% This file is used for ADCHEM

% The number of disturbances can be up to
%   6
% If dec.loop is false, this number is reduced by 1. Additionally, the null
% space method can only be used, if less disturbances are specified. This
% is done automatically.
%
% The H matrices are defined as follows:
%   1   Pure local temperature controllers
%   2   Null space method utilizing all inlet and outlet temperatures as
%       well as the heat exchanger temperature
%   3   Exact local method utilizing all inlet and outlet temperatures as
%       well as the heat exchanger temperature
%   4   Local exact local method utilizing only inlet and outlet
%       temperatures of the respective reactor

% 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'])
import casadi.*

% Definition of decision variables
dec.slave = false;
dec.loss  = true;
dec.save  = false;
dec.plot  = true;
dec.Loop  = false;
dec.LoopL = true*~dec.Loop;
dec.newy  = true*~dec.Loop;
dec.dp    = dec.Loop*true;

% Definition of the chosen H and the measured disturbances
H_var = [4];                % Calculated and considered H matrices
dist_var = [1 2 3 4];       % Considered disturbances for SOC
if numel(dist_var) > 4
    H_var = H_var(H_var ~= 2);
end
if ~dec.loss
    dec.plot = false;
end

% Definition of the save and load name
if dec.loss
    savename = 'SS_opt_SOC_dist';
    savename_H = 'H_dist';
    loadname = 'SS_opt';
    for k = 1:max(dist_var)
        if any(dist_var == k)
            savename = [savename num2str(k)];
            savename_H = [savename_H num2str(k)];
        end
    end
    if dec.slave
        savename   = [savename '_MS'];
        savename_H = [savename_H '_MS'];
    end
    if dec.Loop
        savename   = [savename '_Loop'];
        savename_H = [savename_H '_Loop'];
        loadname   = [loadname '_Loop'];
    end
    if dec.LoopL
        savename   = [savename '_LoopL'];
        savename_H = [savename_H '_LoopL'];
        loadname   = [loadname '_LoopL'];
    end
    if dec.newy
        savename   = [savename '_newy'];
        savename_H = [savename_H '_newy'];
        loadname   = [loadname '_newy'];
    end
    savename   = [savename '.mat'];
    savename_H = [savename_H '.mat'];
    loadname   = [loadname '.mat'];
end    

% Definition of nominal disturbance variables
in.nom = zeros(7,1);
in.nom(1) = 3.0;            % H2/N2 ratio [-]
in.nom(3) = 200;            % Pressure [bar]
in.nom(4) = 250;            % Temperature [C]
in.nom(6) = 1;              % Catalyst activity 
if dec.Loop
    in.nom(2) = 16.2147;    % Feed flow (kg/s)
    in.nom(5) = 0.00;       % Mass fraction ammonia [-]
    in.nom(7) = 0.76;       % Recovery mass fraction of ammonia 
else
    in.nom(2) = 252/3.6;[61.7943997038480];%    % Feed flow (kg/s)
    in.nom(5) = 0.08;[0.0828623888813101];       % Mass fraction ammonia [-]
    in.nom(7) = 1;          % Recovery mass fraction of ammonia 
end    

% Load of the model
[func,dae,d,lb,ub,n,Cineq] = RecLoop_Ver6(1,dec,in.nom(2:7));

% Definition of the initial guesses and redefinition into x and z
if dec.Loop
    load([dirmain '\Data\Init_In_CL.mat'])
else
    load([dirmain '\Data\Init_In_OL.mat'])
end

% Definition of the variation
Var.n = zeros(numel(in.nom),1);
Var.n(1) = 41;41;
Var.n(2) = 41;
Var.n(3) = 41;
Var.n(4) = 21;
Var.n(5) = 61;
Var.lub = [-5 -20 -20 0.00 -0.4 -0.05;...
            5  20  20 0.02  0.2  0.05];
        
% Decision whether the ratio is an optimizing variable or not
dae.u = dae.p(1:n.u);
n.d   = length(d);
n.u   = length(dae.u);
nom_opt.d = zeros(n.d,1);

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

% Definition of the number of measurements in the QP approach
n.y_used = 2;
n.y_var  = 11;

