% TKP4140 Process Control Project
%% Example script showing how to run the Simulink model 
% Example script that:
% - calls the computeTransferFunctions function
% - runs the Simulink model 
% - plots the results

warning off
clc
clear
close all

%% Set default options for plotting. You can change this to your preferences
set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',18,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',18,...
'DefaultLineLineWidth',2,...
'DefaultLineMarkerSize',5)

%% Initializing

%Compute the transfer functions
[G,Gd]=computeTransferFunctions; %Get G and Gd from computeTransferFunctions

%Model parameters
A=4; %area of the tank in m^2
p = [A]; % vertical vector with model parameters. This is used in Simulink

% Nominal  values for u and d.
% These will go to the Workspace and will be used in the "nominal u" and
% "nominal d" blocks.
q1_nom = 1;
q2_nom = 1;
caf_nom = 1;
k_nom = 95;

% The vectors that get sent to the simulink simulation
u_nom=[q1_nom;q2_nom];
d_nom=[caf_nom;k_nom];

% Define the number of variales and model parameters
% These are used in the Interpreted Matlab Fcn block in Simulink
Nx = 2; %--CHANGE HERE-- number of states--NUMBER OF DIFFERENTIAL EQUATIONS
Nu = 2; %--CHANGE HERE-- number of inputs--NUMBER OF MVS
Nd = 2; %--CHANGE HERE-- number of disturbances--NUMBER OF DVS
Np = 1; %--CHANGE HERE-- number of model parameters--CONSTANT VALUES
Ny = Nx;% Do not change--number of outputs---equal to number of states

%INITIAL CONDITIONS FOR THE STATES 
%Set the initial conditions you found for every state.
%These intial conditions will be used to solve the differential equations.
x0 = [1;1/20];  % [h;cA]
  
%% Set the time for giving the steps
u_t_step = 1; % In the "q0 increase" block, the step time will be set to 1.
d_t_step = 1; % In the "q1 increase" block, the step time will be set to 1.

%% Set simulation time
tsim = 20; % same units as your model

%% Step in q1
% Simulating a 10% increase in q1:
q1_step=0.1*q1_nom;
q2_step=0;
k_step=0;
caf_step=0;
% Inserting updated values into step-vectors, which are sent to simulink
u_step = [q1_step;q2_step];
d_step = [caf_step;k_step];

%run Simulink model to test the step
sim('SimulinkPart2') 


figure(1)
subplot(321)
plot(time,u(:,1),'blue')
title('Step on q_1')
ylabel('q_1 [m^3/min]')
xlabel('Time [min]')
ylim([0.995, 1.105]);

figure(1)
subplot(323)
plot(time,x(:,1),'red')
hold on
plot(time,ylinear(:,1),'black', LineStyle='--')
legend('non-linear','linearized','Location','southeast')
ylabel('h [m]')
xlabel('Time [min]')
ylim([0.5,1.5]);
hold off

figure(1)
subplot(325)
plot(time,x(:,2),'red')
hold on
plot(time,ylinear(:,2),'black', LineStyle='--')
legend('non-linear','linearized','Location','northeast')
ylabel('c_A [kmol/m^3]')
ylim([0.04, 0.07]);
xlabel('Time [min]')
hold off

% Finding the slope, integrating process eq(9)
%delta u
du_q1 = 0.1*q1_nom;
%delta h
dh_q1 = x(end,1) - x((end-4),1);
dh_q1_linear = ylinear(end,1) - ylinear((end-4),1);
%delta ca
dcaq1 = x(end,2) - x((end-4),2);
dcaq1_linear = ylinear(end,2) - ylinear((end-4),2);
%delta t
dtq1 = time(end) - time(end-4);

%calculating the gain of h
h_slopeq1 = dh_q1/dtq1;
h_gain_q1 = h_slopeq1/du_q1;
h_slopeq1_linear = dh_q1_linear/dtq1;
h_gain_q1_linear = h_slopeq1_linear/du_q1;

