% Task  : Transformed MVs. Tank. Study static, linear, int transf.
% Date  : July 2020
% Author: CZ
% TODO  : 
% Var   : x=y=T; u=(q or Q); d=[T0; (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)

mygreen = [0 0.5 0]; hotpink = [0.8945    0.4102    0.5898];
clr     = [0 0 1; 1 0 0];
%% 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     = 1500;                                     % [s], simulation time
tstep_T1 = 100;                                 % [s], setpoint change time
tstep_T2 = 3000;                                % [s], setpoint change time
tstep_d1 = 150;                               % [s], disturbance step time
tstep_d2 = 20000;                               % [s], disturbance step time

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

% Disturbances
d1nom   = T0nom;
d1_step = 60;                     % [degC] disturbance in inlet temperature
d1stepcb = [d1_step;T0nom;];                     % [degC] disturbance in inlet temperature

%% Select simulation case
                                         % to count number of figure
ttl = {'Static', 'Linear','Integrating'};                % title of figures
lgd = {'v0', 'vA','T_0','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
tauT  = 400;                             % [s] tracking time constant for v

tau   = V/qnom;      % = -1/A; A =(df/dT)*. tuning param for general linear
k     = V/qnom;
A     = -qnom/V;       % tuning parameter for general linear transformation
tau0  = -1/A;                            % [s] tuning param for linear IFAC
vbias_step = -A*(Tnom+Tstep);
tauq = tau;
tauQ = tau;
nr=1;
for n=1 % 1- no measurment error; 2-use nominal disturbance values
     
     d1_step_cb=d1stepcb(n);
     
for i =1
    for j =0:1
        MV     = i;                                  % MV=1=>u=q; MV=2=>u=Q
        transf = j; % 0=> static; =1=>linear, vL; =2=>linear,vL0;3=> integr
        
        switch MV
            case 1 % u=q
                d2nom = Qnom;
                umax = 30*qnom;
                switch transf                 % select PI controller tuning
                    case 0                          % static transformation
                        Kc=tauq/(tauC+theta); tauI=min(tauq,4*(tauC+theta));
                        vbias= Tnom ; vmin = T0nom+1e-4; % to avoid :0
                        KIS = -1/4000*1/1;                        
                    case 1                          % linear transformation
                        Kc=1/(tauC+theta);  tauI=min(tau,4*(tauC+theta));
                        vbias = -A*Tnom; vmin=-inf;
                    case 2                                       % IFAC vL0
                        Kc=tau/(tauC+theta); tauI=min(tau,4*(tauC+theta));
                        vbias= Tnom ; vmin=-inf;
                    case 3                     % integrating transformation
                        Kc=1/(tauC+theta); tauI=4*(tauC+theta);
                        vbias = 0; vmin=-inf;
                    otherwise
                end
                
            case 2 % u=Q
                d2nom = qnom;
                umax = 3*Qnom;
                switch transf                 % select PI controller tuning
                    case 0                     % static transformation (v0)
                        Kc=tauQ/(tauC+theta);tauI=min(tauQ,4*(tauC+theta));
                        vbias = Tnom ; vmin = T0nom+1e-4; % to avoid :0
                        KIS = 1/0.024*1/1;       % slave time constant = 1;
                        ubias = Qnom; umax = 3*Qnom;
                        
                    case 1                       % linear transformation vL
                        Kc=1/(tauC+theta); tauI=min(tau,4*(tauC+theta));
                        vbias = -A*Tnom; vmin=-inf;
                    case 2                                       % IFAC vL0
                        Kc=tau/(tauC+theta); tauI=min(tau,4*(tauC+theta));
                        vbias= Tnom ;
                    case 3                     % integrating transformation
                        Kc=1/(tauC+theta); tauI=4*(tauC+theta);
                        vbias =0; vmin=-inf;
                    otherwise
                end
                
            otherwise
        end
        
        d2_step = d2nom*0.5;                             % disturbance step
               
        
        p = [V;rho;cp;A;tau0;MV;transf];
        
       
        sim('simTransformedMVs')                                 % simulate
        %sim('simTransformedMVsAW')   
        %sim('simTransformedMVsNum') 
        %sim('simCascade')
        %end
        figure(1)
        %subplot(3,1,1)
        plot(time, T,lintype{j+1},'Color',clr(j+1,:))
        hold on
        xlabel('Time, (s)')
        ylabel('Temperature, (^oC)')
        %ylim([48 53])
        set(gcf,'Position',[100 100 600 300])
        
    
        %subplot(3,1,2)
        figure()
        plot(time, u,lintype{j+1},'Color',clr(j+1,:))
        hold on
        xlabel('Time, (s)')
        ylabel(ylb{i})
        %ylim([0 2000])
        set(gcf,'Position',[100 100 600 300])
        
        
        %subplot(3,1,3)
        figure()
        plot(time,v,lintype{j+1},'Color',clr(j+1,:))
        hold on
        xlabel('Time,(s)')
        ylabel('Transformed input')
        %ylim([0.04 0.06])
        set(gcf,'Position',[100 100 600 300])
        
       
        
    end
end

end
figure(1); plot(time,d(:,1),'Color', mygreen);plot(time,Tsp,'k:'); legend(lgd);
%figure(2); subplot(311); plot(time,Tsp,'k:'); legend(lgd);
%figure(1);  plot(time,Tsp,':','Color','k');legend(lgd);

