Documentation

n4sid

Estimate state-space model using subspace method

Syntax

sys = n4sid(data,nx)
sys = n4sid(data,nx,Name,Value)
sys = n4sid(___,opt)
[sys,x0] = n4sid(___)

Description

sys = n4sid(data,nx) estimates an nx order state-space model, sys, using measured input-output data, data.

sys is an idss model representing the system:

$\begin{array}{l}\stackrel{˙}{x}\left(t\right)=Ax\left(t\right)+Bu\left(t\right)+Ke\left(t\right)\\ y\left(t\right)=Cx\left(t\right)+Du\left(t\right)+e\left(t\right)\end{array}$

A,B,C, and D are state-space matrices. K is the disturbance matrix. u(t) is the input, y(t) is the output, x(t) is the vector of nx states and e(t) is the disturbance.

All the entries of the A, B, C, and K matrices are free estimation parameters by default. D is fixed to zero by default, meaning that there is no feedthrough, except for static systems (nx=0).

sys = n4sid(data,nx,Name,Value) specifies additional attributes of the state-space structure using one or more Name,Value pair arguments. Use the Form, Feedthrough, and DisturbanceModel name-value pair arguments to modify the default behavior of the A, B, C, D, and K matrices.

sys = n4sid(___,opt) specifies estimation options, opt, that configure the initial states, estimation objective, and subspace algorithm related choices to be used for estimation.

[sys,x0] = n4sid(___) also returns the estimated initial state.

Input Arguments

collapse all

 data Estimation data. For time domain estimation, data is an iddata object containing the input and output signal values. For frequency domain estimation, data can be one of the following: Recorded frequency response data (frd or idfrd)iddata object with its properties specified as follows:InputData — Fourier transform of the input signalOutputData — Fourier transform of the output signalDomain — ‘Frequency’ For multiexperiment data, the sample times and intersample behavior of all the experiments must match. You can only estimate continuous-time models using continuous-time frequency domain data. You can estimate both continuous-time and discrete-time models (of sample time matching that of data) using time-domain data and discrete-time frequency domain data. nx Order of estimated model. Specify nx as a positive integer. nx may be a scalar or a vector. If nx is a vector, then n4sid creates a plot which you can use to choose a suitable model order. The plot shows the Hankel singular values for models of different orders. States with relatively small Hankel singular values can be safely discarded. A default choice is suggested in the plot. You can also specify nx as 'best', in which case the optimal order is automatically chosen from nx = 1,..,10. opt Estimation options. opt is an options set, created using n4sidOptions, which specifies options including: Estimation objectiveHandling of initial conditionsSubspace algorithm related choices

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Sample time, specified as a positive scalar. For continuous-time models, use Ts = 0. For discrete-time models, specify Ts as a positive scalar whose value is equal to that of the data sample time.

Type of canonical form of sys, specified as one of the following values:

• 'modal' — Obtain sys in modal form.

• 'companion' — Obtain sys in companion form.

• 'free' — All entries of the A, B, and C matrices are estimated.

• 'canonical' — Obtain sys in observability canonical form [1].

Use the Form, Feedthrough, and DisturbanceModel name-value pair arguments to modify the default behavior of the A, B, C, D, and K matrices.

Logical specifying direct feedthrough from input to output, specified as a logical vector of length Nu, where Nu is the number of inputs.

If Feedthrough is specified as a logical scalar, this value is applied to all the inputs. If the model has no states, then Feedthrough is true(1,Nu).

Specify estimation of the noise component (K matrix), specified as one of the following values:

• 'none' — Noise component is not estimated. The value of the K matrix, is fixed to zero value.

• 'estimate' — The K matrix is treated as a free parameter.

DisturbanceModel must be 'none' when using frequency domain data.

Input delay for each input channel, specified as a numeric vector. For continuous-time systems, specify input delays in the time unit stored in the TimeUnit property. For discrete-time systems, specify input delays in integer multiples of the sample time Ts. For example, InputDelay = 3 means a delay of three sampling periods.

For a system with Nu inputs, set InputDelay to an Nu-by-1 vector. Each entry of this vector is a numerical value that represents the input delay for the corresponding input channel.

You can also set InputDelay to a scalar value to apply the same delay to all channels.

Output Arguments

sys

Identified state-space model, returned as a idss model. This model is created using the specified model orders, delays, and estimation options.

Information about the estimation results and options used is stored in the Report property of the model. Report has the following fields:

Report FieldDescription
Status

Summary of the model status, which indicates whether the model was created by construction or obtained by estimation.

Method

Estimation command used.

InitialState

How initial states were handled during estimation, returned as one of the following values:

• 'zero' — The initial state is set to zero.

• 'estimate' — The initial state is treated as an independent estimation parameter.

This field is especially useful when the InitialState option in the estimation option set is 'auto'.

N4Weight

Weighting scheme used for singular-value decomposition by the N4SID algorithm, returned as one of the following values:

• 'MOESP' — Uses the MOESP algorithm.

• 'CVA' — Uses the Canonical Variable Algorithm.

• 'SSARX' — A subspace identification method that uses an ARX estimation based algorithm to compute the weighting.

This option is especially useful when the N4Weight option in the estimation option set is 'auto'.

N4Horizon

Forward and backward prediction horizons used by the N4SID algorithm, returned as a row vector with three elements —  [r sy su], where r is the maximum forward prediction horizon. sy is the number of past outputs, and su is the number of past inputs that are used for the predictions.

Fit

