API Reference

kf, x_sym, ps, iosys, mats, prob = KalmanFilter(model::System, inputs, outputs; disturbance_inputs, Ts, R1, R2, x0map=[], pmap=[], σ0 = 1e-4, init=false, static=true, split = true, simplify=true, discretize = true, parametric = false, kwargs...)

Construct a Kalman filter for a linear MTK ODESystem. No check is performed to verify that the system is truly linear, if it is nonlinear, it will be linearized.

Returns:

• kf: A Kalman filter. If parametric=true, the A,B,C,D,R1,R2 fields are all functions of (x,u,p,t), otherwise they are matrices that are evaluated at the x0map, pmap values.
• x_sym: The symbolic state variables of the system.
• ps: The symbolic parameters of the system.
• iosys: The simplified MTK System
• mats: A named tuple containing the symbolic system matrices (A,B,C,D,Bw,Dw), where Bw and Dw are the input matrices corresponding to the disturbance inputs.
• prob: A StateEstimationProblem object. This problem object does not play quite the same role as when using Unscented or Extended Kalman filters since the filter is created already by this constructor, but the problem object can still be useful for inspecting the simplified MTK system, to create EstimatedOutput objects and to make use of the symbolic indexing functionality.

Arguments:

• model: An MTK System model, this model must not have undergone structural simplification.
• inputs: The inputs to the dynamical system, a vector of symbolic variables.
• outputs: The outputs of the dynamical system, a vector of symbolic variables.
• disturbance_inputs: The disturbance inputs to the dynamical system, a vector of symbolic variables. These disturbance inputs indicate where dynamics noise $w$ enters the system. The probability distribution with covariance $R_1$ is defined over these variables.
• Ts: The discretization time step.
• R1: The covariance matrix of the dynamics noise (disturbance inputs) $w$.
• R2: The covariance matrix of the measurement noise $e$.
• x0map: A dictionary mapping symbolic variables to their initial values. If a variable is not provided, it is assumed to be initialized to zero. The value can be a scalar number, in which case the covariance of the initial state is set to σ0^2*I(nx), and the value can be a Distributions.Normal, in which case the provided distributions are used as the distribution of the initial state. When passing distributions, all state variables must be provided values.
• pmap: A dictionary mapping symbolic variables to their values.
• σ0: The standard deviation of the initial state. This is used when x0map is not provided.
• init: If true, the initial state is computed using an initialization problem. If false, the initial state is computed using the get_u0 function.
• static: If true, static arrays are used for the state and covariance matrix. This can improve performance for small systems.
• split: Passed to mtkcompile, see the documentation there.
• simplify: Passed to mtkcompile, see the documentation there.
• discretize: If true, the system is discretized using zero-order hold. If false, matrices/functions are generated for the continuous-time system, in which case the user must handle discretization themselves (filtering with a continuous-time system without discretization will yield nonsensical results).
• parametric_A: If true, the A field of the returned filter is a function of (x,u,p,t), otherwise it is a matrix that is evaluated at the x0map, pmap values.
• parametric_B: If true, the B field of the returned filter is a function of (x,u,p,t), otherwise it is a matrix that is evaluated at the x0map, pmap values.
• parametric_C: If true, the C field of the returned filter is a function of (x,u,p,t), otherwise it is a matrix that is evaluated at the x0map, pmap values.
• parametric_D: If true, the D field of the returned filter is a function of (x,u,p,t), otherwise it is a matrix that is evaluated at the x0map, pmap values.
• parametric_R1: If true, the R1 field of the returned filter is a function of (x,u,p,t), otherwise it is a matrix that is evaluated at the x0map, pmap values.
• parametric_R2: If true, the R2 field of the returned filter is a function of (x,u,p,t), otherwise it is a matrix that is evaluated at the x0map, pmap values.
• tuplify: If true, the parameter vector p is returned as a tuple instead of an array. This can improve performance for filters with a small number of parameters of heterogeneous types.
• warn_initialize_determined: Passed on to the internal InitializationProblem construction, see StateEstimationProblem.
• kwargs: Additional keyword arguments passed to mtkcompile.

EstimatedOutput(kf, prob, sym)

Create an output function that can be called like

g(x::Vector,    u, p, t)     # Compute an output
g(xR::MvNormal, u, p, t)     # Compute an output distribution given input distribution xR
g(kf,           u, p, t)     # Compute an output distribution given the current state of an AbstractKalmanFilter

Arguments:

• kf: A Kalman type filter
• prob: A StateEstimationProblem object
• sym: A symbolic expression or vector of symbolic expressions that the function should output.

StateEstimationProblem(model, inputs, outputs; disturbance_inputs, discretization, Ts, df, dg, x0map=[], pmap=[], init=false)

A structure representing a state-estimation problem.

Arguments:

