%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%                                              %%%%%
%%%%%            PERSTORP DYNAMIC STUDY            %%%%%
%%%%%                                              %%%%%
%%%%%          Daniel Greiner Edvardsen            %%%%%
%%%%%                 IKP, NTNU                    %%%%%
%%%%%               June 14, 2011                  %%%%%
%%%%%                                              %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all
load 2010.mat;
k=4000;

% Loading relevant temperatures, mass streams, heat capacities, 
% and densities from 'perstorpdata.m'
[T0 Tp Th11in Th1in Th2in Th3in Th11 Th1 Th2 Th3 T11_ T11 T1 T2 T3 m0 mh11 mh1 mh2 mh3 m11 m1 m2 m3 cp0 cp11 cp1 cp2 cp3 rho0 rho11 rho1 rho2 rho3 w0 w11 w1 w2 w3 wh11 wh1 wh2 wh3 ] = perstorp_data(k);

Th11in  = double(Th11in);
Th1in	= double(Th1in);
Th2in 	= double(Th2in);
Th3in	= double(Th3in);

% ========================================================
%               MAIN HEAT EXCHANGERS
% ========================================================

% Heat transfer areas given from Perstorp [m^2]
A11 = 80;
A1  = 12.3;
A2  = 198;
A3  = 198;

% Volumes given from Perstorp [m^3]
Vh_tot = 0.775;     % Hot side volume [m3] for HE2 and HE3
Vc_tot = 0.760;     % Cold side volume [m3] for HE2 and HE3
Vh11_tot = 0.234;   % Hot side volume [m3] for HE11
Vc11_tot = 0.232;   % Cold side volume [m3] for HE11
Vh1_tot = 0.02226;  % Hot side volume [m3] for HE1
Vc1_tot = 0.02164;  % Cold side volume [m3] for HE1

N = 10;                 % Number of cells (Fixed)

% Estimation of heat transfer coefficients [W/m2*K]
h_h11 = 60.15;
h_h1  = 633;
h_h2  = 5615;
h_h3  = 5800;

h_c11 = h_h11;
h_c1  = h_h1;
h_c2 = h_h2;
h_c3 = h_h3;

% Hot side
Vh11 = Vh11_tot/N;  % Cell volume [m3] for HE11
Ah11 = A11/N;       % Cell area [m2] for HE11
Vh1 = Vh1_tot/N;    % Cell volume [m3] for HE1
Ah1 = A1/N;         % Cell area [m2] for HE1
Vh2 = Vh_tot/N;     % Cell volume [m3] for HE2
Ah2 = A2/N;          % Cell area [m2] for HE2
Vh3 = Vh_tot/N;     % Cell volume [m3] for HE3
Ah3 = A3/N;          % Cell area [m2] for HE3

% Cold side

Vc11 = Vc11_tot/N;
Ac11 = A11/N;
Vc1 = Vc1_tot/N;
Ac1 = A1/N;
Vc2 = Vc_tot/N;	
Ac2 = A2/N;
Vc3 = Vc_tot/N;	
Ac3 = A3/N;		

% Wall
mass_w = 568;           % Weight of heat exchanger [kg] 
rho_w  = 1240;          % Wall density [kg/m3] (TM10-BFG)
V_w    = mass_w/rho_w;  % Wall volume [m3]
cp_w   = 3000;          % Wall heat capacity [J/kg*K]   
V_i_w  = V_w/N;         % Volume of one wall segment [m3]
A_i_w  = 200/N;         % Area of one wall segment [m2]



% ========================================================
%                  INITIAL CONDITIONS [°C]
% ========================================================

% Initial conditions for the heat exchangers, all temperatures:

% First column: Initial conditions, hot side (10)
% Second column: Initial conditions, wall (10)
% Third column: Initial conditions, cold side (10)

HX11init = [124.829214721733,122.027724998537,119.222846071090,116.414573839179,113.602904197629,110.787833036302,107.969356240083,105.147469688884,102.322169257627,99.4934508162486;
            104.126288990191,101.299753148383,98.4697978018694,95.6364188137799,92.7996120422395,89.9593733403622,87.1156985562454,84.2685835329632,81.4180241085612,78.5640161160495;
            57.6345814158503,60.5138789594949,63.3896973770426,66.2620408724073,69.1309136444229,71.9963198868499,74.8582637883812,77.7167495326487,80.5717812982289,83.4233632586492];

HX1init  = [131.401156119446,130.791098624113,130.133243058575,129.423844438905,128.658864361280,127.833948012509,126.944399379317,125.985154515285,124.950752713251,123.835305419060;
            122.180559708185,120.848067494867,119.411174373690,117.861700522323,116.190825228477,114.389036676111,112.446077797376,110.350887882056,108.091539612105,105.655171162826;
            87.4750369065912,91.2323265109586,94.7166212488277,97.9477562154352,100.944125339713,103.722786095673,106.299556605740,108.689105688805,110.905036365621,112.959963296923];
        
