bj

Estimate Box-Jenkins polynomial model using time domain data

Syntax

sys = bj(data, [nb nc nd nf nk])
sys = bj(data,[nb nc nd nf nk], Name,Value)
sys = bj(data, init_sys)
sys = bj(data, ___, opt)
[sys,ic] = bj(___)

Description

sys = bj(data, [nb nc nd nf nk]) estimates a Box-Jenkins polynomial model, sys, using the time-domain data, data. [nb nc nd nf nk] define the orders of the polynomials used for estimation.

sys = bj(data,[nb nc nd nf nk], Name,Value) estimates a polynomial model with additional options specified by one or more Name,Value pair arguments.

sys = bj(data, init_sys) estimates a Box-Jenkins polynomial using the polynomial model init_sys to configure the initial parameterization of sys.

sys = bj(data, ___, opt) estimates a Box-Jenkins polynomial using the option set, opt, to specify estimation behavior.

[sys,ic] = bj(___) returns the estimated initial conditions as an initialCondition object. Use this syntax if you plan to simulate or predict the model response using the same estimation input data and then compare the response with the same estimation output data. Incorporating the initial conditions yields a better match during the first part of the simulation.

Input Arguments

 data Estimation data. data is an iddata object that contains time-domain input and output signal values. You cannot use frequency-domain data for estimating Box-Jenkins models. [nb nc nd nf nk] A vector of matrices containing the orders and delays of the Box-Jenkins model. Matrices must contain nonnegative integers. nb is the order of the B polynomial plus 1 (Ny-by-Nu matrix)nc is the order of the C polynomial plus 1 (Ny-by–1 matrix)nd is the order of the D polynomial plus 1 (Ny-by-1 matrix)nf is the order of the F polynomial plus 1 (Ny-by-Nu matrix)nk is the input delay (in number of samples, Ny-by-Nu matrix) where Nu is the number of inputs and Ny is the number of outputs. opt Estimation options. opt is an options set that configures, among others, the following: estimation objectiveinitial conditionsnumerical search method to be used in estimation Use bjOptions to create the options set. init_sys Polynomial model that configures the initial parameterization of sys. init_sys must be an idpoly model with the Box-Jenkins structure that has only B, C, D and F polynomials active. bj uses the parameters and constraints defined in init_sys as the initial guess for estimating sys. Use the Structure property of init_sys to configure initial guesses and constraints for B(q), F(q), C(q) and D(q). To specify an initial guess for, say, the C(q) term of init_sys, set init_sys.Structure.C.Value as the initial guess. To specify constraints for, say, the B(q) term of init_sys: set init_sys.Structure.B.Minimum to the minimum B(q) coefficient valuesset init_sys.Structure.B.Maximum to the maximum B(q) coefficient valuesset init_sys.Structure.B.Free to indicate which B(q) coefficients are free for estimation You can similarly specify the initial guess and constraints for the other polynomials.

Name-Value 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.

 InputDelay Input delays. InputDelay is a numeric vector specifying a time delay for each input channel. 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, where each entry is a numerical value representing 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. Default: 0 for all input channels IODelay Transport delays. IODelay is a numeric array specifying a separate transport delay for each input/output pair. Specify transport delays as integers denoting delay of a multiple of the sample time Ts. For a MIMO system with Ny outputs and Nu inputs, set IODelay to a Ny-by-Nu array, where each entry is a numerical value representing the transport delay for the corresponding input/output pair. You can also set IODelay to a scalar value to apply the same delay to all input/output pairs. Default: 0 for all input/output pairs IntegrateNoise Logical specifying integrators in the noise channel. IntegrateNoise is a logical vector of length Ny, where Ny is the number of outputs. Setting IntegrateNoise to true for a particular output results in the model: $y\left(t\right)=\frac{B\left(q\right)}{F\left(q\right)}u\left(t-nk\right)+\frac{C\left(q\right)}{D\left(q\right)}\frac{e\left(t\right)}{1-{q}^{-1}}$ Where, $\frac{1}{1-{q}^{-1}}$ is the integrator in the noise channel,e(t). Default: false(Ny,1)(Ny is the number of outputs)

Output Arguments

sys

BJ model that fits the estimation data, returned as a discrete-time idpoly object. 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.

InitialCondition

Handling of initial conditions during model estimation, returned as one of the following values:

• 'zero' — The initial conditions were set to zero.

• 'estimate' — The initial conditions were treated as independent estimation parameters.

• 'backcast' — The initial conditions were estimated using the best least squares fit.

This field is especially useful to view how the initial conditions were handled when the InitialCondition option in the estimation option set is 'auto'.

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 the percentage fitpercent = 100(1-NRMSE).

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 bjOptions 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.

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 [].

Termination

Termination conditions for the iterative search used for prediction error minimization, returned as a structure with the following fields:

FieldDescription
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.

[nb nc nd nf nk]

A vector of matrices containing the orders and delays of the Box-Jenkins model. Matrices must contain nonnegative integers.

• nb is the order of the B polynomial plus 1 (Ny-by-Nu matrix)

• nc is the order of the C polynomial plus 1 (Ny-by–1 matrix)

• nd is the order of the D polynomial plus 1 (Ny-by-1 matrix)

• nf is the order of the F polynomial plus 1 (Ny-by-Nu matrix)

• nk is the input delay (in number of samples, Ny-by-Nu matrix) where Nu is the number of inputs and Ny is the number of outputs.

ic

Estimated initial conditions, returned as an initialCondition object or an object array of initialCondition values.

• For a single-experiment data set, ic represents, in state-space form, the free response of the transfer function model (A and C matrices) to the estimated initial states (x0).

