For multi-experiment data, the sample times and intersample
behavior of all the experiments must match.

np

Number of poles in the estimated transfer function.

np is a nonnegative number.

For systems that are multiple-input, or multiple-output, or
both:

To use the same number of poles for all the input/output
pairs, specify np as a scalar.

To use different number of poles for the input/output
pairs, specify np as an ny-by-nu matrix. ny is
the number of outputs, and nu is the number of
inputs.

nz

Number of zeros in the estimated transfer function.

nz is a nonnegative number.

For systems that are multiple-input, or multiple-output, or
both:

To use the same number of zeros for all the input/output
pairs, specify nz as a scalar.

To use a different number of zeros for the input/output
pairs, specify nz as an ny-by-nu matrix. ny is
the number of outputs, and nu is the number of
inputs.

For a continuous-time model, estimated using discrete-time data,
set nz <= np.

iodelay

Transport delay.

For continuous-time systems, specify transport delays in the
time unit stored in the TimeUnit property of data.
For discrete-time systems, 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 an ny-by-nu array.
Each entry of this array is a numerical value that represents 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.

The specified values are treated as fixed delays.

iodelay must contain either nonnegative
numbers or NaNs. Use NaN in the iodelay matrix
to denote unknown transport delays.

Use [] or 0 to indicate
that there is no transport delay.

opt

Estimation options.

opt is an options set, created using tfestOptions, that specifies estimation
options including:

You obtain init_sys by either performing
an estimation using measured data or by direct construction.

If init_sys is an idtf model, tfest 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 the numerator, denominator,
and transport lag. For example:

To specify an initial guess for the numerator of init_sys,
set init_sys.Structure.Numerator.Value to the initial
guess.

To specify constraints for the numerator of init_sys:

Set init_sys.Structure.Numerator.Minimum to
the minimum numerator coefficient values

Set init_sys.Structure.Numerator.Maximum to
the maximum numerator coefficient values

Set init_sys.Structure.Numerator.Free to
indicate which numerator coefficients are free for estimation

If init_sys is not an idtf model,
the software first converts init_sys to a transfer
function. tfest 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 single quotes (' ').
You can specify several name and value pair
arguments in any order as Name1,Value1,...,NameN,ValueN.

'Ts'

Sample time.

Use the following values for Ts:

0 — Continuous-time model.

data.Ts — Discrete-time
model. In this case, np and nz refer
to the number of roots of z^-1 for the numerator
and denominator polynomials.

Default: 0

'InputDelay'

Input delay for each input channel, specified as a scalar value
or 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 sample times.

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.

Default: 0

'Feedthrough'

Feedthrough for discrete-time transfer function. Must be a Ny-by-Nu logical
matrix. Use a scalar to specify a common value across all channels.

A discrete-time model with 2 poles and 3 zeros takes the following
form:

Hz−1=b0+b1z−1+b2z−2+b3z−31+a1z−1+a2z−2

When the model has direct feedthrough, b0 is
a free parameter whose value is estimated along with the rest of the
model parameters b1, b2, b3, a1, a2.
When the model has no feedthrough, b0 is fixed
to zero.

Default: false (Ny,Nu)

Output Arguments

sys

Identified transfer function, returned as an idtf 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.

InitMethod

Algorithm used to initialize the numerator and denominator
for estimation of continuous-time transfer functions using time-domain
data, returned as one of the following values:

'iv' — Instrument Variable
approach.

'svf' — State Variable Filters
approach.

'gpmf' — Generalized Poisson
Moment Functions approach.

'n4sid' — Subspace state-space
estimation approach.

This field is especially useful to view the algorithm
used when the InitMethod option in the estimation
option set is 'all'.

N4Weight

Weighting matrices used in the singular-value decomposition
step when InitMethod is 'n4sid',
returned as one of the following values:

'MOESP' — Uses the MOESP
algorithm by Verhaegen.

'CVA' — Uses the canonical
variable algorithm (CVA) by Larimore.

'SSARX' — A subspace identification
method that uses an ARX estimation based algorithm to compute the
weighting.

This field is especially useful to view the weighting
matrices used when the N4Weight option in the estimation
option set is 'auto'.

N4Horizon

Forward and backward prediction horizons used when InitMethod is 'n4sid',
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.

InitialCondition

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

Analyze the Origin of Delay in Measured Data

Compare two discrete-time models with and without feedthrough and transport delay.

If there is a delay from the measured input to output, it can be attributed to a lack of feedthrough or to a true transport delay. For discrete-time models, absence of feedthrough corresponds to a lag of 1 sample between the input and output. Estimating a model with Feedthrough = false and IODelay = 0 thus produces a discrete-time system that is equivalent to a system with Feedthrough = true and IODelay = 1. Both systems show the same time- and frequency-domain responses, for example, on step and Bode plots. However, you get different results if you reduce these models using balred or convert them to their continuous-time representation. Therefore, you should check if the observed delay should be attributed to transport delay or to a lack of feedthrough.

