微网多电源容量配置的两阶段鲁棒优化Matlab实现

这两年做微网规划相关项目的人越来越多“微网多电源容量配置”和“两阶段鲁棒优化”几乎成了这类题目里的标配关键词。我自己在某个海岛微网预可研项目和另一个园区微网示范项目中都用 Matlab 完整实现过这套算法。标题拆开看其实就三层微网、多电源、容量配置。微网意味着系统规模不大但结构完整多电源意味着光伏、风电、柴油机、储能甚至燃料电池都要放进优化里容量配置则说明我们最终要回答的是“每种电源装多少”的问题。而“两阶段鲁棒优化”是整篇文章的核心解法它专门用来对付光伏和风电出力不确定带来的风险。这篇文章就把我从建模到代码、从踩坑到调参的完整过程写出来给正在做类似课题或者毕业设计的读者一个可以直接复现的参考。先说清楚这套东西能解决什么实际问题。传统容量配置往往只考虑典型日的确定性出力曲线优化出来的方案一旦遇到连续阴天、风速骤降、负荷尖峰叠加就可能出现大面积切负荷甚至系统失稳。而两阶段鲁棒优化把“投资决策”和“运行调度”分开处理允许光伏、风电在给定的不确定集合内取最坏场景然后在这个最坏场景下仍然保证系统能安全经济运行。也就是说最后得到的容量配置不是“平均意义上最优”而是“最坏情况下也能扛住”。这个特性在微网这种孤岛运行场景下尤其重要因为孤岛微网没有大电网兜底能源安全是第一位的。适合参考这篇文章的读者包括做微网规划方向的研究生、正在做新能源容量配置项目的工程师以及想快速上手 CCG列与约束生成算法的 Matlab 使用者。1. 微网容量配置问题到底卡在哪里不确定性与投资决策的博弈1.1 为什么光算一个典型日远远不够很多入门级的容量配置程序用的是“典型日法”把一整年的光伏、风电、负荷曲线聚合成几个典型日然后建一个确定性优化模型目标函数是投资成本加运行成本最小。这个方法思路简单算得快但有个致命的隐含假设——典型日里用到的出力曲线就是实际运行时的出力曲线。现实根本不是这样。光伏出力受云层遮挡影响五分钟内就能从 80% 额定出力掉到 30%风电更是出了名的“看天吃饭”某些地区夏天连续一周无风也是正常的。如果你只在确定性模型里按典型日照着一条曲线去优化那设备容量大概率是贴着这条曲线的峰值去配的没有冗余。真到了实际运行那天光照不足、风速偏低、负荷又偏高几个因素一叠加储能电量放空之后就只能靠柴油机硬扛扛不住就得切负荷。我在某个海岛项目里做过对比用确定性模型配出来的方案在连续三天阴雨、风速只有预测值 60% 的场景下缺电率达到 7% 以上。这个缺电率对于海岛微网是绝对不可接受的。所以要解决的核心矛盾是投资决策现在就要做但未来的风光出力不确定。容量配置本质上是一个“现在拍板、未来买单”的问题你必须在决策时就考虑到最坏情况而不是等到运行再去补救。1.2 两阶段鲁棒优化如何回答“最坏情况下怎么办”两阶段鲁棒优化的思路一句话概括就是先做投资决策再做最坏场景下的运行调度最后看总成本能不能接受。数学上写作min_x [ C_inv(x) max_u min_y C_op(x, y, u) ]这个式子看起来复杂拆开就清楚了。x 是第一阶段变量代表各电源的装机容量、储能容量u 是不确定性变量也就是光伏、风电的实际出力或者负荷波动y 是第二阶段变量代表每个时段各台设备的出力、储能的充放电功率、购售电量。C_inv 是投资成本C_op 是运行成本。里层 max min 的意思是在所有可能出现的“坏场景”里挑一个让调度成本最高的而外层 min 的意思是我们选择的容量方案要能让这个最坏情况下的成本尽量低。注意这里不是把每个场景都遍历一遍而是直接去寻找那个“最能让你头疼”的场景然后保证在这个场景下系统还能正常运行。这正是鲁棒优化和随机规划最大的区别随机规划讲究概率分布和期望值鲁棒优化讲究边界和最坏情况。微网容量配置这种一次投资、多年运行的问题用鲁棒优化来兜底是完全合理的思路。2. 数学模型这样搭投资决策与运行调度分离2.1 第一阶段目标钱怎么花最值第一阶段只做一件事决定每种电源装多少。这里需要定义决策变量比如光伏装机容量、风电装机容量、柴油机额定功率、储能额定容量和额定功率。对应的成本包括单位容量的购置成本、安装成本、土地或屋顶占用成本等。这一阶段的目标函数就是总投资的等年值折算通常用资金回收系数把一次性投资分摊到每一年便于和后面的年运行成本放在同一个量纲里比较。模型里还应该考虑设备上下限约束比如光伏安装面积限制、柴油机台数限制、储能最大容量限制。这些约束在代码里实现起来很简单就是一堆上下界但它们的意义很大没有这些约束优化器一定会把所有能装的都装到上限去占满屋顶结果完全失真。2.2 第二阶段目标在最坏场景下怎么调度最省第二阶段是在第一阶段容量方案 x 固定的前提下针对某个具体的风光出力场景 u求解一个经济调度问题。调度变量包括每个时段柴油机的出力、光伏和风电的实际出力在不超过可用出力的前提下可以弃光弃风、储能的充放电功率、储能荷电状态、以及与外部电网的交换功率。目标函数是运行成本最小主要包括柴油机燃料成本、启停成本、购电成本必要时还可以加上弃光弃风的惩罚项。同时要满足的约束比较多但每个约束都有物理含义。功率平衡约束是核心它要求每个时段所有电源出力加储能放电功率加上购电功率等于负荷加储能充电功率加上售电功率。储能约束包括充放电功率限值、容量上下限、SOC 动态递推关系以及调度周期始末 SOC 相等的约束。柴油机还有爬坡约束和最小技术出力约束。这些约束不多但全部要写进子问题里因为子问题对偶转换后这些约束会直接影响对偶变量的个数和边界条件。2.3 不确定性集合与预算参数怎么设u 不能随便取必须限定在一个“不确定性集合”里。最常用的是盒式集合加预算约束。盒式集合的意思是每个时段的光伏出力在预测值正负偏差范围内波动风电也一样。预算约束的意思是所有时段的波动量加起来不能超过一个总预算 GammaGamma 越大允许的“集体偏离预测值”程度越高优化结果越保守。这里的 Gamma 就像保险的保额Gamma 等于 0 时等同于确定性优化Gamma 取最大值时就是所有时段全部取最坏值相当于“最悲观专家”替你做决策。实际工程里 Gamma 一般取调度周期时段数的三分之一到二分之一既可以规避大部分极端情况又不至于让投资成本膨胀到没法接受。我后面会专门讲 Gamma 取值对结果的影响。3. CCG 算法原理与 Matlab 实现3.1 从 min-max-min 到主问题-子问题迭代min-max-min 这种三层结构没法直接交给求解器必须通过分解算法把它拆成可以迭代求解的主问题和子问题。主问题 MP 是在一个离散的场景集合里做投资和运行联合优化场景集合一开始只放一个预测场景之后每次迭代加入一个新的最坏场景。子问题 SP 则是在当前投资方案 x 固定的情况下求解内层的 max-min 问题找到那个最坏场景 u_star并返回该场景下的运行成本。迭代过程很简单但每一步都考验细节。先初始化场景集合为预测场景设置下界 LB 为负无穷、上界 UB 为正无穷。每次迭代先求解主问题 MP得到当前投资方案 x 和下界 LB。然后固定 x求解子问题 SP得到最坏场景 u 和对应的运行成本。把投资成本加上这个运行成本得到一个新的总成本用它更新上界 UB。如果上下界相对误差小于收敛阈值就输出结果否则把 u 对应的场景作为一个新的离散场景加入主问题继续迭代。一个很多人忽略的细节是主问题因为是“场景集合”内的松弛问题所以它的最优值给出了原问题的下界而子问题给出了一个可行方案的总成本所以是上界。这个上下界的关系是 CCG 算法收敛性分析的基础也是你在代码里判断“算没算对”的第一道防线。如果代码写反了上下界越迭代相差越大基本就是对偶方向或者更新逻辑出了问题。3.2 子问题对偶转换的关键一步子问题最麻烦的地方在于它是个 max-min 结构求解器不能直接处理。标准的处理方式是把内层 min 问题对偶成一个 max 问题于是 max-min 就变成了一个纯 max 问题。注意这里对偶的方向原问题是最小化运行成本约束是标准形式对偶之后目标函数里会出现不确定性变量 u 与对偶变量的乘积产生双线性项。举个例子功率平衡约束的右侧往往会包含光伏和风电的可用出力也就是 u 的一部分。对偶之后这些 u 会乘上功率平衡约束对应的对偶变量形成一个 u 乘以 λ 的乘积项。这种双线性项是非凸的没办法直接交给 CPLEX 或 Gurobi 求解。工程里最常用的处理办法是 Big-M 线性化把不确定性变量 u 离散到有限个取值点或者引入辅助二进制变量把连续乘积转化为一系列线性约束。工作量大了点但只要 M 选得足够大效果非常稳。我在代码里用的是分段线性化加二进制变量展开的方法。先把每个时段的光伏、风电出力偏差离散成几个档位然后对每个档位引入二进制变量用大 M 把 u 与 λ 的乘积项拆开。实测这种做法的求解时间基本可控而且能够保证子问题求出来的一定是全局最优的最坏场景。3.3 MatlabYalmip 代码骨架与核心逻辑直接说代码结构。整个程序我分成了四个文件数据初始化文件、主问题建模函数、子问题建模函数、主循环文件。数据初始化负责读入负荷曲线、预测出力曲线、设备参数、成本参数主问题函数接收场景集合和当前所有场景下的运行变量定义返回优化模型子问题函数接收投资方案 x返回最坏场景和运行成本主循环负责迭代逻辑和结果汇总。下面是一段核心伪代码真实程序在此基础上扩展场景约束和变量定义即可。% 初始化 u_set {u_forecast}; % 场景集合初始是预测场景 LB -inf; UB inf; tol 1e-3; iter 0; while (UB - LB) / abs(UB) tol iter 20 % 求解主问题 MP build_mp(u_set); % 返回主问题的目标、约束、变量 optimize(MP.cons, MP.obj, ops); LB value(MP.obj); x_val value(MP.x); % 当前最优投资方案 % 求解子问题 SP build_sp(x_val); optimize(SP.cons, SP.obj, ops); SP_cost value(SP.obj); u_star value(SP.u); % 更新上界 UB min(UB, value(invest_cost(x_val)) SP_cost); % 加入新场景 u_set{end1} u_star; iter iter 1; end主问题里最需要注意的是场景变量的组织方式。我习惯把所有场景下的运行变量存成一个 cell 数组每个 cell 对应一个场景的完整变量集合。添加新场景的时候只需要新增一个 cell并把对应的约束也追加进去。这样做的好处是代码结构清晰调试的时候能逐个场景检查变量的值坏处是场景数量多了之后内存占用上升得很快后面我还会说这个问题的解决办法。子问题里最需要注意的是固定 x 的方式。用 Yalmip 写的话不能简单地把 x 替换成数值因为模型中存在 x_val 相关的约束表达式。我的做法是把 x 声明成 sdpvar 变量然后在求解子问题之前用 assign 函数把 x 的值赋进去再调用 optimize。这样既保持了模型的符号表达能力又能保证 x 在子问题里确实是常数。4. 算例结果分析与经济性解读4.1 配置结果对比确定性方案与鲁棒方案的区别我用一个 24 节点时段的测试算例跑了完整迭代。系统包含光伏、风电、柴油机和储能负荷峰值大约 1.2MW。先看确定性优化结果Gamma 取 0 的时候光伏装 500kW、风电装 200kW、柴油机装 1000kW、储能容量 1200kWh总年化成本 1800 万元左右这里面投资折算占了大概 40%。当 Gamma 取到 6也就是允许总共 6 个时段的出力偏差叠加到最坏情况时配置结果发生了明显变化。光伏降到 350kW风电降到 150kW柴油机升到 1200kW储能容量升到 1500kWh总年化成本涨到 2200 万元左右。这个结果非常符合直觉风光出力越不确定优化器就越不愿意把宝押在风光上转而用柴油机和储能来兜底。光伏和风电虽然边际成本接近零但它们的出力不可控在最坏场景下反而成了“靠不住”的电源所以容量被削减。这里要给读者提个醒不要看到鲁棒方案总成本更高就觉得鲁棒优化“不划算”。鲁棒方案多出来的 400 万成本本质上是为“最坏场景下不切负荷”这个保险支付的保费。在孤岛微网里一次非计划停电造成的损失往往远超这 400 万所以这笔账要结合可靠性需求来算不能只看投资数字。4.2 预算参数 Gamma 对容量和总成本的影响Gamma 是最值得做敏感性分析的参数。我把 Gamma 从 0 逐步增加到 12每档都重新跑一遍完整的两阶段迭代。结果是典型的边际递减曲线Gamma 从 0 到 4 那一档总成本上升最快涨幅接近 15%Gamma 从 8 到 12 那一档总成本涨幅不到 3%。这说明不确定性的“边际杀伤力”是在递减的最坏场景里的前几个波动对系统的影响最大。工程上我建议做两步走。第一步先用 Gamma 的中间值比如 6 跑出参考方案第二步把 Gamma 调大调小各跑两档看容量配置结果变化大不大。如果配置结果在相邻两档之间变动很小说明方案对不确定性不敏感可以放心用如果变动很大说明系统本身比较脆弱需要进一步增加储能或者升级线路。这个做法比纠结于某一个精确的 Gamma 值实用得多。4.3 收敛性曲线怎么看主循环的收敛性是一个很容易被忽视但非常关键的检查点。我习惯把每轮迭代的下界和上界都画成两条曲线放在同一张图里。正常情况下下界从很低的值向上爬升上界从很高的值向下收敛两条曲线在 5 到 8 轮迭代之内会靠得很近相对误差低于 1%。如果上界曲线出现震荡也就是某次迭代的总成本比上一轮还高别急着怀疑算法错了先检查是不是子问题返回的最坏场景没有真正被加入主问题或者加入时约束写漏了。如果下界几乎不动、一直是同一个值而上界在缓慢下降说明主问题里新增的场景没有起到收紧可行域的作用。这种情况往往是不确定集合定义得太小极端场景和预测场景区别不大导致每次迭代加进去的新场景形同虚设。解决办法也很直接把不确定集合的偏差范围调大或者强制子问题的目标函数里加入一个最小改善量。5. 常见问题与排查实录跑不通的坑我都踩过了5.1 子问题无界或目标为负这是 CCG 新手最容易碰到的问题。子问题应该是求一个最大值如果你的子问题目标函数在求解器里是负数而且求解器还提示 unbounded那基本可以断定是对偶转换出错了。最常见的原因是把对偶的方向写反或者对偶变量的符号约束写错。我当时排查这个问题花了一整天最后是拿一个固定场景手工代入验算找到的漏洞功率平衡约束写成了小于等于的形式导致对偶变量符号反了。对偶转换的每一行都必须和原约束一一对应我建议在代码里把原问题和对偶问题的约束列表同时打印出来一行一行对着检查比自己盯着屏幕瞎猜快得多。5.2 求解时间爆炸迭代轮数不多但每一轮主问题求解时间越来越长最后甚至单轮就要跑半小时。这个问题的根源在于主问题里的二进制变量太多。储能充放电状态、柴油机启停状态、加上每个场景都要复制一套变量场景一多主问题的规模就指数膨胀。我的解决方案是把“整数变量只保留在第一阶段和柴油机启停”上储能充放电状态先用连续变量近似等主问题算完再在第二阶段精细校验。如果确实需要保留储能二进制变量那就加一个约束相邻时段储能充电状态和放电状态不能同时为 1并给求解器设置一个合理的 MIP gap比如 0.5%省下来的时间远大于多出来的那点误差。5.3 场景一多主问题内存不够跑 20 轮迭代主问题里就有 20 套完整的运行变量和约束内存失控是很正常的。一个比较实用的做法是限制总场景数量比如最多保留 10 个场景超过之后就淘汰对目标函数影响最小的旧场景。淘汰的依据是看每个场景对应的对偶变量值对偶变量越大说明这个场景对目标函数的约束作用越强应该保留反之可以删除。我试过这个办法在保证上下界误差不反弹的前提下内存占用降了差不多一半。如果你的问题场景数特别多还可以考虑用 Benders 分解的思路只保留对偶割平面而不是完整场景但这套方案实现复杂度更高我建议先把 CCG 跑通再考虑。5.4 参数设置经验总结这里把参数设置的经验整理一下方便直接对照使用。陷阱点现象处理经验不确定性偏差范围太小迭代几轮就收敛但结果和确定性优化几乎一样把偏差加大到预测值的 20%-30% 再试Big-M 值取太小子问题求解结果跳变无法稳定收敛M 值至少取系统最大功率的 10 倍收敛阈值太严迭代超过 15 轮还不满足工程场景用 1e-3 就够别追求 1e-6没有设置柴油机最小出力子问题出现不合理的停机状态加最小技术出力约束和启停变量储能 SOC 初值不当第一天调度结果明显不合理强制调度周期始末 SOC 相等6. 个人实操心得与后续扩展这套代码跑通之后我又在几个不同项目上做了调整最大的体会是两阶段鲁棒优化难的不是算法本身而是把物理模型和算法框架合在一起的那一步。你花在写主循环上的时间可能只占两成剩下的时间全在调约束、调参数、排查对偶问题。尤其是当你发现一个结果特别“完美”的时候先别高兴用极端场景去压力测试它一下往往能暴露不少问题。最后再分享一个小技巧迭代完成之后把每个场景下的柴油机出力和储能电量画成曲线和预测场景的曲线放在一起对比。如果最坏场景下柴油机长时间贴着上限跑说明柴油机容量配置偏小应该增加如果储能基本没有放空过说明储能容量偏大可以适当削减。这套可视化的验证手段比单纯看成本数字更能帮你理解配置结果的合理性。后来我也在代码里加了碳交易成本和失负荷惩罚项把单目标扩展成多目标权衡做成另一套工具原项目方案最后验收时也是用这套 CCG 的结果作对比基准。希望这篇总结能帮你少踩几个坑把精力集中在真正需要思考的地方。