企业官网建设流程全解析

配电网里待久了就会发现，同一个节点不同时刻测到的电压根本不是一个固定值：白天屋顶光伏大发、负荷轻，电压往上飘；晚上负荷重、风场出力不足，电压又往下掉。传统确定性潮流在给定一组负荷和出力后给出一个确定解，但这个解并不能回答工程上最关心的问题——电压越限的概率有多大、支路过载的风险有多高。这类问题属于概率潮流（Probabilistic Load Flow）的范畴，而基于半不变量的解析法是求解它的高效手段之一。这篇文章以IEEE 34节点配网测试馈线为对象，讲解基于半不变量的概率潮流计算的数学原理、Matlab实现步骤、灵敏度矩阵生成、概率密度重构以及与蒙特卡洛模拟的对比验证，面向有潮流计算基础、想给配网不确定性分析建立落地工具的研究人员和工程师。

1. 概率潮流到底解决什么问题：一个配网工程的真实场景

1.1 确定性潮流在新能源接入仿真中的局限

我之前在某配电网项目里做分布式光伏接入分析，业主问的第一个问题往往不是"节点电压是多少"，而是"电压会不会越限，一年估计有多少小时会越限"。这个问题用传统潮流是答不出来的。传统潮流把负荷和出力写成一组确定数值，哪怕你按最恶劣工况校核一遍，也只能得到几个离散场景的结果，无法回答"概率"这两个字。

新能源接入之后，系统的不确定性来源明显增加：风速变化导致风电场出力波动、云层遮挡导致光伏出力跳变、电动汽车无序充电导致负荷随机性陡增。这些不确定性最终都会体现在节点注入功率的变化上，再通过潮流方程传播到节点电压和支路功率。要对这样的系统做精细化分析，就需要把输入量当成随机变量，让每个节点的电压和每条支路的潮流也变成有概率分布的随机变量，这种分析框架就是概率潮流。它的输出不是一组单一值，而是电压幅值、支路有功无功的概率密度函数（PDF）和累积分布函数（CDF），有了这些才能定量判断越限风险。

1.2 三类概率潮流算法的取舍

概率潮流的求解思路大体可以分成三类。

第一类是模拟法，最典型的是蒙特卡洛模拟。原理很直接：按输入变量的概率分布大量抽样，对每个样本做一次确定性潮流计算，最后统计状态量的分布。蒙特卡洛的优点是几乎不受系统非线性影响、实现简单、结果可靠；缺点是计算量太大，在含成千上万节点的系统上跑几万次潮流，时间成本很难接受。

第二类是近似法，代表性的是点估计法（Point Estimate Method）。它用少量确定性的采样点来近似输入分布的矩信息，只要做几次潮流就能估算出输出变量的统计矩。点估计法计算量小、也容易实现，但精度受采样点策略影响明显，对强非线性的支路潮流往往不够理想。

第三类是解析法，今天要重点讲的半不变量法就属于这一类。它的核心思路是把输入随机变量的概率信息变换成一串数字特征（各阶半不变量），再利用半不变量在独立随机变量求和下的可加性，直接由输入的半不变量线性组合出状态量的半不变量，最后用Gram-Charlier级数或Edgeworth级数重构概率密度函数。整个过程中确定性潮流只做一次，所以计算速度极快，在配电系统、输电系统的概率分析中都有很成熟的应用。

三类方法的定位差异用一张表可以看得很清楚。

我个人的选择习惯是：先用半不变量法快速扫一遍全网的薄弱节点，再用蒙特卡洛对少量高风险节点做校验，两者结合起来既不牺牲效率，又能兜住解析法尾部的误差。这个思路也是本文标题里把"半不变量"和"IEEE 34节点"放在一起的原因——先用经典算例把方法吃透，再迁移到实际系统。

2. 半不变量法的数学内核：可加性优势与级数展开

2.1 从特征函数到各阶半不变量

要说清楚半不变量法，必须先理解"半不变量"到底是什么。我们常用的数字特征是矩和中心矩：均值、方差、偏度、峰度等。半不变量（Cumulant）是另一套数字特征，它从特征函数中定义而来。

