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方法；而在迭代开发或快速验证阶段，可以考虑矩估计法。

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 end

3.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 end

4. 常见问题与解决方案

4.1 数据全相同导致的拟合错误

当输入数据所有值都相同时(如夜间光伏数据全为0)，fitdist会报错。解决方法：

1. 数据预处理：过滤掉全相同的数据段
2. 异常处理：使用try-catch块捕获错误

try pd = fitdist(data, 'beta'); catch ME if strcmp(ME.identifier, 'stats:probdist:fitdist:AllSameData') disp('所有数据值相同，跳过拟合'); pd = []; else rethrow(ME); end end

4.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'); % 存储结果... end

5.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);