概率潮流计算实战：半不变量法原理与IEEE34节点Matlab实现

搞随机潮流这些年，我最常被问的一句话是：“为什么不能用确定性潮流加一个安全裕度搞定？”说实话，在新能源渗透率不高的时候，这么干确实够用；但等风电、光伏、充电桩都涌进来之后，单一工作点的潮流结果已经没法回答“电压越限概率到底是多少”这种问题了。所以我才认真去做了基于半不变量的概率潮流计算，而且选择在IEEE34节点系统上落地，用Matlab实现一套可复用的代码框架。

这套方案的核心并不复杂：先算一次基态潮流，然后借助灵敏度矩阵把节点注入功率的随机性映射到节点电压和支路潮流上，再用半不变量（也叫累积量）把随机分布的特征值提取出来，最后用Cornish-Fisher级数反推出电压、潮流的概率分布。相比蒙特卡洛那种需要上万次潮流计算的暴力做法，它的计算负担小得多，非常适合规划阶段反复比选方案。这篇文章就把原理、建模、代码实现、踩坑历程一次讲清楚。

1. 为什么我选择半不变量法做概率潮流

1.1 确定性潮流的局限在哪里

确定性潮流解出来的是一组确定的节点电压和支路功率，前提是负荷和发电都给定。可实际系统里，负荷曲线每天都在波动，风电场出力可能在一个小时内从满发跌到零，光伏中午和傍晚完全是两个状态。你拿一个确定性的重负荷场景做校核，只能知道“这个时刻能不能过”，却不知道“一天当中有多少时间会越限”。

我曾经拿某配电网的实际负荷曲线做过统计，馈线首端电流的波动范围能到额定值的±30%，有些节点电压的日波动甚至超过6%。这种情况下，确定性潮流算出来的电压可能是0.97 p.u.，看起来正常，但实际运行中可能每天都有几个小时电压降到0.95 p.u.以下。这个问题只有概率化分析才能暴露出来。

1.2 蒙特卡洛和解析类方法的取舍

概率潮流大致分三条路线：蒙特卡洛模拟、点估计法、解析法。

蒙特卡洛思路最直观，构造几万个随机场景，每个场景都做一次潮流计算，最后把结果统计成概率分布。它的优点是天然兼容非线性、相关性、复杂故障处理；缺点是计算量很大。在IEEE34这种几十节点的系统上，一万次潮流也就几秒钟到几十秒，感觉还好；但放在几百上千节点的实际配电网里，每次潮流都要迭代，一万次就会非常吃力。

点估计法只需要少量确定性计算，然后估算输出变量的各阶矩，它速度快，但高阶矩信息丢失严重，对重尾分布和严重越限场景的刻画不够细腻。

半不变量法属于解析法的一种，它利用随机变量特征函数取对数之后，把卷积运算变成加法运算，将各节点注入功率的随机信息凝聚成几个累积量。配合灵敏度矩阵，一步就能把输入分布的统计特性映射到输出变量上。整个过程只需要一次基态潮流，外加几次向量运算，计算量比蒙特卡洛低两到三个数量级。这就是我选它的核心理由。

1.3 半不变量法适合什么应用场景

做电网规划方案比选时，往往要反复修改电源接入位置、线路容量、无功配置。如果每个方案都跑一遍蒙特卡洛，时间成本太高。半不变量法每一轮计算几乎是瞬时的，适合放在优化循环里做迭代筛选。

它也适合用于运行风险评估，比如评估电压越限概率、支路过载概率、继电保护配合的置信区间。工程上不需要精确重现完整的联合概率密度，只要知道关键变量的期望、方差和若干分位数就足够了，而这恰好是半不变量加Cornish-Fisher展开最擅长的输出形式。

2. 半不变量的数学工具包

2.1 从随机变量到累积量的转化

为什么用半不变量而不用普通矩？因为普通矩有一个麻烦：独立随机变量之和的矩不等于各变量矩的简单相加。而累积量不一样，它具备可加性，这个特性让多节点注入功率的合成变得异常轻松。

累积量的定义来自特征函数的对数展开。设随机变量X的特征函数为φ(t)，则第r阶半不变量为：

[ \kappa_r = \frac{1}{j^r} \frac{d^r \ln \phi(t)}{dt^r}\bigg|_{t=0} ]

前几阶半不变量和中心矩的关系非常直观：

• κ₁ = 数学期望
• κ₂ = 方差
• κ₃ = 三阶中心矩，表征偏度
• κ₄ = 四阶中心矩减去3倍方差平方，表征峰度

实际编程时并不需要从特征函数求导，直接从样本或参数分布计算出均值和各阶中心矩，再换算成累积量就可以。对正态分布来说，三阶及以上的累积量全部是0；对威布尔分布这类偏态分布，高阶累积量则会不断提供偏度和峰度信息。

