基于Matlab与Shapley值法的电力系统调峰成本量化与分摊建模-平芜编程栈

1. 问题重述与模型框架设计

1.1 高比例可再生能源并网后的调峰压力到底从哪来

做电力系统优化这几年，我感触最深的一件事是：风电和光伏装机容量越大的地区，调度员的眉头皱得越紧。过去电网调峰主要靠的是负荷的自然峰谷差，大家心里都有一个大概的经验曲线，什么时候负荷涨、什么时候负荷跌，摸得很准。但当新能源渗透率上来以后，系统净负荷曲线（负荷减去风光出力）开始变得极不稳定，早晨光伏爬坡、傍晚光伏骤降、夜间风电反调峰，这些都是家常便饭。净负荷的剧烈波动直接逼迫火电机组频繁启停、快速升降出力，甚至进入深度调峰区间运行。

调峰压力的本质，其实是系统灵活性资源不足。灵活性资源包括火电的调峰能力、抽水蓄能、储能、需求响应等，它们共同承担“削峰填谷”的任务。当灵活性资源被迫超常规运行，比如火电在50%额定容量以下继续压低出力、一天内启停两次、快速爬坡跨越很大区间，就会带来额外的设备损耗、煤耗升高、寿命折损等一系列成本。这些成本在传统电力系统里数额很小，往往被忽略，但在高比例可再生能源系统里，它们已经成为一笔不可忽视的账。

更麻烦的是，这笔账目前说不清。我们经常听到“谁受益谁承担”的公平原则，但调峰成本涉及多个市场主体：新能源电站给系统带来波动，负荷侧享受了可靠供电，火电厂付出了额外的灵活性服务。到底每个主体该分摊多少？如果分摊机制不合理，火电企业调峰的积极性会越来越低，最终影响整个系统的安全稳定运行。

我做这个项目的初衷，就是想把调峰成本这件事算清楚、分明白。所谓“算清楚”，是通过建立数学模型，把调峰成本的每一个科目精确量化；所谓“分明白”，是设计一套公平、可操作的分摊机制，让新能源企业、负荷主体、火电企业各担其责。Matlab在这整个流程中扮演了关键角色：建模求解、灵敏度分析、结果可视化一条龙完成。

1.2 整体建模思路：两阶段闭环框架

这个项目在架构上采用“先量化、后分摊”的两阶段闭环框架。第一阶段是调峰成本量化，核心是构建一个考虑机组组合（Unit Commitment）和经济调度（Economic Dispatch）的优化模型，把火电调峰的燃料成本、启停成本、深度调峰损耗成本都算出来；第二阶段是成本分摊，核心是把这个总成本按照某种公平原则分配给各市场主体。

两个阶段的关系是：第一阶段算出来的总成本是第二阶段的输入，第二阶段的输出（各主体的分摊费用）又可以为第一阶段的调度决策提供反馈信号。如果分摊结果不合理，可以在下一轮调度中调整火电的报价或新能源的中标电量，形成一个市场化的闭环。这里面其实涉及一个很关键的设计取舍：量化模型要足够精细，但不能过于复杂导致求解困难；分摊模型要足够公平，但不能过于抽象导致无法落地。

我选用Matlab实现还有个现实原因：这个领域的研究者和工程师群体对Matlab非常熟悉，Yalmip工具箱配合CPLEX或Gurobi求解器，几乎是行业标配。相比Python的Pyomo或Julia的JuMP，Matlab在数据预处理、矩阵运算、结果绘图方面更顺手，尤其对电力系统专业出身的人来说，学习成本最低。

2. 调峰成本的构成与量化建模

2.1 火力发电调峰成本的四类构成

调峰成本不是一笔糊涂账，它由多个科目构成，每一项都需要单独建模。我在这套模型里把调峰成本分解为四类：燃料成本、机组启停成本、深度调峰损耗成本和备用容量缩减机会成本。

燃料成本是量化最清晰的一项。机组出力越高，煤耗越大，燃烧效率也会变化。通常用二次函数拟合煤耗特性曲线，再乘以燃料单价得到成本。这一项在传统经济调度里就有，但高比例可再生能源系统里要特别注意低出力工况下的煤耗上升。机组在深度调峰区间运行时，锅炉燃烧稳定性变差，不完全燃烧损失增大，实际煤耗比额定工况高很多，仅仅用常规二次曲线会低估成本。

