做风光联合出力建模这活儿绕不开两个分布风电的Weibull分布和光伏的Beta分布。尤其是做电力系统随机规划、储能容量配置或者可靠性评估的时候输入的风光出力场景不能拍脑袋给个平均值就完事你得把出力的随机特性描述出来。这篇东西我不打算讲大道理就把我实际在Matlab里怎么拟合、怎么组合这两个分布的完整流程、代码和踩过的坑一次性写清楚有直接能抄的代码也有代码背后没写进文档的细节。不管是研究生做课题还是工程师做规划分析按这个思路走一遍基本都能落地。1. 为什么风光输出偏偏要用这两个分布1.1 风速的物理特性决定了Weibull分布的适配性风电出力本质上是风速的函数风速的随机特性直接决定了出力曲线长什么样。我见过不少人上来就直接拿历史出力数据乱拟合其实正确的逻辑链是风速符合某种分布然后经过风力发电机的功率曲线映射到出力。风速本身是正偏态、有下界零的连续随机变量不会出现负值而且绝大多数时候风速集中在一个较低区间偶尔出现高风速的厚尾。两参数的Weibull分布恰好能同时刻画这种偏态和厚尾特性它的概率密度函数长这样[ f(v) \frac{k}{\lambda}\left(\frac{v}{\lambda}\right)^{k-1} e^{-(v/\lambda)^k} ]这里的形状参数k决定了分布曲线的形态k越小尾部越厚风速波动越剧烈尺度参数λ则决定风速的总体量级。实测下来大多数风电场年平均风速对应的k值在2到3之间λ在6到10 m/s之间。为什么不用正态分布因为正态分布是对称的风速不可能对称分布强行用正态拟合会在低风速区出现负值这是物理上不可能存在的后面做抽样就会抽出一堆负出力直接污染整个系统分析结果。1.2 光伏出力的有界性让Beta分布天然占优光伏出力的物理过程是光照强度经过光伏组件转换效率映射到出力。光照强度本身受大气衰减、云层遮挡影响最大值受太阳常数和大气质量限制所以归一化后的光伏出力有一个天然的上界1和一个下界0而且分布形态随着天气状况在极度偏斜和接近对称之间变化。Beta分布正好定义在[0,1]区间上它的概率密度是[ f(x) \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} ]其中α和β是两个形状参数B是Beta函数。当αβ时分布偏向1代表晴朗天气出力偏高当αβ时偏向0代表阴雨天气出力受限。这种有界性是Beta分布最大的优势——它不可能生成超出[0,1]区间的值天然满足光伏出力的物理约束。相比之下如果非要用正态分布或者对数正态分布去拟合归一化光伏出力抽样时必然会产生负值或大于1的值还得做截断处理但一截断概率密度就不再归一化了后面做概率计算全是错的所以Beta分布在这个场景下几乎是唯一合理的选择。1.3 组合研究的核心意义在于捕捉互补特性单独拟合只是一个基础工作把风电和光伏放在一起研究的核心目的是捕捉它们的互补性和相关性。同一个地理位置白天光伏出力高的时候往往风速较弱夜间风速增强但光伏为零这种时间尺度上的互补特性直接关系到系统需要配置多少储能、备用容量如何安排。更关键的是风电和光伏的出力之间并非完全独立大尺度天气系统会同时影响两者所以简单地把两个独立分布相乘得到一个联合分布是不够的需要引入Copula或者场景生成方法来保留它们的相关性结构。这就是组合研究的核心价值所在不是把两个分布摆在那就完事而是要量化它们之间的相依关系。2. Matlab数据准备与预处理这一步错了后面全白搭2.1 数据口径统一和时间分辨率对齐第一步往往是最脏最累的活但恰恰决定整个拟合的成败。风光出力数据来源五花八门有SCADA系统导出的、有气象再分析资料插值过来的还有电网调度系统里面存的这些数据的时间分辨率往往不一致。我用过一个项目数据风电是10分钟一个点光伏是15分钟一个点直接拿来拟合两个分布没问题但做联合场景分析的时候时间轴根本对不上只能重采样。我的建议是统一到1小时分辨率因为电力系统规划和调度分析的标准粒度就是小时级而且小时级数据能滤掉很多短时波动噪声让分布拟合的稳定性更好。时间对齐用Matlab的retime函数配合timetable类型几行代码就搞定不要手工去写循环匹配时间戳又慢又容易错。2.2 风光出力数据的物理量纲归一化拿到数据后必须做归一化而且归一化的基准值要选对。风电出力用风电场额定容量做基准光伏出力用峰值容量做基准这样归一化后的数值区间就在[0,1]之间物理含义是出力率或容量因子。归一化这一步很多人想偷懒不做直接用兆瓦级别的原始数值去拟合Beta分布结果α和β参数会大得离谱后面生成场景再做逆变换的时候数值稳定性极差。实测下来归一化后的数据拟合出来的参数量级通常是零点几到几十数值上处理好很多。代码如下% 读取原始数据假设Data是n行3列的矩阵列依次为时间戳、风电出力、光伏出力 wind_raw Data(:,2); solar_raw Data(:,3); % 额定容量单位MW wind_capacity 150; solar_capacity 100; % 归一化到[0,1]区间 wind_norm wind_raw / wind_capacity; solar_norm solar_raw / solar_capacity; % 剔除超出物理边界的数据传感器异常或检修期间的错误记录 wind_norm(wind_norm 0 | wind_norm 1.05) NaN; solar_norm(solar_norm 0 | solar_norm 1.05) NaN; % 删除NaN行 valid_idx ~isnan(wind_norm) ~isnan(solar_norm); wind_norm wind_norm(valid_idx); solar_norm solar_norm(valid_idx);注意上面我留下了一个1.05的容差因为实际运行中逆变器或变流器在某些工况下短时出力会略微超过额定值超过1就一律杀掉会把一些真实数据误删预留5个百分点的容差是工程实践中的常用做法。2.3 零值堆积问题必须正面处理风电和光伏数据里都有一个让Beta分布拟合很头疼的特征——大量零值。光伏夜间出力严格为零风电无风时段功率也是零如果直接把一大坨零值扔进去做Beta拟合分布参数会被严重拉偏。解决办法有两种第一种是直接剔除零值只对非零区间做Beta拟合第二种是构建混合模型零值单独作为一个离散概率质量点非零部分用Beta分布描述。我建议大多数场景用第一种就够因为风光联合分析关注的核心工况是出力有效的时段夜间光伏零出力是确定性事件不需要用概率分布去描述它硬要建模反而增加复杂度。但注意剔除后要在后面做场景生成的时候把零值时段按时间规律重新补回去否则生成的场景时间序列会和实际不一致。3. Weibull分布拟合的Matlab完整实现3.1 风速数据的两参数极大似然估计拟合Weibull分布最常用的方法是极大似然估计核心思想是找到一组参数λ和k使得观测数据在概率意义下出现的概率最大。两参数Weibull的极大似然估计没有解析解要用数值迭代求解形状参数k的方程[ \hat{k} \left[\frac{\sum_{i1}^{n} x_i^{\hat{k}} \ln(x_i)}{\sum_{i1}^{n} x_i^{\hat{k}}} - \frac{1}{n}\sum_{i1}^{n}\ln(x_i)\right]^{-1} ]这个方程没法直接解出k得靠迭代逼近。Matlab自带的wblfit函数其实就是干这事但作为研究用途我建议自己写一遍迭代过程第一能加深理解第二在遇到异常数据时能更灵活地调整初值和迭代边界。我实际用到的迭代代码如下function [lambda_hat, k_hat] fit_weibull_mle(v) % v是风速数据列向量单位m/s v v(:); % 确保是列向量 n length(v); % 用矩估计给迭代一个靠谱的初值 mu mean(v); sigma std(v); cv sigma / mu; % 经验公式变异系数与形状参数k的关系近似 k0 cv^(-1.086); lambda0 mu / gamma(1 1/k0); % 固定迭代初值 k k0; for iter 1:200 sum_xk sum(v.^k); sum_xk_ln sum(v.^k .* log(v)); sum_ln sum(log(v)); k_new 1 / (sum_xk_ln / sum_xk - sum_ln / n); if abs(k_new - k) 1e-6 k k_new; break; end k k_new; end k_hat k; % 尺度参数的极大似然估计 lambda_hat (mean(v.^k_hat))^(1/k_hat); end这段代码里有两个细节值得说。第一个是初值选取我见过有人直接给k随机初始化结果迭代发散或者收敛到局部极值用变异系数的经验公式算初值之后基本五六次迭代就能收敛到稳定解。第二个是v.^k这种向量化运算Matlab里千万避免用for循环逐元素算幂数据量一特斯拉上去性能差距是数量级的。3.2 最小二乘拟合作为交叉验证手段极大似然估计之外还有一种线性化最小二乘拟合方法它做交叉验证特别好用。Weibull分布的累积分布函数是[ F(v) 1 - e^{-(v/\lambda)^k} ]两边取两次对数可以线性化为[ \ln\left[-\ln(1-F(v))\right] k\ln(v) - k\ln(\lambda) ]这样左边和ln(v)就是线性关系斜率就是k截距结合λ计算。代码如下% 经验累积分布函数值 v_sorted sort(v); n length(v); % 用中位秩公式估计经验CDF比(i/n)更稳健 F_emp (1:n) / (n 1); % 线性化变换 y log(-log(1 - F_emp)); x log(v_sorted); % 线性回归 p polyfit(x, y, 1); k_ls p(1); lambda_ls exp(-p(2) / k_ls); fprintf(MLE结果: k%.3f, lambda%.3f\n, k_hat, lambda_hat); fprintf(LS结果: k%.3f, lambda%.3f\n, k_ls, lambda_ls);注意这里我用了中位秩公式F (i)/(n1)而不是i/n这个细节很关键。如果用i/n最后一个点的经验CDF刚好等于1取对数后就是0再取一次对数变成负无穷线性回归直接崩掉。中位秩是常见的修正法能避开端点处的奇异性。实测两种方法拟合出来的参数差异通常在2%以内如果差异超过5%说明数据质量有问题或者风速分布不适合用两参数Weibull描述需要回头检查数据。3.3 从风速分布到出力分布的功率曲线转换这一步是很多人会跳过的关键步骤。拟合风速的Weibull分布后要得到出力分布必须经过风机功率曲线转换。标准风机功率曲线分三段切入风速V_cin以下出力为零额定风速V_r到切出风速V_co之间出力恒定为额定功率中间段按近似三次方关系映射。我用的转换函数如下function P wind_power_curve(v, P_r, V_cin, V_r, V_co) % v是风速P_r是额定功率 % V_cin切入风速V_r额定风速V_co切出风速 P zeros(size(v)); % 低风速区近似三次方关系 mask_mid (v V_cin) (v V_r); P(mask_mid) P_r * (v(mask_mid).^3 - V_cin^3) / (V_r^3 - V_cin^3); % 额定区 mask_rated (v V_r) (v V_co); P(mask_rated) P_r; % 超出切出风速出力为零保护停机 % 低于切入风速天然为零 end转换后用histogram函数画出力直方图就会发现出力分布往往在零附近和额定功率附近有两个峰值这就是风机工作特性造成的双峰分布。这是很正常的现象不需要强行用某个单一分布去拟合出力数据直接用风速分布加功率曲线转换来描述出力才是物理上正确的做法。4. Beta分布拟合的Matlab完整实现4.1 矩估计法五分钟就能上手光伏出力的Beta分布拟合有两种常用手段矩估计法和极大似然估计法。矩估计的思路是把样本的一阶矩均值和二阶矩方差和Beta分布的参数建立代数关系。Beta分布的均值是α/(αβ)方差是αβ/[(αβ)²(αβ1)]从这两个方程可以反解出参数function [alpha_hat, beta_hat] fit_beta_mom(x) % x是归一化光伏出力数据剔除零值后传入 x x(:); mu_x mean(x); var_x var(x, 1); % 用的是整体方差不是样本方差 % 中间变量 r mu_x * (1 - mu_x) / var_x - 1; alpha_hat mu_x * r; beta_hat (1 - mu_x) * r; % 数值保护参数必须为正 alpha_hat max(alpha_hat, 1e-6); beta_hat max(beta_hat, 1e-6); end矩估计的实现非常直接没有迭代过程不存在收敛性问题在有大量历史数据做快速探索分析时特别方便。它的缺点是估计量的统计效率不如极大似然在样本量小的时候偏差会比较大。实务上我的习惯是先用矩估计快速算一版参数看看数据大致形态再用极大似然做精修。4.2 极大似然估计和Matlab自带函数的使用边界Matlab里Beta分布极大似然拟合直接调betafit函数即可% x是归一化光伏出力剔除零值、剔除夜间数据后 params betafit(x_nonzero); alpha_mle params(1); beta_mle params(2); % 如果要绘制拟合曲线对比 figure; histogram(x_nonzero, 30, Normalization, pdf, FaceColor, [0.8 0.8 0.8]); hold on; x_plot linspace(0.001, 0.999, 200); y_plot betapdf(x_plot, alpha_mle, beta_mle); plot(x_plot, y_plot, r-, LineWidth, 2); xlabel(归一化出力); ylabel(概率密度); legend(实际数据, Beta拟合);betafit函数内部用的是固定点迭代法收敛速度稳定但需要输入数据严格位于开区间(0,1)内。如果数据里有0或者1函数会直接报错或者给出荒谬的参数。这也是为什么预处理阶段必须把零值和满发值单独处理。满发值的情况比零值更隐蔽实际运行中光伏出力长时间稳定在98%、99%附近不会完全等于1但浮点误差可能导致某些点计算出1.0000002这样的值这类越界点也必须清理。4.3 光照强度到出力的链路和分时段拟合策略与风电类似光伏出力的严格物理链路是光照强度经过组件转换模型映射到出力。Beta分布可以直接拟合归一化出力数据但要注意不同时间尺度的数据应该分场景拟合。白天不同时段的光照强度和大气质量差异极大直接拿全年所有时段的出力数据拟合一个Beta分布得到的参数会非常平均无法体现中午出力高、早晚出力低的特征。我实际的项目里会按季节和时间段分桶拟合春夏秋冬各一组每天再按早中晚分三个时段一共12组Beta分布参数。这样做带来的精度提升非常显著缺点是参数数量膨胀后面做场景生成时要设计场景树来调用这些参数。如果嫌复杂折中方案是做日内时段聚合一天分6个时段每个时段一小时级的Beta参数这样既能保留日内变化特征又不至于参数爆炸。5. 风电与光电的分布组合研究怎么把两者桥接起来5.1 联合场景生成的出发点和两种基础方案分布组合研究的落点通常是要生成风光联合出力的随机场景用于后续的随机规划或者可靠性评估。最简单的组合方式是独立抽样假设风电和光伏互不相关各自从自己的分布里独立采样然后按时间戳配对叠加。做法是% 独立抽样生成场景 N_scen 5000; wind_scen_ind wblrnd(lambda_hat, k_hat, N_scen, 1); % 风速分布抽样 wind_power_scen wind_power_curve(wind_scen_ind, P_r, V_cin, V_r, V_co) / P_r; solar_scen_ind betarnd(alpha_mle, beta_mle, N_scen, 1);独立抽样的问题在于它假设风与光完全独立这在大尺度天气系统面前是不成立的。同一个冷锋过境既会带来大风也会带来阴天这种情况下风电出力高而光伏出力低两者之间实际上存在负相关。忽略这种相关性做出来的联合出力场景会高估风光同时出力的概率导致系统备用配置偏低这是有实际安全风险的。5.2 引入Copula模型捕捉相关性结构要正确处理相关性就得用Copula。Copula的本质是把联合分布拆成边缘分布和相关结构两个独立的部分边缘分布你已经用Weibull和Beta拟合好了剩下的问题就是估计相关结构和生成联合样本。具体做法分两步第一步把历史数据转换成均匀分布变量第二步在均匀空间里拟合Copula参数并生成样本。% 第1步用拟合好的边缘分布把数据转化到[0,1]均匀空间 U_wind wblcdf(wind_speed_hist, lambda_hat, k_hat); U_solar betacdf(solar_norm_valid, alpha_mle, beta_mle); % 第2步拟合高斯Copula的相关矩阵 rho copulafit(Gaussian, [U_wind, U_solar]); % 第3步生成联合均匀分布样本 U_sim copularnd(Gaussian, rho, N_scen); % 第4步逆变换回原始物理空间 wind_speed_sim wblinv(U_sim(:,1), lambda_hat, k_hat); solar_sim betainv(U_sim(:,2), alpha_mle, beta_mle);代码里最需要理解的是为什么在Copula里用转换后的U_wind和U_solar而不是直接用原始数据。因为边缘分布已经各自建模完毕相关性分析的唯一目标就是捕捉变量之间的秩相关关系原始数据经过了非线性变换风速到功率、光照到功率Pearson线性相关系数已经失真了但秩相关系数Kendalls tau和Spearman rho对单调变换是不变的。copulafit函数自动处理的就是秩相关所以输入统一的均匀分布变量才是正确的操作。5.3 带约束的风光联合出力场景校正Copula生成的场景在统计特性上和历史数据一致但物理约束还得再检查一遍。风速场景经过功率曲线转换后得到的出力曲线是否与光伏出力协调得看具体的时间尺度。因为风电和光伏的相关性在小时间尺度和大时间尺度上表现完全不同小时级上两者相关性往往较弱日级和季节级上相关性才明显。如果做的是长周期容量规划应该用月平均出力的联合分布来建模时间粒度过细会引入大量噪声掩盖互补特性如果做的是日前调度或实时运行就得用小时级场景并额外叠加时间自相关性否则生成的场景一天之内忽高忽低完全不像真实运行曲线。另外如果生成的联合场景出现极端情况比如风电接近额定且光伏也接近峰值这种场景在物理上不一定不可能但要结合当地天气气候特点判断是否合理。可以加一层场景过滤逻辑计算每个场景的联合出力总和按从小到大的顺序截断到置信水平95%把最极端的5%场景剔除或者标记为特殊场景避免极端场景对优化结果产生过度扰动。6. 常见问题与排查技巧实录6.1 Weibull迭代不收敛或收敛到异常值我自己实际调试的时候遇到过两类情况。一类是数据量太少比如只有几百个小时的风速数据极大似然迭代就容易在少数几次迭代里跳变震荡。解决办法是加迭代上限同时用最小二乘法的结果作为初值再做极大似然如果两者差距很大就说明数据量不足以支撑两参数辨识这时候宁可只用矩估计的结果。另一类是数据含有大量零风速点比如静风时段很长的内陆场址这种情况下Weibull拟合会把k值压得特别小整个分布形态变得极不合理。处理方式是对风速数据做截断处理静风时段单独建模只对有风时段做Weibull拟合。6.2 Beta拟合参数异常与数据预处理排查Beta拟合出现异常参数一般不是算法的问题而是数据的问题。最常见的是alpha和beta都非常小比如小于0.1这说明输入数据的方差接近理论上限分布呈现U型也就是大量数据堆积在0和1两端中间几乎没点。这种情况在光伏出力数据里非常典型全年的日出力数据基本就是夜间0和中午1的混合。解决办法就是按时间区分段拟合白天的时段单独处理不要把所有时段打包。另外还有一个隐蔽问题数据归一化用的容量基准如果小于实际峰值出力会出现大量超过1的数据点Beta分布定义域不包含大于1的值这时估计出来的参数会完全失真。检查方法很简单plot(ecdf(solar_norm))看累积分布有没有在1处出现跳变。6.3 Copula估计不稳定的补救手段高斯Copula的参数估计在高维空间或者数据量有限的时候容易不稳定表现是相关矩阵的某些特征值接近零导致生成的样本各个分量之间有虚假的线性依赖。一个务实的补救方案是改用t-Copula它的尾部相关性比高斯Copula更强能更好地描述极端天气下风光同时偏离均值的情况。在Matlab里实现就是rho_t copulafit(t, [U_wind, U_solar], Tail, 0.5); U_sim_t copularnd(t, rho_t, 0.5, N_scen); % 第二个参数是自由度自由度nu越小尾部相关性越强实际拟合时可以把nu也作为自由参数进行搜索。但注意t-Copula的U_sim_t同样是在[0,1]区间的均匀分布变量逆变换和后面场景生成的流程完全一致。如果连t-Copula都因为样本量太小而不稳定那就退回到独立抽样并显式增加一个相关性修正项把历史数据的秩相关系数作为后验约束叠加到生成的场景上。这个办法虽然理论不够精致但工程上是最稳定的。6.4 代码性能调优的三个实际经验大数据量下拟合迭代的速度直接影响工作效率有几个经验值得分享。一是能用向量化运算的地方坚决不用循环Matlab的向量化优势极大一段v.^k的计算用循环写和向量化写性能差距可能是20倍以上。二是做Monte Carlo场景生成时优先用makedist和random组合而不是直接调用赤裸裸的rand前者在边界处理上更稳健生成速度也没有明显劣势。三是对百万级样本的拟合可以考虑先随机抽样一个子集做参数估计再在全集上验证拟合优度这样既能大幅缩短调试时间又不会明显损失精度。我个人的体会是风光联合出力的分布建模看着像是一个纯粹的统计学问题真正做起来就发现物理机理和数据质量才是决定上限的东西。代码只是把物理过程翻译成数学语言的手段分布选得再漂亮如果数据口径不统一、零值没有处理、时间尺度选错最后的结果照样是一堆没意义的参数。你做完这篇文章里的流程再去跑那些看起来炫酷的规划算法会明显感觉输入数据的底子扎实了结果的可信度自然也就上来了。