% main script for simulating the dynamic tank model, and also the
% control structures. 
% Select file you want to run by in sim('filename')
% 
% Block comment (Ctrl + R); block uncomment (Ctrl + T)
%
% *************************************************************************
% @author: Amos Fang (adapted from Cristina Zotica HEX Temperature SRC)
% @organization: Process Systems Engineering, NTNU
% @project: Master Thesis 2019
% @since: Jan 2019
% @requires: MATLAB R2018b (not tested in other releases)
% @description: Main script for simulating the split-range control
% structure of two tanks
% @estimated run time: N/A
% 
% 
% 
% 
% *************************************************************************
% MATLAB commands
%
% t1       : program start-time
% *************************************************************************
clc
clear

t1=cputime;

%% ************************************************************************
%  Information about the tank
%  
% hmax     : maximum liquid level of the tank [m]
% A        : base tank area [m^2]
% 
% *************************************************************************

hmax = 20;
A = 6000;           % Base Area of Tank

%% ************************************************************************
%  Nominal operating conditions for the tanks
% 
% Fsp      : throughput flow setpoint [-]
% q        : maximum throughput flow rate [m^3/min]
% d1       : nominal hand valve position (z1) [-]
% d2       : nominal hand valve position (z2) [-]
% d3       : nominal hand valve position (z3) [-]
% 
% Define the states. The tank states are assigned x
% x10      : initial height of Tank 1 inventory [m]
% x20      : initial height of Tank 2 inventory [m]
% *************************************************************************

Fsp = 1;
q = 300;

d1 = 1;
d2 = 1; 
d3 = 1;

x10 = 10;
x20 = 10;

%% ************************************************************************
% Simulink Variables
% 
% zmax    : Saturation bias for each input (MVs) of the SRC 
%         :(Typically 0 - fully closed and 1 - fully open) [-]
% d1Tstep : Time when the disturbance is introduced to manual hand valve
%         : z1 [min]
% d2Tstep : Time when the disturbance is introduced to manual hand valve
%         : z2 [min]
% d3Tstep : Time when the disturbance is introduced to manual hand valve
%         : z3 [min]
% interval: Time interval as a multiple of d1Tstep, d2Tstep, d3Tstep
%         : when disturbance step remains [-]
% *************************************************************************
zmax = 1;
d1Tstep = 25;
d2Tstep = 25;
d3Tstep = 25;
interval = 10;

%% ************************************************************************
% Simulink Variables for throughput flow setpoint step change
% 
% FspStep1: 1st throughput flow setpoint step change [% of q]
% FspStep2: 2nd throughput flow setpoint step change [% of q]
% Fsp1Time: Time when 1st throughput flow setpoint change is introduced
% Fsp2Time: Time when 2nd throughput flow setpoint change is introduced
% F1on : Switch for first throughput flow setpoint step change(1:on 2:off)
% F2on : Switch for second throughput flow setpoint step change(1:on 2:off)
% *************************************************************************
            
FspStep1 = -0.5;       
FspStep2 = 0.5;          

Fsp1Time = 18;
Fsp2Time = 37;

F1on = 0;
F2on = 0;

%% ************************************************************************
% Simulink Variables for flow disturbance (moving bottleneck)
%
% d1on     : Switch for disturbance at hand valve z1(1:on 2:off)
% d2on     : Switch for disturbance at hand valve z2(1:on 2:off)
% d3on     : Switch for disturbance at hand valve z3(1:on 2:off)
% dStep    : disturbance step for hand valve closure
% *************************************************************************

d1on = 1;
d2on = 0;
d3on = 0;
dStep = 0.5;

%% ************************************************************************
% Setpoint tracking of the SRC
% r1       : Tank 1 level [m]
% r2       : Tank 2 level [m]

% *************************************************************************

r1 = 10;                 % corresponds to 50% of the Tank 1 level
r2 = 10;                 % corresponds to 50% of the Tank 2 level

%% ************************************************************************
% Simulink Variables for the level split-range controller
% 
% r1delta   : Tank 1 level setpoint step change
% r2delta   : Tank 2 level setpoint step change
% r1on      : Switch for Tank 1 level setpoint change
% r2on      : Switch for Tank 2 level setpoint change
% t1SPstart : Start time for Tank 1 level step change
% t2SPstart : Start time for Tank 2 level step change
% *************************************************************************

r1Delta = 0; %0.1;
r2Delta = 0;
r1on = 0;
r2on = 0;
t1SPstart = 20;
t2SPstart = 20;

%% ************************************************************************
% Controller tunning for SRC
% 
% PI controller tuning parameters calculated with SIMC rules
%
% *************************************************************************

Kc = -1.6;
tauI = 50;

%% ************************************************************************
% MATLAB commands
%
% Run the two-tanks Simulink file; save variables into .mat file
%
% p        : Tank parameters
% t2       : program elapsed-time
% 
% *************************************************************************

p = [A; q];
sim('TwoTanks')
t2 = t1-cputime;

h1 = y1;                 % Tank 1 ODE
h2 = y2;                 % Tank 2 ODE


save('outputTwoTanksSRC.mat','t','h1','h2','q1','q2', 'q3', ...
     'Kc','tauI')

%% ************************************************************************
% Plotting of results
% 
% Figure 1 - subplot of the Tank levels 1 and 2 with time [%]
%          - subplot of hand valve openings with time [%]
% Figure 2 - subplot of SRC signals to control valve q2
%          - subplot of hand valve openings with time [%]
% Figure 3 - subplot of volumetric flow rate with flow disturbance [m^3/min]
%          - subplot of hand valve openings with time [%]
% Figure 4 - subplot of volumetric flow rate with flow throughput setpoint
%            change [m^3/min]
%          - subplot of hand valve openings with time [%]
% 
% *************************************************************************

% *************************************************************************
% Figure 1
%
% This figure displays the plot of the inventory level in both Tank 1
% and 2 when a step change in d1, d2 or d3 is introduced.
% 
% *************************************************************************

figure(1)
set(figure(1),'Color','White')
subplot(211)
plot(t,h1*100/hmax,'r')
hold on
plot(t,h2*100/hmax,'k')
hold off
xlabel('Time,min')
ylabel('Tank Levels, %')
f2 = legend('Tank Level 1','Tank Level 2','Location','Southeast');
set(f2,'FontSize',10)
ax = gca;
ax.FontSize = 12;
grid on
grid minor

subplot(212)
plot(t,z1,'r')
hold on
plot(t,z2,'--k')
hold on
plot(t,z3,'-.b','LineWidth',1)
hold off
xlabel('Time,min')
ylabel('z, hand valve opening, %')
ylim([0.3 1.1]);
legend('z_1','z_2','z_3','Location','Southeast')
ax = gca;
ax.FontSize = 12;
grid on
grid minor

% conditional save
if d1on == 1
    saveas(gcf,'h1-dis','epsc')
elseif d2on == 1
    saveas(gcf,'h2-dis','epsc')
elseif d3on == 1
    saveas(gcf,'h3-dis','epsc')
end

% *************************************************************************
% Figure 2
% 
% This figure displays the plot of the split-range controller signals, x1, 
% x2 or x3 to the control valve q2 when a step change in d1, d2 or d3
% is introduced.
% 
% *************************************************************************

figure(2)
set(figure(2),'Color','White')
plot(t,q2,'Color',[0 0.7 0],'LineWidth',2)
hold on
plot(t,x1,'b')
hold on
plot(t,x2,'--r')
hold on
plot(t,x3,'k')
hold on
ax = gca;
ax.FontSize = 12;
xlabel('Time,min')
ylabel('Valve opening, %')
ylim([0.3 1.1]);
f3=legend('q_2','x_1','x_2','x_3','Location','Southeast');
set(f3,'FontSize',10)
hold off
grid on

