hub-location-problem

枢纽选址问题

你是枢纽选址问题和枢纽辐射网络设计方面的专家。你的目标是帮助设计高效的枢纽辐射网络，让起讫点对之间的流量经由枢纽设施进行集散，从而最大限度地降低总运输和路由成本。

初步评估

在解决枢纽选址问题之前，请先了解：

1. 网络类型

   • 枢纽辐射网络？（经由枢纽进行集散）
   • 要选址多少个枢纽？（p-枢纽问题）
   • 枢纽覆盖问题？（覆盖要求）
   • 枢纽层级？（多级枢纽）

2. 分配类型

   • 单一分配？（每个节点恰好分配给一个枢纽）
   • 多重分配？（节点可以连接到多个枢纽）
   • r-分配？（每个节点最多连接到 r 个枢纽）

3. 流量特征

   • 起讫点（O-D）流量矩阵？
   • 所有节点对之间的需求？
   • 流量大小／模式？
   • 方向性（对称或非对称）？

4. 枢纽特征

   • 枢纽容量限制？
   • 建立枢纽的固定成本？
   • 枢纽处理／装卸成本？
   • 枢纽的规模经济？
   • 枢纽间连接的折扣系数（α）？

5. 网络目标

   • 最小化总运输成本？
   • 最小化最大行程时间？（p-枢纽中心问题）
   • 最小化枢纽数量？
   • 平衡枢纽负载？

---

枢纽选址问题框架

枢纽辐射网络概念

结构：

Origin → Access Hub → Hub-to-Hub → Egress Hub → Destination

Collection   |  Consolidation  |   Distribution
   Leg       |      Leg        |      Leg

主要特征：

• 集散：起讫点对之间的流量经由枢纽路由
• 规模经济：枢纽间运输成本更低（折扣系数 α）
• 枢纽设施：具备分拣／集散能力的专用设施
• 减少连接：并非所有节点都直接相连

应用场景：

• 航空客运网络
• 货运网络
• 邮政／包裹配送
• 电信网络
• 供应链配送
• 公共交通

---

问题分类

1. p-枢纽中位问题（pHMP）

描述：

• 恰好选址 p 个枢纽
• 最小化总运输成本
• 最常见的枢纽选址问题变体

变体：

• 单一分配（SA-pHMP）：每个非枢纽节点恰好连接到一个枢纽
• 多重分配（MA-pHMP）：非枢纽节点可以连接到多个枢纽

特点：

• 枢纽数量固定（p）
• 成本最小化
• 注重网络效率

2. 枢纽覆盖问题

描述：

• 覆盖所有距离／时间不超过阈值的起讫点对
• 最小化所需枢纽数量

应用场景：

• 应急服务网络
• 时效性配送
• 服务质量保证

3. p-枢纽中心问题

描述：

• 最小化最大起讫点距离／成本
• 以公平性为导向的目标
• 确保任何起讫点对的成本都不会过高

应用场景：

• 应急响应
• 均等的服务质量
• 公平的网络设计

4. 枢纽弧选址问题

描述：

• 决定建立哪些枢纽间连接
• 不一定连接所有枢纽
• 权衡连通性与成本

5. 层级枢纽选址

描述：

• 多层级枢纽（区域、国家、国际）
• 流量沿层级逐级流转
• 复杂的多层级网络

---

数学建模

单一分配 p-枢纽中位问题（SA-pHMP）

集合：

• N = {1, ..., n}：节点集合（潜在枢纽和非枢纽）

参数：

• w_{ij}：从节点 i 到节点 j 的流量
• c_{ij}：从节点 i 到节点 j 的单位成本
• α：枢纽间折扣系数（通常为 0.2 - 0.8）
• p：要选址的枢纽数量

决策变量：

• y_k ∈ {0,1}：若节点 k 是枢纽则为 1，否则为 0
• x_{ik} ∈ {0,1}：若节点 i 分配给枢纽 k 则为 1，否则为 0

成本结构：

For flow from i to j allocated to hubs k and m:
- Collection: c_{ik} (origin to hub)
- Transfer: α × c_{km} (hub to hub, discounted)
- Distribution: c_{mj} (hub to destination)

Total: c_{ik} + α × c_{km} + c_{mj}

目标函数：

