% CSTR
% Transformed Inputs
% u = Q; y=c_A
% Chain of transforms
% Callum Kingstree

warning off

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;       %not required, keep in place so that p vector not distrupted

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;                    %value from nominal df/dc_A, where f = dc_Adt
B = -0.1874;                      %value from nominal dg/dT, where g = dTdt

p = [p; A; B];

% 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+A)*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 = 3;                                             %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= 7.32;                              %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)


%% 3. Simulation

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

T_1step=0;% 17.5;                                   %inlet temp step change [K]
c_Afstep= 0;% 0.1;                 %inlet concentration step change [kmol/m3]
q_1step= 0;%0.1;                              %inlet flow step change [m3/min]

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

sim('Model_chain.slx');                                        %run sim


%% 4. Plotting


c_A = out(:,1);                 %export concentration from output [kmol/m3]
T = out(:,2);                           %export temperature from output [K]

% subplot(121)
% plot(t,c_A,t,y_sp,'k--')
% ylim([0.5 1])
% ylabel('Concentration [kmol/m3]')
% xlabel('Time [min]')
% legend('c_A','c_Asp')
% 
% subplot(122)
% plot(t,Q)
% ylim([0 6000])
% ylabel('Heat Added [kJ/min]')
% xlabel('Time [min]')
%legend('Q')

% filename = 'NO_Chain_sp';                         %sa
% save(filename,'t','c_A','T','Q','v','v_sp','d','y_sp')

subplot(121)
plot(t,c_A)
set(gca,'TickLabelInterpreter','latex')
ylabel('c_A')
xlim([0 100])

subplot(122)
plot(t,Q)
ylabel('Q')
xlim([0 100])