设施分散问题：混合整数规划模型
设施分散问题 混合整数规划模型
关于p-分散和最大和模型的全面Python教程
在某些设施位置问题中，需要将设施布置在一起，以免一个设施的影响使其他设施被掩盖或不利地影响。无论是出于风险缓解、避免竞争还是其他考虑，掌握解决这些挑战的方法是运筹学工具包中的关键技能。
Kuby（1987）为这些问题提出了两种不同的混合整数规划（MIP）模型：p-分散问题和最大和问题。在本文中，将使用Python库Pyomo和HiGHS求解器来实现这两个模型。
除了这两个模型，还将介绍一些有用的建模资源。首先，介绍了一种在线性规划上下文中线性化二进制变量乘积的策略，尽管考虑到最大化目标，但在这个问题中不需要显式表示。其次，介绍了一种最大-最小的MIP模型，如果选中该项，则意图最大化小于一组项中的任何参数的某个东西。最后，介绍了一种解决具有明确定义的优先级层次结构的多目标问题的策略，该策略结合了这两个模型的元素。
那些对数值优化尚不熟悉的人在继续阅读本文之前，可以先看一下我之前关于线性规划和分支限界方法的文章。
通常情况下，您可以在此git存储库中找到完整的代码。
您会选择哪些位置来放置5个设施中的哪一个？
二进制变量的乘积
在定义这个问题的基本元素时，可以使用二进制变量，如果选择了一个位置，则二进制变量的值为1，否则为0。让我们将这些变量表示为xᵢ。假设预先计算了两个位置之间的距离（或其他分散度量），让我们将它们表示为dᵢⱼ。我们如何计算所选任意一对设施的分散度？
在这种情况下，使用二进制变量xᵢ和xⱼ的乘积来计算它们同时包含在解决方案中时得到的不相似性是比较直观的。然而，这将导致一个二次规划问题，因此无法应用线性求解器。幸运的是，有一种方法可以在MIP中使用几个约束来表示二进制变量的乘积。
考虑一个有向图G(V, A)，zᵢⱼ是一个新的二进制变量，表示节点i和j都被选中。如果i或j中有一个未被选择，则zᵢⱼ必须为0。这导致了第一组约束：
到目前为止，即使i和j都被选择，zᵢⱼ仍然可以为0。因此，我们必须包括一个额外的约束条件，当i和j被选择时，zᵢⱼ变为1。
当最大化与zᵢⱼ成比例的某个值时，如最大和问题中的情况，第二组约束是不必要的，因为如果它是可行的，那么不计算与zᵢⱼ成比例的增益是没有意义的。然而，在制定其他MIP模型时，这可能会有用。
在下面的部分中，让我们将这些概念应用到最大和问题。
最大和模型
离散最大和问题必须定义在N个离散节点中选择p个设施的位置，以最大化所选设施之间计算的距离之和（或距离之平均值）。因此，考虑将设施视为位于有向图G（V，A）中的节点。从i到j的每个弧的权重是事先已知的距离（离散度）度量dᵢⱼ。此外，考虑使用二进制变量xᵢ指示是否选择位置i，以及zᵢⱼ指示是否同时选择i和j。
目标可以表示为以下内容：
除了在上一节中提出的线性化二进制变量乘积的约束条件外，还需要包括一个约束条件来确保选择了p个位置。
因此，我们可以将问题的约束条件陈述如下：
让我们使用Python将其编写成代码。为此，我们必须首先导入将要使用的库。numpy库将用于线性代数计算和创建随机坐标点；scipy库中的squareform和pdist函数将用于根据坐标矩阵计算距离；matplotlib将是我们的可视化工具；pyomo将是代数建模语言（AML）（使用HiGHS求解器）。
import numpy as npfrom scipy.spatial.distance import squareform, pdistimport matplotlib.pyplot as pltimport pyomo.environ as pyofrom pyomo.contrib.appsi.solvers.highs import Highs我们需要创建N个点，并确定其中有多少个必须被选为设施位置。通过在每次代码执行中固定随机种子（12），可以获得引言中表示的点。
# 固定随机种子np.random.seed(12)# 创建随机点N = 25p = 5coordinates = np.random.random((N, 2))# 计算两两之间的距离distances = squareform(pdist(coordinates))现在我们有了开始构建pyomo模型所需的元素。
在pyomo中，对于建模问题有两种方法：抽象模型和具体模型。在第一种方法中，在提供某些数据值之前定义了问题的代数表达式，而在第二种方法中，模型实例在其元素被定义的同时就被创建。您可以在库文档或Bynum等人的书中找到有关这些方法的更多信息。在本文中，我们将采用具体模型的形式。
model = pyo.ConcreteModel()在这个模型中，我们有两个集合：节点和弧。由于我们考虑的是一个完全有向图，所以除了从节点到自身的节点之外，每对节点之间都存在弧。
# 节点和弧的集合model.V = pyo.Set(initialize=range(N))model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j])预先提供的参数是必须选择的节点数和节点对之间的距离。
# 参数model.d = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A})model.p = pyo.Param(initialize=p)然后，让我们包括决策变量。
# 决策变量model.x = pyo.Var(model.V, within=pyo.Binary)model.z = pyo.Var(model.A, within=pyo.Binary)以及约束条件…
# 选择p个节点def p_selection(model): return sum(model.x[:]) == model.p# 如果起始节点未被选择，则弧为0def dispersion_c1(model, i, j): return model.z[i, j] <= model.x[i]# 如果结束节点未被选择，则弧为0def dispersion_c2(model, i, j): return model.z[i, j] <= model.x[j]# 将约束条件包括在模型中model.p_selection = pyo.Constraint(rule=p_selection)model.dispersion_c1 = pyo.Constraint(model.A, rule=dispersion_c1)model.dispersion_c2 = pyo.Constraint(model.A, rule=dispersion_c2)最后，我们必须创建目标函数。
def disp_obj(model): return sum(model.z[i, j] * model.d[i, j] for (i, j) in model.A)model.obj = pyo.Objective(rule=disp_obj, sense=pyo.maximize)现在最大值模型已准备好进行求解。然后，我们必须实例化与Pyomo兼容的求解器来处理它。 Highs求解器自版本6.4.3以来已在Pyomo中提供（请检查导入），如果安装了highspy软件包。因此，请确保运行pip install highspy。
solver = Highs()solver.solve(model)经过约120秒的计算时间后，我们得到以下结果：
请注意，尽管总离散度被最大化，但左上角的两个设施仍然非常接近，这可能是不可取的。因此，作为最大值模型的替代，我们有p-离散模型，其中我们最大化所选节点之间的最小距离。
p-离散模型
在p-离散模型中，我们需要一个额外的决策变量来计算所选节点之间的最小距离，这是我们要最大化的目标。让我们将此变量记为D。我们必须创建一个大M约束，以确保如果i和j都被选择，则D小于或等于dᵢⱼ；否则，我们必须确保D无限大。如果我们保持zᵢⱼ作为二进制变量的乘积来表达，则可以将此附加约束表达如下。
在此公式中，M是一个任意大的固定参数，用于制定一个分离规则。M的良好选择应足够大，以确保在zᵢⱼ为零时对任何D的值都保持约束（i，j）的可行性，但尽可能小，以使模型的线性松弛与整数版本相似，从而更容易收敛。在这个问题中，dᵢⱼ的最大值可能是一个有趣的替代方案。
尽管此公式对于该模型效果很好，但可以使用更简明的方法，其中甚至不包括变量zᵢⱼ在模型中。
在此公式中，xᵢ或xⱼ等于零是确保不等式对任何D的值都有效的充分条件。目标变为简单地最大化D。
在下面的Python代码中，考虑我们有一个与前一个模型相同的集合和参数以及决策变量x的新模型。
# 最大-最小约束def maxmin_rule(model, i, j): return model.D <= model.d[i, j] + model.M * (1 - model.x[i]) + model.M * (1 - model.x[j])# 新参数 big Mmodel.M = max(model.d[:, :])# 新变量model.D = pyo.Var(within=pyo.NonNegativeReals)# 新约束model.maxmin_rule = pyo.Constraint(model.A, rule=maxmin_rule)# 目标model.obj = pyo.Objective(expr=model.D, sense=pyo.maximize)在调用求解器之后，花费不到1.2秒即可获得以下结果。
这似乎很好，因为位置在空间中分布得很均匀。
有没有办法改善这种分布呢？
多准则方法
请记住p-分散模型的目标函数仅取决于所选节点之间的最小距离。因此，可以使用定义此距离的两个点的组合以及与目标自身大于或等于彼此距离的其他点来获得多个解。我们可以通过选择这些替代方案中的最佳方案来改进我们的解决方案吗？这导致了Kuby（1987）提出的两步“多准则方法”。
第一步，完整解决p-分散模型，并将当前目标存储在参数d_opt中。然后，求解具有附加约束D ≥ d_opt的maxisum模型，以获得在p-分散模型的最优解中，总分散度最大的解之一。
在Python中执行此操作，假设您已经使用maxisum模型的决策变量和约束实例化了p-分散模型。
# D必须是最优的def composed_constr(model): return model.D >= model.d_opt# 解决p-分散问题solver.solve(model)# 新参数model.d_opt = pyo.Param(initialize=model.obj())# 停用旧的目标函数model.obj.deactivate()# 更多变量model.z = pyo.Var(model.A, within=pyo.Binary)# 解决方案不会使当前的D变坏model.composed_constr = pyo.Constraint(rule=composed_constr)# 新目标函数model.obj_disp = pyo.Objective(rule=disp_obj, sense=pyo.maximize)solver.solve(model)这个过程非常快，因为当求解器进入第二个目标函数时，可行解空间大大缩小。不到一秒的时间，就能得到以下结果。
进一步阅读
当客户不均匀分布、设施容量有限或者适当数量的设施事先未知时，可能面临不同的设施位置问题。您可以在Nicolo Cosimo Albanese的这篇令人惊叹的故事中找到使用Python和PuLP实现的Balinski（1965）的公式。
优化：Python中的容量设施位置问题
找到最佳的仓库数量和位置，以降低成本并满足需求
towardsdatascience.com
结论
在本文中，使用Pyomo在Python中实现了两个混合整数规划模型，用于设施分散问题。与文献中先前验证的结果一样，maxisum模型可能导致元素在空间中分布不均匀，而p-分散模型产生的解决方案位置在空间中分布均匀。可以使用maxisum目标函数通过从最优解集合中选择具有最大总分散度的解来优化p-分散模型的解决方案。这些示例中使用的完整代码可供进一步使用。
参考资料
Balinski, M. L. 1965. Integer programming: methods, uses, computations. Management science, 12(3), 253–313.
Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., … & Woodruff, D. L. 2021. Pyomo-optimization modeling in python (Vol. 67). Berlin/Heidelberg, Germany: Springer.
Kuby, M. J. 1987. Programming models for facility dispersion: The p‐dispersion and maxisum dispersion problems. Geographical Analysis, 19(4), 315–329.
