% Cylla Haase Reactor
%Integrating Transform

% Callum Kingstree

clc
clear all
close all

set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',20,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',18,...
'DefaultLineLineWidth',2,...
'DefaultLineMarkerSize',7.75)

%% 1. Initialization

m_M0 = 0; %initial mass of monomer in tank [kg]
m_P0 = 11.227; %initial mass of polymer in tank [kg]
m_W0 = 42.750; %initial mass of water in tank [kg]
rho_M = 900; %density of monomer [kg/m3]
rho_P = 1040; %density of polymer [kg/m3]
rho_W = 1000; %density of water [kg/m3]
c_pM = 1.675; %heat capacity of monomer [kJ/kgK]
c_pP = 3.14; %heat capacity of polymer [kJ/kgK]
c_pW = 4.187; %heat capacity of water [kJ/kgK]
MW_M = 104; %molecular weight of monomer [kg/kmol]
m_c = 21.455; %mass of coolant in system [kg]
mdot_c = 0.9412; %mass flowrate of coolant [kg/s]
c_pc = 4.187; %heat capacity of coolant [kJ/kgK]
k_0 = 55; %pre-exponential rate constant [s-1]
k_1 = 1000; %rate constant correction [ms/kg]
k_2 = 0.4; %rate constant correction [-]
E = 29560.89; %activation energy [kJ/kmol]
c_0 = 0.000052; %batch viscosity correction [-]
c_1 = 16.4; %batch viscosity correction [-]
c_2 = 2.3; %batch viscosity correction [-]
c_3 = 1.563; %batch viscosity correction [-]
a_0 = 555.556; %batch viscosity correction [K]
del_Hp = 70152.16; %heat of polymerization [kJ/kmol]
d_0 = 0.814; %heat transfer coefficient correction [kW/m2K]
d_1 = -5.13; %heat transfer coefficient correction [ms/kg]
mdot_m = 0.00756; %mass flowrate of monomer - subject to scheduling [kg/s]
t_m0_in = 1800; %start time of monomer flow [s]
t_m1_in = 6000; %end time of monmomer flow [s]
T_set = 355.382; %set-point of reactor temp [K]

R = 8.314; %universal gas constant [kJ/kmolK]
UA_loss = 0.00567567; %heat loss to surroundings coefficient [kW/K]
tau_p = 40.2; %heating/cooling time constant [s]
theta_1= 22.8; %time delay 1 [s]
theta_2= 15; %time delay 2 [s]
i_top = 1.2; %highest impurity factor
i_bot = 0.8; %lowest impurity factor
ff_1 = 0; %fouling factor for first batch [m2K/kW]
ff_2 = 0.176; %fouling factor for second batch [m2K/kW]
ff_3 = 0.352; %fouluing factor for third batch [m2K/kW]
ff_4 = 0.528; %fouling factor for fourth batch [m2K/kW]
ff_5 = 0.704; %fouling factor for fifth batch [m2K/kW]
T_ambW = 280.382; %ambient temperature for winter [K]
T_ambS = 305.382; %ambient temperature for summer [K]
T_inlet = 294.26; %inlet temperature of cooling water [K]
T_steam = 449.82; %inlet temperature of steam [K]
a_3 = 10; %temp set-point parameter [-]
a_4 = -15; %temp set-point parameter [-]
a_5 = 6; %temp set-point parameter [-]
B_1 = 0.193; %reactor bottom area [m2]
P = 1.594; %jacket perimeter [m]
B_2 = 0.167; %jacket bottom area [m2]

real_para = [m_M0; m_P0; m_W0; rho_M;rho_P ;rho_W ;c_pM ;c_pP ;c_pW ;MW_M ;m_c ;...
    mdot_c;c_pc ;k_0;k_1 ;k_2 ;E ;c_0;c_1 ;c_2 ;c_3; a_0;del_Hp;d_0 ;d_1;...
    mdot_m;t_m0_in;t_m1_in ;T_set;R;UA_loss;tau_p ;theta_1;...
    theta_2;...
    i_top ;i_bot ;ff_1;ff_2;ff_3;ff_4 ;ff_5;T_ambW;T_ambS;T_inlet;...
    T_steam;a_3;a_4;a_5;B_1;P;B_2];

changed_para = [m_M0; m_P0; m_W0; rho_M;rho_P ;rho_W ;c_pM ;c_pP ;c_pW ;MW_M ;m_c ;...
    mdot_c;c_pc ;k_0;k_1 ;k_2 ;E ;c_0;c_1 ;c_2 ;c_3; a_0;del_Hp;d_0 ;d_1;...
    mdot_m;t_m0_in;t_m1_in ;T_set;R;UA_loss;tau_p ;theta_1;theta_2;...
    i_top ;i_bot ;ff_1;ff_2;ff_3;ff_4 ;ff_5;T_ambW;T_ambS;T_inlet;...
    T_steam;a_3;a_4;a_5;B_1;P;B_2];

x0 = [m_M0; m_P0; T_ambS];
T0 = [T_ambS; T_ambS];


% Main Controller (T error to vset)

k_main = 1;
tauc_main = 30;
theta_main = 0;

Kc_main = 1/(tauc_main+theta_main);
tauI_main = 4*(tauc_main+theta_main);

% Slave controller 
k_slave = 0.00821;                                            %from fitting
tauc_slave = 475;                                        %tuning parameter
theta_slave = 0;                                                %time delay
tau_slave = 475;                                              %from fitting

Kc_slave = 60000; %(1/k_slave)*(tau_slave/(tauc_slave+theta_slave)); %Slave Controller Kc(SIMC)
tauI_slave = 5000; %min(tau_slave,4*(tauc_slave+theta_slave));    %Slave Controller tau_I (SIMC)

sim('CLTI_simulink.slx');

m_M = simout2(:,1);
m_P = simout2(:,2);
T = simout2(:,3);
error = -T+T_set_in;
% 
% subplot(2,2,1)
% plot (t,T)
% ylim([300 400])
% hold on
% plot(t,T_set_in)
% ylabel('Reactor Temp [K]')
% xlabel('Time [s]')
% legend('T','T_{set}^*')
% %print('-depsc2','-r600','ReactorTemp_int.eps')

% subplot(2,2,2)
% plot(t,valve)
% ylabel('Valve position [%]')
% xlabel('Time [s]')
%print('-depsc2','-r600','valve_pos_int.eps')

% subplot(2,2,3)
% plot(t,Tj)
% xlabel('Time [s]')
% ylabel('Jacket Temp [K]')
%print('-depsc2','-r600','Jackettemp_int.eps')

%subplot(2,2,4)
plot(t,error)
ylim([-0.6 0.6])
xlabel('Time [s]')
ylabel ('Error [K]')
%print('-depsc2','-r600','error_int.eps')

n = numel(iae);
iae = iae(n)

% filename = 'Integrating';                         
% save(filename,'t','v','vs','m_M','m_P','T','Tj','error','valve','T_set_in')

