企业官网建设流程全解析

搞风光出力的不确定性建模，有一件事绕不开：你手头要么是一条确定性预测曲线，要么是一堆散布在时序空间里的随机样本，但优化模型真正需要的，是一个规模可控、概率明确、损失函数还能算得动的离散场景集。我自己最早是被一个两阶段随机规划卡住的——确定性出力没法表达不确定性，而直接把2000条蒙特卡洛场景塞进去，求解器直接罢工。后来改用基于概率距离的快速削减法（Fast Forward Selection）做风光场景生成与削减，才把问题理顺。

这套流程现在基本成了我做随机优化的标配：先按风速/光照分布或历史出力特征生成足够多的初始场景，再用概率距离做削减，最后留下十几个带概率的典型场景供优化模型使用。整个过程在MATLAB里可以完整跑通，生成的场景集既能保留原始分布的主要特征，又不会让优化模型爆炸。这篇就按我的实际代码和踩坑记录展开，原理、代码、参数、排错都会涉及，适合正在做电力系统不确定性分析、微网调度、随机规划相关课题的同学参考。

1. 风光场景问题：为什么非得“生成”和“削减”

1.1 随机出力的连续分布怎么塞进优化模型

风光出力的本质是随机过程。风速受地形、气压、季节影响，光照受云层、大气透明度影响，这些因素叠加起来，出力曲线根本没有解析表达式。就算你从历史数据里拟合出了风速的Weibull分布、光照的Beta分布，也只能得到一个连续概率模型。

连续分布虽然数学上漂亮，但优化模型不认。线性规划、混合整数规划处理的是有限维决策变量，你要么把随机性写进约束，要么把随机性离散成有限个“可能世界”。前者对应鲁棒优化——需要构造不确定集合；后者对应随机规划——需要构造场景集。场景这个词听起来玄，说白了就是一组组“某天可能出现的光伏/风电出力时序曲线”，每条曲线配一个发生概率。

所以第一步“生成”解决的是从连续分布到离散样本的转化问题。蒙特卡洛采样最简单，但收敛慢，而且纯随机采样容易在概率密度低的区域浪费样本。拉丁超立方采样（LHS）能让样本更均匀地覆盖分布空间，是我在初始场景生成时优先用的方法。

1.2 削减的本质是分布近似，不是选几条好看的曲线

生成完初始场景之后，问题就来了：样本数量动辄上千，每个场景都是一条24小时或96点的时序曲线。把这些全塞进两阶段随机规划，决策变量会膨胀到求解器无法忍受；如果塞进机会约束规划，场面会更难看。实际操作中，保留5到20个场景是比较舒服的区间。

那怎么从1000个场景里挑出10个？直觉做法是K-means聚类，聚出10个簇，取簇中心。但这里有个关键区别：场景削减不是“找代表性曲线”，而是调整一个离散概率分布，让它尽量接近原始的离散分布。每个初始场景本身带着概率权重（通常等概率1/N），削减后少数场景要继承这些概率。判断削减好坏的标准，不是削减出的曲线“看起来像不像”，而是两个离散分布之间的概率距离有多小。

这个认识上的转变很重要。我见过很多人直接调K-means，把10个簇中心当成削减结果，然后随便给个等概率——分布形状完全变了，随机优化结果自然不可信。用概率距离做指导，每条保留场景吸收周围场景的概率，本质上是在做一个离散分布的最优逼近。

2. 快速前向削减法到底是怎么回事

2.1 场景间的“距离”怎么定义：平方欧氏代价与Wasserstein距离

要衡量一个场景集对另一个场景集的逼近程度，先得定义单条场景之间的代价。因为每条场景都是T维向量（比如24个时刻的归一化出力），最自然的度量就是欧氏距离。实际使用中通常取平方欧氏距离，也就是空间上较远的场景匹配时产生更大的代价。

两个离散分布之间的p阶Wasserstein距离（也叫Kantorovich距离）定义是：W_p = ( Σ_i p_i · min_j c(x_i, y_j) )^(1/p)，其中x_i来自源分布，y_j来自目标分布，c是代价函数。当c取欧氏距离的p次方时，W_2就是最常用的形式。这个公式直观理解是：把源分布中每个点的概率“搬到”目标分布最近的点上，整体搬了多少“土方量”，就是两个分布的距离。

