1. 确定性潮流为什么“不够用”：概率潮流的工程需求

1.1 光伏不确定性给潮流计算带来的实际冲击

以前做潮流计算，输入一张负荷表、确定发电机出力，跑一次牛顿法，拿到一组节点电压和支路功率就可以交差。这套流程在光伏大规模接入之后变得非常别扭：光伏出力不是固定值，晴天、阴天、云层飘过、早晚光照变化，都会让同一节点的注入功率在几分钟内发生大幅波动。你在下午三点用预测均值跑出来的潮流，在下午三点半可能完全不成立，因为几十兆瓦出力消失，电压分布随之整体变化。

这种波动不是偶然事件，而是常态。以常见的标准节点测试系统为例，在某个负荷节点附近接入一个中等容量的光伏电站，若当天气候条件一般，出力可能在额定值的20%到90%之间来回摆动。确定性潮流只能回答“在这个给定的注入下，系统状态是多少”，而工程上真正需要回答的问题是“在这种随机注入下，节点电压有多大的概率越限，支路潮流有多大的概率过载”。这就是概率潮流计算的出发点。

1.2 概率潮流到底要算哪些指标

概率潮流和确定性潮流的核心区别，是把输入的随机性转化为输出的概率分布，而不是只给一个点值。具体要输出的指标大致有下面几类：

举个例子：某节点在确定性潮流里电压是1.04pu，看着很安全，但概率潮流算出来该节点电压高于1.05pu的概率是12%。这么一来，调度员就不会只盯那个“平均值”，而是会关注这个12%的发生条件，比如发生在光伏大发时段，那就要提前安排调峰或无功控制。所以概率潮流对规划的意义，不是给出一个更精确的数字，而是给出一个可以用于风险决策的概率指标。

1.3 两条技术路线的选型思路

现在工程和学术圈里实现概率潮流的方法很多，但最常被拿来组合的，就是蒙特卡洛模拟和半不变量法。蒙特卡洛模拟的思路最直接：反复抽样、反复求解潮流，最后统计结果分布。它的优点是几乎不依赖近似假设，对非线性潮流方程、任意输入分布都能处理，精度高，但代价是计算量大，几千上万次潮流求解在大型系统里非常耗时。

半不变量法则走的是另一条路：在基准运行点把潮流方程线性化，然后用随机变量各阶半不变量在灵敏度矩阵上的叠加，直接算出输出变量的半不变量，最后用级数展开还原概率分布。它的计算量小到可以忽略，但代价是它依赖于线性化和级数截断这两个近似环节。两条路线各有利弊，实际工程中我习惯把蒙特卡洛当成“验收标准”，把半不变量法当成“快速筛选工具”，两种方法互相校验。下面的章节会把两条路线的原理和实现细节都拆开讲清楚。

2. 光伏出力与负荷的不确定性建模：Beta分布是常见起点

2.1 从光照强度到光伏有功出力

处理光伏不确定性，第一步是给光伏出力选择一个合适的概率分布模型。光伏有功出力本质上来自光照强度，光照强度受天气、季节、云层遮挡影响，经验统计显示，在一段固定时间尺度内，归一化的光伏出力往往呈现中间高、两端低的分布特征，Beta分布是业内最常用的近似模型。

设某光伏电站的额定有功为P_rated，归一化出力为p = P_pv / P_rated，则p在[0,1]区间内服从Beta分布，其概率密度函数是：

f(p) = Γ(a+b) / (Γ(a)Γ(b)) · p^(a-1) · (1-p)^(b-1)

其中a和b是形状参数。a和b的取值直接决定分布的形态：a和b都偏大时，分布向中间集中；a较小而b较大时，分布偏向低出力区间。实际建模时，一般根据该电站历史出力的均值和方差，用矩估计法计算这两个参数。设实测归一化出力的均值为μ、方差为σ²，则参数估计公式为：

a = μ · (μ(1-μ)/σ² - 1)

b = (1-μ) · (μ(1-μ)/σ² - 1)

