# Python实战：拉格朗日松弛算法优化共享单车调度问题（附完整代码） 共享单车调度问题是城市交通管理中的经典优化难题。如何在有限的资源下，将单车从过剩区域调运到需求区域，同时最小化总运输成本？本文将带你用Python实现拉格朗日松弛算法，解决这一实际问题。 ## 1. 问题建模与拉格朗日松弛原理 共享单车调度问题可以建模为一个选址问题（p-中值问题）。假设： - **需求点集合** I：需要单车的区域 - **候选投放点集合** J：可以投放单车的站点 - **决策变量**： - x_j：是否选择站点j（0-1变量） - y_ij：需求点i是否由站点j服务（0-1变量） 数学模型如下： ```python # 目标函数：最小化总运输成本 minimize ∑(i∈I)∑(j∈J) d_i * c_ij * y_ij # 约束条件 subject to: ∑(j∈J) y_ij = 1, ∀i∈I # 每个需求点必须被一个站点服务 y_ij ≤ x_j, ∀i∈I, ∀j∈J # 只有被选中的站点才能服务需求点 ∑(j∈J) x_j = p # 总共选择p个站点 x_j ∈ {0,1}, y_ij ∈ {0,1} # 决策变量为二进制 ``` 拉格朗日松弛的核心思想是将复杂约束（y_ij ≤ x_j）松弛到目标函数中，通过乘子λ_ij惩罚违反约束的情况： ```python # 松弛后的拉格朗日函数 L(x,y,λ) = ∑(i∈I)∑(j∈J) d_i*c_ij*y_ij + ∑(i∈I)∑(j∈J) λ_ij*(y_ij - x_j) ``` ## 2. 数据准备与预处理 我们使用摩拜单车上海数据集（2016年8月抽样数据），包含约10万条订单记录。首先进行数据清洗和聚类： ```python import pandas as pd import numpy as np from sklearn.cluster import KMeans # 读取数据 data = pd.read_csv("mobike_shanghai_sample.csv") # 使用K-means聚类生成需求点 cluster_num = 10 # 需求点数量 kmeans = KMeans(n_clusters=cluster_num) X = data[["start_location_x", "start_location_y"]] kmeans.fit(X) centers = kmeans.cluster_centers_[:, :2] # 聚类中心作为需求点 # 计算每个需求点的需求量 labels = kmeans.labels_ d_set = np.zeros(cluster_num) for item in labels: d_set[item] += 1 # 选择需求量最大的5个点作为候选站点 wait_list_num = 5 top_index_list = np.argsort(-d_set)[:wait_list_num] ``` ## 3. 算法实现：次梯度优化与子问题求解 拉格朗日松弛算法包含三个关键部分：子问题求解、次梯度计算和乘子更新。 ### 3.1 子问题求解 **子问题1**（关于y的优化）： ```python def solve_subproblem_1(c_set, d_set, lambda_set): """ 求解y子问题：最小化 ∑(d_i*c_ij + λ_ij)*y_ij 约束：每行y_ij有且仅有一个1 """ A = d_set.reshape(-1,1) * c_set + lambda_set sol_matrix = np.zeros_like(A) sol_value = 0 for i in range(len(A)): j_min = np.argmin(A[i]) # 每行选择最小系数 sol_matrix[i, j_min] = 1 sol_value += A[i, j_min] return sol_matrix, sol_value ``` **子问题2**（关于x的优化）： ```python def solve_subproblem_2(p, lambda_set): """ 求解x子问题：最大化 ∑(∑λ_ij)*x_j 约束：选择p个站点 """ sum_lambda = lambda_set.sum(axis=0) top_p_indices = np.argpartition(-sum_lambda, p)[:p] sol_vector = np.zeros(len(sum_lambda)) sol_vector[top_p_indices] = 1 sol_value = sum_lambda[top_p_indices].sum() return sol_vector, sol_value ``` ### 3.2 次梯度计算与乘子更新 ```python def get_subgradient(y_matrix, x_vector): """计算次梯度：y_ij - x_j""" return y_matrix - x_vector def update_multipliers(lambda_set, subgrad, theta): """更新拉格朗日乘子""" new_lambda = lambda_set + theta * subgrad return np.maximum(new_lambda, 0) # 保持乘子非负 ``` ## 4. 主算法流程与可行化处理 完整的主算法实现如下： ```python def lagrangian_heuristic(c_set, d_set, p, max_iter=100): # 初始化 lambda_set = 4 * np.ones((len(d_set), len(c_set[0]))) theta = 2 rho = 0.9 # 步长衰减系数 best_feasible_val = float('inf') best_x, best_y = None, None for iter in range(max_iter): # 1. 求解子问题 y_matrix, z_lr1 = solve_subproblem_1(c_set, d_set, lambda_set) x_vector, z_lr2 = solve_subproblem_2(p, lambda_set) z_lr = z_lr1 + z_lr2 # 2. 计算次梯度 subgrad = get_subgradient(y_matrix, x_vector) # 3. 更新乘子和步长 theta *= rho**iter lambda_set = update_multipliers(lambda_set, subgrad, theta) # 4. 可行化处理 infeasible = np.any(y_matrix > x_vector) if not infeasible: # 解可行 current_val = np.sum(d_set.reshape(-1,1) * c_set * y_matrix) if current_val < best_feasible_val: best_feasible_val = current_val best_x, best_y = x_vector.copy(), y_matrix.copy() else: # 解不可行，使用启发式修正 y_feasible = y_matrix.copy() for i in range(len(y_feasible)): for j in range(len(y_feasible[0])): if y_feasible[i,j] > x_vector[j]: # 找到该需求点可用的已选站点中距离最近的 feasible_sites = np.where(x_vector == 1)[0] closest = feasible_sites[np.argmin(c_set[i, feasible_sites])] y_feasible[i,:] = 0 y_feasible[i, closest] = 1 current_val = np.sum(d_set.reshape(-1,1) * c_set * y_feasible) if current_val < best_feasible_val: best_feasible_val = current_val best_x, best_y = x_vector.copy(), y_feasible.copy() return best_x, best_y, best_feasible_val ``` ## 5. 结果分析与可视化 运行算法并分析结果： ```python # 计算距离矩阵 def manhattan_distance(a, b): return np.sum(np.abs(a - b)) c_set = np.zeros((cluster_num, wait_list_num)) for i in range(cluster_num): for j in range(wait_list_num): c_set[i,j] = manhattan_distance(centers[i], centers[top_index_list[j]]) # 运行算法 p = 3 # 选择3个站点 best_x, best_y, best_val = lagrangian_heuristic(c_set, d_set, p) # 可视化结果 import matplotlib.pyplot as plt plt.figure(figsize=(10,6)) plt.scatter(centers[:,0], centers[:,1], c='blue', s=d_set*10, alpha=0.5, label='需求点(大小表示需求量)') plt.scatter(centers[top_index_list,0], centers[top_index_list,1], c='red', s=100, label='候选站点') selected_sites = np.where(best_x == 1)[0] plt.scatter(centers[top_index_list[selected_sites],0], centers[top_index_list[selected_sites],1], c='green', s=200, marker='*', label='选中站点') # 绘制服务关系 for i in range(cluster_num): j_selected = np.argmax(best_y[i]) plt.plot([centers[i,0], centers[top_index_list[j_selected],0]], [centers[i,1], centers[top_index_list[j_selected],1]], 'gray', alpha=0.3) plt.legend() plt.title(f'共享单车调度优化结果(总成本:{best_val:.2f})') plt.xlabel('经度') plt.ylabel('纬度') plt.grid(True) plt.show() ``` ## 6. 性能优化与扩展 实际应用中，我们可以通过以下方式提升算法性能： 1. **自适应步长调整**：根据次梯度范数动态调整步长 2. **并行计算**：子问题可以并行求解 3. **混合整数规划求解器**：对子问题使用专业求解器如Gurobi 4. **启发式改进**：在可行化阶段使用更复杂的启发式规则 ```python # 使用Gurobi求解子问题的示例 import gurobipy as gp from gurobipy import GRB def solve_subproblem_2_with_gurobi(p, lambda_set): m = gp.Model("subproblem2") x = m.addVars(len(lambda_set[0]), vtype=GRB.BINARY, name="x") # 目标函数 obj = gp.quicksum(lambda_set[:,j].sum() * x[j] for j in range(len(lambda_set[0]))) m.setObjective(obj, GRB.MAXIMIZE) # 约束 m.addConstr(gp.quicksum(x[j] for j in range(len(lambda_set[0]))) == p, "p_constraint") m.optimize() sol_vector = np.array([x[j].X for j in range(len(lambda_set[0]))]) sol_value = m.ObjVal return sol_vector, sol_value ``` 拉格朗日松弛算法在共享单车调度问题上展现了良好的性能。在我的实际测试中，对于10个需求点和5个候选站点的问题，算法能在秒级内找到接近最优的解。当问题规模扩大到50个需求点和20个候选站点时，与精确求解器相比，拉格朗日松弛在保持解质量的同时，计算时间优势更加明显。