微电网调度中的风光场景生成与削减：蒙特卡洛+概率距离法实战

做微电网调度优化的时候,我大部分时间其实不是在写约束,而是在想办法对付不确定性。前段时间接了一个模拟项目,给定一片风电和一片光伏作为可调度的分布式电源,目标是最小化运行成本。第一步很常规:用蒙特卡洛法把未来24小时的风光出力抽成5000个随机场景。第二步就出问题了——这些场景直接塞进混合整数规划模型后,求解器跑了一个多小时还没收敛,内存占用一度顶到十几GB,最后只能强行掐掉。

我当时的第一个想法是减少抽样数量,但样本太少又会丢掉尾部分布特征,优化结果会过于乐观。后来换成"先生成、后削减"的思路,用概率距离削减法把5000个场景压到20个左右,求解时间降到了几分钟,关键统计特征也保住了。这篇把整套流程完整记录一下,包括蒙特卡洛生成风光场景的建模方式、概率距离削减的数学原理、可用性还不错的MATLAB实现,以及我在实际项目中踩过的几个坑。给同样在折腾风光场景生成与场景削减的朋友做个参考。

1. 被1000个随机场景拖垮的优化模型:为什么需要场景生成与场景削减

1.1 风光出力随机性到底体现在哪里

风电和光伏出力的随机性,本质上来自气象过程。风速受地形、气压、温度影响,呈现出明显的偏态分布;光照则受云层遮挡、大气透明度影响,在一天之内呈现强烈的时段特征。

所以在做随机优化之前,业内一般会对历史数据做概率建模。风速常用双参数Weibull分布,光照强度常用Beta分布。这两个分布不是随便选的——Weibull能很好描述风速"大部分时间不高不低,偶尔有强风"的偏态特征;Beta分布的取值区间是[0,1],和归一化后的光照强度天然匹配。

但要注意,分布模型只描述了"某一时刻的统计特性"。实际问题里还有时序相关性:下午3点的风速和下午2点的风速不会是毫不相关的两个独立样本。如果只按分布独立抽样,生成出来的场景是一个个"切片",而不是一条条"曲线",这个后面会专门讲。

1.2 场景数量怎么影响优化模型

随机优化里最常见的处理方式是场景法:把不确定参数离散成有限个场景,每个场景对应一组决策变量和约束条件。以机组组合为例,如果未来24小时的风光出力有1000个场景,每个场景都要单独刻画机组启停和出力约束,那么模型里的变量数是原来的1000倍。

我那个模拟项目里,优化模型本身有数百个连续变量和数十个二值变量,配合5000个场景后,约束矩阵的行数直接冲到千万级。求解器的运行时间从分钟级变成小时级,内存占用也在不断攀升。这还不是最坏情况——如果做两阶段随机规划,第二阶段变量的数量还和场景数相乘,场景数稍微多一点,模型就不是"跑得慢"的问题,而是根本构建不出来。

所以说,场景削减不是锦上添花,而是随机优化工程落地的一个硬门槛。它的目标非常明确:用尽量少的代表性场景,保持原始概率分布的关键统计特征,让优化模型的规模控制在求解器能接受的范围。

1.3 场景削减的任务定义

从数学上看,原始场景集是一个离散概率分布——每个场景带一个概率权重,比如5000个场景权重相等(各1/5000)。场景削减要做的是找到一个支撑点更少的离散分布,让两个分布之间的"概率距离"尽可能小。

这里的输出有两部分:一是选出来的保留场景本身,二是给每个保留场景重新分配的概率。保留场景不再是等概率的——有些场景吸收了周围被删场景的概率,权重会变大,成为"代表性场景";有些场景权重会变小,对应那些比较偏、但又不能完全忽略的情况。

很多初学者容易把场景削减当成简单的"随机抽几个出来"或者"聚类取中心",这两种做法都能用,但概率距离削减法在理论保证上更严谨,而且保留了原始场景的真实物理特征。聚类取中心产生的场景可能是原始场景中没有的组合状态,分布信息会失真,后面详说。

2. 蒙特卡洛生成风光场景:从概率分布到24时段出力矩阵

