基于粒子群算法的IEEE 33节点配电网重构实现详解

前阵子我在整理配电网优化方向的小项目目标很明确不新增线路、不装无功补偿只靠改变网络里的开关开合状态把运行损耗降下来。这个问题学名叫配电网重构。我选的标准对象是IEEE 33节点测试系统算法用了粒子群算法PSO整套程序用MATLAB实现。从问题建模、编码设计、潮流计算到调参排坑我完整走了一遍今天把这些探索过程整理出来重点讲清楚每一步为什么这么做尤其是一些容易被论文一笔带过、但实际操作时非常容易卡住的细节。如果你正准备用智能算法做配电网重构或者想找一个带潮流的优化算例来练手这篇内容可以直接当作参考路线。1. 重构到底在优化什么目标、约束与可行解1.1 配电网重构的工程本质配电网正常运行时为了保证系统是辐射状的分段开关基本全部闭合联络开关全部打开。这时候电流只能沿着固定的树状路径从变电站流向各个负荷节点。问题在于树状结构决定了潮流分布而有些树状结构会让电能绕远路导致线损偏高、末端电压偏低。重构的核心思路就是改变那些分段开关和联络开关的组合状态把一部分负荷转移到另一条供电路径上让整个网络在仍保持辐射状的前提下实现更优的运行状态。这个思路有点像城市交通的潮汐车道道路数量没变但通过改变通行方向让拥堵路段的车辆分流到空闲路段。配电网里的“车辆”是电流“道路方向”由开关状态决定。重构的收益通常体现在几个方面降低有功网损、改善电压分布、均衡馈线负载以及故障后恢复非故障区供电。科研和工程里最常见的单目标重构就是最小化系统有功损耗。1.2 为什么大家都在用IEEE 33节点“练手”33节点系统来自IEEE标准测试系统库最早作为配电网重构和分布式电源接入研究的公开算例广泛流传。它的规模对算法验证来说非常合适节点数不多不少潮流计算秒级完成又能体现出拓扑变化带来的明显差异结构也足够复杂能验证辐射状约束处理得对不对。系统的基本参数如下不同文献里的节点编号习惯会略有差别有些从0开始编号有些从1开始编号但物理拓扑一致项目典型数值节点数量33含1个电源节点系统基准电压12.66 kV分段开关数量32联络开关数量5系统总有功负荷约3715 kW系统总无功负荷约2300 kvar初始运行网损约202.68 kW不同版本略有差异默认状态下5条联络开关全部打开网络是典型的辐射状单电源系统。把5条联络开关全部合上后网络会形成5个独立的基本环路这5个环路是后面设计粒子编码时最关键的几何结构。标准算例重构后网损通常可以降到140 kW附近降损幅度非常可观这也是初学重构时很有成就感的“验证锚点”。2. 算法选型为什么选了粒子群算法2.1 重构在数学上到底有多难重构的决策变量是开关状态每个开关只有闭合和断开两种状态。如果直接看全系统37个开关就意味着2的37次方种组合规模达到千亿级别。当然实际工程中为了维持辐射状结构不会任意开关组合都行真正有效的候选解远远小于这个数字。但即便我们只考虑闭合几条联络开关并做等量分段开关切换组合数仍然很大而且目标函数无法写成决策变量的显式表达式必须靠潮流计算来得到损耗值。从优化视角来看这是一个典型的离散、非线性、非凸、带强约束的组合优化问题。传统基于梯度的优化方法在这里基本失效因为目标函数对“开关状态”不可导解空间还不连续。工程和学术界的普遍做法是采用元启发式智能算法例如遗传算法、粒子群算法、模拟退火、蚁群算法等。它们不要求梯度信息只要能把候选解映射成一个可评估的适应度值就能搜索。2.2 粒子群与遗传算法的选择考量我早期做这类问题时也纠结过到底用遗传算法还是粒子群算法。遗传算法的思路是通过二进制或十进制编码模拟染色体依靠选择、交叉、变异迭代进化。它的全局搜索能力不错但参数多、收敛速度相对慢而且交叉算子在处理强约束问题时经常生成大量非法个体修起来很麻烦。粒子群算法的机制则完全不同。它模拟鸟群觅食行为每个粒子作为一个候选解在解空间里依靠自身历史最优和群体历史最优不断调整飞行方向和速度。和遗传算法相比粒子群的主要优势是参数少实现简单核心逻辑一套公式就写完前期收敛速度快惯性权重的引入又提供了局部开发和全局探索的平衡能力在处理33节点这种中等规模、维度不高的问题时有天然的高效性。尤其是我们用“基本环路编码”之后粒子维度只有5维PSO的搜索优势会非常明显。如果切换到几百节点的网络PSO同样能跑只是需要针对维度和约束做更多改进。2.3 从连续版本到离散开关决策理解粒子群要把握最核心的速度-位置更新公式v w * v c1 * r1 * (pbest - x) c2 * r2 * (gbest - x)x x v其中w是惯性权重c1和c2是学习因子r1、r2是[0,1]的随机数pbest是粒子自身历史最优位置gbest是整个种群的历史最优位置。式子的物理含义是粒子下一时刻的速度由三部分决定一是维持当前运动方向的惯性二是向自身最优位置学习的认知部分三是向群体最优位置学习的社会部分。但关键问题来了上面的公式工作在连续实数空间中而重构的决策变量是离散的开关编号。如果直接把连续粒子位置强行取整容易产生大量重复或越界的选择后期收敛也会出现震荡。更稳妥的处理方式不是推翻PSO而是把连续变量重新解释为“在某个可选开关集合里的序号”然后通过round取整、越界映射等策略把它转换成可用的整数决策。这部分我在后面第4章会重点展开。3. MATLAB程序框架从原始参数到潮流计算3.1 模块划分与运行环境这个程序不需要安装任何第三方工具箱纯MATLAB脚本就可以运行。我在MATLAB R2022b上测试过基础矩阵运算加几个自定义函数就能跑完整个流程。程序模块我按下面的方式拆分保证每个文件职责单一文件/函数作用main.m主程序完成参数读取、编码初始化、PSO迭代、结果输出calc_powerflow.m前推回代法潮流计算输入开关状态输出节点电压和支路损耗decode_switch.m将粒子向量解码为具体的支路开合状态get_basic_loops.m计算5条联络开关对应的基本环路check_connectivity.m检查网络是否辐射状且全部节点连通objective_func.m调用潮流计算计算目标函数网损 电压越限惩罚文件拆开的直接好处是排错容易。比如怀疑潮流算错了可以单独传一个初始开关状态进去看潮流结果怀疑编码有问题则可以单独检查解码后的开关列表是否合法。3.2 33节点数据怎么组织最顺手很多初学者会把节点负荷和支路阻抗散乱地写在主脚本里一旦要切换算例就要改一堆地方非常容易出错。我的做法是用两个矩阵集中管理数据和程序逻辑分离。支路数据矩阵每行代表一条支路列依次为首端节点编号、末端节点编号、电阻欧姆、电抗欧姆。前面32行是分段开关对应的线路最后5行是联络开关对应的线路。负荷数据矩阵存节点有功和无功单位是kW和kvar。需要注意的是业界流传的IEEE 33节点数据存在0基和1基两套编号版本。我建议在程序里统一使用1基编号节点1是电源节点节点2到节点33是负荷节点这样所有矩阵下标不用做偏移读起来直观。另外潮流计算里我推荐把所有参数归算到标幺值。取基准容量10 MVA、基准电压12.66 kV那么基准阻抗就是Z_base 12.66^2 / 10 16.0276 Ω把原始Ω数据除以基准阻抗得到标幺阻抗把负荷除以基准容量得到标幺功率。标幺处理的好处是潮流计算中电压基本在1.0附近收敛判据设置更直观迭代公式也更简洁。3.3 辐射状配电网为什么要用前推回代法输电网潮流常用牛顿-拉夫逊法因为它收敛快、适用于环网。但配电网有自己的特点辐射状结构、支路数等于节点数减一、线路R/X比值相对较大。直接用牛顿法时雅可比矩阵条件数容易变差甚至不收敛需要各种改进策略。对于辐射状网络前推回代法显然更简单可靠每次迭代只有两个步骤。第一步“前推”假设各节点电压初值已知通常取额定电压1.0从网络末端开始根据节点负荷和下游支路功率逐条支路向上游推算支路功率。第二步“回代”从电源节点开始已知首端电压根据已经求出的支路功率逐条支路向下游推算各节点电压。两步交替迭代直到两次迭代的电压差满足精度要求比如小于1e-6。写成MATLAB核心逻辑大致如下实际运行时需要先根据开关状态生成网络的拓扑顺序保证“前推”是从离电源最远的节点逐层向回走“回代”是从电源逐层向外推% 输入: branchR/branchX 标幺阻抗, Pnode/Qnode 节点注入功率, % order 节点遍历顺序(从末端到根节点的前推顺序) % 输出: V 节点电压向量, Ploss 总有功损耗 V ones(nnode, 1); % 电压初值 S zeros(nnode, 1); % 各支路末端复功率 Snode Pnode 1j * Qnode; % 节点负荷 for iter 1:200 V_old V; % 前推根据电压和下游功率求支路功率 S_branch zeros(nbranch, 1); for k 1:nbranch % 这里用父节点/子节点关系代替从末层往根遍历 % 支路流入功率 下游节点负荷 下游支路损耗 自身损耗 % 为保持简洁具体索引映射省略 end % 回代从根节点向下更新电压 for k 1:nbranch V(child) V(parent) - (conj(S_branch(k)) / conj(V(parent))) * Z(k); end if max(abs(V - V_old)) 1e-6 break; end end Ploss sum(real(S_branch .* conj(S_branch)) .* R_branch);这种分层遍历如果用递归去写会很慢而且容易栈溢出。更工程化的做法是先用广度优先或深度优先搜索建立父子关系表再按照父子关系表从末端到根节点做前推从根节点到末端做回代整个过程就是循环加索引查表非常快。3.4 初学最容易忽略的潮流细节前推回代公式本身不难真正的坑在“节点注入功率的符号”和“支路功率的方向一致性”上。配电网负荷通常看成恒功率模型即PQ节点吸收功率为正电源节点发出功率为负。在前推过程中某支路末端的累积功率应该等于该支路下游所有节点负荷与支路损耗的总和这一项是累加出来的不能用错。回代计算电压时要用该支路的首端电压减去压降压降等于支路电流乘以线路阻抗注意电流要用共轭功率来计算少写一个conj就会出现方向错误。程序里建议把所有的功率相关变量统一用复数表示单位用标幺值这样可以显著减少手工推公式时的符号错误。4. 粒子群重构的四个关键环节4.1 找到5个基本环路是编码的第一步当前面提到“粒子维度只有5维”时很多人的第一反应是5维对应5条联络开关吗如果直接把每条联络开关的开合状态作为一维那等于默认每次只能投切联络开关无法解释分段开关怎么联动。重构的真实决策应该是在保持辐射状的前提下同时改变若干分段开关和联络开关的状态。所以正确的做法是先构建基本环路。方法是把5条联络开关全部合上此时网络从一棵树变成了带5个基本环的连通图。在图论里树的边数为节点数减一每增加一条联络支路就恰好增加一个独立回路。因此5条联络支路对应5个基本环路。每个环路里除了那条联络开关本身其余支路都是初始状态闭合的分段线路。我的MATLAB实现步骤是这样的用基础网络全部联络断开构建邻接表对每条联络支路找到它两端节点在基础网络中的唯一通路把这条通路上的所有普通分段支路收集起来和该联络支路一起组成一个环路集合粒子第k维就决定在第k个环路里选择哪一条分段支路断开对应第k条联络支路强制闭合。由于基础网络是树任意两个节点之间的通路是唯一的所以环路搜索不需要复杂的环基算法用BFS或DFS从联络支路的一端搜到另一端即可。伪代码如下function loop_branches find_path_bfs(adj, startNode, endNode, branchList) % 用BFS在基础树中找 startNode 到 endNode 的路径 % 返回路径上经过的支路编号 parent zeros(size(adj,1),1); visited false(size(adj,1),1); queue startNode; visited(startNode) true; while ~isempty(queue) cur queue(1); queue(1) []; if cur endNode break; end neighbors find(adj(cur,:)); for nb neighbors if ~visited(nb) visited(nb) true; parent(nb) cur; queue(end1) nb; end end end % 回溯得到节点序列再映射为支路编号 loop_branches branchList(path_edges); end对33节点系统来说这5个基本环路在拓扑上是相互独立的同一个普通分段支路不会同时属于两个环路。这一点保证了后面的编码不会出现“两条环路同时选中同一条分段开关”的冲突是编码能用低维整数向量的重要前提。4.2 粒子结构设计与合法开关解码我采用的粒子结构如下x [x1, x2, x3, x4, x5]其中xk表示第k个环路中选中的分段支路在可选集合里的序号。例如第1个环路里包含的可断分段支路有5条那么x1只能在1到5之间取整数。解码时先取出第k个环路对应的分支编号列表然后按xk索引选择要打开的分段支路。与此同时第k条联络支路的开关状态设置为闭合。最终闭合的支路数始终等于32原本的32条分段支路减去5条被打开的分段支路再加上5条被闭合的联络支路刚好满足33节点树的边数要求。这里有个特别值得强调的坑解码时粒子每一维能选择的支路必须是“环内的普通分段支路”不能直接选联络支路。有些教程图省事把所有37条支路都列进可选集合结果显示联络支路闭合后打开的却是联络支路本身等于联络支路一直处于打开状态闭合支路数变化就乱了输出的结构根本不是期望的辐射状网络。粒子更新的速度与位置公式沿用标准PSO但每一维位置更新后必须做整数化处理。我的做法是先用round取整再把越界的索引拉回到该环可选区间的边界。虽然这种方式理论上会降低一点搜索精度但在33节点这个量级上效果良好而且远比复杂的离散PSO变种容易实现。更新约束如下x_new(i) x(i) v_new(i); x_new(i) round(x_new(i)); x_new(i) max(min(x_new(i), length(loop(i))), 1);4.3 辐射状与连通性校验不能省只靠环路编码还不能百分之百保证候选网络一定合法。虽然理论上每个环路断开一条普通分段支路后系统闭合支路数恰好为32且基本环路相互独立但为了保险起见“辐射状 全连通”的校验必须保留。这个校验函数的数据结构很简单从电源节点出发做BFS记录能访问到的节点数如果访问数等于全部33个节点且闭合支路总数等于32则网络判定为可行。function valid check_valid_topology(branchState, nnode) % branchState 是37维逻辑向量1表示闭合0表示断开 if sum(branchState) ~ nnode - 1 valid false; return; end adj build_adjacency(branchState); visited bfs_from_source(adj, 1); valid (length(visited) nnode); end为什么已经用了环路编码还要再做一次校验因为代码在手写过程中容易出边角问题比如环路可选集合生成错误或者粒子解码时索引越界导致某些支路状态漏改。在目标函数入口处统一做一次校验能第一时间暴露编码函数的逻辑错误避免把非法网络送进潮流计算得到莫名其妙的损耗值。4.4 目标函数里的电压惩罚怎么加重构的目标函数通常定义为有功网损最小。但如果只追求网损可能会选出一些让末端电压偏低甚至越限的网络结构因此目标函数里要加入电压约束惩罚项。我使用的适应度表达式为fitness Ploss lambda * sum(max(0, Vmin - V_i)^2)其中Vmin取0.95标幺值lambda是惩罚系数通常取1000到10000之间的值。具体取多大要和网损量纲匹配如果网损是几十千瓦到两百千瓦的有名值惩罚系数就要足够大让越限解在粒子群竞争中快速被淘汰。若用标幺值计算则一般取值1000已经能起到明显效果。目标函数的调用流程如下function fitness objective_func(x, loop, branchParam, nodeParam) openSet decode_switch(x, loop); branchState get_state_from_openSet(openSet); if ~check_valid_topology(branchState, 33) fitness 1e8; % 非法拓扑直接给大罚值 return; end [Ploss, V] calc_powerflow(branchState, branchParam, nodeParam); penalty 1000 * sum(max(0, 0.95 - V).^2); fitness Ploss penalty; end4.5 PSO主循环参数与完整流程主程序的循环框架相对固定。先做参数初始化然后用随机整数初始化粒子的位置和速度迭代里依次做目标函数评估、更新pbest和gbest、更新速度和位置。推荐的参数组合如下种群规模30到60最大迭代次数80到150惯性权重w从0.9线性递减到0.4学习因子c1和c2均取2.0粒子速度限制每维最大速度取该环可选集合长度的三分之一左右。惯性权重从大到小递减是我强烈建议不要省掉的操作。w大时粒子倾向于全局大范围探索w小时粒子更多在局部精细开发。如果用一个固定的小w粒子很容易在前期就全部聚集到某个局部最优附近后面再想跳出来就很难了。5. 常见问题与排坑实录5.1 前推回代一直不收敛怎么办跑PSO时如果频繁出现潮流不收敛多半不是算法问题而是候选网络结构本身就是非法的。比如解码环节把联络开关误处理成打开状态导致某些节点不在任何一个完整供电路径上或者环路中选错了断开支路导致网络出现环网。这时候要在目标函数入口增加拓扑校验并直接惩罚而不是让潮流计算硬跑。我踩过一次坑没做校验直接让前推回代在环网上迭代结果电压震荡完全不收敛程序卡在目标函数里出不来排查