1. 风电与光电联合出力的概率建模思路做新能源发电出力分析的人几乎都绕不开两个分布描述风电出力的Weibull分布以及描述光电出力的Beta分布。把这两个分布组合起来研究本质上是在解决一个非常现实的问题——风、光电源同时接入电网时它们的联合出力特性到底是什么样的。单独看风电或者单独看光电各自都有成熟的概率模型但实际电网调度、储能配置、可靠性评估关心的往往是“风光”这个整体系统的出力波动范围。所以Weibull与Beta分布的组合研究核心不是为了炫技而是为了给后续的随机潮流计算、容量可信度评估、备用容量优化打下一个可量化的概率基础。Matlab在这个场景里几乎是标配工具。它处理概率分布的函数生态非常完善从wblpdf、wblfit到betapdf、betafit一行代码就能完成拟合和绘图。更关键的是Matlab的随机数生成机制便于做蒙特卡洛模拟这对于验证组合模型的输出特性非常友好。我见过不少做新能源方向的研究生前期用Python做数据分析到了要出概率模型、做系统仿真的时候还是会被导师要求切回Matlab——因为后续要接Simulink或者Matpower的话Matlab的接口最少折腾。这篇文章不是要推一套复杂的理论推导而是把我实际跑通过的一套“风电Weibull分布光电Beta分布组合建模”的Matlab实现方案完整拆开。从参数拟合到联合采样从出力时序生成到统计验证每一步都给出可复现的代码和参数选择的依据让还没入门的同学能直接照着跑也让已经做了一段时间的人能对照检查自己的模型有没有踩坑。2. 为什么风电选Weibull、光电选Beta分布特性与机理匹配2.1 风速的偏态特性与Weibull分布的适配逻辑先聊风电。风速的概率分布不是对称的绝大多数地区的中等风速出现概率最高但极端大风比如台风边缘、强对流天气虽然少却会造成出力骤降或切机所以风速分布的右尾相对较长呈明显的正偏态。Weibull分布恰好能描述这种形状它的概率密度函数为f(v) (k / c) * (v / c)^(k-1) * exp(-(v/c)^k) v ≥ 0其中k是形状参数c是尺度参数。k的取值直接决定了分布的形状——k1时退化为指数分布适合描述极端风况k在2到3之间时曲线接近实际风速的典型分布k3时分布逐渐对称偏态减弱。c参数则跟平均风速强相关通常可以通过平均风速和标准差来估计。实际做风电出力建模时不能直接用风速分布代替出力分布因为风机存在切入风速、额定风速、切出风速三个关键点。风速在切入风速以下时出力为零在额定风速和切出风速之间时出力恒定等于额定功率只有中间段才近似按三次方关系变化。所以标准的处理流程是先用实测风速序列拟合Weibull分布参数然后把风速转换成功率转换公式一般为P(u) 0u u_cin 或 u u_coff P(u) P_r * (u - u_cin) / (u_r - u_cin)u_cin ≤ u ≤ u_r P(u) P_ru_r u ≤ u_coff这个转换导致风电出力概率分布在0和额定功率处出现明显的“尖峰”也就是概率堆积。这个特性在后续组合模型中必须保留否则联合出力的结果会失真。2.2 光照辐照度的有界特性与Beta分布的适配逻辑光电的物理机理跟风电差别很大。光伏出力主要取决于水平面辐照度而辐照度天然是“有界”的——不可能为负也不可能超过大气上界的太阳常数约1367 W/m²实际地面辐照度更有上限。因此用无界的正态分布或对数正态分布去拟合辐照度理论上就不严谨。Beta分布定义在[0,1]区间正好跟“归一化辐照度”或者“光伏出力标幺值”的取值域吻合。Beta分布的概率密度为f(x) (1 / B(a, b)) * x^(a-1) * (1-x)^(b-1) 0 ≤ x ≤ 1其中a和b是两个形状参数。a大于b时分布偏向高值区适合晴天比例大的地区a小于b时分布偏向低值区适合多云或阴雨较多的地区。a和b可以通过辐照度样本的均值和方差反推a μ * (μ * (1-μ) / σ² - 1) b (1-μ) * (μ * (1-μ) / σ² - 1)这里μ是归一化辐照度的均值σ²是方差。实际使用中要注意如果样本方差特别小会导致a和b非常大Beta分布变得非常尖这时候要检查是不是数据时间段选得太短比如只取了正午前后两小时的数据那方差肯定偏小拟合出来的分布不能代表全天特性。2.3 组合研究的目标不是求解析解而是生成联合出力样本很多人第一次看到“组合研究”会误以为是要推导风电和光电联合概率密度的解析表达式。直接写出风、光联合概率密度确实可以——如果假设风速和辐照度独立联合密度就是两个边缘密度的乘积。但实际上风速和辐照度并不完全独立它们都受气象条件影响比如阴雨天风速往往偏大晴天风速相对平稳这就产生了相关性。所以更工程化的做法是利用Copula函数把两个边缘分布耦合起来或者干脆用蒙特卡洛生成大量样本直接统计联合出力的分布特性。Matlab实现蒙特卡洛组合建模的路径非常清晰先用wblrnd生成风速样本再用betarnd生成归一化辐照度样本然后同时考虑风速-辐照度的相关系数通过排序法rank correlation对样本进行重新排列让生成的样本对满足目标相关性。这个过程虽然不像解析推导那样优雅但胜在灵活能应对复杂的实际情况。3. Matlab代码实现的完整框架与参数选择3.1 数据准备与预处理我用的示例数据是一组某风电场实测风速序列单位m/s时间分辨率15分钟和同站点的光伏电站辐照度序列单位W/m²时间跨度为一年。数据本身带有典型的新能源实测特征白天辐照度随太阳高度角变化夜间为零风速则呈现明显的季节性和日变化。拿到数据后第一步不是急着拟合分布而是清洗。需要剔除三类异常数据风速为负或者超过50 m/s的明显坏点辐照度夜间非零或超过太阳常数的异常值风、光数据时间戳不对齐的缺失段。清洗完之后再做时间维度的对齐。因为风电出力波动相对平缓而光伏出力突变快如果原始采样频率不一致建议统一插值到相同的时间分辨率我采用的是线性插值。3.2 Weibull分布参数拟合极大似然法与矩估计法的取舍Matlab里拟合Weibull参数最方便的是wblfit函数它默认采用极大似然估计MLE。MLE的优点是渐近无偏且效率高在样本量足够大比如我的一年8760个点条件下估计结果非常稳定。调用方式% v_wind清洗后的风速序列行向量 [param_hat, param_ci] wblfit(v_wind); k_hat param_hat(2); % 形状参数 c_hat param_hat(1); % 尺度参数注意wblfit返回的第一个参数是尺度参数c第二个才是形状参数k这个顺序跟很多教材里的写法相反比较容易搞错。我一开始就栽在这个坑上画出来的拟合曲线跟直方图完全对不上检查半天才发现是参数顺序反了。除了MLE矩估计法在某些工程场景下也很实用。矩估计用风速均值μ和标准差σ来反推参数计算简单稳定性好。Weibull参数的矩估计公式为k (σ / μ)^(-1.086) c μ / Γ(1 1/k)其中Γ是伽马函数。当样本量较小或者数据噪声较大时矩估计比MLE更稳健。我建议在写代码时可以同时输出两种估计结果交叉验证一下拟合优度如果差异过大说明数据本身可能不符合Weibull假设需要重新审视。实测下来在风速数据质量较高的情况下MLE和矩估计得到的k值差异通常在0.1以内c值差异在0.5 m/s以内。3.3 Beta分布参数拟合利用mle函数定制估计Matlab没有直接的betafit函数但可以用betafit的底层——mle函数来实现。mle函数支持自定义分布类型核心代码如下% r_solar清洗后的归一化辐照度序列0~1之间 % 先把辐照度除以一个参考值比如1000 W/m²或者直接除以实测最大值 x r_solar / max(r_solar); % 使用mle函数指定分布为beta phat mle(x, pdf, (x, a, b) betapdf(x, a, b), start, [1, 1]); a_hat phat(1); b_hat phat(2);这里有几个细节必须注意。归一化基准的选择直接影响Beta分布参数的大小。如果选择max(r_solar)作为基准那么样本最大值恰好等于1这会使得Beta分布右端概率密度趋于无穷大拟合出的a参数偏小。更稳妥的做法是用一个物理意义清晰的常数比如标准测试条件STC下的辐照度1000 W/m²这样样本一般不会碰到1的上界拟合更平滑。另一个细节是start初始值要合理否则mle可能迭代不收敛。我通常设成[1,1]对应均匀分布的参数作为中性起点收敛速度还不错。3.4 拟合优度检验不止看图还要看数值指标很多初学者画完拟合曲线叠在直方图上看到大致吻合就觉得没问题。但在论文或工程报告里必须有定量的拟合优度检验。我常用的三个指标均方根误差RMSE衡量拟合概率密度与直方图密度的差异越小越好决定系数R²接近1说明拟合曲线解释了大部分数据变异性K-S检验p值p值大于0.05说明不能拒绝“数据来自该分布”的原假设。Matlab里K-S检验需要用kstest函数但注意它默认针对标准正态分布测试Weibull或Beta分布时需要先通过概率积分变换把数据转换成均匀分布再进行检验。或者直接用kstest2比较经验CDF和理论CDF。类似这样% 对风速做K-S检验 cdf_theory wblcdf(sort(v_wind), c_hat, k_hat); [h, p] kstest(sort(v_wind), CDF, [sort(v_wind), cdf_theory]);其中h0表示在5%显著性水平下不拒绝原假设p就是检验p值。实测一年数据风速的Weibull拟合K-S p值通常在0.1左右辐照度的Beta拟合p值在0.05~0.2之间波动说明模型对实际数据有解释力但并非完美。4. 组合采样与相关性控制从独立到耦合4.1 独立采样与联合出力时序生成最简单的组合模型是假设风速和辐照度相互独立。在这个假设下直接生成两组独立的随机数然后转成出力时序即可。Matlab代码如下% 设定仿真时长与时间分辨率 n_hours 8760; % 一年 dt 1; % 单位小时 % 生成风速样本服从Weibull分布 v_sample wblrnd(c_hat, k_hat, n_hours, 1); % 生成归一化辐照度样本服从Beta分布 r_sample betarnd(a_hat, b_hat, n_hours, 1); % 风速转风电出力 P_r_wind 1.5e6; % 风机额定功率W u_cin 3; u_r 12; u_coff 25; P_wind zeros(n_hours, 1); for i 1:n_hours if v_sample(i) u_cin v_sample(i) u_r P_wind(i) P_r_wind * (v_sample(i) - u_cin) / (u_r - u_cin); elseif v_sample(i) u_r v_sample(i) u_coff P_wind(i) P_r_wind; end end % 辐照度转光伏出力 P_r_pv 1e6; % 光伏额定功率W eta 0.18; % 综合效率 A_pv 5000; % 光伏阵列面积m² P_pv eta * A_pv * r_sample .* 1000; % 辐照度单位W/m²时1000为参考值 P_pv(P_pv P_r_pv) P_r_pv;独立采样模型能快速给出联合出力的边界范围但它的缺陷也很明显完全没有考虑天气过程引起的相关性。例如连续阴雨天风速可能上升而辐照度下降夏季午后风速可能增大而辐照度依然很强。这种独立性假设会让调度部门做出过于乐观的备用容量配置。4.2 基于秩相关的相关性控制方法如果要在模型中引入风速和辐照度的相关性最常见的方法是Iman-Conover方法也叫排序配对法。核心思想是先生成两个独立的随机序列然后通过调节序列的排列顺序使它们的秩相关系数Spearman相关系数逼近目标值。Matlab里实现起来不复杂target_rho 0.3; % 目标秩相关系数根据实测数据估计 % 生成目标相关结构用高斯copula生成中间相关序列 rho_s 2 * sin(pi * target_rho / 6); % 把Spearman秩相关转成Pearson相关近似 Z mvnrnd([0 0], [1 rho_s; rho_s 1], n_hours); % 转换回均匀分布 U normcdf(Z); % 用均匀分布通过逆CDF转成目标边缘分布 v_sample_corr wblinv(U(:,1), c_hat, k_hat); r_sample_corr betainv(U(:,2), a_hat, b_hat);关键在于目标相关性的获取。实际工程中直接把同时间段的风速和辐照度数据求Spearman相关系数作为建模输入。我处理过的那组数据Spearman相关系数约为-0.15也就是风速偏大时辐照度略偏低符合当地“大风往往伴随云系活动”的气候特征。所以组合模型中相关系数一定要体现这种负相关否则会高估联合出力的稳定性。如果你没有实测数据也可以根据地区气候粗略设定沿海多风地区风速与辐照度常呈弱负相关取-0.2到0内陆干旱地区晴热大风天气较多可能呈弱正相关取0到0.2。这个参数对联合出力分布形状影响很大建议做敏感性分析分别跑-0.3、0、0.3三组观察结果差异。4.3 出力时序的平滑与场景削减生成一年8760小时的联合出力序列之后直接用来做随机优化计算计算量太大而且很多相邻时刻出力相似信息冗余。更高效的做法是进行场景削减保留典型场景。Matlab里可以调用kmeans聚类或者K-medoids方法把8760个时段的“风速-辐照度”二维向量聚成若干典型日或典型时段。我习惯先把一天24小时看作一个场景向量那么一年就有365个场景再用K-medoids聚成10个典型日场景。这样可以保留季节与天气类型信息又大幅降低后续优化模型的维度。聚类前数据要做标准化因为风速和辐照度的量纲和数值范围差异很大不标准化的话聚类结果会偏向数值大的辐照度。类似这样% daily_scene365x24的矩阵每行是一天的小时序列 daily_scene_normalized zscore(daily_scene); % 标准化 idx kmedoids(daily_scene_normalized, 10, Distance, sqeuclidean);聚类完成后计算每个簇的质心作为典型日场景簇内样本数除以365作为该场景的概率权重。这样一套流程下来就生成了可用于随机规划的概率场景集。5. 仿真结果分析与可视化如何验证组合模型有效5.1 边缘分布拟合效果对比跑完拟合和采样后第一步验证就是看生成样本的边缘分布是否与原始数据一致。我通常在一张图上画三样东西原始数据的直方图、拟合的理论PDF曲线、以及蒙特卡洛生成样本的直方图。Matlab实现figure(Position, [100 100 900 400]); subplot(1,2,1); histogram(v_wind, 50, Normalization, pdf, FaceColor, [0.8 0.8 0.8]); hold on; v_axis 0:0.1:30; plot(v_axis, wblpdf(v_axis, c_hat, k_hat), r-, LineWidth, 2); histogram(v_sample, 50, Normalization, pdf, EdgeColor, b, FaceColor, none); legend(实测数据, Weibull理论, 采样样本); xlabel(风速 (m/s)); ylabel(概率密度); title(风速分布拟合); subplot(1,2,2); histogram(x, 50, Normalization, pdf, FaceColor, [0.85 0.85 0.85]); hold on; x_axis 0:0.01:1; plot(x_axis, betapdf(x_axis, a_hat, b_hat), r-, LineWidth, 2); histogram(r_sample, 50, Normalization, pdf, EdgeColor, b, FaceColor, none); legend(实测归一化辐照度, Beta理论, 采样样本); xlabel(归一化辐照度); ylabel(概率密度); title(辐照度分布拟合);理想状态下三条线实测直方图、理论曲线、采样直方图应该高度重叠。如果采样直方图明显偏离理论曲线多半是随机数生成时的样本量不够或者逆变换过程中数值溢出导致个别样本异常。5.2 联合出力概率密度与累计分布组合模型的核心输出是风电和光电联合出力的概率特性。先把两组出力加起来得到总出力序列P_total P_wind P_pv然后统计它的经验CDF。这个CDF就是调度部门最关心的东西——它直接回答了“在某个置信水平下系统总出力不低于多少”。比如我计算得到P90出力值含义是90%的时间里总出力超过这个值。这对安排旋转备用容量非常直观。Matlab里计算分位数很简单p_total P_wind P_pv; p90 quantile(p_total, 0.1); % 注意quantile默认是累积概率0.1对应90%分位点下方这里要小心quantile函数的细节。Matlab的quantile(x, p)返回的是下分位数即累计概率为p的值。所以P9090%概率超过的值对应quantile(p_total, 0.1)恰好是累计分布的10%分位点因为超过它的概率是90%。这个方向搞反的话结果完全不同。5.3 相关性对联合出力的影响分析我对比了独立模型和相关模型的联合出力分布发现一个有意思的现象当风速与辐照度呈负相关时联合出力的方差比独立模型更小出力波动更平缓。原因很简单风速大时辐照度小两者出力互补系统总出力不会同时偏低或偏高自然更稳定。而强正相关模型下高风和高辐照同时出现低风和低辐照也同时出现总出力的波动被放大。这个结论对储能配置很有价值。如果在负相关区域做风光互补项目储能容量可以适当降低因为系统天然具备互补性。反之如果正相关明显储能容量就要加大才能平抑联合出力的极端波动。6. 常见问题与排查技巧实录6.1 为什么wblfit拟合出的形状参数偏低我自己遇到过多次这种情况。检查发现要么是数据里混入了大量静风时段风速为0要么是风速序列的采样间隔过短导致风速变化相关性太强有效样本量不足。静风时段太多会让分布形态趋向指数分布k值偏小。解决方法是先做风速数据筛选把持续若干小时风速小于切入风速的时段单独处理或者直接剔除占比过高的极端静风段再重新拟合。6.2 Beta分布拟合时mle报错不收敛怎么处理mle在拟合Beta分布时偶尔会报“迭代超过最大次数”的错误。常见原因是归一化后的数据里包含0或1的边界值而Beta分布概率密度在边界处趋向无穷导致似然函数无穷大。处理办法是把数据略微向中心收缩比如x 0.001 0.998 * x; % 将数据范围从[0,1]移到[0.001, 0.999]这个变换对参数影响极小但能有效避免边界问题。另外增加mle的迭代次数设置或者换用betafit如果有拟合好也是一种选择。6.3 生成的风速样本出现负值或无穷值wblrnd在形状参数k非常小比如小于0.5时会生成极大值甚至超过风机切出风速导致模型失真。这通常是因为原始数据分布太厚尾Weibull假设本身就不适用。这时可以考虑改用混合分布或者先对数据进行截断处理只保留0到切出风速区间的数据来拟合。6.4 联合出力曲线在额定功率处出现“墙”状突变是否正常这是完全正常的。风电出力在额定风速到切出风速区间内恒等于额定功率导致总出力分布在该处出现概率质量集中。如果你在论文中画PDF曲线会看到在额定功率处出现一个尖峰这不是Bug而是物理特性。表述时要注意说明这一点否则审稿人会误以为代码出了错。6.5 如何缩短蒙特卡洛仿真的运行时间一年的时序仿真加上场景削减在普通PC上运行只需要几秒到十几秒瓶颈往往不在采样而在后续的优化计算。如果觉得for循环太慢可以用向量化操作替换。比如风速转出力那段可以直接用逻辑索引P_wind zeros(n_hours, 1); idx_partial v_sample u_cin v_sample u_r; idx_rated v_sample u_r v_sample u_coff; P_wind(idx_partial) P_r_wind * (v_sample(idx_partial) - u_cin) / (u_r - u_cin); P_wind(idx_rated) P_r_wind;这样速度可以提升数倍尤其在仿真时长扩展到多年数据时效果明显。7. 从模型到论文数据分析与成果呈现技巧7.1 图表展示的顺序与逻辑写论文或技术报告时组合模型的图表展示最好按照“数据特性→拟合效果→组合特性→应用示例”的顺序来安排。先展示原始风速和辐照度的统计直方图让人直观看到数据特征再展示拟合的Weibull和Beta分布曲线配合拟合优度指标表格接着展示联合出力样本的散点图或概率密度热力图体现相关性最后给出联合出力的CDF曲线和典型场景聚类结果落到应用端。7.2 关键参数应写成表格一张清晰的参数表比文字描述更高效。我在论文中通常放一个三列表格参数符号、含义、估计值。例如参数含义估计值kWeibull形状参数2.13cWeibull尺度参数 (m/s)7.85aBeta第一形状参数2.36bBeta第二形状参数3.12ρ风速与辐照度Spearman相关系数-0.15表格让审稿人一目了然也方便别人复现你的工作。7.3 代码可复现性的注意事项如果打算把Matlab代码附在论文或博客里一定要保证代码可以一键运行。常见的问题是数据路径写死成绝对路径或者依赖某个工具箱别人没有。我习惯把示例数据生成逻辑直接写在代码里比如先构造一个模拟风速和辐照度序列然后再替换成真实数据。这样读者不需要额外下载数据也能跑通流程。另外Matlab版本的兼容性也要注意。mle、wblrnd、betarnd这些核心函数在2016年以后的所有版本都能用但kmedoids需要Statistic and Machine Learning Toolbox如果没有这个工具箱可以用kmeans替代它们聚类效果接近只是在处理边界点时略有差异。8. 我个人实际操作中的几点体会做风电和光电联合出力的概率模型我最深的感受是分布拟合只是第一步真正的功夫在后续的“组合”环节。很多论文把大量的篇幅花在拟合Weibull和Beta参数上但对组合模型的核心难点——相关性的引入和场景生成——一笔带过这是本末倒置。实际上电网调度关心的不是风速本身服从什么分布而是功率输出和总出力的不确定边界。另一个常常被忽略的点是数据的季节性问题。用全年数据拟合一个统一的Weibull或Beta参数会掩盖不同季节分布形态的差异。比如冬季风大且稳定夏季风小且阵发性强这两季的形状参数k可能差出0.5以上。更精细的做法是分季节拟合参数或者以月为窗口滑动拟合看参数的季节性变化趋势。如果条件允许最好在模型里引入月份作为条件变量让分布参数随月份变化这样生成的风、光出力时序会贴近实际年周期波动。最后提醒一点Matlab里的随机数生成器默认状态是随机的每次跑出来的结果都不一样。如果你的工作涉及对比不同方案务必在开头加一句rng(2024)固定随机种子否则每次运行的结果无法复现往严重的说就是学术诚信问题。固定种子之后所有基于随机采样的结论都可以被他人复核这本身就是科研规范的一部分。这套代码和思路我已经在多个风、光互补项目预研里实际跑过结论具有一定普适性。如果大家想进一步扩展可以考虑把温度和光伏板温度系数加进来或者把风电功率预测误差建模成正态分布去跟风光联合分布叠加那又是另一个深度了。先从这版基础组合模型入手把每一步跑通、看懂、能改后面的路会顺很多。
