调用Simulator求解Trambouze反应Pyomo模型时触发TypeError错误

问题背景 #

使用Pyomo定义Trambouze反应的PFR模型，调用Simulator并指定scipy包模拟时触发TypeError: 'float' object is not callable错误，相关代码及错误信息如下：

模型定义代码 #

# kinetic constants
k1 = 0.0001 # mol/dm3*s
k2 = 0.0015 # s-1
k3 = 0.008 # dm3/mol*s

# create pyomo model
def PFR_trambouze(CAo = 0.112, CBo = 0.132, CXo = 0.0783, CYo = 0.0786):
    m = ConcreteModel()

    # define time variable and bounds
    m.t = ContinuousSet(bounds=(0.0, 300.0))
    # define species concentrations
    m.CA = Var(m.t)
    m.CB = Var(m.t)
    m.CX = Var(m.t)
    m.CY = Var(m.t)

    # define differential variables
    m.dCA = DerivativeVar(m.CA)
    m.dCB = DerivativeVar(m.CB)
    m.dCX = DerivativeVar(m.CX)
    m.dCY = DerivativeVar(m.CY)

    # initial conditions
    m.CA[0.0].fix(CAo)
    m.CB[0.0].fix(CBo)
    m.CX[0.0].fix(CXo)
    m.CY[0.0].fix(CYo)

    # differential equations
    @m.Constraint(m.t)
    def ode_A(m, t):
        return m.dCA[t] == - k1 - k2*m.CA[t] - k3*m.CA[t]**2

    @m.Constraint(m.t)
    def ode_B(m, t):
        return m.dCB[t] == k1*m.CA[t]
    
    @m.Constraint(m.t)
    def ode_X(m, t):
        return m.dCX[t] == k1
    
    @m.Constraint(m.t)
    def ode_Y(m, t):
        return m.dCY[t] == k3*m.CA[t]**2
    
    return m

模拟调用代码 #

sim = Simulator(PFR_trambouze(), package = 'scipy')
tsim, profile = sim.simulate(numpoints= 100)

错误栈信息 #

TypeError                                 Traceback (most recent call last)
<ipython-input-29-00475f40d0ce> in <cell line: 2>()
      1 sim = Simulator(PFR_trambouze(), package = 'scipy')
----> 2 tsim, profile = sim.simulate(numpoints= 100)

4 frames
/usr/local/lib/python3.10/dist-packages/pyomo/dae/simulator.py in simulate(self, numpoints, tstep, integrator, varying_inputs, initcon, integrator_options)
    921                     "The scipy module is not available. Cannot simulate the model."
    922                 )
---> 923             tsim, profile = self._simulate_with_scipy(
    924                 initcon, tsim, switchpts, varying_inputs, integrator, integrator_options
    925             )

/usr/local/lib/python3.10/dist-packages/pyomo/dae/simulator.py in _simulate_with_scipy(self, initcon, tsim, switchpts, varying_inputs, integrator, integrator_options)
    957                         p.set_value(varying_inputs[v][tsim[i - 1]])
    958 
---> 959             profilestep = scipyint.integrate(tsim[i])
    960             profile = np.vstack([profile, profilestep])
    961             i += 1

/usr/local/lib/python3.10/dist-packages/scipy/integrate/_ode.py in integrate(self, t, step, relax)
    429 
    430         try:
---> 431             self._y, self.t = mth(self.f, self.jac or (lambda: None),
    432                                   self._y, self.t, t,
    433                                   self.f_params, self.jac_params)

/usr/local/lib/python3.10/dist-packages/scipy/integrate/_ode.py in run(self, f, jac, y0, t0, t1, f_params, jac_params)
   1342         args = [f, y0, t0, t1] + self.call_args[:-1] + \
   1343                [jac, self.call_args[-1], f_params, 0, jac_params]
-> 1344         y1, t, istate = self.runner(*args)
   1345         self.istate = istate
   1346         if istate < 0:

/usr/local/lib/python3.10/dist-packages/pyomo/dae/simulator.py in _rhsfun(t, x)
    655 
    656                 for d in derivlist:
---> 657                     residual.append(rhsdict[d]())
    658 
    659                 return residual

TypeError: 'float' object is not callable

错误原因 #

问题出在全局定义的动力学常数k1、k2、k3被直接当作浮点数写入约束表达式。Pyomo的Simulator在处理ODE时，需要将约束右侧转换为可调用的符号表达式对象；但对于ode_X这类右侧是纯常数（k1）的约束，Pyomo会直接将其计算为浮点数，而非保留可调用的表达式结构。当Simulator尝试调用这个浮点数时，就触发了类型错误。

解决方法 #

将全局的动力学常数改为Pyomo模型内的Param（参数），让Pyomo将其视为符号化的参数，确保所有约束右侧都是可调用的表达式：

修改后的模型代码：

# create pyomo model
def PFR_trambouze(CAo = 0.112, CBo = 0.132, CXo = 0.0783, CYo = 0.0786):
    m = ConcreteModel()

    # 定义动力学参数（替换全局变量）
    m.k1 = Param(initialize=0.0001, doc='mol/dm3*s')
    m.k2 = Param(initialize=0.0015, doc='s-1')
    m.k3 = Param(initialize=0.008, doc='dm3/mol*s')

    # define time variable and bounds
    m.t = ContinuousSet(bounds=(0.0, 300.0))
    # define species concentrations
    m.CA = Var(m.t)
    m.CB = Var(m.t)
    m.CX = Var(m.t)
    m.CY = Var(m.t)

    # define differential variables
    m.dCA = DerivativeVar(m.CA)
    m.dCB = DerivativeVar(m.CB)
    m.dCX = DerivativeVar(m.CX)
    m.dCY = DerivativeVar(m.CY)

    # initial conditions
    m.CA[0.0].fix(CAo)
    m.CB[0.0].fix(CBo)
    m.CX[0.0].fix(CXo)
    m.CY[0.0].fix(CYo)

    # differential equations（引用模型内的参数m.k1/m.k2/m.k3）
    @m.Constraint(m.t)
    def ode_A(m, t):
        return m.dCA[t] == - m.k1 - m.k2*m.CA[t] - m.k3*m.CA[t]**2

    @m.Constraint(m.t)
    def ode_B(m, t):
        return m.dCB[t] == m.k1*m.CA[t]
    
    @m.Constraint(m.t)
    def ode_X(m, t):
        return m.dCX[t] == m.k1
    
    @m.Constraint(m.t)
    def ode_Y(m, t):
        return m.dCY[t] == m.k3*m.CA[t]**2
    
    return m

修改后，所有约束右侧都是基于Pyomo符号参数的表达式，Simulator可以正确提取并调用这些表达式进行数值积分，不会再触发类型错误。

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