企业官网建设流程全解析

做电力系统随机优化的人，大概率都经历过这样的场景：拿到一批风光出力历史数据，本来想先用蒙特卡洛（MC）生成未来可能会出现的出力路径，再扔进优化模型里去算期望目标、风险约束，结果生成出来的场景经常“不像人样”——上午10点出力接近额定值，11点突然掉到零，12点又拉满，真实的新能源出力哪会这么跳？优化器面对这种伪爬坡，要么给出过度保守的调度方案，要么干脆收敛不下去。问题往往不在MC本身的采样效率上，而在生成场景时没有把时序相关性放进去。

这篇文章想系统梳理一下“考虑时序相关性的MC场景生成与削减”这条技术路线：时序相关性为什么重要、怎么嵌进MC采样、场景削减怎么做才不伤信息，以及完整跑通一遍需要注意哪些坑。内容主要面向正在做新能源出力场景建模、随机优化、时序模拟的研究生和一线工程师，也适合刚接触随机规划但被场景问题困扰的人参考。

1. 先想清楚为什么场景生成必须“看时序”

1.1 随机优化的第一步：把不确定性压缩成场景

随机规划处理不确定性时，通常不是直接面对连续概率分布，而是把分布近似成一组离散场景。每个场景是一条完整的时间序列，比如未来24小时的风电出力曲线，并带一个概率权重。随机期望值模型就是在这些场景上做期望目标优化，相关约束也都要在场景层面满足。

这个思路的核心假设是：只要场景集合能代表原始不确定性，优化结果就有统计意义。以前我在某风电并网项目中做日前调度，最大的体会就是场景质量直接决定优化结果的可用性。场景选得烂，后面调度模型再精巧也救不回来。所以场景生成不是预处理里可有可无的一步，而是整个随机优化链条的地基。

1.2 朴素蒙特卡洛生成法的两个硬伤

最基本的MC做法，是对每个时段单独拟合分布，然后独立采样。假设未来有T个时段，每个时段按自己的均值、方差、分布类型分别抽样，最后拼成一条场景。代码写起来确实简单：

for t in range(T): scene[t] = dist[t].rvs()

但这么做有两个明显的硬伤。

第一个是忽略时段之间的相关性。相邻时刻出力被当成互不相关的独立变量，生成序列的相邻时段相关系数会非常低。真实的风电出力、光伏出力、负荷曲线都带有明显的物理惯性，不会无规则乱跳，独立采样等于把物理规律丢了。

第二个硬伤是序列震荡太剧烈。由于每个时段独立抽样，场景里很容易出现“高—低—高—低”的锯齿状波动。当年我第一次跑风电场景生成，算出来的场景平均每小时爬坡率比真实历史数据高出几倍，调度同事看了一眼就说不能用。这个教训让我后面再也不敢轻视时序结构。

1.3 时序相关性：说你强，你到底强在哪

“时序相关性”说白了就是物理惯性：风电出力不会在短时间内凭空消失又凭空出现，光伏出力从上午到下午是平滑爬升再平滑下降的，负荷也有明显的时段连续性。数学上，这些东西可以用相邻时段协方差、相关系数矩阵或者自相关函数来刻画。

举个例子，某实际风电项目历史出力数据里，相邻1小时间段的相关系数通常在0.85到0.95之间；如果完全不处理相关性，用独立MC采样生成的场景，相邻时段相关系数可能只有0.3甚至更低。差了这么多，生成场景的时序形态当然对不上真实数据。

所以场景生成不能只考虑单点边际分布，必须把时序结构一起放进采样过程。这也是本文题目里“考虑时序相关性”这几个字的来头。

2. 生成阶段怎么把时序相关性嵌进MC采样

2.1 先搞懂相关系数矩阵怎么进采样

时序相关性的核心载体是相关系数矩阵C，维度是T×T，第i行第j列表示时段i与时段j的相关系数。想从一组独立标准正态样本Z出发，得到协方差结构为C的样本X，最常用的工具是Cholesky分解。

假设C = L·Lᵀ，其中L是下三角矩阵。令：

X = Z·Lᵀ

得到的X就是均值为0、协方差矩阵为C的样本。这个思路本质是线性变换，保证样本的二阶矩结构符合要求。但这里有个前提：C必须是对称半正定矩阵，否则Cholesky分解直接报错。我第一次跑这段代码时就被非正定矩阵天天折磨，后面第5章会专门讲怎么处理。

2.2 正态分布下的Cholesky方案

假设每个时段的出力近似服从正态分布，已知均值向量μ、标准差向量σ和时序相关系数矩阵C，那么生成N个场景的代码很简单：