我当初选Wasserstein而不是KL散度，原因很实际。KL散度要求两个分布的支撑集重叠，否则直接算成无穷大；而且它不对称，拿削减后分布去逼近原始分布和反过来算，结果完全不一样。Wasserstein距离对支撑集不敏感，天然适应“场景数变少、支撑点变少”的场景削减场景。它就像搬土方——把A分布的一堆土搬到B分布的位置，最少要搬多远，直观、对称、对离群点也不那么歇斯底里。

2.2 前向选择的迭代逻辑：每一步多保一个最值得保的场景

快速前向削减法（Fast Forward Selection）的思想很直接：从空集开始，每次从还没被保留的候选场景里挑一个，使得“加入它之后，整个原始分布到已选场景集合的Wasserstein距离下降最多”。重复这个过程，直到保留场景数达到目标值。

具体迭代时，代码里维护一个dist_to_selected向量，记录每个候选场景到当前已选集合的最小代价距离。设当前已选集合为S，未选集合为J，任意候选j∈J加入后的整体匹配代价是： Σ_i p_i · min( dist_to_selected(i), c(i,j) ) 其中对已选场景i，dist_to_selected(i)是0，它们不贡献代价；对候选j自己，c(j,j)=0，也不贡献。所以这个式子算的就是：如果保留S∪{j}，所有原始场景把概率转移到最近的保留场景时，产生的加权代价总和。每一轮遍历所有未选场景，取这个代价最小的那个加入S，然后更新dist_to_selected。

第一轮比较特殊，因为dist_to_selected全是无穷大，这时算法等价于挑一个使加权平方距离之和最小的场景，也就是加权“几何中心”附近的场景。所以第一个保留场景往往是主模态，后续轮次才逐步往外围补场景。这是很合理的行为。

2.3 和同步回代、K-means削减的比较，什么时候用哪个

场景削减还有两个常见路线。同步回代削减（Backward Reduction）是从完整集合出发，每次删掉一个场景，使得删除后概率距离增量最小，直到剩下目标数量。K-means/谱聚类则是把场景当成普通样本聚簇，取簇中心并统计簇内概率。

我实际用下来的选择标准是：

如果追求可复现和理论严谨，前向选择是首选。它的每一步都直接对准Wasserstein距离下降，最后拿到的场景集和原始分布的匹配是有量化误差的，而这个误差在后续评估里还能继续用。

3. 初始场景生成：怎么造出足够可信的“母场景”

3.1 从风速Weibull分布和功率曲线出发的生成流程

场景削减再厉害，初始场景太烂也是巧妇难为无米之炊。最经典的风电场景生成路线是：风速用两参数Weibull分布拟合，再通过风机功率曲线把风速转成出力。

Weibull分布的两个参数k（形状）和c（尺度）通常用历史风速数据的极大似然估计求出来。有了分布参数，用蒙特卡洛或LHS抽N组风速时间序列。但要注意，不同时刻的风速不是独立同分布的，白天和夜间、相邻小时之间存在相关性。简单做法是引入AR(1)模型：v(t) = φ·v(t-1) + ε(t)，其中φ是自回归系数，ε是高斯白噪声。这样生成的风速序列在时间上有连续性，不会出现8点风速15m/s、9点风速1m/s这种离谱序列。

风速转功率时用标准功率曲线：

• 低于切入风速v_ci，出力为0；
• 切入风速到额定风速v_r之间，按线性或非线性曲线上升；
• 额定风速到切出风速v_co之间，输出额定功率；
• 高于切出风速，停机保护，出力为0。

光伏场景类似，只是先把光照强度按Beta分布建模，再用光照-功率线性模型加转换效率折算。这里不再展开太多，重点是想说明初始场景生成必须保留时序相关性和物理约束。

3.2 基于历史出力加AR扰动保留时序相关性

上面那种方法需要风速数据和功率曲线参数，有时候手头只有历史出力序列，那就换个思路。先拿历史出力算均值曲线，再用核密度估计或经验分布描述每个时刻的波动范围，最后生成“均值曲线+扰动”的样本。扰动部分我习惯加一阶自回归项，让整条曲线的形态更接近真实日曲线，而不是每个时刻独立抖动出来的脉冲。

