Main Content

This example shows how to manage data and model objects available in the System Identification Toolbox™. System identification is about building models from data. A data set is characterized by several pieces of information: The input and output signals, the sample time, the variable names and units, etc. Similarly, the estimated models contain information of different kinds - estimated parameters, their covariance matrices, model structure and so on.

This means that it is suitable and desirable to package relevant information around data and models into objects. System Identification Toolbox™ contains a number of such objects, and the basic features of these are described in this example.

First create some data:

u = sign(randn(200,2)); % 2 inputs y = randn(200,1); % 1 output ts = 0.1; % The sample time

To collect the input and the output in one object do

z = iddata(y,u,ts);

The data information is displayed by just typing its name:

z

z = Time domain data set with 200 samples. Sample time: 0.1 seconds Outputs Unit (if specified) y1 Inputs Unit (if specified) u1 u2

The data is plotted as iddata by the `plot`

command, as in `plot(z)`

. Press a key to continue and advance between the subplots. Here, we plot the channels separately:

```
plot(z(:,1,1)) % Data subset with Input 1 and Output 1.
```

```
plot(z(:,1,2)) % Data subset with Input 2 and Output 1.
```

To retrieve the outputs and inputs, use

u = z.u; % or, equivalently u = get(z,'u'); y = z.y; % or, equivalently y = get(z,'y');

To select a portion of the data:

zp = z(48:79);

To select the first output and the second input:

```
zs = z(:,1,2); % The ':' refers to all the data time points.
```

The sub-selections can be combined:

```
plot(z(45:54,1,2)) % samples 45 to 54 of response from second input to the first output.
```

The channels are given default names 'y1', 'u2', etc. This can be changed to any values by

set(z,'InputName',{'Voltage';'Current'},'OutputName','Speed');

Equivalently we could write

z.inputn = {'Voltage';'Current'}; % Autofill is used for properties z.outputn = 'Speed'; % Upper and lower cases are also ignored

For bookkeeping and plots, also units can be set:

z.InputUnit = {'Volt';'Ampere'}; z.OutputUnit = 'm/s'; z

z = Time domain data set with 200 samples. Sample time: 0.1 seconds Outputs Unit (if specified) Speed m/s Inputs Unit (if specified) Voltage Volt Current Ampere

All current properties are (as for any object) obtained by get:

get(z)

ans = struct with fields: Domain: 'Time' Name: '' OutputData: [200x1 double] y: 'Same as OutputData' OutputName: {'Speed'} OutputUnit: {'m/s'} InputData: [200x2 double] u: 'Same as InputData' InputName: {2x1 cell} InputUnit: {2x1 cell} Period: [2x1 double] InterSample: {2x1 cell} Ts: 0.1000 Tstart: [] SamplingInstants: [200x0 double] TimeUnit: 'seconds' ExperimentName: 'Exp1' Notes: {} UserData: []

In addition to the properties discussed so far, we have 'Period' which denotes the period of the input if periodic Period = inf means a non-periodic input:

z.Period

ans = Inf Inf

The intersample behavior of the input may be given as 'zoh' (zero-order-hold, i.e. piecewise constant) or 'foh' (first- order-hold, i.e., piecewise linear). The identification routines use this information to adjust the algorithms.

z.InterSample

ans = 2x1 cell array {'zoh'} {'zoh'}

You can add channels (both input and output) by "horizontal concatenation", i.e. z = [z1 z2]:

z2 = iddata(rand(200,1),ones(200,1),0.1,'OutputName','New Output',... 'InputName','New Input'); z3 = [z,z2]

z3 = Time domain data set with 200 samples. Sample time: 0.1 seconds Outputs Unit (if specified) Speed m/s New Output Inputs Unit (if specified) Voltage Volt Current Ampere New Input

Let us plot some of the channels of `z3`

:

```
plot(z3(:,1,1)) % Data subset with Input 2 and Output 1.
```

```
plot(z3(:,2,3)) % Data subset with Input 2 and Output 3.
```

**Generating Inputs**

The command `idinput`

generates typical input signals.

u = idinput([30 1 10],'sine'); % 10 periods of 30 samples u = iddata([],u,1,'Period',30) % Making the input an IDDATA object.

u = Time domain data set with 300 samples. Sample time: 1 seconds Inputs Unit (if specified) u1

