计量、需求分析的交互方式，我在MATLAB里实现时，算了这样一笔账：

给定一天内的风速均值函数，按小时从夜间低风速爬升到午后高风速再回落。如果硬生生从原始序列里直接取96个时刻做场景，低频段占比集中的特征会被切散，很难忠实地保留风速变化规律。改用概率距离消减后，优先保留的是变化幅度最有代表性的时段——典型的风速回升段会被重点保留，夜间平稳段则可能被合并掉一两个相近的小时。这在潮流计算里差别很大：前者算出来的馈线电压波动特征更像实际运行曲线，后者可能把最关键的爬坡时段磨平了。

再比如并网系统的容量配置研究里，如果场景削减后把少数极端低辐照、高负荷的组合丢掉，最终算出的储能配置可能偏小，一旦碰上连阴天就撑不住。这也是为什么同样一批数据，有的论文用随机抽样法做场景削减，结果储能的充放电策略会明显偏离典型场景下的最优策略，就是因为随机抽样对概率密度的保留不如概率距离法精准。

2.4 方法定位：它不是万能的

有一点需要说透：概率距离快速削减法擅长的是“从已有的大量场景里挑出最有代表性的一个子集”，它并不擅长“从无到有生成全新场景”。所以实际工程里通常是把两类方法串联起来用——先用ARMA、Copula或者马尔可夫链等方法生成成百上千个备选场景，再用概率距离快速削减法降维得到几十个典型场景。先扩充再浓缩，这符合“先广撒网后聚焦”的思路。如果数据源本身就比较单薄，比如只有一两周的量测数据，直接套用这套流程生成全年场景，那结果多半不够看，因为备选场景的总量不够，削减出来的东西说服力有限。

3. 代码实现拆解：从关键函数到调用链路

3.1 主流程设计思路

我实现这套MATLAB代码时，参考了一个比较经典的程序结构：主脚本负责数据读取、参数设置和结果可视化，核心削减部分封装成独立函数。这样做的好处是，复现或调整研究场景时，不需要动主逻辑，只需要改输入数据和削减次数，跑出来的结果立即可视化对比。整体流程基本是：

1. 读取或生成原始风光场景序列，数据格式采用矩阵，每一行代表一个时间断面，每一列代表一个备选场景；
2. 设置目标削减数目（比如把500个场景消减成10个）；
3. 调用概率距离快速削减函数，内部迭代计算场景之间的概率距离、逐步剔除最不重要的场景；
4. 输出保留场景的编号、对应的发生概率（权重）以及削减后的场景数据；
5. 绘制削减前后的场景对比图、削减后的典型场景曲线和概率权重柱状图。

提示：很多人刚开始容易把行列方向搞反。MATLAB的惯性习惯是“列向量优先”，如果原始数据里每一列是一个时间序列，需要先转置成“每行一个断面、每列一个场景”的格式，再喂给削减函数，否则算出来的“场景距离”语义完全不对。

3.2 核心削减函数的实现逻辑

概率距离快速削减法的核心函数，我用标准结构实现如下（关键部分保留示意）：

function [keptScenarios, keptProb, reducedData] = probDistReduction(scenarios, probs, targetNum) % scenarios: 矩阵, 每一列是一个备选场景 % probs: 每个场景的初始概率 % targetNum: 期望保留的场景数 numScen = size(scenarios, 2); origData = scenarios; currProbs = probs; removedFlag = false(1, numScen); while (numScen - sum(removedFlag)) > targetNum % 计算当前保留场景之间的概率距离（采用2-范数） activeIdx = find(~removedFlag); numActive = length(activeIdx); distMat = zeros(numActive); for i = 1:numActive for j = 1:numActive if i ~= j diffVec = origData(:, activeIdx(i)) - origData(:, activeIdx(j)); distMat(i, j) = norm(diffVec, 2); else distMat(i, j) = Inf; end end end % 寻找使概率距离损失最小的场景进行剔除 minLoss = Inf; idxToRemove = -1; for i = 1:numActive loss_i = 0; for j = 1:numActive if i ~= j loss_i = loss_i + currProbs(activeIdx(j)) * distMat(i, j); end end if loss_i < minLoss minLoss = loss_i; idxToRemove = i; end end % 剔除场景，并把被剔除场景的概率累加到与之距离最近的保留场景上 removedIdx = activeIdx(idxToRemove); removedFlag(removedIdx) = true; [~, nearestIdx] = min(distMat(idxToRemove, :)); currProbs(activeIdx(nearestIdx)) = currProbs(activeIdx(nearestIdx)) + currProbs(removedIdx); end keptScenarios = activeIdx(find(~removedFlag(activeIdx))); keptProb = currProbs(keptScenarios); reducedData = origData(:, keptScenarios); end

