% Main file to run thermal energy storage system for extended CACE paper
% 08.12.2020
% Cristina Zotica
% time units are in hours

clc;clear;close all
%% Parameters
V    = 5000;                                            % tank total volume
Vmax = V;
rho  = 1000;                                                      % density
p0   = 5;                                              % reference pressure
m0   = 15000;                                 % [kg] reference network mass
eps  = 0.1;                                   % [-] compressibility coeff
qnom = 500;                               % [m3/h] nominal hot water demand
param = [V;rho;p0;m0;eps];
%% Initial conditions differential states
Vh0 = 2500;                                         % [m3] hot layer volume
x0  = [Vh0;p0];                                  % initial condition states
%% Initial conditions disturbances
TS0  = 120;                      % [degC] hot water temperature from boiler
qPT0 = qnom;                      % [m3/h] cold water return from consumers 
qS0  = qnom;                           % [m3/h] hot water from swing boiler

%% Initial conditions inputs
qSP0 = qnom;                         % [m3/h] hot water direct to consumers
qD0  = 0;                                         % [m3/h] hot water dumped
qTP0 = 0;                             % [m3/h] hot water from tank to plant
qEB0 = 0;                           % [m3/h] hot water from electric boiler
qST0 = qS0 - qD0 - qSP0;                % [m3/h] flowrate charging the tank

%% Simulation parameters
nx    = 2; nu = 4; nd = 1; np = length(param);  % variables number

%% Network pressure control. 3 controllers
% dp/dt = p0/m0/eps/rh0*(qSP+qTP+qEB+q). integrating process
kp   = p0/m0/eps*rho;                                       % initial slope
ps   = p0;                          % [bar] pressure setpoint for qSP = MV1
pTPs = ps - 0.5;                 % [bar] pressure setpoint for tank qTP=MV3 
pEPs = pTPs - 0.5;    % [bar] pressure setpoint for electric boiler qEB=MV4

tauCSP = 10/3600; thetaSP=0;        % [h]=10 sec. closed loop time constant
KcSP   = 1/kp*1/(tauCSP+thetaSP);   % proportional gain. qSP (swing boiler)
tauISP = 4*(tauCSP+thetaSP);            % integral time. qSP (swing boiler)
tauTSP = tauISP/2;                                 % tracking time constant

tauCTP = 10/3600; thetaTP=0;        % [h]=10 sec. closed loop time constant
KcTP   = 1/kp*1/(tauCTP+thetaTP);           % proportional gain. qTP (tank)
tauITP = 4*(tauCTP+thetaTP);                    % integral time. qTP (tank)
tauTTP = tauITP/2;                                 % tracking time constant

tauCEP = 10/3600; thetaEP=0;        % [h]=10 sec. closed loop time constant
KcEP   = 1/kp*1/(tauCEP+thetaEP);% proportional gain. qSP (electric boiler)
tauIEP = 4*(tauCEP+thetaEP);         % integral time. qSP (electric boiler)
tauTEP = tauIEP/2;                                 % tracking time constant

%% Network pressure control. Split range control
% dp/dt = p0/m0/eps/rh0*(qSP+qTP+qEB+q). integrating process
qmax  = 1000;                                     % [m3/h] maximum flowrate
N     = 3;                                  % number of inputs given by SRC
u0    = [0;-100;-200];                                         % SRC bias
alpha = [3;3;3];    % slopes in SRC. all MVs have same gain. v=[0;1]
umax  = [qmax;qmax;qmax];
pSR   = [u0;alpha;umax;N];                 % parameters to the split range block

tauCSRC  = 10/3600;%[h]=10sec split range control closed loop time constant
thetaSRC = 0;                                                   % [h] delay
KcSRC    = 1/alpha(1)*1/kp*1/(tauCSRC+thetaSRC);        % proportional gain
tauISRC  = 4*(tauCSRC+thetaSRC);                 %[h] commoun integral time


%% A. Disturbances. constant demand, varying supply
tsim = 30; % [h] simulation time
qsp  = qnom;                               % [m3/h] demand nominal setpoint
tq1  = 12;                          %  [h] step time demand setpoint change
tq2  = 40;
dq   = 500;                                       % [m3/h] step in setpoint

td1  = 1;              % [h] step  time swing boiler availability increases
td2  = 15;             % [h] step  time swing boiler availability decreases
dqS  = 500;          % [m3/h] step in available hot water from swing boiler


%% Simulate split range control
sim('simTES_SRC')
figure;plot(tSRC,ySRC(:,1));legend('Vh')
figure;plot(tSRC,ySRC(:,2));legend('p')
figure; plot(tSRC,qSPSRC); hold on; plot(tSRC,qEPSRC)
% save results
save('ResultsSRC','tSRC','ySRC','qS','qSPSRC','qDSRC','qTPSRC',...
    'qEPSRC','qref','pref')
%% Simulate parallel control
clear time y alg qS qSP qD qTP qEP z qref pref
sim('simTES_3C')
figure;plot(t3C,y3C(:,1));legend('Vh')
figure;plot(t3C,y3C(:,2));legend('p')
figure;plot(t3C,qref);legend('q')

save('Results3C_05','t3C','y3C','qSP3C','qD3C','qTP3C','qEP3C',...,
    'pTPref','pEPref')

save('Results_J_SRC_3C','JSRC','J3C')