% Task  : Transformed MVs. Tank. Study static, linear, int transf.
% Date  : July 2020
% Author: CZ
% TODO  : compare with vL0.
% Var   : x=y=T; u=(q or Q); d=[T0; (Q or q) ]; v=[v0 or vL 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
V    = 10;                                               % [m3] tank volume
rho  = 1000;                                              % [kg/m3] density
cp   = 4.18;                                    % [kJ/kg/deg] specific heat
qnom = 1e-2;                                     % [m3/s] nominal flow rate
Qnom = 1672;                                      % [kW] nominal heat rate
T0nom= 10;                                       % [degC] inlet temperature
%% Set up simulation
tsim    = 1000;                                      % [s], simulation time
tstep_d = 10;                                  % [s], disturbance step time
tstep_T = 500;                                  % [s], setpoint change time

% Stepoint
Tnom  = 50;                                    % [degC] nominal temperature
Tstep = Tnom-10;                      % [degC] setpoint change for temperature

% Disturbances
dnom = zeros(2,1);
dnom(1) = T0nom;

%% Select simulation case
n = 1;                                          % to count number of figure
ttl = {'Static', 'Linear','Integrating'};                % title of figures
lgd = {'Static', 'Linear','Integrating','Setpoint'};          % legend text
ylb = {'Inflow, q (m3/s)';'Heat, Q (kW)'};          % yaxis label for input
lintype = {'-',':','--'};

tauC  = 100;                        % [s] desired closed loop time constant
theta = 0;                                                 % [s] time delay
tau   = V/qnom;      % = -1/A; A =(df/dT)*. tuning param for general linear
k     = V/qnom;
A    = -qnom/V;
tau0 = -1/A;
tauq = tau;1250;
tauQ = tau;1000;

for i =1:2
    for j =0
        MV     = i;                                  % MV=1=>u=q; MV=2=>u=Q
        transf = j;               %0=> static; =1=>linear; =2=> integrating
        
        switch MV
            case 1 % u=q
                dnom(2) = Qnom;
                Kc=tauq/(tauC+theta); tauI=min(tauq,4*(tauC+theta));
                KIS = -1/4000*1/1;               % slave time constant = 1;
                ubias = qnom; vbias = Tnom;
                umax = 3*qnom;
                
            case 2 % u=Q
                dnom(2) = qnom;
                Kc=tauQ/(tauC+theta);tauI=min(tauQ,4*(tauC+theta));
                KIS = 1/0.024*1/1;               % slave time constant = 1;
                ubias = Qnom; vbias = Tnom;
                umax = inf;
            otherwise
        end
        d_step = dnom + dnom.*[0;0];                    % disturbance step
        
        p = [V;rho;cp;A;tau0;MV;transf];
        
        sim('simCascade')                                        % simulate
        figure(i)
        subplot(3,1,1)
        plot(time, T)
        hold on
        xlabel('Time, (s)')
        ylabel('Temperature, (degC)')
        
        subplot(3,1,2)
        plot(time, u,lintype{j+1})
        hold on
        plot(time, uc,lintype{j+1})
        xlabel('Time, (s)')
        ylabel(ylb{i})
        
        subplot(3,1,3)
        plot(time,vc)
        hold on
        plot(time,v)
        xlabel('Time,(s)')
        ylabel('Transformed input')
        
        n = n+1; % count figure number
        
    end
end

figure(1); subplot(311); plot(time,Tsp,'k:'); legend(lgd);
           subplot(313); legend('v calculated','v setpoint');
% figure(2); subplot(311); plot(time,Tsp,'k:'); legend(lgd);
%            subplot(313); legend('v calculated','v setpoint');