随机变量X的特征函数定义为φ(t)=E[e^{jtX}]。如果把对数特征函数展开成幂级数：

ln φ(t) = Σ κ_r (jt)^r / r!

展开系数κ_r就是X的r阶半不变量。实践中最常用的是前四阶，它们与均值μ、方差σ²、三阶中心矩m₃、四阶中心矩m₄的关系是：

• κ₁ = μ
• κ₂ = σ²
• κ₃ = m₃
• κ₄ = m₄ - 3σ⁴

前三阶半不变量和常规统计量是一一对应的，到了四阶开始出现"去正态化"的处理：把正态分布中本来应有的那部分四阶矩扣掉，剩下的才是真正反映非高斯特性的部分。这一点很重要，因为后面Gram-Charlier级数就是把分布看成"标准正态分布加上高阶修正项"来展开的，半不变量正是这些修正项的系数来源。

对于常见分布，半不变量可以直接写出来。正态分布只有前两阶非零：κ₁=μ，κ₂=σ²，三阶以上全部为0。这也是为什么如果一个系统的所有输入都假设成正态分布，输出分布也基本是正态的，Gram-Charlier级数用到二阶就足够了。而风电场出力往往用Weibull分布描述，Beta分布也常用于光伏出力，这些分布的高阶半不变量就要按定义去算或用数值积分求。

2.2 独立随机变量求和的可加性为什么关键

半不变量法之所以比直接的卷积法高效，根本原因在一个非常优美的性质上：如果随机变量X和Y相互独立，那么Z=X+Y的各阶半不变量等于X和Y对应阶半不变量之和，即：

κ_r(Z) = κ_r(X) + κ_r(Y), r=1,2,3,...

这个性质的意义怎么强调都不过分。在没有半不变量之前，求解两个独立随机变量之和的概率密度分布，需要做卷积运算。卷积对任意分布都是成立的，但计算复杂度高，Fortescue当时研究随机潮流时采用的正是这种逐步卷积的思路，数值上有大量近似。半个多世纪后，半不变量法的出现把整个问题从"卷积积分"降维成了"代数加法"。

在配网概率潮流里，一个节点的注入功率往往是多台分布式电源出力和多个负荷之和，它们如果被近似视为独立随机变量，那这个节点总注入功率的r阶半不变量就是每一项注入的r阶半不变量直接相加。这意味着我们可以把输入数据整理成"每个独立随机源"的形式，而不是强行给出节点总注入的具体分布，这在实际工程中非常方便。例如某节点接了一台风机和一个商业负荷，从数据采集的角度，我们更容易获得风机出力分布和负荷波动分布，而不一定拿得到两者之和的完整分布。半不变量的可加性让我们不需要做联合建模，直接加就好。

2.3 Gram-Charlier级数重构分布

知道了输出变量的各阶半不变量，下一步要从这些数字特征还原概率密度函数。这里用到的工具是Gram-Charlier级数展开。

首先把输出变量标准化，令z=(X-μ)/σ，令标准化后的半不变量为λ_r=κ_r/σ^r，其中λ₁=0，λ₂=1。然后利用标准正态分布密度函数φ(z)和它的各阶导数——即Hermite多项式H_n(z)——构造展开式：

f(z) = φ(z)[1 + (λ₃/6)H₃(z) + (λ₄/24)H₄(z) + (λ₅/120)H₅(z) + ((λ₆+10λ₃²)/720)H₆(z) + ...]

对应的累积分布函数为：

F(z) = Φ(z) + φ(z)[(λ₃/6)H₂(z) + (λ₄/24)H₃(z) + (λ₅/120)H₄(z) + ((λ₆+10λ₃²)/720)H₅(z) + ...]

前面几项Hermite多项式是：

• H₂(z)=z²-1
• H₃(z)=z³-3z
• H₄(z)=z⁴-6z²+3
• H₅(z)=z⁵-10z³+15z
• H₆(z)=z⁶-15z⁴+45z²-15

