function OptReconLargeReformer()
%
% Add subdirectory paths
cwd=pwd;
addpath([cwd '\Properties']);
addpath([cwd '\Reactions']);
addpath([cwd '\UnitModels']);
addpath([cwd '\Misc']);
addpath([cwd '\xIndex']);

% Defenition of global variables used in model equations
GlobalVars;

global Vr                           % Matrix of catalyst volumes row : Reactor, column : element in reactor
global SolverFlag                   % Optimization flag (0= Solve model, 2=Optimization, 1=Reconcilliation)

Unscaled=0;
EvalScale=0;

NC   = 5;                           % Number of components in the model
NStr =55;                           % Number of streams in the model
NCSTR=40;                           % Vector of elements in each reactor
NHtr = 4;                           % Number of heaters in the model
NCmp = 1;                           % Number of compressors in the model
NHex = 0;                           % Number of heat exchangers in the model
NCon = 0;                           % Number of condesers in the model
NHex2= 1;                           % Number of type 2 heat exchangers in the model
NHex3= 1;                           % Number of type 3 heat exchangers in the model
NConstr=1;                          % Constraint unit model calc

Vr   = ones(1,10);                  % Reactor element volumes                            [m3]

NVar = (NC+3)*NStr + 1*NCSTR + 1*NHtr + 3*NCmp + 1*NHex + 1*NCon + 2*NHex2 + 5*NHex3 + 7*NConstr;

Sv   = Make_Sv;
if Unscaled
    Sv=speye(size(Sv));
end

%----------------------------------------------------------------
disp(' ')
disp(' ')
disp('-------------------------------------------------------')
disp('Simulation (fmincon)')

% Initial values
SolverFlag=0;                   %Run SSSpecs
Vars=InitVal;

% Set solver options
Opt=optimset('fmincon');
Opt.DerivativeCheck   = 'off';
Opt.Diagnostics       = 'off';
Opt.Display           = 'iter';
Opt.GradConstr        = 'on';
Opt.GradObj           = 'on';
Opt.LargeScale        = 'off';
Opt.TypicalX          = Sv\Vars;

% Scaled initial values
Var0q=Sv\Vars;

% Linear eq. constraints
Aeq=full(SSSpecs(Var0q,'Jacobian')*Sv);
beq=SSSpecs(Var0q,'b');

A=[]; b=[];

% Make upper and lower bounds
lb=[];  ub=[];

SolverFlag=1;                   %Simulation
tic
[Varq,fval,exitflag,output,lambda,grad,hessian] = fmincon(@Sim,Var0q,A,b,Aeq,beq,lb,ub,@NonLinCon,Opt);
toc

% Write output
iter=output.iterations;
SolverFlag=0;                   
[r,J]=fsFun(Varq);
disp('Final solution')
disp(['Rank J : ' num2str(rank(full(J)),3)])
disp(['Cond J : ' num2str(condest(J),3)])
disp(['Max  J : ' num2str(max(max(J)))])
disp(['Min  J : ' num2str(min(min(J)))])

% Save results
Vars=Sv*Varq;
VarSim=Vars;
save VarSim Vars

%----------------------------------------------------------------
disp(' ')
disp(' ')
disp('-------------------------------------------------------')
disp('Reconciliation')

Opt=optimset('fmincon');
Opt.DerivativeCheck   = 'off';
Opt.Diagnostics       = 'off';
Opt.Display           = 'iter';
Opt.GradConstr        = 'on';
Opt.GradObj           = 'on';
Opt.LargeScale        = 'off';
Opt.MaxIter           = 300;
Opt.TypicalX          = Sv\Vars;

% Initial values
% Use simulation solution as initial value
Var0q=Sv\VarSim;

%Linear eq.constraints
Aeq=full(ReconSpecs(Var0q,'Jacobian')*Sv);
beq=ReconSpecs(Var0q,'b');
A=[]; b=[];

%Upper and lower bounds
[ub lb]=UpLowBounds;
ub=Sv\ub;
lb=Sv\lb;

% Make measurement vector
MakeMeas;

SolverFlag=1;                   %Reconciliation
tic
[Varq,fval,exitflag,output,lambda,grad,hessian] = fmincon(@Rec,Var0q,A,b,Aeq,beq,lb,ub,@NonLinCon,Opt);
toc

SolverFlag=0;                   
[r,J]=fsFun(Varq);

% Write output
disp('Final solution')
disp(['Rank J : ' num2str(rank(full(J)),3)])
disp(['Cond J : ' num2str(condest(J),3)])
disp(['Max  J : ' num2str(max(max(J)))])
disp(['Min  J : ' num2str(min(min(J)))])

% Save results
Vars=Sv*Varq;
VarRec=Vars;
save VarRec VarRec

