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

%% Plotting options
set(0,'DefaultTextFontName','Times',...
     'DefaultTextFontSize',14,...
     'DefaultAxesFontName','Times',...
     'DefaultAxesFontSize',14,...
     'DefaultLineLineWidth',4,...
     'DefaultLineMarkerSize',7.75,...
     'DefaultStairLineWidth',2)

close all
clear all
clc
hotpink = [0.8945    0.4102    0.5898];
%% 1. Initialisation

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


% Nominal Values
cAnom  = 2/3;                              %state,  concentration [kmol/m3]
Tnom   = 400;                                         % state tank temp [K]
x0     = [cAnom; Tnom];                 %vector of state nominal conditions
Qnom   = 5000;                           % input nominal heat flow [kJ/min]
Tfnom  = 350;                         % disturbance, nominal inlet temp [K]
cAfnom = 2;            % disturbance, nominal inlet concentration [kmol/m3]
qnom   = 1;                                     % disturbance flow [m3/min]

A     = -(qnom/V+k0);                                    % tuning parameter
vnom  = -A*cAnom;                % Nominal Value for Transformed Input (vA)

p    = [Tref; V; cp; E; rho; k0; R; Hrx;A];           %vector of parameters

%% 2. Controller Parameters
umax = 3*Qnom; umin=0;
ctrl = 1; % outer controller is off
%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= 3.50;                              %time constant of y(s)/v(s) [min]

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

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

k_1 = -3.7e-05;                              %from fitted OL response
tauc_1 = 2;                                               %tuning parameter
theta_1 = 0;                                                    %time delay
tau_1= 4.26;    

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

%% 3. Simulation
tsim=150;

tTf  = 10;                                  %time of inlet temp step change
tcAf = 10;                         %time of inlet concentration step change
tq   = 10;                                       % time of flow step change
tcA  = 10;                                       % time for setpoint change

dTf  =  0; %17.5;                                   %inlet temp step change
dcAf = 0;%  0.1;                                %inlet concentration step change
dq   = 0;% 0.1;                                        % flow step change
dcA  = 0; %2/15;                               % concentration step change

sim('Model_gen.slx');
err(end)

c_A = x(:,1);
T = x(:,2);
Q = u;
%% Save results to .mat file

% filename = 'NO_Gen_sp';
% save(filename,'t','c_A','T','Q','v','vsp','d','ysp') % add other  variables if needed

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

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


