ssest

Estimate state-space model using time or frequency domain data

Syntax

sys = ssest(data,nx)

sys = ssest(data,nx,Name,Value)

sys = ssest(___,opt)

sys = ssest(data,init_sys)

sys = ssest(data,init_sys,opt)

[sys,x0] = ssest(___)

Description

sys = ssest(data,nx) estimates a state-space model, sys, using time- or frequency-domain data, data. sys is a state-space model of order nx and represents:

$$\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 state-space 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 the entries of A, B, C, and K are free estimable parameters by default. D is fixed to zero by default, meaning that there is no feedthrough, except for static systems (nx=0).

sys = ssest(data,nx,Name,Value) estimates the model using the additional options specified by 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 = ssest(___,opt) estimates the model using an option set, opt, that specifies options such as estimation objective, handling of initial conditions and numerical search method used for estimation.

sys = ssest(data,init_sys) estimates a state-space model using the linear system init_sys to configure the initial parameterization.

sys = ssest(data,init_sys,opt) estimates the model using an option set, opt.

[sys,x0] = ssest(___) returns the value of initial states computed during estimation.

Input Arguments

data	Estimation data. For time-domain estimation, data must be 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'
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 ssest 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.
opt	Estimation options. opt is an options set, created using ssestOptions, that specifies options including: Estimation objective Handling of initial conditions Numerical search method used for estimation If opt is not specified and init_sys is a previously estimated idss model, the options from init_sys.Report.OptionsUsed are used.
init_sys	Linear system that configures the initial parameterization of sys. You obtain init_sys by either performing an estimation using measured data or by direct construction. If init_sys is an idss model, ssest uses the parameter values of init_sys as the initial guess for estimating sys. For information on how to specify idss, see Estimate State-Space Models with Structured Parameterization. Constraints on the parameters of init_sys, such as fixed coefficients and minimum/maximum bounds are honored in estimating sys. Use the Structure property of init_sys to configure initial guesses and constraints for the A, B, C, D, and K matrices. For example: To specify an initial guess for the A matrix of init_sys, set init_sys.Structure.A.Value as the initial guess. To specify constraints for the B matrix of init_sys: Set init_sys.Structure.B.Minimum to the minimum B matrix value Set init_sys.Structure.B.Maximum to the maximum B matrix value Set init_sys.Structure.B.Free to indicate if entries of the B matrix are free parameters for estimation If init_sys is not a state-space (idss) model, the software first converts init_sys to an idss model. ssest uses the parameters of the resulting model as the initial guess for estimation. If opt is not specified, and init_sys was obtained by estimation, then the estimation options from init_sys.Report.OptionsUsed are used.

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.

'Ts' — Sample time

0 (continuous-time) (default) | sample time of data (data.Ts) | positive scalar

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 the data sample time.

'InputDelay' — Input delays

0 (default) | scalar | vector

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.

'Form' — Type of canonical form

'free' (default) | 'modal' | 'companion' | 'canonical'

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 treated as free.

• 'canonical' — Obtain sys in the 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.

For more information, see Estimate State-Space Models with Canonical Parameterization.

'Feedthrough' — Direct feedthrough from input to output

0 (default) | 1 | logical vector

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, it is applied to all the inputs.

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

'DisturbanceModel' — Specify whether to estimate the K matrix

'estimate' (default) | 'none'