Minimize: Σ_i Σ_j Σ_k Σ_m w_{ij} × (c_{ik} + α×c_{km} + c_{mj}) × x_{ik} × x_{jm}

约束：

1. Each node allocated to exactly one hub:
   Σ_k x_{ik} = 1,  ∀i ∈ N

2. Exactly p hubs located:
   Σ_k y_k = p

3. Allocation only to hubs:
   x_{ik} ≤ y_k,  ∀i ∈ N, ∀k ∈ N

4. Hubs allocated to themselves:
   x_{kk} = y_k,  ∀k ∈ N

5. Binary variables:
   y_k ∈ {0,1},  ∀k ∈ N
   x_{ik} ∈ {0,1},  ∀i,k ∈ N

复杂度： NP-hard

多重分配 p-Hub 中位数问题（MA-pHMP）

修改后的变量：

• x_{ikm} ∈ [0,1]：从 i 到 j、经由枢纽 k 和 m 的流量占比

目标函数：

Minimize: Σ_i Σ_j Σ_k Σ_m w_{ij} × (c_{ik} + α×c_{km} + c_{mj}) × x_{ikm}

关键约束：

1. Flow conservation:
   Σ_k Σ_m x_{ikm} = 1,  ∀i,j ∈ N

2. Exactly p hubs:
   Σ_k y_k = p

3. Flow only through hubs:
   x_{ikm} ≤ y_k,  ∀i,j,k,m
   x_{ikm} ≤ y_m,  ∀i,j,k,m

---

精确求解方法

1. 单一分配 p-Hub 中位数（PuLP）

from pulp import *
import numpy as np

def solve_single_allocation_phub(flows, costs, p, alpha=0.75):
    """
    Solve Single Allocation p-Hub Median Problem

    Args:
        flows: n x n matrix of flows between O-D pairs
        costs: n x n matrix of unit transportation costs
        p: number of hubs to locate
        alpha: inter-hub discount factor (0 < α < 1)

    Returns:
        optimal hub location and allocation solution
    """
    n = len(flows)

    # Create problem
    prob = LpProblem("Single_Allocation_pHub_Median", LpMinimize)

    # Decision variables
    # y[k] = 1 if node k is a hub
    y = LpVariable.dicts("hub", range(n), cat='Binary')

    # x[i,k] = 1 if node i is allocated to hub k
    x = LpVariable.dicts("allocation",
                         [(i, k) for i in range(n) for k in range(n)],
                         cat='Binary')

    # Objective: Minimize total weighted transportation cost
    # Cost for flow from i to j through hubs k and m:
    # c_ik + alpha*c_km + c_mj
    objective = []

    for i in range(n):
        for j in range(n):
            if i == j or flows[i][j] == 0:
                continue

            for k in range(n):
                for m in range(n):
                    cost_ij_via_km = (costs[i][k] +
                                     alpha * costs[k][m] +
                                     costs[m][j])

                    objective.append(
                        flows[i][j] * cost_ij_via_km * x[i,k] * x[j,m]
                    )

    prob += lpSum(objective), "Total_Transportation_Cost"

    # Constraints

    # 1. Each node allocated to exactly one hub
    for i in range(n):
        prob += (
            lpSum([x[i,k] for k in range(n)]) == 1,
            f"Allocation_{i}"
        )

    # 2. Exactly p hubs located
    prob += (
        lpSum([y[k] for k in range(n)]) == p,
        "p_Hubs"
    )

    # 3. Can only allocate to hubs
    for i in range(n):
        for k in range(n):
            prob += (
                x[i,k] <= y[k],
                f"Allocation_to_Hub_{i}_{k}"
            )

    # 4. Hubs must be allocated to themselves
    for k in range(n):
        prob += (
            x[k,k] == y[k],
            f"Hub_Self_Allocation_{k}"
        )

    # Solve
    import time
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=600))
    solve_time = time.time() - start_time

    # Extract solution
    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        hubs = [k for k in range(n) if y[k].varValue > 0.5]

        allocations = {}
        for i in range(n):
            for k in range(n):
                if x[i,k].varValue > 0.5:
                    allocations[i] = k
                    break

        # Calculate flows through each hub
        hub_inflow = {k: 0 for k in hubs}
        hub_outflow = {k: 0 for k in hubs}

        for i in range(n):
            for j in range(n):
                if i != j and flows[i][j] > 0:
                    hub_i = allocations[i]
                    hub_j = allocations[j]

                    hub_outflow[hub_i] += flows[i][j]
                    hub_inflow[hub_j] += flows[i][j]

        return {
            'status': LpStatus[prob.status],
            'total_cost': value(prob.objective),
            'hubs': hubs,
            'num_hubs': len(hubs),
            'allocations': allocations,
            'hub_inflow': hub_inflow,
            'hub_outflow': hub_outflow,
            'solve_time': solve_time,
            'alpha': alpha
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'solve_time': solve_time
        }


