选址问题模型验证01： The cycle hub location problem

参考文献：

Contreras, I., Tanash, M., & Vidyarthi, N. (2017). Exact and heuristic approaches for the cycle hub location problem. Annals of Operations Research, 258(2), 655-677.

1 背景介绍

Hub location problems (HLPs) arise in the design of hub-and-spoke networks. They have a wide variety of applications in airline transportation, freight transportation, rapid transit systems, trucking industries, postal operations, and telecommunications networks. These systems serve demand for transportation of passengers, commodities, and/or transmission of information (data, voice, video) between multiple origins and destinations. Instead of connecting every origin-destination (O–D) pair directly, hub-and-spoke networks serve customers via a small number of links, where hub facilities consolidate the flows from many origins, transfer them through the hub level network, and eventually distribute them to their final destinations. The use of fewer links in the network concentrates flows at the hub facilities, allowing economies of scale to be applied on routing costs, besides helping to reduce setup costs and to centralize commodity handling and sorting operations. Broadly speaking, HLPs consider the location of a set of hubs and the design of the hub-and-spoke network so as to minimize the total flow cost.

枢纽选址问题大量发生在轴辐式网络中，它在航空运输、快递运输和快速转运系统中中应用广泛。这些系统的一个特点是在多个起点和多个目的地之间运输货物流，但是这些货物流并不是直接运输的，而是通过枢纽节点将一些起点聚合起来，之后通过枢纽间转运，最后分发到目的地节点。这样网络中的链路数据将会大大减少，从而实现系统的规模经济效应，进而避免了大量链路的建设成本。从广义上讲，HLP 考虑了一组枢纽的位置和枢纽辐射型网络的设计，以最大限度地降低总流量成本。

The first one is that flows have to be routed via hubs and thus, paths between O–D nodes must include at least one hub. Second, it is possible to connect hubs with more effective pathways that allow a constant discount factor to be applied to the flow cost between hubs. The third assumption is that hub arcs have no setup cost and thus, hub facilities can be connected at no additional cost. The last one is that distances between nodes satisfy the triangle inequality.

该类问题通常遵循如下四个假设：

• 货物流必须通过枢纽进行转运，也就是每个OD对之间至少包含一个枢纽。
• 通过枢纽将货物流更有效的连接起来，已在枢纽间实现规模经济效益
• 枢纽之间的链路没有设置成本
• 节点间的距离满足三角不等式

可知：后两个假设允许在枢纽间建立完全图；而假设1、3和4导致每个OD对之间最少一个枢纽，最多两个枢纽。

完全互通的网络总是令人望而却步，因此，我们放松了该假设。

2 问题描述与模型构建

2.1 符号描述

我们建立了一个完全图G=(N,A)，其中N={1,2,…,n}表示节点集合和潜在的枢纽点集合；A表示节点间的链路集合。对于任意一个有序点对(i,j)，存在一个从节点i到节点 j 的流量W_{ij}.

cycle hub location problem, CHLP 尝试通过一个回路连接所有枢纽及通过枢纽的链路。每个节点将分配到一个枢纽，且当它被选为枢纽时，就是自己分配给自己。目标是最小化总流成本。对于CHLP的每个可行解，有：

• 存在 p 个枢纽
• 每个枢纽正好与另外两个枢纽相连
• 枢纽派生的路径不包含子回路
• 每个点对恰好有两条路径

当 p ={1, 2, 3} 时，问题将退化为p-hub选址问题。

2.2 模型数学表达式

模型的决策变量说明：

• x_{ikm}：连续型决策变量，表示源子节点i且经过枢纽k和枢纽m的流量
• y_{km}：表示枢纽k和枢纽m之间的弧段被选中
• z_{ik}：表示非枢纽节点i分配给枢纽节点k，当i=k时，表示节点被选做枢纽

 模型解释（ps：中文解释见源码）：

•  Constraints (2) ensure that each node is assigned to exactly one hub.
• Constraint (3) is a cardinality constraint on the number of hubs that must be opened
• constraint (4) state that the number of hub arcs in the cycle is equal to p.
• Constraints (5) guarantee that each hub node must be connected to exactly two other hub nodes.
• Constraints (6) are the flow conservation constraints.
• Constraints (7) and (8) ensure that both end nodes of a hub arc are opened hubs and also, they ensure that non-hub nodes are assigned to an open hub.
• Constraints (9) state that the flow between hubs moves through the hub cycle.
• constraints (10–12) are the standard nonnegativity and integrality constraints.

3 模型验证及代码实现

3.1 模型验证源代码

import matplotlib.pyplot as plt
from gurobipy import GRB, Model, quicksum
import numpy as np
import time

