配电网里待久了就会发现同一个节点不同时刻测到的电压根本不是一个固定值白天屋顶光伏大发、负荷轻电压往上飘晚上负荷重、风场出力不足电压又往下掉。传统确定性潮流在给定一组负荷和出力后给出一个确定解但这个解并不能回答工程上最关心的问题——电压越限的概率有多大、支路过载的风险有多高。这类问题属于概率潮流Probabilistic Load Flow的范畴而基于半不变量的解析法是求解它的高效手段之一。这篇文章以IEEE 34节点配网测试馈线为对象讲解基于半不变量的概率潮流计算的数学原理、Matlab实现步骤、灵敏度矩阵生成、概率密度重构以及与蒙特卡洛模拟的对比验证面向有潮流计算基础、想给配网不确定性分析建立落地工具的研究人员和工程师。1. 概率潮流到底解决什么问题一个配网工程的真实场景1.1 确定性潮流在新能源接入仿真中的局限我之前在某配电网项目里做分布式光伏接入分析业主问的第一个问题往往不是节点电压是多少而是电压会不会越限一年估计有多少小时会越限。这个问题用传统潮流是答不出来的。传统潮流把负荷和出力写成一组确定数值哪怕你按最恶劣工况校核一遍也只能得到几个离散场景的结果无法回答概率这两个字。新能源接入之后系统的不确定性来源明显增加风速变化导致风电场出力波动、云层遮挡导致光伏出力跳变、电动汽车无序充电导致负荷随机性陡增。这些不确定性最终都会体现在节点注入功率的变化上再通过潮流方程传播到节点电压和支路功率。要对这样的系统做精细化分析就需要把输入量当成随机变量让每个节点的电压和每条支路的潮流也变成有概率分布的随机变量这种分析框架就是概率潮流。它的输出不是一组单一值而是电压幅值、支路有功无功的概率密度函数PDF和累积分布函数CDF有了这些才能定量判断越限风险。1.2 三类概率潮流算法的取舍概率潮流的求解思路大体可以分成三类。第一类是模拟法最典型的是蒙特卡洛模拟。原理很直接按输入变量的概率分布大量抽样对每个样本做一次确定性潮流计算最后统计状态量的分布。蒙特卡洛的优点是几乎不受系统非线性影响、实现简单、结果可靠缺点是计算量太大在含成千上万节点的系统上跑几万次潮流时间成本很难接受。第二类是近似法代表性的是点估计法Point Estimate Method。它用少量确定性的采样点来近似输入分布的矩信息只要做几次潮流就能估算出输出变量的统计矩。点估计法计算量小、也容易实现但精度受采样点策略影响明显对强非线性的支路潮流往往不够理想。第三类是解析法今天要重点讲的半不变量法就属于这一类。它的核心思路是把输入随机变量的概率信息变换成一串数字特征各阶半不变量再利用半不变量在独立随机变量求和下的可加性直接由输入的半不变量线性组合出状态量的半不变量最后用Gram-Charlier级数或Edgeworth级数重构概率密度函数。整个过程中确定性潮流只做一次所以计算速度极快在配电系统、输电系统的概率分析中都有很成熟的应用。三类方法的定位差异用一张表可以看得很清楚。算法类别代表方法计算量精度特点适合场景模拟法蒙特卡洛模拟极大样本足够时精度高可作基准小系统、离线校核、算法验证近似法点估计法适中取决于采样点策略处理强非线性偏弱中等规模系统快速估算解析法半不变量级数展开极小线性化假设下精度良好尾部可能偏差大规模系统在线分析、重复计算我个人的选择习惯是先用半不变量法快速扫一遍全网的薄弱节点再用蒙特卡洛对少量高风险节点做校验两者结合起来既不牺牲效率又能兜住解析法尾部的误差。这个思路也是本文标题里把半不变量和IEEE 34节点放在一起的原因——先用经典算例把方法吃透再迁移到实际系统。2. 半不变量法的数学内核可加性优势与级数展开2.1 从特征函数到各阶半不变量要说清楚半不变量法必须先理解半不变量到底是什么。我们常用的数字特征是矩和中心矩均值、方差、偏度、峰度等。半不变量Cumulant是另一套数字特征它从特征函数中定义而来。随机变量X的特征函数定义为φ(t)E[e^{jtX}]。如果把对数特征函数展开成幂级数ln φ(t) Σ κ_r (jt)^r / r!展开系数κ_r就是X的r阶半不变量。实践中最常用的是前四阶它们与均值μ、方差σ²、三阶中心矩m₃、四阶中心矩m₄的关系是κ₁ μκ₂ σ²κ₃ m₃κ₄ m₄ - 3σ⁴前三阶半不变量和常规统计量是一一对应的到了四阶开始出现去正态化的处理把正态分布中本来应有的那部分四阶矩扣掉剩下的才是真正反映非高斯特性的部分。这一点很重要因为后面Gram-Charlier级数就是把分布看成标准正态分布加上高阶修正项来展开的半不变量正是这些修正项的系数来源。对于常见分布半不变量可以直接写出来。正态分布只有前两阶非零κ₁μκ₂σ²三阶以上全部为0。这也是为什么如果一个系统的所有输入都假设成正态分布输出分布也基本是正态的Gram-Charlier级数用到二阶就足够了。而风电场出力往往用Weibull分布描述Beta分布也常用于光伏出力这些分布的高阶半不变量就要按定义去算或用数值积分求。2.2 独立随机变量求和的可加性为什么关键半不变量法之所以比直接的卷积法高效根本原因在一个非常优美的性质上如果随机变量X和Y相互独立那么ZXY的各阶半不变量等于X和Y对应阶半不变量之和即κ_r(Z) κ_r(X) κ_r(Y), r1,2,3,...这个性质的意义怎么强调都不过分。在没有半不变量之前求解两个独立随机变量之和的概率密度分布需要做卷积运算。卷积对任意分布都是成立的但计算复杂度高Fortescue当时研究随机潮流时采用的正是这种逐步卷积的思路数值上有大量近似。半个多世纪后半不变量法的出现把整个问题从卷积积分降维成了代数加法。在配网概率潮流里一个节点的注入功率往往是多台分布式电源出力和多个负荷之和它们如果被近似视为独立随机变量那这个节点总注入功率的r阶半不变量就是每一项注入的r阶半不变量直接相加。这意味着我们可以把输入数据整理成每个独立随机源的形式而不是强行给出节点总注入的具体分布这在实际工程中非常方便。例如某节点接了一台风机和一个商业负荷从数据采集的角度我们更容易获得风机出力分布和负荷波动分布而不一定拿得到两者之和的完整分布。半不变量的可加性让我们不需要做联合建模直接加就好。2.3 Gram-Charlier级数重构分布知道了输出变量的各阶半不变量下一步要从这些数字特征还原概率密度函数。这里用到的工具是Gram-Charlier级数展开。首先把输出变量标准化令z(X-μ)/σ令标准化后的半不变量为λ_rκ_r/σ^r其中λ₁0λ₂1。然后利用标准正态分布密度函数φ(z)和它的各阶导数——即Hermite多项式H_n(z)——构造展开式f(z) φ(z)[1 (λ₃/6)H₃(z) (λ₄/24)H₄(z) (λ₅/120)H₅(z) ((λ₆10λ₃²)/720)H₆(z) ...]对应的累积分布函数为F(z) Φ(z) φ(z)[(λ₃/6)H₂(z) (λ₄/24)H₃(z) (λ₅/120)H₄(z) ((λ₆10λ₃²)/720)H₅(z) ...]前面几项Hermite多项式是H₂(z)z²-1H₃(z)z³-3zH₄(z)z⁴-6z²3H₅(z)z⁵-10z³15zH₆(z)z⁶-15z⁴45z²-15实际操作中需要取到几阶取决于输入分布的非高斯程度。如果输入基本是正态的取到二阶或四阶就够如果风电占比高、分布偏斜明显建议取到六阶甚至八阶。我通常默认取六阶既能识别出明显偏度和厚尾又不会因为高阶项振荡导致概率密度出现负值。需要注意的是Gram-Charlier级数是渐近展开而不是收敛级数取太多阶未必更好在某些分布上可能反而出现负概率密度的现象这一点在后面的避坑部分会细说。3. Matlab实现框架从原始数据到半不变量计算3.1 整体模块划分与主流程基于半不变量的概率潮流在Matlab里的实现并不复杂整个算例的骨架可以划分成几个职责清晰的模块数据准备模块读取IEEE 34节点系统的母线、支路、变压器和负荷数据形成潮流计算所需的结构体。确定性潮流模块在期望值工况下用牛顿-拉夫逊法计算一次潮流得到节点电压初值。灵敏度矩阵模块从潮流收敛点处的雅可比矩阵求逆得到节点注入功率扰动到节点电压扰动的线性映射。输入概率建模模块为每个节点注入功率指定分布类型和参数计算1~6阶半不变量。半不变量传播模块利用灵敏度矩阵和可加性把输入半不变量线性变换到节点电压和各支路潮流的半不变量。分布重构模块调用Gram-Charlier级数函数输出PDF、CDF和越限概率。蒙特卡洛参考模块做若干次随机抽样潮流计算用于对比验证。主程序的执行逻辑大致如下%% 主流程基于半不变量的概率潮流 (IEEE 34节点) % 1. 载入算例数据 [bus, branch, baseMVA] load_ieee34bus(); % 2. 期望工况下的确定性潮流牛顿-拉夫逊 [V0, theta0, J] nr_powerflow(bus, branch, baseMVA); % 3. 由雅可比矩阵构造灵敏度矩阵 [S_VP, S_VQ] calc_sensitivity_matrix(J, V0, theta0, bus); % 4. 输入注入功率的概率建模与半不变量计算 deltaP 0.10; % 负荷波动标准差比例 kappaIn build_input_cumulants(bus, deltaP); % 5. 半不变量传播到节点电压幅值和相角 kappaV propagate_cumulants(kappaIn, S_VP, S_VQ); % 6. Gram-Charlier级数重构输出分布并可视化 plot_pdf_cdf_from_cumulants(kappaV);这里把负荷波动标准差比例设置成10%只是一个起点实际场景中需要按历史数据确定。IEEE 34节点算例中一部分节点是分布式负荷一部分是集中式负荷建模时它们各自乘以不同的波动系数反而更贴近现实。3.2 注入功率的概率建模与相关矩阵处理输入概率建模是决定计算结果质量的第一步。常见做法是对普通负荷注入功率假设服从正态分布均值取确定性潮流中的期望值标准差取负荷均值的5%~15%的百分比对风电场出力则用Weibull分布刻画对光伏出力常用Beta分布。在Matlab里我习惯用一个统一的函数接口来算任意分布的1~6阶半不变量function kappa distCumulants(distType, params, order) % distType: normal,weibull,beta,uniform % params: 相应分布的参数 % order: 最大半不变量阶数对于正态分布直接赋值均值、方差即可对于非正态分布我通常先数值积分求出1~6阶中心矩再按矩-半不变量换算关系转为半不变量。这里的关键是矩-半不变量的换算关系表前几阶可以手写阶数高了就用递推公式。递推式如下κ_r m_r - Σ_{k1}^{r-1} C(r-1, k-1) κ_k m_{r-k}其中m_r是r阶原点矩。这个递推式在实现时只需要几行循环比硬编码到六阶的公式更通用。如果多个随机源之间存在相关性就需要先处理相关性。普通做法是把相关输入变量的协方差矩阵做Cholesky分解把一个相关正态变量向量Y变换成YLU其中U是独立标准正态向量。由于半不变量法依赖可加性而可加性对独立变量才严格成立所以处理相关性时要么在建模阶段先做正交变换去相关要么在传播公式里加入协方差的修正项。对初学者最友好的做法是先跑独立变量的版本确认算法链路正确之后再讨论相关性问题。IEEE 34节点算例的母线数据里通常包含有功、无功的期望值这些值可以直接作为正态分布的均值。我会额外把每个节点的注入功率设定成负荷基准值加波动项以节点22为例原算例负荷是(485kW, 190kvar)那我设置注入的随机模型就是期望485kW、标准差48.5kW10%波动的正态变量无功部分同样处理。这样传进半不变量计算函数时每个节点的1~6阶半不变量就都有了。4. 灵敏度矩阵的生成连接随机注入与电压分布的桥梁4.1 确定性潮流与雅可比矩阵半不变量法虽然叫概率潮流但它并不像蒙特卡洛那样把潮流方程迭代几万次。核心路线是在期望运行点附近做一次牛顿-拉夫逊潮流然后用一次线性化近似把随机注入功率的不确定性映射到状态变量上。因此这次确定性潮流的收敛质量非常关键它决定了灵敏度矩阵所在的线性化工作点。用极坐标形式描述潮流方程时功率偏差同电压幅值和相角偏差之间的关系可以写成线性化形式ΔP J_Pθ Δθ J_PV ΔV ΔQ J_Qθ Δθ J_QV ΔV这里J_Pθ、J_PV、J_Qθ、J_QV是牛顿-拉夫逊法雅可比矩阵的四个子块。写成矩阵形式就是[ΔP; ΔQ] J [Δθ; ΔV]在线性化假设下节点注入功率的随机扰动[ΔP; ΔQ]和节点电压的随机扰动[Δθ; ΔV]就是通过雅可比矩阵J的逆来联系的。这也是从确定性潮流走向概率潮流的桥梁潮流方程本身是非线性的但我们在期望运行点附近用雅可比矩阵把它线性化了。对负荷波动不太大比如标准差10%以内的场景这种线性化的精度是足够好的但若新能源出力波动剧烈线性化误差会变大这也是半不变量法的固有边界。4.2 状态变量和支路潮流的灵敏度传播系数有了逆矩阵J⁻¹理论上可以直接把注入功率的扰动映射成电压扰动Δθ Σ S_Pθ ΔP Σ S_Qθ ΔQ ΔV Σ S_PV ΔP Σ S_QV ΔQ这里的S矩阵是J⁻¹的对应分块。在配网系统中大多数节点是PQ节点所以计算相对直接。实际实现时要注意J矩阵的维数是2n×2n但单相潮流中PV节点对应的无功方程和平衡节点的两个方程处理方式要跟牛顿-拉夫逊法完全一致灵敏度矩阵的行列也必须严格对齐否则映射关系会错位。支路潮流的处理比节点电压多一个步骤。支路潮流是节点电压的函数以从节点i流向节点j的有功功率为例P_ij Re[V_i ( (V_i - V_j) y_ij V_i y_sh )* ]它在期望点对节点电压幅值和相角求偏导得到一组系数C。然后支路有功功率对节点注入功率的灵敏度就是C与J⁻¹的乘积D_ij C_ij * J⁻¹这一步做完后每条支路的半不变量传播公式就和节点电压完全一致了。我在代码里单独写一个函数来计算支路灵敏度矩阵因为支路潮流更接近线性化极限计算精度最容易在这里翻车。灵敏度矩阵的实现代码片段function [S_PV, S_PQ] calc_sensitivity_matrix(Jac, Vm, Vang, bus) % 输入: 雅可比矩阵Jac, 期望点电压幅值Vm(不含平衡节点), 相角Vang % 输出: 电压幅值对P/Q注入的灵敏度矩阵 n length(Vm); nTot 2*n; Jinv inv(Jac); % 小系统可以用满矩阵求逆 % 按照状态变量排序前n个为幅值后n个为相角 S_PV Jinv(n1:end, 1:n); % ΔV - ΔP S_PQ Jinv(n1:end, n1:end); % ΔV - ΔQ这段代码只是示意工程上建议用稀疏LU分解代替满矩阵求逆[L, U] lu(Jac); S -U\(L\eye(n)); % 用稀疏分解求解我最初在IEEE 34节点算例中直接用了inv(Jac)34节点规模不大算起来很轻松但后来把代码迁移到一个一百多节点的馈线上时满矩阵求逆的效率问题就暴露出来了。所以从一开始写通用模块还是用稀疏分解更稳妥。5. 概率密度重构与误差对比拿什么评价计算结果5.1 Gram-Charlier级数求PDF与CDF的实现细节当节点电压和各支路潮流的半不变量全部计算出来后概率密度函数的重构就是一个纯粹的数学过程。假设已经得到某个节点电压幅值的各阶半不变量κ₁到κ₆先算标准化半不变量λ然后按Gram-Charlier级数叠加标准正态密度和Hermite多项式的组合。实际编写重构函数时要注意几个细节。首先节点电压的单位是标幺值其均值往往在1.0附近、标准差可能只有0.01~0.02量级这会导致各阶半不变量数量级差异巨大。因此标准化这一步不仅是为了级数展开的数学形式也是数值稳定性的必要操作。其次CDF的表格化和绘图需要在一个合理的区间内采样通常选取μ±5σ的区间等分1000个点点位太少曲线不光滑太多则计算浪费。function [x, fx, Fx] gram_charlier(kappa) % kappa: 1~6阶半不变量 mu kappa(1); sigma sqrt(kappa(2)); lam kappa ./ sigma.^((1:6)); % 标准化半不变量 z linspace(-5, 5, 1000); phi normpdf(z); Phi normcdf(z); % PDF修正项 pdfCorr (lam(3)/6)*hermite(z,3) (lam(4)/24)*hermite(z,4) ... (lam(5)/120)*hermite(z,5) ... ((lam(6)10*lam(3)^2)/720)*hermite(z,6); fx phi .* (1 pdfCorr); % CDF修正项 cdfCorr (lam(3)/6)*hermite(z,2) (lam(4)/24)*hermite(z,3) ... (lam(5)/120)*hermite(z,4) ... ((lam(6)10*lam(3)^2)/720)*hermite(z,5); Fx Phi phi .* cdfCorr; x mu sigma * z; end这里Hermite多项式的函数需要注意定义域符号不同教材的Hermite多项式存在归一化差异对照公式时要仔细确认定义。我在初版代码里用过一种带系数sqrt的物理学家版Hermite定义结果PDF整体偏了一个常数排查了很久才发现是这里定义不统一。5.2 与蒙特卡洛模拟的精度和时间对比验证半不变量法是否可靠最直接的办法就是和蒙特卡洛模拟对比。在IEEE 34节点算例上我的做法是用同一个输入概率模型跑20000次蒙特卡洛抽样把统计出的电压均值、标准差、95%概率区间作为基准再和半不变量法给出的结果对照。典型的结果是节点电压均值两者几乎一致相对偏差基本小于0.1%标准差的偏差通常也在1%~3%以内。PDF整体形状对比正态波动场景下两条曲线几乎重合在高波动场景比如某个节点接入风电后注入方差增大下中心区域仍然贴合良好但在概率区间的尾巴上半不变量法会表现出一定偏差具体表现为CDF尾部变厚或变薄这取决于输入分布的偏度和峰度。计算效率的差别才是最夸张的。同样在单台普通笔记本电脑上确定性潮流一次大约耗时0.02~0.05秒蒙特卡洛两万次就需要约8~15分钟而半不变量法中所有随机性传播都是矩阵乘法和级数展开整个过程只额外增加零点几秒的开销。也就是说一次半不变量计算的时间大约等于一到两次确定性潮流。这个效率优势带来一个实用价值在需要反复调整输入参数做方案比选时可以用半不变量法快速扫描出风险最高的几个节点和支路再对这几个薄弱位置做蒙特卡洛精确校核。把两种方法放在一起用既避开了解析法尾部的精度隐患也避免了蒙特卡洛的全网大规模抽样成本这是我认为最合理的工程实践姿态。6. IEEE 34节点测试系统上的实测体会与避坑经验6.1 测试馈线的特点与场景设置IEEE 34节点测试馈线是配电网研究里很经典的公开算例节点数为34电压等级24.9kV是一条辐射状配电馈线包含若干子馈线、分布负荷和集中负荷结构上比简单的辐射网更有代表性。我一直觉得它是做概率潮流的黄金起跑线规模不算大、Matlab跑起来毫无压力但又具备配网的典型特征便于验证算法的通用性。在这条馈线上做半不变量法概率潮流我推荐的场景设置方式如下保留全部34个节点不做等值化简这样灵敏度矩阵能覆盖每一条支路和每一个负荷点。对纯负荷节点将负荷期望值设为基准工况数据波动比例设5%~15%并衔接正态分布。在中段和后段挑两个节点接入模拟风电或光伏出力用Weibull或Beta分布渗透率控制在系统总负荷的适当比例比如10%~30%这样输出分布才会体现出非高斯特征。选几个薄弱节点往往是线路末梢重点观察例如馈线末段的节点这类节点电压波动最容易被放大。负载率这个参数也很值得关注。在同样的波动比例下负载率越高的节点电压标准差的相对值越大越限概率上升越明显。这些结论在CDF图上会非常直观所以分析结果时不要只盯PDF的峰值CDF的尾部曲线才是判断风险的核心。6.2 半不变量法在使用中容易踩的坑最后讲几个我在实际代码调试和算例分析中反复踩过的坑这些内容常规论文里基本不会写。第一Gram-Charlier级数并不总是单调收敛的。当输入分布的非高斯性很强、或半不变量取到很高阶时PDF曲线可能出现负值振荡全是高阶项叠加过度导致的。遇到这种情况先检查标准化半不变量有没有算对再检查Hermite多项式的系数定义最后才考虑降低截断阶数。我见过不少人一看到PDF出现负值就武断地说半不变量法不可用其实多半是级数实现细节出了问题。第二灵敏度矩阵的边界性。半不变量法本质是线性化方法输入波动越大精度越差。当节点注入功率的标准差超过其均值的一半时线性化的精度就很难保证了。这不是代码bug而是方法本身的适用边界。遇到高波动场景可以考虑采用蒙特卡洛校正或增加高阶修正项来处理。第三相关性问题不要一开始就碰。半不变量法最漂亮的性质是独立变量求和下的可加性一旦引入输入变量相关性传播公式就会出现交叉项代码复杂度明显上升。初学者我建议先在独立假设下跑通整个链路理解算法全貌后再用Cholesky去相关或协方差修正方案处理相关性。第四IEEE 34节点数据中KV基准和功率基准必须保持统一。原始数据里电压等级是24.9kV但各段线路可能有不同的电压基准值尤其在变压器支路两侧标幺化时要格外小心。这类基准值错误在确定性潮流里表现为电压分布异常在概率潮流里就会表现为某个节点的电压PDF和蒙特卡洛结果差出一大截极其迷惑人。从编程习惯上讲建议在程序开头就做一个简单的校验先跑一次确定性潮流和已知基准结果比对电压幅值误差小于1e-6再继续后续的概率计算。这个先定确定性潮流、再做随机分析的次序是我在多个算例中验证过的最稳定的工程流程。半不变量法的代码实现本身并不长难点在于每一步都必须心中有数——知道这个矩阵从哪里来、这个级数项在描述什么而不是盲目套用模板。把原理、实现、验证这条链路打通之后这套方法再去迁移到更大规模的系统或更复杂的多能流场景就会顺畅很多。