% Task  : Transformed MVs. SISO mixing tank 30-cases, all cases, maybe
% compare
% Date  : 19 February 2020
% Author: SB
% TODO  : 
% Var   : MVs=[q1;q2]; CVs=[M;c]; DV=[q;c1;c2].  

% clc;clear;close all
%% 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];

% colorList = {blue,green,indigo,lilac,{},hotpink,cerulean,burgundy};
colorList = {blue,green,indigo,orange,{},hotpink,cerulean,blackish};
%% 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_list    = [100;20;10];                          % [s], simulation time
% tsim_list    = [1000;20;10];                      % [s], simulation time
% tsim_list    = [500;20;10];                      % [s], simulation time

tstep_c = 1200;         % [s], step time for concentration setpoint changes
tstep_d_list = tsim_list*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;%1.2

cs1 = 1.1;                     % setpoint change for concentration setpoint

d1_list = [2*par.q20;0.5*par.q20];                           % change in DV

ds=[par.q20;par.c10;par.c20;par.M0];

 

theta_list  = [3*par.taur;0.3*par.taur;0];                          % delay
% theta  = 0.3*taur;                                                % delay
% theta  = 0;



%% Set up simulation
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      = 1;                      % step in inputs for open loop simulation


%% Simulation. Get PID tuning, and run simulations
IAE = zeros(6,8);
l=1;
count=1;
progressbar = waitbar(0,'Progress');
filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/Mixingtank30-cases/figures/disturbancegainerror/disturbancerejection';  %tauc>theta,v5=p_control
for k =1%1:3
    for j = 2%1:2
        for i = 3%[0:3,5:7]
            if(i==3)
                f=1;
            end
            par.theta = theta_list(k);
            tsim = tsim_list(k);
            tstep_d=tstep_d_list(k);
            ds(1) = d1_list(j);
            [Kc,tauI,KI,par]=tuning_case(par,i);
            name = ['/',num2str(ds(1)-par.q20),'_q2_',num2str(par.theta/par.taur),'_delay'];
            if i<6
                casenr = ['/v',num2str(i)];
            elseif i==6
                casenr = '/v6_PI';
            else
                casenr = '/v6_P';
            end
                
            out=sim('mixingTank');
            
            waitbar(count/48);
            IAE(l,(i+1)) = out.IAE(end);
            count=count+1;

%% plotting 1
        if (i==0)
            fig2=figure('visible','off');
            fig2.Position = [100,100,450,450];
           %fig2.Position = [left,bottom,width,height];
        end
        if (i<-6 || k>10)
            set(0,'Currentfigure',fig2);
            subplot(2,1,1)
            plot(out.time,out.cm,'color',colorList{i+1});
            hold on
            ylim([0.6,1.5])
            ylabel('y=c_{m}, [kmol/m^3]')
            subplot(2,1,2)
            plot(out.time,out.u,'color',colorList{i+1})
            hold on
            ylim([-1,8])
            ylabel('u=q_1 [m^3/s]')
            xlabel('Time,[s]')
        end
        
        if i == 7
            subplot(211)
            plot(out.time,out.c_setpoint,':','color',red);
            xline(tstep_d,':','color',green,'LineWidth',2);
            if (k>1)
                ldg1=legend('basecase','v_{1}','v_{2}','v_{3}','v_{4}','v_{5PI}','v_{5P}','c_s');
                ldg1.Box = 'off';
                ldg1.Orientation = 'horizontal';
                ldg1.NumColumns = 4;
            else
                ldg1=legend('basecase','v_{1}','v_{2}','v_{3}','v_{4}','c_s');
                ldg1.Box = 'off';
                ldg1.Orientation = 'horizontal';
                ldg1.NumColumns = 4;
            end
            ldg1.Location = 'best';
            fullname = [filepath,name,'.eps'];
            saveas(fig2,fullname,'epsc');
            close(fig2)
        end
        
%% plotting
            fig1=figure('visible','off');
            fig1.Position = [100,100,450,450];
            subplot(2,1,1)
            p1=plot(out.time,out.cm,'color',blue);
            hold on
            p2=plot(out.time,out.c_setpoint,':','color',red);
            ylim([0.6,1.5])
            ylabel('y=c_{m}, [kmol/m^3]')
            xline(tstep_d,':','color',green,'LineWidth',2);
            ldg2=legend([p1,p2],'y = c_{m}','c_s','Box','off','Orientation','horizontal');
            ldg2.Location = 'best';

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

            fullname = [filepath,casenr,name,'.eps'];
            saveas(fig1,fullname,'epsc');
            close(fig1)
        end
        l=l+1;
    end
end
close(progressbar)




%% have all the cases written down here as help functions
function [Kc,tauI,KI,par]=tuning_case(par,i)
switch i
    case 0
        [Kc,tauI,KI,par]=case0(par);
    case 1
        [Kc,tauI,KI,par]=case1(par);
    case 2
        [Kc,tauI,KI,par]=case2(par);
    case 3
        [Kc,tauI,KI,par]=case3(par);
    case 4
        [Kc,tauI,KI,par]=case4(par);
    case 5
        [Kc,tauI,KI,par]=case5(par);
    case 6 
        [Kc,tauI,KI,par]=case6(par,false);
    case 7
        [Kc,tauI,KI,par]=case6(par,true);
    otherwise
        Kc = 1;
        KI = 1;
end
end


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