%calculating the gain of ca
ca_slopeq1 = dcaq1/dtq1;
ca_gain_q1 = ca_slopeq1/du_q1;
ca_slopeq1_linear = dcaq1_linear/dtq1;
ca_gain_q1_linear = ca_slopeq1_linear/du_q1;

%Printing the result
disp('For an increase in q1:')
fprintf('The non-linear gain in h is: %1.4f \n', h_gain_q1)
fprintf('The linearized gain (G_11) in h is: %1.4f \n', h_gain_q1_linear)
fprintf('The non-linear gain in cA is: %1.4f \n', ca_gain_q1)
fprintf('The linearized gain (G_21) in h is: %1.4f \n\n', ca_gain_q1_linear)


%% Step in q2
% Simulating a 10% increase in q2:
q1_step=0;
q2_step=0.1*q2_nom;
k_step=0;
caf_step=0;
% Inserting updated values into step-vectors, which are sent to simulink
u_step = [q1_step;q2_step];
d_step = [caf_step;k_step];

%run Simulink model to test the step
sim('SimulinkPart2') 


figure(1)
subplot(322)
plot(time,u(:,2),'blue')
title('Step on q_2')
ylabel('q_1 [m^3/min]')
xlabel('Time [min]')
ylim([0.995, 1.105]);

figure(1)
subplot(324)
plot(time,x(:,1),'red')
hold on
plot(time,ylinear(:,1),'black', LineStyle='--')
legend('non-linear','linearized','Location','northeast')
ylabel('h [m]')
xlabel('Time [min]')
ylim([0.5,1.5]);
hold off

figure(1)
subplot(326)
plot(time,x(:,2),'red')
hold on
plot(time,ylinear(:,2),'black', LineStyle='--')
legend('non-linear','linearized','Location','southeast')
ylabel('c_A [kmol/m^3]')
ylim([0.04, 0.07]);
xlabel('Time [min]')
hold off

figure(1)
print('-depsc2','-r600','Part2_inputs.eps') % save figure as eps

% Finding the slope, integrating process eq(9)
%delta u
du_q2 = 0.1*q2_nom;
%delta h
dh_q2 = x(end,1) - x((end-4),1);
dh_q2_linear = ylinear(end,1) - ylinear((end-4),1);
%delta ca
dcaq2 = x(end,2) - x((end-4),2);
dcaq2_linear = ylinear(end,2) - ylinear((end-4),2);
%delta t
dtq2 = time(end) - time(end-4);

%calculating the gain of h
h_slopeq2 = dh_q2/dtq2;
h_gain_q2 = h_slopeq2/du_q2;
h_slopeq2_linear = dh_q2_linear/dtq2;
h_gain_q2_linear = h_slopeq2_linear/du_q2;

%calculating the gain of ca
ca_slopeq2 = dcaq2/dtq2;
ca_gain_q2 = ca_slopeq2/du_q2;
ca_slopeq2_linear = dcaq2_linear/dtq2;
ca_gain_q2_linear = ca_slopeq2_linear/du_q2;

%Printing the result
disp('For an increase in q2:')
fprintf('The non-linear gain in h is: %1.4f \n', h_gain_q2)
fprintf('The linearized gain (G_12) in h is: %1.4f \n', h_gain_q2_linear)
fprintf('The non-linear gain in cA is: %1.4f \n', ca_gain_q2)
fprintf('The linearized gain (G_22) in ca is: %1.4f \n\n', ca_gain_q2_linear)

%% Step in cAf
% Simulating a 10% increase in cAf:
q1_step=0;
q2_step=0;
k_step=0;
caf_step=0.1*caf_nom;
% Inserting updated values into step-vectors, which are sent to simulink
u_step = [q1_step;q2_step];
d_step = [caf_step;k_step];

%run Simulink model to test the step
sim('SimulinkPart2') 

figure(2)
subplot(321)
plot(time,d(:,1),'blue')
title('Step on c_{Af}')
ylabel('c_{A,f} [kmol/m^3]')
xlabel('Time [min]')
ylim([0.995, 1.105]);

