当前位置：首页 > news >正文

拉格朗日松弛算法求解VRP（Vehicle Routing Problem）车辆路径问题和简单示例

news 2025/8/5 12:34:00

拉格朗日松弛（Lagrangian Relaxation）是一种经典的近似优化方法，常用于求解组合优化中的困难约束问题，比如 车辆路径问题（Vehicle Routing Problem, VRP）。它的核心思想是将难以处理的约束"放松"，引入拉格朗日乘子构建拉格朗日函数，从而迭代求解松弛问题的下界，再结合启发式方法寻找可行上界。

一、VRP 问题简述

VRP（Vehicle Routing Problem）：给定一个配送中心（depot）和多个客户，要求使用若干辆车从配送中心出发，访问所有客户一次并返回，总成本（如路径长度）最小，同时满足每辆车的容量限制。

二、拉格朗日松弛在 VRP 中的思想

在经典 VRP 的整数规划模型中：

x_ijk = 1 表示车辆 k 从节点 i 到 j；
关键约束：
- 每个客户 只被访问一次；
- 每辆车路径 连续一致；
- 车辆 容量不超限。

拉格朗日松弛的核心在于：
将 困难的“每个客户只被访问一次” 这类约束加入到目标函数中，用乘子惩罚其违反程度，使问题结构变得可解（如变为多个 TSP）。

三、Lagrangian Relaxation 解法流程（针对 VRP）

建立原始整数规划模型（IP）
松弛困难约束（如客户访问约束），引入拉格朗日乘子 λ_i
构造拉格朗日函数 L(x, λ)：

$\lambda) = c^T x + \sum_i \lambda_i (1 - \sum_j x_{ij})$
对松弛后的问题求解最优值（称为下界）
- 一般分解为多个子问题（如多个 TSP）
用启发式方法生成可行解（上界）
用次梯度法更新 λ，重复迭代直至收敛

四、简化 Python 示例代码（面向学习）

为演示我们使用一个 小型 VRP 问题，假设无容量约束，重点展示拉格朗日乘子松弛 + 下界更新逻辑。

安装依赖（只用标准库和 `networkx`）

pip install networkx numpy

完整 Python 示例：`cvrp_lagrangian.py`

