简介：一套基于半不变量法实现概率潮流计算的Matlab源码，面向电力系统分析与新能源并网研究场景，适合需要评估负荷和电源出力随机性对电网运行状态影响的研究生、工程师及科研人员。压缩包共3个文件，全部为.m脚本，整体大小仅9KB，包含核心概率计算算法、IEEE 30节点标准算例数据以及潮流计算主程序，代码结构紧凑，便于直接运行与二次开发。这套源码已有1409人学习下载，验证了其在概率潮流入门与实践中的参考价值。通过研读源码，可掌握半不变量法建模步骤、随机变量概率分布处理以及Matpower接口调用技巧；结合附带算例可快速复现实验，直观理解蒙特卡洛模拟与解析法在概率潮流中的差异。同时，代码中的文件组织方式与变量设计也有助于读者学习如何搭建概率潮流实验框架，为后续扩展分布式能源不确定性分析提供实用工具。

1. 半不变量法概率潮流：不是不用仿真，而是把仿真做成了解析

电力系统规划里最磨人的问题，不是“某一种情况下潮流算不准”，而是“几百上千种情况下，电压和潮流的长尾巴到底长什么样”。光伏出力随风速摆动、负荷随时间波动、电动汽车充电桩接入位置不确定，每一样都让确定性潮流算出来的单点结果显得过于自信。半不变量法概率潮流，就是绕开“把每种场景都跑一遍潮流”的蒙特卡洛式暴力思路，直接用随机变量的统计矩做代数运算，用一次解析计算得到节点电压和支路潮流的概率分布。它解决的问题很具体：在和精确解相差不大的前提下，把计算量从成千上万次潮流削减到一次潮流加若干次矩阵运算。这套方法适合配电网规划、分布式电源接入评估、输电网静态安全分析，也适合手里没有超算、却需要在方案评审前给出概率结论的工程场景。

2. 为什么是半不变量法：原理与数学边界

2.1 概率潮流的两条技术路线：蒙特卡洛与解析法

先给一个判断框架。概率潮流目前就两条路线：模拟法和解析法。模拟法就是蒙特卡洛，把每个随机变量按分布抽样，组成成千上万个场景，每个场景做一次确定性潮流计算，最后对结果做统计。它直观、几乎不受模型复杂度限制，但代价也很直白——一次静态潮流大概几十毫秒到几百毫秒，抽五千次样本就是几十秒到几分钟，如果潮流迭代还不收敛，时间更不可控。在配电网三相不平衡模型、或者含大量逆变器控制策略的模型里，蒙特卡洛的时间成本和收敛问题会被放大得更明显。

解析法的思路完全不同。它不对随机变量抽样，而是把输入随机变量的统计特征（均值、方差、偏度、峰度之类）通过潮流方程传递到输出变量上，再重构出输出变量的概率密度函数或累积分布函数。半不变量法属于解析法里最工程化的一支，核心工具是半不变量（Cumulant，也叫累积量）。半不变量相比普通矩有个关键优势：独立随机变量之和的半不变量等于各自半不变量之和。这个性质让“多个节点注入功率共同作用”这件事变成简单的代数相加，不需要做卷积。加上Gram-Charlier级数或Cornish-Fisher级数做概率密度重构，计算速度通常比蒙特卡洛快两到三个数量级，而且结果是一个连续分布函数，不是一堆散点。

2.2 半不变量的级数展开与Cornish-Fisher改进

半不变量本身不直接给出概率密度，它给出的是分布的“特征参数”。拿到输出变量的各阶半不变量之后，要把它转换回概率密度函数或累积分布函数，这一步叫概率重构。常见做法有两种：Gram-Charlier级数和Cornish-Fisher级数。

Gram-Charlier级数思路是“用正态分布当基底，用高阶半不变量修正尾巴”。它把标准正态分布的概率密度函数和各阶导数做线性组合，前三项对应均值、方差、偏度，偏度修正了单侧尾巴，峰度修正了中心尖峰和两侧厚尾。实际使用中一般取到4阶或6阶，取太少尾部失真，取太多高阶半不变量本身误差会放大，出现负概率密度振荡。Cornish-Fisher级数走的是另一条路：它直接修正随机变量的分位数，从标准正态分布的分位数出发，通过高阶半不变量逐项修正，得到原分布的近似分位数。在重尾场景下，Cornish-Fisher级数给出的尾部分位数通常比Gram-Charlier更稳，因为它在累积分布函数的逆函数上做修正，跳过概率密度振荡问题。