SIM applied to an iddata input delivers an iddata output. Let us use `sim`

to obtain the response of an estimated model `m`

using the input `u`

. We also add noise to the model response in accordance with the noise dynamics of the model. We do this by using the "AddNoise" simulation option:

m = idpoly([1 -1.5 0.7],[0 1 0.5]); % This creates a model; see below. options = simOptions; options.AddNoise = true; y = sim(m,u,options) % simulated response produced as an iddata object

y = Time domain data set with 300 samples. Sample time: 1 seconds Outputs Unit (if specified) y1

The simulation input `u`

and the output `y`

may be combined into a single `iddata`

object as follows:

```
z5 = [y u] % The output-input iddata.
```

z5 = Time domain data set with 300 samples. Sample time: 1 seconds Outputs Unit (if specified) y1 Inputs Unit (if specified) u1

More about the `iddata`

object is found under `help iddata`

.

All models are delivered as MATLAB® objects. There are a few different objects depending on the type of model used, but this is mostly transparent.

load iddata1 m = armax(z1,[2 2 2 1]); % This creates an ARMAX model, delivered as an IDPOLY object

All relevant properties of this model are packaged as one object (here, `idpoly`

). To display it just type its name:

m

m = Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t) A(z) = 1 - 1.531 z^-1 + 0.7293 z^-2 B(z) = 0.943 z^-1 + 0.5224 z^-2 C(z) = 1 - 1.059 z^-1 + 0.1968 z^-2 Sample time: 0.1 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 "z1". Fit to estimation data: 76.38% (prediction focus) FPE: 1.127, MSE: 1.082

Many of the model properties are directly accessible

```
m.a % The A-polynomial
```

ans = 1.0000 -1.5312 0.7293

A list of properties is obtained by get:

get(m)

A: [1 -1.5312 0.7293] B: [0 0.9430 0.5224] C: [1 -1.0587 0.1968] D: 1 F: 1 IntegrateNoise: 0 Variable: 'z^-1' IODelay: 0 Structure: [1x1 pmodel.polynomial] NoiseVariance: 1.1045 InputDelay: 0 OutputDelay: 0 Ts: 0.1000 TimeUnit: 'seconds' InputName: {'u1'} InputUnit: {''} InputGroup: [1x1 struct] OutputName: {'y1'} OutputUnit: {''} OutputGroup: [1x1 struct] Notes: [0x1 string] UserData: [] Name: '' SamplingGrid: [1x1 struct] Report: [1x1 idresults.polyest]

Use `present`

to view the estimated parameter covariance as +/- 1 standard deviation uncertainty values on individual parameters:

present(m)

m = Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t) A(z) = 1 - 1.531 (+/- 0.01801) z^-1 + 0.7293 (+/- 0.01473) z^-2 B(z) = 0.943 (+/- 0.06074) z^-1 + 0.5224 (+/- 0.07818) z^-2 C(z) = 1 - 1.059 (+/- 0.06067) z^-1 + 0.1968 (+/- 0.05957) z^-2 Sample time: 0.1 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: Termination condition: Near (local) minimum, (norm(g) < tol).. Number of iterations: 3, Number of function evaluations: 7 Estimated using ARMAX on time domain data "z1". Fit to estimation data: 76.38% (prediction focus) FPE: 1.127, MSE: 1.082 More information in model's "Report" property.

Use `getpvec`

to fetch a flat list of all the model parameters, or just the free ones, and their uncertainties. Use `getcov`

to fetch the entire covariance matrix.

[par, dpar] = getpvec(m, 'free') CovFree = getcov(m,'value')

par = -1.5312 0.7293 0.9430 0.5224 -1.0587 0.1968 dpar = 0.0180 0.0147 0.0607 0.0782 0.0607 0.0596 CovFree = 0.0003 -0.0003 0.0000 0.0007 0.0004 -0.0003 -0.0003 0.0002 -0.0000 -0.0004 -0.0003 0.0002 0.0000 -0.0000 0.0037 -0.0034 -0.0000 0.0001 0.0007 -0.0004 -0.0034 0.0061 0.0008 -0.0005 0.0004 -0.0003 -0.0000 0.0008 0.0037 -0.0032 -0.0003 0.0002 0.0001 -0.0005 -0.0032 0.0035

nf = 0, nd = 0 denote orders of a general linear model, of which the ARMAX model is a special case.

