% CSTR General transform
% Transformed Inputs
% u = Q; y=c_A
% Callum Kingstree

warning off
clear all
close all
clc

%% 1. Initialisation

% Model Parameters
T_ref = 400;                          %reference temp for rate constant [K]
V = 4;                                                    %tank volume [m3]
c_p = 0.150;                               %specific heat capacity [kJ/kgK]
E = 20000;                                     %activation energy [kJ/kmol]
rho = 1000;                                                %density [kg/m3]
k_ref = 0.5;                               %reference rate constant [min-1]
R = 8.314;                               %universal gas constant [kJ/Kkmol]
H_rx = -1875;                                   %heat of reaction [kJ/kmol]
q_1 = 1; %inlet flow [m3/min]

p = [T_ref; V; c_p; E; rho; k_ref; R; H_rx; q_1];     %vector of parameters
    
% Nominal Values

T_1nom = 350;                                       %nominal inlet temp [K]
Q_nom = 5000;                                   %nominal heat flow [kJ/min]
c_Afnom= 2;                          %nominal inlet concentration [kmol/m3]
A = -0.75;

% Initial Conditions for States (x_0)

c_A0 = 2/3;                         %initial outlet concentration [kmol/m3]
T_0 = 400;                                           %initial tank temp [K]

x_0 = [c_A0; T_0];                            %vector of initial conditions

% Initial Value for Transformed Input (v)

k_nom = k_ref*exp((-E/R)*((1/T_0)-(1/T_ref)));

v_nom = (q_1/V)*c_Afnom - ((q_1/V)+k_nom-0.75)*c_A0;


% Initial Value for Second Transformed Input (v2)

v2_nom = (q_1/V)*(T_1nom-T_0)-(H_rx*k_nom*c_A0)/(rho*c_p)+...
    (Q_nom/(V*rho*c_p))+0.1874*T_0;

% Number of variables and model parameters

Nx = numel(x_0);                                            %number of ODEs
Nu = 1;                                                      %number of MVs
Nv = 1;                                       %number of transformed inputs
Nd = 2;                                             %number of disturbances
Np = numel(p);                                  %number of input parameters

%% 2. Controller Parameters

%Main Controller (2) (y error to v signal)

k_2 = -1/A;                                         %gain of y(s)/v(s)[min]
tauc_2 = 10;                                              %tuning parameter
theta_2 = 0;                                                    %time delay
tau_2= -1/A;                              %time constant of y(s)/v(s) [min]

Kc_2 = (1/k_2)*(tau_2/(tauc_2+theta_2));          %Main Controller Kc(SIMC)
tauI_2 = min(tau_2,4*(tauc_2+theta_2));       %Main Controller tau_I (SIMC)

%Slave Controller (1) (v error to x2 signal)

k_1 = 1/0.1874;                                  %from second linearisation
tauc_1 = 5;                                               %tuning parameter
theta_1 = 0;                                                    %time delay
tau_1= 1/0.1874;                                   %from fitted OL response

Kc_1 = (1/k_1)*(tau_1/(tauc_1+theta_1));         %Slave Controller Kc(SIMC)
tauI_1 = min(tau_1,4*(tauc_1+theta_1));      %Slave Controller tau_I (SIMC)

%% 3. Simulation

T_1tstep= 10;                         %time of inlet temp step change [min]
c_Aftstep= 40;               %time of inlet concentration step change [min]

T_1step= 10;                                    %inlet temp step change [K]
c_Afstep= 0.1;                   %inlet concentration step change [kmol/m3]

c_A_sp = 2/3;                    %Initial concentration set-point [kmol/m3]
c_A_tsp= 70;                               %Time for set-point change [min]
c_A_newsp= 0.8;                       %New concentration set-pomt [kmol/m3]

D_gain =1 ;

sim('CSTR_GL.slx');


%% 4. Plotting

t=ans.t;
c_A = ans.out(:,1);
T = ans.out(:,2);
Q = ans.simout;
y_sp = ans.y;
iae_n = numel(ans.iae);
iae = ans.iae(iae_n)

D_gain =1.5 ;

sim('CSTR_GL.slx');


%% 4. Plotting

t_2=ans.t;
c_A_2 = ans.out(:,1);
T_2 = ans.out(:,2);
Q_2 = ans.simout;
y_sp_2 = ans.y;
iae_n_2 = numel(ans.iae);
iae = ans.iae(iae_n_2)

D_gain =2 ;

sim('CSTR_GL.slx');


%% 4. Plotting

t_3=ans.t;
c_A_3 = ans.out(:,1);
T_3 = ans.out(:,2);
Q_3 = ans.simout;
y_sp_3 = ans.y;
iae_n_3 = numel(ans.iae);
iae = ans.iae(iae_n_3)

subplot(121)
plot(t,c_A,t_2,c_A_2,t_3,c_A_3,t,y_sp,'k--')
ylim([0.5 1])
ylabel('Concentration [kmol/m3]')
xlabel('Time [min]')
legend('c_A; D=1','c_A; D=1.5','c_A; D=2','c_Asp')

subplot(122)
plot(t,Q,t_2,Q_2,t_3,Q_3)
ylim([0 6000])
ylabel('Heat Added [kJ/min]')
xlabel('Time [min]')
legend('Q; D=1', 'Q; D=1.5','Q;D=2')