Estimate a discrete-time model with no feedthrough.

load iddata1z1
np = 2;
nz = 2;
model1 = tfest(z1,np,nz,'Ts',z1.Ts);

model1 has a transport delay of 1 sample and its IODelay property is 0. Its numerator polynomial begins with .

Choice of Feedthrough dictates whether the leading numerator coefficient is zero (no feedthrough) or not (nonzero feedthrough). Delays are expressed separately using InputDelay or IODelay property. This example uses InputDelay to express the delays.

Validate the estimated model. Exclude the data outliers for validation.

Estimate Transfer Function with Unknown, Constrained Transport Delays

Create a transfer function model with the expected numerator and denominator structure and delay constraints.

In this example, the experiment data consists of two inputs and one output. Both transport delays are unknown and have an identical upper bound. Additionally, the transfer functions from both inputs to the output are identical in structure.

init_sys is an idtf model describing the structure of the transfer function from one input to the output. The transfer function consists of one zero, three poles and a transport delay. The use of NaN indicates unknown coefficients.

init_sys.Structure(1).IODelay.Free = true indicates that the transport delay is not fixed.

np specifies the number of poles in the estimated transfer function. The first element of np indicates that the transfer function from the first input to the output contains 3 poles. Similarly, the second element of np indicates that the transfer function from the second input to the output contains 4 poles.

nz specifies the number of zeros in the estimated transfer function. The first element of nz indicates that the transfer function from the first input to the output contains 1 zero. Similarly, the second element of np indicates that the transfer function from the second input to the output does not contain any zeros.

iodelay specifies the transport delay from the first input to the output as 4 seconds. The transport delay from the second input to the output is specified as 0.6 seconds.

data is an idfrd object containing the continuous-time frequency response for G.

Estimate a transfer function for data.

np = [3 4];
nz = [1 0];
sys = tfest(data,np,nz);

np specifies the number of poles in the estimated transfer function. The first element of np indicates that the transfer function from the first input to the output contains 3 poles. Similarly, the second element of np indicates that the transfer function from the second input to the output contains 4 poles.

nz specifies the number of zeros in the estimated transfer function. The first element of nz indicates that the transfer function from the first input to the output contains 1 zero. Similarly, the second element of nz indicates that the transfer function from the second input to the output does not contain any zeros.

More About

Compatibility

Starting in R2016b, a new algorithm is used for performing transfer
function estimation from frequency-domain data. You are likely to
see faster and more accurate results with the new algorithm, particularly
for data with dynamics over a large range of frequencies and amplitudes.
However, the estimation results may not match results from previous
releases. To perform estimation using the previous estimation algorithm,
append '-R2016a' to the syntax.

For example, suppose that you are estimating a transfer function
model with np poles using frequency-domain data, data.

sys = tfest(data,np)

To use the previous estimation algorithm, use the following
syntax.

sys = tfest(data,np,'-R2016a')

Algorithms

The details of the estimation algorithms used by tfest vary
depending on a variety of factors, including the sampling of the estimated
model and the estimation data.

Continuous-Time Transfer Function Estimation Using Time-Domain Data

Parameter Initialization

The estimation algorithm initializes the estimable parameters
using the method specified by the InitMethod estimation
option. The default method is the Instrument Variable (IV) method.

The State-Variable Filters (SVF) approach and the Generalized
Poisson Moment Functions (GPMF) approach to continuous-time parameter
estimation use prefiltered data [1][2].
The constant 1λ in [1] and [2] corresponds
to the initialization option (InitOption) field FilterTimeConstant.
IV is the simplified refined IV method and is called SRIVC in [3].
This method has a prefilter that is the denominator of the current
model, initialized with SVF. This prefilter is iterated up to MaxIter times,
until the model change is less than Tolerance. MaxIter and Tolerance are
options that you can specify using the InitOption structure.
The 'n4sid' initialization option estimates a discrete-time
model, using the N4SID estimation algorithm, that it transforms to
continuous-time using d2c.

You use tfestOptions to
create the option set used to estimate a transfer function.

Parameter Update

The initialized parameters are updated using a nonlinear least-squares
search method, specified by the SearchMethod estimation
option. The objective of the search method is to minimize the weighted
prediction error norm.

Discrete-Time Transfer Function Estimation Using Time-Domain Data

For discrete-time data, tfest uses the
same algorithm as oe to determine
the numerator and denominator polynomial coefficients. In this algorithm,
the initialization is performed using arx,
followed by nonlinear least-squares search based updates to minimize
a weighted prediction error norm.

Perform a bilinear mapping to transform the domain
(frequency grid) of the transfer function. For continuous-time models,
the imaginary axis is transformed to the unit disk. For discrete-time
models, the original domain unit disk is transformed to another unit
disk.

Here W is a frequency-dependent weight that
you specify. D is the denominator of the transfer
function model that is to be estimated, and Ni is
the numerator corresponding to the ith input. y and u are
the measured output and input data, respectively. nf and nu are
the number of frequencies and inputs, and w is
the frequency. Rearranging the terms gives:

