% Task  : Transformed MVs. SISO mixing tank
% compare
% Date  : 14 January 2020
% Author: SB
% TODO  : 
% Var   : MVs=[q1;q2]; CVs=[M;c]; DV=[q;c1;c2].  

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.M0   = 10;                                          % [m3], tank holdup

par.q10  = 1;                                         % [m3/s] inlet flow 1
par.q20  = 9;                                         % [m3/s] inlet flow 2

par.taur  = par.M0/(par.q10+par.q20);             % [s] tank residence time

par.c0   = 1.1;                          % [], nominal outlet concentration
par.c10   = 2;                % [kg/m3], nominal inlet flow 1 concentration
par.c20   = 1;                % [kg/m3], nominal inlet flow 2 concentration

par.A   = -(par.q10+par.q20);

par.v_0 = 0;                             

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

tstep_c = tsim*0.1;     % [s], step time for concentration setpoint changes
tstep_d = tsim*0.1;                       % [s], step time for disturbances


par.gain_error  =  1.0;                      % input error gain proportions
par.disturbance_error = eye(4);       % disturbance gain error, proportions
par.disturbance_error(1,1) = 1.2;
par.disturbance_error(2,2) = 1.2;
par.disturbance_error(3,3) = 1.2;


cs1 =  1.1;                    % setpoint change for concentration setpoint
% cs1 =  1.7;

% d1 = 1.0*par.q20;                                          % change in DV  
% d1 = 2.0*par.q20;                                          % change in DV
d1 = 0.5*par.q20;                                            % change in DV
d2 = 1.0*par.c10;
% d2 = 1.5*par.c10;
d3 = 1.0*par.c20;
% d3 = 0.5*par.c20;
d4 = 1.0*par.M0;
ds = [d1;d2;d3;d4];


par.theta  = 3.0*par.taur;                                          % delay
% par.theta  = 0.3*par.taur;                                        % delay
% par.theta  = 0.0;                                                 % delay  


%% Set up simulation
par.p  = [par.A;0];
x0      = par.c0;                       % initial condition for integration
u0      = par.q10;                                      % physical MVs bias
d0      = [par.q20;par.c10;par.c20;par.M0];                      % DVs bias
dv      = 4;                      % step in inputs for open loop simulation


%% Simulation. Get PID tuning, and run simulations
[Kc,tauI,KI,par]=case3(par);
out=sim('mixingTank');
IAE = out.IAE(end)


%% plotting
fig=figure('visible','on');
subplot(3,1,1)
p1=plot(out.time,out.cm,'color',blue);
hold on
p2=plot(out.time,out.c_setpoint,':','color',red);
% ylim([1,1.2])
ylabel('y=c_{m}, [kmol/m^3]')
% xline(tstep_d,':','color',green,'LineWidth',2);
legend([p1,p2],'y = c_{m}','c_s','Location','southeast')


subplot(3,1,2)
plot(out.time,out.u)
ylabel('u=q_1 [m^3/s]')
xlabel('Time,[s]')
% xline(tstep_d,':','color',green,'LineWidth',2);

subplot(3,1,3)
plot(out.time,out.v)


%% have all the cases written down here as help functions

function [Kc,tauI,KI,par]=case0(par)
k    = (par.c10-par.c0)/(par.q10+par.q20);       % process gain from v to c
tau  = par.M0/(par.q10+par.q20);          % time constant for concentration
par.p(2) = 0;

par.v_0 = par.q10;                          % transformed input bias case 1


tauC = par.theta;         % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end


% PI-controller for concentration
Kc   = 1/k*tau/(tauC+par.theta);                   % first-order plus delay
tauI = min(tau,4*(tauC+par.theta));   
KI   = Kc/tauI;
end

function [Kc,tauI,KI,par]=case1(par)
k    = 1/(par.q10+par.q20);                      % process gain from v to c
tau  = par.M0/(par.q10+par.q20);          % time constant for concentration
par.p(2) = 1;


u_nom = par.q10;
par.v_0 = u_nom*par.c10 + par.q20*par.c20;  % transformed input bias case 2



tauC = par.theta;         % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end


% PI-controller for concentration
Kc   = 1/k*tau/(tauC+par.theta);                   % first-order plus delay
tauI = min(tau,4*(tauC+par.theta)); 
KI   = Kc/tauI;
end

function [Kc,tauI,KI,par]=case2(par)
k    = 1;                                        % process gain from v to c
tau  = par.M0/(par.q10+par.q20);          % time constant for concentration
par.p(2) = 2;


u_nom = par.q10;
par.v_0 = (u_nom*par.c10+par.q20*par.c20)/(u_nom+par.q20);
                                            % transformed input bias case 2

tauC = par.theta;         % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end


% PI-controller for concentration
Kc   = 1/k*tau/(tauC+par.theta);                   % first-order plus delay
tauI = min(tau,4*(tauC+par.theta));   
KI   = Kc/tauI;
end

function [Kc,tauI,KI,par]=case3(par)
k    = -1/par.A;                                 % process gain from v to c
tau  = -par.M0/par.A;                     % time constant for concentration
par.p(2) = 3;


u_nom = par.q10;
par.v_0 = u_nom*(par.c10-par.c0) + par.q20*(par.c20-par.c0) - par.c0*par.A;     
                                            % transformed input bias case 4

tauC = 1.0*par.theta;     % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end


% PI-controller for concentration
Kc   = 1/k*tau/(tauC+1*par.theta);               % first-order plus delay
tauI = min(tau,4*(tauC+1*par.theta));   
KI   = Kc/tauI;
end

function [Kc,tauI,KI,par]=case4(par)
k    = 1;                                        % process gain from v to c
tau  = -par.M0/par.A;                     % time constant for concentration
par.p(2) = 4;


u_nom = par.q10;
par.v_0 = -(u_nom*(par.c10-par.c0)+par.q20*(par.c20-par.c0)-par.c0*par.A)/par.A; 
                                            % transformed input bias case 4

tauC = 1.0*par.theta;     % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end


% PI-controller for concentration
Kc   = 1/k*tau/(tauC+par.theta);                   % first-order plus delay
tauI = min(tau,4*(tauC+par.theta));   
KI   = Kc/tauI;
end

function [Kc,tauI,KI,par]=case5(par)
k    = (par.q20*(par.c10-par.c0))/(par.q10+par.q20);
                                                 % process gain from v to c
tau  = par.M0/(par.q10+par.q20);          % time constant for concentration
par.p(2) = 5;


u_nom = par.q10;
par.v_0 = u_nom/par.q20;                    % transformed input bias case 5


tauC = par.theta;         % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end



% PI-controller for concentration 
Kc   = 1/k*tau/(tauC+par.theta);                   % first-order plus delay
tauI = min(tau,4*(tauC+par.theta));
KI   = Kc/tauI;
end

function [Kc,tauI,KI,par]=case6(par,P)
k    = 1/par.M0;                                 % process gain from v to c

par.p(2) = 6;


u_nom = par.q10;
par.v_0 = u_nom*(par.c10-par.c0) + par.q20*(par.c20-par.c0);
                                            % transformed input bias case 6

tauC = par.theta;         % desired closed loop time constant for flow loop
if tauC==0
    tauC = 0.3*par.taur;
end

% PI-controller for concentration
Kc   = 1/k/(tauC+par.theta);                        % Integrated plus delay
tauI = 4*(tauC+par.theta); 
KI   = Kc/tauI;
if P
    KI   = 0;
end
end
