从组合爆炸到智能搜索:用 Python 求解容量约束 VRP

def distance_callback(from_index, to_index):    from_node = manager.IndexToNode(from_index)    to_node = manager.IndexToNode(to_index)    return distance_matrix[from_node][to_node]

routing.AddDimensionWithVehicleCapacity(    demand_callback_index,    0,                            # 松弛变量，允许超载时为正值    [vehicle_capacity] * num_vehicles,  # 每辆车的容量上界    True,                         # 从起点开始累加    'Capacity')

import randomimport numpy as npimport matplotlib.pyplot as pltfrom ortools.constraint_solver import pywrapcp, routing_enums_pb2# ==================== 1. 数据生成 ====================random.seed(42)np.random.seed(42)num_customers = 15num_vehicles = 2depot = 0vehicle_capacity = 15# 坐标：仓库在原点，客户在[-10, 10]范围内均匀分布points = [(0, 0)] + [(random.uniform(-10, 10), random.uniform(-10, 10)) for _ in range(num_customers)]# 需求量：1~4 均匀分布，仓库为0demands = [0] + [random.randint(1, 4) for _ in range(num_customers)]# ==================== 2. 距离矩阵 ====================def compute_distance_matrix(points):    """欧几里得距离矩阵，乘以1000保留整数精度"""    size = len(points)    matrix = np.zeros((size, size), dtype=int)    for i in range(size):        for j in range(size):            if i != j:                dx = points[i][0] - points[j][0]                dy = points[i][1] - points[j][1]                matrix[i][j] = int(np.hypot(dx, dy) * 1000)    return matrixdistance_matrix = compute_distance_matrix(points)# ==================== 3. OR-Tools 建模 ====================# RoutingIndexManager：管理节点索引与求解器内部索引的映射manager = pywrapcp.RoutingIndexManager(len(points), num_vehicles, depot)routing = pywrapcp.RoutingModel(manager)# 距离回调 → 目标函数def distance_callback(from_index, to_index):    from_node = manager.IndexToNode(from_index)    to_node = manager.IndexToNode(to_index)    return distance_matrix[from_node][to_node]transit_callback_index = routing.RegisterTransitCallback(distance_callback)routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)# 需求回调 → 容量约束维度def demand_callback(from_index):    from_node = manager.IndexToNode(from_index)    return demands[from_node]demand_callback_index = routing.RegisterUnaryTransitCallback(demand_callback)routing.AddDimensionWithVehicleCapacity(    demand_callback_index,    0,                                # 不允许超载的硬约束    [vehicle_capacity] * num_vehicles,    True,    'Capacity')# 搜索参数：贪心初始解 + 引导式局部搜索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 = 5# ==================== 4. 求解并输出 ====================solution = routing.SolveWithParameters(search_parameters)if solution:    print(f'最优解目标值: {solution.ObjectiveValue()/1000:.2f} (总距离)')    routes = []    for vehicle_id in range(num_vehicles):        route = []        index = routing.Start(vehicle_id)        while not routing.IsEnd(index):            node = manager.IndexToNode(index)            route.append(node)            index = solution.Value(routing.NextVar(index))        route.append(depot)        routes.append(route)        load = sum(demands[node] for node in route if node != depot)        dist = sum(distance_matrix[route[i]][route[i+1]] for i in range(len(route)-1)) / 1000        print(f'车辆 {vehicle_id+1} | 路线: {route} | 载重: {load}/{vehicle_capacity} | 距离: {dist:.2f}')else:    print('未找到可行解。请检查约束或增大容量。')    exit()# ==================== 5. 可视化 ====================plt.style.use('dark_background')fig, ax = plt.subplots(figsize=(12, 10), facecolor='#1a1a2e')ax.set_facecolor('#1a1a2e')palette = ['#00f5ff', '#ff6ec7', '#f5ff00', '#00ff7f', '#ffa500', '#af7ac5']for vid, route in enumerate(routes):    color = palette[vid % len(palette)]    xs = [points[n][0] for n in route]    ys = [points[n][1] for n in route]    # 三次样条平滑    t = np.arange(len(xs))    t_dense = np.linspace(0, len(xs) - 1, 500)    xs_smooth = np.interp(t_dense, t, xs)    ys_smooth = np.interp(t_dense, t, ys)    ax.plot(xs_smooth, ys_smooth, color=color, linewidth=2.5, alpha=0.85, zorder=2)    mid = len(route) // 2    ax.annotate(f'车辆{vid+1}', (xs[mid], ys[mid]),                color=color, fontsize=12, fontweight='bold',                bbox=dict(boxstyle='round,pad=0.3', facecolor='#1a1a2e', edgecolor=color, alpha=0.9))# 客户点 (大小 ∝ 需求量)for i in range(1, len(points)):    d = demands[i]    ax.scatter(points[i][0], points[i][1], s=d*80, c='white',               edgecolors='#00f5ff', linewidths=1.5, zorder=5)    ax.annotate(f'{d}', (points[i][0], points[i][1]),                textcoords="offset points", xytext=(0, 8),                ha='center', color='white', fontsize=9, fontweight='bold')# 仓库ax.scatter(points[0][0], points[0][1], s=350, c='#ff073a', marker='s',           edgecolors='white', linewidths=2, zorder=10, label='仓库 (Depot)')# 装饰ax.legend(loc='upper left', frameon=True, facecolor='#2d2d44',          edgecolor='white', fontsize=10)ax.grid(True, linestyle='--', linewidth=0.3, alpha=0.4, color='white')for spine in ax.spines.values():    spine.set_edgecolor('#444466')ax.tick_params(colors='#cccccc')ax.set_xlabel('X', color='white', fontsize=12)ax.set_ylabel('Y', color='white', fontsize=12)ax.set_title('CVRP 优化结果 · 2 辆车 · 容量 Q=15', color='white',             fontsize=18, fontweight='bold', pad=20)plt.tight_layout()plt.savefig('cvrp_solution.png', dpi=200, bbox_inches='tight', facecolor='#1a1a2e')plt.show()