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

clc;clear;close all
%% Parameters
V    = 5000;                                            % tank total volume
rho  = 1000;                                                      % density
m0   = 15*1e3;                                % [kg] reference network mass
eps  = 0.1;                                   % [-] compressibility coeff
p0   = 5;                                              % reference pressure
qnom = 500; 
param = [V;rho;p0;m0;eps;];
%% Initial conditions differential states
Vh0  = 2500;                                        % [m3] hot layer volume
x0   = [Vh0;p0];
%% Initial conditions disturbances
qS0  = qnom;                           % [m3/h] hot water from swing boiler
d0   = [qS0;qnom];                        % vector of disturbances
%% 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
u0   = [qSP0; qTP0; qEB0];                           % vector of inputs

%% Simulation parameters
nx = length(x0); nu = length(u0); nd = length(d0); np = length(param);   % variables number

tsim   = 100/3600;                                          % [hr] simulation time
tustep = 2/3600; % [hr] step change input
tdstep = 2/3600; % [hr] step change disturbance

%% Simulate step on qS (hot water from boiler)
uf     = [qSP0;qTP0+0; qEB0+0];               % step change in inputs
df     = [qS0;qnom*1.5];                  % step change in disturbance

sim('simTES_OL')
%save('ResultsOL_qSP','tOL','yOL','uOL','dOL')
%save('ResultsOL_qS','tOL','yOL','uOL','dOL')
%save('ResultsOL_qTP','tOL','yOL','uOL','dOL')
%save('ResultsOL_q','tOL','yOL','uOL','dOL')

figure;plot(tOL,yOL(:,1)); legend('Vh')
figure;plot(tOL,yOL(:,2)); legend('p')


