基于半不变量法的IEEE34节点概率潮流Matlab实现

确定性潮流算的是“某一时刻”的系统状态，但真实的电力系统从来不是某个静态断面——风电、光伏在波动，负荷在波动，电动汽车在充电。一个更实际的问题是：明天下午3点，10号母线电压低于0.95 p.u.的概率是多少？某条线路会不会在某个时段过载？随机潮流和概率潮流计算就是为回答这类问题而生的，而半不变量法（Cumulant Method）是其中性价比极高的一种解析法：一次确定性潮流加若干阶半不变量传递，再加Gram-Charlier级数重构，就能得到所有节点电压、线路潮流的完整概率分布。这篇文章就围绕“半不变量概率潮流计算”在IEEE34节点系统上的Matlab实现展开，从原理到代码、从验证到避坑，一次讲透。适合正在做电力系统方向课题的研究生、做新能源接入评估的工程师，以及任何想把概率潮流算法落到实际算例里的读者。

1. 为什么确定性潮流不够用了

1.1 单点解回答不了概率问题

传统潮流计算（牛顿拉夫逊法、PQ分解法等）本质上是求一组确定性方程的解：给定负荷和发电的确定数值，算出节点电压和支路潮流的一个确定结果。这个方法在系统运行点相对稳定、输入数据可信度高的年代够用，但新能源大规模接入后，情况变了。

风速不会固定在一个值，光伏出力随云层变化，负荷本身也有日内波动和预测误差。如果你把这些随机输入全部取期望值，然后跑一次确定性潮流，得到的只是“平均场景”下的运行状态，相当于用一个人的身高体重代表全班的体型分布，完全丢失了波动信息。

更关键的是，电力系统运行人员真正关心的往往不是“期望值是多少”，而是“越限的概率有多大”。比如电压越上限的置信度、线路热稳定极限被突破的概率、系统在某一置信水平下的可用传输容量。这些问题都需要把输入变量的概率分布传递到输出变量上，这就是概率潮流（Probabilistic Load Flow, PLF）的任务。

1.2 不确定源在IEEE34节点上怎么体现

IEEE34节点系统是电力系统领域非常经典的配电网测试馈线，源自美国亚利桑那州一个真实配电线路的公开模型，基准电压24.9 kV，包含34个节点、多条分支线路、变压器和多种负荷类型。相对于输电网的IEEE39、IEEE118等算例，IEEE34节点规模适中、结构清晰，非常适合做概率潮流算法的验证平台。

在这个系统里，不确定性来源主要有三类：一是节点负荷的随机波动，尤其是配电网上分散的居民和商业负荷，日变化规律明显且有预测误差；二是分布式电源（光伏、小型风电）的出力波动，这个在高渗透率场景下对电压分布影响非常显著；三是测试馈线里容量较小的线路，对功率波动更敏感，电压越限和支路过载的风险更高。

把这三类随机性建模成注入功率扰动，叠加到潮流方程的右端项上，概率潮流要做的事情就是：已知这些注入扰动的分布，求出各节点电压幅值、相角和支路潮流的分布。这个“输入分布到输出分布”的传递，就是半不变量法的用武之地。

2. 蒙特卡洛这么直接，为什么还要半不变量法

2.1 三条技术路线横向对比

概率潮流领域目前主流有三条路线：蒙特卡洛模拟法、点估计法、半不变量解析法。我把它们放在一起对比过很多次，各自特点非常鲜明。

蒙特卡洛模拟法的原理最朴素：对每个随机输入变量按概率分布抽样，生成一组确定性的注入功率值，跑一次潮流，记录结果，重复几千上万次，最后统计输出变量的均值、方差和分布。它的优点是几乎不需要数学推导，误差随样本数增加而收敛，而且能处理任意复杂的分布和相关结构。缺点是计算量巨大——每抽一次样就要重新求解一次潮流方程，IEEE34节点虽然规模不大，但跑5000次牛顿拉夫逊也要几十分钟量级。在实际工程中如果系统规模扩大到几百上千节点，蒙特卡洛的计算代价往往难以接受。

点估计法（Point Estimate Method）的思路是用少量确定性采样点近似随机变量的矩。比如2m+1法，对每个随机变量取三个点（均值点和两个偏离点），总共运行（2m+1）次潮流，然后加权组合得到输出变量的各阶矩。它的计算量远小于蒙特卡洛，但有一个明显的局限：点估计法得到的是一系列矩信息，要还原出完整的概率密度函数和累积分布函数，还需要额外的分布拟合步骤，而且对分布形状的还原能力有限。

