风电和光伏的随机出力特性一直是新能源建模绕不开的老大难。做微网容量配置、储能调度策略、甚至电网可靠性评估的时候第一步都得把风速和光照的统计规律摸清楚。风速普遍用两参数Weibull分布拟合光伏辐照度或出力则用Beta分布描述这俩在行业里几乎是标配。这篇内容我用Matlab把整套流程完整跑了一遍从参数估计、拟合优度检验到两个分布的组合抽样模拟全部实测可用。适合正在做新能源发电建模、分布式电源规划或者相关课程设计的朋友直接参考。1. 为什么偏偏是Weibull和Beta搞新能源的人大概都听过这俩分布但真要说清楚为什么是它们不少人是含糊的。我先把这个讲透后面代码才有意义。1.1 风速的Weibull模型双参数的物理直觉风速的统计分布靠什么大气边界层里的湍流运动本质上是大量随机扰动叠加的结果但风速有物理下界不低于0且右尾偏长所以正态分布根本不合适。Weibull分布恰好是极值分布的一种它能刻画最小值受限、存在长尾这类自然现象用在风速上既有理论依据又有大量实测支撑。两参数Weibull的概率密度长这样[ f(v) \frac{k}{c}\left(\frac{v}{c}\right)^{k-1}\exp\left[-\left(\frac{v}{c}\right)^k\right], \quad v \ge 0 ]这里形状参数k决定分布形态k1时退化为指数分布k2附近接近Rayleigh分布很多风资源评估教材默认用k2做粗略估算k在3~4之间时接近正态。陆上风电场常见k值在1.5~2.5之间沿海或高原开阔地形可能到3以上。尺度参数c跟平均风速强相关约等于风能资源丰沛程度的标尺。这两个参数的物理意义清楚了后面做估计和结果解读才不会跑偏。还可以换个角度理解风功率公式里风速是三次方关系所以哪怕分布拟合的绝对误差看着不大换算成功率后误差会被放大三倍。这就是为什么风资源评估对拟合精度格外敏感——你偷懒少拟合了长尾可能就把高风速段的发电量低估了。1.2 光伏出力的Beta分布归一化背后的逻辑光伏出力数据本身是0到装机容量之间的连续值辐照度受云层、大气衰减影响呈现出明显的偏态。Beta分布定义在[0,1]区间形状正好匹配这种有界、偏态的物理过程。更重要的是Beta分布只靠两个形状参数α和β就能适配各种峰形——左侧集中阴天多还是右侧集中晴天多都能刻画灵活性很高。实用里第一步是把光伏出力或辐照度归一化[ x \frac{P}{P_{\text{rated}}} ]注意这里分母要用额定装机容量很多人图省事用历史最大值这其实是埋雷一旦数据里出现极端高值归一化后一堆点挤在0.8~0.9区间Beta分布的右尾会被拉得特别离谱后续储能配置跟着偏高。这个坑我后面还会细讲。1.3 组合研究到底在研究什么把Weibull和Beta放一起通常有两层意思。第一层是风电-光伏联合出力概率模型风速和光照天然独立所以联合密度函数可以直接写成二者密度的乘积前提是做完独立性检验然后用蒙特卡洛抽样生成大量风电、光伏出力组合场景。这些场景是后面储能容量优化、并网消纳分析的基础输入。第二层是多能互补特性量化通过组合抽样算风电和光伏出力在不同时段的相关性、互补指数评估风光打捆外送的平滑效应。这类分析在区域级新能源规划里很常见。本文按第一层思路实现全套代码第二层的扩展也会给思路。2. 参数估计理论推导和Matlab实现参数估计是整篇的核心。用自带函数只是一句话的事但你要不知道背后在干什么结果错了都察觉不到。2.1 Weibull参数的极大似然估计极大似然估计MLE是首选方法因为它在大样本下渐进无偏且方差最小。对Weibull分布MLE没有解析解需要迭代。写下对数似然函数[ \ln L n\ln k - nk\ln c (k-1)\sum\ln v_i - \frac{\sum v_i^k}{c^k} ]对k求导置零整理后得到k的迭代式[ \frac{1}{k} \frac{\sum v_i^k \ln v_i}{\sum v_i^k} - \frac{\sum \ln v_i}{n} ]给定一个初值k就能算出新的k反复迭代到收敛。c则直接有显式解[ c \left(\frac{1}{n}\sum v_i^k\right)^{1/k} ]Matlab里我用自己写的迭代实现了一遍也对比了内置函数结果一致。自己实现的好处是能直观看到收敛过程。% 手动实现Weibull MLE迭代估计 % v: 风速数据列向量单位m/s已剔除异常值 v v(:); n length(v); lnV log(v); k 2.0; % 常见初值也可以用矩估计先粗算 for iter 1:100 vk v.^k; sumVk sum(vk); sumVkLnV sum(vk .* lnV); meanLnV mean(lnV); k_new sumVk / (sumVkLnV - meanLnV * sumVk); if abs(k_new - k) 1e-6 k k_new; break; end k k_new; end c (mean(v.^k))^(1/k); fprintf(Weibull: k %.4f, c %.4f\n, k, c);初值选k2一般没问题但如果数据方差很大迭代可能会震荡可以加个阻尼(k_{\text{new}} 0.5(k_{\text{old}} k_{\text{new}}))更稳。样本量少于50时MLE会有偏置工程上建议至少上百个数据点再做估计。直接用内置函数也可以pd_w fitdist(v, Weibull); % 注意a是形状kb是尺度c % 或 [param, ci] wblfit(v); % 返回[尺度c, 形状k]及其置信区间wblfit返回顺序和fitdist正好相反我每次都会踩这个记忆坑。建议代码里注释写清楚避免后续维护的人看晕。2.2 Beta参数的矩估计与极大似然估计Beta分布的密度函数是[ f(x) \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} ]矩估计有解析解算起来飞快。设归一化出力样本的均值为m、方差为s²令中间量 (r \frac{m(1-m)}{s^2} - 1)则[ \alpha m r, \quad \beta (1-m) r ]这组公式是统计书里的标准结果胜在快且稳定适合做初值。MLE则没有显式解Matlab内置的betafit用的就是数值MLE建议直接用。% Beta分布参数估计 % x: 归一化光伏出力范围(0,1) x x(:); x(x 0) 0.001; % 边界处理避免0和1导致的无穷大 x(x 1) 0.999; % 矩估计做初值 m mean(x); s2 var(x); common m * (1 - m) / s2 - 1; alpha0 m * common; beta0 (1 - m) * common; % 内置MLE phat betafit(x); alpha phat(1); beta phat(2); fprintf(Beta: alpha %.4f, beta %.4f\n, alpha, beta);边界处理很关键。归一化后碰上0或1是常见事——多云天气下光伏出力本来就是0深夜检修时也全是0。不做处理log(0)直接给你Inf估计结果整个完蛋。我用的是简单截断法业界还有一种做法是加微小抖动短时间序列里截断法更简单可控。2.3 拟合优度怎么判断才靠谱参数估出来不代表拟合得好。我习惯三个指标一起看避免单一指标骗人。决定系数R²是基础计算经验CDF和理论CDF的差值平方和。Matlab里可以这么算% 计算R² v_sort sort(v); n length(v_sort); emp_cdf (1:n) / n; % 经验CDF theo_cdf wblcdf(v_sort, c, k); % 理论CDF residual emp_cdf - theo_cdf; R2 1 - sum(residual.^2) / sum((emp_cdf - mean(emp_cdf)).^2);偏态分布变异很大R²低于0.95就得重新检查数据或考虑混合分布。对有经验的人来说Q-Q图比R²更直观尾部是否有系统性偏差一眼就能看出来figure; p (1:n) / (n 1); quantile_theo wblinv(p, c, k); quantile_emp sort(v); scatter(quantile_theo, quantile_emp, 10, filled); hold on; plot([min(v), max(v)], [min(v), max(v)], r--, LineWidth, 1.5); xlabel(理论分位数); ylabel(样本分位数); title(Weibull Q-Q图);KS检验是严格的统计检验但我要先给个大警告样本量一大比如超过几千KS检验基本必然拒绝原假设。这不是你的拟合差而是检验本身在大样本下对微小偏差极度敏感。所以大规模数据下我更推荐看Q-Q图和R²KS检验只做辅助参考。[h, pval] kstest(v, CDF, [v_sort, theo_cdf]);记住这里的CDF输入形式需要两列矩阵第一列是数据点第二列是对应的理论CDF值。3. Matlab完整实现从数据预处理到组合抽样这一节给整套可直接跑的流程按顺序操作即可。我用的数据是某风电场一年的逐小时风速和光伏电站同时段的逐小时出力模拟数据也兼容关键是把算法跑通。3.1 数据准备与预处理拿到原始数据第一件事是清洗。风电场测风塔数据经常会有这几个毛病负值传感器故障、连续重复值通信中断补数、飞点雷电干扰。我的处理规则负值和NaN直接删除超过4倍标准差的点视为飞点单独标记并复核不直接删——有可能真是极端天气删了长尾就没了零值保留那是静风是物理真实光伏数据要检查的则是出力超过额定容量这类记录可能是电流互感器标定错了。另外要注意时区光伏出力跟着太阳走北京时间下中午12点的出力峰值才是正常的。% 数据加载 data readtable(wind_pv_data.csv); % 假设有wind和pv两列 v_raw data.wind; p_raw data.pv; % 风速清洗 v v_raw(~isnan(v_raw) v_raw 0); % 飞点检测可选按实际情况决定是否剔除 mu_v mean(v); sigma_v std(v); outlier_idx v mu_v 4 * sigma_v; fprintf(风速飞点数量: %d\n, sum(outlier_idx)); % 建议人工复核后再决定是否删除 % v(outlier_idx) []; % 光伏出力归一化 P_rated 5; % 单位MW额定容量 x p_raw / P_rated; x x(~isnan(x)); x max(0.001, min(0.999, x)); % 边界截断P_rated这个值务必从电站设计资料里查千万别用max(p_raw)替代。我之前接过一个项目数据里面有几天出力异常超过额定值用max归一化后整个Beta分布被拉到荒谬的位置存到配置参数的置信区间全都失真。3.2 分布拟合完整代码块把前面两节的估计代码整合成一个脚本输出参数、置信区间和拟合指标% Weibull拟合 pd_w fitdist(v, Weibull); k_hat pd_w.a; % 形状参数 c_hat pd_w.b; % 尺度参数 ci_w paramci(pd_w); % 参数置信区间 % Beta拟合 phat betafit(x); alpha_hat phat(1); beta_hat phat(2); % 置信区间可以自己用bootstrap跑样本多时可以自动我实测下来fitdist的风速拟合结果非常稳定数值上和手写MLE一致。置信区间方面paramci默认给95%区间样本量小时区间很宽说明参数不确定度大这时候别急着用于下一级计算最好多拿点历史数据。3.3 绘图概率密度叠加直方图图是给评审和论文用的画得规范很重要。我用的是直方图加理论密度曲线的叠加方式figure; histogram(v, Normalization, pdf, NumBins, 40, FaceColor, [0.8 0.8 0.8]); hold on; v_plot linspace(0, max(v)*1.05, 200); plot(v_plot, wblpdf(v_plot, c_hat, k_hat), r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(风速直方图, Weibull拟合); title(风速分布拟合); figure; histogram(x, Normalization, pdf, NumBins, 30, FaceColor, [0.9 0.7 0.7]); hold on; x_plot linspace(0.01, 0.99, 200); plot(x_plot, betapdf(x_plot, alpha_hat, beta_hat), b-, LineWidth, 2); xlabel(归一化光伏出力); ylabel(概率密度); legend(出力直方图, Beta拟合); title(光伏出力分布拟合);直方图NumBins的选择有讲究。数据量2000以上时40个bin能看出长尾形态数据量只有几百的时候bin太多会毛刺严重反而不如20个bin平滑。也可以用Matlab的histogram自动分箱但自动算法偏向保守细节展示不够好。3.4 组合抽样蒙特卡洛场景生成这是最有实用价值的部分。有了两个边缘分布接下来用蒙特卡洛生成大量风速-辐照度场景对。因为风速和光照之间的物理机制相互独立直接抽样再组合就行除非你研究的是气温对光伏效率的耦合影响那是另一回事。% 蒙特卡洛组合抽样 N 50000; % 抽样场景数量 % 独立抽样 v_sample wblrnd(c_hat, k_hat, N, 1); % 注意wblrnd参数顺序是(c, k) x_sample betarnd(alpha_hat, beta_hat, N, 1); % 反归一化得到光伏出力场景MW p_sample x_sample * P_rated; % 风电出力转换基于风机功率曲线简化为分段函数 % v_in3m/s, v_rated12m/s, v_out25m/s, P_w_rated1.5MW每台 P_w_rated 1.5; % MW单台风机额定功率 v_in 3; v_rated 12; v_out 25; P_w_sample zeros(N, 1); idx1 v_sample v_in v_sample v_rated; idx2 v_sample v_rated v_sample v_out; P_w_sample(idx1) P_w_rated .* (v_sample(idx1).^3 - v_in^3) / (v_rated^3 - v_in^3); P_w_sample(idx2) P_w_rated; % 超出切出风速或在切入风速以下出力为0初值 % 组合出力 P_total_sample P_w_sample p_sample; % 统计特征输出 fprintf(组合出力均值: %.4f MW\n, mean(P_total_sample)); fprintf(组合出力标准差: %.4f MW\n, std(P_total_sample)); fprintf(组合出力95%%分位数: %.4f MW\n, quantile(P_total_sample, 0.95));这段代码的功率曲线是简化版实际工程里风机有偏航、变桨控制功率曲线是S形不是严格的三次方关系但用来做规划层面的出力场景模拟足够了。50kW以下的小微系统可以改用一个近似公式风速低于切入风速时出力为0高于额定风速时出力恒定这样更贴合实际运行曲线。为了更好理解组合结果可以画散点图和出力累计分布figure; scatter(v_sample, p_sample, 2, filled, Alpha, 0.3); xlabel(风速 (m/s)); ylabel(光伏出力 (MW)); title(风-光出力场景散点分布); figure; cdfplot(P_total_sample); hold on; cdfplot(P_w_sample); cdfplot(p_sample); legend(组合出力, 风电出力, 光伏出力); xlabel(出力 (MW)); ylabel(CDF); title(出力累积分布对比); grid on;散点图能直观看出两个分布的形状配合组合出力的CDF则是后续做概率潮流、可靠性评估时直接要用的输入。场景数量N建议5万起步少了尾部分位数不稳定多了计算量成倍涨5~10万是性价比最高的区间。4. 常见问题与排查技巧实录实操中遇到的那些坑我一个个列出来都是真实踩过的。4.1 风速数据里静风太多Weibull拟合直接失效某沿海项目测风数据里有30%以上的0值直接丢进fitdist输出的k值小于1概率密度在v0处出现无穷大峰。这明显不符合风电场实际——真正静风时段没有那么长很多0其实是风速低于传感器启动阈值比如低于0.5m/s就记0这是传感器灵敏度问题。处理办法不是简单删0而是做零截断处理清除静风时段后再对非零部分拟合同时单独统计静风概率p0作为混合分布的一部分[ f_{\text{mix}}(v) (1-p_0) f_W(v) p_0 \delta(v) ]这样做出来才是完整的概率模型。实际应用中静风时段往往集中在夜间或特定季节如果只是简单删0拟合出的c值会偏高风能资源评估会乐观10%以上。4.2 Beta拟合发散归一化边界没处理好前面写的max(0.001, min(0.999, x))不是随便拍的。β分布的似然函数里包含(1-x)的β-1次幂如果x1直接算就是0的负几次幂必然Inf。边界截断到0.001和0.999log部分就不会爆炸。但要注意截断幅度太大会改变分布形状。如果你的数据真的有一大堆0值比如夜间出力全为0单纯截断会扭曲估计结果。这时候应该用零膨胀Beta分布Zero-Inflated Beta把0值单独建模而不是混在一起硬拟合。行业里有建议是先用聚类把数据分为白天和夜间时段分别拟合白天时段用Beta夜间单独统计0值和极低出力概率。4.3 kstest大样本必拒绝的陷阱5000个样本以上kstest几乎一定会拒绝原假设p值趋近于0。很多初学者看到这个结果就慌了以为拟合失败其实不然。KS检验的检验功效随样本量增加而增强大样本下任何实际分布与理论分布的微小偏差都会被检测出来。真实数据不可能完美服从理论分布所以不必对这个结果过度解读。我的判断经验1000个样本以内KS检验有参考价值超过5000主要看R²和Q-Q图。R²在0.94以上Q-Q图的点分布在参考线附近没有系统性弯曲这个拟合在实际工程里已经可以用了。4.4 参数的季节漂移问题一年统拟合的结果春夏秋冬分别拟合会差很多。风速的k和c值冬季明显偏高风大且稳定夏季则偏低光伏Beta分布的α和β也随季节显著变化夏天的α和β都大分布更集中冬天则偏态明显。做全年调度策略时用一个年统参数没问题但如果是做季度容量配置或者分时电价下的储能充放策略建议按季度分别拟合四个参数组。这个改动很简单用分组训练循环就行% 按季节分组拟合 seasons {春季, 夏季, 秋季, 冬季}; for s 1:4 idx (month(data.time) (s*3-2)) (month(data.time) s*3); v_s v(idx); x_s x(idx); % 分别做拟合检验 end4.5 抽样结果和实测统计对不上组合抽样得到的出力分布和这一年实际出力分布差距明显先别怀疑蒙特卡洛不对。多半是参数估计阶段的数据处理埋了问题。排查步骤先看拟合阶段R²是否达标再看功率曲线模型是否符合实际风机型号最后核实时段是否匹配——用夜间数据拟合Beta分布自然会偏。还有一种常被忽略的情况光伏数据的额定容量和实际出力单位不对应比如做了单位换算kWh和MWh但忘记除以1000。这种低级错误在项目里非常常见我处理过至少三个类似咨询最后都是单位搞错。5. 后续扩展思路与工具箱推荐做完基础组合研究还有几条明显的扩展路径。考虑时空相关性的Copula扩展。独立假设在天气系统尺度上不一定成立相邻区域的风速往往强相关。用t-Copula或Gaussian Copula把两个边缘分布耦合起来能做出区域风电集群出力模型。Matlab里有copulafit和copularnd系列函数调用起来很方便。从分布到时序用马尔可夫链加分布抽样。把风速状态离散化比如0~25m/s划分成10个状态估计状态转移矩阵每一步从对应状态的Weibull分布里抽一个值就能生成带时序相关性的风速序列比直接蒙特卡洛近真实。置信区间在工程决策中的应用。前面paramci给出的参数置信区间可以直接用于储能容量配置的鲁棒优化——用保守参数组合跑一遍、乐观参数组合跑一遍得到容量配置的上下界这对工程投标报价非常有价值。工具方面Matlab本身很够用。Statistics and Machine Learning Toolbox提供了fitdist、betafit、wblrnd、betarnd、copulafit这些核心函数。我在多个版本上都跑过R2019b到R2024aAPI基本稳定老项目迁移没有遇到大问题。优化工具箱可以用来做风速分布的分段拟合不过一般用不上。如果哪天出了性能问题用parfor做多线程抽样也能很快提速。结合我自己的使用经验用Matlab做这类随机建模的优势是生态完善、出图方便本科甚至专科的学生都能快速上手缺点是大型场景模拟百万级抽样还是有点慢。作为替代可以考虑Python的SciPy和NumPy但Matlab在交互式探索和可视化上更顺手两种工具各有用途。我在实际项目中遇到的一个关键教训是参数估计只是工具数据清洗才是真正耗时的地方。70%的时间花在数据检查上一点不夸张。拿到数据先画出散点图肉眼扫一遍有没有飞点、缺段、单位异常再交给算法处理。算法只是把数据里已有的规律提取出来数据本身脏了再高明的拟合也救不了。最后再分享一个小技巧做完参数估计后用生成的分布抽样生成一组模拟数据把模拟数据和原始数据放在同一张图里做分位数对比。这个方法对检查模型是否忠实于原始数据极其直观比任何统计指标都有说服力。在风电和光伏的联合出力场景分析里这一步做好了后面所有优化和决策才能站得住脚。