import numpy as np def gen_corr_gaussian_scenes(n_scenes, n_periods, mu, std, corr_matrix): L = np.linalg.cholesky(corr_matrix) Z = np.random.normal(size=(n_scenes, n_periods)) X = Z @ L.T # X的协方差矩阵等于corr_matrix scenes = mu + std * X # 叠加各时段均值与标准差 return scenes

这里有个容易混淆的细节：corr_matrix是相关系数矩阵，不是协方差矩阵。如果你直接用协方差矩阵做分解，就不需要再乘std；但如果手里是相关系数矩阵，必须记得在最后乘上各时段的std，才能得到正确量纲的出力场景。

2.3 Nataf变换：让边际分布不再限制你

风电出力通常不是正态的，用Beta分布或者Weibull分布拟合更常见；光伏出力也很少是严格正态。这时如果强行用正态假设，边际分布的范围和尾部特征都会失真。解决办法是走Nataf变换路线。

Nataf变换的步骤是：先用Cholesky生成相关标准正态样本X，再用标准正态累积分布函数Φ(X)把X转成均匀分布样本U，最后用目标分布的分位数函数F⁻¹(U)映射回真实出力。代码逻辑是这样：

from scipy.stats import norm, beta # U = norm.cdf(X) # 标准正态 -> 均匀 # scenes = beta.ppf(U, a, b, loc, scale) # 均匀 -> 目标Beta分布

注意一个关键问题：原始空间的相关系数不能直接照搬到标准正态空间。因为非线性变换会改变相关性。严格做法是用Nataf公式或经验修正系数对相关系数矩阵做一次修正。不过我在工程项目里对比过，T的维度几十个时段时，直接用近似修正，生成场景的误差完全在接受范围内。如果追求严谨，可以先做一次简单的数值模拟校验，看看近似相关矩阵和真实采样相关矩阵差多少。

2.4 Copula与马尔可夫链：两条值得备用的路线

除了Cholesky和Nataf，业界还有两条常见路线。

Copula的核心思想是把边际分布和相关结构拆开。边际分布单独建模，相关结构用Gaussian Copula或t-Copula描述。做法是先生成相关的均匀分布样本，再通过各变量的分位数函数变换回目标分布。优点是理论框架干净、变量之间的相关结构可以显式控制，特别适合风电场和光伏电站联合出力这种“多变量+时序”场景。缺点是实现比Cholesky复杂一点，尤其要选Copula类型和估计参数。

马尔可夫链则适合把出力状态离散化的情况。比如把风电出力分成低、中、高三个状态，用转移概率矩阵控制时序变化。它能刻画非高斯、非线性的状态转换，但状态数量一多，矩阵维度会爆炸，高维场景下用起来很麻烦。

我个人的选择经验是：数据量中等、相关性“规矩”的工程场景，优先Cholesky或Cholesky+Nataf；要在论文里做严格概率建模，再用Copula；马尔可夫链更适合做离散状态下的解析推导，不是场景生成的首要选择。

3. 场景削减：把上千个场景压缩成可用的十几个

3.1 削减问题的数学本质

MC场景动辄几千个，直接丢进优化模型里，尤其模型里带整数变量时，求解器很容易卡死。场景削减的目标，就是在保持概率分布特性尽量不变的前提下，把N个场景压缩到n个，n远小于N，但n个场景的概率加权分布要和原始N个场景的分布尽量接近。

业界最常用的数学框架，是让削减前后概率测度之间的Kantorovich距离（也叫Wasserstein距离，一阶）最小。直觉上就是让削减后的场景集合尽量“贴近”原始场景集合。这个距离既衡量了场景之间距离，又把场景概率权重考虑了进去，比单纯聚类要好解释。

3.2 快速前向选择与同步回代消除

场景削减常用的两种经典启发式算法，一个是同步回代消除，一个是快速前向选择。

同步回代消除（BS）的思路是：反复找一个“最不重要”的场景删掉，并把它的概率合并到最近的保留场景上，一直删到目标数量为止。

快速前向选择（FFS）则是反着来：从空集开始，反复选一个“最能代表”剩余集合的场景加入保留集，一直加到目标数量为止。

我贴一段同步回代消除的关键实现，比较直观：

def backward_reduction(scenes, probs, target): # scenes: (n, T) 原始场景集合 # probs: (n,) 场景概率 remain = list(range(len(scenes))) p = probs.copy() while len(remain) > target: R = np.array(remain) # 计算剩余场景两两之间的距离 D = np.linalg.norm(scenes[R][:, None, :] - scenes[R][None, :, :], axis=-1) np.fill_diagonal(D, np.inf) # 屏蔽自身距离，避免删了自己又合并到自己 # 每个场景被删除造成的概率距离代价 = 概率 * 到最近邻的距离 min_dist = D.min(axis=1) cost = p[R] * min_dist # 删除代价最小的场景，把概率累加到最近保留场景上 rm_idx = np.argmin(cost) nearest_idx = np.argmin(D[rm_idx]) rm = R[rm_idx] near = R[nearest_idx] p[near] += p[rm] remain.remove(rm) return scenes[remain], p[remain]

