Copula+K-means：电力系统风光场景生成与削减实战解析

做电力系统随机优化的同学早晚都要碰上一件事风电场和光伏电站的出力数据都有很强的不确定性建模的时候不能把它们当成两个独立的随机变量来抽样。同一个气象过程驱动下风速和光照强度往往存在隐藏的相关性比如高温无风天光伏出力满发而风电几乎躺平又比如强对流天气时风电出力拉满但云层压低光照。忽略这种相关性做出来的机组组合或储能配置方案可能严重偏离实际。这个项目做的就是两件事用Copula理论把风光出力的相关结构精确刻画出来生成足够大规模的初始场景再用K-means聚类把场景削减到随机规划模型能够接受的数量级。整个过程我都跑通并验证过这篇博文就把完整思路、代码组织和踩过的坑整理出来。1. 项目思路拆解风光不确定性的“两步走”1.1 风光出力不是独立的先说一个很多人容易忽略的事实风电和光伏出力在同一区域内往往不是独立变量。有人会凭直觉觉得“风大就阴天”实际上两者之间的相关关系远没有那么简单它会随着季节、地理条件和天气系统变化甚至可能在某个季节呈现明显的负相关换个季节又接近独立。如果强行假设独立Monte Carlo抽样时会抽出一大堆“风又大、光又强”或“风很小、光很弱”的场景这些场景在真实系统里几乎不会同时出现。把这些伪场景喂给优化模型结果就是要么让系统过度冗余要么对极端场景估计不足。更好的做法是把“边际分布”和“相关结构”分开处理这正是Copula理论的核心价值。1.2 先生成再削减两步走的思路本项目整体思路分两个阶段。第一阶段基于历史数据的边缘分布和相关结构用Copula模型批量生成N个风光联合出力场景N可以取1000到5000甚至更多。第二阶段用K-means聚类把这些场景削减成K个典型场景每个典型场景附带一个概率值K通常取5到20。这两个阶段缺一不可。如果只做生成不做削减几千个场景直接进入混合整数规划模型规模会爆炸如果只做削减不做合理生成直接拿历史数据聚类样本量不够时无法覆盖足够多的极端情况泛化能力很差。先生成后削减等于先把不确定性的“形状”画完整再用聚类去压缩信息损失可控。1.3 方法和工具选型Copula不是某一种具体的分布而是一族连接函数。Gaussian Copula、t Copula、Archimedean类的Clayton、Gumbel、Frank都属于常用选择。项目里我一开始图省事直接用Gaussian Copula结果发现对尾部相关性刻画不够后来改成“多候选模型 AIC自动选择”稳定性明显提升。聚类部分K-means不是唯一选择层次聚类、DBSCAN、K-medoids也能做但K-means在场景削减场景下性价比最高。它计算快、实现简单、对“凸簇”类场景效果好而场景削减恰恰希望把相似性高的场景归成一堆凸簇。K-medoids更稳健但对高维大样本太慢层次聚类要预先算距离矩阵卡在内存上。注意场景削减用K-means时建议用k-means初始化并多次随机重启否则聚类结果可能因为初始点选的不好而出现局部最优。2. Copula相关性建模核心环节的原理2.1 Sklar定理到底说了什么Copula的理论基础是Sklar定理简单说就是任意一个二维联合分布函数都可以拆成两部分——两个边缘分布和一个Copula函数。Copula函数负责描述变量之间的“秩相关结构”边缘分布负责描述每个变量自己的形状。用风功率和光伏功率举个例子。设风功率的边缘分布为F_w(x)光伏功率的边缘分布为F_pv(y)那么联合分布可以写成C(F_w(x), F_pv(y))这里的C就是Copula函数它把两个边缘分布的累积概率值映射到联合概率。实际操作时先分别把原始数据变换成均匀分布U(0,1)上的变量再对这些均匀变量之间的结构建模。这样做的最大优势是灵活风功率可以用Weibull分布光伏功率可以用Beta分布两边互不干扰相关结构单独用一个Copula描述。2.2 常用Copula函数怎么选每个Copula函数描述了不同类型的相关性特征实际使用时要根据数据特点选。Gaussian Copula由相关矩阵决定对称分布尾部渐进独立。它实现最简单适合相关性较弱、没有明显尾部依赖的数据但遇到极端天气场景容易低估“风光同时低出力”或“风光同时高出力”的概率。t Copula比Gaussian多一个自由度参数也能刻画对称的尾部相关性在极端天气样本足够时拟合效果通常更好。代价是参数估计复杂一些自由度越小尾部越厚。Clayton Copula擅长刻画下尾相关性也就是变量在低值区域容易同时出现对应到风光场景里就是“风电低、光伏也低”的时段更频繁。Gumbel Copula则相反擅长刻画上尾相关对应“风电高、光伏也高”的极端情况。Frank Copula整体偏弱相关且对称适合相关性不强的场景。我处理风光数据时通常把这几种Copula全部作为候选计算各自的似然值和AIC选最小的那个。经验上t Copula和Clayton Copula在风光联合场景中胜出的概率较高但具体还是要看数据不能拍脑袋指定。2.3 边缘分布的三种做法边缘分布建模直接影响Copula参数估计和后续采样质量这一步不能应付。常用做法有三种。第一种是参数分布拟合。风速常用两参数Weibull分布光伏功率常用Beta分布。好处是解析形式清楚、采样方便坏处是实际功率数据的分布形态经常和理论分布对不上尤其是风机有切入风速、额定功率限制光伏有逆变器限幅导致功率分布存在大量截断和堆积点硬套参数分布会在尾部产生较大偏差。第二种是经验分布。直接把历史数据的秩变换当作边缘分布即用rankdata除以样本量得到经验累积概率。这个方法在Copula估计阶段很常用因为它不引入分布假设鲁棒性最好。但后续采样时如果只靠经验CDF的逆变换样本范围不会超出历史极值会削弱场景生成覆盖极端情况的能力。第三种是混合分布。把零出力时段单独用Bernoulli变量建模非零出力部分再用参数分布或核密度估计。这种方法最贴近风光功率的真实统计特征但实现复杂度高一些。本项目最终采用的是“经验分布为主关键节点做线性外推修正”的组合方式算是工程精度和实现复杂度的折中。2.4 Copula参数估计与最优选择Copula参数的估计常见有两种思路。如果只考虑Gaussian Copula可以直接用秩相关转换法先对原始数据做经验均匀变换再用正态分位数函数映射到标准正态空间最后求相关系数矩阵。这个方法不需要迭代优化速度极快但只适用于Gaussian Copula。如果候选Copula包含t、Clayton、Gumbel、Frank等就需要用极大似然估计。每一步迭代要把均匀变换后的样本带入Copula密度函数求对数似然再用数值优化器搜索参数。项目里用了两阶段估计先估计边缘分布参数再把边缘分布参数固定只优化Copula参数。这样比一步联合估计稳定得多计算量也小。模型选择方面AIC准则是最常用的AIC -2 * logL 2 * k其中logL是对数似然值k是参数个数。把五个候选Copula都拟合一遍取AIC最小的那个进入采样环节。注意比较Copula模型时AIC必须在“统一的数据变换”下计算也就是所有候选模型都用相同的一套边缘分布参数否则比较没有意义。3. 场景生成实操从历史数据到风光联合场景3.1 数据预处理与时间结构处理很多教程会直接拿历史风功率、光功率数据去拟合Copula然后每时段独立采样。这个做法有个隐患它会打散一天内部的时间相关性。风电功率和光伏功率在一天内都有明显的时序惯性上午的光伏出力曲线是平滑上升的夜里风电也往往呈现持续几小时的高出力或低出力状态。如果每个时段独立抽样生成出来的场景会高频振荡看起来就很假。本项目处理时间结构的方法是“日场景作为基本单元”。先把历史数据整理成一条条完整的日曲线每一条日曲线由24小时的风功率序列和24小时的光伏功率序列拼接而成一共48维。为了提取风光的相关性再对每条日曲线计算日特征比如日平均风功率、日光伏发电量占比等然后用Copula对这些日特征建模。这个做法的好处是日曲线内部的时间结构来自真实历史样本天然保真Copula只负责刻画“日与日之间”风与光的宏观相关性。场景生成出来的每条曲线可能来自历史日曲线的组合或缩放但不会出现每小时乱跳的情况。3.2 Copula采样与场景拼接的完整流程整个生成流程分为六个步骤。第一步数据清洗。剔除数据明显异常的天比如传感器故障导致的长时间零值或越限值。第二部对每天的日特征序列做经验均匀变换得到两个均匀分布变量序列。第三步用前一章的方法拟合多个候选Copula选出AIC最小的模型。第四步从选定的Copula中采样N组日特征组合。第五步对每组采样得到的日特征从历史日曲线库中搜索特征最接近的若干条曲线按相似度加权合成一条新的日场景曲线。第六步把所有合成后的日场景按“风电24点光伏24点”的向量格式保存为场景矩阵。第五步的加权合成可以用最简单的K近邻匹配实现也可以用线性缩放调整日特征值到采样目标。项目里用的是K近邻匹配加功率归一化窗口K取10到30匹配距离用标准化欧氏距离。这样生成的场景既保留了历史曲线的骨架又通过Copula采样引入了足够的新组合。3.3 采样后的相关性校验生成场景后不能直接拿去聚类先做一次相关性校验。这里要特别注意校验Copula建模效果时看的不是Pearson相关系数而是Kendall秩相关系数tau或Spearman秩相关系数。因为Copula本身描述的是秩相关结构在均匀变换下只有秩相关保持不变。校验方法很直接先计算历史数据日特征的Kendall tau再计算生成场景中对应日特征的Kendall tau两者对比。如果偏差超过0.05就要回查边缘分布变换是否有问题或者Copula候选集是否选错了。项目里我遇到过一次偏差特别大的情况后来发现是数据里有大量零出力日没有做特殊处理导致经验均匀变换后样本分布不均匀Copula参数估计被带偏了。注意如果历史数据里“无风日无光日”这类样本占比很高最好先按日特征把样本分成若干天气类型再分别建Copula否则整体模型会被零值主导。4. K-means场景削减降维但不失真4.1 为什么要用聚类做场景削减场景削减的核心目标是用尽量少的代表性场景替代原始大数据场景集同时让削减前后优化目标函数的期望值变化最小。从数学上看这个问题可以转化成某种最优传输问题但直接求解传输矩阵开销很大。聚类做法是工程上的主流近似把相似场景归为一类用一类一中心作为代表。K-means之所以适合场景削减是因为它的目标函数恰好是类内平方和也就是希望每个场景离所属类中心的欧氏距离尽量小。如果把场景向量看作高维空间中的点这个目标函数本质上就是在最小化“用类中心代表所有场景”所带来的期望平方误差和场景削减追求的目标高度契合。4.2 K-means聚类流程与关键设置K-means的流程不多说四步初始化、分配、更新、迭代收敛。真正影响最终效果的往往是几个容易被忽略的设置。第一场景矩阵的标准化方式。一个场景向量是24点风电加24点光伏光伏在夜间全是0风电却可能持续有出力如果不做标准化欧氏距离会被风电维度的绝对数值主导。建议按风电维度集合和光伏维度集合分别做z-score标准化或者干脆用各自的装机容量归一化让两类维度在距离计算中的权重不是由量纲决定。第二初始中心的选择。Scikit-learn的KMeans默认就是k-means效果不错。真正会被坑的是n_init太小默认10在某些数据上仍然不够建议设到50。固定random_state便于复现但在做正式实验时也可以多试几个种子看稳定性。第三收敛判据。默认的tol1e-4基本够用。如果场景量级特别大可以调低迭代次数上限但代价是可能没到最优解就停了。这个参数一般不需要动。4.3 K值选择与效果评估K值指望着“肘部法则”一下子看出来往往不现实场景削减数据的高维特性导致WCSS曲线通常没有明显拐点。我更推荐“目标导向法”从K5开始每隔5个做一次削减把削减后的场景集代入后续的随机优化模型看优化目标值随K的变化。当目标值的改变小于某个阈值时那个K就是工程上够用的值。聚类质量本身可以用轮廓系数辅助判断。轮廓系数范围在-1到1之间越接近1说明簇内紧凑、簇间分离好。但轮廓系数在K值较大时计算成本高48维数据跑起来也不快。项目里的经验是轮廓系数主要用来排除特别差的K值比如轮廓系数低于0.2的K直接排除并不追求最大化它。削减后每个代表场景的概率按类内场景数量占比计算公式很简单P_k count(C_k) / N这里的count(C_k)是第k个聚类中的场景数量N是初始场景总数。代表场景取聚类中心即可。注意如果聚类前做了标准化聚类中心是在标准化空间里的要反标准化回原始功率空间再保存否则后续优化模型读到的是缩放后的假数据。5. 代码组织与关键实现5.1 整体代码结构这个项目的代码整体分成四块数据读取与预处理、Copula拟合与采样、场景合成、K-means削减。建议文件名也按这个结构分避免几个功能塞在一个脚本里互相纠缠。主流程的伪代码大概是这样# 1. 读入历史数据整理成日场景矩阵 wind_day, pv_day load_history_data() # 2. 计算日特征并做经验均匀变换 w_feature rankdata(wind_day.sum(axis1)) / (N 1) p_feature rankdata(pv_day.sum(axis1)) / (N 1) # 3. 拟合Copula并采样N个新特征组合 copula_params fit_selected_copula(w_feature, p_feature) sampled_w_feature, sampled_p_feature sample_copula(copula_params, N) # 4. 基于采样特征从历史日曲线库合成新场景 new_w_scenes, new_p_scenes synthesize_daily_curves( sampled_w_feature, sampled_p_feature, wind_day, pv_day ) # 5. K-means削减 scene_matrix np.hstack([new_w_scenes, new_p_scenes]) kmeans KMeans(n_clustersK, initk-means, n_init50) labels kmeans.fit_predict(scene_matrix) center_scenes kmeans.cluster_centers_ prob np.bincount(labels) / len(labels)这里的fit_selected_copula和sample_copula是两个核心函数下面细说。5.2 Copula拟合与采样的核心代码以Gaussian Copula为例拟合非常简单。先对经验均匀变量做正态分位数变换然后求相关系数矩阵import numpy as np from scipy.stats import norm, rankdata def fit_gaussian_copula(u_w, u_pv): # u_w, u_pv 是经过边缘变换得到的均匀分布变量 z_w norm.ppf(np.clip(u_w, 1e-6, 1 - 1e-6)) z_pv norm.ppf(np.clip(u_pv, 1e-6, 1 - 1e-6)) z np.vstack([z_w, z_pv]) corr np.corrcoef(z) return corr def sample_gaussian_copula(corr, n_samples): mean np.zeros(2) z np.random.multivariate_normal(mean, corr, sizen_samples) u norm.cdf(z) return u[:, 0], u[:, 1]这里的clip很重要因为日特征如果有大量并列值rankdata除以N1后可能出现0.99999之类的值直接norm.ppf会算出很大的数影响相关系数估计。采样得到u_w和u_pv两个均匀分布变量之后再通过这些变量去匹配历史日曲线。对于t Copula、Clayton Copula直接复用scipy.optimize.minimize和对应密度函数自己写极大似然估计。这个工作量并不大但需要熟悉对数似然函数的写法。5.3 K-means削减代码的细节削减部分用Scikit-learn实现很简单但有几个细节容易出错。首先反标准化必须在保存结果之前完成。其次如果后续优化模型要求场景概率之和为1检查一下np.bincount的返回值除以N之后是否为1。还有一个容易被忽视的问题K-means聚类中心的坐标是“均值向量”它可能不代表一条物理上真实的出力曲线。比如某类里光伏曲线都是正午高峰形状均值中心可能把峰值周围时间点的功率摊平了一点看起来比任何历史曲线都平滑。如果后续模型对曲线形状敏感可以改用类内距离质心最近的场景作为代表场景而不是直接用聚类中心。这个改动很小但能避免一些物理上不合理的代表曲线。from sklearn.cluster import KMeans def kmeans_scene_reduction(scene_matrix, n_clusters, random_state42): km KMeans(n_clustersn_clusters, initk-means, n_init50, random_staterandom_state) labels km.fit_predict(scene_matrix) centers km.cluster_centers_ probs np.bincount(labels) / len(labels) return centers, probs, labels6. 常见问题与排错记录6.1 高频问题速查表现象可能原因解决思路生成场景的风光相关系数与历史差距过大边缘变换不彻底、Copula类型选错、样本量不足先校验Kendall tau再用AIC重新选CopulaK-means聚类结果每次跑都不一样初始中心随机、样本量大导致局部最优增大n_init固定random_state用k-means代表场景出现大量负功率或零值过剩采样后逆变换处理不当标准化未还原检查反标准化逻辑用布尔掩码处理零值聚类后轮廓系数极低K值不合适或场景矩阵未按维度块标准化排除K值再检查风电和光伏是否量纲失衡Copula参数估计迭代不收敛样本含大量重复值似然函数存在平台区先做零膨胀处理或改用经验分布做边缘变换削减后优化目标与原始场景差异太大K值偏小丢失极端场景适当增加K值或对极端场景单独保留6.2 实操中的经验教训项目里遇到最多的坑集中在零值处理上。风功率数据里大量0光伏数据更是夜间一整段都是0处理不当会让经验均匀变换后等间距出现一段“平坦区”Copula似然函数在这个区间基本没有辨识度。后来我把每天曲线先分成白天和夜间两段光伏只在白天时段建模风电保留全天24点虽然实现多一点但边缘分布的质量提升非常明显。另一个值得说的问题是K-means的标准差敏感性。我发现直接用原始功率值聚类时风电多的季节场景会主导簇划分光伏的差异几乎体现不出来。后来把风电24维和光伏24维在拼接前各自除以该维度的历史标准差相当于对维度做了加权平衡聚类结果的业务解释性好了很多。还有一个容易犯错的地方是K-means聚类后的代表场景如果直接用来跑随机优化可能导致优化结果对某个代表场景过拟合。项目中我加入了边缘场景保护机制聚类结束后额外保留若干Copula采样中的极端场景比如风光出力最低的5%和最高的5%让它们在场景集里单独占一个低概率槽位。这样既控制了场景数量又防止削减过程把极端但重要的场景抹掉。关于场景数量个人经验是初始生成3000到5000个场景削减到10到20个综合效果最好。太少相关系数的统计误差大太多Copula采样耗时增加削减收益也不明显。最后再分享一个小技巧。Copula采样后的场景矩阵可以直接可视化把风功率和光伏功率的散点图按聚类标签着色。如果削减后的簇边界和Copula等高线形状不匹配说明K-means距离度量和Copula相关结构有冲突这时候可以尝试在削减前先用PCA把48维场景降到20维左右聚类效果往往更稳定。这个项目后期的版本就加了这个处理整体鲁棒性提升了一个台阶。