Transfer function estimation

This page documents how to estimate transfer functions, sometimes called ARX or ARMAX models, i.e. models on any of the forms

\[\begin{aligned} G(z) &= \dfrac{B(z)}{A(z)} = \dfrac{b_m z^m + \dots + b_0}{z^n + a_{n-1} z^{n-1} + \dots + a_0} \\ Ay &= Bu + w \\ Ay &= Bu + Cw \\ Ay &= Bu + 1/D w \end{aligned}\]

This package estimates models in discrete time, but they may be converted to continuous-time models using the function d2c from ControlSystemsBase.jl.

The methods available are:

Functions

  • arx: Transfer-function estimation using least-squares fitting with a closed-form solution. Supports multiple datasets and custom estimator functions.
  • arma: Estimate an ARMA model (no control input).
  • ar: Estimate an AR model (no input).
  • arma_ssa Estimate an ARMA model with estimated noise as input (no control input).
  • plr: Transfer-function estimation (ARMAX model) using pseudo-linear regression. armax is an alias for this function. This method estimates a noise model as well.
  • arxar: Transfer-function estimation using generalized least-squares method. This method estimates a noise model as well.
  • getARXregressor/getARregressor: For low-level control over the estimation

See docstrings for further help.

Note

Most methods for estimation of transfer functions handle SISO, SIMO or MISO systems only. For estimation of MIMO systems, consider using state-space based methods and convert the result to a transfer function using tf after estimation if required.

Usage example:

N  = 2000     # Number of time steps
t  = 1:N
Δt = 1        # Sample time
u  = randn(N) # A random control input
G  = tf(0.8, [1,-0.9], 1)
y  = lsim(G,u,t)[1][:]
yn = y
d  = iddata(y,u,Δt)

na,nb = 1,1   # Number of polynomial coefficients

Gls = arx(d,na,nb,stochastic=false) # set stochastic to true to get a transfer function of MonteCarloMeasurements.Particles
@show Gls
# TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     0.8000000000000005
# --------------------------
# 1.0*z - 0.8999999999999997

As we can see, the model is perfectly recovered. In reality, the measurement signal is often affected by noise, in which case the estimation will suffer. To combat this, a few different options exist:

e  = randn(N)
yn = y + e    # Measurement signal with noise
d  = iddata(yn,u,Δt)

na,nb,nc = 1,1,1

Gls      = arx(d,na,nb, stochastic=true)     # Regular least-squares estimation
Gtls     = arx(d,na,nb, estimator=tls)       # Total least-squares estimation
Gwtls    = arx(d,na,nb, estimator=wtls_estimator(y,na,nb)) # Weighted Total least-squares estimation
Gplr, Gn = plr(d,na,nb,nc, initial_order=20) # Pseudo-linear regression
@show Gls; @show  Gtls; @show  Gwtls; @show  Gplr; @show  Gn;
# TransferFunction{ControlSystemsBase.SisoRational{MonteCarloMeasurements.Particles{Float64,500}}}
#     0.824 ± 0.029
# ---------------------
# 1.0*z - 0.713 ± 0.013

# Gtls = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     1.848908051191616
# -------------------------
# 1.0*z - 0.774385918070221

# Gwtls = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#    0.8180228878106678
# -------------------------
# 1.0*z - 0.891939152690534

# Gplr = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     0.8221837077656046
# --------------------------
# 1.0*z - 0.8896345125395438

# Gn = TransferFunction{ControlSystemsBase.SisoRational{Float64}}
#     0.9347035105826179
# --------------------------
# 1.0*z - 0.8896345125395438

We now see that the estimate using standard least-squares is heavily biased and it is wrongly certain about the estimate (notice the ± in the transfer function coefficients). Regular Total least-squares does not work well in this example, since not all variables in the regressor contain equally much noise. Weighted total least-squares does a reasonable job at recovering the true model. Pseudo-linear regression also fares okay, while simultaneously estimating a noise model. The helper function wtls_estimator(y,na,nb) returns a function that performs wtls using appropriately sized covariance matrices, based on the length of y and the model orders. Weighted total least-squares estimation is provided by TotalLeastSquares.jl. See the example notebooks for more details.

