n4sid
Estimate statespace model using subspace method with timedomain or frequencydomain data
Description
estimates a discretetime statespace model sys
= n4sid(data
,nx
)sys
of order
nx
using data
, which can be timedomain or
frequencydomain data. sys
is a model of the following form:
$$\begin{array}{l}\dot{x}(t)=Ax(t)+Bu(t)+Ke(t)\\ y(t)=Cx(t)+Du(t)+e(t)\end{array}$$
A, B, C, D,
and K are statespace matrices.
u(t) is the input,
y(t) is the output,
e(t) is the disturbance, and
x(t) is the vector of nx
states.
All entries of A, B, C, and
K are free estimable parameters by default. For dynamic systems,
D is fixed to zero by default, meaning that the system has no
feedthrough. For static systems (nx = 0
), D is an
estimable parameter by default.
incorporates additional options specified by one or more namevalue pair arguments. For
example, to estimate a continuoustime model, specify the sample time
sys
= n4sid(data
,nx
,Name,Value
)'Ts'
as 0
. Use the 'Form'
,
'Feedthrough'
, and 'DisturbanceModel'
namevalue
pair arguments to modify the default behavior of the A,
B, C, D, and K
matrices.
Examples
StateSpace Model
Estimate a statespace model and compare its response with the measured output.
Load the inputoutput data z1
, which is stored in an iddata
object.
load iddata1 z1
Estimate a fourthorder statespace model.
nx = 4; sys = n4sid(z1,nx)
sys = Discretetime identified statespace model: x(t+Ts) = A x(t) + B u(t) + K e(t) y(t) = C x(t) + D u(t) + e(t) A = x1 x2 x3 x4 x1 0.8392 0.3129 0.02105 0.03743 x2 0.4768 0.6671 0.1428 0.003757 x3 0.01951 0.08374 0.09761 1.046 x4 0.003885 0.02914 0.8796 0.03171 B = u1 x1 0.02635 x2 0.03301 x3 7.256e05 x4 0.0005861 C = x1 x2 x3 x4 y1 69.08 26.64 2.237 0.5601 D = u1 y1 0 K = y1 x1 0.003282 x2 0.009339 x3 0.003232 x4 0.003809 Sample time: 0.1 seconds Parameterization: FREE form (all coefficients in A, B, C free). Feedthrough: none Disturbance component: estimate Number of free coefficients: 28 Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using N4SID on time domain data "z1". Fit to estimation data: 76.33% (prediction focus) FPE: 1.21, MSE: 1.087
Compare the simulated model response with the measured output.
compare(z1,sys)
The plot shows that the fit percentage between the simulated model and the estimation data is greater than 70%.
You can view more information about the estimation by exploring the idss
property sys.Report
.
sys.Report
ans = Status: 'Estimated using N4SID with prediction focus' Method: 'N4SID' InitialState: 'estimate' N4Weight: 'CVA' N4Horizon: [6 10 10] Fit: [1x1 struct] Parameters: [1x1 struct] OptionsUsed: [1x1 idoptions.n4sid] RandState: [1x1 struct] DataUsed: [1x1 struct]
For example, find out more information about the estimated initial state.
sys.Report.Parameters.X0
ans = 4×1
0.0085
0.0052
0.0193
0.0282
Determine Optimal Estimated Model Order
Load the inputoutput data z1
, which is stored in an
iddata
object.
load iddata1 z1
Determine the optimal model order by specifying argument nx
as a
range from 1 to 10.
nx = 1:10; sys = n4sid(z1,nx);
An automatically generated plot shows the Hankel singular values for models of the
orders specified by nx
.
States with relatively small Hankel singular values can be safely discarded. The
suggested default order choice is 2
.
Select the model order in the Chosen Order list and click Apply.
Specify Estimation Options
Load estimation data.
load iddata2 z2
Specify estimation options. Set the weighting scheme 'N4Weight'
to 'SSARX'
and estimationstatus display option 'Display'
to 'on'
.
opt = n4sidOptions('N4Weight','SSARX','Display','on')
Option set for the n4sid command: InitialState: 'estimate' N4Weight: 'SSARX' N4Horizon: 'auto' Display: 'on' InputOffset: [] OutputOffset: [] EstimateCovariance: 1 OutputWeight: [] Focus: 'prediction' WeightingFilter: [] EnforceStability: 0 Advanced: [1x1 struct]
Estimate a thirdorder statespace model using the updated option set.
nx = 3; sys = n4sid(z2,nx,opt);
Modify Form, Feedthrough, and DisturbanceModel Matrices
Modify the canonical form of the A, B, and C matrices, include a feedthrough term in the D matrix, and eliminate disturbancemodel estimation in the K matrix.
Load inputoutput data and estimate a fourthorder system using the n4sid
default options.
load iddata1 z1 sys1 = n4sid(z1,4);
Specify the modal form and compare the A
matrix with the default A
matrix.
sys2 = n4sid(z1,4,'Form','modal'); A1 = sys1.A
A1 = 4×4
0.8392 0.3129 0.0211 0.0374
0.4768 0.6671 0.1428 0.0038
0.0195 0.0837 0.0976 1.0462
0.0039 0.0291 0.8796 0.0317
A2 = sys2.A
A2 = 4×4
0.7554 0.3779 0 0
0.3779 0.7554 0 0
0 0 0.0669 0.9542
0 0 0.9542 0.0669
Include a feedthrough term and compare the D
matrices.
sys3 = n4sid(z1,4,'Feedthrough',1);
D1 = sys1.D
D1 = 0
D3 = sys3.D
D3 = 0.0487
Eliminate disturbance modeling and compare the K
matrices.
sys4 = n4sid(z1,4,'DisturbanceModel','none'); K1 = sys1.K
K1 = 4×1
0.0033
0.0093
0.0032
0.0038
K4 = sys4.K
K4 = 4×1
0
0
0
0
ContinuousTime CanonicalForm Model
Estimate a continuoustime canonicalform model.
Load estimation data.
load iddata1 z1
Estimate the model. Set Ts
to 0
to specify a continuous model.
nx = 2; sys = n4sid(z1,nx,'Ts',0,'Form','canonical');
sys
is a secondorder continuoustime statespace model in the canonical form.
Estimate StateSpace Model from ClosedLoop Data
Estimate a statespace model from closedloop data using the subspace algorithm SSARX
. This algorithm is better at capturing feedback effects than other weighting algorithms.
Generate closedloop estimation data for a secondorder 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(k1) = K*y(k2) + w(k); u(k1) = K*y(k1) + w(k); z(k) = 1.5*z(k1)  0.7*z(k2) + u(k1) + 0.5*u(k2); y(k) = z(k) + 0.8*v(k); end dat = iddata(y, u, 1);
Specify the weighting scheme 'N4weight'
used by the N4SID algorithm. Create two option sets. For one option set, set 'N4weight'
to 'CVA'
. For the other option set, set the 'N4weight'
to 'SSARX'
.
optCVA = n4sidOptions('N4weight','CVA'); optSSARX = n4sidOptions('N4weight','SSARX');
Estimate statespace models using the option 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);
As the plot shows, the model estimated using the SSARX algorithm produces a better fit than the model estimated using the CVA algorithm.
Input Arguments
data
— Estimation data
iddata
object  frd
object  idfrd
object
Estimation data, specified as an iddata
object, an frd
object, or an idfrd
object.
For timedomain estimation, data
must be an iddata
object containing the input and output signal values.
For frequencydomain estimation, data
can be one of the following:
Estimation data must be uniformly sampled. By default, the software sets the sample time of the model to the sample time of the estimation data.
For multiexperiment data, the sample times and intersample behavior of all the experiments must match.
The domain of your data determines the type of model you can estimate.
Timedomain or discretetime frequencydomain data — Continuoustime and discretetime models
Continuoustime frequencydomain data — Continuoustime models only
nx
— Order of estimated model
1:10
(default)  positive integer scalar  positive integer vector  best
 0
