简介一套基于64种Copula函数与K-means聚类的风光联合场景生成Matlab代码面向电力系统、数据分析等方向的研究者以及进行课程设计或毕业设计的学生。其核心是利用Copula刻画风光出力变量间的依赖结构结合核密度估计生成联合场景再通过K-means聚类提取典型场景为微网优化配置提供输入。资源包共15个文件约4.77MB包含m主源码、xlsx风功率与光伏实测数据、pdf理论资料与案例论文、png结果示意图及txt说明文件结构清晰。代码采用参数化编程支持matlab2014/2019a/2024a注释详细附赠可直接运行的数据方便根据实际需求调整参数并复现实验。已有122人学习适合需要掌握Copula场景生成、K-means聚类分析或开展新能源系统优化研究的人参考。1. 风光联合场景生成为什么“单独建模”会让你的电网规划偏乐观做新能源消纳评估、储能容量配置或者微网日前调度的人迟早会撞上同一个问题光和风不出力不是独立事件。傍晚光伏衰减时风往往起来半夜风大而负荷低这些“联合出力状态”单独看风电曲线、光伏曲线是看不见的。风光联合场景生成就是为了把这层相关性装进模型里给后续优化提供一组“可能同时发生的出力组合”。常见的成熟路线是用 Copula 构造联合分布先采样出一大批初使场景再用 KMeans 聚类削成少数几个典型场景供随机优化使用。我最早做配电网随机潮流时也迷信过独立建模结果算出来的线路负载率明显偏乐观后来切成这套流程才把风险暴露出来。这篇就把 Copula 拟合、场景采样、KMeans 削减的完整走法和踩过的坑一次说清。2. 从历史出力到 Copula 拟合把风光的“相关性脾气”封装进一个模型2.1 为什么非要用 Copula线性相关系数骗了你风电和光伏出力之间的关系并不是一条直线。Pearson 相关系数只能描述线性相关碰上“小出力段相对独立、大出力段同步爬坡”这类尾部相关算出来的系数常常接近 0导致你误以为两者不相关。Copula 的好处是把“边缘分布”和“相关结构”拆开边缘分布负责刻画单变量的“出力形状”很多小时为 0、午间尖峰相关结构单独用一个 Copula 函数描述。Copula 本身就是一个均匀分布空间上的联合分布函数通过它可以把两个变量的分位数对应关系完整保留下来不需要假设线性。我一般先算一下 Kendall 秩相关系数 tau而不是直接看 Pearson。tau 对单调变换不敏感风电出力的非线性爬坡特性不会把它带偏。选 Copula 族时常见做法是先用copulastat把 Frank、Clayton、Gumbel 的理论 tau 算出来再和目标 tau 对比如果尾部相关性明显比如大风时段光也弱Clayton 或 Gumbel 会比 Frank 更像。不确定就多拟合几个族用 AIC 挑。2.2 边缘分布与 Copula 分离估计流程的两步走这一步的顺序不能乱。先用历史样本拟合风电、光伏各自的边缘分布把原始出力值转换为均匀分布变量 U1、U2然后再用这两个 U 去估计 Copula 参数。如果跳过了边缘分布直接拿原始数据拟合 Copula参数会失真后续采样逆变换也会走样。边缘分布常见有两种做法参数化Weibull 拟合风速再转功率和非参数化直接用ksdensity做核密度估计。做场景生成时我更推荐核密度估计它不需要假设具体分布形状对“大量 0 出力”这种非光滑分布更宽容。注意 KDE 得到的 CDF 值两端会非常接近 0 和 1直接取 log 会出 inf所以要加一个很小的裁剪。2.3 Matlab 里拟合并检验copulafit 与 ksdensity 的最小可运行代码下面给出一段能直接跑通的流程骨架。假设data是 N×2 的历史出力矩阵第一列风电、第二列光伏取值范围 0 到装机容量标幺化后 0~1。% 读入数据做标幺化除以装机容量 data data ./ max(data); % 或除以各自的额定容量 wind data(:, 1); solar data(:, 2); % 1) 核密度估计边缘分布得到 CDF 值 u_wind ksdensity(wind, wind, function, cdf); u_solar ksdensity(solar, solar, function, cdf); % 2) 裁剪到 (0,1) 开区间避免 log(0) 或 log(1) eps_val 1e-4; u_wind min(max(u_wind, eps_val), 1 - eps_val); u_solar min(max(u_solar, eps_val), 1 - eps_val); u [u_wind, u_solar]; % 3) 拟合 Frank / Clayton / Gumbel选 AIC 最小的 rho corr(u, Type, Kendall); % 目标 tau [theta_f, ~] copulafit(Frank, u); [theta_c, ~] copulafit(Clayton, u); [theta_g, ~] copulafit(Gumbel, u); % 4) 回验把拟合的 theta 换算成理论 tau和目标对比 tau_f copulastat(Frank, theta_f); tau_c copulastat(Clayton, theta_c); tau_g copulastat(Gumbel, theta_g); fprintf(目标 tau: %.3f\nFrank: %.3f\nClayton: %.3f\nGumbel: %.3f\n, ... rho(1,2), tau_f, tau_c, tau_g);逻辑说明ksdensity的function,cdf返回的是每个样本点对应的经验 CDF 值这一步把原始出力映射到均匀分布。copulafit的第一个参数指定 Copula 族返回该族的参数theta。copulastat则把参数换算成理论 Kendall tau用来检验拟合是否合理。需要特别注意eps_val的裁剪不裁剪的话copulafit内部涉及 log 运算时会出现 inf拟合直接失败。实际项目里我一般会用 AIC 来选族而不是肉眼比 tau% 负对数似然越小越好 [~, nlogl_f] copulafit(Frank, u); [~, nlogl_c] copulafit(Clayton, u); [~, nlogl_g] copulafit(Gumbel, u); aic [2*nlogl_f, 2*nlogl_c, 2*nlogl_g]; [~, best] min(aic); fprintf(最优 Copula: %d (1Frank, 2Clayton, 3Gumbel)\n, best);AIC 比较的好处是它同时惩罚了参数个数避免你为了凑相关性多带参数。很多人到这里就停下来了但别忘了边缘分布本身也参与场景生成后面逆变换要靠它所以边缘分布的质量直接决定生成场景的物理真实性。3. 生成联合场景采样、逆变换与初始场景数的设定3.1 多时段场景怎么组织上一章拟合的是一个静态的“风-光二维联合分布”但实际调度往往需要全天的时序场景。常见做法有三种第一种是忽略时段相关性把每天 24 小时的风光出力拉成一个 48 维向量在第 2 章那种二维 Copula 上按小时各算一遍然后拼接。第二种是先对“日平均出力”做 Copula 采样再按历史条件分配日内形状。第三种是逐小时建立 24 个二维 Copula采样时每个小时独立取一组适用于不考虑自相关的场景。我一般用第一种原因很简单工程上足够用而且第四步 KMeans 聚类正好可以把时段之间的结构一并处理掉。聚类不是只砍数量它还会把“像的日子”合并成一个代表所以只要你生成的初始场景覆盖到了关键相关模式时序上的粗糙可以接受。第三种方法如果每天 24 小时都用独立采样可能出现风电白天猛涨、光伏夜里不为零这种物理上不合理的组合必须在采样后加约束过滤。3.2 copularnd 采样与逆变换还原出力有了第 2 章拟合好的 Copula 参数下一步就是用copularnd在均匀空间生成 N 个样本对再把每个 U 分量通过边缘分布的逆 CDF 还原成出力值。这里的逆 CDF 仍然用ksdensity的function,icdf来实现。% 用选定的 Copula 族生成 N 个均匀空间样本 N 512; % 初始场景数 u_sim copularnd(Frank, theta_f, N); % 逆变换从均匀空间还原到出力空间 wind_sim ksdensity(wind, u_sim(:, 1), function, icdf); solar_sim ksdensity(solar, u_sim(:, 2), function, icdf); % 检查还原后的取值范围 fprintf(风电生成范围: [%.3f, %.3f]\n, min(wind_sim), max(wind_sim)); fprintf(光伏生成范围: [%.3f, %.3f]\n, min(solar_sim), max(solar_sim));逻辑说明copularnd(Frank, theta_f, N)生成 N×2 的矩阵每行是 [0,1] 上均匀分布的样本对并且它们之间带有 Frank Copula 的相关结构。ksdensity(...,function,icdf)是对 KDE 分布做逆变换输入是概率值输出是对应的原始出力。这一步相当于把“相关性已经装好”的均匀随机数翻译回物理量纲。参数说明N越大越好但后续聚类和优化会变慢W 通常取 200~1000标题里的 64 属于偏小的设置样本太少聚类中心的稳定性会差。3.3 “64Copula”的含义与场景数怎么定标题里出现的“64”常见有两种理解一种是初始场景数为 64另一种是生成 64 个典型日。我猜你的压缩包里多半是前者——64 个 Copula 采样场景再聚类成更少的代表。但以我的经验64 作为初始场景太少了。KMeans 聚类需要一个足够稠密的样本空间来保证每个簇都有足够的形状支撑64 个点散布在 48 维空间里会非常稀疏簇中心容易落在样本缝隙里。我一般至少生成 256 个想追求稳健就 512 个。聚类出来的典型场景数量反而比较小通常是 3 到 5 个因为随机优化里每多一个场景计算量就成倍增长。如果你拿到一份代码先确认它的 64 到底是初始采样数还是聚类后的目标数。如果是前者直接把它调大到 256 再跑结果往往会明显改善如果是后者那初始样本数可能已经是 512 或 1024这个参数设置是合理的。4. KMeans 聚类削减从几百个场景到 3~5 个典型场景4.1 场景削减解决的是“计算不可行”随机优化的计算复杂度随场景数量近似线性增长但电力系统的优化模型里带有整数变量机组启停、储能充放电状态场景一多混合整数规划直接算不动。场景削减就是牺牲一点分布精度换计算可行性。KMeans 在这里做的是把 N 个初始场景划分成 K 个簇每个簇的中心就是一个“典型场景”簇内样本的比例就是该场景在优化模型中出现的概率。这样原来 512 个场景的期望值计算就变成了 4 个场景的加权求和。4.2 特征工程归一化、距离度量和特征拼接KMeans 聚类前每个场景要变成一个向量。最简单的是 [风电小时序列(24), 光伏小时序列(24)]即 48 维。但直接拿原始出力向量去算欧氏距离有个问题如果数据是标幺值风电出力范围 0~1 而光伏也是 0~1尺度一致还好但如果你用了实际功率单位MW风电装机 100MW 和光伏 50MW 的数值范围不同欧氏距离会被风电主导。我一般先把风电和光伏分别按各自的装机容量归一化再拼接或者干脆各自单独做标准化。另外如果只关心“日总出力”水平也可以换成特征向量 [风电日电量, 光伏日电量]聚类结果更稳定但会丢失日内爬坡形状。KMeans 默认用平方欧氏距离这个距离对“同时为 0”的相似度非常敏感。如果历史数据里阴雨天和小风日都能让风、光同时落在 0 附近它们在欧氏距离下会合到同一个簇里但实际这两种场景的后续调度策略完全不同一个可能靠储能顶一个可能靠外购电。这种情况我会考虑把距离改成基于形状的度量或者把“出力为 0 的时段数”作为一个额外特征拼进去。4.3 KMeans 聚类的 Matlab 代码与 K 值选择Matlab 自带的kmeans函数直接用即可默认用 kmeans 初始化比随机初始化稳定得多。下面给出一段完整的聚类代码% 构造场景矩阵N_scenes x 48前24列风电后24列光伏 X [wind_sim_mat, solar_sim_mat]; % 每行是一个初始场景 % 可选按列标准化消除量纲影响 X_std (X - mean(X)) ./ std(X); % z-score 标准化 % 跑 KMeansReplicates5 降低随机初始化影响 K 4; rng(42); % 固定随机种子保证可复现 [idx, C, sumd] kmeans(X_std, K, Replicates, 5, MaxIter, 1000); % 统计每个簇的样本数换算成概率 counts histcounts(idx, [1:K1]); probs counts / sum(counts); % 还原聚类中心到原始出力空间把标准化逆回去 C_orig C .* std(X) mean(X);逻辑说明X_std是按列做 z-score 标准化避免风电、光伏量纲不一影响距离。kmeans返回的idx是每个场景的簇标签C是簇中心在标准化空间sumd是每个点到中心距离的平方和。counts算出的probs就是典型场景在随机优化里的权重。特别注意的是rng(42)必须放在kmeans之前否则不固定种子每次运行结果都不一样。K 值选择没有绝对标准。我常用的办法是画“肘部图”横轴 K纵轴总簇内离差平方和SSE找斜率骤降的拐点sse zeros(10, 1); for k 1:10 [~, ~, sumd_k] kmeans(X_std, k, Replicates, 3, MaxIter, 1000); sse(k) sum(sumd_k); end plot(1:10, sse, -o); xlabel(K); ylabel(SSE);如果拐点不明显就换轮廓系数或者直接用业务逻辑来定比如你最多能接受 5 个场景进优化模型那 K 就取 5再看聚类中心是否合理。这类“先定计算约束再选 K”的做法比纯看指标更实用。4.4 聚类后概率分配与结果输出聚类完成后典型场景加上对应的概率输出成一个矩阵或表格供优化模型直接引用典型场景编号风电出力向量24h光伏出力向量24h概率1C_orig(1, 1:24)C_orig(1, 25:48)0.312C_orig(2, 1:24)C_orig(2, 25:48)0.273C_orig(3, 1:24)C_orig(3, 25:48)0.244C_orig(4, 1:24)C_orig(4, 25:48)0.18有个细节KMeans 聚类中心是簇内样本的均值也就意味着典型场景是“平均曲线”它天然会抹掉极端天气。如果你关心的是低概率高风险的场景比如连续阴雨加无风纯 KMeans 会把它跟糙的簇混合掉。常见做法是聚类之后单独保留每个簇内离中心最远的样本或 5% 分位数样本作为“保守典型场景”用来做鲁棒校验。5. 避坑指南Copula 拟合与 KMeans 聚类最容易翻车的 5 个点5.1 0 和 1 的边界ksdensity 算出 0 或 1log 直接 inf现象运行copulafit报错提示输入数据包含 0 或 1或者nlogl出现 NaN。原因ksdensity的 CDF 在两个端点会自然收敛到 0 和 1。Copula 的似然函数里对 U 做对数变换log(0) 就是负无穷优化直接崩。解决在copulafit之前把数据裁剪到安全区间。我在工程里一般用u min(max(u, 1e-4), 1 - 1e-4)裁剪阈值建议别低于 1e-6否则 log 值太大数值稳定性差。裁剪后再检查一次sum(u 1 | u 0)确保没有残留。5.2 单参数 Copula 拟合失败初值与边界怎么给现象copulafit(Clayton, u)返回的 theta 为负或者迭代不收敛。原因copulafit对 Archimedean Copula 的估计用的是单参数 MLE参数空间本身有限制Clayton 的 theta 必须大于 -1 且不为 0Gumbel 必须大于 1。当数据相关性非常弱或者样本量太小时MLE 可能迭代到边界返回一个不合理的值。解决不要直接相信默认输出。我一般在拟合前计算 Kendall tau然后用 tau 与参数的理论关系反解一个初值再交给copulafit。比如 Frank Copula 的 theta 可以从 tau 数值解。如果这样还是不行就用网格搜索在合法参数范围内遍历 theta取负对数似然最小的结果。兜底方案是把 Copula 族换成双参数的 t-Copula自由度参数可以额外吸收尾部特性的差异。5.3 KMeans 的 K 值玄学肘部法判读有歧义现象画出来的 SSE 曲线没有明显拐点K3 和 K5 都说得通聚类结果差异却很大。原因48 维空间里 SSE 的下降比较平滑肘部被高维稀释了。真实的风光数据往往存在多级聚类结构单看一个指标很难定。解决改成两步走——先用轮廓系数筛出候选 K 值比如 2 到 8 都算然后逐个 K 生成典型场景喂给下游优化模型看优化结果比如总成本、切负荷率在哪个 K 附近趋于稳定。优化结果不再变化时就是对你这个问题最合适的 K。这比纯统计指标可靠得多。5.4 聚类的距离度量大量 0 出力让欧氏距离失真现象聚类中心出来的典型场景里总有一个簇的风、光出力同时极小另一个簇同时很大中间过渡的簇几乎不存在。原因欧氏距离在稀疏向量上会把“同时为 0”当成高度相似。零出力时段越多这种失真越严重。这是我做西北某风电场项目时踩过的实坑——聚类结果把“阴天无风”和“深夜无光”混进了一个簇。解决拼接特征时加入一段描述“零出力模式”的辅助特征比如每个小时的 0/1 标记。或者换用动态时间规整DTW距离Matlab 里没有现成轮子自己写一个也就二十行。更简单的是直接对数据做变换先对每个时段的出力开根号或取 log1p把 0 附近的距离拉开。5.5 多次运行结果不一致随机种子与复现现象同一份代码今天跑出的典型场景概率是 0.31/0.27/0.24/0.18明天跑变成 0.35/0.25/0.22/0.18。原因copularnd内部要用随机数kmeans的初始质心也是随机的。没固定种子时每次结果都会有波动。解决在代码最前面统一加一句rng(42)位置必须放在所有随机函数调用之前。我在出报告前还会用一批固定种子42、7、2024各跑一遍确认典型场景差异在一个小范围内。如果不同种子结果差很多说明聚类结构本身不稳要回到 4.2 节调整特征。6. 验证生成的场景三个必做检验与一个时序进阶技巧6.1 相关性回验生成的场景有没有把“脾气”保留下来生成 512 个场景后第一件事是算生成样本的 Kendall tau和原始数据的 tau 对比。常用尺子是偏差不超过 0.05。如果偏差大多半是eps_val裁剪过了头把尾部相关性削掉了或者 Copula 族选错了。这一步一分钟就能跑完别跳过。6.2 分布回验KS 测试与 Q-Q 图对每个单变量把生成场景的 CDF 和原始数据的 KDE 做 KS 检验p 值低于 0.05 就说明边缘分布被 Copula 采样带偏了。常见原因是在逆变换时用了不同的ksdensity带宽或者原始数据标幺化方式不统一。Q-Q 图画出来如果两头翘说明 KDE 尾巴不够重可以换用 t 分布拟合边缘。6.3 典型场景的物理合理性出力包络线把生成场景和典型场景画在同一个图上逐一检查光伏场景的夜间出力是否严格为 0风电场景的爬坡率是否超过物理极限比如 15 分钟内出力变化超过装机容量的一半就要怀疑数据质量。机器不会对物理约束负责但你要在场景生成阶段就把这些约束卡掉。6.4 进阶把独立同分布场景扩展为带自相关的时序场景如果你要做的是 24 小时之内的滚动调度纯独立采样会导致相邻时刻的出力突变无法被储能充电功率限制消化。一个实用的方案是先生成 24 小时的随机数序列再用经验 CDF 把序列的自相关结构调整为目标值。Matlab 里可以用fmincon拟合一个 AR(1) 模型的滞后相关系数再映射到[0,1]均匀空间做 Copula 采样。更简单粗暴的做法是生成大量独立场景后用一个滚动窗口平均滤波再重新排序把自相关“涂抹”进去。不完美但计算量便宜工程上够用。这套“Copula 采样生成 KMeans 聚类削减”的流程我从配电网规划做到微网经济调度超过三年的项目里一直用它做不确定性的入口。最深的体会是Copula 和 KMeans 都不是什么高深算法真正的门槛在数据的边界处理、物理合理性校验和 K 值的业务化选择上。我踩过的这些坑——0/1边界、初值不收敛、零出力距离失真——你大概率也会遇到希望这篇能让你少绕几圈。希望帮到你。本文还有配套的精品资源点击获取