选择依据并不复杂：输出变量接近正态分布时，Gram-Charlier级数取4阶就够；输出分布明显偏斜或重尾时，用Cornish-Fisher级数更稳。如果两种都试过都不理想，问题往往不在展开级数，而在前面的线性化假设，这个坑后面专门说。

2.3 计算主流程：从随机注入到节点电压的八步走

把半不变量法概率潮流的完整计算流程拆开，常见做法可以归纳成下面八步：

步骤一，确定系统的确定性工作点。把各节点负荷和电源出力设成期望值（通常是预测值），做一次确定性潮流计算，得到节点电压幅值、相角的基准值，以及潮流方程雅可比矩阵的逆矩阵。这个工作点是整个概率计算的参考系。

步骤二，建立输入随机变量的概率模型。对有功出力、负荷功率指定分布类型，常见的有正态分布、均匀分布，光照强度用Beta分布，风速用Weibull分布。每个随机变量需要至少前6阶半不变量，这是后面代数运算的原料。

步骤三，用Nataf变换或者Cholesky分解处理输入变量相关性。如果各节点负荷或电源之间存在相关性，这一步必须在半不变量求取之前完成，否则后续所有“半不变量相加”的运算都建立在错误假设上。

步骤四，把输入随机变量的半不变量从直角坐标转换到极坐标分量上。潮流计算里节点注入功率是复数，但概率计算要求半不变量在实部和虚部分别传递，这个细节容易漏。

步骤五，利用灵敏度矩阵（雅可比矩阵的逆）把节点注入功率的半不变量线性映射到节点电压实部和虚部的半不变量。这一步是整个方法的核心，本质是用一阶泰勒展开代替非线性潮流方程。

步骤六，由电压实部、虚部的半不变量求电压幅值和支路潮流的半不变量。电压幅值是实部虚部的非线性组合，需要根据工作点做线性化处理，支路潮流则用支路两端电压和导纳矩阵做相同映射。

步骤七，用Gram-Charlier级数或Cornish-Fisher级数重构输出变量的概率密度函数和累积分布函数。

步骤八，计算工程关心的指标：节点电压越限概率、支路潮流越限概率、期望值、标准差、一定置信度下的区间。这步直接用累积分布函数查值。

整个流程里，步骤一和步骤五牵涉最多工程判断，后面实现部分展开讲。

3. 用Python实现半不变量法概率潮流：最小可运行代码

3.1 输入数据组织：从单点潮流到随机模型

实现半不变量法概率潮流，不需要从零写潮流计算器，常见做法是直接复用成熟的潮流计算工具，取它的雅可比矩阵。下面代码用的是pandaPower做确定性潮流基准计算，然后自行提取节点注入和雅可比矩阵做概率扩展。

import numpy as np import pandapower as pp import scipy.stats as stats from scipy.sparse import csc_matrix, linalg as spla # 创建IEEE 30节点测试系统 net = pp.converter.from_mpc(pp.networks.case30()) pp.run_powerflow(net) # 提取确定性潮流结果：节点注入有功/无功、电压幅值和相角 pg = net.res_ext.loc[:, 'p_mw'].values # 母线有功注入 qg = net.res_ext.loc[:, 'q_mvar'].values # 母线无功注入 v0 = net.res_bus.loc[:, 'vm_pu'].values # 电压幅值标幺值 theta0 = net.res_bus.loc[:, 'va_degree'].values * np.pi / 180.0 # 相角弧度 # 提取雅可比矩阵（直角坐标形式） Ybus = net._ppc['internal']['Ybus'].toarray() V0 = v0 * np.exp(1j * theta0) S = pg + 1j * qg # 极坐标潮流方程：P = Re(V * conj(Y*V))，雅可比由pandapower内部提供 J = net._ppc['internal']['J'].toarray() # 极坐标雅可比矩阵

