Linearization

A nonlinear dynamical system with state (differential and algebraic) $x$ and input signals $u$

\[M \dot x = f(x, u)\]

can be linearized using the function linearize to produce a linear statespace system on the form

\[\begin{aligned} \dot x &= Ax + Bu\\ y &= Cx + Du \end{aligned}\]

The linearize function expects the user to specify the inputs $u$ and the outputs $y$ using the syntax shown in the example below. The system model is not supposed to be simplified before calling linearize:

Example

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
@variables x(t)=0 y(t) u(t) r(t)=0
@parameters kp = 1

eqs = [u ~ kp * (r - y) # P controller
       D(x) ~ -x + u    # First-order plant
       y ~ x]           # Output equation

@named sys = ODESystem(eqs, t) # Do not call @mtkbuild when linearizing
matrices, simplified_sys = linearize(sys, [r], [y]) # Linearize from r to y
matrices

(A = [-2.0;;], B = [1.0;;], C = [1.0;;], D = [0.0;;])

The named tuple matrices contains the matrices of the linear statespace representation, while simplified_sys is an ODESystem that, among other things, indicates the unknown variable order in the linear system through

using ModelingToolkit: inputs, outputs
[unknowns(simplified_sys); inputs(simplified_sys); outputs(simplified_sys)]

3-element Vector{Any}:
 x(t)
 r(t)
 y(t)

Operating point

The operating point to linearize around can be specified with the keyword argument op like this: op = Dict(x => 1, r => 2). The operating point may include specification of unknown variables, input variables and parameters. For variables that are not specified in op, the default value specified in the model will be used if available, if no value is specified, an error is thrown.

Batch linearization and algebraic variables

If linearization is to be performed around multiple operating points, the simplification of the system has to be carried out a single time only. To facilitate this, the lower-level function ModelingToolkit.linearization_function is available. This function further allows you to obtain separate Jacobians for the differential and algebraic parts of the model. For ODE models without algebraic equations, the statespace representation above is available from the output of linearization_function as A, B, C, D = f_x, f_u, h_x, h_u.

All variables that will be fixed by an operating point must be provided in the operating point to linearization_function. For example, if the operating points fix the value of x, y and z then an operating point with constant values for these variables (e.g. Dict(x => 1.0, y => 1.0, z => 1.0)) must be provided. The constant values themselves do not matter and can be changed by subsequent operating points.

One approach to batch linearization would be to call linearize in a loop, providing a different operating point each time. For example:

using ModelingToolkitStandardLibrary
using ModelingToolkitStandardLibrary.Blocks

@parameters k=10 k3=2 c=1
@variables x(t)=0 [bounds = (-0.5, 1.5)]
@variables v(t) = 0

@named y = Blocks.RealOutput()
@named u = Blocks.RealInput()

eqs = [D(x) ~ v
       D(v) ~ -k * x - k3 * x^3 - c * v + 10u.u
       y.u ~ x]

@named duffing = ODESystem(eqs, t, systems = [y, u], defaults = [u.u => 0])

# pass a constant value for `x`, since it is the variable we will change in operating points
linfun, simplified_sys = linearization_function(duffing, [u.u], [y.u]; op = Dict(x => NaN));

println(linearize(simplified_sys, linfun; op = Dict(x => 1.0)))
println(linearize(simplified_sys, linfun; op = Dict(x => 0.0)))

@time linearize(simplified_sys, linfun; op = Dict(x => 0.0))

(A = [-1.0 -16.0; 1.0 0.0], B = [10.0; 0.0;;], C = [0.0 1.0], D = [0.0;;])
(A = [-1.0 -10.0; 1.0 0.0], B = [10.0; 0.0;;], C = [0.0 1.0], D = [0.0;;])
  0.000427 seconds (1.14 k allocations: 49.156 KiB)

However, this route is still expensive since it has to repeatedly process the symbolic map provided to op. linearize is simply a wrapper for creating and solving a ModelingToolkit.LinearizationProblem. This object is symbolically indexable, and can thus integrate with SymbolicIndexingInterface.jl for fast updates.

using SymbolicIndexingInterface

# The second argument is the value of the independent variable `t`.
linprob = LinearizationProblem(linfun, 1.0)
# It can be mutated
linprob.t = 0.0
# create a setter function to update `x` efficiently
setter! = setu(linprob, x)

function fast_linearize!(problem, setter!, value)
    setter!(problem, value)
    solve(problem)
end

println(fast_linearize!(linprob, setter!, 1.0))
println(fast_linearize!(linprob, setter!, 0.0))

