Pyomo 中的慢二次约束创建

编程问答 2022-06-02

问题描述

尝试在 Pyomo 中构建大规模二次约束如下：

import pyomo as pyo
from pyomo.environ import *

scale   = 5000
pyo.n   = Set(initialize=range(scale))
pyo.x   = Var(pyo.n,bounds=(-1.0,1.0))
# Q is a n-by-n matrix in numpy array format,where n equals <scale>
Q_values = dict(zip(list(itertools.product(range(0,scale),range(0,scale))),Q.flatten()))
pyo.Q   = Param(pyo.n,pyo.n,initialize=Q_values)

pyo.xQx = Constraint( expr=sum( pyo.x[i]*pyo.Q[i,j]*pyo.x[j] for i in pyo.n for j in pyo.n ) <= 1.0 )

事实证明，考虑到问题的规模，最后一行慢得令人无法忍受。尝试了 PyPSA、Performance of creating Pyomo constraints 和 pyomo seems very slow to write models 中提到的几件事。但没有运气。

有什么建议（一旦构建模型，Ipopt 求解也很慢。但我猜这与 Pyomo 无关）？

ps：直接构造这样的二次约束也没有帮助（也慢得难以忍受）

pyo.xQx = Constraint( expr=sum( pyo.x[i]*Q[i,j]*pyo.x[j] for i in pyo.n for j in pyo.n ) <= 1.0 )

解决方法

您可以通过使用 quicksum 代替 sum 来获得小幅加速。为了衡量最后一行的性能，我稍微修改了你的代码，如下所示：

import itertools
from pyomo.environ import *
import time
import numpy as np

scale = 5000
m = ConcreteModel()
m.n = Set(initialize=range(scale))
m.x = Var(m.n,bounds=(-1.0,1.0))
# Q is a n-by-n matrix in numpy array format,where n equals <scale>
Q = np.ones([scale,scale])
Q_values = dict(
    zip(list(itertools.product(range(scale),range(scale))),Q.flatten()))
m.Q = Param(m.n,m.n,initialize=Q_values)

t = time.time()
m.xQx = Constraint(expr=sum(m.x[i]*m.Q[i,j]*m.x[j]
                            for i in m.n for j in m.n) <= 1.0)
print("Time to make QuadCon = {}".format(time.time() - t))

我用 sum 测量的时间约为 174.4 秒。使用 quicksum 我有 163.3 秒。

对如此微薄的收益不满意，我尝试重新制定为 SOCP。如果您可以像这样分解 Q：Q= (F^T F)，那么您可以轻松地将您的约束表示为二次圆锥，如下所示：

import itertools
import time
import pyomo.kernel as pmo
from pyomo.environ import *
import numpy as np

scale = 5000
m = pmo.block()
m.n = np.arange(scale)
m.x = pmo.variable_list()
for j in m.n:
    m.x.append(pmo.variable(lb=-1.0,ub=1.0))
# Q = (F^T)F factors (eg.: Cholesky factor)
_F = np.ones([scale,scale])
t = time.time()
F = pmo.parameter_list()
for f in _F:
    _row = pmo.parameter_list(pmo.parameter(_e) for _e in f)
    F.append(_row)
print("Time taken to make parameter F = {}".format(time.time() - t))
t1 = time.time()
x_expr = pmo.expression_tuple(pmo.expression(
    expr=sum_product(f,m.x,index=m.n)) for f in F)
print("Time for evaluating Fx = {}".format(time.time() - t1))
t2 = time.time()
m.xQx = pmo.conic.quadratic.as_domain(1,x_expr)
print("Time for quad constr = {}".format(time.time() - t2))

在同一台机器上运行，我观察到准备传递给锥体的表达式的时间约为 112 秒。实际上准备锥体只需要很少的时间（0.031 秒）。

自然，唯一可以处理 pyomo 中圆锥约束的求解器是 MOSEK。 Pyomo-MOSEK 接口的最新更新也显示出有希望的加速。

您可以通过获取 MOSEK trial license 来免费试用 MOSEK。如果您想阅读更多关于圆锥曲线重构的内容，可以在 MOSEK modelling cookbook 中找到快速而全面的指南。最后，如果您隶属于学术机构，那么我们可以为您提供个人/机构学术许可。希望你觉得这有用。

constantsipoptpyomopythonquadratic