风电、光伏大规模接入之后电网分析的默认假设就变了。以前做潮流计算给一组确定的发电机出力和负荷牛拉法迭代一遍得到一组节点电压和支路功率这套流程在传统火电主导的时代够用但当出力看天吃饭的新能源占比上来输入侧的确定性假设本身就站不住脚。概率潮流计算Probabilistic Load FlowPLF正是用来应对这种不确定性的工具——它把风速、光照、负荷都建模成随机变量通过大规模采样或解析近似得到节点电压和支路功率的概率分布而不是一个孤零零的点值。Matlab凭借矩阵计算优势和自带的统计工具箱一直是做概率潮流最顺手的平台。这篇文章从工程应用的角度把含风光发电的概率潮流计算的数学模型、三种主流算法、Matlab代码框架和调试经验完整过一遍刚接触这个方向的研究生和做新能源接入评估的工程师都能直接拿去参考。1. 概率潮流到底在算什么从单点结果到概率分布的思维转换1.1 确定性潮流的“天花板”在哪里传统的确定性潮流求解的是这样一组方程节点注入功率等于电压与导纳矩阵的乘积方程组给定以后迭代求解出唯一一组节点电压和相角。问题的关键在于“给定”这两个字。系统里有风电和光伏之后注入功率不再是一个确定数值。同一座风电场年平均风速差一点全年发电量可能差出好几个百分点同一天内云层飘过光伏出力能从额定功率跌到零。此时如果你还用单点出力去做潮流分析得到的结果只对应某一种天气场景对系统规划人员来说参考价值非常有限。更麻烦的是电压越限这类风险恰恰容易出现在极端场景里大风天气下风电场满发而本地负荷处于低谷多余功率外送导致局部电压偏高傍晚光伏快速退出负荷却还在高位电压又会往下掉。确定性潮流算出来的“正常工况”电压往往看不出这些边界风险概率潮流正是把这类风险显式地量化出来。1.2 概率潮流的输入输出形态概率潮流做的事情本质上是输入的随机性到输出的随机性的传递。输入侧主要包括三类随机变量风速一般用两参数Weibull分布描述概率密度函数为 f(v) (k/c)(v/c)^(k-1) exp(-(v/c)^k)其中 c 是尺度参数k 是形状参数。光照强度通常用Beta分布描述因为光照强度在0到额定值之间连续变化Beta分布定义在有界区间上形态灵活。负荷功率一般用正态分布或对数正态分布波动范围取均值的3%~8%比较常见。输出侧就是电网运行人员真正关心的东西节点电压幅值的均值、标准差、概率密度函数和累计分布函数支路有功和无功的分布电压越上/下限概率支路过载概率系统网损的期望值。这些指标可以支撑三个层面的决策——规划阶段评估新能源接入容量是否过度激进运行阶段判断当前方式下电压风险水平调度阶段为备用容量和AGC调节留出合理区间。1.3 概率潮流与常规随机分析的差别有读者可能会问这和蒙特卡洛模拟直接撒点有什么区别蒙特卡洛只是概率潮流的实现手段之一。概率潮流本身是一套完整的方法论框架包含输入随机建模、相关性处理、不确定传递、输出统计四个环节。蒙特卡洛是最直观的传递手段点估计法和半不变量法则走解析路线用更少的计算量逼近同样的统计信息。后面会把这三种方法的Matlab实现逐一展开。2. 先把数学模型搭起来风光出力的随机特性描述2.1 风速与风电出力模型风速模型是整个概率潮流里最容易出错的环节因为风速分布和风机出力之间还隔着一道非线性分段函数。Matlab里生成Weibull随机数直接调用 wblrnd(scale, shape) 即可但这里有个经典坑位Matlab的wblrnd第一个参数是尺度参数c单位m/s第二个是形状参数k无量纲我见过不止一次有人把参数填反结果生成的风速样本整体偏小或偏大风电出力分布完全失真。风速到电功率转换用最通用的分段线性模型P_w(v) 0, v v_in 或 v v_out P_w(v) P_rated * (v - v_in)/(v_rated - v_in), v_in v v_rated P_w(v) P_rated, v_rated v v_outv_in 一般取3m/sv_rated 取11~13m/sv_out 取25m/s。这个模型虽然简单但在工程分析里足够用。更精确的风功率曲线可以用厂商实测数据插值不过概率潮流关注的是长期统计分布分段线性模型引入的误差通常可以接受。Weibull参数的获取有两种途径一是用当地测风塔一整年的小时级平均风速数据通过极大似然估计反推二是用平均风速和标准差近似估算。后者在参考资料匮乏的时候很实用c 约等于平均风速的1.12倍k 则通过变异系数查表或数值求解。估算出来的参数用于方案预研足够正式工程评估还是建议用实测数据。2.2 光照强度与光伏出力模型光伏出力的建模同样分两步。第一步用Beta分布描述光照强度的随机性Beta分布的概率密度函数为 f(s) (Γ(ab))/(Γ(a)Γ(b)) * s^(a-1) * (1-s)^(b-1)Matlab里用 betarnd(a, b) 直接生成0到1之间的标幺值。a和b的取值决定了分布的偏斜程度夏季晴天多曲线右偏a取2~3、b取1~1.5比较贴合阴雨多的地方分布相对左偏参数相应调整。第二步做光电转换。理想情况下光伏出力与光照强度近似线性P_pv η * S * A其中η是光电转换效率S是实际光照强度A是光伏阵列面积。实际工程里更常用的是容量标幺法P_pv P_rated * (S / S_ref)S_ref通常取1000W/m²。温度对光伏出力的影响也不能完全忽略组件温度升高导致输出电压下降更精细的模型会在上述公式基础上乘一个温度修正系数(1 - β(T - 25))β一般取0.003~0.005 /°C。对于概率潮流温度修正要不要做得这么细取决于计算目的——如果只是评估年度电压分布简化模型足矣如果要做夏季高温极端场景分析建议把温度项加进来。2.3 负荷随机性与相关性处理负荷波动用正态分布描述最普遍但需要加截断处理。我在Matlab里习惯写成PD_load PD_mean .* (1 0.05 * randn(n, 1)); PD_load(PD_load 0) 0;不做截断的话理论上会生成负负荷样本虽然概率极低但一旦出现就会导致潮流方程出现负注入结果根本没法解释。0.05这个系数对应5%的标准差如果节点负荷本身波动大可以放宽到8%~10%。另一个容易被忽视的问题是相关性。风电场和光伏电站如果处在同一个区域电网风速和光照之间可能存在相关性不同节点的负荷也不会完全独立。忽略相关性的后果是低估系统电压波动的风险区间。Matlab里处理相关高斯随机变量最直接的方式是用Cholesky分解给定相关系数矩阵R计算 L chol(R)然后把独立标准正态样本矩阵 Z 乘以 L 的转置得到带相关性的样本。对于Weibull和Beta这类非高斯变量更严格的做法是引入Copula但工程上如果只是粗略评估Cholesky分解后做概率积分变换也能用。3. 核心方法一蒙特卡洛模拟最稳也最烧钱3.1 三步走采样、潮流计算、统计蒙特卡洛模拟的思路直白到不需要过多解释既然输入是随机变量那就生成大量输入样本逐个做确定性潮流计算最后把所有输出结果汇集起来做统计分析。三个步骤对应三段Matlab代码每一步都可以调优。第一步是采样。样本量太小分布形状出不来样本量太大计算耗时线性增长。第二步是潮流计算每次都调用一次确定性潮流求解器。第三步是统计mean、std、histogram、prctile几个函数就够用。蒙特卡洛最大的优点是“无偏”——只要样本量足够大输出分布一定收敛到真实分布不依赖任何线性化假设。这点在后两种解析方法里是做不到的。代价就是计算量大。如果单次牛顿法潮流耗时0.01秒5000次就是50秒在多节点系统里单次潮流可能到0.1秒5000次就是8分钟级别。3.2 Matlab主程序框架建议用Matpower做潮流计算内核自己手写牛顿法在原理上没问题但处理PV-PQ节点转换、无功越限这些细节时容易踩坑。Matpower自带IEEE 14节点等标准算例loadcase一键加载runpf闭环求解省心很多。下面这套代码是我的常用框架对应风电和光伏分别接入两个不同节点的情况% 蒙特卡洛概率潮流主框架Matpower 统计工具箱 mpc loadcase(case14); base_PD mpc.bus(:, 3); % 保存原始有功负荷 base_QD mpc.bus(:, 4); % 保存原始无功负荷 N 5000; Vr zeros(N, 14); % 电压幅值记录矩阵 Pa zeros(N, 41); % 支路有功记录矩阵case14共41条支路 Sr zeros(N, 14); % 节点注入视在功率记录 % 风速Weibull参数wblrnd(尺度c, 形状k)顺序不要写反 c_w 8.5; k_w 2.2; v_in 3; v_rated 12; v_out 25; Pw_rated 30; % 风电场额定功率 MW % 光照Beta分布参数 a_s 2; b_s 1.5; Pp_rated 15; % 光伏电站额定功率 MW for i 1:N % 1. 采样风速 - 风电出力 v wblrnd(c_w, k_w); if v v_in || v v_out Pw 0; elseif v v_rated Pw Pw_rated * (v - v_in) / (v_rated - v_in); else Pw Pw_rated; end % 2. 采样光照强度 - 光伏出力 s_pu betarnd(a_s, b_s); Pp Pp_rated * s_pu; % 3. 把风光出力等效为相应节点的负负荷注入 PD_node base_PD; PD_node(9) base_PD(9) - Pw / 100; % 风电接入bus9 PD_node(13) base_PD(13) - Pp / 100; % 光伏接入bus13 % 4. 负荷波动对包含负注入的净负荷施加正态扰动 PD_load PD_node .* (1 0.05 * randn(14, 1)); QD_load base_QD .* (1 0.05 * randn(14, 1)); PD_load(PD_load 0) 0; QD_load(QD_load 0) 0; mpc.bus(:, 3) PD_load; mpc.bus(:, 4) QD_load; % 5. 确定性潮流计算 res runpf(mpc, mpoption(out.all, 0)); if res.success 0 warning(第 %d 次潮流不收敛, i); continue; end % 6. 记录输出量 Vr(i, :) res.bus(:, 8); % 电压幅值 Pa(i, :) res.branch(:, 14); % 支路有功 end注意代码第3步和第4步的顺序。先把风电和光伏出力折算成负的净负荷再对这种净负荷施加正态扰动等价于“风光出力随机 负荷随机”的叠加方式。如果你先把原负荷扰动完再把风光注入单独减掉那么在注入很大的节点上净负荷可能出现负值且分布形态被扭歪。另外代码里 Pw/100 是因为case14的基准容量是100MVA把所有功率统一折算到标幺值这一点新手特别容易漏。3.3 收敛性判断与采样规模选择到底采多少组样本才算够经验法则是先看输出量的均值或标准差随样本数的变化曲线当相对波动小于某个阈值时认为收敛。我常用的判据是电压均值的无穷范数误差连续两次采样之间的变化量小于1e-4就停止。更直接的做法是固定样本量5000次起步不好再翻倍到10000次。N_max 10000; V_mean_old zeros(1, 14); for i 1:N_max % …… 上述采样与潮流计算代码 …… V_mean_new mean(Vr(1:i, :), 1); delta max(abs(V_mean_new - V_mean_old)); if delta 1e-4 i 1000 fprintf(均值收敛于第 %d 次采样\n, i); break; end V_mean_old V_mean_new; end这个小循环跑起来有个好处你不用赌样本量机器自己告诉你够了。代价是循环内多了mean运算对整体耗时影响很小。蒙特卡洛还天然支持并行化把for改成parfor前提是循环体内不能有依赖全局变量的操作runpf的mpc结构每次都是基于本次采样数据构造的满足parfor要求。我实测过在8核机器上3000次仿真的耗时能压到原来的三分之一左右。注意parfor里mpc这个变量会被当作广播变量处理数据量不大影响有限。4. 省时方案点估计法与半不变量法的Matlab实现4.1 点估计法用少量确定性潮流逼近统计量点估计法的基本思想很取巧输入随机变量的分布不参与显式采样而是用输入变量的前几阶矩均值、方差、偏度等构造出若干个确定性估计点和对应权重对每个估计点做确定性潮流再对输出加权求和得到统计量。最基础的2点估计法规则如下对每个输入随机变量 x_i取两个估计点x_{i,1} μ_i σ_i x_{i,2} μ_i - σ_i权重各取 1/2。然后把第 i 个输入变量固定在这两个点上其他输入变量固定在均值处分别做两次确定性潮流。如果系统里有 n 个随机输入变量总共需要 2n 次潮流计算。相比蒙特卡洛动辄几千次计算量是天壤之别。输出变量的期望和方差按下面公式聚合E[Y] ≈ Σ_i Σ_k w_{i,k} * Y(x_{i,k}) E[Y^2] ≈ Σ_i Σ_k w_{i,k} * Y^2(x_{i,k}) Var[Y] E[Y^2] - (E[Y])^22点估计只用到均值和方差对线性系统是精确的对非线性系统会有截断误差。想要更高精度可以用3点估计额外引入偏度信息ξ_{i,1} λ3/2 sqrt(λ4 - 3λ3^2/4) ξ_{i,2} λ3/2 - sqrt(λ4 - 3λ3^2/4) ξ_{i,3} 0对应的估计点为 x_{i,k} μ_i ξ_{i,k} * σ_i。3点估计的权重计算链条稍长我在Matlab里建议直接用相关工具箱或者核对Hong在1999年原始文献的公式避免抄错。点估计法最大的短板是它只能给出输出的均值、方差等低阶矩不能直接恢复完整的概率密度分布。如果评估报告里必须画电压概率密度曲线点估计法就帮不上忙了。4.2 半不变量法 Gram-Charlier级数半不变量法走的是“矩-半不变量-级数展开”的分析路线。先利用输入随机变量的概率分布求出各阶半不变量然后在期望运行点做一次确定性潮流得到灵敏度矩阵再把输入半不变量线性映射到输出最后用Gram-Charlier或Edgeworth级数拟合输出分布。Matlab实现的核心步骤是这样第一步由输入随机变量的各阶矩计算半不变量。前四阶半不变量与矩的关系为κ1 μ1 κ2 μ2 - μ1^2 κ3 μ3 - 3μ1μ2 2μ1^3 κ4 μ4 - 4μ1μ3 6μ1^2μ2 - 3μ1^4第二步在系统期望运行点做一次牛顿法潮流取得雅可比矩阵求逆得到灵敏度矩阵 S0 J^(-1)。第三步线性映射。对于第i个输出量Y_i输入随机变量第s阶半不变量的贡献为 (S0(i,j))^s 乘以输入的第s阶半不变量再对j求和。公式为 κ_Y,s Σ_j (S0(i,j))^s * κ_W,s。第四步用Gram-Charlier级数把输出半不变量转换成概率密度。令 z (Y - μ_Y) / σ_Y则f(z) φ(z) * [1 (κ3/6σ^3) * He3(z) (κ4/24σ^4) * He4(z) ...]其中 φ(z) 是标准正态密度函数He3(z)z^3-3zHe4(z)z^4-6z^23。这套方法的计算速度最快适合在线评估场景。代价是灵敏度矩阵来自潮流方程的线性化对非线性强、重尾分布明显的系统展开到四阶筋的精度改善有限偶尔会出现概率密度曲线局部负值的现象这是级数截断本身带来的问题。4.3 三种方法怎么选方法精度计算量实现难度能否重建分布适用场景蒙特卡洛模拟最高无偏极大数千次潮流低能完整直方图标准分析、验证其他方法点估计法中等低阶矩精度高小2~3n次潮流中不能只有矩快速评估均值/标准差半不变量法线性化精度极小1次潮流映射高能近似解析分布在线评估、海量场景遍历从我实际使用的感受来说学术论文里最稳妥的套路是“蒙特卡洛做基准、点估计或半不变量做改进方法”。先用蒙特卡洛给出精确结果再展示改进方法在误差和耗时上的对比。如果直接上来就做半不变量法而缺失基准验证审稿人大概率会追问一句“和蒙特卡洛对比过吗”。5. 案例实操IEEE 14节点系统接入风光电源5.1 系统改造与参数设定这次演示以Matpower自带的case14为基础。IEEE 14节点系统有14个节点、5台发电机系统基准容量100MVA。我在原始算例基础上做了三处改动风电场接入节点9额定功率30MW。节点9原本是纯负荷节点用它接入风电后不需要改变发电机配置直接把注入功率折算成负负荷就行。光伏电站接入节点13额定功率15MW。同样处理为负负荷。所有负荷施加5%标准差的正态扰动截断到非负。风速Weibull参数取 c8.5m/s、k2.2切入风速3m/s、额定风速12m/s、切出风速25m/s。光照Beta分布参数取 a2、b1.5。这些参数偏理想化但演示概率潮流的完整流程足够了。如果要在实际工程中使用参数务必换成现场实测数据。5.2 完整代码实现与运行说明完整代码在第3章的框架基础上增加收敛判断、结果统计和可视化三个环节。我直接贴出循环结束后的统计部分% 剔除不收敛样本假设存于Vr中不收敛行全为0 Vr_valid Vr(all(Vr 1e-8, 2), :); Pa_valid Pa(all(Vr 1e-8, 2), :); % 节点电压统计指标 V_mean mean(Vr_valid, 1); V_std std(Vr_valid, 1); V_p5 prctile(Vr_valid, 5, 1); V_p95 prctile(Vr_valid, 95, 1); % 示例节点4的电压越限概率 prob_low mean(Vr_valid(:, 4) 0.95); prob_high mean(Vr_valid(:, 4) 1.05); fprintf(节点4电压均值 %.4f p.u.标准差 %.4f p.u.\n, V_mean(4), V_std(4)); fprintf(电压低于0.95概率%.4f%%高于1.05概率%.4f%%\n, prob_low*100, prob_high*100); % 支路过载概率有功超过线路容量1.0p.u.基准100MVA overload_prob mean(max(Pa_valid, [], 1) 1.0); fprintf(支路过载概率%.4f%%\n, overload_prob*100); % 绘制节点4电压幅值分布 figure; histogram(Vr_valid(:, 4), 80, Normalization, pdf); xlabel(节点4电压幅值 (p.u.)); ylabel(概率密度); title(节点4电压幅值概率分布5000次蒙特卡洛);运行这段代码需要提前确认Matlab环境具备了统计工具箱wblrnd、betarnd、histogram这些函数都依赖它和Matpower工具箱。Matlab版本我试过R2021b和R2023a都能跑通新版本没有遇到兼容性问题。如果你不想装Matpower也可以自己写牛顿法潮流函数但需要注意几个细节PV节点无功越限时要转换成PQ节点重新迭代平衡节点的相角要固定雅可比矩阵稀疏化用sparse构造不要用满阵否则系统规模一大内存直接爆掉。5.3 结果怎么看分布形态、越限概率与确定性解的差异我这次演示跑出来的典型结果大致是这样参数不同结果会有波动重点看分布形态节点4是系统中比较靠近负荷中心的节点电压均值大约在1.01p.u.标准差在0.012p.u.量级。这看起来波动幅度不大但分布尾部确实会越出 [0.95, 1.05] 的常规运行区间。支路过载概率非常低在千分位以下这符合case14网架结构相对坚强的特点。一个值得注意的现象是蒙特卡洛采样得到的电压均值往往不等于把所有随机变量固定在期望值时做确定性潮流得到的电压值。原因是潮流方程关于注入功率是高度非线性的电压幅值对注入的响应带有凸性期望值变换到了非线性函数内部就不再等价。这也是概率潮流区别于“把期望值代入确定性潮流”的根本原因。如果你在报告里写“风光出力取期望潮流算一遍结果即为系统平均运行状态”这在数学上是站不住脚的。审稿时这个问题是高频质疑点。6. 常见问题与排查技巧实录6.1 潮流不收敛怎么办蒙特卡洛循环里最烦人的就是跑着跑着某一次潮流不收敛。先用if res.success 0 continue把不收敛样本剔掉保证主程序不中断然后回过头排查不收敛的原因。我从实际调试经验看排在前面的原因有三个一是风光注入功率太大。当节点净负荷为负且数值很大时相当于一个功率倒送的发电机节点潮流方程可能走上一条不收敛的迭代路径。解决办法是检查注入功率是否超过系统承受能力适当降低风电场额定容量或者给该节点增加无功补偿设备。二是有功注入过大导致电压偏高触发发电机无功越限PV-PQ转换反复震荡。这种情况可以在潮流计算中打开无功越限处理选项或者调整该节点的无功补偿容量。三是采样到了极端恶化的负荷组合。当多个节点负荷同时处于波动上界时系统运行点可能逼近电压稳定边界。这时需要回溯样本参数看看是不是概率分布参数定得太激进。6.2 计算太慢怎么优化蒙特卡洛的耗时大头在重复潮流计算。同样的网络导纳矩阵结构每次都一样但Matpower每轮都会重新生成和分解。优化手段按收益排序优先把mpoption(out.all, 0)设上关闭MATPOWER的屏幕输出5000次仿真能省掉大约20%的IO时间。用parfor替代for这是最直接的提速手段。需要注意parfor里所有变量都必须符合切片规则我习惯把每次循环需要的数据预先构造成矩阵循环内只做索引切片。如果自己写牛顿法可以把雅可比矩阵中与网络拓扑相关的常数部分离线算好每次迭代只更新与节点注入相关的局部元素。这个方法能压掉不少时间但对代码能力有一定要求前期不建议折腾。对问题规模大、采样次数要求高的场景考虑用点估计法替代蒙特卡洛做快速预筛再用蒙特卡洛对高风险场景重点验算。6.3 结果异常的排查方向碰到概率分布形状诡异、均值偏移明显这类问题我通常按下面这个速查表逐项排查现象常见原因排查与解决电压均值明显偏低或偏高输入随机变量均值参数与基准工况不一致先将所有随机变量固定为期望值跑确定性潮流与Matpower基准结果比对风电功率样本出现负值Weibull函数参数顺序填反检查wblrnd调用正确形式为wblrnd(尺度c, 形状k)分布直方图出现双峰负荷截断过狠导致样本集中在零附近或风光参数组合形成多模态检查输入样本直方图单独绘制风速、光照分布形态概率密度曲线局部负值半不变量法级数截断造成的振荡增加展开阶数或改用Edgeworth级数或直接换蒙特卡洛复核潮流反复不收敛且集中在特定样本段该区间对应高渗透率极端场景检查该样本的风速、光照组合值评估是否超出系统静态稳定约束蒙特卡洛与点估计法的方差结果差异大系统非线性强低阶矩方法截断误差放大以蒙特卡洛为准增加点估计法的估计点数核验还有一个容易忽略的细节如果风光接入节点原本带负荷用负负荷等效后负荷波动生成器会对净负荷做扰动此时同一节点的注入波动和负荷波动被混在一起。严格来说这部分相关性在概率模型中并未分离。想处理干净就把“基础负荷”和“新能源注入”作为两个独立的随机源分开采样后在节点注入方程中相加代码里要预留对应的接口。7. 我踩过的坑和一点个人体会第一次跑通蒙特卡洛概率潮流的时候我用的还是纯手写的牛顿法潮流5000次仿真跑了将近十分钟Matlab风扇嗡嗡转结果电压均值比Matpower基准低了将近2%查了半天才发现是Weibull尺度参数c和形状参数k填反了。从那以后我养成了一个习惯任何随机分布参数进循环之前先单独生成一组样本画直方图目测形态是否合理。分布参数错了后面一切结果都是空中楼阁这一步省不得。另一个体会是方法论选型的顺序。我建议初学者不要一上来就钻研半不变量法和Gram-Charlier级数先用蒙特卡洛把“输入随机到输出随机”的直觉建立起来看懂电压分布是怎么来的再去研究怎么用更少的计算量逼近它。顺序反了的话公式推了一堆结果出了偏差你都不知道该怀疑是哪一步。如果后续想把这套东西扩展到工程应用两个大方向可以考虑一是把风光出力之间的空间相关性特别是同一气候区内多个风电场之间的出力相关性建进去否则风险评估会偏乐观二是结合时序运行模拟把风光出力的时间相关性考虑进来这样得到的电压越限概率才真正对应实际运行中持续时间的累积风险。概率潮流本身解决的是“截面不确定性”问题要和时序信息结合才完整覆盖新能源并网评估的整个拼图。
风电、光伏大规模接入之后电网分析的默认假设就变了。以前做潮流计算给一组确定的发电机出力和负荷牛拉法迭代一遍得到一组节点电压和支路功率这套流程在传统火电主导的时代够用但当出力看天吃饭的新能源占比上来输入侧的确定性假设本身就站不住脚。概率潮流计算Probabilistic Load FlowPLF正是用来应对这种不确定性的工具——它把风速、光照、负荷都建模成随机变量通过大规模采样或解析近似得到节点电压和支路功率的概率分布而不是一个孤零零的点值。Matlab凭借矩阵计算优势和自带的统计工具箱一直是做概率潮流最顺手的平台。这篇文章从工程应用的角度把含风光发电的概率潮流计算的数学模型、三种主流算法、Matlab代码框架和调试经验完整过一遍刚接触这个方向的研究生和做新能源接入评估的工程师都能直接拿去参考。1. 概率潮流到底在算什么从单点结果到概率分布的思维转换1.1 确定性潮流的“天花板”在哪里传统的确定性潮流求解的是这样一组方程节点注入功率等于电压与导纳矩阵的乘积方程组给定以后迭代求解出唯一一组节点电压和相角。问题的关键在于“给定”这两个字。系统里有风电和光伏之后注入功率不再是一个确定数值。同一座风电场年平均风速差一点全年发电量可能差出好几个百分点同一天内云层飘过光伏出力能从额定功率跌到零。此时如果你还用单点出力去做潮流分析得到的结果只对应某一种天气场景对系统规划人员来说参考价值非常有限。更麻烦的是电压越限这类风险恰恰容易出现在极端场景里大风天气下风电场满发而本地负荷处于低谷多余功率外送导致局部电压偏高傍晚光伏快速退出负荷却还在高位电压又会往下掉。确定性潮流算出来的“正常工况”电压往往看不出这些边界风险概率潮流正是把这类风险显式地量化出来。1.2 概率潮流的输入输出形态概率潮流做的事情本质上是输入的随机性到输出的随机性的传递。输入侧主要包括三类随机变量风速一般用两参数Weibull分布描述概率密度函数为 f(v) (k/c)(v/c)^(k-1) exp(-(v/c)^k)其中 c 是尺度参数k 是形状参数。光照强度通常用Beta分布描述因为光照强度在0到额定值之间连续变化Beta分布定义在有界区间上形态灵活。负荷功率一般用正态分布或对数正态分布波动范围取均值的3%~8%比较常见。输出侧就是电网运行人员真正关心的东西节点电压幅值的均值、标准差、概率密度函数和累计分布函数支路有功和无功的分布电压越上/下限概率支路过载概率系统网损的期望值。这些指标可以支撑三个层面的决策——规划阶段评估新能源接入容量是否过度激进运行阶段判断当前方式下电压风险水平调度阶段为备用容量和AGC调节留出合理区间。1.3 概率潮流与常规随机分析的差别有读者可能会问这和蒙特卡洛模拟直接撒点有什么区别蒙特卡洛只是概率潮流的实现手段之一。概率潮流本身是一套完整的方法论框架包含输入随机建模、相关性处理、不确定传递、输出统计四个环节。蒙特卡洛是最直观的传递手段点估计法和半不变量法则走解析路线用更少的计算量逼近同样的统计信息。后面会把这三种方法的Matlab实现逐一展开。2. 先把数学模型搭起来风光出力的随机特性描述2.1 风速与风电出力模型风速模型是整个概率潮流里最容易出错的环节因为风速分布和风机出力之间还隔着一道非线性分段函数。Matlab里生成Weibull随机数直接调用 wblrnd(scale, shape) 即可但这里有个经典坑位Matlab的wblrnd第一个参数是尺度参数c单位m/s第二个是形状参数k无量纲我见过不止一次有人把参数填反结果生成的风速样本整体偏小或偏大风电出力分布完全失真。风速到电功率转换用最通用的分段线性模型P_w(v) 0, v v_in 或 v v_out P_w(v) P_rated * (v - v_in)/(v_rated - v_in), v_in v v_rated P_w(v) P_rated, v_rated v v_outv_in 一般取3m/sv_rated 取11~13m/sv_out 取25m/s。这个模型虽然简单但在工程分析里足够用。更精确的风功率曲线可以用厂商实测数据插值不过概率潮流关注的是长期统计分布分段线性模型引入的误差通常可以接受。Weibull参数的获取有两种途径一是用当地测风塔一整年的小时级平均风速数据通过极大似然估计反推二是用平均风速和标准差近似估算。后者在参考资料匮乏的时候很实用c 约等于平均风速的1.12倍k 则通过变异系数查表或数值求解。估算出来的参数用于方案预研足够正式工程评估还是建议用实测数据。2.2 光照强度与光伏出力模型光伏出力的建模同样分两步。第一步用Beta分布描述光照强度的随机性Beta分布的概率密度函数为 f(s) (Γ(ab))/(Γ(a)Γ(b)) * s^(a-1) * (1-s)^(b-1)Matlab里用 betarnd(a, b) 直接生成0到1之间的标幺值。a和b的取值决定了分布的偏斜程度夏季晴天多曲线右偏a取2~3、b取1~1.5比较贴合阴雨多的地方分布相对左偏参数相应调整。第二步做光电转换。理想情况下光伏出力与光照强度近似线性P_pv η * S * A其中η是光电转换效率S是实际光照强度A是光伏阵列面积。实际工程里更常用的是容量标幺法P_pv P_rated * (S / S_ref)S_ref通常取1000W/m²。温度对光伏出力的影响也不能完全忽略组件温度升高导致输出电压下降更精细的模型会在上述公式基础上乘一个温度修正系数(1 - β(T - 25))β一般取0.003~0.005 /°C。对于概率潮流温度修正要不要做得这么细取决于计算目的——如果只是评估年度电压分布简化模型足矣如果要做夏季高温极端场景分析建议把温度项加进来。2.3 负荷随机性与相关性处理负荷波动用正态分布描述最普遍但需要加截断处理。我在Matlab里习惯写成PD_load PD_mean .* (1 0.05 * randn(n, 1)); PD_load(PD_load 0) 0;不做截断的话理论上会生成负负荷样本虽然概率极低但一旦出现就会导致潮流方程出现负注入结果根本没法解释。0.05这个系数对应5%的标准差如果节点负荷本身波动大可以放宽到8%~10%。另一个容易被忽视的问题是相关性。风电场和光伏电站如果处在同一个区域电网风速和光照之间可能存在相关性不同节点的负荷也不会完全独立。忽略相关性的后果是低估系统电压波动的风险区间。Matlab里处理相关高斯随机变量最直接的方式是用Cholesky分解给定相关系数矩阵R计算 L chol(R)然后把独立标准正态样本矩阵 Z 乘以 L 的转置得到带相关性的样本。对于Weibull和Beta这类非高斯变量更严格的做法是引入Copula但工程上如果只是粗略评估Cholesky分解后做概率积分变换也能用。3. 核心方法一蒙特卡洛模拟最稳也最烧钱3.1 三步走采样、潮流计算、统计蒙特卡洛模拟的思路直白到不需要过多解释既然输入是随机变量那就生成大量输入样本逐个做确定性潮流计算最后把所有输出结果汇集起来做统计分析。三个步骤对应三段Matlab代码每一步都可以调优。第一步是采样。样本量太小分布形状出不来样本量太大计算耗时线性增长。第二步是潮流计算每次都调用一次确定性潮流求解器。第三步是统计mean、std、histogram、prctile几个函数就够用。蒙特卡洛最大的优点是“无偏”——只要样本量足够大输出分布一定收敛到真实分布不依赖任何线性化假设。这点在后两种解析方法里是做不到的。代价就是计算量大。如果单次牛顿法潮流耗时0.01秒5000次就是50秒在多节点系统里单次潮流可能到0.1秒5000次就是8分钟级别。3.2 Matlab主程序框架建议用Matpower做潮流计算内核自己手写牛顿法在原理上没问题但处理PV-PQ节点转换、无功越限这些细节时容易踩坑。Matpower自带IEEE 14节点等标准算例loadcase一键加载runpf闭环求解省心很多。下面这套代码是我的常用框架对应风电和光伏分别接入两个不同节点的情况% 蒙特卡洛概率潮流主框架Matpower 统计工具箱 mpc loadcase(case14); base_PD mpc.bus(:, 3); % 保存原始有功负荷 base_QD mpc.bus(:, 4); % 保存原始无功负荷 N 5000; Vr zeros(N, 14); % 电压幅值记录矩阵 Pa zeros(N, 41); % 支路有功记录矩阵case14共41条支路 Sr zeros(N, 14); % 节点注入视在功率记录 % 风速Weibull参数wblrnd(尺度c, 形状k)顺序不要写反 c_w 8.5; k_w 2.2; v_in 3; v_rated 12; v_out 25; Pw_rated 30; % 风电场额定功率 MW % 光照Beta分布参数 a_s 2; b_s 1.5; Pp_rated 15; % 光伏电站额定功率 MW for i 1:N % 1. 采样风速 - 风电出力 v wblrnd(c_w, k_w); if v v_in || v v_out Pw 0; elseif v v_rated Pw Pw_rated * (v - v_in) / (v_rated - v_in); else Pw Pw_rated; end % 2. 采样光照强度 - 光伏出力 s_pu betarnd(a_s, b_s); Pp Pp_rated * s_pu; % 3. 把风光出力等效为相应节点的负负荷注入 PD_node base_PD; PD_node(9) base_PD(9) - Pw / 100; % 风电接入bus9 PD_node(13) base_PD(13) - Pp / 100; % 光伏接入bus13 % 4. 负荷波动对包含负注入的净负荷施加正态扰动 PD_load PD_node .* (1 0.05 * randn(14, 1)); QD_load base_QD .* (1 0.05 * randn(14, 1)); PD_load(PD_load 0) 0; QD_load(QD_load 0) 0; mpc.bus(:, 3) PD_load; mpc.bus(:, 4) QD_load; % 5. 确定性潮流计算 res runpf(mpc, mpoption(out.all, 0)); if res.success 0 warning(第 %d 次潮流不收敛, i); continue; end % 6. 记录输出量 Vr(i, :) res.bus(:, 8); % 电压幅值 Pa(i, :) res.branch(:, 14); % 支路有功 end注意代码第3步和第4步的顺序。先把风电和光伏出力折算成负的净负荷再对这种净负荷施加正态扰动等价于“风光出力随机 负荷随机”的叠加方式。如果你先把原负荷扰动完再把风光注入单独减掉那么在注入很大的节点上净负荷可能出现负值且分布形态被扭歪。另外代码里 Pw/100 是因为case14的基准容量是100MVA把所有功率统一折算到标幺值这一点新手特别容易漏。3.3 收敛性判断与采样规模选择到底采多少组样本才算够经验法则是先看输出量的均值或标准差随样本数的变化曲线当相对波动小于某个阈值时认为收敛。我常用的判据是电压均值的无穷范数误差连续两次采样之间的变化量小于1e-4就停止。更直接的做法是固定样本量5000次起步不好再翻倍到10000次。N_max 10000; V_mean_old zeros(1, 14); for i 1:N_max % …… 上述采样与潮流计算代码 …… V_mean_new mean(Vr(1:i, :), 1); delta max(abs(V_mean_new - V_mean_old)); if delta 1e-4 i 1000 fprintf(均值收敛于第 %d 次采样\n, i); break; end V_mean_old V_mean_new; end这个小循环跑起来有个好处你不用赌样本量机器自己告诉你够了。代价是循环内多了mean运算对整体耗时影响很小。蒙特卡洛还天然支持并行化把for改成parfor前提是循环体内不能有依赖全局变量的操作runpf的mpc结构每次都是基于本次采样数据构造的满足parfor要求。我实测过在8核机器上3000次仿真的耗时能压到原来的三分之一左右。注意parfor里mpc这个变量会被当作广播变量处理数据量不大影响有限。4. 省时方案点估计法与半不变量法的Matlab实现4.1 点估计法用少量确定性潮流逼近统计量点估计法的基本思想很取巧输入随机变量的分布不参与显式采样而是用输入变量的前几阶矩均值、方差、偏度等构造出若干个确定性估计点和对应权重对每个估计点做确定性潮流再对输出加权求和得到统计量。最基础的2点估计法规则如下对每个输入随机变量 x_i取两个估计点x_{i,1} μ_i σ_i x_{i,2} μ_i - σ_i权重各取 1/2。然后把第 i 个输入变量固定在这两个点上其他输入变量固定在均值处分别做两次确定性潮流。如果系统里有 n 个随机输入变量总共需要 2n 次潮流计算。相比蒙特卡洛动辄几千次计算量是天壤之别。输出变量的期望和方差按下面公式聚合E[Y] ≈ Σ_i Σ_k w_{i,k} * Y(x_{i,k}) E[Y^2] ≈ Σ_i Σ_k w_{i,k} * Y^2(x_{i,k}) Var[Y] E[Y^2] - (E[Y])^22点估计只用到均值和方差对线性系统是精确的对非线性系统会有截断误差。想要更高精度可以用3点估计额外引入偏度信息ξ_{i,1} λ3/2 sqrt(λ4 - 3λ3^2/4) ξ_{i,2} λ3/2 - sqrt(λ4 - 3λ3^2/4) ξ_{i,3} 0对应的估计点为 x_{i,k} μ_i ξ_{i,k} * σ_i。3点估计的权重计算链条稍长我在Matlab里建议直接用相关工具箱或者核对Hong在1999年原始文献的公式避免抄错。点估计法最大的短板是它只能给出输出的均值、方差等低阶矩不能直接恢复完整的概率密度分布。如果评估报告里必须画电压概率密度曲线点估计法就帮不上忙了。4.2 半不变量法 Gram-Charlier级数半不变量法走的是“矩-半不变量-级数展开”的分析路线。先利用输入随机变量的概率分布求出各阶半不变量然后在期望运行点做一次确定性潮流得到灵敏度矩阵再把输入半不变量线性映射到输出最后用Gram-Charlier或Edgeworth级数拟合输出分布。Matlab实现的核心步骤是这样第一步由输入随机变量的各阶矩计算半不变量。前四阶半不变量与矩的关系为κ1 μ1 κ2 μ2 - μ1^2 κ3 μ3 - 3μ1μ2 2μ1^3 κ4 μ4 - 4μ1μ3 6μ1^2μ2 - 3μ1^4第二步在系统期望运行点做一次牛顿法潮流取得雅可比矩阵求逆得到灵敏度矩阵 S0 J^(-1)。第三步线性映射。对于第i个输出量Y_i输入随机变量第s阶半不变量的贡献为 (S0(i,j))^s 乘以输入的第s阶半不变量再对j求和。公式为 κ_Y,s Σ_j (S0(i,j))^s * κ_W,s。第四步用Gram-Charlier级数把输出半不变量转换成概率密度。令 z (Y - μ_Y) / σ_Y则f(z) φ(z) * [1 (κ3/6σ^3) * He3(z) (κ4/24σ^4) * He4(z) ...]其中 φ(z) 是标准正态密度函数He3(z)z^3-3zHe4(z)z^4-6z^23。这套方法的计算速度最快适合在线评估场景。代价是灵敏度矩阵来自潮流方程的线性化对非线性强、重尾分布明显的系统展开到四阶筋的精度改善有限偶尔会出现概率密度曲线局部负值的现象这是级数截断本身带来的问题。4.3 三种方法怎么选方法精度计算量实现难度能否重建分布适用场景蒙特卡洛模拟最高无偏极大数千次潮流低能完整直方图标准分析、验证其他方法点估计法中等低阶矩精度高小2~3n次潮流中不能只有矩快速评估均值/标准差半不变量法线性化精度极小1次潮流映射高能近似解析分布在线评估、海量场景遍历从我实际使用的感受来说学术论文里最稳妥的套路是“蒙特卡洛做基准、点估计或半不变量做改进方法”。先用蒙特卡洛给出精确结果再展示改进方法在误差和耗时上的对比。如果直接上来就做半不变量法而缺失基准验证审稿人大概率会追问一句“和蒙特卡洛对比过吗”。5. 案例实操IEEE 14节点系统接入风光电源5.1 系统改造与参数设定这次演示以Matpower自带的case14为基础。IEEE 14节点系统有14个节点、5台发电机系统基准容量100MVA。我在原始算例基础上做了三处改动风电场接入节点9额定功率30MW。节点9原本是纯负荷节点用它接入风电后不需要改变发电机配置直接把注入功率折算成负负荷就行。光伏电站接入节点13额定功率15MW。同样处理为负负荷。所有负荷施加5%标准差的正态扰动截断到非负。风速Weibull参数取 c8.5m/s、k2.2切入风速3m/s、额定风速12m/s、切出风速25m/s。光照Beta分布参数取 a2、b1.5。这些参数偏理想化但演示概率潮流的完整流程足够了。如果要在实际工程中使用参数务必换成现场实测数据。5.2 完整代码实现与运行说明完整代码在第3章的框架基础上增加收敛判断、结果统计和可视化三个环节。我直接贴出循环结束后的统计部分% 剔除不收敛样本假设存于Vr中不收敛行全为0 Vr_valid Vr(all(Vr 1e-8, 2), :); Pa_valid Pa(all(Vr 1e-8, 2), :); % 节点电压统计指标 V_mean mean(Vr_valid, 1); V_std std(Vr_valid, 1); V_p5 prctile(Vr_valid, 5, 1); V_p95 prctile(Vr_valid, 95, 1); % 示例节点4的电压越限概率 prob_low mean(Vr_valid(:, 4) 0.95); prob_high mean(Vr_valid(:, 4) 1.05); fprintf(节点4电压均值 %.4f p.u.标准差 %.4f p.u.\n, V_mean(4), V_std(4)); fprintf(电压低于0.95概率%.4f%%高于1.05概率%.4f%%\n, prob_low*100, prob_high*100); % 支路过载概率有功超过线路容量1.0p.u.基准100MVA overload_prob mean(max(Pa_valid, [], 1) 1.0); fprintf(支路过载概率%.4f%%\n, overload_prob*100); % 绘制节点4电压幅值分布 figure; histogram(Vr_valid(:, 4), 80, Normalization, pdf); xlabel(节点4电压幅值 (p.u.)); ylabel(概率密度); title(节点4电压幅值概率分布5000次蒙特卡洛);运行这段代码需要提前确认Matlab环境具备了统计工具箱wblrnd、betarnd、histogram这些函数都依赖它和Matpower工具箱。Matlab版本我试过R2021b和R2023a都能跑通新版本没有遇到兼容性问题。如果你不想装Matpower也可以自己写牛顿法潮流函数但需要注意几个细节PV节点无功越限时要转换成PQ节点重新迭代平衡节点的相角要固定雅可比矩阵稀疏化用sparse构造不要用满阵否则系统规模一大内存直接爆掉。5.3 结果怎么看分布形态、越限概率与确定性解的差异我这次演示跑出来的典型结果大致是这样参数不同结果会有波动重点看分布形态节点4是系统中比较靠近负荷中心的节点电压均值大约在1.01p.u.标准差在0.012p.u.量级。这看起来波动幅度不大但分布尾部确实会越出 [0.95, 1.05] 的常规运行区间。支路过载概率非常低在千分位以下这符合case14网架结构相对坚强的特点。一个值得注意的现象是蒙特卡洛采样得到的电压均值往往不等于把所有随机变量固定在期望值时做确定性潮流得到的电压值。原因是潮流方程关于注入功率是高度非线性的电压幅值对注入的响应带有凸性期望值变换到了非线性函数内部就不再等价。这也是概率潮流区别于“把期望值代入确定性潮流”的根本原因。如果你在报告里写“风光出力取期望潮流算一遍结果即为系统平均运行状态”这在数学上是站不住脚的。审稿时这个问题是高频质疑点。6. 常见问题与排查技巧实录6.1 潮流不收敛怎么办蒙特卡洛循环里最烦人的就是跑着跑着某一次潮流不收敛。先用if res.success 0 continue把不收敛样本剔掉保证主程序不中断然后回过头排查不收敛的原因。我从实际调试经验看排在前面的原因有三个一是风光注入功率太大。当节点净负荷为负且数值很大时相当于一个功率倒送的发电机节点潮流方程可能走上一条不收敛的迭代路径。解决办法是检查注入功率是否超过系统承受能力适当降低风电场额定容量或者给该节点增加无功补偿设备。二是有功注入过大导致电压偏高触发发电机无功越限PV-PQ转换反复震荡。这种情况可以在潮流计算中打开无功越限处理选项或者调整该节点的无功补偿容量。三是采样到了极端恶化的负荷组合。当多个节点负荷同时处于波动上界时系统运行点可能逼近电压稳定边界。这时需要回溯样本参数看看是不是概率分布参数定得太激进。6.2 计算太慢怎么优化蒙特卡洛的耗时大头在重复潮流计算。同样的网络导纳矩阵结构每次都一样但Matpower每轮都会重新生成和分解。优化手段按收益排序优先把mpoption(out.all, 0)设上关闭MATPOWER的屏幕输出5000次仿真能省掉大约20%的IO时间。用parfor替代for这是最直接的提速手段。需要注意parfor里所有变量都必须符合切片规则我习惯把每次循环需要的数据预先构造成矩阵循环内只做索引切片。如果自己写牛顿法可以把雅可比矩阵中与网络拓扑相关的常数部分离线算好每次迭代只更新与节点注入相关的局部元素。这个方法能压掉不少时间但对代码能力有一定要求前期不建议折腾。对问题规模大、采样次数要求高的场景考虑用点估计法替代蒙特卡洛做快速预筛再用蒙特卡洛对高风险场景重点验算。6.3 结果异常的排查方向碰到概率分布形状诡异、均值偏移明显这类问题我通常按下面这个速查表逐项排查现象常见原因排查与解决电压均值明显偏低或偏高输入随机变量均值参数与基准工况不一致先将所有随机变量固定为期望值跑确定性潮流与Matpower基准结果比对风电功率样本出现负值Weibull函数参数顺序填反检查wblrnd调用正确形式为wblrnd(尺度c, 形状k)分布直方图出现双峰负荷截断过狠导致样本集中在零附近或风光参数组合形成多模态检查输入样本直方图单独绘制风速、光照分布形态概率密度曲线局部负值半不变量法级数截断造成的振荡增加展开阶数或改用Edgeworth级数或直接换蒙特卡洛复核潮流反复不收敛且集中在特定样本段该区间对应高渗透率极端场景检查该样本的风速、光照组合值评估是否超出系统静态稳定约束蒙特卡洛与点估计法的方差结果差异大系统非线性强低阶矩方法截断误差放大以蒙特卡洛为准增加点估计法的估计点数核验还有一个容易忽略的细节如果风光接入节点原本带负荷用负负荷等效后负荷波动生成器会对净负荷做扰动此时同一节点的注入波动和负荷波动被混在一起。严格来说这部分相关性在概率模型中并未分离。想处理干净就把“基础负荷”和“新能源注入”作为两个独立的随机源分开采样后在节点注入方程中相加代码里要预留对应的接口。7. 我踩过的坑和一点个人体会第一次跑通蒙特卡洛概率潮流的时候我用的还是纯手写的牛顿法潮流5000次仿真跑了将近十分钟Matlab风扇嗡嗡转结果电压均值比Matpower基准低了将近2%查了半天才发现是Weibull尺度参数c和形状参数k填反了。从那以后我养成了一个习惯任何随机分布参数进循环之前先单独生成一组样本画直方图目测形态是否合理。分布参数错了后面一切结果都是空中楼阁这一步省不得。另一个体会是方法论选型的顺序。我建议初学者不要一上来就钻研半不变量法和Gram-Charlier级数先用蒙特卡洛把“输入随机到输出随机”的直觉建立起来看懂电压分布是怎么来的再去研究怎么用更少的计算量逼近它。顺序反了的话公式推了一堆结果出了偏差你都不知道该怀疑是哪一步。如果后续想把这套东西扩展到工程应用两个大方向可以考虑一是把风光出力之间的空间相关性特别是同一气候区内多个风电场之间的出力相关性建进去否则风险评估会偏乐观二是结合时序运行模拟把风光出力的时间相关性考虑进来这样得到的电压越限概率才真正对应实际运行中持续时间的累积风险。概率潮流本身解决的是“截面不确定性”问题要和时序信息结合才完整覆盖新能源并网评估的整个拼图。