2.1 风速与风功率:Weibull分布和功率曲线

风速是场景生成的第一层随机变量。我常用的做法是逐时段建立风速概率模型。对某个站点,先统计一天24个时段的平均风速,再用Weibull分布描述同一时段内的风速波动。

Weibull分布的概率密度函数为:

f(v) = (k/c) * (v/c)^(k-1) * exp(-(v/c)^k)

其中k是形状参数,决定了风速分布的"陡峭程度";c是尺度参数,和平均风速大致成正比。k通常取1.8到2.5之间,c则根据站点历史风速拟合。

有了风速样本之后,还要通过风机功率曲线折算成出力。风机的切入风速、额定风速、切出风速三个节点决定了功率曲线的形状,典型简化模型是:

• 风速小于切入风速或大于切出风速时,出力为0
• 风速在切入风速和额定风速之间时,出力随风速线性上升
• 风速在额定风速和切出风速之间时,出力维持在额定功率

在MATLAB里,这一步用wblrnd函数直接抽样即可,然后用向量化逻辑操作完成功率曲线折算。

2.2 光照与光伏出力:Beta分布和转换模型

光伏出力的随机性来源主要是辐照度。我把每个时段的辐照度先归一化到[0,1],再假设它服从Beta分布。Beta分布有两个参数α和β,改变这两个值可以模拟不同云况下的辐照度特征:晴天时α较大、β较小,出力曲线高而稳定;阴天时α和β的比值接近,均值低,波动大。

辐照度样本出来后,光伏出力用简化线性模型折算:

P_pv = P_rated * (I_t / I_0)

其中I_t是某时段辐照度,I_0是标准辐照度(一般取1000W/m²),P_rated是光伏额定容量。更精细一点可以加温度修正系数,但作为场景生成的初始环节,这个线性模型已经够用。

2.3 MATLAB抽样核心代码:原始场景集怎么造出来

这是生成原始场景集的核心代码。我以一个典型日为例,先构造风速和辐照度的逐时段分布参数,再抽样生成N个场景,最后折算成风/光功率矩阵的拼接形式。

% 场景生成参数 N = 5000; % 原始场景数量 T = 24; % 调度时段数 Pw_rated = 0.8; % 风电额定功率 MW Ppv_rated = 0.5; % 光伏额定功率 MW % 风机功率曲线参数 v_cin = 3; % 切入风速 m/s v_rated = 12; % 额定风速 m/s v_co = 25; % 切出风速 m/s % 逐时段Weibull尺度参数(示意值,实际应由站点风速拟合) c_w = 4.5 + 2.5 * sin(2 * pi * (1:T) / 24 - 0.3); k_w = 2.0; % Weibull形状参数 % 生成风速场景 wind = zeros(N, T); for t = 1:T wind(:, t) = wblrnd(k_w, c_w(t), N, 1); end % 风速转风功率 wind_power = zeros(N, T); for t = 1:T v = wind(:, t); p = zeros(N, 1); linear_zone = v >= v_cin & v < v_rated; rated_zone = v >= v_rated & v <= v_co; p(linear_zone) = Pw_rated * (v(linear_zone) - v_cin) / (v_rated - v_cin); p(rated_zone) = Pw_rated; wind_power(:, t) = p; end % 典型日辐照度曲线(归一化,日出6点日落18点) I_mean = max(0, sin(pi * ((1:T) - 6) / 24)); alpha_b = 2.5; beta_b = 2.5; % 生成光伏辐照度场景 irr = zeros(N, T); sun_hours = 7:18; for t = sun_hours irr(:, t) = betarnd(alpha_b, beta_b, N, 1) .* I_mean(t); end % 辐照度转光伏功率 pv_power = Ppv_rated * irr; % 组合场景矩阵:每行是一个场景,包含24时段风电+24时段光伏 scenes = [wind_power, pv_power];

代码跑完,得到一个5000×48的矩阵。每一行代表一个"日场景",前半段是风电24时段出力,后半段是光伏24时段出力。这个矩阵就是后续场景削减的输入。

注意:这里的分布参数是示意值。真实项目中必须用站点历史风速、辐照度数据做参数估计,不然生成的场景分布会和实际严重偏离。这个问题在第5节专门讲。