机组启停成本包括启动成本和停机成本。启动成本又可以细分为燃料消耗、辅助设备耗电、寿命损耗等，与机组停机时间长短相关：停机时间短，机组处于热态，启动快、成本低；停机时间长，机组冷却到温态或冷态，重启需要更长时间和更多燃料。停机过程看似简单，但涉及吹扫、降温操作，也有一定的成本。在高比例可再生能源系统中，火电机组从“两班制”变成“一天两启停”的情况越来越多，启停成本占比显著上升。

深度调峰损耗成本是最容易被忽略、但在高比例新能源系统中极其重要的一项。当火电机组在低于最小稳定运行出力（通常是50%额定容量以下）继续运行，转子等高温部件会承受额外的交变热应力，形成低周疲劳损耗，直接影响机组寿命。这项成本需要借助材料疲劳理论来建模，通常用Manson-Coffin公式描述寿命损耗与应变幅的关系，再折算成等效成本。

备用容量缩减机会成本的概念稍微绕一点。为了保证新能源出力波动时系统安全，调度中心原本会让火电机组预留一定的上调备用容量。当风电大发时，火电出力被压低，机组实际提供的备用容量相应减少，相当于牺牲了机组在能量市场和备用市场上的获利机会。这部分“少赚的钱”也是调峰的真实成本，在市场化环境下必须量化。

2.2 目标函数与核心约束的数学表达

调峰成本量化模型本质上是一个混合整数线性规划问题（MILP）。目标函数是上述四类成本和的最小化，决策变量包括机组启停状态、出力大小、启停动作、深度调峰状态等。用数学语言描述，目标函数大致是这个形式：

min ∑(燃料成本 + 启动成本 + 停机成本 + 深度调峰损耗成本 + 备用机会成本)

约束条件比较多，但都是电力系统调度的经典约束，包括功率平衡约束、机组出力上下限约束、爬坡约束、最小启停时间约束、备用容量约束和污染物排放约束。功率平衡约束是最基础的等式约束：各机组出力之和与新能源出力之和必须等于负荷需求，这个等式在每个时刻都要成立。爬坡约束限制了机组每个时段出力变化的上限，直接影响系统应对净负荷突变的能力。最小启停时间约束反映机组物理特性，火电机组一旦停机，必须等待足够时间才能重新启动。

在Matlab中求解这个模型，我的做法是先定义所有变量，再逐个添加约束，最后调用Yalmip的optimize函数。关键的一点是，所有常数（比如启动成本热态值、爬坡速率上限）要提前定义为参数，不要直接在约束里写数字，这样之后做灵敏度分析会方便很多。

提示：构建任何优化模型，第一步不是写代码，而是把所有的下标、集合、参数列成一张表。比如火电机组编号集合、时间断面集合、负荷预测序列、新能源出力预测序列、各机组最小出力与技术最小出力等。这一步做得越细，后面写代码越顺。

2.3 深度调峰寿命损耗的简化计算方法

深度调峰损耗是本模型中最难量化的一项，我在这里采用了工程上常用的简化方法。机组进入深度调峰状态后，关键部件承受的热应力变化可以用等效运行小时数来折算。具体做法是，设定一个基准工况（比如75%额定出力），通过寿命损耗计算得到基准工况下每运行一小时对应的损耗值，再按深度调峰出力和基准出力的比值进行修正。

修正系数可以用经验公式拟合。我国电力行业对火电机组深度调峰有比较成熟的试验数据，一般在50%-30%额定出力范围内的损耗修正系数从1.2倍到2.5倍不等。把修正系数乘以该时段内的基准损耗成本和机组电价，就得到这一时段内的深度调峰损耗成本。这个方法的优点是计算简单、参数容易获取，缺点是精度受到经验数据的影响。在实际项目中，如果手头有具体机组的疲劳试验报告，可以直接替换修正系数。

3. 调峰成本分摊模型的设计

3.1 为什么选择Shapley值法作为核心分摊工具

调峰成本算出来之后，第二个核心问题就是：这笔钱怎么分。我在对比了多种分摊方法后，最终选择了合作博弈论中的Shapley值法作为核心分摊工具。

