基于蒙特卡洛的风电光伏功率场景生成：Matlab实现与工程实践

搞新能源并网和电力系统调度的同行对“场景生成”这个词应该都不陌生。风电不像火电你说它明天出力多少它偏偏给你来一段低风速光伏同样一片云飘过来功率曲线当场掉头。传统做法是给一个预测值再给一个预测误差区间但是遇到日内滚动优化、储能充放电策略、备用容量安排这些问题单点预测根本不够用——我们需要一批“可能发生的功率曲线”每一批曲线都意味着一种可能的世界。基于蒙特卡洛的风电功率/光伏功率场景生成方法干的就是这件事用大量随机采样模拟风速、辐照度的随机过程再经过风电机组功率曲线、光伏组件电气模型生成成千上万条24小时出力曲线供后续做场景缩减、随机优化和风险评估。这篇文章不讲教材上的抽象定义直接以Matlab实现为线索从原理、建模、写码到校验完整走一遍把中间踩过的坑也一并说清楚。适合正在做新能源功率预测、微电网调度或者储能容量配置的朋友参考。1. 场景生成要解决的核心问题与整体设计1.1 为什么不是“预测”而是“场景”我先说个例子。假设明天中午的光伏预测功率是80MW调度部门看了这个数觉得并网没问题就让常规机组少发一点。结果实际来了一片云功率从80MW掉到30MW系统频率短时间内就开始波动AGC机组来不及补只能切负荷。这种案子在光伏占比高的区域不是没发生过。问题出在哪儿出在“预测只有一条线”这件事上。预测值本质是条件期望它把所有可能性平均掉了但优化调度关心的不是平均值而是“哪些情况下我会出问题”。蒙特卡洛场景生成的思路就是把这个过程反过来我们先模拟出一大批可能的天气情况每一个情况都得到一条风电或者光伏的功率曲线再把这一批曲线作为输入送进优化模型。这样优化器不仅看到“明天可能出80MW”还会看到“有可能只有30MW也有可能110MW”。它能提前为这些可能性分配备用、安排储能充放电甚至调整检修计划。这就是为什么在随机优化、鲁棒优化、模型预测控制这些方向上场景几乎是标配输入。这里特别提醒一点场景生成和预测不是互斥关系。很多工程落地是先用一个数值天气预报或者统计模型得到一个基准预测曲线再以基准曲线为中心叠加蒙特卡洛随机误差生成一组场景。这种方式看起来和从物理过程直接采样不一样但本质上都在表达“未来可能的出力集合”。如果项目里已经有成熟的预测模型没必要推倒重来把预测误差分布改成蒙特卡洛抽样即可比从零建一个风电物理模型快得多。1.2 蒙特卡洛方法为什么适合干这件事蒙特卡洛这个名词听起来高深本质上就是“大量随机抽样”。你无法准确预测明天每小时的精确风速但有历史数据就能拟合出风速的概率分布从分布里抽样一万次每次当成一个“虚拟的明天”就得到一万种可能的功率曲线。单看其中一条可能都不具备任何意义但这一万条放在一起它们的统计特征比如均值、方差、分位数、极值、时序相关性就近似刻画了真实世界的随机性。做这件事的核心是“生成”而不只是“抽样”。如果我们把每个小时的风速都独立地从分布里抽一次得到的曲线会非常跳一会儿满发一会儿零出力这在现实中不可能。实际的风速有持续性今天10m/s明天大概率不会变成2m/s光伏上空的云也不是每五分钟独立上演一台戏。所以在抽样之外还要把时间相关性、空间相关性这些物理约束加进去。蒙特卡洛只是一个框架质量高低取决于你在这个框架里塞进去的有多少合理的物理模型。我自己比较反感那种“随机数一把梭”的做法先跑几千次再看看效果不行就调随机数。正确姿势是先想清楚不确定性从哪里来。风电的不确定性主要在风速的随机波动光伏的不确定性主要在云遮过程两者对应的概率分布和时间相关结构完全不同混用一套模板必然出问题。这也是下文把风电和光伏拆成两套函数来讲的原因。1.3 标准流程是怎样的任何一个场景生成项目不论风电还是光伏流程大体一致先收集历史数据和基础参数然后估计随机变量的概率分布接着开始随机抽样并叠加时序相关性再通过设备的功率转换模型得到功率曲线最后对生成的场景做质量校验如果场景太多算不动还要做场景缩减。这个流程听起来不复杂但每一步都有细节第4节我会重点讲校验前面这些环节慢慢展开。从工程角度我习惯把整个流程拆成两个阶段离线标定阶段和在线生成阶段。离线标定做一次就够了包括拟合分布参数、估计自回归系数、整理设备功率曲线在线生成阶段则要反复跑配合不同初始状态或不同预报值生成当前场景集。代码结构上也是按这个逻辑拆分能避免每次运行都重新拟合参数省下不少时间。2. 项目准备与Matlab整体架构2.1 需要哪些输入数据和参数场景生成不是凭空跑几千次随机数就算完事建模质量首先取决于输入数据。我在做项目时通常会准备四类东西。第一类是风速或辐照度的历史观测数据。风速一般以10分钟或1小时为间隔辐照度常用GHI水平总辐照度光伏项目还要准备气温数据。第二类是设备的功率转换参数风电机组要切出风速、切入风速、额定风速、额定功率光伏要组件面积、转换效率、温度系数、额定容量。第三类是地点的天文参数光伏必须知道纬度和当天日期用来计算太阳高度角和晴空辐照度。第四类是分布参数可以直接从历史数据拟合也可以用典型值比如风速Weibull分布的尺度参数和形状参数。我整理了一个常用的参数表方便对照以下数值来自常见工程场景实际项目要按对象重新标定参数含义示例值Pr风电机组额定功率2 MWcut_in切入风速3 m/srated_ws额定风速12 m/scut_out切出风速25 m/sscale风速Weibull尺度参数7.5shape风速Weibull形状参数2.2eta光伏组件转换效率0.18A光伏板总面积5000 m²temp_coeff温度系数-0.003 /°Cbeta_alpha云遮系数Beta分布α8beta_beta云遮系数Beta分布β2参数估计算不上复杂但有个坑要提醒直接从历史风速用矩估计或极大似然法拟合Weibull参数时样本里如果有大量静风时段形状参数会被拉得偏小生成的场景会频繁出现低出力。合理做法是先按季节或月份分开拟合同时结合风机实际出力范围把不合理的极端值截断处理。光伏的Beta分布参数也一样分晴天、多云、阴天的典型天气分别拟合效果远比一个全年统一值好。2.2 代码整体分几层我的Matlab实现习惯上分成三层顶层是场景批量生成入口中间是风电场景生成和光伏场景生成两个核心函数底层是晴空模型、功率曲线换算这类小函数。这样做的最大好处是设备参数调整只改底层函数分布参数调整只改中间层批量出图、做校验的代码不用跟着动。顶层入口大致长这样% 生成前先把参数装到结构体里 windPrm.seed 2025; windPrm.P_rate 2e3; % kW windPrm.cut_in 3; windPrm.rated_ws 12; windPrm.cut_out 25; windPrm.scale 7.5; windPrm.shape 2.2; windPrm.ar 0.85; solarPrm.seed 2025; solarPrm.eta 0.18; solarPrm.A 5000; % m2 solarPrm.T_ref 25; % ℃ solarPrm.temp_coeff 0.003; solarPrm.beta_alpha 8; solarPrm.beta_beta 2; solarPrm.c_ar 0.9; nS 1000; % 场景数 T 24; % 预测时长单位小时 windScen gen_wind_scenarios(nS, T, windPrm); solarScen gen_solar_scenarios(nS, T, solarPrm, 75, 32.0, 118.8);注意单位要统一。我在实际项目中吃过亏风速给的单位是m/s功率曲线坐标也是m/s结果代码里某处偷偷乘了3.6换算成km/h所有场景集体失真画图时还看不出来一送进优化模型就发现备用容量全部算偏。工程代码里最重要的不是花哨算法而是单位一致和维度注释。2.3 随机种子与可复现性Monte Carlo的一个陷阱是每次运行结果会变。有人用这个特性来“证明”算法厉害其实很不科学。做研究、写报告、调模型最怕就是前后两次跑的结果对不上。我的习惯是每个函数入口都固定随机种子用一个prm.seed字段传进去。Matlab里写成rng(prm.seed);如果需要并行还要注意parfor worker上的随机流问题。直接用rand或randn每个worker会得到独立子流总体可复现性一般还能保证但在老版本Matlab里不同worker分配方式可能导致结果对不上。稳妥办法是用RandStream设置全局或Worker级别的随机流s RandStream(mlfg6331_64, Seed, prm.seed); RandStream.setGlobalStream(s);这个细节大多数教程不会提但遇到并行运行不一致时真能省下半天排查时间。2.4 可视化与中间量检查场景生成的调试不能光看最终功率曲线中间量非常关键。我每次跑完都会画三张图第一张家图展示场景集中若干条功率曲线第二张画历史数据的平均功率和场景集平均功率的对比第三张画分位数带也就是0.05、0.5、0.95三条曲线形成的区间。如果只画第一张你很可能只看到一堆乱线很难判断统计属性是否正确。figure; hold on; plot(windScen(1:50,:), Color, [0.7 0.7 0.7]); plot(mean(windScen,1), k-, LineWidth, 2); plot(prctile(windScen, 5, 1), b--, LineWidth, 1.5); plot(prctile(windScen, 95, 1), r--, LineWidth, 1.5); legend({单条场景, 均值, P5下边界, P95上边界}, Location, best); xlabel(小时); ylabel(功率/kW); grid on;中间量检查建议把风速序列、云遮系数、晴空辐照度各自画成散点图或曲线图。某个中间量出现负值、突变、超出物理范围多半是模型代码的问题而不是随机数的问题。3. 核心实现风电与光伏场景生成3.1 风电场景Weibull风速加时间相关性风电场景的核心是风速模型。风速长期统计可以用Weibull分布来描述概率密度函数是[ f(v)\frac{k}{\lambda}(\frac{v}{\lambda})^{k-1}e^{-(v/\lambda)^k} ]其中(\lambda)是尺度参数(k)是形状参数。这两个参数决定风速的平均水平、分布胖瘦和极端高风速出现的概率。Matlab里对应的抽样函数就是wblrnd给定λ和k就能直接生成一个风速样本。但如果每个小时都独立抽样曲线会像心电图。解决办法是加一个一阶自回归项让当前时刻风速与上一时刻有关。常用形式是[ v(t) \alpha v(t-1) (1-\alpha) w(t) \epsilon(t) ]其中(w(t))是新抽的Weibull样本(\epsilon(t))是一个小的正态扰动(\alpha)是持续性系数取值一般在0.7到0.95之间。(\alpha)越大风速序列越平滑越小则越毛躁。(\alpha)可以用历史风速数据的一阶自相关系数来估计更简单的做法是直接拿历史时序数据跑一遍var(1)回归。功率则通过风机功率曲线来换算。下面这个函数是我常用的风电场景生成实现function scen gen_wind_scenarios(nS, T, prm) rng(prm.seed); scen zeros(nS, T); for i 1:nS v zeros(1, T); v(1) max(0, wblrnd(prm.scale, prm.shape)); for t 2:T v(t) prm.ar * v(t-1) (1 - prm.ar) * wblrnd(prm.scale, prm.shape) 0.1 * randn; v(t) max(0, v(t)); end for t 1:T if v(t) prm.cut_in || v(t) prm.cut_out p 0; elseif v(t) prm.rated_ws p prm.P_rate * (v(t) - prm.cut_in) / (prm.rated_ws - prm.cut_in); else p prm.P_rate; end scen(i, t) p; end end end这版代码把“风速抽样”和“功率曲线换算”放在同一个循环里便于阅读。风机的功率曲线未必是线性的现实中很多机组的曲线更接近S形或者分段多项式。如果你手头有制造商提供的离散功率曲线点完全可以改成interp1做插值。我遇到过有人直接拿三次样条拟合整条曲线结果在切出风速附近出现过冲低压穿越中间出现负功率这种物理上不可能的数值一定要人工检查。对于超短期风电功率预测场景比如未来4小时每15分钟一个点时间粒度更细风速的持续性系数会更高因为15分钟尺度上风速变化不会太剧烈。这里T不再是24而是16个点AR(1)系数可能需要重新从对应粒度的历史数据里估计不能直接拿小时级参数硬套。3.2 光伏场景晴空模型加云遮系数光伏场景生成比风电多一个环节除了随机性本身还有一个确定性很强的“晴空辐射”基准。光伏功率有天上的太阳几何约束晚上肯定是零中午大概率最大这个趋势是确定的随机性主要体现在云遮挡上。所以我的做法分成两步先算晴空辐照度再乘一个0到1之间的云遮系数。云遮系数可以用Beta分布抽取因为Beta分布定义在0到1区间能刻画多云和晴天的不同形态。典型参数组合晴天多的地方Beta分布的均值会偏高比如参数α8、β2均值0.8阴天多的地方α和β调整到均值低于0.5。晴空模型有简单有复杂。教科书上的Hottel模型以及Duffie和Beckman书里的公式会考虑海拔、大气质量、太阳天顶角等多个因素。我这里给出一个工程常用的简化模型按小时逐点计算太阳高度角function G clear_sky_model(dayOfYear, lat, lon, T) % 简化晴空模型单位W/m2逐小时 G zeros(1, T); for h 1:T decl 23.45 * sin(2 * pi * (284 dayOfYear) / 365); hourAngle 15 * (h - 12); cosZenith sin(decl * pi / 180) * sin(lat * pi / 180) ... cos(decl * pi / 180) * cos(lat * pi / 180) * cos(hourAngle * pi / 180); cosZenith max(0, cosZenith); G(h) 1000 * cosZenith; % 简化值没有做工随海拔修正 end % 夜间、凌晨极低辐照度直接置零避免后续出现微小正功率 G(G 10) 0; end这个模型只用到纬度、日期和小时数不要求气象再分析数据工程演示已经够用。如果你要做高精度光伏并网研究可以再替换成能考虑海拔和大气透明度的模型甚至直接读NASA或当地辐照度表的晴空数据。但注意替换后必须保证输出的G是“理论晴空值”否则后面云遮系数的物理意义就乱了。接下来是光伏场景生成函数。云遮系数同样用AR(1)过程做序列相关模拟云团在一定时间尺度上的连续遮挡。光伏功率用线性简化模型计算[ P(t) \eta \cdot A \cdot G(t) \cdot [1 - \beta (T_c - T_{ref})] ]其中(\eta)是组件转换效率(A)是组件面积(\beta)是温度系数(T_c)是组件温度。组件温度本身也有不确定性这里用一个围绕参考温度的随机过程来模拟。完整函数function scen gen_solar_scenarios(nS, T, prm, dayOfYear, lat, lon) rng(prm.seed); G0 clear_sky_model(dayOfYear, lat, lon, T); scen zeros(nS, T); for i 1:nS c zeros(1, T); c(1) betarnd(prm.beta_alpha, prm.beta_beta); for t 2:T c(t) prm.c_ar * c(t-1) (1 - prm.c_ar) * betarnd(prm.beta_alpha, prm.beta_beta) 0.03 * randn; c(t) min(1, max(0, c(t))); end Tmp prm.T_ref 8 * randn(1, T); % 简化的组件温度波动 for t 1:T if G0(t) 0 scen(i, t) 0; else p prm.eta * prm.A * G0(t) * c(t) * (1 - prm.temp_coeff * (Tmp(t) - prm.T_ref)); scen(i, t) max(0, p); % 功率不可能为负 end end end endBeta分布抽样出来的云遮系数是随机的裸用AR(1)可能让系数跑出[0,1]区间所以代码里做了min和max裁剪。这里的裁剪要小心如果裁剪量过大会破坏预设的Beta分布形状场景功率分布和真实偏差会越来越大。解决方向要么增大AR(1)的新息扰动权重要么直接把随机过程建在Logit变换后的空间上生成完再映射回[0,1]。超短期光伏功率预测的常见应用场景里云遮过程的时效性很强可能一块云的遮挡周期只有十几分钟。如果T是4小时、15分钟粒度AR(1)系数需要根据场地历史实测重新估计。我踩过的坑是照搬小时级参数到15分钟级生成出来的场景要么云遮变化慢得像静止要么极端频繁地闪变跟实际完全不搭。3.3 批量生成与并行化场景数量一上来比如一次生成5000条、每条96个时段两层循环跑下来Matlab可能一小段时间出不了结果。这时候有两个提速方向。第一个是向量化。把nS维循环尽量改成矩阵运算。光伏场景尤其适合向量化因为晴空模型是固定的随机扰动可以生成nS×T的矩阵一次性算完。第二个是并行化这是大场景最有效的方案。把普通for改成parfor前先确认Parallel Computing Toolbox已安装并执行parpool;这里有个从坑里爬出来的经验parfor和rng的搭配要格外小心。我在一次项目里用parfor跑5000场景笔记本16核全部拉满速度确实快但连续两次运行生成的场景居然不一样。后来发现是随机数流的问题。彻底解决的办法是用RandStream设定可复现随机流再配合parfor的worker设置保证运行结果可复现。如果只是做参数粗筛临时用parfor没问题但要出报告、提交论文还是建议串行加固定种子或者把随机数流管理代码一并写进去。提速效果可以用简单表格评估场景数串行耗时参考8核parfor耗时参考1000约18秒约4秒5000约85秒约15秒10000约170秒约30秒具体数据因机器和代码实现差异很大但至少说明并行化在场景数多时是值得的。我见过不少人在100条场景时也开并行其实启动并行池本身就花时间小场景下得不偿失。4. 场景质量校验与后处理4.1 场景数量选多少合适场景数不是越大越好。太少极端情况没覆盖到太多后续随机优化模型解起来慢。工程上有个经验判断先跑场景数为50、100、200、500、1000的几组实验对比生成的功率均值曲线和分位数包络。当组间统计量变化已经很小就说明场景数够了。具体操作可以在Matlab里循环生成不同规模场景集然后检查某几个关键时段的分位数qLevels [0.05, 0.5, 0.95]; for n [50, 100, 200, 500, 1000] s gen_wind_scenarios(n, T, windPrm); q prctile(s, qLevels * 100, 1); fprintf(n%4d, 第12小时: P50%.1f, P05%.1f, P95%.1f\n, ... n, q(2,12), q(1,12), q(3,12)); end我自己项目的常见规律是风速场景用到300到800条已经能稳定光伏场景因为云遮的极端事件稀疏偶尔需要到1000条以上尤其是研究极端天气下的备用缺口时。再多的时间优化模型规模就会几何级上涨此时应该用场景缩减而不是继续加码。4.2 用统计指标验质量肉眼看好不好看是一回事统计上过不过关是另一回事。场景集的分布必须和历史数据或者高精度预测结果的统计特征一致不然优化模型会学到一套“假世界”。我常用四个校验点均值曲线是否接近历史日均值分位数包络是否包裹了大部分历史曲线每个时刻的方差是否在一个量级时序相关性比如相邻时刻功率差的一阶自相关是否与历史接近。写成代码也就是prctile、mean、std和autocorr几个函数的事。如果想再做严格一点可以算连续等级概率评分CRPS对比场景预测相对历史观测的累积分布表现。CRPS在气象和新能源领域很常用它同时惩罚了不准确和过度自信的场景集。只是Matlab自带没有现成函数需要自己按公式写或者下载社区实现。我这里重点说一个容易忽略的问题只看均值不叫校验。我见过一个案例生成的场景集均值曲线和真实数据几乎完全重叠但是每个场景内部时序完全是白噪声没有任何持续性。送进储能优化模型后模型发现每时每刻都用最便宜策略就行最终方案和真实运行成本差了30%。所以校验时序相关性特别重要比均值对齐更关键。4.3 场景缩减从1000条到10条一个随机优化模型如果塞进1000个场景每个场景又有96个时段变量规模会非常庞大求解时间可能从小时变成天。此时需要场景缩减目标是用更少的场景近似原始场景集的概率分布。最常用的是快速前向选择法也称后向削减法基本思路是迭代地删掉对分布“贡献最小”的场景并把被删场景的概率权重累加到距离它最近的保留场景上。实际工程里不少人也直接用K-means聚类先把所有场景聚成10到20类再用每类中心场景代表那一类聚类占比作为场景权重。K-means实现代码很短Matlab有内置kmeans函数[idx, C] kmeans(scen, 10); weights histcounts(idx, 10) / nS;这里有个要记住的点场景缩减应该在所有场景生成验证之后再做而不是一上来就缩减。缩减本质是“有损压缩”如果最初场景集质量就不行缩减完的少数场景会把误差放大。先保证原始集合质量再考虑少化顺序不能颠倒。5. 常见问题与排查技巧实录5.1 生成的场景曲线过于平滑或过于激烈如果场景和真实曲线相比太平均化多半是持续性系数α设置太高或者分布参数被平均掉了。如果曲线剧烈跳变则可能是α太低或者AR(1)模型的新息方差太大。排查时先别改模型把场景画出来和历史曲线放一起看再写一小段代码算出相邻时段功率差的标准差和历史数据对比调整就有方向了。还有一种情况是生成的场景全部集中在一个很小的区间里看不出不确定性。这往往是分布参数拟合得过窄比如用了一个月的数据而这个月恰好天气稳定得到的Beta分布方差极小。解决办法是拉长历史数据窗口至少覆盖一个季度或者直接按不同天气类型分层处理。5.2 随机数种子影响实验结果场景生成的结果本质上带有随机性不能用一次运行就下结论。跑实验前固定一个主种子再在参数变化时保持种子一致。报告里明确记录“种子2025场景数500”。如果发现不同种子下结果差异巨大说明场景数不够不是运行出错了要回到第4节检查收敛性。我不止一次看到有人用不同种子对比两种调度策略结果差异其实是噪声差点得出错误结论。5.3 parfor运行慢或结果不一致如果parfor没有加速查一下是否开了并行池。Matlab里parpool后任务管理器里看CPU使用率如果一直只有一两个核在工作多半是parfor循环内部有全局变量或者文件写入冲突导致Matlab退化成串行执行。结果不一致的问题大概率出在随机数流用RandStream锁worker的随机流可以解决。另外提醒一句场景生成函数里如果调用的是自定义函数parfor会要求这些函数都在当前路径可见或者正常运行不会报错但在远程worker里找不到函数。提前把所有函数放到同一个文件夹并addpath设置好能省很多麻烦。5.4 中文注释乱码和路径问题Matlab版本中文注释乱码是老问题。Matlab 2023之后默认编码变化打开旧脚本时中文变乱码多半是文件编码和编辑器默认编码不一致。在编辑器里改UTF-8或者GBK统一一下就行别在脚本里混存两种编码。另外有人把生成的.mat数据直接存在中文路径下到别的机器上跑不出来了建议项目路径用纯英文数据文件名不要带空格和特殊符号。5.5 功率出现负值或夜间不归零风电场景要检查切出风速以上功率是否置零光伏场景要检查夜间G00时是否有微小正功率。出现负值基本可以猜到是温度修正项在低温时算出负值或者风机功率曲线的多项式拟合出了问题。处理办法是在功率输出函数出口统一加一句max(0,p)但更重要的是搞懂物理来源不能只靠裁剪掩盖模型错误。5.6 数据质量差导致场景失真这个最隐蔽。新建场站没有长期历史数据只能用周围气象站的数据替代但海拔和地形差异会导致风速、辐照度的分布都不同。如果直接拿来拟合数据分布偏移场景生成再花哨也没有意义。我的建议是先用一小段实测数据做分布校准至少把Weibull参数的均值误差控制在可接受范围再进入场景生成流程。另外历史数据里要处理缺测和异常跳变否则拟合出来的分布会混入噪声。做这个项目的过程中我最大的体会是蒙特卡洛场景生成不是把随机数往函数里一丢就完事。它最花精力的地方其实在参数标定和质量校验风速的Weibull分布到底该按季节拆还是按季度拆云遮系数的Beta参数怎么从历史数据反推场景集和真实数据的相关性差异怎么修正这些问题在一个标准教程里往往找不到统一答案只能是每个项目重新过一遍。我现在的固定工作流是先花半天把数据处理干净再花半天把几个关键统计量跑出来然后才开始真正调模型算下来反而比一开始就闷头写代码快得多。分享一个小技巧把所有场景生成的中间变量比如风速样本、云遮系数、晴空辐照度都单独存成变量并做内存检查。哪怕最终功率曲线看起来正常只要中间变量有一点点物理异常你也能顺藤摸瓜找到源头。这些中间量就像医生的化验单别只信最终那张诊断结论。