3. 概率距离削减法的数学直觉:Kantorovich距离与两类主流算法

3.1 为什么是概率距离而不是聚类

很多人会问:场景削减直接用k-means聚类不就行了吗?k-means确实能压缩场景数量,但它有两个问题。

第一,k-means的聚类中心是簇内样本的平均,这个"平均场景"在原始场景集中往往不存在。风功率曲线经过平均之后,高峰被抹平,低谷被抬高,最后得到一条"温和得不像话"的代表曲线。这种平均化效应会低估系统的调峰压力,进而影响储能配置和机组组合结果。

第二,k-means优化的是"簇内平方距离",没有显式考虑概率分布之间的差异。两个场景簇即便内部距离很小,如果簇间的概率权重分配和原始分布不一致,整体统计特征还是会失真。

概率距离削减法走的是另一条路:保留原始场景中的真实样本作为代表,只调整这些代表场景的概率权重。削减后的分布Q中的每个支撑点都是原始场景P中的某个样本,Q通过吸收被删场景的概率来逼近P。由于支撑点不会超出原始样本集合,削减结果保留了真实的物理相关性,不会出现"风功率持续稳定在中间值"这种现实中很少见的状态。

3.2 Kantorovich距离怎么算

概率距离削减法使用的度量是Kantorovich距离,也叫Wasserstein距离。对两个离散分布P和Q,Kantorovich距离的一种方便计算形式是:

W(P, Q) = Σ p_i * min_j c(s_i, z_j), 对所有被削减的场景 s_i

这里的c(s_i, z_j)是场景s_i和保留场景z_j之间的距离度量,一般用欧氏距离;pj是被删场景s_i对应的概率。直观理解:每个被删的场景都要找到一个最近的保留场景"投靠",它贡献的距离代价等于自身概率乘以最近距离,所有被删场景的代价之和就是削减后的概率距离。

距离越小,说明被删场景距离保留场景越近,信息损失越小。这个方法把"场景削减"转化成一个最小化Kantorovich距离的组合优化问题。实际操作中,精确最优解很难求,所以工程上普遍使用两种启发式算法:快速前向选择和同步回代消除。

3.3 快速前向选择与同步回代消除的算法步骤

快速前向选择是"从零开始往上加"。过程是:保留集初始为空,每一步从剩余场景中选一个加入保留集,使得"当前保留集对全体场景的Kantorovich距离最小"。重复这个步骤,直到保留集达到指定数量K。

这个算法的特点是:每一步都是基于全局距离增量做的贪心选择,场景一旦选进来就不会再被删掉。前向选择生成的保留集倾向于覆盖原始分布的"骨架",比较适合需要保留极端场景的场合。

同步回代消除则是"从全集开始往下删"。过程是:初始把N个场景全部放入保留集,每一轮计算"如果把某个场景删除,它对当前保留集的Kantorovich距离会增加多少",选代价增量最小的那个场景删除,并把它的概率转移给最近保留场景。重复删除,直到剩余K个场景。

这个算法的特点是"全局出发、逐步收缩",计算速度通常比前向选择快,但对极端场景的保留效果不如前向选择,更倾向于保留"主流区域"的场景。

两种算法不是互斥的。我在实际项目里的经验是:初始场景数在几千以内时,前向选择的保留效果更稳定;初始场景数上万时,回代消除的矩阵运算优势更明显,可以先用回代消除粗削到200个,再用前向选择精削到20个。

4. 手写前向选择与回代消除:MATLAB代码走读与典型输出

4.1 前向选择的核心实现

先说前向选择的MATLAB实现。核心是维护一个距离矩阵D,第(i,j)个元素是场景i和场景j的欧氏距离。因为这个矩阵是N×N的,场景数5000时占内存约190MB,还能接受;如果场景数过万,建议分块计算。