Report contains information about the estimation process:

m.Report m.Report.DataUsed % record of data used for estimation m.Report.Fit % quantitative measures of model quality m.Report.Termination % search termination conditions

ans = Status: 'Estimated using ARMAX with prediction focus' Method: 'ARMAX' InitialCondition: 'zero' Fit: [1x1 struct] Parameters: [1x1 struct] OptionsUsed: [1x1 idoptions.polyest] RandState: [1x1 struct] DataUsed: [1x1 struct] Termination: [1x1 struct] ans = struct with fields: Name: 'z1' Type: 'Time domain data' Length: 300 Ts: 0.1000 InterSample: 'zoh' InputOffset: [] OutputOffset: [] ans = struct with fields: FitPercent: 76.3807 LossFcn: 1.0824 MSE: 1.0824 FPE: 1.1266 AIC: 887.1256 AICc: 887.4123 nAIC: 0.1192 BIC: 909.3483 ans = struct with fields: WhyStop: 'Near (local) minimum, (norm(g) < tol).' Iterations: 3 FirstOrderOptimality: 7.2436 FcnCount: 7 UpdateNorm: 0.0067 LastImprovement: 0.0067

To obtain on-line information about the minimization, use the 'Display' estimation option with possible values 'off', 'on', and 'full'. This launches a progress viewer that shows information on the model estimation progress.

Opt = armaxOptions('Display','on'); m1 = armax(z1,[2 2 2 1],Opt);

There are several types of linear models. The one above is an example of the `idpoly`

version for polynomial type models. Different variants of polynomial-type models, such as Box-Jenkins models, Output Error models, ARMAX models etc are obtained using the corresponding estimators - `bj, oe, armax, arx`

etc. All of these are presented as `idpoly`

objects.

Other variants are `idss`

for state-space models; `idgrey`

for user-defined structured state-space models; `idtf`

for transfer function models, and `idproc`

for process models (gain+delay+static gain).

The commands to evaluate the model: `bode, step, iopzmap, compare`

, etc, all operate directly on the model objects, for example:

compare(z1,m1)

Transformations to state-space, transfer function and zeros/poles are obtained by `idssdata`

, `tfdata`

and `zpkdata`

:

```
[num,den] = tfdata(m1,'v')
```

num = 0 0.9430 0.5224 den = 1.0000 -1.5312 0.7293

The 'v' means that num and den are returned as vectors and not as cell arrays. The cell arrays are useful to handle multivariable systems. To also retrieve the 1 standard deviation uncertainties in the values of `num`

and `den`

use:

```
[num, den, ~, dnum, dden] = tfdata(m1,'v')
```

num = 0 0.9430 0.5224 den = 1.0000 -1.5312 0.7293 dnum = 0 0.0607 0.0782 dden = 0 0.0180 0.0147

The objects also connect directly to the Control System Toolbox™ model objects, like `tf`

, `ss`

, and `zpk`

, and can be converted into these LTI objects, if Control System Toolbox is available. For example, `tf`

converts an `idpoly`

object into a `tf`

object.

CSTBInstalled = exist('tf','class')==8; if CSTBInstalled % check if Control System Toolbox is installed tfm = tf(m1) % convert IDPOLY model m1 into a TF object end

tfm = From input "u1" to output "y1": 0.943 z^-1 + 0.5224 z^-2 ---------------------------- 1 - 1.531 z^-1 + 0.7293 z^-2 Sample time: 0.1 seconds Discrete-time transfer function.

When converting an IDLTI model to the Control Systems Toolbox's LTI models, noise component is not retained. To also include the noise channels as regular inputs of the LTI model, use the 'augmented' flag:

if CSTBInstalled tfm2 = tf(m1,'augmented') end

tfm2 = From input "u1" to output "y1": 0.943 z^-1 + 0.5224 z^-2 ---------------------------- 1 - 1.531 z^-1 + 0.7293 z^-2 From input "v@y1" to output "y1": 1.051 - 1.113 z^-1 + 0.2069 z^-2 -------------------------------- 1 - 1.531 z^-1 + 0.7293 z^-2 Input groups: Name Channels Measured 1 Noise 2 Sample time: 0.1 seconds Discrete-time transfer function.

The noise channel is named `v@y1`

in the model `tfm2`

.