save ReconLargeReformerData % for use in excel

%----------------------------------------------------------------
disp(' ')
disp(' ')
disp('-------------------------------------------------------')
disp('Optimization')

% Used reconcilliation result as initial value
Var=VarRec;
Var0q=Sv\Var;

Opt=optimset('fmincon');
Opt.DerivativeCheck   = 'off';
Opt.Diagnostics       = 'off';
Opt.Display           = 'iter';
Opt.GradConstr        = 'on';
Opt.GradObj           = 'on';
Opt.LargeScale        = 'off';
Opt.TypicalX          = Sv\Var;

%Linear eq. constraints
Aeq=full(OptSpecs(Var0q,'Jacobian')*Sv);
beq=OptSpecs(Var0q,'b');
A=[]; b=[];

%Upper and lower bounds
[ub lb]=UpLowBounds;

% Optimization constraints
ix=IxStr(5);
ub(ix(NC+2))=900;        % Heater 1 Outlet temperature
lb(ix(NC+2))=773;        % Heater 1 Outlet temperature
ix=IxStr(16);
ub(ix(NC+2))=900;        % Heater 2 Outlet temperature
lb(ix(NC+2))=773;        % Heater 2 Outlet temperature
ix=IxStr(27);
ub(ix(NC+2))=900;        % Heater 3 Outlet temperature
lb(ix(NC+2))=773;        % Heater 3 Outlet temperature
ix=IxStr(38);
ub(ix(NC+2))=900;        % Heater 4 Outlet temperature
lb(ix(NC+2))=773;        % Heater 4 Outlet temperature

ix=IxHtr(1);
ub(ix(1))=   9500;       % Heater 1 duty
ix=IxHtr(2);
ub(ix(1))=  13000;       % Heater 2 duty
ix=IxHtr(3);
ub(ix(1))=   13000;      % Heater 3 duty
ix=IxHtr(4);
ub(ix(1))=   10000;      % Heater 4 duty

ix=IxConstr(1);
ub(ix(1))= 150;          % RON
lb(ix(1))= 103;          % RON
ub(ix(2))=   4;          % H2/HC
lb(ix(2))=   3;          % H2/HC
ub(ix(3))=   100;        % Feed Massflow
lb(ix(3))=   50;         % Feed Massflow
lb(ix(7))=    3.5;       % Hydrogen product massflow

ix=IxStr(53);
ub(ix(NC+3))=10;         % D-1405 Pressure
lb(ix(NC+3))=7;          % D-1405 Pressure

ub=Sv\ub;
lb=Sv\lb;

SolverFlag=2;            %Optimization
tic
[Varq,fval,exitflag,output,lambda,grad,hessian] = fmincon(@Obj,Var0q,A,b,Aeq,beq,lb,ub,@NonLinCon,Opt);
toc

% Write output
SolverFlag=0;                        
[r,J]=fsFun(Varq);
disp('Final solution')
disp(['Rank J : ' num2str(rank(full(J)),3)])
disp(['Cond J : ' num2str(condest(J),3)])
disp(['Max  J : ' num2str(max(max(J)))])
disp(['Min  J : ' num2str(min(min(J)))])

% Save the results
Vars=Sv*Varq;
VarOpt=Vars;
save VarOpt VarOpt

save OptLargeReformerData % for use in excel
% Remove paths
cwd=pwd;
rmpath([cwd '\Properties']);
rmpath([cwd '\Reactions']);
rmpath([cwd '\UnitModels']);
rmpath([cwd '\Misc']);
rmpath([cwd '\xIndex']);

return
% --------------------------------------------------------------------------

function [r,J]=fsFun(Var)
  % Steady state model 
  % Returns residuals, r=f(x), and Jacobian, J=df/dx.
  r=f(Var);
  if nargout > 1
      J=full(Jacobian(Var));
  end
return         

function [g,dg]=Obj(Var)
  % Optimization objective 
  % g  :objective to be minimized
  % dg :gradient of g at var
  
  if nargout == 1
      [g]   =Objective(Var);
  else
      [g dg]=Objective(Var);
  end
return         

function [g,dg]=Sim(Var)
  % Optimization objective 
  % g  :objective to be minimized
  % dg :gradient of g at var
  
  if nargout == 1
      g=0;
  else
      g=0; dg=zeros(size(Var))';
  end
return         

function [g,dg]=Rec(Var)
  % Optimization objective 
  % g  :objective to be minimized
  % dg :gradient of g at var
  
  if nargout == 1
      [g]   =Recon(Var);
  else
      [g dg]=Recon(Var);
  end
return         