注意里面两个细节：一是距离矩阵的对角线要设成无穷大，否则每个场景的最近邻都是自己，删了等于没删；二是概率合并一定要用原始场景数组的索引，别在删除过程中错位。这个实现时间复杂度偏高，但几千个场景完全跑得动，够用。

3.3 聚类削减：另一种直观解法

除了SBR（基于概率距离的削减），聚类削减也很常见。做法是把每个场景当成T维空间里的一个点，跑K-means或K-medoids，把聚类中心当作代表场景，聚类内所有场景的概率之和就是这个代表场景的概率。

聚类削减的优点是用现成库就能做，非常快。缺点也很明显：K-means用的是欧氏距离，对时序错位很敏感。两条波形相同但相位错一个小时的场景，欧氏距离会非常大，聚类中心也会把两条场景的形态“平均”得面目全非。如果要用聚类削减，可以试试用动态时间规整（DTW）距离替代欧氏距离做K-medoids，效果会好不少，但计算量也上去了。

我自己常用的策略是：如果场景量特别大（上万），先用聚类粗筛到几百个，再用SBR精削到目标数量，速度和效果都觉得可以。

3.4 怎么判断削减是否“削坏了”

削减效果不能只看削减后的几条曲线好不好看。我做场景削减，必看的指标有三个：

第一，削减前后期望出力曲线的相对误差，一般要控制在3%以内。第二，极端分位数误差，比如 5% 和 95% 分位点出力值，削减不应该把极端场景都洗掉。第三，削减前后概率距离本身，保持在一个可接受范围。

经验值上，把5000个场景削减到20个到50个，期望误差控制在3%以内常有的事；如果削减到10个以下，极端场景丢失的风险会明显上升。建议做一条“削减数量-误差”曲线，看曲线什么时候拐平，再决定保留多少个场景。

4. 完整实操：从历史出力数据到最终场景集

4.1 数据清洗与边际分布拟合

第一步，先把历史出力数据按时间尺度整理成矩阵。每一行是一条历史出力序列，每一列是一个时段。比如做日前调度，时段取24，每条历史序列就是一天24个点的出力；如果做滚动优化，可以按15分钟间隔取96个点。

数据清洗要处理缺失值、异常值和跳变点。我遇到最多的异常是出力超过装机容量的记录、连续零点突然冒出的尖峰，这些多半来自数据采集错误，直接剔除或用前后时段插值替换就好。

清洗完之后做边际分布拟合。先画直方图和QQ图，判断数据大概服从什么分布。风电功率在接近0和接近额定值的位置概率密度通常很高，Beta分布拟合效果普遍不错；有些地区数据也可能更接近Weibull分布。不要迷信某一种分布，对每个时段分别估计参数，对比一下拟合优度再选。

4.2 时序相关性矩阵估计与修正

边际分布搞定后，计算任意两个时段间的样本相关系数，得到T×T时序相关矩阵。理论上的相关系数矩阵必须对称半正定，但样本估计出来的矩阵经常不满足，尤其在时段数多、历史样本少的情况下，总会出现一两个负特征值，直接导致Cholesky分解失败。

修正方法我常用特征值裁剪：对相关矩阵做特征值分解，把所有负特征值置为0（或者一个很小的正数），用剩余特征值重构矩阵，然后重新归一化，让对角线元素变成1。简单实现如下：

def fix_corr_matrix(C): vals, vecs = np.linalg.eigh(C) vals = np.clip(vals, a_min=1e-8, a_max=None) D = np.diag(vals) C_fixed = vecs @ D @ vecs.T # 归一化对角线为1 inv_sd = 1.0 / np.sqrt(np.diag(C_fixed)) C_fixed = C_fixed * np.outer(inv_sd, inv_sd) return C_fixed

这个方法效果很稳，我也试过直接缩小非对角元素来强行走通Cholesky，但相关性保留度会损失很大，不推荐。

4.3 生成带时序相关性的MC场景（代码示例）

下面给一段基于Cholesky+Nataf思路的完整示例。这里假设选择Beta分布作为边际分布，先生成相关标准正态样本，再转换到Beta边际：

import numpy as np from scipy.stats import norm, beta def gen_corr_beta_scenes(n_scenes, n_periods, beta_params, corr_matrix): # beta_params: list of (a_t, b_t, loc_t, scale_t)，每个时段一组参数 L = np.linalg.cholesky(fix_corr_matrix(corr_matrix)) Z = np.random.normal(size=(n_scenes, n_periods)) X = Z @ L.T # 相关标准正态样本 U = norm.cdf(X) # 等概率变换到均匀分布 scenes = np.zeros_like(U) for t in range(n_periods): a, b, loc, scale = beta_params[t] scenes[:, t] = beta.ppf(U[:, t], a, b, loc, scale) return scenes

