

FileName = mfilename('fullpath');
[directory,~,~] = fileparts(FileName);
loadname.SOC = [directory(1:42) loadname.SOC(3:end)];
loadname.opt = [directory(1:42) loadname.opt(3:end)];
%%

clc, close all, clear fig

% Initialization of the system
% Load of the values

load(loadname.SOC)
load(loadname.opt)
% Initialization of the variables
if true
    opt.u = SS_opt{1}.u(:,16);
    opt.x = SS_opt{1}.x(:,16);
    opt.z = SS_opt{1}.z(:,16);
    opt.y = SS_opt{1}.y(:,16);
    opt.f = SS_opt{1}.f(:,16);
    SOC.u = SS_SOC{1,dec.H_comp}.u(:,16);
    SOC.x = SS_SOC{1,dec.H_comp}.x(:,16);
    SOC.z = SS_SOC{1,dec.H_comp}.z(:,16);
    SOC.f = SS_SOC{1,dec.H_comp}.f(:,16);
    opt.t = 0;
    opt.loss_int = 0;
    SOC.loss_int = 0;
    Ver1.opt.loss_int = 0;
    Ver2.opt.loss_int = 0;
    Ver3.opt.loss_int = 0;
    Ver4.opt.loss_int = 0;
    Ver1.SOC.loss_int = 0;
    Ver2.SOC.loss_int = 0;
    Ver3.SOC.loss_int = 0;
    Ver4.SOC.loss_int = 0;
end

for k = 2:t.N+1
    if any(int.sim(:,k)-int.sim(:,k-1))
        opt.t = [opt.t; t.samp(k); t.samp(k)];
        for l = 1:n.d
            if int.sim(l,k)-int.sim(l,k-1) ~= 0 && int.sim(l,k) ~=0
                opt.u = [opt.u opt.u(:,end) SS_opt{l}.u(:,SS_opt{l}.Var == int.sim(l,k)+d.nom(l))];
                opt.x = [opt.x opt.x(:,end) SS_opt{l}.x(:,SS_opt{l}.Var == int.sim(l,k)+d.nom(l))];
                opt.y = [opt.y opt.y(:,end) SS_opt{l}.y(:,SS_opt{l}.Var == int.sim(l,k)+d.nom(l))];
                opt.z = [opt.z opt.z(:,end) SS_opt{l}.z(:,SS_opt{l}.Var == int.sim(l,k)+d.nom(l))];
                opt.f = [opt.f opt.f(:,end) SS_opt{l}.f(:,SS_opt{l}.Var == int.sim(l,k)+d.nom(l))];

                SOC.u = [SOC.u SOC.u(:,end) SS_SOC{l,dec.H_comp}.u(:,SS_SOC{l,dec.H_comp}.Var == int.sim(l,k)+d.nom(l))];
                SOC.x = [SOC.x SOC.x(:,end) SS_SOC{l,dec.H_comp}.x(:,SS_SOC{l,dec.H_comp}.Var == int.sim(l,k)+d.nom(l))];
                SOC.z = [SOC.z SOC.z(:,end) SS_SOC{l,dec.H_comp}.z(:,SS_SOC{l,dec.H_comp}.Var == int.sim(l,k)+d.nom(l))];
                SOC.f = [SOC.f SOC.f(:,end) SS_SOC{l,dec.H_comp}.f(:,SS_SOC{l,dec.H_comp}.Var == int.sim(l,k)+d.nom(l))];
            end
        end
    end
    % Calculation and integration of the loss for Version 1
    if ~isempty(dec.Ver{1})
        Ver1.SOC.loss(k) = Ver1.EoR.sim(k)-SOC.f(end);
        Ver1.opt.loss(k) = Ver1.EoR.sim(k)-opt.f(end);
        Ver1.SOC.loss_int(k) = Ver1.SOC.loss_int(k-1)+t.int*(Ver1.SOC.loss(k)+Ver1.SOC.loss(k-1))/2;
        Ver1.opt.loss_int(k) = Ver1.opt.loss_int(k-1)+t.int*(Ver1.opt.loss(k)+Ver1.opt.loss(k-1))/2;
    end

    % Calculation and integration of the loss for Version 2
    if ~isempty(dec.Ver{2})
        Ver2.SOC.loss(k) = Ver2.EoR.sim(k)-SOC.f(end);
        Ver2.opt.loss(k) = Ver2.EoR.sim(k)-opt.f(end);
        Ver2.SOC.loss_int(k) = Ver2.SOC.loss_int(k-1)+t.int*(Ver2.SOC.loss(k)+Ver2.SOC.loss(k-1))/2;
        Ver2.opt.loss_int(k) = Ver2.opt.loss_int(k-1)+t.int*(Ver2.opt.loss(k)+Ver2.opt.loss(k-1))/2;
    end

    % Calculation and integration of the loss for Version 3
    if ~isempty(dec.Ver{3})
        Ver3.SOC.loss(k) = Ver3.EoR.sim(k)-SOC.f(end);
        Ver3.opt.loss(k) = Ver3.EoR.sim(k)-opt.f(end);
        Ver3.SOC.loss_int(k) = Ver3.SOC.loss_int(k-1)+t.int*(Ver3.SOC.loss(k)+Ver3.SOC.loss(k-1))/2;
        Ver3.opt.loss_int(k) = Ver3.opt.loss_int(k-1)+t.int*(Ver3.opt.loss(k)+Ver3.opt.loss(k-1))/2;
    end

    % Calculation and integration of the loss for self-optimizing control
    if ~isempty(dec.Ver{4})
        Ver4.SOC.loss(k) = Ver4.EoR.sim(k)-SOC.f(end);
        Ver4.opt.loss(k) = Ver4.EoR.sim(k)-opt.f(end);
        Ver4.SOC.loss_int(k) = Ver4.SOC.loss_int(k-1)+t.int*(Ver4.SOC.loss(k)+Ver4.SOC.loss(k-1))/2;
        Ver4.opt.loss_int(k) = Ver4.opt.loss_int(k-1)+t.int*(Ver4.opt.loss(k)+Ver4.opt.loss(k-1))/2;
    end