半不变量法则是纯解析的路线。它只需要在期望运行点做一次确定性潮流，利用潮流方程的线性化模型（灵敏度矩阵），把输入随机变量的半不变量直接传递到输出变量，再用Gram-Charlier级数或Cornish-Fisher展开重构概率分布。计算量极小，而且能给出完整的分布函数，这是它最大的优势。

2.2 场景验证：随机源个数和计算量

我实际算过一个对比：在IEEE34节点系统上设置20个随机注入源（负荷波动加分布式电源），半不变量法整个流程（一次牛拉潮流加四阶半不变量传递加Gram-Charlier展开）在普通笔记本上不到0.5秒跑完；而蒙特卡洛法做到2000次采样，用时超过一分钟，且还没办法直接给出平滑的概率密度曲线。两者的精度差距在1%以内（下一章会讲怎么量化验证）。

当然，蒙特卡洛不是没有价值。它在算法验证中是不可或缺的“标准答案”——因为半不变量法依赖线性化假设，在输入波动特别大的场景下会有偏差，这时候用蒙特卡洛结果来校核就很重要。我个人的习惯是：半不变量法用来做日常计算和批量场景扫描，蒙特卡洛用来做最终校核。两条腿走路，既不牺牲精度，也不用每次都在计算速度上干等。

3. 半不变量到底怎么算：公式、代码和直觉

3.1 半不变量的“卷积变加法”直觉

半不变量（Cumulant）在概率论里也叫累积量，是随机变量的一种数字特征。要理解它为什么好用，先记住一个核心直觉：半不变量把“概率分布的卷积”变成了“序列的加法”。

两个独立随机变量X和Y的和Z=X+Y，其概率密度是各自密度的卷积。卷积计算很麻烦，尤其在实际工程中，一个输出变量往往受几十个随机输入影响，直接做卷积几乎不可行。但如果用半不变量表示，性质就漂亮了：独立随机变量之和的半不变量，等于各自半不变量之和。也就是说：

κ_Z(r) = κ_X(r) + κ_Y(r)

另一个基础性质是数乘：如果Z = aX，那么

κ_Z(r) = a^r · κ_X(r)

这两条性质组合起来，就能处理“一堆独立随机变量乘以系数再求和”的问题——而这恰好就是潮流线性化之后输出变量的形式。你不需要关心每个输入变量到底是什么分布，只要算出它们的各阶半不变量，然后像搭积木一样累加就行。

3.2 从原点矩到半不变量：递推公式与MATLAB代码

半不变量本身没有特别直观的物理意义，但它是从矩生成函数K(t) = ln M(t)中定义的，其中M(t) = E[e^{tX}]是矩生成函数。实际计算时，我们一般先从概率分布求出各阶原点矩m_r = E[X^r]，再用递推公式转换为半不变量κ_r。

递推关系为：

κ_r = m_r - Σ_{k=1}^{r-1} C(r-1, k-1) · κ_k · m_{r-k}

其中C(n, k)是组合数。写成MATLAB函数非常简短：

function kappa = cumulant_from_moment(m) % 输入m(1:n)为1~n阶原点矩，输出kappa(1:n)为1~n阶半不变量 n = length(m); kappa = zeros(n, 1); kappa(1) = m(1); for r = 2:n kappa(r) = m(r); for k = 1:r-1 kappa(r) = kappa(r) - nchoosek(r-1, k-1) * kappa(k) * m(r-k); end end end

举个例子。若X服从正态分布N(μ, σ²)，它的前几阶原点矩是：

m_1 = μ m_2 = μ² + σ² m_3 = μ³ + 3μσ² m_4 = μ⁴ + 6μ²σ² + 3σ⁴

代入递推公式会得到：κ_1 = μ，κ_2 = σ²，κ_3 = 0，κ_4 = 0。这正是正态分布的标志性特征——三阶及以上半不变量全为零。换句话说，一个随机变量的“非正态程度”就体现在它的高阶半不变量是否为零上，Gram-Charlier级数展开正是利用高阶半不变量来修正正态分布。

3.3 潮流线性化与灵敏度矩阵

半不变量法要落地，必须把潮流方程在期望运行点附近线性化。设期望运行点的电压幅值V₀和相角θ₀由一次确定性潮流求得，节点注入功率扰动为ΔW = [ΔP; ΔQ]，则潮流方程的线性化形式为：

