% Task: SRC for tanks in series with moving bottlnecks
% Date: April 2021

clc; clear; close all

%% Initial conditions
q0   = 1;                                      % nominal flowrates [m3/min]
zOP0 = [1;1;1;1];             % initial operator setpoints for the 4 valves
hmax = [2.3; 2.8; 3.2];              % [m] max L=[2.3; 2.8; 3.2]. from h=2d
hH   = 0.5;                                     % level high setpoint = 90%
hL   = 0.5;                                     % level low setpoint = 10 %
d0   = 0;                                                    % disturbances
%% Select TPM
TPM = 3;                                    % select where to place the TPM
switch TPM
     case 0 % TPM = z0
          z0 = [1;   0.8; 0.7; 0.6];                   % operator set flows
          x0 = [hL; hL; hL].*hmax;                     % initial tank level
          load SRCparamTPMz0
     case 1 % TPM = z1
          z0 = [0.8; 1;   0.7; 0.6];                  % operator set valves
          x0 = [hH; hL; hL].*hmax;                     % initial tank level
          load SRCparamTPMz1
     case 2 % TPM = z2
          z0 = [0.8; 0.7; 1; 0.6];                    % operator set valves
          x0 = [hH; hH; hL].*hmax;                     % initial tank level
          load SRCparamTPMz2
     case 3 % TPM = z3
          z0 = [0.6; 0.7; 0.8; 1];                    % operator set valves
          x0 = [hH; hH; hH].*hmax;                     % initial tank level
          load SRCparamTPMz3
     otherwise
          disp('Incorrect TPM')
end

%% Parameters
A     = [1; 1.5; 2];                        % [m3] tanks cross section area
rho   = 1000;                                       % [kg/m3] water density
N     = 3;                                                % number of tanks
DP    = [0.5; 0.6; 0.7; 1]*1e5;                                       %[Pa]
Cv    = q0./(z0.*sqrt(DP/rho));                  % linear valve coefficient
param = [A; Cv; DP; rho; N];% vector of parameters passed to the ode system

%% Set-up simulation
nx     = length(x0); nu = length(z0); nd = length(d0); np = length(param);
td     = 0;           % [min] disturbances step time. given for consistency
df     = 0;                                               % no disturbances
th     = ones(6,1);              % [min] time for setpoint change. not used
tsim   = 120;                                       % [min] simulation time
tdelay = 1;

%% Change in TPM
zOP      = zOP0-[0.7;0;0;0];                  % operator sets valve positions
%duration = 'short';                                       % short bottlneck
duration = 'long';                                   % prolonged bottlneck
switch duration
     case 'short'
          tBend = 19;
     case 'long'
          tBend = 400000;
     otherwise
          disp('error duration TPM')
end
tz0OP = [10; tBend];  % [min] time for operator decrease and increase in z0
tz1OP = [10; tBend];  % [min] time for operator decrease and increase in z1
tz2OP = [10; tBend];  % [min] time for operator decrease and increase in z2
tz3OP = [10; tBend];  % [min] time for operator decrease and increase in z3

%% Simulate SRC
Kc     = Kc*100; 
tauI   = tauI*1;
pr     = 50;                                                % setpoint in %
x0     = pr/100*hmax;                                   % initial inventory
h0     = x0;                                               %  setpoint in m
hf     = h0;                                               % final setpoint
MV     = 2;                                    % number of MVs for each SRC
pSR1   = [u01, alpha1, MV];             % parameters to SR block for tank 1
pSR2   = [u02, alpha2, MV];            % parameters to SRC block for tank 2
pSR3   = [u03, alpha3, MV];            % parameters to SRC block for tank 3

sim('simSRC')

%% Save results for controllers with different setpoint and SRC
idx         = find(zOP~=1)-1; % find which valve was changed
filenameSRC = strcat('SRC_','TPM_z',num2str(TPM),'_h',num2str(pr),'_',duration,'_z',num2str(idx));   % create filename
save(filenameSRC,'tSRC','uSRC','ySRC','h0','zOP','z0');

%% Plot results
FSRC = uSRC'.*Cv.*sqrt(DP/rho); % flows
figure('Name','SRC')
subplot(3,1,1)
plot(tSRC,ySRC./hmax'*100)
legend('h1', 'h2', 'h3')
xlabel('Time [min]')
ylabel('Level [%]')

subplot(3,1,2)
plot(tSRC,uSRC'*100)
legend('z0','z1', 'z2', 'z3')
xlabel('Time [min]')
ylabel('Valve opening [%]')

subplot(3,1,3)
plot(tSRC, FSRC)
legend('F0','F1','F2','F3')
xlabel('Time [min]')
ylabel('Flowrate [m3/min]')
