Pyomo模型评估耗时优化求助：两种约束实现均无改善

我有一个ConcreteModel，采用完全向量化的Pyomo ConstraintList构建约束，IPOPT求解仅需约8秒，但模型评估耗时约40秒。以下是report_timing输出及约束代码片段：

[+   0.20] Processed 1 objectives
[+   0.20] Processed 1 objectives
[+   0.20] Processed 1 objectives
[+   0.20] Processed 1 objectives
[+  23.26] Processed 12410 constraints
[+  23.26] Processed 12410 constraints
[+  23.26] Processed 12410 constraints
[+  23.26] Processed 12410 constraints
[   24.22] Generated NL representation
[   24.22] Generated NL representation
[   24.22] Generated NL representation
[   24.22] Generated NL representation
24.29 seconds required to write file
24.29 seconds required for presolve

ConstraintList构建约束的代码片段 #

# Define the constraints
pmo.constraints = pyo.ConstraintList()

for i in pmo.i:
    force = list()
    for mi in pmo.mi:
        # (#1) geometry constraint
        pmo.constraints.add( pmo.MLAT[i,mi]**2 + pmo.height[mi]**2 == (pmo.NML[i,mi]*pmo.MLopt[mi])**2 )

        # (#2) length velocity constraint
        if i == 0:     # second order forward difference
            _NMLDot = ((-3*pmo.NML[i,mi]+4*pmo.NML[i+1,mi]-pmo.NML[i+2,mi])/2/iDt)
        elif i == 1:    # first order forward difference
            _NMLDot = ((pmo.NML[i+1,mi]-pmo.NML[i,mi])/1/iDt)
        elif i == iN-2:  # first order backward difference
            _NMLDot = ((pmo.NML[i,mi]-pmo.NML[i-1,mi])/1/iDt)
        elif i == iN-1:  # second order backward difference
            _NMLDot = ((3*pmo.NML[i,mi]-4*pmo.NML[i-1,mi]+pmo.NML[i-2,mi])/2/iDt)
        else:           # central difference
            _NMLDot = ((pmo.NML[i+1,mi]-pmo.NML[i-1,mi])/2/iDt)
        pmo.constraints.add( _NMLDot / pmo.CVmax[mi] == pmo.NMV[i,mi])

        # (#3) force equilibrium [external function]
        TF, MFAT = utils.calcForcePyomo(
                pmo.act[i,mi], pmo.NML[i,mi], pmo.NMV[i,mi], pmo.MLAT[i,mi], \
                pmo.MTL[i,mi], pmo.IFmax[mi], pmo.MLopt[mi], pmo.TLslk[mi] )
        pmo.constraints.add( TF == MFAT )
        force.append(MFAT)

    # (#4) moment equilibrium
    for cName in pmo.cName:
        _tau = iTauT.getDependentColumn(cName)[i] # float
        _moment = pyo.quicksum( iMomentArm[cName][i,mi] * force[mi] for mi in pmo.mi)
        pmo.constraints.add( _moment == _tau )

我还改写了等效的规则式约束，但计算时间并未得到改善：

规则式约束的代码片段 #

def constraint_1(model, i, mi):
    return model.MLAT[i,mi]**2 + model.height[mi]**2 == (model.NML[i,mi]*model.MLopt[mi])**2

def constraint_2(model, i, mi):
    if i == 0:     # second order forward difference
        _NMLDot = ((-3*model.NML[i,mi]+4*model.NML[i+1,mi]-model.NML[i+2,mi])/2/iDt)
    elif i == 1:    # first order forward difference
        _NMLDot = ((model.NML[i+1,mi]-model.NML[i,mi])/1/iDt)
    elif i == iN-2:  # first order backward difference
        _NMLDot = ((model.NML[i,mi]-model.NML[i-1,mi])/1/iDt)
    elif i == iN-1:  # second order backward difference
        _NMLDot = ((3*model.NML[i,mi]-4*model.NML[i-1,mi]+model.NML[i-2,mi])/2/iDt)
    else:           # central difference
        _NMLDot = ((model.NML[i+1,mi]-model.NML[i-1,mi])/2/iDt)
    
    return (_NMLDot / model.CVmax[mi] == model.NMV[i,mi])

pmo.eq_1 = pyo.Constraint(pmo.i, pmo.mi, rule=constraint_1)
pmo.eq_2 = pyo.Constraint(pmo.i, pmo.mi, rule=constraint_2)


# store the output of the external function
pmo.TF, pmo.MFAT = {},{}

def buildForces(model, i, mi):
    model.TF[i,mi], model.MFAT[i,mi] = utils.calcForcePyomo(
            model.act[i,mi], model.NML[i,mi], model.NMV[i,mi], model.MLAT[i,mi], \
            model.MTL[i,mi], model.MLopt[mi], pmo.TLslk[mi] )

pmo.updForces = pyo.BuildAction(pmo.i, pmo.mi, rule=buildForces)


def constraint_3(model, i, mi):
    return model.TF[i,mi] == model.MFAT[i,mi]

def constraint_4(model, i, cName):
    _tau       = iTauT.getDependentColumn(cName)[i] # float
    moment     = pyo.quicksum( iMomentArm[cName][i,mi] * model.MFAT[i,mi] for mi in model.mi)
    return moment == _tau

pmo.eq_3 = pyo.Constraint(pmo.i, pmo.mi,    rule=constraint_3)
pmo.eq_4 = pyo.Constraint(pmo.i, pmo.cName, rule=constraint_4)

优化建议 #

从timing报告看，大部分时间消耗在约束处理和NL文件生成环节，以下是针对性优化方向：

• 优化外部函数calcForcePyomo
  这是核心瓶颈：当前逐点调用的方式会产生大量Python循环开销。

  • 改为向量化计算：将整个维度的变量（如pmo.act[:,mi]）一次性传入函数，用NumPy向量化操作替代逐元素循环，大幅减少函数调用次数。
  • 预计算结果：提前计算所有(i,mi)点的TF和MFAT，存储为Pyomo的Param或Expression，避免在约束构建阶段重复计算。

• 提升NL文件生成效率
  Pyomo生成NL文件的耗时直接关联约束复杂度和数量：

  • 替换quicksum：大规模模型中，Python原生sum的性能优于quicksum，或提前将求和项定义为Expression复用。
  • 清理冗余约束：检查是否存在重复定义或可合并的约束项，减少约束总数。

• 采用向量化约束定义
  对于约束1、约束2这类结构化强的约束，直接利用Pyomo集合操作实现向量化定义，避免Python循环：

  pmo.eq_1 = pyo.Constraint(pmo.i, pmo.mi, expr=pmo.MLAT**2 + pmo.height**2 == (pmo.NML*pmo.MLopt)**2)
  

• 预计算固定参数
  像iTauT.getDependentColumn(cName)[i]这类固定值，提前计算并存储为Pyomo的Param，避免在约束构建或评估时重复调用外部接口。

• 启用Pyomo快速模式
  创建模型时开启快速模式，减少内部运行开销：

  pmo = pyo.ConcreteModel(fast=True)
  

  或设置环境变量PYOMO_FAST=1。

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