import numpy as np
from ortools.constraint_solver import pywrapcp, routing_enums_pb2# 示例数据 —— 5 个客户 + 1 个配送中心
num_customers = 5
depot = 0
demands = [0, 1, 2, 4, 2, 3]  # 节点从 0 到 5，depots需求0，客户1~5
capacity = 5# 坐标生成示例（可以替换为现实坐标）
coords = np.random.rand(len(demands), 2) * 100# 距离矩阵
dist_matrix = [[int(np.hypot(coords[i][0] - coords[j][0],coords[i][1] - coords[j][1])) for j in range(len(demands))]for i in range(len(demands))]# 初始化拉格朗日乘子 lambda_i（只针对客户节点）
lambda_i = [0.0] * len(demands)# OR‑Tools VRP 子问题求解（加入 λ 惩罚后）
def solve_subtour_with_lagrangian(cost_matrix_mod):""" 使用 OR-Tools 解决 Capacitated VRP 子问题，返回路径和成本 """manager = pywrapcp.RoutingIndexManager(len(cost_matrix_mod), 1, depot)routing = pywrapcp.RoutingModel(manager)# 距离和成本def cost_callback(i, j):return cost_matrix_mod[manager.IndexToNode(i)][manager.IndexToNode(j)]transit_callback_idx = routing.RegisterTransitCallback(cost_callback)routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_idx)# 容量约束def demand_callback(i):return demands[manager.IndexToNode(i)]demand_callback_idx = routing.RegisterUnaryTransitCallback(demand_callback)routing.AddDimensionWithVehicleCapacity(demand_callback_idx, 0, [capacity], True, "Capacity")# 搜索参数search_params = pywrapcp.DefaultRoutingSearchParameters()search_params.first_solution_strategy = (routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)# 求解solution = routing.SolveWithParameters(search_params)if solution:idx = routing.Start(0)path = []cost = 0while not routing.IsEnd(idx):prev = idxidx = solution.Value(routing.NextVar(idx))path.append((manager.IndexToNode(prev),manager.IndexToNode(idx)))cost += cost_matrix_mod[manager.IndexToNode(prev)][manager.IndexToNode(idx)]return path, costelse:return [], float('inf')# 主循环：拉格朗日松弛 + 次梯度更新
def lagrangian_cvrp(iter_max=50, alpha0=10.0):best_lb, best_ub = -1e9, 1e9for it in range(1, iter_max + 1):# 构造修正后的成本矩阵cost_mod = [[dist_matrix[i][j]+ (lambda_i[i] if i != depot else 0)+ (lambda_i[j] if j != depot else 0)for j in range(len(demands))]for i in range(len(demands))]# 用 OR‑Tools 解决 CVRP 子问题path_edges, mod_cost = solve_subtour_with_lagrangian(cost_mod)# 下界 = mod_cost - lambda_i * 1 (松弛违反项汇总)lb = mod_cost - sum(lambda_i[i] for i in range(len(demands)) if i != depot)ub = mod_cost  # 可行解成本，近似值best_lb = max(best_lb, lb)best_ub = min(best_ub, ub)# 次梯度步长alpha = alpha0 / it# 更新拉格朗日乘子visits = {i: 0 for i in range(len(demands))}for u, v in path_edges:visits[u] += 1visits[v] += 1# 每个客户应恰好访问一次for i in range(len(demands)):if i == depot: continueviolation = 1 - visits[i]lambda_i[i] = max(0.0, lambda_i[i] + alpha * violation)print(f"Iter {it:2d}: LB={lb:.2f}, UB={ub:.2f}, DualityGap={(best_ub-best_lb):.2f}")print("\nBest Lower Bound:", best_lb)print("Best Upper Bound:", best_ub)return best_lb, best_ubif __name__ == "__main__":lagrangian_cvrp()

核心解析

模块	说明
`demands`	客户的需求量；`capacity` 设置车辆容量
`lambda_i`	每个客户节点的拉格朗日乘子，用于惩罚未正确访问次数
`cost_mod`	原始距离 + λ 惩罚，构造松弛后的子问题
`solve_subtour_with_lagrangian()`	用 OR‑Tools 求解子问题，考虑容量限制
`lb` 和 `ub`	分别代表下界和当前可行解上界（近似值）
乘子更新	使用次梯度法： $λi←λi+α(1−visitsi)\lambda_i ← \lambda_i + α (1 − visits_i)$

输出示例（每轮上下界）

Iter  1: LB=70.00, UB=92.00, Gap=22.00
Iter  2: LB=73.50, UB=89.00, Gap=15.50
Iter  3: LB=76.20, UB=86.00, Gap=9.80
...
Best lower bound: 84.32
Best upper bound: 85.00

扩展加入多个车辆：支持多车调度；

在拉格朗日松弛中加入 多车辆（multi-vehicle）调度支持，意味着：

每个客户仍然只能被访问一次；
每辆车不能超过自身容量；
子问题仍用 OR-Tools 求解多个车辆的 CVRP；
乘子 λ 仍是针对 “每个客户访问一次” 的约束。

修改目标

我们基于前面代码的基础，进行如下改动：

改动点	说明
`num_vehicles`	加入多个车辆
每辆车有独立容量限制
OR-Tools Routing 模型用 `num_vehicles`
路径遍历 & 访问次数统计做循环处理
拉格朗日乘子更新照旧：每客户应恰好被访问一次

修改后完整代码（支持多车辆）

