改进二进制粒子群算法在IEEE33节点配电网重构中的Matlab复现实践

1. 项目定位与核心价值

1.1 这个项目解决什么问题

做配电网重构的朋友，应该绕不开IEEE33节点和二进制粒子群算法这两个关键词。我最近完整复现了一篇改进二进制粒子群算法用于配电网重构的核心论文，目标很直接：在Matlab环境下，用IEEE33节点作为标准算例，把网络重构抽象成0-1开关组合优化问题，然后用改进后的BPSO去搜索低损耗、电压质量更好的运行方式。这个工程几乎覆盖了配电网重构入门要碰到的所有关键环节，包括拓扑连通性、前推回代潮流、二进制编码、约束处理、收敛性分析和算法对比。不管你是电气工程的研究生、做课程设计的本科生，还是刚接触电网规划的工程师，这套复现思路和代码框架都值得完整过一遍。

为什么它值得写成一篇文章？因为很多论文只给公式和结论，真正落实到Matlab时会遇到一堆隐性问题。比如编码后粒子对应的开关状态不满足辐射状约束、潮流迭代半天不收敛、改进算法和标准算法收敛曲线差异不明显、IEEE33节点的支路编号和联络开关对不上等等。这些问题单纯看论文是看不出来的，只有亲手跑一遍才能理解。下面我按从数学模型到代码实现，再到实验对比的顺序，把整个复现过程拆开讲。

1.2 为什么选IEEE33节点当“试验田”

IEEE33节点系统在配电网重构算法论文里出现频率极高，几乎就是“试验田”。它包含33个节点、37条支路，其中32条分段开关支路、5条联络开关支路，规模不大不小。节点太多，像IEEE123节点，算法迭代一次潮流计算的时间会明显增加，调试不方便；节点太少，又体现不出配电网重构中辐射状约束和电压变化的基本特征。IEEE33节点在常规额定负荷下，原始网损大约202.68kW，最低电压接近0.9131pu，重构后有明显的降损和电压提升空间，很多文献都把它的结果当作基准，方便我们对照。

另一个关键原因是IEEE33节点的参数公开稳定，不同论文里差别极小。这就给复现提供了“锚点”：如果你的算法跑下来的结果和文献对不上，可以先怀疑编码或潮流计算有问题，而不是怀疑算例本身。对新手来说，先在IEEE33节点上把基准跑通，再扩展到复杂网络或分布式电源场景，是效率最高的一条路。我自己的习惯也是每次研究新算法，都先在IEEE33上重跑一遍原始版本，作为“对照组”。

2. 配电网重构的数学模型与评价指标

2.1 目标函数：让网损成为唯一指挥棒

配电网重构最经典的目标是让有功网损最小。对于一个闭合支路集合确定的辐射网，网损可以写成支路潮流和电压的函数：

min f = Σ_{i=1}^{N_b} [ (P_i^2 + Q_i^2) / V_i^2 ] * R_i

其中P_i、Q_i是第i条支路首端流过的有功和无功功率，V_i是首端节点电压，R_i是支路电阻。这个公式把重构问题“翻译”成可以计算的适应度：网损小的开关组合得分高。实际处理时，也可以直接用支路电流平方乘电阻，但两种公式在电压偏离1.0pu的时候会有微小差异，复现论文时要看清原文用哪一种。

很多核心论文第一步先用单目标网损最小，后面再扩展成“网损+电压偏差+开关动作次数”的多目标。我在复现时也建议先只做单目标。原因很简单：多目标涉及到Pareto前沿、权重系数调节，变量多了容易把问题搞复杂。先把网损降到文献那个水平，再逐步加目标，每加一个目标都记录对解集分布的影响，这样后面写论文或者改代码都有据可查。

2.2 约束条件不能只做摆设

IEEE33节点重构问题看着简单，真正的难点在约束。至少有三条硬约束必须处理：

1. 拓扑约束：重构后网络必须是辐射状，既要连通，又不能出现环网。用支路数量来看，33个节点的树必须正好有32条闭合支路，同时从根节点能到达所有节点。
2. 电压约束：通常要求各节点电压在0.90pu到1.05pu之间，低于下限的要施加惩罚。
3. 支路容量约束：每条支路电流或视在功率不能超过允许上限，过载也需要惩罚。