时间相关性为什么重要？削减算法用的是整条时序曲线的欧氏距离，如果原始场景里相邻时刻是随机跳变的，距离度量就会被高频噪声主导，削减出的场景也会跟着毛躁。实际做微网调度时，大家关心的是日内爬坡、净负荷峰值、储能充放电窗口，这些都需要时序曲线平滑可解释。用AR扰动生成初始场景，比逐时刻独立采样的效果好一个档次。

3.3 初始场景数量怎么定

初始场景太少，削减算法等于在矮子里拔将军；太多，距离矩阵内存和计算时间就上来了。我自己的经验分几个档位：

• 初步验证算法：N=500，够用，跑得快；
• 常规论文/项目场景：N=1000到2000，比较稳妥；
• 追求分布精度：N=5000，这时候必须用矩阵化代码和单精度距离矩阵，否则MATLAB会卡死。

初始场景数N和最终保留场景数R没有严格公式，但经验法则是R不要超过N的5%，也不要少于3个。低于3个场景，随机优化基本失去意义。

4. MATLAB代码拆解：核心函数我这样写

4.1 文件组织与主脚本流程

整个流程我分成三个文件：主脚本、场景生成函数、快速前向削减函数。另加一个Wasserstein距离评估函数，用来量化削减质量。

先给主脚本，输入输出都写在注释里，跑一遍就能看到削减前后的对比图。

%% main_scenario_reduction.m % 基于概率距离快速削减法的风光场景生成与削减 % 流程: 生成初始场景 -> 快速前向削减 -> 概率再分配 -> 质量评估 clear; clc; rng(42); N = 1000; % 初始场景数 T = 24; % 时间点数 N_red = 10; % 保留场景数 % 用一条基准日出力曲线加AR扰动生成初始场景 % 这里模拟的是归一化后的出力序列，范围[0,1] base = 0.55 + 0.30 * sin((0:T-1)/T*2*pi - pi/2); base = max(base, 0.05); phi = 0.85; % 自回归系数，控制时序平滑性 sigma_noise = 0.08; % 扰动强度 scen = zeros(N, T); for n = 1:N noise = zeros(1, T); for t = 2:T noise(t) = phi * noise(t-1) + sigma_noise * randn(); end scen(n, :) = base + noise; end scen = min(max(scen, 0), 1); % 钳位到[0,1]，对应归一化出力约束 prob = ones(N, 1) / N; % 初始场景等概率 % 快速前向削减 [red_scen, red_prob, selected_idx, hist_wass] = fast_forward_reduction(scen, prob, N_red); % 削减前后概率距离评估（用平方欧氏代价的平方根，对应W_2） wass_before = wasserstein_distance(scen, prob, scen(1,:), 1.0, 'one_scene'); wass_after = wasserstein_distance(scen, prob, red_scen, red_prob, 'multi'); % 可视化 figure('Color','w','Position',[100 100 900 420]); subplot(1,2,1); plot(scen(1:100,:)', 'Color',[0.7 0.7 0.7]); hold on; plot(red_scen', 'LineWidth',1.5); title('初始场景与削减后场景'); xlabel('时刻'); ylabel('归一化出力'); subplot(1,2,2); bar(1:N_red, red_prob); title('削减后场景概率'); xlabel('场景编号'); ylabel('概率'); disp('削减前后Wasserstein距离对比:'); disp([' 削减前(用单场景近似全分布): ', num2str(wass_before)]); disp([' 削减后(用10场景近似全分布): ', num2str(wass_after)]);

这里为了演示，我直接用基准曲线加AR噪声代替了Weibull采样流程。真实项目里只要把scen换成generate_wind_scenarios的返回值就行，后面的削减逻辑完全不用改。

4.2 fast_forward_reduction逐行解析

核心削减函数是全篇的关键。我写的时候充分考虑了可读性和性能的平衡，没有用特别花哨的MATLAB技巧，但每一步都有明确注释。

