基于半不变量法的概率潮流Matlab实现与IEEE34节点验证

做配电系统或者新能源接入研究的人应该都遇到过这个问题负荷和风机出力明明都在波动可手头只有一套确定性潮流程序算出来的结果永远是一个点根本回答不了“电压越限概率有多大”“线路满载率超过80%的可能性是多少”这类问题。这几年随机潮流Probabilistic Load Flow, PLF被反复提起本质上就是把潮流计算从“单点快照”升级成“概率分布”。我这次在IEEE34节点系统上用基于半不变量Cumulant的方法做了一版概率潮流计算的Matlab实现算下来效果不错代码结构也相对清晰适合拿来当模板改。先说结论相比蒙特卡洛动辄几千次上万次重复潮流求解半不变量法只需要一次基准潮流加少量线性化计算就能把节点电压和支路潮流的概率分布给还原出来精度在工程可接受范围内速度却快了两三个数量级。这篇文章把原理、实现、验证和坑全部铺开讲研究随机潮流、做配电网分析、或者正在用Matlab写电力系统算法的人可以直接拿里面的思路和代码框架用。1. 为什么选半不变量法做概率潮流1.1 确定性潮流到概率潮流的本质转变传统潮流计算求解的是[ P_i V_i\sum_{j}V_j(G_{ij}\cos\theta_{ij}B_{ij}\sin\theta_{ij}) ]给定一组确定的负荷和发电通过牛顿-拉夫逊法迭代得到节点电压幅值和相角然后计算出支路潮流。问题在于实际系统的负荷永远在随机波动分布式光伏和风电的出力更是强随机性变量。一个确定的注入功率值本身就不存在那算出来的“确定解”意义自然有限。概率潮流把输入功率从定值改成随机变量希望得到的是输出量节点电压、支路功率的概率分布包括期望、方差、偏度甚至完整的PDF/CDF曲线。这样就能回答节点电压低于0.95 p.u.的概率是百分之几某条支路输送功率超过限额的概率有多大。对配电网规划、新能源消纳评估、静态安全分析来说这个信息密度比单个数字有用得多。1.2 三大概率潮流方法的取舍随机潮流的主流实现路线有三条我列个表对比一下方便你判断什么场景用哪种。方法核心思想优点缺点蒙特卡洛模拟法对输入随机变量反复抽样逐次做确定性潮流精度最高几乎无假设能处理任意分布和相关性计算量巨大几千次潮流才能收敛到稳定统计量半不变量法Cumulant法利用半不变量的可加性结合灵敏度矩阵线性化传递统计特性再用级数展开还原分布速度极快一次潮流加少量矩阵运算适合在线分析依赖线性化假设强非线性场景下精度会下降点估计法PEM取少量代表性采样点做确定性潮流再由点值估计输出统计矩不用线性化实现简单对非线性有一定适应对分布尾部刻画较弱高阶矩精度不稳定蒙特卡洛胜在“稳”但做参数扫描和方案比较时太慢点估计法实现简单可尾部精度不够半不变量法速度快、能给出完整分布需要付出的代价是理解门槛——你至少要弄懂半不变量和级数展开是什么。我做这个项目时综合考量后选了半不变量法主要原因是它特别适合IEEE34节点这类中等规模配电网随机扰动相对负荷水平不算极端线性化误差可控而且算出来的CDF曲线非常平滑后续做风险指标评估很方便。1.3 半不变量法核心思路一句话版半不变量法的本质是把“求两个随机变量之和的分布”这个复杂卷积问题转换成“两个变量的半不变量相加”这个简单代数问题。因为节点电压对注入功率波动的响应在基准点附近可以近似看成线性关系所以注入功率的随机性经过灵敏度矩阵传递后输出变量的半不变量可以直接由输入变量的半不变量线性组合得到最后用Gram-Charlier级数或Cornish-Fisher级数把分布函数还原出来。一句话总结输入变量的统计特性 → 半不变量 → 线性变换 → 输出变量的半不变量 → 级数展开 → 输出变量的分布函数。整个过程核心计算量只在一次牛顿法潮流和几个矩阵乘法。2. 半不变量的数学原理与级数展开2.1 矩与半不变量的关系先交代一个基础概念。随机变量 (X) 的 (k) 阶原点矩定义为 (m_k E[X^k])中心矩定义为 (\mu_k E[(X - E[X])^k])。半不变量 (\kappa_k) 是矩的另一种等价描述它的定义方式比较间接——通常用累积量生成函数 (K(t) \ln E[e^{tX}]) 在零点的泰勒展开系数来定义。实际工程里没必要每次都从生成函数积分求半不变量直接用矩到半不变量的换算公式即可。前几阶换算关系如下[ \kappa_1 m_1 \mu ][ \kappa_2 m_2 - m_1^2 \sigma^2 ][ \kappa_3 m_3 - 3m_1m_2 2m_1^3 ][ \kappa_4 m_4 - 4m_1m_3 6m_1^2m_2 - 3m_1^4 ]如果是零均值化的变量公式还能再简化但通用做法就是先算各阶矩再套换算关系。Matlab里写一个从中心矩到半不变量的通用函数也不难我后面会给出代码思路。2.2 半不变量的可加性是整个方法的地基半不变量最迷人的性质是如果随机变量 (X_1) 和 (X_2) 相互独立那么 (Y X_1 X_2) 的各阶半不变量等于两者各阶半不变量直接相加[ \kappa_k(Y) \kappa_k(X_1) \kappa_k(X_2) ]这个性质对矩是不成立的——两个随机变量之和的 (E[Y^2]) 里还带着交叉项 (2E[X_1]E[X_2])算起来要额外处理协方差。半不变量则完全免掉了交叉项独立变量求和时只要把各阶半不变量累加即可。这就是半不变量法能替代卷积计算的根本原因。在概率潮流里节点注入功率 (S_i P_i jQ_i) 通常被建模成基准值加随机扰动(P_i P_{i0} \Delta P_i)。如果各节点负荷扰动相互独立那么全网所有节点注入扰动的“总和效应”在各阶半不变量的意义上可以直接叠加再配合线性化潮流方程输出量电压幅值、支路功率的半不变量就变成输入半不变量的加权组合。2.3 从半不变量恢复PDF/CDF的两种展开有了输出变量的各阶半不变量还差最后一步怎么把半不变量还原成熟悉的概率密度函数或累积分布函数工程上常用两种级数展开。第一种是Gram-Charlier展开。核心思想是以标准正态分布为基础用输出变量的偏度、峰度等标准化半不变量去修正正态分布。设标准化的随机变量 (y (V - \mu_V)/\sigma_V)其PDF可以写成[ f(y) \varphi(y)\left[1 \frac{\kappa_3}{6}H_3(y) \frac{\kappa_4}{24}H_4(y) \dots\right] ]其中 (\varphi(y)) 是标准正态密度(H_3(y) y^3 - 3y)(H_4(y) y^4 - 6y^2 3) 是Hermite多项式。Gram-Charlier的优势是实现简单劣势是当随机变量分布偏离正态较多时截断误差会明显增大甚至出现密度函数局部为负的尴尬情况。第二种是Cornish-Fisher展开。它不是直接展开密度函数而是展开分位数。如果我们知道标准正态分布的分位数 (z_\alpha)那么随机变量 (y) 的 (\alpha) 分位数可以近似写成[ y_\alpha \approx z_\alpha \frac{\kappa_3}{6}(z_\alpha^2 - 1) \frac{\kappa_4}{24}(z_\alpha^3 - 3z_\alpha) - \frac{\kappa_3^2}{36}(2z_\alpha^3 - 5z_\alpha) ]这个展开能直接给出累计概率对应值在做电压越限概率、支路过载概率这类风险评估时非常好用。我在IEEE34节点的测试中发现Cornish-Fisher展开对偏度较大的负荷分布比Gram-Charlier更稳健尾部精度更靠谱。所以我的代码中默认用Cornish-Fisher来计算分位数和CDFGram-Charlier作为校核选项保留。3. IEEE34节点系统与输入随机变量建模3.1 系统概况与数据准备IEEE34节点测试馈线是北美配电网研究常用的标准算例原始拓扑来自一个实际运行的配电系统包含34个节点、两条主馈线和若干分支系统电压等级主要是24.9kV和4.16kV负荷分布很不均匀末端节点电压偏低是一个很典型的辐射状配网结构。原版数据还包含三相不平衡的细节如果用的是三相潮流模型处理起来会繁琐不少。我在实现概率潮流时做了一步简化先把系统转换成单相正序等值模型三相负荷累加成单相负荷线路参数取正序参数。这样做的原因是半不变量法本质上是基于交流潮流方程在基准点做灵敏度线性化单相模型已经完全能验证算法有效性再叠加三相不平衡会让灵敏度矩阵推导复杂好几倍但结论并不会发生质变。数据准备阶段最需要仔细的是节点和支路编号。IEEE34节点原始数据的节点编号从800开始到848结束中间有跳跃不是连续的1到34。我在Matlab里新建了一个映射表把原始编号映射成内部连续编号免得索引出错。这块看似简单实际容易栽跟头建议一开始就做编号映射不要偷懒直接在原始编号上操作。3.2 输入随机变量建模负荷功率和分布式电源出力是概率潮流输入随机性的两大来源。负荷建模我用的是正态分布第 (i) 个节点的有功和无功负荷分别设成 (P_i \sim N(P_{i0}, (\sigma_{Pi})^2))、(Q_i \sim N(Q_{i0}, (\sigma_{Qi})^2))基准值 (P_{i0})、(Q_{i0}) 直接取IEEE34节点数据里给定的负荷值标准差按基准值的5%~10%取。这个比例参考了电力系统负荷短期预测误差的典型水平实际工程里可以根据负荷类型调整。如果有风电或光伏出力通常不是正态分布。典型的风速服从Weibull分布风机有功出力又是风速的非线性函数所以风机出力的分布往往带有明显的偏态。对于非正态输入处理方式和正态一样——依然是先求各阶矩、再换算半不变量只不过一阶、三阶、四阶都不再是简单的那几个值。Matlab里可以直接用数值积分方法从概率密度函数求矩也可以对历史数据样本直接用矩估计。我在代码里提供了一个通用的“从样本或PDF计算前四阶矩和半不变量”的函数接入新分布时不用改主流程。3.3 输入相关性问题的取舍如果不同节点间的负荷扰动不是独立的半不变量法直接用可加性会引入误差。处理相关性有成熟套路先把相关系数矩阵做Cholesky分解再用分解矩阵把独立随机变量线性变换成具有目标相关性的随机变量然后对这个变换后的线性组合求半不变量。由于线性变换不改变半不变量的“累加结构”只需要在传递公式里多套一层变换矩阵即可。IEEE34节点的原始数据没有给出负荷间的相关系数我这次默认各节点负荷独立。但代码架构里预留了相关性处理的接口把相关系数矩阵作为参数传进去就行。实际配电网中同一馈线下的负荷受气温和用电习惯影响相关性是不可忽略的做工程评估时建议把相关系数设为0.3~0.5再算一遍看看结果差异。4. Matlab代码实现细节4.1 整体架构与模块划分我把代码拆成五个模块层次关系很明确改起来不费劲数据预处理模块load_ieee34.m读取系统数据做节点编号映射组织成潮流计算需要的结构体。确定性潮流模块run_pf.m牛顿-拉夫逊法求解基准运行点返回电压、相角、支路潮流以及雅可比矩阵。半不变量计算模块calc_cumulant.m根据输入随机变量的分布类型或样本计算各阶矩并换算成半不变量。灵敏度传递模块calc_sensitivity.m从雅可比矩阵构造灵敏度矩阵把输入半不变量线性传递给输出变量。分布还原与统计模块fit_distribution.m用Cornish-Fisher展开计算输出变量的分位数和CDF并输出期望、方差、越限概率等指标。这样划分的好处是以后换测试系统只需要改模块1换输入分布只需要改模块3换展开方法只需要改模块5主流程几乎不动。4.2 确定性潮流求解模块牛顿-拉夫逊潮流是整个算法的基础它的收敛精度直接影响后续灵敏度矩阵的质量。IEEE34节点规模不大极坐标形式的牛顿法迭代20到30次就能收敛到 (10^{-8}) 的精度。核心代码框架function [V, theta, J] run_pf(bus, branch) % bus: 节点数据 [编号, 类型, P, Q, V0, theta0, ...] % branch: 支路数据 [首端, 末端, R, X, ...] % J: 最后一次迭代的雅可比矩阵 % 初始化 V bus(:,5); theta bus(:,6); max_iter 50; tol 1e-8; for iter 1:max_iter % 计算失配功率 [dP, dQ, J] mismatch_and_jacobian(V, theta, bus, branch); if max(abs([dP; dQ])) tol break; end % 修正方程 dx J \ [-dP; -dQ]; % 更新状态 theta theta dx(1:nbus-1); V V dx(nbus:end); end end这里要注意雅可比矩阵一定要在收敛点重新计算一次并输出因为灵敏度矩阵要从收敛点的雅可比矩阵求逆用迭代中途的雅可比矩阵会引入明显误差。我第一次实现时偷懒直接用了最后一次迭代的中间矩阵结果后面的电压方差误差到了15%以上改成收敛点重算后误差立刻降下来了。4.3 注入功率半不变量计算各节点注入功率的随机扰动 (\Delta P_i)、(\Delta Q_i) 都是随机变量。对正态分布 (N(\mu, \sigma^2))前四阶半不变量非常简洁[ \kappa_1 \mu,\quad \kappa_2 \sigma^2,\quad \kappa_3 0,\quad \kappa_4 0 ]这意味着如果所有输入都是正态分布输出变量的前三阶以上半不变量全部为零分布还原时Gram-Charlier展开直接退化为正态分布——这显然是丢信息的。所以实际做概率潮流时即使输入取正态分布建议至少把负荷扰动建模成带有一定偏度的分布比如Gamma分布或者直接把三阶、四阶半不变量按照典型负荷曲线的统计特性填入非零值这样还原出的分布才有偏度和尾部信息。通用情况下我的calc_cumulant.m支持两种输入方式一种是传入概率密度函数的函数句柄用数值积分求矩另一种是直接传入历史样本用样本矩估计半不变量。核心换算代码如下function kappa cumulant_from_moments(m) % m: 前4阶原点矩向量 kappa zeros(4,1); kappa(1) m(1); kappa(2) m(2) - m(1)^2; kappa(3) m(3) - 3*m(1)*m(2) 2*m(1)^3; kappa(4) m(4) - 4*m(1)*m(3) 6*m(1)^2*m(2) - 3*m(1)^4; end用样本算矩时注意用无偏修正还是用原始矩对结果影响不大因为后级还有标准化操作但建议保持一致性要么全部用原始矩要么全部用无偏矩不要混用。4.4 输出量半不变量与分布拟合这是整个方法最核心的一步。潮流方程在基准点附近做一阶泰勒展开[ \Delta V S_0 \cdot \Delta W ]其中 (\Delta W) 是所有节点注入功率扰动组成的向量(S_0) 是灵敏度矩阵由收敛点雅可比矩阵求逆得到。对于 (Y \sum_j a_j X_j) 这种线性组合半不变量的传递公式是[ \kappa_k(Y) \sum_j a_j^k \cdot \kappa_k(X_j) ]注意是 (a_j^k) 而不是 (a_j) 的普通乘法——半不变量的阶数不同权重的幂次也不同。这是初学者最容易写错的地方。我一开始就漏掉了这个 (k) 次幂导致三阶以上的结果完全不对CDF曲线扭曲得很厉害。Matlab里实现这段传递直接用矩阵逐行操作% S0: 灵敏度矩阵size为n_out × n_in % cum_in: 输入变量的半不变量矩阵size为n_in × n_order % 对每阶k输出半不变量 S0.^k * cum_in(:, k) for k 1:4 S_power S0 .^ k; % 每个元素取k次幂 cum_out(:, k) S_power * cum_in(:, k); end得到输出变量的半不变量后就可以用Cornish-Fisher展开求分位数。CDF曲线的绘制方式是对一系列概率值 (\alpha [0.01, 0.02, ..., 0.99])用Cornish-Fisher公式求对应的分位数 (y_\alpha)描点连成曲线就是经验CDF。实测下来这个方式比直接展开密度函数再积分稳定得多。4.5 结果统计与性能对比除了绘制CDF程序最后会输出几个关键风险指标节点电压幅值低于0.95 p.u.的概率支路潮流超过额定容量80%的概率节点电压幅值的期望和标准差与蒙特卡洛基准抽样10000次相比的最大绝对误差性能数据方面我在同一台机器上对比过半不变量法总耗时约0.3秒其中基准潮流0.15秒其余是矩阵运算和Cornish-Fisher展开而10000次蒙特卡洛模拟需要跑大约3分钟。这个速度差距在配电网在线安全分析和多场景扫描场景下是决定性的。5. 实测结果与误差分析5.1 期望值与方差对比在IEEE34节点上我设置了12个负荷节点为随机注入标准差取基准值的10%用半不变量法计算得到电压幅值的期望和蒙特卡洛10000次仿真的平均值对比最大偏差出现在末端节点842附近大约0.0008 p.u.这个误差在工程上完全可以接受。电压标准差的对比误差稍大大约在3%~6%之间波动原因主要是灵敏度矩阵的线性化假设在高渗透率末端节点上偏离实际较多——末端节点重负荷时电压对注入功率的响应有轻微非线性。如果想进一步降低误差可以考虑二阶灵敏度修正即在泰勒展开中增加二阶项但计算复杂度也会明显上升。5.2 PDF与CDF曲线效果用Cornish-Fisher展开画出来的电压幅值CDF曲线与蒙特卡洛结果吻合度很高在中段5%~95%概率区间几乎重合尾部有轻微偏离。偏度较大时比如节点负荷用Weibull分布Gram-Charlier展开的尾部会出现波浪形抖动Cornish-Fisher则没有这个问题两者对比非常明显。这印证了我之前说的尾部精度优先选Cornish-Fisher。5.3 误差的来源与边界要清楚一点半不变量法和蒙特卡洛的偏差不是“算错了”而是模型的近似程度不同。蒙特卡洛把潮流方程当成黑箱反复计算保留了完整的非线性半不变量法则在基准点做了线性化截断天然牺牲了一部分非线性响应信息。负荷波动小、系统运行点接近基准点时两者差异就小如果波动大到电压偏离基准值5%以上建议用二阶近似或者干脆切到蒙特卡洛做校核。我这里给一个经验参考值当负荷标准差占基准值10%以内时半不变量法的期望误差在1%以内标准差误差在6%以内当标准差达到20%时标准差误差会上升到10%左右。负荷波动水平电压期望误差电压标准差误差建议5% 0.3% 2%放心使用10% 0.8% 6%正常使用20%约 2%约 10%结合蒙特卡洛校核6. 常见问题与排查技巧实录6.1 基准潮流不收敛概率潮流跑不出来八成是基准潮流就没收敛。IEEE34节点原始数据里部分负荷节点的有功无功比例比较极端直接套用标准牛顿法可能导致迭代发散。我的处理技巧是先用平启动flat start全部节点电压取1.0 p.u.、相角取0计算一遍如果迭代发散就把负荷按比例下调到80%算出一个中间解再以这个解为初值带全负荷重新迭代。这个方法在处理配电网重负荷算例时非常好用。另外注意PQ节点的无功上下限约束。IEEE34节点有些节点类型是PV节点或者带无功限制牛顿法中如果不处理无功越限迭代过程会出现振荡或者解出错误的运行点。我建议在每次迭代后检查无功出力越限就把节点类型切回PQ重新计算。6.2 半不变量数值溢出高阶矩的计算在数值上很不稳定。假如直接用 (E[X^k]) 的定义去算四阶矩在变量数值很大比如功率以kW为单位达到几千时四阶矩的数量级会到 (10^{12}) 以上直接带入半不变量换算公式会出现严重的浮点精度损失。解决办法有两个一是把输入随机变量做标准化后再算半不变量算完再逆标准化回去二是所有计算都用标幺值体系功率基准值取系统总容量电压取额定电压这样所有变量都在(1\times10^{-2})到(10^2)的量级内数值稳定性好非常多。我强烈建议整个程序统一用标幺值不仅算法稳定代码的通用性也好。6.3 输出CDF曲线异常扭曲如果画出的CDF曲线在中段出现S形褶皱或者尾部出现负概率多半是半不变量传递时 (k) 次幂写错了或者是级数展开截断阶数不够。先检查S0.^k这一步再检查半不变量输入是否合理尤其三阶四阶符号是否反了。如果检查无误但曲线还是扭曲就改成Cornish-Fisher展开并保留到四阶项不要用Gram-Charlier后者对偏度大的分布确实不稳定。6.4 与蒙特卡洛基准对比时对不上这里有一个很多人忽略的细节半不变量法的输入随机变量是按“节点注入功率”建模的而蒙特卡洛模拟时每次抽样都要重新求解潮流两者对负荷随机性的传播路径是一致的。但如果蒙特卡洛抽样时忽略了负荷无功和有功的联动关系比如有功按正态抽样、无功却独立抽了另一个分布两种方法的结果会系统性偏离。正确的做法是保持相关系数一致有功扰动和无功扰动之间按功率因数关联起来也就是说 (Q) 的扰动从 (P) 的扰动乘一个固定系数得到而不是独立抽样。6.5 常见问题速查表问题现象可能原因解决动作潮流迭代发散初值不当或无功越限未处理平启动/降负荷过渡/PV-PQ切换电压方差误差大于10%灵敏度矩阵用了迭代中间值在收敛点重新计算雅可比高阶半不变量数值畸变功率单位过大造成浮点溢出全程序切换标幺值体系CDF曲线中段褶皱半不变量权重缺失(k)次幂检查S0.^k传递公式与蒙特卡洛对比系统偏移有功无功抽样相关性不一致按功率因数关联抽样Cornish-Fisher尾部异常输入偏度过大或截断阶数不够增加前四阶是否保留必要时提高展开阶数到5~6阶7. 一些扩展方向和最后想说的事这个项目做完之后我顺手做了两个方向的扩展一是把输入变量从正态扩展到Weibull和Beta分布用来模拟光伏出力的随机性二是加入了节点负荷间相关系数矩阵的处理用Cholesky分解变换独立半不变量。扩展之后算法主流程一行没改只动了输入建模模块这说明按模块划分的架构确实省心。我个人经验里最值得提醒的就一件事半不变量法和蒙特卡洛不是替代关系而是互补关系。做机理分析、在线校核、多场景快速扫描时用半不变量法递进筛选对重点工况或者高精度要求的场合再用蒙特卡洛或混合法验证。把两者组合起来既保速度又保精度这才是工程上最务实的做法。如果你也在做随机潮流相关的工作建议先把IEEE34节点这个算例完整跑通再迁移到实际馈线数据上遇到的坑会少很多。