Mathworks SYSTEM IDENTIFICATION TOOLBOX 7 Reference

System Identification Toolbox™ 7

Reference

Lennart Ljung

How to Contact The MathWorks

www.mathworks.
comp.soft-sys.matlab Newsgroup
www.mathworks.com/contact_TS.html T echnical Support

[e-mail] Product enhancement suggestions
bugs@mathwo
[e-mail] Documentation error reports
[e-mail] Order status, license renewals, passcodes
info@mathwo

com

rks.com

rks.com

Web

Bug reports

Sales, prici

ng, and general information

508-647-7000 (Phone)

508-647-7001 (Fax)

The MathWorks, Inc.
3 Apple Hill Drive
Natick, MA 01760-2098

For contact information about worldwide offices, see the MathWorks Web site.

System Identification Toolbox™ Reference

© COPYRIGHT 1988–20 10 by The MathWorks, Inc.

The software described in this document is furnished under a license agreement. The software may be used
or copied only under the terms of the license agreement. No part of this manual may be photocopied or
reproduced in any form without prior written consent from The MathW orks, Inc.

FEDERAL ACQUISITION: This provision applies to all acquisitions of the Program and Documentation
by, for, or through the federal government of the United States. By accepting delivery of the Program
or Documentation, the government hereby agrees that this software or documentation qualifies as
commercial computer software or commercial computer software documentation as such terms are used
or defined in FAR 12.212, DFARS Part 227.72, and DFARS 252.227-7014. Accordingly, the terms and
conditions of this Agreement and only those rights specified in this Agreement, shall pertain to and govern
theuse,modification,reproduction,release,performance,display,anddisclosureoftheProgramand
Documentation by the federal government (or other entity acquiring for or through the federal government)
and shall supersede any conflicting contractual terms or conditions. If this License fails to meet the
government’s needs or is inconsistent in any respect with federal procurement law, the government agrees
to return the Program and Docu mentation, unused, to The MathWorks, Inc.

Trademarks

MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See

www.mathworks.com/trademarks for a list of additional trademarks. Other product or brand

names may be trademarks or registered trademarks of their respective holders.

Patents

The MathWorks products are protected by one or more U.S. patents. Please see

www.mathworks.com/patents for more information.

Revision History

September 2007 Online only Revised for Version 7.1 (Release 2007b)
March 2008 Online only Revised for Version 7.2 (Release 2008a)
October 2008 Online only Revised for Version 7.2.1 (Release 2008b)
March 2009 Online only Revised for Version 7.3 (Release 2009a)
September 2009 Online only Revised for Version 7.3.1 (Release 2009b)
March 2010 Online only Revised for Version 7.4 (Release 2010a)

Function Reference

1

Data Import and Processing ........................ 1-3

Contents

Linear Model Identification

Nonlinear Black -Box Model Identification

ODE Param eter Estimation

Recursive Model Identification

Model Analysis

Simulation and Prediction

System Identification Tool GUI

.................................... 1-13

........................ 1-5

......................... 1-12

......................... 1-16

Functions – Alphabetical List

2

........... 1-8

..................... 1-12

..................... 1-17

Block Reference

3

Data Import and Processing ........................ 3-2

Linear Model Identification

........................ 3-3

v

Simulation and Prediction ......................... 3-4

Blocks — Alphabetical List

4

Index

vi Contents

Function Reference

1

Data Import and Processing (p. 1-3)

Linear Model Identification (p. 1-5) Estimate time response, frequency

Nonlinear Black-Box Model
Identification (p. 1-8)

ODE Parameter Estimation (p. 1-12) Estimate parameters of linear and

Recursive Model Identification
(p. 1-12)

Model Analysis (p. 1-13)

Represent, process, analyze, and
manipulate data

response, transfer function,
input-output polynomial, and
state-space models from time and
frequency domain data

Estimate nonlinear ARX and
Hammerstein-Wiener models

nonlinear ordinary differential
or difference equations (grey-box
models)

Recursively estimate input-output
linear models, such as AR,
ARX, ARMAX, Box-Jenkins, and
Output-Error models

Validate and analyze model s by
comparing model output, computing
parameter confidence intervals and
prediction errors, and getting advice
on estimated models

1 Function Reference

Simulation and Prediction (p. 1-16) Simulate and predict linear and

nonlinear model output, and
estimate initial states

System Identification Tool GUI
(p. 1-17)

Start System Identification
Toolbox™ GUI and customize
preferences

1-2

Data Import and Processing

Data Import and Processing

advice

covf

delayest

detrend

diff

fcat

feedback

fft

fselect

get

getexp

getTrend

iddata

idfilt

idfrd

idinput

idresamp

Analysis and recommendations for
data or estimated linear polynomial
and state-space models

Estimate covariance functions for
time-domain iddata object

Estimate time delay (dead time)
from data

Subtract offset or trend from data
signals

Difference sig nals in iddata objects

Concatenate frequency-domain
signals in data objects

Identify possible feedback data

Transform iddata object to
frequency domain data

Frequencies from frequency response
data

Query p roperties of data and model
objects

Specific experiments from
multiple-experiment data set

Data offset and trend information

Time- or frequency-domain d ata

Filter data using user-defined
passbands, g eneral filters, or
Butterworth filters

Frequency-response data or model

Generate input signals

Resample time-domain data by
decimation or interpolation

1-3

1 Function Reference

ifft

isreal

merge (iddata)

misdata

nkshift

pexcit

plot

realdata

resample

set

size

timestamp

TrendInfo

Transform iddata objects from
frequency to time domain

Determine whether model
parameters or data values are
real

Merge data sets into iddata object

Reconstruct missing input and
output data

Shift data sequences

Level of excitation of input signals

Plot iddata or model objects

Determine w hether iddata is based
on real-valued signals

Resample time-domain data by
decimation or interp olatio n (requires
Signal P rocess ing Toolbox™
software)

Set properties of data and model
objects

Dimensions of data and model
objects

Return date and time when object
was created or last modified

Offset and linear trend slope values
for detrending data

1-4

Linear Model Identification

Linear Model Identification

ar

armax

arx

arxdata

arxstruc

balred

bj

c2d

cra

d2c

delayest

etfe

feedback

frd

freqresp

Estimate parameters of AR model
for scalar time series

Estimate parameters of ARMAX or
ARMA model

Estimate parameters of ARX or AR
model using least squares

ARX parameters from
multiple-output models with
variance information

Compute and compare loss functions
for single-output ARX models

Reduce model order (requires
Control System Toolbox™ product)

Box-Jenkins (BJ) model estimation

Transform linear model from
continuous to discrete time

Estimate impulse response using
prewhitened-based correlation
analysis

Transform linear model from
discrete to continuous time

Estimate time delay (dead time)
from data

Estimate empirical transfer
functions and periodograms

Identify possible feedback data

Convert idfrd objects to Control
System Toolbox frequency-response
LTI model

Frequency response data from linear
models

1-5

1 Function Reference

get

idarx

idfrd

idgrey

idmodel

idpoly

idproc

idss

impulse

init

iv4

ivar

ivstruc

ivx

LTI Commands

merge

Query p roperties of data and model
objects

Multiple-output ARX polynomials,
impulse response, or step response
model

Frequency-response data or model

Linear ODE (grey-box model) with
known and unknown parameters

Superclass for linear models

Linear polynomial input-output
model

Linear, low-order, continuous-time
transfer function

State-space model

Plot impulse response with
confidence interval

Set or randomize initial parameter
values

Estimate ARX model using
four-stage instrumental variable
method

Estimate AR model using
instrumental variable method

Loss functions for sets of ARX model
structures

Estimate parameters of ARX model
using instrumental variable method
with arbitrary instruments

Apply Control System Toolbox
commands to linear model

Merge estimated models

1-6

Linear Model Identification

n4sid

nuderst

oe

pem

pexcit

polydata

selstruc

set

setpname

setPolyFormat

setstruc

size

spa

spafdr

Estimate state-space models using
subspace method

Set step size for numerical
differentiation

Output-error (OE) model parameter
estimation

Estimate model parameters
using iterative prediction-error
minimization method

Level of excitation of input signals

Parameters from single-input and
single-output polynomial model

Select model order for single-output
ARX models

Set properties of data and model
objects

Set mnemonic parameter names for
linear black-box model structures

Specify format for B and F
polynomials of m ulti-input
polynomial model for backward
compatibility

Set matrix structure for idss model
objects

Dimensions of data and model
objects

Estimate frequency response with
fixed frequency resolution using
spectral analysis

Estimate frequency response and
spectrum using spectral analysis
with frequency-dependent resolution

1-7

1 Function Reference

ss

ssdata

step

struc

tf

tfdata

timestamp

zpk

zpkdata

Convert linear models to Control
System Toolbox LTI models

State-space matrices from
parametric linear model

Plot step response with confidence
interval

Generate model-order combinations
for single-output ARX model
estimation

Convert linear models to
transfer-function Control System
Toolbox LTI models

Numerator and denominator of
transfer function from linear model

Return date and time when object
was created or last modified

Convert linear model to Control
System Toolbox state-space LTI
models

Zeros, poles, and gains of transfer
function from linear model

Nonlinear Black-Box Model Identification

addreg

customnet

customreg

1-8

Add custom regressors to nonlinear
ARX model

Custom nonlinearity estimator
for nonlinear ARX and
Hammerstein-Wiener models

Custom regre ssor for nonlinear ARX
models

Nonlinear Black-Box Model Identification

data2state(idnlarx)

deadzone

evaluate

findop(idnlarx)

findop(idnlhw)

get

getDelayInfo

getreg

idnlarx

idnlhw

idnlmodel

init

linapp

linear

linearize(idnlarx)

linearize(idnlhw)

Map past input/output data to
current states of nonlinear ARX
model

Class representing dead-zone
nonlinearity estimator for
Hammerstein-Wiener models

Value of nonlinearity estimator at
given input

Compute operating point for
nonlinear ARX model

Compute operating point for
Hammerstein-Wiener model

Query p roperties of data and model
objects

Get input/output delay information
for

idnlarx model structure

Regressor expressions and numerical
values in nonlinear ARX model

Nonlinear ARX model

Hammerstein-Wiener model

Superclass for nonlinear m odels

Set or randomize initial parameter
values

Linear approximation of nonlinear
ARX and Hammerstein -Wiener
models for given input

Specify to estimate no n linear ARX
model that is linear in (nonlinear)
custom regressors

Linearize nonlinear ARX model

Linearize Hammerstein-Wiener
model

1-9

1 Function Reference

neuralnet

nlarx

nlhw

operspec(idnlarx)

operspec(idnlhw)

pem

poly1d

polyreg

pwlinear

saturation

set

sigmoidnet

Class representing neural
network object created in Neural
Network Toolbox™ product for
estimating nonlinear ARX and
Hammerstein-Wiener models

Estimate nonlinear ARX model

Estimate Hammerstein-Wiener
model

Construct operating point
specification object for

idnlarx

model

Construct operating point
specification object for

idnlhw

model

Estimate model parameters
using iterative prediction-error
minimization method

Class representing single-variable
polynomial nonlinear estimator for
Hammerstein-Wiener models

Powers and products of standard
regressors

Class representing piecewise-linear
nonlinear estimator for
Hammerstein-Wiener models

Class representing saturation
nonlinearity estimator for
Hammerstein-Wiener models

Set properties of data and model
objects

Class representing sigmoid network
nonlinearity estimator for nonlinear
ARX and Hammerstein -Wiener
models

1-10

Nonlinear Black-Box Model Identification

treepartition

unitgain

wavenet

Class representing binary-tree
nonlinearity estimator for nonlinear
ARX models

Specify absence of nonlinearities for
specific input or output channels in
Hammerstein-Wiener models

Class representing wavelet network
nonlinearity estimator for nonlinear
ARX and Hammerstein -Wiener
models

1-11

1 Function Reference

