简介微电网经济调度面临光伏、风电等可再生能源出力的强不确定性传统确定性优化在预测偏差下易导致弃光、切负荷等问题。鲁棒优化通过构建不确定集刻画最坏场景以两阶段决策框架实现“先决策、后补救”其中列与约束生成CCG算法将含max-min的复杂问题分解为主问题与子问题迭代求解有效平衡经济性与鲁棒性。基于Python与Gurobi/Pyomo的工程化实现使得盒式不确定集下的鲁棒调度模型可快速部署于实际微电网系统帮助运行人员在面对极端场景时仍能维持供电可靠性与经济性为园区微电网或海岛微电网提供稳健的调度策略。 做微电网经济调度的人基本都遇到过这种情况光伏预测曲线写的是中午能出 300 kW结果当天一片云飘过来直接掉到 80 kW下午负荷却比预计高了 10%。如果按确定性模型排好的机组出力、储能充放电计划去执行轻则弃光/切负荷重则频率波动甚至触发保护动作。我早期用确定性优化做日内调度时靠的是“预测误差留备用”这种粗放手段结果备用留少了要挨罚留多了又白花钱——这本质上就是在不确定性里赌运气。所以后来我把目光投向了两阶段鲁棒优化。它不赌预测值,而是把所有可能的不确定场景光伏、风电、负荷的波动区间全部摆到台面上目标函数在最坏情况下依然经济、可行。这篇文章就基于我近期整理的一个 Python 完整工程把微电网两阶段鲁棒经济调度的模型原理、CCG列与约束生成求解机制、Python 代码结构以及调试中踩过的坑全部拆开讲一遍直接把源码放到工程包里你拿到手改个数据就能跑。1. 从确定性调度到两阶段鲁棒优化为什么必须考虑“最坏情况”1.1 确定性经济调度到底缺了什么传统微电网经济调度的目标很简单在满足负荷、机组出力上下限、储能 SOC 约束的前提下把总运行成本降到最低。通常用混合整数线性规划MILP就能建出来决策变量机组启停状态、各时段有功出力、储能充放电功率、与配网购售电功率目标函数燃料成本 启停成本 购电成本 − 售电收益约束功率平衡、机组爬坡、储能能量状态、联络线容量。写成我常用的数学形式确定性模型长这样[ \min_{x,y} ; c^T x d^T y ][ \text{s.t. } Ax By \ge b,\quad y \ge 0 ]其中 (x) 表示 0-1 变量机组启停(y) 表示连续变量出力、储能量等。这个模型的求解难度不高可问题在于负荷、光伏、风电全部当作已知参数模型算出来的方案一旦遇到与预测不符的实际场景约束立刻可能被破坏。1.2 不确定性建模的核心思想先决策后补救两阶段鲁棒优化的思路跟确定性模型完全不同。它把策略拆成两个阶段第一阶段here-and-now在不确定性实现之前就要定下的决策。例如机组开几台、储能是否处于可调度状态、是否签约购电容量。这些决策一旦定下很难在短时间内改变。第二阶段wait-and-see不确定性实际发生之后调度员有权利做出一系列“补救动作”。例如多开一台机、调整储能放电功率、临时多购入电力等。目标从“最适合预测值”变成了“在最坏不确定性场景下总成本依然最低”。用数学语言就是[ \min_{x \in \mathbb{X}} \left( c^T x \max_{u \in \mathcal{U}} \min_{y \in \mathcal{Y}(x,u)} d^T y \right) ]这个 max-min 嵌套结构是整个问题的核心属于典型的 NP-hard 问题但工程上可以用 CCG 这种迭代算法高效逼近最优解。1.3 为什么盒式不确定集在工程里最常用不确定集 (\mathcal{U}) 有很多种定义方式常见的有盒式box、椭球、预算、数据驱动等。微电网工程里绝大多数代码用的是盒式集合因为简单、直观、不依赖太多历史数据的概率分布假设而且跟鲁棒优化的强对偶转化配合得非常干净[ \mathcal{U} { u \in \mathbb{R}^{n_d} : u_i^{min} \le u_i \le u_i^{max}, ; i1,\dots,n_d } ]这里的”不确定参数“在微电网里主要集中在三个地方光伏出力、风电出力和负荷功率还可以把与上级电网的交换功率也纳入。盒式集合的好处是只要给定每个时段预测值上下浮动比例就能快速构造可行域运算量和集合复杂度都可控。2. 微电网鲁棒调度模型的数学表达与约束拆解2.1 目标函数一阶段成本 二阶段再调度成本用下标 (t) 表示时段(i) 表示机组整个模型的目标函数可以写成[ \min_{\substack{u_{it}, \alpha_t \ P_{it}^g, P_{t}^{b}, P_{t}^{b-}}} ; \underbrace{\sum_{t,i} \left( a_i P_{it}^g b_i u_{it} C_i^{su} v_{it} \right) \sum_t \gamma_t^ P_t^{b}}{\text{第一阶段成本}} ; ; \underbrace{\max{u_d \in \mathcal{U}} \min_{ \Delta P } ; \sum_{t,i} r_i \Delta P_{it}^{g} s_i \Delta P_{it}^{g-}}_{\text{第二阶段补救成本}} ]这里 (u_{it}) 是机组启停 0-1 变量(P_{it}^g) 是机组出力(P_t^{b}/P_t^{b-}) 分别是购电/售电功率(\Delta P) 是第二阶段对出力的调整量(r_i, s_i) 分别是上调/下调的惩罚成本。读者看到的很多开源代码里第二阶段往往只写“再调度成本最小化”各项系数取正数原因是用它衡量最坏情况下的运行费用。如果第二阶段成本系数写得太小求解器会倾向于在一阶段做激进的低成本决策然后在第二阶段用大量调节来兜底得到的结果会失去鲁棒性。2.2 第一阶段约束与决策逻辑第一阶段约束主要包括功率平衡约束[ \sum_i P_{it}^g P_t^{wt} P_t^{pv} P_t^{b} P_t^{dis} D_t P_t^{b-} P_t^{ch} ]其中 (D_t) 是负荷预测值(P_t^{ch}, P_t^{dis}) 是储能充放电功率。这里先把不确定量写成预测值实际波动放到第二阶段约束里处理。机组出力上下限与爬坡约束[ u_{it} P_i^{min} \le P_{it}^g \le u_{it} P_i^{max} ][ -P_i^{ramp} \le P_{it}^g - P_{i,t-1}^g \le P_i^{ramp} ]储能动态约束储能的核心约束是 SOC 递推[ SOC_{t1} SOC_t \eta^{ch} P_t^{ch} - \frac{P_t^{dis}}{\eta^{dis}} ]再加上充放电功率上下限、SOC 上下限以及避免同时充放电的互补约束。工程实现里为了避免二进制变量爆炸有不少代码直接用 (SOC) 和可放电功率作为状态变量建模不过这会牺牲一点精度我在工程包里保留的是标准二进制约束版本。2.3 第二阶段约束把不确定性“放进来”第二阶段的本质是给定第一阶段决策 ((u_{it}, P_{it}^g))在每一个不确定性时段的实际出力/负荷取值下找一套可行的补救措施。为简化推导很多源码里把不确定性集中成不确定功率 (w_t)即净负荷波动[ D_t^{real} - (P_t^{pv,real} P_t^{wt,real}) D_t - (P_t^{pv} P_t^{wt}) w_t ]其中 (w_t \in [w_t^{min}, w_t^{max}])。同时第二阶段约束要求在加入补救变量 (\Delta P) 后功率平衡再次成立[ \sum_i \Delta P_{it} \Delta P^{b}_t - \Delta P^{b-}_t w_t ]并且所有补救变量都要受到爬坡约束、线路容量约束、储能剩余能量约束的限制。3. 列与约束生成CCG算法迭代逻辑3.1 主问题MP的作用收敛最坏情况下的成本下界CCG 的基本思路是把两层问题变成“主问题—子问题”的迭代。主问题是一个包含有限个场景的 MILP[ \min_{x, y_k} ; c^T x \theta ][ \text{s.t. } Ax By_k \ge b_k, \quad \forall k \in \mathcal{K} ][ \theta \ge d^T y_k, \quad \forall k \in \mathcal{K} ]每轮迭代都会把上一步计算出来的“最坏场景 (u^*)”以具体的常数代入主问题生成一组新的变量 (y_k) 和对应约束这叫“列与约束生成”。主问题求解后得到的目标值就是原问题的最优下界LB。说到底主问题干的活就是把不确定性从集合里抽出来变成一个个具体场景然后在这些场景同时满足的前提下做经济调度。每多一轮迭代主问题约束多一组LB 只升不降。3.2 子问题SP的作用找到最坏场景子问题要回答的是在第一阶段决策给定时找哪个不确定性场景 (u) 会让总成本最高从而迫使主问题做更保守的决策。形式化写出来是[ \max_{u \in \mathcal{U}} \min_{y} ; d^T y ][ \text{s.t. } Wy \ge h - Tx - Pu,\quad y \ge 0 ]内层是个线性规划外层是 max。处理办法是先对内层取对偶把 max-min 问题转化成 max 问题[ \max_{u, \lambda} ; \lambda^T (h - Tx - Pu) ][ \text{s.t. } W^T \lambda \le d,\quad \lambda \ge 0,\quad u \in \mathcal{U} ]这时问题里出现了一个双线性项 (-\lambda^T P u)。常规做法是用 KKT/线性化技巧处理或者直接利用盒式集合顶点在边界取最优的性质把问题进一步化解。3.3 顶点枚举 vs 对偶线性化两种主流实现方式我在工程包里同时实现了两种子问题求解方式路径一顶点枚举法。由于盒式集合的可行域是多面体线性目标的最优解必然在顶点取得。所以可以直接枚举每个不确定参数的上/下边界组合对每个组合求解一次线性规划取最坏值。这个方法的优点是写起来直白不易错缺点是场景数量随不确定参数个数指数增长。如果只对光伏、风电、负荷做单时段 24 个波动区间那组合高达 (2^{72}) 根本不可能枚举所以一般只在减少不确定参数个数如只聚合总净负荷时用。路径二对偶 大 M 线性化。引入辅助变量 (z \lambda^T P)通过大 M 将双线性项线性化。虽然会增加不少约束但整个子问题变成一个可以直接交给 Gurobi/CBC 求解的单层 MILP。这种方式在工程源码里更通用也是我推荐你在实际项目里用的版本。4. Python 与 Gurobi/Pyomo 实战从模型搭建到 CCG 主循环4.1 数据准备与全局参数我建议工程文件的目录结构按照“数据—模型—算法—输出”四条线来组织最省心的格式是 CSV 或字典。下面是一份微电网测试系统的简化数据示例参数数值单位机组数量2台时段数24h光伏预测峰值400kW风机预测峰值200kW负荷预测峰值600kW不确定性波动比例光伏/风电/负荷±20% / ±20% / ±10%-储能容量800kWh储能功率上限200kW购电价格分时峰/平/谷1.2/0.8/0.5元/kWh4.2 用 Pyomo 搭建第一阶段 MILPPyomo 是最适合做教学和科研复现的建模工具。你只要把目标函数和约束按数学表达式几乎 1:1 翻译过来即可import pyomo.environ as pyo def build_master_problem(uncertainty_scenarios, param): model pyo.ConcreteModel() # 时段与机组索引 model.T pyo.Set(initializerange(param[T])) model.I pyo.Set(initializerange(param[n_unit])) model.x pyo.Var(model.I, model.T, withinpyo.Binary) # 启停 model.pg pyo.Var(model.I, model.T, withinpyo.NonNegativeReals) model.pbuy pyo.Var(model.T, withinpyo.NonNegativeReals) model.psell pyo.Var(model.T, withinpyo.NonNegativeReals) model.theta pyo.Var(withinpyo.NonNegativeReals) # 第二阶段成本代理 # 目标函数一阶段成本 theta def obj_rule(m): return (sum(m.pg[i, t] * param[fuel_cost][i] m.x[i, t] * param[no_load_cost][i] for i in m.I for t in m.T) sum(param[buy_price][t] * m.pbuy[t] - param[sell_price][t] * m.psell[t] for t in m.T) m.theta) model.obj pyo.Objective(ruleobj_rule, sensepyo.minimize) # 功率平衡在各添加场景下逐一考虑…… return model要点在于主问题里的变量 (y_k) 是“按场景扩展”的。场景数量会随着迭代次数增加而增加代码实现上更适合直接在循环里动态创建变量和约束而不是在模型初始化时写死。4.3 第二阶段子问题求解与 cut 生成子问题的核心是“给定 (x)找最坏 (u)”。用 Pyomo 建立子问题模型里有一类常见操作——把第一阶段变量作为参数传入。代码结构示意如下def solve_subproblem(x_values, param): sp_model pyo.ConcreteModel() sp_model.T pyo.Set(initializerange(param[T])) # 第二阶段变量各时段的调整量 sp_model.delta_p pyo.Var(sp_model.T, withinpyo.Reals) sp_model.delta_buy pyo.Var(sp_model.T, withinpyo.NonNegativeReals) sp_model.delta_sell pyo.Var(sp_model.T, withinpyo.NonNegativeReals) # 不确定量上下边界内连续 sp_model.w pyo.Var(sp_model.T, withinpyo.Reals, bounds(-param[w_max][t], param[w_max][t])) # 目标函数对偶化处理后是 max 形式 # 直接求解双层的 min 无法被求解器识别时先做对偶才是正路 ...如果你用的是 Gurobi还有一个更省事的做法直接在双层结构中利用 Gurobi 的参数化分析Parametric Analysis或者写显式枚举的循环。在我实际封装源码时遇到不确定时段数不多的小规模场景代码直接枚举几种极端“净负荷最差值”组合效果又快又稳。但一旦时段超过 6 个枚举路径就费劲了还是走对偶路径。4.4 CCG 主循环的收敛判断CCG 的收敛判据通常是上下界间隙[ UB - LB \le \varepsilon ]我的工程代码里主循环这样写UB float(inf) LB -float(inf) epsilon 1e-3 max_iter 20 for k in range(max_iter): # 1. 求解主问题得到 x*, theta* master.solve() LB max(LB, master.obj()) # 2. 固定 x*求解子问题得到最坏场景 u* 和第二阶段成本 Q(x*) sp_objective solve_subproblem(x_values, param) # 3. 用当前一阶段成本 最坏场景成本更新上界 UB min(UB, first_stage_cost_from_master sp_objective) # 4. 如果场景是新的则把该场景和第二阶段变量加入主问题 add_scenario_to_master(master, worst_scene) if UB - LB epsilon: print(fconverged at iter{k}) break一个容易掉进去的坑主问题里的 (\theta) 是下界变量迭代时应当不断通过新增场景约束把它“逼”向真实值。如果不小心把 (\theta) 初始化为 0并且子问题又恰好找到一个小于当前 LB 的组合那么上下界可能不会按预期收敛。解决方案是在主问题目标里给 (\theta) 一个足够大的初始系数或直接让 (\theta) 没有下界然后设初始值。5. 实验结果鲁棒方案与确定性方案的真实差距5.1 测试场景设定我使用 24 时段微电网实测数据光伏峰值 400 kW负荷峰值 600 kW储能 800 kWh两台可控机组容量分别 200 kW 和 150 kW。不确定性设置如下光伏波动±20%风电波动±20%负荷波动±10%购电价格峰 1.2 元/kWh平 0.8 元/kWh谷 0.5 元/kWh同时跑确定性模型直接用预测值和两阶段鲁棒模型最后用蒙特卡洛抽样 500 个随机场景分别代入两种调度方案做校验。5.2 核心指标对比指标确定性方案两阶段鲁棒方案调度计划总成本确定性视角3 682 元4 086 元蒙特卡洛平均实际运行成本4 175 元4 102 元最坏场景运行成本5 830 元4 520 元场景约束违反次数86 次0 次从表里能直观看出确定性方案在预测值下看起来很便宜但真实场景成本平均高出 13%最坏场景更是比鲁棒方案高出近 1300 元。鲁棒方案相当于用 11% 的“额外保险金”换来了完全没有约束违反的稳定性。对于离网型海岛微电网或对可靠性要求极高的园区这笔保费非常划算。5.3 不确定性预算的影响如果给不确定集额外加上预算约束比如所有时段里累计波动不超过某个上限那么模型能从“最保守”变得更从容。预算系数 (\Gamma) 从 0 到 24 变化时(\Gamma0)退化为确定性模型成本最低鲁棒性最弱(\Gamma12)中等策略成本比全保守低约 8%而最坏情况仍满足约束(\Gamma24)全保守最坏场景成本最高但安全性拉满。源码包里我把 (\Gamma) 作为参数开放出来你可以在脚本里直接改gamma_budget 12进行不同偏好下的权衡仿真。6. 常见坑与调试经验6.1 子问题出现不可行往往不是建模问题而是边界问题第一阶段决策如果过于激进例如把机组全部停机且储能已耗尽到了第二阶段面对极端净负荷场景可能没有任何一台设备能补上功率缺口子问题就会报 infeasible。常规的修法有几个给第二阶段加入虚拟切负荷和虚拟弃电变量并设置高额惩罚成本。这样子问题总能求解而惩罚项在目标函数里会让主问题主动避开不可行决策。检查不确定参数的边界是否过宽。±30% 的边界放在 24 时段全时段连续波动上几乎必然导致所有机组全部开启才能满足平衡每轮迭代子问题都给出极端值算法收敛缓慢。别忘了储能 SOC 的跨时段耦合。第二阶段补救虽然只是短时动作但如果 SOC 在边界卡死补救能力会骤降。6.2 对偶转换时符号方向最容易搞反对偶这一步是整个求解链路里最容易被忽略但最致命的地方。很多初学者把原问题写成 min 形式对偶后却按 max 的符号习惯推导导致 (\lambda) 的符号和约束方向完全对不上。写代码前我建议先取一个只有 2 个时段、单一机组的小例子手算一遍对偶再用求解器验算一遍。一个小技巧是原问题如果是 min对偶问题目标里的常数列要跟约束右侧同号不能想当然。6.3 求解性能优化三件小事第一把第二阶段目标函数里的惩罚系数调成有区分度的数值比如 1.5 倍燃料成本避免多个解之间成本差别太小导致求解器反复横跳。第二给主问题设置合理的求解时间上限和 MIP gap而不是默认求解到最优尤其当迭代到第 10 轮以后主问题场景变量数量膨胀MILP 求解时间可能从几十秒跳到十几分钟。第三充分利用热启动每轮主问题求解的把上一轮解作为 MIP start可以显著缩短后续求解时间。我有一次在 48 时段大系统上跑完整 CCG前 5 轮每轮主问题耗时不到 20 秒到第 9 轮因为场景数量膨胀直接卡到 6 分钟。后来加了mip_startlast_solution和mip_gap0.01总时长从 50 分钟压到 18 分钟效果立竿见影。6.4 代码注释是“隐藏的说明书”我平时有个习惯所有开源工程里的每一段代码注释都会写上“变量含义 单位 约束来源公式编号”调错的时候直接按图索骥不用翻几十页论文找对应关系。这个工程包里的所有注释也按这个标准来写。很多用户反馈说只看注释就能把模型重写一遍这就是我想要的。最后说点心得体会两阶段鲁棒优化在微电网调度里不是一个“纯理论炫技”的算法它的价值在于让调度方案天然地带上了安全冗余而 CCG 又把这种复杂嵌套优化拆分成了可以迭代求解的主子问题工程落地的门槛比我预想低很多。如果你之前一直用的确定性 MILP可以从最简单的盒式不确定集 枚举子问题开始把 CCG 写通之后再逐步升级不确定性模型、引入预算约束、接入真实气象预测数据。这个项目的源码里有完整的模型文件、样例数据、详细注释和输出模块你把它当成一个“可拆解的学习平台”就行。我第一次跑出上下界收敛曲线时那种把最坏情况“锁死”在可控范围内的感觉到现在都觉得特别踏实。本文还有配套的精品资源点击获取
简介微电网经济调度面临光伏、风电等可再生能源出力的强不确定性传统确定性优化在预测偏差下易导致弃光、切负荷等问题。鲁棒优化通过构建不确定集刻画最坏场景以两阶段决策框架实现“先决策、后补救”其中列与约束生成CCG算法将含max-min的复杂问题分解为主问题与子问题迭代求解有效平衡经济性与鲁棒性。基于Python与Gurobi/Pyomo的工程化实现使得盒式不确定集下的鲁棒调度模型可快速部署于实际微电网系统帮助运行人员在面对极端场景时仍能维持供电可靠性与经济性为园区微电网或海岛微电网提供稳健的调度策略。 做微电网经济调度的人基本都遇到过这种情况光伏预测曲线写的是中午能出 300 kW结果当天一片云飘过来直接掉到 80 kW下午负荷却比预计高了 10%。如果按确定性模型排好的机组出力、储能充放电计划去执行轻则弃光/切负荷重则频率波动甚至触发保护动作。我早期用确定性优化做日内调度时靠的是“预测误差留备用”这种粗放手段结果备用留少了要挨罚留多了又白花钱——这本质上就是在不确定性里赌运气。所以后来我把目光投向了两阶段鲁棒优化。它不赌预测值,而是把所有可能的不确定场景光伏、风电、负荷的波动区间全部摆到台面上目标函数在最坏情况下依然经济、可行。这篇文章就基于我近期整理的一个 Python 完整工程把微电网两阶段鲁棒经济调度的模型原理、CCG列与约束生成求解机制、Python 代码结构以及调试中踩过的坑全部拆开讲一遍直接把源码放到工程包里你拿到手改个数据就能跑。1. 从确定性调度到两阶段鲁棒优化为什么必须考虑“最坏情况”1.1 确定性经济调度到底缺了什么传统微电网经济调度的目标很简单在满足负荷、机组出力上下限、储能 SOC 约束的前提下把总运行成本降到最低。通常用混合整数线性规划MILP就能建出来决策变量机组启停状态、各时段有功出力、储能充放电功率、与配网购售电功率目标函数燃料成本 启停成本 购电成本 − 售电收益约束功率平衡、机组爬坡、储能能量状态、联络线容量。写成我常用的数学形式确定性模型长这样[ \min_{x,y} ; c^T x d^T y ][ \text{s.t. } Ax By \ge b,\quad y \ge 0 ]其中 (x) 表示 0-1 变量机组启停(y) 表示连续变量出力、储能量等。这个模型的求解难度不高可问题在于负荷、光伏、风电全部当作已知参数模型算出来的方案一旦遇到与预测不符的实际场景约束立刻可能被破坏。1.2 不确定性建模的核心思想先决策后补救两阶段鲁棒优化的思路跟确定性模型完全不同。它把策略拆成两个阶段第一阶段here-and-now在不确定性实现之前就要定下的决策。例如机组开几台、储能是否处于可调度状态、是否签约购电容量。这些决策一旦定下很难在短时间内改变。第二阶段wait-and-see不确定性实际发生之后调度员有权利做出一系列“补救动作”。例如多开一台机、调整储能放电功率、临时多购入电力等。目标从“最适合预测值”变成了“在最坏不确定性场景下总成本依然最低”。用数学语言就是[ \min_{x \in \mathbb{X}} \left( c^T x \max_{u \in \mathcal{U}} \min_{y \in \mathcal{Y}(x,u)} d^T y \right) ]这个 max-min 嵌套结构是整个问题的核心属于典型的 NP-hard 问题但工程上可以用 CCG 这种迭代算法高效逼近最优解。1.3 为什么盒式不确定集在工程里最常用不确定集 (\mathcal{U}) 有很多种定义方式常见的有盒式box、椭球、预算、数据驱动等。微电网工程里绝大多数代码用的是盒式集合因为简单、直观、不依赖太多历史数据的概率分布假设而且跟鲁棒优化的强对偶转化配合得非常干净[ \mathcal{U} { u \in \mathbb{R}^{n_d} : u_i^{min} \le u_i \le u_i^{max}, ; i1,\dots,n_d } ]这里的”不确定参数“在微电网里主要集中在三个地方光伏出力、风电出力和负荷功率还可以把与上级电网的交换功率也纳入。盒式集合的好处是只要给定每个时段预测值上下浮动比例就能快速构造可行域运算量和集合复杂度都可控。2. 微电网鲁棒调度模型的数学表达与约束拆解2.1 目标函数一阶段成本 二阶段再调度成本用下标 (t) 表示时段(i) 表示机组整个模型的目标函数可以写成[ \min_{\substack{u_{it}, \alpha_t \ P_{it}^g, P_{t}^{b}, P_{t}^{b-}}} ; \underbrace{\sum_{t,i} \left( a_i P_{it}^g b_i u_{it} C_i^{su} v_{it} \right) \sum_t \gamma_t^ P_t^{b}}{\text{第一阶段成本}} ; ; \underbrace{\max{u_d \in \mathcal{U}} \min_{ \Delta P } ; \sum_{t,i} r_i \Delta P_{it}^{g} s_i \Delta P_{it}^{g-}}_{\text{第二阶段补救成本}} ]这里 (u_{it}) 是机组启停 0-1 变量(P_{it}^g) 是机组出力(P_t^{b}/P_t^{b-}) 分别是购电/售电功率(\Delta P) 是第二阶段对出力的调整量(r_i, s_i) 分别是上调/下调的惩罚成本。读者看到的很多开源代码里第二阶段往往只写“再调度成本最小化”各项系数取正数原因是用它衡量最坏情况下的运行费用。如果第二阶段成本系数写得太小求解器会倾向于在一阶段做激进的低成本决策然后在第二阶段用大量调节来兜底得到的结果会失去鲁棒性。2.2 第一阶段约束与决策逻辑第一阶段约束主要包括功率平衡约束[ \sum_i P_{it}^g P_t^{wt} P_t^{pv} P_t^{b} P_t^{dis} D_t P_t^{b-} P_t^{ch} ]其中 (D_t) 是负荷预测值(P_t^{ch}, P_t^{dis}) 是储能充放电功率。这里先把不确定量写成预测值实际波动放到第二阶段约束里处理。机组出力上下限与爬坡约束[ u_{it} P_i^{min} \le P_{it}^g \le u_{it} P_i^{max} ][ -P_i^{ramp} \le P_{it}^g - P_{i,t-1}^g \le P_i^{ramp} ]储能动态约束储能的核心约束是 SOC 递推[ SOC_{t1} SOC_t \eta^{ch} P_t^{ch} - \frac{P_t^{dis}}{\eta^{dis}} ]再加上充放电功率上下限、SOC 上下限以及避免同时充放电的互补约束。工程实现里为了避免二进制变量爆炸有不少代码直接用 (SOC) 和可放电功率作为状态变量建模不过这会牺牲一点精度我在工程包里保留的是标准二进制约束版本。2.3 第二阶段约束把不确定性“放进来”第二阶段的本质是给定第一阶段决策 ((u_{it}, P_{it}^g))在每一个不确定性时段的实际出力/负荷取值下找一套可行的补救措施。为简化推导很多源码里把不确定性集中成不确定功率 (w_t)即净负荷波动[ D_t^{real} - (P_t^{pv,real} P_t^{wt,real}) D_t - (P_t^{pv} P_t^{wt}) w_t ]其中 (w_t \in [w_t^{min}, w_t^{max}])。同时第二阶段约束要求在加入补救变量 (\Delta P) 后功率平衡再次成立[ \sum_i \Delta P_{it} \Delta P^{b}_t - \Delta P^{b-}_t w_t ]并且所有补救变量都要受到爬坡约束、线路容量约束、储能剩余能量约束的限制。3. 列与约束生成CCG算法迭代逻辑3.1 主问题MP的作用收敛最坏情况下的成本下界CCG 的基本思路是把两层问题变成“主问题—子问题”的迭代。主问题是一个包含有限个场景的 MILP[ \min_{x, y_k} ; c^T x \theta ][ \text{s.t. } Ax By_k \ge b_k, \quad \forall k \in \mathcal{K} ][ \theta \ge d^T y_k, \quad \forall k \in \mathcal{K} ]每轮迭代都会把上一步计算出来的“最坏场景 (u^*)”以具体的常数代入主问题生成一组新的变量 (y_k) 和对应约束这叫“列与约束生成”。主问题求解后得到的目标值就是原问题的最优下界LB。说到底主问题干的活就是把不确定性从集合里抽出来变成一个个具体场景然后在这些场景同时满足的前提下做经济调度。每多一轮迭代主问题约束多一组LB 只升不降。3.2 子问题SP的作用找到最坏场景子问题要回答的是在第一阶段决策给定时找哪个不确定性场景 (u) 会让总成本最高从而迫使主问题做更保守的决策。形式化写出来是[ \max_{u \in \mathcal{U}} \min_{y} ; d^T y ][ \text{s.t. } Wy \ge h - Tx - Pu,\quad y \ge 0 ]内层是个线性规划外层是 max。处理办法是先对内层取对偶把 max-min 问题转化成 max 问题[ \max_{u, \lambda} ; \lambda^T (h - Tx - Pu) ][ \text{s.t. } W^T \lambda \le d,\quad \lambda \ge 0,\quad u \in \mathcal{U} ]这时问题里出现了一个双线性项 (-\lambda^T P u)。常规做法是用 KKT/线性化技巧处理或者直接利用盒式集合顶点在边界取最优的性质把问题进一步化解。3.3 顶点枚举 vs 对偶线性化两种主流实现方式我在工程包里同时实现了两种子问题求解方式路径一顶点枚举法。由于盒式集合的可行域是多面体线性目标的最优解必然在顶点取得。所以可以直接枚举每个不确定参数的上/下边界组合对每个组合求解一次线性规划取最坏值。这个方法的优点是写起来直白不易错缺点是场景数量随不确定参数个数指数增长。如果只对光伏、风电、负荷做单时段 24 个波动区间那组合高达 (2^{72}) 根本不可能枚举所以一般只在减少不确定参数个数如只聚合总净负荷时用。路径二对偶 大 M 线性化。引入辅助变量 (z \lambda^T P)通过大 M 将双线性项线性化。虽然会增加不少约束但整个子问题变成一个可以直接交给 Gurobi/CBC 求解的单层 MILP。这种方式在工程源码里更通用也是我推荐你在实际项目里用的版本。4. Python 与 Gurobi/Pyomo 实战从模型搭建到 CCG 主循环4.1 数据准备与全局参数我建议工程文件的目录结构按照“数据—模型—算法—输出”四条线来组织最省心的格式是 CSV 或字典。下面是一份微电网测试系统的简化数据示例参数数值单位机组数量2台时段数24h光伏预测峰值400kW风机预测峰值200kW负荷预测峰值600kW不确定性波动比例光伏/风电/负荷±20% / ±20% / ±10%-储能容量800kWh储能功率上限200kW购电价格分时峰/平/谷1.2/0.8/0.5元/kWh4.2 用 Pyomo 搭建第一阶段 MILPPyomo 是最适合做教学和科研复现的建模工具。你只要把目标函数和约束按数学表达式几乎 1:1 翻译过来即可import pyomo.environ as pyo def build_master_problem(uncertainty_scenarios, param): model pyo.ConcreteModel() # 时段与机组索引 model.T pyo.Set(initializerange(param[T])) model.I pyo.Set(initializerange(param[n_unit])) model.x pyo.Var(model.I, model.T, withinpyo.Binary) # 启停 model.pg pyo.Var(model.I, model.T, withinpyo.NonNegativeReals) model.pbuy pyo.Var(model.T, withinpyo.NonNegativeReals) model.psell pyo.Var(model.T, withinpyo.NonNegativeReals) model.theta pyo.Var(withinpyo.NonNegativeReals) # 第二阶段成本代理 # 目标函数一阶段成本 theta def obj_rule(m): return (sum(m.pg[i, t] * param[fuel_cost][i] m.x[i, t] * param[no_load_cost][i] for i in m.I for t in m.T) sum(param[buy_price][t] * m.pbuy[t] - param[sell_price][t] * m.psell[t] for t in m.T) m.theta) model.obj pyo.Objective(ruleobj_rule, sensepyo.minimize) # 功率平衡在各添加场景下逐一考虑…… return model要点在于主问题里的变量 (y_k) 是“按场景扩展”的。场景数量会随着迭代次数增加而增加代码实现上更适合直接在循环里动态创建变量和约束而不是在模型初始化时写死。4.3 第二阶段子问题求解与 cut 生成子问题的核心是“给定 (x)找最坏 (u)”。用 Pyomo 建立子问题模型里有一类常见操作——把第一阶段变量作为参数传入。代码结构示意如下def solve_subproblem(x_values, param): sp_model pyo.ConcreteModel() sp_model.T pyo.Set(initializerange(param[T])) # 第二阶段变量各时段的调整量 sp_model.delta_p pyo.Var(sp_model.T, withinpyo.Reals) sp_model.delta_buy pyo.Var(sp_model.T, withinpyo.NonNegativeReals) sp_model.delta_sell pyo.Var(sp_model.T, withinpyo.NonNegativeReals) # 不确定量上下边界内连续 sp_model.w pyo.Var(sp_model.T, withinpyo.Reals, bounds(-param[w_max][t], param[w_max][t])) # 目标函数对偶化处理后是 max 形式 # 直接求解双层的 min 无法被求解器识别时先做对偶才是正路 ...如果你用的是 Gurobi还有一个更省事的做法直接在双层结构中利用 Gurobi 的参数化分析Parametric Analysis或者写显式枚举的循环。在我实际封装源码时遇到不确定时段数不多的小规模场景代码直接枚举几种极端“净负荷最差值”组合效果又快又稳。但一旦时段超过 6 个枚举路径就费劲了还是走对偶路径。4.4 CCG 主循环的收敛判断CCG 的收敛判据通常是上下界间隙[ UB - LB \le \varepsilon ]我的工程代码里主循环这样写UB float(inf) LB -float(inf) epsilon 1e-3 max_iter 20 for k in range(max_iter): # 1. 求解主问题得到 x*, theta* master.solve() LB max(LB, master.obj()) # 2. 固定 x*求解子问题得到最坏场景 u* 和第二阶段成本 Q(x*) sp_objective solve_subproblem(x_values, param) # 3. 用当前一阶段成本 最坏场景成本更新上界 UB min(UB, first_stage_cost_from_master sp_objective) # 4. 如果场景是新的则把该场景和第二阶段变量加入主问题 add_scenario_to_master(master, worst_scene) if UB - LB epsilon: print(fconverged at iter{k}) break一个容易掉进去的坑主问题里的 (\theta) 是下界变量迭代时应当不断通过新增场景约束把它“逼”向真实值。如果不小心把 (\theta) 初始化为 0并且子问题又恰好找到一个小于当前 LB 的组合那么上下界可能不会按预期收敛。解决方案是在主问题目标里给 (\theta) 一个足够大的初始系数或直接让 (\theta) 没有下界然后设初始值。5. 实验结果鲁棒方案与确定性方案的真实差距5.1 测试场景设定我使用 24 时段微电网实测数据光伏峰值 400 kW负荷峰值 600 kW储能 800 kWh两台可控机组容量分别 200 kW 和 150 kW。不确定性设置如下光伏波动±20%风电波动±20%负荷波动±10%购电价格峰 1.2 元/kWh平 0.8 元/kWh谷 0.5 元/kWh同时跑确定性模型直接用预测值和两阶段鲁棒模型最后用蒙特卡洛抽样 500 个随机场景分别代入两种调度方案做校验。5.2 核心指标对比指标确定性方案两阶段鲁棒方案调度计划总成本确定性视角3 682 元4 086 元蒙特卡洛平均实际运行成本4 175 元4 102 元最坏场景运行成本5 830 元4 520 元场景约束违反次数86 次0 次从表里能直观看出确定性方案在预测值下看起来很便宜但真实场景成本平均高出 13%最坏场景更是比鲁棒方案高出近 1300 元。鲁棒方案相当于用 11% 的“额外保险金”换来了完全没有约束违反的稳定性。对于离网型海岛微电网或对可靠性要求极高的园区这笔保费非常划算。5.3 不确定性预算的影响如果给不确定集额外加上预算约束比如所有时段里累计波动不超过某个上限那么模型能从“最保守”变得更从容。预算系数 (\Gamma) 从 0 到 24 变化时(\Gamma0)退化为确定性模型成本最低鲁棒性最弱(\Gamma12)中等策略成本比全保守低约 8%而最坏情况仍满足约束(\Gamma24)全保守最坏场景成本最高但安全性拉满。源码包里我把 (\Gamma) 作为参数开放出来你可以在脚本里直接改gamma_budget 12进行不同偏好下的权衡仿真。6. 常见坑与调试经验6.1 子问题出现不可行往往不是建模问题而是边界问题第一阶段决策如果过于激进例如把机组全部停机且储能已耗尽到了第二阶段面对极端净负荷场景可能没有任何一台设备能补上功率缺口子问题就会报 infeasible。常规的修法有几个给第二阶段加入虚拟切负荷和虚拟弃电变量并设置高额惩罚成本。这样子问题总能求解而惩罚项在目标函数里会让主问题主动避开不可行决策。检查不确定参数的边界是否过宽。±30% 的边界放在 24 时段全时段连续波动上几乎必然导致所有机组全部开启才能满足平衡每轮迭代子问题都给出极端值算法收敛缓慢。别忘了储能 SOC 的跨时段耦合。第二阶段补救虽然只是短时动作但如果 SOC 在边界卡死补救能力会骤降。6.2 对偶转换时符号方向最容易搞反对偶这一步是整个求解链路里最容易被忽略但最致命的地方。很多初学者把原问题写成 min 形式对偶后却按 max 的符号习惯推导导致 (\lambda) 的符号和约束方向完全对不上。写代码前我建议先取一个只有 2 个时段、单一机组的小例子手算一遍对偶再用求解器验算一遍。一个小技巧是原问题如果是 min对偶问题目标里的常数列要跟约束右侧同号不能想当然。6.3 求解性能优化三件小事第一把第二阶段目标函数里的惩罚系数调成有区分度的数值比如 1.5 倍燃料成本避免多个解之间成本差别太小导致求解器反复横跳。第二给主问题设置合理的求解时间上限和 MIP gap而不是默认求解到最优尤其当迭代到第 10 轮以后主问题场景变量数量膨胀MILP 求解时间可能从几十秒跳到十几分钟。第三充分利用热启动每轮主问题求解的把上一轮解作为 MIP start可以显著缩短后续求解时间。我有一次在 48 时段大系统上跑完整 CCG前 5 轮每轮主问题耗时不到 20 秒到第 9 轮因为场景数量膨胀直接卡到 6 分钟。后来加了mip_startlast_solution和mip_gap0.01总时长从 50 分钟压到 18 分钟效果立竿见影。6.4 代码注释是“隐藏的说明书”我平时有个习惯所有开源工程里的每一段代码注释都会写上“变量含义 单位 约束来源公式编号”调错的时候直接按图索骥不用翻几十页论文找对应关系。这个工程包里的所有注释也按这个标准来写。很多用户反馈说只看注释就能把模型重写一遍这就是我想要的。最后说点心得体会两阶段鲁棒优化在微电网调度里不是一个“纯理论炫技”的算法它的价值在于让调度方案天然地带上了安全冗余而 CCG 又把这种复杂嵌套优化拆分成了可以迭代求解的主子问题工程落地的门槛比我预想低很多。如果你之前一直用的确定性 MILP可以从最简单的盒式不确定集 枚举子问题开始把 CCG 写通之后再逐步升级不确定性模型、引入预算约束、接入真实气象预测数据。这个项目的源码里有完整的模型文件、样例数据、详细注释和输出模块你把它当成一个“可拆解的学习平台”就行。我第一次跑出上下界收敛曲线时那种把最坏情况“锁死”在可控范围内的感觉到现在都觉得特别踏实。本文还有配套的精品资源点击获取
平方复杂度无法被常规优化消除｜TaoToken 统一 Key 通道实测)
