搞新能源建模的朋友，估计都绕不开这两件事：怎么描述风力的随机性，又怎么描述光伏的随机性。风电功率和光伏功率的出力曲线都不是线性稳定的，但如果我们把风速、光照强度这类原始气象要素用合适的概率分布刻画清楚，后面做容量配置、可靠性分析、储能优化，全都变得有据可依。风电场出力建模，业内最常用的就是两参数Weibull分布；光伏出力建模，最经典的则是Beta分布。把这两个分布拿到一起做组合研究，再用Matlab实现整套拟合、绘图与统计分析流程，正是这篇博文要解决的完整问题。
这篇内容适合谁？一个是正在做新能源并网、微电网规划类课题的研究生，另一个是从事风光资源评估的前期设计人员。你不需要把数理统计从头学一遍，只要跟着我把分布参数估计的原理弄明白、把代码框架拿走去改数据，基本就能直接上手出结果。我会把Weibull和Beta两个分布从物理背景、参数推导到Matlab代码实现全部拆开讲透，包括拟合后怎么评估好坏、组合研究到底在做什么，以及我实测踩过的那些坑。
1. 内容整体设计与思路拆解
1.1 风速为什么选Weibull而非正态——先把物理直觉讲清楚
很多刚接触风资源分析的人会先入为主地认为风速数据取个平均，画个直方图，长得像钟形就可以用正态分布近似。但实际的风速观测数据几乎全是右偏的，也就是小风速出现的频率很高，大风速出现的频率很低，拖一条长长的右尾。正态分布是对称的，它既不能刻画“微风日数多”这种常见形态，又无法处理风速永远大于等于0这个天然边界。瑞利分布虽然也能描述风速，但它本质上是Weibull分布形状参数k=2的特例，灵活性太差。
Weibull分布之所以成为IEC标准推荐的风速拟合模型，关键在于它的两个参数各有明确物理意义。形状参数k控制分布曲线的形态：k小于1时曲线呈指数型衰减，k等于2时退化为瑞利分布，k大于3时开始接近正态分布，实际风电场风速的k值通常在1.5到3之间。尺度参数c则直接与平均风速挂钩，c越大，整体风速水平越高，分布曲线的横轴范围越宽。这两个参数一旦估计出来，就可以直接计算平均风速、最可能风速、最大能量风速，以及任意风速区间出现的概率。
在Matlab里做这件事的工具链非常成熟，只需要fitdist或者wblfit函数就能完成参数估计。但如果你只会调用函数而不理解背后的极大似然原理，一旦遇到拟合结果异常或者数据形态很特殊的场景，就很难判断是数据问题还是算法问题。所以下面我会把参数估计的公式和代码逻辑一起讲。
1.2 光伏出力为什么用Beta分布——区间受限数据的天然选择
光伏出力本身受到光照强度、组件温度、转换效率等多重因素影响，建模难度比风电更大。但研究者很早就发现一个关键规律：归一化后的光伏出力数据始终落在0到1的区间内，而且分布形态可以千变万化——有的地区晴天多，出力集中在01附近；有的地区多云天气多，出力偏向0附近的低值；有些场景还会出现两端高中间低的U型分布。
Beta分布恰好是定义在[0,1]区间上的连续分布，它有两个形状参数α和β，通过调整这两个参数，几乎可以覆盖前面提到的所有形态。α大β小时，分布向1偏斜，代表光照条件好的光伏电站；α小β大时，分布向0偏斜，代表光照条件差或者遮挡严重的场景；α和β都很小时，曲线两端高中间低，适合描述多云-晴天交替频繁的地区。
物理上还有一个关键点：Beta分布是贝叶斯统计里的共轭先验，后续如果要做光伏出力的贝叶斯推断，用Beta分布会让后验始终保持在同一个分布族内，计算上非常友好。而在Matlab里，fitdist(x, 'Beta')或betafit都能直接估计参数，但要注意Beta分布要求数据严格在开区间(0,1)内，这就涉及数据预处理的边界处理，我在实操部分会专门给出一套处理方案。
1.3 组合研究到底在解决什么问题
把Weibull和Beta分开做拟合，是很多论文的基础操作，只算“单变量统计建模”。真正的组合研究，是把风速和光照两个随机过程放到同一个分析框架下，回答三类工程问题。
第一类是互补性量化。风能和光能天然存在时间互补，白天光照好但往往风速较小，夜间风速大但无光照。通过组合分布的蒙特卡洛模拟，可以计算风光联合出力的波动率相对单一电源降低多少，这个数值是储能容量配置的重要输入。第二类是置信区间评估。基于拟合出的两个分布，可以模拟出成千上万组风光出力场景，从而给出系统出力在一定置信水平下的上下界，直接服务于可靠性分析。第三类是相关性影响分析。风速和光照强度之间并非完全独立，比如同一个天气系统会同时影响两者，这就需要研究独立假设下的组合结果和存在相关性时的差异。
我在博文里给出的Matlab代码，会把这几个层面的组合分析都覆盖到，这样你拿到的就不只是一个拟合工具，而是一个面向工程评估的完整分析框架。
2. 核心细节解析与实操要点
2.1 Weibull分布参数估计的三种方法
Weibull分布的密度函数是：
f(v; k, c) = (k/c)(v/c)^(k-1) exp(-(v/c)^k)
其中v表示风速，k为形状参数，c为尺度参数。估计这两个参数，工程上常见三种方法。
极大似然估计是精度最高也最推荐的方法。对数似然函数对k和c求偏导并令其为零，化简后得到：
1/k = (∑v_i^k ln v_i)(∑v_i^k)^(-1) - (∑ln v_i)/n
c = ((1/n)∑v_i^k)^(1/k)
注意第一个方程里有未知的k出现在v_i^k中，没法直接解出，需要借助牛顿迭代法或直接让Matlab的mle函数内部完成数值求解。
线性回归法利用Weibull分布的累积函数变换。对F(v) = 1 - exp(-(v/c)^k)取两次对数得到：
ln(-ln(1-F(v))) = k ln v - k ln c
令y = ln(-ln(1-F(v)))，x = ln v，这就是一条斜率为k、截距为-k ln c的直线，用最小二乘拟合即可得到参数。这个方法直观、快速，但精度受经验累积概率F(v)的计算方式影响，在尾部数据上偏差较大。
Matlab里还有一种不用写迭代公式的便捷方法，直接调用wblfit函数或者fitdist：
% 风速数据为列向量 wind_speed pd = fitdist(wind_speed, 'Weibull'); k = pd.Params(1); % 形状参数 c = pd.Params(2); % 尺度参数对比三种方法：极大似然最可靠，适合最终报告和论文使用；线性回归适合快速估算和教学演示；fitdist本质上调用的就是极大似然，只是把迭代细节封装好了。我的建议是正式分析一律以fitdist结果为准，但至少要知道它内部在跑什么，否则遇到不收敛时会无从下手。
2.2 Beta分布参数估计与边界处理
Beta分布的密度函数是：
f(x; α, β) = (Γ(α+β))/(Γ(α)Γ(β)) x^(α-1)(1-x)^(β-1)
其中x是归一化后的光伏出力，取值在0到1之间。估计α和β最简单的做法是矩估计，利用样本均值和样本方差反解：
μ = α/(α+β)
σ² = αβ / ((α+β)²(α+β+1))
反解得到：
α = μ²(1-μ)/σ² - μ
β = α(1-μ)/μ
这个公式看起来简单，但用矩估计在样本量小或数据集中在边界时结果很不稳定，可能出现α、β为负的荒谬结果。更稳健的还是极大似然估计，Matlab中直接用：
pd = fitdist(pv_norm, 'Beta'); alpha = pd.Params(1); beta = pd.Params(2);这里有一个容易踩坑的细节：Beta分布的极大似然要求样本数据严格大于0且严格小于1。但实际光伏出力数据里，夜间出力等于0、中午可能出现1.0（归一化后），这会直接导致估计失败。处理办法是在归一化时留出缓冲区间：把数据压缩到[0.001, 0.999]之间，或者对等于0和等于1的极少数样本做微小扰动。我个人更推荐前一种做法，因为扰动会改变原始数据的概率质量，压缩则更可控。
具体归一化公式可以用：
pv_norm = (pv_data - min(pv_data)) / (max(pv_data) - min(pv_data)); pv_norm = pv_norm * 0.998 + 0.001; % 映射到[0.001, 0.999]这样既保留了原始数据的分布形态，又避开了Beta分布的边界奇点。
2.3 拟合优度怎么量化评估
分布拟合会不会，光靠眼睛看直方图和拟合曲线重叠得好不好是不够的，必须上统计指标。我常用四个指标。
决定系数R²：先对方差直方图的高度做归一化，把每个柱子的概率密度和拟合曲线在柱中点处的理论密度比较，R²越接近1说明拟合越好。均方根误差RMSE：计算所有采样点处经验密度与理论密度的误差平方和均值再开方，RMSE越小越好。KS检验：这是最严格的分布拟合检验，比较经验累积分布函数和理论累积分布函数之间的最大距离，p值小于0.05时说明拟合显著不佳。AIC和BIC：这两个指标用于比较不同分布族拟合同一组数据的优劣，值越小越好，也适合后面做混合分布模型的比较。
Matlab里计算这些指标非常方便，KS检验有现成的kstest，AIC可以从negloglik函数拿负对数似然值后手动计算。关于这些指标的解读，我的经验是：不要只看单一指标。R²高但KS检验不通过的情况经常发生，说明整体形态贴合但局部偏差大，往往数据里有离群风速段。综合使用R²和KS检验，才能对拟合质量有完整把握。
3. 实操过程与核心环节实现
这一部分我会给出一套完整可运行的Matlab代码框架，数据我用了模拟数据来演示，你把它替换成自己的实测风速、光照序列就能直接跑。代码分四个阶段：数据准备、Weibull拟合、Beta拟合、组合分析。
3.1 数据准备与预处理
实际工作中你拿到的风速数据可能来自测风塔，光照数据来自气象站或光伏电站的SCADA系统。两者时间分辨率可能不一致，风电场往往是10分钟平均风速，光伏逆变器可能只有小时级出力记录。我的建议是先统一时间尺度，再对齐时间戳。如果要做四季或月度分析，还需要把数据按季节或月份分组。
数据清洗这块，我列出几个必做步骤：
% 假设已导入风速数据 wind_raw 和光伏出力数据 pv_raw % 1. 剔除异常值：风速为负或超过40m/s以上的记录直接剔除 wind_clean = wind_raw(wind_raw >= 0 & wind_raw <= 40); % 2. 光伏出力非负，且不超过额定装机对应的出力上限 pv_clean = pv_raw(pv_raw >= 0 & pv_raw <= rated_capacity); % 3. 缺失值处理：线性插值补全时间序列中的NaN wind_clean = fillmissing(wind_clean, 'linear');光伏出力归一化用的基准，不同研究里口径不同。有的用理论峰值出力（STC条件下的功率），有的用当月实测最大出力。我推荐用STC理论峰值，因为这样归一化后的数据分布形态更符合Beta分布的建模假设，实测最大值作为基准会把所有数据压缩到更小的区间，导致α、β估计失真。
最后提醒一点：做分布拟合的数据要打乱时间顺序吗？不需要。分布拟合只关心数据的统计特征，不关心时序关系，但后续做组合模拟时需要保留时间相关性，所以建议同时保留原始时序数据和清洗后的数据集。
3.2 Weibull拟合与绘图完整代码
以下代码实现风速直方图叠加Weibull拟合曲线，并输出参数和拟合指标。
% 风速数据 wind_speed = wind_clean; % 单位m/s % 分布拟合 pd_w = fitdist(wind_speed, 'Weibull'); k_w = pd_w.Params(1); c_w = pd_w.Params(2); % 计算拟合优度指标 [~, p_ks_w] = kstest(wind_speed, 'CDF', pd_w); nll_w = negloglik(pd_w); n_w = length(wind_speed); aic_w = 2*nll_w + 2*2; % Weibull有2个参数 bic_w = 2*nll_w + 2*log(n_w); % 绘制对比图 figure('Color', 'w'); histogram(wind_speed, 'Normalization', 'pdf', 'FaceAlpha', 0.3, 'EdgeColor', 'k'); hold on; v = linspace(min(wind_speed), max(wind_speed), 200); plot(v, pdf(pd_w, v), 'r-', 'LineWidth', 2); xlabel('风速 (m/s)'); ylabel('概率密度'); legend('实测频率直方图', 'Weibull拟合曲线'); title(['Weibull拟合 k=', num2str(k_w, '%.3f'), ', c=', num2str(c_w, '%.3f')]); grid on;这里有几个参数计算的说明。k和c的物理意义在风资源评估里可以直接落地：平均风速与c成正比，近似关系是E(v) = c·Γ(1+1/k)，比如k=2.2、c=7.5时，平均风速约等于7.5×Γ(1.455)≈7.5×0.886≈6.6m/s。最大能量风速v_max = c(1+1/k)^(1/k)，这个风速对应风功率密度最大的点，风电机组选型时看这个参数比看平均风速更直接。KS检验的p值如果大于0.05，说明在95%置信水平下认为数据服从Weibull分布的假设成立。
3.3 Beta拟合与绘图完整代码
光伏出力归一化后做Beta拟合，代码如下：
% 光伏出力归一化到(0,1)开区间 pv_norm = (pv_clean - min(pv_clean)) / (max(pv_clean) - min(pv_clean)); pv_norm = pv_norm * 0.998 + 0.001; % Beta分布拟合 pd_p = fitdist(pv_norm, 'Beta'); alpha_p = pd_p.Params(1); beta_p = pd_p.Params(2); % 拟合优度 [~, p_ks_p] = kstest(pv_norm, 'CDF', pd_p); nll_p = negloglik(pd_p); n_p = length(pv_norm); aic_p = 2*nll_p + 2*2; bic_p = 2*nll_p + 2*log(n_p); % 绘图 figure('Color', 'w'); histogram(pv_norm, 'Normalization', 'pdf', 'FaceAlpha', 0.3, 'EdgeColor', 'k'); hold on; x = linspace(0.001, 0.999, 200); plot(x, pdf(pd_p, x), 'b-', 'LineWidth', 2); xlabel('归一化光伏出力'); ylabel('概率密度'); legend('实测频率直方图', 'Beta拟合曲线'); title(['Beta拟合 \alpha=', num2str(alpha_p, '%.3f'), ', \beta=', num2str(beta_p, '%.3f')]); grid on;Beta分布的α、β参数解读有个简单规律：α/(α+β)就是分布均值，α+β越大，分布越集中、方差越小。一个实测案例中我算过某光伏电站的α=3.1、β=2.4，均值0.56，说明这个电站整体出力水平略偏上，而且由于β不算大，低出力时段占比也不低，对应多云天气较多的气候特征。α和β都小于1时分布呈U型，这种情况在天气剧烈交替的地区会出现，拟合结果也合理。
绘制完单分布拟合图之后，建议再加一张Q-Q图来判断尾部拟合情况。Matlab里可以用probplot函数，如果Q-Q图中的散点沿45度角直线分布，说明拟合非常好，特别是在两端的点偏离直线明显时，就要考虑数据是否存在混合分布特征。
3.4 风光组合模型的叠加与互补性计算
组合分析的核心是蒙特卡洛模拟。利用拟合好的Weibull和Beta分布，分别生成大量风速和光照样本，再通过风机和光伏的功率转换模型得到出力序列，合成联合出力后做统计分析。
风机功率转换模型我用分段函数来近似：
function P = wind_power(v, P_rated, v_in, v_out, v_r) % v_in切入风速, v_out切出风速, v_r额定风速 P = zeros(size(v)); idx = v >= v_in & v <= v_r; P(idx) = P_rated * (v(idx) - v_in) / (v_r - v_in); idx = v > v_r & v <= v_out; P(idx) = P_rated; end光伏出力模型简化处理，假设线性关系：
P_pv = pv_norm .* P_pv_rated;组合模拟的完整代码框架：
n_sim = 10000; % 从拟合分布中生成随机样本 sim_wind = random(pd_w, n_sim, 1); sim_pv = random(pd_p, n_sim, 1); % 风功率转换 P_w_sim = wind_power(sim_wind, P_w_rated, 3, 25, 12); % 光伏出力 P_p_sim = sim_pv .* P_p_rated; % 如果考虑风光装机容量比例 % 例如风80MW + 光20MW P_w_rated = 80; P_p_rated = 20; P_total = P_w_sim + P_p_sim; % 互补性量化 std_w_only = std(P_w_sim); std_total = std(P_total); reduction = (std_w_only - std_total) / std_w_only * 100; % 出力不确定性区间 lower_5 = prctile(P_total, 5); upper_95 = prctile(P_total, 95); % 绘制联合出力分布 figure('Color', 'w'); histogram(P_total, 50, 'Normalization', 'pdf', 'FaceAlpha', 0.4); xlabel('风光联合出力 (MW)'); ylabel('概率密度'); title(['互补性削峰效果: 标准差降低', num2str(reduction, '%.1f') '%']); grid on;这个模拟有个重要的前提假设：风速和光照被当作相互独立的随机变量。现实中两者往往负相关，白天风小、晚上风大，独立假设会低估互补性削峰效果。如果手里有同步观测数据，可以用Copula函数建模相关性再模拟，那是一个更大的话题，但先跑通独立模型是第一步。
需要特别说明的是，random函数从pd对象中生成随机数时，Weibull和Beta分布的尾部行为会被忠实保留。也就是说，模拟序列里会出现极小概率的极端大风天和极端低光照日，这恰好是可靠性分析最关心的场景。模拟次数方面，10000次是最低要求，想得到稳定的5%分位数结果，建议至少跑50000次。
4. 常见问题与排查技巧实录
4.1 分布拟合不收敛或参数异常
拟合不收敛的第一典型场景是Beta分布拟合时数据里包含0或1的样本。前面已经提到解决办法，但还有一个容易忽略的问题：归一化时分母用max-min作为基准，如果数据中出现个别极端大的毛刺值，会把所有正常数据压缩到很小范围，导致α、β估计值巨大而且拟合曲线震荡。遇到这种情况先检查归一化前的数据是否有离群点，比如光伏出力偶尔出现超过额定功率的异常记录。
Weibull拟合出现k值异常偏大或偏小也要警惕。k大于5时说明风速数据非常集中，这通常不是全年数据而可能是某一小段时间的数据；k小于1时说明小风速占绝对主导，要检查是否包含了大量静风记录。如果数据确实包含大量零风速，考虑用零膨胀Weibull模型，而不是单纯Weibull拟合。
4.2 KS检验始终不通过怎么办
实测数据经常出现KS检验p值小于0.05的情况，很多人会非常焦虑，觉得模型失效。我的经验是分三步排查。第一步，看Q-Q图确认偏差发生在哪一段。如果偏差集中在右尾，说明大风速段有混合分布特征，考虑用双峰混合Weibull模型。第二步，增加样本量再试。KS检验对样本量非常敏感，数据量上万时细微偏差都会被放大，这时更应关注RMSE和R²，或者把显著性水平放宽到0.01。第三步，按时间分段重新拟合。全年数据往往包含了不同季节的风速特征，分季节拟合后各段的KS检验往往能通过。
Beta分布拟合后发现α和β都小于1，KS检验如果还不通过，就要考虑是不是Beta分布根本不适用。有两种情况：一是数据里存在一个强峰值接近0.5且围绕峰集中，这时峰度比Beta分布更尖锐，换用混合Beta分布效果更好；二是数据有很多极端接近0或1的点，Beta分布虽然形态上能拟合，但局部概率密度与实际偏差大，需要先做滤波或重采样。
4.3 组合模拟结果明显偏离实际
这个问题比较隐蔽。有时候单分布拟合各项指标都很好，合并模拟出来的联合出力期望值却和实测同期联合出力对不上。我排查过几次，最常见的原因是装机容量比例设置错误。风光容量配比必须和模拟研究的情境一致，有的研究做成50:50，有的则反映实际场站比例，两者结果差异巨大。
其次是功率转换模型的参数问题。风机额定风速、切入切出风速如果不按实际机型设置，模拟的功率曲线会在高风速段大幅失真。我建议实测的风速-功率散点图拿出来先拟合一遍，确认分段模型的参数正确后再做组合模拟。
最后注意时间尺度的匹配。如果你的Weibull分布是从10分钟平均风速数据拟合出来的，而Beta分布用的是小时级光伏出力数据，两者混用会引入尺度不一致的系统偏差。统一时间尺度是组合研究的基本前提。
4.4 代码层面的几个细节坑
画图时中文乱码是Matlab老生常谈的问题，用set(gcf, 'Color', 'w')搭配保存时指定字体可以解决大部分显示问题：
set(0, 'DefaultAxesFontName', 'SimHei'); set(0, 'DefaultTextFontName', 'SimHei');fitdist与mle函数对参数顺序的定义不一致也是一个容易踩的坑。fitdist返回的Weibull对象Params顺序是[k, c]，而mle(wind_speed, 'pdf', @(x,k,c) wblpdf(x,k,c), 'start', [1,1])的输入顺序必须和wblpdf一致。Beta分布同样，betafit返回[a, b]，如果直接用mle需确认参数顺序。
我在实际跑数万次模拟时还遇到过random函数生成大量随机数占用内存过高的问题。解决方法是一次生成一个大的随机矩阵，不要循环调用。比如random(pd_w, 50000, 1)一次性生成，后续所有分析都基于这个矩阵切片完成。
5. 一些实操体会
用Matlab把Weibull和Beta分布的组合研究完整跑通之后，再看风光出力数据的感觉会完全不同。你不再只盯着均值、最大值这些描述统计量，而是能回答更复杂的问题：风速在切入风速附近的概率有多大，归一化光伏出力在0.3以下的时段占多少比例，风光互补后出力低于某个阈值的天数期望是多少。这些答案直接决定了储能配置的容量和调度策略的保守程度。
我个人在做这类分析时，最看重的是拟合参数的可解释性。不管数据换了几批、地区换了几个，k、c、α、β的变化总能对应上气候和场址的物理特征。比如内陆和沿海的k值差异、阴雨地区和干旱地区的α、β差异，这些规律反过来还能用于数据缺失地区的参数估算，这算是分布拟合研究最有工程价值的部分。
最后分享一个小技巧：把整套拟合流程封装成函数后，配合generatePDF批量导出不同月份、不同季节的分析报告，能大幅提升做风光资源评估的效率。后续如果数据量大了或者要分析多个场站，可以把这套代码并行化跑，每个场站的拟合互不干扰，输出结果统一汇总到一个结构体里，再做横向对比分析。这套流程我已经在多个实际项目中验证过稳定性，希望能帮你少走一些弯路。
