高比例可再生能源系统调峰成本量化与分摊模型Matlab实现

在电力系统分析领域待久了你会发现一个现象系统里真实的成本往往藏在那些没人愿意核算的“隐性代价”里。我去年参与一个区域电网的新能源并网评估项目调度同事拿着净负荷曲线问我“风电半夜满发的时候火电得压到50%以下出力多耗的煤、多磨损的机组、早上再爬坡的损耗谁来认账”这个问题直接催生了这套《高比例可再生能源电力系统的调峰成本量化与分摊模型Matlab代码实现》。今天这篇文章就是把这套模型的思路、公式、Matlab实现路线和调试中踩过的坑全部梳理出来供做电力系统规划、运行分析或者写论文的同学参考。这套东西能解决什么具体问题简单说当风电光伏占比升高以后系统为了平衡净负荷负荷减去新能源出力的剧烈波动需要火电频繁启停、深度压出力、快速爬坡这些动作都会产生额外成本。传统电力市场或调度核算里这部分成本要么被摊进平均电价要么直接没人算。模型要做两件事——先量化出这部分成本到底有多少再用一个可解释的分摊规则把成本分配到风电、光伏、负荷、火电各主体头上。代码层面基于Matlab实现核心求解用到linprog/intlinprog适合电力专业研究生、新能源企业分析工程师、设计院规划人员参考。1. 调峰成本从哪来高比例可再生能源改变了系统“调峰”的底层逻辑想准确量化一个东西先得搞清楚它产生的机理。调峰成本不是新能源接入以后才有的但确实是新能源比例高了以后才变得不可忽略的。这一节把底层逻辑理清楚后面照着公式建模才不会跑偏。1.1 传统电力系统怎么调峰边际成本排序与显性成本在传统火电为主、负荷曲线相对平滑的系统里调峰其实就是“跟着负荷走”。负荷白天高、夜间低调度员安排机组出力峰荷时段让效率高的机组多发谷荷时段让机组压低出力。这种调节产生的成本是燃料消耗的变化而燃料消耗本身可以通过煤耗特性曲线算出来——煤耗曲线上每一点的斜率就是边际煤耗率调度优化本质上是按边际成本由低到高去排机组带负荷。这套逻辑下调峰成本是显性的跟着燃料价格走谁多发谁多承担争议不大。但传统调峰有一个隐含前提峰谷差在系统可调能力范围内火电机组即便压低出力也还维持在一个经济工况附近。常规煤电机组的最小技术出力一般在40%~50%额定容量也就是说你让它压到50%出力它还能保持不太离谱的煤耗率。这时候的调峰成本大体可以用单位煤耗的升高来估算量级很小常常被归入“运行成本”而不单独核算。1.2 新能源接入后的“鸭子曲线”调峰需求从峰谷问题变成爬坡与低出力问题高比例风电光伏接入以后调峰问题的性质完全变了。净负荷曲线不再是一条平滑的“白天高、晚上低”而是出现典型的“鸭子曲线”白天光伏大发时净负荷被压出一个深谷傍晚光伏骤减、负荷攀升净负荷又急速拉高形成非常陡的爬坡段。这意味着系统需要的不再是“缓慢跟随负荷”的能力而是两种更苛刻的能力一是深度调峰能力。午间光伏大发时段净负荷可能跌到很低火电必须压到最小技术出力以下甚至进入不经济的低负荷运行区间。很多机组的深度调峰技术出力下限是30%左右低于这个值要么投油稳燃要么停机。二是快速爬坡能力。傍晚净负荷在2~3小时内可能陡增几千兆瓦常规火电爬坡速率往往跟不上调度只能提前多开几台机组“空转”备用这些机组在等待时段几乎不发电却仍然耗能。所以我一直跟人强调高比例可再生能源下的调峰本质上是“净负荷曲线的极端化”——峰谷差更大、爬坡更陡、不确定性更高。成本就藏在这些“更”字里。1.3 调峰成本的本质一项长期被隐性化、如今必须显性化的系统成本为什么说调峰成本必须单独量化因为它长期混在系统总运行成本里没人说得清新能源每发一度电到底给系统带来了多少额外调峰负担。举一个直观的例子某电网原本10台火电机组带负荷每台都在60%以上出力运行系统整体煤耗率很低。光伏接入后午间净负荷低到只需要6台机组出力另外4台要么停机要么压到深调状态——停机再启动有启停成本深度压出力时煤耗率急剧升高早高峰再让它们爬坡回来又产生爬坡损耗和寿命损耗。这些成本如果只算“系统总燃料费用”根本没法看出来是“谁”造成的一旦要分摊就必须把每一类的机理、每一笔的钱数拆开算。把这件事想明白以后量化模型的设计目标就很明确了把总运行成本里属于“调峰相关”的组成部分识别出来逐项公式化再通过一个公平的规则分摊到各个引起调峰需求的主体身上。下面两节分别讲这两步的具体做法。2. 调峰成本的量化模型四个成本项与两条实现路径量化模型的核心是“拆项”。我把调峰成本拆成四个主要组成部分每一部分都有明确的物理含义和计算公式。实际算例中可以按需选择计入哪些项但为了让结果有说服力建议四项都包含。2.1 四类必须计入的成本及其公式**第一类机组启停成本。**新能源出力波动导致的机组频繁启停是调峰成本里最直观的一块。计算公式为 [ C_{su} \sum_{t} \sum_{i} (SC_i \cdot y_{i,t}) ] 其中 (SC_i) 是机组i的单次启动成本(y_{i,t}) 是0-1变量表示机组i在时段t是否发生了“由停转启”的状态切换。停机成本同理用 (SD_i) 和切换变量计算。这里注意启停成本往往不是一个固定值冷态、温态、热态启动成本差别很大工程上多数模型会简化成与停机时长相关的分段函数但初版代码用固定值就够跑通流程。**第二类深度调峰运行成本。**机组压到最小技术出力以下运行时煤耗率显著上升甚至需要投入助燃燃料。简化计算是设定一个“经济出力下限”低于该下限的出力部分单位发电成本乘以一个惩罚系数。公式 [ C_{deep} \sum_{t}\sum_{i} \beta_i \cdot (P_{i,t}^{eco} - P_{i,t}) \cdot \Delta t, \quad P_{i,t} P_{i,t}^{eco} ] (\beta_i) 是深度调峰的单位惩罚成本元/MWh需要根据机组煤耗特性曲线拟合。**第三类爬坡损耗成本。**机组频繁大幅爬坡会加剧热应力影响寿命。工程上常用“爬坡里程”折算成本也就是每个时段出力变化的绝对值乘以一个单位爬坡损耗系数 [ C_{ramp} \sum_{t}\sum_{i} \lambda_i \cdot |P_{i,t} - P_{i,t-1}| ] 这个系数不容易从教科书上查到实际项目里多半参考厂家提供的寿命损耗曲线折算。没有资料时可以先按燃煤成本的2%~5%估算后面做敏感性分析时再看结果对系数敏不敏感。**第四类弃风弃光机会成本。**当系统调峰能力用尽时只能弃掉一部分新能源电量。这部分电量的“机会成本”可以按新能源上网电价扣除运行成本后的净值估算 [ C_{curtail} \sum_{t} (P_{curtail,t} \cdot c_{new}) ] 弃电本身不是“花钱”是“少赚钱”但在系统总成本视角下它是调峰能力不足的真实代价必须计入。2.2 净负荷调节需求量化从曲线到系统参数量化成本之前先要把“调峰需求”翻译成可计算的系统参数。我做这一步时用净负荷曲线 (P_{net,t} P_{load,t} - P_{wind,t} - P_{pv,t}) 作为基础输入然后计算几个关键指标净负荷峰谷差(\max(P_{net}) - \min(P_{net}))衡量系统需要应对的最大调节幅度。最大爬坡需求(\max(|P_{net,t} - P_{net,t-1}|))衡量系统需要的最快爬坡速率。低出力持续时间净负荷低于某阈值的连续时段数决定火电需要压出力的深度和时长。这些指标直接决定了机组组合优化的边界条件。例如净负荷峰谷差超过系统可调容量时必须允许机组启停而不是只调出力最大爬坡需求超过单台机组爬坡速率时必须预开多台机组分摊爬坡。2.3 量化路径选择先复现实际调度再做全局最小成本优化量化调峰成本有两条实现路径我的代码里提供了两种模式但默认推荐第二种。**路径A模拟实际调度过程。**按调度员实际采取的操作给定机组组合、给定运行方式做时序仿真统计各项成本。优点是不用求解优化问题速度快、贴近现实缺点是“实际调度”未必经济最优算出的成本可能偏高而且机组组合一变结果就得重算。**路径B机组组合优化。**以系统总运行成本最小为目标求解一个混合整数规划问题得到最优的机组启停计划和出力计划再从中拆分调峰相关成本。优点是结果具备“最小成本”标杆意义可解释性强缺点是建模复杂求解耗时。工程上我建议以路径B为主、路径A做交叉验证。两者结果差距如果很大说明实际调度偏离经济最优太远调度优化本身就有改进空间如果差距很小说明当前运行方式已经比较接近最优量化结果可以直接用来做成本分摊。3. 分摊模型设计责任追溯法与Shapley值法的取舍成本算出来以后最难的不是技术而是“谁来承担”。分摊模型设计得好不好直接决定这个工具能不能在实际项目里落地。我试过两种主流思路各有优劣下面展开讲。3.1 责任追溯法按净负荷边际贡献分摊简单但容易引来争议责任追溯法的思路是净负荷曲线的每一次变化都是由负荷变化、风电出力变化、光伏出力变化共同叠加出来的。那么某类电源对调峰成本的责任就是它出力变化对净负荷波动的“贡献份额”。具体做法先按时间点把净负荷变化量分解(\Delta P_{net,t} \Delta P_{load,t} - \Delta P_{wind,t} - \Delta P_{pv,t})。对每一个时段把系统总调峰成本按三类主体的调峰需求量用各自的出力变化绝对值度量占比来分摊。这个方法的优点是直观、好算、可解释性强Excel都能复算。缺点是它把一个系统性的问题所有主体共同造成的系统波动当作线性叠加来分摊忽视了机组组合的耦合效应。实际做项目时遇到过一个经典争议光伏在午间大发客观上造成了火电深度压出力追责时光伏要承担很大一块成本——但光伏本身是清洁能源被要求承担高额调峰费用在商务上很难被接受。3.2 合作博弈视角Shapley值分摊的公平性与计算代价责任追溯法的争议促使我认真研究合作博弈分摊。Shapley值的核心思想是把风电、光伏、负荷、甚至火电都看作合作博弈的“参与人”总调峰成本是各参与人共同作用的结果。每个参与人应分摊的成本等于它对所有可能联盟组合的“边际贡献期望值”。公式形式为 [ \phi_i(v) \sum_{S \subseteq N \setminus {i}} \frac{|S|! (n-|S|-1)!}{n!} \left( v(S \cup {i}) - v(S) \right) ] 其中 (v(S)) 是联盟S的系统总调峰成本函数需要为每个联盟重新运行一次优化或仿真才能得到。它满足三个公理有效性所有分摊的和等于总成本、对称性相同贡献者分摊相同、零贡献者不分摊。理论上是最“公平”的分摊方案。但代价也很明显计算复杂度随参与人数量指数增长(2^n) 次联盟成本计算在参与人稍多时根本跑不动。我的做法是第一把同类机组聚合为少量“联盟参与人”比如所有火电作为一个参与人风电算一个光伏算一个负荷算一个这样联盟数大幅减少第二用蒙特卡洛采样近似Shapley值对参与人排列随机抽样用几百上千次抽样逼近精确值。实测下来结果稳定性和精确度对工程定价足够了。3.3 工程推荐的两段式分摊流程单用责任追溯法或Shapley值法都有各自的硬伤。实际项目里我推荐一个折中的两段式分摊流程**第一段系统级分摊。**把总调峰成本在“波动性资源风电光伏负荷波动”与“常规电源火电”之间划分。怎么划看总调峰成本与系统不装新能源时的基准调峰成本之差差额部分归因于新能源接入由风电和光伏承担基准部分由负荷侧承担。**第二段新能源内部分摊。**在第一段确定风电和光伏合起来承担的总金额后再在风电与光伏之间细分。细分依据可以是各自的净负荷贡献率责任追溯法、Shapley值、或者实际弃电率——按项目谈判最容易达成共识的口径来选。这套流程的好处总盘子用“增量成本”思路定话说得清楚细分规则灵活商务谈判空间大。代码实现上只需要在成本计算模块之外另写一个分摊计算模块两种算法都预留接口就行。4. Matlab实现路线数据准备、优化建模、结果输出理论框架讲完之后进入代码实现部分。很多同学拿到这类题目最容易卡住的地方不是算法本身而是“数据怎么组织、约束怎么写成矩阵、结果怎么拆出来”。下面按我的实现习惯一步步说清楚。4.1 数据准备与参数对象设计我习惯用结构体数组管理所有机组参数不用凌散的变量堆满workspace。一个典型的机组参数结构如下% 机组参数10台火电机组 gen struct(); gen.name {G1,G2,G3,G4,G5,G6,G7,G8,G9,G10}; gen.Pmax [800, 800, 600, 600, 400, 400, 300, 300, 200, 200]; % 最大出力 MW gen.Pmin [320, 320, 240, 240, 160, 160, 120, 120, 80, 80]; % 最小技术出力 MW gen.ramp [100, 100, 80, 80, 60, 60, 40, 40, 30, 30]; % 爬坡速率 MW/h gen.SC [10, 10, 8, 8, 5, 5, 3, 3, 2, 2] * 1e4; % 启动成本 元/次 gen.SD [8, 8, 6, 6, 4, 4, 2, 2, 1.5, 1.5] * 1e4; % 停机成本 元/次 gen.coal [0.22, 0.22, 0.24, 0.24, 0.26, 0.26, 0.28, 0.28, 0.30, 0.30]; % 煤耗 t/MWh关于时间尺度强烈建议用15分钟一个时段、一天96个点做基本算例。15分钟能反映爬坡过程又不至于像5分钟粒度那样把优化问题撑得太大。第一版代码不要直接上全年8760小时先跑通一天96点的算例结果合理再扩展。负荷曲线和新能源出力曲线可以用readtable从Excel读入也可以先用Matlab内置函数生成模拟曲线。初学者注意新能源出力曲线一定要做归一化处理用装机容量乘以归一化出力序列避免单位和量纲混乱。我习惯把所有数据统一成“MW”和“元”两个单位前后一致排查错误时能省一半时间。4.2 用intlinprog求解机组组合目标函数、约束与关键代码骨架机组组合问题是一个混合整数线性规划MILPMatlab里对应intlinprog求解器。为什么不用linprog因为机组启停状态是0-1整数变量linprog只支持连续变量。核心变量定义方式(P_{i,t})连续变量机组i在时段t的出力MW。(u_{i,t})0-1变量机组i在时段t是否开机。(z_{i,t})0-1变量机组i在时段t是否发生启动从0变1。(w_{i,t})0-1变量机组i在时段t是否发生停机从1变0。目标函数写成Matlab向量形式时需要把二维变量展成一维决策变量x。我的做法是用reshape把矩阵按列展开然后组装目标系数向量f顺序保持一致nG length(gen.Pmax); T 96; % 时段数 % 决策变量顺序: [P(1..nG,1..T), u(1..nG,1..T), z(1..nG,1..T), w(1..nG,1..T)] nVar nG * T * 4; f zeros(1, nVar); % 填入出力成本系数煤耗 × 煤价 % 依次填 z 的启动成本、 w 的停机成本、以及深度调峰惩罚项约束条件里功率平衡约束是等式约束写成Aeq*x beq机组出力上下限、爬坡约束、启停逻辑约束是不等式约束写成A*x b。最关键的约束之一是启停逻辑约束(z_{i,t} - w_{i,t} u_{i,t} - u_{i,t-1})这个约束把“状态切换”和“启停动作”联系起来没有它模型里会出现“启停了但不花钱”的错误结果。爬坡约束的写法也要当心。正确的数学逻辑是 [ P_{i,t} - P_{i,t-1} \leq R_i \cdot u_{i,t-1} P_{i}^{max} \cdot (1 - u_{i,t-1}) ] 这个约束表达的是机组从开机状态爬坡时受爬坡速率限制如果上一时段本来停机那本时段开机可以满出力启动忽略热启动初期出力限制时。初学者最容易踩的坑是写成 (P_{i,t} - P_{i,t-1} \leq R_i)这样会让停机后启动的机组也受爬坡限制导致模型无解或结果保守。调用求解器的骨架options optimoptions(intlinprog, Display, iter, ... RelativeGapTolerance, 1e-3, MaxTime, 600); [x, fval, exitflag] intlinprog(f, intcon, A, b, Aeq, beq, lb, ub, options);其中intcon指定哪些变量是整数z、w和u变量对应编号全部列入lb/ub设置变量上下限连续变量出力范围按gen.Pmin和gen.Pmax设置0-1变量范围设0到1。4.3 从优化结果中剥离“调峰相关成本”求解完成以后fval是系统总运行成本但调峰相关成本需要从结果里分离出来。我的做法是写一个独立函数splitCost(x, gen, T)按时间维度和出力水平做判断。剥离逻辑分三层第一层把启动成本、停机成本直接从目标函数系数里提取——每个启动和停机动作都是一个固定金额累加即得。第二层深度调峰成本遍历每个机组每个时段判断出力是否低于经济出力下限Peco低于部分乘以惩罚系数。第三层爬坡损耗成本遍历相邻时段用出力变化绝对值乘以单位爬坡损耗系数。这里有个容易被忽略的细节如果启停成本和爬坡成本都已经包含在总燃料成本里因为目标函数整体最小化那拆分时必须先去重避免重复计算。我的做法是优化时目标函数里就区分好“常规运行成本”和“调峰附加成本”两组系数求解完成后分别汇总。也就是说目标函数写成 (f_{total} f_{base} f_{peak})其中f_peak部分专门对应上面四类调峰成本项这样拆出来的数天然是对得上的。4.4 Shapley值分摊模块与可视化输出Shapley值计算模块核心是一个联盟价值函数v(S)。对每一组参与人联盟需要重新设定“哪些参与者存在、其出力归零或保留”再调用优化求解器算一次联盟成本。参与人合并以后联盟数控制在8~16之间暴力枚举或者蒙特卡洛都能接受。解析解的代码思路蒙特卡洛采样版更简单function phi shapleyBySampling(v_func, nPlayer, nSample) phi zeros(1, nPlayer); for s 1:nSample perm randperm(nPlayer); for i 1:nPlayer idx perm(1:i); idx_prev perm(1:i-1); phi(perm(i)) phi(perm(i)) v_func(idx) - v_func(idx_prev); end end phi phi / nSample; end可视化输出部分我标配三张图第一张是净负荷曲线各机组出力堆叠图展示机组组合的时序演化第二张是四类调峰成本的堆叠柱状图按时段展示成本分布第三张是各类主体分摊金额的饼图或条形图。作图时注意统一色系我用parula色图出图比较专业。5. 仿真实验与调参避坑记录代码写完之后真正决定成果质量的是实验设计和调试经验。这一节把我实际跑算例的过程、结果规律和踩过的坑完整记录下来方便你直接参考。5.1 三组典型场景的实验结果规律我用同一套负荷曲线只改变新能源渗透率风光装机占比分别为10%、30%、50%做了三组96时段算例。结果规律非常清晰**新能源渗透率越高总调峰成本增长是非线性的。**10%渗透率时系统总调峰成本约占总运行成本的3%~5%50%渗透率时这个比例跳到12%~18%。原因在于一旦净负荷深谷跌破火电最小技术出力总和系统就从“多压一点出力”变成“必须整台停机”启停成本台阶式上升。**爬坡损耗成本在傍晚时段集中爆发。**下午16点到20点净负荷陡增所有火电同时爬坡这一段的爬坡损耗成本能占到全天爬坡成本的70%以上。这个规律提醒你如果做季度或年度分析重点盯住傍晚这个“魔鬼时段”。**Shapley值分摊和线性追责法的结果差异明显。**在50%渗透率场景里线性追责法会分给光伏约40%的调峰成本因为午间深谷主要是光伏造成的而Shapley值只分给光伏约26%风电扛了更多。原因并不复杂——光伏出力集中在午间午间恰好是系统调节压力最大的时段但它和负荷高峰是错开的单看“谁改变了净负荷”会低估这种错峰效益。5.2 intlinprog求解中常见的坑与解决方案这部分的坑我都是一步步踩出来的每个都值得单独提醒。第一个坑模型无解却不知道约束哪里冲突。intlinprog返回exitflag为负值时先别急着改参数用linprog先跑一个放松整数变量的版本看连续松弛是否可行。如果松弛后仍然不可行说明约束矩阵本身有问题如果松弛可行但整数解不可行多数是启停逻辑约束写错了。**第二个坑爬坡约束导致启动机组无法满出力。**前面讲过爬坡约束必须区分“连续运行”和“启动”两种状态。我见过太多人把爬坡约束写成纯差分形式结果模型永远无解——因为一台停着的机组你要让它启动如果不加(u_{i,t-1}0)的豁免逻辑它每时段只能爬(R_i)MW根本开不起来。**第三个坑求解时间太长。**96时段、10台机组、4类变量耦合intlinprog硬跑可能要十几分钟甚至更久。我的经验是按顺序做三件事一是把机组聚合相似的机组合并成等效机组二是时间尺度先粗后细——先用24时段跑通确认结果合理再切到96时段三是设置RelativeGapTolerance到1e-3对于工程分析这个精度完全够用能显著缩短求解时间。**第四个坑结果里出现“启停抖动”。**相邻时段机组反复启停这在物理上不现实但模型为了微小的成本优势可能这么做。解决办法是给启停成本加一个最小运行/最小停机时间约束或者干脆给启停次数设置上限。工程上建议加入最小停机时间约束因为新能源波动下热态机组频繁启停对寿命的影响不可忽视。5.3 判断结果合理性的几个经验指标模型算完以后不要急着把数字写进报告。先对照几个经验指标做合理性校验能挡掉90%的低级错误**指标一单位电量调峰成本。**把总调峰成本除以新能源发电量得到的单位值一般在几元/MWh到几十元/MWh范围。如果你算出几百元/MWh基本可以断定有成本重复计算或者机组参数设置不合理。**指标二启停成本与机组利用小时数的匹配度。**如果模型优化结果里某台机组一天启停3次以上而其启动成本又特别高那大概率是约束条件缺失比如没加最小运行时间需要回查建模逻辑。**指标三深度调峰时段与净负荷低谷时段的对应关系。**深度调峰必然发生在净负荷低位时段如果结果里出现净负荷很高时仍有大量机组深调出力说明经济出力下限或目标函数权重设置有问题。这三个指标都不需要额外数据直接从结果矩阵里统计就能算出来。养成每次跑完先做校验的习惯比反复调参数有效得多。做这套模型最深的体会是调峰成本分摊这个问题的难点大部分不在数学和Matlab代码里而在你对系统物理过程的理解和对“公平”的判断上。我劝所有要拿这个模型做项目或写论文的朋友——第一先花三天把净负荷曲线的“形状”看透比直接写代码重要得多第二分摊规则的选型要预留可解释性选的方案要能对不懂优化的人讲明白第三结果出来以后一定要拿实际调度数据做一次交叉验证哪怕只是粗略对比总量都能帮你发现很多隐藏的建模偏差。后期想扩展还可以在这个框架上叠加储能调峰的动态过程、需求响应的可中断时序以及多区域联络线功率约束调峰成本的种类和分摊维度都可以更细。