System Identification Toolbox™

Algorithm Properties

Algorithm properties affecting estimation process for linear models

Syntax

idprops idmodel algorithm
Model.algorithm.PropertyName='PropertyValue'

Description

Algorithm is a property of the idmodel class that specifies the estimation algorithm. The idmodel subclasses are used to define various linear models, including idarx, idss, idpoly, idproc, and idgrey. These models inherit the idmodel Algorithm property.

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

Note You can use the get function or dot notation to fetch fields of Algorithm as if they were the properties of the model itself. Similarly you can use set or dot notation to set a particular field. This shortcut access is available for the fields of Algorithm only (not for deeper level struct fields such as Search or Threshold options). For example:

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

is equivalent to

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

When you create a new model by refining an existing model m, the algorithm properties of m are inherited by the new model.

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

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

You can also specify the algorithm properties (except advanced properties) as property-value pairs when creating the linear model (using idpoly, idss, etc.) or when estimating them (using pem, n4sid, armax, oe etc.).

Algorithm Properties

Criterion: Specifies criterion used during minimization. Criterion can have the following values:

'Det': Minimize , where E represents the prediction error. This is the optimal choice in a statistical sense and leads to the maximum likelihood estimates in case nothing is known about the variance of the noise. It uses the inverse of the estimated noise variance as the weighting function. This is the default criterion used for all models, except idnlgrey which uses 'Trace' by default.
'Trace': Minimize the trace of the weighted prediction error matrix trace(E'*E*W), where E is the matrix of prediction errors, with one column for each output, and W is a positive semi-definite symmetric matrix of size equal to the number of outputs. By default, W is an identity matrix of size equal to the number of model outputs (so the minimization criterion becomes trace(E'*E), or the traditional least-sum-of-squared-errors criterion. You can specify the relative weighting of prediction errors for each output using the Weighting field of the Algorithm property.

Note The difference between the two criteria is meaningful in multiple-output cases only. In single-output models, the two criteria are equivalent. Both the Det and Trace criteria are derived from a general requirement of minimizing a weighted sum of squares of prediction errors. The Det criterion can be interpreted as estimating the covariance matrix of the noise source and using the inverse of that matrix as the weighting. When using the Trace criterion, you must specify the weighting using the Weighting property.

If you want to achieve better accuracy for a particular channel in multiple-input multiple-output models, you should use Trace with weighting that favors that channel. Otherwise it is natural to use the Det criterion. When using Det, you can check cond(model.NoiseVariance) after estimation. If the matrix is ill-conditioned, it may be more robust to use the Trace criterion. You can also use compare on validation data to check whether the relative error for different channels corresponds you your needs or expectations. Use the Trace criterion if you need to modify the relative errors, and check model.NoiseVariance to determine what weighting modifications to specify.

The search method of lsqnonlin supports the Trace criterion only.

Focus: Defines how the errors e between the measured and the modeled outputs are weighed at specific frequencies during the minimization of the following loss function:

Higher weighting at specific frequencies emphasizes the requirement for a good fit at these frequencies. Focus can have the following values:
- 'Prediction': (Default) Automatically calculates the weighting function as a product of the input spectrum and the inverse of the noise model. This minimizes the one-step-ahead prediction, which typically favors fitting small time intervals (higher frequency range). From a statistical-variance point of view, this is the optimal weighting function. However, this method neglects the approximation aspects (bias) of the fit. Might not result in a stable model. Use 'Stability' when you want to ensure a stable model.
- 'Simulation': Estimates the model using the frequency weighting of the transfer function that is given by the input spectrum. Typically, this method favors the frequency range where the input spectrum has the most power. In other words, the resulting model will produce good simulations for inputs that have the same spectra as used for estimation. For models that have no disturbance model, there is no difference between 'Simulation' and 'Prediction'. In this case, with H=1, which is equivalent to A=C=D=1 for idpoly models and K = 0 for idss models.
  For models that have a disturbance model, G is first estimated with H=1, and then H is estimated using a prediction-error method with a fixed estimated transfer function . This guarantees a stable transfer function G.
- 'Stability': This weighting is the same as for 'Prediction', but the model is forced to be stable. Use only when you know the system is stable. In some cases, forcing the model to be stable can result in a bad model.
- Enter a row vector or a matrix containing frequency values that define desired passbands. For example:
```
[wl,wh]
[w1l,w1h;w2l,w2h;w3l,w3h;...]
```
  where wl and wh represent upper and lower limits of a passband. For a matrix with several rows defining frequency passbands, the algorithm uses union of frequency ranges to define the estimation passband.
