电网扩展改造是我这几年接触最多的一块业务之一。去年有个配电网项目甲方一口气定了十来台PMU结果安装完成以后调度大屏上仍然有一片区域的状态量是“灰色”的——数据上不来。问题不在装置本身而在安装位置没选对。这个项目后来逼着我重新把最优PMU配置OPP问题完整梳理了一遍最后用二进制粒子群优化BPSO在Matlab里做了全套实现才把选点方案彻底定了下来。今天把这段研究过程和代码思路整理出来给做电网规划、算法仿真或者电力系统毕业论文的朋友做个参考。文章会从OPP问题的数学本质讲起解释为什么BPSO这类离散智能算法在这类组合爆炸问题上比穷举法和连续PSO更好用然后逐步拆解Matlab建模、可观性分析函数、BPSO核心迭代逻辑最后用IEEE标准节点做算例验证并补充几个实际工程中经常遇到但论文里很少提的坑。1. 先说清楚OPP问题你装的不是PMU是“电网的眼睛”1.1 PMU和SCADA的本质区别很多刚接触PMU配置的同学容易把PMU当成SCADA监控与数据采集系统的“升级版遥测终端”但实际上两者完全不是一个量级的东西。SCADA采集的是稳态电气量扫描周期是秒级不同变电站之间的数据在时间轴上根本不是同步的。而PMU相量测量单元依赖GPS授时可以在微秒级别对电压、电流相量做同步采样。也就是说全网的PMU在同一时刻给系统“拍了一张照片”有了这张照片动态过程、扰动传播、功角变化才能被完整还原。问题就出在这里PMU不是白菜价变电站也不是想装就能装。一套三相电压加若干支路电流的PMU装置加上通信、服务器改造和后期运维单点成本不低。所以实际工程里就出现了一个很典型的规划问题——在给定拓扑结构下最少装几台PMU、装在哪些节点才能让全网状态量完全可观测。这就是OPP问题Optimal PMU Placement要解决的核心。1.2 可观性约束的数学表达判断装一组PMU之后系统是否可观测有一个非常简单实用的判据任何一个总线母线节点要么自己装了PMU要么至少有一个相邻节点装了PMU。因为PMU测量的是安装点的电压相量和流出电流相量通过线路的电流和两端电压之间的关系欧姆定律完全可以把相邻节点的电压相量推算出来。用数学语言说假设电网拓扑可以用一个n×n的邻接矩阵A表示A(i,j)1表示节点i和节点j之间存在输电线路并且A不包含自环那么定义位置决策向量x∈{0,1}^nx_i1表示在节点i安装PMU。可观性向量c就可以写成c x A·x按布尔逻辑理解大于0即为可观其中x那一项代表“自己装了PMU”A·x那一项代表“至少有一个邻居装了PMU”。OPP的目标函数和约束条件就变成了min ∑x_is.t. c_i ≥ 1, i 1,2,…,n这个形式非常干净也是一个标准的0-1整数规划问题。后面你写Matlab适应度函数的时候就完全围绕c_i是否全部大于0来做判断。1.3 零注入节点的“免费午餐”上面那个判据有一个被很多教材忽略的例外——零注入节点。如果一个节点没有发电机、没有负荷、没有无功补偿设备仅仅是一个线路汇集点比如某些中间开关站那么它的注入电流为零。根据基尔霍夫电流定律KCL只要它所有相邻节点中有两个及以上可观测就可以通过电流平衡关系把该零注入节点的电压相量间接推算出来不需要它本身装PMU也不算作“不可观测”。零注入节点是OPP问题里的一个经典优化空间。很多IEEE标准系统里都存在这种节点利用好它们可以把最优PMU数量再压下一两个。但说实话这个处理会显著增加代码复杂度如果只是做算法研究第一步可以先不考虑零注入节点先把基本框架跑通后面我在第7章会讲如何用两步法把它补上。2. 算法选型思路探讨为什么在OPP问题上选BPSO而非其他算法2.1 2^N的组合爆炸问题OPP问题本质上是组合优化问题指数爆炸不可避免。以IEEE 57节点系统为例穷举所有可能的位置组合就是2^57约1.44×10^17个方案就算把约束条件剪枝剪掉一部分这个规模也远超一台普通PC能承受的范围。而工程实际中我还遇到过几百个节点的区域电网模型穷举基本属于天方夜谭。确定性方法比如整数线性规划ILP在小规模系统上可以直接调用CPLEX或Gurobi求解速度很快。但问题是很多研究机构没有商业求解器授权纯Matlab环境下的整数规划函数intlinprog虽然也能算但遇到大系统或者要加入各种工程约束比如某些节点禁止安装、单PMU通道数量限制时改模型和调参数都非常痛苦。这种情况下启发式算法就成了更灵活的选择。2.2 连续PSO为什么不能直接用标准PSO是为连续空间优化设计的。粒子的位置是一个连续实数向量通过速度更新在超平面上飞行寻找适应度函数的极小值。但OPP的决策变量是0/1离散值你不能直接把连续PSO跑出来的x0.732这样的小数当成装0.732台PMU。有人说可以加阈值判断连续值大于0.5就置1否则置0。听着像那么回事但实际效果很差——因为标准的连续速度更新公式产生的速度和位置信息是带方向和大小的硬量化之后这个信息就丢掉了粒子群很容易在某个局部解附近来回震荡收敛性非常不稳定。我早期试过这种变通做法十个粒子里有七八个最后收敛到的是次优解。2.3 BPSO的概率编码逻辑BPSO二进制粒子群优化算法最早是Kennedy和Eberhart在1997年提出的思路很巧妙速度向量保留连续值但不再直接代表位置移动量而是代表位置取1的概率。具体来说用sigmoid函数把速度映射到(0,1)区间映射结果就是该维度取1的概率sigmoid(v) 1 / (1 exp(-v))位置更新变成x 1如果 rand sigmoid(v)否则 x 0这个设计保留了连续PSO的搜索方向信息又天然适配0/1编码。粒子速度大则取1的概率高速度负得厉害则取0的概率高。这样每一次迭代都相当于在01空间里按概率抽样既保证了搜索的随机性又保留了群体经验导向。回到选型这个问题上遗传算法GA也能做0/1优化但GA需要设计交叉算子、变异算子、选择策略参数更多调参周期更长。BPSO的迭代公式相对简洁实现起来代码量小而且我在多个算例上的对比结果是对于OPP这种约束相对规整的组合问题BPSO在收敛速度和稳定性上不输GA。这也是我最终选BPSO的直接原因。3. MATLAB建模第一步邻接矩阵、可观性判断与目标函数3.1 数据准备别在邻接矩阵上偷懒OPP问题的所有输入数据归根结底就是一张反映电网拓扑的邻接矩阵。所以项目开始阶段建议把时间花在把拓扑关系搞清楚上这个数据错了后面全白算。以IEEE 14节点系统为例14条母线之间的连接关系可以用下面这种形式表示。我习惯把邻接矩阵写成一个函数文件方便多个算例共用function A ieee14_adjacency() % IEEE 14-bus system adjacency matrix (without self-loop) % 节点编号 1~14 A zeros(14, 14); lines [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(lines,1) i lines(k,1); j lines(k,2); A(i,j) 1; A(j,i) 1; end end如果是导入外部数据——比如从MATPOWER的case文件转换可以用mpc loadcase(case14.m)然后从mpc.branch矩阵中提取首末端节点编号自动生成邻接矩阵。这里有一个细节需要留意mpc.branch的节点编号是1到n_bus的原始编号如果系统里有节点编号不连续的情况最好先做一个重映射否则邻接矩阵维度对不上。3.2 可观性判断的向量化写法判断一组PMU配置是否让全网可观测用循环写很直观但系统大了以后效率低下。向量化写法其实只有一行核心逻辑function cObs check_observability(A, x) % x: 0/1 行向量, 1表示该节点安装PMU n length(x); cObs (A * x(:) 0) | x(:); cObs cObs; % 保持行向量 end这里的逻辑就是之前说的节点i可观测的条件是“x_i1”或者“存在邻居j使x_j1”。用A乘以x向量得到的就是每个节点被多少个装PMU的邻居覆盖的计数大于0就说明被覆盖了再或上自身安装标志即可。需要注意A * x(:)在Matlab里是数值矩阵乘法如果A是稀疏矩阵也没问题结果里存的是被覆盖的邻居数量。有人会问为什么不直接在邻接矩阵对角线加1然后用一次乘法算c (AI)x也可以两种写法等价但是加单位阵的做法在后续做冗余度分析统计每个节点被几台PMU覆盖时反而更清晰因为对角线加1之后矩阵乘出来的数值就是“自身安装数邻居安装数”。3.3 目标函数与罚函数设计别让适应度变成两座大山OPP的目标函数首先是“安装数量最小”。但是直接拿sum(x)做适应度会有一个严重问题成千上万个粒子中大量方案是不可观的如果只看安装数量不可观的那些方案反而“成本更低”优化过程就被带偏了。标准的处理手段是引入罚函数。我会把适应度函数写成这样function cost opp_fitness(x, A, nBus) % 第一部分PMU安装数量成本 cost sum(x); % 第二部分不可观节点惩罚 x x(:); cObs (A * x 0) | x; nUnobs sum(~cObs); if nUnobs 0 cost cost 1e6 * nUnobs; % 惩罚项系数取远大于安装成本的值 end end惩罚系数为什么取1e6而不是10、100因为安装一台PMU的成本在目标函数里是1而一个不可观节点的代价必须大到任何“少装一台PMU导致一个节点不可观”的方案都不可能成为最优解。我用1e6是为了在不同规模系统里都不用改这个参数。实测下来14节点系统、30节点系统、57节点系统用这个系数都没出过问题。如果你想让搜索结果更偏向“部分可观测但成本低”的工程妥协方案可以把这个系数调小但那是另一个优化目标了。4. BPSO核心代码拆解从速度更新到概率翻转4.1 算法参数怎么定BPSO的参数设置继承了标准PSO的经验但有几个地方需要针对0/1问题单独调整。我常用的参数表如下参数推荐值说明粒子数20~40节点数越大粒子数适当增加但一般40足够迭代次数100~200OPP问题拓扑规模有限200代以内基本收敛惯性权重w0.9线性递减到0.4前期全局搜索后期局部精调学习因子c1,c21.49445Clerc经典参数稳定性好速度钳位[-6, 6]防止sigmoid饱和这是关键细节变异概率pm0.01~0.05增强跳出局部最优的能力关于速度钳位这里多说两句。sigmoid函数当v6时输出已经接近1约0.9975v-6时输出接近0约0.0025。如果速度值不加限制地增大到十几甚至几十粒子所有维度的位置概率都会极端化几乎永远取1或者永远取0搜索能力大幅下降。这就是所谓的sigmoid饱和问题。把速度钳位在[-6,6]内既保留了概率的多样性也让粒子在接近收敛时仍然有概率翻转一些位保持探索能力。4.2 完整的BPSO迭代框架我在Matlab里实现的BPSO主循环逻辑核心代码可以浓缩成下面这一段function [gbest, gbestCost, convCurve] bpso_opp(A, nBus, opts) % 参数初始化 nP opts.nParticles; % 粒子数 maxIter opts.maxIter; % 迭代次数 wMax 0.9; wMin 0.4; c1 1.49445; c2 1.49445; vClamp 6; pm 0.03; % 种群初始化每个粒子是一个1×nBus的0/1向量 X randi([0 1], nP, nBus); V zeros(nP, nBus); % 速度初始化为0 pbest X; % 个体历史最优 pbestCost inf(nP, 1); for i 1:nP pbestCost(i) opp_fitness(X(i,:), A, nBus); end [gbestCost, idx] min(pbestCost); gbest pbest(idx, :); convCurve zeros(maxIter, 1); for t 1:maxIter w wMax - (wMax - wMin) * t / maxIter; for i 1:nP % 更新速度 r1 rand(1, nBus); r2 rand(1, nBus); V(i,:) w * V(i,:) ... c1 * r1 .* (pbest(i,:) - X(i,:)) ... c2 * r2 .* (gbest - X(i,:)); % 速度钳位 V(i,:) max(min(V(i,:), vClamp), -vClamp); % 用sigmoid概率更新位置 X(i,:) double(rand(1, nBus) 1 ./ (1 exp(-V(i,:)))); % 变异操作 if rand pm flipIdx randi(nBus); X(i, flipIdx) 1 - X(i, flipIdx); end % 评估新的粒子 newCost opp_fitness(X(i,:), A, nBus); if newCost pbestCost(i) pbestCost(i) newCost; pbest(i,:) X(i,:); end end % 更新全局最优 [bestCostNow, idx] min(pbestCost); if bestCostNow gbestCost gbestCost bestCostNow; gbest pbest(idx, :); end convCurve(t) gbestCost; end end有几个细节值得解释一下。首先是rand(1, nBus) 1 ./ (1 exp(-V(i,:)))这行。Matlab里exp(-V)对向量计算没问题但需要注意当V是负的大数时exp(-V)可能溢出为Inf导致sigmoid逼近1。这正是为什么要做速度钳位——把V限制在[-6,6]exp(-6)0.0025exp(6)约403都不会溢出。其次是变异操作。标准BPSO没有变异算子容易早熟尤其OPP问题的搜索空间是离散的超立方体粒子可能在某个局部最优配置附近反复横跳。我加了pm0.03的小概率变异每个粒子每代有3%的概率随机翻转一个位置。这个变异率不宜过大否则算法退化成了随机搜索收敛速度会明显下降。4.3 为什么初始速度设为零而不是随机数很多人写PSO时喜欢把初始速度设成随机小量这在连续PSO中没问题。但在BPSO中初始速度如果是个比较大的正数或负数会导致第一代粒子位置直接偏向全1或全0种群多样性被破坏。初始速度归零让粒子第一代只依赖随机位置和群体引力进行搜索反而更稳定。我的经验是初始速度设为0初始位置完全随机这样种群覆盖了搜索空间的各个角落。5. IEEE标准节点算例跑了哪些系统结果怎么验证5.1 三个典型系统的测试结果我用这套BPSO代码在IEEE 14、IEEE 30和IEEE 57三个标准系统上做了测试。这里先说明一下下面的数量结果是在“不考虑零注入节点、不要求N-1冗余”的基本条件下的结论系统节点数线路数BPSO找到的最少PMU数文献参考值IEEE 14142044不考虑零注入时IEEE 3030411010IEEE 5757801717~18比如IEEE 14节点系统BPSO数次运行得到的典型最优位置是{4, 7, 9, 13}这样一组节点不同次运行可能会有等价解比如{2, 7, 10, 13}因为对称性。判断结果对不对我的验证方法是先用穷举法在小系统上验证BPSO是否真的找到了全局最优。14节点系统穷举2^14是16384种组合一分钟之内就能跑完。结果BPSO找到的最优解和穷举结果吻合这让我对算法实现有信心了——很多算法“跑出来一个结果”不等于“结果正确”务必回归到小规模验证。30节点和57节点系统无法穷举就用两个办法验证一是把BPSO结果代入可观性检查函数确认所有节点c_i≥1二是和已知文献中的结果对照。我的测试里BPSO每次都能在和文献一致的数量级上收敛。5.2 收敛曲线和参数敏感性在IEEE 30节点系统上我统计了收敛曲线前15代目标值快速下降从初始几十台PMU的不可观高成本方案迅速降到十几台30代之后基本稳定在10台附近。如果后期连续50代没有改进我会适当调高变异概率重新跑一次看能否冲出局部最优。粒子数的影响也很明显。我用10个粒子和40个粒子分别跑10个粒子的情况下有几次收敛到11台局部最优而40个粒子的情况下基本每次都能稳定找到10台。建议在算力允许时粒子数不要低于30尤其是系统规模超过60个节点之后粒子数太少很容易丢失搜索多样性。5.3 对算法结果的“复现检查”这里分享一个我踩过不少次的坑某些节点系统里有“孤岛”或断开的连通分量如果邻接矩阵用错了数据源比如把某个变压器支路漏掉算出来的最优解可能是一个“在断开子网络中装了一堆PMU、主网却完全不可观”的荒谬方案。所以每次跑完最优解我建议把结果可视化出来画一张拓扑图把PMU安装节点标红逐个节点检查是不是真的被覆盖了。可视化代码思路也很简单g graph(A); figure; plot(g, Layout, force); highlight(g, find(gbest1), NodeColor, r, MarkerSize, 8);看着图检查一遍比单纯打印数据靠谱得多。这个问题也是很多纯做算法优化的朋友容易忽略的——毕竟对电力系统拓扑不熟悉很难从结果数组里一眼看出“节点8和节点12之间其实根本没有线路”。6. 工程中会踩的坑可观测不等于可通信通道数不等于母线数6.1 PMU通道容量限制很多研究论文里做的OPP假设“一台PMU可以测量所有相邻支路”这在仿真里没问题但实际装置不是这样的。一台PMU的采集通道数量是有限的典型配置是8通道或16通道。在变电站里每个通道可以接入一组电压互感器或电流互感器的测量信号。如果一个变电站有多条出线而每台PMU最多只能接有限的电流通道那就需要考虑——一个节点装一台PMU可能覆盖不了它所有的邻居节点。举个实际例子某个220kV枢纽站有6条出线PMU需要同时测量母线电压和6条出线的电流如果装置只有8通道刚好够用但如果有10条出线一台PMU就不够了要么装两台要么选测一部分关键线路。反映到OPP建模上就是在可观性判断函数里增加通道约束或者把“一个PMU覆盖多少支路”作为设备参数加进去。6.2 通信网络和同步精度约束PMU的定位是“同步测量”同步依赖GPS/北斗授时信号但数据传输还需要通信网络。OPP问题通常只解决“装在哪个节点”可实际工程里还有一个配套问题PMU离通信主站太远或者光纤链路带宽不够数据传不回来装了等于白装。我们在一个项目中就遇到过这样的问题某山区变电站拓扑上确实是最优PMU安装点但那个变电站的光纤通道老旧升级成本比PMU本身还高。最后被迫把PMU挪到相邻变电站虽然从纯拓扑角度看多装了一台但整体成本反而更低。所以工程上的OPP优化目标函数里最好加入通信链路改造费用的权重不能只看PMU数量。6.3 N-1可靠性别让一个PMU坏了全系统瘫痪最后一个容易被忽略的工程需求是N-1可靠性。OPP基本版本保证的是“所有PMU都正常时全网可观测”但现实中PMU也有故障率。区域电网对可靠性的要求是任意一台PMU故障退出后全网仍然可观测。这个约束在建模上并不复杂只需要在目标函数里增加遍历对每个安装位置x逐一假设某台PMU故障把x对应位临时置0检查剩余配置是否仍全网可观测。如果任何一个位置故障后不可观就给该方案加惩罚。计算量会比基本版大n倍但BPSO本身计算量不大完全可以承受。我在30节点系统上跑过带N-1约束的版本最优PMU数量从10台增加到了15台左右——这是冗余的代价。需要提醒的是这个N-1约束和零注入节点的KCL推算会互相影响你假设某台PMU故障时要考虑是否还能依靠零注入节点推算处理起来比较繁琐。工程实践中如果不需要那么高的可靠性建议先跑不带N-1约束的版本作为成本下限参考。6.4 “禁止安装节点”的掩码处理配电系统里还有一种常见约束某些节点由于变电站空间限制、绝缘配合或者产权原因物理上不允许新增PMU。处理方式是在编码上增加一个掩码向量maskmask(i)0表示节点i禁止安装PMU。在位置更新和初始化时把禁止节点的位置强制设为0X(:, mask 0) 0;并确保速度更新后在强制归零。这个改动的成本几乎为零但在实际项目中特别实用。我经常和同行说BPSO的优势之一就是这种工程约束加进去非常容易而改用整数规划时每加一条约束都要重新构造数学模型两个方案在灵活度上的差距非常明显。7. 把BPSO代码封装成可复用的研究工具7.1 零注入节点的两步法前面说的基本版BPSO不考虑零注入节点想要把结果做得更贴近文献最优值这里说一个可靠的两步法。第一步全系统不考虑零注入节点用BPSO跑出一个基础最优配置。第二步检查每个零注入节点看它相邻的节点中有没有至少两个已经被PMU覆盖或已可观测。如果满足条件就把这个零注入节点标记为“通过KCL可观测”然后尝试把被它“替代覆盖”的那些PMU中的一个移除再验证全网是否依然可观测。这个两步法的好处是不需要改动BPSO的核心代码只需要在BPSO结果基础之上做一个局部精化。我在IEEE 14节点系统上验证过基本版最优是4台PMU经典考虑零注入节点的最优解是3台用两步法确实能把4台压到3台。7.2 多目标扩展不只是最少还要最稳OPP问题在工程上往往是多目标的既要安装数量少又要冗余度高每个节点被尽量多的PMU覆盖有的项目还要求充分考虑通信跳数或者地震、台风等灾害下的脆弱性。这些目标不一定完全冲突但很难用一个加权函数搞定。我后来把代码扩展成了多目标BPSO版本每个粒子的适应度是一个向量[sum(x), -min_coverage, ...]采用多目标粒子群的支配关系思想维护Pareto前沿代码框架和单目标版本是共用的。如果你有这方面的需求建议在单目标版本跑通之后再做扩展不要一步到位写多目标否则排错会很痛苦。7.3 最终输出什么我最终把这套代码封装成了一个接口比较简单的研究工具[gbest, gbestCost, convCurve] bpso_opp(A, nBus, struct(nParticles,40,maxIter,200));输入邻接矩阵和系统规模输出最优安装位置、最小安装数量和收敛曲线。再加上一节可视化和算例验证代码整个研究闭环就完整了。如果把所有细节邻接矩阵生成、可观性检查、BPSO主程序、结果可视化都写成一个脚本整体代码量也就两百行左右维护成本很低。根据自己的实际操作经验再给大家几个小建议第一代码里所有用到rand的地方如果想复现结果建议用rng固定随机种子第二写论文的话要按照“不同系统下的多次运行统计”来呈现结果而不是只给一次运行的最优值审稿人一般看重这个第三如果后续要扩展到更大规模的电网建议把邻接矩阵改用稀疏矩阵存储并把有可能改成向量化的循环全部向量化否则大系统的运行时间会让人等得失去耐心。