做微网动态经济调度的人，基本都会撞上同一个拦路虎：风机出力、光伏波动、负荷预测偏差这些不确定性，到底该怎么塞进优化模型里。如果你只取期望值，调度方案往往过于乐观，遇到极端天气直接抓瞎；如果你把所有可能场景都枚举进去，模型维度爆炸，求解器算到天亮也出不来。场景生成与场景削减这套组合拳，就是为了解决这个矛盾——先用概率分布把不确定性描述清楚，再通过抽样和聚类把它压缩成有限几个代表性场景，最后放进随机优化模型里求一个兼顾经济性与鲁棒性的调度方案。
这篇内容适合正在做微电网、主动配电网经济调度、或者接触随机规划但总觉得理论到代码隔着一层纱的人。我会把从概率分布建模、场景生成、场景削减到随机优化建模的完整链路拆开讲，附上Matlab里的实操代码和参数设置，包括那些文献里不会写、但你在复现时一定会踩的坑。这里不讨论具体某篇论文的案例抄作业，重点是把这套方法论吃透——换任何系统、任何概率分布，你都能搭出自己的调度框架。
1. 核心思路：为什么微网调度绕不开“场景”
1.1 场景法在整个调度链条里的位置
微网动态经济调度，本质上是一个多时段的资源分配问题。传统确定性调度把光伏、负荷都当成固定值，模型里全是等式和不等式约束，求解快，结果也漂亮。问题是现实世界里没有“固定值”——云飘过来，光伏出力两分钟内掉一半，这个出力值在调度时刻到来之前，是不可能精确预知的。
处理不确定性的主流思路有三条：鲁棒优化、机会约束规划、以及本文要展开的随机优化（场景法）。鲁棒优化追求最坏情况下的可用性，但往往过于保守，经济性很差；机会约束允许一定概率的越限，处理起来数学上比较绕。场景法的思路更直接：既然不确定量服从某个概率分布，那就从这个分布里抽出一批代表性样本，每个样本带一个发生概率，然后把所有样本推进优化模型，目标函数是“所有场景下的期望成本最小化”。这样一来，不确定性就被自然地转化成了确定性的、但规模更大的优化问题。
这套“抽样本—算概率—建模型”的流程，就是整个技术链条的骨架。场景生成解决“样本从哪里来”，场景削减解决“样本太多怎么办”，随机优化解决“样本怎么用”。在实际的微网调度中，这套流程通常和模型预测控制（MPC）的滚动优化配合使用——每次滚动窗口开始前，重新生成最新的场景集，再求解一次优化问题，实现闭环校正。
1.2 五种不确定性来源的建模取舍
微网里值得建模的不确定性来源，我归纳下来主要是五类：风速、光照辐照度、负荷预测误差、电价（如果参与市场化交易）、以及电动汽车充电等随机行为。其中前四类在文献里最常见，处理手段也相对成熟。
风速通常用两参数威布尔分布描述，形状参数k和尺度参数c可以通过历史数据的极大似然估计得到。光照辐照度在晴天时呈Beta分布特征，但遇到多云天气就变成多峰，简单Beta往往不够用，实践里更多是直接用历史样本的核密度估计。负荷预测误差一般假设为零均值正态分布，标准差取预测值的5%到10%，这个假设在工程上足够用了。电价序列则复杂得多，通常需要时间序列模型，比如ARIMA或者场景树生成。
我在实操中的建模取舍原则是这样的：能查到历史数据的，优先用统计分布拟合；数据量不足的，宁可用带偏置的保守分布，也不强行套标准分布。很多新手一上来就用标准正态分布描述所有不确定性，结果生成出来的场景和实际出力规律严重脱节，后面再怎么削减都救不回来。
2. 场景生成实操：从概率分布到出力时序
2.1 风速、光照与负荷误差的分布建模
先看风电出力。风机的功率转换曲线通常是一个分段函数：切入风速以下出力为零，额定风速到切出风速之间出力恒定，中间区域按三次方关系近似。因此要生成风电出力场景，第一步是生成风速场景，再通过功率曲线映射到出力序列。
风速的威布尔分布概率密度函数是：
f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k决定分布形状，c决定特征尺度。很多地区的k值在1.8到2.3之间，c值大约等于该地平均风速的1.12倍左右。这个参数如果懒得从历史数据里拟合，先用k=2、c按平均风速估算也能跑通流程，但做实际项目时，强烈建议拿至少一整年的实测风速数据去做极大似然估计。Matlab里可以用wblfit函数直接拟合，传入历史风速序列即可得到k和c。
光伏出力和风速逻辑类似但模型不同。太阳辐照度在小时尺度上更接近Beta分布，其概率密度函数为：
f(x) = (x^(a-1) * (1-x)^(b-1)) / B(a,b)Beta分布的两个形状参数a和b，可以由辐照度历史数据的均值和方差反推。但这里有个关键点：Beta分布的定义域是[0,1]，当辐照度转化为标幺值后可以直接套用，要记得乘以安装容量和转换效率才能得到有功出力。实际项目里还有一个更省事的方案——直接对历史光伏出力数据做拟合，跳过辐照度这个中间变量，因为辐照度到出力的转换效率本身就有不确定性，绕开它反而能减少一层建模误差。
负荷预测误差相对简单，直接假设服从均值为零、标准差为预测值一定百分比的正态分布。实际调度系统里预测误差的标准差一般取5%左右，保守场景可以取到10%。如果在做的是配电网层面的调度，负荷的波动性会比单一馈线小，误差标准差可以适当压缩。
2.2 拉丁超立方采样：比蒙特卡洛更聪明的抽样
有了概率分布函数，下一步就是从分布里抽样本生成场景。最朴素的做法是蒙特卡洛随机抽样——按分布随机生成几千个风速样本，每个时刻都独立抽样。这个方法简单，但有个致命的低效问题：要想让样本空间覆盖均匀，尤其是尾部分布（也就是极端大风、极端低负荷这类小概率但高影响场景），需要的抽样次数大得惊人。几千次抽样可能只有个位数落到尾部区间，这些极端场景对调度鲁棒性恰恰又是最关键的。
拉丁超立方采样（LHS）是更聪明的做法。核心思想是分层抽样——把概率分布函数的值域[0,1]均匀切成N个等概率小区间，从每个区间里随机取一个分位数，再通过累积分布函数的逆变换映射回变量取值。这样无论样本量多大，都能保证变量在整个值域范围内被“填满”，不会出现大量样本聚在均值附近、尾部稀疏的情况。
两种方法对比下来：
| 对比项 | 蒙特卡洛抽样 | 拉丁超立方采样 |
|---|---|---|
| 基本思路 | 完全随机抽样 | 分层均匀抽样 |
| 覆盖均匀性 | 依赖样本量，尾部容易稀疏 | 天然覆盖整个分布区间 |
| 相同样本量下的精度 | 较低 | 较高，方差缩小30%-50% |
| 实现复杂度 | 极低 | 低，核心是分层和逆变换 |
| 适用场景 | 快速验证 | 正式场景生成 |
LHS实现起来也就十几行代码，先对每个随机变量生成分层均匀数，再用wblinv、norminv这类逆累积分布函数映射回物理量。这里有一个细节：如果多个随机变量独立抽样，分层会失去意义，因为它们的排列是随机的。这个时候需要对分层序列做一次随机排序（shuffle），确保风速、光照、负荷之间的相关性是随机的，而不是被分层结构锁死。
2.3 Matlab代码实现与时间相关性处理
下面是我在Matlab里常用的场景生成代码骨架，基于LHS从风速分布生成风速场景，再映射为风电出力：
% 场景生成：拉丁超立方采样 + 逆变换 rng(42); % 固定随机种子，便于复现 Nscenarios = 2000; % 生成场景数 Nperiods = 24; % 调度周期（小时） k = 2.0; % 威布尔形状参数 c = 7.5; % 威布尔尺度参数 % 生成LHS分层样本：每列是分层后的均匀随机数 U = lhsdesign(Nscenarios, Nperiods); WindSpeed = wblinv(U, c, k); % 从威布尔分布逆变换得到风速 % 通过功率曲线映射为出力（简化线性模型） v_in = 3; v_rate = 12; P_rate = 1.0; % 切入、额定风速与额定功率 WindPower = zeros(size(WindSpeed)); WindPower(WindSpeed < v_in) = 0; idx = (WindSpeed >= v_in) & (WindSpeed < v_rate); WindPower(idx) = P_rate * (WindSpeed(idx) - v_in) / (v_rate - v_in); WindPower(WindSpeed >= v_rate) = P_rate;这里lhsdesign是Matlab自带的函数，默认返回0到1之间、按分层方式均匀分布的设计矩阵。wblinv是威布尔逆累积分布函数，把分层均匀数映射为风速值。这套代码跑完，WindPower这个矩阵就是2000行乘24列的出力场景集合，每个场景代表24小时内风电出力的一条完整时序轨迹。
这里有一个我在实际项目中反复强调的点：如果你直接把24个时段独立抽样，生成的风速序列会非常“毛躁”——前一小时风速15m/s，下一小时直接掉到3m/s，这在物理上几乎不可能发生。风电出力有显著的时间连续性，必须引入时间相关性。最简单的做法是给风速序列加一个自回归AR(1)模型：
% 对每个场景施加AR(1)时间相关性 phi = 0.85; % 自回归系数，代表风速的惯性 for s = 1:Nscenarios for t = 2:Nperiods WindSpeed(s,t) = phi * WindSpeed(s,t-1) + (1-phi) * WindSpeed(s,t); end end自回归系数phi的取值决定了出力时序的平滑程度。phi越高，相邻时段的风速越接近，时序曲线越平滑；phi越低，序列越跳变。实测中，风速自相关系数通常在0.8到0.95之间，我一般取0.85左右。注意AR(1)处理后风速的分布会略微偏离威布尔分布——这是一个取舍，工程上这个偏差可以接受，因为物理上的时间相关性比精确的分布形态更重要。
如果不做时间相关性处理，后面的场景削减就会出问题：聚类结果会得到很多“锯齿形”场景，削减出来的典型场景看起来极不自然，调度结果也没有实际参考意义。这个坑我踩过，第一次做出来的场景集每个时段都是独立抽样，削减结果里的“典型场景”在现实中根本不会出现。
光伏出力场景的生成逻辑完全一样，只是把威布尔分布换成Beta分布，然后乘上光照转出力的效率系数即可。负荷场景则直接用正态分布加AR(1)相关性，这里不再赘述。
3. 场景削减实操：把2000个场景压到5个
3.1 削减动机：求解复杂度指数级上升
场景生成出来以后，直接放进优化模型是不现实的。假设你有2000个场景，每个场景24个时段，加上机组启停的整数变量，这个优化问题的规模会急剧膨胀。以一个含5台机组、24小时、2000场景的动态经济调度问题为例，光功率平衡约束就有24×2000=48000条，加上各机组的爬坡约束、启停约束，变量数量轻松突破几十万，其中还有大量二进制变量。直接拿分支定界法求解，计算时间是以小时计的，在实际的日前调度场景中根本无法接受。
所以场景削减的目标很明确：在尽量保留原始场景集概率分布信息的前提下，把2000个场景压缩到5-10个具有代表性的典型场景，并给每个场景重新赋一个概率值。学术界对这个问题有明确的衡量标准——削减前后的两个概率分布之间的Wasserstein距离（也称作推土机距离）最小。通俗地说，就是让削减后的场景集“长得”和原始场景集“尽量像”，概率分布形态不要出现大的畸变。
3.2 两种主流削减方法：聚类与快前向选择
实践中用得最多的两类方法是聚类法和快速前向选择法。
聚类法很简单：把N个场景看作N维空间中的N个点，其中每个点是一个24维向量（代表一整天24小时的风电出力序列），然后用k-means或k-medoids算法把它们聚成K类，每类的质心就是一个典型场景。这一步看似就是数据挖掘里的标准操作，但有几个电力系统特有的坑：一是场景向量的维度是24，直接用欧氏距离在多维空间里算，各时段会等权处理，实际上调度更关心出力高峰时段的准确性；二是k-means只考虑几何距离，没有考虑概率分布信息，聚类后每个场景的概率取所属类的场景数量占比，这在场景本身就等概率生成时是对的，但如果场景带有不同的权重，就需要注意。
我用的更多是改进后的k-means变体——在距离度量里对关键时段加权重。比如峰荷时段误差权重设为2，谷荷时段设为1，这样削减出来的场景会更侧重于保存价格敏感时段的精度：
% k-means场景削减：带时段权重的距离 function [reduced, prob] = scenario_reduction_kmeans(scenarios, K, weights) % scenarios: N x T 原始场景矩阵 % K: 削减后场景数 % weights: 1 x T 时段权重 N = size(scenarios, 1); % 对场景加权 W = sqrt(weights); scen_w = scenarios .* W; % 逐元素加权 [idx, C] = kmeans(scen_w, K, 'MaxIter', 500, 'Replicates', 10); reduced = C ./ W; % 还原到原尺度 prob = histcounts(idx, 1:K+1)' / N; % 场景概率按频次统计 end快前向选择法的思路完全不同，它是从原始场景集合中“挑”场景，每轮挑选一个使原分布与当前选中集合分布之间距离下降最多的场景加入集合，直到选满K个。这种方法不产生“合成场景”，每个保留场景都是真实历史场景，物理上更可信，但计算量比k-means大。对于2000个场景削减到5个，快前向选择在Matlab里跑一轮大约需要数秒到数十秒，对离线场景分析完全没有压力，但对在线滚动优化则偏慢。
两种方法的取舍，我的建议是：做离线场景库分析用快前向选择，物理保真度高；做在线滚动优化用k-means，速度快、结果足够稳定。这里贴一段快前向选择的简化实现，核心逻辑是维护一个“场景距离矩阵”，每轮挑选使总体距离最小的候选场景：
function [reduced, prob] = scenario_reduction_fastforward(scen, K) % 基于Wasserstein距离的快速前向选择 % 计算场景两两之间的欧氏距离 N = size(scen, 1); D = pdist2(scen, scen); selected = []; remaining = 1:N; for k = 1:K best_score = inf; best_idx = -1; for i = remaining % 第i个场景加入后，其他剩余场景到选中集合的最短距离 temp_set = [selected, i]; dist_to_set = min(D(remaining, temp_set), [], 2); score = sum(dist_to_set .* scen_prob(remaining)); % 概率加权 if score < best_score best_score = score; best_idx = i; end end selected = [selected, best_idx]; remaining(remaining == best_idx) = []; end reduced = scen(selected, :); % 概率需要重新分配，参见下节 end这段代码是教学级示例，做项目时还需要处理概率重分配和效率优化。但核心思想很清晰：每次挑一个场景加入，目标是让剩余场景集合到已选场景集合的总距离最小。
3.3 削减后的场景概率：更新与归一化
3.3.1 概率更新规则
场景削减之后，概率的重新赋权比削减本身更容易出错。原始2000个场景如果是等概率采样，每个概率是1/2000；但削减后剩下的K个场景不能简单地把概率设为1/K，因为每个保留场景“代表”的原始场景数量不同。
我用的标准做法是这样的：对每个原始场景，找它距离最近的保留场景，然后把这个原始场景的概率累加到该保留场景上。用公式表达就是：保留场景j的概率 = 所有被分配到以j为最近保留场景的原始场景概率之和。在k-means里这恰好等于聚类成员数除以总场景数；在快前向选择里则需要按最近邻规则重新分配。
% 概率重分配：基于最近邻 function prob = assign_probability(original_scen, reduced) N = size(original_scen, 1); K = size(reduced, 1); D = pdist2(original_scen, reduced); [~, assign] = min(D, [], 2); prob = accumarray(assign, 1/N, [K 1]); end3.3.2 归一化细节
概率分配完之后必须检查一件事：所有保留场景的概率之和是否等于1。在某些实现中，由于初始场景不是等概率的（比如用重要性采样生成），或者削减算法在迭代中有微小的舍入误差，概率之和可能不等于1。
处理方式很简单，在把概率送进优化模型之前做一次归一化：
prob = prob / sum(prob);看上去这是个小到不能再小的操作，但在实际求解随机优化时，概率不归一化会导致目标函数值系统性偏大或偏小——比如所有概率都乘了0.99，最终费用期望值就凭空缩小1%，你拿这个结果跟确定性调度的结果做对比，得出“随机优化更省钱”的错误结论。这种错误在文献里没见过，但在项目复现里非常容易踩到。
3.4 削减质量评价：怎么知道削减得好不好
场景削减的质量评价是个容易被忽略的环节。很多同学做完削减，看一眼削减出来的曲线还算平滑，就直接扔进优化模型了。实际上，评价削减质量有一个量化指标：削减前后经验分布的Wasserstein距离。具体做法是，画出原始2000个场景在某个时段（比如正午12点）出力的经验累积分布函数，然后再画削减后5个场景在同一时段的累积分布函数概率加权结果，两条曲线之间的面积差就是分布误差。
实际操作中，我会分别对比峰荷时段、腰荷时段、谷荷时段的累积分布曲线。如果某个时段的分布偏差特别大，说明这个时段在削减过程中被当成了“不重要的距离”给抹掉了，需要在距离权重里提高这个时段的权重，重新做削减。
还有一个经验法则：削减后的典型场景数量K，不是越大越好，也不是固定的5个或10个。K的大小取决于你的优化模型可接受的求解时间与精度权衡。一般做24小时日前经济调度，K取5-10就够了；如果做的是更精细的日内滚动调度，K取3-5即可，因为滚动窗口会频繁重新求解，新信息会不断修正不确定性。
4. 把削减后的场景落到随机优化模型里
4.1 两阶段随机优化的目标函数与约束
场景削减出来的5个典型场景，每个都带一个概率值。接下来把这套场景集放进微网经济调度模型里。标准的模型结构是两阶段随机规划：第一阶段在“看到”实际出力场景之前做决策，对应于机组启停、备用安排等需要在日前确定的物理决策；第二阶段在第一阶段决策已知的前提下，针对每个具体场景做出经济调度决策，比如各机组的实际出力、与主网的交换功率、负荷削减量等。
以一台柴油机、一台燃气轮机、一组储能电池、光伏风电并网、与主网可交换功率的典型微网为例，目标函数写成：
min sum_t (启动成本 + 运行成本) + sum_s prob_s * sum_t (场景s下购电成本 + 弃风惩罚 + 失负荷惩罚)第一项是所有开机时段的基础启停成本和空载成本，与场景无关；第二项是各场景下运行费用的期望值。弃风惩罚和失负荷惩罚是“软约束”手段——给弃风和切负荷一个很高但不至于离谱的惩罚单价（通常是电价的几倍到几十倍），让优化器在遇到极端场景时宁可少量弃风或切负荷，也不违反物理硬约束。这种方法在工程实施中很有效，因为完全靠硬约束去保证所有场景下的功率平衡，往往会导致可行域为空。
约束条件包括：各场景下的有功功率平衡约束、机组出力上下限约束、爬坡约束、储能荷电状态递推约束与容量约束、联络线功率约束、系统备用约束等。其中功率平衡约束要按场景展开：
对每个时段 t、每个场景 s： 柴油机出力 + 燃机出力 + 风电出力(s) + 光伏出力(s) + 储能放电 - 储能充电 + 主网购电 - 主网售电 - 负荷(s) - 失负荷(s) + 弃风(s) = 0注意这里风电和光伏出力是场景给出的已知参数，柴发和燃机出力是决策变量，储能充放电和主网交换也是决策变量。弃风和失负荷是非负的松弛变量，仅在极端场景下有非零值。
备用约束通常是所有场景共用的，表示如果某个场景的出力低于预期，系统的向上/向下可调能力要足够。这个约束在随机规划里是耦合约束——它把所有场景的系统状态联结在一起，增加了求解难度，但却是保证鲁棒性的关键。
4.2 场景编号法与决策变量组织
把场景和时段耦合进模型时，一个常见的建模技巧是“场景编号法”：先把所有决策变量加上场景下标和时段下标，然后按“大M矩阵”一次性把所有约束写进Matlab的优化问题结构里。以Yalmip工具箱为例，决策变量定义方式如下：
% 决策变量定义 sdpvar Pg(Ngen, Nperiods); % 第一阶段：机组出力基准 binvar u(Ngen, Nperiods); % 第一阶段：机组启停状态 sdpvar Pg_s(Ngen, Nperiods, K); % 第二阶段：各场景下机组出力 sdpvar Pbuy_s(Nperiods, K); % 第二阶段：向主网购电 sdpvar Psell_s(Nperiods, K); % 第二阶段：向主网售电 sdpvar shed_s(Nperiods, K); % 第二阶段：失负荷 sdpvar curtail_s(Nperiods, K); % 第二阶段：弃风关键的设计决策是：机组启停变量u不带场景下标，因为启停决策在日前已经确定；而实际出力变量带场景下标，它们是第二阶段根据具体场景调整的变量。这样设定保证了“两个阶段”的逻辑一致性——启动成本按确定的启停方案发生，而运行成本和购电成本按每个场景的概率加权。这个设定和实际情况一致：调度员在日前决定哪些机组开、哪些机组关，日内根据实际风光出力决定每台机组具体发多少。
约束条件的场景展开在Yalmip里用一个双层循环即可：
Constraints = []; for t = 1:Nperiods for s = 1:K Constraints = [Constraints, sum(Pg_s(:,t,s)) + Pwind(s,t) + Ppv(s,t) + Pbuy_s(t,s) ... - Psell_s(t,s) + Pdis_s(t,s) - Pch_s(t,s) + shed_s(t,s) ... - curtail_s(t,s) == Pload(s,t)]; end end4.3 求解器选择与Yalmip建模要点
求解随机优化模型，Matlab原生自带的linprog只支持线性规划，处理混合整数线性规划要用intlinprog。如果模型里目标函数是线性的，约束也是线性的（包括二进制启停变量），intlinprog就能解决中小规模问题。但一旦场景数超过10个，机组数超过5台，整数变量超过数百个，intlinprog的求解速度会明显下降。这种情况下我会用Gurobi或CPLEX这类商用求解器，它们在处理大规模混合整数规划时性能要好得多，配合Yalmip接口使用体验很顺滑。
下面是我在Matlab中用Yalmip构建随机优化模型的标准结构：
% 目标函数：启动成本 + 期望运行成本 Objective = sum(sum(StartCost .* u)) + sum(sum(FuelCost .* Pg)); for s = 1:K Objective = Objective + prob(s) * (sum(Pbuy_cost(:,s)) + ... penalty_shed * sum(sum(shed_s(:,:,s))) + ... penalty_curtail * sum(sum(curtail_s(:,:,s)))); end % 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 1, 'mipgap', 0.001); optimize(Constraints, Objective, ops);这里有几个关键参数要说明一下。mipgap是混合整数规划的相对对偶间隙，设置成0.001表示求解器找到的解与理论最优解的差距在0.1%以内就停止。这个参数很有用，因为微网经济调度问题的规模决定了对偶间隙降到很低很花时间，而0.1%的精度在工程上绰绰有余，可以大幅缩短求解时间。verbose控制求解器的输出详细程度，设置为1能让你看到每次迭代的间隙变化和上下界，对调试很有帮助。
还有一个非常实用的小技巧：把约束写成矩阵形式而不是循环展开形式，求解速度会快很多。Yalmip内置的repmat、kron函数可以帮助做约束展开。比如备用约束在所有时段所有场景都成立，完全可以用维度广播的方式一次性写入，避免几十个循环嵌套拖慢建模速度。
5. 我踩过的坑：场景法的常见问题与排查
5.1 场景集失真导致调度结果过于乐观
我在做第一个微网随机调度项目时，犯过一个很典型的错误：场景生成时只考虑了风速的威布尔分布，没有考虑风机出力的截断效应。威布尔分布可以生成任意大的风速值，但风机出力在额定风速以上就封顶了，在切入风速以下就是零。生成出来的风速场景经过功率曲线映射后，许多场景的出力被“压”在0和额定值两个点上，导致风电出力的实际分布呈双峰形态，而不是光滑的连续分布。这个失真信息直接传递到优化模型里，调度结果高估了风电的调节能力，某些时段的备用容量严重不足。
排查思路很简单：生成场景之后先画分布直方图，对比原始风速分布映射后的理论出力分布。如果两个分布形状差异很大，找原因。这类问题在文献中不太显眼，但在实际数据面前表露无疑。
5.2 k-means聚类结果的随机性问题
k-means算法的初始聚类中心是随机选择的，意味着同一份场景数据跑两次，削减结果可能不同。如果你在写论文或者做项目报告，这个问题是致命的——评审人要求你的结果可复现，而代码两次跑出来的典型场景不一样，数据对不上。
解决方法有两个层面。第一，固定随机种子，Matlab里用rng(42)这样固定种子后再跑k-means；第二，更稳健的做法是设置Replicates参数，让k-means从多组不同的初始中心出发，最终取目标函数值最小的一组结果。我通常设置Replicates = 20，这样聚类结果的随机性基本可以忽略。在科研项目中，我还会把固定种子的值写进代码注释和报告的方法学部分，确保任何人重新运行都得到一致结果。
5.3 场景数K的选择：不是越大越好
场景削减数K的确定，没有放之四海而皆准的公式。我见过有人把K设成50、100，理由仅仅是“怕丢失信息”，结果模型规模爆炸，求解时间从几分钟涨到几小时。我也见过K设成3，结果调度结果过于保守，运行成本比确定性模型还高，完全失去了随机优化的意义。
我用过一个实用的经验法则：K的取值应当让削减前后每个典型时段出力分布的Wasserstein距离缩小到某个可接受的阈值以下。具体操作时，我通常从K=3开始，以2为步长逐步增加，每次削减后计算削减前后的分布距离误差，画一条“K-误差”曲线。曲线明显拐点的位置就是合理的K取值。对于24时段微网调度，这个拐点通常出现在K=5到K=10之间。
5.4 别忘了做场景生成与削减的联动校验
最后一个忠告：场景生成和场景削减永远不是两个孤立的步骤。削减算法对原始场景集的特征非常敏感，原始场景集如果本身相关性太弱（比如24个时段独立抽样），削减结果大概率是一些在物理上没有意义的“伪典型场景”。正确流程是先做完时间相关性处理，检查生成场景的统计特征（均值、方差、自相关系数），确认合理后再进入削减。
我自己的项目流程中，场景生成后一定会先保存一份“统计特征自检表”，包括：每个时段出力的均值、标准差与历史数据的偏差、相邻时段出力的自相关系数、极端场景（出力接近0或额定值）的占比。这份自检表会伴随整个项目交付，作为场景质量的可审计记录。很多项目做完后报告写不出来，就是因为缺了这一步——没有这份统计记录，你拿什么证明你的场景库能代表真实的不确定性？
把场景生成、场景削减和随机优化串在一起后，这套方法还可以往多个方向扩展。比如把场景生成里的静态分布换成时序场景树，用于多阶段决策；把削减算法换成基于深度生成的对抗网络来生成高维场景；或者在优化模型里引入条件风险价值（CVaR）来控制尾部风险。这些扩展都是在同一个框架上的演进，核心逻辑不变——先用概率描述不确定性，再压缩成可计算的场景，最后在优化模型里显式地、而不是拍脑袋地处理风险。我在实际项目中验证过的最有效路径，是先跑通5个场景的小模型，确认整个链路正确、结果合理，再去扩展场景数和模型复杂度。直接上大规模模型，出了问题连排查都无从下手。
