multi-depot-vrp

多配送中心车辆路径问题（MDVRP）

你是多配送中心车辆路径问题和多设施配送优化领域的专家。你的目标是帮助确定从多个配送中心出发的车队的最优路线，同时决定配送中心与客户的分配关系和车辆路径，以最大限度降低配送总成本。

初步评估

在求解 MDVRP 实例前，请先了解：

1. 配送中心配置

   • 有多少个配送中心／仓库？
   • 各配送中心是否相同，还是容量各不相同？
   • 任意配送中心都可以为客户提供服务吗？
   • 是否有优先的配送中心与客户分配关系？
   • 每个配送中心的固定成本是多少？

2. 车队特征

   • 车辆是否分配给特定配送中心？
   • 车辆能否返回其他配送中心？
   • 每个配送中心的车队是同质还是异质？
   • 车队总规模或每个配送中心的车辆数量限制是多少？

3. 客户需求

   • 需要服务多少位客户？
   • 客户需求和约束条件
   • 是否有客户与配送中心之间的限制？
   • 是否有服务时间要求？

4. 问题目标

   • 最小化总距离？
   • 最小化车辆数量？
   • 平衡各配送中心的工作量？
   • 最小化最长路线的长度？

5. 问题规模

   • 小型（< 50 位客户，2-3 个配送中心）：可使用精确方法
   • 中型（50-200 位客户）：使用高级启发式算法
   • 大型（200+ 位客户）：需要使用元启发式算法

---

数学建模

MDVRP 数学模型

集合：

• D = {1, ..., m}：配送中心集合
• C = {1, ..., n}：客户集合
• V = D ∪ C：所有节点
• K_d：配送中心 d 的车辆集合

参数：

• c_{ij}：节点 i 到 j 的成本／距离
• q_i：客户 i 的需求量
• Q_k：车辆 k 的容量
• M_d：配送中心 d 可用车辆的最大数量

决策变量：

• x_{ijk} ∈ {0,1}：若车辆 k 从 i 行驶到 j，则为 1
• y_{dk} ∈ {0,1}：若车辆 k 从配送中心 d 出发并投入使用，则为 1

目标函数：

Minimize: Σ_{d∈D} Σ_{k∈K_d} Σ_{i∈V} Σ_{j∈V} c_{ij} * x_{ijk}

约束条件：

1. Each customer visited exactly once:
   Σ_{d∈D} Σ_{k∈K_d} Σ_{i∈V} x_{ijk} = 1,  ∀j ∈ C

2. Vehicle starts and ends at same depot:
   Σ_{j∈C} x_{djk} = y_{dk},  ∀d ∈ D, k ∈ K_d
   Σ_{i∈C} x_{idk} = y_{dk},  ∀d ∈ D, k ∈ K_d

3. Flow conservation:
   Σ_{i∈V} x_{ihk} = Σ_{j∈V} x_{hjk},  ∀h ∈ C, ∀d ∈ D, k ∈ K_d

4. Capacity constraint:
   Σ_{i∈C} Σ_{j∈V} q_i * x_{ijk} ≤ Q_k,  ∀d ∈ D, k ∈ K_d

5. Maximum vehicles per depot:
   Σ_{k∈K_d} y_{dk} ≤ M_d,  ∀d ∈ D

6. Subtour elimination constraints

7. Binary variables:
   x_{ijk}, y_{dk} ∈ {0,1}

---

精确算法和启发式算法

1. 先聚类，再规划路线

import numpy as np
from sklearn.cluster import KMeans
import random