这里把相关性修正函数直接挂进去了，避免矩阵非正定导致程序中断。用Beta分布做边际，出力会被限制在0到装机容量范围内，更贴近实际。生成5000个场景之后，下一步就是削减。

4.4 场景削减落地（代码示例）

沿用第3章的backward_reduction，把5000个场景削减到20个，代码就一行：

scenes_red, probs_red = backward_reduction(scenes_all, probs_all, target=20)

但target选多少，不要拍脑袋。我的做法是先分别削减到5、10、20、30、50个场景，跑一遍下游优化模型，画一条目标函数随场景数变化的曲线。如果20个场景和50个场景算出来的目标函数差不了多少，那就取20个，省计算量；如果误差还是很大，继续往上加。

概率合并完成之后，检查一下：

print("sum probs:", probs_red.sum())

正常情况下必须等于1。如果概率和不对，说明合并概率时漏了场景，回头排查删除索引。

4.5 削减结果验证

削减完成后，别急着进优化模型，先做几个快速验证。

第一个是可视化对比。把原始5000个场景的包络带画出来，也就是每个时段取5%到95%分位数，再叠加削减后的20条场景曲线。削减后的场景应该基本落在包络带内部，并且分布比较均匀，不能全挤在上沿或下沿。

第二个是统计指标对比。计算削减前后期望出力曲线和标准差曲线，误差控制在可接受范围。第三个是自相关系数对比。削减后场景集合的相邻时段相关性，应该和原始场景集以及历史数据比较接近。这个指标常被忽略，但恰好在验证时序相关性的环节特别重要。

5. 常见问题与排查技巧实录

5.1 相关矩阵非正定，一运行就报错

这是最常遇到的坑，几乎每个做Cholesky方案的人都会撞上。原因无非两种：样本相关矩阵估计噪声导致出现负特征值，或者时段数T大于历史样本数N。

处理办法前面已经给了：特征值裁剪加对角线归一化。如果修完矩阵之后相关性损失太大，可以试试增加历史样本量或者对相关矩阵做收缩估计，把非对角元素往0方向收缩一点，整体数值稳定性会提升。

5.2 削减后分布指标不对

削减之后发现期望误差很大，大概率不是算法本身的问题，而是操作细节出错。重点排查三件事：

一是距离度量是不是用了原始量纲。不同时段出力数值跨度大，建议先对场景做标准化处理再算距离，否则某些时段的距离会被放大。二是概率合并时有没有漏项。同步回代里每次删除一个场景都要把概率累加到最近保留场景上，漏一步概率和就不对，后面指标自然全错。三是有没有做归一化。削完后概率总和重新检查一次，归一化到1。

5.3 削减数量怎么选

没有绝对标准，但有一个足够实用的经验法。

场景数量少，计算快，但会丢掉极端场景；场景数量多，分布保真度好，但优化模型会变慢。实践中可以按比例试：5000个原始场景，分别削减到10%、5%、1%的数量，看目标函数变化。小规模测试场景取10到20个就够，出论文或者做精细调度，取30到50个更稳妥。

我在某风电场项目中做过对比，5000个场景削减到20个，目标函数值和用5000个场景算出来的几乎一致，但求解时间从两小时降到了十分钟以内，性价比极高。

5.4 时序相关性处理过度的表现

相关性处理过头，也会出问题。

表现就是所有场景几乎平行，相邻时段变化太“顺滑”，场景之间高度相似，集合整体方差被压缩，极端情况覆盖不足。原因多半是把相关性全塞进线性变换后，随机样本的自由度被相关结构吃掉太多。

解决办法是在生成时保留一部分独立扰动项，也就是在相关结构样本基础上叠加一个小的独立噪声项。这样既能保住时序连续性，又不会让所有场景都长得差不多。检查方法很直接：对比削减后场景集合的5%分位带与原始历史数据的5%分位带，如果新场景带太窄，就要把独立扰动调大一点。

我个人在实际操作中的体会是，MC加Cholesky（或者Cholesky+Nataf）生成，配合同步回代消除完成削减，是一条性价比极高的经典路线。风电、光伏这类出力随机性“有规律可循”的场景，完全够用，没必要一上来就上复杂的深度生成模型。尤其当你需要生成结果可解释、可复现、能写进论文或工程报告时，这条路线的优势会非常明显。

最后再分享一个小技巧：场景削减完成之后，后续优化里要同时利用“场景抽稀”和“概率加权”，让每个场景标明权重后再代入目标函数，而不要简单地把20个场景做平均后当成一个确定性场景。很多人在这里丢掉概率信息，前面花的工夫就白费了。