刚拿到这篇《基于两阶段鲁棒优化算法的微网多电源容量配置》的论文解读和配套程序资源时我第一反应是终于有人把微电网规划的硬骨头啃下来还愿意分享源码了。搞过微电网规划的人都知道容量配置这事儿看似简单——光伏、风电、储能、柴油机各装多少容量能让投资最省、运行最稳——但一旦把风光出力的不确定性考虑进去问题性质立刻从“一个优化问题”变成“一堆优化问题”。这篇论文能发到EI、CSCD、北大核心三检索核心卖点就是用两阶段鲁棒优化Two-Stage Robust Optimization框架处理风光出力的最恶劣场景同时通过CCG列与约束生成算法把“min-max-min”这个没法直接求解的三层结构拆成主问题-子问题迭代求解。我在复现程序的过程中踩了不少坑也摸清了不少门道。这篇文章就围绕这套方法本身和程序实现细节把我的理解和实操经验完整写出来。1. 这篇论文到底解决了什么问题项目背景与核心思路拆解1.1 微网容量配置的行业痛点微电网多电源容量配置简单说就是回答一个问题在一个包含光伏、风电、储能、柴油发电机或者燃气轮机的微电网里每种电源该装多少容量才能在满足负荷需求的前提下让全生命周期的总成本最低。这问题的难度不在“算一个最优解”而在“不确定性”上。光伏出力看天吃饭风电出力看风吃饭负荷曲线本身也有波动。传统的确定性优化方法把风光出力当成一个固定值来算比如“取全年平均出力”算出来的容量配置在实际运行中很容易出现两种情况要么容量装少了遇到连续阴天或无风天微网供不上电导致切负荷要么容量装多了大部分时间设备闲置投资浪费严重。而传统的随机规划方法需要假设风光出力服从某种概率分布然后用蒙特卡洛抽样生成大量场景。但问题在于真实的天气过程很难用精确的概率分布去刻画而且抽样场景再多也可能漏掉那种“十年一遇”的极端连续恶劣天气。一旦这种没被抽样到的恶劣场景发生配置方案就可能直接崩掉。这就是鲁棒优化的用武之地。它的核心思想是我不预测你具体什么时候出多少力我直接假设一个“出力最恶劣的取值范围”只要能扛住这个范围内最恶劣的组合那实际运行一定稳。代价是配置结果会偏保守但换来了极高的可靠性。1.2 为什么是“两阶段”鲁棒而不是单阶段单阶段鲁棒优化Single-Stage Robust Optimization的思路是在所有不确定参数最恶劣的情况下同时做出“建多大容量”和“怎么运行”的决策。相当于你在买设备之前就必须把所有可能发生的运行情况全部考虑进去而且所有决策都是“事前”拍板的。这样做的麻烦在于现实中的决策是有先后顺序的。你建设光伏电站、买储能电池是“今天”就要拍板的事属于长期投资决策。但微网建成以后具体某台柴油机今天发多少电、储能今天充还是放那是“明天”根据实际天气和负荷情况才能定的属于短期运行决策。这两层决策的性质完全不同。两阶段鲁棒优化精准对应了这种决策结构。第一阶段Here-and-Now做容量投资决策——装多少光伏、多少风电、多少储能、多少柴油机这个决策一旦做了就改不了。第二阶段Wait-and-See做运行调度决策——在第一阶段决定的装机容量下针对每一种可能出现的风光出力场景找到最经济的运行方式。而鲁棒优化的“恶意”之处在于第二阶段的外部环境风光出力不会乖乖按你的预测来它会“故意”往最让你的运行成本变高的方向去变。所以整个问题变成了一个三层结构外层是投资决策min中层是找一个让运行成本最大的恶劣出力场景max内层是在这个恶劣场景下最小化运行成本min。即经典的“min-max-min”形式。这个结构不是随便选的它精准刻画了“先投资决策再面对恶劣环境做最优应对”的现实逻辑同时也让模型保有了数学上的可解性。1.3 程序包里到底有什么资源概览我拿到的这套配套程序核心就是两阶段鲁棒优化的求解实现。程序主体用MATLAB编写求解器采用Yalmip工具箱Gurobi或CPLEX的组合。这种搭配在电力系统优化领域非常主流Yalmip负责把数学模型“翻译”成求解器能识别的格式Gurobi/CPLEX负责真正把优化问题算出来。程序包的结构大概分为几个模块主程序负责搭建两阶段鲁棒模型、调用CCG迭代框架、不确定性参数生成模块负责生成风光出力的不确定集合、以及数据文件包含典型日负荷曲线、风光出力预测值、各类电源的经济性参数等。这套资源最值钱的地方在于它完整走通了“论文公式→可运行代码”的链路。论文里那些密密麻麻的约束条件每一步在代码里都有对应实现CCG算法的迭代主循环、子问题的对偶转换全部是能直接跑起来看效果的。这是很多论文做不到的——不少论文公式写得漂亮但代码复现起来各种细节缺失。2. 核心模型架构目标函数、约束与不确定性建模的精髓2.1 目标函数怎么构建投资成本与运行成本的双层博弈两阶段鲁棒优化模型的目标函数是第一阶段投资成本和第二阶段运行成本的总和。第一阶段的投资成本包括光伏、风电、储能按容量和单位投资成本折算、柴油机的建设成本这里面储能比较特殊除了功率成本还有容量成本电池本身的价格按kWh算。第二阶段运行成本主要包括柴油发电机的燃料成本和运维成本、储能的充放电损耗成本以及与主网交互的购售电成本。如果微网是并网型的那么向主网购电要花钱向主网售电有收入如果系统允许切负荷即无法满足全部负荷还要加上切负荷惩罚成本。目标函数的表达式可以写成min_x { C_invest(x) max_u min_y C_operate(y, u, x) }其中x代表第一阶段的容量配置决策u代表第二阶段的不确定参数风光出力、负荷y代表运行调度决策各机组出力、储能充放电功率、与主网交互功率。一个关键细节是投资成本前通常需要乘一个“资金回收系数”来把一次性投资分摊到规划年限的每一年或者把运行成本按年份折现到投资时点。论文里采用的是等年值法就是把总投资额乘以一个系数折算成每年等额支付的金额然后和年运行成本相加。这个系数CRF的计算公式为CRF r * (1r)^T / ((1r)^T - 1)其中r为贴现率T为规划年限。程序里默认T取20年贴现率r取0.088%这个参数可以根据实际项目调整。这个小细节非常影响最终结果——贴现率取5%和取8%光伏容量的最优配置可能差出20%以上。2.2 两阶段模型的数学表达式与约束拆解第一阶段模型的约束比较简单直接——各类电源的配置容量不能超过上限比如场地限制、电网接纳能力限制。写成标准形式就是min a^T x Q(x) s.t. x_min x x_max各电源容量上下限 x为整数或连续变量其中Q(x)是第二阶段价值的函数表示在给定第一阶段决策x的情况下最恶劣场景下的最小运行成本。这个Q(x)没法直接求出来只能通过CCG算法迭代逼近。第二阶段的运行约束就复杂多了包括功率平衡约束、各机组出力上下限约束、爬坡约束、储能充放电约束、储能SOC状态约束以及与主网交互功率约束等。用标准形式表达为对任意u ∈ U: min b^T y s.t. D y d - E x - F u功率平衡与运行约束 y 0这里D、E、F是系数矩阵体现了运行约束中决策变量y、第一阶段变量x、不确定参数u之间的耦合关系。功率平衡约束是最核心的任意时刻光伏出力风电出力柴油机出力储能放电功率购电功率 负荷功率储能充电功率售电功率。这是在程序里最容易出错的地方——因为储能既可以充电相当于负荷也可以放电相当于电源它的变量方向会改变。程序里通常用“储能SOC”作为状态变量把充放电功率统一处理成SOC(k1) SOC(k) η_c * P_ch(k) * Δt - P_dis(k) * Δt / η_d这样处理能避免同一时刻既充电又放电的逻辑矛盾。2.3 不确定性集合鲁棒优化的灵魂所在鲁棒优化相对随机规划的核心差别在于不确定性集合U的建模方式。常见的类型有三种不确定性集类型核心思想保守度计算复杂度盒式不确定集认为出力在[预测值-偏差, 预测值偏差]区间内任意取值最高最低椭球不确定集用椭球刻画不确定参数的联合分布范围中等较高需要二阶锥约束预算不确定集Budget限制所有时段中偏移预测值的时段总数上限可控较低论文采用的是预算不确定集也叫Γ-鲁棒。它的思想是虽然每个时段的光伏/风电出力都可能偏离预测值但不能指望它每个时段都同时达到最大偏差——现实里不会有这么坏的天气。引入一个“不确定预算”Γ限制所有时段里允许出现最大偏差的总时段数上限。数学上预算不确定集可以写成U { u | u_t u_t^f δ_t * Δu_t, Σ|δ_t| Γ, 0 |δ_t| 1 }其中u_t^f是预测出力Δu_t是最大偏差幅度δ_t是偏差比例变量Γ是预算参数。Γ0时退化为确定性模型Γ24时表示24个小时都可以达到最大偏差即盒式集合保守度拉满。我从程序调试中得到的经验是Γ的取值不是越大越好。Γ取值增大系统配置的容量会趋向保守容量变大投资成本上升但Γ太小又会低估风险。论文里通常取Γ为预测点数的50%~75%作为一个合理区间。具体工程上还需要根据项目所在地的历史天气数据做敏感性分析来最终标定。2.4 为什么选择CCG算法而不是Benders分解两阶段鲁棒优化的经典求解方法有两种Benders分解和列与约束生成CCG。程序选择CCG是有道理的。Benders分解的思路是把第二阶段问题对偶后把“最恶劣场景”和“运行成本”的信息通过割平面Cut一层层反馈给主问题。它的特点是需要不断往主问题里加割但加了割的约束只是原来的线性对偶关系不需要新增变量。CCG的思路不同它在迭代过程中每轮找到一个最恶劣的场景后不是仅仅加一条割而是直接把x, y联合变量和这个场景下的运行约束整块塞进主问题。也就是说每迭代一轮主问题就长胖一圈——多一组运行决策变量和多一组约束。从收敛速度来看CCG比Benders快得多。原因是Benders用对偶信息“间接”逼近最恶劣场景而CCG直接生成对应最恶劣场景的变量和约束相当于把第二阶段的实际运行情况“原样”建模进了主问题。大量数值实验表明CCG在有限次迭代内就能达到很高的精度特别适合微电网这种规模适中、但约束关系复杂的规划问题。3. 程序复现全流程从数学模型到可执行代码3.1 环境配置与求解器选择建议的环境配置是MATLAB R2020a及以上版本Yalmip工具箱用R20210331版本求解器用Gurobi 9.x如果你有学术许可或CPLEX 12.x。Gurobi和CPLEX对线性规划和混合整数线性规划的支持都极其稳定这两个随便选一个都可以。需要注意的一个坑是Yalmip版本过旧对新版MATLAB的兼容性会很差而Gurobi版本过新Yalmip可能无法识别。我踩过的组合是MATLAB R2022bYalmip旧版Gurobi 10.x直接报“No suitable solver found”——Yalmip不认识Gurobi 10的接口。解决办法是升级Yalmip到最新版或者干脆把Gurobi降到9.5版本。安装好之后在MATLAB里运行yalmiptest命令可以测试求解器是否配置成功。看到Gurobi/CPLEX后面显示“Found”状态就说明环境没问题了。这一关不过后面的所有步骤都白搭。3.2 主程序架构CCG迭代怎么实现整个程序实现的CCG算法主逻辑如下这是伪代码理解流程用初始化LB -infUB inf迭代次数k 1 设初始恶劣场景 u^(0) 为预测场景即δ_t 0 循环迭代 1. 求解主问题MP min a^T x θ s.t. x的容量约束 θ b^T y^(k)对每个已生成的场景k D y^(k) d - E x - F u^(k)对每个已生成的场景k 得到最优解(x^(k), θ^(k))更新LB max{LB, a^T x^(k) θ^(k)} 2. 将x^(k)代入第二阶段子问题SP 求解 max_{u ∈ U} min_{y} { b^T y | D y d - E x^(k) - F u } 得到最恶劣场景u^(k1)和最优目标值f(x^(k)) 更新UB min{UB, a^T x^(k) f(x^(k))} 3. 如果 (UB - LB) / UB εε取0.01或0.001则收敛退出循环 4. 否则把新场景u^(k1)对应的运行变量y^(k1)和约束加入主问题k k1回到步骤1这里最核心也最费力的一步是步骤2——子问题是两层结构max-min无法直接求解需要对内层的min进行对偶转换把“max-min”变成“max-max”也就是一个单层的max问题。对偶转换处理得当子问题就变成了一个普通的最大化线性规划问题可以直接用Gurobi求解。具体来说第二阶段子问题的对偶形式可以写成max_{u, λ} λ^T (d - E x^(k) - F u) s.t. D^T λ b λ 0 u ∈ U其中λ是第二阶段原问题的对偶变量。注意这个对偶问题里λ和u都是变量目标里有一个λ^T F u的双线性项这导致它不是线性规划。论文中通常假设不确定性集合U是盒式/预算式也就是u_t是线性约束此时利用u_t的取值范围上下界约束可进一步消去双线性项把它转化成一个可解的线性规划或混合整数规划MILP问题。如果子问题中存在二进制变量比如储能充放电状态对偶的处理会更加复杂要么用大M法把互补约束线性化要么用KKT条件替换后做线性化处理。这套程序里针对储能模型做的是连续变量处理也就是忽略充放电状态的整数特性但用SOC约束和充放电功率的上下限天然限制了储能不能同时充放。如果非要精确处理储能的二进制变量则需要引入大M法程序规模会膨胀不少。3.3 不确定性建模的代码实现细节程序里对不确定性的建模重点在两个模块预测场景的生成以及不确定集合的参数设置。预测场景数据通常由典型日曲线给定。程序里默认给出了夏季、冬季两个典型日每个典型日24个小时包含光伏出力标幺值序列、风电出力标幺值序列和负荷标幺值序列。光伏出力的预测偏差系数默认设为±20%阴雨天偏差更大晴朗天偏差小风电出力的预测偏差系数设为±15%风的大小波动剧烈负荷的预测偏差设为±5%负荷相对稳定变化有规律。这个偏差系数就是Δu_t的数值基础。不确定集合构建的核心代码如下简化描述逻辑% 不确定参数的标称值与偏差 u_nominal PV_forecast; % 光伏预测出力 u_dev 0.2 * PV_rated; % 最大偏差20%的额定容量 Gamma_PV 12; % 光伏不确定预算默认12个时段可以同时达到最大偏差 % 用Yalmip构造不确定变量 u_PV sdpvar(24, 1); for t 1:24 F [F, u_nominal(t) - u_dev(t) u_PV(t) u_nominal(t) u_dev(t)]; end F [F, sum(abs(u_PV - u_nominal) ./ u_dev) Gamma_PV];这段代码的意思是光伏每个时段的出力可以在预测值±20%内浮动但所有时段浮动总量加起来不能超过Γ对应的总预算。注意这里的Σ|δ_t|是一个绝对值约束实际建模通常引入辅助变量将绝对值线性化Yalmip里可以直接写norm或者写两对不等式来完成。我实际调试时发现Γ的取值对迭代次数影响很大。Γ12左右时CCG通常在5~8轮内收敛但如果Γ设到24盒式集合子问题每次都会找到“所有时段全部最大偏差”的极端场景主问题反复加重约束迭代次数会明显增加而且最终配置的容量会显著偏大——这是鲁棒优化“保守性”的代价。论文里用了不同Γ下的对比实验来展示这种权衡程序里也保留了参数调节的接口方便做敏感性分析。3.4 储能模型的建模要点储能是微网容量配置里最容易建模出错的环节。程序采用的方式是简化储能模型SOC作为连续变量充放电功率均为连续变量且通过以下约束限制在同一时刻只能充电或放电SOC(t1) SOC(t) η_ch * P_ch(t) * Δt - P_dis(t) * Δt / η_dis 0 P_ch(t) P_ch_max 0 P_dis(t) P_dis_max SOC_min SOC(t) SOC_max SOC(1) SOC(25)调度周期始末SOC一致注意到这里没有用二进制变量来强制“P_ch × P_dis 0”而是依靠充电效率η_ch和放电效率η_dis低于1这一损耗特性让优化问题自动倾向于不同时充放因为同时充放会产生能量损耗增加成本。这是个很聪明的建模技巧避免了二进制变量带来的计算复杂度爆炸代价是最优性会有极微小的影响。储能容量配置的另一个关键参数是“储能时长”即储能电池容量kWh对应多少小时的额定功率。程序里默认取2小时也就是一个100kW/200kWh的储能系统。这个参数直接影响储能的单位投资成本功率成本kW和容量成本kWh是分开算的程序里用两个变量分别表达。3.5 Yalmip代码骨架展示Yalmip的建模方式非常直观几乎和数学公式一一对应。我把第一阶段的容量投资问题的核心代码逻辑写出来简化版完整程序以资源包为准% 决策变量 x_PV sdpvar(1, 1); % 光伏配置容量kW x_WT sdpvar(1, 1); % 风电配置容量kW x_ES sdpvar(1, 1); % 储能功率配置kW x_ES_cap sdpvar(1, 1); % 储能容量配置kWh x_DE sdpvar(1, 1); % 柴油机配置容量kW % 目标函数投资成本等年值 第二阶段期望 CRF 0.08 * (10.08)^20 / ((10.08)^20 - 1); InvestCost CRF * (c_PV * x_PV c_WT * x_WT ... c_ES_p * x_ES c_ES_c * x_ES_cap c_DE * x_DE); % 约束容量上下限 Constraints [x_PV_min x_PV x_PV_max, x_WT_min x_WT x_WT_max, ...]; % 两阶段鲁棒目标 Objective InvestCost theta; % theta是第二阶段成本的替代变量在CCG迭代中动态更新 % 求解第一步先求解主问题 optimize(Constraints, Objective, sdpsettings(solver, gurobi));这段代码的逻辑就是先确定第一阶段的容量决策x和第二阶段成本代理变量θ然后在CCG迭代中不断往Constraints里添加新场景对应的运行变量和约束同时更新θ的取值约束直到上下界收敛。我建议拿到程序后先在原始参数下跑一遍确认结果能复现论文里的Case 1数据然后再手动修改Γ参数、贴现率和单位投资成本观察结果变化趋势。这套验证流程能帮你确认自己对模型的理解是否准确。4. 程序调试心得与微网规划实操经验4.1 影响CCG收敛速度的几个因素实际跑程序时最关心的就是迭代多少次能收敛。我实验发现除了Γ参数外以下几个因素也显著影响收敛表现松弛变量的引入如果在第二阶段加入了松弛变量比如切负荷量子问题会更容易找到“让切负荷最严重”的场景从而快速提升下界加速收敛。但如果松弛变量的惩罚系数设得太小子问题会优先选择切负荷导致配置结果不真实设得太大又可能引发数值病态问题。我建议惩罚系数至少设置为最高运行成本的10倍以上。求解器数值容差Gurobi默认的MIP gap是1e-4但对于CCG这种嵌套求解结构内层求解的微小误差会通过主问题累积放大。如果发现上下界长期无法收敛或出现LBUB这种违背理论的情况把求解器的MIPGap适当地放宽到1e-3或者把FeasibilityTol放宽到1e-6往往能解决问题。初始场景的选择CCG算法需要一个初始的恶劣场景。程序默认用预测场景即不确定参数全取标称值作为第一轮的主问题输入。事实上用极端场景光伏出力全零、负荷最大作为初始场景可以明显减少迭代次数——因为主问题在一开始就能生成一个相对接近最优的容量配置。这点小技巧值得一试。4.2 常见报错与排查速查表整理这份速查表都是我实际调试中遇到过的问题相信很多复现程序的同行也遇到过类似情况报错现象原因分析解决方法No suitable solver foundYalmip未识别Gurobi/CPLEX运行yalmiptest检查升级Yalmip或重装求解器检查求解器路径是否在MATLAB搜索路径中Infeasible problem主问题无解通常是容量上限过小或功率平衡约束过紧检查储能初始SOC约束是否导致无可行解检查风电穿透功率约束是否过严迭代次数过多几十轮不收敛子问题对偶转换出错检查对偶变量维度是否正确检查拉格朗日乘子和原问题的约束是否一一对应用一个小规模随机算例人工验证对偶间隙为零Out of memory主问题规模爆炸CCG每轮都新增变量和约束迭代二三十轮后问题规模剧增可在循环中定期“清理”冗余约束或使用optimize后的recover功能复用KKT信息结果中的储能容量为0储能经济性参数不匹配储能充放电效率太低或储能单位成本过高导致配置不上对比文献常用参数确认数量级检查贴现率是否过大导致储能收益被过度折现风光容量优于直觉偏大/偏小不确定性集合参数设置不合理检查Γ取值检查Δu_t的百分比是否过高检查负荷曲线的峰值和能量需求水平4.3 从论文到工程这套程序能怎么扩展论文和程序自带的是标准测试算例但真实工程项目比这复杂得多。我自己做过的扩展方向包括以下几个第一把单目标扩展为多目标。容量配置往往是成本与可靠性的博弈可以在目标函数中同时考虑年化总成本和供电可靠性如电量不足期望值EENS然后用ε-约束法或加权法生成Pareto前沿。程序中的第二阶段子问题已有切负荷惩罚项稍微改造惩罚系数就能把EENS作为第二优化目标输出。第二把独立微网扩展为并网微网。独立微网没有主网支撑必须满足严格的功率平衡和备用容量约束并网微网可以买卖电不确定性场景下甚至可以依靠主网购电兜底。程序里的并网模型核心只需要加两个变量购电功率、售电功率和一条联络线功率约束改造量不大但实用价值很高。第三考虑季节性差异和长期退化。光伏组件会老化、储能电池会衰减、柴油机效率会下降这些影响在20年规划期内不容忽视。可以把一个规划年细分为多个典型季分别建立不确定性集合并在储能容量约束中增加年衰减系数。这会显著增加计算量但工程结果更真实。第四接入实测数据。把典型日曲线替换成项目地实测的风光资源与负荷数据再用统计方法标定不确定集合参数。实测数据往往比预测曲线更“不光滑”对CCG子问题的求解会带来新的数值挑战但这是让程序从论文走向和项目落地的必由之路。4.4 审稿人视角的加分项如果你准备沿着这个方向发论文有几个细节能明显提升审稿人的观感。一是画出CCG的收敛曲线上下界随迭代次数的变化并解释为什么在第几轮之后收敛加速——这比单纯给最终结果更有说服力。二是做不确定性预算参数的敏感性分析说明不同鲁棒水平下容量配置的结构性变化风向、光伏、储能、柴油机的装机比例如何随Γ改变。三是增加算法对比用Benders分解、采样平均近似方法与CCG做对比给出收敛速度与解质量的定量结论。结语一些务实建议整套程序跑通、结果复现以后我最大的感受是两阶段鲁棒优化并不神秘它本质上是把“决策时间轴”和“不确定性对抗性”写进了数学模型里。主问题-子问题迭代求解的CCG框架就像一个经验丰富的谈判专家——主问题先报一个投资方案子问题立刻抬出最恶劣场景反驳主问题被反驳后修订方案、增加底线约束直到双方达成一致。这个过程看似繁琐但每一步都有明确的经济学含义。如果你刚拿到这套程序我的建议是先别急着改代码把第2章的数学模型用纸笔画一遍搞懂每个变量和每个约束在说什么再对照代码一行行看。第一次跑通之后故意改错几个参数比如把Γ设成0、把储能效率设成1看看结果如何反常——这比任何教程都能加深你对模型的理解。运行资源里给出的算例数据、输出结果和主程序、子程序都整合在一起了。如果你在复现过程中卡在某个报错上或者对某个约束的建模有疑问欢迎交流。代码能跑通只是起点真正把它用到你自己的微网规划项目里你会发现问题背后的工程判断力远比算法本身更有价值。