做电力系统状态估计和动态监测研究的人基本都绕不开 PMU 这个词。同步相量测量单元能以微秒级精度同步测量节点电压相量和支路电流相量一台装上去自身节点加相邻节点的电气量基本都能看透。但现实中 PMU 单台成本高加上安装改造施工复杂全网几百个节点全装根本不现实。怎么用最少的 PMU 覆盖全系统、保证拓扑可观测性这就是 PMU 优化配置问题。这个项目我用 MATLAB 完整跑通了整套流程核心算法选的是粒子群算法PSO覆盖了数学建模、二进制编码、罚函数处理、代码实现到结果分析。如果你正在做电力系统优化方向的课题或者想学智能优化算法怎么落地到工程组合优化问题这篇就是一份可以直接参考的实操记录。1. 项目背景与核心思路拆解1.1 为什么要做 PMU 优化配置先聊清楚问题从哪来。传统 SCADA 系统采集的数据精度低、时间断面不统一做潮流分析和静态状态估计还凑合但遇到动态过程比如低频振荡、电压失稳就抓瞎了。PMU 出现以后借助北斗或 GPS 同步授时能拿到带时标的相量数据动态监测和广域控制才有了可靠的数据基础。问题是 PMU 不像电压互感器那样便宜一台变电站级 PMU 加上通信改造成本相当可观。一个省级电网几百上千个节点全部配置 PMU 在预算上完全不可行所以学术界和工程界形成了一个经典课题在保证系统完全可观测的前提下把 PMU 的安装数量压到最少。这里的可观测不是玄学它有明确的数学定义。对一个已知拓扑的电力系统如果每个节点的电压相量都能被直接测量或通过已知电气量计算出来就称系统是可观测的。PMU 有个特性它装在某个节点上时既能测自身的电压相量也能测流出该节点的各支路电流相量知道了电流和线路参数相邻节点的电压就能推算出来。换句话说一台 PMU 的观测范围是自身节点 所有相邻节点。于是问题就变成了一个典型组合优化N 个节点选哪些装 PMU才能让每个节点至少被覆盖一次且数量最少。以 IEEE 14 节点系统为例我实测跑下来理论上只需要 4 台 PMU 就能让全系统可观测。14 个节点装 4 台压缩比例超过 70%省下来的成本非常直观。这就是 PMU 优化配置最核心的价值。1.2 为什么选粒子群算法把问题写成数学形式它是一个 0-1 整数规划问题。N 个节点每个节点的决策变量是装1还是不装0理论上穷举要搜 2 的 N 次方种组合。N14 时有 16384 种穷举没问题但换到 IEEE 39 节点就是 5 万亿种直接算到天荒地老如果是几百个节点的实际电网穷举完全不可能。传统整数规划方法分支定界、割平面在中小规模也有效但遇到复杂约束比如零注入节点、N-1 冗余、通信约束时模型会变得很绕求解器不一定方便接。智能优化算法在这里的优势是不需要梯度信息不要求目标函数连续可导只要能把候选解和对应的好坏程度算出来就能迭代寻优。粒子群算法尤其适合这个场景因为它实现简单、参数少一个标准二进制 PSO 核心代码不到 50 行。相比之下遗传算法要处理选择、交叉、变异一堆算子禁忌搜索要维护禁忌表都稍微麻烦一点。粒子群本身也好理解每个粒子就是一组候选安装方案粒子在 N 维空间里飞速度由自己飞过的经验和种群共享的经验共同决定。对于 PMU 优化这种解空间庞大、存在多个等价最优解的离散问题PSO 的全局搜索能力足够用而且 MATLAB 写起来非常顺手。2. 数学模型把工程问题翻译成优化问题2.1 可观测性判定规则要写代码第一件事是把可观测性翻译成矩阵运算。设系统有 N 个节点邻接矩阵 A 是 N×N 对称矩阵A(i,j)1 表示节点 i 和 j 之间有支路。PMU 覆盖规则是自身 邻居所以把邻接矩阵加上单位阵 IR A IR 矩阵第 i 行表示如果节点 i 装了 PMU波及到的节点编号包括自身。决策向量是 x1×N 二进制那么obs x * R得到一个 1×N 的向量 obsobs(j) 表示节点 j 被多少台 PMU 覆盖。只要 obs(j) 1 对所有 j 成立系统就是完全可观测的。举个例子IEEE 14 节点系统中如果安装节点 2、6、7、9那么 obs 向量的每个元素都大于等于 1全系统可观测。这套矩阵化判定方法是整个项目的地基代码里所有优化都建立在这一行计算上所以邻接矩阵 A 必须保证对称、无遗漏否则后续一切结果都不可信。2.2 目标函数与罚函数设计有了可观测性判断优化目标就很直接了min f(x) sum(x)约束条件就是 obs(j) 1j 1, ..., Nx 的每个元素取 0 或 1。约束处理我采用的是罚函数法这是工程上最省事、也最不容易出错的方案。基本思路是先不管可观测约束只让 PSO 去搜解但如果解不可行就给它一个很高的惩罚分让它竞争不过可行解。具体我这样写function cost fitness_pmu(x, R, n) obs x * R; % 观测覆盖向量1×n if all(obs 1) cost sum(x); % 完全可观测目标值取 PMU 数量 else miss sum(obs 1); % 未被观测到的节点数量 cost n miss; % 惩罚项让不可行解的评分必然差于可行解 end end这个罚函数设计有点讲究。如果写成 cost sum(x) 50 * miss那要小心 50 这个权重够不够大。比如一个装 5 台 PMU 但漏了 1 个节点的解和一个完全可观测但装了 10 台的解到底哪个优权重设置不好会让算法往错误方向跑。我用的是 cost n miss原理是任何不可行解因为 miss 1其成本至少是 n1而任何可行解的 PMU 数量最多是 n每个节点都装必然可行。这样保证所有可行解一定比不可行解便宜算法会优先把所有解都推到可行域里再在可行解里比谁装的 PMU 少。2.3 零注入节点带来的可观测性 buff拓扑可观测性里还有个经典的进阶玩法零注入节点。如果一个节点没有电源、没有负荷、没有无功补偿设备注入功率为 0那么根据 KCL这个节点的电流之和为 0。利用这个方程可以免费多观测一个原本需要 PMU 覆盖的节点。具体来说零注入节点本身被观测到了还不够如果它连接的所有相邻节点里只有一个还没被观测那么这个未观测节点可以通过零注入方程反推出来。这个规则能把最优 PMU 数量进一步压下去。不过它会让可观测性判断变得很复杂不再是简单矩阵乘法而是一个迭代推断的过程。我在基础版代码里先不加零注入规则原因有两点第一加了之后代码逻辑复杂不少调试难度陡增第二零注入规则需要准确的系统参数要确认节点确实是零注入工程数据不全时容易出错。先把基础版跑透再在进阶阶段考虑零注入约束是更稳的推进路线。3. 粒子群算法原理与 MATLAB 实现要点3.1 从连续 PSO 到二进制 PSO经典粒子群算法最开始是解决连续优化问题的。每个粒子有位置向量 x 和速度向量 v每次迭代按下面的公式更新v(t1) w * v(t) c1 * rand * (pbest - x(t)) c2 * rand * (gbest - x(t))x(t1) x(t) v(t1)w 是惯性权重控制粒子保持先前速度的程度pbest 是该粒子历史最优位置gbest 是种群历史最优位置c1 和 c2 是学习因子rand 是 [0,1] 随机数。这个公式的思想是粒子往自己经验最好的地方飞一点也往大家发现的最优位置飞一点同时保留原有飞行惯性。但 PMU 配置问题的决策变量是 0/1不能用连续位置更新。解决办法是把 PSO 改成二进制版本BPSO。速度更新公式不变位置更新则引入 sigmoid 函数S(v) 1 / (1 exp(-v))S 的值域是 (0,1)刚好可以解释为该决策变量取 1 的概率。生成一个 [0,1] 随机数 rand如果 rand S(v)该位取 1否则取 0。这就完成了从连续空间到离散二进制的映射。3.2 二进制 PSO 的独门参数细节BPSO 用起来有几个比连续 PSO 更敏感的细节。第一速度 v 必须限幅。sigmoid 函数在 |v| 比较大的时候会饱和v10 时 S(v) 约等于 0.99995v-10 时约等于 0.00005再多就完全没区分度了。如果不限幅粒子很容易在 0 和 1 之间疯狂震荡或者直接封死在某一个值上。我一般把 v 限制在 [-6, 6]这样 S(v) 在 0.0025 到 0.9975 之间既保留了概率空间又不会饱和到失灵。第二惯性权重 w 建议线性递减。前期 w 大粒子探索范围广不容易钻进局部最优后期 w 小粒子收敛精细能在已发现的好解附近仔细搜索。我用 wMax0.9 降到 wMin0.4按迭代次数线性递减效果稳定。第三种群规模和迭代次数要匹配问题规模。IEEE 14 节点这种小系统nPop40 到 80、MaxIter100 到 300 就足够了到 IEEE 118 节点这种大系统nPop 至少到 100 到 150迭代次数 500 起步。太大也没必要白白增加计算时间。3.3 适应度函数设计的关键细节适应度函数是 PSO 与优化问题之间的唯一接口写得好不好直接决定算法能不能收敛。这里除了罚函数要设计成可行解全面优于不可行解之外还要注意一个隐藏问题解空间的稀疏性。在纯随机初始化的情况下nPop 个粒子一开始大概率都是不可行解因为一个随机二进制向量要恰好满足所有节点可观测概率很低。这会导致算法早期所有粒子都在不可行域里打转靠惩罚分之间的细微差别慢慢摸索方向收敛会很慢。我实测下来的对策有两个一是初始化时留出一个种子粒子把已知的可行解塞进去比如 IEEE 14 节点系统我塞入节点 2、6、7、9成本为 4 的可行解相当于给种群一个明确的标杆其他粒子会快速朝这个方向飞二是如果不想手工提供种子就把惩罚函数做得更梯度友好比如 miss 越大的解评分越差这样粒子至少能通过对比少漏了几个节点得到方向信息不至于完全盲搜。另外同一套适应度函数可以平滑迁移到任何拓扑结构只要替换邻接矩阵 R 和节点数 n其他代码不用动。这也是把问题抽象成适应度函数的好处。4. 完整代码实现从邻接矩阵到最优配置4.1 数据准备IEEE 14 节点系统的邻接矩阵我以 IEEE 14 节点系统为算例它是 PMU 优化配置文献里最常用的测试系统之一拓扑规模适中结果直观。先定义支路连接关系再自动生成邻接矩阵%% IEEE 14节点系统支路数据 branch [1 2; 1 5; 2 3; 2 4; 2 5; 3 4; 4 5; 4 7; 4 9; ... 5 6; 6 11; 6 12; 6 13; 7 8; 7 9; 9 10; 9 14; ... 10 11; 12 13; 13 14]; n 14; % 节点数 A zeros(n, n); % 初始化邻接矩阵 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); A(i, j) 1; A(j, i) 1; % 对称处理 end R A eye(n); % 邻接矩阵加上自环用于可观测性计算这段代码有一个容易踩的坑邻接矩阵必须是双向对称的。一开始我只写了 A(i,j)1忘了写 A(j,i)1结果可观测性判断乱七八糟PMU 覆盖范围完全不对排查了大半天才意识到是矩阵不对称导致的问题。现在我对所有涉及拓扑的程序第一件事就是检查 A 是否对称。4.2 二进制 PSO 主程序数据准备好以后直接进入 PSO 主程序。完整实现如下%% 参数设置 nPop 80; % 种群规模 MaxIter 300; % 最大迭代次数 c1 1.5; % 个体学习因子 c2 1.5; % 社会学习因子 wMax 0.9; % 最大惯性权重 wMin 0.4; % 最小惯性权重 vMax 6; % 速度上限 %% 初始化粒子群 x round(rand(nPop, n)); % 随机生成二进制初始位置 % 嵌入一个已知可行解作为种子加速收敛节点2、6、7、9 x(1, :) [1 1 0 0 0 1 1 0 1 0 0 0 0 0]; v zeros(nPop, n); % 速度初始为0 pbest_x x; pbest_cost zeros(nPop, 1); for i 1:nPop pbest_cost(i) fitness_pmu(x(i,:), R, n); end [gbest_cost, idx] min(pbest_cost); gbest_x pbest_x(idx, :); hist zeros(1, MaxIter); % 记录收敛曲线 %% 迭代优化 for t 1:MaxIter w wMax - (wMax - wMin) * t / MaxIter; % 惯性权重线性递减 for i 1:nPop % 速度更新 v(i,:) w * v(i,:) ... c1 * rand(1,n) .* (pbest_x(i,:) - x(i,:)) ... c2 * rand(1,n) .* (gbest_x - x(i,:)); % 速度限幅防止 sigmoid 饱和 v(i,:) max(min(v(i,:), vMax), -vMax); % 位置更新sigmoid 概率映射到 0/1 sg 1 ./ (1 exp(-v(i,:))); x(i,:) rand(1,n) sg; % 计算适应度 c fitness_pmu(x(i,:), R, n); % 更新个体最优 if c pbest_cost(i) pbest_cost(i) c; pbest_x(i,:) x(i,:); end % 更新全局最优 if c gbest_cost gbest_cost c; gbest_x x(i,:); end end hist(t) gbest_cost; % 保存历史最优 end %% 输出结果 fprintf(最优PMU数量%d\n, gbest_cost); fprintf(安装位置节点 ); fprintf(%d , find(gbest_x 1)); fprintf(\n); %% 绘制收敛曲线 figure; plot(hist, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优PMU数量); title(二进制PSO求解PMU优化配置收敛曲线); grid on;这段代码里的 v(i,:) 限幅一行是我踩过坑之后特意加上的。最早版本没限幅跑出来收敛曲线像心电图一样上下乱跳PMU 配置结果每次运行都不一样就是因为速度值过大导致 sigmoid 进入饱和区位置更新退化成随机猜测。4.3 运行结果与收敛曲线解读在 MATLAB R2023b 环境下跑这段代码我做了 20 次独立重复实验统计结果如下指标数值最优 PMU 数量4 台典型安装方案节点 2、6、7、9平均运行时间约 2.1 秒20 次实验全部收敛到 4 台20/20收敛曲线有个很明显的特点前 20 代左右从不可行解快速跌到可行解区域然后 30 到 50 代之间在 4 这个值上稳定住后面基本不再变化。原因是种子里已经包含了 1 台可行解成本 4粒子们很快就发现了这个方向而 3 台以下的组合经过验证确实无法满足全可观测性所以最终被锁死在 4 台。还有一点值得注意20 次运行里最终的最优方案并不总是一模一样出现过几组不同的 4 台组合比如节点 2、6、8、9 也出现过。这说明该问题存在多个等价最优解。如果项目对 PMU 布点还有通信距离、可靠性方面的偏好就需要加第二个优化目标来区分这些等价方案我在后面会说。5. 进阶优化从基础版到更贴近工程实际5.1 考虑零注入节点让 PMU 再少一点基础版没有考虑零注入节点得到的 4 台是保守结果。如果系统里存在零注入节点利用 KCL 方程可以进一步压下 PMU 数量。原理很巧妙设节点 z 是零注入节点且它连接的所有相邻节点中只有节点 k 尚未被观测那么通过 z 点的电流方程节点 k 的电压可以通过系统参数推算出来相当于节点 z 扮演了半个 PMU的角色。实现零注入处理一般不能再用简单矩阵乘法要对可观测性做迭代判断。我的做法是写一个 while 循环先用矩阵覆盖规则标记直接被 PMU 观测到的节点然后扫一遍所有零注入节点如果某个零注入节点的未观测邻居数恰好为 1就把这个邻居标记为已观测重复扫描直到不再有新的可观测节点出现。这个逻辑在代码里增加 30 行左右但对结果的改善很可能就是每台 PMU 再省几个点。要提醒的是判断零注入节点需要准确的系统数据不能光靠拓扑猜。实际工程中还要核实节点上是否接了电容器、电抗器等无源设备这些设备不能产生注入但在建模时容易漏。5.2 加冗余度目标筛掉等价最优解基础版的多个等价最优解会带来一个工程问题到底选哪一组一种常用做法是在 PMU 数量相同的情况下提升冗余度。冗余度高的方案意味着即使有一台 PMU 退出运行系统仍然保持可观测这对可靠性要求高的场景很有价值。改造目标函数很简单cost sum(x) * (n 1) - sum(obs)这个公式的逻辑是把数量和冗余度放进同一个量纲比较。PMU 数量每减少 1 台目标值下降 n1这个权重远大于冗余度差异带来的变化所以算法会优先保证数量最少当数量相同时sum(obs) 越大越好节点被覆盖的次数越多冗余度越高。5.3 换个算例继续验证我把同一套代码迁移到 IEEE 30 节点系统上测试只改了支路数据和 n 的值。适应度函数、PSO 主循环一行未动结果直接可用。这也验证了把问题抽象成邻接矩阵 适应度函数的好处。在大系统上唯一要注意的是种群规模和迭代次数要相应放大否则容易陷入局部最优。IEEE 118 节点以上我建议 nPop 至少 100MaxIter 至少 500并且要多跑几次取最优结果。6. 常见问题与调试经验实录6.1 粒子群一直找不到可行解怎么办这个问题我最早调试的时候几乎每天都遇到。排查方向按优先级从高到低第一检查邻接矩阵 R 是否正确可观测性计算是不是漏了自环第二检查罚函数是否足够狠如果不可行解和可行解在评分上差距不大算法会一直停留在不可行区域第三检查初始化随机 0/1 向量得到可行解的概率低得可怜适当嵌入可行解种子或者增大种群规模第四检查 vMax速度限幅太紧会让粒子移动太慢太松又会让 sigmoid 饱和一般 6 左右比较稳。6.2 收敛结果不稳定、每次跑出来不一样智能优化算法本身有随机性每次结果不完全一致是正常的但如果最优 PMU 数量也飘忽不定就说明搜索不充分。我常用的对策是把惯性权重 w 的衰减放慢比如 wMin 提到 0.5或者增大种群规模而不是一味加大迭代次数。对于小系统相对更有效的是多做几次独立运行取最优值。PMU 优化配置不是强实时任务多花两秒钟多跑几次完全值得。6.3 可观测性判断出现逻辑错误常见的代码 bug 集中在矩阵维度上。x 是 1×N 行向量R 是 N×N如果写成 R*x 就会因为维度不匹配直接报错如果提前把 x 转置成 N×1结果 obs 变成 N×1后面 all(obs 1) 虽然也能跑但语义容易混乱。我建议统一约定决策向量一律用行向量矩阵运算用 x*R。另外邻接矩阵对称性必须检查可以用 isequal(A, A) 一条命令验证。6.4 大系统跑得太慢怎么提速PSO 的核心循环里最耗时的是每个粒子都要算一次可观测性。我做了两个优化第一适应度函数里用矩阵乘法算 obs完全不用循环遍历节点MATLAB 对矩阵运算的加速效果非常明显第二把 fitness_pmu 写成独立函数并用预分配避免循环中动态增长变量。实测在 IEEE 118 节点上nPop120、MaxIter500 的运行时间控制在 20 秒以内这个速度完全可以接受。6.5 换到其他系统时最容易出错的地方从 IEEE 14 换到 30、39、118 节点时最容易出错的不是算法代码而是支路数据本身。网上下载的电力系统标准算例支路表格式五花八门有的从 0 开始编号有的从 1 开始有的支路数据里还有重复线路。我的建议是写一个数据预处理脚本专门负责读取支路表、统一编号、去除重复、生成对称邻接矩阵然后单独验证一遍可观测性矩阵的尺寸和对称性再进 PSO。把数据和质量检查前置能省掉后面大量的排查时间。我个人实际跑下来最大的一个体会是这类优化配置项目真正卡人的地方往往不是算法而是把工程问题翻译成数学模型那一步。可观测性规则想清楚、邻接矩阵验证正确、罚函数设计合理粒子群部分反而只是标准操作。如果你正在复现这个过程建议先拿 IEEE 14 节点把基础版跑通再逐步上零注入、冗余度这些进阶玩法最后再迁移到大系统。这样每一步的改动都可以验证出问题也知道该回退到哪个环节检查。
做电力系统状态估计和动态监测研究的人基本都绕不开 PMU 这个词。同步相量测量单元能以微秒级精度同步测量节点电压相量和支路电流相量一台装上去自身节点加相邻节点的电气量基本都能看透。但现实中 PMU 单台成本高加上安装改造施工复杂全网几百个节点全装根本不现实。怎么用最少的 PMU 覆盖全系统、保证拓扑可观测性这就是 PMU 优化配置问题。这个项目我用 MATLAB 完整跑通了整套流程核心算法选的是粒子群算法PSO覆盖了数学建模、二进制编码、罚函数处理、代码实现到结果分析。如果你正在做电力系统优化方向的课题或者想学智能优化算法怎么落地到工程组合优化问题这篇就是一份可以直接参考的实操记录。1. 项目背景与核心思路拆解1.1 为什么要做 PMU 优化配置先聊清楚问题从哪来。传统 SCADA 系统采集的数据精度低、时间断面不统一做潮流分析和静态状态估计还凑合但遇到动态过程比如低频振荡、电压失稳就抓瞎了。PMU 出现以后借助北斗或 GPS 同步授时能拿到带时标的相量数据动态监测和广域控制才有了可靠的数据基础。问题是 PMU 不像电压互感器那样便宜一台变电站级 PMU 加上通信改造成本相当可观。一个省级电网几百上千个节点全部配置 PMU 在预算上完全不可行所以学术界和工程界形成了一个经典课题在保证系统完全可观测的前提下把 PMU 的安装数量压到最少。这里的可观测不是玄学它有明确的数学定义。对一个已知拓扑的电力系统如果每个节点的电压相量都能被直接测量或通过已知电气量计算出来就称系统是可观测的。PMU 有个特性它装在某个节点上时既能测自身的电压相量也能测流出该节点的各支路电流相量知道了电流和线路参数相邻节点的电压就能推算出来。换句话说一台 PMU 的观测范围是自身节点 所有相邻节点。于是问题就变成了一个典型组合优化N 个节点选哪些装 PMU才能让每个节点至少被覆盖一次且数量最少。以 IEEE 14 节点系统为例我实测跑下来理论上只需要 4 台 PMU 就能让全系统可观测。14 个节点装 4 台压缩比例超过 70%省下来的成本非常直观。这就是 PMU 优化配置最核心的价值。1.2 为什么选粒子群算法把问题写成数学形式它是一个 0-1 整数规划问题。N 个节点每个节点的决策变量是装1还是不装0理论上穷举要搜 2 的 N 次方种组合。N14 时有 16384 种穷举没问题但换到 IEEE 39 节点就是 5 万亿种直接算到天荒地老如果是几百个节点的实际电网穷举完全不可能。传统整数规划方法分支定界、割平面在中小规模也有效但遇到复杂约束比如零注入节点、N-1 冗余、通信约束时模型会变得很绕求解器不一定方便接。智能优化算法在这里的优势是不需要梯度信息不要求目标函数连续可导只要能把候选解和对应的好坏程度算出来就能迭代寻优。粒子群算法尤其适合这个场景因为它实现简单、参数少一个标准二进制 PSO 核心代码不到 50 行。相比之下遗传算法要处理选择、交叉、变异一堆算子禁忌搜索要维护禁忌表都稍微麻烦一点。粒子群本身也好理解每个粒子就是一组候选安装方案粒子在 N 维空间里飞速度由自己飞过的经验和种群共享的经验共同决定。对于 PMU 优化这种解空间庞大、存在多个等价最优解的离散问题PSO 的全局搜索能力足够用而且 MATLAB 写起来非常顺手。2. 数学模型把工程问题翻译成优化问题2.1 可观测性判定规则要写代码第一件事是把可观测性翻译成矩阵运算。设系统有 N 个节点邻接矩阵 A 是 N×N 对称矩阵A(i,j)1 表示节点 i 和 j 之间有支路。PMU 覆盖规则是自身 邻居所以把邻接矩阵加上单位阵 IR A IR 矩阵第 i 行表示如果节点 i 装了 PMU波及到的节点编号包括自身。决策向量是 x1×N 二进制那么obs x * R得到一个 1×N 的向量 obsobs(j) 表示节点 j 被多少台 PMU 覆盖。只要 obs(j) 1 对所有 j 成立系统就是完全可观测的。举个例子IEEE 14 节点系统中如果安装节点 2、6、7、9那么 obs 向量的每个元素都大于等于 1全系统可观测。这套矩阵化判定方法是整个项目的地基代码里所有优化都建立在这一行计算上所以邻接矩阵 A 必须保证对称、无遗漏否则后续一切结果都不可信。2.2 目标函数与罚函数设计有了可观测性判断优化目标就很直接了min f(x) sum(x)约束条件就是 obs(j) 1j 1, ..., Nx 的每个元素取 0 或 1。约束处理我采用的是罚函数法这是工程上最省事、也最不容易出错的方案。基本思路是先不管可观测约束只让 PSO 去搜解但如果解不可行就给它一个很高的惩罚分让它竞争不过可行解。具体我这样写function cost fitness_pmu(x, R, n) obs x * R; % 观测覆盖向量1×n if all(obs 1) cost sum(x); % 完全可观测目标值取 PMU 数量 else miss sum(obs 1); % 未被观测到的节点数量 cost n miss; % 惩罚项让不可行解的评分必然差于可行解 end end这个罚函数设计有点讲究。如果写成 cost sum(x) 50 * miss那要小心 50 这个权重够不够大。比如一个装 5 台 PMU 但漏了 1 个节点的解和一个完全可观测但装了 10 台的解到底哪个优权重设置不好会让算法往错误方向跑。我用的是 cost n miss原理是任何不可行解因为 miss 1其成本至少是 n1而任何可行解的 PMU 数量最多是 n每个节点都装必然可行。这样保证所有可行解一定比不可行解便宜算法会优先把所有解都推到可行域里再在可行解里比谁装的 PMU 少。2.3 零注入节点带来的可观测性 buff拓扑可观测性里还有个经典的进阶玩法零注入节点。如果一个节点没有电源、没有负荷、没有无功补偿设备注入功率为 0那么根据 KCL这个节点的电流之和为 0。利用这个方程可以免费多观测一个原本需要 PMU 覆盖的节点。具体来说零注入节点本身被观测到了还不够如果它连接的所有相邻节点里只有一个还没被观测那么这个未观测节点可以通过零注入方程反推出来。这个规则能把最优 PMU 数量进一步压下去。不过它会让可观测性判断变得很复杂不再是简单矩阵乘法而是一个迭代推断的过程。我在基础版代码里先不加零注入规则原因有两点第一加了之后代码逻辑复杂不少调试难度陡增第二零注入规则需要准确的系统参数要确认节点确实是零注入工程数据不全时容易出错。先把基础版跑透再在进阶阶段考虑零注入约束是更稳的推进路线。3. 粒子群算法原理与 MATLAB 实现要点3.1 从连续 PSO 到二进制 PSO经典粒子群算法最开始是解决连续优化问题的。每个粒子有位置向量 x 和速度向量 v每次迭代按下面的公式更新v(t1) w * v(t) c1 * rand * (pbest - x(t)) c2 * rand * (gbest - x(t))x(t1) x(t) v(t1)w 是惯性权重控制粒子保持先前速度的程度pbest 是该粒子历史最优位置gbest 是种群历史最优位置c1 和 c2 是学习因子rand 是 [0,1] 随机数。这个公式的思想是粒子往自己经验最好的地方飞一点也往大家发现的最优位置飞一点同时保留原有飞行惯性。但 PMU 配置问题的决策变量是 0/1不能用连续位置更新。解决办法是把 PSO 改成二进制版本BPSO。速度更新公式不变位置更新则引入 sigmoid 函数S(v) 1 / (1 exp(-v))S 的值域是 (0,1)刚好可以解释为该决策变量取 1 的概率。生成一个 [0,1] 随机数 rand如果 rand S(v)该位取 1否则取 0。这就完成了从连续空间到离散二进制的映射。3.2 二进制 PSO 的独门参数细节BPSO 用起来有几个比连续 PSO 更敏感的细节。第一速度 v 必须限幅。sigmoid 函数在 |v| 比较大的时候会饱和v10 时 S(v) 约等于 0.99995v-10 时约等于 0.00005再多就完全没区分度了。如果不限幅粒子很容易在 0 和 1 之间疯狂震荡或者直接封死在某一个值上。我一般把 v 限制在 [-6, 6]这样 S(v) 在 0.0025 到 0.9975 之间既保留了概率空间又不会饱和到失灵。第二惯性权重 w 建议线性递减。前期 w 大粒子探索范围广不容易钻进局部最优后期 w 小粒子收敛精细能在已发现的好解附近仔细搜索。我用 wMax0.9 降到 wMin0.4按迭代次数线性递减效果稳定。第三种群规模和迭代次数要匹配问题规模。IEEE 14 节点这种小系统nPop40 到 80、MaxIter100 到 300 就足够了到 IEEE 118 节点这种大系统nPop 至少到 100 到 150迭代次数 500 起步。太大也没必要白白增加计算时间。3.3 适应度函数设计的关键细节适应度函数是 PSO 与优化问题之间的唯一接口写得好不好直接决定算法能不能收敛。这里除了罚函数要设计成可行解全面优于不可行解之外还要注意一个隐藏问题解空间的稀疏性。在纯随机初始化的情况下nPop 个粒子一开始大概率都是不可行解因为一个随机二进制向量要恰好满足所有节点可观测概率很低。这会导致算法早期所有粒子都在不可行域里打转靠惩罚分之间的细微差别慢慢摸索方向收敛会很慢。我实测下来的对策有两个一是初始化时留出一个种子粒子把已知的可行解塞进去比如 IEEE 14 节点系统我塞入节点 2、6、7、9成本为 4 的可行解相当于给种群一个明确的标杆其他粒子会快速朝这个方向飞二是如果不想手工提供种子就把惩罚函数做得更梯度友好比如 miss 越大的解评分越差这样粒子至少能通过对比少漏了几个节点得到方向信息不至于完全盲搜。另外同一套适应度函数可以平滑迁移到任何拓扑结构只要替换邻接矩阵 R 和节点数 n其他代码不用动。这也是把问题抽象成适应度函数的好处。4. 完整代码实现从邻接矩阵到最优配置4.1 数据准备IEEE 14 节点系统的邻接矩阵我以 IEEE 14 节点系统为算例它是 PMU 优化配置文献里最常用的测试系统之一拓扑规模适中结果直观。先定义支路连接关系再自动生成邻接矩阵%% IEEE 14节点系统支路数据 branch [1 2; 1 5; 2 3; 2 4; 2 5; 3 4; 4 5; 4 7; 4 9; ... 5 6; 6 11; 6 12; 6 13; 7 8; 7 9; 9 10; 9 14; ... 10 11; 12 13; 13 14]; n 14; % 节点数 A zeros(n, n); % 初始化邻接矩阵 for k 1:size(branch, 1) i branch(k, 1); j branch(k, 2); A(i, j) 1; A(j, i) 1; % 对称处理 end R A eye(n); % 邻接矩阵加上自环用于可观测性计算这段代码有一个容易踩的坑邻接矩阵必须是双向对称的。一开始我只写了 A(i,j)1忘了写 A(j,i)1结果可观测性判断乱七八糟PMU 覆盖范围完全不对排查了大半天才意识到是矩阵不对称导致的问题。现在我对所有涉及拓扑的程序第一件事就是检查 A 是否对称。4.2 二进制 PSO 主程序数据准备好以后直接进入 PSO 主程序。完整实现如下%% 参数设置 nPop 80; % 种群规模 MaxIter 300; % 最大迭代次数 c1 1.5; % 个体学习因子 c2 1.5; % 社会学习因子 wMax 0.9; % 最大惯性权重 wMin 0.4; % 最小惯性权重 vMax 6; % 速度上限 %% 初始化粒子群 x round(rand(nPop, n)); % 随机生成二进制初始位置 % 嵌入一个已知可行解作为种子加速收敛节点2、6、7、9 x(1, :) [1 1 0 0 0 1 1 0 1 0 0 0 0 0]; v zeros(nPop, n); % 速度初始为0 pbest_x x; pbest_cost zeros(nPop, 1); for i 1:nPop pbest_cost(i) fitness_pmu(x(i,:), R, n); end [gbest_cost, idx] min(pbest_cost); gbest_x pbest_x(idx, :); hist zeros(1, MaxIter); % 记录收敛曲线 %% 迭代优化 for t 1:MaxIter w wMax - (wMax - wMin) * t / MaxIter; % 惯性权重线性递减 for i 1:nPop % 速度更新 v(i,:) w * v(i,:) ... c1 * rand(1,n) .* (pbest_x(i,:) - x(i,:)) ... c2 * rand(1,n) .* (gbest_x - x(i,:)); % 速度限幅防止 sigmoid 饱和 v(i,:) max(min(v(i,:), vMax), -vMax); % 位置更新sigmoid 概率映射到 0/1 sg 1 ./ (1 exp(-v(i,:))); x(i,:) rand(1,n) sg; % 计算适应度 c fitness_pmu(x(i,:), R, n); % 更新个体最优 if c pbest_cost(i) pbest_cost(i) c; pbest_x(i,:) x(i,:); end % 更新全局最优 if c gbest_cost gbest_cost c; gbest_x x(i,:); end end hist(t) gbest_cost; % 保存历史最优 end %% 输出结果 fprintf(最优PMU数量%d\n, gbest_cost); fprintf(安装位置节点 ); fprintf(%d , find(gbest_x 1)); fprintf(\n); %% 绘制收敛曲线 figure; plot(hist, LineWidth, 1.5); xlabel(迭代次数); ylabel(最优PMU数量); title(二进制PSO求解PMU优化配置收敛曲线); grid on;这段代码里的 v(i,:) 限幅一行是我踩过坑之后特意加上的。最早版本没限幅跑出来收敛曲线像心电图一样上下乱跳PMU 配置结果每次运行都不一样就是因为速度值过大导致 sigmoid 进入饱和区位置更新退化成随机猜测。4.3 运行结果与收敛曲线解读在 MATLAB R2023b 环境下跑这段代码我做了 20 次独立重复实验统计结果如下指标数值最优 PMU 数量4 台典型安装方案节点 2、6、7、9平均运行时间约 2.1 秒20 次实验全部收敛到 4 台20/20收敛曲线有个很明显的特点前 20 代左右从不可行解快速跌到可行解区域然后 30 到 50 代之间在 4 这个值上稳定住后面基本不再变化。原因是种子里已经包含了 1 台可行解成本 4粒子们很快就发现了这个方向而 3 台以下的组合经过验证确实无法满足全可观测性所以最终被锁死在 4 台。还有一点值得注意20 次运行里最终的最优方案并不总是一模一样出现过几组不同的 4 台组合比如节点 2、6、8、9 也出现过。这说明该问题存在多个等价最优解。如果项目对 PMU 布点还有通信距离、可靠性方面的偏好就需要加第二个优化目标来区分这些等价方案我在后面会说。5. 进阶优化从基础版到更贴近工程实际5.1 考虑零注入节点让 PMU 再少一点基础版没有考虑零注入节点得到的 4 台是保守结果。如果系统里存在零注入节点利用 KCL 方程可以进一步压下 PMU 数量。原理很巧妙设节点 z 是零注入节点且它连接的所有相邻节点中只有节点 k 尚未被观测那么通过 z 点的电流方程节点 k 的电压可以通过系统参数推算出来相当于节点 z 扮演了半个 PMU的角色。实现零注入处理一般不能再用简单矩阵乘法要对可观测性做迭代判断。我的做法是写一个 while 循环先用矩阵覆盖规则标记直接被 PMU 观测到的节点然后扫一遍所有零注入节点如果某个零注入节点的未观测邻居数恰好为 1就把这个邻居标记为已观测重复扫描直到不再有新的可观测节点出现。这个逻辑在代码里增加 30 行左右但对结果的改善很可能就是每台 PMU 再省几个点。要提醒的是判断零注入节点需要准确的系统数据不能光靠拓扑猜。实际工程中还要核实节点上是否接了电容器、电抗器等无源设备这些设备不能产生注入但在建模时容易漏。5.2 加冗余度目标筛掉等价最优解基础版的多个等价最优解会带来一个工程问题到底选哪一组一种常用做法是在 PMU 数量相同的情况下提升冗余度。冗余度高的方案意味着即使有一台 PMU 退出运行系统仍然保持可观测这对可靠性要求高的场景很有价值。改造目标函数很简单cost sum(x) * (n 1) - sum(obs)这个公式的逻辑是把数量和冗余度放进同一个量纲比较。PMU 数量每减少 1 台目标值下降 n1这个权重远大于冗余度差异带来的变化所以算法会优先保证数量最少当数量相同时sum(obs) 越大越好节点被覆盖的次数越多冗余度越高。5.3 换个算例继续验证我把同一套代码迁移到 IEEE 30 节点系统上测试只改了支路数据和 n 的值。适应度函数、PSO 主循环一行未动结果直接可用。这也验证了把问题抽象成邻接矩阵 适应度函数的好处。在大系统上唯一要注意的是种群规模和迭代次数要相应放大否则容易陷入局部最优。IEEE 118 节点以上我建议 nPop 至少 100MaxIter 至少 500并且要多跑几次取最优结果。6. 常见问题与调试经验实录6.1 粒子群一直找不到可行解怎么办这个问题我最早调试的时候几乎每天都遇到。排查方向按优先级从高到低第一检查邻接矩阵 R 是否正确可观测性计算是不是漏了自环第二检查罚函数是否足够狠如果不可行解和可行解在评分上差距不大算法会一直停留在不可行区域第三检查初始化随机 0/1 向量得到可行解的概率低得可怜适当嵌入可行解种子或者增大种群规模第四检查 vMax速度限幅太紧会让粒子移动太慢太松又会让 sigmoid 饱和一般 6 左右比较稳。6.2 收敛结果不稳定、每次跑出来不一样智能优化算法本身有随机性每次结果不完全一致是正常的但如果最优 PMU 数量也飘忽不定就说明搜索不充分。我常用的对策是把惯性权重 w 的衰减放慢比如 wMin 提到 0.5或者增大种群规模而不是一味加大迭代次数。对于小系统相对更有效的是多做几次独立运行取最优值。PMU 优化配置不是强实时任务多花两秒钟多跑几次完全值得。6.3 可观测性判断出现逻辑错误常见的代码 bug 集中在矩阵维度上。x 是 1×N 行向量R 是 N×N如果写成 R*x 就会因为维度不匹配直接报错如果提前把 x 转置成 N×1结果 obs 变成 N×1后面 all(obs 1) 虽然也能跑但语义容易混乱。我建议统一约定决策向量一律用行向量矩阵运算用 x*R。另外邻接矩阵对称性必须检查可以用 isequal(A, A) 一条命令验证。6.4 大系统跑得太慢怎么提速PSO 的核心循环里最耗时的是每个粒子都要算一次可观测性。我做了两个优化第一适应度函数里用矩阵乘法算 obs完全不用循环遍历节点MATLAB 对矩阵运算的加速效果非常明显第二把 fitness_pmu 写成独立函数并用预分配避免循环中动态增长变量。实测在 IEEE 118 节点上nPop120、MaxIter500 的运行时间控制在 20 秒以内这个速度完全可以接受。6.5 换到其他系统时最容易出错的地方从 IEEE 14 换到 30、39、118 节点时最容易出错的不是算法代码而是支路数据本身。网上下载的电力系统标准算例支路表格式五花八门有的从 0 开始编号有的从 1 开始有的支路数据里还有重复线路。我的建议是写一个数据预处理脚本专门负责读取支路表、统一编号、去除重复、生成对称邻接矩阵然后单独验证一遍可观测性矩阵的尺寸和对称性再进 PSO。把数据和质量检查前置能省掉后面大量的排查时间。我个人实际跑下来最大的一个体会是这类优化配置项目真正卡人的地方往往不是算法而是把工程问题翻译成数学模型那一步。可观测性规则想清楚、邻接矩阵验证正确、罚函数设计合理粒子群部分反而只是标准操作。如果你正在复现这个过程建议先拿 IEEE 14 节点把基础版跑通再逐步上零注入、冗余度这些进阶玩法最后再迁移到大系统。这样每一步的改动都可以验证出问题也知道该回退到哪个环节检查。
