% Task  : Transformed MVs. Tank. Study static, linear, int transf.
% Date  : August 2020
% Author: CZ
% TODO  :
% Var   : x=[T1;T2]; y = T2; w=T1; u=Q; ...
%         d=[T0; TD; qD; q]; v=[v0 or vL or vL0 or vF]

clc;clear;close all; warning off
%% Plotting options
set(0,'DefaultTextFontName','Times',...
     'DefaultTextFontSize',14,...
     'DefaultAxesFontName','Times',...
     'DefaultAxesFontSize',14,...
     'DefaultLineLineWidth',4,...
     'DefaultLineMarkerSize',7.75,...
     'DefaultStairLineWidth',2)

%% Parameters and nominal conditions
V1    = 10;                                            % [m3] tank 1 volume
V2    = 10;                                            % [m3] tank 2 volume
rho   = 1000;                                             % [kg/m3] density
cp    = 4.18;                                   % [kJ/kg/deg] specific heat
qnom  = 1e-2;                                    % [m3/s] nominal flow rate
qDnom = 2*qnom;                  % [m3/s] nominal flow rate added to tank 2
T0nom = 10;                                      % [degC] inlet temperature
TDnom = 20;                  % [degC] % temperature of flow added to tank 2
T1nom = 50;                             % [degC] tank 1 nominal temperature
Qnom  = qnom*(T1nom-T0nom)*rho*cp;                 % [kW] nominal heat rate
T2nom = 30;                             % [degC] tank 2 nominal temperature
Tnom = [T1nom;T2nom];                     % [degC] nominal temperature cond
%% Set up simulation
tsim     = 2e4;                                     % [s], simulation time

% Setpoint step change
t0sp = 1000; tfsp = t0sp+14e3;         % time when the step starts and ends
t    = [0 t0sp t0sp tfsp];                % time vector for setpoint change
dsp  = 10; % [degC] change in setpoint
sp   = [0 0 dsp dsp];                         % change in setpoint - vector
T2sp = timeseries(sp',t);                     % setpoint step as timeseries

% Disturbances step change
dnom  = [T0nom; TDnom; qDnom; qnom];
dstep = [15;-10 ;qDnom*(-0.5); qnom*(-0.5)];        % steps for disturbances
t0    = [5;7;9;11]*1e3;     % time when the step starts for each disturbance
tf    = ones(length(dnom),1)*tsim;            % time when disturbance stops
nd    = numel(dstep);
for i = 1:nd                        % construct timeseries for disturbances
     t         = [0 t0(i) t0(i) tf(i)];
     d         = [0 0 dstep(i) dstep(i)];
     dist {i}  = timeseries(d',t);
end

%% Select simulation case

lgd     = {'v0', 'vL','vL0','vI','Setpoint'};                 % legend text
lintype = {'-',':','--',':'};
ttl     = {'Tank 1-v0','Tank 1-vL','Tank 1-vL0','Tank 1-vI'};  % title text

A1     = -qnom/V1;      % tuning parameter for general linear transformation
A2     = -(qnom+qDnom)/V2;
tau01  = -1/A1*1;                          % [s] tuning param for linear IFAC
tau02  = -1/A2*1;


n        = 1;
umax     = 3*Qnom; ubias = Qnom;

vbias1   = [T1nom;-A1*T1nom; T1nom;0];  % biasses for tank 1 transformation
vbias2   = [T2nom;-A2*T2nom; T2nom;0];   % biases for tank 2 transformation

k1     = [1;-1/A1;1;1e6]; k2     = [1;-1/A2;1;0.6*1e6]; % process gain
tau1   = [ones(3,1)*tau01; 1e6];           % open loop time constant tank 1
tauC1  = ones(4,1)*250;    tauC2  = ones(4,1)*1000;  %[s] closed loop const 
theta1 = zeros(4,1);                                % [s] time delay tank 1

tau2   = [ones(3,1).*[1740 610 610 1e6]; 1330 530 530 1e6];
         %row:transf for tank 1, column:transf for tank 2
theta2 = [20 10 10 10; 20 10 10 20; 20 10 10 30; 10 5 5 5];
         %row:transf for tank 1, column:transf for tank 2
         

for i = 0:3 % for tank1=>0:static;1:linear; 2:lin k=1;3:int
     transf1 = i;
      
     for j =0:3
          transf2 = j;      % for tank2=>0:static;1:linear; 2:lin k=1;3:int
          vbias = [vbias1(i+1);vbias2(j+1)];              % controller bias
          
          p1 = [V1;V2;rho;cp;A1;A2;tau01; tau02;transf1];        % parameters for tank 1
          p2 = [V1;V2;rho;cp;A1;A2;tau01; tau02;transf2];        % parameters for tank 2
          
          k     = [k1(i+1);k2(j+1)]; 
          tau   = [tau1(i+1);tau2(i+1,j+1)];
          theta = [theta1(i+1);theta2(i+1,j+1)];
          tauC  = [tauC1(i+1);tauC2(j+1)];
          
          Kc = 1./k.*tau./(tauC+theta);
          tauI = min(tau,4.*(tauC+theta));
          
          sim('sim_uHeat')                                       % simulate
          
          figure(i+1)
          
          subplot(3,2,1)
          plot(time, T(:,1),lintype{j+1})  
          title(ttl{i+1})
          hold on
          xlabel('Time, (s)')
          ylabel('w=T1, (degC)')
          xlim([0 tsim])  
          
          subplot(3,2,2)
          plot(time, T(:,2),lintype{j+1})          
          hold on       
          xlabel('Time, (s)')
          ylabel('y=T2, (degC)')
          xlim([0 tsim])  
          
          subplot(3,2,3)
          plot(time, v1,lintype{j+1})           
          hold on
          xlabel('Time, (s)')
          ylabel('v1')          
          xlim([0 tsim])
          
          subplot(3,2,4)
          plot(time, v2,lintype{j+1})          
          hold on
          xlabel('Time, (s)')
          ylabel('v2')          
          xlim([0 tsim])
          
          subplot(3,2,5)
          plot(time, u,lintype{j+1})          
          hold on
          xlabel('Time, (s)')
          ylabel('u=Q, ')
          xlim([0 tsim] )      

% figure(n)
% subplot(1,2,1)
% plot(time,v2); xlabel('Time, (s)'); ylabel('v2'); legend(lgd{j+1})
% subplot(1,2,2)
% plot(time,T(:,2));xlabel('Time, (s)'); ylabel('v2');
% set(gcf,'Position',[100 100 400 200])  
%           n=n+1;
      
     end
     
     subplot(3,2,2); plot(time,Tsp,'k--');
     subplot(3,2,5);
     lg=legend(lgd); lg.Position = [0.7430 0.0812 0.1210 0.2317];
     set(gcf,'Position',[100*(j+1) 100*(j+1) 600 300])  
     
end


