企业官网建设流程全解析

1. 为什么风光场景生成非得做"削减"这一刀

做新能源电力系统规划或调度优化的朋友，对"场景法"应该都不陌生。风电和光伏出力受天气影响，具有强烈的随机性和波动性，我们没办法用一个确定的曲线去描述"明天风多大、太阳多强"，只能通过历史数据或概率模型，生成大量可能出现的出力场景。理论上场景数量越多，对不确定性刻画越精确，但随之而来的问题是：计算量呈指数级膨胀。

举个例子，一台风机的出力场景如果生成1000个，那么一个含10台风机、5个光伏电站的系统，组合起来的场景总数就是天文数字。调度模型里每个场景都对应一组约束条件，求解器面对成千上万的场景约束，轻则求解时间从分钟级变成小时级，重则直接内存爆炸无法求解。所以业界普遍的做法是：先生成大量场景，再用某种方法削减到少量有代表性的场景，用这几十个场景去逼近原始上千场景的统计特性。

这个思路听起来简单，但削减方法选得不好，结果会非常难看。我见过有人直接用K-means聚类做场景削减，聚类中心倒是选出来了，但聚类算法本质上是把场景当成"点"来划分，它压根不考虑场景之间的转移概率和时序相关性，削减后的典型场景集在概率分布上往往失真，后续优化结果偏差很大。

而概率距离快速削减法（也叫基于Kantorovich距离的场景削减法，简称SD场景削减）是目前学术界和工业界都认可的主流方法。它的核心逻辑不是"找聚类中心"，而是通过计算场景之间的概率距离矩阵，逐步剔除对整体概率分布影响最小的场景，并把被剔除场景的概率累加到距离它最近的保留场景上。每一步削减都保证原始场景集与削减后场景集的概率距离增量最小，因此能在给定场景数条件下，最大程度保留原始分布特性。

MATLAB实现这套方法，代码量不大，但有几个坑非常隐蔽。本文我用一个完整的风电+光伏联合场景生成与削减案例，把原理、代码、调参经验和验证方法一次讲透。适合正在做电力系统随机优化、微电网容量配置、风电光伏消纳分析的研究生和工程师参考。

2. 风电光伏联合场景生成：先弄懂我们削减的对象是什么

2.1 风力和光伏出力的数学模型怎么建

要削减场景，第一步是得有场景。场景生成的方法很多，常见的有：基于历史数据的直接采样、基于ARIMA时间序列模型、基于马尔可夫链、基于Copula函数、基于概率密度函数的蒙特卡洛采样等。不同的方法适合不同的数据条件。

在实际项目中，我通常采用"概率分布模型+时序相关性还原"的方式。风电出力近似服从Weibull分布，光伏出力近似服从Beta分布，但直接按这两个分布独立采样出来的场景毫无时序逻辑——风不会上一秒满发下一秒突然归零，光伏也不会中午出力呈白噪声式跳动。所以完整的场景生成链路是：

1. 对风速历史数据拟合Weibull分布，通过风速-功率转换曲线得到风电出力序列。
2. 对光照强度历史数据拟合Beta分布，通过光照-功率转换模型得到光伏出力序列。
3. 引入自相关系数和互相关性，让生成的场景在时间维度上平滑、在风-光之间具有合理的相关性。
4. 用蒙特卡洛采样生成N个初始场景，每个场景包含24个时段的双维出力（风电+光伏）。

需要说明的是，场景削减算法本身并不关心场景是怎么生成的，它只负责在给定初始场景集和概率权重（通常假设等概率）的前提下，找到最优的保留子集。所以哪怕你用别的方法生成场景，只要格式是"场景数×时段数×变量数"，后面的削减代码可以直接复用。

2.2 初始场景规模如何确定

初始场景规模N的选取直接影响生成精度和削减效果。太少了，削减后的典型场景缺乏代表性；太多了，计算概率距离矩阵时的时间和内存开销会显著增加。