% Definition of indices used in the solution
ind.T    = [2*N+1 1:N/3       2*N+2 N/3+1:2*N/3   2*N+3  2*N/3+1:N]; % 3*N+4 
ind.wNH3 = [2*N+6 N+1:4*N/3   2*N+7 4*N/3+1:5*N/3 2*N+8  5*N/3+1:2*N];
ind.wN2  = [];

ind.Rea = (2*N+9):(2*N+12);
ind.Rec = (2*N+13):(2*N+16);
ind.W   = (2*N+17);

ind.x = 1:N;
ind.y = [2*N+1 N/3 2*N+2 N*2/3 2*N+3 N];
ind.z = (N+1):(2*N+17);

ind.u = (2*N+17+1):(2*N+17+n.u);

z0 = [z0([1:30 61:64]); 420; z0([65:67 71:74 76:79 81])];

if n.y_var == 11
    ind.y = ind.T;
elseif n.y_var == 22
    ind.y = [ind.T(1:33) ind.wNH3(1:33)];
end
n.y = length(ind.y);
    
%% Part a) Root finding for the modified 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.u,p0(1:n.u));
J_root = substitute(dae.quad    ,d,zeros(n.d,1));
J_root = substitute(J_root,dae.u,p0(1:n.u));
cost = Function('cost',{w_root},{J_root});

% 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);

% 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.wNH3 = full(root.comb.o0(ind.wNH3));
root.wN2  = full(root.comb.o0(ind.wN2));

root.Rea = full(root.comb.o0(ind.Rea));
root.Rec = full(root.comb.o0(ind.Rec));
root.W   = full(root.comb.o0(ind.W));

root.f   = -full(cost(root.comb.o0));
root.u  = p0(1:n.u);

x0 = root.x;
z0 = root.z;

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; dae.u(4); Cineq];
lbg  = [zeros(size(gadd,1)-1,1); -1];
ubg  = [zeros(size(gadd))];

% Substitution of the inlet values
g = substitute(gadd,d,     zeros(n.d,1));
J = substitute(dae.quad,d, zeros(n.d,1));

% Definition of the NLP function
nlp = struct('x',w_opt,'f',J,'g',g);
opts = struct('warn_initial_bounds',false	...
          ,'print_time',false               ...
          ,'ipopt',struct('print_level',3)  ...
          );
solver = nlpsol('solver','ipopt',nlp,opts);
nom_opt.comb = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',lbg,'ubg',ubg);

% Reassignment of variables
nom_opt.x = full(nom_opt.comb.x(ind.x));
nom_opt.z = full(nom_opt.comb.x(ind.z));
nom_opt.y = full(nom_opt.comb.x(ind.y));
    
nom_opt.T    = full(nom_opt.comb.x(ind.T));
nom_opt.wNH3 = full(nom_opt.comb.x(ind.wNH3));
nom_opt.wN2  = full(nom_opt.comb.x(ind.wN2));

nom_opt.Rea = full(nom_opt.comb.x(ind.Rea));
nom_opt.Rec = full(nom_opt.comb.x(ind.Rec));
nom_opt.W   = full(nom_opt.comb.x(ind.W));

nom_opt.u = full(nom_opt.comb.x(ind.u));
nom_opt.p = [full(nom_opt.u); zeros(4,1)];
nom_opt.f = full(nom_opt.comb.f);

fprintf('\n\n\nThe steady-state optimization is completed. \n\n\n')

%% Part c) Calculation of the sensitivity matrix F

% Definition of the variation
pert_opt.dist = [-1 -1 -1 -0.001 -0.05 -0.005]*-1;
n.pd = numel(dist_var);

% Predefinition of the solution vectors
pert_opt.x = zeros(n.x,n.pd);
pert_opt.z = zeros(n.z,n.pd);
pert_opt.y = zeros(n.y,n.pd);

pert_opt.u = zeros(n.u,n.pd);
pert_opt.f = zeros(1,n.pd);

% Definition of the varying vector
in.sens = diag(pert_opt.dist);
pert_opt.count = 1:n.d;
k = 1;

