搞新能源的人应该都听过这句话风电看分布光伏看天气。这里的“分布”不是玄学而是决定电场出力特性的概率密度函数。国内做风资源评估、光伏电站设计、储能容量配置的工程师十有八九都绕不开两个数学模型——风速的Weibull分布和光伏出力的Beta分布。这次我把两者在Matlab里完整实现了一遍包括参数估计、拟合优度检验、联合场景生成以及相关性控制把整个过程和踩过的坑一条一条理清楚。这篇内容适合正在做风-光互补研究、概率潮流、储能容量配置或者毕业设计需要做新能源出力建模的读者尤其是想用Matlab快速落地又不想只会调函数、不理解背后逻辑的人。先说结论这套组合建模方法解决的核心问题是把“风”和“光”两种随机性用数学模型描述出来再通过蒙特卡洛采样生成工程可用的联合场景。它不是简单的两个分布相加而是要解决“风大光少、光大风小”这类天然互补特性如何在数值仿真中被正确还原。下面从选型逻辑开始逐步拆到代码实现和工程落地的每一个细节。1. 先想清楚Weibull和Beta为什么能凑成一对黄金搭档1.1 风速的物理特性决定了Weibull的统治地位风速数据的统计特性有几个硬约束风速恒为非负值、分布具有明显的右偏特征、低风速时段占比通常较高且不同地域的分布形态差异极大。Weibull分布恰好具备这些特性它的形状参数k控制曲线的偏斜程度当k接近2时分布近似Rayleigh分布当k在2到3之间时曲线峰值明显且右尾较厚这个形态和绝大多数风电场实测风速直方图是非常吻合的。尺度参数λ则直接和平均风速挂钩平均风速约等于λ乘以Γ(11/k)所以从概率密度曲线就能推断出该地点的风资源整体水平。从行业规范角度看IEC 61400-12等风资源评估相关标准中Weibull分布是描述风速长期分布的主流推荐模型。实际项目中我也对比过其他候选分布比如对数正态分布和Gamma分布这两种分布虽然也能拟合偏态数据但面对风速数据中常见的“长右尾近零频次较低”的组合形态时往往在尾部拟合精度上吃亏。在Matlab中做拟合时Weibull的参数估计过程也很稳定极少出现数值震荡这一点在后文实测代码中会看到。1.2 Beta分布恰好落在光伏出力的定义域上光伏出力和风速不同物理上它天然被限制在一个区间内。如果我们把某时刻的实际出力除以额定容量得到一个0到1之间的标幺值那么这个值在夜间是0在阴天是低值在晴天高值附近频繁波动。历史出力标幺值的直方图通常表现为单峰、偏斜、贴近边界这三种特性之一Beta分布正好定义在[0,1]区间上且通过调节两个形状参数α和β可以拟合出左侧偏斜、右侧偏斜、近似对称甚至U型曲线灵活性很强。Beta分布的均值公式是α/(αβ)方差公式是αβ/((αβ)²(αβ1))。这两个公式在实际工程中非常有用因为当我拿到实测出力数据后可以先计算样本均值和方差代入公式反推出α和β的初值再用最大似然估计做精修。这种两步走的方式大大降低了数值迭代不收敛的概率。另外需要说明的是严格意义上辐照度本身也常用Beta分布描述工程中之所以直接用光伏出力标幺值是因为出力与辐照度之间经过温度、效率等因素转换后仍然保持0到1的区间特征Beta分布对出力数据的拟合效果在多数情况下足够好。1.3 所谓“组合研究”到底在组合什么很多人看到“组合”两个字第一反应是把Weibull概率密度和Beta概率密度画在同一张图上或者把两个分布函数乘一乘。这些都不对。真正的组合研究核心是构建风-光联合出力场景。风电场和光伏电站的出力在时间尺度上存在复杂的相关关系。晴天的中午光伏出力接近满发但风速往往偏低夜间光伏出力为零但此时风速可能达到日峰值。如果我们要做储能容量配置或者概率潮流分析就必须知道“风光同时段的典型出力组合”长什么样而不是只知道各自边缘分布。组合的方法就是先对风速和光伏出力分别建立精确的边缘分布再通过相关性建模把两者关联起来最后用采样生成联合场景。这里面涉及到概率积分变换、Copula函数、蒙特卡洛采样、场景削减等一整套工具也是整篇内容的重点。Matlab在这条技术路线上的支持非常完善统计工具箱中的wblfit、betafit、copulafit、copularnd都是为这类需求设计的。2. Matlab环境准备与参数估计算法拆解2.1 工具箱版本与数据准备规范做这套实现需要Matlab的Statistics and Machine Learning Toolbox也就是统计与机器学习工具箱。新版本R2023b、R2024a自然没问题老版本R2016a之后的统计工具箱里wblfit、betafit、copulafit这些核心函数都存在接口也基本没变过所以即使是十年前的老环境也能跑通。需要提醒一句脚本里如果写了中文注释务必把Matlab编辑器默认编码设置成UTF-8否则在部分中文本地化版本里注释和字符串会出现乱码耽误排查时间。数据准备是整个流程里最容易被低估的环节。风速数据建议取测风塔10米或轮毂高度处的10分钟平均风速序列周期至少覆盖一年才能包含四季特征。光伏出力数据则按电站实际并网功率曲线采集时间分辨率同样建议10分钟或1小时。拿到原始数据后需要先剔除异常值和停机时段风速数据中零值如果来自测风设备故障应该标记为缺失而不是直接参与拟合光伏数据的处理更讲究夜间出力为0是正常物理现象但如果直接把大量0值塞进Beta分布的最大似然估计会得到荒谬的参数后文我会单独讲这个边界问题。2.2 最大似然估计wblfit和betafit背后在干什么Matlab的wblfit函数默认使用最大似然估计MLE来求解Weibull分布的两个参数。Weibull的似然函数涉及形状参数k的隐式方程没有解析解函数内部通过数值迭代方式求解。实际使用时wblfit返回两个输出第一个是参数估计向量parmhat元素顺序依次是形状参数k和尺度参数λ第二个是参数置信区间parmci这个区间在做工程报告时非常有用比如k的95%置信区间能帮我判断拟合结果是否稳定。如果你不想用内置函数也可以自己写目标函数交给fminsearch优化但实际测试下来wblfit的迭代速度和稳定性都更好没必要重复造轮子。betafit函数的逻辑类似但有一个重要的实现差异Beta分布的定义域是[0,1]闭区间而最大似然估计在数据恰好等于0或1时对数似然函数会趋向无穷直接导致迭代崩溃。betafit内部会对边界数据做轻微调整但现实中如果光伏出力数据包含大量精确满发或精确零出力点仍然会得到偏差较大的结果。所以我在实际代码里会先做数据清洗把接近0和接近1的数据做截断处理再交给betafit。2.3 拟合优度检验画图别只看个热闹拟合完成之后直接画一条概率密度曲线叠在直方图上确实最直观但图上的“差不多”并不能替代定量检验。我自己做项目时的标准动作是三件套概率图、KS检验、卡方检验。概率图用Matlab的wblplot它会按Weibull分布的理论分位数和样本分位数画散点如果数据点大致落在一条直线上说明拟合良好。对于Beta分布没有专门的plot函数可以用probplot加自定义分布句柄实现或者用经验CDF和理论CDF画对比散点。KS检验用kstest直接把样本和指定分布对象传进去即可输出p值大于显著性水平比如0.05表示不能拒绝“样本来自该分布”的原假设。卡方检验用chi2gof需要提前把数据分箱然后比较频数。这里有个经典的坑当样本量很大比如8760个小时点KS检验和卡方检验都会非常敏感轻微的尾部偏差就会导致p值极小但这并不代表拟合失败。我的处理方式是同时看拟合误差指标RMSE和概率图综合判断而不是机械地看p值。3. 核心代码实现风速Weibull与辐照Beta的完整拟合流程3.1 风速Weibull拟合三行核心代码加一张诊断图以下代码是我在一个实际风电场项目里的标准写法数据来自某测风塔一年的10分钟平均风速序列共52560个点。先剔除掉小于0的异常数据和明显的设备故障时段然后调用wblfit%% 风电风速 Weibull 分布拟合 load wind_data.mat; % wind_speed: 52560×1, 单位 m/s speed wind_speed(wind_speed 0); % 剔除静风和异常 [parmhat, parmci] wblfit(speed); k parmhat(1); lambda parmhat(2); figure; histogram(speed, 40, Normalization, pdf, FaceAlpha, 0.35); hold on; x_axis 0:0.1:max(speed); y_fit wblpdf(x_axis, k, lambda); plot(x_axis, y_fit, r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); title(sprintf(风速 Weibull 拟合结果: k%.2f, lambda%.2f, k, lambda)); grid on;这段代码运行后我实测的一组数据得到的拟合结果是k2.13λ7.45。k值略大于2说明风速分布比标准的Rayleigh分布更尖锐峰值风速大约在5到6米每秒附近右尾存在较大的高风速概率。λ7.45意味着该测风塔的年平均风速大约在6.6米每秒左右换算结果和实测平均风速吻合度很高。这里看直方图和拟合曲线时重点看峰值位置和右尾衰减趋势是否一致很多人只盯着峰值忽略了右尾但在风资源评估中高风速段的拟合精度直接影响发电量估算和大风工况判断。3.2 光伏出力Beta拟合边界处理才是真正的坑光伏出力数据的Beta拟合比风速要麻烦得多因为数据边界问题直接决定拟合成败。我的一个实际光伏电站数据是额定容量1兆瓦按小时采集全年8760个点。首先把出力折算成标幺值也就是除以电站额定容量得到0到1之间的序列。然后面临的第一个问题是夜间零值占比接近一半直接把所有零值丢给betafit会出现严重的数值异常第二个问题是晴天的中午出力可能精确等于1这也会导致MLE迭代崩溃。我的处理方式是先剔除边界数据再进行拟合。严格来说这会损失一部分边界概率信息但工程上可接受。如果想要更精细可以采用零膨胀Beta分布或截断Beta分布但在常规的容量配置研究中简单的边界剔除配合良好的中段拟合已经足够。代码如下%% 光伏出力 Beta 分布拟合 load pv_data.mat; % pv_output: 8760×1, 单位 kW rated_power 1000; % 额定容量 kW pv_pu pv_output / rated_power; % 标幺化 % 边界处理剔除夜间零出力与满发时段 pv_clean pv_pu(pv_pu 0.01 pv_pu 0.99); phat betafit(pv_clean); alpha phat(1); beta phat(2); figure; histogram(pv_clean, 40, Normalization, pdf, FaceAlpha, 0.35); hold on; x_axis 0:0.005:1; y_fit betapdf(x_axis, alpha, beta); plot(x_axis, y_fit, b-, LineWidth, 2); xlabel(光伏出力标幺值); ylabel(概率密度); title(sprintf(光伏出力 Beta 拟合结果: alpha%.2f, beta%.2f, alpha, beta)); grid on;这组数据的拟合结果是α0.86β1.94。α小于β说明分布形态右偏也就是低出力时段的概率密度更大这个特点符合实际光伏电站的特性夜间零出力被剔除后剩下的大量时段集中在低辐照区间晴天高出力时段的概率密度相对较低但尾部较长。如果某个电站处于光照资源极好的地区α可能大于β分布会向左偏高出力时段占比增大这也是Beta分布灵活性的体现。3.3 拟合参数怎么用从概率密度到工程指标有了参数以后不能只满足于画图好看要把它转化成工程上可用的指标。第一类指标是期望出力Weibull分布的期望是λ乘以Γ(11/k)Beta分布的期望是α除以αβ。以我们前面的拟合结果为例风速期望约6.6米每秒光伏出力的条件期望约0.31标幺值。第二类指标是置信区间风电出力P90和P10对应的风速值可以通过wblinv和betainv函数直接计算这在风电场可发电量评估中经常用到。第三类指标是出力波动区间如果储能容量配置需要覆盖95%的场景那就需要把风光联合出力的5%分位数和95%分位数算出来这就要用到下一章的联合采样方法了。参数本身还有一个自检功能。当k值小于1.5时说明当地风速分布极度偏斜小风速频次极高这种场址往往不具备开发价值当α和β都小于1时Beta分布会变成U型曲线意味着光伏出力大部分时间要么接近0要么接近满发这对储能系统来说是非常严苛的运行工况设计容量要留足余量。这些判断用输出了几秒钟就能得到比反复试凑场景要高效得多。4. 风光联合场景生成从独立采样到相关性控制4.1 独立蒙特卡洛采样先看随机数生成说了什么建立完边缘分布后最直接的做法是假设风速和光伏出力相互独立用wblrnd和betarnd分别生成大量随机样本再拼在一起。代码如下N 10000; wind_sample_ind wblrnd(k, lambda, N, 1); pv_sample_ind betarnd(alpha, beta, N, 1); figure; scatter(wind_sample_ind, pv_sample_ind, 5, k, filled, Alpha, 0.2); xlabel(风速 (m/s)); ylabel(光伏出力标幺值); title(独立采样风-光联合散点图);这个散点图跑出来是一个近似椭圆形的云团风速从0到30均匀分布光伏出力从0到1基本随机搭配风速大小和光伏出力高低之间没有明显的规律。独立假设在数学上很简单但在实际风-光资源评估里是有问题的。同一地区的风速和辐照度往往由相似的气象系统驱动晴天可能同时伴随大风也可能大风出现在阴雨天这种相关性如果用独立采样来模拟会高估“大风高光伏”这类极端场景的出现概率也会低估“无风无光”的极端枯竭场景。做储能容量配置时如果输入场景集本身就失真得到的容量结论自然不可靠。4.2 用Copula管住相关性把两个边缘分布粘起来Copula方法的核心想法是先把每个变量的原始观测值通过概率积分变换转换成0到1之间的均匀分布变量再用一个连接函数描述这些均匀变量之间的相关性最后再逆变换回原始物理量。Matlab中实现非常直接%% 数据预处理确保两个序列按时间戳对齐 valid wind_speed 0 pv_pu 0.01 pv_pu 0.99; s_w wind_speed(valid); s_pv pv_pu(valid); %% 经验CDF变换 u1 tiedrank(s_w) / (length(s_w) 1); u2 tiedrank(s_pv) / (length(s_pv) 1); %% 拟合Gaussian Copula rho copulafit(Gaussian, [u1, u2]); %% 生成带相关性的联合场景 U copularnd(Gaussian, rho, N); wind_sample wblinv(U(:,1), k, lambda); pv_sample betainv(U(:,2), alpha, beta); figure; scatter(wind_sample, pv_sample, 5, r, filled, Alpha, 0.2); xlabel(风速 (m/s)); ylabel(光伏出力标幺值); title(sprintf(Gaussian Copula联合采样秩相关系数%.2f, rho(1,2)));我这里用的tiedrank是秩变换分母加1是为了避免变换后的值取到0或1防止后续copulafit在边界处出现数值异常。这个细节是我早期踩过坑后总结出来的直接用ecdf函数得到的经验CDF在样本最大值处会等于1而Copula模型对边界值非常敏感。在真实项目中实测数据的秩相关系数通常在-0.3到0.2之间波动负相关说明一个资源上升另一个资源下降风-光互补特性明显如果r接近0.5以上说明该区域风光资源同涨同跌这对储能容量设计影响极大需要重点评估。如果高斯Copula拟合效果不理想尤其是风速和光伏出力存在厚尾相关性时可以考虑改用t Copula只需要把copulafit的第一个参数改成t它会额外估计一个自由度参数ν输出变成[rho_t, nu]。我在实际项目中测试过当样本量少于2000时t Copula的参数估计方差较大高斯Copula反而更稳定样本量充足时t Copula的尾部拟合更精准。选择依据不是越复杂越好而是看数据特征和后续用途。4.3 场景削减从10000个样本到50个典型日Matlab生成的10000个联合场景直接传给后续的时序仿真程序计算代价太大。储能容量优化或者概率潮流分析往往需要在优化迭代中反复调用场景集必须做场景削减。最常用的做法是K-means聚类把10000个样本聚成50类每一类的聚类中心代表一个典型场景该类样本数量占总样本数量的比例作为该场景的概率。num_scenarios 50; [idx, C] kmeans([wind_sample, pv_sample], num_scenarios, Replicates, 5); scenario_prob accumarray(idx, 1, [num_scenarios, 1]) / N; figure; scatter(C(:,1), C(:,2), scenario_prob * 1000, filled); xlabel(风速 (m/s)); ylabel(光伏出力标幺值); title(削减后的50个典型风光联合场景);这段代码用scatter的第三个参数控制点的大小点的面积正比于场景概率直观展示哪些组合场景出现频率高。设置Replicates参数为5表示用5组不同的初始聚类中心分别计算取全局最优避免落入局部极值。聚类数选择没有绝对标准通常的做法是从20到100每隔10个做一次肘部法则检验看总离差平方和的下降趋势在哪个点开始变缓。50个场景是我在储能容量配置项目中使用较多的一组配置既能控制优化模型求解时间又能保留风光出力的典型波动特征。5. 常见问题与排查技巧实录5.1 数据预处理里的三条教训第一条教训是风速数据中的零值不要一刀切删除。零值可能来自静风时段也可能来自测风设备冻结或故障。如果静风在气候特征中确实存在那么零值本身是有物理意义的直接删除会让k值偏大、λ值偏小如果是设备故障导致的假零值不删除则会严重拉低λ。我的区分方法是结合温度、湿度等辅助气象数据做合理性校验批量查看零值出现的时段分布。第二条教训是光伏出力数据的满发点处理。晴天中午出力精确等于1的点如果占比超过5%简单的剔除会让尾部概率失真更合理的做法是采用混合Beta模型也就是在Beta分布基础上叠加两个离散点分别对应0和1的概率。第三条教训是和联合建模强相关的风速序列和光伏序列在做秩变换之前必须严格按时间戳对齐不能一个序列剔除了停机数据、另一个序列剔除了夜间数据否则两个向量长度对不上tiedrank和copulafit直接报维度错误就算勉强对齐了秩相关的含义也变了。5.2 参数估计不稳定的排查思路betafit报错最常见的原因是输入向量里有0或1。我先检查pv_clean的范围一旦发现极值就按前文方法剔除。wblfit不收敛的情况比较少见但我在处理样本量极小比如少于30个点或者风速数据全部集中在某个狭窄区间时遇到过参数估计方差过大的问题这时parmci的区间会非常宽。解决方法是补充数据量或者采用矩估计给出初值再交给MLE精修。另一个容易被忽略的问题是数据的时间分辨率如果输入的是逐小时平均风速而你的仿真步长是10分钟直接把逐小时拟合参数应用到10分钟尺度会产生偏差。风速的采样间隔越小k值通常会略有减小尾部更厚所以在报告参数时需要标注清楚数据的时间分辨率。5.3 联合采样合理性检查清单生成联合场景之后不要直接拿去用先做三个检查。第一个检查是边缘分布还原度把采样的风速和光伏出力各自画直方图和原始的Weibull、Beta理论密度曲线做对比偏差过大说明逆变换环节写错了。第二个检查是相关性保持度计算联合样本的Kendall或Spearman秩相关系数和原始数据的秩相关系数对比差距应该在0.02以内。第三个检查是物理合理性风速样本不应该出现负值光伏出力样本应该在0到1范围内联合散点图中不应该出现在极端不合理区域的高密度簇。我见过不少初学者的代码跑完copularnd后直接忘记逆变换结果风速样本变成了一堆0到1之间的小数画出来的散点图完全看不出物理含义这类低级错误在自检清单下无处遁形。5.4 问题排查速查表下面这个表格是我在实际项目答疑时经常用到的速查列表整理给读者参考。问题现象可能原因排查与解决方案betafit报错或输出异常输入数据含0或1剔除边界数据或改用零膨胀Beta模型wblfit的parmci区间过宽样本量太小或风速分布过窄增加数据覆盖时长检查异常值剔除是否过度联合采样散点图无物理意义忘记对Copula样本做逆变换检查是否调用wblinv和betainv转换U矩阵联合样本秩相关系数严重偏离原始值两个序列未按时间戳对齐先做valid索引统一筛选再进行秩变换场景削减后极端场景消失聚类数过少或K-means未收敛增大聚类数设置Replicates并比较多次优化结果中文注释在Matlab中乱码编码设置不统一统一使用UTF-8保存脚本检查系统区域设置我在实际项目里比较偏爱把整套流程封装成两个独立函数一个处理风速Weibull拟合一个处理光伏Beta拟合输入历史序列后直接返回参数、置信区间和诊断图。第一次做的时候我也直接调内置函数一把梭直到有一次评审专家问“拟合优度检验结果怎么样”才发现自己拿不出定量证据只能指着图说大概可以。后来补上了概率图和KS检验的输出再遇到质疑时直接把p值和RMSE一摆问题就解决了。这套流程的另一个延伸方向是把它推广到多个风电场和光伏电站的联合出力建模只需要把Copula的维度从二维扩展到多维其余流程完全一致。如果你正在做新能源出力不确定性建模建议先把本文的代码用自己的数据跑一遍再去调整分布族和Copula类型这样踩坑的成本会小很多。