% 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',2,...
'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];

%% 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*1;                                    % [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.2;                       % [s], step time for disturbances
% tsim = 45;


cas =  2.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 = 1.0*par.taur;

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



%% Simulation. Get PID tuning, and run simulations
tic
[Kc,tauI,KI,par]=transformedcase(par);
out1=sim('CSTR');
toc
tic
[Kc,tauI,KI,par]=basecase(par);
out2=sim('CSTR');
toc
tic
[Kc,tauI,KI,par]=basecase_decouple2(par);
out3=sim('CSTR_decouple2'); %this is the best decoupling
toc




%% plotting
fig=figure('visible','on');
subplot(321)
p1=plot(out1.time,out1.V,'color',blue);
hold on
p2=plot(out2.time,out2.V,'-.','color',hotpink);
p3=plot(out3.time,out3.V,':','color',green);
% p4=plot(out4.time,out4.V,'-.','color',indigo);
% hold on
% p2=plot(out.time,out.cbs,':','color',red);
% ylim([1,1.2])
ylabel('x=V')
% xline(tstep_d,':','color',green,'LineWidth',2);

subplot(322)
p5=plot(out1.time,out1.ca,'color',blue);
hold on
p6=plot(out2.time,out2.ca,'-.','color',hotpink);
p7=plot(out3.time,out3.ca,':','color',green);
% p8=plot(out4.time,out4.ca,'-.','color',indigo);
ylabel('x=c_{a}')

subplot(323)
plot(out1.time,out1.u(:,1))
hold on 
plot(out2.time,out2.u(:,1),'-.','color',hotpink)
plot(out3.time,out3.u(:,1),':','color',green)
% plot(out4.time,out4.u(:,1),'-.','color',indigo)
ylabel('u1=q1 [-]')
xlabel('Time,[s]')
subplot(324)
plot(out1.time,out1.u(:,2))
hold on 
plot(out2.time,out2.u(:,2),'-.','color',hotpink)
plot(out3.time,out3.u(:,2),':','color',green)
% plot(out4.time,out4.u(:,2),'-.','color',indigo)
ylabel('u2=q2 [-]')
xlabel('Time,[s]')


subplot(325)
plot(out1.time,out1.v(:,1))
ylabel('v1 [-]')
subplot(326)
plot(out1.time,out1.v(:,2))
ylabel('v2 [-]')

% subplot(427)
% plot(out3.time,out3.u1)
% subplot(428)
% plot(out3.time,out3.u2)




