Pyomo.DAE中添加含解变量的PDE项时循环引用错误的解决方法

问题描述 #

我尝试用Pyomo.DAE求解某PDE方程，当仅保留方程前三项（$Re\frac{\partial C}{\partial t} + v\frac{\partial C}{\partial x} - D\frac{\partial^2 C}{\partial x^2} = 0$）时，代码可正常运行；但添加包含解变量C的第四项$-Re·v_{var}·C$后，运行SolverFactory('ipopt').solve(m, tee=True)时返回循环引用错误。已知该方程存在解析解，需解决如何在PDE中正确整合解变量C的问题。

可正常运行的前三项代码 #

m = ConcreteModel()

# 定义参数
m.D= Param(initialize = D)
m.v = Param(initialize = v)
m.Re = Param(initialize = Re)
m.v_var = Param(initialize = v_var)
m.u = Param(initialize = u)

# 定义空间/时间网格
m.x = ContinuousSet(bounds=(0,L))
m.t = ContinuousSet(bounds=(0,t_end))

# 定义解空间变量
m.C = Var(m.t, m.x)

# 定义偏导数
m.dCdt   = DerivativeVar(m.C, wrt=m.t)
m.dCdx   = DerivativeVar(m.C, wrt=m.x)
m.d2Cdx2 = DerivativeVar(m.C, wrt=(m.x, m.x))

@m.Constraint(m.t, m.x)
def pde(m, t, x):
    if t == 0:
        return Constraint.Skip
    if x == 0 or x == L:
        return Constraint.Skip
    return m.Re * m.dCdt[t,x] + m.v * m.dCdx[t,x] - m.D * m.d2Cdx2[t,x] == 0

添加第四项后的问题代码 #

修改后的PDE约束代码：

return m.Re * m.dCdt[t,x] + m.v * m.dCdx[t,x] - m.D * m.d2Cdx2[t,x] - m.Re*m.v_var*m.C[t,x] == 0

解决方法 #

循环引用错误通常源于Pyomo处理DAE时变量依赖关系解析异常，或是边界/初始条件定义不完整，可通过以下步骤修复：

1. 补充完整的初始与边界条件
   原代码仅跳过了t=0、x=0/x=L处的PDE约束，但未单独定义这些位置的条件，导致变量依赖关系混乱。需添加：

   # 初始条件（t=0时的C值，替换为你的实际初始条件）
   @m.Constraint(m.x)
   def initial_condition(m, x):
       return m.C[0, x] == 1.0  # 示例值，按需修改
   
   # 左边界条件（x=0处的C值）
   @m.Constraint(m.t)
   def boundary_left(m, t):
       return m.C[t, 0] == 1.0  # 示例值，按需修改
   
   # 右边界条件（x=L处的C值）
   @m.Constraint(m.t)
   def boundary_right(m, t):
       return m.C[t, L] == 0.0  # 示例值，按需修改
   

2. 确保执行DAE离散化步骤
   Pyomo.DAE模型必须先离散化才能用IPOPT求解，需在求解前添加离散化代码：

   # 采用有限差分法离散化
   discretizer = TransformationFactory('dae.finite_difference')
   # 离散化空间变量x，nfe为有限元步数，可按需调整
   discretizer.apply_to(m, nfe=20, wrt=m.x, scheme='BACKWARD')
   # 离散化时间变量t
   discretizer.apply_to(m, nfe=10, wrt=m.t, scheme='BACKWARD')
   

3. 检查参数与变量初始化
   确保所有参数（D、v、Re等）均已正确赋值，变量C可设置初始值辅助求解器收敛：

   m.C = Var(m.t, m.x, initialize=0.0)  # 添加initialize参数
   

完整可运行代码 #

from pyomo.environ import ConcreteModel, Param, Var, Constraint, SolverFactory
from pyomo.dae import ContinuousSet, DerivativeVar, TransformationFactory

# 替换为你的实际参数值
D = 0.1
v = 0.5
Re = 100
v_var = 0.01
u = 0.2
L = 1.0
t_end = 5.0

m = ConcreteModel()

# 定义参数
m.D= Param(initialize = D)
m.v = Param(initialize = v)
m.Re = Param(initialize = Re)
m.v_var = Param(initialize = v_var)
m.u = Param(initialize = u)

# 定义空间/时间网格
m.x = ContinuousSet(bounds=(0,L))
m.t = ContinuousSet(bounds=(0,t_end))

# 定义解空间变量并初始化
m.C = Var(m.t, m.x, initialize=0.0)

# 定义偏导数
m.dCdt   = DerivativeVar(m.C, wrt=m.t)
m.dCdx   = DerivativeVar(m.C, wrt=m.x)
m.d2Cdx2 = DerivativeVar(m.C, wrt=(m.x, m.x))

# 初始条件
@m.Constraint(m.x)
def initial_condition(m, x):
    return m.C[0, x] == 1.0 if x < 0.5 else 0.0  # 示例初始分布

# 边界条件
@m.Constraint(m.t)
def boundary_left(m, t):
    return m.C[t, 0] == 1.0

@m.Constraint(m.t)
def boundary_right(m, t):
    return m.C[t, L] == 0.0

# 完整PDE约束
@m.Constraint(m.t, m.x)
def pde(m, t, x):
    if t == 0 or x == 0 or x == L:
        return Constraint.Skip
    return m.Re * m.dCdt[t,x] + m.v * m.dCdx[t,x] - m.D * m.d2Cdx2[t,x] - m.Re*m.v_var*m.C[t,x] == 0

# 离散化DAE模型
discretizer = TransformationFactory('dae.finite_difference')
discretizer.apply_to(m, nfe=20, wrt=m.x, scheme='BACKWARD')
discretizer.apply_to(m, nfe=10, wrt=m.t, scheme='BACKWARD')

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

# 输出求解状态
print("求解状态:", result.solver.status)
print("终止原因:", result.solver.termination_condition)

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