GitHub

ControlSystemsMTK.jl

ControlSystemsMTK provides an interface between ControlSystems.jl and ModelingToolkit.jl.

See the videos below for examples of using ControlSystems and ModelingToolkit together.

Installation

pkg> add ControlSystemsMTK

From ControlSystems to ModelingToolkit

Simply calling ODESystem(sys) converts a StateSpace object from ControlSystems into the corresponding ModelingToolkitStandardLibrary.Blocks.StateSpace. If sys is a named statespace object, the names will be retained in the ODESystem.

Example:

julia> using ControlSystemsMTK, ControlSystemsBase, ModelingToolkit, RobustAndOptimalControl

julia> P0 = tf(1.0, [1, 1])  |> ss
StateSpace{Continuous, Float64}
A = 
 -1.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

julia> @named P = ODESystem(P0)
Model P with 2 equations
States (3):
  x[1](t) [defaults to 0.0]
  input₊u(t) [defaults to 0.0]
  output₊u(t) [defaults to 0.0]
Parameters (0):

julia> equations(P)
2-element Vector{Equation}:
 Differential(t)(x[1](t)) ~ input₊u(t) - x[1](t)
 output₊u(t) ~ x[1](t)

From ModelingToolkit to ControlSystems

An ODESystem can be converted to a named statespace object from RobustAndOptimalControl.jl by calling named_ss

named_ss(ode_sys, inputs, outputs; op)

this performs a linearization of ode_sys around the operating point op (defaults to the default values of all variables in ode_sys).

Example:

Using P from above:

julia> @unpack input, output = P;

julia> P02_named = named_ss(P, [input.u], [output.u])
NamedStateSpace{Continuous, Float64}
A = 
 -1.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model
With state  names: x[1](t)
     input  names: input₊u(t)
     output names: output₊u(t)

julia> using Plots;

julia> bodeplot(P02_named)

plot

julia> ss(P02_named) # Convert to a statespace system without names
StateSpace{Continuous, Float64}
A = 
 -1.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

ModelingToolkit tends to give weird names to inputs and outputs etc., to access variables easily, named_ss implements prefix matching, so that you can access the mapping from input₊u(t) to output₊u(t) by

P02_named[:out, :in]

To learn more about linearization of ModelingToolkit models, see the video below

Symbolic linearization and code generation

ModelingToolkit has facilities for symbolic linearization that can be used on sufficiently simple systems. The function linearize_symbolic behaves similarly to linearize but returns symbolic matrices A,B,C,D rather than numeric. A StateSpace system with such symbolic coefficients can be used to generate a function that takes parameter values and outputs a statically sized statespace system with numeric matrices. An example follows

System model

We start by building a system mode, we'll use a model of two masses connected by a flexible transmission

using ControlSystemsMTK, ControlSystemsBase
using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
using ModelingToolkitStandardLibrary.Mechanical.Rotational
using ModelingToolkitStandardLibrary.Blocks: Sine
using ModelingToolkit: connect
import ModelingToolkitStandardLibrary.Blocks
t = Blocks.t

# Parameters
m1 = 1
m2 = 1
k = 1000 # Spring stiffness
c = 10   # Damping coefficient
@named inertia1 = Inertia(; J = m1)
@named inertia2 = Inertia(; J = m2)
@named spring = Spring(; c = k)
@named damper = Damper(; d = c)
@named torque = Torque(use_support=false)

function SystemModel(u=nothing; name=:model)
    @named sens = Rotational.AngleSensor()
    eqs = [
        connect(torque.flange, inertia1.flange_a)
        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
        connect(inertia2.flange_b, sens.flange)
    ]
    if u !== nothing
        push!(eqs, connect(u.output, :u, torque.tau))
        return @named model = ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper, u])
    end
    ODESystem(eqs, t; systems = [sens, torque, inertia1, inertia2, spring, damper], name)
end

model = SystemModel() |> complete

