★ FEATURED ARTICLE 在新能源项目里摸爬滚打久了你会发现一个特别有意思的现象无论做风电功率预测、光伏容量配置还是微电网调度策略大家最终都在跟一个问题较劲——风光出力到底怎么描述才靠谱。你说用平均出力调度方案直接失真极端天气一来就抓瞎。你说用典型日曲线一年8760小时那真是成千上万种可能性典型日根本代表不了全貌。所以业内普遍的做法就是用蒙特卡洛模拟配合场景削减把“说不清”的风光出力转成一堆看得见、摸得着、能参与计算的代表性场景。这篇文章我就结合自己做过的模拟项目把这个套路从头到尾拆一遍包括场景是怎么生成的、削减算法是干什么吃的、参数怎么设置以及那些最容易踩的坑。1. 内容整体设计与思路拆解1.1 从“确定性思维”到“随机场景思维”的转变初做风光出力模拟的人最容易犯的毛病是拿着历史平均风速、平均光照强度直接推功率。这样算出来的结果是一条光滑的、毫无脾气的出力曲线。真拿去算经济调度问题就来了实际风电场出力波动极大一条平均曲线把所有风险都给抹平了算出来的方案激进得吓人。蒙特卡洛方法的思路完全不同。它不追求用一条曲线代表所有情况而是承认风光出力的随机本质然后用概率分布描述这种随机性通过大量随机采样去逼近未来可能出现的各种运行状态。换句话说一次采样对应一个“如果”成千上万次采样你就得到了成千上万个可能性。这些可能性就是原始场景集。但这里有个现实的矛盾采样次数少了生成不了足够的可能性采样次数多了比如抽了5000个场景每一个都扔进优化模型里算一遍计算量直接爆炸。所以另一个关键技术就出场了场景削减。它的目标是在保留原始场景集统计特征的前提下用尽量少的代表性场景代替大规模场景集让计算量降下来还不损失太多精度。1.2 解决方案的整体架构我用的整体方案可以拆成四层第一层输入数据层包括风速、光照强度的历史统计参数通常是威布尔分布或正态分布的拟合结果。第二层场景生成层利用蒙特卡洛采样生成N个24小时的出力场景含时序相关性处理。第三层场景削减层用聚类或快速前向选择算法把N个场景削减到用户指定的K个。第四层应用层削减后的典型场景和对应概率直接作为优化模型如机组组合、储能容量配置的输入。整个流程一句话概括采样生成大量可能 → 削减保留少量精华 → 带着概率进优化。这个思路不仅适用于风光出力只要你能给出随机变量的概率分布这个框架几乎都能套。2. 核心细节解析场景生成的关键步骤与参数设计2.1 风速与光照的概率建模场景生成的第一步是给随机变量选一个合适的概率分布。风电出力最常用的风速建模工具是两参数威布尔分布概率密度函数长这样[ f(v) \frac{k}{c}\left(\frac{v}{c}\right)^{k-1} e^{-(v/c)^k} ]其中 ( k ) 是形状参数控制风速分布的“胖瘦”( c ) 是尺度参数控制平均风速的高低。我有个项目用的实测数据显示某风电场全年的拟合结果是 ( k2.1 )( c7.8 )意味着平均风速大概在6.9米/秒附近风速波动比较集中极端大风占比不高。光伏出力则完全换一套逻辑它的随机性主要来自光照强度业内常用Beta分布建模。光照强度归一化处理后Beta分布的概率密度为[ f(x) \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)} ]其中的 ( \alpha ) 和 ( \beta ) 是形状参数从历史辐照度数据的均值和方差反推得到。实际操作中我常用矩估计法先算历史光照数据的均值 ( \mu ) 和方差 ( \sigma^2 )然后套这两个式子推算[ \alpha \mu \left(\frac{\mu(1-\mu)}{\sigma^2} - 1\right) ] [ \beta (1-\mu) \left(\frac{\mu(1-\mu)}{\sigma^2} - 1\right) ]这里要特别注意一个细节光照强度有很强的日周期性——白天有光、夜里归零。如果只做单纯的分布采样很容易采出“半夜光照充足”这种荒谬场景。所以必须按小时划分时段对每一小时单独建模。比如上午10点这个时段光照均值可能已经升到0.55方差也大深夜2点则基本趋近于0。把24小时每个时段各拟合一组Beta分布参数生成的时候逐时段采样时序特征就保留下来了。2.2 时序相关性处理不能让场景“各说各话”蒙特卡洛采样的一个天然缺陷是每次采样互相独立。如果你逐小时独立抽样生成出来的风速曲线会在时间轴上疯狂跳变——前一小时风速8米/秒后一小时直接掉到2米/秒现实中基本不会发生。这种“伪随机”场景拿去做调度计算结果必然失真。解决这个问题的办法很多我常用的是AR(1)一阶自回归模型加上逆变换采样。思路是这样的先利用历史数据建立时间相关性模型生成具有时序相关性的随机数再通过逆变换映射到目标分布上。简单来说我先用历史风速序列估计自相关系数 ( \phi )通常一小时尺度的风速自相关系数在0.8到0.95之间。然后生成一个服从标准正态分布的白噪声序列 ( \epsilon )按下面的方式构造一个标准正态的时序序列[ z_t \phi z_{t-1} \epsilon_t ]最后利用逆变换采样的核心结论任意随机变量 ( X ) 的累积分布函数 ( F_X(x) ) 是区间(0,1)上的均匀分布所以将 ( z_t ) 经过 ( \Phi(z_t) )标准正态的累积分布转成均匀分布数再带入威布尔分布的逆函数 ( F_W^{-1}(\cdot) )就能得到带相关性的威布尔风速序列[ v_t F_W^{-1}(\Phi(z_t)) ]这套组合拳下来生成的场景既能保持历史风速的分布特征又在时间上像真实的天气过程一样平滑连续。光伏光照序列的处理方式同理只是把目标分布换成Beta分布。这是我试过最稳健的方案没有之一。2.3 功率转换模型与出力曲线生成风速和光照序列只是中间产物最终要的是出力。风电功率转换常采用分段函数模型这是行业最普遍的做法具体形式如下当 ( v v_{in} ) 或 ( v \ge v_{out} ) 时出力为0。当 ( v_{in} \le v v_r ) 时出力按线性或三次方关系上升。当 ( v_r \le v v_{out} ) 时出力恒定等于额定功率。[ P_w \begin{cases} 0, v v_{in} \text{ 或 } v \ge v_{out} \ \dfrac{P_r}{v_r - v_{in}}(v - v_{in}), v_{in} \le v v_r \ P_r, v_r \le v v_{out} \end{cases} ]我用的风电机组切入风速3米/秒、额定风速12米/秒、切出风速25米/秒额定功率1.5兆瓦。光伏出力则简化为功率与光照强度的线性关系加上通过效率[ P_s \eta \cdot S \cdot A \cdot G_t ]其中 ( \eta ) 是综合转换效率( S ) 是光伏板面积( A ) 是光照强度归一化值( G_t ) 是参考辐照度。这套模型简单、参数少、够用做场景分析完全没问题。生成场景的时候直接对24小时逐时段风速/光照采样、转换成功率就得到了一个24维的出力向量 ( [P_1, P_2, ..., P_{24}] )。多采样几轮场景集就出来了。2.4 生成规模的确定场景生成规模 ( N ) 怎么定我一直遵循一个原则样本量要达到“统计特征稳定”的水平。我做过一个简单实验分别生成500、1000、2000、5000个场景统计均值、方差和各分位数。结果发现2000个场景以上时均值和方差的抖动已经很小到5000个场景时基本完全稳定。所以我个人建议常规项目取2000到3000个初始场景就够追求更极端情况的预测可以取到5000个。场景数量再多一方面计算时间明显增加另一方面削减算法也会更吃力边际收益其实很低。从工程角度来说N取值3000是个性价比非常高的选择。3. 场景削减算法从“海量可能”到“少数代表”3.1 为什么必须削减计算代价与信息保留的博弈你可能会想问蒙特卡洛生成了5000个场景每个场景都是独立可能的未来凭什么把它们硬砍到10个或者20个信息不就丢了吗这个问题的关键在于优化模型的计算代价。以机组组合问题为例目标函数涉及到每个场景下的运行成本和启停状态约束模型规模是原本的 ( N ) 倍。当 ( N5000 )哪怕是一个中等规模的微网系统求解时间也从分钟级直接飙升到小时甚至天数。更麻烦的是大量概率相近、形状相似的场景会让矩阵病态求解器更容易陷入数值问题。所以场景削减的真正目标不是“丢信息”而是“去冗余”把分布相近的场景合并成一个代表场景并赋予一个能反映原场景集权重的概率让削减后的场景集与原始集的统计特征尽可能一致。常见的衡量指标是K-L散度或Wasserstein距离前者大家了解得多一点后者在实际计算中效果往往更好。3.2 方案一K-means聚类削减K-means做场景削减的思路很直观把N个场景看作N个24维空间中的点聚类成K类每一类的质心就是代表场景类中场景数量占比就是该场景的概率。具体步骤是这样的将N个场景向量组成矩阵 ( X \in R^{N \times 24} )。随机初始化K个聚类中心。计算每个场景到各中心的欧氏距离或马氏距离分配到距离最近的中心。更新中心每个类内所有场景向量的均值成为新中心。重复3和4直到中心变化小于阈值或达到指定迭代次数。统计每个类中的场景数量 ( n_k )概率 ( p_k n_k / N )。K值的选取我试过用肘部法则看SSE随K下降的趋势拐点也试过直接用求解器算最终目标值来对比。如果从最终优化结果角度看K在10到20时削减前后的目标函数误差通常能控制在3%以内。所以一般结合项目需求取K10到20这之间的差异对结果影响不大。做K-means有几个细节容易坑人。第一一定要先做归一化。24维向量里如果直接做欧氏距离某些时段的数值天然偏大比如白天光伏出力接近1深夜接近0会无形中给这些时段加权重导致聚类结果偏向高值时段。把每个维度都归一化到0到1之后各时段的权重才是均衡平等的。第二距离度量不要死磕欧氏距离。如果你关心曲线的“形状”相似性而非“幅值”相似性动态时间规整距离会更合适但计算量会上好几个台阶小规模场景还行规模一大就有点吃力。第三K-means对初始中心敏感我习惯跑10次随机初始化保留SSE最小的那一次避免掉进局部最优的坑里。3.3 方案二基于Wasserstein距离的快速前向削减如果说K-means是“先聚类再取中心”那快速前向选择削减就是“从头到尾按相似度逐个挑代表”。它的算法原理是这样的假设你有初始场景集合 ( \xi {\xi_1, \xi_2, ..., \xi_N} )各自有初始概率 ( p_i )。目标是从中选出一个子集 ( S )让子集到原集合的Wasserstein距离最小。贪心策略的实现流程是初始时所有场景都在待选集合 ( J ) 中。找到一个场景索引 ( j )使得把 ( j ) 从待选集合移入保留子集后整体距离增量最小。将被 ( j ) “替代”的场景的概率累加到 ( j ) 上。重复步骤2和3直到保留场景数达到K。这里的“替代”关系通俗讲就是对于原始集合中每一个不在保留子集里的场景 ( \xi_i )找一个保留子集里距离它最近的场景 ( \xi_j ) 来代表它并把 ( \xi_i ) 的概率转移到 ( \xi_j ) 上。这样每次削减都能保证被删掉的场景概率被合理地“托付”给了最像它的那个代表场景整体的统计特征损失最小。相比K-means基于Wasserstein距离的快速前向选择的优势在于它能更好地保持原场景集的边界特征和尾部特征。比如风速场景里的极端大风场景、光伏场景里的连续阴雨天场景这些在K-means的质心里很可能会被稀释掉但在快速前向选择里它们更有可能被保留为独立的代表场景。所以在做极端事件分析或者可靠性评估的项目里我强烈建议用这种方法。3.4 两种削减算法的对比与选型建议对比维度K-means聚类快速前向选择实现难度低各类库直接调用中需一定算法实现计算耗时低收敛快高K个代表需要循环K次边界场景保留较弱质心会稀释极端较强代表场景更贴近原始场景概率分配方式类内占比最近邻概率转移适用场景常规规划问题、场景数大极端事件分析、可靠性评估常用场景数5000以上很轻松建议3000以内根据我的经验常规的容量配置、经济调度问题K-means完全够用而且效率极高涉及可靠性校验、备用容量决策的最好选择快速前向选择边界信息更重要。3.5 削减效果评估场景集质量怎么看削减完之后不能拍拍脑袋说“感觉差不多”你得拿出实实在在的指标来证明削减后的场景集还能代表原始场景集。我最常检查三个东西第一个是均值曲线对比。分别画出原始场景集和削减后场景集的平均出力曲线两者应该高度重合这是最直观的检验指标。第二个是方差对比。逐时段对比原始集和保留集的出力方差误差控制在5%以内就算合格。方差对后续优化问题的风险度量影响很大方差严重缩水的场景集算出来的方案会过于乐观。第三个是分位数对比。选取几个关键时段早晚高峰、午间光伏大发时段对比原始集和削减后场景集的5%、50%、95%分位数。如果极端场景被削减掉了分位数会明显偏离。我给自己定的经验准则是均值曲线误差小于1%方差误差小于5%分位数误差小于3%这个削减就算成功。4. 完整实操过程从原始数据到场景集的落地4.1 项目背景与数据准备我搭建的模拟项目设定一个含有风电场和光伏电站的区域微网。风电装机100兆瓦光伏装机50兆瓦。需要生成模拟未来24小时的风光联合出力场景用于后续的储能容量优化配置。数据准备分两块。风电部分通过当地测风塔历史数据分析得到了威布尔分布参数 ( k2.1 )、( c7.8 )光伏部分采集了一年的辐照度数据按小时时段分别拟合Beta分布参数上午10点对应的 ( \alpha 2.3 )( \beta 2.7 )下午2点对应的 ( \alpha 2.8 )( \beta 2.2 )。4.2 场景生成的实现细节含代码思路场景生成过程原则上可以直接用Python的NumPy完成完全不需要依赖重型科学计算平台。核心步骤是先生成带相关性的标准正态序列再逆变换到威布尔或Beta分布。下面是我在这个项目里写过的生成代码的核心片段伪代码风格换任何环境都能复现import numpy as np # 参数设置 n_scenarios 3000 n_hours 24 phi_wind 0.85 # 风速自相关系数 k_w, c_w 2.1, 7.8 # 威布尔参数 # 生成带时序相关性的风速场景 wind_scenarios np.zeros((n_scenarios, n_hours)) for i in range(n_scenarios): z np.zeros(n_hours) z[0] np.random.normal(0, 1) for t in range(1, n_hours): z[t] phi_wind * z[t-1] np.random.normal(0, 1) # 逆变换到威布尔分布 u 0.5 * (1 np.vectorize(np.math.erf)(z / np.sqrt(2))) wind_scenarios[i] c_w * (-np.log(1 - u)) ** (1 / k_w)光照的生成逻辑完全类似只是把逆变换的目标分布换成各小时自己的Beta分布并增加“夜间时段出力为0”的条件。这里有个小细节Beta分布没有解析形式的逆函数所以实际实现中我用的是scipy.stats.beta.ppf这类数值逆函数或者直接用scipy.stats.beta.rvs配合相关性处理。4.3 从风速/光照到功率的映射生成风速场景后代入第二节的风电功率转换模型得到风电出力场景矩阵光照场景则通过线性功率模型换算为光伏出力矩阵。把两套矩阵对应时段相加就得到风光联合出力场景集。实际操作中我习惯不直接相加而是分开保留两套矩阵。原因很简单后续优化任务可能是“风电储能”或“光伏储能”分开建模的合在一起的联合出力矩阵虽然方便但会丢掉各电源的独立特征。所以我是存两套矩阵需要联合场景时再临时组合。4.4 快速前向削减的核心实现我在实际项目中对比过两种削减算法这里重点展示快速前向选择的核心步骤实现因为它比K-means更需要自己动手。核心代码如下def fast_forward_selection(scenarios, probs, k): n scenarios.shape[0] selected [] remaining list(range(n)) while len(selected) k: best_j None best_dist np.inf for j in remaining: # 计算以j为代表的总距离 total_dist 0 for i in remaining: if i j: continue d np.linalg.norm(scenarios[i] - scenarios[j]) total_dist probs[i] * d if total_dist best_dist: best_dist total_dist best_j j selected.append(best_j) remaining.remove(best_j) return selected当然这个朴素的循环实现效率并不高场景数到几千个时每一步都要算全量距离耗时会拉满。实际工程项目中我会预计算全部场景对的距离矩阵并缓存这样每一步削减只需要查表和更新速度能提升一到两个数量级。4.5 削减参数选择与最终场景集输出我把初始场景数设为3000削减目标设为15个代表场景。削减后的场景概率分布不一定均匀有的场景概率可能高达0.2有的只有0.02这些概率要作为权重带进后续的优化模型。最终输出的场景集包含三个文件或表格削减后的15个风电出力场景曲线、15个光伏出力场景曲线、15个场景对应的概率值。另外我会附上一份削减质量报告包含均值误差、方差误差、分位数误差三个指标方便后续查看这组场景的可信度。5. 常见问题与排查技巧实录5.1 生成的风速曲线为什么在时间轴上疯狂振荡这是初学者最常遇到的问题原因就是忽略了时序相关性。如果你直接对每一小时独立采样得到的序列第1小时风速12米/秒、第2小时突然掉到3米/秒这显然不符合气象学常识。解决办法就是我前面说的引入AR(1)模型生成相关正态序列再逆变换到目标分布。风速的时间自相关系数越高曲线越平滑。风力大的地区夏季和冬季的 ( \phi ) 值略有差异应针对季节分别标定。5.2 削减后的概率加起来不等于1这个问题的本质是概率归一化没做好。初始采样时每个场景概率通常是 ( 1/N )概率转移过程中可能出现浮点累加误差。解决办法很简单削减后强制归一化一次让所有场景概率之和等于1。如果做完归一化后某个场景的概率是负值那说明概率转移时出了逻辑错误检查一下“替代关系”的归属判断。5.3 K值到底怎么定既要算得快又要算得准我常用的一个经验做法是从K5开始按步长5逐步增加到K50分别计算削减前后目标函数比如总购电成本或弃风率的差值。当K增加到某个值后差值变化率趋于平缓那个拐点就是合适的K值。在我经手的几个项目里这个拐点普遍出现在K10到20之间。有的文章写K可以取到几十甚至上百但这要看具体模型的可解范围。K30以上时计算时间的增长已经很可观了边际收益却很小。5.4 削减后极端场景丢失可靠性评价过于乐观用K-means就特别容易出现这个问题。大风速和高光照的极端场景在质心计算中会被“平均”掉。这种情况下我的建议有两个一是改用快速前向削减它对边界场景更友好二是采用混合策略先预留少量极端场景作为“强制保留场景”余下的场景再用聚类削减。这个方法在一些实际工程中效果显著避免了备用容量配置因为极端场景被磨平而偏低的风险。5.5 不同随机种子生成的结果差异很大怎么判断是否正常蒙特卡洛方法本来就是随机抽样不同随机种子之间的结果差异是正常的。判断标准是看统计量是否收敛而不要试图让两次运行结果完全相等。多次运行后的均值曲线应当越来越平稳各分位数也应趋于稳定。建议在报告中注明使用的随机种子和采样规模方便他人复现。如果两次运行差异巨大通常意味着样本量不够或者分布参数拟合有问题。5.6 场景削减效率太低先算距离矩阵如果场景数上到数千甚至上万快速前向削减最耗时的部分就是反复计算距离。我一开始没注意这个问题跑了3000个场景花了近半小时后来改成预计算距离矩阵时间压缩到两分钟以内。距离矩阵的存储对内存也有要求3000×3000的矩阵大概是70多MB还在可接受范围内。如果场景达到上万就要考虑更高效的数据结构或分块处理。6. 实操过程中的几点体会这套蒙特卡洛加场景削减的组合拳我在不同项目里用了很多次总结几条个人的心得。分布参数拟合别偷懒。风速的威布尔参数和光照的Beta参数直接用原始数据拟合后再采样和用文献经验值采样结果差异很大。尤其是Beta参数的矩估计数据样本少的时候波动大最好能积累至少一年的逐小时数据。做场景削减永远保留一份原始场景集。我见过有人直接把削减后的场景集当成最终分析结果不保留原始集后面要改K值或改算法时还得重新生成一遍场景非常被动。原始场景集是中间资产建议落盘存好后面做灵敏度分析、换削减算法时直接复用。分析场景集的时候重点关注代表场景对应的概率是否出现“极端占大头”的现象。如果快速前向削减后某个场景概率超过0.3说明原场景集中很多场景都长得差不多或者削减过程中某个最近邻“吸”走了太多场景的概率。这个场景在后来的优化里权重会被拉得很高如果它同时是一个比较极端的场景结果可能对决策产生不成比例的影响。项目的可复现性也别忘了。生成场景时用到的随机种子、AR(1)的自相关系数、削减的K值这些参数都建议写进项目文档。蒙特卡洛方法天然随机没有参数记录过三个月回来再看结果自己都说不清楚当时是怎么算出来的。最后还想分享一个小技巧。如果场景削减后仍然觉得K太大可以先削减到中期规模比如30个再人工合并那些形状相似且概率都较小的场景。这种“先机器削减再人工微调”的做法在高可靠性的项目里很有价值。毕竟算法再聪明也无法完全替代人对现场物理过程的理解。
阅读完成 · 觉得有帮助？
