概率潮流这个东西说实话我刚接触的时候心里是有点打鼓的。做电力系统潮流分析常规做法是把风电出力和负荷当成确定值牛拉法一跑出来一组节点电压和支路功率完事了。但现实中风电和负荷是天生不确定的风大风小负荷高峰低谷这些波动会在电网里传播最终体现在节点电压越限概率、支路过载概率这些指标上。所以就有了概率潮流——把输入的不确定性考虑进去输出的不再是单一数值而是一个分布。这篇文章主要想聊聊我实际做过的这类Matlab程序基于蒙特卡洛模拟和半不变量Gram-Charlier级数展开的概率潮流计算。内容包括风电和负荷的不确定性建模、两种计算方法的原理、完整的程序架构和关键代码、以及我踩过的坑和调通之后的经验。适合正在做相关课题的电力系统研究生、刚接触概率潮流的工程师以及想找一个可靠参考实现的人。1. 为什么非要算概率潮流确定性潮流到底缺了什么1.1 确定性方法解决不了“概率”问题先摆一个最简单的场景。某个节点接了一座风电场风速今天是8m/s明天可能是12m/s风电场出力自然不同。传统潮流里你把风电出力固定成一个数比如额定容量的60%然后计算潮流。问题是这个60%是拍脑袋定的风速实际分布可能让出力在20%到95%之间波动那节点电压会不会越限支路功率会不会超载确定性潮流给不了答案因为你只算了一个点而不是整个空间。电网规划人员和调度人员需要的不是“某一种场景下电压是1.03pu”而是“电压超过1.05pu的概率有多大”。概率潮流要解决的就是这个问题给定风电场出力和负荷的概率分布求输出电气量节点电压幅值、相角、支路有功无功的概率分布或者至少求出一阶矩、二阶矩以及越限概率。这里有个很关键的认知概率潮流不是把多个确定性潮流结果简单拼起来。它需要一套系统的方法处理输入随机变量的分布、变量之间的相关性、以及非线性潮流方程对这些随机变量的传播。1.2 蒙特卡洛和半不变量法的路线选择目前主流的概率潮流方法可以粗略分成三类模拟法、解析法、近似法。蒙特卡洛属于模拟法半不变量加Gram-Charlier级数展开属于解析法还有比如点估计法算是近似法。我在这篇文章里重点讲蒙特卡洛和半不变量法因为它们是课题里最常被拿来对比验证的一对搭档。蒙特卡洛的思路非常简单粗暴根据风电场风速分布和负荷分布抽样几万次每次抽样都跑一次确定性潮流最后把几万个节点电压、支路功率结果做统计分析得到它们的概率密度函数、累积分布函数和越限概率。它的优点是不需要对潮流方程做任何线性化假设模型适应性强各种非线性环节都能塞进去缺点是计算量巨大几万次潮流计算哪怕是用牛拉法数据量一大也相当耗时。半不变量法的思路则完全相反。它基于随机变量的半不变量累积量以及潮流方程在一定精度下的线性化假设把复杂的卷积运算变换成简单的代数加法。也就是说不需要跑成千上万次潮流只需要在基准运行点做一次潮流计算得到灵敏度矩阵然后通过半不变量和Gram-Charlier级数快速重构输出量的概率分布。计算效率比蒙特卡洛高出好几个数量级代价是精度受线性化假设影响。我在实际做程序时习惯把两种方法都写了蒙特卡洛作为“标准答案”半不变量法作为“快速算法”拿前者验证后者的精度这也正是很多论文的标准做法。下面两张图程序里画的通常放在一起对比一条CDF曲线是蒙特卡洛的另一条是半不变量法的看它们重不重合。这里有一张简单的对比表方便理解两条技术路线的差异对比维度蒙特卡洛模拟半不变量法Gram-Charlier基本原理大量抽样确定性潮流线性化半不变量代数运算计算开销极高上万次潮流低一次潮流多次简单运算精度高理论上可收敛到真值中高受线性化假设限制非线性适应能力强弱适用场景验证、精细分析在线计算、快速评估2. 不确定性建模风电和负荷的“性格”要对得上2.1 风速的Weibull分布与风功率曲线做概率潮流的第一步不是写潮流程序而是先把不确定性源建模搞清楚。风电出力不确定性的源头是风速工程上最常用的风速概率模型是两参数Weibull分布概率密度函数长这样f(v) (k / λ) * (v / λ)^(k-1) * exp(-(v/λ)^k)其中v是风速k是形状参数λ是尺度参数。k一般在2到3之间取2.5左右的时候曲线形态和实际风速统计比较接近λ和平均风速相关。在Matlab里直接用wblrnd函数就能生成服从Weibull分布的随机风速样本% 风速采样 k 2.5; % 形状参数根据当地风资源统计得到 lambda 7.5; % 尺度参数单位 m/s和平均风速有关系 v wblrnd(lambda, k, 1, N_samples); % N_samples 是抽样次数拿到风速之后需要通过风功率曲线得到风电场有功出力。标准的风功率曲线是一个分段函数低于切入风速不发电在切入和额定风速之间近似按三次方关系增加达到额定风速后输出额定功率超过切出风速则停机保护。% 风功率曲线换算P_rated 为风电场额定功率 v_in 3; v_rated 12; v_out 25; P zeros(size(v)); idx (v v_in) (v v_rated); P(idx) P_rated * (v(idx).^3 - v_in^3) / (v_rated^3 - v_in^3); idx (v v_rated) (v v_out); P(idx) P_rated;这一块看起来简单但里面有个容易忽略的点实际工程里风电场还有尾流效应、风机之间的相互遮挡、风电场的功率控制策略等这些都会让实际出力和“单台风机功率曲线乘以台数”不一样。如果你做的是大电网级别的概率潮流风电场往往整体建模成一个PQ节点或者PV节点用简化功率曲线已经完全够用。如果研究的是风电场内部那要上更精细的模型那就是另一个课题了。2.2 负荷的正态分布假设负荷的不确定性标准做法是假设各节点负荷服从正态分布均值为预测值标准差取均值的1%到10%之间看负荷预测精度的实际情况。负荷预测越准标准差比例越小。对于正态分布采样Matlab里用normrnd函数% 负荷采样P_load0 为负荷母线有功预测值 sigma_ratio 0.05; % 5% 的标准差比例 P_load P_load0 P_load0 * sigma_ratio * randn(1, N_samples);这里有个细节需要注意无功负荷通常和有功负荷是耦合的工程上一般假设负荷的功率因数保持不变即无功和有功同比例波动。所以采样的时候最好是先采有功无功用功率因数乘出来而不是独立采两个正态分布。我在做程序的时候踩过这个坑刚开始独立采样有功和无功结果导致某些节点的功率因数和物理实际差得很远算出来的电压分布明显失真。2.3 变量之间的相关性不能假装它们互不认识风电和风电之间、负荷和负荷之间实际上是有相关性的。比如同一区域内的两座风电场如果地理位置接近风速会高度相关同一个区域内的负荷受气温和用电习惯影响也会同步波动。但半不变量法中随机变量的独立性假设很重要否则半不变量的可加性就不成立了。处理相关性的常用手段是通过Cholesky分解把相关的正态变量转化成不相关的标准正态变量。具体思路是先构造各随机变量之间的相关系数矩阵R对它做Cholesky分解得到下三角矩阵L使得R L * L。然后对独立标准正态采样矩阵Z做变换X L * Z得到的X就具有指定的相关系数矩阵。之后再通过等概率变换比如把正态分布的分位数映射到Weibull分布得到相关风速样本。需要说明的是这种做法在处理正态变量时比较标准处理非正态变量时需要在相关正态空间和原始分布空间之间做等概率转换过程稍复杂一些。如果只是做最简单的验证一般可以忽略相关性先假定所有输入变量独立。但在实际工程汇报中评委或者导师很容易问一句“你们的输入变量有没有考虑相关性”所以程序里最好预留相关性模块的接口。3. 蒙特卡洛模拟最“笨”但最可靠的办法3.1 采样流程和样本量怎么定蒙特卡洛的程序流程可以说没有任何智力上的难度但工程细节反而最考验人。我第一次写的时候直接把抽样、潮流计算、统计塞在一个大循环里结果跑了一个多小时才结束数据存得还乱七八糟的。后来重构了程序把流程拆成了清晰的模块。整体流程是这样的读入电网数据包括母线参数、支路参数、发电机参数、负荷参数。设置抽样次数N比如10000次。对每一次采样根据风速Weibull分布采样经风功率曲线计算风电场有功出力根据负荷正态分布采样得到各节点有功和无功更新潮流计算中的PQ节点注入功率调用牛拉法潮流计算记录该次采样的节点电压幅值、相角、支路功率循环结束后对所有记录做统计分析求均值、标准差画概率密度估计、经验CDF算越限概率。样本量怎么定理论上蒙特卡洛的收敛速度正比于1/sqrt(N)也就是说想提高一位小数精度样本量要增加到100倍。实际操作中如果只是算均值一般5000次就有比较好的效果了如果要看分布尾部比如越限概率建议至少20000次以上。我在验证半不变量法的精度时通常取20000到50000次蒙特卡洛作为基准。这里很关键的一点是不要在Matlab里用for循环写牛拉法求解器和数据记录全部挤在一起。实测下来直接在循环内反复构造稀疏矩阵并调用求解函数速度慢得让人想砸电脑。更好用的是预分配存储空间把每次潮流计算的结果存入一个预先分配好的大矩阵避免动态扩展内存。3.2 牛拉法潮流计算的时间和空间优化概率潮流里的潮流求解器本质就是确定性潮流你用Matpower的runpf也行自己写一个牛拉法也可以。但如果你自己写有几点可以优化用稀疏矩阵存储导纳矩阵Y节点导纳矩阵在电力系统里极度稀疏全矩阵存储会浪费大量内存并拖慢速度预先计算导纳矩阵并因子分解因为每次采样只是改变注入功率导纳矩阵不变所以Y矩阵和LU分解可以在循环外算好不过注意牛拉法每次迭代的雅可比矩阵是会变的这个没法完全复用。我用的是自己实现的一段基于PQ分解和牛拉法混合的求解器对于IEEE 30节点系统单次潮流大约在毫秒级别20000次采样跑下来几分钟能完成可以接受。如果你用Matpower代码会简洁很多但每次runpf内部的求解开销较大20000次可能需要近一个小时需要提前有心理准备。3.3 收敛判据和无效样本处理蒙特卡洛采样过程中有一个很容易被忽略的问题潮流计算不收敛。某些极端的风电出力和负荷组合会导致潮流根本无法收敛这在实际项目中并不罕见。处理方法不是简单的continue跳过而是要做记录。如果你跳过的样本数比例很高说明这个运行点在工程上本身就接近静态电压稳定极限这是很有价值的信息不是应该悄悄丢掉的数据。我在程序里一般设定牛拉法收敛精度为1e-8的功率偏差最大迭代次数30次。如果30次还没收敛我会把这组输入记录下来统计不收敛次数和相应的风电场出力、负荷水平。这样最后写报告的时候可以说明“在多少比例的场景下潮流不收敛对应的边界条件是什么”。4. 半不变量法把卷积变成加法的魔法4.1 先理解半不变量是什么半不变量也叫累积量统计里用希腊字母κ表示。它的好处在于有非常漂亮的可加性独立随机变量之和的半不变量等于各自半不变量之和。注意这种性质一阶矩、二阶矩也都满足但三阶中心矩不满足加法而半不变量每一阶都满足。这正是概率潮流里使用半不变量的根本原因。具体来说如果随机变量Y是若干个独立输入随机变量的线性组合比如Y a1X1 a2X2 ... an*Xn那么Y的各阶半不变量可以直接由Xi的半不变量加权求和得到不需要做卷积不需要做复杂的积分。而潮流方程在运行点附近做泰勒展开忽略二阶及以上项之后节点电压和支路功率的随机偏差恰好可以表达成输入随机变量偏差的线性组合。两个条件一拍即合这就是半不变量法能大幅简化概率潮流计算的底层逻辑。4.2 由矩求半不变量怎么算程序实现上需要从输入随机变量的分布出发求各阶半不变量。最直接的方式是先求随机变量的各阶原点矩再通过矩和半不变量之间的递推关系计算半不变量。n阶原点矩mn的定义是E[X^n]对离散样本就是均值对连续分布就是概率密度积分的展形。半不变量和原点矩之间的递推关系可以用下面的公式表达κ1 m1 κ2 m2 - m1^2 κ3 m3 - 3m1m2 2m1^3 κ4 m4 - 4m1m3 12m1^2m2 - 3m1^4 - 3m2^2 κ5 m5 - 5m1m4 20m1^2m3 - 60m1^3m2 24m1^5 - 10m2m3 30m1m2^2这些公式看起来吓人但是代码里用一个已知的递推算法就能搞定不需要手算到五阶以上。一般取到六阶半不变量就够用了再高阶的稳定性反而变差。对于常用的正态分布、Weibull分布可以直接用概率论里的公式计算矩。Matlab里可以用数值积分求矩比如用integral函数或者对理论分布直接计算% 以Weibull分布前四阶原点矩为例 k 2.5; lambda 7.5; m1 lambda * gamma(1 1/k); m2 lambda^2 * gamma(1 2/k); m3 lambda^3 * gamma(1 3/k); m4 lambda^4 * gamma(1 4/k);这里gamma是伽马函数Matlab里就是gamma函数。再代入上述递推式求出各阶半不变量。4.3 Gram-Charlier级数从半不变量到概率密度函数有了输出变量的各阶半不变量之后我们还需要还原成概率密度函数。这就是Gram-Charlier级数登场的地方。核心思想是任何一条分布曲线都可以用标准正态分布的概率密度函数φ(x)及其各阶导数来逼近。经过推导概率密度函数可以表示为f(x) φ(z) * [1 (κ3/6) * H3(z) (κ4/24) * H4(z) ...]其中z是标准化之后的变量H3、H4是Hermite多项式。实际写程序时需要的Hermite多项式H3(z) z^3 - 3z H4(z) z^4 - 6z^2 3 H5(z) z^5 - 10z^3 15z这里κ3和κ4分别是标准化后的三阶、四阶半不变量。注意κ3和κ4其实对应着统计里的偏度和峰度所以Gram-Charlier级数的本质就是在正态分布基础上做偏度修正和峰度修正。实际程序实现中标准化方式是把随机变量Y减去均值再除以标准差z (y - mean_y) / std_y;然后把上面级数加进去就得到了概率密度。累积分布函数同样可以通过对标准正态累积分布积分得到。Matlab代码不复杂% 已知输出半不变量 kappa求Gram-Charlier展开的PDF值 mu kappa(1); sigma sqrt(kappa(2)); z (x - mu) / sigma; pdf_val normpdf(z) .* (1 ... (kappa(3)/sigma^3)/6 * (z.^3 - 3*z) ... (kappa(4)/sigma^4)/24 * (z.^4 - 6*z.^2 3));就这一段代码是整篇文章的核心之一。实测下来对于节点电压幅值这种近似对称的分布Gram-Charlier展开到四阶已经能给出相当不错的拟合效果对于支路功率这种偏态比较明显的量可能需要六阶展开才能把尾部分布形态描出来。5. Matlab程序实现从零搭一个能跑的概率潮流工具5.1 程序模块怎么划分一个完整的概率潮流Matlab程序我建议按模块组织文件结构大致如下prob_pf/ ├── data/ │ ├── bus_data.m % 母线数据 │ ├── branch_data.m % 支路数据 │ └── gen_data.m % 发电机数据 ├── models/ │ ├── wind_sampling.m % 风速采样与风电出力计算 │ ├── load_sampling.m % 负荷采样 │ └── cumulant_calc.m % 半不变量计算 ├── pf/ │ ├── newton_pf.m % 牛拉法潮流 │ └── jacobian.m % 雅可比矩阵形成 ├── analysis/ │ ├── gram_charlier.m % Gram-Charlier级数还原 │ └── monte_carlo.m % 蒙特卡洛主循环 └── main.m % 主控脚本这种模块化组织不是形式主义调试的时候好处非常明显。我之前把采样和潮流写在一起出问题之后根本分不清是采样出错了还是潮流求解器报错只能到处断点调试。拆分模块之后每个环节都可以单独验证。5.2 关键代码实现思路先说半不变量法的主流程。这一段的逻辑其实是比较直的在基准运行点风电出力和负荷都取期望值计算一次确定性潮流从潮流结果中提取雅可比矩阵或者灵敏度矩阵S它描述了输入功率变化如何影响输出电压和支路功率对每个输入随机变量风电场出力、负荷计算其各阶半不变量利用灵敏度矩阵将输入半不变量线性组合成输出半不变量用Gram-Charlier级数重构输出概率密度函数。以节点注入功率到节点电压的灵敏度为例核心代码如下% 基准潮流得到 Jacobi 矩阵分解 [V0, S] base_power_flow(bus, branch, gen); % S 是灵敏度矩阵线性化系数 % 输入随机变量的半不变量矩阵C_input 维度k阶数 x 变量数 kappa_out zeros(max_order, n_bus); for k 1:max_order kappa_out(k, :) abs(S) * C_input(k, :) .* sign(mean(C_input(k, :))); end这中间灵敏度矩阵的推导是程序里最麻烦的部分。具体来说牛拉法最后一次迭代的雅可比矩阵实际上是用最后一次修正量构造的矩阵它的逆就是输入功率变化到节点电压变化的线性映射。这个细节很多资料没写清楚我当时也是翻了很久文献才搞清楚自己实现时务必注意。蒙特卡洛那部分的代码更直接核心循环如下N 20000; Vmc zeros(N, n_bus); Pflow zeros(N, n_branch); for i 1:N % 采样风速和输出 v wblrnd(lambda, k); P_wind wind_power_curve(v, P_rated); % 采样负荷 P_load P_load0 P_load0 * 0.05 * randn; % 构造该次采样的注入功率 bus(:, 3) P_load; % 有功注入 bus(:, 4) Q_load; % 无功注入 % 潮流计算 [V, success] newton_pf(bus, branch, P_wind); % 记录 Vmc(i, :) abs(V).; Pflow(i, :) branch_power_flow(V, branch); end5.3 输出结果与验证指标程序跑完之后我通常会输出三类结果第一类是概率密度对比图。选择某个风电场接入节点画蒙特卡洛直方图拟合的概率密度曲线和Gram-Charlier级数还原的概率密度曲线放在同一张图里对比。视觉上就能直观看到匹配程度。第二类是累计分布对比图。横坐标是节点电压纵坐标是累积概率。这张图的价值在于可以直接读出电压越限概率比如P(V 0.95pu)是多少P(V 1.05pu)是多少。第三类是数字指标。我习惯计算三个误差指标均值绝对误差、标准差绝对误差、CDF上最大垂直距离类似Kolmogorov-Smirnov统计量。用这几项指标量化半不变量法的精度写报告的时候有用。下面是我跑IEEE 30节点系统接入两座风电场后的典型结果仅供示意参考电压均值误差小于0.0005 pu基本可以忽略电压标准差误差1%左右CDF最大偏差通常出现在分布的左尾大约在0.02到0.05之间。这个精度在工程上是可以接受的但你必须清楚它是有代价的半不变量法基于线性化假设当风电场渗透率很高或者系统运行点接近电压稳定极限时误差会明显放大。6. 常见问题与排查技巧实录6.1 Gram-Charlier展开出现负概率密度怎么办这是用Gram-Charlier级数最大的坑没有之一。理论上概率密度函数处处非负但Gram-Charlier展开作为一种截断级数在某些区域通常是偏离均值较远的尾部会出现负值。我在测试支路功率的PDF时就遇到过概率密度在尾部变成负数的情形当时第一反应是程序写错了反复检查索引和系数后来查阅文献才明白这是级数展开本身固有的振荡问题。解决方法有几种将展开阶数从四阶降到三阶有时候低阶反而比高阶稳定改用Cornish-Fisher展开它直接作用于分位数函数天然避开了概率密度为负的问题另一种常用手段是不直接用Gram-Charlier还原概率密度而是用Edgeworth级数区别在于系数排列方式不同。实际处理上如果负概率出现在尾部远离关注区的地方也可以选择忽略只展示主要区间。但如果是出现在均值附近那说明输出分布的偏度或峰度太大线性化假设已经不太成立了这时候更建议回到蒙特卡洛结果分析原因而不是硬调级数。6.2 蒙特卡洛到底要跑多少次才够这个问题我每次都会跟同学讨论。原则是要看你想分析什么指标。算均值和标准差几千次就够想算小概率的越限事件比如0.5%的电压越限概率你需要更多样本不然统计波动太大。简单估算一下如果真实越限概率是p0.005那么N次采样中越限次数近似服从二项分布标准差是sqrt(Np(1-p))。为了让相对误差控制在20%以内需要sqrt(Np(1-p))/(Np) 0.2也就是N 25(1-p)/p大概需要5000次左右。但这只是估算。工程上我一般取20000次起步敏感场景取50000次。另外有人会用拉丁超立方采样或者重要采样来减少方差从而降低所需样本量。在概率潮流里拉丁超立方是个不错的折中方案用同样的样本量估计精度能提升一个档次。如果你论文里对计算时间有要求可以考虑在蒙特卡洛部分加一个拉丁超立方采样的选项。6.3 线性化误差在哪些场景会“爆表”半不变量法最怕非线性强的运行点。具体来说以下几种情况会使半不变量法的结果和蒙特卡洛产生明显偏差第一系统重负荷。当系统运行点接近潮流可行域边界时电压对注入功率的响应呈现出很强的非线性泰勒展开只取一阶项的假设失效Gram-Charlier还原出的分布就会严重偏离实际。第二风电场渗透率特别高。风电出力波动范围大比如从0到80%额定容量这么大的输入偏差已经远远超出“小扰动”的范围线性化自然会出现误差。第三负荷分布方差过大。如果某些节点负荷标准差设到均值的20%以上输出的电压分布尾部会被拉得很长级数展开很难完美拟合。遇到这些情况时我的实际建议是拿蒙特卡洛结果做一次对比验证。如果偏差确实大可以考虑半不变量法结合二阶敏感度分析把潮流方程的泰勒展开保留到二阶项。只是那样推导和程序实现的工作量会成倍增加适合作为进阶研究方向。6.4 程序调试的顺序和技巧写这套程序的时候我的调试顺序非常重要按照这个排查能少走弯路第一步先单独验证不确定性建模模块。把Weibull采样结果画直方图和理论概率密度曲线对照确认分布参数设置正确负荷采样同理。第二步验证半不变量计算模块。用标准正态分布做测试因为它的各阶半不变量有解析结果一阶是均值二阶是方差三阶以上全为0。如果程序输出的高阶半不变量不为0那递推公式一定有问题。第三步验证半不变量法的线性组合逻辑。用一个小型测试系统手动设置一两个随机变量手动计算期望输出再和程序结果对比。第四步最后才跑整个系统做蒙特卡洛对比。这套顺序帮我避开了无数低级错误。很多人一上来就把所有模块连起来跑结果出错后哪里都像有问题排查成本极高。最后再分享一个小技巧做了这么多概率潮流程序我最想提醒后来者的一点是不要迷信“最先进”的方法先吃透蒙特卡洛。蒙特卡洛虽然笨但它给出的结果是你判断其他一切算法好坏的标尺。我见过有人直接跳到半不变量法却连误差来源都说不清楚就是因为没有蒙特卡洛结果做基准。如果你刚开始做这个方向建议先完整地跑通蒙特卡洛把不确定性建模、潮流计算、统计分析整个链条走一遍然后再加入半不变量法模块逐步对比优化。这样虽然前期进度慢一点但你对每个环节的理解都会扎实很多。程序框架搭好之后后面无论是换更大的测试系统、加相关性模型还是扩展成计及储能和柔性负荷的概率潮流都只是在现有骨架上增加模块而已。