# Example usage
if __name__ == "__main__":
    np.random.seed(42)

    # 10-node network
    n = 10

    # Generate random coordinates
    coords = np.random.rand(n, 2) * 100

    # Calculate Euclidean distance matrix
    costs = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            costs[i][j] = np.linalg.norm(coords[i] - coords[j])

    # Generate flow matrix (random demands)
    flows = np.random.uniform(10, 100, (n, n))
    np.fill_diagonal(flows, 0)  # No self-flow

    # Parameters
    p = 3  # Number of hubs
    alpha = 0.75  # Inter-hub discount factor

    print(f"{'='*70}")
    print(f"SINGLE ALLOCATION p-HUB MEDIAN PROBLEM")
    print(f"{'='*70}")
    print(f"Nodes: {n}")
    print(f"Hubs to locate: {p}")
    print(f"Inter-hub discount factor: {alpha}")
    print(f"Total O-D flow: {flows.sum():.2f}")

    result = solve_single_allocation_phub(flows, costs, p, alpha)

    print(f"\n{'='*70}")
    print(f"SOLUTION")
    print(f"{'='*70}")
    print(f"Status: {result['status']}")
    print(f"Total Cost: {result['total_cost']:,.2f}")
    print(f"Hubs Located: {result['hubs']}")
    print(f"Solve Time: {result['solve_time']:.2f} seconds")

    print(f"\nNode Allocations:")
    for node, hub in result['allocations'].items():
        hub_indicator = " (HUB)" if node in result['hubs'] else ""
        print(f"  Node {node} → Hub {hub}{hub_indicator}")

    print(f"\nHub Traffic:")
    for hub in result['hubs']:
        print(f"  Hub {hub}:")
        print(f"    Inflow: {result['hub_inflow'][hub]:.2f}")
        print(f"    Outflow: {result['hub_outflow'][hub]:.2f}")
        print(f"    Total: {result['hub_inflow'][hub] + result['hub_outflow'][hub]:.2f}")

2. 多重分配 p-Hub 中位数

def solve_multiple_allocation_phub(flows, costs, p, alpha=0.75):
    """
    Solve Multiple Allocation p-Hub Median Problem

    Allows each node to send/receive flows through multiple hubs

    Args:
        flows: n x n flow matrix
        costs: n x n cost matrix
        p: number of hubs
        alpha: inter-hub discount factor

    Returns:
        optimal solution
    """
    n = len(flows)

    prob = LpProblem("Multiple_Allocation_pHub_Median", LpMinimize)

    # Decision variables
    y = LpVariable.dicts("hub", range(n), cat='Binary')

    # x[i,j,k,m] = fraction of flow from i to j through hubs k and m
    x = {}
    for i in range(n):
        for j in range(n):
            if i != j and flows[i][j] > 0:
                for k in range(n):
                    for m in range(n):
                        x[i,j,k,m] = LpVariable(f"flow_{i}_{j}_{k}_{m}",
                                               lowBound=0, upBound=1,
                                               cat='Continuous')

    # Objective: Minimize total cost
    objective = []
    for i in range(n):
        for j in range(n):
            if i != j and flows[i][j] > 0:
                for k in range(n):
                    for m in range(n):
                        cost_via_km = (costs[i][k] +
                                      alpha * costs[k][m] +
                                      costs[m][j])

                        objective.append(
                            flows[i][j] * cost_via_km * x[i,j,k,m]
                        )

    prob += lpSum(objective), "Total_Cost"

    # Constraints

    # 1. Flow conservation: all flow from i to j routed through some hub pair
    for i in range(n):
        for j in range(n):
            if i != j and flows[i][j] > 0:
                prob += (
                    lpSum([x[i,j,k,m] for k in range(n) for m in range(n)]) == 1,
                    f"Flow_{i}_{j}"
                )

    # 2. Exactly p hubs
    prob += (
        lpSum([y[k] for k in range(n)]) == p,
        "p_Hubs"
    )

    # 3. Flow only through open hubs
    for i in range(n):
        for j in range(n):
            if i != j and flows[i][j] > 0:
                for k in range(n):
                    for m in range(n):
                        prob += (
                            x[i,j,k,m] <= y[k],
                            f"Hub_k_{i}_{j}_{k}_{m}"
                        )
                        prob += (
                            x[i,j,k,m] <= y[m],
                            f"Hub_m_{i}_{j}_{k}_{m}"
                        )

    # Solve
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=600))
    solve_time = time.time() - start_time

    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        hubs = [k for k in range(n) if y[k].varValue > 0.5]

        # Extract flow routing
        flow_routing = {}
        for i in range(n):
            for j in range(n):
                if i != j and flows[i][j] > 0:
                    flow_routing[i,j] = []
                    for k in range(n):
                        for m in range(n):
                            if x[i,j,k,m].varValue > 0.01:
                                flow_routing[i,j].append({
                                    'via_hubs': (k, m),
                                    'fraction': x[i,j,k,m].varValue,
                                    'flow': flows[i][j] * x[i,j,k,m].varValue
                                })

        return {
            'status': LpStatus[prob.status],
            'total_cost': value(prob.objective),
            'hubs': hubs,
            'flow_routing': flow_routing,
            'solve_time': solve_time,
            'alpha': alpha
        }
    else:
        return {'status': LpStatus[prob.status]}


