如何用SymPy消去中间变量并化简键合图方程组为ODE

我有一个键合图程序，它生成了一组方程关系，但最终应化简为仅含非中间变量的几个常微分方程（ODE）。通过替换中间关系，给定的方程组应化简为指定的ODE形式。用sym.solve硬编码解决可行，但无法泛化，请问如何用SymPy实现这类替换操作？

原方程组 #

Eq(SE_01, e_01)
Eq(f_02, p_02/I_02)
Eq(e_05, q_05/C_05)
Eq(e_03, R_03*f_03)
Eq(f_06, e_06/R_06)
Eq(e_02, pdot_02)
Eq(f_05, qdot_05)
Eq(f_01, f_02)
Eq(f_02, f_03)
Eq(f_02, f_04)
Eq(e_01, e_02 + e_03 + e_04)
Eq(e_04, e_05)
Eq(e_05, e_06)
Eq(f_04, f_05 + f_06)

目标化简结果 #

pdot_02 = -R_03*p_02/I_02 - q_05/C_05 + SE_01
qdot_05 = p_02/I_02 - q_05/(C_05*R_06)

已尝试的硬编码代码 #

solution_explicit = sym.solve(equations, [pdot_02, qdot_05, e_01, e_02, e_03, e_04, e_05, e_06, f_01, f_02, f_03, f_04, f_05, f_06], dict=True)

print("\nExplicit solution for pdot_02 and qdot_05:")
if solution_explicit:
    sol = solution_explicit[0]
    pdot_02_sol = sol.get(pdot_02)
    qdot_05_sol = sol.get(qdot_05)

    if pdot_02_sol is not None:
        sym.pprint(sym.simplify(sym.Eq(pdot_02, pdot_02_sol)))
    else:
        print("pdot_02 not found in this solution set.")

    if qdot_05_sol is not None:
        sym.pprint(sym.simplify(sym.Eq(qdot_05, qdot_05_sol)))

泛化解法实现 #

要实现泛化的替换化简，无需硬编码所有中间变量，可按以下步骤操作：

核心思路 #

1. 区分变量类型：明确目标变量（状态导数，如pdot_02、qdot_05）、保留变量（状态、输入、参数，如SE_01、p_02）、中间变量（待消去的e_*、f_*类变量）。
2. 自动提取中间变量：通过集合运算从方程组的所有符号中筛选出待消去的中间变量，无需手动枚举。
3. 分步消元替换：先解出中间变量的显式表达式，再代入含目标导数的方程，最终化简得到目标ODE。

泛化代码示例 #

import sympy as sym

# 定义所有符号
SE_01, e_01, f_02, p_02, I_02, e_05, q_05, C_05, e_03, R_03, f_03, f_06, e_06, R_06, e_02, pdot_02, f_05, qdot_05, f_01, f_04 = sym.symbols('SE_01 e_01 f_02 p_02 I_02 e_05 q_05 C_05 e_03 R_03 f_03 f_06 e_06 R_06 e_02 pdot_02 f_05 qdot_05 f_01 f_04')

# 构建方程组
equations = [
    sym.Eq(SE_01, e_01),
    sym.Eq(f_02, p_02/I_02),
    sym.Eq(e_05, q_05/C_05),
    sym.Eq(e_03, R_03*f_03),
    sym.Eq(f_06, e_06/R_06),
    sym.Eq(e_02, pdot_02),
    sym.Eq(f_05, qdot_05),
    sym.Eq(f_01, f_02),
    sym.Eq(f_02, f_03),
    sym.Eq(f_02, f_04),
    sym.Eq(e_01, e_02 + e_03 + e_04),
    sym.Eq(e_04, e_05),
    sym.Eq(e_05, e_06),
    sym.Eq(f_04, f_05 + f_06)
]

# 定义目标变量（ODE的导数）
target_derivatives = [pdot_02, qdot_05]
# 定义需要保留的非导数变量（状态、输入、参数）
retained_vars = [SE_01, p_02, q_05, I_02, C_05, R_03, R_06]

# 自动提取中间变量：所有不在目标导数和保留变量里的符号
all_symbols = set().union(*[eq.free_symbols for eq in equations])
intermediate_vars = [sym for sym in all_symbols if sym not in target_derivatives and sym not in retained_vars]

# 第一步：解出所有中间变量的表达式
intermediate_solutions = sym.solve(equations, intermediate_vars, dict=True)[0]

# 第二步：将中间变量代入包含目标导数的方程，并整理出目标导数的表达式
target_eq = [eq for eq in equations if any(deriv in eq.free_symbols for deriv in target_derivatives)]
final_odes = []

for eq in target_eq:
    substituted_eq = eq.subs(intermediate_solutions)
    for deriv in target_derivatives:
        if deriv in substituted_eq.free_symbols:
            sol = sym.solve(substituted_eq, deriv)[0]
            final_odes.append(sym.Eq(deriv, sym.simplify(sol)))

# 输出结果
for ode in final_odes:
    sym.pprint(ode)

代码说明 #

• 无需硬编码中间变量，通过集合运算自动适配不同的键合图方程组；
• 分步替换逻辑清晰，先解中间变量再代入目标方程，可扩展至线性/非线性方程组；
• 最终输出的ODE会自动化简为仅含保留变量的形式。

内容的提问来源于stack exchange，提问作者Austin Flynt