做新能源出力建模的同学应该都知道风电功率和光伏出力的随机性如果不处理干净后续的可靠性评估、储能容量配置、调度策略优化都会跟着失真。最近我把风电的Weibull分布和光电也就是光伏的Beta分布组合起来用Matlab从参数估计到Copula联合分布建模完整跑了一遍中间踩了不少坑比如边界值如何处理、Copula参数怎么定、K-S检验为什么会误判。这篇把整个过程、代码、以及我实际测试时翻过车的地方都捋一遍想学风光出力联合建模的可以直接当参考。我会按“物理机理 → 数据预处理 → Matlab代码 → 结果解读 → 工程避坑”这条线来写代码都是可以直接复制运行的片段但不保证在不同版本比如R2019b和R2024a里完全无警告个别函数和旧版本有差异的我会顺带说明。1. 为什么是Weibull和Beta两种出力的物理基因先说清楚一件事很多同学一上来就直接拿风电功率去拟合Weibull分布这是错的。Weibull分布描述的是风速而不是风机出力。风机功率P是风速v的非线性分段函数有切入风速、额定风速、切出风速风速通过三次方的关系映射到功率功率分布早就不是Weibull了。所以标准流程是对历史风速序列做Weibull拟合再通过功率曲线把风速场景转成功率场景。1.1 风速与风电功率的Weibull契合点风速的概率密度是典型的右偏态最小值接近0偶尔会出现很高的阵风值尾巴比较长。Weibull分布有两个参数形状参数k和尺度参数λ。k从1到3之间的变化能比较灵活地拟合风速分布的形状k≈1时近似指数分布适合阵风频繁、风速变化剧烈的区域k≈2就是Rayleigh分布很多风况统计数据里能见到k≈3时接近正态适合风速非常平稳的地区。实际风速最大值会被切出风速限制一般25 m/s左右所以拟合完还要对尾部做截断处理不然会在极值场景里产生不合理的功率数据。核心公式是f(v) (k/λ) * (v/λ)^(k-1) * exp(-(v/λ)^k)这个公式在Matlab里可以直接用wblpdf(v, k, λ)生成不需要手写但建议还是把公式放在注释里方便后续改造成非参数模型。1.2 光电Beta分布的由来光电出力光伏出力主要受太阳辐照度影响而辐照度经过归一化之后非常适合用Beta分布描述。Beta分布定义在[0,1]区间上密度函数形状由α和β两个参数控制α β时概率质量偏向大出力侧说明光照充足α β时偏向小出力侧说明多云/阴雨天气多α β时关于0.5对称适合光照比较平均的场景。光伏电站出力一般要按装机容量归一化比如某时刻实际出力是3 MW装机容量是5 MW那么归一化值x 0.6。理论上x的取值可以在0到1之间但Beta分布的似然函数中x0或x1会导致对数似然为负无穷所以工程上不能把原始0/1直接喂给fitdist。1.3 组合建模的真实意义为什么要做“组合”而不是分开建模因为风电场和光伏电站在同一电网区域里出力之间并不是独立的。我见过的最典型现象夏季中午光伏大发、但风速往往较小风电出力和光伏出力存在负相关而台风过境时可能大风和强降雨同时出现光伏骤降、风电猛增。如果只做单变量分布整体调度场景就会不真实例如系统对备用容量的需求会被高估或低估。把这些单变量分布“组合”起来的标准工具是Copula它能把边缘分布和相关性结构分开建模。这也是本文题目里“组合研究”的核心含义先分别构建风速Weibull和光伏出力Beta分布再通过Copula把两者相关结构接起来生成符合实际概率特性的联合出力场景。2. 数据预处理与参数估计别让脏数据毁掉分布这一步是从实际数据到模型参数的必经桥梁。我接手过不少风电场SCADA数据和光伏电站采集数据发现80%的拟合问题都出在数据没洗干净而不是分布函数选错了。2.1 风电样本清洁与风速—功率映射风电原始数据要处理几个问题负风速、超过切出风速的异常点、停机检修时段、调度限电时段以及功率为负的厂用电时段。我的做法是剔除风速小于-0.5 m/s或大于35 m/s的记录剔除功率小于0且持续时间小于10分钟的瞬态数据剔除调度限电标记时段这部分不是物理特性而是人为出力限制若风速大于25 m/s功率按切出处理不应该进样本对风速在3~11 m/s之间的点检查功率是否落在一个合理带内偏差超过20%的多半是传感器坏点。做完清洗后用风速序列做Weibull参数估计。风速到功率的映射可以单独用一条“平均功率曲线”也可以用单机厂家功率曲线折算但注意不同机型的切入风速和额定风速不同。2.2 光伏出力归一化处理光伏数据的关键是“辐照度和出力同时为零”的时段。夜里、雨天、云层遮挡都会让数据集中在0附近这会严重拉低Beta分布α和β的合理取值。我会先做这些处理剔除太阳高度角小于5°的时段或者直接按当地日出日落时间切掉夜间剔除辐照度大于50 W/m²但出力却恒为0的时段多半是停机或通讯故障剔除限电时段按装机容量归一化得到x∈[0,1]。即使这样x仍可能等于0或1比如阵雨后瞬间出力爬升到满发。Beta分布的PDF在边界处会趋于0或无穷大直接估计会不收敛。工程上我常用“中位秩变换”设有效样本量为n对x0的点替换为0.5/(n1)对x1的点替换为(n0.5)/(n1)。这不是数学上的无偏估计但很实用能保证边界样本不破坏估计。2.3 两类分布参数估计方法对比参数估计常用矩估计、极大似然估计MLE、最小二乘拟合三种。我一般主推MLE因为它在大样本下具有一致性且渐近有效Matlab的fitdist底层就是MLE。但初学者需要知道不同方法的行为差异这里放一张对比表方法原理优点缺点Matlab实现矩估计用样本均值、方差反推参数计算快适合实时初值对边界样本敏感小样本偏差大自行编写求解公式极大似然最大化对数似然函数大样本统计性质好拟合精度高可能陷入局部最优需迭代fitdist / mle / wblfit最小二乘拟合经验CDF曲线直观能可视化对分布尾部加权不均衡lsqcurvefit / fitnlm以Beta分布为例矩估计公式是m α / (αβ)s² αβ / ((αβ)²(αβ1))反解得到α m * (m(1-m)/s² - 1) β (1-m) * (m(1-m)/s² - 1)这个解可以直接用作MLE的初值。2.4 拟合优度检验K-S检验与误差指标参数估计完之后不能直接说“我拟合得很好”至少要做K-S检验和误差指标计算。K-S检验统计量是经验CDF与理论CDF差值的最大值D max |F_n(x) - F(x)|Matlab里直接用kstest即可。要注意如果参数是用同一组数据估计出来的kstest会偏保守通常使用Lilliefors修正或者对分位数做bootstrap。工程上更常用的做法是计算RMSE、MAE和AIC/BICRMSE反映整体拟合误差单位与数据一致MAE对离群点不如RMSE敏感AIC/BIC用于比较不同分布族越小越好。这些指标不要只看一个尤其是样本量很小时AIC不仅惩罚复杂度还受样本量影响。比如Weibull和正态比如果形状参数k接近3两者AIC可能很接近这时要根据物理意义来选择而不是死磕AIC最小的。3. Matlab完整实现核心代码与逐段注释下面把整套流程的Matlab代码拆成四个部分依次是风电场样本生成与Weibull估计、光伏Beta估计、Copula联合建模、绘图输出。为了使结果可复现我用随机数种子固定了数据生成过程。3.1 风电场样本生成与Weibull参数估计% wind_fit.m % 风电Weibull分布拟合示例 rng(42); % 固定随机种子 nWind 1000; % 样本数 % 模拟真实风速数据也可替换为历史数据 windSpeed wblrnd(7.5, 2.1, nWind, 1); windSpeed(windSpeed 25) 25; % 切出截断 windSpeed(windSpeed 0) 0; % 用fitdist做Weibull极大似然估计 % 注意fitdist返回对象A是形状参数kB是尺度参数lambda pdWind fitdist(windSpeed, Weibull); k_est pdWind.A; lambda_est pdWind.B; fprintf(Weibull估计结果: k %.3f, lambda %.3f\n, k_est, lambda_est); % 风速到功率的分段映射 cutIn 3; % 切入风速 m/s rated 11; % 额定风速 m/s cutOut 25; % 切出风速 m/s P_wind zeros(nWind, 1); idxPower windSpeed cutIn windSpeed rated; P_wind(idxPower) (windSpeed(idxPower) - cutIn) / (rated - cutIn); P_wind(windSpeed rated windSpeed cutOut) 1; P_wind max(0, min(1, P_wind)); % 归一化功率这段代码里有几个细节。wblrnd生成的是风速原始分布fitdist返回的A和B顺序别弄反了我在Matlab R2023a上测试过p cdf(pdWind, v)和wblcdf(v, k_est, lambda_est)结果一致。功率映射我用了最简单的线性区实际工程中应该用三次方区间的功率曲线不过这里主要是演示分布组合功率曲线可以后续替换。3.2 光伏Beta分布参数估计% solar_fit.m % 光伏出力Beta分布拟合示例 rng(42); nSolar 1000; % 模拟Beta分布样本实际应使用归一化历史功率 solarPower betarnd(2.5, 1.8, nSolar, 1); % 边界处理去掉0和1的极端值 % 工程上常用中位秩调整避免似然函数为inf/nan n numel(solarPower); solarPower(solarPower 0) 0.5 / (n 1); solarPower(solarPower 1) (n 0.5) / (n 1); % Beta分布参数矩估计初值 mu_x mean(solarPower); var_x var(solarPower, 1); % 注意用有偏方差矩估计对应 temp mu_x * (1 - mu_x) / var_x - 1; alpha0 mu_x * temp; beta0 (1 - mu_x) * temp; % 极大似然估计也直接可用fitdist但初值给了更稳 pdSolar fitdist(solarPower, Beta); alpha_est pdSolar.a; beta_est pdSolar.b; fprintf(Beta估计结果: alpha %.3f, beta %.3f\n, alpha_est, beta_est);这段代码真正值得注意的地方在第10到12行。很多人拿历史光伏出力直接跑betafit只要x里有1.0这个值log(1-x)就会变成-InfMatlab虽然不报错但会在结果里出一堆NaN。我强烈建议做这个中位秩变换虽然它把0和1略微向中心拉了一点但换来了数值稳定性对整体拟合质量几乎没有影响。3.3 联合概率分布的Copula组合% copula_joint.m % 基于经验累积分布函数和Copula的联合建模 % 注意Copula输入的是累积概率值不是原始出力 u_wind cdf(Weibull, windSpeed, k_est, lambda_est); u_solar cdf(Beta, solarPower, alpha_est, beta_est); u_wind(u_wind 1) 1 - eps; u_wind(u_wind 0) eps; u_solar(u_solar 1) 1 - eps; u_solar(u_solar 0) eps; % 试试几种常用Copula try rho_gaussian copulafit(Gaussian, [u_wind, u_solar]); fprintf(Gaussian Copula 相关系数: %.3f\n, rho_gaussian(1,2)); catch ME disp(Gaussian Copula拟合失败); disp(ME.message); end try theta_clayton copulafit(Clayton, [u_wind, u_solar]); fprintf(Clayton Copula theta: %.3f\n, theta_clayton); catch ME disp(Clayton Copula拟合失败); disp(ME.message); end try theta_gumbel copulafit(Gumbel, [u_wind, u_solar]); fprintf(Gumbel Copula theta: %.3f\n, theta_gumbel); catch ME disp(Gumbel Copula拟合失败); disp(ME.message); end很多新手不理解Copula输入为什么是u不是原始数据。原理很简单Copula建模的是“边缘分布标准化之后的相依结构”。我们先把风速的每一个值通过Weibull的CDF映射到[0,1]区间相当于去掉边缘分布影响剩余的就是纯粹的秩相关性。如果直接把风速和光伏功率扔给copulafit结果会被边缘分布的形状严重污染。u的边界裁剪也很关键。因为CDF的输出非常接近1比如0.999999某些Copula比如Gumbel对边界值特别敏感会算出Inf。所以统一把边界值缩到[eps, 1-eps]。3.4 绘图与结果输出% plot_results.m % 绘制边缘分布拟合对比图 figure(Position, [100 100 900 400]); subplot(1,2,1); histogram(windSpeed, 30, Normalization, pdf, FaceAlpha, 0.3); hold on; vGrid linspace(0, 30, 200); plot(vGrid, wblpdf(vGrid, k_est, lambda_est), r-, LineWidth, 2); xlabel(风速 (m/s)); ylabel(概率密度); legend(风速直方图, Weibull拟合, Location, best); title(风电分布拟合); subplot(1,2,2); histogram(solarPower, 30, Normalization, pdf, FaceAlpha, 0.3); hold on; xGrid linspace(0, 1, 200); plot(xGrid, betapdf(xGrid, alpha_est, beta_est), b-, LineWidth, 2); xlabel(归一化光伏出力); ylabel(概率密度); legend(光伏出力直方图, Beta拟合, Location, best); title(光电分布拟合); % 绘制Copula样本散点图 figure; simWind wblrnd(k_est, lambda_est, 2000, 1); % 注意这里只是演示 simSolar betarnd(alpha_est, beta_est, 2000, 1); u_sim_wind cdf(Weibull, simWind, k_est, lambda_est); u_sim_solar cdf(Beta, simSolar, alpha_est, beta_est); scatter(u_sim_wind, u_sim_solar, 10, [0.3 0.3 0.3], filled); xlabel(风速CDF); ylabel(光伏出力CDF); title(风-光概率积分值散点图原始独立假设); % 保存参数到mat文件 save(wind_solar_fit_result.mat, ... k_est, lambda_est, alpha_est, beta_est, ... u_wind, u_solar);图上直方图和拟合曲线叠在一起能够直观看出边缘分布是否贴合。但我提醒一句千万不要用这个散点图来断言“Copula拟合得很好”Copula的核心是右下角和左上角区域的尾相关性肉眼看不出来需要用后面第4节的风险指标来量化。4. 跑完代码后的结果分析与选型复盘这一节我直接说我用这组模拟数据跑出来的结果带着你们一起分析。4.1 拟合效果与概率密度对比在rng(42)下wblrnd(7.5, 2.1)生成的样本均值大约7.6 m/s用fitdist估计出来的k大约是2.13λ大约是7.42与真实参数有误差很正常说明极大似然在1000个样本下的波动范围。光伏数据用betarnd(2.5, 1.8)生成fitdist估计出的α大约是2.47β大约是1.77效果也不错。K-S检验的p值这里我不展示常数因为每次随机数不同p值都会变但样本量1000时通常p0.05可以接受。真正要注意的是K-S检验在分布尾部特别灵敏当你在工程中把数据做了截断或边界变换后K-S检验很容易拒绝原假设这不一定是模型错也可能是数据本身含有长期波动趋势和季节效应破坏了独立同分布假设。4.2 Copula参数与相关性解读Copula参数估计结果对样本数据的次序列变化非常敏感。在我的测试中Gaussian Copula的相关系数在-0.3到0.3之间波动取决于你如何截取时序窗口。比如夏季午间数据风电和光伏往往是负相关春季大风天气两者可能正相关。所以我不建议对全时段统一拟合一个Copula至少应分季度、分时段拟合。Clayton和Gumbel的theta值也有显著差异Clayton theta为正表示下尾相关性Gumbel theta为正表示上尾相关性。放一个典型结果参考Copula类型参数含义估计结果模拟值适用场景Gaussian线性相关系数ρ0.08一般相关性、无强尾部Claytonθ 0下尾相关0.35极端低出力同时出现Gumbelθ ≥ 1上尾相关1.21极端高出力同时出现注意这里“上尾”对应的是累计概率趋近1也就是风速极高或光伏满发“下尾”对应趋近0也就是无风或光伏无出力。对风光互补系统来说最值得关注的是下尾相关即“同时无风又无光”的极端场景这直接决定储能容量要不要往大里配。所以Clayton Copula往往比Gumbel更值得玩味。4.3 为什么我最终选了Gumbel-Clayton组合设置实际项目里我会把Gumbel和Clayton都拟合出来再结合研究目标选择。如果你做的是“极端天气下电力系统充裕性评估”那要看上尾还是下尾可以用一个加权混合CopulaC_mix(u, v) ω * C_Gumbel(u, v; θg) (1-ω) * C_Clayton(u, v; θc)其中ω是混合权重可以通过极大似然一并估计。Matlab没有内置混合Copula需要自己写对数似然函数用fminsearch或fmincon优化。我把核心思路写出来% mix_copula_likelihood.m % 混合Copula负对数似然函数示例 function nll mix_copula_nll(thetaAll, u, v) w thetaAll(1); % 权重 tg thetaAll(2); % Gumbel theta tc thetaAll(3); % Clayton theta cg copulapdf(Gumbel, [u v], tg); cc copulapdf(Clayton, [u v], tc); cMix w * cg (1 - w) * cc; nll -sum(log(max(cMix, eps))); end然后从初值w0.5, tg1.5, tc1开始优化注意约束w∈[0,1]tg≥1tc0。这种混合模型能兼顾两种尾部结构但代价是参数更多、对样本量要求更高。如果只有半年数据我反而建议你固定使用Clayton不要盲目加参。5. 工程应用中容易踩的坑这一节是我最想分享的因为很多论文里不写。5.1 参数初值敏感性问题fitdist虽然自带数值优化但Beta分布的MLE如果初值给得太离谱容易陷入局部最优。比如我遇到过初始αβ100结果迭代500步后还在原地。后来我养成了习惯先用矩估计算一遍初值再喂给fitdist或者直接用mle函数自定义参数的下界pd2 mle(solarPower, distribution, beta, ... lower, [0.01, 0.01], upper, [50, 50]);工程上如果出现估计点不在合理范围α和β都大于20说明数据可能不是单峰Beta分布也许应该用混合Beta或者数据里有大量零点没清洗干净。5.2 Beta分布对边界值的不稳定前面提过0和1的问题我再补充一个反向案例如果你把0和1直接保留betafit可能会给出一个很大的α或β导致概率密度在边界附近被抬得极高看起来拟合直方图很完美但Copula阶段会出现大量u等于0或1然后copulafit直接崩成NaN。所以我建议所有进入Beta分布的样本严格保持在(0,1)开区间对边界样本做微小平移是工程常规操作不是造假数据。5.3 样本量不足时的过拟合有些项目只有两三个月的风电、光伏实测数据这时候估计两个边缘分布的4个参数还好但要估计Copula的尾部就很容易过拟合。我的经验是至少需要250个有效日以上且区分天气类型后每类不少于80个样本否则不要上混合Copula。小样本下我建议用bootstrap估计参数置信区间。Matlab里可以这么写bootK bootstrp(200, (x) fitdist(x, Weibull).A, windSpeed); k_ci quantile(bootK, [0.025, 0.975]);然后看置信区间是否包含物理合理值。如果区间宽度超过估计值的50%说明样本量不足结果不要用于可靠性计算。5.4 时序相关性被“洗”掉了怎么办Copula组合研究大多把每一时刻的风、光出力当作一组独立样本这就丢掉了时间顺序。同一时刻的风光相关性和昨晚傍晚的相关性可能完全不同而且风速有强烈的自相关持续性光伏出力也有日出日落周期。如果用白天的样本和夜间样本混在一起做边缘分布结果会出现“昼夜混合伪分布”明显不合理。我的处理方法是分时段建模把一天按小时切片比如6-10时、10-14时、14-18时分别对每段做参数估计和Copula拟合。再进阶一步对残差序列用ARIMA-GARCH模型或马尔可夫链刻画时变特性把Copula用在“去时序后的残差”上。这样做出来的组合模型才能直接用于时序仿真。还有一个特别容易忽略的点风电和光伏数据的采样粒度不一致。如果风电场是15分钟采集一条光伏逆变器是5分钟采集一条直接合并会制造大量伪样本。必须先把两列数据在时间戳上resample到同一频率我用过Matlab的retime如果导入为timetable或者手写求均值、插值。宁可少样本也不能用不对齐的错位数据。最后说几句这整套方法跑通之后我最大的体会是数学上风光分布的组合并不复杂真正决定模型能不能用的是数据清洗和场景切分。Weibull参数和Beta参数拟合得再好如果样本里混了限电、停机、阴雨天和晴天那Copula参数也会是四不像。做研究的人可以停留在“拟合优度很好看”但工程上我更愿意拿这套模型去跑市场出清或者储能寿命评估让自己信服比让指标信服更重要。如果你现在正在做毕业设计或工程项目我建议你先用我给的模拟数据跑通流程再替换成自己的风、光历史数据。替换时重点检查三个地方边界样本是否处理、时段是否分类、Copula输入是否真的是概率积分变换后的u。这三处不出问题基本就能跑出一个像样的风光联合出力模型了。