为什么不用简单的“按发电量比例分摊”或者“按峰值负荷分摊”？这两种方法虽然计算简单、容易解释，但存在明显的公平性缺陷。按发电量比例分摊的话，承担调峰责任的火电企业分到的成本比例反而不高，因为它们在新能源大发时段出力很低；受益于调峰的新能源企业却因为自身发电量高而承担了最多成本，这显然不符合“谁受益谁承担”的原则。按峰值负荷分摊则完全忽略了波动性的贡献，只要用户峰值负荷相同，哪怕其中一个用户负荷平稳、另一个用户负荷剧烈波动，分摊的费用也一模一样，这也不公平。

Shapley值法来自合作博弈理论，它的核心思想是：每位参与者应该获得的成本分摊额，等于该参与者对所有可能的合作联盟边际贡献的期望值。放在调峰成本分摊这个情景下，“参与者”可以是新能源电站、负荷聚合商、火电企业的组合。某个参与者给联盟带来的调峰成本增量越大，它分摊的调峰成本就越高。

3.2 Shapley值计算的数学原理与实现步骤

Shapley值法的计算公式看起来有点抽象：对包含N个参与者的集合，每个参与者i的Shapley值定义为对所有不包含i的子集S，计算(S的规模×(N-S的规模-1)!/N!)?(v(S∪{i})-v(S))的累加和。其中v(S)是联盟S的总调峰成本，称为特征函数。

用大白话说，这个公式就是在模拟这样一个过程：把所有参与者按随机顺序排成一队，逐个加入联盟，每加入一个参与者，联盟成本就会增加一个边际量。每个参与者取自己所有可能的“加入顺序”下的边际成本的加权平均，就是它应该承担的调峰成本。这个机制同时保证了效率性（所有参与者分摊额之和等于总成本）和公平性（对称参与者分摊相同成本、不参与调峰的参与者分摊为零）。

在Matlab中实现Shapley值计算，核心是枚举所有可能的参与者子集。当参与者数量小于等于10时，2的N次方减1个子集可以快速遍历；参与者数量超过10后，计算量指数增长，这时候需要用蒙特卡洛模拟近似计算Shapley值，以随机采样联盟顺序的方式估计边际贡献均值。我在实际项目中设置了两种模式，默认用精确计算，在参与者数量大于12时自动切换为蒙特卡洛近似。

项目案例里，我设置了三个典型参与者：风电场A、光伏电站B、和负荷中心C。通过运行调峰成本优化模型，得到“只有风电+负荷”“只有光伏+负荷”“三者都在”等情景下的系统总调峰成本，以此作为特征函数值，最终算出风电场A分摊30%、光伏电站B分摊25%、负荷中心C分摊45%。这个结果很好解释：负荷中心虽然总额不变，但它的日内波动特性是调峰需求的根源，所以承担的比例最高；光伏因为集中在白天出力，与负荷曲线相对匹配，分摊比例反而低于风电。

3.3 不同分摊机制的对比分析

Shapley值法算出来的结果是理论上最公平的，但在实际应用中，它的计算复杂度和解释成本都偏高。行业里还有其他几种常用分摊机制，我在项目中做了对比分析。

按电量比例分摊是操作最简单的方式，只需要知道各主体一年的发电量或用电量，按比例分账即可。这种方式的优点是数据要求低、规则透明，缺点是完全没有体现各主体对调峰的真实贡献或责任，容易导致“劣币驱逐良币”：负荷曲线平稳的用户为波动性用户买单，出力稳定的电厂为出力随机的电厂买单。

按最大净负荷波动贡献分摊是另一种思路，先计算各主体单独接入系统时造成的净负荷峰谷差增量，再按增量比例分摊调峰成本。这种方法能体现波动性的影响，但忽略了各主体之间的相互抵消效应。风电场和光伏电站的波动方向可能相反，单独计算再相加会高估总波动成本。Shapley值法的优势恰好在这里：它考虑所有参与者组合的边际效果，能准确捕捉主体之间的交互作用。

三种方法各有优劣，我对它们的结论是：如果系统规模小、数据质量高、各方对公平性要求较高，优先用Shapley值法；如果系统规模庞大、实时性要求高、解释成本敏感，可以用Shapley值法计算结果作为基准，再拟合出一个简化的分摊公式来逼近。

4. Matlab代码实现过程与算例演示

4.1 代码整体结构与数据准备

整套Matlab程序我按功能划分成五个模块：数据输入模块、参数设置模块、优化模型构建模块、Shapley值计算模块、结果输出与可视化模块。这样分手的好处是可以单独调试，数据换了不用改模型，模型改了不用动数据，后面维护非常省心。

