%Prico_Find_Constraints

%Find where constraints change from active to inactive for the PRICO
%process


clear all
clc

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

%Load heat exchanger area, scaling matrices and nominal point
load nominal.mat
load hx_area.mat

load opt_nom.mat

Reset_Prico_Temp ;


par.what = 'optim' ;
par.ua = UA ;
par.d0 = d0 ;

%Scaling matrices for constraints and variables

scalex = diag([1e-1 1e-1 1e-3 1e-2 1e-1 1e-1]) ;  %Scaling matrix for inputs
par.scalex = scalex ;
scalec = diag([1e-2 1e-2 1e-3 1e-2 1e-5]) ;  %Scaling matrix for inequality constraints
par.scalec = scalec ;
scaleceq = diag([1 1e-3 1e-3]) ; %Scaling matrix for equality constraints. If solutions are noisy, increase these.
par.scaleceq = scaleceq;


%Prepare scaling matrices, starting point and upper and lower bounds

par.scalec = scalec ; par.scaleceq=scaleceq ;
lb = [ -15 -175 220 120 250/60 9]' ;
ub = [45 -150 2500 700 1200/60 22]' ;
par.x0 = scalex*x_opt_prico ;
par.lb = scalex*lb ;
par.ub = scalex*ub ;


% Constraints: 

% 1: Margin to saturation at turbine outlet
% 2: Superheating at compressor inlet
% 3: Compressor surge margin
% 4: Minimum delta T in heat exchanger
% 5: Maximum work constraint
% 6: Maximum compressor speed constraint

%Define which constraint we are seeking to find whether is active or not

par.constraintID = 2 ; %Number from 1 to 6 picking from the list above

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

% runwhat = 'nominal' ;
runwhat = 'optimize' ;
% runwhat = 'find' ;

d = d0 ;
%If we want to change the other disturbance, do it here:

d(1) = 30; %Change the temperature
d(3) = 1.0 ; %Change feed rate

%Define which disturbance we are changing and which one is constant
par.var = 'feed' ; % 'tamb' or 'feed'

par.d = d ;
par.runwhat = runwhat ;

switch runwhat
    case('nominal')
        nominal_optimization ;
    case('optimize')
        %close all

        dis = [1:0.005:1.075 1.0775:0.0025:1.1025] ;

        
        XOPT = zeros(6,length(dis));
        JOPT = zeros(1,length(dis));
                 ETA = zeros(1,length(dis)) ;
                PH = zeros(size(ETA));
        COPT = zeros(5,length(dis));
        LOPT = zeros(6,length(dis));
        ROPT = zeros(6,length(dis));
        flags = zeros(1,length(dis));
        NoOfSuccesses = 0 ;
        
        exitflag = 0 ; par.exitflag = exitflag ;
        tic
        for i=1:length(dis)
            if exitflag >=1
                par.x0 = xopt ;
            end
            switch(par.var)
                case('tamb')
                    disp(['Optimization number ' num2str(i) ' of ' num2str(length(dis))])
                    disp(['T_amb = ' num2str(dis(i)) ' degrees C'])
                case('feed')
                    disp(['Optimization number ' num2str(i) ' of ' num2str(length(dis))])
                    disp(['Feed rate = ' num2str(dis(i)) ' times the nominal value'])
            end

            [res,xopt,jopt,con,l,exitflag,ETA(i),PH(i)] = prico_optres(dis(i),par) ;
            XOPT(:,i) = scalex\xopt ;
            JOPT(i) = 1e5*jopt ;
            COPT(:,i) = scalec\con(1:5) ;
            ROPT(:,i) = res ;
            LOPT(:,i) = l ;
            flags(i) = exitflag ;
            par.exitflag = exitflag ;
            
        end
        opttime = toc  %#ok<NOPTS>
        avgtime = opttime/length(dis) %#ok<NOPTS>
        ind_2 = find(flags==1) ;
        dis_short = dis(ind_2) ;
        XOPT_short = XOPT(:,ind_2)     ;
        JOPT_short = JOPT(ind_2)       ;
        COPT_short = COPT(:,ind_2)     ;
        LOPT_short = LOPT(:,ind_2)     ;
        ROPT_short = ROPT(:,ind_2)     ;
        ETA_short = ETA(ind_2) ;
        ph_short = PH(ind_2) ;
%         if length(dis) > 1 && length(dis_short) > 1
%             plot_stuff
%         end
        
        disp([flags ; dis]') ; % Show exitflag together with disturbance value to see when the optimizer starts failing
    case('find')
        par.options.Display = 'final' ;
        exitflag = 0 ; par.exitflag = exitflag ;
        options = optimset('TolX',1e-4,'Display','iter','MaxIter',15) ; % Options for fsolve
        solution= fzero(@prico_optres,[26.1 30],options,par)
end

Reset_Prico_Temp ;