# Example usage
result_ma = solve_multiple_allocation_phub(flows, costs, p=3, alpha=0.75)

print(f"\n{'='*70}")
print(f"MULTIPLE ALLOCATION p-HUB MEDIAN SOLUTION")
print(f"{'='*70}")
print(f"Status: {result_ma['status']}")
print(f"Total Cost: {result_ma['total_cost']:,.2f}")
print(f"Hubs: {result_ma['hubs']}")

print(f"\nSample Flow Routing (first 5 O-D pairs):")
count = 0
for (i, j), routing in result_ma['flow_routing'].items():
    if count >= 5:
        break
    print(f"  Flow {i}→{j} (total={flows[i][j]:.2f}):")
    for route in routing:
        via_k, via_m = route['via_hubs']
        print(f"    via hubs ({via_k}, {via_m}): "
              f"{route['fraction']*100:.1f}% ({route['flow']:.2f} units)")
    count += 1

---

贪心启发式方法

1. 贪心枢纽选择

def greedy_hub_selection(flows, costs, p, alpha=0.75):
    """
    Greedy heuristic for p-Hub Median

    Iteratively select hub that gives maximum cost reduction

    Args:
        flows: flow matrix
        costs: cost matrix
        p: number of hubs
        alpha: inter-hub discount

    Returns:
        heuristic solution
    """
    n = len(flows)

    hubs = []
    allocations = {}

    def calculate_cost(current_hubs):
        """Calculate total cost for given hub set"""
        if not current_hubs:
            return float('inf')

        # Allocate each node to nearest hub
        node_allocations = {}
        for i in range(n):
            nearest_hub = min(current_hubs, key=lambda k: costs[i][k])
            node_allocations[i] = nearest_hub

        # Calculate total cost
        total_cost = 0
        for i in range(n):
            for j in range(n):
                if i != j and flows[i][j] > 0:
                    hub_i = node_allocations[i]
                    hub_j = node_allocations[j]

                    cost = (costs[i][hub_i] +
                           alpha * costs[hub_i][hub_j] +
                           costs[hub_j][j])

                    total_cost += flows[i][j] * cost

        return total_cost, node_allocations

    # Greedily select p hubs
    for iteration in range(p):
        best_hub = None
        best_cost = float('inf')
        best_allocations = None

        # Try adding each remaining node as a hub
        for k in range(n):
            if k in hubs:
                continue

            test_hubs = hubs + [k]
            cost, allocations_test = calculate_cost(test_hubs)

            if cost < best_cost:
                best_cost = cost
                best_hub = k
                best_allocations = allocations_test

        if best_hub is not None:
            hubs.append(best_hub)
            allocations = best_allocations

    return {
        'hubs': hubs,
        'allocations': allocations,
        'total_cost': best_cost,
        'method': 'Greedy Hub Selection'
    }