这段代码的前半部分是标准的潮流计算初始化，重点在最后两行：Ybus是全网络的节点导纳矩阵，J是潮流迭代收敛后的雅可比矩阵。半不变量法的线性化映射全靠这个矩阵，它的规模是2n乘2n（n为节点数），前n行对应有功方程对相角和电压幅值的偏导，后n行对应无功方程。

提示：如果用的不是pandaPower而是自己写的牛顿法潮流，雅可比矩阵一定要取最后一次迭代收敛时的值，不要用初始值。初始雅可比和收敛雅可比在某些系统里差别很大，直接影响概率结果的精度。

3.2 核心计算：线性化潮流与半不变量传递

拿到雅可比矩阵和基准工作点后，下一步是把输入随机功率波动的半不变量通过逆雅可比映射到状态变量上。这里有个关键点：潮流方程在极坐标下对电压幅值和相角求偏导，得到的雅可比矩阵是极坐标形式的，但半不变量法输入输出的随机变量一般以注入功率的实部和虚部表示，需要做坐标转换。

# 输入随机变量建模：以10个节点的负荷有功为例，设为正态分布 n = net.bus.shape[0] mu_p = pg.copy() # 各节点有功期望值 sigma_p = 0.03 * np.abs(mu_p) + 0.5 # 标准差：3%波动+0.5MW基础噪声 mu_q = qg.copy() sigma_q = 0.03 * np.abs(mu_q) + 0.2 # 用矩法求半不变量：先求中心矩，再转半不变量 def cumulants_from_normal(mu, sigma, order=6): """正态分布半不变量：前两阶为均值和方差，更高阶为0""" res = np.zeros(order) res[0] = mu res[1] = sigma ** 2 return res cum_p = np.array([cumulants_from_normal(mu_p[i], sigma_p[i]) for i in range(n)]) cum_q = np.array([cumulants_from_normal(mu_q[i], sigma_q[i]) for i in range(n)]) # 构造输入功率波动向量（有功+无功按节点顺序拼接） cum_in = np.vstack([cum_p, cum_q]) # 形状 (2n, order) # 灵敏度矩阵：逆雅可比矩阵取负（注入功率方程F(V)=S - Y*V中S为已知量） J_inv = np.linalg.inv(J) # 半不变量线性传递：输出半不变量 = (J_inv)^k 对输入半不变量的加权组合 # 一阶半不变量（均值）直接通过灵敏度映射 delta_cum_1 = J_inv @ cum_in[:, 0] # 二阶及以上的半不变量需考虑灵敏度矩阵的Hadamard幂 H = np.abs(J_inv) ** 2 # 对二阶取平方，三阶取立方，以此类推 cum_out = np.zeros((2 * n, 6)) cum_out[:, 0] = delta_cum_1 # 均值：线性映射 cum_out[:, 1] = H @ cum_in[:, 1] # 方差：平方叠加 cum_out[:, 2] = (np.abs(J_inv) ** 3) @ cum_in[:, 2] # 三阶半不变量 cum_out[:, 3] = (np.abs(J_inv) ** 4) @ cum_in[:, 3] # 四阶 cum_out[:, 4] = (np.abs(J_inv) ** 5) @ cum_in[:, 4] # 五阶 cum_out[:, 5] = (np.abs(J_inv) ** 6) @ cum_in[:, 5] # 六阶

这段代码的逻辑分三块。第一块是输入随机变量建模，工程上用期望值加波动幅度描述负荷和电源出力的不确定性。标准差设成“基准值的3%加常数项”，是为了避免零负荷节点出现零方差导致后续矩阵运算奇异。第二块是半不变量初始化，对正态分布来说前两阶是均值和方差，三阶以上全为零，但代码仍然保留到六阶，是因为后面级数展开需要统一格式。第三块是核心传递公式，半不变量线性传输遵循“阶数对应幂次”的规律：第k阶输出半不变量等于灵敏度矩阵元素的k次方与第k阶输入半不变量的加权和。

注意：这里用的是np.abs(J_inv) ** k而不是J_inv ** k，因为半不变量传递用的是泰勒展开系数的绝对值幂次。雅可比矩阵元素的正负号已经在潮流方程的方向性中体现，这里只取幅值贡献。

