% TKP4140 Process Control Project
%% Example script showing how to run the Simulink model 
%  and how you should plot the results
% This is the script that you should run.
% Update this script with your model
%
warning off
clc
clear 
close all

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

%% Initializing

%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.
% The nominal values are used for 
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.

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

% Setpoints = nominal values of states
setpoint_h = 1;
setpoint_cA = 1/20;



%run Simulink model to test the step in q1
tsim = 10;                  % set simulation time.
sim('exampleSimulink') 

%Plotting the results
figure(1)
subplot(341)
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(345)
plot(time,x(:,1),'red')
ylabel('h [m]')
xlabel('Time [min]')
ylim([0.77,1.23]);

figure(1)
subplot(349)
plot(time,x(:,2),'red')
ylabel('c_A [kmol/m^3]')
ylim([0.047, 0.057])
xlabel('Time [min]')

% 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);
%delta ca
dcaq1 = x(end,2) - x((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;

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

%Printing the result
disp('The gain of h for increase in q1 is:')
disp(h_gain_q1)
disp('The gain of ca for increase in q1 is:')
disp(ca_gain_q1)


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

%Simulating again
sim('exampleSimulink') 

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

figure(1)
subplot(346)
plot(time,x(:,1),'red')
ylabel('h [m]')
xlabel('Time [min]')
ylim([0.77,1.23]);

figure(1)
subplot(3,4,10)
plot(time,x(:,2),'red')
ylabel('c_A [kmol/m^3]')
ylim([0.047, 0.057])
xlabel('Time [min]')

% 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);
%delta ca
dcaq2 = x(end,2) - x((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;

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

%Printing the result
disp('The gain of h for increase in q2 is:')
disp(h_gain_q2)
disp('The gain of ca for increase in q2 is:')
disp(ca_gain_q2)


%% Step in caf
% Simulating a 10% increase in caf:
q1_step=0;
q2_step=0;
caf_step=0.1*caf_nom;
k_step=0;
% Inserting updated values into step-vectors
u_step = [q1_step;q2_step];
d_step = [caf_step;k_step];

%Simulating again
sim('exampleSimulink') 

%Plotting
figure(1)
subplot(343)
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(1)
subplot(347)
plot(time,x(:,1),'red')
ylabel('h [m]')
xlabel('Time [min]')
ylim([0.77,1.23]);
%ylim([9.9 11.1])

figure(1)
subplot(3,4,11)
plot(time,x(:,2),'red')
ylabel('c_A [kmol/m^3]')
ylim([0.047, 0.057])
xlabel('Time [min]')

%Finding the gain of cA (h has zero gain)
ducaf = 0.1*caf_nom;
dycaf = max(x(:,2)) - min(x(:,2));
disp('The SS gain of cA for increase in cAf is:')
disp(dycaf/ducaf)


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

%Simulating again
sim('exampleSimulink') 

%Plotting
figure(1)
subplot(344)
plot(time,d(:,2),'blue')
title('Step on k')
ylabel('k [m^3/kmol*min]')
xlabel('Time [min]')
ylim([94.525, 104.975])

figure(1)
subplot(348)
plot(time,x(:,1),'red')
ylabel('h [m]')
xlabel('Time [min]')
ylim([0.77,1.23]);
%ylim([9.9 11.1])

figure(1)
subplot(3,4,12)
plot(time,x(:,2),'red')
ylabel('c_A [kmol/m^3]')
ylim([0.047, 0.057])
xlabel('Time [min]')

%Calculating the gain
duk = 0.1*k_nom;
dyk = max(x(:,2)) - min(x(:,2));
var = dyk/duk;
disp('The SS gain of cA for increase in k is:')
disp(var)

%% Export figure in .eps format. This gives optimal results in Latex

set(gcf,'Position',[100 100 1200 800])                 % set the figure size 
                       % first 2 numbers are the coordinates on your screen 
                                           % 3rd number is the figure width 
                                           % 4th number is the figure hight
                                           
print('-depsc2','-r600','Resulting_Graph.eps') % save eps.