\[ \begin{equation} \left[ \begin{array}{c} \mathrm{connect}\left( torque_{+}flange, inertia1_{+}flange_{a} \right) \\ \mathrm{connect}\left( inertia1_{+}flange_{b}, spring_{+}flange_{a}, damper_{+}flange_{a} \right) \\ \mathrm{connect}\left( inertia2_{+}flange_{a}, spring_{+}flange_{b}, damper_{+}flange_{b} \right) \\ \mathrm{connect}\left( inertia2_{+}flange_{b}, sens_{+}flange \right) \\ sens_{+}phi_{+}u\left( t \right) = sens_{+}flange_{+}phi\left( t \right) \\ sens_{+}flange_{+}tau\left( t \right) = 0 \\ torque_{+}phi_{support}\left( t \right) = 0 \\ torque_{+}flange_{+}tau\left( t \right) = - torque_{+}tau_{+}u\left( t \right) \\ inertia1_{+}phi\left( t \right) = inertia1_{+}flange_{a_{+}phi}\left( t \right) \\ inertia1_{+}phi\left( t \right) = inertia1_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} inertia1_{+}phi\left( t \right)}{\mathrm{d}t} = inertia1_{+}w\left( t \right) \\ \frac{\mathrm{d} inertia1_{+}w\left( t \right)}{\mathrm{d}t} = inertia1_{+}a\left( t \right) \\ inertia1_{+}J inertia1_{+}a\left( t \right) = inertia1_{+}flange_{a_{+}tau}\left( t \right) + inertia1_{+}flange_{b_{+}tau}\left( t \right) \\ inertia2_{+}phi\left( t \right) = inertia2_{+}flange_{a_{+}phi}\left( t \right) \\ inertia2_{+}phi\left( t \right) = inertia2_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} inertia2_{+}phi\left( t \right)}{\mathrm{d}t} = inertia2_{+}w\left( t \right) \\ \frac{\mathrm{d} inertia2_{+}w\left( t \right)}{\mathrm{d}t} = inertia2_{+}a\left( t \right) \\ inertia2_{+}J inertia2_{+}a\left( t \right) = inertia2_{+}flange_{b_{+}tau}\left( t \right) + inertia2_{+}flange_{a_{+}tau}\left( t \right) \\ spring_{+}phi_{rel}\left( t \right) = spring_{+}flange_{b_{+}phi}\left( t \right) - spring_{+}flange_{a_{+}phi}\left( t \right) \\ spring_{+}flange_{b_{+}tau}\left( t \right) = spring_{+}tau\left( t \right) \\ spring_{+}flange_{a_{+}tau}\left( t \right) = - spring_{+}tau\left( t \right) \\ spring_{+}tau\left( t \right) = spring_{+}c \left( - spring_{+}phi_{rel0} + spring_{+}phi_{rel}\left( t \right) \right) \\ damper_{+}phi_{rel}\left( t \right) = - damper_{+}flange_{a_{+}phi}\left( t \right) + damper_{+}flange_{b_{+}phi}\left( t \right) \\ \frac{\mathrm{d} damper_{+}phi_{rel}\left( t \right)}{\mathrm{d}t} = damper_{+}w_{rel}\left( t \right) \\ \frac{\mathrm{d} damper_{+}w_{rel}\left( t \right)}{\mathrm{d}t} = damper_{+}a_{rel}\left( t \right) \\ damper_{+}flange_{b_{+}tau}\left( t \right) = damper_{+}tau\left( t \right) \\ damper_{+}flange_{a_{+}tau}\left( t \right) = - damper_{+}tau\left( t \right) \\ damper_{+}tau\left( t \right) = damper_{+}d damper_{+}w_{rel}\left( t \right) \\ \end{array} \right] \end{equation} \]

Numeric linearization

We can linearize this model numerically using named_ss, this produces a NamedStateSpace{Continuous, Float64}

lsys = named_ss(model, [model.torque.tau.u], [model.inertia1.phi, model.inertia2.phi])
model: NamedStateSpace{Continuous, Float64}
A = 
     0.0      0.0    1.0    0.0
     0.0      0.0    0.0    1.0
 -1000.0   1000.0  -10.0   10.0
  1000.0  -1000.0   10.0  -10.0
B = 
  0.0
  0.0
  1.0
 -0.0
C = 
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
D = 
 0.0
 0.0

