% 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; 
%% 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     = 1e4;                                     % [s], simulation time

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

% Disturbances step change
dnom  = [T0nom; TDnom; qDnom; Qnom];
dstep = [15; -10;qDnom*(-0.5); Qnom*(-0.5)];            % disturbances step
t0d   = [5;7;9;11]*1e3;                                % disturbances start
tfsp  = ones(length(dnom),1)*30e4;                      % disturbances stop
nd    = numel(dstep);
for i = 1:nd                        % construct timeseries for disturbances
     td        = [0 t0d(i) t0d(i) tfsp(i)];
     d         = [0 0 dstep(i) dstep(i)];
     dist {i}  = timeseries(d',td);
end
                 

%% Select simulation case

lgd = {'v0', 'vL','vL0','vI','Setpoint'};          % legend text
ylb = {'Inflow, q (m3/s)';'Heat, Q (kW)'};          % yaxis label for input
lintype = {'-',':','--',':'};


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;

lgd = {'v0', 'vL','vL0','vI','Setpoint'};          % legend text
n   = 1;

for j =0
     MV     = 1;                                % MV=1=>u=q; MV=2=>u=Q
     transf = j; % 0=> static; =1=>linear, vL; =2=>linear,vL0;3=> integr
        
     d2nom   = Qnom;
     umax    = 3*qnom; umin=1e-4; ubias = qnom;
     thetauv = 630;    tauCS = 2*thetauv;                  % slave controller
     switch transf                            % select PI controller tuning
          case 0                                    % static transformation
               k     = 1; tau = 744; theta = 655;
               vbias = T2nom;
               kuv   = -645; tuv=816;
          case 1                                    % linear transformation
               k     = 333; tau = 1202; theta = 610;
               vbias = -A2*T2nom;
               kuv   = -1.83; tuv=815;
          case 2                                                 % IFAC vL0
               k     = 1; tau = 760; theta = 600;
               vbias = T2nom;
               kuv   = -645; tuv = 892;
               
          case 3                               % integrating transformation
               k = 1e10; tau = 1e10;
               vbias = 0;
               kuv   = -645; tuv=815;
               
          otherwise
     end
     tauC  = 5*tauCS;
     Kc    = 1/k*tau/(tauC+theta);      % master PI-ctrl
     tauI  = min(tau,4*(tauC+theta));
     
     KcS   = 1/kuv*tuv/(tauCS+thetauv); % slave PI-ctrl
     tauIS = min(tuv,4*(tauCS+thetauv));
     KIS   = 1/kuv*1/(tauCS+thetauv); % slave I-ctrl
     
     
     p = [V1;V2;rho;cp;A1;A2;tau01; tau02;transf];        % parameters
     
     sim('sim_uFlow')                                       % simulate
     
     figure(i)
     subplot(3,1,1)
     plot(time, T(:,2))
     hold on
     xlabel('Time, (s)')
     ylabel('Temperature, (degC)')
     legend(lgd{j+1})
     
     subplot(3,1,2)
     plot(time, u)
     hold on
     xlabel('Time, (s)')
     
     legend(lgd{j+1})
     
     subplot(3,1,3)
     plot(time, vsp)
     hold on
     xlabel('Time, (s)')
     ylabel('v')
     legend(lgd{j+1})
     
     n=n+1;
end


