
% This file is used to calculate the inlet flowrate to the reactor in the
% case of a given feed and hence extent of reaction. It still contains
% redundant parts which can be removed, but these parts do not increase the
% calculation drastically.

% Clearing of the workspace and closing of all windows
clc; clear; close all;

% Definition of the current path folder
FileName = mfilename('fullpath');
[directory,~,~] = fileparts(FileName);
[dirmain,~,~] = fileparts(directory);

% Import of CASADI
addpath([dirmain '\Util1. Casadi3.1.0'])
addpath([dirmain '\Util2. Functions'])
addpath(genpath([dirmain '\Util3. MATLAB2tikz']))
addpath([dirmain '\Data'])
import casadi.*

% Definition of decision variables
dec.plot = true;
dec.scan = false;

% Definition of the save name
savename = 'SS_opt.mat';

% Definition of the total integration time and the number of time steps for
% the integrator
t.tot  = 100*60;
t.int = 10;
t.N  = t.tot/t.int;
t.samp = (0:t.int:t.tot)/60;

% Definition of feed variables
in.nom = zeros(4,1);
in.nom(1) = 252/3.6;        % Feed flow (kg/s)
in.nom(2) = 200;            % Pressure [bar]
in.nom(3) = 250;            % Temperature [C]
in.nom(4) = 0.08;           % Mass fraction ammonia [-]

% Definition of the initial guesses and redefinition into x and z
load([dirmain '\Data\Init_In.mat'])

% Load of the model
[func,dae,d,lb,ub,n,Cineq] = AmmonReac_NewDist(t.int);

%% Part a) Root finding for the original system

% Combination of the variables
w_root = [dae.x;dae.z];
w0     = [x0;z0];

% Combination of the equality constraints and substitution of inlet values
g = vertcat(dae.ode,dae.alg);
g_root = substitute(g,d,zeros(n.d,1));
g_root = substitute(g_root,dae.p(1:n.u),p0(1:n.u));

% Definition and solving of the system
g_fun = Function('g_fun',{w_root},{g_root});
G = rootfinder('G','newton',g_fun);
[root.comb] = G('i0',w0);

% Calculation of the number of CSTRs in total
N = size(dae.x,1);

% Definition of indices used in the solution
ind.T = [3*N+1 1:N/3 3*N+2 N/3+1:2*N/3 3*N+3 2*N/3+1:3*N/3]; % 2*N+4 
ind.w = [3*N+5 N+1:4*N/3 3*N+6 4*N/3+1:5*N/3 3*N+7 5*N/3+1:2*N];
ind.x = 1:N;
ind.z = (N+1):(3*N+10);
ind.u = (3*N+11):(3*N+13);

% Reassignment of variables
root.x = full(root.comb.o0(ind.x));
root.z = full(root.comb.o0(ind.z));
root.T = full(root.comb.o0(ind.T));
root.w = full(root.comb.o0(ind.w));
root.f = -(root.z(N)-(in.nom(4)))*(in.nom(1));
x0 = root.x;
z0 = root.z;

% Print of the found results
fprintf('\n\n\nThe root-finding is completed. \n\n\n')

%% Part b) Steady state optimization of the system

% Combination of the variables
w_opt = [dae.x; dae.z; dae.p(1:n.u)];
lbw   = [lb.x;  lb.z;  lb.p(1:n.u)];
ubw   = [ub.x;  ub.z;  ub.p(1:n.u)];
w0    = [x0;    z0;    p0(1:n.u);];

% Combination of the equality and inequality constraints
gadd = [g; Cineq];
lbg  = [zeros(size(gadd,1)-1,1); -1];
ubg  = [zeros(size(gadd))];

% Predefinition of the solution vector
nom_opt.x = zeros(n.x,100);
nom_opt.z = zeros(n.z,100);
nom_opt.T = zeros(N+3,100);
nom_opt.w = zeros(N+3,100);
nom_opt.p = zeros(n.p,100);
nom_opt.u = zeros(n.u,100);
nom_opt.f = zeros(1,100);

l = 0;
nom_opt.f = 0;

