简介：这份资源面向电力系统分析与智能优化算法的学习者，聚焦粒子群算法在无功优化中的工程实现。包内仅含1个m文件，压缩包约1KB，核心脚本围绕粒子群优化流程展开，涵盖种群规模、速度与位置上下限等参数设定，以及以电压稳定性和网络损耗最小化为目标的适应度函数设计，并涉及节点平衡方程、线路阻抗模型与无功设备约束等数学建模环节。读者可借助MATLAB环境直接运行与调试，理解粒子如何依据个体经验与全局最优经验迭代更新，逐步逼近满足运行约束的无功补偿方案，从而掌握最优潮流与无功配置的求解思路。目前已有374人学习，适合具备一定电力系统分析基础、希望用智能算法解决多目标非线性优化问题的读者参考。
1. 粒子群无功优化：为什么随机搜索能驯服电网的“电压焦虑”
电网里最让人头疼的问题之一，不是发电够不够，而是无功功率怎么分配。有功功率管的是“能不能干活”，无功功率管的是“电压稳不稳”。一台发电机发出的无功多了，节点电压就往上抬；少了，电压就往下掉。更麻烦的是，网损会随着无功潮流的分布而变化——同样一批负荷，无功分配方案不同，线路上烧掉的功率可能差出好几个百分点。传统做法靠灵敏度矩阵和线性规划，在局部线性化点附近找最优解，可一旦系统运行点偏移、负荷波动，线性假设就崩了。粒子群算法之所以被拿来做无功优化，恰恰因为它不依赖梯度信息，能在非凸、离散、多约束的空间里“撒一把粒子”去试探。这份资源把整套流程做成了一个可运行的 MATLAB 工程，从潮流计算到粒子更新再到约束处理，链路是完整的。搞懂它，等于搞懂了一类随机优化算法怎么跟电力系统分析对接。
2. 无功优化的数学本质：为什么梯度法会在这里翻车
2.1 目标函数不是凸的，梯度方向可能指向更差的解
无功优化的标准写法通常是：在满足潮流方程和运行约束的前提下，最小化网损，或者最小化电压偏差，或者两者加权。控制变量包括发电机端电压、变压器抽头、无功补偿装置容量。问题在于，潮流方程本身是非线性的，网损对控制变量的二阶导数矩阵并不保证正定。换句话说，目标函数的等高线不是一个个同心圆，而是像山谷里套着小山包。梯度法从某个初始点出发，沿着最陡下降方向走，很容易一头扎进局部最优——那个点的梯度确实是零，但全局看它只是个小坑。
常见的实现思路是：先用潮流计算求出当前控制变量下的状态变量，再通过数值扰动或灵敏度公式求梯度。但数值扰动的步长很难选，步长大了精度不够，步长小了被潮流收敛误差淹没。粒子群算法绕开了这个问题，它不关心梯度，只关心每个粒子当前位置的目标函数值。粒子之间共享历史最优信息，整个群体在解空间里既保留个体记忆，又向全局最优靠拢。这种机制天然适合处理无功优化这种“目标函数写得出解析式但求不出解析梯度”的场景。
2.2 离散变量和连续变量混在一起，传统优化器要分两层
变压器抽头是离散的，通常有 17 个档位或 19 个档位；无功补偿电容也是分组投切的，每组容量固定。发电机端电压是连续的。如果直接用连续优化算法，得到的抽头位置可能是 8.37 档，实际装置根本没法执行。传统做法是外层枚举离散变量、内层优化连续变量，计算量随离散变量数量指数增长。粒子群算法处理这个问题的方式很直接：把离散变量也编码进粒子位置，在更新后做取整或就近映射。虽然理论上不保证收敛到整数最优，但工程上通常能得到可执行的满意解。
注意：取整操作会破坏粒子速度的连续性，常见的补救是限制取整后粒子的速度更新幅度，或者在适应度函数里对偏离整数档位的行为加惩罚项。
2.3 约束处理是真正的分水岭
无功优化的约束分两类：等式约束是潮流方程，每个节点都要满足有功和无功平衡；不等式约束包括发电机无功出力上下限、节点电压上下限、线路传输功率上限。粒子群算法本身是无约束搜索，必须外挂约束处理机制。常见的做法有三种：罚函数法、可行解优先法、修复法。罚函数法把约束违反量乘以惩罚系数加到目标函数上，实现简单但惩罚系数难调；可行解优先法在比较粒子优劣时，先看是否可行，再看目标值；修复法直接把越界粒子拉回边界。这份资源里用的是罚函数加边界截断的组合，后面会展开说。
3. 粒子群算法的内核：速度更新公式到底在干什么
3.1 三个分量各自承担什么角色
粒子群的速度更新公式通常写成这样：
% 速度更新：惯性项 + 个体认知项 + 社会学习项 v(i,:) = w * v(i,:) ... + c1 * rand(1, nVar) .* (pbest(i,:) - x(i,:)) ... + c2 * rand(1, nVar) .* (gbest - x(i,:)); % 位置更新 x(i,:) = x(i,:) + v(i,:); % 边界处理：越界截断 x(i,:) = max(x(i,:), lb); x(i,:) = min(x(i,:), ub);第一项w * v是惯性项，保留粒子上一时刻的运动方向。w越大，粒子越“固执”，全局探索能力强但局部收敛慢；w越小，粒子越容易被拉向最优解，收敛快但可能早熟。第二项c1 * rand .* (pbest - x)是个体认知项，把粒子拉向它自己历史上到达过的最好位置。第三项c2 * rand .* (gbest - x)是社会学习项，把粒子拉向整个群体当前找到的最好位置。c1和c2分别控制“相信自己”和“相信群众”的权重，常见取值在 1.5 到 2.0 之间。
rand(1, nVar)是每个维度独立的随机数，这一点很关键。如果所有维度共用同一个随机数，粒子就只能沿着从当前位置指向 pbest 或 gbest 的直线方向运动，搜索轨迹会被限制在低维流形上。逐维随机化让粒子在每个维度上有独立的步长，搜索覆盖更充分。
3.2 惯性权重为什么要从大到小线性递减
固定惯性权重有个两难：前期需要大w来广泛探索，后期需要小w来精细收敛。常见的改进是线性递减：
% 惯性权重从 w_max 线性降到 w_min w = w_max - (w_max - w_min) * (iter / maxIter);iter是当前迭代次数，maxIter是最大迭代次数。典型取值w_max = 0.9，w_min = 0.4。这个策略背后的直觉是：搜索初期粒子应该像撒出去的沙子，散得越开越好；搜索后期应该像收网，逐步向最优区域聚拢。如果w一直很大，粒子会在最优解附近来回震荡，永远停不下来；如果w一直很小，粒子很快失去探索能力，整个群体挤在一个局部最优里出不来。
3.3 适应度函数怎么把潮流计算嵌进去
粒子群每迭代一次，每个粒子都要算一次适应度。在无功优化里，适应度函数的核心是调用潮流计算：
function fitness = evaluate_fitness(x, baseMVA, bus, branch, gen, load) % 把粒子位置映射到控制变量 gen(:, 6) = x(1:nGen); % 发电机端电压 tap = round(x(nGen+1:nGen+nTap)); % 变压器抽头取整 shunt = round(x(nGen+nTap+1:end)); % 补偿容量取整 % 调用潮流计算 [busOut, genOut, branchOut, success] = run_powerflow(baseMVA, bus, branch, gen); if ~success fitness = 1e6; % 潮流不收敛，给一个极大的惩罚值 return; end % 计算网损 loss = sum(real(get_branch_flows(branchOut))); % 计算电压越限惩罚 vViolation = sum(max(0, abs(busOut(:, 8)) - 1.05) + max(0, 0.95 - abs(busOut(:, 8)))); % 计算发电机无功越限惩罚 qViolation = sum(max(0, genOut(:, 3) - genOut(:, 9)) + max(0, genOut(:, 10) - genOut(:, 3))); % 综合适应度 penalty = 100 * vViolation + 100 * qViolation; fitness = loss + penalty; end这段代码有几个关键点。第一，潮流不收敛时直接给一个极大的适应度值，相当于把这个粒子判死刑，不让它参与后续的 pbest 和 gbest 竞争。第二，电压越限惩罚用的是abs(busOut(:, 8))，因为潮流计算输出的电压可能是标幺值，需要看幅值。第三，惩罚系数100不是随便定的，它要足够大，让可行解的适应度明显优于不可行解，但又不能大到让目标函数完全失去区分度。实际调参时，可以先跑几轮，观察可行解和不可行解的适应度分布，再决定惩罚系数。
4. 从粒子到电网：控制变量编码与潮流计算的对接细节
4.1 粒子维度怎么确定
假设系统有 5 台发电机、4 台有载调压变压器、3 组并联电容器。粒子位置向量的维度就是 5 + 4 + 3 = 12 维。前 5 维是发电机端电压，中间 4 维是变压器抽头档位，后 3 维是电容器投切组数。每一维的上下界由设备物理限制决定：发电机端电压通常在 0.95 到 1.05 标幺值之间，变压器抽头在 0.9 到 1.1 之间（对应不同档位），电容器组数在 0 到最大组数之间。
编码方式有两种：一种是直接编码物理量，比如电压值、抽头标幺值；另一种是编码档位索引，比如第 3 档、第 7 档。直接编码的好处是粒子位置直接对应控制变量，不需要额外映射；坏处是离散变量的取整操作会引入不连续。档位索引编码的好处是离散变量天然整数，坏处是连续变量和离散变量混在一起时，速度更新的物理意义不直观。这份资源用的是直接编码加取整，工程上更常见。
4.2 潮流计算是黑匣子，调用前要保证数据格式一致
MATLAB 里做潮流计算常用的工具是 MATPOWER 风格的数据结构。bus矩阵每一行是一个节点，列的含义包括节点编号、类型、有功负荷、无功负荷、电压幅值、相角等。gen矩阵每一行是一台发电机，列包括节点编号、有功出力、无功出力、电压设定值、无功上下限等。粒子群在更新位置后，要把新的控制变量写回gen和bus矩阵，再调用潮流计算。
这里有个容易翻车的地方：潮流计算会修改bus和gen矩阵的内容，如果每次评估都在原矩阵上改，下一次评估就会拿到被污染的数据。常见的做法是在评估函数内部先做深拷贝：
busCopy = bus; genCopy = gen; branchCopy = branch; % 然后在 busCopy、genCopy 上修改控制变量另一个坑是节点类型。平衡节点的电压和相角是给定的，不参与优化；PV 节点的有功出力和电压幅值给定，无功出力由潮流计算得出；PQ 节点的有功和无功负荷给定，电压幅值和相角待求。粒子群在编码时，只能把 PV 节点的电压设定值作为控制变量，不能把平衡节点的电压也放进去，否则潮流方程会过约束。
4.3 网损计算不能只看线路电阻
网损的精确计算需要用到线路两端的电压和相角：
% 计算线路损耗 for k = 1:nBranch from = branch(k, 1); to = branch(k, 2); Vf = busOut(from, 8) * exp(1j * busOut(from, 9) * pi / 180); Vt = busOut(to, 8) * exp(1j * busOut(to, 9) * pi / 180); If = (Vf - Vt) * branch(k, 4) + Vf * branch(k, 5); % 支路电流 It = (Vt - Vf) * branch(k, 4) + Vt * branch(k, 6); Sf = Vf * conj(If); % 首端视在功率 St = Vt * conj(It); % 末端视在功率 loss = loss + real(Sf + St); % 有功损耗 endbranch(k, 4)是串联导纳，branch(k, 5)和branch(k, 6)分别是首端和末端的并联导纳。如果只算I^2 * R，会忽略并联导纳上的损耗，在长线路或高压系统中误差不可忽略。另外，real(Sf + St)得到的是有功损耗，无功损耗是imag(Sf + St)，无功优化主要关心有功损耗，但无功损耗会影响电压分布，间接影响有功损耗。
5. 避坑与排查：粒子群无功优化里最容易翻车的几个地方
5.1 潮流不收敛导致整个种群“假死”
现象：迭代初期还有适应度值在变化，跑了几十代之后所有粒子的适应度都变成同一个大数，gbest 不再更新。
原因：某个粒子被随机初始化到一个极端位置，比如所有发电机电压都拉到上限、所有电容器都投入，导致潮流方程无解。这个粒子的适应度被赋值为 1e6，如果它恰好是初始 gbest，整个种群会被拉向这个不可行区域。更糟的是，如果惩罚系数设置不当，不可行解的适应度可能比某些可行解还小，粒子会主动往不可行区域跑。
解决：初始化时对每个粒子做一次潮流校验，不收敛的粒子重新随机生成。迭代过程中，如果连续多代 gbest 不变，检查当前 gbest 对应的控制变量是否导致潮流不收敛。如果是，把 gbest 重置为历史可行最优解。
5.2 离散变量取整后粒子“卡死”在边界
现象：变压器抽头或电容器组数在某一代之后不再变化，粒子位置在边界附近来回震荡。
原因：取整操作把连续位置映射到离散档位后，速度更新仍然基于连续位置计算。如果粒子位置在 8.4 和 8.6 之间震荡，取整后都是 8 档，但速度更新看到的是连续值的变化，可能把粒子推向 9 档或 7 档，下一轮又拉回来。这种震荡消耗迭代次数，但不产生有效搜索。
解决：对离散维度，在速度更新后直接做取整，然后用取整后的位置计算 pbest 和 gbest。或者对离散维度单独设置较小的速度上限，限制单步跳档幅度。
5.3 惩罚系数过大导致可行解“淹没”在噪声里
现象：所有粒子的适应度值都很大，可行解和不可行解的适应度差异被惩罚项主导，目标函数（网损）的优化方向丢失。
原因：惩罚系数设得太大，比如 1e4，导致任何微小的电压越限都会产生巨大的惩罚值。可行解之间的网损差异可能只有 0.01，但惩罚项的波动可能有 0.1，粒子群实际上在优化惩罚项而不是网损。
解决：惩罚系数应该与目标函数的量级匹配。网损通常在 0.01 到 0.1 标幺值之间，惩罚系数取 10 到 100 比较合适。更好的做法是采用自适应惩罚：迭代初期惩罚系数小一些，让粒子先探索可行域；迭代后期惩罚系数增大，迫使粒子留在可行域内。
5.4 种群多样性过早丧失
现象：迭代不到 50 代，所有粒子的位置就聚成一团，gbest 不再明显更新，但距离理论最优还有差距。
原因：c2设得太大，或者w下降太快，粒子过早向 gbest 靠拢。如果初始种群本身多样性就不够，比如所有粒子都在同一个局部区域随机生成，问题会更严重。
解决：初始化时用拉丁超立方采样或均匀设计，保证粒子在解空间里分散。c1和c2可以设成动态的：前期c1大一些鼓励个体探索，后期c2大一些鼓励群体收敛。另一种做法是引入变异算子，每隔若干代对部分粒子做随机扰动。
5.5 目标函数量纲不统一导致加权系数难调
现象：如果目标函数同时包含网损和电压偏差，网损的量级是 0.01 到 0.1，电压偏差的量级是 0.001 到 0.01，直接相加的话电压偏差项几乎不起作用。
原因：两个目标的物理量纲不同，数值范围差一个数量级。如果不做归一化，优化结果会偏向量级大的那个目标。
解决：对每个目标做归一化，比如除以初始值或基准值。或者用权重系数显式调整：fitness = alpha * loss / loss0 + beta * vDev / vDev0，其中loss0和vDev0是初始潮流下的值。alpha和beta根据运行需求调整，如果更关心电压质量，就增大beta。
6. 参数调优的取舍：没有免费午餐，只有场景匹配
粒子群算法在无功优化里的表现，很大程度上取决于参数怎么选。种群规模、最大迭代次数、惯性权重范围、学习因子、惩罚系数，这些参数之间不是独立的。种群规模大，每代计算量大，但需要的迭代次数少；种群规模小，每代快，但容易早熟。我一般会这么理解：如果系统规模小（比如 IEEE 14 节点），种群取 30 到 50 就够；如果系统规模大（比如 IEEE 118 节点），种群至少取 100，否则解空间覆盖不足。
最大迭代次数不是越多越好。无功优化的控制变量维度通常在 10 到 50 之间，粒子群在 200 到 500 代内通常能收敛到满意解。超过 500 代之后，gbest 的变化往往只是在小数点后第四位，工程上已经没有意义。与其把迭代次数设成 2000，不如多跑几次独立实验，看结果的方差。如果 10 次独立运行得到的网损差异在 1% 以内，说明算法稳定；如果差异超过 5%，说明参数需要调整，或者问题本身有多个局部最优。
惯性权重的线性递减策略在大多数场景下够用，但有一种情况需要特别注意：如果目标函数在迭代初期就快速下降，然后长时间停滞，可能是w下降太快，粒子过早失去探索能力。这时候可以试试非线性递减，比如用余弦函数或者指数函数，让w在前期保持较大值更长时间。
学习因子c1和c2的经典取值是 2.0 和 2.0，但这个组合偏向社会学习，容易早熟。我通常会把c1设成 2.5，c2设成 1.5，让粒子更相信自己的历史经验。如果问题对全局最优要求很高，可以反过来，c1取 1.5，c2取 2.5，加快群体收敛。还有一种做法是让c1从大到小、c2从小到大动态变化，模拟“先探索后开发”的过程。
惩罚系数是最需要根据实际系统调整的参数。如果系统电压水平本来就很好，电压越限很少发生，惩罚系数可以小一些，让网损优化占主导。如果系统电压经常越限，惩罚系数必须足够大，否则优化结果不可行。一个实用的技巧是：先跑一次不带惩罚的优化，看看电压越限的分布，再根据越限的严重程度设定惩罚系数。如果最大越限量是 0.02 标幺值，惩罚系数取 50 到 100 就能把越限压到 0.001 以内。
从那以后我每次做无功优化，都会先跑一次基础潮流，记录初始网损和电压分布，再根据越限情况设定惩罚系数，最后用 10 次独立运行验证稳定性。这套流程走下来，基本不会出现“跑了一晚上结果不可行”的情况。希望帮到你。
本文还有配套的精品资源，点击获取
