如何解决Pyomo 中的慢二次约束创建
尝试在 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 中找到快速而全面的指南。最后,如果您隶属于学术机构,那么我们可以为您提供个人/机构学术许可。希望你觉得这有用。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。