搞电力系统不确定性的朋友应该都有个共同体会风光一多传统确定性潮流算完就得出一个“标准答案”但现场负荷和出力天天在变那个单点结果基本只能当安慰剂用。我最近拿IEEE33节点系统当试验对象把蒙特卡洛法跟风光出力模型绑在一起做概率潮流计算跑完一堆场景之后发现这活儿看起来简单实际坑不少但用好了信息量非常可观。这篇文章就把我整个折腾过程、代码逻辑、参数选择和踩过的坑完整记录下来希望给正在搞不确定性分析的同学一点可参考的操作路径。蒙特卡洛法在概率潮流里的定位其实很直接对风光出力、负荷波动这些随机变量进行大量采样然后反复求解确定性潮流最后对节点电压、支路功率做统计分析得到概率分布和越限风险。IEEE33节点系统结构经典、参数公开用来做算法验证再合适不过。搭配风光出力模型后不仅能模拟源端随机性还能把光伏和风电的相关性、时序特性一并考虑进去算是入门概率潮流最不折腾的组合。1. 为什么拿IEEE33节点当试验田1.1 不确定性分析到底在解决什么问题传统潮流计算的前提是运行工况固定负荷和出力都取确定值。但实际电网里分布式光伏随风速和光照波动、负荷随着季节和时段变化任何一个输入量的偏离都会导致节点电压或支路潮流偏离计算值。一旦偏离超出设备限值就可能出现过电压、线路过载甚至保护误动。所以我们需要回答的不是“某个时刻系统安不安全”而是“在一年或者一天内系统有多大概率越限”。概率潮流的本质就是把输入随机变量的分布映射到输出量分布量化风险。这里蒙特卡洛法是数理统计里最朴素也最鲁棒的思路只要采样次数足够多计算结果就会依概率收敛到真实分布。它不像解析法那样要做大量线性化假设也不像点估计法那样丢失高阶信息很适合前期探索和精度验证。1.2 IEEE33节点的优势与仿真意义IEEE33节点系统是配电领域特别常用的算例拓扑是单电源辐射网包含33个节点、32条支路首端是平衡节点负荷有有功和无功还有联络开关可以处理重构场景。这个规模很讨巧潮流收敛快、计算量低但包含了支路末端电压偏低、轻载重载切换这类特性足够展示概率潮流结果里的细节差异。拿它做蒙特卡洛试验还有一个好处公开参数里负荷基准值是3715kW和2300kVar恰好可以在此基础上叠加风光模型。我们在节点18接入风电、节点22接入光伏容量按渗透率设定就能清楚看到随机源接入后末端电压和馈线潮流的分布变化。这种算例复现门槛低非常适合先跑通流程再迁移到更大系统。2. 风光出力模型该怎么搭2.1 风速与光照的概率特性风光出力模型的核心是描述随机出力背后的概率分布。风速常用两参数威布尔分布概率密度函数记作[ f(v)\frac{k}{\lambda}\left(\frac{v}{\lambda}\right)^{k-1}\exp\left(-\left(\frac{v}{\lambda}\right)^k\right) ]其中形状参数k和尺度参数λ决定风速的峰态和平均量级。实际应用中如果这里采集的是某地区风电场风速数据可以用最大似然估计或者经验公式拟合出参数。风电出力与风速之间通常用分段线性曲线表示PWind在切入风速以下或切出风速以上为0额定风速以上封顶中间段按线性或二次插值。光伏出力的随机性来自光照辐照度一般用Beta分布描述概率密度为[ f(r)\frac{\Gamma(\alpha\beta)}{\Gamma(\alpha)\Gamma(\beta)}\left(\frac{r}{r_{max}}\right)^{\alpha-1}\left(1-\frac{r}{r_{max}}\right)^{\beta-1} ]Alpha和Beta可以根据历史辐照数据的均值和方差反推。在温度波动不大时光伏有功出力可以简化为标准光照条件下的线性折算。这样每个采样时刻就能生成一组与风速和光照对应的风光出力值。2.2 相关性怎么处理单独采样风速和光照会忽略它们之间的天然联系。典型城区中风速白天偏大、夜晚偏小而光照中午最大、傍晚为零二者存在时段上的负相关。如果采样时完全独立得到的场景组合可能过于极端或脱离实际天气过程。处理相关性比较主流的方法是使用Copula函数把边缘分布和相关性结构分开建模。常用有Gaussian Copula和Frank Copula前者参数容易根据秩相关矩阵估算后者对尾部依赖描述更精细。实操中我先对历史风速、光照数据计算秩相关系数然后用该系数构造相关矩阵再生成服从该相关结构的标准正态随机数最后做逆变换得到对应的风速和光照样本。这样生成的场景序列和原始数据统计特性高度一致。2.3 参数选择与场景生成风光出力模型里最容易翻车的是参数选择不靠数据而靠拍脑袋。我建议第一步先获取本地风电场和光伏电站的小时级历史数据至少覆盖一年8760个小时然后分场景如按月份或按天气类型拟合参数。如果没有真实数据也可以参照公开研究中的典型参数风速Weibull分布的k取2.0到2.5λ取6到8光伏Beta分布的α、β取4到6。关键是要在后续蒙特卡洛采样前先对模型做统计检验比如把采样生成的出力分布与原始数据直方图对比偏差大就要重新拟合。场景生成数量取决于蒙特卡洛规模。我一般先抽500次画直方图看分布形态再逐步增大到5000或10000次观察均值和方差的收敛情况。代码里风光场景采样通常放在主循环外层一次生成全部样本再逐个代入潮流比每次循环里重复调用分布函数更快。3. 蒙特卡洛概率潮流计算实操3.1 整体流程拆解概率潮流分析的完整流程分四步定义输入随机变量、生成样本场景、求解确定性潮流、统计输出结果。第一步要明确哪些节点是风电、哪些是光伏负荷在哪个节点什么范围波动。第二步上一个章节已经说了用Copula采样获得一组组风速光照样本再换算成风机和光伏的有功功率输出。第三步是核心计算环节每组样本中把风光出力替代对应节点的注入功率同时给负荷节点叠加正态扰动然后调用牛顿拉夫逊法或前推回代法求解潮流。辐射配电网用前推回代法宰得更稳但为保持通用性我用Matpower里的runpf函数它可以指定每种潮流方法。第四步统计结果对所有样本求得节点电压幅值、支路有功/无功潮流后计算均值、标准差、最大值、最小值以及节点电压越上限/下限概率、线路过载概率。输出如果做成概率密度曲线或者累积分布曲线会比单纯列数字更直观。3.2 核心代码逻辑逐段讲解实际工程中我用Matlab做全流程代码规模不大。下面这段是主循环的核心骨架我加了注释说明每一步作用。先定义随机变量维数和样本矩阵% 初始化 N 10000; % 蒙特卡洛样本数 wind_nodes [18]; % 节点18接风电 pv_nodes [22]; % 节点22接光伏 Pg zeros(33, N); % 存储每个节点的注入有功样本 for i 1:N % 1. 抽取Copula相关风速和光照 U mvnrnd([0,0], [1, c; c, 1], 1); U normcdf(U); % 转换到均匀分布 v_sample wblinv(U(1), lambda_w, k_w); % 风速 r_sample betainv(U(2), alpha_pv, beta_pv); % 光照 % 2. 根据风速-出力转化曲线计算风电功率 P_wind wind_curve(v_sample); % 自定义函数 P_pv pv_curve(r_sample); % 自定义函数 % 3. 叠加负荷波动正态分布 load_noise normrnd(1, 0.05, 33, 1); P_load base_load .* load_noise; % 4. 构造各节点注入功率 P_net -P_load; % 负荷为负注入 P_net(wind_nodes) P_net(wind_nodes) P_wind; P_net(pv_nodes) P_net(pv_nodes) P_pv; % 5. 保存该样本的节点注入 Pg(:, i) P_net; end这里的核心细节是把“随机数生成”和“潮流计算”解耦。先生成全部输入样本再用矩阵运算或并行循环依次求解潮流这样可以利用Matlab的parfor并行加速。如果按照边采样边潮流的老写法每次循环里都要初始化潮流参数实际耗时能差一个数量级。风速转功率的wind_curve函数在我这里用的是典型三段式曲线function Pw wind_curve(v) v_in 3; v_out 25; v_rated 12; P_rated 500; % 单位kW if v v_in || v v_out Pw 0; elseif v v_rated Pw P_rated * (v - v_in) / (v_rated - v_in); else Pw P_rated; end end光伏输出采用线性折算光照超过额定辐照度就限幅。这些函数都要单独测试一下边界情况尤其是风速在切入风速附近以及光照接近零的情况否则会产出奇怪的负功率或超大功率。潮流计算部分我对每列样本调用一次runpf并记录电压和潮流结果。值得注意Matpower默认情况下如果潮流不收敛会直接报错中断。我采用try-catch结构记录不收敛场景数量并把不收敛场景置为上一次成功结果或丢弃。如果采样场景导致潮流不收敛次数太多大概率是生成的出力或负荷波动范围设置过宽需要回头检查分布参数。V_amp zeros(33, N); S_branch zeros(32, N); conv_flags zeros(1, N); parfor i 1:N mpc loadcase(case33bus); % IEEE33节点的标准潮流数据 mpc.bus(:, 2) mpc.bus(:, 2) - Pg(:, i) / 1000; % 注意单位转换Matpower基准功率单位是MW mpc.bus(:, 3) mpc.bus(:, 3) - Qg(:, i) / 1000; % 无功类似处理 try opt mpoption(pf.alg, NR, out.all, 0); result runpf(mpc, opt); V_amp(:, i) result.bus(:, 8); S_branch(:, i) abs(result.branch(:, 14) 1i * result.branch(:, 15)); conv_flags(i) 1; catch conv_flags(i) 0; end end3.3 关键参数设置与收敛性判断蒙特卡洛法的精度只受到两个因素影响一个是样本量另一个是输入分布和实际物理过程的匹配程度。样本量不是越大越好因为当前推回代法在这种配电网规模下10000次计算可能在几分钟内完成但很多情况下5000次已经足够稳定。判断样本量是否足够的常用指标是输出均值随样本数的收敛曲线。我一般会设定一个滑窗计算前200个样本的平均电压之后每增加100个样本更新一次看数值波动是否低于0.0001pu。如果低于这个阈值说明再多样本也只是把尾部分布描得更细不会改变整体结论。如果想精确预测极端场景比如概率5%以下的越限样本量需要明显增大至少到20000次以上采用重要抽样或拉丁超立方采样会更划算。参数设置里还有个容易忽略的点负荷波动标准差。很多人习惯给所有节点统一加5%波动但实际配电网里不同节点的负荷性质差异很大居民负荷波动大商业负荷波动小。如果能拿到节点类型信息最好分节点设置不同标准差。实在没有数据至少也要让负荷波动均值归一化后保持原有基准值平衡否则会造成系统总负荷偏差过大让潮流结果整体偏移。4. 结果分析与坑位记录4.1 典型输出与分析指标跑完10000个场景后我习惯先输出节点电压幅值的概率分布。可以画成盒须图横轴是33个节点纵轴是电压标幺值每个节点展示最小、四分之一分位、中位数、四分之三分位和最大值。通常在末端节点18和22附近能明显看到电压分布更宽这是接入风光后的典型特征末端电压不仅受负荷波动影响还受随机出力变化影响区间范围比首段节点大得多。另一个值得分析的是支路潮流越限概率。分支潮流越限概率定义为P(Si Simax)其中Simax是支路容量。当风电大发、负荷低谷时潮流反向或达到峰值更容易触发越限。我计算出的结果往往能直观显示是哪几条联络支路风险最大这对接下来的网络重构很有参考价值。分析层面还可以进一步做灵敏度分析比如计算电压越限概率对风电接入容量的灵敏度用来评估并网容量上限。4.2 蒙特卡洛法常见坑与替代思路第一个坑是伪随机数种子不一致。同一套代码如果不同人运行没有固定随机种子结果差异可能很大。所以在发布代码或做对比实验时务必用rng(固定数字)初始化保证可复现。第二个坑是忽略相关性。独立采样导致风电和光伏同时出高功率的场景发生频率偏高会让系统高渗透率场景被高估越限概率被放大。我试过Copula修正前后末端电压越上限概率能差出近3倍这个差异足以影响并网决策。第三个坑是样本量不足却强行看极值。想比较节点电压的最大值和最小值500次采样得到的极值根本不可信因为极值依赖稀有的极端场景。要稳定估计尾部真得使用分层采样技术。我后来把普通蒙特卡洛换成拉丁超立方采样在同一样本量下电压标准差和越限概率的方差明显减小代码改动不大就是在采样前对均匀随机数做了一次排序。第四个坑是潮流求解器设置不对。有些开源潮流工具对配电网多平衡节点的处理并不好尤其当风电场容量占比高时如果只把大电网平衡节点设为唯一平衡源系统很容易因为无功不足而迭代不收敛。实际处理方法要么换成分布式的松弛节点要么在风光场站配置无功补偿模型。4.3 实用建议如果只是想快速评估方案优劣没必要一上来就跑10000次。先跑200次看趋势再针对高风险场景加密采样这就是所谓序贯蒙特卡洛思路。具体做法是第一轮粗采样后用KDE估计越限区域的概率密度然后在越限区域增加权重重新采样得到的估计方差最小。工程实践中我还建议把概率潮流结果和时序生产模拟结合。单纯随机采样得到的是“不考虑时间顺序”的概率结果但储能调度、爬坡需求都要看时序性。我在项目里的做法是先按小时划分场景在每个小时段内做蒙特卡洛这样得到的逐时概率分布能直接给调度员提供风险预警窗口。结果展示时别只给均值。把“5%分位数”“95%分位数”和“越限概率”三件套做成表格或热力图运维人员看到才知道风险集中在哪条馈线、哪个时段。我这里整理了自己用的一组典型结果格式给大家参考节点电压均值(pu)电压5%分位电压95%分位越上限概率越下限概率180.9750.9631.0212.1%3.2%220.9820.9691.0181.8%2.5%330.9580.9440.9890.4%8.9%支路过载概率可以单独整理成另一张表标明支路编号、线路容量、最大负载率和过载概率。通过这些表格规划人员能直接找准薄弱环节而不是大海捞针。最后再分享一点个人经验做这类仿真最消耗精力的不是潮流算法而是输入数据的清洗和概率模型的校验。我前几版跑出来结果特别夸张后来一查是光照辐照度没有做昼夜归属把夜间零辐照和白天高辐照混在一起导致Beta分布拟合结果失真。花点时间把历史数据按时间切片、按天气聚类再拟合参数出来的概率潮流才真正有工程参考价值。概率潮流不是炫技舞台它解决的是“系统风险能不能被看到”的问题。把数据模型做扎实蒙特卡洛这把牛刀才不至于切碎了豆腐。