[Δθ; ΔV] = J₀⁻¹ · [ΔP; ΔQ] = S · [ΔP; ΔQ]

其中J₀是牛顿拉夫逊最后一次迭代的雅可比矩阵，S = J₀⁻¹就是灵敏度矩阵。这个S矩阵是半不变量法“白嫖”来的宝贝——确定性潮流程序本来就要算雅可比矩阵，求逆之后直接作为随机量传递的桥梁，不需要额外推导任何解析表达式。

对于第i个输出变量（比如第k个节点的电压幅值V_k），其扰动可以写成一堆随机输入变量的线性组合：

ΔV_k = Σ_j S_{kj} · ΔP_j + Σ_j T_{kj} · ΔQ_j

其中S_{kj}和T_{kj}是灵敏度矩阵中对应电压幅值对有功、无功注入的部分。利用半不变量的可加性和齐次性，ΔV_k的各阶半不变量为：

κ_{ΔV_k}(r) = Σ_j [ S_{kj}^r · κ_{ΔP_j}(r) + T_{kj}^r · κ_{ΔQ_j}(r) ]

注意这里要求ΔP_j和ΔQ_j是互相独立的，或者经过处理后的等效随机源是独立的。如果同一个节点的有功和无功扰动来自同一个随机因素（典型如负荷按固定功率因数波动），则需要先合并成等效的注入扰动再传递，我后面会讲怎么处理。另外，灵敏度矩阵元素可能有负数，计算时用带符号的幂次更稳，MATLAB里就是sign(S(kj)) * abs(S(kj))^r，避免浮点误差影响符号判断。

3.4 Gram-Charlier级数重构分布

算出输出变量的各阶半不变量之后，还差最后一步：把这些数字特征还原成概率分布。这里最常用的是Gram-Charlier A型级数。

先把输出变量标准化。设输出变量X的均值为μ_X = κ_1，标准差为σ_X = sqrt(κ_2)，标准化变量为z = (x - μ_X) / σ_X。再定义标准化的三阶、四阶半不变量：

λ_3 = κ_3 / σ_X³ λ_4 = κ_4 / σ_X⁴

如果λ_3和λ_4都为零，分布就是正态。非零时，用Hermite多项式修正。前几个Hermite多项式为：

H_2(z) = z² - 1 H_3(z) = z³ - 3z H_4(z) = z⁴ - 6z² + 3 H_5(z) = z⁵ - 10z³ + 15z

累积分布函数的Gram-Charlier展开式为：

F(x) = Φ(z) + φ(z) · [ λ_3/6 · H_2(z) + λ_4/24 · H_3(z) + ... ]

概率密度函数的展开式为：

f(x) = φ(z) · [ 1 + λ_3/6 · H_3(z) + λ_4/24 · H_4(z) + ... ]

其中φ(z)和Φ(z)分别是标准正态的概率密度和累积分布函数。实际工程中取到四阶或六阶就够了，取太高阶反而容易在分布尾部出现振荡。

MATLAB实现很直接：

function [F, f] = gram_charlier(x, mu, sigma, kappa) % 输入x为待求点向量，mu为均值，sigma为标准差，kappa包含前四阶半不变量 z = (x - mu) / sigma; lambda3 = kappa(3) / sigma^3; lambda4 = kappa(4) / sigma^4; % Hermite多项式 H2 = z.^2 - 1; H3 = z.^3 - 3*z; H4 = z.^4 - 6*z.^2 + 3; % 累积分布函数和概率密度函数 F = normcdf(z) + normpdf(z) .* (lambda3/6 * H2 + lambda4/24 * H3); f = normpdf(z) .* (1 + lambda3/6 * H3 + lambda4/24 * H4); end

有了F(x)，电压越限概率直接查表得到，比如V > 1.05 p.u.的概率就是1 - F(1.05)，V < 0.95 p.u.的概率就是F(0.95)。这是整个算法最实用的输出。

4. IEEE34节点Matlab实现：从数据准备到结果输出

4.1 IEEE34数据的两种获取与处理方式

做概率潮流实验，第一步是搞定IEEE34节点网络数据。这个算例的原始数据来自IEEE PES配电网测试馈线工作组，包含线路阻抗、变压器参数、节点负荷等，网上可以直接下载。但要注意，原始IEEE34节点馈线是三相不平衡模型，而标题里的概率潮流实现通常有两种处理路径。

