% Task  : Transformed MVs. extraction process casestudy
% compare
% Date  : 14. April 2020
% Author: Simen Bjorvand
% TODO  : 
% Var   : MVs=[]; CVs=[]; States=[]  DV=[].

clc;clear;%close all
% profile on
%% Plotting options
set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',16,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',16,...
'DefaultLineLineWidth',2,...
'DefaultLineMarkerSize',7.75,...
'DefaultStairLineWidth',2)
%% Colors
blue    = [0 0.4470 0.7410];
red     = [0.83 0 0.33];
green   = [0 0.5 0];
indigo  = [0.2930 0 0.5078];
lilac   = [0.7383    0.5898    0.7930];
hotpink = [0.8945    0.4102    0.5898];
cerulean = [0.1094    0.6641    0.8359];

%% Parameters
%initial values for states
par.Tinit = 350;                             % [K] temperature of cold side

%initial values for input and disturbance variables
par.zinit = 0.1;                           % [-] steam inlet valve position

par.T0init = 300;                         % [K] cold side inlet temperature
par.winit = 10;                                     % [kg/s] cold side flow
par.Td0init = 390;                       % [K] steam uppstreams temperature


%other parameters
par.m = 1000;                                   % [kg] hold up og cold side
par.cp = 4;                           % [kJ/kgK] heat capacity of cold side
par.U = 1;                             % [kW/m2K] heat transfer coefficient 
par.A = 100;                                  % [m2] Area of heat transfer
par.lambda = 2200; 
               % [kJ/kg] The enthalpy of vaporization of medium in hot side
par.R = 8.314;                              % [kJ/kmolK] ideal gas constant
par.pdref = 1.0;                % [Bar] saturated steam refference pressure
par.Tdref = 370;               % [K] saturated steam refference temperature
par.Mm = 18;                                % [Kg/kmol] Molar mass of steam

Tdinit = 370;                                 % [K] Temperature of hot side
pdinit = 1.0;                                  % [Bar] Pressure of hot side

%parameters derived from the other parameters
Rs = par.R/par.Mm;
par.pd0init = par.pdref*exp(-par.lambda/Rs*(1/par.Td0init-1/par.Tdref));                           
                                          % [Bar] steam uppstreams pressure
                                          
Qinit = par.U*par.A*(Tdinit - par.Tinit);
wdinit = Qinit/par.lambda;
par.k = wdinit/(par.zinit*(par.pd0init-pdinit));
                                             % [kg/Bar*s] valve coefficient
                                             
% find different taus, non are "correct"
par.tau1 = par.m/par.winit;
par.tau2 = par.cp*par.m/(par.cp*par.winit+par.U*par.A);

%calculate As = df/dx|* and B = df/dz|*    
k1 = par.k*par.zinit*par.lambda^2*pdinit/(Rs*Tdinit^2);
par.k2 = k1/(par.U*par.A);
par.k3 = par.lambda*par.k/(par.U*par.A)*(par.pd0init-pdinit);

par.As = -(1/par.m)*(par.winit + (1/par.cp)*(k1/(1+par.k2))); 

par.B = (par.lambda*par.k/(par.m*par.cp)) * (par.pd0init-pdinit)/(1+par.k2);


v0 = 0;                              %update v0 in some help function later

%% Set up simulation
par.p   = [par.m       %1
           par.cp      %2                                
           par.U       %3
           par.A       %4
           par.lambda  %5
           par.k       %6
           par.R       %7
           par.pdref   %8 
           par.Tdref   %9
           par.Mm      %10    
           par.As      %11
           0           %12 indicates the structure which is used        
           0];         %13 indicates cascade and which variable is measured
                                                           % Parameter list
       
x0      = par.Tinit;                                        % initial state
       
u0      = par.zinit;                                    % physical MVs bias
  
d0      = [par.T0init
           par.winit
           par.pd0init];                                         % DVs bias

y0      = x0;

dv      = 0.1;
%% variables which will be varied to make up the 30 cases
tsim_0    = 500;                                     % [s], simulation time

tstep_y = tsim_0*0.02;                 % [s], step time for setpoint change