• model: An MTK ODESystem model, this model must not have undergone structural simplification.
• inputs: The inputs to the dynamical system, a vector of symbolic variables that must be of type @variables.
• outputs: The outputs of the dynamical system, a vector of symbolic variables that must be of type @variables.
• disturbance_inputs: The disturbance inputs to the dynamical system, a vector of symbolic variables that must be of type @variables. These disturbance inputs indicate where dynamics noise $w$ enters the system. The probability distribution $d_f$ is defined over these variables.
• discretization: A function discretization(f_cont, Ts, x_inds, alg_inds, nu) = f_disc that takes a continuous-time dynamics function f_cont(x,u,p,t) and returns a discrete-time dynamics function f_disc(x,u,p,t). x_inds is the indices of differential state variables, alg_inds is the indices of algebraic variables, and nu is the number of inputs.
• Ts: The discretization time step.
• df: The probability distribution of the dynamics noise $w$. When using Kalman-type estimators, this must be a MvNormal or SimpleMvNormal distribution.
• dg: The probability distribution of the measurement noise $e$. When using Kalman-type estimators, this must be a MvNormal or SimpleMvNormal distribution.
• x0map: A dictionary mapping symbolic variables to their initial values. If a variable is not provided, it is assumed to be initialized to zero. The value can be a scalar number, in which case the covariance of the initial state is set to σ0^2*I(nx), and the value can be a Distributions.Normal, in which case the provided distributions are used as the distribution of the initial state. When passing distributions, all state variables must be provided values.
• σ0: The standard deviation of the initial state. This is used when x0map is not provided or when the values in x0map are scalars.
• pmap: A dictionary mapping symbolic variables to their values. If a variable is not provided, it is assumed to be initialized to zero.
• init: If true, the initial state is computed using an initialization problem. If false, the initial state is computed using the get_u0 function.
• warn_initialize_determined: Passed on to the internal InitializationProblem/ODEProblem construction, default matches MTK's true.
• xscalemap: A dictionary mapping state variables to scaling factors. This is used to scale the state variables during integration to improve numerical stability. If a variable is not provided, it is assumed to have a scaling factor of 1.0. If provided, discretization is a function with signature discretization(f_cont, Ts, x_inds, alg_inds, nu, scale_x) where scale_x is a vector of scaling factors for the state variables.

Usage:

Pseudocode

prob      = StateEstimationProblem(...)
kf        = get_filter(prob, ExtendedKalmanFilter)      # or UnscentedKalmanFilter
filtersol = forward_trajectory(kf, u, y)
sol       = StateEstimationSolution(filtersol, prob)   # Package into higher-level solution object
plot(sol, idxs=[prob.state; prob.outputs; prob.inputs]) # Plot the solution

StateEstimationSolution(kfsol, prob)

A solution object that provides symbolic indexing to a KalmanFilteringSolution object.

Fields:

• sol: a KalmanFilteringSolution object.
• prob: a StateEstimationProblem object.

Example

sol = StateEstimationSolution(kfsol, prob)
sol[model.x]                 # Index with a variable
sol[model.y^2]               # Index with an expression
sol[model.y^2, dist=true]    # Obtain the posterior probability distribution of the provided expression
sol[model.y^2, Nsamples=100] # Draw 100 samples from the posterior distribution of the provided expression

get_filter(prob::StateEstimationProblem, ::Type{DAEUnscentedKalmanFilter}; constraint_solver, p = prob.p, kwargs...)

Instantiate a DAEUnscentedKalmanFilter from a state-estimation problem built around an MTK model with algebraic equations. The package auto-generates get_x_z, build_xz, and the algebraic residual callback from the equation/unknown splitting that MTK produced during mtkcompile.

The user must supply:

• constraint_solver: a callable (f, z0) -> z that solves f(z) ≈ 0. Typically LowLevelParticleFilters.scimlbase_solver(SimpleNewtonRaphson(); reltol=1e-12).

The discretization callback passed to StateEstimationProblem must be a DAE-aware integrator such as SeeToDee.Trapezoidal(f, Ts, x_inds, a_inds, nu) or SeeToDee.SimpleColloc(...) — the resulting prob.f is forwarded directly as the DAE UKF's dynamics.

The parameter object p defaults to prob.p. It is stored in the filter and used to compute Bw below. See remake for replacement of the parameters of the problem.

The number of disturbance inputs nw may differ from the differential-state count nx_diff:

• When nw == nx_diff, prob.df.Σ is used directly as the process-noise covariance on the differential state (treated as variance-per-step). This is the natural convention when each disturbance input maps 1-to-1 onto one differential state.
• When nw != nx_diff, the wrapper linearizes the continuous-time RHS prob.f_cont w.r.t. the disturbance inputs to obtain Bw = ∂(f_cont)/∂w (size nx_diff × nw), and uses R1_diff = Bw · prob.df.Σ · Bwᵀ as the effective process-noise covariance on the differential state. Bw is a "select" matrix in the typical case where each w_i enters one force equation directly: a row of zeros for kinematic-relation states like D(x) = vx, and a 1 on whichever input drives a given force equation. So prob.df.Σ keeps its variance-per-step interpretation — it just gets routed onto the diff states whose RHS actually contains a disturbance input. This is necessary for index-3 mechanical systems with dim Q ≥ 2, where only the lowest-order (force) equations admit Pantelides-safe noise placement (typically nw = dim Q + 1 < nx_diff).

Initial state should be on the constraint manifold; pass init=true to StateEstimationProblem or provide a consistent x0map.