2.2 累积量的可加性原则

假设某个母线的注入功率由多个独立随机分量组成，例如风电场出力ΔPg、负荷波动ΔPl。那么注入功率ΔPi = ΔPg - ΔPl的前几阶累积量可以写成：

[ \kappa_{r}^{(\Delta P_i)} = \kappa_{r}^{(Pg)} + (-1)^r \kappa_{r}^{(Pl)} ]

这个公式就是整个半不变量法的基础。它表明，每个随机注入分量不需要做卷积，不需要通过采样合成，只需要把各自的累积量按阶数相加即可。累积量加法对计算效率的提升是决定性的。

不过有一点要提醒：可加性成立的前提是各分量相互独立。如果两个风电场距离很近、共享同一风区，或者多个负荷节点同步受气温影响，那累积量相加前必须先做相关性解耦。工程上常用的做法是先用Nataf变换把相关正态变量转换成独立正态变量，再做累积量叠加；相关性的处理会让代码复杂度上一个台阶，所以很多初版实现都不处理它。

2.3 Cornish-Fisher展开：把累积量变回概率分布

有了输入节点的累积量，经过灵敏度矩阵传递后，可以得到输出变量（如电压幅值）的累积量。但累积量本身不是概率分布，我们必须把它还原成累计分布函数或分位数。

Cornish-Fisher展开的思路是，从标准正态分布的分位数出发，用偏度、峰度等修正项逼近实际分布的分位数。四阶修正公式如下：

[ y_p \approx x_p + \frac{1}{6}(x_p^2 - 1)\frac{\kappa_3}{\sigma^3}

• \frac{1}{24}(x_p^3 - 3x_p)\frac{\kappa_4}{\sigma^4}

• \frac{1}{36}(2x_p^3 - 5x_p)\frac{\kappa_3^2}{\sigma^6} ]

其中x_p是标准正态分布对应概率p的分位数。如果只算到二阶修正，得到的就是正态近似；把三阶、四阶项加进去后，对偏态分布的效果明显改善。

我也试过Gram-Charlier级数，它在中心区域和Cornish-Fisher差别不大，但在尾巴位置容易振荡，出现概率小于0或大于1的非物理值。Cornish-Fisher在重尾情况下同样存在这个问题，但可以通过限制展开阶数和后期修正缓解。实际工程中一般取四阶到六阶就够用，取到八阶以上对数据噪声会很敏感，反而得不偿失。

3. IEEE34节点系统与计算参数

3.1 为什么选IEEE34作为测试对象

IEEE34节点测试馈线来自IEEE PES配电系统委员会，是一个接近真实工程的配电网络，基础电压等级24.9 kV，基准容量约2.5 MVA。它包含长距离线路、单相和三相混合馈线、分布式负荷和集中负荷、电压调节器、变压器等环节，比常见的IEEE30、IEEE118这类输电系统更贴近配电网的物理特征。

在随机潮流研究里选IEEE34有几个天然优势。一是网络不大不小，既能体现多节点注入对电压分布的交互影响，又不会因为节点太多导致灵敏度矩阵逆运算过于庞大。二是线路较长，电压对功率波动更敏感，概率潮流能看出明显的尾部差异。三是它自带负荷模型种类多，方便我们测试不同负荷波动假设对结果的影响。

3.2 负荷和分布式电源的不确定性建模

随机潮流的输入随机量主要是负荷和新能源出力，我在这里用最常用的两组模型：

负荷有功和无功功率采用正态分布，均值取潮流计算中的基准值，标准差按均值的5%到10%设置。假设负荷节点24的有功基准为0.184 MW，均值0.184、标准差0.012时，在95%置信区间下负荷会在0.160到0.208 MW之间波动。无功功率类似，取功率因数0.9左右。

风机出力采用两参数威布尔分布。风速v的概率密度为：

[ f(v)=\frac{k}{c}\left(\frac{v}{c}\right)^{k-1}\exp\left(-\left(\frac{v}{c}\right)^k\right) ]

风速的k阶原点矩有闭式解：

[ E[v^k]=c^k\Gamma\left(1+\frac{k}{\lambda}\right) ]

有了风速矩就能换算出风机出力的分布矩，再由矩转累积量，最终进入半不变量计算流程。我一般取形状参数λ=2.2，尺度参数c=7.5 m/s，并设定切入风速3 m/s、额定风速12 m/s、切出风速25 m/s。

3.3 灵敏度矩阵怎么取

半不变量法有一个隐含假设：输出变量对输入注入功率的变化是近似线性的。这个假设在潮流方程小扰动范围内是合理的。

设潮流方程写成F(x)=0，x包含节点电压幅值和相角。在基态运行点附近线性化得到：

