我读过Haase 和 Muller (2014)的论文，他们介绍了多项式 logit 选择概率的线性重新表述，并比较了不同的方法。我尝试了一个练习来编写他们的模型，特别是Haase (2009)的变体，并用 Gurobi 6.5 来解决它。问题是我获得了非常长的运行时间。

在他们的数值研究中，他们在具有 1 个 Intel Xeon 2.4 GHz 处理器和 24 GiB RAM 的 64 位 Windows Server 2008 上使用了 CPLEX 12。另一方面，我在具有 2 个六核 AMD Opteron 2.00 GHz 处理器和 32 GiB RAM 的 64 位 Windows 8 上使用 Python + Gurobi。

我了解 Gurobi 和 CPLEX 在不同模型/程序上的运行时间可能存在差异。在某些情况下，Gurobi 更好，在某些情况下，CPLEX 更好。另外，我知道我的机器不是很有代表性，我希望 0.4 GHz 不会产生如此大的差异。但是差别真的很大。对于一些小问题，哈斯和穆勒需要一分钟多一点的时间来解决它们，在我的机器上我需要两个多小时。求解器和机器的选择有可能那么重要吗？有人可以检查我的代码，因为我是这个游戏的新手。

这是用 LaTeX 编写的数学模型：

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

英石

$\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$

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

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