def mdvrp_cluster_first(customer_coords, depot_coords, demands,
                       vehicle_capacity, vehicles_per_depot):
    """
    Cluster-first, route-second heuristic for MDVRP

    Phase 1: Assign customers to depots (clustering)
    Phase 2: Solve VRP for each depot

    Args:
        customer_coords: n x 2 array of customer coordinates
        depot_coords: m x 2 array of depot coordinates
        demands: customer demands
        vehicle_capacity: vehicle capacity
        vehicles_per_depot: vehicles available at each depot

    Returns:
        solution dictionary
    """
    n_customers = len(customer_coords)
    n_depots = len(depot_coords)

    # Phase 1: Cluster customers to depots using K-means
    all_coords = np.vstack([depot_coords, customer_coords])

    # Initialize K-means with depot locations
    kmeans = KMeans(n_clusters=n_depots, init=depot_coords, n_init=1)
    customer_clusters = kmeans.fit_predict(customer_coords)

    # Phase 2: Solve VRP for each depot
    depot_routes = []
    total_distance = 0

    for depot_id in range(n_depots):
        # Get customers assigned to this depot
        cluster_customers = [i for i in range(n_customers)
                           if customer_clusters[i] == depot_id]

        if not cluster_customers:
            continue

        # Build distance matrix for this depot and its customers
        depot_coord = depot_coords[depot_id]
        cluster_coords = [depot_coord] + [customer_coords[i]
                                         for i in cluster_customers]

        sub_n = len(cluster_coords)
        sub_dist_matrix = np.zeros((sub_n, sub_n))

        for i in range(sub_n):
            for j in range(sub_n):
                sub_dist_matrix[i][j] = np.linalg.norm(
                    cluster_coords[i] - cluster_coords[j])

        # Demands for this cluster
        cluster_demands = [0] + [demands[i] for i in cluster_customers]

        # Solve VRP for this depot using Clarke-Wright
        routes = clarke_wright_for_depot(
            sub_dist_matrix, cluster_demands,
            vehicle_capacity, vehicles_per_depot)

        # Convert local indices back to global customer indices
        for route in routes:
            global_route = []
            for node in route:
                if node == 0:
                    global_route.append(('depot', depot_id))
                else:
                    global_route.append(('customer', cluster_customers[node-1]))
            depot_routes.append(global_route)

            # Calculate route distance
            route_distance = sum(sub_dist_matrix[route[i]][route[i+1]]
                               for i in range(len(route)-1))
            total_distance += route_distance

    return {
        'routes': depot_routes,
        'customer_to_depot': customer_clusters,
        'total_distance': total_distance,
        'num_vehicles': len(depot_routes)
    }


def clarke_wright_for_depot(dist_matrix, demands, vehicle_capacity, max_vehicles):
    """Clarke-Wright savings for single depot"""
    n = len(dist_matrix)
    customers = list(range(1, n))
    depot = 0

    # Calculate savings
    savings = []
    for i in customers:
        for j in customers:
            if i < j:
                saving = (dist_matrix[depot][i] +
                         dist_matrix[depot][j] -
                         dist_matrix[i][j])
                savings.append((saving, i, j))

    savings.sort(reverse=True)

    # Initialize routes
    routes = [[depot, c, depot] for c in customers]
    route_loads = [demands[c] for c in customers]

    # Merge routes
    for saving_val, i, j in savings:
        if len(routes) <= max_vehicles:
            break

        route_i_idx = next((idx for idx, r in enumerate(routes) if i in r), None)
        route_j_idx = next((idx for idx, r in enumerate(routes) if j in r), None)

        if route_i_idx is None or route_j_idx is None:
            continue
        if route_i_idx == route_j_idx:
            continue

        route_i = routes[route_i_idx]
        route_j = routes[route_j_idx]

        i_at_end = (route_i[1] == i or route_i[-2] == i)
        j_at_end = (route_j[1] == j or route_j[-2] == j)

        if not (i_at_end and j_at_end):
            continue

        combined_load = route_loads[route_i_idx] + route_loads[route_j_idx]
        if combined_load > vehicle_capacity:
            continue

        # Merge
        route_i_interior = route_i[1:-1]
        route_j_interior = route_j[1:-1]

        if route_i_interior[-1] == i and route_j_interior[0] == j:
            new_route = [depot] + route_i_interior + route_j_interior + [depot]
        elif route_i_interior[-1] == i and route_j_interior[-1] == j:
            new_route = [depot] + route_i_interior + route_j_interior[::-1] + [depot]
        elif route_i_interior[0] == i and route_j_interior[0] == j:
            new_route = [depot] + route_i_interior[::-1] + route_j_interior + [depot]
        elif route_i_interior[0] == i and route_j_interior[-1] == j:
            new_route = [depot] + route_j_interior + route_i_interior + [depot]
        else:
            continue

        routes[route_i_idx] = new_route
        route_loads[route_i_idx] = combined_load
        del routes[route_j_idx]
        del route_loads[route_j_idx]

    return routes

2. 多配送中心 OR-Tools 实现

from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp

def solve_mdvrp_ortools(customer_coords, depot_coords, demands,
                       vehicle_capacity, vehicles_per_depot, time_limit=60):
    """
    Solve MDVRP using Google OR-Tools

    Args:
        customer_coords: customer coordinates
        depot_coords: depot coordinates
        demands: customer demands
        vehicle_capacity: vehicle capacity
        vehicles_per_depot: list of vehicles per depot
        time_limit: time limit in seconds

    Returns:
        solution dictionary
    """
    n_customers = len(customer_coords)
    n_depots = len(depot_coords)
    total_vehicles = sum(vehicles_per_depot)

    # Combine all coordinates (depots first, then customers)
    all_coords = np.vstack([depot_coords, customer_coords])
    n_locations = len(all_coords)

    # Build distance matrix
    dist_matrix = np.zeros((n_locations, n_locations))
    for i in range(n_locations):
        for j in range(n_locations):
            dist_matrix[i][j] = np.linalg.norm(all_coords[i] - all_coords[j])

    # Demands (depots have 0 demand)
    all_demands = [0] * n_depots + list(demands)

    # Create starts and ends for each vehicle
    # Each vehicle starts and ends at its assigned depot
    starts = []
    ends = []
    vehicle_to_depot = []

    for depot_id in range(n_depots):
        for _ in range(vehicles_per_depot[depot_id]):
            starts.append(depot_id)
            ends.append(depot_id)
            vehicle_to_depot.append(depot_id)

    # Create routing manager
    manager = pywrapcp.RoutingIndexManager(n_locations, total_vehicles,
                                          starts, ends)

    # Create routing model
    routing = pywrapcp.RoutingModel(manager)

    # Distance callback
    def distance_callback(from_index, to_index):
        from_node = manager.IndexToNode(from_index)
        to_node = manager.IndexToNode(to_index)
        return int(dist_matrix[from_node][to_node] * 100)

    transit_callback_index = routing.RegisterTransitCallback(distance_callback)
    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)

    # Demand callback
    def demand_callback(from_index):
        from_node = manager.IndexToNode(from_index)
        return int(all_demands[from_node])

    demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)

    # Add capacity dimension
    routing.AddDimensionWithVehicleCapacity(
        demand_callback_index,
        0,  # null capacity slack
        [int(vehicle_capacity)] * total_vehicles,
        True,  # start cumul to zero
        'Capacity')

    # Prevent visiting depots as intermediate stops
    for depot_id in range(n_depots):
        depot_index = manager.NodeToIndex(depot_id)
        routing.solver().Add(routing.ActiveVar(depot_index) == 0)

    # Search parameters
    search_parameters = pywrapcp.DefaultRoutingSearchParameters()
    search_parameters.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
    search_parameters.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    search_parameters.time_limit.seconds = time_limit

    # Solve
    solution = routing.SolveWithParameters(search_parameters)

    if solution:
        routes = []
        depot_assignments = {}
        total_distance = 0

        for vehicle_id in range(total_vehicles):
            index = routing.Start(vehicle_id)
            route = []
            route_distance = 0

            while not routing.IsEnd(index):
                node = manager.IndexToNode(index)
                route.append(node)
                previous_index = index
                index = solution.Value(routing.NextVar(index))
                route_distance += routing.GetArcCostForVehicle(
                    previous_index, index, vehicle_id) / 100.0

            route.append(manager.IndexToNode(index))

            if len(route) > 2:  # Route has customers
                depot_id = vehicle_to_depot[vehicle_id]
                routes.append({
                    'vehicle_id': vehicle_id,
                    'depot_id': depot_id,
                    'route': route,
                    'distance': route_distance
                })
                total_distance += route_distance

                # Track customer-depot assignments
                for node in route[1:-1]:  # Skip start/end depot
                    if node >= n_depots:  # Is customer
                        customer_id = node - n_depots
                        depot_assignments[customer_id] = depot_id

        return {
            'status': 'Optimal',
            'routes': routes,
            'depot_assignments': depot_assignments,
            'total_distance': total_distance,
            'num_vehicles_used': len(routes)
        }
    else:
        return {
            'status': 'No solution found',
            'routes': None
        }


# Example usage with visualization
def visualize_mdvrp_solution(customer_coords, depot_coords, routes):
    """Visualize MDVRP solution"""
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots(figsize=(12, 8))

    n_depots = len(depot_coords)
    colors = plt.cm.tab10(np.linspace(0, 1, 10))

    # Plot routes by depot
    for route_info in routes:
        depot_id = route_info['depot_id']
        route = route_info['route']
        vehicle_id = route_info['vehicle_id']

        # Get coordinates for this route
        route_coords = []
        for node in route:
            if node < n_depots:
                route_coords.append(depot_coords[node])
            else:
                customer_id = node - n_depots
                route_coords.append(customer_coords[customer_id])

        xs = [c[0] for c in route_coords]
        ys = [c[1] for c in route_coords]

        ax.plot(xs, ys, 'o-', color=colors[depot_id],
               linewidth=2, markersize=6,
               label=f'Depot {depot_id+1}, Vehicle {vehicle_id+1}',
               alpha=0.7)

    # Plot depots
    for i, depot in enumerate(depot_coords):
        ax.plot(depot[0], depot[1], 's', color=colors[i],
               markersize=20, label=f'Depot {i+1}',
               markeredgecolor='black', markeredgewidth=2, zorder=10)

    # Plot customers
    for customer in customer_coords:
        ax.plot(customer[0], customer[1], 'o', color='lightgray',
               markersize=8, markeredgecolor='black', markeredgewidth=1)

    ax.set_xlabel('X Coordinate')
    ax.set_ylabel('Y Coordinate')
    ax.set_title('Multi-Depot VRP Solution')
    ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()