这些约束不能只写在论文里，在算法实现中要通过适应度函数或编码方式体现。常见做法是“编码约束+惩罚函数”：编码约束尽量让粒子生成的开关组合满足闭合支路数为32，连通性则在适应度函数里检查，不连通或过电压的粒子直接给一个很大适应度。如果惩罚太轻，不可行解可能比可行解得分还高，算法会跑偏；惩罚太重，又会压缩搜索空间。一般我取比最大可行网损高两个数量级，比如可行网损在0.14MW左右，那么惩罚项取10^4或更高。

2.3 潮流计算怎么嵌进来

配电网重构需要反复评估大量开关组合，一次优化可能要跑几千次甚至上万次潮流。IEEE33节点是典型的放射状网络，用牛顿-拉夫逊法也能算，但每次都要形成雅可比矩阵，速度偏慢，而且对辐射状网络的节点编号顺序比较敏感。我更推荐前推回代法：从末端节点开始向根节点逐段累加支路功率，再从根节点向末端逐段计算节点电压，反复迭代几次就能收敛。它的计算量小，对配电网的R/X比值大、电压等级不高的场景适应得也很好。

在Matlab里实现前推回代的核心是先把网络按支路关系整理好，知道每一条支路对应的“首端节点”和“末端节点”。前推时从网络末端开始，把下层节点的负荷和支路功率逐层往根节点累计；回代时从根节点已知电压（通常1.0pu）开始，沿支路推算下一层节点电压。迭代到最大电压偏差小于1e-6就认为收敛。这套逻辑说白了就是把“树的叶子往根部汇总、再从根部往叶子分配”重复几次，非常符合配电网单电源辐射状的结构。

3. 二进制粒子群算法与改进方向

3.1 标准BPSO是怎么工作的

标准粒子群算法针对连续变量，粒子速度直接加到位置上。但配电网重构的开关状态只有开/合两种，需要二进制粒子群算法。BPSO的核心思路是保留粒子的速度概念，但位置不再是连续值，而是通过Sigmoid函数把速度映射成“取1的概率”：

v = wv + c1r1*(pbest - x) + c2r2(gbest - x) S(v) = 1 / (1 + exp(-v)) if rand < S(v), then x = 1, else x = 0

这里的x是0/1向量，每一位代表一条可操作支路的开关状态。速度v越大，S(v)越接近1，粒子更倾向选择闭合；v越负，则更倾向打开。标准BPSO实现简单，但在配电网重构这种离散组合优化问题上容易出两个毛病。一是速度更新后期粒子位置变化很小，容易早熟；二是随机生成的二进制串很难满足辐射状约束，大量粒子无效，搜索效率低。所以核心论文里的“改进”基本都围绕这两个痛点展开，要么改惯性权重，要么加变异，要么改编码方式。

3.2 我做了哪几个改进

我复现的论文里，改进BPSO大致包含四个动作。

第一，惯性权重非线性递减。标准BPSO的w是常数或线性递减，但配电网重构解空间存在大量不可行区域，前期需要强探索能力去发现不同连通区域，后期需要细致搜索。我用的是：

w = w_max - (w_max - w_min) * (t / T)^2

在迭代初期w接近0.9，全局搜索能力强；后期w接近0.4，局部搜索更精细。相比线性递减，这种凹形曲线让前期保持更久的探索时间，收敛曲线也更平滑。

第二，加入变异重置算子。当全局最优连续多代没有更新时，对部分粒子按变异概率随机翻转若干位，同时重置这些粒子的历史最优位置。这本质上就是在BPSO里引入遗传算法的突变思想，让陷入停滞的种群重新活跃起来。方差太小的种群，就算粒子位置还有变化，也只是在小邻域里反复试探，这时候变异比调参数更直接。

第三，用环路编码代替全支路编码。标准BPSO如果直接用37位0/1表示所有支路，随机产生恰好32条闭合支路的概率极低，几乎不可能得到可行解。我复现时把37条支路按5个基本环路分组，每个粒子只维护5个整数，分别表示5个环路里打开哪条支路。这样“打开支路总数=5”变成隐式约束，天然满足闭合支路数为32。剩下只需要检查连通性，不可行解数量少了一个数量级。