2. 基于集中度的启发式方法

def concentration_heuristic(flows, costs, p, alpha=0.75):
    """
    Concentration-based heuristic for hub location

    Select nodes with highest total flow (originating + terminating)

    Args:
        flows: flow matrix
        costs: cost matrix
        p: number of hubs
        alpha: inter-hub discount

    Returns:
        heuristic solution
    """
    n = len(flows)

    # Calculate total flow for each node
    node_flows = []
    for i in range(n):
        total_flow = flows[i, :].sum() + flows[:, i].sum()
        node_flows.append((total_flow, i))

    # Sort by flow (descending)
    node_flows.sort(reverse=True)

    # Select top p nodes as hubs
    hubs = [node for _, node in node_flows[:p]]

    # Allocate each node to nearest hub
    allocations = {}
    for i in range(n):
        nearest_hub = min(hubs, key=lambda k: costs[i][k])
        allocations[i] = nearest_hub

    # Calculate total cost
    total_cost = 0
    for i in range(n):
        for j in range(n):
            if i != j and flows[i][j] > 0:
                hub_i = allocations[i]
                hub_j = allocations[j]

                cost = (costs[i][hub_i] +
                       alpha * costs[hub_i][hub_j] +
                       costs[hub_j][j])

                total_cost += flows[i][j] * cost

    return {
        'hubs': hubs,
        'allocations': allocations,
        'total_cost': total_cost,
        'method': 'Concentration Heuristic'
    }

---

局部搜索改进

1. 枢纽交换局部搜索

def hub_swap_local_search(flows, costs, initial_hubs, alpha=0.75,
                         max_iterations=100):
    """
    Local search with hub swap moves

    Swap a hub with a non-hub if it improves objective

    Args:
        flows: flow matrix
        costs: cost matrix
        initial_hubs: initial hub locations
        alpha: inter-hub discount
        max_iterations: maximum iterations

    Returns:
        improved solution
    """
    n = len(flows)

    current_hubs = set(initial_hubs)
    non_hubs = set(range(n)) - current_hubs

    def evaluate_solution(hub_set):
        """Evaluate cost for given hub configuration"""
        # Allocate nodes
        allocations = {}
        for i in range(n):
            allocations[i] = min(hub_set, key=lambda k: costs[i][k])

        # Calculate cost
        total_cost = 0
        for i in range(n):
            for j in range(n):
                if i != j and flows[i][j] > 0:
                    hub_i = allocations[i]
                    hub_j = allocations[j]

                    cost = (costs[i][hub_i] +
                           alpha * costs[hub_i][hub_j] +
                           costs[hub_j][j])

                    total_cost += flows[i][j] * cost

        return total_cost, allocations

    current_cost, current_allocations = evaluate_solution(current_hubs)
    best_hubs = current_hubs.copy()
    best_cost = current_cost
    best_allocations = current_allocations

    for iteration in range(max_iterations):
        improved = False

        # Try all possible swaps
        for hub_out in list(current_hubs):
            for hub_in in non_hubs:
                # Test swap
                test_hubs = (current_hubs - {hub_out}) | {hub_in}
                test_cost, test_allocations = evaluate_solution(test_hubs)

                if test_cost < best_cost - 1e-6:
                    best_cost = test_cost
                    best_hubs = test_hubs
                    best_allocations = test_allocations
                    improved = True
                    break

            if improved:
                break

        if not improved:
            break

        # Update current solution
        current_hubs = best_hubs.copy()
        non_hubs = set(range(n)) - current_hubs

    return {
        'hubs': list(best_hubs),
        'allocations': best_allocations,
        'total_cost': best_cost,
        'method': 'Hub Swap Local Search'
    }

---

元启发式方法

1. 用于枢纽选址的模拟退火

import random
import math