# Complete example
if __name__ == "__main__":
    np.random.seed(42)
    random.seed(42)

    # Generate problem
    n_depots = 3
    n_customers = 40

    depot_coords = np.random.rand(n_depots, 2) * 100
    customer_coords = np.random.rand(n_customers, 2) * 100

    demands = [random.randint(5, 20) for _ in range(n_customers)]
    vehicle_capacity = 100
    vehicles_per_depot = [3, 3, 2]  # 8 total vehicles

    print(f"Problem: {n_customers} customers, {n_depots} depots")
    print(f"Total demand: {sum(demands)}")
    print(f"Vehicles per depot: {vehicles_per_depot}")

    print("\nSolving MDVRP with OR-Tools...")
    result = solve_mdvrp_ortools(customer_coords, depot_coords, demands,
                                vehicle_capacity, vehicles_per_depot,
                                time_limit=60)

    print(f"\nStatus: {result['status']}")
    print(f"Total Distance: {result['total_distance']:.2f}")
    print(f"Vehicles Used: {result['num_vehicles_used']}")

    # Count customers per depot
    depot_counts = {}
    for customer_id, depot_id in result['depot_assignments'].items():
        depot_counts[depot_id] = depot_counts.get(depot_id, 0) + 1

    print("\nCustomers assigned to each depot:")
    for depot_id in range(n_depots):
        count = depot_counts.get(depot_id, 0)
        print(f"  Depot {depot_id+1}: {count} customers")

    print("\nRoute Details:")
    for route_info in result['routes']:
        depot_id = route_info['depot_id']
        route = route_info['route']
        distance = route_info['distance']
        n_cust = len(route) - 2

        print(f"  Depot {depot_id+1}, Vehicle {route_info['vehicle_id']+1}:")
        print(f"    Customers: {n_cust}, Distance: {distance:.2f}")

    # Visualize
    visualize_mdvrp_solution(customer_coords, depot_coords, result['routes'])

---

工具和库

Python 库

• OR-Tools（Google）：最适合 MDVRP（推荐）
• PuLP/Pyomo：混合整数规划建模
• scikit-learn：用于配送中心分配的 K-means 聚类

方法

• 先聚类，再规划路线：速度快，适合地理位置聚类
• 先规划路线，再聚类：集成效果更好，但更复杂
• 统一方法：同时求解（OR-Tools）

---

常见挑战和解决方案

挑战：配送中心分配不平衡

问题：

• 部分配送中心负荷过重，其他配送中心利用不足
• 地理位置分布不均衡

解决方案：

• 添加工作量平衡约束
• 使用加权聚类
• 进行优化后再平衡

挑战：不允许跨配送中心调配

问题：

• 车辆必须返回出发的配送中心
• 车辆不能为不同于其出发配送中心的客户提供服务

解决方案：

• OR-Tools 可通过设置单独的起点／终点索引自然处理
• 在聚类阶段强制执行

挑战：配送中心容量限制

问题：

• 配送中心存在吞吐量限制
• 每个配送中心的车辆数或客户数有限制

解决方案：

• 添加配送中心容量约束
• 使用带约束的聚类
• 采用多阶段分配

---

输出格式

MDVRP 求解报告

问题：

• 客户：60 位
• 配送中心：3 个
• 车辆总数：10 辆（各配送中心分别为 4、3、3 辆）
• 车辆容量：100 个单位

解：

指标	数值
总距离	1,845 km
使用车辆数	9 / 10
平均路线长度	205 km

配送中心分配：

配送中心	客户数	车辆数	总距离	平均负载
1	24	4	782 km	96%
2	20	3	623 km	94%
3	16	2	440 km	92%

---

需要提出的问题

1. 有多少个车场／仓库？
2. 是否可以从任意车场为客户提供服务？
3. 车辆是否固定分配给特定车场，还是可以灵活调配？
4. 是否应在各车场之间均衡工作量？
5. 车场是否有容量限制？
6. 车辆能否返回与出发时不同的车场？

---

相关技能

• vehicle-routing-problem：用于单车场 VRP
• hub-location-problem：用于车场选址优化
• facility-location-problem：用于车场选址
• capacitated-vrp：用于侧重容量的路径规划
• network-design：用于配送网络设计