1. 为什么风功率用Weibull、光伏用Beta——分布选型的底层逻辑做新能源出力建模的朋友对Weibull分布和Beta分布这两个名字应该都不陌生。风速基本都用两参数Weibull光伏出力归一化之后习惯用Beta但真正要把它们组合起来做风光互补分析中间其实隔着一整套完整的参数估计、拟合检验和场景生成流程。很多人手头有数据、有Matlab但卡就卡在不知道每一步怎么落地尤其是“两个分布组合”这个需求——不是简单把两张图拼在一起而是要围绕它们建立一套能复现、能扩展的建模框架。先解决一个最基础的问题为什么偏偏是这两个分布这不是拍脑袋选的背后有物理机制和数据形态两方面的原因。1.1 风速的物理机制与Weibull分布的适配原因风速本质上是一个非负随机变量受大气环流、地形、地表摩擦等多因素叠加影响实测风速序列的概率密度通常呈现明显的右偏特征——低速时段占比高偶尔出现强风拖出一条长尾。两参数Weibull分布的概率密度函数长这样f(v) (k/A) · (v/A)^(k-1) · exp(-(v/A)^k), v ≥ 0其中A是尺度参数量纲与风速一致决定分布整体的“位置”k是形状参数决定分布曲线的形态。k值越大曲线越集中k越接近1.52.5时密度曲线和实测风速直方图贴合得非常好。这个分布能用一个数学式子同时描述“小风居多、大风偶发”的统计规律又只带两个参数估计起来不麻烦所以在风资源评估、风电出力模拟、可靠性计算中被广泛采用。我自己最常用也最推荐的拟合方式是最大似然估计MLE因为MLE的估计值在大样本下具备一致性和渐近有效性。对一组独立风速样本 v1, v2, …, vn写出对数似然函数后对A和k求偏导置零可以推导出形状参数k满足的隐式方程1/k - (Σ vᵢ^k · ln vᵢ) / Σ vᵢ^k (1/n) · Σ ln vᵢ 0这个方程没有解析解需要用数值迭代牛顿法或者Matlab自带的wblfit函数内部处理。尺度参数A则可以直接显式计算A ( (1/n) · Σ vᵢ^k )^(1/k)1.2 光伏出力归一到[0,1]后为什么妥帖地跟着Beta走光伏出力或者说归一化辐照度是一个取值严格落在[0,1]区间的随机变量。Beta分布恰好定义在[0,1]上概率密度为f(x) x^(α-1) · (1-x)^(β-1) / B(α, β), 0 ≤ x ≤ 1, α 0, β 0其中B(α, β)是Beta函数起到归一化作用。α和β两个形状参数一个控制左侧的“抬升”程度一个控制右侧的“衰减”程度。晴天时光伏出力会集中在较高区间数据分布偏向右端此时α显著大于β阴雨天出力普遍偏低分布偏向左端β会更大。这种灵活的形态变化让Beta分布在描述不同气候区、不同季节的光伏出力时都很从容。使用Beta分布前必须先做归一化——把实测出力除以装机容量变成标幺值。这一步看起来简单但很多人在这里栽跟头如果直接用原始兆瓦数值去拟合矩估计公式会直接失效因为Beta分布不支持超出[0,1]区间的数据。1.3 一个容易被混淆的边界风速用Weibull风电功率不要直接套Weibull这里必须多说一句它关系到整个研究的严谨性。多数文献说的是风速服从Weibull分布而不是风电功率服从Weibull分布。风电功率是由风速经过风机功率曲线映射得到的而功率曲线带着切入风速、额定风速、切出风速三段非线性特征映射之后功率分布的形状和风速分布已经完全不同了。正确的做法是先用Weibull拟合风速得到参数后用蒙特卡洛抽样生成风速场景再经过功率曲线换算成风电功率场景。这样物理过程就是自洽的。如果你手头只有现成的风电功率数据也别直接拿Weibull去硬套可以考虑用核密度估计或者分段分布去描述硬套的结果往往是拟合检验被拒得很难看。2. 参数估计方法拟合出靠谱参数的关键细节分布选型定下来之后核心工作就转到参数估计上。你当然可以直接调用Matlab的wblfit和betafit但如果不理解这些函数背后的逻辑遇到拟合结果异常时只能干瞪眼。这一节把两套估计方法都讲透包括初值怎么选、数据怎么预处理、结果怎么验证。2.1 Weibull双参数的最大似然估计与初值选择用Matlab自带的wblfit最省事一行代码就能得到[尺度参数A, 形状参数k]的估计值% 假设v是一列风速数据m/s单位要统一 [param_hat, param_ci] wblfit(v); A_hat param_hat(1); % 尺度参数 k_hat param_hat(2); % 形状参数但wblfit内部用的迭代算法需要初值它会自动用矩估计思想选初始点。绝大多数情况下这没问题可一旦风速序列里包含大量静风时段数值接近0或者数据量特别少迭代就可能不收敛或者收敛到不合理的参数上。所以我习惯自己实现一遍MLE主要图两点一是能完全控制初值二是能借此充分理解数据形态。自己实现的思路是用fminsearch最小化负对数似然函数。负对数似然写成NLL(A, k) -n·ln(k) n·k·ln(A) - (k-1)·Σ ln(vᵢ) Σ (vᵢ/A)^k直接对A和k两维优化即可。初值的选择建议用矩估计的经验公式先算出风速样本的均值和标准差A的初值取均值附近对大多数风况尺度参数和平均风速数值接近k的初值取1.52.5之间的值。这样优化过程又快又稳。2.2 Beta分布的矩估计与极大似然估计以及数据预处理Beta分布的参数估计有两种常用路线矩估计和最大似然估计。我建议先用矩估计给初值再用最大似然精修。矩估计的核心是用样本均值 μ 和样本方差 σ² 反推参数α μ · ( μ·(1-μ)/σ² - 1 ) β (1-μ) · ( μ·(1-μ)/σ² - 1 )注意中间那项 μ·(1-μ)/σ² 必须大于1。如果算出来小于等于1说明样本方差太大数据形态不适合Beta分布这时候要回头检查数据里是不是混合了夜间零出力时段和白天非零出力时段。数据预处理我总结了一套固定流程把原始功率序列除以装机容量得到[0,1]区间标幺值去掉夜间必然为零的时段一般按辐照度阈值或者太阳高度角判断对剩余数据再做一次截断把极小值比如小于0.001和极接近1的值做边界处理避免Beta密度函数在端点处奇异用处理后的数据做矩估计求初值再用fmincon或betafit做最大似然精修。% 光伏出力归一化与Beta拟合p_raw是原始功率 p_pu p_raw / Pn; % 标幺化 p_pu p_pu(p_pu 0.001); % 去掉零值和接近零的时段 p_pu min(max(p_pu, 0.001), 0.999); % 避免端点奇异 [alpha_hat, beta_hat] betafit(p_pu);2.3 拟合质量怎么判断KS检验、Q-Q图、误差指标参数算出来不代表工作完成了还必须验证拟合质量。最常用的定量方法是Kolmogorov-Smirnov检验KS检验它比较经验累积分布函数和理论累积分布函数的最大差距。Matlab里这样写% Weibull拟合的KS检验 [h1, p1] kstest(v, CDF, (x) wblcdf(x, A_hat, k_hat)); % Beta拟合的KS检验 [h2, p2] kstest(p_pu, CDF, (x) betacdf(x, alpha_hat, beta_hat));h0表示在默认显著性水平0.05下不拒绝原假设说明拟合效果可以接受h1则说明拟合被拒绝。但我要提醒一句当样本量特别大比如一年8760个小时的数据KS检验会非常敏感一个很小的偏差都会导致拒绝。这时候别急着改分布先画Q-Q图看一眼偏差到底在哪个区间。如果偏差集中在尾部比如极端大风时段对后续场景抽样影响有限可以接受如果偏差发生在高概率密度的核心区间就得认真考虑是不是数据预处理出了问题。除了KS检验我还会计算均方根误差RMSE和决定系数 R²把经验直方图对应的概率密度与拟合分布的理论密度做对比定量衡量拟合精度。这几个指标放在一起才算一道完整的检验工序。项目Weibull拟合Beta拟合核心参数尺度A、形状k形状α、形状β主要估计方法最大似然MLE矩估计定初值 MLE精修数据要求非负连续风速归一化到[0,1]的光伏出力初值建议A取均值附近k取1.52.5矩估计公式直接计算常用检验KS检验、Q-Q图、RMSEKS检验、Q-Q图、RMSE3. Matlab代码实现从零跑通风光分布拟合与组合场景生成下面这套代码我按可复现的标准来写。你只需要把你的实测数据替换掉对应变量名整个流程可以直接跑通。我用模拟数据做演示这样你可以在没有实测数据的情况下先把流程验证一遍。3.1 风速数据的Weibull拟合完整代码与绘图%% 第1步生成模拟风速数据一年8760小时 rng(2026); v_true_A 6.8; % 真实尺度参数约等于平均风速 v_true_k 2.1; % 真实形状参数 v wblrnd(v_true_A, v_true_k, 8760, 1); v min(v, 30); % 模拟风机切出保护截断到30m/s %% 第2步Weibull参数拟合 [param_hat, param_ci] wblfit(v); A_hat param_hat(1); k_hat param_hat(2); fprintf(拟合尺度参数 A %.4f形状参数 k %.4f\n, A_hat, k_hat); %% 第3步绘制概率密度拟合对比图 figure(Color, w, Position, [100 100 900 400]); subplot(1,2,1); histogram(v, Normalization, pdf, NumBins, 30, ... FaceColor, [0.7 0.7 0.7], EdgeColor, none); hold on; x_axis 0:0.1:30; y_fit wblpdf(x_axis, A_hat, k_hat); plot(x_axis, y_fit, r-, LineWidth, 2); xlabel(风速 / (m/s)); ylabel(概率密度); title(风速直方图与Weibull拟合曲线); legend(实测直方图, Weibull拟合, Location, northeast); grid on; %% 第4步KS检验 [h1, p1] kstest(v, CDF, (x) wblcdf(x, A_hat, k_hat)); fprintf(Weibull KS检验: h %d, p %.4f\n, h1, p1);跑完这段代码你会看到拟合出的A和k与真实模拟参数非常接近。用真实风速数据时直方图和拟合曲线重叠度高的区域一目了然。有一点要特别注意wblrnd生成的模拟数据是分布已知的代码里用于验证流程闭环没问题但换成实测数据后如果直方图和拟合曲线在风速接近0的区间严重偏差就要考虑数据采集设备的启动风速限制——很多测风设备在低于0.5m/s时测不准会造成该区间数据异常。3.2 光伏出力的Beta拟合完整代码与绘图%% 第1步生成模拟光伏出力数据归一化标幺值 % 用Beta分布模拟白天时段出力叠加96点日曲线 time_steps 96; % 一天96个15分钟点 days 365; n_total time_steps * days; % 构造太阳位置系数简化为正弦形状白天才能有出力 solar_env repmat(max(sin(linspace(0, pi, time_steps)), 0), days, 1); alpha_true 5.0; % Beta真实参数 beta_true 2.5; base_pv betarnd(alpha_true, beta_true, n_total, 1); p_pv base_pv .* solar_env; % 夜间时段乘0 p_pv p_pv(p_pv 0.01); % 去掉夜间零值 %% 第2步Beta参数拟合矩估计初值 betafit精修 mu mean(p_pv); var_pv var(p_pv); % 矩估计初值 alpha_mom mu * (mu * (1 - mu) / var_pv - 1); beta_mom (1 - mu) * (mu * (1 - mu) / var_pv - 1); % 最大似然精修如果不放心矩估计结果可直接用betafit [alpha_hat, beta_hat] betafit(p_pv); fprintf(矩估计初值: alpha %.4f, beta %.4f\n, alpha_mom, beta_mom); fprintf(MLE拟合结果: alpha %.4f, beta %.4f\n, alpha_hat, beta_hat); %% 第3步绘制拟合对比图 subplot(1,2,2); histogram(p_pv, Normalization, pdf, NumBins, 30, ... FaceColor, [0.9 0.8 0.6], EdgeColor, none); hold on; x_beta 0.001:0.001:0.999; y_beta_fit betapdf(x_beta, alpha_hat, beta_hat); plot(x_beta, y_beta_fit, b-, LineWidth, 2); xlabel(归一化光伏出力); ylabel(概率密度); title(光伏出力直方图与Beta拟合曲线); legend(实测直方图, Beta拟合, Location, northeast); grid on; %% 第4步KS检验 [h2, p2] kstest(p_pv, CDF, (x) betacdf(x, alpha_hat, beta_hat)); fprintf(Beta KS检验: h %d, p %.4f\n, h2, p2);这段代码里有一个关键的仿真设计光伏出力不是直接抽一个Beta随机数就完事而是先采样Beta随机数再乘一个“太阳位置”系数。这样能模拟夜间零出力现象使数据更接近真实情况。处理真实数据时同样要遵循这个逻辑建模的目标是白天的非零出力部分零出力时段要单独统计概率不要和连续分布混在一起。3.3 组合研究从分布抽样到风光联合出力场景拟合出两个分布参数后组合研究的真正重头戏才开始。所谓“组合”就是把风速的Weibull分布和光伏出力的Beta分布接进同一个场景生成框架然后对生成的风光联合出力序列做各种分析。先看核心抽样代码%% 组合场景生成参数 N_scenario 20000; % 场景数量越大分布越稳定 P_w_rated 1500; % 风机额定功率 kW P_pv_rated 1000; % 光伏装机容量 kW v_ci 3; v_r 12; v_co 25; % 切入、额定、切出风速 %% 第1步从Weibull分布抽样风速 v_sample wblrnd(A_hat, k_hat, N_scenario, 1); %% 第2步简化功率曲线风速转功率 % 线性段切入风速到额定风速之间线性上升 % 额定段额定风速到切出风速之间保持额定 % 零功率段低于切入风速或高于切出风速 P_w zeros(N_scenario, 1); idx_linear v_sample v_ci v_sample v_r; idx_rated v_sample v_r v_sample v_co; P_w(idx_linear) P_w_rated .* (v_sample(idx_linear) - v_ci) ./ (v_r - v_ci); P_w(idx_rated) P_w_rated; % 归一化让风功率标幺值在[0,1]区间与光伏一致 P_w_norm P_w / P_w_rated; %% 第3步从Beta分布抽样光伏出力 % 先生成白天非零出力再按日照概率叠加夜间零时段 solar_prob 0.55; % 白天时段占比可按季节调整 pv_nonzero betarnd(alpha_hat, beta_hat, N_scenario, 1); pv_mask rand(N_scenario, 1) solar_prob; % 判断是否白天 P_pv_norm zeros(N_scenario, 1); P_pv_norm(pv_mask) pv_nonzero(pv_mask); %% 第4步组合出力风电光伏 % 这里假设风电场、光伏电站同时接入同一母线忽略线路损耗 P_total_norm P_w_norm P_pv_norm; P_total_norm min(P_total_norm, 1.0); % 限幅模拟逆变器/变流器限功率 %% 第5步组合结果分析 figure(Color, w, Position, [120 120 1200 400]); subplot(1,3,1); histogram(P_w_norm, Normalization, pdf, NumBins, 40); xlim([0 1]); xlabel(风电标幺出力); ylabel(概率密度); title(风电出力分布); subplot(1,3,2); histogram(P_pv_norm, Normalization, pdf, NumBins, 40); xlim([0 1]); xlabel(光伏标幺出力); ylabel(概率密度); title(光伏出力分布); subplot(1,3,3); histogram(P_total_norm, Normalization, pdf, NumBins, 40); xlim([0 1]); xlabel(组合标幺出力); ylabel(概率密度); title(风光组合出力分布); % 输出统计指标 fprintf(风电出力均值: %.4f方差: %.4f\n, mean(P_w_norm), var(P_w_norm)); fprintf(光伏出力均值: %.4f方差: %.4f\n, mean(P_pv_norm), var(P_pv_norm)); fprintf(组合出力均值: %.4f方差: %.4f\n, mean(P_total_norm), var(P_total_norm));这段代码生成的组合出力序列后续能做很多延展分析把P_total_norm和负荷曲线做差可以算电量不足期望值LOLE把组合出力曲线按不同光伏装机配比重复抽样可以找到方差最小的最优配比把结果输入储能容量优化模型可以得到不同置信水平下的储能配置建议。3.4 结果解读与延伸思路跑完上面代码你会看到三个明显的现象。第一个现象是风电出力分布呈现“两头大中间小”的形态大量场景里风速不在有效发电区间内出力为0一部分场景出力达到额定值上限真正连续分布的中间段反而不多。这就是功率曲线映射造成的效果直接用Weibull拟合功率分布根本复现不了这个特征。第二个现象是光伏出力分布中零出力场景占比基本被solar_prob参数控制住了剩下的非零出力服从Beta分布。这说明光伏出力的随机性应该拆成两个部分来理解日升日落的确定性过程决定零出力概率云层遮挡和大气透明度决定非零阶段的连续波动。第三个现象是组合出力方差。大多数情况下风电和光伏的出力峰值时段是错开的——白天光伏强、夜间风电往往更稳组合后的总出力方差会明显小于两者单独出力方差之和。这个“互补平滑效应”是风光组合研究的核心价值所在。你可以在此基础上做一层简单延伸遍历不同的光伏装机比例系数比如从0到2.5按0.1步长变化重新计算P_total_norm P_w_norm ratio * P_pv_norm记录每种配比下组合出力方差、低于某个阈值比如0.05装机容量的概率。画出一条曲线后最优风光配比基本一目了然。这个分析用到的新代码不多但产出非常有说服力。4. 组合研究到底能解决什么问题一个小的容量配置示例单纯把两个分布拟合出来学术上可以发一篇笔记但工程价值有限。组合研究的价值在于它为后续的系统分析提供了大量的随机场景。这里用一个具体的容量配置示例说明这套方法怎么落地。4.1 为什么要做风光分布的组合研究风光电站规划时最核心的一个问题是在保证供电可靠性的前提下风电和光伏的装机容量怎么配最经济这类问题没法用“平均出力”算因为风光出力是强随机变量不同时段的相关性和互补性直接决定系统可信容量。先从两个分布独立抽样再组合成大量场景本质上是在用蒙特卡洛方法做“数字试验”——每一条抽样结果都是一次“可能发生的天气状态”跑几千上万次就能系统性地评估在不同装机配比下系统失负荷的概率和电量。我的经验是这个框架特别适合以下三类应用场景风电场和光伏电站共建基地的容量配比规划比如在同一升压站下风电装多少、光伏装多少才能让送出线路利用率最高含高比例新能源的配电网可靠性评估用组合出力序列和负荷曲线做时序仿真统计年缺电时间储能容量配置的前期估算组合出力的波动幅度直接决定储能是偏向能量型还是功率型需要多大功率和容量。4.2 基于场景抽样的容量可信度与不足电量分析容量可信度衡量的是风光电站能替代多少常规机组容量。计算思路很简单把组合出力序列按升序排列看它能稳住到多高的出力水平。比如取P90值意思是90%的场景下组合出力不低于这个值那么从可靠性角度看这套风光系统至少可以“可信地”提供这么大体量的容量。% 基于前述P_total_norm场景结果 p90 prctile(P_total_norm, 10); % 注意低百分位数对应高可靠性 p50 prctile(P_total_norm, 50); fprintf(组合出力P90 %.4f, P50 %.4f\n, p90, p50); % 用电量不足期望值分析示例假设负荷恒定在组合出力均值的80% load_level 0.8 * mean(P_total_norm); shortage max(0, load_level - P_total_norm); EENS mean(shortage); % 电量不足期望值 fprintf(电量不足期望值 EENS %.4f标幺值\n, EENS);EENS这个指标非常有意思它把场景抽样和系统可靠性联系起来了。EENS单位是“标幺电量”乘以风光总装机容量和年小时数就得到实际电量缺口。做风光配比优化时把不同配比下的EENS画成曲线再叠加成本信息最优解就浮出水面了。4.3 从独立抽样到相关性建模的进阶方向前面的组合场景假设风电和光伏出力相互独立这个假设在多数情况下可以接受但严格来说同一个地区风电和光伏受同一天气系统影响两者存在一定的负相关或正相关性——多云天气往往伴随大风晴空天气光伏出力高但风速可能偏弱。如果研究精读要求更高可以考虑引入Copula模型把风速和辐照度的尾部相关性刻画出来。从技术栈演进的角度看WeibullBeta组合属于“边缘分布明确、相关性独立”的第一阶段模型第二阶段是使用经验Copula或Gaussian Copula来捕捉相关性第三阶段是引入时间序列模型如ARIMA-GARCH加上分布假设刻画逐时出力动态过程。但所有进阶模型的边缘分布仍然绕不开本文讲的Weibull和Beta拟合这一层。参数拟合的可靠性决定上层模型的天花板。5. 实操中踩过的坑常见问题与排查心得这部分是我最想分享的内容。理论方法和代码框架大家都能搜到但真正跑数据时那些让人抓狂的细节通常只有自己踩过坑才能体会。5.1 数据零值太多导致拟合失败最典型的坑来自数据里的零值。风速序列里的静风时段、光伏序列里的夜间时段、还有停机检修时段都会产生大量精确为0的数据。直接把这样的数据丢给wblfit或者betafit问题立刻暴露。对Weibull拟合零值本身不是大问题但如果静风时段占比太高拟合出来的形状参数k会被“压小”导致分布曲线把大量概率质量推向接近0的区间中间段的密度反而被低估。处理办法是单独统计静风占比把它当作一个离散概率点在正风速子集上拟合连续的Weibull部分。对Beta拟合0和1端点的处理更加微妙。Beta分布在端点处密度可能趋向无穷数据里一旦混入精确的0或1最大似然估计就会向端点“塌缩”。我的做法是先物理判断0值来源夜间时段直接剔除白天时段如果出现0大概率是辐照度不足或者设备死区同样剔除。然后把剩余数据压缩到[0.001, 0.999]给端点留出安全余量。5.2 拟合参数与实际物理量对不上有一种情况特别迷惑人拟合结果、统计检验都通过了但参数数值怎么看怎么不对劲。比如风速均值10m/s的地方Weibull形状参数k拟合出来却是0.8尺度参数A高达15。这时候别急着信结果先画直方图看看数据形态。形状参数k小于1意味着密度函数在风速接近0处趋向无穷大部分概率堆在低速区。风速均值10m/s的数据不可能拟合出这种参数除非数据里有大量记录错误的超低速时段或者时间断档。这是数据质量问题的信号不是分布模型的问题。遇到这种情形先回去检查原始数据的时间戳连续性、传感器量程、以及是否有单位混用问题m/s和knots混在一起是常见事故。5.3 工具箱版本与随机数差异的坑Matlab不同版本之间随机数生成器默认算法有差异。以前用老版本写的数据模拟代码在新版本里跑出的随机序列可能完全不同。直接影响就是同样一份代码两次运行出来的拟合参数有细小差异KS检验的p值也可能在不同版本间波动。这不是bug是随机数种子和算法变更造成的正常现象。解决办法有两个一是写代码时统一用rng(固定种子)显式设置随机数种子保证实验可复现二是报告结果时同时记录Matlab版本号方便后续追溯。如果做学术研究这个细节很影响结果可复现性审稿人问到的时候你能拿出证据。5.4 不要过度拟合分布最后一个建议可能听起来有点反直觉。拟合优度指标越好不代表模型越实用。我在项目里见过有人为了提高R²、降低RMSE把简单两参数Weibull换成三参数Weibull、混合Weibull甚至非参数核密度模型。拟合精度确实提升了但代价是参数数量暴增、物理可解释性下降场景生成时出现大量极端值后续系统分析反而被带偏。做工程的时候我给自己定了一条规矩分布模型的复杂度不超过“能准确描述核心区间的概率分布特征并保留合理的尾部行为”这个限度。两参数Weibull和两参数Beta在这个框架下通常已经够用。先跑通这条链路再考虑扩展。坦白讲风光出力建模这条路难点从来不在某个单独环节而在整条链路的串联。每一个参数都需要问到“它怎么来的、怎么验证、怎么用起来”。把Weibull和Beta这一关打通后面的Copula、时序模型、优化算法才有立足点。希望这份实操笔记能帮你少走一些弯路把时间花在真正有价值的分析上。