数据输入模块读取的是系统基础数据，包括各火电机组的额定容量、最小出力、爬坡速率、煤耗曲线系数、启停成本参数，以及典型日24个时段的负荷预测数据和新能源出力预测数据。我习惯把这些数据统一保存在Excel文件里，程序开头用readtable函数批量读取，避免在代码里硬编码。参数设置模块集中定义所有常量：燃料单价、备用容量价格、深度调峰损耗修正系数等。

下面这段代码是核心模型构建部分的骨架：

% 定义优化变量 u = binvar(nGen, nTime, 'full'); % 机组启停状态，1表示运行 p = sdpvar(nGen, nTime, 'full'); % 机组出力，单位MW startUp = binvar(nGen, nTime, 'full');% 启动动作 shutDown = binvar(nGen, nTime, 'full');% 停机动作 % 功率平衡约束 Constraints = []; for t = 1:nTime Constraints = [Constraints, sum(p(:,t)) + sum(Pwind(:,t)) == Pload(t)]; end % 机组出力上下限约束 for g = 1:nGen for t = 1:nTime Constraints = [Constraints, Pmin(g)*u(g,t) <= p(g,t) <= Pmax(g)*u(g,t)]; end end % 目标函数：燃料成本 + 启停成本 + 深度调峰损耗成本 Objective = sum(sum(fuelCostCoeff1 .* p.^2 + fuelCostCoeff2 .* p + fuelCostCoeff3 .* u))... + sum(sum(startUpCost .* startUp)) + sum(sum(shutDownCost .* shutDown))... + sum(sum(deepPeakCost .* deepPeakFlag));

4.2 一个24时段的实际算例演示

为了让读者更直观地理解模型工作流程，我设计了一个简化的算例：系统含2台火电机组（一台600 MW、一台300 MW）、一个200 MW风电场和一个50 MW光伏电站，负荷数据采用某地区典型夏季日负荷曲线，风电出力按照“夜间大、白天小”的反调峰特性设置，光伏出力按照“中午大、早晚零”的典型曲线设置。

首先运行调峰成本优化模型，计算出24个时段的机组组合结果。结果很典型：夜间风电大发时段，两台火电机组中的一台被迫深度调峰到40%额定出力附近，另一台直接停机；白天中午光伏大发时段，300 MW机组短暂停机；傍晚光伏骤降时段，停机的机组需要快速热启动，承担晚高峰负荷。整个优化过程求解时间不到5秒（使用Gurobi求解器），总调峰成本为42.6万元。

接着运行Shapley值分摊模块，把风电场、光伏电站、负荷中心三个主体作为参与者，计算各自应分摊的成本。特征函数的计算需要跑7次优化模型（2的3次方减1），每次求解2秒左右，整个过程不到20秒。计算结果是：风电场分摊15.2万元（占比35.7%）、光伏电站分摊9.8万元（占比23.0%）、负荷中心分摊17.6万元（占比41.3%）。负荷中心占比最高，原因是该地区典型日负荷峰谷差本身较大，即便没有新能源接入，也需要火电深度调峰来平衡负荷。

4.3 结果可视化与灵敏度分析

Matlab在结果可视化方面的优势在这个项目里体现得很充分。我绘制了三类图：第一类是系统运行结果图，包含负荷曲线、风电光伏出力曲线、各火电机组出力曲线和启停状态时序图，这种图可以直观展示整个系统一天的调度过程；第二类是调峰成本构成饼图，清晰对比燃料成本、启停成本、深度调峰损耗成本和备用机会成本的占比；第三类是各主体分摊比例的横向柱状图。

灵敏度分析是做这个项目不可跳过的一步。我重点考察了两个关键参数的影响：新能源装机容量的变化和天然气价格（火电燃料成本）的变化。结果显示，随着新能源渗透率从20%提高到40%，总调峰成本几乎呈指数增长，但单位新能源发电量承担的调峰成本却先降后升，原因是低渗透率时系统还有足够闲置调峰容量，高渗透率时系统被迫频繁启停机组，边际调峰成本快速上升。这个结论对政策制定有参考价值：不能简单地说“新能源越多越好”，还要看系统调峰资源的匹配程度。

5. 常见问题排查与实操经验分享

