%This file is for running optimization of the PRICO process, using the
%PRICO_TEMPSPECS model.

close all
%clear all
clc

global h
h=actxserver('UnisimDesign.Application');

% Load the nominal values of x and d
load nominal.mat
load hx_area.mat
load scaling.mat
load opt_nom.mat

Reset_Prico_Temp ;

% Define parameters
par.scalec = scalec ; 
par.scaleceq = scaleceq; 
par.scalex = scalex ; 
par.scalex_small = scalex_small ;
par.ua = UA ;
par.what = 'optim' ;
par.composition = 'no' ;

switch par.composition 
    case('yes')
        lb = [ -2 -159 1000 400 800/60 0.1 0.1 0 0.1 17]' ;
        ub = [20 -155 1200 500 1200/60 0.4 0.4 0.1 0.4 21]' ;
        %Scale x0 and upper and lower bounds
        x_0 = scalex*x0 ;
        lb = scalex*lb ;
        ub = scalex*ub ;
    case('no')
%         scalec = removerows(scalec,4) ;
%         ind = [1 2 3 5 6];
%         scalec = scalec(:,ind) ;
        par.scalec = scalec ;
        lb = [ -10 -175 250 150 250/60 9]' ;
        ub = [35 -153 2000 500 1200/60 22]' ;
        x_0 = scalex_small*x_opt_prico ;
        lb = scalex_small*lb ;
        ub = scalex_small*ub ;
end

opt_options = optimset('Algorithm','interior-point','TolCon',1e-4,'TolFun',1e-4,'Display','iter','MaxIter',125,'MaxFunEvals',1e4,'FinDiffType','central') ; % Options for fmincon

% Optimization at nominal disturbances with randomly generated starting
% points x0 inside [lb,ub]:

%nominal_optimization ;

% Do optimization

pert = [1:1/15:4/3]; %#ok<NBRAK> %

%  load start.mat ;
% newx0 = x_start ;

disturb = zeros(size(pert));

XOPT = zeros(6,length(pert));
JOPT = zeros(1,length(pert));
ETA = zeros(1,length(pert)) ;
PH = zeros(size(ETA));
COPT = zeros(5,length(pert));
flags = zeros(1,length(pert));
NoOfSuccesses = 0 ;

%Disturbances are:

% d(1): Ambient temperature (used to set temperature of the two hot streams
% entering the HX)
% d(2): Natural gas feed pressure
% d(3): Natural gas feed flow rate

j = 1 ; % Which one to change

iterations = [];
optimtime = [];

for i=1:length(pert)
    d = d0 ;
    d(3) = d0(3)*1 ; 
    
    d(j) = d0(j)*pert(i) ;
    disturb(i) = d(j) ; 
    disp(['Optimization number ' num2str(i) ' of ' num2str(length(pert))])
    switch j
        case(1)
            disp(['Ambient temperature: ' num2str(d(j)) ' deg C'])
        case(2)
            disp(['Feed pressure: ' num2str(d(j)/100) ' bar'])
        case(3)
            disp(['Feed flow rate: ' num2str(d(j)*3600) ' kmol/h'])
    end
    tic
     [x_opt, fval, exitflag,output] = fmincon(@obj_prico_temp,x_0,[],[],[],[],lb,ub,@constraints_prico_temp,opt_options,d,par);
%	[x_opt, fval, exitflag,output] = fmincon(@obj_prico_temp,newx0,[],[],[],[],lb,ub,@constraints_prico_temp,opt_options,d,par);
    opttime = toc ;
    [j_opt eta ph] = obj_prico_temp(x_opt,d,par) ;
    c_opt = constraints_prico_temp(x_opt,d,par) ;

    XOPT(:,i) = scalex_small\x_opt ;
    JOPT(i) = 1e5*j_opt ;
    COPT(:,i) = c_opt ; %scalec\
    ETA(i) = eta/100 ; 
    PH(i) = ph/100 ; % Go from kPa to bar
    flags(i) = exitflag ;
    
    if exitflag > 0
        NoOfSuccesses = NoOfSuccesses + 1;
        iterations = [iterations output.iterations] ; %#ok<AGROW> % Count iterations of successful optimizations
        optimtime = [optimtime opttime] ; %#ok<AGROW>
        x_0 = x_opt ; % Use the last successful optimization result as starting point! Has to be better than using the nominal point.
    end
    
end

NoOfSuccesses %#ok<NOPTS>
avg_iter = mean(iterations) %#ok<NOPTS>
avt_time = mean(optimtime) %#ok<NOPTS>
%Plot results of the converged optimizations
ind_2 = find(flags==1) ; % Single out the converged optimizations

disturb_short = disturb(ind_2) ;
XOPT_short = XOPT(:,ind_2)     ;
JOPT_short = JOPT(ind_2)       ;
COPT_short = COPT(:,ind_2)     ;
ETA_short = ETA(ind_2)         ;
PH_short = PH(ind_2)           ;

PL = XOPT_short(4,:)/100 ;

if length(pert) > 1 && length(disturb_short) > 1

%Plot constraint functions
figure('Name','Constraints')
subplot(5,1,1),plot(disturb_short,COPT_short(1,:),'LineWidth',2),title('DPsat')
subplot(5,1,2),plot(disturb_short,COPT_short(2,:),'LineWidth',2),title('DTsup')
subplot(5,1,3),plot(disturb_short,COPT_short(3,:),'LineWidth',2),title('DMsurge')
subplot(5,1,4),plot(disturb_short,COPT_short(4,:),'LineWidth',2),title('DTmin')
subplot(5,1,5),plot(disturb_short,COPT_short(5,:),'LineWidth',2),title('W-Wmax')

figure('Name','Objective function'),plot(disturb_short,JOPT_short) ; %Plot the objective function

opt_speed = 100*XOPT_short(5,:)*scalex_small(5,5)/ub(5) ;
% max_speed = ub(5)/scalex_small(5,5)*ones(size(opt_speed)) ;
% min_speed = lb(5)/scalex_small(5,5)*ones(size(opt_speed)) ;

% speed = 60*[min_speed ; opt_speed ; max_speed] ;

figure('Name','CompressorVariables')
subplot(4,1,1),plot(disturb_short,opt_speed,'LineWidth',2),legend('Optimal speed (%)') ;
subplot(4,1,2),plot(disturb_short,ETA_short,'LineWidth',2),legend('Adiabatic efficiency')
subplot(4,1,3),plot(disturb_short,PL,'LineWidth',2),legend('Compressor inlet pressure')
subplot(4,1,4),plot(disturb_short,PH_short,'LineWidth',2),legend('Compressor outlet pressure')


lb6 = 3600*lb(6)/scalex_small(6,6);
min_flow = lb6*ones(size(XOPT_short(6,:)));
ub6=3600*ub(6)/scalex_small(6,6);
max_flow = ub6*ones(size(XOPT_short(6,:)));
flow = 3600*XOPT_short(6,:) ;

flow = [min_flow ; flow ; max_flow];

figure('Name','MR flow rate'),plot(disturb_short,flow,'LineWidth',2),legend('Minimum flow rate','Refrigerant flow rate','Maximum flow rate')
%figure(6),plot(disturb(ind_2),PH(ind_2),'LineWidth',2),legend('Compressor outlet pressure')

end

Reset_Prico_Temp ;