风电和光电作为当前新能源装机的两大主力各自的波动性和随机性一直是规划与调度绕不开的难题。这个项目把风速的Weibull分布和光照的Beta分布放进同一个Matlab框架里做组合研究本质上是用两把“概率尺子”去量化风、光资源的统计特征再通过蒙特卡洛模拟把两种随机过程叠加起来评估联合出力特性。我早期做风资源评估的时候只拟合单一分布看不出互补性后来把光资源也拉进来一起建模很多东西才真正对上了。这篇博文我会从理论公式、参数估计讲到代码实现和踩坑记录适合电力系统方向的研究生、做新能源规划的工程师以及所有想用Matlab复现风-光组合分析的读者。代码都是实际跑过的直接复制到你的工程里就能用。1. 项目整体思路拆解两把尺子量同一个问题1.1 为什么非要用概率分布刻画风、光资源风电场和光伏电站的运行数据摆在那里最直观的感受就是“忽高忽低”。今天平均风速6m/s明天可能变成3m/s同一块光伏板上午出力200kW下午云一过来直接掉到50kW。如果只拿平均值去做发电量预测或者容量配置误差会大到让你怀疑人生。概率分布的引入就是为了把这些波动性“收编”成一组可以计算的参数。Weibull分布在风速拟合上几乎是行业标准原因在于它的形状参数k非常灵活——k小的时候分布偏斜能刻画多静风或者多强风的地区k大的时候分布收窄对应风速比较稳定的场址。Beta分布则完全不同它定义在[0,1]这个有界区间上天然适合处理归一化之后的光照强度或者光伏功率既不会出现负值也不会突破物理上限。组合研究的第一步就是先把这两种分布分别拟合好让它们各自成为描述风、光资源的两把可靠标尺。1.2 两个概率模型的数学基础Weibull分布的概率密度函数表达式为f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中v是风速k是形状参数c是尺度参数。实际工程里经常提到的“平均风速”和“有效风速区间”其实都可以从k和c里推导出来。累积分布函数F(v) 1 - exp(-(v/c)^k)在后续做KS检验和场景生成时都会用到。这里特别强调尺度参数c和平均风速不是一回事c通常略大于平均值两者之间还差着一个gamma函数的系数。Beta分布的表达式为f(x) x^(a-1) * (1-x)^(b-1) / B(a, b)其中x落在[0,1]区间a和b是形状参数B(a,b)是Beta函数。这个分布最妙的地方在于当a和b都小于1时它是U形的对应“要么没光、要么满发”的极端天气当a和b都大于1时它是钟形的对应光照相对平稳的场景。做光资源拟合时Beta分布的形态自适应能力非常强。1.3 组合研究的核心目标是什么单看风或者单看光都只能回答“这个资源怎么样”。组合研究要回答的则是“风光在一起到底靠不靠谱”。比如某个地区风大的季节正好光照弱那风光的联合出力可能比单独看任何一种资源都稳定反过来如果风和光总是同强同弱那联合出力的波动反而可能被放大。这个项目通过构建联合概率密度和生成大量随机场景把风速和辐照度映射到功率空间再统计联合出力的概率分布从而量化互补性、评估电力不足风险。这才是整个项目的灵魂所在。2. 从数据到参数Weibull拟合的完整过程2.1 风速数据的清洗与预处理我见过太多人拿到测风塔数据就直接往拟合函数里丢结果拟合出来的曲线惨不忍睹。真实的风速数据里通常混着三类“脏数据”一是缺测值时间戳上就是NaN二是超出物理范围的值比如风速82m/s这种明显是传感器故障的记录三是大量零风速段这些静风数据如果占比过高会严重压低形状参数k使拟合曲线低估高风速段的概率。我的清洗流程是先把风速列单独提取出来用isnan筛选缺测值超过物理范围的值直接置为NaN再统一插值对于静风记录如果占比超过总数据的5%我会考虑用零均值的小噪声去替代避免拟合时出现奇异行为。插值我习惯用线性插值fillmissing简单稳定。清洗之后一定要再画一次时序图确认数据形态合理再进入参数估计环节这一步能帮你省掉后面80%的麻烦。2.2 极大似然估计、矩估计与最小二乘的取舍Weibull参数估计最主流的方法是极大似然估计MLE。Matlab中提供了两个入口一个是直接调用wblfit函数返回形状参数和尺度参数另一个是把数据封装成fitdist对象再取对象的a和b字段。两者的数学本质一样但fitdist在后续配合pdf、cdf等函数时更方便。MLE的迭代公式看起来很简单实际求解时要用数值方法。Matlab内部处理了这部分但默认初值在数据量较小时偶尔会失效。我自己的习惯是先用矩估计算一个初值k ≈ (std/mean)^(-1.086)超参数c mean / gamma(1 1/k)然后把这个初值传给mle函数手动迭代。这样做的原因很实际——MLE的本质是求解非线性方程初值离真实值太远时牛顿迭代法可能发散或者收敛到局部解。矩估计虽然效率不如MLE但胜在稳定用它做初值几乎没出过问题。最小二乘拟合是另一种思路对累积分布函数做双对数变换把Weibull的拟合问题变成线性回归问题。具体做法是对ln(-ln(1-F))关于ln(v)做直线拟合斜率就是k截距和c有关。这个方法在早期文献里很常见现在的工程实践中更多用来做交叉验证。如果MLE拟合出来的结果和最小二乘的结果差距很大那大概率是数据本身有问题而不是方法的问题。2.3 拟合优度检验怎么才算过拟合完参数不能直接就看图说“差不多”。我在项目中至少做两件事第一用KS检验判断拟合分布和经验分布是否来自同一总体第二对比高风速段的尾部拟合效果。KS检验在Matlab里的用法有一个容易踩的坑。kstest函数要求传入CDF函数句柄和参数很多人直接写成kstest(data, CDF, pd)结果是错的正确的写法是kstest(data, CDF, {wblcdf, k, c})。这里的参数cell数组结构必须严格匹配否则Matlab会报错或者给出错误结果。检验输出的p值如果大于0.05说明在95%置信水平下不能拒绝原假设拟合通过。但KS检验也不是万能的它对分布中部的偏差敏感对尾部的偏差反而不敏感。而风速分布尾部恰恰对应着高风速、大出力的区间是发电量评估的关键。所以我每次都会把累积概率在0.9以上的数据点单独画出来看如果尾部拟合偏差超过10%我会考虑改用混合Weibull分布或者截尾Weibull而不是硬着头皮用一个分布去拟合到底。3. 光照强度与Beta分布的建模细节3.1 为什么Beta分布比正态分布更适合光资源很多人习惯性地对任何数据都用正态分布去拟合但光资源数据有一个非常硬性的约束物理上有上下界。辐照度不可能是负值也不可能无限大。如果用正态分布去拟合光强数据拟合结果的左尾会有一部分概率落在负值区间这在物理上是无意义的。Beta分布定义在[0,1]区间只要先做归一化就不会出现越界问题。另一个微妙的地方是光资源的分布形态会随天气类型剧烈变化。晴天时归一化辐照度集中在一个偏高的值附近分布呈尖峰状阴天时分布变得扁平甚至偏向0。Beta分布通过a、b两个参数可以覆盖这两种极端形态这是正态分布做不到的。Beta分布也有两个参数但这两个参数刻画的是“有界区间内的形状”和光资源的物理特征完全对齐。3.2 归一化处理理论最大值还是实测最大值用Beta分布拟合光资源数据之前必须先把原始辐照度数据归一到[0,1]之间。这一步看似简单实际暗藏一个关键决策归一化因子到底取什么。如果取理论最大值比如太阳常数对应的天文辐射上限这样得到的归一化数据几乎很难达到1Beta拟合出来的均值会整体偏低参数形态失真。如果取实测历史最大值拟合精度会明显提升但模型迁移到另一个地区时就不适用了。我在项目中通常采用折中方案先看数据质量如果实测序列足够长至少一年就用实测最大值的某个百分位数比如99.5%分位数做归一化因子既能保留拟合精度又能避免极端峰值对整体归一化的干扰。还有一类特殊情况——夜间零辐照度数据。Beta分布对x0和x1这两端的概率密度是发散的如果在拟合前不处理这些零值会导致参数a趋近于0拟合失败。我的做法是把辐照度低于某阈值的记录单独剔除或者在归一化时给零值加一个非常小的正数偏移。3.3 Beta参数估计与Matlab代码实现Beta分布的参数估计在Matlab里直接调用betafit就能完成函数内部也是用MLE框架。但为了给读者一个更直观的理解我通常先用矩估计算初值。矩估计的公式很简洁a μ * (μ*(1-μ)/σ² - 1) b (1-μ) * (μ*(1-μ)/σ² - 1)其中μ是样本均值σ²是样本方差。如果你发现计算出来的a或者b小于等于0说明这组数据可能不是一个单峰的Beta分布能描述的需要检查数据处理环节。对于大多数实际数据矩估计给出的初值已经足够接近真实值。下面是完整的Matlab拟合代码包含数据归一化和参数估计% 导入辐照度原始数据单位W/m^2 irradiance ...; % 剔除夜间零值保留白天有效数据 valid irradiance 10; % 阈值为10W/m^2 data_day irradiance(valid); % 归一化处理使用99.5%分位数作为归一化因子 norm_factor prctile(data_day, 99.5); x data_day / norm_factor; x(x 0) 1e-4; % 避免零值导致边界发散 % 矩估计初值 mu mean(x); sigma2 var(x); temp mu * (1 - mu) / sigma2 - 1; a_init mu * temp; b_init (1 - mu) * temp; % MLE拟合 phat betafit(x, [], [a_init, b_init]); a phat(1); b phat(2); % 拟合效果画图 figure; histogram(x, 30, Normalization, pdf, FaceColor, [0.8 0.8 0.8]); hold on; xx linspace(0, 1, 200); plot(xx, betapdf(xx, a, b), r-, LineWidth, 2); xlabel(归一化辐照度); ylabel(概率密度); legend(实测直方图, Beta拟合曲线);注意betafit的第二个参数是初始值第三个参数是优化选项。如果直接不给初始值Matlab会内部自动选取初值但实测下来给一个合理的初值能让收敛更快、结果更稳定。3.4 从辐照度到光伏出力的转换模型Beta分布拟合出来的是辐照度的统计特征但最终评价光伏出力还需要一个转换环节。工程中常用的是两参数分段模型P Pmax * (R - Rc) / (Rmax - Rc)当 R Rc P 0当 R ≤ Rc其中R是辐照度Rc是启动阈值一般取120-200W/m²Rmax是参考辐照度Pmax是装机容量。这个模型比线性模型多了一个“死区”能更真实地反映光伏组件在低辐照度下近乎不出力的物理特性。在组合仿真时先用betarnd生成归一化辐照度样本乘以归一化因子还原成辐照度再通过这个转换模型得到功率序列最终统计的光伏出力分布就会非常贴近实际。4. 风光组合仿真的Matlab代码实现4.1 蒙特卡洛场景生成的基本流程组合研究最核心的代码实现是蒙特卡洛模拟。整体思路是这样的先用拟合好的Weibull参数生成风速样本用Beta参数生成归一化辐照度样本再分别转换成功率最后把两种功率叠加起来形成联合出力序列。这里有一个重要的前提假设要说明我在基础版实现里假定风速和辐照度在小时尺度上是相互独立的。这个假设并不完全成立尤其在某些特定气候区风和光可能存在弱相关。但作为第一版组合研究独立假设已经能给出很多有价值的结论。想考虑相关性的话可以引入Copula函数做联合抽样相当于给两个边缘分布加一个相关结构这个在后续进阶版本里再详细展开。由于蒙特卡洛需要大量样本我建议用向量化操作而不是循环。生成10万个样本在Matlab里就是一条命令的事运行时间只需要零点几秒% 设置随机种子保证实验可复现 rng(42); % 生成风速场景样本 n_samples 100000; wind_speed wblrnd(k, c, n_samples, 1); % 生成归一化辐照度场景样本 norm_irradiance betarnd(a, b, n_samples, 1); % 反归一化得到真实辐照度 irradiance norm_irradiance * norm_factor;wblrnd的第一个参数是形状参数k第二个是尺度参数c和wblfit的输出顺序一致这个细节写代码时很容易搞反。betarnd的两个参数就是Beta分布的形状参数a和b。4.2 风功率和光伏功率的转换计算风速样本生成之后要借助风功率曲线把风速映射到风机出力。风电功率曲线一般是S形切入风速以下出力为0额定风速以上出力被限制在额定功率。用一个分段函数来表示% 风功率曲线参数 v_in 3; % 切入风速单位m/s v_out 25; % 切出风速 v_rated 12; % 额定风速 P_w_rated 1.0; % 标幺值取额定功率为1 % 计算风功率 P_wind zeros(n_samples, 1); idx1 wind_speed v_in wind_speed v_rated; idx2 wind_speed v_rated wind_speed v_out; P_wind(idx1) (wind_speed(idx1)^3) / (v_rated^3 - v_in^3) - (v_in^3) / (v_rated^3 - v_in^3); P_wind(idx2) P_w_rated;这里的功率曲线模型是简化版实际项目里应该用风机厂商提供的实测功率曲线表。但简化模型对于算法验证和分布组合研究已经足够了因为核心目的是展示随机变量通过非线性变换后的分布变化规律。光伏功率转换就用前面提到的分段模型% 光伏转换模型 R_c 150; % 启动阈值W/m^2 R_max norm_factor; % 参考辐照度 P_pv zeros(n_samples, 1); valid_idx irradiance R_c; P_pv(valid_idx) (irradiance(valid_idx) - R_c) / (R_max - R_c); P_pv min(P_pv, 1); % 限制最大出力不超过额定值两组功率都算出来之后定义联合出力P_total P_wind P_pv; % 风光联合出力标幺值这里的P_total就是每次随机场景对应的联合出力。联合出力的上限是2风、光各1下限是0。有了这10万个联合出力样本后面所有的统计分析都顺手了。4.3 联合概率密度与互补性分析的可视化联合出力样本生成之后最直观的展示方式是二维直方图。用Matlab的histogram2函数可以把风速和辐照度的联合分布画成一个三维地形图颜色深浅代表概率密度的高低。如果数据集中在对角线附近说明风、光存在同涨同落的相关性如果数据分布呈现“反对角”形态则说明互补性强——风速大的时候光照小风速小的时候光照大。另一个更实用的输出是联合出力的累积分布函数。用ecdf函数可以直接画出经验CDF曲线从这条曲线上可以读出工程上非常有价值的指标比如P_total小于某个阈值的概率。假设电力系统要求联合出力必须大于0.4标幺值如果从CDF上读出这个概率大于某个可接受值那就说明需要配置储能或者扩大装机容量。这就是组合研究从理论走向工程决策的落点。互补性指标我一般用Pearson相关系数和Spearman秩相关系数同时计算corr_pearson corr(wind_speed, irradiance, Type, Pearson); corr_spearman corr(wind_speed, irradiance, Type, Spearman);Pearson相关系数描述线性相关Spearman描述单调相关性。如果两者符号一致都能反映互补趋势如果符号相反说明两个变量之间存在非线性关系这时候单纯用相关系数做判断就不够了需要进一步借助Copula模型的尾部相关系数来识别极端场景下的互补性。5. 踩坑记录与排查技巧实录5.1 参数初值不当导致的拟合失败我在做Beta分布拟合时遇到的问题最典型。有一组光照数据来自某高海拔地区晴天极多归一化辐照度集中在0.8-0.95之间方差非常小。直接用betafit拟合结果返回了离谱的参数a50b2.5拟合曲线尖得几乎像一根线。后来排查发现是矩估计初值中的σ²出现了下溢导致a_init计算异常。解决办法是对原始数据先做一次标准化检查如果方差过小给矩估计公式中的方差项设置一个下限或者直接改用极大似然估计的默认初值。这个坑给我的教训是不要想当然地认为可调函数总能自适应一切数据理解初值的作用机制才能在异常结果出现时快速定位问题。5.2 零值数据和极端值对拟合的干扰风光数据的共同特点是大量零值。风速数据中的静风段表现在时域上就是一段平直为零的曲线光资源数据中的夜间时段则是另一个零值来源。如果不加处理直接喂给拟合函数Beta分布会因为x0处的概率密度发散而完全崩溃Weibull分布中的零值则会让形状参数k被拉到极低值。处理零值我总结了三条经验第一风资源分析中保留静风段但拟合前先计算静风占比占比过高时要考虑两参数Weibull是否还适用第二光资源分析中直接剔除夜间零值只拟合白天的有效数据第三对归一化之后出现的零值统一加上一个1e-4到1e-3量级的小偏移避免数值问题。这个小偏移看似不起眼却能拯救整个拟合过程。5.3 随机种子的设置与结果复现问题蒙特卡洛仿真如果每次运行结果都不一样实验就失去了可复现性。这个问题在论文审稿和工程评审时尤其致命。第一次跑出功率不足概率是8.5%第二次跑变成9.1%评审直接质疑数据的可靠性。解决方法是严格使用rng函数设置随机种子而且是每次生成随机样本之前都设置一次。另外要特别注意wblrnd和betarnd如果每次调用前都重置种子确实能保证序列完全一致但同时意味着统计上每次都是同一组伪随机数在某些极端情况下可能引入偏差。在实践上我通常生成一次大规模样本后保存到.mat文件后续所有分析都从文件中读取样本既保证可复现又方便不同参数方案的横向对比。5.4 大样本仿真时的内存与速度优化10万样本对Matlab来说是小菜一碟但如果你做的是8760小时的全年场景再加上几十年的历史数据回测内存占用就会变成一个大问题。一个常见的错误是频繁在循环里调用wblrnd和betarnd每次生成一个样本再拼接。正确的做法是提前一次性生成所有样本用向量化计算代替循环。如果样本量实在太大考虑分批生成并增量累加统计量而不是把全部样本都存在内存里。我在处理百万级样本量时还会利用parfor做并行计算配合数据分块策略运行时间可以减少60%以上。我个人做风光组合研究走下来最大的体会是拟合分布只是整个项目的地基真正出成果的地方在组合分析和场景生成。Weibull和Beta的公式本身并不复杂Matlab里对应的拟合函数也都是现成的但如果忽略了数据清洗、参数初值、零值处理这些细节后面每一步都会带着偏差走。建议读者拿到项目后先花一半时间把两份真实数据彻底清理干净再把参数估计和组合仿真的代码跑通最后用一组已知结果的数据做验证。这样下来你的风光联合出力分析模型才算真正立得住。
风光联合出力概率建模：Weibull与Beta分布组合的Matlab实现
风电和光电作为当前新能源装机的两大主力各自的波动性和随机性一直是规划与调度绕不开的难题。这个项目把风速的Weibull分布和光照的Beta分布放进同一个Matlab框架里做组合研究本质上是用两把“概率尺子”去量化风、光资源的统计特征再通过蒙特卡洛模拟把两种随机过程叠加起来评估联合出力特性。我早期做风资源评估的时候只拟合单一分布看不出互补性后来把光资源也拉进来一起建模很多东西才真正对上了。这篇博文我会从理论公式、参数估计讲到代码实现和踩坑记录适合电力系统方向的研究生、做新能源规划的工程师以及所有想用Matlab复现风-光组合分析的读者。代码都是实际跑过的直接复制到你的工程里就能用。1. 项目整体思路拆解两把尺子量同一个问题1.1 为什么非要用概率分布刻画风、光资源风电场和光伏电站的运行数据摆在那里最直观的感受就是“忽高忽低”。今天平均风速6m/s明天可能变成3m/s同一块光伏板上午出力200kW下午云一过来直接掉到50kW。如果只拿平均值去做发电量预测或者容量配置误差会大到让你怀疑人生。概率分布的引入就是为了把这些波动性“收编”成一组可以计算的参数。Weibull分布在风速拟合上几乎是行业标准原因在于它的形状参数k非常灵活——k小的时候分布偏斜能刻画多静风或者多强风的地区k大的时候分布收窄对应风速比较稳定的场址。Beta分布则完全不同它定义在[0,1]这个有界区间上天然适合处理归一化之后的光照强度或者光伏功率既不会出现负值也不会突破物理上限。组合研究的第一步就是先把这两种分布分别拟合好让它们各自成为描述风、光资源的两把可靠标尺。1.2 两个概率模型的数学基础Weibull分布的概率密度函数表达式为f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中v是风速k是形状参数c是尺度参数。实际工程里经常提到的“平均风速”和“有效风速区间”其实都可以从k和c里推导出来。累积分布函数F(v) 1 - exp(-(v/c)^k)在后续做KS检验和场景生成时都会用到。这里特别强调尺度参数c和平均风速不是一回事c通常略大于平均值两者之间还差着一个gamma函数的系数。Beta分布的表达式为f(x) x^(a-1) * (1-x)^(b-1) / B(a, b)其中x落在[0,1]区间a和b是形状参数B(a,b)是Beta函数。这个分布最妙的地方在于当a和b都小于1时它是U形的对应“要么没光、要么满发”的极端天气当a和b都大于1时它是钟形的对应光照相对平稳的场景。做光资源拟合时Beta分布的形态自适应能力非常强。1.3 组合研究的核心目标是什么单看风或者单看光都只能回答“这个资源怎么样”。组合研究要回答的则是“风光在一起到底靠不靠谱”。比如某个地区风大的季节正好光照弱那风光的联合出力可能比单独看任何一种资源都稳定反过来如果风和光总是同强同弱那联合出力的波动反而可能被放大。这个项目通过构建联合概率密度和生成大量随机场景把风速和辐照度映射到功率空间再统计联合出力的概率分布从而量化互补性、评估电力不足风险。这才是整个项目的灵魂所在。2. 从数据到参数Weibull拟合的完整过程2.1 风速数据的清洗与预处理我见过太多人拿到测风塔数据就直接往拟合函数里丢结果拟合出来的曲线惨不忍睹。真实的风速数据里通常混着三类“脏数据”一是缺测值时间戳上就是NaN二是超出物理范围的值比如风速82m/s这种明显是传感器故障的记录三是大量零风速段这些静风数据如果占比过高会严重压低形状参数k使拟合曲线低估高风速段的概率。我的清洗流程是先把风速列单独提取出来用isnan筛选缺测值超过物理范围的值直接置为NaN再统一插值对于静风记录如果占比超过总数据的5%我会考虑用零均值的小噪声去替代避免拟合时出现奇异行为。插值我习惯用线性插值fillmissing简单稳定。清洗之后一定要再画一次时序图确认数据形态合理再进入参数估计环节这一步能帮你省掉后面80%的麻烦。2.2 极大似然估计、矩估计与最小二乘的取舍Weibull参数估计最主流的方法是极大似然估计MLE。Matlab中提供了两个入口一个是直接调用wblfit函数返回形状参数和尺度参数另一个是把数据封装成fitdist对象再取对象的a和b字段。两者的数学本质一样但fitdist在后续配合pdf、cdf等函数时更方便。MLE的迭代公式看起来很简单实际求解时要用数值方法。Matlab内部处理了这部分但默认初值在数据量较小时偶尔会失效。我自己的习惯是先用矩估计算一个初值k ≈ (std/mean)^(-1.086)超参数c mean / gamma(1 1/k)然后把这个初值传给mle函数手动迭代。这样做的原因很实际——MLE的本质是求解非线性方程初值离真实值太远时牛顿迭代法可能发散或者收敛到局部解。矩估计虽然效率不如MLE但胜在稳定用它做初值几乎没出过问题。最小二乘拟合是另一种思路对累积分布函数做双对数变换把Weibull的拟合问题变成线性回归问题。具体做法是对ln(-ln(1-F))关于ln(v)做直线拟合斜率就是k截距和c有关。这个方法在早期文献里很常见现在的工程实践中更多用来做交叉验证。如果MLE拟合出来的结果和最小二乘的结果差距很大那大概率是数据本身有问题而不是方法的问题。2.3 拟合优度检验怎么才算过拟合完参数不能直接就看图说“差不多”。我在项目中至少做两件事第一用KS检验判断拟合分布和经验分布是否来自同一总体第二对比高风速段的尾部拟合效果。KS检验在Matlab里的用法有一个容易踩的坑。kstest函数要求传入CDF函数句柄和参数很多人直接写成kstest(data, CDF, pd)结果是错的正确的写法是kstest(data, CDF, {wblcdf, k, c})。这里的参数cell数组结构必须严格匹配否则Matlab会报错或者给出错误结果。检验输出的p值如果大于0.05说明在95%置信水平下不能拒绝原假设拟合通过。但KS检验也不是万能的它对分布中部的偏差敏感对尾部的偏差反而不敏感。而风速分布尾部恰恰对应着高风速、大出力的区间是发电量评估的关键。所以我每次都会把累积概率在0.9以上的数据点单独画出来看如果尾部拟合偏差超过10%我会考虑改用混合Weibull分布或者截尾Weibull而不是硬着头皮用一个分布去拟合到底。3. 光照强度与Beta分布的建模细节3.1 为什么Beta分布比正态分布更适合光资源很多人习惯性地对任何数据都用正态分布去拟合但光资源数据有一个非常硬性的约束物理上有上下界。辐照度不可能是负值也不可能无限大。如果用正态分布去拟合光强数据拟合结果的左尾会有一部分概率落在负值区间这在物理上是无意义的。Beta分布定义在[0,1]区间只要先做归一化就不会出现越界问题。另一个微妙的地方是光资源的分布形态会随天气类型剧烈变化。晴天时归一化辐照度集中在一个偏高的值附近分布呈尖峰状阴天时分布变得扁平甚至偏向0。Beta分布通过a、b两个参数可以覆盖这两种极端形态这是正态分布做不到的。Beta分布也有两个参数但这两个参数刻画的是“有界区间内的形状”和光资源的物理特征完全对齐。3.2 归一化处理理论最大值还是实测最大值用Beta分布拟合光资源数据之前必须先把原始辐照度数据归一到[0,1]之间。这一步看似简单实际暗藏一个关键决策归一化因子到底取什么。如果取理论最大值比如太阳常数对应的天文辐射上限这样得到的归一化数据几乎很难达到1Beta拟合出来的均值会整体偏低参数形态失真。如果取实测历史最大值拟合精度会明显提升但模型迁移到另一个地区时就不适用了。我在项目中通常采用折中方案先看数据质量如果实测序列足够长至少一年就用实测最大值的某个百分位数比如99.5%分位数做归一化因子既能保留拟合精度又能避免极端峰值对整体归一化的干扰。还有一类特殊情况——夜间零辐照度数据。Beta分布对x0和x1这两端的概率密度是发散的如果在拟合前不处理这些零值会导致参数a趋近于0拟合失败。我的做法是把辐照度低于某阈值的记录单独剔除或者在归一化时给零值加一个非常小的正数偏移。3.3 Beta参数估计与Matlab代码实现Beta分布的参数估计在Matlab里直接调用betafit就能完成函数内部也是用MLE框架。但为了给读者一个更直观的理解我通常先用矩估计算初值。矩估计的公式很简洁a μ * (μ*(1-μ)/σ² - 1) b (1-μ) * (μ*(1-μ)/σ² - 1)其中μ是样本均值σ²是样本方差。如果你发现计算出来的a或者b小于等于0说明这组数据可能不是一个单峰的Beta分布能描述的需要检查数据处理环节。对于大多数实际数据矩估计给出的初值已经足够接近真实值。下面是完整的Matlab拟合代码包含数据归一化和参数估计% 导入辐照度原始数据单位W/m^2 irradiance ...; % 剔除夜间零值保留白天有效数据 valid irradiance 10; % 阈值为10W/m^2 data_day irradiance(valid); % 归一化处理使用99.5%分位数作为归一化因子 norm_factor prctile(data_day, 99.5); x data_day / norm_factor; x(x 0) 1e-4; % 避免零值导致边界发散 % 矩估计初值 mu mean(x); sigma2 var(x); temp mu * (1 - mu) / sigma2 - 1; a_init mu * temp; b_init (1 - mu) * temp; % MLE拟合 phat betafit(x, [], [a_init, b_init]); a phat(1); b phat(2); % 拟合效果画图 figure; histogram(x, 30, Normalization, pdf, FaceColor, [0.8 0.8 0.8]); hold on; xx linspace(0, 1, 200); plot(xx, betapdf(xx, a, b), r-, LineWidth, 2); xlabel(归一化辐照度); ylabel(概率密度); legend(实测直方图, Beta拟合曲线);注意betafit的第二个参数是初始值第三个参数是优化选项。如果直接不给初始值Matlab会内部自动选取初值但实测下来给一个合理的初值能让收敛更快、结果更稳定。3.4 从辐照度到光伏出力的转换模型Beta分布拟合出来的是辐照度的统计特征但最终评价光伏出力还需要一个转换环节。工程中常用的是两参数分段模型P Pmax * (R - Rc) / (Rmax - Rc)当 R Rc P 0当 R ≤ Rc其中R是辐照度Rc是启动阈值一般取120-200W/m²Rmax是参考辐照度Pmax是装机容量。这个模型比线性模型多了一个“死区”能更真实地反映光伏组件在低辐照度下近乎不出力的物理特性。在组合仿真时先用betarnd生成归一化辐照度样本乘以归一化因子还原成辐照度再通过这个转换模型得到功率序列最终统计的光伏出力分布就会非常贴近实际。4. 风光组合仿真的Matlab代码实现4.1 蒙特卡洛场景生成的基本流程组合研究最核心的代码实现是蒙特卡洛模拟。整体思路是这样的先用拟合好的Weibull参数生成风速样本用Beta参数生成归一化辐照度样本再分别转换成功率最后把两种功率叠加起来形成联合出力序列。这里有一个重要的前提假设要说明我在基础版实现里假定风速和辐照度在小时尺度上是相互独立的。这个假设并不完全成立尤其在某些特定气候区风和光可能存在弱相关。但作为第一版组合研究独立假设已经能给出很多有价值的结论。想考虑相关性的话可以引入Copula函数做联合抽样相当于给两个边缘分布加一个相关结构这个在后续进阶版本里再详细展开。由于蒙特卡洛需要大量样本我建议用向量化操作而不是循环。生成10万个样本在Matlab里就是一条命令的事运行时间只需要零点几秒% 设置随机种子保证实验可复现 rng(42); % 生成风速场景样本 n_samples 100000; wind_speed wblrnd(k, c, n_samples, 1); % 生成归一化辐照度场景样本 norm_irradiance betarnd(a, b, n_samples, 1); % 反归一化得到真实辐照度 irradiance norm_irradiance * norm_factor;wblrnd的第一个参数是形状参数k第二个是尺度参数c和wblfit的输出顺序一致这个细节写代码时很容易搞反。betarnd的两个参数就是Beta分布的形状参数a和b。4.2 风功率和光伏功率的转换计算风速样本生成之后要借助风功率曲线把风速映射到风机出力。风电功率曲线一般是S形切入风速以下出力为0额定风速以上出力被限制在额定功率。用一个分段函数来表示% 风功率曲线参数 v_in 3; % 切入风速单位m/s v_out 25; % 切出风速 v_rated 12; % 额定风速 P_w_rated 1.0; % 标幺值取额定功率为1 % 计算风功率 P_wind zeros(n_samples, 1); idx1 wind_speed v_in wind_speed v_rated; idx2 wind_speed v_rated wind_speed v_out; P_wind(idx1) (wind_speed(idx1)^3) / (v_rated^3 - v_in^3) - (v_in^3) / (v_rated^3 - v_in^3); P_wind(idx2) P_w_rated;这里的功率曲线模型是简化版实际项目里应该用风机厂商提供的实测功率曲线表。但简化模型对于算法验证和分布组合研究已经足够了因为核心目的是展示随机变量通过非线性变换后的分布变化规律。光伏功率转换就用前面提到的分段模型% 光伏转换模型 R_c 150; % 启动阈值W/m^2 R_max norm_factor; % 参考辐照度 P_pv zeros(n_samples, 1); valid_idx irradiance R_c; P_pv(valid_idx) (irradiance(valid_idx) - R_c) / (R_max - R_c); P_pv min(P_pv, 1); % 限制最大出力不超过额定值两组功率都算出来之后定义联合出力P_total P_wind P_pv; % 风光联合出力标幺值这里的P_total就是每次随机场景对应的联合出力。联合出力的上限是2风、光各1下限是0。有了这10万个联合出力样本后面所有的统计分析都顺手了。4.3 联合概率密度与互补性分析的可视化联合出力样本生成之后最直观的展示方式是二维直方图。用Matlab的histogram2函数可以把风速和辐照度的联合分布画成一个三维地形图颜色深浅代表概率密度的高低。如果数据集中在对角线附近说明风、光存在同涨同落的相关性如果数据分布呈现“反对角”形态则说明互补性强——风速大的时候光照小风速小的时候光照大。另一个更实用的输出是联合出力的累积分布函数。用ecdf函数可以直接画出经验CDF曲线从这条曲线上可以读出工程上非常有价值的指标比如P_total小于某个阈值的概率。假设电力系统要求联合出力必须大于0.4标幺值如果从CDF上读出这个概率大于某个可接受值那就说明需要配置储能或者扩大装机容量。这就是组合研究从理论走向工程决策的落点。互补性指标我一般用Pearson相关系数和Spearman秩相关系数同时计算corr_pearson corr(wind_speed, irradiance, Type, Pearson); corr_spearman corr(wind_speed, irradiance, Type, Spearman);Pearson相关系数描述线性相关Spearman描述单调相关性。如果两者符号一致都能反映互补趋势如果符号相反说明两个变量之间存在非线性关系这时候单纯用相关系数做判断就不够了需要进一步借助Copula模型的尾部相关系数来识别极端场景下的互补性。5. 踩坑记录与排查技巧实录5.1 参数初值不当导致的拟合失败我在做Beta分布拟合时遇到的问题最典型。有一组光照数据来自某高海拔地区晴天极多归一化辐照度集中在0.8-0.95之间方差非常小。直接用betafit拟合结果返回了离谱的参数a50b2.5拟合曲线尖得几乎像一根线。后来排查发现是矩估计初值中的σ²出现了下溢导致a_init计算异常。解决办法是对原始数据先做一次标准化检查如果方差过小给矩估计公式中的方差项设置一个下限或者直接改用极大似然估计的默认初值。这个坑给我的教训是不要想当然地认为可调函数总能自适应一切数据理解初值的作用机制才能在异常结果出现时快速定位问题。5.2 零值数据和极端值对拟合的干扰风光数据的共同特点是大量零值。风速数据中的静风段表现在时域上就是一段平直为零的曲线光资源数据中的夜间时段则是另一个零值来源。如果不加处理直接喂给拟合函数Beta分布会因为x0处的概率密度发散而完全崩溃Weibull分布中的零值则会让形状参数k被拉到极低值。处理零值我总结了三条经验第一风资源分析中保留静风段但拟合前先计算静风占比占比过高时要考虑两参数Weibull是否还适用第二光资源分析中直接剔除夜间零值只拟合白天的有效数据第三对归一化之后出现的零值统一加上一个1e-4到1e-3量级的小偏移避免数值问题。这个小偏移看似不起眼却能拯救整个拟合过程。5.3 随机种子的设置与结果复现问题蒙特卡洛仿真如果每次运行结果都不一样实验就失去了可复现性。这个问题在论文审稿和工程评审时尤其致命。第一次跑出功率不足概率是8.5%第二次跑变成9.1%评审直接质疑数据的可靠性。解决方法是严格使用rng函数设置随机种子而且是每次生成随机样本之前都设置一次。另外要特别注意wblrnd和betarnd如果每次调用前都重置种子确实能保证序列完全一致但同时意味着统计上每次都是同一组伪随机数在某些极端情况下可能引入偏差。在实践上我通常生成一次大规模样本后保存到.mat文件后续所有分析都从文件中读取样本既保证可复现又方便不同参数方案的横向对比。5.4 大样本仿真时的内存与速度优化10万样本对Matlab来说是小菜一碟但如果你做的是8760小时的全年场景再加上几十年的历史数据回测内存占用就会变成一个大问题。一个常见的错误是频繁在循环里调用wblrnd和betarnd每次生成一个样本再拼接。正确的做法是提前一次性生成所有样本用向量化计算代替循环。如果样本量实在太大考虑分批生成并增量累加统计量而不是把全部样本都存在内存里。我在处理百万级样本量时还会利用parfor做并行计算配合数据分块策略运行时间可以减少60%以上。我个人做风光组合研究走下来最大的体会是拟合分布只是整个项目的地基真正出成果的地方在组合分析和场景生成。Weibull和Beta的公式本身并不复杂Matlab里对应的拟合函数也都是现成的但如果忽略了数据清洗、参数初值、零值处理这些细节后面每一步都会带着偏差走。建议读者拿到项目后先花一半时间把两份真实数据彻底清理干净再把参数估计和组合仿真的代码跑通最后用一组已知结果的数据做验证。这样下来你的风光联合出力分析模型才算真正立得住。