ODE Parameter Estimation

get

getinit

getpar

idgrey

idnlgrey

idnlmodel

init

pem

set

setinit

setpar

Query p roperties of data and model
objects

Values of idnlgrey model initial
states

Parameter values and properties of

idnlgrey model parameters

Linear ODE (grey-box model) with
known and unknown parameters

Nonlinear ODE (grey-box model)
with unknown parameters

Superclass for nonlinear m odels

Set or randomize initial parameter
values

Estimate model parameters
using iterative prediction-error
minimization method

Set properties of data and model
objects

Set initial states of idnlgrey model
object

Set initial parameter values of

idnlgrey model object

Recursive Model Identification

armax

r

rarx

1-12

timate recursively parameters of

Es

RMAX or ARMA models

A

stimate parameters of ARX or AR

E
models recursively

Model Analysis

rbj

roe

rpem

rplr

segment

Model Analysis

advice

aic

arxdata

balred

bode

compare

Estimate recursively parameters of
Box-Jenkins models

Estimate recursively output-error
models (IIR-filters)

Estimate general input-output
models using recursive
prediction-error minimization
method

Estimate general input-output
models using recursive pseudolinear
regression method

Segment data and estimate models
for each segment

Analysis and recommendations for
data or estimated linear polynomial
and state-space models

Akaike Information Criterion for
estimated model

ARX parameters from
multiple-output models with
variance information

Reduce model order (requires
Control System Toolbox product)

Compute and plot frequency
response magnitude and phase for
logarithmic frequencies

Compare model output and
measured output

1-13

1 Function Reference

ffplot

fpe

freqresp

fselect

impulse

isreal

ivstruc

noisecnv

nyquist

pe

plot

polydata

predict

predict(idnlarx)

predict(idnlgrey)

predict(idnlhw)

Compute and plot frequency
response magnitude and phase for
linear frequencies

Akaike Final Prediction Error for
estimated model

Frequency response data from linear
models

Frequencies from frequency response
data

Plot impulse response with
confidence interval

Determine whether model
parameters or data values are
real

Loss functions for sets of ARX model
structures

Transform idmodel object with noise
channels to model with measured
channels only

Plot Nyquist curve of frequency
response with confidence interval

Prediction errors associated with
model and data set

Plot iddata or model objects

Parameters from single-input and
single-output polynomial model

Predict output k steps ahead

Predict output k steps ahead for
nonlinear ARX model

Predict output k steps ahead for
nonlinear O DE model

Predict output k steps ahead for
Hammerstein-Wiener model

1-14

Model Analysis

present

pzmap

resid

selstruc

sim

sim(idnlarx)

sim(idnlgrey)

sim(idnlhw)

simsd

ssdata

step

tfdata

view

zpkdata

Display model information, including
estimated uncertainty

Plotzerosandpoleswithconfidence
interval

Compute and test model residuals
(prediction errors)

Select model order for single-output
ARX models

Simulate linear models with
confidence interval

Simulate nonlinear ARX model

Simulate nonlinear ODE model

Simulate Hammerstein-Wiener
model

Simulate models with uncertainty
using Monte Carlo method

State-space matrices from
parametric linear model

Plot step response with confidence
interval

Numerator and denominator of
transfer function from linear model

Plot model characteristics using
Control Sy stem Toolbox LTI Viewer
GUI

Zeros, poles, and gains of transfer
function from linear model

1-15

1 Function Reference

Simulation and Prediction

findstates(idmodel)

findstates(idnlarx)

findstates(idnlgrey)

findstates(idnlhw)

predict

predict(idnlarx)

predict(idnlgrey)

predict(idnlhw)

retrend

sim

sim(idnlarx)

sim(idnlgrey)

sim(idnlhw)

simsd

Estimate initial states of linear
model from data

Estimate initial states of nonlinear
ARX model from data

Estimate initial states of nonlinear
grey-box model from data

Estimate initial states of nonlinear
Hammerstein-Wiener model from
data

Predict output k steps ahead

Predict output k steps ahead for
nonlinear ARX model

Predict output k steps ahead for
nonlinear O DE model

Predict output k steps ahead for
Hammerstein-Wiener model

Add offsets or trends to data signals

Simulate linear models with
confidence interval

Simulate nonlinear ARX model

Simulate nonlinear ODE model

Simulate Hammerstein-Wiener
model

Simulate models with uncertainty
using Monte Carlo method

1-16

System Identification Tool GUI

System Identification Tool GUI

ident

midprefs

Open System Identification Tool
GUI

Set folder for storing idprefs.mat
containing GUI startup information

1-17

1 Function Reference

1-18

Functions – Alphabetical
List

2

addreg

Purpose Add custom regressors to nonlinear ARX model

Syntax m = addreg(model,regressors)

m = addreg(model,regressors,output)

Description m = addreg(model,regressors) adds custom regressors to a nonlinear

ARX model by appending the
and m are idnalrx objects. For single-output models, regresso rs is
an object array of regressors you create using
or a cell array of string expressions. For multiple-output models,

regressors is 1-by-ny cell array of customreg objects or 1-by-ny cell

array of cell arrays of string expressions.
the

ny cells to the corresponding model output channel. If regressors

is a single regressor, addreg adds this regressor to all output channels.

m = addreg(model,regressors,output) adds regressors regressors

to specific output channels output of a multiple-output model. output
is a scalar integer or vector of integers, where each integer is the index
of a model output channel. Specify several pairs of

output values to add different regresso r variables to the corresponding

output channels.

CustomRegressors model property. model

customreg or polyreg,

addreg adds each element of

regressors and

Examples Add regressors to a nonlinear ARX model as a cell array of strings:

% Create nonlinear ARX model with standard regressors:

m1 = idnlarx([4 2 1],'wavenet','nlr',[1:3]);

% Create model with additional custom regressors:

m2 = addreg(m1,{'y1(t-2)^2';'u1(t)*y1(t-7)'})

% List all standard and custom regressors of m2:

getreg(m2)

Add regressors to a nonlinear ARX model as customreg objects:

% Create nonlinear ARX model with standard regressors:

m1 = idnlarx([4 2 1],'wavenet','nlr',[1:3]);

% Create a model based on m1 with custom regressors:

2-2

r1 = customreg(@(x)x^2, {'y1'}, 2)
r2 = customreg(@(x,y)x*y, {'u1','y1'}, [0 7])
m2 = addreg(m1,[r1 r2]);

See Also customreg | getreg | nlarx | polyreg

How To • “Identifying Nonlinear A RX Models”

addreg

2-3

advice

Purpose Analysis and recommendations for data or estimated linear polyn omial

and state-space models

Syntax advice(model)

advice(data)

Inputs model

Name of the idarx, idgrey, idpoly, idproc,oridss model
object. These model objects belong to the
representing linear polynomial and state-space models.

data

Name of the iddata object.

idmodel abstract class,

Description advice(model) displays the following inform a tion about the estimated

model in the MATLAB

• Does the model capture essential dynamics of the system and the

disturbance characteristics?

• Is the model order higher than necessary?

• Is there potential output feedback in the validation data?

• Would a nonlinear ARX model perform better than a linear ARX

model?

advice(data) displays the following information about the data in

the MATLAB Command Window:

• What are the excitation levels of the signals and how does this affects

the model orders? See also

• Doesitmakesensetoremoveconstantoffsetsandlineartrendsfrom

thedata? Seealso

• Isthereanindicationofoutputfeedbackinthedata? Seealso

feedback.

®

Command Window:

pexcit.

detrend.

2-4

advice

See Also

detrend

feedback

iddata

pexcit

2-5

aic

Purpose Akaike Information Criterion for estimated model

Syntax am = aic(model)

am = aic(model1,model2,...)

Arguments model

Name of an idarx, idgrey, idpoly, idproc, idss, idnlarx,

idnlhw,oridnlgrey model object.

Description am = aic(model) returns a scalar value of the Akaike’s Information

Criterion (AIC) for the estimated

am = aic(model1,model2,...) returns a row vector containing AIC

values for the estimated models

Remarks Akaike’s Information Criterion (AIC) provides a measure of model

quality by simulating the situation where the model is tested on a
different data set. After computing several different models, you can
compare them using this criterion. According to Akaike’s theory, the
most accurate model has the smallest AIC.

model.

model1,model2,....

2-6

Note If you use the same data set for both model estimation and
validation, the fit always improves as you increase the model order and,
therefore, the flexibility of the model structure.

Akaike’s Information Criterion (AIC) is defined by the following
equation:

d

AIC V

=+log

where V is the loss function, d is the num ber of e stimated parameters,
and N is the number of values in the estimation data set.

The loss function V is defined by the following equation:

2

N

N

⎝

⎠

⎛

Vtt

=

1

det , ,

εθ εθ

⎜
⎜

()()

∑

N

1

()

NN

⎞

T

⎟
⎟

aic

where

For d<<N:

Note AIC is approximately equal to log(FPE).

AIC is formally defined as the negative log-likelihood functionΛ,
evaluated at the estimated parameters, plus the number of estimated
parameters. You can derive AIC from this definition, as follows:

If the disturbance source is Gaussian with the covariance matrix
logarithm of the likelihood function is

Maximizing this analytically with respect toΛ, and then maximizing
the result with respect to

θ

represents the estimated parameters.

N

d

N

∑

2

N

T

⎞
⎟

⎠

⎞
⎟

⎠

−

θ

1

,gives

N

()

2

Λ

+

⎛
⎜

⎝

1

AIC V

Ltt const

( , ) ( , ) ( , ) log detθεθεθΛΛ Λ=− −

⎛

=+

log 1

⎜
⎝

1
2

,the

Np

L const V

(, ) log( )θΛ = + +

where p is the number of outputs.

N

22

References Ljung, L. System Identification: Theory for the User, Upper Saddle

River, NJ, Prentice-Hal PTR, 1 999. Seesectionsaboutthestatistical
framework for parameter estimation and maximum likelihood method
and comparing model structures.

2-7

aic

See Also

EstimationInfo

fpe

2-8

Algorithm Properties

Purpose Algorithm properties affecting estimation process for linear models

Syntax idprops idmodel algorithm

Model.algorithm.PropertyName='PropertyValue '

Description Algorithm is a property of the idmodel class that specifies the

estimation algorithm. The
various linear models, including

idgrey. These models inherit the idmodel Algorithm property.

Note For a description of nonlinear model Algorithm property, see the
corresponding nonlinear model reference page.

Property names are not case sensitive. When you type a property
name, you only need to enter enough characters to uniquely identify
the property. The Algorithm fields can be accessed and modified as
for any structure using dot syntax. For example, you can access the

SearchMethod field by typing Model.Algorithm.SearchMethod.

idmodel subclasses are used to define

idarx, idss, idpoly, idproc,and

Note You can use the get function or dot notation to fetch fields of

Algorithm as if they were the properties of the model itself. Similarly

you can use
access is available for the fields of
struct fields such as Search or Threshold options). For example:

method = Model.SearchMethod;
Model.MaxIter = 100;

is equivalent to

get(Model, 'SearchMethod')
set(Model, 'maxiter', 100);

set or dot notation to set a particular field. This shortcut

