Transfer function estimation

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.

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 

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

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

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

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

yt,A = getARregressor(y::AbstractVector, na)

Returns values such that x = A\yt. See getARXregressor for more details.