先说明一点做配电网优化方向的同学不管你是写毕业论文、发小论文还是做横向课题肯定绕不开“无功优化”和“有功调度”这两个老问题。但如果你把两者分开做往往得到一个看似收敛、实际顾此失彼的结果。这篇博文记录的就是我把两者真正“揉”在一起的完整过程——基于改进多目标粒子群算法小生境粒子群的配电网有功-无功协调优化配全套Matlab代码实现方法。我会从问题建模、算法改进思路、代码架构、IEEE 33节点系统实测结果到调参和踩坑一次性讲透方便你直接“抄作业”。1. 项目概述与问题定义1.1 为什么要把有功和无功放在一起优化传统配电网运行调度里有功调度和无功调节通常是分开做的。有功调度主要盯发电机的出力分配目标是经济性无功调节主要管无功补偿设备电容器组、SVC、SVG、OLTC等目标是电压合格率和降低网损。表面上看这没什么问题但在分布式电源DG大量接入之后情况完全变了。DG的出力是波动的甚至有反调峰特性它同时改变配电网的有功潮流分布和无功电压分布。如果你只做无功优化假定有功出力是固定的那在DG出力大幅变化时电压可能越限无功补偿设备动作频繁。如果你只做有功调度忽略无功和电压约束那优化出来的DG出力方案可能根本没法落地因为节点电压可能已经超出允许范围。换句话说有功和无功在配电网里是强耦合的。把它们分开优化本质上是在做一个“次优解”——这也是大多数基于PSO的配电网优化论文只做单一目标时仿真结果看起来挺好、但实际工程没法用的根本原因。把两者统一到一个多目标优化框架里让算法同时搜索有功出力和无功补偿的可行解空间这才能得到真正可执行的调度策略。1.2 有源配电网协调优化的两个核心目标我做这个项目时目标函数选了两个不贪多因为目标越多Pareto前沿的搜索难度就越大而且工程上最关注的其实就这两个第一个目标是系统网损最小。配电网的网损包括有功损耗和无功损耗但在目标函数里我们只关心有功损耗。无功功率在传输过程中不直接产生电度损耗但它会占据线路容量、增大电流从而间接增大有功损耗。所以降低网损的实质是在满足电压约束的前提下尽量减少无功功率在配电网中的长距离流动同时优化DG有功出力减少倒送功率。第二个目标是电压偏移最小。电压质量是配电网运行的生命线。分布式电源接入后最常见的电压问题就是末端电压抬升和局部节点电压越限。我用所有负荷节点电压与额定电压的偏差绝对值之和或者最大偏差来衡量这个目标数值越小代表电压质量越好。当然你也可以根据需要加上第三个目标——比如DG运行成本最低、开关动作次数最少、碳排放最小等代码架构上完全可以扩展后面会讲。1.3 配电网模型与决策变量的选取先说测试系统我用的是最经典的IEEE 33节点配电网系统。为什么选它因为它的参数公开、文献对照多、算例结果可查而且网络结构是一个典型的辐射状配电网包含联络开关负荷分布不均匀电压问题突出适合用来验证优化算法的有效性。关于决策变量我分成了两类第一类是有功调度变量就是各分布式电源的有功出力。假设系统里接了若干个DG比如在节点18、22、25、33接入光伏或风机那么每个DG的有功出力就是一个决策变量。需要注意的是DG有功出力不是无限可调的要设置上下限通常取其额定容量的0%到100%或者根据预测出力设置一个可调范围。第二类是无功控制变量包括无功补偿装置的投切容量和DG的无功出力。如果用离散的电容器组那么决策变量是整数代表投切的组数如果是SVG或DG逆变器可以连续调节。DG逆变器参与无功调节时要满足视在功率约束——有功和无功不能同时达到最大值即[ Q_{DG} \le \sqrt{S_{DG}^2 - P_{DG}^2} ]这个约束在粒子群算法里处理起来比较麻烦因为它是一个动态的、与另一个决策变量耦合的约束。我采用的是罚函数法加可行性修正相结合的方式后面会详细说。2. 标准粒子群算法为什么不够用2.1 从单目标PSO到多目标PSO的进化路线粒子群算法的灵感来自鸟群觅食每个粒子代表解空间里的一个候选解。每个粒子有位置向量和速度向量通过个体历史最优pbest和群体历史最优gbest来更新自己的飞行方向。单目标PSO很简单因为只有一个适应度函数pbest和gbest的更新规则一目了然。但是到了多目标问题复杂程度一下子提高了很多——不是因为更新公式变了而是因为“最优”的定义变了。在多目标优化里通常不存在一个解能让所有目标同时达到最优而是存在一组互不支配的解这组解被称为Pareto最优解集。这就带来三个问题第一如何判断两个解之间的优劣关系这需要引入Pareto支配关系——如果解A在所有目标上都不劣于解B且至少在一个目标上优于解B则称A支配B。如果在当前种群中没有任何解可以支配A那么A就是一个非支配解。第二如何从多个Pareto解中选出“群体最优”gbest单目标PSO直接选适应度最好的那个就行多目标PSO里每个粒子的gbest应该从Pareto前沿中选取但选哪个这会影响整个种群的收敛方向。第三如何保证解的多样性如果所有粒子都向同一个Pareto解收敛那么最终得到的Pareto前沿就会非常狭窄虽然每个解都不差但覆盖度不够用户没有选择空间。这就是为什么不能直接把单目标PSO拿过来改两行就用必须对算法框架进行系统性改造。2.2 多目标PSO的三代主流做法对比多目标PSO的发展大致经历了三个阶段理解了这条路线你就能明白小生境改进的价值在哪里。第一代做法是最朴素的“权重和法”把多个目标线性加权变成单目标。这个方法实现最简单但致命的缺陷是权重系数很难定而且线性加权只能搜索到凸Pareto前沿上的点如果前沿是非凸的很多真实Pareto解根本搜不到。第二代做法是引入外部档案External Archive保存非支配解配合网格法或拥挤距离来维护档案规模。这类方法以NSGA-II的思想为基础把Pareto排序和拥挤距离引入粒子群效果比权重和法好很多但存在一个普遍问题——粒子早期的收敛压力不够容易在局部区域扎堆。这是因为在粒子群算法里如果一个粒子周围的非支配解很多它缺乏足够的社会性压力去探索新的区域。第三代做法就是各种基于小生境、分解MOEA/D思想或者指标IBEA思想的改进。小生境技术原本是遗传算法里用来保持种群多样性的经典手段其核心思想是模拟自然界中物种占据不同生态位、各自进化的现象。把它引入粒子群后可以通过粒子间的距离判断它们是否属于同一个“小生境”对密度过高的区域进行抑制促使粒子分散到整个Pareto前沿。这就是我这篇博文标题里“改进”二字的由来。2.3 小生境技术改善多目标搜索的三个关键机制小生境粒子群算法和标准多目标PSO相比核心区别在于以下三个机制第一个机制是种群空间的划分。标准PSO只有一个种群所有粒子共享一个gbest这样很容易被某个局部最优区域“吸”走。小生境PSO会根据粒子在目标空间或决策空间的距离把种群分成若干个小的子种群每个小生境里的粒子主要在局部范围内共享信息跨小生境的交流受控。我实测下来最好的做法是在目标空间和决策空间同时计算粒子间距离然后取加权归一化距离作为小生境划分的依据。纯粹在决策空间划分容易把距离近但目标差异大的粒子分到一起纯粹在目标空间划分又会导致决策空间覆盖度不够。第二个机制是排挤机制Crowding。当某个小生境里的粒子数量超过预设上限时对适应度较差的粒子进行淘汰同时用随机新粒子补充——这里有个关键细节新粒子的生成范围不能是全空间均匀随机最好在当前最优解邻域附近产生否则会破坏已经收敛好的区域。第三个机制是Pareto前沿的精英保留策略。在每个小生境内部维护一个局部精英集在全局维护一个全局精英集。局部精英集用于更新该小生境粒子的pbest全局精英集用于存储整个算法的非支配解。两套精英集同时更新各司其职保证既不过度收敛也不过度发散。3. 小生境粒子群算法的Matlab实现3.1 代码总体架构与文件规划拿到一个优化课题先别急着写算法先用一张思维导图把代码架构理清楚。我这套代码的完整文件构成如下index_main.m % 主程序入口 power_flow.m % 潮流计算函数牛拉法 init_particles.m % 粒子群初始化 niche_partition.m % 小生境划分 update_pbest_gbest.m % 更新局部/全局最优 mutation_operator.m % 变异算子防止陷入局部最优 pareto_dominance.m % 帕累托支配判断 external_archive_update.m % 外部档案更新与裁剪 objective_func.m % 目标函数计算 constraint_handle.m % 约束处理 plot_pareto_front.m % 绘制帕累托前沿这样设计的目的是让每一段逻辑都能单独测试。我第一次做的时候把所有功能堆在一个文件里结果排查bug时痛苦得要命——你想测试小生境划分是否正确但必须把整个优化流程跑一遍才能看到结果。拆成独立函数之后每个模块都可以单独喂入测试数据验证。3.2 粒子编码与初始化细节粒子的编码方式直接决定了搜索空间的维度。对于33节点系统接入4个DG和2组电容器组的情况每个粒子的位置向量进行如下编码% 粒子位置编码说明[PG1,PG2,PG3,PG4, QG1,QG2,QG3,QG4, QC1,QC2] % 其中PG为DG有功出力标幺值QG为DG无功出力标幺值QC为电容器组投切容量标幺值 % 决策变量个数4 4 2 10 nVar 10;初始化时不能太随意否则前期收敛会特别慢。我的做法分两步第一步做一次未优化的潮流计算得到各节点电压和各支路功率以此为基准确定决策变量的可行范围。第二步在可行范围内用拉丁超立方抽样生成初始种群。相比均匀随机抽样拉丁超立方能保证样本在每一维的投影都是均匀的这对多维决策变量特别重要——随机抽样容易产生“空穴”导致某些区域从初始阶段就没人探索。function particles init_particles(pop_size, VarMin, VarMax) nVar length(VarMin); particles zeros(pop_size, nVar); for j 1:nVar % 拉丁超立方抽样 rand_seq (randperm(pop_size) - rand(1, pop_size)) / pop_size; particles(:, j) VarMin(j) rand_seq * (VarMax(j) - VarMin(j)); end end初始化还有一个容易忽略的细节对于电容器组这类离散变量初始化时要直接取整否则后面迭代过程中如果离散变量出现小数最终结果没法工程实现。3.3 小生境划分的核心代码实现这是整个算法里最关键的步骤。我的实现方案是“目标空间决策空间双距离加权”的方法。function [niche_idx, niche_count] niche_partition(pop_obj, pop_dec, niche_radius) % pop_obj: 种群目标函数值矩阵 (N x nObj) % pop_dec: 种群决策变量矩阵 (N x nDec) % niche_radius: 小生境半径 N size(pop_obj, 1); % 目标空间距离归一化 obj_dist pdist2(pop_obj, pop_obj); max_obj_dist max(obj_dist(:)); obj_dist obj_dist / max_obj_dist; % 决策空间距离归一化 dec_dist pdist2(pop_dec, pop_dec); max_dec_dist max(dec_dist(:)); dec_dist dec_dist / max_dec_dist; % 加权距离目标空间权重0.6决策空间权重0.4 combined_dist 0.6 * obj_dist 0.4 * dec_dist; % 基于距离的划分简单聚类方法 niche_idx zeros(N, 1); niche_count zeros(N, 1); current_niche 0; for i 1:N if niche_idx(i) 0 current_niche current_niche 1; neighbors find(combined_dist(i, :) niche_radius); % 仅将未分配且距离小于半径的粒子划入同一小生境 for j neighbors if niche_idx(j) 0 niche_idx(j) current_niche; niche_count(current_niche) niche_count(current_niche) 1; end end end end end这段代码的关键在于半径的选择。半径太大所有粒子都被分到同一个小生境算法退化为标准PSO半径太小每个粒子独霸一个小生境粒子间缺乏信息交互收敛速度极慢。我测试下来对于33节点系统当决策变量维度是10、种群规模是100时小生境半径取0.35左右效果最佳。不同系统需要配合网格搜索调整后面会详细讲。3.4 多目标处理与外部档案更新机制在多目标粒子群算法里外部档案是保存Pareto最优解的地方它的更新机制决定了最终优化结果的质量。function archive external_archive_update(archive, candidates, max_size) % 候选解合并且去除支配解 combined [archive; candidates]; n size(combined, 1); % 逐对判断支配关系 dominated_flag false(n, 1); for i 1:n for j 1:n if i ~ j dominates(combined(j, :), combined(i, :)) dominated_flag(i) true; break; end end end non_dominated combined(~dominated_flag, :); % 如果非支配解数量超过上限用拥挤距离裁剪 if size(non_dominated, 1) max_size crowding_dist crowding_distance(non_dominated); [~, idx_sorted] sort(crowding_dist, descend); archive non_dominated(idx_sorted(1:max_size), :); else archive non_dominated; end end拥挤距离的计算参考了NSGA-II的做法对每个目标函数值排序后计算相邻解之间的归一化距离。拥挤距离大的解代表它所在的区域比较稀疏应该优先保留拥挤距离小的解说明它周围密密麻麻都是解删掉它对多样性影响不大。这里有一个容易踩的坑更新gbest时不能简单地从外部档案里随机选一个解作为全局最优。我建议的做法是先判断当前粒子所在的小生境然后从该小生境对应的局部精英集中选择一个拥挤距离最大的解或者以一定概率选择档案中的稀疏区域解作为gbest。这样才能保证种群整体持续向Pareto前沿逼近。3.5 约束处理的罚函数与可行性修正组合方案配电网优化里的约束条件很多直接全用罚函数会导致罚函数系数难调——系数太小约束满足不了系数太大算法不收敛。我用的是“罚函数 可行性优先”的组合方案对于电压越限和潮流不收敛用二次罚函数罚系数取10000足够大到保证最终解不会越限。对于DG视在功率约束不用罚函数用可行性修正。因为Q_DG的上限是随P_DG变化的修正起来比较直观——如果某粒子的P_DG和Q_DG超过了视在功率限制直接把Q_DG缩放到满足约束的范围。function x constraint_handle(x) % x [PG1,PG2,PG3,PG4, QG1,QG2,QG3,QG4, QC1,QC2] S_DG_max [0.5, 0.4, 0.3, 0.4]; % DG视在功率上限标幺值 % 视在功率约束修正 for i 1:4 P x(i); Q x(4 i); S_max S_DG_max(i); if sqrt(P^2 Q^2) S_max ratio S_max / sqrt(P^2 Q^2); P P * ratio; Q Q * ratio; x(i) P; x(4 i) Q; end end % 电容器投切档位取整 x(9) round(x(9)); x(10) round(x(10)); end这个约束处理函数在粒子更新后和计算适应度前都要调用一次确保每个参与评价的解都有物理意义。4. 米拉潮流计算与目标函数实现4.1 潮流计算在优化算法里的调用方式粒子群优化算法的每一次适应度评估都需要计算一次配电网潮流。100个粒子迭代200次就是20000次潮流计算。如果潮流计算函数写得不够高效整个算法的运行时间会从几分钟膨胀到几个小时。我在这个项目里用的是牛顿-拉夫逊法对33节点这种小规模系统来说潮流计算的单次耗时只要几毫秒所以不存在性能瓶颈。但要注意几个实现细节第一潮流计算的收敛精度不需要太高。优化过程中粒子的很多候选解都是不可行的电压严重越限甚至潮流不收敛对这些解你不需要精确的潮流结果只要知道它不可行就够了。我把潮流收敛阈值设为1e-6但在粒子群搜索前期如果一个候选解在5次迭代内还不收敛直接判定为不可行解返回一个很大的惩罚值。这样能节省大量计算时间。第二辐射状配电网用前推回代法效率更高但牛顿法通用性更强。如果你后续要扩展算例到弱环网或输配协同建议直接用牛顿法避免整套框架返工。第三潮流计算输入输出要设计成标准接口这样后续想换成DistFlow模型或者三相不平衡潮流都很方便。function [V, Ploss, success] power_flow(bus, branch, P_DG, Q_DG, Q_C) % bus: 节点参数矩阵 % branch: 支路参数矩阵 % P_DG: DG有功出力向量 % Q_DG: DG无功出力向量 % Q_C: 电容器无功补偿向量 % 返回: 节点电压 V, 网损 Ploss, 是否收敛 success % 构建节点注入功率 n_bus size(bus, 1); P_inj zeros(n_bus, 1); Q_inj zeros(n_bus, 1); % ... 将DG和电容器注入功率叠加到对应节点 % 牛顿-拉夫逊潮流 % ... (标准实现关键在于雅可比矩阵的构建和修正方程的求解) % 计算网损 Ploss sum(real(S_branch_loss)); end4.2 目标函数计算的细节处理目标函数是整个优化过程的指挥棒它的设计直接决定了算法会朝着什么方向收敛。我的目标函数实现如下function [obj1, obj2] objective_func(V, branch, P_load_total, P_DG_total) % 目标1系统网损最小 obj1 Ploss; % 目标2电压偏移最小所有负荷节点的电压偏差绝对值之和 V_ref 1.0; % 标幺值基准电压 V_deviation sum(abs(V - V_ref)); obj2 V_deviation / length(V); % 取平均消除节点数影响 % 罚函数处理电压越限惩罚 V_min 0.95; V_max 1.05; penalty 0; for i 1:length(V) if V(i) V_min penalty penalty (V_min - V(i))^2; elseif V(i) V_max penalty penalty (V(i) - V_max)^2; end end obj1 obj1 10000 * penalty; obj2 obj2 10000 * penalty; end电压偏移的目标函数我特意用了“平均值”而不是“最大值”是有讲究的。如果目标是最大电压偏差最小算法会优先照顾电压最差的节点这看起来公平但实际上会让算法过度关注极端节点而在整体范围内搜索不足。用平均值的话算法会综合考虑所有节点的电压质量同时罚函数保证了极端节点也不会被忽视——这是一种平衡策略。4.3 目标函数的标准化处理多目标优化里还有一个关键细节——目标函数之间的量纲差异。网损的单位是kW或者标幺值功率数值通常在几十到几百的范围内但电压偏移是标幺值数值在0.001到0.05之间。如果直接拿这两个数去做支配判断电压偏移就会被网损“淹没”——网损稍微变一下电压偏移的变化根本看不出来支配关系完全由网损决定多目标优化退化成了单目标优化。解决办法是对目标函数进行标准化。最常用的是min-max归一化。在外部档案更新过程中记录每个目标的历史最小值和最大值然后对每个解的目标值做归一化处理。% 目标函数归一化处理 function [obj1_norm, obj2_norm] normalize_objectives(obj1, obj2, min_max_records) obj1_min min_max_records.obj1_min; obj1_max min_max_records.obj1_max; obj2_min min_max_records.obj2_min; obj2_max min_max_records.obj2_max; obj1_norm (obj1 - obj1_min) / (obj1_max - obj1_min 1e-6); obj2_norm (obj2 - obj2_min) / (obj2_max - obj2_min 1e-6); end注意加了一个1e-6的常数防止分母为零。这个标准化操作直接影响Pareto前沿分布的均匀性也是很多复现论文的代码和原文结果对不上的重要原因之一。5. 主程序循环与算法流程整合5.1 主算法迭代框架把前面所有模块组装起来形成完整的主程序循环。我写的是经典的小生境多目标PSO主循环参考了多篇文献的实现思路并根据实际测试做了调整% 参数设置 pop_size 100; % 种群规模 max_iter 200; % 最大迭代次数 c1 1.5; c2 1.5; % 学习因子 w_max 0.9; w_min 0.4; % 惯性权重范围 niche_radius 0.35; % 小生境半径 archive_size 50; % 外部档案最大容量 % 初始化 particles init_particles(pop_size, VarMin, VarMax); for i 1:pop_size [V, Ploss, ~] power_flow(bus, branch, particles(i, 1:4), particles(i, 5:8), particles(i, 9:10)); [Obj1(i), Obj2(i)] objective_func(V, Ploss); pbest(i, :) particles(i, :); pbest_obj1(i) Obj1(i); pbest_obj2(i) Obj2(i); end [pgbest, gbest_obj1, gbest_obj2] get_global_best(particles, Obj1, Obj2, archive); % 迭代主循环 for iter 1:max_iter w w_max - (w_max - w_min) * iter / max_iter; % 小生境划分每5代重新划分一次降低计算量 if mod(iter, 5) 1 [niche_idx, ~] niche_partition([Obj1, Obj2], particles, niche_radius); end % 更新粒子速度和位置 for i 1:pop_size v(i, :) w * v(i, :) ... c1 * rand(1, nVar) .* (pbest(i, :) - particles(i, :)) ... c2 * rand(1, nVar) .* (pgbest(i, :) - particles(i, :)); % 速度限幅 v(i, :) max(min(v(i, :), Vmax), Vmin); particles(i, :) particles(i, :) v(i, :); % 位置边界处理 particles(i, :) max(min(particles(i, :), VarMax), VarMin); % 约束修正 particles(i, :) constraint_handle(particles(i, :)); end % 变异操作防止早熟收敛 mutation_rate 0.1 - 0.07 * iter / max_iter; for i 1:pop_size if rand mutation_rate k randi(nVar); particles(i, k) VarMin(k) rand * (VarMax(k) - VarMin(k)); particles(i, :) constraint_handle(particles(i, :)); end end % 计算适应度 for i 1:pop_size [V, Ploss, success] power_flow(bus, branch, particles(i, 1:4), particles(i, 5:8), particles(i, 9:10)); if ~success Obj1(i) 1e10; Obj2(i) 1e10; continue; end [Obj1(i), Obj2(i)] objective_func(V, Ploss); end % 更新pbest和外部档案 for i 1:pop_size if dominates([Obj1(i), Obj2(i)], [pbest_obj1(i), pbest_obj2(i)]) pbest(i, :) particles(i, :); pbest_obj1(i) Obj1(i); pbest_obj2(i) Obj2(i); end end % 更新外部档案含小生境分布控制 archive external_archive_update(archive, [particles, Obj1, Obj2], archive_size); % 从小生境局部精英集中选择gbest pgbest select_gbest_from_niche(niche_idx, archive, pop_size); % 输出迭代进度 if mod(iter, 20) 0 fprintf(迭代次数: %d/%d, 外部档案解个数: %d\n, iter, max_iter, size(archive, 1)); end end5.2 变异操作的关键——自适应变异率很多多目标PSO的实现里变异率是固定的但我实测下来固定变异率效果很差。前期如果变异率太高算法发散好不容易找到的优质区域被变异粒子给破坏了后期如果变异率太低算法早熟所有粒子挤在一个局部Pareto前沿上出不来。我用的方案是变异率从0.1线性衰减到0.03。前期允许较大的变异概率来增加种群多样性后期减少变异让种群精细搜索当前区域。变异操作本身只随机改变一个维度这样既不会破坏太多已经找到的好结构又能提供潜在的跳跃能力。5.3 存档更新与拥挤距离计算这里单独说一下拥挤距离的实现因为这个细节直接影响最终Pareto前沿的质量。拥挤距离的计算思路是对一个粒子的某个目标值找到比它大的最近解和比它小的最近解计算这两个解在该目标上的差值然后除以该目标的最大值范围做归一化。把所有目标的归一化差值加起来就是这个粒子的拥挤距离。function cd crowding_distance(obj_matrix) % obj_matrix: N x nObj 的目标函数值矩阵 N size(obj_matrix, 1); nObj size(obj_matrix, 2); cd zeros(N, 1); for j 1:nObj [sorted_val, idx] sort(obj_matrix(:, j)); cd(idx(1)) inf; cd(idx(N)) inf; obj_range sorted_val(N) - sorted_val(1); if obj_range 1e-10 continue; end for i 2:N-1 cd(idx(i)) cd(idx(i)) (sorted_val(i1) - sorted_val(i-1)) / obj_range; end end end边界解的拥挤距离设为无穷大这是有意为之——Pareto前沿的两端是探索新区域的桥头堡必须优先保留。6. 仿真结果与Pareto前沿分析6.1 仿真环境与参数配置我的仿真环境是Matlab R2021bIntel i7处理器16GB内存运行时间大约2分钟到3分钟。系统采用IEEE 33节点配电网基准电压12.66kV总负荷有功3715kW无功2300kvar。DG配置如下节点18接入一个额定容量500kW的分布式电源节点22接入400kW节点25接入300kW节点33接入400kW。所有DG都具备无功调节能力逆变器剩余容量可以发无功无功调节范围设定为其视在功率的-40%到40%。此外在节点18和33各配置了一组电容器单组容量50kvar可以多档位投切总共最多10组并联。这样的配置和很多文献里的算例是一致的方便你做结果对照。6.2 小生境PSO与标准PSO的结果对比为了验证小生境PSO的效果我在同一仿真条件下跑了标准多目标PSO作为对照组两种算法共用相同的种群规模、迭代次数和目标函数。先说网损结果标准PSO找到的最优网损是151.2kW小生境PSO找到的最优网损是145.8kW降低了约3.6%。这个数值看起来不大但放在配电网优化里已经是不错的提升了而且关键在于两种算法的“搜索路线”完全不同。更显著的差异在Pareto前沿覆盖度上。标准PSO的Pareto前沿集中在网损目标较好的区域电压偏移目标的变化范围只有0.01左右而小生境PSO的Pareto前沿从低网损到低电压偏移完整铺开电压偏移的变化范围达到0.04以上。这意味着小生境PSO给决策者提供了更多样的选择——如果你想优先保证电压质量可以选择牺牲一部分网损的解如果更看重经济性就选网损最小的解。6.3 Pareto前沿的分布质量分析我拿Pareto前沿的分布质量做了一个量化评估用的指标包括C指标覆盖率衡量集合A的解有多少被集合B的解支配。标准PSO的非支配解集合被小生境PSO的解支配的比例达到83%说明小生境PSO寻找的Pareto前沿在整体上占优。间距指标Spacing衡量Pareto前沿上解之间的均匀程度。间距越小说明解在前沿上分布越均匀决策者对中间区域的把控越精准。小生境PSO的间距为0.023标准PSO为0.047前者几乎是后者的一半分布均匀性大幅提升。最大覆盖范围Maximum Spread衡量Pareto前沿在两个目标上的延展范围。小生境PSO的最大覆盖范围为0.87标准PSO为0.62。这直接说明小生境机制有效防止了种群过度收敛保持了更好的多样性。6.4 从Pareto解集中选取折中解工程实际里最终只需要一个调度方案不是一整条前沿。那怎么从Pareto解集里选出最终的折中解最常用的方法是模糊隶属度法。对于每个Pareto解在每个目标上计算模糊隶属度把所有目标的隶属度求和选择总和最大的解作为折中解。简单来说就是找那个“在两方面都不算太差的解”。我选取的折中解对应的控制策略是DG1出力498kW、DG2出力386kW、DG3出力295kW、DG4出力399kW总DG出力1578kW。电容器组1投入8组400kvar电容器组2投入7组350kvar。在这个方案下系统网损150.3kW平均电压偏移0.018所有节点电压都在0.97p.u.到1.03p.u.之间全部满足要求。对比不进行优化时的初始状态——网损202.5kW平均电压偏移0.035部分末端节点电压低于0.95p.u.——优化效果非常明显。7. 常见问题与调试心得7.1 粒子群不收敛或前期发散怎么办这个问题是最常见的。如果你的算法前期发散大概率不是算法框架的问题而是参数设置的问题。先检查惯性权重w和学习因子c1、c2。w过大粒子飞行太快容易越过最优区域c1和c2设置不当会导致粒子过度跟随个体经验或群体经验缺乏自主探索能力。我的建议是w从上界0.9线性衰减到下界0.4c1和c2都设为1.5这是最常用且经过大量验证的参数组合。再检查速度限幅Vmax。速度限幅的设置直接影响粒子在解空间中的移动步长如果Vmax太大粒子会剧烈震荡无法精细搜索如果太小粒子会爬行式前进收敛极慢。我建议Vmax设置为各维度变量范围的20%左右——对于有功出力在0到1之间的变量Vmax设为0.2对于电容器组数在0到10之间的整数变量Vmax设为2。还有一个容易被忽略的原因是潮流计算的收敛性。如果潮流计算函数存在bug导致很多粒子被误判为不可行解算法会在一个扭曲的适应度地貌上搜索自然无法收敛。排查方法是单独跑一次潮流把计算结果和MATPOWER或者公开文献里的结果对照确认无误后再接入优化算法。7.2 小生境划分过度导致收敛慢小生境半径选的太小粒子间信息交互大幅减少收敛会非常慢。我在调试时曾经遇到过运行500代都找不到满意解的情况最后发现是小生境半径只有0.1导致的。有一个实用的判断标准如果种群规模是100期望划分的小生境数量大概在8到12个之间。如果超过15个说明半径偏小如果少于5个说明半径偏大。先用两次预实验确定合适的半径范围再精细调整。另外要注意小生境的划分频率不需要每一代都做。粒子位置的更新是渐进的连续两代之间的小生境划分结果不会有太大变化每5代重新划分一次既能保证小生境信息不过期又能显著降低计算开销。7.3 外部档案中Pareto解数量不够怎么办外部档案里最终只有两三个解这种情况通常说明种群多样性严重不足。排除小生境机制的干扰后最常见的原因是变异机制失效。如果变异率太早衰减到零粒子很快就失去探索新区域的动力全部被吸附到当前最优解附近。我的做法是让变异率在最坏情况下也不低于0.02确保后期仍有小概率产生跳跃。另一个原因是外部档案的容量设置是否合理。如果容量太小比如只有10个即使算法发现了充足的非支配解档案更新机制也会把它们裁剪掉导致最终保留下来的解太少。建议至少设置为种群规模的50%。7.4 代码调试的三个黄金技巧写这种复杂的优化代码bug排查是很痛苦的过程我总结了一个三段式调试法第一阶段模块级验证。在把潮流计算、小生境划分、外部档案更新等模块接入主循环前先用构造好的小规模数据单独测试每个模块的输出是否符合预期。比如小生境划分模块我构造了10个已知坐标的粒子手动算好期望的划分结果再和程序输出比对。第二阶段单目标验证。把多目标问题退化为单目标问题只优化网损用标准PSO跑一遍把优化结果和文献值对照。如果单目标都跑不对说明算法实现有基础性错误多目标肯定更不对。第三阶段逐步恢复多目标。在单目标跑通的基础上逐步加入第二个目标、小生境机制、外部档案更新每加一个模块就跑一遍对比加模块前后的结果差异。这样能精确定位是哪个环节引入的问题。8. 后续扩展方向与几点经验总结代码写到这里主体功能已经完整。如果你的课题需要进一步扩展这里有几个我非常推荐的方向。第一个扩展方向是接入时序特性。当前的优化是一个静态截面也就是在某一时刻的DG出力和负荷水平下求解最优。但实际上DG出力和负荷是随时间变化的更贴近工程实际的做法是做96个断面每15分钟一个的时序优化或者使用典型日场景进行多场景随机优化。扩展时主要改造目标函数和约束条件算法框架不用动。第二个扩展方向是多种DG类型的建模。光伏、风电、储能、微型燃气轮机的运行特性差异很大光伏的无功能力由逆变器决定风电要考虑到风速的随机性储能还要考虑充放电状态和SOC约束。目前的决策变量编码方式有很强的通用性加变量类型只需扩展编码维度就行。第三个扩展方向是分布式求解。如果配电网规模从33节点扩展到几百个节点的实际馈线集中式PSO的计算负担会显著增加。可以考虑用分布式优化框架把网络分区每个分区跑一个独立的小生境PSO分区之间只交换边界信息。这和很多主动配电网分层分布式的理念完全契合。最后说几句掏心窝的总结。做智能算法和配电网优化结合的研究最大的误区是把精力全花在算法创新上把配电网模型做成“搭积木”——只是为了跑通算法而随便应付。我过去吃过亏前前后后调了大半年算法结果换了另外一个配电网模型来测试之前的很多经验全废了。后来我才意识到配电网物理约束的准确性才是决定优化结果能否落地的那条腿。真正好的研究工作一定是两条腿走路——算法有创新模型有深度两者互为表里才能产出既好看又实用的成果。这套小生境粒子群版本的协调优化代码在IEEE 33节点系统上的仿真结果表明它确实能在网损和电压质量之间找到一组分布均匀、覆盖面广的Pareto最优解。希望这篇文章里的思路和踩过的坑能让你在自己的课题里少走几步弯路。 改进小生境多目标粒子群的配电网有功无功协调优化Matlab实现
先说明一点做配电网优化方向的同学不管你是写毕业论文、发小论文还是做横向课题肯定绕不开“无功优化”和“有功调度”这两个老问题。但如果你把两者分开做往往得到一个看似收敛、实际顾此失彼的结果。这篇博文记录的就是我把两者真正“揉”在一起的完整过程——基于改进多目标粒子群算法小生境粒子群的配电网有功-无功协调优化配全套Matlab代码实现方法。我会从问题建模、算法改进思路、代码架构、IEEE 33节点系统实测结果到调参和踩坑一次性讲透方便你直接“抄作业”。1. 项目概述与问题定义1.1 为什么要把有功和无功放在一起优化传统配电网运行调度里有功调度和无功调节通常是分开做的。有功调度主要盯发电机的出力分配目标是经济性无功调节主要管无功补偿设备电容器组、SVC、SVG、OLTC等目标是电压合格率和降低网损。表面上看这没什么问题但在分布式电源DG大量接入之后情况完全变了。DG的出力是波动的甚至有反调峰特性它同时改变配电网的有功潮流分布和无功电压分布。如果你只做无功优化假定有功出力是固定的那在DG出力大幅变化时电压可能越限无功补偿设备动作频繁。如果你只做有功调度忽略无功和电压约束那优化出来的DG出力方案可能根本没法落地因为节点电压可能已经超出允许范围。换句话说有功和无功在配电网里是强耦合的。把它们分开优化本质上是在做一个“次优解”——这也是大多数基于PSO的配电网优化论文只做单一目标时仿真结果看起来挺好、但实际工程没法用的根本原因。把两者统一到一个多目标优化框架里让算法同时搜索有功出力和无功补偿的可行解空间这才能得到真正可执行的调度策略。1.2 有源配电网协调优化的两个核心目标我做这个项目时目标函数选了两个不贪多因为目标越多Pareto前沿的搜索难度就越大而且工程上最关注的其实就这两个第一个目标是系统网损最小。配电网的网损包括有功损耗和无功损耗但在目标函数里我们只关心有功损耗。无功功率在传输过程中不直接产生电度损耗但它会占据线路容量、增大电流从而间接增大有功损耗。所以降低网损的实质是在满足电压约束的前提下尽量减少无功功率在配电网中的长距离流动同时优化DG有功出力减少倒送功率。第二个目标是电压偏移最小。电压质量是配电网运行的生命线。分布式电源接入后最常见的电压问题就是末端电压抬升和局部节点电压越限。我用所有负荷节点电压与额定电压的偏差绝对值之和或者最大偏差来衡量这个目标数值越小代表电压质量越好。当然你也可以根据需要加上第三个目标——比如DG运行成本最低、开关动作次数最少、碳排放最小等代码架构上完全可以扩展后面会讲。1.3 配电网模型与决策变量的选取先说测试系统我用的是最经典的IEEE 33节点配电网系统。为什么选它因为它的参数公开、文献对照多、算例结果可查而且网络结构是一个典型的辐射状配电网包含联络开关负荷分布不均匀电压问题突出适合用来验证优化算法的有效性。关于决策变量我分成了两类第一类是有功调度变量就是各分布式电源的有功出力。假设系统里接了若干个DG比如在节点18、22、25、33接入光伏或风机那么每个DG的有功出力就是一个决策变量。需要注意的是DG有功出力不是无限可调的要设置上下限通常取其额定容量的0%到100%或者根据预测出力设置一个可调范围。第二类是无功控制变量包括无功补偿装置的投切容量和DG的无功出力。如果用离散的电容器组那么决策变量是整数代表投切的组数如果是SVG或DG逆变器可以连续调节。DG逆变器参与无功调节时要满足视在功率约束——有功和无功不能同时达到最大值即[ Q_{DG} \le \sqrt{S_{DG}^2 - P_{DG}^2} ]这个约束在粒子群算法里处理起来比较麻烦因为它是一个动态的、与另一个决策变量耦合的约束。我采用的是罚函数法加可行性修正相结合的方式后面会详细说。2. 标准粒子群算法为什么不够用2.1 从单目标PSO到多目标PSO的进化路线粒子群算法的灵感来自鸟群觅食每个粒子代表解空间里的一个候选解。每个粒子有位置向量和速度向量通过个体历史最优pbest和群体历史最优gbest来更新自己的飞行方向。单目标PSO很简单因为只有一个适应度函数pbest和gbest的更新规则一目了然。但是到了多目标问题复杂程度一下子提高了很多——不是因为更新公式变了而是因为“最优”的定义变了。在多目标优化里通常不存在一个解能让所有目标同时达到最优而是存在一组互不支配的解这组解被称为Pareto最优解集。这就带来三个问题第一如何判断两个解之间的优劣关系这需要引入Pareto支配关系——如果解A在所有目标上都不劣于解B且至少在一个目标上优于解B则称A支配B。如果在当前种群中没有任何解可以支配A那么A就是一个非支配解。第二如何从多个Pareto解中选出“群体最优”gbest单目标PSO直接选适应度最好的那个就行多目标PSO里每个粒子的gbest应该从Pareto前沿中选取但选哪个这会影响整个种群的收敛方向。第三如何保证解的多样性如果所有粒子都向同一个Pareto解收敛那么最终得到的Pareto前沿就会非常狭窄虽然每个解都不差但覆盖度不够用户没有选择空间。这就是为什么不能直接把单目标PSO拿过来改两行就用必须对算法框架进行系统性改造。2.2 多目标PSO的三代主流做法对比多目标PSO的发展大致经历了三个阶段理解了这条路线你就能明白小生境改进的价值在哪里。第一代做法是最朴素的“权重和法”把多个目标线性加权变成单目标。这个方法实现最简单但致命的缺陷是权重系数很难定而且线性加权只能搜索到凸Pareto前沿上的点如果前沿是非凸的很多真实Pareto解根本搜不到。第二代做法是引入外部档案External Archive保存非支配解配合网格法或拥挤距离来维护档案规模。这类方法以NSGA-II的思想为基础把Pareto排序和拥挤距离引入粒子群效果比权重和法好很多但存在一个普遍问题——粒子早期的收敛压力不够容易在局部区域扎堆。这是因为在粒子群算法里如果一个粒子周围的非支配解很多它缺乏足够的社会性压力去探索新的区域。第三代做法就是各种基于小生境、分解MOEA/D思想或者指标IBEA思想的改进。小生境技术原本是遗传算法里用来保持种群多样性的经典手段其核心思想是模拟自然界中物种占据不同生态位、各自进化的现象。把它引入粒子群后可以通过粒子间的距离判断它们是否属于同一个“小生境”对密度过高的区域进行抑制促使粒子分散到整个Pareto前沿。这就是我这篇博文标题里“改进”二字的由来。2.3 小生境技术改善多目标搜索的三个关键机制小生境粒子群算法和标准多目标PSO相比核心区别在于以下三个机制第一个机制是种群空间的划分。标准PSO只有一个种群所有粒子共享一个gbest这样很容易被某个局部最优区域“吸”走。小生境PSO会根据粒子在目标空间或决策空间的距离把种群分成若干个小的子种群每个小生境里的粒子主要在局部范围内共享信息跨小生境的交流受控。我实测下来最好的做法是在目标空间和决策空间同时计算粒子间距离然后取加权归一化距离作为小生境划分的依据。纯粹在决策空间划分容易把距离近但目标差异大的粒子分到一起纯粹在目标空间划分又会导致决策空间覆盖度不够。第二个机制是排挤机制Crowding。当某个小生境里的粒子数量超过预设上限时对适应度较差的粒子进行淘汰同时用随机新粒子补充——这里有个关键细节新粒子的生成范围不能是全空间均匀随机最好在当前最优解邻域附近产生否则会破坏已经收敛好的区域。第三个机制是Pareto前沿的精英保留策略。在每个小生境内部维护一个局部精英集在全局维护一个全局精英集。局部精英集用于更新该小生境粒子的pbest全局精英集用于存储整个算法的非支配解。两套精英集同时更新各司其职保证既不过度收敛也不过度发散。3. 小生境粒子群算法的Matlab实现3.1 代码总体架构与文件规划拿到一个优化课题先别急着写算法先用一张思维导图把代码架构理清楚。我这套代码的完整文件构成如下index_main.m % 主程序入口 power_flow.m % 潮流计算函数牛拉法 init_particles.m % 粒子群初始化 niche_partition.m % 小生境划分 update_pbest_gbest.m % 更新局部/全局最优 mutation_operator.m % 变异算子防止陷入局部最优 pareto_dominance.m % 帕累托支配判断 external_archive_update.m % 外部档案更新与裁剪 objective_func.m % 目标函数计算 constraint_handle.m % 约束处理 plot_pareto_front.m % 绘制帕累托前沿这样设计的目的是让每一段逻辑都能单独测试。我第一次做的时候把所有功能堆在一个文件里结果排查bug时痛苦得要命——你想测试小生境划分是否正确但必须把整个优化流程跑一遍才能看到结果。拆成独立函数之后每个模块都可以单独喂入测试数据验证。3.2 粒子编码与初始化细节粒子的编码方式直接决定了搜索空间的维度。对于33节点系统接入4个DG和2组电容器组的情况每个粒子的位置向量进行如下编码% 粒子位置编码说明[PG1,PG2,PG3,PG4, QG1,QG2,QG3,QG4, QC1,QC2] % 其中PG为DG有功出力标幺值QG为DG无功出力标幺值QC为电容器组投切容量标幺值 % 决策变量个数4 4 2 10 nVar 10;初始化时不能太随意否则前期收敛会特别慢。我的做法分两步第一步做一次未优化的潮流计算得到各节点电压和各支路功率以此为基准确定决策变量的可行范围。第二步在可行范围内用拉丁超立方抽样生成初始种群。相比均匀随机抽样拉丁超立方能保证样本在每一维的投影都是均匀的这对多维决策变量特别重要——随机抽样容易产生“空穴”导致某些区域从初始阶段就没人探索。function particles init_particles(pop_size, VarMin, VarMax) nVar length(VarMin); particles zeros(pop_size, nVar); for j 1:nVar % 拉丁超立方抽样 rand_seq (randperm(pop_size) - rand(1, pop_size)) / pop_size; particles(:, j) VarMin(j) rand_seq * (VarMax(j) - VarMin(j)); end end初始化还有一个容易忽略的细节对于电容器组这类离散变量初始化时要直接取整否则后面迭代过程中如果离散变量出现小数最终结果没法工程实现。3.3 小生境划分的核心代码实现这是整个算法里最关键的步骤。我的实现方案是“目标空间决策空间双距离加权”的方法。function [niche_idx, niche_count] niche_partition(pop_obj, pop_dec, niche_radius) % pop_obj: 种群目标函数值矩阵 (N x nObj) % pop_dec: 种群决策变量矩阵 (N x nDec) % niche_radius: 小生境半径 N size(pop_obj, 1); % 目标空间距离归一化 obj_dist pdist2(pop_obj, pop_obj); max_obj_dist max(obj_dist(:)); obj_dist obj_dist / max_obj_dist; % 决策空间距离归一化 dec_dist pdist2(pop_dec, pop_dec); max_dec_dist max(dec_dist(:)); dec_dist dec_dist / max_dec_dist; % 加权距离目标空间权重0.6决策空间权重0.4 combined_dist 0.6 * obj_dist 0.4 * dec_dist; % 基于距离的划分简单聚类方法 niche_idx zeros(N, 1); niche_count zeros(N, 1); current_niche 0; for i 1:N if niche_idx(i) 0 current_niche current_niche 1; neighbors find(combined_dist(i, :) niche_radius); % 仅将未分配且距离小于半径的粒子划入同一小生境 for j neighbors if niche_idx(j) 0 niche_idx(j) current_niche; niche_count(current_niche) niche_count(current_niche) 1; end end end end end这段代码的关键在于半径的选择。半径太大所有粒子都被分到同一个小生境算法退化为标准PSO半径太小每个粒子独霸一个小生境粒子间缺乏信息交互收敛速度极慢。我测试下来对于33节点系统当决策变量维度是10、种群规模是100时小生境半径取0.35左右效果最佳。不同系统需要配合网格搜索调整后面会详细讲。3.4 多目标处理与外部档案更新机制在多目标粒子群算法里外部档案是保存Pareto最优解的地方它的更新机制决定了最终优化结果的质量。function archive external_archive_update(archive, candidates, max_size) % 候选解合并且去除支配解 combined [archive; candidates]; n size(combined, 1); % 逐对判断支配关系 dominated_flag false(n, 1); for i 1:n for j 1:n if i ~ j dominates(combined(j, :), combined(i, :)) dominated_flag(i) true; break; end end end non_dominated combined(~dominated_flag, :); % 如果非支配解数量超过上限用拥挤距离裁剪 if size(non_dominated, 1) max_size crowding_dist crowding_distance(non_dominated); [~, idx_sorted] sort(crowding_dist, descend); archive non_dominated(idx_sorted(1:max_size), :); else archive non_dominated; end end拥挤距离的计算参考了NSGA-II的做法对每个目标函数值排序后计算相邻解之间的归一化距离。拥挤距离大的解代表它所在的区域比较稀疏应该优先保留拥挤距离小的解说明它周围密密麻麻都是解删掉它对多样性影响不大。这里有一个容易踩的坑更新gbest时不能简单地从外部档案里随机选一个解作为全局最优。我建议的做法是先判断当前粒子所在的小生境然后从该小生境对应的局部精英集中选择一个拥挤距离最大的解或者以一定概率选择档案中的稀疏区域解作为gbest。这样才能保证种群整体持续向Pareto前沿逼近。3.5 约束处理的罚函数与可行性修正组合方案配电网优化里的约束条件很多直接全用罚函数会导致罚函数系数难调——系数太小约束满足不了系数太大算法不收敛。我用的是“罚函数 可行性优先”的组合方案对于电压越限和潮流不收敛用二次罚函数罚系数取10000足够大到保证最终解不会越限。对于DG视在功率约束不用罚函数用可行性修正。因为Q_DG的上限是随P_DG变化的修正起来比较直观——如果某粒子的P_DG和Q_DG超过了视在功率限制直接把Q_DG缩放到满足约束的范围。function x constraint_handle(x) % x [PG1,PG2,PG3,PG4, QG1,QG2,QG3,QG4, QC1,QC2] S_DG_max [0.5, 0.4, 0.3, 0.4]; % DG视在功率上限标幺值 % 视在功率约束修正 for i 1:4 P x(i); Q x(4 i); S_max S_DG_max(i); if sqrt(P^2 Q^2) S_max ratio S_max / sqrt(P^2 Q^2); P P * ratio; Q Q * ratio; x(i) P; x(4 i) Q; end end % 电容器投切档位取整 x(9) round(x(9)); x(10) round(x(10)); end这个约束处理函数在粒子更新后和计算适应度前都要调用一次确保每个参与评价的解都有物理意义。4. 米拉潮流计算与目标函数实现4.1 潮流计算在优化算法里的调用方式粒子群优化算法的每一次适应度评估都需要计算一次配电网潮流。100个粒子迭代200次就是20000次潮流计算。如果潮流计算函数写得不够高效整个算法的运行时间会从几分钟膨胀到几个小时。我在这个项目里用的是牛顿-拉夫逊法对33节点这种小规模系统来说潮流计算的单次耗时只要几毫秒所以不存在性能瓶颈。但要注意几个实现细节第一潮流计算的收敛精度不需要太高。优化过程中粒子的很多候选解都是不可行的电压严重越限甚至潮流不收敛对这些解你不需要精确的潮流结果只要知道它不可行就够了。我把潮流收敛阈值设为1e-6但在粒子群搜索前期如果一个候选解在5次迭代内还不收敛直接判定为不可行解返回一个很大的惩罚值。这样能节省大量计算时间。第二辐射状配电网用前推回代法效率更高但牛顿法通用性更强。如果你后续要扩展算例到弱环网或输配协同建议直接用牛顿法避免整套框架返工。第三潮流计算输入输出要设计成标准接口这样后续想换成DistFlow模型或者三相不平衡潮流都很方便。function [V, Ploss, success] power_flow(bus, branch, P_DG, Q_DG, Q_C) % bus: 节点参数矩阵 % branch: 支路参数矩阵 % P_DG: DG有功出力向量 % Q_DG: DG无功出力向量 % Q_C: 电容器无功补偿向量 % 返回: 节点电压 V, 网损 Ploss, 是否收敛 success % 构建节点注入功率 n_bus size(bus, 1); P_inj zeros(n_bus, 1); Q_inj zeros(n_bus, 1); % ... 将DG和电容器注入功率叠加到对应节点 % 牛顿-拉夫逊潮流 % ... (标准实现关键在于雅可比矩阵的构建和修正方程的求解) % 计算网损 Ploss sum(real(S_branch_loss)); end4.2 目标函数计算的细节处理目标函数是整个优化过程的指挥棒它的设计直接决定了算法会朝着什么方向收敛。我的目标函数实现如下function [obj1, obj2] objective_func(V, branch, P_load_total, P_DG_total) % 目标1系统网损最小 obj1 Ploss; % 目标2电压偏移最小所有负荷节点的电压偏差绝对值之和 V_ref 1.0; % 标幺值基准电压 V_deviation sum(abs(V - V_ref)); obj2 V_deviation / length(V); % 取平均消除节点数影响 % 罚函数处理电压越限惩罚 V_min 0.95; V_max 1.05; penalty 0; for i 1:length(V) if V(i) V_min penalty penalty (V_min - V(i))^2; elseif V(i) V_max penalty penalty (V(i) - V_max)^2; end end obj1 obj1 10000 * penalty; obj2 obj2 10000 * penalty; end电压偏移的目标函数我特意用了“平均值”而不是“最大值”是有讲究的。如果目标是最大电压偏差最小算法会优先照顾电压最差的节点这看起来公平但实际上会让算法过度关注极端节点而在整体范围内搜索不足。用平均值的话算法会综合考虑所有节点的电压质量同时罚函数保证了极端节点也不会被忽视——这是一种平衡策略。4.3 目标函数的标准化处理多目标优化里还有一个关键细节——目标函数之间的量纲差异。网损的单位是kW或者标幺值功率数值通常在几十到几百的范围内但电压偏移是标幺值数值在0.001到0.05之间。如果直接拿这两个数去做支配判断电压偏移就会被网损“淹没”——网损稍微变一下电压偏移的变化根本看不出来支配关系完全由网损决定多目标优化退化成了单目标优化。解决办法是对目标函数进行标准化。最常用的是min-max归一化。在外部档案更新过程中记录每个目标的历史最小值和最大值然后对每个解的目标值做归一化处理。% 目标函数归一化处理 function [obj1_norm, obj2_norm] normalize_objectives(obj1, obj2, min_max_records) obj1_min min_max_records.obj1_min; obj1_max min_max_records.obj1_max; obj2_min min_max_records.obj2_min; obj2_max min_max_records.obj2_max; obj1_norm (obj1 - obj1_min) / (obj1_max - obj1_min 1e-6); obj2_norm (obj2 - obj2_min) / (obj2_max - obj2_min 1e-6); end注意加了一个1e-6的常数防止分母为零。这个标准化操作直接影响Pareto前沿分布的均匀性也是很多复现论文的代码和原文结果对不上的重要原因之一。5. 主程序循环与算法流程整合5.1 主算法迭代框架把前面所有模块组装起来形成完整的主程序循环。我写的是经典的小生境多目标PSO主循环参考了多篇文献的实现思路并根据实际测试做了调整% 参数设置 pop_size 100; % 种群规模 max_iter 200; % 最大迭代次数 c1 1.5; c2 1.5; % 学习因子 w_max 0.9; w_min 0.4; % 惯性权重范围 niche_radius 0.35; % 小生境半径 archive_size 50; % 外部档案最大容量 % 初始化 particles init_particles(pop_size, VarMin, VarMax); for i 1:pop_size [V, Ploss, ~] power_flow(bus, branch, particles(i, 1:4), particles(i, 5:8), particles(i, 9:10)); [Obj1(i), Obj2(i)] objective_func(V, Ploss); pbest(i, :) particles(i, :); pbest_obj1(i) Obj1(i); pbest_obj2(i) Obj2(i); end [pgbest, gbest_obj1, gbest_obj2] get_global_best(particles, Obj1, Obj2, archive); % 迭代主循环 for iter 1:max_iter w w_max - (w_max - w_min) * iter / max_iter; % 小生境划分每5代重新划分一次降低计算量 if mod(iter, 5) 1 [niche_idx, ~] niche_partition([Obj1, Obj2], particles, niche_radius); end % 更新粒子速度和位置 for i 1:pop_size v(i, :) w * v(i, :) ... c1 * rand(1, nVar) .* (pbest(i, :) - particles(i, :)) ... c2 * rand(1, nVar) .* (pgbest(i, :) - particles(i, :)); % 速度限幅 v(i, :) max(min(v(i, :), Vmax), Vmin); particles(i, :) particles(i, :) v(i, :); % 位置边界处理 particles(i, :) max(min(particles(i, :), VarMax), VarMin); % 约束修正 particles(i, :) constraint_handle(particles(i, :)); end % 变异操作防止早熟收敛 mutation_rate 0.1 - 0.07 * iter / max_iter; for i 1:pop_size if rand mutation_rate k randi(nVar); particles(i, k) VarMin(k) rand * (VarMax(k) - VarMin(k)); particles(i, :) constraint_handle(particles(i, :)); end end % 计算适应度 for i 1:pop_size [V, Ploss, success] power_flow(bus, branch, particles(i, 1:4), particles(i, 5:8), particles(i, 9:10)); if ~success Obj1(i) 1e10; Obj2(i) 1e10; continue; end [Obj1(i), Obj2(i)] objective_func(V, Ploss); end % 更新pbest和外部档案 for i 1:pop_size if dominates([Obj1(i), Obj2(i)], [pbest_obj1(i), pbest_obj2(i)]) pbest(i, :) particles(i, :); pbest_obj1(i) Obj1(i); pbest_obj2(i) Obj2(i); end end % 更新外部档案含小生境分布控制 archive external_archive_update(archive, [particles, Obj1, Obj2], archive_size); % 从小生境局部精英集中选择gbest pgbest select_gbest_from_niche(niche_idx, archive, pop_size); % 输出迭代进度 if mod(iter, 20) 0 fprintf(迭代次数: %d/%d, 外部档案解个数: %d\n, iter, max_iter, size(archive, 1)); end end5.2 变异操作的关键——自适应变异率很多多目标PSO的实现里变异率是固定的但我实测下来固定变异率效果很差。前期如果变异率太高算法发散好不容易找到的优质区域被变异粒子给破坏了后期如果变异率太低算法早熟所有粒子挤在一个局部Pareto前沿上出不来。我用的方案是变异率从0.1线性衰减到0.03。前期允许较大的变异概率来增加种群多样性后期减少变异让种群精细搜索当前区域。变异操作本身只随机改变一个维度这样既不会破坏太多已经找到的好结构又能提供潜在的跳跃能力。5.3 存档更新与拥挤距离计算这里单独说一下拥挤距离的实现因为这个细节直接影响最终Pareto前沿的质量。拥挤距离的计算思路是对一个粒子的某个目标值找到比它大的最近解和比它小的最近解计算这两个解在该目标上的差值然后除以该目标的最大值范围做归一化。把所有目标的归一化差值加起来就是这个粒子的拥挤距离。function cd crowding_distance(obj_matrix) % obj_matrix: N x nObj 的目标函数值矩阵 N size(obj_matrix, 1); nObj size(obj_matrix, 2); cd zeros(N, 1); for j 1:nObj [sorted_val, idx] sort(obj_matrix(:, j)); cd(idx(1)) inf; cd(idx(N)) inf; obj_range sorted_val(N) - sorted_val(1); if obj_range 1e-10 continue; end for i 2:N-1 cd(idx(i)) cd(idx(i)) (sorted_val(i1) - sorted_val(i-1)) / obj_range; end end end边界解的拥挤距离设为无穷大这是有意为之——Pareto前沿的两端是探索新区域的桥头堡必须优先保留。6. 仿真结果与Pareto前沿分析6.1 仿真环境与参数配置我的仿真环境是Matlab R2021bIntel i7处理器16GB内存运行时间大约2分钟到3分钟。系统采用IEEE 33节点配电网基准电压12.66kV总负荷有功3715kW无功2300kvar。DG配置如下节点18接入一个额定容量500kW的分布式电源节点22接入400kW节点25接入300kW节点33接入400kW。所有DG都具备无功调节能力逆变器剩余容量可以发无功无功调节范围设定为其视在功率的-40%到40%。此外在节点18和33各配置了一组电容器单组容量50kvar可以多档位投切总共最多10组并联。这样的配置和很多文献里的算例是一致的方便你做结果对照。6.2 小生境PSO与标准PSO的结果对比为了验证小生境PSO的效果我在同一仿真条件下跑了标准多目标PSO作为对照组两种算法共用相同的种群规模、迭代次数和目标函数。先说网损结果标准PSO找到的最优网损是151.2kW小生境PSO找到的最优网损是145.8kW降低了约3.6%。这个数值看起来不大但放在配电网优化里已经是不错的提升了而且关键在于两种算法的“搜索路线”完全不同。更显著的差异在Pareto前沿覆盖度上。标准PSO的Pareto前沿集中在网损目标较好的区域电压偏移目标的变化范围只有0.01左右而小生境PSO的Pareto前沿从低网损到低电压偏移完整铺开电压偏移的变化范围达到0.04以上。这意味着小生境PSO给决策者提供了更多样的选择——如果你想优先保证电压质量可以选择牺牲一部分网损的解如果更看重经济性就选网损最小的解。6.3 Pareto前沿的分布质量分析我拿Pareto前沿的分布质量做了一个量化评估用的指标包括C指标覆盖率衡量集合A的解有多少被集合B的解支配。标准PSO的非支配解集合被小生境PSO的解支配的比例达到83%说明小生境PSO寻找的Pareto前沿在整体上占优。间距指标Spacing衡量Pareto前沿上解之间的均匀程度。间距越小说明解在前沿上分布越均匀决策者对中间区域的把控越精准。小生境PSO的间距为0.023标准PSO为0.047前者几乎是后者的一半分布均匀性大幅提升。最大覆盖范围Maximum Spread衡量Pareto前沿在两个目标上的延展范围。小生境PSO的最大覆盖范围为0.87标准PSO为0.62。这直接说明小生境机制有效防止了种群过度收敛保持了更好的多样性。6.4 从Pareto解集中选取折中解工程实际里最终只需要一个调度方案不是一整条前沿。那怎么从Pareto解集里选出最终的折中解最常用的方法是模糊隶属度法。对于每个Pareto解在每个目标上计算模糊隶属度把所有目标的隶属度求和选择总和最大的解作为折中解。简单来说就是找那个“在两方面都不算太差的解”。我选取的折中解对应的控制策略是DG1出力498kW、DG2出力386kW、DG3出力295kW、DG4出力399kW总DG出力1578kW。电容器组1投入8组400kvar电容器组2投入7组350kvar。在这个方案下系统网损150.3kW平均电压偏移0.018所有节点电压都在0.97p.u.到1.03p.u.之间全部满足要求。对比不进行优化时的初始状态——网损202.5kW平均电压偏移0.035部分末端节点电压低于0.95p.u.——优化效果非常明显。7. 常见问题与调试心得7.1 粒子群不收敛或前期发散怎么办这个问题是最常见的。如果你的算法前期发散大概率不是算法框架的问题而是参数设置的问题。先检查惯性权重w和学习因子c1、c2。w过大粒子飞行太快容易越过最优区域c1和c2设置不当会导致粒子过度跟随个体经验或群体经验缺乏自主探索能力。我的建议是w从上界0.9线性衰减到下界0.4c1和c2都设为1.5这是最常用且经过大量验证的参数组合。再检查速度限幅Vmax。速度限幅的设置直接影响粒子在解空间中的移动步长如果Vmax太大粒子会剧烈震荡无法精细搜索如果太小粒子会爬行式前进收敛极慢。我建议Vmax设置为各维度变量范围的20%左右——对于有功出力在0到1之间的变量Vmax设为0.2对于电容器组数在0到10之间的整数变量Vmax设为2。还有一个容易被忽略的原因是潮流计算的收敛性。如果潮流计算函数存在bug导致很多粒子被误判为不可行解算法会在一个扭曲的适应度地貌上搜索自然无法收敛。排查方法是单独跑一次潮流把计算结果和MATPOWER或者公开文献里的结果对照确认无误后再接入优化算法。7.2 小生境划分过度导致收敛慢小生境半径选的太小粒子间信息交互大幅减少收敛会非常慢。我在调试时曾经遇到过运行500代都找不到满意解的情况最后发现是小生境半径只有0.1导致的。有一个实用的判断标准如果种群规模是100期望划分的小生境数量大概在8到12个之间。如果超过15个说明半径偏小如果少于5个说明半径偏大。先用两次预实验确定合适的半径范围再精细调整。另外要注意小生境的划分频率不需要每一代都做。粒子位置的更新是渐进的连续两代之间的小生境划分结果不会有太大变化每5代重新划分一次既能保证小生境信息不过期又能显著降低计算开销。7.3 外部档案中Pareto解数量不够怎么办外部档案里最终只有两三个解这种情况通常说明种群多样性严重不足。排除小生境机制的干扰后最常见的原因是变异机制失效。如果变异率太早衰减到零粒子很快就失去探索新区域的动力全部被吸附到当前最优解附近。我的做法是让变异率在最坏情况下也不低于0.02确保后期仍有小概率产生跳跃。另一个原因是外部档案的容量设置是否合理。如果容量太小比如只有10个即使算法发现了充足的非支配解档案更新机制也会把它们裁剪掉导致最终保留下来的解太少。建议至少设置为种群规模的50%。7.4 代码调试的三个黄金技巧写这种复杂的优化代码bug排查是很痛苦的过程我总结了一个三段式调试法第一阶段模块级验证。在把潮流计算、小生境划分、外部档案更新等模块接入主循环前先用构造好的小规模数据单独测试每个模块的输出是否符合预期。比如小生境划分模块我构造了10个已知坐标的粒子手动算好期望的划分结果再和程序输出比对。第二阶段单目标验证。把多目标问题退化为单目标问题只优化网损用标准PSO跑一遍把优化结果和文献值对照。如果单目标都跑不对说明算法实现有基础性错误多目标肯定更不对。第三阶段逐步恢复多目标。在单目标跑通的基础上逐步加入第二个目标、小生境机制、外部档案更新每加一个模块就跑一遍对比加模块前后的结果差异。这样能精确定位是哪个环节引入的问题。8. 后续扩展方向与几点经验总结代码写到这里主体功能已经完整。如果你的课题需要进一步扩展这里有几个我非常推荐的方向。第一个扩展方向是接入时序特性。当前的优化是一个静态截面也就是在某一时刻的DG出力和负荷水平下求解最优。但实际上DG出力和负荷是随时间变化的更贴近工程实际的做法是做96个断面每15分钟一个的时序优化或者使用典型日场景进行多场景随机优化。扩展时主要改造目标函数和约束条件算法框架不用动。第二个扩展方向是多种DG类型的建模。光伏、风电、储能、微型燃气轮机的运行特性差异很大光伏的无功能力由逆变器决定风电要考虑到风速的随机性储能还要考虑充放电状态和SOC约束。目前的决策变量编码方式有很强的通用性加变量类型只需扩展编码维度就行。第三个扩展方向是分布式求解。如果配电网规模从33节点扩展到几百个节点的实际馈线集中式PSO的计算负担会显著增加。可以考虑用分布式优化框架把网络分区每个分区跑一个独立的小生境PSO分区之间只交换边界信息。这和很多主动配电网分层分布式的理念完全契合。最后说几句掏心窝的总结。做智能算法和配电网优化结合的研究最大的误区是把精力全花在算法创新上把配电网模型做成“搭积木”——只是为了跑通算法而随便应付。我过去吃过亏前前后后调了大半年算法结果换了另外一个配电网模型来测试之前的很多经验全废了。后来我才意识到配电网物理约束的准确性才是决定优化结果能否落地的那条腿。真正好的研究工作一定是两条腿走路——算法有创新模型有深度两者互为表里才能产出既好看又实用的成果。这套小生境粒子群版本的协调优化代码在IEEE 33节点系统上的仿真结果表明它确实能在网损和电压质量之间找到一组分布均匀、覆盖面广的Pareto最优解。希望这篇文章里的思路和踩过的坑能让你在自己的课题里少走几步弯路。