Order of the estimated model, specified as a nonnegative integer or as a vector containing a range of positive integers.
If you already know what order you want your estimated model to have, specify
nx
as a scalar.If you want to compare a range of potential orders to choose the most effective order for your estimated model, specify that range for
nx
.n4sid
creates a Hankel singularvalue plot that shows the relative energy contributions of each state in the system. States with relatively small Hankel singular values contribute little to the accuracy of the model and can be discarded with little impact. The index of the highest state you retain is the model order. The plot window includes a suggestion for the order to use. You can accept this suggestion or enter a different order. For an example, see Determine Optimal Estimated Model Order.If you do not specify
nx
, or if you specifynx
asbest
, the software automatically choosesnx
from the range1:10
.If you are identifying a static system, set
nx
to0
.
opt
— Estimation options
n4sidOptions
option set
Estimation options, specified as an n4sidOptions
option set. Options specified by opt
include:
Estimation objective
Handling of initial conditions
Subspace algorithmrelated choices
For examples showing how to specify options, see Specify Estimation Options and ContinuousTime CanonicalForm Model.
NameValue Arguments
Specify optional pairs of arguments as
Name1=Value1,...,NameN=ValueN
, where Name
is
the argument name and Value
is the corresponding value.
Namevalue arguments must appear after other arguments, but the order of the
pairs does not matter.
Before R2021a, use commas to separate each name and value, and enclose
Name
in quotes.
Example:
sys = n4sid(data,nx,'Form','modal')
Ts
— Sample time of estimated model
sample time of data (data.Ts
) (default)  0
 positive scalar
