风光出力组合建模：Weibull与Beta分布联合分析及Matlab实现

风电的Weibull分布及光电的Beta分布组合研究Matlab代码实现做新能源规划或者微电网仿真的朋友肯定绕不开一个问题风电场和光伏电站的出力到底怎么用数学描述搞并网仿真、容量配置、储能调度都得先把电源侧的随机性刻画清楚。这个领域里最常用的两把尺子就是——风电出力用Weibull分布光伏出力用Beta分布。最近我把两者组合起来做了一套完整的Matlab代码从参数拟合到联合概率分析一步到位这篇就把思路和踩过的坑一起梳理出来。先说清楚这套东西解决什么问题单看风电Weibull分布描述风速的统计规律单看光伏Beta分布描述光照强度或归一化出力的概率特征。但实际工程里风光往往是同一个园区、同一个配网节点你需要知道风大时光照如何光伏出力高时风电状态怎样这就必须把两个分布组合起来研究。代码实现的核心包括三块分布参数估计、随机样本生成、联合概率分析框架。适合做新能源出力建模、微电网可靠性分析、储能容量优化的朋友参考。1. 为什么要把风电和光电放在同一个模型里算1.1 风光互补的前提是能联合描述随机性很多人都知道风光互补这个说法但真正做仿真时才发现互补性不是一句口号你得用一个数学框架把两边的随机性同时装进去。单独建模风电、单独建模光伏这两件事都不难。难的是回答这样的工程问题当风电出力处于低水平时光伏出力的概率分布是什么样的两者同时处于高出力区间的概率是多少如果要配置一条联络线或一组储能需要知道风光联合出力的最坏情况这根本无法从两个边缘分布独立推算出来。我记得第一次做微电网可靠性评估时想当然地把风电和光伏当作独立变量分别抽样然后相加。结果算出来的系统缺电概率明显偏乐观后来才发现问题出在忽略了两个电源在时间尺度上的关联性——比如同一个地区的风和光往往受同一天气系统支配夜间风大但无光白天光照强但风速可能下降。这种时间上的此消彼长单靠两个独立的边缘分布根本表达不出来。1.2 Weibull和Beta各自承担的角色这个组合模型里Weibull和Beta各有各的阵地。Weibull分布描述风速的概率密度它有两个参数形状参数k和尺度参数c。k决定了风速分布的偏斜程度k越小曲线越右偏说明小风速出现频率更高c决定了平均风速的量级。实际风电场的风速数据用Weibull分布拟合效果通常非常好这有气象统计的理论支撑也有大量工程验证。Beta分布则是描述区间[0,1]内连续随机变量的利器。光伏出力归一化后实际出力除以额定容量天然落在0到1之间Beta分布的两个形状参数α和β非常灵活——α和β都大于1时曲线呈单峰α小于1、β大于1时曲线呈单调递减反过来则递增。这正好能刻画光伏出力在不同天气条件下的分布形态晴天出力集中在高位阴天则集中在低位甚至接近零。不过要提醒的是Beta分布的参数初值如果给得不好拟合时很容易发散这点后面专门讲。2. Weibull分布建模风电从风速实测到出力曲线2.1 分布参数的物理含义别搞反了先上一段最基础的Matlab代码用极大似然估计拟合风速数据的Weibull参数% 风速实测数据单位m/s可以是某风电场一年的小时级数据 wind_speed [2.3, 3.1, 4.5, 5.2, 6.8, 7.1, 8.4, 9.2, 7.8, 6.3, ...]; % 方法一用wblfit直接拟合Statistics Toolbox [parmhat, parmci] wblfit(wind_speed); k parmhat(1); % 形状参数 c parmhat(2); % 尺度参数 % 方法二手动实现极大似然估计便于理解原理 % Weibull的对数似然函数 negloglik (params) -sum(log(wblpdf(wind_speed, params(1), params(2)))); params0 [2, mean(wind_speed)/gamma(1.5)]; % 初值估计 options optimset(Display, off); [params_ml] fminsearch(negloglik, params0, options); k_ml params_ml(1); c_ml params_ml(2);先说结论直接用wblfit最省事底层就是极大似然估计数值稳定性有保证。手动实现那套主要是为了让大家看懂原理实际项目里我建议直接用内置函数。但这里有个关键点拟合风速分布和拟合风电出力分布是两回事。很多初学者把风速的Weibull分布直接当成出力分布这是错的。风速通过功率曲线转换为出力后分布形态会发生显著变化——切入风速以下出力为零额定风速以上出力被截断为额定值这两处会产生概率堆积。所以如果你要做的是出力层面的组合分析建议直接对归一化出力数据建模或者用风速分布叠加功率曲线转换。2.2 形状参数为什么重要三个典型场景形状参数k的现实意义常被低估。k2时就是Rayleigh分布很多教材默认用这个值但实际风电场的数据拟合出来k通常落在1.5到3之间。我做过几个不同风场的数据对比风场类型典型k值对出力的影响季风影响明显风速变化剧烈1.5~1.8中等风速占比高出力波动大沿海稳定风场2.0~2.4出力曲线较平滑高原或山口地形2.5~3.0风速集中在中高区间利用小时数高这里有一个实操中的坑如果直接用wblfit拟合有时k的估计值会异常偏大比如超过4这往往是因为数据里有大量重复值或极端值。我建议拟合前先做数据清洗剔除异常记录例如风速超过风机的切出风速通常是25m/s的记录。% 数据清洗示例 valid_idx wind_speed 0 wind_speed 25; wind_speed_clean wind_speed(valid_idx);2.3 从风速到功率分段函数的处理细节风速的Weibull分布拟合完成后下一步是转换为出力。标准的风功率曲线是分段函数% 风机参数 v_cut_in 3; % 切入风速 m/s v_rated 12; % 额定风速 m/s v_cut_out 25; % 切出风速 m/s P_rated 1500; % 额定功率 kW % 功率曲线转换为归一化出力 function p_norm wind_power_curve(v, v_ci, v_r, v_co) p_norm zeros(size(v)); % 区域1低于切入风速出力为0 % 区域2切入风速到额定风速之间按三次方关系近似 idx (v v_ci) (v v_r); p_norm(idx) ((v(idx) - v_ci) ./ (v_r - v_ci)).^3; % 区域3额定风速到切出风速之间出力为1 idx (v v_r) (v v_co); p_norm(idx) 1; % 区域4超过切出风速停机出力为0 end需要注意的是这个三次方关系的近似处理在不同文献里有不同版本。有的用平方关系有的用线性关系还有的用更精细的分段函数。对于组合研究来说三次方近似已经足够但如果你要精确模拟某型号风机最好从厂家数据手册中获取实际功率曲线。3. Beta分布建模光伏光照随机性的数学刻画3.1 为什么偏偏是Beta分布而不是正态分布光伏出力的核心驱动因素是光照强度而光照强度受云层遮挡、大气透明度、昼夜变化等因素影响。用归一化出力P实际出力/额定容量作为随机变量时P的取值范围是[0,1]Beta分布恰恰定义在[0,1]区间上这是它最根本的优势。正态分布虽然用起来方便但它定义在整个实数轴上如果硬要用正态分布拟合归一化光伏出力会在0和1附近产生不合理的概率密度。另一方面Beta分布的形状极其灵活——α1且β1时呈U型适合描述要么大晴天要么阴天的极端天气模式α1且β1时呈钟形适合描述过渡季节的温和天气。这种灵活性是其他分布很难替代的。给一个直观的代码示例用Beta分布拟合归一化光伏出力数据% 归一化光伏出力数据取值范围[0,1] pv_output [0.1, 0.15, 0.2, 0.35, 0.5, 0.65, 0.72, 0.8, 0.85, 0.9, 0.88, ...]; % Beta分布参数估计 % 方法一内置函数betafit [phat, pci] betafit(pv_output); alpha_hat phat(1); beta_hat phat(2); % 方法二极大似然估计数据量小时更稳定 negloglik_beta (params) -sum(log(betapdf(pv_output, params(1), params(2)))); params0_beta [2, 2]; % 初值 [params_beta] fminsearch(negloglik_beta, params0_beta, options); alpha_ml params_beta(1); beta_ml params_beta(2);3.2 数据归一化时最容易出错的边界问题Beta分布定义在开区间(0,1)上但实测光伏出力数据经常出现恰好等于0夜间或阴天或恰好等于1满功率输出的情况。这些边界值在Beta分布下概率密度为0直接拟合会出问题。常见的处理办法是把0和1附近的极值做微小偏移% 处理边界值 epsilon 1e-6; pv_output_clipped pv_output; pv_output_clipped(pv_output_clipped 0) epsilon; pv_output_clipped(pv_output_clipped 1) 1 - epsilon;这个处理看起来简单但epsilon的取值会影响参数估计结果。我试过不同的epsilon发现取1e-3到1e-6之间对拟合结果的影响很小只要不取太大会歪曲分布形状即可。另外一个经验是如果0或1的数据占比很大比如超过20%单纯用Beta分布拟合就不太合适了可能需要零膨胀Beta分布Zero-Inflated Beta这是进阶方向代码里可以留一个扩展接口。3.3 拟合效果如何量化检验光看拟合曲线的形状差不多远远不够。在组合研究中参数估计的精度会直接影响后续联合概率计算的可信度。我推荐至少做两个检验。第一个是Kolmogorov-Smirnov检验直接用Matlab的kstest% KS检验判断数据是否来自指定的Beta分布 pd_beta makedist(Beta, a, alpha_hat, b, beta_hat); [h, p, ksstat] kstest(pv_output, CDF, pd_beta);第二个是Q-Q图分位数-分位数图目视检查尾部拟合情况% Q-Q图 figure; qqplot(pv_output, pd_beta); title(Beta分布拟合Q-Q图);需要说明的是KS检验对样本量很敏感样本量大的时候很容易拒绝原假设。我的经验是不要只看p值是否大于0.05而是结合Q-Q图判断偏差是否在可接受范围内工程建模追求的是足够好不是统计学上完美。4. 两种分布组合研究的三个实用层级4.1 第一层独立假设下的出力叠加最基础的组合方式是假设风速和光照相互独立分别从Weibull分布和Beta分布抽样然后叠加得到系统总出力。% 独立抽样蒙特卡洛 N 10000; % 从Weibull分布抽样风速 wind_samples wblrnd(k, c, [N, 1]); % 风速转出力 pv_wind wind_power_curve(wind_samples, v_cut_in, v_rated, v_cut_out); % 从Beta分布抽样光伏出力 pv_solar betarnd(alpha_hat, beta_hat, [N, 1]); % 总出力 total_output pv_wind pv_solar; % 统计总出力分布 figure; histogram(total_output, 50, Normalization, pdf);这种方法的优点是简单直接Matlab里两行抽样代码就能跑。缺点是忽略了风光的实际关联性结果会低估某些极端场景比如连续阴雨天且风速也低的发生概率。实际工程中把出力等于0的情况单独统计是很有必要的——风小到停机、光照弱到为零两者同时发生的概率就是系统绝对失电的风险。独立假设下这个概率是两者零出力概率的乘积但真实数据往往比这高这就是独立假设的局限。4.2 第二层联合概率密度与互补性量化更深入的做法是构造两变量的联合概率密度。这里有一个很实用的工程近似把风速和光照的关联性通过时间序列相关性引入用Copula函数连接两个边缘分布。但在Matlab代码实现层面可以先从更直观的分箱统计入手用历史数据构造经验联合分布% 假设有历史同步数据wind_power_history 和 pv_power_history % 均为归一化出力序列长度相同 edges 0:0.05:1; [N_joint, ~, ~] histcounts2(wind_power_history, pv_power_history, edges, edges); % 转换为联合概率 P_joint N_joint / sum(N_joint(:)); % 可视化联合概率热力图 figure; imagesc(edges, edges, P_joint); colorbar; xlabel(风电归一化出力); ylabel(光伏归一化出力);这个热力图能直观显示出力组合的高发区和空白区。如果高发区沿对角线分布说明两电源正相关互补性差如果高发区分布在两个轴附近说明互补性好。这个分析不需要复杂数学但工程上非常实用。4.3 第三层Copula函数建模进阶扩展当历史数据不够多或者你希望生成更灵活的联合场景时Copula是标准解法。核心思路是既然Weibull和Beta分别刻画了边缘分布那关联结构就单独交给Copula函数例如Gumbel Copula适合刻画上尾相关极端高出力同时发生的风险。Matlab的Statistics Toolbox里没有直接内置Copula拟合函数但可以用copulafit。需要说明的是系统自带的Copula族主要是椭圆族和阿基米德族。我在实际使用中踩过坑——对于风光这类存在尾部负相关的数据一端高时另一端往往低Clayton Copula有时比Gaussian更贴合。选型逻辑可以简单记为Gaussian是万金油Clayton能刻画一损俱损的下尾相关Gumbel能刻画一荣俱荣的上尾相关。5. Matlab实现中的关键避坑和实测经验5.1 参数估计不收敛先查这三件事我用这套组合模型处理过不少真实数据参数估计不收敛或结果明显不合理时90%是以下三个原因。第一数据没清洗。风速数据里混入负值、异常尖峰、停机时段的数据都会干扰MLE估计。光伏数据里的0和1边界值必须做偏移处理否则betafit会直接报错。第二初值给得太随意。优化算法对初值敏感。Weibull分布的初值可以用矩估计给一个粗略值风速均值μ和标准差σ根据公式k(σ/μ)^(-1.086)、cμ/Γ(11/k)计算。Beta分布的初值也同理% 矩估计法给Beta分布提供初值 mu_hat mean(pv_output); var_hat var(pv_output); common mu_hat * (1 - mu_hat) / var_hat - 1; alpha_init mu_hat * common; beta_init (1 - mu_hat) * common;第三样本量不足或数据时间段单一。只用一个月的冬季数据去拟合全年分布参数肯定偏离。至少用一年的数据且要覆盖四季。5.2 随机抽样时容易忽略的时间相关性在做Monte Carlo模拟时很多人直接用wblrnd和betarnd生成独立随机数。但实际的风速和光照存在显著的自相关性——今天是晴天明天大概率也是晴天这一小时风速大下一小时风速大概率也不小。如果你做的是小时级时序仿真我建议用马尔可夫链或多状态转移矩阵来引入时间相关性而不是直接从静止分布抽样。这个扩展在Matlab里实现也不困难核心是构造转移概率矩阵% 简化示例将光伏出力划分为5个状态 states 1:5; trans_matrix zeros(5); % 转移计数矩阵 for t 1:length(pv_output)-1 s_now discretize(pv_output(t), 0:0.2:1); % 离散化到状态 s_next discretize(pv_output(t1), 0:0.2:1); trans_matrix(s_now, s_next) trans_matrix(s_now, s_next) 1; end % 归一化为概率转移矩阵 trans_matrix trans_matrix ./ sum(trans_matrix, 2);5.3 代码框架的完整搭配建议这里给一个参考的代码目录结构实际跑下来工作效率会高很多wind_pv_combination/ ├── data/ │ ├── wind_speed.csv % 风速时间序列 │ └── pv_output.csv % 光伏归一化出力时间序列 ├── main.m % 主脚本流程编排 ├── fit_distributions.m % 拟合两个分布参数 ├── power_curve.m % 风速转功率曲线函数 ├── joint_analysis.m % 联合概率分析 └── utils/ ├── data_cleaning.m └── plot_save.m主脚本的逻辑就是读数据 → 清洗 → 拟合 → 检验 → 抽样 → 联合分析 → 出图。每个模块单独函数化方便后续替换成别的分布比如用Log-normal替代Weibull做对比研究或接进微电网优化模型。5.4 可视化出图的一点心得组合研究必然涉及大量出图图表的呈现质量直接影响结论的传达。三个图我认为是必须的第一个是风速Weibull拟合曲线与直方图叠加直观展示拟合效果第二个是光伏Beta拟合曲线与直方图叠加第三个是联合概率热力图或散点密度图。出图时推荐统一设置字体和坐标轴范围比如X轴固定0到1归一化出力域方便两图对比。散点图在做大数据量可视化时容易过度绘制可以改用hist3或histcounts2的热力图信息密度更高。6. 这套组合方法还能往哪个方向扩展我目前在做的事情是把这套WeibullBeta组合模型嵌入到微电网储能容量优化里。基本思路是用组合模型生成风光联合出力的蒙特卡洛场景然后对每个场景做储能充放电的时序仿真以失负荷率最小为目标优化储能容量。相比直接用历史数据做场景生成分布模型的好处是可以生成历史数据里没出现过的极端场景让优化结果更稳健。另一个方向是用于风电场和光伏电站的协同选址。不同地理位置的风光资源相关性不同用Copula函数量化这种空间相关性可以指导园区规划——尽量让两个电源的出力时间序列相关性低一些互补性更强一些。回到Matlab实现本身有几个细节最后再叮嘱一下。fitdist函数用起来也很顺手适用于需要更多分布类型对比的场景makedist可以方便地创建分布对象用于后续各种运算。如果你用的是旧版本Matlab比如2016年之前的R2016a部分概率分布对象函数可能不支持建议用wblfit、betafit这类独立的拟合函数兼容性更好。我自己在跑这套代码时有一个很深的体会分布的参数拟合本身只是第一步真正花时间的是理解数据背后的物理过程和工程约束。Weibull分布和Beta分布参数拟合成功的关键是要先理解风电、光电出力的物理特性字面上是数学问题实际上是个工程问题把这点想通了代码反而简单。