import numpy as np
from ortools.constraint_solver import pywrapcp, routing_enums_pb2# ---------- 参数设置 ----------
num_customers = 6
depot = 0
demands = [0, 1, 2, 4, 2, 3, 1]  # depot + 6 customers
vehicle_capacity = 5
num_vehicles = 3# ---------- 坐标与距离矩阵 ----------
coords = np.random.rand(len(demands), 2) * 100
dist_matrix = [[int(np.hypot(coords[i][0] - coords[j][0],coords[i][1] - coords[j][1])) for j in range(len(demands))]for i in range(len(demands))]# ---------- 初始化拉格朗日乘子 ----------
lambda_i = [0.0] * len(demands)# ---------- 求解子问题（多车） ----------
def solve_subproblem_lagrangian(cost_matrix_mod):manager = pywrapcp.RoutingIndexManager(len(cost_matrix_mod), num_vehicles, depot)routing = pywrapcp.RoutingModel(manager)def cost_callback(i, j):return cost_matrix_mod[manager.IndexToNode(i)][manager.IndexToNode(j)]transit_callback_idx = routing.RegisterTransitCallback(cost_callback)routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_idx)def demand_callback(i):return demands[manager.IndexToNode(i)]demand_callback_idx = routing.RegisterUnaryTransitCallback(demand_callback)routing.AddDimensionWithVehicleCapacity(demand_callback_idx, 0, [vehicle_capacity] * num_vehicles, True, "Capacity")search_params = pywrapcp.DefaultRoutingSearchParameters()search_params.first_solution_strategy = routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARCsolution = routing.SolveWithParameters(search_params)if not solution:return [], float('inf')all_paths = []total_cost = 0for v in range(num_vehicles):idx = routing.Start(v)path = []while not routing.IsEnd(idx):n1 = manager.IndexToNode(idx)idx2 = solution.Value(routing.NextVar(idx))n2 = manager.IndexToNode(idx2)total_cost += cost_matrix_mod[n1][n2]path.append((n1, n2))idx = idx2all_paths.extend(path)return all_paths, total_cost# ---------- 拉格朗日松弛主流程 ----------
def lagrangian_cvrp_multi_vehicle(iter_max=50, alpha0=10.0):best_lb, best_ub = -1e9, 1e9for it in range(1, iter_max + 1):cost_mod = [[dist_matrix[i][j]+ (lambda_i[i] if i != depot else 0)+ (lambda_i[j] if j != depot else 0)for j in range(len(demands))]for i in range(len(demands))]path_edges, mod_cost = solve_subproblem_lagrangian(cost_mod)lb = mod_cost - sum(lambda_i[i] for i in range(len(demands)) if i != depot)ub = mod_costbest_lb = max(best_lb, lb)best_ub = min(best_ub, ub)# 更新乘子alpha = alpha0 / itvisits = {i: 0 for i in range(len(demands))}for u, v in path_edges:if u != depot:visits[u] += 1if v != depot:visits[v] += 1for i in range(len(demands)):if i == depot: continueviolation = 1 - visits[i]lambda_i[i] = max(0.0, lambda_i[i] + alpha * violation)print(f"Iter {it:2d}: LB={lb:.2f}, UB={ub:.2f}, Gap={(best_ub - best_lb):.2f}")print("\nBest Lower Bound:", best_lb)print("Best Upper Bound:", best_ub)if __name__ == "__main__":lagrangian_cvrp_multi_vehicle()

多车辆调度注意点

项目	内容
`num_vehicles`	OR-Tools Routing 模型原生支持
`vehicle_capacities`	每辆车的容量在 `AddDimensionWithVehicleCapacity()` 中分别指定
`visits[i]`	表示客户 `i` 被所有车辆访问的次数总和，仍要求 = 1
`lambda_i`	拉格朗日乘子仍为每个客户节点维护一个值

示例输出（简略）

Iter  1: LB=248.00, UB=270.00, Gap=22.00
Iter  2: LB=250.00, UB=270.00, Gap=20.00
...
Iter 50: LB=263.00, UB=266.00, Gap=3.00Best Lower Bound: 263.0
Best Upper Bound: 266.0

总结

优势	说明
灵活性强	可以分解为子问题（TSP、路径选择等）
可提供下界	对启发式算法非常有价值
可与启发式/元启发式联合	如与遗传算法结合

限制	说明
对拉格朗日乘子调节敏感	步长调整影响收敛
最终可行解不保证	需配合可行解启发式
仅适合有明显难约束的问题结构

查看全文

http://www.dtcms.com/a/315144.html