% TKP4140 Process Control Project part 3
warning off
clc
clear 
close all

%% Set default options for plotting
set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',18,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',18,...
'DefaultLineLineWidth',2,...
'DefaultLineMarkerSize',7.75)

%% Tuning parameters for the controllers
% Collecting the gains from the output of part 2
h_gain_q1 = 0.25;
cA_gain_q2 = 0.006089743589744;
% The delay in the transfer functions as shown in the report
theta11 = 0;
theta22 = 2/34;
% Determining the closed loop time constant, as explained in the report
tauc2 = theta22*25;
tauc1 = tauc2/5;
% Calculating Kc and tauI
Kc1 = (1/h_gain_q1)*(1/(tauc1+theta11));
Kc2 = (1/cA_gain_q2)*(1/(tauc2+theta22));
tauI1 = 4*(tauc1+theta11);
tauI2 = 4*(tauc2+theta22);

%% Initializing the model

%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
setpoint_t_step = 1; % Sets the time for changes in the setpoint
d_t_step = 1; % Sets the time for disturbance

%% Setting the setpoints to be nominal values of states
s_h = 1;
s_cA = 1/20;

%% Simulation time
tsim = 20;

%% Testing the model

%% 10% increase in setpoint of h
s_h_step = 0.1*s_h;
s_cA_step = 0;
k_step=0;
caf_step=0;
% Inserting updated disturbance values into step-vectors, which are sent to simulink
d_step = [caf_step;k_step];

%% Running the model
sim('part3sim.mdl')

%% Plotting the result
figure(1)
subplot(421)
plot(time,x(:,1),'red', 'DisplayName','h')
hold on;
plot(time, h_set, '--', 'color', 'black', 'DisplayName', 'h_s')
hold off;
title('Setpoint increase of h')
ylabel('h [m]')
legend()

figure(1)
subplot(423)
plot(time,x(:,2),'red', 'Displayname', 'c_A')
hold on;
plot(time, cA_set, '--', 'color', 'black', 'DisplayName', 'c_{As}')
hold off;
ylabel('cA [kmol/m3]')
legend()

figure(1)
subplot(425)
plot(time,u(:,1),'blue')
hold on;
ylabel('q_1 [m^3/min]')

figure(1)
subplot(427)
plot(time,u(:,2),'blue')
hold on;
ylabel('q_2 [m^3/min]')
xlabel('Time [min]')

IAE_h_h = IAE_h(end);
IAE_cA_h = IAE_cA(end);

%% 10% increase in setpoint of cA
s_h_step = 0;
s_cA_step = 0.1*s_cA;
k_step=0;
caf_step=0;
% Inserting updated disturbance values into step-vectors, which are sent to simulink
d_step = [caf_step;k_step];

%% Running the model
sim('part3sim.mdl')

%% Plotting the result
figure(1)
subplot(422)
plot(time,x(:,1),'red', 'DisplayName','h')
hold on;
plot(time, h_set, '--', 'color', 'black', 'DisplayName', 'h_s')
hold off;
title('Setpoint increase of c_A')
ylabel('h [m]')
legend()

figure(1)
subplot(424)
plot(time,x(:,2),'red', 'Displayname', 'c_A')
hold on;
plot(time, cA_set, '--', 'color', 'black', 'DisplayName', 'c_{As}')
hold off;
ylabel('cA [kmol/m3]')
legend()

figure(1)
subplot(426)
plot(time,u(:,1),'blue')
hold on;
ylabel('q_1 [m^3/min]')

figure(1)
subplot(428)
plot(time,u(:,2),'blue')
hold on;
ylabel('q_2 [m^3/min]')
xlabel('Time [min]')

IAE_h_cA = IAE_h(end);
IAE_cA_cA = IAE_cA(end);

%% 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','SetPointIncrease.eps') % save eps.

%% 10% increase in setpoint of cAf
s_h_step = 0;
s_cA_step = 0;
k_step=0;
caf_step=0.1*caf_nom;
% Inserting updated disturbance values into step-vectors, which are sent to simulink
d_step = [caf_step;k_step];

%% Running the model
sim('part3sim.mdl')

%% Plotting the result
figure(2)
subplot(521)
plot(time,d(:,1),'green')
title('Increase in c_{Af}')
ylabel('c_{Af} [kmol/m^3]')

figure(2)
subplot(523)
plot(time,x(:,1),'red', 'DisplayName','h')
hold on;
plot(time, h_set, '--', 'color', 'black', 'DisplayName', 'h_s')
hold off;
ylabel('h [m]')
legend()

figure(2)
subplot(525)
plot(time,x(:,2),'red', 'Displayname', 'c_A')
hold on;
plot(time, cA_set, '--', 'color', 'black', 'DisplayName', 'c_{As}')
hold off;
ylabel('cA [kmol/m3]')
legend()

figure(2)
subplot(527)
plot(time,u(:,1),'blue')
hold on;
ylabel('q_1 [m^3/min]')

figure(2)
subplot(529)
plot(time,u(:,2),'blue')
hold on;
ylabel('q_2 [m^3/min]')
xlabel('Time [min]')

IAE_h_caf = IAE_h(end);
IAE_cA_caf = IAE_cA(end);

%% 10% increase in setpoint of k
s_h_step = 0;
s_cA_step = 0;
k_step=0.1*k_nom;
caf_step=0;
% Inserting updated disturbance values into step-vectors, which are sent to simulink
d_step = [caf_step;k_step];

%% Running the model
sim('part3sim.mdl')

%% Plotting the result
figure(2)
subplot(522)
plot(time,d(:,2),'green')
ylabel('k_{Af} [m^3/kmol*min]')
title('Increase in k')

figure(2)
subplot(524)
plot(time,x(:,1),'red', 'DisplayName','h')
hold on;
plot(time, h_set, '--', 'color', 'black', 'DisplayName', 'h_s')
hold off;
ylabel('h [m]')
legend()

figure(2)
subplot(526)
plot(time,x(:,2),'red', 'Displayname', 'c_A')
hold on;
plot(time, cA_set, '--', 'color', 'black', 'DisplayName', 'c_{As}')
hold off;
ylabel('cA [kmol/m3]')
legend()

figure(2)
subplot(528)
plot(time,u(:,1),'blue')
hold on;
ylabel('q_1 [m^3/min]')

figure(2)
subplot(5,2,10)
plot(time,u(:,2),'blue')
hold on;
ylabel('q_2 [m^3/min]')
xlabel('Time [min]')

IAE_h_k = IAE_h(end);
IAE_cA_k = IAE_cA(end);

%% 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','DisturbanceControlled.eps') % save eps.

%% Printing the IAE-values
fprintf(['For step change the setpoint of h\n' ...
    ' the IAE_h = %f, and IAE_cA = %f \n'], IAE_h_h, IAE_cA_h)
fprintf(['For step change the setpoint of cA\n' ...
    ' the IAE_h = %f, and IAE_cA = %f \n'], IAE_h_cA, IAE_cA_cA)
fprintf(['For step change in cAf\n' ...
    ' the IAE_h = %f, and IAE_cA = %f \n'], IAE_h_caf, IAE_cA_caf)
fprintf(['For step change in k\n' ...
    ' the IAE_h = %f, and IAE_cA = %f \n'], IAE_h_k, IAE_cA_k)