处理含风光出力的随机优化问题时最让人头疼的往往不是优化算法本身而是怎么把“不确定性”这件看不见摸不着的事翻译成模型能读懂的输入。我之前在做一个微电网日前调度方案时图省事直接拿预测曲线喂给优化器结果调度计划拿到现场一跑就露馅了光伏午间的尖峰、夜间风电的随机波动根本没被模型覆盖到。后来老老实实改用场景法先把风光不确定性生成成大量可能出现的场景再用概率距离削减法压缩到可计算的规模问题才算真正解决。这篇文章就围绕“MATLAB蒙特卡洛法概率距离削减法”这套组合拳展开把风光场景生成与削减从原理到代码一次讲透。适合刚接触随机优化、鲁棒优化的研究生或者是做新能源并网、微电网/主动配电网调度、电力市场出清的工程师参考。内容不绕弯子直接给思路、给代码、给参数、给避坑经验。1. 场景法解决的核心问题为什么单点预测不够用1.1 从一次调度方案翻车说起做含风电光伏的优化调度最基本的做法是用预测曲线作为输入跑一个确定性的优化模型。但这里有个致命前提——预测是准的。实际中风电和光伏的出力波动远比想象中剧烈天气预报给的一条曲线只是众多可能中的一条你要是只拿中线去优化遇到实际偏差大的日子要么弃风弃光严重要么电量不平衡被考核。我在项目里吃过一次亏模型算出来的“最优”购电计划在实际运行中连续三天出现偏差原因是那几天实际辐照比预测值高了将近30%光伏实际出力把计划打乱了。回头做分析时才意识到问题的本质不是预测不准而是我压根没在模型里给“不确定性”留位置。从那以后凡涉及风光出力的优化问题我都会先做场景生成与削减这步预处理把随机性显式地表达成一组带概率的离散场景。1.2 场景法的基本逻辑与完整流程场景法的核心思想其实很朴素用一堆“有可能发生”的出力曲线去逼近真实的风光出力随机过程。每一条曲线就是一个场景每个场景附带一个发生概率。这样做的好处是场景可以直接写进优化模型的目标函数和约束里比如目标函数变成“所有场景下成本的期望最小化”约束变成“在绝大多数场景下都满足”这样优化出来的方案对不确定性天然具有适应能力。一个标准的场景分析流程是先基于风速、光照的历史统计特性用蒙特卡洛法随机抽样生成大量初始场景比如1000条风电曲线、1000条光伏曲线然后把其中“长得差不多”的场景合并删减用概率距离削减法压缩成少量具有代表性的场景比如10条、20条最后把这组削减后的场景连同各自概率一起喂给上层优化模型。这里生成和削减是一对配合动作生成阶段追求的是“覆盖充分”宁可多不能漏把所有可能的天气形态都抽样出来削减阶段追求的是“保留精髓”在尽量不损失概率分布信息的前提下把场景数压到优化模型能承受的规模。1.3 为什么选蒙特卡洛概率距离削减这对组合做场景生成有很多路线比如基于历史数据直接聚类、用马尔可夫链生成时序、用生成对抗网络学分布。蒙特卡洛法能成为最主流的方案核心原因是它简单可靠只要你能给出风速和光照的概率分布模型就能在几分钟内生成几千条出力序列数学基础扎实工程上也好解释。场景削减的路线也不少最简单的就是K-means聚类。但聚类方法有一个问题它默认所有场景等概率聚类中心在概率意义上是“平均”的容易把低概率但影响大的极端天气场景磨平。概率距离削减法则不一样它显式地把场景的概率和场景间的距离结合在一起削减过程中会考虑哪些场景该合并、概率该怎么转移在保持概率分布形态方面明显更有优势。这也是为什么它在电力系统的随机规划里几乎是标准配置。2. 蒙特卡洛生成初始风光场景从概率分布到出力序列2.1 风速的Weibull分布与风电出力换算风速概率建模行业内用两参数Weibull分布最多。它的概率密度函数长这样[ f(v) \frac{k}{\lambda}\left(\frac{v}{\lambda}\right)^{k-1} e^{-(v/\lambda)^k} ]其中λ是尺度参数和平均风速直接相关k是形状参数控制风速分布的“胖瘦”。国内比如三北地区k大概在1.8到2.2之间λ在6到9之间南方风电场因为风况更不稳定k有时候会到1.5以下。你可以直接用历史风速数据去拟合这两个参数MATLAB里有现成工具也可以按经验取值测试。采样用MATLAB的wblrnd(λ, k, 行, 列)一行代码就能生成一批风速样本。需要注意的是wblrnd的输入参数顺序是先λ后k这和很多教材里写先k后λ的习惯不一样我第一次写就搞反了生成的分布完全不对调试了半天才发现是参数顺序的问题。拿到风速样本后不是直接当作风电出力还要经过风机的“输出特性曲线”转换。典型的风机有三个关键风速切入风速通常3m/s、额定风速通常12m/s、切出风速通常25m/s再加上一个风电额定功率。出力转换规则是% windParam [lambda, k, v_ci, v_r, v_co] % 输出风电出力归一化到 [0,1] p_wind zeros(T, N); p_wind(v v_ci v v_r) (v(v v_ci v v_r) - v_ci) / (v_r - v_ci); p_wind(v v_r v v_co) 1; % v v_ci 或 v v_co 的场景出力自动保持为 0这实际上是一个三段式函数风速太低不发电、中等风速时出力随风速线性上升、额定风速以上出力恒定、切出风速后出于安全保护停机。这个非线性转换是蒙特卡洛场景生成里最关键的一步因为风速分布经过这个“歪曲”之后出力分布形态会完全变掉——中间风速区间的出力概率被拉平这个特性会直接影响后面的场景削减结果。2.2 光照的Beta分布与光伏出力换算光伏出力的源头是光照辐照度工程上常用Beta分布来描述辐照度的随机性。之所以不用正态分布是因为光照强度天然是非负的而且有物理上限日最大辐照度Beta分布定义在[0,1]区间只要乘一个最大辐照度系数就能很好地贴合实际。Beta分布有两个形状参数α和β它们的取值决定了分布形态α比较大时光照容易偏高β比较大时光照容易偏低。实际拟合时晴天的α、β和阴雨天差别非常大所以一般按天气类型分开建模或者用混合Beta分布。入门阶段用单组参数也够用。光伏出力换算比风电简单通常认为光伏出力与辐照度近似成正比% solarParam [alpha, beta, r_max] r betarnd(alpha, beta, [T, N]) * r_max; % 辐照度 p_solar r / r_max; % 光伏出力归一化到 [0,1]如果想把温度和转换效率也考虑进去可以在上面公式基础上加一个温度修正系数但那样要引入环境温度的概率模型复杂度高了不少。对大多数场景分析来说辐照度正比模型已经能抓住光伏波动的主要特征可以先把基础版本跑通再加细节。2.3 MATLAB场景生成函数实现把上面两部分组装起来就是一个完整的场景生成函数。我实际用的版本长这样function scen gen_wind_solar_scenarios(N, T, windParam, solarParam, seed) % 蒙特卡洛生成风光联合出力场景 % 输入: % N 初始场景数 % T 调度时段数 % windParam [lambda, k, v_ci, v_r, v_co] % solarParam [alpha, beta, r_max] % seed 随机种子便于结果复现 % 输出: % scen 维度为 2*T x N 的矩阵 % 前 T 行是风电出力后 T 行是光伏出力 % 每列代表一个可能的出力场景 if nargin 4 ~isempty(seed) rng(seed); end % --- 风速采样与风电出力 --- lambda windParam(1); k windParam(2); v_ci windParam(3); v_r windParam(4); v_co windParam(5); v wblrnd(lambda, k, [T, N]); p_wind zeros(T, N); p_wind(v v_ci v v_r) (v(v v_ci v v_r) - v_ci) / (v_r - v_ci); p_wind(v v_r v v_co) 1; % --- 辐照度采样与光伏出力 --- alpha solarParam(1); beta solarParam(2); r_max solarParam(3); r betarnd(alpha, beta, [T, N]) * r_max; p_solar r / r_max; % --- 合并 --- scen [p_wind; p_solar]; end生成场景时有个很容易忽略的点场景矩阵的维度设计。建议把矩阵做成“维度×场景数”的形式也就是每一列是一个完整场景。这样后续做削减、算距离都是对列操作逻辑非常顺。如果做成“场景数×维度”pdist2那块还得转置来转置去代码可读性差还容易出维度错误。2.4 时序相关性和采样改进LHS、AR(1)上面的基础版本假设每个时段的采样互相独立这隐含了一个不太真实的假设某一时刻风速大下一时刻风速小完全是随机跳跃的。实际的风速过程有明显的时序相关性白天到晚上的变化是渐进的不是随机打点。处理时序相关最简单的办法是用一阶自回归模型AR(1)对风速做时间上的平滑[ v_{t1} \mu \phi (v_t - \mu) \varepsilon_t ]其中φ是自回归系数通常在0.7~0.95之间φ越大说明风速越“懒”前后时刻关联越强ε_t是白噪声。这样生成的风速曲线就会带上自然的惯性特征风电出力曲线也不再是高频跳变的锯齿状。代价是代码量增加不少而且要对风速先做正态化处理入门阶段可以先不搞但心里要有这根弦。另一个性价比很高的改进是拉丁超立方采样LHS。普通蒙特卡洛是完全随机撒点容易在某个区间堆一堆样本、另一个区间空荡荡LHS先把概率空间均匀分成N个区间再从每个区间各抽一个点这样样本对分布的覆盖均匀得多。对随机优化来说LHS通常能以更少的场景数达到同样的概率覆盖效果。MATLAB实现也不复杂先对[0,1]均匀分层抽样再调Weibull和Beta分布的分位数函数逆变换回去。不过要提醒一句LHS适合改善样本覆盖但不会自动解决时序相关性问题。想全面一点的话LHS和AR(1)可以一起上生成质量会明显上一个台阶。3. 概率距离削减把1000个场景压缩到10个还不出戏3.1 概率距离Kantorovich距离到底在衡量什么假设你生成了1000个场景那优化模型里每个随机变量就得有1000个可能值约束规模直接膨胀1000倍绝大多数情况下模型根本解不动。场景削减就是为了把这个数量压下来同时让削减后的场景集合在概率意义上仍然“接近”原始场景集合。“接近”这个词怎么量化概率距离削减法用的是Kantorovich距离也叫最优传输距离。直观理解是要把原始场景集合的“概率质量”搬到削减后场景集合上需要付出的最小“搬运成本”。场景A和场景B之间的距离度量了搬运单个概率单位的成本Kantorovich距离就是所有场景搬运成本的总和。它同时考虑了场景间的几何距离和各自的概率这和K-means这类纯几何聚类有本质区别。在实际削减操作里核心思想是每一步找到一个“最不值得保留”的场景把它删掉顺便把这个场景的概率合并到离它最近的保留场景上。这样场景总数逐步减少而概率总和始终是1损失的Kantorovich距离被控制在最小。两个最常用的实现分别是快速前向选择法和同步回代削减法下面分别讲。3.2 快速前向选择算法的原理与MATLAB实现快速前向选择的思路是“从N个场景里挑K个作为保留集合”每一步从当前场景集合里找出“到其他场景距离最近的那个”既然它和别的场景长得太像留着的价值最低就把它删掉。重复这个操作直到只剩K个场景。删掉的场景概率怎么处理合并到离它最近的那个保留场景上。这个过程保持了场景集合的总概率始终为1。function [keep_idx, keep_prob] fast_forward_selection(scen, prob, K) % 快速前向选择法场景削减 % scen: nDim x N 矩阵N个初始场景 % prob: N x 1 向量场景初始概率通常为 1/N % K: 目标保留场景数 % 输出: % keep_idx: K x 1 向量保留场景在原矩阵中的列索引 % keep_prob: K x 1 向量保留场景的更新后概率 N size(scen, 2); idx (1:N); % 当前场景在原矩阵中的索引 P prob(:); % 当前概率 D pdist2(scen, scen, euclidean); while length(idx) K D_self D; D_self(1:size(D_self,1)1:end) Inf; % 对角线置为Inf % minDist_i 场景i到最近其他场景的距离 % nearestIdx_i 最近场景在当前D矩阵中的行号 [minDist, nearestIdx] min(D_self, [], 2); % 快速前向选择删除“最近距离最小”的场景 [~, delPos] min(minDist); keepPos nearestIdx(delPos); % 概率合并到最近场景 P(keepPos) P(keepPos) P(delPos); % 删除对应行、列 D(delPos, :) []; D(:, delPos) []; idx(delPos) []; P(delPos) []; end keep_idx idx; keep_prob P / sum(P); % 归一化保证概率和为1 end这段代码是教学优化版本每一步重新计算最近距离逻辑清楚方便理解。如果初始场景N到几千循环几百次耗时也在可接受范围内。要是追求性能可以维护一个最近距离矩阵每次只更新受影响的部分但代码复杂度会明显上升建议先把基础版本跑通再考虑优化。3.3 同步回代削减算法的原理与MATLAB实现同步回代削减的删除标准和前向选择正好相反。它不是找“最冗余”的场景删而是算每个场景的“删除代价”——这个场景的概率乘以它到最近其他场景的距离。删除代价最小的场景意味着去掉它对整体概率分布影响最小所以优先删除它。这个实现特别适合处理概率不均匀的场景集。如果某个场景的概率很低即使它离其他场景很远删除它的代价也可能很小反过来一个高概率场景哪怕和别的场景很近删除代价也会很大于是更倾向于保留。function [keep_idx, keep_prob] backward_reduction(scen, prob, K) % 同步回代削减法场景削减 % 输入输出与 fast_forward_selection 一致 N size(scen, 2); idx (1:N); P prob(:); D pdist2(scen, scen, euclidean); while length(idx) K D_self D; D_self(1:size(D_self,1)1:end) Inf; [minDist, nearestIdx] min(D_self, [], 2); % 同步回代删除“概率×最近距离”最小的场景 J P .* minDist; [~, delPos] min(J); keepPos nearestIdx(delPos); P(keepPos) P(keepPos) P(delPos); D(delPos, :) []; D(:, delPos) []; idx(delPos) []; P(delPos) []; end keep_idx idx; keep_prob P / sum(P); end这两种算法的迭代框架几乎一样差别就在删除准则那一行——前向选择看minDist同步回代看P .* minDist。这个差别看起来小实际削减效果差异不小下面展开说。3.4 两种削减方法怎么选从我实际测试的体验来说不存在绝对优劣要看你的下游任务是什么。前向选择更倾向于保留“空间上分散”的场景。因为每次删的是离所有场景最近的那个那些孤零零待在角落里的极端场景反而不容易被删。这在高比例新能源系统里特别重要——极端低风、持续阴天这种场景虽然概率低但对电力电量平衡冲击最大如果被削减掉系统就会低估风险。我做过一个案例用前向选择保留的场景集合最差场景下的失负荷量明显比同步回代高但这恰恰更接近原始1000个场景里真实存在的风险水平。同步回代更偏向保留“概率高”的场景。它按概率加权删除代价低概率场景更容易被优先合并掉。这在概率分布比较集中、极值不重要的场景下性能很好削减后期望值误差更小。如果是做优化调度追求期望成本最小同步回代通常更好用如果是做风险评估、可靠性分析我会优先用前向选择。当然也可以两种都跑一遍对比结果后再定多花的时间成本很低。4. 削减效果怎么看评价指标与调参经验4.1 三个必看的数字概率和、期望出力、距离指标削减完别急着往优化模型里塞先过三个体检指标。第一个是概率和。削减后的所有场景概率加起来必须是1这是硬指标。我见过不少人在循环里合并概率时忘了删概率向量里的对应项或者最后少归一化一步导致概率和变成0.98、1.05这种数优化模型的目标函数直接被带偏。代码里我最后都加了归一化但自己写的时候还是要养成检查习惯。第二个是期望出力误差。对比削减前所有场景的平均出力和削减后按概率加权的平均出力这个误差一般控制在个位数百分比以内就算不错。具体算哪段出力要看你关心的是总发电量还是某个时段的出力一般先看总出力期望% 原始场景期望总出力所有场景等概率 orig_exp mean(sum(scen, 1)); % 削减后期望总出力按概率加权 reduced_exp sum(sum(reduced_scen) .* keep_prob); fprintf(原始期望总出力: %.4f\n, orig_exp); fprintf(削减后期望总出力: %.4f\n, reduced_exp); fprintf(相对误差: %.2f%%\n, abs(orig_exp - reduced_exp) / orig_exp * 100);我测试过几组参数N1000削减到K20时期望出力误差通常能压在1%以内效果相当好。第三个是削减前后的Kantorovich距离。这个数值越小说明削减对概率分布的扭曲越小。MATLAB里没现成函数但可以自己写对原始每个场景找它到最近保留场景的距离乘以它的概率后求和。实现也不复杂D_reduced pdist2(scen, reduced_scen, euclidean); [minD, ~] min(D_reduced, [], 2); kd sum(prob .* minD); % 近似Kantorovich距离注意这个算的是原始场景集到削减场景集的单向距离严格意义上不是完整对称的Kantorovich距离但在工程上做相对比较完全够用。4.2 K选多大合适从优化模型复杂度反推削减后的场景数K没有绝对标准更多看下游优化模型能承受多大规模。场景数直接影响决策变量的数量一个含风光不确定性的机组组合模型如果每个时段都有随机约束场景数K等于把约束或决策变量放大K倍。K20时几分钟能解完的模型K100时可能要跑几个小时甚至解不动。我的经验是分三步处理先跑K5、10、20、50四个档位分别记录优化目标值然后画一条“场景数-目标值”曲线找到目标值变化趋于平缓的拐点最后取拐点对应的K作为正式配置。比如我做过的一个微电网日前经济调度案例K从5增到20时最优成本明显下降但K从20增到50时变化不到0.5%就果断定在K20。初始生成场景数N一般取500到2000就够用了。N太小削减出的场景缺乏多样性N超过2000pdist2计算N×N的距离矩阵开始吃内存削减循环的耗时也明显上升收益却很小。综合考虑准确性和计算负担N1000、K20是我最常用的配置。4.3 画图检查削减后的场景带是否还认识指标归指标我每次削减完还会画一张图用眼睛检查一下削减效果。具体做法是把原始场景的“场景带”和削减后场景叠在一张图里看figure(Color, w); % 原始场景前50条画浅色细线 plot(scen(1:T, 1:50), Color, [0.8 0.8 1], LineWidth, 0.3); hold on; % 削减后的场景画红色粗线 plot(reduced_scen(1:T, :), r, LineWidth, 1.5); xlabel(时段); ylabel(风电出力 (pu)); title(削减前后风电场景对比);观察几个点削减后的场景是否还覆盖原始场景带的主要形态是否存在明显偏离所有原始场景的“幽灵场景”场景之间是不是过于集中丢失了大范围波动的信息。正常削减结果应该是红线和蓝色区域的主要包络形态一致同时K条红线之间保持一定差异。如果削减后的场景挤成一团说明K取得太大或者削减准则和你的分布不匹配需要调整。5. 常见报错与排查技巧实录5.1 快速排查表做场景生成与削减代码本身不长但报错的点不少。我把这半年被问得最多的几个问题整理成一张表基本覆盖了新手会遇到的大部分坑。报错/现象可能原因解决办法wblrnd生成的场景出力全为0参数λ、k顺序写反或风机参数设置不合理确认wblrnd(lambda, k)检查v_ci、v_r、v_co之间大小关系pdist2报错或内存不足场景数N过大或场景矩阵维度放错了方向N控制在2000以内确认scen是nDim x N方向削减后概率和不是1删除场景时概率没正确转移或没归一化循环结束后统一除以sum(P)削减结果每次都不一样蒙特卡洛采样没固定随机种子生成场景前用rng(固定值)固定种子光伏出力场景恒为0Beta分布参数α、β设得不合理采样值接近0尝试α5、β3这类偏高的参数检查r_max是否太小削减后场景太集中失去极端场景削减准则过于偏向高概率场景改用快速前向选择法或手动保留极端场景再削减没有Statistics工具箱跑不了pdist2工具箱依赖用下面给出的替代代码关于pdist2依赖这里给一个不依赖工具箱的欧氏距离实现可以直接替换函数名function D euclidean_dist_custom(X, Y) % X: nDim x N1, Y: nDim x N2 % 返回 N1 x N2 距离矩阵 D sqrt(max(sum(X.^2,1) sum(Y.^2,1) - 2 * X * Y, 0)); end这个版本用矩阵乘法展开欧氏距离速度不差而且不依赖任何额外工具箱。如果你的MATLAB没有安装Statistics and Machine Learning Toolbox把削减函数里的pdist2换成这个函数就行。5.2 两个容易被坑的细节第一个坑是场景矩阵的“展平方式”不统一。如果你把每个场景存成nPeriod x nVar的矩阵削减函数却按向量来处理距离计算就会错乱。我的建议很简单从生成函数开始就统一用“维度×场景数”的形式风电和光伏按T行堆叠后面所有环节都不改省去一堆维度转换的麻烦。第二个坑是削减后概率的“信息流失”。削减完场景概率会发生转移原来每个场景都是1/N削减后保留场景的概率差别很大。这时候如果你把削减后的场景当作等概率直接用等于白做了削减。我把概率存在keep_prob里传到优化模型时目标函数会写成sum(keep_prob .* 每个场景的成本)而不是简单平均。很多初学者在这步翻车优化结果和原始场景法的结果对不上查了半天才发现概率没带进去。5.3 一个实用的调试小技巧场景生成和削减是随机算法出了bug难定位。我的习惯是先固定随机种子再用小规模跑通逻辑N50、K5削减完把保留场景索引、概率全部打印出来手算核对一遍。比如两个几乎一样的场景削减后应该合并成一个概率翻倍完全相同的两个场景合并后概率是原来两个之和。确认这些基本规律都对再放大到N1000正式跑。这样能过滤掉大部分逻辑错误而不是在1000个场景的矩阵里大海捞针。结尾聊聊我的一点实际体会场景生成与削减这套流程在我做过的微电网优化、配电网可靠性评估、售电公司购电策略好几个项目里都派上了大用场。刚开始总觉得“多生成点场景、多保留点场景”更保险后来实测下来发现关键不在于场景多而在于削减方法合不合适、概率有没有用对。快速前向选择和同步回代削减我都会跑一遍对比一下关键指标再定最终版本花不了多少时间但结果会稳不少。最后再分享一个小技巧削减后的场景集不要一次性丢到复杂优化模型里。先挑一两个典型场景单独跑通确定性模型确认边界条件、约束格式都对再切换到全场景的随机优化版本。这样排查问题快也不会在模型报错时分不清是场景的问题还是优化的问题。这套方法和代码我已经在多个项目里反复用直接拿去改改参数就能跑通。