%% Estimating Transfer Function Models for a Heat Exchanger
%
% The measured data is stored in an excel file heat_exchanger.xlsx and includes measurements of changes in coolant temperature around a nominal and changes in product temperature around a nominal. Estimate a transfer function for the heat exchanger.

dat = readtable("heat_exchanger.xlsx");
head(dat)
t = dat.Time(:);                 % minutes
u = dat.CoolantTemperature(:);   % input (coolant temp change)
y = dat.ProductTemperature(:);   % output (product temp change)

%

figure
plot(t,u,'-','DisplayName','Coolant'), hold on
plot(t,y,'-','DisplayName','Product')
xlabel('Time [min]'), ylabel('Deviation from nominal')
legend, grid on

%

% Sample time --> simple
Ts  = dat.Time(2) - dat.Time(1);  

data = iddata(y,u,Ts,'TimeUnit','min')

sysTF = tfest(data,1,0,nan)

% fit FOPDT
% nd = max(0, round(delayest(dat)));      % delay in samples
% m  = tfest(dat, 1, 0, nd);              % np=1, nz=0, ioDelay=nd*Ts

%% Report the FOPDT parameters
K     = dcgain(sysTF);                      % process gain
tau   = 1/abs(pole(sysTF));                 % time constant [s]
theta = sysTF.IODelay;                      % input–output delay [s]
fprintf('Estimated FOPDT:  G(s) = %g * exp(-%g s) / ( %g s + 1 )\n', K, theta, tau);

%

%% Fit quality
figure
compare(data, sysTF), grid on

%

%% Step and residual analysis
figure, step(sysTF), grid on, title('Estimated step response')
figure, resid(data, sysTF), grid on, title('Residual analysis')

%
%% Transfer Function Estimation from an Initial System

sysInit = idtf(NaN,[1 NaN],'ioDelay',NaN);
sysInit.TimeUnit = 'minutes';

sysInit.Structure.num.Value   = 1;
sysInit.Structure.num.Minimum = 0;
sysInit.Structure.den.Value   = [1 1];
sysInit.Structure.den.Minimum = [0 0];

sysInit.Structure.ioDelay.Value   = 0.2;
sysInit.Structure.ioDelay.Minimum = 0;
sysInit.Structure.ioDelay.Maximum = 1;

sysTF_initialized = tfest(data,sysInit)

compare(data,sysTF,sysTF_initialized)
resid(sysTF_initialized,data);

%
%% Process Model Estimation

sysP1D  = procest(data,'P1D')
compare(data,sysTF,sysTF_initialized, sysP1D)
resid(sysP1D,data);

%
%% Process Model Estimation with Disturbance Model

sysInit = idproc('P1D','TimeUnit','minutes');
sysInit.Structure.Kp.Value    = 1;
sysInit.Structure.Kp.Minimum  = 0;
sysInit.Structure.Tp1.Value   = 1;
sysInit.Structure.Tp1.Maximum = 10;
sysInit.Structure.Td.Value    = 0.2;
sysInit.Structure.Td.Minimum  = 0;
sysInit.Structure.Td.Maximum  = 1;

% specify ARMA1 disturbance model
opt = procestOptions('DisturbanceModel','ARMA1');
sysP1D_noise = procest(data,sysInit,opt)

resid(sysP1D_noise,data);
compare(data,sysTF,sysTF_initialized,sysP1D,sysP1D_noise)