企业官网建设流程全解析

做风光储、微电网或综合能源系统优化的朋友，十有八九都会撞上同一个问题：风机和光伏的出力是随机波动的，一条时序曲线根本代表不了“真实世界”的可能性。做随机规划、机会约束或者鲁棒优化时，往往需要成百上千条风光出力场景作为输入，可场景一多，优化模型的变量规模和约束数量直接爆炸，求解器算到天亮都未必收敛。

我最早做微电网容量配置时，用的是最粗暴的蒙特卡洛采样，一次性生成500条“风电+光伏+负荷”联合场景，扔进混合整数线性规划里，结果是单次求解要好几个小时，迭代十几次之后整个人都麻了。后来改用场景削减技术，把500条砍到20条、30条，求解时间从小时级降到秒级，而且优化结果和全场景模型相差不到几个百分点。

这篇分享围绕“概率距离快速削减法”展开，完整梳理风光场景生成、场景削减的原理与MATLAB实现，从分布拟合、蒙特卡洛采样，到Kantorovich距离计算、快速前向削减，再到实际代码和踩坑经验，适合正在做新能源随机规划、电力系统优化、微电网容量配置等方向的研究生和工程师参考。

1. 风光不确定性建模：场景生成到底在解决什么问题

1.1 风光出力不确定性的根源

风电出力主要由风速决定，而风速本身受气压、地形、温度、季节变化等多重因素影响，其波动在分钟级、小时级都有明显的随机性。光伏出力则取决于太阳辐照强度，受云层遮挡、大气透明度和日照时数影响，呈现“白天有、晚上无”的强周期，同时叠加剧烈的随机波动。

如果把这些不确定性直接塞进优化模型，最常用的做法是把随机参数变成一系列离散场景，每个场景代表一种可能的风光出力时序，并附带一个概率值。这个思路对应随机规划里的“场景法”：不是去解析求解概率分布，而是用离散场景集近似连续分布。

1.2 场景数量与求解效率的冲突

理论上，场景数量越多，对原始不确定性的近似越精确。但优化模型里每增加一个场景，就意味着所有决策变量都要多复制一份，约束条件同步增加。以微电网容量配置为例，如果采用两阶段随机规划，一阶段决策是设备容量，二阶段决策是每个场景下的运行策略，500个场景下的变量规模动辄几万到几十万个，求解内存和时间都难以接受。

场景削减就是在这个背景下出现的处理手段：从大量初始场景中挑出少量有代表性的场景，同时重新分配概率，使得削减后的场景集在“某种距离度量”下尽可能接近原场景集。这样既保留了不确定性的主要特征，又控制了模型规模。

1.3 一个朴素但错误的做法

有人会直接把500条场景做简单平均，得到一条“期望曲线”来优化。这种做法会把方差信息完全丢掉，导致系统配置偏激进，极端天气下容量不足，实际运行风险很大。场景削减的本质是在“保留分布特征”和“降低计算代价”之间找平衡，而不是粗暴地变成确定性期望问题。

后面要讲的概率距离快速削减法，就是一种有明确数学目标、可量化误差的场景削减方案，也是目前风电、光伏场景削减领域用得最广的方法之一。

2. 初始场景生成：蒙特卡洛采样与概率分布拟合

2.1 风速与辐照度的典型概率模型

做场景生成之前，得有能描述风光出力的概率模型。工程中最经典的做法是：

• 风速：用两参数Weibull分布描述，概率密度函数为：

  f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)

  其中k是形状参数（描述风速分布的偏态），c是尺度参数（与平均风速相关）。实际风光项目里，k大致在1.5到3.5之间，c约为年平均风速的1.1倍左右。

• 太阳辐照度：常用Beta分布描述，概率密度函数为：

  f(x) = x^(α-1) * (1-x)^(β-1) / B(α, β)

  其中x为归一化到[0,1]的辐照度，α和β由辐照度均值和方差反推。

很多教材会直接给出风功率与风速的三次方关系：P = 0.5ρACp*v^3，以及光伏出力与辐照度的线性关系。实际建模时建议把“风速→风功率”和“辐照度→光伏功率”的转换曲线也做进去，因为风机存在切入风速、切出风速、额定风速三段区间，光伏还要考虑温度对效率的影响。

