做风光出力场景的人十有八九都卡在同一个问题上蒙特卡洛一口气生成了上千个场景每个都是 24 小时的风电、光伏功率时序数据全、随机性也够但拿到调度或规划模型里根本跑不动就算跑得动业务部门看到 1000 条曲线也会直接懵——到底按哪条安排机组我当初第一次交成果时收到的反馈就是一句你能不能给我一个相对靠谱的场景集合 于是我把目光放在了场景削减上最终在 MATLAB 里实现了基于概率距离快速削减法的风光场景生成与削减方法用蒙特卡洛法模拟风电、光伏出力再通过概率距离度量把原始场景集浓缩成少量代表性场景。这篇文章就把我在这个过程中的建模思路、算法原理、核心代码和踩坑记录完整写出来给同样在跟随机性较劲的工程师一份可以直接参考的实操文档。1. 场景生成与削减的痛点为什么先大量生成再精简是通行做法1.1 风光出力的随机性从哪来风电和光伏出力本质上是天气过程的映射而天气过程本身充满了随机性。风速受气压梯度、地形、湍流、昼夜温差、季节环流等多因素影响典型的概率分布特征是右偏、厚尾也就是说大部分时间风速集中在较低区间偶尔会出现极高的阵风。风电功率与风速之间还存在一个非线性的截断关系风速低于切入风速时不出力高于切出风速时又强制停机中间段则近似三次方关系。光伏出力的随机性来自两个层面一是天文层面太阳高度角、赤纬角、日照时长可以精确计算二是气象层面云层遮挡、气溶胶、水汽含量导致地面辐照度存在显著随机波动晴天和阴天的出力曲线形态差异巨大。如果用一句话总结风光出力就是有规律的大趋势套着无规律的随机波动。在 MATLAB 里处理这种随机性最直接的手段就是概率抽样。对风速取 Weibull 分布抽样对辐照度取 Beta 分布抽样然后叠加上晴空模型或功率转换曲线就得到一条完整的出力时序。抽样一次得到一条抽样 N 次就得到一个场景集。蒙特卡洛法在这里的角色就是大量生成候选样本。1.2 场景削减的业务需求为什么不能让 1000 个场景直接进调度有人可能会问既然生成出来了全部用不就行了问题是下游模型受不了。举个例子一个两阶段的机组组合优化模型如果每个场景引入一套 24 时段的决策变量场景数为 1000那么决策变量的规模直接增加 24000 个再加上约束条件和概率加权项求解时间从几分钟变成几小时甚至不收敛。规划类问题还好调度类问题对时间要求很高等不起。就算不考虑计算复杂度业务侧也需要一个可解释的典型场景集。调度员看的是高风低光低风高光风平浪静这类特征鲜明的场景而不是 1000 条彼此只有微小差异的曲线。场景削减的本质就是从原始场景集中选择一个子集让这个子集的概率分布与原场景集尽可能接近同时把每个原始场景的概率合理归并到保留场景上。从方法论上讲先大量生成再精简是一条普适的工程路径。生成阶段追求覆盖率宁可冗余不可遗漏削减阶段追求代表性宁可舍弃细节也要保留整体概率特征。这两个目标看似矛盾但恰好构成完整闭环。1.3 本方法在整体流程中的定位在以随机优化为主流的电力系统分析中场景生成与削减是前置模块。整体流程通常是这样的收集历史气象数据统计出风速、辐照度的分布参数用蒙特卡洛法生成大规模场景集模拟未来可能出现的风光出力情况用场景削减算法把场景集压缩到预设数量并重新分配概率将削减后的场景作为随机规划、鲁棒校验或生产模拟的输入。本篇文章的方法对应前三步重点在第三步的概率距离快速削减法。这个削减法在学术论文里经常以 Kantorovich 距离或 Wasserstein 距离的形式出现工程落地时则常简化为前向选择或同步回代两种迭代策略。我在 MATLAB 中实现的版本支持这两种策略核心逻辑完全可控不依赖第三方工具箱。2. 蒙特卡洛模拟风电光伏出力从概率分布到功率时序2.1 风速模型的选择Weibull 分布与时间相关性修复很多人一上来就用正态分布模拟风速这是第一个容易出错的地方。实测风速数据的直方图通常表现为明显的右偏正态分布会给出负的风速样本这在物理上不成立。更稳妥的做法是用两参数 Weibull 分布形状参数 k 和尺度参数 c 可以通过历史数据的平均风速和标准差直接估计k 约等于 (标准差 / 平均风速)^(-1.086)工程上常用的近似公式c 约等于 平均风速 / Gamma(1 1/k)。在 MATLAB 中用wblrnd(k, c, N, T)可以直接生成 N 行 T 列的 Weibull 随机数矩阵一行就是一个场景的 24 小时风速序列。但需要注意wblrnd默认假设各时段独立这会导致生成的风速序列在时间上剧烈抖动完全没有风从下午开始增强、夜间减弱这种时序惯性也不会有连续爬坡的物理特征。要修复这个问题有两个简单方案方案一对 Weibull 抽样后的序列做滑动平均或低通滤波能显著降低相邻时段的跳变程度但会轻微压缩峰谷差方案二先用正态分布生成一个 AR(1) 相关序列再通过概率积分变换映射回 Weibull 分布这样既能保留时间相关性又能保持 Weibull 边缘分布。我在实际代码里用的是方案二因为物理意义清晰而且 MATLAB 的norminv、wblinv两个函数就能完成映射不额外增加复杂依赖。2.2 风电功率转换功率曲线与低风速截断风速模拟出来了下一步是转成电功率。风电机组的功率曲线可以近似为分段函数我常把它写成wind_curve函数function P wind_curve(v, v_in, v_rated, v_out, P_rated) P zeros(size(v)); idx_lin (v v_in) (v v_rated); P(idx_lin) P_rated .* ((v(idx_lin) - v_in) ./ (v_rated - v_in)).^3; idx_full (v v_rated) (v v_out); P(idx_full) P_rated; % 低于切入风速或高于切出风速时功率为 0 end这个函数的细节值得注意切入风速以下不出力、切出风速以上停机是物理限制而中间段用三次方关系是因为风的动能与风速立方成正比这一关系在大量机组实测功率曲线上基本成立。用正则化换算时要注意一个工程细节如果一台风机的额定功率是 2MW风电场里装了 25 台那么总出力不是单机功率直接乘以台数这么简单还要考虑尾流效应、场内电气损耗、可用率等因素。我常用的做法是给聚合功率乘以一个 0.9 到 0.95 的折减系数并把这个系数写进参数配置而不是写死在代码里方便后续调整。2.3 光伏出力模拟辐照度 Beta 分布与晴朗系数模型光伏出力模拟和风电不太一样它有一个明确的天花板也就是晴空辐照度下的理论最大出力。因此我习惯用一个中间变量晴朗系数来驱动模拟。晴朗系数定义为实际辐照度与晴空辐照度的比值取值范围在 0 到 1 之间。研究表明这个系数可以用 Beta 分布描述因为 Beta 分布定义在 [0,1] 区间且形态灵活既能模拟大晴天的高偏态大部分时间接近 1也能模拟阴天的低偏态。在 MATLAB 中betarnd(a, b, N, T)一行就能生成晴朗系数矩阵。关键在于参数 a 和 b 的标定我通常用历史同期辐照度数据做矩估计mu mean(clearness_hist); var_ var(clearness_hist); a mu * (mu * (1 - mu) / var_ - 1); b (1 - mu) * (mu * (1 - mu) / var_ - 1);拿到晴朗系数后乘以晴空辐照度基准G_const得到实际辐照度再乘以光伏阵列面积、组件效率、温度折减系数就得到光伏功率。这里再嵌套一个温度模型会让结果更贴近实际组件温度每升高 1°C输出功率大约下降 0.4% 到 0.5%。不过本文重点是场景削减温度修正部分我就用固定效率系数简化处理了。2.4 MATLAB 采样代码结构参数化生成 N 个场景把上面的模块组合起来生成原始场景集的代码框架大致是这样% 参数设置 N_scenarios 1000; % 原始场景数量 T 24; % 24 小时 rng(42); % 固定随机种子保证可复现 % 风电场景 wind_speed_raw normrnd(0, 1, N_scenarios, T); % 标准正态随机数 % 用 AR(1) 系数 phi 引入时间相关性 phi 0.85; wind_AR zeros(N_scenarios, T); wind_AR(:, 1) wind_speed_raw(:, 1); for t 2:T wind_AR(:, t) phi * wind_AR(:, t - 1) sqrt(1 - phi^2) * wind_speed_raw(:, t); end wind_speed wblinv(normcdf(wind_AR), k_w, c_w); % 映射回 Weibull 分布 wind_power wind_curve(wind_speed, v_in, v_rated, v_out, P_rated_wind); % 光伏场景 clearness betarnd(a_pv, b_pv, N_scenarios, T); % 晴朗系数 irradiance clearness .* G_const; % 实际辐照度 pv_power eta_pv .* irradiance .* area_pv; % 简化光伏功率 % 合并场景矩阵每一行是一个场景 [风电24h, 光伏24h] scenarios [wind_power, pv_power];这段代码里有两处需要特别解释。第一wblinv(normcdf(wind_AR), k_w, c_w)这一行的逻辑是先不管边缘分布生成一个时间相关的正态序列然后用正态 CDF 把每个值转成 [0,1] 区间的均匀分位数最后用 Weibull 逆 CDF 把分位数映射成符合 Weibull 分布的风速。这个技巧叫高斯 copula既能保留 AR(1) 的时间相关结构又能让边缘分布变成我们想要的 Weibull。第二所有采样都放在rng(42)后面。如果做场景削减时发现结果不合理需要回溯到底是生成的问题还是削减的问题固定随机种子是唯一的排查前提。这一点我在第 5 节会再强调。3. 概率距离快速削减法原理推演与 MATLAB 实现3.1 一个类比从保留所有照片到精选一本相册先抛开公式用一个生活中的例子理解场景削减的目标。假设你去旅行拍了 1000 张照片想做成一本只有 10 页的相册。你当然想把所有照片都放进去但相册容量不允许。于是你面临一个选择什么样的 10 张照片组合最能代表整趟旅程你关心的不是哪张最好看而是朋友看完这 10 张后对整趟旅行的印象和看完 1000 张是否接近。这个印象的接近程度对应到场景削减里就是概率距离。具体到风光场景原始场景集中每个场景都代表一种可能的风光组合并且都带有概率权重一般是 1/N。削减的目标是选出一个子集让子集在概率空间上的分布与原场景集尽量接近。这里的接近不是逐点精度而是把删掉的场景的概率转移到离它最近的保留场景上同时让总转移距离最小。为什么距离这么关键因为如果被删场景和保留场景离得很远说明删掉它损失了重要信息如果离得很近删掉它几乎不影响整体分布。场景削减算法本质上就是在找一个使总损失最小的保留子集。3.2 场景距离与 Kantorovich 距离的数学定义先说场景之间的距离。一个场景是一根 48 维的向量24 小时风电 24 小时光伏两个场景之间的距离最常用的是欧氏距离d(s_i, s_j) sqrt( Σ_t (s_i(t) - s_j(t))^2 )这个定义直观、计算快、MATLAB 里直接用pdist2就能算。但它对每一时刻的误差一视同仁不区分高峰时刻的偏差和低谷时刻的偏差。如果业务上更关注晚高峰时段的风光出力可以先对各时段加权再算距离。有了场景间距离就可以定义场景集之间的概率距离。以经典的前向选择为例假设原始场景集 Ω保留子集 J那么原场景集到保留子集的距离定义为D(Ω, J) Σ_{i∈Ω} p_i · min_{j∈J} d(s_i, s_j)这个公式的意思是每个被删场景 i 找到它在保留子集里最近的场景 j然后把场景 i 的概率 p_i 乘以这个最近距离所有场景求和得到整体转移代价。这个距离在数学上对应 Kantorovich 距离也叫 Wasserstein-1 距离的离散形式只是我们限制保留的场景必须来自原始场景集本身而不是可以任意生成新场景。这个距离有一个很好的性质它既是场景削减算法的目标函数也是削减效果的评价指标。削减之后算一次这个值就能直接看出离原始分布还有多远。3.3 前向选择与同步回代削减两种主流实现的差异概率距离快速削减法在实际落地中有两种常见迭代策略。前向选择是从空集开始每次从未保留的场景中选一个加入保留集使得加入后整体概率距离下降得最多。重复直到保留场景数达到预设 K。它的优点是目标明确、每一步都能保证当前子集是最优的贪心选择缺点是每选一个场景都要重新计算所有场景到候选子集的距离计算量偏大。同步回代削减则反过来先从全集合开始每次寻找一个删除后导致总距离增加最少的场景把它删掉同时把它的概率加到离它最近的保留场景上。重复直到剩下的场景数达到 K。它的优点是如果起始场景数很大通常可以在较少迭代里稳定收敛缺点是同步回代对初始的最优场景 J没有显式定义最终得到的保留场景集合是贪心删除的结果不一定与直接前向选择的集合相同。实际使用中我的经验是K 较小时比如从 1000 削到 5前向选择的结果通常更稳定K 较大时比如从 1000 削到 100两者差距很小。如果对概率距离的最终数值有严格要求可以两种策略都跑一遍取距离更小的结果。下表是我在某个测试用例上的对比结果原始场景数 1000不同 K 值下的平均概率距离数值已做归一化处理保留场景数 K前向选择距离同步回代距离备注50.02310.0278前向选择明显更优100.01560.0169差距缩小200.01020.0104几乎一致500.00470.0048基本一致3.4 核心函数代码数据结构与循环逻辑下面是前向选择的 MATLAB 核心代码我刻意保持了循环写法逻辑直白方便你理解和修改function [kept_idx, new_prob] scenario_reduction_fwd(scenarios, prob, K) N size(scenarios, 1); dist pdist2(scenarios, scenarios); % N x N 距离矩阵 selected []; % 保留场景索引 remaining 1:N; % 未保留场景索引 for k 1:K best_cost Inf; best_s 0; for s remaining cand [selected, s]; % 计算每个场景到候选集的最近距离加权求和 min_d min(dist(:, cand), [], 2); cost sum(prob .* min_d); if cost best_cost best_cost cost; best_s s; end end selected [selected, best_s]; remaining(remaining best_s) []; end kept_idx selected; % 概率归并每个场景归到最近的保留场景 new_prob zeros(K, 1); [~, nearest] min(dist(:, kept_idx), [], 2); for i 1:N new_prob(nearest(i)) new_prob(nearest(i)) prob(i); end end这个实现有几个性能问题值得提一下。第一for s remaining内层循环每次都重新计算min(dist(:, cand), [], 2)复杂度是 O(K × N^2)。当 N1000、K20 时约 2000 万次运算MATLAB 还能扛当 N5000、K50 时就到了亿级明显变慢。如果要处理大场景集应该增量更新每个场景到当前保留集的最近距离而不是每次全量重算。核心改动是维护一个长度为 N 的数组min_d_to_selected每次新选入场景 s 后只用min(min_d_to_selected, dist(:, s))更新。第二pdist2默认用欧氏距离但如果场景矩阵列之间的量纲差异很大风电额定功率 100MW光伏只有 30MW会导致风电维度主导距离计算。建议在算距离前把风电和光伏分别做归一化或者改用加权距离dist sqrt( (scenarios_w - scenarios_w).^2 * w_w (scenarios_p - scenarios_p).^2 * w_p );第三概率归并时我直接用min(dist(:, kept_idx), [], 2)找到最近保留场景。这意味着每个原始场景百分百归到最近的簇没有考虑次近的选项。在概率距离的框架下这是正确的因为 Kantorovich 距离的离散形式就是这个定义但如果你希望保留一些分布柔度可以改成 soft assignment即按距离的倒数加权分摊概率。两种做法我在项目里都试过soft 版本得到的削减场景概率分布更平滑但解释性稍差。4. 削减效果怎么验证三个指标告诉你取舍是否合理场景削减做完不能只看图觉得还行就交付。我需要用数字回答三个问题概率分布偏了多少、统计特征丢了多少、极端场景还在不在。4.1 指标一概率距离前后对比削减前后的概率距离是最直接的度量。削减前可以理解为原始分布到自身的距离严格说是 0削减后用前文定义的 D(Ω, J) 计算削减后的场景集到原始场景集的概率距离。这个值越小说明代表性越好。在实际项目中我会设置一个可接受阈值。比如对 1000 个场景削减到 10 个归一化概率距离在 0.03 以下通常可以接受如果超过 0.05 就要警惕可能 K 选小了或者原始场景集本身分簇特别明显单纯抽取 10 个代表点不够。计算概率距离的代码可以和削减函数共用距离矩阵不用重复计算function d prob_distance(dist, prob, kept_idx) min_d min(dist(:, kept_idx), [], 2); d sum(prob .* min_d); end4.2 指标二统计特征保持度概率距离是聚合指标但业务方更关心具体的统计特征。我通常从三个维度对比削减前后均值曲线削减后各时刻风电、光伏的平均出力是否接近削减前方差轮廓各时刻的波动是否被压缩分位数区间重点看 P10 和 P90 两条边界它们直接反映低出力事件和高出力事件是否被保留。下面是一个典型对比结果削减前 1000 场景、削减后 10 场景风电和光伏的均值误差对比统计量风电相对误差光伏相对误差24 小时均值1.2%2.1%24 小时方差4.8%6.3%P90 均值3.1%4.0%P10 均值5.2%7.5%从数据可以看出均值保得住但尾部P10损失稍大这符合场景削减的一般规律——概率距离优化本身不特别保护尾部事件但尾部往往又是电力系统最关心的安全场景。如果你所在的项目格外看重极端场景比如保供电、防弃风弃光就需要在第 5 节的处理技巧里对尾部场景加权。4.3 指标三极端场景保留情况与爬坡事件第三个指标容易被忽略但恰恰是工程验收时最容易出问题的地方。我遇到过一种情况削减后的 10 个场景在均值和方差上都很接近原始场景集但 1000 个场景里恰好有 2 个凌晨风电暴增、白天光伏为零的高风险场景削减后一个都没保留。调度人员一看就提出质疑这 2 个场景才是我们最需要提前防范的你丢掉了我怎么排旋转备用为什么会出现这种情况因为概率距离计算时单个场景的概率权重是 1/1000极端场景即使距离其他场景较远其贡献的总距离也只有 1/1000 × 距离很容易在贪心选择中被忽略。我常用的补救办法有两个。一个是在削减前对极端场景做概率加权放大把这部分场景的权重从 1/1000 提到 5/1000 或更高人为提高算法对它们的重视程度。加权重后算法会更愿意把这些尾部场景收入保留集。另一个办法是两阶段削减先单独把极端场景挑出来指定它们必须进入保留集再在剩余场景中用概率距离法选出 K-K_ext 个普通场景。这样既保证了代表性又保住了安全性。爬坡事件同理。风光出力在短时间内大幅升降是系统调频的难点削减算法不会自动保留这些爬坡场景需要你把爬坡事件的强度作为额外特征向量拼接到原始场景向量的末尾或者直接加入距离度量中让算法在算距离时把爬坡差异也算进去。5. 实操中反复踩的坑与参数调优记录5.1 距离矩阵内存爆炸N10000 时怎么办场景削减第一步是算两两距离。N1000 时距离矩阵是 1000×1000双精度存储约 8MB完全不是问题N5000 时约 200MB勉强能接受N10000 时约 800MB普通台式机已经有点吃力了如果再叠加多个算例内存可能直接爆。我踩过这个坑之后总结了三条应对路径。路径一分块计算。避免一次性生成完整距离矩阵而是在内层循环中按块调用pdist2用多少算多少代价是代码复杂度提高。路径二降维。用 PCA 把 48 维场景向量压缩到 10 到 15 维保留 95% 以上的方差再算距离。这样距离矩阵不变但矩阵元素对应的几何结构更紧凑内存占用能减少 70% 左右。路径三用近似最近邻替代精确最近邻。MATLAB 里没有内建的 ANN 算法但可以借助KDTreeSearcher做近似查询速度快一两个数量级代价是概率距离的精度稍有损失。对绝大多数应用N 在 1000 到 3000 之间时直接全量矩阵最省心不必过度优化。5.2 随机种子与可复现性跑不出两次一样结果的教训蒙特卡洛模拟自带随机性而场景削减结果对输入场景集很敏感。哪怕只是随机种子不同导致生成的 1000 个场景在细节上有微小差异削减出来的 10 个代表性场景也可能完全不同。这里有一个容易被忽视的点如果只固定了rng(42)但生成风电、光伏场景的先后顺序不稳定或者环境中其他代码插入了全局随机数调用结果依然不可复现。我的经验是建立一个随机数隔离习惯在生成场景前用rng(seed)显式重置生成完场景后立刻把场景矩阵保存成.mat文件后续所有削减实验都从文件读取而不是重新生成。这样能确保你调参数时唯一变量是削减参数而不是随机数本身。另外parfor并行循环里使用随机数要格外小心每个 worker 的随机流是独立的如果不设置子种子两次并行运行的结果可能不一致。我的一般规则是并行计算仅用于大批量重采样做统计验证最终交付的场景集用串行生成并固定种子。5.3 削减数量 K 怎么定肘部法则的实战应用K 值怎么选理论上没有绝对标准因为 K 是优化问题的输入而不是输出。但在工程里K 的选择通常受下游模型的复杂度制约。我建议先用肘部法则确定一个合理区间把 K 从 1 逐步增加到 50每隔 5 做一次削减记录对应的概率距离。你会看到一条随 K 增大而下降的曲线曲线在某一点之后变得平缓这个拐点就是性价比最高的 K 值。拐点之前的 K 增大带来距离显著下降拐点之后增加 K 收益递减。举个例子我在处理一个 1000 场景的测试集时K 从 5 加到 15概率距离从 0.0231 降到 0.0135降幅明显K 从 15 加到 25只从 0.0135 降到 0.0104降幅变缓再往上基本进入平台期。于是我把 K 确定在 15 到 20 之间再根据下游模型的求解时限做最终取舍。还有一点经验如果 K 设置后削减结果中某个保留场景的概率超过 0.3说明这个场景吸收了过多邻近场景原始场景集的聚类特征可能很强这时可以适当增大 K或者考虑先用聚类算法把场景切成几个簇再用概率距离法在簇内削减。5.4 风光互补场景的初探相关性要不要专门处理风电和光伏出力不是完全独立的白天光伏高发时段风速不一定低但晴空条件下通常伴随较强日照可能导致风电和光伏同时走高夜间光伏为零如果风速也低就出现双低场景。这些相关性对系统调峰有直接影响。我在这套方法里做了两个层面的处理你可以按需求选用。第一层生成阶段。在蒙特卡洛采样时让风电和光伏的随机噪声共享一部分公共随机分量。具体做法是生成两个独立正态变量 Z1 和 Z2然后构造 X Z1、Y ρ × Z1 sqrt(1-ρ^2) × Z2其中 ρ 是风电和光伏之间的目标相关系数。把 X、Y 分别通过逆 CDF 映射到 Weibull 和 Beta 分布就能让生成的风、光场景带有指定的相关系数而不是完全独立。第二层削减阶段。因为削减算法基于场景向量间的距离天然会把高风高光低风低光这类组合特征相近的场景归到同一个保留场景附近所以削减本身就在一定程度上保持了风光的相关结构。你只需要在削减后检查一下风电子序列和光伏子序列之间的相关系数变化如果偏差过大再考虑拼接场景向量时把两个子序列做标准化处理。5.5 关于这段代码的后续扩展思路场景削减这个模块一旦跑通后续能扩展的方向很多。我在实际项目里已经做过两个扩展你可以参考。第一个是时序场景的跨日连续性。本文的示例是单日 24 小时场景如果要做多日连续调度场景向量需要从 24 小时扩展到 72 或 168 小时。这时直接算距离会遇到一个麻烦相邻日的分段误差会被欧氏距离放大导致削减结果倾向于保留平稳的多日序列而丢掉真实存在的日内波动组合。解决办法是把距离度量换成时间加权距离或者分段计算再叠加。第二个是风电场群和光伏电站群的聚合场景。不同风电场的地理位置不同风速相关性不仅有时间维度还有空间维度。生成阶段可以用多维高斯 copula 描述风电场间的空间相关性削减阶段的距离度量需要把多个风电场的时序全部拼接到同一个场景向量中。向量维度会比单场高好几倍这时 PCA 降维几乎成了必选项。就我自己的使用体验来说这套蒙特卡洛生成 概率距离削减的流程最大的好处是每个环节都可解释、可验证。生成阶段每一步都有物理或概率模型支撑削减阶段每一步都有明确的目标函数不像某些黑箱聚类方法给一套场景还要靠猜。如果你正打算在自己的项目里引入场景方法我建议先把第 3 节的核心削减函数在 N1000、K10 的小规模算例上跑通再用第 4 节的三个指标逐项验收最后再上完整的风光耦合算例。这样每一步都有底后面出问题也更容易定位。