基于多目标粒子群优化算法的计及光伏波动性的主动配电网有功无功协调优化（Matlab代码实现）
搞配电网优化这块的同行，对“主动配电网”这个词应该都不陌生。传统配电网是被动接收电能，但随着分布式光伏大规模接入，电压越限、潮流倒送、网损攀升这些问题变得很棘手，配电网必须从“被动”转向“主动”，也就是要主动去管理配电网内部的分布式电源、储能、无功补偿装置等资源。而我今天要分享的这个项目，就是围绕主动配电网调度中的一个核心难题——有功无功协调优化，在计及光伏波动性这一实际约束下，采用多目标粒子群优化算法（MOPSO）进行求解，并给出完整的Matlab实现方案。
这段时间我在做相关课题，发现很多人一上来就套标准PSO，忽略了一个关键点：光伏出力是波动的，不是恒定值。如果不把这个波动性建模进去，你算出来的调度策略可能在实际运行中根本不可行。而这个项目恰恰把光伏波动性作为约束条件之一纳入优化模型，同时用多目标框架处理网损和电压偏差这两个相互冲突的指标，思路很清晰，也非常贴合工程实际。这篇文章我就把整个项目的设计思路、数学模型、算法原理、Matlab实现框架以及我在实际调试中踩过的坑都整理出来，希望能帮到正在做配电网优化调度或者打算用多目标粒子群算法做研究的朋友。
1. 项目整体思路与研究价值
1.1 为什么一定要考虑光伏的波动性
光伏出力受光照强度、云层遮挡、温度等因素影响，具有明显的随机性和间歇性。很多研究者在做优化时，直接把某个时段的光伏出力当作固定值来处理，这种做法在理论上是简化了，但工程上往往站不住脚。
打个比方，就好比你规划一条公交线路，假设每个站点上下车人数是固定的，那班次安排特别容易；但真实情况是早晚高峰和平峰的人流量差异很大，你必须把这个波动考虑进去，否则要么高峰期运力不足、要么平峰期空跑浪费。光伏波动性的引入也是这个道理，它让配电网调度从“理想化的确定性优化”走向“贴近实际的鲁棒性优化”。
在本项目中，光伏波动性主要通过两方面体现：
在约束条件中，引入节点电压的随机波动区间，光伏出力的变化会改变潮流分布，导致节点电压随之波动，若优化结果没有预留足够的调节裕度，实际运行时极易出现电压越限。
在目标函数中体现为网损和电压偏差这两个指标的累积统计特性，不是在单一场景下优化，而是在典型光伏波动场景下综合优化。
项目中采用的常见做法是设定多个典型光伏出力场景，比如晴天平稳出力、多云波动出力、阴天低出力等，然后让优化算法在所有这些场景下都能给出可行的调度结果，这就叫“场景鲁棒”。
1.2 有功无功协调优化的本质
传统配电网调度中，有功和无功往往分开处理：有功靠发电机调整、无功靠电容器投切和有载调压变压器分接头调整。但在有高比例分布式光伏接入的主动配电网中，有功潮流和无功电压是强耦合的。光伏逆变器既可以发有功，也可以吸收或发出无功，所以有功和无功必须联合优化，才能达到更好的效果。
举个例子，某个节点光伏大发时，如果只看有功调度，只要保证光伏出力被消纳就行，但此时线路可能因为潮流倒送导致末端节点电压升高。如果这时候合理调节逆变器的无功输出，吸收一部分无功，就能把电压压住，同时还可能降低线路中的无功流动，减少网损。这就是有功无功协调优化的价值所在。
本项目建立的优化模型中，控制变量一般包括：
- 分布式光伏的有功削减量（如果允许）
- 光伏逆变器的无功输出量
- 无功补偿装置（如电容器组）的投入容量
- 有载调压变压器（OLTC）的分接头位置
目标函数则通常包括两个：系统网损最小，节点电压偏差最小。这两个目标在数学上往往是矛盾的：降低网损可能要求增加无功补偿，这会影响电压分布；而追求电压质量可能要以一定的网损为代价。所以这个问题天然适合用多目标优化算法求解，得到一组Pareto最优解，供调度人员根据实际需求灵活选择。
1.3 为什么选多目标粒子群优化算法
配电网有功无功协调优化是一个非线性的、多约束的、混合整数优化问题，传统数学规划方法在处理这类问题时比较吃力，尤其是控制变量中含有离散变量（分接头档位、电容器组投切组数）时，求解难度更大。
粒子群优化算法（PSO）的优势在于实现简单、不依赖目标函数的梯度信息、全局搜索能力强，而且对连续和离散混合变量有很好的适应性。多目标粒子群优化算法（MOPSO）则是在PSO基础上引入Pareto支配关系、外部档案集和拥挤距离机制，使得算法能够同时寻优多个相互冲突的目标，最终得到一组均匀分布的Pareto前沿解。
相比遗传算法（NSGA-II），MOPSO收敛速度更快，参数更少，在配电网这类中等规模优化问题上表现非常出色。这也是我选择MOPSO作为求解算法的核心原因。
2. 数学模型与约束条件详解
2.1 目标函数设计
先说两个目标函数，这是后面所有代码的核心，也是评判优化结果好坏的标尺。
第一个目标——系统网损最小。配电网中的有功网损可以表示为所有支路有功损耗之和：
[ f_1 = \sum_{i,j \in \Omega} G_{ij}(U_i^2 + U_j^2 - 2U_iU_j\cos\theta_{ij}) ]
其中 (U_i)、(U_j) 是节点电压幅值，(\theta_{ij}) 是节点 (i)、(j) 之间的相角差，(G_{ij}) 是支路电导。
这个公式看着复杂，但潮流计算软件或者Matpower一算就能出来，不用手算。实际代码中直接调用潮流计算函数，输入控制变量，输出各支路潮流，进而累加得到网损。
第二个目标——节点电压偏差最小。电压偏差用所有负荷节点电压偏移参考电压值的累计量来衡量：
[ f_2 = \sum_{k=1}^{N} \left| \frac{U_k - U_{ref}}{U_{ref}} \right| ]
(U_{ref}) 通常取1.0（标幺值），(N) 是节点数。这个指标反映了系统电压质量的优劣，数值越小代表电压越接近额定值。
2.2 等式约束与不等式约束
等式约束就是潮流平衡方程，这个不用多说，任何潮流计算都必须满足。在Matlab中一般通过调用Matpower或自己写的牛顿-拉夫逊法潮流程序来保证。
不等式约束分几类：
- 节点电压幅值约束：(U_{min} \leq U_i \leq U_{max})，一般是0.95~1.05 p.u.
- 支路潮流约束：(P_{ij} \leq P_{ij,max})
- 光伏逆变器无功出力约束：(-\sqrt{S_{PV}^2 - P_{PV}^2} \leq Q_{PV} \leq \sqrt{S_{PV}^2 - P_{PV}^2})，也就是逆变器剩余容量约束，这也是有功无功耦合的体现
- 有载调压变压器分接头约束：分接头档位在上下限之间，且只能整数调节
- 电容器组投切容量约束：投切容量是离散的，按组整数投切
2.3 光伏波动性的场景刻画
光伏波动性在数学上通常用场景集来描述。具体做法是：
第一步，根据历史光照数据拟合光伏出力的概率分布。常用的有Beta分布，因为Beta分布定义在[0,1]区间，非常适合描述光伏出力的标幺值。
第二步，通过蒙特卡洛采样或拉丁超立方采样，生成大量光伏出力场景。
第三步，用场景缩减技术（比如聚类方法，K-means或者同步回代消除法）将场景数量缩减到几个典型场景，同时保证典型场景能代表原始场景的概率特征。
每个典型场景都有对应的概率权重，优化目标可以转化为期望值形式：
[ \min \sum_{s=1}^{S} \pi_s \cdot f_i^{(s)} ]
其中 (\pi_s) 是场景 (s) 的概率，(f_i^{(s)}) 是该场景下的目标函数值。这种处理方式叫作随机优化，相比确定性优化更能真实反映光伏波动带来的影响。
本项目的参数选择上，我一般设置3~5个典型场景：一个高光照平稳场景、一个中等光照波动场景、一个低光照场景，再配一个极端波动场景（如云层快速移动导致出力骤升骤降）。极端场景虽然发生概率小，但对电压影响大，不纳入优化容易导致实际运行中的电压越限事故。
3. 多目标粒子群优化算法的实现细节
3.1 标准PSO回顾
标准粒子群算法源于鸟群觅食行为的模拟。每个粒子代表解空间中的一个候选解，粒子通过追踪个体历史最优（pbest）和群体历史最优（gbest）来更新自己的速度和位置：
[ v_{i}^{(t+1)} = w v_{i}^{(t)} + c_1 r_1 (pbest_i - x_i^{(t)}) + c_2 r_2 (gbest - x_i^{(t)}) ]
[ x_{i}^{(t+1)} = x_{i}^{(t)} + v_{i}^{(t+1)} ]
(w) 是惯性权重，(c_1)、(c_2) 是学习因子，(r_1)、(r_2) 是[0,1]之间的随机数。
惯性权重的设置很大程度上决定算法的搜索性能。我一开始直接用固定值0.8，发现前期收敛还行，后期容易陷入局部最优。后来改为线性递减策略：迭代初期权重0.9，强调全局搜索；迭代后期权重0.4，强调局部精细搜索。实测下来效果提升明显，Pareto前沿的分布均匀性也更好。
3.2 MOPSO的核心机制：Pareto支配、外部档案与拥挤距离
MOPSO在标准PSO的基础上增加了几项关键机制：
Pareto支配关系。对于最小化问题，解A支配解B，当且仅当A在所有目标上都不劣于B，且至少在一个目标上严格优于B。如果两个解互不支配，它们就是非支配关系。所有当前种群中的非支配解组成了Pareto前沿的近似。
外部档案集。MOPSO维护一个外部档案，用来存放搜索过程中找到的非支配解。每轮迭代后，将当前种群中的非支配解放入档案，并删除其中被新解支配的旧解。档案有容量上限，超出容量时需要用拥挤距离来裁剪。
拥挤距离。这个机制保证Pareto前沿的多样性。简单来说，计算档案中每个非支配解在目标空间中周围解的密度，密度越小（拥挤距离越大）的解越值得保留。在更新gbest时，我通常从档案集中拥挤距离较大的解里随机选一个，这样能让搜索向稀疏区域拓展，避免前沿集中在某个局部。
另外，还有一个需要注意的点是网格划分策略。MOPSO算法中还有一种用自适应网格来保持多样性的做法，把目标空间划分为若干网格，统计每个网格内解的数量，解越少的网格优先级越高。两种策略各有优劣，拥挤距离简单直接，网格划分计算量稍大但在前沿形状复杂时表现更好。我自己的实现里用的是拥挤距离法，原因是配电网问题的目标空间维度通常不高（就网损和电压偏差两个），拥挤距离的计算开销完全可控。
3.3 离散变量的处理方法
配电网优化里的控制变量并不全是连续的。有载调压变压器分接头是整数变量，电容器投切容量也是离散的。
处理离散变量的常见做法有三种：取整法、罚函数法、混合编码法。
取整法最简单，粒子飞行后在每次计算适应度之前对离散变量做四舍五入，但在某些情况下容易导致搜索方向失真，因为梯度信息被打乱了。
罚函数法是在目标函数中加入对离散变量越界的惩罚项，但需要反复调惩罚系数，比较麻烦。
混合编码法是我更推荐的做法，在速度更新时连续飞行，在位置更新时对指定维度做离散映射。具体来说，电容器投切组数这个维度，粒子位置记为连续值，但实际物理量是 (round(x)) 乘以单组容量。这样做既保持了PSO的搜索连续性，又保证了实际控制变量的物理可行性。
3.4 约束处理策略
配电网问题约束很多，处理不好容易出现不可行解。我采用的方法是外点罚函数法加可行性修复相结合。
对于电压越限、支路过载这类不等式约束，在目标函数后面加上惩罚项：
[ F_1 = f_1 + \lambda_1 \sum_{k} [\max(0, U_k - U_{max})]^2 + \lambda_2 \sum_{k} [\max(0, U_{min} - U_k)]^2 ]
简单解释就是，如果某个节点电压超出了上下限，就统计所有越限节点的越限量的平方和，乘上一个很大的惩罚系数，加到目标函数里。这样算法自然会避开那些电压越限的解。
对于潮流不收敛的情况（有时候粒子的位置对应着从数学上无法收敛的潮流状态），直接给这个粒子赋一个极大的适应度值，相当于宣判死刑，让它在后续迭代中自然被淘汰。
4. Matlab代码实现框架与步骤
4.1 整体代码结构
我实现的Matlab代码，整体模块划分如下：
main.m：主脚本，定义系统参数、算法参数，调用各模块执行优化，并输出结果initialization.m：粒子群初始化opf_objective.m：目标函数计算（网损+电压偏差）power_flow.m：潮流计算模块mo pso_update.m：粒子群速度和位置更新non_dominated_sort.m：非支配排序crowding_distance.m：拥挤距离计算update_archive.m：外部档案更新scene_generate.m：光伏出力典型场景生成
其中power_flow.m是计算瓶颈，每次迭代每个粒子都要调用一次潮流计算，所以这里的效率直接决定整个程序跑得快不快。
4.2 光伏场景生成的Matlab实现
场景生成的代码逻辑分为两步：采样和聚类。
先用Beta分布拟合历史光伏出力数据，通过betafit函数估计出形状参数 (a) 和 (b)，然后用betarnd(a,b,N,M)生成 (N) 个采样场景、每个场景 (M) 个时段的光伏出力序列。最后用K-means聚类，把N个场景缩减到S个典型场景。
% 场景生成示例代码 rng(42); % 固定随机种子，保证结果可复现 N_samples = 500; % 初始采样场景数量 S_scenes = 4; % 缩减后的典型场景数量 % 假设历史数据拟合得到Beta分布参数 alpha = 2.1; beta = 1.8; % 生成光伏出力标幺值样本 pv_samples = betarnd(alpha, beta, N_samples, 24); % K-means聚类缩减场景 [idx, centers] = kmeans(pv_samples, S_scenes, 'MaxIter', 500); % 计算每个典型场景的概率 weights = histcounts(idx, S_scenes) / N_samples; % 保存结果 save('pv_scenes.mat', 'centers', 'weights');这里有个细节：光伏出力在夜间为0，白天的Beta分布参数应该分段拟合。我通常的做法是把24小时分成几个时段，比如白天时段和夜间时段分开建模，而不是全天用一个分布，否则容易出现晚间采样出很大的光伏出力这种不合实际的场景。
4.3 MOPSO核心代码框架
粒子编码设计上，一个粒子的位置向量由四部分组成：
[ x = [P_{PV,curtail}(1:n_{PV}), Q_{PV}(1:n_{PV}), T_{OLTC}(1:n_{OLTC}), Q_{CB}(1:n_{CB})] ]
第一部分是光伏有功削减量，第二部分是光伏逆变器无功出力，第三部分是变压器分接头位置，第四部分是电容器组投切容量。因为分接头和电容器都是离散的，编码时单独区分处理。
粒子更新的核心代码如下：
% 粒子群速度和位置更新（含离散变量处理） for i = 1:nPop % 速度更新 inertia = w * pop(i).Velocity; cognitive = c1 * rand(size(pop(i).Position)) .* (pop(i).Best.Position - pop(i).Position); social = c2 * rand(size(pop(i).Position)) .* (global_best.Position - pop(i).Position); pop(i).Velocity = inertia + cognitive + social; % 位置更新 pop(i).Position = pop(i).Position + pop(i).Velocity; % 处理离散变量：变压器分接头档位（第nPV*2+1到nPV*2+nOLTC维） discrete_idx = nPV*2+1 : nPV*2+nOLTC+nCB; pop(i).Position(discrete_idx) = round(pop(i).Position(discrete_idx)); % 越界处理 pop(i).Position = min(max(pop(i).Position, lb), ub); end有几个细节需要特别注意：
惯性权重不是固定的，我采用线性递减的策略：
w = w_max - (w_max - w_min) * iter / maxIter;这样一来，迭代刚开始时粒子有较强的探索能力，可以在解空间大范围搜索；到后期权重变小，粒子会倾向于在当前最优解附近精细搜索，收敛性明显改善。
学习因子 (c_1) 和 (c_2) 我建议都取2.0。如果 (c_1) 过大，粒子容易过度偏向个体历史最优，导致探索不足；(c_2) 过大则会过早聚集到当前全局最优附近，容易陷入局部Pareto前沿。这两个值在绝大多数问题上取1.5~2.5之间都能有不错的效果，不需要过度调参。
4.4 外部档案更新与gbest选择
外部档案更新是本算法的核心。每迭代一步，把当前所有非支配解并入档案，然后去掉被支配的解，如果档案容量超出预设值，则用拥挤距离最小的几个解淘汰掉。
% 外部档案更新 function archive = update_archive(archive, pop) % 合并当前种群和档案中的解 all_solutions = [archive.solutions, pop]; % 找出所有非支配解 non_dominated = find_non_dominated(all_solutions); % 如果数量超过档案容量，按拥挤距离裁剪 if length(non_dominated) > archive_size cd = crowding_distance(non_dominated); [~, idx_sorted] = sort(cd, 'ascend'); non_dominated = non_dominated(idx_sorted(end-archive_size+1:end)); end archive.solutions = non_dominated; endgbest的选取直接关系到算法的探索方向。我选择从档案中拥挤距离较大的解里随机挑一个作为gbest，因为拥挤距离大的位置代表目标空间中的稀疏区域，向这里搜索可以让Pareto前沿更均匀地铺开，而不是堆在某一小段。
初始粒子群的产生也很关键。不能纯随机初始化，最好把一批粒子按零无功补偿、全部无功补偿、分接头中位等几个典型的工程场景摆上去，这样初始种群中就有一些性能尚可的可行解，算法收敛速度会快很多。
5. 实操中的问题与排查经验
5.1 潮流计算不收敛
这是我最常遇到的问题。粒子在搜索空间中飞到某些位置时，对应的控制参数组合可能导致潮流计算不收敛，尤其当光伏无功出力与电压支撑方向相反，或者变压器分接头挡位极端不合理的时候。
解决办法有两个方向。一个是尽量给初始种群中多放一些通过潮流验证的可行解，让粒子大部分时间在“安全区域”附近搜索。另一个是增加罚函数值，对潮流不收敛的粒子直接罚到无穷大，算法会逐步把搜索方向拉回可行域，但代价是收敛速度变慢。
实测下来，把潮流计算模块加上“最大迭代次数保护”很关键，如果牛顿-拉夫逊迭代超过一定次数就判定不收敛退出，而不是让它陷入死循环。这个保护可以大幅提升整体运行效率。
5.2 Pareto前沿分布不均匀
有时候跑出来的Pareto前沿集中在一小段，网损很小的解和电压偏差很小的解之间有空档。这种情况下，调度人员就很难选择折中方案。
排查后发现，问题往往出在拥挤距离更新机制上。最初我是每个迭代周期末尾统一更新一次拥挤距离，但种群在迭代过程中不断变化，这个值应该动态调整。改成每轮更新档案后立刻重算拥挤距离，并在选择gbest时基于最新的拥挤距离，前沿分布问题就解决了。
另外，增大外部档案容量也能改善分布均匀度，但代价是每一轮迭代耗时增加，需要权衡。一般取100左右比较合适。
5.3 离散变量导致搜索效率低
变压器分接头和电容器组这种离散变量的处理，一开始我直接对位置取整，结果发现算法经常在两个整数之间来回震荡，收敛速度很慢。
后来我换了一种思路：位置更新时不取整，只在计算目标函数时取整。这样粒子的“数学位置”是连续的，但“物理位置”是离散的。这样既不会失去PSO的连续搜索特性，又保证了实际可行性。注意最终的输出结果需要把离散变量重新取整还原成实际可操作的档位值。
这种方法在文献中经常被称为“连续松弛+离散映射”，在处理混合整数问题时相当好使。
5.4 光伏场景数量对结果的影响
光伏场景设置太多，计算量成倍增加；设置太少，又无法充分反映光伏波动性对结果的影响。我做了几组对比实验。
场景数为3时，计算速度快，但优化结果在某些极端天气下不具备足够的电压调节能力；场景数为7时，结果更鲁棒，但运行时间增加了近一倍。最后我选了5个场景作为折中：两个正常天气场景、一个低出力场景、一个高波动场景、一个极端场景。
实际做的时候可以通过肘部法则确定场景数，核心是看增加场景数带来的目标函数变化量是否明显减少。如果增加场景数对结果影响不大了，说明场景数已经足够。
5.5 运行时间优化
这个项目涉及的粒子数量多、迭代次数多，每个粒子每轮迭代都要跑一次潮流，确实很耗时。我做了几个优化：
一是尽量减小潮流计算中的冗余计算。如果某个控制变量变化前后，网络拓扑没有改变，可以利用上一次计算的电压初值作为本次迭代的初值，牛顿-拉夫逊法会更快收敛。
二是利用Matlab的向量化运算和并行计算。粒子群中每个粒子的潮流计算相互独立，完全可以用parfor替换for，我的四核机器实测大约能提速3倍左右。
三是给粒子群设置最大无改进迭代次数，如果连续多代Pareto前沿没有改进，就提前终止迭代，省出一大块时间。
6. 仿真结果解读与工程应用思考
6.1 典型案例结果分析
我在IEEE 33节点配网系统上做了仿真测试，把光伏接入节点设置在负荷较重、位置偏末端的地方，设定了5个光伏出力场景，用MOPSO进行求解。
优化后的Pareto前沿呈现出明显的下凸形态：左侧是电压偏差很小的解，右侧是网损很小的解，中间区域存在一系列折中方案。相比优化前，网损可以降低约18%到25%，电压偏差指标也能显著改善，具体数值取决于选中的Pareto解在曲线上的位置。
这个结果说明，在不增加额外硬件投资的前提下，通过合理协调光伏逆变器无功出力、变压器分接头档位和电容器组投切，主动配电网的运行效率和安全裕度都可以得到可观的提升。
6.2 光伏波动性对调度策略的影响
对比“计及光伏波动性”和“不考虑波动性”两种情况，可以发现重要差异：不考虑波动性的优化结果在典型场景下看似很好，网损很低、电压偏差小，但一旦光伏出力发生较大波动，电压越限和支路过载的风险显著上升；而计及波动性的方案虽然在常规场景下的网损略高一点，但在全部场景下都能保持电压在安全范围内。
这就好比一个人只按晴天规划出行，不准备雨具，遇到下雨就狼狈；而考虑了各种天气的人虽然背包重了一点，但在任何天气下都不慌。对于实际电网运行来说，安全裕度远比微小网损更有价值。
6.3 如何扩展和深入
这个项目的代码框架具备很好的扩展性。如果手头有实际配电网数据，只需要把IEEE 33节点的拓扑结构替换成实际网络，把负荷曲线和光伏历史数据换成实际采集数据，再重新设定节点参数，剩下的算法框架基本可以复用。
在这基础上，还可以做几个方向的延展。碳排放约束可以考虑进模型，变成三目标优化问题；储能系统的充放电策略也可以纳入控制变量；更进一步，可以用深度强化学习来替代元启发式算法，实现在线快速决策。
我个人的体会是，这类优化项目最大的价值不只是跑出一组漂亮的Pareto曲线，而是深刻理解模型假设和实际工程之间的差距。比如光伏波动性到底用什么分布来建模、哪些场景必须纳入考虑、离散变量如何处理才能既符合物理又方便求解——这些经验往往是论文和代码注释里不会写清楚的。把这些问题想透，改造成实际问题也就轻车熟路了。
