1. 为什么纯蒙特卡洛会在“时序相关”面前失灵先交代一个背景MCMonte Carlo蒙特卡洛方法做场景生成在电力系统、能源调度、金融风险这些领域里已经算常规操作了。思路也不复杂——对随机变量的概率分布做大量采样用成千上万条样本序列去逼近真实世界里可能出现的各种情况。光伏出力、风电功率、负荷曲线这些不确定性比较强的时序数据都能用MC方法生成一堆“如果明天是这样的天气/负荷”的场景集合。但这里有个很微妙的问题。如果你只是对每个时间点的分布单独采样然后把它们硬拼成一条时序曲线生成出来的场景在“形状”上会很怪一会儿冲到峰值、一会儿跌到谷底相邻时间点之间没有任何连贯性。真实世界不是这样的——风速不会从0瞬间跳到15m/s再瞬间掉回2m/s光伏出力也不会上午十点还是满发、十点零一分就变成零。相邻时刻之间是有惯性、有记忆效应的这个在专业术语里叫自相关两个风电场之间、或者风电和负荷之间也可能存在空间上的联动关系这叫互相关。两者合在一起就是标题里说的“时序相关性”。纯MC采样的致命缺陷就在这里。它只刻画了每个时间点各自的边缘分布完全没有刻画时间维度和空间维度上的相关结构。用这种场景去做随机优化算出来的调度策略会过于保守或者过于激进因为模型看到的“极端场景”很多在物理上根本不可能出现——同一时刻所有风电场都满发、下一时刻全部归零这种事情现实中不会发生。所以场景生成这个方向核心难点从来不是“怎么采样”而是“怎么把相关性结构正确嵌入到采样过程里”。后面我会先讲清楚生成端怎么做再讲削减端怎么做。削减也是个同样关键的话题MC采样动不动就是几千上万条场景直接塞进随机规划模型里根本算不动必须用算法把它们压缩到十几个、几十个有代表性的场景同时尽量不损失原始概率分布的信息。这两件事串起来才是完整的技术链路缺一环都没法落地。2. 场景生成实操从独立采样到相关性植入的完整链路2.1 第一步把样本频率定在“够用”的粒度上时序场景生成的第一步不是写采样代码而是先想清楚时间分辨率和场景长度。这两个参数直接决定后续所有计算量。以电力系统为例常见的做法是1小时一个点、取24个点代表一天或者15分钟一个点、取96个点。分辨率越细能刻画的波动细节越多但需要估计和存储的相关参数也会变多——96个时间点的相关矩阵就是96×96接近一万个参数样本量不够的话很容易估计出噪声来。我的建议是先想明白下游任务到底需要多细的粒度。如果是做日前调度1小时分辨率足够了如果是做日内滚动优化可能要15分钟甚至5分钟。不要一上来就追求最细的分辨率因为采样时的方差、削减时的计算量都是随着维数增长的。先用粗粒度把整个流程跑通确认相关性处理逻辑没问题再逐步加细。样本规模也有讲究。很多论文里生成2000条、5000条甚至10000条场景看起来越多越好但要注意MC的收敛速度是O(1/√N)想让精度提高一倍样本量得翻四倍。实际经验是如果最后要削减到10个场景初始生成500到1000条就够了再多对最终结果的影响非常有限纯粹浪费计算资源。这里有个更合理的做法是先估算下游模型对场景数目的敏感性用50、100、200条分别试跑看目标函数值变化幅度。变化不大就用小的变化明显就加大。后面我会专门讲这个。2.2 第二步用Cholesky分解把相关结构塞进白噪声处理时序相关性的经典手段是协方差矩阵Cholesky分解。它的原理是这样的假设你有m个随机变量比如m个风电场的出力它们之间的相关性用一个m×m的相关系数矩阵R来描述。你先从标准正态分布里独立采样得到一组互不相关的向量Z然后通过矩阵分解得到下三角矩阵L使得R L·Lᵀ用L去线性变换Z得到的新向量X L·Z就自动带上了R描述的相关系数结构。放到时序场景生成里的具体操作是把每个时间点当作一个变量先根据历史数据估计出整个时间跨度上的相关矩阵包括同一风电场不同时刻的自相关和不同风电场同一时刻的互相关然后对它做Cholesky分解再用分解得到的下三角矩阵去变换独立正态采样序列。核心代码如下import numpy as np # 假设有3个风电场24个时间点生成24*372维的相关场景 n_times 24 n_sites 3 n_total n_times * n_sites # 1. 根据历史数据估计相关矩阵简化示例用随机矩阵代替 # 实际使用中应该用历史出力序列直接计算 np.corrcoef R np.eye(n_total) # 这里应该替换为真实估计的相关矩阵 # 2. Cholesky分解 L np.linalg.cholesky(R) # 3. 采样独立标准正态随机向量 Z np.random.normal(size(n_total, n_samples)) # n_samples是场景数 # 4. 线性变换得到带相关性的样本 X L Z # 形状: (n_total, n_samples) # 5. 将结果重塑为 (n_sites, n_times, n_samples) scenarios X.reshape(n_sites, n_times, n_samples)这里面有一个工程上的关键点Cholesky分解要求相关矩阵必须是正定矩阵。但实际用历史数据估计出来的相关矩阵经常不是正定的——数据有缺失、时间序列长度不够、或者某些变量之间线性相关度极高都会导致矩阵半正定甚至不定。遇到这种情况直接做Cholesky分解会报错。我的处理办法有三个一是检查数据质量把缺失过多的时段剔除二是对估计出来的相关矩阵做特征值修正把所有小于某个阈值的特征值拉到一个正数三是改用特征分解替代Cholesky分解原理类似但容错性更好。具体怎么判断这个矩阵的病态程度我在后面问题排查那节会展开讲。2.3 第三步用Copula方法处理更复杂的尾部依赖Cholesky分解能处理线性相关性但真实世界的时序数据经常有更复杂的依赖结构。最典型的是风电出力——小风天和大风天都容易出现“一片区域的风电场同时不出力”或者“同时满发”的极端情况这种在统计上叫尾部依赖。用线性相关系数正态假设去刻画会低估极端场景发生的概率而电力系统最怕的恰恰就是极端的低出力高负荷场景。进阶做法是用Copula函数。它的核心思想是把“边缘分布”和“相关结构”分开建模先用历史数据分别拟合每个风电场出力的边缘分布可以用非参数核密度估计再选一个Copula函数来描述变量之间的依赖结构。常用的有Gaussian Copula和t-Copula——前者计算方便后者尾部依赖更强适合描述风电场出力这种有“同时极值”倾向的数据。具体操作流程是把历史数据通过各自的边缘分布变换成均匀分布序列再用这些均匀分布序列去拟合Copula参数最后从Copula里采样得到带上相关结构且服从真实边缘分布的样本。我在实际项目里的体会是如果数据量够用至少一两年的历史出力数据用t-Copula的效果会比纯Cholesky好不少尤其是在刻画极端场景分布形态方面。但代价是计算复杂度上去了参数估计也更敏感。所以我的建议是先用简单方法把流程跑通再逐步升级模型复杂度不要一上来就上Copula否则出了问题根本分不清是相关性建模错了还是边缘分布拟合错了。3. 场景削减到底在“削”什么从概率距离到代表性场景3.1 同步回代消除法的完整操作流程生成几千条场景以后不可能直接把原始样本丢进优化模型里因为计算量完全不可控。前面提过一句——8G内存的机器跑十万条场景的随机规划基本上会卡到怀疑人生。削减的本质是造一个规模更小的场景集使得新旧场景集之间的“概率分布距离”最小。这里的距离一般指Kantorovich距离通俗理解就是两个概率分布之间需要“搬运多少概率质量”才能变得一致。同步回代消除法Scenario Backward Reduction是目前用得最多的算法之一。它的核心逻辑很朴素反复找到“最没有存在感”的场景把它的概率叠加到离它最近的另一个场景上然后淘汰它。每一步都让场景总数减一直到数量满足要求。具体步骤拆开来看是四个环节为每个场景i找到它与所有其他场景j的最小距离min_dist(i)以及对应最近的场景编号在所有场景里找到min_dist最小的那个场景k它就是当前“最容易被替代”的场景把场景k的概率加到它的最近邻场景上然后从场景集合里删掉k更新删除后受影响的邻居关系重复1-3直到场景数达标这个算法的好处是每一步都保证场景集在Kantorovich距离意义上损失最小所以理论性质比较干净。实现上有一点需要注意距离的定义。更细致的做法是逐时间点计算欧氏距离再累加形成两条时序曲线之间的总距离。代码逻辑大概长这样def backward_reduction(scenarios, probs, target_num): scenarios: (n_scenarios, n_times) probs: 每个场景的概率和为1 target_num: 削减后的目标场景数 indices list(range(len(scenarios))) probs np.array(probs) while len(indices) target_num: # 计算所有场景两两之间的距离 dist_matrix np.zeros((len(indices), len(indices))) for i in range(len(indices)): for j in range(len(indices)): dist_matrix[i, j] np.linalg.norm( scenarios[indices[i]] - scenarios[indices[j]] ) np.fill_diagonal(dist_matrix, np.inf) # 找每个场景最近邻的距离 min_dist dist_matrix.min(axis1) nearest_idx dist_matrix.argmin(axis1) # 找到全局最小距离的场景 k min_dist.argmin() j nearest_idx[k] # 把概率合并给最近邻删除场景k probs[indices[j]] probs[indices[k]] indices.pop(k) new_scenarios scenarios[indices] new_probs probs[indices] new_probs new_probs / new_probs.sum() # 重新归一化 return new_scenarios, new_probs大厂和小研究团队在工程实现上的一个差别是距离矩阵的计算方式——用双重循环在场景数超过两三千时会有明显的性能瓶颈。我建议用scipy.spatial.distance.cdist或者直接矩阵广播运算能快一个数量级。3.2 聚类式削减与回代式削减的取舍同步回代消除法虽然经典但也不是唯一的选择。K-means和K-medoids聚类也经常被用于场景削减。聚类思路是把几千条场景按照相似度分成k组用每一组的质心均值曲线或中心点离所有样本最近的那条实际场景来代表这个组组的样本占比就是概率。两类方法放在一起对比各自的优劣比较明显。同步回代消除法得到的结果一定是原始场景里的某一条不会出现“平均出来但物理上不真实”的曲线这一点在新能源出力场景里很关键——把三条完全不同形态的出力曲线平均成一条可能就会得到一条现实中不可能出现的“中间态”曲线。我在实际做风电场景削减时特别在意这个因为调度模型对出力下限和上限很敏感平均值可能把顶峰削平、把低谷抬高这种场景看着合理其实会误导优化结果。K-means的优点是快非常适合场景数非常大比如上万条时做第一轮粗削减缺点是它对初始聚类中心敏感如果初始点选得不好容易收敛到局部最优导致削减后的场景多样性不足。我的常规做法是组合拳先用K-means把一万条粗削减到两百条再用同步回代消除法精削到十条以内。这样兼顾了速度和场景质量而且K-means那步产生的“质心平均”问题在第二轮同步回代时会被消除掉影响可控。削减规模怎么定也有经验可循。我见过很多文章不管三七二十一就削到10个但目标场景数和数据本身的结构复杂度是直接相关的。如果历史数据呈现出明显的季节性差异10个场景可能根本不够表达冬夏两种完全不同的出力形态。我的判断标准是削减之后按照新场景集重新算一遍关键统计量日均出力、最大出力、出力波动方差与原始场景集偏差不超过5%就可以接受。这个校验步骤比任何理论推导都实在。4. 实操中踩过的坑与排查清单4.1 相关矩阵不是正定矩阵Cholesky直接报错做Cholesky分解时最常遇到的报错就是相关矩阵非正定。我刚接触这个方向时踩过一次很深的坑用三四年历史出力数据估计相关系数中间有一段时间数据质量特别差很多时刻的出力值直接填了0结果算出来的相关矩阵行列式接近0一分解就崩。排查思路分两步走。先查数据——把历史出力序列画出来看看有没有异常的大段零值或者缺失值零值过多的时段会人为压低相关系数造成矩阵病态。我后来写了一套简单的数据清洗逻辑任何一天里出力为0的时刻超过全部时刻的30%这条样本直接丢掉宁可数据少一点也要保证相关结构是从真实物理过程里估计出来的。如果确认数据没问题矩阵仍然不正定就做特征值修正。思路是把相关矩阵做特征值分解把所有小于阈值的特征值推到一个正值再重新合成矩阵这一步在文献里叫“矩阵收缩”或“特征值裁剪”。还有一种更实用的方法是只截取跟风电场出力【原始描述】中“时序相关性”密切对应的变量子集去做Cholesky避免把所有变量无差别一锅炖相关性弱的变量会自动让矩阵变得很臃肿。我后来发现把96维的大矩阵拆成4个24维的小矩阵对各自时段内的相关结构分别建模问题会少很多。4.2 削减后概率分布“跑偏”了很多人做完削减只检查场景数对不对、曲线形不相似却忘了检查概率分布的统计特性。有一次我用K-means削完场景目测几条曲线都跟原始场景很接近但拿削减后的场景去跑机组组合模型结果目标函数值和全场景模型差了将近20%。排查下来发现是概率分配出了问题——K-means是均匀采样思想它默认每个簇的样本权重就是样本数占比但这种处理对极端场景不友好某些簇里只有一两条极端场景被K-means硬编码进了质心导致极端情况在削减后的场景集里被放大或丢失。后来我在这类场景削减任务中固定了一个审查习惯削减前后必须逐项比对这些统计量——各类场景的均值出力曲线、出力标准差、概率加权后的极端出力值比如5%分位点和95%分位点。发现问题就调削减方法比如改用K-medoids代替K-means或者用同步回代法。另外还有一个非常隐蔽的坑削减完以后概率一定要重新归一化。因为概率合并过程中舍入误差会积累不归一化的话后续算期望值时所有结果都会系统性偏小。4.3 场景数量到底取多少计算资源与精度的拉扯前面讲了场景数量是计算资源和精度之间的折中。我在做实际项目时给过一个比较通用的方法画出“场景数量-目标函数值”的收敛曲线。选一组目标场景数比如5、10、20、50、100分别跑同一个优化模型记录目标函数值的变化幅度。理论上场景越多结果越接近“真实随机优化”的极限值但这个边际收益是递减的。经验阈值是当场景数量翻倍但目标函数变化不到1%时就说明当前数量已经足够了。这个方法在小规模算例上很快就能跑出来大算例可以先削减到10个场景跑一次再把20个场景跑一次对比两次结果的差异。差异大就翻倍继续差异小就用少的。另外还要看你对“风险偏好”的态度求解的是期望成本最小化场景少点问题不大求解的是CVaR这类关注极端风险的指标场景少了根本覆盖不到尾部此时保守一点反而更稳妥。我个人在这个问题上吃过不少亏现在对极端风险类问题一般会比常规情况多留一倍的场景数量。5. 一个小型可复现案例两风电场出力场景从生成到削减理论和坑讲完了我拿一个完整的迷你案例把整条链路串起来。假设有两个风电场历史出力数据是按小时记录的一共365天每天24个点。目标是生成1000个日场景再削减到10个供一个简单的两机组经济调度模型使用。第一步是估计相关矩阵。把两个场各自24小时的数据拼接成48维向量用历史数据逐日计算相关系数得到48×48的相关矩阵R。注意这里的历史样本是365天等于365组48维样本样本量大约为维数的7倍勉强够用但不算富余。我用的是现场记录中的真实出力数据处理时要格外仔细地剔除检修日和极端停机日否则这些非典型样本会把相关结构拉偏。第二步是采样生成1000条带相关性的正态样本再通过历史数据拟合的逆累积分布函数转换成出力值。这里有个容易被忽视的操作风电出力是0到1之间的比例值直接从标准正态反变换会出现超出边界的情况我一般用截断正态分布做边缘分布或者用Beta分布拟合后再做概率积分变换。第三步是削减。按照前面说的方法先跑K-means粗削到100条再跑同步回代消除法精削到10条。原始1000条场景的总计算时间包括相关性采样、削减和统计校验在普通办公笔记本上大约是两三分钟完全可接受。最后做校验。我对比了削减前后两个风电场的日出力期望值、方差和95%分位值偏差基本稳定在2%-4%。下表是这个过程的典型结果记录数据做了脱敏处理但数量级和真实项目接近校验指标原始1000场景削减后10场景偏差风电场1日均出力0.4120.4051.7%风电场2日均出力0.3860.3972.9%两场联合出力方差0.0210.0204.8%联合出力95%分位低出力侧0.1680.1614.2%有一年我在实际项目中反复用这套流程后来总结出一条心得无论采样多么精细、削减算法多么先进场景生成只是整个不确定性分析链路上的一个环节下游优化模型的形态、目标函数的凸性、约束条件的严格程度都会反过来影响你对场景质量的评判标准。不要试图一次性做出完美的场景集——先把整个链路跑通再回头针对薄弱环节做迭代优化这才是最省时的路线。另外分享一个最后的小技巧做完削减以后不要只把10条场景扔给优化模型就完事了。可以把原始1000条场景按削减得到的10条场景的“所属簇”画在一张图上看看每个簇内部的形态差异大不大。如果某个簇内部曲线形态特别分散说明这个簇的代表性不足需要增加场景数量或者换用更细腻的距离度量比如用“出力变化速率”加上“出力水平”一起定义距离。这套视觉化检查方法我一直在用比任何指标都直观。 时序场景生成与削减：从蒙特卡洛采样到相关性建模的完整实践
1. 为什么纯蒙特卡洛会在“时序相关”面前失灵先交代一个背景MCMonte Carlo蒙特卡洛方法做场景生成在电力系统、能源调度、金融风险这些领域里已经算常规操作了。思路也不复杂——对随机变量的概率分布做大量采样用成千上万条样本序列去逼近真实世界里可能出现的各种情况。光伏出力、风电功率、负荷曲线这些不确定性比较强的时序数据都能用MC方法生成一堆“如果明天是这样的天气/负荷”的场景集合。但这里有个很微妙的问题。如果你只是对每个时间点的分布单独采样然后把它们硬拼成一条时序曲线生成出来的场景在“形状”上会很怪一会儿冲到峰值、一会儿跌到谷底相邻时间点之间没有任何连贯性。真实世界不是这样的——风速不会从0瞬间跳到15m/s再瞬间掉回2m/s光伏出力也不会上午十点还是满发、十点零一分就变成零。相邻时刻之间是有惯性、有记忆效应的这个在专业术语里叫自相关两个风电场之间、或者风电和负荷之间也可能存在空间上的联动关系这叫互相关。两者合在一起就是标题里说的“时序相关性”。纯MC采样的致命缺陷就在这里。它只刻画了每个时间点各自的边缘分布完全没有刻画时间维度和空间维度上的相关结构。用这种场景去做随机优化算出来的调度策略会过于保守或者过于激进因为模型看到的“极端场景”很多在物理上根本不可能出现——同一时刻所有风电场都满发、下一时刻全部归零这种事情现实中不会发生。所以场景生成这个方向核心难点从来不是“怎么采样”而是“怎么把相关性结构正确嵌入到采样过程里”。后面我会先讲清楚生成端怎么做再讲削减端怎么做。削减也是个同样关键的话题MC采样动不动就是几千上万条场景直接塞进随机规划模型里根本算不动必须用算法把它们压缩到十几个、几十个有代表性的场景同时尽量不损失原始概率分布的信息。这两件事串起来才是完整的技术链路缺一环都没法落地。2. 场景生成实操从独立采样到相关性植入的完整链路2.1 第一步把样本频率定在“够用”的粒度上时序场景生成的第一步不是写采样代码而是先想清楚时间分辨率和场景长度。这两个参数直接决定后续所有计算量。以电力系统为例常见的做法是1小时一个点、取24个点代表一天或者15分钟一个点、取96个点。分辨率越细能刻画的波动细节越多但需要估计和存储的相关参数也会变多——96个时间点的相关矩阵就是96×96接近一万个参数样本量不够的话很容易估计出噪声来。我的建议是先想明白下游任务到底需要多细的粒度。如果是做日前调度1小时分辨率足够了如果是做日内滚动优化可能要15分钟甚至5分钟。不要一上来就追求最细的分辨率因为采样时的方差、削减时的计算量都是随着维数增长的。先用粗粒度把整个流程跑通确认相关性处理逻辑没问题再逐步加细。样本规模也有讲究。很多论文里生成2000条、5000条甚至10000条场景看起来越多越好但要注意MC的收敛速度是O(1/√N)想让精度提高一倍样本量得翻四倍。实际经验是如果最后要削减到10个场景初始生成500到1000条就够了再多对最终结果的影响非常有限纯粹浪费计算资源。这里有个更合理的做法是先估算下游模型对场景数目的敏感性用50、100、200条分别试跑看目标函数值变化幅度。变化不大就用小的变化明显就加大。后面我会专门讲这个。2.2 第二步用Cholesky分解把相关结构塞进白噪声处理时序相关性的经典手段是协方差矩阵Cholesky分解。它的原理是这样的假设你有m个随机变量比如m个风电场的出力它们之间的相关性用一个m×m的相关系数矩阵R来描述。你先从标准正态分布里独立采样得到一组互不相关的向量Z然后通过矩阵分解得到下三角矩阵L使得R L·Lᵀ用L去线性变换Z得到的新向量X L·Z就自动带上了R描述的相关系数结构。放到时序场景生成里的具体操作是把每个时间点当作一个变量先根据历史数据估计出整个时间跨度上的相关矩阵包括同一风电场不同时刻的自相关和不同风电场同一时刻的互相关然后对它做Cholesky分解再用分解得到的下三角矩阵去变换独立正态采样序列。核心代码如下import numpy as np # 假设有3个风电场24个时间点生成24*372维的相关场景 n_times 24 n_sites 3 n_total n_times * n_sites # 1. 根据历史数据估计相关矩阵简化示例用随机矩阵代替 # 实际使用中应该用历史出力序列直接计算 np.corrcoef R np.eye(n_total) # 这里应该替换为真实估计的相关矩阵 # 2. Cholesky分解 L np.linalg.cholesky(R) # 3. 采样独立标准正态随机向量 Z np.random.normal(size(n_total, n_samples)) # n_samples是场景数 # 4. 线性变换得到带相关性的样本 X L Z # 形状: (n_total, n_samples) # 5. 将结果重塑为 (n_sites, n_times, n_samples) scenarios X.reshape(n_sites, n_times, n_samples)这里面有一个工程上的关键点Cholesky分解要求相关矩阵必须是正定矩阵。但实际用历史数据估计出来的相关矩阵经常不是正定的——数据有缺失、时间序列长度不够、或者某些变量之间线性相关度极高都会导致矩阵半正定甚至不定。遇到这种情况直接做Cholesky分解会报错。我的处理办法有三个一是检查数据质量把缺失过多的时段剔除二是对估计出来的相关矩阵做特征值修正把所有小于某个阈值的特征值拉到一个正数三是改用特征分解替代Cholesky分解原理类似但容错性更好。具体怎么判断这个矩阵的病态程度我在后面问题排查那节会展开讲。2.3 第三步用Copula方法处理更复杂的尾部依赖Cholesky分解能处理线性相关性但真实世界的时序数据经常有更复杂的依赖结构。最典型的是风电出力——小风天和大风天都容易出现“一片区域的风电场同时不出力”或者“同时满发”的极端情况这种在统计上叫尾部依赖。用线性相关系数正态假设去刻画会低估极端场景发生的概率而电力系统最怕的恰恰就是极端的低出力高负荷场景。进阶做法是用Copula函数。它的核心思想是把“边缘分布”和“相关结构”分开建模先用历史数据分别拟合每个风电场出力的边缘分布可以用非参数核密度估计再选一个Copula函数来描述变量之间的依赖结构。常用的有Gaussian Copula和t-Copula——前者计算方便后者尾部依赖更强适合描述风电场出力这种有“同时极值”倾向的数据。具体操作流程是把历史数据通过各自的边缘分布变换成均匀分布序列再用这些均匀分布序列去拟合Copula参数最后从Copula里采样得到带上相关结构且服从真实边缘分布的样本。我在实际项目里的体会是如果数据量够用至少一两年的历史出力数据用t-Copula的效果会比纯Cholesky好不少尤其是在刻画极端场景分布形态方面。但代价是计算复杂度上去了参数估计也更敏感。所以我的建议是先用简单方法把流程跑通再逐步升级模型复杂度不要一上来就上Copula否则出了问题根本分不清是相关性建模错了还是边缘分布拟合错了。3. 场景削减到底在“削”什么从概率距离到代表性场景3.1 同步回代消除法的完整操作流程生成几千条场景以后不可能直接把原始样本丢进优化模型里因为计算量完全不可控。前面提过一句——8G内存的机器跑十万条场景的随机规划基本上会卡到怀疑人生。削减的本质是造一个规模更小的场景集使得新旧场景集之间的“概率分布距离”最小。这里的距离一般指Kantorovich距离通俗理解就是两个概率分布之间需要“搬运多少概率质量”才能变得一致。同步回代消除法Scenario Backward Reduction是目前用得最多的算法之一。它的核心逻辑很朴素反复找到“最没有存在感”的场景把它的概率叠加到离它最近的另一个场景上然后淘汰它。每一步都让场景总数减一直到数量满足要求。具体步骤拆开来看是四个环节为每个场景i找到它与所有其他场景j的最小距离min_dist(i)以及对应最近的场景编号在所有场景里找到min_dist最小的那个场景k它就是当前“最容易被替代”的场景把场景k的概率加到它的最近邻场景上然后从场景集合里删掉k更新删除后受影响的邻居关系重复1-3直到场景数达标这个算法的好处是每一步都保证场景集在Kantorovich距离意义上损失最小所以理论性质比较干净。实现上有一点需要注意距离的定义。更细致的做法是逐时间点计算欧氏距离再累加形成两条时序曲线之间的总距离。代码逻辑大概长这样def backward_reduction(scenarios, probs, target_num): scenarios: (n_scenarios, n_times) probs: 每个场景的概率和为1 target_num: 削减后的目标场景数 indices list(range(len(scenarios))) probs np.array(probs) while len(indices) target_num: # 计算所有场景两两之间的距离 dist_matrix np.zeros((len(indices), len(indices))) for i in range(len(indices)): for j in range(len(indices)): dist_matrix[i, j] np.linalg.norm( scenarios[indices[i]] - scenarios[indices[j]] ) np.fill_diagonal(dist_matrix, np.inf) # 找每个场景最近邻的距离 min_dist dist_matrix.min(axis1) nearest_idx dist_matrix.argmin(axis1) # 找到全局最小距离的场景 k min_dist.argmin() j nearest_idx[k] # 把概率合并给最近邻删除场景k probs[indices[j]] probs[indices[k]] indices.pop(k) new_scenarios scenarios[indices] new_probs probs[indices] new_probs new_probs / new_probs.sum() # 重新归一化 return new_scenarios, new_probs大厂和小研究团队在工程实现上的一个差别是距离矩阵的计算方式——用双重循环在场景数超过两三千时会有明显的性能瓶颈。我建议用scipy.spatial.distance.cdist或者直接矩阵广播运算能快一个数量级。3.2 聚类式削减与回代式削减的取舍同步回代消除法虽然经典但也不是唯一的选择。K-means和K-medoids聚类也经常被用于场景削减。聚类思路是把几千条场景按照相似度分成k组用每一组的质心均值曲线或中心点离所有样本最近的那条实际场景来代表这个组组的样本占比就是概率。两类方法放在一起对比各自的优劣比较明显。同步回代消除法得到的结果一定是原始场景里的某一条不会出现“平均出来但物理上不真实”的曲线这一点在新能源出力场景里很关键——把三条完全不同形态的出力曲线平均成一条可能就会得到一条现实中不可能出现的“中间态”曲线。我在实际做风电场景削减时特别在意这个因为调度模型对出力下限和上限很敏感平均值可能把顶峰削平、把低谷抬高这种场景看着合理其实会误导优化结果。K-means的优点是快非常适合场景数非常大比如上万条时做第一轮粗削减缺点是它对初始聚类中心敏感如果初始点选得不好容易收敛到局部最优导致削减后的场景多样性不足。我的常规做法是组合拳先用K-means把一万条粗削减到两百条再用同步回代消除法精削到十条以内。这样兼顾了速度和场景质量而且K-means那步产生的“质心平均”问题在第二轮同步回代时会被消除掉影响可控。削减规模怎么定也有经验可循。我见过很多文章不管三七二十一就削到10个但目标场景数和数据本身的结构复杂度是直接相关的。如果历史数据呈现出明显的季节性差异10个场景可能根本不够表达冬夏两种完全不同的出力形态。我的判断标准是削减之后按照新场景集重新算一遍关键统计量日均出力、最大出力、出力波动方差与原始场景集偏差不超过5%就可以接受。这个校验步骤比任何理论推导都实在。4. 实操中踩过的坑与排查清单4.1 相关矩阵不是正定矩阵Cholesky直接报错做Cholesky分解时最常遇到的报错就是相关矩阵非正定。我刚接触这个方向时踩过一次很深的坑用三四年历史出力数据估计相关系数中间有一段时间数据质量特别差很多时刻的出力值直接填了0结果算出来的相关矩阵行列式接近0一分解就崩。排查思路分两步走。先查数据——把历史出力序列画出来看看有没有异常的大段零值或者缺失值零值过多的时段会人为压低相关系数造成矩阵病态。我后来写了一套简单的数据清洗逻辑任何一天里出力为0的时刻超过全部时刻的30%这条样本直接丢掉宁可数据少一点也要保证相关结构是从真实物理过程里估计出来的。如果确认数据没问题矩阵仍然不正定就做特征值修正。思路是把相关矩阵做特征值分解把所有小于阈值的特征值推到一个正值再重新合成矩阵这一步在文献里叫“矩阵收缩”或“特征值裁剪”。还有一种更实用的方法是只截取跟风电场出力【原始描述】中“时序相关性”密切对应的变量子集去做Cholesky避免把所有变量无差别一锅炖相关性弱的变量会自动让矩阵变得很臃肿。我后来发现把96维的大矩阵拆成4个24维的小矩阵对各自时段内的相关结构分别建模问题会少很多。4.2 削减后概率分布“跑偏”了很多人做完削减只检查场景数对不对、曲线形不相似却忘了检查概率分布的统计特性。有一次我用K-means削完场景目测几条曲线都跟原始场景很接近但拿削减后的场景去跑机组组合模型结果目标函数值和全场景模型差了将近20%。排查下来发现是概率分配出了问题——K-means是均匀采样思想它默认每个簇的样本权重就是样本数占比但这种处理对极端场景不友好某些簇里只有一两条极端场景被K-means硬编码进了质心导致极端情况在削减后的场景集里被放大或丢失。后来我在这类场景削减任务中固定了一个审查习惯削减前后必须逐项比对这些统计量——各类场景的均值出力曲线、出力标准差、概率加权后的极端出力值比如5%分位点和95%分位点。发现问题就调削减方法比如改用K-medoids代替K-means或者用同步回代法。另外还有一个非常隐蔽的坑削减完以后概率一定要重新归一化。因为概率合并过程中舍入误差会积累不归一化的话后续算期望值时所有结果都会系统性偏小。4.3 场景数量到底取多少计算资源与精度的拉扯前面讲了场景数量是计算资源和精度之间的折中。我在做实际项目时给过一个比较通用的方法画出“场景数量-目标函数值”的收敛曲线。选一组目标场景数比如5、10、20、50、100分别跑同一个优化模型记录目标函数值的变化幅度。理论上场景越多结果越接近“真实随机优化”的极限值但这个边际收益是递减的。经验阈值是当场景数量翻倍但目标函数变化不到1%时就说明当前数量已经足够了。这个方法在小规模算例上很快就能跑出来大算例可以先削减到10个场景跑一次再把20个场景跑一次对比两次结果的差异。差异大就翻倍继续差异小就用少的。另外还要看你对“风险偏好”的态度求解的是期望成本最小化场景少点问题不大求解的是CVaR这类关注极端风险的指标场景少了根本覆盖不到尾部此时保守一点反而更稳妥。我个人在这个问题上吃过不少亏现在对极端风险类问题一般会比常规情况多留一倍的场景数量。5. 一个小型可复现案例两风电场出力场景从生成到削减理论和坑讲完了我拿一个完整的迷你案例把整条链路串起来。假设有两个风电场历史出力数据是按小时记录的一共365天每天24个点。目标是生成1000个日场景再削减到10个供一个简单的两机组经济调度模型使用。第一步是估计相关矩阵。把两个场各自24小时的数据拼接成48维向量用历史数据逐日计算相关系数得到48×48的相关矩阵R。注意这里的历史样本是365天等于365组48维样本样本量大约为维数的7倍勉强够用但不算富余。我用的是现场记录中的真实出力数据处理时要格外仔细地剔除检修日和极端停机日否则这些非典型样本会把相关结构拉偏。第二步是采样生成1000条带相关性的正态样本再通过历史数据拟合的逆累积分布函数转换成出力值。这里有个容易被忽视的操作风电出力是0到1之间的比例值直接从标准正态反变换会出现超出边界的情况我一般用截断正态分布做边缘分布或者用Beta分布拟合后再做概率积分变换。第三步是削减。按照前面说的方法先跑K-means粗削到100条再跑同步回代消除法精削到10条。原始1000条场景的总计算时间包括相关性采样、削减和统计校验在普通办公笔记本上大约是两三分钟完全可接受。最后做校验。我对比了削减前后两个风电场的日出力期望值、方差和95%分位值偏差基本稳定在2%-4%。下表是这个过程的典型结果记录数据做了脱敏处理但数量级和真实项目接近校验指标原始1000场景削减后10场景偏差风电场1日均出力0.4120.4051.7%风电场2日均出力0.3860.3972.9%两场联合出力方差0.0210.0204.8%联合出力95%分位低出力侧0.1680.1614.2%有一年我在实际项目中反复用这套流程后来总结出一条心得无论采样多么精细、削减算法多么先进场景生成只是整个不确定性分析链路上的一个环节下游优化模型的形态、目标函数的凸性、约束条件的严格程度都会反过来影响你对场景质量的评判标准。不要试图一次性做出完美的场景集——先把整个链路跑通再回头针对薄弱环节做迭代优化这才是最省时的路线。另外分享一个最后的小技巧做完削减以后不要只把10条场景扔给优化模型就完事了。可以把原始1000条场景按削减得到的10条场景的“所属簇”画在一张图上看看每个簇内部的形态差异大不大。如果某个簇内部曲线形态特别分散说明这个簇的代表性不足需要增加场景数量或者换用更细腻的距离度量比如用“出力变化速率”加上“出力水平”一起定义距离。这套视觉化检查方法我一直在用比任何指标都直观。