做微电网调度这块的人多多少少都被同一个问题折磨过光伏和风电出力说不准。今天预测明天中午光照不错结果一片云飘过来出力直接掉三成夜间风电大发负荷却只有那么一点电送不出去又存不下只能眼睁睁看着弃风。确定性优化模型算出来的“最优”结果在实际运行里往往大打折扣。而含可再生能源和储能的区域微电网最优运行问题恰恰就是要在这种充满不确定性的环境里找到一个既经济又安全、还不至于过度保守的调度方案。这篇文章我基于一个可完全复现的多阶段鲁棒调度模型把建模、求解器配置、CCG算法迭代、MATLABYALMIPCplex代码实现以及各种坑从头到尾捋一遍。适合正在做微电网优化调度、鲁棒优化入门、或者被“不确定集”和“min-max-min”结构搞得头疼的研究生和工程师参考。1. 项目背景与核心问题拆解1.1 微电网最优运行到底在优化什么区域微电网通常包含光伏、风电、储能、负荷以及和上级电网相连的公共连接点有些还会带微型燃气轮机。所谓“最优运行”简单说就是在满足负荷供电、设备运行约束的前提下让一天的运行总成本最低。成本项包括向上级电网购电的费用、微型燃气轮机的燃料成本、储能充放电带来的损耗成本以及弃风和切负荷的惩罚成本。但这里有一个关键点这个“最优”是在什么信息条件下求出来的。如果是预测值全部已知、且完全准确的确定性优化那问题就是一个标准线性规划或混合整数线性规划用求解器一把梭就能跑完。然而实际运行中光伏和风电的预测误差是客观存在且不可消除的。你要是完全按预测值做计划到了实际时刻很可能出现两个问题一是可再生能源实际出力低于预测导致功率缺额不得不临时高价购电或者切负荷二是实际出力高于预测系统又不敢多消纳被迫弃风弃光。无论哪种都会让“最优”偏离理想值。所以研究微电网鲁棒调度的核心动机就是要找到一条调度策略无论风光实际出力在合理预测区间内如何波动这条策略都能保证系统安全运行并且运行成本的“最坏情况”尽可能低。1.2 不确定性来源与鲁棒优化的引入微电网的不确定性来源比大电网更复杂因为系统规模小、惯性小、调节手段有限。主要不确定项包括光伏出力的间歇性、风电出力的随机波动、负荷预测偏差甚至是实时电价的波动。在多阶段鲁棒调度模型里通常把这些不确定量建模成“不确定集”假设它们在一个有界集合内变化而不需要精确的概率分布。这里就有个很有意思的对比。随机优化需要用场景树或蒙特卡洛抽样来描述不确定性计算量大而且概率分布假设不准时结果会偏鲁棒优化则是“保险思维”——我不赌概率我就假设最坏情况会出现然后保证在最坏情况下方案仍然可行且成本可控。代价是结果会比随机优化保守一些但换来的是工程上特别看重的鲁棒可行性。这个项目标题里特意强调“考虑鲁棒性和不确定性”说白了就是拒绝拍脑袋加备用容量而是用一套严谨的数学框架把不确定性显式地放进优化模型让调度计划自带“抗扰动”能力。注意鲁棒优化不是越保守越好。过分保守的结果是储能永远满充待命、购电合同永远按最大负荷签成本高得离谱在工程上没有实用价值。所以模型里通常会引入不确定预算参数鲁棒度用来调节保守程度。2. 多阶段鲁棒调度模型的数学构建2.1 确定性基础模型先搞清楚要优化的对象在引入不确定性之前先把确定性模型写清楚。这个项目里的微电网系统拓扑简化后包括光伏机组、风电机组、储能系统、微型燃气轮机、本地负荷、PCC点与主网交换功率。决策变量按时间尺度展开一个典型调度周期是24小时步长1小时也可以细化到15分钟。目标函数一般写成min Σ_t [ c_buy(t)·P_buy(t) - c_sell(t)·P_sell(t) c_gas·P_mt(t) c_bess·(P_ch(t)P_dis(t)) λ_curtail·P_curtail(t) λ_load·P_shed(t) ]其中购电和售电价格是分时电价燃气轮机成本用线性近似储能成本可以理解成充放电的折旧成本弃风和切负荷则用惩罚系数表示数值设得很大逼着优化模型优先保证供电和消纳。约束条件包括功率平衡约束P_pv(t) P_wt(t) P_mt(t) P_dis(t) P_buy(t) P_load(t) P_ch(t) P_sell(t) P_curtail(t) P_shed(t)储能动态约束SOC(t1) SOC(t) η_ch·P_ch(t)·Δt - P_dis(t)·Δt/η_dis储能容量与充放电功率上下限约束联络线传输功率约束燃气轮机出力上下限和爬坡约束弃风弃光与切负荷的非负约束这些约束构成了一个典型的混合整数线性规划。为啥会有整数变量因为储能充放电状态、购售电状态、燃气轮机启停状态都是0-1变量否则模型会在同一个时刻既充电又放电、既买电又卖电产生无意义的“环流”。2.2 盒式不确定集与鲁棒对应的引入现在把不确定性放进来。光伏和风电的实际出力可以写成P_pv(t) P_pv_forecast(t) ξ_pv(t)·P_pv_error_max(t)这里的ξ_pv(t)是标准化扰动取值范围[-1,1]。类似地风电和负荷也有同样的表达。把这些扰动向量收集成一个不确定向量ξ然后定义盒式不确定集U { ξ | ξ_min ≤ ξ ≤ ξ_max, Σ_t |ξ(t)| ≤ Γ }这个Γ就是不确定预算。它限制的是“同一时间区间内最多有多少个点同时达到最坏偏差”本质上是给保守度加一个上限。Γ越大系统要对抗的极端场景越多方案越保守Γ0时模型退化成确定性模型。用盒式不确定集而不是更复杂的多面体不确定集主要是出于两个考虑第一实际工程里对风光出力的误差范围比较容易从历史数据中统计出来上下界是现成的第二盒式集在CCG算法里处理起来简单子问题的对偶变换和线性化都方便求解效率高。2.3 两阶段鲁棒优化的min-max-min结构拆解这是整个模型最核心也最容易劝退人的地方。多阶段鲁棒调度模型本质上是两阶段决策第一阶段也就是“日前决策阶段”需要在不确定性实现之前作出决策。这些决策包括燃气轮机的启停状态、与主网的购售电状态、储能的充放电状态等“0-1开关量”以及机组的基础出力计划。这些决策一旦定了短期内改不了所以叫“here-and-now”决策。第二阶段是“实时调整阶段”等到风光实际出力确定之后系统在第一阶段决策给定的前提下通过调整储能功率、机组出力、弃风切负荷量等连续变量来保证功率平衡和运行安全。这一层决策可以随场景变化而灵活调整叫“wait-and-see”决策。整个鲁棒模型可以写成min_x { c^T·x max_{ξ∈U} min_{y∈Ω(x,ξ)} d^T·y }外层的min对应第一阶段决策成本内层的max-min是核心表示在所有可能的不确定场景里找到一个让第二阶段运行成本最高的“最坏场景”并在这个最坏场景下求最小的调整成本。两层合起来就是让第一阶段方案在“最坏情况下”的总成本最小。这个结构在工程上可以这样理解你在前一天晚上制定计划时不是问“明天最可能的场景是什么”而是问“如果明天风光出差错哪种差错组合对系统最不利我能否扛得住”。第二阶段的min就是当最坏场景真的来了系统通过储能、机组、切负荷等手段做出的最优响应。需要注意这个模型不是单纯的多阶段决策而是两阶段鲁棒优化框架。所谓“多阶段”体现在储能SOC等时序耦合变量跨时段传递子问题内部本身是含时序的动态优化。3. 求解思路与CCG算法完全解析3.1 为什么不能直接用简单鲁棒方法解决这类问题最直观的思路是把不确定集所有顶点依次带入确定性模型取最坏结果对应的方案。听起来不难但组合爆炸如果24个时段每个时段有两个不确定参数风电、光伏每个参数取上下两个顶点组合数就是4的24次方根本算不完。而且即便算完了得到的结果往往是最极端组合下的方案保守得吓人成本和实际运行需求严重脱节。另一种思路是用强对偶把内层min问题转换成max问题把整个min-max-min转成单层优化。办法可行但问题规模会迅速膨胀而且遇到这种含0-1变量的外层问题直接对偶并不容易处理。所以实际工程中最常用的是Benders分解和CCG列与约束生成Column-and-Constraint Generation两类分解算法。CCG因为收敛速度快、生成割的质量高在这类问题里几乎是首选。3.2 CCG迭代主问题与子问题的配合CCG算法说起来并不神秘核心就是“先猜一个最坏场景求解带这个场景约束的主问题然后固定主问题得到的第一阶段决策去子问题里找真正的最坏场景如果找到的场景让目标值变大了就把它作为新的约束加回主问题继续迭代”。具体迭代流程如下初始化设置下界LB-∞上界UB∞设定初始不确定场景ξ_0通常取预测值或某个典型极端场景迭代次数k1收敛精度ε。求解主问题主问题包含第一阶段决策x以及对应已知场景集{ξ_1, ξ_2, ..., ξ_k}的第二阶段决策变量y_i。由于每个场景对应一套第二阶段变量主问题规模会随迭代次数增加而变大。主问题的最优值记为下界LB_k。求解子问题固定主问题求出的第一阶段决策x_k代入第二阶段模型求解f(x_k) max_{ξ∈U} min_{y∈Ω(x_k,ξ)} d^T·y得到最坏场景ξ_{k1}和对应的最坏运行成本f(x_k)。此时更新上界UB_k min(UB_{k-1}, c^T·x_k f(x_k))收敛判断如果UB_k - LB_k ≤ ε迭代终止输出当前方案否则将ξ_{k1}作为新的场景加入主问题的场景集kk1返回第2步。这个算法的巧妙之处在于主问题每添加一个场景约束其实就是在逼近真实的最坏情况。好比你去谈判对方每次提出一个新的无理要求你就把这个要求写进合同里最后合同条款越来越多但每一版都比上一版更接近对方的真实底线。3.3 子问题中max-min的处理技巧强对偶与大M法子问题是个max-min问题求解它需要先把内层min转成max从而把整个问题统一成一个单层max问题。具体做法是写出内层LP的对偶问题利用强对偶定理把min换成对偶max与外层max合并。这里有个关键难点转换后的目标函数里会出现不确定变量ξ与对偶变量π的乘积项也就是双线性项模型瞬间变成非线性的没法直接用MILP求解器处理。工程上的标准解法是大M法引入辅助变量把双线性项线性化。大M的取值不能太大也不能太小太小可能把可行域压缩掉太大又会在数值上引发病态问题。实际测试下来取约束相关量级上限的10到100倍通常比较稳。如果第二阶段模型不含整数变量那么LP对偶是完全有效且严格的。一旦把储能充放电状态等整数变量放进第二阶段问题就变成混合整数规划对偶性不再成立处理难度陡增。所以在标准CCG框架中通常会让第一阶段承担全部0-1决策第二阶段只保留连续变量。实操提示第二阶段子问题如果遇到数值不稳定或对偶无界优先检查第二阶段约束是否冗余以及不确定集预算约束写没写对。这类问题八成出在建模细节上而不是算法本身。4. MATLABYALMIPCplex完整复现流程4.1 工具箱与求解器环境配置这个项目的完整复现我推荐环境组合是MATLAB R2022a或更新版本 YALMIP R2021版以上 Cplex 12.10或Gurobi 9.x。YALMIP作为建模层屏蔽了底层求解器的差异写约束跟写数学公式差不多直观。安装时有一个老生常谈但特别容易翻车的点YALMIP的路径必须在MATLAB搜索路径中排在Cplex工具箱之前否则调用Cplex时可能被MATLAB自带的优化工具箱截胡导致求解器识别异常。装完后在命令行跑一下yalmiptest它会列出每个求解器的状态看到Cplex和Gurobi都显示可用再开始建模。4.2 系统参数与数据准备以某区域微电网为例为了方便复现我用一个典型的区域微电网算例参数全部列在下面。这个拓扑包含一台光伏机组、一台风电机组、一组储能电池、一台微型燃气轮机和主网联络线。设备参数项数值光伏额定装机容量500 kW风电额定装机容量800 kW储能额定功率500 kW储能额定容量1000 kWh储能初始SOC50%燃气轮机额定功率300 kW联络线最大传输功率600 kW负荷峰值800 kW调度周期时段数24小时光伏和风电的日前预测出力曲线按典型日实测数据归一化后缩放得到预测误差上界设为预测值的15%~20%。分时电价设置成峰、平、谷三段价格分别是1.2元/kWh、0.7元/kWh、0.35元/kWh。储能充放电效率取0.95SOC上下限取20%和90%。4.3 核心代码实现主问题与子问题的关键片段下面给出的是复现过程中的核心代码骨架我用YALMIP语法实现。首先是搭建基础模型参数% 定义时段 T 24; % 光伏/风电预测 Ppv_fc pv_forecast; % 1x24 Pwt_fc wt_forecast; % 误差上界 delta_pv 0.2 * Ppv_fc; delta_wt 0.15 * Pwt_fc; % 不确定预算 Gamma_pv 8; % 24个时段里最多8个时段达到偏差上限然后是主问题的YALMIP建模。注意CCG迭代过程中场景集合会不断扩充所以主问题的第二阶段变量需要定义为元胞数组每轮迭代添加一组% 第一阶段变量 x_mt binvar(1, T); % 燃气轮机启停 z_buy binvar(1, T); % 购电状态 z_sell binvar(1, T); % 售电状态 z_ch binvar(1, T); % 储能充电状态 z_dis binvar(1, T); % 储能放电状态 P_mt sdpvar(1, T); P_buy sdpvar(1, T); P_sell sdpvar(1, T); P_ch sdpvar(1, T); P_dis sdpvar(1, T); SOC sdpvar(1, T1);第二阶段的调整变量放元胞数组里for k 1:Kmax y{k}.P_buy2 sdpvar(1, T); y{k}.P_sell2 sdpvar(1, T); y{k}.P_ch2 sdpvar(1, T); y{k}.P_dis2 sdpvar(1, T); y{k}.P_shed sdpvar(1, T); y{k}.P_curtail sdpvar(1, T); % 针对场景 xi{k} 构建约束详见下一段 end关键的功率平衡约束针对每个已知场景都要加一套for k 1:iter Constraints [Constraints, Ppv_fc xi_pv{k}.*delta_pv Pwt_fc xi_wt{k}.*delta_wt ... P_mt y{k}.P_dis2 - y{k}.P_ch2 y{k}.P_buy2 - y{k}.P_sell2 ... P_load - y{k}.P_shed y{k}.P_curtail]; % 储能SOC递推 Constraints [Constraints, SOC(2:T1) SOC(1:T) 0.95*y{k}.P_ch2*dt - y{k}.P_dis2*dt/0.95]; end子问题的核心是处理max-min。我通常不手动写对偶而是借助YALMIP内部机制先用yalmip的min表达式建模内层再用dualize或recover手段对偶化。如果版本支持不方便就需要手工推导对偶表达式并把双线性项用大M线性化。这里给出手工处理后的子问题目标函数片段% 子问题固定 x_star 后求解最坏场景 % 变量不确定量 xi_pv, xi_wt第二阶段调整量 y % 经过对偶变换后目标函数含有 xi * lambda 项 % 用大M法引入辅助变量 z_xi_lambda等价线性化 MBig 1e6; Constraints [Constraints, z1 lambda - MBig*(1 - abs(xi_up)), ... z1 -lambda - MBig*(1 - abs(xi_up)), ... z1 lambda MBig*(1 - abs(xi_up)), ... z1 -lambda MBig*(1 - abs(xi_up))];这段代码的细节不值得完全照抄因为它依赖具体对偶推导。复现时更推荐的做法是手写对偶问题然后用YALMIP直接建模对偶问题避免依赖大M线性化。对偶问题的推导虽然繁琐但每一步都是确定的不容易出错。4.4 确定性模型与鲁棒模型结果对比分析模型跑通后的结果对比是整个复现里最有价值的部分。我以某天的数据为例确定性优化Γ0的总运行成本是约1.83万元当Γ8时鲁棒方案总成本上升到约2.21万元多出的3800元就是“鲁棒性溢价”代价换来的是系统在风光大幅偏差下不切负荷、不违反联络线功率约束。看储能SOC曲线会更直观。确定性方案里储能倾向于在低谷电价时段多充电、高峰时段多放电SOC曲线几乎是贴着边界来回跑经济性拉满但几乎没有安全裕度。鲁棒方案的SOC曲线则留出了更多余量不会轻易充到90%上限或放到20%下限因为最坏场景里需要储能随时有能力补偿功率缺额。风电和光伏场景下鲁棒方案还会主动把燃气轮机的出力抬高一些相当于花钱买“旋转备用”。这些结果都在验证一个结论鲁棒优化不是靠运气而是靠“预支”一部分经济性换取安全性。5. 复现过程中最常踩的坑与排查技巧5.1 求解器报错与模型调试复现这类模型时我最常遇到的报错就是YALMIP提示“Unable to prove infeasibility”或“NaN in objective”。前者多半是约束写重了或者某些变量没定义完整后者则常见于SOC初始值没给、功率平衡约束里符号搞反、或者价格参数中出现未初始化变量。调试技巧上我习惯把模型拆开测先跑确定性模型确保所有约束和求解器交互正常再单独验证子问题固定一组x后看子问题能否求解最后才跑完整CCG。这样出了错能快速定位是外层迭代的问题还是内层模型的问题。另外YALMIP里可以用check(Constraints)检查约束残差如果某条约束残差很大逐条打印出来看数值很快能找到是哪个约束写错了。5.2 CCG收敛慢与最优性间隙问题CCG算法理论上收敛很快但实际复现时常常遇到迭代十几轮还不收敛的情况。排查方向有三个第一不确定集定义不合理。比如 Γ 设得过大导致子问题最坏场景特别极端主问题不断添加场景约束却始终追不上最坏值。解决办法是把 Γ 从0开始逐步增加观察UB和LB间隙的变化确定一个适中值。第二子问题对偶推导有误导致生成的割太弱甚至错误。这种问题最隐蔽表现为主问题目标值在迭代过程中出现回退。我建议子问题求解后手动校验把求出的最坏场景代回第二阶段原问题非对偶形式看目标值是否一致不一致就是对偶变换或线性化出了问题。第三大M参数不当。M值太小会导致子问题可行域被错误收紧求出的“最坏场景”实际上不是最坏M值太大会让求解器数值稳定性下降。建议对M值做敏感性测试选一个让子问题结果稳定不变的最小数量级。5.3 不确定预算Γ调节的工程建议既然 Γ 直接决定方案的保守程度那它怎么取才合理我的经验是如果历史预测误差数据充足可以统计每天实测误差超过预测上限的时段数分布取90%分位数作为Γ的参考值。如果数据不足采用“试凑法”以0.5为步长从0往上递增每档跑一次CCG记录成本和最坏场景下的切负荷量。成本曲线会出现一个明显拐点拐点对应的Γ通常就是兼顾经济性与鲁棒性的选择。负荷重要程度高的园区微电网Γ适当取大普通商业微电网Γ可以取小一些没必要为小概率极端场景付出过高的成本代价。还有一个容易被忽略的操作细节光伏和风电的不确定预算应该分开设置。光照和风速的随机特性差异很大混在一起设置Γ会掩盖各自的风险特征。比如光伏可以设Γ_pv6风电设Γ_wt10这样更贴近实际波动规律。6. 一点实际操作体会这套模型我完整复现过不止一次最深的感触是多阶段鲁棒调度的难点不在数学公式本身而在“把公式变成可求解的模型”这个过程。YALMIPCplex的组合让建模门槛降低了很多但CCG迭代、子问题对偶、大M线性化这些环节仍然需要耐心调试才能稳定跑出结果。如果你第一次跑这个模型就遇到不收敛或者结果离奇不用怀疑人生大概率是某个约束符号写反了或者大M取值不当。最后再分享一个小技巧把CCG迭代过程中每一轮的UB和LB打印出来画成收敛曲线。这条曲线不仅能直观看出算法的收敛速度还能帮你判断模型是不是出了问题——正常情况下LB和UB应该单调逼近如果出现回退说明生成割的过程中有严重错误赶紧回头检查子问题。这个习惯我建议所有做分解算法的人都培养起来。