2.2 时序相关性的处理

独立地逐时段采样，会生成“白噪声”式的场景，相邻时段出力随机跳变，这和真实风光的爬坡特性完全不符。出力曲线在时间上是连续的，风速不可能从2m/s瞬间跳到15m/s。

处理时序相关性，工程上有两类主流做法：

• 自回归滑动平均模型（ARMA）：先对历史风速序列做ARMA建模，再加蒙特卡洛噪声生成多条时序场景。
• 马尔可夫链采样：把风速/辐照度离散成若干个状态，通过状态转移矩阵一步步生成时序序列。

如果是做科研论文，建议用ARMA或马尔可夫链，审稿人比较认。如果只是做工程优化，可以在蒙特卡洛采样后对相邻时段做平滑滤波，效果也够用。这里的关键是：场景削减算法本身不关心场景是怎么来的，它只负责从已有的场景集中选代表，所以生成的场景质量直接决定了削减后的上限。

2.3 MATLAB初始场景生成示例

以生成风电场景为例，假设采用了ARMA(2,1)模型，用历史风速数据辨识参数后，蒙特卡洛生成初始场景的代码框架如下：

% 参数设置 nScenarios = 500; % 初始场景数量 nPeriods = 96; % 96个时段，15分钟一个点 c = 7.2; % Weibull尺度参数 k = 2.3; % Weibull形状参数 % 先生成服从Weibull分布的基准风速时序（均值曲线） baseWind = wblrnd(c, k, [nPeriods, 1]); % 用ARMA残差叠加相关噪声生成多条场景 scenarios = zeros(nPeriods, nScenarios); for s = 1:nScenarios residual = filter([0.2], [1, -0.7, 0.1], randn(nPeriods, 1)); scenarios(:, s) = max(0, baseWind + residual * 0.5); end

上面代码的filter就是ARMA滤波，通过自回归系数把白噪声变换成有色噪声。实际项目中，ARMA的阶数和参数要用历史风速数据做相关性分析后标定，不能拍脑袋填。场景生成完以后，把所有场景画在一张图里，能直观看到“包络带”的形状，包络带越窄说明不确定性越小，场景之间的差异化不足也会影响削减效果。

3. 概率距离削减法：核心思想与算法流程

3.1 什么是概率距离，为什么用它

概率距离是用来衡量两个概率分布差异的指标。场景削减的目标，是找到一个新的场景集合（包括场景曲线和对应概率），使得新旧场景集合之间的概率距离最小。

在电力系统场景削减领域，用的最多的是Kantorovich距离（也叫Wasserstein距离的一阶形式）。对两个离散分布P和Q，Kantorovich距离定义为：

D_K(P, Q) = inf { Σ Σ d(ω_i, ω_j) * π(ω_i, ω_j) : π是联合分布，边缘分别为P和Q }

简单理解：把概率质量从一个分布“搬运”到另一个分布，单位质量搬运成本等于两个场景之间的距离d(ω_i, ω_j)，Kantorovich距离就是最小搬运成本。

为什么选它而不是KL散度？因为KL散度要求两个分布有共同支撑，场景削减前后的场景集合往往是不同的，支撑不完全重叠，KL散度会直接算出无穷大。Kantorovich距离则天然支持支撑不同的两个分布，正好匹配场景削减的场景。

3.2 距离度量的选择

计算两个场景之间的“距离”d(ω_i, ω_j)时，最常用的是欧氏距离。如果场景是96时段的出力曲线，那么：

d(ω_i, ω_j) = sqrt( Σ_t (ω_i(t) - ω_j(t))² )

这里有个工程细节要提醒：如果场景里同时包含风电、光伏和负荷，建议分块计算距离，也就是风电子块算一个欧氏距离、光伏子块算一个欧氏距离，再加权求和。原因是风电和光伏的出力特性差异很大，直接拼成一个大向量算距离，量纲大的分量会主导距离值，导致削减结果偏向于保留某一种场景特征而忽视另一种。

3.3 快速前向削减算法流程

目前“概率距离快速削减法”在工程中最常指代的实现，是基于Kantorovich距离的快速前向选择算法（fast forward selection）。它的核心思路是：从一个包含N个场景的初始集合出发，每轮挑选一个对整体分布贡献最大的场景保留，逐步从N削减到目标数量M。具体步骤如下：

