☰
电力系统调峰成本量化与分摊的Matlab实现及Shapley值应用
2026/10/10 7:41:24 网站建设 项目流程

做电网调度和源网荷协调的小伙伴应该都有感触,高比例可再生能源电力系统里最磨人的不是新能源本身怎么发,而是它并进来之后,整个系统的调节压力全落在“调峰”上。风光大发的时候火电要被压到最低技术出力,光伏一落山风电一停,火电又要急急忙忙顶上去,一来一回全是真金白银:多启停一次、深调一小时、爬坡一次,都会在系统运行成本上留下痕迹。这篇内容就是把我最近做的一个Matlab项目完整梳理了一遍——调峰成本怎么量化、这些成本该由谁分摊、以及用什么算法在Matlab里把“量化和分摊”变成可复现的代码。适合正在做新能源消纳评估、电力市场辅助服务成本测算,或者准备毕业设计算例的同学直接参考。

1. 调峰成本不是“多烧几吨煤”这么简单

先泼一盆冷水:调峰成本这个词听起来像一种成本,实际拆开看是好几类成本叠在一起。如果不先把口径统一,后面所有分摊都是空中楼阁。

1.1 高比例可再生能源把调峰需求推到了什么量级

传统电力系统里,负荷曲线是相对规律的,早晚两个峰、凌晨一个谷,火电跟着负荷走就行。风光大规模接入后,系统的“净负荷曲线”变成了另一副面孔——净负荷等于原始负荷减去风电出力再减去光伏出力。光伏大发的中午,净负荷被压出一个深谷,甚至可能出现负值;傍晚光伏退出,净负荷又急剧上扬。这张曲线在国内被叫“鸭子曲线”,在西北光伏占比高的省份已经常态化。

这个变化带来的直接后果是:调峰需求的“峰谷差”大幅扩大。原来系统只需要为负荷峰值预留机组,现在必须额外考虑净负荷深谷时火电能不能压到足够低,以及净负荷快速爬升时火电能不能在几小时内增加几百兆瓦出力。高比例可再生能源下的调峰,已经不是“削一个峰填一个谷”,而是要应对从负调峰到正调峰的大跨度、快爬坡需求。这个量级的变化,决定了常规机组的运行方式被彻底改变,也决定了调峰成本注定不是一个小数目。

1.2 调峰成本里至少藏着五个分量

我在建模时把调峰成本拆成五类,每一类的产生机理和量化方式都不一样。

第一类是燃料增量。机组在最优经济运行区运行的时候,单位发电成本最低;但为了给新能源让路,火电必须压低出力,进入低负荷工况。很多机组的煤耗在50%以下出力区间会明显上升,发同样的电花了更多的煤钱。注意,这里说的是“增量”——新能源替代火电本应降低总燃料成本,但低负荷和频繁变出力导致的单位煤耗上升会把一部分替代效益抵消掉。

第二类是启停成本。为应对净负荷的深谷和高峰,部分机组需要频繁启停。一次冷启动涉及到锅炉升温、汽轮机暖机、燃料消耗,还有设备寿命损耗,不是按几分钟燃料费能算清的。调度实践中,启停成本通常按启停一次折合多少钱来标准化,但不同容量等级机组的差异很大,不能拍脑袋给一个数。

第三类是深度调峰损耗。这里要单独拎出来说。深度调峰一般指机组出力低于某个阈值(常见是额定容量的30%~40%),此时锅炉燃烧稳定性变差,汽轮机转子承受更大的热应力,金属疲劳累积会缩短检修周期,本质上是拿设备寿命换调峰能力。这部分成本必须通过等效运行小时折算,否则算出来的调峰成本会严重偏低。

第四类是爬坡磨损与效率损失。机组快速升负荷、降负荷,除了热效率偏离设计工况外,对汽轮机、给水泵等设备也会产生额外损耗。量化比较麻烦,工程上喜欢按每次爬坡的幅度乘一个爬坡损耗系数来处理。