第四，精英保留策略。每代在更新位置前，先把当前全局最优解保存一份，种群更新后如果新解没有超过它，就用保存的精英直接替换种群中最差个体。这一步成本极低，但能保证算法在非常恶劣的初始种群下也不容易丢失全局最优。很多工程复现忽略这一点，导致改进算法在几次随机初始化下结果波动很大，其实不一定是算法问题，是缺少精英保护。

3.3 改进后的算法循环体

改进BPSO的整体流程可以简化为：

1. 初始化种群。根据环路编码，每个粒子生成一组打开的支路开关，再映射成37位开关状态，检查连通性并计算初始适应度。
2. 把每个粒子的当前位置记录为个体最优pbest，全局最优位置记作gbest。
3. 按非线性公式更新惯性权重w。
4. 对每个粒子更新速度v，并对速度做限幅处理。
5. 用Sigmoid函数把v映射成概率，再随机决定每个环路内打开开关的更新方向。这里需要注意：环路编码里每个粒子是5个整数，更新时要对每个整数位做二进制化，不是简单加1减1。
6. 对更新后的粒子做变异重置，判断是否改善。
7. 计算适应度，更新pbest和gbest。
8. 精英替换差个体。
9. 判断是否达到最大迭代次数或连续未更新代数，若满足则输出gbest对应的开关组合和潮流结果。

整个流程看起来不复杂，但真正写代码时，步骤5和步骤6是最容易出bug的地方。尤其环路编码在做速度更新时，不能简单把5个整数当连续变量直接取整，否则很多粒子会生成相同的开关组合，种群多样性快速下降。我用的办法是把“每个环里选择哪条支路打开”拆成几个二进制的“热编码”，用BPSO的标准更新方式处理，然后再解码回整数。具体代码逻辑我在下一节展开。

4. Matlab代码实现关键细节

4.1 编码和解码：最容易翻车的一步

先说清楚IEEE33节点的支路构成。33个节点之间共有37条支路，其中前32条是分段开关支路，第33到第37条是联络开关支路。原始运行状态是分段开关全部闭合、联络开关全部打开，形成一个辐射状网络。重构时我们需要打开若干分段开关、闭合对应的联络开关，最终仍然满足33个节点32条闭合支路的辐射树。

如果你直接用一个长度为37的0/1向量，1表示闭合，那初始种群随机生成时，大多数粒子闭合支路数不是32，根本过不了拓扑约束，这也是我在6.2节会展开讲的“全是不可行解”坑。改进版的复现中，我强烈建议采用“5个环路编码”：把5条联络开关对应的基本环分别找出来，每个环内只能打开一条支路，粒子用5个整数表示这5个打开位置。这样打开支路总数固定为5，闭合支路数自动等于32，问题规模也从2^37缩小到每个环内可选支路数的乘积，搜索空间大幅减小。

但环路编码也有代价：它要求你先把每个环包含哪些支路列准确。IEEE33节点的标准环路信息可以从前推回代需要的网络拓扑中得到，也可以参考文献直接整理成常量表。实际复现时，我会用Matlab的结构体保存每个环路的支路编号，比如loop{1}=[7 8 33 20 21 19 18]，loop{2}=[9 10 34 21 22 20]，不同论文的环路划分可能有差异，但只要保证5个环覆盖所有可打开支路，并且任意一个粒子在每个环中选一条打开后，网络有可能连通，就可以用于搜索。解码时先得到打开的5条支路，再和37位全闭合状态取差集，得到最终开关状态。

4.2 潮流计算函数：前推回代法的Matlab实现

前推回代虽然原理简单，但Matlab代码的坑在于支路数据要按“层次”排好，或者每次迭代都从末端往前扫描。IEEE33节点结构固定，最简单的是用邻接表+后序遍历，也可以直接用节点支路关联矩阵做前推。下面是我复现时常用的核心片段，去掉了完整节点数据，只保留逻辑：

