简介基于MATLAB的粒子群优化在电力系统无功优化中的应用代码包面向电气工程研究生、电力系统分析相关技术人员以及需要解决无功补偿配置问题的工程实践者。核心代码为单个m脚本完整实现了粒子群算法流程包含种群规模、速度上下限、位置边界等参数定义适应度函数围绕电压稳定性与网络损耗最小化两个目标设计迭代过程中粒子依据个体经验和全局最优经验持续更新速度和位置直至满足收敛条件可辅助理解全局优化算法与最优潮流计算的结合方式。代码注释清晰结构精简压缩包仅含1个m文件大小约1KB便于直接导入MATLAB运行或对参数、目标函数做二次改造。已有374人学习下载适合正在学习智能优化算法、需要搭建无功优化仿真框架的读者作入门参考也可作为电力系统课程设计或学术研究的算法基线实现。1. 粒子群无功优化这名字听着刻板跑通之后是真有用粒子群无功优化说直白点就是让 MATLAB 里的粒子群算法PSO去替一个电力系统重新分配“无功”。电网电压偏低、网损偏高很多时候不是线路或变压器硬件出了问题而是发电机无功出力、变压器分接头档位、并联电容器组投切这三样东西配合没到位。无功优化就是把这种配合变成一个带约束的数学问题目标函数是网损最小或电压偏差最小约束是潮流方程、节点电压上下限和新能源与调压设备各自的调节范围。PSO 的好处在于它对变量类型不挑食——连续变量和离散档位能一起搜不用算海森矩阵不容易陷在初始点上。这套粒子群无功优化资源实际就是一份 MATLAB 工程包Matpower 做潮流计算自写 PSO 做搜索拿 IEEE 14 或 30 节点算例跑一遍几十行代码就能看到网损下降和节点电压整体抬升。适合正在做电力系统分析课程设计、本科毕业设计或者刚进课题组需要一套基准算法的同学。2. 无功优化的数学模型先搞清楚在算什么再让粒子去飞2.1 目标函数与状态变量的理论边界无功优化在前人的文献里已经写成很标准的范式但第一次接触的人很容易把控制变量和状态变量搅在一起。简单拆一下节点电压幅值 V 和相角 θ 是状态变量由潮流计算解得控制变量是我们能动的设备包括发电机无功出力 QG、变压器分接头变比 tap、并联电容器无功投切量 QC。目标函数一般写成网损最小min P_loss Σ P_ij(θ, V)其中 P_ij 是节点 i 到节点 j 的线路有功损耗它本质上是电压幅值和相角的函数。之所以说这个优化是混合整数非线性规划是因为 tap 和 QC 是离散量——变压器分接头按档位调电容器组按“投/切”决定投几组粒子连续搜索出来的位置不能直接用需要映射到最近的有效档位。发电机无功 QG 虽然是连续量但又有上下限。这种连续和离散混在一起的可行域用梯度类算法很别扭因为一旦对离散变量求导雅可比矩阵就失去意义。PSO 这类元启发式算法没有这个问题它只关心目标函数输入输出的数值关系离散变量通过取整或舍入映射来适配。这里有一个容易想偏的点不是所有节点电压都能自由优化。平衡节点的电压幅值和相角由潮流计算兜底PV 节点的电压幅值由发电机无功出力支撑PQ 节点的电压才完全受潮流方程支配。写目标函数和约束时这几类节点要分开处理否则罚函数会惩罚到不该惩罚的变量让粒子往错误方向飞。2.2 约束条件的矩阵化表达与罚函数无功优化的约束一般分三块。等式约束是潮流方程本身这部分由 Matpower 的 runpf 自动满足只要潮流收敛等式约束就成立。不等式约束包括节点电压幅值上下限、发电机无功上下限、变压器分接头档位范围、电容器组投切组数范围。把不等式约束写进粒子群搜索的常见做法有两种一是把约束作为边界限制粒子位置更新后直接回拉到合法区间二是把越限程度叠加到目标函数上构成罚函数。我一般在初版代码里用第二种因为返工少。罚函数的形式可以这样构造penalty k_V * sum(max(0, V - Vmax).^2 max(0, Vmin - V).^2) k_Q * sum(max(0, QG - QGmax).^2 max(0, QGmin - QG).^2)其中 k_V 和 k_Q 是两个惩罚系数取值通常在 1000 到 10000 之间视目标函数量级而定。为什么取这么大的数因为网损的量级可能只有几兆瓦而电压越限 0.05标幺值的平方是 0.0025如果不放大惩罚系数粒子会发现越限带来的“收益”比约束惩罚还大最后收敛出来的结果电压普遍越界。这个系数可以从 1000 开始试观察收敛后的最大电压越限量如果还有越限就继续加大直到越限量为零。注意罚函数是给 PSO 当“向导”用的不是潮流计算的一部分。计算适应度时先跑一次潮流如果潮流不收敛直接给一个大值比如 1e6省去后续所有罚函数计算。2.3 Matpower 的 loadcase 数据结构与潮流接口这套代码离不开 Matpower它不是独立写的潮流求解器。Matpower 的数据结构是三个核心矩阵mpc.bus、mpc.gen、mpc.branch。bus 矩阵第 8 列存节点电压幅值初值第 9 列存电压相角初值VM 列第 8 列是初始电压猜测值gen 矩阵第 2 列是有功出力 PG第 3 列是 QG第 4 和第 5 列是 QG 的上下限branch 矩阵从第 6 列到第 9 列是变压器的变比和相移参数其中第 9 列是变比 tap 的位置。在写 PSO 之前得先确认 mpc 里哪些发电机节点是要优化的。方法很简单看 gen 矩阵里 QGmax - QGmin 大于某个阈值比如 0.1 Mvar的发电机它们有调节空间才作为控制变量。那些 QGmax 和 QGmin 几乎相等的发电机说明它是固定无功运行不能动。同理branch 里只有存在变压器TAP 列非零的线路才需要把 tap 纳入粒子编码普通线路没有可调分接头不需要编码进去。Matpower 跑潮流的接口是mpopt mpoption(verbose, 0, out.all, 0, pf.alg, NR); results runpf(mpc, mpopt);这段代码的含义mpoption 创建一组潮流计算选项verbose 置 0 表示不打印逐次迭代日志out.all 置 0 表示不输出详细结果文件pf.alg 指定用牛顿-拉夫逊法。runpf 是 Matpower 的潮流计算入口传入 mpc 和 mpopt返回 results 结构。results.bus 的第 8 列是收敛后的节点电压幅值results.branch 的第 14 列是线路有功损耗加起来就是总网损。每次粒子更新位置后要改 mpc 里对应控制变量的值重新调 runpf才能得到这一组控制变量下的潮流结果和网损。有一点实测经验值得提Matpower 的版本不同runpf 返回结果的结构略有差异。老版本里 results.branch 第 14 列是损耗有些新版本里第 16 列才是。稳妥的做法是跑完用 size(results.branch, 2) 查看列数再确认列含义。这个检查习惯能省下后面一大半排查时间。3. 粒子群算法落地编码、适应度与边界处理3.1 粒子的编码方式与解码过程PSO 里一个粒子就是一组控制变量的集合。假设系统里有 3 台可调发电机、2 台可调变压器、2 组并联电容器那一个粒子的位置向量就是 7 维x [QG1, QG2, QG3, tap1, tap2, QC1, QC2]粒子位置向量的每一位都连续取值但 tap 和 QC 实际是离散的这就需要在计算适应度之前加一个解码步骤。常见做法是就近取整到合法档位。比如变压器分接头范围是 0.9 到 1.1步长 0.025那么连续粒子值 1.0337 应该被解码成 1.025因为 1.025 离它更近。电容器组如果最多投 4 组每组 5 Mvar那么连续值需要先从 [0, 4] 取整再乘以 5 得到实际无功。解码这一步必须放在修改 mpc 之前不能把连续值直接写进 mpc。原因很简单Matpower 的潮流计算要求 tap 必须是有效档位值如果填入 1.0337 这种数变压器模型会认为变比是连续可调算出来的结果虽然能收敛但不符合实际设备约束而且最优解往往是那些“量变不连续”的边界点取整后性能对不上。解码函数的核心逻辑大致是这样的function vars decode(x, gen_idx, tap_idx, cap_idx, step) vars x; if ~isempty(tap_idx) base vars(tap_idx); % 粒子给出的连续变比 tap_levels 0.9 : 0.025 : 1.1; % 分接头候选档位 [~, nearest] min(abs(tap_levels - base)); % 找最近档位 vars(tap_idx) tap_levels(nearest); % 回写解码结果 end if ~isempty(cap_idx) vars(cap_idx) round(x(cap_idx)); % 电容器组数取整 vars(cap_idx) vars(cap_idx) * step; % 组数乘单组容量 end end这个函数里 tap_idx 是粒子位置中变压器分接头对应的索引cap_idx 是电容器组对应的索引step 是单组电容器容量。min 函数配合 abs 距离找到距离当前值最近的候选档位round 把电容器连续值取整到最近整数。逻辑上解码不影响粒子的原始位置每次更新后重新解码保证粒子在连续空间里自由飞只在评估时才约束成实际设备值。这种“连续搜索、离散评估”的方式是元启发式算法处理混合整数问题的通用做法比直接在离散空间里穷举省事得多。边界的问题也要在这里处理好。解码后必须检查每个变量是否还在合法范围内因为 decode 只负责离散化不负责越界判断。比如粒子位置可能飞成 1.3这已经超过变压器档位上限了。我习惯在解码函数里把所有变量统一做一次 min/max 截断而不是放在主循环里逐个判断。3.2 适应度函数网损与罚函数的拼接策略适应度函数是整个 PSO 搜索的“裁判”它的设计直接决定粒子往哪飞。一个工程上验证可用的适应度函数长这样function fitness objective(x, mpc_base, mpopt, decode_info) mpc mpc_base; vars decode(x, decode_info.gen_idx, decode_info.tap_idx, decode_info.cap_idx, decode_info.cap_step); mpc.gen(decode_info.gen_idx, 3) vars(decode_info.gen_idx); % 写入发电机无功 mpc.branch(decode_info.tap_idx, 9) vars(decode_info.tap_idx); mpc.bus(decode_info.bus_shunt_idx, 5) mpc.bus(decode_info.bus_shunt_idx, 5) vars(decode_info.cap_idx); % 并联无功 try results runpf(mpc, mpopt); if ~results.success fitness 1e6; % 潮流不收敛直接给大罚值 return; end catch fitness 1e6; return; end Ploss sum(results.branch(:, 14)); % 线路有功损耗之和 V results.bus(:, 8); % 收敛后的节点电压幅值 QG results.gen(:, 3); % 收敛后的发电机无功 V_up max(0, V - 1.05); % 电压上越限量 V_down max(0, 0.95 - V); % 电压下越限量 Q_up max(0, QG - mpc.gen(:, 5)); % 无功上越限量 Q_down max(0, mpc.gen(:, 4) - QG); % 无功下越限量 penalty 10000 * (sum(V_up.^2) sum(V_down.^2)) 500 * (sum(Q_up.^2) sum(Q_down.^2)); fitness Ploss penalty; end这段代码的关键点有三个。第一写 QG 时用的是 mpc.gen 第 3 列对应 Matpower 数据格式里的 QG 字段而 QG 上下限在第 4 和第 5 列不能写反。第二电容器组并联到节点时修改的是 bus 矩阵第 5 列并联电纳初值要取原系统已有的并联值再加上新投入量直接覆盖会丢掉原始并联参数。第三罚函数里电压罚系数取 10000无功罚系数取 500这个比例是我在多个 IEEE 算例上试出来的经验值——电压越限对系统危害远大于无功越限所以加权更大。如果只关心网损不管电压可以把惩罚系数调小但那样收敛出来的结果很可能低压侧电压在 0.90 以下显然是不行的。3.3 速度-位置更新公式的参数含义与越限处理PSO 的标准更新公式是这个资源的核心引擎MATLAB 里的实现如下w 0.9 - 0.5 * iter / maxIter; % 惯性权重线性递减 v w .* v ... c1 * rand(size(v)) .* (pbest - x) ... c2 * rand(size(v)) .* (gbest - x); v max(vmax_min, min(vmax_max, v)); % 速度截断 x x v; x min(ub, max(lb, x)); % 位置越限回拉参数含义分别说明w 是惯性权重控制当前速度对下一步的影响程度初始 0.9、终值 0.4 是工程上最常用的区间前半段大惯性保证全局探索后半段小惯性利于局部精化。c1 是认知学习因子让粒子朝自己的历史最优飞c2 是社会学习因子让粒子朝全局最优飞两者通常都设 2。rand(size(v)) 产生与速度矩阵等大、0 到 1 之间的随机数矩阵给更新公式引入随机性。越限处理有两种流派截断法与反射法。截断法是直接 min(ub, max(lb, x)) 把位置拉回边界简单但会让粒子“堆积”在边界上如果最优解正好在边界附近且边界有平坦区粒子容易卡住。反射法是如果某个变量越过上界就把它从边界向内部反弹一段距里比如 x(i) 超出 ub则 x(i) 2*ub - x(i)这样缓冲区的粒子还能继续参与搜索。我在这套包里默认用截断法因为反射法在高维问题上要额外处理多变量同时越界的情况代码容易出边界 bug。速度上限 vmax 也很关键。vmax 设太小粒子搜索范围受限可能绕过最优解设太大会导致粒子飞过整个搜索空间方向和位置都失去意义。经验是 vmax 取每个变量搜索范围的 10% 左右。比如 QG 范围是 [-50, 50] Mvarvmax 就设 10tap 范围宽度是 0.2vmax 就设 0.02。在代码里 vmax 应该是一个与粒子维度同长的向量而不是单一标量因为不同变量量纲差异很大。4. 代码框架逐块拆解主程序、参数表与运行结果解读4.1 主程序结构与调用关系这套代码的目录结构很清爽四个文件各司其职main.m 负责初始化、循环、画图objective.m 是计算适应度的独立函数decode.m 负责离散变量解码case30.m 是 Matpower 自带的 IEEE 30 节点数据文件。实际项目中我建议把粒子群参数集中放在 main.m 开头改起来方便。%% main.m 主循环骨架 clc; clear; close all; rng(42); % 固定随机种子保证可复现 mpc loadcase(case30); % 载入IEEE 30节点系统 mpopt mpoption(verbose, 0, out.all, 0, pf.alg, NR); % 粒子群参数 nPop 50; % 种群大小 maxIter 100; % 最大迭代次数 dim 12; % 控制变量维度 lb zeros(1, dim); % 变量下界 ub ones(1, dim); % 变量上界后面按实际映射 vmax 0.1 * (ub - lb); % 速度上限 % 初始化种群 x repmat(lb, nPop, 1) rand(nPop, dim) .* repmat(ub - lb, nPop, 1); v zeros(nPop, dim); pbest x; gbest x(1, :); % 计算初始适应度 for i 1:nPop fitness(i) objective(x(i, :), mpc, mpopt, decode_info); if fitness(i) fitness(gbest_idx) gbest_idx i; end end pbest_fit fitness; ...这个骨架里 rng(42) 是很多人忽略但极关键的一行。PSO 的随机初始化意味着每次运行结果都不同如果不固定种子论文或报告里的数据就无法复现答辩时导师让你再跑一遍结果对不上非常尴尬。rng 函数的参数可以是任意整数42 只是随手写的重要的是它让后续所有 rand 调用都生成同一序列的随机数。loadcase(case30) 从 Matpower 的数据目录加载 IEEE 30 节点系统。dim 的取值不是随便定的它等于可调发电机数 可调变压器数 电容器组数具体多少个要打开 case30.m 数一下不要照抄我这里的 12实际系统不同维度不同。lb 和 ub 这一段我在骨架里用零一区间示意真实代码里要根据每个变量的实际范围组装成向量否则粒子初始化就全乱套了。4.2 参数设置参考表与推荐取值整理一份在 IEEE 14、30、57 节点系统上实测过有效的参数表新手可以直接照抄快点点调参数名称推荐取值调整方向说明种群规模 nPop50 100节点数多或变量维度高时取大值迭代次数 maxIter100 200看迭代曲线是否在后期仍有下降惯性权重初值 w_start0.9大于 0.95 时收敛慢震荡严重惯性权重终值 w_end0.4小于 0.3 时容易早熟停滞学习因子 c12.0过大会过度依赖自身历史学习因子 c22.0过大会过早聚集到 gbest速度上限系数0.1 × 变量范围过大导致位置飞散不收敛电压罚因子 k_V10000有电压越限就加大无功罚因子 k_Q500有 QG 越限就加大这套参数的一个替换逻辑是如果看到迭代曲线在整个迭代区间里一直平稳下降、没有平台期说明迭代次数偏少可以从 100 调到 150 或者 200如果曲线在 30 代就完全走平、后面 70 代毫无变化说明种群多样性不足优先把 nPop 提到 80 或 100而不是去动 c1/c2。这个判断依据是惯性权重线性递减到 0.4 后搜索步伐已经很小如果此时还没有找到好的区域加大迭代次数只会浪费算力。还有一个容易忽略的参数是 Marxer——不是 Matpower 都有我自己比较常用的是把 pf.tol潮流收敛容差设置为 1e-8。容差太大会导致不同的粒子算出来的网损精度不一致迭代曲线会出现毛刺容差太小会增加每次潮流的计算量对 30 节点这种规模影响不大但对 118 节点就会明显变慢。4.3 迭代曲线与优化结果的判读方式跑完主循环后代码末尾会画出一条网损随代数下降的曲线。判读这条曲线时不要只看最后收敛到什么值要看它的下降形态。标准的健康形态是前 20% 的迭代区间内网损快速下降中间 50% 缓慢下降且伴随小幅波动最后 30% 基本平缓、波动幅度很小。如果曲线在前几代就直线跌到底后面完全平走说明初始种群质量太高或粒子太快聚集到局部最优可以重新调整随机种子或减小 w_start。如果整条曲线都在锯齿形震荡说明 w 衰减太快、后期惯性仍然偏大尝试把 w_end 从 0.4 降到 0.3。结果判读除了看网损数字还要看电压分布。优化后最直观的验证方式是分别打印优化前和优化后的电压最低点fprintf(优化前网损: %.4f MW最低电压: %.4f p.u.\n, base_loss, min(base_V)); fprintf(优化后网损: %.4f MW最低电压: %.4f p.u.\n, best_loss, min(best_V));如果优化后网损降了但最低电压反而更低说明罚函数的电压权重不足粒子牺牲电压去换网损这在工程上是不可接受的。电压是电能质量的核心指标网损降低的收益在幅度上远不如一台变压器普遍处于欠电压状态下运行带来的隐患大。这时要果断把 k_V 从 10000 提到 20000 或 50000重新运行。反过来如果网损没怎么变但电压被抬得很高说明罚函数权重过大粒子把主要精力放在调整电压上但网损搜索不够充分。5. 常见问题排查收敛慢、越限、结果不稳定的四个坑5.1 现象一runpf 报错或收敛标志为假现象运行过程中某个粒子的潮流计算失败程序直接中断或在结果里看到大量 success0 的粒子。原因粒子初始化时某些控制变量取了极端组合导致潮流方程无解。最常见的触发点是发电机无功被取到接近下限而系统有功负载又高节点电压被拉垮牛顿-拉夫逊法迭代发散。解决第一个办法是初始化时就把粒子位置限制在变量的实际可行区间内而不是用全范围随机。第二个办法是在 objective.m 里捕获 runpf 异常遇到不收敛直接返回 1e6不中断主程序。第三个办法是检查 mpc 里 QGmax、QGmin 的单位是不是 MVAr与负荷单位 MVA 是否一致单位混乱会导致潮流误判。我自己被坑过最严重的一次是某个数据文件的 QGmin 写的是标幺值实际范围全部偏大粒子飞出去后一大批解都不可行。5.2 现象二迭代曲线震荡剧烈长时间不下降现象网损曲线不是平滑下降而是每次迭代上下跳动甚至有个别代网损比上一代高出一大截。原因惯性权重 w 起始值太高或者速度上限 vmax 太大。粒子以过高的速度飞越最优区域无法稳定落入极值点。另外罚函数系数如果突然让某个不可行解的罚值巨大也会导致曲线出现尖峰。解决把 w_start 从 0.95 降到 0.85同时把 vmax 缩小到变量范围的 5%。还有一个处理技巧是给 gbest 增加“停滞触发判断”——如果连续 10 代 gbest 没有变化说明粒子已经聚集此时把 w 调回 0.6用较大速度冲出局部区域重新搜索比坐等迭代次数耗尽有效得多。5.3 现象三每次运行结果差异很大复现困难现象同一套参数每次运行网损值可能差 0.5 MW 以上电压分布也不同。原因PSO 是随机算法初始化位置和速度由随机数决定不同运行会生成不同的初始种群。这个问题在种群规模较小时尤其明显比如 nPop20 时初始种群的随机性对结果影响非常大。解决在 main.m 第一行加 rng(固定种子)保证每次运行生成完全相同的随机序列。如果需要统计意义上的最优值就是固定种子跑一遍作为报告数据或者改用 10 个不同种子各跑一遍取最好值。注意 rng 必须放在初始化种群之前放在粒子循环内部没有意义。另外把 nPop 提高到 80 以上也能显著降低随机种子对结果的影响。5.4 现象四早熟收敛结果明显不如其他算法现象迭代曲线前 20 代就完全走平最终网损比文献里报道的高出一截且电压越限依旧存在。原因PSO 在搜索后期粒子多样性严重丢失所有个体都聚集到某个局部最优解附近无法跳出去。这种情况在变量维度高、可行域非凸时更容易发生。解决给粒子添加变异操作每代按 5% 的概率随机重置某个维度的取值或者采用多种群策略把种群分成 3 个子群每 10 代交换一次信息。这两种改法都会让代码复杂度上一个台阶但对 IEEE 57、118 节点这类中大规模算例来说是值得的。还有一种快速验证办法把 PSO 的跑结果和 Matpower 内置的 SQP 优化器结果做对比如果 PSO 长期差于 SQP说明搜索机制需要改而不是参数没调好。6. 让结果更进一步自适应惯性权重与多种群协同参数表里固定 w 线性递减有两种改法实测有效的是按粒子聚集程度自适应调节。思路是每一代计算所有粒子的平均适应度与 gbest 适应度的相对差距差距大说明种群还没收敛w 保持较大值维持探索差距小说明粒子都快挤在一起了w 减小以精化搜索。实现上把原来那行 w 0.9 - 0.5 * iter / maxIter 替换成fitness_std std(fitness) / max(abs(gbest_fit), 1e-10); w 0.9 - 0.5 * fitness_std; w max(0.4, min(0.9, w));其中 std(fitness) 是所有粒子适应度的标准差如果所有粒子都聚集在同一个解周围标准差趋近于零w 就自动趋近 0.4如果粒子分散w 保持高位。这个做法的好处是不需要额外调迭代区间比例参数适应度分布本身就是种群多样性的天然指标。多种群协同的做法更简单粗暴把 nPop 分成 3 个子群每个子群独立演化 10 代后把所有子群的 gbest 放到一起比较把最差的子群淘汰用最好 gbest 附近的随机粒子补齐。这种做法对避免早熟收敛的效果很直接但会让目标函数调用次数增加大约 20%在 118 节点系统上计算时间会从几分钟涨到十几分钟需要评估你的算力耐受度。关于结果的验证方法我个人的习惯是用 PSO 的结果和 fmincon 里的 SQP 算法做一次交叉验算拿 PSO 优化得到的最优控制变量代入同样的 mpc用 runopf 跑一遍确认网损和电压都匹配避免 PSO 那边因为罚函数误差给出一个“看起来最优但实际不可行”的解。这个交叉验算不用写进代码包作为手工检查步骤跑一遍几十秒就够。回头说我自己最初改这个包的经历——第一次跑的时候没注意 rng结果前后两次运行网损差得离谱我还以为是代码逻辑写错了花了一晚上逐段排查最后才发现只是随机种子的问题。从那以后我每次跑新算例都强制走一遍固定种子、检查越限、验证潮流收敛三步流程再也没被随机性坑过。粒子群无功优化这种东西参数和代码都只是敲门砖真正的门槛在于你能不能稳定复现一个可信的结果。希望帮到你。本文还有配套的精品资源点击获取