# Comparison of Various Model Identification Methods

This example shows several identification methods available in System Identification Toolbox™. We begin by simulating experimental data and use several estimation techniques to estimate models from the data. The following estimation routines are illustrated in this example: `spa`, `ssest`, `tfest`, `arx`, `oe`, `armax` and `bj`.

### System Description

A noise corrupted linear system can be described by:

`y = G u + H e`

where `y` is the output and `u` is the input, while `e` denotes the unmeasured (white) noise source. `G` is the system's transfer function and `H` gives the description of the noise disturbance.

There are many ways to estimate `G` and `H`. Essentially they correspond to different ways of parameterizing these functions.

### Defining a Model

System Identification Toolbox provides users with the option of simulating data as would have been obtained from a physical process. Let us simulate the data from an IDPOLY model with a given set of coefficients.

```B = [0 1 0.5]; A = [1 -1.5 0.7]; m0 = idpoly(A,B,[1 -1 0.2],'Ts',0.25,'Variable','q^-1'); % The sample time is 0.25 s.```

The third argument `[1 -1 0.2]` gives a characterization of the disturbances that act on the system. Execution of these commands produces the following discrete-time idpoly model:

`m0`
```m0 = Discrete-time ARMAX model: A(q)y(t) = B(q)u(t) + C(q)e(t) A(q) = 1 - 1.5 q^-1 + 0.7 q^-2 B(q) = q^-1 + 0.5 q^-2 C(q) = 1 - q^-1 + 0.2 q^-2 Sample time: 0.25 seconds Parameterization: Polynomial orders: na=2 nb=2 nc=2 nk=1 Number of free coefficients: 6 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Created by direct construction or transformation. Not estimated. ```

Here `q` denotes the shift operator so that A(q)y(t) is really short for y(t) - 1.5 y(t-1) + 0.7 y(t-2). The display of model `m0` shows that it is an ARMAX model.

This particular model structure is known as an ARMAX model, where AR (Autoregressive) refers to the A-polynomial, MA (Moving average) to the noise C-polynomial and X to the "eXtra" input B(q)u(t).

In terms of the general transfer functions `G` and `H`, this model corresponds to a parameterization:

`G(q) = B(q)/A(q)` and `H(q) = C(q)/A(q)`, with common denominators

### Simulating the Model

We generate an input signal `u` and simulate the response of the model to these inputs. The `idinput` command can be used to create an input signal to the model and the `iddata` command will package the signal into a suitable format. The `sim` command can then be used to simulate the output as shown below:

```prevRng = rng(12,'v5normal'); u = idinput(350,'rbs'); %Generates a random binary signal of length 350 u = iddata([],u,0.25); %Creates an IDDATA object. The sample time is 0.25 sec. y = sim(m0,u,'noise') %Simulates the model's response for this```
```y = Time domain data set with 350 samples. Sample time: 0.25 seconds Outputs Unit (if specified) y1 ```
``` %input with added white Gaussian noise e %according to the model description m0 rng(prevRng);```

One aspect to note is that the signals `u` and `y` are both IDDATA objects. Next we collect this input-output data to form a single iddata object.

`z = [y,u];`

To get information on the data object (which now incorporates both the input and output data samples), just enter its name at the MATLAB® command window:

`z`
```z = Time domain data set with 350 samples. Sample time: 0.25 seconds Outputs Unit (if specified) y1 Inputs Unit (if specified) u1 ```

To plot the first 100 values of the input `u` and output `y`, use the `plot` command on the iddata object:

```h = gcf; set(h,'DefaultLegendLocation','best') h.Position = [100 100 780 520]; plot(z(1:100));``` It is good practice to use only a portion of the data for estimation purposes, `ze` and save another part to validate the estimated models later on:

```ze = z(1:200); zv = z(201:350);```

### Performing Spectral Analysis

Now that we have obtained the simulated data, we can estimate models and make comparisons. Let us start with spectral analysis. We invoke the `spa` command which returns a spectral analysis estimate of the frequency function and the noise spectrum.

`GS = spa(ze);`