这段代码里最关键的几个点：

• 概率距离的计算：用的是两个场景向量之间的欧氏距离，本质上是衡量“这两个场景在时序曲线形状上差多大”。如果想进一步考虑不同时刻的重要性差异，可以在这个距离公式上叠加一个权重向量，比如对负荷尖峰时段加大权重，算出来的削减结果会更贴合调度关心的时段。
• 剔除标准的推导：对候选场景 i 计算剔除后“所有其他场景到它”的加权距离总和，选总和最小的那个剔除。这个总和越小，说明被剔除场景和其他场景越相似，去掉它对整体概率分布的影响越小。
• 概率再分配：剔除场景后，它的概率要累加到离它最近的保留场景上。这个动作相当关键，它就是“概率距离削减”这个名字里“概率”二字的落点——削减之后不是简单扔掉一些曲线，而是把它们的代表性折到幸存的场景里，保证总概率质量守恒。

注意：这个函数里我用双重循环求距离矩阵，在场景数量达到千级以上时会有性能压力。实际工程中可以先预分配distMat矩阵，再考虑用向量化或并行计算优化。但在教学和中小规模研究中，双重循环的写法更直观，也更方便读者把算法逻辑“翻译”成论文里的公式。

3.3 主脚本与可视化验证

主脚本部分，我按“生成原始场景 -> 削减 -> 对比展示”的顺序组织：

clear; clc; rng(42); % 生成模拟场景: 以典型日风速曲线为基准叠加上噪声 hours = 0:24; baseWind = 4 + 3 * sin(pi * (hours - 6) / 12); numScen = 500; windData = zeros(length(hours), numScen); probs = ones(numScen, 1) / numScen; % 初始等概率 for k = 1:numScen windData(:, k) = baseWind' + 1.2 * randn(length(hours), 1); end % 调用削减函数, 保留10个典型场景 targetKeep = 10; [keptIdx, keptProb, reducedWind] = probDistReduction(windData, probs, targetKeep); % 绘制对比图 figure; subplot(1,2,1); plot(hours, windData, 'Color', [0.8 0.8 0.8]); hold on; plot(hours, reducedWind, 'LineWidth', 1.5); xlabel('时刻/h'); ylabel('风速/(m/s)'); title('削减前后场景对比'); subplot(1,2,2); bar(keptIdx, keptProb); xlabel('保留场景编号'); ylabel('概率权重'); title('典型场景概率权重');

这里有一个实操时很有用的细节：rng(42) 固定随机种子。为什么强调这一点？因为场景生成本身带有随机性，如果你不固定种子，每次运行主脚本得到的原始场景都不一样，削减结果自然也不同。在调试阶段这会非常痛苦——你可能改了一行代码，以为结果变化是代码引起的，实际上只是随机数变了。固定随机种子能确保实验的可复现性，这也是很多论文对算法结果做敏感性分析的基础。

另外，绘制削减前后的场景对比图时，原始场景用浅灰色细线绘制，削减后的典型场景用深色实线加粗显示，两者叠加后能直观看出削减后的场景是否很好地“笼罩”在原始场景的变化范围内。如果削减后的曲线形态明显偏离原始场景包络线，说明削减数目设置太少或者备选场景生成质量不足，需要回头检查。

4. 风景场景生成与削减的完整实操流程

4.1 从原始数据到备选场景：一个通用流程

方法和框架说清楚了，还得聊聊具体的实操闭环。我习惯把整个流程分成四大步：数据准备、场景生成、场景削减、结果应用。

数据准备阶段，需要确定研究区域和时间粒度。时间粒度一般取1小时，一天就是24个断面；如果做日内滚动优化，可能用到15分钟粒度，一天就是96个断面。粒度越细，场景向量的维度越高，削减时的计算量也越大，但保留的时序细节越丰富。实际上做全年8760小时的场景削减，计算压力会非常大，所以更常见的做法是按“典型日”来做——把全年数据按季节或月份聚类成若干类典型日，再在每一类内部做场景生成和削减。每个典型日只需处理24或48个时间断面的场景，计算量小很多。

场景生成阶段，根据可用数据情况选择方法：

• 如果只有历史风速/辐照度的统计分布参数，最省事的做法是用蒙特卡洛抽样，从边缘分布里抽数，再通过相关性矩阵生成风速-辐照度的联合场景。
• 如果有历史序列的时序相关性，可以用ARMA模型或一阶马尔可夫链来生成，能更好地保留“今天风速高明天风速也高”这种持续性特征。
• 如果有风电场和光伏电站的地理相关性数据，用Copula函数把两个边际分布“粘”起来，效果更好。