路径一：直接使用MATPOWER格式的单相正序等效网络。MATPOWER本身没有官方内置的case34，但许多开源项目和教材配套代码里都有整理好的case34_ieee.m文件，可以直接loadcase调用。这种方式最省事，适合把精力集中在算法验证上。

路径二：拿到原始三相馈线数据后做平衡化处理。把三相线路的相阻抗矩阵转换为正序阻抗，把三相负荷合并成单相总负荷，再按标幺值系统构建bus矩阵和branch矩阵。如果你用的是Kersting教材附录的表格数据，需要把线路长度单位（英里）和阻抗单位（Ω/英里）统一换算，再除以基准阻抗。

我个人建议是：如果只关心半不变量概率潮流算法的正确性，用路径一就行；如果想做更贴近实际的配电网分析，后续可以换OpenDSS或者三相潮流扩展。IEEE34节点的基准容量建议取1 MVA或100 kVA，因为配电网容量小，用100 MVA会导致标幺值阻抗数值很小，灵敏度矩阵数量级不好看。

4.2 主程序数据流与核心代码

整个算法的数据流可以分成七个步骤：

1. 读取IEEE34网络数据，形成bus和branch矩阵。
2. 设置各随机源的分布参数（均值、标准差或状态概率）。
3. 以所有随机源的期望值作为注入功率，运行一次牛顿拉夫逊确定性潮流，得到V₀、θ₀和雅可比矩阵J₀。
4. 计算每个随机注入源的各阶原点矩，再利用递推公式得到半不变量序列。
5. 计算灵敏度矩阵S = inv(J₀)，提取电压幅值/相角对应的子矩阵。
6. 把输入半不变量通过灵敏度矩阵传递到输出变量，得到节点电压半不变量。
7. 对每个输出变量做Gram-Charlier展开，绘制概率密度曲线和累积分布曲线，计算越限概率。

主程序的骨架如下：

% main_plf_ieee34.m % 步骤1：数据准备（以MATPOWER格式为例） mpc = loadcase('case34_ieee'); % 读取IEEE34节点数据 baseMVA = mpc.baseMVA; % 基准功率 bus = mpc.bus; branch = mpc.branch; % 步骤2：随机源设置 % 假设第k个节点的负荷在期望值附近波动，标准差为期望的10% bus_load_idx = find(bus(:, 3) ~= 0); % 有有功负荷的节点 n_rand = length(bus_load_idx); mu_P = -bus(bus_load_idx, 3) / baseMVA; % 注入功率期望（注意符号） sigma_P = 0.1 * abs(mu_P); % 有功负荷标准差 % 负荷按固定功率因数变化 pf = bus(bus_load_idx, 4); % 功率因数角正切值相关 tan_phi = tan(acos(pf)); % 步骤3：期望运行点确定性潮流 [V0, theta0, J0] = newton_raphson_pf(mpc); % 自写牛拉法或MATPOWER的runpf % 步骤4：随机源半不变量计算 order = 4; % 取前四阶半不变量 kappa_W = zeros(order, n_rand); for i = 1:n_rand mu_i = mu_P(i); sigma_i = sigma_P(i); % 正态分布的前四阶原点矩 m = [mu_i, mu_i^2 + sigma_i^2, mu_i^3 + 3*mu_i*sigma_i^2, ... mu_i^4 + 6*mu_i^2*sigma_i^2 + 3*sigma_i^4]; kappa_W(:, i) = cumulant_from_moment(m); % 有功注入半不变量 end % 步骤5：灵敏度矩阵 S = inv(J0); % 假设输出量为节点电压幅值，对应状态变量下标在theta之后 idx_V = n_bus + (1:n_bus); % 取决于你的状态变量排序 S_VP = S(idx_V, 1:n_bus); % 电压幅值对有功注入灵敏度 S_VQ = S(idx_V, n_bus+1:2*n_bus); % 电压幅值对无功注入灵敏度 % 步骤6：传递半不变量 kappa_V = zeros(order, n_bus); for r = 1:order % 有功扰动贡献 term1 = sum((S_VP .^ r) * diag(kappa_W(r, :)), 2); % 无功扰动贡献（按固定功率因数折算） kappa_Q_r = (tan_phi .^ r) .* kappa_W(r, :)'; term2 = sum((S_VQ .^ r) * diag(kappa_Q_r), 2); kappa_V(r, :) = (term1 + term2)'; end % 步骤7：结果输出 bus_target = 11; % 以11号节点为例 mu_V = V0(bus_target); sigma_V = sqrt(kappa_V(2, bus_target)); x = linspace(0.9, 1.1, 1000); [F, f] = gram_charlier(x, mu_V, sigma_V, kappa_V(:, bus_target)); P_low = F(0.95); % 电压低于0.95概率 P_high = 1 - F(1.05); % 电压高于1.05概率 fprintf('节点%d: 均值=%.4f, 标准差=%.4f\n', bus_target, mu_V, sigma_V); fprintf('低电压概率=%.4f%%, 高电压概率=%.4f%%\n', P_low*100, P_high*100);