Sample time of the estimated model, specified as the commaseparated pair
consisting of 'Ts'
and either 0
or a positive scalar.
For continuoustime models, specify
'Ts'
as0
. In the frequency domain, using continuoustime frequencydomain data results in a continuoustime model.For discretetime models, the software sets
'Ts'
to the sample time of the data in the units stored in theTimeUnit
property.
InputDelay
— Input delays
0
(default)  scalar  vector
Input delay for each input channel, specified as the commaseparated pair
consisting of 'InputDelay'
and a numeric vector.
For continuoustime models, specify
'InputDelay'
in the time units stored in theTimeUnit
property.For discretetime models, specify
'InputDelay'
in integer multiples of the sample timeTs
. For example, setting'InputDelay'
to3
specifies a delay of three sampling periods.
For a system with N_{u} inputs, set
InputDelay
to an
N_{u}by1 vector. Each entry of this vector is
a numerical value that represents the input delay for the corresponding input
channel.
To apply the same delay to all channels, specify 'InputDelay'
as a scalar.
Form
— Type of canonical form
'free'
(default)  'modal'
 'companion'
 'canonical'
Type of canonical form of sys
, specified as the
commaseparated pair consisting of 'Form'
and one of the following values:
'free'
— Treat all entries of the matrices A, B, C, D, and K as free.'modal'
— Obtainsys
in modal form.'companion'
— Obtainsys
in companion form.'canonical'
— Obtainsys
in the observability canonical form.
For definitions of the canonical forms, see Canonical StateSpace Realizations.
For more information about using these forms for identification, see Estimate StateSpace Models with Canonical Parameterization.
For an example, see Modify Form, Feedthrough, and DisturbanceModel Matrices.
Feedthrough
— Direct feedthrough from input to output
0
(default)  1
 logical vector
Direct feedthrough from input to output, specified as the commaseparated pair
consisting of 'Feedthrough'
and a logical vector of length
N_{u}, where
N_{u} is the number of inputs. If you
specify 'Feedthrough'
as a logical scalar, that value is applied to
all the inputs. For static systems, the software always assumes
'Feedthrough'
is 1
.
For an example, see Modify Form, Feedthrough, and DisturbanceModel Matrices.
DisturbanceModel
— Option to estimate timedomain noise component parameters
'estimate'
 'none'
Option to estimate timedomain noise component parameters in the K matrix,
specified as the commaseparated pair consisting of
'DisturbanceModel'
and one of the following values:
'estimate'
— Estimate the noise component. The K matrix is treated as a free parameter. For timedomain data,'estimate'
is the default.'none'
— Do not estimate the noise component. The elements of the K matrix are fixed at zero. For frequencydomain data,'none'
is the default and the only acceptable value.
For an example, see Modify Form, Feedthrough, and DisturbanceModel Matrices.
Output Arguments
sys
— Identified statespace model
idss
model
Identified statespace model, returned as an 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 Field  Description  

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:
This field is especially useful when the
 
N4Weight  Weighting scheme used for singularvalue decomposition by the N4SID algorithm, returned as one of the following values:
This option is especially useful when the
 
N4Horizon  Forward and backward prediction horizons used by the N4SID algorithm,
returned as a row vector with three elements
 
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:
 
Parameters  Estimated values of model parameters.  
OptionsUsed  Option set used for estimation. If you did not configure any custom
options,  
RandState  State of the random number stream at the start of estimation. Empty,
 
DataUsed  Attributes of the data used for estimation, returned as a structure with the following fields.

For more information on using Report
, see Estimation Report.
x0
— Initial states computed during estimation
column vector  matrix
Initial states computed during the estimation, returned as an array containing a column vector corresponding to each experiment.
This array is also stored in the Parameters
field of the model
Report
property.
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.
[5] McKelvey, T., H. Akcay, and L. Ljung. "Subspacebased multivariable system identification from frequency response data." IEEE Transactions on Automatic Control, 1996, Vol. 41, pp. 960–979.
Version History
Open Example
You have a modified version of this example. Do you want to open this example with your edits?
MATLAB Command
You clicked a link that corresponds to this MATLAB command:
Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
 América Latina (Español)
 Canada (English)
 United States (English)
Europe
 Belgium (English)
 Denmark (English)
 Deutschland (Deutsch)
 España (Español)
 Finland (English)
 France (Français)
 Ireland (English)
 Italia (Italiano)
 Luxembourg (English)
 Netherlands (English)
 Norway (English)
 Österreich (Deutsch)
 Portugal (English)
 Sweden (English)
 Switzerland
 United Kingdom (English)