Harden your SymPy control algebra: compute the resolvent (sI - A)^-1 exactly to avoid order inflation from floats

Learn why SymPy matrix inverses of (sI - A) inflate model order with floats, and how exact rationals and nsimplify yield clean resolvent and transfer functions.

SymPy and control design often meet at the same place: computing the resolvent (sI − A)⁻¹, the matrix many engineers casually call the state transition matrix in the Laplace domain. A surprising trap appears when symbolic algebra meets floating-point arithmetic. Depending on the inversion method, identical code can yield expressions of very different apparent order, which derails simplification and transfer-function inspection. Let’s walk through what happens and how to make the algebra robust.

Problem setup

Consider a 5×5 system matrix A with physically motivated parameters and the goal of forming ϕ(s) = (sI − A)⁻¹. With floating-point inputs, the same inverse can look second order with one method and fifth order with another.

from sympy import symbols, Matrix, eye
import numpy as np
s = symbols('s')
vol_tank = 20000  # [m^3]
flow_k = 20000/(35*60)
disc_area = np.pi * 20.25**2
coef0 = flow_k/(0.06*disc_area)
coef1 = flow_k/vol_tank
tau = 0.02
A_mat = Matrix([
    [-coef0, 0,      0,      0,      0],
    [0,     -coef1,  1,      1,      1],
    [0,      0,     -1/tau,  0,      0],
    [0,      0,      0,     -1/tau,  0],
    [0,      0,      0,      0,     -1/tau]
])
# Same matrix, different inversion paths
Phi_lu = (s*eye(5) - A_mat).inv(method='LU')
Phi_def = (s*eye(5) - A_mat).inv()

One of these inverses presents a compact, low-order rational structure; the other inflates the order. That discrepancy is not about the underlying system but about arithmetic choices.

Why the same inverse looks different

The key is mixing floats and symbols. When floats enter polynomial expressions, coefficients become inexact. Even if the true numerator and denominator share a large common divisor, tiny floating-point perturbations prevent exact GCD cancellation. The result is a fraction that looks higher order than it really is.

Internally, SymPy can convert floats to rationals to do exact polynomial algebra. In this case, however, the conversion back to floats happens too early for the GCD cancellation to fully simplify the result. That’s why one method seems to “work” while another produces an unexpectedly high-order expression. A particular algorithm like LU can occasionally escape the worst of it, but that is incidental and not a sound strategy to rely on.

A stable way to fix it

The most reliable approach is to make the expressions exact before inversion. Either construct the model with exact rationals and sympy.pi from the start, or coerce the matrix to exact form via nsimplify prior to inversion. Both routes keep polynomial arithmetic exact so that cancellations occur as intended. You can always convert to floats at the end for evaluation or display.

If you already built A from floats, coerce to rationals on the fly:

from sympy import nsimplify, eye
As = s*eye(5) - A_mat
Phi_exactlike = nsimplify(As).inv()

Even better, construct everything as exact quantities and only approximate at the end if needed:

import sympy as sp
from sympy import Rational as Q
s = sp.symbols('s')
vol_tank = 20000
flow_k = Q(20000)/(35*60)
disc_area = sp.pi * Q('20.25')**2
coef0 = flow_k/(Q('0.06')*disc_area)
coef1 = flow_k/vol_tank
tau = Q('0.02')
A_mat = sp.Matrix([
    [-coef0, 0,      0,      0,      0],
    [0,     -coef1,  1,      1,      1],
    [0,      0,     -1/tau,  0,      0],
    [0,      0,      0,     -1/tau,  0],
    [0,      0,      0,      0,     -1/tau]
])
Phi_sym = (s*sp.eye(5) - A_mat).inv().applyfunc(sp.factor)
Phi_num = sp.N(Phi_sym)  # optional numeric view at the end

This pattern avoids spurious degree inflation and yields the clean factorized structure that matches expectations. If you prefer a single-shot correction without changing how A is built, the nsimplify approach above provides that.

Additional practical notes

If your pipeline builds transfer functions, discrepancies can show up in functions that convert state-space to tf. It can be helpful to cross-check with control.dcgain, control.poles, and control.zeros, which in one workflow produced values consistent with manually derived transfer functions even when a conversion step did not. When a matrix has a clean structure, exploiting it may also help; here the first row and column are a simple diagonal block, so focusing on the 4×4 submatrix can reduce algebraic clutter.

There is also active work to improve the underlying behavior: a pull request was opened to address the premature float conversion during inverse simplification at https://github.com/sympy/sympy/pull/28168.

Why this matters for control engineers

Model order, pole-zero locations, and DC gain all depend on exact cancellations in the algebra. If spurious higher-order terms sneak in, you can end up with misleading expressions, poor numerical conditioning, or incorrect intuition about system dynamics. This is especially relevant when auditing ss2tf outputs or matching symbolic derivations to numeric checks.

At the same time, exact arithmetic is not a substitute for real-world variability. Even a perfectly simplified symbolic transfer function won’t capture tolerances and approximations in physical components. The point here is to separate algebraic correctness from physical modeling uncertainty and treat them appropriately.

Conclusion

Keep symbolic control algebra exact whenever possible. Use sympy.pi instead of numpy.pi, prefer Rational to binary floats, or coerce mixed expressions with nsimplify before inverting. Choose inversion methods for robustness, not because they accidentally hide floating-point noise. If needed, verify key properties with dcgain, poles, and zeros as an independent sanity check. Only convert to floats at the end, after symbolic simplification has done its job. That small discipline prevents inflated polynomial degrees, preserves intended cancellations, and keeps your control design calculations faithful and readable.