3.3 概率密度重构：Gram-Charlier与Cornish-Fisher

得到状态变量（电压实部、虚部）的半不变量后，需要把它们合成电压幅值和相角的半不变量，再用级数展开重构概率分布。电压幅值的平方等于实部平方加虚部平方，是非线性关系，工程上通常在工作点附近做一阶线性化处理，用复变函数的微分关系把实部虚部的半不变量合成幅值的近似半不变量。

# 合成节点电压幅值的半不变量（基于线性化） # 电压幅值增量 |dV| ≈ (Vr/Vm)*dVr + (Vi/Vm)*dVi Vr0 = v0 * np.cos(theta0) Vi0 = v0 * np.sin(theta0) coef_r = Vr0 / v0 # 实部系数 coef_i = Vi0 / v0 # 虚部系数 # 实际项目中由dP/dQ到dVr/dVi的映射已完成，此处直接合成 # 假设cum_out前n行为电压相角半不变量，后n行为电压幅值半不变量 cum_vm = cum_out[n:2*n, :] # 电压幅值半不变量（极坐标形式） # Gram-Charlier级数重构概率密度 def gram_charlier_pdf(x, cum, n_order=4): """基于半不变量重构概率密度，x为标准化的变量""" mu = cum[0] sigma = np.sqrt(cum[1]) gamma1 = cum[2] / sigma ** 3 # 偏度 gamma2 = cum[3] / sigma ** 4 - 3 # 超值峰度 std_x = (x - mu) / sigma phi = stats.norm.pdf(std_x) pdf = phi * (1 + gamma1 / 6 * (std_x ** 3 - 3 * std_x) + gamma2 / 24 * (std_x ** 4 - 6 * std_x ** 2 + 3)) return pdf # 用Cornish-Fisher求95%分位数（用于越限判断） def cornish_fisher_ppf(q, cum): """Cornish-Fisher分位数修正，q为概率值""" mu = cum[0] sigma = np.sqrt(cum[1]) z = stats.norm.ppf(q) gamma1 = cum[2] / sigma ** 3 gamma2 = cum[3] / sigma ** 4 - 3 # 一阶修正 z_cf = z + (gamma1 / 6) * (z ** 2 - 1) + \ (gamma2 / 24) * (z ** 3 - 3 * z) - \ (gamma1 ** 2 / 36) * (2 * z ** 3 - 5 * z) return mu + sigma * z_cf # 计算第7个节点的电压越限概率上限值 node_idx = 6 q95 = cornish_fisher_ppf(0.95, cum_vm[node_idx]) print(f"节点{net.bus.name[node_idx]}电压95%分位值为: {q95:.4f} p.u.")

这段代码展示了两件关键事。Gram-Charlier级数在标准化变量上展开，前四项对应正态基底、偏度修正、峰度修正，输出概率密度。Cornish-Fisher分位数公式则直接给出累积概率对应的电压值，这在工程上比概率密度更常用，因为越限判断本质是“P(V > Vmax) < 5%还是> 5%”的分位数问题。

注意：代码里cum_out[n:2*n]假设了状态变量排序是“相角在前、幅值在后”，实际工程中这个顺序取决于你的雅可比矩阵怎么排的。我踩过这个坑，排反了之后算出来的电压分布方差对不上蒙特卡洛，排查了半天才发现是索引顺序问题。

3.4 结果验证：与蒙特卡洛对拍

半不变量法算完之后，第一件事不是出报告，而是和蒙特卡洛结果对比验证。验证方法很简单：抽出同样的随机变量分布，跑2000次确定性潮流，统计电压幅值的均值和标准差，与半不变量法结果做相对误差对比。

