前几年做新能源随机场景时我踩过一个特别深的坑风电场景单独生成了一套光伏场景单独生成了一套最后按蒙特卡洛的索引硬拼在一起。当时觉得没什么直到拿实际运行数据复盘发现模型预测的“风光同时低出力”日子比真实情况少了一半备用容量安排得非常被动。原因后来很清楚——风光出力不是独立变量它们在同一片天气系统下会同步变化这种相关结构如果不在生成阶段保住后面做再多优化都是白搭。后来我把整套流程改成基于Copula理论构造联合分布再用K-means做场景削减才算把这个问题理顺。这篇文章就把这套方法的完整思路、实操步骤和踩坑记录摊开讲。核心是两件事第一怎么用Copula把风电、光伏的边缘分布“捏”成一个尊重相关性的联合分布从而生成大批量初选场景第二怎么用K-means把这些成千上万的场景削减成几十个带概率权重的代表场景让下游的随机调度、可靠性评估跑得动。适合电力系统方向的研究生、刚入行的算法工程师以及任何需要给随机优化提供风光联合场景的开发者。整套流程并不神秘数学门槛也不算高只要理解“边缘分布”和“相依结构”可以分开建模后面就顺了。我会从问题本质一直讲到代码骨架再把我实际遇到的几个坑放在最后希望能帮你少折腾几个晚上。1. 风光出力不是“11”联合建模到底解决什么问题1.1 分开生成再拼接为什么总是对不上很多初学场景生成的人第一步都会先做一个“独立采样”风电用一组时序模型生成场景光伏用另一组模型生成场景最后按场景编号一对一对齐。这样做的隐含假设是风、光出力条件独立。但实际气象过程里风电和光伏共享同一个大气环流背景冷锋过境时风速加大同时云量增多导致辐照下降夏季午后局地热对流可能同时带来大风和云层遮挡连续阴雨天则可能让风电和光伏同时处于低出力区间。这些现象都说明二者之间存在复杂且非线性的相依关系。如果无视这种关系把独立生成的两批场景任意组合生成的联合分布与真实分布会有明显偏差。最典型的问题就是“同时低出力”和“同时高出力”这类极值事件的频率失真。我的经验是独立拼接出来的场景往往让调度结果偏乐观系统看起来有很多冗余电源真到极端天气来的时候却不够用。1.2 相关系数描述不了“一起糟糕”的日子有人可能会说要反映相关性算一个皮尔逊相关系数不就行了不行。皮尔逊相关系数衡量的是线性相关而风光出力之间经常不是线性关系低出力区间的收敛、高出力区间的发散、局部的厚尾行为都不是一个标量能刻画的。更重要的是皮尔逊相关系数对极端值敏感样本里几个异常点就能把系数拉得很高或很低。这时候需要用秩相关来辅助判断。Kendall tau和Spearman rho不关心变量的绝对数值只关心排序的一致性对非线性单调变换是稳健的。但即使知道了秩相关也仍然不足以完整描述联合分布——两只变量的尾巴到底是一起翘还是分开翘秩相关回答不了。Copula理论解决的正是这个问题先把每个变量的边缘分布单独建出来再用一个Copula函数把它们的相依结构完整地连接起来相当于把“每个变量自己的脾性”和“变量之间的合拍程度”分开了各建各的最后再拼在一起。1.3 两阶段思路先扩充场景再压缩场景有了Copula之后还需要回答一个工程问题下游优化到底需要多少场景随机规划里几十个场景通常就是极限几百个就会让求解时间不可接受但统计学上又需要足够多样本来覆盖可能的出力组合。这是一个天然矛盾。所以标准做法是两步走。第一步用Copula从联合分布里抽出大量初选场景比如5000个这一步目的是保证样本对真实相关结构的覆盖面足够广第二步用K-means把5000个场景聚成K类每类取聚类中心作为代表场景并用类内样本占比作为概率权重。这样既保留了联合分布里的主要特征又把场景数量压缩到优化问题能承受的范围。后面讲的所有内容都是围绕这两步展开的。2. Copula选型和拟合我把Sklar定理拆成了可落地的四步2.1 边缘分布先别套正态尊重数据真实形状Copula建模的第一步不是选Copula而是先搞定边缘分布。Sklar定理说设风电出力W和光伏出力S的联合分布函数为F(w,s)那么一定存在一个Copula函数C使得F(w,s) C(F_w(w), F_s(s))其中F_w和F_s就是两个边缘分布函数。这个等式最大的价值在于建模时可以独立估计边缘分布然后单独拟合相依结构。实操中第一步就是把历史出力的原始数值转换成均匀分布U即 u F_w(w)v F_s(s)。边缘分布怎么选很多教材喜欢用Weibull拟合风速、Beta拟合光伏但到了出力层面数据特征往往比教材复杂光伏日平均出力可能带有零点、风电出力的分布可能厚尾、限电导致的历史弃风弃光也会改变分布形状。我的建议是不要硬套参数分布直接用核密度估计或经验CDF。经验CDF最简单u rankdata(w) / (n 1)这里分母用n1而不是n是为了避免样本最大值对应的u变成1否则后面Copula采样时逆变换会取到样本边界上容易失真。用核密度估计更平滑但需要处理边界因为出力是有界的[0, 额定功率]普通高斯核会在边界外漏概率要使用带边界修正的核密度估计。对大多数工程场景经验CDF已经够用计算快还不容易出错。提示如果光伏数据里混入大量夜间零出力或者有显著的多峰分布经验CDF依然可以计算但Copula的唯一性要求边缘分布连续这时要么只取正出力时段要么改成日尺度聚合数据否则后面会出问题。2.2 用Kendall tau而不是Pearson来估计相关参数边缘分布处理完之后数据变成了两个在[0,1]区间上的均匀分布序列。接下来要估计Copula的参数。最常用的入参是秩相关Kendall tau因为它与Copula的参数往往有解析关系。以Gaussian Copula为例它的相关矩阵参数rho和Kendall tau之间的关系是rho sin(pi * tau / 2)这个公式在工程中太常用了。先算历史数据的Kendall tau再转成rho就可以初始化Gaussian Copula的参数。t Copula的散度矩阵也常用同一个转换公式。不同tau值对应的rho如下表秩相关tau对应rho0.000.0000.200.3090.400.5880.600.8090.800.951注意Kendall tau对单调变换不变所以即使风电和光伏的出力量纲不同、分布形态不同也不影响相关结构的估计。这一点比直接用皮尔逊相关稳健得多。如果数据里含有异常值Kendall tau给出的估计也更耐脏。2.3 常用Copula族的性格与选型有了tau和rho还需要回答用哪种Copula。不同Copula相当于不同的“性格设定”适合不同的相依结构。下面是我常用的选型对比Copula族允许负相关尾部相关典型适用场景Gaussian是无基准场景对称依赖尾部不突出t是对称尾部风光极端出力同时出现的频率高于正态假设Clayton基本为正下尾相关低出力同时出现、联合极端低谷Gumbel基本为正上尾相关高出力同时出现Frank是无允许负相关的复杂组合对风电和光伏来说最值得关注的是“同时低出力”的尾部风险。若某基地的天气特征是静稳天气下风和光一起萎靡那么Clayton的下尾相关可能更合适如果高低两头都同时出现极端t Copula通常比Gaussian更贴合。实际项目中不会靠肉眼拍板而是在备选族里分别估计参数计算对数似然和AIC/BIC选出拟合更好的族。补充一个细节t Copula比Gaussian多一个自由度参数nunu越小尾部相关越强nu很大时t Copula趋近Gaussian。所以如果AIC显示t Copula显著优于Gaussian往往说明数据存在尾部相关这个时候不要只用Gaussian。2.4 从Copula采样到还原出力场景拟合完Copula之后生成场景就是一条流水线先从Copula抽样得到均匀分布数对再用边缘分布的逆函数还原成原始出力。Gaussian Copula采样示例写成Python大致如下import numpy as np from scipy import stats # w, p 分别是历史风电、光伏日平均出力数组 # 1. 边缘分布 - 均匀分布 w_u stats.rankdata(w) / (len(w) 1) p_u stats.rankdata(p) / (len(p) 1) # 2. 用Kendall tau估计rho tau, _ stats.kendalltau(w, p) rho np.sin(np.pi * tau / 2) # 3. 从Gaussian Copula中采样 z np.random.multivariate_normal( mean[0, 0], cov[[1.0, rho], [rho, 1.0]], size5000 ) u_s stats.norm.cdf(z) # 采样得到的均匀分布对 # 4. 逆变换回出力量纲 w_scen np.quantile(w, u_s[:, 0]) p_scen np.quantile(p, u_s[:, 1])t Copula采样基本一样只是把多元正态换成多元t分布再把结果按t分布的CDF转成均匀分布nu 6 # 实际中通过极大似然或网格搜索估计 t_z stats.multivariate_t.rvs( locnp.zeros(2), shape[[1.0, rho], [rho, 1.0]], dfnu, size5000 ) u_s stats.t.cdf(t_z, dfnu)逆变换这一环最容易出错如果边缘分布是用经验CDF估计的逆变换可以用np.quantile如果用了带参数的分布就要用该分布的ppf函数。两者结果会有细微差别但量级一致。采样得到的5000个[w_scen, p_scen]对就是初选场景集。3. K-means场景削减的细节聚类中心、概率权重和K值选法3.1 场景进入聚类前先想清楚向量维度Copula生成的初选场景可以是二维点日平均风、光出力对也可以是高维时序曲线比如24小时风电加上24小时光伏拼成48维向量。K-means本身不关心维度它只在你给的距离定义下做划分。但维度不同距离的含义不同结果差异很大。如果场景只是二维的风光出力组合用欧氏距离很自然如果是48维的日曲线直接算欧氏距离时幅度大的时段会主导聚类结果。举例说风电在凌晨出力普遍偏高光伏在午间才高如果两个场站的量级差很多欧氏距离会主要反映风电时段的差异光伏的特征容易被淹没。这时建议先确认所有维度是否在同一量纲下要么统一归一化到各自额定功率要么在聚类特征里把每个时段的权重显式做出来不要糊里糊涂地丢进去。注意K-means聚类时常用欧氏距离它隐含地假设每个维度重要性相同。如果对某些极端时段更关心比如晚高峰时段的风光互补性可以考虑给这些维度加权或者改用加权的欧氏距离。3.2 K-means、多次重启和“代表场景”的真实含义K-means的目标是把N个场景分成K组使得组内平方误差和WSS最小。算法本身不复杂初始化K个质心反复迭代分配样本和更新质心直到目标函数不再明显下降。但它对初始质心很敏感随机初始化可能收敛到局部最优。所以实操中要做三件事。第一初始化必须用K-means它让初始质心尽量分散避免一开始就扎堆。第二要多次重启比如n_init20每次用不同的随机种子最后挑目标函数最小的那次结果。第三固定随机种子保证实验可复现。很多项目里场景削减结果不稳定根因往往只是没有固定种子。K-means得到的聚类中心是簇内所有样本的均值向量它不一定是历史上真实出现过的场景。这个点很关键聚类中心是一个“合成场景”可能介于某些实际场景之间但物理上不一定是某一天真实发生的出力曲线。如果下游优化对场景的可实现性要求高比如要求每个场景都必须是可观测的典型日那可以在每个簇里挑一个最接近聚类中心的真实历史场景作为代表而不用聚类中心本身。3.3 K值怎么定才不至于拍脑袋K值太小代表场景丢失太多分布信息K值太大又违背场景削减压缩场景数量的初衷。K值选择常用三个指标肘部图看WSS下降趋势、轮廓系数、Davies-Bouldin指数。肘部图最简单把K从2试到30画WSS与K的关系拐弯最明显的位置就是合理K。轮廓系数则综合了簇内紧密度和簇间分离度取值越高越好Davies-Bouldin指数越低越好。我通常在项目里把K控制在5到20之间。用日平均风光出力做场景集时K取8到15都比较常见如果是小时级时序曲线K可以稍大一些因为维度高、形态组合更多。一个典型的评估思路是先跑一组K记录指标变化再结合下游优化的求解成本做最终决定。理论上K越大削减后保留的信息越多但随机规划求解时间可能成倍上涨所以K值本质是精度与速度的平衡。3.4 概率权重不是摆设K-means削减完之后每个代表场景都要带一个概率权重否则下游优化不知道哪些场景更常见。权重计算很简单p_k n_k / N即第k个簇的样本数n_k除以总样本数N。这K个权重之和恒为1。需要注意两个细节一是权重必须在聚类完成后再统计不能在聚类前就按等概率假设二是如果后续要对代表场景做二次筛选比如剔除物理上不可行的场景那么剩余场景的权重要重新归一化。另外削减后的场景数量虽然变少但场景集仍然应该保持秩相关结构。验证方法很直接算削减前5000个点里风、光的Kendall tau再算削减后12个代表场景按权重展开后的加权相关指标两者不要差太多。如果削减后相关性明显变小往往意味着K取小了或者聚类距离偏好出了问题。4. 从5000组初选场景到12个代表场景的完整实测4.1 演示数据与参数设置我用一份虚构的演示数据走了一遍完整流程方便展示每个环节。假设某风电场额定容量100MW某光伏电站额定容量80MW历史数据是全年365天的日平均出力清洗掉限电日和异常数据后剩350天。参数如下阶段取值边缘分布经验CDFrankdata / (n1)备选CopulaGaussian与t Copula按AIC选择初始场景数5000场景向量风电日平均出力、光伏日平均出力二维聚类算法KMeansk-meansn_init20代表场景数12权重簇内样本数占比这组数据里历史风光的Kendall tau约为0.35说明确实存在中等强度的正相关。用Gaussian Copula和t Copula分别拟合后t Copula的对数似然更高AIC更小所以最终选t Copula自由度nu估计出来大约是5。这代表数据里尾部相关不可忽视。4.2 核心脚本骨架下面的代码把“估计Copula 采样 K-means削减”串成一条完整流水线方便你直接改成自己的数据格式。import numpy as np from scipy import stats from sklearn.cluster import KMeans # 假设 w, p 是长度相同的历史出力数组已清洗 # 1. 边缘分布 - 均匀分布 w_u stats.rankdata(w) / (len(w) 1) p_u stats.rankdata(p) / (len(p) 1) # 2. 估计依赖参数 tau, _ stats.kendalltau(w, p) rho np.sin(np.pi * tau / 2) # 3. 从t Copula采样5000个初选场景 nu 5 # 项目中用AIC选定 t_z stats.multivariate_t.rvs( locnp.zeros(2), shape[[1.0, rho], [rho, 1.0]], dfnu, size5000 ) u_s stats.t.cdf(t_z, dfnu) w_scen np.quantile(w, u_s[:, 0]) p_scen np.quantile(p, u_s[:, 1]) # 4. 组装场景矩阵并做K-means削减 X np.column_stack([w_scen, p_scen]) model KMeans(n_clusters12, initk-means, n_init20, random_state42) labels model.fit_predict(X) weights np.bincount(labels) / len(labels) centers model.cluster_centers_这个骨架故意保持简单方便逐段验证。如果你要扩展成24小时曲线场景只需要把X换成由风、光时序拼接成的多维向量并把前面的边缘分布和Copula按逐时段处理。骨架的整体逻辑不变。4.3 看结果相关性骨架保留了多少跑完流程后我对比了削减前后的几个关键指标如下表所示演示数据结果仅作参考指标削减前5000点削减后12点加权风电平均出力归一化值0.4120.406光伏平均出力归一化值0.2830.291风光的Kendall tau0.350.31风光的Pearson相关系数0.440.40可以看到削减后风光出力的均值和相关性都有轻微缩小但方向没变、骨架还在。这说明K-means在保留联合分布主要趋势上是有效的。如果你想进一步压缩误差有两种做法一是增大K比如从12增到20相关性指标会更接近初选场景二是在聚类时把这些代表场景的权重纳入距离计算让相关性指标本身参与优化而不是事后检验。4.4 代表场景如何进入随机优化得到12个代表场景之后下游随机规划里常见的是场景树形式每个代表场景对应一组风光出力并带一个概率权重。比如用两阶段随机优化目标函数里就会出现对场景概率加权的期望项sum_k p_k * f(x, y_k)其中x是第一阶段的决策变量y_k是第k个场景下的风光出力p_k是权重。这样求解出来的方案既考虑了风电和光伏的联合相关结构又不至于被成千上万个场景拖垮求解速度。还有一个工程细节如果12个代表场景中某个场景在物理上明显不合理比如风、光同时达到额定功率的极端样本出现概率过高要在削减后做一次合理性校验必要时把该场景剔除并把权重重新归一化。不要把这个校验放在聚类前因为聚类前剔除会影响所有簇的结构。5. 我没写进代码注释里的坑Copula拟合与K-means削减的实战避坑5.1 边缘分布“过拟合”会毁掉Copula边缘分布看似简单却是整套流程里最容易出错的地方。不要一上来就用高斯核密度估计去拟合有界出力数据常见后果是在0和额定功率附近产生边界泄漏拟合出的分布会出现明显的高峰或负值区域。更隐蔽的问题是当你选的边缘分布和真实分布偏差很大时后续Copula拟合会把这个偏差传导到相依结构里本来无关的变量也可能被“拟合”出虚假相关性。我的建议是先用经验CDF做基准方案因为它不引入额外参数、不会过拟合。如果追求光滑和样本外性能再用带边界修正的核密度估计并做交叉验证。没有把握的时候经验CDF永远是安全的选择。5.2 只选Gaussian Copula尾部风险容易丢很多现成代码默认用Gaussian Copula因为它实现简单、参数好估计。但对新能源出力这种容易出现“同时极端”的场景Gaussian Copula的尾部相关系数为0也就是说它认为两种变量不会同时取极端值。这会导致生成的场景里风、光同时低出力的极端事件被系统性低估。对于可靠性要求高的调度问题这种低估是致命的。所以至少要把Gaussian和t Copula都拟合一遍用AIC或交叉验证选型。如果数据里显现出明显的非对称尾部比如低出力尾部更厚还要试试Clayton。选型这一步多花半小时后面能省下大量补救成本。5.3 光伏出力含大量零时Copula会失灵光伏电站的原始数据如果包含夜间时段会有大量零出力样本这会让边缘分布出现一个很大的“质量堆积”在0处。Copula理论严格来说要求边缘分布连续否则Sklar定理里的Copula不唯一拟合出来的相依结构会变得不稳定。我在实际项目中撞到过这个问题把24小时逐时光伏序列直接丢进去边缘CDF在0附近出现平台Copula参数估计结果随样本变化很大生成场景也明显不合理。处理办法有三种按推荐程度排序一是把分析尺度放到日平均或日总发电量避开夜间零值二是只取日出到日落时段把“是否出力”和“出力多少”分开建模比如先用一个0-1变量模拟是否发电再用截断分布模拟出力大小三是使用零膨胀Copula的扩展形式但这需要额外实现一般工程上不太值当。我的经验是大多数规划问题用日尺度就能解决没必要硬啃小时级零膨胀建模。5.4 K-means对极端场景不友好可靠性分析要小心K-means本质上是在平方误差意义下找代表性场景它天然偏好把样本归到“平均形态”附近因此孤立且极端的场景很容易被平均掉。举个例子5000个初选场景里可能只有几十个是“风小、光也小”的极端场景它们聚成一簇后聚类中心里的光伏出力可能从0.05被拉高到0.12看起来不算极端了。但恰恰是这个场景决定了备用容量的多少。如果你做的是可靠性或安全校核类分析不能只依赖K-means。建议在削减后的场景集之外把历史数据中或采样样本中排名靠前的极端联合事件单独挑出来作为确定性约束场景加入优化。或者在K-means之前对那些尾部场景做特殊加权提高它们被保留的概率。记住一句话削减是用来压缩“常态”不是用来消灭“风险”。5.5 高维场景的维度诅咒与距离度量选择当场景从二维扩展到24小时、96小时甚至更高维度时K-means的欧氏距离会面临维度诅咒。高维空间里所有点的距离会趋向接近聚类结构变得模糊。Wind和Solar的曲线又各有不同的峰谷形态直接用原始曲线距离聚类结果往往是一堆“平均形态”丢掉了真正的形状差异。我建议在进入K-means前先做两件事。第一对曲线做特征压缩比如用主成分分析保留解释方差95%以上的主成分或者提取每天的风电利用小时数、光伏利用小时数、峰值时刻、爬坡速率等业务特征第二如果仍然想用原始曲线距离可以考虑动态时间规整DTW距离但DTW计算成本高5000个场景上做聚类会比较慢。工程上我更推荐“先降维特征再K-means”速度快结果也更好解释。做了这么多年场景生成我现在的习惯是先画风、光出力的成对散点图和秩相关热力图再决定用哪种Copula。如果秩相关在高低出力区间不对称优先考虑Clayton这类非对称Copula如果尾部同时极端事件概率明显偏高就直接上t Copula。等场景削减完成后再做一次相关性校验而不是只看聚类误差。这套流程不敢说适合所有项目但至少能让你在风光联合场景上少走很多弯路。 
长河君 · 2026/10/21 16:09:39 · 阅读 1.2k 基于Copula与K-means的风光出力场景生成与削减实战
前几年做新能源随机场景时我踩过一个特别深的坑风电场景单独生成了一套光伏场景单独生成了一套最后按蒙特卡洛的索引硬拼在一起。当时觉得没什么直到拿实际运行数据复盘发现模型预测的“风光同时低出力”日子比真实情况少了一半备用容量安排得非常被动。原因后来很清楚——风光出力不是独立变量它们在同一片天气系统下会同步变化这种相关结构如果不在生成阶段保住后面做再多优化都是白搭。后来我把整套流程改成基于Copula理论构造联合分布再用K-means做场景削减才算把这个问题理顺。这篇文章就把这套方法的完整思路、实操步骤和踩坑记录摊开讲。核心是两件事第一怎么用Copula把风电、光伏的边缘分布“捏”成一个尊重相关性的联合分布从而生成大批量初选场景第二怎么用K-means把这些成千上万的场景削减成几十个带概率权重的代表场景让下游的随机调度、可靠性评估跑得动。适合电力系统方向的研究生、刚入行的算法工程师以及任何需要给随机优化提供风光联合场景的开发者。整套流程并不神秘数学门槛也不算高只要理解“边缘分布”和“相依结构”可以分开建模后面就顺了。我会从问题本质一直讲到代码骨架再把我实际遇到的几个坑放在最后希望能帮你少折腾几个晚上。1. 风光出力不是“11”联合建模到底解决什么问题1.1 分开生成再拼接为什么总是对不上很多初学场景生成的人第一步都会先做一个“独立采样”风电用一组时序模型生成场景光伏用另一组模型生成场景最后按场景编号一对一对齐。这样做的隐含假设是风、光出力条件独立。但实际气象过程里风电和光伏共享同一个大气环流背景冷锋过境时风速加大同时云量增多导致辐照下降夏季午后局地热对流可能同时带来大风和云层遮挡连续阴雨天则可能让风电和光伏同时处于低出力区间。这些现象都说明二者之间存在复杂且非线性的相依关系。如果无视这种关系把独立生成的两批场景任意组合生成的联合分布与真实分布会有明显偏差。最典型的问题就是“同时低出力”和“同时高出力”这类极值事件的频率失真。我的经验是独立拼接出来的场景往往让调度结果偏乐观系统看起来有很多冗余电源真到极端天气来的时候却不够用。1.2 相关系数描述不了“一起糟糕”的日子有人可能会说要反映相关性算一个皮尔逊相关系数不就行了不行。皮尔逊相关系数衡量的是线性相关而风光出力之间经常不是线性关系低出力区间的收敛、高出力区间的发散、局部的厚尾行为都不是一个标量能刻画的。更重要的是皮尔逊相关系数对极端值敏感样本里几个异常点就能把系数拉得很高或很低。这时候需要用秩相关来辅助判断。Kendall tau和Spearman rho不关心变量的绝对数值只关心排序的一致性对非线性单调变换是稳健的。但即使知道了秩相关也仍然不足以完整描述联合分布——两只变量的尾巴到底是一起翘还是分开翘秩相关回答不了。Copula理论解决的正是这个问题先把每个变量的边缘分布单独建出来再用一个Copula函数把它们的相依结构完整地连接起来相当于把“每个变量自己的脾性”和“变量之间的合拍程度”分开了各建各的最后再拼在一起。1.3 两阶段思路先扩充场景再压缩场景有了Copula之后还需要回答一个工程问题下游优化到底需要多少场景随机规划里几十个场景通常就是极限几百个就会让求解时间不可接受但统计学上又需要足够多样本来覆盖可能的出力组合。这是一个天然矛盾。所以标准做法是两步走。第一步用Copula从联合分布里抽出大量初选场景比如5000个这一步目的是保证样本对真实相关结构的覆盖面足够广第二步用K-means把5000个场景聚成K类每类取聚类中心作为代表场景并用类内样本占比作为概率权重。这样既保留了联合分布里的主要特征又把场景数量压缩到优化问题能承受的范围。后面讲的所有内容都是围绕这两步展开的。2. Copula选型和拟合我把Sklar定理拆成了可落地的四步2.1 边缘分布先别套正态尊重数据真实形状Copula建模的第一步不是选Copula而是先搞定边缘分布。Sklar定理说设风电出力W和光伏出力S的联合分布函数为F(w,s)那么一定存在一个Copula函数C使得F(w,s) C(F_w(w), F_s(s))其中F_w和F_s就是两个边缘分布函数。这个等式最大的价值在于建模时可以独立估计边缘分布然后单独拟合相依结构。实操中第一步就是把历史出力的原始数值转换成均匀分布U即 u F_w(w)v F_s(s)。边缘分布怎么选很多教材喜欢用Weibull拟合风速、Beta拟合光伏但到了出力层面数据特征往往比教材复杂光伏日平均出力可能带有零点、风电出力的分布可能厚尾、限电导致的历史弃风弃光也会改变分布形状。我的建议是不要硬套参数分布直接用核密度估计或经验CDF。经验CDF最简单u rankdata(w) / (n 1)这里分母用n1而不是n是为了避免样本最大值对应的u变成1否则后面Copula采样时逆变换会取到样本边界上容易失真。用核密度估计更平滑但需要处理边界因为出力是有界的[0, 额定功率]普通高斯核会在边界外漏概率要使用带边界修正的核密度估计。对大多数工程场景经验CDF已经够用计算快还不容易出错。提示如果光伏数据里混入大量夜间零出力或者有显著的多峰分布经验CDF依然可以计算但Copula的唯一性要求边缘分布连续这时要么只取正出力时段要么改成日尺度聚合数据否则后面会出问题。2.2 用Kendall tau而不是Pearson来估计相关参数边缘分布处理完之后数据变成了两个在[0,1]区间上的均匀分布序列。接下来要估计Copula的参数。最常用的入参是秩相关Kendall tau因为它与Copula的参数往往有解析关系。以Gaussian Copula为例它的相关矩阵参数rho和Kendall tau之间的关系是rho sin(pi * tau / 2)这个公式在工程中太常用了。先算历史数据的Kendall tau再转成rho就可以初始化Gaussian Copula的参数。t Copula的散度矩阵也常用同一个转换公式。不同tau值对应的rho如下表秩相关tau对应rho0.000.0000.200.3090.400.5880.600.8090.800.951注意Kendall tau对单调变换不变所以即使风电和光伏的出力量纲不同、分布形态不同也不影响相关结构的估计。这一点比直接用皮尔逊相关稳健得多。如果数据里含有异常值Kendall tau给出的估计也更耐脏。2.3 常用Copula族的性格与选型有了tau和rho还需要回答用哪种Copula。不同Copula相当于不同的“性格设定”适合不同的相依结构。下面是我常用的选型对比Copula族允许负相关尾部相关典型适用场景Gaussian是无基准场景对称依赖尾部不突出t是对称尾部风光极端出力同时出现的频率高于正态假设Clayton基本为正下尾相关低出力同时出现、联合极端低谷Gumbel基本为正上尾相关高出力同时出现Frank是无允许负相关的复杂组合对风电和光伏来说最值得关注的是“同时低出力”的尾部风险。若某基地的天气特征是静稳天气下风和光一起萎靡那么Clayton的下尾相关可能更合适如果高低两头都同时出现极端t Copula通常比Gaussian更贴合。实际项目中不会靠肉眼拍板而是在备选族里分别估计参数计算对数似然和AIC/BIC选出拟合更好的族。补充一个细节t Copula比Gaussian多一个自由度参数nunu越小尾部相关越强nu很大时t Copula趋近Gaussian。所以如果AIC显示t Copula显著优于Gaussian往往说明数据存在尾部相关这个时候不要只用Gaussian。2.4 从Copula采样到还原出力场景拟合完Copula之后生成场景就是一条流水线先从Copula抽样得到均匀分布数对再用边缘分布的逆函数还原成原始出力。Gaussian Copula采样示例写成Python大致如下import numpy as np from scipy import stats # w, p 分别是历史风电、光伏日平均出力数组 # 1. 边缘分布 - 均匀分布 w_u stats.rankdata(w) / (len(w) 1) p_u stats.rankdata(p) / (len(p) 1) # 2. 用Kendall tau估计rho tau, _ stats.kendalltau(w, p) rho np.sin(np.pi * tau / 2) # 3. 从Gaussian Copula中采样 z np.random.multivariate_normal( mean[0, 0], cov[[1.0, rho], [rho, 1.0]], size5000 ) u_s stats.norm.cdf(z) # 采样得到的均匀分布对 # 4. 逆变换回出力量纲 w_scen np.quantile(w, u_s[:, 0]) p_scen np.quantile(p, u_s[:, 1])t Copula采样基本一样只是把多元正态换成多元t分布再把结果按t分布的CDF转成均匀分布nu 6 # 实际中通过极大似然或网格搜索估计 t_z stats.multivariate_t.rvs( locnp.zeros(2), shape[[1.0, rho], [rho, 1.0]], dfnu, size5000 ) u_s stats.t.cdf(t_z, dfnu)逆变换这一环最容易出错如果边缘分布是用经验CDF估计的逆变换可以用np.quantile如果用了带参数的分布就要用该分布的ppf函数。两者结果会有细微差别但量级一致。采样得到的5000个[w_scen, p_scen]对就是初选场景集。3. K-means场景削减的细节聚类中心、概率权重和K值选法3.1 场景进入聚类前先想清楚向量维度Copula生成的初选场景可以是二维点日平均风、光出力对也可以是高维时序曲线比如24小时风电加上24小时光伏拼成48维向量。K-means本身不关心维度它只在你给的距离定义下做划分。但维度不同距离的含义不同结果差异很大。如果场景只是二维的风光出力组合用欧氏距离很自然如果是48维的日曲线直接算欧氏距离时幅度大的时段会主导聚类结果。举例说风电在凌晨出力普遍偏高光伏在午间才高如果两个场站的量级差很多欧氏距离会主要反映风电时段的差异光伏的特征容易被淹没。这时建议先确认所有维度是否在同一量纲下要么统一归一化到各自额定功率要么在聚类特征里把每个时段的权重显式做出来不要糊里糊涂地丢进去。注意K-means聚类时常用欧氏距离它隐含地假设每个维度重要性相同。如果对某些极端时段更关心比如晚高峰时段的风光互补性可以考虑给这些维度加权或者改用加权的欧氏距离。3.2 K-means、多次重启和“代表场景”的真实含义K-means的目标是把N个场景分成K组使得组内平方误差和WSS最小。算法本身不复杂初始化K个质心反复迭代分配样本和更新质心直到目标函数不再明显下降。但它对初始质心很敏感随机初始化可能收敛到局部最优。所以实操中要做三件事。第一初始化必须用K-means它让初始质心尽量分散避免一开始就扎堆。第二要多次重启比如n_init20每次用不同的随机种子最后挑目标函数最小的那次结果。第三固定随机种子保证实验可复现。很多项目里场景削减结果不稳定根因往往只是没有固定种子。K-means得到的聚类中心是簇内所有样本的均值向量它不一定是历史上真实出现过的场景。这个点很关键聚类中心是一个“合成场景”可能介于某些实际场景之间但物理上不一定是某一天真实发生的出力曲线。如果下游优化对场景的可实现性要求高比如要求每个场景都必须是可观测的典型日那可以在每个簇里挑一个最接近聚类中心的真实历史场景作为代表而不用聚类中心本身。3.3 K值怎么定才不至于拍脑袋K值太小代表场景丢失太多分布信息K值太大又违背场景削减压缩场景数量的初衷。K值选择常用三个指标肘部图看WSS下降趋势、轮廓系数、Davies-Bouldin指数。肘部图最简单把K从2试到30画WSS与K的关系拐弯最明显的位置就是合理K。轮廓系数则综合了簇内紧密度和簇间分离度取值越高越好Davies-Bouldin指数越低越好。我通常在项目里把K控制在5到20之间。用日平均风光出力做场景集时K取8到15都比较常见如果是小时级时序曲线K可以稍大一些因为维度高、形态组合更多。一个典型的评估思路是先跑一组K记录指标变化再结合下游优化的求解成本做最终决定。理论上K越大削减后保留的信息越多但随机规划求解时间可能成倍上涨所以K值本质是精度与速度的平衡。3.4 概率权重不是摆设K-means削减完之后每个代表场景都要带一个概率权重否则下游优化不知道哪些场景更常见。权重计算很简单p_k n_k / N即第k个簇的样本数n_k除以总样本数N。这K个权重之和恒为1。需要注意两个细节一是权重必须在聚类完成后再统计不能在聚类前就按等概率假设二是如果后续要对代表场景做二次筛选比如剔除物理上不可行的场景那么剩余场景的权重要重新归一化。另外削减后的场景数量虽然变少但场景集仍然应该保持秩相关结构。验证方法很直接算削减前5000个点里风、光的Kendall tau再算削减后12个代表场景按权重展开后的加权相关指标两者不要差太多。如果削减后相关性明显变小往往意味着K取小了或者聚类距离偏好出了问题。4. 从5000组初选场景到12个代表场景的完整实测4.1 演示数据与参数设置我用一份虚构的演示数据走了一遍完整流程方便展示每个环节。假设某风电场额定容量100MW某光伏电站额定容量80MW历史数据是全年365天的日平均出力清洗掉限电日和异常数据后剩350天。参数如下阶段取值边缘分布经验CDFrankdata / (n1)备选CopulaGaussian与t Copula按AIC选择初始场景数5000场景向量风电日平均出力、光伏日平均出力二维聚类算法KMeansk-meansn_init20代表场景数12权重簇内样本数占比这组数据里历史风光的Kendall tau约为0.35说明确实存在中等强度的正相关。用Gaussian Copula和t Copula分别拟合后t Copula的对数似然更高AIC更小所以最终选t Copula自由度nu估计出来大约是5。这代表数据里尾部相关不可忽视。4.2 核心脚本骨架下面的代码把“估计Copula 采样 K-means削减”串成一条完整流水线方便你直接改成自己的数据格式。import numpy as np from scipy import stats from sklearn.cluster import KMeans # 假设 w, p 是长度相同的历史出力数组已清洗 # 1. 边缘分布 - 均匀分布 w_u stats.rankdata(w) / (len(w) 1) p_u stats.rankdata(p) / (len(p) 1) # 2. 估计依赖参数 tau, _ stats.kendalltau(w, p) rho np.sin(np.pi * tau / 2) # 3. 从t Copula采样5000个初选场景 nu 5 # 项目中用AIC选定 t_z stats.multivariate_t.rvs( locnp.zeros(2), shape[[1.0, rho], [rho, 1.0]], dfnu, size5000 ) u_s stats.t.cdf(t_z, dfnu) w_scen np.quantile(w, u_s[:, 0]) p_scen np.quantile(p, u_s[:, 1]) # 4. 组装场景矩阵并做K-means削减 X np.column_stack([w_scen, p_scen]) model KMeans(n_clusters12, initk-means, n_init20, random_state42) labels model.fit_predict(X) weights np.bincount(labels) / len(labels) centers model.cluster_centers_这个骨架故意保持简单方便逐段验证。如果你要扩展成24小时曲线场景只需要把X换成由风、光时序拼接成的多维向量并把前面的边缘分布和Copula按逐时段处理。骨架的整体逻辑不变。4.3 看结果相关性骨架保留了多少跑完流程后我对比了削减前后的几个关键指标如下表所示演示数据结果仅作参考指标削减前5000点削减后12点加权风电平均出力归一化值0.4120.406光伏平均出力归一化值0.2830.291风光的Kendall tau0.350.31风光的Pearson相关系数0.440.40可以看到削减后风光出力的均值和相关性都有轻微缩小但方向没变、骨架还在。这说明K-means在保留联合分布主要趋势上是有效的。如果你想进一步压缩误差有两种做法一是增大K比如从12增到20相关性指标会更接近初选场景二是在聚类时把这些代表场景的权重纳入距离计算让相关性指标本身参与优化而不是事后检验。4.4 代表场景如何进入随机优化得到12个代表场景之后下游随机规划里常见的是场景树形式每个代表场景对应一组风光出力并带一个概率权重。比如用两阶段随机优化目标函数里就会出现对场景概率加权的期望项sum_k p_k * f(x, y_k)其中x是第一阶段的决策变量y_k是第k个场景下的风光出力p_k是权重。这样求解出来的方案既考虑了风电和光伏的联合相关结构又不至于被成千上万个场景拖垮求解速度。还有一个工程细节如果12个代表场景中某个场景在物理上明显不合理比如风、光同时达到额定功率的极端样本出现概率过高要在削减后做一次合理性校验必要时把该场景剔除并把权重重新归一化。不要把这个校验放在聚类前因为聚类前剔除会影响所有簇的结构。5. 我没写进代码注释里的坑Copula拟合与K-means削减的实战避坑5.1 边缘分布“过拟合”会毁掉Copula边缘分布看似简单却是整套流程里最容易出错的地方。不要一上来就用高斯核密度估计去拟合有界出力数据常见后果是在0和额定功率附近产生边界泄漏拟合出的分布会出现明显的高峰或负值区域。更隐蔽的问题是当你选的边缘分布和真实分布偏差很大时后续Copula拟合会把这个偏差传导到相依结构里本来无关的变量也可能被“拟合”出虚假相关性。我的建议是先用经验CDF做基准方案因为它不引入额外参数、不会过拟合。如果追求光滑和样本外性能再用带边界修正的核密度估计并做交叉验证。没有把握的时候经验CDF永远是安全的选择。5.2 只选Gaussian Copula尾部风险容易丢很多现成代码默认用Gaussian Copula因为它实现简单、参数好估计。但对新能源出力这种容易出现“同时极端”的场景Gaussian Copula的尾部相关系数为0也就是说它认为两种变量不会同时取极端值。这会导致生成的场景里风、光同时低出力的极端事件被系统性低估。对于可靠性要求高的调度问题这种低估是致命的。所以至少要把Gaussian和t Copula都拟合一遍用AIC或交叉验证选型。如果数据里显现出明显的非对称尾部比如低出力尾部更厚还要试试Clayton。选型这一步多花半小时后面能省下大量补救成本。5.3 光伏出力含大量零时Copula会失灵光伏电站的原始数据如果包含夜间时段会有大量零出力样本这会让边缘分布出现一个很大的“质量堆积”在0处。Copula理论严格来说要求边缘分布连续否则Sklar定理里的Copula不唯一拟合出来的相依结构会变得不稳定。我在实际项目中撞到过这个问题把24小时逐时光伏序列直接丢进去边缘CDF在0附近出现平台Copula参数估计结果随样本变化很大生成场景也明显不合理。处理办法有三种按推荐程度排序一是把分析尺度放到日平均或日总发电量避开夜间零值二是只取日出到日落时段把“是否出力”和“出力多少”分开建模比如先用一个0-1变量模拟是否发电再用截断分布模拟出力大小三是使用零膨胀Copula的扩展形式但这需要额外实现一般工程上不太值当。我的经验是大多数规划问题用日尺度就能解决没必要硬啃小时级零膨胀建模。5.4 K-means对极端场景不友好可靠性分析要小心K-means本质上是在平方误差意义下找代表性场景它天然偏好把样本归到“平均形态”附近因此孤立且极端的场景很容易被平均掉。举个例子5000个初选场景里可能只有几十个是“风小、光也小”的极端场景它们聚成一簇后聚类中心里的光伏出力可能从0.05被拉高到0.12看起来不算极端了。但恰恰是这个场景决定了备用容量的多少。如果你做的是可靠性或安全校核类分析不能只依赖K-means。建议在削减后的场景集之外把历史数据中或采样样本中排名靠前的极端联合事件单独挑出来作为确定性约束场景加入优化。或者在K-means之前对那些尾部场景做特殊加权提高它们被保留的概率。记住一句话削减是用来压缩“常态”不是用来消灭“风险”。5.5 高维场景的维度诅咒与距离度量选择当场景从二维扩展到24小时、96小时甚至更高维度时K-means的欧氏距离会面临维度诅咒。高维空间里所有点的距离会趋向接近聚类结构变得模糊。Wind和Solar的曲线又各有不同的峰谷形态直接用原始曲线距离聚类结果往往是一堆“平均形态”丢掉了真正的形状差异。我建议在进入K-means前先做两件事。第一对曲线做特征压缩比如用主成分分析保留解释方差95%以上的主成分或者提取每天的风电利用小时数、光伏利用小时数、峰值时刻、爬坡速率等业务特征第二如果仍然想用原始曲线距离可以考虑动态时间规整DTW距离但DTW计算成本高5000个场景上做聚类会比较慢。工程上我更推荐“先降维特征再K-means”速度快结果也更好解释。做了这么多年场景生成我现在的习惯是先画风、光出力的成对散点图和秩相关热力图再决定用哪种Copula。如果秩相关在高低出力区间不对称优先考虑Clayton这类非对称Copula如果尾部同时极端事件概率明显偏高就直接上t Copula。等场景削减完成后再做一次相关性校验而不是只看聚类误差。这套流程不敢说适合所有项目但至少能让你在风光联合场景上少走很多弯路。