Symbolic Linearization with SymPy

The SymPy backend preserves the full algebraic structure of a Modelica model. With it you can derive linearizations by hand — analytically — instead of finite-differencing a black-box simulator. The Jacobians come out as exact rational expressions in the model's physical parameters.

The notebook imports the rumoca-generated SymPy module, builds the symbolic right-hand side, differentiates to get A and B at an operating point, and inspects controllability from the symbolic matrices.

Symbolic linearization of the PVTOL with SymPy

Take the Rumoca-emitted SymPy model of the PVTOL, build its symbolic right-hand side, and linearize it at hover. The result is LaTeX-ready A, B matrices small enough to read by eye, plus a symbolic controllability rank.

Export step — Rumoca compiles the Modelica source into the SymPy model this notebook imports (run once, already done for you):

rumoca compile model/pvtol.mo -m PVTOL --target sympy -o _generated/sympy

_generated/sympy/PVTOL_sympy.py is machine-generated and never hand-edited.

In [1]:

import sys
from pathlib import Path
import numpy as np
import sympy as sp

ROOT = Path.cwd()
sys.path.insert(0, str(ROOT / "_generated" / "sympy"))

import PVTOL_sympy as sym_plant

STATE_NAMES = ["x", "z", "theta", "x_dot", "z_dot", "theta_dot"]
INPUT_NAMES = ["T", "M"]

Build the symbolic right-hand side from `Model.explicit_solution`

In [2]:

mdl = sym_plant.Model()
mdl.solve_explicit()
sol = mdl.explicit_solution

x_dot, z_dot, theta_dot = sp.symbols("x_dot z_dot theta_dot")
m_sym, J_sym = sp.symbols("m J")

der_x_dot     = sol[m_sym * sp.Symbol("x_dot_dot")]     / m_sym
der_z_dot     = sol[m_sym * sp.Symbol("z_dot_dot")]     / m_sym
der_theta_dot = sol[J_sym * sp.Symbol("theta_dot_dot")] / J_sym

f_sym = sp.Matrix([x_dot, z_dot, theta_dot, der_x_dot, der_z_dot, der_theta_dot])
for n, expr in zip(STATE_NAMES, f_sym):
    print(f"  d{n}/dt = {expr}")

  dx/dt = x_dot
  dz/dt = z_dot
  dtheta/dt = theta_dot
  dx_dot/dt = -T*sin(theta)/m
  dz_dot/dt = (T*cos(theta) - g*m)/m
  dtheta_dot/dt = M/J

Symbolic Jacobians A and B

In [3]:

x_syms = sp.Matrix([sp.Symbol(n) for n in STATE_NAMES])
u_syms = sp.Matrix([sp.Symbol(n) for n in INPUT_NAMES])

A_sym = sp.simplify(f_sym.jacobian(x_syms))
B_sym = sp.simplify(f_sym.jacobian(u_syms))

from IPython.display import display, Math
display(Math("A(x, u, p) = " + sp.latex(A_sym)))
display(Math("B(x, u, p) = " + sp.latex(B_sym)))

A (x, u, p) = 000000000000 000 - \frac{T c o s ( θ )}{m} - \frac{T s i n ( θ )}{m} 0 100000010000001000

B (x, u, p) = 000 - \frac{s i n ( θ )}{m} \frac{c o s ( θ )}{m} 0 00000 \frac{1}{J}

Substitute hover trim

x = 0, θ = 0, T = m g, M = 0.

In [4]:

hover_subs = {
    sp.Symbol("x"): 0, sp.Symbol("z"): 0, sp.Symbol("theta"): 0,
    sp.Symbol("x_dot"): 0, sp.Symbol("z_dot"): 0, sp.Symbol("theta_dot"): 0,
    sp.Symbol("T"): m_sym * sp.Symbol("g"),
    sp.Symbol("M"): 0,
}
A_hover = sp.simplify(A_sym.subs(hover_subs))
B_hover = sp.simplify(B_sym.subs(hover_subs))

display(Math("A_{\\rm hover} = " + sp.latex(A_hover)))
display(Math("B_{\\rm hover} = " + sp.latex(B_hover)))

A_{hover} = 000000000000 000 - g 00 100000010000001000

B_{hover} = 0000 \frac{1}{m} 0 00000 \frac{1}{J}

Symbolic controllability

Matrix.rank() computes the rank of the controllability matrix $,[B,\ AB,\ A^2 B,\ \dots,\ A^{5} B],$ symbolically and parameter-independently — no numerical pipeline can make this statement.

In [5]:

n = A_hover.shape[0]
ctrl_cols = [B_hover]
Ak = sp.eye(n)
for _ in range(1, n):
    Ak = Ak * A_hover
    ctrl_cols.append(Ak * B_hover)
C_sym = sp.Matrix.hstack(*ctrl_cols)
print(f"Controllability rank = {C_sym.rank()} / {n}")

Controllability rank = 6 / 6

← Parameter Identification Nonlinear MPC →