实际操作中需要取到几阶，取决于输入分布的非高斯程度。如果输入基本是正态的，取到二阶或四阶就够；如果风电占比高、分布偏斜明显，建议取到六阶甚至八阶。我通常默认取六阶，既能识别出明显偏度和厚尾，又不会因为高阶项振荡导致概率密度出现负值。需要注意的是，Gram-Charlier级数是渐近展开而不是收敛级数，取太多阶未必更好，在某些分布上可能反而出现负概率密度的现象，这一点在后面的避坑部分会细说。

3. Matlab实现框架：从原始数据到半不变量计算

3.1 整体模块划分与主流程

基于半不变量的概率潮流在Matlab里的实现并不复杂，整个算例的骨架可以划分成几个职责清晰的模块：

1. 数据准备模块：读取IEEE 34节点系统的母线、支路、变压器和负荷数据，形成潮流计算所需的结构体。
2. 确定性潮流模块：在期望值工况下用牛顿-拉夫逊法计算一次潮流，得到节点电压初值。
3. 灵敏度矩阵模块：从潮流收敛点处的雅可比矩阵求逆，得到节点注入功率扰动到节点电压扰动的线性映射。
4. 输入概率建模模块：为每个节点注入功率指定分布类型和参数，计算1~6阶半不变量。
5. 半不变量传播模块：利用灵敏度矩阵和可加性，把输入半不变量线性变换到节点电压和各支路潮流的半不变量。
6. 分布重构模块：调用Gram-Charlier级数函数，输出PDF、CDF和越限概率。
7. 蒙特卡洛参考模块：做若干次随机抽样潮流计算，用于对比验证。

主程序的执行逻辑大致如下：

%% 主流程：基于半不变量的概率潮流 (IEEE 34节点) % 1. 载入算例数据 [bus, branch, baseMVA] = load_ieee34bus(); % 2. 期望工况下的确定性潮流（牛顿-拉夫逊） [V0, theta0, J] = nr_powerflow(bus, branch, baseMVA); % 3. 由雅可比矩阵构造灵敏度矩阵 [S_VP, S_VQ] = calc_sensitivity_matrix(J, V0, theta0, bus); % 4. 输入注入功率的概率建模与半不变量计算 deltaP = 0.10; % 负荷波动标准差比例 kappaIn = build_input_cumulants(bus, deltaP); % 5. 半不变量传播到节点电压（幅值和相角） kappaV = propagate_cumulants(kappaIn, S_VP, S_VQ); % 6. Gram-Charlier级数重构输出分布并可视化 plot_pdf_cdf_from_cumulants(kappaV);

这里把"负荷波动标准差比例"设置成10%只是一个起点，实际场景中需要按历史数据确定。IEEE 34节点算例中一部分节点是分布式负荷，一部分是集中式负荷，建模时它们各自乘以不同的波动系数反而更贴近现实。

3.2 注入功率的概率建模与相关矩阵处理

输入概率建模是决定计算结果质量的第一步。常见做法是：对普通负荷注入功率，假设服从正态分布，均值取确定性潮流中的期望值，标准差取负荷均值的5%~15%的百分比；对风电场出力，则用Weibull分布刻画；对光伏出力，常用Beta分布。

在Matlab里我习惯用一个统一的函数接口来算任意分布的1~6阶半不变量：

function kappa = distCumulants(distType, params, order) % distType: 'normal','weibull','beta','uniform' % params: 相应分布的参数 % order: 最大半不变量阶数

对于正态分布，直接赋值均值、方差即可；对于非正态分布，我通常先数值积分求出1~6阶中心矩，再按矩-半不变量换算关系转为半不变量。这里的关键是矩-半不变量的换算关系表，前几阶可以手写，阶数高了就用递推公式。递推式如下：

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

其中m_r是r阶原点矩。这个递推式在实现时只需要几行循环，比硬编码到六阶的公式更通用。

如果多个随机源之间存在相关性，就需要先处理相关性。普通做法是：把相关输入变量的协方差矩阵做Cholesky分解，把一个相关正态变量向量Y变换成Y=LU，其中U是独立标准正态向量。由于半不变量法依赖"可加性"，而可加性对独立变量才严格成立，所以处理相关性时要么在建模阶段先做正交变换去相关，要么在传播公式里加入协方差的修正项。对初学者最友好的做法是先跑独立变量的版本，确认算法链路正确之后，再讨论相关性问题。

