I've been trying to use cplex to solve an optimal transportation problem. The problem model is normally very large (in my description below the total number of variables are 1048576 (= 1024^2), and the number of constraints is 2048). My problem is that the process of adding constraints is far too slow to be practical ( the time spent in solving the model is fine though). I googled this issue, there are some tips, but still I couldn't find a feasible solution.
The problem is as follows: Given two nonnegative vectors a and b of the same length 1024, and a 1024-by-1024 nonnegative matrix C. Assuming that the sum over all elements of a is the same as that of b (np.sum(a) == np.sum(b)). I want to find a 1024-by-1024 nonnegative matrix X such that the sum of C[i,j] * X[i,j] is minimized, subject to the constraints that the sum of all elements of the i-th row of X equals to the i-th element of a and the sum of all elements of the j-th column of X equals to the j-th element of b, for all possible i and j, i.e.
Minimize:
C[0,0] * X[0,0] + C[0,1] * X[0,1] + ... + C[1023,1023] * X[1023,1023]
Subject to:
All X[i,j] >= 0
X[0,0] + X[0,1] + ... + X[0,1023] == a[0]
X[1,0] + X[1,1] + ... + X[1,1023] == a[1]
...
X[1023,0] + X[1023,1] + ... X[1023,1023] == a[1023]
X[0,0] + X[1,0] + ... + X[1023,0] == b[0]
X[0,1] + X[1,1] + ... + X[1023,1] == b[1]
...
X[0,1023] + X[1,1023] + ... X[1023,1023] == b[1023]
My code is roughly like: (in the following code DOT is the transportation model; a and b are lists with length 1024, C is a list with length 1048576(= 1024 ** 2).
from __future__ import print_function
import cplex
DOT = cplex.Cplex()
DOT.objective.set_sense(DOT.objective.sense.minimize)
size = 1024
# set up names of variables
nms = ["x{0}".format(i) for i in range(size * size)]
# add variables to the model
DOT.variables.add(obj = C, names = nms) # C is a nonnegative list with length 1048576
constraints = list()
for i in range(size):
constraints.append([nms[i * size : (i + 1) * size], [1.0] * size])
for i in range(size):
constraints.append(cplex.SparsePair(nms[i : size * size : size], [1.0] * size))
rhs = a + b # a, b are nonnegative lists with the same length and sum
constraint_senses = ["E"] * (size * 2)
# the following line: adding constraints to model is too slow
DOT.linear_constraints.add(lin_expr = constraints, senses = constraint_senses, rhs = rhs)
# solve the model
DOT.solve()
# print some information
print("Solution status :", DOT.solution.get_status())
print("Cost : {0:.5f}".format(DOT.solution.get_objective_value()))
print()
As I write in the comment, the process of adding constraints to the model is too slow. Is there any way to speed up it?
Any help will be appreciated. Thanks in advance!