end

% Addition of the final time step
opt.t = [opt.t; t.samp(end)];
opt.u = [opt.u opt.u(:,end)];
opt.x = [opt.x opt.x(:,end)];
opt.y = [opt.y opt.y(:,end)];
opt.z = [opt.z opt.z(:,end)];
opt.f = [opt.f opt.f(:,end)];
opt.dy = opt.y-opt.y(:,1)*ones(1,size(opt.y,2));

SOC.u = [SOC.u SOC.u(:,end)];
SOC.x = [SOC.x SOC.x(:,end)];
SOC.z = [SOC.z SOC.z(:,end)];
SOC.f = [SOC.f SOC.f(:,end)];

% Definition of integer values
itg.fig = 1;
itg.lab = 1;
itg.fda = 1;
itg.axe = 1;
itg.ann = 1;
if t.tot <= 6000
    itg.xlab = 'Time $t$ [min]';
elseif  t.tot <= 3600*24
    itg.xlab = 'Time $t$ [h]';
else
    itg.xlab = 'Time $t$ [d]';
end
tStr = '\color[rgb]{1, 1, 1}titleString';

% Definition of figure and values
fig.dimw21 = 8.4;
fig.dimh21 = 12;
fig.dimw = 8.4;
fig.dimh = 6;
fig.Font = [12 10];
fig.label = [];
fig.ann   = [];

% Plot of the total extent of reaction with the optimal given as well
fig.in(itg.fig) = figure(itg.fig);
set(fig.in(itg.fig),'units','centimeters','position', [5+itg.fig,5+itg.fig, fig.dimw21, fig.dimh21]);
subplot(2,1,1)
fig.axes(itg.axe) = gca;
fig.axes(itg.axe).Position = [.15 .55 .775 .4];
fig.dat(itg.fda) = plot(opt.t,-opt.f*3.6,'Color',[0.93 0.69 0.13]); itg.fda = itg.fda+1;
hold on
if ~isempty(dec.Ver{1})
    fig.dat(itg.fda) = plot(t.samp(2:end),-Ver1.EoR.sim(2:end)*3.6,'r'); itg.fda = itg.fda+1;
end
if ~isempty(dec.Ver{2})
    fig.dat(itg.fda) = plot(t.samp(2:end),-Ver2.EoR.sim(2:end)*3.6,'b--'); itg.fda = itg.fda+1;
end
if ~isempty(dec.Ver{3})
    fig.dat(itg.fda) = plot(t.samp(2:end),-Ver3.EoR.sim(2:end)*3.6,'k'); itg.fda = itg.fda+1;