• For a multiple-experiment data set with Ne experiments, ic is an object array of length Ne that contains one set of initialCondition values for each experiment.

If bj returns ic values of 0 and the you know that you have non-zero initial conditions, set the 'InitialCondition' option in bjOptions to 'estimate' and pass the updated option set to bj. For example:

opt = bjOptions('InitialCondition','estimate')
[sys,ic] = bj(data,[nb nc nd nf nk],opt)
The default 'auto' setting of 'InitialCondition' uses the 'zero' method when the initial conditions have a negligible effect on the overall estimation-error minimization process. Specifying 'estimate' ensures that the software estimates values for ic.

For more information, see initialCondition. For an example of using this argument, see Obtain Initial Conditions.

Examples

collapse all

Estimate the parameters of a single-input, single-output Box-Jenkins model from measured data.

nb = 2;
nc = 2;
nd = 2;
nf = 2;
nk = 1;
sys = bj(z1,[nb nc nd nf nk]);

sys is a discrete-time idpoly model with estimated coefficients. The order of sys is as described by nb, nc, nd, nf, and nk.

Use getpvec to obtain the estimated parameters and getcov to obtain the covariance associated with the estimated parameters.

Estimate the parameters of a multi-input, single-output Box-Jenkins model from measured data.

nb = [2 1 1];
nc = 1;
nd = 1;
nf = [2 1 2];
nk = [5 10 15];
sys = bj(z8,[nb nc nd nf nk]);

sys estimates the parameters of a model with three inputs and one output. Each of the inputs has a delay associated with it.

Estimate a regularized BJ model by converting a regularized ARX model.

Estimate an unregularized BJ model of order 30.

m1 = bj(m0simdata(1:150),[15 15 15 15 1]);

Estimate a regularized BJ model by determining Lambda value by trial and error.

opt = bjOptions;
opt.Regularization.Lambda = 1;
m2 = bj(m0simdata(1:150),[15 15 15 15 1],opt);

Obtain a lower-order BJ model by converting a regularized ARX model followed by order reduction.

opt1 = arxOptions;
[L,R] = arxRegul(m0simdata(1:150),[30 30 1]);
opt1.Regularization.Lambda = L;
opt1.Regularization.R = R;
m0 = arx(m0simdata(1:150),[30 30 1],opt1);
mr = idpoly(balred(idss(m0),7));

Compare the model outputs against data.

opt2 = compareOptions('InitialCondition','z');
compare(m0simdata(150:end),m1,m2,mr,opt2); Estimate the parameters of a single-input, single-output Box-Jenkins model while configuring some estimation options.

Generate estimation data.

B = [0 1 0.5];
C = [1 -1 0.2];
D = [1 1.5 0.7];
F = [1 -1.5 0.7];
sys0 = idpoly(1,B,C,D,F,0.1);
e = iddata([],randn(200,1));
u = iddata([],idinput(200));
y = sim(sys0,[u e]);
data = [y u];

data is a single-input, single-output data set created by simulating a known model.

Estimate initial Box-Jenkins model.

nb = 2;
nc = 2;
nd = 2;
nf = 2;
nk = 1;
init_sys = bj(data,[2 2 2 2 1]);

Create an estimation option set to refine the parameters of the estimated model.

opt = bjOptions;
opt.Display = 'on';
opt.SearchOptions.MaxIterations = 50;

opt is an estimation option set that configures the estimation to iterate 50 times at most and display the estimation progress.

Reestimate the model parameters using the estimation option set.

sys = bj(data,init_sys,opt);

sys is estimated using init_sys for the initial parameterization for the polynomial coefficients.

To view the estimation result, enter sys.Report.

Estimate a multi-input, multi-output Box-Jenkins model from estimated data.

data = [z1 z2(1:300)];

data contains the measured data for two inputs and two outputs.

Estimate the model.

nb = [2 2; 3 4];
nc = [2;2];
nd = [2;2];
nf = [1 0; 2 2];
nk = [1 1; 0 0];
sys = bj(data,[nb nc nd nf nk]);

The polynomial order coefficients contain one row for each output.

sys is a discrete-time idpoly model with two inputs and two outputs.

Estimate a second-order Box-Jenkins model sys and return the initial conditions in ic.

nb = 2;
nc = 2;
nd = 2;
nf = 2;
nk = 1;
[sys,ic] = bj(z1i,[nb nc nd nf nk]);
ic
ic =
initialCondition with properties:

A: [4x4 double]
X0: [4x1 double]
C: [0.8744 0.5426 0.4647 -0.5285]
Ts: 0.1000

ic is an initialCondition object that encapsulates the free response of sys, in state-space form, to the initial state vector in X0. You can incorporate ic when you simulate sys with the z1i input signal and compare the response with the z1i output signal.

collapse all

Box-Jenkins Model Structure

The general Box-Jenkins model structure is:

$y\left(t\right)=\sum _{i=1}^{nu}\frac{{B}_{i}\left(q\right)}{{F}_{i}\left(q\right)}{u}_{i}\left(t-n{k}_{i}\right)+\frac{C\left(q\right)}{D\left(q\right)}e\left(t\right)$

where nu is the number of input channels.

The orders of Box-Jenkins model are defined as follows:

Alternatives

To estimate a continuous-time model, use:

• tfest — returns a transfer function model

• ssest — returns a state-space model

• bj to first estimate a discrete-time model and convert it a continuous-time model using d2c.

References

 Ljung, L. System Identification: Theory for the User, Upper Saddle River, NJ, Prentice-Hall PTR, 1999.

Extended Capabilities 