function [red_scen, red_prob, selected_idx, hist_wass] = fast_forward_reduction(scen, prob, N_red) % 基于概率距离的快速前向削减法 % 输入: % scen : N x T 初始场景矩阵，每行一条时序曲线 % prob : N x 1 场景概率向量，和为1 % N_red : 保留场景数 % 输出: % red_scen : N_red x T 削减后的场景 % red_prob : N_red x 1 削减后的概率 % selected_idx : N_red x 1 保留场景在原始矩阵中的行号 % hist_wass : N_red x 1 每一轮迭代后的Wasserstein代价 [N, T] = size(scen); if N_red >= N error('削减后的场景数必须小于初始场景数'); end if abs(sum(prob) - 1) > 1e-6 warning('概率之和不为1，建议检查输入概率向量'); prob = prob / sum(prob); % 防御性归一化 end % 预计算距离矩阵。N=5000时是200MB double，可接受。 % 若N>5000，建议改用单精度或分块计算，见常见问题 D = pdist2(scen, scen, 'squaredeuclidean'); % 平方欧氏距离 D(1:N+1:end) = 0; % 对角线置零 selected = false(N,1); % 是否已选中 dist_to_sel = inf(N,1); % 每个场景到已选集合的最小距离 hist_wass = zeros(N_red,1); for r = 1:N_red cand_val = inf(N,1); for j = 1:N if selected(j) continue; end % 如果选中j，每个场景取“到已选集合”和“到j”的较小距离 d_new = min(dist_to_sel, D(:, j)); cand_val(j) = sum(prob .* d_new); % 加权Wasserstein代价 end [best_val, best_k] = min(cand_val); % 更新选中集合和距离结构 selected(best_k) = true; dist_to_sel = min(dist_to_sel, D(:, best_k)); hist_wass(r) = best_val; end % 提取保留场景 selected_idx = find(selected); red_scen = scen(selected_idx, :); % 概率再分配：每个原始场景的概率转移到距其最近的保留场景 D_red = pdist2(scen, red_scen, 'squaredeuclidean'); [~, nearest] = min(D_red, [], 2); red_prob = accumarray(nearest, prob, [N_red 1]); red_prob = red_prob / sum(red_prob); % 防御数值误差，再归一化 end

有几个实现细节值得单独说明。

第一，距离矩阵用了平方欧氏距离。如果这里用普通欧氏距离，Wasserstein代价的尺度会偏小，削减结果对远距离场景的惩罚不够，场景分布容易挤在中心附近。平方欧氏距离把远处的代价放大，保留场景会更均匀地铺满整个分布空间。

第二，D(1:N+1:end)=0这一行把对角线清零。如果不处理，候选j到自己的距离是0，但不影响min运算；处理之后逻辑更清晰。

第三，概率再分配没有用复杂的最优传输，直接按最近保留场景转移。因为前向选择本来就是基于最近匹配的，最近邻转移和Wasserstein距离的匹配方式一致，实现简单且结果正确。

4.3 概率再分配与最终评估函数

评估函数用来量化削减质量，我单独写成wasserstein_distance，方便对比不同保留数量下的削减误差：

function wass = wasserstein_distance(scen, prob, ref_scen, ref_prob, mode) % 计算两个离散分布之间的Wasserstein距离(W_2) % 输入: % scen : N x T 源分布场景 % prob : N x 1 源分布概率 % ref_scen : M x T 参考分布场景，M可等于1 % ref_prob : M x 1 参考分布概率，模式为one_scene时可传1 % mode : one_scene 表示单场景近似，multi表示多场景 % 输出: % wass : W_2距离 [N, ~] = size(scen); [M, ~] = size(ref_scen); % 计算源场景到参考场景的最小代价 D = pdist2(ref_scen, scen, 'squaredeuclidean'); % M x N [dmin, ~] = min(D, [], 1); % 每个源场景的最近匹配代价 if strcmp(mode, 'one_scene') % 单场景近似时，参考分布概率就是1，直接加权 wass = sqrt(sum(prob(:) .* dmin(:))); else % 多场景近似时，参考分布的概率在削减函数里已经算好， % 这里通过ref_prob对匹配方向做对称校验（可选）。 % 实际使用时可以只算单向距离，因为削减优化目标就是单向匹配 wass = sqrt(sum(prob(:) .* dmin(:))); end end

注意一个哲学问题：Wasserstein距离严格来说是对称的，但场景削减的优化目标是“原始分布到削减分布的匹配代价”，是单向距离。实际评估时用单向就够了。如果非要算对称距离，还需要反过来把削减分布的概率搬到原始分布上，那样会增加不少计算量，而且和削减目标不一致。

