LTI state-space models

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.

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 Identication 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)
• y: Measurements N×ny
• u: Control signal N×nu
• 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

sys, x0, opt = pem(data; nx, kwargs...)

System identification using the prediction-error method.

Arguments:

• data: iddata object containing y and u.
  • y: Measurements, either a matrix with time along dim 2, or a vector of vectors
  • u: Control signals, same structure as y
• nx: Number of poles in the estimated system. Thus number should be chosen as number of system poles plus number of poles in noise models for measurement noise and load disturbances.
• focus: Either :prediction or :simulation. If :simulation is chosen, a two stage problem is solved with prediction focus first, followed by a refinement for simulation focus.
• metric: A Function determining how the size of the residuals is measured, default sse (e'e), but any Function such as norm, e->sum(abs,e) or e -> e'Q*e could be used.
• regularizer(p)=0: function for regularization. The structure of p is detailed below
• solver Defaults to Optim.BFGS()
• stabilize_predictor=true: Modifies the estimated Kalman gain K in case A-KC is not stable by moving all unstable eigenvalues to the unit circle.
• difficult=false: If the identification problem appears to be difficult and ends up in a local minimum, set this flag to true to solve an initial global optimization problem to supply a good initial guess. This is expected to take some time.
• kwargs: additional keyword arguments are sent to Optim.Options.

Return values

• sys::StateSpaceNoise: identified system. Can be converted to StateSpace by convert(StateSpace, sys) or ss(sys), but this will discard the Kalman gain matrix, see innovation_form to obtain a predictor system.
• x0: Estimated initial state
• opt: Optimization problem structure. Contains info of the result of the optimization problem

Structure of parameter vector p

The parameter vector is of type ComponentVector and the fields A,B,K,x0 can be accessed as p.A etc. The internal storage is according to

A = size(nx,nx)
B = size(nx,nu)
K = size(nx,ny)
x0 = size(nx)
p = [A[:]; B[:]; K[:]; x0]

newpem(
    d,
    nx;
    zeroD = true,
    sys0 = subspaceid(d, nx; zeroD),
    focus = :prediction,
    optimizer = BFGS(
        alphaguess = LineSearches.InitialStatic(alpha = 1),
        linesearch = LineSearches.HagerZhang(),
    ),
    zerox0 = false,
    initx0 = false,
    store_trace = true,
    show_trace = true,
    show_every = 50,
    iterations = 10000,
    allow_f_increases = false,
    time_limit = 100,
    x_tol = 0,
    f_abstol = 0,
    g_tol = 1e-12,
    f_calls_limit = 0,
    g_calls_limit = 0,
)

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

Arguments:

• d: iddata
• nx: Model order
• zeroD: Force zero D matrix
• 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
• zerox0: Force initial state to zero.
• initx0: Estimate initial state once, otherwise at each iteration

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

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

era(d::AbstractIdData, r, m = 2r, n = 2r, l = 5r; λ=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

H = okid(d::AbstractIdData, r, l = 5r; λ=0)

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

Arguments:

• r: Model order
• l: Number of Markov parameters to estimate.
• λ: Regularization parameter