% Task  : Transformed MVs. level and concentration example
% Date  : 05 July 2019
% Author: CZ
% TODO  : 
% Var   : MVs=[q1;q2]; CVs=[M;c]; DV=[q].  

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

%% Parameters

M0   = 100;                                             % [kg], tank holdup

q0   = 10;                                   % [kg/s], nominal outlet flow
q10  = 5;                                             % [kg/s] inlet flow 1
q20  = 5;                                             % [kg/s] inlet flow 2

taur  = M0/q0;                                    % [s] tank residence time
tau02 = taur*1;                                          % tuning parameter

c0   = 10;                       % [kg/kg, %], nominal outlet concentration
c1   = 5;                  % [kg/kg, %], nominal inlet flow 1 concentration
c2   = 15;                 % [kg/kg, %], nominal inlet flow 2 concentration

v1_0 = q10+q20-q0;                            % =0 transformed input 1 bias
v2_0 = tau02/M0*(q10*c1+q20*c2)+(1-tau02/taur)*c0;%transformed input 2 bias

%% Set up simulation
p       = [c1; c2; tau02];                              % system parameters
x0      = [M0;c0];                      % initial condition for integration
u0      = [q10;q20];                                    % physical MVs bias
du      = [0;0.5];                % step in inputs for open loop simulation
tsim    = 500;                                       % [s], simulation time

tstep_M = 20;                  % [s], step time for holdup setpoint changes
tstep_c = 10;           % [s], step time for concentration setpoint changes
tstep_d = 100;                            % [s], step time for disturbances

dM = 0*M0;                                 % step change for flow setpoint
dc = 0.1*c0;                          % step change for concentration setpoint
dq = 0*q0;                       % step change for outflow

% Model parameters
theta  = 3*taur;                                       % no delay for both loops
k.M    = 1;                                    % process gain from v1 to M
k.c    = 1;                                    % process gain from v2 to c
tau.c  = tau02;                           % time constant for concentration

tauC.M = 1;               % desired closed loop time constant for flow loop
tauC.c = 1;               % desired closed loop time constant for flow loop

% PI-controller for holdup
Kc_M   = 1/k.M/(tauC.M+theta); % integrating process
tauI_M = 4*(tauC.M+theta); 

% PI-controller for concentration
Kc_c   = 1/k.c*tau.c/(tauC.c+theta);               % first-order plus delay
tauI_c = min(tau.c,4*(tauC.c+theta));          

%% Simulation. Find u from a static calculation block
sim('simTransformedMVs')
figure
subplot(2,2,1)
plot(time,M)
ylabel('y_1=M, [kg]')
subplot(2,2,3)
plot(time,u(:,1))
ylabel('u_1=q_1')
xlabel('Time,[s]')

%figure(2)
subplot(2,2,2)
plot(time,c)
xlabel('Time,[s]')
ylabel('y_2=c, [kg/kg]')
subplot(2,2,4)
plot(time,u(:,2))
ylabel('u_2=F_2')
xlabel('Time,[s]')