The result of spectral analysis is a complex-valued function of frequency: the frequency response. It is packaged into an object called `IDFRD` object (Identified Frequency Response Data):

`GS`
```GS = IDFRD model. Contains Frequency Response Data for 1 output(s) and 1 input(s), and the spectra for disturbances at the outputs. Response data and disturbance spectra are available at 128 frequency points, ranging from 0.09817 rad/s to 12.57 rad/s. Sample time: 0.25 seconds Output channels: 'y1' Input channels: 'u1' Status: Estimated using SPA on time domain data "ze". ```

To plot the frequency response it is customary to use the Bode plot, with the `bode` or `bodeplot` command:

```clf h = bodeplot(GS); % bodeplot returns a plot handle, which bode does not ax = axis; axis([0.1 10 ax(3:4)])``` The estimate `GS` is uncertain, since it is formed from noise corrupted data. The spectral analysis method provides also its own assessment of this uncertainty. To display the uncertainty (say for example 3 standard deviations) we can use the `showConfidence` command on the plot handle `h` returned by the previous `bodeplot` command.

`showConfidence(h,3)` The plot says that although the nominal estimate of the frequency response (blue line) is not necessarily accurate, there is a 99.7% probability (3 standard deviations of the normal distribution) that the true response is within the shaded (light-blue) region.

### Estimating Parametric State Space Models

Next let us estimate a default (without us specifying a particular model structure) linear model. It will be computed as a state-space model using a prediction error method. We use the `ssest` function in this case:

`m = ssest(ze) % The order of the model will be chosen automatically`
```m = Continuous-time identified state-space model: dx/dt = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x1 -0.007167 1.743 x2 -2.181 -1.337 B = u1 x1 0.09388 x2 0.2408 C = x1 x2 y1 47.34 -14.4 D = u1 y1 0 K = y1 x1 0.04108 x2 -0.03751 Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: none Disturbance component: estimate Number of free coefficients: 10 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using SSEST on time domain data "ze". Fit to estimation data: 75.08% (prediction focus) FPE: 0.9835, MSE: 0.9262 ```

At this point the state-space model matrices do not tell us very much. We could compare the frequency response of the model `m` with that of the spectral analysis estimate `GS`. Note that `GS` is a discrete-time model while `m` is a continuous-time model. It is possible to use `ssest` to compute discrete time models too.

```bodeplot(m,GS) ax = axis; axis([0.1 10 ax(3:4)])``` We note that the two frequency responses are close, despite the fact that they have been estimated in very different ways.

### Estimating Simple Transfer Functions

There is a wide variety of linear model structures, all corresponding to linear difference or differential equations describing the relation between u and y. The different structures all correspond to various ways of modeling the noise influence. The simplest of these models are the transfer function and ARX models. A transfer function takes the form:

`Y(s) = NUM(s)/DEN(s) U(s) + E(s)`

where NUM(s) and DEN(s) are the numerator and denominator polynomials, and Y(s), U(s) and E(s) are the Laplace Transforms of the output, input and error signals (y(t), u(t) and e(t)) respectively. NUM and DEN can be parameterized by their orders which represent the number of zeros and poles. For a given number of poles and zeros, we can estimate a transfer function using the `tfest` command:

`mtf = tfest(ze, 2, 2) % transfer function with 2 zeros and 2 poles`
```mtf = From input "u1" to output "y1": -0.05428 s^2 - 0.02386 s + 29.6 ------------------------------- s^2 + 1.361 s + 4.102 Continuous-time identified transfer function. Parameterization: Number of poles: 2 Number of zeros: 2 Number of free coefficients: 5 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on time domain data "ze". Fit to estimation data: 71.69% FPE: 1.282, MSE: 1.195 ```

AN ARX model is represented as: A(q)y(t) = B(q)u(t) + e(t),

or in long-hand,