[ J\Delta x = \Delta W ]

其中J就是牛顿拉夫逊法最后一次迭代得到的雅可比矩阵，ΔW是节点注入功率的随机扰动量。那么节点电压的随机扰动为：

[ \Delta x = -J^{-1}\Delta W ]

这里的灵敏度矩阵就是-J^{-1}。如果再需要支路潮流的概率分布，只需要把支路潮流对节点电压的偏导乘上去即可。

这里有个工程细节：雅可比矩阵中包含松弛节点对应的行和列，求逆之前我曾直接对整个矩阵求逆，结果在电压幅值分量上出现数值扰动。后来我改成先消去松弛节点相关的行列，再对缩减后的矩阵求逆，数值稳定性好很多。另外，如果系统包含PV节点，雅可比矩阵中对应的电压幅值行要替换为无功功率行，别弄混。

4. Matlab代码实现与核心模块

4.1 代码整体框架

我建议把整个实现拆成以下模块，方便后续扩展和复用：

prob_pf_main.m % 主程序：读数据、算潮流、算累积量、画图 load_ieee34.m % 读取IEEE34节点网络数据 newton_pf.m % 基态潮流 + 返回雅可比矩阵 calc_input_cumulants.m % 计算输入变量累积量 propagate_cumulants.m % 灵敏度矩阵传递累积量 cornish_fisher.m % Cornish-Fisher分位数还原 plot_pf_results.m % 绘制电压概率分布图

这里面最关键的是输入数据的组织。Matpower自带案例里没有IEEE34配电馈线，我最初用OpenDSS导出了TXT数据，再在Matlab里组装成bus和branch结构体；如果只想跑通算法，也可以把IEEE34的线路数据手工整理成表格，只保留正序单相等值参数。注意，IEEE34本身是三相不平衡馈线，严格做工程分析应该用三相潮流。本文为了聚焦半不变量算法，我把网络等效成单相正序模型，结果趋势有参考价值，但真实配网工程务必回到三相模型校验。

4.2 基态潮流与雅可比矩阵提取

基态潮流我用自己写的牛顿法，避免引入额外工具箱的接口限制。核心流程如下：

function [V, theta, J] = newton_pf(bus, branch, Ybus) n = size(bus, 1); V = bus(:, 8); % 初始电压幅值 theta = bus(:, 9) * pi / 180; tol = 1e-8; maxiter = 20; for k = 1:maxiter [F, J] = power_mismatch(V, theta, Ybus, bus); if norm(F, inf) < tol break; end dX = -J \ F; theta = theta + dX(1:n); V = V + dX(n+1:2*n); end % 返回最后一次迭代的雅可比矩阵 end

IEEE34线路的R/X比值普遍较高，初值选择对收敛影响很大。我用平启动时某些节点电压会震荡，后来把初始相角设为0、电压幅值设为1.02 p.u.，并限制每次迭代步长，收敛就稳定了。实际算例从开始到收敛约5到8次迭代。

雅可比矩阵提取时，要同时返回PQ节点对应的子块。对于随机潮流，我会把松弛节点的行和列消掉，因为松弛节点电压固定，不参与概率扰动。

4.3 注入功率半不变量计算

每个随机注入节点的注入功率定义为电源注入与负荷吸收的代数和。以负荷正态分布为例，直接用公式换算：

function kappa = calc_cumulants_normal(mu, sigma, norder) kappa = zeros(norder, 1); kappa(1) = mu; kappa(2) = sigma^2; % 其余高阶累积量为0 end

对于威布尔分布的风电出力，我先算出各阶原点矩，再转累积量。这里给出一个直接从样本构造累积量的通用函数，它适合任何分布：

function kappa = cumulants_from_samples(x, norder) mu = mean(x); xc = x - mu; kappa = zeros(norder, 1); kappa(1) = mu; kappa(2) = mean(xc.^2); kappa(3) = mean(xc.^3); kappa(4) = mean(xc.^4) - 3 * kappa(2)^2; kappa(5) = mean(xc.^5) - 10 * kappa(2) * kappa(3); kappa(6) = mean(xc.^6) - 15 * kappa(2) * kappa(4) ... - 10 * kappa(3)^2 + 30 * kappa(2)^3; end

如果需要更严格的原始矩转累积量公式，可以查任意数理统计教材里的连接多项式。实操中直接用这些公式配合数千次抽样的样本，就能得到稳定的累积量估计，比推导参数分布矩再转换更省事。

4.4 响应变量半不变量反推

假设输入随机变量ΔW经过灵敏度矩阵S映射到输出变量ΔX。在线性化假设下，输出变量的第r阶累积量为：

[ \kappa_r^{(\Delta X)} = \sum_{i=1}^{m} \left| S_{X,i} \right|^{r} \kappa_r^{(\Delta W_i)} ]

代码实现起来非常直接：

function kappa_out = propagate_cumulants(S, kappa_in, norder) n_out = size(S, 1); kappa_out = zeros(n_out, norder); for r = 1:norder % 对灵敏度取模值，再按阶数作幂运算 Sr = abs(S) .^ r; kappa_out(:, r) = Sr * kappa_in(:, r); end end

注意两个关键点。第一，灵敏度矩阵元素必须取绝对值，因为电压幅值对功率注入的灵敏度可能为负，若保留符号，高阶累积量的物理意义会错乱。第二，这个分段线性映射隐含了各输入随机变量相互独立的假设；如果输入间存在相关性，需要先做相关性变换后再进入这个环节。

4.5 最终概率结果输出

求得各电压幅值累积量后，Cornish-Fisher展开就可以输出任意分位数。我在代码中按以下方式实现电压越限概率计算：

function [prob_lo, prob_hi] = calc_violation(kappa, v_lo, v_hi) mu = kappa(1); sigma = sqrt(kappa(2)); g1 = kappa(3) / sigma^3; g2 = kappa(4) / sigma^4; % 求解分位点 z 对应的概率，用二分法或fzero p_lo = normcdf((v_lo - mu) / sigma); % 简单近似 % 考虑偏度修正时，把Cornish-Fisher反函数代入 ... end

更严格的做法是先对累积量做反Cornish-Fisher变换，求出累计概率。实际计算时，我会画出每个关键节点的电压概率密度曲线，曲线横轴为电压幅值p.u.，纵轴为概率密度。配电网关心的是电压在0.95 p.u.到1.05 p.u.之间的概率，低于下限或高于上限的概率就是越限风险，这也是规划人员最想要的单值指标。

5. 结果验证与常见坑

5.1 与蒙特卡洛对拍

半不变量法本质是近似算法，所以我始终保留蒙特卡洛作为校验基准。在IEEE34系统上，我用10000次蒙特卡洛模拟做对照，比较某个关键节点的电压分布。下表是一个典型算例的输出对比：

从结果看，前四阶累积量配合Cornish-Fisher展开已经能把越限概率估得很准。偏差主要来自潮流方程的线性化误差，还有高阶累积量截断误差。如果对尾部概率要求非常高，比如要评估万分之一的小概率事件，半不变量法会明显需要更高阶累积量和分段线性化修正。

5.2 结果出现非物理值怎么办

我最初跑出来的电压概率密度曲线在左尾出现了概率密度，这一点让同事一度以为代码有bug。排查后发现是两个原因叠加：一是Cornish-Fisher展开在极深尾部会出现振荡，导致概率密度变负；二是某些节点电压对注入功率的灵敏度系数被高估，线性化误差在重尾分布下被放大。

常用的处理手段是：限制输出分位数范围，将计算得到的分位数小于0.9 p.u.的部分用0.9 p.u.截断，或者改用Gram-Charlier级数并只保留四阶。更严谨的方案是对输入分布做截断，避免风功率出现负值等非物理采样，这样尾部振荡会缓解很多。

5.3 条件受限下的实用经验

先说相关性。我最早把所有负荷当成独立的，结果在多个负荷同步增长场景下严重低估了电压越限概率。后来我改用相关系数矩阵描述相邻负荷的同步性，并先用Cholesky分解生成相关正态向量，再转换成样本和累积量。这一步虽然增加了代码量，但结果可信度提升明显。

第二点是雅可比矩阵的质量。牛顿法迭代步数太少时，把未收敛的雅可比拿来做概率传播，误差会被后续步骤放大。建议至少让潮流残差降到1e-8以下，再提取雅可比矩阵。也别直接对全矩阵求逆，稀疏矩阵用LU分解代替inv，速度和内存都好得多。

6. 在实际项目里沉淀下来的细节

代码写完不等于算法能落地。真正让我记住这套方法价值的是一个分布式光伏接入方案的对比：两个候选接入节点，确定性潮流算出来电压都在合格范围，但概率潮流显示其中一个方案在午后时段电压越上限概率高达6%，另一个方案只有0.5%。如果没有随机潮流，这个问题在可研阶段很难被量化发现。

最后分享一个小技巧，做灵敏度矩阵传播时，不要一次性把所有节点的所有阶累积量都算出来，先只算电压幅值这个最重要的响应量，等验证通过后再扩展到相角和支路潮流。这个习惯能让你在调试阶段少掉一半的查错时间。半不变量法的核心优势是快，但它更依赖对工程模型的理解；把线性化条件、相关性处理、尾部修正这三个点想明白，在IEEE34上跑通之后，迁移到实际配电网模型只是数据替换的问题。