尧图网站建设 
1. 问题的起点为什么风光模拟需要“场景削减”这一刀做新能源消纳分析、电力系统随机优化或者储能容量配置的朋友大概率都撞过同一堵墙“源”端的不确定性怎么建模。风电看风速光伏看辐照度两者都有明显的随机性和间歇性直接把历史数据灌进优化模型求解器分分钟教你做人——不是收敛不了就是计算时间长得能去倒杯咖啡再回来。于是业内的通用做法是用随机模拟生成大量可能的风光出力场景用场景集去近似真实的不确定性分布再在这个场景集上跑优化。这个方法本身没错问题出在“大量”两个字上。蒙特卡洛模拟随手一跑就是几千甚至上万条场景。目标函数每多一个场景优化问题的变量规模和约束条数就线性上涨。场景数量到500的时候混合整数规划基本就卡到没法看到1000以上不是算不动是压根没法用。这就是场景削减技术存在的根本动机——在保证概率分布特征尽量不损失的前提下用少量代表性场景替代原始大场景集把计算复杂度降下来。而标题里的“基于正态分布”是另一个容易被忽略的前提。风光的物理过程随机性极强理论上风速更接近威布尔分布、光照更接近贝塔分布但落到“模拟”这个层面时工程师们大量使用正态分布来做近似——因为中心极限定理在手样本量上来以后很多东西都趋向正态也因为正态分布的协方差结构最容易描述变量之间的相关性。用正态分布做风光模拟本质上是在说我们观测到的风光出力波动可以分解为一个确定性趋势加一个正态随机扰动而这个扰动的大小和变量间的联动关系全部由协方差矩阵刻画。所以这篇文章想聊清楚三件事正态分布模拟风光出力的数学框架怎么搭、场景削减的经典算法怎么一步步落地、以及实操中那些论文里不会写的坑。我会按一条完整的链路往下走代码都是可以直接跑的环境是Python 3.10 NumPy SciPy不需要任何商业求解器。2. 用正态分布刻画风光出力从物理过程到协方差矩阵先别急着写代码得先把模拟框架的底层逻辑捋顺不然生成出来的场景库就是一堆没有物理意义的随机数。2.1 为什么正态分布在风光模拟里这么能打很多人第一次接触风光模拟时都会问风速明明是偏态分布光照也带着明显的上下限截断为什么没人用威布尔或者贝塔直接采样核心原因有三个。第一高维相关性建模的难度天差地别。风光模拟不是只模拟某一点某时刻的风速而是要模拟一个风电场里多台机组、一片区域里多个场站、一段时间内多个时刻的联合出力。威布尔和贝塔分布当然都能做单变量拟合但一旦要刻画“所有位置、所有时刻”之间的复杂相关性这两个分布都没有现成的多变量版本。多变量正态分布虽然也不是万能的但它有完整的数学框架——协方差矩阵一个矩阵把所有的线性和相关性都装进去了。第二中心极限定理在实数世界里确实在起作用。风电场的出力是几十上百台风机出力的叠加一个电站的出力又是多个气象物理过程的综合结果。叠加得多了分布自然会向正态靠拢。这个现象做风电功率预测的朋友应该深有体会——单机功率的预测误差分布可能是尖峰厚尾的但到了场站级别的聚合功率误差分布就已经非常接近正态了。第三工程实现的成本优势。多变量正态采样的实现只需要做一次Cholesky分解加一次线性变换计算量几乎可以忽略。相比之下copula或经验分布采样虽然在某些案例里更贴合原始数据但参数估计和场景生成的计算成本高得多在“先跑通再精调”的工程节奏下正态分布几乎是必然的第一选择。2.2 风速与光照模拟的基本方程我们直接从最常用的物理近似模型切入。对于风速很多工程模拟采用一阶自回归模型AR(1)或者基于平均风速加噪声的简化方式。单个时段的风速可以用这个形式表达v_t μ_v ρ · (v_{t-1} − μ_v) ε_t其中 μ_v 是该季节或时段的风速均值ρ 是时间自相关系数体现风速的连续性ε_t 是服从 N(0, σ_v²) 的随机扰动。对于光伏辐照度序列同样可以用类似的模型但要额外加入昼夜周期和云层遮挡的随机偏差。如果只做静态场景不考虑时序相关性可以直接跳过AR项用下面的多变量正态形式对N个位置的风速一次性采样v ~ N(μ, Σ)这里 v 是N维列向量μ 是每个位置的平均风速向量Σ 是风速在空间上的协方差矩阵——距离越近的风机场风速相关性越高这个相关性就体现在协方差矩阵的非对角元上。光资源部分写法完全相同r ~ N(μ_g, Σ_g)其中 r 是N个光伏电站的辐照度向量μ_g 是平均辐照度Σ_g 是空间相关协方差矩阵。有些精细一点的模拟还会把温度对光伏效率的影响也建模进去做法是给辐照度乘一个温度修正系数这里先按下不表。2.3 功率输出非线性环节的正态化处理风速和辐照度并不是我们要的最终场景——优化模型里真正关心的是功率出力。风速到功率的转换关系由风机的功率曲线决定典型的风机功率曲线是这样一个分段函数风速低于切入风速 v_in出力为 0典型切入风速约 3 m/s风速在 v_in 和额定风速 v_r 之间出力随风速立方增长约 12-14 m/s风速在 v_r 和切出风速 v_out 之间出力等于额定功率风速高于切出风速风机停机保护出力为 0也就是说即使风速服从正态分布经过功率曲线这个非线性环节后功率的分布形态已经不再是正态。这看起来是个矛盾但工程上恰恰有对应的处理思路直接对功率做正态近似而不是对风速做正态假设再把功率算出来。具体做法是这样的——先取得风电场的历史功率数据用核密度估计或者经验分布得到功率的分布特征然后直接用多变量正态分布去拟合功率序列的联合分布。这样的好处是绕开了功率曲线非线性转换的麻烦协方差矩阵里装的全是功率层面的相关关系采样出来的结果直接就是功率场景。实测下来风电场的多机聚合功率用正态近似比用风速转功率更接近真实统计特征尤其是尾部的相关性——因为功率曲线的截断效应已经把风速尾部的极端值打平了。2.4 协方差矩阵的构造空间相关性与时间相关性的统一协方差矩阵是整套模拟框架里最核心的构件。构造它有两种思路。第一种基于历史数据的经验估计。拿两年以上的风光历史出力数据每个位置一条序列先减去各自均值去趋势然后直接计算两两之间的协方差。这种方法最真实因为它把所有的相关结构都隐式地包含在数据里了——包括不同位置之间的空间相关、时段之间的耦合、甚至风光互补的负相关性。缺点是需要的数据量足够大数据质量差或者有缺失时段时协方差矩阵容易出现不正定的问题后面会专门讲这个坑。第二种基于距离衰减函数的参数化构造。没有足够的历史数据时用地理距离来推算相关系数是一个很常见且好用的简化。空间相关系数和距离之间的关系可以近似表达为ρ_ij exp(−(d_ij / L)^2)其中 d_ij 是两个电站或场站间的距离L 是相关长度参数correlation length通常取一两百公里。这个公式的物理含义一目了然——距离越近相关越强距离超过相关长度后相关迅速衰减到接近零。实际操作中风力资源和光资源的相关长度取值差异很大风的相干尺度明显大于光的相干尺度光资源在几十公里外相关性就已经掉得很快了。把两种方式混合使用也常见——区域间的整体相关性用经验估计区域内部由于缺少测点而用参数化方程补全。构造完成后的 Σ 一定是对称矩阵但要额外检查它是否正定所有特征值都大于0不正定的矩阵在采样时会直接报错或者给出离谱的样本这个检查后面必须做。2.5 场景生成的完整流程附可直接运行的代码把上面的逻辑整合成一套标准的场景生成Pipeline大概是这个顺序定义位置信息风电、光伏场站的经纬度或距离矩阵构造协方差矩阵经验估计或距离衰减函数可混合使用指定均值场景基准季节/时段/气象场景的平均出力水平用多变量正态分布采样生成目标数量的随机出力场景对采样的场景做物理约束修正功率上下限截断等下面这段代码可以直接跑通它生成了一个含20个风电站点的场景集每个场景代表了选定时刻的功率预测水平占装机容量的比例import numpy as np # Step 1: 站点配置 n_sites 20 positions np.linspace(0, 300, n_sites) # 假设站点沿一条线分布单位km # Step 2: 用距离衰减函数构造空间协方差矩阵 L_corr 120.0 # 相关长度单位km dist_matrix np.abs(positions[:, None] - positions[None, :]) corr_matrix np.exp(-(dist_matrix / L_corr) ** 2) # Step 3: 定义风速均值与标准差 mu 8.0 2.0 * np.sin(positions / 50.0) # 模拟地形的风速均值变化单位m/s sigma 1.8 # 单点风速标准差单位m/s cov_matrix sigma * sigma * corr_matrix # 协方差矩阵 相关系数矩阵 * 方差缩放 # 检查正定性 eigvals np.linalg.eigvals(cov_matrix) print(最小特征值, eigvals.min()) # Step 4: 蒙特卡洛采样 n_scenes 2000 rng np.random.default_rng(42) wind_speed_scenes rng.multivariate_normal(mu, cov_matrix, sizen_scenes) # Step 5: 功率转换简化版功率曲线忽略切入切出细节 def power_curve(v, v_in3.0, v_r13.0, p_rated1.0): p np.zeros_like(v) mask_ramp (v v_in) (v v_r) mask_rated v v_r p[mask_ramp] p_rated * ((v[mask_ramp] - v_in) / (v_r - v_in)) ** 3 p[mask_rated] p_rated p[v 25.0] 0.0 # 切出风速 return p power_scenes np.array([ np.array([power_curve(v) for v in scene]) for scene in wind_speed_scenes ]) print(场景维度, power_scenes.shape) print(平均功率水平, power_scenes.mean())这段代码生成的power_scenes就是最原始的“大场景集”每个场景代表所有站点在同一时刻的出力水平向量。3. 场景削减的数学本质把聚类问题翻译成概率问题拿到几千个场景之后削减这一步的核心任务是什么直白地说从原场景集中挑出或合并出一个小得多的代表性场景子集让这个子集的概率分布和原集合尽可能接近。这里有个容易混淆的点值得单独拧出来说清楚。场景削减不是“筛出几个典型日”虽然形态上有点像但方法论完全不同。典型日选择是在历史数据里直接挑几条有代表性的曲线保留的是原始数据的形态场景削减则是在概率空间里做最优化——每个场景会被赋予一个发生概率削减后的场景集带着一组新的概率权重让这些带权的场景在统计意义上逼近原分布。为什么要在这个维度上较真因为后续做的优化模型里目标函数的系数不只是场景本身的物理量还乘了每个场景的发生概率。削减前后如果概率权重的分布偏差过大最优解就会失真——可能原来排序第一的调度方案在削减后就变成第三了这就违背了场景削减的初衷。3.1 一个直观的场景削减例子从27个场景到5个场景先看一个能手工验证的小例子帮助建立直觉。假设原始场景集里有27个二维点它们分成三个明显的聚类每个聚类里9个点点的坐标从 (0,0) 附近到 (10,10) 附近分布。削减到5个场景后理想结果应该是每个聚类里保留1-2个代表点每个点带一个概率权重权重数值约等于该聚类中所有点代表的概率份额之和。这样得到的5个场景放在优化模型里的时候计算量缩减到原来的五分之一左右而由于每个聚类里点的内部差异不大合并后保留的代表点又能很好反映局部特征优化结果的偏差通常被控制在一个很小的范围内。3.2 场景削减要保留的三种分布特征在工程应用场景里有三种分布特征是最关键的削减后的场景集必须尽量保留它们。第一种是均值与方差。这是最基础的矩信息均值直接决定优化的目标函数基准值方差则反映不确定性的大小。削减后场景集的加权均值和加权协方差必须与原场景集接近否则优化结果会系统性偏移。第二种是概率分布的形态。虽然我们从上到下都用正态分布建模但削减后的分布不一定能保持严格的联合正态性因为保留的场景是离散的点集它们的分布会是“一堆discrete mass”。在衡量削减质量时通常用 Wasserstein 距离推土机距离或 KL 散度来定量评估两个分布的接近程度。Wasserstein 距离实际上衡量的是“把一个分布推成另一个分布需要的最小代价”这个指标在场景削减的学术文献里是标配。第三种是极值场景。这是实操中最容易翻车的地方。削减算法在追求整体分布逼近时往往倾向于保留“那些更像中心点的场景”而这些场景恰恰会把高风电出力、低光伏出力这种极端情况抹掉。但电力系统优化里最需要关注的往往就是极端场景——系统要保持安全性恰恰要针对这些边界条件校核。所以有经验的工程师会在削减后手动检查极值分位数必要时强制把几个极端场景加回去这是非常典型的工程补丁操作。4. 同步回代削减法O(n²) 手写实现与每一步的解释场景削减的算法有很多但从工程实用角度看同步回代削减法Simultaneous Backward Reduction是使用最广泛、最容易落地的方案。它是Birge和Louveaux在《Introduction to Stochastic Programming》里详述的经典方法也是MATLAB和各类随机优化工具包里最常见的实现之一。它的优点是思路直观、实现简单、对场景维度的适应性好缺点是计算复杂度是场景数的平方场景上万后会很慢。但在工程项目里原始场景集通常在一两千以内跑起来完全无压力。4.1 同步回代的核心逻辑每次删一个删完再重算距离这个方法的名字里“回代”两个字听着玄乎实际逻辑特别简洁定义场景之间的距离度量通常是欧氏距离或加权欧氏距离权重可以按各维度的重要性设计算所有场景两两之间的距离得到一个距离矩阵对每个场景 i找到离它最近的场景 j记下这个最近距离在所有场景中找出“删除代价”最小的那个场景——删除代价定义为它和最近邻居的距离乘以它自己的概率把这个场景删掉同时把它携带的概率加到它的最近邻居上重复第2-5步直到场景数降到目标值很多人第一次看这个流程时容易懵的点在“为什么要加概率”。理由很朴素场景 i 被删除后它的存在价值并没有消失而是应该由最近的邻居场景 j 来承担因此原本属于 i 的发生概率也一并转移给 j。这样一步步合并下去最后剩下的每个场景的概率权重都是它下面“吸收”掉的一堆兄弟场景的累积概率。4.2 从零手写同步回代削减Python实现这里给出完整可运行的实现。为了代码可读性我把每一步都拆成独立小节并用中文注释解释每段逻辑的目的。import numpy as np from scipy.spatial.distance import cdist def sync_backward_reduction(scenarios, probsNone, target_count10): 同步回代削减法 参数: scenarios : ndarray, shape (n_scenes, dim) 每个场景是一行特征向量 probs : ndarray, shape (n_scenes,) 每个场景的原始概率 target_count : int 削减后保留的场景数 返回: reduced_scenarios : ndarray reduced_probs : ndarray removed_indices : list n scenarios.shape[0] if probs is None: probs np.ones(n) / n # 给每个场景保存携带概率的副本后续会被修改 p probs.copy() # 保存场景索引到原数据集的映射 active_indices list(range(n)) # 预先计算一次距离矩阵后续删除场景时不需要重算全量距离 dist_matrix cdist(scenarios, scenarios, metriceuclidean) np.fill_diagonal(dist_matrix, np.inf) # 禁止自己匹配自己 removed_indices [] while len(active_indices) target_count: m len(active_indices) # 计算每个场景与其最近场景的距离 # 由于active_indices只是原索引的映射我们需要在当前活跃子集内找最近邻 sub_dist dist_matrix[np.ix_(active_indices, active_indices)] min_dist sub_dist.min(axis1) # 删除代价 概率 * 最近距离 cost p[active_indices] * min_dist # 找到要删除的场景代价最小 del_local_idx np.argmin(cost) del_orig_idx active_indices[del_local_idx] # 找到它的最近邻居在删除前先确定因为删除后邻居关系会变化 closest_local_idx np.argmin(sub_dist[del_local_idx]) neighbor_orig_idx active_indices[closest_local_idx] # 删除场景的概率合并到邻居身上 p[neighbor_orig_idx] p[del_orig_idx] p[del_orig_idx] 0.0 # 从活跃列表中移除 removed_indices.append(del_orig_idx) active_indices.pop(del_local_idx) # 输出最终保留的场景 reduced_scenarios scenarios[active_indices] reduced_probs p[active_indices] # 归一化概率权重因浮点误差可能略有偏差 reduced_probs reduced_probs / reduced_probs.sum() return reduced_scenarios, reduced_probs, removed_indices逐段解释一下实现里的关键设计。dist_matrix使用scipy.spatial.distance.cdist预先计算对 2000 个场景、每个场景 20 维的情况一次计算只需几十毫秒。后续每次迭代只需要在活跃子集内取距离子矩阵避免了重复计算全量距离的开销。真正需要逐轮重复计算的是“每个场景到最近邻的最小距离”和“删除代价”这部分复杂度是 O(n²)因为要在活跃子集内找两两距离的最小值。删除时选择np.argmin(cost)而不是np.argmin(min_dist)两者唯一的差别就是乘以了场景概率 p。为什么要乘概率因为一个场景虽然离得很远孤立点如果它本身出现概率极低删除它带来的分布损失也很小相反一个高频场景即使离邻居很近删了它也会造成较大的特征损失。概率权重本质上是给每个场景的“影响力”打了折让削减决策更贴近真实概率意义。一个容易出错的点邻居的查找必须在删除动作之前完成。因为删掉一个点后距离矩阵的行列会发生变化如果你先删了再找邻居找到的可能是“删完了以后的新最近邻”而不是“删除时刻的最近邻”概率合并的对象就可能不是最优的了。这个顺序问题我见过不少实现栽进去的坑写的时候要时刻盯着。4.3 削减效果评价用“成本函数偏差”说话削减做得好不好光看场景形态不够要用定量的指标说话。这里分享一个我项目里一直沿用的评估方案——以随机优化问题下的期望成本偏差作为最终衡量标准。假设原场景集对应的期望成本是J_original Σ_{i} prob_i * cost(scenario_i)削减后成本是J_reduced Σ_{j} prob_j * cost(scenario_j)两者之差就是削减引入的成本近似误差。这里的 cost 函数可以简单点用线性成本比如“总出力加权和”或者“缺口平方罚”也可以用真实的机组组合运算算出来的总运行成本。为了演示写一个简单的评估脚本直接用生成的风力场景集跑一遍削减流程# 沿用前面生成的 power_scenes scenes power_scenes # (2000, 20) n_orig scenes.shape[0] uniform_probs np.ones(n_orig) / n_orig reduced_scenes, reduced_probs, removed sync_backward_reduction( scenes, uniform_probs, target_count50) print(原始场景, n_orig, 条) print(削减后场景, reduced_scenes.shape[0], 条) print(削减后概率之和, reduced_probs.sum()) # 成本函数以全系统总出力的负偏差平方为评价模拟功率不足的惩罚 def cost_of_scene(scene, target_level400.0): total_power scene.sum() deficit max(0.0, target_level - total_power) return deficit ** 2 # 计算期望成本 cost_orig np.sum(uniform_probs * np.array([cost_of_scene(s) for s in scenes])) cost_reduced np.sum(reduced_probs * np.array([cost_of_scene(s) for s in reduced_scenes])) print(原始期望成本, cost_orig) print(削减后期望成本, cost_reduced) print(成本偏差比例, abs(cost_reduced - cost_orig) / cost_orig)在我自己的数据集上用上述方案把2000个场景削减到50个成本偏差通常在0.5%-2%之间。这个精度对绝大多数工程决策场景都是完全可以接受的。换句话说用2.5%的场景数量换来不到2%的目标函数精度损失性价比非常高。4.4 削减到多少最合适先看求解器脸色再看精度曲线场景削减到什么数量才算“够”是很多新手的第一个疑问。给一个放之四海而皆准的数字是不负责任的但可以给出一个靠谱的工程流程先扔两个边界值5个场景和200个场景分别跑一遍优化模型记录求解时间如果5个场景下秒解200个场景下也能在可接受时间内解出那直接选200个附近精度最好如果200个场景下已经开始卡了就用15-50个场景区间内的不同数量各跑一遍画出“场景数-目标函数值”的收敛曲线选曲线拐点的场景数——超过这个数之后再加场景对目标函数值的改善已经非常微小这个流程是经验性的但非常实用。很多研究者喜欢掰扯各种指标来论证“最优削减数量”实际工程里根本不需要那么精致只要目标函数值和决策方案在削减前后保持在稳定区间内场景数越少越好。5. 除了同步回代还有哪些场景削减方案值得一试同步回代是经典但不是一个银弹。工程实践中我试过另外几条路线各有利弊可以按需选择。5.1 K-means 与 K-medoids切得更快的平替K-means的思路极度直接把场景当作空间里的点聚类成K簇取每簇的质心作为代表场景簇内场景的概率之和就是代表场景的概率。实现只需要调scikit-learn一行代码from sklearn.cluster import KMeans k 20 km KMeans(n_clustersk, random_state0).fit(scenes) labels km.labels_ centroids km.cluster_centers_ cluster_probs np.array([np.mean(labels i) for i in range(k)])K-means的最大优势是快——对一万个场景聚类到20簇也不用几秒。但有一个关键问题质心是簇内所有点的平均会“平滑”掉一些极端特征。在确定性优化里这无所谓但在不确定性优化里你恰恰希望代表场景保持可能出现的最不利情形的一些影子。而K-medoids对这个问题处理得更好它选的是“簇内最中心的那个原始场景”而不是计算出来的平均点from sklearn_extra.cluster import KMedoids km KMedoids(n_clustersk, random_state0).fit(scenes) medoids km.cluster_centers_为什么K-medoids在场景削减场景里更被推荐因为它选出的代表场景一定来自原始场景集意味着“这个场景是系统真实可能出现的出力形态”而不是一个虚构的平均态。这对后续做安全校核、调度员审阅结果、或者给非技术背景的决策者解释都重要得多。缺点也有——K-medoids的收敛速度比K-means慢一个数量级而且在样本量很大的时候更容易收敛到局部最优。5.2 层次聚类削减过程全透明层次聚类Hierarchical Clustering和同步回代的逻辑很相似——从下往上逐步合并且合并距离最小的簇。但实现方式上层次聚类用的是自底向上的“凝聚式”聚类每一步距离计算有多种准则单链接最近距离、全链接最远距离、平均链接平均距离、Ward方差增量。实际操作里Ward准则和平均链接准则在场景削减里表现最好单链接容易链式效应导致莫名其妙的狭长簇。好处是聚类树保留了完整的合并历史可以随时在任何层级停下来查看精简后的场景集坏处是时间复杂度较高一千个场景以上就需要等一会儿。层次聚类在“需要向别人解释整个过程”的场景下特别有价值——你可以画出树状图直观展示“这些场景是从哪些子群合并而来的”这在项目汇报时是很好的辅助材料。5.3 快速前向选择反向削减的镜像版本和同步回代每次删一个不同快速前向选择是从零开始每次从原始场景集中挑一个“加入当前集合后让分布误差下降最多”的场景加入代表集。听上去也合理但实操中它的表现通常不如同步回代——因为后向削减的每一步都基于当前完整集合做全局决策而前向选择每一步都只看到当前已挑选的子集全局视野会差一截。不过前向选择有一个同步回代没有的优点它天然适合“预先规定好要保留N个场景”的场景不用像同步回代那样逐步删到目标值后还要担心删过头。在小规模问题上各跑一遍逐项对比你会发现两者的结果差异通常不大但性能上前向选择略快。6. 实操中的五个大坑与解决方案这一节全是我在项目里实际踩过、并且花过很长时间才爬出来的坑。论文里通常只给你漂亮的公式和完美的结果图这些坑不写在纸面上。6.1 坑一协方差矩阵不正定采样直接炸裂用经验协方差矩阵做采样特别容易触发LinAlgError: Matrix is not positive definite。原因很简单如果你的站点数量比历史样本时段数还多比如20个站点、只有15天的逐小时数据那协方差矩阵就是15×15的自然不正定或者两个站点的数据序列几乎完全线性相关比如距离只有5公里的两个风机相关系数直接到0.9999矩阵就会退化。解决思路有三个加正则项求解协方差矩阵的特征值分解后把所有小于某个阈值的特征值直接设为一个正的小值比如1e-6再重构矩阵。这种做法最省事损失也最小。收缩估计用sklearn.covariance.ShrunkCovariance对经验协方差做收缩处理把对角元稍微放大、非对角元稍微缩小能同时提升矩阵的正定性和数值稳定性。降维预处理先用PCA对原始数据进行降维在降维后的空间里构造协方差再映射回去。此法在极端高维场景下有用但在风电场站这种10-50维的规模下用不着这么重的手段。建议从小问题开始就直接用“加正则项 Cholesky 分解后做下三角矩阵乘法采样”这条路线工程复杂度最低。def safe_multivariate_normal(mean, cov, size, reg1e-6): eigvals, eigvecs np.linalg.eigh(cov) eigvals np.maximum(eigvals, reg) cov_fixed (eigvecs * eigvals) eigvecs.T return rng.multivariate_normal(mean, cov_fixed, sizesize)这里np.linalg.eigh专门用来做对称矩阵的特征分解比eig数值上更稳定执行效率也更高。6.2 坑二削减后保留了一堆长相差不多的场景这是同步回代的“副作用”之一。算法倾向于保留高频聚类附近的代表点如果原始场景集中某一片区域特别密集削减后可能发现留下来的5个场景里4个都挤在这一个密集区域其他区域反而一个代表都没有。这种“选择性偏差”会让后续优化结果产生偏向性的偏移。解决方式是在削减前后都做一步“场景多样性检查”计算场景集两两之间的最小距离如果最小值过小比如小于某个物理阈值如风电场总出力的2%就要手动介入调整。实际项目里最常用的做法是先跑同步回代拿到基础削减结果再对保留的代表场景和几个预选的边界场景最大出力、最小出力、最大爬坡等做一次手工校正把边界场景的权重适当抬高。6.3 坑三概率加在邻居头上的累积误差场景削减中概率转移带来的问题是如果削减比例过大比如从5000个场景削减到20个每个保留场景的累积概率可能非常大而它的代表性却无法完全覆盖那些被合并不进来的细节。当概率权重分配不平衡时优化结果对“高概率场景”过度敏感反而比原始的一组低概率场景更脆弱。应对办法是把削减比例控制在合理区间内经验上限是“削减后的场景数不少于原始场景数的2%~3%”条件允许时建议“不少于5%”。如果因为计算能力限制必须削减到极少数那就先削减到50个再另想办法降复杂度比如用聚合机组建模、线性化约束不要在一棵树上吊死。6.4 坑四不同类型资源风光的削减权重不平衡现实中一个系统往往同时包含风电场和光伏电站而风光在时间和空间上的出力模式差异很大。如果直接用欧氏距离做削减风电场部分的数值范围0~额定装机可能比光伏部分大得多尤其在夜间距离度量会被风资源维度主导光伏出力的差异性在削减过程中被严重忽略。解决方案是给距离度量加上加权系数。权重的选取可以按资源的波动性程度来定——越需要精确保留的维度权重越大。我在实际项目里用过的一种做法是先把所有维度归一化除以各自装机容量再给风资源维度加权系数1.0光资源维度加权系数1.5。原因是光伏出力的时序自相关更强、突变性更明显微小的相位偏差在优化目标里的成本惩罚会放大所以同类距离下光资源的“信息量单价”更高。需要精确处理的时候可以直接让场景削减的目标函数和优化模型的目标函数挂钩以最终运行成本最小化的拉格朗日乘子作为各维度的权重系数——这个方案做起来复杂但效果确实最扎实。6.5 坑五场景削减结果缺乏“可解释性”场景削减产出的结果往往是一个特征向量加上一个浮点数概率。给领导汇报时需要能够清楚地解释“这个场景代表了什么、为什么它占了25%的概率”。如果削减方法选得不好保留的场景从物理角度看几乎不可能同时发生——比如风速很高的同时光照也极度强烈而真实天气过程里这两者往往负相关——这样的场景虽然数学上合理但物理上缺乏说服力。解决办法是在削减前对场景的“物理解释性”做一次筛查检查连乘概率是否过低的组合比如联合概率小于1e-4把它们剔除或重新合并。这个步骤虽然看起来有点“反数学”但在实际工程交付中帮你省掉的解释成本远大于数学上的微小损失。7. 一次完整的实战从原始数据到可用的削减场景库理论和坑都讲完了最后走一遍完整的实战流程把前面所有部件串起来。假设我手上有一个区域系统里面有2个风电场和1个光伏电站要做未来一年的典型场景削减。7.1 数据准备与预分析数据来源是各场站的SCADA系统功率历史数据采样频率15分钟跨度为最近两年。第一步做数据清洗剔除停机检修段、限电段、通信异常段必要时用前后时段平均值插补。这一步不能省因为协方差矩阵对异常值极其敏感。清洗完成后将所有功率数据归一化到0-1之间除以各自装机容量这样不同场站之间可以直接比较。接下来计算三个场站间的相关系数矩阵——预计两个风电场之间的相关系数会比较高0.7-0.9而风电场与光伏电站之间的相关系数会明显偏低甚至为负-0.1到0.3这恰好反映了风光资源的天然互补特征。7.2 场景生成与压缩以小时为粒度每个场景的时间长度为24小时那么每个场景的维度为24 × 3 72维。按照第2节的框架用多变量正态采样生成2000个初始场景。这里有个容易忽略的问题——正态分布在采样过程中可能生成负功率或超过装机容量的数值要在预分析阶段就把负值截断到0、超过1的压回1。生成完后用同步回代削减到30个场景。因为场景维度较高72维距离矩阵会受“维数灾难”的影响欧氏距离在高维空间区分度会下降。折中做法是先对72维场景做PCA降维到10-20维把累计方差贡献率保持到95%以上在高维空间做削减再把场景映射回原空间。我实测过这种先降维再削减的方案和直接在高维空间削减相比不仅在计算时间上快一个量级削减质量以目标函数偏差衡量也更好。7.3 削减结果的工程验证削减完成后我做三重验证矩检验计算原场景集与削减后场景集的均值差和标准差差确保均值偏差不超过原标准差的2%标准差偏差不超过5%分位数检验比较原场景集和削减后场景集在 5%、50%、95% 三个分位点上的出力水平差异经济性验证把原场景和削减场景分别代入一个简化的机组组合模型比较两者算出的总运行成本差距偏差控制在2%以内才算过关这套验证体系并不能保证万无一失但作为工程验收标准已经足够。真正到了学术研究级别需要更精细的Wasserstein距离评估和灵敏度分析但通用工程项目的投资决策对2%以内的目标函数误差是完全能接受的。最后分享一点个人体会。场景削减这个环节在整套风光模拟和优化体系里占据的位置看起来不大但它直接决定了后续所有分析的输入质量。我见过不少项目在预测模型、优化算法上投入大量精力却因为场景削减做得粗糙导致一整条分析链的结论偏掉。与其追求更花哨的不确定性模型不如多花点时间把场景削减这一步做扎实——把距离度量选对、把概率权重算清、把极端场景保住。这步地基稳固了上面盖什么楼都不慌。如果后续有精力我打算再写一篇把“风电光伏联合出力场景模拟”和“源荷互动下的场景削减”结合起来的内容涉及负荷相关性的部分会更复杂一些。先到这里大家在实操中遇到有意思的场景削减案例也欢迎一起交流。 基于正态分布的风光场景削减：从协方差矩阵到同步回代实现
1. 问题的起点为什么风光模拟需要“场景削减”这一刀做新能源消纳分析、电力系统随机优化或者储能容量配置的朋友大概率都撞过同一堵墙“源”端的不确定性怎么建模。风电看风速光伏看辐照度两者都有明显的随机性和间歇性直接把历史数据灌进优化模型求解器分分钟教你做人——不是收敛不了就是计算时间长得能去倒杯咖啡再回来。于是业内的通用做法是用随机模拟生成大量可能的风光出力场景用场景集去近似真实的不确定性分布再在这个场景集上跑优化。这个方法本身没错问题出在“大量”两个字上。蒙特卡洛模拟随手一跑就是几千甚至上万条场景。目标函数每多一个场景优化问题的变量规模和约束条数就线性上涨。场景数量到500的时候混合整数规划基本就卡到没法看到1000以上不是算不动是压根没法用。这就是场景削减技术存在的根本动机——在保证概率分布特征尽量不损失的前提下用少量代表性场景替代原始大场景集把计算复杂度降下来。而标题里的“基于正态分布”是另一个容易被忽略的前提。风光的物理过程随机性极强理论上风速更接近威布尔分布、光照更接近贝塔分布但落到“模拟”这个层面时工程师们大量使用正态分布来做近似——因为中心极限定理在手样本量上来以后很多东西都趋向正态也因为正态分布的协方差结构最容易描述变量之间的相关性。用正态分布做风光模拟本质上是在说我们观测到的风光出力波动可以分解为一个确定性趋势加一个正态随机扰动而这个扰动的大小和变量间的联动关系全部由协方差矩阵刻画。所以这篇文章想聊清楚三件事正态分布模拟风光出力的数学框架怎么搭、场景削减的经典算法怎么一步步落地、以及实操中那些论文里不会写的坑。我会按一条完整的链路往下走代码都是可以直接跑的环境是Python 3.10 NumPy SciPy不需要任何商业求解器。2. 用正态分布刻画风光出力从物理过程到协方差矩阵先别急着写代码得先把模拟框架的底层逻辑捋顺不然生成出来的场景库就是一堆没有物理意义的随机数。2.1 为什么正态分布在风光模拟里这么能打很多人第一次接触风光模拟时都会问风速明明是偏态分布光照也带着明显的上下限截断为什么没人用威布尔或者贝塔直接采样核心原因有三个。第一高维相关性建模的难度天差地别。风光模拟不是只模拟某一点某时刻的风速而是要模拟一个风电场里多台机组、一片区域里多个场站、一段时间内多个时刻的联合出力。威布尔和贝塔分布当然都能做单变量拟合但一旦要刻画“所有位置、所有时刻”之间的复杂相关性这两个分布都没有现成的多变量版本。多变量正态分布虽然也不是万能的但它有完整的数学框架——协方差矩阵一个矩阵把所有的线性和相关性都装进去了。第二中心极限定理在实数世界里确实在起作用。风电场的出力是几十上百台风机出力的叠加一个电站的出力又是多个气象物理过程的综合结果。叠加得多了分布自然会向正态靠拢。这个现象做风电功率预测的朋友应该深有体会——单机功率的预测误差分布可能是尖峰厚尾的但到了场站级别的聚合功率误差分布就已经非常接近正态了。第三工程实现的成本优势。多变量正态采样的实现只需要做一次Cholesky分解加一次线性变换计算量几乎可以忽略。相比之下copula或经验分布采样虽然在某些案例里更贴合原始数据但参数估计和场景生成的计算成本高得多在“先跑通再精调”的工程节奏下正态分布几乎是必然的第一选择。2.2 风速与光照模拟的基本方程我们直接从最常用的物理近似模型切入。对于风速很多工程模拟采用一阶自回归模型AR(1)或者基于平均风速加噪声的简化方式。单个时段的风速可以用这个形式表达v_t μ_v ρ · (v_{t-1} − μ_v) ε_t其中 μ_v 是该季节或时段的风速均值ρ 是时间自相关系数体现风速的连续性ε_t 是服从 N(0, σ_v²) 的随机扰动。对于光伏辐照度序列同样可以用类似的模型但要额外加入昼夜周期和云层遮挡的随机偏差。如果只做静态场景不考虑时序相关性可以直接跳过AR项用下面的多变量正态形式对N个位置的风速一次性采样v ~ N(μ, Σ)这里 v 是N维列向量μ 是每个位置的平均风速向量Σ 是风速在空间上的协方差矩阵——距离越近的风机场风速相关性越高这个相关性就体现在协方差矩阵的非对角元上。光资源部分写法完全相同r ~ N(μ_g, Σ_g)其中 r 是N个光伏电站的辐照度向量μ_g 是平均辐照度Σ_g 是空间相关协方差矩阵。有些精细一点的模拟还会把温度对光伏效率的影响也建模进去做法是给辐照度乘一个温度修正系数这里先按下不表。2.3 功率输出非线性环节的正态化处理风速和辐照度并不是我们要的最终场景——优化模型里真正关心的是功率出力。风速到功率的转换关系由风机的功率曲线决定典型的风机功率曲线是这样一个分段函数风速低于切入风速 v_in出力为 0典型切入风速约 3 m/s风速在 v_in 和额定风速 v_r 之间出力随风速立方增长约 12-14 m/s风速在 v_r 和切出风速 v_out 之间出力等于额定功率风速高于切出风速风机停机保护出力为 0也就是说即使风速服从正态分布经过功率曲线这个非线性环节后功率的分布形态已经不再是正态。这看起来是个矛盾但工程上恰恰有对应的处理思路直接对功率做正态近似而不是对风速做正态假设再把功率算出来。具体做法是这样的——先取得风电场的历史功率数据用核密度估计或者经验分布得到功率的分布特征然后直接用多变量正态分布去拟合功率序列的联合分布。这样的好处是绕开了功率曲线非线性转换的麻烦协方差矩阵里装的全是功率层面的相关关系采样出来的结果直接就是功率场景。实测下来风电场的多机聚合功率用正态近似比用风速转功率更接近真实统计特征尤其是尾部的相关性——因为功率曲线的截断效应已经把风速尾部的极端值打平了。2.4 协方差矩阵的构造空间相关性与时间相关性的统一协方差矩阵是整套模拟框架里最核心的构件。构造它有两种思路。第一种基于历史数据的经验估计。拿两年以上的风光历史出力数据每个位置一条序列先减去各自均值去趋势然后直接计算两两之间的协方差。这种方法最真实因为它把所有的相关结构都隐式地包含在数据里了——包括不同位置之间的空间相关、时段之间的耦合、甚至风光互补的负相关性。缺点是需要的数据量足够大数据质量差或者有缺失时段时协方差矩阵容易出现不正定的问题后面会专门讲这个坑。第二种基于距离衰减函数的参数化构造。没有足够的历史数据时用地理距离来推算相关系数是一个很常见且好用的简化。空间相关系数和距离之间的关系可以近似表达为ρ_ij exp(−(d_ij / L)^2)其中 d_ij 是两个电站或场站间的距离L 是相关长度参数correlation length通常取一两百公里。这个公式的物理含义一目了然——距离越近相关越强距离超过相关长度后相关迅速衰减到接近零。实际操作中风力资源和光资源的相关长度取值差异很大风的相干尺度明显大于光的相干尺度光资源在几十公里外相关性就已经掉得很快了。把两种方式混合使用也常见——区域间的整体相关性用经验估计区域内部由于缺少测点而用参数化方程补全。构造完成后的 Σ 一定是对称矩阵但要额外检查它是否正定所有特征值都大于0不正定的矩阵在采样时会直接报错或者给出离谱的样本这个检查后面必须做。2.5 场景生成的完整流程附可直接运行的代码把上面的逻辑整合成一套标准的场景生成Pipeline大概是这个顺序定义位置信息风电、光伏场站的经纬度或距离矩阵构造协方差矩阵经验估计或距离衰减函数可混合使用指定均值场景基准季节/时段/气象场景的平均出力水平用多变量正态分布采样生成目标数量的随机出力场景对采样的场景做物理约束修正功率上下限截断等下面这段代码可以直接跑通它生成了一个含20个风电站点的场景集每个场景代表了选定时刻的功率预测水平占装机容量的比例import numpy as np # Step 1: 站点配置 n_sites 20 positions np.linspace(0, 300, n_sites) # 假设站点沿一条线分布单位km # Step 2: 用距离衰减函数构造空间协方差矩阵 L_corr 120.0 # 相关长度单位km dist_matrix np.abs(positions[:, None] - positions[None, :]) corr_matrix np.exp(-(dist_matrix / L_corr) ** 2) # Step 3: 定义风速均值与标准差 mu 8.0 2.0 * np.sin(positions / 50.0) # 模拟地形的风速均值变化单位m/s sigma 1.8 # 单点风速标准差单位m/s cov_matrix sigma * sigma * corr_matrix # 协方差矩阵 相关系数矩阵 * 方差缩放 # 检查正定性 eigvals np.linalg.eigvals(cov_matrix) print(最小特征值, eigvals.min()) # Step 4: 蒙特卡洛采样 n_scenes 2000 rng np.random.default_rng(42) wind_speed_scenes rng.multivariate_normal(mu, cov_matrix, sizen_scenes) # Step 5: 功率转换简化版功率曲线忽略切入切出细节 def power_curve(v, v_in3.0, v_r13.0, p_rated1.0): p np.zeros_like(v) mask_ramp (v v_in) (v v_r) mask_rated v v_r p[mask_ramp] p_rated * ((v[mask_ramp] - v_in) / (v_r - v_in)) ** 3 p[mask_rated] p_rated p[v 25.0] 0.0 # 切出风速 return p power_scenes np.array([ np.array([power_curve(v) for v in scene]) for scene in wind_speed_scenes ]) print(场景维度, power_scenes.shape) print(平均功率水平, power_scenes.mean())这段代码生成的power_scenes就是最原始的“大场景集”每个场景代表所有站点在同一时刻的出力水平向量。3. 场景削减的数学本质把聚类问题翻译成概率问题拿到几千个场景之后削减这一步的核心任务是什么直白地说从原场景集中挑出或合并出一个小得多的代表性场景子集让这个子集的概率分布和原集合尽可能接近。这里有个容易混淆的点值得单独拧出来说清楚。场景削减不是“筛出几个典型日”虽然形态上有点像但方法论完全不同。典型日选择是在历史数据里直接挑几条有代表性的曲线保留的是原始数据的形态场景削减则是在概率空间里做最优化——每个场景会被赋予一个发生概率削减后的场景集带着一组新的概率权重让这些带权的场景在统计意义上逼近原分布。为什么要在这个维度上较真因为后续做的优化模型里目标函数的系数不只是场景本身的物理量还乘了每个场景的发生概率。削减前后如果概率权重的分布偏差过大最优解就会失真——可能原来排序第一的调度方案在削减后就变成第三了这就违背了场景削减的初衷。3.1 一个直观的场景削减例子从27个场景到5个场景先看一个能手工验证的小例子帮助建立直觉。假设原始场景集里有27个二维点它们分成三个明显的聚类每个聚类里9个点点的坐标从 (0,0) 附近到 (10,10) 附近分布。削减到5个场景后理想结果应该是每个聚类里保留1-2个代表点每个点带一个概率权重权重数值约等于该聚类中所有点代表的概率份额之和。这样得到的5个场景放在优化模型里的时候计算量缩减到原来的五分之一左右而由于每个聚类里点的内部差异不大合并后保留的代表点又能很好反映局部特征优化结果的偏差通常被控制在一个很小的范围内。3.2 场景削减要保留的三种分布特征在工程应用场景里有三种分布特征是最关键的削减后的场景集必须尽量保留它们。第一种是均值与方差。这是最基础的矩信息均值直接决定优化的目标函数基准值方差则反映不确定性的大小。削减后场景集的加权均值和加权协方差必须与原场景集接近否则优化结果会系统性偏移。第二种是概率分布的形态。虽然我们从上到下都用正态分布建模但削减后的分布不一定能保持严格的联合正态性因为保留的场景是离散的点集它们的分布会是“一堆discrete mass”。在衡量削减质量时通常用 Wasserstein 距离推土机距离或 KL 散度来定量评估两个分布的接近程度。Wasserstein 距离实际上衡量的是“把一个分布推成另一个分布需要的最小代价”这个指标在场景削减的学术文献里是标配。第三种是极值场景。这是实操中最容易翻车的地方。削减算法在追求整体分布逼近时往往倾向于保留“那些更像中心点的场景”而这些场景恰恰会把高风电出力、低光伏出力这种极端情况抹掉。但电力系统优化里最需要关注的往往就是极端场景——系统要保持安全性恰恰要针对这些边界条件校核。所以有经验的工程师会在削减后手动检查极值分位数必要时强制把几个极端场景加回去这是非常典型的工程补丁操作。4. 同步回代削减法O(n²) 手写实现与每一步的解释场景削减的算法有很多但从工程实用角度看同步回代削减法Simultaneous Backward Reduction是使用最广泛、最容易落地的方案。它是Birge和Louveaux在《Introduction to Stochastic Programming》里详述的经典方法也是MATLAB和各类随机优化工具包里最常见的实现之一。它的优点是思路直观、实现简单、对场景维度的适应性好缺点是计算复杂度是场景数的平方场景上万后会很慢。但在工程项目里原始场景集通常在一两千以内跑起来完全无压力。4.1 同步回代的核心逻辑每次删一个删完再重算距离这个方法的名字里“回代”两个字听着玄乎实际逻辑特别简洁定义场景之间的距离度量通常是欧氏距离或加权欧氏距离权重可以按各维度的重要性设计算所有场景两两之间的距离得到一个距离矩阵对每个场景 i找到离它最近的场景 j记下这个最近距离在所有场景中找出“删除代价”最小的那个场景——删除代价定义为它和最近邻居的距离乘以它自己的概率把这个场景删掉同时把它携带的概率加到它的最近邻居上重复第2-5步直到场景数降到目标值很多人第一次看这个流程时容易懵的点在“为什么要加概率”。理由很朴素场景 i 被删除后它的存在价值并没有消失而是应该由最近的邻居场景 j 来承担因此原本属于 i 的发生概率也一并转移给 j。这样一步步合并下去最后剩下的每个场景的概率权重都是它下面“吸收”掉的一堆兄弟场景的累积概率。4.2 从零手写同步回代削减Python实现这里给出完整可运行的实现。为了代码可读性我把每一步都拆成独立小节并用中文注释解释每段逻辑的目的。import numpy as np from scipy.spatial.distance import cdist def sync_backward_reduction(scenarios, probsNone, target_count10): 同步回代削减法 参数: scenarios : ndarray, shape (n_scenes, dim) 每个场景是一行特征向量 probs : ndarray, shape (n_scenes,) 每个场景的原始概率 target_count : int 削减后保留的场景数 返回: reduced_scenarios : ndarray reduced_probs : ndarray removed_indices : list n scenarios.shape[0] if probs is None: probs np.ones(n) / n # 给每个场景保存携带概率的副本后续会被修改 p probs.copy() # 保存场景索引到原数据集的映射 active_indices list(range(n)) # 预先计算一次距离矩阵后续删除场景时不需要重算全量距离 dist_matrix cdist(scenarios, scenarios, metriceuclidean) np.fill_diagonal(dist_matrix, np.inf) # 禁止自己匹配自己 removed_indices [] while len(active_indices) target_count: m len(active_indices) # 计算每个场景与其最近场景的距离 # 由于active_indices只是原索引的映射我们需要在当前活跃子集内找最近邻 sub_dist dist_matrix[np.ix_(active_indices, active_indices)] min_dist sub_dist.min(axis1) # 删除代价 概率 * 最近距离 cost p[active_indices] * min_dist # 找到要删除的场景代价最小 del_local_idx np.argmin(cost) del_orig_idx active_indices[del_local_idx] # 找到它的最近邻居在删除前先确定因为删除后邻居关系会变化 closest_local_idx np.argmin(sub_dist[del_local_idx]) neighbor_orig_idx active_indices[closest_local_idx] # 删除场景的概率合并到邻居身上 p[neighbor_orig_idx] p[del_orig_idx] p[del_orig_idx] 0.0 # 从活跃列表中移除 removed_indices.append(del_orig_idx) active_indices.pop(del_local_idx) # 输出最终保留的场景 reduced_scenarios scenarios[active_indices] reduced_probs p[active_indices] # 归一化概率权重因浮点误差可能略有偏差 reduced_probs reduced_probs / reduced_probs.sum() return reduced_scenarios, reduced_probs, removed_indices逐段解释一下实现里的关键设计。dist_matrix使用scipy.spatial.distance.cdist预先计算对 2000 个场景、每个场景 20 维的情况一次计算只需几十毫秒。后续每次迭代只需要在活跃子集内取距离子矩阵避免了重复计算全量距离的开销。真正需要逐轮重复计算的是“每个场景到最近邻的最小距离”和“删除代价”这部分复杂度是 O(n²)因为要在活跃子集内找两两距离的最小值。删除时选择np.argmin(cost)而不是np.argmin(min_dist)两者唯一的差别就是乘以了场景概率 p。为什么要乘概率因为一个场景虽然离得很远孤立点如果它本身出现概率极低删除它带来的分布损失也很小相反一个高频场景即使离邻居很近删了它也会造成较大的特征损失。概率权重本质上是给每个场景的“影响力”打了折让削减决策更贴近真实概率意义。一个容易出错的点邻居的查找必须在删除动作之前完成。因为删掉一个点后距离矩阵的行列会发生变化如果你先删了再找邻居找到的可能是“删完了以后的新最近邻”而不是“删除时刻的最近邻”概率合并的对象就可能不是最优的了。这个顺序问题我见过不少实现栽进去的坑写的时候要时刻盯着。4.3 削减效果评价用“成本函数偏差”说话削减做得好不好光看场景形态不够要用定量的指标说话。这里分享一个我项目里一直沿用的评估方案——以随机优化问题下的期望成本偏差作为最终衡量标准。假设原场景集对应的期望成本是J_original Σ_{i} prob_i * cost(scenario_i)削减后成本是J_reduced Σ_{j} prob_j * cost(scenario_j)两者之差就是削减引入的成本近似误差。这里的 cost 函数可以简单点用线性成本比如“总出力加权和”或者“缺口平方罚”也可以用真实的机组组合运算算出来的总运行成本。为了演示写一个简单的评估脚本直接用生成的风力场景集跑一遍削减流程# 沿用前面生成的 power_scenes scenes power_scenes # (2000, 20) n_orig scenes.shape[0] uniform_probs np.ones(n_orig) / n_orig reduced_scenes, reduced_probs, removed sync_backward_reduction( scenes, uniform_probs, target_count50) print(原始场景, n_orig, 条) print(削减后场景, reduced_scenes.shape[0], 条) print(削减后概率之和, reduced_probs.sum()) # 成本函数以全系统总出力的负偏差平方为评价模拟功率不足的惩罚 def cost_of_scene(scene, target_level400.0): total_power scene.sum() deficit max(0.0, target_level - total_power) return deficit ** 2 # 计算期望成本 cost_orig np.sum(uniform_probs * np.array([cost_of_scene(s) for s in scenes])) cost_reduced np.sum(reduced_probs * np.array([cost_of_scene(s) for s in reduced_scenes])) print(原始期望成本, cost_orig) print(削减后期望成本, cost_reduced) print(成本偏差比例, abs(cost_reduced - cost_orig) / cost_orig)在我自己的数据集上用上述方案把2000个场景削减到50个成本偏差通常在0.5%-2%之间。这个精度对绝大多数工程决策场景都是完全可以接受的。换句话说用2.5%的场景数量换来不到2%的目标函数精度损失性价比非常高。4.4 削减到多少最合适先看求解器脸色再看精度曲线场景削减到什么数量才算“够”是很多新手的第一个疑问。给一个放之四海而皆准的数字是不负责任的但可以给出一个靠谱的工程流程先扔两个边界值5个场景和200个场景分别跑一遍优化模型记录求解时间如果5个场景下秒解200个场景下也能在可接受时间内解出那直接选200个附近精度最好如果200个场景下已经开始卡了就用15-50个场景区间内的不同数量各跑一遍画出“场景数-目标函数值”的收敛曲线选曲线拐点的场景数——超过这个数之后再加场景对目标函数值的改善已经非常微小这个流程是经验性的但非常实用。很多研究者喜欢掰扯各种指标来论证“最优削减数量”实际工程里根本不需要那么精致只要目标函数值和决策方案在削减前后保持在稳定区间内场景数越少越好。5. 除了同步回代还有哪些场景削减方案值得一试同步回代是经典但不是一个银弹。工程实践中我试过另外几条路线各有利弊可以按需选择。5.1 K-means 与 K-medoids切得更快的平替K-means的思路极度直接把场景当作空间里的点聚类成K簇取每簇的质心作为代表场景簇内场景的概率之和就是代表场景的概率。实现只需要调scikit-learn一行代码from sklearn.cluster import KMeans k 20 km KMeans(n_clustersk, random_state0).fit(scenes) labels km.labels_ centroids km.cluster_centers_ cluster_probs np.array([np.mean(labels i) for i in range(k)])K-means的最大优势是快——对一万个场景聚类到20簇也不用几秒。但有一个关键问题质心是簇内所有点的平均会“平滑”掉一些极端特征。在确定性优化里这无所谓但在不确定性优化里你恰恰希望代表场景保持可能出现的最不利情形的一些影子。而K-medoids对这个问题处理得更好它选的是“簇内最中心的那个原始场景”而不是计算出来的平均点from sklearn_extra.cluster import KMedoids km KMedoids(n_clustersk, random_state0).fit(scenes) medoids km.cluster_centers_为什么K-medoids在场景削减场景里更被推荐因为它选出的代表场景一定来自原始场景集意味着“这个场景是系统真实可能出现的出力形态”而不是一个虚构的平均态。这对后续做安全校核、调度员审阅结果、或者给非技术背景的决策者解释都重要得多。缺点也有——K-medoids的收敛速度比K-means慢一个数量级而且在样本量很大的时候更容易收敛到局部最优。5.2 层次聚类削减过程全透明层次聚类Hierarchical Clustering和同步回代的逻辑很相似——从下往上逐步合并且合并距离最小的簇。但实现方式上层次聚类用的是自底向上的“凝聚式”聚类每一步距离计算有多种准则单链接最近距离、全链接最远距离、平均链接平均距离、Ward方差增量。实际操作里Ward准则和平均链接准则在场景削减里表现最好单链接容易链式效应导致莫名其妙的狭长簇。好处是聚类树保留了完整的合并历史可以随时在任何层级停下来查看精简后的场景集坏处是时间复杂度较高一千个场景以上就需要等一会儿。层次聚类在“需要向别人解释整个过程”的场景下特别有价值——你可以画出树状图直观展示“这些场景是从哪些子群合并而来的”这在项目汇报时是很好的辅助材料。5.3 快速前向选择反向削减的镜像版本和同步回代每次删一个不同快速前向选择是从零开始每次从原始场景集中挑一个“加入当前集合后让分布误差下降最多”的场景加入代表集。听上去也合理但实操中它的表现通常不如同步回代——因为后向削减的每一步都基于当前完整集合做全局决策而前向选择每一步都只看到当前已挑选的子集全局视野会差一截。不过前向选择有一个同步回代没有的优点它天然适合“预先规定好要保留N个场景”的场景不用像同步回代那样逐步删到目标值后还要担心删过头。在小规模问题上各跑一遍逐项对比你会发现两者的结果差异通常不大但性能上前向选择略快。6. 实操中的五个大坑与解决方案这一节全是我在项目里实际踩过、并且花过很长时间才爬出来的坑。论文里通常只给你漂亮的公式和完美的结果图这些坑不写在纸面上。6.1 坑一协方差矩阵不正定采样直接炸裂用经验协方差矩阵做采样特别容易触发LinAlgError: Matrix is not positive definite。原因很简单如果你的站点数量比历史样本时段数还多比如20个站点、只有15天的逐小时数据那协方差矩阵就是15×15的自然不正定或者两个站点的数据序列几乎完全线性相关比如距离只有5公里的两个风机相关系数直接到0.9999矩阵就会退化。解决思路有三个加正则项求解协方差矩阵的特征值分解后把所有小于某个阈值的特征值直接设为一个正的小值比如1e-6再重构矩阵。这种做法最省事损失也最小。收缩估计用sklearn.covariance.ShrunkCovariance对经验协方差做收缩处理把对角元稍微放大、非对角元稍微缩小能同时提升矩阵的正定性和数值稳定性。降维预处理先用PCA对原始数据进行降维在降维后的空间里构造协方差再映射回去。此法在极端高维场景下有用但在风电场站这种10-50维的规模下用不着这么重的手段。建议从小问题开始就直接用“加正则项 Cholesky 分解后做下三角矩阵乘法采样”这条路线工程复杂度最低。def safe_multivariate_normal(mean, cov, size, reg1e-6): eigvals, eigvecs np.linalg.eigh(cov) eigvals np.maximum(eigvals, reg) cov_fixed (eigvecs * eigvals) eigvecs.T return rng.multivariate_normal(mean, cov_fixed, sizesize)这里np.linalg.eigh专门用来做对称矩阵的特征分解比eig数值上更稳定执行效率也更高。6.2 坑二削减后保留了一堆长相差不多的场景这是同步回代的“副作用”之一。算法倾向于保留高频聚类附近的代表点如果原始场景集中某一片区域特别密集削减后可能发现留下来的5个场景里4个都挤在这一个密集区域其他区域反而一个代表都没有。这种“选择性偏差”会让后续优化结果产生偏向性的偏移。解决方式是在削减前后都做一步“场景多样性检查”计算场景集两两之间的最小距离如果最小值过小比如小于某个物理阈值如风电场总出力的2%就要手动介入调整。实际项目里最常用的做法是先跑同步回代拿到基础削减结果再对保留的代表场景和几个预选的边界场景最大出力、最小出力、最大爬坡等做一次手工校正把边界场景的权重适当抬高。6.3 坑三概率加在邻居头上的累积误差场景削减中概率转移带来的问题是如果削减比例过大比如从5000个场景削减到20个每个保留场景的累积概率可能非常大而它的代表性却无法完全覆盖那些被合并不进来的细节。当概率权重分配不平衡时优化结果对“高概率场景”过度敏感反而比原始的一组低概率场景更脆弱。应对办法是把削减比例控制在合理区间内经验上限是“削减后的场景数不少于原始场景数的2%~3%”条件允许时建议“不少于5%”。如果因为计算能力限制必须削减到极少数那就先削减到50个再另想办法降复杂度比如用聚合机组建模、线性化约束不要在一棵树上吊死。6.4 坑四不同类型资源风光的削减权重不平衡现实中一个系统往往同时包含风电场和光伏电站而风光在时间和空间上的出力模式差异很大。如果直接用欧氏距离做削减风电场部分的数值范围0~额定装机可能比光伏部分大得多尤其在夜间距离度量会被风资源维度主导光伏出力的差异性在削减过程中被严重忽略。解决方案是给距离度量加上加权系数。权重的选取可以按资源的波动性程度来定——越需要精确保留的维度权重越大。我在实际项目里用过的一种做法是先把所有维度归一化除以各自装机容量再给风资源维度加权系数1.0光资源维度加权系数1.5。原因是光伏出力的时序自相关更强、突变性更明显微小的相位偏差在优化目标里的成本惩罚会放大所以同类距离下光资源的“信息量单价”更高。需要精确处理的时候可以直接让场景削减的目标函数和优化模型的目标函数挂钩以最终运行成本最小化的拉格朗日乘子作为各维度的权重系数——这个方案做起来复杂但效果确实最扎实。6.5 坑五场景削减结果缺乏“可解释性”场景削减产出的结果往往是一个特征向量加上一个浮点数概率。给领导汇报时需要能够清楚地解释“这个场景代表了什么、为什么它占了25%的概率”。如果削减方法选得不好保留的场景从物理角度看几乎不可能同时发生——比如风速很高的同时光照也极度强烈而真实天气过程里这两者往往负相关——这样的场景虽然数学上合理但物理上缺乏说服力。解决办法是在削减前对场景的“物理解释性”做一次筛查检查连乘概率是否过低的组合比如联合概率小于1e-4把它们剔除或重新合并。这个步骤虽然看起来有点“反数学”但在实际工程交付中帮你省掉的解释成本远大于数学上的微小损失。7. 一次完整的实战从原始数据到可用的削减场景库理论和坑都讲完了最后走一遍完整的实战流程把前面所有部件串起来。假设我手上有一个区域系统里面有2个风电场和1个光伏电站要做未来一年的典型场景削减。7.1 数据准备与预分析数据来源是各场站的SCADA系统功率历史数据采样频率15分钟跨度为最近两年。第一步做数据清洗剔除停机检修段、限电段、通信异常段必要时用前后时段平均值插补。这一步不能省因为协方差矩阵对异常值极其敏感。清洗完成后将所有功率数据归一化到0-1之间除以各自装机容量这样不同场站之间可以直接比较。接下来计算三个场站间的相关系数矩阵——预计两个风电场之间的相关系数会比较高0.7-0.9而风电场与光伏电站之间的相关系数会明显偏低甚至为负-0.1到0.3这恰好反映了风光资源的天然互补特征。7.2 场景生成与压缩以小时为粒度每个场景的时间长度为24小时那么每个场景的维度为24 × 3 72维。按照第2节的框架用多变量正态采样生成2000个初始场景。这里有个容易忽略的问题——正态分布在采样过程中可能生成负功率或超过装机容量的数值要在预分析阶段就把负值截断到0、超过1的压回1。生成完后用同步回代削减到30个场景。因为场景维度较高72维距离矩阵会受“维数灾难”的影响欧氏距离在高维空间区分度会下降。折中做法是先对72维场景做PCA降维到10-20维把累计方差贡献率保持到95%以上在高维空间做削减再把场景映射回原空间。我实测过这种先降维再削减的方案和直接在高维空间削减相比不仅在计算时间上快一个量级削减质量以目标函数偏差衡量也更好。7.3 削减结果的工程验证削减完成后我做三重验证矩检验计算原场景集与削减后场景集的均值差和标准差差确保均值偏差不超过原标准差的2%标准差偏差不超过5%分位数检验比较原场景集和削减后场景集在 5%、50%、95% 三个分位点上的出力水平差异经济性验证把原场景和削减场景分别代入一个简化的机组组合模型比较两者算出的总运行成本差距偏差控制在2%以内才算过关这套验证体系并不能保证万无一失但作为工程验收标准已经足够。真正到了学术研究级别需要更精细的Wasserstein距离评估和灵敏度分析但通用工程项目的投资决策对2%以内的目标函数误差是完全能接受的。最后分享一点个人体会。场景削减这个环节在整套风光模拟和优化体系里占据的位置看起来不大但它直接决定了后续所有分析的输入质量。我见过不少项目在预测模型、优化算法上投入大量精力却因为场景削减做得粗糙导致一整条分析链的结论偏掉。与其追求更花哨的不确定性模型不如多花点时间把场景削减这一步做扎实——把距离度量选对、把概率权重算清、把极端场景保住。这步地基稳固了上面盖什么楼都不慌。如果后续有精力我打算再写一篇把“风电光伏联合出力场景模拟”和“源荷互动下的场景削减”结合起来的内容涉及负荷相关性的部分会更复杂一些。先到这里大家在实操中遇到有意思的场景削减案例也欢迎一起交流。