JoseleMG JoseleMG - 1 year ago 200
Python Question

Using multiprocess to build Gurobi expression in Python

First of all, I am working with jupyter-notebook, python version 3.5, and Gurobi 7.0.2 with its python interface, all on Red Hat.

This is the context of my problem: I want to solve a Quadratic problem, which has a huge number of variables. It takes over 1-2 hours, to build the objective function.

I thought about using NumPy GPU acceleration but the expression is a little bit tricky, and so, this can't be a solution.

Therefore, I am trying to build the objective function using several threads. However, I am getting an error that I don't know how to handle it.

I simplified my code, so it can be more readable (The error is still the same).

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())


Output:

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


I guess from this thread https://groups.google.com/forum/#!topic/gurobi/fwLRrWLLJqo, that it is something related with pickling the gurobi expression that I output from the queue, but I am not really sure.

Do you know how can I solve this? Is there any other way to "return" results from process that could work here?. I would like to avoid writing in the disk, due to it slowness (maybe the last resource :S).

Thank you in advance :).

P.D. The slow part of my code is this one, which I tried to split it in several process:

# 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

Answer Source

First, Gurobi Optimizer does not support multiple threads for model building. Even if we did, model building is almost never the bottleneck in properly written applications.

In this case, you have a lot of expressions that look like x + x + x. While correct, that is very inefficient - it is much better to write 3*x. Here is a quick rewrite of 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

And a quick workaround for the issue raised by Erwin Kavelagen is to set PreQLinearize=1.

Edit: We can make penalty_function slightly more efficient by combining the two var[i]*var[j] terms:

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
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download