`$\mathit{y}\left(\mathit{t}\right)+{\mathit{a}}_{1}^{}\mathit{y}\left(\mathit{t}-1\right)+{\cdots \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\mathit{a}}_{\mathrm{na}}^{}\mathit{y}\left(\mathit{t}-\mathit{na}\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\mathit{b}}_{1}^{}\mathit{u}\left(\mathit{t}-\mathit{nk}\right)+\text{\hspace{0.17em}}\cdots \text{\hspace{0.17em}}{\mathit{b}}_{\mathrm{nb}}^{}\mathit{u}\left(\mathit{t}-\mathit{nb}\text{\hspace{0.17em}}-\mathit{nk}+1\right)+\mathit{e}\left(\mathit{t}\right)\text{\hspace{0.17em}}$`

This model is linear in coefficients. The coefficients ${\mathit{a}}_{1}^{}$, ${\mathit{a}}_{2}^{}$, ..., ${\mathit{b}}_{1}^{}$, ${\mathit{b}}_{2}^{}$, .. can be computed efficiently using a least squares estimation technique. An estimate of an ARX model with a prescribed structure - two poles, one zero and a single lag on the input is obtained as ([na nb nk] = [2 2 1]):

`mx = arx(ze,[2 2 1])`
```mx = Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t) A(z) = 1 - 1.32 z^-1 + 0.5393 z^-2 B(z) = 0.9817 z^-1 + 0.4049 z^-2 Sample time: 0.25 seconds Parameterization: Polynomial orders: na=2 nb=2 nk=1 Number of free coefficients: 4 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using ARX on time domain data "ze". Fit to estimation data: 68.18% (prediction focus) FPE: 1.603, MSE: 1.509 ```

### Comparing the Performance of Estimated Models

Are the models (`mtf`, `mx`) better or worse than the default state-space model `m`? One way to find out is to simulate each model (noise free) with the input from the validation data `zv` and compare the corresponding simulated outputs with the measured output in the set `zv`:

`compare(zv,m,mtf,mx)` The fit is the percent of the output variation that is explained by the model. Clearly `m` is the better model, although `mtf` comes close. A discrete-time transfer function can be estimate using either the `tfest` or the `oe` command. The former creates an IDTF model while the latter creates an IDPOLY model of Output-Error structure (y(t) = B(q)/F(q)u(t) + e(t)), but they are mathematically equivalent.

`md1 = tfest(ze,2,2,'Ts',0.25) % two poles and 2 zeros (counted as roots of polynomials in z^-1)`
```md1 = From input "u1" to output "y1": 0.8383 z^-1 + 0.7199 z^-2 ---------------------------- 1 - 1.497 z^-1 + 0.7099 z^-2 Sample time: 0.25 seconds Discrete-time identified transfer function. Parameterization: Number of poles: 2 Number of zeros: 2 Number of free coefficients: 4 Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using TFEST on time domain data "ze". Fit to estimation data: 71.46% FPE: 1.264, MSE: 1.214 ```
`md2 = oe(ze,[2 2 1]) `
```md2 = Discrete-time OE model: y(t) = [B(z)/F(z)]u(t) + e(t) B(z) = 0.8383 z^-1 + 0.7199 z^-2 F(z) = 1 - 1.497 z^-1 + 0.7099 z^-2 Sample time: 0.25 seconds Parameterization: Polynomial orders: nb=2 nf=2 nk=1 Number of free coefficients: 4 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using OE on time domain data "ze". Fit to estimation data: 71.46% FPE: 1.264, MSE: 1.214 ```
`compare(zv, md1, md2)` The models `md1` and `md2` deliver identical fits to the data.

### Residual Analysis

A further way to gain insight into the quality of a model is to compute the "residuals": i.e. that part `e` in `y = Gu + He` that could not be explained by the model. Ideally, these should be uncorrelated with the input and also mutually uncorrelated. The residuals are computed and their correlation properties are displayed by the `resid` command. This function can be used to evaluate the residues both in the time and frequency domains. First let us obtain the residuals for the Output-Error model in the time domain:

`resid(zv,md2) % plots the result of residual analysis` We see that the cross correlation between residuals and input lies in the confidence region, indicating that there is no significant correlation. The estimate of the dynamics `G` should then be considered as adequate. However, the (auto) correlation of `e` is significant, so `e` cannot be seen as white noise. This means that the noise model `H` is not adequate.