Specify whether to estimate the K matrix which specifies the noise component, 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.

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

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 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: 'zero' — The initial state is set to zero. 'estimate' — The initial state is treated as an independent estimation parameter. 'backcast' — The initial state is estimated using the best least squares fit. Column vector of length Nx, where Nx is the number of states. For multi-experiment data, a matrix with Ne columns, where Ne is the number of experiments. Parametric initial condition object (x0obj) created using idpar. Only for discrete-time state-space models. 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: Field Description 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 ssestOptions 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. Structure with the following fields: Field Description Name Name of the data set. Type Data type. For ARX models, this is set to 'Time domain data'. Length Number of data samples. Ts Sample time. This is equivalent to Data.Ts. InterSample Input intersample behavior. 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. The value of Intersample has no effect on estimation results for discrete-time models. InputOffset Offset removed from time-domain input data during estimation. OutputOffset Offset removed from time-domain output data during estimation. Termination Termination conditions for the iterative search used for prediction error minimization. Structure with the following fields: Field Description WhyStop Reason for terminating the numerical search. Iterations Number of search iterations performed by the estimation algorithm. FirstOrderOptimality $$\infty$$-norm of the gradient search vector when the search algorithm terminates. FcnCount Number of times the objective function was called. UpdateNorm Norm of the gradient search vector in the last iteration. Omitted when the search method is 'lsqnonlin' or 'fmincon'. LastImprovement Criterion improvement in the last iteration, expressed as a percentage. Omitted when the search method is 'lsqnonlin' or 'fmincon'. Algorithm Algorithm used by 'lsqnonlin' or 'fmincon' search method. Omitted when other search methods are used. For estimation methods that do not require numerical search optimization, the Termination field is omitted. For more information on using Report, see Estimation Report.

x0

Initial states computed during the estimation. If data contains multiple experiments, then x0 is an array with each column corresponding to an experiment. This value is also stored in the Parameters field of the model’s Report property.

Examples

Determine Optimal Estimated Model Order

Obtain measured input-output data.

load icEngine.mat;
data = iddata(y,u,0.04);

data is an iddata object containing 1500 input-output data samples. The data sample time is 0.04 seconds.

Estimate a state-space model for measured input-output data. Determine the optimal model order within a given model order range.

nx = 1:10;
sys = ssest(data,nx);

A plot that shows the Hankel singular values (SVD) for models of the orders specified by nx appears.

States with relatively small Hankel singular values can be safely discarded. The suggested default order choice is 3.

Select the model order in the Model Order drop-down list and click Apply.

Identify State-Space Model With Input Delay

Load time-domain system response data.

load iddata7 z7;

Identify a fourth-order state-space model of the data. Specify a known delay of 2 seconds for the first input and 0 seconds for the second input.

nx = 4;
sys = ssest(z7(1:300),nx,'InputDelay',[2;0]);

Estimate State-Space Model Using Regularization

Obtain a regularized 5th order state-space model for a 2nd order system from a narrow bandwidth signal.

Load estimation data.

load regularizationExampleData eData;

Create the transfer function model used for generating the estimation data (true system).

trueSys = idtf([0.02008 0.04017 0.02008],[1 -1.561 0.6414],1);

Estimate an unregularized state-space model.

opt = ssestOptions('SearchMethod','lm');
m = ssest(eData,5,'form','modal','DisturbanceModel','none','Ts',eData.Ts,opt);

Estimate a regularized state-space model.

opt.Regularization.Lambda = 10;
mr = ssest(eData,5,'form','modal','DisturbanceModel','none','Ts',eData.Ts,opt);

Compare the model outputs with the estimation data.

compare(eData,m,mr);

Compare the model impulse responses.

impulse(trueSys,m,mr,50);
legend('trueSys','m','mr');

Estimate Partially Known State-Space Model Using Structured Estimation

Estimate a state-space model of measured input-output data. Configure the parameter constraints and initial values for estimation using a state-space model.

Create an idss model to specify the initial parameterization for estimation.

A = blkdiag([-0.1 0.4; -0.4 -0.1],[-1 5; -5 -1]);
B = [1; zeros(3,1)]; 
C = [1 1 1 1]; 
D = 0; 
K = zeros(4,1);
x0 = [0.1 0.1 0.1 0.1];
Ts = 0;
init_sys = idss(A,B,C,D,K,x0,Ts);

Setting all entries of K to 0 creates an idss model with no state disturbance element.

Use the Structure property to fix the values of some of the model parameters. Configure the model so that B and K are fixed, and only the nonzero entries of A are estimable.

