Exported functions and types

Index

ControlSystemIdentification.FRDType
FRD(w, r)

Represents frequency-response data. w holds the frequency vector and r the response. Methods defined on this type include

  • +-*
  • length, vec, sqrt
  • plot
  • feedback
  • freqvec
  • tfest to estimate a rational model
  • Indexing in the frequency domain using, e.g., G[1Hz : 5Hz], G[1rad : 5rad]

If r represents a MIMO frequency response, the dimensions are ny × nu × nω. freqresp returns a PermutedDimsArray whose .parent field follows this convention.

source
ControlSystemIdentification.N4SIDStateSpaceType
N4SIDStateSpace <: AbstractPredictionStateSpace

The result of statespace model estimation using the n4sid method.

Fields:

  • sys: estimated model in the form of a StateSpace object
  • Q: estimated covariance matrix of the states
  • R: estimated covariance matrix of the measurements
  • S: estimated cross covariance matrix between states and measurements
  • K: kalman observer gain
  • P: solution to the Riccatti equation
  • x: estimated state trajectory
  • s: singular values
  • fve: Fraction of variance explained by singular values
source
ControlSystemIdentification.PredictionStateSpaceType
PredictionStateSpace{T, ST <: AbstractStateSpace{T}, KT, QT, RT, ST2} <: AbstractPredictionStateSpace{T}
PredictionStateSpace(sys, K, Q=nothing, R=nothing, S=nothing)

A statespace type that contains an additional Kalman-filter model for prediction purposes.

Arguments:

  • sys: DESCRIPTION
  • K: Infinite-horizon Kalman gain
  • Q = nothing: Dynamics covariance
  • R = nothing: Measurement covariance
  • S = nothing: Cross-covariance
source
ControlSystemIdentification.arMethod
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}, ControlSystems.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}, ControlSystems.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.armaMethod
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.arma_ssaMethod
arma_ssa(d::AbstractIdData, na, nc; L=nothing, estimator=\, robust=false)

Fit arma models using Singular Spectrum Analysis (SSA). A low-rank factorization (svd or robust svd) is performed on the data in order to decompose the signal and the noise. The noise is then used as input to fit an arma model.

Arguments:

  • d: iddata
  • na: number of denominator parameters
  • nc: number of numerator parameters
  • L: length of the lag-embedding used to separate signal and noise. nothing corresponds to automatic selection.
  • estimator: The function to solve the least squares problem. Examples \,tls,irls,rtls.
  • robust: Use robust PCA to be resistant to outliers.
source
ControlSystemIdentification.arxMethod
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₂...]

source
ControlSystemIdentification.arxarMethod
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

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}, ControlSystems.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}, ControlSystems.SisoRational{Int64}}
1
-
z

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

julia> G = minreal(B / A)
TransferFunction{Discrete{Int64}, ControlSystems.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}, ControlSystems.SisoRational{Float64}}
1.0z + 0.7
----------
   1.0z

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

julia> H = 1 / D
TransferFunction{Discrete{Int64}, ControlSystems.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}, ControlSystems.SisoRational{Float64}}
   0.9987917259291642
-------------------------
1.0z - 0.7937837464682017

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

Sample Time: 1 (seconds)
Discrete-time transfer function model, e = [...]
source
ControlSystemIdentification.coherenceMethod
κ = coherence(d; n = length(d)÷10, noverlap = n÷2, window=hamming)

Calculates the magnitude-squared coherence Function. κ close to 1 indicates a good explainability of energy in the output signal by energy in the input signal. κ << 1 indicates that either the system is nonlinear, or a strong noise contributes to the output energy. κ: Coherence function (not squared) N: Noise model