4.4 向量化与内存处理的取舍

上面代码里最重的部分是预计算D矩阵，复杂度是O(N^2·T)。N=1000，T=24时，pdist2计算飞快；N=5000时就需要一点时间，但还在可接受范围。真正的问题在内存：N=5000时D是5000×5000的double，约200MB，普通笔记本没问题；N=10000时逼近1GB，MATLAB会很吃力。

如果场景规模特别大，我一般做三件事。一是用single类型存距离矩阵，内存直接减半，代价是极小概率的数值误差。二是分块计算距离，不要一次pdist2全部求完，每块500行，边算边用。三是干脆不存完整距离矩阵，每轮迭代动态算D(:,j)的某几列，牺牲速度换内存。

代码里给的还是完整矩阵版本，因为对大多数课题场景来说，N=1000到2000已经够用，没必要为极端规模牺牲可读性。

5. 一次完整实操：从1000个场景削减到10个

5.1 参数设定与跑通效果

我在实际跑流程时，最常用的一组参数是：N=1000，T=24，N_red=10。这个规模既有说服力，也不会让等待时间超过几秒。

用第4节主脚本直接跑，第一次输出会看到迭代过程中Wasserstein代价的变化。初始时单场景近似全分布的代价大约在一个较大数值（具体取决于基准曲线和噪声强度），第一轮选中一个近似中心场景后，代价立刻大幅下降，之后每多保留一个场景，代价都会继续下降，但下降幅度逐渐趋缓。

实际操作时要注意，MATLAB输出里如果出现削减后某条场景的概率特别小，比如0.005以下，说明这个场景只代表极少数极端样本。这种场景对期望值目标贡献很小，但在风险规避型优化里可能恰恰是尾部风险的代表，要不要保留取决于你的优化目标。我通常不会因为概率小就删掉它，除非总概率低于0.1%。

5.2 削减过程的可视化与迭代观察

主脚本第一张图画的是前100条初始场景和削减后的10条场景叠加。这里能看到一个很有意思的现象：削减出的10条曲线不是简单地把100条曲线叠在一起取包络，而是每条都落在初始场景分布的不同“密度峰”附近。有的贴近基准线（大概率中心），有的明显偏离（代表高风/低风日），还有一条贴近上界附近的极端场景。

第二张概率条形图更直观。按降序排列后，概率最大的场景往往占30%-40%的总概率，中间几个场景各有5%-15%，尾巴上有一两个场景概率低于5%。这个分布形态和原始数据本身的离散程度强相关：如果基准曲线附近聚集了60%的样本，削减后大概率会把大部分概率合到一到两个中心场景上。

5.3 削减结果质量怎么量化评估

我习惯做一个消融实验：把N_red从1扫到30，记录每个数量对应的Wasserstein距离，画一条下降曲线。曲线会在某个位置出现拐点，再增加场景数，距离下降已经不明显。这个拐点就是适合你问题的保留场景数。

用第4节的评估函数算出削减前后距离，你会发现削减到10个场景时，距离可能只有单场景近似时的20%-30%，这就是场景削减带来的信息压缩效率。如果下降曲线很平缓，说明原始分布本身比较集中；如果下降很慢，说明分布很分散，可能需要保留15到20个场景才够用。

6. 常见坑位与排查记录

6.1 削减出来的概率总和不对

概率总和不对基本是输入prob没归一化，或者accumarray的索引和概率维度对不上。我在代码里加了防御性归一化和warning，但如果你自己改代码，最容易翻车的是nearest的数量不等于N，或者prob里有NaN。

排查方法很简单：削减后直接sum(red_prob)，如果不是1，先查prob，再查D_red计算是否正确。另一个隐蔽问题是场景矩阵里存在NaN值，pdist2遇到NaN会直接返回NaN，导致min运算结果全是NaN。

提示：生成场景时如果做了插值或缺失值修补，务必确认没有NaN残留。我踩过一次这个坑，排查了半天才发现是某条风速数据缺测导致整行场景污染。

6.2 距离矩阵内存爆了怎么办