IEEE 34节点算例的母线数据里通常包含有功、无功的期望值，这些值可以直接作为正态分布的均值。我会额外把每个节点的注入功率设定成"负荷基准值加波动项"，以节点22为例，原算例负荷是(485kW, 190kvar)，那我设置注入的随机模型就是期望485kW、标准差48.5kW（10%波动）的正态变量，无功部分同样处理。这样传进半不变量计算函数时，每个节点的1~6阶半不变量就都有了。

4. 灵敏度矩阵的生成：连接随机注入与电压分布的桥梁

4.1 确定性潮流与雅可比矩阵

半不变量法虽然叫"概率潮流"，但它并不像蒙特卡洛那样把潮流方程迭代几万次。核心路线是：在期望运行点附近做一次牛顿-拉夫逊潮流，然后用一次线性化近似把随机注入功率的不确定性映射到状态变量上。因此这次确定性潮流的收敛质量非常关键，它决定了灵敏度矩阵所在的线性化工作点。

用极坐标形式描述潮流方程时，功率偏差同电压幅值和相角偏差之间的关系可以写成线性化形式：

ΔP = J_Pθ Δθ + J_PV ΔV ΔQ = J_Qθ Δθ + J_QV ΔV

这里J_Pθ、J_PV、J_Qθ、J_QV是牛顿-拉夫逊法雅可比矩阵的四个子块。写成矩阵形式就是：

[ΔP; ΔQ] = J [Δθ; ΔV]

在线性化假设下，节点注入功率的随机扰动[ΔP; ΔQ]和节点电压的随机扰动[Δθ; ΔV]就是通过雅可比矩阵J的逆来联系的。这也是从确定性潮流走向概率潮流的桥梁：潮流方程本身是非线性的，但我们在期望运行点附近用雅可比矩阵把它线性化了。对负荷波动不太大（比如标准差10%以内）的场景，这种线性化的精度是足够好的；但若新能源出力波动剧烈，线性化误差会变大，这也是半不变量法的固有边界。

4.2 状态变量和支路潮流的灵敏度传播系数

有了逆矩阵J⁻¹，理论上可以直接把注入功率的扰动映射成电压扰动：

Δθ = Σ S_Pθ ΔP + Σ S_Qθ ΔQ ΔV = Σ S_PV ΔP + Σ S_QV ΔQ

这里的S矩阵是J⁻¹的对应分块。在配网系统中大多数节点是PQ节点，所以计算相对直接。实际实现时要注意：J矩阵的维数是2n×2n，但单相潮流中PV节点对应的无功方程和平衡节点的两个方程，处理方式要跟牛顿-拉夫逊法完全一致，灵敏度矩阵的行列也必须严格对齐，否则映射关系会错位。

支路潮流的处理比节点电压多一个步骤。支路潮流是节点电压的函数，以从节点i流向节点j的有功功率为例：

P_ij = Re[V_i ( (V_i - V_j) y_ij + V_i y_sh )* ]

它在期望点对节点电压幅值和相角求偏导，得到一组系数C。然后支路有功功率对节点注入功率的灵敏度就是C与J⁻¹的乘积：

D_ij = C_ij * J⁻¹

这一步做完后，每条支路的半不变量传播公式就和节点电压完全一致了。我在代码里单独写一个函数来计算支路灵敏度矩阵，因为支路潮流更接近线性化极限，计算精度最容易在这里翻车。

灵敏度矩阵的实现代码片段：

function [S_PV, S_PQ] = calc_sensitivity_matrix(Jac, Vm, Vang, bus) % 输入: 雅可比矩阵Jac, 期望点电压幅值Vm(不含平衡节点), 相角Vang % 输出: 电压幅值对P/Q注入的灵敏度矩阵 n = length(Vm); nTot = 2*n; Jinv = inv(Jac); % 小系统可以用满矩阵求逆 % 按照状态变量排序，前n个为幅值，后n个为相角 S_PV = Jinv(n+1:end, 1:n); % ΔV -> ΔP S_PQ = Jinv(n+1:end, n+1:end); % ΔV -> ΔQ