figure(2)
subplot(323)
plot(time,x(:,1),'red')
hold on
plot(time,ylinear(:,1),'black', LineStyle='--')
legend('non-linear','linearized','Location','northwest')
ylabel('h [m]')
xlabel('Time [min]')
hold off

figure(2)
subplot(325)
plot(time,x(:,2),'red')
hold on
plot(time,ylinear(:,2),'black', LineStyle='--')
legend('non-linear','linearized','Location','southeast')
ylabel('c_A [kmol/m^3]')
ylim([0.047, 0.053])
xlabel('Time [min]')
hold off

%Finding the gains, SS-gain, eq(8)
ducaf = 0.1*caf_nom;

%calculating the gains of h
dycaf_h = max(x(:,1)) - min(x(:,1));
dycaf_linear_h = max(ylinear(:,1)) - min(ylinear(:,1));
h_gain_caf = dycaf_h/ducaf;
h_gain_caf_linear = dycaf_linear_h/ducaf;

%calculating the gains of ca
dycaf_ca = max(x(:,2)) - min(x(:,2));
dycaf_linear_ca = max(ylinear(:,2)) - min(ylinear(:,2));
ca_gain_caf = dycaf_ca/ducaf;
ca_gain_caf_linear = dycaf_linear_ca/ducaf;

%Printing the result
disp('For an increase in caf:')
fprintf('The non-linear gain in h is: %1.4f \n', h_gain_caf)
fprintf('The linearized gain (Gd_11) in h is: %1.4f \n', h_gain_caf_linear)
fprintf('The non-linear gain in cA is: %1.4f \n', ca_gain_caf)
fprintf('The linearized gain (Gd_12) in ca is: %1.4f \n\n', ca_gain_caf_linear)

%% Step in k
% Simulating a 10% increase in k:
q1_step=0;
q2_step=0;
k_step=0.1*k_nom;
caf_step=0;
% Inserting updated values into step-vectors, which are sent to simulink
u_step = [q1_step;q2_step];
d_step = [caf_step;k_step];

%run Simulink model to test the step
sim('SimulinkPart2') 


figure(2)
subplot(322)
plot(time,d(:,2),'blue')
title('Step on k')
ylabel('k [m^3/kmol*min]')
xlabel('Time [min]')

figure(2)
subplot(324)
plot(time,x(:,1),'red')
hold on
plot(time,ylinear(:,1),'black', LineStyle='--')
legend('non-linear','linearized','Location','northwest')
ylabel('h [m]')
xlabel('Time [min]')
hold off

figure(2)
subplot(326)
plot(time,x(:,2),'red')
hold on
plot(time,ylinear(:,2),'black', LineStyle='--')
legend('non-linear','linearized','Location','northeast')
ylabel('c_A [kmol/m^3]')
ylim([0.047, 0.053])
xlabel('Time [min]')
hold off

figure(2)
print('-depsc2','-r600','Part2_disturbances.eps') % save figure as eps

%Finding the gains, SS-gain, eq(8)
duk = 0.1*k_nom;

%calculating the gains of h
dyk_h = max(x(:,1)) - min(x(:,1));
dyk_linear_h = max(ylinear(:,1)) - min(ylinear(:,1));
h_gain_k = dyk_h/duk;
h_gain_k_linear = dyk_linear_h/duk;

%calculating the gains of ca
dyk_ca = min(x(:,2)) - max(x(:,2));
dyk_linear_ca = min(ylinear(:,2)) - max(ylinear(:,2));
ca_gain_k = dyk_ca/duk;
ca_gain_k_linear = dyk_linear_ca/duk;

%Printing the result
disp('For an increase in k:')
fprintf('The non-linear gain in h is: %1.4f \n', h_gain_k)
fprintf('The linearized gain (Gd_21) in h is: %1.4f \n', h_gain_k_linear)
fprintf('The non-linear gain in cA is: %1.4e \n', ca_gain_k)
fprintf('The linearized gain (Gd_22) in ca is: %1.4e \n\n', ca_gain_k_linear)