如果没有历史数据，也可以直接取一组经验参数，比如a=5、b=2，表示光伏出力多数时间偏高的晴朗天气场景；或者a=2、b=6，表示多云天气场景。在程序中，我用MATLAB的betarnd或者betainv函数来生成对应分布样本，前者出的是随机数序列，后者配合rand可以生成指定分位数。个人更推荐直接约定a、b参数，因为后续需要解析计算半不变量时，参数化的Beta分布比纯经验分布好处理得多。

2.2 负荷波动的建模与参数选择

除了光伏出力，负荷本身也在实时变化。工程上做概率潮流时，负荷通常被建模为以预测值为期望、以一定标准差波动的正态分布，即：

P_load = μ_load + ε，ε ~ N(0, σ_load²)

需要说明的是，这里直接用正态分布有一个隐患：如果σ_load太大，抽样时可能出现负荷为负的荒谬样本。所以我在程序里做了截断处理，把低于μ_load - 3σ_load或高于μ_load + 3σ_load的样本强行截掉，不让极端值进入潮流计算。标准差一般取预测负荷的5%~20%，具体取决于系统规模和分析场景：规划分析可以取大一点，运行校核则取小一点。

负荷的有功和无功不是独立波动的。实际测量中，无功的波动往往和有功保持一定的功率因数关联，因此在建模时最好保持每个负荷节点的功率因数不变，或者给无功也设一个与有功同源的随机乘数，否则会出现“有功波动但无功纹丝不动”这种明显不符合物理的输入。这一点在蒙特卡洛抽样时不难实现，但在半不变量法中容易被忽略，后面会详细说。

3. 蒙特卡洛路线：为什么说它笨但可靠

3.1 蒙特卡洛概率潮流的总体流程

蒙特卡洛模拟的思路一句话就能讲清楚：把随机因素抽样成千上万次，每次得到一个确定的潮流解，最后把这些解的期望、方差、分布统计出来。整个流程可以拆成下面五步：

1. 设置抽样次数N，并设定输入随机变量的概率模型；
2. 当前抽样中，生成光伏出力样本和负荷波动样本，把它们叠加到节点注入功率上；
3. 对修改后的注入功率执行一次确定性潮流计算；
4. 记录本次计算的节点电压幅值、相角、支路潮流等输出量；
5. 重复步骤2到4，直到N次完成，最后对记录结果做统计分析。

这个流程几乎没有数学上的门槛，唯一的要求是你已经有一个可靠的确定性潮流求解器。我在初始版本里直接用自己的牛顿法函数，后来换到标准节点测试系统的数据格式时，也只是把数据解析部分做了适配，核心循环不变。

3.2 MATLAB实现细节：抽样、循环与随机数种子

蒙特卡洛在MATLAB里的实现，最大问题不是算法本身，而是随机数种子管理。实际写代码时，一定要在循环前用rng(固定种子)设定随机数生成器，否则每次运行结果都不同，后面验证半不变量法时连对照都没法做。

一个典型的蒙特卡洛核心循环大概长这样：

rng(2024); % 固定随机种子，确保结果可复现 N = 5000; nPV = 3; % 光伏电站数量 a = 5; b = 2; % Beta分布形状参数 V_store = zeros(30, N); Pflow_store = zeros(41, N); for k = 1:N % 抽样光伏归一化出力并换算到系统基准功率 p_pv_sample = betarnd(a, b, nPV, 1); P_pv = P_pv_rated .* p_pv_sample / Sbase; % 抽样负荷波动 loadFactor = 1 + 0.1 * randn(30, 1); % 10%标准差 P_load = P_load_base .* loadFactor; % 更新节点有功/无功注入 Pbus = P_gen_base - P_load + P_pv_inj; Qbus = Q_gen_base - Q_load; % 求解潮流 [V_k, theta_k] = myNewtonPF(Pbus, Qbus, Ybus, bus_type); V_store(:, k) = abs(V_k); Pflow_store(:, k) = computeBranchFlow(V_k, theta_k, Yf, Yt); end % 统计结果 V_mean = mean(V_store, 2); V_std = std(V_store, 0, 2);

