
% This script calculates the maximum deviation from the selected setpoints
% for the outlet concentration from each tank for a pH neutralization
% process with three tanks in series.

% The disturbances are the following:
%       Tank 1:     cin1,        qin1, cr1, qr1
%       Tank 2:     cin2 = c1max,      cr2, qr2
%       Tank 3:     cin3 = c2max,      cr3, qr3
% Changes in the flowrate q are neglected after the first tank.

% It passes the parameters to the function nomcalc, which
% calculates the nominal values for the inlet and outlet streams to each of
% the three tanks, e.g. flowrate of base and acid, concentration of acid,
% flowrate out of the tanks and the residence time tau in each tank.
% The function then uses the nominal values returned by nomcalc to
% calculate the steady-state gains for different disturbances and the
% control input gains for each tank.

clear all
clc

% Setpoints for the outlet concentration from each tank (excess
% H+ [mol/l] (Selected arbitrarily):

sp = [0.1 0.01 0];   % Tank 1, Tank 2, Tank 3

% Time delay in each tank [s]:
theta = 10;

% Volume of each tank [l]
V = 2*40.7/3*1000;

% Controller settings:
Kc = -133;
tau_i = 80;

% Disturbances :
d_cin1 = 0.5;       % First tank (50%)
d_qin1 = 0.5;       % First tank (50%)
d_cr  = 0.1;        % All tanks  (10%)
d_qr  = 0.1;        % All tanks  (10%)

% Nominal values::
qin1_nom = 5;       % Nominal Flow rate of acid to Tank 1, [l/s]
pH_acid = -1;       % pH of feed stream (pH=-1 --> c = 10 mol/l excess H+)
pH_base = 15;       % pH of all reagent streams (pH=15 --> c = -10 mol/l excess H+)

% 
cin1_nom = 10^(-pH_acid)-10^(-14+pH_acid); % Concentration (excess H+) in acid stream to Tank 1. (can also call conc_calc)
cr_nom = 10^(-pH_base)-10^(-14+pH_base);   % Concentration (excess H+) of base added to each tank.

% Collecting the parameters in an array:
par_nom = [V cin1_nom qin1_nom cr_nom];

% Calculating nominal values
N = nomcalc(sp,par_nom);  % Returns nom = [qr_nom ; q_nom ; cin_nom ; qin_nom ; tau];

qr_nom = N(1,:);
q_nom = N(2,:);
cin_nom = N(3,:);
qin_nom = N(4,:);
tau = N(5,:);

% Calculating the steady state gains for different disturbances and the gain from 
% the control input:
K = kcalc(N,cr_nom,sp);

% Extracting the arrays from the matrix K:
ktank1 = K(1,:);
ktank2 = K(2,:);
ktank3 = K(3,:);

% Process gains
Kp = [ktank1(5) ; ktank2(5) ; ktank3(5)];

% Defining s as the Laplace variable s
s = tf('s');

% Controller TF:
Gc = Kc*(1+1/(tau_i*s));

% Process model TF and sensitivity for each tank:
for i=1:3
    
    % Process model TF
    G(i) = Kp(i)/(tau(i)*s+1);             
    
    % Sensitivity;    y(i) = sens(i)*gd(i)*d(i)
    sens(i) = 1/(1 + G(i)*Gc*exp(-theta*s)/(1+G(i)*Gc*(1-exp(-theta*s))));
    
end

% SGd transfer functions for each disturbance in Tank 1:
SGd_tank1_c  = sgd(ktank1(1),sens(1),tau(1));
SGd_tank1_q  = sgd(ktank1(2),sens(1),tau(1));
SGd_tank1_cr = sgd(ktank1(3),sens(1),tau(1));
SGd_tank1_qr = sgd(ktank1(4),sens(1),tau(1));

% SGd transfer functions for each disturbance in Tank 2:
SGd_tank2_c  = sgd(ktank2(1),sens(2),tau(2));
SGd_tank2_q  = sgd(ktank2(2),sens(2),tau(2));
SGd_tank2_cr = sgd(ktank2(3),sens(2),tau(2));
SGd_tank2_qr = sgd(ktank2(4),sens(2),tau(2));

% SGd transfer functions for each disturbance in Tank 3:
SGd_tank3_c  = sgd(ktank3(1),sens(3),tau(3));
SGd_tank3_q  = sgd(ktank3(2),sens(3),tau(3));
SGd_tank3_cr = sgd(ktank3(3),sens(3),tau(3));
SGd_tank3_qr = sgd(ktank3(4),sens(3),tau(3));
 
% Multiplying the norm(max values) of sens(i)*Gd(i) with the maximum expected
% disturbance to obtain the maximum deviation from the setpoint for the
% outlet concentration from each tank:

%__________________________________________________________________________
%/////////////////////////////////TANK 1: /////////////////////////////////

% Disturbances, Tank 1:
dc1  = d_cin1*cin1_nom;
dq1  = d_qin1*qin1_nom; 
dcr1 = d_cr*cr_nom;
dqr1 = d_qr*qr_nom(1);


% _________________________Single disturbances:____________________________

% When d = c_in:
 c1(1) = norm(SGd_tank1_c,inf)*dc1;
 
% When d = q_in:
 c1(2) = norm(SGd_tank1_q,inf)*dq1;
 
% When d = cr1:
 c1(3) = norm(SGd_tank1_cr,inf)*dcr1;
 
% When d = +qr1:
 c1(4) = norm(SGd_tank1_qr,inf)*dqr1;
 
%__________________________Multiple disturbances:__________________________

% Four disturbances at a time:

c1m(1) = + c1(1) + c1(2) + c1(3) + c1(4);
c1m(2) = + c1(1) + c1(2) + c1(3) - c1(4);
c1m(3) = + c1(1) + c1(2) - c1(3) + c1(4);
c1m(4) = + c1(1) + c1(2) - c1(3) - c1(4); 
c1m(5) = + c1(1) - c1(2) + c1(3) + c1(4);
c1m(6) = + c1(1) - c1(2) + c1(3) - c1(4);
c1m(7) = + c1(1) - c1(2) - c1(3) + c1(4);
c1m(8) = + c1(1) - c1(2) - c1(3) - c1(4);

%________________________Max. output from Tank 1:__________________________

c1max = max(abs(c1m));


%__________________________________________________________________________
%/////////////////////////////////TANK 2: /////////////////////////////////

dc2  = c1max;
dq2  = 0;
dcr2 = dcr1;
dqr2 = d_qr*qr_nom(2);

% _________________________Single disturbances:____________________________
% When d = c1max:
 c2(1) = norm(SGd_tank2_c,inf)*dc2;
 
% When d = cr2:
 c2(2) = norm(SGd_tank2_cr,inf)*dcr2;
  
% When d = qr2:
 c2(3) = norm(SGd_tank2_qr,inf)*dqr2;
 
%_________________________Multiple disturbances:___________________________

c2m(1) = + c2(1) + c2(2) + c2(3);
c2m(2) = + c2(1) + c2(2) - c2(3);
c2m(3) = + c2(1) - c2(2) + c2(3);
c2m(4) = + c2(1) - c2(2) - c2(2);

%________________________Max. output from Tank 2:__________________________ 

c2max = max(abs(c2m));

%__________________________________________________________________________
%///////////////////////////////// TANK 3: ////////////////////////////////

dc3  = c2max;
dq3  = 0;
dcr3 = dcr1;
dqr3 = d_qr*qr_nom(3);

% _________________________Single disturbances:____________________________
% When d = c2max:
 c3(1) = norm(SGd_tank3_c,inf)*dc3;
 
% When d = cr2:
 c3(2) = norm(SGd_tank3_cr,inf)*dcr3;
  
% When d = qr2:
 c3(3) = norm(SGd_tank3_qr,inf)*dqr3;
 
%________________________Multiple disturbances:____________________________

c3m(1) = c3(1) + c3(2) + c3(3);
c3m(2) = c3(1) + c3(2) - c3(3);
c3m(3) = c3(1) - c3(2) + c3(3);
c3m(4) = c3(1) - c3(2) - c3(3);
   
%________________________Max. output from Tank 3:__________________________

c3max = max(abs(c3m));

%__________________________________________________________________________

% Finding the indices of the worst-case combinations:

index1 = find(c1m==(max((max(c1m)))));
index2 = find(c2m==(max((max(c2m)))));
index3 = find(c3m==(max((max(c3m)))));

% -------------------Displaying the results:-------------------------------

format shortEng

if index1 == 1
    comb1 = 'cin1, qin1, cr1, qr1';
elseif index1 == 2
    comb1 = 'cin1, qin1, cr1, -qr1';
elseif index1 == 3
    comb1 = 'cin1, qin1, -cr1, qr1';
elseif index1 == 4
    comb1 = 'cin1, qin1, -cr1, -qr1';
elseif index1 == 5
    comb1 = 'cin1, -qin1, cr1, qr1';
elseif index1 == 6
    comb1 = 'cin1, -qin1, cr1, -qr1';
elseif index1 == 7
    comb1 = 'cin1, -qin1, -cr1, qr1';
else
    comb1 = 'cin1, -qin1, -cr1, -qr1';
end
            
            
if index2 == 1
    comb2 = 'cin2, cr2, qr2';
elseif index2 == 2
    comb2 = 'cin2, cr2, -qr2';
elseif index2 == 3
    comb2 = 'cin2, -cr2, qr2';
else
    comb2 = 'cin2, -cr2, -qr2';
end
    

if index3 == 1
    comb3 = 'cin3, cr3, qr3';
elseif index2 == 2
    comb3 = 'cin3, cr3, -qr3';
elseif index2 == 3
    comb3 = 'cin3, -cr3, qr3';
else
    comb3 = 'cin3, -cr3, -qr3';
end

text1 = 'Tank 1: ';
text2 = 'Tank 2: ';
text3 = 'Tank 3: ';


disp('Maximum allowed dev. from sp, Tank 3: +-1e-06')

disp(' ')

disp('Selected setpoints:')
disp([text1, num2str(sp(1))])
disp([text2, num2str(sp(2))])
disp([text3, num2str(sp(3))])

disp(' ')

disp('Calculated max. deviation from setpoint:')
disp([text1, num2str(c1max)])
disp([text2, num2str(c2max)])
disp([text3, num2str(c3max)])

disp(' ')

disp('Worst-case disturbance combinations:')
disp([text1, comb1])
disp([text2, comb2])
disp([text3, comb3])
