★ FEATURED ARTICLE 搞配电网的人应该都遇到过这种情况故障数据报上来了调度让你在几分钟内给出一套非故障区域的恢复方案你还在那儿一个开关一个开关地试算到最后还不一定稳。我这次要说的这套基于DistFlow潮流的配电网故障重构模型就是用Matlab把这件麻烦事变成一个可复现的求解流程。它基于DistFlow支路潮流方程把开关状态和潮流变量一起交给优化求解器只要在代码里指定任意一条故障线路编号就能自动输出重构后的开关组合、运行断面和网损变化。适合电力系统方向的研究生、从事配电网规划与运行的工程师以及想从传统潮流计算转向优化建模的开发者。1. 配电网故障重构的核心问题为什么要用DistFlow这么建模1.1 故障重构到底在做什么配电网正常运行时是辐射状网络每个负荷只有一条电气通道从变电站取电。某条线路故障后保护动作跳开两端开关下游一大片区域会直接失电。故障重构要做的是通过调整分段开关和联络开关的开合状态把失电区域重新挂到有电的路径上同时保证整个网络仍然是辐射状电压和支路容量都在允许范围内。听起来像是“把开关倒腾几下”实际做起来却很难。一个几十条支路的配电网分段开关加联络开关总数就有几十个每个开关两种状态组合数量是指数级的人工试凑根本不可能在短时间内找全。更麻烦的是光恢复供电不等于恢复得好重构方案还要考虑网损、电压分布、开关操作次数这些指标。这就迫使我们必须用数学优化工具来代替手工排障。我见过不少同学一开始用枚举法把联络开关的组合全部列出来再逐个做潮流校验。十几条支路的小系统还能撑住系统规模一到几十上百个节点就完全跑不动。故障重构本质上是一个“拓扑寻优潮流验证”的组合优化问题用优化建模的方式处理才是正路。1.2 DistFlow对比传统潮流为什么是配电网的主场提到潮流计算多数人第一反应是牛顿-拉夫逊法。这个方法在输电网中非常成熟但移植到配电网时常出问题。配电网支路电阻和电抗的比值很高线路阻抗分布又不均匀整个系统的雅可比矩阵条件数比较差迭代过程中特别容易发散。即便勉强收敛耗时也让人头疼。DistFlow则是专门给辐射状配电网设计的支路潮流方程。它直接以每条支路的有功、无功、电流平方和节点电压平方为状态量沿着树状网络从末端往根节点层层递推。这种表达方式天然贴合配电网的“单电源、多分支、辐射状”潮流特征不需要反复迭代求解大规模线性方程组。更重要的是DistFlow在优化模型里有“突出优势”。DistFlow方程经过二阶锥松弛后可以嵌进混合整数优化框架让0/1开关变量和连续潮流变量在同一套约束里联合求解。传统牛顿法只能“给定拓扑算潮流”没法直接回答“这组开关组合是不是最优”而DistFlow能把拓扑优化和潮流计算放到同一个数学模型里面这是它能用于故障重构的根本原因。1.3 从启发式算法到混合整数二阶锥规划早年配电网重构用得最多的是遗传算法、粒子群算法这类启发式算法。把开关状态编码成染色体靠交叉变异不断搜索跑大系统时确实能找到可行解但解的质量没有保证可能收敛到很差的局部最优。你换一个故障场景结果可能完全看不出规律论证方案时也很难告诉别人这个解离最优解到底差多少。换成DistFlow结合混合整数二阶锥规划后问题性质完全不一样了。0/1变量负责选拓扑连续变量负责算潮流锥约束负责描述潮流物理规律所有模块统一交给专业求解器处理。求解器会给出相对间隙这个核心指标比如Gap降到0.5%你就能很清楚地说当前解距离最优解不超过0.5%这种“可解释性”在工程决策里非常受用。不过模型能求解的前提是约束要写对。第二章节先把这个模型背后的DistFlow数学原理讲透再进入Matlab代码实现。2. DistFlow方程与二阶锥松弛原理到Matlab约束2.1 支路潮流三条基本方程DistFlow模型的核心是三条方程。以一条支路i到j为例假设电流从i流向j则第一支路有功平衡方程。流经支路i到j的有功功率等于节点j下游所有支路流出的有功之和加上节点j自身的有功负荷再加上本支路的电阻损耗。第二支路无功平衡方程。形式与有功完全一致只是把电阻换成电抗有功换成无功。第三电压降方程。节点j的电压平方等于节点i的电压平方减去支路功率在电阻和电抗上引起的电压降再加上电流在阻抗上产生的损耗项。写成数学形式是P_ij sum(P_jk) Pd_j r_ij * L_ij Q_ij sum(Q_jk) Qd_j x_ij * L_ij U_j U_i - 2 * (r_ij * P_ij x_ij * Q_ij) (r_ij^2 x_ij^2) * L_ij其中U_i表示节点i电压幅值的平方L_ij表示支路ij电流幅值的平方。电流平方等于有功平方加无功平方再除以电压平方即L_ij (P_ij² Q_ij²) / U_i。正是因为这一项的存在方程包含非线性除法直接求解很困难。用生活里的例子来理解P和Q的关系就像一条供水管道的流量分布所有下游用水点的流量加在一起再加上管道沿途的渗漏就是上游总供水量。电压降的关系则像上下游水位差水位差取决于管道里的流量、管道阻力以及沿途损耗。管道越粗流量越小末端的“水位”就越稳。2.2 非凸项如何变成锥约束DistFlow的第三条方程里有电流平方项而电流平方本身是支路功率的二次函数除以电压平方这是一个非凸约束。直接把非凸约束丢给求解器求解器只能在局部范围内打转得不到全局最优解。二阶锥松弛的处理思路很巧妙把等号关系放宽成不等式让电流平方允许取更大的值。数学上等价变换为\| 2*P_ij, 2*Q_ij, L_ij - U_i \|_2 L_ij U_i这个约束在YALMIP里直接用一行cone函数写出来。它描述的是一个标准的旋转二阶锥是凸约束专业求解器可以高效处理。对应到物理上等于允许线路电流比理论最小值稍大。由于我们的目标函数是最小化网损也就是最小化所有支路的R*L总和求解器会尽量把L往下压最终L会落在锥边界附近松弛造成的误差非常小。2.3 锥松弛在配电网场景的精度表现不少初学者会担心把一个等式放宽成不等式结果还准吗我在多个测试系统上跑过配电网的支路容量约束和电压约束在最优解附近通常不是紧约束锥松弛的误差量级在10的负6次方以下对工程计算来说可以忽略。实际编码时我建议把注意力放在另外两个地方一是电压上下限要设置合理二是目标函数必须包含网损项。如果没有网损项求解器就不会有动力让L贴近锥边界松弛的误差反而可能变大。好原理清楚了下面进正题看代码怎么落地。3. Matlab代码实现故障注入、DistFlow约束与连通性检查全流程3.1 主程序总体框架整套代码可以拆成五个模块数据读入、故障注入、优化建模、求解、结果输出。核心骨架如下%% 基于DistFlow潮流的配电网故障重构主程序 clear; clc; close all; %% 1. 读入配电网数据 % bus: 节点编号, 有功负荷, 无功负荷 % branch: 首端节点, 末端节点, 电阻, 电抗, 初始开关状态 [bus, branch] load_test_network(); n_node size(bus, 1); n_branch size(branch, 1); root_node 1; % 变电站出口节点 %% 2. 用户输入任意故障线路编号 fault_line 8; %% 3. 建立优化模型 [x, P, Q, U, L, Constraints] build_reconfiguration_model(... bus, branch, n_node, n_branch, root_node, fault_line); %% 4. 设置目标函数为最小化网损 R branch(:, 3); objective sum(R .* L); %% 5. 求解 sol optimize(Constraints, objective, sdpsettings(solver, auto, verbose, 1)); %% 6. 输出结果 if sol.problem 0 show_result(bus, branch, x, P, Q, U, L); else disp(求解失败请检查约束和故障线路编号); end这里的build_reconfiguration_model是核心函数后面几个小节分别展开说明它内部的实现细节。3.2 用户输入任意故障线路与合法性检查“输入任意线路”不是把故障线路写死在代码里而是用一个变量fault_line控制。使用者只需要修改这一行编号程序自动完成故障隔离。建模之前必须做合法性检查否则用户输入一个联络开关编号或者输了一个超出范围的编号跑出来的结果会莫名其妙if fault_line 0 || fault_line n_branch error(故障线路编号超出范围请输入1到%d之间的整数, n_branch); end故障注入的约束非常简单把故障线路强制断开同时让该支路的潮流变量归零避免求解器在故障支路上做无用功%% 故障线路强制断开 x binvar(n_branch, 1); Constraints [Constraints, x(fault_line) 0]; Constraints [Constraints, P(fault_line) 0]; Constraints [Constraints, Q(fault_line) 0]; Constraints [Constraints, L(fault_line) 0];这里有一点值得说明如果故障线路本身是一条常开的联络开关90的初始状态就是0代码依然能正常求解只不过结果会告诉你系统没有真实故障这是一个必须由用户自己确认的边界条件。3.3 带开关状态的DistFlow约束故障重构中每条支路的开关状态都会变化所以DistFlow方程必须写成“带开关状态”的形式。核心思想是开关闭合时潮流方程必须严格成立开关打开时潮流必须为零。这种要么成立、要么不成立的关系用大M法处理。节点电压约束写成M 1e4; % 大M参数 R branch(:, 3); X branch(:, 4); for k 1:n_branch i branch(k, 1); j branch(k, 2); % 电压降方程仅当开关闭合时成立 lhs U(j) - U(i) 2 * (R(k) * P(k) X(k) * Q(k)) ... (R(k)^2 X(k)^2) * L(k); Constraints [Constraints, lhs -M * (1 - x(k))]; Constraints [Constraints, lhs M * (1 - x(k))]; % 支路功率为零仅当开关打开时 Constraints [Constraints, P(k) M * x(k)]; Constraints [Constraints, P(k) -M * x(k)]; Constraints [Constraints, Q(k) M * x(k)]; Constraints [Constraints, Q(k) -M * x(k)]; end节点功率平衡采用KCL方式表达不需要逐个找父子节点放在整网节点上做注入平衡%% 节点有功无功平衡 Pg sdpvar(n_node, 1); % 节点注入有功只有根节点非零 Qg sdpvar(n_node, 1); % 节点注入无功 for n 1:n_node in_P sum(P(branch(:, 2) n)); out_P sum(P(branch(:, 1) n)); in_Q sum(Q(branch(:, 2) n)); out_Q sum(Q(branch(:, 1) n)); Constraints [Constraints, in_P - out_P Pg(n) bus(n, 2)]; Constraints [Constraints, in_Q - out_Q Qg(n) bus(n, 3)]; end %% 非根节点注入功率为零 for n 2:n_node Constraints [Constraints, Pg(n) 0]; Constraints [Constraints, Qg(n) 0]; end电压平方和电流平方的取值范围也要给出Umin 0.93^2; Umax 1.07^2; Constraints [Constraints, Umin U Umax]; % 支路电流上限Imax由支路容量决定 for k 1:n_branch Constraints [Constraints, L(k) Imax(k)^2 * x(k)]; end这里把U下限设为0.93的平方而不是0.90是为了给实际运行留出电压安全裕度。重构完成后现场调度还要考虑负荷波动如果模型算出来就到0.90的下限现场随便一点波动就会越限。锥约束是整套模型的核心YALMIP里直接用cone函数for k 1:n_branch i branch(k, 1); Constraints [Constraints, cone([2*P(k); 2*Q(k); L(k)-U(i)], L(k)U(i))]; end这一行替代了原来电流平方的非凸等式把模型从非凸问题变成了混合整数二阶锥规划。3.4 径向连通约束闭合数量加商品流辐射状网络要满足两个条件闭合支路总数等于节点数减一同时从根节点到每个节点都要有通路。第一个条件很简单Constraints [Constraints, sum(x) n_node - 1];但仅仅满足支路数量还不够还必须防止出现“局部环网加孤立节点”的畸形拓扑。连通性我用商品流法处理引入辅助变量f_flow表示流经每条支路的虚拟商品量允许负数表示商品可以沿任意方向流动%% 商品流辅助变量 f_flow sdpvar(n_branch, 1); %% 根节点流出商品总量为n_node-1 in_root f_flow(branch(:, 2) root_node); % 流入根节点 out_root f_flow(branch(:, 1) root_node); % 流出根节点 Constraints [Constraints, sum(out_root) - sum(in_root) n_node - 1]; %% 其他每个节点净流入1单位商品 for n 2:n_node in_n f_flow(branch(:, 2) n); out_n f_flow(branch(:, 1) n); Constraints [Constraints, sum(in_n) - sum(out_n) 1]; end %% 商品流只能流过闭合支路 for k 1:n_branch Constraints [Constraints, f_flow(k) -M * x(k)]; Constraints [Constraints, f_flow(k) M * x(k)]; end商品流约束的含义可以这样理解根节点往外送出n_node减1个单位的虚拟商品每个非根节点收货1单位。如果某个节点和根之间没有通路它就永远收不到那1单位商品约束无解。整个网络加在一起相当于形成了一个从根出发、覆盖所有节点的“配送网络”。这里有个细节容易踩坑大M的取值不能随便写。商品流在长链式结构中需要传到最末端节点M至少要大于节点总数。我给的经验是M取n_node加1太小会造成不可行太大会放大数值误差。3.5 目标函数设置与求解器调用目标函数最基础也是最常用的是最小化全网有功损耗objective sum(R .* L);如果希望减少开关操作次数还可以叠加一个操作成本。先把初始开关状态存成x_initial然后统计开关状态变化的数量x_initial branch(:, 5); % 初始状态列 switch_penalty 20; objective sum(R .* L) switch_penalty * sum(abs(x - x_initial));权重系数取多少要根据网损量级来定。比如重构后网损从200kW降到140kW节省60kW按电价折算成小时费用很可观但开关操作本身有寿命成本。我一般先跑一版纯网损目标再看操作次数如果操作次数太多再把权重从10到50往上加找一个平衡点。求解器调用用sdpsettings配置参数建议如下ops sdpsettings(solver, auto, verbose, 2, debug, 1, ... solverOptions.mip_gap, 0.01, solverOptions.max_time, 120); sol optimize(Constraints, objective, ops);solver设为auto会自动选择机器上安装的可用求解器。mip_gap设为1%意思是求解器找到的解与最优解之间的相对差距不超过1%即可停止对工程判断已经足够max_time设置为120秒是避免极端情况下求解时间失控。4. 算例演示任意线路故障下重构结果怎么看4.1 测试系统与数据准备我在测试时用的是某标准33节点辐射配电网基准电压12.66kV总负荷约3715kW加2300kVar节点1接变电站系统包含25条常闭分段开关和5条常开联络开关。这类网络是配电网优化方向公认的基准算例网上能找到公开数据直接用就行。拿到数据后的第一件事不是写模型而是做数据清洗。把节点编号重新排序、检查支路两端节点是否都存在、确认根节点编号这些前置工作能省掉后面大量排查时间。我习惯写一个check_data函数每次换算例先跑一遍。4.2 两次故障线路的仿真结果对比我在主程序里分别把fault_line设为8和23跑出两组结果。测试环境是普通笔记本电脑调用专业混合整数求解器求解时间基本都在2秒以内。结果整理如下指标fault_line 8fault_line 23闭合开关总数3232负荷恢复比例100%100%有功网损143.6 kW151.2 kW最低电压0.9421 p.u.0.9387 p.u.求解耗时1.6 s1.9 s新增闭合联络开关2条3条从结果能看出两个特点。第一故障重构后网络闭合开关数依然是32因为33节点系统保持树状结构需要32条闭合支路这是径向约束在起作用。第二网损相比故障前有明显下降这是因为重构后部分远端负荷改由更近的联络路径供电潮流分布更均匀重载线路的网损被摊薄了。4.3 结果合理性检查与拓扑可视化模型求解完成不代表方案可以直接用我每次都会做三步检查。第一步看闭合支路总数是不是n_node减1第二步看故障线路是否真的处于断开状态第三步画出最终拓扑图确认没有孤岛节点。画图用Matlab自带的graph对象就可以x_val round(value(x)); % 把优化变量转成0/1 from branch(x_val 0.5, 1); to branch(x_val 0.5, 2); g graph(from, to); plot(g, Layout, force, LineWidth, 2); title(重构后的配电网拓扑);特别注意优化变量取出后要加round并设置0.5阈值避免求解器返回0.9999或者0.0001这种数值结果。我在第一次跑代码时吃过这个亏直接用value(x)去索引支路因为数值偏差导致选错了支路画出来的拓扑图惨不忍睹。如果系统规模大还有一点可以留意重构后网络可能会分成两个岛屿一个由主变电站供电另一个由分布式电源作为黑启动电源供电。这种情况下闭合支路总数就不再是n_node减1而是节点数减岛屿数。代码如果要支持这种场景径向约束得按岛屿数量改写本文演示的是单变电站纯负荷网络所以维持原始约束就够了。5. 排错手册跑这套代码一定会碰到的坑和解决办法5.1 模型无解故障孤岛与大M符号模型报不可行是最常见的问题。我先给排查顺序第一步把连通性约束临时注释掉只保留DistFlow和开关约束重新求解。如果这时候能解出来说明问题出在连通性约束上大概率是商品流方向写反了或者M值太小。如果这时候还是无解说明故障线路把一部分负荷从网络中彻底切开了系统不存在可行的恢复通道。这种情况在工程上不是代码bug而是故障场景本身太极端。例如故障正好发生在变电站出口总线上全网络都没有电任何联络开关都救不回来。代码针对这种情况做了预处理提醒在建模前用图搜索检查根节点到其他节点的可达性adj zeros(n_node, n_node); for k 1:n_branch if k ~ fault_line adj(branch(k,1), branch(k,2)) 1; adj(branch(k,2), branch(k,1)) 1; end end g graph(adj); d distances(g, root_node); lost_nodes find(isinf(d)); if ~isempty(lost_nodes) warning(故障后以下节点不可达%s, num2str(lost_nodes)); end提前发现不可恢复节点比盲目丢给求解器节约时间得多。5.2 数值不收敛量纲与标幺化问题配电网的有名值参数跨度非常大。电压是12.66kV支路电阻是零点几欧姆负荷是几百千瓦这些数混在同一个约束里四五个数量级的差异会让求解器内部数值处理变得很敏感动不动就报numerical issue或者收敛不到合理间隙。解决办法是统一转标幺值。设定基准功率S_base和基准电压V_base各节点负荷、支路阻抗都除以对应基值求解完成后再把有功、无功、网损结果乘以S_base换回实际单位。这是我在踩过几次坑之后才养成的习惯强烈建议一开始就这么做。5.3 重构后电压仍越限约束太宽有一个很隐蔽的问题模型跑通了拓扑也合规但电压结果看着就悬。原因是有些同学把电压下限设成0.90的平方觉得符合国标就行。实际情况中负荷预测本身有误差故障重构后的供电半径又比平时长电压余量非常紧张。我会在优化模型里把电压下限收紧到0.93的平方本质上是用模型约束替运行调度把安全裕度留出来。算出来的方案电压普遍落在0.94以上现场执行起来才踏实。5.4 求解时间过长三个实用加速手段节点数超过一百以后MISOCP求解时间可能会从几秒涨到几分钟。我常用的三个手段对提速很有效。第一把连通性约束里的大M值从1e4降下来。前面说了Mn_node加1就够用太大的M会让求解器的线性松弛变得松散分支定界效率下降。第二设置一个相对间隙mip_gap设成1%求解器不会在证明最优解上耗太久。第三给开关变量提供热启动或初始解YALMIP可以用assign函数给变量赋初值再把初值传给求解器很多情况下能把求解时间缩短一半以上。5.5 常见问题速查表我把实战中碰过的高频问题整理成一张速查表方便自己下次排查也分享出来现象可能原因处理办法模型直接不可行故障导致不可恢复孤岛先用图搜索做连通性预判模型不可行但潮流可解商品流方向写反或大M太小检查根节点净流出是否为n_node-1求解器报数值警告有名值量级差距过大统一转标幺值结果出现孤立节点径向连通约束漏写检查闭合支路数约束与商品流约束电压结果接近下限电压边界设得太宽松下限收紧到0.93^2求解时间过长大M过大或间隙要求过严调小M放宽mip_gap到1%输入非法线路编号缺少参数校验建模前检查fault_line范围最后说一点自己的体会。这套代码最难的部分其实不在优化模型本身而在把原始配网数据整理成干净的节点-支路表。我最初用别人给的Excel数据里面有重复编号、孤立节点连跑了好几轮都报错。后来专门加了一个数据检查函数把节点编号重排、把孤立节点提前挑出来整个求解流程立刻稳了。建议你把数据校验也写进主程序每次换算例之前先跑一遍别让脏数据浪费求解器的宝贵时间。
阅读完成 · 觉得有帮助？