这一步输出的是一大批备选场景，比如500到2000个，数量多多益善，因为削减算法本质上是在一个高维空间里做代表点选取，备选集越稠密，选出来的代表场景越有说服力。

场景削减阶段，调用我们第三部分写好的核心函数，把1000个场景消减成10到20个。削减后的场景集合就构成了后续仿真优化的输入。

结果应用阶段，最常见的是把削减后的典型场景及其概率权重代入日前调度模型，目标函数是“期望运行成本最小”——每个典型场景对应一组出力曲线和概率，把各个场景下的成本按概率加权求和，得到的就是考虑了不确定性的期望成本。相比只取确定性预测场景（单场景）的做法，这种多场景期望模型明显更贴合真实运行需求。

4.2 风光联合场景削减的特殊细节

当题目里的“风光场景”是风速和辐照度联合场景时，需要注意一个实际操作细节——量纲。

风速的量纲是m/s，典型值在0到15之间；辐照度的量纲是W/m²，典型值在0到1000之间。如果把风速和辐照度直接拼成一个向量去算欧氏距离，辐照度的数值会完全主导距离计算，风速的差异几乎可以忽略。这样削减出来的场景，很可能风速特征保留得不好，而辐照度特征保留得“过好”，达不到两种能源特性均衡保留的目标。

解决办法是在拼装场景向量之前做归一化或标准化处理。我常用的做法是 min-max 归一化，把风速和辐照度分别压到 [0,1] 区间再拼接，等效于给两种变量相同的初始权重。如果研究更关注某一种资源，也可以故意把该变量的缩放比例调大，实现加权效果。

% 风光联合场景拼装前的归一化示例 windNorm = (windData - min(windData)) ./ (max(windData) - min(windData)); irrNorm = (irrData - min(irrData)) ./ (max(irrData) - min(irrData)); jointData = [windNorm; irrNorm];

还有一点值得注意：风速和辐照度在很多地区存在天然的负相关——白天辐照度高的时候，风速可能偏低或偏高，需要用Copula函数保留这种相关结构。如果直接在拼装时独立抽样，就会把相关性抹掉，生成的联合场景在物理上就不自洽。

4.3 削减数目的选择技巧

目标削减数目不是拍脑袋定的。它直接影响后续优化问题的规模和精度。

从优化角度看，保留的场景数越多，期望成本模型的精度越高，但决策变量的规模和求解时间也越大。对典型日前调度模型来说，如果保留场景数从5个增加到20个，模型里的变量数量会成倍增加，混合整数规划求解时间可能从几十秒膨胀到几十分钟。所以实际中需要做一个权衡：先用10个场景跑一遍，再逐步增加到20、30个，观察目标函数值的变化。如果目标函数值基本稳定（比如变化幅度小于1%），就没必要保留更多场景了。

从统计角度看，有一种经验法则是保留场景数大致取备选场景数的平方根。比如从400个备选场景里削减，保留20个左右比较合适。这个经验法则来自聚类分析中的类数选择直觉，不一定严格最优，但作为初始值非常好用。之后再结合目标函数值做敏感性分析，找到一个“性价比最高”的削减数目。

我在实操中常用这样一组对比：分别保留5/10/20/30个场景，每个配置跑同一组优化问题，记录求解时间和目标函数值，做成一张二维对比表。目标函数从5到10个场景通常变化明显，从10到20个场景变化收窄，从20到30个场景几乎不动，这时候就锁定10到20个作为最终配置。

5. 常见问题与排查技巧

5.1 削减结果出现负值或超界

我自己上手时遇到的第一个问题，就是削减后的场景曲线出现负值或超过历史最大值的点。排查下来根源在场景生成阶段：如果直接用蒙特卡洛从正态分布抽样，抽到的风速值可能是负的，或者光伏辐照度超过1000W/m²。这在物理上不合理。

解决办法有两个层次。

第一层是生成阶段就限制：抽样后把小于零的值截断为0，把大于辐照度上限的值截断为上限。这虽然简单粗暴，但对后续削减影响不大，因为截断后的分布依然保留了主要的概率质量。

第二层是削减后处理：如果削减出来的典型场景里还有非物理值，可以再做一次边界钳制。需要注意的是，钳制操作会轻微改变场景间的距离结构，但对于工程研究来说影响完全可以接受。

