两阶段鲁棒优化在微电网调度中的建模与CCG算法实现

简介：本资源是一套面向电力系统优化方向研究生与科研人员的微电网两阶段鲁棒经济调度完整实现方案，聚焦解决含不确定性（如风电出力波动）下的调度保守性与经济性平衡问题。压缩包共13个文件，含4个核心MATLAB脚本（MP.m、SP.m、main_1.m、MP2.m）、3份PDF文献（含经典论文《微电网两阶段鲁棒优化经济调度方法》及拓展研究）、6张关键结果图（如目标函数收敛曲线、场景调度对比图等），总大小1.57MB；代码基于YALMIP建模，调用CPLEX/Gurobi求解，注释详尽，目标函数与约束均以紧凑矩阵形式表达，支持快速修改鲁棒调节系数以调控保守程度。已有3182人学习下载，代码非通用模板，为原创高复现度实现，收敛稳定、结构清晰、模块解耦良好，特别适合作为两阶段鲁棒优化入门范例，并可便捷迁移至机组组合、综合能源系统等同类问题建模。 做微电网调度的人，十有八九都遇到过这种情况：光伏预测明天中午功率特别高，结果当天一片云过来，发电量直接砍半；负荷预测跑了一堆深度学习模型，实际运行还是偏差好几个百分点。如果调度方案是拿预测值硬算出来的，遇到这种偏差就只能临时切负荷或者弃风弃光，运行成本蹭蹭往上涨。之前我做的确定性经济调度在理想条件下跑得确实飞快，但拿到实际运行数据一对比，问题就出来了——方案没有为不确定性留任何裕量。后来我把两阶段鲁棒优化方法完整复现了一遍，用matlab做平台、yalmip搭模型，cplex和gurobi都试过，把从模型推导到代码实现的整条链路都走通了。这篇博文就按复现的顺序，把两阶段鲁棒优化的建模思路、CCG列与约束生成算法的实现细节、yalmip调用求解器的常见坑，以及最后的算例分析一次性讲清楚。不管是正在做毕业设计的同学，还是刚接触鲁棒优化的工程师，这篇内容都能帮你少走不少弯路。

1. 项目整体设计与思路拆解

1.1 为什么确定性优化在微电网调度中不够用

先明确一个前提：微电网经济调度是一个“先决策、后运行”的问题。日前调度要在前一天就定下机组的启停、联络线功率、储能充放电计划，但真正运行的时候，光伏和风电的实际出力跟预测值总会有偏差。确定性优化把预测值当成真值来求解，得到一组“最优”方案，可一旦实际值偏离预测值，方案就可能不可行——比如计划里光伏要出100kW，实际只有50kW，剩下50kW没着落，只能快速切负荷。

有人可能会想：把预测值乘一个安全系数不就行了。这确实是工程上最常见的粗暴做法，但问题在于安全系数是拍脑袋定的，系统不同、季节不同、预测精度不同，一刀切的裕量要么太保守浪费钱，要么不够用还是出问题。确定性优化的本质是“赌预测会准”，而微电网场景下这个赌注风险太高。

1.2 两阶段鲁棒优化的核心思路与场景定位

鲁棒优化的思路完全不同：它不赌预测准，而是假设不确定量一定会在某个范围内波动，并保证即便最坏情况发生，调度方案也能满足运行约束。这样得到的方案不一定是最经济的，但一定是“安全”的。

但单阶段鲁棒优化有个臭名昭著的问题——太保守。因为它要求所有决策都能容忍最坏情况，而实际上很多决策是可以等不确定性实现后再调整的。这就引出两阶段结构：

• 第一阶段决策：必须在知道实际风光出力之前确定，比如机组启停状态、联络线购电协议；
• 第二阶段决策：等不确定量实现后，在运行中实时调整，比如机组出力、储能充放电、需求响应。

用生活化的话来说，这就像提前订酒店和当天灵活改签。第一阶段把房间订好，第二阶段看到现场情况后做调整。两阶段鲁棒优化要解决的正是“如何在保留调整空间的前提下，让总成本最低，同时保证最坏情况下不违约”。这个结构能覆盖微电网日前-日内两级调度、虚拟电厂聚合调度、园区综合能源系统等多类场景，适用范围很广。

1.3 为什么推荐matlab + yalmip + cplex/gurobi这套工具链

求解两阶段鲁棒优化不是直接扔给求解器就能出结果的，算起来需要反复迭代，每一步都有大量线性规划或混合整数线性规划要解。这时候工具链的选型就很重要了。

我用matlab做平台，原因是数据处理、画图、结果分析能一站式搞定，对电力系统的研究者来说也最熟悉。yalmip负责建模，它最大的好处是不锁定求解器——同样一套模型代码，想用cplex就用cplex，想用gurobi就用gurobi，底层切换只需要改一行设置。而cplex和gurobi是当前性能第一梯队的商用求解器，处理大规模MILP时比matlab自带的intlinprog快一个数量级以上。这套组合，基本是学术界和工业界做鲁棒优化的标配。