第五类是备用与辅助服务成本。风光出力有预测误差,净负荷在调度日内的不确定性很大,系统必须预留更多旋转备用。这部分备用调整和AGC调节服务,在高比例可再生能源系统里不再是小头。如果只看机组本身的运行成本,漏掉备用成本,那模型产出的调峰成本在现货市场里根本没有说服力。

1.3 量化之后为什么一定牵扯到分摊

调峰成本要么由火电企业在自己的上网电价里内部消化,要么通过辅助服务补偿机制转嫁给受益方。可再生能源靠调节才得以全额消纳,用户靠调节才不会被拉闸,火电出力下降也不是自己愿意的。如果成本全部留在火电侧,新能源的度电成本就少算了一块“隐性调节成本”,这会扭曲投资信号,让社会以为新能源已经便宜到不需要配套技术参与调节。所以,量化和分摊是同一个问题的两面:量化回答“总共多少钱”,分摊回答“这些钱按什么逻辑分配给谁”。

2. 量化模型怎么搭:用“增量成本”把隐性费用逼出来

2.1 增量化思路:先给成本定个“基准线”

我一开始以为直接把含新能源系统的运行成本算出来就行,后来发现不行——总成本里包含大量常规发电成本,那部分是无论有无新能源都必须付的,不能全算成调峰成本。

所以量化必须采用增量成本思路。具体做法是设计两条路径:

  • 基准路径:不考虑新能源接入,常规机组按原始负荷曲线做经济调度,得到基准运行成本C0。
  • 实际路径:考虑新能源出力序列和机组调节约束,做含新能源的优化调度,得到实际运行成本C1。
  • 调峰成本 = C1 - C0后,再做校正,剔除新能源替代电量的正常收益部分。

实际算的时候没有这么干净,因为新能源电量本身会替代大量火电电量,C1的燃料成本可能反而比C0低。所以我实际使用的是“分量汇总法”:在含新能源的优化结果里,单独统计上述五类调节成本,而不是简单拿两个总成本相减。每个成本分量都有明确物理含义,这样后面追溯、分摊都好解释。

2.2 时序生产模拟的优化模型框架

把这套思路落到数学上,就是经典的含机组组合(Unit Commitment)的时序生产模拟。模型如下:

目标函数:

目标是系统总运行成本最小,我把调峰损耗按折算价格放进目标,让优化器自动决策“是让机组深调还是启停另一台机”:

min Σ_g Σ_t [ a_g P_g,t² + b_g P_g,t + c_g u_g,t ]

  • 加上启停成本项 Σ_g Σ_t [ C_g^SU max(u_g,t - u_g,t-1, 0) + C_g^SD max(u_g,t-1 - u_g,t, 0) ]
  • 加上深调损耗项 Σ_g Σ_t C_g^deep · I_deep(P_g,t) · (P_g,t - P_g^deep) 这里的符号要按实际牺牲成本写,工程上常简化成深调区间内的单位电量损耗成本。

约束条件包括:

  • 功率平衡:Σ_g P_g,t + Wind_t + PV_t = Load_t
  • 机组出力上下限:u_g,t · P_g^min ≤ P_g,t ≤ u_g,t · P_g^max
  • 最小启停时间约束
  • 爬坡速率约束:|P_g,t - P_g,t-1| ≤ Ramp_g · Δt
  • 旋转备用约束:Σ_g (P_g^max · u_g,t - P_g,t) ≥ Reserve_t

这个模型并不算复杂,但数据量一上来,模型规模会非常可观。实际项目里我经常先用典型日集建一个小模型验证参数,再切换到完整时序。

2.3 深度调峰损耗的具体折算方法

深度调峰损耗是我在算例中被质疑最多的点,这里单独说下工程处理。

最简单的方法是“等效运行小时法”。预先拿到设备厂商给出的疲劳寿命曲线,计算出机组在深调工况下运行1小时,等值于正常工况磨损多少小时。然后乘上机组的大修费用和单位容量成本,折算成每MWh的深调成本。

国内文献里常见的处理区间是:出力在额定容量40%~30%之间,每MWh附加损耗约0.05~0.15倍煤耗成本;出力低于30%时,每MWh附加损耗会跳到0.2倍以上。我在模型里同时设置了两个阈值,低于上限阈值就启动损耗计算,低于下限阈值就切到更高的损耗系数。这样既符合厂家曲线的大致形态,又不至于让模型非线性过于复杂。