get_filter(prob::StateEstimationProblem, ::Type{ExtendedKalmanFilter}; p = prob.p, constant_R1=true, kwargs)
get_filter(prob::StateEstimationProblem, ::Type{UnscentedKalmanFilter}; p = prob.p, kwargs)

Instantiate a filter from a state-estimation problem. kwargs are sent to the filter constructor.

The parameter object p is stored in the filter and defaults to prob.p. For the ExtendedKalmanFilter with constant_R1=true, p is also used to compute the constant R1. The keyword p does not affect the initial-state distribution prob.d0, from which the filter takes the numeric type of its internal state. If p contains, e.g., ForwardDiff.Dual numbers, use get_filter(remake(prob; p), ...) instead, see remake.

If constant_R1=true, the dynamics noise covariance matrix R1 is assumed to be constant and is computed as $R_1 = B_w Σ_w B_w^T$ with $B_w = ∂f/∂w$ evaluated at the mean of prob.d0 and u = 0, where $Σ_w$ is the covariance of prob.df. Otherwise, R1 is computed at each time step through repeated linearization w.r.t. the disturbance inputs w.

The UnscentedKalmanFilter propagates the disturbance inputs through the dynamics using augmented sigma points, its R1 is the covariance of prob.df, of size nw × nw.

get_parameters(prob::StateEstimationProblem, syms)

Return a vector containing the current values of the parameters syms in prob.p. The result is, e.g., useful as an initial guess for the parameter vector θ passed to a setter obtained from parameter_setter. The same restrictions on syms as for parameter_setter apply.

setter = parameter_setter(prob::StateEstimationProblem, syms)
p = setter(θ)

Return a callable object that maps a parameter vector θ to a parameter object for prob. The returned parameter object is a copy of prob.p with the entries corresponding to the symbolic parameters syms replaced by the corresponding entries of θ, such that length(θ) == length(syms). The result is intended to be passed to remake or to the p keyword of get_filter, e.g.,

setter = parameter_setter(prob, [cmodel.c, cmodel.k])
filt = get_filter(remake(prob; p = setter(θ)), UnscentedKalmanFilter)

The returned parameter object has the same container type as prob.p. If prob.p is a Tuple (the default), the entries that are not replaced are kept as they are, and the replaced entries take the types of the entries of θ, which makes the setter suitable for use with automatic differentiation, e.g., using ForwardDiff. The construction of the tuple is type stable. If prob.p is an AbstractVector, a copy with the element type promoted to accommodate the entries of θ is returned.

The symbols in syms are given as symbolic variables of the completed model, e.g., complete(model).k. Bound parameters (parameters defined as expressions of other parameters), inputs, disturbance inputs, state variables, observed variables and array-valued parameters are not supported and result in an ArgumentError. Bound parameters are recomputed from the parameters they depend on by the generated code.

See also get_parameters to obtain the current values of the parameters, e.g., as an initial guess for θ.

propagate_distribution(f, kf, dist, args...; kwargs...)

Propagate a probability distribution dist through a nonlinear function f using the covariance-propagation method of filter kf.

Arguments:

• f: A nonlinear function f(x, args...; kwargs...) that takes a vector x and returns a vector.
• kf: A state estimator, such as an ExtendedKalmanFilter or UnscentedKalmanFilter.
• dist: A probability distribution, such as a MvNormal or SimpleMvNormal.
• args: Additional arguments to f.
• kwargs: Additional keyword arguments to f.

remake(prob::StateEstimationProblem; p = prob.p, df = prob.df, dg = prob.dg, d0 = nothing)

Return a copy of prob with the parameter object p, the dynamics-noise distribution df, the measurement-noise distribution dg and the initial-state distribution d0 replaced. No symbolic processing or code generation is performed, the generated functions prob.f, prob.f_cont and prob.g are reused. This makes remake suitable for use inside cost functions for parameter and covariance estimation, see parameter_setter for construction of p.

If d0 === nothing (the default), the initial-state distribution prob.d0 is reused, with its mean and covariance converted to the element type obtained by promoting the element types of prob.d0, df.Σ, dg.Σ and the numeric entries of p. This conversion is required when, e.g., p contains ForwardDiff.Dual numbers, since the filters take the numeric type of their internal state from d0. Static arrays are preserved. An explicitly provided d0 is used as given.

The sample interval Ts cannot be changed by remake since it is part of the discretized dynamics prob.f; construct a new problem to change it.

Pitfalls

• The mean of prob.d0 is not recomputed when p changes. If the problem was constructed with init = true, the initial state solved for during initialization corresponds to the original parameters. Provide a new d0, widen the initial covariance, or construct a new problem if the initial state depends on the parameters.
• Parameters solved for during initialization (parameters with binding missing) are not recomputed by remake, their values in p are used as given.
• Bound parameters, i.e., parameters defined as expressions of other parameters, are not part of p. They are computed by the generated code and are thus recomputed from the entries of p.

ParameterSetter

Callable object returned by parameter_setter. setter(θ) returns a parameter object with the entries corresponding to the symbols setter.syms replaced by the entries of θ.