function [c,ceq,Gc,Gceq] = NonLinCon(Var)
  % Steady state model 
  % c    :nonlinear constraints
  % ceq  :nonlinear equality constraints
  % GC   :gradient of c
  % GCeq :gradient of ceq
  t=0;
  
  c =[];
  Gc=[];
  ceq=f(Var);
  if nargout > 2
      Gceq=full(Jacobian(Var)');
  end
return         


function Var0=InitVal
  % Find initial values
  global NVar
  Var0=zeros(NVar,1);
  r=Model(Var0,'InitVal');
  
  for i=1:NVar:size(r,1)
      Var0=Var0+r(i:i+NVar-1);
  end
return
  
function dxdt = f(Var)
  % Compute xdot in M*xdot = f(t,x)
  global Sv
  dxdt=Model(Sv*Var,'dxdt');
return

function M=Mmtx(Var)
  % Compute mass matrix
  global NVar
  global Sv
  x=zeros(NVar,1);
  M=Model(Sv*Var,'Mass')*Sv;
return

function J=Jacobian(Var)
   %Compute racobian matrix d(f(t,x))/dx
   global Sv
   J=Model(Sv*Var,'Jacobian')*Sv;
return
   
% --------------------------------------------------------------------------
% Model
function res=Model(Var,Mode)
    Kca=0.2;    % Mass of catalyst tuning factor
    V1=1.00*Kca;
    V2=1.07*Kca;
    V3=1.78*Kca;
    V4=3.28*Kca;
    res=[ ...
         Specs(Var,Mode)     
         VapLiqMix(Var,1,55,1,2,3,Mode)
         E1401(Var,1,2,3,4,48,49,50,Mode)
         Heater(Var, 1, 4, 5   ,Mode)
         CSTR(  Var, 1, 5, 6,V1,Mode)
         CSTR(  Var, 2, 6, 7,V1,Mode)
         CSTR(  Var, 3, 7, 8,V1,Mode)
         CSTR(  Var, 4, 8, 9,V1,Mode)
         CSTR(  Var, 5, 9,10,V1,Mode)
         CSTR(  Var, 6,10,11,V1,Mode)
         CSTR(  Var, 7,11,12,V1,Mode)
         CSTR(  Var, 8,12,13,V1,Mode)
         CSTR(  Var, 9,13,14,V1,Mode)
         CSTR(  Var,10,14,15,V1,Mode)
         Heater(Var, 2,15,16   ,Mode)
         CSTR(  Var,11,16,17,V2,Mode)
         CSTR(  Var,12,17,18,V2,Mode)
         CSTR(  Var,13,18,19,V2,Mode)
         CSTR(  Var,14,19,20,V2,Mode)
         CSTR(  Var,15,20,21,V2,Mode)
         CSTR(  Var,16,21,22,V2,Mode)
         CSTR(  Var,17,22,23,V2,Mode)
         CSTR(  Var,18,23,24,V2,Mode)
         CSTR(  Var,19,24,25,V2,Mode)
         CSTR(  Var,20,25,26,V2,Mode)
         Heater(Var, 3,26,27   ,Mode)
         CSTR(  Var,21,27,28,V3,Mode)
         CSTR(  Var,22,28,29,V3,Mode)
         CSTR(  Var,23,29,30,V3,Mode)
         CSTR(  Var,24,30,31,V3,Mode)
         CSTR(  Var,25,31,32,V3,Mode)
         CSTR(  Var,26,32,33,V3,Mode)
         CSTR(  Var,27,33,34,V3,Mode)
         CSTR(  Var,28,34,35,V3,Mode)
         CSTR(  Var,29,35,36,V3,Mode)
         CSTR(  Var,30,36,37,V3,Mode)
         Heater(Var, 4,37,38   ,Mode)
         CSTR(  Var,31,38,39,V4,Mode)
         CSTR(  Var,32,39,40,V4,Mode)
         CSTR(  Var,33,40,41,V4,Mode)
         CSTR(  Var,34,41,42,V4,Mode)
         CSTR(  Var,35,42,43,V4,Mode)
         CSTR(  Var,36,43,44,V4,Mode)
         CSTR(  Var,37,44,45,V4,Mode)
         CSTR(  Var,38,45,46,V4,Mode)
         CSTR(  Var,39,46,47,V4,Mode)
         CSTR(  Var,40,47,48,V4,Mode)
         E1402_D1405(Var,1,49,50,51,52,Mode)
         Split(Var,51,53,54,Mode)
         Compressor(Var,1,54,55,Mode)
         ConstraintCalc(Var,Mode)
        ];
    if imag(res)~=0
        warning('Complex numbers in residual vector')
    end
return

function r=Specs(Var,Mode)
global SolverFlag
  if SolverFlag==0
    r=SSSpecs(Var,Mode);
  else
    r=[];    
  end
return

