Learn more about System Identification Toolbox

pem - Estimate model parameters using iterative prediction-error minimization method

Syntax

m = pem(data)
m = pem(data,mi)
m = pem(data,mi,'Property1',Value1,...,'PropertyN',ValueN)
m = pem(data,orders)
m = pem(data,'P1D')
m = pem(data,'nx',ssorder)
m = pem(data,'na',na,'nb',nb,'nc',nc,'nd',nd,'nf',nf,'nk',nk)
m = pem(data,orders,'Property1',Value1,...,'PropertyN',ValueN)

How to Use

If you are using the System Identification Tool GUI, you can specify PEM for low-order continuous-time process models, linear state-space, and polynomial models. If you are working in the MATLAB Command Window, you can use the pem command to both construct and estimate these linear models and to also estimate linear and nonlinear grey-box models.

Alternatively, you can use PEM to try to refine initial parameter estimates for all linear and nonlinear parametric models. For more information about refining initial model estimates, see Refining Linear Parametric Models.

Description

pem is the basic estimation command in the toolbox and covers a variety of situations.

data is always an iddata object that contains the input/output data. Both time-domain and frequency-domain signals are supported. data can also be an frd or idfrd frequency-response data object. Estimation of noise models (K in state-space models and A, C, and D in polynomial models) is not supported for frequency-domain data.

With Initial Model

mi is any idmodel or idnlmodel object. It could be a result of another estimation routine, or constructed and modified by the constructors (idarx, idpoly, idss, idgrey, idproc) and set. The properties of mi can also be changed by any property name/property value pairs in pem as indicated in the syntax.

m is then returned as the best fitting model in the model structure defined by mi. The iterative search is initialized at the parameters of the initial/nominal model mi. m will be of the same class as mi.

Black-Box State-Space Models

With m = pem(data,n), where n is a positive integer, or m = pem(data,'nx',n), a state-space model of order n is estimated.

If data is continuous-time (frequency-domain) data, a corresponding continuous-time state space model is estimated.

The default is that it is estimated in a 'Free' parameterization that can be further modified by the properties 'nk', 'DisturbanceModel', and 'InitialState' (see the corresponding reference pages for idss and n4sid). The model is initialized by n4sid and then further adjusted by optimizing the prediction error fit.

You can choose among several different orders by

m = pem(data,'nx',[n1,n2,...nN])

and you are then prompted for the "best" order. By

m = pem(data,'best')

an automatic choice of order among 1:10 is made.

m = pem(data)

is short for m = pem(data,'best'). To work with other delays, use, for example, m = pem(data,'best','nk',[0,...0]).

In this case m is returned as an idss model.

Estimating the D, K, and X0 Matrices

Whether the D matrix is estimated or not is governed by the property nk, which is further described below. The default is that D is not estimated. By setting the kth entry of nk to 0, the kth column of D (corresponding to the kth input) is estimated. To estimate a full D matrix, let nk = zeros(1,nu), as in

m = pem(data,order,'nk',[0 .. 0])

This holds for both discrete- and continuous-time models.

For frequency-domain data, K is always fixed to 0. For time-domain data, K is estimated by default. To fix K to 0 in this case, use

m = pem(data,order,'DisturbanceModel','none')

Similarily, X0 is estimated if 'InitialState' is set to 'Estimate', and fixed to 0 if 'InitialState' is set to 'Zero'.

Black-Box Multiple-Input-Single-Output Models

The function pem also handles the general multiple-input-single-output structure

The orders of this general model are given either as

orders = [na nb nc nd nf nk]

or with (...'na',na,'nb',nb,...) as shown in the syntax. Here na, nb, nc, nd, and nf are the orders of the model, and nk is the delay(s). For multiple-input systems, nb, nf, and nk are row vectors giving the orders and delays of each input. (See What Are Black-Box Polynomial Models? in the User's Guide for a definition of the orders.) When the orders are specified with separate entries, those not given are taken as zero.

For frequency-domain data, only estimation of B and F is supported. It is simpler to use oe in that case.

In this case, m is returned as an idpoly object.

Continuous-Time Process Models

Entering for the initial model an acronym for a process model, as in

m = pem(data,'P2UI')

will estimate a continuous-time process model of the indicated type. See the reference page for idproc for details of possible model types and associated property name/property value pairs.

In this case, m is returned as an idproc model.

Properties

In all cases the algorithm is affected by the properties (see Algorithm Properties for details):

Focus can be set to 'Prediction' (default), 'Simulation', or a passband range.
MaxIter and Tolerance govern the stopping criteria for the iterative search.
LimitError deals with how the criterion can be made less sensitive to outliers and bad data.
MaxSize determines the largest matrix ever formed by the algorithm. The algorithm goes into for loops to avoid larger matrices, which can be more efficient than using virtual memory.
Display, with possible values 'Off', 'On', and 'Full', governs the information sent to the MATLAB Command Window.

For black-box state-space models, 'N4Weight' and 'N4Horizon' will also affect the result, since these models are initialized with an n4sid estimate. See the reference page for n4sid.

Typical idmodel properties are (see idmodel properties for more details):

Ts is the sampling interval. Set 'Ts'= 0 to obtain a continuous-time state-space model. For discrete-time models, 'Ts' is automatically set to the 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 using d2c.
nk is the time delays from the inputs (not applicable to structured state-space models). Time delays specified by 'nk' will be included in the model.
DisturbanceModel determines the parameterization of K for free and canonical state-space parameterizations, as well as for idgrey models. It also determines whether a noise model should be included for idproc models.
InitialState: The initial state can have a substantial influence on the estimation result for systems 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'.

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')

Building an ARMAX model for the response to output 2,

ma = pem(z(:,2,:),'na',2,'nb',[2 3 1],'nc',2,'nk',[1 2 0])

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.

PEM uses optimization to minimize the cost function, defined as follows for scalar outputs:

where e(t) is the difference between the measured output and the predicted output of the model. For a linear model, this error is defined by the following equation:

e(t) is a vector and the cost function is a scalar value. The subscript N indicates that the cost function is a function of the number of data samples and becomes more accurate for larger values of N. For multiple-output models, the previous equation is more complex.

For black-box models, PEM estimates an initial model and then varies the parameter values along a specific direction to decrease the cost function. As with any nonlinear optimization algorithm, there is a chance that the model might find a local minimum that is not accurate for a specific system.