2. 数学模型拆解：两阶段鲁棒优化的核心推导

2.1 不确定性集合的选择与预算约束

两阶段鲁棒优化里，第一步就是定义“不确定性往哪个范围波动”。这里最常用的是“箱式 + 预算”的不确定集合：

U = { u | u_i^min ≤ u_i ≤ u_i^max, ∑ |(u_i - u_i^pre) / (u_i^max - u_i^pre)| ≤ Γ }

用大白话讲，第一组约束限定了每个不确定量的波动上下限；第二组约束限制了“总共有多少个变量可以同时达到最大偏差”——这就是预算参数Γ的作用。

为什么要加这个Γ？因为如果所有光伏、负荷都同时取最坏值，概率极低，方案会过于保守。Γ值越小，不确定性集合越小，方案越经济；Γ值越大，集合越大，方案越保守。实际使用中Γ是一个可调节旋钮，用来平衡经济性和鲁棒性，比如Γ取2表示最多两个不确定量同时到边界。

我推荐入门先选这种线性不确定性集合，因为它是线性的，对偶转换后不会引入非线性项，求解简单。相比之下，椭球不确定集合虽然更精确，但会引入二阶锥约束，对初学者来说调试难度陡增。

2.2 两阶段模型的一般化表示

两阶段鲁棒经济调度模型可以写成下面这个简洁的“min-max-min”结构：

min_x c^T x + max_{u∈U} min_{y∈Ω(x,u)} d^T y

• 第一阶段：min_x c^T x，x 是第一阶段的“待定决策”，对应机组启停、事前签约量等；
• 第二阶段：max_{u∈U} min_{y} d^T y，意思是“先找到最坏的不确定场景u，再在u下做最优的经济调度y”。

约束部分涉及功率平衡、机组出力上下限、爬坡约束、储能SOC递推、联络线容量等。这些约束本质上是线性约束，这也是我们能用对偶方法求解的前提。

一个常见的疑问是：第二阶段问题如果包含整数变量怎么办？比如储能充放电状态、需求响应分段状态。答案是：如果第二阶段必须含0-1变量，对偶方法就不再严格适用，需要改用KKT条件或者组合Benders等更复杂的手段。所以在入门复现时，建议先把储能建模成连续变量（充放电功率连续可调，SOC连续），先跑通整体流程，再逐步增加复杂度。

2.3 max-min对偶转换：从三层问题变成两层问题

现在的问题是内层的 min_y d^T y 是个线性规划，对于LP来说强对偶定理成立。于是我们可以把内层LP写成其对偶形式 max_λ(...)，整个max-min-min结构就变成了：

max_{u∈U} max_{λ∈Λ} (对偶目标函数)

两个max可以合并成一个max，问题从三层嵌套变成两层。但这个过程中会出现一个麻烦——双线性项。具体地说，不确定性变量u和对偶变量λ会相乘，比如u^T λ。这一项不是线性的，没法直接扔给cplex/gurobi求解。

解决办法是big-M线性化：引入辅助变量，把u^T λ这种双线性项拆成带大M的线性不等式组。M的取值不能取得太大，太大会造成数值病态；但也不能太小，太小会错误地截断可行域。实操中我会先用确定性情况试算一遍，观测对偶变量的量级，再定M的数量级，通常设成对偶变量量级的10到100倍。

另一个可选的路线是用KKT条件替换第二阶段。相比对偶方法，KKT的推导对初学者更直观，但互补松弛条件同样需要big-M线性化，而且需要补充大量辅助变量。我在复现时两种都试过，体感上对偶方法在CCG迭代中实现更简单，求解速度也更稳。

3. CCG算法求解流程与matlab/yalmip代码实现

3.1 列与约束生成（CCG）的核心迭代逻辑

模型建完之后，最常用的求解算法是列与约束生成（CCG）。它的思想是“不要一次性考虑所有不确定场景，而是在迭代中逐步找到最关键的场景”。CCG维护两个值：

• 下界值LB：来自主问题。主问题只包含“已经找到的那些最坏场景”，所以它得到的目标值只会更低，是全局最优的下界；
• 上界值UB：来自子问题。子问题固定了第一阶段决策，求解真实的最坏场景，所以得到的成本代表实际可达到的上界。

每次迭代过程：

1. 求解主问题，得到第一阶段决策x*和目标值，更新LB；
2. 固定x*，求解子问题（也就是max-min问题），找到当前最坏场景u*和目标值，更新UB；
3. 如果UB - LB小于给定阈值，算法收敛；
4. 否则把新找到的场景u*作为新约束加入主问题，回到步骤1。