### Estimating ARMAX and Box-Jenkins Models

Let us now compute a second order ARMAX model and a second order Box-Jenkins model. The ARMAX model has the same noise characteristics as the simulated model `m0` and the Box-Jenkins model allows a more general noise description.

`am2 = armax(ze,[2 2 2 1]) % 2nd order ARMAX model`
```am2 = Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t) A(z) = 1 - 1.516 z^-1 + 0.7145 z^-2 B(z) = 0.982 z^-1 + 0.5091 z^-2 C(z) = 1 - 0.9762 z^-1 + 0.218 z^-2 Sample time: 0.25 seconds Parameterization: Polynomial orders: na=2 nb=2 nc=2 nk=1 Number of free coefficients: 6 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using ARMAX on time domain data "ze". Fit to estimation data: 75.08% (prediction focus) FPE: 0.9837, MSE: 0.9264 ```
`bj2 = bj(ze,[2 2 2 2 1]) % 2nd order BOX-JENKINS model`
```bj2 = Discrete-time BJ model: y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t) B(z) = 0.9922 z^-1 + 0.4701 z^-2 C(z) = 1 - 0.6283 z^-1 - 0.1221 z^-2 D(z) = 1 - 1.221 z^-1 + 0.3798 z^-2 F(z) = 1 - 1.522 z^-1 + 0.7243 z^-2 Sample time: 0.25 seconds Parameterization: Polynomial orders: nb=2 nc=2 nd=2 nf=2 nk=1 Number of free coefficients: 8 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using BJ on time domain data "ze". Fit to estimation data: 75.47% (prediction focus) FPE: 0.9722, MSE: 0.8974 ```

### Comparing Estimated Models - Simulation and Prediction Behavior

Now that we have estimated so many different models, let us make some model comparisons. This can be done by comparing the simulated outputs as before:

```clf compare(zv,am2,md2,bj2,mx,mtf,m)``` We can also compare how well the various estimated models are able to predict the output, say, 1 step ahead:

`compare(zv,am2,md2,bj2,mx,mtf,m,1)` The residual analysis of model `bj2` shows that the prediction error is devoid of any information - it is not auto-correlated and also uncorrelated with the input. This shows that the extra dynamic elements of a Box-Jenkins model (the C and D polynomials) were able to capture the noise dynamics.

`resid(zv,bj2)` ### Comparing Frequency Functions

In order to compare the frequency functions for the generated models we use again the `bodeplot` command:

```clf opt = bodeoptions; opt.PhaseMatching = 'on'; w = linspace(0.1,4*pi,100); bodeplot(GS,m,mx,mtf,md2,am2,bj2,w,opt); legend({'SPA','SS','ARX','TF','OE','ARMAX','BJ'},'Location','West');``` The noise spectra of the estimated models can also be analyzed. For example here we compare the noise spectra the ARMAX and the Box-Jenkins models with the state-space and the Spectral Analysis models. For this, we use the `spectrum` command:

```spectrum(GS,m,mx,am2,bj2,w) legend('SPA','SS','ARX','ARMAX','BJ');``` ### Comparing Estimated Models with the True System

Here we validate the estimated models against the true system `m0` that we used to simulate the input and output data. Let us compare the frequency functions of the ARMAX model with the true system.

`bode(m0,am2,w)` The responses practically coincide. The noise spectra can also be compared.

`spectrum(m0,am2,w)` Let us also examine the pole-zero plot:

`h = iopzplot(m0,am2);` It can be seen that the poles and the zeros of the true system (blue) and the ARMAX model (green) are very close.

We can also evaluate the uncertainty of the zeros and poles. To plot confidence regions around the estimated poles and zeros corresponding to 3 standard deviations, use `showConfidence` or turn on the "Confidence Region" characteristic from the plot's context (right-click) menu.

`showConfidence(h,3)` `%`

We see that the true, blue zeros and poles are well inside the green uncertainty regions.