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

用Python+Gurobi复现Benders分解算法：一个供应链优化问题的完整建模与求解过程

news 2026/5/26 1:20:10

用Python+Gurobi实战Benders分解：供应链网络优化建模与求解全流程

在供应链管理中，如何以最低成本满足客户需求始终是核心挑战。想象一个跨国零售商需要决定在亚洲地区建立多少个区域仓库，每个仓库服务哪些城市，以及如何规划运输路线。这类问题往往涉及数百万美元的运营成本差异，而Benders分解算法正是解决此类大规模混合整数规划问题的利器。本文将带您从零实现一个完整的仓库-配送中心选址优化方案，通过Python和Gurobi的商业求解器组合，揭开Benders分解在真实业务场景中的应用面纱。

1. 问题建模：多级供应链网络优化

我们考虑一个典型的两级供应链网络：

决策变量：
- 二进制变量$y_j$表示是否在位置$j$建立仓库（1=建设，0=不建）
- 连续变量$x_{ij}$表示从仓库$j$到客户$i$的配送量
目标函数：
```
\min \sum_j f_j y_j + \sum_{i,j} c_{ij} x_{ij}
```
其中$f_j$是仓库$j$的固定建设成本，$c_{ij}$是单位运输成本
约束条件：
- 每个客户需求必须被满足：$\sum_j x_{ij} \geq d_i \quad \forall i$
- 仓库流量平衡：$\sum_i x_{ij} \leq M y_j \quad \forall j$（$M$为足够大的常数）

关键建模技巧：

# Gurobi建模片段 model = gp.Model('supply_chain') y = model.addVars(warehouses, vtype=GRB.BINARY, name='open') x = model.addVars(customers, warehouses, name='ship') model.setObjective( gp.quicksum(f[j]*y[j] for j in warehouses) + gp.quicksum(c[i,j]*x[i,j] for i in customers for j in warehouses), GRB.MINIMIZE )

2. Benders分解实现框架

将原问题分解为：

主问题：处理整数变量（仓库选址决策）
子问题：处理连续变量（最优运输方案）

2.1 主问题构建

初始主问题仅包含选址决策：

master = gp.Model('master') y = master.addVars(warehouses, vtype=GRB.BINARY, name='open') q = master.addVar(lb=-GRB.INFINITY, name='auxiliary') # 初始目标函数 master.setObjective( gp.quicksum(f[j]*y[j] for j in warehouses) + q, GRB.MINIMIZE )

2.2 子问题对偶形式

固定$y$值后，子问题的对偶变量$\alpha_i$（客户需求对偶）和$\beta_j$（仓库容量对偶）构成关键切割：

def solve_subproblem(y_fixed): sub = gp.Model('subproblem') alpha = sub.addVars(customers, name='alpha') beta = sub.addVars(warehouses, ub=0, name='beta') # 注意上界约束 sub.setObjective( gp.quicksum(d[i]*alpha[i] for i in customers) + gp.quicksum(M*y_fixed[j]*beta[j] for j in warehouses), GRB.MAXIMIZE ) # 对偶约束 sub.addConstrs( (alpha[i] + beta[j] <= c[i,j] for i in customers for j in warehouses), name='dual_constraint' ) sub.optimize() return sub, alpha, beta

3. 切割平面生成与迭代

Benders算法的核心在于动态生成两种切割：

3.1 可行性切割

当子问题无界时，添加极射线切割：

# 获取极射线后添加约束 feas_cut = master.addConstr( gp.quicksum(d[i]*alpha_ray[i] for i in customers) + gp.quicksum(M*y[j]*beta_ray[j] for j in warehouses) <= 0, name=f'feas_cut_{iter}' )

3.2 最优性切割

当子问题有解时，添加最优性切割：

opt_cut = master.addConstr( q >= gp.quicksum(d[i]*alpha_val[i] for i in customers) + gp.quicksum(M*y[j]*beta_val[j] for j in warehouses), name=f'opt_cut_{iter}' )

迭代控制逻辑：

while UB - LB > tolerance: # 求解主问题 master.optimize() LB = master.objVal # 求解子问题 sub, alpha, beta = solve_subproblem([y[j].X for j in warehouses]) if sub.status == GRB.UNBOUNDED: # 处理无界情况 add_feasibility_cut(alpha.UnbdRay, beta.UnbdRay) else: # 更新上界并添加最优性切割 current_UB = ... # 计算当前上界 if current_UB < UB: UB = current_UB add_optimality_cut(alpha.X, beta.X)

4. 实战技巧与性能优化

4.1 初始解加速

提供合理的初始解可显著减少迭代次数：

# 启发式初始解：选择运输成本最低的仓库组合 def heuristic_initial(): transport_cost = { j: sum(min(c[i,j] for j in warehouses) for i in customers) for j in warehouses } return sorted(warehouses, key=lambda j: f[j] + transport_cost[j])[:3]

4.2 切割管理策略

策略类型	实现方法	适用场景
惰性约束	使用Gurobi的LazyConstraint	大规模问题
池化切割	维护切割池定期筛选	内存有限时
聚合切割	合并相似切割	减少约束数量

4.3 收敛监控可视化

import matplotlib.pyplot as plt def plot_convergence(history): plt.plot(history['LB'], label='Lower Bound') plt.plot(history['UB'], label='Upper Bound') plt.xlabel('Iteration') plt.ylabel('Objective Value') plt.legend() plt.show()