做主动配电网优化这个方向的人十有八九会在“怎么把光伏波动性放进去”这件事上卡一阵子。光伏不像火电那样出力稳定它像一路随光照强度上下跳的电源中午可能猛烈反送傍晚又断崖式掉功率你要是拿一组固定出力去做优化投到实际里基本没法用。以前我刚开始做这块时习惯性用线性加权把网损和电压偏差变成一个单目标再丢给标准PSO去跑结果跑出来的解很“拧巴”——低压台区电压合格了网损却高得离谱网损压下来了末端电压又越限。权重的调整特别主观20次实验换6组权重能给你6种完全不同风格的方案。后来换成多目标粒子群优化算法MOPSO把网损最小和电压偏差最小当作两个并列的目标用Pareto支配关系去寻优问题的思路一下就顺了。这篇笔记就是围绕这个完整项目展开的工程总结从光伏波动性建模、MOPSO算法设计一直写到Matlab代码怎么组织、常见坑怎么排查。适合正在做主动配电网方向研究的研究生、做分布式光伏接入规划的工程师也适合刚接触多目标优化想快速上手的读者。1. 项目背景与问题拆解光伏波动性为什么必须“计及”1.1 高渗透率光伏带来的“双高”挑战传统配电网的设计前提是单向潮流变电站流向负荷电压沿馈线逐渐降低末端低电压是主要矛盾。光伏大量接入后这个前提被打破了。馈线上出现了分布式电源潮流可能反向流动末端电压在中午光伏大发时反而会抬升。我见过一个10kV馈线的案例轻载时段接入2MW光伏后末端电压从0.98标幺值直接抬到1.04以上如果再叠加次日多云天气下的功率波动电压会在短时间内反复越限这已经不是靠传统人工调压能解决的问题了。这类问题通常被概括为“双高”——高渗透率、高波动性。高渗透率意味着光伏在馈线功率中的占比大影响不再是“局部小扰动”高波动性意味着场景切换快中午晴天和云遮时段的光伏出力可能相差70%以上。如果你在优化模型里只用一组确定性的光伏出力数据那么最优解只对那一个时刻有效换一个光照条件可能完全失效。所以“计及光伏波动性”不是锦上添花而是主动配电网优化必须考虑的基础条件。这里有个容易混淆的点光伏波动性到底在优化模型里怎么体现一种是把它放进场景集里用多个典型场景描述可能的光照情况然后优化所有场景下的期望目标另一种是把它写进鲁棒约束比如要求电压越限概率小于5%。后者更严格但计算量大做研究时通常用前者也就是场景期望模型。这也是本项目的核心思路。1.2 有功和无功协调优化的内在逻辑配电网里的电压调节设备并不少有载调压变压器OLTC、并联电容器组、静止无功补偿器SVC、分布式光伏逆变器的无功能力。这些都是无功手段。无功补偿解决低电压和高电压问题的“性价比”很高因为无功功率就地平衡可以显著改善电压分布。但无功手段有一个天然局限它改变的是电压幅值和无功潮流对有功功率在线路上的流动影响有限。馈线重载时有功潮流带来的网损大头仍然存在光靠投电容器、调变比是压不下来的。有功手段则不同。光伏限功率、储能充放电直接改变馈线上的有功流动方向。比如中午光伏倒送导致末端电压偏高适当削减光伏出力或者让储能充电潮流方向改变电压很快就降下来。但有功调节的代价也很直观光伏限功率等于放弃发电量储能充放电则涉及电池寿命和调度成本。所以这里有了一组矛盾无功调节成本低但能力有限有功调节效果好但有代价两者必须一起考虑。把它比喻成开车下山可能不太贴切我更习惯把它想成调节供水管网无功手段像调水压阀有功手段像调进水总量。坡度大了光靠调压阀会憋管光调总量又可能水压不足。主动配电网的有功无功协调优化就是要把这两类手段放进同一个决策框架让它们互相配合。这也是为什么题目特别强调“协调”——如果分两步做先优化变比和无功再优化有功那么第二步的结果往往会推翻第一步的最优方案因为电压和潮流是耦合的。1.3 多目标优化模型的具体数学形式既然要协调就要有一个明确的数学模型。本项目把优化目标定义为两个网络损耗最小、电压偏差最小。用公式写出来是f1 Σ_s ω_s · ( Σ_t P_loss(s,t) ) f2 Σ_s ω_s · ( Σ_t Σ_i (V_i(s,t) - V_ref)^2 )其中 s 表示光伏场景ω_s 是场景概率t 是时段i 是节点P_loss 是网络有功损耗V_ref 一般取 1.0 标幺值。两个目标函数都按所有场景的期望值来算这样就把“光伏波动性”纳入了优化目标——一组控制变量好不好要看它在所有典型场景下的平均表现。决策变量主要分两类一类是离散控制变量包括OLTC变比档位、电容器组投切组数一类是连续控制变量包括各分布式光伏逆变器的无功出力以及可选的光伏有功削减率。有些研究还会加入储能充放电功率这时决策变量就更多了。下表列一个常见的决策变量集合决策变量类型取值范围示例说明OLTC变比档位离散整数-8 ~ 8 档每档0.0125调节变电站侧电压电容器投切组数离散整数0 ~ 10 组每组容量约50kvar光伏逆变器无功出力连续[-0.8S, 0.8S]功率因数范围内可调光伏有功削减率连续[0, 0.2]极端场景下限制出力储能充放电功率连续[-P_max, P_max]可选视系统配置约束条件包括节点功率平衡、节点电压上下限、支路电流容量、OLTC变比范围、电容器投切上下限、光伏逆变器视在功率约束等。其中功率平衡方程由潮流计算自然满足其余不等式约束要在优化过程中通过罚函数或可行性修复来处理。这一步做好了后面的算法才有意义。2. 光伏波动性的建模与场景削减2.1 光伏出力模型从光照强度到并网功率光伏出力波动性的根源是光照强度的随机变化。工程上常用Beta分布来描述一段时间内光照强度的概率分布因为光照强度是一个有上下界的随机变量Beta分布的形态灵活和实测光照数据拟合得也不错。Beta分布的概率密度函数是f(G) Γ(αβ) / (Γ(α)Γ(β)) · (G/G_max)^(α-1) · (1 - G/G_max)^(β-1)其中 G 是实际光照强度G_max 是最大光照强度α 和 β 是分布的形状参数。这两个参数可以从历史光照数据的平均值 μ 和标准差 σ 反推α μ²(1-μ)/σ² - μ β μ(1-μ)²/σ² - (1-μ)实际项目里我不会直接用光照强度作为决策输入而是把它转成光伏并网功率。典型公式是P_pv η_pv · S_pv · G · (1 - β_T · (T_c - T_ref))η_pv 是光伏板转换效率S_pv 是光伏阵列面积β_T 是温度系数T_c 是电池板工作温度T_ref 取25℃。温度对出力的影响在夏季中午尤其明显板温可能到50℃以上效率掉得比较厉害。如果你手头没有温度数据可以作简化处理只保留光照强度的影响但论文或报告里最好说明这个简化假设。做仿真时光照强度的时序数据通常会按小时或15分钟一个时段来离散。每个时段按照Beta分布抽样就可以生成一整天的光伏出力曲线。这里要注意的是不同时段之间往往存在相关性严格的做法是用时序采样或马尔可夫链而不是每个时段独立抽样。不过对于以潮流优化为主的论文和初步工程分析独立时段抽样加场景削减已经足够误差在可接受范围内。2.2 蒙特卡洛场景生成与K-means削减有了光伏出力模型下一步是生成场景集。最直接的方法是蒙特卡洛抽样对每个时段的Beta分布抽取大量样本组合成N个完整的光伏出力场景。比如生成1000个场景每个场景包含24个时段的光照强度。但1000个场景意味着每次目标函数评估要跑1000次潮流计算粒子群算法跑50个粒子迭代100代那就是500万次潮流计算量完全不可接受。所以场景必须削减。削减的思路是从1000个原始场景中选出M个有代表性的场景M通常取515并给每个场景赋一个概率。常用算法有两种同步回代削减法backward reduction和K-means聚类法。同步回代的核心思想是贪心删除——每次找与其他场景距离之和最小的场景删掉把它的概率累加到最近邻场景上直到剩下M个。K-means更简单直接把所有原始场景聚成M类用聚类中心作为典型场景用类内的样本占比作为场景概率。Matlab里K-means实现非常直接% 生成 1000 个光伏出力场景每个场景 24 个时段 N 1000; M 10; S betarnd(alpha, beta, N, 24); % 按Beta分布抽样 % 归一化到光伏额定容量 S S * P_pv_rated; % K-means削减到10个典型场景 [cluster_index, centers] kmeans(S, M, Distance, sqeuclidean, Replicates, 20); prob accumarray(cluster_index, 1) / N;Replicates设为20是为了防止K-means陷入局部最优初始化跑20次取最优聚类结果。削减后centers就是典型场景的光伏出力曲线prob是对应的概率。这样计算量就从1000次潮流降到了10次目标函数评估速度快了100倍。这里有一个经验场景数不是越多越好。10个场景和15个场景的优化结果差异通常很小但计算时间差了50%。如果做灵敏度分析或参数扫描我一般把场景数压到812个既能保留光伏出力的分布特性又不至于让仿真等得心焦。2.3 波动性进入目标函数的方式期望目标与越限惩罚场景生成好之后难点是怎么在目标函数里把波动性体现出来。最简单的方式是用所有场景的期望指标。比如网损期望值f1 Σ_s prob(s) · P_loss(s)也就是说对一组给定的控制变量先逐个场景调用潮流计算算出每个场景下的网损和电压偏差再按场景概率加权求和。这组控制变量在所有典型场景下的平均表现就得到了光伏波动性自然被“计及”了进去。但只算期望值有一个问题期望值会把极端场景“平均掉”。比如10个场景里9个电压正常1个场景末端电压到了1.07这个越限情况对期望电压偏差的影响可能并不大实际运行却可能触发保护跳闸。因此我会在目标函数或约束条件里加一层“越限惩罚”统计所有场景中的电压越限情况并乘一个大系数加入目标。做法也很直接viol_prob 0; for s 1:M V runPowerFlow(x, scenario(s)); viol_s max(0, V - V_max) max(0, V_min - V); if any(viol_s 0) viol_prob viol_prob prob(s); end end f1 f1 lambda_penalty * viol_prob;这样算法在寻优时就会主动避开那些会让电压在某些场景下越限的方案优化结果更贴近实际运行要求。如果你做的是机会约束规划可以把viol_prob大于给定阈值比如5%的判断直接写成约束但罚函数法实现起来更简单调参也直观。3. 多目标粒子群优化算法的核心设计3.1 为什么单目标加权法在这里不够用很多第一次接触多目标优化的人会问为什么不把网损和电压偏差乘上权重相加变成一个目标函数我之前也这么干过但实际做下来至少有三个问题。第一个是权重的主观性。网损单位是千瓦电压偏差单位是标幺值两者的数量级差很多不归一化的话权重几乎没有物理意义即便归一化权重怎么取还是靠试。第二个更致命——加权法遇到非凸Pareto前沿时有些最优解无论怎么取权重都找不到。第三配电网里OLTC和电容器是离散变量目标函数本身不光滑加权法很容易卡在局部解附近。所以这个项目用多目标粒子群优化算法MOPSO不是炫技而是因为它能一次性给出整条Pareto前沿。Pareto前沿上的每个点都代表一个“不能在不牺牲一个目标的前提下改善另一个目标”的折中方案。做研究也好做工程方案也罢拿到的是一组候选方案而不是被一个权重锁死的唯一解后续决策空间大很多。下面简单对比一下两种思路对比项单目标加权法MOPSO结果单个解一组Pareto解集权重需人工设定敏感不需要解集中自动覆盖折中非凸前沿可能丢失解可以覆盖计算量较小较大但可控结果决策无选择空间可结合偏好选折中解3.2 Pareto支配与外部档案机制再解释一下Pareto支配用大白话说就是方案A在所有目标上都不比方案B差并且至少有一个目标严格比B好那A就支配B。所有不被任何其他方案支配的解构成Pareto前沿。MOPSO里有一个专门存放非支配解的地方叫外部档案Archive。每次迭代结束把当前种群中的非支配解塞进外部档案再剔除掉已经被其他解支配的旧解。但外部档案容量有限默认100个装满了怎么办就需要裁剪。裁剪的原则是保留“稀疏区域”的解删掉“拥挤区域”的解这样Pareto前沿才能分布均匀而不是堆在一小块区域。计算拥挤度的方法有很多我常用网格法把目标空间划分成网格统计每个网格里的粒子数粒子越多的网格表示这一带的解越拥挤优先从这些网格里删。这相当于在目标空间里做了一次密度估计实现简单效果也不错。3.3 速度更新、离散变量处理与变异机制MOPSO的粒子更新公式和标准PSO一样v(t1) w·v(t) c1·r1·(pbest - x) c2·r2·(gbest - x) x(t1) x(t) v(t1)w是惯性权重c1和c2是学习因子r1和r2是[0,1]均匀随机数。区别在于gbest怎么选单目标PSO里gbest是全局最优解MOPSO里没有单一的最优解所以gbest从外部档案里选而且倾向于从粒子密度最低的网格里选目的是引导粒子去探索Pareto前沿上还比较“空旷”的区域。配电网优化里有离散变量OLTC档位和电容器投切组数本身是整数。粒子位置更新后是连续值必须做取整处理。但直接四舍五入会导致粒子在某个档位附近反复跳动影响收敛因此我通常还加一个边界检查和邻域扰动。如果粒子越界把它拉回边界如果连续多次收敛到同一个离散值附近就随机探索相邻档位避免死锁。变异机制也是MOPSO保持多样性的关键。标准PSO没有变异MOPSO在早期容易聚到局部Pareto前沿。我习惯让变异率随时间递减p_m(t) 0.2 - 0.15 * (t / T_max)迭代初期变异率高粒子探索范围大后期变异率降低收敛到局部精细搜索。这个机制对配电网这种多局部极值的问题非常有效。核心更新代码骨架如下for i 1:pop 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, discreteIdx) round(x(i, discreteIdx)); % 边界检查 x(i,:) max(min(x(i,:), ub), lb); % 变异 if rand p_m idx randi(nVar); x(i, idx) lb(idx) rand * (ub(idx) - lb(idx)); if ismember(idx, discreteIdx) x(i, idx) round(x(i, idx)); end end end这段代码是整个MOPSO的核心后面接上潮流计算和目标函数评估就是一个可以跑通的闭环。4. Matlab代码实现从主程序到潮流计算4.1 项目目录与主程序流程Matlab代码的组织直接影响调试效率。我现在的习惯是把功能拆开每个文件只做一件事目录结构大致是这样moPso_ADN/ ├─ main.m % 主程序设置参数并调用优化 ├─ case33.m % 配电网基础数据33节点系统为例 ├─ loadData.m % 负荷数据、光伏场景数据 ├─ genScenarios.m % 生成光伏场景并削减 ├─ moPSO_Core.m % MOPSO主算法 ├─ powerflow.m % 前推回代潮流计算 ├─ objectiveFun.m % 目标函数评估 ├─ plotResults.m % 绘制Pareto前沿和节点电压分布main.m只做五件事加载数据、生成场景、初始化MOPSO参数、调用优化循环、输出和画图。结构清晰之后调参和加功能都很方便。主程序的核心流程大致是% main.m clc; clear; close all; load(case33.mat); % 网络阻抗、拓扑 [Scenario, Prob] genScenarios(1000, 10); % 生成10个典型场景 % M0PSO 参数 pop 50; maxGen 100; archiveSize 100; nGrid 20; [front, bestSolution] moPSO_Core(...); plotResults(front, bestSolution);4.2 前推回代潮流计算配电网版本的“迭代收敛”配电网优化里潮流计算是内层循环每个粒子每次迭代都要调用很多次因此它的速度决定了整个算法能不能跑完。33节点这类辐射状配电网我强烈推荐前推回代法而不是牛顿拉夫逊法。原因有两点前推回代不需要求导和形成雅可比矩阵每步只是简单的复数运算速度快辐射状网结构天然适合前推回代收敛性也很稳定。前推回代的基本思路分两步回代从末端节点向首端节点根据当前节点电压和负荷/电源注入功率计算每条支路的电流或功率。前推从首端节点向末端节点根据支路电流和线路阻抗更新各节点电压。两步骤反复迭代直到前后两次节点电压幅值差的最大值小于收敛阈值比如1e-6。Matlab代码骨架function [V, P_loss] powerflow(branch, bus, load_power, gen_power) n size(bus, 1); V ones(n, 1); I_branch zeros(size(branch, 1), 1); for iter 1:100 V_old V; % 回代从末端往首端算支路电流 for k size(branch,1):-1:1 from branch(k,1); to branch(k,2); S_inj load_power(to) - gen_power(to); I_branch(k) conj(S_inj / V(to)); % 累加支路下游电流略实际需拓扑排序 end % 前推从首端往末端算节点电压 for k 1:size(branch,1) from branch(k,1); to branch(k,2); V(to) V(from) - branch(k,3) * I_branch(k); end if max(abs(V - V_old)) 1e-6 break; end end P_loss sum(abs(I_branch).^2 .* branch(:,3)); end实际工程代码要考虑拓扑顺序先做广度优先搜索给支路编号然后按逆拓扑序回代、按拓扑序前推。前推回代对纯辐射状网收敛很快一般迭代几十次就达到1e-6。如果网络里有弱环需要先开环再补偿会稍微复杂一点但主流配电网优化算例基本都用辐射状网。光伏逆变器有剩余无功能力时我会让光伏节点在潮流中按“恒无功注入”处理无功值由决策变量给出如果某些场景下节点电压越限再把无功修正为维持电压的PV控制模式。两块逻辑加在一起仿真结果会更贴近真实逆变器的控制策略。4.3 目标函数与罚函数让搜索“听话”objectiveFun.m是核心评估函数输入是一个粒子的决策变量输出是两个目标函数值和违约量。大致的流程是function [f1, f2] objectiveFun(x, Scenario, Prob, data) % 解析决策变量 tap x(1); cap x(2); dg_q x(3:end); % 初始化目标 f1 0; f2 0; for s 1:length(Prob) % 根据场景修改光伏出力 data.gen(:,2) Scenario(s,:); % 潮流计算 [V, P_loss] powerflow(data, tap, cap, dg_q); % 电压偏差 v_dev sum((V - 1).^2); % 场景加权累加 f1 f1 Prob(s) * P_loss; f2 f2 Prob(s) * v_dev; end % 罚函数处理电压越限 f1 f1 1000 * max(0, max(abs(V-1) - 0.05)); f2 f2 1000 * max(0, max(abs(V-1) - 0.05));罚函数系数取1000是我常用的值。太小了惩罚不够越限方案可能留在Pareto前沿上太大了会压制其他目标的梯度信息导致粒子早早丢失多样性。1000在大多数量纲下都够用如果你发现优化结果中有大量越限解可以把它提到1e4。4.4 参数设置速查表与归一化细节MOPSO参数对收敛速度和Pareto前沿质量影响很大给一张我调通这类型项目的速查表参数推荐值说明种群规模 pop50 ~ 100决策变量多时取大最大迭代次数 maxGen100 ~ 300看收敛曲线再调整惯性权重 w0.9 → 0.4 线性递减前期全局探索后期局部搜索学习因子 c1、c22.0、2.0常见配置也可略做调整外部档案容量100太大拥挤度计算变慢网格划分数20 × 20按目标个数调整变异率初值/终值0.2 → 0.05线性递减这里特别要强调归一化。外部档案的拥挤度计算依赖目标值之间的距离网损可能是几千瓦电压偏差可能是零点几甚至几点几如果不归一化拥挤度全被网损主导电压偏差的多样性就丢了。我在moPSO_Core里对每个目标都做了一遍min-max归一化再算网格和拥挤度Pareto前沿立刻就均匀了。这一步非常关键算是踩了坑才得出的经验。5. 调试实录常见问题与排查技巧5.1 结果不收敛、Pareto前沿稀烂时先排查这三件事这个项目调试过程里最折腾的问题是算法跑完Pareto前沿稀稀疏疏甚至只有几个点。后来总结出三个排查顺序。第一先跑单场景再跑多场景。光伏场景一多每个粒子都要对所有场景做潮流计算量大不说出了问题根本不知道是潮流算错了还是场景概率算错了。我的一般做法是先把场景数临时改成1跑确定性优化确认潮流、罚函数、粒子更新这些环节在单场景下都能正常收敛再切回多场景模型。第二检查潮流迭代次数和残差。前推回代如果不收敛通常是因为设定的迭代上限太小或负荷/电源数据有误导致功率严重不平衡。在powerflow里加一行disp(iter, max(abs(V - V_old)))观察残差变化一秒就能定位。第三检查目标函数量纲和归一化。如果f1和f2的量级差一个数量级以上Pareto前沿就会偏向某一侧。先把两个目标都归一化到[0,1]再跑MOPSO前沿分布会正常很多。5.2 离散变量导致振荡换挡永远在犹豫OLTC变比和电容器投切组数是离散变量粒子速度更新后四舍五入取整这是最简单的处理方式但有个副作用粒子可能在两个相邻档位之间反复横跳尤其是迭代后期速度已经很小的时候四舍五入的取整操作会把粒子“踢”到相邻档位目标函数值随之跳变外部档案里的非支配解也跟着震荡。我的解决办法是分两步走。第一步在粒子更新流程里离散变量取整后加一个邻域概率扰动。代码里我常写if rand 0.1 x(i, discreteIdx) x(i, discreteIdx) randsample([-1 0 1], 1); x(i, discreteIdx) round(max(min(x(i, discreteIdx), ub(discreteIdx)), lb(discreteIdx))); end第二步外部档案里保存的解全部保证是合法取整后的结果绝不允许把未取整的连续值放进档案。否则最后输出的“最优解”可能对应不存在的OLTC档位工程上完全不可用。做完这两步离散变量振荡的问题基本就消失了。5.3 电压越限概率与折中解选取多目标优化得到Pareto前沿之后怎么从100个解里选一个如果做工程报告我喜欢用模糊隶属度法。先找出每个目标在所有非支配解中的最大值和最小值然后计算每个解的隶属度μ_i (f_i_max - f_i) / (f_i_max - f_i_min)把每个解的所有目标隶属度累加取总隶属度最大的那个作为折中解。这个方法的直观含义是挑一个在网损和电压偏差两个目标上相对最接近“各目标最好值”的方案。代码很简单fmax max(front); fmin min(front); mu (repmat(fmax, size(front,1), 1) - front) ./ (repmat(fmax - fmin, size(front,1), 1)); total_mu sum(mu, 2); [~, bestIdx] max(total_mu); bestSolution front(bestIdx, :);选完折中解后我还会在全部场景下重新运行一遍潮流统计每个节点的电压越限概率画成箱线图或者“电压区间包络线”。这个结果放到论文或报告里非常有说服力比单纯贴一张Pareto前沿散点图更能说明“计及光伏波动性”的效果。画图用Matlab的fill函数画包络带效果不错。% 统计各节点电压分布 V_all zeros(numStatic, numScenarios); for s 1:numScenarios V_all(:, s) runPowerFlow(bestSolution, Scenario(s)); end V_min min(V_all, [], 2); V_max max(V_all, [], 2); fill([1:n; n:-1:1], [V_min; flip(V_max)], [0.9 0.9 0.9]);5.4 场景削减与计算时间的平衡最后聊一个容易被忽略的细节场景削减比例和计算时间的平衡。如果你的算例是IEEE 33节点系统10个场景100个粒子跑100代单次潮流计算大概几毫秒总体时间还能接受。但如果换成大规模配电网比如300节点以上就必须考虑优化计算效率了。我有两个实操建议。第一把场景削减后的场景数控制在58个优先用K-means不要用同步回代——同步回代在样本量大时非常慢K-means在Matlab里有内置函数向量化运算快得多。第二潮流计算做矩阵化预处理把支路阻抗、拓扑结构提前展开成矩阵前推回代改成矩阵运算而不是for循环速度能提升好几倍。本质上做工程仿真精度和耗时的平衡一直是重点。最后说一点个人经验这类项目最容易出现的认知误区是把算法当主角把配电网模型当配角。事实上模型建模对了场景生成合理了哪怕用标准MOPSO也能得到不错的结果反过来模型和场景一团糟再高级的算法也救不回来。我自己的习惯是先把单场景确定性问题调通然后再引入波动性接着再上多目标算法每一步过了再往后走。这个顺序看起来“慢”但实际是整个项目里省时间最多的做法。如果你也在折腾类似的方向建议试试这个思路应该能少走不少弯路。
主动配电网优化：计及光伏波动性的MOPSO算法与Matlab实践
做主动配电网优化这个方向的人十有八九会在“怎么把光伏波动性放进去”这件事上卡一阵子。光伏不像火电那样出力稳定它像一路随光照强度上下跳的电源中午可能猛烈反送傍晚又断崖式掉功率你要是拿一组固定出力去做优化投到实际里基本没法用。以前我刚开始做这块时习惯性用线性加权把网损和电压偏差变成一个单目标再丢给标准PSO去跑结果跑出来的解很“拧巴”——低压台区电压合格了网损却高得离谱网损压下来了末端电压又越限。权重的调整特别主观20次实验换6组权重能给你6种完全不同风格的方案。后来换成多目标粒子群优化算法MOPSO把网损最小和电压偏差最小当作两个并列的目标用Pareto支配关系去寻优问题的思路一下就顺了。这篇笔记就是围绕这个完整项目展开的工程总结从光伏波动性建模、MOPSO算法设计一直写到Matlab代码怎么组织、常见坑怎么排查。适合正在做主动配电网方向研究的研究生、做分布式光伏接入规划的工程师也适合刚接触多目标优化想快速上手的读者。1. 项目背景与问题拆解光伏波动性为什么必须“计及”1.1 高渗透率光伏带来的“双高”挑战传统配电网的设计前提是单向潮流变电站流向负荷电压沿馈线逐渐降低末端低电压是主要矛盾。光伏大量接入后这个前提被打破了。馈线上出现了分布式电源潮流可能反向流动末端电压在中午光伏大发时反而会抬升。我见过一个10kV馈线的案例轻载时段接入2MW光伏后末端电压从0.98标幺值直接抬到1.04以上如果再叠加次日多云天气下的功率波动电压会在短时间内反复越限这已经不是靠传统人工调压能解决的问题了。这类问题通常被概括为“双高”——高渗透率、高波动性。高渗透率意味着光伏在馈线功率中的占比大影响不再是“局部小扰动”高波动性意味着场景切换快中午晴天和云遮时段的光伏出力可能相差70%以上。如果你在优化模型里只用一组确定性的光伏出力数据那么最优解只对那一个时刻有效换一个光照条件可能完全失效。所以“计及光伏波动性”不是锦上添花而是主动配电网优化必须考虑的基础条件。这里有个容易混淆的点光伏波动性到底在优化模型里怎么体现一种是把它放进场景集里用多个典型场景描述可能的光照情况然后优化所有场景下的期望目标另一种是把它写进鲁棒约束比如要求电压越限概率小于5%。后者更严格但计算量大做研究时通常用前者也就是场景期望模型。这也是本项目的核心思路。1.2 有功和无功协调优化的内在逻辑配电网里的电压调节设备并不少有载调压变压器OLTC、并联电容器组、静止无功补偿器SVC、分布式光伏逆变器的无功能力。这些都是无功手段。无功补偿解决低电压和高电压问题的“性价比”很高因为无功功率就地平衡可以显著改善电压分布。但无功手段有一个天然局限它改变的是电压幅值和无功潮流对有功功率在线路上的流动影响有限。馈线重载时有功潮流带来的网损大头仍然存在光靠投电容器、调变比是压不下来的。有功手段则不同。光伏限功率、储能充放电直接改变馈线上的有功流动方向。比如中午光伏倒送导致末端电压偏高适当削减光伏出力或者让储能充电潮流方向改变电压很快就降下来。但有功调节的代价也很直观光伏限功率等于放弃发电量储能充放电则涉及电池寿命和调度成本。所以这里有了一组矛盾无功调节成本低但能力有限有功调节效果好但有代价两者必须一起考虑。把它比喻成开车下山可能不太贴切我更习惯把它想成调节供水管网无功手段像调水压阀有功手段像调进水总量。坡度大了光靠调压阀会憋管光调总量又可能水压不足。主动配电网的有功无功协调优化就是要把这两类手段放进同一个决策框架让它们互相配合。这也是为什么题目特别强调“协调”——如果分两步做先优化变比和无功再优化有功那么第二步的结果往往会推翻第一步的最优方案因为电压和潮流是耦合的。1.3 多目标优化模型的具体数学形式既然要协调就要有一个明确的数学模型。本项目把优化目标定义为两个网络损耗最小、电压偏差最小。用公式写出来是f1 Σ_s ω_s · ( Σ_t P_loss(s,t) ) f2 Σ_s ω_s · ( Σ_t Σ_i (V_i(s,t) - V_ref)^2 )其中 s 表示光伏场景ω_s 是场景概率t 是时段i 是节点P_loss 是网络有功损耗V_ref 一般取 1.0 标幺值。两个目标函数都按所有场景的期望值来算这样就把“光伏波动性”纳入了优化目标——一组控制变量好不好要看它在所有典型场景下的平均表现。决策变量主要分两类一类是离散控制变量包括OLTC变比档位、电容器组投切组数一类是连续控制变量包括各分布式光伏逆变器的无功出力以及可选的光伏有功削减率。有些研究还会加入储能充放电功率这时决策变量就更多了。下表列一个常见的决策变量集合决策变量类型取值范围示例说明OLTC变比档位离散整数-8 ~ 8 档每档0.0125调节变电站侧电压电容器投切组数离散整数0 ~ 10 组每组容量约50kvar光伏逆变器无功出力连续[-0.8S, 0.8S]功率因数范围内可调光伏有功削减率连续[0, 0.2]极端场景下限制出力储能充放电功率连续[-P_max, P_max]可选视系统配置约束条件包括节点功率平衡、节点电压上下限、支路电流容量、OLTC变比范围、电容器投切上下限、光伏逆变器视在功率约束等。其中功率平衡方程由潮流计算自然满足其余不等式约束要在优化过程中通过罚函数或可行性修复来处理。这一步做好了后面的算法才有意义。2. 光伏波动性的建模与场景削减2.1 光伏出力模型从光照强度到并网功率光伏出力波动性的根源是光照强度的随机变化。工程上常用Beta分布来描述一段时间内光照强度的概率分布因为光照强度是一个有上下界的随机变量Beta分布的形态灵活和实测光照数据拟合得也不错。Beta分布的概率密度函数是f(G) Γ(αβ) / (Γ(α)Γ(β)) · (G/G_max)^(α-1) · (1 - G/G_max)^(β-1)其中 G 是实际光照强度G_max 是最大光照强度α 和 β 是分布的形状参数。这两个参数可以从历史光照数据的平均值 μ 和标准差 σ 反推α μ²(1-μ)/σ² - μ β μ(1-μ)²/σ² - (1-μ)实际项目里我不会直接用光照强度作为决策输入而是把它转成光伏并网功率。典型公式是P_pv η_pv · S_pv · G · (1 - β_T · (T_c - T_ref))η_pv 是光伏板转换效率S_pv 是光伏阵列面积β_T 是温度系数T_c 是电池板工作温度T_ref 取25℃。温度对出力的影响在夏季中午尤其明显板温可能到50℃以上效率掉得比较厉害。如果你手头没有温度数据可以作简化处理只保留光照强度的影响但论文或报告里最好说明这个简化假设。做仿真时光照强度的时序数据通常会按小时或15分钟一个时段来离散。每个时段按照Beta分布抽样就可以生成一整天的光伏出力曲线。这里要注意的是不同时段之间往往存在相关性严格的做法是用时序采样或马尔可夫链而不是每个时段独立抽样。不过对于以潮流优化为主的论文和初步工程分析独立时段抽样加场景削减已经足够误差在可接受范围内。2.2 蒙特卡洛场景生成与K-means削减有了光伏出力模型下一步是生成场景集。最直接的方法是蒙特卡洛抽样对每个时段的Beta分布抽取大量样本组合成N个完整的光伏出力场景。比如生成1000个场景每个场景包含24个时段的光照强度。但1000个场景意味着每次目标函数评估要跑1000次潮流计算粒子群算法跑50个粒子迭代100代那就是500万次潮流计算量完全不可接受。所以场景必须削减。削减的思路是从1000个原始场景中选出M个有代表性的场景M通常取515并给每个场景赋一个概率。常用算法有两种同步回代削减法backward reduction和K-means聚类法。同步回代的核心思想是贪心删除——每次找与其他场景距离之和最小的场景删掉把它的概率累加到最近邻场景上直到剩下M个。K-means更简单直接把所有原始场景聚成M类用聚类中心作为典型场景用类内的样本占比作为场景概率。Matlab里K-means实现非常直接% 生成 1000 个光伏出力场景每个场景 24 个时段 N 1000; M 10; S betarnd(alpha, beta, N, 24); % 按Beta分布抽样 % 归一化到光伏额定容量 S S * P_pv_rated; % K-means削减到10个典型场景 [cluster_index, centers] kmeans(S, M, Distance, sqeuclidean, Replicates, 20); prob accumarray(cluster_index, 1) / N;Replicates设为20是为了防止K-means陷入局部最优初始化跑20次取最优聚类结果。削减后centers就是典型场景的光伏出力曲线prob是对应的概率。这样计算量就从1000次潮流降到了10次目标函数评估速度快了100倍。这里有一个经验场景数不是越多越好。10个场景和15个场景的优化结果差异通常很小但计算时间差了50%。如果做灵敏度分析或参数扫描我一般把场景数压到812个既能保留光伏出力的分布特性又不至于让仿真等得心焦。2.3 波动性进入目标函数的方式期望目标与越限惩罚场景生成好之后难点是怎么在目标函数里把波动性体现出来。最简单的方式是用所有场景的期望指标。比如网损期望值f1 Σ_s prob(s) · P_loss(s)也就是说对一组给定的控制变量先逐个场景调用潮流计算算出每个场景下的网损和电压偏差再按场景概率加权求和。这组控制变量在所有典型场景下的平均表现就得到了光伏波动性自然被“计及”了进去。但只算期望值有一个问题期望值会把极端场景“平均掉”。比如10个场景里9个电压正常1个场景末端电压到了1.07这个越限情况对期望电压偏差的影响可能并不大实际运行却可能触发保护跳闸。因此我会在目标函数或约束条件里加一层“越限惩罚”统计所有场景中的电压越限情况并乘一个大系数加入目标。做法也很直接viol_prob 0; for s 1:M V runPowerFlow(x, scenario(s)); viol_s max(0, V - V_max) max(0, V_min - V); if any(viol_s 0) viol_prob viol_prob prob(s); end end f1 f1 lambda_penalty * viol_prob;这样算法在寻优时就会主动避开那些会让电压在某些场景下越限的方案优化结果更贴近实际运行要求。如果你做的是机会约束规划可以把viol_prob大于给定阈值比如5%的判断直接写成约束但罚函数法实现起来更简单调参也直观。3. 多目标粒子群优化算法的核心设计3.1 为什么单目标加权法在这里不够用很多第一次接触多目标优化的人会问为什么不把网损和电压偏差乘上权重相加变成一个目标函数我之前也这么干过但实际做下来至少有三个问题。第一个是权重的主观性。网损单位是千瓦电压偏差单位是标幺值两者的数量级差很多不归一化的话权重几乎没有物理意义即便归一化权重怎么取还是靠试。第二个更致命——加权法遇到非凸Pareto前沿时有些最优解无论怎么取权重都找不到。第三配电网里OLTC和电容器是离散变量目标函数本身不光滑加权法很容易卡在局部解附近。所以这个项目用多目标粒子群优化算法MOPSO不是炫技而是因为它能一次性给出整条Pareto前沿。Pareto前沿上的每个点都代表一个“不能在不牺牲一个目标的前提下改善另一个目标”的折中方案。做研究也好做工程方案也罢拿到的是一组候选方案而不是被一个权重锁死的唯一解后续决策空间大很多。下面简单对比一下两种思路对比项单目标加权法MOPSO结果单个解一组Pareto解集权重需人工设定敏感不需要解集中自动覆盖折中非凸前沿可能丢失解可以覆盖计算量较小较大但可控结果决策无选择空间可结合偏好选折中解3.2 Pareto支配与外部档案机制再解释一下Pareto支配用大白话说就是方案A在所有目标上都不比方案B差并且至少有一个目标严格比B好那A就支配B。所有不被任何其他方案支配的解构成Pareto前沿。MOPSO里有一个专门存放非支配解的地方叫外部档案Archive。每次迭代结束把当前种群中的非支配解塞进外部档案再剔除掉已经被其他解支配的旧解。但外部档案容量有限默认100个装满了怎么办就需要裁剪。裁剪的原则是保留“稀疏区域”的解删掉“拥挤区域”的解这样Pareto前沿才能分布均匀而不是堆在一小块区域。计算拥挤度的方法有很多我常用网格法把目标空间划分成网格统计每个网格里的粒子数粒子越多的网格表示这一带的解越拥挤优先从这些网格里删。这相当于在目标空间里做了一次密度估计实现简单效果也不错。3.3 速度更新、离散变量处理与变异机制MOPSO的粒子更新公式和标准PSO一样v(t1) w·v(t) c1·r1·(pbest - x) c2·r2·(gbest - x) x(t1) x(t) v(t1)w是惯性权重c1和c2是学习因子r1和r2是[0,1]均匀随机数。区别在于gbest怎么选单目标PSO里gbest是全局最优解MOPSO里没有单一的最优解所以gbest从外部档案里选而且倾向于从粒子密度最低的网格里选目的是引导粒子去探索Pareto前沿上还比较“空旷”的区域。配电网优化里有离散变量OLTC档位和电容器投切组数本身是整数。粒子位置更新后是连续值必须做取整处理。但直接四舍五入会导致粒子在某个档位附近反复跳动影响收敛因此我通常还加一个边界检查和邻域扰动。如果粒子越界把它拉回边界如果连续多次收敛到同一个离散值附近就随机探索相邻档位避免死锁。变异机制也是MOPSO保持多样性的关键。标准PSO没有变异MOPSO在早期容易聚到局部Pareto前沿。我习惯让变异率随时间递减p_m(t) 0.2 - 0.15 * (t / T_max)迭代初期变异率高粒子探索范围大后期变异率降低收敛到局部精细搜索。这个机制对配电网这种多局部极值的问题非常有效。核心更新代码骨架如下for i 1:pop 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, discreteIdx) round(x(i, discreteIdx)); % 边界检查 x(i,:) max(min(x(i,:), ub), lb); % 变异 if rand p_m idx randi(nVar); x(i, idx) lb(idx) rand * (ub(idx) - lb(idx)); if ismember(idx, discreteIdx) x(i, idx) round(x(i, idx)); end end end这段代码是整个MOPSO的核心后面接上潮流计算和目标函数评估就是一个可以跑通的闭环。4. Matlab代码实现从主程序到潮流计算4.1 项目目录与主程序流程Matlab代码的组织直接影响调试效率。我现在的习惯是把功能拆开每个文件只做一件事目录结构大致是这样moPso_ADN/ ├─ main.m % 主程序设置参数并调用优化 ├─ case33.m % 配电网基础数据33节点系统为例 ├─ loadData.m % 负荷数据、光伏场景数据 ├─ genScenarios.m % 生成光伏场景并削减 ├─ moPSO_Core.m % MOPSO主算法 ├─ powerflow.m % 前推回代潮流计算 ├─ objectiveFun.m % 目标函数评估 ├─ plotResults.m % 绘制Pareto前沿和节点电压分布main.m只做五件事加载数据、生成场景、初始化MOPSO参数、调用优化循环、输出和画图。结构清晰之后调参和加功能都很方便。主程序的核心流程大致是% main.m clc; clear; close all; load(case33.mat); % 网络阻抗、拓扑 [Scenario, Prob] genScenarios(1000, 10); % 生成10个典型场景 % M0PSO 参数 pop 50; maxGen 100; archiveSize 100; nGrid 20; [front, bestSolution] moPSO_Core(...); plotResults(front, bestSolution);4.2 前推回代潮流计算配电网版本的“迭代收敛”配电网优化里潮流计算是内层循环每个粒子每次迭代都要调用很多次因此它的速度决定了整个算法能不能跑完。33节点这类辐射状配电网我强烈推荐前推回代法而不是牛顿拉夫逊法。原因有两点前推回代不需要求导和形成雅可比矩阵每步只是简单的复数运算速度快辐射状网结构天然适合前推回代收敛性也很稳定。前推回代的基本思路分两步回代从末端节点向首端节点根据当前节点电压和负荷/电源注入功率计算每条支路的电流或功率。前推从首端节点向末端节点根据支路电流和线路阻抗更新各节点电压。两步骤反复迭代直到前后两次节点电压幅值差的最大值小于收敛阈值比如1e-6。Matlab代码骨架function [V, P_loss] powerflow(branch, bus, load_power, gen_power) n size(bus, 1); V ones(n, 1); I_branch zeros(size(branch, 1), 1); for iter 1:100 V_old V; % 回代从末端往首端算支路电流 for k size(branch,1):-1:1 from branch(k,1); to branch(k,2); S_inj load_power(to) - gen_power(to); I_branch(k) conj(S_inj / V(to)); % 累加支路下游电流略实际需拓扑排序 end % 前推从首端往末端算节点电压 for k 1:size(branch,1) from branch(k,1); to branch(k,2); V(to) V(from) - branch(k,3) * I_branch(k); end if max(abs(V - V_old)) 1e-6 break; end end P_loss sum(abs(I_branch).^2 .* branch(:,3)); end实际工程代码要考虑拓扑顺序先做广度优先搜索给支路编号然后按逆拓扑序回代、按拓扑序前推。前推回代对纯辐射状网收敛很快一般迭代几十次就达到1e-6。如果网络里有弱环需要先开环再补偿会稍微复杂一点但主流配电网优化算例基本都用辐射状网。光伏逆变器有剩余无功能力时我会让光伏节点在潮流中按“恒无功注入”处理无功值由决策变量给出如果某些场景下节点电压越限再把无功修正为维持电压的PV控制模式。两块逻辑加在一起仿真结果会更贴近真实逆变器的控制策略。4.3 目标函数与罚函数让搜索“听话”objectiveFun.m是核心评估函数输入是一个粒子的决策变量输出是两个目标函数值和违约量。大致的流程是function [f1, f2] objectiveFun(x, Scenario, Prob, data) % 解析决策变量 tap x(1); cap x(2); dg_q x(3:end); % 初始化目标 f1 0; f2 0; for s 1:length(Prob) % 根据场景修改光伏出力 data.gen(:,2) Scenario(s,:); % 潮流计算 [V, P_loss] powerflow(data, tap, cap, dg_q); % 电压偏差 v_dev sum((V - 1).^2); % 场景加权累加 f1 f1 Prob(s) * P_loss; f2 f2 Prob(s) * v_dev; end % 罚函数处理电压越限 f1 f1 1000 * max(0, max(abs(V-1) - 0.05)); f2 f2 1000 * max(0, max(abs(V-1) - 0.05));罚函数系数取1000是我常用的值。太小了惩罚不够越限方案可能留在Pareto前沿上太大了会压制其他目标的梯度信息导致粒子早早丢失多样性。1000在大多数量纲下都够用如果你发现优化结果中有大量越限解可以把它提到1e4。4.4 参数设置速查表与归一化细节MOPSO参数对收敛速度和Pareto前沿质量影响很大给一张我调通这类型项目的速查表参数推荐值说明种群规模 pop50 ~ 100决策变量多时取大最大迭代次数 maxGen100 ~ 300看收敛曲线再调整惯性权重 w0.9 → 0.4 线性递减前期全局探索后期局部搜索学习因子 c1、c22.0、2.0常见配置也可略做调整外部档案容量100太大拥挤度计算变慢网格划分数20 × 20按目标个数调整变异率初值/终值0.2 → 0.05线性递减这里特别要强调归一化。外部档案的拥挤度计算依赖目标值之间的距离网损可能是几千瓦电压偏差可能是零点几甚至几点几如果不归一化拥挤度全被网损主导电压偏差的多样性就丢了。我在moPSO_Core里对每个目标都做了一遍min-max归一化再算网格和拥挤度Pareto前沿立刻就均匀了。这一步非常关键算是踩了坑才得出的经验。5. 调试实录常见问题与排查技巧5.1 结果不收敛、Pareto前沿稀烂时先排查这三件事这个项目调试过程里最折腾的问题是算法跑完Pareto前沿稀稀疏疏甚至只有几个点。后来总结出三个排查顺序。第一先跑单场景再跑多场景。光伏场景一多每个粒子都要对所有场景做潮流计算量大不说出了问题根本不知道是潮流算错了还是场景概率算错了。我的一般做法是先把场景数临时改成1跑确定性优化确认潮流、罚函数、粒子更新这些环节在单场景下都能正常收敛再切回多场景模型。第二检查潮流迭代次数和残差。前推回代如果不收敛通常是因为设定的迭代上限太小或负荷/电源数据有误导致功率严重不平衡。在powerflow里加一行disp(iter, max(abs(V - V_old)))观察残差变化一秒就能定位。第三检查目标函数量纲和归一化。如果f1和f2的量级差一个数量级以上Pareto前沿就会偏向某一侧。先把两个目标都归一化到[0,1]再跑MOPSO前沿分布会正常很多。5.2 离散变量导致振荡换挡永远在犹豫OLTC变比和电容器投切组数是离散变量粒子速度更新后四舍五入取整这是最简单的处理方式但有个副作用粒子可能在两个相邻档位之间反复横跳尤其是迭代后期速度已经很小的时候四舍五入的取整操作会把粒子“踢”到相邻档位目标函数值随之跳变外部档案里的非支配解也跟着震荡。我的解决办法是分两步走。第一步在粒子更新流程里离散变量取整后加一个邻域概率扰动。代码里我常写if rand 0.1 x(i, discreteIdx) x(i, discreteIdx) randsample([-1 0 1], 1); x(i, discreteIdx) round(max(min(x(i, discreteIdx), ub(discreteIdx)), lb(discreteIdx))); end第二步外部档案里保存的解全部保证是合法取整后的结果绝不允许把未取整的连续值放进档案。否则最后输出的“最优解”可能对应不存在的OLTC档位工程上完全不可用。做完这两步离散变量振荡的问题基本就消失了。5.3 电压越限概率与折中解选取多目标优化得到Pareto前沿之后怎么从100个解里选一个如果做工程报告我喜欢用模糊隶属度法。先找出每个目标在所有非支配解中的最大值和最小值然后计算每个解的隶属度μ_i (f_i_max - f_i) / (f_i_max - f_i_min)把每个解的所有目标隶属度累加取总隶属度最大的那个作为折中解。这个方法的直观含义是挑一个在网损和电压偏差两个目标上相对最接近“各目标最好值”的方案。代码很简单fmax max(front); fmin min(front); mu (repmat(fmax, size(front,1), 1) - front) ./ (repmat(fmax - fmin, size(front,1), 1)); total_mu sum(mu, 2); [~, bestIdx] max(total_mu); bestSolution front(bestIdx, :);选完折中解后我还会在全部场景下重新运行一遍潮流统计每个节点的电压越限概率画成箱线图或者“电压区间包络线”。这个结果放到论文或报告里非常有说服力比单纯贴一张Pareto前沿散点图更能说明“计及光伏波动性”的效果。画图用Matlab的fill函数画包络带效果不错。% 统计各节点电压分布 V_all zeros(numStatic, numScenarios); for s 1:numScenarios V_all(:, s) runPowerFlow(bestSolution, Scenario(s)); end V_min min(V_all, [], 2); V_max max(V_all, [], 2); fill([1:n; n:-1:1], [V_min; flip(V_max)], [0.9 0.9 0.9]);5.4 场景削减与计算时间的平衡最后聊一个容易被忽略的细节场景削减比例和计算时间的平衡。如果你的算例是IEEE 33节点系统10个场景100个粒子跑100代单次潮流计算大概几毫秒总体时间还能接受。但如果换成大规模配电网比如300节点以上就必须考虑优化计算效率了。我有两个实操建议。第一把场景削减后的场景数控制在58个优先用K-means不要用同步回代——同步回代在样本量大时非常慢K-means在Matlab里有内置函数向量化运算快得多。第二潮流计算做矩阵化预处理把支路阻抗、拓扑结构提前展开成矩阵前推回代改成矩阵运算而不是for循环速度能提升好几倍。本质上做工程仿真精度和耗时的平衡一直是重点。最后说一点个人经验这类项目最容易出现的认知误区是把算法当主角把配电网模型当配角。事实上模型建模对了场景生成合理了哪怕用标准MOPSO也能得到不错的结果反过来模型和场景一团糟再高级的算法也救不回来。我自己的习惯是先把单场景确定性问题调通然后再引入波动性接着再上多目标算法每一步过了再往后走。这个顺序看起来“慢”但实际是整个项目里省时间最多的做法。如果你也在折腾类似的方向建议试试这个思路应该能少走不少弯路。