@time fast_linearize!(linprob, setter!, 1.0)

(A = [-1.0 -10.0; 1.0 0.0], B = [10.0; 0.0;;], C = [0.0 1.0], D = [0.0;;])
(A = [-1.0 -10.0; 1.0 0.0], B = [10.0; 0.0;;], C = [0.0 1.0], D = [0.0;;])
  0.000113 seconds (139 allocations: 17.016 KiB)

Note that linprob above can be interacted with similar to a normal ODEProblem.

julia> prob[x]ERROR: UndefVarError: `prob` not defined
julia> prob[x] = 1.5ERROR: UndefVarError: `prob` not defined
julia> prob[x]ERROR: UndefVarError: `prob` not defined

Symbolic linearization

The function ModelingToolkit.linearize_symbolic works similar to ModelingToolkit.linearize but returns symbolic rather than numeric Jacobians. Symbolic linearization have several limitations and no all systems that can be linearized numerically can be linearized symbolically.

Input derivatives

Physical systems are always proper, i.e., they do not differentiate causal inputs. However, ModelingToolkit allows you to model non-proper systems, such as inverse models, and may sometimes fail to find a realization of a proper system on proper form. In these situations, linearize may throw an error mentioning

Input derivatives appeared in expressions (-g_z\g_u != 0)

This means that to simulate this system, some order of derivatives of the input is required. To allow linearize to proceed in this situation, one may pass the keyword argument allow_input_derivatives = true, in which case the resulting model will have twice as many inputs, $2n_u$, where the last $n_u$ inputs correspond to $\dot u$.

If the modeled system is actually proper (but MTK failed to find a proper realization), further numerical simplification can be applied to the resulting statespace system to obtain a proper form. Such simplification is currently available in the package ControlSystemsMTK.

Tools for linear analysis

ModelingToolkit contains a set of tools for more advanced linear analysis. These can be used to make it easier to work with and analyze causal models, such as control and signal-processing systems.

Also see ControlSystemsMTK.jl for an interface to ControlSystems.jl that contains tools for linear analysis and frequency-domain analysis.

Docstrings

• ModelingToolkit.LinearizationProblem
• ModelingToolkit.linearization_function
• ModelingToolkit.linearize
• ModelingToolkit.linearize_symbolic

(; A, B, C, D), simplified_sys = linearize(sys, inputs, outputs;    t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
(; A, B, C, D)                 = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)

  r ┌─────┐       ┌─────┐     ┌─────┐
───►│     ├──────►│     │  u  │     │
    │  F  │       │  C  ├────►│  P  │ y
    └─────┘     ┌►│     │     │     ├─┬─►
                │ └─────┘     └─────┘ │
                │                     │
                └─────────────────────┘

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
function plant(; name)
    @variables x(t) = 1
    @variables u(t)=0 y(t)=0
    eqs = [D(x) ~ -x + u
           y ~ x]
    ODESystem(eqs, t; name = name)
end

function ref_filt(; name)
    @variables x(t)=0 y(t)=0
    @variables u(t)=0 [input = true]
    eqs = [D(x) ~ -2 * x + u
           y ~ x]
    ODESystem(eqs, t, name = name)
end

function controller(kp; name)
    @variables y(t)=0 r(t)=0 u(t)=0
    @parameters kp = kp
    eqs = [
        u ~ kp * (r - y),
    ]
    ODESystem(eqs, t; name = name)
end

@named f = ref_filt()
@named c = controller(1)
@named p = plant()

connections = [f.y ~ c.r # filtered reference to controller reference
               c.u ~ p.u # controller output to plant input
               p.y ~ c.y]

@named cl = ODESystem(connections, t, systems = [f, c, p])

lsys0, ssys = linearize(cl, [f.u], [p.x])
desired_order = [f.x, p.x]
lsys = ModelingToolkit.reorder_unknowns(lsys0, unknowns(ssys), desired_order)

@assert lsys.A == [-2 0; 1 -2]
@assert lsys.B == [1; 0;;]
@assert lsys.C == [0 1]
@assert lsys.D[] == 0

## Symbolic linearization
lsys_sym, _ = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x])

@assert substitute(lsys_sym.A, ModelingToolkit.defaults(cl)) == lsys.A

(; A, B, C, D), simplified_sys = linearize_symbolic(sys::AbstractSystem, inputs, outputs; simplify = false, allow_input_derivatives = false, kwargs...)

lin_fun, simplified_sys = linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, initialize = true, initialization_solver_alg = TrustRegion(), kwargs...)

mutable struct LinearizationProblem{F<:ModelingToolkit.LinearizationFunction, T}