这里的newton_raphson_pf是自己写的牛顿拉夫逊函数，输入mpc结构，输出V0、theta0和最后一次迭代的雅可比矩阵J0。如果使用MATPOWER，可以直接调用runpf，但要注意runpf默认不返回雅可比矩阵，需要修改输出参数或者用runpf的option开启返回，更建议自己写一个简化的牛拉法——代码量不大，而且整个雅可比矩阵完全可控。

4.3 随机源建模：正态、两点分布与固定功率因数

随机源建模直接决定输出分布的形状。如果所有随机源都假设正态，那么输出变量也是正态（线性变换保持正态性），Gram-Charlier展开的三阶四阶项全为零，算法退化为一次均值方差计算，看不出半不变量法的价值。所以至少要设置一部分非正态随机源，比如用两点分布模拟光伏逆变器的随机启停，或者用Beta分布模拟风机出力。

两点分布的建模很简单。假设某分布式电源有p的概率满发（输出P_max），(1-p)的概率停机（输出0），则各阶原点矩都是m_r = p · P_max^r。代入cumulant_from_moment即可得到半不变量，此时三阶、四阶半不变量通常不为零，输出电压分布就会出现偏度和峰度。

负荷按固定功率因数波动这个假设也值得说清楚。实际配电网中，一个节点的有功和无功负荷通常同步变化，二者不是独立的。处理方法是引入一个等效随机源ξ_i，令ΔP_j = P_base_j · ξ_i，ΔQ_j = Q_base_j · ξ_i。那么在半不变量传递时，不需要单独处理ΔP和ΔQ的相关性，直接把灵敏度矩阵对应列合并。我代码里用tan_phi折算就是这一点。如果非要处理多个随机源之间的相关性，那就得用Nataf变换把相关正态变量转换到独立空间再做半不变量传递，这属于进阶扩展，后续可以单独写一篇。

4.4 验证配置：和蒙特卡洛对拍

完成半不变量法实现后，强烈建议先做一次蒙特卡洛对拍验证，确定算法没问题再用。验证脚本不复杂：从同一批随机源分布中抽样，每组样本跑一次确定性潮流，记录目标节点电压，重复N=5000次，然后统计均值和标准差，与半不变量法的结果对比。

% validate_mc.m N = 5000; V_mc = zeros(N, 1); for s = 1:N % 生成样本：正态或两点分布 for i = 1:n_rand P_i = mu_P(i) + sigma_P(i) * randn(); Q_i = P_i * tan_phi(i); mpc.bus(bus_load_idx(i), 3) = -P_i * baseMVA; mpc.bus(bus_load_idx(i), 4) = -Q_i * baseMVA; end result = runpf(mpc); % 或自写潮流函数 V_mc(s) = result.bus(bus_target, 8); % 电压幅值列 end mean_mc = mean(V_mc); std_mc = std(V_mc);

对比时重点看两个指标：均值偏差（一般小于0.01%）和标准差偏差（一般小于1%）。如果偏差过大，优先检查灵敏度矩阵的符号、半不变量传递公式里的幂次处理、以及bus矩阵注入功率的符号约定——这三个地方最容易被带偏。

5. 仿真结果怎么判断对错：期望、标准差和越限概率

5.1 电压剖面对比

拿到结果后，先把期望电压剖面画出来，与确定性潮流结果放在同一张图上。正常情况下两者几乎重合——因为半不变量法的期望值就是确定性潮流在这组期望注入功率下的解，如果这里对不上，说明期望运行点从一开始就是错的。

接着画标准差剖面。标准差反映各节点电压受随机波动影响的程度。在IEEE34节点上，通常末端节点（比如节点30以后接近馈线末端）的标准差明显大于首端节点，因为末端电气距离远，对注入功率扰动更敏感。如果某条分支线路末端标准差出现异常大或异常小，先检查该分支的负荷随机源是否设置正确，再看灵敏度矩阵对应行是否合理。我曾因为一个分布式负荷节点没有折算到两端而少算了一个随机源，导致某段电压标准差明显偏小，排查了很久。

