MATLAB实战：用fitdist函数搞定风光数据Weibull和Beta分布拟合（附完整代码）
MATLAB实战：用fitdist函数搞定风光数据Weibull和Beta分布拟合（附完整代码）
在可再生能源领域，风速和太阳辐照度的概率分布建模是电力系统规划、风机性能评估和光伏发电预测的基础工作。对于工程师和研究人员来说，掌握快速准确的数据拟合方法，意味着能够更高效地完成不确定性分析、场景生成和系统优化等任务。本文将深入探讨MATLAB中fitdist函数在Weibull和Beta分布拟合中的实战应用，提供可直接运行的代码示例，并分享处理实际工程数据时的实用技巧。
1. 理解Weibull和Beta分布在风光数据中的应用
1.1 Weibull分布与风速特性
Weibull分布因其形状灵活而被广泛用于描述风速分布。它的概率密度函数(PDF)为：
f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中：
- v表示风速
- k是形状参数(决定分布形状)
- c是尺度参数(影响分布范围)
实际应用中发现：当形状参数k≈2时，Weibull分布特别适合描述中等风速区域的风况，这也是为什么它成为风电行业的标准模型。
1.2 Beta分布与光伏特性
Beta分布则常用于描述归一化后的太阳辐照度，其PDF为：
f(x) = (x^(a-1) * (1-x)^(b-1)) / beta(a,b)其中：
- x是归一化的辐照度(0到1之间)
- a, b是形状参数
- beta(a,b)是Beta函数
注意：Beta分布要求输入数据必须在(0,1)区间，因此实际应用中需要先对原始辐照度数据进行归一化处理。
2. fitdist函数基础与参数拟合方法对比
2.1 fitdist函数的基本用法
MATLAB的fitdist函数提供了最大似然估计(MLE)方法来拟合分布参数，基本语法为：
pd = fitdist(data, 'distribution_name')对于Weibull和Beta分布，分别使用：
% Weibull分布拟合 weibull_pd = fitdist(wind_data, 'weibull'); % Beta分布拟合 beta_pd = fitdist(solar_data, 'beta');2.2 不同拟合方法的对比
| 方法 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| fitdist(MLE) | 统计性质最优，结果精确 | 计算量较大，可能不收敛 | 数据质量好，样本量大 |
| 矩估计法 | 计算简单快速 | 对小样本不稳健 | 初步分析，快速估算 |
| 标准差法 | 实现简单 | 精度较低 | 需要快速结果时 |
工程经验建议：对于正式分析报告，推荐使用fitdist的MLE方法；而在迭代开发或快速验证阶段，可以考虑矩估计法。
3. 实战案例：完整的风光数据拟合流程
3.1 风速数据Weibull拟合
假设我们有一个月的风速数据，采样间隔为10分钟，数据存储在变量wind_speed中(维度：N×24，N为天数)：
% 初始化参数存储 k_values = zeros(1, 24); c_values = zeros(1, 24); % 逐小时拟合 for hour = 1:24 % 提取当前小时的所有风速数据 current_hour_data = wind_speed(:, hour); % 移除无效数据(如NaN) valid_data = current_hour_data(~isnan(current_hour_data)); % 拟合Weibull分布 if ~isempty(valid_data) pd = fitdist(valid_data, 'weibull'); k_values(hour) = pd.B; % 形状参数 c_values(hour) = pd.A; % 尺度参数 else k_values(hour) = NaN; c_values(hour) = NaN; end end3.2 光伏数据Beta拟合
对于光伏数据，需要特别注意夜间零值处理：
% 假设solar_irradiance是原始辐照度数据(维度：N×24) max_irradiance = max(solar_irradiance(:)); % 获取最大辐照度用于归一化 solar_normalized = solar_irradiance / max_irradiance; a_values = zeros(1, 24); b_values = zeros(1, 24); for hour = 1:24 hour_data = solar_normalized(:, hour); % 移除零值和无效数据 non_zero_data = hour_data(hour_data > 0 & ~isnan(hour_data)); if length(non_zero_data) > 10 % 确保有足够样本 try pd = fitdist(non_zero_data, 'beta'); a_values(hour) = pd.a; b_values(hour) = pd.b; catch % 处理拟合失败情况 a_values(hour) = NaN; b_values(hour) = NaN; end else a_values(hour) = NaN; b_values(hour) = NaN; end end4. 常见问题与解决方案
4.1 数据全相同导致的拟合错误
当输入数据所有值都相同时(如夜间光伏数据全为0)，fitdist会报错。解决方法：
- 数据预处理：过滤掉全相同的数据段
- 异常处理：使用try-catch块捕获错误
try pd = fitdist(data, 'beta'); catch ME if strcmp(ME.identifier, 'stats:probdist:fitdist:AllSameData') disp('所有数据值相同，跳过拟合'); pd = []; else rethrow(ME); end end4.2 拟合结果可视化验证
良好的实践是总是可视化拟合结果：
% Weibull拟合可视化示例 figure; histogram(wind_data, 'Normalization', 'pdf'); hold on; x = linspace(min(wind_data), max(wind_data), 100); pdf_values = pdf(pd, x); plot(x, pdf_values, 'LineWidth', 2); title('Weibull分布拟合验证'); xlabel('风速(m/s)'); ylabel('概率密度'); legend('实际数据', '拟合曲线');4.3 提高拟合稳定性的技巧
- 数据量不足时：考虑使用移动窗口增加样本量
- 参数初始值：对于复杂数据，可以指定初始参数值
- 拟合选项：调整最大迭代次数等优化参数
% 指定初始参数示例(Weibull) options = statset('MaxIter', 1000, 'Robust', 'on'); pd = fitdist(data, 'weibull', 'Options', options, 'StartPoint', [2, 5]);5. 高级应用与性能优化
5.1 并行计算加速大规模数据拟合
对于多年的高频采样数据，可以使用并行计算：
parfor hour = 1:24 % 并行化的拟合代码 pd = fitdist(wind_data(:, hour), 'weibull'); % 存储结果... end5.2 分布参数的时间序列分析
将拟合得到的参数按时间序列分析，可发现潜在模式：
% 分析形状参数k的日变化规律 figure; plot(1:24, k_values, '-o'); title('Weibull形状参数k的日变化'); xlabel('小时'); ylabel('k值'); grid on;5.3 生成符合拟合分布的随机数据
基于拟合结果生成场景数据：
% 生成1000个随机样本 simulated_wind = random(pd, [1000, 1]); % 或者使用更专业的场景生成函数 wind_scenarios = mvnrnd(mean_wind, cov_wind, 1000);