function [idxKeep, pKeep] = fastForwardSelect(data, K) % data: N x T 原始场景矩阵 % K: 目标保留场景数 N = size(data, 1); p = ones(N, 1) / N; D = pdist2(data, data); % 欧氏距离矩阵 idxKeep = zeros(1, K); for m = 1:K if m == 1 % 第一个保留点:选到所有场景加权距离最小的场景 bestCost = inf; for j = 1:N cost = sum(p .* D(:, j)); if cost < bestCost bestCost = cost; pick = j; end end else % 已有保留集到全场景的最小距离 distToKeep = min(D(:, idxKeep(1:m-1)), [], 2); bestCost = inf; for j = 1:N if any(idxKeep == j) continue; end % 把候选j并入保留集后,全场景新增的总距离代价 distCand = min(distToKeep, D(:, j)); cost = sum(p .* distCand); if cost < bestCost bestCost = cost; pick = j; end end end idxKeep(m) = pick; end % 概率分配:每个被删场景找最近保留场景,把概率累加过去 [~, nearest] = min(D(:, idxKeep), [], 2); pKeep = p(idxKeep); for i = 1:N if any(idxKeep == i) continue; end j = nearest(i); pKeep(j) = pKeep(j) + p(i); end end

这段代码的复杂度是O(K·N·M),这里的M是当前保留集大小。K一般不超过50,N是5000或1万,实际跑起来也就几秒钟到半分钟,完全可以接受。如果嫌慢,可以用向量化手段把内层j循环改成矩阵计算,但代码可读性会下降,我一般只在正式工程里做这种优化。

4.2 回代消除的核心实现

回代消除的逻辑正好反过来。初始所有场景都在保留集里,迭代删除"代价增量最小"的场景,直到剩下K个。

function [idxKeep, pKeep] = backwardReduction(data, K) N = size(data, 1); p = ones(N, 1) / N; D = pdist2(data, data); keep = true(N, 1); while sum(keep) > K keepIdx = find(keep); bestCost = inf; delIdx = 0; for j = keepIdx.' % 当前保留集中除j外,距离j最近的点 others = setdiff(keepIdx, j); [minDist, ~] = min(D(j, others)); cost = p(j) * minDist; if cost < bestCost bestCost = cost; delIdx = j; end end % 删除代价最小的场景,概率转移给最近保留场景 keep(delIdx) = false; keepIdx = find(keep); [~, pos] = min(D(delIdx, keepIdx)); nearestIdx = keepIdx(pos); p(nearestIdx) = p(nearestIdx) + p(delIdx); p(delIdx) = 0; end idxKeep = find(keep); pKeep = p(idxKeep); end

这个实现每轮循环都要重算一次最近距离,时间复杂度略高,但好处是逻辑非常清楚,不容易出错。场景数在1万以内、目标K在50以内时,性能不是瓶颈。

4.3 削减效果到底怎么样:一组典型对比数据

用上面两个函数跑5000个场景削减到20个时,我记录过一组典型数据,可以给大家一个直观参考。不同站点数据、不同分布参数下结果会有差异,但量级基本一致。

从这张表能看出两件事:前向选择在保留尾部特征上确实比回代消除好一点;而无论哪种方法,削减到20个场景后,优化模型的计算规模都在工程可接受的范围内。

注意:均值误差控制在1%左右,不等于优化结果误差也在1%左右。随机优化关心的是累积分布尾部的表现,尤其是系统最不利时的成本和功率平衡,所以一定要看分位数误差,不能只看均值。

5. 实际项目中绕不开的五个坑:参数拟合、时序相关性、削减数量与验证流程

5.1 分布参数别照抄文献,先用站点数据拟合

这是很多初学者最容易犯的错。文献里写了Weibull形状参数k取2.0,尺度参数c取8,就照搬进自己的模型。但每个站点的风资源特性差很多——同一天内陆的强对流天气和沿海的稳定海风,在Weibull参数上的体现完全不同。

正确的做法是拉取站点至少一年的历史风速数据,用MATLAB的fitdist函数做最大似然估计:

pd = fitdist(wind_history, 'Weibull'); k = pd.A; % 形状参数 c = pd.B; % 尺度参数

光照的Beta分布参数也可以用类似思路,先从历史辐照度数据算出每个时段的均值和方差,再用矩估计反推α和β。拟合用的是"实测数据",而不是"经验参数",生成出来的场景才具备真实的统计特征。

