clc
clear all
warning off

%% Null space method
open('CSTRandColumnAConnected')
load HelS
H=HelS';                    % Sending Hns to Simulink
csns=-1.0e4;                % Sending csns (setpoint of c) to Simulink

%% Step test in L for finding the tuning parameters for the L controller 
% Making a step in L
set_param('CSTRandColumnAConnected/LT - step','time','500','before','0','after','778/60*0.1')
% Running the simulation
sim('CSTRandColumnAConnected')
% Removing the step in L
set_param('CSTRandColumnAConnected/LT - step','time','500','before','0','after','0')

% Plotting the step response in c versus time
figure (1)
plot(c.time,c.signals.values,'b','LineWidth',2.5) 
% The step was made at time=500 min, and the time interval 500-510 min 
% was considered and assumed as an integating process.
xlim([500 510])
set(gca,'XTick',500:2:510);
xlabel('Time [min]','FontSize',12,'FontWeight','bold')
ylabel('c [-]','FontSize',12,'FontWeight','bold')
title('Step Test in L - Exact Local Method','FontSize',12,'FontWeight','bold')

% Calculating the SIMC tuning parameters for an integrating process
dy=c.signals.values(1378)-c.signals.values(522);
dt=c.time(1378)-c.time(522);
du=778/60*0.1;

kprime=dy/(dt*du)                           % kprime
theta_L=0                                   % time delay in L
tauc_L=5                                    % time constant for L
tauI_L=4*(tauc_L+theta_L)                   % integral time
Kc_L=(1/kprime)*(1/(tauc_L+theta_L))        % Controller gain for L
P_L=Kc_L                                    % Tuning parameter P for controlling L
I_L=Kc_L/tauI_L                             % Tuning parameter I for controlling L
% P=0.0336                                  % Calculated value of P
% I=0.0017                                  % Calculated value of I
break 

%% Controlling L and finding V for different disturbances in F0

% F0 is set to different values in the range F0opt +- 20% and the simulaton
% is runned, followed by storing the value of V when using the L controller 
% to keep (csns-c) constant.  

set_param('CSTRandColumnAConnected/LT - step','time','500','before','0','after','0')
set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.2')
sim('CSTRandColumnAConnected') 
V1=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.175')
sim('CSTRandColumnAConnected') 
V2=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.15')
sim('CSTRandColumnAConnected') 
V3=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.125')
sim('CSTRandColumnAConnected') 
V4=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.10')
sim('CSTRandColumnAConnected') 
V5=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.075')
sim('CSTRandColumnAConnected') 
V6=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.05')
sim('CSTRandColumnAConnected') 
V7=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','-460/60*0.025')
sim('CSTRandColumnAConnected') 
V8=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.0')
sim('CSTRandColumnAConnected') 
V9=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.025')
sim('CSTRandColumnAConnected') 
V10=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.05')
sim('CSTRandColumnAConnected') 
V11=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.075')
sim('CSTRandColumnAConnected') 
V12=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.10')
sim('CSTRandColumnAConnected') 
V13=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.125')
sim('CSTRandColumnAConnected') 
V14=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.15')
sim('CSTRandColumnAConnected') 
V15=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.175')
sim('CSTRandColumnAConnected') 
V16=ws.signals.values(end,11)*60;

set_param('CSTRandColumnAConnected/F0_step','time','200','before','0','after','460/60*0.20')
sim('CSTRandColumnAConnected') 
V17=ws.signals.values(end,11)*60;

VelSim=[V1;V2;V3;V4;V5;V6;V7;V8;V9;V10;V11;V12;V13;V14;V15;V16;V17];
