Estimate output-error polynomial model using time-domain or frequency-domain data
Output-error (OE) models are a special configuration of polynomial models, having
only two active polynomials—B and F. OE models represent
conventional transfer functions that relate measured inputs to outputs while also including
white noise as an additive output disturbance. You can estimate OE models using time- and
frequency-domain data. The tfest
command offers the same functionality as
oe
. For tfest
, you specify the model orders using
number of poles and zeros rather than polynomial degrees. For continuous-time estimation,
tfest
provides faster and more accurate results, and is
recommended.
estimates an OE model sys
= oe(data
,[nb
nf nk]
)sys
, represented by
y(t) is the output, u(t) is the input, and e(t) is the error.
oe
estimates sys
using the measured
input-output data data
, which can be in the time or the frequency
domain. The orders [nb nf nk]
define the number of parameters in each
component of the estimated polynomial.
specifies model structure attributes using additional options specified by one or more
name-value pair arguments.sys
= oe(data
,[nb
nf nk]
,Name,Value
)
[
returns the estimated initial conditions as an sys
,ic
] = oe(___)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.
Estimate an OE polynomial from time-domain data using two methods to specify input delay.
Load the estimation data.
load iddata1 z1
Set the orders of the B and F polynomials nb
and nf
. Set the input delay nk
to one sample. Compute the model sys
.
nb = 2; nf = 2; nk = 1; sys = oe(z1,[nb nf nk]);
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%.
Instead of using nk
, you can also use the name-value pair argument 'InputDelay'
to specify the one-sample delay.
nk = 0;
sys1 = oe(z1,[nb nf nk],'InputDelay',1);
figure
compare(z1,sys1)
The results are identical.
You can view more information about the estimation by exploring the idpoly
property sys.Report
.
sys.Report
ans = Status: 'Estimated using OE' Method: 'OE' InitialCondition: 'zero' Fit: [1x1 struct] Parameters: [1x1 struct] OptionsUsed: [1x1 idoptions.polyest] RandState: [1x1 struct] DataUsed: [1x1 struct] Termination: [1x1 struct]
For example, find out more information about the termination conditions.
sys.Report.Termination
ans = struct with fields:
WhyStop: 'Near (local) minimum, (norm(g) < tol).'
Iterations: 3
FirstOrderOptimality: 0.0708
FcnCount: 7
UpdateNorm: 1.4809e-05
LastImprovement: 5.1744e-06
The report includes information on the number of iterations and the reason the estimation stopped iterating.
Load the estimation data.
load oe_data1 data;
The idfrd
object data
contains the continuous-time frequency response for the following model:
Estimate the model.
nb = 2; nf = 3; sys = oe(data,[nb nf]);
Evaluate the goodness of fit.
compare(data,sys);
Estimate a high-order OE model from data collected by simulating a high-order system. Determine the regularization constants by trial and error and use the values for model estimation.
Load the data.
load regularizationExampleData.mat m0simdata
Estimate an unregularized OE model of order 30.
m1 = oe(m0simdata,[30 30 1]);
Obtain a regularized OE model by determining the Lambda value using trial and error.
opt = oeOptions; opt.Regularization.Lambda = 1; m2 = oe(m0simdata,[30 30 1],opt);
Compare the model outputs with the estimation data.
opt = compareOptions('InitialCondition','z'); compare(m0simdata,m1,m2,opt);
The regularized model m2
produces a better fit than the unregularized model m1
.
Compare the variance in the model responses.
h = bodeplot(m1,m2); opt = getoptions(h); opt.PhaseMatching = 'on'; opt.ConfidenceRegionNumberSD = 3; opt.PhaseMatching = 'on'; setoptions(h,opt); showConfidence(h);
The regularized model m2
has a reduced variance compared to the unregularized model m1
.
Load the estimation data data
and sample time Ts
.
load oe_data2.mat data Ts
An iddata
object data
contains the discrete-time frequency response for the following model:
View the estimation sample time Ts
that you loaded.
Ts
Ts = 1.0000e-03
This value matches the property data.Ts
.
data.Ts
ans = 1.0000e-03
You can estimate a continuous model from data
by limiting the input and output frequency bands to the Nyquist frequency. To do so, specify the estimation prefilter option 'WeightingFilter
' to define a passband from 0
to 0.5*pi/Ts
rad/s. The software ignores any response values with frequencies outside of that passband.
opt = oeOptions('WeightingFilter',[0 0.5*pi/Ts]);
Set the Ts
property to 0
to treat data
as continuous-time data.
data.Ts = 0;
Estimate the continuous model.
nb = 1; nf = 3; sys = oe(data,[nb nf],opt);
Load the data.
load iddata1ic z1i
Estimate an OE polynomial model sys
and return the initial conditions in ic
.
nb = 2; nf = 2; nk = 1; [sys,ic] = oe(z1i,[nb,nf,nk]); ic
ic = initialCondition with properties: A: [2x2 double] X0: [2x1 double] C: [0.9428 0.4824] 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.
data
— Estimation dataiddata
object | frd
object | idfrd
objectEstimation data, specified as an iddata
object, an frd
object, or an idfrd
object.
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:
Time-domain 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.
You can compute discrete-time models from time-domain data or discrete-time
frequency-domain data. Use tfest
to compute continuous-time
models.
[nb nf nk]
— OE model ordersOE model orders, specified as a 1-by-3 vector or a vector of integer matrices.
For a system represented by
where y(t) is the output,
u(t) is the input, and
e(t) is the error, the elements of [nb
nf nk]
are as follows:
nb
— Order of the B(q)
polynomial + 1, which is equivalent to the length of the
B(q) polynomial. nb
is an
Ny-by-Nu
matrix. Ny is the number of outputs and
Nu is the number of inputs.
nf
— Order of the F polynomial.
nf
is an
Ny-by-Nu
matrix.
nk
— Input delay, expressed as the number of samples.
nk
is an
Ny-by-Nu
matrix. The delay appears as leading zeros of the B
polynomial.
For estimation using continuous-time frequency-domain data, specify only
[nb nf]
and omit nk
. For an example, see Estimate Continuous-Time OE Model Using Frequency Response.
init_sys
— Linear system idpoly
model | linear model | structureLinear system that configures the initial parameterization of
sys
, specified as an idpoly
model, another
linear model, or a structure. You obtain init_sys
either by
performing an estimation using measured data or by direct construction.
If init_sys
is an idpoly
model of the OE
structure, oe
uses the parameter values of
init_sys
as the initial guess for estimating
sys
. The sample time of init_sys
must match
the sample time of the data.
Use the Structure
property of init_sys
to
configure initial guesses and constraints for B(q)
and F(q). For example:
To specify an initial guess for the F(q)
term of init_sys
, set
init_sys.Structure.F.Value
as the initial guess.
To specify constraints for the B(q) term
of init_sys
:
Set init_sys.Structure.B.Minimum
to the minimum
B(q) coefficient values.
Set init_sys.Structure.B.Maximum
to the maximum
B(q) coefficient values.
Set init_sys.Structure.B.Free
to indicate which
B(q) coefficients are free for
estimation.
If init_sys
is not a polynomial model of the OE structure, the
software first converts init_sys
to an OE structure model.
oe
uses the parameters of the resulting model as the initial
guess for estimating
sys
.
If you do not specify opt
and init_sys
was
obtained by estimation, then the software uses estimation options from
init_sys.Report.OptionsUsed
.
opt
— Estimation optionsoeOptions
option setEstimation options, specified as an oeOptions
option set. Options specified by opt
include:
Estimation objective
Handling of initial conditions
Numerical search method and the associated options
For examples of specifying estimation options, see Estimate Continuous Model Using Band-Limited Discrete-Time Frequency-Domain Data.
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',1
'InputDelay'
— Input delaysInput delays for each input channel, specified as the comma-separated pair
consisting of 'InputDelay'
and a numeric vector.
For continuous-time models, specify 'InputDelay'
in the
time units stored in the TimeUnit
property.
For discrete-time models, specify 'InputDelay'
in integer
multiples of the sample time Ts
. For example, setting
'InputDelay'
to 3
specifies 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.
To apply the same delay to all channels, specify 'InputDelay'
as a scalar.
For an example, see Estimate OE Polynomial Model.
'IODelay'
— Transport delaysTransport delays for each input-output pair, specified as the comma-separated pair
consisting of 'IODelay'
and a numeric array.
For continuous-time models, specify 'IODelay'
in the
time units stored in the TimeUnit
property.
For discrete-time models, specify 'IODelay'
in integer
multiples of the sample time Ts
. For example, setting
'IODelay'
to 4
specifies a transport
delay of four sampling periods.
For a system with Nu inputs and
Ny outputs, set
'IODelay'
to an
Ny-by-Nu
matrix. Each entry is an integer value representing the transport delay for the
corresponding input-output pair.
To apply the same delay to all channels, specify 'IODelay'
as
a scalar.
You can specify 'IODelay'
as an alternative to the
nk
value. Doing so simplifies the model structure by reducing the
number of leading zeros in the B polynomial. In particular, you can
represent max(nk-1,0)
leading zeros as input-output delays using
'IODelay'
instead.
sys
— OE polynomial modelidpoly
objectOE polynomial model that fits the estimation data, returned as an idpoly
model object. This model is created using the specified model
orders, delays, and estimation options. The sample time of sys
matches the sample time of the estimation data. Therefore, sys
is
always a discrete-time model when estimated from time-domain data. For continuous-time
model identification using time-domain data, use tfest
.
The Report
property of the model stores information about the
estimation results and options used. 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. | ||||||||||||||||||
InitialCondition | Handling of initial conditions during model estimation, returned as one of the following values:
This field is especially useful to view
how the initial conditions were handled when the | ||||||||||||||||||
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 no custom options were
configured, this is a set of default options. See | ||||||||||||||||||
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:
| ||||||||||||||||||
Termination | Termination conditions for the iterative search used for prediction error minimization, returned as a structure with the following fields:
For estimation methods that do not require numerical search optimization,
the |
For more information on using Report
, see Estimation Report.
ic
— Initial conditionsinitialCondition
object | object array of initialCondition
valuesEstimated 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 oe
returns ic
values of
0
and the you know that you have non-zero initial conditions, set
the 'InitialCondition'
option in oeOptions
to 'estimate'
and pass the updated option set
to oe
. For
example:
opt = oeOptions('InitialCondition,'estimate')
[sys,ic] = oe(data,np,nz,opt)
'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.
The general output-error model structure is:
The orders of the output-error model are:
If data
is continuous-time frequency-domain data,
oe
estimates a continuous-time model with the following transfer
function:
The orders of the numerator and denominator are nb
and
nf
, similar to the discrete-time case. However, the sample delay
nk
does not exist in the continuous case, and you should not specify
nk
when you command the estimation. Instead, express any system delay
using the name-value pair argument 'IODelay'
along with the system delay
in the time units that are stored in the property TimeUnit
. For example,
suppose that your continuous system has a delay of iod
seconds. Use
model = oe(data,[nb nf],'IODelay',iod)
.
Parallel computing support is available for estimation using the
lsqnonlin
search method (requires Optimization Toolbox™). To enable parallel computing, use oeOptions
, set SearchMethod
to
'lsqnonlin'
, and set
SearchOptions.Advanced.UseParallel
to true
.
For example:
opt = oeOptions;
opt.SearchMethod = 'lsqnonlin';
opt.SearchOptions.Advanced.UseParallel = true;
armax
| arx
| bj
| compare
| iddata
| idfrd
| idpoly
| iv4
| n4sid
| oeOptions
| polyest
| sim
| tfest
You have a modified version of this example. Do you want to open this example with your edits?