Continuous-time state-space model
With state  names: inertia1₊phi(t) inertia2₊phi(t) inertia1₊w(t) inertia2₊w(t)
     input  names: torque₊tau₊u(t)
     output names: inertia1₊phi(t) inertia2₊phi(t)

Symbolic linearization

If we instead call linearize_symbolic and pass the jacobians into ss, we get a StateSpace{Continuous, Num}

mats, simplified_sys = ModelingToolkit.linearize_symbolic(model, [model.torque.tau.u], [model.inertia1.phi, model.inertia2.phi])
symbolic_sys = ss(mats.A, mats.B, mats.C, mats.D)
StateSpace{Continuous, Num}
A = 
                0                         0               1                         0
                0                         0               0                         1
 (-spring₊c) / inertia1₊J     spring₊c / inertia1₊J  (-damper₊d) / inertia1₊J    damper₊d / inertia1₊J
    spring₊c / inertia2₊J  (-spring₊c) / inertia2₊J    damper₊d / inertia2₊J   (-damper₊d) / inertia2₊J
B = 
      0
      0
 1 / inertia1₊J
      0
C = 
 1  0  0  0
 0  1  0  0
D = 
 0
 0

Continuous-time state-space model

Code generation

That's pretty cool, but even nicer is to generate some code for this symbolic system. Below, we use build_function to generate a function that takes a numeric vector x representing the values of the state, and a vector of parameters, and returns a StaticStateSpace{Continuous, Float64}. We pass the keyword argument force_SA=true to build_function to get an allocation-free function. 

defs = ModelingToolkit.defaults(simplified_sys)
x, pars = ModelingToolkit.get_u0_p(simplified_sys, defs, defs) # Extract the default state and parameter values

fun = Symbolics.build_function(symbolic_sys, states(simplified_sys), ModelingToolkit.parameters(simplified_sys);
    expression = Val{false}, # Generate a compiled function rather than a Julia expression
    force_SA   = true,       # Use static arrays instead of regular arrays, for higher performance 🚀
)

static_lsys = fun(x, pars) # We now have a function that takes state and parameters and returns the linearized system!
StaticStateSpace{Continuous, 2, 1, 4, StaticArraysCore.SMatrix{4, 4, Float64, 16}, StaticArraysCore.SMatrix{4, 1, Float64, 4}, StaticArraysCore.SMatrix{2, 4, Int64, 8}, StaticArraysCore.SMatrix{2, 1, Int64, 2}}
A = 
     0.0      0.0    1.0    0.0
     0.0      0.0    0.0    1.0
 -1000.0   1000.0  -10.0   10.0
  1000.0  -1000.0   10.0  -10.0
B = 
 0.0
 0.0
 1.0
 0.0
C = 
 1  0  0  0
 0  1  0  0
D = 
 0
 0

Continuous-time state-space model

It's pretty fast

using BenchmarkTools
@btime $fun($x, $pars)
8.484 ns (0 allocations: 0 bytes)

faster than multiplying two integers in python.

C-code generation

If you prefer to get C-code for deployment onto an embedded target, the types from Symbolics can be converted to SymPy symbols using symbolics_to_sympy. After this, the function SymbolicControlSystems.ccode is called to generate the C-code. The symbols that are present in the system will be considered input arguments in the generated code.

using SymbolicControlSystems
Asp = SymbolicControlSystems.Sym.(Symbolics.symbolics_to_sympy.(mats.A))
Bsp = SymbolicControlSystems.Sym.(Symbolics.symbolics_to_sympy.(mats.B))
Csp = SymbolicControlSystems.Sym.(Symbolics.symbolics_to_sympy.(mats.C))
Dsp = SymbolicControlSystems.Sym.(Symbolics.symbolics_to_sympy.(mats.D))
sys_sp = ss(Asp, Bsp, Csp, Dsp)
sympars_sp = Symbolics.symbolics_to_sympy.(sympars)

discrete_sys_sp = c2d(sys_sp, 0.01, :tustin) # We can only generate C-code for discrete systems

code = SymbolicControlSystems.ccode(discrete_sys_sp; function_name="perfectly_grilled_hotdogs")

This produces the following code, 

#include <stdio.h>
#include <math.h>