tstep_d = tsim_0*0.02;                    % [s], step time for disturbances



ys = y0;

 
d1_list = par.T0init *  [0.9, 1.1];                        % change in DV 1
d2_list = par.winit *   [0.9, 1.1];                        % change in DV 2
d3_list = par.pd0init * [0.9, 1.1];                        % change in DV 3

d1 = 1.0*par.T0init;                                       % change in DV 1
d2 = 1.0*par.winit;                                        % change in DV 2
d3 = 1.0*par.pd0init;                                      % change in DV 3



ds = [d1;d2;d3];


tauC = 50;
theta_list = [0,50];
vlist = [2,4,5,6];
name_case_list ={'v2','v4','v5','v6'};
colorList = {green,cerulean,hotpink,lilac,indigo};

%input gain error
Ku_list = [1.0 2.0];
Ku = 1.0;


IAE_list = zeros(4*1,6*4);
% out = cell(4,1);


filepath = '/Users/simenbjorvand/Documents/OneDrive - NTNU/master/heatExchanger/figures/disturbance/cascadePd'; 
h = waitbar(0,'Running simulation');
%% Simulation. Get PID tuning, and run simulation
count = 1;
for l = 1:2
    for f = 1:2
        for i = 1:1
            for k = 1:6
                for j = 1:4
                    if (l == 3 && f == 3)
                    else
                        Ku = Ku_list(l);
                        theta = theta_list(f);
                        [Kc,KI,Kci,KIi,par] = cascade_tunings(par,tauC,theta,vlist(j),2);
                        tsim = tsim_0+theta*10;
                        if (k == 1 || k == 2)
                            ds(1) = d1_list(k);
                            prename = ['T0disturbance',num2str(d1_list(k))];
                        elseif (k == 3 || k == 4)
                            ds(2) = d2_list(k-2);
                            prename = ['wdisturbance',num2str(d2_list(k-2))];
                        elseif (k == 5 || k == 6)
                            ds(3) = d3_list(k-4);
                            prename = ['pddisturbance',num2str(d3_list(k-4))];
                        end
                        
                        out = sim('HEXcascade');
                        IAE = out.IAE(end);
                        IAE_list(j,k+6*(f-1)+12*(l-1)) = IAE;
                        casenr = name_case_list{j};
                        name = ['theta',num2str(theta),'Ku',num2str(Ku)];

                        %plot some shite down here
                        if j == 1
                            fig1=figure('visible','off');
                        end

                        set(0,'Currentfigure',fig1);

                        subplot(211)
                        plot(out.time,out.T,'color',colorList{j})
                        hold on
                        ylabel('y = T [K]')
                        subplot(212)
                        plot(out.time,out.u,'color',colorList{j})
                        ylabel('u = z [-]')
                        xlabel('time [s]')
                        hold on 
                        if j == 4
                            subplot(211)
                            plot(out.time,out.ys,'--','color',red)
                            legend('v_{2}','v_{4}','v_{5}','v_{6}')

                            fullname = [filepath,'/compare/',prename,name,'.eps'];
                            saveas(fig1,fullname,'epsc');
                            close(fig1)
                        end

                        fig2 = figure('visible','off');
                        subplot(311)
                        plot(out.time,out.T)
                        hold on
                        plot(out.time,out.ys,'--','color',red)
                        ylabel('y = T [K]')
                        legend('y = T','y_{s}')
                        subplot(312)
%                         plot(out.time,out.wd)
                        plot(out.time,out.Td)
                        hold on
                        plot(out.time,out.yis,'--','color',red)
%                         ylabel('yi = w_{d} [kg/s]')
                        ylabel('yi = T_{d} [K]')
                        legend('y_{i} = T_{d}','y_{is}')
                        subplot(313)
                        plot(out.time,out.u)
                        ylabel('u = z [-]')
                        xlabel('time [s]')

                        fullname = [filepath,'/individual/',casenr,prename,name,'.eps'];
                        saveas(fig2,fullname,'epsc');
                        close(fig2)
                    end
                    
                    waitbar(count/(4*1*6*4));
                    count = count+1;
                end
            end
        end
    end
end
close(h)