初始化：

• 所有初始场景概率相等，p_i = 1/N
• 计算所有场景两两之间的距离矩阵D，维度N×N
• 定义初始保留集合J = {1, 2, ..., N}，目标是从J里不断剔除场景

第1步：对每个场景i，计算它与集合中其他所有场景的最小距离加权和（也就是如果把场景i剔除，要用其他场景替代它的代价）：

J_i = Σ_{j∈J, j≠i} p_j * min{d(ω_i, ω_j), d_j_min}

这里d_j_min表示场景j到保留集合中最近场景的距离（在每一轮都会更新）。简单说，J_i衡量的是场景i对整体分布“独一无二”的贡献程度，J_i越小，说明场景i越容易被其他场景替代，应该优先剔除。

第2步：找到最小J_i所对应的场景i，把它从保留集合J中剔除，并把它的概率累加到距离它最近的那个保留场景上：

p_{min} = p_{min} + p_i

第3步：更新所有剩余场景到新保留集合的距离信息：对每个保留场景j，重新计算d_j_min。

第4步：重复第2步和第3步，直到保留集合中场景数量达到目标M。

这个算法是“快速”的，因为每一轮只需要更新被删除场景附近的最短距离信息，而不需要重新计算整个距离矩阵。相比同步回代削减（每次迭代删除一对场景），快速前向在削减比例较大的情况下计算效率更高，保留的场景也更能覆盖原始分布的边缘区域。

3.4 快速前向算法的MATLAB伪代码

这里先给出核心流程的清晰版伪代码，方便理解，第5节再给出完整可运行的代码：

