如果你平时只做确定性潮流计算第一次听到“概率潮流计算”这个概念脑子里大概率会冒出一串问题潮流结果还能有概率分布多跑几次仿真不就行了半不变量又是什么东西跟蒙特卡洛比到底强在哪这套基于半不变量的随机潮流方法本质上是把“输入随机性”通过潮流方程映射成“输出随机性”的一种解析近似手段计算速度快、占用资源少特别适合用在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的思路是取输入变量的若干特征点每个特征点跑一次潮流然后用加权组合还原输出矩信息。计算量是2n1次潮流n为随机输入维数比蒙特卡洛快很多但精度对特征点选取敏感。半不变量法走的是另一条路不做大量采样而是把输入随机变量的各阶半不变量求出来通过潮流方程的线性化灵敏度映射到输出变量再用级数展开重构概率密度函数和累积分布函数。整个过程只需要一次基态潮流加几次矩阵运算毫秒级出结果。下表是三种方法的对比。方法计算代价精度特征适用场景蒙特卡洛数千到数万次潮流精度随样本数提升无线性化假设基准校验、高精度分析点估计法2n1次潮流前三阶矩精度较好中等规模快速评估半不变量法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⁻¹ · ΔSJ⁻¹就是灵敏度矩阵。它把输入随机变量的随机性“传导”到输出状态变量上。对于第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:nPVnPQ); V V dX(nPVnPQ1: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里实现这套算法建议按文中的模块顺序一步步搭建每完成一个模块就用蒙特卡洛交叉验证一下这样出了问题时能很快定位到具体环节。