function [V, Ploss, converged] = backward_forward_33(branch_ij, loadNode) % branch_ij: 闭合支路列表，每行[首端节点, 末端节点, R, X] % loadNode: 每个节点的有功无功负荷 nNode = 33; V = ones(nNode, 1); % 初始电压标幺值 Sload = loadNode(:,1) + 1j*loadNode(:,2); % 标幺值 for iter = 1:100 Vold = V; % 前推：从末端向首端累加注入功率 branchS = zeros(size(branch_ij,1), 1); for i = size(branch_ij,1):-1:1 n1 = branch_ij(i,1); n2 = branch_ij(i,2); % 需要根据节点深度判断n2是否为末端节点 % 这里省略了子支路功率累加过程 branchS(i) = conj(Sload(n2)/conj(V(n2))); end % 回代：从首端向末端更新电压 for i = 1:size(branch_ij,1) n1 = branch_ij(i,1); n2 = branch_ij(i,2); Sij = branchS(i); Iij = conj(Sij / V(n1)); V(n2) = V(n1) - (branch_ij(i,3) + 1j*branch_ij(i,4)) * Iij; end if max(abs(V - Vold)) < 1e-6 break; end end Ploss = 0; for i = 1:size(branch_ij,1) Iij = conj(branchS(i) / V(branch_ij(i,1))); Ploss = Ploss + abs(Iij)^2 * branch_ij(i,3); end converged = (iter < 100); end

这段代码表现了前推回代的核心，但直接照抄不一定能跑通，因为你必须把节点分层信息和“子支路累计功率”处理对。我的经验是先用节点深度（从节点1开始的支路跳数）做排序，确定每个节点的子节点集合，再按深度从大到小前推、从小到大回代。如果只针对IEEE33节点，甚至可以手工把所有支路按末端节点编号映射到前推顺序，代码写死也能接受。关键点是所有功率、电压都用标幺值，节点负荷的单位别和基准容量弄混。

4.3 主程序框架与参数设置

主程序需要做的事情很明确：加载网络参数、设置算法参数、初始化种群、循环迭代、输出结果。下面是一个骨架：

% IEEE33节点标准数据，branch = [支路编号 首节点 末节点 R X] % loadData = [节点编号 P(kW) Q(kvar)] rng(2024); N = 40; % 种群规模 T = 100; % 最大迭代次数 wMax = 0.9; wMin = 0.4; c1 = 1.5; c2 = 1.5; vMax = 6; % 速度上限，控制位置翻转概率 mutProb = 0.1; % 变异概率 % 初始化5环路编码的位置（每个环路打开一条支路） % 这里简化为直接随机选中每条环路的支路编号 pop = zeros(N, 5); for k = 1:N for loopIdx = 1:5 pop(k, loopIdx) = loopList{loopIdx}(randi(length(loopList{loopIdx}))); end end % 主循环省略，具体见流程 for t = 1:T w = wMax - (wMax - wMin) * (t/T)^2; for k = 1:N % 更新速度、sigmoid映射、解码到37位开关状态 % 计算适应度 end % 更新个体最优和全局最优 % 变异重置、精英替换 end

参数方面，我建议种群规模N取30到60，太大收敛慢，IEEE33节点没必要盲目扩大种群；迭代次数80到150次足够。速度上限vMax不能太小，否则Sigmoid函数中间区域概率变化不明显，粒子容易一直集中在某些位置上；也不能太大，否则取1概率长期接近0或1，容易失去二进制翻转能力。vMax取4到6是比较稳的区间。

4.4 适应度函数怎么写

适应度函数是整个代码的“翻译官”。它输入一个粒子的开关状态，输出一个标量适应度。我的实现顺序是先检查是否满足连通性和闭合支路数，再调用潮流计算，最后计算网损和惩罚。

function fit = evalFitness(switchState, branchData, loadData) % switchState: 37维0/1，1表示闭合 closedIdx = find(switchState == 1); if length(closedIdx) ~= 32 || ~isRadial(closedIdx) fit = 1e6; % 不可行解给大惩罚 return; end % 根据closedIdx生成闭合支路矩阵 branch_ij = branchData(closedIdx, :); [V, Ploss, converged] = backward_forward_33(branch_ij, loadData); if ~converged fit = 1e6; return; end % 加入电压越限惩罚 vPenalty = sum(max(0.9 - V, 0))*1000 + sum(max(V - 1.05, 0))*1000; fit = Ploss + vPenalty; end