如果安装了并行计算工具箱，可以把循环改成parfor。但要注意，parfor下每个worker的随机数序列与单核运行不一样，必须在循环体内显式设置独立的随机数流，否则并行前后结果无法对齐。更简单的做法是先在循环外生成全部所需的随机样本矩阵，再用parfor只做潮流计算。这种“样本先备好、并行只计算”的写法可复现性最好。

3.3 抽样次数怎么定

抽样次数直接决定蒙特卡洛结果的稳定性和计算耗时。理论上N取得越大越好，但实际工程中一定有个合理范围。以IEEE 30节点这样的规模，我一般取5000次，已经能稳定分辨电压期望的千分位变化；如果只关心越限概率这种尾部指标，可能要增加到10000次以上，因为尾部事件本身出现概率低，样本太少时统计波动很大。

可以做一个简单验证：把N从500逐步增加到10000，观察某个关键节点电压期望值的变化。当N增加后均值不再明显漂移，就说明样本量足够了。经验法则：误差大致按1/√N衰减，想提高一位有效精度，样本量要放大100倍，所以蒙特卡洛的精度靠堆样本，成本确实不低。

3.4 耗时优化的一个经验

刚开始写蒙特卡洛循环时，我每个样本都从统一的平坦初值启动牛顿法。后来发现，把上一个样本的潮流解作为当前样本的初始值，能明显减少牛顿迭代次数，因为相邻两个样本的注入功率变化很小，前一组的电压相角分布已经很接近新解了。算下来，每个样本平均少迭代两到三次。在小系统里这个优化看起来无关痛痒，但一旦系统规模上到几百节点、N上一万，累计节省的时间相当可观。

4. 半不变量法：把概率计算降到线性代数级

4.1 一次潮流 + 灵敏度矩阵 + 半不变量叠加

半不变量法的核心思想，是把潮流方程在基准运行点附近做一阶泰勒展开，把非线性的随机传播问题，转化为线性的灵敏度叠加问题。设系统基准潮流解为X_0和注入功率W_0，当注入功率出现随机扰动ΔW时，状态变量增量近似为：

ΔX ≈ J^(-1) · ΔW

其中J是潮流计算的雅可比矩阵。令H = J^(-1)，则H就是灵敏度矩阵，它的第i行第j列元素表示第j个注入随机变量发生单位变化时，第i个状态变量的变化量。

接下来的关键利用了半不变量的可加性。对于两个相互独立的随机变量U和V，有：

κ_k(U + V) = κ_k(U) + κ_k(V)

如果随机变量X是多个独立随机变量ΔW_j的线性组合，即X_i = Σ_j H(i,j)·ΔW_j，那么X_i的各阶半不变量可以逐项累加，但要注意线性系数的幂次：

κ_k(X_i) = Σ_j [H(i,j)]^k · κ_k(ΔW_j)

第2节提到过，光伏出力用Beta分布、负荷用正态分布建模，各自的半不变量我们是可以求出来的。把它们按上述公式叠加到状态变量上，就得到了节点电压、支路潮流等输出量的各阶半不变量。

4.2 让输入随机变量的半不变量“对齐”

半不变量法在程序实现中最容易出错的地方，是输入随机变量向量与灵敏度矩阵列的对应关系。框架上需要做一个对齐动作：把每个随机注入源按其所在的节点和类型（有功还是无功）排序，然后构造对应的灵敏度矩阵列。

以一个PQ节点接入光伏为例，光伏只向该节点注入有功，所以它的随机变量只对应注入功率向量中的P那一行；而同一个节点上的负荷如果按固定功率因数建模，负荷波动会同时引起有功和无功变化，此时应该把负荷有功和负荷无功组合成一个随机变量，而不是拆成两个独立变量。我在程序里用了这样一个简化方案：每个负荷节点只定义一个随机乘子，乘子同时作用于有功和无功两个分量，这样既符合物理规律，也让半不变量的叠加个数大幅减少。

代数表达式上，如果负荷节点j的有功扰动ΔP_j = L_j·cosφ、无功扰动ΔQ_j = L_j·sinφ，其中L_j是负荷功率波动的随机变量，那么状态变量增量的半不变量累加式应该写成：

κ_k(X_i) = Σ_j [A_ij]^k · κ_k(L_j)