我的经验法则是：单变量场景至少1000个起步，风光联合场景建议2000到5000个。1000个以上场景经过削减到10-20个，统计特性保留效果较好；如果初始只有200个场景，削减到10个之后误差会明显偏大。

这里贴一段生成初始风电场景的MATLAB核心代码，采用Weibull分布采样并加入时序自相关处理：

% 参数设置 nScenarios = 2000; % 初始场景数 nPeriods = 24; % 调度时段数（小时） k = 2.1; % Weibull形状参数，根据历史风速拟合 c = 8.5; % Weibull尺度参数，m/s vRated = 12; % 额定风速 vCutIn = 3; % 切入风速 vCutOut = 25; % 切出风速 Prated = 1.5; % 额定功率 MW % 生成独立Weibull风速样本 windSpeed = wblrnd(c, k, nScenarios, nPeriods); % 加入时序平滑（简单一阶自回归，让相邻时段风速平滑过渡） rho = 0.85; % 自相关系数 for t = 2:nPeriods windSpeed(:, t) = rho * windSpeed(:, t-1) + ... (1-rho) * windSpeed(:, t) + randn(nScenarios, 1) * 0.3; windSpeed(windSpeed < 0) = 0; end % 风速→功率转换 windPower = zeros(nScenarios, nPeriods); for s = 1:nScenarios for t = 1:nPeriods v = windSpeed(s, t); if v < vCutIn || v > vCutOut windPower(s, t) = 0; elseif v >= vRated windPower(s, t) = Prated; else % 线性区间的近似功率曲线 windPower(s, t) = Prated * (v - vCutIn) / (vRated - vCutIn); end end end

光伏场景的生成逻辑类似，区别在于Beta分布拟合的是光照强度的归一化值，并且时序上要体现"白天出力、晚上归零"的强周期特征。联合场景就是把两个矩阵横向拼接，每个场景行向量维度变成48（24时段风电+24时段光伏）。

3. 概率距离快速削减法原理：每一步剔除都让分布损失最小

3.1 Kantorovich距离到底是什么

很多教程一上来就堆公式，把人看懵。我用一句话概括概率距离削减法的本质：它在反复问一个问题——"如果把某个场景删掉，把它的概率转移给离它最近的另一个场景，整体概率分布变化有多大？如果这个变化是所有候补删除项中最小的，那就删它。"

衡量"变化有多大"的标尺就是Kantorovich距离，也叫Wasserstein距离。两个场景集合之间的Kantorovich距离，可以理解为"把一个概率分布搬到另一个概率分布所需的最小代价"，这里的代价就是场景差异乘以转移概率。

对于离散场景集，Kantorovich距离的计算公式可以转化为一个最优运输问题。但在场景削减这个特定情境下，我们不需要求解完整的最优运输问题，因为削减过程是逐步进行的，每一步从一个场景集中删除一个场景，并把其概率累加到最近邻场景上。这种迭代贪婪策略得到的削减结果，在每一步都是局部最优的，实际应用中已经足够逼近全局最优。

3.2 算法的完整流程拆解

概率距离快速削减法的具体步骤如下：

1. 初始化：有N个场景，每个场景概率权重p_i = 1/N（也可根据历史频率赋权）。
2. 计算所有场景两两之间的欧氏距离d(i, j)，得到N×N对称距离矩阵。距离定义可以是逐时段绝对误差之和（L1范数），也可以是欧氏距离（L2范数），常用的是L2范数。
3. 对于每个场景i，找到与其距离最近的场景j（j ≠ i），记录距离d_min(i) = min d(i, j)。
4. 计算每个场景的"可删除代价"：cost(i) = p_i × d_min(i)。这个值的物理含义是：删除场景i后，概率损失（p_i）乘以它不得不转移给最近邻居的"距离代价"（d_min(i)）。
5. 选出cost最小的场景k，把它删除，并将其概率p_k累加到最近邻居j上：p_j = p_k + p_j。
6. 更新距离矩阵（删除第k行第k列），重复步骤3-5，直到保留的场景数量达到预设值M。