# 蒙特卡洛验证：2000次抽样 np.random.seed(42) n_mc = 2000 vm_samples = np.zeros((n_mc, n)) for i in range(n_mc): # 对负荷波动抽样，重算潮流 net_mc = net.deepcopy() # 按正态分布给负荷添加波动 load_p = net_mc.load.p_mw.values * (1 + np.random.normal(0, 0.03, len(net_mc.load))) net_mc.load.p_mw = load_p try: pp.run_powerflow(net_mc) vm_samples[i, :] = net_mc.res_bus.vm_pu.values except: vm_samples[i, :] = np.nan # 剔除不收敛样本 valid = ~np.isnan(vm_samples[:, node_idx]) mc_mean = np.mean(vm_samples[valid, node_idx]) mc_std = np.std(vm_samples[valid, node_idx]) # 半不变量法结果 cum_vm_node = cum_vm[node_idx] scn_mean = cum_vm_node[0] scn_std = np.sqrt(cum_vm_node[1]) print(f"节点电压均值相对误差: {abs(scn_mean - mc_mean) / mc_mean * 100:.2f}%") print(f"节点电压标准差相对误差: {abs(scn_std - mc_std) / mc_std * 100:.2f}%")

对比的指标主要是均值和标准差，误差在1%以内说明线性化假设成立，误差超过5%基本可以判定系统运行点非线性过强或者波动范围太大，需要调整。这里有一个重要细节：蒙特卡洛抽样里要捕捉不收敛样本，半不变量法不会告诉你系统在哪些场景下潮流不收敛，这是它的盲区。如果蒙特卡洛里的不收敛率超过1%，说明系统静态安全储备不足，此时半不变量法的概率结果已经不可信，需要先做静态安全分析找出失稳场景。

4. 参数设置与工程边界：哪些值一变就翻车

4.1 随机变量建模：正态假设的适用范围

半不变量法本身并不要求随机变量服从正态分布，但很多工程实现默认套正态分布，原因有两个：正态分布的半不变量只需要均值方差，三阶以上全为零，计算量最小；节点负荷由大量独立小用户叠加，中心极限定理保证其近似正态。但现实中有三类场景正态假设撑不住：光伏出力的Beta分布、风速的Weibull分布、电动汽车充电负荷的混合分布。这些分布的特征是偏度明显、有厚尾，三阶以上的半不变量不可忽略。

处理办法不是换算法，而是把非正态分布的高阶半不变量求出来再代入流程。常见的做法是先用矩法求中心矩，再通过半不变量与中心矩的递推关系转换。对Beta分布和Weibull分布，可以直接查统计手册里的矩公式，或者用数值积分算前6阶矩。我在实际项目里遇到最多的是Beta分布的光伏模型，这个分布的形状参数需要从历史辐照度数据用极大似然估计拟合。拟合结果的精度对概率潮流的影响非常大，很多“算出来不对味”的案例，根源不在半不变量算法，而在输入分布参数估偏了。

4.2 节点注入相关性：Nataf变换与相关系数矩阵

概率潮流最容易被忽略的非理想因素是节点间的注入功率相关性。相邻光伏电站共用一片云层遮挡，出力强正相关；同一区域电动汽车同时充电，负荷也强相关。如果假设各节点独立，半不变量传递公式里“半不变量相加”直接成立。但相关性存在时，简单相加会低估输出变量的方差，尤其低估支路潮流的尾部风险。

工程上处理相关性的标准做法是Nataf变换。核心思路：输入变量的边缘分布可以是任意分布，先把每个变量用等概率变换转成标准正态分布，把原始相关系数矩阵修正为标准正态空间的相关系数矩阵，然后用Cholesky分解生成相关标准正态样本，最后逆变换回原始分布空间。在半不变量法里，这个过程不做抽样，而是把修正后的相关系数矩阵用于生成相关半不变量，本质上是对中心矩做修正。

需要注意的坑：Nataf变换要求输入的相关系数矩阵必须是正定的。实际工程里直接从历史数据算出的相关系数矩阵经常非正定，原因包括数据不完全同步、部分节点数据缺失、四舍五入误差。处理手段是特征值修正：把非正定矩阵的负特征值截断到零，再用原特征向量重构矩阵，或者直接在原始数据上做奇异值分解重构。千万别直接对非正定矩阵做Cholesky分解，会直接报错。

4.3 雅可比矩阵与系统状态：选工作点比选展开阶数更重要

半不变量法在数学上等价于“在工作点处用一阶泰勒展开代替潮流方程”，所以工作点选取直接决定精度。一个直觉判断：如果蒙特卡洛抽样得到的电压分布范围在±5%以内，一阶展开够用；如果电压波动超过±10%，线性化误差会显著增大，概率密度曲线出现变形甚至负值。