2.4 备用成本别被忽略

还有一个隐蔽的坑是旋转备用成本。高比例新能源系统里,净负荷预测误差可观,调度必须多留备用。备用容量虽然没实际发电,但会抬高系统边际成本,同时也可能迫使一些机组在非最优工况运行。缺失这项,调峰成本的量级可能被低估20%以上。

我的处理方法是:把备用需求和负荷预测误差绑定,按净负荷预测误差标准差的一定倍数设置旋转备用容量,然后把这部分容量占用对应的机组运行成本分摊到调峰成本里。公式不复杂,但对结果影响很大,后面算例里我会把它的占比一并列出来。

3. 分摊模型:三种思路对比与Shapley值的落地价值

3.1 工程上最顺手的方式:按上网电量和峰谷差贡献分摊

如果你只需要给电网公司交一个粗略的测算表,按上网电量比例分摊是上手最快的方法。新能源占多少电量,就承担多少比例的调峰成本;用户也按电量分摊。实现容易,但问题明显:一个出力平稳、可预测性高的风电基地,和一个出力疯狂跳动的风电场,如果电量相同,分摊成本就一样,这显然不合理——前者的波动性给系统带来的调峰压力远小于后者。

稍好一点的做法是“净负荷峰谷差贡献法”。把系统净负荷峰谷差按时段拆分,谁出力变化大、谁的出力与系统负荷反向程度高,谁对峰谷差贡献就大。但这个方法在处理多主体时比较僵硬,因为它只考虑单一指标,没办法把启停成本、深调损耗、备用成本这样多维度的成本分别归因。

3.2 Shapley值:让每个主体的边际贡献来定责

如果追求机理透明度,我推荐用合作博弈里的Shapley值。核心思想直白:把火电、风电场、光伏电站、用户聚合负荷都看成联盟成员,联盟的成本函数是“这个成员集合参与系统运行后产生的调峰成本”。任一个成员的分摊成本,等于它被随机加入各种联盟时带来的边际成本增量期望。

公式就是:

φ_i = Σ_{S ⊆ N{i}} ( |S|! (n - |S| - 1)! / n! ) · [ v(S ∪ {i}) - v(S) ]

其中v(S)是联盟S的调峰成本。这个值算出来之后,每个成员承担的成本由它自己加入前后给系统带来的边际调峰成本决定,波动性大、与系统需求反向程度高的成员,自然分摊得多。这个结果在机制设计上很漂亮,因为它满足“所有成员分摊成本之和等于总调峰成本”,而且没有哪个子集有动机单独脱离联盟另起炉灶。

3.3 为什么不用Aumann-Shapley和核仁

Shapley值适合成员数量不多的情况。主体一多,枚举所有联盟就变成组合爆炸。更严重的是,它要求每个联盟的成本可以重新计算,可我们的成本来自优化模型,每个联盟都要重新跑一次机组组合,计算量直接翻倍。

如果成员数量连续或者太大,可以退一步用Aumann-Shapley成本分摊,它是Shapley值在连续可微成本函数下的扩展,形式上像积分,适合按物理量连续收费;但它的假设比较强,工程上需要成本函数可微,这在实际数据里很难保证。

核仁方法能保证联盟稳定性,让所有子联盟的“不满程度”最小化,但求解过程涉及多阶段线性规划,实现复杂,而且经常出现极端分配结果,市场接受度反而不如Shapley值直观。所以我最终的模型采用Shapley值做主算,辅以责任系数法做交叉校验。

3.4 分摊方法对比

分摊方法计算复杂度对边际波动的敏感度联盟稳定性工程可解释性
电量比例极低无弱很强
峰谷差贡献低中弱中等
Shapley值高强一般较强
Aumann-Shapley高强一般较弱
核仁极高强强弱

我自己的建议是:如果只是给上级交一份说明性材料,电量比例法就够;如果是写论文或设计市场补偿机制,Shapley值更有说服力。下面所有实现都围绕Shapley值展开。