source
ControlSystemIdentification.coherenceplotFunction
coherenceplot(d, [(;n=..., noverlap=...); hz=false)

Calculates and plots the (squared) coherence Function κ. κ close to 1 indicates a good explainability of energy in the output signal by energy in the input signal. κ << 1 indicates that either the system is nonlinear, or a strong noise contributes to the output energy.

hz indicates Hertz instead of rad/s

Keyword arguments to coherence are supplied as a named tuple as a second positional argument .

source
ControlSystemIdentification.crosscorplotFunction
crosscorplot(data, [lags])
crosscorplot(u, y, Ts, [lags])

Plot the cross correlation between input and output for lags that default to 10% of the length of the dataset on the negative side and 50% on the positive side but no more than 100 on each side.

source
ControlSystemIdentification.eraFunction
era(d::AbstractIdData, r, m = 2r, n = 2r, l = 5r; p = l, λ=0)

Eigenvalue realization algorithm. Uses okid to find the Markov parameters as an initial step.

Arguments:

  • r: Model order
  • l: Number of Markov parameters to estimate.
  • λ: Regularization parameter
  • p: Optionally, delete the first p columns in the internal Hankel matrices to account for initial conditions != 0. If x0 != 0, for era, p defaults to l, while when calling okid directly, p defaults to 0.
source
ControlSystemIdentification.eraMethod
era(YY::AbstractArray{<:Any, 3}, Ts, r::Int, m::Int, n::Int)

Eigenvalue realization algorithm.

Arguments:

  • YY: Markov parameters (impulse response) size n_out×n_in×n_time
  • Ts: Sample time
  • r: Model order
  • m: Number of rows in Hankel matrix
  • n: Number of columns in Hankel matrix
source
ControlSystemIdentification.estimate_x0Function
estimate_x0(sys, d, n = min(length(d), 3 * slowest_time_constant(sys)); fixed = fill(NaN, sys.nx)

Estimate the initial state of the system

Arguments:

  • d: iddata
  • n: Number of samples to use.
  • fixed: If a vector of the same length as x0 is provided, finite values indicate fixed values that are not to be estimated, while nonfinite values are free.

Example

sys   = ssrand(2,3,4, Ts=1)
x0    = randn(sys.nx)
u     = randn(sys.nu, 100)
y,t,x = lsim(sys, u; x0)
d     = iddata(y, u, 1)
x0h   = estimate_x0(sys, d, 8, fixed=[Inf, x0[2], Inf, Inf])
x0h[2] == x0[2] # Should be exact equality
norm(x0-x0h)    # Should be small
source
ControlSystemIdentification.find_similarity_transformFunction
find_similarity_transform(sys1, sys2)

Find T such that ControlSystems.similarity_transform(sys1, T) == sys2

Ref: Minimal state-space realization in linear system theory: an overview, B. De Schutter

julia> T = randn(3,3);

julia> sys1 = ssrand(1,1,3);

julia> sys2 = ControlSystems.similarity_transform(sys1, T);

julia> T2 = find_similarity_transform(sys1, sys2);

julia> T2 ≈ T
true
source
ControlSystemIdentification.getARXregressorMethod
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
ControlSystemIdentification.iddataFunction
iddata(y,       Ts = nothing)
iddata(y, u,    Ts = nothing)
iddata(y, u, x, Ts = nothing)

Create a time-domain identification data object.

Arguments

  • y::AbstractArray: output data (required)
  • u::AbstractArray: input data (if available)
  • x::AbstractArray: state data (if available)
  • Ts::Union{Real,Nothing} = nothing: optional sample time

If the time-series are multivariate, time is in the last dimension.

Operations on iddata

  • detrend
  • prefilter
  • resample
  • append two along the time dimension [d1 d2]
  • index time series d[output_index, input_index]
  • index the time axis with indices d[time_indices]
  • index the time axis with seconds d[3Sec:12Sec] (using ControlSystemIdentification: Sec)
  • access number of inputs, outputs and sample time: d.nu, d.ny, d.Ts
  • access the time time vector d.t
  • premultiply to scale outputs C * d
  • postmultiply to scale inputs d * B
  • writedlm
  • ramp_in, ramp_out
  • plot
  • specplot
  • crosscorplot

Examples

julia> iddata(randn(10))
Output data of length 10 with 1 outputs

julia> iddata(randn(10), randn(10), 1)
InputOutput data of length 10 with 1 outputs and 1 inputs

julia> d = iddata(randn(2, 10), randn(3, 10), 0.1)
InputOutput data of length 10 with 2 outputs and 3 inputs

julia> [d d] # Concatenate along time
InputOutput data of length 20 with 2 outputs and 3 inputs

julia> d[1:3]
InputOutput data of length 3 with 2 outputs and 3 inputs

julia> d.nu
3

julia> d.t # access time vector
0.0:0.1:0.9
source
ControlSystemIdentification.impulseestMethod
ir, t, Σ = impulseest(d::AbstractIdData, n; λ=0, estimator=ls)

Estimates the system impulse response by fitting an n:th order FIR model. Returns impulse-response estimate, time vector and covariance matrix. See also impulseestplot

source
ControlSystemIdentification.minimum_phaseMethod
minimum_phase(G)

Move zeros and poles of G from the unstable half plane to the stable. If G is a statespace system, it's converted to a transfer function first. This can incur loss of precision.

source
ControlSystemIdentification.model_spectrumMethod
model_spectrum(f, h, args...; kwargs...)

Arguments:

  • f: the model-estimation function, e.g., ar,arma
  • h: The sample time
  • args: arguments to f
  • kwargs: keyword arguments to f

Example:

using ControlSystemIdentification, DSP
T = 1000
s = sin.((1:T) .* 2pi/10)
S1 = spectrogram(s,window=hanning)
estimator = model_spectrum(ar,1,2)
S2 = spectrogram(s,estimator,window=rect)
plot(plot(S1),plot(S2)) # Requires the package LPVSpectral.jl
source
ControlSystemIdentification.n4sidMethod
res = n4sid(data, r=:auto; verbose=false)

Estimate a statespace model using the n4sid method. Returns an object of type N4SIDStateSpace where the model is accessed as res.sys.

Implements the simplified algorithm (alg 2) from "N4SID: Subspace Algorithms for the Identification of Combined Deterministic Stochastic Systems" PETER VAN OVERSCHEE and BART DE MOOR

The frequency weighting is borrowing ideas from "Frequency Weighted Subspace Based System Identification in the Frequency Domain", Tomas McKelvey 1996. In particular, we apply the output frequency weight matrix (Fy) as it appears in eqs. (16)-(18).

Arguments:

  • data: Identification data data = iddata(y,u)
  • r: Rank of the model (model order)
  • verbose: Print stuff?
  • Wf: A frequency-domain model of measurement disturbances. To focus the attention of the model on a narrow frequency band, try something like Wf = Bandstop(lower, upper, fs=1/Ts) to indicate that there are disturbances outside this band.
  • i: Algorithm parameter, generally no need to tune this
  • γ: Set this to a value between (0,1) to stabilize unstable models such that the largest eigenvalue has magnitude γ.
  • zeroD: defaults to false

See also the newer implementation subspaceid which allows you to choose between different weightings (n4sid being one of them). A more accurate prediciton model can sometimes be obtained using newpem, which is also unbiased for closed-loop data.

source
ControlSystemIdentification.newpemMethod
sys, x0, res = newpem(
    d,
    nx;
    zeroD  = true,
    focus  = :prediction,
    stable = true,
    sys0   = subspaceid(d, nx; zeroD, focus, stable),
    metric = abs2,
    regularizer = (p, P) -> 0,
    optimizer = BFGS(
        linesearch = LineSearches.BackTracking(),
    ),
    store_trace = true,
    show_trace  = true,
    show_every  = 50,
    iterations  = 10000,
    time_limit  = 100,
    x_tol       = 0,
    f_abstol    = 0,
    g_tol       = 1e-12,
    f_calls_limit = 0,
    g_calls_limit = 0,
    allow_f_increases = false,
)

A new implementation of the prediction-error method (PEM). Note that this is an experimental implementation and subject to breaking changes not respecting semver.

The prediction-error method is an iterative, gradient-based optimization problem, as such, it can be extra sensitive to signal scaling, and it's recommended to perform scaling to d before estimation, e.g., by pre and post-multiplying with diagonal matrices d̃ = Dy*d*Du, and apply the inverse scaling to the resulting system. In this case, we have

\[D_y y = G̃ D_u u ↔ y = D_y^{-1} G̃ D_u u\]

hence G = Dy \ G̃ * Du where $ G̃ $ is the plant estimated for the scaled iddata.

Arguments:

  • d: iddata
  • nx: Model order
  • zeroD: Force zero D matrix
  • stable if true, stability of the estimated system will be enforced by eigenvalue reflection using schur_stab with ϵ=1/100 (default). If stable is a real value, the value is used instead of the default ϵ.
  • sys0: Initial guess, if non provided, subspaceid is used as initial guess.
  • focus: prediction or :simulation. If :simulation, hte K matrix will be zero.
  • optimizer: One of Optim's optimizers
  • metric: The metric used to measure residuals. Try, e.g., abs for better resistance to outliers.

The rest of the arguments are related to Optim.Options.

  • regularizer: A function of the parameter vector and the corresponding PredictionStateSpace/StateSpace system that can be used to regularize the estimate.

Example

using ControlSystemIdentification, ControlSystems, Plots
G = DemoSystems.doylesat()
T = 1000  # Number of time steps
Ts = 0.01 # Sample time
sys = c2d(G, Ts)
nx = sys.nx
nu = sys.nu
ny = sys.ny
x0 = zeros(nx) # actual initial state
sim(sys, u, x0 = x0) = lsim(sys, u; x0)[1]

σy = 1e-1 # Noise covariance

u  = randn(nu, T)
y  = sim(sys, u, x0)
yn = y .+ σy .* randn.() # Add measurement noise
d  = iddata(yn, u, Ts)

sysh, x0h, opt = ControlSystemIdentification.newpem(d, nx, show_every=10)

plot(
    bodeplot([sys, sysh]),
    predplot(sysh, d, x0h), # Include the estimated initial state in the prediction
)

Extended help

This implementation uses a tridiagonal parametrization of the A-matrix that has been shown to be favourable from an optimization perspective.¹ The initial guess sys0 is automatically transformed to a special tridiagonal modal form. [1]: Mckelvey, Tomas & Helmersson, Anders. (1997). State-space parametrizations of multivariable linear systems using tridiagonal matrix forms.

The parameter vector used in the optimizaiton takes the following form

p = [trivec(A); vec(B); vec(C); vec(D); vec(K); vec(x0)]

Where ControlSystemIdentification.trivec vectorizes the -1,0,1 diagonals of A. If focus = :simulation, K is omitted, and if zeroD = true, D is omitted.

source
ControlSystemIdentification.noise_modelMethod
noise_model(sys::AbstractPredictionStateSpace)

Return a model of the noise driving the system, v, in

\[x' = Ax + Bu + Kv\\ y = Cx + Du + v\]

The model neglects u and is given by

\[x' = Ax + Kv\\ y = Cx + v\]

Also called the "innovation form". This function calls ControlSystems.innovation_form.

source
ControlSystemIdentification.okidFunction
H = okid(d::AbstractIdData, nx, l = 5nx; p = 1, λ=0, estimator = /)

Observer Kalman filter identification. Returns the Markov parameters H size n_out×n_in×l+1

Arguments:

  • nx: Model order
  • l: Number of Markov parameters to estimate.
  • λ: Regularization parameter
  • p: Optionally, delete the first p columns in the internal Hankel matrices to account for initial conditions != 0. If x0 != 0, try setting p around the same value as l.
source
ControlSystemIdentification.plrMethod
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

source
ControlSystemIdentification.predplotFunction
predplot(sys, data, x0=nothing; ploty=true, plote=false, h=1, sysname="")

Plot system simulation and measured output to compare them. ploty determines whether or not to plot the measured signal plote determines whether or not to plot the residual h is the prediction horizon.

source
ControlSystemIdentification.prefilterMethod
prefilter(d::AbstractIdData, responsetype::FilterType)

Filter both input and output of the identification data using zero-phase filtering (filtfilt). Since both input and output is filtered, linear identification will not be affected in any other way than to focus the fit on the selected frequency range, i.e. the range that has high gain in the provided filter. Note, if the system that generated d is nonlinear, identification might be severely impacted by this transformation. Verify linearity with, e.g., coherenceplot.

source
ControlSystemIdentification.prefilterMethod
prefilter(d::AbstractIdData, l::Number, u::Number)

Filter input and output with a bandpass filter between l and u Hz. If l = 0 a lowpass filter will be used, and if u = Inf a highpass filter will be used.

source
ControlSystemIdentification.schur_stabMethod
schur_stab(A::AbstractMatrix{T}, ϵ = 0.01)

Stabilize the eigenvalues of discrete-time matrix A by transforming A to complex Schur form and projecting unstable eigenvalues 1-ϵ < λ ≤ 2 into the unit disc. Eigenvalues > 2 are set to 0.

source
ControlSystemIdentification.simplotFunction
simplot(sys, data, x0=nothing; ploty=true, plote=false, sysname="")

Plot system simulation and measured output to compare them. ploty determines whether or not to plot the measured signal plote determines whether or not to plot the residual

source
ControlSystemIdentification.subspaceidMethod
subspaceid(frd::FRD, args...; estimate_x0 = false, bilinear_transform = false, kwargs...)

If a frequency-reponse data object is supplied

  • The FRD will be automatically converted to an InputOutputFreqData
  • estimate_x0 is by default set to 0.
  • bilinear_transform transform the frequency vector to discrete time, see note below.

Note: if the frequency-response data comes from a frequency-response analysis, a bilinear transform of the data is required before estimation. This transform will be applied if bilinear_transform = true.

source
ControlSystemIdentification.subspaceidMethod
subspaceid(
    data::InputOutputData,
    nx = :auto;
    verbose = false,
    r = nx === :auto ? min(length(data) ÷ 20, 20) : nx + 10, # the maximal prediction horizon used
    s1 = r, # number of past outputs
    s2 = r, # number of past inputs
    W = :MOESP,
    zeroD = false,
    stable = true, 
    focus = :prediction,
    svd::F1 = svd!,
    scaleU = true,
    Aestimator::F2 = \,
    Bestimator::F3 = \,
    weights = nothing,
)

Estimate a state-space model using subspace-based identification.

Ref: Ljung, Theory for the user.

Arguments:

  • data: Identification data iddata
  • nx: Rank of the model (model order)
  • verbose: Print stuff?
  • r: Prediction horizon. The model may perform better on simulation if this is made longer, at the expense of more computation time.
  • s1: past horizon of outputs
  • s2: past horizon of inputs
  • W: Weight type, choose between :MOESP, :CVA, :N4SID, :IVM
  • zeroD: Force the D matrix to be zero.
  • stable: Stabilize unstable system using eigenvalue reflection.
  • focus: :prediction or simulation
  • svd: The function to use for svd
  • scaleU: Rescale the input channels to have the same energy.
  • Aestimator: Estimator function used to estimate A,C.
  • Bestimator: Estimator function used to estimate B,D.
  • weights: A vector of weights can be provided if the Bestimator is wls.

Extended help

A more accurate prediciton model can sometimes be obtained using newpem, which is also unbiased for closed-loop data (subspaceid is biased for closed-loop data, see example in the docs). The prediction-error method is iterative and generally more expensive than subspaceid, and uses this function (by default) to form the initial guess for the optimization.

source
ControlSystemIdentification.subspaceidMethod
subspaceid(data::InputOutputFreqData,
    Ts = data.Ts,
    nx = :auto;
    cont = false,
    verbose = false,
    r = nx === :auto ? min(length(data) ÷ 20, 20) : 2nx, # Internal model order
    zeroD = false,
    estimate_x0 = true,
    stable = true, 
    svd = svd!,
    Aestimator = \,
    Bestimator = \,
    weights = nothing
)

Estimate a state-space model using subspace-based identification in the frequency domain.

See the docs for an example.

Arguments:

  • data: A frequency-domain identification data object.
  • Ts: Sample time at which the data was collected
  • nx: Desired model order, an interer or :auto.
  • cont: Return a continuous-time model? A bilinear transformation is used to convert the estimated discrete-time model, see function d2c.
  • verbose: Print stuff?
  • r: Internal model order, must be ≥ nx.
  • zeroD: Force the D matrix to be zero.
  • estimate_x0: Esimation of extra parameters to account for initial conditions. This may be required if the data comes from the fft of time-domain data, but may not be required if the data is collected using frequency-response analysis with exactly periodic input and proper handling of transients.
  • stable: For the model to be stable (uses schur_stab).
  • svd: The svd function to use.
  • Aestimator: The estimator of the A matrix (and initial C-matrix).
  • Bestimator: The estimator of B/D and C/D matrices.
  • weights: An optional vector of frequency weights of the same length as the number of frequencies in `data.
source
ControlSystemIdentification.tfestFunction
H, N = tfest(data, σ = 0.05)

Estimate a transfer function model using the Correlogram approach. Both H and N are of type FRD (frequency-response data).

  • σ determines the width of the Gaussian window applied to the estimated correlation functions before FFT. A larger σ implies less smoothing.
  • H = Syu/Suu Process transfer function
  • N = Sy - |Syu|²/Suu Noise PSD
source
ControlSystemIdentification.tfestFunction
tfest(
    data::FRD,
    p0,
    link = log ∘ abs;
    freq_weight = sqrt(data.w[1]*data.w[end]),
    refine = true,
    opt = BFGS(),
    opts = Optim.Options(
        store_trace       = true,
        show_trace        = true,
        show_every        = 1,
        iterations        = 100,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 0,
        f_tol             = 0,
        g_tol             = 1e-8,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
    ),
)

Fit a parametric transfer function to frequency-domain data.

The initial pahse of the optimization solves

\[\operatorname{minimize}_{B,A}{|| B/l - A||}\]

and the second stage (if refine=true) solves

\[\operatorname{minimize}_{B,A}{|| \text{link}\left(\dfrac{B}{A}\right) - \text{link}\left(l\right)||}\]

(abs2(link(B/A) - link(l)))

Arguments:

  • data: An FRD onbject with frequency domain data.
  • p0: Initial parameter guess. Can be a NamedTuple or ComponentVector with fields b,a specifying numerator and denominator as they appear in the call to tf, i.e., (b = [1.0], a = [1.0,1.0,1.0]). Can also be an instace of TransferFunction.
  • link: By default, phase information is discarded in the fitting. To include phase, change to link = log.
  • freq_weight: Apply weighting with the inverse frequency. The value determines the cutoff frequency before which the weight is constant, after which the weight decreases linearly. Defaults to the geometric mean of the smallest and largest frequency.
  • refine: Indicate whether or not a second optimization stage is performed to refine the results of the first.
  • opt: The Optim optimizer to use.
  • opts: Optim.Options controlling the solver options.

See also minimum_phase to transform a possibly non-minimum phase system to minimum phase.

source
ControlSystemIdentification.tfestMethod
tfest(data::FRD, basis::AbstractStateSpace; 
    freq_weight = 1 ./ (data.w .+ data.w[2]),
    opt = BFGS(),
    metric::M = abs2,
    opts = Optim.Options(
        store_trace       = true,
        show_trace        = true,
        show_every        = 50,
        iterations        = 1000000,
        allow_f_increases = false,
        time_limit        = 100,
        x_tol             = 1e-5,
        f_tol             = 0,
        g_tol             = 1e-8,
        f_calls_limit     = 0,
        g_calls_limit     = 0,
)

Fit a parametric transfer function to frequency-domain data using a pre-specified basis.

Arguments:

  • data: An FRD onbject with frequency domain data.

function kautz(a::AbstractVector)

  • basis: A basis for the estimation. See, e.g., laguerre, laguerre_oo, kautz
  • freq_weight: A vector of weights per frequency. The default is approximately 1/f.
  • opt: The Optim optimizer to use.
  • opts: Optim.Options controlling the solver options.
source
ControlSystemIdentification.wtls_estimatorFunction
wtls_estimator(y,na,nb, σu=0)

Create an estimator function for estimation of arx models in the presence of measurement noise. If the noise variance on the input σu (model errors) is known, this can be specified for increased accuracy.

source
DSP.Filters.resampleMethod
DSP.resample(sys::AbstractStateSpace{<:Discrete}, Qd::AbstractMatrix, newh::Real)

Change sample time of covariance matrix Qd beloning to sys to newh. This function does not handle the measurement covariance, how to do this depends on context. If the faster sampled signal has the same measurement noise, no change should be made. If the slower sampled signal was downsampled with filtering, the measurement covariance should be increased if the system is changed to a faster sample rate. To maintain the frequency response of the system, the measurement covariance should be modified accordinly.

Arguments:

  • sys: A discrete-time system that has dynamics noise covariance matric Qd.
  • Qd: Covariance matrix of dynamics noise.
  • newh: The new sample time.
source
DSP.Filters.resampleMethod
dr = resample(d::InputOutputData, f)

Resample iddata d with fraction f, e.g., f = fs_new / fs_original.

source
StatsAPI.predictMethod
predict(ARX::TransferFunction, d::InputOutputData)

One step ahead prediction for an ARX process. The length of the returned prediction is length(d) - max(na, nb)

Example:

julia> predict(tf(1, [1, -1], 1), iddata(1:10, 1:10))
9-element Vector{Int64}:
  2
  4
  6
  8
 10
 12
 14
 16
 18
source
StatsAPI.residualsMethod
residuals(ARX::TransferFunction, d::InputOutputData)

Calculates the residuals v = Ay - Bu of an ARX process and InputOutputData d. The length of the returned residuals is length(d) - max(na, nb)

Example:

julia> ARX = tf(1, [1, -1], 1)
TransferFunction{Discrete{Int64}, ControlSystems.SisoRational{Int64}}
  1
-----
z - 1

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

julia> u = 1:5
1:5

julia> y = lsim(ARX, u, 1:5)[1][:]
5-element Vector{Float64}:
  0.0
  1.0
  3.0
  6.0
 10.0

julia> d = iddata(y, u)
InputOutput data of length 5 with 1 outputs and 1 inputs

julia> residuals(ARX, d)
4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
source
DelimitedFiles.writedlmFunction
DelimitedFiles.writedlm(io::IO, d::AbstractIdData, args...; kwargs...)

Write identification data to disk.

source
ControlSystems.c2dFunction
c2d(w::AbstractVector{<:Real}, Ts; w_prewarp = 0)
c2d(frd::FRD, Ts; w_prewarp = 0)

Transform continuous-time frequency vector w or frequency-response data frd from continuous to discrete time using a bilinear (Tustin) transform. This is useful in cases where a frequency response is obtained through frequency-response analysis, and the function subspaceidf is to be used.

source
Qd = ControlSystems.c2d(sys::StateSpace{Discrete}, Q::Matrix)
Qd, Rd = ControlSystems.c2d(sys::StateSpace{Discrete}, Q::Matrix, R::Matrix)

Sample a continuous-time covariance matrix to fit the provided discrete-time system. The measurement covariance R may also be provided

The method used comes from theorem 5 in the reference below.

Ref: "Discrete-time Solutions to the Continuous-time Differential Lyapunov Equation With Applications to Kalman Filtering", Patrik Axelsson and Fredrik Gustafsson

On singular covariance matrices: The traditional double integrator with covariance matrix Q = diagm([0,σ²]) can not be sampled with this method. Instead, the input matrix ("Cholesky factor") of Q must be manually kept track of, e.g., the noise of variance σ² enters like N = [0, 1] which is sampled using ZoH and becomes Nd = [1/2 Ts^2; Ts] which results in the covariance matrix σ² * Nd * Nd'.

source