新能源出力建模：Weibull与Beta分布参数估计及Matlab实现

做新能源出力建模这行，绕不开“风电用Weibull，光伏用Beta”这对组合。我第一次在某个微电网规划项目里被要求把风光出力用概率模型描述出来时，第一反应是拿正态分布硬套，结果拟合出来的曲线和实测数据差得离谱，后来才慢慢理解，这两个分布背后对应的是两种完全不同的物理过程——风速的正偏态长尾特性和辐照度在单位区间内的随机波动。这篇博文就围绕这个“组合”展开：为什么偏偏是这两个分布、参数怎么估计、Matlab代码怎么落地、常见坑怎么排，以及这套模型到底能接哪些实际应用。适合正在做新能源并网研究、微电网容量配置、储能调度或者概率潮流的学生和工程师参考，不用重新翻一堆论文，照着思路走就能出结果。

1. 为什么风电选Weibull而光伏选Beta：先从物理意义说起

1.1 风速的正偏态和长尾特性，决定了Weibull是天然选择

风速数据有一个很明显的特点：大量时间风速集中在较低区间，偶尔出现阵风或强风，分布曲线呈现出明显的右偏和长尾。均值附近并不是最常出现风速的地方，而是略低于均值的某个风速段出现频率最高。正态分布描述这种数据时会出两个问题：一是正态分布允许负值，而风速物理上不可能为负；二是正态分布的尾部衰减太快，难以刻画强风事件的概率。Weibull分布的定义域恰好是[0, +∞)，形状参数k可以灵活调整偏态程度，尺度参数c则决定了风速量级。特别地，当k=2时退化为瑞利分布，这也是风电行业最常用的简化模型。

从概率密度表达式看，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 < 1

a和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)是尺度A，param(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); % 防止归一化后恰好出现1（Beta密度在端点奇异） 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或1，Beta拟合崩塌

光伏辐照度数据很大比例集中在0附近（阴天）或1附近（晴空），归一化后直方图两端翘得老高。这种情况下矩估计算出来的a或b可能小于1甚至小于0，Beta密度形态完全失控。我的处理经验是：对有效辐照度样本做“先归一化、后截断”，把大于0.995和小于0.005的值强制收敛到边界内，再重新估计。虽然会损失一点点尾部信息，但参数稳定性和后续优化计算的收敛性提升明显。

4.4 风速转功率的“平台期”问题：别忽略额定风速到切出风速段

很多文章的风速-功率转换用简化的三段式，但处理细节很容易出错。额定风速到切出风速这一段，风机出力不再是三次方关系，而是保持在额定功率的平台期。有同学直接把v^3的关系式作用到所有大于切入风速的段，导致高风速段出力无限上涨，场景里出现2倍、3倍额定功率的荒唐结果。正确做法是先对风速分区间标记，再分段计算，变量命名上建议用idx_region1和idx_region2这种带区域编号的名字，避免逻辑混乱。

4.5 实测问题速查表

5. 进阶扩展：这套组合模型还能往哪走

5.1 从独立组合到相关组合：用Copula构造联合分布

独立组合的最大缺陷是忽视了风资源和光照之间可能存在的相关性。如果想更精细地刻画“大风天往往多云、辐照偏低”这类客观规律，可以在两个边缘分布之上加一个Copula结构，常用的有Frank Copula和Clayton Copula，分别对应正相关和负相关的尾部特征。实现上不需要重新拟合边缘分布，只需要把历史风速和辐照度的样本投影到均匀空间，估计Copula参数，再反过来生成相关均匀样本，最后用Weibull和Beta的逆CDF映射回原始空间。

5.2 从概率模型到优化模型：给储能配置和容量规划用

拟合出参数不是终点，把参数化的概率模型接入优化问题才是价值所在。储能容量配置里最常见的方法是：用蒙特卡洛生成全年8760小时的风光出力场景，叠加负荷曲线后，对每个时段计算功率缺额和富余量，再统计亏电概率和弃电率，最后用这些指标做容量优化。我实际做过的项目里，Weibull+Beta组合模型生成的场景集参与优化，得到的储能容量配置比用平均日出力估算的方法节省约15%的装机成本，这还是在忽略时序相关性的前提下。

5.3 从静态分布到时序建模：日尺度波动模式的引入

Weibull和Beta描述的是长期统计特性，丢掉了时间先后顺序。如果想研究风光出力的爬坡率、日内波动特征，就需要引入自相关结构，比如加一个Markov状态切换，把天气划分为晴、多云、阴雨几种状态，每种状态下用不同的Beta参数建模。这样从“静态分布组合”升级为“隐状态时序组合”，复杂度上升一个台阶，但模型的表现力也强得多。这个方向是后续比较值得投入的研究点。

最后再分享一点我个人的体会。做这种概率建模，千万别一上来就把数据丢给函数库跑拟合。先花一半时间做数据清洗、画直方图、看物理规律，再用代码验证三个数：形状参数是否在常识范围、均值与实测均值差多少、极端场景生成是否符合直觉。这三关过了，模型基本就能用了。代码写得多花哨不如数据理解得透，拟合这活儿啊，七分靠清洗，三分靠算法。