用伪代码表达就是：

while nScenarios > nTarget % 对每个场景找最近邻居距离 minDist = zeros(nScenarios, 1); minIdx = zeros(nScenarios, 1); for i = 1:nScenarios dist_i = distMatrix(i, :); dist_i(i) = inf; % 排除自身 [minDist(i), minIdx(i)] = min(dist_i); end % 计算删除代价 deleteCost = probWeights .* minDist; % 找到代价最小的场景 [~, delIdx] = min(deleteCost); % 概率转移给最近邻居 probWeights(minIdx(delIdx)) = probWeights(minIdx(delIdx)) + probWeights(delIdx); % 删除场景和对应行列 distMatrix(delIdx, :) = []; distMatrix(:, delIdx) = []; probWeights(delIdx) = []; nScenarios = nScenarios - 1; end

注意这里有个细节：距离矩阵每删除一个场景就要更新一次，因为删除后某些场景的最近邻居可能发生变化，必须重新计算最近邻距离。这也是该算法时间复杂度为O(N³)的原因——每轮删除都要重新扫描一遍距离矩阵。不过对于几千个场景来说，MATLAB跑起来也就几秒到几十秒，完全可以接受。

3.3 为什么不是一次删一批

有一种快速近似做法是预先计算所有场景的距离矩阵，然后一次性选出多个待删除对象，再统一做概率转移。这样做速度更快，但存在一个问题：第一批被判定为"可删"的场景中，可能会有两个场景互为最近邻居。如果同时删除它们，各自的概率都要转移到对方身上，就出现了"概率无处安放"的矛盾。

更合理的做法是每删除一个场景立即更新距离矩阵并重新计算最近邻，保证每一步操作都是基于最新场景集的。虽然计算量更大，但结果更可靠。对于超大规模场景（上万个），可以考虑用KD-Tree等近似最近邻搜索加速，但对于电力系统场景削减的典型规模，暴力计算完全够用。

4. MATLAB完整实现：从距离矩阵到削减主循环

4.1 场景切削函数封装

实际项目中我不会把主脚本堆成一坨，而是封装成一个通用函数，输入初始场景矩阵和期望保留数，输出削减后的场景、对应概率和削减过程中的误差记录。以下是完整的MATLAB函数代码：

function [reducedScenarios, reducedProb, distHistory] = ... scenarioReduction(scenarios, probs, targetCount) % 基于概率距离快速削减法的场景削减函数 % 输入： % scenarios - 初始场景矩阵，维度 [N, T]，N为场景数，T为时段数（含多变量拼接） % probs - 场景概率向量，维度 [N, 1]，缺省时默认等概率 % targetCount - 目标保留场景数 % 输出： % reducedScenarios - 削减后场景矩阵，维度 [targetCount, T] % reducedProb - 削减后场景概率向量，维度 [targetCount, 1] % distHistory - 每次删减后的Kantorovich距离记录，用于观察削减误差 N = size(scenarios, 1); if nargin < 2 || isempty(probs) probs = ones(N, 1) / N; end if nargin < 3 error('必须指定目标场景数'); end if targetCount >= N error('目标场景数必须小于初始场景数'); end % 计算初始距离矩阵（欧氏距离） distMatrix = squareform(pdist(scenarios, 'euclidean')); % 避免对角线为0，置为无穷大 distMatrix(1:N+1:end) = inf; distHistory = zeros(N - targetCount, 1); iter = 0; while size(distMatrix, 1) > targetCount iter = iter + 1; currentN = size(distMatrix, 1); % 找每个场景的最近邻居 [minDist, minIdx] = min(distMatrix, [], 2); % 计算每个场景被删除时造成的概率距离增量 deleteCost = probs(:) .* minDist(:); % 选择代价最小的场景删除 [~, delIdx] = min(deleteCost); % 概率累加到最近邻居 neighborIdx = minIdx(delIdx); probs(neighborIdx) = probs(neighborIdx) + probs(delIdx); % 记录当前的累积削减误差（用当前最小代价近似单步损失） distHistory(iter) = deleteCost(delIdx); % 从距离矩阵和概率向量中删除该场景 distMatrix(delIdx, :) = []; distMatrix(:, delIdx) = []; probs(delIdx) = []; end % 最后提取保留的场景：需要知道原场景索引 % 由于不断删行，距离矩阵的行序和场景矩阵的行序在删除后不一致， % 这里需要用逻辑索引追踪原位置，简化方式是通过在削减过程中同步删除场景矩阵的行来实现。 % 由于函数输入未返回索引，这里采用同步删除方式重建场景矩阵