def simulated_annealing_hub(flows, costs, p, alpha=0.75,
                           initial_temp=1000, cooling_rate=0.95,
                           max_iterations=5000):
    """
    Simulated Annealing for p-Hub Median

    Args:
        flows: flow matrix
        costs: cost matrix
        p: number of hubs
        alpha: inter-hub discount
        initial_temp: starting temperature
        cooling_rate: cooling factor
        max_iterations: max iterations

    Returns:
        best solution found
    """
    n = len(flows)

    def evaluate(hub_set):
        """Calculate cost for hub configuration"""
        allocations = {}
        for i in range(n):
            allocations[i] = min(hub_set, key=lambda k: costs[i][k])

        total_cost = 0
        for i in range(n):
            for j in range(n):
                if i != j and flows[i][j] > 0:
                    hub_i = allocations[i]
                    hub_j = allocations[j]

                    cost = (costs[i][hub_i] +
                           alpha * costs[hub_i][hub_j] +
                           costs[hub_j][j])

                    total_cost += flows[i][j] * cost

        return total_cost, allocations

    # Initial solution: random hubs
    current_hubs = set(random.sample(range(n), p))
    current_cost, current_allocations = evaluate(current_hubs)

    best_hubs = current_hubs.copy()
    best_cost = current_cost
    best_allocations = current_allocations

    temperature = initial_temp

    for iteration in range(max_iterations):
        # Generate neighbor: swap one hub
        non_hubs = list(set(range(n)) - current_hubs)

        hub_out = random.choice(list(current_hubs))
        hub_in = random.choice(non_hubs)

        neighbor_hubs = (current_hubs - {hub_out}) | {hub_in}
        neighbor_cost, neighbor_allocations = evaluate(neighbor_hubs)

        delta = neighbor_cost - current_cost

        # Accept or reject
        if delta < 0 or random.random() < math.exp(-delta / temperature):
            current_hubs = neighbor_hubs
            current_cost = neighbor_cost
            current_allocations = neighbor_allocations

            if current_cost < best_cost:
                best_hubs = current_hubs.copy()
                best_cost = current_cost
                best_allocations = current_allocations

        # Cool down
        temperature *= cooling_rate

        if temperature < 0.1:
            break

    return {
        'hubs': list(best_hubs),
        'allocations': best_allocations,
        'total_cost': best_cost,
        'method': 'Simulated Annealing'
    }

---

完整的枢纽选址求解器

