Identifiaction With Matlab PDF
Identifiaction With Matlab PDF
Identifiaction With Matlab PDF
Toolbox
For Use with MATLAB
Lennart Ljung
Computation
Visualization
Programming
Users Guide
Version 5
Web
Newsgroup
Technical support
Product enhancement suggestions
Bug reports
Documentation error reports
Order status, license renewals, passcodes
Sales, pricing, and general information
508-647-7000
Phone
508-647-7001
Fax
[email protected]
[email protected]
[email protected]
[email protected]
[email protected]
For contact information about worldwide offices, see the MathWorks Web site.
System Identification Toolbox Users Guide
COPYRIGHT 1988 - 2002 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 MathWorks, Inc.
FEDERAL ACQUISITION: This provision applies to all acquisitions of the Program and Documentation by
or for the federal government of the United States. By accepting delivery of the Program, the government
hereby agrees that this software qualifies as "commercial" computer software within the meaning of FAR
Part 12.212, DFARS Part 227.7202-1, DFARS Part 227.7202-3, DFARS Part 252.227-7013, and DFARS Part
252.227-7014. The terms and conditions of The MathWorks, Inc. Software License Agreement shall pertain
to the governments use and disclosure of the Program and Documentation, and shall supersede any
conflicting contractual terms or conditions. If this license fails to meet the governments minimum needs or
is inconsistent in any respect with federal procurement law, the government agrees to return the Program
and Documentation, unused, to MathWorks.
MATLAB, Simulink, Stateflow, Handle Graphics, and Real-Time Workshop are registered trademarks, and
TargetBox is a trademark of The MathWorks, Inc.
Other product or brand names are trademarks or registered trademarks of their respective holders.
First printing
Second printing
Third printing
Fourth printing
Fifth printing
Online only
Contents
Preface
1
Basic Questions About System Identification . . . . . . . . . . . . 1-2
Common Terms Used in System Identification . . . . . . . . . . 1-4
Basic Information About Dynamic Models . . . . . . . . . . . . . . 1-6
The Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-6
The Basic Dynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-7
Variants of Model Descriptions . . . . . . . . . . . . . . . . . . . . . . . . . 1-7
How to Interpret the Noise Source . . . . . . . . . . . . . . . . . . . . . . . 1-8
Terms to Characterize the Model Properties . . . . . . . . . . . . . . 1-10
The Basic Steps of System Identification . . . . . . . . . . . . . . . 1-12
A Startup Identification Procedure . . . . . . . . . . . . . . . . . . . . 1-14
Step 1: Looking at the Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 1-14
Step 2: Getting a Feel for the Difficulties . . . . . . . . . . . . . . . . 1-14
2
The Big Picture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Model and Data Boards . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Working Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Views . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Validation Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Work Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Management Aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Workspace Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Help Texts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2-2
2-2
2-3
2-3
2-4
2-4
2-4
2-5
2-6
ii
Contents
2-15
2-15
2-15
2-16
2-17
2-21
2-23
2-25
2-26
Examining Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Views and Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Plot Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Frequency Response and Disturbance Spectra . . . . . . . . . . . .
Transient Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Poles and Zeros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Compare Measured and Model Output . . . . . . . . . . . . . . . . . . .
Residual Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Text Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
LTI Viewer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Further Analysis in the MATLAB Workspace . . . . . . . . . . . . .
2-28
2-28
2-29
2-30
2-30
2-31
2-31
2-32
2-33
2-33
2-34
2-35
2-35
2-36
2-36
2-37
2-37
Tutorial
3
Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-2
The Toolbox Commands . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-3
An Introductory Example to Command Mode . . . . . . . . . . . . 3-5
The System Identification Problem . . . . . . . . . . . . . . . . . . . . . 3-9
Impulse Responses, Frequency Functions, and Spectra . . . . . . 3-9
Polynomial Representation of Transfer Functions . . . . . . . . . 3-11
State-Space Representation of Transfer Functions . . . . . . . . . 3-13
Continuous-Time State-Space Models . . . . . . . . . . . . . . . . . . . 3-14
Estimating Impulse Responses . . . . . . . . . . . . . . . . . . . . . . . . . 3-15
Estimating Spectra and Frequency Functions . . . . . . . . . . . . . 3-15
Estimating Parametric Models . . . . . . . . . . . . . . . . . . . . . . . . . 3-16
Subspace Methods for Estimating State-Space Models . . . . . . 3-17
iii
iv
Contents
3-19
3-19
3-20
3-20
3-22
3-26
3-27
3-27
3-28
3-30
3-31
3-37
3-38
3-39
3-41
Examining Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Parametric Models: idmodel and its children . . . . . . . . . . . . . .
Frequency Function Format: the idfrd model . . . . . . . . . . . . .
Graphs of Model Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Transformations to Other Model Representations . . . . . . . . .
Discrete and Continuous Time Models . . . . . . . . . . . . . . . . . . .
3-51
3-51
3-57
3-58
3-61
3-62
3-65
3-65
3-68
3-68
3-68
3-69
3-70
3-70
3-44
3-46
3-49
3-76
3-76
3-76
3-77
3-77
3-78
3-79
3-80
3-80
3-81
3-83
3-85
Function Reference
4
Functions By Category . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Help Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The Graphical User Interface . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulation and Prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Data Manipulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Nonparametric Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Model Structure Creation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Manipulating Model Structures . . . . . . . . . . . . . . . . . . . . . . . . .
Model Conversions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Model Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Assessing Model Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . .
Model Structure Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Recursive Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . .
General . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vi
Contents
4-2
4-2
4-2
4-2
4-2
4-3
4-3
4-3
4-5
4-5
4-5
4-7
4-7
4-7
4-8
4-8
Preface
Preface
xi
Preface
Typographical Conventions
This manual uses some or all of these conventions.
Item
Convention
Example
Example code
Monospace font
Monospace font
f = freqspace(n,'whole')
Mathematical
expressions
MATLAB output
Monospace font
xii
Italics
An array is an ordered
collection of information.
[c,ia,ib] = union(...)
Monospace italics
sysc = d2c(sysd,'method')
Related Products
Related Products
The MathWorks provides several products that are especially relevant to the
kinds of tasks you can perform with the System Identification Toolbox. In
particular, the Systems Identification Toolbox requires these products:
MATLAB
For more information about any of these products, see either:
The online documentation for that product, if it is installed or if you are
reading the documentation from the CD
The MathWorks Web site, at http://www.mathworks.com; see the products
section
Note The products listed below complement the functionality of the System
Identification Toolbox.
Product
Description
Simulink
Financial Toolbox
xiii
Preface
xiv
Product
Description
Optimization Toolbox
Signal Processing
Toolbox
Statistics Toolbox
xv
Preface
xvi
1
The System Identification
Problem
1-2
1-3
1-4
1-5
The Signals
Models describe relationships between measured signals. It is convenient to
distinguish between input signals and output signals. The outputs are then
partly determined by the inputs. Think for example of an airplane where the
inputs would be the different control surfaces, ailerons, elevators, and the like,
while the outputs would be the airplanes orientation and position. In most
cases, the outputs are also affected by more signals than the measured inputs.
In the airplane example it would be wind gusts and turbulence effects. Such
unmeasured inputs will be called disturbance signals or noise. If we denote
inputs, outputs, and disturbances by u, y, and e, respectively, the relationship
can be depicted in the following figure.
e
All these signals are functions of time, and the value of the input at time t will
be denoted by u(t). Often, in the identification context, only discrete-time points
are considered, since the measurement equipment typically records the signals
just at discrete-time instants, often equally spread in time with a sampling
interval of T time units. The modeling problem is then to describe how the
three signals relate to each other.
1-6
( ARX )
Such a relationship tells us, for example, how to compute the output y(t) if the
input is known and the disturbance can be ignored:
y ( t ) = 1.5y ( t T ) 0.7y ( t 2T ) + 0.9u ( t 2T ) + 0.5u ( t 3T )
The output at time t is thus computed as a linear combination of past outputs
and past inputs. It follows, for example, that the output at time t depends on
the input signal at many previous time instants. This is what the word
dynamic refers to. The identification problem is then to use measurements of
u and y to figure out:
The coefficients in this equation (i.e., -1.5, 0.7, etc.).
How many delayed outputs to use in the description (two in the example:
y (t-T) and y (t-2T) ) .
The time delay in the system is (2T in the example: you see from the second
equation that it takes 2T time units before a change in u will affect y).
How many delayed inputs to use (two in the example: u(t-2T) and u(t-3T)).
The number of delayed inputs and outputs are usually referred to as the
model order(s).
y=G u+He
1-7
which says that the measured output y(t) is a sum of one contribution that
comes from the measured input u(t) and one contribution that comes from the
noise He. The symbol G then denotes the dynamic properties of the system,
that is, how the output is formed from the input. For linear systems it is called
the transfer function from input to output. The symbol H refers to the noise
properties, and is called the disturbance model. It describes how the
disturbances at the output are formed from some standardized noise source
e(t).
State-space models are common representations of dynamical models. They
describe the same type of linear difference relationship between the inputs and
the outputs as in the ARX model, but they are rearranged so that only one
delay is used in the expressions. To achieve this, some extra variables, the
state variables, are introduced. They are not measured, but can be
reconstructed from the measured input-output data. This is especially useful
when there are several output signals, i.e., when y(t) is a vector. Chapter 3,
Tutorial, gives more details about this. For basic use of the toolbox it is
sufficient to know that the order of the state-space model relates to the
number of delayed inputs and outputs used in the corresponding linear
difference equation. The state-space representation looks like
x(t+1)=Ax( t) +Bu(t)+Ke(t)
y(t)=Cx(t)+D u(t)+e(t)
Here x(t) is the vector of state variables. The model order is the dimension of
this vector. The matrix K determines the disturbance properties. Notice that if
K = 0, then the noise source e(t) affects only the output, and no specific model
of the noise properties is built. This corresponds to H = 1 in the general
description above, and is usually referred to as an Output-Error model. Notice
also that D = 0 means that there is no direct influence from u(t) to y(t). Thus
the effect of the input on the output all passes via x(t) and will thus be delayed
at least one sample. The first value of the state variable vector x(0) reflects the
initial conditions for the system at the beginning of the data record. When
dealing with models in state-space form, a typical option is whether to estimate
D, K, and x(0) or to let them be zero.
1-8
1-9
The need and use of the noise model can be summarized as follows:
It is, in most cases, required to obtain a better estimate for the dynamics, G.
It indicates how reliable noise-free simulations are.
It is required for reliable predictions and stochastic control design.
Impulse Response
The impulse response of a dynamical model is the output signal that results
when the input is an impulse, i.e., u(t) is zero for all values of t except t=0,
where u(0)=1. It can be computed as in the equation following (ARX), by letting
t be equal to 0, 1, 2, ... and taking y(-T)=y(-2T)=0 and u(0)=1.
Step Response
The step response is the output signal that results from a step input, i.e., u(t)
is zero for negative values of t and equal to one for positive values of t. The
impulse and step responses together are called the models transient
response.
Frequency Response
The frequency response of a linear dynamic model describes how the model
reacts to sinusoidal inputs. If we let the input u(t) be a sinusoid of a certain
frequency, then the output y(t) will also be a sinusoid of this frequency. The
amplitude and the phase (relative to the input) will however be different. This
frequency response is most often depicted by two plots; one that shows the
amplitude change as a function of the sinusoids frequency and one that shows
the phase shift as function of frequency. This is known as a Bode plot.
1-10
1-11
identified.
2 Examine the data. Polish it so as to remove trends and outliers, and select
another model set. Possibly also try other estimation methods (Step 4) or
work further on the input-output data (Steps 1 and 2).
1-12
The System Identification Toolbox offers several functions for each of these
steps.
For Step 2 there are routines to plot data, filter data, and remove trends in
data, as well as to resample and reconstruct missing data.
For Step 3 the System Identification Toolbox offers a variety of nonparametric
models, as well as all the most common black-box input-output and state-space
structures, and also general tailor-made linear state-space models in discrete
and continuous time.
For Step 4 general prediction error (maximum likelihood) methods, as well as
instrumental variable methods and sub-space methods are offered for
parametric models, while basic correlation and spectral analysis methods are
used for nonparametric model structures.
To examine models in Step 5, many functions allow the computation and
presentation of frequency functions and poles and zeros, as well as simulation
and prediction using the model. Functions are also included for
transformations between continuous-time and discrete-time model
descriptions and to formats that are used in other MATLAB toolboxes, like the
Control System Toolbox and the Signal Processing Toolbox.
1-13
1-14
Measured Validation Data output and the ARX and state-space models
simulated outputs
If these agreements are reasonable, the problem is not so difficult, and a
relatively simple linear model will do a good job. Some fine tuning of model
orders, and noise models have to be made and you can proceed to Step 4.
Otherwise go to Step 3.
Model Unstable
The ARX or state-space model may turn out to be unstable, but could still be
useful for control purposes. Change to a 5- or 10-step ahead prediction instead
of simulation in the Model Output View.
Feedback in Data
If there is feedback from the output to the input, due to some regulator, then
the spectral and correlations analysis estimates are not reliable. Discrepancies
between these estimates and the ARX and state-space models can therefore be
disregarded in this case. In the Model Residuals View of the parametric
models, feedback in data can also be visible as correlation between residuals
and input for negative lags.
Disturbance Model
If the state-space model is clearly better than the ARX model at reproducing
the measured output, this is an indication that the disturbances have a
substantial influence, and it will be necessary to model them carefully.
Model Order
If a fourth order model does not give a good Model Output plot, try eighth
order. If the fit clearly improves, it follows that higher order models will be
required, but that linear models could be sufficient.
Additional Inputs
If the Model Output fit has not significantly improved by the tests so far, think
over the physics of the application. Are there more signals that have been, or
1-15
could be, measured that might influence the output? If so, include these among
the inputs and try again a fourth order ARX model from all the inputs. (Note
that the inputs need not at all be control signals, anything measurable,
including disturbances, should be treated as inputs).
Nonlinear Effects
If the fit between measured and model output is still bad, consider the physics
of the application. Are there nonlinear effects in the system? In that case, form
the nonlinearities from the measured data and add those transformed
measurements as extra inputs. This could be as simple as forming the product
of voltage and current measurements, if you realize that it is the electrical
power that is the driving stimulus in, say, a heating process, and temperature
is the output. This is of course application dependent. It does not take very
much work, however, to form a number of additional inputs by reasonable
nonlinear transformations of the measured ones, and just test if inclusion of
them improves the fit.
Still Problems?
If none of these tests leads to a model that is able to reproduce the Validation
Data reasonably well, the conclusion might be that a sufficiently good model
cannot be produced from the data. There may be many reasons for this. It may
be that the system has some quite complicated nonlinearities, which cannot be
realized on physical grounds. In such cases, nonlinear, black-box models could
be a solution. Among the most used models of this character are the Artificial
Neural Networks (ANN).
Another important reason is that the data simply do not contain sufficient
information, e.g., due to bad signal to noise ratios, large and nonstationary
disturbances, varying system properties, etc.
Otherwise, use the insights of which inputs to use and which model orders to
expect and proceed to Step 4.
1-16
properties of the obtained models. There are a few things to look for in these
comparisons.
1-17
using the colon (:) notation. You can also press the Order Selection button.
When you select Estimate, models for all combinations (easily several
hundreds) are computed and their (prediction error) fit to Validation Data is
shown in a special plot. By clicking in this plot the best models with any chosen
number of parameters will be inserted into the Model Board, and evaluated as
desired.
Many State-space models: A similar feature is also available for black-box
state-space models, estimated using n4sid. When a good order has been found,
try the PEM estimation method, which often improves on the accuracy.
ARMAX, OE, and BJ models: Once you have a feel for suitable delays and
dynamics orders, if is often useful to try out ARMAX, OE, and/or BJ with these
orders and test some different orders for the disturbance transfer functions (C
and D). Especially for poorly damped systems, the OE structure is suitable.
Multivariable Systems
Systems with many input signals and/or many output signals are called
multivariable. Such systems are often more challenging to model. In particular
systems with several outputs could be difficult. A basic reason for the
difficulties is that the couplings between several inputs and outputs lead to
more complex models. The structures involved are richer and more parameters
will be required to obtain a good fit.
Available Models
The System Identification Toolbox as well as the GUI handle general, linear
multivariable models. All earlier mentioned models are supported in the single
output, multiple input case. For multiple outputs, ARX models and state-space
models are covered. Multi-output ARMAX and OE models are covered via
state-space representations: ARMAX corresponds to estimating the K-matrix,
while OE corresponds to fixing K to zero. (These are pop-up options in the GUI
model order editor.)
Generally speaking, it is preferable to work with state-space models in the
multivariable case, since the model structure complexity is easier to deal with.
It is essentially just a matter of choosing the model order.
1-18
1-19
menu Preprocess > Select Channels for this. Dont change the Validation
Data. The GUI will keep track of the input and output channels. It will do
the right thing when evaluating the channel-restricted models using the
Validation Data. It might also be appropriate to see if improvements in the
fit are obtained for various model types, built for one output at a time.
If you decide for a multi-output model, it is often easiest to use state-space
models. Use n4sid as a primary tool and try out pem when a good order has
been found.
1-20
1-21
1-22
2
The Graphical User
Interface
2-2
The Views
Below the Data and Model Boards are buttons for different views. These
control what aspects of the data sets and models you would like to examine, and
are described in more detail in Handling Data on page 2-7 and in Examining
Models on page 2-28. To select a data set or a model, so that its properties are
displayed, click on its icon. A selected object is marked by a thicker line in the
icon. To deselect, click again. An arbitrary number of data/model objects can be
examined simultaneously. To have more information about an object,
double-click (or right-click or Ctrl-click) on its icon.
2-3
Management Aspects
Diary: It is easy to forget what you have been doing. By double-clicking on a
data/model icon, a complete diary will be given of how this object was created,
along with other key information. At this point you can also add comments and
change the name of the object and its color.
Layout: To have a good overview of the created models and data sets, it is good
practice to try rearranging the icons by dragging and dropping. In this way
models corresponding to a particular data set can be grouped together, etc. You
can also open new boards (Options menu Extra model/data boards) to further
rearrange the icons. These can be dragged across the screen between different
windows. The extra boards are also equipped with notepads for your comments.
2-4
Sessions: The Model and Data Boards with all models and data sets together
with their diaries can be saved (under menu item File) at any point, and
reloaded later. This is the counterpart of save/load workspace in the
command-driven MATLAB. The four most recent sessions are listed under File
for immediate open.
Cleanliness: The boards will hold an arbitrary number of models and data sets
(by creating clones of the board when necessary). It is however advisable to
clear (delete) models and data sets that no longer are of interest. Do that by
dragging the object to the Trash Can. (Double-clicking on the trash can will
open it up, and its contents can be retrieved.) Empty the can if you run into
memory problems.
Window Culture: Dialog and plot windows are best managed by the GUIs
close function (submenu item under File menu, or select Close, or
Workspace Variables
The models and data sets created within the GUI are normally not available in
the MATLAB workspace. Indeed, the workspace is not at all littered with
variables during the sessions with the GUI. The variables can however be
exported at any time to the workspace, by dragging and dropping the object
icon on the To Workspace box. They will then carry the name in the workspace
that marked the object icon at the time of export. You can work with the
variables in the workspace, using any MATLAB commands, and then perhaps
import modified versions back into the GUI. Note that models and data are
exported as the toolboxs objects idmodel, idfrd, and iddata. For how to
extract information and work with these objects, see Chapter 3, Tutorial and
Model Conversions on page 4-5 of the Command Reference chapter.
The GUIs names of data sets and models are suggested by default procedures.
Normally, you can enter any other name of your choice at the time of creation
of the variable. Names can be changed (after double-clicking on the icon) at any
time. Unlike the workspace situation, two GUI objects can carry the same
name (i.e., the same string in their icons).
2-5
Help Texts
The GUI contains some 100 help texts that are accessible in a nested fashion,
when required. The main ident window contains general help topics under the
Help menu. This is also the case for the various plot windows. In addition,
every dialog box has a Help push button for current help and advice.
2-6
Handling Data
Handling Data
Data Representation
In the System Identification Toolbox, signals and observed data are
represented as column vectors, e.g.,
u(1)
u(2)
u =
u(N )
The entry in row number k, i.e., u(k), will then be the signals value at sampling
instant number k. It is generally assumed in the toolbox that data are sampled
at equidistant sampling times, and the sampling interval T is supplied as a
specific argument.
We generally denote the input to a system by the letter u and the output by y.
If the system has several input channels, the input data is represented by a
matrix, where the columns are the input signals in the different channels:
u = u1 u2 um
The same holds for systems with several output channels.
The observed input-output data record is represented in the System
Identification Toolbox by the iddata object, that is created from the input and
output signals by
Data = iddata(y,u,Ts)
2-7
In addition to this mandatory information, you may add further properties that
will help in the bookkeeping:
4 The starting time for the sampling
5 Input and output channel names
6 Input and output channel units
7 Periodicity and intersample behavior of the input
8 Data notes: These are notes for your own information and bookkeeping that
will follow the data and all models created from them.
2-8
Handling Data
As you select the pop-up menu Data and choose the item Import, a dialog box
will open, where you can enter the information items 1 - 8, just listed. This box
has five fields for you to fill in.
Figure 2-2: The Dialog for Importing Data into the GUI
Actually, you can enter any MATLAB expressions in these fields, and they will
be evaluated to compute the input and the output before inserting the data into
the GUI.
Data name: Enter the name of the data set to be used by the GUI. This name
can be changed later on.
Starting time and Sampling interval: Fill these out for correct time and
frequency scales in the plots.
2-9
Channel names: Enter strings for the different input and output channels
names. Separate the strings by comma. The number of names must be equal to
the number of channels. If these entries are not filled out, default names, y1,y2,
..., u1, u2 ..., will be used.
Channel units: Enter, in analogous format, the units in which the
measurements are made. These will follow to all models built from data, but
are used only for plot information.
Period: If the input is periodic, enter here the period length. Inf means a
non-periodic input, which is default.
Intersample: Choose the intersample behavior of the input as one of ZOH
(zero-order hold, i.e., the input signal piecewise constant between the samples)
or FOH (first-order hold, i.e., the input signal is piecewise linear between the
samples) or BL (Band-limited, i.e., the continuous time input signal has no
power above the Nyquist frequency). ZOH is default.
The box at the bottom is for Notes, where you can enter any text you want to
accompany the data for bookkeeping purposes.
Finally, select Import to insert the data into the GUI. When no more data sets
are to be inserted, select Close to close the dialog box. Reset will empty all the
fields of the box.
The procedure just described will create an iddata object, with all its
properties. If you already have an iddata object available in the workspace,
you can import that directly by selecting the data format Iddata Object in the
pop-up menu at the top of the Import Data dialog.
2-10
Handling Data
shown instead. By default the periodograms of the data are shown, i.e., the
absolute square of the Fourier transforms of the data. The plot can be changed
to any chosen frequency range and a number of different ways of estimating
spectra, by the Options menu item in the spectra window.
The purpose of examining the data in these ways is to find out if there are
portions of the data that are not suitable for identification, if the information
contents of the data is suitable in the interesting frequency regions, and if the
data have to be preprocessed in some way, before using them for estimation.
Preprocessing Data
Detrending
Detrending the data involves removing the mean values or linear trends from
the signals (the means and the linear trends are then computed and removed
from each signal individually). This function is accessed under the pop-up
menu Preprocess, by selecting item Remove Means or Remove Trends. More
advanced detrending, such as removing piecewise linear trends or seasonal
variations cannot be accessed within the GUI. It is generally recommended
that you remove at least the mean values of the data before the estimation
phase. There are however situations when it is not advisable to remove the
sample means. It could for example be that the physical levels are built into the
underlying model, or that integrations in the system must be handled with the
right level of the input being integrated.
2-11
when you evaluate data and model properties, for models covering different
subsets of the data.
Prefiltering
By filtering the input and output signals through a linear filter (the same filter
for all signals) you can, e.g., remove drift and high frequency disturbances in
the data, that should not affect the model estimation. This is done by selecting
the pop-up menu item Preprocess > Filter... in the main window. The dialog
is quite analogous to that of selecting data ranges in the time domain. You
mark with a rectangle in the spectral plots the intended passband or stop band
of the filter, you select a button to check if the filtering has the desired effect,
and then you insert the filtered data into the GUIs Data Board.
Prefiltering is a good way of removing high frequency noise in the data, and
also a good alternative to detrending (by cutting out low frequencies from the
pass band). Depending on the intended model use, you can also make sure that
the model concentrates on the important frequency ranges. For a model that
will be used for control design, for example, the frequency band around the
intended closed-loop bandwidth is of special importance.
If you intend to use the data to build models both of the system dynamics and
the disturbance properties, it is recommended to do the filtering at the
estimation phase. That is achieved by selection the pop-up menu item
Estimate > Parametric Models, and then select the estimation Focus to be
Filter. This opens the same filter dialog as above. The prefiltering will
however apply only for estimating the dynamics from input to output. The
disturbance model is determined from the original data.
Resampling
If the data turn out to be sampled too fast, they can be decimated, i.e., every
k-th value is picked, after proper prefiltering (antialias filtering). This is
obtained from menu item Preprocess > Resample.
You can also resample at a faster sampling rate by interpolation, using the
same command, and giving a resampling factor less than one.
Quickstart
The pop-up menu item Preprocess > Quickstart performs the following
sequence of actions: It opens the Time plot Data view, removes the means
from the signals, and it splits these detrended data into two halves. The first
2-12
Handling Data
one is made Working Data and the second one becomes Validation Data. All the
three created data sets are inserted into the Data Board.
Multi-Experiment Data
The toolbox allows the handling of data sets that contain several different
experiments. Both estimation and validation can be applied to such data sets.
This is quite useful to deal with experiments that have been conducted at
different occasions but describe the same system. It is also useful to be able to
keep together pieces of data that have been obtained by cutting out
informative pieces from a long data set. Multi-experiment data can be
imported and used in the GUI as any iddata object. Selecting specific part of a
multi-experiment data set is done from the pop-up menu item Preprocess >
Select Experiment. To merge several data sets in the Data board (obtained,
e.g., by cutting out nice portions from other data sets) use the pop-up menu
item Preprocess > Merge Experiment.
Simulating Data
The GUI is intended primarily for working with real data sets, and does not
itself provide functions for simulating synthetic data. That has to be done in
command mode, and you can use your favorite procedure in Simulink, the
Signal Processing Toolbox, or any other toolbox for simulation and then insert
the simulated data into the GUI as described above.
The System Identification Toolbox also has several commands for simulation.
For example, should check idinput and sim in the Command Reference
chapter for details. The following example shows how the ARMAX model
y ( t ) 1.5y ( t 1 ) + 0.7y ( t 2 ) =
u ( t 1 ) + 0.5u ( t 2 ) + e ( t ) e ( t 1 ) + 0.2e ( t 1 )
2-13
The input, u, and the output, y, can now be imported into the GUI as data, and
the various estimation routines can be applied to them. By also importing the
simulation model, model1, into the GUI, its properties can be compared to those
of the different estimated models.
To simulate a continuous-time state-space model
x = Ax + Bu + Ke
y = Cx + e
with the same input, and a sampling interval of 0.1 seconds, do the following
in the System Identification Toolbox.
A = [-1 1;-0.5 0]; B = [1; 0.5]; C = [1 0]; D = 0; K = [0.5;0.5];
Model2 = idss(A,B,C,D,K,'Ts', 0) % Ts = 0 means continuous time
Data = iddata([],[u e]);
Data.Ts = 0.1
y=sim(Model2,Data);
2-14
Estimating Models
Estimating Models
The Basics
Estimating models from data is the central activity in the System
Identification Toolbox. It is also the one that offers the most variety of
possibilities and thus is the most demanding one for the user.
All estimation routines are accessed from the pop-up menu Estimate in the
ident window. The models are always estimated using the data set that is
currently in the Working Data box.
One can distinguish between two different types of estimation methods:
Direct estimation of the Impulse or the Frequency Response of the system.
These methods are often also called nonparametric estimation methods, and
do not impose any structure assumptions about the system, other than that
it is linear.
Parametric methods. A specific model structure is assumed, and the
parameters in this structure are estimated using data. This opens up a large
variety of possibilities, corresponding to different ways of describing the
system. Dominating ways are state-space and several variants of difference
equation descriptions.
y(t) =
gk u ( t k )
k=1
The name derives from the fact that if the input u(t) is an impulse, i.e., u(t)=1
when t=0 and 0 when t>0, then the output y(t) will be y ( t ) = g t . For a
multivariable system, the impulse response g k will be a ny-by-nu matrix,
where ny is the number of outputs and nu is the number of inputs. Its i-j
element thus describes the behavior of the i-th output after an impulse in the
j-th input.
By choosing menu item Estimate > Correlation Model impulse response
coefficients are estimated directly from the input/output data using so called
2-15
2-16
Estimating Models
different options on how to graph the curves. The frequencies for which to
estimate the response can also be selected as an option under the Options
menu in this View window.
The Spectral Analysis command also estimates the spectrum of the additive
disturbance v(t) in the system description. This estimated disturbance
spectrum is examined under the Model View item Noise Spectrum.
The Spectral Analysis estimate is stored as an idfrd object. If you need to
further work with the estimates, you can export the model to the MATLAB
workspace and retrieve the responses directly from this object or by Nyquist
and Bode. See idfrd, bode, and nyquist in the Command Reference chapter
for more information. (A model is exported by dragging and dropping it over the
To Workspace icon.)
Two options that affect the spectral analysis estimate can be set in the dialog
box. The most important choice is a number, M, (the size of the lag window)
that affects the frequency resolution of the estimates. Essentially, the
frequency resolution is about 2 /M radians/(sampling interval). The choice of
M is a trade-off between frequency resolution and variance (fluctuations). A
large value of M gives good resolution but fluctuating and less reliable
estimates. The default choice of M is good for systems that do not have very
sharp resonances and may have to be adjusted for more resonant systems.
The options also offer a choice between the Blackman-Tukey windowing
method spa (which is default) and a method based on smoothing direct Fourier
transforms, etfe. etfe has an advantage for highly resonant systems, in that
it is more efficient for large values of M. It however has the drawbacks that it
requires linearly spaced frequency values, does not estimate the disturbance
spectrum, and does not provide confidence intervals. The actual methods are
described in more detail in the Command Reference chapter under spa and
etfe. To obtain the spectral analysis model for the current settings of the
options, you can just type the hotkey s in the ident window.
2-17
Parametric Models, which contains the basic dialog for all parametric
2-18
Estimating Models
You can fill out the Order box yourself at any time, but for assistance you can
select Order Editor. This will open up another dialog box, depending on the
chosen Structure, in which the desired model order and structure information
can be entered in a simpler fashion.
You can also enter a name of a MATLAB workspace variable in the order edit
box. This variable should then have a value that is consistent with the
necessary orders for the chosen structure.
Note For the state-space structure and the ARX structure, several orders
and combination of orders can be entered. Then all corresponding models will
be compared and displayed in a special dialog window for you to select
suitable ones. This could be a useful tool to select good model orders. This
option is described in more detail later in this section. When it is available, a
button Order selection is visible.
Estimation Method
A common and general method of estimating the parameters is the prediction
error approach, where simply the parameters of the model are chosen so that
the difference between the models (predicted) output and the measured output
is minimized. This method is available for all model structures. Except for the
ARX case, the estimation involves an iterative, numerical search for the best
fit.
To obtain information from and interact with this search, select Iteration
control... This is a button which is visible when an iterative estimation process
has been selected. This also gives access to a number of options that govern the
search process. (See Algorithm Properties on page 4-14 in the Command
Reference chapter.)
For some model structures (the ARX model, and black-box state-space models)
methods based on correlation are also available: Instrumental Variable (IV)
and Sub-space (N4SID) methods. The choice between methods is made in the
Parametric Models dialog box.
The dialog box also has three pop-up menus that offer further options: Focus
allows you to choose between a frequency weighting that concentrates on the
models prediction or simulation performance. Another alternative is
prefiltering, which was described in Prefiltering. Moreover, the pop-up menu
2-19
InitialState gives options to estimate the initial state or to fix it to zero. The
value Auto makes an automatic choice among these options. Finally, the menu
Covariance allows the choice between Estimate and None. The normal
situation is that the covariance of the model is estimated, so that various
uncertainty measures can be displayed in the plots. However, for high order
state-space models estimated by N4SID, or large multivariable ARX models,
the computation of the covariance matrix may take quite a long time. Choosing
Covariance: None will then greatly reduce the computation time.
Resulting Models
The estimated model is inserted into the GUIs Model Board. You can then
examine its various properties and compare them with other models properties
using the Model View plots. More about that in Examining Models on
page 2-28.
To take a look at the model itself, double-click on the models icon (or right-click
or Ctrl-click). The Data/Model Info window that then opens gives you
information about how the model was estimated. You can then also select the
Present button, which will list the model, and its parameters with estimated
standard deviations in the MATLAB command window.
If you need to work further with the model, you can export it by dragging and
dropping it over the To Workspace icon, and then apply any MATLAB and
toolbox commands to it. (See, in particular, the commands ssdata, tfdata, d2c,
and get in the Command Reference chapter.)
2-20
Estimating Models
ARX Models
The Structure
The most used model structure is the simple linear difference equation
y ( t ) + a 1 y ( t 1 ) + + a na y ( t na ) =
b 1 u ( t nk ) + + b nb u ( t nk nb + 1 )
which relates the current output y(t) to a finite number of past outputs y(t-k)
and inputs u(t-k).
The structure is thus entirely defined by the three integers na, nb, and nk. na
is equal to the number of poles and nb 1 is the number of zeros, while nk is the
pure time-delay (the dead-time) in the system. For a system under
sampled-data control, typically nk is equal to 1 if there is no dead-time.
For multi-input systems nb and nk are row vectors, where the i-th element
gives the order/delay associated with the i-th input.
the input orders and delays as a vector. The number of models resulting from
all combinations of orders and delays can however be very large. As an
alternative, you may enter one vector (like nb=1:10) for all inputs and one
vector for all delays. Then only such models are computed that have the same
orders and delays from all inputs.
2-21
Estimation Methods
There are two methods to estimate the coefficients a and b in the ARX model
structure:
Least Squares: Minimizes the sum of squares of the right-hand side minus the
left-hand side of the expression above, with respect to a and b. This is obtained
by selecting ARX as the Method.
Instrumental Variables: Determines a and b so that the error between the
right- and left- hand sides becomes uncorrelated with certain linear
combinations of the inputs. This is obtained by selecting IV in the Method box.
The methods are described in more detail in the Command Reference chapter
under arx and iv4.
Multi-Output Models
For a multi-output ARX structure with ny outputs and nu inputs, the
difference equation above is still valid. The only change is that the coefficients
a are ny-by-ny matrices and the coefficients b are ny-by-nu matrices.
The orders [NA NB NK] define the model structure as follows:
NA: an ny-by-ny matrix whose i-j entry is the order of the polynomial (in the
delay operator) that relates the j-th output to the i-th output
NB: an ny-by-nu matrix whose i-j entry is the order of the polynomial that
2-22
Estimating Models
A ( q )y ( t ) =
[ B i ( q ) Fi ( q ) ] ui ( t nki ) + [ C ( q ) D ( q ) ] e ( t )
i=1
Here ui denotes input #i, and A, Bi, C, D, and Fi, are polynomials in the shift
operator (z or q). (Dont get intimidated by this: It is just a compact way of
writing difference equations; see below.)
The general structure is defined by giving the time-delays nk and the orders of
these polynomials (i.e., the number of poles and zeros of the dynamic model
from u to y, as well as of the disturbance model from e to y).
y ( t ) = [ B ( q ) F ( q ) ]u ( t nk ) + e ( t ) (Output-Error)
y ( t ) = [ B ( q ) F ( q ) ]u ( t nk ) + [ C ( q ) D ( q ) ]e ( t ) (Box-Jenkins)
The shift operator polynomials are just compact ways of writing difference
equations. For example the ARMAX model in longhand would be
y ( t ) + a 1 y ( t 1 ) + + a na y ( t na ) = b 1 u ( t nk ) + +
b nb u ( t nk nb + 1 ) + e ( t ) + c 1 e ( t 1 ) + + c nc e ( t nc )
2-23
Note that A(q) corresponds to poles that are common between the dynamic
model and the disturbance model (useful if disturbances enter the system
close to the input). Likewise F i ( q ) determines the poles that are unique for
the dynamics from input # i, and D(q) the poles that are unique for the
disturbances.
The reason for introducing all these model variants is to provide for flexibility
in the disturbance description and to allow for common or different poles
(dynamics) for the different inputs.
Estimation Method
The coefficients of the polynomials are estimated using a prediction
error/Maximum Likelihood method, by minimizing the size of the error term
e in the expression above. Several options govern the minimization
procedure. These are accessed by activating Iteration Control in the
Parametric Models window, and selecting Options.
The algorithms are further described in Function Reference under armax,
Algorithm Properties, bj, oe, and pem. See also Parametric Model
Estimation and Defining Model Structures.
2-24
Estimating Models
Note These model structures are available only for the scalar output case.
For multi-output models, the state-space structures offer the same flexibility.
Also note that it is not possible to estimate many different structures
simultaneously for the input-output models.
State-Space Models
The Model Structure
The basic state-space model in innovations form can be written
x(t+1) = A x( t) + B u(t) + K e( t)
y (t) = C x (t) + D u(t) + e( t)
The System Identification Toolbox supports two kinds of parametrizations of
state-space models: black-box, free parametrizations, and parametrizations
tailor-made to the application. The latter is discussed below under the heading
User Defined Model Structures. First we will discuss the black-box case.
2-25
Estimation Methods
There are two basic methods for the estimation:
PEM: Is a standard prediction error/maximum likelihood method, based on
quality of the resulting estimates may significantly depend some options called
N4Weight and N4Horizon. These options can be chosen in the Order Editor
window. If N4Horizon is entered with several rows, the models corresponding
to the horizons in each row are examined separately using the Working data.
The best model in terms of prediction (or simulation, if K = 0) performance is
selected. A figure is shown that illustrates the fit as a function of the horizon.
If the N4Horizon box is left empty, a default choice is made.
See n4sid and pem in the Command Reference chapter for more information.
2-26
Estimating Models
To use these structures in conjunction with the GUI, just define the
appropriate structure in the MATLAB command window. Then use the
Structure pop-up menu to select By Initial Model and enter the variable
name of the structure in the edit box Initial Model in the Parametric Models
window and select Estimate.
2-27
Examining Models
Having estimated a model is just a first step. It must now be examined,
compared with other models, and tested with new data sets. This is primarily
done using the six Model View functions, at the bottom of the main ident
window:
Frequency response
Transient response
Zeros and poles
Noise spectrum
Model output
Model residuals
In addition, you can double-click on the models icon to get Text Information
about the model. Finally, you can export the model to the MATLAB workspace
and use any commands for further analysis and model use.
2-28
Examining Models
File
File allows you to copy the current figure to another, standard MATLAB figure
window. This might be useful, e.g., when you intend to print a customized plot.
Other File items cover printing the current plot and closing the plot window.
Options
Options first of all cover actions for setting the axes scaling. This menu item
also gives a number of choices that are specific for the plot window in question,
like a choice between step response or impulse response in the Transient
response window.
A most important option is the possibility to show confidence intervals. Each
estimated model property has some uncertainty. This uncertainty can be
estimated from data. By checking Show confidence intervals, a confidence
region around the nominal curve (model property) will be marked (by
dash-dotted lines). The level of confidence can also be set under this menu item.
Note Confidence intervals are supported for most models and properties,
except models estimated using etfe, and the k-step ahead prediction-property.
For n4sid, the covariance properties are actually not fully known. The
Cramer-Rao lower limit for the covariance matrix is then used instead.
2-29
Style
The style menu gives access to various ways of affecting the plot. You can add
gridlines, turn the zoom on and off, and change the linestyles. The menu also
covers a number of other options, like choice of units and scale for the axis.
Channel
For multivariate systems, you can choose which input-output channel to
examine. The current choice is marked in the figure title.
Help
The Help menu has a number of items, which explain the plot and its options.
Transient Response
Good and simple insight into a models dynamic properties is obtained by
looking at its step response or impulse response. This is the output of the model
when the input is a step or an impulse. These responses are plotted when the
Model View Transient Response is checked.
2-30
Examining Models
2-31
The fit will also be displayed. This is computed as the percentage of the output
variation that is reproduced by the model. So, a model that has a fit of 0% gives
the same mean square error as just setting the model output to be the mean of
the measured output.
If the model is unstable, or has integration or very slow time constants, the
levels of the simulated and the measured output may drift apart, even for a
model that is quite good (at least for control purposes). It is then a good idea to
evaluate the models predicted output rather than the simulated one. With a
prediction horizon of k, the k-step ahead predicted output is then obtained as
follows:
The predicted value y(t) is computed from all available inputs u ( s ) ( s t ) (used
according to the model) and all available outputs up to time t-k, y ( s ) ( s t k ) .
The simulation case, where no past outputs at all are used, thus formally
corresponds to k=. To check if the model has picked up interesting dynamic
properties, it is wise to let the predicted time horizon (kT, T being the sampling
interval) be larger than the important time constants.
Note here that different models use the information in past output data in their
predictors in different ways. This depends on the disturbance model. For
example, so called Output-Error models (obtained by fixing K to zero for
state-space models and setting na=nc=nd=0 for input output models, see the
previous section) do not use past outputs at all. The simulated and the
predicted outputs, for any value of k, thus coincide.
Residual Analysis
In a model
y ( t ) = G ( z )u ( t ) + H ( z )e ( t )
the noise source e(t) represents that part of the output that the model could not
reproduce. It gives the left-overs or, in Latin, the residuals. For a good model,
the residuals should be independent of the input. Otherwise, there would be
more in the output that originates from the input and that the model has not
picked up.
To test this independence, the cross-correlation function between input and
residuals is computed by checking the Model View Model Residuals. It is wise
to also display the confidence region for this function. For an ideal model the
correlation function should lie entirely between the confidence lines for positive
2-32
Examining Models
lags. If, for example, there is a peak outside the confidence region for lag k, this
means that there is something in the output y(t) that originates from u(t-k) and
that has not been properly described by the model. The test is carried out using
the Validation Data. If these were not used to estimate the model, the test is
quite tough. See also Model Structure Selection and Validation.
For a model also to give a correct description of the disturbance properties (i.e.,
the transfer function H), the residuals should be mutually independent. This
test is also carried out by the view Model Residuals, by displaying the
auto-correlation function of the residuals (excluding lag zero, for which this
function by definition is 1). For an ideal model, the correlation function should
be entirely inside the confidence region.
Text Information
By double-clicking (right mouse button or Ctrl-click) on the model icon, a
Data/model Info dialog box opens, which contains some basic information
about the model. It also gives a diary of how the model was created, along with
the notes that originally were associated with the estimation data set. At this
point you can do a number of things:
Present
Selecting the Present button displays details of the model in the MATLAB
command window. The models parameters along with estimated standard
deviations are displayed, as well as some other notes.
Modify
You can simply type in any text you want anywhere in the Diary and Notes
editable text field of the dialog box. You can also change the name of the model
just by editing the text field with the model name. The color, which the model
is associated with in all plots, can also be edited. Enter RGB-values or a color
name (like 'y') in the corresponding box.
LTI Viewer
If you have the Control System Toolbox, you will see an icon To LTI Viewer in
the main window. By dragging and dropping a model onto this icon you will
open the LTI Viewer. This viewer handles an arbitrary amount of models, but
it requires all of them to have the same number of inputs and outputs.
2-33
ss, idss,
ssdata
tf, tfdata
zpk, zpkdata
Note that the commands ss, tf, and zkp transform the model to the Control
System Toolboxs LTI models. Moreover, if you have that toolbox many of its
LTI-commands can be applied directly to the model objects of the Identification
Toolbox. See Connections Between the Control System Toolbox and the
System Identification Toolbox
Also, if you need to prepare specialized plots that are not covered by the Views,
all the System Identification Toolbox commands for computing and extracting
simulations, frequency functions, zeros and poles, etc., are available. See the
Tutorial and Function Reference sections.
2-34
Plot Windows
In the various plot windows the action of the mouse buttons depends on
whether the zoom is activated or not:
Zoom Active: Then the left and middle mouse buttons are associated with the
zoom functions as in the standard MATLAB zoom. Left button zooms in and the
2-35
middle one zooms out. In addition, you can draw rectangles with the left
button, to define the area to be zoomed. Double-clicking restores the original
plot. The right mouse button is associated with special GUI actions that depend
on the window. In the View plots, the right mouse button is used to identify the
curves. Point and click on a curve, and a box will display the name of the
model/data set that the curve is associated with, and also the current
coordinate values for the curve. In the Model Selection plots the right mouse
button is used to inspect and select the various models. In the Prefilter and
Data Range plots, rectangles are drawn with this mouse button down, to
define the selected range.
Zoom not active: The special GUI functions just mentioned are obtained by
any mouse button.
The zoom is activated and deactivated under the menu item Style. The default
setting differs between the plots. Dont activate the zoom from the command
line! That will destroy the special GUI functions. (If you happen to do so
anyway, quit the window and open it again.)
Troubleshooting in Plots
The function Auto-range, which is found under the menu item Options, sets
automatic scales to the plots. It is also a good function to invoke when you think
that you have lost control over the curves in the plot. (This may happen, for
example, if you have zoomed in a portion of a plot and then change the data of
the plot).
If the view plots dont respond the way you expect them to, you can always
quit the window and open it again. By quit here we mean using the
underlying window systems own quitting mechanism, which is called different
things in the different platforms. The normal way to close a window is to use
the Close function under the menu item File, or to uncheck the corresponding
check box.
2-36
GUI to start with your current window layout and current plot options, select
menu item
Options > Save preferences
in the main ident window. This saves the information in a file idprefs.mat.
This file also stores information about the four most recent sessions with ident.
This allows the session File menu to be correctly initialized. The session
information is automatically stored upon exit. The layout and preference
information is only saved when the indicated option is selected.
The file idprefs.mat is created the first time you close the GUI. It is by default
stored in the same directory as your startup.m file. If this default does not
work, you are prompted for a directory where to store the file. This can be
ignored, but then session and preference information cannot be saved.
To change or select a directory for idprefs.mat, use the command midprefs.
See the Command Reference chapter for details.
To change model colors and default options to your own customized choice,
make a copy of the M-file idlayout.m to your own directory (which should be
before the basic ident directory in the MATLABPATH), and edit it according to its
instructions.
Customized Plots
If you need to prepare hardcopies of your plots with specialized texts, titles and
so on, make a copy of the figure first, using Copy Figure under the File menu
item. This produces a copy of the current figure in a standard MATLAB figure
format.
For plots that are not covered by the View windows, (e.g., Nyquist plots), you
have to export the model to the MATLAB workspace and construct the plots
there.
2-37
2-38
3
Tutorial
Tutorial
Overview
This chapter has three purposes:
It gives an overview of System Identification theory, the basic models and
disturbance descriptions used, and the character of the basic algorithms. It
also provides some practical advice for a number of issues that are essential
for a successful application.
It describes the commands and objects of the System Identification Toolbox,
their syntax and use. If you primarily use the graphical user interface (GUI),
you will not have to bother about these aspects.
It also describes the commands that are not reached from the GUI; i.e.,
simulation, the recursive algorithms and more advanced model structure
definitions.
3-2
3-3
Tutorial
Layer 3: Model Structure Selection. The third layer of the toolbox contains some
useful techniques to select orders and delays.
arxstruc, selstruc
Layer 4: Structured Models and Further Model Conversions. The fourth layer contains
transformations between continuous and discrete time, and functions for
estimating completely general model structures for linear systems. The
commands are
c2d, d2c, idss, idgrey, pe, predict
ss, tf, zp, frd (to be used with the Control System Toolbox)
3-4
It contains the input vector u2, the output vector y2. First form the data object.
dry = iddata(y2,u2,0.08);
To get an overview of all the information contained in the iddata object dry,
type
get(dry)
3-5
Tutorial
Let us first estimate the impulse response of the system by correlation analysis
to get some idea of time constants and the like.
impulse(ze,'sd',3)
When the calculations are finished, a display of the basic information about m1
is shown. Anytime m1 is typed, this display is shown. Typing present(m1) will
give some more information about the model, including uncertainties.
To retrieve the properties of this model we could, e.g., find the A matrix of the
state space representation by
A = m1.a
m1 is a model object, and
get(m1)
3-6
How good is this model? One way to find out is to simulate it and compare the
model output with measured output. We then select a portion of the original
data that was not used to build the model, e.g., from sample 800 to 900.
zv = dry(800:900);
zv = detrend(zv);
compare(zv,m1);
We can also compare the step response of the model with one that is directly
computed from data (ze) in a nonparametric way.
step(m1,ze)
3-7
Tutorial
3-8
e
y
t = 1, 2, , N
t = 1, 2, , N
Assuming the signals are related by a linear system, the relationship can be
written
y ( t ) = G ( q )u ( t ) + v ( t )
(3-1)
G ( q )u ( t ) =
g ( k )u ( t k )
(3-2)
k=1
and
3-9
Tutorial
g ( k )q
G(q) =
q u(t) = u(t 1)
(3-3)
k=1
The numbers { g ( k ) } are called the impulse response of the system. Clearly,
g ( k ) is the output of the system at time k if the input is a single (im)pulse at
time zero. The function G ( q ) is called the transfer function of the system. This
i
function evaluated on the unit circle ( q = e ) gives the frequency function
i
G( e )
(3-4)
(3-5)
which is defined by
v ( ) =
R v ( )e
(3-6)
(3-7)
(3-8)
v ( ) = H ( e )
(3-9)
Equations (Equation 3-1) and (Equation 3-8) together give a time domain
description of the system
y ( t ) = G ( q )u ( t ) + H ( q )e ( t )
(3-10)
where G is the transfer function of the system. Equations (Equation 3-4) and
(Equation 3-5) constitute a frequency domain description.
3-10
G ( e );
v ( )
(3-11)
The impulse response (Equation 3-3) and the frequency domain description
(Equation 3-11) are called nonparametric model descriptions since they are not
defined in terms of a finite number of parameters. The basic description
(Equation 3-10) also applies to the multivariable case; i.e., to systems with
several (say nu) input signals and several (say ny) output signals. In that
case G ( q ) is an ny-by-nu matrix while H ( q ) and v ( ) are ny-by-ny matrices.
nk
B(q)
------------ ;
A(q)
1
H ( q ) = -----------A(q)
(3-12)
1
B ( q ) = b 1 + b2 q
+ + a na q
na
+ + b nb q
nb + 1
(3-13)
Here, the numbers na and nb are the orders of the respective polynomials. The
number nk is the number of delays from input to output. The model is usually
written
A ( q )y ( t ) = B ( q )u ( t nk ) + e ( t )
(3-14)
or explicitly
y ( t ) + a 1 y ( t 1 ) + + a na y ( t na ) =
b 1 u ( t nk ) + b 2 u ( t nk 1 ) + + b nb u ( t nk nb + 1 ) + e ( t )
(3-15)
Note that (Equation 3-14) - (Equation 3-15) apply also to the multivariable
case, with ny output channels and nu input channels. Then A ( q ) and the
coefficients a i become ny-by-ny matrices, B ( q ) and the coefficients b i become
ny-by-nu matrices.
3-11
Tutorial
Another very common, and more general, model structure is the ARMAX
structure
A ( q )y ( t ) = B ( q )u ( t nk ) + C ( q )e ( t )
(3-16)
+ + c nc q
nc
(3-17)
with
F ( q ) = 1 + f1 q
+ + fnf q
nf
(3-18)
with
D(q ) = 1 + d1q
+ + d nd q
nd
All these models are special cases of the general parametric model structure.
B(q )
C(q)
A ( q )y ( t ) = ------------ u ( t nk ) + ------------- e ( t )
F(q)
D(q)
(3-19)
3-12
The same type of models can be defined for systems with an arbitrary number
of inputs. They have the form
B1 ( q )
B nu ( q )
C(q)
A ( q )y ( t ) = --------------- u 1 ( t nk 1 ) + ...+ ------------------- u nu ( t nk nu ) + ------------- e ( t )
D(q)
F1 ( q )
F nu ( q )
(3-20)
(3-21)
Here the relationship between the input u ( t ) and the output y ( t ) is defined
via the nx-dimensional state vector x ( t ) . In transfer function form
(Equation 3-21) corresponds to (Equation 3-1) with
1
G ( q ) = C ( qI nx A ) B + D
(3-22)
(3-23)
H ( q ) = C ( qI nx A ) K + I ny
(3-24)
3-13
Tutorial
x ( t + 1 ) = Ax ( t ) + Bu ( t ) + w ( t )
(3-25)
y ( t ) = Cx ( t ) + Du ( t ) + e ( t )
(3-26)
y ( t ) = Hx ( t ) + Du ( t ) + v ( t )
Here, x means the time derivative of x . If the input is piece-wise constant over
time intervals kT t < ( k + 1 )T , then the relationship between u [ k ] = u ( kT )
and y [ k ] = y ( kT ) can be exactly expressed by (Equation 3-21) by taking
T
A = e
FT
B =
G d;
C = H
(3-27)
(3-28)
K =
e
0
3-14
d;
K
(3-29)
(3-30)
and estimate the g-coefficients by the linear least squares method. In fact, to
check if there are non-causal effects from input to output, e.g., due to feedback
from y in the generation of u (closed loop data), g for negative lags can also be
estimated.
y ( t ) = g ( m )u ( t + m ) + + g ( 1 )u ( t + 1 ) + g ( 0 )u ( t ) +
g ( 1 )u ( t 1 ) + + g ( n )u ( t n )
(3-31)
y ( ) = G ( e ) u ( ) + v ( )
i
(3-32)
yu ( ) = G ( e ) u ( )
By estimating the various spectra involved, the frequency function and the
disturbance spectrum can be estimated as follows.
3-15
Tutorial
1
( ) = ---R
yu
N
y ( t + )u ( t )
(3-33)
t=1
and analog expressions for the others. Then, form estimates of the
corresponding spectra
M
y ( ) =
i
Ry ( )W M ( )e
(3-34)
= M
and analogously for u and yu . Here W M ( ) is the so-called lag window and
M is the width of the lag window. The estimates are then formed as
yu ( )
( e i ) = ------------------;
G
N
u ( )
yu ( )
v ( ) = y ( ) -----------------------
u( )
(3-35)
e ( t ) = H ( q ) [ y ( t ) G ( q )u ( t ) ]
(3-36)
These errors are, for given data y and u, functions of G and H. These in turn
are parametrized by the polynomials in (Equation 3-14)-(Equation 3-19) or by
entries in the state-space matrices defined in (Equation 3-26)(Equation 3-29).
The most common parametric identification method is to determine estimates
of G and H by minimizing
N
V N ( G, H ) =
e
t=1
that is
3-16
(t)
(3-37)
,H
] = argmin
[G
N
N
(t)
(3-38)
t=1
(3-39)
3-17
Tutorial
observed data sequences. See Sections 7.3 and 10.6 in Ljung (1999). For more
details, see the references under n4sid in the Command Reference chapter.
3-18
Data Representation
The observed output and input signals, y ( t ) and u ( t ) , are represented as
column vectors y and u. Row k corresponds to sample number k. For
multivariable systems, each input (output) component is represented as a
column vector, so that u becomes an N-by-nu matrix (N = number of sampled
observations, nu = number of input channels). The output-input data is
collectively represented in the iddata format. This is the basic object for
dealing with signals in the toolbox. It is used by most of the commands. It is
created by
Data = iddata(y,u,Ts)
More details about the iddata object is given at the end of this section.
3-19
Tutorial
Correlation Analysis
The correlation analysis procedure described in Estimating Impulse
Responses is implemented in the function impulse.
impulse(Data)
This function plots the estimated impulse response. Adding an argument 'sd'
as in
impulse(Data,'sd',3)
it also marks a confidence region corresponding to (in this case) three standard
deviations. The result can be stored and replotted.
ir = impulse(Data)
impulse(ir,'sd',3)
An alternative is the command step that plots the step response, calculated
from the impulse estimate.
step(Data)
Spectral Analysis
The function spa performs spectral analysis according to the procedure in
(Equation 3-35)(Equation 3-37).
g = spa(Data)
Here Data contains the output-input data in the iddata object as above. g is
returned as an idfrd (Identified frequency domain) model object, that contains
the estimated
frequency function G N and the estimated disturbance
3-20
performs the spectral analysis, and plots first G and then v . bode gives
logarithmic amplitude and frequency scales (in rad/sec) and linear phase scale,
while ffplot gives linear frequency scales (in Hz). The uncertainty of the
estimates is displayed by adding the argument 'sd' as in
bode(g,'sd',3)
Similarly
nyquist(g)
gives a Nyquist plot of the frequency function, i.e. a plot of the real part versus
the imaginary part of G
If Data = y is a time series, that is Data has no input channel, spa returns an
estimate of the spectrum of that signal:
g= spa(y)
ffplot(g)
The rule is that as M increases, the estimated frequency functions show sharper
details, but are also more affected by random disturbances. A typical sequence
of commands that test different window sizes is
g10 = spa(Data,10)
g25 = spa(Data,25)
g50 = spa(Data,50)
bode(g10, g25, g50)
3-21
Tutorial
g = etfe(Data)
This can also be interpreted as the spectral analysis estimate for a window size
that is equal to the data length. For time series, etfe gives the periodogram as
a spectral estimate. The function also allows some smoothing of the crude
estimate; it can be a good alternative for signals and systems with sharp
resonances. See Function Reference for more information.
(two inputs and one output is this example) and these names will then follow
the object and appear in all plots. The names will also be inherited by models
that are estimated from the data.
Similarly, channel units can be specified using the properties OutputUnit and
InputUnit. These units, when specified, will be used in plots.
The timepoints associated with the data samples are determined by the
sampling interval Ts and the time of the first sample, Tstart.
Data.Tstart = 24
The actual time point values are given by the property, SamplingInstants, as
in
plot(Data.sa,Data.u)
for a plot of the input with correct time points. Autofill is used for all properties,
and they are case insensitive. For easy writing 'u' is synonymous to 'Input'
and 'y' to 'Output' when referring to the properties.
Manipulating Channels
An easy way to set and retrieve channel properties is to use subscripting. The
subscripts are defined as
Data(samples,outputs,inputs),
3-22
so Dat(:,3,:) is the data object obtained from Dat by keeping all input
channels, but only output channel 3. (Trailing ':'s can be omitted so
Dat(:,3,:)= Dat(:,3).)
The channels can also be retrieved by their names, so that
Dat(:,{'speed','flow'},[ ])
is the data object where the indicated output channels have been selected and
no input channels are selected.
Moreover
Dat1(101:200,[3 4],[1 3]) = Dat2(1001:1100,[1 2],[6 7])
will change samples 101 to 200 of output channels 3 and 4 and input channels
1 and 3 in the iddata object Dat1 to the indicated values from iddata object
Dat2. The names and units of these channels will then also be changed
accordingly.
To add new channels, use horizontal concatenation of iddata objects
Dat =[Dat1, Dat2];
Nonequal Sampling
The property SamplingInstants gives the sampling instants of the data points.
It can always be retrieved by get(Dat,'SamplingInstants') (or Dat.s) and is
then computed from Dat.Ts and Dat.Tstart. SamplingInstants can also be
set to an arbitrary vector of the same length as the data, so that nonequal
sampling can be handled. Ts is then automatically set to [ ]. Most of the
estimation routines, though, do not handle unequally sampled data.
Multiple Experiments
The iddata object can also store data from separate experiments. The property
ExperimentName is used to separate the experiments. The number of data as
well as the sampling properties can vary from experiment to experiment, but
the input and output channels must be the same. (Use NaN to fill unmeasured
3-23
Tutorial
channels in certain experiments.) The data records will be cell arrays, where
the cells contain data from each experiment.
Multiple experiments can be defined directly by letting the 'y' and 'u'
properties as well as 'Ts' and 'Tstart' be cell arrays.
It is normally easier to create multi-experiment data by merging experiments
as in
Dat = merge(Dat1,Dat2)
adds the data in Dat2 as a new experiment with name 'Run4'. See iddemo # 8
for an illustration of how multiple experiments can be used.
iddata Properties
Type get(Dat) or see iddata in the Command Reference for a complete list
of iddata properties.
Subreferencing
The samples, outputs and input channels can be referenced according to
Data(samples,outputs,inputs)
3-24
Use colon (:) to denote all samples/channels and the empty matrix ([ ]) to denote
no samples/channels. The channels can be referenced by number or by name.
For several names a cell array must be used.
Dat2 = Dat(:,'y3',{'u1','u4'})
Dat2 = Dat(:,3,[1 4])
will select the samples with time marks between 1.27 and 9.3.
Subreferencing with curly brackets refers to the experiment.
Data{Experiment}(samples,outputs,inputs)
Adding Channels
Dat = [Dat1,Dat2,...,DatN]
creates an iddata object Dat, consisting of the input and output channels in
Dat1,... DatN. Default channel names ('u1', 'u2', 'y1' , 'y2' etc) will
be changed so that overlaps in names are avoided, and the new channels will
be added.
If Datk contains channels with user-specified names that are already present
in the channels of Datj, j<k, these new channels will be ignored.
Adding Samples
Dat = [Dat1;Dat2;... ;DatN]
creates an iddata object Dat whose signals are obtained by stacking those of
Datk on top of each other. That is.
Dat.y = [Dat1.y;Dat2.y; ... DatN.y]
and similarly for the inputs. The Datk objects must all have the same number
of channels and experiments.
3-25
Tutorial
The argument Data is an iddata object that contains the output and input data
sequences, while modstruc specifies the particular structure of the model to be
estimated. The resulting estimated model is contained in m. It is a model object
that stores various information. The model objects will be described in
Defining Model Structures, but for most use of the toolbox, you do not have to
consider the details of these objects. Just typing the model name
m
gives a complete list of the model s properties. The property values can be
easily retrieved just by dot-referencing; for example,
m.par
3-26
In the following, Data denotes an iddata object that contains the input output
data as described in the previous section. It may also just contain an output
signal, i.e., a time series.
ARX Models
To estimate the parameters a i and b i of the ARX model (Equation 3-14), use
the function arx.
m
= arx(Data,[na nb nk])
Here na, nb, and nk are the corresponding orders and delays in (Equation 3-15)
on page 3-11 that define the exact model structure. The function arx
implements the least squares estimation method, using QR-factorization for
overdetermined linear equations.
An alternative is to use the Instrumental Variable (IV) method described in
connection with (Equation 3-39). This is obtained with
m = iv4(Data,[na nb nk])
AR Models
For a single output signal y ( t ) the counterpart of the ARX model is the AR
model.
A ( q )y ( t ) = e ( t )
(3-40)
3-27
Tutorial
but for scalar signals more options are offered by the command
m = ar(y,na)
which has an option that allows you to choose the algorithm from a group of
several popular techniques for computing the least squares AR model. Among
these are Burgs method, a geometric lattice method, the Yule-Walker
approach, and a modified covariance method. (See Function Reference for
details.) The counterpart of the iv4 command is
m = ivar(y,na)
The nonzero orders of the model can also be defined as property name/property
value pairs as in
m = pem(Data,'na',na,'nb',nb,'nc',nc,'nk',nk)
3-28
B1 ( q )
B nu ( q )
C(q)
A ( q )y ( t ) = --------------- u 1 ( t nk 1 ) + + ------------------- u nu ( t nk nu ) + ------------- e ( t ) (3-41)
D(q)
F1 ( q )
F nu ( q )
where nb, nf, and nk are row vectors of the same lengths as the number of input
channels containing each of the orders and delays
nb = [nb1 ...
nf = [nf1 ...
nk = [nk1 ...
nbnu];
nfnu];
nknu];
While the search is typically initialized using the built-in start-up procedure
giving just orders and delays (as described above), the ability to force a specific
initial condition is useful in several contexts. Some examples are mentioned in
Initial Parameter Values.
Information about how the minimization progresses can be supplied to the
MATLAB command window by the property trace. See the list at the end of
this section.
3-29
Tutorial
State-Space Models
Black-Box, Discrete Time Parametrizations
Suppose first that there is no particular knowledge about the internal
structure of the discrete-time state-space model (Equation 3-15). Any linear
model is sought. A simple approach then is to use
m = pem(Data)
To get a plot, from which the order can be determined among a list of orders
nn = [n1,n2,...,nN], use
m = pem(Data,'nx',nn)
All these black-box models are initialized by the subspace method n4sid. To
obtain the estimate from this routine, use
m = n4sid(Data,n)
When the systems are multi-output, the following criterion is used for the
minimization,
N
det
e ( t )e
(t)
(3-42)
t=1
which is the maximum likelihood criterion for Gaussian noise with unknown
covariance matrix.
3-30
Optional Variables
The estimation functions accept a list of property name/property value pairs
that may affect both the model structure and the estimation algorithm. For
complete lists of these properties, see algorithm properties, idarx, idmodel,
idpoly, idss, and idgrey in the Function Reference. Some of them are listed
here. Note that any property, as well as values that are strings, can be entered
as any unambiguous, case-insensitive abbreviation, as is
m = pem(Data,mi,'fo','si').
Note The algorithm properties, like all other model properties, will be
inherited by the resulting model m. If you continue the estimation using m as
initial model, all previously set algorithm features will thus apply, unless you
specify otherwise.
3-31
Tutorial
Prediction: This is the default and means that the model is determined by
minimizing the prediction errors. It corresponds to a frequency weighting
that is given by the input spectrum times the inverse noise model. Typically,
this favors a good fit at high frequencies. From a statistical variance point of
view, this is the optimal weighting, but then the approximation aspects
(bias) of the fit are neglected.
Simulation: This means that frequency weighting of the transfer function fit
is given by the input spectrum. Frequency ranges where the input has
considerable power will thus be better described by the model. In other
words, the model approximation is such that the model will produce as good
simulations as possible, when applied to inputs with the same spectra as
used for the estimation. For models that have no disturbance model (A=C=D
for idpoly models and K=0 for idss models) there is no difference between
the Simulation and Prediction values. For models with a disturbance
description, this is estimated by a prediction error method, keeping the
estimated transfer function from input to output fixed. The resulting model
is guaranteed to be stable.
Stability: The algorithm is modified so that a stable model is guaranteed, but
the weighing still correspond to prediction.
Any SISO linear filter. Then the transfer function from input to output is
determined by a frequency fit with this filter times the input spectrum as
weighting function. The noise model is determined by a prediction error
method, keeping the transfer function estimate fixed. To obtain a good model
fit over a special frequency range, the filter should thus be chosen with a
passband over this range. For a model with no disturbance model, the result
3-32
is the same as first applying prefiltering to data using idfilt. The filter can
be specified in a few different ways:
- as any single-input-single-output idmodel
- as a ss, tf, or zpk model from the Control System Toolbox
- as {A,B,C,D} with the state-space matrices for the filter
- as {numerator, denominator} with the transfer function
numerator/denominator of the filter
MaxSize: No matrix with more than MaxSize elements is formed by the
algorithm. Instead, for loops will be used. MaxSize will thus decide the
memory/speed trade-off, and can prevent slow use of virtual memory. MaxSize
can be any positive integer, but it is of course required that the input-output
data themselves contain less than MaxSize elements. The default value of
MaxSize is Auto, which means that the value is determined in the M-file
idmsize. This file may be edited by the user to optimize speed on a particular
computer. See also Memory/Speed Trade-Offs.
FixedParameter: A list of parameters that will be kept fixed to the
nominal/initial values and not estimated. This is a vector of integers containing
the indices of the fixed parameters, or a cell array of parameter names. If
names are used, wildcard entries apply, which may be convenient if you have
groups of parameters in your model. See the reference page of Algorithm
properties.
3-33
Tutorial
3-34
quadratic to one that gives linear weight to large errors. Errors larger than
LimitError times the estimated standard deviation will carry a linear weight
in the criteria. The default value of LimitError is 1.6. LimitError =0 disables
gn: The Gauss-Newton direction (inverse of the Hessian times the gradient
direction) If no improvement is found along this direction, the gradient
direction is also tried out.
gns: A regularized version of the Gauss-Newton direction. Eigenvalues less
than pinvtol of the Hessian are neglected, and the Gauss-Newton direction
is computed in the remaining subspace. (pinvtol is part of the 'advanced'
field: See Algorithm Properties in the Command Reference.
lm: The Levenberg-Marquard method is used. This means that the next
parameter value is -pinv(H+d*I)*grad from the previous one, where H is the
Hessian, I is the identity matrix, grad is the gradient. d is a number that is
increased until a lower value of the criterion is found.
Auto: A choice between the above is made in the algorithm. This is the
default choice.
One property of the returned model is EstimationInfo. That is a structure that
contains useful information about the estimation process. See EstimationInfo
in Command Reference for a list of fields.
3-35
Tutorial
3-36
to see what the assignable values are. See get and/or set in the Command
Reference chapter. Each property value can easily also be retrieved by
subreferencing as in
m.A
and set as in
m.b(3) = 27
See Function Reference for complete property lists. Here only examples are
given. Note that it is sufficient to use any case insensitive, unambiguous
abbreviation of the property names.
3-37
Tutorial
(3-43)
is defined by the five polynomials A(q), B(q), C(q), D(q), and F(q). These are
represented in the standard MATLAB format for polynomials. Polynomial
coefficients are stored as row vectors ordered by descending powers. For
example, the polynomial
A ( q ) = 1 + a1 q
+ a2 q
+ + an q
is represented as
A = [1 a1 a2 ...
an]
Delays in the system are indicated by leading zeros in the B ( q ) polynomial. For
example, the ARX model
y ( t ) 1.5y ( t 1 ) + 0.7y ( t 2 ) = 2.5u ( t 2 ) + 0.9u ( t 3 )
(3-44)
In the multi-input case (Equation 3-41) B and F are matrices, whose row
number k corresponds to the k-th input. The command idpoly can also be used
to define time-continuous systems. See Function Reference for details.
3-38
When m is defined, the polynomials and their orders can be easily retrieved and
changed, as in
m.a % for the A-polynomial
roots(m.a)
m.a(3)=0.95
(3-45)
Here A(q) is an ny-by-ny matrix whose entries are polynomials in the delay
operator q-1. You can represent it as
A ( q ) = I ny + A 1 q
+ + A na q
na
(3-46)
a 21 ( q ) a 22 ( q ) a 2ny ( q )
(3-47)
a ny1 ( q ) a ny 2 ( q ) a nyny ( q )
a kj ( q ) = kj + a kj q
na k j na kj
+ + a kj q
:
(3-48)
This polynomial describes how old values of output number j affect output
number k. Here kj is the Kronecker-delta; it equals 1 when k = j , otherwise,
it is 0. Similarly, B ( q ) is an ny-by-nu matrix
B( q ) = B0 + B1 q
+ B nb q
nb
(3-49)
or
3-39
Tutorial
b 11 ( q ) b 12 ( q ) b 1nu ( q )
B(q ) =
b 21 ( q ) b 22 ( q ) b 2nu ( q )
(3-50)
with
1
b kj ( q ) = b kj q
nk kj
nb kj nk kj nb kj + 1
+ + b kj q
The delay from input number j to output number k is nk kj . To link with the
structure definition in terms of na, nb, nk in the arx and iv4 commands, note
that na is a matrix whose kj-element is nakj , while the kj-elements of nb and
nk are nb kj and nk kj respectively.
The idarx representation of the model (Equation 3-45) can be created by
m = idarx(A,B)
Note that A(:,:,1) is always the identity matrix, and that leading zero
coefficients in B matrices define the delays.
Consider the following system with two outputs and three inputs.
y 1 ( t ) 1.5y 1 ( t 1 ) + 0.4y 2 ( t 1 ) + 0.7y 1 ( t 2 ) =
0.2u 1 ( t 4 ) + 0.3u 1 ( t 5 ) + 0.4u 2 ( t ) 0.1u 2 ( t 1 ) + 0.15u 2 ( t 2 ) + e 1 ( t )
y 2 ( t ) 0.2y 1 ( t 1 ) 0.7y 2 ( t 2 ) + 0.01y 1 ( t 2 ) =
u 1 ( t ) + 2u 2 ( t 4 ) + 3u 3 ( t 1 ) + 4u 3 ( t 2 ) + e 2 ( t )
which in matrix notation can be written as
3-40
3-41
Tutorial
(a)
(b)
x ( 0 ) = x0
(3-51)
(c)
(3-52)
x ( 0 ) = x0
(See Ljung (1999), page 93.) It is often easier to define a parameterized
state-space model in continuous time because physical laws are most often
described in terms of differential equations. The matrices F, G, H, and D
contain elements with physical significance (for example, material constants).
The numerical values of these may or may not be known. To estimate unknown
parameters based on sampled data (assuming T is the sampling interval) first
transform (Equation 3-52) to (Equation 3-51) using the formulas
(Equation 3-27). The value of the Kalman gain matrix K in (Equation 3-51) or
~
K in (Equation 3-52) depends on the assumed character of the additive noises
w ( t ) and e ( t ) in (Equation 3-25), and its continuous-time counterpart.
~
Disregard that link and view K in (Equation 3-51) (or K in (Equation 3-52)) as
the basic tool to model the disturbance properties. This gives the directly
parametrized innovations form. (See Ljung (1999) page 99.) If the internal
noise structure is important, you could use user-defined greybox structures
(the idgrey object) as in the example.
The discrete time model (Equation 3-51) can be put into the idss model by
m = idss(A,B,C,D,K,X0,'Ts',T)
3-42
any elements in the matrices are freely adjustable by the estimation routines.
The parameters will be adjusted to data by
me = pem(Data,m)
The iterative search for the best fit is then initialized in the nominal matrices
A, B, C, D, K, X0. Note that the command me = pem(Data,4), which just defines
the model order, first estimates (using n4sid) a starting model m, from which
the search is initialized.
In this free parameterization, you can decide how to deal with the disturbance
model matrix K. Letting
m.DisturbanceModel = 'None'
(no delays) means that all elements of the D-matrix should be estimated, while
m.nk = [1,1,..,1]
(rather than 'Free'). This is still a black-box model, since the canonical form
covers all models of a certain order. The structure modifications can all be
combined at the estimation call
me = pem(Data,m,'sspar','can','dist','none','ini','z')
3-43
Tutorial
sets the structure matrix for A, called As, to a diagonal matrix, where the
diagonal elements are freely adjustable. Defining
m.A = [2 0; 0 3]
x( t + 1 ) =
1 1
0 1
x( t) +
2
3
u( t) +
4
5
e(t)
y( t) = 1 0 x( t) + e( t)
x(0) = 0
0
with five unknown parameters i , i=1,...,5. Suppose the nominal/initial values
of these parameters are -1, 2, 3, 4 and 5. This structure is then defined by
3-44
The definition thus follows in two steps. First the nominal model is defined.
Then the structure (known and unknown parameter values) is defined by the
structure matrices, As, Bs, etc.
Example 3.2: A Continuous-Time Model Structure. Consider the following model
structure
0 1
0
x ( t ) =
x(t) +
u(t)
0 1
2
y( t) = 1 0 x( t) + e(t )
01
x( 0) =
3
0
3-45
Tutorial
m.x0s = [NaN;0]
m.noisevar = [0.01 0; 0 0.1]
The structure m can now be used to estimate the unknown parameters i from
observed data
Data = iddata([y1 y2], u, 0.1)
by
model = pem(Data,m)
The iterative search for minimum is then initialized at the parameters in the
nominal model m. The continuous time model is automatically sampled to agree
with the sampling interval of the data. The structure can also be used to
simulate the system above with sampling interval T=0.1 for input u and noise
realization e.
e = randn(300,2)
u = idinput(300);
simdat = iddata([],[u e],'Ts',0.1);
y = sim(m,simdat) % The continuous system will automatically be
% sampled using Ts =0.1
The nominal parameter values are then used, and the noise sequence is scaled
according to the matrix m.noisevar.
When estimating models, you can try a number of neighboring structures,
such as what happens if I fix this parameter to a certain value or what
happens if I let loose these parameters. This is easily handled by the structure
matrices As, Bs, etc. For example, to free the parameter x2(0) (perhaps the
motor wasnt at rest after all), you can use
model = pem(Data,m,'x0s',[NaN;NaN])
3-46
greybox modeling and write an M-file that describes the structure. The
format is
[A,B,C,D,K,x0] = mymfile(par,T,aux);
where mymfile is the user-defined name for the M-file, par contains the
parameters as a column vector, T is the sampling interval, and aux contains
auxiliary variables. The latter variables are used to house options, so that some
different cases can be tried out without having to edit the M-file. The matrices
A, B, C, D, K, and x0 refer either to the continuous time description
(Equation 3-52) or to the discrete-time description (Equation 3-51). When a
continuous time description is fitted to sampled data, the estimation routines
perform the necessary sampling of the model. To obtain the same structure as
in the Example 3.2, you can do the following.
function [A,B,C,D,K,x0] = mymfile(par,T,aux)
A = [0 1; 0 par(1)];
B = [0;par(2)];
C = eye(2);
D = zeros(2,2);
K = zeros(2,1);
x0 =[par(3);0];
Once the M-file has been written, the idgrey model m is defined by the command
m = idgrey('mymfile',par,'c',aux);
where par contains the nominal (initial) values of the corresponding entries in
the structure. 'c' signals that the underlying parametrization is continuous
time. aux contains the values of the auxiliary parameters. Note that T and aux
must be given as input arguments, even if they are not used by the code.
From here on, estimate models and evaluate results as for any other model
structure. Some further examples of user-defined model structures are given
below.
3-47
Tutorial
3-48
for a 10th order approximation of a heat rod one meter in length, with an initial
temperature of 22 degrees. The initial estimate of the heat conductivity is 0.27,
and of the heat transfer coefficient is 1.
The model parameters are estimated by
me = pem(Data,m)
If you would like to try a finer grid, that is, take Ngrid larger, you can do this
easily with
me = pem(Data,m,'Filearg',[20,1,22])
Example 3.4: Parametrized Disturbance Models. Consider a discrete-time model
x ( t + 1 ) = Ax ( t ) + Bu ( t ) + w ( t )
y ( t ) = Cx ( t ) + e ( t )
where w and e are independent white noises with covariance matrices R1 and
R2, respectively. Suppose that you know the variance of the measurement
noise R2, and that only the first component of w ( t ) is nonzero. This can be
handled by the following M-file.
function [A,B,C,D,K,x0] = mynoise(par,T,aux)
R2 = aux(1); % The assumed known measurement noise variance
A = [par(1) par(2);1 0];
B = [1;0];
C = [par(3) par(4)];
D = 0;
R1 = [par(5) 0;0 0];
K = Adlqe(A,eye(2),C,R1,R2); % from the Control System Toolbox
x0 = [0;0];
3-49
Tutorial
In the search for the minimum, the gradient of the prediction errors e ( t ) is
computed by numerical differentiation. The step-size is determined by the
M-file nuderst. In its default version the step-size is simply 10 4 times the
absolute value of the parameter in question (or the number 10 7 if this is
larger). When the model structure contains parameters with different orders of
magnitude, try to scale the variables so that the parameters are all roughly the
same magnitude. You may need to edit the M-file nuderst to address the
problem of suitable step sizes for numerical differentiation.
3-50
Examining Models
Examining Models
Once you have estimated a model, you need to investigate its properties. You
have to simulate it, test its predictions, and compute its poles and zeros and so
on. You thus have to transform the model to various ways of representing and
presenting it. This section deals with how this is done. The following topics will
be covered:
Parametric models: basic use, accessing properties, simulation and
prediction. Also manipulating channels, in particular the noise input
channels.
Frequency domain models
Graphing model properties
Transformations to other representations
Transformations between continuous and discrete time
Basic Use
If you just estimate models from data, the model objects should be transparent.
All parametric estimation routines return idmodel results.
m = arx(Data,[2 2 1])
The model m contains all relevant information. Just typing m will give a brief
account of the model. present(m) also gives information about the
uncertainties of the estimated parameters. get(m) gives a complete list of
model properties.
Most of the interesting properties can be directly accessed by subreferencing.
m.a
m.da
See the property list obtained by get(m), as well as the property lists of
idgrey, idarx, idpoly, and idss in the Command Reference for more details
on this.
3-51
Tutorial
and by idfrd and freqresp, the frequency response of the model can be
computed.
The number of input channels must either be equal to the number of measured
channels in m, in which case a noise free simulation is obtained, or equal to the
sum of the number of input and output channels in m. In the latter case the last
input signals (v) are interpreted as white noise. They are then scaled by the
NoiseVariance matrix of m and added to the output via the disturbance model
y = Gu + He
e = Lv
T
The output is returned as an iddata object with just output channels. Here is
a typical string of commands.
A = [1 -1.5 0.7];
B = [0 1 0.5];
m0 = idpoly(A,B,[1 -1 0.2]);
u = iddata([],idinput(400,'rbs',[0 0.3]));
v= iddata([],randn(400,1));
3-52
Examining Models
y = sim(m0, [u v]);
plot(y)
The inverse model (Equation 3-38), which computes the prediction errors
from given input-output data, is simulated with
e = pe(m,[y u])
To compute the k-step ahead prediction of the output signal based on a model
m, the procedure is as follows.
yhat = predict(m,[y u],k)
say, monthly sales figures. A model is estimated based on the first half, and
then its ability to predict half a year ahead is checked out on the second half of
the observations.
plot(y)
y1 = y(1:48), y2 = y(49:96)
m4 = ar(y1,4)
yhat = predict(m4,y2,6)
plot(y2,yhat)
The command compare is useful for any comparisons involving sim and
predict.
3-53
Tutorial
Use colon (:) to denote all channels and the empty matrix ([ ]) to denote no
channels. The channels can be referenced by number or by name. For several
names, a cell array must be used.
m3 = m('position',{'power','speed'})
or
m3 = m(3,[1 4])
Thus m3 is the model obtained from m by considering the transfer functions from
input numbers 1 and 4 (with input names 'power' and 'speed') to output
number 3 (with name 'position')
For a single output model m
m4 = m(inputs)
will select the corresponding input channels, and for a single input model
m5 = m(outputs)
(3-53)
3-54
Examining Models
(3-54)
For the model m in (Equation 3-53), the restriction to the transfer function
matrix G is obtained by
m1 = m('measured') or just m1 = m('m')
will plot the additive noise spectra according to the model m, while
bode(m)
This creates a model m3 with all input channels, both measured u and noise
sources e, being treated as measured signals,. That is, m3 is a model from u and
e to y, describing the transfer functions G and H. The information about the
variance of the innovations e is then lost. For example, studying the step
response from the noise channels, will then not take into consideration how
large the noise contributions actually are.
3-55
Tutorial
ny inputs)
fcn(noisecnv(m)) returns the properties of the transfer function [G H] (ny
returns the properties of the transfer function HL. (ny outputs and ny inputs).
If m is a time series model, fcn(m) returns the properties of H , while
fcn(noisecnv(m,'Norm'))
3-56
Examining Models
idmodel Properties
See Function Reference for a complete list of idmodel properties.
Adding Channels
m = [m1,m2,...,mN]
creates an idmodel object m, consisting of all the input channels in m1,... mN.
The output channels of mk must be the same. Analogously
m = [m1;m2;... ;mN]
creates an idmodel object m consisting of all the output channels in m1, m2,...,
mN. The input channels of mk must all be the same.
If you have the Control System Toolbox, interconnections between idmodels,
like G1+G2, G1*G2, append(G1,G2), feedback(G1,G2), etc, can be performed
just as for LTI-objects. However, covariance information is typically lost.
This gives G and v in (Equation 3-11) along with their estimated covariances,
which are translated from the covariance matrix of the estimated parameters.
If m corresponds to a time-continuous model, the frequency functions are
computed accordingly.
3-57
Tutorial
shows a comparison of several models. Modk can be any idmodel models. They
can be used with any of the Control System Toolboxs LTI models. For some
commands Modk can also be idfrd and iddata objects. For multivariable
models, the plots are grouped so that each input/output channel (for all models)
are plotted together. The InputName and OutputName properties of the models
are used for this. The number of channels need not be the same in the different
3-58
Examining Models
allows you to define colors, linestyles and markers associated with the different
models. PlotStyle takes values such as 'b' (for blue), 'b:' (for a blue dotted
line) or 'b*-' (for a blue solid line with the points marked by a star). This is
the same as for the usual plot command.
To also show the uncertainty of the model characteristics, use
command(Mod1,...,ModN,'sd',SD)
This will mark, using dash-dotted lines, a confidence region around the
nominal model corresponding to SD standard deviations (for the Gaussian
distribution). This region is computed using the estimated covariance matrix
for the estimated parameters.
command(Mod1,...,ModN,'sd',SD,'fill')
For the frequency response graphs, this shows the additive disturbance
spectra, i.e., the spectra of the signal H(q)e(t) in Equation (Equation 3-53), so
that the properties of the noise source e are included in the plot.
For the other graphs, the properties of the transfer function H are shown, i.e.,
no noise normalization is done. The same is true if Model is a time series and
has no measured input channels. That means that, for example, step shows
the step response of the transfer function H, without accounting for the size
(covariance matrix) of e. To include such effects, the disturbances should first
be converted to normalized noise sources, using the command noisecnv. See
The Noise Channels.
3-59
Tutorial
Model Output
An important and visually suggestive plot, is to compare the measured output
signal with the models simulated or predicted outputs. This is achieved by
compare(Data,model)
The input signal in Data is used by the model(s) to simulate the output. This
simulated output is shown together with the measured output, which reveals
what features in the data the model can and cannot reproduce. Also a legend
shows the fit between the signals, in terms of how much of the output variation
is reproduced by the model(s).
Frequency Response
Three functions offer graphic display of the frequency functions and spectra:
bode, ffplot, and nyquist.
bode(G)
Transient Response
The impulse and step responses of the models are shown by
impulse(Model) and step(Model)
impulse and step follow the general syntax, but can also accept iddata objects
as arguments. For direct estimation of step and impulse responses from data,
the procedure described in Estimating Impulse Responses is used.
3-60
Examining Models
This gives a plot with x marking poles and o marking zeros. Otherwise, pzmap
follows the general syntax.
General
If you have the Control System Toolbox
view(Model)
The two last commands were described previously. The three first commands
clearly transform to the state-space, the polynomial, and the multivariable
3-61
Tutorial
or directly by
mc = pem(Data,nx,'Ts',0,'ss','can')
3-62
Examining Models
The search for the continuous-time model must be carried out in a canonical (or
any other structured) parameterization. The fit is still made to the sampled
signals in Data. The model is sampled with the datas sampling interval for the
fit. The information about the input intersample behavior in Data.
InterSample is used to determine whether the sampling should be
zero-order-hold (zoh, piecewise constant input) or first-order-hold (foh,
piecewise linear input). All this gives black-box state-space models without any
prescribed internal structure. In these cases, and for a zoh input, it may be
easier to first estimate a black-box model in discrete time and then transform
it to continuous time with d2c as described below. For a foh input it might be
better to directly estimate the continuous-time model, since the mapping from
discrete to continuous under a foh assumption is somewhat tricky.
The major reason for identifying continuous-time model is to secure a
particular structure of the continuous-time state-space matrices. This would
typically reflect a physical interpretation or some greybox modeling work done.
This situation is handled by defining the structure as a continuous time idss
or idgrey model, as described inBlack-Box State-Space Models: the idss
Model and onwards. The resulting structure mi is fitted to data in the usual
way.
m = pem(Data,mi)
Transformations
Transformations between continuous-time and discrete-time model
representations are achieved by c2d and d2c. Note that it is not sufficient to
just assign a new value of Ts to the model object. The corresponding
uncertainty measure (the estimated covariance matrix of the internal
parameters) is also transformed in most cases. The syntax is
modc = d2c(modd)
modd = c2d(mc,T)
If the discrete-time model has some pure time delays (nk > 1 ) the default
command removes them before forming the continuous-time model, and
appends them using the property InputDelay in model modc. This property is
used to add appropriate phase lag and shift the data, whenever the model is
used. d2c also offers as an option to approximate the dead time by a finite
dimensional system. Note that the disturbance properties are translated by the
somewhat questionable formula (Equation 3-29). The covariance matrix is
translated by Gauss approximation formula using numerical derivatives. The
3-63
Tutorial
M-file nuderst is then invoked. You may want to edit it for applications where
the parameters have very different orders of magnitude. See the comments in
State-Space Structures: Initial Values and Numerical Derivatives.
Here is an example that compares the Bode plots of an estimated model and its
continuous-time counterpart.
m= armax(Data,[2 3 1 2]);
mc = d2c(m); bode(m,mc)
3-64
This invokes a plot from which a best order is chosen. Just omitting the order
argument, m = n4sid(Data) or pem(Data) makes a default choice of the best
order.
For models of ARX type, various orders and delays can be efficiently studied
with the command arxstruc. Collect in a matrix NN all of the ARX structures
you want to investigate, so that each row of NN is of the type
[na nb nk]
With
V = arxstruc(Date,Datv,NN)
an ARX model is fitted to the data set Date for each of the structures in NN.
Next, for each of these models, the sum of squared prediction errors is
computed, as they are applied to the data set Datv. The resulting loss functions
are stored in V together with the corresponding structures.
To select the structure that has the smallest loss function for the validation set
Datv, use
3-65
Tutorial
nn = selstruc(V,0)
where you first establish a suitable value of the delay nk by testing second
order models with delays between one and ten. The best fit selects the delay,
and then all combinations of ARX models with up to five a and b parameters
are tested with delays around the chosen value (a total of 75 models).
If the model is validated on the same data set from which it was estimated; i.e.,
if Date = Datv, the fit always improves as the flexibility of the model structure
increases. You need to compensate for this automatic decrease of the loss
functions. There are several approaches for this. Probably the best known
technique is Akaikes Final Prediction Error (FPE) criterion and his closely
related Information Theoretic Criterion (AIC). Both simulate the cross
validation situation, where the model is tested on another data set.
3-66
If substantial noise is present, the ARX models may need to be of high order to
describe simultaneously the noise characteristics and the system dynamics.
(For ARX models the disturbance model 1/A(q) is directly coupled to the
dynamics model B(q)/A(q).)
3-67
Tutorial
shows a nonparametric estimate of the impulse response. In the call above, also
a confidence region around zero is shown, corresponding to three standard
deviations (ca 99.9%). Any part of the impulse response that is outside this
region is thus significant. The first sample after t=0, at which the impulse
response estimate crosses the confidence band is thus a good estimate of the
delay in the channel in question.
Significant impulse response estimates for negative time lags are indications
of feedback in the data.
Residual Analysis
The residuals associated with the data and a given model, as in
(Equation 3-38), are ideally white and independent of the input for the model
to correctly describe the system. The function
resid(Model,Data)
3-68
computes the residuals (prediction errors) e from the model when applied to
Data, and performs whiteness and independence analyses. The auto
correlation function of e and the cross-correlation function between e and u are
computed and displayed for up to lag 25. Also displayed are 99% confidence
intervals for these variables, assuming that e is indeed white and independent
of u.
The rule is that if the correlation functions go significantly outside these
confidence intervals, do not accept the corresponding model as a good
description of the system. Some qualifications of this statement are necessary:
Model structures like the OE structure (Equation 3-17) and methods like the
IV method (Equation 3-41) focus on the dynamics G and less about the
disturbance properties H. If you are interested primarily in G, focus on the
independence of e and u rather than the whiteness of e.
Correlation between e and u for negative lags, or current e ( t ) affecting
future u ( t ) , is an indication of output feedback. This is not a reason to reject
the model. Correlation at negative lags is of interest, since certain methods
do not work well when feedback is present in the input-output data, (see
Feedback in Data), but concentrate on the positive lags in the
cross-correlation plot for model validation purposes.
When using the ARX model (Equation 3-14), the least squares procedure
automatically makes the correlation between e ( t ) and u ( t k ) zero for
k = nk , nk + 1 , nk + nb 1 , for the data used for the estimation.
The residuals e together with the input u are returned by
E= resid(Model,Data)
as an iddata object. As part of the validation process, you can graph the
residuals using
plot(E)
for a simple visual inspection of irregularities and outliers. (See also Outliers
and Bad Data; Multi-Experiment Data.)
3-69
Tutorial
returns the iddata object e which has the inputs in Data as inputs and the
prediction errors (residuals) as outputs. Building models using e will thus
reveal if there is any significant influence from u to e left in the data. Such
models are called Model Error Models, and examining them is a good
complement to traditional residual analysis.
E= resid(Model,Data)
impulse(E,'sd',3) % An alternative to residual analysis
bode(spa(E),'sd',3) % Shows the frequency ranges
% with significant model errors
m = arx(E,[0 10 0])
bode(m,'sd',3)
Note that the resid command has several options to display model error
properties rather than correlation functions.
Noise-Free Simulations
To check whether a model is capable of reproducing the observed output when
driven by the actual input, you can run a simulation.
u = Data(:,[],:) % Extracting the input from the data
yh = sim(Model,u)
y = Data(:,:,[]) % extracting the output from the data
plot(y,yh)
3-70
sequence. It does not account for systematic errors due to an inadequate choice
of model structure. There is no guarantee that the true system lies in the
confidence interval.
The uncertainty in the different model views is displayed if the argument 'sd'
is included in the argument list
command(Model,'sd',sd)
Ten possible models are drawn from the asymptotic distribution of the model
Model. The response of each of them to the input u is graphed on the same
diagram.
The uncertainty of these responses concerns the external, input-output
properties of the model. It reflects the effects of inadequate excitation and the
presence of disturbances.
You can also directly get the standard deviation of the simulated result by
[ysim,ysimsd] = sim(Model,u)
which is used to give the standard deviations of all model characteristics. The
parametric uncertainty is directly available as
Model.da for the standard deviations of Model.a
3-71
Tutorial
Modelc = Model
Modelc.ss = 'canon'
Modelc.da
All routines for computing and displaying model characteristics also offer to
calculate and show the uncertainties. See Transformations to Other Model
Representations.
Large uncertainties in these representations are caused by excessively high
model orders, inadequate excitation, or bad signal-to-noise ratios.
3-72
Model Structures
Multivariable system offer a potentially richer internal structure. The easiest
approach, in the black-box situation, is to think just in terms of input delays
and state-space model order. A recommended approach is to get an idea of
input delays from the nonparametric impulse response estimate, and
determine the vector nk = [nk1,nk2,...,nkm] where nkj is the minimal delay
from input j to any of the output channels, and then try state-space models
with several orders and with these delays.
impulse(Data,'sd',3)
Model = n4sid(Data(1:500),'nx',1:10,'nk',nk)
compare(Data(501:1000),Model)
The compare plot will reveal which output channels are easy and which are
difficult to reproduce.
An alternative to find the delays is to first estimate a parametric model with
delays 1, and then examine the impulse responses of this model and determine
the delays.
Model = pem(Data) % This uses 'best' model order
impulse(Model,'sd',3)
Model = pem(Data,'nx',1:10,'nk',nk)
3-73
Tutorial
Channel Selection
A particular aspect of multivariable models regards the selection of channels.
Models for subselections of input-output channels may be quite useful and
informative. Generally speaking the models become better when more input
channels are used, and worse when more output channels are used. The
latter observation is due to the fact that such models have more to explain.
If you build models with several outputs and find, using compare, a certain
output channel to be difficult to reproduce, then try to build model of this
channel alone. This will reveal if there are inherent difficulties with this
output, or that it is just too difficult to handle it together with other outputs.
Analogously, if you see that, using, e.g., step or impulse, a certain input
channel seems to have an insignificant influence on the outputs, then remove
that channel, and examine if the corresponding model becomes any worse, e.g.,
in the compare plots.
A main feature of the toolboxs data and model objects is that it gives full
support for the bookkeeping required for these channel subselections.
Channels are selected by direct subreferencing, and the InputName and
OutputName properties form the basis for a correct combination of channels.
The subreferencing follows
Data(Samples,Outputs,Inputs)
Model(Outputs,Inputs)
3-74
Date = Data(1:500)
Datv = Data(501:1000)
m = pem(Date)
compare(Datv,m)
m1 = pem(Date(:,3,4))
compare(Datv,m,m1)
bode(m,m1)
compare(Datv,m(:,4),m1)
3-75
Tutorial
Offset Levels
When the data have been collected from a physical plant, they are typically
measured in physical units. The levels in these raw input and output
measurements may not match in any consistent way. This will force the models
to waste some parameters correcting the levels.
Typically, linearized models are sought around some physical equilibrium. In
such cases offsets are easily dealt with: subtract the mean levels from the input
and output sequences before the estimation. It is best if the mean levels
correspond to the physical equilibrium, but if such values are not known, use
the sample means.
Data = detrend(Data);
Section 14.1 in Ljung (1999) discusses this in more detail. There are situations
when it is not advisable to remove the sample means. It could for example be
that the physical levels are built into the underlying model, or that
integrations in the system must be handled with the right level of the input
being integrated.
With the detrend command, you can also remove piece-wise linear trends.
3-76
Often the data has portions with bad behavior. This may, e.g., be due to big
disturbances or sensor failures over a period of time. It can also be that there
are time periods where nothing happens, the input is not exciting, etc. Then
the best alternative is to break up the data into pieces of informative portions.
By merging the pieces into a multiexperiment iddata object, they can still be
used together to build models. Another situation when multiexperiment data
is useful is when several different experiments have been performed on the
same plant. The estimation routines take proper action to handle the different
pieces. All estimation, simulation, and validation routines in the toolbox
handle multi-experiment data in a transparent fashion. A typical string of
commands could be
plot(Data)
Datam = merge(Data(1:340),Data(500:897), ...
Data(1001:1200),Data(1550:2000))
m =pem(Datam{[1,2,4]}) % Portions 1,2 and 4 for estimation
compare(Datam{3},m) % Portion 3 for validation
Missing Data
In practice it is often the case that certain measurement samples are missing.
The reason may be sensor failures or data acquisition failures. It may be that
the data are directly reported as missing, or that plots reveal that some values
are obviously in error. This may apply both to inputs and outputs. In these
cases, replace the missing data by NaN when forming the signal matrices and
the iddata object. The routine misdata can then be applied to reconstruct the
missing data in a reasonable way.
dat = iddata(y,u,0.2) % y and/or u contain NaN for missing data.
dat1 = misdata(dat);
plot(dat,dat1) % Checking how the missing data
% have been estimated in dat1
m = pem(dat1) % Model estimated using reconstructed missing data
3-77
Tutorial
This computes and uses a fifth order Butterworth bandpass filter with
passband between the indicated frequencies. The data is filtered through this
filter before fitting the transfer function from the measured inputs (G in
Equation (Equation 3-53)) to the outputs. The disturbance model (H) is
however estimated using the unfiltered data. Chapter 14 in Ljung (1999)
discusses the role of filtering in more detail.
The command butter is from the Signal Processing Toolbox. If you do not have
that toolbox, the filter can be computed using idfilt from the System
Identification Toolbox.
[Df,Mfilt] = idfilt(Data,5,[0.02 0.1])
m = pem(Data,3,'Foc',Mfilt)
For a model that does not use a disturbance description (that is, H = 1 in
(Equation 3-53), which corresponds to K = 0 for state-space, and na=nc=nd=0
for polynomial models), the Focus effect is the same as applying the routine to
filtered data. That is,
m = pem(Data,3,'Foc',Mfilt,'dist','none')
m = pem(Df,3,'dist','none')
Feedback in Data
If the system was operating in closed loop (feedback from the past outputs to
the current input) when the data were collected, some care has to be exercised.
3-78
Basically, all the prediction error methods work equally well for closed-loop
data. Note, however, that the Output-Error model (Equation 3-17) and the
Box-Jenkins model (Equation 3-18) are normally capable of giving a correct
description of the dynamics G, even if H (which equals 1 for the output error
model) does not allow a correct description of the disturbance properties. This
is no longer true for closed-loop data. You then need to model the disturbance
properties more carefully. Another thing to be cautious about is that impulse
response effects at delay 0 very well could be traced to the feedback mechanism
and not to the system itself.
The spectral analysis method and the instrumental variable techniques (with
default instruments) as well as n4sid may give unreliable results when
applied to closed-loop data. These techniques should be avoided when feedback
is present.
To detect if feedback is present, use the basic method of applying impulse to
estimate the impulse response. Significant values of the impulse response at
negative lags is a clear indication of feedback. When a parametric model has
been estimated and the resid command is applied, a graph of the correlation
between residuals and inputs is given. Significant correlation at negative lags
again indicates output feedback in the generation of the input. Testing for
feedback is, therefore, a natural part of model validation.
Delays
The selection of the delay nk in the model structure is a very important step in
obtaining good identification results. You can get an idea about the delays in
the system by the impulse response estimate from impulse.
Incorrect delays are also visible in parametric models. Underestimated delays
(nk too small) show up as small values of leading b k estimates, compared to
their standard deviations. Overestimated delays (nk too large) are usually
visible as a significant correlation between the residuals and the input at the
lags corresponding to the missing b k terms in the resid plot.
A good procedure is to start by using arxstruc to test all feasible delays
together with a second-order model. Use the delay that gives the best fit for
further modeling. When you have found an otherwise satisfactory structure,
vary nk around the nominal value within the structure, and evaluate the
results.
3-79
Tutorial
(3-55)
(3-56)
y ( t ) = ( t ) 0 ( t ) + e ( t )
(3-57)
where the regression vector ( t ) contains old values of observed inputs and
outputs, and 0 ( t ) represents the true description of the system. Moreover,
e ( t ) is the noise source (the innovations). Compare with (Equation 3-14). The
T
natural prediction is y ( t ) = ( t ) ( t 1 ) and its gradient with respect to
becomes exactly ( t ) .
3-80
For models that cannot be written as linear regressions, you cannot recursively
compute the exact prediction and its gradient for the current estimate ( t 1 ) .
Then approximations y ( t ) and ( t ) must be used instead. Section 11.4 in
Ljung (1999) describes suitable ways of computing such approximations for
general model structures.
The matrix Q ( t ) that affects both the adaptation gain and the direction in
which the updates are made, can be chosen in several different ways. This is
discussed in the following.
(3-58)
Ew ( t )w ( t ) = R 1
(3-59)
(3-60)
T
( t 1 ) ( t ) ( t ) P ( t 1 )
P ( t ) = P ( t 1 ) + R1 P
-------------------------------------------------------------------T
R 2 + ( t ) P ( t 1 ) ( t )
Here R 2 is the variance of the innovations e ( t ) in (Equation 3-57):
2
R 2 = Ee ( t ) (a scalar). The algorithm (Equation 3-60) will be called the
Kalman filter (KF) approach to adaptation, with drift matrix R 1 . See eq
3-81
Tutorial
tk 2
e (k)
(3-61)
k=1
+ ( t ) P ( t 1 ) ( t )
(3-62)
3-82
(3-63)
Q ( t ) = ----------------- I
2
(t )
(3-64)
See eqs (11.45) and (11.46), respectively in Ljung (1999). These choices of Q ( t )
move the updates of in (Equation 3-55) in the negative gradient direction
(with respect to ) of the criterion (Equation 3-39). Therefore, (Equation 3-63)
will be called the Unnormalized Gradient (UG) approach and
(Equation 3-64) the Normalized Gradient (NG) approach to adaptation, with
gain . The gradient methods are also known as least mean squares (LMS) in
the linear regression case.
Available Algorithms
The System Identification Toolbox provides the following functions that
implement all common recursive identification algorithms for model structures
in the family (Equation 3-43): rarmax, rarx, rbj, rpem, rplr, and roe. They all
share the following basic syntax.
[thm,yh] = rfcn(z,nn,adm,adg)
gives the forgetting factor algorithm (Equation 3-62), with forgetting factor
lam.
adm = 'ug'; adg = gam
gives the unnormalized gradient approach (Equation 3-63) with gain gam.
Similarly,
adm = 'ng'; adg = gam
gives the normalized gradient approach (Equation 3-64). To obtain the Kalman
filter approach (Equation 3-60) with drift matrix R1, enter
adm = 'kf'; adg = R1
The value of R 2 is always 1. Note that the estimates in (Equation 3-60) are
not affected if all the matrices R 1, R 2 and P ( 0 ) are scaled by the same
3-83
Tutorial
number. You can therefore always scale the original problem so that R 2
becomes 1.
The output argument thm is a matrix that contains the current models at the
different samples. Row k of thm contains the model parameters, in
alphabetical order at sample time k, corresponding to row k in the data matrix
z. The ordering of the parameters is the same as m.par would give when
applied to a corresponding offline model.
The output argument yh is a column vector that contains, in row k, the
predicted value of y ( k ) , based on past observations and current model. The
vector yh thus contains the adaptive predictions of the outputs, and is useful
also for noise cancelling and other adaptive filtering applications.
The functions also have optional input arguments that allow the specification
of ( 0 ), P ( 0 ) , and ( 0 ) . Optional output arguments include the last value of
the matrix P and of the vector .
Now, rarx is a recursive variant of arx; similarly rarmax is the recursive
counterpart of armax and so on. Note, though that rarx does not handle
multi-output systems, and rpem does not handle state-space structures.
The function rplr is a variant of rpem, and uses a different approximation of
the gradient . It is known as the recursive pseudo-linear regression approach,
and contains some well known special cases. See Equation (11.57) in Ljung
(1999). When applied to the output error model (nn=[0 nb 0 0 nf nk]) it
results in methods known as HARF ('ff'-case) and SHARF ('ng'-case). The
common extended least squares (ELS) method is an rplr algorithm for the
ARMAX model (nn=[na nb nc 0 0 nk]).
The following example shows a second order output error model, which is built
recursively and its time varying parameter estimates plotted as functions of
time.
thm = roe(z,[2 2 1],'ff',0.98);
plot(thm)
The next example shows how a second order ARMAX model is recursively
estimated by the ELS method, using Kalman filter adaptation. The resulting
static gains of the estimated models are then plotted as a function of time.
[N,dum]=size(z);
thm = rplr(z,[2 2 2 0 0 1],'kf',0.01eye(6));
nums = sum(thm(:,3:4)')';
3-84
dens = ones(N,1)+sum(thm(:,1:2)')';
stg = nums./dens;
plot(stg)
Segmentation of Data
Sometimes the system or signal exhibits abrupt changes during the time when
the data is collected. It may be important in certain applications to find the
time instants when the changes occur and to develop models for the different
segments during which the system does not change. This is the segmentation
problem. Fault detection in systems and detection of trend breaks in time
series can serve as two examples of typical problems.
The System Identification Toolbox offers the function segment to deal with the
segmentation problem. The basic syntax is
thm = segment(z,nn)
3-85
Tutorial
with a format like rarx or rarmax. The matrix thm contains the piecewise
constant models in the same format as for the algorithms described earlier in
this section.
The algorithm that is implemented in segment is based on a model description
like (Equation 3-58), where the change term w ( t ) is zero most of the time, but
now and then it abruptly changes the system parameters 0 ( t ) . Several
Kalman filters that estimate these parameters are run in parallel, each of them
corresponding to a particular assumption about when the system actually
changed. The relative reliability of these assumed system behaviors is
constantly judged, and unlikely hypotheses are replaced by new ones. Optional
arguments allow the specification of the measurement noise variance R 2 in
(Equation 3-57), of the probability of a jump, of the number of parallel models
in use, and also of the guaranteed lifespan of each hypothesis. See segment in
the Command Reference chapter.
3-86
Spectral analysis (etfe and spa) returns results in the idfrd model format,
that now just contains SpectrumData and its variance. bode will only plot these
signal spectra and, if required, the confidence intervals.
g = spa(y)
p= etfe(y)
bode(g,p,'sd',3)
3-87
Tutorial
Note that arx also handles multivariable signals, and so do n4sid and pem.
m = n4sid(y) %default order
bode(m)
compare(y,m,10) % 10-step ahead predictions being evaluated.
In addition there are two commands that are specifically constructed for
building scalar AR models of time series. One is
m = ar(y,na)
which has an option that allows you to choose the algorithm from a group of
several popular techniques for computing the least squares AR model. Among
these are Burgs method, a geometric lattice method, the Yule-Walker
approach, and a modified covariance method. See Function Reference for
details. The other command is
m = ivar(y,na)
3-88
Periodic Inputs
It is often an advantage to use a periodic input for identification, whenever
possible. See Section 13.3 in Ljung(1999). If you import or create a periodic
input, as in
u = idinput([300 2 5]) % Period 300, 2 inputs, 5 periods
3-89
Tutorial
Function Calls
The function calls are the same for many essential functions. bode, nyquist,
step, impulse, ssdata, tfdata, zpkdata, freqresp, minreal, etc., all do the
same things with essentially the same syntax. The System Identification
Toolbox commands however also handle model uncertainty. The System
Identification Toolbox commands will be used whenever at least one of the
objects in the argument list is an idmodel or idfrd object.
Also, subreferencing channels and concatenations follow the same syntax.
Moreover, most of the LTI-commands for model manipulation, like G1+G2,
G1*G2, feedback, append, balreal, augstate, canon, etc, will work (using the
Control System Toolbox) in the expected way, returning idmodel objects.
However, covariance information is in most cases lost.
Object Relations
Since the System Identification Toolbox can be run without the Control System
Toolbox, there are no formal parent/child relations between the objects in the
two toolboxes. There are however easy transformations between them. The
command that creates idmodel, idss, and idpoly will accept any LTI object,
zpk, tf or ss. idfrd can similarly be created from frd objects. If the LTI object
has an InputGroup named 'noise' these input will be treated as normalized
white noise, when creating the idmodel object with correct disturbance model
information.
Analogously, ss, zpk, tf, and frd accept any idmodel or idfrd (in case of frd)
object. The covariance information is then not stored in the LTI objects, but all
disturbance information will be translated to a group of extra input channels
with the group name 'noise'. If these are interpreted as normalized white
noise, the LTI objects have the same disturbance properties as the original
imdmodel object.
These simple relations also mean that it is easy to use any LTI command in the
Control System Toolbox and return to System Identification Toolbox objects.
Mb = idss(balreal(ss(M)))
Plot Relations
Although the calls bode, step etc., have essentially the same syntax, the plots
look different. The System Identification Toolbox commands show, when
required, confidence regions, and typically show the different input/output
3-90
Memory/Speed Trade-Offs
On machines with no formal memory limitations, it is still of interest to
monitor the sizes of the matrices that are formed. The typical situation is when
an overdetermined set of linear equations is solved for the least squares
solution. The solution time depends, of course, on the dimensions of the
corresponding matrix. The number of rows corresponds to the number of
observed data, while the number of columns corresponds to the number of
estimated parameters. The property MaxSize used with all the relevant
M-files, guarantees that no matrix with more than MaxSize elements is
formed. Larger data sets and/or higher order models are handled by for loops.
for loops give linear increase in time when the data record is increased, plus
some overhead.
If you regularly work with large data sets and/or high order models, it is
advisable to tailor the memory and speed trade-off to your machine by choosing
MaxSize carefully. You could also change the default value of MaxSize in the
M-file idmsize. Then the default value of MaxSize (that is 'Auto') will be
tailored to your needs. Note that this value is allowed to depend on the number
of rows and columns of the matrices formed.
3-91
Tutorial
Local Minima
The iterative search procedures in pem, armax, oe, and bj lead to models
corresponding to a local minimum of the criterion function (Equation 3-39).
Nothing guarantees that this local minimum is also a global minimum. The
start-up procedure for black-box models in these routines is, however,
reasonably efficient in giving initial estimates that lead to the global minimum.
If there is an indication that a minimum is not as good as you expected, try
starting the minimization at several different initial conditions, to see if a
smaller value of the loss function can be found. The function init can be used
for that.
Another example is when you want to try a model with one more delay (for
example, three instead of two) because the leading b-coefficient is quite small.
m1 = armax(Data,[3 3 2 2]);
m1.b(3) = 0
m2 = armax(Data,m1);
If you decrease the number of delays, remember that leading zeros in the
B-polynomial are treated as delays. Suppose you go from three to two delays
in the above example.
m1 = armax(z,[3 3 2 3]);
m1.b(3) = 0.00001;
m2 = armax(Data,m1);
3-92
Note that when constructing homemade initial conditions, the conditions must
correspond to a stable predictor (C and F being Hurwitz polynomials), and they
should not contain any exact pole-zero cancellations.
For user defined structured state-space and multi-output models, you must
provide the initial parameter values (initial model) when defining the
structure in idss or idgrey. The basic approach is to use physical insight to
choose initial values of the parameters with physical significance, and try some
different (randomized) initial values for the others. The routine init can be
used for that.
Initial State
The filter that computes the prediction errors in (Equation 3-36) needs to be
properly initialized. For input-output (polynomial) models, values of inputs,
outputs and predictions prior to time t = 0 are required, and state-space models
need the initial state x(0). There are several ways to handle these unknown
states. A simple one is to take all unknown values as zero. If the model
predictor has slow dynamics (i.e. the poles of CF, or the eigenvalues of A-KC
are close to the unit circle), this could have a very bad effect on the parameter
estimates. It is particularly pronounced for output-error models, where the
noise model cannot be adjusted to handle slow transients form initial
conditions.
The toolbox offers a number of options how to deal with the initial state of the
predictor. They are handled by the model property InitialState. The
unknown state can be treated as a vector of unknown parameters
(InitialState = 'Estimate'), they can be set to zero (InitialState =
'Zero'), or estimated by a backwards prediction method (InitialState =
'Backcast'). It can also be fixed to any user defined value. The default value
is InitialState = 'Auto', which makes an automatic choice between the
options, guided by the estimation data. For details, see idss and idpoly in the
Command Reference chapter. Basically, the effect of the initial conditions on
the prediction errors are tested and if they seem to be negligible, 'zero' is
chosen, which gives a fast and efficient algorithm. Otherwise the initial state
is estimated or backcasted. EstimationInfo will contain information about
which method was chosen in this case.
Proper handling of the initial state is necessary both when models are
estimated, and when predictions and simulations are compared. The
3-93
Tutorial
commands predict, pe, sim, and compare all offer options for how to deal with
this.
No Covariance
Evaluating and visualizing the uncertainty of the estimated models is a very
important aspect of system identification. Handling, and translating
covariance information takes a major part of the time in many of the routines
of the System Identification Toolbox. For example, in n4sid, calculating the
Cramer-Rao bound (which in this case is used and an indication of the
covariance properties) takes much longer than estimating the actual model. In
d2c and c2d, most of the time is spent on covariance handling. If you build
models that are of a preliminary nature, and you would like to speed up the
calculations, you can add the property name/property value pair
'Covariance'/'None' to the list of arguments in most relevant routines. This
will prevent covariance calculations and set a flag not to spend time on this in
future use of the model. This flag can also be set in the model at any time by
3-94
Model.cov = 'no'
nk and InputDelay
Whats the difference between the properties nk and InputDelay? InputDelay
is defined for all idmodel and idfrd objects, while nk is defined for idarx,
idpoly as well as for 'Free' and 'Canonical' idss models. Both properties
indicate a delay from the input channels to the outputs. For idarx, nk is a
matrix, describing the delays in the different input/output channels, but
otherwise both nk and InputDelay describe the delay from a certain input
channel to all the output channels.
InputDelay is really a flag that tells the model to append the input delays as
time lags, when the model is simulated, or as phase lags when the frequency
functions are computed. The InputDelay does not show up when the model is
represented in state-space form, nor as transfer functions, nor in the
input-output polynomials. InputDelay can be used both for continuous and
discrete time models. In the latter case, the InputDelay is measured in number
of samples. Moreover InputDelay may assume negative values, in order to
handle noncausal models.
The property nk, on the other hand, is a model structure property, requiring the
model to contain the indicated number of delays whatever the parameter
values. This means that the state-space matrices, the transfer functions, etc.,
will show these delays in an explicit manner. Consequently, nk is not defined
for continuous-time models.
Otherwise the two properties can be used in the same way
m1 = pem(Data,4,'InputDelay',[3 2 0])
m2 = pem(Data,4,'nk',[3 2 0])
bode(m1,m2)
A1 = m1.A
A2 = m2.A
give identical bode-plots (up to minor variations due to end-effects in the data
records), while A1 and A2 are different. In fact while A1 is of size 4-by-4, the
matrix A2 is of size 7-by-7, since three extra states are required to accommodate
the extra 2+1 input delays.
Note that setting nk to a certain value for a given model gives a model structure
that has the indicated delay for any parameter values. The impulse response
of the model may however change (not only be shifted) by this assignment.
3-95
Tutorial
(3-65)
y ( t ) = b0 + b1 u ( t ) + b 2 u ( t ) + b3 u ( t )
(3-66)
let
Data = iddata(y,[ones(size(u)), u, u.^2, u.^3]);
m= arx(Data,'na',0,'nb',[1 1 1 1],'nk',[ 0 0 0 0])
This is formally a model with one output and four inputs, but all the model
testing in terms of compare, sim, and resid operate in the natural way for the
model (Equation 3-65), once the data set Data is defined as above.
Note that when pem is applied to linear regression structures, by default a
robustified quadratic criterion is used. The search for a minimum of the
criterion function is carried out by iterative search. Normally, use this
robustified criterion. If you insist on a quadratic criterion, then set the
argument LimitError in pem to zero. Then pem also converges in one step.
y ( ) = T
k = M
where
3-96
R y ( kT )e
iT
WM( k )
(3-67)
1
( kT ) = ---R
y
N
y ( lT kT )y ( lT )
l=1
1
Ey ( t ) = -----2
2
y ( ) d
(3-68)
= T
spa(y);
= squeeze(sp.spec) % squeeze takes out the spurios dimensions
sum(phiy)/length(phiy)/T;
sum(y.^2)/size(y,1);
3-97
Tutorial
Ee ( t ) =
iT 2
(3-69)
3-98
and these functions. Modify the C-polynomial accordingly. Make the degree of
the monic C-polynomial in continuous time equal to the sum of the degrees of
the monic A- and D-polynomials; i.e., in continuous time.
length(C)-1 = (length(A)-1)+(length(D)-1)
3-99
Tutorial
nu
nu
nu
f 1, f nf1, , f 1 , , f nfnu ]
Note that the property PName (for Parameter Name) may be useful to help with
the bookkeeping in these cases, and when fixing certain parameters using
FixedParameter. The routine setpname may be helpful in automatically setting
mnemonic parameter names for black-box models.
Complex-Valued Data
Some applications of system identification work with complex-valued data, and
thus create complex-valued models. Most of the routines in the System
Identification Toolbox support complex data and models. This is true for the
estimation routines ar, armax, arx, bj, covf, ivar, iv4, oe, pem, spa, and n4sid.
3-100
The transformation routines, like freqresp, zpkdata, etc., also work for
complex-valued models, but no pole-zero confidence regions are given. Note
also that the parameter variance-covariance information then refers to the
complex valued parameters, so no separate information about the accuracy of
the real and imaginary parts will be given. Some display functions like compare
and plot do not work for the complex case. Use sim and plot real and
imaginary parts separately.
Strange Results
Strange results can of course be obtained in any number of ways. We only
point out two cautions: It is tempting in identification applications to call the
residuals eps. Dont do that. This changes the machine , which certainly will
give you strange results.
It is also natural to use names like step, phase, etc., for certain variables. Note
though that these variables take precedence over M-files with the same name
so be sure you dont use variable names that also are names of M-files.
3-101
Tutorial
3-102
4
Function Reference
This section contains detailed descriptions of all of the functions of user interest in the System
Identification Toolbox. It begins with a list of functions grouped by subject area and continues with
the entries in alphabetical order.
Information is also available through the on-line Help facility. By typing a function name without
arguments, you also get immediate syntax help about its arguments for most functions
For ease of use, most functions have several default arguments. The Syntax first lists the function
with the necessary input arguments and then with all the possible input arguments. The functions
can be used with any number of arguments between these extremes. The rule is that missing, trailing
arguments are given default values, as defined in the manual. Default values are also obtained by
entering the arguments as the empty matrix [ ].
MATLAB does not require that you specify all of the output arguments; those not specified are not
returned. For functions with several output arguments in the System Identification Toolbox, missing
arguments are, as a rule, not computed, in order to save time.
Function Reference
Functions By Category
Help Functions
help ident
idhelp
A micro-manual.
idprops,
help idprops
midprefs
pe
predict
sim
Data Manipulation
4-2
detrend
get/set
iddata
idfilt
Filter data.
merge (iddata)
misdata
plot (iddata)
Plot data.
resample
Resample data.
Functions By Category
Nonparametric Estimation
covf
cra
impulse, step
etfe
spa
Parameter Estimation
ar
Estimate AR model.
armax
arx
bj
ivar
iv4
oe
n4sid
pem
idfrd
idgrey
4-3
Function Reference
4-4
idpoly
idss
Functions By Category
init
merge (idmodel)
Model Conversions
arxdata
idmodred
c2d
d2c
freqresp
idfrd
noisecnv
polydata
ssdata
tfdata
zpkdata
Model Analysis
bode
compare
ffplot
impulse, step
nyquist
4-5
Function Reference
4-6
present
pzmap
view
Functions By Category
Model Validation
aic, fpe
arxstruc, selstruc
Select ARX-structure
compare
pe
predict
resid
sim
Simulate a model
bode, nyquist
pzmap
arxdata, polydata,
ssdata, tfdata,
zpkdata
ivstruc
n4sid, pem
selstruc
Select structure.
struc
4-7
Function Reference
rarx
rbj
roe
rpem
rplr
segment
General
4-8
get
set
setpname
size
timestamp
aic
Purpose
4aic
Syntax
am = aic(Model)
Description
2d
AIC = log ( V ) + ------N
where V is the loss function, d is the number of estimated parameters, and N
is the number of estimation data values.
See Also
EstimationInfo, fpe
Reference
4-13
Algorithm Properties
Purpose
4Algorithm Properties
Syntax
idprops algorithm
m.algorithm
Description
All the idmodel objects in the toolbox, idarx, idss, idpoly, and idgrey, have a
property Algorithm, which is a structure that contains a number of options
that govern the estimation algorithms. The fields of this structure can be
individually set and retrieved in the usual way, such as get(m,'MaxIter') or
m.SearchDirection = 'gn'. Also, autofill applies and the names are case
insensitive.
Note The algorithm properties, like all other model properties, will be
inherited by the resulting model m. If you continue the estimation using m as
initial model, all previously set algorithm features will thus apply, unless you
specify otherwise.
4-14
Algorithm Properties
statistical variance point of view, this is the optimal weighting, but then
the approximation aspects (bias) of the fit are neglected.
- 'Simulation': This means that frequency weighting of the transfer
function fit is given by the input spectrum. Frequency ranges where the
input has considerable power will thus be better described by the model.
In other words, the model approximation is such that the model will
produce as good simulations as possible, when applied to inputs with the
same spectra as used for the estimation. For models that have no
disturbance model, that is y = G u + e, (A=C=D=1 for idpoly models and K=0
for idss models) there is no difference between 'Simulation' and
'Prediction'. For models with a disturbance description, i.e. y = Gu + H
e, G is first estimated with H = 1 and then H is estimated by a prediction
fixed. This
error method, keeping the estimated transfer function G
option will also guarantee a stable transfer function G.
- 'Stability': The resulting model is guaranteed to be stable, but a
prediction weighing is still maintained. Note that forcing the model to be
stable could mean that a bad model is obtained. Use only when you know
the system to be stable.
- Any SISO linear filter. Then the transfer function from input to output is
determined by a frequency fit with this filter times the input spectrum as
weighting function. The disturbance model is determined by a prediction
error method, keeping the transfer function estimate fixed, as in the
simulation case. To obtain a good model fit over a special frequency range,
the filter should thus be chosen with a passband over this range. For a
model with no disturbance model, the result is the same as first applying
prefiltering to data using idfilt. The filter can be specified in a few
different ways as:
- Any single-input-single-output idmodel
- An ss, tf or zpk model from the Control System Toolbox
- {A,B,C,D} with the state-space matrices for the filter
- {numerator, denominator} with the transfer function
numerator/denominator of the filter
MaxSize: No matrix with more than MaxSize elements is formed by the
algorithm. Instead, for-loops will be used. MaxSize will thus decide the
memory/speed trade-off, and can prevent slow use of virtual memory.
MaxSize can be any positive integer, but it is required that the input-output
4-15
Algorithm Properties
data contain less than MaxSize elements. The default value of MaxSize is
'Auto', which means that the value is determined in the M-file idmsize. You
can edit this file to optimize speed on a particular computer.
FixedParameter: A list of parameters that will be kept fixed to the
nominal/initial values and not estimated. This is a vector of integers
containing the indices of the fixed parameters. The numbering of the
parameters is the same as in the model property 'ParameterVector'. The
parameter names from the property 'PName' can also be used. For
structured state-space models, it is easier to fix/unfix parameters by the
structure matrices, As, Bs, etc. See idss. When using parameter names to
specify the fixed parameters, Fixedparameter is a cell array of strings. The
strings may contain the wildcards * (meaning any continuation of the given
string) and ? (meaning any character). For example, if all disturbance model
parameters start with k, FixedParameter = {'k*'} will fix all these
parameters. The function setpname may be useful in this context.
4-16
Algorithm Properties
4-17
Algorithm Properties
StepReduction: In the line search used for other directions than LM, the
step is reduced by the factor stepred in each try. Default:
StepReduction = 2.
considered stable if all poles are within the distance Zstability from the
origin. Default 1.01.
4-18
Algorithm Properties
See Also
Reference
For the iterative minimization, see Dennis, J.E. Jr. and R.B. Schnabel.
Numerical Methods for Unconstrained Optimization and Nonlinear Equations,
Prentice Hall, Englewood Cliffs, N.J. 1983.
For a general reference to the identification algorithms, see Ljung (1999),
Chapter 10.
4-19
ar
Purpose
4ar
Syntax
m = ar(y,n)
[m ,refl] = ar(y,n,approach,window,maxsize)
Description
approach. The sum of a least-squares criterion for a forward model and the
analogous criterion for a time-reversed model is minimized.
approach = 'ls': The least-squares approach. The standard sum of squared
4-20
ar
approach.
The combinations of approaches and windowing have a variety of names. The
least-squares approach with no windowing is also known as the covariance
method. This is the same method that is used in the arx routine. The MATLAB
default method, forward-backward with no windowing, is often called the
modified covariance method. The Yule-Walker approach, least-squares plus
pre- and post-windowing, is also known as the correlation method.
See Algorithm Properties for an explanation of the input argument maxsize.
Examples
Compare the spectral estimates of Burg s method with those found from the
forward-backward nonwindowed method, given a sinusoid in noise signal.
y = sin([1:300]') + 0.5*randn(300,1);
y = iddata(y);
mb = ar(y,4,'burg');
mfb = ar(y,4);
bode(mb,mfb)
4-21
ar
See Also
References
4-22
armax
Purpose
4armax
Syntax
m = armax(data,orders)
m = armax(data,'na',na,'nb',nb,'nc',nc,'nk',nk)
m = armax(data,orders,'Property1',Value1,...,'PropertyN',ValueN)
Description
A ( q )y ( t ) = B ( q )u ( t nk ) + C ( q )e ( t )
using a prediction error method.
data is an iddata object containing the output-input data. The model orders
can be specified as (...,'na',na,'nb',nb,...) or by setting the argument
orders to
orders = [na nb nc nk]
The parameters na, nb, and nc are the orders of the ARMAX model, and nk is
the delay. Specifically,
1
na:
A ( q ) = 1 + a1 q
nb:
B ( q ) = b1 + b2 q
nc:
C ( q ) = 1 + c1 q
+ + a na q
na
+ + b nb q
+ + c nc q
nb + 1
nc
where mi is an initial guess at the ARMAX model given in idpoly format. See
Polynomial Representation of Transfer Functions in the Tutorial for more
information.
For multi-input systems, nb and nk are row vectors, such that the k-th entry
corresponds to the order and delay associated with the k-th input.
If data has no input channels and just one output channel (i.e., it is a time
series), then
4-23
armax
Algorithm
4-24
armax
See Also
References
4-25
arx
Purpose
4arx
Syntax
m = arx(data,orders)
m = arx(data,'na',na,'nb',nb,'nk',nk)
m= arx(data,orders,'Property1',Value1,...,'PropertyN',ValueN)
Description
as
orders = [na nb nk]
na:
A ( q ) = 1 + a1 q
nb:
B ( q ) = b1 + b2 q
+ + a na q
na
+ + b nb q
nb + 1
For a time series, data contains no input channels and orders = na. Then an
AR model of order na for y is computed.
A ( q )y ( t ) = e ( t )
Models with several inputs
A ( q )y ( t ) = B 1 ( q )u 1 ( t nk 1 ) + B nu u nu ( t nk nu ) + e ( t )
are handled by allowing nb and nk to be row vectors defining the orders and
delays associated with each input.
4-26
arx
Models with several inputs and several outputs are handled by allowing na,
nb, and nk to contain one row for each output number. See Multivariable
ARX Models: The idarx Model in the Tutorial for exact definitions.
The algorithm and model structure are affected by the property name/property
value list in the input argument.
Useful options are reached by the properties 'Focus', 'InputDelay', and
'MaxSize'.
See Algorithm Properties for details of these properties and possible values
When the true noise term e ( t ) in the ARX model structure is not white noise
and na is nonzero, the estimate does not give a correct model. It is then better
to use armax, bj, iv4, or oe.
Examples
Algorithm
See Also
4-27
arxdata
Purpose
4arxdata
Syntax
[A,B] = arxdata(m)
[A,B,dA,dB] = arxdata(m)
Description
m is the model as an idarx or idpoly model object. arxdata will work on any
idarx model. For idpoly it will give an error unless the underlying model is an
ARX model, i.e., the orders nc=nd=nf=0. (See the reference page for idpoly.)
A and B are returned in the standard multivariable ARX format (see idarx),
describing the model.
y ( t ) + A 1 y ( t 1 ) + A 2 y ( t 2 ) + + A na y ( t na ) =
B 0 u ( t ) + B 1 u ( t 1 ) + + B nb u ( t nb ) + e ( t )
Here A k and B k are matrices of dimensions ny-by-ny and ny-by-nu, respectively
(ny is the number of outputs, i.e., the dimension of the vector y ( t ) and nu is the
number of inputs). See Multivariable ARX Models: The idarx Model in the
Tutorial.
The arguments A and B are 3-D arrays that contain the A matrices and the B
matrices of the model in the following way:
A is an ny-by-ny-by-(na+1) array such that
A(:,:,k+1) = Ak
A(:,:,1) = eye(ny)
Note that A always starts with the identity matrix, and that leading entries in
B equal to zero means delays in the model. For a time series B = [].
dA and dB are the estimated standard deviations of A and B.
See Also
4-28
idarx
arxstruc
Purpose
4arxstruc
Syntax
V = arxstruc(ze,zv,NN)
V = arxstruc(ze,zv,NN,maxsize)
Description
with the same interpretation as described for arx. See struc for easy
generation of typical NN matrices for single-input systems.
Each of ze and zv are iddata objects containing output-input data. Models for
each of the model structures defined by NN are estimated using the 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.
Note that arxstruc is intended for single-output systems only.
The output argument V is best analyzed using selstruc. It contains the loss
functions in its first row. The remaining rows of V contain the transpose of NN,
so that the orders and delays are given just below the corresponding loss
functions. The last column of V contains the number of data points in ze. The
selection of a suitable model structure based on the information in v is
normally done using selstruc. See Model Structure Selection and Validation
in the Tutorial for advice on model structure selection and cross-validation.
See Algorithm Properties for an explanation of maxsize.
4-29
arxstruc
Examples
Compare first to fifth order models with one delay using cross-validation on the
second half of the data set. Then select the order that gives the best fit to the
validation data set.
NN = struc(1:5,1:5,1);
V = arxstruc(z(1:200),z(201:400),NN);
nn = selstruc(V,0);
m = arx(z,nn);
See Also
4-30
bj
Purpose
4bj
Syntax
m = bj(data,orders)
m = bj(data,'nb',nb,'nc',nc,'nd',nd,'nf',nf,'nk',nk)
m = bj(data,orders,'Property1',Value1,'Property2',Value2,...)
Description
B(q)
C(q )
y ( t ) = ------------ u ( t nk ) + ------------- e ( t )
F( q)
D(q)
using a prediction error method.
data is an iddata object containing the output-input data. The model orders
can be specified by setting the argument orders to
orders = [ nb nc nd nf nk]
The parameters nb, nc, nd, and nf are the orders of the Box-Jenkins model and
nk is the delay. Specifically,
1
nf:
F ( q ) = 1 + f1q
nb:
B ( q ) = b1 + b2 q
nc:
C ( q ) = 1 + c1 q
nd:
D ( q ) = 1 + d1 q
+ + f nf q
1
nf
+ + b nb q
+ + c nc q
+ + d nd q
nb + 1
nc
nd
4-31
bj
For multi-input systems, nb, nf, and nk are row vectors with as many entries
as there are input channels. Entry number i then describes the orders and
delays associated with the i-th input.
The structure and the estimation algorithm are affected by any property
name/property value pairs that are set in the input argument list. Useful
properties are 'Focus', 'InitialState', 'Trace', 'MaxIter', 'Tolerance',
'LimitError', and 'FixedParameter'.
See Algorithm Properties and the reference pages for idmodel and idpoly for
details of these properties and their possible values.
bj does not support multi-output models. Use state-space model for this case
(see n4sid and pem).
Examples
Here is an example that generates data and stores the results of the startup
procedure separately.
B = [0 1 0.5];
C = [1 -1 0.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 = [y u];
mi = bj(z,[2 2 2 2 1],'MaxIter',0)
m = bj(z,mi)
m.EstimationInfo
m = bj(z,m); % Continue if m.es.WhyStop shows that maxiter has
% been reached.
compare(z,m,mi)
Algorithm
See Also
4-32
bode
Purpose
Syntax
4bode
bode(m1,m2,m3,..'sd',sd,'mode',mode,'ap',ap,'fill')
Description
bode computes the magnitude and phase of the frequency response of idmodel
and idfrd models. When invoked without left-hand arguments, bode produces
a Bode plot on the screen.
bode(m) plots the Bode response of an arbitrary idmodel or idfrd model m. This
model can be continuous or discrete, and SISO or MIMO. The InputNames and
OuputNames of the models are used to plot the responses for different I/O
channels in separate plots. Pressing the Enter key advances the plot from one
input-output pair to the next one.
If m contains information about both I/O channels and output noise spectra,
only the I/O channels are shown. To show the output noise spectra enter
m('n') ('n' for 'noise') in the model list. Analogously, specific I/O channels can
be selected by the normal subreferencing: m(ky,ku).
4-33
bode
Note that the frequencies cannot be specified for idfrd objects. For those the
plot and response are calculated for the internally stored frequencies.
Several Models
bode(m1,m2,...,mN) or bode(m1,m2,...mN,w) plots the Bode response of
several idmodel or idfrd models on a single figure. The models may be mixes
of different sizes and continuous/discrete. The sorting of the plots is made
based on the InputNames and OutputNames. If the frequencies w are specified,
these will apply to all non-idfrd models in the list. If you want different
frequencies for different models, you should thus first convert them to idfrd
objects using the idfrd command.
bode(m1,'PlotStyle1',...,mN,'PlotStyleN') further specifies which color,
Arguments
The output argument w contains the frequencies for which the response is
given, whether specified among the input arguments or not. The output
arguments mag and phase are 3-D arrays with dimensions
(number of outputs)x(number of inputs)x(length of w)
For SISO systems mag(1,1,k) and phase(1,1,k) give the magnitude and
phase (in degrees) at the frequency k = w(k). To obtain the result as a normal
vector of responses use mag = mag(:) and phase = phase(:).
For MIMO systems mag(i,j,k) is the magnitude of the frequency response at
frequency w(k) from input j to output i, and similairly for phase(i,j,k).
4-34
bode
If sdmag and sdphase are specified, the standard deviations of the magnitude
and phase are also computed. Then sdmag is an array of the same size as mag,
containing the estimated standard deviations of the response, and analogously
for sdphase.
See Also
4-35
compare
Purpose
4compare
Syntax
compare(data,m);
compare(data,m,k,sampnr,init)
compare(data,m1,m2,...,mN,Yplots)
compare(data,m1,'PlotStyle1',...,mN,'PlotStyleN',k,sampnr,init)
[yh,fit] =
compare(data,m1,'PlotStyle1',...,mN,'PlotStyleN',k,sampnr,init)
Description
4-36
compare
init = x0, where x0 is a column vector of the same size as the state vector of
the models, uses x0 as the initial state.
init = 'e' is default.
If data contains several experiments, separate plots are given for the different
experiments. In this case sampnr, if specified, must be a cell array with as many
entries as there are experiments.
Arguments
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);
See Also
sim, predict
4-37
covf
Purpose
4covf
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.
y(t)
u(t)
1
R ( i + ( j 1 )nz ,k + 1 ) = ---N
zi ( t )zj ( t + k )
(k)
= R
ij
t=1
where z j is the j-th row of z and missing values in the sum are replaced by
zero.
The optional argument maxsize controls the memory size as explained under
Algorithm Properties.
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
are handled. The expectation symbol E corresponds to the sample means.
Algorithm
See Also
spa
4-38
cra
Purpose
4cra
Syntax
cra(data);
[ir,R,cl] = cra(data,M,na,plot);
cra(R);
Description
4-39
cra
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')
See Also
4-40
impulse, step
c2d
Purpose
4c2d
Syntax
md = c2d(mc,T)
md = c2d(mc,T,method)
Description
Examples
Define a continuous-time system and study the poles and zeros of the sampled
counterpart.
mc = idpoly(1,1,1,1,[1 1 0],'Ts',0);
md = c2d(mc,0.5);
pzmap(md)
See Also
d2c
4-41
detrend
Purpose
4detrend
Syntax
zd = detrend(z)
zd = detrend(z,o,brkp)
Description
The default (o = 0) removes the zero-th order trends, i.e., the sample means
are subtracted.
With o = 1, linear trends are removed, after a least-squares fit. With brkp not
specified, one single line is subtracted from the entire data record. A
continuous piecewise linear trend is subtracted if brkp contains breakpoints at
sample numbers given in a row vector.
Note that detrend for iddata objects differs somewhat from detrend in the
Signal Processing Toolbox.
Examples
Remove a V-shaped trend from the output with its peak at sample number 119,
and remove the sample mean from the input.
zd(:,1,[]) = detrend(z(:,1,[]),1,119);
zd(:,[],1) = detrend(z(:,[],1));
4-42
d2c
Purpose
4d2c
Syntax
mc = d2c(md)
mc = d2c(md,'CovarianceMatrix',cov,'InputDelay',inpd)
Description
Note The transformation from discrete to continuous time is not unique. d2c
selects the continuous-time counterpart with the slowest time constants,
consistent with the discrete-time model. The lack of uniqueness also means
that the transformation may be ill-conditioned or even singular. In particular,
poles on the negative real axis, in the origin, or in the point 1, are likely to
cause problems. Interpret the results with care.
4-43
d2c
Examples
Note that the transformation to continuous time can be included in the n4sid
command by specifying the model to be continuos time.
mc = n4sid(data,3,'Ts',0)
See Also
c2d, nuderst
References
See Discrete and Continuous Time Models and Spectrum Normalization and
the Sampling Interval in the Tutorial.
4-44
EstimationInfo
Purpose
4EstimationInfo
Syntax
m.EstimationInfo
m.es
m.es.DataLength, etc
Description
4-45
EstimationInfo
vector' or 'Near (local) minimum'. The latter means that the expected
improvement is less than Tolerance (see Algorithm Properties).
See Also
4-46
idfrd
etfe
Purpose
4etfe
Syntax
g = etfe(data)
g = etfe(data,M,N)
Description
y ( t ) = G ( q )u ( t ) + v ( t )
data contains the output-input data and is an iddata object.
g is given as an idfrd object with the estimate of G ( e
) at the frequencies
w = [1:N]/Npi/T
Examples
4-47
etfe
Generate a periodic input, simulate a system with it, and compare the
frequency response of the estimated model with the true system at the excited
frequency points.
m = idpoly([1 -1.5 0.7],[0 1 0.5]);
u = iddata([],idinput([50,1,10],'sine'));
u.Period = 50;
y = sim(m,u);
me = etfe([y u])
bode(me,'b*',m)
Algorithm
The empirical transfer function estimate is computed as the ratio of the output
Fourier transform to the input Fourier transform, using fft. The periodogram
is computed as the normalized absolute square of the Fourier transform of the
time series.
The smoothed versions (M less than the length of z) are obtained by applying a
Hamming window to the output fast Fourier transform (FFT) times the
conjugate of the input FFT, and to the absolute square of the input FFT,
respectively, and subsequently forming the ratio of the results. The length of
this Hamming window is equal to the number of data points in z divided by M,
plus one.
See Also
4-48
spa
ffplot
Purpose
4ffplot
Syntax
ffplot(m)
[mag,phase,w] = ffplot(m)
[mag,phase,w,sdmag,sdphase] = ffplot(m)
ffplot(m1,m2,m3,...,w)
ffplot(m1,'PlotStyle1',m2,'PlotStyle2',...)
ffplot(m1,m2,m3,..'sd',sd,'mode',mode,'ap',ap)
Description
This function has exactly the same syntax as bode. The only difference is that
it gives graphs with linear frequency scales and Hz as the frequency unit.
See Also
bode, nyquist
4-49
freqresp
Purpose
4freqresp
Syntax
H = freqresp(m)
[H,w,covH] = freqresp(m,w)
Description
and in radians/second.
If m has ny outputs and nu inputs, and w contains Nw frequencies, the output
H is a ny-by-nu-by-Nw array such that H(:,:,k) gives the complex valued
response at the frequency w(k).
For a SISO model, H(:) to obtain a vector of the frequency response.
If w is not specified, a default choice is made. For a discrete-time model w will
be 128 equally spaced frequency points from 0 (excluded) to the Nyquist
frequency. For a continuous-time model, the default is
w=logspace(log10(pi/abs(Ts)/100),log10(10*pi/abs(Ts)),128)'
where Ts is the sampling interval of the data from which the model was
estimated. If the model is not estimated, Ts is taken as 1, which may make it
necessary to specify w as in input argument in this case.
[H,w,covH] = freqresp(M,w)
also returns the frequencies w and the covariance covH of the response. covH is
a 5-D array where covH(ky,ku,k,:,:) is the 2-by-2 covariance matrix of the
response from input ku to output ky at frequency w(k). The 1,1 element is the
variance of the real part, the 2,2 element the variance of the imaginary part
and the 1,2 and 2,1 elements the covariance between the real and imaginary
parts. squeeze(covH(ky,ku,k,:,:)) gives the covariance matrix of the
corresponding response.
If m is a time series (no input channels), H is returned as the (power) spectrum
of the outputs; an ny-by-ny-by-Nw array. Hence H(:,:,k) is the spectrum
matrix at frequency w(k). The element H(k1,k2,k) is the cross spectrum
between outputs k1 and k2 at frequency w(k). When k1=k2, this is the
real-valued power spectrum of output k1.
4-50
freqresp
See Also
4-51
fpe
Purpose
4fpe
Syntax
am = fpe(Model)
Description
1+dN
FPE = V ---------------------1dN
where V is the loss function, d is the number of estimated parameters, and N
is the number of estimation data.
See Also
EstimationInfo, aic
Reference
4-52
get
Purpose
4get
Syntax
Value = get(m,'PropertyName')
get(m)
Struct = get(m)
Description
accessed directly).
Struct = get(m) converts the object m into a standard MATLAB structure
with the property names as field names and the property values as field values.
Remark
See Also
4-53
getexp
4getexp
Purpose
Syntax
Description
See merge (iddata) and iddata for how to create multi-experiment data
objects.
The experiments can also be retrieved by a fourth subscript as in
d1 = data(:,:,:,ExperimentNumber). See help iddata/subsref for details
on this.
4-54
idarx
4idarx
Purpose
Syntax
m = idarx(A,B,Ts)
m = idarx(A,B,Ts,'Property1',Value1,...,,'PropertyN',ValueN)
Description
y ( t ) + A 1 y ( t 1 ) + A 2 y ( t 2 ) + + A na y ( t na ) =
B 0 u ( t ) + B 1 u ( t 1 ) + + B nb u ( t nb ) + e ( t )
Here A k and B k are matrices of dimensions ny-by-ny and ny-by-nu, respectively
(ny is the number of outputs, i.e., the dimension of the vector y ( t ) and nu is the
number of inputs). See Multivariable ARX Models: The idarx Model in the
Tutorial chapter.
The arguments A and B are 3-D arrays that contain the A matrices and the B
matrices of the model in the following way.
A is an ny-by-ny-by-(na+1) array such that
A(:,:,k+1) = Ak
A(:,:,1) = eye(ny)
Note that A always starts with the identity matrix, and that delays in the model
are defined by setting the corresponding leading entries in B to zero. For a
multivariate time series take B = [].
The optional property NoiseVariance sets the covariance matrix of the driving
noise source e ( t ) in the model above. The default value is the identity matrix.
4-55
idarx
idarx
Properties
4-56
idarx
Examples
Simulate a second order ARX model with one input and two outputs, and then
estimate a model using the simulated data.
A = zeros(2,2,3);
B = zeros(2,1,3)
A(:,:,1) =eye(2);
A(:,:,2) = [-1.5 0.1;-0.2 1.5];
A(:,:,3) = [0.7 -0.3;0.1 0.7];
B(:,:,2) = [1;-1];
B(:,:,3) = [0.5;1.2];
m0 = idarx(A,B,1);
u = iddata([],idinput(300));
e = iddata([],randn(300,2));
y = sim(m0,[u e]);
m = arx([y u],[[2 2;2 2],[2;2],[1;1]]);
See Also
4-57
iddata
Purpose
4iddata
Syntax
data = iddata(y,u)
data = iddata(y,u,Ts,'Property1',Value1,...,'PropertyN',ValueN)
Description
iddata is the basic object for dealing with signals in the toolbox. It is used by
most of the commands.
Basic Use
Let y be a column vector or an N-by-ny matrix. The columns of y correspond to
the different output channels. Similarly, u is a column vector or an N-by-nu
matrix containing the signals of the input channels.
data = iddata(y,u,Ts)
creates an iddata object containing these output and input channels. Ts is the
sampling interval. This construction is sufficient for most purposes.
The data is then plotted by plot(data) (see plot), and portions of the data
record are selected as in ze = data(1:300) or zv = data(501:700).
The signals in the output channels are retrieved by data.OutputData, or for
short data.y. Similarly the input signals are obtained by data.InputData or
data.u.
For a time series (no input channels) use data = iddata(y), or let u = [].
An iddata object can also contain just an input, by letting y = [].
The sampling interval can be changed by set(data,'Ts',0.3) or, simpler, by
data.Ts = 0.3.
The input and output channels are given default names like 'y1', 'y2',
'u1','u2', etc. The channel names can be set by
set(data,'InputName',{'Voltage','Current'},'OutputName','Tempera
ture')
(two inputs and one output is this example) and these names will then follow
the object and appear in all plots. The names will also be inherited by models
that are estimated from the data.
4-58
iddata
The actual time point values are given by the property 'SamplingInstants',
as in
plot(data.sa,data.u)
for a plot of the input with correct time points. Autofill is used for all properties,
and they are case insensitive.
Manipulating Channels
An easy way to set and retrieve channel properties is to use subscripting. The
subscripts are defined as
data(Samples,Outputs,Inputs)
so dat(:,3,:) is the data object obtained from dat by keeping all input
channels, but only output channel 3. (Trailing :s can be omitted so
dat(:,3,:)= dat(:,3).).
The channels can also be retrieved by their names, so that
dat(:,{'speed','flow'},[])
is the data object where the indicated output channels have been selected and
no input channels are selected.
Moreover
dat1(101:200,[3 4],[1 3]) = dat2(1001:1100,[1 2],[6 7])
will change samples 101 to 200 of output channels 3 and 4 and input channels
1 and 3 in the iddata object dat1 to the indicated values from iddata object
dat2. The names and units of these channels will then also be changed
accordingly.
To add new channels, use horizontal concatenation of iddata objects.
dat =[dat1, dat2];
4-59
iddata
(see Horizontal Concatenation below) or add the data record directly. Thus
dat.u(:,5) = U
Nonequal Sampling
The property 'SamplingInstants' gives the sampling instants of the data
points. It can always be retrieved by get(DAT,'SamplingInstants') (or dat.s)
and is then computed from dat.Ts and dat.Tstart. 'SamplingInstants' can
also be set to an arbitrary vector of the same length as the data, so that
nonequal sampling can be handled. Ts is then automatically set to [ ]. Most of
the estimation routines, though, do not handle unequally sampled data.
Multiple Experiments
The iddata object can also store data from separate experiments. The property
'ExperimentName' is used to separate the experiments. The number of data as
well as the sampling properties can vary from experiment to experiment, but
the input and output channels must be the same. (Use NaN to fill possibly
unmeasured channels in certain experiments.) The data records will be cell
arrays, where the cells contain data from each experiment.
Multiple experiments can be defined directly by letting the 'y' and 'u'
properties as well as 'Ts' and 'Tstart' be cell arrays.
It is normally easier to create multi-experiment data by merging experiments
as in
dat = merge(dat1,dat2)
See the reference page for merge (data). Storing multiple experiments as one
iddata object may be very useful to handle experimental data that has been
collected on different occasions, or when a data set has been split up to remove
bad portions of the data. All the toolbox s routines accept multiple
experiment data.
Experiments can be retrieved by the command getexp. They can also be
retrieved by subscripting with a fourth index: dat(:,:,:,3) is experiment
number 3 and dat(:,:,:,{'Day1','Day4'}) retrieves the two experiments
with the indicated names.
4-60
iddata
which adds the data in dat2 as a new experiment with name 'Run4'. See
iddemo number 8 for an illustration of how multiple experiments can be used.
iddata
Properties
In the list below, N denotes the number of samples of the signals, ny the number
of output channels, nu the number of input channels, and Ne the number of
experiments:
Domain: Assumes the values 'Time' or 'Frequency' and denotes whether the
data are time domain or frequency domain data.
Name: An optional name for the data set. An arbitrary string.
OutputData, InputData: The data matrices y and u. In the single
experiment case y is an N-by-ny matrix and u is an N-by-nu matrix. For
multiple experiments y and u are 1-by-Ne cell arrays, with each cell
containing the data for the different experiments.
OutputName, InputName: Cell arrays of length ny-by-1 and nu-by-1
containing the names of the output and input channels. If not specified,
default names, {'y1';'y2';...} and {'u1';'u2';...} are given.
OutputUnit, InputUnit: Cell arrays of length ny-by-1 and nu-by-1
containing the units of the output and input channels.
TimeUnit: The unit for the sampling instants.
Ts: Sampling interval. A positive scalar. For multiexperiment data, Ts is a
1-by-Ne cell array, with each cell containing the sampling interval of the
corresponding experiment. For nonequally sampled data, Ts = [].
Tstart: The starting time of the data record. A scalar. For multiexperiment
data, Tstart is a 1-by-Ne cell array, with each cell containing the starting
time for the corresponding experiment
SamplingInstants: The time values of the sample points. A N-by-1 vector.
For multiple experiment data, SamplingInstants is a 1-by-Ne cell array,
with each cell containing the sampling instants of the corresponding
experiment. For equally sampled data, SamplingInstants is generated from
Ts and Tstart.
4-61
iddata
Period: The period of the input. A nu-by-1 vector, where the k-th entry
contains the period of the k-th input. Period = inf means nonperiodic data.
For multiexperiment data, Period is a 1-by-Ne cell array with each cell
containing the period(s) for the input of the corresponding experiment.
InterSample: Describes the intersample behavior of the input channels. An
nu-by-1 cell array where the (k,1) element is 'zoh', 'foh', or 'bl',
denoting that input number k is piecewise constant, piecewise linear, or
band limited. For multiple experiment data, InterSample is an nu-by-Ne cell
array.
ExperimentName: A string containing the name of the experiment. For
multiple experiment data ExperimentName is a 1-by-Ne cell array with each
cell containing the name of the corresponding experiment. It can be freely
set, and is by default given names {'Exp1', 'Exp2',...}.
Notes: An arbitrary field to store extra information and notes about the
object.
UserData: An arbitrary field for any possible use.
Note that all properties can be set or retrieved either by set/get or by
subscripts. Autofill applies to all properties and values, and are case
insensitive:. 'y' and 'u' can be used as short for 'OutputData' and
'InputData'. 'y' and 'u' can also replace 'Output' and 'Input' in the other
properties.
data.y=randn(100,2)
data.una = 'Voltage'
set(data,'tim','minute')
p = data.per
For a complete list of property values, use get(data). To see possible value
assignments, use set(data).
Subreferencing
Use a colon (:) to denote all samples/channels and the empty matrix ([ ]) to
denote no samples/channels. The channels can be referenced by number or by
name. For several names, a cell array must be used.
4-62
iddata
dat2 = dat(:,'y3',{'u1','u4'})
dat2 = dat(:,3,[1 4])
will select the samples with time marks between 1.27 and 9.3.
Subreferencing with curly brackets refers to the experiment.
data{Experiment}(samples,outputs,inputs)
Horizontal
Concatenation
dat = [dat1,dat2,...,datN]
creates an iddata object dat, consisting of the input and output channels in
dat1,... datN. Default channel names ('u1', 'u2', 'y1', 'y2', etc.) will
be changed so that overlaps in names are avoided, and the new channels will
be added.
If datk contains channels with user specified names, that are already present
in the channels of Datj, j<k, these new channels will be ignored.
Vertical
Concatenation
creates an iddata object dat whose signals are obtained by stacking those of
datk on top of each other. That is
dat.OutputData = [dat1.Ouputdata;dat2.OutputData; ...
datN.OutputData]
and similarly for the inputs. The datk objects must all have the same number
of channels and experiments.
Online Help
Functions
See Also
4-63
ident
Purpose
4ident
Syntax
ident
ident(session,directory)
Description
When the session is specified, the interface will open with this session active.
Typing ident(session,directory) on the MATLAB command line, when the
interface is active, will load and open the session in question.
For more information about the graphical user interface, see Chapter 2, The
Graphical User Interface.
Examples
4-64
ident('iddata1.sid')
ident('mydata.sid','\matlab\data\cdplayer\')
idfilt
Purpose
4idfilt
Syntax
zf = idfilt(z,filter)
zf = idfilt(z,ord,Wn)
zf = idfilt(z,ord,causality)
[zf,mf] = idfilt(z,ord,Wn,hs)
Description
As an explicit argument filter. This in turn can be given either as any SISO
idmodel or LTI model object, or as a cell array {A,B,C,D} of SISO state-space
matrices or as a cell array {num,den} of numerator/denominator filter
coefficients.
As a triplet of arguments ...,ord,Wn,hs,..., which defines a Butterworth
filter of order ord. If hs is not specified and Wn contains just one element, a
low pass filter with cutoff frequency Wn (measured as a fraction of the
Nyquist frequency) is obtained. If hs =' high' a high pass filter with this
cutoff frequency is obtained instead. If Wn = [Wnl Wnh] is a vector with two
elements, a filter (of order 2*ord) with passband between Wnl and Wnh is
obtained is hs is not specified. If hs = 'stop' a bandstop filter with stop band
between these two frequencies is obtained instead.
The output argument mf is the filter given as an idmodel object.
With causality = 'causal' (default) causal filtering is used. With
causality = 'noncausal', a noncausal, zero-phase filter is used for the
filtering.
It is common practice in identification to select a frequency band where the fit
between model and data is concentrated. Often this corresponds to bandpass
filtering with a pass band over the interesting breakpoints in a Bode diagram.
For identification where a disturbance model also is estimated, it is better to
achieve the desired estimation result by using the property 'Focus' (see
Algorithm Properties) than just to prefilter the data.
4-65
idfilt
Algorithm
The Butterworth filter is the same as butter in the Signal Processing Toolbox.
Also, the zero-phase filter is equivalent to filtfilt in that toolbox.
References
4-66
idfrd
Purpose
Create the idfrd (Identified Frequency Response Data) object that stores
frequency function and spectrum data along with covariance information.
Syntax
h = idfrd(Response,Freqs,Ts)
h = idfrd(Response,Freqs,Ts,'CovarianceData',Covariance, ...
'SpectrumData',Spec,'NoiseCovariance',Speccov,'property1', ...
Value1,'PropertyN',ValueN)
h = idfrd(mod)
h = idfrd(mod,Freqs)
Description
4idfrd
For a model
y ( t ) = G ( q )u ( t ) + H ( q )e ( t )
stores the transfer function estimate G (see equation (Equation 3-4) in the
Tutorial chapter)
i
G(e )
as well as the spectrum of the additive noise ( v ) at the outputs
v ( ) = T H ( e
iT 2
4-67
idfrd
element is the variance of the imaginary part and the 1-2 and 2-1 elements is
the covariance between the real and imaginary parts.
squeeze(Covariance(ky,ku,kf,:,:)) thus gives the covariance matrix of the
corresponding response.
The information about spectrum is optional. The format is as follows:
spec is a 3-D array of dimension ny-by-ny-by-Nf, such that spec(ky1,ky2,kf)
is the cross spectrum between the noise at output ky1 and the noise at output
ky2, at frequency Freqs(kf). When ky1=ky2 the (power) spectrum of the noise
at output ky1 is thus obtained. For a single output model, spec can be given as
a vector.
speccov is a 3-D array of dimension ny-by-ny-by-Nf, such that
speccov(ky1,ky1,kf) is the variance of the corresponding power spectrum.
where Ts is the sampling interval specified by mod and for the continuous-time
case
Freqs = logspace(log10(pi/Ts/100),log10(10*pi/ Ts),128)
where Ts is the sampling interval of the data from which the model was
estimated. If the model is not estimated, a simple default choice of Freqs is
made. In this case it may be necessary to supply the argument Freqs explicitly.
If mod has InputDelay different from zero, these are appended as phase lags,
and h will then have an InputDelay 0.
4-68
idfrd
idfrd
Properties
4-69
idfrd
InputUnit: The units in which the input channels are measured. It has the
same format as 'InputName'.
OutputUnit: Correspondingly for the output channels.
InputDelay: A row vector of length equal to the number of input channels.
Contains the delays from the input channels. These should thus be appended
as phase lags when the response is calculated. This is done automatically by
freqresp, bode, ffplot, and nyquist. Note that if the idfrd is calculated
form an idmodel, possible input delays in that model are converted to phase
lags, and InputDelay of the idfrd model is set to zero.
Notes: An arbitrary field to store extra information and notes about the
object.
UserData: An arbitrary field for any possible use.
EstimationInfo: A structure that contains information about the estimation
process that is behind the frequency data. It contains the following fields:
- Status: Gives the status of the model, e.g., 'Not estimated'.
- Method: The identification routine that created the model.
- WindowSize: If the model was estimated by spa or etfe, the size of window
(input argument M) that was used.
- DataName: The name of the data set from which the model was estimated.
- DataLength: The length of this data set.
Note that all properties can be set or retrieved either by set/get or by
subscripts. Autofill applies to all properties and values, and these are case
insensitive:
h.ts = 0
loglog(h.fre,squeeze(h.spe(2,2,:)))
For a complete list of property values, use get(m). To see possible value
assignments, use set(m). See also idprops idfrd.
4-70
idfrd
Subreferencing
contains the information for measured inputs only, that is, just ResponseData,
while
h('n')
Horizontal
Concatenation
Vertical
Concatenation
4-71
idfrd
Examples
Compute separate idfrd models, one containing g and the other the noise
spectrum.
g = idfrd(m('m'))
phi = idfrd(m('n'))
See Also
4-72
idgrey
4idgrey
Purpose
Syntax
m = idgrey(MfileName,ParameterVector,CDmfile)
m = idgrey(MfileName,ParameterVector,CDmfile,FileArgument,Ts,...
'Property1',Value1,...,'PropertyN',ValueN)
Description
4-73
idgrey
There are also two properties, DisturbanceModel and InitialState that can
be used to affect the parameterizations of K and X0, thus overriding the
outputs from the M-file. See below.
idgrey
Properties
4-74
idgrey
For a complete list of property values, use get(m). To see possible value
assignments, use set(m). See also idprops and idgrey.
M-File Details
4-75
idgrey
x ( t ) = A ( )x ( t ) + B ( )u ( t ) + K ( )e ( t )
x( 0 ) = x0 ( )
y ( t ) = C ( )x ( t ) + D ( )u ( t ) + e ( t )
Here x ( t ) is the time derivative x ( t ) for a continuous time model and x ( t + Ts )
for a discrete time model.
The matrices in this time-discrete model can be parameterized in an arbitrary
way by the vector . Write the format for the M-file as follows.
[A,B,C,D,K,x0] = mymfile(pars,T,Auxarg)
Here the vector pars contains the parameters , and the output arguments A,
B, C, D, K, and x0 are the matrices in the model description that correspond to
this value of the parameters and this value of the sampling interval T.
T is the sampling interval, and Auxarg is any variable of auxiliary variables
with which you want to work. (In that way you can change certain constants
and other aspects in the model structure without having to edit the M-file.)
Note that the two arguments T and Auxarg must be included in the function
head of the M-file, even if they are not utilized within the M-file.
4-76
idgrey
Examples
4-77
idinput
Purpose
4idinput
Syntax
u = idinput(N)
u = idinput(N,type,band,levels)
[u,freqs] = idinput(N,'sine',band,levels,sinedata)
Description
idinput generates input signals of different kinds, that are typically used for
that determine the lower and upper bound of the pass-band. The frequencies
wlow and whigh are expressed in fractions of the Nyquist frequency. A white
noise character input is thus obtained for band = [0 1], which also is the
default value.
4-78
idinput
where B is such that the signal is constant over intervals of length 1/B (the
clock period). Also in this case the default is band = [0 1].
The argument levels defines the input level. It is a row vector
levels = [minu, maxu]
such that the signal u will always be between the values minu and maxu for the
choices type = 'rbs', 'prbs' and 'sine'. For type = 'rgs', the signal level is such
that minu is the mean value of the signal, minus one standard deviation, while
maxu is the mean value plus one standard deviation. Gaussian white noise with
zero mean and variance one is thus obtained for levels = [-1, 1], which is
also the default value.
In the 'sine' case, the sinusoids are chosen from the frequency grid
freq = 2*pi*[1:Grid_Skip:fix(P/2)]/P intersected with pi*[band(1)
band(2)]
(for Grid_Skip, see below.) For multi-input signals, the different inputs use
different frequencies from this grid. An integer number of full periods is always
delivered. The selected frequencies are obtained as the second output
argument, freqs, where row ku of freqs contains the frequencies of input
number ku. The resulting signal is affected by a fifth input argument sinedata
sinedata = [No_of_Sinusoids, No_of_Trials, Grid_Skip]
Algorithm
Very simple algorithms are used. The frequency contents is achieved for 'rgs'
by an eighth order Butterworth, noncausal filter, using idfilt. This is quite
reliable. The same filter is used for the 'rbs' case, before making the signal
4-79
idinput
binary. This means that the frequency contents is not guaranteed to be precise
in this case.
For the 'sine' case, the frequencies are selected to be equally spread over the
chosen grid, and each sinusoid is given a random phase. A number of trials are
made, and the phases that give the smallest signal amplitude are selected. The
amplitude is then scaled so as to satisfy the specifications of levels.
Reference
4-80
See Sderstrm and Stoica (1989), Chapter C5.3. For a general discussion of
input signals, see Ljung (1999), Section 13.3.
idmodel
Purpose
4idmodel
Description
idmodel is an object that the user does not deal with directly. It contains all the
common properties of the model objects idarx, idgrey, idpoly, and idss,
Basic Use
If you just estimate models from data, the model objects should be transparent.
All parametric estimation routines return idmodel results.
m = arx(Data,[2 2 1])
The model m contains all relevant information. Just typing m will give a brief
account of the model. present(m) also gives information about the
uncertainties of the estimated parameters. get(m) gives a complete list of
model properties.
Most of the interesting properties can be directly accessed by subreferencing.
m.a
m.da
See the property list obtained by get(m), as well as the property lists of
idgrey, idarx, idpoly, and idss in the Command Reference for more details
on this. See also idprops.
The characteristics of the model m can be directly examined and displayed by
commands like impulse, step, bode, nyquist, pzmap. The quality of the model
is assessed by commands like compare, and resid. If you have the Control
System Toolbox, just typing view(m) gives access to various display functions.
To extract state-space matrices, transfer function polynomials, etc., use the
commands
arxdata, polydata, tfdata, ssdata, zpkdata
and by idfrd and freqresp, the frequency response of the model can be
computed.
4-81
idmodel
Use colon (:) to denote all channels and the empty matrix ([ ]) to denote no
channels. The channels can be referenced by number or by name. For several
names, a cell array must be used.
m3 = m('position',{'power','speed'})
or
m3 = m(3,[1 4])
Thus m3 is the model obtained from m by looking at the transfer functions from
input numbers 1 and 4 (with input names 'power' and 'speed') to output
number 3 (with name position).
For a single output model m
m4 = m(inputs)
will select the corresponding input channels, and for a single input model
m5 = m(outputs)
4-82
idmodel
(4-1)
creates a time-series model m2 from m by ignoring the measured input. That is,
m2 describes the signal He.
For a system with measured inputs, bode, step, and many other
transformation and display functions just deal with the transfer function
matrix G. To obtain or graph the properties of the disturbance model H, it is
therefore important to make the transformations m('n'). For example,
bode(m('n'))
will plot the additive noise spectra according to the model m, while
4-83
idmodel
bode(m)
This creates a model m3 with all input channels, both measured u and noise
sources e, being treated as measured signals,. That is, m3 is a model from u and
e to y, describing the transfer functions G and H. The information about the
variance of the innovations e is then lost. For example, studying the step
response from the noise channels, will then not take into consideration how
large the noise contributions actually are.
To include that information, e should first be normalized e = Lv , so that v
becomes white noise with an identity covariance matrix.
m4 = noisecnv(m,'Norm')
ny inputs)
4-84
idmodel
returns the properties of the transfer function HL (ny outputs and ny inputs).
If m is a time series model, fcn(m) returns the properties of H, while
fcn(noisecnv(m,'Norm'))
idmodel
Properties
In the list below, ny is the number of output channels, and nu is the number of
input channels:
Name: An optional name for the data set. An arbitrary string.
OutputName, InputName: Cell arrays of length ny-by-1 and nu-by-1
containing the names of the output and input channels. For estimated
models, these are inherited from the data. If not specified, they will be given
default names: {'y1','y2',...} and {'u1','u2',...}.
OutputUnit, InputUnit: Cell arrays of length ny-by-1 and nu-by-1 containing
the units of the output and input channels. Inherited from data for estimated
models.
TimeUnit: The unit for the sampling interval.
Ts: Sampling interval. A non-negative scalar. Ts = 0 denotes a
continuous-time model. Note that changing just Ts will not recompute the
model parameters. Use c2d and d2c for recomputing the model to other
sampling intervals.
ParameterVector: The vector of adjustable parameters in the model
structure. Initial/nominal values or estimated values, depending on the
status of the model. A column vector.
PName: The names of the parameters. A cell array of the length of the
parameter vector. If not specified, it will contain empty strings. See also
setpname.
4-85
idmodel
4-86
idmodel
For a complete list of property values, use get(m). To see possible value
assignments, use set(m).
Subreferencing
Use colon (:) to denote all channels and the empty matrix ([ ]) to denote no
channels. The channels can be referenced by number or by name. For several
names, a cell array must be used.
m2 = m('y3',{'u1','u4'})
m3 = m(3,[1 4])
will select the corresponding input channels, and for a single input model
m5 = m(outputs)
Similarily the string 'noise' (or any abbreviation) refers to the noise input
channels. See above under The Noise Channels for more details.
Horizontal
Concatenation
creates an idmodel object m, consisting of all the input channels in m1,... mN.
The output channels of mk must be the same.
Vertical
Concatenation
4-87
idmodel
creates an idmodel object m consisting of all the output channels in m1, m2,
..mN. The input channels of mk must all be the same.
Online Help
Functions
See Also
4-88
idmodred
Purpose
4idmodred
Syntax
MRED = idmodred(M)
MRED = idmodred(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 is an idss model.
The function requires several routines in the Control System Toolbox.
ORDER: The desired order (dimension of the state-space representation). If
ORDER = [], which is the default, a plot will show how the diagonal elements
Algorithm
The functions balreal and modred from the Control System Toolbox are used.
The plot, in case ORDER = [], shows the vector g as returned from 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,[4ones(3,3),4ones(3,2),ones(3,2)]);
MRED = idmodred(M,3);
bode(M,MRED)
Use the reduced order model as initial condition for a third order state-space
model.
M2= pem(data,MRED);
4-89
idpoly
Purpose
Syntax
4idpoly
Description
idpoly creates a model object containing parameters that describe the general
multi-input-single-output model structure.
B1 ( q )
B nu ( q )
C(q)
A ( q )y ( t ) = --------------- u 1 ( t nk 1 ) + ------------------- u nu ( t nk nu ) + ------------- e ( t )
D(q)
F1 ( q )
F nu ( q )
A, B, C, D, and F specify the polynomial coefficients.
For single-input systems, these are all row vectors in the standard MATLAB
format.
A = [1 a1 a2 ...
ana]
consequently describes
A ( q ) = 1 + a1 q
+ + a na q
na
A, C, D, and F all start with 1, while B contains leading zeros to indicate the
delays. See Polynomial Representation of Transfer Functions in the
Tutorial chapter.
For multi-input systems, B and F are matrices with one row for each input.
For time series, B and F are entered as empty matrices.
B = [];
F = [];
4-90
idpoly
Properties
4-91
idpoly
For a complete list of property values, use get(m). To see possible value
assignments, use set(m). See also idprops idpoly.
Examples
4-92
idpoly
Note that the continuous time model will automatically be sampled to the
sampling interval of the data, when simulated, so the above is also achieved by
u = iddata([],[u1 u2],0.1)
y = sim(m,u)
See Also
sim, idss
References
See Ljung (1999) Section 4.2 for the model structure family.
See
T. Knudsen (1994): A new method for estimating ARMAX models. In Proc. 10th
IFAC Symposium on System Identification, pp 611-617. Copenhagen, Denmark
for the backcast method.
4-93
idss
Purpose
4idss
Syntax
m = idss(A,B,C,D)
m = idss(A,B,C,D,K,x0,Ts,'Property1',Value1,...,'PropertyN',ValueN)
mss = idss(m1)
Description
4-94
idss
Parameterizations
There are several different ways to define the parameterization of the
state-space matrices. The parameterization will determine which parameters
can be adjusted to data by the parameter estimation routine pem:
Free Black-Box parameterizations: This is the default situation and
corresponds to letting all parameters in A, B, and C be freely adjustable. This
is obtained by setting the property 'SSParameterization' = 'Free'. The
parameterizations of D, K, and X0 are then determined by the following
properties:
- 'nk': A row vector of the same length as the number of inputs. The ku-th
element is the delay from input channel no ku. nk =[0,...,0], thus means
that there is no delay from any of the inputs, and that consequently all
elements of the D matrix should be estimated. nk =[1,...,1] means that
there is a delay of 1 from each input, so that the D matrix is fixed to be zero.
- 'DisturbanceModel': This property affects the parametrization of K and
can assume the values:
'Estimate' which means that all elements of the K matrix are to be
estimated.
'None': all elements of K are fixed to zero.
'Fixed': all elements of K are fixed to their nominal/initial values.
4-95
idss
idss Properties
4-96
idss
- 'Free': Means that all parameters in A,B and C are freely adjustable and
the parameterizations of D, K and X0 depend on the properties 'nk',
'DisturbanceModel' and 'InitialState'
- 'Canonical': Means that A and C are parameterized as an observer
canonical form. The details of this parameterization depends on the
property 'CanonicalIndices'. The B-matrix is always fully
parameterized, and the parameterizations of D, K, and X0 depend on the
properties 'nk', 'DisturbanceModel', and 'InitialState'.
- 'Structured': Means that the parametrization is determined by the
properties (the structure matrices) 'As', 'Bs', 'Cs', 'Ds', 'Ks', and
'X0s'. A NaN in any position in these matrices denotes a freely adjustable
parameter and a numeric value denotes a fixed and nonadjustable
parameter.
nk: A row vector with as many entries as the number of input channels. The
entry number k denotes the time delay from input number k to y(t). This
property is relevant only for 'Free' and 'Canonical' parameterizations. If
any delay is larger than 1, the structure of the A, B, and C matrices will accommodate this delay, at the price of a higher order model.
DisturbanceModel with possible values:
- 'Estimate': Means that the K matrix is fully parameterized.
- 'None': Means that the K matrix is fixed to zero. This gives a so-called
Output-Error model, since the model output depends on past inputs only.
- 'Fixed': Means that the K matrix is fixed to the current nominal values
InitialState with possible values:
- 'Estimate': Means that X0 is fully parameterized.
- 'Zero': Means that X0 is fixed to zero.
- 'Fixed': Means that X0 is fixed to the current nominal value.
- 'Backcast': The value of X0 is estimated by the identification routines as
the best fit to data, but it is not stored.
- 'Auto': Gives an automatic, and data-dependent choice between
'Estimate', 'Zero' and 'Backcast'.
A, B, C, D, K, and X0: The state-space matrices that can be set and
retrieved at any time. These contain both fixed values and
estimated/nominal values.
4-97
idss
dA, dB, dC, dD, dK, and dX0: The estimated standard deviations of the
state-space matrices. These cannot be set, only retrieved. Note that these are
not defined for an idss model with 'Free' SSParameterization. You can
then convert the parameterization to 'Canonical' and study the
uncertainties of the matrix elements in that form.
As, Bs, Cs, Ds, Ks, and X0s: These are the structure matrices that have
the same sizes as A, B, C etc. and show the freely adjustable parameters as
NaN in the corresponding position. These properties are used to define the
model structure for 'SSParameterization' = 'Structured'. They are
however always defined and can be studied also for the other
parametrizations.
CanonicalIndices: Determines the details of the canonical
parameterization. It is a row vector of integers with as many entries as there
are outputs. They sum up to the system order. This is the so-called
pseudo-canonical multi-index, with an exact definition, e.g., on page 132 in
Ljung (1999). A good default choice is 'Auto'. This property is relevant only
for the canonical parameterization case. Note however, that for 'Free'
parameterizations, the estimation algorithms also store a canonically
parameterized model, to handle the model uncertainty.
In addition to these properties, idss objects also have all the properties of the
idmodel object. See idmodel properties, Algorithm Properties, and
EstimationInfo.
Note that all properties can be set and retrieved either by set/get or by
subscripts. Autofill applies to all properties and values, and are case
insensitive.
m.ss='can'
set(m,'ini','z')
p = eig(m.a)
For a complete list of property values, use get(m). To see possible value
assignments, use set(m). See also idprops idss.
Examples
4-98
idss
1 0
x =
x+ 3 u
4
0 2
y = 1 1 x+e
4-99
idss
See Also
4-100
n4sid, pem
impulse
Purpose
Syntax
Description
4impulse
mixture.
For a discrete time idmodel m, the impulse response y and, when required, its
estimated standard deviation ysd, is computed using sim. When called with
output arguments, y, ysd and the time vector t are returned. When impulse is
called without output arguments, a plot of the impulse response is shown. If sd
is given a value larger than zero, a confidence region around zero is drawn. It
corresponds to the confidence of sd standard deviations. In the plots, the
impulse is inversely scaled with the sampling interval, so that it has the same
energy regardless of the sampling.
Adding an argument 'fill' among the input arguments gives an uncertainty
region marked by a filled area, rather than by dash-dotted lines.
The start time T1 and the end time T2 can be specified by Time= [T1 T2]. If T1
is not given, it is set to -T2/4. The negative time lags (the impulse is always
assumed to occur at time 0) show possible feedback effects in the data, when
the impulse is estimated directly from data. If Time is not specified, a default
value is used.
For an iddata set data, impulse(data) estimates a high order, noncausal FIR
model after first having prefiltered the data so that the input is as white as
possible. The impulse response of this FIR model and, when asked for, its
confidence region is then plotted. When called with an output argument,
impulse, in the iddata case, returns this FIR model, stored as an idarx
model.The order of the prewhitening filter can be specified as na. The default
value is na = 10.
4-101
impulse
Any number and any mixture of models and data sets can be used as input
arguments. The responses are plotted with each input/output channel (as
defined by the models and data sets InputName and OutputName) as a
separate plot. Colors, linestyles, and marks can be defined by PlotStyle
values. These are the same as for the regular plot command, like
impulse(m1,'b-*',m2,'y--',m3,'g')
The noise input channels in m are treated as follows: Consider a model m with
both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that m.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v:
y = Gu + HLv
cov ( v ) = I
impulse(m) plots the impulse response of the transfer function G.
impulse(m('n')) plots the impulse response of the transfer function H. (ny
inputs and ny outputs).The input channels have names e@yname, where
yname is the name of the corresponding output.
If m is a time series, that is nu = 0, impulse(m) plots the impulse response of
the transfer function H.
impulse(noisecnv(m)) plots the impulse response of the transfer function
[G H] (nu+ny inputs and ny outputs). The noise input channels have names
e@yname, where yname is the name of the corresponding output.
impulse(noisecnv(m,'norm')) plots the impulse response of the transfer
function [G HL] (nu+ny inputs and ny outputs. The noise input channels
have names v@yname, where yname is the name of the corresponding output.
Arguments
4-102
impulse
ysd has the same dimensions as y and contains the standard deviations of y.
If impulse is called with an output argument and a single data set in the input
arguments, the output is returned as an idarx model mod containing the high
order FIR model, and its uncertainty. By calling impulse with mod, the
responses can be displayed and returned without having to redo the
estimation.
Example
See Also
cra, step
4-103
init
Purpose
4init
Syntax
m = init(m0)
m = init(m0,R,pars,sp)
Description
See Also
4-104
ivar
Purpose
4ivar
Syntax
m = ivar(y,na)
m = ivar(y,na,nc,maxsize)
Description
Examples
See Also
References
4-105
ivstruc
Purpose
4ivstruc
Compute fit between simulated and measured output for a group of model
structures.
Syntax
v = ivstruc(ze,zv,NN)
v = ivstruc(ze,zv,NN,p,maxsize)
Description
with the same interpretation as described for arx. See struc for easy
generation of typical NN matrices for single-input systems.
Each of ze and zv are iddata objects containing output-input data. Models for
each model structure defined in NN are estimated using the instrumental
variable (IV) method on data set ze. The estimated models are simulated using
the inputs from data set zv. The normalized quadratic fit between the
simulated output and the measured output in zv is formed and returned in v.
The rows below the first row in v are the transpose of NN, and the last row
contains the logarithms of the condition numbers of the IV matrix
( t )
(t)
Note The IV method used does not guarantee that the obtained models are
stable. The output-error fit calculated in v may then be misleading.
4-106
ivstruc
Examples
Compare the effect of different orders and delays, using the same data set for
both the estimation and validation.
v = ivstruc(z,z,struc(1:3,1:2,2:4));
nn = selstruc(v)
m = iv4(z,nn);
Algorithm
See Also
4-107
ivx
Purpose
4ivx
Estimate the parameters of an ARX model using the instrumental variable (IV)
method with arbitrary instruments.
Syntax
m = ivx(data,orders,x)
m = ivx(data,orders,x,maxsize)
Description
ivx is a routine analogous to the iv4 routine, except that you can use arbitrary
instruments. These are contained in the matrix x. Make this the same size as
the output, data.y. In particular, if data contains several experiments x must
be a cell array with one matrix/vector for each experiment. The instruments
used are then analogous to the regression vector itself, except that y is replaced
by x.
Note that ivx does not return any estimated covariance matrix for m, since that
requires additional information. m is returned as an idpoly object for single
output systems and as an idarx object for multi-output systems.
Use iv4 as the basic IV routine for ARX model structures. The main interest
in ivx lies in its use for nonstandard situations; for example when there is
feedback present in the data, or when other instruments need to be tried out.
Note that there is also an IV version that automatically generates instruments
from certain filters you define (type help iv).
See Also
iv4, ivar
References
4-108
iv4
Purpose
4iv4
Syntax
m = iv4(data,orders)
m = iv4(data,'na',na,'nb',nb,'nk',nk)
m= iv4(data,orders,'Property1',Value1,...,'PropertyN',ValueN)
Description
This routine is an alternative to arx and the use of the arguments are entirely
analogous to the arx function. The main difference is that the procedure is not
sensitive to the color of the noise term e ( t ) in the model equation.
For an interpretation of the loss function (innovations covariance matrix),
consult Interpretation of the Loss Function in the Tutorial chapter.
Examples
Algorithm
The first stage uses the arx function. The resulting model generates the
instruments for a second-stage IV estimate. The residuals obtained from this
model are modeled as a high-order AR model. At the fourth stage, the
input-output data are filtered through this AR model and then subjected to the
IV function with the same instrument-filters as in the second stage.
For the multi-output case, optimal instruments are obtained only if the noise
sources at the different outputs have the same color. The estimates obtained
with the routine are reasonably accurate though even in other cases.
See Also
arx, oe
References
4-109
LTI Commands
4LTI Commands
Purpose
To allow direct calls to typical LTI commands from idmodel objects. The
Control System Toolbox is required for these commands.
Syntax
Description
If you have the Control System Toolbox, most of the relevant LTI-commands,
as listed under Syntax, can be directly applied to any idmodel (idarx, idgrey,
idpoly, idss). You can also use the overloaded operations +, -, and *. The same
operations are performed and the result is delivered as an idmodel. The
original covariance information is however lost, most of the time.
Example
You have two more or less identical processes connected in series. Estimate a
model for one of them, and use that to form an initial estimate for a model of
the connected process.
m = pem(data) % data concerns one of the processes
m2 = pem(data2,m*m) % data2 are from the whole connected process
4-110
merge (iddata)
Purpose
4merge (iddata)
Syntax
dat = merge(dat1,dat2,....,datN)
Description
dat collects the data sets in dat1,.. datN into one iddata object, with several
experiments. The number of experiments in dat will the sum of the number of
experiments in datk. For the merging to be allowed a number of conditions
must be satisfied:
All of datk must have the same number of input channels, and the
InputNames must be the same.
All of datk must have the same number of output channels, and the
OutputNames must be the same. If some input or output channel is lacking in
one experiment, it can be replaced by a vector of NaN s to conform with these
rules.
If the ExperimentNames of datk have been specified to something else than
the default 'Exp1', 'Exp2', etc., they must all be unique. If default names
overlap, they will be modified, so that dat will have a unique list of
ExperimentNames.
The sampling intervals, the number of observations, and the input properties
(Period, InterSample) may be different in the different experiments.
The individual experiments can be retrieved by the command getexp. The can
also be retrieved by subreferencing with a fourth index:.
dat1 = dat(:,:,:,ExperimentNumber) or dat1 =
dat(:,:,:,ExperimentName)
Example
Bad portions of data have been detected around sample 500 and between
samples 720 - 730. Cut out these bad portions and form a multiple experiment
data set that can be used to estimate models, without the bad data destroying
the estimate.
dat = merge(dat(1:498),dat(502:719),dat(719:1000))
4-111
merge (iddata)
m = pem(dat)
Use the first two parts to estimate the model and the third one for validation.
m = pem(dat{1:2});
compare(dat{3},m)
See Also
4-112
iddata, getexp
merge (idmodel)
Purpose
4merge (idmodel)
Syntax
m = merge(m1,m2,....,mN)
[m,tv] = merge(m1,m2)
Description
The models m1,m2,...,mN must all be of the same structure, just differing in
parameter values and covariance matrices. m is then the merged model, where
the parameter vector is a statistically weighted mean (using the covariance
matrices to determine the weights) of the parameters of mk.
When two models are merged,
[m, tv] = merge(m1,m2)
2
and
mb = arx(merge(z1,z2),[2 3 4]);
lead to models ma and mb that are related and should be close. The difference is
that merging the data sets assumes that the signal-to-noise ratios are about
the same in the two experiments. Merging the models allows one model to
much more uncertain, e.g, due to more disturbances in that experiment. If the
conditions are about the same, it is recommended to merge data rather than
models, since this is more efficient and typically involves better conditioned
calculations.
4-113
midprefs
Purpose
Select a directory for idprefs.mat, a file that stores the graphical user
interface s start-up information.
Syntax
midprefs
midprefs(path)
Description
The graphical user interface ident allows a large number of variables for
customized choices. These include the window layout, the default choices of
plot options, and names and directories of the four most recent sessions with
ident. This information is stored in the file idprefs.mat, which should be
placed on the user s MATLABPATH. The default, automatic location for this file is
in the same directory as the user s startup.m file.
4midprefs
4-114
misdata
4misdata
Purpose
Syntax
Datae = misdata(Data)
Datae = misdata(Data,Model)
Datae = misdata(Data,Maxiter,Tol)
Description
Data is input-output data in the iddata object format. Missing data samples
(both in inputs and in outputs) are entered as NaN.
Datae is an iddata object where the missing data has been replaced by
reasonable estimates.
Model is any idmodel (idarx, idgrey, idpoly, idss) used for the
reconstruction of missing data.
iterations will be terminated when the difference between two consecutive data
estimates differ by less than tol%. The default value of tol is 1.
Algorithm
4-115
nkshift
4nkshift
Purpose
Syntax
Datas = nkshift(Data,nk)
Description
is related to
m2 = pem(nkshift(dat,nk),4);
such that m1 and m2 are the same models, but m1 stores the delay information
for use when frequency responses etc. are computed.
Note the difference with the idss and idpoly property nk
m3 = pem(dat,4,'nk',nk)
See Also
4-116
noisecnv
4noisecnv
Purpose
Syntax
mod1 = noisecnv(mod)
mod2 = noisecnv(mod,'norm')
Description
The noise input channels in mod are converted as follows: Consider a model
with both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that mod.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v:
y = Gu + HLv
cov ( v ) = I
mod1 = noisecnv(mod) converts the model to a representation of the system
[G H] with nu+ny inputs and ny outputs. All input are treated as measured,
and mod1 does not have any noise model. The former noise input channels
have names e@yname, where yname is the name of the corresponding output.
mod2 = noisecnv(mod,'norm') converts the model to a representation of the
system [G HL] with nu+ny inputs and ny outputs. All input are treated as
measured, and mod2 does not have any noise model. The former noise input
channels have names v@yname, where yname is the name of the
corresponding output. Note that the noise variance matrix factor L typically
is uncertain (has a non-zero covariance). This is taken into account in the
uncertainty description of mod2.
If mod is a time series, that is nu = 0, mod1 is a model that describes the
transfer function H with measured input channels. Analogously, mod2
describes the transfer function HL.
4-117
noisecnv
4-118
nuderst
Purpose
4nuderst
Syntax
nds = nuderst(pars)
Description
The function pem uses numerical differentiation with respect to the model
parameters when applied to state-space structures. The same is true for many
functions that transform model uncertainties to other representations.
The step-size used in these numerical derivatives is determined by the M-file
nuderst. The output argument nds is a row vector whose k-th entry gives the
increment to be used when differentiating with respect the k-th element of the
parameter vector pars.
The default version of nuderst uses a very simple method. The step-size is the
maximum of 10 4 times the absolute value of the current parameter and 10 7 .
You can adjust this to the actual value of the corresponding parameter by
editing nuderst. Note that the nominal value, for example 0, of a parameter
may not reflect its normal size.
4-119
nyquist
Purpose
4nyquist
Syntax
nyquist(m)
[fr,w] = nyquist(m)
[fr,w,covfr] = nyquist(m)
nyquist(m1,m2,m3,...,w)
nyquist(m1,'PlotStyle1',m2,'PlotStyle2',...)
nyquist(m1,m2,m3,..'sd*5',sd,'mode',mode)
Description
When sd is specified as a number larger than zero, confidence regions will also
be plotted. These are ellipses in the complex plane and correspond to the region
4-120
nyquist
Arguments
no plot is drawn on the screen. If m has ny outputs and nu inputs, and w contains
nw frequencies, then fr is an ny-by-nu-by-Nw array such that fr(ky,ku,k) gives
the complex-valued frequency response from input ku to output ky at the
frequency w(k). For a SISO model, use fr(:) to obtain a vector of the frequency
response. The uncertainty information covfr is a 5-D array where
covfr(ky,ku,k,:,:)) is the 2-by-2 covariance matrix of the response from
input ku to output ky at frequency w(k). The 1,1 element is the variance of the
real part, the 2,2 element the variance of the imaginary part and the 1,2 and
2,1 elements the covariance between the real and imaginary parts.
squeeze(covfr(ky,ku,k,:,:)) gives the covariance matrix of the
corresponding response.
If m is a time series (no input), fr is returned as the (power) spectrum of the
outputs; an ny-by-ny-by-Nw array. Hence fr(:,:,k) is the spectrum matrix at
frequency w(k). The element fr(k1,k2,k) is the cross spectrum between
outputs k1 and k2 at frequency w(k). When k1=k2, this is the real-valued power
spectrum of output k1. covfr is then the covariance of the spectrum fr, so that
covfr(k1,k1,k) is the variance of the power spectrum of output k1 at
frequency w(k). No information about the variance of the cross spectra is
normally given. (That is, covfr(k1,k2,k) = 0 for k1 not equal to k2.)
4-121
nyquist
Examples
See Also
4-122
g = spa(data)
m = n4sid(data,3)
nyquist(g,m,3)
bode, etfe, ffplot, idfrd, spa
n4sid
Purpose
4n4sid
Syntax
m = n4sid(data)
m = n4sid(data,order,'Property1',Value1,...,'PropertyN',ValueN)
Description
The function n4sid estimates models in state-space form, and returns them as
an idss object m. It handles an arbitrary number of input and outputs,
including the time-series case (no input). The state-space model is in the
innovations form
x ( t + Ts ) = Ax ( t ) + Bu ( t ) + Ke ( t )
y ( t ) = Cx ( t ) + Du ( t ) + e ( t )
m: The resulting model as an idss object.
data: An iddata object containing the output-input data.
order: The desired order of the state-space model. If order is entered as a row
vector (like order = [1:10]), preliminary calculations for all the indicated
orders are carried out. A plot will then be given that shows the relative
importance of the dimension of the state vector. More precisely, the singular
values of the Hankel matrices of the impulse response for different orders are
graphed. You will be prompted to select the order, based on this plot. The idea
is to choose an order such that the singular values for higher orders are
comparatively small. If order = 'best', a model of best (default choice) order
is computed, among the orders 1:10. This is the default choice of order.
The list of property name/property value pairs may contain any idss and
algorithm properties. See idss and Algorithm Properties.
idss properties that are of particular interest for n4sid are:
nk: The delays from the inputs to the outputs, a row vector with the same
number of entries as the number of input channels. Default is nk = [1 1 ...
1]. Note that delays being 0 or 1 show up as zeros or estimated parameters in
the D matrix. Delays larger than 1 means that a special structure of the A, B
and C matrices are used to accommodate the delays. This also means that the
actual order of the state-space model will be larger than order.
CovarianceMatrix (can be abbreviated to 'co'): Setting CovarianceMatrix
to 'None' will block all calculations of uncertainty measures. These may
4-123
n4sid
take the major part of the computation time. Note that, for a 'Free'
parameterization, the individual matrix elements cannot be associated with
any variance (these parameters are not identifiable). Instead, the resulting
model m stores a hidden state-space model in canonical form, that contains
covariance information. This is used when the uncertainty of various
input-output properties are calculated. It can also be retrieved by m.ss =
'can'. The actual covariance properties of n4sid estimates are not known
today. Instead the Cramer-Rao bound is computed and stored as an
indication of the uncertainty.
DisturbanceModel: Setting DisturbanceModel to None will deliver a model
with K = 0. This will have no direct effect on the dynamics model, other that
that the default choice of N4Horizon will not involve past outputs.
InitialState: The initial state is always estimated for better accuracy.
However. it is returned with m only if InitialState = Estimate .
Algorithm properties that are special interest are:
Focus: Assumes the values 'Prediction' (default), 'Simulation',
Stability , or any SISO linear filter (given as an LTI or idmodel object, or
as filter coefficients. See Algorithm Properties.) Setting 'Focus' to
'Simulation' chooses weights that should give a better simulation
performance for the model. In particular, a stable model is guaranteed.
Selecting a linear filter will focus the fit to the frequency ranges that are
emphasized by this filter.
N4Weight: This property determines some weighting matrices used in the
singular-value decomposition that is a central step in the algorithm. Two
choices are offered: 'MOESP' that corresponds to the MOESP algorithm by
Verhaegen, and 'CVA' which is the canonical variable algorithm by
Larimore. The default value is 'N4Weight' = 'Auto', which gives an
automatic choice between the two options. m.EstimationInfo.N4Weight
tells you what the actual choice turned out to be.
N4Horizon: Determines the prediction horizons forward and backward, used
by the algorithm. This is a row vector with three elements: N4Horizon =[r
sy su], where r is the maximum forward prediction horizon, i.e., the
algorithms uses up to r-step ahead predictors. sy is the number of past
outputs, and su is the number of past inputs that are used for the
predictions. See pages 209-210 in Ljung(1999) for the exact meaning of this.
These numbers may have a substantial influence of the quality of the
4-124
n4sid
resulting model, and there are no simple rules for choosing them. Making
'N4Horizon' a k-by-3 matrix, means that each row of 'N4Horizon' will be
tried out, and the value that gives the best (prediction) fit to data will be
selected. (This option cannot be combined with selection of model order.) If
the property 'Trace' is 'On', information about the results will be given in
the MATLAB command window.
4-125
n4sid
Algorithm
Examples
Build a fifth order model from data with three inputs and two outputs. Try
several choices of auxiliary orders. Look at the frequency response of the model.
z = iddata([y1 y2],[ u1 u2 u3]);
m = n4sid(z,5,'n4h',[7:15]','trace','on');
bode(m,'sd',3)
See Also
References
4-126
oe
Purpose
4oe
Syntax
m = oe(data,orders)
m = oe(data,'nb',nb,'nf',nf,'nk',nk)
m = oe(data,orders,'Property1',Value1,'Property2',Value2,...)
Description
B( q)
y ( t ) = ------------ u ( t nk ) + e ( t )
F(q)
are estimated using a prediction error method.
data is an iddata object containing the output-input data. The structure
information can either be given explicitly as
(...,'nb',nb,'nf',nf,'nk',nk,...), or in the argument orders given as
orders = [nb nf nk]
The parameters nb and nf are the orders of the Output-Error model and nk is
the delay. Specifically,
nb:
B ( q ) = b 1 + b2 q
nf:
F ( q ) = 1 + f1 q
+ + b nb q
+ + f nf q
nb + 1
nf
4-127
oe
The structure and the estimation algorithm are affected by any property
name/property value pairs that are set in the input argument list. Useful
properties are 'Focus', 'InitialState', 'InputDelay', 'SearchDirection',
'MaxIter', 'Tolerance', 'LimitError', 'FixedParameter', and 'Trace'.
See Algorithm Properties, idpoly, and idmodel for details of these properties
and their possible values.
oe does not support multi-output models. Use state-space model for this case
(see n4sid and pem).
Algorithm
See Also
4-128
pe
Purpose
4pe
Compute the prediction errors associated with a model and a data set.
Syntax
e = pe(m,data)
[e,x0] = pe(m,data,init)
Description
data is the output-input data set, given as an iddata object, and m is any
idmodel object.
e is returned as an iddata object, so that e.OutputData contains the prediction
errors that result when model m is applied to the data:
1
e ( t ) = H ( q ) [ y ( t ) G ( q )u ( t ) ]
The argument init determines how to deal with the initial conditions:
init ='estimate' means that the initial state is chosen so that the norm of
prediction error is minimized. This initial state is returned as x0.
init = 'zero' sets the initial state to zero.
init = 'model' used the model s internally stored initial state.
init = x0, where x0 is a column vector of appropriate dimension uses that
value as initial state.
If init is not specified, the model property m.InitialState is used, so that
'Estimate', 'Backcast' and 'Auto' set init = 'Estimate', while
m.InitialState = 'Zero' sets init = 'zero', and 'Fixed' and 'Model'
set init = 'model'.
The output argument x0 is the used value of the initial state. If data contains
several experiments, x0 will be a matrix, containing the initial states from each
experiment.
See Also
idmodel, resid
4-129
pem
Purpose
4pem
Syntax
m
m
m
m
m
m
m
Description
pem is the basic estimation command in the toolbox and covers a variety of
situations.
=
=
=
=
=
=
=
pem(data)
pem(data,mi)
pem(data,mi,'Property1',Value1,...,'PropertyN',ValueN)
pem(data,orders)
pem(data,'nx',ssorder)
pem(data,'na',na,'nb',nb,'nc',nc,'nd',nd,'nf',nf,'nk',nk)
pem(data,orders,'Property1',Value1,...,'PropertyN',ValueN)
4-130
pem
Properties
4-131
pem
Trace, with possible values 'Off', 'On', 'Full', that governs the
information sent to the MATLAB command window.
For black-box state-space models, also 'N4Weight' and 'N4Horizon' will affect
the result, since these models are initialized with n4sid estimate. See the
reference page for n4sid.
Typical idmodel properties to affect are (see idmodel properties for more
details):
Ts, the sampling interval. Set 'Ts'=0 to obtain a continuous-time state-space
model. For discrete-time models, 'Ts' is automatically set to sampling
interval of the data. Note that, in the black box case, it is usually better to
first estimate a discrete-time model, and then convert that to continuous
time by d2c.
nk, the time delays from the inputs (not applicable to structured state-space
models). Time delays specified by 'nk', will be included in the model.
DisturbanceModes determines the parameterization of K for free and
canonical state-space parameterizations, as well as for idgrey models.
InitialState. The initial state may have a substantial influence on the
estimation result for system with slow responses. It is most pronounced for
Output-Error models (K=0 for state-space, and na=nc=nd=0 for input/output
models). The default value 'Auto', estimates the influence of the initial state
and sets the value to 'Estimate', 'Backcast' or 'Zero', based on this
effect. Possible values of 'InitialState' are 'Auto', 'Estimate',
'Backcast', 'Zero' and 'Fixed'. See Initial State in the Tutorial.
Examples
Here is an example of a system with three inputs and two outputs. A canonical
form state-space model of order 5 is sought.
z = iddata([y1 y2],[ u1 u2 u3]);
m = pem(z,5,'ss','can')
4-132
pem
Comparing the models (compare automatically matches the channels using the
channel names).
compare(z,m,ma)
Algorithm
pem uses essentially the same algorithm as armax with modifications to the
computation of prediction errors and gradients.
See Also
4-133
plot (iddata)
Purpose
4plot (iddata)
Syntax
plot(data)
plot(d1,...,dN)
plot(d1,PlotStyle1,...,dN,PlotStyleN)
Description
One plot for each I/O channel combination is produced. Pressing the Return
key advances the plot.
To plot a specific interval, use plot(data(200:300)). To plot specific
input/output channels, use plot(data(:,ky,ku)), consistent with the
subreferencing of iddata objects. (See iddata).
If data.intersample = 'zoh', the input is piecewise constant between
sampling points, and it is then graphed accordingly.
To plot several iddata sets d1,...,dN, use plot(d1,...,dN). I/O channels
with the same experiment name, input name, and output name, will always be
plotted in the same plot.
With PlotStyle, the color, linestyle, and marker of each data set can be
specified
plot(d1,'y:*',d2,'b')
See Also
4-134
iddata
plot (idmodepolydata)
Purpose
4plot (idmodepolydata)
Syntax
[A,B,C,D,F] = polydata(m)
[A,B,C,D,F,dA,dB,dC,dD,dF] = polydata(m)
Description
See Also
4-135
predict
Purpose
4predict
Syntax
yp = predict(m,data)
[yp,mpred] = predict(m,data,k,init)
Description
data is the output-input data as an iddata object, and m is any idmodel object
(idpoly, idss, idgrey, or idarx)
The argument k indicates that the k step ahead prediction of y according to the
model m is computed. In the calculation of yp(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 k is 1.
The output yp is an iddata object containing the predicted values as
OutputData. The output argument mpred contains the k-step ahead predictor.
This is given as a cell array, whose k:th entry is an idpoly model for the
predictor of output number k.
init determines how to deal with the initial state:
init ='estimate': The initial state is set to value that minimizes the norm
of the prediction error associated with the model and the data.
init = 'zero' sets the initial state to zero.
init = 'model' used the model's internally stored initial state
init = x0, where x0 is a column vector of appropriate dimension uses that
value as initial state
If init is not specified, the model property m.InitialState is used, so that
'Estimate', 'Backcast' and 'Auto' set init = 'Estimate', while
m.InitialState = 'Zero' sets init = 'zero', and 'Model' and 'Fixed' set
init = 'model'.
An important use of predict is to evaluate a models properties in the
mid-frequency range. Simulation with sim (which conceptually corresponds to
k = inf) can lead to levels that drift apart, since the low frequency behavior is
emphasized. One step ahead prediction is not a powerful test of the model s
properties, since the high frequency behavior is stressed. The trivial predictor
4-136
predict
y ( t ) = y ( t 1 ) can give good predictions in case the sampling of the data is
fast.
Another important use of predict is to evaluate models of time series. The
natural way of studying a time-series model s ability to reproduce
observations is to compare its k-step ahead predictions with actual data.
Note that for Output-Error models, there is no difference between the k-step
ahead predictions and the simulated output, since, by definition, Output-Error
models only use past inputs to predict future outputs.
Algorithm
The model is evaluated in state-space form, and the state equations are
simulated k-steps ahead with initial value x ( t k ) = x ( t k ) , where x ( t k )
is the Kalman filter state estimate.
Examples
Simulate a time series, estimate a model based on the first half of the data, and
evaluate the four step ahead predictions on the second half.
m0 = idpoly([1 -0.99],[],[1 -1 0.2]);
e = iddata([],randn(400,1));
y = sim(m0,e);
m = armax(y(1:200),[1 2]);
yp = predict(m,y,4);
plot(y(201:400),yp(201:400))
See Also
compare, sim, pe
4-137
present
Purpose
4present
Syntax
present(m)
Description
This function displays the model m, together with the estimated standard
deviations of the parameters, loss function, and Akaike s Final Prediction
Error (FPE) Criterion (which essentially equals the AIC). It also displays
information about how m was created.
present thus gives more detailed information about the model than the
4-138
pzmap
Purpose
4pzmap
Syntax
pzmap(m)
pzmap(m,'sd',sd)
pzmap(m1,m2,m3,...)
pzmap(m1,'PlotStyle1',m2,'PlotStyle2',...,'sd',sd)
pzmap(m1,m2,m3,..,'sd',sd,mode,axis)
Description
The zeros and poles of m are graphed, with o denoting zeros and x denoting
poles. Poles and zeros at infinity are ignored. For discrete-time models, zeros
and poles at the origin are also ignored.
If sd has a value larger than zero, confidence regions around the poles and
zeros are also graphed. The regions corresponding to sd standard deviations
are marked. The default value is sd = 0. Note that the confidence regions may
sometimes stretch outside the plot, but they are always symmetric around the
indicated zero or pole.
If the poles and zeros are associated with a discrete-time model, a unit circle is
also drawn. For continuous-time models the real and imaginary axes are
drawn
When mi contain information about several different input/output channels,
there are some options:
mode = 'sub' splits the screen into several plots, one for each input/output
channel. These are based on the InputName and OutputName properties
associated with the different models.
mode = 'same' gives all plots in the same diagram. Pressing the Return key
advances the plots.
mode = 'sep' erases the previous plot before the next channel pair is treated.
as
axis = [-s s -s s]
4-139
pzmap
The colors associated with the different models can be selected by the
arguments PlotStyle. Use PlotStyle = 'b', 'g', etc. Markers and line styles
are not used.
The noise input channels in m are treated as follows: Consider a model m with
both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that m.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v
y = Gu + HLv
cov ( v ) = I
Then:
pzmap(m) plots the zeros and poles of the transfer function G.
pzmap(m('n')) plots the zeros and poles of the transfer function H. (ny
inputs and ny outputs). The input channels have names e@yname, where
yname is the name of the corresponding output.
If m is a time series, that is nu = 0, pzmap(m) plots the zeros and poles of the
the transfer function H.
pzmap(noisecnv(m)) plots the zeros and poles of the transfer function [G H]
(nu+ny inputs and ny outputs). The noise input channels have names
e@yname, where yname is the name of the corresponding output.
pzmap(noisecnv(m,'norm')) plots the zeros and poles of the transfer
function [G HL] (nu+ny inputs and ny outputs. The noise input channels
have names v@yname, where yname is the name of the corresponding output.
Examples
shows all zeros and poles of two models along with the confidence regions
corresponding to three standard deviations.
See Also
4-140
idmodel, zpkdata
rarmax
Purpose
4rarmax
Syntax
thm = rarmax(z,nn,adm,adg)
[thm,yhat,P,phi,psi] = rarmax(z,nn,adm,adg,th0,P0,phi0,psi0)
Description
where na, nb, and nc are the orders of the ARMAX model, and nk is the delay.
Specifically,
1
na:
A ( q ) = 1 + a1 q
nb:
B ( q ) = b1 + b 2 q
nc:
C ( q ) = 1 + c1 q
+ + a na q
na
+ + b nb q
+ + c nc q
nb + 1
nc
4-141
rarmax
thm(k,:) = [a1,a2,...,ana,b1,...,bnb,c1,...,cnc]
yhat is the predicted value of the output, according to the current model, i.e.,
row k of yhat contains the predicted value of y(k) based on all past data.
The actual algorithm is selected with the two arguments adm and adg. These
are described under rarx.
The input argument th0 contains the initial value of the parameters, a row
vector, consistent with the rows of thm. The default value of th0 is all zeros.
The arguments P0 and P are the initial and final values, respectively, of the
scaled covariance matrix of the parameters. See rarx. The default value of P0
is 104 times the unit matrix. The arguments phi0, psi0, phi, and psi contain
initial and final values of the data vector and the gradient vector, respectively.
The sizes of these depend in a rather complicated way on the chosen model
orders. The normal choice of phi0 and psi0 is to use the outputs from a
previous call to rarmax with the same model orders. (This call could of course
be a dummy call with default input arguments.) The default values of phi0 and
psi0 are all zeros.
Note that the function requires that the delay nk be larger than 0. If you want
nk=0, shift the input sequence appropriately and use nk=1.
Algorithm
Examples
Compute and plot, as functions of time, the four parameters in a second order
ARMA model of a time series given in the vector y. The forgetting factor
algorithm with a forgetting factor of 0.98 is applied.
thm = rarmax(y,[2 2],'ff',0.98);
plot(thm)
4-142
rarx
Purpose
4rarx
Syntax
thm = rarx(z,nn,adm,adg)
[thm,yhat,P,phi] = rarx(z,nn,adm,adg,th0,P0,phi0)
Description
where na and nb are the orders of the ARX model, and nk is the delay.
Specifically,
1
na:
A ( q ) = 1 + a1 q
nb:
B ( q ) = b 1 + b2 q
+ + a na q
na
+ + b nb q
nb + 1
unu]
and by allowing nb and nk to be row vectors defining the orders and delays
associated with each input.
Only single-output models are handled by rarx.
4-143
rarx
The estimated parameters are returned in the matrix thm. The k-th row of thm
contains the parameters associated with time k, i.e., they are based on the data
in the rows up to and including row k in z. Each row of thm contains the
estimated parameters in the following order.
thm(k,:) = [a1,a2,...,ana,b1,...,bnb]
In the case of a multi-input model, all the b parameters associated with input
number 1 are given first, and then all the b parameters associated with input
number 2, and so on.
yhat is the predicted value of the output, according to the current model, i.e.,
row k of yhat contains the predicted value of y(k) based on all past data.
The actual algorithm is selected with the two arguments adg and adm. These
are described in Recursive Parameter Estimation in the Tutorial chapter.
The options are as follows.
With adm ='ff' and adg = lam the forgetting factor algorithm
(Equation 3-60abd)+(Equation 3-62) is obtained with forgetting factor = lam.
This is what is often referred to as Recursive Least Squares, RLS. In this case
the matrix P (see below) has the following interpretation: R 2 /2 P is
approximately equal to the covariance matrix of the estimated parameters.
Here R 2 is the variance of the innovations (the true prediction errors e(t) in
(Equation 3-57)
With adm ='ug' and adg = gam, the unnormalized gradient algorithm
(Equation 3-60abc) + (Equation 3-63) is obtained with gain gamma= gam. This
algorithm is commonly known as normalized Least Mean Squares, LMS.
Similarly, adm ='ng' and adg = gam give the normalized gradient or
Normalized Least Mean Squares, NLMS algorithm (Equation 3-60abc) +
(Equation 3-64). In these cases, P is not defined or applicable.
With adm ='kf' and adg = R1 the Kalman Filter Based algorithm
(Equation 3-60) is obtained with R2= 1 and R1 = R1. If the variance of the
innovations e(t) is not unity but R 2 ; then R 2 P is the covariance matrix of the
parameter estimates, while R 1 =R1 / R 2 is the covariance matrix of the
parameter changes in (Equation 3-58).
4-144
rarx
The input argument th0 contains the initial value of the parameters; a row
vector, consistent with the rows of thm. (See above.) The default value of th0 is
all zeros.
The arguments P0 and P are the initial and final values, respectively, of the
scaled covariance matrix of the parameters. The default value of P0 is 104
times the identity matrix.
The arguments phi0 and phi contain initial and final values, respectively, of
the data vector.
( t ) = [ y ( t 1 ), , y ( t na ), u ( t 1 ), u ( t nb nk + 1 ) ]
Then, if
z = [y(1),u(1); ...
;y(N),u(N)]
you have phi0 = ( 1 ) and phi = ( N ) . The default value of phi0 is all zeros.
For online use of rarx, use phi0, th0, and P0 as the previous outputs phi, thm
(last row), and P.
Note that the function requires that the delay nk be larger than 0. If you want
nk=0, shift the input sequence appropriately and use nk=1. See nkshift.
Examples
Adaptive noise canceling: The signal y contains a component that has its origin
in a known signal r. Remove this component by estimating, recursively, the
system that relates r to y using a sixth order FIR model together with the
NLMS algorithm.
z = [y r];
[thm,noise] = rarx(z,[0 6 1],'ng',0.1);
% noise is the adaptive estimate of the noise
% component of y
plot(y-noise)
If the above application is a true online one, so that you want to plot the best
estimate of the signal y - noise at the same time as the data y and u become
available, proceed as follows.
phi = zeros(6,1); P=1000eye(6);
th = zeros(1,6); axis([0 100 -2 2]);
plot(0,0,''), hold on
% The loop:
4-145
rarx
while ~abort
[y,r,abort] = readAD(time);
[th,ns,P,phi] = rarx([y r],'ff',0.98,th,P,phi);
plot(time,y-ns,'')
time = time +Dt
end
This example uses a forgetting factor algorithm with a forgetting factor of 0.98.
readAD represents an M-file that reads the value of an A/D converter at the
indicated time instant.
4-146
rbj
Purpose
4rbj
Syntax
thm = rbj(z,nn,adm,adg)
[thm,yhat,P,phi,psi] = ... rbj(z,nn,adm,adg,th0,P0,phi0,psi0)
Description
where nb, nc, nd, and nf are the orders of the Box-Jenkins model, and nk is the
delay. Specifically,
nb:
B ( q ) = b1 + b 2 q
nc:
C ( q ) = 1 + c1 q
nd:
D ( q ) = 1 + d1 q
nf:
F ( q ) = 1 + f1 q
+ + c nc q
+ + b nb q
nc
+ + d nd q
+ + f nf q
nb + 1
nd
nf
4-147
rbj
yhat is the predicted value of the output, according to the current model, i.e.,
row k of yhat contains the predicted value of y(k) based on all past data.
The actual algorithm is selected with the two arguments adm and adg. These
are described under rarx.
The input argument th0 contains the initial value of the parameters, a row
vector, consistent with the rows of thm. (See above.) The default value of th0 is
all zeros.
The arguments P0 and P are the initial and final values, respectively of the
scaled covariance matrix of the parameters. See rarx. The default value of P0
is 104 times the unit matrix. The arguments phi0, psi0, phi, and psi contain
initial and final values of the data vector and the gradient vector, respectively.
The sizes of these depend in a rather complicated way on the chosen model
orders. The normal choice of phi0 and psi0 is to use the outputs from a
previous call to rbj with the same model orders. (This call could, of course, be
a dummy call with default input arguments.) The default values of phi0 and
psi0 are all zeros.
Note that the function requires that the delay nk is larger than 0. If you want
nk=0, shift the input sequence appropriately and use nk=1.
Algorithm
4-148
resample
Purpose
4resample
Syntax
datar = resample(data,P,Q)
datar = resample(data,P,Q,,filter_order)
Description
Algorithm
Examples
4-149
resid
Purpose
4resid
Syntax
resid(m,data)
resid(m,data,Type)
resid(m,data,Type,M)
e = resid(m,data);
Description
In all cases the residuals e associated with the data and the model are
computed. This is done as in the command pe with a default choice of init.
When called without output arguments, resid produces a plot. The plot can be
of three kinds depending on the argument Type:
Type = 'Corr' (default): The autocorrelation function of e and the cross
correlation between e and the input(s) u are computed and displayed. The
99% confidence intervals for these values are also computed and shown as a
yellow region. The computation of the confidence region is done assuming e
to be white and independent of u. The functions are displayed up to lag M,
which is 25 by default.
Type = 'ir': The impulse response (up to lag M, which is 25 by default) from
the input to the residuals is plotted with a 99% confidence region around zero
marked as a yellow area. Negative lags up tp M/4 are also included to
investigate feedback effects. (The result is the same as
impulse(e,'sd',2.58,'fill',M).)
Type = 'fr': The frequency response from the input to the residuals (based
on a high order FIR model) is shown as a Bode plot. A 99% confidence region
around zero is also marked as a yellow area.
With an output argument, no plot is produced, and e is returned with the
residuals (prediction errors) associated with the model and the data. It is an
iddata object with the residuals as outputs and the input in data as inputs.
That means that e can be directly used to build model error models, i.e., models
that describe the dynamics from the input to the residuals (which should be
negligible if m is a good description of the system).
4-150
resid
See Model Structure Selection and Validation in the Tutorial chapter for
more information.
Examples
To compute a model error model, that is, a model to input to the residuals to
see if any essential unmodeled dynamics are left,
e = resid(m,data);
me = arx(e,[10 10 0]);
bode(me,'sd',3,fill')
See Also
References
4-151
roe
Purpose
4roe
Syntax
thm = roe(z,nn,adm,adg)
[thm,yhat,P,phi,psi] = roe(z,nn,adm,adg,th0,P0,phi0,psi0)
Description
where nb and nf are the orders of the Output-Error model, and nk is the delay.
Specifically,
nb:
B ( q ) = b1 + b2 q
nf:
F ( q ) = 1 + f1 q
+ + b nb q
+ + f nf q
nb + 1
nf
4-152
roe
The actual algorithm is selected with the two arguments adg and adm. These
are described under rarx.
The input argument th0 contains the initial value of the parameters, a row
vector, consistent with the rows of thm. (See above.) The default value of th0 is
all zeros.
The arguments P0 and P are the initial and final values, respectively, of the
scaled covariance matrix of the parameters. See rarx. The default value of P0
is 104 times the unit matrix. The arguments phi0, psi0, phi, and psi contain
initial and final values of the data vector and the gradient vector, respectively.
The sizes of these depend in a rather complicated way on the chosen model
orders. The normal choice of phi0 and psi0 is to use the outputs from a
previous call to roe with the same model orders. (This call could be a dummy
call with default input arguments.) The default values of phi0 and psi0 are all
zeros.
Note that the function requires that the delay nk is larger than 0. If you want
nk=0, shift the input sequence appropriately and use nk=1.
Algorithm
See Also
4-153
rpem
Purpose
4rpem
Syntax
thm = rpem(z,nn,adm,adg)
[thm,yhat,P,phi,psi] = rpem(z,nn,adm,adg,th0,P0,phi0,psi0)
Description
where na, nb, nc, nd, and nf are the orders of the model, and nk is the delay.
For multi-input systems nb, nf, and nk are row vectors giving the orders and
delays of each input. See Polynomial Representation of Transfer Functions in
the Tutorial for an exact definition of the orders.
The estimated parameters are returned in the matrix thm. The k-th row of thm
contains the parameters associated with time k, i.e., they are based on the data
in the rows up to and including row k in z. Each row of thm contains the
estimated parameters in the following order.
thm(k,:) = [a1,a2,...,ana,b1,...,bnb,...
c1,...,cnc,d1,...,dnd,f1,...,fnf]
For multi-input systems the B part in the above expression is repeated for each
input before the C part begins, and also the F part is repeated for each input.
This is the same ordering as in m.par.
yhat is the predicted value of the output, according to the current model, i.e.,
row k of yhat contains the predicted value of y(k) based on all past data.
The actual algorithm is selected with the two arguments adg and adm. These
are described under rarx.
4-154
rpem
The input argument th0 contains the initial value of the parameters;, a row
vector, consistent with the rows of thm. (See above.) The default value of th0 is
all zeros.
The arguments P0 and P are the initial and final values, respectively, of the
scaled covariance matrix of the parameters. See rarx. The default value of P0
is 104 times the unit matrix. The arguments phi0, psi0, phi, and psi contain
initial and final values of the data vector and the gradient vector, respectively.
The sizes of these depend in a rather complicated way on the chosen model
orders. The normal choice of phi0 and psi0 is to use the outputs from a
previous call to rpem with the same model orders. (This call could be a dummy
call with default input arguments.) The default values of phi0 and psi0 are all
zeros.
Note that the function requires that the delay nk is larger than 0. If you want
nk=0 , shift the input sequence appropriately and use nk=1.
Algorithm
See Also
4-155
rplr
Purpose
4rplr
Syntax
thm = rplr(z,nn,adm,adg)
[thm,yhat,P,phi] = rplr(z,nn,adm,adg,th0,P0,phi0)
Description
This is a direct alternative to rpem and has essentially the same syntax. See
rpem for an explanation of the input and output arguments.
rplr differs from rpem in that it uses another gradient approximation. See
Section 11.5 in Ljung (1999). Several of the special cases are well known
algorithms.
When applied to ARMAX models, (nn = [na nb nc 0 0 nk]), rplr gives the
Extended Least Squares (ELS) method. When applied to the output error
structure (nn = [0 nb 0 0 nf nk]) the method is known as HARF in the
adm = 'ff' case and SHARF in the adm = 'ng' case.
rplr can also be applied to the time-series case in which an ARMA model is
estimated with
z = y
nn = [na nc]
You can thus use rplr as an alternative to rarmax for this case.
See Also
4-156
segment
Purpose
4segment
Syntax
segm = segment(z,nn)
[segm,V,thm,R2e] = segment(z,nn,R2,q,R1,M,th0,P0,ll,mu)
Description
A ( q )y ( t ) = B ( q )u ( t nk ) + C ( q )e ( t )
assuming that the model parameters are piece-wise constant over time. It
results in a model that has split the data record into segments over which the
model remains constant. The function models signals and systems that may
undergo abrupt changes.
The input-output data are contained in z, which is either an iddata object or a
matrix z = [y u] where y and u are column vectors. If the system has several
inputs, u has the corresponding number of columns.
The argument nn defines the model order. For the ARMAX model
nn = [na nb nc nk]
where na, nb, and nc are the orders of the corresponding polynomials. See
Polynomial Representation of Transfer Functions in the Tutorial. Moreover
nk is the delay. If the model has several inputs, nb and nk are row vectors,
giving the orders and delays for each input.
For an ARX model (nc = 0) enter
nn = [na nb nk]
The output argument segm is a matrix, whose k-row contains the parameters
corresponding to time k. This is analogous to the output argument thm in rarx
and rarmax. The output argument thm of segment contains the corresponding
model parameters that have not yet been segmented. That is, they are not
4-157
segment
value is 5.
th0 is the initial value of the parameters. Its default is zero. P0 is the initial
The most critical parameter for you to choose is R2. It is usually more robust to
have a reasonable guess of R2 than to estimate it. Typically, you need to try
different values of R2 and evaluate the results. (See the example below.)
sqrt(R2) corresponds to a change in the value y(t) that is normal, giving no
indication that the system or the input might have changed.
4-158
segment
Algorithm
Examples
Check how the algorithm segments a sinusoid into segments of constant levels.
Then use a very simple model y(t) = b1 * 1, where 1 is a fake input and b 1
describes the piecewise constant level of the signal y(t) (which is simulated as
a sinusoid).
y = sin([1:50]/3)';
thm = segment([y,ones(y)],[0 1 1],0.1);
plot([thm,y])
By trying various values of R2 (0.1 in the above example), more levels are
created as R2 decreases.
4-159
selstruc
Purpose
4selstruc
Syntax
nn = selstruc(v)
[nn,vmod] = selstruc(v,c)
Description
The default value of c is 'plot'. The plot shows the percentage of the output
variance that is not explained by the model, as a function of the number of
parameters used. Each value shows the best fit for that number of parameters.
By clicking in the plot you can examine which orders are of interest. By clicking
on Select the variable nn is returned in the workspace as the optimal model
structure for your choice of number of parameters. Several choices can be
made.
c = 'aic' gives no plots, but returns in nn the structure that minimizes
Akaike s Information Criterion (AIC),
2d
V mod = V 1 + -------
N
where V is the loss function, d is the total number of parameters in the
structure in question, and N is the number of data points used for the
estimation.
c = 'mdl' returns in nn the structure that minimizes Rissanen s Minimum
Description Length (MDL) criterion.
log ( N )
V mod = V 1 + d
----------------------
N
When c equals a numerical value, the structure that minimizes
cd
V mod = V 1 + ------
N
is selected.
The output argument vmod has the same format as v, but it contains the
logarithms of the accordingly modified criteria.
4-160
selstruc
Example
V = arxstruc(data(1:200),data(201:400),struc(1:10,1:10,1:10))
nn = selstruc(V,0); %best fit to validation data
m = arx(data,nn)
4-161
set
Purpose
4set
Syntax
set(m,'Property',Value)
set(m,'Property1',Value1,...'PropertyN',ValueN)
set(m,'Property')
set(m)
Description
set is used to set or modify the properties of any of the objects in the toolbox
(iddata, idmodel, idgrey, idarx, idpoly, idss). See the corresponding
reference entries for a complete list of properties.
set(m,'Property',Value) assigns the value Value to the property of the
object m, specified by the string 'Property'. This string can be the full property
name (e.g., 'SSParameterization') or any unambiguous case-insensitive
abbreviation (e.g., 'ss').
set(m,'Property1',Value1,...'PropertyN',ValueN) sets multiple properties
with a single statement. In certain cases this may be necessary, since the model
m must, e.g., have state-space matrices of consistent dimensions after each set
statement.
set(m,'Property') displays admissible values for the property specified by
'Property'.
set(m) displays all assignable values of m and their admissible values.
4-162
setpname
Purpose
4setpname
Syntax
model = setpname(model)
Description
model is an idmodel object of idarx, idpoly or idss type. The returned model
has the PName property set to a cell array of strings that correspond to the
symbols used in this manual to describe the parameters.
For single input idpoly models, the parameters are called 'a1',
'a2', ...,'fn', just as defined in Section Polynomial Representation of
Transfer Functions.
For multi-input idpoly models, the b- and f-parameters have the output/input
channel number in parenthesis as in 'b1(1,2)', 'f3(1,2)' etc.
For idarx models, the parameter names are as in 'Ar(ky,ku)' for the ky-ku
entry of the matrix in (Equation 3-46) and similarly for the B-parameters.
For idss models the parameters are named for the matrix entries they
represent, like 'A(4,5)', 'K(2,3)' etc.
This function is particularly useful when certain parameters are to be fixed.
See the property 'FixedParameter' under Algorithm Properties.
4-163
sim
Purpose
4sim
Syntax
y = sim(m,ue)
[y, ysd] = sim(m,ue,init)
Description
the sum of the number of inputs and noise sources (= number of outputs). In
the latter case the last inputs in ue are regarded as noise sources and a
noise-corrupted simulation is obtained. The noise is scaled according to the
property m.NoiseVariance in m, so in order to obtain the right noise level
according to the model, the noise inputs should be white noise with zero mean
and unit covariance matrix. If no noise sources are contained in ue, a noise-free
simulation is obtained.
sim returns y containing the simulated output, as an iddata object.
init gives access to the initial states:
init = 'm' (default) uses the model m's internally stored initial state.
init = 'z' uses zero initial state.
init = x0, where x0 is a column vector of appropriate length uses this value
as the initial state.
The second output argument ysd is the standard deviation of the simulated
output.
If m is a continuous-time model, it is first converted to discrete time with the
sampling interval given by ue taking into account the intersample behavior of
the input (ue.InterSample). See the sectionDiscrete and Continuous Time
Models in the Tutorial.
Examples
4-164
=
=
=
=
iddata([],randn(500,1));
iddata([],idinput(500,'prbs'));
sim(m0,[u e]);
[y u]; % An iddata object with y as output and u as input.
sim
See Also
4-165
simsd
Purpose
4simsd
Syntax
simsd(m,u)
simsd(m,u,N,noise,Ky)
Description
Examples
Plot the step response of the model m and evaluate how it varies in view of the
model's uncertainty.
step1 = [zeros(5,1); ones(20,1)];
simsd(m,step1)
See Also
4-166
sim
size
Purpose
4size
Syntax
d = size(m)
[ny,nu,Npar,Nx] = size(model)
[N, ny, nu, Nexp] = size(data)
ny = size(data,2)
ny = size(data,'ny')
size(model)
Description
To access just one of these sizes use size(data,k) for the k-th dimension or
size(data,'N'), size(data,'ny'), etc.
When called with one output argument d = size(data) returns:
d = [N ny nu] if the number of experiments is 1
d = [sum(N) ny nu Ne] if the number of experiments is Ne>1
4-167
size
[ny nu Npar]
For idfrd models, the sizes are:
ny = number of output channels
nu = number of input channels
Nf = number of frequencies
Ns = number of spectrum channels
When size is called with no output arguments in any of these cases, the
information is displayed in the MATLAB command window.
4-168
spa
Purpose
4spa
Syntax
g = spa(data)
g = spa(data,M,w,maxsize)
[g,phi,spe] = spa(data)
Description
spa estimates the transfer function g and the noise spectrum v of the general
linear model
y ( t ) = G ( q )u ( t ) + v ( t )
where v ( ) is the spectrum of v ( t ) .
data contains the output-input data as an iddata object. The data may be
complex-valued.
g is returned as an idfrd object (see idfrd) with the estimate of G ( ei ) at the
frequencies specified by row vector w. The default value of w is
w = [1:128]/128pi/Ts
Changing the value of M exchanges bias for variance in the spectral estimate.
As M is increased, the estimated functions show more detail, but are more
corrupted by noise. The sharper peaks a true frequency function has, the
higher M it needs. See etfe as an alternative for narrowband signals and
systems. See also Estimating Spectra and Frequency Functions in the
Tutorial.
maxsize controls the memory-speed trade-off (see Algorithm Properties).
For time series, where data contains no input channels, g is returned with the
estimated output spectrum and its estimated standard deviation.
4-169
spa
S =
m = M
Examples
4-170
spa
Algorithm
The covariance function estimates are computed using covf. These are
multiplied by a Hamming window of lag size M and then Fourier transformed.
The relevant ratios and differences are then formed. For the default
frequencies, this is done using FFT, which is more efficient than for
user-defined frequencies. For multi-variable systems, a straightforward
for-loop is used.
Note that M = is in Table 6.1 of Ljung (1999). The standard deviations are
computed as on pages 184 and 188 in Ljung (1999).
See Also
4-171
Purpose
Transform the model objects of the System Identification Toolbox to the LTI
models of the Control System Toolbox.
Syntax
sys
sys
sys
sys
Description
mod is any idmodel object: idgrey, idarx, idpoly, idss, idmodel, or idfrd.
However, idfrd objects can only be converted to frd objects.
=
=
=
=
ss(mod)
tf(mod)
zpk(mod)
frd(mod)
sys is returned as the indicated LTI model. The noise input channels in mod are
treated as follows.
Consider a model mod with both measured input channels u (nu channels) and
noise channels e (ny channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that mod.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v.
y = Gu + HLv
cov ( v ) = I
For ss, tf, and zpk, both measured input channels u and normalized noise
input channels v in mod will be input channels in sys. The noise input channels
will belong to the InputGroup 'Noise', while the others belong to the
InputGroup 'Measured'. The names of the noise input channels will be
v@yname, where yname is the name of the corresponding output channel. This
means that the LTI object realizes the transfer function [G HL].
For frd, no noise input information is included.
To transform only the measured input channels in sys, use
sys = ss(mod('m'))
and analogously for tf and zpk. This will give a representation of just G.
4-172
For a time series, (no measured input channels, nu = 0), the LTI
representations in ss, tf, zpk, and frd contain the transfer functions from the
normalized noise sources v to the outputs., that is HL. If the model m has both
measured and noise inputs, sys = ss(mod('n')) will give a representation of
the additive noise. For frd no output is given for a time-series model.
In addition, the normal subreferencing can be used.
sys = ss(mod(1,[3 4]))
to convert the noise channels e to regular input channels. These channels will
have the names e@yname.
4-173
ssdata
Purpose
4ssdata
Syntax
[A,B,C,D,K,X0] = ssdata(m)
[A,B,C,D,K,X0,dA,dB,dC,dD,dK,dX0] = ssdata(m)
Description
m is the model given as any idmodel object. A, B, C, D,K, and X0 are the matrices
in the state-space description
x ( t ) = Ax ( t ) + Bu ( t ) + Ke ( t )
x ( 0 ) = x0
y ( t ) = Cx ( t ) + Dx ( t ) + e ( t )
where x ( t ) is x ( t ) or x ( t + Ts ) depending on whether m is a continuous or
discrete-time model.
dA, dB, dC, dD, dK, and dX0 are the standard deviations of the state-space
matrices.
gives the response from the measured inputs. This is a model without a
disturbance description. Moreover
[A,B,C,D,K] = ssdata(m('n'))
4-174
ssdata
To obtain state-space descriptions that treat all input channels, both u and e
as measured inputs, first apply the conversion
m = noisecnv(m)
or
m = noisecnv(m,'norm')
where the latter case first normalizes e to v, where v has a unit covariance
matrix. See the reference page for noisecnv.
Algorithm
See Also
4-175
step
Purpose
4step
Syntax
step(m)
step(data)
step(data,'sd',sd,'PW',na,Time)
step(m,'sd',sd,Time)
step(m1,m2,...,dat1, ...,mN,Time,'sd',sd)
step(m1,'PlotStyle1',m2,'PlotStyle2',...,dat1,'PlotStylek',...,mN,
'PlotStyleN',Time,'sd',sd)
[y,t,ysd] = step(m)
mod = step(data)
Description
step can be applied both to idmodels and to iddata sets, as well as to any
mixture.
For a discrete-time idmodel m, the step response y and, when required, its
estimated standard deviation ysd, is computed using sim. When called with
output arguments, y, ysd and the time vector t are returned. When step is
called without output arguments, a plot of the step response is shown. If sd is
given a value larger than zero, a confidence region around the response is
drawn. It corresponds to the confidence of sd standard deviations. If the input
argument list contains 'fill', this region is plotted as a filled area.
The start time T1 and the end time T2 can be specified by Time= [T1 T2]. If T1
is not given, it is set to -T2/4. The negative time lags (the step is always
assumed to occur at time 0) show possible feedback effects in the data, when
the step is estimated directly from data. If Time is not specified, a default value
is used.
For an iddata set data, step(data) estimates a high order, noncausal FIR
model after first having prefiltered the data so that the input is as white as
possible. The step response of this FIR model and, when asked for, its
confidence region is then plotted. When called with an output argument, step,
in the iddata case, returns this FIR model, stored as an idarx model.The order
of the prewhitening filter can be specified as na. The default value is na = 10.
Any number and any mixture of models and data sets can be used as input
arguments. The responses are plotted with each input/output channel (as
defined by the models and data sets InputName and OutputName) as a separate
plot. Colors, linestyles, and marks can be defined by PlotStyle values, as in
4-176
step
step(m1,'b-*',m2,'y--',m3,'g')
The noise input channels in m are treated as follows: Consider a model m with
both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that m.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v:
y = Gu + HLv
cov ( v ) = I
step(m) plots the step response of the transfer function G.
step(m('n')) plots the step response of the transfer function H. (ny inputs
and ny outputs).The input channels have names e@yname, where yname is the
name of the corresponding output.
If m is a time series, that is nu = 0, step(m) plots the step response of the
transfer function H.
step(noisecnv(m)) plots the step response of the transfer function [G H]
(nu+ny inputs and ny outputs). The noise input channels have names
e@yname, where yname is the name of the corresponding output.
step(noisecnv(m,'norm')) plots the step response of the transfer function
[G HL] (nu+ny inputs and ny outputs). The noise input channels have names
v@yname, where yname is the name of the corresponding output.
Arguments
If step is called with a single idmodel m, the output argument y is a 3-D array
of dimension Nt-by-ny-by-nu. Here Nt is the length of the time vector t, ny is the
number of output channels and nu is the number of input channels. Thus
y(:,ky,ku) is the response in the ky-th output channel to a step in the ku-th
input channel. No plot is produced when output arguments are used.
ysd has the same dimensions as y and contains the standard deviations of y.
If step is called with an output argument and a single data set in the input
arguments, the output is returned as an idarx model mod containing the high
order FIR model, and its uncertainty. By calling step with mod, the responses
can be displayed and returned without having to redo the estimation.
4-177
step
Example
See Also
4-178
cra, impulse
struc
Purpose
4struc
Syntax
NN = struc(NA,NB,NK)
Description
Examples
The statement
NN = struc(1:2,1:2,4:5);
produces
NN =
1 1
1 1
1 2
1 2
2 1
2 1
2 2
4
5
4
5
4
5
5
4-179
timestamp
Purpose
4timestamp
Syntax
timestamp(obj)
ts = timestamp(obj)
Description
4-180
tfdata
Purpose
4tfdata
Syntax
[num,den] = tfdata(m)
[num,den,sdnum,sdden] = tfdata(m)
[num,den,sdnum,sdden] = tfdata(m,'v')
Description
m is a model given as any idmodel object with ny output channels and nu input
channels
num is a a cell array of dimension ny-by-nu. num{ky,ku} (note the curly
brackets) contains the numerator of the transfer function from input ku to
output ky. This numerator is a row vector whose interpretation is described
below.
If m is a SISO model, adding an extra input argument 'v' (for vector) will
return num and den as vectors rather than cell arrays.
The formats of num and den are the same as those used by the Signal Processing
Toolbox and the Control System Toolbox, both for continuous-time and
discrete-time models. See Examining Models in the Tutorial chapter and
the examples below.
The noise input channels in m are treated as follows: Consider a model m with
both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that m.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v:
y = Gu + HLv
cov ( v ) = I
4-181
tfdata
Examples
2z + 4z
H ( z ) = -------------------------------------------3
2
z + 2z + 3z + 5
which is the same as
1
2q + 4q
H ( q ) = ----------------------------------------------------------1
2
3
1 + 2q + 3q + 5q
Note that for discrete time, idpoly and polydata has a different interpretation
of the numerator vector, in case it does not have the same length as the
denominator vector. To avoid confusion, it is advised to fill out with zeros to
make numerator and denominator vectors the same length. That is done by
tfdata.
See Also
4-182
idpoly, noisecnv
view
Purpose
4view
Syntax
view(m)
view(m('n'))
view(m1,...,mN,Plottype)
view(m1,PlotStyle1,...,mN,PlotStyleN)
Description
Adding Plottype as a last argument specifies the type of plot in which view is
initialized. Plottype is any of 'impulse', 'step', 'bode', 'nyquist',
'nichols', 'sigma', or 'pzmap'. It can also be given as a cell array containing
any collection of these strings (up to 6) in which case a multiplot is shown.
view will not display confidence regions. For that use bode, nyquist, impulse,
step, and pzmap.
The noise input channels in m are treated as follows: Consider a model m with
both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
4-183
view
See Also
4-184
zpkdata
Purpose
4zpkdata
Syntax
[z,p,k] = zpkdata(m)
[z,p,k,dz,dp,dk] = zpkdata(m)
[z,p,k,dz,dp,dk] = zpkdata(m,'v')
Description
m is a model given as any idmodel object with ny output channels and nu input
channels.
z is a a cell array of dimension ny-by-nu. z{ky,ku} (note the curly brackets)
contains the zeros of the transfer function from input ku to output ky. This is a
column vector of possibly complex numbers.
variance of the real part, the 2-2 element is the variance of the imaginary part
and the 1-2 and 2-1 elements contain the covariance between the real and
imaginary parts.
dp contains the covariance matrices of the poles in the same way.
dk is a matrix containing the variances of the elements of k.
If m is a SISO model, adding an extra input argument 'v' (for vector) will
return z and p as vectors rather than cell arrays.
Note that the zeros and the poles are associated with the different channels
combinations. To obtain the so-called transmission zeros, use tzero.
4-185
zpkdata
The noise input channels in m are treated as follows: Consider a model m with
both measured input channels u (nu channels) and noise channels e (ny
channels) with covariance matrix
y = Gu + He
cov ( e ) = = LL
where L is a lower triangular matrix. Note that m.NoiseVariance = . The
model can also be described with unit variance, normalized noise source v.
y = Gu + HLv
cov ( v ) = I
Then:
zpkdata(m) returns the zeros and poles of G.
zpkdata(m('n')) returns the zeros and poles of H. (ny inputs and ny
outputs)
If m is a time series, that is nu = 0, zpkdata(m) returns the zeros and poles of
H.
zpkdata(noisecnv(m)) returns the zeros and poles of the transfer function
[G H] (nu+ny inputs and ny outputs)
zpkdata(noisecnv(m,'norm')) returns the zeros and poles of the transfer
function [G HL] (nu+ny inputs and ny outputs.
The procedure handles both models in continuous and discrete time.
Note that you cannot rely on information about zeros and poles at the origin
and at infinity for discrete time models. (This is a somewhat confusing issue
anyway.)
Algorithm
4-186
The poles and zeros are computed using ss2zp. The covariance information is
computed using Gauss s approximation formula, using the parameter
covariance matrix contained in m. When the transfer function depends on the
parameters in a complicated way, numerical differentiation is applied. The
step-sizes for the differentiation are determined in the M-file nuderst.
Index
A
adaptive noise cancelling 4-141
Advanced 4-14
AIC, the Akaike Information Criterion 4-9
Akaikes Final Prediction Error (FPE) 3-66
AR model 3-27
ARARMAX structure 3-12
ARMAX 2-23
ARMAX model 2-23
ARMAX structure 3-12
ARX 2-23
ARX model 2-21
B
basic tools 3-3
BJ 2-23
Bode diagram 2-30
Bode plot 1-10
Box-Jenkins (BJ) structure 3-12
Box-Jenkins model 2-23
Burgs method 3-28
C
communication window ident 2-2
comparisons using compare 3-53
complex-valued data 3-100
confidence interval 2-29
correlation analysis 1-4
covariance matrix 3-29
suppressing calculation 3-94
covariance method 3-28
creating models from data 2-2
cross correlation function 1-17
cross spectrum 3-15
customized plots 2-37
D
data
channels 3-22
feedback 3-68
multiple experiments 3-23
Data Board 2-2
data handling checklist 2-13
data representation 2-7
data views 1-4
delays 3-11
detrending the data 2-11
difference equation 1-7
disturbance 1-6
disturbance spectra 2-30
drift matrix 3-81
dynamic models, introduction 1-6
E
empirical transfer function estimate 3-21
estimation
parametric 3-26
estimation data 1-4
estimation method
instrumental variables 3-17
nonparametric 3-19
parametric 3-26
prediction error approach 2-19
estimation methods
direct 2-15
parametric 2-15
exporting to the MATLAB workspace 2-34
extended least squares (ELS) 3-84
I-1
Index
F
fault detection 3-85
feedback 1-15
feedback in data 3-68
focus 2-12
frequency
function 2-16
functions 3-9
plots 3-9
range 3-9
response 2-16
scales 3-9
frequency domain description 3-10
frequency response 1-10
G
Gauss-Newton direction 4-14
Gauss-Newton minimization 3-29
geometric lattice method 3-28
graphical user interface (GUI) 2-2
greybox-modeling 3-47
GUI 2-2
topics 2-35
H
Hamming window 3-21
I
idarx model object 3-39
ident window 2-35
identification method
subspace 2-26
identification process, basic steps 1-12
idfrd model object 3-20
I-2
K
Kalman gain 3-14
L
lag widow 3-16
lag window 3-21
layout 2-36
least squares 2-22
Levenberg-Marquard 4-14
LimitError 4-13
M
main ident window 2-35
maximum likelihood
criterion 3-30
method 3-17
MaxIter 4-13
Index
MaxSize 4-11
memory horizon 3-82
merge experiments 3-24
model
nonparametric 3-11
output-error 1-8
parametric 2-18
properties 1-10
set 1-4
state-space 1-8
structure 2-17
uncertainty 3-70
view functions 2-28
views 1-8
Model Board 2-2
model order 1-7
model structure 1-4
model uncertainty 2-29
model validation 1-5
model views 1-4
multi-output models
criterion 3-30
Multiple experiments 3-23
multivariable ARX model 3-39
multivariable systems 1-18
N
N4Horizon 3-33
N4Weight 3-33
na,nb,nc,nd,nf
parameter definitions 3-11
noise 1-6
noise model 1-8
noise source 1-8
noise-free simulation 1-9
noise-free simulations 3-70
O
OE 2-23
offsets 3-76
order editor 2-19
outliers
signals 1-4
output error model 2-23
output signals 1-6
Output-Error model 1-8
output-error model 2-23
state space model 3-43
P
parametric identification 1-4
parametric model 2-18
parametric model estimation 3-26
periodogram 3-22
pole 2-31
poles 1-11
poorly damped systems 1-18
prediction
error identification 2-19
error method 3-17
prediction error 3-16
preferences
in GUI 2-37
prefiltering 2-12
prefiltering signals 2-11
I-3
Index
R
recursive
identification 3-4
references list 1-21
resampling 2-12
residual analysis 1-17
residuals 1-2
robustification 4-13
S
sampling interval 1-6
SearchDirection 4-13
selecting data ranges 2-11
Sessions 2-5
shift operator 2-23
simulating data 2-13
spectral analysis 1-4
spectrum 3-9
startup identification procedure 1-14
state space model
continuous time 3-14
output-error model 3-43
stochastic 3-13
state variables 1-8
state vector 3-13
state-space
model 2-25
state-space model 3-13
state-space models 1-7
step response 1-10
structure 1-4
structure matrices 3-44
I-4
T
time delay 1-7
time domain description 3-10
time series model 3-21
time-continuous systems 3-38
Tolerance 4-13
trace 3-29
transfer function 1-8
transient response 1-10
typographical conventions (table) x
U
uncertainty
suppressing calculation 3-94
Unnormalized Gradient (UG) Approach 3-83
V
validation data 1-2
W
white noise 1-9
window sizes 3-21
working data 1-4
Working Data set 2-3
Y
Yule-Walker approach 3-28
Z
zero 2-31
Index
zeros 1-11
I-5
Index
I-6