HX2init  = [138.846162306895,126.020623397690,114.568023658816,104.341393674112,95.2094967049858,87.0551445499288,79.7736936867086,73.2717023984603,67.4657316507897,62.2812743317322;
            134.347502135017,122.003532601391,110.980951541971,101.138307902861,92.3492925820153,84.5011175243392,77.4930683302683,71.2352118036758,65.6472418532726,60.6574489371328;
            59.0336235425334,63.8287520557554,69.1987212088914,75.2124429738280,81.9470904987497,89.4890884590447,97.9352221316091,107.393879425125,117.986441805093,129.848841963138];
        
HX3init  = [128.956874151343,119.294028842280,110.285462669462,101.886873858987,94.0569603461932,86.7572166629025,79.9517445775879,73.6070765572558,67.6920111828653,62.1774597089048;
            126.654079683043,117.147158465433,108.283958636675,100.020893296551,92.3173268494088,85.1353751689136,78.4397192938219,72.1974317385598,66.3778145644549,60.9522484152974;
            59.7270371216900,65.0636179460445,70.7877869198639,76.9276940100560,83.5135336749246,90.5776933526244,98.1549127341161,106.282454603887,115.000288088586,124.351285214742];


        
% Controllers

tc10 = -0.634632447104129;
tc20 = 0.278170600005430;

%% Tuning parameters

% Controller 1
kc1 = 0.026;
tc1 = 62;
I1  = kc1/tc1;

% Controller 2
kc2 = 0.102;
tc2 = 174;
I2  = kc2/tc2;

%% Simulation

%Simulation time
simtime=3600;

sim('perstorp_sim',simtime);
systemHandle = get_param(gcs, 'Handle');
saveas(systemHandle , 'simulink.pdf', 'pdf')
%% PLOTS
close all;

%Convert from seconds to minutes
time=t/60;

figure
plot(time,T11(:,1),time,T1,time,T2,time,T3,time,Tend,'LineWidth',2)
% Title('Cold side outlet temperatures');
legend('T11','T1','T2','T3','Tend')
xlabel('Time [min]');
ylabel('Temperature [\circC]')
axis([0 15 80 140]);

figure
plot(time,Th11,time,Th1,time,Th2,time,Th3)
% Title('Hot side outlet temperatures');
legend('Th11','Th1','Th2','Th3')
xlabel('Time [min]');
ylabel('Temperature [\circC]')

figure
plot(time,u1,time,u2,time,u3,'LineWidth',2);
% Title('Mass fractions');
axis([0 10 0 0.7]);
xlabel('Time [min]'); ylabel('Splits');
legend('u_1', 'u_2', 'u_3');
grid on;

figure
plot(time,c1,time,c2,'LineWidth',2);
xlabel('Time [min]'); ylabel('Controlled variables');
% axis([0 time(end) -1 1]);
legend('c_1', 'c_2');
grid on;

%% Tuning
% close all;
% tstart=28;
% tslutt=50;
% y1start=30;
% y1slutt=46;
% sprangy1=5;
% y2start=1.47;
% y2slutt=3;
% sprangy2=0.2;
% 
% figure
% [AX,H1,H2]=plotyy(time,c1,time,u1*m0);
% set(get(AX(1),'Ylabel'),'String','Controller variable c_1 [\circC]','fontsize',12) 
% set(get(AX(2),'Ylabel'),'String', 'Mass flow [kg/s]','fontsize',12)
% axis(AX(1),[tstart tslutt y1start y1slutt]);
% axis(AX(2),[tstart tslutt y2start y2slutt]);
% set(AX(1),'YLim',[y1start y1slutt])
% set(AX(1),'YTick',y1start:sprangy1:y1slutt)
% set(AX(2),'YLim',[y2start y2slutt])
% set(AX(2),'YTick',y2start:sprangy2:y2slutt)
% set(H1,'linewidth',2)
% set(H2,'linewidth',2)
% xlabel('Time [min]','fontsize',12);
% grid on;

% tstart=25;
% y1start=-6;
% y1slutt=8;
% y2start=7.2;
% y2slutt=8.66;
% figure
% [AX,H1,H2]=plotyy(time,c2,time,u3*m0);
% set(get(AX(1),'Ylabel'),'String','Controller variable c_2 [\circC]','fontsize',12) 
% set(get(AX(2),'Ylabel'),'String', 'Mass flow [kg/s]','fontsize',12)
% axis(AX(1),[tstart time(end) y1start y1slutt]);
% axis(AX(2),[tstart time(end) y2start y2slutt]);
% set(AX(1),'YLim',[y1start y1slutt])
% set(AX(1),'YTick',y1start:2:y1slutt)
% set(AX(2),'YLim',[y2start y2slutt])
% set(AX(2),'YTick',y2start:0.2:y2slutt)
% set(H1,'linewidth',2)
% set(H2,'linewidth',2)
% xlabel('Time [min]','fontsize',12);
% grid on;