有一点要提醒：连通性判断不能只靠闭合支路数等于32。IEEE33节点上，可能出现一个粒子恰好有32条闭合支路但支路分成两个连通块的情况，比如节点1到10孤立成一个小树，节点11到33形成另一个树。这种情况潮流计算无法进行，必须直接用深度优先搜索或Matlab的graph函数检查根节点是否能到达所有节点。我在调试阶段吃过这个亏，当时觉得“支路数对了肯定没问题”，结果一跑潮流全是NaN。

5. 实验结果与复现心得

5.1 IEEE33节点基准算例的典型结果

复现完成以后，我的结果和主流文献基本一致。原始状态是32条分段开关全部闭合、5条联络开关全部打开，网损约202.68kW，最低电压约0.9131pu。经过改进BPSO重构后，找到的最优开关组合为：打开7-8、9-10、14-15、28-29、32-33这5条分段支路，同时闭合5条联络支路8-21、9-15、12-22、18-33、25-29。这个组合下网损约139.55kW，最低电压约0.9378pu，降损幅度约31.2%。

我把结果整理成对照表：

需要说明的是，不同文献给出的最优网损可能在139.4到139.6kW之间波动，这很正常，因为编码方式、潮流公式、电压约束的惩罚系数以及开关环路划分都可能造成微小差异。复现时只要看到收敛后的网损落在该区间，开关组合也能对应到合理的辐射树，基本可以判定程序正确。

5.2 改进算法的收敛表现

我用标准BPSO和改进BPSO各跑了20次，统计结果如下：

改进BPSO的平均收敛代数明显更早，稳定性也更好。标准BPSO经常在迭代后期陷入局部解，尤其是当种群里的粒子都汇聚到同一个不可行区域附近的连续空间时，很难再跳出来。变异重置算子在这里起到了关键作用，它让粒子在连续几代没有更新后重新开始探索，而不是在同一个坑里反复打转。

不过也要强调，收敛代数提前并不完全代表算法更好，还需要看最终解质量和多次运行稳定性。有些改进算法加了强局部搜索，收敛很快，但容易把搜索过程变成“爬山”，本质上是用局部搜索替代全局搜索，对更大规模网络未必有效。所以我在复现时并不是只看收敛曲线，会把20次运行的最优解分布、最差值、平均值全部记录下来。

5.3 复现时最容易误判的“陷阱”

第一是单位混乱。IEEE33节点典型基准容量是100MVA，基准电压12.66kV，负荷数据常用kW和kvar。如果直接将kW除以100后当作标幺值，数值倒是没错，但电压、阻抗和功率之间的配合容易出错。我建议在写数据文件时就完成标幺转换，无论是节点负荷还是支路阻抗，统一用标幺值。

第二是网损的单位容易写错。前推回代返回的Ploss是基于标幺功率的结果，如果基准容量是100MVA，0.0013955标幺网损对应139.55kW。看到这个数字不要激动，也不要直接把0.001395当kW写进论文。要在输出函数里乘Sbase。

第三是开关编号错位。支路编号和节点编号不是一回事。比如支路7-8可能对应支路编号7，但联络开关8-21对应支路编号33。编码矩阵里的索引一旦错位，结果可能是一个看起来很合理但完全无法物理实现的开关组合。

6. 常见问题与排查技巧实录

6.1 潮流计算不收敛怎么办

如果前推回代函数返回不收敛，我一般按顺序排查四件事。第一，闭合支路数是否为32；第二，网络是否连通。很多情况下不收敛不是因为方程问题，而是拓扑根本不可行，前推回代在遇到环路或孤岛时会计算出错误功率。第三，节点负荷是放在节点上还是支路上，有没有把支路充电功率也混进去；IEEE33节点是三相平衡模型，没有充电电容，但如果换其他算例就可能出错。第四，迭代初值是否全为1.0pu。只要拓扑正确，前推回代在辐射状配电网中很少真不收敛，如果反复出现NaN，多半是除零，而在配电潮流里除零通常是因为某个节点被孤立，电压为0。

还有一种情况是算法搜索到的树虽然连通，但某些支路方向判定错了。前推回代需要知道电流从根节点到末端的方向，如果支路方向列表和节点编号不符，回代时电压会算成负值或复数。遇到这种现象，先把最小电压打印出来，如果出现负电压，基本就是支路方向反了。

6.2 种群全是不可行解