其中A_ij = cosφ·H(i,P_j) + sinφ·H(i,Q_j)。这个细节不处理好，无功和有功的权重就配错了，结果会差得很离谱。

4.3 Gram-Charlier级数：从半不变量还原PDF/CDF

求出输出状态变量的前几阶半不变量之后，还要把它们还原为概率密度函数和累积分布函数。最常用的做法是Gram-Charlier级数展开：以正态分布为基础，用Hermite多项式去修正偏度、峰度带来的偏差。

将半不变量标准化，定义：

α_3 = κ_3 / σ³，α_4 = κ_4 / σ⁴

其中σ是标准差。随机变量标准化形式为z = (X - μ) / σ，则概率密度函数的四阶Gram-Charlier近似为：

f(z) ≈ φ(z) · [1 + (α_3/6)·H_3(z) + (α_4/24)·H_4(z)]

其中φ(z)是标准正态密度函数，H_3和H_4是Hermite多项式：

H_3(z) = z³ - 3z

H_4(z) = z⁴ - 6z² + 3

累积分布函数则对密度积分得到，或者直接用标准正态分布函数Φ(z)加上修正项近似。工程上一般取前四到六阶半不变量就够用了，取到更高阶时，级数修正项不仅不能保证收敛，还可能让概率密度在局部出现负值，这一点后续会详细讲。

4.4 计算量对比：为什么说差两个数量级

蒙特卡洛和半不变量法的计算资源消耗，根本不在一档。仍然以IEEE 30节点为例：

半不变量法最划算的地方在于，它把成千上万次非线性求解压缩成一次潮流求解和一次矩阵求逆。矩阵求逆的规模等于状态变量维度，对30节点系统来说就是约50乘50的矩阵，毫秒级就能完成。所以半不变量法特别适合做方案批量扫描，比如不同光伏容量、不同接入点排列组合下的风险评估，蒙特卡洛在这种场景下会算到怀疑人生。

5. IEEE 30节点系统上的MATLAB实现细节

5.1 测试系统与光伏接入方案

本文程序以IEEE 30节点标准测试系统为对象。这个系统的数据很规范：30个母线、41条支路、6台发电机，基准功率一般取100MVA。母线1是平衡节点，节点2、5、8、11、13是PV节点，其余是PQ节点。程序中我直接载入了标准测试系统的bus和branch数据格式，没有额外导入其他工具包。

光伏接入方案参考了工程上的典型做法：在三个不同的PQ节点各接一个光伏电站，容量分别取15MW、20MW、25MW，总容量60MW，相当于系统总负荷的20%左右。这个渗透率下电压波动已经比较明显，适合观察概率潮流的区分度。节点选择上也考虑了电气距离差异：一个靠近发电机区域，一个在系统中部，一个在远端负荷密集区，这样能看出不同位置光伏出力不确定性传播路径的差异。

所有光伏电站统一用Beta分布建模，取a=5、b=2，意味着归一化出力期望约为0.71，对应额定容量乘以期望后，三个电站的期望出力合计约42.6MW。

5.2 基准潮流和雅可比矩阵怎么处理

半不变量法必须用到雅可比矩阵，但直接调用现成潮流工具包时，通常拿不到矩阵本身。我在实现时没有绕弯子，直接把一个标准牛顿法潮流函数改写成了返回雅可比矩阵的版本。函数签名大致是：

function [V, theta, J] = myNewtonPF(Pbus, Qbus, Ybus, pv, pq)

函数内部先对PQ节点和PV节点编号，构造未压缩的功率不平衡方程，形成完整雅可比。返回的J矩阵维度等于未知量个数，即PQ节点数×2 + PV节点数。需要特别说明的是，这里返回的J是最终迭代收敛时的雅可比矩阵，不是初始点矩阵。因为半不变量法的线性化点就是基准潮流运行点，只有收敛点处的灵敏度才反映真实运行状态的局部动态。

5.3 灵敏度矩阵构造的关键细节

得到雅可比矩阵J后，直接求逆得到灵敏度矩阵H。但在构造H与注入随机变量的对应关系之前，必须完成两类节点剔除：

