搞配电网规划或者分布式电源接入评估的朋友大概率都遇到过同一个尴尬你手里有一套确定的负荷数据算完潮流报告里写着“XX节点电压为0.98 p.u.”但实际上现场负荷一直在波动光伏风电更是看天吃饭这个单一的数值说白了也就是一个样本点。真正需要回答的问题是电压落在0.95到1.05之间的概率有多大线路过载的风险有多高。概率潮流Probabilistic Load Flow, PLF就是干这个的而基于半不变量Cumulant的解析法是我在IEEE34节点配电系统上反复对比过之后觉得性价比最高的一套方案。这篇文章我打算把整套思路完整拆开讲先从概率潮流的基本问题说起说明为什么我选半不变量而不是蒙特卡洛再讲IEEE34节点系统做算例的特点然后推导半不变量的核心原理给出Matlab代码框架和关键实现细节最后把我实际调试中踩过的坑整理成速查表。内容偏工程实用代码思路直接可复现适合正在做配电网规划、新能源接入评估、或者正在写相关课程设计和论文的同学。1. 概率潮流到底在算什么为什么用半不变量1.1 从确定性潮流到概率潮流传统潮流计算大家都很熟给定网络拓扑、发电机出力和负荷功率用牛顿-拉夫逊或者PQ分解法求解节点电压幅值和相角。这个计算的前提是所有输入是确定的。但实际系统里负荷随时间波动风电、光伏出力受气象条件影响电动汽车充电行为更是随机性拉满。这些不确定性会通过潮流方程传播到节点电压和支路功率上导致最终的结果也是一个随机量。概率潮流解决的就是这个问题。它把输入功率描述成随机变量经过潮流方程的“传播”输出也是带概率分布的随机变量。最终得到的不再是“电压0.98 p.u.”而是“电压落在0.95~1.05 p.u.的概率为97%越限概率为3%”。这种信息对于风险评估、规划决策来说价值极高因为你不仅知道系统在某个典型工况下能不能运行还知道它在全年负荷和天气变化下有多大概率“出事”。1.2 蒙特卡洛模拟和解析法怎么选概率潮流的主流实现思路分成两大类模拟法和解析法。模拟法的代表是蒙特卡洛模拟Monte Carlo Simulation, MCS它的思路极其简单粗暴——生成大量输入样本比如10000组符合正态分布的负荷功率然后每组都调一次潮流计算把得到的10000组潮流结果做统计分析画出电压或支路功率的频数分布就能近似得到概率密度函数。MCS的优点是准确、实现简单、几乎不受模型复杂度和非线性程度限制你甚至不需要做任何数学推导直接套个mpower的runpf循环就行。但缺点同样明显计算量太大。一个IEEE34节点的算例跑5000次潮流在普通电脑上可能就得好几分钟如果网络规模扩大到几百上千节点或者输入变量维度很高计算时间会让人崩溃。而在工程上尤其是在设计阶段需要反复调整参数、做多场景对比的时候MCS的耗时往往是不可接受的。解析法的思路则是用数学工具“推导”出输出的概率分布不需要大规模采样计算潮流。主要包括卷积法、点估计法、半不变量法、Gram-Charlier级数展开等。其中半不变量法因为计算效率极高、原理清晰、实现难度适中是我最推荐的方案。1.3 半不变量法的核心优势半不变量法为什么快核心在于它把“求随机变量和的分布”这个问题简化成了“求各阶半不变量相加”。在概率论中半不变量有一个非常漂亮的性质如果两个随机变量独立那么它们和的半不变量等于各自半不变量之和。这个线性可加性让原本复杂的多维卷积运算变成了一堆代数加法省掉了大量计算。在潮流方程中如果我们把输入功率的波动看成随机扰动那么节点电压偏离基准值的小扰动可以近似通过雅可比矩阵的逆来线性映射。这样输出的半不变量就等于输入半不变量乘以灵敏度矩阵对应元素的幂次再做加权叠加。整个过程不涉及任何循环迭代只需要一次基态潮流加上一系列矩阵运算。实测在IEEE34节点系统上半不变量法配合Gram-Charlier级数展开计算速度比蒙特卡洛快两到三个数量级而精度在常见的负荷波动范围标准差为均值10%~20%内完全够用。2. IEEE34节点系统与输入不确定性建模2.1 IEEE34节点系统有什么特点IEEE34节点测试馈线是北美配电系统比较经典的算例原始模型是三相不平衡系统包含两条主馈线、大量单相和两相支路、调压器、电容器组线路总长度比较长而且负荷沿线路分布很不均匀。相比IEEE30节点、IEEE118节点这类输电网算例它有两个对概率潮流方法影响很大的特征。第一配电网的电阻和电抗比值R/X比大线路阻抗中电阻占比高导致有功功率和无功功率对电压的影响耦合更强牛顿-拉夫逊迭代时雅可比矩阵的条件数更大。这意味着我们在做线性化时灵敏度矩阵的元素量级差异可能会很大尤其对于靠近馈线末端的节点电压对末端注入功率的灵敏度会特别高。在做半不变量法时这种高灵敏度会放大方差容易出现电压波动范围偏大的结果需要特别注意校验。第二IEEE34原始数据是三相互不独立的。很多时候我们做研究为了方便会把它整理成单相正序等效模型用Matpower的case格式直接跑。但如果你是做三相概率潮流那就需要在每一相上都建立概率模型输入变量维度翻三倍半不变量法的计算量会明显增加。这篇文章我以单相等效模型为例来演示因为这是多数论文和课程设计采用的处理方式三相扩展的原理是一致的只是维度和数据准备更繁琐。2.2 负荷随机性建模正态分布假设在概率潮流里负荷功率是最基本的随机输入。工程上最常见的简化假设是每个节点的有功和无功负荷服从正态分布均值就是基态潮流里的负荷值标准差一般是均值的10%~20%比如有一篇文章会写“负荷标准差取均值的10%”。这个假设的合理性在于大量的独立小负荷叠加后根据中心极限定理总体近似服从正态分布。需要注意的是正态分布只取正数区间因为负荷不可能为负所以实际使用时要对样本做修剪或直接用截断正态分布这个细节虽然简单但经常被人忽略。在Matlab的代码实现中生成负荷样本非常简单可以用normrnd(mu, sigma, N, 1)批量生成。但对于半不变量法正态分布有一个非常友好的数学性质它的半不变量只有前两阶非零一阶半不变量等于均值二阶半不变量等于方差三阶及以上的半不变量全为零。这意味着如果系统中的随机源全部是正态分布那输出电压的半不变量也只需要算到二阶Gram-Charlier级数展开就退化成纯正态分布。这种情况下解析法的精度会非常好但也正因为这样很多人做半不变量法时会加入非正态的分布式电源模型让方法真正发挥出处理非正态分布的优势。2.3 分布式电源与风光的随机性模型现在的配电网算例不接入点分布式电源感觉都不完整。光伏的出力模型通常可以用Beta分布来拟合某个时段内的光照强度风电的出力则一般用Weibull分布描述风速再通过风速-功率转换曲线映射到输出功率。相比于正态分布Beta分布和Weibull分布的高阶半不变量不为零这正好是半不变量法发挥价值的地方。举个例子风速v如果服从Weibull分布概率密度函数中有两个参数形状参数k和尺度参数c。根据风速-输出功率的关系风电机组的输出功率是一个分段函数在切入风速以下为0在额定风速以上为额定功率中间段近似线性。要直接用半不变量法计算最方便的做法是先用解析公式求风速的原始矩再转换出半不变量最后结合功率曲线的分段关系算出输出功率的半不变量。如果觉得这个推导麻烦工程上还有一个偷懒但有效的办法先按分布大量采样生成风电场出力样本从样本里直接估算各阶半不变量然后把这个半不变量作为输入带入解析计算。这种方法本质上是在用“局部蒙特卡洛”获取输入模型参数但后续的潮流传播仍然是解析的计算精度和效率的平衡很好。3. 半不变量法的数学原理与核心公式3.1 半不变量的定义和可加性原理先说说半不变量是什么。随机变量X的特征函数φ(t)在原点的k阶对数导数值就对应X的k阶半不变量κ_k。听起来很抽象但实际计算时根本不需要真的去求特征函数导数因为半不变量和中心矩有固定的换算关系。比如前三阶的关系是κ1 μ均值κ2 σ²方差κ3 μ3三阶中心矩。四阶及以上的换算公式会复杂一点但Matlab代码里可以直接用现成的cumulant函数计算。这个公式的本质逻辑是随机变量的各阶矩描述了分布的形状而半不变量是把这些信息重新组织成另一种形式让它具备“可加性”。假如我告诉你节点5的负荷波动是一个随机变量X节点8的负荷波动是另一个独立随机变量Y那这两个节点对潮流方程的联合影响从概率分布角度看就是XY的分布。如果你用卷积直接算需要算积分很麻烦但如果你用半不变量只需要算κ_k(XY) κ_k(X) κ_k(Y)每一阶都是加法。这一条性质直接决定了解析法的计算效率。3.2 注入功率到状态变量的灵敏度映射有了输入的半不变量下一步就是考虑它们如何通过潮流方程传播到输出。潮流方程在数学上是一组非线性方程f(V, θ) 0其中注入功率向量W是输入电压幅值V和相角θ是状态变量。严格来说输入概率分布经过非线性函数后输出分布做解析求解非常困难。但如果我们考虑的是围绕基态运行点附近的小扰动就可以把潮流方程在基态解处做泰勒展开保留一阶项忽略高阶项得到一个线性映射关系ΔX J⁻¹ · ΔW其中ΔX是状态变量增量ΔW是注入功率增量J是牛顿-拉夫逊迭代收敛时的雅可比矩阵。J⁻¹就是灵敏度矩阵S。这个线性化假设是半不变量法最核心的近似。在这个线性映射下状态变量X的半不变量和输入W的半不变量之间有一个简洁的关系。第k阶半不变量满足κ_k(ΔX) |S|^k · κ_k(ΔW)这里的|S|^k并不是指矩阵乘法而是指矩阵中每个元素都取k次幂后再乘以对应各阶的系数。严格地说对于多维输入输出节点i的k阶半不变量是所有输入节点的第k阶半不变量乘以S(i, j)^k之后求和。因为半不变量具有可加性每个输入节点独立贡献最终的结果就是累加。实现时用矩阵的逐元素幂运算Matlab里一个点乘幂就能搞定。3.3 Gram-Charlier级数从半不变量到概率密度算出了输出的各阶半不变量离概率密度函数还差一步。半不变量本身只是分布的数字特征不是分布本身。要从有限阶半不变量恢复概率密度函数目前主流做法是用Gram-Charlier级数展开。Gram-Charlier级数的基本思路是以标准正态分布为基准函数用一个多项式乘以标准正态密度函数来逼近真实的概率密度。具体来说先将状态变量X标准化令z (X - κ1)/√κ2那么标准化的概率密度函数可以展开为f(z) φ(z) c1·φ(z) c2·φ(z) c3·φ(z) ...其中φ(z)是标准正态分布的密度函数φ(z)、φ(z)是它的导数系数c_i由标准化后的半不变量决定。实际使用中通常截断到4阶或6阶因为阶数太高反而会出现数值振荡尾部可能发散。前几阶系数的物理意义很明确三阶半不变量对应偏度四阶对应峰度它们分别描述了分布相对于正态分布的扭曲程度和尖峰程度。Matlab里实现Gram-Charlier展开时比较方便的做法是利用标准正态分布密度函数的导数与Hermite多项式之间的关系。你不需要手算高阶导数公式可以直接调用Hermite多项式函数利用递推关系自己写一个然后逐项累加。最终得到的f(z)就是标准化的概率密度函数还原成原变量X的密度函数时还要做一个变量替换把标准化时的缩放系数乘回去。4. Matlab实操全过程4.1 数据准备与运行环境我通常用Matpower作为基态潮流的求解工具。Matpower自带多个IEEE标准算例数据但需要注意IEEE34节点原始数据是三相不平衡的配电网模型Matpower自带的case格式并不直接包含这个算例需要自行转换或者从第三方下载整理好的单相等效版本。我一般是把IEEE34的母线数据、支路数据、变压器数据整理成Matpower的case结构然后用runpf函数跑一遍基态潮流确认结果收敛、节点电压在合理范围内。数据准备阶段最容易出的问题有两个一是母线编号的连续性Matpower要求母线编号是1到34连续整数但IEEE34原始数据里母线编号本来就是1到34基本不用调整二是支路参数的单位Matpower统一用p.u.而IEEE34原始数据很多给的是欧姆和千伏需要先算出统一的基准容量和基准电压然后做归算。这里我习惯把基准容量设为1 MVA基准电压取馈线标称电压24.9 kV然后把每条支路的阻抗值先归算到该基准下。4.2 半不变量计算的核心代码框架基态潮流跑完之后就可以开始半不变量法的核心计算了。整个流程分为四步算输入半不变量、算灵敏度矩阵、算输出半不变量、Gram-Charlier展开。先看第一步构造输入随机变量的半不变量矩阵。假设系统有nb个节点每阶半不变量是一个nb维向量。% 半不变量阶数 maxOrder 6; % 每个节点的负荷标准差比例 sigmaRatio 0.1; % 初始化半不变量矩阵 K_input(order, nb) K_input zeros(maxOrder, nb); for i 1:nb baseLoad bus(i, PD); % 基态负荷有功 mu baseLoad; sigma sigmaRatio * baseLoad; % 正态分布半不变量一阶为均值二阶为方差 K_input(1, i) mu; K_input(2, i) sigma^2; % 三阶及以上为0 end如果还要加风电在节点k那就在K_input矩阵的第k列加上风机出力对应的半不变量。风机出力的半不变量可以通过采样估算也可以用解析公式。下面给出一个用样本估算半不变量的函数function cum sampleCumulant(data, maxOrder) % data是服从该分布的大量样本 % 先求中心矩再换算成半不变量 m1 mean(data); m2 mean((data - m1).^2); m3 mean((data - m1).^3); m4 mean((data - m1).^4); cum zeros(1, maxOrder); cum(1) m1; cum(2) m2; cum(3) m3; cum(4) m4 - 3*m2^2; % 更高阶可依照中心矩-半不变量公式换算 end第二步是灵敏度矩阵。在Matpower中跑完基态潮流后可以用牛顿-拉夫逊迭代的最终雅可比矩阵来构造% 算例基态求解 mpc loadcase(case34); result runpf(mpc); % 提取雅可比矩阵Matpower的runpf在result里并没有直接给出J % 常见做法是用nargout2的牛顿法接口或者自己重新组装 % 这里演示直接调mpower内部函数获取 [J, ~] opf_jacobian(result); % 如果版本支持 % 如果上面不行用pf的函数接口 [V, converged, iterations, J] newtonpf(result.Ybus, ...);注意Matpower不同版本获取雅可比矩阵的方式不同。有些版本里runpf的结果结构体会有internal字段里面存放了雅可比矩阵。如果实在拿不到有一个笨办法自己写一个简单的牛顿-拉夫逊潮流顺手把迭代收敛时的雅可比矩阵输出。这在IEEE34节点这种规模下完全可行代码量也不大。第三步用输入半不变量和灵敏度矩阵计算输出半不变量S inv(J); % 灵敏度矩阵 K_output zeros(maxOrder, 2*nb); % 每个节点有幅值和相角两个量 for k 1:maxOrder % 对每个节点求和 S(i,j)^k * K_input(k,j) Sk S.^k; K_output(k, :) Sk * K_input(k, :); end这里需要注意一点雅可比矩阵对应的状态变量是电压幅值和相角所以J的行数和状态变量数一致是2*nb-2平衡节点相角除外。如果直接把J求逆得到的S矩阵维度和输入向量不一致要么在输入向量里也把平衡节点对应的元素排除掉要么在结果里补上平衡节点的零上下波动。实际代码里我习惯把平衡节点相角波动统一设为0电压幅值波动也设为0在矩阵运算时单独处理避免维度混乱。第四步Gram-Charlier展开。这里给出一个常用的实现用Hermite多项式递推得到f(z)function pdf gramCharlier(z, cum) % z为标准化后的变量取值 % cum为输出的1到maxOrder阶半不变量 sigma sqrt(cum(2)); % 标准化半不变量 s3 cum(3) / sigma^3; s4 cum(4) / sigma^4; if length(cum) 5 s5 cum(5) / sigma^5; else s5 0; end if length(cum) 6 s6 cum(6) / sigma^6; else s6 0; end pdf (1/sigma) * exp(-z.^2/2) / sqrt(2*pi) .* (1 ... s3/6 * hermitepoly(3, z) ... s4/24 * hermitepoly(4, z) ... s5/120 * hermitepoly(5, z) ... s6/720 * hermitepoly(6, z)); end function H hermitepoly(n, x) % Hermite多项式递推计算 H0 ones(size(x)); H1 x; if n 0 H H0; elseif n 1 H H1; else for k 2:n H x .* H1 - (k-1) * H0; H0 H1; H1 H; end end end到这里你已经有办法在任意z值处计算概率密度了。画图时把电压范围从最小值到最大值划分成若干点依次求PDF值然后plot就行。累积分布函数CDF可以对PDF做积分得到或者在输出半不变量基础上用Gram-Charlier的累积分布展开公式直接算。4.3 结果对比与精度检验代码写完必须验证。我验证的标准方法是把半不变量法的结果和蒙特卡洛模拟做对比。以IEEE34节点中某个末端节点比如节点848的电压幅值为例蒙特卡洛做5000次采样得到电压的经验分布半不变量法用Gram-Charlier级数展开得到解析的PDF和CDF。把两条曲线画在同一张图上观察趋势是否一致。为了定量评估误差可以用两个指标一是期望值的相对误差二是方差的相对误差。两者的计算公式分别是误差_e |μ_ML - μ_CU| / μ_ML * 100% 误差_var |σ²_ML - σ²_CU| / σ²_ML * 100%在我实测的IEEE34算例中负荷标准差取均值的10%负荷全部为独立正态分布时期望值误差在0.1%以内方差误差在3%以内这个精度对于工程评估来说已经足够好。如果加入Weibull分布的风电模型高阶级数的影响会更明显Gram-Charlier截断到6阶时偏度和峰度特征也能较好地刻画出来误差会略大一些但趋势依然正确。精度校验的另一个实用技巧是画“概率密度曲线叠加频数直方图”。把蒙特卡洛模拟得到的5000个电压幅值画成直方图再把半不变量法求出的PDF曲线画在上面如果曲线和直方图的包络吻合说明解析法很好地捕捉了分布的形态。这样可以直观地检查出是否存在截断阶数不足、尾部振荡等问题。5. 实战中常见的坑与排查方法5.1 IEEE34数据与Matpower接口的坑IEEE34原始数据是配电网格式里面的变压器是三相三线制还带有调压器站直接塞进Matpower的case35格式里会有几个绕不过去的问题。首先是变比和阻抗基准。IEEE34给的变压器阻抗是标幺值但基准容量和基准电压和Matpower默认不同直接填进去会导致潮流结果异常最常见的就是电压全超限、甚至迭代不收敛。解决办法是先仔细阅读用例数据说明把变压器支路的阻抗归算到统一基准下再填入buses和branches矩阵。其次是调压器。IEEE34系统里有两处调压器它们的作用是自动调整变压器变比来维持下游电压。Matpower的branch结构本身也带变比字段但不会自动模拟调压器动作。所以在概率潮流里我一般采取两种策略要么把调压器固定在中性档位不做调节要么在潮流迭代中增加一个变比调整的逻辑。如果你做的是静态概率潮流建议直接固定变比因为半不变量法只围绕一个基态运行点做线性化变比变化本身会改变雅可比矩阵结构处理起来太复杂。5.2 半不变量法的适用边界半不变量法不是万能的最大的限制来自线性化假设。当输入功率波动范围太大时比如某节点负荷标准差达到均值的50%以上或者风电接入容量占系统总负荷比例很大线性化的误差会显著增大基于一阶泰勒展开的灵敏度映射就会失真概率密度曲线可能和蒙特卡洛结果出现明显偏差。遇到这种场景我的建议是先调整基态运行点把大扰动看成多个小扰动叠加或者引入二阶项修正。但引入二阶项会让代码复杂度上一个台阶工程上如果只是初步评估不如直接退回到蒙特卡洛模拟更省心。其次半不变量法要求输入变量独立。如果节点负荷之间存在显著相关性简单的半不变量相加就不成立了。处理相关性的方案是做Cholesky分解或Copula变换先把相关的输入变量转换成独立的标准正态变量再代入半不变量计算。这个扩展我在文章里不展开但你需要知道它存在的必要性。5.3 实际调试过程实录我在调试这个算例时遇到过一个很典型的bugGram-Charlier展开后的概率密度函数在某个区域出现了负值。概率密度为负这在数学上当然不合理但出现的原因并不难理解因为Gram-Charlier级数是一个带截断的展开在分布的尾部截断误差可能导致负的贡献值。查到问题后我换了两种处理方式。一是增大阶数从4阶截断换到6阶截断尾部振荡明显改善二是利用CDF单调不减的特性先算CDF再对CDF做数值微分得到PDF这样即使CDF本身有微小振荡PDF也不会出现负的极端值。这两种方法配合使用后曲线光滑度改善很多。遇到概率密度为负的情况千万别直接认为半不变量法失效多数时候只是展开阶数和数值稳定性问题。另一个调试坑是雅可比矩阵元素量级差异过大导致的数值灵敏度问题。在配电网中相角对有功灵敏度往往是10⁻²量级电压幅值对无功灵敏度可能是10⁻¹量级两者相差悬殊。在计算高阶半不变量时S^k会对这些差异做指数级放大导致线性方程组出现数值病态。我的处理办法是在计算时把灵敏度矩阵按列归一化算完结果再还原或者直接把某些灵敏度特别大的元素做岭值截断防止某个节点单独主导整个结果。6. 概率潮流结果的工程应用6.1 配电网规划中的电压越限风险评估概率潮流做出来之后能干什么最直接的应用就是电压越限风险评估。在配电网规划阶段你要判断新增一条线路或者一个分布式电源接入点后全年各节点电压会不会有越限风险。用半不变量法可以一次性算出每个节点电压的完整概率分布然后通过CDF快速读取越限概率。比如IEEE34节点系统中节点848的电压一旦低于0.95 p.u.就认为是低电压风险事件。算完累积分布函数F(v)在v0.95处求值如果F(0.95)0.03意味着有3%的概率电压越下限。这个概率如果超过了规划允许的阈值比如2%那就需要加强线路或者调整调压器设定。这种定量评估比单纯看基态潮流结果要可靠得多因为基态只能告诉你典型工况不越限但无法回答全年随机场景下的风险有多大。6.2 分布式电源接入容量评估另一个典型应用场景是分布式电源的接入容量评估。假设要在IEEE34节点的某个末端节点接入一座光伏电站接入容量多大合适传统做法是不断调大光伏容量运行潮流看电压是否越限。这种确定性分析方法给出的答案是边界值比如最大接入5 MW但实际运行时光伏出力波动可能瞬间突破这个边界导致电压越限。用概率潮流的话你可以把光伏出力建模为Beta分布然后扫描不同接入容量每个容量下都算电压越限概率。最终得到一条“接入容量-越限概率”的曲线规划人员可以根据概率风险承受能力去选择合适容量。这种思路在现在新能源大量接入配电网的背景下尤其有参考价值。6.3 方法扩展方向半不变量法本身也有不少扩展空间。比如除了电压和支路功率还可以把网损的概率分布也算出来在输入相关性建模基础上还能扩展到多维相关的负荷场景如果结合时序数据和负荷预测的误差分布还能做短时间尺度的概率安全评估。对于正在做相关研究的朋友我建议在掌握本文基础代码后优先往相关性处理和二阶线性化两个方向深入这两个方向都是论文和实际工程比较认可的价值点。回到代码本身我个人在实际使用中有一个习惯始终保留一套蒙特卡洛模拟的低采样版本比如500次的MCS作为快速校验基准每改动一次解析法代码就跑一遍对比。这个习惯虽然多花一点时间但能极快地暴露线性化假设、数据接口和数值实现上隐藏的问题。做概率潮流这种计算数学推导错了你很难从结果里一眼看出来但和基准结果一对比偏差立马现形。这篇文章提到的代码框架和坑基本都是从这套对比流程里反复磨出来的直接拿去用能帮你少走不少弯路。