用Python和Gurobi复现L-Shape方法：一个两阶段随机规划的实战案例

用Python和Gurobi复现L-Shape方法一个两阶段随机规划的实战案例当你在供应链管理中面临需求不确定的库存决策或在能源系统中需要应对可再生能源出力波动时两阶段随机规划提供了将不确定性纳入数学模型的优雅框架。而L-Shape方法正是将这类先决策-后观察-再调整的决策过程转化为可计算算法的关键钥匙。本文将带你从零开始用Python和Gurobi实现一个完整的L-Shape算法解决一个简化的电力系统调度问题——这个案例中的数据和方法可以轻松迁移到生产计划、物流优化等场景。1. 问题建模电力调度中的两阶段决策假设我们管理一个包含火电机组和风电场的微型电网需要制定24小时的发电计划。第一阶段决策日前市场必须提前确定火电机组的启停和基础出力而第二阶段决策实时市场可以根据实际风电出力调整燃气机组的灵活出力。这正是典型的两阶段随机规划问题第一阶段变量火电机组启停状态$u_t$、基础出力$p_t^{base}$第二阶段变量燃气机组调节量$p_t^{adjust}$、弃风量$w_t^{curtail}$随机参数风电实际出力$W_t(\omega)$考虑三种典型场景高/中/低# 场景定义 (单位MW) scenarios { high_wind: [45, 50, 55, ..., 60], # 24小时风电预测 medium_wind: [35, 40, 45, ..., 50], low_wind: [25, 30, 35, ..., 40], prob: [0.3, 0.5, 0.2] # 场景概率 }注意实际应用中应通过历史数据生成更多场景本文为演示简化为3个代表性场景2. L-Shape方法的核心实现步骤2.1 主问题建模主问题包含第一阶段决策变量和初始近似的最优性割平面开始时可能为空def build_master_problem(): m gp.Model(Master_Problem) # 第一阶段变量 u m.addVars(24, vtypeGRB.BINARY, namecommitment) p_base m.addVars(24, lb0, ub100, namebase_generation) # 辅助变量θ近似第二阶段的期望成本 theta m.addVar(lb-GRB.INFINITY, nametheta) # 第一阶段目标 m.setObjective( gp.quicksum(500*u[t] 10*p_base[t] for t in range(24)) theta, GRB.MINIMIZE ) # 初始约束可根据问题添加 m.addConstrs(p_base[t] 30*u[t] for t in range(24)) # 最小出力约束 return m2.2 子问题求解与割平面生成对每个场景求解子问题收集对偶信息生成最优性割def solve_subproblem(master_solution, scenario): sub_m gp.Model(Sub_Problem) # 第二阶段变量 p_adjust sub_m.addVars(24, lb-50, ub50, nameadjustment) w_curtail sub_m.addVars(24, lb0, namecurtailment) # 目标函数调节成本 弃风惩罚 sub_m.setObjective( gp.quicksum(15*abs_(p_adjust[t]) 50*w_curtail[t] for t in range(24)), GRB.MINIMIZE ) # 平衡约束关键连接约束 demand [55, 50, ..., 60] # 24小时负荷 sub_m.addConstrs( master_solution.p_base[t] p_adjust[t] min(scenario[t], master_solution.wind_cap) - w_curtail[t] demand[t] for t in range(24) ) sub_m.optimize() # 提取对偶变量 duals [c.Pi for c in sub_m.getConstrs()] return sub_m.ObjVal, duals生成最优性割平面的数学形式为$$ \theta \geq \sum_{k1}^K p_k \left[ Q(x^k,\xi_k) \pi_k^T T_k (x - x^k) \right] $$对应代码实现def add_optimality_cut(master, x_hat, sub_results): # sub_results包含各场景的目标值和对偶变量 expr gp.LinExpr() for k in range(len(scenarios)): Q_k, pi_k sub_results[k] # 计算割平面系数 coeff sum(pi_k[t] * scenarios[k][t] for t in range(24)) expr scenarios[prob][k] * (Q_k - coeff) # 添加割平面约束 master.addConstr(master.getVarByName(theta) expr)3. 完整算法流程实现将各组件整合为完整的L-Shape算法def L_shape_algorithm(max_iter100, tol1e-4): # 初始化 master build_master_problem() UB float(inf) LB -float(inf) iteration 0 while iteration max_iter and UB - LB tol: # 求解主问题 master.optimize() LB master.ObjVal # 固定第一阶段变量求解各场景子问题 x_hat {p_base: [v.X for v in master.getVars() if base in v.VarName]} sub_results [] for s in [high_wind, medium_wind, low_wind]: obj, duals solve_subproblem(x_hat, scenarios[s]) sub_results.append((obj, duals)) # 计算上界 current_UB (sum(500*x_hat[u][t] 10*x_hat[p_base][t] for t in range(24)) sum(p*sub_results[k][0] for k, p in enumerate(scenarios[prob]))) UB min(UB, current_UB) # 添加最优性割 add_optimality_cut(master, x_hat, sub_results) iteration 1 print(fIter {iteration}: LB{LB:.2f}, UB{UB:.2f}, Gap{(UB-LB)/UB*100:.2f}%) return master4. 结果分析与算法改进运行算法后我们通常会观察到收敛过程呈现以下典型特征迭代次数下界 (LB)上界 (UB)相对间隙112,45038,72067.8%524,38028,15013.4%1026,85027,2301.4%1527,11027,1250.06%关键改进技巧初始割平面添加基于场景平均解的初始割加速早期收敛割平面选择保留活跃割平面避免模型膨胀并行计算各场景子问题可并行求解# 并行求解子问题示例 from concurrent.futures import ThreadPoolExecutor def parallel_solve_subproblems(x_hat): with ThreadPoolExecutor() as executor: futures [] for s in scenarios: future executor.submit(solve_subproblem, x_hat, scenarios[s]) futures.append(future) return [f.result() for f in futures]5. 工程实践中的挑战与解决方案在实际应用中我们往往会遇到以下典型问题内存爆炸当场景数增加到数百个时存储所有子问题模型会消耗大量内存。解决方案是使用延迟生成策略只在需要时构建子问题采用基于Benders分解的分布式计算框架数值不稳定对偶值可能因问题退化出现异常。应对措施包括在Gurobi中设置参数DualReductions0添加正则化项保证子问题强对偶性成立收敛震荡上下界出现振荡时可以引入信任域限制主问题解的变化范围采用割平面聚合技术# 信任域实现示例 def add_trust_region(master, x_prev, delta0.1): for t in range(24): p master.getVarByName(fbase_generation[{t}]) master.addConstr(p x_prev[t] - delta) master.addConstr(p x_prev[t] delta)6. 扩展到更复杂的应用场景本文展示的基础框架可以进一步扩展多阶段问题通过嵌套L-Shape方法处理多阶段决策class MultistageProblem: def __init__(self): self.stages [] # 各阶段模型 self.scenario_tree {} # 情景树结构整数追索当第二阶段包含整数变量时需要使用分支定界框架def solve_integer_subproblem(): sub_m gp.Model(Integer_Sub) sub_m.Params.MIPGap 0.01 # 控制求解精度 # ... 整数变量定义 ... return sub_m风险规避通过条件风险价值(CVaR)等指标替代期望值def add_cvar_constraints(master, alpha0.9): # 添加CVaR相关变量和约束 eta master.addVar(nameVaR) zeta master.addVars(len(scenarios), namezeta) master.addConstrs(zeta[k] scenarios[k].obj - eta for k in scenarios) master.addConstr(theta eta 1/(1-alpha)*sum(p[k]*zeta[k] for k in scenarios))