Pyomo中缓慢的二次约束创建

Question

尝试在 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 )

Answer 1

您可以使用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 中找到快速而全面的指南。最后，如果您隶属于学术机构，那么我们可以为您提供个人/机构学术许可证。希望你觉得这很有用。