%% Main optimization 
close all
clear all
clc

%% Set up initial point and boundaries
load('EOSS.mat','w0','lb','ub')
filename='SSresults.mat';
%[w0,lb,ub] = Initialp(filename);

%% Set up results structure
ActCon=struct('Flag',[],'T', [],'MBOG',[],'QEnv',[], 'Act', [],...
    'Objective',[],'QBOG', [],'Cross',[],'Flow',[],'SplitFlow',[],...
    'LowP',[],'HighP',[],'T3',[],'Output',[],'w', [],'Time',[]);
tol = 1e-3; % Tolerance for a constraint to be active
m=1;

% Set up perturbations order based on QBOG
QB = [];
for p1= 0:5:10
    for p2= 0.9:0.1:1.1
        Tin_bog = -32 + 273.15+p1; %Inlet temperature [K]
        m_bog=610*p2; %Gas mass flow [Kg/h] (Neskaa, 2010)
        QB(m,1) = -BOGC( Tin_bog, m_bog );
        QB(m,2)=p1;
        QB(m,3)=p2;
        m=m+1;
    end
end
QB = sortrows(QB,1);

%% Solve perturbed optimization problems
m =1;
for j= 1:size(QB,1) % Vary on QBOG
    for k= 0.9:0.05:1.1 % Vary on QEnv
        options = optimoptions('fmincon','Algorithm', 'interior-point',...
            'Diagnostics','on','PlotFcns',{@optimplotfval,... 
            @optimplotfirstorderopt,@optimplotconstrviolation,...
            @optimplotstepsize},'FunValCheck','on' ,'TolCon',...
            1e-4,'TolFun', 1e-5, 'MaxIter', 1000,'MaxFunEvals', 150000,...
            'TolX', 1e-8, 'ScaleProblem','obj-and-constr',...
            'UseParallel',true,'FinDiffType','forward',...
            'OutputFcn',@outfun);
        
        %Solve optimization
        tic
        [w,f,eflag,outpt,w0] = runobjconst(w0,lb,ub,options,...
            QB(j,2),QB(j,3),k,1);
        t=toc        
        
        % Save optimal and feasible points
        if eflag == 1 || eflag == 2 
            [~,cons,~,Wc,QBOG,Cross,fil,THA_01]=MainEO(w,...
                QB(j,2),QB(j,3),k);
            A = [0];
            for n=drange(31:1:38)
                if abs(cons(n))<= tol
                    A = [A n];
                end
            end
        else
            A = [-1];
            Wc= [-1];
        end
        
        % Save information in a structure
        ActCon(m).Flag=eflag;
        ActCon(m).T=-32 + 273.15+QB(j,2);
        ActCon(m).MBOG=610*QB(j,3);
        ActCon(m).QEnv=80.5852*k;
        ActCon(m).Act=[A];
        ActCon(m).Objective=Wc;
        ActCon(m).QBOG=-QBOG;
        ActCon(m).Cross=Cross;
        ActCon(m).Output=outpt;
        ActCon(m).Flow=w(130);
        ActCon(m).SplitFlow=w(131);
        ActCon(m).LowP=w(132);
        ActCon(m).HighP=w(133);
        ActCon(m).T3=THA_01;
        ActCon(m).w=w;
        ActCon(m).Time=t;
        m=m+1
    end
end

%% Save results to file
file='Actcons.mat';
save(file)