这段代码只是示意，工程上建议用稀疏LU分解代替满矩阵求逆：

[L, U] = lu(Jac); S = -U\(L\eye(n)); % 用稀疏分解求解

我最初在IEEE 34节点算例中直接用了inv(Jac)，34节点规模不大，算起来很轻松，但后来把代码迁移到一个一百多节点的馈线上时，满矩阵求逆的效率问题就暴露出来了。所以从一开始写通用模块，还是用稀疏分解更稳妥。

5. 概率密度重构与误差对比：拿什么评价计算结果

5.1 Gram-Charlier级数求PDF与CDF的实现细节

当节点电压和各支路潮流的半不变量全部计算出来后，概率密度函数的重构就是一个纯粹的数学过程。假设已经得到某个节点电压幅值的各阶半不变量κ₁到κ₆，先算标准化半不变量λ，然后按Gram-Charlier级数叠加标准正态密度和Hermite多项式的组合。

实际编写重构函数时要注意几个细节。首先，节点电压的单位是标幺值，其均值往往在1.0附近、标准差可能只有0.01~0.02量级，这会导致各阶半不变量数量级差异巨大。因此标准化这一步不仅是为了级数展开的数学形式，也是数值稳定性的必要操作。其次，CDF的表格化和绘图需要在一个合理的区间内采样，通常选取μ±5σ的区间等分1000个点，点位太少曲线不光滑，太多则计算浪费。