seedt = time.time()
rnd = np.random
rnd.seed(int(seedt))
# 1. 数据初始化
nodeNum = 20
alpha = 0.6
p = 4

N = range(0,nodeNum)
A = [(i,j) for i in range(0,nodeNum) for j in range(0,nodeNum)]
O = np.zeros(nodeNum)
# 定义节点间距离矩阵
pts = rnd.rand(nodeNum, 2) * 200
cost = {(i, j): np.hypot(pts[i, 0] - pts[j, 0], pts[i, 1] - pts[j, 1]) for i, j in A}  # 计算距离向量
# 定义节点间流量矩阵
w = {(i,j): np.random.randint(30,50) for i,j in A}

# 从节点i流出的流量
O = [i for i in N]
D = [i for i in N]
for i in N:
    sumO = 0
    for j in N:
        sumO = sumO+w[i,j]
    O[i] = sumO
print(f'源流量汇总：{O}')
for j in N:
    sumD = 0
    for i in N:
        sumD = sumD + w[i,j]
    D[j] = sumD
print(f'汇流量汇总：{D}')



# 定义决策变量需要的集合
IK = [(i,k) for i in N for k in N]
IKM = [(i,k,m) for i in N for k in N for m in N]
KM = [(k,m) for k in N for m in N]

chlp = Model('CHLP')
z = chlp.addVars(IK, vtype='B', name='z') # 二进制决策变量，将节点i分配到枢纽k
x = chlp.addVars(IKM, lb=0, vtype='C', name='x') # 连续型决策变量，起点为i的流经过枢纽k和枢纽m的流量
y = chlp.addVars(KM, vtype='B', name='y') # 二进制行决策变量，枢纽k和枢纽m被选中

chlp.modelSense = GRB.MINIMIZE
chlp.setObjective(quicksum((cost[i,k]*O[i] + cost[k,i]*D[i])*z[i,k] for i,k in IK)
                  + quicksum(cost[k,m]*alpha*x[i,k,m] for i,k,m in IKM))

# 约束：每个节点只分到一个枢纽节点
chlp.addConstrs(quicksum(z[i,k] for k in N) == 1 for i in N)
# 约束：选出p个枢纽
chlp.addConstr(quicksum(z[k,k] for k in N) == p)
# 约束：回路中的枢纽枢纽等于p
chlp.addConstr(quicksum(y[k,m] for k in N for m in N) == p)
# 约束：每个枢纽节点精确连接其他两个枢纽
chlp.addConstrs(quicksum(y[k,m] for m in N if m>k) + quicksum(y[m,k] for m in N if k>m)
                == 2*z[k,k] for k in N)
# 约束：流平衡约束
chlp.addConstrs(O[i]*z[i,k] + quicksum(x[i,m,k] for m in N if m != k)
                == quicksum(x[i,k,m] for m in N if m!=k) + quicksum(w[i,m]*z[m,k] for m in N)
                for i in N for k in N if k != i)

# 约束：开放枢纽的两端必定是枢纽
chlp.addConstrs(z[k,m] + y[k,m] <= z[m,m] for k,m in KM if m>k)
chlp.addConstrs(z[m,k] + y[k,m] <= z[k,k] for k,m in KM if m>k)
# 约束：枢纽之间的流在回路上进行
chlp.addConstrs(x[i,k,m] + x[i,m,k] <= O[i]*y[k,m] for i,k,m in IKM if m>k)

chlp.optimize()


# for i, m, k in IKM:
#     if x[i, m, k].x > 0:
#         plt.plot([pts[i, 0], pts[m, 0]], [pts[i, 1], pts[m, 1]], c='g')
#         plt.plot([pts[m, 0], pts[k, 0]], [pts[m, 1], pts[k, 1]], c='b')

# 显示hub和spoke之间的分配连接
for i, k in IK:
    if z[i, k].x > 0.99:
        plt.plot([pts[i, 0], pts[k, 0]], [pts[i, 1], pts[k, 1]], c='g')
# 显示hub之间的连接回路
for k,m in KM:
    if y[k,m].x > 0.99:
        plt.plot([pts[m, 0], pts[k, 0]], [pts[m, 1], pts[k, 1]], c='b')
# 为每个点标注数字记号
for i, txt in enumerate(N):
    plt.scatter(pts[i, 0], pts[i, 1], c='y', marker='s')
    plt.annotate(txt, (pts[i, 0] + 2, pts[i, 1] + 2))
# 可视化绘制枢纽节点
for k in N:
    if z[k,k].x >= 0.99:
        plt.scatter(pts[k, 0], pts[k, 1], c='r', marker='s')
plt.show()  # 最后必须要显示图片

3.2 验证结果可视化