高比例可再生能源电力系统调峰成本量化与分摊：模型、Matlab实现与Shapley值应用

“这600万到底贵在哪”两年前我第一次接触“高比例可再生能源电力系统的调峰成本量化与分摊模型”这个课题导师看完我第一版结果只问了这么一句。我当时手里的数据是算出了带新能源和不带新能源两套运行成本差值600万但我既说不清这600万里有多少是机组启停造成的、多少是深度调峰煤耗上升造成的也说不清这笔钱该记在风电头上还是光伏头上还是该由火电自己承担。那一刻我才意识到调峰成本不只是一个统计口径问题它本质上是一条完整的链路先把“调峰到底动用了哪些资源”拆开再把每一项成本量化进优化模型最后按博弈规则或物理责任分摊到各个市场主体身上。这篇博文就把我这套流程完整写出来调峰成本的构成与数学建模、MatlabYalmip实现机组组合与经济调度的量化代码、基于调峰深度/调峰电量/Shapley值的三种分摊方法以及我做敏感性分析和结果可视化时踩过的坑。适合电力系统优化方向的研究生、做新能源并网评估的工程师如果你想在一个标准24节点或IEEE 30节点系统上复现“可再生能源渗透率提高→调峰成本上升→谁该承担多少”这个分析链条这篇文章可以直接当脚手架用。1. 高比例可再生能源的“净负荷陷阱”为什么传统成本口径失效1.1 从负荷曲线到净负荷曲线调峰压力的真正来源先明确一个核心概念我们常说的“负荷每天有峰有谷”其实只是传统电网的情况。当风电、光伏大规模接入后系统实际要平衡的对象不再是负荷曲线而是净负荷曲线净负荷 负荷 − 风光出力风电出力的随机性和光伏出力的昼夜特性叠加在负荷上之后净负荷曲线的形态会比原始负荷更加极端。最典型的是光伏大发的中午时段负荷还在高位但净负荷被压低甚至出现负值而傍晚光伏快速退坡、负荷继续攀升净负荷又急剧抬高。这个过程在业内有个形象的说法叫“鸭子曲线”。我做过一组典型日数据对比可以很直观地展示这个问题时刻原始负荷(MW)风光出力(MW)净负荷(MW)08:0072001200600013:0080005000300019:00105001800870023:0065008005700原始负荷的峰谷差只有4000MW左右但净负荷峰谷差可以拉到5700MW而且13:00到19:00之间的爬坡速率远超传统负荷变化速率。这就意味着系统里必须有足够多的机组能够快速降出力、快速提出力、频繁启停——这件事本身就值钱并不是免费的。1.2 传统电能量成本为什么不能反映调峰传统经济调度模型的成本口径非常单一只算机组的燃料耗量成本谁发的电多谁花的煤钱就多。这个口径在基荷运行的场景下问题不大到了高比例可再生能源场景就失效了原因有三第一机组为了给新能源让路被迫压到很低的技术出力以下运行甚至进入“深度调峰”区间。在这个区间里锅炉热效率明显下降单位电量的煤耗率不是线性上升而是曲线上升。这部分其实是“额外成本”。第二启停一台大型火电机组绝不是按一下开关那么简单。从热态启动到并网要烧掉大量燃料停机过程还有管道工质损耗和热应力导致的寿命损耗。传统只按发电量算成本就相当于把出租车司机接单送人的油钱算了但从空车赶去接人的油钱没算。第三为了应对净负荷的剧烈波动系统还要预留大量备用容量调用储能、需求响应等灵活资源。这些资源本身的折旧、损耗和补偿费用同样不在传统电量成本里体现。所以我在这套模型里采用的口径是增量成本法调峰成本 计及高比例可再生能源后的系统运行总成本包括各项灵活性成本 − 同一系统在基准场景下不考虑调峰需求的电量成本。这个差额才是真正由“调峰需求”引发的成本后期做分摊时也以这个增量为对象避免把本来就该发电的那部分煤耗也摊到新能源头上。2. 调峰成本的构成拆解与数学模型建立2.1 成本构成全景量化模型的第一步不是写代码而是先把“调峰成本”这个抽象概念拆成一张可以计量的账单。不同文献口径不同我在这套模型里选择了下面五个部分每一部分都对应一组明确的决策变量和成本系数成本项驱动因素模型中的体现方式机组启停成本为跟踪净负荷变化而启停机组启动事件变量 × 单次启动成本深度调峰煤耗增量机组在低负荷区间效率下降分段煤耗函数在低负荷段取更高斜率储能运行折旧成本储能循环充放、寿命损耗充放电功率 × 单循环折算成本弃风弃光惩罚成本调峰能力不足被迫削减新能源弃电变量 × 惩罚电价需求响应补偿成本调用柔性负荷参与调峰可调负荷功率 × 单位补偿价格这套构成里前两项是传统电源侧的调峰成本中间两项是新型灵活性资源成本最后一项是负荷侧成本。实际项目中可以根据数据情况增删比如加入碳排放惩罚、备用容量成本等但核心框架基本一致。2.2 火电调峰成本的关键建模启动、深调与爬坡火电机组是系统中最大的调峰资源也是最需要精细建模的对象。我总结了三个必须模块化的成本来源。**启动成本。**工程上常用指数形式描述机组启动后的热状态带来的启动燃料变化SUC_i a_i b_i × (1 − exp(−T_off_i / τ_i))其中T_off是机组停机时长τ_i是时间常数。简化到调度模型中通常会按“热启动/冷启动”分两档取固定值。我做的系统规模不大直接给每台机组一个线性化的启停成本系数也不影响趋势判断。**深度调峰煤耗增量。**传统经济调度默认机组出力在[P_min, P_max]区间内运行煤耗函数近似为F_i(P) a_i × u_i b_i × P_i c_i × P_i²但深度调峰意味着允许机组出力低于常规最小技术出力进入P_extreme到P_min的区间这个区间内煤耗曲线斜率明显变陡。做法是分段线性化每段用不同的斜率系数b在求解MILP时直接通过额外变量引入分段区间。这一步非常关键因为如果没有深调区间模型会倾向于让机组“不降出力”而多弃风从而低估系统的实际调峰成本。**爬坡成本/备用成本。**爬坡本身不消耗燃料但它挤占了机组的调节空间导致系统需要额外的向上/向下备用。这部分在机组组合模型中通常体现为旋转备用约束并预留一定比例的容量成本则反映在启停和煤耗的联动结果里不必单独设置惩罚项。2.3 新能源弃电、储能与需求响应的成本建模新能源机组的运行边际成本几乎为零但出现弃电时意味着系统“没消纳掉本来可以发的便宜电”。我用一个可调变量curtail表示弃电功率并在目标函数中加入惩罚项C_curtail λ × Σ curtailλ的数值取多少取决于想表达什么如果按市场化口径可以取新能源上网电价或者绿证价值如果按社会福利口径可以取一个足够大的惩罚系数保证模型先调用一切可调节资源、最后才弃电。储能模型在这里不是简单的能量守恒需要同时描述充放电功率限制、SOC递推和循环损耗成本。核心约束是SOC(t1) SOC(t) (η_ch × P_ch(t) − P_di(t) / η_dis) × Δt目标函数里的储能成本项用C_sto × (P_ch P_dis)近似相当于把折旧和损耗折算到单次充放电功率上。需要注意SOC可以直接取连续变量行为约束通过上下限实现充放电状态要用一个二进制变量避免同时充放。需求响应建模我用了最简单的可调负荷形式设定每个时段可削减/可平移的负荷量和对应的单位补偿价格然后让优化模型自己决定“调用多少柔性负荷去填谷”比“让火电深度调峰”更划算。2.4 完整的目标函数与约束集把上面的模块拼起来整条优化模型的骨架是这样的目标函数min Σ_i (启动成本 停机成本) Σ_i 煤耗成本 Σ弃电惩罚 Σ储能损耗成本 Σ需求响应补偿约束集合功率平衡Σ火电出力 Σ风电出力 − 弃电 Σ储能放电 − Σ储能充电 Σ可调负荷 原始负荷上下备用旋转备用 ≥ 负荷预测误差与新能源出力预测误差的覆盖需求机组出力上下限u … 需要区分常规区间和深调区间爬坡约束机组相邻时段出力变化量不超过爬坡速率最小启停时间约束储能SOC与充放电状态约束新能源弃电变量范围约束这是一个典型的混合整数线性规划MILP问题。变量规模大概是机组数 × 时段数 × 状态类型。用Matlab的Yalmip工具箱可以非常方便地把这套约束“翻译”成求解器能吃的模型接下来进入工程实现环节。3. MatlabYalmip工程化实现从公式到可跑通的代码3.1 环境准备与求解器选型我的开发环境是Matlab R2021b Yalmip工具箱 CPLEX 12.10求解器。如果你是新用户我建议至少装两个求解器CPLEX和Gurobi任选其一再留一个Gurobi姊妹版或SCIP作为免费备用。Yalmip本身不是求解器它只是建模语言层真正干苦力的是底层的MILP solver。这里有一个必须提前提醒的坑Yalmip需要匹配Matlab的Java版本和求解器的版本。我最早装的时候Matlab R2020b配新版Gurobi直接报“Unable to load DLL”原因是Gurobi的JAR包和Matlab自带JDK不兼容。解决办法是去MathWorks官方文档查到当前Matlab版本对应的JDK版本然后下载与该JDK匹配的Gurobi历史版本。别小看这一步它浪费了我整整一个下午。如果只需要做小规模验证或者不想折腾商业licenseMatlab自带的intlinprog也能解中小规模的MILP。24节点、6台机组、24时段的规模intlinprog硬解也能在两三分钟内出结果只是在大规模场景下速度明显不如CPLEX/Gurobi。我的习惯是模型调试阶段用intlinprog正式跑敏感性分析时切到Gurobi。3.2 输入数据结构设计建模前的数据准备我统一放在一个CaseData.m脚本里返回一个structfunction caseData CaseData() % 机组参数6台火电机组 caseData.gen.pmax [200; 200; 150; 150; 100; 80]; % MW caseData.gen.pmin [80; 80; 50; 50; 30; 20]; % MW caseData.gen.rampUp [40; 40; 30; 30; 25; 20]; % MW/h caseData.gen.rampDown [40; 40; 30; 30; 25; 20]; caseData.gen.a [100; 100; 80; 80; 60; 40]; % 空载成本 caseData.gen.b [0.15; 0.15; 0.18; 0.18; 0.22; 0.25]; % 增量煤耗系数 caseData.gen.startCost [5000; 5000; 3000; 3000; 2000; 1500]; caseData.gen.minUp [6; 6; 4; 4; 2; 2]; % 最小开机时间 caseData.gen.minDown [4; 4; 3; 3; 2; 2]; % 最小停机时间 caseData.gen.initStatus ones(6,1); % 初始运行状态 % 24h负荷曲线与风光出力曲线 caseData.load.P [ ... % 24行1列数据 5800 5600 5350 5200 5250 5400 5900 6800 ... 7800 8400 8700 8600 8100 7900 7800 7600 ... 7900 8500 8900 9200 9000 8200 7000 6100] ; caseData.re.wind [ ... % 24行1列风电数据 800 850 900 950 900 800 700 600 500 450 400 ... 380 360 350 400 500 600 650 700 720 680 750 780 810] ; caseData.re.solar [ ... 0 0 0 50 200 500 900 1200 1400 1500 1600 ... 1650 1600 1400 1100 800 500 200 50 0 0 0 0 0] ; caseData.re.curtailPenalty 50; % 元/MWh end数据组织的原则是“批量、向量化”所有机组参数尽量写成列向量避免在约束循环里内嵌太多标量变量。参数单位统一用MW和MWh成本用元方便后面直接出图和汇报。3.3 Yalmip建模与求解核心代码定义决策变量是整个建模里最容易出错的地方我用的命名规范是u表示机组启停状态、su/sd表示启动和停机事件、p表示出力、curtail表示弃电、soc表示储能电量、pch/pdis表示储能充放电功率。nG 6; T 24; u binvar(nG, T, full); su binvar(nG, T, full); sd binvar(nG, T, full); p sdpvar(nG, T, full); % 机组出力 % 新能源弃电变量 nRe 2; curtail sdpvar(nRe, T, full); % 储能变量可选 nS 1; soc sdpvar(nS, T1, full); pch sdpvar(nS, T, full); pdis sdpvar(nS, T, full); chState binvar(nS, T, full); % 1表示充电状态 Constraints []; % 功率平衡 Constraints [Constraints, sum(p,1) sum(pch,1)*(0) ... ... % 实际按 netLoad 写 ];写到这里我要特别强调一个逻辑功率平衡方程里火电出力、储能放电、新能源可用出力减去弃电量之和要等于净负荷加上充电功率。正确写法是把原始负荷、风光出力放在一起先算netLoad再让决策变量去平衡它netLoad caseData.load.P - caseData.re.wind - caseData.re.solar; % 功率平衡含储能 Constraints [Constraints, ... sum(p,1) sum(pdis,1) - sum(pch,1) sum(caseData.re.wind caseData.re.solar - curtail, 1) ... caseData.load.P ];机组状态逻辑约束% 启停状态一致性u(t) - u(t-1) su(t) - sd(t) for t 2:T Constraints [Constraints, u(:,t) - u(:,t-1) su(:,t) - sd(:,t)]; Constraints [Constraints, su(:,t) sd(:,t) 1]; end最小启停时间约束是最容易写错的。比较可靠的线性化写法是“事件触发式”如果在t时刻启动则接下来minUp个时段必须保持开机for t 1:T for i 1:nG if t minUp(i) - 1 T Constraints [Constraints, ... sum(u(i, t:min(tminUp(i)-1, T))) minUp(i) * su(i, t)]; end end end目标函数Objective sum(sum(caseData.gen.a .* u)) ... % 空载成本 sum(sum((caseData.gen.b) .* p)) ... % 边际煤耗成本 sum(sum(caseData.gen.startCost .* su)) ... % 启动成本 sum(sum(caseData.re.curtailPenalty * curtail)); % 弃电惩罚求解放入ops结构ops sdpsettings(solver, gurobi, verbose, 2, ... savesolveroutput, 1, solveroutputfile, out.json); % 设置MIP间隙防止卡在太严格的终止条件上 ops.gurobi.MIPGap 0.01; result optimize(Constraints, Objective, ops);3.4 求解过程中踩过的四个典型坑**坑一大M参数过大导致数值病态。**很多约束需要引入大M把非线性关系线性化比如“机组出力必须小于Pmax×u”。M值不要直接用99999那会让求解器在预处理阶段产生大量数值误差。我的做法是把M设置成对应约束里物理上限的1.1~1.5倍。比如Pmax是200MWM就取220既保证约束有效性又不会破坏LP松弛质量。**坑二二进制变量维度过高导致求解时间爆炸。**完整机组组合问题在数百节点系统里会有几千个二进制变量。如果直接对24h逐时段建模Gurobi可能在几分钟内出结果但如果扩展到8760h还保持逐时二进制建模内存就直接爆了。解决思路有两个一是用典型日聚类压缩时间维度二是把机组组合分成“日前开停机决策”和“实时经济调度”两个阶段第二阶段不再引入启停变量。**坑三SOC递推约束里的充放电同时性问题。**初次建模时只写SOC递推公式没有限制pch和pdis互斥结果目标函数为了最大化放电收益让储能在同一个时段既充电又放电相当于白赚了两倍的损耗成本。加入chState二元变量后约束变成Constraints [Constraints, pch chState .* chMax]; Constraints [Constraints, pdis (1 - chState) .* disMax];**坑四初始状态对结果影响巨大。**尤其是机组最小启停时间约束的累加段如果t1时刻的初始状态设定较随意优化结果会出现第一小时频繁启停的伪影。务必将initStatus、初始出力、初始SOC都从数据脚本里读出来而不是在建模脚本中写死为0。4. 成本分摊从总账到分账的博弈逻辑与实现4.1 为什么要做分摊以及公平性底线量化模型算出的调峰总成本是一个“总账”。但在实际项目中最后一定要回答“这笔钱谁来出”的问题是风电光伏企业按出力造成净负荷波动的程度补偿给火电还是火电因自身灵活性不足要买单又或者是用户因为用电时段分布不合理而承担调峰压力。分摊需要遵守三条底线第一是可解释性分摊结果要能讲清楚为什么某个主体承担这个比例第二是激励相容性不能让分摊机制反过来鼓励新能源故意不预测、或者鼓励火电故意把成本做大第三是可操作性分摊方法在现有结算和计量体系里要能落地不能公式漂亮但算不出来。4.2 基于调峰深度与调峰电量的简化分摊法最实用的简化方案是调峰深度法把每个调峰资源实际承担的调峰深度作为分摊权重。定义某台火电机组在时段t的调峰深度为它在基准出力P_base_i通常是Pmax附近与实际出力P_i(t)之间的差depth_i(t) P_base_i − P_i(t)系统总调峰深度是各机组之和。那么新能源需要承担的调峰成本份额可以按“系统净负荷峰谷差中被新能源出力波动放大的部分占比”来计算。这种方法计算量极小里面不含任何优化求解适合做快速评估和初步汇报。还有一种更细的调峰电量法按各主体在峰时段的“贡献电量”来分摊。峰时段里新能源出力越小、负荷越高说明调峰压力越大谁造成了这个状态谁就按比例承担。这个方法的好处是只依赖电量计量数据和结算系统天然兼容。但这两个方法都有一个共同的弱点它们只看物理量不看经济边际。可能出现的情况是某个新能源场站出力波动完全平滑但因为并网点位置刚好在负荷中心被“平均”分摊了调峰成本。想要更公平还得上互动博弈的方法。4.3 基于Shapley值的边际贡献分摊模型Shapley值来自合作博弈论思想很简单一个主体的贡献不取决于它单独做了什么而取决于它加入不同联盟时给联盟带来的边际变化。“谁让系统多花了钱谁就多承担”。把这套逻辑映射到调峰成本分摊上我把参与调峰的四个主体看成联盟成员火电机组群T、新能源群R、储能S、需求响应D。定义特征函数v(K) “只有联盟K中的资源参与调峰时系统需要承担的调峰增量成本”。注意这里所谓的“只有联盟K参与调峰”需要在优化模型里把联盟外资源的灵活性置零比如不允许新能源联盟外的机组参与深调、不允许联盟外的储能充放电、不允许联盟外的负荷做灵活调节。然后主体i的Shapley值为φ_i Σ_{S ⊆ N \ {i}} [ |S|! (n − |S| − 1)! / n! ] × ( v(S ∪ {i}) − v(S) )n4时只有16个联盟实际需要考虑的非空联盟15个每个联盟跑一次优化模型算力完全可控。Matlab实现框架如下nPlayers 4; % 1:火电 2:新能源 3:储能 4:需求响应 v zeros(2^nPlayers, 1); for s 1:2^nPlayers-1 idx find(bitget(s-1, 1:nPlayers)); % 联盟成员列表 v(s) computeCoalitionCost(idx); % 调用模型计算该联盟下调峰增量成本 end phi zeros(1, nPlayers); for i 1:nPlayers for s 1:2^nPlayers-1 % 只处理不包含i的联盟 if ~bitget(s-1, i) S find(bitget(s-1, 1:nPlayers)); % 联盟S∪{i}的索引 sIdx s 2^(i-1); w factorial(length(S)) * factorial(nPlayers - length(S) - 1) / factorial(nPlayers); phi(i) phi(i) w * (v(sIdx) - v(s)); end end endcomputeCoalitionCost函数内部就是第3节那套Yalmip模型只是把参数掩码传给CaseData屏蔽掉联盟外的灵活性资源。这一步是整个分摊模型的核心工作量所在——我建议为每个联盟单独生成一个CaseData副本而不是在建模脚本里写一堆if分支不然debug时会非常痛苦。Shapley值在这里体现的经济含义是某个主体的调峰成本责任等于它对系统总调峰成本的“平均边际贡献”。火电机组在联盟里贡献大量调峰能力如果它的存在显著降低了弃电惩罚它的Shapley值可能是正的“贡献值”新能源如果只增加波动而不提供灵活性它的值就显著为正意味着它应当承担较高成本份额。这个方法最大的价值是能给出一个完备且可解释的分摊结果缺点是模型求解次数多而且特征函数的计算口径需要提前与各方对齐。4.4 三种分摊结果的对比与使用场景分摊方法计算开销优点缺点适用场景调峰深度法极低物理意义清晰数据易得只看深度不看边际可能平均主义快速评估、政府对账、初步方案调峰电量法极低与结算系统兼容性好无法体现储能/DR的灵活性价值月/季度结算、场站补贴核算Shapley值法中15次求解公平性最强、可解释性高需要统一各方的特征函数口径机制设计、争议处理、年度清算我的实际建议是先用调峰深度法跑一版结果用于汇报再用Shapley值法做一版精细结果用于机制谈判或论文分析两者差异本身就是很好的研究素材。5. 结果可视化与模型扩展方向5.1 用Matlab画好五张关键图的实操模板模型跑通之后结果可视化决定了这份工作能不能“说服人”。我每次会固定输出五张图第一张是调度结果堆叠面积图横轴是24小时纵轴是功率依次堆叠火电出力、风电出力、光伏出力、储能放电/充电、负荷曲线。这张图能让评审一眼看到“谁在什么时候顶上了缺口”。figure; area(t_hours, [p_thermal p_wind p_solar p_dis], LineStyle, -); hold on, plot(t_hours, caseData.load.P, k-, LineWidth, 2); ylabel(功率/MW); xlabel(时刻/h); legend(火电,风电,光伏,储能放电,负荷, Location,best);第二张是净负荷与火电最小技术出力对比图把净负荷曲线和系统开机机组允许的最小出力之和画在一起两线之间的面积就是系统“必须动用额外调峰手段”的压力量化图这张图比任何文字都有说服力。第三张是成本构成柱状图把启停成本、深调煤耗增量、储能损耗、弃电惩罚、DR补偿分列为五个柱子观察新能源渗透率从20%升到40%时哪根柱子长得最快。第四张是分摊比例的横向条形图Shapley结果用barh绘制每个主体一条可以按承担份额降序排列。注意给图配上“基础运行成本不算在内图中仅统计调峰增量成本”的说明避免误读。第五张是敏感性分析曲线图横轴是可再生能源渗透率纵轴是单位调峰成本元/MWh画2~3条曲线分别对应不同峰谷差或燃料价格场景。5.2 敏感性分析设计渗透率、峰谷差与燃料价格我把敏感性分析当成分摊模型的“压力测试”来做。最常做的是三组第一组可再生能源渗透率10%、20%、30%、40%、50%看调峰总成本和单位调峰成本的变化趋势。你会发现在渗透率超过某个阈值后成本曲线会出现明显加速上升的拐点——这个拐点就是系统“免费消纳能力”的上限也是规划层面最值得关注的数字。第二组把原始负荷曲线的峰谷差按比例拉大或缩小模拟用电结构变化。峰谷差变大时通常火电深度调峰和储能的调用量同时上升两笔钱叠加调峰成本上升较明显。第三组燃料价格上下浮动20%观察分摊结果是否会出现“角色反转”。比如煤价高企时火电深度调峰的机会成本更高分摊到新能源头上的比例会明显加大。每组敏感性分析建议固定其他条件只改一个参数保证因果链路清晰。我把所有场景的优化求解和Shapley计算封装成一个两层循环外面跑频次里面跑联盟一个典型场景的全套计算在Gurobi下大约10分钟出完。5.3 从24h教学版到8760h工程版时序聚合与滚动求解第3节代码是标准的24h单时段断面模型适合跑通逻辑和汇报演示。但真实工程评估里一个“高比例可再生能源系统”需要看的不只是一天而是全年8760小时的运行效果。直接扩展会面临二进制变量爆炸问题我自己的过渡路线是第一步先做典型日聚类。用k-means对全年的净负荷曲线、风光出力曲线聚类得到春夏秋冬各一个或者六个典型日。这个做法能保留季节特性和天气模式同时把问题规模压到可控范围。第二步做周滚动优化。保留48小时重叠窗口每次只优化未来一周的机组组合滚动更新开机计划。Matlab实现时每次滚动求解都复用上一次末段的机组状态作为初始状态这个做法能自然处理跨日启停约束。第三步如果课题需要接现货市场可以把调峰成本模型改写成两阶段市场出清日前阶段做机组组合和启停实时阶段做经济调度和平衡成本计算。调峰成本量化模块保持不变只是从“事后统计”变成“机制设计工具”。做完这三步你已经从一个“调峰成本计算脚本”升级到一套可以支撑规划分析和政策研究的模型体系了。我在这条路上踩过很多坑写出来也是希望后来者能把精力花在系数标定和结果解读上而不是卡在Yalmip语法和大M参数这种本不该耽误时间的地方。