上面代码有个需要注意的地方：我在函数开头提到需要追踪原场景索引。实际上最稳妥的做法是在删除距离矩阵行列时，同步删除场景矩阵中的对应行，否则最后不知道哪些原始场景被保留了。修正后的核心循环如下：

reducedScenarios = scenarios; idxMap = (1:N)'; while size(reducedScenarios, 1) > targetCount currentN = size(reducedScenarios, 1); distMat = squareform(pdist(reducedScenarios, 'euclidean')); distMat(1:currentN+1:end) = inf; [minDist, minIdx] = min(distMat, [], 2); deleteCost = probs(:) .* minDist(:); [~, delIdx] = min(deleteCost); neighborIdx = minIdx(delIdx); probs(neighborIdx) = probs(neighborIdx) + probs(delIdx); % 同步删除场景行和概率 reducedScenarios(delIdx, :) = []; probs(delIdx) = []; end

每轮循环重新调用pdist计算距离矩阵，比维护一个动态更新的距离矩阵更直观。虽然理论上多次调用pdist会增加计算量，但对几千个场景、48维的场景矩阵来说，MATLAB性能完全扛得住，代码可读性却高了很多。

4.2 主脚本：生成场景→削减→可视化

主脚本的逻辑很清晰，生成2000个风光联合场景，削减到10个，然后画出削减前后对比图。代码如下：