% Solving of the system in the case of disturbances
for j = pert_opt.count(dist_var)
    
    % Substitution of the inlet values
    g = substitute(gadd,d,in.sens(:,j));
    J = substitute(dae.quad,d,in.sens(:,j));

    % Definition of the NLP function
    nlp = struct('x',w_opt,'f',J,'g',g);
    solver = nlpsol('solver','ipopt',nlp,opts);
    pert_opt.comb = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',lbg,'ubg',ubg);

    % Reassignment of variables
    pert_opt.x(:,k) = full(pert_opt.comb.x(ind.x));
    pert_opt.y(:,k) = full(pert_opt.comb.x(ind.y));
    pert_opt.z(:,k) = full(pert_opt.comb.x(ind.z));
    
    pert_opt.u(:,k) = full(pert_opt.comb.x(ind.u));
    pert_opt.f(:,k) = -full(pert_opt.comb.f);
    k = k+1;
end

% Calculation of the disturbance
pert_opt.dd = ones(n.y,1)*pert_opt.dist(dist_var);

% Calculation of the sensitivity matrix
F = (pert_opt.y - nom_opt.y*ones(1,n.pd))./pert_opt.dd;

%% Part d) Calculation of the plant gain matrix Gy_n by manipulating one DOF

% Predefinition of the solution vectors
pert_con_n.x = zeros(n.x,n.u-1);
pert_con_n.y = zeros(n.y,n.u-1);
pert_con_n.z = zeros(n.z,n.u-1);
pert_con_n.u = zeros(n.u,n.u-1);
clc
w0 = full(nom_opt.comb.x);
aux = 1:n.u-1;

% Solving of the system in the case of disturbances
for j = 1:n.u-1
    
    % Switch in betwee having a low level slave controller or not
    if dec.slave
        g = vertcat(gadd, (dae.z(N+j) - (nom_opt.z(N+j)+1)));
        g = vertcat(g,    (dae.z(N+aux(aux~=j)) - (nom_opt.z(N++aux(aux~=j)))));
    else
        g = vertcat(gadd, (dae.p(j) - nom_opt.u(j)*0.9999999));
        g = vertcat(g, (dae.p(aux(aux~=j)) - nom_opt.u(aux(aux~=j))));
    end
    
    % Substitution of the inlet values
    g = substitute(g,d,zeros(n.d,1));
    J = substitute(dae.quad,d,zeros(n.d,1));

    % Definition of the NLP function
    nlp = struct('x',w_opt,'f',J,'g',g);
    solver = nlpsol('solver','ipopt',nlp,opts);
    pert_con_n.comb = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',[lbg;zeros(n.u-1,1)],'ubg',[ubg;zeros(n.u-1,1)]);

    % Reassignment of variables
    pert_con_n.x(:,j) = full(pert_con_n.comb.x(ind.x));
    pert_con_n.y(:,j) = full(pert_con_n.comb.x(ind.y));
    pert_con_n.z(:,j) = full(pert_con_n.comb.x(ind.z));
    pert_con_n.u(:,j) = full(pert_con_n.comb.x(ind.u));
    pert_con_n.f(:,j) = full(pert_con_n.comb.f);
end

% Calculation of the setpoint change
if dec.slave
    pert_con_n.du = pert_con_n.y(2:2:end,:)-nom_opt.y(2:2:end)*ones(1,n.u-1);
else
    pert_con_n.du = pert_con_n.u-nom_opt.u*ones(1,n.u-1);
end
pert_con_n.du = ones(n.y,1)*diag(pert_con_n.du)';

% Calculation of the plant gain matrix
Gy_n = (pert_con_n.y - nom_opt.y*ones(1,n.u-1))./pert_con_n.du;

% Calculation of the relative gain array
RGA.nom = Gy_n.*pinv(Gy_n)';


%% Part e) Alternative method for calculation of Gy and F
% Redefinition of the constraints, cost function and variables
g = [dae.ode; dae.alg];
J = dae.quad;
u = dae.p(1:n.u);
z = [dae.x; dae.z];

% Definition of the rootfinder function
rf_function = casadi.Function('rf_function',{z,u,d},...
                    {g,J,z},{'z','u','d'},{'g','J','z'});
rf  = casadi.rootfinder('rf','newton',rf_function);

% Definition of the casadi functions
fun.Gy  = rf.jacobian('u','z');
fun.Gd  = rf.jacobian('d','z');
fun.Juu = rf.hessian('u','J');
fun.Jud = jacobian(rf.jacobian('u','J'),'d');
fun.Jdu = jacobian(rf.jacobian('d','J'),'u');