5.1 建模和求解中的典型坑

我在开发这套模型的过程中踩了不少坑，总结几个典型的希望对大家有帮助。

第一个坑是启停变量和状态变量之间的逻辑约束缺失。很多初学者只定义了机组的运行状态变量u，却没有定义启停动作变量，导致优化结果里出现“机组没有任何动作就改变状态”的荒谬结论。正确做法是添加这样的约束：启动动作大于等于当前时刻状态减上一时刻状态，停机动作大于等于上一时刻状态减当前时刻状态。没有这组约束，计算出的启停成本就是错的。

第二个坑是爬坡约束的循环索引错误。爬坡约束描述的是相邻时段出力的最大变化量，在Matlab里用循环写的时候很容易粗心把t和t-1写反，或者第一时段约束没有单独处理。我习惯先把第一个时段的约束单独写出，再从第二时段开始循环，这样可以避免索引错误。

第三个坑是Yalmip参数传递的错误使用。在计算Shapley值时，需要多次修改参数重新求解优化模型，如果用assign函数直接给sdpvar变量赋值，会造成旧的变量定义残留在求解器内存中，影响后续求解，甚至导致结果错误。正确做法是对每个情景重新定义相关的约束和目标函数，或者直接用optimize函数传入新的参数。

第四个坑是求解器选择问题。默认的求解器（如sedumi）处理MILP问题时效率很低，几分钟都算不完。一定要安装并配置CPLEX或Gurobi，并在代码里通过Yalmip的sdpsettings指定求解器。配置的关键是确保求解器路径正确添加到Matlab中，这一步虽然简单，但很多人在这里卡住。

5.2 数据质量与单位一致性检查清单

电力系统模型对数据质量和单位一致性要求极高。我整理了一份自检清单，每次运行模型前都会过一遍：

所有功率数据的单位统一为MW，能量单位统一为MWh，费用单位统一为元或万元；
负荷曲线和新能源出力曲线的时间分辨率一致，如果负荷是15分钟一个点，新能源是1小时一个点，需要先做插值或聚合处理；
火电机组的煤耗曲线系数单位是吨/兆瓦时还是千克/兆瓦时，要换算清楚；
爬坡速率单位是兆瓦/小时，如果是兆瓦/分钟需要乘以60；
启停时间约束的单位是小时，如果时间分辨率是15分钟，minUp和minDown要换算成对应的时间步数。

5.3 提高求解速度的几个实用技巧

算例规模一大，MILP求解就会变得很慢。我在实践中总结出几个提速技巧。

第一个技巧是给变量提供好的初始解。先用简化模型（去掉整数约束的连续松弛问题）求解，得到的连续解作为MILP问题整数变量的初始可行解，配合CPLEX的MIP start功能，可以显著减少分枝定界搜索时间。第二个技巧是添加对称破缺约束。多台同类型机组在模型中是完全对称的，求解器会重复探索大量对称解。通过按机组编号强制排序的约束（比如机组1的出力不低于机组2的出力），可以大幅压缩搜索空间。第三个技巧是合理设置求解器参数，比如开启CPLEX的并行计算、设置合理的MIP gap容忍度（通常设置为0.5%-1%就够了，没必要追求0%最优解），求解速度能提高5到10倍。

5.4 从研究到落地的扩展思路

这个模型目前是基于日前调度的框架，实际系统运行中还有日内滚动调度和实时调度的层级结构，调峰成本在不同时间尺度上的量化方法并不相同。后续可以在模型中加入多时间尺度协调机制，让调峰成本分摊更贴合实际运行方式。

另一个值得探索的方向是引入储能系统。储能和火电联合调峰时，调峰成本的结构会发生变化：火电的启停次数减少，但储能电池的循环寿命损耗成为新的成本项。Shapley值框架同样适用，只是需要把储能电池作为新的参与者加入联盟，重新计算各参与者的边际贡献。

我对实际落地还有一个体会：调峰成本分摊模型算出来的结果再好，最终都要面对“谁来买单”的现实问题。在实际工程中，分摊方案需要与市场规则、结算体系衔接，还要兼顾各方可接受度。这让我想起在一个项目验收会上一位老调度员的评价：“模型是冷的，但分钱是热的问题。”所以在输出结果时，我都会附上完整的分摊依据和可解释性说明，让各方能看到自己的成本构成来源，这样的模型才能真正有生命力。