function [x, fx, Fx] = gram_charlier(kappa) % kappa: 1~6阶半不变量 mu = kappa(1); sigma = sqrt(kappa(2)); lam = kappa ./ sigma.^((1:6)'); % 标准化半不变量 z = linspace(-5, 5, 1000)'; phi = normpdf(z); Phi = normcdf(z); % PDF修正项 pdfCorr = (lam(3)/6)*hermite(z,3) + (lam(4)/24)*hermite(z,4) + ... (lam(5)/120)*hermite(z,5) + ... ((lam(6)+10*lam(3)^2)/720)*hermite(z,6); fx = phi .* (1 + pdfCorr); % CDF修正项 cdfCorr = (lam(3)/6)*hermite(z,2) + (lam(4)/24)*hermite(z,3) + ... (lam(5)/120)*hermite(z,4) + ... ((lam(6)+10*lam(3)^2)/720)*hermite(z,5); Fx = Phi + phi .* cdfCorr; x = mu + sigma * z; end

这里Hermite多项式的函数需要注意定义域符号：不同教材的Hermite多项式存在归一化差异，对照公式时要仔细确认定义。我在初版代码里用过一种带系数sqrt的物理学家版Hermite定义，结果PDF整体偏了一个常数，排查了很久才发现是这里定义不统一。

5.2 与蒙特卡洛模拟的精度和时间对比

验证半不变量法是否可靠，最直接的办法就是和蒙特卡洛模拟对比。在IEEE 34节点算例上，我的做法是：用同一个输入概率模型跑20000次蒙特卡洛抽样，把统计出的电压均值、标准差、95%概率区间作为基准，再和半不变量法给出的结果对照。

典型的结果是：节点电压均值两者几乎一致，相对偏差基本小于0.1%；标准差的偏差通常也在1%~3%以内。PDF整体形状对比，正态波动场景下两条曲线几乎重合；在高波动场景（比如某个节点接入风电后注入方差增大）下，中心区域仍然贴合良好，但在概率区间的尾巴上，半不变量法会表现出一定偏差，具体表现为CDF尾部变厚或变薄，这取决于输入分布的偏度和峰度。

计算效率的差别才是最夸张的。同样在单台普通笔记本电脑上，确定性潮流一次大约耗时0.02~0.05秒，蒙特卡洛两万次就需要约8~15分钟；而半不变量法中所有随机性传播都是矩阵乘法和级数展开，整个过程只额外增加零点几秒的开销。也就是说，一次半不变量计算的时间大约等于一到两次确定性潮流。

这个效率优势带来一个实用价值：在需要反复调整输入参数做方案比选时，可以用半不变量法快速扫描出风险最高的几个节点和支路，再对这几个薄弱位置做蒙特卡洛精确校核。把两种方法放在一起用，既避开了解析法尾部的精度隐患，也避免了蒙特卡洛的全网大规模抽样成本，这是我认为最合理的工程实践姿态。

6. IEEE 34节点测试系统上的实测体会与避坑经验

6.1 测试馈线的特点与场景设置

IEEE 34节点测试馈线是配电网研究里很经典的公开算例，节点数为34，电压等级24.9kV，是一条辐射状配电馈线，包含若干子馈线、分布负荷和集中负荷，结构上比简单的辐射网更有代表性。我一直觉得它是做概率潮流的"黄金起跑线"：规模不算大、Matlab跑起来毫无压力，但又具备配网的典型特征，便于验证算法的通用性。

在这条馈线上做半不变量法概率潮流，我推荐的场景设置方式如下：

1. 保留全部34个节点，不做等值化简，这样灵敏度矩阵能覆盖每一条支路和每一个负荷点。
2. 对纯负荷节点，将负荷期望值设为基准工况数据，波动比例设5%~15%并衔接正态分布。
3. 在中段和后段挑两个节点接入模拟风电或光伏出力，用Weibull或Beta分布，渗透率控制在系统总负荷的适当比例（比如10%~30%），这样输出分布才会体现出非高斯特征。
4. 选几个薄弱节点（往往是线路末梢）重点观察，例如馈线末段的节点，这类节点电压波动最容易被放大。

负载率这个参数也很值得关注。在同样的波动比例下，负载率越高的节点，电压标准差的相对值越大，越限概率上升越明显。这些结论在CDF图上会非常直观，所以分析结果时不要只盯PDF的峰值，CDF的尾部曲线才是判断风险的核心。

6.2 半不变量法在使用中容易踩的坑

最后讲几个我在实际代码调试和算例分析中反复踩过的坑，这些内容常规论文里基本不会写。

第一，Gram-Charlier级数并不总是单调收敛的。当输入分布的非高斯性很强、或半不变量取到很高阶时，PDF曲线可能出现负值振荡，全是高阶项叠加过度导致的。遇到这种情况，先检查标准化半不变量有没有算对，再检查Hermite多项式的系数定义，最后才考虑降低截断阶数。我见过不少人一看到PDF出现负值就武断地说半不变量法不可用，其实多半是级数实现细节出了问题。

第二，灵敏度矩阵的边界性。半不变量法本质是线性化方法，输入波动越大，精度越差。当节点注入功率的标准差超过其均值的一半时，线性化的精度就很难保证了。这不是代码bug，而是方法本身的适用边界。遇到高波动场景，可以考虑采用蒙特卡洛校正或增加高阶修正项来处理。

第三，相关性问题不要一开始就碰。半不变量法最漂亮的性质是独立变量求和下的可加性，一旦引入输入变量相关性，传播公式就会出现交叉项，代码复杂度明显上升。初学者我建议先在独立假设下跑通整个链路，理解算法全貌后再用Cholesky去相关或协方差修正方案处理相关性。

第四，IEEE 34节点数据中KV基准和功率基准必须保持统一。原始数据里电压等级是24.9kV，但各段线路可能有不同的电压基准值，尤其在变压器支路两侧标幺化时要格外小心。这类基准值错误在确定性潮流里表现为电压分布异常，在概率潮流里就会表现为某个节点的电压PDF和蒙特卡洛结果差出一大截，极其迷惑人。

从编程习惯上讲，建议在程序开头就做一个简单的校验：先跑一次确定性潮流，和已知基准结果比对电压幅值，误差小于1e-6再继续后续的概率计算。这个"先定确定性潮流、再做随机分析"的次序，是我在多个算例中验证过的最稳定的工程流程。半不变量法的代码实现本身并不长，难点在于每一步都必须心中有数——知道这个矩阵从哪里来、这个级数项在描述什么，而不是盲目套用模板。把原理、实现、验证这条链路打通之后，这套方法再去迁移到更大规模的系统或更复杂的多能流场景，就会顺畅很多。