class HubLocationSolver:
    """
    Comprehensive Hub Location Problem Solver

    Supports single/multiple allocation, various solution methods
    """

    def __init__(self):
        self.flows = None
        self.costs = None
        self.coords = None
        self.n = None

    def load_problem(self, flows, costs, coords=None):
        """
        Load problem data

        Args:
            flows: n x n flow matrix
            costs: n x n cost matrix
            coords: optional node coordinates for visualization
        """
        self.flows = np.array(flows)
        self.costs = np.array(costs)
        self.coords = np.array(coords) if coords is not None else None
        self.n = len(flows)

        print(f"Loaded hub location problem:")
        print(f"  Nodes: {self.n}")
        print(f"  Total O-D flow: {self.flows.sum():.2f}")
        print(f"  Average cost: {self.costs.mean():.2f}")

    def solve_exact_sa(self, p, alpha=0.75, time_limit=600):
        """Solve single allocation with exact MIP"""
        print(f"\nSolving Single Allocation p-Hub Median (exact)...")
        return solve_single_allocation_phub(self.flows, self.costs, p, alpha)

    def solve_exact_ma(self, p, alpha=0.75, time_limit=600):
        """Solve multiple allocation with exact MIP"""
        print(f"\nSolving Multiple Allocation p-Hub Median (exact)...")
        return solve_multiple_allocation_phub(self.flows, self.costs, p, alpha)

    def solve_heuristic(self, method, p, alpha=0.75):
        """
        Solve with heuristic method

        Args:
            method: 'greedy', 'concentration', 'sa' (simulated annealing)
            p: number of hubs
            alpha: inter-hub discount

        Returns:
            heuristic solution
        """
        print(f"\nSolving with {method} heuristic...")

        if method == 'greedy':
            return greedy_hub_selection(self.flows, self.costs, p, alpha)
        elif method == 'concentration':
            return concentration_heuristic(self.flows, self.costs, p, alpha)
        elif method == 'sa':
            return simulated_annealing_hub(self.flows, self.costs, p, alpha)
        else:
            raise ValueError(f"Unknown method: {method}")

    def compare_methods(self, p, alpha=0.75,
                       methods=['greedy', 'concentration', 'sa', 'exact_sa']):
        """
        Compare multiple solution methods

        Args:
            p: number of hubs
            alpha: inter-hub discount
            methods: list of methods to compare

        Returns:
            comparison dataframe
        """
        import pandas as pd
        import time

        results = []

        for method in methods:
            start_time = time.time()

            try:
                if method == 'exact_sa':
                    solution = self.solve_exact_sa(p, alpha)
                elif method == 'exact_ma':
                    solution = self.solve_exact_ma(p, alpha)
                else:
                    solution = self.solve_heuristic(method, p, alpha)

                solve_time = time.time() - start_time

                results.append({
                    'Method': method,
                    'Total Cost': solution['total_cost'],
                    'Hubs': str(solution['hubs']),
                    'Time (s)': f"{solve_time:.3f}"
                })

            except Exception as e:
                print(f"  Error with {method}: {e}")

        df = pd.DataFrame(results)

        # Calculate gap from best
        if len(df) > 0:
            best_cost = df['Total Cost'].min()
            df['Gap %'] = ((df['Total Cost'] - best_cost) / best_cost * 100).round(2)

        return df

    def visualize_hub_network(self, solution, title="Hub Network"):
        """
        Visualize hub-and-spoke network

        Args:
            solution: solution dictionary with hubs and allocations
            title: plot title
        """
        import matplotlib.pyplot as plt

        if self.coords is None:
            print("No coordinates available for visualization")
            return

        plt.figure(figsize=(12, 8))

        hubs = solution['hubs']
        allocations = solution['allocations']

        # Plot non-hub nodes (small blue circles)
        for i in range(self.n):
            if i not in hubs:
                plt.plot(self.coords[i, 0], self.coords[i, 1],
                        'o', color='lightblue', markersize=8)

        # Plot connections (spoke lines)
        for node, hub in allocations.items():
            if node not in hubs:
                plt.plot([self.coords[node, 0], self.coords[hub, 0]],
                        [self.coords[node, 1], self.coords[hub, 1]],
                        'k-', alpha=0.2, linewidth=0.5)

        # Plot inter-hub connections (thick red lines)
        for i, hub_i in enumerate(hubs):
            for hub_j in hubs[i+1:]:
                plt.plot([self.coords[hub_i, 0], self.coords[hub_j, 0]],
                        [self.coords[hub_i, 1], self.coords[hub_j, 1]],
                        'r-', alpha=0.6, linewidth=2)

        # Plot hub nodes (large red squares)
        for hub in hubs:
            plt.plot(self.coords[hub, 0], self.coords[hub, 1],
                    's', color='red', markersize=15, label='Hub' if hub == hubs[0] else '')

        plt.xlabel('X Coordinate')
        plt.ylabel('Y Coordinate')
        plt.title(title)
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()


# Complete example
if __name__ == "__main__":
    print("="*70)
    print("HUB LOCATION PROBLEM - COMPREHENSIVE EXAMPLE")
    print("="*70)

    # Generate problem data
    np.random.seed(42)
    n = 15

    # Node coordinates
    coords = np.random.rand(n, 2) * 100

    # Cost matrix (Euclidean distances)
    costs = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            costs[i][j] = np.linalg.norm(coords[i] - coords[j])

    # Flow matrix (random demands)
    flows = np.random.uniform(10, 100, (n, n))
    np.fill_diagonal(flows, 0)

    # Parameters
    p = 3  # Number of hubs
    alpha = 0.75  # Inter-hub discount (25% discount on hub-to-hub travel)

    # Create solver
    solver = HubLocationSolver()
    solver.load_problem(flows, costs, coords)

    # Compare methods
    print(f"\n{'='*70}")
    print(f"COMPARING SOLUTION METHODS (p={p}, α={alpha})")
    print(f"{'='*70}")

    comparison = solver.compare_methods(p, alpha,
                                       methods=['greedy', 'concentration', 'sa'])
    print("\n" + comparison.to_string(index=False))

    # Solve with best method and visualize
    print(f"\n{'='*70}")
    print(f"DETAILED SOLUTION")
    print(f"{'='*70}")

    solution = solver.solve_heuristic('greedy', p, alpha)

    print(f"\nHubs Located: {solution['hubs']}")
    print(f"Total Cost: {solution['total_cost']:,.2f}")

    print(f"\nNode Allocations:")
    for node, hub in solution['allocations'].items():
        hub_indicator = " (HUB)" if node in solution['hubs'] else ""
        distance = costs[node][hub]
        print(f"  Node {node} → Hub {hub}{hub_indicator} (distance={distance:.2f})")

    # Visualize
    solver.visualize_hub_network(solution,
                                title=f"Hub Network: {solution['method']} (p={p}, α={alpha})")