void perfectly_grilled_hotdogs(double *y, double u, double damper_d, double inertia1_J, double inertia2_J, double spring_c) {
    static double x[4] = {0};  // Current state
    double xp[4] = {0};        // Next state
    int i;

    // Common sub expressions. These are all called xi, but are unrelated to the state x
    double x0 = inertia2_J*spring_c;
    double x1 = 2.0*x0;
    double x2 = 200.0*damper_d;
    double x3 = inertia2_J*x2;
    double x4 = x0 + x3;
    double x5 = 40000.0*inertia2_J;
    double x6 = inertia1_J*x5;
    double x7 = inertia1_J*spring_c;
    double x8 = inertia1_J*x2;
    double x9 = x7 + x8;
    double x10 = x6 + x9;
    double x11 = 1.0/(x10 + x4);
    double x12 = x11*x[2];
    double x13 = damper_d*inertia1_J;
    double x14 = inertia1_J*inertia2_J;
    double x15 = 1.0/(damper_d*x5 + 200.0*x0 + 40000.0*x13 + 8000000.0*x14 + 200.0*x7);
    double x16 = spring_c + x2;
    double x17 = 0.01*u;
    double x18 = x17*(x16 + x5);
    double x19 = 20000.0*damper_d;
    double x20 = 1.0/(inertia1_J*x19 + inertia2_J*x19 + 100.0*x0 + 4000000.0*x14 + 100.0*x7);
    double x21 = x20*x[3];
    double x22 = x20*x[1];
    double x23 = -x0;
    double x24 = x10 + x3;
    double x25 = x11*x[0];
    double x26 = x0*x12;
    double x27 = x11*x[3];
    double x28 = x11*x[1];
    double x29 = 2.0*x7;
    double x30 = x16*x17;
    double x31 = x4 + x6;
    double x32 = -x7;
    double x33 = x31 + x8;
    double x34 = x25*x7;
    double x35 = x15*x[3];
    double x36 = x15*x[1];
    double x37 = 0.005*u*x15/inertia1_J;

    // Advance the state xp = Ax + Bu
    xp[0] = (x1*x12 + x10*x22 + x15*x18 + x21*x4 + x25*(x23 + x24));
    xp[1] = (-400.0*x0*x25 + x11*x18 + 400.0*x26 + x27*(400.0*damper_d*inertia2_J + x1) + x28*(x10 + x23 - x3));
    xp[2] = (x12*(x32 + x33) + x15*x30 + x21*x31 + x22*x9 + x25*x29);
    xp[3] = (x11*x30 - 400.0*x12*x7 + x27*(x31 + x32 - x8) + x28*(400.0*x13 + x29) + 400.0*x34);

    // Accumulate the output y = C*x + D*u
    y[0] = (x10*x36 + x10*x37 + x24*x25 + x26 + x35*x4);
    y[1] = (x12*x33 + x31*x35 + x34 + x36*x9 + x37*x9);

    // Make the predicted state the current state
    for (i=0; i < 4; ++i) {
        x[i] = xp[i];
    }

}

Additional resources

Internals: Transformation of non-proper models to proper statespace form

For some models, ModelingToolkit will fail to produce a proper statespace model (a non-proper model is differentiating the inputs, i.e., it has a numerator degree higher than the denominator degree if represented as a transfer function) when calling linearize. For such models, given on the form $\dot x = Ax + Bu + \bar B \dot u$ we create the following augmented descriptor model

\[\begin{aligned} sX &= Ax + BU + s\bar B U \\ [X_u &= U]\\ s(X - \bar B X_u) &= AX + BU \\ s \begin{bmatrix}I & -\bar B \\ 0 & 0 \end{bmatrix} &= \begin{bmatrix} A & 0 \\ 0 & -I\end{bmatrix} \begin{bmatrix}X \\ X_u \end{bmatrix} + \begin{bmatrix} B \\ I_u\end{bmatrix} U \\ sE &= A_e x_e + B_e u \end{aligned}\]

where $X_u$ is a new algebraic state variable and $I_u$ is a selector matrix that picks out the differentiated inputs appearing in $\dot u$ (if all inputs appear, $I_u = I$).

This model may be converted to a proper statespace model (if the system is indeed proper) using DescriptorSystems.dss2ss. All of this is handled automatically by named_ss.