- Enter any SISO linear filter in any of the following ways:
  A single-input-single-output (SISO) idmodel object.
  An ss, tf, or zpk model from the Control System Toolbox™ product.
  Using the format {A,B,C,D}, which specifies the state-space matrices of the filter.
  Using the format {numerator, denominator}, which specifies the numerator and denominator of the filter transfer function
  This calculates the weighting function as a product of the filter and the input spectrum to estimate the transfer function from input to output, G. To obtain a good model fit for a specific frequency range, you must choose the filter with a passband in this range. After estimating G, the algorithm computes the disturbance model using a prediction-error method and keeping the estimated transfer function fixed (similar to the 'Simulation' case). For a model that has no disturbance model, the estimation result is the same if you first prefilter the data using idfilt.
- For frequency-domain data only, enter a column vector of weights for 'Focus'. This vector must have the same size as length of the frequency vector of the data set, Data.Frequency. Each input and output response in the data is multiplied by the corresponding weight at that frequency.
Maxsize: MaxSize: A positive integer, specified such that the input-output data is split into segments where each contains fewer than MaxSize elements. Setting MaxSize can improve computational performance. The default value of MaxSize is 'Auto', which uses the M-file idmsize to set the value. You can edit this file to optimize computational speed on a particular computer. MaxSize does not affect the numerical properties of the estimate except when you use InitialState = 'backcast' for frequency-domain data. In this case, the frequency ranges where backcasting takes place might depend on MaxSize and affects estimates.
FixedParameter: Vector of integers containing the indices of parameters that are not estimated and remain fixed at nominal or initial values. Parameter indices refer to their position in the list, stored in the property 'ParameterVector'. You can also specify parameter names as values from the property 'PName'. To specify fixed parameters using parameter names, enter Fixedparameter as a cell array of strings. For example, to fix parameters with names 'a' and 'b', type m.FixedParameter = {'a','b','c'}. Strings can contain wildcards, such as '*' to specify the automatic completion of a string, or '?' to represent an arbitrary character. For example, to fix three parameters in a disturbance model that start with 'k', such as 'k1', 'k2','k3', use FixedParameter = {'k*'}. For structured state-space models, you can fix and free parameters by specifying structure matrices, such as As and Bs (see idss).
Weighting: Positive semi-definite symmetric matrix W to use as weighting for minimization of the trace criterion trace(E'*E*W). Weighting can be used to specify relative importance of outputs in multiple-input multiple-output models (or reliability of corresponding data) by specifying W as a diagonal matrix of non-negative values. Weighting is not useful in single-input single-output models. By default, Weighting is the identity matrix of size equal to the number of model outputs, assigning equal importance to each output during estimation.

Trace: Specifies what information displays in the MATLAB^® Command Window about the iterative search during estimation.

'Trace'='Off': Displays no information.
'Trace'='On': Displays the loss-function values for each iteration.
'Trace'='Full': Displays the same information as 'On' and also include the current parameter values and the search direction (except when the Advanced SSParameterization model property is set to 'Free' for idss models and the list of parameters can change between iterations).

LimitError: Specifies when to adjust the weight of large errors from quadratic to linear. Default value is 0. Errors larger than LimitError times the estimated standard deviation have a linear weight in the criteria. The standard deviation is estimated robustly as the median of the absolute deviations from the median and divided by 0.7. (See the section about choosing a robust norm in [2].) LimitError = 0 disables the robustification and leads to a purely quadratic criterion. When estimating with frequency-domain data, LimitError is set to zero.

Note You can estimate the model with the default value of LimitError (zero) and plot the prediction errors using pe(data.model). If the resulting plot shows occasional large values, repeat the estimation with model.Algorithm.LimitError set to a value between 1 and 2.

MaxIter: Specifies the maximum number of iterations during loss-function minimization. The iterations stop when MaxIter is reached or another stopping criterion is satisfied, such as the Tolerance. The default value of MaxIter is 20. Setting MaxIter = 0 returns the result of the startup procedure. Use EstimationInfo.Iterations to get the actual number of iterations during an estimation.
Tolerance: Specifies the minimum percentage difference (divided by 100) between the current value of the loss function and its expected improvement after the next iteration: When the percentage of expected improvement is less than Tolerance, the iterations are stopped. Default value is 0.01. The estimate of the expected loss-function improvement at the next iteration is made based on the Gauss-Newton vector computed for the current parameter value.
SearchMethod: The search method used for iterative parameter estimation. It can take one of the following values:
- 'gn': The subspace Gauss-Newton direction. Singular values of the Jacobian matrix less than GnPinvConst*eps*max(size(J))*norm(J) are discarded when computing the search direction, where J is the Jacobian matrix. The Hessian matrix is approximated by J^TJ. If there is no improvement along this direction, the gradient direction is also tried.
- 'gna': An adaptive version of subspace Gauss-Newton approach, suggested by Wills and Ninness (IFAC World congress, Prague 2005). Eigenvalues less than gamma*max(sv) of the Hessian are neglected , where sv are the singular values of the Hessian. The Gauss-Newton direction is computed in the remaining subspace. gamma has the initial value InitGnaTol (see below) and is increased by the factor LmStep each time the search fails to find a lower value of the criterion in less than 5 bisections. It is decreased by the factor 2*LmStep each time a search is successful without any bisections.
- 'lm': Uses the Levenberg-Marquardt method. 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, and grad is the gradient. d is a number that is increased until a lower value of the criterion is found.
- 'Auto': A choice among the above is made in the algorithm. This is the default choice.
- 'lsqnonlin': Uses lsqnonlin optimizer from Optimization Toolbox software. You must have Optimization Toolbox installed to use this option. This search method can only handle the Trace criterion.
Advanced: Structure that specifies advanced algorithm options and has the following fields:
- Search: Uses the following fields to specify options for the iterative search:
1. GnPinvConst: Must be a positive real value. Specifies that singular values of the Jacobian that are smaller than GnPinvConst*max(size(J)*norm(J)*eps are discarded when computing the search direction using the 'gn' method. Default value is 1e4.
2. InitGnaTol: The initial value of gamma in the gna search algorithm. See SearchMethod for description of gna. Default is 10^-4.
3. LmStep: The size of the Levenberg-Marquardt step. The next value of the search-direction length d in the Levenberg-Marquardt method is LmStep times the previous one. Default is 2.
4. StepReduction: For search directions other than the Levenberg-Marquardt direction, the step is reduced by the factor StepReduction after each iteration. Default is 2.
5. MaxBisection: The maximum number of bisections used by the line search along the search direction. Default is 25.
6. LmStartValue: The starting value of search-direction length d in the Levenberg-Marquardt method. Default is 0.001.
7. RelImprovement: The iterations are stopped if the relative improvement of the criterion is less than RelImprovement. Default is RelImprovement = 0. This property is different from Tolerance in that it uses the actual improvement of the loss function, as opposed to the expected improvement.
- Threshold: Contains fields with thresholds for several tests:
1. Sstability: Specifies the location of the rightmost pole to test the stability of continuous-time models. A model is considered stable when its rightmost pole is to the left of Sstability. Default is 0.
2. Zstability: Specifies the maximum distance of all poles from the origin to test stability of discrete-time models. A model is considered stable if all poles are within the distance Zstability from the origin. Default is 1+sqrt(eps).
- AutoInitialState: Specifies when to automatically estimate the initial state. When InitialState = 'Auto', the initial state is estimated when the ratio of the prediction-error norm with a zero initial state to the norm with an estimated initial state exceeds AutoInitialState. Default is 1.05.

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

Note These properties apply to n4sid. Since pem commonly uses n4sid to initialize a model for iterative estimation, these properties affect the results of pem too.

N4Weight: Calculates the weighting matrices used in the singular-value decomposition step of the algorithm and has three possible values:
- 'Auto': (Default) Automatically chooses between 'MOESP' and 'CVA' .
- 'MOESP': Uses the MOESP algorithm by Verhaegen.
- 'CVA': Uses the canonical variable algorithm by Larimore.
For more information about setting this property, see the n4sid reference page.
N4Horizon: Determines the forward and backward prediction horizons used by the algorithm. Enter a row vector with three elements: N4Horizon=[r sy su], where r is the maximum forward prediction horizon; that is, the algorithm uses up to r step-ahead predictors. sy is the number of past outputs, and su is the number of past inputs used for predictions. For an exact definition of these integers, see the section about subspace methods in [2], where they are called r, s1, and s2. These numbers can have a substantial influence on the quality of the resulting model and there are no simple rules for choosing them. Making 'N4Horizon' a k-by-3 matrix means that the algorithm tries each row of 'N4Horizon' and selects the value that gives the best (prediction) fit to the data. Choosing the best row is not available when you also specify to select the best model order. When you specify one column in 'N4Horizon', the interpretation is r=sy=su. The default choice is 'N4Horizon' = 'Auto', which uses an Akaike Information Criterion (AIC) for the selection of sy and su.

Note For algorithm properties of nonlinear models, see the reference pages for idnlarx, idnlhw, and idnlgrey.

References

[1] Dennis, J.E., Jr., and R.B. Schnabel, Numerical Methods for Unconstrained Optimization and Nonlinear Equations, Prentice Hall, Englewood Cliffs, N.J., 1983. See the chapter about iterative minimization.

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

Algorithm Properties

Syntax

Description

Algorithm Properties

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

References

See Also