%% Simple Ammonia Loop - Optimization of the combined flowsheet.
% This file is used to evaluate the performance of the surrogate-based 
% optimization flowsheet. The surrogate models are defined in the folder
% "Prog1. Adaptive Smapling". There are 3 different surrogate models safed
% in the folder "LocalSave". If desired, the surrogate models can be
% changed in the file "ParOpt.m".
%
% Written by Julian Straus, Trondheim 09.12.2018
%
%  The components in this system are defined as 
%   x(1) - Hydrogen
%   x(2) - Nitrogen
%   x(3) - Ammonia
%   x(4) - Argon
%   x(5) - Methane

% Clearing of the workspace and the memory, closing of all open
% windows/figures
clc, clear variables, close all

% Start of the timer for the overall function file
time.in.overall = tic;

% Definition of the number of reactors in the CSTR cascade. Currently, the
% initialization requires 10 CSTRs for each reactor. If this should be
% changed, then the initialization in InitSamp has to be adjusted.
Reac1_n  = 10;
Reac2_n  = 10;

% Loading of the needed parameters
ParOpt

% Definition of the inlet values for the makeup section
def.M.FMake = [1809.17; 27.32e5; 7.495;...                  Feed Stream Q [mol/s], p [Pa], and T [C]
              [73.02; 24.28; 1.62; 0.49; 0.59]/100;...      Feed Stream mole fractions [-]
              par.M{1}.S];...                               Feed Compressor RPM
def.M.HEx1  = [9274.61; 10e5; -10;...                       Water Stream HEx1 Q [mol/s], p [Pa], and T [C]
               ones(5,1)];...                               Water Stream HEx1 mole fractions [-]
def.M.FMix  = [6561.38; 127.50e5; 25.82;...              	Mixer Stream Q [mol/s], p [Pa], and T [C]
              [64.77; 20.25; 2.66; 6.98; 5.34]/100];...     Mixer Stream Q mole fractions [-]
                        
% Definition of the inlet values for the reaction section
def.R.FReact = [8369.7; 140e5; 51.62;...                    Feed Stream Q [mol/s], p [Pa], and T [C]
                0.6658; 0.2113; 0.0243; 0.0557; 0.0429];... Feed Stream mole fractions [-]
def.R.Val1   = 832.4;...                                    Split Valve Q [mol/s]
def.R.Temp   = [379.1; 410.9];...                           Temperature Reactor Inlets [C]
def.R.Tref   = par.R{17}.Tref+5;...                           Reference Temperature HEx5
def.R.Sratio = {par.R{3}.ratio+[-5; -5; 0]};...                         Split ratios Split 1

% Calculation of the inlet values for the separation section
def.S.FSep = [7539.40; 131.2e5; 53.73;...                   Feed Stream Q [mol/s], p [P a], and T [C]
              0.5728; 0.1791; 0.1385; 0.0619; 0.0477];...   Feed Stream Mole fractions [-]
def.S.HEx6 = [11184.22; 3.828e5; 10;...                     Water Stream HEx6 Q [mol/s], p [Pa], and T [C]
              ones(5,1)];...                                Water Stream HEx6 mole fractions
def.S.Tref = [par.S{9}.Tref+4 ...                             Reference Temperature HEx8
              par.S{11}.Tref+4 ...                            Reference Temperature HEx9
              par.S{15}.Tref+4];...                           Reference Temperature HEx10
def.S.Sratio = {par.S{3}.ratio-[0;-10];...                          Split ratios Split 2
                par.S{7}.ratio+[0;0.5]};...                         Split ratios Split 3
                
% Calculation of version histories
Ver.SM    = [1;...
             1];

% Values for the cost function
pr.feed = 8.8*0.08*8e3/1000;...                 Feed cost
pr.prod = -17*0.2*1e4/1000;...                  Product cost
pr.purg = -11.2*0.01*1e4/1000;...               Purge cost
pr.comp = 0.08*1e4/4000;...                     Compressor cost
pr.stea = -0.1;...                              Steam production cost
pr.cool = .5*0.075;...                             Cooling cost
% Calculations
%   Flow: Molar mass*cost Antonio*scaling factor to mol*scaling factor to kmol
%   Comp: Cost Antonio*scaling factor to kW, COst adjustment for wrong
%   value in Antonio's code

% Definition of the initial conditions and the anonymous functions
[initComb,      ConAnComb,...
 initComb_NN,   ConAnComb_NN]...
= InitOpt(par,scl,con,ther,def,itg,dirmain);
    
load('InitCombOpt.mat')

% Definition of the opt.fsol used in each solver 'iter-detailed'
opt.fmin = optimoptions('fmincon','Algorithm','sqp','Display','off', ...
                        'SpecifyConstraintGradient',true,'SpecifyObjectiveGradient',true, ...
                        'MaxIter',1e3,'MaxFunEvals',1e10,...
                        'UseParallel',true);
opt.fmin.Display = 'iter-detailed';

% Addition of the  Path direction of the unit operations
UO = [dirmain,'\Util4. Unit Operations\Optimization'];
addpath(UO)

% Addition of an additional extra line break
fprintf('\n\n')
time.in.overall = tic;
    
%% Code for calculating the combined flowsheet using the initial models