end
if ~isempty(dec.Ver{4})
    fig.dat(itg.fda) = plot(t.samp(2:end),-Ver4.EoR.sim(2:end)*3.6,'k--'); itg.fda = itg.fda+1;
end
fig.axes(itg.axe).XTickLabel = [];
fig.label(itg.lab) = ylabel('$J$ [t/h]'); itg.lab = itg.lab+1;
fig.ann(itg.ann) = annotation('textbox',[0 .925 .05 .05],'String','a)'); itg.ann = itg.ann+1;
% legend('$\xi_{opt,SS}$','T+SOC+ESC','T+SOC','T+ESC','T')
[~,fig.leg1] = legend('$\xi_{opt,SS}$','bounding','SS detection','without correction','T');
% set(fig.axes(itg.axe).Legend,'Location','northeast','FontSize',fig.Font(2),'Interpreter','latex','Box','off')
set(fig.axes(itg.axe).Legend,'Location','southeast','FontSize',fig.Font(2),'Interpreter','latex','Box','off')
fig.axes(itg.axe).XLim(2) = 10;
itg.axe = itg.axe+1;

% Plot of the loss with respect to the optimum
subplot(2,1,2)
fig.axes(itg.axe) = gca;
fig.axes(itg.axe).Position = [.15 .10 .775 .4];
hold on
if ~isempty(dec.Ver{1})
    fig.dat(itg.fda) = plot(t.samp(2:end),Ver1.opt.loss_int(2:end)/1e3,'r'); itg.fda = itg.fda+1;
end
if ~isempty(dec.Ver{2})
    fig.dat(itg.fda) = plot(t.samp(2:end),Ver2.opt.loss_int(2:end)/1e3,'b--'); itg.fda = itg.fda+1;
end
if ~isempty(dec.Ver{3})
    fig.dat(itg.fda) = plot(t.samp(2:end),Ver3.opt.loss_int(2:end)/1e3,'k'); itg.fda = itg.fda+1;
end
if ~isempty(dec.Ver{4})
    fig.dat(itg.fda) = plot(t.samp(2:end),Ver4.opt.loss_int(2:end)/1e3,'k--'); itg.fda = itg.fda+1;
end
fig.label(itg.lab) = xlabel(itg.xlab); itg.lab = itg.lab+1;
% fig.label(itg.lab) = ylabel('Integrated cost difference $\Delta J_{int}$ [t]'); itg.lab = itg.lab+1;
fig.label(itg.lab) = ylabel('Loss, $J_{int}$ [t]'); itg.lab = itg.lab+1;
fig.ann(itg.ann) = annotation('textbox',[0 .475 .05 .05],'String','b)'); itg.ann = itg.ann+1;
fig.axes(itg.axe).YLim = [0 10];
fig.axes(itg.axe).XLim(2) = 10;
itg.axe = itg.axe+1;
itg.fig = itg.fig + 1;  

set(fig.ann, 'Interpreter', 'Latex');
set(fig.ann(1:2), 'Fontsize', 10);
set(fig.ann(1:2), 'LineStyle', 'none');
set(fig.label, 'Fontsize', fig.Font(1));
set(fig.label, 'Interpreter', 'Latex');
set(fig.axes,'FontSize', fig.Font(2));
% set(fig.axes,'Legend.FontSize', fig.Font(2));
set(fig.axes,'Box','on');
set(fig.axes,'TickLabelInterpreter', 'Latex');
set(fig.axes,'LineWidth',.75);
set(fig.axes,'XMinorTick', 'on');
set(fig.axes,'YMinorTick', 'on');
set(fig.dat,'LineWidth',.75)
% fig.axes(2).FontSize = 6;
% fig.axes(end).Legend.FontSize = 8;

fig.axes(2).YTick = 0:10;
fig.axes(2).YTickLabel = {'0';'1';'2';'3';'4';'5';'6';'7';'8';'9';'10'};

%% Only to be used for evaluation of the detection
% for k = 5:2:12
%     fig.leg1(k).XData = [fig.leg1(k).XData(1) 0.15];
% end
% % fig.axes(1).Legend.Position = [0.85 0.4 0.2 0.5];

for k = 1:4
    fig.leg1(k).Interpreter = 'latex';
    fig.leg1(k).FontSize = fig.Font(2);
end