N=5000以上时，完整D矩阵可能让你的MATLAB直接内存不足。我建议按顺序尝试：先转single，内存减半；再分块算距离，参考这个思路：

nblocks = ceil(N / 1000); D = zeros(N, N, 'single'); for bi = 1:nblocks idx_i = (bi-1)*1000+1 : min(bi*1000, N); D(:, idx_i) = single(pdist2(scen, scen(idx_i,:), 'squaredeuclidean')); end D(1:N+1:end) = 0;

分块不会改变结果，只是把峰值内存降下来。如果连single分块都撑不住，那就只能在迭代里动态计算，每轮只算候选j对应的列。

6.3 场景削减后形态失真

如果削减出的场景曲线出现高频抖动，或者和原始场景形态差异过大，问题一般出在初始场景生成阶段。最典型的原因是逐时刻独立采样导致时间相关性丢失，这时削减算法会把那些“锯齿形”场景当成独立模态保留下来。

处理方法是在生成阶段加入AR项或滑动平均，保证相邻时刻的出力平滑。另一个办法是在削减时把距离代价从平方欧氏改成带时间权重的高阶距离，但那是特殊需求，常规场景用平方欧氏就够。

6.4 保留场景数N_red怎么选

没有绝对答案，但有参考流程。我的习惯是跑一遍距离扫描曲线，找拐点；同时结合优化模型规模，优先保证求解器能跑动。对两阶段随机规划，10个场景通常是甜点区；对机会约束或CVaR类风险约束，15到20个场景更稳，因为尾部事件的分布覆盖更充分。

另一个提示：N_red不是越大越好。随机规划求解时间和场景数近似线性增长，20个场景比10个场景慢了不止一倍，而Wasserstein距离可能只下降百分之几。花50%的时间换取3%的精度，性价比很低。

7. 把场景集用到随机优化里的一些经验

7.1 两阶段随机规划中的接入方式

削减后的场景集最直接的用途是两阶段随机规划。典型形式是：第一阶段决策在不确定实现前做出（比如机组启停），第二阶段根据场景决策调整（比如实际出力、储能出力）。场景数和概率直接进入目标函数：

min c^T x + Σ_s p_s · Q(x, ξ_s)

其中ξ_s就是削减后的场景，p_s是它的概率。这块模型不复杂，但有两个经验值得分享：一是场景概率最好在模型外归一化好再传入，不要在YALMIP或求解器里再除一遍；二是场景数量不要超过求解器能承受的上限，否则每次迭代求解都会变成灾难。

7.2 场景树和多阶段问题扩展

如果做的是多阶段决策，单层削减不够，需要构造场景树。思路是前几个时刻用少量场景，后几个时刻在每个分支下再生成子场景。快速前向削减法天然支持这个扩展：先对整条时序的前T1时段做削减，再在保留场景分支下对后续时段做条件削减。这样构造出的场景树既控制了总节点数，又保留了跨时段的相关性。

我做过一个模拟项目X的多阶段机组组合，三层场景树（根节点分3个、二层各分3个、三层各分4个）总叶子数36个，比直接生成1000个完整场景再整体削减更能体现时序决策的序贯特性。

7.3 鲁棒与分布鲁棒方向的后续玩法

场景削减产出的不只是场景，还顺带定义了“概率距离”这个度量。这个度量在分布鲁棒优化里可以直接用——拿Wasserstein球构造模糊集，半径就是削减前后距离。这个方向的好处是，你不需要知道真实分布，只需要一个名义分布（历史数据或削减后场景集）和一个球半径，就能在覆盖所有近邻分布的意义下做最坏情况优化。

这个扩展现在很热，核心思想就是把场景削减的误差边界显式写进优化模型。如果你的项目刚好卡在“场景数太多不好解、太少怕不准”的困境，这条路线值得研究。

我个人在实际操作中的体会是：场景生成与削减这套流程，参数不用一上来就追求完美。先用N=500、N_red=10跑通全流程，把代码和可视化框架搭好，再回来调Weibull参数、AR系数、保留数量。顺序反了会非常折磨人，因为问题的源头往往在生成阶段而不是削减阶段。另外，削减后一定要看概率分布和距离下降曲线，别只盯着曲线形状。场景集的价值在于它能代表原始不确定性，而不在于它看起来多漂亮。