Uncertain transfer function with Particles coefficients can be used like any other model. Try, e.g., nyquistplot(Gls) to get a Nyquist plot with confidence bounds.

See also function arma for estimation of signal models without inputs.

A video using the function arx is available here:

Time-series modeling

Time-series modeling can be seen as special cases of transfer-function modeling where there are no control inputs. This package is primarily focused on control system identification, but we nevertheless provide two methods aimed at time-series estimation:

  • ar: Estimate an AR model (no input).
  • arma_ssa Estimate an ARMA model with estimated noise as input (no control input).

Code generation

To generate C-code for, e.g., simulating a system or filtering through an estimated transfer function, see SymbolicControlSystems.jl.

Transfer-function API

ControlSystemIdentification.arxFunction
Gtf = arx(d::AbstractIdData, na, nb; inputdelay = ones(Int, size(nb)), λ = 0, estimator=\, stochastic=false)

Fit a transfer Function to data using an ARX model and equation error minimization.

  • nb and na are the number of coefficients of the numerator and denominator polynomials.

Input delay can be added via inputdelay = d, which corresponds to an additional delay of z^-d. An inputdelay = 0 results in a direct term. The highest order of the B polynomial is given by nb + inputdelay - 1. λ > 0 can be provided for L₂ regularization. estimator defaults to \ (least squares), alternatives are estimator = tls for total least-squares estimation. arx(Δt,yn,u,na,nb, estimator=wtls_estimator(y,na,nb) is potentially more robust in the presence of heavy measurement noise. The number of free parameters is na+nb

  • stochastic: if true, returns a transfer function with uncertain parameters represented by MonteCarloMeasurements.Particles.

Supports MISO estimation by supplying an iddata with a matrix u, with nb = [nb₁, nb₂...] and optional inputdelay = [d₁, d₂...]

This function supports multiple datasets, provided as a vector of iddata objects.

source
ControlSystemIdentification.arFunction
ar(d::AbstractIdData, na; λ=0, estimator=\, scaleB=false, stochastic=false)

Estimate an AR transfer function G = 1/A, the AR process is defined as A(z⁻¹)y(t) = e(t)

Arguments:

  • d: IdData, see iddata
  • na: order of the model
  • λ: λ > 0 can be provided for L₂ regularization
  • estimator: e.g. \,tls,irls,rtls
  • scaleB: Whether or not to scale the numerator using the variance of the prediction error.
  • stochastic: if true, returns a transfer function with uncertain parameters represented by MonteCarloMeasurements.Particles.

Estimation of AR models using least-squares is known to struggle with heavy measurement noise, using estimator = tls can improve the result in this case.

Example

julia> N = 10000
10000

julia> e = [-0.2; zeros(N-1)] # noise e
10000-element Vector{Float64}:
[...]

julia> G = tf([1, 0], [1, -0.9], 1) # AR transfer function
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   1.0z
----------
1.0z - 0.9

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> y = lsim(G, e, 1:N)[1][:] # Get output of AR transfer function from input noise e
10000-element Vector{Float64}:
[...]

julia> Gest = ar(iddata(y), 1) # Estimate AR transfer function from output y
TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Float64}}
          1.0z
-------------------------
1.0z - 0.8999999999999998

Sample Time: 1.0 (seconds)
Discrete-time transfer function model

julia> G ≈ Gest # Test if estimation was correct
true

julia> eest = lsim(1/Gest, y, 1:N)[1][:] # recover the input noise e from output y and estimated transfer function Gest
10000-element Vector{Float64}:
[...]

julia> isapprox(eest, e, atol = eps()) # input noise correct recovered
true 
source
ControlSystemIdentification.armaFunction
model = arma(d::AbstractIdData, na, nc; initial_order=20, method=:ls)

Estimate a Autoregressive Moving Average model with na coefficients in the denominator and nc coefficients in the numerator. Returns the model and the estimated noise sequence driving the system.

Arguments:

  • d: iddata
  • initial_order: An initial AR model of this order is used to estimate the residuals
  • estimator: A function (A,y)->minimizeₓ(Ax-y) default is \ but another option is wtls_estimator(1:length(y)-initial_order,na,nc,ones(nc))