% Saving MATLAB plots to .eps
if d1on == 1
    saveas(gcf,'h1-dis-flows','epsc')
elseif d2on == 1
    saveas(gcf,'h2-dis-flows','epsc')
elseif d3on == 1
    saveas(gcf,'h3-dis-flows','epsc')
end

% *************************************************************************
% Figure 3
%
% This figure displays the plot of the flow rate through each control
% valve, q1, q2 or q3 when a step change in d1, d2 or d3 is introduced.
% 
% *************************************************************************

figure(3)
set(figure(3),'Color','White')
subplot(211)
plot(t,v1flow.*z1,'b')
hold on
plot(t,v2flow.*z2,'--r')
hold on
plot(t,v3flow.*z3,'k')
hold on
ax = gca;
ax.FontSize = 12;
xlabel('Time,min')
ylabel('Flow rates m^3/min')
f4 = legend('q_1flow','q_2flow','q_3flow','Location','Southeast');
set(f4,'FontSize',10)
hold off
grid on

subplot(212)
plot(t,z1,'r')
hold on
plot(t,z2,'--k')
hold on
plot(t,z3,'-.b','LineWidth',1)
hold off
xlabel('Time,min')
ylabel('z, hand valve opening, %')
ylim([0.3 1.1]);
legend('z_1','z_2','z_3','Location','Southeast')
ax = gca;
ax.FontSize = 12;
grid on
grid minor

% Saving MATLAB plots to .eps
if d1on == 1
    saveas(gcf,'h1-dis-flows','epsc')
elseif d2on == 1
    saveas(gcf,'h2-dis-flows','epsc')
elseif d3on == 1
    saveas(gcf,'h3-dis-flows','epsc')
end

% *************************************************************************
% Figure 4
%
% Integral wind-up test
% 1) Introduce step flow disturbance to d1 so that the level will deviate
%    from 50%.
% 2) Before Tank 2 level can recover to 50%, introduce a step change to
%    set F1on = 1 (denoted as the variable x3 in the Simulink model
% 3) Leave Fsp = 0 for a period of time
% 4) Reinstate the Fsp to its original setpoint and note the rise in Tank 2
%    ; set F2on = 1
% 5) Apply the anti-windup on the PI controller.
% 
% This figure displays the plot of the inventory level of Tank 1
% and 2 when a step change in flow throughput is introduced, before Tank 2
% level can recover to nominal level. The purpose of this Figure is to
% investigate the integral wind-up effect when PI controllers are used
%
% *************************************************************************

if F1on == 1 || F2on == 1 
    figure(4)
    set(figure(4),'Color','White')
    subplot(211)
    plot(t,h1*100/hmax,'r')
    hold on
    plot(t,h2*100/hmax,'k')
    plot([198 198],[48 52],'r--','LineWidth',1)
    plot([400 400],[48 52],'r--','LineWidth',1)
    text(210,51.5,'Integral wind-up','Interpreter','latex')
    hold off
    xlabel('Time,min')
    ylabel('Tank Levels, %')
    f2 = legend('Tank Level 1','Tank Level 2','Location','Southeast');
    set(f2,'FontSize',10)
    ax = gca;
    ax.FontSize = 12;
    grid on
    grid minor

    subplot(212)
    plot(t,x3*300,'r')
    hold on
    plot(t,z1*300,'b')
    plot([198 198],[0 300],'r--','LineWidth',1)
    plot([408 408],[0 300],'r--','LineWidth',1)
    text(210,270,'Integral wind-up','Interpreter','latex')
    xlabel('Time,min')
    ylabel('Flow rates m^3/min')
    legend('F_{sp}','z_1','Location','Northeast')
    ax = gca;
    ax.FontSize = 12;
    grid on
    grid minor

    % Saving MATLAB plots to .eps
    saveas(gcf,'anti-windup','epsc')
end