我最早用全37位二进制编码跑标准BPSO时，出现了一个让人崩溃的现象：前几十代所有粒子适应度都是1e6，也就是说一个可行解都没搜到。原因很好理解：从37条支路里恰好选中32条闭合支路的概率是C(37,32)/2^37，而C(37,32)=C(37,5)=435897，分母是137438953472，算下来只有约0.0003%，这还没算连通性。所以随机二进制编码的初始种群，绝大多数粒子闭合支路数都不对。

解决办法有两个。一是换成环路编码，我在4.1节说过，把每条基本环路看成一个“候选槽”，粒子只需要在每个槽里选一条支路打开，打开支路总数自动等于5。二是继续用37位编码，但要增加修复算子：初始化时先随机选出5条支路打开，其它支路置为1，保证闭合支路数为32；后续每次粒子更新后，同样检查1的数量，如果多了就随机把一部分1置0，少了就随机把一部分0置1。修复以后再用DFS检查连通性。在实际复现中，我更推荐修复算子，因为它的通用性强，后面换IEEE123节点时不需要重新定义环路。

6.3 Matlab代码的细节坑

有几个Matlab层面的坑值得单独说。第一个是老生常谈的rand和randi混用。BPSO更新位置时，rand是用来和概率值比较的，randi是用来随机选择整数支路编号的。如果你忘了初始化随机种子，每次跑结果都不一样，复现论文时会很难判断参数好坏。建议在main函数开头固定rng(2024)，需要统计多次运行结果时再在循环里重新设置不同种子。

第二个坑是矩阵运算维度不匹配。IEEE33节点是33个节点、37条支路，有些数据文件把节点编号从0开始，也有文件把支路编号从0开始，导致索引偏移一位。我的习惯是节点1作为电源根节点，支路编号从1到37，所有函数统一用这个约定。

第三个坑是保存结果时直接把复数功率写入表格。如果V和I都是用复数计算，Ploss是实数，但中间变量branchS是复数。输出时要用real()和abs()转换，否则命令窗口里一切正常，导出Excel时出现“x + yi”格式，后面处理数据又得重新清洗。

7. 从复现走向改进

7.1 后续扩展方向

IEEE33节点上的改进BPSO复现只是第一步，很多人会把结果继续往不同方向扩展。一是多目标化，把网损最小、电压偏差最小、开关动作次数最少三个目标放进同一个框架，用多目标BPSO求Pareto解集，再用模糊隶属度选折中解。二是分布式电源接入，风机、光伏接入后会改变潮流分布和辐射状约束的边界，算法需要处理DG出力的不确定性，常见做法是场景生成加机会约束规划。三是动态重构，负荷随时间变化时，不能只做一个时间断面的开关组合，需要在时序维度上优化开关动作次数，改进BPSO这时候通常要和动态规划或滑窗策略结合。四是算法对比，多跑几种群智能算法，比如灰狼优化、樽海鞘群、鲸鱼优化、遗传算法，在IEEE33和IEEE69节点上做统计对比。

这些扩展的方向都建立在基础算例和潮流程序正确的前提上。如果你已经把IEEE33节点的潮流、BPSO编码、适应度函数和收敛统计全部调通，换成其他算法时只需要替换种群更新部分，潮流和评价函数完全复用。这也是复现核心论文带来的最大收益：你不仅复现了一个算法，还搭出了一个可复用的配电网重构测试平台。

7.2 几条经验

最后说几条个人体会。第一，改进算法一定要和标准算法对比，且对比时参数规模要一致。不能改进版用40个粒子跑150代，标准版用20个粒子跑50代，这种对比没有意义。第二，多运行几次，不要拿单次运行结果写结论。群智能算法本质是随机优化，一次运行得到139.55kW可能是运气好，也可能只是碰巧。我看论文时特别关注作者是否给出多次运行的标准差，如果只有一条漂亮收敛曲线，可信度要打个问号。第三，保存好每一次实验的随机种子、参数、开关组合，方便回溯。复现过程中我有一半时间不是在改算法，而是在找“上次那个结果是怎么跑出来的”。有完整的实验记录，这些问题都能避免。

我现在拿到任何配电网重构论文，第一件事就是把IEEE33的基准复现出来，再套论文里的算法。这个流程看起来很慢，实际是最省时间的：算例数据对了、潮流函数对了，后面的算法改进才站得住。希望这篇记录能帮你少踩几个坑。