5.2 削减后概率权重总和不为1

这是概率再分配实现错误时的典型症状。正常情况下，初始概率和为1，每次剔除场景时把它的概率累加到最近场景上，总和应该恒为1。如果发现削减后概率总和偏离1，大概率是代码里“最近邻居”的查找逻辑写错了，比如找最近邻居时把被剔除场景自己也包含了进去，导致概率累加时加到了0上。

调试手段很简单：在每次剔除场景后打印所有保留场景的概率总和，观察是否恒为1。如果某一轮开始偏离，说明问题出现在那一轮的最近邻查找或概率更新逻辑上。这类问题用断点调试比空想快得多。

5.3 削减结果不够“典型”

这个问题的表现是：削减后的场景曲线五花八门，很难看出共性，或者每个场景都像是加了很多噪声的毛刺曲线，而不是光滑的典型日曲线。出现这种情况，多半是备选场景生成时的噪声太大，削减算法被高维噪声带偏了。

解决方案是降低噪声方差，或者在生成备选场景之前先做一次平滑滤波。实际操作中，我建议把噪声强度的选择与后续研究需求绑定：如果是研究极端场景对系统的影响，噪声可以适当大一些，保留尾部风险；如果是研究典型日运行方式，噪声应该小一些，让削减后的场景更平滑、更接近理想典型曲线。

5.4 性能瓶颈与仿真卡顿

当备选场景数量超过2000、时间断面超过48时，双重循环计算距离矩阵的耗时可能达到几分钟量级，这时候主脚本调用削减函数会明显卡顿。

几个优化思路：

• 先按时间维度降采样，比如从15分钟粒度改成30分钟或1小时，维度降一半，计算量减半；
• 在MATLAB里预先分配distMat，避免循环内动态扩容；
• 使用向量化计算距离矩阵，比如利用矩阵乘法和广播特性替代双重循环。

% 向量化计算欧氏距离矩阵的示例 numActive = size(activeData, 2); a2 = sum(activeData.^2, 1); distMat2 = sqrt(max(bsxfun(@plus, a2', a2) - 2 * activeData' * activeData, 0));

这段代码里的核心技巧是：任意两列向量之间的欧氏距离平方，等于两列向量自身模长的和减去两倍的点积。使用bsxfun或隐式扩展，可以一次性得到整个距离矩阵，避免显式循环。我实测下来，500×24的数据量，向量化版本比双重循环快上几十倍，场景规模越大越明显。

5.5 削减结果与预期场景差异大

最后一种常见情况是，削减出来的场景时间分布与直观经验不符。比如冬季典型日的风速场景削减结果，居然没有保留明显的冬季大风过程，反而凑了一堆过渡季节曲线。这种时候首先要怀疑的是备选场景生成的季节性偏差，而不是削减算法本身。检查生成备选场景时使用的分布参数是不是按季节分别统计的；如果全年混在一起生成，季节性规律就会被平均掉，削减自然也找不回冬季特征。

解决办法是按季节或月份分别生成备选场景，再分别削减，最后按各季节的天数占比加权组合成全年典型场景集。

6. 一个小技巧与个人体会

最后分享一个我在多次实操中总结出来的小技巧：不要只用削减后的场景做一次优化就收工，一定要追加一个回代校验。

回代校验的做法是：把削减后的每个典型场景当作确定性输入，算出一组调度决策；然后把所有备选场景（包括被削减掉的那些）代入这批决策，计算它们各自的运行成本和约束越限情况，最后按原始概率加权统计。如果加权越限率接近零、期望成本与基于削减场景计算的期望成本相差不大，说明削减后的场景集很好地代表了原始场景集；如果差异明显，说明削减数目太少或者削减算法丢失了重要尾部风险，需要调整。

这个方法从理论上对应了“场景削减误差评估”的思路——削减本身只是手段，核心目的是让削减后的场景集合在后续优化中产生和原始场景集足够接近的决策质量。多花五分钟做一次回代校验，比反复调算法参数要靠谱得多。

我在用这套流程处理一个模拟项目X的风光容量配置研究时，最开始削减到10个场景，回代校验发现极端高风速场景下的出力越限率达到4%，后来把保留场景数提高到20个，越限率降到1%以内，同时求解时间增加了大约一半。这个代价换取的风险控制收益，我认为非常值得。

概率距离快速削减法作为一个场景降维工具，真正吃透之后，你会发现它不只是论文里的一个公式，而是可以反复使用、稳定输出的一把好用的尺子。希望这篇文章能帮你在自己的研究或项目中，把风光场景生成与削减这件事做得又快又好。