Pyomo模拟一阶ODE未呈现预期指数衰减问题排查

我正在测试Pyomo的仿真能力，编写了如下脚本模拟一阶ODE（dx/dt + 10*x = 0，时间常数100ms），但运行后输出的曲线未呈现预期的指数衰减，初始点后数值直接归零，调整nfe参数也无改善，请问这是什么原因？

from pyomo.environ import *
from pyomo.dae import *
import matplotlib.pyplot as plt
plt.clf()
# Define a model
m = ConcreteModel()

# Define time and discretize it
m.t = ContinuousSet(bounds=(0, 1))
m.tau = TransformationFactory('dae.finite_difference')
m.tau.apply_to(m, nfe=500, scheme='BACKWARD')

# Define variables
m.x = Var(m.t)
m.dxdt = DerivativeVar(m.x)

# Initial condition
m.ic = Constraint(expr=m.x[0] == 1)
m.x[0].fix(1)  # or any other non-zero value

# DAE Equation (a simple differential equation: dx/dt + 10*x = 0)
#time constant should be 100ms
m.ode = Constraint(m.t, rule=lambda m, t: m.dxdt[t] + 10*m.x[t] == 0)

# Solve the model
solver = SolverFactory('ipopt')
results = solver.solve(m, tee=True)

# Extract the results and plot
time_points = [t for t in m.t]
x_values = [m.x[t]() for t in m.t]
%matplotlib ipympl
plt.plot(time_points, x_values, label='x(t)')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.title('Time Evolution of x(t)')
plt.show()

• 核心原因：DAE离散化变换的顺序错误
  Pyomo的dae.finite_difference变换需要在所有微分变量（包括DerivativeVar）定义完成后再执行。你的代码中先对空模型应用了离散化变换，之后才定义m.x和m.dxdt，导致变换没有处理这些变量，微分约束m.ode未被正确离散化，求解器无法识别有效的微分方程约束，因此出现不符合预期的结果。

• 修复步骤：

  1. 调整代码顺序，先定义变量和导数，再执行离散化变换：

     # 先定义变量和导数
     m.x = Var(m.t)
     m.dxdt = DerivativeVar(m.x)
     # 再应用离散化变换
     m.tau = TransformationFactory('dae.finite_difference')
     m.tau.apply_to(m, nfe=500, scheme='BACKWARD')
     

  2. 移除重复的初始条件约束：m.x[0].fix(1)已经足够固定初始值，无需再添加m.ic = Constraint(expr=m.x[0] == 1)，避免冗余约束干扰求解。

• 修复后的完整代码

from pyomo.environ import *
from pyomo.dae import *
import matplotlib.pyplot as plt
plt.clf()
# Define a model
m = ConcreteModel()

# 先定义时间集合
m.t = ContinuousSet(bounds=(0, 1))

# 先定义变量和导数
m.x = Var(m.t)
m.dxdt = DerivativeVar(m.x)

# 变量定义完成后再执行离散化
m.tau = TransformationFactory('dae.finite_difference')
m.tau.apply_to(m, nfe=500, scheme='BACKWARD')

# 初始条件（仅需一种方式）
m.x[0].fix(1)

# DAE方程约束
m.ode = Constraint(m.t, rule=lambda m, t: m.dxdt[t] + 10*m.x[t] == 0)

# 求解模型
solver = SolverFactory('ipopt')
results = solver.solve(m, tee=True)

# 提取结果并绘图
time_points = [t for t in m.t]
x_values = [m.x[t]() for t in m.t]
%matplotlib ipympl
plt.plot(time_points, x_values, label='x(t)')
plt.xlabel('Time')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.title('Time Evolution of x(t)')
plt.show()

修复后运行代码，即可得到符合预期的指数衰减曲线x(t) = e^(-10t)。

内容的提问来源于stack exchange，提问作者Erik Iverson