搞新能源建模的朋友估计都绕不开这两件事怎么描述风力的随机性又怎么描述光伏的随机性。风电功率和光伏功率的出力曲线都不是线性稳定的但如果我们把风速、光照强度这类原始气象要素用合适的概率分布刻画清楚后面做容量配置、可靠性分析、储能优化全都变得有据可依。风电场出力建模业内最常用的就是两参数Weibull分布光伏出力建模最经典的则是Beta分布。把这两个分布拿到一起做组合研究再用Matlab实现整套拟合、绘图与统计分析流程正是这篇博文要解决的完整问题。这篇内容适合谁一个是正在做新能源并网、微电网规划类课题的研究生另一个是从事风光资源评估的前期设计人员。你不需要把数理统计从头学一遍只要跟着我把分布参数估计的原理弄明白、把代码框架拿走去改数据基本就能直接上手出结果。我会把Weibull和Beta两个分布从物理背景、参数推导到Matlab代码实现全部拆开讲透包括拟合后怎么评估好坏、组合研究到底在做什么以及我实测踩过的那些坑。1. 内容整体设计与思路拆解1.1 风速为什么选Weibull而非正态——先把物理直觉讲清楚很多刚接触风资源分析的人会先入为主地认为风速数据取个平均画个直方图长得像钟形就可以用正态分布近似。但实际的风速观测数据几乎全是右偏的也就是小风速出现的频率很高大风速出现的频率很低拖一条长长的右尾。正态分布是对称的它既不能刻画“微风日数多”这种常见形态又无法处理风速永远大于等于0这个天然边界。瑞利分布虽然也能描述风速但它本质上是Weibull分布形状参数k2的特例灵活性太差。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)/nc ((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检验有现成的kstestAIC可以从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·Γ(11/k)比如k2.2、c7.5时平均风速约等于7.5×Γ(1.455)≈7.5×0.886≈6.6m/s。最大能量风速v_max c(11/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分布拟合后发现α和β都小于1KS检验如果还不通过就要考虑是不是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批量导出不同月份、不同季节的分析报告能大幅提升做风光资源评估的效率。后续如果数据量大了或者要分析多个场站可以把这套代码并行化跑每个场站的拟合互不干扰输出结果统一汇总到一个结构体里再做横向对比分析。这套流程我已经在多个实际项目中验证过稳定性希望能帮你少走一些弯路。