做新能源出力建模这行绕不开“风电用Weibull光伏用Beta”这对组合。我第一次在某个微电网规划项目里被要求把风光出力用概率模型描述出来时第一反应是拿正态分布硬套结果拟合出来的曲线和实测数据差得离谱后来才慢慢理解这两个分布背后对应的是两种完全不同的物理过程——风速的正偏态长尾特性和辐照度在单位区间内的随机波动。这篇博文就围绕这个“组合”展开为什么偏偏是这两个分布、参数怎么估计、Matlab代码怎么落地、常见坑怎么排以及这套模型到底能接哪些实际应用。适合正在做新能源并网研究、微电网容量配置、储能调度或者概率潮流的学生和工程师参考不用重新翻一堆论文照着思路走就能出结果。1. 为什么风电选Weibull而光伏选Beta先从物理意义说起1.1 风速的正偏态和长尾特性决定了Weibull是天然选择风速数据有一个很明显的特点大量时间风速集中在较低区间偶尔出现阵风或强风分布曲线呈现出明显的右偏和长尾。均值附近并不是最常出现风速的地方而是略低于均值的某个风速段出现频率最高。正态分布描述这种数据时会出两个问题一是正态分布允许负值而风速物理上不可能为负二是正态分布的尾部衰减太快难以刻画强风事件的概率。Weibull分布的定义域恰好是[0, ∞)形状参数k可以灵活调整偏态程度尺度参数c则决定了风速量级。特别地当k2时退化为瑞利分布这也是风电行业最常用的简化模型。从概率密度表达式看Weibull分布为f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)v ≥ 0其中k无量纲c的单位与风速一致m/s。k值越小分布越偏低风速占比越高k值越大风速分布越集中曲线越接近正态。实际风电场的数据拟合下来k通常在1.8到2.5之间c在6到9 m/s之间这意味着大多数时候风速集中在3到10 m/s区间恰好是风机切入风速到额定风速附近的工作区间。1.2 辐照度的单位区间特性让Beta分布成为光伏建模首选光伏出力建模和风电最大的不同在于物理量纲。风速可以无限增大但辐照度受大气层顶太阳常数限制经过大气衰减后到达地面的值天然落在一个有界范围内。工程上将实测辐照度除以基准值一般取1000 W/m²或历史最大值后得到一个0到1之间的归一化变量x。这时问题就变成了如何选择一个定义在[0,1]上的分布来描述x的变化规律。Beta分布恰好满足这个要求。其概率密度为f(x) x^(a-1)(1-x)^(b-1) / B(a,b)0 x 1a和b两个形状参数给了它极大的灵活性。a和b都大于1时是单峰钟形a小于1且b大于1时呈反J形适合描述多云天气下低辐照度频发的场景a大于1且b小于1时呈J形对应晴空少云、高辐照度占比高的场景a和b都小于1时则是U形描述极端天气分布。光伏电站的实际辐照度数据往往是“晴空高值”和“阴天低值”两个状态交替出现Beta分布能比正态分布更好地包住这种双峰偏向的形态。1.3 “组合研究”究竟在组合什么独立边缘分布与联合场景生成标题里的“组合”有两层含义。第一层是工程仿真层面最常见的用法把风电出力用Weibull抽样得到光伏出力用Beta抽样得到两者作为两个独立的随机变量一起输入到概率潮流、容量配置或储能调度模型中。这种组合不关心风资源和太阳辐射之间的相关性实现简单适用于规划阶段的粗评估。第二层是统计建模层面的“组合”即考虑风资源和辐照度之间的相关结构。一个地区风速大时往往伴随云层移动加快辐照度波动也会加剧两者并非完全独立。如果想精细刻画风光互补特性需要在两个边缘分布之上引入Copula函数或者联合分布构造。本文主体以第一种独立组合为主——这也是大多数场景下够用的方案第二种进阶做法放在第5章讨论。2. 参数估计的两种核心路径极大似然与矩估计2.1 Weibull参数的极大似然估计公式推导与Matlab代码Weibull参数估计最常用的方法是极大似然估计。对一组观测风速v1, v2, ..., vn构造似然函数后形状参数k满足如下隐式方程1/k Σ(vi^k * ln(vi)) / Σ(vi^k) - Σln(vi)/n尺度参数c则直接由k决定c (Σ(vi^k) / n)^(1/k)实际写代码时没人手算这个迭代Matlab里两条路都可以走。最快的是直接用wblfit函数% 假设wind_data是列向量单位m/s [param, pci] wblfit(wind_data); c_fit param(1); % 尺度参数注意顺序 k_fit param(2); % 形状参数这里必须提醒一下很多教程把wblfit返回顺序搞反param(1)是尺度Aparam(2)是形状B不是我们口语里常说的k在前。如果你用fitdist用法是这样的pd fitdist(wind_data, Weibull); k_fit pd.B; % 形状 c_fit pd.A; % 尺度fitdist返回的对象属性是A对应尺度、B对应形状代码里注释写清楚省得下次翻车。还有一种工程上更实用的拟合思路不是直接对风速密度做拟合而是对风速的三次方加权做拟合。因为风机功率近似与风速立方成正比低风速区的拟合误差对发电量影响远小于高风速区。做法是构造目标函数让拟合曲线在风速立方空间中与历史数据的偏差最小。我实际项目中试过这种方法得到的k、c算出来的年发电量比普通极大似然结果误差小3%到5%。对小风机会差异更明显。2.2 Beta参数矩估计四两拨千斤的闭式解Beta分布的参数估计可以用极大似然但Matlab自带的betafit有时候迭代不稳定而且返回值顺序容易搞混。更直接的办法是矩估计——用样本均值和样本方差反解a和b。假设归一化后的辐照度样本为x1, x2, ..., xn计算样本均值m和样本方差var那么a m * (m*(1-m)/var - 1) b (1-m) * (m*(1-m)/var - 1)代码写出来只有三行m mean(x_data); v var(x_data, 1); % 注意用总体方差 a_fit m * (m*(1-m)/v - 1); b_fit (1-m) * (m*(1-m)/v - 1);矩估计的物理意义很直观均值决定了分布峰值的位置方差决定了分布展宽。辐照度数据如果晴天多均值偏大a会大于b曲线峰值靠近1阴天多则反之。需要注意的是当样本均值接近0或1时m*(1-m)/v会变得非常大可能导致a、b出现负值或者数量级爆炸。这种情况通常说明数据里混入了大量零值或满值样本需要先做数据清洗而不是直接怪公式。2.3 数据预处理静风、夜间零值和缺失值怎么处理才不出错数据预处理这一步决定了整个拟合工作的成败却恰恰是很多文章一笔带过的部分。风速数据里经常有大量静风时段风速为0如果直接把全序列丢给wblfit曲线在低风速区会严重偏离实际因为Weibull分布理论上可以取任何非负值但大量0值会让似然函数变形。我现在的处理办法是先区分“全时段建模”和“有风时段建模”。做年发电量评估时用全时段但额外引入一个静风概率p0最终风速分布看成“p0概率为0 (1-p0)概率服从Weibull”的混合模型做短期出力分布时直接剔除静风样本只对有效风速区间拟合Weibull。辐照度数据处理更讲究。光伏发电只在白天有意义如果把夜间零值一起纳入Beta拟合a、b会被拉向极端值曲线没法看。我一般先筛选出辐照度大于阈值的样本比如大于10 W/m²再对这些“有效光照”样本做归一化和拟合。归一化基准也别总用1000 W/m²有些高纬度或阴雨地区全年辐照度很少接近1000用这个基准会导致归一化数据全部堆积在0.3以下Beta拟合出来的a、b极为不对称。更稳妥的做法是用历史实测最大值的1.05倍或95分位数的1.1倍做基准让归一化数据在(0,1)区间内铺得开。3. Matlab代码实现从单分布拟合到联合场景生成3.1 风速Weibull拟合与可视化直方图和概率密度曲线对齐完整的风速拟合流程我习惯这么写。先在脚本里读取风速数据剔除缺失值然后分全时段和有风时段分别拟合最后把实测直方图和拟合密度曲线画在同一张图上对比。% 风速Weibull拟合完整示例 rng(2024); % 真实工程中这里用 load(wind_data.mat) 读取数据 v_data 30 * wblrnd(2.1, 7.5, 2000, 1); v_data(v_data 0) []; % 剔除静风保留有效风速 v_eff v_data(v_data 0.5); fprintf(总样本 %d有效风速样本 %d\n, length(v_data), length(v_eff)); % 方式一wblfit phat wblfit(v_eff); c_fit phat(1); k_fit phat(2); % 方式二自定义MLE迭代验证 fun (k) (mean(v_eff.^k .* log(v_eff)) / mean(v_eff.^k) - mean(log(v_eff)) - 1/k); k_mle fzero(fun, 2); c_mle mean(v_eff.^k_mle)^(1/k_mle); % 绘图对比 figure; histogram(v_eff, 40, Normalization, pdf, FaceAlpha, 0.3, EdgeColor, none); hold on; x_axis linspace(0, max(v_eff), 200); plot(x_axis, wblpdf(x_axis, c_fit, k_fit), r-, LineWidth, 2); plot(x_axis, wblpdf(x_axis, c_mle, k_mle), b--, LineWidth, 2); legend(实测直方图,wblfit拟合,自定义MLE拟合); xlabel(风速 (m/s)); ylabel(概率密度); title(风速Weibull分布拟合对比); grid on;我实际跑下来自定义MLE和wblfit结果基本一致差异在小数点后三位所以平时直接用wblfit就行。写自定义迭代主要是为了理解计算原理面试或回答审稿意见时用得上。3.2 辐照度Beta拟合归一化、矩估计和边界保护光伏部分的核心是归一化和边界保护。直接看代码% 辐照度Beta拟合完整示例 % 真实工程中这里用 load(radiation_data.mat) 读取小时级辐照度 s_data 60 * betarnd(1.9, 3.2, 2000, 1); % 模拟数据 % 剔除夜间和极低值 s_valid s_data(s_data 10); % 归一化基准实测最大值1.05倍留出白空间 s_norm s_valid / (max(s_valid) * 1.05); % 防止归一化后恰好出现1Beta密度在端点奇异 s_norm min(s_norm, 0.9999); % 矩估计 mu mean(s_norm); variance var(s_norm, 1); a_fit mu * (mu*(1-mu)/variance - 1); b_fit (1-mu) * (mu*(1-mu)/variance - 1); % 直接用betafit对照 beta_params betafit(s_norm); a_params beta_params(1); b_params beta_params(2); % 绘图 figure; histogram(s_norm, 40, Normalization, pdf, FaceAlpha, 0.3, EdgeColor, none); hold on; x_axis linspace(0.001, 0.999, 200); plot(x_axis, betapdf(x_axis, a_fit, b_fit), r-, LineWidth, 2); plot(x_axis, betapdf(x_axis, a_params, b_params), g--, LineWidth, 2); legend(归一化辐照度直方图,矩估计拟合,betafit拟合); xlabel(归一化辐照度); ylabel(概率密度); title(辐照度Beta分布拟合对比); grid on;这段代码里有两个细节特别关键。第一个是归一化基准选择max(s_valid)*1.05这样所有样本严格落在(0,1)内不会出现落在端点导致Beta密度无穷大的情况。第二个是担保逻辑s_norm min(s_norm, 0.9999)防止个别样本归一化后等于1把a、b参数带偏。这两个坑我在初版代码里都踩过出来的Beta曲线要么右端点翘上天要么参数负得没法解释。3.3 蒙特卡洛联合出力场景生成把两个分布“组合”起来模型参数都拟合好之后真正的应用在于生成风光联合出力的场景集。所谓场景集就是一组风电出力pw和光伏出力pp的组合样本每个样本代表一种可能的运行状态。生成方法很简单——从两个分布里独立抽样然后通过出力转换公式把风速映射成风机出力、把归一化辐照度映射成光伏出力。% 蒙特卡洛生成风光联合出力场景 N 10000; rng(2024); % 独立抽样生成风速和归一化辐照度 v_sample wblrnd(c_fit, k_fit, N, 1); s_sample betarnd(a_fit, b_fit, N, 1); % 风机出力模型vci切入、vr额定、vco切出 P_wind_rated 1.5; % MW vci 3; vr 12; vco 25; P_wind zeros(N, 1); idx_region1 v_sample vci v_sample vr; idx_region2 v_sample vr v_sample vco; P_wind(idx_region1) P_wind_rated * (v_sample(idx_region1).^3 - vci^3) / (vr^3 - vci^3); P_wind(idx_region2) P_wind_rated; % 低于切入或高于切出保持0满足 idx_region1 和 idx_region2 之外的样例 % 光伏出力模型简化为线性忽略温度影响 P_pv_rated 1.0; % MW P_pv P_pv_rated * s_sample; P_pv(P_pv P_pv_rated) P_pv_rated; % 合并为联合出力场景矩阵 scenarios [P_wind, P_pv]; % 查看统计特征 fprintf(风电平均出力: %.3f MW\n, mean(P_wind)); fprintf(光伏平均出力: %.3f MW\n, mean(P_pv)); fprintf(联合出力平均: %.3f MW\n, mean(P_wind P_pv));这段是整套流程里最实用的一块生成出来的场景矩阵可以直接喂给后续的优化模型。比如做储能容量配置时把每一行看成一天的出力基线结合负荷曲线就能算出缺电概率和弃风弃光率。做概率潮流时则把这些场景作为输入统计节点电压和支路功率的分布。3.4 场景验证别让模型输出变成“自说自话”场景生成完必须做验证不然模型的输出只是自说自话。我常做三件事。第一件是边缘分布验证。拿生成的10000个风速样本再画直方图和原始实测风速直方图叠加肉眼判断形状是否吻合。更严谨一点可以用K-S检验[h_w, p_w] kstest(v_sample, CDF, makedist(Weibull, A, c_fit, B, k_fit)); [h_s, p_s] kstest(s_sample, CDF, makedist(Beta, a, a_fit, b, b_fit));p值大于0.05说明生成样本与目标分布在统计上没有显著差异模型可信。第二件是联合场景的合理性检查。把风电出力和光伏出力画成散点图检查是否存在不合理的组合——比如风电满载同时光伏也为零的场景这种组合概率极低如果出现比例过高说明两个分布的抽样逻辑有问题需要考虑引入相关结构。第三件是极端场景的压力测试。人为设置极端气象条件——风速达到切出风速、辐照度长期低于阈值看模型生成的场景是否符合物理预期。这个测试往往能发现风机功率曲线分段函数的边界错误。4. 实操中常见的五个坑与排查速查表4.1 拟合不收敛或参数异常大先查数据而不是查函数有次我帮一个课题组看拟合结果他们用wind数据跑wblfit出来的k值高达6.8完全不符合物理常识。排查了一圈发现数据文件里混入了风速传感器的故障常值——连续几百个小时都显示同一个数。这种“平台式”数据会让似然函数误以为存在一个巨大的概率集中点把分布形状拉得特别尖。处理办法很简单预处理时先看一眼数据的时间序列图把持续不变的异常段剔掉。参数异常大或者异常小的时候九成是数据问题不是算法问题。4.2 直方图和密度曲线对不上问题在归一化方式直方图和拟合曲线不重合是新手最容易困惑的地方。Matlab里histogram默认纵轴是频数而拟合的pdf纵轴是概率密度两者量纲都不一样直接叠图当然对不上。解决办法是画直方图时设置Normalization,pdf参数让直方图总面积归一化成1再叠加概率密度曲线。这一步代码就一行但效果天差地别。4.3 归一化后样本贴近0或1Beta拟合崩塌光伏辐照度数据很大比例集中在0附近阴天或1附近晴空归一化后直方图两端翘得老高。这种情况下矩估计算出来的a或b可能小于1甚至小于0Beta密度形态完全失控。我的处理经验是对有效辐照度样本做“先归一化、后截断”把大于0.995和小于0.005的值强制收敛到边界内再重新估计。虽然会损失一点点尾部信息但参数稳定性和后续优化计算的收敛性提升明显。4.4 风速转功率的“平台期”问题别忽略额定风速到切出风速段很多文章的风速-功率转换用简化的三段式但处理细节很容易出错。额定风速到切出风速这一段风机出力不再是三次方关系而是保持在额定功率的平台期。有同学直接把v^3的关系式作用到所有大于切入风速的段导致高风速段出力无限上涨场景里出现2倍、3倍额定功率的荒唐结果。正确做法是先对风速分区间标记再分段计算变量命名上建议用idx_region1和idx_region2这种带区域编号的名字避免逻辑混乱。4.5 实测问题速查表现象可能原因排查与解决wblfit返回的k值极大5数据混入平台常值或大量重复值检查时间序列剔除非物理的恒定段Beta拟合出的a或b为负归一化数据大量堆在0或1附近调整归一化基准截断端点增加白空间直方图与pdf曲线错位histogram未设置概率密度归一化加Normalization,pdf参数重画场景中出现风电出力远超额定功率风速-功率三次方公式未做分段限制先分段标记再计算不同区间出力蒙特卡洛结果每次跑都不一样随机数种子未固定加rng(2024)固定种子保证结果可复现全时段风速拟合效果差静风时段占比过高拉低形状参数按全时段/有风时段分开建模额外引入静风概率5. 进阶扩展这套组合模型还能往哪走5.1 从独立组合到相关组合用Copula构造联合分布独立组合的最大缺陷是忽视了风资源和光照之间可能存在的相关性。如果想更精细地刻画“大风天往往多云、辐照偏低”这类客观规律可以在两个边缘分布之上加一个Copula结构常用的有Frank Copula和Clayton Copula分别对应正相关和负相关的尾部特征。实现上不需要重新拟合边缘分布只需要把历史风速和辐照度的样本投影到均匀空间估计Copula参数再反过来生成相关均匀样本最后用Weibull和Beta的逆CDF映射回原始空间。5.2 从概率模型到优化模型给储能配置和容量规划用拟合出参数不是终点把参数化的概率模型接入优化问题才是价值所在。储能容量配置里最常见的方法是用蒙特卡洛生成全年8760小时的风光出力场景叠加负荷曲线后对每个时段计算功率缺额和富余量再统计亏电概率和弃电率最后用这些指标做容量优化。我实际做过的项目里WeibullBeta组合模型生成的场景集参与优化得到的储能容量配置比用平均日出力估算的方法节省约15%的装机成本这还是在忽略时序相关性的前提下。5.3 从静态分布到时序建模日尺度波动模式的引入Weibull和Beta描述的是长期统计特性丢掉了时间先后顺序。如果想研究风光出力的爬坡率、日内波动特征就需要引入自相关结构比如加一个Markov状态切换把天气划分为晴、多云、阴雨几种状态每种状态下用不同的Beta参数建模。这样从“静态分布组合”升级为“隐状态时序组合”复杂度上升一个台阶但模型的表现力也强得多。这个方向是后续比较值得投入的研究点。最后再分享一点我个人的体会。做这种概率建模千万别一上来就把数据丢给函数库跑拟合。先花一半时间做数据清洗、画直方图、看物理规律再用代码验证三个数形状参数是否在常识范围、均值与实测均值差多少、极端场景生成是否符合直觉。这三关过了模型基本就能用了。代码写得多花哨不如数据理解得透拟合这活儿啊七分靠清洗三分靠算法。