% Task  : Transformed MVs. Shower example
% Date  : 20 Aug 2019
% Author: CZ
% TODO  : change v2 to include T.  
% Var   : MVs=[F1;F2]; CVs=[F;T]; DVs=[T1;T2].  

clc;clear;close all; warning off
%% Plotting options
set(0,'DefaultTextFontName','Times',...
'DefaultTextFontSize',14,...
'DefaultAxesFontName','Times',...
'DefaultAxesFontSize',14,...
'DefaultLineLineWidth',4,...
'DefaultLineMarkerSize',7.75,...
'DefaultStairLineWidth',2)

%% Colorsgreen = [0 0.5 0];
green    = [0 0.5 0];
pistachio = [166 180 1]/256;
gold     = [0.9805, 0.7773,  0.25];
brick    = [108, 1, 2]/256;
cherry   = [222, 58, 58]/256;%[213, 1, 2]/256;
orange   = [255, 144, 51]/256;
indigo   = [0.2930 0 0.5078];
lilac    = [0.7383    0.5898    0.7930];
hotpink  = [0.8945    0.4102    0.5898];
cerulean = [0.1094    0.6641    0.8359];
plum     = [194, 147, 141]/256;

%% Parameters

m     = 100;                                           % [kg], tank holdup

F0    = 10;                                   % [kg/s], nominal outlet flow
F10   = 5;                                            % [kg/s] inlet flow 1
F20   = 5;                                            % [kg/s] inlet flow 2

T0    = 35;                           % [K], nominal outlet temperature
T10   = 20;                     % [K], nominal inlet flow 1 temperature
T20   = 50;                     % [K], nominal inlet flow 2 temperature

taur0 = m/(F10+F20);                           % [s] nominal residence time
n     = 1; 
tau0  = taur0*n;                                         % tuning parameter

v1_0 = F10+F20;                                       % transf input 1 bias
v2_0 = tau0/m*(F10*T10+F20*T20)+(1-tau0/taur0)*T0;    % transf input 2 bias

u0 = [F10;F20];                                                % input bias
p  = [m; tau0];                                                % parameters
%% Set up simulation
tsim    = 500;                                       % [s], simulation time

% Model parameters
theta  = 0;                                       % no delay for both loops
k.F    = 1;                                    % process gain from u1t to F
k.T    = 1;                                    % process gain from u2t to T
tau.T  = tau0;                               % time constant for temprature

tauC.F = 5;               % desired closed loop time constant for flow loop
tauC.T = 5;               % desired closed loop time constant for flow loop

% I-controller for flow
KI   = 1/k.F/(tauC.F+theta);      

% PI-controller for temprature
Kc   = 1/k.T*tau.T/(tauC.T+theta);      % proportional gain flow controller
tauI = min(tau.T,4*(tauC.T+theta));       % integral time flow controller

%% Simulation no delay.

tstep_F  = 200;                   % [s], step time for flow setpoint changes
tstep_T  = 150;            % [s], step time for temperature setpoint changes
tstep_d1 = 50;                            % [s], step time for disturbances
tstep_d2 = 100;                            % [s], step time for disturbances

dF  = 1;0*F0;                                 % step change for flow setpoint
dT  = 1;0*T0;                          % step change for temperature setpoint
dT1 = 0.15*T10;                       % step change for inflow 1 temperature
dT2 = 0.1*T20;                         % step change for inflow 2 temperature

delay = 0; % no measurment delay

sim('simTransformedMVs')

%% Plot individual figures
figure(1)
plot(time,F,'red')
hold on
plot(time, Fsp,'k:')
ylim([9 11.5])
xlim([0 300])
ylabel('Outlet flow, [kg/s]')
xlabel('Time,[s]')
legend('y_1=q', 'y_1^s=q^s ', 'Location','east')
set(gcf,'Position',[100 100 500 250])
%print('-depsc2','-r600','y1.eps')

figure(2)
plot(time,T,'blue')
hold on 
plot(time, Tsp, ':','Color',hotpink) 
%ylim([34.5 36.5])
ylim([34.9 36.1])
xlim([0 300])
ylabel('Temperature [degC]')
xlabel('Time,[s]')
legend('y_2=T','y_2^s=T^s', 'Location','east') 
set(gcf,'Position',[100 100 500 250]) 
%print('-depsc2','-r600','y2.eps')

figure(3)
plot(time, d(:,1),'Color',green) 
hold on 
plot(time,d(:,2),'Color',lilac) 
ylim([19 56])
xlim([0 300])
ylabel('Temperature [degC]')
xlabel('Time,[s]') 
legend('d_1=T_1','d_2=T_2','Location','east')
set(gcf,'Position',[100 100 500 250])
print('-depsc2','-r600','d.eps')

figure(4)
plot(time, u(:,1),'Color',indigo)
hold on 
plot(time,u(:,2),'Color',cerulean)
ylim([3 7]) 
xlim([0 300])
ylabel('Inflow (input) [kg/s]')
xlabel('Time,[s]')
legend('u_1=q_1','u_2=q_2','Location','northeast')
set(gcf,'Position',[100 100 500 250]) 
%print('-depsc2','-r600','u.eps')

%% Small plots
figure()
subplot(1,2,2)
plot(time,T,'blue')
hold on
plot(time, Tsp, ':','Color',hotpink)
ylim([34.9 36.1])
xlim([0 300])
ylabel('Temperature (^\circ C)')
xlabel('Time (s)')
legend('y_2=T','y_2^s','Location','northwest')
legend boxoff

subplot(1,2,1)
plot(time,F,'Color',brick)
hold on
plot(time, Fsp,':','Color',cherry)
plot(time, u(:,1),'Color',indigo)
plot(time, u(:,2),'Color',cerulean)
ylim([3 11.5])
xlim([0 300])
ylabel('Flow (kg/s)')
xlabel('Time (s)')
lgd=legend('y_1=F','y_1^s','u_1=F_1','u_2=F_2');%,'Location','north');
lgd.Position = [0.1641    0.493    0.2823    0.1];
lgd.NumColumns=2;
legend boxoff
set(gcf,'Position',[100 100 600 250])
%print('-depsc2','-r600','Shower.eps')
