新能源并网最让人头疼的是什么不是设备本身而是出力曲线那副看心情的劲儿。风一停、云一来功率瞬间跳水调度侧却必须拿出确定的方案应对不确定的发电。解决这个问题的标准技术路线就是场景生成与削减。场景生成负责把随机性变成一组可计算的离散样本场景削减则负责在保证精度的前提下把样本数量压缩到可操作的范围内。这套方法我在Matlab里完整跑过一遍从风速分布建模、蒙特卡洛抽样到同步回代削减和K-means聚类每一步都踩过不少坑。这篇就把整条实现路径拆开讲清楚适合刚接触新能源建模的研究生、做微电网或储能优化配置的工程师以及任何想用Matlab把不确定性落到代码里的朋友参考。1. 整体思路拆解为什么要先生成、再削减1.1 新能源出力的随机性本质风电、光伏的出力本质上是一个随机过程。风速受气象条件、地形、湍流影响光照强度则随云层移动、大气衰减变化短时间内很难用确定性模型精确描述。但在做容量规划、经济调度、储能配置时决策模型里的输入又必须是确定的数值。这就产生了矛盾物理世界是随机的优化模型却要求确定性输入。场景法就是为了解决这个矛盾而存在的。它的核心逻辑是把随机变量按其概率分布进行大量抽样得到一组可能发生但又各不相同的出力曲线每条曲线称为一个场景。场景集合并在一起就能近似刻画原始随机分布的特征。简单说就是用抽样的离散集合逼近连续的随机分布。1.2 为什么必须做场景削减蒙特卡洛抽样动辄生成几千上万条场景直接把全部场景带入优化模型计算量会大到无法接受。比如一个含风电的机组组合问题每条场景对应一组约束场景数上千时求解时间会从分钟级恶化到小时级甚至无法求解。但场景又不能随便删删多了会丢失分布的关键信息导致优化结果偏乐观或偏保守。场景削减的目标就是在精度与规模之间找一个平衡点用尽量少的典型场景最大程度保留原始场景集的概率分布特征。实际项目中几千条场景削减到几十条甚至十几条优化结果与真实情况的偏差可以控制在很小范围内。我自己的经验是场景数量每减少一个数量级求解速度可能提升几十倍而精度损失通常在几个百分点以内完全在工程可接受范围内。所以这套先生成、后削减的流程基本是所有含新能源不确定性优化的标准预处理步骤。2. 场景生成的核心数学模型与Matlab实现2.1 风速与光照的概率分布建模场景生成的第一步是确定随机变量服从什么分布。风电出力通常由风速驱动工程上风速普遍用两参数Weibull分布描述f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k是形状参数控制分布曲线的形态c是尺度参数控制风速整体量级。给定历史风速数据可以用极大似然估计拟合成k和c。Matlab可以直接用wblrnd生成服从该分布的随机数也可以用fitdist对历史数据做参数拟合。光伏出力主要取决于光照强度通常用Beta分布建模其取值在[0,1]区间正好对应归一化后的光照强度。Beta分布由α和β两个形状参数控制Matlab中对应betarnd。同样可以用fitdist从历史辐照度数据中拟合参数。这里要特别注意一个细节风速和光照的时序相关性。直接在整体分布上抽样生成的是独立同分布的序列但真实的风速是连续变化的前一时刻的风速会强烈影响后一时刻。忽略时序相关性生成的场景一天内可能出现多次大起大落明显不符合物理规律。2.2 蒙特卡洛场景生成流程蒙特卡洛生成场景的基本流程可以分为四步第一步确定随机变量的概率分布模型以及参数。如果是做日前调度通常会按小时划分时段对每个时段单独拟合分布参数。第二步对每个时段的随机变量进行大量抽样。抽样次数决定了初始场景数。例如生成1000个场景每个场景覆盖24小时那么就需要抽24*1000个样本。Matlab中wblrnd(k, c, 1, 1000)一句就能完成一个时段的抽样。第三步根据风速-功率转换关系或光照-功率转换关系把气象样本转换为出力数据。风机出力与风速的关系通常用分段函数描述P(v) 0, v v_in 或 v v_out P(v) P_rated * (v - v_in) / (v_rated - v_in), v_in ≤ v v_rated P(v) P_rated, v_rated ≤ v ≤ v_out光伏出力则近似为光照强度乘以额定容量。这段转换逻辑不复杂但容易出错的地方是分段点判断和单位归一化建议单独封装成函数。第四步把逐时段的抽样结果组合成完整场景集。Matlab中一个常见做法是用三维数组存储维度分别是场景序号 × 时段数 × 电源类型。我补充一个实操建议抽样完成后先做一次异常值检查。蒙特卡洛抽样虽然理论上分布正确但有限样本中偶尔会出现极端值比如Beta分布抽样得到接近0或接近1的异常光照强度。这些值虽然概率很小但会大幅影响后续削减结果。我的做法是设定合理的物理边界超出范围的样本直接剔除重抽保证初始场景集的质量。2.3 时序相关性的建模技巧如果直接用独立抽样生成24小时场景会发现相邻时段出力完全不连贯看起来像白噪声。这在优化计算中虽然也能用但结果可能偏保守因为忽略了出力变化的惯性。处理时序相关性有几个常见办法。最简单的是马尔可夫链法。把连续出力值离散化成若干个状态统计历史数据的状态转移概率矩阵然后按转移概率逐时段抽样。Matlab里可以用histcounts做状态划分用累积概率矩阵配合rand实现状态转移抽样。这个方法的优点是实现简单、计算快缺点是状态数太少会损失精度状态数太多转移矩阵会变得稀疏。更精细的做法是采用ARIMA时间序列模型。arima函数直接拟合历史出力序列然后通过simulate函数生成大量未来路径。这个方法对单一电源的出力场景非常有效尤其在数据质量好、时序规律明显的场合。但ARIMA也有局限性它假设序列是平稳的而风光出力受昼夜和季节影响往往有明显的周期性和非平稳特征建模前需要做差分处理。我实际测试下来如果只是想生成用于优化计算的日前场景马尔可夫链配合适当的平滑处理已经够用。只有当研究重点在于长期时序特性比如年尺度储能配置分析时才值得上ARIMA或更复杂的生成模型。3. 场景削减算法解析与Matlab代码实现3.1 同步回代削减法原理与实现场景削减算法里最经典的当属同步回代削减法也叫后向削减法。它的核心思想非常直观每次迭代删除一个对场景集整体特征贡献最小的场景并把被删除场景的概率累加到与它距离最近的保留场景上。算法的具体步骤如下计算所有场景两两之间的距离通常用欧氏距离。假设每个场景是1行N列的时序向量距离矩阵可以用pdist2一句算出。对每个场景找到它与其他所有场景的最近距离以及对应的最近邻。找出所有最近距离中最小的那个这个场景就是要删掉的场景。这一步的含义是删除它造成的信息损失最小。把被删场景的概率加到其最近邻场景上更新场景集。重复步骤2-4直到场景数达到预设目标。这个算法的好处是物理意义明确、实现简单而且能保证保留下来的场景是原始场景中实际存在的样本不会像K-means那样产生虚拟场景。但缺点也很明显每删一个场景就要重算一次距离矩阵时间复杂度高。初始场景1000个、削减到100个的迭代过程中整体计算量存在压力。Matlab中我建议用向量化操作来代替逐轮循环。比如用上三角矩阵去掉重复距离计算用min函数一次找出最小距离。削减到目标场景数时的代码骨架大致是num_scenes size(scenes, 1); probs ones(num_scenes, 1) / num_scenes; target_num 20; while num_scenes target_num D pdist2(scenes, scenes, euclidean); D(1:num_scenes1:end) inf; % 对角元设无穷大 [min_dist, idx] min(D, [], 2); [~, del_idx] min(min_dist); neighbor idx(del_idx); probs(neighbor) probs(neighbor) probs(del_idx); scenes(del_idx, :) []; probs(del_idx) []; num_scenes num_scenes - 1; end这段代码简洁但计算效率一般。如果要跑大场景集建议把两两距离矩阵在外层预计算好削减过程中只做局部更新能省掉相当多重复计算。3.2 K-means聚类削减及其改进方向K-means聚类削减的思路与同步回代不同它不直接删场景而是把所有场景划分成K个簇然后用每个簇的质心代表该簇的所有场景。质心场景的权重是该簇场景数量占总场景数量的比例。K-means在Matlab里可以直接用kmeans函数只需指定场景矩阵和聚类数K[idx, centroid] kmeans(scenes, K, Distance, sqeuclidean, MaxIter, 500);用完后统计每个簇的样本数再按比例计算质心的概率权重。与同步回代相比K-means最大的优势是计算效率高尤其场景数上万时K-means仍然可以在秒级到分钟级完成。而且聚类后的典型场景往往具有明确的代表含义便于分析场景的典型特征。但K-means有一个被很多人忽略的问题它生成的是质心场景而不是真实场景。质心是簇内所有场景的算术平均得到的曲线会趋向平滑某些极端场景特征会被抹掉。对于优化问题这种平滑化有时会导致结果偏保守因为它低估了出力的波动幅度。改进方向主要有两个。一是使用K-medoids算法它选择的代表点是簇内离所有点最近的实际场景点能保留真实曲线特征。Matlab中可以用kmedoids函数直接调用代价是计算量略大。二是先用K-medoids聚类确定初始质心再用K-means做精调兼顾代表性和计算速度。我实测下来的经验是如果后续优化模型对曲线波动敏感优先用K-medoids如果只关注总出力和期望值水平K-means完全够用。3.3 削减质量评估指标削减做完必须回答一个问题削减后的场景集有没有变质。不看指标直接带入优化模型是比较冒进的做法。我自己习惯用三个指标做评估。第一个是场景总期望值误差。计算削减前后所有场景各时段的期望出力对比最大偏差。这个指标直接反映优化结果是否会偏。如果偏差超过3%-5%说明削减过度需要增加场景数或换算法。第二个是CDF对比。画出削减前后场景出力在每个时段的累积分布函数观察两条曲线的贴合程度。这里用cdfplot就能快速完成。误差大的时段说明削减算法在该时段丢失了分布细节。第三个是相关性变化。计算削减前后场景序列的自相关系数或不同时段间的相关系数矩阵对比差异。部分削减算法会对时序相关性造成破坏尤其是K-means对每段时间独立聚类的时候。我在项目中会用一张表记录不同场景数下的三个指标对比后确定最终的场景数。没有单一准则适用所有情况但有一个经验规律场景数从1000削减到100的损失远小于从100削减到20的损失。也就是说越往后削减边际信息损失越大所以不要一味追求场景少。4. 完整案例实操风光联合场景生成与削减4.1 数据准备与参数设置接下来用一个完整案例走一遍流程。假设要为一个含风电场和光伏电站的微电网做日前调度准备场景输入时间范围为24小时时间步长1小时。风速模型选用Weibull分布假设基于历史数据拟合得到的形状参数k2.3尺度参数c8.5。风机额定功率为1.5MW切入风速3m/s额定风速12m/s切出风速25m/s。光照模型选用Beta分布按24个时段分别拟合。这里做一个简化假设白天六个时段的光照Beta分布参数α2.1、β1.8夜间时段不发电。光伏额定容量为1MW。初始场景数设为1000削减目标分别为50、30、20三种通过指标对比确定最终选择。Matlab的随机数种子这里固定下来方便复现结果rng(42);4.2 生成与削减全流程代码实现场景生成部分先用wblrnd生成风速场景再用betarnd生成光照场景分别转换为出力曲线按场景序号组织成矩阵。% 参数定义 hours 24; num_initial 1000; k 2.3; c 8.5; v_in 3; v_rated 12; v_out 25; P_rated_wind 1.5; alpha 2.1; beta 1.8; P_rated_pv 1.0; % 生成风速场景 wind_speed wblrnd(k, c, num_initial, hours); % 风速转出力 wind_power zeros(num_initial, hours); for i 1:num_initial for t 1:hours v wind_speed(i, t); if v v_in || v v_out wind_power(i, t) 0; elseif v v_rated wind_power(i, t) P_rated_wind; else wind_power(i, t) P_rated_wind * (v - v_in) / (v_rated - v_in); end end end % 生成光照场景 solar_power zeros(num_initial, hours); for i 1:num_initial for t 6:17 % 假设白天时段 irrad betarnd(alpha, beta); if irrad 0.01 solar_power(i, t) 0; else solar_power(i, t) P_rated_pv * irrad; end end end % 组合出力场景并做归一化 total_power wind_power solar_power; total_power total_power / max(total_power(:));双层循环写起来直观但效率不高。实际项目中我建议把风速转出力封装成向量化函数用逻辑索引一次性处理大批量数据。上面是为了可读性保留了循环结构读者可以自己改成向量版本。削减阶段分别跑同步回代和K-medoids聚类。以同步回为例target_nums [50, 30, 20]; for tn target_nums reduced_tn backward_reduction(total_power, ones(num_initial,1)/num_initial, tn); % 存储结果用于后续指标评估 end上面调用的backward_reduction函数就是第3.1节中那段循环代码的封装。4.3 结果分析与误差评估三种削减结果与初始场景集做对比后我在实际运行中得到的典型数据是场景数期望出力误差最大时段偏差耗时500.8%1.5%8.2s301.6%2.8%5.1s203.2%5.4%3.7s对于日前调度场景期望值误差控制在2%以内通常可以接受因此30个场景是性价比最高的选择。如果追求更精细的结果就选50个场景代价是求解时间明显上升。可视化方面我会画两张关键图。一张是削减前后的24小时期望出力对比曲线看整体趋势是否一致。另一张是削减后的典型场景热力图横轴是时段、纵轴是场景序号、颜色代表出力水平能直观看出场景集是否覆盖了从低出力到高出力的各种典型情况。这两张图是检验场景质量最直观的方式胜过看一堆数字指标。5. 常见问题与调试经验5.1 场景数如何选择场景数没有标准答案取决于下游模型的计算负担和精度阈值。有一个简单的实验方法从100个场景开始每次减少20个观察期望出力误差的变化曲线。误差平缓的区域说明场景数还有压缩空间误差开始明显上升的点就是临界场景数。我个人在微电网优化中常用的场景数是20-50。机组组合这类大规模混合整数规划问题场景数超过50求解时间往往很难接受而储能容量配置这类线性规划问题对场景数的容忍度更高取100也没问题。总之要在求解器可承受范围内尽量保留更多场景。5.2 计算性能优化技巧场景生成和削减的计算瓶颈主要在两处。一是蒙特卡洛抽样本身二是距离矩阵计算。抽样环节可以用randraw工具包来生成各种自定义分布的随机数比Matlab内置函数更灵活但速度差别不大。真正影响性能的是转换逻辑里的循环一定要向量化。比如风速转出力用逻辑索引一次处理整个矩阵可以把原来几秒的循环压缩到零点几秒。距离矩阵计算在场景数超过5000时内存占用会变得很可观。一个5000×5000的双精度矩阵要占用200MB内存迭代过程反复重算会非常慢。建议预先算一次距离矩阵之后用索引更新或者改用K-means这种不需要完整距离矩阵的算法。如果场景数量级上万优先考虑K-means路线不要硬跑同步回代。5.3 我踩过的几个坑Beta分布拟合在光照数据非零的情况下才有意义。如果历史数据里包含大量零值时段整体拟合一个Beta分布会得到非常奇怪的参数。我踩过这个坑之后改为对非零时段单独拟合零值时段直接按0处理。另一个坑是K-means聚类结果的随机性。K-means初始质心是随机选择的同样的数据跑两次结果可能不同。解决办法是固定随机种子或者使用多次运行取最优结果。kmeans函数可以通过Replicates, 10参数让Matlab自动跑10次并返回最优聚类结果强烈建议加上。还有一点容易被忽略场景生成时坐标和单位的归一化。风速用m/s、出力用MW、光照用百分比直接混在一起算距离的话量级差异会导致距离矩阵几乎由量级最大的变量主导。需要先把所有变量归一化到同一尺度比如统一除以各自最大值或标准差。否则削减结果会偏向于照顾量级大的电源另一个电源的场景特征可能被悄悄丢掉。5.4 排查问题的基本流程场景生成结果异常时先不要急着调算法参数按顺序排查这三件事。第一检查输入数据。历史风速和光照数据是否包含异常值、缺测值时间序列是否有明显跳变。数据质量差是场景生成结果离谱的最常见原因。第二检查概率分布拟合结果。用plot把拟合曲线和实际数据直方图叠在一起看直观判断拟合得好不好。仅靠拟合出的数值参数很难发现分布形态不对。第三检查随机数生成过程。固定随机种子重复生成两次初始场景对比结果是否一致。如果不一致说明代码某处存在未固定的随机源比如K-means初始质心或抽样函数。固定rng后结果仍然不一致就要检查是否有循环内重新设置随机种子的操作。调试期间我习惯把中间变量逐一保存成mat文件每一步结束后都画图检查。场景生成这类随机性很强的代码一步错往往要到最后才发现结果异常提前可视化能省下大量排查时间。结尾分享场景生成与削减这套流程看起来是简单的统计抽样加聚类真正跑通后才发现细节远比想象的多。分布参数拟合不好后续场景质量全线崩塌削减算法选错精度和速度双双受损连距离计算里量纲归一化这种小事都可能让削减结果偷偷跑偏。我在多次项目迭代中的一个体会是场景削减不是一个一步到位的步骤而是一个需要反复试算、不断与下游优化结果对照反馈的过程留出足够的调参时间比追求一次跑通更重要。如果你也在用Matlab做新能源不确定性建模建议先拿一套自己熟悉的历史数据把生成、削减、评价这三步完整跑通再逐步换成新数据和新场景。希望这篇分享能帮你少走几步弯路。