5.2 时序相关性缺失时,场景削减会把"噪声"当特征

前面提到过,如果模型对每个时段独立抽样,生成的风速场景会像锯齿一样剧烈跳动:上午10点风速8m/s,11点突然掉到1m/s,12点又飙到9m/s。这种剧烈爬坡在实际风电场很少出现,但蒙特卡洛抽样会频繁产生。

问题是,场景削减算法并不认识"这是不合理时序",它只按距离度量选代表场景。如果原始场景集里充满了锯齿状噪声,削减出来的保留场景依然带锯齿,后面做机组组合时,爬坡约束会松得离谱,储能调度策略也会被带偏。

我常用的解决办法是给风速过程加一个AR(1)自回归结构:

v_t = μ_t + φ * (v_{t-1} - μ_{t-1}) + ε_t

其中μ_t是t时段的平均风速,φ是一阶自回归系数(通常取0.7到0.95),ε_t是Weibull噪声。这样生成的风速曲线既有统计分布特征,又保持了时间上的连续性。光伏出力则更简单——夜间出力本来都是0,白天按辐照度Beta分布走,时序相关性主要来自太阳几何位置的变化,只要把时段平均辐照度曲线建模好,问题不大。

5.3 削减数量别拍脑袋,用拐点法加统计校验

削减到多少个场景合适?我见过有人固定用10个,有人用50个,都说"凭经验"。更稳妥的做法是画一条拐点曲线:横轴是K,纵轴是削减后的Kantorovich距离(也就是削减带来的信息损失)。随着K从1增加到100,Kantorovich距离通常是先快速下降,然后趋于平缓。曲线的"拐弯处"就是均衡点——再增加场景,信息损失改善很小,但计算量继续上升。

具体操作是:先用前向选择分别算K=5、10、15、20、30、50时的最小概率距离,画出曲线,找一个斜率开始明显变小的K值。我那个项目里,拐点大概在K=15附近,为了保险起见选了20。如果模型规模允许,选拐点值的1.5倍左右比较稳妥。

5.4 距离度量里的量纲问题:风速和光功率不能直接比

场景矩阵是48维的——24个风电时段加24个光伏时段。直接对这48维算欧氏距离,风速维度(单位m/s)和光伏功率维度(单位MW)不在一个数量级。比如风速数值在0到15之间波动,光伏功率在0到0.5之间波动,距离计算会被风速部分主导,光伏的随机特性几乎被淹没。

解决方法是在构造场景矩阵时先做标准化。我习惯对风电和光伏出力分别做z-score标准化,再拼成场景矩阵做削减。标准化带来的问题是削减后的场景概率分布会变化,所以在统计校验环节要回到原始物理量纲下重新评估误差。

如果不做标准化,至少也要给风电段和光伏段分别加权重。这个细节看着小,但对削减质量的影响非常大。我最初没处理量纲问题的时候,削减出来的20个光伏场景几乎长得一模一样,等于浪费了10个场景名额。

5.5 削减之后必须做的统计校验

最后一步,也是我每次必做的一步,是对削减结果的统计校验。校验内容包括三项:

• 削减前后各时段风/光出力的均值误差,控制在2%以内
• 削减前后各时段风/光出力的标准差误差,控制在5%以内
• 削减前后总出力的5%、50%、95%分位数偏差,分别记录对比

如果某项误差超标,就要回到参数估计或削减算法里找原因,而不是硬着头皮往下算。通常问题出在分布参数拟合不准或者削减数量太少。做过一轮校验之后,优化结果才具备可信度。

最后再分享一个我固定的操作习惯:削减完成并校验通过后,我不会直接把削减场景集扔进优化模型,而是先用它跑一次单场景确定性优化,再把几个典型的极端场景(如低风低光、高风低光)单独拿出来跑一遍,对比两种结果对系统备用容量的要求。如果极端场景的结果明显更恶劣,就说明削减场景虽然统计上达标,但还是漏掉了对决策影响大的尾部情况,这时候我会适当增加K值或者调整距离权重,重新做一轮削减。这套流程走下来,优化方案才真正让人放心。