% Task  : Transformed MVs. SISO flash tank casestudy
% compare
% Date  : 5. March 2020
% Author: SB
% TODO  : 
% Var   : MVs=[zd]; CVs=[p]; States=[x]  DV=[].

clc;clear;%close all
% profile on
%% Plotting options
set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',16,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',16,...
'DefaultLineLineWidth',4,...
'DefaultLineMarkerSize',7.75,...
'DefaultStairLineWidth',2)
%% Colors
blue    = [0 0.4470 0.7410];
red     = [0.83 0 0.33];
green   = [0 0.5 0];
indigo  = [0.2930 0 0.5078];
lilac   = [0.7383    0.5898    0.7930];
hotpink = [0.8945    0.4102    0.5898];
cerulean = [0.1094    0.6641    0.8359];
orange  = [255,196,31]/255; 
bluegreen = [37,235,175]/255;
burgundy = [84,20,20]/255;
orange2  = [255,96+50,31]/255; 
blackish = [0,0,0];

%% Parameters
par.ca0 = 0.05;                    % [kmol/m3] Concentration of component A
par.V0 = 4;                                     % [m3] Hold up CSTR reactor
par.q10 = 1;                                        % [m3/min] Inlet stream 
par.q20 = 1;                                       % [m3/min] Outlet stream 
par.caf0 = 1;              % [kmol/m3] Concentration of component A in feed
par.k0 = 95;                         % [m3/kmol/min] Reaction rate constant   


par.A22 = -(par.q10/par.V0+2*par.k0*par.ca0);             % A22 = df2/dca|*
par.A21 = -(par.q10/par.V0^2)*(par.caf0-par.ca0);          % A21 = df2/dV|*
par.B21 = (par.caf0-par.ca0)/par.V0;                      % B21 = df2/dq1|*

par.taur = par.V0/par.q10;

par.v10 = 0;                                        % transformed MV 1 bias
par.v20 = 0;                                        % transformed MV 2 bias

%% calculate A = df/dx (x*,u*,d*)



%% Set up simulation
par.p   = [par.A22   %1
           0];       %2                                    % Parameter list
x0      = [par.V0;par.ca0];             % initial condition for integration
y0      = [par.V0;par.ca0];                                % initial output
u0      = [par.q10;par.q20];                            % physical MVs bias
d0      = [par.caf0;par.k0];                                     % DVs bias


%% variables which will be varied to make up the 30 cases
tsim    = 400;                                       % [s], simulation time

tstep_y1 = tsim*0.2;    % [s], step time for concentration setpoint changes
tstep_y2 = tsim*0.2;    % [s], step time for concentration setpoint changes
tstep_d = tsim*0.1;                       % [s], step time for disturbances
% tsim = 45;


cas =  1.0*par.ca0;            % setpoint change for concentration setpoint
Vs  =  1.0*par.V0;             % setpoint change for concentration setpoint
 
d1 = 1.0*par.caf0;                                           % change in DV
d2 = 1.0*par.k0;                                             % change in DV

ds = [d1;d2];

par.thetaca = 0.0*par.taur;
thetaList = [0,1,4]*par.taur;
tsim_list = [1,2,5]*tsim;
tstep_y1_list = [1,2,5]*tstep_y1;
tstep_y2_list = [1,2,5]*tstep_y2;

%input gain error
Ku = [2.0;2.0];

Vslist = [1.1 2 0.9 0.5]*Vs;
caslist = [1.1 2 0.9 0.5]*cas;



filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/CSTR/figures/bothGainError';
timeDelayFolder = {'/noTimeDelay','/1TimeDelay','/4timeDelay'};

IAEV_list = zeros(9,8);
IAEcA_list = zeros(9,8);

h = waitbar(0,'Running simulation');
%% Simulation. Get PID tuning, and run simulations
count = 1;
for i = 2%1:3
    for k = 1%1:2
        for j = 4%1:4
            if k == 1
                Vs = Vslist(j);
                cas = 1.0*par.ca0;
                name = ['/V',num2str(Vs/par.V0)];
            else
                Vs = 1.0*par.V0;
                cas = caslist(j);
                name = ['/cA',num2str(cas/par.ca0)];
            end
            tsim = tsim_list(i);
            tstep_y1 = tstep_y1_list(i);
            tstep_y2 = tstep_y2_list(i);
            
            par.thetaca = thetaList(i);
            [Kc,tauI,KI,par]=transformedcase(par);
            out1=sim('CSTR');
            [Kc,tauI,KI,par]=basecase(par);
            out2=sim('CSTR');
            [Kc,tauI,KI,par]=basecase_decouple2(par);
            out3=sim('CSTR_decouple2'); %this is the best decoupling
            
            IAEV_list((i-1)*3+1,(k-1)*4+j) = out1.IAE1(end);
            IAEV_list((i-1)*3+2,(k-1)*4+j) = out2.IAE1(end);
            IAEV_list((i-1)*3+3,(k-1)*4+j) = out3.IAE1(end);
            
            IAEcA_list((i-1)*3+1,(k-1)*4+j) = out1.IAE2(end);
            IAEcA_list((i-1)*3+2,(k-1)*4+j) = out2.IAE2(end);
            IAEcA_list((i-1)*3+3,(k-1)*4+j) = out3.IAE2(end);



            %% plotting
            fig=figure('visible','off');
            fig.Position = [100,100,600,450];

            subplot(221)
            plot(out1.time,out1.V,'color',blue);
            hold on
            plot(out2.time,out2.V,'--','color',hotpink);
            plot(out3.time,out3.V,'-.','color',green);
             plot(out1.time,out1.Vs,':','color',red);
            ylabel('y_1=V [m^{3}]')
            
%             ldg2=legend('Transformed','Basecase','Decoupled','ys');
% %             ldg2.Location='east';
%             ldg2.Location = 'northoutside';
%             ldg2.Box='off';

            subplot(222)
            plot(out1.time,out1.ca,'color',blue);
            hold on
            plot(out2.time,out2.ca,'--','color',hotpink);
            plot(out3.time,out3.ca,'-.','color',green);
            plot(out1.time,out1.cas,':','color',red);
            ylabel('y_2=c_{A} [kmol/m^{3}]')
            
            ldg2=legend('Transformed','Basecase','Decoupled','ys');
            ldg2.Location='northeast';
%             ldg2.Location = 'northoutside';
            ldg2.Box='off';

            subplot(223)
            plot(out1.time,out1.u(:,1),'color',blue)
            hold on 
            plot(out2.time,out2.u(:,1),'--','color',hotpink)
            plot(out3.time,out3.u(:,1),'-.','color',green)
            ylabel('u1=q_f [m^{3}/min]')
            xlabel('Time,[min]')

            subplot(224)
            plot(out1.time,out1.u(:,2),'color',blue)
            hold on 
            plot(out2.time,out2.u(:,2),'--','color',hotpink)
            plot(out3.time,out3.u(:,2),'-.','color',green)
            ylabel('u2=q [m^{3}/min]')
            xlabel('Time,[min]')

            fullname = [filepath,timeDelayFolder{i},name,'.eps'];
            saveas(fig,fullname,'epsc')
            
            waitbar(count/(3*8));
            count = count+1;
        end
    end
end
close(h)





