% Task  : Transformed MVs. Tank. Study static, linear, int transf.
% Date  : August 2020
% Author: CZ
% TODO  :
% Var   : x=[T1;T2]; y = T2; w0 T1; u=(q or Q); ...
%         d=[T0; TD; qD; (Q or 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] nominal temperature
Qnom  = qnom*(T1nom-T0nom)*rho*cp;                 % [kW] nominal heat rate

T2nom = 30;
Tnom = [T1nom;T2nom];
%% Set up simulation
tsim     = 7000;                                     % [s], simulation time
tstep_u  = 1000;

% Disturbances step change
dstep = [T0nom*0; TDnom*0;qDnom*0];                % steps for first 3 disturbances
t0    = [10;20;30;40];                          % time when the step starts
tf    = t0+[1 1 1 1];                            % time when the step stops
nd    = numel(dstep);
for i = 1:nd % construct
     t         = [0 t0(i) t0(i) tf(i)];
     d         = [0 0 dstep(i) dstep(i)];
     dist {i}  = timeseries(d',t);
end

%% Open loop simulation
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;

t = [0 t0(4) t0(4) tf(4)];

lgd = {'v0', 'vL','vL0','vI','Setpoint'};          % legend text
ylb = {'Inflow, q (m3/s)';'Heat, Q (kW)'};          % yaxis label for input
lintype = {'-',':','--',':'};
n=1;
for i = 1
     for j = 2
          MV     = i;
          transf = j;
          
          switch MV
               case 1
                    unom    = qnom;
                    d4nom   = Qnom;
                    dist{4} = timeseries([0 0 Qnom*0 Qnom*0]',t);% set dist 
                    
               case 2
                    unom    = Qnom;
                    d4nom   = qnom;
                    dist{4} = timeseries([0 0 qnom*0 qnom*0]',t);% set dist 
          end
          dnom  = [T0nom; TDnom; qDnom; d4nom];%nominal values for disturbances
          
          p = [V1;V2;rho;cp;A1;A2;tau01; tau02;transf;MV];     % parameters
          
          % Input step change
          ustep = unom*1.1;                                      % [m3/sec] 
          d4_step = d4nom*1.1;
          sim('simOL')                                           % simulate
          
          figure(n)
          subplot(2,2,1)
          plot(time,u)
          hold on
          ylabel(ylb{i})
          xlabel('time,s')
          legend(lgd{j+1});
          xlim([0 tsim])
          
          subplot(2,2,2)
          plot(time,vc)
          hold on
          ylabel('v')
          xlabel('time,s')
          xlim([0 tsim])
          
          subplot(2,2,3)
          plot(time,T(:,1))
          hold on
          ylabel('w=T1')
          xlabel('time,s')
          xlim([0 tsim])
          
          subplot(2,2,4)
          plot(time,T(:,2))
          hold on
          ylabel('y=T2')
          xlabel('time,s')
          xlim([0 tsim])
          n=n+1;
          kuv(j+1,i) = (vc(1)-vc(end))/(u(1)-u(end)) 
     end
   
end