for l = 1:100
    % Substitution of the disturbance values
    in.dist = [-(8.6+l*.01);zeros(n.d-1,1)];
    J = substitute(dae.quad,d,in.dist);
    g_opt = substitute(gadd,d,in.dist);

    % Definition of the NLP function
    nlp = struct('x',w_opt,'f',J,'g',g_opt);
    solver = nlpsol('solver','ipopt',nlp);
    nom_opt.comb = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',lbg,'ubg',ubg);
    
    % Reassignment of variables
    nom_opt.x(:,l) = full(nom_opt.comb.x(ind.x));
    nom_opt.z(:,l) = full(nom_opt.comb.x(ind.z));
    nom_opt.T(:,l) = full(nom_opt.comb.x(ind.T));
    nom_opt.w(:,l) = full(nom_opt.comb.x(ind.w));
    nom_opt.u(:,l) = full(nom_opt.comb.x(ind.u));
    nom_opt.p(1:n.u,l) = full(nom_opt.u(:,l));
    nom_opt.f(:,l) = full(nom_opt.comb.f);
    
    % Redefinition of the new initial values
    w0 = full(nom_opt.comb.x);

    if (root.f-nom_opt.f(:,l)) < 0.0001 && nom_opt.f(:,l) < -10
        break
    end
end
if l < 100
    fprintf('\n\n\nThe inlet flowrate of the reactor for a fixed outlet flowrate was found. \n\n\n')
    nom_opt.x = nom_opt.x(:,1:l);
    nom_opt.z = nom_opt.z(:,1:l);
    nom_opt.T = nom_opt.T(:,1:l);
    nom_opt.w = nom_opt.w(:,1:l);
    nom_opt.u = nom_opt.u(:,1:l);
    nom_opt.p = nom_opt.p(:,1:l);
    nom_opt.f = nom_opt.f(:,1:l);
else
    fprintf('\n\n\nThe inlet flowrate was not found within 100 iterations.')
end


%% Part c) Plotting of the results

% Definition of integer values
itg.fig = 1;
itg.lab = 1;
itg.fda = 1;

N1 = N/3+1;
N2 = N/3+1;
N3 = N/3+1;

% Definition of figure and values
fig.dimw = 8.6;
fig.dimh = 6;
fig.Font = [10 8];
fig.label = [];

if dec.plot
    
% Plotting of the steady state reactor profiles
fig.in(itg.fig) = figure(itg.fig);
set(fig.in(itg.fig),'units','centimeters');
set(gcf, 'position', [5+itg.fig,5+itg.fig, fig.dimw, fig.dimh]);
hold on
yyaxis left
fig.axes(itg.fig) = gca;
fig.axes(itg.fig).YAxis(1).Color = 'k';
fig.axes(itg.fig).YAxis(2).Color = 'r';
fig.dat(itg.fda) = plot([linspace(0,1,N1) linspace(1,2,N2) linspace(2,3,N3)],...
                        full(root.T),'k'); itg.fda = itg.fda+1;
fig.dat(itg.fda:itg.fda+l-1) = plot([linspace(0,1,N1) linspace(1,2,N2) linspace(2,3,N3)],...
                                     full(nom_opt.T),'k--'); itg.fda = itg.fda+l;
fig.label(itg.lab) = ylabel('Temperature $T$ [$^{\circ}$C]'); itg.lab = itg.lab+1;
fig.label(itg.lab) = xlabel('Reactor Coordinate [-]'); itg.lab = itg.lab+1;

yyaxis right
fig.dat(itg.fda) = plot([linspace(0,1,N1) linspace(1,2,N2) linspace(2,3,N3)],...
                        full(root.w)*100,'r'); itg.fda = itg.fda+1;
hold on
fig.dat(itg.fda:itg.fda+l-1) = plot([linspace(0,1,N1) linspace(1,2,N2) linspace(2,3,N3)],...
                             full(nom_opt.w)*100,'r--'); itg.fda = itg.fda+l;
fig.label(itg.lab) = ylabel('Mass fraction NH$_3$ $w$ [wt.\%]'); itg.lab = itg.lab+1;
fig.label(itg.lab) = legend('Original $\mathbf{u}$','Optimized $\mathbf{u}$','Location','SE'); itg.lab = itg.lab+1;
legend BOXOFF 
itg.fig = itg.fig + 1;

% Redefinition of the labels
set(fig.label, 'Interpreter', 'Latex');
set(fig.label, 'Fontsize', fig.Font(1));
set(fig.axes,'TickLabelInterpreter', 'Latex');
set(fig.axes,'Fontsize', fig.Font(2));
end