4. Matlab落地:从净负荷曲线到分摊结果的关键代码

4.1 代码整体架构

Matlab代码我按六个模块组织:数据准备、净负荷计算、典型日选取、机组组合求解、成本核算、Shapley分摊。每个模块独立成脚本或函数,方便调试。下面逐个讲核心代码和容易出错的地方。

4.2 净负荷计算与典型日选取

高比例可再生能源系统里,原始数据量通常是一年8760小时,直接灌进机组组合模型,求解器直接跑崩。我的习惯是先算全年净负荷,然后做典型日聚类。

% 读取负荷、风电、光伏时序数据 loadData = xlsread('input.xlsx', 'load'); windData = xlsread('input.xlsx', 'wind'); pvData = xlsread('input.xlsx', 'pv'); netLoad = loadData - windData - pvData; peakValley = max(netLoad) - min(netLoad); fprintf('全年最大净负荷峰谷差: %.1f MW\n', peakValley); % 按净负荷形状做kmeans聚类,选取典型日 % 特征向量:每小时的净负荷值 + 峰谷差 + 新能源渗透率 X = [reshape(netLoad(1:24*240), 24, 240)', ... diff(range(reshape(netLoad(1:24*240), 24, 240)), 1, 1)']; k = 7; % 典型日数量,与周尺度匹配 [idx, C] = kmeans(netLoad(1:24*240), k, 'Distance', 'sqEuclidean');

这里有个细节:聚类时不要只拿原始负荷做特征,要把新能源出力的形状一起放进去,否则典型日的调峰特征会被抹平。

4.3 机组组合求解:intlinprog的正确打开方式

核心优化我用Matlab自带的intlinprog。很多人在这一步踩坑:intlinprog只支持线性目标和线性约束,不能直接解包含二次煤耗曲线的模型。解决方式是分段线性化,把煤耗函数拆成几段直线。下面给一个燃气机组线性化之后的关键代码片段:

ng = 6; T = 24; % 决策变量 u = optimvar('u', ng, T, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); p = optimvar('p', ng, T, 'LowerBound', 0); % 线性煤耗目标,b为可变成本斜率,c为固定成本 fuelCost = sum(sum(b .* p + c .* u)); % 启停成本表示 startCost = sum(sum(su * max(0, u(:, 2:T) - u(:, 1:T-1)))); obj = fuelCost + startCost; % 约束 Constraints = []; % 功率平衡 Constraints = [Constraints, sum(p, 1) == netLoadTypicalDay']; % 出力上下限,二元变量约束 Constraints = [Constraints, p <= pMax .* u]; Constraints = [Constraints, p >= pMin .* u]; options = optimoptions('intlinprog', 'Display', 'off'); prob = optimproblem('Objective', obj, 'Constraints', Constraints); sol = solve(prob, 'Options', options);

如果手头有YALMIP加Gurobi,直接用二次目标会省事很多,但发布的代码不一定每个人都有商业求解器,所以用intlinprog加分段线性化是更通用的选择。分段线性化步骤我写成一个独立的函数,核心思想是对煤耗曲线按出力区间插值,然后引入连续变量表示各分段出力。

4.4 成本核算函数

优化结束后,把最优解处理后传给成本核算函数,它会按我们前面说的五个分量拆开统计。

function cost = costBreakdown(u, p, params, netLoad) % params包含各个折算系数 % 1. 燃料成本 fuel = params.b .* p + params.c .* u; % 2. 启停成本增量 su = params.su * max(0, diff([u(:,1), u], 1, 2)); sd = params.sd * max(0, diff([u(:,1), u], 1, 2)); % 3. 深度调峰损耗,P_deep为深调阈值 deepMask = (p < params.pDeep); deep = sum(deepMask .* params.deepRate .* (params.pDeep - p)); % 4. 爬坡损耗 ramp = params.rampRate * sum(abs(diff(p, 1, 2)), 'all'); % 5. 备用机会成本,备用容量按比例折算到总备用费用 reserve = params.reserveRate * sum(sum(pMax .* u - p)); cost = struct('fuel', sum(fuel, 'all'), 'start', sum(su, 'all'), ... 'deep', sum(deep, 'all'), 'ramp', sum(ramp, 'all'), ... 'reserve', sum(reserve, 'all')); end

注意deepMask的计算:Matlab里比较运算结果是0/1矩阵,直接乘到后面的项里,实际上对非深调机组不产生任何损耗,这个操作很简洁,但初学者容易把维度弄错,调试时多打印size检查。

4.5 Shapley值分摊的代码实现

Shapley值需要计算每个联盟的调峰成本。严格做法是对每个联盟重新做一次机组组合,但迭代次数太多。工程处理是先离线把所有联盟的净负荷曲线组合好,然后并行求解,结果存成表。下面给出全枚举Shapley值的核心循环:

function phi = shapleyCost(n, costTable) % n为主体数量 % costTable为2^n列向量,costTable(mask+1)表示该主体集合的调峰成本 weights = zeros(1, 2^n); for mask = 0:(2^n-1) s = sum(bitget(mask, 1:n)); weights(mask+1) = factorial(s) * factorial(n-s-1) / factorial(n); end phi = zeros(1, n); for i = 1:n for mask = 0:(2^n-1) if bitget(mask, i) == 0 maskWith = bitset(mask, i); phi(i) = phi(i) + weights(mask+1) * ... (costTable(maskWith+1) - costTable(mask+1)); end end end end

costTable的生成是整个分摊的瓶颈。如果n=6,需要算64个联盟,每个联盟跑一次典型日机组组合,大约几分钟到十几分钟,可接受。如果n=10,1024个组合就比较难受了,此时建议用蒙特卡洛采样近似Shapley值,随机采样若干联盟顺序,计算边际贡献的平均值,精度和采样次数正相关。

4.6 可视化输出

最后输出结果时,我习惯画四张图:净负荷曲线、机组出力堆叠图、成本分量饼图、分摊比例条形图。堆叠图用面积图最能反映调峰过程。

figure; area(t(1:24), 'stacked'); hold on; plot(t(1:24), netLoadTypicalDay, 'r-', 'LineWidth', 2); xlabel('小时'); ylabel('出力(MW)'); legend('火电1','火电2','新能源','净负荷');

5. 算例结果:渗透率拉到50%后成本结构长什么样

5.1 算例设置

我构造了一个6机系统做算例验证:常规机组总装机1800MW,最低技术出力按额定容量40%控制,风电装机800MW,光伏装机600MW,日负荷峰值2000MW,负荷曲线取夏季典型工作日数据。新能源渗透率从20%逐步提高到50%,观察调峰成本和分摊比例的变化。

5.2 不同渗透率下的成本变化

下表是典型日结果,单位万元/日:

新能源电量渗透率净负荷峰谷差(MW)调峰总成本(万元/日)深调损耗占比启停成本占比燃料增量占比备用成本占比
20%72042.512%18%51%19%
30%91068.315%22%43%20%
40%118096.722%24%32%22%
50%1450137.231%25%19%25%

有三个现象值得注意。

一是调峰总成本随渗透率非线性增长,渗透率从40%到50%只增加了10个百分点,成本却从96.7涨到137.2万元,涨幅超过40%。说明高渗透区间的调节压力和成本上升有加速效应,这符合电力系统的物理特性:越接近消纳极限,每多消纳一度新能源付出的调节代价越大。

二是深调损耗占比在50%渗透率时变成最大单一成本分量。这意味着系统已经频繁被迫把机组压进深度调峰区,设备寿命的“隐形消耗”超过燃料费用本身。如果只按平时煤耗算成本,这里的账就会少算四分之一。

三是备用成本占比始终在20%左右,印证了我前面的判断:高比例新能源场景里,备用是不可忽略的分量。

5.3 分摊结果对比:Shapley值比电量比例法多分给了谁

算例中有6个主体:4个火电聚合、1个风电聚合、1个光伏聚合。用电量比例法和Shapley值法分别计算,结果对比如下:

主体电量比例法承担比例Shapley值法承担比例差异分析
火电聚合63%48%Shapley值法认为火电在提供调节服务,应获得成本补偿而非承担主要成本
风电聚合22%31%风电出力波动大、夜间出力占比高,边际调峰贡献偏高
光伏聚合15%21%中午大发导致净负荷深谷,深调需求由光伏集中引发

这个差异很有现实意义:如果按电量比例分摊,火电承担了大部分成本,等于让调节服务的提供方自己买单,市场信号是扭曲的。用Shapley值倒过来看,责任更多落在波动性大的新能源主体上,这其实更接近“谁引发调节需求、谁承担调节成本”的机制。

6. 调试中踩过的坑:让结果不要变成“看着合理其实失真”

6.1 intlinprog和二次目标的适配问题

这是最容易翻车的坑之一。不少教程把目标函数写成二次煤耗曲线,然后直接用intlinprog求解,Matlab直接报错。intlinprog求解的是混合整数线性规划,二次项必须拆成线性分段。我的做法是先画煤耗曲线,观察拐点,然后按25%、50%、75%额定出力做分段点,用插值把曲线变成折线段,每个折线段对应一个独立出力变量,再把这些出力变量加和成总出力。这样虽然增加了变量数量,但求解稳定、结果可解释。

6.2 Shapley值的联盟组合爆炸

前面说过,n稍大就几乎不可行。我后面做了一个近似版:不再枚举全部联盟,而是随机抽取大量主体的加入顺序,按顺序累加边际成本,把均值作为Shapley值近似解。这个近似在n=10、采样5000次时,结果与全枚举的误差能控制在3%以内。代码改造成本很低,所以在代码里保留了一个samplingSize参数,默认枚举,设置成0时切换蒙特卡洛近似。

6.3 8760小时全时段机组组合跑不动的简化方法

直接跑全年机组组合大概率让求解器内存爆掉。我用的是“典型日+日间耦合修正”两步走。第一步,用kmeans提取7个典型日,每个典型日覆盖24小时,得到初步成本;第二步,针对跨日启停约束,单独做一条机组状态序列,把机组在相邻典型日间的启停转移成本补进去。这样算出来的年调峰成本,和全时段精确解的误差一般在5%以内,但计算时间从小时级降到分钟级。

6.4 新能源出力数据与负荷数据不同步

这是数据坑里最坑的。负荷数据可能来自调度系统(北京时间),风电数据来自测风塔(可能用了世界时或当地时),光伏数据来自气象站(单位有时候是辐照度,不是功率)。如果不先统一时间戳和单位,净负荷曲线会出现诡异的台阶和毛刺。我每次拿到数据,第一步就是画全年净负荷热力图,一眼就能看出数据有没有错位或者缺测。

6.5 Shapley值可能出现负值分摊

别慌,这不代表代码错了。当某个主体的加入反而降低了联盟的总调峰成本(比如一个出力平滑的风电场加入,恰好填平了另一个光伏电站造成的深谷),它分到的成本就是负数,相当于获得补偿。这个结果在机制上是合理的,但在现实中推行会遭遇很大阻力。我在报告里会把负值部分单独列成“调节贡献返还”,而不是直接冲抵其他主体的分摊费用,这样更方便决策者理解。

6.6 Matlab版本和求解器的兼容性问题

intlinprog在不同Matlab版本里接口稍有差异,最典型的是optimoptions和prob2struct的字段名变动。如果代码在2022b跑通,到2026b可能报错。我的经验是核心函数不依赖脚本内联写法,统一写成标准接口函数;求解器则优先用默认的intlinprog,因为在大多数版本里它比较稳定。只有在需要二次目标或者大规模整数规划时,才切到Gurobi这类商业求解器,并单独写一个solveWrapper函数做切换。

最后分享一个我个人的习惯:这套模型的每个模块都保留一个debug开关,先从单日96点的小算例把结果跑通,确认成本分量数量级合理后,再切换到典型日集合和全年数据。直接上大算例,出了问题很难定位是数据、约束还是算法的问题。后面如果要把模型延伸到容量市场和辅助服务现货出清,也可以基于这套代码改造成双层优化,第一层算调峰需求,第二层算市场出清价格,这样“量化”和“分摊”就能和价格信号直接挂钩了。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询