• 平衡节点不参与灵敏度计算，因为它的电压和相角是给定的，不随注入波动改变；
• PV节点的无功注入能力视为充足，光伏随机注入只加到PQ节点的有功上，不受PV节点电压控制逻辑干扰。

举例来说，30节点系统中PQ节点数目约为20个，除去平衡节点的相角方程后，状态变量总数约49个，随机注入源包括3个光伏随机变量和多个负荷随机变量。实际程序里我不会把随机注入源的维度固定死，而是设计成动态拼接，这样换系统时不用大改代码。

状态变量半不变量的累加是三层循环完成的：

kmax = 6; kappaX = zeros(nState, kmax); for i = 1:nState for j = 1:nInj for k = 1:kmax kappaX(i,k) = kappaX(i,k) + H(i,j)^k * kappaInj(j,k); end end end

外面两重循环的规模是nState乘nInj，对30节点系统来说，这个矩阵运算在一次矩阵乘法里就能完成，后面的kmax循环只是逐阶累加。真正耗时的反而是求输入随机变量的半不变量，因为Beta分布的解析矩计算涉及到Gamma函数，数值上比较繁琐。

5.4 输入分布半不变量的计算程序

如果不想费劲推导Beta分布的解析矩公式，一个更通用的办法是用“从样本估计半不变量”的数值路线。先生成大量分布样本，然后由样本的中心矩换算半不变量。

function cum = cumulantFromSample(x, order) x = x(:); mx = mean(x); m = @(r) mean((x - mx).^r); cum = zeros(order, 1); cum(1) = mx; cum(2) = m(2); % 方差 cum(3) = m(3); % 偏度相关 if order >= 4 cum(4) = m(4) - 3 * m(2)^2; end if order >= 5 cum(5) = m(5) - 10 * m(3) * m(2); end if order >= 6 cum(6) = m(6) - 15 * m(4) * m(2) - 10 * m(3)^2 + 30 * m(2)^3; end end

这个数值方法有个低调的好处：它不仅适用于Beta分布，也适用于任何你能生成样本的分布，比如带截断的正态分布。缺省情况下程序只计算前六阶，因为Gram-Charlier展开用到六阶已经足够，继续往高阶走不仅矩阵运算量增大，数值稳定性反而会变差。

5.5 结果输出：画图比看数字更有效

跑完两种方法，程序最后会输出两类图表：节点电压幅值的概率密度曲线和累积分布曲线，以及指定支路的有功潮流概率分布。我个人的实操建议是，CDF图优先于PDF图。工程决策人员看到PDF曲线往往一脸茫然，但看到CDF曲线，可以直接读出“电压低于0.95pu的概率是多少”“支路潮流超过限额的概率是多少”这类明确结论。程序里我专门封装了一个输出函数，给定半不变量后直接返回到分位数表：

function q = quantileByCumulants(cum, p) % 用Cornish-Fisher展开求百分位数 mu = cum(1); sigma = sqrt(cum(2)); z = norminv(p); q = mu + sigma * (z + (cum(3)/sigma^3)/6 * (z^2 - 1) ... + (cum(4)/sigma^4)/24 * (z^3 - 3*z) ... - (cum(3)/sigma^3)^2/36 * (2*z^3 - 5*z)); end

Cornish-Fisher展开比Gram-Charlier更适合算分位数，尤其在求95%这种尾部分位数时，不会出现概率密度为负的尴尬问题。

6. 两种方法实测对比：误差从哪里来，又该怎么排查

6.1 低渗透率场景：结果几乎一致

当光伏渗透率较低（比如总容量只有系统负荷的5%~10%）时，半不变量法和蒙特卡洛的结果非常接近。我对节点25的电压做了对照，蒙特卡洛5000次的期望电压是1.0316pu，半不变量法算出来是1.0314pu，误差不到0.02%；标准差上，蒙特卡洛是0.0023pu，半不变量法是0.0023pu，几乎没有差别。越限概率方面，因为电压偏离限值较远，两种方法都给出接近0的结果。

