做电力系统规划的同学对PMU（同步相量测量单元）应该不陌生。单台PMU价格不低，安装位置又决定广域测量系统（WAMS）的“视野”边界，所以PMU站址选在哪里，和选多少台一样重要。最佳PMU位置配置（OPP）就是从全网N个候选节点里找出一个最小集合，装上PMU后让整个电网满足完全可观测性。我最近用二进制粒子群优化（BPSO）把整套OPP流程在Matlab里完整实现了一遍，从IEEE 14节点小算例一路测试到118节点大网，这篇就把原理、代码和踩过的坑一次性说透。对正在做WAMS布点、电网状态估计或者智能优化算法对比的同行，应该有一定参考价值。
1. 先把OPP问题看明白：为什么“放哪里”比“放几个”更烧脑
1.1 PMU、可观测性与广域测量系统
传统SCADA系统的远方终端装置（RTU）采集的是慢速、不同步的稳态数据，而PMU靠GPS时钟同步，能在同一时间断面测出所在节点的电压相量和相连支路的电流相量，时间戳一致是它最值钱的地方。WAMS要做全网动态监视、低频振荡辨识、故障定位，基础就是先把这些“同步相量”铺到一组合理的网络位置上去。但一台PMU连同配套的通信、授时和屏柜设备并不便宜，而且变电站的二次仓位、光纤通道资源都有限，不可能在所有节点都装。于是问题变成：哪些节点装了PMU之后，全网每个节点的电压幅值和相角都能被算出来，同时数量最少。
“可观测”这个词，头一次接触的朋友容易想复杂。其实在这个问题里就是一句大白话：每个节点的电压相量，要么被PMU直接量到，要么能通过已知的电气关系推出来。PMU装在某个节点上，能测这个节点的电压、从这个节点出去的支路电流，而线路阻抗是已知参数，根据欧姆定律，相邻节点的电压相量就能算出来。所以一台PMU的实际效果等于“照亮”它自己，以及所有直接连接的邻居节点。系统完全可观测，就是每个节点至少被这样“照亮”一次。
这里有个特别贴切的类比：城市道路监控探头的选点。一个探头装在路口，能看到自身路口和紧邻的路口；目标是找最少的路口装探头，让所有路口都在至少一个探头的视野里。PMU的OPP问题，就是在电力系统拓扑图上做同样的事，只不过节点之间靠输电线路连接，“视野”由电气定律保证。理解了这个类比，后面再看约束矩阵就不会懵。
1.2 OPP问题的数学刻画与约束条件
OPP可以写成0-1整数规划。记决策变量xi，节点i装PMU取1，不装取0。候选集合就是全部N个节点。目标函数很简单：min sum(x)。难点全在约束上。把电网的节点连接关系写成扩展邻接矩阵A，A对角线全为1，A(i,j)=1表示节点i和j相邻。这样Ax的结果是一个N维列向量，第k行的数值表示节点k自己以及它所有邻居中装PMU的总数。所以“每个节点至少被覆盖一次”的约束，就是Ax >= 1，这里的“1”是N维全1向量。
为什么很多教科书用A*x而不是别的写法？因为矩阵乘法的本质是“行与列的点积”。扩展邻接矩阵第i行只有i自己和邻居这些位置是1，与x做点积，统计的正是在节点i自身和邻居中PMU的数量。这个表达式一次矩阵乘法就完成了全网覆盖情况检查，比循环遍历节点要快得多。后面代码里我用向量化写法，也是基于这一点。
还有一个行业里常提的概念叫零注入节点，指既没有发电机也没有负荷的纯联络节点，注入电流为零。它能额外提供一个KCL功率平衡方程，帮助推导相邻节点的电压，因此等效“视野”更大。考虑零注入后，IEEE 14节点系统这样的算例，最优PMU数量还能再压一两台。本文实现的是最常用的基础一阶可观测模型，零注入属于扩展分支，我在第3.2节会给出扩展思路。
到这里你应该也看出来了，OPP本质上是一个NP难组合优化问题。N个节点每个都有装和不装两种状态，解空间是2^N。IEEE 14节点候选解是16384个，暴力枚举勉强可以；到IEEE 118节点就是2^118，超过10^35种方案，常规枚举根本走不到头。这就是为什么需要智能优化算法。
2. 为什么选BPSO：二进制粒子群求解OPP的改造逻辑
2.1 从连续PSO到离散BPSO
粒子群算法（PSO）的老用户都知道经典速度-位置更新公式：v(t+1)=wv(t)+c1r1*(pbest-x)+c2r2(gbest-x)，然后x(t+1)=x(t)+v(t+1)。连续优化里这一套很顺，但OPP的决策变量是“装/不装”的0/1逻辑，不能直接把连续速度加到离散位置上。BPSO的解决办法很巧妙，把速度从“位移”的含义，改成“取1的概率倾向”。
具体做法是，对速度v做一次sigmoid变换：S(v)=1/(1+exp(-v))，得到0到1之间的概率值，再生成一个0到1均匀随机数rand，如果rand小于S(v)，该维度取1，否则取0。v越大，S越接近1，取1概率越高；v越负，取0概率越高。这就是BPSO和PSO的核心差异：连续PSO里位置是确定的数值，BPSO里位置是随机采样出来的0/1序列。
这里必须提醒一个常见误区。很多刚上手的人会把BPSO的速度更新理解成“pbest减x的二进制意义”，其实这个减法在二进制空间里没有严格的数学定义，它只是提供一种数值倾向：当前这一维是1而pbest是0时，差值为-1，把速度往负方向拉，概率式更新时就更倾向取0；反过来也一样。工程上这么用没问题，但心里要清楚，BPSO本质上是用逻辑回归式的概率模型在做搜索，不是真的在连续空间里“移动粒子”。
2.2 粒子编码与适应度函数设计
用BPSO解OPP，粒子编码最自然的方式就是“一个粒子对应一种PMU布置方案”。每个粒子是长度为N的0/1数组，第k位为1表示在节点k装PMU，0表示不装。初始种群建议随机生成，保证既有比较稀疏的个体也有比较稠密的个体。经验上看，如果初始化全部偏向全0，搜索前期都在给惩罚项“还债”；全部偏向全1，种群多样性差，后期很难把多余PMU减下来。
适应度函数我直接给一个通用结构：fitness = numPMU + penalty * unobservedCount。numPMU是当前方案装PMU的台数，unobservedCount是未覆盖节点数，penalty是惩罚系数。只有全网可观测时，unobservedCount为0，适应度退化为真正的目标值。penalty建议至少设为N，稳妥一点取2N。原因很简单：一台PMU最多给目标函数省下的量级是N，如果penalty小于N，算法会发现“把PMU全拆掉、漏掉所有节点”也是低适应度，搜索方向就毁了。我在初版代码里把penalty设为0.5N，结果粒子全部往稀疏解跑，适应度曲线很漂亮，实际每个解都不可行。这个坑印象非常深。
这里还有一个细节：gbest更新时，不要拿“当前迭代的速度采样结果”去覆盖，而要在适应度评估之后，把历史最优快照（比如pbest和gbest）保存下来。BPSO的gbest本质是一组0/1序列，不是连续空间里的点，很多人因为沿用连续PSO的写法，在全局最优更新时直接把某个中间值放进去，导致收敛方向漂移。代码里我建议用pbest(i,:)这样的二进制行向量快照来存。
3. Matlab代码实现：从邻接矩阵到主循环
3.1 输入数据准备与扩展邻接矩阵
下面用IEEE 14节点系统做示范。这个测试系统有14个母线（bus）和20条支路，常见于OPP论文的算例。先把支路表写出来，再用它生成邻接矩阵。我习惯先用zeros(n)建全零方阵，再逐条支路把两个方向都置1，最后把对角线全部置1，表示节点能被自身安装的PMU覆盖。
n = 14; adj = zeros(n); edges = [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]; for k = 1:size(edges, 1) i = edges(k, 1); j = edges(k, 2); adj(i, j) = 1; adj(j, i) = 1; end % 对角线置1：节点自身可被自己安装的PMU直接观测 adj(1:n+1:end) = 1;有两点必须提醒。第一，如果从外部数据文件读支路表，有些公开数据集的节点编号从0开始，读进来之后要整体加1，否则索引越界，或者生成一个莫名其妙的矩阵。第二，邻接矩阵必须是对称的，因为电网输电线路是可逆的，i到j能推电压，j到i也能推。有的参考代码只写了adj(i,j)=1，忘了写反向，覆盖约束就会漏掉一半邻居，最优解会整体偏大。这个错误很隐蔽，因为矩阵不全为零，程序不报错，但结果就是不对。
3.2 适应度与可观测性检验的向量化写法
可观测性检查的核心就一句：obsMap = adj * x，然后判断 all(obsMap >= 1)。obsMap的第r个元素是节点r自身和邻居中PMU总数，只要每个元素都大于等于1，系统就是完全可观测的。
function [obsFlag, unObs] = calcObservability(adj, x) obsMap = adj * x(:); % 向量化统计每个节点的被覆盖次数 obsFlag = all(obsMap >= 1); % 所有节点都被覆盖才算满足 unObs = sum(obsMap < 1); % 未覆盖节点数，用于惩罚项 end这里用矩阵乘法一次算出全网覆盖情况，比for循环逐个节点判断快得多。在小系统上区别不明显，但到了IEEE 118节点、甚至几千节点的区域电网，这种向量化写法能把单次评估时间从毫秒级降到几十微秒，累计到几百代就是从几分钟变成几秒钟。适应度函数在可观测性函数外层包一层即可。
function f = fitness(adj, x, N) penalty = 2 * N; [obsFlag, unObs] = calcObservability(adj, x); if obsFlag f = sum(x); else f = sum(x) + penalty * unObs; end end如果要做零注入节点的扩展，思路是这样：把零注入节点的电流基尔霍夫方程也纳入判断，当一个零注入节点的邻居数量与已知量测数满足一定条件时，这个零注入节点可以提供额外的等式约束，从而把更多未知节点“推”出来。零注入处理通常不改变基础框架，只是在适应度评估里增加一层潮流方程统一性检查。很多论文会把“考虑零注入”的OPP最优解压到更小，就是这个原因。
3.3 BPSO主程序：参数、初始化与迭代
主循环我按通用BPSO流程写。种群规模pop，迭代次数maxIter，学习因子c1和c2都取2.0，惯性权重w从0.9线性下降到0.4，速度上限vmax取4。vmax这一点特别关键：sigmoid函数在v超过正负4之后已经非常接近1或0，再大的速度值只会让位置更新完全确定化，粒子失去随机探索能力。这就是为什么很多BPSO代码把vmax设在4到6之间，不是随便写的。
rng(42); popSize = 30; maxIter = 100; c1 = 2.0; c2 = 2.0; wMax = 0.9; wMin = 0.4; vmax = 4; X = double(rand(popSize, n) > 0.5); % 初始种群 V = zeros(popSize, n); pbestX = X; pbestFit = inf(popSize, 1); gbestX = X(1, :); gbestFit = inf; fitHistory = zeros(maxIter, 1); for t = 1:maxIter w = wMax - (wMax - wMin) * t / maxIter; for i = 1:popSize r1 = rand(1, n); r2 = rand(1, n); V(i, :) = w * V(i, :) + c1 * r1 .* (pbestX(i, :) - X(i, :)) ... + c2 * r2 .* (gbestX - X(i, :)); V(i, :) = max(min(V(i, :), vmax), -vmax); S = 1 ./ (1 + exp(-V(i, :))); X(i, :) = double(rand(1, n) < S); f = fitness(adj, X(i, :), n); if f < pbestFit(i) pbestFit(i) = f; pbestX(i, :) = X(i, :); end if f < gbestFit gbestFit = f; gbestX = X(i, :); end end fitHistory(t) = gbestFit; end简单解读这段循环。速度更新里(pbestX-X)和(gbestX-X)的作用是产生方向性：当前位为1而个体最优位为0，该维速度会偏负，sigmoid后取0的概率更大，粒子就倾向于向优秀方案看齐。随机数r1、r2让每个维度的更新带独立随机性，避免所有粒子同步收敛。位置生成用的是概率式采样，而不是简单的x+v，这一点和连续PSO完全不同。
惯性权重w在循环内线性递减，前期w接近0.9，粒子惯性大，偏向全局探索；后期w接近0.4，粒子慢下来做局部精细调整。这种从粗到细的过渡，和模拟退火的温度下降思路类似。运行结束后，gbestX就是算法给出的最优PMU位置方案，gbestFit是相应PMU数量。
注意：BPSO里gbest和pbest保存的是历史0/1快照，不要用当前代的临时解去覆盖。我在调试时因为在这里写错，收敛曲线一路“假优化”，白烧了一晚上笔记本。
4. 在IEEE标准节点系统上实测：配置结果与调参记录
4.1 基准算例的已知最优与BPSO表现
OPP这个方向有个好处，就是公开文献里有很多验证过的基准结果。不计零注入节点时，IEEE 14节点系统的最优PMU数量是4台，IEEE 30节点系统约10台。我在Matlab R2023b环境下用上面的代码跑IEEE 14，不同随机种子结果略有波动，但绝大部分独立运行都能收敛到4台。比如一组常见优质解是节点2、6、7、9，覆盖范围正好铺满全系统。
| 算例系统 | 节点规模 | 不考虑零注入的文献最优值 | BPSO多次运行通常结果 |
|---|---|---|---|
| IEEE 14 | 14 | 4台 | 4台 |
| IEEE 30 | 30 | 约10台 | 10~11台 |
| IEEE 118 | 118 | 依赖约束假设 | 高于14节点算例，需加大种群和迭代 |
这个解不是唯一解，但可以作为测试算法的基准点。如果你跑出来的配置数量超过4，先不要怀疑算法，检查一下可观测性矩阵或惩罚系数。IEEE 30节点系统规模大一些，粒子长度变成30，我通常把种群提高到40、迭代增加到150。多次独立测试里，BPSO基本能拿到10台或偶尔11台。这里提醒一句：不同论文对算例的约束假设不完全一致，比如有的把电流幅值限制考虑进可观测条件，有的在目标里加通信冗余度，这些都会让最优数量发生变化。比对文献时先确认对方用的是哪套约束，不要只盯数字。
4.2 参数如何微调：种群、迭代、vmax与惩罚系数
给一套我实测下来比较稳的起始参数：
| 参数 | 推荐值 | 说明 |
|---|---|---|
| popSize | 30~50 | 规模越大，种群越大 |
| maxIter | 100~200 | 先加迭代，再加种群 |
| c1, c2 | 2.0 | 经典粒子群参数 |
| w | 0.9降到0.4 | 前期探索，后期收敛 |
| vmax | 4 | 避免sigmoid饱和 |
| penalty | 2N | 必须大于单台PMU可能省下的最大量 |
小算例完全可以跑得动。规模变大时优先加迭代次数，其次是种群规模，最后才考虑调c1/c2。很多人一上来就动学习因子，结果往往不是收敛慢，就是振荡大，因为真正的问题根本不在那里。
vmax的调整要小心。vmax太大，比如设为10，sigmoid饱和严重，粒子维度经常直接定死为0或1，搜索失去随机性；vmax太小，比如1，速度都集中在sigmoid中段，粒子每次都很随缘，收敛会很慢。如果你发现多轮结果差异极大，先看vmax；如果结果稳定但普遍偏大，很可能w衰减太快，后期粒子没有足够扰动跳出局部极值。我自己的做法是先用默认参数跑一遍，再把gbestFit和已知最优解做对比，有差距再针对性调。
4.3 收敛曲线怎么看：早熟与停滞的处理
收敛曲线的形状能透露很多信息。理想情况下，前20代适应度快速下降，后面形成长长的平台期，平台期的高度就是当前找到的最优PMU数量。如果你看到曲线“断崖式”掉到很低，然后一直平坦，要小心：大概率是惩罚系数太小，算法找到了一个“漏覆盖但适应度低”的假最优。
早熟收敛是启发式算法的通病，BPSO也不例外。我常用的三个补救办法：一是重随机重启，连续跑20次独立实验，取所有gbest里的最优；二是对gbest做局部扰动，选几个维度随机翻转0/1再评估，相当于小范围邻域搜索；三是把种群拓扑从全连接改成环形结构，让信息传播慢一点，降低全体同步陷入局部极值的概率。三个方法都不复杂，但实操效果明显，尤其是重启配合多轮试验，基本能把小系统的结果稳定到文献最优。
5. 实操中常见错误与排查技巧
5.1 代码改完结果不对？照这个清单查
我在调试这段代码时踩过的坑，整理成一个自查清单。第一，邻接矩阵对角线是不是置1了？忘了这一步，可观测性判断会漏掉“节点自身被PMU直接观测”的情况，最优解通常偏大。第二，可观测性判断是all(obs>=1)而不是sum(obs)==N？sum判等会错误地把“被覆盖3次的节点”和“漏掉的节点”抵消掉，实际是不可行解被当成可行解。第三，惩罚系数是否大于N？前面强调过，小于N时算法宁愿漏点也要少装PMU。第四，pbest和gbest更新是否存的是历史快照？如果用当前代临时值覆盖，整个收敛方向会乱。第五，位置更新有没有写成X=X+V？那是在连续PSO里的写法，在BPSO里应该用概率采样。
另外有个很现实的问题：Matlab环境本身。这代码用到的基础操作只有矩阵乘法、比较和循环，R2018之后的版本都能直接跑，不需要所谓“新版本专用函数”。如果遇到license manager error -8或者启动闪退，先检查许可证服务，别急着重装。文件路径里有中文或者全角字符，也会让脚本读数据时静默出错。把路径全部改成英文，能省掉很多莫名其妙的报错。
5.2 从14节点到千节点：提速与规模扩展
如果要把算法用到更大规模系统，我建议做三件事。第一，可观测性检查保持向量化，不要因为图方便改成循环内逐节点判断。第二，在适应度函数里提前剪枝，比如sum(x)已经大于当前gbestFit时，可以直接返回一个很大的值，省掉一次完整的覆盖判断。第三，对邻接矩阵用稀疏矩阵存储，IEEE 118节点网架稀疏度很高，稀疏矩阵乘法可以再压掉一截计算时间。实测在118节点算例、种群50、200代的情况下，核心循环时间在几十秒量级，做实验完全够用。
规模上去之后还会遇到一个新问题：不是所有节点都可安装PMU。有些区域电网的通信条件、二次设备仓位不允许新增测点，这时只需要把对应维度固定为0，在初始化时强制置零并在每次位置更新后再次钳制。这个约束加得很自然，因为BPSO的0/1编码天然支持掩码操作，比整数规划改约束来得灵活。
5.3 和精确求解器配合：BPSO不是替代品
最后说一点我自己的体会。在仅有“最小数量+完全可观测”这种线性目标时，ILP精确求解器才是王道，Matlab自带的intlinprog几行就能搞定IEEE 14和30，严格证明最优。BPSO真正的价值在于目标函数会变，比如同时考虑N-1冗余、通信通道容量、故障后网络重构，这时ILP建模会变得相当笨重，而BPSO只需改适应度函数。
所以我现在做项目的基本流程是：先用一个简单的ILP求出理论下界，建立一个“标尺”，再用BPSO在多种运行方式、多种约束组合下快速重算PMU调整方案。二者不是替代关系，是互相验证的关系。如果你也打算在OPP上继续往下做，我建议把零注入、N-1冗余和通信约束这三个扩展方向都留好接口，BPSO的灵活性会让后面的工作省很多力。
