% Main script for simulating the dynamic heat exchanger model, and also the
% control structures. 
% Select file you want to run by uncommenting sim('filename')
% Select (comment;uncomment) the correct tuning parameters. 

clc
clear
%close all

N=10;                % Number of cells

F_c = 5;             % cold fluid mass flow, kg/s
F_h = 3;             % hot fluid mass flow,  kg/s
F_c_max = 10;        % maximum cold fluid mass flow, kg/s
F_h_max = 3;         % maximum hot fluid mass flow, kg/s 

z_c = F_c/F_c_max;   % cold valve 
z_h = F_h/F_h_max;   % hot valve

rho_c = 980 ;        % cold fluid density, kg/m3
rho_h = 890;         % hot fluid density,  kg/m3

Tin_h = 70+273.15 ;  % hot fluid inlet temperature, K
U = 150*1e-3;        % heat transfer coefficient,   kW/m2/K
A = 90;              % heat transfer area, m2
cp_c = 1.5;          % specific heat for cold fluid, kJ/kg/K
cp_h = 1.2;          % specific heat for hot fluid,  kJ/kg/K

Tin_c = 20 + 273.15; % cold fluid inlet temperature

u = [z_c; z_h];      % input

p = [rho_c; rho_h; Tin_h; U; A; cp_c; cp_h;N;F_c_max;F_h_max]; % parameters

d = Tin_c;           % disturbances

Tc0 = [314.1411  310.3426  307.0824  304.2841  301.8823  299.8208  298.0514  296.5326  295.2291  294.1102]';
Th0 = [335.2388  328.4495  322.6229  317.6224  313.3307  309.6473  306.4859  303.7724  301.4434  299.4442]';
x0 = [Tc0; Th0]; % Initial conditions
ThSP = x0(end);  % set point for cold valve
ThSP2 = ThSP+1;  % set point for hot valve in the case of two controllers structure

%% Controller tunings for the case of 2 controllers
Kc_c = -0.055;     % controller gain for cold valve. it is the same for input resetting
tauI_c = 74;       % tau I for cold valve. it is the same for input reseting

Kc_h = 0.0803;    % controller gain for hot valve, 2 cotrollers. Set to 0 if this case is not simulated
tauI_h = 86;      % tau I for hot valve, 2 controllers

%% Contrller tunings for input resetting. 
Kc_h_IR = 0.1;        % controller gain for valve position control. 
tauI_h_IR = 77;     % tau I for valve position control

%% Controller tunning for SRC
Kc = -0.06;
tauI = 80;

sim('SRC')                % split range control
%sim('InputResetting')     % input resetting
%sim('TestController')     % test controllers tunings
%sim('Model')              % heat exchanger dynamic model
%sim('TwoControllers')     % separate controllers for each valve

Tc = y(:,1:N);             % cold fluid temperature in each cell
Th = y(:,N+1:2*N);         % hot fluid temperature in each cell

%save('outputSRC.mat','t','Tc','Th','u1','u2','Kc','tauI','ThSP')  % results from SRC sim
%save('outputIR.mat','t','Tc','Th','u1','u2','Kc_c','Kc_h','tauI_c','tauI_h','ThSP')   % results from input resetting
%save('output2C.mat','t','Tc','Th','u1','u2','Kc_c','Kc_h','tauI_c','tauI_h','ThSP','ThSP2')   % results from 2 controllers
%save('outputOL-cold.mat','t','Tc','Th','u1','u2')   % results from open loop simulation
%save('outputOL-hot.mat','t','Tc','Th','u1','u2')   % results from open loop simulation
%save('outputOL-hot-cold.mat','t','u1','u2','Kc_c','tauI_c') % results 


figure 
plot(t,Tc(:,1)-273.15)
hold on
plot(t,Th(:,N)-273.15)
grid on
ax = gca;
ax.FontSize = 16;
xlabel('Time,s')
ylabel('Outlet temperature, ^oC','Interpreter','none')
plot(t,ones(1,length(t))*ThSP-273.15)
legend('cold','hot','Ts')
hold off

figure
plot(t,Th(:,N)-273.15,'r')
xlabel('Time,s')
ylabel('Th_out, ^oC')
ax = gca;
ax.FontSize = 16;
grid on
grid minor

figure
plot(t,u1)
hold on
plot(t,u2)
ax = gca;
ax.FontSize = 16;
xlabel('Time,s')
ylabel('Valve opening')
legend('cold','hot')
hold off
grid on

%% Controller trunig
% % Fitting a first order + delay. (given delay)
% diff = Th(:,N)-Th(1,N); % Difference in the hot outlet temperature
% tol = 0.1;
% idx = find(abs(diff)>tol,1,'first'); % find index when the temperature start changing with a given tolerance
% tstep = 10;
% theta = t(idx)-tstep % delay
% k = (Th(end,N)-Th(1,N))/(u1(end)-u1(1)); % steady-state gain, for change in cold fluid
% %k = (Th(end,N)-Th(1,N))/(u2(end)-u2(1)); % steady-state gain, for change in hot fluid
% ytau = (Th(end,N)-Th(1,N))*0.63+Th(1,N);
% 
% n = [1 2 3 4 5 10 15 20 25]; % tauc = n*theta
% 
% for i = 1:length(n)
%     tauc(i) = n(i)*theta;
%     tau1 = 128.6-theta-tstep;% +10% F_c, tol =0.1
%     %tau1 = 102.8-theta-tstep;  % +10% F_h
%     Kc(i) = 1/k*tau1/(tauc(i)+theta);
%     tauI(i) = min (tau1,4*(theta+tauc(i)));    
% end
% 


%% Input resetting controller (VPC) tuning. y = u1(cold valve); u = u1(hot valve).

% diff = u1-u1(1);
% tol = 0.01;
% idx = find(abs(diff)>tol,1,'first');
% tstep = 100;
% theta = t(idx)-tstep; % delay
% k = (u1(end)-u1(1))/(u2(end)-u2(1));
% dy = u1(end)-u1(1);
% ytau = dy*0.63+u1(1);
% 
% n = [1 2 3 4 5 10 15 20 25 30 35 40 45 50 55 60 65 70]; % tauc = n*theta
% for i = 1:length(n)
%     tauc(i) = n(i)*theta;
%     tau1 = 188.9454-theta-tstep;% tol=0.001
%     Kc(i) = 1/k*tau1/(tauc(i)+theta);
%     tauI(i) = min (tau1,4*(theta+tauc(i)));    
% end