init_sys.Structure.A.Free = (A~=0);
init_sys.Structure.B.Free = false;
init_sys.Structure.K.Free = false;

The entries in init_sys.Structure.A.Free determine whether the corresponding entries in init_sys.A are free (true) or fixed (false).

Load the measured data and estimate a state-space model using the parameter constraints and initial values specified by init_sys.

load iddata2 z2;
sys = ssest(z2,init_sys);

The estimated parameters of sys satisfy the constraints specified by init_sys.

Model Order Reduction by Estimation

Consider the Simulink model idF14Model. Linearizing this model gives a ninth-order model. However, the dynamics of the model can be captured, without compromising the fit quality too much, using a lower-order model.

Obtain the linearized model.

load_system('idF14Model');
io = getlinio('idF14Model');
sys_lin = linearize('idF14Model',io);

sys_lin is a ninth-order state-space model with two outputs and one input.

Simulate the step response of the linearized model, and use the data to create an iddata object.

Ts = 0.0444;
t = (0:Ts:4.44)';
y = step(sys_lin,t);

data = iddata([zeros(20,2);y],[zeros(20,1); ones(101,1)],Ts);

data is an iddata object that encapsulates the step response of sys_lin.

Compare the data to the model linearization.

compare(data,sys_lin);

Because the data was obtained by simulating the linearized model, there is a 100% match between the data and model linearization response.

Identify a state-space model with a reduced order that adequately fits the data.

Determine an optimal model order.

nx = 1:9;
sys1 = ssest(data,nx,'DisturbanceModel','none');

A plot showing the Hankel singular values (SVD) for models of the orders specified by nx appears.

States with relatively small Hankel singular values can be safely discarded. The plot suggests using a fifth-order model.

At the MATLAB command prompt, select the model order for the estimated state-space model. Specify the model order as 5, or press Enter to use the default order value.

Compare the data to the estimated model.

compare(data,sys1);

sys1 provides a 98.4% fit for the first output and a 97.7% fit for the second output.

Examine the stopping condition for the search algorithm.

sys1.Report.Termination.WhyStop

ans = 
'Maximum number of iterations reached.'

Create an estimation options set that specifies the 'lm' search method and allows a maximum of 50 search iterations.

opt = ssestOptions('SearchMethod','lm');
opt.SearchOptions.MaxIterations = 50;
opt.Display = 'on';

Identify a state-space model using the estimation option set and sys1 as the estimation initialization model.

sys2 = ssest(data,sys1,opt);

Compare the response of the linearized and the estimated models.

compare(data,sys_lin,sys2);

sys2 provides a 99% fit for the first output and a 98% fit for the second output while using 4 less states than sys_lin .

More About

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 $$({\lambda }_1,\sigma \pm j\omega ,{\lambda }_2)$$, 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(s)=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& \mathrm{..}& \mathrm{..}& 0& -{\alpha }_n\\ 1& 0& 0& \mathrm{..}& 0& -{\alpha }_n-1\\ 0& 1& 0& .& :& :\\ :& 0& .& .& :& :\\ 0& .& .& 1& 0& -{\alpha }_2\\ 0& \mathrm{..}& \mathrm{..}& 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.

Algorithms

ssest initializes the parameter estimates using either a noniterative subspace approach or an iterative rational function estimation approach. It then refines the parameter values using the prediction error minimization approach. See pem and ssestOptions for more information.

References

[1] Ljung, L. System Identification: Theory For the User, Second Edition, Upper Saddle River, N.J: Prentice Hall, 1999.

Extended Capabilities

Automatic Parallel Support

Accelerate code by automatically running computation in parallel using Parallel Computing Toolbox™.

Parallel computing support is available for estimation using the lsqnonlin search method (requires Optimization Toolbox™). To enable parallel computing, use ssestOptions, set SearchMethod to 'lsqnonlin', and set SearchOptions.Advanced.UseParallel to true.

For example:

opt = ssestOptions;
opt.SearchMethod = 'lsqnonlin';
opt.SearchOptions.Advanced.UseParallel = true;