See also estimate_residuals

source
ControlSystemIdentification.plrFunction
G, Gn = plr(d::AbstractIdData,na,nb,nc; initial_order = 20)

Perform pseudo-linear regression to estimate a model on the form Ay = Bu + Cw The residual sequence is estimated by first estimating a high-order arx model, whereafter the estimated residual sequence is included in a second estimation problem. The return values are the estimated system model, and the estimated noise model. G and Gn will always have the same denominator polynomial.

armax is an alias for plr. See also pem, ar, arx and arxar

source
ControlSystemIdentification.arxarFunction
G, H, e = arxar(d::InputOutputData, na::Int, nb::Union{Int, Vector{Int}}, nd::Int)

Estimate the ARXAR model Ay = Bu + v, where v = He and H = 1/D, using generalized least-squares method. For more information see Söderström - Convergence properties of the generalized least squares identification method, 1974.

Arguments:

  • d: iddata
  • na: order of A
  • nb: number of coefficients in B, the order is determined by nb + inputdelay - 1. In MISO estimation it takes the form nb = [nb₁, nb₂...].
  • nd: order of D

Keyword Arguments:

  • H = nothing: prior knowledge about the AR noise model
  • inputdelay = ones(Int, size(nb)): optional delay of input, inputdelay = 0 results in a direct term, takes the form inputdelay = [d₁, d₂...] in MISO estimation
  • λ = 0: λ > 0 can be provided for L₂ regularization
  • estimator = \: e.g. \,tls,irls,rtls, the latter three require using TotalLeastSquares
  • δmin = 10e-4: Minimal change in the power of e, that specifies convergence.
  • iterations = 10: maximum number of iterations.
  • verbose = false: if true, more information is printed

See also plr, arx.

Example:

julia> N = 500 
500

julia> sim(G, u) = lsim(G, u, 1:N)[1][:]
sim (generic function with 1 method)

julia> A = tf([1, -0.8], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
1.0z - 0.8
----------
   1.0z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> B = tf([0, 1], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Int64}}
1
-
z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> G = minreal(B / A)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   1.0
----------
1.0z - 0.8

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> D = tf([1, 0.7], [1, 0], 1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
1.0z + 0.7
----------
   1.0z

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> H = 1 / D
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   1.0z
----------
1.0z + 0.7

Sample Time: 1 (seconds)
Discrete-time transfer function model

julia> u, e = randn(1, N), randn(1, N)
[...]

julia> y, v = sim(G, u), sim(H * (1/A), e) # simulate process
[...]

julia> d = iddata(y .+ v, u, 1)
InputOutput data of length 500 with 1 outputs and 1 inputs

julia> na, nb , nd = 1, 1, 1
(1, 1, 1)

julia> Gest, Hest, res = arxar(d, na, nb, nd)
(G = TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
   0.9987917259291642
-------------------------
1.0z - 0.7937837464682017

Sample Time: 1 (seconds)
Discrete-time transfer function model, H = TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
          1.0z
-------------------------
1.0z + 0.7019519225937721

Sample Time: 1 (seconds)
Discrete-time transfer function model, e = [...]
source
ControlSystemIdentification.getARXregressorFunction
getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb; inputdelay = ones(Int, size(nb)))

Returns a shortened output signal y and a regressor matrix A such that the least-squares ARX model estimate of order na,nb is y\A

Return a regressor matrix used to fit an ARX model on, e.g., the form A(z)y = B(z)u with output y and input u where the order of autoregression is na, the order of input moving average is nb and an optional input delay inputdelay. Caution, changing the input delay changes the order to nb + inputdelay - 1. An inputdelay = 0 results in a direct term.

Example

A     = [1,2*0.7*1,1] # A(z) coeffs
B     = [10,5] # B(z) coeffs
u     = randn(100) # Simulate 100 time steps with Gaussian input
y     = filt(B,A,u)
yr,A  = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
x     = A\yr # Estimate model polynomials
plot([yr A*x], lab=["Signal" "Prediction"])

For nonlinear ARX-models, see BasisFunctionExpansions.jl. See also arx

source

Video tutorials

Video tutorials performing transfer-function estimation are available here: