半不变量法概率潮流计算：IEEE34节点系统Matlab实现与验证

做潮流计算的朋友应该都遇到过这种尴尬：明明把负荷数据调得再准，算出来的节点电压也只是“某个工况下的值”。真正的电网随时都在波动，负荷在变，光伏在变，风电场出力也在变，单点算出来的潮流结果，很难回答一个最关键的问题：电压越限的概率到底有多大？

要回答这个问题，就得从确定性潮流走向概率潮流。这篇文章我用 Matlab 在 IEEE34 节点标准测试系统上实现了基于半不变量的概率潮流计算，通过 Gram-Charlier 级数展开还原节点电压的概率分布，并在同一测试系统下与蒙特卡洛模拟做了交叉验证。整个思路不难，但里面有不少容易踩坑的细节，我用实际踩过的坑把完整链路写出来，希望对刚接触概率潮流的朋友有帮助。

1. 为什么要用概率潮流：确定性潮流的盲区与半不变量法的切入点

1.1 一个节点电压到底是多少？这个问题本身就不完整

传统的潮流计算，比如最常用的 Newton-Raphson 法，输入是一组固定的负荷、固定的发电机出力，输出是一组确定的节点电压和支路功率。实际系统里根本没有“固定”这回事：用户负荷随早晚峰谷变化，屋顶光伏随云量抖动，电动汽车充电行为更是随机性极强。你把早上 10 点的负荷算得再准，也预测不了下午 3 点的电压。

这时候如果只看确定性潮流结果，往往会得到一种虚假的安全感。比如某节点电压算出来是 0.96 pu，距离下限 0.95 pu 还有一点裕度，你可能觉得没事。但负荷波动一大，电压低于 0.95 pu 的概率可能是 15%，也就是说每 7 次抽样就有一次越限。这个概率信息是确定性潮流给不了你的，而概率潮流给的就是这个：电压不是一个点，而是一条分布曲线。

1.2 概率潮流的几种主流路线，为什么单独挑半不变量法

概率潮流的实现路线大致有四类：

• 蒙特卡洛模拟（Monte Carlo Simulation）：最直观，随机生成大量负荷场景，逐个做潮流计算，最后统计电压样本的均值、标准差、越限概率。缺点是计算量大，尤其在大规模系统里跑一万次潮流很要命。
• 点估计法（Point Estimate Method）：用少量确定性的“代表点”代替整个分布，计算量比 MC 小很多，但只能给出前几阶矩，恢复完整分布时比较粗糙。
• 解析法（Analytical Method）：基于潮流方程线性化，直接把输入随机变量的统计特征传播到输出，速度快到可以忽略不计。半不变量法就是解析法里最主流的一个。
• 场景法 / 拟蒙特卡洛：介于 MC 与解析法之间，适合处理强非线性问题，但实现复杂度高。

半不变量法的核心优势是计算效率极高。整个流程只需要做一次确定性潮流得到线性化工作点，再做几阶矩阵运算，就能得到全系统所有节点电压的统计矩，再用 Gram-Charlier 或 Edgeworth 级数恢复概率密度。在 IEEE34 这种规模的系统上，跑一次解析法概率潮流的时间基本在 0.1 秒量级，而五千次蒙特卡洛模拟要几十秒。如果扩展到几百个节点的大规模输电网，差距会更明显。

1.3 半不变量到底是什么？一个类比帮你建立直觉

半不变量（Cumulant）是随机变量的一种数字特征，和矩（Moment）类似，但有一个矩不具备的优良性质：相互独立随机变量之和的半不变量可以直接相加。

这句话怎么理解？想象你手里有三堆豆子，每堆豆子的重量分布都不一样，你想知道三堆豆子倒在一起后的总重量分布。用矩做这件事，需要做卷积运算，很麻烦；用半不变量做，直接把每一阶半不变量加起来就行。

正是这个“线性可加性”，让半不变量特别适合处理潮流问题。因为潮流方程线性化之后，节点电压的随机偏差可以看成是各个节点注入功率随机偏差的加权和，而注入功率来自许多独立的随机源。有了半不变量，就能避开卷积，通过简单的代数运算完成随机变量的叠加。

2. 随机输入建模：负荷与光伏怎么变成统计参数

2.1 负荷波动：正态分布是最常用的假设

配电网负荷由大量独立小用户汇集而成，根据中心极限定理，节点净负荷通常可以用正态分布近似。在 IEEE34 节点系统中，每个负荷节点的有功期望值取基准潮流场景下的值，标准差取期望值的一定比例，比如 10% 到 15%。无功功率的处理我习惯于让它跟随有功变化，保持基准功率因数。

这样做的好处是每个节点只需要一个随机源“有功功率”，无功 Q 通过功率因数角换算成有功的线性函数，后面传播半不变量时只需构造一个等效灵敏度系数，实现起来简洁，物理上也符合实际：一个台区电压波动时，有功和无功往往是同涨同跌的。

如果确实需要让 P 和 Q 独立波动，也很简单，把 Q 当成第二个独立随机源，在传播公式里对 P 和 Q 分别算一遍半不变量相加即可。但要注意这样会丢失 P、Q 之间的相关性，实际负荷的有功和无功很少完全独立。

2.2 光伏出力：Beta 分布比正态分布更合理

光伏出力受辐照度影响，出力区间是 0 到额定容量之间，而且大概率不会满发也不会零出力，分布呈现明显的截断特性。这种情况下用正态分布会出现负出力的荒谬结果，所以通常用 Beta 分布拟合。

Beta 分布有两个形状参数 a 和 b，期望和方差分别满足：

• E[X] = a / (a + b)
• Var(X) = ab / ((a+b)^2 (a+b+1))

我在算例里把接入光伏的节点有功出力归一化到额定容量后，取 a = 4、b = 3，这样期望出力约为 0.57 倍额定容量，标准差约为 0.175 倍额定容量，比较接近实际光伏出力的波动形态。如果你手头有实测出力曲线，可以先用极大似然估计拟合出 a 和 b，再把 Beta 分布的矩转换成半不变量，流程是完全一样的。

2.3 从原始矩到半不变量：四阶以内就够用

得到随机变量的分布模型之后，下一步要把它转成半不变量。对于正态分布，前两阶半不变量就是均值和方差，偏度为 0，峰度修正也是 0。对 Beta 分布或其他分布，需要从中心矩换算。

设中心矩 m_r 为：

• m2 = E[(X - μ)^2]
• m3 = E[(X - μ)^3]
• m4 = E[(X - μ)^4]

则前四阶半不变量可以通过以下公式得到：

• κ1 = μ
• κ2 = m2
• κ3 = m3
• κ4 = m4 - 3 m2^2

如果只有离散样本，也可以先用历史数据计算中心矩，再套上面的公式。只要样本量足够大，半不变量的估计一致性是有保证的。我做这个项目时一开始直接把样本矩当半不变量用，结果四阶项明显偏大，后来才发现漏了 -3 m2^2 的修正项，画出来的概率密度尾部翘得离谱。

3. 半不变量法概率潮流推导：从卷积到 Gram-Charlier 级数

3.1 潮流方程线性化：把非线性问题变成线性映射

半不变量法的基础是潮流方程在期望工作点附近的线性化。设确定性潮流解为 V0 和 θ0，节点注入功率偏差为 ΔP、ΔQ，则节点电压幅值和相角的偏差可以写成：

[ \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix}

\begin{bmatrix} \frac{\partial P}{\partial \theta} & \frac{\partial P}{\partial V} \ \frac{\partial Q}{\partial \theta} & \frac{\partial Q}{\partial V} \end{bmatrix} \begin{bmatrix} \Delta \theta \ \Delta V \end{bmatrix} ]

这里的矩阵就是潮流计算的 Jacobian 矩阵 J。两边求逆，就能得到灵敏度矩阵 S = J^{-1}。也就是说，电压偏差是注入功率偏差的线性函数：

[ \Delta V_i = \sum_j S_{VP}(i,j) \Delta P_j + \sum_j S_{VQ}(i,j) \Delta Q_j ]

线性化看着简单，但有一个关键点容易被忽略：线性化点必须选在期望值所在的工况。如果在重载工况下做线性化，而期望负荷其实是轻载，算出来的电压分布会整体偏移。更稳妥的做法是先以各随机变量的期望值做一次确定性潮流，再在潮流解处取 Jacobian。

3.2 半不变量的线性传播：为什么系数要取 r 次方

这是整个方法最核心的一步。半不变量有一个重要性质：如果 Y = cX，其中 c 是常数，那么：

[ κ_r(Y) = c^r κ_r(X) ]

注意是 c 的 r 次方，不是 c 本身。这个性质直接决定了注入功率半不变量如何传播到电压：

[ κ_r(V_i) = \sum_j [S_{VP}(i,j)]^r κ_r(P_j) + \sum_j [S_{VQ}(i,j)]^r κ_r(Q_j) ]

如果 Q 通过基准功率因数跟随 P，即 ΔQ_j = k_j ΔP_j，可以定义等效灵敏度：

