Устойчивая алгебра в SymPy: резольвента (sI − A)⁻¹ без ложного роста порядка

Почему резольвента (sI − A)⁻¹ в SymPy «раздувается» из‑за плавающей точки, и как это исправить: рациональная арифметика и nsimplify для передаточных функций.

SymPy и проектирование систем управления часто сходятся в одной точке: вычислении резольвенты (sI − A)⁻¹ — матрицы, которую многие инженеры неформально называют матрицей перехода состояния в области Лапласа. Когда символическая алгебра сталкивается с арифметикой плавающей точкой, возникает неожиданная ловушка. В зависимости от выбранного способа обращения одна и та же программа может породить выражения с резко различным видимым порядком, что мешает упрощению и анализу передаточных функций. Разберёмся, что происходит и как сделать алгебру устойчивой.

Постановка задачи

Рассмотрим системную матрицу A размером 5×5 с физически обоснованными параметрами и задачей получить ϕ(s) = (sI − A)⁻¹. При входных данных с плавающей точкой одна и та же обратная по одному методу выглядит как второго порядка, а по другому — как пятого.

from sympy import symbols, Matrix, eye
import numpy as np
s = symbols('s')
vol_tank = 20000  # [м^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]
])
# Одна и та же матрица, разные пути обращения
Phi_lu = (s*eye(5) - A_mat).inv(method='LU')
Phi_def = (s*eye(5) - A_mat).inv()

Один из вариантов обратной даёт компактную рациональную структуру низкого порядка; другой искусственно раздувает порядок. Это расхождение связано не с системой как таковой, а с выбором арифметики.

Почему одна и та же обратная выглядит по‑разному

Суть в смешении чисел с плавающей точкой и символов. Как только флоты попадают в полиномы, коэффициенты становятся неточными. Даже если у истинных числителя и знаменателя есть большой общий делитель, малые погрешности плавающей точки мешают точному сокращению по НОД. В итоге получается дробь, выглядящая как более высокого порядка, чем есть на самом деле.

Внутри SymPy умеет переводить флоты в рациональные числа, чтобы выполнять точную полиномиальную алгебру. Но здесь обратное преобразование к флотам происходит слишком рано, и сокращение по НОД не успевает полностью упростить выражение. Поэтому один способ будто бы «срабатывает», а другой возвращает неожиданно высокопорядковую форму. Отдельные алгоритмы, вроде LU, порой случайно обходят худшие эффекты, но это побочный эффект, на который нельзя опираться как на стратегию.

Надёжный способ исправить

Самый надёжный подход — сделать выражения точными до обращения. Либо сразу соберите модель в рациональной арифметике, используя sympy.pi, либо принудительно приведите матрицу к точному виду с помощью nsimplify перед инвертированием. Оба пути сохраняют точность полиномиальной арифметики, и необходимые сокращения выполняются как задумано. Для вычислений или отображения всегда можно перейти к флотам в самом конце.

Если A уже построена из флот‑коэффициентов, можно на лету привести её к рациональному виду:

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

Ещё лучше — изначально задать все величины точно и лишь при необходимости приблизить результат в конце:

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)  # по желанию — числовой вид в конце

Такая схема предотвращает ложное завышение степени и даёт аккуратную факторизованную структуру, соответствующую ожиданиям. Если хочется обойтись одним шагом, не меняя способ построения A, используйте приём с nsimplify, показанный выше.

Дополнительные практические замечания

Если ваш конвейер формирует передаточные функции, расхождения могут проявляться в функциях преобразования из пространственного представления состояния в tf. Полезно перепроверять результаты с помощью control.dcgain, control.poles и control.zeros: в одном рабочем процессе они давали значения, согласующиеся с вручную выведенными передаточными функциями, даже если шаг преобразования — нет. Когда у матрицы чистая структура, ею тоже стоит воспользоваться; здесь первая строка и первый столбец образуют простой диагональный блок, поэтому фокус на подматрице 4×4 уменьшает алгебраический «шум».

Над улучшением поведения тоже ведётся работа: открыт pull request, устраняющий преждевременное преобразование к флотам при упрощении обратной матрицы: https://github.com/sympy/sympy/pull/28168.

Почему это важно для инженеров по управлению

Порядок модели, расположение полюсов и нулей и DC‑усиление зависят от точных сокращений в алгебре. Если в выражения просачиваются ложные высокие степени, вы рискуете получить вводящую в заблуждение запись, плохую численную обусловленность или неверную интуицию о динамике системы. Это особенно актуально при проверке результатов ss2tf или сопоставлении символических выводов с численными проверками.

При этом точная арифметика не заменяет реальную изменчивость. Даже идеально упрощённая символическая передаточная функция не отражает допуски и аппроксимации физических компонентов. Важно отделять алгебраическую корректность от неопределённости физического моделирования и относиться к ним по‑разному.

Выводы

Старайтесь по возможности сохранять символическую алгебру управления точной. Используйте sympy.pi вместо numpy.pi, предпочитайте Rational двоичным флотам или приводите смешанные выражения к точному виду через nsimplify до обращения. Выбирайте методы инвертирования из соображений устойчивости, а не потому, что они случайно прячут шум плавающей точки. При необходимости проверяйте ключевые свойства через dcgain, poles и zeros как независимую sanity‑проверку. Переходите к флотам только в самом конце, когда символическое упрощение уже выполнено. Эта небольшая дисциплина предотвращает раздувание степеней, сохраняет задуманные сокращения и делает расчёты по управлению корректными и читаемыми.