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 signal`OutputData` — Fourier transform of the output signal`Domain` — ‘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 `[]`.

For more information on using `Report`, see Estimation Report.

`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

`load iddata2 z2`

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.

`load iddata1 z1`

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 right-most 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

`$A=\left[\begin{array}{cccccc}0& 0& ..& ..& 0& -{\alpha }_{n}\\ 1& 0& 0& ..& 0& -{\alpha }_{n}-1\\ 0& 1& 0& .& :& :\\ :& 0& .& .& :& :\\ 0& .& .& 1& 0& -{\alpha }_{2}\\ 0& ..& ..& 0& 1& -{\alpha }_{1}\end{array}\right]$`

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.