MATLAB实战：用fitdist函数搞定风速与光伏数据的Weibull和Beta分布拟合

MATLAB实战：用fitdist函数搞定风速与光伏数据的Weibull和Beta分布拟合

在可再生能源领域，风速和太阳辐照度的概率分布建模是电力系统规划、风机性能评估和光伏出力预测的基础工作。对于工程师和研究人员来说，掌握MATLAB中fitdist函数的实战应用技巧，能够快速实现Weibull和Beta分布的参数拟合，大幅提升数据分析效率。本文将深入解析这两个分布的特性、fitdist函数的核心参数配置、常见报错解决方案以及结果可视化技巧，帮助读者避开实际应用中的各种"坑"。

1. 理解Weibull与Beta分布在能源领域的应用

1.1 Weibull分布与风速特性

Weibull分布因其形状参数的灵活性，能够很好地描述风速的随机特性。其概率密度函数为：

f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)

其中：

• k为形状参数（shape parameter），决定分布曲线的形态
• c为尺度参数（scale parameter），反映特征风速大小

典型参数范围：

• 陆地风场：k通常介于1.5-2.5之间
• 海上风场：k值可能接近3.0
• 尺度参数c与平均风速正相关

提示：实际项目中，建议先对风速数据进行质量控制，剔除异常值和故障时段记录

1.2 Beta分布与光伏特性

Beta分布特别适合描述太阳辐照度的概率特性，因为：

• 定义域为[0,1]，与归一化辐照度匹配
• 双参数可调，适应不同天气模式

其概率密度函数为：

f(x) = x^(a-1)*(1-x)^(b-1)/beta(a,b)

关键特性：

• 参数a,b > 0控制分布形状
• 当a=b=1时退化为均匀分布
• 日出日落时段易出现数据全零情况

2. fitdist函数深度解析

2.1 基础语法与参数配置

fitdist的标准调用格式：

pd = fitdist(data, distname, 'Name', value)

关键参数选项：

对于Weibull拟合的完整示例：

wind_data = xlsread('wind_2022.xlsx'); pd_weibull = fitdist(wind_data(:,2), 'Weibull',... 'Options', statset('MaxIter',1000,'TolX',1e-6));

2.2 数据预处理技巧

风速数据准备：

1. 清洗异常值（通常风速>25m/s需验证）
2. 处理缺失数据：

   wind_data(isnan(wind_data)) = mean(wind_data,'omitnan');

3. 按时间维度分组（如每小时一组）

光伏数据特殊处理：

% 处理全零数据情况 solar_data(solar_data==0) = eps; % 归一化到(0,1)区间 solar_norm = solar_data / max(solar_data);

3. 实战案例：24小时时序数据拟合

3.1 风速数据批量处理

假设有全年逐小时风速数据矩阵（365×24）：

% 初始化参数存储 k_24h = zeros(1,24); c_24h = zeros(1,24); for hour = 1:24 current_data = wind_matrix(:,hour); current_data = current_data(current_data > 0.5); % 过滤低风速 pd = fitdist(current_data, 'Weibull'); k_24h(hour) = pd.B; c_24h(hour) = pd.A; % 异常值处理 if pd.B > 5 || pd.A > 30 k_24h(hour) = median(k_24h); c_24h(hour) = median(c_24h); end end

3.2 光伏数据特殊处理方案

针对夜间零值问题，推荐分段处理策略：

for hour = 1:24 current_solar = solar_matrix(:,hour); if all(current_solar == 0) a_params(hour) = 1; b_params(hour) = 1; % 设为均匀分布 else current_solar(current_solar==0) = min(current_solar(current_solar>0))/2; pd = fitdist(current_solar/max(current_solar), 'beta'); a_params(hour) = pd.a; b_params(hour) = pd.b; end end

4. 结果验证与可视化

4.1 拟合优度检验

Kolmogorov-Smirnov检验：

[h_wind,p_wind] = kstest(wind_data, 'CDF', pd_weibull); if h_wind == 0 disp('Weibull拟合通过KS检验 (p>0.05)'); else disp(['拟合可能存在问题，p值为',num2str(p_wind)]); end

4.2 专业可视化技巧

多子图对比展示：

figure('Position',[100,100,1200,600]) subplot(1,2,1) histfit(wind_data,50,'weibull') title('风速Weibull分布拟合') subplot(1,2,2) pd_beta = fitdist(solar_norm, 'beta'); x = 0:0.01:1; y = pdf(pd_beta,x); bar(histcounts(solar_norm,20)/length(solar_norm),'FaceAlpha',0.5) hold on plot(x,y,'LineWidth',2) title('光伏Beta分布拟合')

参数时序图：

figure yyaxis left plot(1:24, k_24h, '-o') ylabel('形状参数k') yyaxis right plot(1:24, c_24h, '-s') ylabel('尺度参数c') title('24小时Weibull参数变化') grid on

5. 高级技巧与异常处理

5.1 常见报错解决方案

错误1："Data must be in the open interval (0,1) for the beta distribution"

解决方法：

solar_data = (solar_data - min(solar_data)) / (max(solar_data) - min(solar_data) + eps); solar_data = max(min(solar_data,1-eps),eps);

错误2："Weibull fitting failed to converge"

优化方案：

options = statset('MaxIter',10000,'TolX',1e-8); pd = fitdist(data, 'Weibull','Options',options);

5.2 并行计算加速

对于大规模数据集：

parfor hour = 1:24 % 拟合代码保持不变 end

5.3 结果存储与报告生成

建议将拟合结果保存为结构体数组：

results(1).DistName = 'Weibull'; results(1).Parameters = [k_24h; c_24h]; results(2).DistName = 'Beta'; results(2).Parameters = [a_params; b_params]; save('fit_results.mat','results');