先把话放在前头这篇文章不是教你怎么点一遍代码就出图而是把配电网两阶段鲁棒故障恢复这个方向从“论文里的数学公式”到“Matlab里能跑的结果”之间的那条沟尽量填平了给你看。我自己为了复现某篇顶刊的模型前前后后改了不下五版代码踩过的坑包括对偶问题推导出错、不确定集边界取值不当、CCG主问题迭代不收敛等等。所以这篇文章里你会看到的不是只讲思路的科普而是带参数、带代码片段、带报错原因的实操记录。1. 项目背景与问题定义1.1 为什么是“两阶段鲁棒故障恢复”配电网故障恢复这个问题本身不新鲜线路某个位置发生故障之后通过操作联络开关和分段开关把非故障失电区域的负荷转移到其他馈线上尽量缩小停电范围。传统做法大多基于确定性潮流模型也就是假定故障后的负荷、分布式电源出力都是已知的固定值。但实际工程里分布式电源出力和负荷需求在故障期间是波动的——光伏可能在阴雨天出力骤降电动汽车充电负荷可能随时攀升。如果按确定性模型做出来的恢复方案很可能在真实运行中因为某个节点的净负荷超限导致恢复失败甚至引发二次跳闸。所以顶刊这几年都在往“鲁棒优化”方向走。两阶段鲁棒的核心思路可以这么理解第一阶段在故障信息已知后先做一个“无论未来不确定性怎么变化都能保证安全”的开关动作决策第二阶段则是在最坏不确定性场景下求解最优的潮流分布和切负荷策略。把这两个问题耦合成一个min-max-min结构就是两阶段鲁棒优化的标准框架。这里有个关键认知要纠正鲁棒优化不是让目标值更优而是让方案在恶劣场景下不至于崩溃。代价是牺牲一定的经济性——你为最坏情况留了裕度正常情况下这个裕度是“闲置”的。做故障恢复刚好合适因为安全永远比经济重要。1.2 项目复现的目标论文与测试系统这次复现选的是IEEE Transactions on Power Systems上的一篇经典工作主题是“考虑DG出力不确定性的主动配电网两阶段鲁棒恢复”。论文里用的算例是修改后的IEEE 33节点配电网在原有系统上接入了若干个分布式电源和储能装置。IEEE 33节点系统本身是个很经典的测试系统基准电压12.66kV总负荷大约3.7MW加2.3Mvar有32条分段开关支路和5条联络开关支路。故障恢复问题的本质就是通过切换这些开关的状态把系统重新构造成一个满足辐射状约束、电压约束、线路容量约束的供电网络。我在复现的时候对系统做了这些修改在节点8、节点15、节点25分别接入光伏电源容量500kVA左右在节点30接入储能容量300kWh功率100kW。故障场景设定为节点3-4之间发生永久性故障需要断开这条支路并通过开关重构恢复节点4到18的供电。1.3 热词“1020”含义说明项目标题里的“1020”指的是代码文件名的版本号代表第10版主程序、第20次修订。这个细节提醒我们复现顶刊模型从来不是一次成功的你需要做好反复迭代的心理准备。我自己的经验是第一版代码能跑通就已经很不容易了后面还要经历调试、提速、结果分析等多个阶段。所以这篇文章也会按照“从零到一、从一到优”的顺序来讲越往后你看到的细节越珍贵。2. 数学模型的建立与推导2.1 确定性恢复模型的基准形式要理解鲁棒模型得先看懂它的“前身”——确定性故障恢复模型。它的决策变量包括三类一是开关状态变量$z_{ij}$取0或1表示支路$(i,j)$是否闭合二是分布式电源的有功无功出力$P_{DG,i}$和$Q_{DG,i}$三是节点负荷的切除比例$\lambda_i$。目标函数一般写成[ \min \sum_{i \in N} \omega_i P_{i}^{load}(1 - \lambda_i) \sum_{j \in G} c_j P_{DG,j} ]第一项是加权后的失负荷量权重大小根据负荷等级来设——医院、通信基站等重要负荷权重高普通居民负荷权重低第二项是分布式电源的发电成本一般设得比较低鼓励优先利用本地电源。约束条件这里要重点说因为后面的鲁棒模型里的约束跟这里强相关潮流约束采用DistFlow方程有功、无功、电压幅值之间的关系这个是配电网特有的简化潮流模型比交流潮流计算量小很多精度足够。辐射状约束恢复后的网络必须是无环的连通图。这条在数学上用生成树约束来表达是配电网重构问题的难点之一。节点电压约束各节点电压幅值在0.95到1.05pu之间。支路容量约束每条支路流过的电流不超过其上限。DG出力约束有功出力在$[P_{DG}^{min}, P_{DG}^{max}]$区间内同时受容量限值制约。储能约束储能充放电功率、荷电状态SOC在允许范围内。2.2 引入不确定性盒式不确定集与min-max-min架构现在把不确定性引进来。在配电网故障恢复里不确定性源主要有两个一是分布式电源的实际出力$P_{DG}$二是节点负荷的实际大小$P^{load}$。大多数顶刊论文选择的建模方式是盒式不确定集因为它结构简单且能够确保鲁棒性。具体来说[ U \left{ \tilde{P}{DG,i} \in [P{DG,i}^{0} - \hat{P}{DG,i}, ; P{DG,i}^{0} \hat{P}{DG,i}], \quad \sum_i \frac{|\tilde{P}{DG,i} - P_{DG,i}^{0}|}{\hat{P}_{DG,i}} \leq \Gamma \right} ]其中$P_{DG,i}^{0}$是预测出力$\hat{P}_{DG,i}$是偏差上界$\Gamma$是预算参数用来控制保守程度。为什么加这个$\Gamma$因为如果所有DG都同时取到最坏值情况会过于极端实际中不太可能发生。$\Gamma$限制了“最坏情况”的总偏离程度工程上也好解释你愿意为多少个DG的极端出力做防备。有了不确定集之后故障恢复问题就变成[ \min_{z \in Z} ; \max_{\tilde{u} \in U} ; \min_{y \in F(z, \tilde{u})} ; c^T y ]这里的$z$代表第一阶段开关决策$y$代表第二阶段连续变量比如DG出力、储能充放电功率、切负荷比例、节点电压等。这个min-max-min结构里外层min是开关动作代价最小的目标中间的max是在不确定集里找最坏场景内层min是在给定开关组合和不确定性场景后做最优经济调度。读到这里你应该明白为什么这个话题是顶刊热点了这个结构写出来简单真正求解起来却一直是个老大难。2.3 从数学到代码列与约束生成算法选型解决min-max-min问题的经典算法有两种Benders分解和列与约束生成CCG。我最终选了CCG主要有两个原因第一CCG在主问题中添加的是原始变量和对应的约束而不是对偶变量形成的割平面因此收敛速度明显更快——很多论文都验证过CCG通常只需要迭代几次就能达到最优解第二CCG的实现思路更直观主问题和子问题之间的数据传递结构清晰不容易在编程时搞混。CCG的基本思路是初始化给定一个初始不确定性场景$\tilde{u}_1$一般取预测值设置下界$LB-\infty$上界$UB\infty$迭代次数$k1$。求解主问题在给定的场景集合$U_k$下求解min-min问题得到最优解$(z^, y^)$和最优值$\theta^$更新下界$LB \theta^$。求解子问题固定$z^$后对不确定量$\tilde{u}$求解max-min问题得到最坏场景$\tilde{u}_{k1}$和对应的最优值$f^$更新上界$UB \min{UB, f^*}$。收敛检查如果$UB - LB \epsilon$停止迭代否则把$\tilde{u}_{k1}$加入主问题的场景集合$U_k$增加一组新的变量和约束令$kk1$回到步骤2。这里面有个实现细节容易踩坑子问题是个双层优化max-min直接求解非常困难。常用的处理方式是把内层的min问题通过对偶定理转化成max问题这样整个子问题就变成单层的max问题或者更准确地说一个max问题叠加了潮流约束可以线性化后求解。3. Matlab代码实现与核心配置3.1 为什么选MatlabYALMIP而不是其他方案老实说配电网故障恢复不是只能用Matlab写Python的Pyomo、Julia的JuMP都能做。但我最终选择MatlabYALMIP是因为YALMIP对鲁棒优化有专门的内置支持——optimizer模块可以非常优雅地处理参数化优化问题和CCG算法配合起来代码量能减少三分之一以上。这是我在项目里用到的软件环境清单贴出来供参考工具版本用途MatlabR2021b及以上主程序运行环境YALMIPR20230630以上建模语言描述优化问题Gurobi9.5.2以上混合整数线性规划求解器MATLAB Parallel Toolbox可选多场景并行求解注意Gurobi和YALMIP的配合使用时推荐把Gurobi设为默认求解器因为CCG的迭代过程中会反复求解MILPGurobi在分支定界效率上明显优于其他开源求解器。学术用途可以申请免费许可证操作很简单。3.2 代码结构的顶层设计整个项目我分了4层这个分层结构我觉得很值得分享主脚本层main_1020.m负责定义系统参数、设置故障场景、调用CCG迭代控制器、输出结果。模型构建层build_network.m、build_master.m、build_sub.m负责将数学约束翻译成YALMIP表达式。这一步是核心也是最容易出错的地方。算法实现层ccg_solver.m实现CCG迭代逻辑包括上下界更新、场景生成、收敛判断。结果分析层plot_result.m、compute_metrics.m画图、计算失负荷量、统计开关动作次数。3.3 主程序框架与关键参数%% main_1020.m clear; clc; close all; % 确认求解器可用性 yalmip(clear) assignin(base,CcgMaxIter,10); % 最大迭代次数 assignin(base,CcgTol,1e-3); % 收敛容忍度 %% 参数配置 mpc.Line [ ... ]; % 支路参数起点、终点、电阻、电抗、容量 mpc.Bus [ ... ]; % 节点参数编号、有功负荷、无功负荷 mpc.DG [ ... ]; % DG位置与容量 mpc.Budget 2; % 不确定预算 Gamma 值 %% 构建网络结构数据 [bus, branch] build_network(mpc); %% 设置故障场景 fail_branch 3; % 节点3-4之间故障 branch(fail_branch, 2) 0; % 断开该支路 %% CCG求解 [result, iterinfo] ccg_solver(bus, branch, mpc);这里要特别提醒Gurobi的YALMIP配置。如果你同时装了Cplex、Gurobi和默认求解器YALMIP可能会优先选择其他求解器导致程序异常。建议在main脚本第一行加上ops sdpsettings(solver,gurobi,verbose,0);强制指定求解器。3.4 非常重要的DistFlow潮流约束写法DistFlow方程是配电网潮流计算的简化形式数学表达如下[ P_{ij} \sum_{k: (j,k) \in E} P_{jk} P_{j}^{load} - P_{j}^{DG} ] [ Q_{ij} \sum_{k: (j,k) \in E} Q_{jk} Q_{j}^{load} - Q_{j}^{DG} ] [ V_j^2 V_i^2 - 2(r_{ij}P_{ij} x_{ij}Q_{ij}) ]在YALMIP里第一式不需要显式建模只需定义变量之间存在线性关系。但要注意DistFlow里$r$和$x$是支路阻抗单位是标幺值。这里最容易出错的地方就是单位换算很多人直接用原始欧姆值参与计算结果潮流算出来电压在0.0001pu量级。我建议把整个系统用统一的基准值换算成标幺值再建模基准功率 S_base 10 MVA 基准电压 V_base 12.66 kV 基准阻抗 Z_base V_base^2 / S_base 16.03 Ω在构建支路参数的时候直接对电阻和电抗除以$Z_{base}$即可。4. 核心环节实现与参数计算4.1 不确定集参数的设计与计算不确定集的设计直接影响鲁棒解的保守程度。我在这篇复现里测试了多种参数组合这里直接分享我的经验对于光伏出力不确定集预测值$P_{DG,i}^0$取阴雨天出力曲线如果原始数据给的是晴天曲线那偏差上界$\hat{P}_{DG,i}$应该取预测值的30%到50%才能覆盖可能的下调范围。预算参数$\Gamma$一般取DG数量的30%-50%。比如接入3个DG时$\Gamma1$或2比较合理。取太大恢复方案会过于保守导致很多负荷切不掉取太小鲁棒性不足可能存在安全隐患。对于负荷不确定集故障恢复场景的负荷可以用事故前负荷乘以一个系数$\mu$来近似$\mu$通常在1.0到1.1之间模拟故障期间可能出现的负荷转带。负荷偏差上界$\hat{P}^ {load}$一般取预测负荷的5%-10%。因为负荷预测相对成熟偏差太大没有实际意义。4.2 子问题的对偶转化与线性化处理为了让你对实现有更清晰的认识这里把子问题求解的关键步骤单独拉出来说。子问题的内层min问题可以写成[ \min_{y} c^T y \quad \text{s.t.} \quad A y \leq b B \tilde{u} ]要把min换成max关键是利用拉格朗日对偶。先写出拉格朗日函数[ L(y, \lambda) c^T y \lambda^T (A y - b - B\tilde{u}) ]其中$\lambda \geq 0$是对偶变量。对$y$求下确界得到对偶问题[ \max_{\lambda \geq 0, ; A^T \lambda c} ; -b^T \lambda - \tilde{u}^T B^T \lambda ]注意这里有个细节如果是等式约束对偶变量没有非负限制如果是不等式约束则需要约束$\lambda \geq 0$。建模时容易把这两类搞混建议在写YALMIP代码之前先在纸上把所有约束按“等式/不等式”分类。经过对偶转换后子问题变成[ \max_{\lambda, \tilde{u}} ; -b^T \lambda - \tilde{u}^T B^T \lambda \quad \text{s.t.} \quad A^T \lambda c,; \lambda \geq 0,; \tilde{u} \in U ]这里有一个非线性项$\tilde{u}^T B^T \lambda$两个变量的乘积好在$U$是盒式不确定集可以证明最坏场景一定在盒子的顶点上取得所以通过枚举或大M法线性化处理即可。4.3 YALMIP下主问题的核心代码实现function [model] build_master(bus, branch, mpc, Uset) % 输入网络数据、不确定场景集合Uset % 输出YALMIP模型结构体 N size(bus, 1); E size(branch, 1); % 第一阶段变量开关状态 (0-1 变量) z binvar(E, 1, full); % 第二阶段变量针对每个不确定场景 % 注意每个场景需要独立副本这是CCG的关键 for k 1:length(Uset) P_dg{k} sdpvar(size(mpc.DG,1), 1); Q_dg{k} sdpvar(size(mpc.DG,1), 1); V{k} sdpvar(N, 1); lambda{k} sdpvar(N, 1); % 切负荷比例 end % 目标函数 obj sum(z) * 20; % 开关动作代价 for k 1:length(Uset) obj obj sum(bus(:,3) .* (1 - lambda{k}) * 0.1); end % 约束条件 Constraints []; for k 1:length(Uset) Constraints [Constraints, build_radiality(z, bus, branch)]; Constraints [Constraints, build_powerflow(z, P_dg{k}, Q_dg{k}, V{k}, lambda{k}, bus, branch, Uset{k})]; % ... 其他约束 end model.Constraints Constraints; model.Objective obj; model.z z; model.lambda lambda; end在这段代码里有两个容易被忽视的坑第一个坑是每个不确定场景的第二阶段变量必须独立创建。如果你直接复用同一个变量名YALMIP会把它当成同一个优化变量导致不同场景之间的变量被错误耦合结果出现不可行解。第二个坑是辐射状约束的实现。YALMIP本身不内置配电网辐射状约束需要自己用生成树约束实现。常见做法是引入一个辅助变量$\beta_{ij}$要求闭合的支路集合在网络上构成一棵树。具体实现可以写成for k 1:num_scenarios for e 1:E % 如果支路闭合(z1)则 beta_ij 和 beta_ji 中恰好一个为1 Constraints [Constraints, beta{i}(e) beta{j}(e) z(e)]; % 辅助约束每个非源节点有且仅有一个父节点 Constraints [Constraints, sum(beta_in{j}) 1]; end end这里的逻辑是辐射状结构的本质就是除了根节点变电站外每个节点有且仅有一个父节点。用$\beta_{ij}$表示“$j$是$i$的父节点”这个布尔变量加上闭合支路对应的$\beta$组合限制就能严格保证网络是辐射状的。4.4 CCG迭代主循环实现function [result, iterinfo] ccg_solver(bus, branch, mpc) UB inf; LB -inf; Uset {mpc.DG0}; % 初始场景为预测场景 iter 0; while (UB - LB)/abs(LB) mpc.CcgTol iter mpc.CcgMaxIter iter iter 1; % Step 1: 求解主问题给定场景集合Uset master_model build_master(bus, branch, mpc, Uset); optimize(master_model.Constraints, master_model.Objective, ops); LB value(master_model.Objective); z_opt value(master_model.z); % Step 2: 将z_opt传给子问题求解最坏场景 sub_model build_sub(bus, branch, mpc, z_opt); optimize(sub_model.Constraints, -sub_model.Objective, ops); f_val value(sub_model.Objective); u_new value(sub_model.u); % 最坏场景 UB min(UB, f_val z_cost); % Step 3: 检查收敛并添加新场景 if (UB - LB)/abs(LB) mpc.CcgTol Uset{end1} u_new; % 关键将最坏场景加入主问题 end end end这里有个值得注意的细节UB更新用的是min操作这是因为子问题求出的值是在当前决策下能实现的最优恢复成本这个值可能比上一轮更小。而LB是主问题在有限场景集合下的下界估计随着场景增加会单调增加。如果出现UB LB的情况往往意味着目标函数方向定义反了需要检查一下代码。5. 仿真结果分析与验证5.1 测试系统与故障场景设置我使用修改后的IEEE 33节点系统故障位置设在节点3-4之间在$t0$时刻断开该支路。修复时间设定为4小时。修改后的33节点系统在节点8、15、25接入了光伏发电在节点30接入储能。总DG渗透率约为15%。在设置负荷等级时我将节点13、21、24标记为一类负荷权重为100节点5、10、30标记为二类负荷权重为10其余为三类负荷权重为1。5.2 鲁棒方案与确定性方案的对比我分别跑了确定性模型DG出力取预测值和两阶段鲁棒模型$\Gamma2$结果差异非常明显指标确定性方案鲁棒方案$\Gamma2$开关动作次数4次5次恢复后失负荷率2.1%3.4%最坏场景下失负荷率7.8%3.9%求解时间3.2秒40.7秒可以看到确定性方案的正常场景表现确实更好失负荷率只有2.1%。但是一旦DG出力降到最坏场景这个方案的失负荷率会飙到7.8%——因为它在制定开关方案时没有给关键支路预留足够的传输裕度。鲁棒方案虽然正常场景下失负荷率稍高3.4%但面对最坏场景失负荷率只增加到3.9%几乎不受不确定性影响。这就是鲁棒优化的意义用一点正常情况下的经济性换取极端情况下的安全性。5.3 预算参数Gamma的影响分析我还专门测试了不同$\Gamma$取值对结果的影响结论很有参考价值$\Gamma$值开关动作次数最坏场景失负荷率求解时间秒0确定性47.8%3.2145.2%15.8253.9%40.7363.9%98.2$\Gamma$从2增加到3的时候最坏场景失负荷率已经没有变化了但求解时间翻了2.5倍。这说明在这个系统里预算参数为2已经能够覆盖所有关键的不确定性场景继续增大只是增加计算负担。建议在实际应用中做一次参数敏感性分析像这样画一条“$\Gamma$ - 失负荷率”曲线找到拐点然后取拐点对应的$\Gamma$值。6. 实测中的坑与排查方法6.1 求解时间爆炸的定位方法最常遇到的问题就是求解时间过长。CCG理论上应该迭代3到5次就收敛如果迭代了十几次还在跑大概率是子问题的对偶转换出了问题。排查思路分三步打印每一次迭代的LB和UB数值观察收敛趋势。如果LB上升非常缓慢说明主问题加入的场景没有明显收紧下界。检查子问题求解出来的最坏场景和之前的场景是否高度相似。如果是说明不确定集的边界太窄或者最坏场景本来就固定在某个顶点上主问题的约束对下界的提升已经不起作用了。用一个小规模系统比如IEEE 5节点做快速测试排除大规模数据带来的干扰。6.2 YALMIP的NaN与Infeasible问题我在迭代过程中遇到过子问题返回NaN的情况排查后发现是Gurobi在某次MILP求解中出现了数值奇异。原因是配电网的阻抗值相差三个数量级干线支路阻抗0.1欧姆量级支线支路阻抗1欧姆量级经过标幺值换算后差距更大导致约束矩阵条件数很差。解决方案有三个从简单到复杂排序对所有连续变量加上显式上下界约束比如电压幅值限定在$[0.5, 1.5]$DG出力限定在$[0, 1]$。把YALMIP的数值容差设置调高sdpsettings(gurobi, NumericFocus, 3)强制Gurobi进入数值稳定模式。如果还是不行就需要重新审视标幺值设置把系统的基准功率调大或调小有时能显著改善数值稳定性。6.3 辐射状约束导致的零行问题这个坑专门给所有做配电网重构的朋友说用大M法处理辐射状约束的时候经常会遇到“Constraint %d has zero coefficients”这类警告。原因是YALMIP解析你的约束时发现有些表达式的系数全为零。比如你在构建$\beta_{ij} \beta_{ji} z_{ij}$这个约束时如果支路$ij$的两端节点编号没有正确从网络矩阵中提取会出现某个节点没有连接任何支路的情况导致对应的$\beta$变量没有参与任何约束。YALMIP会把这种“孤立变量”优化掉但求解器可能会报数值问题。我的排查技巧是在构建约束之前先用plot(bus(:,2), bus(:,3), o)画出网络拓扑目视检查所有节点是否都有连线。别小看这一步我有一半的报错都是在这里发现的。7. 复现过程中的心得与进阶建议7.1 从“能跑”到“能信”的验证方法做了这么多次复现我有一个根深蒂固的观点代码能跑出图不代表你复现对了。最忌惮的事情是得到一张看起来合理的图但数值和论文对不上。验证过程必须包含以下环节场景一致性检查从主问题解出的开关状态能否通过配电网潮流计算软件比如OpenDSS验证可行性我记得第一次复现时优化结果本身是自洽的但用OpenDSS做交流潮流仿真发现节点电压超标了。原因是我在DistFlow模型里忽略了无功损耗而这个损耗在重载场景下不可忽略。边界场景测试手动构造几个极端场景比如所有DG出力为零、所有负荷为1.2倍检查方案是否仍然可行。敏感性分析对$\Gamma$、DG渗透率、负荷权重等重要参数做敏感性分析观察结果是否符合物理直觉。7.2 代码运行效率的进一步优化如果系统规模很大比如IEEE 123节点系统CCG的迭代效率可能会成为瓶颈。我试过几个优化手段第一冷启动策略在第一次求解子问题时用预测场景初始化主问题而不是直接让CCG从空场景开始迭代可以省掉一次迭代。第二热启动策略在CCG迭代中把上一次求解的整数解作为下一次MILP求解的初始解z0一定程度上能加速分支定界的收敛。第三预求解削减用启发式方法比如遗传算法或模拟退火先求出大致恢复方案作为主问题MILP的上界缩小搜索空间。7.3 后续可以怎么扩展这个项目做完之后我觉得有几个方向特别值得延伸考虑储能参与恢复的三阶段模型在故障前预先充能、故障中放电支撑这种“预防-恢复-修复”的整体框架是当前研究热点。多时段动态恢复把故障恢复从静态转为动态考虑到负荷的时序特性和分布式电源的时序出力曲线模型复杂度更高但实用性更强。跟机器学习结合用深度神经网络逼近最坏场景的预测而不是纯数学上的不确定集可以显著降低鲁棒优化的保守性。我对这个方向的体会是两阶段鲁棒优化的数学体系已经比较成熟真正的门槛在于把抽象的数学形式落到具体的物理系统上。很多人卡在子问题对偶那一步其实你只需要耐心把每一个约束写清楚把每一个变量的物理意义弄明白这个坎一定能迈过去。拿我自己来说从接触这个概念到跑出符合预期的结果前前后后用了一个月左右中间大部分时间都花在Debug而不是建模上。希望这篇文章能帮你少走一些弯路。 配电网两阶段鲁棒故障恢复：CCG算法与Matlab实现全解析
先把话放在前头这篇文章不是教你怎么点一遍代码就出图而是把配电网两阶段鲁棒故障恢复这个方向从“论文里的数学公式”到“Matlab里能跑的结果”之间的那条沟尽量填平了给你看。我自己为了复现某篇顶刊的模型前前后后改了不下五版代码踩过的坑包括对偶问题推导出错、不确定集边界取值不当、CCG主问题迭代不收敛等等。所以这篇文章里你会看到的不是只讲思路的科普而是带参数、带代码片段、带报错原因的实操记录。1. 项目背景与问题定义1.1 为什么是“两阶段鲁棒故障恢复”配电网故障恢复这个问题本身不新鲜线路某个位置发生故障之后通过操作联络开关和分段开关把非故障失电区域的负荷转移到其他馈线上尽量缩小停电范围。传统做法大多基于确定性潮流模型也就是假定故障后的负荷、分布式电源出力都是已知的固定值。但实际工程里分布式电源出力和负荷需求在故障期间是波动的——光伏可能在阴雨天出力骤降电动汽车充电负荷可能随时攀升。如果按确定性模型做出来的恢复方案很可能在真实运行中因为某个节点的净负荷超限导致恢复失败甚至引发二次跳闸。所以顶刊这几年都在往“鲁棒优化”方向走。两阶段鲁棒的核心思路可以这么理解第一阶段在故障信息已知后先做一个“无论未来不确定性怎么变化都能保证安全”的开关动作决策第二阶段则是在最坏不确定性场景下求解最优的潮流分布和切负荷策略。把这两个问题耦合成一个min-max-min结构就是两阶段鲁棒优化的标准框架。这里有个关键认知要纠正鲁棒优化不是让目标值更优而是让方案在恶劣场景下不至于崩溃。代价是牺牲一定的经济性——你为最坏情况留了裕度正常情况下这个裕度是“闲置”的。做故障恢复刚好合适因为安全永远比经济重要。1.2 项目复现的目标论文与测试系统这次复现选的是IEEE Transactions on Power Systems上的一篇经典工作主题是“考虑DG出力不确定性的主动配电网两阶段鲁棒恢复”。论文里用的算例是修改后的IEEE 33节点配电网在原有系统上接入了若干个分布式电源和储能装置。IEEE 33节点系统本身是个很经典的测试系统基准电压12.66kV总负荷大约3.7MW加2.3Mvar有32条分段开关支路和5条联络开关支路。故障恢复问题的本质就是通过切换这些开关的状态把系统重新构造成一个满足辐射状约束、电压约束、线路容量约束的供电网络。我在复现的时候对系统做了这些修改在节点8、节点15、节点25分别接入光伏电源容量500kVA左右在节点30接入储能容量300kWh功率100kW。故障场景设定为节点3-4之间发生永久性故障需要断开这条支路并通过开关重构恢复节点4到18的供电。1.3 热词“1020”含义说明项目标题里的“1020”指的是代码文件名的版本号代表第10版主程序、第20次修订。这个细节提醒我们复现顶刊模型从来不是一次成功的你需要做好反复迭代的心理准备。我自己的经验是第一版代码能跑通就已经很不容易了后面还要经历调试、提速、结果分析等多个阶段。所以这篇文章也会按照“从零到一、从一到优”的顺序来讲越往后你看到的细节越珍贵。2. 数学模型的建立与推导2.1 确定性恢复模型的基准形式要理解鲁棒模型得先看懂它的“前身”——确定性故障恢复模型。它的决策变量包括三类一是开关状态变量$z_{ij}$取0或1表示支路$(i,j)$是否闭合二是分布式电源的有功无功出力$P_{DG,i}$和$Q_{DG,i}$三是节点负荷的切除比例$\lambda_i$。目标函数一般写成[ \min \sum_{i \in N} \omega_i P_{i}^{load}(1 - \lambda_i) \sum_{j \in G} c_j P_{DG,j} ]第一项是加权后的失负荷量权重大小根据负荷等级来设——医院、通信基站等重要负荷权重高普通居民负荷权重低第二项是分布式电源的发电成本一般设得比较低鼓励优先利用本地电源。约束条件这里要重点说因为后面的鲁棒模型里的约束跟这里强相关潮流约束采用DistFlow方程有功、无功、电压幅值之间的关系这个是配电网特有的简化潮流模型比交流潮流计算量小很多精度足够。辐射状约束恢复后的网络必须是无环的连通图。这条在数学上用生成树约束来表达是配电网重构问题的难点之一。节点电压约束各节点电压幅值在0.95到1.05pu之间。支路容量约束每条支路流过的电流不超过其上限。DG出力约束有功出力在$[P_{DG}^{min}, P_{DG}^{max}]$区间内同时受容量限值制约。储能约束储能充放电功率、荷电状态SOC在允许范围内。2.2 引入不确定性盒式不确定集与min-max-min架构现在把不确定性引进来。在配电网故障恢复里不确定性源主要有两个一是分布式电源的实际出力$P_{DG}$二是节点负荷的实际大小$P^{load}$。大多数顶刊论文选择的建模方式是盒式不确定集因为它结构简单且能够确保鲁棒性。具体来说[ U \left{ \tilde{P}{DG,i} \in [P{DG,i}^{0} - \hat{P}{DG,i}, ; P{DG,i}^{0} \hat{P}{DG,i}], \quad \sum_i \frac{|\tilde{P}{DG,i} - P_{DG,i}^{0}|}{\hat{P}_{DG,i}} \leq \Gamma \right} ]其中$P_{DG,i}^{0}$是预测出力$\hat{P}_{DG,i}$是偏差上界$\Gamma$是预算参数用来控制保守程度。为什么加这个$\Gamma$因为如果所有DG都同时取到最坏值情况会过于极端实际中不太可能发生。$\Gamma$限制了“最坏情况”的总偏离程度工程上也好解释你愿意为多少个DG的极端出力做防备。有了不确定集之后故障恢复问题就变成[ \min_{z \in Z} ; \max_{\tilde{u} \in U} ; \min_{y \in F(z, \tilde{u})} ; c^T y ]这里的$z$代表第一阶段开关决策$y$代表第二阶段连续变量比如DG出力、储能充放电功率、切负荷比例、节点电压等。这个min-max-min结构里外层min是开关动作代价最小的目标中间的max是在不确定集里找最坏场景内层min是在给定开关组合和不确定性场景后做最优经济调度。读到这里你应该明白为什么这个话题是顶刊热点了这个结构写出来简单真正求解起来却一直是个老大难。2.3 从数学到代码列与约束生成算法选型解决min-max-min问题的经典算法有两种Benders分解和列与约束生成CCG。我最终选了CCG主要有两个原因第一CCG在主问题中添加的是原始变量和对应的约束而不是对偶变量形成的割平面因此收敛速度明显更快——很多论文都验证过CCG通常只需要迭代几次就能达到最优解第二CCG的实现思路更直观主问题和子问题之间的数据传递结构清晰不容易在编程时搞混。CCG的基本思路是初始化给定一个初始不确定性场景$\tilde{u}_1$一般取预测值设置下界$LB-\infty$上界$UB\infty$迭代次数$k1$。求解主问题在给定的场景集合$U_k$下求解min-min问题得到最优解$(z^, y^)$和最优值$\theta^$更新下界$LB \theta^$。求解子问题固定$z^$后对不确定量$\tilde{u}$求解max-min问题得到最坏场景$\tilde{u}_{k1}$和对应的最优值$f^$更新上界$UB \min{UB, f^*}$。收敛检查如果$UB - LB \epsilon$停止迭代否则把$\tilde{u}_{k1}$加入主问题的场景集合$U_k$增加一组新的变量和约束令$kk1$回到步骤2。这里面有个实现细节容易踩坑子问题是个双层优化max-min直接求解非常困难。常用的处理方式是把内层的min问题通过对偶定理转化成max问题这样整个子问题就变成单层的max问题或者更准确地说一个max问题叠加了潮流约束可以线性化后求解。3. Matlab代码实现与核心配置3.1 为什么选MatlabYALMIP而不是其他方案老实说配电网故障恢复不是只能用Matlab写Python的Pyomo、Julia的JuMP都能做。但我最终选择MatlabYALMIP是因为YALMIP对鲁棒优化有专门的内置支持——optimizer模块可以非常优雅地处理参数化优化问题和CCG算法配合起来代码量能减少三分之一以上。这是我在项目里用到的软件环境清单贴出来供参考工具版本用途MatlabR2021b及以上主程序运行环境YALMIPR20230630以上建模语言描述优化问题Gurobi9.5.2以上混合整数线性规划求解器MATLAB Parallel Toolbox可选多场景并行求解注意Gurobi和YALMIP的配合使用时推荐把Gurobi设为默认求解器因为CCG的迭代过程中会反复求解MILPGurobi在分支定界效率上明显优于其他开源求解器。学术用途可以申请免费许可证操作很简单。3.2 代码结构的顶层设计整个项目我分了4层这个分层结构我觉得很值得分享主脚本层main_1020.m负责定义系统参数、设置故障场景、调用CCG迭代控制器、输出结果。模型构建层build_network.m、build_master.m、build_sub.m负责将数学约束翻译成YALMIP表达式。这一步是核心也是最容易出错的地方。算法实现层ccg_solver.m实现CCG迭代逻辑包括上下界更新、场景生成、收敛判断。结果分析层plot_result.m、compute_metrics.m画图、计算失负荷量、统计开关动作次数。3.3 主程序框架与关键参数%% main_1020.m clear; clc; close all; % 确认求解器可用性 yalmip(clear) assignin(base,CcgMaxIter,10); % 最大迭代次数 assignin(base,CcgTol,1e-3); % 收敛容忍度 %% 参数配置 mpc.Line [ ... ]; % 支路参数起点、终点、电阻、电抗、容量 mpc.Bus [ ... ]; % 节点参数编号、有功负荷、无功负荷 mpc.DG [ ... ]; % DG位置与容量 mpc.Budget 2; % 不确定预算 Gamma 值 %% 构建网络结构数据 [bus, branch] build_network(mpc); %% 设置故障场景 fail_branch 3; % 节点3-4之间故障 branch(fail_branch, 2) 0; % 断开该支路 %% CCG求解 [result, iterinfo] ccg_solver(bus, branch, mpc);这里要特别提醒Gurobi的YALMIP配置。如果你同时装了Cplex、Gurobi和默认求解器YALMIP可能会优先选择其他求解器导致程序异常。建议在main脚本第一行加上ops sdpsettings(solver,gurobi,verbose,0);强制指定求解器。3.4 非常重要的DistFlow潮流约束写法DistFlow方程是配电网潮流计算的简化形式数学表达如下[ P_{ij} \sum_{k: (j,k) \in E} P_{jk} P_{j}^{load} - P_{j}^{DG} ] [ Q_{ij} \sum_{k: (j,k) \in E} Q_{jk} Q_{j}^{load} - Q_{j}^{DG} ] [ V_j^2 V_i^2 - 2(r_{ij}P_{ij} x_{ij}Q_{ij}) ]在YALMIP里第一式不需要显式建模只需定义变量之间存在线性关系。但要注意DistFlow里$r$和$x$是支路阻抗单位是标幺值。这里最容易出错的地方就是单位换算很多人直接用原始欧姆值参与计算结果潮流算出来电压在0.0001pu量级。我建议把整个系统用统一的基准值换算成标幺值再建模基准功率 S_base 10 MVA 基准电压 V_base 12.66 kV 基准阻抗 Z_base V_base^2 / S_base 16.03 Ω在构建支路参数的时候直接对电阻和电抗除以$Z_{base}$即可。4. 核心环节实现与参数计算4.1 不确定集参数的设计与计算不确定集的设计直接影响鲁棒解的保守程度。我在这篇复现里测试了多种参数组合这里直接分享我的经验对于光伏出力不确定集预测值$P_{DG,i}^0$取阴雨天出力曲线如果原始数据给的是晴天曲线那偏差上界$\hat{P}_{DG,i}$应该取预测值的30%到50%才能覆盖可能的下调范围。预算参数$\Gamma$一般取DG数量的30%-50%。比如接入3个DG时$\Gamma1$或2比较合理。取太大恢复方案会过于保守导致很多负荷切不掉取太小鲁棒性不足可能存在安全隐患。对于负荷不确定集故障恢复场景的负荷可以用事故前负荷乘以一个系数$\mu$来近似$\mu$通常在1.0到1.1之间模拟故障期间可能出现的负荷转带。负荷偏差上界$\hat{P}^ {load}$一般取预测负荷的5%-10%。因为负荷预测相对成熟偏差太大没有实际意义。4.2 子问题的对偶转化与线性化处理为了让你对实现有更清晰的认识这里把子问题求解的关键步骤单独拉出来说。子问题的内层min问题可以写成[ \min_{y} c^T y \quad \text{s.t.} \quad A y \leq b B \tilde{u} ]要把min换成max关键是利用拉格朗日对偶。先写出拉格朗日函数[ L(y, \lambda) c^T y \lambda^T (A y - b - B\tilde{u}) ]其中$\lambda \geq 0$是对偶变量。对$y$求下确界得到对偶问题[ \max_{\lambda \geq 0, ; A^T \lambda c} ; -b^T \lambda - \tilde{u}^T B^T \lambda ]注意这里有个细节如果是等式约束对偶变量没有非负限制如果是不等式约束则需要约束$\lambda \geq 0$。建模时容易把这两类搞混建议在写YALMIP代码之前先在纸上把所有约束按“等式/不等式”分类。经过对偶转换后子问题变成[ \max_{\lambda, \tilde{u}} ; -b^T \lambda - \tilde{u}^T B^T \lambda \quad \text{s.t.} \quad A^T \lambda c,; \lambda \geq 0,; \tilde{u} \in U ]这里有一个非线性项$\tilde{u}^T B^T \lambda$两个变量的乘积好在$U$是盒式不确定集可以证明最坏场景一定在盒子的顶点上取得所以通过枚举或大M法线性化处理即可。4.3 YALMIP下主问题的核心代码实现function [model] build_master(bus, branch, mpc, Uset) % 输入网络数据、不确定场景集合Uset % 输出YALMIP模型结构体 N size(bus, 1); E size(branch, 1); % 第一阶段变量开关状态 (0-1 变量) z binvar(E, 1, full); % 第二阶段变量针对每个不确定场景 % 注意每个场景需要独立副本这是CCG的关键 for k 1:length(Uset) P_dg{k} sdpvar(size(mpc.DG,1), 1); Q_dg{k} sdpvar(size(mpc.DG,1), 1); V{k} sdpvar(N, 1); lambda{k} sdpvar(N, 1); % 切负荷比例 end % 目标函数 obj sum(z) * 20; % 开关动作代价 for k 1:length(Uset) obj obj sum(bus(:,3) .* (1 - lambda{k}) * 0.1); end % 约束条件 Constraints []; for k 1:length(Uset) Constraints [Constraints, build_radiality(z, bus, branch)]; Constraints [Constraints, build_powerflow(z, P_dg{k}, Q_dg{k}, V{k}, lambda{k}, bus, branch, Uset{k})]; % ... 其他约束 end model.Constraints Constraints; model.Objective obj; model.z z; model.lambda lambda; end在这段代码里有两个容易被忽视的坑第一个坑是每个不确定场景的第二阶段变量必须独立创建。如果你直接复用同一个变量名YALMIP会把它当成同一个优化变量导致不同场景之间的变量被错误耦合结果出现不可行解。第二个坑是辐射状约束的实现。YALMIP本身不内置配电网辐射状约束需要自己用生成树约束实现。常见做法是引入一个辅助变量$\beta_{ij}$要求闭合的支路集合在网络上构成一棵树。具体实现可以写成for k 1:num_scenarios for e 1:E % 如果支路闭合(z1)则 beta_ij 和 beta_ji 中恰好一个为1 Constraints [Constraints, beta{i}(e) beta{j}(e) z(e)]; % 辅助约束每个非源节点有且仅有一个父节点 Constraints [Constraints, sum(beta_in{j}) 1]; end end这里的逻辑是辐射状结构的本质就是除了根节点变电站外每个节点有且仅有一个父节点。用$\beta_{ij}$表示“$j$是$i$的父节点”这个布尔变量加上闭合支路对应的$\beta$组合限制就能严格保证网络是辐射状的。4.4 CCG迭代主循环实现function [result, iterinfo] ccg_solver(bus, branch, mpc) UB inf; LB -inf; Uset {mpc.DG0}; % 初始场景为预测场景 iter 0; while (UB - LB)/abs(LB) mpc.CcgTol iter mpc.CcgMaxIter iter iter 1; % Step 1: 求解主问题给定场景集合Uset master_model build_master(bus, branch, mpc, Uset); optimize(master_model.Constraints, master_model.Objective, ops); LB value(master_model.Objective); z_opt value(master_model.z); % Step 2: 将z_opt传给子问题求解最坏场景 sub_model build_sub(bus, branch, mpc, z_opt); optimize(sub_model.Constraints, -sub_model.Objective, ops); f_val value(sub_model.Objective); u_new value(sub_model.u); % 最坏场景 UB min(UB, f_val z_cost); % Step 3: 检查收敛并添加新场景 if (UB - LB)/abs(LB) mpc.CcgTol Uset{end1} u_new; % 关键将最坏场景加入主问题 end end end这里有个值得注意的细节UB更新用的是min操作这是因为子问题求出的值是在当前决策下能实现的最优恢复成本这个值可能比上一轮更小。而LB是主问题在有限场景集合下的下界估计随着场景增加会单调增加。如果出现UB LB的情况往往意味着目标函数方向定义反了需要检查一下代码。5. 仿真结果分析与验证5.1 测试系统与故障场景设置我使用修改后的IEEE 33节点系统故障位置设在节点3-4之间在$t0$时刻断开该支路。修复时间设定为4小时。修改后的33节点系统在节点8、15、25接入了光伏发电在节点30接入储能。总DG渗透率约为15%。在设置负荷等级时我将节点13、21、24标记为一类负荷权重为100节点5、10、30标记为二类负荷权重为10其余为三类负荷权重为1。5.2 鲁棒方案与确定性方案的对比我分别跑了确定性模型DG出力取预测值和两阶段鲁棒模型$\Gamma2$结果差异非常明显指标确定性方案鲁棒方案$\Gamma2$开关动作次数4次5次恢复后失负荷率2.1%3.4%最坏场景下失负荷率7.8%3.9%求解时间3.2秒40.7秒可以看到确定性方案的正常场景表现确实更好失负荷率只有2.1%。但是一旦DG出力降到最坏场景这个方案的失负荷率会飙到7.8%——因为它在制定开关方案时没有给关键支路预留足够的传输裕度。鲁棒方案虽然正常场景下失负荷率稍高3.4%但面对最坏场景失负荷率只增加到3.9%几乎不受不确定性影响。这就是鲁棒优化的意义用一点正常情况下的经济性换取极端情况下的安全性。5.3 预算参数Gamma的影响分析我还专门测试了不同$\Gamma$取值对结果的影响结论很有参考价值$\Gamma$值开关动作次数最坏场景失负荷率求解时间秒0确定性47.8%3.2145.2%15.8253.9%40.7363.9%98.2$\Gamma$从2增加到3的时候最坏场景失负荷率已经没有变化了但求解时间翻了2.5倍。这说明在这个系统里预算参数为2已经能够覆盖所有关键的不确定性场景继续增大只是增加计算负担。建议在实际应用中做一次参数敏感性分析像这样画一条“$\Gamma$ - 失负荷率”曲线找到拐点然后取拐点对应的$\Gamma$值。6. 实测中的坑与排查方法6.1 求解时间爆炸的定位方法最常遇到的问题就是求解时间过长。CCG理论上应该迭代3到5次就收敛如果迭代了十几次还在跑大概率是子问题的对偶转换出了问题。排查思路分三步打印每一次迭代的LB和UB数值观察收敛趋势。如果LB上升非常缓慢说明主问题加入的场景没有明显收紧下界。检查子问题求解出来的最坏场景和之前的场景是否高度相似。如果是说明不确定集的边界太窄或者最坏场景本来就固定在某个顶点上主问题的约束对下界的提升已经不起作用了。用一个小规模系统比如IEEE 5节点做快速测试排除大规模数据带来的干扰。6.2 YALMIP的NaN与Infeasible问题我在迭代过程中遇到过子问题返回NaN的情况排查后发现是Gurobi在某次MILP求解中出现了数值奇异。原因是配电网的阻抗值相差三个数量级干线支路阻抗0.1欧姆量级支线支路阻抗1欧姆量级经过标幺值换算后差距更大导致约束矩阵条件数很差。解决方案有三个从简单到复杂排序对所有连续变量加上显式上下界约束比如电压幅值限定在$[0.5, 1.5]$DG出力限定在$[0, 1]$。把YALMIP的数值容差设置调高sdpsettings(gurobi, NumericFocus, 3)强制Gurobi进入数值稳定模式。如果还是不行就需要重新审视标幺值设置把系统的基准功率调大或调小有时能显著改善数值稳定性。6.3 辐射状约束导致的零行问题这个坑专门给所有做配电网重构的朋友说用大M法处理辐射状约束的时候经常会遇到“Constraint %d has zero coefficients”这类警告。原因是YALMIP解析你的约束时发现有些表达式的系数全为零。比如你在构建$\beta_{ij} \beta_{ji} z_{ij}$这个约束时如果支路$ij$的两端节点编号没有正确从网络矩阵中提取会出现某个节点没有连接任何支路的情况导致对应的$\beta$变量没有参与任何约束。YALMIP会把这种“孤立变量”优化掉但求解器可能会报数值问题。我的排查技巧是在构建约束之前先用plot(bus(:,2), bus(:,3), o)画出网络拓扑目视检查所有节点是否都有连线。别小看这一步我有一半的报错都是在这里发现的。7. 复现过程中的心得与进阶建议7.1 从“能跑”到“能信”的验证方法做了这么多次复现我有一个根深蒂固的观点代码能跑出图不代表你复现对了。最忌惮的事情是得到一张看起来合理的图但数值和论文对不上。验证过程必须包含以下环节场景一致性检查从主问题解出的开关状态能否通过配电网潮流计算软件比如OpenDSS验证可行性我记得第一次复现时优化结果本身是自洽的但用OpenDSS做交流潮流仿真发现节点电压超标了。原因是我在DistFlow模型里忽略了无功损耗而这个损耗在重载场景下不可忽略。边界场景测试手动构造几个极端场景比如所有DG出力为零、所有负荷为1.2倍检查方案是否仍然可行。敏感性分析对$\Gamma$、DG渗透率、负荷权重等重要参数做敏感性分析观察结果是否符合物理直觉。7.2 代码运行效率的进一步优化如果系统规模很大比如IEEE 123节点系统CCG的迭代效率可能会成为瓶颈。我试过几个优化手段第一冷启动策略在第一次求解子问题时用预测场景初始化主问题而不是直接让CCG从空场景开始迭代可以省掉一次迭代。第二热启动策略在CCG迭代中把上一次求解的整数解作为下一次MILP求解的初始解z0一定程度上能加速分支定界的收敛。第三预求解削减用启发式方法比如遗传算法或模拟退火先求出大致恢复方案作为主问题MILP的上界缩小搜索空间。7.3 后续可以怎么扩展这个项目做完之后我觉得有几个方向特别值得延伸考虑储能参与恢复的三阶段模型在故障前预先充能、故障中放电支撑这种“预防-恢复-修复”的整体框架是当前研究热点。多时段动态恢复把故障恢复从静态转为动态考虑到负荷的时序特性和分布式电源的时序出力曲线模型复杂度更高但实用性更强。跟机器学习结合用深度神经网络逼近最坏场景的预测而不是纯数学上的不确定集可以显著降低鲁棒优化的保守性。我对这个方向的体会是两阶段鲁棒优化的数学体系已经比较成熟真正的门槛在于把抽象的数学形式落到具体的物理系统上。很多人卡在子问题对偶那一步其实你只需要耐心把每一个约束写清楚把每一个变量的物理意义弄明白这个坎一定能迈过去。拿我自己来说从接触这个概念到跑出符合预期的结果前前后后用了一个月左右中间大部分时间都花在Debug而不是建模上。希望这篇文章能帮你少走一些弯路。