where m is the current iteration, and Dm-1(ω) is
the denominator response identified at the previous iteration. Now
each step of the iteration is a linear least-squares problem, where
the identified parameters capture the responses Dm(ω) and Ni,m(ω) for i =
1,2,...nu. The iteration
is initialized by choosing D0(ω) =
1.

The first iteration of the algorithm identifies D1(ω).
The D1(ω) and Ni,1(ω) polynomials
are expressed in monomial basis.

The second and following iterations express the polynomials Dm(ω) and Ni,m(ω) in
terms of orthogonal rational basis functions on the unit disk. These
basis functions have the form:

Bj,m(ω)=(1−|λj,m−1|2q−λj,m−1)∏r=0j−11−(λj,m−1)*q(ω)q(ω)−λr,m−1

Here λj,m-1 is
the jth pole that is identified at the previous
step, m-1, of the iteration. λj,m-1* is
the complex conjugate of λj,m-1,
and q is the frequency-domain variable on the unit
disk.

The algorithm runs for a maximum of 20 iterations.
The iterations are terminated early if the relative change in the
value of the loss function is less than 0.001 in the last three iterations.

If you specify bounds on transfer function coefficients, these
bounds correspond to affine constraints on the identified parameters.
If you only have equality constraints (fixed transfer function coefficients),
the corresponding equality constrained least-squares problem is solved
algebraically. To do so, the software computes an orthogonal basis
for the null space of the equality constraint matrix, and then solves
least-squares problem within this null space. If you have upper or
lower bounds on transfer function coefficients, the corresponding
inequality constrained least-squares problem is solved using interior-point
methods.

Perform linear refinements — The S-K iterations,
even when they converge, do not always yield a locally optimal solution.
To find a critical point of the optimization problem that may yield
a locally optimal solution, a second set of iterations are performed.
The critical points are solutions to a set of nonlinear equations.
The algorithm searches for a critical point by successively constructing
a linear approximation to the nonlinear equations and solving the
resulting linear equations in the least-squares sense. The equations
are:

The first iteration is started with the best solution found
for the numerators Ni and
denominator D parameters during S-K iterations.
Unlike S-K iterations, the basis functions Bj(ω) are
not changed at each iteration, the iterations are performed with the
basis functions that yielded the best solution in the S-K iterations.
As before, the algorithm runs for a maximum of 20 iterations. The
iterations are terminated early if the relative change in the value
of the loss function is less than 0.001 in the last three iterations.

If you specify bounds on transfer function coefficients, these
bounds are incorporated into the necessary optimality
conditions via generalized Lagrange multipliers. The resulting constrained
linear least-squares problems are solved using the same methods explained
in the S-K iterations step.

Return the transfer function parameters corresponding
to the optimal solution — Both the S-K and linear refinement
iteration steps do not guarantee an improvement in the loss function
value. The algorithm tracks the best parameter value observed during
these steps, and returns these values.

Invert the bilinear mapping performed in step 1.

Perform an iterative refinement of the transfer function
parameters using the nonlinear least-squares search method specified
in the SearchMethod estimation option. This step
is implemented in the following situations:

When you specify the EnforceStability estimation
option as true (stability is requested), and the
result of step 5 of this algorithm is an unstable model. The unstable
poles are reflected inside the stability boundary and the resulting
parameters are iteratively refined. For information about estimation
options, see tfestOptions.

You use frequency domain input-output data to identify
a multi-input model.

If you are using the estimation algorithm from R2016a or earlier
(see Compatibility) for estimating a continuous-time model
using continuous-time frequency-domain data, then for continuous-time
data and fixed delays, the Output-Error algorithm is used for model
estimation. For continuous-time data and free delays, the state-space
estimation algorithm is used. In this algorithm, the model coefficients
are initialized using the N4SID estimation method. This initialization
is followed by nonlinear least-squares search based updates to minimize
a weighted prediction error norm.

tfest command first estimates a discrete-time
model from the discrete-time data. The estimated model is then converted
to a continuous-time model using the d2c command.
The frequency response of the resulting continuous-time model is then
computed over the frequency grid of the estimation data. A continuous-time
model of the desired (user-specified) structure is then fit to this
frequency response. The estimation algorithm for using the frequency-response
data to obtain the continuous-time model is the same as that for continuous-time
transfer function estimation using continuous-time data.

If you are using the estimation algorithm from R2016a or earlier
(see Compatibility), the state-space estimation algorithm
is used for estimating continuous-time models from discrete-time data.
In this algorithm, the model coefficients are initialized using the
N4SID estimation method. This initialization is followed by nonlinear
least-squares search based updates to minimize a weighted prediction
error norm.

Delay Estimation

When delay values are specified as NaN,
they are estimated separate from the model numerator and denominator
coefficients, using delayest.
The delay values thus determined are treated as fixed values during
the iterative update of the model using a nonlinear least-squares
search method. Thus, the delay values are not iteratively updated.