Quantitative assessment of the estimation, returned as a structure. See Loss Function and Model Quality Metrics for more information on these quality metrics. The structure has the following fields:

FieldDescription
FitPercent

Normalized root mean squared error (NRMSE) measure of how well the response of the model fits the estimation data, expressed as a percentage.

LossFcn

Value of the loss function when the estimation completes.

MSE

Mean squared error (MSE) measure of how well the response of the model fits the estimation data.

FPE

Final prediction error for the model.

AIC

Raw Akaike Information Criteria (AIC) measure of model quality.

AICc

Small sample-size corrected AIC.

nAIC

Normalized AIC.

BIC

Bayesian Information Criteria (BIC).

Parameters

Estimated values of model parameters.

OptionsUsed

Option set used for estimation. If no custom options were configured, this is a set of default options. See n4sidOptions for more information.

RandState

State of the random number stream at the start of estimation. Empty, [], if randomization was not used during estimation. For more information, see rng in the MATLAB® documentation.

DataUsed

Attributes of the data used for estimation, returned as a structure with the following fields:

FieldDescription
Name

Name of the data set.

Type

Data type.

Length

Number of data samples.

Ts

Sample time.

InterSample

Input intersample behavior, returned as one of the following values:

• 'zoh' — Zero-order hold maintains a piecewise-constant input signal between samples.

• 'foh' — First-order hold maintains a piecewise-linear input signal between samples.

• 'bl' — Band-limited behavior specifies that the continuous-time input signal has zero power above the Nyquist frequency.

InputOffset

Offset removed from time-domain input data during estimation. For nonlinear models, it is [].

OutputOffset

Offset removed from time-domain output data during estimation. For nonlinear models, it is [].

x0

Initial states computed during the estimator of sys.

If data contains multiple experiments, then x0 is an array with each column corresponding to an experiment.

Examples

collapse all

Specify the estimation options.

opt = n4sidOptions('Focus','simulation','Display','on');

Estimate the model.

nx = 3;
sys = n4sid(z2,nx,opt);

sys is a third-order, state-space model.

Estimate a state-space model from closed-loop data using the subspace algorithm SSARX. This algorithm is better at capturing feedback effects than other weighting algorithms.

Generate closed-loop estimation data for a second-order system corrupted by white noise.

N = 1000;
K = 0.5;
rng('default');
w = randn(N,1);
z = zeros(N,1);
u = zeros(N,1);
y = zeros(N,1);
e = randn(N,1);
v = filter([1 0.5],[1 1.5 0.7],e);
for k = 3:N
u(k-1) = -K*y(k-2) + w(k);
u(k-1) = -K*y(k-1) + w(k);
z(k) = 1.5*z(k-1) - 0.7*z(k-2) + u(k-1) + 0.5*u(k-2);
y(k) = z(k) + 0.8*v(k);
end
dat = iddata(y, u, 1);

Specify the weighting scheme used by the N4SID algorithm. In one options set, specify the algorithm as CVA and in the other, specify as SSARX.

optCVA = n4sidOptions('N4weight','CVA');
optSSARX = n4sidOptions('N4weight','SSARX');

Estimate state-space models using the options sets.

sysCVA = n4sid(dat,2,optCVA);
sysSSARX = n4sid(dat,2,optSSARX);

Compare the fit of the two models with the estimation data.

compare(dat,sysCVA,sysSSARX);

From the plot, you see that the model estimated using the SSARX algorithm produces a better fit than the CVA algorithm.

Estimate a continuous-time, canonical-form model.

Specify the estimation options.

opt = n4sidOptions('Focus','simulation','Display','on');

Estimate the model.

nx = 2;
sys = n4sid(z1,nx,'Ts',0,'Form','canonical',opt);

sys is a second-order, continuous-time, state-space model in the canonical form.

collapse all

Modal Form

In modal form, A is a block-diagonal matrix. The block size is typically 1-by-1 for real eigenvalues and 2-by-2 for complex eigenvalues. However, if there are repeated eigenvalues or clusters of nearby eigenvalues, the block size can be larger.

For example, for a system with eigenvalues $\left({\lambda }_{1},\sigma ±j\omega ,{\lambda }_{2}\right)$, the modal A matrix is of the form

$\left[\begin{array}{cccc}{\lambda }_{1}& 0& 0& 0\\ 0& \sigma & \omega & 0\\ 0& -\omega & \sigma & 0\\ 0& 0& 0& {\lambda }_{2}\end{array}\right]$

Companion Form

In the companion realization, the characteristic polynomial of the system appears explicitly in the rightmost column of the A matrix.

For a system with characteristic polynomial

$P\left(s\right)={s}^{n}+{\alpha }_{1}{s}^{n-1}+\dots +{\alpha }_{n-1}s+{\alpha }_{n}$

the corresponding companion A matrix is

The companion transformation requires that the system be controllable from the first input. The companion form is poorly conditioned for most state-space computations; avoid using it when possible.

References

[1] Ljung, L. System Identification: Theory for the User, Appendix 4A, Second Edition, pp. 132–134. Upper Saddle River, NJ: Prentice Hall PTR, 1999.

[2] van Overschee, P., and B. De Moor. Subspace Identification of Linear Systems: Theory, Implementation, Applications. Springer Publishing: 1996.

[3] Verhaegen, M. "Identification of the deterministic part of MIMO state space models." Automatica, 1994, Vol. 30, pp. 61—74.

[4] Larimore, W.E. "Canonical variate analysis in identification, filtering and adaptive control." Proceedings of the 29th IEEE Conference on Decision and Control, 1990, pp. 596–604.