Algorithm only (not for deeper level

2-9

Algorithm Properties

When you create a new model by refining an existing model m,the
algorithm p roperties of

Note You can estimate a model with specific algorithm settings and
define a structure variable to store the algorithm values. For example:

model = n4sid(data,order)
myalg = model.Algorithm;
myalg.Focus='Simulation';
m = pem(data,model,'alg',myalg)

You can also specify the algorithm properties (except advanced
properties) as property-value pairs when creating the linear model
(using

armax, oe etc.).

Algorithm Properties

idpoly, idss, etc.) or when estimating them (using pem, n4sid,

m are inherited by the new model.

2-10

• Criterion: Specifies criterion used during minimization. Criterion

can have the fo llowing values:

- 'Det': Minimize

error. This is the optimal choice in a statistical sense and leads to
the maximum likelihood estimates in case nothing is known about
the variance of the noise. It uses the inverse of the estimated
noise variance as the weighting function. This is the default
criterion used for all models, except
by default.

det( ’* )EE

,whereE represents the prediction

idnlgrey which uses 'Trace'

- 'Trace': Minimize the trace of the weighted prediction error

matrix
with one column for each output, and
symmetric matrix of size equal to the number of outputs. By
default,
outputs (so the minimization criterion becomes
the traditional least-sum-of-squared-errors criterion. You can

trace(E'*E*W),whereE is the matrix of prediction errors,

W is a positive semi-definite

W is an identity matrix of size equal to the number of model

trace(E'*E),or

Algorithm Properties

specify the relative weighting of prediction errors for each output
using the

Note The difference between the two criteria is meaningful in
multiple-output cases only. In single-output models, the two criteria
are equivalent. Both the
general requirement of minimizing a weighted sum of squares of
prediction errors. The
the covariance matrix of the noise source and using the inverse of
that matrix as the weighting. When using the
must specify the weighting using the

If you want to achieve better accuracy for a particular channel
in multiple-input multiple-output models, you should use

Trace with weighting that favors that channel. Otherwise it is

natural to use the

cond(model.NoiseVariance) after estimation. If the matrix is

ill-conditioned, it may be more robust to use the
You can also use
relative error for different channels corres ponds to your needs or
expectations. Use the
relative errors, and check
weighting modifications to specify.

Weighting field of the Algorithm property.

Det and Trace criteria are derived f rom a

Det criterion can be interpreted as estimating

Trace criterion, you

Weighting property.

Det criterion. When using Det, you can check

Trace criterion.

compare on validation data to check whether the

Trace criterion if you need to modify the

model.NoiseVariance to determine what

The search method of

lsqnonlin supports the Trace criterion only.

• Focus: Defines how the errors e between the measured and the

modeled outputs are weighed at specific frequencies during the
minimization of the following loss function:

Ve

=∑λ

iii

2

2-11

Algorithm Properties

Higher weighting at specific frequencies emphasizes the requirement
for a good fit at these frequencies.
values:

- 'Prediction': (Default) Automatically calculates the weighting

function as a product of the input spectrum and the inverse of the
noise model. This minimizes the one-step-ahead prediction, which
typically favors f itting small time intervals (higher frequency
range). From a statistical-variance point of view, this is the
optimal weighting function. However, this method neglects the
approximation aspects (bias) of the fit. Might not result in a stable
model. Use

- 'Simulation': Estimates the model using the frequency

weighting of the transfer function that is given by the input
spectrum. Typically, this method favors the frequency range
where the input spectrum has the most power. In other words, the
resulting model will produce good simulations for inputs that have
the same spectra as used for estimation. For models that have no
disturbance model, there is no difference between

and 'Prediction'.Inthiscase,
is equivalent to
models.

Focus can have th e following

'Stability' when you want to ensure a stable model.

'Simulation'

yGuHe=+

A=C=D=1 for idpoly models and K=0for idss

with H=1,which

2-12

For models that have a disturbance model, G is first estimated
with H=

with a fixed estimated transfer function
stable transfer function G.

1,andthenH is estimated using a prediction-error method

ˆ

G

. This guarantees a

- 'Stability': This weighting is the same as for 'Prediction',

but the model is forced to be stable. Use only when you know the
system is stable. In some cases, forcing the model to be stable
can result in a bad model.

- Enter a row vector or a matrix containing frequency values that

define desired passbands. For example:

[wl,wh]
[w1l,w1h;w2l,w2h;w3l,w3h;...]

Algorithm Properties

where wl and wh represent upper and lower limits of a passband.
For a matrix with several rows defining frequency passbands, the
algorithm uses union of frequency ranges to define the estimation
passband.

- Enter any SISO linear filter in any of the following ways:

A single-input-single-output (SISO)

ss, tf,orzpk model from the Control System Toolbox product.

An

Using the format
matrices of the filter.

Using the format
numerator and denominator of the filter transfer function

This calculates the weighting function as a product of the filter and
the input spectrum to estimate the transfer function from input
to output, G. To obtain a good model fit for a specific frequency
range, you must choose the filter with a passband in this range.
After estimating G, the algorithm computes the disturbance
model using a prediction-error method and keeping the estimated
transfer function fixed (similar to the
model that has no disturbance model, the estimation result is the
same if you first prefilter the data using

{A,B,C,D}, which specifies the state-space

{numerator, denominator}, which specifie s the

idmodel object.

'Simulation' case). For a

idfilt.

- For frequency-domain data only, enter a column vector of weights

for

'Focus'. This vector must have the same size as length of the

frequency vector of the data set,
output response in the data is multiplied by the corresponding
weight at that frequency.

•

Maxsize: A positive integer, specified such that the input-output

data is split into segments where each contains fewer than
elements. Setting MaxSize can improve computational performance.
The default value of
set the value. You can edit this file to optimize computational speed
on a particular computer.
properties of the estimate except when you use

'backcast'

for frequency-domain data. In this case, the frequency

MaxSize is 'Auto', which uses the file idmsize to

MaxSize does not affect the numerical

Data.Frequency. Each input and

MaxSize

InitialState =

2-13

Algorithm Properties

ranges where backcasting takes place might depend on MaxSize and
affects estimates.

•

FixedParameter: Vector of integers containing the indices of

parameters that are not estimated and remain fixed a t nominal or
initial values. Parameter indices refer to their position in the list,
stored in the property
parameter names as values from the property
fixed parameters using parameter names, enter
a cell array of strings. For example, to fix parameters with names

'a' and 'b',typem.FixedParameter = {'a','b','c'}.Stringscan

contain wildcards, such as
astring,or
fix three parameters in a disturbance model that start with
as

'k1', 'k2','k3',useFixedParameter = {'k*'}. For structured

state-space models, you can fix and free parameters by specifying
structure matrices, such as

Note By default, the property 'PName' is e mpty. Use setpname to
assign default parameter names. For example,

'ParameterVector'. You can also specify

'PName'.Tospecify

Fixedparameter as

'*' to specify t he a utomatic completion of

'?' to represent an arbitrary character. For example, to

'k',such

As and Bs (see idss).

m = setpname(m).

2-14

• Weighting: Positive semi-definite symmetric matrix W to use as

weighting for minimization of the trace criterion

Weighting can be used to specify relative importance of outputs in

multiple-input multiple-output models (or reliability of corresponding
data) by specifying

Weighting is not useful in single-input single-output models. By

default,
of model outputs, assigning equal importance to each output during
estimation.

•

Display: S pecifies what information displays in the MATLAB

Command Window about the iterative search during estimation.

Weighting is the identity matrix of size equal to the number

W as a diagonal matrix of non-negative values.

trace(E'*E*W).

- 'Off': D isplays no information.

- 'On': Displays the loss-function values for each iteration.

Algorithm Properties

- 'Full': Displays the same information as 'On' and also i nclude

the current parameter values and the search direction (except
when the
to

'Free' for idss models and the list of parameters can change

between iterations).

•

LimitError: Specifies when to adjust the weight of large errors

from quadratic to linear. Default value is

LimitError times the estimated standard deviation have a linear

weight in the criteria. The standard deviation is estimated robustly
as the median of the absolute deviations from the median and
divided by

LimitError = 0 disables the robustification and leads to a purely

quadratic criterion. When estimating with frequency-domain data,

LimitError is set to zero.

Note You can estimate the model with the default value
of

LimitError (zero) and plot the prediction errors using

pe(data.model). If the resulting plot shows occasional large values,

repeat the estimat ion with
value between

Advanced SSParameterization model property is set

0. Errors larger than

0.7. (See the section about choosing a robust norm in [2].)

model.Algorithm.LimitError set to a

1 and 2.

• MaxIter: Specifies the maximum number of iterations during

loss-function minimization. The iterations stop when

MaxIter

is reached or another stopping criterion is satisfied, such as
the

Tolerance. The default value of MaxIter is 20. Setting
MaxIter = 0 returns the result of the startup procedure. Use
EstimationInfo.Iterations to get the actual number of iterations

during an estimation.

•

Tolerance: Specifies the minimum percentage difference (divided

by 100) between the current value of the loss function and its
expected improvem ent after the next iteration: When the percentage
of expected improvement is less than
are stopped. Default value is

0.01. The estimate of the expected

Tolerance, the iterations

2-15

Algorithm Properties

loss-function improvement at the next iteration is made based on the
Gauss-Newton vector computed for the current parameter value.

•

SearchMethod: The search method used for iterative parameter

estimation. It can take one of the following values:

- 'gn': The subspace Gauss-Newton direction.

Singular values of the Jacobian matrix less than

GnPinvConst*eps*max(size(J))*norm(J) are discarded when

computing the search direction, where
The Hessian matrix is approximated by
improvement along this direction, the gradient direction is also
tried.

- 'gna': An adaptive version of subspace Gauss-Newton approach,

suggested by Wills and Ninness (IFAC World congress, Prague

2005). Eigenvalues less than
neglected , where
The Gauss-Newton direction is compute d in the rem aining
subspace.
is increased by the factor
find a lower value of t h e criterion in less than 5 bisections. It is
decreased by the factor
without any bisections.

J is the Jacobian matrix.

JTJ. If there is no

gamma*max(sv) of the Hessian are

sv are the singular values of the Hessian.

gamma has the initial value InitGnaTol (see below) and

LmStep each time the search fails to

2*LmStep each time a search is successful

2-16

- 'lm': Uses the Levenberg-Marquardt method. This means that

the next parameter value is
one, where
the gradient.
of the criterion is found.

H is the Hessian, I is the identity m atrix, and grad is

d is a number that is increased until a lower value

-pinv(H+d*I)*grad fro m the previous

- 'Auto': A choice among the above is made in the algorithm. This

is the default choice .

- 'lsqnonlin':Useslsqnonlin optim izer from Optimization

Toolbox™ software. You must have O ptimizatio n Toolbox installed
to use this option. This search method can only handle the
criterion.

Trace

Algorithm Properties

• Advanced: Structure that specifies advanced algorithm options and

has the following fields:

- Search: Uses the following fields to specify options for the iterative

search:

1

GnPinvConst: Must be a positive real value. Specifies

that singular values of the Jacobian that are smaller than

GnPinvConst*max(size(J)*norm(J)*eps are discarded when

computing the search direction u sing the
value is

InitGnaTol: The initial value of gamma in the gna search

2

algorithm. See

10^-4.

LmStep: The size of the Levenberg-Marquardt step. The next value

3

of the search-direction length
method is

StepReduction: For search directions other than the

4

Levenberg-Marquardt direction, the step is reduced by the factor

StepReduction after each iteration. Default is 2.

1e4.

SearchMethod for description of gna.Defaultis

d in the Levenberg-Marquardt

LmStep times the previous one. Default is 2.

'gn' method. Default

MaxBisection: The maximum number of bisections used b y the

5

line search along the search direction. Default is

LmStartValue: The starting value of search-direction length d in

6

the Levenberg-Marquardt method. Default is

RelImprovement: The iterations are stopped if the relative

7

improvement of the criterion is less than
is

RelImprovement = 0. This property is different from Tolerance

in that it uses the actual improvement of the loss function, as
opposed to the expected improvement.

RelImprovement.Default

25.

0.001.

- Threshold: Contains fields wi th thresholds for several tests:

Sstability: Specifies the location of the rightmost pole to test the

1

stability of continuous-time models. A model is considered stable
when its rightmost pole is to the left of

Sstability. Default is 0.

2-17

Algorithm Properties

2 Zstability: Specifies the maximum distance of all poles from

the origin to test stability of discrete-time models. A model is
considered stable if all poles are within the distance
from the origin. Default is 1+sqrt(eps).

- AutoInitialState: Specifies when to automatically estimate the

initial state. When
estimated when the ratio of the prediction-error norm with a zero
initial state to the norm with an estimated initial state exceeds

AutoInitialState. Default is 1.05.

Properties Relevant to Estimation of n4sid, State-Space (idss)
Models

Note These properties apply to n4sid.Sincepem com m only uses n4sid

to initialize a model for iterative estimation, these properties affect
the results of

Zstability

InitialState = 'Auto', the initial state is

pem too.

2-18

• N4Weight: Calculates the w eighting matrices used in the

singular-value decomposition step of the algorithm and has three
possible values:

- 'Auto': (Default) Automatically chooses between 'MOESP' and

'CVA' .

- 'MOESP': Uses the MOESP algorithm by Verhaegen.

- 'CVA': Uses the canonical variable algorithm by Larimore.

For more information about setting this property, see the
reference page.

N4Horizon: Determines the forward and backward p rediction

•

horizons used by the algorithm. Enter a row vector with three
elements:
prediction horizon; that is, the algorithm uses up to
predictors.
of past inputs used for predictions. For an exact definition of these
integers, see the section about subspace methods in [2], where they

N4Horizon=[r sy su],wherer is the maximum forward

sy is the number of past outputs, and su is the number

n4sid

r step-ahead

Algorithm Properties

are called r, s1,ands2. These numbers can have a substantial
influence on the quality of the resulting model and there are no
simple rules for choosing them. Making
means that the algorithm tries each row of
the value that gives the best (prediction) fit to the data. Choosing
the best row is not available when you also specify to select the best
model order. When you specify one co lumn in
interpretation is

'Auto'

, which uses an Akaike Information Criterion (AIC) for the

selection of

r=sy=su. The default choice is 'N4Horizon' =

sy and su.

Note For algorithm properties of nonlinear models, see the reference
pages for

idnlarx, idnlhw,andidnlgrey.

References [1] Dennis, J.E., Jr., and R.B. Schnabel, Numerical Methods for

Unconstrained Optimization and Nonlinear Equations,Prentice

Hall, Englewood Cliffs, N.J., 1983. See the chapter about iterative
minimization.

'N4Horizon' a k-by-3 matrix

'N4Horizon' and selects

'N4Horizon',the

See Also

[2] Ljung, L. System Identification: Theory for the User, Upper Saddle
River, NJ, Prentice-Hal PTR, 1999. See the chapter about computing
the estimate.

armax

bj

EstimationInfo

idpoly

idss

n4sid

oe

pem

2-19

ar

Purpose Estimate parameters of AR model for scalar time series

Syntax m = ar(y,n)

[m,refl] = ar(y,n,approach,window)
[m,refl] = ar(y,n,approach,window,'P1',V1,...,'PN',VN)

Arguments y

iddata

channel).

n

Scalar that specifies the order of the model you want to estimate
(the number of A parameters in the AR model).

approach

Lets you choose the algorithm for computing the least squares AR
model from the following options:

•

object that contains the time-series data (one output

'burg': Burg’s lattice-based method. Solves the lattice filter

equations using the harmonic mean of forward and backward
squared prediction e rrors.

2-20

•

'fb': (Default) Forward-backward approach. Minimizes the

sum of a least- squares criterion for a forward model, and the
analogous criterion for a time-reversed model.

•

'gl': Geometric lattice approach. Similar to Burg’s method,

but uses the geometric mean instead of the harmonic me an
during minimization.

•

'ls': Least-squares approach. Minimizes the standard sum of

squared forward-prediction errors.

•

'yw': Yule-Walker approach. Solves the Yule-Walker

equations, formed from sample covariances.

window

Lets you specify how to use information about the data outside the
measured time interval (past and future values). The following
windowing options are available:

• 'now': (Default) No windowing. This value is the default

except when the
data is used to form regression vectors. The summation in the
criteria starts at the sample index equal to

'pow': Postwindowing. M issing end values are replaced with

•

zeros and the summation is extended to time
number of observations).

•

'ppw': Pre- and postwindowing. Us e d in the Yule-Walker

approach.

•

'prw': Prewindowing. Missing past values are replaced with

zeros so that the summ ation in the criteria can start at time
equal to zero.

'P1',V1,...,'PN',VN

Pairs of property names and property values can include any of
the following.

approach argument is 'yw'.Onlymeasured

n+1.

N+n (N is the

ar

Property Name

'Covariance' • 'None'

Property Value

Suppresses the
calculation of
the covariance
matrix.

• []

Empty.

• Square matrix

containing
covariance
values of size
equal to the
length of the
parameter vector

Description

Specifies
calculation of
uncertainties
in parameter
estimates.

2-21

ar

Description

Property Name

'MaxSize'

'Ts'

Note Use for scalar time series only. For multivariate data, use arx.

m = ar(y,n) returns an idpoly model m.

[m,refl] = ar(y,n,approach,window) returns an idpoly model m

and the variable refl. For the two lattice-based approaches, 'burg'
and 'gl', refl stores the reflection coefficients in the first row, and the
corresponding loss function values in the second row. The first column
of

refl is the zeroth-order model, and the (2,1) element of refl is

the norm of the time series itself.

[m,refl] = ar(y,n,approach,window,'P1',V1,...,'PN',VN)

returns an idpoly model m and the variable refl using additional
windowing criteria.

Property Value

Integer

Real positive
number

Description

Algorithm

See

Properties

description.

Sets the sampling
time and overrides
the sampling time
of

y.

for the

Remarks The AR model structure is given by the following equation:

Aqyt et( ) () ()=

AR model parameters are estimated using variants of the least-squares
method. The following table summarizes the common names for
methods with a specific combination of
values.

2-22

approach and window argument

Method Approach and Windowing

Modified Covariance Method (Default) Forward-backward

approach and no windowing.

ar

Correlation Method

Covariance Method

Yule-Walker approach, which
corresponds to least squares plus
pre- and postwindowing.

Least squares approach with no
windowing. arx uses this routine.

Examples Given a sinusoidal signal with noise, compare the spectral estimates of

Burg’s method with those found from the forward-backward approach
and no-windowing method on a Bode plot.

y = sin([1:300]') + 0.5*randn(300,1);
y = iddata(y);
mb = ar(y,4,'burg');
mfb = ar(y,4);
bode(mb,mfb)

References Marple, Jr., S.L., Digital Spectral Analysis with Applications,Prentice

Hall, Englewood Cliffs, 1987, Chapter 8.

See Also

Algorithm Properties

arx

EstimationInfo

etfe

idpoly

ivar

pem

2-23

ar

spa

step

2-24

Purpose Estimate parameters of ARMAX or ARMA model

Syntax m = armax(data,orders)

m = armax(data,orders,'P1',V1,...,'PN',VN)
m = armax(data,'na',na,'nb',nb,'nc',nc,'nk',nk)

Arguments data

iddata

orders

Vector of integers, specified using the format

For m ultiple-input systems, nb and nk are row vectors where the

ith element corresponds to the order and delay associated with

the

When
then

object that contains the input-output data.

orders = [na nb nc nk]

ith input.

data is a time series, which has no input and one output,

armax

orders = [na nc]

Tip When refining an estimated model mi, set the model orders
as follows:

orders = mi

'na',na,'nb',nb,'nc',nc,'nk',nk

, 'nb',and'nc' are orders of the ARMAX model. nk is the

'na'

delay.

'P1',V1,...,'PN',VN

Pairs of property names and property values can include any of
the following

na, nb, nc,andnk are the corresponding integer values.

idmodel properties:

2-25

armax

Description

'Focus', 'InitialState', 'Display', 'MaxIter', 'Tolerance',
'LimitError',and'FixedParameter'.

See

Algorithm Properties, idpoly,andidmodel for more

information.

Note armax only supports time-domain data w ith single or multiple
inputs and single output. For frequency-domain data, use
multiple-output case, use ARX or a state-space model (see

pem).

m = armax(data,orders) returns an idpoly model m with estimated

oe.Forthe

n4sid and

parameters and covariances (parameter uncertainties). Estimates the
parameters using the prediction-error method and specified

m = armax(data,orders,'P1',V1,...,'PN',VN) returns an idpoly

orders.

model m. Use additional property-value pairs to specify the estimation
algorithm properties.

m = armax(data,'na',na,'nb',nb,'nc',nc,'nk',nk) returns an
idpoly model m with orders and delays specified as parameter-value

pairs.

Remarks The ARMAX model structure is

yt a yt a yt n

() ( ) ( )

2-26

+−++ −=

butn butn n

1

Amorecompactwaytowritethedifferenceequationis

Aqyt Bqut n Cqet

()() ()( ) ()()=−+

where

11…

1

() (

−++ −−+

… ))

knkb

1 …

() ( )()

−++ − + cet c et n et

1

na

a

b

nc

c

k

+

yt()

•

•

•

•

•

—Outputattimet.

—Numberofpoles.

n

a

n

—Numberofzeroesplus1.

b

n

—NumberofC coefficients.

c

— Number of input samples that occur before the input affects

n

k

the output, also called the dead time in the system. For discrete
systems with no dead time, there is a minimum 1–sample delay

because the output depends on the previous input and

yt yt n

()( )−−1 …

•

— Previous outputs on which the current output

a

nk= 1

depends.

armax

.

•

ut n ut n n

()( )−−−+… 1

kkb

— Previous and delayed inputs on which

the current output depends.

•

et et n

()( )−−1 …

The parameters

— White-noise disturbance value.

c

na, nb,andnc are the orders of the ARMAX model, and

nk is the delay. q is the delay operator. Specifically,

n

−

Aq aq a q

()=+ + +

1

Bq b bq b q

()=+ ++

12

Cq cq c q

()=+ + +

1

1

…

1

−

1

…

−

1

…

1

−

a

n

a

n

−+

1

b

n

b

n

−

c

n

c

If data is a time series, which ha s no input channels and one output
channel, then

armax calculates an ARMA model for the time series

Aqyt et( ) () ()=

2-27

armax

In this case

orders = [na nc]

Algorithm An iterative search algorithm with the properties 'SearchMethod',

'MaxIter', 'Tolerance',and'Advanced' minimizes a robustified

quadratic prediction error criterion. The iterations are terminated
either when
is less than
be found. You can get information about the search criteria using

m.EstimationInfo.

When you do not specify initial parameter values for the iterative search
in

orders, they are constructed in a special four-stage LS-IV algorithm.

The cutoff value for the robustification is based on the property

LimitError and on the estimated standard deviation of the residuals

from the initial parameter estimate. It is not recalculated during the
minimization.

To ensure that only models corresponding to stable predictors are tested,
the algorithm performs a stability test of the predictor. Generally, both

and

Cq()

MaxIter is reached, or when the expected improvement

Tolerance, or when a lower value of the criterion cannot

(if applicable) must have all zeros inside the unit circle.

Fq()

Minimization information is displayed on the screen when the property

'Display' is 'On' or 'Full'.With'Display' ='Full',boththe

current and the previous parameter estimates are displayed in
column-vector form, listing parameters in alphabetical order. Also,
the values of the criterion function are given and the Gauss-Newton
vector and its norm are also displayed. With
the criterion values are displayed.

'Display' = 'On' only

References Ljung, L. System Identification: Theory for the User, Upper Saddle

River, NJ, Prentice-Hal PTR, 1999. See chapter about computing the
estimate.

2-28

armax

See Also

Algorithm Properties

EstimationInfo

idpoly

pem

2-29

arx

Purpose Estimate parame te rs of A RX or AR model using least squares

Syntax m = arx(data,orders)

m = arx(data,orders,'P1',V1,...,'PN',VN)
m = arx(data,'na',na,'nb',nb,'nc',nc,'nk',nk)

Arguments data

An iddata object, an frd object, or an idfrd
frequency-response-data object.

orders

Vector of integers, specified using the format

orders = [na nb nk]

For m ultiple-input systems, nb and nk are row vectors where the

ith element corresponds to the order and delay associated with

the

ith input.

2-30

When
then

Note When refining an estimated model mi, set the model orde rs
as follows:

'na',na,'nb',nb,'nc',nc,'nk',nk

'na'

delay.

data is a time series, which has no input and one output,

orders = [na]

orders = mi

, 'nb',and'nc' are orders of the ARMAX model. nk is the

na, nb, nc,andnk are the corresponding integer values.

Description

'P1',V1,...,'PN',VN

Pairs of property names and property values can include any of
the following

'Focus', 'InitialState', 'Display', 'MaxIter', 'Tolerance',
'LimitError',and'FixedParameter'.

See

Algorithm Properties, idpoly,andidmodel for more

idmodel properties:

information.

Note arx does not support multiple-output continuous-time models.
Use state-space model structure instead. When the true noise term

et et n

()( )−−1 …

in the ARX model structure is not white noise and na

c

is nonzero, the model estimate is incorrect. In this case, use armax,

bj, iv4,oroe.

m = arx(data,orders) returns a model m as an idpoly or idarx

object with estimated parameters and covariances (parameter
uncertainties). For single-output data, the model is an
For multiple-output data, the model is an
least-squares method and specified

orders.

idarx object. Uses the

idpoly object.

arx

m = arx(data,orders,'P1',V1,...,'PN',VN) returns a model m.Use

additional property-value pairs to specify the estimation algorithm
properties.

m = arx(data,'na',na,'nb',nb,'nc',nc,'nk',nk) returns a model
m with orders and delays specified as parameter-value pairs.

Remarks arx estimates the parameters of the ARX model structure:

yt a yt a yt na

() ( ) ... ( )

+ −++ − =

b u t nk b u t nb nk

( ) ... ( )

−++ −−+

1

1

1

nb

na

1 ++ et()

2-31

arx

The parameters na and nb are the orders of the ARX model, and nk
is the delay.

yt()

•

•

•

•

—Outputattimet.

n

—Numberofpoles.

a

—Numberofzeroesplus1.

n

b

n

— Number of input samples that occur before the input affects

k

the output, also called the dead time in the system. For discrete
systems with no dead time, there is a minimum 1–sample delay

because the output depends on the previous input and

nk= 1

.

•

yt yt n

()( )−−1 …

— Previous outputs on which the current output

a

depends.

ut n ut n n

()( )−−−+… 1

•

kkb

— Previous and delayed inputs on which

the current output depends.

et et n

•

()( )−−1 …

— White-noise disturbance value.

c

Amorecompactwaytowritethedifferenceequationis

Aqyt Bqut n et

()() ()( ) ()=−+

k

q is the delay operator. Specifically,

n

−

Aq aq a q

()=+ + +

1

Bq b bq b q

()=+ ++

12

1

…

1

−

1

…

n

a

n

b

−

a

n

−+

b

1

Time Series Models

For time-series data that contains no inputs, one output and orders =

, the model has AR structure of order na.

na

2-32

The AR model structure is

Aqyt et( ) () ()=

Multiple Inputs and Single-Output Models

For multiple-input systems, nb and nk are row v ectors where the ith
element corresponds to the order and delay associated with the

ith

input.

yt A yt A yt A yt na

() ( ) ( ) ( )

+−+−++ −=

But But

−−++ − +1) ( ) ( )… Butnb et

12…

12

() (

+

01

na

nb

Multi-Output Models

For m odels with multiple inputs and multiple outputs, na, nb,andnk
contain one row for each output signal.

arx

In the multiple-output case,

arx minimizes the trace of the prediction

error covariance matrix, or the norm

N

T

etet

() ()

t

=∑1

To transform this to an arbitrary quadratic norm using a weighting
matrix

Lambda

N

T

et et

() ()Λ

∑

t

=

1

−

1

use the syntax

m = arx(data,orders,'NoiseVariance', Lambda)

You can use arx to refine an existing model m_initial as an argum

m = arx(data,m_initial)

ent.

2-33

arx

The new model m uses the orders and the weighting matrix for the
prediction errors from

m_initial.Youcanfurthermodifym_initial

by adding a list of property name and value p airs as arguments.
This is especially useful when some parameters must be fixed using

'FixedParameter' property.

Continuous-Time Models

For models with one output, continuous-time models can be estimated
from continuous-time frequency-domain data. In this case,
number of estimated denominator coefficients and

nb is number of

na is the

estimated numerator coefficients.

Note For continuous-time models, omit the delay parameter nk because
it has no meaning in this case. Because estimating continuous-time
ARX models often produces bias, you might get better results by using
the

oe method.

2-34

For example, when na = 4, nb = 2, the model structure is:

bs b

+

Gs

()=

4

sasasasa

++++

3

1

12

2

2

34

Tip When using continuous-time data, limit the fit to a smaller
frequency range using the

m = arx(datac,[na nb],'focus',[0 wh])

'Focus' idmodel property:

Estimating Initial Conditions

For time-domain data, the signals are shifted such that unmeasured
signals are never required in the predictors. Therefore, there is no need
to estimate initial conditions.

For frequency-domain data, it might be necessary to adjust the data by
initial conditions that support circular convolution.

arx

You can set the property

•

'zero' — No adjustment.

'estimate' — Perform adjustment to the data by initial conditions

•

that support circular convolution.

•

'auto' — Auto matically choose between 'zero' and 'estimate'

based on the data.

See

Algorithm Properties for more information on model properties.

'InitialState' to one of the following values:

Algorithm QR factorization solves the overdetermined set of linear equations that

constitutes the least-squares estimation problem.

The regression matrix is formed so that only measured quantities are
used (no fill-out with zeros). When the regression matrix is larger
than

MaxSize, data is segmented and QR factorization is performed

iteratively on these data segments.

Examples This example generates input data based on a specified ARX model, and

then uses this data to estimate an ARX model.

A = [1 -1.5 0.7]; B = [0 1 0.5];
m0 = idpoly(A,B);
u = iddata([],idinput(300,'rbs'));
e = iddata([],randn(300,1));
y = sim(m0, [u e]);
z = [y,u];
m = arx(z,[2 2 1]);

See Also Algorithm Properties | EstimationInfo | ar | idarx | idpoly

| iv4 | ivar | ivx | pem

How To • “Using Linear Model for Nonlinear ARX Estimation”

2-35

arxdata

Purpose ARX parameters from multiple-output models with variance

information

Syntax [A,B] = arxdata(m)

[A,B,dA,dB] = arxdata(m)

Arguments m

An idarx model object.

Also accepts single-output
ARX structure with orders

idpoly models with an underlying

nc=nd=nf=0.

Description [A,B] = arxdata(m) returns A and B as 3-D arrays.

Suppose ny is the number of outputs (the dimension of the vector y(t))
and nu is the number of inputs.

A is an ny-by-ny-by-(na+1) array such that

A(:,:,k+1) = Ak
A(:,:,1) = eye(ny)

where k=0,1,...,na.

B is an ny-by-nu-by-(nb+1) array with

B(:,:,k+1) = Bk

A(0)

is always the identity matrix. The leading entries in B equal to

zero, which means there are no delays in the model.

Note For a time series, B=[].

2-36

[A,B,dA,dB] = arxdata(m) returns A and B matrices, and dA and dB

as the estimated standard deviations of A and B,respectively.

arxdata

Remarks A and B are 2-D or 3-D arrays and are returned in the standard

multivariable ARX format (see

yt A yt A yt A yt na

() ( ) ( ) ( )

+−+−++ −=

But But

−−++ − +1) ( ) ( )… Butnb et

where Akand Bkmatrices have dimensions ny-by-ny and ny-by-nu,
respectively. ny is the number of outputs (the dimension of the vector
y(t))andnu is the number of inputs.

12…

12

() (

+

01

idarx), describing the model.

na

nb

See Also

idarx

idpoly

2-37

arxstruc

Purpose Compute and compare loss functions for single-output ARX models

Syntax V = arxstruc(ze,zv,NN)

V = arxstruc(ze,zv,NN,maxsize)

Arguments ze

Estimation data set can be iddata or idfrd object.

zv

Validation data set can be iddata or idfrd object.

NN

Matrix defines the number of different ARX-model structures.
Each row of

nn = [na nb nk]

maxsize

Specified maximum data size. See Algorithm Properties for
an explanation.

NN is o f the form:

Description

2-38

Note Use arxstruc for single-output systems only. arxstruc supports
both single-input and multiple-input systems.

V = arxstruc(ze,zv,NN) returns V, which contains the loss functions

in its first row. The remaining rows of
that the orders and delays are given just below the corresponding loss
functions. The last column of

V = arxstruc(ze,zv,NN,maxsize) uses the additiona l specification

of the maximum data size.

with the same interpretation as described for
generation of typical

NN matrices.

V contains the n umber of da t a points in ze.

V contain the transpose of NN,so

arx.Seestruc for easy

arxstruc

The output argument V is best analyzed using selstruc.Theselection
of a suitable m odel structure based on the information in
done using

selstruc.

Remarks Each of ze and zv is an iddata object containing output-input data.

Frequency-domain data and
for each of the model structures defined by
data set

ze. The loss functions (normalized sum of squared prediction

errors) are then computed for these models when applied to the
validation data set

zv. The data sets ze and zv need not be of equal size.

They could, however, be the same sets, in which case the computation
is faster.

idfrd objects are also supported. Models

NN are estimated using the

Examples This example uses the simulation data from a second-order idpoly

model w ith additive noise. The data is split into two parts, where one
part is the estimation data and the other is the validation data. You
select the best model by comparing the output of models with orders
ranging between 1 and 5 with the validating data. All models have an
input-to-output delay of 1.

v is normally

% Create an ARX model for generaing data:
A = [1 -1.5 0.7]; B = [0 1 0.5];
m0 = idpoly(A,B);
% Generate a random input signal:
u = iddata([],idinput(400,'rbs'));
e = iddata([],0.1*randn(400,1));
% Simulate the output signal from the model m0:
y = sim(m0, [u e]);
z = [y,u]; % analysis data
NN = struc(1:5,1:5,1);
V = arxstruc(z(1:200),z(201:400),NN);
nn = selstruc(V,0);
m = arx(z,nn);

2-39

arxstruc

See Also

Algorithm Properties

arx

idpoly

ivstruc

selstruc

struc

2-40

balred

Purpose Reduce model order (requires Control System Toolbox product)

Syntax MRED = balred(M)

MRED = balred(M,ORDER,'DisturbanceModel','None')

Description This function reduces the order of any model M given as an idmodel

object. The resulting reduced-order model, MRED,isanidss model.

The function requires routines from the Control System Toolbox
product.

ORDER: The desired order (dimension of the state-space representation).

If

ORDER = [], which is the default, a plot shows how the diagonal

elements of the observability and controllability Gramians of a balanced
realization decay with the order of the representation. You are then
prompted to select an orde r based on this plot. The idea is that such a
small element has a negligible influence on the input-output behavior
of the model. We recommend that you choose an order s uch that only
large elements in these matrices are retained.

'DisturbanceModel': If the property DisturbanceModel is set to
'None',thenanoutput-errormodelMRED is produced: that is, one with

the Kalman gain K equal to zero. Otherwise (default), the disturbance
model is also reduced.

The function recognizes whether
model and performs the reduction accordingly. The resulting model,

MRED,issimilartoM in this respect.

There are several options for how the reduction is performed:

RelTol, Offset, Elimination

M is a continuous- or discrete-time

AbsTol,

.

Algorithm The function uses the balred algorithm in Control System Toolbox. The

plot, in case

ORDER = [], shows the vector g returned by balreal.

Examples Build a high-order multivariable ARX model, reduce its order to 3,and

compare the frequency responses of the original and reduced models:

M = arx(data,[4*ones(3,3),4*ones(3,2),ones(3,2)]);

2-41

balred

MRED = balred(M,3);
bode(M,MRED)

Use the reduced-order model as an initial condition for a third-order
state-space model.

M2 = pem(data,MRED);

See Also

balreal

2-42

Purpose Box-Jenkins (BJ) model estimation

Syntax m = bj(data,[nb nc nd nf nk])

m = bj(data,[nb nc nd nf nk],'PropertyName',PropertyValue)
m = bj(data,m_initial)

Description m = bj(data,[nb nc nd nf nk]) estimates Box-Jenkins model

parameters and their covariances from input-output data.

idpoly object. data is a time-domain, single-output iddata object.
nb, nc, nd, and nf are orders of the B, C, D,andF polynomials,

respectively.

nk is the input delay, specified as the number of samples.

Orders and delay are scalar for single-input data, and row vectors for
multiple-input data with the same size as the number of input channels.

m = bj(data,[nb nc nd nf nk],'PropertyName',PropertyValue)

estimates Box-Jenkins model using algorithm options specified by

idpoly property name-value pairs. See Algorithm Properties.

m = bj(data,m_initial) refines previouslyestimatedmodel
m_initial,whichisanidpoly object.

m is an

bj

bj does not support frequency-domain and multiple-output data.

Definitions The general Box-Jenkins m odel structure is:

nu

Bq

()

i

=∑1

i

ut nk

ii

Fq

()

i

()

yt

()

where nu is the number of input channels.

The orders of Box-Jenkins model are defined as follows:

nbBqbbq bq

:

() ...

=+ ++

12

nc C q c q c

:

:

nd D q d q d q

:

nf F q f q

1

() ...

=+ + +
=+ + +

() ...

1

=+ +

() ...

1

−−+

−

1

1

−−

1

1

−

1

1

Cq

()

et

Dq

()

nb

q

n

cc

nd

nf

()=−

nb

nc

−

nd

nf

+

11

++−fq

2-43

bj

Examples Estimate parameters of a sing le-input single-output Box-Jenkins model:

% Load SISO data.

load iddata1;

% Estimate model parameters

mbj = bj(z1,[2 2 2 2 1])

Estimate parameters of a multi-input single-output Box-Jenkins model:

% Load MISO data.

load iddata8;

% Estimate model parameters

mbj = bj(z8,[[2 1 1] 1 1 [2 1 2] [5 10 15]])

Estimate parameters of a single-input single-output Box-Jenkins model
using estimation algorithm properties:

2-44

% Generate estimation data using simulation.

B = [0 1 0.5];
C=[1-10.2];
D = [1 1.5 0.7];
F = [1 -1.5 0.7];
m0 = idpoly(1,B,C,D,F,0.1);
e = iddata([],randn(200,1));
u = iddata([],idinput(200));
y = sim(m0,[u e]);
z=[yu];

% Estimate model parameters.

mbj_i = bj(z,[2 2 2 2 1]);

% Repeat the estimation with more iterations.

mbj = bj(z,mbj_i,'MaxIter',50)

% View the estimation results.

mbj.EstimationInfo

% Compare initial and refined model parameters.

compare(z,mbj,mbj_i)

References Ljung, L. System Identification: Theory for the User,2nded.,Upper

Saddle River, NJ, Prentice-Hall, 1999. See the chapter on computing
the estimate.

See Also idmodel | oe | idpoly | n4sid | pem | Algorithm Properties |

EstimationInfo

Tutorials • “Tutorial – Identifying Linear Models Using the Command Line”

How To • “Identifying Input-Output Polynomial Models”

• “Algorithms for Estimating Polynomial Models”

bj

2-45

bode

Purpose Compute and plot frequency response magnitude and phase for

logarithmic frequencies

Syntax bode(m)

bode(m,w)
bode(m('noise')
bode(m1,...,mN,'sd',sd,'mode','same','ap',ap,'fill')
[mag,phase,w] = bode(m)
[mag,phase,w,sdmag,sdphase] = bode(m)

Description bode(m) plots a Bode plot for the model m, which can be an idpoly,

idss, idarx, idgrey,oridfrd object. This frequency response is a

function of logarithmic frequencies in radians per unit time (stored as
the

TimeUnit model property). D efault frequency values are computed

from the model dynamics. For time series spectra, phase plots are
omitted. For MIMO models, press Enter to view the next plot in the
sequence of different I/O channel pairs, annotated using the
and OuputNames model properties.

bode(m,w) plots a Bode plot at specified frequencies w in radians per

unit time, which can be

InputNames

2-46

• Avectorofvalues.

•

{wmin,wmax}, which specifie s 100 logarithmically spaced frequency

values ranging from a minimum value

wmax.

{wmin,wmax,np}, which specifies np logarithmically spaced frequency

•

values.

Note For idfrd models, you cannot specify individual frequencies
and can only limit the frequencies range for the internally stored
frequencies using

{wmin,wmax}.

wmin and a maximum value

bode(m('noise') plots a Bode plot of the output noise spectra when

the model contains noise spectrum information.

bode(m1,...,mN,'sd',sd,'mode','same','ap',ap,'fill') plots a

Bode plot for several models.

sd specifies the confidence region as a

positive number that represents the number of standard deviations.
The argument

mode = 'same' displays all I/O channels in the sam e plot. Set ap =

to show only amplitude plots, or ap = 'P' to show only phase plots.

'A'

[mag,phase,w] = bode(m) computes the magnitude mag and
phase values of the frequency response, w hich are 3-D arrays with

'fill' indicates that the confidence region is color filled.

dimensions (number of outputs)-by-(number of inputs)-by-(length of
w).

w specifies the frequency values for computing the response even

if you did not specify it as an input. For SISO systems,

mag(1,1,k)

and phase(1,1,k) are the magnitude and phase (in degrees) at the
frequency w(k). For MIMO s ystems,
the frequency response at frequency
similarly for
spectrum and

phase(i,j,k).Whenm is a time series, mag is its power

phase is zero.

mag(i,j,k) is the magnitude of

w(k) from input j to output i,and

bode

See Also

[mag,phase,w,sdmag,sdphase] = bode(m) computes the standard

deviations of the magnitude
array of the same size as
size as

phase.

etfe

ffplot

freqresp

idfrd

nyquist

spa

spafdr

sdmag and the phase sdphase. sdmag is an

mag,andsdphase is an array of the same

2-47

compare

Purpose Compare model output and measured output

Syntax compare(data,m);

compare(data,m,k)
compare(data,m,k,'Samples',sampnr,'InitialState',init,'OutputPlots
',Yplots)
compare(data,m1,m2,...,mN)
compare(data,m1,'PlotStyle1',...,mN,'PlotStyleN')
[yh,fit,x0] = compare(data,m1,'PlotStyle1',...,mN,'PlotStyleN',k)

Description data is the output-input data in the usual iddata object format. data

can a lso b e an idfrd object with frequenc y-re sponse data.

compare computes the output yh that results when the model m is

simulated with the input
object. The result is plotted together with the corresponding measured
output

y. The percentage of the output variation that is explained by

the model

fit = 100*(1 - norm(yh - y)/norm(y-mean(y)))

u. m can be any idmodel or idnlmodel model

2-48

is also computed and displayed. For multiple-output systems, this is
done separately for each output. For frequency-domain data (or in
general for complex valued data) the
onlytheabsolutevaluesof

When the argument

y and yh are plotted.

k is specified, the k step-ahead prediction of y

fit is still calculated as above, but

according to the model m are computed instead of the simulated output.
In the calculation of yh(t), the model can use outputs up to time t–k:
y(s),s = t –k,t –k –1, … (and inputs up to the current time t). The
default value of
only. Note that for frequency-domain data, only simulation (
is allowed, and for time-series data (no input) only prediction (

inf) is possible.

k is inf, which gives a pure simulation from the input

k = inf)

k not

compare

Property Name/Property Value Pairs

The optional property name/property value pairs 'Samples'/sampnr,

'InitialState'/init,and'OutputPlots'/Yplots can be given in

any order.

The argument
with

OutputName in this array are plotted, while all are used for the

necessary computations. If

Yplots can be a cell array of strings. Only the outputs

Yplots is not specified, all outputs are

plotted.

The argument

sampnr indicates that only the sample numbers in this

row vector are plotted and used for the calculation of the fit. The whole
data record is used for the simulation/prediction.

The argument

init determines how to handle initial conditions in the

models:

•

init = 'e' (for 'estimate') estimates the initial conditions for

best fit.

•

init = 'm' (for 'model') uses internally stored initial state of the

model.

•

init = 'z' (for 'zero') uses zero initial conditions.

init = x0,wherex0 is a column vector of the same size as the state

•

vector of the models, uses

init = 'e' is the default.

•

x0 as the initial state.

Several Models

When several models are specified, as in compare(data,m1,m2,...,mN),
the plots sho w responses and fits for all models. In that case

data

should contain all inputs and outputs that are required for the different
models. However, some models might correspond to subselections of
channels and might not need all channels in
proper handling of signals is based on the

data. In that case the

InputNames and OutputNames

of data and the models.

2-49

compare

With compare(data,m1,'PlotStyle1',...mN,'PlotStyle2'),the
color, line style, and/or marker can be specified for the curves associated
with the different m odels . The markers are the same as for the regular

plot command. For example,

compare(data,m1,'g_*',m2,'r:')

If data contains several experiments, separate plots are given for the
different experiments. In this case
array w ith as many entries as there are experiments.

Arguments When output arguments [yh,fit,x0] = compare(data,m1,..,mN)

are specified, no plots are produced.

yh is a cell array of length equal to the number of models. Each cell

contains the corresponding model output as an

fit is, in the general case, a 3-D array with fit(kexp,kmod,ky)

containing the fit (computed as above) for output ky, model kmod,and
experiment

kexp.

sampnr,ifspecified,mustbeacell

iddata object.

x0 is a cell array, such that x0{kmod} is the estimated initial state for

model number
whose column number
number

kmod. If data is multiexperiment, X0{kmod} isamatrix

kexp is the initial state vector for experiment

kexp.

Examples Split the data record into two parts. Use the first one for estimating a

model and the second one to check the model’s ability to predict six
steps ahead.

ze = z(1:250);
zv = z(251:500);
m= armax(ze,[2 3 1 0]);
compare(zv,m,6);
compare(zv,m,6,'Init','z') % No estimation of

% the initial state.

2-50

compare

See Also

findstates(idmodel)

pe

predict

sim

2-51

covf

Purpose Estimate covariance functions for time-domain iddata object

Syntax R = covf(data,M)

R = covf(data,M,maxsize)

Description data is an iddata object and M is the maximum delay -1 for which

the covariance function is estimated. The routine is intended for
time-domain data only.

Let

z contain the output and input channels

yt

()

⎡

zt

()

where y and u are the rows of data.OutputData and data.InputData,
respectively, with a total of

R is returned as an nz

Ri j nzk

(( ), ) ()( ) ()+− += +=

⎤

=

⎢

⎥

ut

()

⎣

⎦

11

nz channels.

2

-by- M matrix with entries

N

1

ztz t k R k

ij

∑

N

t

=

1



ij

2-52

where zjis the jth row of z, and missing values in the sum are replaced
by zero.

The optional argument
under

Algorithm Properties.

maxsize controls the memory size as explained

The easiest way to describe and unpack the result is to use

reshape(R(:,k+1),nz,nz) = E z(t)*z'(t+k)

Here ' is complex conjugate transpose, which also explains how complex
data is handled. The expectation symbol

E corresponds to the sample

means.

covf

Algorithm When nz is at most two, and w hen permitted by maxsize,afastFourier

transform technique is applied. Otherwise, straightforward summing
is used.

See Also

iddata

spa

2-53

cra

Purpose Estimate impulse response using prewhitened-based correlation

analysis

Syntax cra(data);

[ir,R,cl] = cra(data,M,na,plot);
cra(R);

Description data is the output-input data given as an iddata object. The routine is

intended for time-domain data only.

The routine only handles single-input-single-output data pairs. (For the
multivariate case, apply

cra prewhitens the input sequence; that is, cra filters u through a

filter chosen so that the result is as uncorrelated (white) as possible.
The output
functions of the filtered
cross correlation function between (prewhitened) input and output is
also computed and graphed. Positive values of the lag variable then
correspond to an influence from
significant correlation for negativelagsisanindicationoffeedback
from

y to u in the data.

y is subjected to the same filter, and then the covariance

cra to two signals at a time, or use impulse.)

y and u are computed and graphed. The

u to later values of y.Inotherwords,

2-54

A properly scaled version of this correlation function is also an estimate
of the system’s impulse response
99% confidence levels. The output argument
estimate, so that its first entry corresponds to lag zero. (Ne gative lags
are excluded in
it corresponds to an impulse of height 1
the sampling interval of the data.

The output argument
as follows:

• The first column of

• The second column contains the covariance function of the (possibly

filtered) output.

ir.) In the plot, the impulse response is scaled so that

R contains the covariance/correlation information

R contains the lag indices.

ir. This is also graphed along with

ir is this impulse response

/T and duration T,whereT is

• The third column contains the covariance function of the (possibly

prewhitened) input.

• The fourth column contains the correlation function. The plots can

be redisplayed by

cra(R).

cra

The output argument

cl is the 99% confidence level for the impulse

response estimate.

The optional argument
covariance/correlation functions are computed. These are from
so that the length of

0 to M. The default value of M is 20.

For the prewhitening, the input is fitted to an A R model of order
The third argument of

na = 10.Withna = 0 the covariance and correlation functions of the

M defines the number of lags for which the

-M to M,

R is 2M+1. The impulse response is computed from

na.

cra can change this order from its default value

original data sequences are obtained.

plot: plot = 0 gives no plots. plot = 1 (the default) gives a plot of

the estimated impulse response together with a 99% confidence region.

plot = 2 gives a plot of all the covariance functions.

An often better alternative to

cra is the functions impulse and step,

which use a high-order FIR model to estimate the impulse response.

Examples Compare a second-order ARX model’s impulse response with the one

obtained by correlation analysis.

ir = cra(z);
m = arx(z,[2 2 1]);
imp = [1;zeros(19,1)];
irth = sim(m,imp);
subplot(211)
plot([ir irth])
title('impulse responses')
subplot(212)
plot([cumsum(ir),cumsum(irth)])
title('step responses')

2-55

cra

See Also

impulse

step

2-56

customnet

Purpose Custom nonlinearity estimator for nonlinear ARX and

Hammerstein-Wiener models

Syntax C=customnet(H)

C=customnet(H,PropertyName,PropertyValu e)

Description customnet is an object that stores a custom nonlinear estimator with a

user-defined unit function. This custom unit function uses a weighted
sum of inputs to compute a scalar output.

Construction C=customnet(H) creates a nonlinearity estimator object with a

user-defined unit function using the function handle
to a function of the form
value of the function,
range. Name the function
the interval
custom nonlinearity have only one input and one output.

C=customnet(H,PropertyName,PropertyValu e) creates a nonlinearity

estimator using property-value pairs defined in “customnet Propertie s”
on page 2-58.

[-a a]. Hammerstein-Wiener models require that your

[f,g,a] = gaussunit(x),wheref is the

g=df/dx,anda indicates the unit function active

gaussunit.m. g is significantly nonzero in

H. H must point

Remarks Use customnet to define a nonlinear function

scalar and x is an
on the following function expansion with a possible linear term L:

Fx xrPLaf xrQb c() ( )=− + −

aaf x rQb c d

where f is a unit function that you define using the function handle H.

P and Q are

projection matrices P and Q are determined by principal component
analysis of estimation data. Usually,
the estimation data are linearly dependent, then
columns of Q,
in the unit function.

m-dimensional row vector. The unit function is based

+

+

()

()

111

+

m-by-p and m-by-q projection matrices, respectively. The

q, corresponds to the number of components of x used

−

()

()

nnn

…

+

+

p=m. If the components of x in

yFx= ()

,wherey is

p<m. The number of

2-57

customnet

When used to estimate nonlinear ARX models, q is equal to the size of
the

NonlinearRegressors property of the idnlarx object. When used

to estimate Hammerstein-Wiener models,

m=q=1 and Q is a scalar.

customnet
Properties

r is a 1-by-

m vector and represents the mean value of the regressor

vector computed from estimation data.

d, a,andc are scalars.

p-by-1 vector.

L is a

b represents

q-by-1 vectors.

The function handle of the unit function of the custom net must have the
form

[f,g,a] = function_name(x). This function must be vectorized,

which means that for a vector or matrix

g musthavethesamesizeasx and be computed element-by-element.

x, the output arguments f and

You can include property-value pairs in the constructor to specify the
object.

After creating the object, you can use

get or dot notation to access the

object property values. For example:

% List all property values
get(C)
% Get value of NumberOfUnits property
C.NumberOfUnits

You can also use the set function to set the value of particular
properties. For example:

2-58

set(C, 'LinearTerm', 'on')

The first argument to set must be the name of a MATLAB variable.

customnet

Property Name

NumberOfUnits

LinearTerm

Parameters

Description

Integer specifies the number of nonlinearity units in the
expansion.
Default=

10.

For example:

customnet(H,'NumberOfUnits',5)

Can have the following values:

•

'on'—Estimates the vector L in the expansion.

'off'—Fixes the vecto r L to zero.

•

For example:

customnet(H,'LinearTerm','on')

A structure containing the parameters in the nonlinear
expansion, as follows:

•

RegressorMean:1-by-m vector containing the means of x

in estimation data, r.

NonLinearSubspace: m-by-q matrix containing Q.

•

LinearSubspace: m-by-p matrix containing P.

•

UnitFcn

LinearCoef: p-by-1 vector L.

•

Dilation: q-by-1 matrix containing the values b_k.

•

Translation:1-by-n vector containing the values c_k.

•

OutputCoef: n-by-1 vector containing the values a_k.

•

OutputOffset:scalard.

•

Typically, the values of this structure are set by estimating a
model with a

customnet nonlinearity.

Stores the function handle that points to the unit function.

2-59

customnet

Algorithm customnet uses an iterative search technique for estimating

parameters.

Examples Define custom unit function and save it in gaussunit.m:

function [f, g, a] = GAUSSUNIT(x)
% x: unit function variable
% f: unit function value
%g:df/dx
% a: unit active range (g(x) is significantly
% nonzero in the interval [-a a])

% The unit function must be "vectorized": for
% a vector or matrix x, the output arguments f and g
% must have the same size as x,
% computed element-by-element.

% GAUSSUNIT customnet unit function example
[f, g, a] = gaussunit(x)
f = exp(-x.*x);
if nargout>1

g = - 2*x.*f;
a = 0.2;

end

Use custom networks in nlarx and nlhw model estimation commands:

% Define handle to example unit function.
H = @gaussunit;
% Estimate nonlinear ARX model using
% Gauss unit function with 5 units.
m = nlarx(Data,Orders,customnet(H,'NumberOfUnits',5));

See Also evaluate | nlarx | nlhw

How To • “Identifying Nonlinear A RX Models”

2-60

• “Identifying Hammerstein-Wiener Models”

customnet

2-61

customreg

Purpose Custom regressor for nonlinear ARX models

Syntax C=customreg(Function,Variables)

C=customreg(Function,Variables,Dela ys,Vectorized)

Description customreg class represents arbitrary functions of past inputs and

outputs, such as products, powers, and other MATLAB expressions of
input and output variables.

You can specify custom regressors in addition to or instead of standard
regressors for greater flexibility in modeling your data using nonlinear
ARX models. Fo r example, you can define regressors like tan(u(t-1)),

2

u(t-1)

,andu(t-1)*y(t-3).

For simpler regressor expressions, specify custom regressors directly
in the GUI or in the
expressions, create a
specify these objects as inputs to the estimation. Regardless of how you
specify custom regressors, the toolbox represents these regressors as

customreg objects. Use getreg to list the expressions of all standard

and custom regressors in your model.

nlarx estimation command. For more complex

customreg object for each custom regressor and

2-62

A special case of custom regressors involves polynomial combinations
of past inputs and outputs. For example, it is common to capture
nonlinearities in the system using polynomial expressions like y(t−1)2,
u(t−1)2, y(t−2)2, y(t−1)*y(t−2), y(t−1)*u(t−1), y(t− 2)*u(t−1). At the
command line, use the
regressors automatically by computing all combinations of input and
output variables up to a specified degree.
objects that you specify as inputs to the estimation.

The nonlinear ARX model (
as the
using
MIMO models, to retrieve the

m.CustomRegressors{ky}(r).

Use the
regressors using vectorized form during estimation. If you know

CustomRegressors property. You can list all custom regressors

m.CustomRegressors,wherem is a nonlinear ARX model. For

Vectorized property to specify whether to compute custom

polyreg command to generate polynomial-type

polyreg produces customreg

idnlarx object) stores all custom regressors

rth custom regressor for output ky,use

customreg

that your regressor formulas can be vecto r ize d, set Vectorized to 1
to achieve better performance. To better understand vectorization,
consider the custom regressor function handle
and y are vectors and each variable is evaluated over a time grid.
Therefore,
are concatenated to produce a

for k = 1:length(x)

end

z must be evaluated for each (xi,yi) pair, and the results

z vector:

z(k) = x(k)^2*y(k)

The above expression is a nonvectorized computation and tends
to be slow. Specifying a

Vectorized computation uses MATLAB

vectorization rules to evaluate the regressor expression using matrices
instead of the

% ".*" indicates element-wise operation
z=(x.^2).*y

FOR-loop and results in faster computation:

z=@(x,y)x^2*y. x

Construction C=customreg(Function,Variables) specifies a custom regressor for

a nonlinear ARX model.
regressor.

Function is a handle or string representing a function

of input and output variables.
that represent the names of model inputs and outputs in the function

Function. Each input and output name must coincide with the strings

in the

idnlarx object. The size of Variables must match the number of
Function inputs. For multiple-output models with p outputs, the

InputName and OutputName properties of the corresponding

custom regressor is a
object, where the kyth entry defines the custom regressor for output

ky. You must add these re gres sors to the model by assigning the
CustomRegressors model property or by using addreg.

C=customreg(Function,Variables,Dela ys,Vectorized) create a

custom regressor that includes the delays corresponding to inputs or
outputs in

Arguments. Delays is a vector of positive integers that

represent the delays of
vector element). The size of

C is a customreg object that stores custom

Variables is a cel l array of strings

p-by-1 cell array or an array of customreg

Variables variables (default is 1 for each

Delays must match the size of Variables.

2-63

customreg

Vectorized value of 1 uses MATLAB vectorization rules to evaluate

the regressor expression
(false).

Properties After creating the object, you can use get or dot notation to access the

object property values. For example:

% List all property values
get(C)
% Get value of Arguments property
C.Arguments

You can also use the set function to set the value of particular
properties. For example:

set(C,'Vectorized',1)

Function.Bydefault,Vectorized value is 0

Property Name

Function

Variables

Description

Function handle or string representing a function of standards
regressors.

For example:

cr = @(x,y) x*y

Cell array of strings that represent the names of model input
and output variables in the function
output name must coincide with the strings in the

Function. Each input and

InputName

and OutputName properties of the idnlarx object—the model
for which you define custom regressors. The size of

Variables

must match the number of Function inputs.

For example,

C = customreg(cr,{'y1','u1'},[2 3])

Variables correspond to {'y1','u1'} in:

2-64

customreg

Property Name

Delays

Vectorized

Description

Vector of positive integers representing the delays of

Variables.ThesizeofDelays must match the size of
Arguments.

Default:

For example,

1 foreachvectorelement.

Delays are [2 3] in:

C = customreg(cr,{'y1','u1'},[2 3])

Assignable values:

•

0 (default)—Function is not computed in vectorized form.

1—Function is computed in vectorized form when called

•

with vector arguments.

Examples Define custom regressors as a cell array of strings:

load iddata1
m = nlarx(z1,[2 2 1]);
C={'u1(t-1)*sin(y1(t-3))','u1(t-2)^3'};
% u1 and y1 are system input and output

m.CustomRegressors = C;
m=pem(z1,m)

Define custom regressors directly in the estimation command nlarx:

m = nlarx(data,[na nb nk],'linear',...

'CustomRegressors',...
{'u1(t-1)*sin(y1(t-3))','u1(t-2)^3'});

Define custom regressors as an object array of customreg objects:

2-65

customreg

cr1=@(x,y) x*sin(y);
cr2=@(x) x^3;
C=[customreg(cr1,{'u' 'y'},[1 3]),...

customreg(cr2,{'u'},2)];

m=addreg(m,C);

Use vectorization rules to evaluate regressor expression during
estimation:

C = customreg(@(x,y) x*sin(y),{'u' 'y'},[1 3])
set(C,'Vectorized',1)
m = nlarx(data,[na nb nk],'sigmoidnet','CustomReg',C)

See Also addreg | getreg | idnlarx | nlarx | polyreg

How To • “Identifying Nonlinear A RX Models”

2-66

Purpose Transform linear model from continuous to discrete time

Syntax md = c2d(mc,T)

md = c2d(mc,T,method)
[md,G] = c2d(mc,T,method)

Description mc is a continuous-time model such as any idmodel object (idgrey,

idproc, idpoly,oridss). md is the model that i s obtained when it is

sampled with sampling interval

method = 'zoh' (default) makes the translation to discrete time under

the assumption that the input is piecewise constant (zero-order hold).

method = 'foh' assumes the input to be piecewise linear between the

sampling instants (first-order hold).

When you have the Control System Toolbox product installed, you can
also use the following methods:
cases, no translation of the covariance matrix takes place.

Note that the innovations variance λ of the continuous-time model
is interpreted as the intensity of the spectral density of the noise
spectrum. The noise variance in

T.

'tustin',and'matched'.Inthese

md is thus given as λ/T.

c2d

idpoly and idss models are returned in the same format. idgrey

structures are preserved if their CDMfile property is equal to 'cd'.
Otherwise they are transformed to
returned as

idpoly models, the covariance matrix is translated by the use of

For
numerical derivatives. The step sizes used for the differentiation are
given by the function
thecovariancematrixisnottranslated, but covariance information
about the input-output properties is included in
translation of covariance information (which may take some time), use

c2d(mc,T,'covariance','none').

The output argument G is a matrix that transforms the initial state x0
of mc to the initial state of md as

idgrey objects.

nuderst.Foridss, idproc, and idgrey models,

idss objects. idproc models are

md. To inhibit the

2-67

c2d

X0d=G * [X0; u(0)],

where u(0) is the input at time 0. For idproc models, the state
variables correspond to those of
is returned as the empty matrix.

Examples Define a continuous-time system and study the poles and zeros o f the

sampled counterpart.

mc = idpoly(1,1,1,1,[1 1 0],'Ts',0);
md = c2d(mc,0.5);
pzmap(md)

idgrey(mc).Foridpoly models, G

See Also

d2c

2-68

data2state(idnlarx)

Purpose Map past input/output data to current states of nonlinear ARX model

Syntax X = data2state(MODEL,IOSTRUCT)

X = data2state(MODEL,DATA)

Description X = data2state(MODEL,IOSTRUCT) maps the input and output samples

in

IOSTRUCT to th e current states of MODEL, X. For a definition of the

states of
2-189. The data in
that the returned state values must be interpreted as values at the
time immediately after the time corresponding to the last (most recent)
sample in the data.

X = data2state(MODEL,DATA) maps the input and output samples

from

Input • MODEL: idnlarx model.

IOSTRUCT: Structure with fields Input and Output.Samplesin

•

IOSTRUCT mustbeintheorderofincreasingtime(thelastrowof

values corresponds to the most recent time). Each field contains
data samples corresponding to the past input and output of
respectively.

idnlarx models, see “Defi n ition of idnl arx States” on page

IOSTRUCT is interpreted as past samples of data, so

DATA to the current states, X,ofthemodel.

MODEL

- Input:MatrixofNU columns, where NU is the number of inputs.

The number of rows can be equal to either of the following:

• Maximum input delay in MODEL (maximum across all input

variables).

• 1 to specify steady-state (constant) input values.

- Output:MatrixofNY columns, where NY is the number of outputs.

The number of rows can be equal to either of the following:

• Maximum input delay in MODEL (max imum across all output

variables).

• 1 to specify steady-state (constant) output values.

2-69

data2state(idnlarx)

• DATA: iddata object containing data samples. Samples in DATA must

be in the order of increasing time (the last row of values corresponds
to the most recent time). The number of samples in
greater than or equal to the maximum delay in the model across
all input and output variables.

Note To determine maximum delay in each input and output channel
of

MODEL,usethegetDelayInfo command. For more information, see

the

getDelayInfo reference page.

Output X is the state vector of MODEL correspondingtothetimeafterthemost

recent sample in the input data (

IOSTRUCT or DATA).

Examples In this example you determine the current state of an idnlarx model.

1 Load your data and create a data object.

DATA must be

2-70

load motorizedcamera;
z = iddata(y,u,0.02,'Name','Motorized Camera', ...

'TimeUnit','s');

2 Estimate an idnlarx model from the data. The model has 6 inputs

and 2 outputs.

mw1 = nlarx(z,[ones(2,2),ones(2,6),ones(2,6)],wavenet);

3 Compute the maximum delays across all output variables in mw1.

MaxDelays = getDelayInfo(mw1);

4 Represent the past input and output samples:

IOData = struct('Input', ...

rand(max(MaxDelays(3+1:end)),6),...

'Output', ...

data2state(idnlarx)

rand(max(MaxDelays(1:3)),2));

5 Compute the current states of mw1 based on the past data in

IOSTRUCT.

X = data2state(mw1,IOData)

The previous command computes the state vector.

Note You can specify constant input levels with scalar values

(10,20,30,40,50,60) for the input variables by setting
IOSTRUCT.Input = [10, 20, 30, 40, 50, 60] instead of a matrix of

values.

See Also

findop(idnlarx)

findstates(idnlarx)

getDelayInfo

2-71

deadzone

Purpose Class representing dead-zone nonlinearity estimator for

Hammerstein-Wiener models

Syntax s=deadzone(ZeroInterval,I)

Description deadzone is an object that stores the dead-zone nonlinearity estimator

for estimating Hammerstein-Wiener models.

You can use the constructor to create the nonlinearity object, as follows:

s=deadzone(ZeroInterval,I) creates a dead-zone nonlinearity

estimator object, initialized with the zero interval

evaluate(d,x) to compute the value of the function defined by

Use
the

deadzone object d at x.

I.

Remarks Use deadzone to define a nonlinear function

function of x and has the following characteristics:

axb Fx

≤< =

xa Fx xa

<=−

xb

≥

y and x are scalars.

()()
Fx x b()=−

0

yFx= ()

,whereF is a

Properties You can specify the property value as an argument in the constructor

to specify the object.

After creating the object, you can use
object property values. For example:

% List ZeroInterval property value
get(d)
d.ZeroInterval

You can also use the set function to set the value of particular
properties. For example:

get or dot notation to access the

2-72

deadzone

set(d, 'ZeroInterval', [-1.5 1.5])

The first argument to set must be the name of a MATLAB variable.

Property Name

ZeroInterval

Description

1-by-2 row vector that specifies the initial zero interval of the
nonlinearity.
Default=

For example:

[NaN NaN].

deadzone('ZeroInterval',[-1.5 1.5])

Examples Use deadzone to specify the dead-zone nonlinearity estimator in

Hammerstein-Wiener models. For example:

m=nlhw(Data,Orders,deadzone([-1 1]),[]);

The dead-zone nonlinearity is initialized at the interval [-1 1].The

See Also

interval values are adjusted to the estimation data by

nlhw

nlhw.

2-73

delayest

Purpose Estimate time delay (dead time) from data

Syntax nk = delayest(Data)

nk = delayest(Data,na,nb,nkmin,nkmax,maxtest)

Description Data is an iddata object containing the input-output data. It can also be

an

idfrd object defining frequency-response data. Only single-output

data can be handled.

nk is returned as an integer or a row vector of integers, containing the

estimated time delay in samples from the input(s) to the output in

The estimate is based on a comparison of ARX models with different
delays:

yt a yt a yt na

() ( ) ... ( )

+ −++ − =

b u t nk b u t nb nk

( ) ... ( )

−++ −−+

1

The integer na is the order of the A polynomial (default 2). nb is a row
vector of length equal to the number of inputs, containing the order(s) of
the B polynomial(s) (default all 2).

1

1

nb

na

1 ++ et()

Data.

2-74

nkmin and nkmax are row vectors of the same length as the number of

inputs, contain ing the smallest and largest delays to be tested. Defaults
are

nkmin = 0 and nkmax = nkmin+20.

nb, nkmax,and/ornkmin are entered as scalars in the multiple-input

If
case, all inputs will be assigned the same values.

maxtest is the largest number of tests allowed (default 10,000).

Purpose Subtract offset or trend from data signals

Syntax data_d = detrend(data)

data_d = detrend(data,Type)
[data_d,T] = detrend(data,Type)
data_d = detrend(data,1,brkp)

Description data_d = detrend(data) subtracts the mean value from each

time-domain or time-series signal
objects.

data_d = detrend(data,Type) subtracts a mean value from each

signal when

,oratrendspecifiedbyaTrendInfo object when Type = T.

1

[data_d,T] = detrend(data,Type) stores the trend information as a
TrendInfo object T.

data_d = detrend(data,1,brkp) subtracts a piecewise linear

Type = 0, a linear trend (least-squares fit) when Type =

trend at one or more breakpoints
discontinuities between successive linear trends occur. When
contains breakpoints that match the time vector, detrend subtracts a
continuous piecewise linear trend. You cannot store piecewise linear
trend information as an output argument.

data. data_d and data are iddata

brkp. brkp is a data index where

detrend

brkp

Examples Subtract mean values from input and output signals and store the

trend information:

% Load SISO data containing vectors u2 and y2.
load dryer2
% Create data object with sampling interval of 0.08 sec.
data=iddata(y2,u2,0.08)
% Plot data on a time plot. Data has a nonzero mean.
plot(data)
% Remove the mean from the data.
[data_d,T] = detrend(data,0)
% Plot detrended data on the same plot.
hold on

2-75

detrend

plot(data_d)

Remove specified offset from input and output signals:

% Load SISO data containing vectors u2 and y2.
load dryer2
% Create data object with sampling time of 0.08 sec.
data=iddata(y2,u2,0.08)
plot(data)
% Create a TrendInfo object for storing offsets and trends.
T = getTrend(data)
% Assign offset values to the TrendInfo object.
T.InputOffset=5;
T.OutputOffset=5;
% Subtract offset from the data.
data_d = detrend(data,T)
% Plot detrended data on the same plot.
hold on
plot(data_d)

Subtract several linear trends that connect at three breakpoints [30

60 90]

:

data = detrend(data,1,[30 60 90]);
% [30 60 90] are data indexes where breakpoints occur.

Subtract a mean value from the input signal and a V-shaped trend
from the output signal, such that the V peak occurs at the breakpoint
value of

119:

zd1 = z(:,:,[]); zd2 = z(:,[],:);
zd1(:,1,[]) = detrend(z(:,1,[]),1,119);
zd2(:,[],1) = detrend(z(:,[],1));
zd = [zd1,zd2];

See Also getTrend | retrend | TrendInfo

How To • “Handling Offsets and Trends in Data”

2-76