% Calculation of the matrices
nom_opt.w = [nom_opt.x; nom_opt.z];
ad.Gy  = full(fun.Gy( nom_opt.w, nom_opt.u, nom_opt.d));
ad.Gd  = full(fun.Gd( nom_opt.w, nom_opt.u, nom_opt.d));
ad.Juu = full(fun.Juu(nom_opt.w, nom_opt.u, nom_opt.d));
ad.Jud = full(fun.Jud(nom_opt.w, nom_opt.u, nom_opt.d));
ad.Jdu = full(fun.Jdu(nom_opt.w, nom_opt.u, nom_opt.d));
ad.Gy    = ad.Gy(:,1:3);
ad.Gd    = ad.Gd(:,1:4);
ad.Juu   = ad.Juu(1:3,1:3);
ad.Jud   = ad.Jud(1:3,1:4);
ad.Gyout = ad.Gy([65 60 69],1:3);
ad.F = ad.Gd-ad.Gy/ad.Juu*ad.Jud;

%% Part f) Application of the exact local method on the system

% Definition of the scaling matrices
if dec.Loop
    noise.dist = [5 20 20 0.01 0.3000 0.1000];%mean(abs(Var.lub));
else
    noise.dist = [5 20 20 0.01 0.3000 0.1000];%mean(abs(Var.lub))
end

noise.dist = noise.dist(dist_var);
Wd = diag(noise.dist);

% F = ad.F(ind.y,:);
% Gy_n = ad.Gy(ind.y,:);
noise.meas = ones(1,3*n.y)*4;
Wn = diag(noise.meas);

ind.y_used = ind.y([1 11 12 22 23 33]);
ind.y_red = [1 11 12 22 23 33];

% Definition of the measurement indices
ind.m1 = 1:2;
ind.m2 = 3:4;
ind.m3 = 5:6;
ind.m4 = 2:7;

% Self-optimizing control based on direct selection
H{1} = zeros(3,7);
H{1}(1,3) = 1;
H{1}(2,5) = 1;
H{1}(3,6) = 1;

