% Task  : Transformed MVs. SISO mixing tank 30-cases, all cases, maybe
% compare
% Date  : 5 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    =  [300,60,30];                  % [s], simulation time, small
% tsim_list    =  [60,30,30];                 % [s], simulation time, small

tstep_c_list = tsim_list/10;% [s], step time for concentration setpoint changes

tstep_d = 10000;                          % [s], step time for disturbances
par.gain_error  =  2.0;                      % input error gain proportions
par.disturbance_error = eye(4);       % disturbance gain error, proportions
par.disturbance_error(2,2) = 1.0;

% cs1 = 1.2;              % setpoint change for concentration setpoint, small
cs1 = 1.7;                                                          % large
% limity = [1.05,1.30];
limity = [0.9,2.0];
% limitu = [0,8]; 
limitu = [-2,120];


theta_list  = [3*par.taur;0.3*par.taur;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      = 1;                      % step in inputs for open loop simulation
ds=d0;

%% Simulation. Get PID tuning, and run simulations
IAE = zeros(3,8);
l=1;
count=1;
progressbar = waitbar(0,'Progress');

filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/Mixingtank30-cases/figures/gainerror/setpointchange';  %tauc>theta,v5=p_control
for k =1%1:3
    for i = 6%[0:3,5:7]
        par.theta = theta_list(k);
        tsim = tsim_list(k);
        tstep_c=tstep_c_list(k);
        [Kc,tauI,KI,par]=tuning_case(par,i);
        name = ['/',num2str(cs1-par.c0),'_ys_change_',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/24);
        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];
        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(limity)
            ylabel('y=c_{m}, [kmol/m^3]')
            subplot(2,1,2)
            plot(out.time,out.u,'color',colorList{i+1})
            hold on
            ylim(limitu)
            ylabel('u=q_1 [m^3/s]')
            xlabel('Time,[s]')
        end
        
        if i == 7
            subplot(211)
            plot(out.time,out.c_setpoint,':','color',red)
            if (k>0)
                 ldg=legend('basecase','v_{1}','v_{2}','v_{3}','v_{4}','v_{5PI}','v_{5P}','c_s');
                 ldg.NumColumns = 4;
            else
                ldg=legend('basecase','v_{1}','v_{2}','v_{3}','v_{4}','c_s');
                ldg.NumColumns = 3;
            end
            ldg.Box = 'off';
            ldg.Location = 'best';
            ldg.Orientation = 'horizontal';
            fullname = [filepath,name,'.eps'];
            saveas(fig2,fullname,'epsc');
            close(fig2)
        end
        


%% plotting 2
        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(limity)
        ylabel('y=c_{m}, [kmol/m^3]')
        ldg2 = legend([p1,p2],'y = c_{m}','c_s','Location','southeast');
        ldg2.Box = 'off';
        ldg2.Location = 'best';

        subplot(2,1,2)
        plot(out.time,out.u)
%         ylim(limitu)
        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
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