% CSTR
%no transform
% Callum Kingstree
% u = Q; y=c_A

warning off
close all
clear 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;

p = [p; A];

% 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;


% 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 = -7.4*10^(-3);                                         %gain of y(s)/v(s)[min]
tauc_2 = 10;                                               %tuning parameter
theta_2 = 0;                                                    %time delay
tau_2= 3.65;                              %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 = 0.0073;                              %from fitted OL response
tauc_1 = 2;                                               %tuning parameter
theta_1 = 0;                                                    %time delay
tau_1= 4.27;                                         %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
c_Aftstep= 10; %time of inlet concentration step change
q_1tstep = 10;

T_1step=  0;%17.5; %inlet temp step change
c_Afstep= 0; %0.1; %inlet concentration step change
q_1step= 0; % 0.1;

c_A_sp = 2/3;
c_A_tsp= 10;
c_A_newsp=0.7; % 0.8;

sim('Model_notrans.slx');


%% 4. Plotting


c_A = out(:,1);
T = out(:,2);

subplot(121)
plot(t,c_A)

subplot(122)
plot(t,Q)


% subplot(121)
% plot(t,c_A,t,y_sp,'k--')
% ylabel('Concentration [kmol/m3]')
% xlim([0 100])
% 
% xlabel('Time [min]')
% legend('cA','c_Asp')

filename = 'PI_sp_plot';
save(filename,'t','c_A','T','Q','y_sp','d') % add other  variables if needed