% SOC using the nullspace method
H{2} = null(F')';

% SOC using the Exact local method for all inputs
Y  = [F(ind.y_red,:)*Wd Wn(ind.y_red,ind.y_red)];
H{3} = ((Y*Y')\Gy_n(ind.y_red,1:3))';

% SOC using the exact local method for different inputs
Y1 = [F(ind.y_red(ind.m1),:)*Wd Wn(ind.y_red(ind.m1),ind.y_red(ind.m1))];
Y2 = [F(ind.y_red(ind.m2),:)*Wd Wn(ind.y_red(ind.m2),ind.y_red(ind.m2))];
Y3 = [F(ind.y_red(ind.m3),:)*Wd Wn(ind.y_red(ind.m3),ind.y_red(ind.m3))];
H1 = ((Y1*Y1')\Gy_n(ind.y_red(ind.m1),1))';   
H2 = ((Y2*Y2')\Gy_n(ind.y_red(ind.m2),2))';
H3 = ((Y3*Y3')\Gy_n(ind.y_red(ind.m3),3))';
H{4}  = [blkdiag(H1,H2,H3)];

% Scaling of the H matrix
for l = H_var
    H{l}(1,:) = H{l}(1,:)/max(abs(H{l}(1,:)));
    H{l}(2,:) = H{l}(2,:)/max(abs(H{l}(2,:)));
    H{l}(3,:) = H{l}(3,:)/max(abs(H{l}(3,:)));
end

% Calculation of the RGA for the control structures
for j = H_var
    RGA.H{j} = (H{j}*Gy_n(ind.y_red,:)).*pinv(H{j}*Gy_n(ind.y_red,:))';
end

% Save of the results
if dec.save
    save([dirmain '\Data\' savename_H],'H','noise')
end


%% Part g) Root finding for the new system

if dec.loss
% Clear of aux
clear aux

% Redefinition of the variables
if dec.LoopL
    
    Var.lub(:,1) = [-5 5]';
    Var.lub(1,4) = 0;
    Var.n(4)     = 21;
    Var.n(7)     = 0;
    
    dec.Loop = dec.LoopL;
    dec.dp   = true;
    
    % Definition of nominal disturbance variables
    in.nom = zeros(7,1);
    in.nom(1) = 3.00;            % H2/N2 ratio [-]
    in.nom(3) = 200;            % Pressure [bar]
    in.nom(4) = 250;            % Temperature [C]
    in.nom(6) = 1;              % Catalyst activity 
    if dec.Loop
        in.nom(2) = 16.2147;    % Feed flow (kg/s)
        in.nom(5) = 0.00;       % Mass fraction ammonia [-]
        in.nom(7) = 0.76;       % Recovery mass fraction of ammonia 
    else
        in.nom(2) = 252/3.6;    % Feed flow (kg/s)
        in.nom(5) = 0.08;       % Mass fraction ammonia [-]
        in.nom(7) = 1;          % Recovery mass fraction of ammonia 
    end    
    
    load([dirmain '\Data\Init_In_CL.mat'])
    [func,dae,d,lb,ub,n,Cineq] = RecLoop_Ver6(1,dec,in.nom(2:7));
    dae.u = dae.p(1:n.u);
    n.d   = length(d);
    n.u   = length(dae.u);

    z0 = [z0([1:30 61:64]); 420; z0([65:67 71:74 76:79 81])];
    
    g = vertcat(dae.ode,dae.alg);

    % 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; dae.u(4); Cineq];
    lbg  = [zeros(size(gadd,1)-1,1); -1];
    ubg  = [zeros(size(gadd))];

    % Substitution of the inlet values
    g = substitute(gadd,d,     zeros(n.d,1));
    J = substitute(dae.quad,d, zeros(n.d,1));

    % Definition of the NLP function
    nlp = struct('x',w_opt,'f',J,'g',g);
    solver = nlpsol('solver','ipopt',nlp,opts);
    nom_opt.comb = solver('x0',w0,'lbx',lbw,'ubx',ubw,'lbg',lbg,'ubg',ubg);
    if dec.newy
        nom_opt.y = full(nom_opt.comb.x(ind.y));
    end
end

% Redefinition of the initial conditions
w0 = full(nom_opt.comb.x);
    
% Combination of the variables
w = [dae.x; dae.z; dae.u];
Var.n(5) = 0;

% Definition of the number of different H matrices
n.H = numel(H);

% Predefinition of variables and load of the optimal values
SS_SOC = cell(n.d,n.H);
in.var = zeros(n.d,1);
loadname = 'SS_opt';
if dec.Loop
    loadname   = [loadname '_Loop'];
end
loadname   = [loadname '.mat'];
load([dirmain '\Data\' loadname])

for l = H_var

% if dec.Loop && ~dec.LoopL
%     if l == 1
%         Var.n(2) = 99;
%         Var.lub(1,2) = -48;
%     elseif l == 3
%         Var.n(2) = 84;
%         Var.lub(1,2) = -33;
%         Var.n(3) = 93;
%         Var.lub(1,3) = -42;
%     elseif l == 4
%         Var.n(2) = 77;
%         Var.lub(1,2) = -26;
%         Var.n(3) = 91;
%         Var.lub(1,3) = -40;
%     end
% end

% Combination of the equality constraints to keep the setpoint constant
g_org = vertcat(dae.ode,dae.alg,H{l}*w(ind.y_used)-H{l}*nom_opt.y(ind.y_red),dae.u(4)-0);

aux2 = 1:n.d;

for j = aux2(Var.n ~= 0)
    in.SOC{j,1} = linspace(Var.lub(1,j),Var.lub(2,j),Var.n(j));
    for k = 1:Var.n(j)
        % Substitution of the inlet values
        in.dist = in.var;
        in.dist(j) = in.dist(j)+in.SOC{j}(k);
        g = substitute(g_org,d,in.dist);
        J_root = substitute(dae.quad    ,d,in.dist);
        cost = Function('cost',{w},{J_root});

        % Definition and solving of the system
        g_fun = Function('g_fun',{w},{g});
        G = rootfinder('G','newton',g_fun);
        aux.comb = G('i0',w0);
        aux.x(:,k) = full(aux.comb.o0(ind.x));
        aux.y(:,k) = full(aux.comb.o0(ind.y));
        aux.z(:,k) = full(aux.comb.o0(ind.z));
        aux.u(:,k) = full(aux.comb.o0(ind.u));
        aux.f(k)   = -full(cost(aux.comb.o0));
%         w0 = full(aux.comb.o0);
    end
    w0  = full(nom_opt.comb.x);
    
    % Combination into a structure
    SS_SOC{j,l} = struct('Var',in.SOC{j}+in.nom(j+1),'z',aux.z,...
                         'u',aux.u,'x',aux.x,'y',aux.y,'f',aux.f);
                     
    % Calculation of the loss
    if numel(SS_opt{j}.Var) == numel(SS_SOC{j,l}.Var)
    if sum(abs(diff([SS_opt{j}.Var;SS_SOC{j,l}.Var]))) < 1e-6
        SS_SOC{j,l}.Var = SS_opt{j}.Var;
    end
    end
    SS_SOC{j,l}.loss = abs((SS_opt{j}.f(end,ismember(SS_opt{j}.Var,SS_SOC{j,l}.Var))-SS_SOC{j,l}.f(end,:))'...
                          ./SS_opt{j}.f(end,ismember(SS_opt{j}.Var,SS_SOC{j,l}.Var))')*100;
%     SS_SOC{j,l}.loss = (SS_opt{j}.f(end,:)-SS_SOC{j,l}.f(end,:))'...
%                         ./SS_opt{j}.f(end,:)'*100;

    clear aux
    fprintf('The steady-state simulation for the disturbance %1i and the H matrix %1i is finished. \n',j,l)
end
end

% Save of the results
if dec.save
    save([dirmain '\Data\' savename],'SS_SOC','H','noise')
end
end


%% Part g) Plotting of the results

if dec.plot
    % Definition of integer values
    itg.fig = 1;
    itg.lab = 1;
    itg.arr = 1;
    itg.sub = 1;

    % Definition of figure and values
    fig.xlab = {'Ratio H$_2$/N$_2$ [-]','Flowrate [kg/s]','Pressure [bar]','Temperature [C]',...
                'Mass fraction [-]','Catalyst activity [-]','Recovery NH$_3$ [-]'};
    fig.xlim = [];
    fig.d     = [0.1 5 25 25 0.02 0.2];
    for k = aux2(Var.n ~= 0)
        fig.xlim = [fig.xlim;
                    floor(round(SS_opt{k}.Var(1)/fig.d(k+1)))*fig.d(k+1) ...
                    fig.d(k+1) ...
                    ceil(round(SS_opt{k}.Var(end)/fig.d(k+1)))*fig.d(k+1)];
    end
    fig.dimw  = 12;
    fig.dimh  = 7;
    fig.label = [];
    fig.Font  = [10 8];
    fig.ann   = [];
    fig.sub   = [];

    % Definition of the colours
    rk = {'k','k--','r','r--'};
    H_cell = {'T','NSM','ELM$_1$','ELM$_2$'};

    % Initialization of the figure
    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]);

    % Plotting of the loss with flowrate change
    for k = aux2(Var.n ~= 0)
        fig.sub(k) = subplot(2,ceil(n.d/2),k);
        fig.axes(k) = gca;
        xlim(fig.xlim(k,[1 3]))
%         ylim([0 3])
        hold on
        for l = H_var
            plot(SS_SOC{k,l}.Var,SS_SOC{k,l}.loss,rk{l})
        end
        fig.axes(k).XTick = fig.xlim(k,1):fig.xlim(k,2):fig.xlim(k,3);
%         fig.axes(k).YTick = 0:1:3;
        fig.label(itg.lab) = xlabel(fig.xlab{k+1}); itg.lab = itg.lab+1;
        fig.axes(k).LineWidth = 1;
        if mod(k-1,n.d/2) == 0
            fig.label(itg.lab) = ylabel('Loss [\%]'); itg.lab = itg.lab+1;
        end
    end

    % Modificiation of the plotting properties
    set(fig.ann, 'Interpreter', 'Latex');
    set(fig.ann, 'LineStyle', 'none');
    set(fig.label, 'Interpreter', 'Latex');
    set(fig.axes,'TickLabelInterpreter', 'Latex');
    set(fig.axes,'XMinorTick', 'on');
    set(fig.ann, 'Fontsize', fig.Font(1));
    set(fig.label, 'Fontsize', fig.Font(1));
    set(fig.label(3), 'Fontsize', fig.Font(2));
    set(fig.axes,'Fontsize', fig.Font(2));
end