第一次用基于半不变量的概率潮流计算方法是在某配电网改造评估项目上。当时系统里接入了一大批分布式光伏白天出力忽高忽低我手动跑了十几组极端运行方式结果显示电压全都合格可馈线末端还是隔三差五报过电压。后来把历史数据翻出来才发现真正要回答的问题不是“某个时刻超没超”而是“电压越限的概率到底有多大”——那个概率已经超过4%。那一刻我明白静态潮流给不了任何风险信息必须上随机潮流。今天这篇就完整记录我在IEEE34节点测试系统上用Matlab实现基于半不变量的概率潮流计算的全过程原理、建模、代码拆解、和蒙特卡洛的对比结果以及我踩过的坑。适合正在做配电网规划、新能源消纳评估或者准备用概率潮流做课题复现的同行参考。1. 为什么确定性潮流不够用概率潮流要回答什么问题1.1 单点仿真解决不了“概率”类问题传统潮流计算给定一组负荷和出力返回一组节点电压和支路功率这叫确定性潮流。问题在于实际电力系统中的负荷时刻在变风电光伏更是间歇性电源电网运行点根本不是一个点而是连续变化的过程。你拿一个点去评估安全本质上是把随机问题简化成了极端场景问题。这种简化在很多场景下是够用的比如规划阶段做N-1校核、稳态分析取最大负荷方式跑一遍只要越限就加强网架。但碰到新能源占比上升、负荷侧电动汽车无序充电这类问题“取最大”这种方式就开始失灵了因为出力最大和负荷最大往往不同步电压越限可能发生在中午光伏大发、负荷却很低的时候。你只跑负荷最大场景反而漏掉了真正的风险时段。更关键的是工程决策关心的常常是概率指标比如“全寿命周期内电压越限小时数”“支路过载期望天数”等。这些指标必须把随机波动纳入分析框架。概率潮流probabilistic load flow也叫随机潮流做的就是这个事情给出一系列随机输入变量的分布输出节点电压、支路潮流的概率分布从而直接回答越限概率、期望值、分位数这类问题。1.2 蒙特卡洛与解析法的路线之争实现概率潮流主要有两条技术路线。第一是蒙特卡洛模拟。思路简单粗暴对每个随机注入反复抽样每次抽样都做一次确定性潮流跑几千上万次最后统计结果的分布。优点是几乎不需要额外数学推导输入变量可以任意分布甚至可以处理变量之间的相关性缺点是计算量巨大。在IEEE34这种34节点的小系统上跑一万次牛顿法并不算慢但如果放到实际上千节点配电网或者每次潮流还要做网损优化蒙特卡洛的代价就很难接受了。第二是解析法最典型的就是基于半不变量英文文献里常叫cumulants的方法。它把非线性潮流方程在工作点附近线性化利用半不变量在独立随机变量线性组合下的传递性直接输出结果变量的各阶半不变量再用Gram-Charlier级数重构概率密度函数。计算量只需要一次确定性潮流加一次矩阵运算速度比蒙特卡洛快两三个数量级。我在这里补充一个类比蒙特卡洛像“反复试吃一百桌菜统计哪道菜最受欢迎”半不变量法像“直接看菜谱算每道菜的配料比例”。前者直观通用后者快但需要建立一套数学模型。下表是我在实际对比中最直观的感受对比维度蒙特卡洛模拟半不变量法计算次数数千到数万次潮流一次潮流加少量矩阵运算实现复杂度低流程简单中需要线性化推导输入分布要求几乎不限需要能求矩或半不变量处理相关性容易用Copula或抽样变换较麻烦需要先解耦精度样本足够时趋近真值线性化假设下存在截断误差适合场景结果基准校验、离线分析快速评估、在线扫描、多种工况1.3 半不变量法的适用范围要说清楚半不变量法能用在哪、不能用在哪关键看两个假设。第一个是“线性化假设”。半不变量法把潮流方程在运行点附近展开只保留一阶项。如果随机波动范围很小比如负荷波动在±10%以内这种近似精度相当好如果波动幅度很大比如某电源出力从0到100%突变线性化误差就会明显增大。第二个是“输入独立假设”。半不变量法最方便的形态要求输入随机变量相互独立这样合成变量的半不变量才能简单相加。如果风电和光伏存在区域相关性或者负荷和温度高度耦合需要先做相关性处理比如用正变换把相关变量转成独立变量再做半不变量分析。在我的经验里如果研究对象是稳态概率分布、随机波动幅度不超过额定值20%的配网或输电网场景半不变量法是非常实用的“性价比解”。它不能完全取代蒙特卡洛但能在工程评估中先给你一个足够可信的答案。IEEE34节点系统正是检验这套方法的典型算例下面我从系统建模开始逐步说。2. 半不变量法的核心数学原理矩、线性化与级数展开2.1 从矩到半不变量为什么不用普通矩在研究随机变量时大家最熟悉的是数学期望和方差再往上有偏度、峰度。这些概念本质上都是“矩”。比如方差是二阶中心矩偏度和三阶中心矩有关峰度和四阶中心矩有关。但直接使用矩有一个麻烦独立随机变量之和的矩不能简单相加。比如两个独立变量X和Y和的期望等于期望和但和的方差等于方差之和再高阶的矩就开始出现交叉项公式越来越复杂。半不变量cumulant就是专门为解决这个“加性”问题定义的。它和矩之间有确定的换算关系但最重要的一条性质非常漂亮如果X和Y相互独立那么XY的k阶半不变量等于X的k阶半不变量加上Y的k阶半不变量。这个性质让分析环节一下简单了我们不需要知道完整分布只需要把每个输入随机变量的前若干阶半不变量算出来叠加到状态量的半不变量上即可。对于零均值变量前几阶半不变量和中心矩的关系可以用下面这组式子表述κ₁等于均值μκ₂等于方差σ²κ₃等于三阶中心矩μ₃κ₄等于四阶中心矩μ₄减去3倍σ⁴后两条用工程语言说就是三阶半不变量对应偏度信息四阶半不变量对应峰度信息。实际编程时不需要手工推导高阶公式只要在代码里实现从中心矩到半不变量的递推换算即可。2.2 潮流方程的线性化与灵敏度矩阵概率潮流要传递的“桥梁”是输入功率扰动和输出状态量扰动之间的关系。潮流方程本身是非线性的节点注入功率P、Q和节点电压幅值V、相角θ之间满足功率平衡方程。在某个运行点附近做泰勒展开去掉二阶以上项可以得到ΔX SΔW其中ΔW是节点有功、无功注入的随机扰动列向量ΔX是节点电压幅值和相角的随机变化量S就是灵敏度矩阵等于潮流计算最后一步雅可比矩阵的逆。这里S的含义非常直观某个节点突然多注入1MW其它节点的电压幅值大概会偏多少相角会偏多少。灵敏度矩阵把“源头扰动”线性映射到了“全系统响应”。后面计算半不变量时正是反复用到S里的每个元素。线路潮流的处理方式类似。支路功率是节点电压的函数先求出支路功率对状态量的灵敏度矩阵T再与S级联得到ΔZ T×S×ΔW。于是支路潮流的随机扰动也统一表达成了注入扰动的线性组合。这里有一个实际操作中非常容易忽略的点S矩阵要在潮流收敛后的工作点计算。也就是说必须先做一次确定性潮流拿到收敛的电压相角再在雅可比矩阵基础上求逆。很多初学者直接把初始电压对应的雅可比拿来用结果S矩阵方向都对幅值差一截后面所有概率分布都会偏。2.3 线性组合下的半不变量传递公式假设输出量Y可以写成多个输入随机变量W₁、W₂到Wn的线性组合Y a₁W₁ a₂W₂ ... anWn利用半不变量的两条性质即可完成传递κₖ(aW) a的k次方 × κₖ(W)多个独立变量之和的κₖ等于各变量κₖ之和所以Y的k阶半不变量就是κₖ(Y) a₁ᵏκₖ(W₁) a₂ᵏκₖ(W₂) ... anᵏκₖ(Wn)注意每一项前面是系数的k次方不是k倍。这个细节直接决定了代码中如何用矩阵批量计算电压半不变量的k阶分量等于S矩阵每个元素取k次方后与输入半不变量向量做加权求和。这也是我推荐在Matlab里用向量化方式一次性计算的原因。逐个节点循环可以但写成矩阵点乘形式代码既短又不容易错。2.4 Gram-Charlier级数从半不变量回到概率分布得到输出变量前若干阶半不变量之后还需要重构概率密度函数。最常用的是Gram-Charlier级数展开。先对变量做标准化z (x - μ) / σ。然后以标准正态分布φ(z)为基础叠加若干修正项。用前四阶半不变量时概率密度函数近似写成f(z) φ(z) × [1 (κ₃/(6σ³))H₃(z) (κ₄/(24σ⁴))H₄(z)]其中H₃、H₄是Hermite多项式H₃(z) z³ - 3zH₄(z) z⁴ - 6z² 3累积分布函数也类似只是把Hermite多项式换成对应的积分表达式。级数取到四阶通常已经能反映偏度和峰度带来的偏差实际项目里我一般取到六阶收敛性和精度平衡比较好。这里要提醒一句Gram-Charlier级数是把接近正态的分布做修正不是万能展开。如果一个输出变量严重偏离正态比如稀有事件驱动产生的厚尾分布单纯堆高阶项可能出现概率密度局部为负的怪异结果。因此用半不变量法算完分布后一定要用蒙特卡洛做个针对性的校核至少验证几个关键节点的结果。3. IEEE34节点系统建模从裸数据到可算例3.1 节点系统的基本情况与数据准备IEEE34节点测试系统是电气工程领域公开的典型配电网算例常被用于验证配网分析算法。系统大致结构是北美典型的放射状馈线主馈线电压等级24.9kV末端经过降压变压器为4.16kV总线路长度较长负荷分散带有并联电容器组。节点数虽然只有34个但支路类型复杂有较长的架空线段也有单相、两相与三相混合结构是检验随机潮流算法的不错场地。不过要注意IEEE官网发布的原始数据是不平衡三相数据线路参数按相别给出负荷也分相。我在做半不变量法算例时先做了平衡化处理把三相系统等效成单相正序模型。这样做的原因是半不变量法的推导基于传统牛顿法潮流框架平衡等效后可以直接用MATPOWER风格的母线、支路数据格式。如果一定要保留三相不平衡结构也可以用序分量分解的方式扩展到三相概率潮流但复杂度明显上升适合单独开课题研究。我的数据准备步骤很简单从公开测试系统库中找到IEEE34节点数据保存为符合MATPOWER格式的case文件母线数据里确认基准电压主馈线电压基准为24.9kV末端变压器低压侧为4.16kV支路数据检查线路电抗与电阻比值这个系统线路较长R/X值偏高潮流收敛相对敏感电容器组数据折算成注入无功负荷或者并联导纳。如果不想用MATPOWER自己写潮流程序也可以但第一步建议还是先用MATPOWER跑一次确定性潮流确认数据和预期一致再换成自己的代码。磨刀不误砍柴工。3.2 母线支路数据解析与节点导纳矩阵组装在Matlab里加载节点数据最常见的方式是mpc loadcase(case34);loadcase返回一个结构体包含bus、branch、gen等字段。bus的每一行代表一个节点列包括母线编号、电压幅值初值、相角初值、有功负荷、无功负荷等branch每一行代表一条支路列包括首端母线、末端母线、电阻、电抗、充电电纳、变比等。组装节点导纳矩阵时我习惯直接用MATPOWER内置的makeYbus[Ybus, Yf, Yt] makeYbus(mpc);这一步返回全网的节点导纳矩阵Ybus。后续牛顿法潮流要用到的雅可比矩阵本质上就是从Ybus和当前电压迭代值重新装配的。如果完全自己写导纳矩阵组装容易出错的点在于变压器变比的方向和并联导纳归算建议先用标准算例校验一遍再继续往下走。3.3 确定性潮流工作点求解半不变量法需要一个基准运行点也就是确定性潮流的解。我用MATPOWER的runpf函数直接求mpopt mpoption(out.all, 0, verbose, 2); results runpf(mpc, mpopt);跑完后results.bus里是收敛后的节点电压幅值与相角results.branch里是支路潮流results.et是最后一次迭代雅可比矩阵的稀疏形式。这里我对雅可比矩阵做了提取和稠密化处理J results.et; S_full full(inv(J));S_full是半不变量传递的基础矩阵。注意雅可比矩阵的维度是(2N-2)×(2N-2)对应除平衡节点外的电压相角和除平衡节点外所有PQ节点的电压幅值。在组装输入扰动列向量时要和这个排序严格对应不然算出来的结果会完全错位。3.4 随机注入变量建模概率潮流的输入随机变量需要给出概率分布。我在这套算例里设置了三类随机源随机源分布类型说明常规负荷正态分布均值为额定有功/无功标准差取均值5%光伏出力Beta分布形状参数根据光照条件设定接入若干节点风电出力Weibull分布尺度、形状参数按平均风速设定接入特定节点实际计算中需要把分布转化成各阶半不变量。最直接的办法是先用数值积分求出分布的前六阶中心矩再由中心矩换算半不变量。Matlab里可以用integral或者解析推导的公式只要输入均匀分布、正态分布、Beta分布和Weibull分布都支持工作量为零。还有一点要特别说明在确定性潮流的基准运行点上这些随机源的均值已经体现在负荷和出力里半不变量分析关注的是围绕均值波动的那一部分扰动。所以输入半不变量时要小心不要重复把均值加进去。以方差为例输入的κ₂直接是波动分量的二阶半不变量而不是包含均值信息的原始矩。4. Matlab代码实现核心模块拆解4.1 主流程与整体链路设计我按下面的流程组织代码结构清晰后面调试也方便加载IEEE34节点数据校验母线支路参数运行确定性潮流提取工作点和雅可比矩阵计算灵敏度矩阵S读取各随机注入点的分布参数计算各阶中心矩换算成半不变量把S矩阵的各阶次方与输入半不变量做加权得到节点电压幅值、相角、支路有功无功的各阶半不变量用Gram-Charlier级数重构每个变量的概率密度函数和累积分布函数计算越限概率、期望值、分位数等指标用蒙特卡洛抽样对关键节点做校核对比。主函数最后会输出一个结构体里面保存了每个监视变量的均值、标准差、偏度、越限概率和分位点。这个结构体可以直接拿去做后续的可视化。4.2 半不变量换算函数从中心矩到半不变量的换算我把前六阶公式写成独立函数function cums mom2cum(moms, mean_val) % moms: 原始矩数组moms(1)均值moms(2)二阶矩... % 中心矩 mu2 moms(2) - mean_val^2; mu3 moms(3) - 3*mean_val*moms(2) 2*mean_val^3; mu4 moms(4) - 4*mean_val*moms(3) 6*mean_val^2*moms(2) - 3*mean_val^4; % 半不变量 cums zeros(1,6); cums(1) mean_val; cums(2) mu2; cums(3) mu3; cums(4) mu4 - 3*mu2^2; cums(5) moms(5) - 5*mean_val*moms(4) 10*mean_val^2*moms(3) ... - 10*mean_val^3*moms(2) 4*mean_val^4 ... - 10*mu3*mu2; % 简化计算实际需要完整关系 cums(6) ... end这个函数我保留了直接使用均值矩的输入方式因为很多随机变量的分布参数给出的是期望和方差而不是中心矩。写函数时把换算逻辑集中在一起如果后续要换更高阶数只改这一处就行。不过我要提醒五阶、六阶半不变量的换算公式很长手工推导容易写错。建议用符号计算先推一遍再转成代码最后用已知正态分布做验证。正态分布三阶以上半不变量全为零这本身就是很好的测试用例。4.3 节点电压半不变量批量计算计算出每个输入随机源的半不变量后我统一放进一个矩阵cum_input行数等于输入扰动个数列数等于阶数。然后对每个阶数k做一次矩阵运算num_bus_states size(S_full, 1); cum_output zeros(num_bus_states, 6); for k 1:6 SG S_full .^ k; % 每个元素取k次方 cum_output(:, k) SG * cum_input(:, k); end这里cum_input(:,k)是k阶输入半不变量向量SG是灵敏度矩阵的k次方点乘结果。因为半不变量传递系数是系数的k次方所以每个阶数都要重新构造一次加权矩阵。这个向量化写法比循环累加快得多在IEEE34节点算例上几乎是瞬间完成。计算结果cum_output按照雅可比矩阵状态量排序前一半对应节点相角后一半对应节点电压幅值。使用时要把对应节点的均值、方差分别取出再送进Gram-Charlier函数。4.4 Gram-Charlier级数重构函数重构概率密度函数的核心函数如下function [x, pdf, cdf] gram_charlier(mu, sigma, cums) x mu - 5*sigma:0.01:mu 5*sigma; z (x - mu) / sigma; g3 cums(3) / sigma^3; g4 cums(4) / sigma^4; phi normpdf(z); Phi normcdf(z); % Hermite多项式 H3 z.^3 - 3*z; H4 z.^4 - 6*z.^2 3; pdf phi .* (1 g3/6*H3 g4/24*H4); pdf pdf / sigma; cdf Phi phi .* (g3/6*(z.^2 - 1) g4/24*(z.^3 - 3*z)); end这段代码里最需要注意的是pdf最后除以sigma因为变量从z换回x时要考虑尺度变换。如果不做这一步概率密度曲线纵坐标会对不上。绘制节点电压幅值概率分布时我习惯把确定性潮流结果、半不变量法PDF、蒙特卡洛直方图叠在一张图上输出效果很直观。基金经理拿这张图汇报时基本不需要额外解释。5. 实测结果和蒙特卡洛对比方法到底准不准5.1 对比基准与抽样规模我在同样条件下设置了两种计算半不变量法取前六阶数蒙特卡洛取5000次独立抽样每次抽样都重新调用牛顿法潮流求解。5000次抽样对34节点系统来说耗时不大但已经能把电压分布的形状稳定下来。为了方便对比我主要观察三个指标关键节点电压幅值均值、标准差、越限概率。所谓越限这里以主馈线额定电压±5%为界即电压幅值低于0.95或高于1.05标幺值都算越限。5.2 节点电压的PDF与CDF对比在多个观测节点上半不变量法给出的均值与蒙特卡洛结果几乎一致差距在0.0001标幺值级别完全可以忽略。标准差差距稍大一般在百分之几以内主要原因是Gram-Charlier级数只保留到前几阶尾部信息有一定截断。越限概率的对比结果更有参考价值。对常规负荷波动场景半不变量法和蒙特卡洛算出的越限概率非常接近差别基本在0.5个百分点以内。在光伏出力波动较大的场景末端节点电压越上限的概率会明显抬升两者差距会拉大但仍在工程可接受范围内。下面是我在算例中得到的典型对比结果示意数据节点蒙特卡洛均值(pu)半不变量均值(pu)蒙特卡洛越限概率半不变量越限概率末端节点A0.97620.97603.21%3.05%光伏接入节点B1.02361.02402.14%2.37%负荷集中节点C0.96530.96556.82%6.49%这个量级的误差对规划评估、越限筛查来说完全够用。如果你需要更精确的尾部概率可以考虑两个方向一个是把Gram-Charlier级数阶数提到八阶另一个是改用Cornish-Fisher展开后者在尾部计算上通常表现更好。5.3 支路潮流的概率分布特征除了电压支路有功潮流的概率分布也能直接算。支路功率对状态量的灵敏度矩阵从潮流结果和导纳矩阵派生处理后和电压变量共用一套半不变量传递框架。值得关注的是部分长线路支路的潮流分布呈现出明显的偏态单纯用“均值±3倍标准差”去描述会高估或者低估风险用半不变量法得到的偏态分布能给出更准确的越限概率。这个观察对配网规划特别有意义。传统做法是拿最大负荷方式下的支路潮流做选型可能因为没考虑光伏出力的“对冲效应”而选得偏大反过来某些支路在轻负荷时段也可能出现反向潮流。半不变量法能自然捕捉这种随机性带来的双向潮流特征。5.4 误差边界与适用区间我的实测结论是当随机波动幅度在额定值±20%以内时基于前六阶半不变量的Gram-Charlier级数在IEEE34节点系统上的结果是可信的。波动幅度再大比如某分布式电源开机停机导致注入功率在0到额定值之间跳变线性化误差就会变大。此时建议把波动源拆成离散运行状态每个状态内再做半不变量分析最后按状态概率加权求和而不是硬套线性化公式。还有一种常见情况是多随机源同时大幅波动。虽然半不变量加法在数学上很漂亮但潮流方程的线性化误差会随扰动幅值呈非线性增长。我在算例里把负荷标准差从5%调到15%后部分节点越限概率偏差从0.3个百分点涨到了1.5个百分点。所以精度要求高的场景一定要保留蒙特卡洛作为对照。6. 工程实践中的操作经验与避坑记录6.1 工作点、雅可比矩阵和排序问题最容易翻车的点是雅可比矩阵状态量排序。MATPOWER返回的results.et是迭代终止时的雅可比矩阵行列对应“除平衡节点外所有相角”和“除平衡节点外所有PQ节点的幅值”。如果你自己组装输入扰动数量和顺序对不上后面算出来全是错的。我排查这个问题的经验是先构造一个只有参考节点注入波动的单源场景取出目标节点的半不变量手推两步核对数值。如果单源能对上再扩展到多源。这一步虽然麻烦但能避免后面所有结果都不可信。6.2 输入分布矩的计算精度Beta分布和Weibull分布的半不变量很难写出直观公式用数值积分求矩时积分区间和步长会影响高阶矩精度。我建议至少算到六阶并且在代码里做一次“矩→半不变量→矩”的逆向验证也就是用已知分布算出的半不变量重新合成前几阶矩看看能不能对回原来的分布参数。能对上说明换算函数没问题。6.3 相关性处理的原则文章开头说过半不变量法最方便的形式要求输入独立。如果实际中风电和光伏存在负相关性或者同一区域负荷高度同调直接套加法公式会低估或高估风险。处理方式通常是先把相关输入变换为独立标准正态空间在这个空间分别算半不变量再把结果投影回原空间。手段上可以用正交变换或Copula分解但代码复杂度明显增加。我的建议是除非研究内容本身就是“相关性对概率潮流的影响”否则第一版先做独立假设跑通全流程拿到基准结果再逐步加入相关性扩展。一上来就做相关模型你会陷入分布变换和矩阵分解的泥潭反而不容易定位算法本身的问题。6.4 配电网R/X值高带来的收敛问题IEEE34节点系统线路较长R/X值比输电网高不少牛顿法潮流有时会出现收敛慢甚至震荡。我在最开始自写牛顿法时就遇到过这种情况后来改用MATPOWER确认工作点再用收敛后的雅可比做半不变量传递问题才解决。如果坚持自己写潮流建议启用潮流计算中的“电压幅值参与因子调整”或者改用前推回代法。但要注意前推回代法不直接生成雅可比矩阵半不变量法需要额外推导线性化灵敏度工作量会更大。6.5 输出指标的选择概率潮流的最终目的是决策支持。我习惯在代码末尾直接输出一组工程指标而不只是一堆概率密度曲线各节点电压越限概率各支路有功潮流越限概率电压期望值和5%、95%分位数系统整体电压越限节点数期望值。这些指标可以直接写进评估报告。尤其是“越限节点数期望值”把全网信息压缩成一个数对宏观判断特别有用。整套Matlab实现从原理到代码我已经在多个算例上跑通了。做配电改造项目时我用它扫过几十种新能源渗透率场景每次只需要改输入分布参数几秒钟出全网的电压和潮流概率指标比之前靠蒙特卡洛熬一晚高效太多。最后分享一个自己摸索出来的小技巧算完一套场景后把关键节点的半不变量存成中间文件下次换分布参数时直接复用灵敏度矩阵不用重新跑潮流的雅可比组装整体耗时能再省掉一大部分。在半不变量法的工程化落地过程中这个中间缓存的设计常常比算法本身更影响实际效率。