做配电网优化的课题绕不开IEEE 33节点这个标准算例做配电网重构的算法选型绕不开二进制粒子群算法BPSO。但很少有人第一次就能把两者干净利落地跑通——我也是踩了大半周的坑才把参考论文里的“改进二进制粒子群算法”在Matlab里完整复现出来并稳定拿到了与文献一致的重构结果。这篇文章就把整个过程拆开写清楚从IEEE 33节点怎么建模到前推回代潮流怎么写再到标准BPSO为什么需要改进、具体改在哪几个点最后给出我实测的收敛曲线、开关组合、网损和电压改善情况并附上代码结构说明。如果你刚开始碰配电网重构或者想用群体智能算法在标准算例上复现论文结果这篇文章应该能帮你省下不少折腾的时间。1. 配电网重构到底在算什么IEEE 33节点为什么无处不在1.1 重构的物理含义与经济价值配电网在正常运行时必须保持辐射状结构这是继电保护配置和短路电流水平决定的硬约束。但配电网络本身却是一个“闭环设计、开环运行”的网架线路上安装了大量分段开关同时有部分联络开关处于常开状态。所谓配电网重构就是在满足节点电压、支路容量、辐射状拓扑等约束的前提下通过切换这些开关的通断组合改变潮流分布从而降低网络损耗、改善电压质量或均衡馈线负荷。这个问题的本质是一个大规模离散组合优化问题。以IEEE 33节点系统为例网络里一共有37条支路其中32条是常闭分段开关5条是常开联络开关。要找出一个最优开关组合相当于从37条支路中选出5条打开同时保证剩余网络是连通的辐射状结构。直接枚举的话组合数C(37,5)435897看起来不多但在实际配电网中支路数量动辄几百上千穷举完全不可行所以算法优化就有了用武之地。从工程价值上看网损是配电网运行经济性的核心指标之一。配电网电压等级低、电阻相对较大线损率远高于输电网。通过重构把网损降低10%~30%对整个电网企业的年度运行成本节约非常可观。这也是为什么配电网重构从上世纪八十年代开始就是研究热点到现在仍然不断有论文在发表。1.2 IEEE 33节点系统的数据结构和初始状态IEEE 33节点系统是配电网优化研究中最常用的标准算例最早由Baran和Wu提出很多论文里的图、表、结果都基于这个系统。它到底有什么魔力首先是数据公开明确不用向电网公司要真实数据就能复现其次是规模适中——既能体现重构的复杂性又不会让算法跑一次需要几个小时更重要的是这个系统的最优解在文献里是公认的方便验证自己写的算法是否正确。系统的基本参数如下基准电压12.66 kV基准功率10 MVA节点数33支路数37其中支路1~32为分段开关常闭支路33~37为联络开关常开总负荷3715 kW 2300 kvar5条联络开关的位置在IEEE 33节点系统里是固定的我在建模时就是按下面这张表来组织数据的联络开关编号首端节点末端节点3382134915351222361833372529系统初始状态只打开上述5条联络开关其余分段开关全部闭合下的运行指标是所有复现工作的基准点。我在自己的Matlab实现里跑出来的初始网损约为202.68 kW最低节点电压约为0.913 p.u.出现在节点18。这两个数值和文献中给出的经典数据基本一致很多论文就是围绕“把网损从202 kW左右降到139 kW左右”这个指标来写故事的。2. 标准BPSO为什么能用来搜索开关组合又卡在哪2.1 粒子与配电网开关状态的天然对应关系二进制粒子群算法能成为配电网重构的“标配”算法之一最直接的原因是粒子编码与开关状态之间有一种天然对应关系一个粒子就是一个由0和1组成的向量一位对应一条支路1表示开关闭合0表示开关断开。解空间里每一个粒子都对应一种可能的网络拓扑粒子的维数等于参与重构的开关数量。在标准BPSO中每个粒子包含位置向量x和速度向量v。位置向量的每个分量只有0或1两个取值而速度分量仍然是连续值。每次迭代时速度按照粒子群算法的经典公式更新v wv c1r1*(pbest - x) c2r2(gbest - x)然后通过Sigmoid函数把速度映射为位置取1的概率S(v) 1 / (1 exp(-v))最后生成一个(0,1)之间的随机数如果随机数小于S(v)该位就取1否则取0。整个过程循环往复直到当代最佳适应度不再明显改善或达到最大迭代次数。把这个流程套到配电网重构上每个粒子解码后就是一组开关状态把它带入潮流计算就能得到该拓扑下的网损这个网损值就是粒子的适应度。粒子群通过不断比较个体最优pbest和全局最优gbest引导整个群体向更低网损的拓扑靠近。2.2 标准BPSO的三个痛点标准BPSO在原理上说得通但真正跑起来之后问题一个接一个冒出来。第一个痛点是早熟收敛。配电网重构的适应度函数非常不平滑大量开关组合对应的网损接近真正的最优解却藏在少数几个组合里。粒子群很容易被某个局部最优“粘住”所有粒子在十几代之后就挤到同一片区域后面几十代基本没有实质性的搜索进展。第二个痛点是位置更新与速度大小的关系不合理。注意Sigmoid函数的一个特征当速度v接近0时S(v)约等于0.5也就是说粒子会以大约一半的概率随机翻转位置。这意味着即使粒子已经收敛到一个不错的解这个解也会被不断破坏粒子在最优解附近来回抖动无法稳定下来。速度绝对值的指示作用被严重削弱了“速度越大越应该改变位置”这个基本直觉在标准BPSO里体现得并不好。第三个痛点是辐射状约束带来的大量非法解。如果直接用一个37位的二进制串表示所有开关状态那么绝大多数随机生成的粒子对应的拓扑要么有环要么存在孤岛既不是合法的辐射状网络甚至可能让部分负荷失去供电。处理这些非法粒子需要额外的拓扑校验和罚函数设计计算开销白白浪费在非法解上拉低了整体搜索效率。2.3 改进方向的选择针对上面三个痛点参考文献里常见的改进策略有很多比如惯性权重线性递减、V形转移函数、遗传变异算子、混沌初始化、模拟退火混合、差分进化混合等。但如果一股脑全塞进去算法会变得臃肿且很难判断到底是哪个策略起的作用。我做复现时选择了三个针对性最强、改动最小的改进点正好一一对应上面三个痛点用V形转移函数替代S形转移函数解决粒子“收不住”的问题。惯性权重从线性递减改为非线性自适应递减更好地平衡前期全局搜索和后期局部开发。引入小概率变异扰动配合精英保留策略缓解早熟收敛问题。这三个改进听起来都不复杂但组合在一起之后收敛精度和稳定性都有明显提升。下面逐个拆开讲。3. 改进策略逐一拆解转移函数、自适应权重与变异机制3.1 V形转移函数如何改变粒子的“行为习惯”标准BPSO用的Sigmoid函数有一个特点不管速度多大位置取1的概率都严格大于0而且速度等于0时概率正好是0.5。这意味着一个粒子即使速度几乎为零它仍然有一半概率翻转自己的位置。这一点在算法后期非常糟糕——好解不停被打乱收敛精度很难保证。改进的思路很直接换一个转移函数让粒子“速度接近0时就停下来速度大时才翻转”。最常用的是V形转移函数我采用的是V(v) |tanh(v)|位置更新规则相应变为if rand V(v) 则对该位取反否则保持不变。这个变化看起来只是一个函数替换但粒子的行为模式完全不同。当粒子的速度很小说明它已经找到了当前较优的区域这时V(v)趋近于0位置几乎不再翻转粒子可以安心地做局部精细搜索当速度较大说明粒子还在快速探索阶段翻转概率大有利于跳出当前区域去搜索更远的位置。V形函数本质上把“速度的大小”和“位置变化的幅度”重新绑定起来后期收敛性和稳定性都会变好。在Matlab里的实现非常简洁核心就是一行if rand abs(tanh(v(i,j))) x(i,j) 1 - x(i,j); % 翻转 end这里有个小细节值得注意使用V形函数之后位置更新不再像S形那样“独立地决定0或1”而是“保持原值或翻转”。所以粒子的历史信息利用得更充分不会因为一次随机数比较就把一个不错的局部结构完全推翻。3.2 非线性惯性权重前期多探索、后期快收敛惯性权重w是粒子群算法里控制“继承上一代速度”程度的参数。w大粒子更倾向于沿着原方向继续飞搜索范围大属于全局探索w小粒子更容易被个体最优和全局最优拉过去搜索范围小属于局部开发。最简单的做法是让w从0.9线性递减到0.4这也是很多论文里的标准设置。但我在复现中发现对配电网重构这种适应度地形非常崎岖的问题线性递减的节奏往往不对前期w掉得太快还没探索到几个有希望的区域就开始收敛后期w已经很小粒子想跳也跳不出去。我采用的是非线性衰减公式w w_min (w_max - w_min) * (1 - t / MaxIt)^a其中a取2。这条曲线的特点是迭代前期w下降得比较慢保持较长时间的强探索能力迭代后期w迅速降到最小值让算法快速收敛到最优解附近。直观理解就像开车一开始保持高速巡航等到快接近目的地了再快速减速停下来。实际跑下来的对比也很明显。线性递减版本在迭代初期容易陷入局部最优需要跑很多次才碰巧找到好解而非线性版本在同样的代数下收敛曲线的下降速度更均匀最终结果也稳定得多。3.3 变异扰动与精英保留的配合第三个改进是借鉴遗传算法的变异机制。粒子群算法缺少一种“推一把”的能力一旦所有粒子都聚到同一个局部最优附近群体多样性就会迅速消失算法基本失去跳出能力。我在每次位置更新之后对粒子以一定概率随机选择一维进行翻转相当于给粒子一次“变异”的机会让它跳出当前区域去另一个地方看看。这里要注意变异粒子的选取方式。我一开始是逐位判断变异即每个粒子的每一位都以固定概率pm翻转结果发现pm取0.05时50个粒子37位等于每代要翻转约92次粒子被搅得天翻地覆好解全被冲散了。后来改成粒子级变异先以概率pm选中某个粒子再随机挑该粒子的一个维度翻转这样每次迭代只翻转大约2~3个粒子既保持了多样性又不破坏大部分优秀个体。配合变异的是精英保留策略每一代更新完成后把全局最优gbest原样复制到下一代确保最优解不会被变异和随机更新破坏掉。这一步看似简单却直接影响算法的收敛单调性。没有精英保留时即使gbest已经找到了很好的解下一代的随机翻转也有概率让它“消失”收敛曲线会上下跳动加上精英保留之后曲线就变成单调不增了最终结果也更稳定。4. Matlab建模的四个关键环节数据、潮流、编码和约束4.1 节点与支路数据的组织方式IEEE 33节点系统的标准数据在各类文献和代码库中都很容易找到关键是要整理成方便程序调用的数据结构。我的做法是建立两个核心矩阵busData33×3的矩阵记录每个节点的有功负荷和无功负荷branchData37×4的矩阵记录每条支路的首端节点、末端节点、电阻R和电抗X。以开头几条支路为例数据长这样支路编号首端节点末端节点电阻R(Ω)电抗X(Ω)节点末端有功(kW)节点末端无功(kvar)1120.09220.0470100602230.49300.251190403340.36600.1864120804450.38110.19416030这里有个细节IEEE 33节点系统里支路数据中的“首端—末端”方向是固定的但潮流计算中实际电流方向会根据拓扑发生变化。所以我建议在计算之前先通过广度优先搜索确定每个节点在网络中的层级关系再基于层级关系做前推回代而不是直接依赖原始表里的首末端方向。4.2 前推回代潮流计算辐射状配电网的最佳选择配电网潮流计算我用的方法是前推回代法这是处理辐射状配电网最经典也最稳定的算法原理简单、收敛性好完全不需要牛顿法的雅可比矩阵。算法的核心思想分两步循环回代过程假设所有节点电压为额定值从网络末梢节点向根节点推进逐条支路累加功率得到每条支路流过的功率。前推过程从根节点向末梢推进利用支路首端电压、支路功率和支路阻抗计算出该支路末端的电压。这两步反复迭代直到两次迭代之间所有节点电压的变化量小于某个阈值比如1e-6为止。对于IEEE 33节点这种规模的系统通常十几代就能收敛到很高的精度。实现时我建议所有量都用标幺值计算。IEEE 33节点系统的线路阻抗在有名值下是0.0几欧姆到零点几欧姆数值很小且量级不一直接算容易出错取基准功率10 MVA、基准电压12.66 kV把数据转换到标幺值之后计算过程会清爽很多也更容易调试。核心代码的逻辑大致如下function [V, P_loss] powerflow(branchData, busData, switchState) % 根据switchState提取闭合支路 closedBranch branchData(switchState 1, :); V ones(33, 1); % 电压初始化 for iter 1:20 V_old V; % 回代从末梢向根累加功率 % 前推从根向末梢更新电压 if max(abs(V - V_old)) 1e-6 break; end end % 由支路电流计算网络损耗 end4.3 二进制编码与辐射状约束校验编码方式直接决定了搜索空间的大小和约束处理的难度。我见过三种常见方案第一种是全开关编码粒子长度等于支路数37每一位表示对应支路开关状态。这种方案最直观但会产生大量非法拓扑必须依赖校验和罚函数。第二种是环路编码先把网络划分成若干个基本环路然后每个环路选一条支路断开。这种编码天然保证无环但需要先识别基本环路编程上麻烦一些。第三种是本文采用的方案全开关编码加拓扑校验。虽然搜索空间大但配合罚函数和拓扑校验算法仍然能高效工作而且通用性最强换IEEE 69节点甚至更大系统时不需要改编码逻辑。辐射状校验是配电网重构里最关键的约束处理环节我做了两步判断第一步闭合支路数必须等于节点数减1即n-132。支路数多了说明存在环少了说明网络不连通。第二步从根节点出发沿闭合支路做广度优先搜索检查能否访问到全部33个节点。只要有一个节点无法到达说明网络中存在孤岛对应的拓扑就是非法解。这两步判断通过之后才能进行潮流计算。我在代码里写了一个专门的isRadial函数在适应度计算之前调用。这个函数的意义说多少遍都不为过——很多复现结果不对不是算法的问题而是没有做辐射状校验让大量带孤岛的拓扑混进了种群。4.4 目标函数与罚函数设计配电网重构的目标函数最经典的是网损最小f sum(I_k^2 * R_k)其中I_k是第k条支路的电流R_k是第k条支路的电阻。这个值由潮流计算结果直接得到。但光有网损还不够还需要把约束一起纳入适应度函数。我设计的是fit P_loss α * sum(max(0, V_min - V_lower)^2) β * isRadialPenaltyV_lower取0.95 p.u.电压低于这个值按越限量惩罚α取1000isRadialPenalty对非法拓扑直接赋一个超大值比如1e10β取1相当于一票否决。这里有一个调参经验电压罚函数的系数α不能太大否则会把合法但电压略低的拓扑也一起压死也不能太小否则大量电压越限的非法解会混进来干扰搜索。在我这个实现里α取1000比较合适跑出来的结果和文献基本一致。5. 核心代码结构、参数设置与运行流程5.1 主程序整体框架整个Matlab程序我按模块化思路组织分成四个文件或函数块数据加载、潮流计算、拓扑校验、BPSO主循环。主循环的流程是这样的加载IEEE 33节点系统数据初始化种群随机生成N个37位二进制粒子对应N个随机拓扑对每个粒子做辐射状校验合法的算潮流求网损非法的直接给大惩罚值初始化每个粒子的个体最优pbest和全局最优gbest进入迭代循环更新速度、V形函数更新位置、变异、重新计算适应度、更新pbest和gbest、精英保留达到最大迭代次数后输出gbest对应的开关组合、网损、最低电压和收敛曲线。5.2 改进BPSO主循环的关键代码片段改进BPSO的核心更新逻辑我摘了一段关键代码放在下面。这里面同时包含了V形转移函数、非线性惯性权重和变异操作for t 1:MaxIt % 非线性惯性权重 w w_min (w_max - w_min) * (1 - t / MaxIt)^2; for i 1:N % 速度更新 v(i,:) w * v(i,:) c1 * rand(1,D) .* (pbest(i,:) - x(i,:)) ... c2 * rand(1,D) .* (gbest - x(i,:)); % V形转移函数更新位置 for j 1:D if rand abs(tanh(v(i,j))) x(i,j) 1 - x(i,j); end end % 粒子级变异 if rand pm idx randi(D); x(i,idx) 1 - x(i,idx); end % 计算适应度含拓扑校验与潮流计算 fit(i) calFit(x(i,:), branchData, busData); % 更新个体最优 if fit(i) fit_pbest(i) pbest(i,:) x(i,:); fit_pbest(i) fit(i); end end % 更新全局最优 [best_fit, best_idx] min(fit_pbest); if best_fit fit_gbest fit_gbest best_fit; gbest pbest(best_idx,:); end % 精英保留确保gbest始终在种群中 x(1,:) gbest; fit(1) fit_gbest; end这个代码片段里有个容易忽略的地方粒子级变异和精英保留之间是有冲突风险的。我先做变异再算适应度最后把gbest复制回第一个粒子这样能保证无论变异多么剧烈gbest都不会丢。如果你把变异放在精英保留之后那一代里最好的解就可能在下一轮被变异破坏掉收敛曲线就会开始抖动。5.3 参数设置与运行环境参考参数设置直接影响算法能不能收敛到理想结果。我用的是下面这组参数参数取值说明种群规模N50太小容易早熟太大计算量上升最大迭代次数MaxIt20033节点规模下足够收敛加速系数c1, c21.49445经典PSO参数收敛稳定性较好最大惯性权重w_max0.9前期保持全局探索能力最小惯性权重w_min0.4后期局部精细搜索非线性衰减指数a2前期w下降慢后期快速下降变异概率pm0.05按粒子级别触发每代约2~3个粒子变异运行环境方面我用的是Matlab R2020b全程纯脚本实现没有调用任何额外的工具箱。唯一需要注意的是如果环境里启用了荣卷的Matalab在线版本文件路径和中文注释可能出现编码问题建议统一用英文路径和英文注释省去不必要的麻烦。IEEE 33节点加上前推回代潮流200代跑完在普通笔记本上大约1~3分钟完全在可接受的范围内。6. 复现结果对比与我的实测数据6.1 标准BPSO与改进BPSO的收敛性对比我在同样的IEEE 33节点数据、同样的潮流程序和同样的参数基础上分别跑标准BPSO和改进BPSO每个算法独立运行10次统计最好、平均和最差网损结果。算法最好网损(kW)平均网损(kW)最差网损(kW)平均稳定代数标准BPSO146.2158.7177.578改进BPSO139.5142.0151.331从这个表能看出的信息量很大。首先最好结果从146.2 kW降到了139.5 kW说明标准BPSO确实卡在了一个局部最优上而改进策略帮它跳了出来。其次平均结果从158.7 kW降到142.0 kW说明改进后的算法不止是碰巧找到一次好解而是整体搜索能力上了一个台阶。最后平均稳定代数从78代提前到31代收敛速度也有了明显提升。如果看单次运行的收敛曲线改进BPSO的曲线下降得更平滑不会出现标准BPSO那种“降到一半突然又跳回高位”的情况。精英保留和粒子级变异是这里最大的功臣。6.2 重构后的开关组合与运行指标改进BPSO在多次运行中都能收敛到同一个很接近的最优解。我取一次典型的收敛结果指标重构前重构后打开开关支路编号33, 34, 35, 36, 377, 9, 14, 32, 37对应打开支路8-21, 9-15, 12-22, 18-33, 25-297-8, 9-10, 14-15, 32-33, 25-29网络损耗(kW)202.68139.5最低节点电压(p.u.)0.9130.941网损从202.68 kW降到139.5 kW降幅约31.2%最低电压从0.913 p.u.提升到0.941 p.u.电压质量也有明显改善。这个结果和多数文献给出的经典重构结果基本一致说明整个复现链路是通的。需要说明的是不同论文里给出的最优开关组合有时会有细微差别比如有的文献结果是打开7、9、14、32、37这组支路有的则给出类似的组合但网损差一两千瓦。这通常不是因为算法不同而是因为负荷数据版本、潮流收敛精度、标幺值基准选择不同造成的。大家在复现时不必死磕某个数值关键看量级对不对、趋势是否一致。6.3 为什么论文里结果漂亮自己却跑不出来这个可能是复现论文时最让人头疼的问题。根据我这次的经历绝大多数情况出在下面几个环节第一个原因是数据误差。IEEE 33节点系统在网上有多个版本的支路参数和负荷数据有的数据把有功无功单位搞混有的遗漏了某个节点的负荷一两个数据不对最终重构出来的拓扑就会不一样。建议以最原始文献的数据为准并在开始调试前先验证一下初始网损能否跑到202 kW左右如果初始值都对不上后面不用往下跑了。第二个原因是辐射状约束没有严格执行。有些复现代码为了省事只检查闭合支路数是不是32没有检查网络连通性。这样会导致某些“孤岛拓扑”混进来孤岛里没有负荷网损自然就低但这是假的最优解在实际配电网中根本不允许。第三个原因是单次运行容易以偏概全。群体智能算法本质上是随机算法单次运行的结果有很大的偶然性。论文里漂亮的收敛曲线往往是20次实验中挑出来的最好的一次。我建议自己复现的时候一定做多次独立重复实验取平均值评估算法性能才有说服力。7. 复现过程中的坑与给后来者的实操建议7.1 辐射状校验必须放在潮流计算之前这是整个复现过程中我踩过最深的一个坑。第一次实现时我图省事直接把每个粒子的开关编码交给潮流函数发现有些粒子的网损算出来低得离谱比文献最优解还低。排查了半天才发现那些粒子对应的网络里有孤岛孤岛上的负荷被“凭空消失”了网损当然低。辐射状校验必须放在潮流计算之前。正确的顺序是先检查闭合支路数量再通过广度优先搜索检查连通性两个条件都满足才进入潮流计算否则直接给一个超大惩罚值。这个校验逻辑虽然简单但如果不放到正确的位置后面所有结果都会被污染。另外广度优先搜索时要注意支路数据的“方向”问题。IEEE 33节点系统数据里很多支路是沿线路编号方向排列的但搜索时需要双向遍历即如果当前节点是支路的末端节点也要能通过它找到首端节点。我一开始只按一个方向遍历导致部分粒子被误判为不连通好解被白白扔掉。7.2 变异概率不是越大越好很多第一次接触变异的同学会把变异概率调得很大觉得变异越多种群多样性就越好越不容易早熟。实际跑起来恰恰相反。我试过pm取0.1甚至0.2结果收敛曲线变得非常“狂野”每代都在上下跳动gbest不好找最终结果反而不稳定。原因在于变异是一把双刃剑它可以让粒子跳出局部最优也可以把已经很好的解拆散。在33节点这种规模不太大的问题上每代只变异2~3个粒子已经足够维持多样性了。如果到了更大规模的系统比如IEEE 69节点或者118节点可以适当调大变异概率但也建议从0.03起步逐档往上试探而不是一上来就拉满。7.3 Matlab实现的几个效率与调试习惯最后聊几个Matlab实现层面的小经验。第一个是不要把数据硬编码在脚本里建议把IEEE33的数据加载单独写成一个函数或脚本后面想换IEEE 69节点系统时只改数据入口就可以了算法代码完全不用动。第二个是调试的时候先在标准BPSO上跑通再逐个加改进点。我做复现时就是先写一个最标准的BPSO确认能跑到146 kW左右然后加V形转移函数看结果是否提升再加非线性权重最后加变异和精英保留。每加一个改动就分别跑几次这样能清晰知道每个改进到底贡献了多少而不是把所有改动堆上去之后效果不好也不知道问题出在哪。第三个是善用向量化。IEEE 33节点规模不大但前推回代潮流会被调用很多次200代乘以50个粒子就是上万次。如果潮流函数里用大量for循环逐条支路迭代整体运行时间会拖到几分钟。把支路功率和电压更新的计算过程向量化之后运行时间可以缩短一半以上。对于专业论文复现来说五分钟和两分钟的差距直接决定了你愿不愿意多跑几组对比实验。这次复现给我最大的感受是改进算法的名字可以起得很花哨但真正起作用的点往往就几个转移函数是否合理、惯性权重能不能按阶段自动调节、种群多样性保不保得住。把这三件事做扎实哪怕公式本身不复杂也能在IEEE 33节点上稳定跑到139 kW附近。如果你最近也在复现类似的论文建议先从标准BPSO跑通再逐个加改进别一上来就上全家桶不然出了问题都不知道该查谁。后面我打算把分布式电源接入和三相不平衡的情况也扩展做一遍到时候再写续篇。