[ S_{eff}(i,j) = S_{VP}(i,j) + k_j S_{VQ}(i,j) ]

然后一步到位：

[ κ_r(V_i) = \sum_j [S_{eff}(i,j)]^r κ_r(P_j) ]

这样每个节点电压的四阶半不变量都能一次算出来。这里请一定小心矩阵元素的符号和数量级：灵敏度矩阵的元素可能很小，但取三次方、四次方之后数量级会剧烈变化，如果直接用全矩阵乘法而不是逐元素求幂，结果会完全错误。

3.3 Gram-Charlier 级数：从半不变量还原概率密度

半不变量只是数字特征，最终用户还是想看概率密度和越限概率。Gram-Charlier 级数做的事情，是用标准正态分布作为基底，加上 Hermite 多项式修正项去逼近真实分布。

标准化变量 z = (V - μ) / σ，其中 μ = κ1，σ = sqrt(κ2)。定义归一化半不变量：

• λ3 = κ3 / σ^3
• λ4 = κ4 / σ^4

则节点电压的概率密度函数可以展开为：

[ f(z) = \phi(z) \left[ 1 + \frac{λ3}{6} H_3(z) + \frac{λ4}{24} H_4(z) \right] ]

其中 φ(z) 是标准正态密度，Hermite 多项式前三项为：

• H3(z) = z^3 - 3z
• H4(z) = z^4 - 6z^2 + 3

有了概率密度，越限概率就通过数值积分计算。比如电压上限越限概率：

[ P(V > V_{max}) = \int_{z_{max}}^{\infty} f(z) dz ]

在 Matlab 里用 integral 函数直接算就行。不过有个经验要提前说：Gram-Charlier 级数展开到四阶时，在分布尾部偶尔会出现轻微负概率密度，这是级数截断的固有现象。遇到这种情况，优先检查输入半不变量的四阶项是否估计得过大，如果没问题，可以考虑取到六阶展开，甚至改用 Edgeworth 展开。

3.4 完整计算流程：一张图看懂主框架

整个计算流程可以归纳成六步，我实际写代码时也是按这个顺序拆函数的：

1. 从随机输入模型计算每个节点注入功率的前四阶半不变量。
2. 以各随机变量的期望值做一次确定性潮流，得到期望工作点。
3. 在期望工作点提取潮流 Jacobian，求逆得到灵敏度矩阵。
4. 按半不变量线性传播公式，把注入功率半不变量映射为节点电压半不变量。
5. 用 Gram-Charlier 级数把电压半不变量还原为概率密度函数。
6. 对每个节点做数值积分，得到电压越限概率以及均值、标准差等统计量。

4. IEEE34 节点系统搭建与 Matlab 实现细节

4.1 IEEE34 测试系统的基本情况

IEEE34 节点系统是一个经典的配电网测试馈线，额定电压 24.9 kV，包含 34 个节点、33 条支路，既有三相段也有单相段，负荷类型包括恒功率、恒阻抗、恒电流等多种。这个系统常用来验证配电网潮流、重构、分布式电源接入等算法。

用它来做概率潮流有个好处：规模不大不小，既能体现概率算法相对蒙特卡洛的效率优势，又不会让调试复杂度失控。我在实现时做了简化处理，先按三相平衡等效模型计算，忽略线路不平衡带来的中性点电压偏移。如果你要做更精细的分析，需要换成三相四线模型，但半不变量法的整体框架是不变的。

4.2 代码结构设计

整个项目的 Matlab 代码我拆成了五个文件，职责非常清晰：

• load_system_data.m：读取 IEEE34 节点系统参数，构造导纳矩阵和初始潮流数据。
• deterministic_pf.m：实现 Newton-Raphson 确定性潮流，返回电压解和最后一次迭代的 Jacobian 矩阵。
• compute_input_cumulants.m：根据设置的随机输入模型（正态负荷、Beta 光伏）计算各节点注入功率的前四阶半不变量。
• propagate_cumulants.m：传入灵敏度矩阵和输入半不变量，输出每个节点电压的四阶半不变量。
• gram_charlier_pdf.m：根据电压半不变量还原概率密度，并计算越限概率。

这样的拆分方式让我在调试时能单独验证每一步。比如先跑一次确定性潮流，核对结果是不是与已知潮流解一致，再单独检查灵敏度矩阵的元素，看符号有没有反。

4.3 核心函数实现：灵敏度传播与 Gram-Charlier 恢复

