在 Python 中使用多进程构建 Gurobi 表达式

Using multiprocess to build Gurobi expression in Python

首先,我正在使用 jupyter-notebook,python 3.5 版,以及具有 python 界面的 Gurobi 7.0.2,所有这些都在 Red Hat 上。

这是我的问题的上下文:我想求解一个包含大量变量的二次方程式问题。构建 objective 函数需要 1-2 个小时。

我考虑过使用 NumPy GPU 加速,但表达式有点棘手,因此,这不是解决方案。

因此,我尝试使用多个线程构建 objective 函数。但是,我收到一个错误,我不知道如何处理它。

我简化了我的代码,所以它可以更具可读性(错误仍然是一样的)。

from gurobipy import *
import multiprocessing as mp
import queue
mp.set_start_method('fork')

def function(obj,q):

    print('We enter')
    obj = x*x + x*y + y*y + y*z + z*z + 2*x
    q.put(obj)
    print('We end')

m = Model("qp")
obj = QuadExpr()

x = m.addVar(ub=1.0, name="x")
y = m.addVar(ub=1.0, name="y")
z = m.addVar(ub=1.0, name="z")
q = mp.Queue()
if __name__ == '__main__':
    for k in range (4): 
        p = mp.Process(target=function, args=(obj,q,))
        p.start()
        obj+=q.get()
        p.join()
m.setObjective(obj)

# Add constraint: x + 2 y + 3 z <= 4
m.addConstr(x + 2 * y + 3 * z >= 4, "c0")

# Add constraint: x + y >= 1
m.addConstr(x + y >= 1, "c1")

m.optimize()

for v in m.getVars():
    print('%s %g' % (v.varName, v.x))

print('Obj: %g' % obj.getValue())

输出:

We enter
We end

---------------------------------------------------------------------------
GurobiError                               Traceback (most recent call last)
<ipython-input-31-c71f0667f39b> in <module>()
     17         p = mp.Process(target=function, args=(obj,q,))
     18         p.start()
---> 19         obj+=q.get()
     20         p.join()
     21 m.setObjective(obj)

quadexpr.pxi in gurobipy.QuadExpr.__iadd__ (../../src/python/gurobipy.c:39916)()

quadexpr.pxi in gurobipy.QuadExpr.add (../../src/python/gurobipy.c:35739)()

linexpr.pxi in gurobipy.LinExpr.add (../../src/python/gurobipy.c:29245)()

GurobiError: Unsupported type (<class 'NoneType'>) for LinExpr addition argument

我从这个帖子 https://groups.google.com/forum/#!topic/gurobi/fwLRrWLLJqo 猜想这与酸洗我从队列输出的 gurobi 表达式有关,但我不太确定。

你知道我该如何解决这个问题吗?有没有其他方法可以从可以在这里工作的过程中获得 "return" 结果?。我想避免写入磁盘,因为它很慢(可能是最后一个资源:S)。

提前谢谢你:)。

P.D。我的代码中最慢的部分是这个,我试图将它分成几个过程:

# var is an array of GRB.BINARY
# D=edge_costs
def penalty_function(var,obj,D):
    num_nodes = len(var)
    for i,fil in enumerate(D):
        for j,val in enumerate(fil):
            # -x*D*x
            if val!=0:
                obj+=var[i]*var[j]*val
            # -x(i)x(j)*min(Ds)
            if(j>i):
                for k in range(num_nodes):
                    if(not j==i):
                    minval=min(D[j][k],D[i][k])
                        if (minval!=0 ):
                            obj+=var[i]*var[j]*minval
    return obj

首先,Gurobi优化器不支持多线程建模。即使我们这样做了,模型构建也几乎从来不是正确编写应用程序的瓶颈。

在这种情况下,您有很多看起来像 x + x + x 的表达式。虽然正确,但效率很低——最好写成 3*x。这是 penalty_function 的快速重写:

# var is an array of GRB.BINARY
# D=edge_costs
def penalty_function(var,obj,D):
    num_nodes = len(var)
    for i,fil in enumerate(D):
        for j,val in enumerate(fil):
            # -x*D*x
            if val!=0:
                obj+=var[i]*var[j]*val
            # -x(i)x(j)*min(Ds)
            if(j>i):
                minval = sum(min(D[j][k],D[i][k]) for k in range(num_nodes))
                if minval != 0:
                    obj += var[i]*var[j]*minval
    return obj

Erwin Kavelagen 提出的问题的快速解决方法是 set PreQLinearize=1

编辑:我们可以通过组合两个 var[i]*var[j] 项使 penalty_function 稍微更有效率:

def penalty_function(var,obj,D):
    num_nodes = len(var)
    for i,fil in enumerate(D):
        for j,val in enumerate(fil):
            # -x(i)x(j)*min(Ds)
            if(j>i):
                val += sum(min(D[j][k],D[i][k]) for k in range(num_nodes))
            # -x*D*x
            if val!=0:
                obj += var[i]*var[j]*val
    return obj

添加 Greg Glockner 改进后,我意识到我可以以不同的方式使用 numpy 和多线程。我在这里添加它以防有人对此感兴趣。

def parallel_function(D,num_nodes,start,end):     
    minD=[]
    if(end>num_nodes-1):
        end=num_nodes
    for i in range(start,end,1):
        # fill with 0 (triangular matrix)
        minval=[0]*(i+1)
        for j,val in enumerate(D[i]):
            if(j>i):
                minval.append(sum(min(D[j][k],D[i][k]) for k in range(num_nodes)))
        minD.append(minval)
    return minD

def penalty_function(var,obj,D,num_nodes):
    num_process=8
    inc=math.ceil(num_nodes/num_process)
    if __name__ == '__main__':
        pool = mp.Pool(processes=num_process)              # start 8 worker processes
        minval=pool.starmap(parallel_function, [(D,num_nodes,u*inc,u*inc+inc) for u in range(8)])
        pool.close()
    minD=[]
    for x in minval:
        minD=minD+x
    obj+=np.array(var).dot((np.array(minD)+np.array(D))).dot(var)
    return obj