if Ver.SM(1)
    fprintf('Combined Flowsheet - Original Model:\n')
    
    % Solving of the system
    time.in.C = tic;
    for i = 1:pow
    [a.C,fval.C,exitflag.C,output.C] ...
                = fmincon(@(x) objfun_Comb(x,pr),initComb,[],[],[],[],[],[],ConAnComb,opt.fmin);
    end
    time.C = toc(time.in.C);
    if exitflag.C > 0
        fprintf('It took %4.4f s to solve the problem.\n',time.C)
    elseif exitflag.C == 0
        fprintf('The problem exceeded the number of predefined iterations.')
    else
        fprintf('The problem did not converge.')
    end

    % Reassignment of the values of each UO
    [into.C,UOVec.C,Par.C,Q.C] = ReassigFlex(a.C,ConAnComb,name.Comb,scl,itg);
    [~,~,Par.C_init,~]         = ReassigFlex(initComb,ConAnComb,name.Comb,scl,itg);

    fprintf('\n')
end
     
%% Code for calculating the combined flowsheet using the PLS + neural network model

if Ver.SM(2)
    fprintf('Combined Flowsheet - Neural Network:\n')
        
    % Solving of the system
    time.in.C = tic;
    for i = 1:pow
    [a.CNN,fval.CNN,exitflag.CNN,output.CNN] ...
                = fmincon(@(x) objfun_Comb_NN(x,pr),initComb_NN,[],[],[],[],[],[],ConAnComb_NN,opt.fmin);
    end
    time.CNN = toc(time.in.C);
    if exitflag.CNN > 0
        fprintf('It took %4.4f s to solve the problem.\n',time.CNN)
    elseif exitflag.CNN == 0
        fprintf('The problem exceeded the number of predefined iterations.')
    else
        fprintf('The problem did not converge.')
    end

    % Reassignment of the values of each UO
    [into.CNN,UOVec.CNN,Par.CNN,Q.CNN] = ReassigFlex(a.CNN,ConAnComb_NN,name.CombNN,scl,itg);
    [~,~,Par.CNN_init,~]               = ReassigFlex(a.CNN,ConAnComb_NN,name.CombNN,scl,itg);

    fprintf('\n')

end    
time.overall = toc(time.in.overall);

%% Calculation of the comparison in-between the models

% Comparison between stream values and duties
if Ver.SM
    
    % Makeup Section
    Comp.MaCAbs = abs([Q.C.Comp1; Q.C.Comp2; Q.C.Comp3] - [Q.CNN.Comp1; Q.CNN.Comp2; Q.CNN.Comp3]);
    Comp.MaCRel = Comp.MaCAbs./[Q.C.Comp1; Q.C.Comp2; Q.C.Comp3];
    
    % Reactor section
    Comp.ReaAbs = abs((UOVec.C.HEx2(2:3,:)-UOVec.CNN.NN_Reac));
    Comp.ReaRel = abs(Comp.ReaAbs./UOVec.C.HEx2(2:3,:))*100;
    Comp.ReaQAbs = abs(Q.C.HEx5-Q.CNN.NN_Reac);
    Comp.ReaQRel = Comp.ReaQAbs./Q.C.HEx5*100;

    % Separation Section
    Comp.SepAbs = abs([(UOVec.C.HEx6(4,:)-UOVec.CNN.NN_Sep(5,:));...
                       (UOVec.C.Sep(3,:) -UOVec.CNN.NN_Sep(2,:));...
                       (UOVec.C.HEx7(4,:)-UOVec.CNN.NN_Sep(3,:))]);
    Comp.SepRel = abs(Comp.SepAbs./[UOVec.C.HEx6(4,:);...
                                    UOVec.C.Sep(3,:);...
                                    UOVec.C.HEx7(4,:)])*100;
    Comp.SepQAbs = abs([Q.C.HEx8;Q.C.HEx9;Q.C.HEx10]-Q.CNN.NN_Sep);
    Comp.SepQRel = Comp.SepQAbs./[Q.C.HEx8;Q.C.HEx9;Q.C.HEx10]*100;

    Comp.MaCRel
    Comp.ReaRel
    Comp.ReaQRel
    Comp.SepRel
    Comp.SepQRel
    (sum(Q.CNN.NN_Sep)-sum([Q.C.HEx8;Q.C.HEx9;Q.C.HEx10]))/sum([Q.C.HEx8;Q.C.HEx9;Q.C.HEx10])*100
end


% Comparion of the parameters
fprintf('Compressor Speed')
[Par.C.Comp1 Par.CNN.Comp1]*1e4
fprintf('Compressor Speed Error')
abs(Par.C.Comp1-Par.CNN.Comp1)*1e4/(10493-9980)*100

fprintf('Split Ratio Reactor')
[Par.C.Split1(1:2) Par.CNN.NN_Reac(1:2)]*100

fprintf('Reference Temperature Reactor')
[Par.C.HEx5 Par.CNN.NN_Reac(3)]*1000

fprintf('Split Ratio pre Cooling Train')
[Par.C.Split2(2) Par.CNN.NN_Sep(1)]*100

fprintf('Reference Temperature Cooling Train')
[[Par.C.HEx8; Par.C.HEx9; Par.C.HEx10] Par.CNN.NN_Sep(2:4)]*100

fprintf('Purge Split')
[Par.C.Split3(2) Par.CNN.Split3(2)]*100

fprintf('Cost function [%%]')
(fval.C-fval.CNN)/fval.C*100

