如果你平时只做确定性潮流计算，第一次听到“概率潮流计算”这个概念，脑子里大概率会冒出一串问题：潮流结果还能有概率分布？多跑几次仿真不就行了？半不变量又是什么东西，跟蒙特卡洛比到底强在哪？这套基于半不变量的随机潮流方法，本质上是把“输入随机性”通过潮流方程映射成“输出随机性”的一种解析近似手段，计算速度快、占用资源少，特别适合用在IEEE34节点这类中等规模系统上做电压越限概率评估、新能源出力波动影响分析。
这篇文章我按自己从理论推导到Matlab代码落地、再到IEEE34节点算例实测的完整过程来写。适合正在做随机潮流方向课题、毕业论文涉及概率潮流计算、或者工作中需要评估配电网/输电网电压风险的工程师。我会把半不变量的数学原理、潮流方程线性化处理、Cornish-Fisher级数重构分布、关键代码模块、以及对标蒙特卡洛的精度校验一次讲透，顺便把我在调试中踩过的坑也一并交代。
1. 确定性潮流只给一个点，随机注入下你需要一条曲线
1.1 确定性潮流能回答什么，不能回答什么
传统潮流计算给定一组确定的负荷和发电出力，求解出节点电压幅值、相角和支路功率，得到的是一个确定性的运行点。比如IEEE34节点系统在某组负荷下，节点20的电压是0.9825 pu，这就是一个点。这个点有没有用？有，但它回答不了工程上几个很现实的问题：未来半小时负荷波动，节点电压低于0.95 pu的概率是多少？光伏出力从200 kW跳到800 kW，馈线末端电压越上限的风险有多大？
这些问题本质上都涉及随机性。电网里的负荷在变，新能源出力在变，故障和机组停运也是随机的。确定性的潮流解只是把随机输入取了一组期望值或典型值，算出来一个“平均场景”的结果，但平均场景往往掩盖了极端风险。
举个生活化的例子，一个城市的日用电量均值是1000 MWh，但峰值日可能到1300 MWh，谷值日只有700 MWh。如果你只看均值，配电变压器的容量规划就会出大问题。同样的道理，电力系统运行决策需要的是“分布”而非“单点”。
1.2 概率潮流的常规解法与效率对比
概率潮流（Probabilistic Load Flow, PLF）就是为了解决这个问题出现的。它把节点注入功率视作随机变量，计算节点电压、支路潮流的概率分布。目前主流做法大概分三类。
蒙特卡洛模拟（MCS）是最直观的思路：按输入随机变量的分布抽样成千上万组，每组都做一次确定性潮流，最后统计输出量的分布。好处是原理简单、几乎没有线性化假设，坏处是慢。IEEE34节点做一次牛顿-拉夫逊潮流大约几十毫秒，抽一万次就是几百秒，用来做在线评估肯定不现实，但作为验证基准非常合适。
点估计法（Point Estimate Method）的思路是取输入变量的若干特征点，每个特征点跑一次潮流，然后用加权组合还原输出矩信息。计算量是2n+1次潮流（n为随机输入维数），比蒙特卡洛快很多，但精度对特征点选取敏感。
半不变量法走的是另一条路：不做大量采样，而是把输入随机变量的各阶半不变量求出来，通过潮流方程的线性化灵敏度映射到输出变量，再用级数展开重构概率密度函数和累积分布函数。整个过程只需要一次基态潮流加几次矩阵运算，毫秒级出结果。下表是三种方法的对比。
| 方法 | 计算代价 | 精度特征 | 适用场景 |
|---|---|---|---|
| 蒙特卡洛 | 数千到数万次潮流 | 精度随样本数提升，无线性化假设 | 基准校验、高精度分析 |
| 点估计法 | 2n+1次潮流 | 前三阶矩精度较好 | 中等规模快速评估 |
| 半不变量法 | 1次潮流+解析计算 | 小扰动下精度高，重载/大波动下精度下降 | 在线评估、海量场景遍历 |
我做IEEE34节点的随机潮流时，首选半不变量法，然后用蒙特卡洛结果做交叉验证。这样既保证了速度，又能知道近似方法到底损失了多少精度。
2. 半不变量法的数学内核：从随机注入到输出分布的两步映射
2.1 半不变量的定义和“加法友好”特性
很多资料一上来就抛公式，容易把人劝退。我尽量说得直白一些。半不变量（Cumulant）是随机变量的一种数字特征，跟矩类似，但有一个非常好的性质：如果两个随机变量相互独立，那么它们之和的各阶半不变量等于各自同阶半不变量之和。
这跟均值和方差的性质是统一的。一阶半不变量就是均值，二阶半不变量就是方差。均值和方差确实满足“加法可加性”，但三阶以上中心矩不满足。半不变量就把这套性质推广到了任意阶。定义上，它是特征函数对数展开后的系数，计算上可以由中心矩递推得到。
几个常用的对应关系：
- κ₁ = μ（均值）
- κ₂ = σ²（方差）
- κ₃ = E[(X-μ)³]（三阶中心矩，反映偏斜）
- κ₄ = E[(X-μ)⁴] - 3σ⁴（四阶累积量，反映峰度）
工程计算中一般取到四阶就够用，再高阶的数值稳定性差，而且对分布形状的改善有限。
这里要强调一下为什么“半不变量”适合做概率潮流。潮流计算中，某个输出量（比如节点电压幅值）是大量输入随机变量共同作用的结果，它可以近似看成输入变量的线性组合。如果输入变量相互独立，输出变量的半不变量就能用输入变量半不变量的加权和直接算出来，不需要先求联合分布，也不需要做卷积。这就是半不变量法计算效率的核心来源。
2.2 潮流方程线性化与灵敏度映射
潮流方程本身是非线性的，但在基态运行点附近可以线性化。极坐标形式的潮流方程为：
P_i = V_i Σ V_j (G_ij cosθ_ij + B_ij sinθ_ij) Q_i = V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij)
写成矩阵形式，线性化后有：
[ΔP; ΔQ] = J · [Δθ; ΔV]
其中J是牛顿-拉夫逊法的雅可比矩阵。做概率潮流时，我们关心的是相反的方向：输入注入功率扰动ΔS = [ΔP; ΔQ]，会引起状态变量变化ΔX = [Δθ; ΔV]。两边求逆：
ΔX = J⁻¹ · ΔS
J⁻¹就是灵敏度矩阵。它把输入随机变量的随机性“传导”到输出状态变量上。对于第i个输出状态变量，有：
ΔX_i = Σ_j S_ij · ΔS_j
其中S_ij是灵敏度矩阵第i行第j列的元素。这就是半不变量法实现“从输入到输出”映射的关键一步。如果输入变量相互独立，输出的r阶半不变量就是：
κ_r(X_i) = Σ_j (S_ij)^r · κ_r(S_j)
注意这里有个r次方。一阶时对应均值映射，二阶时对应方差映射（平方项），三阶以上同理。这个式子看着简单，却是整个算法的心脏。
2.3 用Cornish-Fisher级数重构概率分布
半不变量不是分布本身，用它还原概率密度函数（PDF）和累积分布函数（CDF）还需要一步。常见做法有两种：Gram-Charlier级数和Cornish-Fisher级数。
Gram-Charlier级数以标准正态分布为基础，用Hermite多项式叠加修正项逼近真实密度函数。优点是概念直观，缺点是某些情况下密度函数会出现负值或尾部振荡，阶数越高越不稳定。
Cornish-Fisher级数是基于标准正态分布分位数的一种展开，直接修正分位数，进而获得CDF的反函数。工程上我更喜欢用它，因为求越限概率本质上就是在查CDF的反函数，Cornish-Fisher给的结果更稳定。
Cornish-Fisher三阶截断公式大致形式：
y(α) = z(α) + (γ₁/6)(z(α)² - 1) + (γ₂/24)(z(α)³ - 3z(α)) - (γ₁²/36)(2z(α)³ - 5z(α))
其中γ₁是三阶标准化半不变量（偏度），γ₂是四阶标准化半不变量（超额峰度），z(α)是标准正态分布的α分位数。得到标准化的y后，再反变换回实际变量：
X(α) = μ + σ · y(α)
这样就能算出任意分位点对应的电压值。把分位点α在[0.01, 0.99]区间内密集取值，就能画出CDF曲线。如果想要PDF，可以对CDF做数值微分，或者直接用Gram-Charlier展开密度。
3. 从IEEE34节点系统开始：数据准备与随机模型设定
3.1 为什么要选IEEE34节点作为验证平台
IEEE34节点系统是一个经典的测试算例，规模适中，既能体现算法在中等规模系统中的表现，又不会像上百节点系统那样调试困难。MATPOWER工具包中自带的case34是34节点输电网算例，数据文件可以直接用loadcase命令载入，包含34个节点、33条支路和4台发电机。对于做随机潮流算法验证来说，这个规模恰到好处。
IEEE34节点还有一个更知名的版本是34节点配电网测试馈线（IEEE 34-bus test feeder），属于美国亚利桑那州一条实际馈线的简化模型，电压等级24.9 kV/4.16 kV，特点是供电半径长、负荷分散、线路阻抗大，末端电压问题突出。如果你想专门做配电网电压越限分析，用配电网版本更有代表性。但无论是输电网版本还是配电网版本，半不变量法的算法框架完全一致，区别只在于输入数据和随机变量的设置方式。
我自己的做法是用MATPOWER的case34打好算法框架，再导入配电馈线数据跑配电网场景。这样一套代码两种系统都能测，效率很高。
3.2 负荷随机模型和参数选择
在随机潮流中，负荷是最基础的随机源。工程上最常用的假设是节点有功和无功负荷服从正态分布，均值取基态潮流中的确定值，标准差取均值的5%~10%。
比如节点k的基础有功负荷是1.5 MW，取标准差为5%，则：
P_k ~ N(1.5, 0.075²) MW
无功负荷Q_k类似，可以取相同的相对标准差，也可以用功率因数关联。需要注意两点。第一，P和Q通常不是独立的，实际中常用恒功率因数假设，即Q_k = P_k · tanφ，这样Q的随机性由P决定，两个输入变量之间引入了相关性。第二，PV节点（发电机节点）的注入有功一般设为恒定，无功出力是为了维持电压，它的随机性来自对端负荷的随机波动，不能简单当作独立随机变量处理。
最简单可靠的建模方式是把所有PQ节点的P和Q都视为独立正态随机变量，标准差取5%，这样能充分利用半不变量法的可加性，避免相关性处理的额外复杂度。在此基础上再逐步引入相关性扩展，是比较稳妥的学习路径。
对于新能源出力，可以用Beta分布或Weibull分布建模，不一定要求正态。半不变量法的一大优势就是它不要求输入变量服从正态分布，只要你能算出一到四阶半不变量就能参与运算。Beta分布和Weibull分布的半不变量可以通过数值积分或者先求矩再转换得到。
3.3 数据单位换算：最容易翻车的细节
这里必须单独说一个坑。MATPOWER中所有数据默认使用标幺值，基准功率baseMVA在case34中通常是100 MVA。case34中某节点负荷显示为1.5 MW，在数据文件里实际上是0.015 pu。如果做随机模型时直接用1.5作为均值、0.075作为标准差，半不变量法算出来的电压波动会大得离谱，因为电压幅值本来就在1.0 pu附近，输入注入的波动单位搞错整个结果直接报废。
统一的做法是先把负荷的有名值除以baseMVA变成标幺值，再做均值、标准差的设定。我一开始自己写代码时就吃过这个亏，后来把所有参数都放在一个结构体里统一换算，就再没出过问题。
4. Matlab逐模块实现：从潮流内核到分布重构
4.1 模块一：确定性潮流与雅可比矩阵获取
半不变量法的第一步是获取基态潮流结果和雅可比矩阵。我的实现没有直接依赖MATPOWER的runpf黑盒，而是自己写了一个精简的牛顿-拉夫逊潮流内核。原因很简单：runpf的返回结果里不直接提供雅可比矩阵，虽然可以通过修改MATPOWER内部函数或利用option打入补丁获得，但总归不够干净。自己实现NR潮流对IEEE34节点这种规模并没有难度，而且雅可比矩阵的获取完全在掌控之中。
核心流程如下，首先是构建导纳矩阵：
% 基于bus和branch数据构造导纳矩阵Y Y = makeYbus(mpc.bus, mpc.branch);然后是牛顿迭代主循环：
for iter = 1:maxIter [dP, dQ] = powerMismatch(V, theta, Y, bus, gen); J = computeJacobian(V, theta, Y, bus); dX = J \ [dP; dQ]; theta = theta + dX(1:nPV+nPQ); V = V + dX(nPV+nPQ+1:end); if max(abs([dP; dQ])) < tol break; end end这里的关键是节点类型的编号重排。平衡节点（slack bus）的相角是参考基准，不参与迭代；PV节点的电压幅值给定，只参与相角迭代；PQ节点的电压幅值待求。在组装雅可比矩阵时，要把这些索引关系理清楚，不然矩阵维度对不上。
得到收敛结果后，J就是灵敏度计算所需的基础矩阵。如果你不想自己写完整的NR潮流，也可以在MATPOWER安装目录中找到内部潮流函数，在迭代结束处把J返回出来，效果一样，但要注意版本兼容性。
4.2 模块二：输入随机变量的各阶半不变量
假设系统中有NLoad个PQ节点，每个节点的有功、无功负荷都视为正态随机变量。各阶半不变量计算非常直接：
% 输入注入随机变量个数 nInput = 2 * NLoad; % 半不变量矩阵：4行，对应1~4阶 kappaX = zeros(4, nInput); for k = 1:NLoad % 有功注入，一阶半不变量=均值(标幺值)，二阶=方差 kappaX(1, 2*k-1) = P_load_pu(k); kappaX(2, 2*k-1) = (sigma_pu(k))^2; % 三阶、四阶为0（正态分布） % 无功注入 kappaX(1, 2*k) = Q_load_pu(k); kappaX(2, 2*k) = (sigma_q_pu(k))^2; end如果负荷模型用的是Beta分布或其他分布，就需要先计算中心矩再转换成半不变量。矩阵形式里每列对应一个输入随机变量，每行对应一阶到四阶。这个组织结构在后面聚合运算时非常方便。
4.3 模块三：输出状态变量的半不变量聚合
在基态潮流收敛后，从雅可比矩阵J求逆得到灵敏度矩阵S_inv。由于状态变量和输入注入的排列顺序需要在代码中对齐，建议把S_inv拆成两部分：S_theta对应相角输出，S_V对应电压幅值输出。各自的第r阶输出半不变量计算方式为：
% kappaX_r：所有输入变量的第r阶半不变量列向量 % S_V_abs：电压幅值对输入注入的灵敏度矩阵绝对值 kappaV_r = (abs(S_V_abs).^r) * kappaX_r;用矩阵运算一次完成所有节点电压第r阶半不变量的聚合，代码非常简洁。四个阶数分别算一遍，就得到了全部34个节点的电压幅值一阶到四阶半不变量。
值得注意的一个工程细节是，灵敏度矩阵的量级会直接影响结果。如果某行元素特别大，说明对应节点对某个注入的随机扰动特别敏感，往往意味着该节点电气距离远、网络支撑弱。在IEEE34节点配电网版本中，馈线末端节点的电压灵敏度通常明显高于靠近电源侧的节点，这也是馈线末端电压波动大的数学体现。
4.4 模块四：Cornish-Fisher重构与可视化
有了节点电压的前四阶半不变量，就可以用Cornish-Fisher级数重构CDF。以节点20为例，核心代码如下：
mu = kappaV(1, i); % 均值 sigma = sqrt(kappaV(2, i)); % 标准差 gamma1 = kappaV(3, i) / sigma^3; % 偏度 gamma2 = kappaV(4, i) / sigma^4; % 超额峰度 alpha = (0.001:0.001:0.999)'; % 分位点序列 z_a = norminv(alpha); % 标准正态分位数 y = z_a + (gamma1/6) .* (z_a.^2 - 1) ... + (gamma2/24) .* (z_a.^3 - 3*z_a) ... - (gamma1^2/36) .* (2*z_a.^3 - 5*z_a); V_cdf = mu + sigma * y; % 累计分布函数曲线这里V_cdf和alpha一一对应，plot(V_cdf, alpha)就是电压幅值的CDF曲线。如果想看PDF，对V_cdf的差分取倒数即可：
dV = diff(V_cdf); pdf_approx = 1 ./ dV / (length(alpha) - 1);把34个节点的计算结果循环跑一遍，整个系统的概率潮流信息就全出来了。最让人惊艳的是，从负荷建模到CDF重构，整个流程在普通笔记本上运行时间不超过1秒，而同样精度的蒙特卡洛基准试验要跑小十几分钟。
5. 仿真结果怎么看：输出、校验与误差分析
5.1 从分布曲线到越限概率
半不变量法跑通后，最直观的输出是各节点电压幅值的概率密度曲线和累计分布曲线。以IEEE34节点中靠近馈线末端的节点为例，在负荷标准差取5%时，电压幅值均值大约在0.98 pu附近，标准差在0.002~0.005 pu这个数量级。看到这个数值要注意，电压标幺值的标准差很小，说明5%的负荷波动对电压造成的绝对影响也就千分之几pu，符合实际物理规律。
更值得关注的是分布的“形状”。如果电压概率密度曲线明显偏向一侧，说明运行点附近电压对负荷扰动的响应是非线性的，偏离正态分布。这就是半不变量法相对简单均值方差分析的优势：三阶半不变量反映了偏斜方向，四阶半不变量反映了尾部厚度。重负荷场景下，电压分布会呈左偏形态，即低压尾部被拉长，电压越下限的风险不能用对称分布来估计。
越限概率的计算也很直接。比如想知道节点20电压低于0.95 pu的概率，只需要在CDF曲线上找到0.95 pu对应的CDF值：
P_low = interp1(V_cdf, alpha, 0.95, 'linear');这个数字就是“电压越下限风险”。做运行方式调整或新能源接入容量分析时，把不同场景下的越限概率画成柱状图，哪些节点危险、哪几类场景风险大，一目了然。
5.2 与蒙特卡洛对比的校验指标
只跑半不变量法不跟蒙特卡洛对比，很难判断算法实现和参数设置是否正确。我的校验流程是固定的：对同样的随机模型，用蒙特卡洛抽5000~10000组样本，每组做一次确定性潮流，统计电压幅值的均值、标准差和CDF，然后与半不变量法结果对比。
在负荷标准差为5%时，两者CDF曲线几乎重合，最大偏差通常在1e-3 pu以内，完全满足工程分析需求。负荷标准差增加到10%时，偏差会有所上升，尤其是分布的尾部。原因不难理解，半不变量法依赖潮流方程在基态点附近的线性化，负荷波动越大，线性化误差越明显。记住这个误差规律，以后用半不变量法就知道它的适用边界了。
我习惯用两个定量指标来评估精度：
- 平均绝对误差（MAE）：所有分位点上CDF差值的平均值。
- 最大绝对误差（MaxAE）：所有分位点上CDF差值的最大值，通常出现在分布尾部。
实测下来，小波动场景下MaxAE数量级为10⁻³，大波动场景会到10⁻²。如果出现这个量级的退化，优先检查是不是灵敏度矩阵在重载点出现了病态，再考虑改用更高阶的非线性映射或者切换到蒙特卡洛。
6. 我在调试中踩过的坑和最后几点建议
6.1 雅可比矩阵病态与灵敏度矩阵爆炸
第一次跑通半不变量法后，我遇到的最隐蔽的问题来自雅可比矩阵的数值状态。在接近重载的运行条件下，雅可比矩阵可能病态，求逆后灵敏度矩阵中某些元素数值异常大，导致电压方差被严重高估。
排查方法其实不难：就算完成基态潮流后，用cond(J)看一眼条件数。条件数超过1e12就要警惕。针对这种情况，我做了两件事。第一，检查基态潮流本身是否已经接近收敛极限，如果NR迭代次数很多或残差下降缓慢，先解决潮流收敛问题再说概率计算。第二，将对角占优性较差的节点负荷适当降低波动幅度，让运行点远离电压崩溃边界。
还有一个小技巧，不要直接对完整的雅可比矩阵求逆，而是用稀疏矩阵的分解结果来解线性方程组。Matlab中只需要写成S_full = J \ eye(n)，得到的灵敏度矩阵就是J⁻¹。这样数值稳定性更好，也比inv(J)快不少。
6.2 Cornish-Fisher级数的阶数选择
半不变量法重构分布时，展开阶数不是越高越好。我做过一个对比实验，用三阶截断、四阶截断和六阶截断分别重构同一个算例的电压分布。结果很有意思：四阶截断的尾部明显比三阶更好，但六阶截断出现了小幅振荡，在一些分位点上的CDF值甚至超过了[0,1]范围。
原因在于高阶半不变量的数值精度有限，随着阶数升高，估计误差会被放大。工程上建议最多用到四阶，即只保留偏度和峰度修正。如果四阶截断仍然不满足精度要求，问题大概率不在展开阶数上，而是半不变量法本身对本文场景的线性化假设已经不成立了。
6.3 扩展到含新能源随机出力的思路
最后聊聊扩展。我在IEEE34节点上跑通负荷随机模型后，又加入了风机和光伏出力随机性。新能源出力的分布明显不是正态的，风机出力常用Weibull分布描述，光伏出力常用Beta分布描述。但半不变量法依然适用，只需把新能源节点的注入功率半不变量算出来，追加到输入随机变量矩阵kappaX中即可。
计算Weibull分布或Beta分布的半不变量，最省事的办法是先数值积分求原点矩，再通过原点矩与半不变量之间的递推关系转换。Matlab中有自带的wblstat、betastat函数可以直接获得均值和方差，更高阶矩需要自己积分。我自己写了一个小函数，输入分布类型和参数，输出前四阶半不变量，几十行代码就搞定了。
整个IEEE34节点随机潮流项目做完，我最深的体会是：概率潮流的价值不在于算得多花哨，而在于把“风险”这个东西量化出来了。确定性潮流告诉我们电压是多少，概率潮流告诉我们电压超出安全范围的概率是多少。基于半不变量的方法虽然有一些线性化近似，但胜在速度快、结果稳定，在需要遍历大量场景的在线分析和规划评估中，实用性远高于蒙特卡洛。如果你也想在Matlab里实现这套算法，建议按文中的模块顺序一步步搭建，每完成一个模块就用蒙特卡洛交叉验证一下，这样出了问题时能很快定位到具体环节。