这里有个工程上常犯的错误：把基准工作点设成“当前时刻的运行状态”，而不是“未来场景的期望状态”。概率潮流的输入是预测分布，工作点必须与预测期望值一致。比如评估明天中午光伏大规模接入的电压分布，工作点应该设成中午时刻的光伏出力和负荷期望值，而不是当前时刻的实测值。否则均值偏移直接让方差和偏度全部失真。

另外，雅可比矩阵在重负荷节点附近可能接近奇异，逆矩阵元素非常大，半不变量传递会把微小的输入波动放大成巨大的输出电压波动，出现离谱的方差。这种情况下不能直接信任逆雅可比的结果，需要对系统做静态稳定性分析。实际项目里，若发现某个节点的方差比其他节点大一个数量级，优先检查该节点附近是否有接近电压崩溃的运行点。

4.4 输出变量选择：节点电压还是支路潮流

半不变量法的输出并不限于节点电压幅值，支路潮流、网损、变压器负载率都能算。但不同输出变量的线性化精度差别很大。支路有功潮流的传输方程是P_ij = V_i V_j (G_ij cos(theta_ij) + B_ij sin(theta_ij))，是节点电压幅值和相角的非线性函数，且支路两端节点电压相关性很强，在线性化过程中需要把两个节点的半不变量联合传递，而不是分别独立传递再相乘。

工程实践里的常见做法：先把节点电压实部虚部的半不变量全部算出来，再对每条支路构造“传输函数”，用链式法则把电压的半不变量映射为支路潮流的半不变量。这里的高阶交叉项（实部与虚部的乘积）非常多，处理起来比节点电压麻烦不少。如果项目里只需判断几条关键支路是否过载，建议只算这几条支路，不要全网络支路一起算，既减少计算量，也避开无关支路的线性化误差放大问题。

5. 半不变量法的常见问题与排查：现象、原因、解决

5.1 概率密度曲线出现负值

现象：Gram-Charlier级数重构出的概率密度函数在分布尾部变成负值，或者出现多峰振荡。

原因：级数截断或者高阶半不变量取值异常。三阶和四阶半不变量偏大时，修正项会在尾部分子上超过正态基底项，密度变负。另一个常见原因是输入随机变量取的正态分布假设与实际偏差较大，高阶半不变量不为零且数值过大，级数对非线性形状“修正过头”。

解决：先降展开阶数，从4阶退回2阶，看是否恢复为纯正态结果；若2阶也振荡，说明线性化出的半不变量本身有问题。再检查输入方差是否过大，把标准差从3%降到1%试算，若振荡消失，确定是大波动导致的线性化失效，需要改用分段线性化或直接上蒙特卡洛。也可以用Cornish-Fisher级数做交叉验证，分位数方法不显式构造概率密度函数，尾部表现更稳健。

5.2 高阶半不变量计算溢出或不准

现象：六阶半不变量达到1e6量级，重构出的分布宽得离谱，分位数超出物理极限。

原因：半不变量随时间阶数升高呈阶乘级增长，尤其对有偏分布，高阶半不变量剧烈放大。数值计算中用双精度浮点保存六阶以上结果，精度会严重丢失。

解决：工程上把处理重心放在四阶以内，六阶只用于尾部分位数计算且不做概率密度还原。改用法：所有半不变量在传递过程中统一做归一化处理，先行标准化再做级数展开，避免大数运算。如果使用的是电力系统计算平台，优先检查平台默认的半不变量阶数设置，很多商业软件默认取到四阶，强行调高不会带来精度收益。

5.3 计算结果与蒙特卡洛对不上：系统非线性过强

现象：节点电压均值对得上，但标准差的误差达到10%以上，或者95%分位数偏差明显。

原因：系统运行点附近潮流方程非线性强，或者输入波动范围太大，一阶泰勒展开不再成立。典型场景：重负荷线路接近热稳定极限、光伏出力波动大导致节点电压偏离基准值太多、含大量恒功率负荷的低压配电网。