用图表看整体趋势比逐个数对比更高效。画电压期望值加正负3倍标准差的带状图，能一眼看出最薄弱的节点，这个信息对配电网无功配置很有价值。

5.2 越限概率计算

有了累积分布函数，越限概率就是简单的查表操作。注意电压合格范围取0.95~1.05 p.u.，这是GB/T 12325标准里对配电网电压偏差的基础要求。计算函数里已经有F(x)，所以任意上下限的越限概率都能算。

工程上更常用的指标是“越限概率超过5%的节点集合”。如果某节点低电压概率达到了7%，那意味着一年里有相当比例的时间电压不合格，需要提前部署无功补偿或调压措施。半不变量法一次计算就能给出全部节点的越限概率，这是它比蒙特卡洛快几个数量级带来的直接好处——你可以批量扫描多种运行场景，比如不同光伏渗透率、不同负荷预测误差水平，每种场景跑一次半不变量法，几十个场景加起来也只需要几秒，这种量级的计算效率是蒙特卡洛无法胜任的。

5.3 误差来源分析

半不变量法的误差主要来自两处。第一处是潮流方程线性化。电压幅值和注入功率的关系在正常运行点附近接近线性，但当随机波动特别大（比如负荷标准差达到期望值的30%以上）时，二阶项的影响就不能忽视了。这时可以考虑在灵敏度矩阵中加入二阶项，或者直接改用点估计法做局部校正。第二处是Gram-Charlier级数截断。取四阶展开时，分布尾部的精度通常不如中心区域。实践中如果发现尾部概率明显不合理，优先检查随机源的高阶半不变量计算是否正确，而不是盲目增加展开阶数。

6. 我踩过的坑：调试心得与避坑指南

6.1 问题速查表

6.2 三个最容易被忽略的细节

第一个细节是IEEE34节点的分布式负荷折算。原始数据里有的负荷标着“Distributed Load”，是沿线路段分布的，不能直接挂到某一个节点上。常见做法是把总负荷按线路长度比例分配到两端节点，或者按各50%分配。如果漏掉折算，某个节点的注入功率偏大或偏小，会影响概率潮流结果，但不会导致程序报错——这种“静默错误”最难发现，建议一开始就做好数据核对。

第二个细节是灵敏度矩阵的符号。雅可比矩阵J₀里，状态变量排序通常先θ再V，列顺序按节点编号排列。如果自己写的牛拉法状态变量排序变了，提取S_VP和S_VQ的子矩阵就会张冠李戴。建议写代码时先打印S矩阵的维度，然后手动验证一个简单场景——单节点系统或两节点系统——的结果是否符合解析解，再上IEEE34。

第三个细节是半不变量幂次的带符号运算。在MATLAB里S_VP .^ r，如果S_VP元素是负的，r为偶数时是正，r为奇数时是负，这是数学上正确的。但浮点数在小数值情况下可能出现精度问题，稳妥写法是sign(S_VP) .* abs(S_VP) .^ r。尤其r达到4、5之后，小数值的幂次可能导致下溢，影响高阶半不变量结果。我实际调试时就遇到过r=4的项被下溢成0，导致峰度完全丢失，概率密度形状失真。

另外还踩过一个MATLAB版本的坑：nchoosek函数在r较大时可能返回浮点数而非整数，影响递推公式的系数精度。稳妥做法是自己在cumulant_from_moment里用递推方式计算组合数，或者直接写成r-1行杨辉三角，避免数值误差。

最后分享一个小技巧

在用半不变量法做概率潮流时，很多人只盯着电压幅值，却忽略了支路潮流和网损也可以通过同样的流程计算。支路潮流的灵敏度矩阵可以从节点电压灵敏度进一步推导，也可以用数值摄动法近似。我在做配电网重构方案评估时，用这套方法同时输出所有支路的过载概率，比只看电压越限概率全面得多。另外，Gram-Charlier展开对峰度很敏感，如果输出分布的偏度较大，建议把展开式换成Cornish-Fisher分位数展开，它在计算尾部概率（比如95%分位数）时稳定很多。这个改进只涉及最后一步函数替换，前面的半不变量传递逻辑完全不用动，性价比很高。