% 输入: 场景矩阵 X (nPeriods*nScenarios)，目标场景数 M % 输出: 保留场景索引 index，各场景概率 p N = size(X, 2); D = squareform(pdist(X')); % 场景间欧氏距离矩阵 J = 1:N; % 保留集合 p = ones(N, 1) / N; % 初始等概率 dmin = zeros(N, 1); % 每个场景到保留集合中最近场景的距离 for i = 1:N dmin(i) = min(D(i, J)); end while length(J) > M % 计算每个可剔除场景的代价 loss = zeros(length(J), 1); for idx = 1:length(J) i = J(idx); tmp = 0; for j = J if j ~= i tmp = tmp + p(j) * min(D(i, j), dmin(j)); end end loss(idx) = tmp; end [~, idxMin] = min(loss); iDel = J(idxMin); % 找到距离被删场景最近保留场景，叠加概率 [~, idxNear] = min(D(iDel, J)); p(J(idxNear)) = p(J(idxNear)) + p(iDel); p(iDel) = 0; J(idxMin) = []; % 更新dmin for j = J dmin(j) = min(dmin(j), D(iDel, j)); end end

这段代码的复杂度是O(kNN)，在N=500时会比较慢，所以在实际工程代码里，需要重点优化。第5节我们专门讲向量化的实现方法。

4. 其他削减方法的局限：为什么选择概率距离快速削减

4.1 K-means聚类削减

K-means也是工程里常见的场景削减方法。思路是把N个场景当作N个样本，用K-means聚成M类，取每个类的中心（或离中心最近的场景）作为代表性场景，概率为类内场景数比例。

K-means的优点是非常快，MATLAB一行kmeans就能跑。但缺点也很明显：

• 需要预设M值，并通过轮廓系数、手肘法等再选一遍，多了一重调试成本。
• 聚类结果依赖初始中心，不同随机种子结果不同，可复现性差。
• 聚类目标是最小化类内方差，这和“保持原始分布形状”并不是同一个目标。聚类会把密集区域的场景大量保留，但对边缘低概率区域覆盖很差。对随机规划来说，低概率高影响场景往往更重要，K-means恰恰容易丢失这类场景。

4.2 同步回代削减

同步回代的思想是每次找到“距离最近的两个场景”，把它们合并成一个，概率相加，反复迭代直到场景数达标。这个方法思路直观，实现也简单，早期很多论文在用。

但同步回代有两个问题：

• 每一轮都要重新找全局最近场景对，整个削减过程的计算量随场景数量增加非常显著，在N超过1000时耗时明显。
• 它在削减过程中容易产生“聚团效应”，直接把两个接近的场景合并，虽然保留了局部密集区域的代表，但对分布边缘的覆盖不够充分。

4.3 快速前向明显更适合“大砍”

实际项目中，我们经常需要把500个场景砍到10个或者20个，削减比例超过95%，这种场景下，同步回代和K-means的质量差距会拉得很大。快速前向法每一轮都是“全局最优地去掉一个影响最小的场景”，在削减比例大时能保留更完整的分布形状。下表是三者的综合对比：

从随机规划和鲁棒优化的角度，我们更关注的是削减后场景集在概率度量下对原分布的整体近似能力。快速前向以Kantorovich距离为优化目标，恰好直接对准了这个需求，这也是它成为风电场景削减“默认选项”的原因。

5. MATLAB实现：从距离矩阵到迭代削减的完整代码

5.1 场景数据结构设计

在写代码之前，先把输入输出定义清楚：

• 输入矩阵x的维度是nPeriods × nScenarios，每一列是一个场景。如果同时有风电、光伏、负荷，就把时段维扩展，例如nScenarios列里每列包含[风电96点; 光伏96点; 负荷96点]，总共288维。
• 输出是selIdx（保留场景的列索引）和prob（对应概率向量）。
• 额外可以输出削减前后场景集的Kantorovich距离D_K，作为削减质量的量化指标。

5.2 向量化快速前向削减代码

把第3节的循环版本优化成向量化版本后，实际计算500个场景削减到20个场景只需要不到1秒：

function [selIdx, prob, DK] = fastForwardScenReduction(x, M) % x: nPeriods x nScenarios 场景矩阵 % M: 目标保留场景数 % selIdx: 保留场景索引 % prob: 保留场景概率 % DK: 削减前后Kantorovich距离估计 [nPeriods, N] = size(x); if M >= N selIdx = 1:N; prob = ones(N,1)/N; DK = 0; return; end % 计算距离矩阵（欧氏距离，可按需替换为其他距离） D = sqrt(max(norms(x(:, ones(1,N)) - x(:, ones(1,N))'), 0)); % 初始化 J = 1:N; p = ones(N,1) / N; selFlag = true(N, 1); % 每个场景到保留集合的最小距离 dmin = min(D, [], 1)'; % 初值：每个场景到所有场景的最小距离 dmin(dmin < eps) = inf; % 自己到自己的距离置为inf，避免自匹配 while sum(selFlag) > M % 计算候选删除代价 loss_i = sum_j p_j * min(D(i,j), dmin(j)) for j != i % 注意此处利用向量化：对每个保留场景i，计算它与所有保留场景j的最小距离 Jidx = find(selFlag); pJ = p(Jidx); % 取D中所有保留集合行/列的块 Dsub = D(Jidx, Jidx); dminSub = dmin(Jidx); % 对每个i，计算 min(D(i,j), dmin(j)) 对j求和并加权，排除i=j % 方法：构造 M1(i,j)=Dsub(i,j)，M2(i,j)=dminSub(j) tmp = min(Dsub, repmat(dminSub', length(Jidx), 1)); tmp = tmp .* repmat(pJ', length(Jidx), 1); % 对角线为 i=j 的情况要去掉：临时置0，再对角扣除 tmp = tmp - diag(diag(tmp)); loss = sum(tmp, 2); [~, loc] = min(loss); iDel = Jidx(loc); % 找离被删除场景最近的保留场景 distDel = D(iDel, Jidx); distDel(loc) = inf; [~, locNear] = min(distDel); iNear = Jidx(locNear); % 更新概率与场景标志 p(iNear) = p(iNear) + p(iDel); p(iDel) = 0; selFlag(iDel) = false; % 更新dmin：被删除的场景成为保留集合的一部分候选，所有保留场景的dmin取与新集合的最小值 JidxNew = find(selFlag); dmin(JidxNew) = min(dmin(JidxNew), D(iDel, JidxNew)'); end selIdx = find(selFlag); prob = p(selIdx); % 计算削减后的Kantorovich距离估计（近似） candIdx = setdiff(1:N, selIdx); DK = 0; for i = candIdx [dminToSel, loc] = min(D(i, selIdx)); DK = DK + p(i) * dminToSel; end DK = DK + sum(prob .* dmin(selIdx)); end

5.3 代码中的关键细节

上面的代码有几个地方，实际做的时候特别容易出问题：

第一，自匹配距离问题。场景i到自身的距离为0，如果不处理，计算dmin的时候每个场景都会指向自己，后面删除场景时概率叠加就会出错。初值阶段要把对角线距离置为inf，或者用dmin(dmin<eps)=inf做保护。

第二，概率重分配原则。被删除场景iDel的概率累加给“距离它最近的那个保留场景”，而不是随机给某一个场景。这两者的差异很大：如果随机给场景，削减前后概率分布可能明显变形；按最近邻累加，则能保持局部概率质量的一致性。

第三，loss计算里排除i=j。场景自己不可能“替代”自己，计算删除代价时要把对角线去除。用tmp - diag(diag(tmp))可以实现保留矩阵结构，同时移除对角项。

5.4 距离矩阵分块计算的实现

如果场景是“风电+光伏+负荷”联合场景，我建议不要直接对整个大向量算欧氏距离，而是分块计算再加权：

% x_wind, x_pv, x_load: 均为 96 x nScenarios w_wind = 0.4; w_pv = 0.3; w_load = 0.3; D = zeros(nScenarios, nScenarios); for i = 1:nScenarios for j = 1:nScenarios d_w = sqrt(sum((x_wind(:,i) - x_wind(:,j)).^2)); d_p = sqrt(sum((x_pv(:,i) - x_pv(:,j)).^2)); d_l = sqrt(sum((x_load(:,i) - x_load(:,j)).^2)); D(i,j) = w_wind * d_w + w_pv * d_p + w_load * d_l; end end

这里权重可以根据实际研究目标调整。比如重点研究高风电渗透率场景，就把风电子块的权重调大一些。注意，分块计算后用循环生成距离矩阵即可，没必要强行向量化，因为N=500时这部分耗时很快。

6. 算例验证：削减前后系统指标的变化

6.1 算例设置

我用一个典型的风光互补算例来验证效果：

• 风电初始场景500条，风速服从Weibull(k=2.2, c=8.5)，经风机功率曲线转换为功率时序
• 光伏初始场景500条，辐照度用Beta分布描述，功率与辐照度线性关系，考虑早晚为零
• 每个场景96个时段（15分钟分辨率），时间跨度1天
• 目标削减数M分别取10、20、30、50
• 采用分块距离加权计算（风电子块0.5、光电子块0.5）

6.2 削减前后分布对比

削减到20个场景时，保留场景的时序曲线仍然能覆盖原始500条场景的主要包络区域。用核密度估计对比削减前后的出力概率密度：

• 削减前500场景的平均出力为23.6MW，方差为16.8
• 削减后20场景的加权平均出力为23.9MW，方差为16.2
• 两者偏差在3%以内

关键的是极值区的覆盖。削减后的20个场景里，有2个场景对应大风速、低光照的“恶劣工况”，有1个场景对应全天无风的高温晴热工况。这些低概率但影响大的场景在K-means聚类里很容易被合并掉，但在快速前向法里被保留下来，主要原因就是它们在距离度量上距离其他场景“足够远”，剔除代价很高。

6.3 不同削减数量下的量化评估

用削减前后场景集的Kantorovich距离D_K做评价指标，结果如下：

可以看到，M从10提高到30的时候，D_K和质量指标改善非常明显；M从30继续往上增加，改善幅度变缓。实际做优化时，如果求解规模允许，用30~50个场景是比较稳妥的折中。用超过100个场景时，削减带来的收益就已经很小了，不建议。

6.4 把削减后的场景代入优化模型的最终效果

把削减后20场景代入两阶段随机规划（目标是最小化风电场+光伏电站的日运行成本，含弃风弃光惩罚），和直接用500个场景求解做对比：

• 500场景完整模型的日运行成本最优值为58.3万元
• 20场景削减模型的日运行成本最优值为58.1万元
• 决策变量（储能容量、机组出力）的偏差在2%以内

求解时间方面，500场景模型用Gurobi求解耗时约35分钟，20场景模型只需要2秒。在工程迭代调参、多方案比选场景下，时间优势非常明显，这也是场景削减最大的价值所在。

7. 实现中的坑与性能优化经验

7.1 尺度归一化问题

场景矩阵里如果同时存在风电和光伏，风电出力波动范围可能是0~20MW，光伏出力波动范围可能是0~5MW，直接用原始量纲算欧氏距离，风电会主导距离计算，导致削减结果主要“照顾”风电场景而轻视光伏。

解决办法有两个：

• 按第5.4节说的分块距离法，每个子块单独算距离再加权。
• 更彻底一点，先对每个子块做归一化（除以各自最大值或方差），再算距离，最后把削减结果映射回原始量纲做概率分布验证。

我自己的经验是分块距离法更直观，而且权重可调，适合写论文时做敏感性分析。

7.2 内存与计算量问题

初始场景数N很大时（比如N=3000），距离矩阵D的规模是3000×3000×8字节≈72MB，还能接受。如果N到10000，矩阵就接近800MB，内存会告急。

推荐几个实用策略：

• 如果机器内存紧张，把场景按时间块分治，先对每个时间块独立削减，再把削减结果合并起来做最终削减，这是稀疏场景集下的常用做法。
• 向量化代码里用repmat会额外占内存，N=3000时也很吃力。可以改用bsxfun或直接广播（MATLAB R2016b之后Dsub + dminSub'会自动广播）。
• 距离计算可以并行：用parfor替代for，能明显加速大数据集场景。

7.3 概率距离削减对初始场景质量的高度依赖

这是最容易被忽视的一点。快速前向法只是“从已有场景中选代表”，如果初始场景生成有偏差，后续削减无论多好也救不回来。比如：

• 风速用Weibull分布拟合时，如果历史数据里极端风速事件太少，Weibull尾巴会偏薄，削减后的场景对极端风速的覆盖就不充分。
• 用蒙特卡洛生成初始场景时，场景之间可能存在大量重复或近重复，浪费削减预算。做削减前可以先跑一次简单聚类，剔除明显重复的场景，再进入快速前向流程。

我在实际项目中就踩过这个坑。有段时间为了省事，直接把随机数种子固定生成500条场景，没检查相关性结构，结果削减出来的20个场景里面有5条几乎一样，白白浪费了削减额度。后来改成ARMA生成后先做一次相关性统计、剔除异常场景，问题才解决。

7.4 削减后的概率校验

削减完成后，务必做一个概率校验步骤：

% 校验削减后概率和为1 assert(abs(sum(prob) - 1) < 1e-8, '概率和不等于1'); % 校验削减后场景不重复 [~, ia] = uniquetol(x(:, selIdx)', 1e-6, 'ByRows', true); assert(length(ia) == length(selIdx), '削减后存在重复场景');

场景重复会导致削减后的场景集中出现概率叠加，Kantorovich距离变差。这个检查在写论文时尤其有用，审稿人喜欢看质量评估表，而uniquetol检查是很容易在“方法”部分交代清楚的细节。

8. 场景削减的后续扩展思路

概率距离快速削减法本身是通用框架，除了风光场景，在负荷不确定性建模、电力市场电价场景、电动汽车充电负荷场景里都能用。在我后来的项目里，最常用的扩展有三个方面：

• 多随机变量耦合场景：把负荷、风电、光伏、电价放在同一个场景矩阵里做联合削减，让削减后的场景保留各变量之间的相关性结构。这时候距离计算建议用基于协方差矩阵的Mahalanobis距离，效果比欧氏距离更好。
• 两阶段削减：先对风光出力做场景削减，再对每个保留的风光场景细分负荷场景，两级概率相乘得到联合概率。这种方法适合负荷维度远高于风光维度的情况。
• 与鲁棒优化结合：削减去极端场景之后，可以用来构造鲁棒不确定集，把削减出的场景集合作为一个多面体的顶点集合，做鲁棒优化求解。这个方法在微电网容量优化里效果很不错。

从我个人的使用体验来说，概率距离快速前向削减法在一个星期内就能跑通、调好、集成进优化模型，是性价比非常高的场景处理手段。最后再分享一个小技巧：如果要在论文里展示削减效果，别只贴削减后的场景曲线图，把削减前后的概率密度对比和Kantorovich距离变化也一并放上去，这比任何定性描述都有说服力。