解决：先检查蒙特卡洛样本里的电压分布范围，若电压偏移超过基准值10%，考虑用二阶半不变量方法，引入潮流方程的海森矩阵做二阶修正。但二阶方法在实现上需要构造三阶张量，代码量和调试难度翻倍。工程上还有个退路——把这个场景切分成几个子区间，每个子区间单独做半不变量法，再按子区间概率加权合并，这种方法虽然没有二阶方法理论优雅，但实施快很多。

5.4 相关系数矩阵不合法

现象：Nataf变换时Cholesky分解报错，提示矩阵不是正定矩阵。

原因：历史数据中部分节点量测缺失、时间序列不同步、或者多个节点之间存在接近完全的线性相关（比如两个光伏电站容量相同且共享辐照度数据），导致相关系数矩阵的特征值接近零或为负。

解决：先做特征值分解，把小于1e-6的特征值设为1e-6，重构矩阵再分解。若重构后矩阵仍然病态，检查数据是否包含重复测点或者强共线节点，直接剔除冗余节点。实际操作中，我一般对相关系数矩阵做一次L2正则化，给对角线加一个很小的常数，再检查条件数，如果条件数超过1e4级别，说明数据质量本身不适合做相关性分析，该回去补数据。

5.5 尾部分布误差大：换Cornish-Fisher后仍不行

现象：50%分位数附近结果极准，但1%和99%分位数的误差达到20%以上，“中间准、两边飞”。

原因：Cornish-Fisher级数本质上是对正态分位数做多项式修正，修正项越多，尾部收敛性反而越差，尤其对重尾分布会出现振荡。Gram-Charlier在尾部的表现类似，本质都是多项式逼近的Radon-Nikodym导数问题。

解决：把评估重点从“精确尾部”换成“越限概率级别”，也就是把问题从“99%分位数是多少”改成“越限概率是否低于1%”。这个切换在工程上意义重大，因为后者不需要精确尾部，只需要判断累积概率是否跨过阈值。若客户坚持要尾部精确值，没有捷径，老老实实上重要性采样蒙特卡洛，把抽样集中在越限区域。半不变量法在这个场景的真正价值是快速筛选，蒙特卡洛负责精算尾部。

6. 验证方法：用蒙特卡洛做基准校验与精度判定

半不变量法算完之后，总要回答“准不准”。最可靠的答案来自与蒙特卡洛的对比，但不是随便抽几百次样本就对拍，抽样规模直接决定“准”的结论可信度。

校验抽样规模有个经验公式：若目标是对比均值，500次样本就够；对比标准差，至少1000次；对比1%尾部分位数，5000次起步。原因是分位数估计的收敛速度远慢于均值和方差。验证时不光看误差均值，还要记录抽样波动区间，两次不同种子的蒙特卡洛结果本身就有±0.5%级别的标准差波动。如果半不变量法落在蒙特卡洛的置信区间内，才算真正对得上。

评价精度我习惯用三个指标：均值相对误差（小于0.5%算优秀）、标准差相对误差（小于2%算合格）、95%分位数绝对误差（小于额定电压的0.5%或者支路载流量的2%算合格）。这三个指标按项目需求分配权重，做调度运行方案用分位数指标为主，做规划评估用均值和分位数并重。

最后一个技巧：把半不变量法当作蒙特卡洛的前置预筛器。先用半不变量法快速识别越限风险节点，对低风险节点直接采信半不变量法结果，对高风险节点单独做蒙特卡洛精细计算。这样既保住计算速度，又把非线性影响最大的区域用高精度方法覆盖了。我自己做含高比例光伏接入的配电网项目时，通常用这种“解析初筛加蒙特卡洛精算”的组合，几十个节点系统几秒就能给出可信概率结论。

这套方法用顺了之后，回头看确定性潮流反而觉得可疑——一个数怎么能代表全部场景呢。但也要记住它的边界：系统运行于强非线性区时、含大量离散控制策略切换时，半不变量法的解析解会失真。判断是否该用它的标准很简单：先跑500次蒙特卡洛看电压波动范围，若波动小于10%且无潮流不收敛样本，放心用半不变量法；反之，老老实实加大蒙特卡洛规模。这个先验证后计算的习惯，帮我避掉过很多次“结果很漂亮但落不了地”的尴尬，希望帮到你。

本文还有配套的精品资源，点击获取