这里有个容易误解的点：主问题的变量会随着迭代不断增多。每识别出一个新的最坏场景，主问题就追加一组第二阶段变量和约束。这也是CCG名字里“列生成”的由来——每一轮都在往主问题里加一列（一组新的决策变量）。

3.2 初始化与收敛判据的选取

初始化上界UB一般取正无穷，下界LB取负无穷。第一轮求解主问题时如果没加任何场景，模型可能无界，所以第一个场景一般用不确定量的预测值（或者让所有不确定量取区间中点）来初始化。

收敛判据建议用相对间隙而不是绝对值：abs(UB-LB)/abs(LB) ≤ 1e-3。因为如果成本量级是几百万，绝对间隙1块钱就非常严格，求解器很难达到，而相对间隙1e-3通常已经足够工程精度。实际复现中CCG迭代次数一般在3到8次之间就能收敛，这个收敛速度比Benders对偶切割要快不少。

3.3 matlab + yalmip的核心代码框架

下面给出核心代码框架，去掉细节数据定义，重点看逻辑：

% 定义不确定集合参数 u_pred = [...]; % 预测值 u_max = ...; u_min = ...; Gamma = 2; % 预算参数 % 主问题变量 x = binvar(n_x,1); % 第一阶段变量（如机组启停） y_mp = sdpvar(n_y, K, n_day); % 第二阶段变量，K为场景数 % ... 其他变量定义 Constraints_mp = []; obj_mp = c'*x; for k = 1:K Constraints_mp = [Constraints_mp, A*x + B*y_mp(:,k) >= b + E*u_scenario(:,k)]; obj_mp = obj_mp + d'*y_mp(:,k); end ops = sdpsettings('solver','gurobi','verbose',2,'gurobi.mipgap',1e-4); optimize(Constraints_mp, obj_mp, ops); x_hat = value(x); % 子问题：经过对偶和big-M线性化后的形式 lambda = sdpvar(n_con,1); % 对偶变量 % 双线性项 u'*lambda 用辅助变量和big-M约束线性化 Constraints_sp = [对偶可行域约束, big-M线性化约束]; obj_sp = -(对偶目标函数表达式); % 把max问题转成min求解 optimize(Constraints_sp, obj_sp, ops);

这段代码是简化的逻辑示意，实际的子问题还有big-M线性化的辅助变量、对偶可行域约束等。但整体框架就是这样：主问题不断累积场景，子问题用对偶形式解出最坏场景。

写代码时有三个细节值得注意：

1. 每次迭代主问题要完全重新求解，而不是在上一次解的基础上热启动。虽然热启动理论上更快，但新增场景会改变整个约束结构，热启动效果不稳定，我实测下来完整重解更省心。
2. 子问题里固定x_hat的方式是把x_hat作为数值代入，而不是作为等式约束。如果作为约束加入，会把二阶段问题变成两阶段耦合问题，求解器求解方式完全不同，容易出错。
3. 子问题求解结束后要检查求解器的退出标志，碰到“infeasible”要先查对偶转换是否有误，而不是盲目调参。这一步是调试鲁棒优化代码最高频的坑。

4. 求解器配置实操：yalmip调用cplex与gurobi的关键细节

4.1 yalmip的安装与型号设置

先把yalmip安装讲清楚：从官网下载最新zip包后解压，在matlab里addpath(genpath('yalmip目录'))，然后savepath保存路径。验证安装是否成功，在命令行跑一下yalmiptest，它会自动检测当前机器上可用的求解器。

yalmip调用求解器的设置核心是sdpsettings：

ops = sdpsettings('solver', 'gurobi'); ops = sdpsettings('solver', 'cplex');

切换求解器只需要改这一行。但要提醒的是，gurobi和cplex的MILP参数名略有不同，比如gurobi的mipgap在yalmip中要写成‘gurobi.mipgap’，cplex写成‘cplex.mip.tolerances.mipgap’。这些细节在yalmip的文档里都能查到，但新手经常卡在这一步。

4.2 cplex与gurobi的选择建议

两者都是商业求解器，但使用体验上差别不小。简单给个对比：

我的建议是：学术用途优先申请gurobi的免费license，使用体验更顺滑；但如果你所在课题组或公司已经买了cplex，那用cplex完全没问题，因为代码层面你只需要改设置那一行，不影响其他逻辑。

4.3 求解器参数调优：MIPGap、TimeLimit与数值稳定性

CCG迭代中主问题是一个混合整数线性规划，随着场景不断加入，变量和约束规模快速增长，求解时间会明显上升。调参的关键是：