clear; clc; close all; rng(2024); % 固定随机种子，保证结果可复现 % 生成初始场景（以2.2节的windPower为基础，增加光伏部分） % 此处直接载入预生成的场景矩阵，实际使用时替换为自己的场景生成代码 % scenarios: [2000, 48]，前24列风电，后24列光伏 load('initialScenarios.mat', 'scenarios'); N = size(scenarios, 1); probs = ones(N, 1) / N; % 设置削减目标数 M = 10; % 调用削减函数 [redScen, redProb, histErr] = scenarioReduction(scenarios, probs, M); % 输出结果 fprintf('初始场景数: %d，削减后场景数: %d\n', N, M); fprintf('削减后概率分布: \n'); disp(redProb(:)'); % 可视化：削减前后某几个场景的对比 figure; subplot(2,1,1); plot(1:24, scenarios(1:50, 1:24)', 'Color', [0.7 0.7 0.7]); hold on; plot(1:24, redScen(:, 1:24)', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('风电出力/MW'); title('风电场景：削减前后对比'); legend('初始场景', '削减后典型场景'); grid on; subplot(2,1,2); plot(1:24, scenarios(1:50, 25:48)', 'Color', [0.9 0.7 0.7]); hold on; plot(1:24, redScen(:, 25:48)', 'LineWidth', 1.5); xlabel('时段/h'); ylabel('光伏出力/MW'); title('光伏场景：削减前后对比'); legend('初始场景', '削减后典型场景'); grid on;

运行后，削减后的10个典型场景会在2000个原始场景的"云图"中清晰勾勒出主要分布模式。我截取一张典型输出效果描述：风电场景削减后的曲线基本覆盖了从低出力到高出力的各个典型水平，光伏场景则呈现出明显的午间高峰形态，概率最大的场景是整个分布中"最普通"的那种——白天光伏高发、晚间风电接力。

5. 削减效果怎么评估：三种Kantorovich距离验证方法

5.1 削减前后概率距离的直接计算

最直观的评估方法是计算削减前场景集与削减后场景集之间的Kantorovich距离。这个距离越小，说明削减后的场景集对原始分布的还原度越高。MATLAB中可以通过调用Wasserstein距离的相关实现，或直接按定义用线性规划求解最优运输问题。对于中等规模问题，用transport函数或者CVX求解即可。

但实际中我们更关心的是：随着目标场景数M从1增加到50，削减误差怎样变化？把这个变化曲线画出来，就能直观看到"用多少个场景能在精度和计算量之间取得平衡"。我的经验是：对于风力发电场景，M=10时距离通常能降到初始距离的5%以内，继续增加M，边际收益开始大幅递减。

M_values = 1:10:50; errorValues = zeros(size(M_values)); for i = 1:length(M_values) [~, ~, ~] = deal([]); [redScen, redProb, histErr] = scenarioReduction(scenarios, probs, M_values(i)); errorValues(i) = sum(histErr); % 累计削减误差 end figure; plot(M_values, errorValues, 'o-', 'LineWidth', 1.5); xlabel('保留场景数'); ylabel('累计削减误差'); title('削减误差随目标场景数的变化'); grid on;

5.2 与K-means和K-medoids聚类的对比

很多场景削减文献会对比概率距离法与聚类法的差异。我建议你在自己的算例上也做一次对比，因为在某些数据结构下K-medoids和概率距离法的结果很接近，但概率距离法在概率分配上天然更合理。

K-means聚类削减的问题在于：聚类得到的簇中心是一个虚拟场景，它可能不代表任何真实存在的出力模式，而且它无法保留场景的出现概率概念。K-medoids虽然选取的是真实场景作为簇中心，但每个簇内场景被视为等权重，忽略了簇内场景的密度信息。

概率距离法则通过概率转移机制，让概率密度高的区域保留更多场景，概率密度低的区域被合并。直观理解：如果某个区域有400个相似场景，削减后这个区域可能保留4个场景，每个概率约0.1；如果另一区域只有50个场景，削减后可能只保留1个场景。这个特性对后续优化计算特别友好，因为概率权重可以被直接用作随机规划中的场景概率。

5.3 进一步检验：用削减后场景求解优化问题对比

这是最实际也最有说服力的验证方式。假设你正在做微电网容量配置，先用2000个场景求解一版结果，再用削减后的10个场景求解，对比两类结果的目标函数值偏差。如果偏差在可接受范围内（比如5%以内），说明削减方案有效，可以放心在更大规模求解中使用。

具体操作上，你需要一个两阶段随机规划模型（如投资决策+运行模拟），场景不确定性体现在风电光伏出力上。用削减后场景求解得到的投资方案，再拿到全场景下做运行模拟回代检验，计算期望运行成本。这个回代过程基本涵盖了随机优化中最关键的一步——场景削减与模型求解的闭环验证。

6. 实操中的关键坑与调参经验

6.1 归一化处理：风电和光伏的量纲差异是隐形杀手

这是我最想强调的一点。风电出力单位是MW，数值可能在0到1500之间；光伏出力如果做归一化处理，数值在0到1之间。如果你把这两个量直接拼接到一个48维向量里计算欧氏距离，距离计算会被风电的数值尺度主导，光伏的差异在距离中几乎不起作用。最终削减出来的场景，风电部分很精确，光伏部分则非常粗糙。

正确做法是先归一化。把风电和光伏都标幺到各自的装机容量下（即除以额定功率），然后再拼接计算距离。这样每个维度的量纲一致，距离度量才能均衡反映两类电源的差异。削减完成后，如果需要还原到实际功率，再乘以各自的装机容量即可。

% 归一化后再削减 windNormalized = windPower / windCapacity; pvNormalized = pvPower / pvCapacity; scenariosNormalized = [windNormalized, pvNormalized]; % 削减... % 还原 windReduced = reducedScenarios(:, 1:24) * windCapacity; pvReduced = reducedScenarios(:, 25:48) * pvCapacity;

6.2 距离度量选择：L1还是L2

默认的欧氏距离（L2）对异常场景更敏感，会让削减结果更偏向于方差较大的区域；曼哈顿距离（L1）则更稳健，在一些风光场景数据中表现更好。我建议两种都跑一遍，对比削减后的累积距离误差。如果数据中有明显的离群场景，L1通常更稳。

另外，有一种进阶做法是在距离计算中引入时段权重。比如在调度问题中，峰时段的出力偏差对系统运行的影响远大于谷时段，可以给峰时段赋予更高权重。具体做法是把原始场景向量乘上一个权重向量，再计算欧氏距离。这个技巧在工程实践中很实用，能让削减结果更贴合约束的关键区域。

6.3 目标场景数M怎么定

M的取值没有绝对标准，但有几个参考维度：

• 随机规划求解能力：如果用的是Cplex/Gurobi求解混合整数规划，M在10-20个时求解速度通常可接受；如果模型简单（如线性规划），M可以放到50个。
• 分布还原精度需求：用5.1节的误差曲线，观察误差下降的拐点，取拐点处的M值。
• 对比文献常用值：多数期刊论文中M取10、20、50三档做敏感性分析。

我个人的习惯是：先取M=20作为基准，然后做M=10和M=50的敏感性分析，把三条结果曲线放在同一张图里，审稿人或评审专家看到这种分析会觉得你考虑得很周全。

6.4 运行效率优化技巧

当初始场景数达到5000以上时，每次循环都重算距离矩阵会变得有些吃力。两个优化方向：

第一，用向量化操作替代循环。MATLAB的pdist本身是C语言实现的，效率已经很高，关键是尽量避免在循环里做额外的矩阵复制。第二，可以在每轮削减前先检查最近邻矩阵是否可以增量更新，只有当被删除场景的邻居关系受影响时才重算对应行。不过这个优化实现复杂度较高，一般场景规模下收益不大。

如果场景数量过万，建议先用K-means粗削减到500个，再用概率距离法精确削减到目标数。这种两级削减策略兼顾速度和精度，是工程中很实用的做法。

7. 案例：把风光场景削减用到微电网容量规划中

7.1 问题设定

我做一个简化版的案例演示。假设一个独立微电网，含风力发电、光伏发电和储能系统，需要确定风电、光伏和储能的配置容量，目标是在满足负荷需求的前提下，年综合成本最小。风速和光照的不确定性用场景描述。

初始生成2000个风光联合场景，用本文方法削减到20个典型场景。

7.2 求解结果对比

用削减后的20个场景求解容量配置，再将配置结果放回2000个全场景中做运行模拟，验证实际期望成本。结果是这样的（数据仅供演示，不同算例会浮动）：

从结果看，概率距离法削减后求解的决策变量与全场景求解结果非常接近，回代成本偏差仅0.36%，而K-means削减后求解的偏差达到4.8%。这个案例很直观地说明了削减方法选择对优化决策的真实影响。

7.3 从案例中得出的经验

第一，场景削减的价值不仅在于压缩计算量，更在于保留概率分布的关键特征。选对方法，20个场景能扛住2000个场景的统计精度；选错方法，50个场景也照样失真。

第二，评估削减质量不能只看削减前后的概率距离，必须结合下游具体问题（如容量配置结果）做闭环检验。有时候概率距离根小，但关键尾部分布被削没了，导致规划方案在极端场景下失去鲁棒性。

第三，风光联合削减不是简单的"分别削减再把结果拼起来"。联合削减能保留风-光之间的时序相关性和互补性，这对含高比例新能源的微电网规划至关重要。如果你分别对风电和光伏做削减，得到的典型场景组合可能根本不存在于原始联合分布中。

8. MATLAB实现中的代码细节与扩展思路

8.1 函数化封装的设计心得

我在写scenarioReduction函数时，刻意把输入输出设计成纯矩阵操作，不依赖全局变量。这样好处很明显：可以直接放进parfor循环做批量实验，比如同时测试不同的M值、不同的距离度量、不同的初始场景集，互不干扰。

推荐你把函数文件保存为scenarioReduction.m，放在工作目录下，然后写一个runExperiments.m主脚本，集中管理参数设置、调用和结果输出。我所有的MATLAB项目都遵循这个习惯：数据生成、核心算法、实验分析三个层次分离，改一个模块不影响其他模块。

8.2 可视化场景质量的小技巧

除了简单的曲线图，我还会画一种"削减效果热力图"：横轴是24小时时段，纵轴是不同的削减后场景，颜色表示出力水平。这种图可以快速看出削减后的场景是否覆盖了全时段的各种出力水平。另外，把削减后场景的概率权重画成一个饼图或条形图，能直观看到概率分布是否过于集中——如果出现某个场景权重超过0.3的情况，说明削减后代表性可能不足，建议增加场景数。

8.3 扩展到多变量和时序耦合场景

上面的例子只有风、光两个变量，实际中可能还需要考虑负荷的不确定性、电价的不确定性等。实现上只需把更多变量的时序数据拼接到场景矩阵中，维度变为[T × 变量数]。归一化时注意每个变量单独处理，避免某一变量主导距离计算。

对于需要考虑季节差异的场景集（比如四季各有不同特性），可以分季节分别生成和削减，再合并成包含季节标签的完整场景集。在随机规划中，这种做法还能同时表达"季节概率"和"场景条件概率"的层次结构。

8.4 与Matlab工具箱的配合

如果你的MATLAB安装了Optimization Toolbox，可以把削减后的场景直接接入随机规划求解框架。场景概率向量对应随机规划中的场景权重，削减后的场景矩阵作为不确定参数输入到约束中。另外，Statistics and Machine Learning Toolbox中的pdist函数是计算距离矩阵的基础工具；如果你要做K-means对比实验，kmeans函数也是现成的。

整套代码不需要任何额外的第三方工具箱，基础的MATLAB环境加Statistics工具箱就完全够用。这点对很多用学校正版License或者教育版的朋友很友好，不用担心依赖问题。

9. 写在最后：场景削减的实际项目经验小结

概率距离快速削减法是我在多个新能源规划项目中验证过的最可靠场景削减方案。它的优势不是某个特定算例的巧合，而是有严格的最优运输理论支撑——每一步削减都保证了概率距离损失最小，整个削减过程可解释、可验证、可复现。

分享几条我在实际项目中总结的经验：

1. 初始场景质量决定削减效果上限。生成场景时务必保证概率分布拟合的准确性，时序相关性处理不到位，削减再精确也没用。
2. 归一化是所有场景削减的第一步，不要跳过去。不同量纲的变量拼在一起，距离计算就失真了。
3. 削减后的验证必须闭环。算一下削减前后场景的均值、方差、相关系数等统计指标不够，最好能代入你的实际优化模型看决策差异。
4. MATLAB代码注意随机种子设置。场景生成有随机性，不固定rng种子，前后两次结果可能明显不同，论文复现性会受质疑。
5. 如果读者把场景数从10加到50，计算时间翻倍但误差改善甚微，说明关键信息已经被10个场景捕捉到了，没必要继续增加。

最后再分享一个小技巧：在削减过程中可以记录每一步的单步削减误差（即每次被删除场景的p_i × d_min(i)），画成误差曲线观察衰减趋势。如果发现中间某一步误差突然跳升，说明当前场景集中存在一个"孤立场景"——它与所有其他场景距离都很远，删它会造成很大分布损失。这种情况下保留它比多保留两个相似场景更有价值，你甚至可以把这种现象作为自适应选择目标场景数的依据。

这套方法我已经在多个项目里跑顺了，代码量不大，原理也不复杂，但每次用都能体会到"用数学优化替代拍脑袋削减"的踏实感。如果你正在为风光场景数量发愁，不妨把文中的代码跑一遍，应该能帮你少走不少弯路。