基于多目标粒子群优化算法的计及光伏波动性的主动配电网有功无功协调优化Matlab代码实现搞配电网优化这块的同行对“主动配电网”这个词应该都不陌生。传统配电网是被动接收电能但随着分布式光伏大规模接入电压越限、潮流倒送、网损攀升这些问题变得很棘手配电网必须从“被动”转向“主动”也就是要主动去管理配电网内部的分布式电源、储能、无功补偿装置等资源。而我今天要分享的这个项目就是围绕主动配电网调度中的一个核心难题——有功无功协调优化在计及光伏波动性这一实际约束下采用多目标粒子群优化算法MOPSO进行求解并给出完整的Matlab实现方案。这段时间我在做相关课题发现很多人一上来就套标准PSO忽略了一个关键点光伏出力是波动的不是恒定值。如果不把这个波动性建模进去你算出来的调度策略可能在实际运行中根本不可行。而这个项目恰恰把光伏波动性作为约束条件之一纳入优化模型同时用多目标框架处理网损和电压偏差这两个相互冲突的指标思路很清晰也非常贴合工程实际。这篇文章我就把整个项目的设计思路、数学模型、算法原理、Matlab实现框架以及我在实际调试中踩过的坑都整理出来希望能帮到正在做配电网优化调度或者打算用多目标粒子群算法做研究的朋友。1. 项目整体思路与研究价值1.1 为什么一定要考虑光伏的波动性光伏出力受光照强度、云层遮挡、温度等因素影响具有明显的随机性和间歇性。很多研究者在做优化时直接把某个时段的光伏出力当作固定值来处理这种做法在理论上是简化了但工程上往往站不住脚。打个比方就好比你规划一条公交线路假设每个站点上下车人数是固定的那班次安排特别容易但真实情况是早晚高峰和平峰的人流量差异很大你必须把这个波动考虑进去否则要么高峰期运力不足、要么平峰期空跑浪费。光伏波动性的引入也是这个道理它让配电网调度从“理想化的确定性优化”走向“贴近实际的鲁棒性优化”。在本项目中光伏波动性主要通过两方面体现在约束条件中引入节点电压的随机波动区间光伏出力的变化会改变潮流分布导致节点电压随之波动若优化结果没有预留足够的调节裕度实际运行时极易出现电压越限。在目标函数中体现为网损和电压偏差这两个指标的累积统计特性不是在单一场景下优化而是在典型光伏波动场景下综合优化。项目中采用的常见做法是设定多个典型光伏出力场景比如晴天平稳出力、多云波动出力、阴天低出力等然后让优化算法在所有这些场景下都能给出可行的调度结果这就叫“场景鲁棒”。1.2 有功无功协调优化的本质传统配电网调度中有功和无功往往分开处理有功靠发电机调整、无功靠电容器投切和有载调压变压器分接头调整。但在有高比例分布式光伏接入的主动配电网中有功潮流和无功电压是强耦合的。光伏逆变器既可以发有功也可以吸收或发出无功所以有功和无功必须联合优化才能达到更好的效果。举个例子某个节点光伏大发时如果只看有功调度只要保证光伏出力被消纳就行但此时线路可能因为潮流倒送导致末端节点电压升高。如果这时候合理调节逆变器的无功输出吸收一部分无功就能把电压压住同时还可能降低线路中的无功流动减少网损。这就是有功无功协调优化的价值所在。本项目建立的优化模型中控制变量一般包括分布式光伏的有功削减量如果允许光伏逆变器的无功输出量无功补偿装置如电容器组的投入容量有载调压变压器OLTC的分接头位置目标函数则通常包括两个系统网损最小节点电压偏差最小。这两个目标在数学上往往是矛盾的降低网损可能要求增加无功补偿这会影响电压分布而追求电压质量可能要以一定的网损为代价。所以这个问题天然适合用多目标优化算法求解得到一组Pareto最优解供调度人员根据实际需求灵活选择。1.3 为什么选多目标粒子群优化算法配电网有功无功协调优化是一个非线性的、多约束的、混合整数优化问题传统数学规划方法在处理这类问题时比较吃力尤其是控制变量中含有离散变量分接头档位、电容器组投切组数时求解难度更大。粒子群优化算法PSO的优势在于实现简单、不依赖目标函数的梯度信息、全局搜索能力强而且对连续和离散混合变量有很好的适应性。多目标粒子群优化算法MOPSO则是在PSO基础上引入Pareto支配关系、外部档案集和拥挤距离机制使得算法能够同时寻优多个相互冲突的目标最终得到一组均匀分布的Pareto前沿解。相比遗传算法NSGA-IIMOPSO收敛速度更快参数更少在配电网这类中等规模优化问题上表现非常出色。这也是我选择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_{k1}^{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_{s1}^{S} \pi_s \cdot f_i^{(s)} ]其中 (\pi_s) 是场景 (s) 的概率(f_i^{(s)}) 是该场景下的目标函数值。这种处理方式叫作随机优化相比确定性优化更能真实反映光伏波动带来的影响。本项目的参数选择上我一般设置3~5个典型场景一个高光照平稳场景、一个中等光照波动场景、一个低光照场景再配一个极端波动场景如云层快速移动导致出力骤升骤降。极端场景虽然发生概率小但对电压影响大不纳入优化容易导致实际运行中的电压越限事故。3. 多目标粒子群优化算法的实现细节3.1 标准PSO回顾标准粒子群算法源于鸟群觅食行为的模拟。每个粒子代表解空间中的一个候选解粒子通过追踪个体历史最优pbest和群体历史最优gbest来更新自己的速度和位置[ v_{i}^{(t1)} w v_{i}^{(t)} c_1 r_1 (pbest_i - x_i^{(t)}) c_2 r_2 (gbest - x_i^{(t)}) ][ x_{i}^{(t1)} x_{i}^{(t)} v_{i}^{(t1)} ](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*21到nPV*2nOLTC维 discrete_idx nPV*21 : nPV*2nOLTCnCB; 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_size1: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曲线而是深刻理解模型假设和实际工程之间的差距。比如光伏波动性到底用什么分布来建模、哪些场景必须纳入考虑、离散变量如何处理才能既符合物理又方便求解——这些经验往往是论文和代码注释里不会写清楚的。把这些问题想透改造成实际问题也就轻车熟路了。