clear all;
close all;
clc;
%% Load data
load('Modeling_Test.mat'); % Load modeling test data
load('Parameters.mat'); % Load values of the parameters of the model (Table 1)

%% Filtering
% In the paper, signals are represented without being filtered (see Fig 12). However, in
% order to obtain a more accurate value of the average errors, they were filtered by means of
% the following code. So, uncomment this to get exactly the same values
% as the ones shown in Table 2:

% setM.Tw_out=Filter(setM.Tw_out);
% setM.Tw_in=Filter(setM.Tw_in);
% setM.Ts_out=Filter(setM.Ts_out);
% setM.Ts_in=Filter(setM.Ts_in);
% setM.Tamb=Filter(setM.Tamb);
% setM.Ta_in=Filter(setM.Ta_in);
% setM.Ta_out=Filter(setM.Ta_out);
% setM.Tt2=Filter(setM.Tt2);

%% Simulation parameters
Ts=0.1; % Sample time
t_sim_i=setM.t(1); % Initial time
t_sim_f=setM.t(end); % Final time

%% Initial contitions for Submodel 3 (SM3)
Tw_out_0=mean(setM.Tw_out(1:10));
Ts_out_0=mean(setM.Ts_out(1:10));
Fw_2_0=mean(setM.Fw_2(1:10))/60000;
Fw_1_0=mean(setM.Fw_1(1:10))/60000;
Ts_in_0=mean(setM.Ts_in(1:10));
Tt_in_0=Tw_out_0;
Tt_out_0=(Fw_2_0/Fw_1_0)*(Ts_in_0-Ts_out_0)+Tt_in_0;
Tt1_0=mean(setM.Tw_in(1:10))+0.6;
Tp1_out_0=mean(setM.Tw_in(1:10));
%% Additional initial contitions for Submodel 1 (SM1)
Ta_out_0=mean(setM.Ta_out(1:10));
Tamb_0=mean(setM.Tamb(1:10));
%% Additional initial contitions for Submodel 2 (SM2)
Tp4_out=setM.Ts_in;
Tp4_out_0=mean(Tp4_out(1:10));
%% Additional initial contitions for Submodel 4 (SM4)
Tt2_0=mean(setM.Tt2(1:10));
Tr_out_0=Ts_out_0;
Tp2_out_0=Ts_out_0;
Tp3_out_0=Ts_out_0;

%% Constants
kw=4112120.21; % [J/(m^3*K)]
Vs=0.18/1000; % [m3]
Vt=0.18/1000; % [m3]
k_plate=15; % [W/m*K]
A_plate=0.172*0.042; % [m^2]
L_plate=0.0005; % [m]
nch=9; % [-]
he=k_plate*A_plate/L_plate; % [W/K]
Fw_1min=min(setM.Fw_1); % [l/min]
Fw_1max=max(setM.Fw_1); % [l/min]
Fw_2min=min(setM.Fw_2); % [l/min]
Fw_2max=max(setM.Fw_2); % [l/min]
nc=16; % [-]
Cfc=12941.5; % [J/K]
Fa_max=164; % [l/min]
Fa_min=120; % [l/min]
D_pipes=0.015; % [m]
Area=pi*(D_pipes/2)^2; % [m^2]
L_p2=4; % [m]
Vp2=L_p2*Area; % [m]
L_p3=2; % [m]
Vp3=L_p3*Area; % [m^3]

%% Simulation
sim('Complete_Model_M.slx');

%% Results of the adjustment of the complete model with the modeling test data.

set(0,'DefaultTextFontSize', 7.5)
set(0,'DefaultAxesFontSize', 7.5)
figure1 = figure;
set(gcf,'Units', 'centimeters')
set(gcf,'Position', [14 6 21 18.1]);

%% subplot1
subplot1 = subplot(6,1,1,'Parent',figure1);
hold(subplot1,'on');
plot1=plot(setM.t, [setM.Tw_out Tw_out]);
set(plot1(2),'LineWidth',0.8);
axis([0 t_sim_f 60 70])
ylabel('T_{w_{out}} (ºC)','FontAngle','italic','FontName','Times New Roman',...
    'Rotation',90);
box(subplot1,'on');
grid(subplot1,'on');
legend1=legend('Plant output', 'Model output');
set(legend1,...
    'Position',[0.738115212996716 0.826022821576622 0.156432746098056 0.044590642142017],...
    'Box','Off',...
    'FontSize',7.5);

%% subplot2
subplot2 = subplot(6,1,2,'Parent',figure1);
hold(subplot2,'on');
plot2=plot(setM.t, [setM.Tw_in Tw_in]);
set(plot2(2),'LineWidth',0.8);
axis([0 t_sim_f 55 65])
ylabel('T_{w_{in}} (ºC)','FontAngle','italic','FontName','Times New Roman',...
    'Rotation',90);
box(subplot2,'on');
grid(subplot2,'on');

%% subplot3
subplot3 = subplot(6,1,3,'Parent',figure1);
hold(subplot3,'on');
plot3=plot(setM.t, [setM.Tt2 Tt2]);
set(plot3(2),'LineWidth',0.8);
axis([0 t_sim_f 50 60])
ylabel('T_{t2} (ºC)','FontAngle','italic','FontName','Times New Roman',...
    'Rotation',90);
box(subplot3,'on');
grid(subplot3,'on');

%% subplot4
subplot4 = subplot(6,1,4,'Parent',figure1);
hold(subplot4,'on');
plot4=plot(setM.t, [setM.Ta_out Ta_out]);
set(plot4(2),'LineWidth',0.8);
axis([0 t_sim_f 55 65])
ylabel('T_{a_{out}} (ºC)','FontAngle','italic','FontName','Times New Roman',...
    'Rotation',90);
box(subplot4,'on');
grid(subplot4,'on');

%% subplot5
subplot5 = subplot(6,1,5,'Parent',figure1);
hold(subplot5,'on');
plot5=plot(setM.t, [setM.Ts_in Ts_in]);
set(plot5(2),'LineWidth',0.8);
axis([0 t_sim_f 50 60])
ylabel('T_{s_{in}} (ºC)','FontAngle','italic','FontName','Times New Roman',...
    'Rotation',90);
box(subplot5,'on');
grid(subplot5,'on');

%% subplot6
subplot6 = subplot(6,1,6,'Parent',figure1);
hold(subplot6,'on');
plot6=plot(setM.t, [setM.Ts_out Ts_out]);
set(plot6(2),'LineWidth',0.8);
axis([0 t_sim_f 55 65])
ylabel('T_{s_{out}} (ºC)','FontAngle','italic','FontName','Times New Roman',...
    'Rotation',90);
xlabel('Time (s)')
box(subplot6,'on');
grid(subplot6,'on');

%% Average error of each output of the model with respect to the outputs of the process
% Table 2
disp(mean(Errors)');


