做电力系统随机优化的人,很难绕开“场景生成与场景削减”这两个词。当你在机组组合、经济调度或者电网规划模型里引入风光负荷不确定性时,多场景技术是当下跑得最成熟的一条路径:先用基于拉丁超立方抽样的方法生成成百上千条风、光、负荷时间序列,再用后向场景削减把这些序列压缩成十几个“典型日”,既保留不确定性信息,又让优化模型真正算得动。这篇文章是我在这类项目里的完整实操记录,从抽样原理讲到削减评估,顺手把踩过的坑也一并列出来,适合正在做风光负荷预测、随机优化和电力系统不确定性分析的同行参考。
1. 思路拆解:为什么场景生成必须分层抽样
先说结论:场景生成不是简单地从历史数据里复制几段曲线,而是要把“未来可能发生什么”用可计算的方式表达出来。风光负荷的波动性来自气象和用电行为,它们有随机性,也有统计规律。场景生成要做的就是抽取足够多样、分布合理的时间序列样本,让后面的优化模型既不会因为只看到一种情况而盲目乐观,也不会因为考虑了极端情况而过于保守。
1.1 拉丁超立方抽样到底解决了一个什么问题
最开始的初学者往往用直接随机抽样做这件事,然后发现结果很不稳定。普通蒙特卡洛随机抽样,相当于闭着眼睛往概率分布曲线上抓点,只要总样本量足够大,比如跑到上万条,分布自然会被摸清楚。但电力系统场景生成的工程场景里,样本量通常被压在2000以内,因为后面还有削减和优化计算,场景太多优化模型根本收敛不了。
样本量一旦降下来,普通随机抽样的缺陷就很明显:它可能让样本整体偏向分布的某一段,把低概率的“大风天”“连续阴雨天”“夏季午后的强负荷”漏掉,而这些恰好是影响系统安全性的关键场景。
拉丁超立方抽样做了个关键动作:把每个输入变量的概率区间平均切成N份,在每一份里强制取一个代表点,最后再把顺序随机打乱。它保证样本的累积分布函数和理论分布函数高度贴近,样本量能比蒙特卡洛随机抽样显著减少,同时覆盖分布的整个范围,尤其是尾部区域。一个直观类比:普通随机抽样像食堂师傅凭手感打菜,运气不好一勺全是肉一勺全是菜;拉丁超立方抽样则像个细心的营养师,把每道菜分成N份,保证每勺里荤素比例都均衡。
1.2 风光负荷三个输入量的建模起点
场景生成的目标不是一个单一的物理量,而是把风速、太阳辐照度、负荷预测误差同时作为三个输入变量来处理。风速在统计上常用两参数威布尔分布拟合,形状参数通常在1.5到3.5之间,尺度参数与平均风速有关。太阳辐照度由于受到云层遮挡等因素影响,通常用Beta分布描述,取值落在0到1之间,再乘上光照理论峰值得到实际辐照度。负荷的不确定性则更多来自预测误差,误差项通常近似服从正态分布,或者用t分布来加厚尾部。
这里有个容易搞错的细节:真正需要抽样的不是风速和负荷绝对值的分布,而是围绕预测值的“误差分布”。也就是说,先有确定性预测曲线,在它上面叠加从误差分布里抽出的扰动,这样生成出来的场景才符合实际调度场景。直接用历史绝对值的分布去抽样,场景会产生全局偏置,优化结果会出现虚假的“安全”或“紧张”。
1.3 相关性假设:场景生成里最容易丢的环节
拉丁超立方抽样是对每个变量分别分层抽样,但“分别抽样”不等于可以让它们互相独立。真实天气系统里,风速偏大的白天,云层往往被吹散,辐照度偏高;静稳天气下,辐照和风速都可能偏低。到了午后用电高峰,光伏出力和空调负荷又会叠加。如果抽样时不考虑这些相关性,生成的场景会过度乐观或过度保守,优化结果失真。
实操中,人们会在抽样之后加一个相关性调整步骤,最常用的是Iman-Conover方法:通过Cholesky分解构造带有目标相关系数的正态样本,再用排序的方式把目标相关结构“灌”进拉丁超立方采样矩阵里。这个步骤的成本很低,但效果非常明显,后面我会给出可直接复用的代码。
2. 场景生成实操:从数据准备到产出场景集
有了前面这些铺垫,下面进入具体环节。我以生成500个场景、每个场景覆盖24小时为例,走一遍完整流程。这套流程我在多个项目里用过,输入数据稍有差别时只需要调整分布参数和相关系数矩阵。
2.1 关键第一步:把历史数据拆解成可抽样的误差模型
拿到原始数据后,第一件事不是直接拟合分布,而是做残差分析。假设你手头有历史的风速预测序列、辐照度预测序列和负荷预测序列,以及对应的实际观测值。对每个时段,计算出相对误差或者绝对误差,把误差序列收集起来。
接着,把误差序列按时段分组,因为凌晨和下午的负荷预测误差结构差异很大,正午和傍晚的辐照度波动也不一样。常规做法是对每个时段分别拟合分布参数,得到24组误差分布,也可以按季节聚类后分簇拟合。风资源建议先对风速建模,再通过风机功率曲线把风速换算成出力;光伏则分“晴空辐射模型”和“随机波动项”两部分处理,这样比直接拟合功率分布更接近物理过程。
这一步虽然不起眼,但直接决定后续场景的质量。我见过不少项目跳过残差分析,直接对原始功率数据拟合,生成出来的场景均值偏移明显,削减后典型日怎么看都不对劲。
2.2 拉丁超立方抽样核心代码实现
下面这段代码实现拉丁超立方抽样的核心步骤,包括均匀分层和随机排列。之后用逆累积分布函数把均匀值映射到目标分布上。
import numpy as np from scipy.stats import weibull_min, beta, norm def lhs_sample(n_dim, n_samples, random_state=None): """生成 n_samples 行 n_dim 列的拉丁超立方均匀样本""" rng = np.random.default_rng(random_state) sample = np.zeros((n_samples, n_dim)) for d in range(n_dim): # 每个维度都把 [0,1] 区间切成 n_samples 份 edges = np.linspace(0.0, 1.0, n_samples + 1) for i in range(n_samples): # 在每个小区间内部均匀采样,保证全覆盖 sample[i, d] = rng.uniform(edges[i], edges[i + 1]) # 打乱该维度的顺序,避免各维度分层点固定对应 rng.shuffle(sample[:, d]) return sample # 用统一的均匀样本,分三列映射到风速、辐照度、负荷误差 unit_samples = lhs_sample(n_dim=3, n_samples=500, random_state=42) # 风速:威布尔分布,shape=2.2, scale=6.8 wind_speed = weibull_min.ppf(unit_samples[:, 0], c=2.2, scale=6.8) # 辐照度指数:Beta分布,参数取 alpha=2.0, beta=2.5 solar_index = beta.ppf(unit_samples[:, 1], a=2.0, b=2.5) # 负荷预测误差:均值为0,标准差为给定值 load_error = norm.ppf(unit_samples[:, 2], loc=0.0, scale=0.06)代码里的关键点,是均匀采样那一步用了小区间内部的随机数,而不是直接取区间中点。这样做保留了随机性,同时因为每个区间最多只出一个点,所以分布覆盖仍然有保证。
生成完整场景时,还需要做一步时间序列展开:对每个小时,用分布参数生成当日的随机变量,再叠加一个时变分量。比如风速场景可以写成基础风速乘以时变系数,再加一个自回归扰动;光伏场景则按太阳高度角算出理论辐照度曲线,再乘上抽样得到的辐照度指数;负荷场景在预测负荷曲线上叠加误差项。这样每个场景都是一条带日内波动形态的24小时序列,而不是一组互不相关的散点。
2.3 用Cholesky和排序恢复场景间相关性
如果直接把上面的三列抽样结果当作独立变量用,风速大的场景里辐照度仍然可能很大,这与实际气象规律相悖。处理办法是构造一个目标相关系数矩阵,然后重塑抽样结果的排列顺序。
假设目标相关系数矩阵如下:风速与辐照度相关系数为-0.35,风速与负荷误差为0.1,辐照度与负荷误差为0.45。实际操作时,这些数应当根据你所用地区的历史数据估算。
from scipy.stats import rankdata def apply_correlation(X, target_corr, random_state=None): """用 Iman-Conover 方法,按目标相关系数矩阵调整样本排列""" n, d = X.shape rng = np.random.default_rng(random_state) # 先造一个正态样本,再注入目标相关结构 Z = rng.standard_normal((n, d)) L = np.linalg.cholesky(target_corr) Z_corr = Z @ L.T # 对每一列,按 Z_corr 中该列的排名来重排 X 的对应列 X_ranked = np.empty_like(X) for j in range(d): ranks = rankdata(Z_corr[:, j], method='ordinal') - 1 X_ranked[:, j] = np.sort(X[:, j])[ranks] return X_ranked target_corr = np.array([ [1.0, -0.35, 0.10], [-0.35, 1.0, 0.45], [0.10, 0.45, 1.0] ]) scene = np.column_stack([wind_speed, solar_index, load_error]) scene_corr = apply_correlation(scene, target_corr, random_state=7)这样处理完之后,场景间的秩相关关系就能与目标矩阵基本一致。注意这里调整的是排列顺序,不改变每列边缘分布的值,所以不会破坏拉丁超立方分层带来的分布覆盖特性。这个细节是场景生成项目里最容易被忽略但又最影响结果质量的一环。
3. 后向场景削减:从500个场景压缩到十几个典型日
场景生成完毕,手里有了500条时间序列,每条还带一组概率权重。但这个体量对优化模型来说仍然过重,必须削减到几十条以内。后向削减就是一种从大集合出发,逐步删除冗余场景的方法。
3.1 削减的核心逻辑:不是随便删,而是最小化概率运输距离
要砍掉一批场景,最简单的方法是随机抽,或者按场景出现概率排序后保留概率最高的几个。这两种方法都不可靠,因为它们忽略了一个关键问题:场景之间的“距离”有多大。比如500个场景里可能有300个非常相似,都代表“多云微风晴热”,它们其实可以合并成几个典型场景;而“台风天”“寒潮夜”虽然出现概率低,离其他场景很远,一旦删掉,整个场景集的不确定性信息就丢了。
后向削减法由Growe-Kuska等人提出,核心是用Kantorovich距离来衡量两个场景集合之间的差异。通俗地讲,Kantorovich距离计算的是“把一个概率分布改造成另一个概率分布,最少需要搬运多少概率质量”。每个场景都有一个权重,当我们要删掉场景i时,它的概率质量要转移到剩余场景里,最优转移方式是转移到离它最近的另一个场景上。
于是后向削减的逻辑就很清晰了:每一轮找到那个“删掉它以后,整个场景集的距离增加最少”的场景,删除它,把它的权重加到离它最近的场景上,然后继续下一轮,直到场景数量达到目标值。这个过程保留了概率分布的总体形态,也保住了离群场景,因为离群场景一旦删除,距离代价会非常大,通常会被留到最后。
3.2 经典的迭代削减实现,附可直接改用的代码
下面给出一版简明的后向削减实现。它先算好场景间的距离矩阵,然后循环执行“找最小代价删除”和“权重转移”。
from scipy.spatial.distance import cdist def backward_reduction(scenarios, weights, n_target): """ scenarios: 三维数组,形状 [N, T, dims] weights: 初始概率权重,长度为 N n_target: 削减目标场景数 """ N = scenarios.shape[0] # 把场景展平成向量,计算两两欧氏距离 flat = scenarios.reshape(N, -1) dist = cdist(flat, flat, metric='euclidean') alive = np.ones(N, dtype=bool) w = weights.copy() while alive.sum() > n_target: idx = np.where(alive)[0] min_dist = np.full(len(idx), np.inf) for a, i in enumerate(idx): for b, j in enumerate(idx): if i != j and dist[i, j] < min_dist[a]: min_dist[a] = dist[i, j] # 删除代价 = 被删场景权重 * 到最近存活场景的距离 cost = w[idx] * min_dist rm_pos = np.argmin(cost) rm_idx = idx[rm_pos] # 找到离 rm_idx 最近的存活场景,承接被删场景的权重 neighbor_pos = np.argmin([dist[rm_idx, j] for j in idx if j != rm_idx]) neighbor_idx = idx[neighbor_pos] w[neighbor_idx] += w[rm_idx] alive[rm_idx] = False keep_idx = np.where(alive)[0] new_weights = w[keep_idx] / w[keep_idx].sum() return scenarios[keep_idx], new_weights, keep_idx使用时要留意距离度量的选择。如果场景里包含风速、辐照度、负荷多维数据,直接展平求欧氏距离会把三个量纲混在一起。建议先做归一化,或者分维度加权后再合并,比如给负荷误差乘一个常数比例,让它在距离计算中的权重符合你对不同风险的重视程度。另外,距离矩阵的计算复杂度是O(N²),N为2000时还可以接受,如果场景数上万,就先做一次快速预削减,比如用最近邻法先粗略去掉明显冗余的场景,再跑正式后向削减。
权重转移完成后,必须做一次归一化,把所有权重视为1。很多人在这一步忘了重新归一化,后面带入优化模型时,约束条件里的概率和会变成非1的值,结果全错。
3.3 削减完之后,怎么判断这批典型日能用了
削减不是结束,验证才是。一个合格的削减结果,要让削减后的场景集在统计指标上与原始场景集尽量接近,同时不能丢掉有工程意义的极值场景。
我通常做三组检查:
- 均值与标准差对比。把原始场景集和削减后场景集的逐小时均值、标准差画在同一张图里,看曲线是否贴合。如果削减后某个时段标准差明显变小,说明该时段的不确定性被低估了。
- 分位数对比。重点看5%、50%和95%分位数,尤其是风速和负荷的尾部区域。分位数差异超过10%时,需要检查距离度量是否合理。
- 极端场景保留情况。削减后场景里是否还包含“全年最大风速日”“光伏连续低出力日”。如果这些场景被削光了,优化模型算出来的方案可能在真实事故中失效。
为了更直观,这里放一张典型的指标对比逻辑:
| 指标 | 原始500场景 | 削减后15场景 | 说明 |
|---|---|---|---|
| 负荷均值偏差 | - | ≤1% | 偏差越小越好 |
| 风速标准差偏差 | - | ≤5% | 需要训练调参数 |
| 95%分位数偏差 | - | ≤8% | 关注尾部覆盖 |
| 极端场景个数 | 5个 | ≥2个 | 保留关键风险 |
削减数量也不是越少越好。我做过尝试,把15个场景压到5个,均值还贴着,但95%分位数系统性偏低,说明削减过头了。工程上常用的削减目标是10到20个场景,既能覆盖不确定性,又不会让随机优化问题陷入组合爆炸。
4. 踩坑实录与排查清单
这一节没有理论,全是项目实战里反复遇到过的具体问题。每个问题我都给出一线排查思路,算是这段时间下来最有价值的一部分。
4.1 高频问题速查表
| 问题现象 | 可能原因 | 解决办法 |
|---|---|---|
| 生成场景均值与历史均值明显偏移 | 没做残差分析,直接抽样绝对值 | 先构建误差分布,再叠加确定性预测 |
| 削减后的场景过度集中,缺少离群场景 | 距离度量的权重设置不合理,或削减目标数过小 | 归一化各维度,调整加权系数,适当增加场景数 |
| 风速与辐照度场景相关性不明显 | 抽样后没有做相关性调整 | 增加Iman-Conover步骤,按历史相关系数重排 |
| 负荷场景尾部过薄,优化结果过于乐观 | 负荷误差用了正态分布,未考虑厚尾 | 改成t分布或混合分布,加厚尾部 |
| 削减后场景权重出现负数或总和不为1 | 权重转移后没有重新归一化 | 每次迭代结束统一除以权重和 |
| 模型结果对削减后场景非常敏感 | 削减掉的某类关键场景恰好影响约束识别 | 检查极端场景是否保留,必要时人工锁存关键场景 |
4.2 实操心法:先小样本跑通,再大样本调优
第一次搭建整套流程时,别直接冲2000个场景、50万行数据。我先用200个场景、削减到10个场景快速跑通全链路,确认代码逻辑无误后,再把场景量拉上去。这样能极大节省调试时间。还有一个习惯是每次修改相关系数矩阵或分布参数前,先把上次生成结果存档,方便回滚对比。
后向削减的代码逻辑比较绕,但它本质上就是一个“删除 + 权重转移”的循环。调试时可以在循环里打印本轮删除的场景编号、删除代价和最近邻居编号,看到哪一步出现了不合逻辑的转移,问题通常就集中在距离矩阵计算那几行。
4.3 关于场景削减后优化收敛的一点个人经验
场景削减后接入随机优化模型时,容易被忽视的是“场景数量”与“优化变量规模”的联动效应。即使削减到了15个场景,机组组合问题里二进制变量也会成倍增加,求解时间可能从分钟级跳到小时级。一个有效做法是先用削减到5个场景的模型做快速预演,确认约束和边界条件没有明显冲突,再切换到完整15个场景做正式求解。这样跑出来的结果,既保留了削减精度,又不会让求解器卡死在你服务器上。
另外,我习惯在削减后人工检查一两个关键场景的时间序列曲线,用肉眼看它们是不是“正常天气”。有时候统计指标全部合格,画出来却发现场景出现了风速夜间突增、负荷时段错位这类物理上不可能的情况。这类问题往往出在残差模型没有考虑时段相关性,及时发现能省很多弯路。
这套方法不只是风电场和光伏电站规划能用。凡是碰到“不确定性输入过多、必须压缩成有限场景”的随机规划问题,拉丁超立方抽样加后向削减这个组合都值得试一遍。我最近还在想,能不能把削减环节做得更智能一些,直接根据下游优化模型的目标函数来动态决定保留哪些场景,而不只是做纯统计距离削减,这个方向后续项目里会继续推。先把眼前这套流程跑稳,比追任何新概念都管用。