简介本资源为基于Copula理论与K-means聚类的风光联合场景生成MATLAB代码包面向计算机、电子信息工程、数学等专业的学生及科研人员可用于课程设计、期末大作业与毕业设计中的新能源出力场景建模与聚类分析。压缩包共15个文件约4.77MB包含1个main.m主程序、2个xlsx风光功率数据文件、3份Copula与场景生成相关PDF学习资料、8张png运行结果图及说明文本结构清晰、注释明细支持参数化编程与灵活改参。代码兼容MATLAB 2014、2019a与2024a版本附赠案例数据可直接运行帮助读者快速复现从边缘分布拟合、Copula联合分布构建到K-means场景聚类的完整流程理解变量间复杂依赖关系的建模思路。目前已有122人学习下载适合需要掌握场景生成方法并完成相关课题的中高级学习者参考。1. 风光联合场景生成到底在算什么从 Copula 到 Kmeans 的一条完整链路做新能源出力建模的人迟早会撞上一个问题风电场和光伏电站的出力曲线单独看都还好一旦要放在同一个系统里做联合调度或可靠性评估两条曲线之间的相关性就成了绕不过去的坎。你如果直接把风电和光伏的历史出力序列各自拟合一个边缘分布再独立采样拼在一起得到的场景集合会严重失真——因为现实中太阳落山前后光伏骤降的那段时间风速往往也在变化两者不是各走各的。这就是风光联合场景生成要解决的核心问题在保留每个变量自身统计特性的同时把变量之间的相关结构也复现出来。标题里的 64Copula 指的是一种具体的 Copula 建模方案Kmeans 聚类则负责把生成出来的海量场景压缩成少量典型场景供后续优化模型使用。整条链路是历史出力数据 → 边缘分布拟合 → Copula 相关结构建模 → 采样生成大量场景 → Kmeans 聚类削减 → 输出典型场景及其概率。适合做电力系统规划、微网容量配置、储能优化调度的工程师也适合正在写相关方向论文、需要一套能跑通的 MATLAB 代码做实验的研究生。下面我把这条链路拆开从原理到代码到踩坑一步步讲清楚。2. Copula 建模风光相关性的原理与 MATLAB 实现2.1 为什么不能用协方差矩阵直接描述风光相关性很多人第一反应是用多元正态分布来建模风电和光伏的联合分布协方差矩阵直接描述相关性。这个思路在变量服从正态分布时没问题但风电出力典型分布是 Beta 或 Weibull光伏出力在白天呈 Beta 分布、夜间恒为零两者都不是正态的。硬套多元正态边缘分布就被扭曲了尾部相关性也完全对不上。Copula 的思路是把联合分布拆成两部分边缘分布各自拟合变量之间的依赖结构单独用一个 Copula 函数描述。Sklar 定理保证了这种拆分的数学合法性——任何联合分布都可以写成边缘分布和 Copula 函数的组合。这样一来风电用 Weibull 拟合、光伏用 Beta 拟合两者之间的相关结构用 Gaussian Copula 或 t-Copula 来刻画各管各的互不干扰。选 Gaussian Copula 还是 t-Copula如果风光出力序列在极端天气下表现出明显的尾部依赖比如连续阴天伴随低风速t-Copula 更合适因为它有更厚的尾部。如果数据量不大、尾部特征不明显Gaussian Copula 参数少、估计稳定是更稳妥的起点。我一般先用 Gaussian 跑通再看尾部拟合效果决定要不要换。2.2 边缘分布拟合风电 Weibull、光伏 Beta 的参数估计先看风电。风速通常用 Weibull 分布拟合但注意这里要建模的是风电出力而非风速。如果只有出力数据可以直接对出力做 Beta 拟合或者用风速-功率曲线反推。常见做法是对归一化后的风电出力序列做 Beta 拟合形状参数用最大似然估计。光伏出力白天归一化后同样适合 Beta 分布。但夜间出力为零这部分不能直接丢——丢了会改变日内的时序结构。处理方式有两种一是只对白天时段建模夜间单独标记为零二是对全时段做零膨胀 Beta 拟合。代码里通常采用第一种简单且不容易翻车。% 假设 wind_pu 和 solar_pu 是归一化后的风电和光伏出力列向量 % 去除夜间零值只对光伏白天时段拟合 day_idx solar_pu 0; solar_day solar_pu(day_idx); % 风电 Beta 拟合 wind_beta fitdist(wind_pu, Beta); a_w wind_beta.a; b_w wind_beta.b; % 光伏 Beta 拟合仅白天 solar_beta fitdist(solar_day, Beta); a_s solar_beta.a; b_s solar_beta.b; % 输出参数供后续 Copula 使用 fprintf(Wind Beta: a%.4f, b%.4f\n, a_w, b_w); fprintf(Solar Beta: a%.4f, b%.4f\n, a_s, b_s);这段代码的关键点fitdist是 MATLAB 统计工具箱的函数返回的分布对象可以直接调cdf和pdf。参数a和b是 Beta 分布的两个形状参数控制分布的偏度和峰度。风电出力的 Beta 拟合通常 a 较小、b 较大说明出力偏低的时间居多光伏白天出力的 a 和 b 则取决于当地日照条件。拟合完一定要画 QQ 图或 PP 图检查光看参数数值没有直观判断力。2.3 Gaussian Copula 的参数估计与采样步骤边缘分布拟合好之后下一步是把原始数据变换到均匀分布空间。具体做法是对每个变量做概率积分变换U F_wind(x_wind)V F_solar(x_solar)。变换后 U 和 V 都服从 [0,1] 均匀分布它们之间的相关性就完全由 Copula 函数承载了。Gaussian Copula 的参数就是一个相关系数矩阵。对于二元情况就是一个标量 ρ。估计方法用 Kendall 秩相关系数 τ 换算ρ sin(πτ/2)。这个换算关系是 Gaussian Copula 的解析性质比直接算 Pearson 相关系数更稳健因为秩相关不受边缘分布形式影响。% 概率积分变换 U cdf(wind_beta, wind_pu); V cdf(solar_beta, solar_day); % 注意长度对齐 % 用 Kendall tau 估计 Gaussian Copula 参数 tau corr(U, V, Type, Kendall); rho sin(pi * tau / 2); % 生成 N 个联合场景 N 10000; rng(42); % 固定随机种子保证可复现 Z mvnrnd([0 0], [1 rho; rho 1], N); U_sim normcdf(Z(:,1)); V_sim normcdf(Z(:,2)); % 逆变换回原始空间 wind_sim icdf(wind_beta, U_sim); solar_sim icdf(solar_beta, V_sim);逻辑说明mvnrnd生成二元标准正态样本协方差矩阵里的 ρ 就是 Copula 参数。normcdf把正态样本变成均匀样本这一步等价于从 Gaussian Copula 中采样。最后用icdf做逆变换把均匀样本映射回风电和光伏的原始出力空间。rng(42)是固定随机种子做实验时务必加上否则每次跑出来的场景集合都不一样没法对比。参数方面N 取 10000 是经验值——太少则聚类效果不稳定太多则后续 Kmeans 计算量大。ρ 的值一般在 0.2 到 0.6 之间取决于当地风资源和光照的日内互补程度。如果算出来 ρ 是负数说明风电和光伏在日内呈互补关系白天光照好、风速低这在很多地区是正常现象不用觉得模型出了问题。3. Kmeans 场景削减从一万条曲线到五个典型日3.1 场景削减为什么必须做以及 Kmeans 的适用边界Copula 采样出来一万个场景直接丢进优化模型是不现实的——计算量扛不住。场景削减的目标是用少量典型场景加对应概率尽可能保留原始场景集合的统计特征。Kmeans 是最常用的手段因为它实现简单、收敛快而且聚类中心本身就是一条完整的出力曲线物理意义明确。但 Kmeans 有它的边界。它假设簇是球形的、各簇方差相近对于风光出力这种时序曲线如果不同季节的曲线形态差异很大比如夏季和冬季的日内模式完全不同Kmeans 可能把形态差异大的曲线硬分到同一簇。这时候可以考虑先用层次聚类做初步划分再在每簇内用 Kmeans 细化。不过对于大多数工程应用直接对归一化后的曲线做 Kmeans 已经够用。另一个关键决策是聚类数 K 怎么定。常用的方法有肘部法、轮廓系数法、以及基于场景削减误差的指标。我一般先用肘部法粗选再结合后续优化模型对场景数的容忍度来定。如果优化模型是线性规划K 取 5 到 10 都行如果是混合整数规划K 最好不超过 5否则求解时间会明显拉长。3.2 把 24 小时出力曲线喂给 Kmeans特征构造与归一化Kmeans 的输入是特征矩阵每一行是一个场景每一列是一个特征。对于 24 小时出力曲线最直接的做法是把 24 个时间点的出力值作为 24 个特征。但这样有个问题不同场景的出力幅值差异可能很大Kmeans 基于欧氏距离幅值大的特征会主导距离计算。解决办法是先对每条曲线做归一化让所有曲线在同一尺度下比较形态。归一化方式有两种一是最大最小值归一化到 [0,1]二是 Z-score 标准化。对于风光出力曲线我倾向于用最大最小值归一化因为出力本身有物理上界额定容量归一化后保留了这个信息。% 假设 wind_sim 和 solar_sim 是 N x 24 的矩阵每行一个场景 % 拼接成联合场景矩阵 joint_scenes [wind_sim, solar_sim]; % N x 48 % 对每条曲线做最大最小值归一化 scene_norm zeros(size(joint_scenes)); for i 1:size(joint_scenes, 1) row joint_scenes(i, :); scene_norm(i, :) (row - min(row)) / (max(row) - min(row) eps); end % 或者直接对每一列做归一化保留时间维度的形态 scene_norm (joint_scenes - min(joint_scenes)) ./ ... (max(joint_scenes) - min(joint_scenes) eps);这里有个容易搞混的地方按行归一化和按列归一化效果完全不同。按行归一化是把每条曲线自己的幅值信息抹掉只保留形态按列归一化是保留每个时间点的相对大小。做场景削减通常用按行归一化因为我们关心的是曲线形态的相似性而不是绝对幅值。但如果你后续优化模型需要绝对出力值那聚类完还要把幅值信息还原回去。3.3 聚类数 K 的选取与典型场景概率计算K 的选取没有一刀切的标准。我通常跑一个循环K 从 2 到 10每次记录簇内平方和然后画肘部图。肘部拐点对应的 K 就是候选值。但肘部法有时候拐点不明显这时候结合轮廓系数一起看。K_range 2:10; sse zeros(length(K_range), 1); sil zeros(length(K_range), 1); for k K_range [idx, C, sumd] kmeans(scene_norm, k, Replicates, 5, MaxIter, 500); sse(k-1) sum(sumd); sil(k-1) mean(silhouette(scene_norm, idx)); end % 画肘部图 figure; yyaxis left; plot(K_range, sse, -o); ylabel(SSE); yyaxis right; plot(K_range, sil, -s); ylabel(Silhouette); xlabel(K);Replicates设为 5 表示用 5 个不同的初始中心各跑一次取最优结果避免陷入局部最优。MaxIter是最大迭代次数500 对于 24 维特征通常够用。silhouette计算轮廓系数值越接近 1 说明聚类效果越好。聚类完成后每个簇的典型场景就是聚类中心。该场景的概率等于簇内样本数除以总样本数。注意聚类中心是在归一化空间里的如果要还原到原始出力空间需要用对应簇内样本的最大最小值做逆变换。% 假设最终选定 K5 K 5; [idx, C] kmeans(scene_norm, K, Replicates, 10); % 计算每个典型场景的概率 prob zeros(K, 1); for k 1:K prob(k) sum(idx k) / length(idx); end % 还原到原始空间以风电为例 wind_typical zeros(K, 24); for k 1:K cluster_samples joint_scenes(idx k, 1:24); wind_typical(k, :) mean(cluster_samples, 1); end这里用簇内样本的均值作为典型场景比直接用归一化空间的聚类中心更合理因为均值保留了原始幅值信息。概率计算就是简单的频率统计所有簇的概率之和为 1。4. 避坑与排查风光联合场景生成里最容易翻车的五个地方4.1 概率积分变换后出现 Inf 或 NaN现象调用cdf做概率积分变换时结果里出现 Inf 或 NaN后续 Copula 参数估计直接报错。原因边缘分布拟合的参数不合理导致某些数据点的 CDF 值趋近于 0 或 1取对数或做逆变换时溢出。常见于 Beta 分布的边界值处理——如果原始数据里有 0 或 1Beta 分布的 CDF 在边界处可能返回极端值。解决在拟合前对数据做微小的截断比如把 0 替换为 1e-6把 1 替换为 1-1e-6。另外检查fitdist的拟合结果是否收敛如果参数估计时出现警告换用矩估计法重新拟合。4.2 Kendall 相关系数算出来是负数以为代码写错了现象风电和光伏的 Kendall τ 是负值ρ 也是负的采样出来的场景里风电高的时候光伏低看起来“反了”。原因这不是代码错误而是数据本身的特征。很多地区风光在日内确实呈互补关系——白天光照强、风速低夜间风速高、无光照。负相关是合理的。解决不用改代码但要在论文或报告里解释清楚这个负相关的物理含义。如果确实预期是正相关检查数据对齐是否正确——比如风电和光伏的时间戳是否错位或者单位是否统一。4.3 Kmeans 聚类结果每次跑都不一样现象同样的数据同样的 K每次运行kmeans得到的聚类中心和概率都不同。原因Kmeans 对初始中心敏感默认只跑一次容易陷入局部最优。另外如果没有固定随机种子初始化的随机性也会导致结果波动。解决设置Replicates参数为 5 到 10让算法多跑几次取最优。同时在脚本开头加rng(42)固定随机种子。如果数据量特别大可以用kmeans的Start参数指定初始中心比如用层次聚类的中心作为起点。4.4 典型场景的日出力曲线出现“锯齿”或不平滑现象聚类还原出来的典型场景曲线在某些时间点突然跳变看起来不像真实出力曲线。原因Kmeans 在归一化空间聚类还原时用簇内均值如果簇内样本的形态差异较大均值曲线就会出现不自然的波动。另外如果归一化时用了按列归一化还原后时间维度上的连续性会被破坏。解决改用按行归一化保留每条曲线的时序形态。如果锯齿仍然明显增加 K 值让每个簇更紧凑或者改用 K-medoids 算法它选的是实际样本作为中心曲线更自然。4.5 场景概率之和不为 1 或某个簇概率为 0现象计算完概率发现总和不是 1或者某个典型场景的概率是 0。原因概率为 0 说明该簇没有分配到任何样本通常是 K 设得太大或者初始中心选得不好导致空簇。概率之和不为 1 一般是浮点误差累积或者计算时用了错误的样本总数。解决概率为 0 时减少 K 值或者用kmeans的EmptyAction参数设为singleton让空簇重新分配。概率之和的误差在 1e-10 量级以内可以忽略如果偏差大检查idx的长度是否和样本总数一致。5. 从典型场景到优化模型一个验证场景质量的实用技巧场景生成完怎么判断这组典型场景到底能不能用最直接的办法是回代验证把典型场景和对应概率输入你的优化模型看目标函数值和用全量场景或历史数据算出来的差多少。如果偏差在可接受范围内说明场景削减没有丢失关键信息。但更细一步的验证是看统计特征保留度。具体做法是对原始场景集合和典型场景集合分别计算风电出力的均值、方差、以及风光之间的 Kendall τ然后对比。典型场景的统计量应该和原始场景接近如果偏差大说明 K 选小了或者聚类特征构造有问题。% 原始场景的统计特征 wind_mean_orig mean(wind_sim(:)); wind_std_orig std(wind_sim(:)); tau_orig corr(U_sim, V_sim, Type, Kendall); % 典型场景的统计特征加权 wind_mean_typ sum(prob .* mean(wind_typical, 2)); wind_std_typ sqrt(sum(prob .* var(wind_typical, 0, 2))); % 典型场景的 Kendall tau 需要把典型场景展开成加权样本后计算这个验证步骤花不了几分钟但能帮你避免“场景看起来没问题、一放进优化模型结果就偏”的尴尬。我自己的习惯是每次调整 K 或 Copula 参数后先跑这个验证通过了再往下做。还有一个实用技巧是可视化对比。把原始场景集合画成灰色半透明曲线典型场景画成粗实线叠在一起看。如果典型场景能覆盖原始场景的主要形态区域说明削减合理如果某些区域的原始曲线完全没有被典型场景代表说明 K 不够或者聚类特征需要调整。这个图在论文里也很好用审稿人一眼就能看出场景削减的效果。最后说一个我踩过的坑有一次用 t-Copula 生成场景自由度参数设得太小导致尾部样本过多聚类出来的典型场景里出现了好几个极端场景概率还不低。后来把自由度调到 5 以上尾部就正常了。t-Copula 的自由度参数 ν 控制尾部厚度ν 越小尾部越厚一般取 3 到 10 之间具体看数据的尾部特征。如果不确定先用 Gaussian Copula 跑一遍做基准再换 t-Copula 对比差异明显再考虑用 t-Copula。希望帮到你。本文还有配套的精品资源点击获取 风光联合场景生成：Copula建模与Kmeans削减的MATLAB实战
简介本资源为基于Copula理论与K-means聚类的风光联合场景生成MATLAB代码包面向计算机、电子信息工程、数学等专业的学生及科研人员可用于课程设计、期末大作业与毕业设计中的新能源出力场景建模与聚类分析。压缩包共15个文件约4.77MB包含1个main.m主程序、2个xlsx风光功率数据文件、3份Copula与场景生成相关PDF学习资料、8张png运行结果图及说明文本结构清晰、注释明细支持参数化编程与灵活改参。代码兼容MATLAB 2014、2019a与2024a版本附赠案例数据可直接运行帮助读者快速复现从边缘分布拟合、Copula联合分布构建到K-means场景聚类的完整流程理解变量间复杂依赖关系的建模思路。目前已有122人学习下载适合需要掌握场景生成方法并完成相关课题的中高级学习者参考。1. 风光联合场景生成到底在算什么从 Copula 到 Kmeans 的一条完整链路做新能源出力建模的人迟早会撞上一个问题风电场和光伏电站的出力曲线单独看都还好一旦要放在同一个系统里做联合调度或可靠性评估两条曲线之间的相关性就成了绕不过去的坎。你如果直接把风电和光伏的历史出力序列各自拟合一个边缘分布再独立采样拼在一起得到的场景集合会严重失真——因为现实中太阳落山前后光伏骤降的那段时间风速往往也在变化两者不是各走各的。这就是风光联合场景生成要解决的核心问题在保留每个变量自身统计特性的同时把变量之间的相关结构也复现出来。标题里的 64Copula 指的是一种具体的 Copula 建模方案Kmeans 聚类则负责把生成出来的海量场景压缩成少量典型场景供后续优化模型使用。整条链路是历史出力数据 → 边缘分布拟合 → Copula 相关结构建模 → 采样生成大量场景 → Kmeans 聚类削减 → 输出典型场景及其概率。适合做电力系统规划、微网容量配置、储能优化调度的工程师也适合正在写相关方向论文、需要一套能跑通的 MATLAB 代码做实验的研究生。下面我把这条链路拆开从原理到代码到踩坑一步步讲清楚。2. Copula 建模风光相关性的原理与 MATLAB 实现2.1 为什么不能用协方差矩阵直接描述风光相关性很多人第一反应是用多元正态分布来建模风电和光伏的联合分布协方差矩阵直接描述相关性。这个思路在变量服从正态分布时没问题但风电出力典型分布是 Beta 或 Weibull光伏出力在白天呈 Beta 分布、夜间恒为零两者都不是正态的。硬套多元正态边缘分布就被扭曲了尾部相关性也完全对不上。Copula 的思路是把联合分布拆成两部分边缘分布各自拟合变量之间的依赖结构单独用一个 Copula 函数描述。Sklar 定理保证了这种拆分的数学合法性——任何联合分布都可以写成边缘分布和 Copula 函数的组合。这样一来风电用 Weibull 拟合、光伏用 Beta 拟合两者之间的相关结构用 Gaussian Copula 或 t-Copula 来刻画各管各的互不干扰。选 Gaussian Copula 还是 t-Copula如果风光出力序列在极端天气下表现出明显的尾部依赖比如连续阴天伴随低风速t-Copula 更合适因为它有更厚的尾部。如果数据量不大、尾部特征不明显Gaussian Copula 参数少、估计稳定是更稳妥的起点。我一般先用 Gaussian 跑通再看尾部拟合效果决定要不要换。2.2 边缘分布拟合风电 Weibull、光伏 Beta 的参数估计先看风电。风速通常用 Weibull 分布拟合但注意这里要建模的是风电出力而非风速。如果只有出力数据可以直接对出力做 Beta 拟合或者用风速-功率曲线反推。常见做法是对归一化后的风电出力序列做 Beta 拟合形状参数用最大似然估计。光伏出力白天归一化后同样适合 Beta 分布。但夜间出力为零这部分不能直接丢——丢了会改变日内的时序结构。处理方式有两种一是只对白天时段建模夜间单独标记为零二是对全时段做零膨胀 Beta 拟合。代码里通常采用第一种简单且不容易翻车。% 假设 wind_pu 和 solar_pu 是归一化后的风电和光伏出力列向量 % 去除夜间零值只对光伏白天时段拟合 day_idx solar_pu 0; solar_day solar_pu(day_idx); % 风电 Beta 拟合 wind_beta fitdist(wind_pu, Beta); a_w wind_beta.a; b_w wind_beta.b; % 光伏 Beta 拟合仅白天 solar_beta fitdist(solar_day, Beta); a_s solar_beta.a; b_s solar_beta.b; % 输出参数供后续 Copula 使用 fprintf(Wind Beta: a%.4f, b%.4f\n, a_w, b_w); fprintf(Solar Beta: a%.4f, b%.4f\n, a_s, b_s);这段代码的关键点fitdist是 MATLAB 统计工具箱的函数返回的分布对象可以直接调cdf和pdf。参数a和b是 Beta 分布的两个形状参数控制分布的偏度和峰度。风电出力的 Beta 拟合通常 a 较小、b 较大说明出力偏低的时间居多光伏白天出力的 a 和 b 则取决于当地日照条件。拟合完一定要画 QQ 图或 PP 图检查光看参数数值没有直观判断力。2.3 Gaussian Copula 的参数估计与采样步骤边缘分布拟合好之后下一步是把原始数据变换到均匀分布空间。具体做法是对每个变量做概率积分变换U F_wind(x_wind)V F_solar(x_solar)。变换后 U 和 V 都服从 [0,1] 均匀分布它们之间的相关性就完全由 Copula 函数承载了。Gaussian Copula 的参数就是一个相关系数矩阵。对于二元情况就是一个标量 ρ。估计方法用 Kendall 秩相关系数 τ 换算ρ sin(πτ/2)。这个换算关系是 Gaussian Copula 的解析性质比直接算 Pearson 相关系数更稳健因为秩相关不受边缘分布形式影响。% 概率积分变换 U cdf(wind_beta, wind_pu); V cdf(solar_beta, solar_day); % 注意长度对齐 % 用 Kendall tau 估计 Gaussian Copula 参数 tau corr(U, V, Type, Kendall); rho sin(pi * tau / 2); % 生成 N 个联合场景 N 10000; rng(42); % 固定随机种子保证可复现 Z mvnrnd([0 0], [1 rho; rho 1], N); U_sim normcdf(Z(:,1)); V_sim normcdf(Z(:,2)); % 逆变换回原始空间 wind_sim icdf(wind_beta, U_sim); solar_sim icdf(solar_beta, V_sim);逻辑说明mvnrnd生成二元标准正态样本协方差矩阵里的 ρ 就是 Copula 参数。normcdf把正态样本变成均匀样本这一步等价于从 Gaussian Copula 中采样。最后用icdf做逆变换把均匀样本映射回风电和光伏的原始出力空间。rng(42)是固定随机种子做实验时务必加上否则每次跑出来的场景集合都不一样没法对比。参数方面N 取 10000 是经验值——太少则聚类效果不稳定太多则后续 Kmeans 计算量大。ρ 的值一般在 0.2 到 0.6 之间取决于当地风资源和光照的日内互补程度。如果算出来 ρ 是负数说明风电和光伏在日内呈互补关系白天光照好、风速低这在很多地区是正常现象不用觉得模型出了问题。3. Kmeans 场景削减从一万条曲线到五个典型日3.1 场景削减为什么必须做以及 Kmeans 的适用边界Copula 采样出来一万个场景直接丢进优化模型是不现实的——计算量扛不住。场景削减的目标是用少量典型场景加对应概率尽可能保留原始场景集合的统计特征。Kmeans 是最常用的手段因为它实现简单、收敛快而且聚类中心本身就是一条完整的出力曲线物理意义明确。但 Kmeans 有它的边界。它假设簇是球形的、各簇方差相近对于风光出力这种时序曲线如果不同季节的曲线形态差异很大比如夏季和冬季的日内模式完全不同Kmeans 可能把形态差异大的曲线硬分到同一簇。这时候可以考虑先用层次聚类做初步划分再在每簇内用 Kmeans 细化。不过对于大多数工程应用直接对归一化后的曲线做 Kmeans 已经够用。另一个关键决策是聚类数 K 怎么定。常用的方法有肘部法、轮廓系数法、以及基于场景削减误差的指标。我一般先用肘部法粗选再结合后续优化模型对场景数的容忍度来定。如果优化模型是线性规划K 取 5 到 10 都行如果是混合整数规划K 最好不超过 5否则求解时间会明显拉长。3.2 把 24 小时出力曲线喂给 Kmeans特征构造与归一化Kmeans 的输入是特征矩阵每一行是一个场景每一列是一个特征。对于 24 小时出力曲线最直接的做法是把 24 个时间点的出力值作为 24 个特征。但这样有个问题不同场景的出力幅值差异可能很大Kmeans 基于欧氏距离幅值大的特征会主导距离计算。解决办法是先对每条曲线做归一化让所有曲线在同一尺度下比较形态。归一化方式有两种一是最大最小值归一化到 [0,1]二是 Z-score 标准化。对于风光出力曲线我倾向于用最大最小值归一化因为出力本身有物理上界额定容量归一化后保留了这个信息。% 假设 wind_sim 和 solar_sim 是 N x 24 的矩阵每行一个场景 % 拼接成联合场景矩阵 joint_scenes [wind_sim, solar_sim]; % N x 48 % 对每条曲线做最大最小值归一化 scene_norm zeros(size(joint_scenes)); for i 1:size(joint_scenes, 1) row joint_scenes(i, :); scene_norm(i, :) (row - min(row)) / (max(row) - min(row) eps); end % 或者直接对每一列做归一化保留时间维度的形态 scene_norm (joint_scenes - min(joint_scenes)) ./ ... (max(joint_scenes) - min(joint_scenes) eps);这里有个容易搞混的地方按行归一化和按列归一化效果完全不同。按行归一化是把每条曲线自己的幅值信息抹掉只保留形态按列归一化是保留每个时间点的相对大小。做场景削减通常用按行归一化因为我们关心的是曲线形态的相似性而不是绝对幅值。但如果你后续优化模型需要绝对出力值那聚类完还要把幅值信息还原回去。3.3 聚类数 K 的选取与典型场景概率计算K 的选取没有一刀切的标准。我通常跑一个循环K 从 2 到 10每次记录簇内平方和然后画肘部图。肘部拐点对应的 K 就是候选值。但肘部法有时候拐点不明显这时候结合轮廓系数一起看。K_range 2:10; sse zeros(length(K_range), 1); sil zeros(length(K_range), 1); for k K_range [idx, C, sumd] kmeans(scene_norm, k, Replicates, 5, MaxIter, 500); sse(k-1) sum(sumd); sil(k-1) mean(silhouette(scene_norm, idx)); end % 画肘部图 figure; yyaxis left; plot(K_range, sse, -o); ylabel(SSE); yyaxis right; plot(K_range, sil, -s); ylabel(Silhouette); xlabel(K);Replicates设为 5 表示用 5 个不同的初始中心各跑一次取最优结果避免陷入局部最优。MaxIter是最大迭代次数500 对于 24 维特征通常够用。silhouette计算轮廓系数值越接近 1 说明聚类效果越好。聚类完成后每个簇的典型场景就是聚类中心。该场景的概率等于簇内样本数除以总样本数。注意聚类中心是在归一化空间里的如果要还原到原始出力空间需要用对应簇内样本的最大最小值做逆变换。% 假设最终选定 K5 K 5; [idx, C] kmeans(scene_norm, K, Replicates, 10); % 计算每个典型场景的概率 prob zeros(K, 1); for k 1:K prob(k) sum(idx k) / length(idx); end % 还原到原始空间以风电为例 wind_typical zeros(K, 24); for k 1:K cluster_samples joint_scenes(idx k, 1:24); wind_typical(k, :) mean(cluster_samples, 1); end这里用簇内样本的均值作为典型场景比直接用归一化空间的聚类中心更合理因为均值保留了原始幅值信息。概率计算就是简单的频率统计所有簇的概率之和为 1。4. 避坑与排查风光联合场景生成里最容易翻车的五个地方4.1 概率积分变换后出现 Inf 或 NaN现象调用cdf做概率积分变换时结果里出现 Inf 或 NaN后续 Copula 参数估计直接报错。原因边缘分布拟合的参数不合理导致某些数据点的 CDF 值趋近于 0 或 1取对数或做逆变换时溢出。常见于 Beta 分布的边界值处理——如果原始数据里有 0 或 1Beta 分布的 CDF 在边界处可能返回极端值。解决在拟合前对数据做微小的截断比如把 0 替换为 1e-6把 1 替换为 1-1e-6。另外检查fitdist的拟合结果是否收敛如果参数估计时出现警告换用矩估计法重新拟合。4.2 Kendall 相关系数算出来是负数以为代码写错了现象风电和光伏的 Kendall τ 是负值ρ 也是负的采样出来的场景里风电高的时候光伏低看起来“反了”。原因这不是代码错误而是数据本身的特征。很多地区风光在日内确实呈互补关系——白天光照强、风速低夜间风速高、无光照。负相关是合理的。解决不用改代码但要在论文或报告里解释清楚这个负相关的物理含义。如果确实预期是正相关检查数据对齐是否正确——比如风电和光伏的时间戳是否错位或者单位是否统一。4.3 Kmeans 聚类结果每次跑都不一样现象同样的数据同样的 K每次运行kmeans得到的聚类中心和概率都不同。原因Kmeans 对初始中心敏感默认只跑一次容易陷入局部最优。另外如果没有固定随机种子初始化的随机性也会导致结果波动。解决设置Replicates参数为 5 到 10让算法多跑几次取最优。同时在脚本开头加rng(42)固定随机种子。如果数据量特别大可以用kmeans的Start参数指定初始中心比如用层次聚类的中心作为起点。4.4 典型场景的日出力曲线出现“锯齿”或不平滑现象聚类还原出来的典型场景曲线在某些时间点突然跳变看起来不像真实出力曲线。原因Kmeans 在归一化空间聚类还原时用簇内均值如果簇内样本的形态差异较大均值曲线就会出现不自然的波动。另外如果归一化时用了按列归一化还原后时间维度上的连续性会被破坏。解决改用按行归一化保留每条曲线的时序形态。如果锯齿仍然明显增加 K 值让每个簇更紧凑或者改用 K-medoids 算法它选的是实际样本作为中心曲线更自然。4.5 场景概率之和不为 1 或某个簇概率为 0现象计算完概率发现总和不是 1或者某个典型场景的概率是 0。原因概率为 0 说明该簇没有分配到任何样本通常是 K 设得太大或者初始中心选得不好导致空簇。概率之和不为 1 一般是浮点误差累积或者计算时用了错误的样本总数。解决概率为 0 时减少 K 值或者用kmeans的EmptyAction参数设为singleton让空簇重新分配。概率之和的误差在 1e-10 量级以内可以忽略如果偏差大检查idx的长度是否和样本总数一致。5. 从典型场景到优化模型一个验证场景质量的实用技巧场景生成完怎么判断这组典型场景到底能不能用最直接的办法是回代验证把典型场景和对应概率输入你的优化模型看目标函数值和用全量场景或历史数据算出来的差多少。如果偏差在可接受范围内说明场景削减没有丢失关键信息。但更细一步的验证是看统计特征保留度。具体做法是对原始场景集合和典型场景集合分别计算风电出力的均值、方差、以及风光之间的 Kendall τ然后对比。典型场景的统计量应该和原始场景接近如果偏差大说明 K 选小了或者聚类特征构造有问题。% 原始场景的统计特征 wind_mean_orig mean(wind_sim(:)); wind_std_orig std(wind_sim(:)); tau_orig corr(U_sim, V_sim, Type, Kendall); % 典型场景的统计特征加权 wind_mean_typ sum(prob .* mean(wind_typical, 2)); wind_std_typ sqrt(sum(prob .* var(wind_typical, 0, 2))); % 典型场景的 Kendall tau 需要把典型场景展开成加权样本后计算这个验证步骤花不了几分钟但能帮你避免“场景看起来没问题、一放进优化模型结果就偏”的尴尬。我自己的习惯是每次调整 K 或 Copula 参数后先跑这个验证通过了再往下做。还有一个实用技巧是可视化对比。把原始场景集合画成灰色半透明曲线典型场景画成粗实线叠在一起看。如果典型场景能覆盖原始场景的主要形态区域说明削减合理如果某些区域的原始曲线完全没有被典型场景代表说明 K 不够或者聚类特征需要调整。这个图在论文里也很好用审稿人一眼就能看出场景削减的效果。最后说一个我踩过的坑有一次用 t-Copula 生成场景自由度参数设得太小导致尾部样本过多聚类出来的典型场景里出现了好几个极端场景概率还不低。后来把自由度调到 5 以上尾部就正常了。t-Copula 的自由度参数 ν 控制尾部厚度ν 越小尾部越厚一般取 3 到 10 之间具体看数据的尾部特征。如果不确定先用 Gaussian Copula 跑一遍做基准再换 t-Copula 对比差异明显再考虑用 t-Copula。希望帮到你。本文还有配套的精品资源点击获取