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 $$\frac{1}{\lambda}$$ 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.
Continuous-Time Transfer Function Estimation Using Continuous-Time Frequency-Domain Data
The estimation algorithm performs the following tasks:
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.
Perform S-K iterations [4] to
solve a nonlinear least-squares problem — Consider a multi-input
single-output system. The nonlinear least-squares problem is to minimize
the loss function:
$$\underset{D,{N}_{i}}{\text{minimize}}{\displaystyle \sum _{k=1}^{{n}_{f}}{\left|W({\omega}_{k})\left(y({\omega}_{k})-{\displaystyle \sum _{i=1}^{{n}_{u}}\frac{{N}_{i}({\omega}_{k})}{D({\omega}_{k})}{u}_{i}({\omega}_{k})}\right)\right|}^{2}}$$
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 N_{i} is
the numerator corresponding to the ith input. y and u are
the measured output and input data, respectively. n_{f} and n_{u} are
the number of frequencies and inputs, and w is
the frequency. Rearranging the terms gives:
$$\underset{D,{N}_{i}}{\text{minimize}}{\displaystyle \sum _{k=1}^{{n}_{f}}{\left|\frac{W({\omega}_{k})}{D({\omega}_{k})}\left(D({\omega}_{k})y({\omega}_{k})-{\displaystyle \sum _{i=1}^{{n}_{u}}{N}_{i}({\omega}_{k}){u}_{i}({\omega}_{k})}\right)\right|}^{2}}$$
To perform the S-K iterations, the algorithm iteratively solves
$$\underset{{D}_{m},{N}_{i,m}}{\text{minimize}}{\displaystyle \sum _{k=1}^{{n}_{f}}{\left|\frac{W({\omega}_{k})}{{D}_{m-1}({\omega}_{k})}\left({D}_{m}({\omega}_{k})y({\omega}_{k})-{\displaystyle \sum _{i=1}^{{n}_{u}}{N}_{i,m}({\omega}_{k}){u}_{i}({\omega}_{k})}\right)\right|}^{2}}$$
where m is the current iteration, and D_{m-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 D_{m}(ω) and N_{i,m}(ω) for i =
1,2,...n_{u}. The iteration
is initialized by choosing D_{0}(ω) =
1.
The first iteration of the algorithm identifies D_{1}(ω).
The D_{1}(ω) and N_{i,1}(ω) polynomials
are expressed in monomial basis.
The second and following iterations express the polynomials D_{m}(ω) and N_{i,m}(ω) in
terms of orthogonal rational basis functions on the unit disk. These
basis functions have the form:
$${B}_{j,m}(\omega )=\left(\frac{\sqrt{1-{\left|{\lambda}_{j,m-1}\right|}^{2}}}{q-{\lambda}_{j,m-1}}\right){\displaystyle \prod _{r=0}^{j-1}\frac{1-{({\lambda}_{j,m-1})}^{*}q(\omega )}{q(\omega )-{\lambda}_{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 N_{i} and
denominator D parameters during S-K iterations.
Unlike S-K iterations, the basis functions B_{j}(ω) 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
.
When you add a regularization penalty to the loss
function using the Regularization
estimation option.
For an example about regularization, see Regularized Identification of Dynamic Systems.
You estimate a continuous-time model using discrete-time
data (see Discrete-Time Transfer Function Estimation Using Discrete-Time
Frequency-Domain Data).
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.
Discrete-Time Transfer Function Estimation Using Discrete-Time Frequency-Domain Data
The estimation algorithm is the same as for continuous-time
transfer function estimation using continuous-time frequency-domain
data, except discrete-time data is used.
If you are using the estimation algorithm from R2016a or earlier
(see Compatibility), the algorithm is the same as the algorithm
for discrete-time transfer function estimation using time-domain data.
Note:
The software does not support estimation of a discrete-time
transfer function using continuous-time frequency-domain data. |
Continuous-Time Transfer Function Estimation Using Discrete-Time Frequency-Domain Data
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.
For an initial model, init_sys
,
with:
the initial delay values are left unchanged.
Estimation of delays is often a difficult problem. You should
assess the presence and the value of a delay. To do so, use physical
insight of the process being modeled and functions such as arxstruc
, delayest
,
and impulseest
. For an example
of determining input delay, see Model Structure Selection: Determining Model Order and Input Delay.