1. MIPGap：默认情况gurobi会追求非常高的精度，但CCG迭代中其实不需要每一轮都精确到1e-6。把MIPGap设置为1e-4甚至5e-4，三到五轮迭代就能节省大量时间。
2. TimeLimit：给每一轮求解设置时间上限，比如60秒或120秒。防止某轮主问题卡死导致整个迭代死循环。
3. 数值容忍度：big-M处理过的模型往往带有数值病态问题，适当放宽求解器的可行容差（如gurobi的feasibilitytol设置到1e-6，不要默认的1e-7），能减少很多“fake infeasible”问题。

实测下来，同样一个18节点微电网算例，默认参数下主问题一轮可能要跑两分钟，把MIPGap放宽到1e-4后，一轮只要四十秒左右，最终结果几乎没差别。

5. 常见问题与排查技巧实录

在复现过程中我踩了不少坑，这里整理成表格，每个都很典型。

这里单独说两个值得展开的点。

第一，big-M大小要反复试。M取太小，双线性项被错误限制，求解出来的所谓“最坏场景”不是真正最坏；M取太大，数值条件数恶化，主问题求解会变得奇慢甚至出现病态结果。我的经验是先跑一次确定性调度，把所有对偶变量的量级打印出来，然后M取最大的对偶变量绝对值乘20左右，再微调。

第二，求解器的版本匹配是稳定性大敌。matlab 2022b这类新版本对旧版mex文件兼容性差，如果cplex或gurobi版本太旧，常常出现一调用就崩溃或报error 9。解决方案是去官网下载与matlab版本匹配的最新版求解器，然后用求解器自带的matlab接口重新编译。遇到这类问题不要反复重装matlab，先查求解器版本。

6. 复现结果分析与参数敏感性讨论

6.1 典型算例设置

为了验证方法，我搭了一个典型的小型微电网场景：

• 1台微型燃气轮机，额定功率100kW，爬坡速率10kW/min；
• 光伏电站50kW，风机40kW；
• 储能系统80kWh/50kW，SOC范围0.1-0.9；
• 基础负荷约80kW，允许5%的需求响应；
• 与配电网的联络线功率上限50kW。

风光和负荷的预测值由历史曲线生成，设置了±20%的波动区间。不确定性预算Γ取2。分别跑确定性调度、两阶段鲁棒优化两个方案。

6.2 关键结果对比

确定性调度在预测场景下总成本最低，比如一天成本约1800元。但把实际场景拉到最坏情况（光伏降20%、负荷升20%）后，确定性方案的功率平衡约束直接不满足，要么切负荷要么弃风弃光，额外惩罚成本可能超过300元，实际总成本冲到2100元以上。而两阶段鲁棒优化把最坏情况一开始就纳入考量，调度计划在全体可行场景中保持可行，总成本约1950元，比确定性方案在最坏场景下的实际成本低。

这就是鲁棒优化的真正价值：它不是让“理想情况”更省钱，而是让“糟糕情况”不至于失控。鲁棒方案在预测场景下会比确定性方案贵5%-10%，这部分“额外成本”可以理解为买保险的保费。

6.3 Γ参数的影响：保守性-经济性权衡

Γ从0（完全忽略相关性）逐渐增大到5，系统总成本呈现明显的单调上升趋势。Γ=0时等价于确定性优化，成本最低但鲁棒性最差；Γ=5时所有不确定量都同时取最坏值，成本最高但绝对安全。

实操建议是：不要直接把Γ拉满，而是画一条“成本-Γ”曲线，结合历史运行数据中实际出现的最坏偏差次数来选Γ。如果统计显示同时出现两个以上大幅度偏差的概率很低，就取Γ=2或3，这样既不会太保守，又能覆盖绝大多数情况。

我在复现时还发现一个有意思的现象：当Γ超过一定值后，成本上升的斜率会变缓。这是因为系统里燃气轮机和储能等调节资源有限，在极端场景下已经达到出力上限，裕量无法继续增加，多出来的不确定性只能用切负荷惩罚成本来兜底。这个拐点可以作为系统扩容决策的参考信号——如果实际运行中Γ经常需要取得很大，说明系统的调节资源确实不够。

整套复现下来，我最大的体会是三句话。第一，两阶段鲁棒优化的核心难点不在数学公式，而在CCG迭代里每一个子模块的衔接，任何一个环节的bug都会表现为“不收敛”或者“结果不合理”，所以一定要先跑通确定性模型，再加不确定性。第二，yalmip这套工具链是真的省心，建模和求解器解耦，换求解器只改一行配置，但前提是把版本匹配和许可证问题提前处理好。第三，不确定集合和预算参数的选取不是纯数学问题，而是对实际运行场景的建模决策，脱离数据拍脑袋选参数，跑出来的结果再漂亮也没有工程意义。最后再分享一个小技巧：把主问题和子问题分别封装成两个函数，加一段日志打印每一轮的LB、UB、间隙和最坏场景信息，调试效率立刻上一个台阶。希望这篇内容能帮正在做微电网调度和鲁棒优化的你省下几周时间。

本文还有配套的精品资源，点击获取