基于Copula与蒙特卡洛的风光出力场景生成及Kmeans缩减方法

我最早被风光出力场景分析折磨是在做一个含风电、光伏、储能的容量配置项目。项目本身不复杂在给定负荷的前提下回答装多少风、装多少光、配多少储能才能让全生命周期成本最低、供电可靠性达标。但风光出力的不确定性让这个问题变得很棘手——今天和明天的出力曲线完全不同不同季节更是天差地别。摸着均值曲线做容量配置结果大概率要么过度投资要么容量不足。当时团队的技术路线最终收敛到一句话用蒙特卡洛生成海量风光出力场景用Copula函数刻画风电与光伏出力的相关性再用Kmeans做场景缩减提取少量有代表性的发电场景。这条路线我在后来的多个项目中反复验证值得把这套方法论完整记录下来。1. 风光出力场景这条路是怎么被逼出来的1.1 只拿均值曲线做规划会得到什么结果很多初学者上来就画一条风光出力的典型日曲线然后拿去算电量平衡。这个做法在项目汇报时看起来很直观但做规划时隐患很大均值会抹平极端情况而极端情况恰恰决定了系统的可靠性短板。举个具体例子。某天傍晚恰好无风风电出力接近0光伏也已经落日归零此时晚高峰负荷还在高位这个时段全靠储能和备用电源扛。均值曲线上这个时段可能看起来还有0.2左右的出力直接套进优化模型后储能容量会被显著低估。反过来盛夏午后光伏满发而风电几乎为零如果只看均值可能会错误地认为风光能同时提供出力导致系统在无风晴天的午间出现容量缺口。容量配置项目中储能的设计尤其敏感。储能的本质就是对冲时间维度上的波动用均值曲线做储能容量优化算出来的结果一定是电量过剩但功率不足实际运行中会频繁触发削负荷。这个教训我体会很深不确定性分析不是为了显示自己手段多而是因为这玩意儿不分析方案根本落不了地。1.2 场景分析在容量配置中的两种用法场景法在规划模型中有两种典型用法。第一种是把场景当作约束条件逐一校验生成N个场景后系统在每一个场景下都不能失负荷以此校核容量配置的可靠性。这种方式逻辑简单、责任清晰缺点是如果场景数量太大逐一约束会让优化模型变成巨型问题。第二种是把场景纳入期望值目标函数容量投资成本加上所有场景运行成本的期望整体最小化。在这种范式下每个场景都带一个概率权重目标函数写成投资成本加 Σ(场景概率 × 场景运行成本)。这个方式能自然地反映大概率场景影响经济性、小概率极端场景影响可靠性的折中关系也是我实际项目中最常用的方法。但无论哪种用法直接扔2000个场景给优化求解器都不现实。这就需要一个前置环节先把2000个场景缩减到10到20个带概率的典型场景。场景缩减的意义就在这里——它不追求保留每一条曲线而是保留整个场景集合的统计特征。所以整条技术链路实际上是生成-缩减-优化三段式。2. 数据准备与边际分布场景生成的地基2.1 数据清洗与时间分辨率的选择风光出力场景分析的第一步是数据这一步做不好后面全是白费。数据来源一般是场站实测或者气象再分析数据我通常要求至少1年、推荐2到3年的历史出力数据时间分辨率按项目需求选。容量配置级别的规划1小时分辨率足够如果做日内调度或储能充放电策略建议搞到15分钟。数据清洗有几个容易忽略的点。第一是限电时段风电场或者光伏电站被调度限出力时SCADA记录的出力并不是真实可发能力直接保留会让低频段的概率密度虚高。第二是检修和通讯故障时段出力会出现长时间跳零这类数据要剔除或插补。第三是异常尖峰比如某条曲线突然出现超过装机容量的数值多半是采样错误。清洗之后要做归一化按各场站装机容量把出力折算到[0,1]区间。这一步很多人觉得理所当然但它有一个隐藏作用归一化之后风电和光伏可以用同一套概率分布框架来分析后面构建Copula时才不会因为量纲差异出问题。2.2 边际分布选型为什么我倾向Beta分布归一化后的风光出力呈现出典型的两头高、中间低形态大量时刻出力接近0比如夜晚的光伏、静风期的风电另一个尖峰出现在满发附近比如大风天的深夜风电、晴天的正午光伏。中间出力反而相对稀疏。风电出力常用的分布有Weibull、Lognormal等光伏出力有人用Beta分布。但在实际拟合中我发现对归一化后的出力数据Beta分布是最省心的选择——它定义在[0,1]区间能自然处理有界数据形状参数灵活能同时拟合零出力堆积和满发堆积这两个尖峰。用Weibull拟合光伏这类非对称、双峰特征明显的数据往往效果欠佳中间段的概率密度会显著偏离经验分布。核密度估计KDE理论上拟合得最好但它在边界处会泄漏到[0,1]之外做分位数逆变换时还容易产生不合理的尾部波动实现成本偏高。所以我的一般建议是首选Beta分布拟合结束后用KS检验和Q-Q图确认效果如不理想再改用混合分布或KDE。2.3 用概率积分变换把历史数据压到均匀区间Copula函数只处理[0,1]区间上的均匀变量所以拿到边际分布之后必须把历史数据做一次概率积分变换对每个出力值x计算F(x)——也就是该值在已拟合边际分布下的累积概率。变换后的U序列理论上服从[0,1]均匀分布。这一步做完一定要做检验别跳过。最简单的方式是画直方图看U序列是不是均匀覆盖更严谨一点做KS检验判断均匀性。我踩过一次坑某次项目里数据没清洗干净限电时段让Beta分布在零出力附近拟合得极差变换后的U序列在0附近堆成了一堵墙后续Copula参数估计全被带偏。后来加了数据清洗和分布诊断U序列的直方图才变得平直。理论上U序列是所有Copula建模的原料原料不行后面的高斯空间变换和蒙特卡洛采样都无从谈起。3. Copula函数建模风与光的依赖关系怎么刻画3.1 为什么不能直接拍一个相关系数了事风电场和光伏电站如果地理位置接近出力之间往往存在一定的相关性同一个天气系统既影响风速又影响云量无风的静夜可能同时导致风电出力低和光伏出力为零白天的强对流可能让光伏骤降的同时带来大风。如果假设风电和光伏独立生成场景时可能会把高风高光同时出现的频率放大让规划方案过度乐观如果只用皮尔逊相关系数做一个线性近似又无法刻画尾部联动。Copula函数解决的是边际分布与依赖结构的分离问题。边际分布描述单个变量的概率形态Copula描述多个变量之间如何纠缠。用Sklar定理的语言说联合分布可以被拆成边际分布加Copula函数。这样做的好处非常明显我可以给风电配Beta分布、给光伏配Beta分布然后单独用一个结构来刻画风弱光弱同时发生这类极端事件的相关性而不必强行构造一个多维参数分布。Copula还有个直观解释它像一把锁把两个变量边缘上各自独立的概率值锁在一起。两个U序列都是均匀分布如果没有相关性散点图就是一片均匀的云如果存在相关性散点图就会呈现出朝对角线聚拢、或朝左下角/右上角聚拢的结构。3.2 Copula族的选择与参数估计步骤实际项目中常见的Copula族就这么几个选型时主要看对称性、尾部依赖特征和参数估计的稳定性。Copula族尾部依赖适合描述的风光场景Gaussian无尾部依赖风电和光伏相关性较弱、没有明显极端联动t-Copula对称尾部依赖极端天气下风电和光伏可能同时大幅波动Clayton下尾相关无风无光同时发生即系统最危险的时段Gumbel上尾相关大风暴可能让风电满发、但光伏因云层骤降Frank无尾部依赖温和的正负相关不考虑极端事件选型的依据一部分来自经验一部分来自数据。先画出风电和光伏历史出力变换到均匀空间后的散点图看尾部是偏向左下角、右上角还是两个角都偏。如果左下角堆积明显优先试Clayton如果两个角都堆积优先试t-Copula如果整体温和无明显极端联动Gaussian就行。不要一上来就套t-Copula它的自由度参数在小样本下估计波动很大。参数估计用两步极大似然法IFM。第一步用历史出力估计Beta分布参数得到边际分布第二步用概率积分变换得到U序列然后建立想要测试的Copula对U序列做极大似然估计。以Gaussian Copula为例实际操作就是先把U序列做标准正态逆变换成Z序列然后直接算Z序列的相关系数矩阵t-Copula则在相关系数基础上多一个自由度参数可以用数值优化估计。3.3 拟合结果怎么看Copula拟合得好不好不能只看参数估计是否收敛要看生成数据能不能复现历史数据的依赖结构。我常用的验证手段有三个。第一画对比散点图把历史U序列的散点图和从拟合Copula中抽样的散点图并排看眼睛就能判断依赖结构像不像。第二计算经验Copula与理论Copula之间的偏差比如Cramer-von Mises距离数值越小越好。第三比较几个Copula族的AIC或BIC值这个指标直接体现拟合质量和参数数量的平衡。这里有一个容易犯的错只看相关系数。Gaussian Copula和t-Copula可能估计出非常接近的相关系数但它们对极端天气场景的刻画差别很大。在容量配置项目里系统失负荷率往往由尾部场景决定所以尾部依赖拟合比整体相关系数更重要。如果只是为了报告里写一个相关系数那真不如不做Copula直接历史数据重采样更省事。4. 蒙特卡洛生成风光出力场景从概率模型到一批可用的曲线4.1 采样流程与逆变换的完整实现Copula拟合好之后蒙特卡洛生成场景的流程其实很机械先从Copula中抽取相关随机数再做边际逆变换。以Gaussian Copula为例完整流程可以这样写。import numpy as np from scipy import stats # 假设 wind_hist, solar_hist 是清洗并归一化后的历史出力序列 # 1. 拟合Beta边际分布 wind_beta stats.beta.fit(wind_hist, floc0, fscale1) solar_beta stats.beta.fit(solar_hist, floc0, fscale1) # 2. 概率积分变换到均匀空间 u_wind stats.beta.cdf(wind_hist, *wind_beta) u_solar stats.beta.cdf(solar_hist, *solar_beta) # 3. 高斯空间变换并估计相关系数 z_wind stats.norm.ppf(u_wind) z_solar stats.norm.ppf(u_solar) rho np.corrcoef(z_wind, z_solar)[0, 1] # 4. 蒙特卡洛采样生成2000个场景 np.random.seed(42) # 固定种子保证可复现 n_scenarios 2000 mean [0, 0] cov [[1, rho], [rho, 1]] z_sample np.random.multivariate_normal(mean, cov, n_scenarios) # 5. 逆变换回[0,1]均匀空间再经过各自边际分布的逆函数得到出力场景 u_sample stats.norm.cdf(z_sample) wind_scenarios stats.beta.ppf(u_sample[:, 0], *wind_beta) solar_scenarios stats.beta.ppf(u_sample[:, 1], *solar_beta)如果你的场景是24小时曲线而不是单个时间断面的出力操作方式略有差异把全天24个时刻分别建立出动力的联合分布模型或者更常见的做法是在多个典型时段上分别构建Copula再按时间维度拼接。对于容量配置项目我通常按季节或者按典型日类型先分组每一组单独拟合数据和采样这样能保留日出力和季节差异。4.2 场景数量设置2000个够不够很多第一次做场景分析的人会问到底生成多少个场景才够我的答案是1000个起步2000个比较稳如果条件允许可以到5000个但5000个以上对大多数规划问题的边际收益就开始下降了。判断够不够不能拍脑袋要做收敛性检验。分别生成500、1000、2000、5000个场景对比各组的出力期望值、方差、5%分位数和95%分位数。当指标随场景数增加而趋于平稳时说明样本量已经足够。比如某次项目中2000个场景的95%分位数比1000个场景提升了0.3个百分点但3000、5000个场景几乎没变那就说明2000个已经到位。场景数量还取决于后续Kmeans缩减的最小数。如果你的目标是缩减到10个典型场景那原始场景太少会导致聚类结构不稳所以宁可多生成一些。但过多也没必要——生成20000个场景后面聚类时计算距离矩阵的内存和时间开销都会急剧上升收益却微乎其微。4.3 采样优化的一个小改进拉丁超立方简单随机采样有一个缺点稀疏区域容易出现扎堆极端场景的覆盖不稳定。今天生成2000个场景和明天重新生成2000个可能极端分位数会有可感知的差异这在工程上很恼火。一个低成本改进是拉丁超立方采样LHS。思想是把每个变量的[0,1]区间均匀切分成N个小区间在每个区间内各抽一个点然后随机重排来保留变量间的相关性。这样保证每次采样都能覆盖全部分布范围尤其是尾部的极端区域而且蒙特卡洛方差显著降低。如果已经走了Gaussian Copula的路径LHS可以直接作用在高斯空间的Z序列上先在两个维度上分别做分层采样再通过相关性矩阵施加结构最后逆变换回出力空间。代码量不大但效果非常明显。我现在的项目里标准随机采样基本只做基线对比正式计算一律LHS。5. K-means场景缩减2000个场景压缩到10个5.1 特征设计与距离度量的坑蒙特卡洛生成完2000个场景后下一步就是Kmeans缩减。但很多人直接在这一步踩坑问题大多出在特征设计上。先说最简单的形式如果生成的是单时段的风电和光伏出力每个场景是一个二维向量[风电出力, 光伏出力]直接在二维空间做Kmeans就行距离用欧氏距离概念清晰。但如果场景是24小时出力曲线事情就复杂一点。我的做法是将风电24小时曲线和光伏24小时曲线首尾拼接构造一个48维特征向量再在这个向量空间里聚类。这里有两个坑。第一个坑是量纲不统一。如果除了出力还塞进了温度、辐照度等物理量量纲大的特征会主导欧氏距离聚类结果完全失真。所以要么只使用归一化后的出力要么对全部特征做标准化再聚类。第二个坑是Kmeans聚类中心是簇内平均得到的典型场景可能不是一条实际出现过的出力曲线。对容量配置类的优化模型这个特性通常不影响使用但如果你需要拿一条真实曲线去跟电网调度核对可以考虑K-medoids它会从原始场景中挑一条真实存在的曲线作为代表。5.2 K值怎么选肘部法则不是唯一依据K值选择是场景缩减里最艺术的部分。理论派会说用肘部法则加轮廓系数实践派会告诉你最终看求解器能不能扛得住。肘部法则的做法很简单对不同的K值分别计算簇内平方和SSE画一条随K递增而递减的曲线找拐点。轮廓系数衡量每个样本与簇内其他样本的紧密度和与相邻簇的分离度数值越大越好。但我的实际经验是K值往往不是纯统计数据说了算。规划模型里动辄加几十个整数变量每多一个场景运算量就上一个大台阶。我记得某次项目里轮廓系数在K8附近最好但优化模型在K6时还能在两个小时内解完K8时就跑到了十个小时。最终我们取K6并且把误差控制在了可接受范围内。场景缩减的目的本来就是服务决策不是追求聚类指标的完美。另外K值还要兼顾极端场景的保留。通过Kmeans得到的10个场景大概率会落在高概率区域极端低风低光的危险场景可能被平均掉。如果你关心的是系统可靠性建议在缩减后的场景集之外额外手动保留2到3个历史极端场景作为约束校验场景。这个操作成本低但对结果可靠性提升非常大。5.3 缩减后场景的概率权重与质量验证Kmeans聚类完成后每个簇的代表场景就是质心但这个质心本身不携带概率信息。概率必须从原始场景集合中继承每个簇的概率等于该簇包含的场景数除以总场景数。from sklearn.cluster import KMeans features np.hstack([wind_scenarios, solar_scenarios]) k 10 kmeans KMeans(n_clustersk, random_state42, n_init10) labels kmeans.fit_predict(features) # 每个簇的代表场景质心 representatives kmeans.cluster_centers_ # 每个簇的概率权重 probabilities np.bincount(labels) / len(labels)缩减完别急着丢进优化模型先做质量验证。我一般比较三样东西缩减前后场景集合的均值出力、标准差、以及风力与光伏出力的相关系数。如果缩减前后这三个统计量的偏差超过5%说明K值太小或者特征设计有问题需要回头调整。散点图对比也很直观把2000个原始场景画成灰色散点再把10个缩减后的代表场景用红色标记叠加在同一张图上。要是红色点分布的趋势与灰色云团明显不一致那就说明缩减过程把依赖结构破坏了这时候习惯性增加K值或者检查特征是否需要配合标准化。6. 实测踩坑记录与经验总结6.1 Copula拟合中最容易翻车的三个细节第一条历史数据不干净边际分布全盘皆输。我前面提过限电时段堆积的案例这是Copula建模范式中最经典的翻车点。一旦U序列没有呈现出均匀分布后续所有相关性估计都会失真。建议每次完成边际分布拟合后先看一眼U序列的直方图再往前走。第二条样本量不足时硬上复杂Copula。t-Copula自由度参数在小样本下极不稳定有时一次估计出10另一次估计出60结果毫无意义。如果历史数据只有半年我的建议是老老实实用Gaussian或者Frank哪怕明知道有尾部相关也不要为了追求精确而引入高方差估计。第三条用错误的对象做拟合。Copula拟合的输入必须是均匀分布边际变换后的U序列而不是原始出力。有人图省事直接拿归一化后的出力数据去拼Copula得到的结果没有任何数学依据。这条看着基础但我在复查别人代码时真的见过好几次。6.2 聚类缩减与优化模型衔接要注意什么场景缩减的结果要正确地进入优化模型需要明确每个场景的概率权重并且保证所有场景概率之和为1。这个看似简单但如果在代码里做了标准化却忘记把权重带进目标函数优化结果会被概率权重严重的场景主导计算出的期望成本会偏差很大。另一个衔接问题是时间连续性。如果在规划模型中需要模拟连续多日的运行而Kmeans缩减出来的是独立日场景那么场景之间就没有时序关联储能电量的跨日转移会算不准。遇到这种情况我倾向于采用典型日日间转场约束的改良方案或者对多日场景进行聚类而不是简单对单日场景做缩减。6.3 我现在的推荐操作流程这么多年做下来我在风光出力场景分析上基本固定了这样一套操作流程获取并清洗至少一年的场站级历史出力数据剔除限电、检修、异常数据归一化到[0,1]。对风电和光伏分别拟合Beta边际分布通过KS检验和Q-Q图确认拟合质量。做概率积分变换到均匀区间用直方图确认U序列均匀性。观察历史U序列散点图的尾部形态初选Copula族用极大似然估计参数对比AIC/BIC。用LHS方法从Copula中抽取1000到2000个随机场景做收敛性检验确认场景数量合理。对场景做Kmeans缩减用肘部法则初步定K结合模型计算能力敲定最终K值。验证缩减前后场景集合的均值、方差与相关性确认质量合格后再进入优化模型。这套流程每一步都有对应的检查点不会让错误扩散到下一环节。我在实际项目中体会到这套链条里最容易出问题的往往不是算法本身而是数据质量。边际分布歪了Copula拟合再精细也是白搭Kmeans再快如果特征设计把不同特性的场景揉在一起缩减出来的典型场景也代表不了真实情况。做这类项目宁可每个环节多停下来看图也不要追求一口气跑完。先用散点图看清尾部形态再决定选哪种Copula先用K值灵敏度测试摸清边界再做最终缩减。慢即是快。