两阶段鲁棒优化这几年在微电网调度里确实是个绕不开的热点方向。我看不少朋友卡在“模型懂了、代码跑不通”这一步尤其是“关键场景辨别算法”这个切入点——它对应的是两阶段鲁棒求解中怎么加速找到最恶劣场景的问题。这篇就把我实际调试这套Matlab代码的过程、模型怎么搭、关键场景怎么识别、以及那些文档里不写但特别影响结果的细节一次性梳理清楚。1. 项目整体设计与思路拆解1.1 为什么选两阶段鲁棒而不是传统确定性优化微电网调度里最头疼的从来不是“已知明天负荷曲线然后排机组”而是风光出力、负荷波动这些不确定性怎么处理。我最早做微电网调度用的是确定性优化——给定一组预测值直接求解机组启停和出力计划。结果一到实际运行就露馅光伏预测偏差一大要么弃光弃得肉疼要么功率不平衡导致联络线越限。后来转向鲁棒优化。鲁棒优化的核心思想是“我不赌一个具体预测值而是设一个不确定集模型在这个集里所有可能场景下都必须可行”。这个思路很对但单阶段鲁棒往往太保守——它要求所有不确定参数同时取最坏值而现实里风电大、光伏未必小这些极端情况根本不会同时发生。两阶段鲁棒就是来修正这个问题的第一阶段做“现在必须定下来”的决策比如机组启停、与主网交互的购售电计划第二阶段在不确定性逐步揭示后做“可以调整”的决策比如蓄电池充放电、微燃机出力修正。用一句大白话说第一天先把房间隔断搭好第二天根据实际天气再决定每个房间怎么用。这样既保证了最坏情况下的安全性又不会因为过度保守白白浪费经济性。1.2 关键场景辨别算法在这个框架里的定位两阶段鲁棒模型的求解最常用的思路是CCG列与约束生成算法也就是把原问题拆成主问题和子问题不停迭代主问题给一个第一阶段的决策子问题在这个决策下寻找到最恶劣的不确定场景再把该场景对应的约束“喂”回主问题重新优化直到上下界收敛。这里有个工程痛点子问题搜最恶劣场景时不确定集里可能有成百上千个候选场景如果每个主问题迭代都把所有场景跑一遍计算量会非常感人。我实测过一个小型微网6节点、3台机组、含储能直接用穷举方式遍历场景池单次CCG迭代要6到8分钟整套算法收敛要十几轮跑完一次算例等得怀疑人生。关键场景辨别算法解决的正是这个问题它先通过某种规则或灵敏度分析从场景池里挑出“对目标函数影响最大、最能代表极端情形”的少数几个关键场景主问题只把这些关键场景对应的约束加进去。场景池是动态更新的——每次子问题求解后实际找到的最恶劣场景会被加入“已识别关键场景”集合下一轮迭代就不用再全量扫描。这个“先辨识、后加约束”的思路本质上是用很小的精度损失换来了数量级的求解速度提升工程实用价值很高。1.3 这套方案的适用范围和预期效果这套方法适合含风电/光伏、储能、微燃机并且与主网有功率交换的交流或直流微电网。你只需要修改基础数据负荷曲线、风光预测出力、机组参数即可迁移到不同算例。在我自己复现的案例里6节点微网48个时段的调度问题接入关键场景辨别后整体求解时间从原来的40多分钟压缩到了6分钟以内目标函数值和穷举法相比偏差控制在0.5%以内。对于需要反复跑敏感性分析的研究场景来说这个收益非常明显。2. 核心细节解析与实操要点2.1 不确定性集的选择与数学建模两阶段鲁棒的第一步就是把你关心的不确定参数构造成一个“集合”而不是一个点。最常见的是盒式不确定集[ U { \tilde{p} \in R^T \mid p_t^{fore} - \hat{p}_t \le \tilde{p}_t \le p_t^{fore} \hat{p}_t } ]其中 ( p_t^{fore} ) 是t时段预测值( \hat{p}_t ) 是偏差上限。你别小看这个简单的盒子它直接决定了鲁棒模型长什么样。如果你只用盒式集子问题会退化成“每个时段都取最坏值”这样两阶段模型又变得过度保守所以实际工程中几乎都会引入一个预算约束budget of uncertainty[ \sum_t \frac{|\tilde{p}_t - p_t^{fore}|}{\hat{p}_t} \le \Gamma ](\Gamma) 的含义是“最多允许几个时段同时达到极端偏差”。它的取值直接影响保守度(\Gamma0) 退化为确定性模型(\GammaT) 变回纯盒式鲁棒。我在项目里用蒙特卡罗模拟不同 (\Gamma) 下的成本和违反概率最终选了一个违约概率不超过5%的最大 (\Gamma)。这里要特别提醒一个新手容易踩的坑两个不确定源比如光伏和风电如果各自都建一个 (\Gamma)子问题求解时它们的取值是耦合的。实际项目里我采用的是统一的预算约束加权重系数把不同源的偏差归一化后再共享同一个(\Gamma)。这一处理让模型在数学上更干净求解时也不用担心两个“最坏”叠加导致过度保守。2.2 两阶段模型的目标函数与约束体系目标函数分两阶段写。第一阶段目标主要包含机组启停成本、与主网的购电成本/售电收益这些是在不确定性揭示前就要做好的决策。第二阶段目标是在给定第一阶段决策和最坏场景下调整微燃机出力、储能充放电使得调整成本最小。写成数学形式就是[ \min_{x} (c^T x \max_{u \in U} \min_{y} d^T y) ]其中 x 是一阶段决策u 是不确定参数y 是二阶段调整决策。注意这个 min-max-min 结构是整个模型的核心也是区别于单阶段优化的标志。约束体系我按这几类来列功率平衡约束这是必须严格满足的等式约束注意如果含直流潮流或交流潮流线性化处理方式不同。微燃机出力约束不仅要有上下限还要考虑爬坡约束。爬坡约束我一开始漏掉了导致输出的调度方案在实际执行时飞车结果做动态仿真直接失稳。储能约束包括充放电功率限制和SOC范围。这里有一个关键的建模技巧充电效率和放电效率如果不一致必须分开写约束不能简单用一个效率系数。否则模型会利用“同时充放电”来白嫖能量我见过不少代码栽在这上面。联络线功率约束与主网的交换功率限制注意购电和售电要分别处理不能用一个可正可负的变量代替否则目标函数里购电价和售电价不一样时模型会“买贵卖贱”来套利。2.3 两阶段模型的求解架构主问题与子问题的交互采用CCG算法来求解整体框架如下主问题MP——在已知部分关键场景的前提下最小化第一阶段成本加上第二阶段成本的近似同时满足第一阶段约束和“已经被识别出来的关键场景”对应的第二阶段可行性约束。主问题的解给出了原问题的下界如果是min问题。子问题SP——给定主问题求出的第一阶段决策 (x^)寻找使第二阶段成本最大的那个不确定场景 (u^)同时计算对应的第二阶段最小成本。子问题的解给出了原问题的上界。算法的停止判据是上下界之差小于某个阈值比如0.01。如果你发现上下界收敛很慢多半不是算法实现的问题而是子问题求解精度设置不当或者割平面添加方式有误——这个我在第5节会详细讲。3. 实操过程与核心环节实现3.1 系统建模与基础数据准备为了这篇讲解我搭了一个简化但不失代表性的微网算例包含1台微型燃气轮机额定功率1.6MW、1台柴油发电机1MW、一组储能容量2MWh功率0.5MW、一个光伏电站装机2MW、一个风电场装机1.5MW通过联络线与上级电网相连调度周期为24小时分辨率为1小时。基础数据我放在Excel里用Matlab的 readtable 读取这样方便你换算例。三个比较关键的数据集一是负荷预测曲线二是光伏出力预测和风电出力预测三是分时电价。注意这里的负荷预测不是确定值同样要叠加偏差否则构建不确定集时数据来源不全。负荷曲线我用了一个典型的夏季日负荷峰值出现在19:00-21:00大约1.8MW光伏预测取正午12:00达到峰值1.6MW夜间为零风电预测则是夜间出力大、白天出力小典型的反调峰特性。这样设置能让算例呈现出“白天光伏多了要储能充电晚间风电不足要靠燃气轮机顶上”的有趣调度逻辑。3.2 关键场景辨别算法的实现步骤关键场景辨别算法是这套代码的核心我把它拆成5步第一步场景池生成根据预测值和偏差范围随机生成N个候选场景。这里我用的是拉丁超立方采样而不是简单蒙特卡洛因为LHS能保证采样点均匀覆盖整个不确定空间用较少的样本数就能获得较好的代表性。第二步场景聚类预筛选对生成的场景做K-means聚类聚成K簇。聚类时特征向量就是各时段的不确定参数值对光伏场景就是24个时段的光伏出力值。每个簇的中心场景作为“候选关键场景”进入下一步筛选。这一步能挡掉大量重复的相似场景大幅降低后续辨识的计算负担。第三步灵敏度排序这一步是算法的精髓。对每个候选场景计算它对应的第二阶段成本对一阶段决策的灵敏度。怎么理解呢如果一个场景本身就很恶劣比如光伏出力远低预期、负荷又远高预期那它一定对成本影响巨大反之一个接近预测值的“温和”场景不会给系统带来额外压力。通过这个灵敏度指标排序我们只需要关注排名前m个场景。具体实现上我没有直接求梯度因为子问题求解本身就耗时而是用一个更聪明的办法先对所有候选场景跑一次子问题求解记录各自的目标函数值然后按目标函数值降序排列取前m个作为“预选关键场景”。这个做法的理论依据是在凸优化框架下目标函数值最大的场景就是在当前一阶段决策下最危险的场景。第四步主问题迭代修正主问题求解时只把“预选关键场景”所对应的约束加进去。求解完主问题后把新的一阶段决策代回子问题重新评估所有候选场景。如果某个原本排序靠后的场景在新决策下的目标函数值升到了前列就把它加入关键场景集合。这样形成一个“辨识—求解—再辨识”的动态闭环。第五步收敛判据子问题找到的“当前最恶劣场景”和主问题目标函数值之间的gap小于阈值时停止迭代。我用的判据是相对gap 1e-3同时严格监控最大迭代次数设置50次上限防止因为数值问题导致死循环。3.3 Matlab代码实现与Yalmip建模整体代码架构我分成五个模块主程序main.m、数据读取模块load_data.m、主问题建模build_MP.m、子问题建模build_SP.m、场景生成与辨别scenario_identification.m。这样分层的好处是便于调试——每个模块跑完可以单独检查不必每次从头跑。我用的求解器组合是Yalmip Cplex。Cplex求解混合整数线性规划是公认的稳而且Yalmip提供了非常友好的Matlab语法接口。如果你没有Cplex许可证Gurobi也可以替代只需把代码里的cplex改成gurobi即可。主问题的核心代码逻辑如下% 主问题最小化第一阶段成本 第二阶段成本基于已识别关键场景 x binvar(2, 24); % 机组启停状态2台机组 x 24小时 p sdpvar(2, 24); % 机组出力 soc sdpvar(1, 25); % 储能SOC25个点因为要包含初始和结束 u_grid sdpvar(1, 24); % 联络线功率正为购电负为售电 % 第二阶段变量针对每个关键场景都需要一套 y cell(1, n_key_scenarios); for k 1:n_key_scenarios y{k}.p2 sdpvar(2, 24); % 第二阶段机组调整出力 y{k}.soc2 sdpvar(1, 25); % 第二阶段SOC修正 end % 目标函数启动成本 燃料成本 购电成本 第二阶段调整成本 objective sum(sum(start_cost .* max(0, diff([zeros(2,1), x], 1, 2)))) ... sum(sum(fuel_cost .* p)) ... sum(price_buy .* max(0, u_grid)) - sum(price_sell .* max(0, -u_grid)) ... sum_over_scenarios(adjust_cost);注意这里用了一个Matlab的小技巧用 max(0, u_grid) 和 max(0, -u_grid) 来分别处理购电和售电这样在Yalmip里会自动引入辅助变量不需要手动构建额外的二进制变量。子问题的核心代码是求解 max-min 问题。标准做法是利用强对偶将内层min问题转化为max问题使整个子问题变成单层max问题function [worst_cost, worst_scenario] solve_SP(x_star) % x_star 是主问题传来的一阶段决策 % 内层min问题的对偶变量 lambda sdpvar(1, 24); % 功率平衡约束的对偶变量 mu_lb sdpvar(2, 24); % 机组出力下限约束的对偶变量 mu_ub sdpvar(2, 24); % 机组出力上限约束的对偶变量 % 对偶问题目标函数最大化 dual_obj sum(lambda .* (load - scenario_pv - scenario_wt)) ... - sum(mu_ub .* p_max) sum(mu_lb .* p_min); % 对偶约束 constraints []; for t 1:24 for i 1:2 constraints [constraints, lambda(t) - mu_ub(i,t) mu_lb(i,t) fuel_cost(i)]; end constraints [constraints, mu_lb(:,t) 0, mu_ub(:,t) 0]; end % 外层max还要遍历不确定参数 optimize(constraints, -dual_obj); % 注意Yalmip默认是min所以取负号 end这里有个特别容易出错的地方对偶问题里不确定参数和x_star是线性耦合的但它们不是优化变量而是参数。因此子问题其实是一个带不确定参数的最大化问题。实际求解时我用的是枚举关键候选场景 每个场景下求解一个线性规划的做法先固定不确定参数为某个候选场景然后求解对偶LP得到对应的目标值最后取最大的那个场景。这样比直接用Yalmip处理双线性项要稳健得多效率也高。3.4 参数设置与结果分析关键参数设置如下参数数值说明预测偏差比例光伏15%风电20%负荷10%偏差越大不确定性越强预算参数 (\Gamma)824时段中最多8个时段同时达到极端偏差场景池大小500初始随机场景数量K-means簇数20聚类后候选场景数预选关键场景数5每轮迭代加入主问题的场景数收敛阈值1e-3上下界相对gap仿真结果总成本约12,847元。其中第一阶段固定成本启停基础购电约9,623元第二阶段调整成本约3,224元。相比不考虑不确定性的确定性模型总成本约10,158元鲁棒模型付出了约26.5%的“保险溢价”——这个溢价就是系统在最坏场景下保持安全运行所需付出的代价。而相比纯盒式鲁棒总成本约15,632元关键场景辨别方法节省了约17.8%的成本充分说明它有效规避了过度保守的问题。4. 常见问题与排查技巧实录4.1 子问题对偶推导出现符号错误这是所有两阶段鲁棒代码里最隐蔽的坑。对偶问题推导时如果原问题是min且约束是≤那么对偶变量应满足≥0约束是≥对偶变量≤0等式约束的对偶变量无符号限制。我在第一次实现时就栽在这里——把等式约束对偶变量也加上了非负限制结果导致子问题上界始终偏低主问题和子问题之间的gap永远收敛不了。排查技巧先用一个很小的确定性用例手算一遍对偶推导再对照代码验证。如果发现上界异常偏低优先检查对偶变量的符号约束这个问题的出现频率占我调试时间的一半以上。4.2 CCG迭代不收敛或收敛极慢遇到的第二个典型问题是迭代了20多轮gap还在5%附近晃悠。排查后发现根因是主问题添加关键场景约束时遗漏了不确定参数的变化对功率平衡的影响。具体来说主问题中的功率平衡约束是针对“名义场景”建立的加入关键场景约束后没有同步更新功率平衡等式里的不确定参数项。解决方案在添加场景约束时把该场景对应的不确定参数值作为一个参数传入主问题重新生成该场景下的功率平衡约束和其他相关约束。这一步代码量不大但很容易漏。4.3 场景辨别算法误删关键场景还有一个工程细节在聚类预筛选阶段如果K值设得太小比如5某些“单独不恶劣但在组合下很恶劣”的场景可能被聚类中心抹平。我实测发现K20时结果最稳定K10时偶尔会出现收敛到次优解的情况。解决办法在灵敏度排序阶段我把“临近聚类中心的边界场景”也加入候选池而不是只取聚类中心。这样虽然增加了预选场景的计算量从K个增加到大约3K个但能有效避免漏掉极端组合。4.4 Yalmip求解器报错“Infeasible problem”这个报错多数时候不是真的不可行而是建模时变量定义域出了问题。最常见的两种一是SOC变量没有定义上下限求解器找不到解二是储能结束SOC约束设得太死比如必须回到初始值而系统实际能力达不到。我在项目里把结束SOC改成柔性约束——允许偏离初始值但目标函数里加一个惩罚项。这样既保证了储能的长期可用性又不会因为硬约束导致整体不可行。4.5 常见问题速查表现象可能原因解决方案子问题上界过低对偶变量符号错误手推小算例验证对偶推导主问题不可行一阶段决策约束过紧检查启停逻辑和功率平衡约束收敛慢关键场景集合太小增大预选场景数或聚类数K目标函数异常偏小储能同时充放电检查充放电效率是否分开建模结果波动大场景池随机性不足改用拉丁超立方采样增加场景数求解时间过长场景池过大未聚类调整K-means簇数先粗筛再精筛5. 代码运行环境与迁移复现建议5.1 运行环境配置我在Matlab R2022b环境下测试的这套代码需要安装以下工具箱和外部求解器Yalmip推荐从GitHub拉取最新版这不是官方工具箱但维护活跃Cplex 12.10 或 Gurobi 9.5学术许可免费Matlab Optimization ToolboxYalmip运行的基础依赖安装顺序建议先装Matlab再装Cplex或Gurobi最后把Yalmip添加到路径。注意Cplex安装时要选择“添加到Matlab路径”选项否则Yalmip找不到求解器。5.2 迁移到自己的微网算例换算例时需要修改的地方基本集中在 load_data.m 和参数设置区。具体来说修改机组参数表把微燃机和柴油机的容量、爬坡率、燃料成本系数换成你自己的参数。修改储能参数容量、最大充放电功率、充放电效率、SOC上下限。修改负荷和风光预测曲线Excel里直接替换即可注意列顺序要和读取代码匹配。调整不确定集参数预测偏差比例和预算参数(\Gamma)需要根据实际数据重新标定。有个细节要提醒如果你把调度周期改成96时段15分钟分辨率场景池大小和聚类簇数都需要相应增加。我测试过96时段案例500个场景30个簇比较合适如果还用24时段的参数配置聚类效果会比较差影响场景辨别的准确度。5.3 扩展方向这套代码后续可以往几个方向扩展都是当前研究的热点一是把单微网扩展成多微网互联的架构此时关键场景辨别算法需要处理的不确定源更多维度爆炸的问题会更明显但也更能体现算法的加速优势。二是引入电动汽车充放电行为的不确定性。EV的充放电行为不仅有时段的不确定性还有电量需求的不确定性建模难度更大但模型价值更高。三是在目标函数中加入碳排放约束这对应双碳背景下微电网低碳调度的需求也是现在各类基金项目喜欢的方向。6. 写在最后的调试心得这套代码我前前后后调试了两周多最大的体会是两阶段鲁棒优化的难点不在“鲁棒”而在“两阶段”。很多人把主问题、子问题分别写对就以为完事了实际上两个问题之间的变量传递、场景添加、上下界计算才是真正需要反复调试的地方。建议你自己复现时先跑一个3节点微型算例比如只有1台机组储能固定负荷把CCG的每一轮迭代结果打印出来手动验算一遍上下界的更新逻辑是否正确。确认小算例跑通后再替换成完整微网算例。别一上来就对着24时段大系统调试那样问题定位会非常痛苦。另外建议把 gap 随迭代次数的变化曲线打印出来。如果gap稳步下降说明算法逻辑正确只是收敛速度问题如果gap出现“回弹”甚至发散那一定是约束添加或场景辨识环节有逻辑错误别急着调参数先回头检查建模逻辑。最后分享一个小技巧主问题每次迭代后保存所有已识别关键场景对应的约束不要重新生成整个主问题模型而是用Yalmip的“增量添加约束”功能。这样不仅能减少建模时间还能让你随时检查某些场景的约束是否正确添加对调试非常有帮助。这套代码的正常运行只是第一步真正吃透它才能在你自己的研究里灵活改造和落地。有具体问题欢迎交流我也还在不断迭代这个框架。