简介:面向电力系统自动化与智能电网方向的技术人员及研究人员,这份PDF文档围绕二阶锥松弛(SOCP)在主动配电网故障重构中的应用展开,将非凸潮流方程转化为可高效求解的凸优化形式,并结合改进的禁忌克隆遗传算法搜索最优供电恢复拓扑,配合颜色编码可视化模块直观呈现网络状态,适合需要快速制定故障恢复方案、提升电网可靠性的中高级学习场景。资源包共1个PDF文件,约125KB,内容涵盖核心算法实现、故障模拟、拓扑校验与可视化等完整MATLAB代码示例,并讨论了目标函数设计、约束处理及含分布式电源的主动配电网注意事项。已有62人学习,读者可从中获取从故障检测到最优方案输出的全流程实现思路与可复用脚本,便于对照调试与二次开发。
1. 配电网故障重构:当线路跳闸后,8 秒内给出恢复方案
凌晨两点,配电房警报响了,某条 10kV 线路跳闸。运维值班的人最怕这种时刻——不是怕停电本身,而是怕在几十条备选线路里找不到那个既满足辐射状约束、又能把负荷损失压到最低的恢复方案。传统做法要么靠经验手动倒闸,要么跑遗传算法等上几十秒还未必收敛到全局最优。这份 MATLAB 代码干的事很直接:把非凸的配电网潮流方程用二阶锥松弛(SOCP)转成凸优化问题,故障发生后自动搜索最优拓扑,再用颜色编码把结果画出来。适合正在做配电网自动化、智能电网方向的研究生和工程师,尤其是需要快速验证重构算法、又不想从零搭潮流模型的人。代码基于 IEEE 33 节点系统,主线清晰,改起来不费劲。
2. 架构拆解:从 case33bw 到求解器的四层结构
2.1 整体分层与数据流
这份代码不是一个大脚本从头跑到尾,拆开看是四层结构,每层职责很明确。第一层是数据层,case33bw()负责生成 IEEE 33 节点系统的节点-支路关系矩阵,包含首末端节点编号、电阻、电抗、开关状态这些字段。第二层是建模层,build_socp_model()把物理约束翻译成 YALMIP 能识别的约束表达式,核心工作就是把非线性潮流方程做二阶锥松弛。第三层是求解层,调用 CPLEX 或 GUROBI 解这个凸优化问题,拿到各支路功率和节点电压。第四层是展示层,visualize_topology()把结果画成拓扑图,故障线路标红、联络开关标绿、电压偏低节点变黄。
数据流是这样的:main()先读系统参数,设置故障支路,把故障线路阻抗拉到极大值模拟断开,然后构建 SOCP 模型丢给求解器,最后拿求解结果刷新可视化。整个链路没有隐藏的状态传递,每一层的输入输出都是显式的结构体或矩阵,调试的时候在任意一层打断点都能看清数据长什么样。
我一般看这类代码先找main(),因为它把调用顺序和参数流向全暴露了。这份代码的main()写得很直白:
function main() % 初始化电网参数:节点、支路、负荷 [bus, branch] = case33bw(); % 设置故障线路编号,这里以 12 号支路为例 fault_branch = 12; branch(fault_branch, 4) = 1e6; % 故障阻抗设为极大值,等效断开 % 构建二阶锥松弛模型 model = build_socp_model(bus, branch); % 调用求解器(默认 CPLEX,可切换) result = solve_with_cplex(model); % 可视化:拓扑 + 故障标记 + 重构结果 visualize_topology(branch, fault_branch, result.x); end这段代码里branch矩阵的第 4 列存的是支路阻抗,把故障支路阻抗设成1e6是一种“软断开”手法——比直接删行更省事,因为拓扑矩阵的维度不变,后续索引不用重新映射。result.x是求解器返回的决策变量向量,里面包含各支路有功、无功和节点电压平方值,可视化函数靠它来着色。
2.2 二阶锥松弛的数学转化:为什么 norm 不等式是关键
配电网潮流的原始约束里,支路功率和电压的关系是非线性的,直接扔给求解器会卡在非凸性上。二阶锥松弛的做法是引入辅助变量,把电压平方和电流平方单独拎出来当变量,然后用一个范数不等式把非线性项包进锥里。代码里这段是核心中的核心:
% 支路潮流约束:对每条支路 k 构建 SOC 约束 for k = 1:nbranch i = branch(k, 1); % 支路首端节点 j = branch(k, 2); % 支路末端节点 Pij = sdpvar(1); % 支路有功功率 Qij = sdpvar(1); % 支路无功功率 constraints = [constraints, ... % 节点功率平衡:注入功率 = 流出功率 - 损耗 P(k) == Pij - r(k) * (I(k)^2), ... Q(k) == Qij - x(k) * (I(k)^2), ... % 电压降方程:用电压平方差表示 u(j) == u(i) - 2*(r(k)*Pij + x(k)*Qij) + (r(k)^2 + x(k)^2)*I(k)^2, ... % 二阶锥约束:把 (P,Q,I,u) 的关系压进锥里 norm([2*Pij; 2*Qij; u(i)-I(k)^2], 2) <= u(i)+I(k)^2]; end逐行看。前两行是节点功率平衡,P(k)是节点 k 的净注入有功,Pij是支路首端流出的有功,r(k)*(I(k)^2)是支路电阻上的有功损耗。第三行是电压降方程,注意这里用的是电压平方u而不是电压本身,因为平方之后欧姆定律里的二次项能写成线性组合。最后一行norm([...], 2) <= ...就是二阶锥约束,它等价于(2Pij)^2 + (2Qij)^2 + (u(i)-I(k)^2)^2 <= (u(i)+I(k)^2)^2,展开后正好还原出Pij^2 + Qij^2 <= u(i)*I(k)^2这个原始的非线性关系。
参数上要留意u和I的定义域。u是电压平方,物理上必须非负,代码里通过u >= 0的约束保证;I是电流平方,同样非负。如果求解器报 infeasible,第一件事就是检查这两个变量有没有被意外赋了负值。另外norm的第二个参数2不能省,虽然默认就是 2-范数,但显式写出来在排查时更清楚。
2.3 故障模拟与拓扑约束:软断开比硬删除更稳
故障模拟这块,代码用了两种手法。第一种是前面看到的把阻抗拉到1e6,第二种是直接改状态矩阵:
% 方式一:阻抗软断开(推荐,拓扑维度不变) fault_branch = 15; branch(fault_branch, 4) = 1e8; % 方式二:状态位硬标记(用于可视化区分) branch_status = ones(nbranch, 1); branch_status(fault_branch) = 0;为什么推荐方式一?因为配电网重构的约束里有一堆跟支路索引绑定的表达式,如果直接branch(fault_branch,:) = []删行,后面所有for k = 1:nbranch的循环索引全得重算,联络开关的映射关系也会乱。用大阻抗等效断开,物理意义是“这条支路流过的功率趋近于零”,求解器能自然处理,代码也不用改结构。
拓扑约束里最容易被忽略的是辐射状校验。配电网正常运行时必须是辐射状(没有环),重构方案如果搞出环,保护装置会误动。代码里check_radial()函数干的就是这事:
function feasible = check_radial(status) % 基于图论:闭合支路数 = 节点数 - 1 且连通 n_closed = sum(status == 1); n_node = size(bus, 1); if n_closed ~= n_node - 1 feasible = false; return; end % 用邻接矩阵做连通性检查 adj = build_adjacency(status); feasible = (rank(adj) == n_node - 1); end这个函数在遗传算法的每一代都会调用,用来筛掉不满足辐射状的个体。n_closed ~= n_node - 1是必要条件,rank(adj)检查连通性是充分条件,两个都过才算可行。实测在 33 节点系统里,这个校验单次耗时不到 1 毫秒,比把辐射状约束写进优化模型里让求解器处理要快得多。
3. 跑起来:环境配置、求解器切换与第一次运行
3.1 依赖清单与版本踩坑
这份代码跑起来需要三样东西:MATLAB 本体(R2018b 及以上,因为用到了sdpvar的某些重载)、YALMIP 工具箱(建模层)、一个能解二阶锥的求解器(CPLEX / GUROBI / MOSEK 三选一)。YALMIP 直接去官网下最新版解压到 MATLAB 路径下就行,不用编译。求解器这块,CPLEX 对学术用途有免费授权,配置好之后在 MATLAB 里yalmiptest能检测到就算通了。
版本上有个血泪经验:YALMIP 2020 之前的版本对norm约束的解析有 bug,会把二阶锥约束错误地转成一般的非线性约束,导致求解器报“不支持该约束类型”。我第一次跑的时候就卡在这,yalmiptest显示 CPLEX 可用,但一跑主程序就报错。后来把 YALMIP 换到 2021 之后的版本才解决。所以如果你跑出来报 SOC 约束不支持,先查 YALMIP 版本,别急着怀疑求解器。
环境配好之后,在 MATLAB 命令行里跑一遍yalmiptest,看输出里 CPLEX 那一行是不是found。如果是not found,检查求解器路径有没有加到 MATLAB 的path里,或者savepath有没有执行成功。
3.2 主程序运行与求解器切换
主程序入口就是main(),直接跑就行。默认走 CPLEX,如果你想换 GUROBI 或 MOSEK,改solve_with_cplex()里的 solver 配置:
function result = solve_with_cplex(model) % 求解器配置:可切换 cplex / gurobi / mosek options = sdpsettings('solver', 'cplex', ... 'verbose', 1, ... 'cplex.timelimit', 60); % 求解并检查状态 sol = optimize(model.constraints, model.objective, options); if sol.problem ~= 0 warning('求解失败:%s', sol.info); end result.x = value(model.variables); endsdpsettings里verbose设 1 会打印求解过程,第一次跑建议开着,能看到求解器迭代了多少轮、gap 收敛到多少。cplex.timelimit是超时限制,单位秒,33 节点系统一般 10 秒内就收敛,设 60 是留余量。如果换成 GUROBI,把'cplex.timelimit'改成'gurobi.timelimit',其余不变。
跑完之后命令行会输出类似“第 3 轮迭代,间隙值 0.0021”的日志,这是while gap > 1e-4循环在打印收敛过程。最终 gap 小于1e-4就退出,弹出拓扑图。如果图没弹出来,检查visualize_topology()里有没有figure被注释掉,或者 MATLAB 的图形后端是不是设成了-nodisplay。
3.3 结果解读:红色虚线、绿色节点和电压色阶
跑通之后弹出的图里,红色粗虚线是故障线路,绿色实心点是重构后闭合的联络开关,节点颜色按电压标幺值渐变——低于 0.95 p.u. 的偏黄,正常范围偏蓝。图标题会显示“辐射状拓扑校验通过”或“校验失败”,后者说明重构方案有环或者不连通,得回去查约束。
有个细节值得注意:故障线路在图上画的是虚线而不是直接消失,这是故意的。运维人员看的时候需要知道“哪条线断了”,如果直接不画,反而丢失了故障位置信息。代码里plot(branch(fault_branch,1:2), 'r--', 'LineWidth', 3)这行就是干这个的,'r--'是红色虚线,LineWidth设 3 让它比正常线路粗一圈,视觉上更突出。
电压色阶用的是colormap自定义的红蓝配色,不是 MATLAB 默认的jet。默认jet在中间段颜色跳变太剧烈,看电压分布的时候容易误判。自定义的[0.8 0.2 0.2; 0.2 0.6 1]从红到蓝过渡平滑,低电压区域一眼就能认出来。
4. 核心算法拆解:禁忌克隆遗传算法怎么搜拓扑
4.1 种群初始化与克隆扩增
重构搜索用的是改进的禁忌克隆遗传算法,跟标准遗传算法比,多了克隆扩增和禁忌筛选两步。初始化的时候随机生成一批满足辐射状约束的拓扑个体,每个个体是一个二进制向量,1 表示支路闭合、0 表示断开。克隆扩增是把当前种群里的精英个体复制多份,扩大搜索半径:
% 克隆扩增:精英个体复制 clone_num 份 top_indices = find(fitness < prctile(fitness, 30)); % 前 30% 作为精英 clones = repmat(pop(top_indices, :), clone_num, 1); % 自适应变异:变异强度随代数衰减 mutation_strength = 0.1 * (1 - gen / max_gen); clones = clones + mutation_strength .* randn(size(clones)); % 禁忌筛选:剔除不满足辐射状的个体 feasible_clones = check_radial(clones); % 精英保留:原精英 + 可行克隆组成新种群 new_pop = [pop(top_indices, :); feasible_clones];prctile(fitness, 30)取适应度排名前 30% 的个体当精英,这个比例不是拍脑袋定的——太低(比如 10%)会导致种群多样性不足,容易早熟收敛;太高(比如 50%)又会让克隆扩增失去聚焦效果。实测 30% 在 33 节点系统上平衡得比较好。
mutation_strength随代数线性衰减是关键设计。早期gen小,变异强度接近 0.1,搜索范围大,利于全局探索;后期接近 0,变异强度小,利于局部精细搜索。如果变异强度固定不变,后期会在最优解附近来回震荡,收敛不了。
check_radial(clones)这一步是禁忌策略的体现——不满足辐射状的克隆直接淘汰,不进入下一代。这比在适应度函数里加惩罚项更直接,因为惩罚项需要调权重,调不好会让可行解和不可行解的适应度混在一起,选择压力失真。
4.2 目标函数设计:网损、负荷削减与开关惩罚的权重平衡
目标函数是三个目标的加权和:网损最小、负荷削减最小、开关动作次数最少。代码里写成:
function obj = create_objective(Pg, Pd, sw_status) % 网损项:支路电流平方乘电阻 loss_term = sum(I.^2 .* r); % 负荷削减惩罚:系数 1000 是经验值 shed_term = 1000 * sum(Pd); % 开关动作惩罚:每次操作计 1000 sw_term = 1000 * sum(sw_status); obj = loss_term + shed_term + sw_term; end这里1000这个系数是调试出来的经验值,不是理论推导的。系数低于 500 的时候,优化器会倾向于“宁可断电也不投备用电源”,因为投备用电源要动开关,开关惩罚相对负荷削减惩罚太便宜了。系数高于 2000 又会过度保守,明明有更优的拓扑也不敢切,因为怕多动一次开关。实际调的时候建议先固定网损项权重为 1,然后单独调shed_term和sw_term的系数,看重构方案在“保供电”和“少操作”之间的取舍是否符合预期。
有个坑要注意:三个目标的量纲不一样。网损是 kW 级别,负荷削减也是 kW,但开关动作次数是无量纲的整数。直接加权求和会让量纲大的项主导优化方向。代码里没做归一化,靠系数硬调,这在 33 节点系统上能跑通,但换到更大的系统(比如 69 节点或 118 节点)就得重新调系数。更稳的做法是先做 min-max 归一化再加权,不过那会增加代码复杂度,看你的精度要求。
4.3 收敛判据与迭代日志
迭代终止条件是两个:gap < 1e-4或者iter >= max_iter。gap是当前最优解和上一代最优解的相对变化量,小于1e-4认为收敛了。日志打印在循环里:
while gap > 1e-4 && iter < max_iter fprintf('第%d轮迭代,间隙值%.4f\n', iter, gap); % ... 克隆、变异、筛选、评估 ... gap = abs(best_obj - prev_best) / abs(prev_best); iter = iter + 1; endmax_iter默认设 100,33 节点系统一般 20 到 30 轮就收敛了。如果跑到 100 轮还没收敛,大概率是变异强度衰减太慢或者种群规模太小,可以试着把clone_num从默认值调大,或者把mutation_strength的初始值从 0.1 降到 0.05。
日志里gap的下降曲线能反映算法行为。正常情况是前几轮快速下降,中间平缓,最后在1e-4附近震荡几下退出。如果gap一直不降或者来回跳,检查check_radial是不是把太多可行解误判成不可行,导致种群多样性枯竭。
5. 避坑与排查:五个真实踩过的坑
5.1 求解器报“SOC 约束不支持”
现象:optimize返回problem = 4,提示二阶锥约束类型不被当前求解器识别。
原因:YALMIP 版本过旧,或者sdpsettings里 solver 设成了不支持 SOCP 的求解器(比如fmincon)。fmincon只能解非线性规划,不认锥约束。
解决:确认 YALMIP 是 2021 之后的版本,solver 设成cplex、gurobi或mosek之一。跑yalmiptest看这三个里哪个显示found,用那个。
5.2 重构结果出现环,但check_radial返回 true
现象:拓扑图上明显有个环,但标题显示“校验通过”。
原因:check_radial里用rank(adj) == n_node - 1判连通性,如果邻接矩阵构建时把故障支路也算进去了(故障支路阻抗大但状态位还是 1),rank会虚高,误判成连通。
解决:构建邻接矩阵前先把故障支路的状态位置 0,或者用branch(:,4) < 1e5作为闭合判据,别只看状态位。代码里build_adjacency函数要同步改。
5.3 电压越限但目标函数没惩罚
现象:重构方案里某些节点电压低于 0.9 p.u.,但优化器认为这是最优解。
原因:目标函数只加了网损、负荷削减、开关惩罚三项,没加电压越限惩罚。约束里虽然有u >= 0.81(0.9 的平方),但如果松弛后求解器在边界附近找到可行解,电压可能贴着下限走。
解决:在约束里加u >= 0.81和u <= 1.21(对应 0.9 到 1.1 p.u.),或者在目标函数里加电压偏差惩罚项。前者更硬,后者更软,看你对电压质量的容忍度。
5.4 切换求解器后结果差异大
现象:CPLEX 跑出来的重构方案和 GUROBI 跑出来的不一样,网损差了好几个 kW。
原因:不同求解器对二阶锥松弛的紧度处理不同,CPLEX 默认的锥松弛可能比 GUROBI 松,导致松弛后的解偏离原始非凸问题的最优解。
解决:在sdpsettings里加'cplex.preind'或'gurobi.method'参数调紧松弛,或者用'solver'同时指定两个求解器做交叉验证。如果差异超过 5%,说明松弛不够紧,得检查 SOC 约束有没有写漏。
5.5 含 DG 的主动配电网潮流不收敛
现象:接入分布式电源后,build_socp_model报潮流方程无解。
原因:DG 的注入功率是负的负荷,如果代码里Pd只允许非负,DG 节点会被当成“负负荷”处理,但约束里没加反向功率项,导致功率平衡方程矛盾。
解决:在节点功率平衡里加Pg项,允许Pg取负值表示 DG 注入。同时检查u的下界,DG 接入点电压可能被抬升,u <= 1.21的上界要留够余量。
6. 二次开发:从 33 节点到自定义网架,以及一个调试习惯
拿到这份代码之后,最常见的二次开发需求是换网架。case33bw()里的节点坐标和支路参数是硬编码的,换成 69 节点或者你们自己的实际配电网,得改三处:节点数n_node、支路数nbranch、以及branch矩阵的每一行。改完之后第一件事是跑check_radial看初始拓扑是不是辐射状,如果不是,说明支路连接关系录错了。
第二常见的需求是加约束。比如限制开关动作次数不超过 3 次,在约束里加sum(sw_status) <= 3就行。但要注意sw_status的定义——如果它是连续变量(松弛后的),得先做取整或者加sw_status <= 1和sw_status >= 0的边界。我一般会在build_socp_model末尾统一加这类工程约束,不散落在各个函数里,方便后面查。
第三个是可视化定制。visualize_topology里的颜色和线型都是写死的,想改成你们单位的标准配色,直接改colormap和plot的参数就行。但有个坑:scatter画联络开关的时候,MarkerSize设太大(比如 200)会盖住节点编号,设太小又看不清。实测 100 到 150 之间比较合适,具体看你的节点密度。
验证改造是否成功,我习惯用“对照法”:先跑原始 33 节点系统,记下网损和收敛轮数作为基线;改完之后再跑,如果网损变化在 10% 以内、收敛轮数没有翻倍,基本说明改造没引入结构性错误。如果网损突然翻了好几倍,大概率是某条支路的阻抗参数录错了,或者辐射状约束被破坏导致求解器在可行域边缘找解。
有个调试习惯我保持了挺久:每次改完build_socp_model,先在optimize之前把model.constraints打印出来,数一下约束条数。33 节点系统正常是 200 条左右(每条支路 4 条约束加节点平衡),如果突然变成 150 条或者 300 条,说明某处循环写错了或者约束被重复添加。这个习惯帮我省了不少排查时间——约束条数不对,后面求解结果再“好看”也是假的。
从那以后我每次换网架或者加约束,都强制先跑一遍约束条数核对和辐射状校验,再去看优化结果。希望帮到你。
本文还有配套的精品资源,点击获取