Very long running time for Haase and Muller (2014) coded in Python + Gurobi

Question

I have read the paper of Haase and Muller (2014) where they present the linear reformulations of the multinomial logit choice probabilities and compare different approaches. I have tried for an exercise to code their models, particularly the variant of Haase (2009), and to solve it with Gurobi 6.5. The issue is that I obtain really long running times.

In their numerical investigation they used CPLEX 12 on a 64-bit Windows Server 2008 with 1 Intel Xeon 2.4 GHz processor and 24 GiB of RAM. On the other hand I am using Python + Gurobi on a 64-bit Windows 8 with 2 Six-Core AMD Opteron 2.00 GHz processors and 32 GiB of RAM.

I understand that there can be differences between Gurobi and CPLEX in running times on different models/programs. In some cases the Gurobi is better, and in some cases the CPLEX is better. Also, I understand that my machine is not very representative and I hope that 0.4 GHz doesn't make such a big difference. But the difference is really huge. For some small problems, that Haase and Muller needed a little bit more than a minute to solve them, on my machine I need more than two hours. Is it possible that choice of the solver and machine can matter that much? Can someone, please, check my code, because I am newbie in this game.

This is the mathematical model written in LaTeX:

$ \max F = \sum_{i\in I} \sum_{j\in J} x_{ij} $

s.t.

$ \sum_{j\in J} y_j = r $

$ \tilde{x}_{i} + \sum_{j\in J} \leq 1, \forall i \in I $

$ x_{ij} - \frac{\phi_{ij}}{1 + \phi_{ij}} y_j \leq 0, \forall i \in I, j \in J $

$ x_{ij} - \phi_{ij}\tilde{x}_i \leq 0, \forall i \in I, j \in J $

$ y_i \in \{0, 1\}, \forall j \in J $

$ x_{ij} \geq 0, \forall i \in I, j \in J $

$ \tilde{x}_i \geq 0, \forall i \in I $

where $ \phi_{ij} = \frac{ e^{v_{ij}} }{\sum_{k\in M\setminus J} e^{v_{ik}}} $.

Here is the code in Python.

from gurobipy import *
import math
import numpy as np

n = 200 # The number of customers
m = 60 # The number of possible facility locations.
x_coords = 30 * np.random.random_sample(n + m)
y_coords = 30 * np.random.random_sample(n + m)
I = { i: (float(x_coords[i]), float(y_coords[i])) for i in range(n) } # The customers' locations.
M = { j: (float(x_coords[n + j]), float(y_coords[n + j])) for j in range(m) } # The possible facility locations.
J = {i: M[i] for i in range(m-10)} # J is subset of M. The competitor has located facilities in M\J.
MdiffJ = {i: M[m-10+i] for i in range(10)} # Competitor's facility locations. |M\J| = 10 as in the paper.
d = {}
v = {}
for i in I:
    for j in J:
        d[i,j] = abs(I[i][0] - J[j][0]) + abs(I[i][1] - J[j][1]) # Rectangular distances.
        v[i,j] = -0.2 * d[i,j]
r = round(0.2 * len(J)) # == 0.2 * |J|, as in the paper.


model = Model()

# Creating variables
y = {}
x_tilda = {}
x = {}
for j in J:
    y[j] = model.addVar(vtype=GRB.BINARY, name="y_%i" % j)
for i in I:
    x_tilda[i] = model.addVar(lb=0, name="x_tilda_%i" % i)
    for j in J:
        x[i,j] = model.addVar(lb=0, name="x_%i%i" % (i,j))

# Integrating variables
model.update()

# Setting the objective.
obj = LinExpr()
obj += quicksum(quicksum(x[i,j] for j in J) for i in I)

# Setting the objectiv as a maximization problem.
model.setObjective(obj, GRB.MAXIMIZE)

# Computing phi (constants)
phi = {}
for i in I:
    for j in J:
        phi[i,j] = math.e**v[i,j] / (sum(math.e**v[i,k] for k in MdiffJ))

# Adding constraints.
model.addConstr(quicksum(y[j] for j in J) - r == 0, "sum(y_*) == r")
for i in I:
    model.addConstr(x_tilda[i] + quicksum(x[i,j] for j in J) <= 1, "x_tilda_%i+sum(x_%i*) <= 1" % (i,i))
    for j in J:
        model.addConstr(x[i,j] - (phi[i,j] / (1 + phi[i,j])) * y[j] <= 0, "x_%i%i - frac * y_%i <= 0" % (i,j,j))
        model.addConstr(x[i,j] - phi[i,j] * x_tilda[i] <= 0, "x_%i%i - phi_%i%i * x_tilda_%i <= 0" % (i,j,i,j,i))

# Computing the optimal solution.
model.optimize()