1. 这个项目到底在解决什么问题1.1 电力系统研究里的不确定性从哪来做电力系统规划、微电网调度或者配电网分析的工程师和研究生大概率都碰到过这个场景风电出力时大时小光伏白天晴雨不定负荷曲线跟着季节和作息波动偏偏这些不确定性又直接影响调度策略的可行性、备用容量的设置、储能配置的合理性。你要是只用一条典型日曲线去做优化结果往往偏乐观实际运行一兑现就出问题。那怎么把不确定性考虑进去业界主流思路之一就是多场景法把风光和负荷可能出现的各种情况用大量随机生成的场景表示出来再通过缩减算法挑出少量有代表性的典型场景用这几十个场景去逼近原始上千个场景的统计特性。这样一来随机优化问题就转化成了确定性场景集合上的优化问题计算复杂度可控模型也足够贴近真实情况。这也是风光及负荷多场景随机生成与缩减Matlab代码这类项目如此普遍的根源。它本质上是一套不确定性建模的预处理工具箱不管你后续做两阶段鲁棒优化、随机规划还是做基于场景的调度策略评估这套代码都能直接复用。1.2 多场景方法的核心思想把连续问题变成离散问题我用一个简单的类比来解释多场景法的价值。想象你要评估一个城市未来一年的天气对出行的影响理论上天气有无穷多种组合。但实际决策时你不需要这一年每一天的天气你只需要知道晴天、雨天、台风天、高温天各有几成概率每种天气下大家的出行方式如何变化。把无穷种天气聚类成几种典型天气就是场景缩减的基本思想。放到电力系统里风电、光伏、负荷的出力都是连续随机变量理论上它们的联合分布是无穷维的。多场景方法分两步走第一根据风、光、负荷的历史统计规律用蒙特卡洛采样生成大量随机场景比如1000个每个场景包含一组某时刻风电出力光伏出力负荷大小的时序数据第二用场景缩减算法常见的有同步回代缩减、K-means聚类等把1000个场景精简成10~20个典型场景同时给每个典型场景算出一个概率权重。缩减后的场景集合在概率分布意义上尽量接近原始场景集合但规模大幅缩小可以直接嵌入优化模型。这套流程的价值在于它把不确定性从抽象的统计分布变成了优化模型里实实在在可以写进约束的一组组系数和权重。这正是随机规划、机会约束规划、鲁棒优化等高级建模方法的地基。2. 场景生成与缩减的数学原理2.1 风光负荷的概率分布模型怎么选在做场景生成之前第一步是确定每个随机变量的概率分布。这块没有统一标准但工程实践里有一些默认的老规矩。风电出力通常先用两参数Weibull分布描述风速再做风机功率曲线转换。风速概率密度函数为f(v) (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)其中k是形状参数c是尺度参数。但直接用出力数据做分布拟合时很多人也会直接用Beta分布拟合风电出力因为出力范围是0到额定功率之间Beta分布定义在[0,1]区间天然匹配。光伏出力主流做法是用Beta分布拟合光照强度或出力水平同样是因为光伏出力归一化后在0到1之间Beta分布的两个形状参数alpha和beta可以灵活拟合不同偏态。也有不少论文直接用正态分布近似精度略低但操作简单。负荷呢一般假设服从正态分布或对数正态分布均值和方差从历史数据中统计得到。负荷的特性是相关性比较强——同一区域内的不同负荷节点往往同涨同跌所以做多节点负荷场景时还要考虑相关性。![注意] 分布假设并不是越复杂越好。我见过不少论文把分布建模写得很花哨但实际用下来Beta分布拟合光伏、Weibull或Beta拟合风电、正态拟合负荷加上相关性修正对绝大多数规划调度问题已经够用。花哨模型带来的收益微乎其微反而让参数估计和代码调试变得复杂。2.2 蒙特卡洛采样与相关性处理确定了边缘分布之后场景生成的核心就是采样。最简单的做法是直接对每个随机变量独立采样然后组合成场景。但这样会产生一个问题风、光、负荷之间往往存在相关性比如同一区域的风电和光伏可能负相关白天风小、晚上风大不同节点的负荷强正相关。忽略这些相关性生成出来的场景虽然单个变量分布正确但联合分布完全失真后续优化结果也就失去意义。处理相关性的标准做法是Cholesky分解配合Nataf变换。大致步骤是根据历史数据计算各随机变量之间的相关系数矩阵R对R做Cholesky分解得到下三角矩阵L满足R L * L生成一组独立的标准正态随机向量Z做线性变换得到相关标准正态向量Y L * Z用等概率变换把相关标准正态分布映射回各自的目标分布。这个流程在Matlab里实现很顺手核心只是几个矩阵运算。等概率变换的原理是标准正态分布的累计分布函数值服从[0,1]均匀分布再通过目标分布的逆累计分布函数就能把相关性保留下来。需要注意Cholesky分解要求相关系数矩阵是正定的。实际操作中如果矩阵里存在完全线性相关的变量或者采样数量太少导致经验相关矩阵不正定分解会直接报错。我的建议是在生成前对相关矩阵做一次特征值修正把所有负特征值往零方向做极小偏移。2.3 场景缩减同步回代法的核心逻辑缩减这一步是整个流程的灵魂。做缩减的方法很多工程里最常见的是两种基于概率距离的同步回代缩减和K-means聚类缩减。同步回代法的思路非常直观原始有N个场景每个场景有自身的概率权重等概率时就是1/N现在要删掉一部分场景把被删场景的概率累加到保留场景上使得删减前后场景集合的概率分布距离最小。具体步骤是迭代进行的每一轮删除一个场景对当前场景集合中的每一对场景(i, j)计算它们之间的距离d(i,j)常用的是 Kantorovich 距离本质上就是在某范数如2范数下求两个场景时序数据的距离找到距离最近的那一对删掉其中一个把被删场景的概率累加到与之距离最近的那个保留场景上重复上述过程直到场景数达到预设的目标值。这个算法实现起来不复杂但直接暴力循环在初始场景数很大时会很慢。1000个场景缩减到10个理论上最多迭代990轮每轮都要遍历当前所有场景对复杂度近似O(N^3)Matlab纯循环跑起来可能要几十秒甚至几分钟。不过考虑到这是离线预处理一次性计算时间通常可以接受。K-means聚类的思路则是把N个场景看成N个样本点每个样本点是T维向量T是时间段数然后直接聚成K类每个类的质心就是典型场景类内样本占比就是场景概率。K-means实现简单、速度快但缺点也很明显K值要事先指定聚类结果依赖初始中心的选择每次运行可能不同而且K-means用欧氏距离聚类本质上假设每个时段独立同权对峰谷差异明显的数据会偏向平均化。我在后面的代码实现里两种方法都会给出来默认用同步回代法因为它在论文里认可度更高、结果更稳定。3. Matlab代码实现与实操流程3.1 整体架构与关键函数设计整个代码我建议拆成四个模块每个模块一个脚本或函数方便独立调试和复用数据输入模块读入风光负荷的历史数据或者直接设置分布参数对于没有实测数据的研究场景用分布参数合成数据是常见的做法场景生成模块实现蒙特卡洛采样、相关性处理输出初始场景矩阵场景缩减模块实现同步回代缩减或K-means缩减输出典型场景和对应概率评估与可视化模块对比缩减前后场景的均值、方差、累积分布画出场景集和典型场景曲线。模块化的好处是后面你换一套数据只需要改数据输入模块想换缩减算法只替换场景缩减模块其他部分完全不用动。3.2 核心代码实现含参数计算过程我给出一个可以直接运行的完整示例场景设定为两个风电场、一个光伏电站、一个负荷节点一共4个随机变量24小时时段初始生成1000个场景缩减到10个。%% 参数设置 num_scenarios 1000; % 初始场景数 num_typ_scenarios 10; % 缩减后的典型场景数 num_vars 4; % 随机变量数风电2个 光伏1个 负荷1个 num_periods 24; % 时间维度24小时 %% 分布参数设置以Beta分布为例 % 风电1Beta分布参数 a1, b1每个时段一组这里为简洁只给全天均值示例 % 实际使用时应按24个时段分别拟合参数 beta_params_wind1 [2.0, 3.5]; % 形状参数 beta_params_wind2 [1.8, 4.0]; beta_params_pv [3.0, 5.0]; % 光伏出力夜间时段应近似为0 % 负荷正态分布 mu_load 100; % 均值单位MW sigma_load 10; % 标准差单位MW %% 相关系数矩阵4x4 R [1.0, 0.3, -0.2, 0.1; 0.3, 1.0, -0.1, 0.2; -0.2, -0.1, 1.0, 0.0; 0.1, 0.2, 0.0, 1.0]; %% 场景生成 % 步骤1Cholesky分解 % 先做一次特征值修正确保正定 [V, D] eig(R); d diag(D); d(d 1e-6) 1e-6; R_fixed V * diag(d) * V; L chol(R_fixed, lower); % 步骤2生成独立标准正态随机数 Z randn(num_vars, num_periods, num_scenarios); % 步骤3施加相关性每个时段独立处理 Y zeros(size(Z)); for t 1:num_periods Y(:, t, :) L * squeeze(Z(:, t, :)); end % 步骤4等概率变换到目标分布 % 标准正态CDF - 均匀分布 U normcdf(Y); % 均匀分布 - Beta分布风电和光伏 scenarios_wind1 betainv(U(1, :, :), beta_params_wind1(1), beta_params_wind1(2)); scenarios_wind2 betainv(U(2, :, :), beta_params_wind2(1), beta_params_wind2(2)); scenarios_pv betainv(U(3, :, :), beta_params_pv(1), beta_params_pv(2)); % 均匀分布 - 正态分布负荷 scenarios_load norminv(U(4, :, :), mu_load, sigma_load); % 组装成场景矩阵每行一个场景每列对应 (变量, 时段) % 存储为 struct 更清晰 for k 1:num_scenarios scen(k).wind1 squeeze(scenarios_wind1(1, :, k)); scen(k).wind2 squeeze(scenarios_wind2(1, :, k)); scen(k).pv squeeze(scenarios_pv(1, :, k)); scen(k).load squeeze(scenarios_load(1, :, k)); scen(k).prob 1 / num_scenarios; end上面这段代码的关键点在Cholesky分解和等概率变换。我在调试时最容易出错的地方是维度搞混这里统一逻辑每个场景是一个24×4的矩阵即时段×变量但Matlab里用三维数组存储时(变量, 时段, 场景)这种排布更利于矩阵运算。后面读取某个场景时要小心别把维度弄反。3.3 同步回代缩减的Matlab实现%% 同步回代缩减 % 输入scen结构体数组目标场景数 K % 输出缩减后的场景结构体数组每个场景带概率 function scen_reduced sbr_reduction(scen, K) N length(scen); probs ones(N, 1) / N; % 记录场景编号对应的数据矩阵便于距离计算 data zeros(N, 24 * 4); % 每个场景拉平成一行 for i 1:N data(i, :) [scen(i).wind1, scen(i).wind2, scen(i).pv, scen(i).load]; end active true(N, 1); % 标记场景是否存活 num_active N; while num_active K % 计算当前存活场景两两之间的距离 idx_active find(active); m length(idx_active); % 距离矩阵 dist_mat zeros(m, m); for i 1:m for j i1:m d norm(data(idx_active(i), :) - data(idx_active(j), :)); dist_mat(i, j) d; dist_mat(j, i) d; end end % 对每个场景i找到距离它最近的场景j以及对应的距离和概率乘权值 min_dist inf(m, 1); min_dist_idx zeros(m, 1); for i 1:m row dist_mat(i, :); row(i) inf; [val, j] min(row); min_dist(i) val; min_dist_idx(i) j; end % 加权距离概率 * 距离 weighted probs(idx_active) .* min_dist; % 删除加权距离最小的那个场景 [~, del_idx] min(weighted); del_global idx_active(del_idx); % 找到距离该场景最近且存活的邻居 neighbor_idx idx_active(min_dist_idx(del_idx)); % 概率转移被删场景的概率累加到邻居 probs(neighbor_idx) probs(neighbor_idx) probs(del_global); % 标记删除 active(del_global) false; num_active num_active - 1; end % 输出 idx_keep find(active); scen_reduced scen(idx_keep); for i 1:length(idx_keep) scen_reduced(i).prob probs(idx_keep(i)); end end这段代码是同步回代最直白的实现版本每一步都对应前面讲的原理。实际运行时如果初始场景数特别大比如5000个这个双重循环会很吃力。我的经验是初始场景在1000以内时问题不大超过2000建议先把数据降维一下再算距离例如用主成分分析把24×4维降到10~20维能显著提速。缩减完成之后正规做法是做一个评估看看缩减前后的均值曲线、标准差曲线、累积分布曲线是否吻合。我习惯画两张图%% 可视化对比 figure; % 原始场景的均值曲线 mean_orig mean(data_orig, 1); plot(1:96, mean_orig, LineWidth, 1.5); hold on; % 缩减后场景的加权均值 mean_red sum(data_red .* prob_red, 1); plot(1:96, mean_red, LineWidth, 1.5); legend(缩减前均值, 缩减后均值);这里的96列是24时段×4个变量画图时用垂直虚线分隔每个变量会更直观。缩减结果若均值偏差在3%以内工程上完全可以接受。4. 调试经验与参数调优建议4.1 常见问题速查表我在把这套代码用在不同的项目场景时踩过不少坑整理成表格分享出来按出现频率排序问题现象根本原因解决方案Cholesky分解报错Matrix must be positive definite相关系数矩阵非正定变量间相关性设置不合理对相关矩阵做特征值修正或检查是否出现完全线性相关的变量组合缩减后典型场景全部挤在一起多样性差初始场景数太少或初始场景本身分布过于集中增大初始场景数建议至少500检查分布参数是否合理场景均值和原始数据均值偏差过大等概率变换时Beta逆函数对边界值0/1处理不当使用betainv时注意返回值可能为0或1加入微小裁剪区间比如限制在[0.001, 0.999]同步回代运行速度极慢距离矩阵反复计算复杂度高减少初始场景数对数据降维或改用K-means方法负荷出现负值正态分布采样时尾部产生负值电力系统负荷不可能为负对负荷采样做截断处理生成后把负值置为0或取绝对值更推荐用对数正态分布光伏夜间时段出力不为零Beta分布采样没有区分昼夜夜间时段仍然采到非零出力事先将夜间时段的出力直接置零只在白天时段做随机采样缩减概率之和不为1场景概率累加过程出现数值误差缩减结束后统一归一化probs probs / sum(probs)4.2 参数设置的实践经验与坑分布参数怎么估计是整个流程里最容易被低估的环节。很多人直接从论文里抄一套参数就用但不同地区、不同季节的风光出力特性差别很大。我建议如果手头有历史数据哪怕只有一年的小时级数据也值得按分时段拟合来做——每个时段单独拟合Beta分布或Weibull分布的参数这样能保留日内波动特性。比如光伏在上午10点和下午2点的Beta分布参数差异很大直接用全天统一参数会抹平这种差别。相关系数矩阵怎么定。如果数据充分直接拿历史出力数据计算Pearson相关系数即可。如果数据不足凭经验设置也行但要注意相关系数的绝对值不宜超过0.8否则Cholesky分解后生成的场景会表现出过强的跟随性看起来不自然。比如风电出力相关系数设到0.9两个风电场场景几乎同步波动调度模型会严重低估系统灵活性需求。缩减到多少场景合适这也是高频问题。我的经验是初步研究用5~10个场景够用论文级的结果建议分别做10、20、30个场景的对比实验看均值和分位数的变化趋势选取边际收益开始递减的那个数目。缩减场景数太少典型场景丢失极端情况优化结果偏乐观太多则场景规模带来的计算负担反而得不偿失。关于K-means与同步回代的选择我的建议是追求结果稳定、准备写论文用同步回代法做快速预实验、只是看看大致场景形态用K-means。K-means在Matlab里可以直接调kmeans函数代码量少一个数量级缺点是要多跑几次选最优聚类结果因为初始中心随机。另外K-means聚类前最好对数据做标准化否则负荷量级远大于风电光伏时聚类结果基本只由负荷变量主导风电光伏的差异被忽略。%% K-Means 缩减的极简实现 % data_orig: N行每行一个场景拉平后的向量 % K: 聚类数 [idx, C] kmeans(data_orig, K, Replicates, 10); % 重复10次取最优 prob_red histcounts(idx, 1:K1) / length(idx); % C就是典型场景矩阵每行对应一个典型场景4.3 相关性不匹配时的排查思路还有一个值得单独说的场景生成的数据单变量分布是对的但变量之间相关性显著偏离了你设定的相关系数矩阵。这个问题我在早期调试时经常遇到后来排查发现是采样量不足导致的统计噪声——1000个场景下发样本相关矩阵就能大致收敛如果只生成100个估计值会波动得比较厉害。另一个隐蔽原因是等概率变换环节的离散化截断比如Beta分布采样后把输出裁剪到[0.001, 0.999]会轻微改变线性相关性但通常影响不大。如果对相关性精度要求极高可以采用更严格的方法采样完成后计算实际相关矩阵然后做迭代修正调整初始Cholesky分解用的相关矩阵使输出相关矩阵逼近目标值。不过说句实在话对绝大多数电力系统规划调度问题初始相关的轻微偏差±0.05以内对后端优化结果几乎没影响没必要上迭代修正这种复杂方案。5. 代码的扩展方向与真实项目中的注意事项5.1 从单时段独立采样到时序相关性前面展示的代码是每个时段独立采样没有考虑时间维度上的自相关性——也就是上一小时风电出力高下一小时大概率也不会太低。这种逐时段独立采样的做法会在场景里产生锯齿状的出力曲线看起来很不真实也直接影响调度优化里爬坡约束的合理性。改进思路有两条路线。一条是用马尔可夫链模型把风速或出力的状态离散成若干区间通过历史数据统计状态转移概率矩阵再按转移概率逐时段生成场景。另一条是用多元正态分布直接生成整个时间序列时间相关性体现在协方差矩阵里——但协方差矩阵规模是96×9624小时×4变量需要大量历史数据才能可靠估计。我在实际项目中常用的是第一种先按马尔可夫链生成风速序列再经过功率曲线转换得到风电出力序列。这样做物理意义清晰而且能自然保证时序连续性。5.2 场景缩减之后怎么用缩减完成得到典型场景和概率后下一步通常是嵌入优化模型。不管你是用Yalmip还是直接用Matlab的linprog/intlinprog典型场景的用法都是一样的把每个场景的出力曲线作为确定性参数代入约束条件目标函数用场景概率加权求和。举个例子随机经济调度问题的最简形式是目标min sum_k prob_k * (火电成本 弃风惩罚 失负荷惩罚) 约束每个场景k下都要满足功率平衡、机组出力上下限、爬坡约束等这里特别注意功率平衡约束必须在每个场景下都成立而不是只在期望场景下成立。很多初学者在这里犯错——只用了期望值做调度得到的最优解在大多数场景下根本不可行。场景缩减的价值就在于此它对每个典型场景都单独建模约束条件更贴近真实运行情况。5.3 真实项目中一定要避开的三个坑第一个坑是数据口径不统一。风电、光伏、负荷的历史数据可能来自不同数据库时间分辨率不一致有的15分钟、有的1小时量纲不一致有的MW、有的kW采集时间有偏差。做场景生成前必须统一插值、统一量纲、对齐时间戳否则后面的相关性分析全是错的。这一步工作琐碎耗时但做了十年数据处理的工程师都会告诉你数据清洗占整个项目60%的时间但决定80%的结果质量。第二个坑是场景缩减目标数设得太少。有人为了省计算时间一上来就缩到3个场景结果极端场景大风无光、无风高温被合并掉调度方案在极端情况下面临切负荷风险。我的建议是至少要保证缩减后场景覆盖原始数据的最小值情景和最大值情景——也就是说缩减前先把样本中风电出力最小的一天、光伏出力最大的那天找出来单独保留再跑缩减算法防止极端情况被平均掉。第三个坑是忘记做缩减效果的定量评估。很多人在论文里只放一张缩减前后场景对比图看起来很直观但说服力不够。规范的评估指标至少有三个缩减前后各时段均值的最大偏差、标准差的最大偏差、累积分布函数的Kolmogorov-Smirnov统计量。把这些数值列出来审稿人和答辩老师都会觉得你的工作做得很扎实。%% 缩减效果评估KS检验 % 对每个时段、每个变量做两样本KS检验 % 原假设缩减前后场景集合来自同一分布 h_ks zeros(24, 4); p_ks zeros(24, 4); for t 1:24 for v 1:4 [h_ks(t, v), p_ks(t, v)] kstest2(data_orig(:, (v-1)*24t), ... expand_scen_reduced(:, (v-1)*24t)); end end % h1表示拒绝原假设说明该时段缩减前后分布差异显著6. 我个人的实操体会这套多场景随机生成与缩减代码我在不同项目里迭代过四五个版本。最初版只是简单地对每个变量独立采样、用K-means缩减那时候做完的结果被合作方一句场景怎么都在平均曲线附近就给打回来了。后来重构成现在这套相关性处理等概率变换同步回代缩减的架构效果才稳定下来。中间最大的体会是场景生成不是目的分布保真才是目的——不要追求生成过程的数学华丽要反复确认缩减后的场景集合在统计特性上经得起对比。另外一个小技巧调试的时候不要把初始场景数直接设1000先用100个把整个流程跑通检查生成的曲线形态是否符合物理直觉——风电曲线有无跳变、光伏是否只在白天有出力、负荷是否出现负值。确认代码逻辑无误后再把场景数调大做正式实验。这样能省下大量等待循环的时间。最后再分享一个我在展示结果时的做法把10个典型场景画在同一张图上不标颜色顺序而是按概率从大到小排列用线条粗细区分概率大小。这样一眼就能看出系统最可能落在什么运行状态比摆一摞曲线图有说服力得多。这套代码本身改造成本很低不管是加储能、加多节点负荷还是改成季度场景分析你都能直接在现有框架上扩展。