尝试在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中提到的几件事，创建Pyomo约束的性能和Pyomo编写模型的速度似乎非常慢。但是没有运气。
任何建议（一旦模型被构建，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.031s）。
当然，在pyomo中，唯一能够处理圆锥约束的解算器是MOSEK。Pyomo-MOSEK界面的最新更新也显示出了良好的加速效果。
您可以通过获得MOSEK试用许可证免费试用MOSEK。如果你想阅读更多关于圆锥曲线重新配方的信息，可以在MOSEK建模食谱中找到一个快速而全面的指南。最后，如果您隶属于学术机构，那么我们可以为您提供个人/机构学术许可证。希望你觉得这有用。