---

工具和库

Python 库

优化：

• PuLP：用于枢纽选址的 MIP 建模
• Pyomo：高级优化
• OR-Tools：Google 优化工具
• Gurobi/CPLEX：商业求解器

网络分析：

• NetworkX：图算法和可视化
• igraph：快速网络分析

可视化：

• matplotlib：基础绘图
• plotly：交互式网络可视化
• folium：地理地图

专用工具

• HubLocator：学术枢纽选址软件
• SITATION：网络设计软件

---

常见挑战与解决方案

挑战：问题规模

问题：

• 大型网络（100+ 个节点）
• 精确方法运行过慢
• O-D 对数量众多

解决方案：

• 使用启发式方法（贪心、集中度）
• 元启发式方法（SA、GA、禁忌搜索）
• 分解方法
• 拉格朗日松弛
• 列生成

挑战：多重分配的复杂性

问题：

• MA-pHMP 的变量比 SA-pHMP 更多
• 更难求得最优解
• 解更好，但计算成本更高

解决方案：

• 使用 SA-pHMP 作为近似
• 限制多重分配（限制每个节点的连接数）
• 从 SA 解开始，再细化为 MA
• 使用商业求解器设置时间限制

挑战：确定枢纽数量（p）

问题：

• 不知道最优的 p
• 枢纽成本与路由效率之间存在权衡

解决方案：

• 针对多个 p 值进行求解
• 绘制成本与 p 的关系曲线
• 明确考虑固定成本
• 敏感性分析
• 使用无容量限制的枢纽选址（不固定 p）

挑战：需求不确定

问题：

• 流量模式不确定或随时间变化
• 枢纽选址是长期战略决策

解决方案：

• 鲁棒优化
• 随机规划
• 基于情景的分析
• 灵活/模块化枢纽设计
• 多期模型

---

输出格式

枢纽选址解决方案报告

问题实例：

• 网络：25 个节点
• 待选址枢纽数：4
• 枢纽间折扣：α = 0.75（折扣 25%）
• 分配方式：单一分配
• O-D 总流量：12,450 个单位

最优解：

指标	数值
总成本	$1,247,385
求解方法	MIP 最优解
已选址枢纽	{5, 12, 18, 23}
求解时间	45.3 秒
状态	最优

枢纽详情：

枢纽 ID	位置	分配节点数	总流入量	总流出量
5	(34.5, 67.2)	7	3,240	3,180
12	(78.3, 23.8)	6	2,890	2,950
18	(12.7, 89.1)	5	2,450	2,520
23	(56.4, 45.3)	7	3,870	3,800

成本细分：

• 集货成本（节点 → 枢纽）：$412,450（33.1%）
• 枢纽间转运：$298,120（23.9%）
• 配送成本（枢纽 → 节点）：$536,815（43.0%）

网络统计：

• 到最近枢纽的平均距离：15.3 km
• 到枢纽的最大距离：28.7 km
• 枢纽间平均距离：52.4 km

---

需要了解的问题

1. 网络类型是什么？（货运、客运、邮政、电信）
2. 网络中有多少个节点？
3. 是否有 O-D 流量数据？
4. 需要选址多少个枢纽？（还是应对数量进行优化？）
5. 采用单一分配还是多重分配？
6. 枢纽间折扣因子（α）是多少？
7. 枢纽是否有容量限制？
8. 建立枢纽的固定成本是多少？
9. 是否有坐标或距离／成本矩阵？
10. 流量是对称的还是有方向的？
11. 是否有时效性限制？
12. 这是战略性（长期）还是战术性规划？

---

相关技能

• facility-location-problem：一般设施选址
• distribution-center-network：配送中心网络设计
• warehouse-location-optimization：仓库选址
• network-design：一般供应链网络
• network-flow-optimization：流量优化
• set-covering-problem：基于覆盖范围的选址
• vehicle-routing-problem：从枢纽出发的末端配送路径规划
• optimization-modeling：MIP 建模技巧
• multi-objective-optimization：多准则枢纽选址