做主动配电网故障恢复研究的同学,十有八九都要跟两个词打交道:重构和孤岛划分。传统做法是先把网络重构问题单独跑一遍,再把孤岛划分拎出来算,结果一到联动仿真就露馅——重构后的拓扑把本来可以孤岛运行的分布式电源区域切断了,或者孤岛划分的结果跟开关状态完全对不上。这个升级版的统一模型,就是把重构和孤岛划分放进同一个混合整数优化框架里同时决策,在Matlab里用YALMIP建模、Gurobi求解,我实测下来,33节点系统单时段恢复方案十几秒就能出一个可行解,多时段模型也能在可接受时间内收敛。这篇文章把模型原理、代码架构和踩过的坑都摊开讲,适合正在做配网恢复算法设计、或者想从分离求解切换到统一模型的同学参考。
1. 统一模型的设计逻辑与升级方向
1.1 重构和孤岛划分不统一,恢复方案就会“打架”
主动配电网发生故障后,恢复供电主要有两个手段:一个是重构,通过调整分段开关和联络开关改变网络拓扑,绕开故障支路;另一个是孤岛划分,让分布式电源在脱离主网的情况下带起一部分负荷,形成一个个独立的供电岛。
问题在于,这俩决策是强耦合的。开关怎么闭合,直接决定了哪些区域还有可能跟变电站连通、哪些区域只能靠DG带;反过来,你期不期望某个区域形成孤岛,也会影响开关动作方案。分开建模的时候,重构模型默认所有负荷都能由主网供应,孤岛划分模型又假设孤岛边界是预先固定的,两套结果经常互相矛盾。我之前做过一个对比实验:同样一个故障场景,先重构再固定孤岛边界,恢复的负荷量比分离开算平均少了10%以上,而且开关倒切次数更多。这不是算法不够好,是模型本身在结构上就漏掉了决策之间的关联。
统一模型的核心思路,是把“开关状态”和“孤岛电源选择”放到同一个决策空间里,共享同一套潮流约束和连通性约束。目标还是最小化失负荷量,但开关动作和DG出力是同时求出来的。这样得到的恢复方案才是工程上可以直接执行的:哪些开关闭合、哪些DG带哪些负荷,是一张完整的操作票。
1.2 目标函数与决策变量怎么定
统一模型的决策变量分三类:支路开关状态、DG出力、实际供应的负荷量。开关状态是二进制变量,决定网络拓扑;DG出力是连续变量,决定孤岛内的功率平衡;实际供应负荷量可以小于原始负荷,表示切除了一部分非关键负荷。
目标函数一般写成多目标加权的形式:
min sum_t [ sum_i alpha_i * (P_load_i(t) - pl_i(t)) + beta * Ploss(t) ] + gamma * sum_t sum_k |z_k(t) - z_k(t-1)|其中alpha_i是节点负荷的权重系数,这个很重要。一级负荷、二级负荷、三级负荷的权重不能随便取,我试过用[10, 1, 0.1]这种量级的分级权重,效果比较符合电网调度的实际逻辑。beta是网损惩罚系数,网损在恢复问题里一般不是主目标,所以权重相对小一些。gamma是开关动作惩罚,多时段模型里不加上这一项,模型会让开关来回倒切,实际现场根本不允许。
约束条件里最有挑战性的是网络辐射状约束和孤岛连通性约束。配网正常运行时要求辐射状拓扑,故障恢复也要求恢复后的网络不成环。传统做法是枚举所有可能的辐射状拓扑,但节点一多就组合爆炸。统一模型里一般用“单商品流”约束或者生成树约束来隐式保证连通性,这个后面展开讲。
1.3 升级版本到底改了什么
这个升级版本对比旧版,最核心的改动有三个维度:
第一,从单时段变成了多时段。故障恢复不是瞬间完成的,DG出力会变化,负荷也在变。升级版模型按时间断面展开,每个断面都满足潮流约束,同时通过储能SOC状态方程把时段耦合起来。旧版一个静态断面算完就结束,根本不考虑DG中午出力大、傍晚出力小这些时序特性,恢复方案在现场执行时经常出现“算出来能带负荷,实际带不动”的情况。
第二,孤岛边界从固定变成动态。旧版孤岛划分把孤岛范围看成固定集合,升级版模型里孤岛边界是由开关状态和DG出力共同决定的,每个时段可以微调。比如同一个DG,午间出力足的时候可以带大一点的孤岛,到了傍晚出力下降,自动切掉部分三级负荷。这种动态孤岛才是主动配电网最需要的。
第三,建模和求解层面做了工程化处理。变量按矩阵批量声明,约束用向量化方式组装,减少了YALMIP内部的开销;求解器配置、初始解构造、Big-M取值这些都调过一轮。升级版代码里我把数据、模型、求解器重构成几个独立模块,换算例时只需要改数据文件,不用动核心建模代码。
2. 核心建模原理与关键公式拆解
2.1 潮流约束:DistFlow方程与二阶锥松弛
故障恢复的决策模型没法用牛顿-拉夫逊法那种交流潮流迭代求解,因为拓扑是变量,雅可比矩阵根本构建不起来。解决办法是用DistFlow支路潮流方程,把潮流关系表达成一组代数约束,放进数学规划框架里。
DistFlow的基本形式是:
P_ij - r_ij * l_ij - sum(P_jk) + p_load_j - p_dg_j = 0 Q_ij - x_ij * l_ij - sum(Q_jk) + q_load_j - q_dg_j = 0 v_j = v_i - 2 * (r_ij * P_ij + x_ij * Q_ij) + (r_ij^2 + x_ij^2) * l_ij其中P_ij、Q_ij是支路首端有功和无功,l_ij是支路电流的平方,v_i是节点电压的平方。这套方程本身是非线性的,直接扔给求解器会很痛苦。实际建模时做一个标准处理:把第三个方程里的等号关系保留,但把电流平方和电压平方的乘积关系用旋转锥约束替换:
|| [2*P_ij, 2*Q_ij, l_ij - v_i] ||_2 <= l_ij + v_i这就是二阶锥松弛。它把原本等式约束的非线性项松弛成一个凸锥,整个问题从非凸的MINLP变成混合整数二阶锥规划(MISOCP)。对辐射状配电网,只要网络本身是连通的,这个松弛在绝大多数实际场景下是紧的,也就是说松弛解跟原问题最优解一致。我做了大量测试,IEEE 33节点系统各种故障场景下没有出现过松弛不紧导致结果离谱的情况,但如果某些支路重载或者电压偏低,还是建议加一个小的惩罚项检查一下对偶间隙。
2.2 开关动作与系统连通性的二进制建模
支路开关状态用二进制变量z_ij表示,1是闭合,0是断开。支路功率只要开关断开就必须为0,这个用Big-M约束实现:
-P_max * z_ij <= P_ij <= P_max * z_ij -Q_max * z_ij <= Q_ij <= Q_max * z_ijBig-M的取值非常影响数值稳定性。我踩过坑:取太大会让求解器误判可行性,取太小又会把本来可行的方案切掉。我的经验是P_max取支路容量上限的1.1倍,电压偏移的Big-M取0.1 p.u.量级,不要为了追求“通用”而设一个特别大的值。
连通性约束是最容易写错的地方。孤岛划分的核心就是保证每个通电节点要么由变电站供电,要么由某个投运的DG供电,而且路径不能环路。我用的是“虚拟功率流”方法:给每个节点注入一个单位的虚拟电流需求,开关闭合的支路才能传输虚拟电流,变电站节点和投入运行的DG节点作为虚拟电流源。约束可以写成:
sum_{k in delta_in(i)} f_k - sum_{k in delta_out(i)} f_k = u_i其中u_i是节点是否恢复供电的二进制变量,f_k是支路虚拟功率,只有z_k=1时f_k才能非零。这个方法的好处是把“连通性”这种图论性质转化成了线性约束,求解器能直接处理。另一个做法是用生成树约束,限制闭合支路数正好等于通电节点数减孤岛数,但虚拟流方法实现起来更直观,调试也更方便。
2.3 孤岛功率平衡和DG出力约束
孤岛运行的核心约束是岛内的有功、无功必须平衡。这个约束不需要显式区分“哪个DG带哪个负荷”,只要通过支路潮流约束自然满足即可。因为每个孤岛都是连通的子图,孤岛内所有节点的注入功率加起来一定为0,这个功率平衡就会自动约束DG容量必须大于岛内负荷。
DG出力本身有上下限约束:
P_dg_min <= P_dg <= P_dg_max Q_dg_min <= Q_dg <= Q_dg_max如果DG是光伏,还要考虑功率因数或者逆变器容量约束。我经常在升级版模型里把DG按PQ节点建模,也就是给定有功出力和功率因数范围,让无功在允许范围内自由调节,这样孤岛的电压质量会好很多。
还有一个活生生的教训:DG的容量约束写成P_dg <= P_dg_capacity * s_dg,这个s_dg是DG是否投入孤岛运行的二进制变量。如果不加这个变量,模型会让一个关着的DG也出力,形成“僵尸DG”。这个问题我刚写模型时犯过一次,查了半天才发现是忘了把DG状态和出力耦合。
2.4 多时段协同和储能SOC建模
升级版最核心的改动是多时段的协同优化。多时段模型把时间轴离散成T个断面,每个断面都有一组潮流约束和开关状态变量,储能SOC方程把相邻时段串起来:
SOC(t+1) = SOC(t) + eta_c * P_ch(t) * dt / E_cap - P_dis(t) * dt / (eta_d * E_cap) SOC_min <= SOC(t) <= SOC_max P_ch(t) * P_dis(t) = 0 % 线性化后变成同一时段只能充或放充放电互斥约束是典型的非线性且非凸约束。实际建模用二进制变量u_ch和u_dis,加约束:
P_ch(t) <= M * u_ch(t) P_dis(t) <= M * u_dis(t) u_ch(t) + u_dis(t) <= 1这样就把互斥关系线性化了。储能SOC初值要按实际运行状态给,不能默认满电——我见过不少论文直接设SOC(1)=1,算出来的恢复方案在实际中根本不可行。升级版代码里我按故障发生前的调度结果给SOC初值,或者直接作为输入参数暴露出来,让使用者自己填。
多时段模型的开关动作次数也要控制。如果不加限制,求解器会给出一种方案:每个时段开关状态都不一样,频繁倒切,现场没法执行。我一般在目标函数里加开关动作惩罚项,再额外加总动作次数上限约束,比如一天内最多动作5次,这样结果工程上更可接受。
3. Matlab代码框架与实操流程
3.1 算例准备:先用IEEE 33节点把流程跑通
做故障恢复研究,IEEE 33节点系统是绕不开的经典算例。这个系统是配网恢复领域的“Hello World”:33个节点、37条支路,其中5条是常开的联络开关,基准电压12.66 kV,总负荷大约3.7 MW加2.3 Mvar。系统不大,但辐射状拓扑、联络开关、负荷分区这些特征全都有,验证统一模型再合适不过。
数据准备这一步最容易出乱子。我建议直接从MATPOWER里读case33bw,或者用网上流传的经典数据脚本,手搓数据容易出错。升级版代码里我把支路数据整理成一个表格,每行是[首端节点, 末端节点, 电阻, 电抗, 初始状态, 是否联络开关],节点负荷单独一张表。DG接在18、22、25、33这几个节点,容量在300到600 kW之间,故障后能形成至少两个有效孤岛。
故障场景可以这样设置:假设节点23到24之间的支路发生三相短路故障,保护动作切除故障支路,这时候节点24到33这一大片区域失电。算例验证时要把故障支路从候选闭合支路集合里剔除,否则模型会自作主张把故障支路“恢复”掉,得出一个完全错误的方案。
3.2 YALMIP建模:变量定义、约束组装与目标函数
Matlab里我推荐用YALMIP作为建模层。YALMIP的语法和数学公式一一对应,排查约束错误比直接调Gurobi的C API或Matlab接口方便得多。下面是建模的关键片段:
% 变量定义 z = binvar(nbranch, T, 'full'); % 支路通断状态 P = sdpvar(nbranch, T, 'full'); % 支路有功 Q = sdpvar(nbranch, T, 'full'); % 支路无功 l = sdpvar(nbranch, T, 'full'); % 支路电流平方 v = sdpvar(nbus, T, 'full'); % 节点电压平方 pg = sdpvar(ndg, T, 'full'); % DG有功出力 qg = sdpvar(ndg, T, 'full'); % DG无功出力 pl = sdpvar(nbus, T, 'full'); % 实际供货负荷约束组装要用循环处理每条支路和每个时段:
Constraints = []; for t = 1:T for k = 1:nbranch i = br_from(k); j = br_to(k); % DistFlow 等式约束 Constraints = [Constraints, ... P(k,t) - R(k)*l(k,t) - sum(P(out_br{j}, t)) + pl(j,t) - pg_at_bus(j,t) == 0]; Constraints = [Constraints, ... Q(k,t) - X(k)*l(k,t) - sum(Q(out_br{j}, t)) + ql(j,t) - qg_at_bus(j,t) == 0]; % 二阶锥松弛 Constraints = [Constraints, ... norm([2*P(k,t); 2*Q(k,t); l(k,t) - v(i,t)]) <= l(k,t) + v(i,t)]; % 开关状态耦合支路功率 Constraints = [Constraints, ... -M_P(k) * z(k,t) <= P(k,t) <= M_P(k) * z(k,t)]; Constraints = [Constraints, ... -M_Q(k) * z(k,t) <= Q(k,t) <= M_Q(k) * z(k,t)]; end end这段代码里的out_br{j}表示所有以节点j为首端的支路集合,需要提前构建好。pg_at_bus(j,t)表示注入到节点j的DG有功,也要做一个节点到DG的映射。目标函数部分:
Objective = 0; for t = 1:T Objective = Objective + sum(alpha .* (P_load - pl(:,t))); Objective = Objective + beta * sum(R .* l(:,t)); end for t = 2:T Objective = Objective + gamma * sum(sum(abs(z(:,t) - z(:,t-1)))); endabs用在二进制变量上没问题,YALMIP会把绝对值线性化。
3.3 求解器配置:Gurobi的安装与调用
YALMIP只是建模层,真正求解MISOCP需要底层求解器。我用的是Gurobi,学术版免费申请,经常做优化研究的人值得备一个。CPLEX和MOSEK也可以,但Gurobi对二阶锥的支持更顺手,数值稳定性也更好。
调用方式:
ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'debug', 1); sol = optimize(Constraints, Objective, ops); if sol.problem == 0 z_opt = value(z); pl_opt = value(pl); disp('求解成功'); else disp(sol.info); end第一次跑的时候务必把debug设为1。YALMIP内置了约束诊断功能,一旦模型里有矛盾约束会直接告诉你哪行出了问题,省掉大量手查的时间。我见过太多人上来就ops全默认,模型不可行了也不知道去哪定位,折腾一晚上不如一个debug开关。
性能方面,升级版模型在IEEE 33节点、24时段、带4个DG和1个储能的规模下,Gurobi 9.x以上版本大概一二十秒收敛到1%的MIP gap。如果跑半小时还没收敛,大概率是模型里有冗余约束或者Big-M取值问题,不是求解器不行。
3.4 结果输出:恢复方案表、电压分布和孤岛划分可视化
算完之后的输出环节很重要,不然一堆二进制数值根本看不出方案好不好。我的习惯是生成三样东西:
第一,开关动作方案表。把每个时段的开关状态变化整理成表格,哪些开关从断开变成闭合、哪些从闭合变成断开,一目了然。这个表本质上就是现场操作票。
第二,电压分布图。把每个恢复节点的电压幅值画出来,重点看有没有节点电压越限。恢复方案不仅要求能满足负荷,还得电压合格。我通常把电压上限设成1.05 p.u.,下限0.95 p.u.,如果是二级以上重要负荷,下限可以放宽到0.90 p.u.,但要在方案说明里写清楚。
第三,孤岛划分示意图。用Matlab画网络拓扑图,恢复区域用一种颜色,失电区域用另一种颜色,孤岛边界高亮出来。我用过plot(G, 'Layout', 'force')加highlight函数,效果足够直观,省得再导入到第三方工具里。这张图拿去汇报或者写论文都是直接能用的。
4. 常见问题与排查技巧实录
4.1 求解报错或返回不可行解
升级版模型跑出来不可行,90%的情况不是求解器的问题,而是约束写错了。我的排查顺序是固定的:先把孤岛连通性约束停掉,只跑重构部分,看是否可行;然后单独打开孤岛约束,关掉多时段,看单时段是否可行;最后才放开多时段。这样做能快速定位是哪个模块引入了不可行。
另一个常见原因是DG容量和孤岛负荷不匹配。比如DG只有300 kW,孤岛里却挂了800 kW的负荷,功率平衡约束永远满足不了。不要想着靠目标函数“惩罚”,不可行问题直接让求解器崩溃。解决办法是检查DG容量曲线,或者把负荷按可切性分组,允许切除部分三级负荷。
电压上下限约束也有可能造成不可行。33节点系统里末端节点电压本来就偏低,故障恢复后线路路径变长,电压下降更严重。我一般先按0.90到1.05 p.u.设置一个较宽的范围跑通模型,再逐渐收紧到合格范围,看哪些节点需要加无功补偿或者调整DG出力。
4.2 孤岛约束失效,恢复区域出现环路
环路的典型特征是闭合支路数大于通电节点数减孤岛数。虚拟流方法失效的常见原因,是虚拟流的Big-M取值不对。约束里f_k <= M * z_k,M取太小会把合法路径也截断,取太大又让求解器数值发虚。我的经验值是M等于节点总数的2倍,对33节点系统取66,实测很稳。
还有一种隐蔽问题:虚拟流模型把变电站和DG都设为电源,但如果某个DG虽然投入但没有实际出力,模型会给它设置虚拟流为0,那么这个孤岛实际上没有电源,节点却可能因为开关闭合而“通电”。要避免这个问题,必须保证DG出力下限约束和虚拟流电源状态绑定,也就是s_dg=1时DG有功出力必须大于某个下限,才允许这个节点作为虚拟流源。
4.3 多时段模型越跑越慢,怎么调优
多时段模型规模是单时段的T倍,但如果运行时间不是线性增长而是爆炸式增长,说明模型里有非凸约束或者不必要的二进制变量。我见过最典型的问题是储能模型里写了P_ch * P_dis == 0这种非线性等式,YALMIP会把它当MINLP处理,求解时间直接翻几十倍。要严格按照上面写的二进制变量互斥方法建模,始终保持MISOCP结构。
另外可以加初始解。Gurobi支持传入MIP start,我把单时段模型求出来的解按时间复制一份作为多时段模型的初始解,收敛速度有显著提升。具体做法是先用assign给YALMIP变量赋初值,YALMIP会自动传给求解器。
4.4 MATLAB版本、YALMIP和求解器的环境坑
环境问题虽然不算模型算法本身,但真的会卡人。我在MATLAB R2023b上遇到过一个很典型的问题:YALMIP升级到最新版后,sdpsettings里的某些Gurobi选项被更新了,旧代码直接报错找不到选项名。解决方案是先运行yalmiptest检查兼容性,再逐个试配置项。
中文注释乱码问题也值得说一下。R2023b以后MATLAB默认编码改成了UTF-8,但很多老脚本是用GBK存的,打开后中文注释全是乱码。解决方法是把文件另存为UTF-8编码,或者改MATLAB编辑器的默认编码设置。这个跟模型本身无关,但真要遇上了很窝火。
4.5 排障速查表
| 现象 | 可能原因 | 处理建议 |
|---|---|---|
| 模型不可行 | 孤岛内负荷超过DG容量 | 允许切除三级负荷或调整DG位置 |
| 模型不可行 | 电压约束过严 | 先放宽到0.90~1.05 p.u.再逐步收紧 |
| 出现环路 | 连通性约束Big-M不合适 | 虚拟流Big-M取节点数的2倍 |
| 求解极慢 | 存在P_ch * P_dis == 0非线性项 | 改用二进制变量表示充放电互斥 |
| 求解极慢 | Big-M取值太大 | M取支路容量上限的1.1倍,电压M取0.1 p.u. |
| 结果有“僵尸DG” | DG状态变量没与出力耦合 | 加P_dg <= P_dg_capacity * s_dg约束 |
| 中文注释乱码 | 文件编码GBK与MATLAB默认UTF-8冲突 | 另存为UTF-8编码再重新打开 |
| YALMIP报错 | 求解器版本和YALMIP不兼容 | 运行yalmiptest检查各项支持 |
最后再分享一个我实测下来的小技巧:升级版模型里,把目标函数中的三级负荷权重设成非常小但不为0,会显著降低搜索难度。原因是模型需要判断哪些负荷优先切除,如果权重一样,求解器要在对称解里反复试探,浪费大量时间。拉开权重差距相当于给搜索树加了方向,收敛速度提升有时候非常明显。做这行的都清楚,算法研究最怕的不是模型复杂,而是模型复杂到连最优解都找不到;统一模型本身已经把问题规模压下来了,剩下的就是这些建模细节的打磨。