做配电网规划、微电网调度或者储能容量配置的人对“风光及负荷多场景随机生成与缩减”这个说法应该不会陌生。很多时候优化模型本身并不复杂真正让人头疼的是喂给模型的那组场景——风电出力、光伏出力、负荷曲线三者叠加不确定性不是简单相加而是一层摞一层根本没法用一条“典型日曲线”概括。风大可能同时伴随着光伏出力波动负荷低谷时段的光伏为零这些耦合关系如果处理不好模型算出来的结果基本没法用。我最初处理这类问题时直接拿全年8760小时的数据去跑随机优化模型规模直接爆炸解一次要几个小时而且结果在重复运行时还不稳定。后来才意识到问题不在求解器而在数据前处理用多场景随机生成与缩减技术把海量历史数据转化为一组数量少、概率明确、统计特性接近原始分布的代表性场景这才是随机优化落地的前提。这套基于Matlab的实现覆盖了从历史数据输入、概率场景生成、聚类缩减到结果验证的完整链路不依赖付费工具箱核心逻辑清晰适合做规划研究的工程师和学生参考。1. 随机优化真正需要的不是“典型日”而是一组带概率的场景1.1 单条曲线无法表达风光荷的三重不确定性很多刚接触这个方向的人会问我能不能取一个平均风速日、一个平均光照日、一个平均负荷日组合成“典型场景”我的回答是可以但结果多半没有参考价值。原因很直接——取平均值的过程把波动抹掉了而波动恰恰是风光荷不确定性的核心。风电出力呈现强随机性风速的Weibull分布偏斜明显一个高风速场景和一个低风速场景对优化结果的影响完全不是“平均风速”能代表的光伏出力受云层遮挡、天气状态跳变影响峰值时段可能出现剧烈跌落负荷则叠加了工作日/休息日、季节、气温等多种周期性因素。三者组合在一起真实的联合分布可能是多模态的平均值场景落在概率密度很低的区域代表不了任何典型情况。所以随机优化需要的是一组覆盖关键概率区间的场景集合每个场景带一个发生概率。这组场景既要包含大概率的高频状态也要包含概率较低但影响极大的极端状态两者在优化目标中都会被加权。1.2 场景缩减不是偷懒而是概率分布意义上的逼近这里必须澄清一个常见误解场景缩减不是“随便挑几个日子凑数”而是在概率分布层面做近似。理论上场景越多对原始不确定性的刻画越精细但场景一多优化模型的变量和约束成倍增加随机规划的一个场景就要复制一套决策变量。以两阶段随机规划为例如果第一阶段有100个决策变量每个场景第二阶段有200个变量300个场景就会产生6万个第二阶段变量求解规模直接失控。数学上缩减的本质是给定原始场景集S和缩减场景集S寻找S上的概率分布使两者之间的Kantorovich概率距离最小。通俗解释就是用“搬移概率质量”的成本来衡量两个分布的距离——把原始场景的概率质量移动到缩减场景上移动距离越大损失越大缩减算法的目标就是让这个总移动成本最小。1.3 从历史数据到可用场景集的完整链路在我做的Matlab实现里完整流程是四步走数据准备整理历史的风速、辐照度、气温、负荷数据统一时间分辨率和单位处理缺失值按季节或时段分段。随机生成基于历史数据的统计特征利用参数分布拟合、时间序列模型和相关矩阵生成N个原始随机场景通常300到1000个。场景缩减用聚类或快速前向选择算法把N个场景缩减为K个代表场景通常5到20个并重新计算每个缩减场景的概率。结果验证对比缩减前后场景集的均值、方差、极值分布和时序相关性确认缩减损失在可接受范围内。这四步看着简单但每一步都有不少细节坑。接下来我按这几个环节展开说。2. 场景生成的三条路线以及我在Matlab代码里最终怎么取舍2.1 参数分布直接抽样数学简单但容易生成“一眼假”的场景最直觉的做法是风速用Weibull分布拟合辐照度用Beta分布拟合负荷用正态分布拟合然后分别独立抽样拼成一个场景。代码三五行就能写完但生成的场景几乎一定会出问题——独立性假设不符合物理规律。实际的风速时间序列有强自相关性上一小时风速高下一小时风速往往也不低光伏出力在晴天条件下是平滑的钟形曲线不会出现逐小时剧烈跳变独立抽样生成的序列会像白噪声一样高频振荡这种场景放进调度模型会诱导模型做出过度保守或明显不合实际的决策。如果你只是想快速测试一个算法框架独立抽样可以应付但要做正经的规划分析这条路走不通。2.2 时间序列模型把风速的惯性、负荷的周期性强行加回来解决独立抽样问题的自然思路是引入时序相关性。常用方法有两种我都试过第一种是自回归滑动平均模型。以风速为例先拟合确定性部分比如用日平均风速的周期性曲线然后对残差序列建立ARMA(p,q)模型。ARMA的好处是能保留短时相关性生成的风速曲线看起来“像真的”——有持续的爬升和下降段不会乱跳。但ARMA只能捕捉线性相关对风速的非线性波动比如阵风尖峰表现一般。第二种是马尔可夫链。把风速划分成若干个状态区间比如按分位数分成5到10个状态统计状态转移概率矩阵然后按转移概率生成状态序列。马尔可夫链天然适合处理离散状态转移比如晴转多云再转阴的天气过程对光伏出力的模拟效果比ARMA更好。缺点是要构造状态转移矩阵需要较多历史数据否则转移概率估计不稳。实际项目中我见的比较多的是ARMA和马尔可夫链混合使用风速用ARMA天气状态用马尔可夫链负荷用周期项加ARMA残差。2.3 多变量相关性处理风光荷不是各自独立的随机数这是很多人容易忽略的一环。同一时刻的风、光、负荷之间是存在统计相关性的——大风天往往云层厚、光伏出力低极端高温天负荷飙升但风速可能偏低早晚时段负荷爬坡和光伏出力曲线存在天然的时序错位。处理多变量相关性最严谨的工具是Copula函数可以用一个相关结构把多个边缘分布粘在一起。Copula理论本身不难难在Matlab里没有官方的高斯Copula拟合函数需要自己写经验Copula或者用辅助函数代码量不小。对大多数规划场景我采用更简单的折中方案估计原始数据的相关系数矩阵然后用Cholesky分解对独立抽样序列施加线性相关性。具体做法是先生成互不相关的随机矩阵乘以相关系数矩阵的Cholesky因子再做逆变换映射回目标分布。这个方法只能近似线性相关无法刻画尾部相依性但对风-光-荷这种线性相关为主的场景已经够用。2.4 我在代码里采用的生成组合确定性预测随机扰动时序耦合综合考虑拟合精度、计算成本和实现复杂度我这套代码最终采用的组合方案是把历史数据按季节和工作日/休息日分组分别求每个时刻的均值得到确定性基准曲线。对基准曲线的残差按变量分别拟合分布同时估计残差的相关矩阵。利用ARMA(1,1)模型为风速和负荷生成时序相关随机残差光伏残差让马尔可夫天气状态控制。通过Cholesky分解引入风光荷之间的同期相关。把确定性基准、自回归残差、相关性扰动叠加起来生成N个原始场景。这套组合的优点是物理意义清晰生成场景的波动形态贴近真实数据。缺点是参数估计依赖历史数据质量如果历史数据只有半年季节分组后每段样本很少ARMA系数会估计得不稳。我的处理办法是加大季节性分组的粗粒度先按两季制冷季和供暖季跑通后续再细化。3. 场景缩减原理从上百个原始场景到十几个代表场景到底“缩”掉了什么3.1 缩减的数学目标最小化与原场景集的概率距离生成阶段产生的原始场景数量通常会设到300个以上为的是充分覆盖不确定性空间。缩减阶段要考虑的是选K个场景重新分配概率使缩减后场景集和原始场景集在概率分布上尽可能接近。前面提到这个“接近程度”用Kantorovich距离度量。直观理解可以这样想原始场景是散在地上的N堆沙子每个场景有对应的“质量”概率缩减后要用K个桶去装这些沙子桶的位置可以移动移动每堆沙子的成本等于沙子质量乘以移动距离最优缩减就是找到K个桶的位置和每个桶的沙子质量让总搬运成本最低。这就是缩减的核心数学模型。所有具体算法本质上都是这个最小化问题的近似求解。3.2 三种主流缩减手段的效果与选择我在Matlab里实现和测试过三种缩减方法各有明显倾向性。先列个表缩减方法输出场景类型计算复杂度适用场景主要缺陷K-means聚类类均值中心可能是合成场景低适合大规模数据离散场景较多、对精度要求不极端苛刻时中心不是真实历史场景可能过度平滑K-medoids聚类类内真实样本保证可回溯较高需要计算全距离矩阵要求缩减场景必须来自原始数据时时间序列数据量大时内存和耗时都很高快速前向选择原始场景子集概率重新分配中等K值较小时很快概率距离优先、保留极端场景时算法是贪心的不能保证全局最优K-means三千个场景几十毫秒就能缩完但它输出的是簇中心的均值序列。时序数据做均值后峰会变矮、谷会变浅典型的“平均化失真”缩减场景会丢掉极端值。K-medoids每个簇里挑一个最接近中心的真实样本能保证场景是历史真实出现的形态但矩阵是N×N的距离矩阵N500时就要存12万个距离值内存和时间压力陡然上升。我在代码里把后二种都开放出来了用户可以在函数参数里指定。对大多数场景我推荐K-medoids因为缩减场景保留了真实的时序形态后续做优化分析时更容易向非专业人士解释。3.3 时序场景不能被当成普通多维点来处理这里有个非常典型的坑。每个场景的长度是T比如24小时或168小时场景本身就是一个长度为T的向量。T越大维度越高而高维空间中的欧氏距离区分度会下降——所有点之间的距离都差不多聚类算法很难有效工作。所以在做聚类缩减前我习惯先做一步标准化和主成分压缩。先用z-score把风速、光伏、负荷三个量的量纲统一到同尺度然后拼成一个3T维的向量再对全部场景做主成分分析只保留累计方差贡献率超过95%的前几个主成分最后在主成分空间里做K-means或K-medoids。这样既避免了量纲问题又压缩了无效维度。如果你不想用主成分压缩至少也要做标准化。不然负荷是几十到几百千瓦的量级光伏是零到几千瓦的量级欧氏距离基本会被负荷主导风光的波动对聚类结果的影响会非常小。3.4 用Matlab实现K-means缩减的完整步骤与关键代码以Matlab自带的kmeans函数为例核心流程如下% scenarios: N x T 矩阵N个原始场景每个场景长度为T % 先做z-score标准化 scenStd (scenarios - mean(scenarios, 1)) ./ std(scenarios, [], 1); scenStd(isnan(scenStd)) 0; % 主成分压缩保留95%方差 [coeff, score, ~, ~, explained] pca(scenStd); cumVar cumsum(explained); idxKeep find(cumVar 95, 1, first); scenPca score(:, 1:idxKeep); % K-means聚类用kmeans初始化提高稳定性 k 10; rng(42); [clusterId, clusterCenter] kmeans(scenPca, k, Start, plus, ... MaxIter, 500, Replicates, 20); % 重新计算每个簇的概率 prob accumarray(clusterId, 1, [k, 1]) / N; % 找到每个簇中离中心最近的原始场景作为缩减场景 reducedIdx zeros(k, 1); for i 1:k pointsInCluster find(clusterId i); [~, nearest] min(sum((scenPca(pointsInCluster, :) ... - clusterCenter(i, :)).^2, 2)); reducedIdx(i) pointsInCluster(nearest); end reducedScenarios scenarios(reducedIdx, :);注意第7行我把随机数种子固定成了42这很重要。如果不固定每次运行聚类结果都会变下游优化结果也就无法复现。第12行用20次重复跑取最优解能有效规避初值敏感问题。第23到27行做的是把K-means的簇心结果转成K-medoids式的真实样本输出——这是我从K-means里“抢救”真实场景的手段很实用。4. 代码结构设计与参数选择哪些值值得调哪些不该乱动4.1 输入数据格式与接口约定我这套代码的输入约定很直接——一个CSV文件三列时间序列第一列风速第二列辐照度第三列负荷行数是研究期内的时段数。如果原始数据是15分钟分辨率的但你要做的是小时级日前调度记得先做聚合平均否则生成场景的时间尺度和优化模型对不上这个错很隐蔽。数据准备的三个硬性要求时间序列必须等间隔缺失值要提前处理不能留NaN否则标准化和PCA步骤直接报错。量纲差异不要在路上纠结代码内部会做标准化处理但原始单位别混用风速不要一会儿m/s一会儿km/h。尽量保留至少一年的完整数据。少于一年季节性拟合没意义生成出的场景可能在冬夏之间随机漂移物理性很差。4.2 参数面板里每个字段的含义和推荐取值我把参数集中在一个结构体里方便批量调参参数名含义推荐取值说明timeRes时间分辨率小时1或0.25必须和优化模型一致seasonSplit季节分组数2或4数据量不足时选2nOriginal原始生成场景数300到1000越多越稳定耗时也越长nReduced缩减目标场景数5到20越小越激进损失越大armaOrderARMA模型阶数[1, 1]对风速和负荷够用corrMethod相关性处理方法chol 或 copula默认用Cholesky速度快randomSeed随机数种子固定整数必须固定否则结果不可复现pcaVariancePCA保留方差比例95降低维度去除噪声clusterMethod缩减方法kmedoids也支持 kmeans 和 forward新手最容易犯的错误是把nReduced设得太小。我见过有人直接把1000个场景缩减成2个然后拿两个场景去做储能规划结果几乎是拍脑袋。两个场景的概率再准确也覆盖不了风光的季节性差异和日间波动。通常我先按5、10、15、20几档做灵敏度分析看优化目标对场景数什么时候收敛再选拐点附近的K值。4.3 核心函数的调用关系与输出结果组织代码的主体是3个函数加1个主脚本% 主脚本示例 histData readtable(history.csv); params defaultParams(); scenSet generateScenarios(histData, params); [reducedScen, prob] reduceScenarios(scenSet, params); plotScenario(reducedScen, prob, load(typicalDay.mat)); checkReductionQuality(scenSet, reducedScen, prob);generateScenarios负责做数据分组、分布拟合、ARMA残差生成和相关性注入输出是一个N×T×3的三维数组三个通道分别对应风、光、荷。reduceScenarios内部调用标准化的PCA和缩减方法输出缩减后的场景集和对应的概率向量。checkReductionQuality是验证函数计算覆盖度和概率偏差下面第5章会细说。我特别提醒一点输出场景矩阵的排列顺序要统一。我见过有人生成的时候场景维度是N×T×3可下游优化模型读进来的时候按N×3×T解析结果曲线全乱了还很不容易察觉。代码里最好在输出前加一句断言明确矩阵的维度定义。4.4 随机数种子最容易忽略的可复现性问题做随机生成与聚类最怕结果不可复现。今天跑出来的场景和昨天不一样优化结果整整差了几个百分点然后你会怀疑模型写错了查半天发现只是随机数种子没固定。这是这类代码最常见的假问题。我的建议是随机数种子必须作为参数显式暴露出来而且在脚本开头就执行rng(params.randomSeed)。如果你在主循环里用了并行计算还要注意每个worker的随机数流也要独立设置否则并行下的结果和串行下的不一致。简化做法是固定用单线程跑生成步骤缩减步骤再并行。一个额外的技巧把生成场景和缩减场景分别保存成独立的.mat文件文件名带上参数哈希或者时间戳。这样即使后面改了参数旧结果也能追溯。做科研和写报告时可追溯性比计算速度更值钱。5. 实操中反复踩到的坑以及一套管用的结果验证方法5.1 场景缩减后丢失极端工况优化结果明显偏乐观我最早用纯K-means缩减时踩过一个大坑。缩减后的场景数量只有10个聚类算法倾向于把样本密集的区域拆得更细而概率极低的极端场景——比如持续3天无风同时光伏为零的冬季夜晚——会被合并到附近的大类里。这些极端工况在概率上占比极小但一旦发生对储能容量配置和可靠性分析的影响极大。观察结果就是用缩减场景算出来的期望运行成本偏低系统看起来“很稳”实际一到极端天气就兜不住。解决思路有两个层面。第一把极端场景单独保护起来在聚类之前先识别出极端事件样本比如风速低于阈值的场景让它们强制占据几个缩减名额。第二在覆盖度和稀疏度指标上做硬性约束防止缩减场景全部挤在概率中心区域。5.2 聚类初值敏感导致结果抖动用kmeans和多次启动解决K-means本身对初始簇中心很敏感初始点选得不好很容易陷入局部最优。第一次缩减到10个场景A场景被保留了第二次同样的输入A场景没了B场景顶上。如果你拿着这两组场景分别去跑优化得到的储能容量配置可能差20%。这是聚类算法固有缺陷规避手段就是第3.4节代码里的两行设置Start, plus和Replicates, 20。前者用K-means初始化大幅提升初始点质量后者用20组不同的初始点跑20次保留所有运行中组内总距离最小的一次。实测下来这个配置基本消除了缩减结果在重复运行时的抖动。代价是计算时间增加大约10倍但原始场景只有几百上千个这点时间完全可以接受。5.3 用覆盖度、稀疏度两个指标客观评价缩减质量缩减质量不能靠“看起来差不多”来评价要有量化指标。我长期用两个指标组合判断一个管“准”一个管“全”覆盖度衡量缩减后场景对原始场景的覆盖能力计算方法是每个原始场景到最近缩减场景的距离的平均值。这个值越小说明缩减场景整体越接近原始分布。稀疏度衡量缩减场景之间的差异程度计算缩减场景两两之间距离的平均值。这个值越大说明缩减场景相互之间越分散重复覆盖越少。实际操作中两个指标要配合着看。如果覆盖度小但稀疏度也小说明缩减场景扎堆了虽然每一块覆盖得都准但场景间的多样性不足。如果覆盖度大而稀疏度大说明缩减场景太散对原始分布的代表性差。理想的缩减结果是覆盖度尽量小同时稀疏度尽量大两者的平衡点通常出现在K从5到15的区间内。5.4 把缩减场景放进优化模型前先做一个偏差体检最后一个建议是缩减结果别直接用先做一个偏差体检。体检方法很简单——拿原始场景全集和缩减场景集分别去跑同一个简化优化模型对比两者在最优解上的目标函数值和关键决策变量。偏差在3%以内说明缩减质量可接受超过5%说明信息损失过多需要调整K值或缩减方法。我在实际项目里习惯用储能容量配置作为体检模型。这类优化问题对场景特性很敏感而且目标函数的经济含义直观决策者和评审人都容易理解。如果缩减后场景集的容量配置和原始场景集相差无几那么这套缩减结果就可以放心用于更复杂的规划模型。做这个体检时要注意两个优化模型的结构必须完全一致只有输入场景不同否则偏差来源说不清楚。第一遍跑的时候先不固定随机数种子跑三次看波动范围第二遍再固定种子确保每次体检结果一致。这样既验证了缩减算法又验证了优化模型本身的稳定性。最后再分享一个实操技巧如果你做的是多节点配电网规划场景缩减不能在系统层面一次性做完。我建议按关键节点分别生成和缩减场景再检查节点之间的场景相关性是否被破坏。只做全局缩减很容易把局部节点的负荷尖峰细节抹平。这个坑花了我一周时间才发现希望对你有帮助。