下面这段代码是半不变量传播的核心逻辑。输入是灵敏度矩阵 SVP、SVQ，功率因数斜率数组 k，以及每个节点的前四阶半不变量矩阵 cumulants_in（每一列对应一阶，每一行对应一个节点）。

% k: 每个节点的无功/有功比率 k_j % cumulants_in: nNode x 4 矩阵，依次为 kappa1..kappa4 % cumulants_v: nNode x 4 矩阵，输出节点电压幅值的前四阶半不变量 function cumulants_v = propagate_cumulants(SVP, SVQ, k, cumulants_in) nNode = size(SVP, 1); nSource = size(SVP, 2); % 等效灵敏度：考虑 Q 跟随 P Seff = SVP + SVQ .* repmat(k(:)', nNode, 1); cumulants_v = zeros(nNode, 4); % 一阶：均值直接通过灵敏度传播 cumulants_v(:, 1) = Seff * cumulants_in(:, 1); % 二阶到四阶：系数取 r 次方 for r = 2:4 Seff_r = abs(Seff).^r; cumulants_v(:, r) = Seff_r * cumulants_in(:, r); end end

注意这里我对 Seff 取了绝对值再求 r 次方。因为灵敏度矩阵里某些元素可能是负数，但半不变量的传播系数是常数的 r 次方，负数的偶数次方为正、奇数次方保持负号，而绝对值写法只适用于偶数阶。严格来说应该写成 sign(Seff).^r .* abs(Seff).^r，或者直接 Seff.^r。上面代码里用 abs 处理有简化嫌疑，更严谨的写法是：

Seff_r = Seff.^r;

当时我为了保险起见先取绝对值，是因为电压幅值对注入功率的灵敏度大多是正的，负元素数量不多，但严谨起见建议直接用 .^r，代数上才是对的。

Gram-Charlier 概率密度恢复的代码如下：

% mu, sigma: 电压估值的期望和标准差 % k3, k4: 电压的三阶、四阶半不变量 % vline: 待求概率密度的电压点向量 function p = gram_charlier_pdf(vline, mu, sigma, k3, k4) z = (vline - mu) / sigma; lambda3 = k3 / sigma^3; lambda4 = k4 / sigma^4; H3 = z.^3 - 3*z; H4 = z.^4 - 6*z.^2 + 3; phi = (1 / sqrt(2*pi)) * exp(-0.5 * z.^2); p = phi .* (1 + lambda3/6 .* H3 + lambda4/24 .* H4); end

算越限概率时就对这个密度从 1.05 pu 或 0.95 pu 积分：

p_upper = integral(@(v) gram_charlier_pdf(v, mu, sigma, k3, k4), 1.05, Inf); p_lower = integral(@(v) gram_charlier_pdf(v, mu, sigma, k3, k4), -Inf, 0.95);

5. 算例结果解读：电压分布、越限概率与蒙特卡洛对比

5.1 算例设置

我在 IEEE34 节点系统的负荷节点上设置了两种随机场景：

• 场景 A：所有节点有功负荷服从正态分布，标准差取期望值的 10%，无功按基准功率因数跟随有功。
• 场景 B：在节点 822 和节点 858 分别接入额定容量 0.5 MW 和 0.8 MW 的光伏，光伏出力用 Beta(a=4, b=3) 分布模拟，负荷波动参数与场景 A 一致。

蒙特卡洛对比组的样本数为 5000，每次抽样都用完整的 Newton-Raphson 潮流求解。

5.2 电压期望和标准差的对比

以场景 B 中两个光伏接入点附近的电压为例，半不变量法与蒙特卡洛的结果对比如下：

从这个表能看出两件事。第一，解析法得到的电压期望与蒙特卡洛结果非常接近，相对偏差基本在 0.1% 以内，说明线性化点选在期望值处是合理的。第二，标准差越接近电压下限的节点，偏差越大，节点 858 的标准差相对偏差到了 5% 左右，这是因为该节点距离电源点远、灵敏度系数大，高阶半不变量在传播过程中被放大。

另外要注意，确定性潮流得到的 0.9497 pu 已经低于 0.95 pu，但这个值只是一个工况点。如果光伏大发，该节点电压会有所抬升，所以实际越限情况并没有确定性结果显示的那么悲观。

5.3 概率密度和越限概率的对比

节点 858 的电压概率密度曲线，用 Gram-Charlier 级数恢复出来的形状和蒙特卡洛直方图对比如下：均值附近拟合得非常好，但左侧尾部比蒙特卡洛直方图略薄一些。这直接反映在越限概率上：

越限概率对分布尾部非常敏感，所以解析法和蒙特卡洛的误差会比均值和标准差大一些。6.8% 对 5.9%，绝对偏差不到一个百分点，对于规划阶段的评估来说是可以接受的。如果你的场景里负荷波动特别大，或者有强非线性设备，建议把 Gram-Charlier 展开到第六阶甚至第八阶，能显著改善尾部精度。

5.4 计算效率对比

这是半不变量法最让人满意的地方。在我自己的电脑上，场景 B 下：

• 半不变量法整个过程（一次确定性潮流 + 半不变量传播 + 概率密度恢复）耗时约 0.11 秒。
• 蒙特卡洛模拟 5000 次完整 Newton-Raphson 潮流耗时约 27 秒。
• 点估计法如果能收敛，大约在 0.5 秒左右。

三千多倍的效率差异，放在实时调度或在线风险评估场景里就是能不能用的区别。当然，解析法的精度受线性化假设约束，并不适合在所有场景下替代蒙特卡洛，但它非常适合做“快速筛查”：先用半不变量法把所有节点扫一遍，找出越限风险高的薄弱节点，再对这些节点单独跑更精确的蒙特卡洛验证。

6. 踩坑与调参经验：比理论推导更重要的实战细节

6.1 四阶半不变量的负值陷阱

算节点电压的四阶半不变量时，即使输入 Beta 分布的四阶半不变量是正的，经过灵敏度矩阵传播后，某些节点的四阶半不变量可能变成负值，甚至数值大到不合理。这通常不是程序算错了，而是输入模型的四阶中心矩与正态分布偏离太大，导致四阶半不变量在数值上不稳定。

我调试时习惯先把所有节点的 κ4 画出来看一遍，如果某个节点数值比 κ2^2 还大几个量级，就要怀疑输入数据的问题。也可以用公式 κ4 = m4 - 3m2^2 反推，如果 m4 和 3m2^2 几乎相等，说明输入分布接近正态，四阶修正项很小，这是最理想的情况。

6.2 Gram-Charlier 尾部的负概率密度怎么处理

Gram-Charlier 级数展开不是万能的，在偏度或峰度偏离正态比较大的节点上，概率密度尾部可能出现负值，导致越限概率积分出现负数。我的处理办法是：先检查 λ4 的绝对值是否大于 3，如果 λ4 超过约 4，四阶展开基本不可靠，需要增加展开阶数或改用 Edgeworth 展开。

另一个实用技巧是：解析法得到概率密度后，在计算越限积分前，把负密度的部分用周围正密度的插值顶掉，或者直接把负值截断后重新归一化。这个做法不严谨，但工程上能避免出现“越限概率为负”这种尴尬结果，我通常只在快速排查阶段用，正式报告里还是以蒙特卡洛复核为准。

6.3 线性化误差的来源与控制

半不变量法最大的理论短板是潮流方程线性化。负荷波动越大，线性化误差越大。我的经验界限是：节点注入功率的标准差不超过期望值的 20% 时，解析法的均值和标准差精度都非常好，越限概率也有参考价值；超过 30% 以后，建议把潮流方程在期望点处重新线性化，做一两次迭代修正，或者直接用点估计法兜底。

另外一点，IEEE34 节点系统里有不少恒阻抗和恒电流负荷成分，这些负荷的功率在电压变化时本身会跟着变。简单地把所有负荷都当成恒功率正态随机源，会低估节点电压的波动范围。我的做法是：对于恒阻抗负荷，先把阻抗按期望电压折算成等效功率，再做随机建模，这样线性化点更贴近真实工况。

6.4 一点个人体会

做这个项目之前，我一直觉得概率潮流是高门槛的学术概念，真上手后发现核心思想并不复杂：用半不变量把随机变量之间的卷积变成加法，用灵敏度矩阵把注入功率的随机性传播到电压，再用 Gram-Charlier 把统计特征还原成直观的分布曲线。难点恰恰在那些理论推导里不会写的细节，比如灵敏度矩阵元素的次方处理、四阶矩修正项、分布尾部的数值稳定性。

如果让我给刚接触这个方向的朋友一个建议，我会说：先别急着在 IEEE34 这种标准系统上跑完整代码，拿一个三节点手算系统，把每一步的矩阵、半不变量、概率密度都手推一遍，再用 Matlab 核对。把这个小系统彻底搞通之后，再上 IEEE34 节点系统，你会有一种“原来如此”的通透感。概率潮流这东西，理论半小时，调参三天，但一旦把链路跑通，后面扩展到风机、储能、电动汽车充电负荷，都只是换输入分布模型的事。