This example shows how to estimate a transfer function from measured signal data.

**Heat Exchanger**

In this example we estimate the transfer function for a heat exchanger. The heat exchanger consists of a coolant temperature, product temperature, and disturbance ambient temperature. We will estimate the coolant to product temperature transfer function.

**Configuring the Measured Data**

The measured data is stored in a MATLAB file and includes measurements of changes in coolant temperature around a nominal and changes in product temperature around a nominal.

```
load iddemo_heatexchanger_data
```

Collect the measured data using the `iddata`

command and plot it.

data = iddata(pt,ct,Ts); data.InputName = '\Delta CTemp'; data.InputUnit = 'C'; data.OutputName = '\Delta PTemp'; data.OutputUnit = 'C'; data.TimeUnit = 'minutes'; plot(data)

**Transfer Function Estimation**

From the physics of the problem we know that the heat exchanger can be described by a first order system with delay. Use the `tfest`

command specifying one pole, no zeros, and an unknown input/output delay to estimate a transfer function

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

sysTF = From input "\Delta CTemp" to output "\Delta PTemp": 1.468 exp(-0.0483*s) * --------- s + 1.561 Continuous-time identified transfer function. Parameterization: Number of poles: 1 Number of zeros: 0 Number of free coefficients: 2 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on time domain data "data". Fit to estimation data: 36% (simulation focus) FPE: 0.0264, MSE: 0.02622

The `compare`

and `resid`

commands allow us to investigate how well the estimated model matches the measured data.

set(gcf,'DefaultAxesTitleFontSizeMultiplier',1,... 'DefaultAxesTitleFontWeight','normal',... 'Position',[100 100 780 520]); resid(sysTF,data);

compare(data,sysTF)

The figure shows that residuals are strongly correlated implying there is information in the measured data that has not be adequately captured by the estimated model.

**Transfer Function Estimation from an Initial System**

Previously we estimated a transfer function from data but apart for the system order did not include much a priori knowledge. From the physics of the problem we know that the system is stable and has positive gain. Inspecting the measured data we also know that the input/output delay is around 1/5 of a minute. We use this information to create an initial system and estimate a transfer function using this system as the initial guess.

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

Restrict the transfer function numerator and denominator terms so that the system is stable with positive gain.

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

Restrict the input/output delay to the range [1/5 1] minute.

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

Use the system as an initial guess for the estimation problem

sysTF_initialized = tfest(data,sysInit)

sysTF_initialized = From input "\Delta CTemp" to output "\Delta PTemp": 2.94 exp(-0.287*s) * --------- s + 3.302 Continuous-time identified transfer function. Parameterization: Number of poles: 1 Number of zeros: 0 Number of free coefficients: 2 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on time domain data "data". Fit to estimation data: 68.28% (simulation focus) FPE: 0.006168, MSE: 0.006442

resid(sysTF_initialized,data);

compare(data,sysTF,sysTF_initialized)

**Process Model Estimation**

In the above we treated the estimation problem as a transfer function estimation problem, but we know that there is some additional structure we can impose. Specifically the heat exchanger system is known to be a 1st order process with delay, or 'P1D' model of the form:

Use the `procest`

command to further impose this structure on the problem

```
sysP1D = procest(data,'P1D')
```

sysP1D = Process model with transfer function: Kp G(s) = ---------- * exp(-Td*s) 1+Tp1*s Kp = 0.90548 Tp1 = 0.32153 Td = 0.25435 Parameterization: 'P1D' Number of free coefficients: 3 Use "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using PROCEST on time domain data "data". Fit to estimation data: 70.4% FPE: 0.005614, MSE: 0.005607

resid(sysP1D,data);

compare(data,sysTF,sysTF_initialized,sysP1D)

**Process Model Estimation with Disturbance Model**

The residual plots of all the estimations performed so far show that the residual correlation is still high, implying that the model is not rich enough to explain all the information in the measured data. The key missing piece is the disturbance ambient temperature which we have not yet included in the model.

Create a 'P1D' process model with restriction on delay and time constant values and use this as the initial guess for the estimation problem.

sysInit = idproc('P1D','TimeUnit','minutes');

Restrict the model to have positive gain, and delay in the range [1/5 1] minute.

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 the option to use a first-order model ('ARMA1') for the disturbance component. Use the template model sysInit along with the option set to estimate the model.

opt = procestOptions('DisturbanceModel','ARMA1'); sysP1D_noise = procest(data,sysInit,opt)

sysP1D_noise = Process model with transfer function: Kp G(s) = ---------- * exp(-Td*s) 1+Tp1*s Kp = 0.91001 Tp1 = 0.3356 Td = 0.24833 An additive ARMA disturbance model has been estimated y = G u + (C/D)e C(s) = s + 591.6 D(s) = s + 3.217 Parameterization: 'P1D' Number of free coefficients: 5 Use "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using PROCEST on time domain data "data". Fit to estimation data: 96.86% (prediction focus) FPE: 6.307e-05, MSE: 6.294e-05

resid(sysP1D_noise,data);

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

The residual plot clearly shows that the residuals are uncorrelated implying we have a model that explains the measured data. The 'ARMA1' disturbance component we estimated is stored as numerator and denominator values in the "NoiseTF" property of the model.

sysP1D_noise.NoiseTF

ans = num: {[1 591.6038]} den: {[1 3.2172]}

**Compare Different Models**

Although we have identified a model that explains the measured data we note that the model fit to measured data is around 70%. The loss in fit value is a consequence of the strong effect of the ambient temperature disturbance which is illustrated as follows.

The measured data was obtained from a simulink model with the following exact values (d is a gaussian nose disturbance)

y = (1-pi/100)*exp(-15s)/(21.3s+1)*u + 1/(25s+1)*d

Create a 'P1D' model with these values and see how well that model fits the measured data.

sysReal = idproc('P1D','TimeUnit','minutes'); sysReal.Structure.Kp.Value = 1-pi/100; sysReal.Structure.Td.Value = 15/60; sysReal.Structure.Tp1.Value = 21.3/60; sysReal.NoiseTF = struct('num',{[1 10000]},'den',{[1 0.04]});

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

The comparison plot shows that the true system fit to measured data is also around 70% confirming that our estimated models are no better than the real model in fitting measured data - a consequence of the strong effect of the ambient temperature disturbance.

**Additional Information**

For more information on identification of dynamic systems with System Identification Toolbox visit the System Identification Toolbox product information page.

Was this topic helpful?