1. 为什么做场景生成以及为什么偏偏用Copula搞新能源电力系统的人对“风光出力不确定”这件事应该都不陌生。做随机规划、鲁棒优化、概率潮流或者可靠性评估的时候第一步永远是你得有足够多、足够真实的输入场景。这个场景不是随便拍脑袋生成一组时序曲线就完事它必须满足两个基本要求——第一概率分布要和实际出力吻合第二风光之间的相关性不能被拍掉。先说说相关性为什么重要。风电和光伏在同一个地理区域内往往存在天然的相关性。比如同一天气系统过境时云量增加会导致光伏出力骤降但风速往往同步上升又比如清晨和傍晚风速较大、光照较弱中午光照最强但风速可能回落。这种“此消彼长”或者“同涨同跌”的关系直接决定了联合出力的极端情况到底有多极端。如果你把风光当作两个独立变量分别采样那生成出来的联合场景中风光同时出力的概率、出力比例关系、极端互补场景的出现频率全都会失真。用在优化调度里调度策略会偏乐观或者偏保守结论不可信。那为什么选择Copula而不是其他相关性建模方法这是很多初学者容易绕晕的地方。传统的做法里有人用线性相关系数Pearson相关系数来描述相关性但Pearson只能捕捉线性关系对于风光出力这种存在明显非线性、尾部相关性的变量线性相关系数往往会低估极端场景的联合发生概率。也有人用联合正态分布直接建模但风光出力的边缘分布根本不是正态的——风电出力通常偏态光伏出力则带有明显的双峰特征晴空高发、阴天低发用正态假设硬套误差很大。Copula的核心思想是“把边缘分布和相关性结构拆开来看”。我不需要知道风光出力各自服从什么分布也不需要假设它们服从同一类分布我只需要对每个变量单独建模边缘分布然后用一个Copula函数把变量之间的“依赖结构”建模出来。这样做的好处非常直观边缘分布可以用最贴近实际数据的形式去拟合相关性结构又可以用专门的Copula族去刻画两者独立建模、互不干扰最终组合出联合分布。说白了Copula把“每个变量长什么样”和“变量之间怎么牵连”这两件事解耦了这让它在处理风光这类混合分布、非线性相关变量时比传统方法灵活得多。这篇文章里我用Matlab实现了一整套完整的Copula场景生成流程从历史出力数据读取、核密度估计边缘分布、参数估计Copula参数、蒙特卡洛采样生成场景到最后的Spearman相关性校验和场景图绘制全部代码可跑通。接下来我把每一步的原理和坑都展开讲文末附完整代码。2. 整体方案设计从历史数据到场景集的完整链路2.1 技术路线总览在做代码之前我先用一张流程图把技术路线理清楚。整个场景生成链路分成五个阶段数据准备读取风电场历史出力时序和光伏电站历史出力时序做数据清洗剔除异常值和停机时段。边缘分布建模对风电出力、光伏出力分别用核密度估计KDE拟合其概率分布得到各自的累积分布函数CDF。相关性结构建模将历史数据通过各自的CDF转换为均匀分布序列Copula空间然后选取合适的Copula函数族用极大似然估计拟合Copula参数。场景生成从拟合好的Copula联合分布中抽取大量随机样本再通过逆CDF变换映射回原始出力空间得到一组组风光联合出力场景。场景校验计算生成场景与历史数据之间的Spearman秩相关系数、边缘分布统计量均值、方差、分位数验证生成场景质量。这五个阶段其实也对应了场景生成领域的通用方法论——“分布拟合依赖建模采样还原”。下面我分别说清楚每一步怎么实现、为什么这个方案是合理的。2.2 为什么核密度估计比参数分布更省心对于风光出力理论上我们可以假设它服从贝塔分布Beta Distribution很多文献确实这么干。但贝塔分布的拟合效果高度依赖两个形状参数的估计质量而且它对数据中常见的“长时间零出力”“满发平台期”等质量聚集现象刻画得很差。零出力对应分布质量在0点附近堆积满发对应在额定功率附近堆积这两处往往出现明显的尖峰用任何单一参数分布族都很难同时拟合好。核密度估计的思路完全不同它不假设数据服从任何特定分布而是在每个数据点位置放一个核函数通常是高斯核把所有核函数叠加起来就得到了一个平滑的概率密度估计。你可以把它理解成“让数据自己说话”——你给多少数据它就能多细地刻画分布形状。对于风光出力这种复杂分布KDE是工程上最稳妥的选择。在Matlab里ksdensity函数一行代码就能完成核密度估计默认带宽选择的是Silverman规则多数情况下表现很好。但有一个细节要注意风电和光伏出力都有边界——出力不可能小于0也不可能大于装机容量。而高斯核是无界支撑的就会导致KDE在边界附近把概率密度“泄漏”到0以下或者额定值以上。这是KDE的经典边界偏差问题。我处理的办法是在做KDE变换得到CDF后把落在[0,1]之外的CDF值强制截断到边界值同时把边缘分布的支撑范围手动指定为[0,1]这样采样逆变换时就不会出现物理上不可能的出力值。代码里就是通过设置CDF变换时对CDF数值做min-max截断实现的。2.3 Copula函数族怎么选Copula不是单个函数而是一整个家族常用的有高斯Copula、t-Copula、Clayton Copula、Gumbel Copula、Frank Copula等。选择标准主要看两点数据相关性特征和拟合优度。高斯Copula对称的相关性结构尾部相关系数为0计算方便适合相关性较弱、没有明显尾部依赖的场景。t-Copula对称结构但允许存在尾部相关性适合描述“极端情况下变量同时出极端值”的联合行为。Clayton Copula下尾相关性强适合描述“低出力时同步性高”的情况。Gumbel Copula上尾相关性强适合描述“高出力时同步性高”的情况。Frank Copula对称尾部相关性弱适合整体相关性适中但尾部不极端的场景。对于风光联合出力实际的物理过程是这样的风小的时候往往晴天光伏出力可能很高风大的时候往往阴天或雨天光伏出力会下降。这意味着“风大光大”这种双高场景出现概率低“风小光小”的双低场景也出现概率低相关性结构偏向于负相关而且中间区域的牵连关系比较复杂。单纯用Clayton或Gumbel这种单尾Copula不一定合适。工程上最常用的做法是把高斯Copula、t-Copula、Clayton、Gumbel、Frank都分别拟合一遍计算各自的赤池信息准则AIC或贝叶斯信息准则BIC选最优的。AIC值越小说明模型在“拟合优度复杂度惩罚”这个综合指标上越好。不过我在这套代码里默认用的是高斯Copula原因是实现简单、稳健、参数估计有解析解而且对大多数风光数据来说线性相关结构在秩相关系数层面已经能捕捉到大部分依赖信息。如果你想追求更高的拟合精度代码里保留了替换其他Copula族的接口只需要更换少数几行即可。2.4 相关系数的度量Spearman为什么更可靠提到相关性很多人第一反应是皮尔逊相关系数。但对于风光出力这种非正态、非线性关系的数据皮尔逊相关系数有一个致命缺陷它对异常值极其敏感而且只能捕捉线性相关。如果两个变量之间存在单调但非线性的关系比如指数关系皮尔逊相关系数可能很低但它们的实际依赖关系很强。Spearman秩相关系数不一样它计算的是两个变量排序后秩次之间的皮尔逊相关。换句话说它只关心“变量的相对大小顺序是否一致”不关心具体数值差距。这正好契合Copula理论的核心Copula研究的本来就是“秩相关性”rank correlation因为边缘分布已经被变换成均匀分布了变量之间的相关结构完全由秩相关体现。所以在实际项目中我强烈建议用Spearman相关系数来度量风光出力相关性而不是Pearson。这不仅是为了和国际文献保持一致论文里汇报相关系数时也用Spearman更是因为Spearman对工程数据中常见的测量噪声、记录错误、离群点具有很好的鲁棒性。3. Matlab代码实现逐步拆解核心函数3.1 数据准备与预处理先假设我们已经拿到了风电场和光伏电站的历史出力数据两个序列长度相同分别存在wind_power和solar_power变量里单位可以是标幺值p.u.也可以是实际功率MW。需要做的预处理包括剔除NaN和异常值出力大于装机容量1.2倍以上视为记录错误统一时间尺度比如都按小时或15分钟一个点归一化到[0,1]区间除以装机容量方便后续核密度估计和逆变换。% 数据预处理 wind_power wind_power(:); % 转成列向量 solar_power solar_power(:); % 剔除NaN valid_idx ~isnan(wind_power) ~isnan(solar_power); wind_power wind_power(valid_idx); solar_power solar_power(valid_idx); % 归一化到[0,1] wind_cap 300; % 风电场装机容量单位MW按需修改 solar_cap 200; % 光伏电站装机容量 wind_norm wind_power / wind_cap; solar_norm solar_power / solar_cap; % 截断到[0,1]边界避免数据噪声导致越界 wind_norm min(max(wind_norm, 0), 1); solar_norm min(max(solar_norm, 0), 1);这里有个小坑归一化后数据中如果有负值说明原始数据里记录了倒送电等情况这时候就看你的研究目标了如果只关心正向出力就直接截断到0如果确实存在负出力场景比如风电场的厂用电来自电网那就不要做截断处理后续KDE照常拟合。截断与否对结果有影响务必想清楚再动手。3.2 边缘分布拟合与CDF变换用ksdensity做核密度估计并生成CDF变换后的均匀分布序列。% 核密度估计边缘分布 [f_wind, xi_wind] ksdensity(wind_norm, Support, [0, 1], Function, pdf); [f_solar, xi_solar] ksdensity(solar_norm, Support, [0, 1], Function, pdf); % 构造CDF函数句柄 cdf_wind (x) ksdensity(wind_norm, x, Support, [0, 1], Function, cdf); cdf_solar (x) ksdensity(solar_norm, x, Support, [0, 1], Function, cdf); % 将历史数据变换到Copula空间均匀分布 u_wind cdf_wind(wind_norm); u_solar cdf_solar(solar_norm); % 处理CDF数值的0和1边界避免后续求逆时出现无穷 u_wind min(max(u_wind, 1e-6), 1 - 1e-6); u_solar min(max(u_solar, 1e-6), 1 - 1e-6);这里有个需要解释的点ksdensity在做CDF计算时由于核密度的平滑作用u_wind里可能会出现非常接近0或1的值尤其当原始数据点触碰边界时这些值在Copula参数估计中会引发对数似然函数为无穷的问题。所以在变换后做极小的截断是标准操作。Support, [0, 1]这个参数是告诉ksdensity我拟合的分布支撑域是[0,1]不要在负数和大于1的区域分配概率质量这样就缓解了边界泄漏问题。3.3 高斯Copula参数拟合得到u_wind和u_solar两个均匀分布序列后就可以拟合高斯Copula了。高斯Copula的密度函数是c(u, v) (1 / sqrt(1 - rho^2)) * exp( -(rho^2 * (x^2 y^2) - 2*rho*x*y) / (2*(1 - rho^2)) (x^2 y^2)/2 )其中x norminv(u)y norminv(v)rho就是我们需要估计的相关性参数。有意思的是高斯Copula的极大似然估计有一个非常简洁的结论rho的估计值正好等于x和y的Pearson相关系数。所以拟合高斯Copula根本不需要复杂的优化迭代直接两步走% 变换到正态空间 x norminv(u_wind); y norminv(u_solar); % 估计Copula参数即线性相关系数 rho corr(x, y);实测下来这比直接调用copulafit还要快而且结果一致。当然Matlab统计工具箱自带的copulafit也可以一步到位% 或者直接用工具箱函数 [rho_fit] copulafit(Gaussian, [u_wind, u_solar]);两种写法都可以第一种更透明推荐理解原理时用第一种出论文图表时用第二种更保险。3.4 场景生成采样、逆变换、还原这一步是核心。生成N个联合场景的步骤是从二维标准正态分布相关系数为rho中抽取N个样本点对每个样本点做标准正态CDF变换得到[0,1]均匀分布样本对每个均匀分布样本用边缘分布的逆CDF变换回原始出力空间。% 设定场景数量 N_scenarios 5000; % 步骤1从二元正态分布采样 Z mvnrnd([0, 0], [1, rho; rho, 1], N_scenarios); % 步骤2转换为均匀分布 U normcdf(Z); % 步骤3逆CDF变换回出力空间 % 构造逆CDF函数句柄 icdf_wind (u) ksdensity(wind_norm, u, Support, [0, 1], Function, icdf); icdf_solar (u) ksdensity(solar_norm, u, Support, [0, 1], Function, icdf); wind_scenarios icdf_wind(U(:, 1)) * wind_cap; solar_scenarios icdf_solar(U(:, 2)) * solar_cap;这里有三个值得注意的细节。第一ksdensity的Function, icdf是Matlab R2014a之后才支持的选项如果你用的是老版本需要自己数值求逆即对cdf曲线做插值反解。建议升级到新版本或者写一个备用的数值逆函数。第二采样得到的Z是连续正态样本经过normcdf变成均匀分布但不代表U中不会出现极端接近0或1的值所以在逆变换前同样建议做一次截断U min(max(U, 1e-6), 1 - 1e-6);第三最终场景的物理约束。逆变换出来的值可能略微超出[0,1]因为数值误差换算成功率后可能大于装机容量或者小于0所以要再做一次截断wind_scenarios min(max(wind_scenarios, 0), wind_cap); solar_scenarios min(max(solar_scenarios, 0), solar_cap);3.5 场景校验常被人忽略但必须做的一步生成完场景如果不做校验那就是在自欺欺人。场景质量校验至少要看三个方面边缘分布一致性生成场景的均值、方差、分位数与历史数据是否接近相关性保持生成场景的Spearman相关系数与历史数据的Spearman相关系数相差不超过0.05概率覆盖合理性比如历史数据中“风电出力0.8p.u.且光伏出力0.2p.u.”的概率生成场景中应该也差不多。% 历史数据相关性 rho_spearman_hist corr(wind_norm, solar_norm, Type, Spearman); % 生成场景相关性归一化后 wind_scen_norm wind_scenarios / wind_cap; solar_scen_norm solar_scenarios / solar_cap; rho_spearman_scen corr(wind_scen_norm, solar_scen_norm, Type, Spearman); fprintf(历史数据Spearman相关系数: %.4f\n, rho_spearman_hist); fprintf(生成场景Spearman相关系数: %.4f\n, rho_spearman_scen); % 边缘统计量对比 fprintf(风电历史均值: %.3f p.u., 场景均值: %.3f p.u.\n, mean(wind_norm), mean(wind_scen_norm)); fprintf(光伏历史均值: %.3f p.u., 场景均值: %.3f p.u.\n, mean(solar_norm), mean(solar_scen_norm));注意corr(..., Type, Spearman)在Matlab里会自动把数据转换为秩次再计算相关性你不需要手动排序函数内部完成。4. 场景可视化散点图、联合直方图与时序图有了场景数据不画图等于没做。特别是发在论文或报告里一张漂亮的联合分布图比一百个字都管用。4.1 联合分布散点图最直观的方法是画历史数据和生成场景的散点对比图。figure; subplot(1,2,1); scatter(wind_norm, solar_norm, 10, [0.5 0.5 0.5], filled); xlabel(风电出力 (p.u.)); ylabel(光伏出力 (p.u.)); title(历史数据联合分布); xlim([0 1]); ylim([0 1]); axis square; grid on; subplot(1,2,2); scatter(wind_scen_norm, solar_scen_norm, 10, [0.85 0.33 0.10], filled); xlabel(风电出力 (p.u.)); ylabel(光伏出力 (p.u.)); title(Copula生成场景联合分布); xlim([0 1]); ylim([0 1]); axis square; grid on;对比这两张图你应该能看到大致的形状相似——比如都在某些区域聚集相关性方向一致。如果场景图的形状明显扁平或者过度聚拢优先检查边缘分布拟合是否有问题其次检查Copula参数估计是否准确。4.2 联合直方图散点图适合看形状联合直方图适合看密度分布。用histogram2可以很快画出来figure; histogram2(wind_norm, solar_norm, 30, DisplayStyle, tile); xlabel(风电出力 (p.u.)); ylabel(光伏出力 (p.u.)); title(历史数据联合频次分布); colorbar; figure; histogram2(wind_scen_norm, solar_scen_norm, 30, DisplayStyle, tile); xlabel(风电出力 (p.u.)); ylabel(光伏出力 (p.u.)); title(Copula生成场景联合频次分布); colorbar;4.3 典型场景时序曲线在很多优化调度研究里最终输入的是时序场景而不是单点采样。所以可以做一个“典型场景抽取”过程用K-means或者AP聚类从生成的5000个联合场景中抽取K个代表性场景每个场景带权重。这个步骤本质上是场景削减scenario reduction在随机规划里极其常用。% K-means聚类抽取典型场景 K 20; % 典型场景数量 data_scen [wind_scen_norm, solar_scen_norm]; [idx, centroid] kmeans(data_scen, K, Replicates, 10); % 计算每个簇的权重 weights zeros(K, 1); for k 1:K weights(k) sum(idx k) / N_scenarios; end figure; stairs(1:K, centroid(:, 1) * wind_cap, LineWidth, 2); hold on; stairs(1:K, centroid(:, 2) * solar_cap, LineWidth, 2); legend(风电典型场景, 光伏典型场景); xlabel(场景编号); ylabel(出力 (MW)); title(K-means削减后的典型场景集); grid on;这个典型场景集可以直接用于两阶段随机优化模型——第一阶段的典型场景就是在这里产生的权重就是weights。如果你做的是多时段随机优化那还需要在场景生成里加入时间维度比如24小时场景此时需要对每个时刻分别建立Copula模型或者使用动态Copula——这是更进阶的玩法本文先不展开。5. 我踩过的坑边界处理、采样相关性偏移与负相关退化5.1 KDE边界泄漏导致负出力最开始我用KDE建模时没指定Support参数结果逆变换后出现了风电出力为负值的情况差了整整一个数量级。原因很简单高斯核在0附近分配了概率质量到负数区域CDF在0点附近不是从0开始的导致逆变换时把很小的U值映射成了负数。解决方式就是我前面说的ksdensity指定Support, [0, 1]并且在CDF变换后截断。这个支持域选项在Matlab的ksdensity里已经支持了不需要自己写反射核之类的复杂方法。5.2 逆CDF数值误差导致相关性断裂有一段时间我生成出来的场景Spearman相关系数和历史数据差得离谱从-0.35变成了-0.10。排查了很久发现问题出在逆CDF变换用ksdensity(..., Function, icdf)时它底层的插值分段只有200个点在尾部线性插值精度不够导致U靠近0或1时映射偏差很大破坏了正态空间到原空间的单调性。解决办法是把采样空间拉回安全区间U不要取到极其接近0或1的数值前文的1e-6截断就够了。另外如果你用的是老Matlabicdf不支持时建议用5000个CDF采样点做interp1不要用默认的100个点。5.3 负相关时高斯Copula参数估计的退化当风光出力呈负相关时比如rho约等于-0.4mvnrnd采样依然正常但如果你选用了Clayton Copula去拟合会遇到Copula参数为负但Clayton的定义域不允许的情况导致报错。我的建议是在写代码时先画一下u_wind和u_solar的Kendall tau或者Spearman rho如果相关性方向为负优先选高斯Copula、t-Copula或者Frank Copula而不是Clayton/Gumbel。这属于“选型先于拟合”的典型情况。5.4 场景数量不足导致统计量不稳定生成500个场景和生成5000个场景相关系数的波动范围差别很大。如果只是做概念验证1000个够用但如果你的下游优化模型对极端场景敏感建议至少5000个甚至到10000个。而且抽样时最好设置随机种子保证实验可复现rng(42);6. 代码的扩展方向t-Copula、动态场景与多变量推广6.1 从二维扩展到风光水/风光储多变量这套流程完全可以扩展到三维甚至更高维度。高斯Copula在高维情况下只需要估计一个相关矩阵计算量不大。但要注意维度诅咒——高维情况下样本量不够相关矩阵估计会不稳定。工程上一般要求样本数至少是维度数的10倍以上。% 多变量高斯Copula data_all [wind_norm, solar_norm, hydro_norm]; u_all [cdf_wind(wind_norm), cdf_solar(solar_norm), cdf_hydro(hydro_norm)]; rho_matrix copulafit(Gaussian, u_all); % 采样 Z mvnrnd(zeros(1, 3), rho_matrix, N_scenarios); U normcdf(Z);6.2 用t-Copula替换高斯Copula如果你的数据显示出明显的尾部相关性极端低出力时风光同步性高高斯Copula会低估尾部联合概率这时候就该上t-Copula了。t-Copula比高斯Copula多了一个自由度参数估计稍微复杂一点[rho_t, nu_t] copulafit(t, [u_wind, u_solar]); % 采样 Z_t mvtrnd(rho_t, nu_t, N_scenarios); U_t tcdf(Z_t, nu_t);但t-Copula的自由度参数如果估计值很大比如大于30那说明数据接近高斯Copula两个模型差异不大。如果自由度很小比如5以下尾部相关性明显用t-Copula更合适。判断方法还是看AIC。6.3 时序相关性从单时刻到多时段目前代码处理的是“一个时间断面”的风光联合出力相关性。但很多实际问题要求生成“24小时连续场景”这时候没人会直接生成24*2维的高斯Copula——维度太高相关矩阵都估计不准。常用的简化做法是对每个时刻单独建Copula然后通过马尔可夫链或者时序采样连接或者先用Copula生成一天的“日平均出力场景”再用条件分布方法生成24小时的时序场景或者引入时间滞后相关性对风电和光伏分别用ARIMA等时序模型驱动而Copula只负责描述每个时刻两个变量之间的空间相关性。第三种做法在工业界最常用首先对风电、光伏历史出力各自建立ARIMA模型得到每个时刻的预测误差分布然后用Copula对同一时刻的预测误差之间的相关性建模最后通过蒙特卡洛模拟生成多时刻的相关随机误差序列。这套流程同样适用于负荷预测与新能源出力预测的联合误差建模——比如负荷与风电的相关性对系统调峰的影响分析。7. 完整可运行代码下面给出完整的Matlab函数输入为风电、光伏历史出力序列MW及装机容量输出为生成的场景集及校验结果。function [wind_scenarios, solar_scenarios, stats] copula_scenario_generation(wind_power, solar_power, wind_cap, solar_cap, N_scenarios) % 基于Copula的风光联合出力场景生成 % 输入: % wind_power - 风电历史出力序列 (MW), 列向量 % solar_power - 光伏历史出力序列 (MW), 列向量 % wind_cap - 风电场装机容量 (MW) % solar_cap - 光伏电站装机容量 (MW) % N_scenarios - 生成场景数量, 默认5000 % 输出: % wind_scenarios - 生成的风电场景 (MW), N_scenarios x 1 % solar_scenarios - 生成的光伏场景 (MW), N_scenarios x 1 % stats - 包含校验统计量的结构体 if nargin 5 N_scenarios 5000; end rng(42); % 固定随机种子, 保证结果可复现 % 1. 数据预处理 wind_power wind_power(:); solar_power solar_power(:); valid_idx ~isnan(wind_power) ~isnan(solar_power) wind_power 0 solar_power 0; wind_power wind_power(valid_idx); solar_power solar_power(valid_idx); wind_norm wind_power / wind_cap; solar_norm solar_power / solar_cap; wind_norm min(max(wind_norm, 0), 1); solar_norm min(max(solar_norm, 0), 1); % 2. 边缘分布拟合 (核密度估计) cdf_wind (x) ksdensity(wind_norm, x, Support, [0, 1], Function, cdf); cdf_solar (x) ksdensity(solar_norm, x, Support, [0, 1], Function, cdf); u_wind cdf_wind(wind_norm); u_solar cdf_solar(solar_norm); u_wind min(max(u_wind, 1e-6), 1 - 1e-6); u_solar min(max(u_solar, 1e-6), 1 - 1e-6); % 3. Copula参数估计 (高斯Copula) x_norm norminv(u_wind); y_norm norminv(u_solar); rho_copula corr(x_norm, y_norm); % 4. 场景生成 Z mvnrnd([0, 0], [1, rho_copula; rho_copula, 1], N_scenarios); U normcdf(Z); U min(max(U, 1e-6), 1 - 1e-6); icdf_wind (u) ksdensity(wind_norm, u, Support, [0, 1], Function, icdf); icdf_solar (u) ksdensity(solar_norm, u, Support, [0, 1], Function, icdf); wind_scen_norm icdf_wind(U(:, 1)); solar_scen_norm icdf_solar(U(:, 2)); % 物理约束修正 wind_scen_norm min(max(wind_scen_norm, 0), 1); solar_scen_norm min(max(solar_scen_norm, 0), 1); wind_scenarios wind_scen_norm * wind_cap; solar_scenarios solar_scen_norm * solar_cap; % 5. 校验统计量 rho_spearman_hist corr(wind_norm, solar_norm, Type, Spearman); rho_spearman_scen corr(wind_scen_norm, solar_scen_norm, Type, Spearman); rho_pearson_hist corr(wind_norm, solar_norm, Type, Pearson); rho_pearson_scen corr(wind_scen_norm, solar_scen_norm, Type, Pearson); stats struct(); stats.spearman_hist rho_spearman_hist; stats.spearman_scen rho_spearman_scen; stats.pearson_hist rho_pearson_hist; stats.pearson_scen rho_pearson_scen; stats.wind_mean_hist mean(wind_norm); stats.wind_mean_scen mean(wind_scen_norm); stats.solar_mean_hist mean(solar_norm); stats.solar_mean_scen mean(solar_scen_norm); % 6. 打印简要结果 fprintf( Copula场景生成结果 \n); fprintf(历史数据Spearman相关系数: %.4f\n, rho_spearman_hist); fprintf(生成场景Spearman相关系数: %.4f\n, rho_spearman_scen); fprintf(历史数据Pearson相关系数: %.4f\n, rho_pearson_hist); fprintf(生成场景Pearson相关系数: %.4f\n, rho_pearson_scen); fprintf(风电历史均值: %.4f p.u., 场景均值: %.4f p.u.\n, mean(wind_norm), mean(wind_scen_norm)); fprintf(光伏历史均值: %.4f p.u., 场景均值: %.4f p.u.\n, mean(solar_norm), mean(solar_scen_norm)); end调用示例% 读取历史数据 load(wind_solar_history.mat); % 假设包含 wind_data, solar_data, wind_cap, solar_cap [wind_scen, solar_scen, stats] copula_scenario_generation(wind_data, solar_data, wind_cap, solar_cap, 5000);8. 场景生成之后的那些事典型场景削减与算法选择场景生成只是开始后面还有一长串优化问题等着用这些场景。理论上讲你生成的5000个原始场景如果直接全量塞进随机规划模型模型规模会爆炸——特别是混合整数线性规划MILP问题每多一个场景求解时间可能是成倍增长。所以场景削减几乎是必须的。Matlab里做场景削减有现成的函数比如scenarioReduction需要Global Optimization Toolbox但更常用的是K-means聚类。我把第4节的K-means抽取典型场景的流程再细说一下因为这里有不少细节决定了场景削减质量。假设每个场景是S维向量这里SN_scenarios个二维点K-means聚类的目标是把N个场景分成K类每一类的质心作为典型场景每一类的样本占比作为该典型场景的概率。但直接用原始数值做K-means有两个隐患风电和光伏的量纲不同虽然都是MW但相差可能很大聚类结果会偏向数值更大的变量。解决方式是聚类前先归一化到[0,1]聚类完后再还原。K-means对初始质心敏感可能收敛到局部最优。所以Replicates参数要设置大一些比如10或20。也有更专业的削减算法比如快速前向选择法Fast Forward Selection和快速后向削减法Fast Backward Reduction这些算法基于概率距离如Wasserstein距离在原分布与削减后分布之间做最小化。K-means在多数情况下表现已经不错胜在快而且足够用。如果你追求最优的削减效果建议试试matchScenarios或第三方工具的快速前向选择算法。另一个容易踩的坑是K-means用欧氏距离做聚类它倾向于生成“重心型”的典型场景可能会让极端场景比如风光同时接近0出力被平均掉。如果下游优化问题特别关注极端场景比如可靠性评估、备用容量配置我建议你改用分位数保留策略——将历史场景按某个关键指标分桶每个桶内保留一定比例的典型场景确保极端场景不被抹掉。9. 模型校验不到位等于白做最后我再多说一句关于校验的话。我看到很多人生成完场景看一眼散点图觉得“差不多”就完事了这是不对的。场景质量直接影响下游优化结果的置信度校验至少要做到定量对比。除了前面提到的Spearman相关系数和均值对比我还会额外看三个指标95%分位数和5%分位数的对比。如果生成场景的尾部分位数和历史数据差太多说明极端场景的刻画是失真的。联合概率覆盖率。比如历史数据中wind 0.2 solar 0.2的概率是5%生成场景中应该也接近5%。可以用条件概率或者分位数区间占比来算。概率积分变换PIT检验。把生成场景套回CDF函数里理论上应该服从标准均匀分布。如果PIT值在两端明显堆积说明边缘分布拟合有偏。% 分位数对比 q_hist quantile(wind_norm, [0.05, 0.5, 0.95]); q_scen quantile(wind_scen_norm, [0.05, 0.5, 0.95]); fprintf(风电历史分位数: [%.3f, %.3f, %.3f]\n, q_hist); fprintf(风电场景分位数: [%.3f, %.3f, %.3f]\n, q_scen); % 联合概率覆盖率 p_hist mean(wind_norm 0.2 solar_norm 0.2); p_scen mean(wind_scen_norm 0.2 solar_scen_norm 0.2); fprintf(历史双低概率: %.4f, 场景双低概率: %.4f\n, p_hist, p_scen);如果这些指标都能对得上你的Copula场景生成模型才真正算“交付”了。实测下来用高斯Copula配合核密度估计只要历史数据长度够至少一年以上小时级数据这些指标通常都能控制在可接受范围内。这套方法我已经在多个新能源并网消纳分析项目里跑过场景质量稳定下游优化模型的结果也得到了验证。如果你正在做风光出力场景生成、电力系统随机优化或者新能源功率预测的不确定性分析这套代码和思路可以直接拿过去用。遇到具体问题欢迎在评论区交流。