这说明在小波动场景下，一阶线性化近似造成的误差完全可以忽略，半不变量法完全能替代蒙特卡洛。这也是为什么很多规划工具愿意用它：几十个接入方案排下来，结果稳定且速度极快。

6.2 高光伏渗透率场景：线性化假设开始吃力

当光伏总容量提高到系统负荷的30%以上时，差距开始显现。蒙特卡洛的电压分布往往带有明显右偏特征，而半不变量法在峰度上经常估不准，导致CDF曲线在尾部出现偏差。有一次我在节点21附近接入40MW光伏并取a=2、b=6这种偏低压出力的分布，蒙特卡洛算出该节点电压高于1.05pu的概率约3.5%，半不变量法却只算出2.1%。这个误差在工程决策上已经不可忽视了。

误差来源可以拆成三个部分：

1. 一阶线性化误差：潮流方程本身非线性，注入波动大时，二阶及更高阶项不能被忽略；
2. 级数截断误差：Gram-Charlier展开只取前几阶Hermite多项式，分布偏斜越明显误差越大；
3. 输入分布假设误差：Beta分布和正态分布只是对真实随机过程的近似，如果真实分布和假设不符，任何方法都救不回来。

6.3 我踩过的三个坑与修复思路

第一个坑是Gram-Charlier展开算出概率密度为负。出现这种问题的原因通常是取到了六阶以上展开，且分布的偏度较大。我的修复方法是把展开截断到四阶，同时对偏度超过一定阈值的节点改用Cornish-Fisher求分位数，不再用Gram-Charlier生成PDF曲线。

第二个坑是灵敏度矩阵与随机注入变量错位。这个问题在程序里藏得很深：明明单个输入变量的半不变量算对了，但叠加后结果量级全乱了。排查链路是先取一个只有单一随机源的简化算例，手动验证某个状态变量的半不变量是否满足κ_k(X)=H^k·κ_k(W)，确认没有错位后，再逐步加回多个随机源。这个排查过程很枯燥，但只要怀疑灵敏度矩阵，就要坚决回头核对节点编号顺序，而不是去调后面的级数展开。

第三个坑最隐蔽：蒙特卡洛和半不变量法用的输入分布不一致。比如蒙特卡洛用Betarnd(5,2)生成样本，但半不变量法计算Beta分布矩时手写的参数却是(2,5)，两种方法结果当然对不上。后来我在程序里加了一处断言：先分别从两种方法得到的输入半不变量做一致性检验，不一致则直接报错。这算一个很值得分享的小技巧，能省掉大量排查时间。

6.4 工程配合：先快筛，后验算

基于实测对比，我给这套程序的定位是“半不变量法批量筛选、蒙特卡洛重点验算”。具体操作是，先用半不变量法扫描几十种光伏接入方案，把电压越限概率超过1%的方案挑出来；再对这些风险方案跑5000次蒙特卡洛，做最终确认。这样既避免了蒙特卡洛在全方案扫描时的巨大计算量，也不会因为半不变量法的近似误差导致风险漏判。

在IEEE 30节点系统上，这套配合流程单台普通电脑就能跑完。如果用全蒙特卡洛扫描50个方案，每个方案5000次潮流，总耗时可能是数小时；而半不变量法扫描50个方案只需要几分钟，再对其中五六个高风险方案做蒙特卡洛，总耗时控制在半小时以内，结论的可靠性也得到保证了。

最后再说一点实际体会

做概率潮流这几年，我越来越觉得程序本身不是最难的部分，难的是让分析结果真正被用起来。跑完两种方法之后，习惯别急着去看那一大堆期望和标准差，先画一条关键节点的电压越限概率随光伏渗透率变化的曲线，再画一条CDF曲线给非专业的人看，这两种图比任何均值方差表都有说服力。

更实用的建议是，写程序时一定要把输入随机变量建模和潮流求解解耦。我最初把Beta采样直接写进潮流循环里，导致后来想换负荷模型时几乎重写了一遍。现在代码里三段式结构很清楚：随机变量建模模块、确定性潮流求解模块、概率统计模块，哪个模块出了问题单独排查，互不牵连。以后再遇到馈线、微网或者含储能的系统，也能快速套用这套框架。