高比例可再生能源电力系统调峰成本量化与分摊模型Matlab实现
“高比例可再生能源电力系统”“调峰成本量化”“分摊模型”这三个词放在一起的时候很多做电力系统优化的朋友心里应该已经有画面了风电光伏占比一高净负荷曲线变成“鸭子”火电机组一天到晚启停、深调、爬坡系统运行成本飙升然后大家开会争论——这笔钱到底该算谁的这个题目我拆过很多次也帮人复现过好几版Matlab代码。说实话调峰成本量化与分摊这事核心难点不在数学建模本身而在于“算账口径”和“分钱规则”怎么定。模型框架对不对、求解器选没选对都是执行问题真正费心思的是怎么把“调峰成本”这笔账算得既合理又可落地。这篇文章我就从实际操作的角度把整个建模思路、Matlab实现方案、以及我在调试过程中踩过的坑完整捋一遍。不管你是做毕业论文、工程项目还是科研课题这套方法都能直接拿来改。1. 调峰成本到底要算哪些账1.1 先搞清楚“调峰”的实质在传统电力系统里负荷曲线是相对可控的调峰压力主要来自用户侧的变化——早晚高峰、季节差异。这个时候调峰说白了就是“发电跟着负荷走”。但高比例可再生能源介入以后情况完全变了风电靠天吃饭、光伏日出夜伏电源侧的出力不确定性和反调峰特性让系统的净负荷负荷减去风电光伏出力曲线变得陡峭得多。我习惯用“净负荷波动率”这个指标来判断系统调峰压力它可以直观反映可再生接入前后系统需要调整出力的幅度变化。净负荷峰谷差越大意味着系统需要快速增减出力的范围越大。这就像一辆车传统系统开的是平路只需要偶尔踩刹车高比例可再生系统开的是连续急上急下的山路过弯对发动机扭矩和变速箱换挡速度的要求完全不是一回事。这块多出来的“驾驶压力”就是系统额外付出的调峰成本。量化调峰成本的第一性原理就是把“山路的消耗”和“平路的消耗”做差。也就是说调峰成本并不是某一种燃料成本本身而是可再生能源接入之后系统为了维持功率平衡而额外付出的运行代价。这是整个量化模型的逻辑起点。1.2 调峰成本的五大构成实际操作中我把调峰成本拆成五个部分每一部分都有明确的物理含义和可计算的数学表达式。这五块是成本类别产生机理典型量化方式机组启停成本净负荷峰谷差增大火电机组被迫频繁启停启动耗量折算费用 停机维护费用深度调峰煤耗增量机组低于最低稳燃负荷运行煤耗率显著上升分段煤耗曲线低负荷区间二次项拟合机组寿命损耗成本深调状态下转子热应力、低周疲劳加剧按转子低周疲劳曲线折算停机当量运行备用成本为应对风电光伏预测误差需要预留额外旋转备用备用容量单价 × 备用容量需求增量弃电惩罚/机会成本调峰能力不足时的弃风弃光弃电量 × 上网电价或环境价值系数注意这里的每一项都不是拍脑袋定的。比如深度调煤耗增量我看过不少论文直接用常数煤耗率这个简化在小规模算例里问题不大但如果你的模型要跟实际运行数据对标就必须把煤耗率曲线做成分段函数额定负荷区间用线性段50%以下负荷进入深调区间煤耗率按二次曲线上升。机组寿命损耗这块最容易被人忽略。很多初学者算调峰成本只算煤耗和启停结果总成本偏低后续分摊的时候各主体分到的小份额看起来很美实际运行中根本无法解释。火电机组参与深度调峰转子承受的交变热应力会显著加速低周疲劳损耗这个成本在真实市场里是被计入机组补偿的。在学术模型里我一般用“折算寿命损耗小时”来量化每参与一次深度调峰等效减少若干小时的设备寿命然后按机组单位容量投资成本折算成金额。2. 调峰成本量化的数学模型2.1 目标函数怎么定标准的调峰成本最小化模型本质上是一个带机组组合的优化调度问题。目标函数我是这样写的[ \min \sum_{t1}^{T} \sum_{i1}^{N} \left[ C_i^{fuel}(P_{i,t}) C_i^{su} \cdot y_{i,t} C_i^{sd} \cdot z_{i,t} C_i^{deep}(P_{i,t}) C_i^{loss}(P_{i,t}) \right]\sum_{t1}^{T} \left[ C^{res} \cdot R_t^{req} C^{curt} \cdot \left( P_{t}^{wind,curt} P_{t}^{pv,curt} \right) \right] ]解释一下每一项的含义第一项是常规机组燃料成本( P_{i,t} ) 是机组 i 在时段 t 的出力燃料成本用二次函数拟合。第二项和第三项是启动成本和停机成本( y_{i,t} ) 和 ( z_{i,t} ) 是启停动作的0-1变量。第四项是深度调峰附加煤耗第五项是寿命损耗折算成本。后面就是旋转备用购买成本和弃风弃光惩罚成本。为什么把深调成本、寿命损耗成本和启停成本单列而不并入燃料成本从数学上当然可以合并但单列对后续的分摊有巨大好处——因为不同可再生能源场站对不同成本项的“贡献”是不一样的风电的波动性主要增加启停成本光伏的“鸭子曲线”主要拉高峰谷差进而增加深调成本。如果你把所有成本糊成一团分摊时就完全没有抓手了。2.2 约束条件不可少模型要能真实反映运行物理过程约束条件必须完整。我每次复现这种模型下面几条约束是跑不掉的机组出力上下限约束[ P_i^{min} \cdot u_{i,t} \le P_{i,t} \le P_i^{max} \cdot u_{i,t} ]( u_{i,t} ) 是机组运行状态0是停机、1是运行。这里要注意如果允许深度调峰那么 ( P_i^{min} ) 其实不是固定值我通常把运行状态细分——正常状态、深度调峰状态两种状态的出力下限不一样。深度调峰状态下机组出力可以下探到额定值的40%甚至30%但煤耗率也更高。爬坡约束[ -P_i^{ramp} \le P_{i,t} - P_{i,t-1} \le P_i^{ramp} ]高比例可再生能源场景下爬坡约束经常成为瓶颈。很多时候一个算例无解问题就出在爬坡率设置过紧。我的习惯是先算一下净负荷的最大爬坡需求反推系统需要的爬坡能力再回头校核模型参数。最小启停时间约束[ \sum_{kt}^{tT_i^{on}-1} u_{i,k} \ge T_i^{on} \cdot y_{i,t} ]这个约束是Big-M处理过的形式。最小运行时间和最小停机时间必须建模否则求解器为了省启动成本会给出“频繁启停”这种现实中不可能的操作方案。功率平衡约束和备用约束[ \sum_{i1}^{N} P_{i,t} P_{t}^{wind} P_{t}^{pv} - P_{t}^{curt} L_t ][ \sum_{i1}^{N} \min(P_i^{max} - P_{i,t}, P_i^{ramp}) \ge R_t^{req} ]备用约束里的 ( R_t^{req} ) 可以设成固定比例如最大负荷的5%也可以考虑风电光伏预测误差动态生成。动态设置更准确但需要额外的预测误差数据。我在初版模型里通常先用固定比例跑通后再升级。2.3 关键参数的经验设定参数设得好不好直接决定模型求解速度和结果质量。分享几个我的设定经验机组启动成本燃煤机组一次冷启动成本大约是其满负荷发电2-4小时的燃料费用热启动减半。不要用对称的启停成本现实中停机成本远低于启动成本。深度调峰分段点一般取额定出力的50%以下为深调区。比如一台300MW机组150MW以下算深调。这个设定要参考实际机组的最低稳燃负荷别拍脑袋定。寿命损耗折算我常用“等效运行小时”法每有一次深度调峰低于50%额定出力等效损耗0.1%-0.5%的转子寿命。再根据机组单位造价折算成钱。这个比例各家文献取值差异很大你只要保证算例里量级合理、趋势正确即可。弃电惩罚成本建议设置在可再生能源上网电价的0.5-1.5倍之间既能让模型优先消纳新能源又不至于让惩罚项完全主导目标函数。3. 成本分摊钱怎么分才服众3.1 分摊的三个基本原则成本算完只是第一步真正麻烦的是分摊。调峰成本总额摆在那里怎么在不同电源、不同主体之间分配直接关系到利益格局。我在做分摊方案时始终守三条原则谁引发、谁承担。分摊对象是对调峰需求增量有贡献的主体最典型的就是引入可再生能源后产生的净负荷峰谷差扩大。谁受益、谁分担。调峰成本本质上是为了整个系统的安全稳定运行而付出的代价所有受益者都有责任承担。这个原则在引入储能时特别关键因为储能既是调峰服务的提供者本身也可能因为充放电策略而增加其他机组压力。激励相容。分摊方案要能让各主体“心服口服”。如果一个可再生能源场站发现自己出力越平稳、预测越准分摊的费用越少那这个方案就是成功的。3.2 从责任系数法到Shapley值法实际项目里我用过三种分摊方法从简单到复杂各有适用场景。方法一按调峰需求贡献比例分摊把每个可再生能源场站对净负荷峰谷差的边际贡献作为分摊权重。假设不接入某个场站i时的系统净负荷峰谷差是 ( D_i^{base} )接入后的峰谷差是 ( D^{with} )那么该场站引发的峰谷差增量为[ \Delta D_i D^{with} - D_i^{base} ]然后将总调峰成本按 ( \Delta D_i ) 的比例分摊。这个方法的优点是直观、易解释缺点是只考虑了单一指标的边际变化对负荷时序特征比如爬坡速度不敏感。我的实践经验是这种方法适合做初步估算不适合出正式结果。方法二按运行场景差值分摊思路是在同一个系统里分别求解“含高比例可再生”和“不含高比例可再生”两种场景的运行成本差值就是可再生能源引发的总调峰成本再按各场站的容量占比或电量占比分摊。这种方法实现最简单也是很多论文的基础模型。缺点是完全忽略了场站之间的差异性和交互作用两个出力特性截然不同的风电场分摊到的成本可能完全相同。方法三合作博弈Shapley值法理论上最完善的是Shapley值。把所有可再生能源场站看成一个博弈联盟的参与人联盟的特征函数 ( v(S) ) 定义为联盟S内场站接入系统时的调峰成本变化量。每个场站的分摊成本按所有可能联盟中的边际贡献平均加权[ \varphi_i(v) \sum_{S \subseteq N \setminus {i}} \frac{|S|! (n-|S|-1)!}{n!} \left[ v(S \cup {i}) - v(S) \right] ]Shapley值法的优势是同时考虑了一个场站单独接入和与多个场站共同接入时的边际影响分摊结果满足效率性、对称性等良好性质。代价是计算量巨大——n个场站需要求解 ( 2^n ) 个场景的运行优化模型。n超过8个场景计算时间基本不可接受了。所以在工程落地时我通常采用一个折中方案——随机抽样子集计算近似Shapley值。用蒙特卡洛抽样的方式从全部联盟中抽取一定数量的子集用这些子集的平均边际贡献近似Shapley值。我实测过在5-6个场站的算例中抽取2000次联盟即可让分摊结果的波动控制在3%以内对大多数分析场景足够用了。3.3 分摊结果如何检验分摊结果出来之后别急着写报告。我用两个指标检验结果是否合理个体理性检验任何一个场站单独承担的分摊成本不能高于它与所有其他场站共同承担时的总成本减去其他场站的成本即联盟偏离收益为负。这个条件不满足说明分摊方案存在“反激励”某些场站会倾向于拉帮结派或者退出系统。稳定性检验Shapley值法算出来的结果天然满足核分配条件但责任系数法不满足。检验方法是看能不能找到某个子联盟其成员分摊成本之和大于该子联盟单独接入系统时的调峰成本。如果存在这个分摊方案就是“不稳定”的。这两条检验过关分摊结果才算有说服力。4. Matlab代码实现全流程4.1 整体架构设计我自己写这套模型的Matlab代码从来不用“一个大脚本从头跑到尾”的方式那样调试起来非常痛苦。统一的分层架构是case_study.m % 主脚本数据录入与结果汇总 data_processing.m % 数据处理模块 build_model.m % 构建调度优化模型Yalmip solve_model.m % 调用求解器求解 alloc_shapley.m % 分摊计算模块 plot_result.m % 结果可视化把数据、建模、求解、分摊、出图分离开最大的好处是当你需要改算例规模比如从3台机组改成10台机组只需要动数据文件和主脚本模型构建代码完全不用碰。我早期贪图方便把一堆代码写在一起后面光debug就花了一半时间。4.2 核心代码模块讲解模型构建部分我使用Yalmip工具它可以把线性规划和混合整数规划模型描述得非常简洁。核心代码如下% 定义决策变量 P sdpvar(n_gen, T, full); % 机组出力 u binvar(n_gen, T, full); % 机组启停状态 y binvar(n_gen, T, full); % 启动动作 z binvar(n_gen, T, full); % 停机动作 P_w_curt sdpvar(n_wind, T, full); % 弃风量 P_pv_curt sdpvar(n_pv, T, full); % 弃光量 % 目标函数 objective 0; for i 1:n_gen for t 1:T objective objective fuel_cost(P(i,t), gen.a(i), gen.b(i), gen.c(i)) ... gen.su(i) * y(i,t) gen.sd(i) * z(i,t) ... deep_peak_cost(P(i,t), gen) ... loss_cost(P(i,t), gen); end end objective objective reserve_price * sum(R_res) ... curtail_penalty * (sum(P_w_curt(:)) sum(P_pv_curt(:))); % 约束条件 Constraints []; for t 1:T % 功率平衡 Constraints [Constraints, sum(P(:,t)) ... sum(P_w(:,t) - P_w_curt(:,t)) ... sum(P_pv(:,t) - P_pv_curt(:,t)) L(t)]; % 出力上下限 Constraints [Constraints, gen.Pmin .* u(:,t) P(:,t) gen.Pmax .* u(:,t)]; % 爬坡约束 if t 1 Constraints [Constraints, -gen.ramp P(:,t) - P(:,t-1) gen.ramp]; end % 启动停机逻辑 Constraints [Constraints, u(:,t) - u(:,t-1) y(:,t) - z(:,t)]; end % 求解 options sdpsettings(solver, gurobi, verbose, 1); optimize(Constraints, objective, options);这段代码看起来简单但里面有几个关键细节我要特别提醒第一个坑是燃料成本函数不能直接用二次函数加进去。Yalmip处理 ( P^2 ) 没问题但 ( P^2 ) 的非线性会让求解器变成求解MIQP混合整数二次规划速度比MILP慢很多倍。我的做法是做拉格朗日松弛把二次函数分段线性化。代码里用fuel_cost这个函数接口把二次函数转成线性逼近。第二个坑是启停动作约束的索引对齐。 ( y_t u_t - u_{t-1} ) 只在 ( t \ge 2 ) 时有定义t1时段要用初始状态 ( u_0 ) 处理。很多人漏掉这个边界条件结果第一时段的机组组合结果完全不合理。第三个坑是深度调峰状态的定义。如果你的模型真的有“常规运行”和“深调运行”两档就要把每台机组扩展成两个状态变量分别建模。这个会让模型规模翻倍所以初版模型建议先用单档出力下限跑通后再上深度调峰逻辑。4.3 求解器配置与调试求解器选择上Yalmip默认的sedumi和linprog都很慢我用的是gurobi——在MILP问题上的求解速度有代差级别的优势。以下是配置建议% 求解器参数配置 options sdpsettings(solver, gurobi); options.gurobi.MIPGap 0.01; % 设置1%的允许间隙大幅提速 options.gurobi.TimeLimit 600; % 10分钟求解上限 options.gurobi.MIPFocus 1; % 1寻找可行解优先, 2寻找最优解优先 options.verbose 1; % 1显示求解日志MIPGap这个参数我要重点讲一下。电力系统优化里很多时候不需要证明“严格最优”只需要在最优解1%-2%的邻域内给出可执行方案即可。把MIPGap从默认的1e-4放宽到0.01求解时间常常能缩短一个数量级。我跑过的一个36机组两阶段机组组合模型默认间隙跑了40多分钟没收敛放宽到1%后不到3分钟就出结果了成本差异不到0.5%。还有一点很实用MIPFocus1在模型难找可行解的时候特别好使。它会让求解器优先寻找可行解而不是下边界提升特别适合机组组合这种有严格启停时间约束、容易卡在“无可行解”的模型。5. 案例实测一个典型日的完整算例5.1 算例参数设置为了验证模型的可行性我搭了一个缩微算例。这个算例的灵感来自一个实际风电场周边的简化数据规模控制得比较小方便初学者跑通类型参数数值火电机组G1容量300 MW火电机组G2容量200 MW火电机组G3容量100 MW风电场W1装机150 MW风电场W2装机100 MW光伏电站S1装机120 MW系统峰值负荷—480 MW各机组爬坡率50%额定容量/h—24小时负荷曲线取典型冬季负荷风电和光伏出力曲线用实测数据的时间序列。这个配比下可再生渗透率大约40%净负荷曲线有明显反调峰特征。5.2 求解结果与成本分析跑完模型之后第一件事是看机组组合和出力曲线。结果很典型负荷高峰和光伏大发时段火电机组被压到较低出力甚至停机夜间风电大发时段G3机组基本处于深度调峰状态120MW额定出力降到45MW左右G1机组维持60%以上负荷。具体成本数据如下表所示成本项金额万元燃料成本正常运行段86.4深度调峰附加煤耗12.7机组启动成本8.2机组停机成本3.5寿命损耗折算2.1旋转备用成本4.8合计117.7作为对比把可再生能源出力清零、只保留负荷的“基准场景”跑一遍系统运行成本是74.6万元。也就是说可再生能源接入带来的可量化调峰成本增量约为43.1万元。有意思的地方在于这部分增量中燃料成本只占不到一半启动/停机/深调/寿命损耗这些“隐性成本”合计超过一半。如果你只按燃料成本算调峰成本结果会严重偏低后面分到各方头上的金额也会缺乏说服力。5.3 分摊结果展示接下来用三种方法分别做分摊计算。分摊对象是W1、W2和S1三个可再生能源场站。结果如下分摊方法W1分摊成本万元W2分摊成本万元S1分摊成本万元调峰需求贡献比例法18.215.69.3运行场景差值法16.816.89.5Shapley值法14.717.910.5差异很明显三条曲线各说各话。我的解读是Shapley值法的结果最“讲理”。W1虽然装机容量大但它位于风资源较好且出力平稳的区域调峰压力贡献小所以分摊率反而低于容量占比W2虽然只有100MW但地处风速波动剧烈区域单位电量的出力波动显著分摊了最多成本。光伏在白天削峰、夜间反调峰的特征决定了它比风电更需要系统调峰支持所以在单位容量口径下分摊成本高于W1。这个结果给我们一个重要启发调峰成本分摊不是按装机容量“摊大饼”而是按出力特性“看行为”。同样的装机规模、不同的出力特性对系统调峰资源的消耗可能差出2-3倍。6. 常见问题与排查技巧实录6.1 模型不可行的快速定位这是初学者问得最多的问题也是我最想让你们记住的一段经验。“模型不可行”一句话背后的潜台词是你的约束之间相互矛盾可能没有真实物理意义。我的排查顺序是先检查功率平衡约束。看看负荷数据最高峰是否超过了机组出力上限总和加上可再生能源最大出力。我见过不少算例负荷曲线是480MW但机组总容量只有450MW这当然无解。再检查爬坡约束。当净负荷在某一个小时内从100MW跳到350MW而所有在线机组的爬坡能力之和不到200MW/h必然无解。这时候要么放宽爬坡率要么允许更大弃电。最后检查最小启停时间约束。这个约束最阴险——表面上看每个时段数据都合理但跨时段组合之后某台机组根本无法满足“至少要开机6小时”的要求。我常用的调试方法是逐个打开调节器先注释掉最小启停时间加了功率平衡和爬坡约束之后看有没有可行解再逐步加硬约束找到“压垮骆驼的最后一根稻草”。6.2 求解效率优化模型规模一大求解时间就从秒级跳到分钟级甚至小时级。我的优化手段按优先级排列设置合理的MIPGap前面已经讲过这是性价比最高的提速手段。给变量赋初始值。比如上一轮调度方案的机组组合结果通过assign函数传给优化变量作为备选初始解。Gurobi会从初始解开始搜极大缩短找可行解的时间。减少Big-M系数的量级。比如最小启停时间约束里的常数因子 ( T_i^{on} )不要用T这种全局大数用具体的机组最小运行小时数。Big-M值越紧求解器分支定界的上下界越紧剪枝效率越高。我实测过一个36机组96时段的大算例同一套模型默认设置跑了62分钟优化三步走之后跑进9分钟成本只差0.3%。6.3 结果合理性检验模型能出结果不代表结果一定对。我最后一遍例行检查通常包括看机组组合是否符合“经济顺序”边际成本低的机组负载率是否更高、是否更优先保持在运行状态。如果出现大煤耗机组开机、小煤耗机组停机大概率是模型的启动成本参数设定有问题。看弃电量的时间分布正常情况下弃风弃光应该集中在负荷低谷和风电大发时段。如果随机出现白天弃光就要检查光伏出力和负荷数据是否对齐。看备用约束是否收紧如果备用约束在所有时段都没有起作用可以尝试把备用需求调高20%看系统成本是否显著变化。如果完全没有变化说明备用约束形同虚设需要重新审视备用需求参数的合理性。这三项检查各花不了10分钟但能帮你少被审稿人或者评审专家问倒。6.4 一个容易忽略的成本项最后再分享一个很多人不知道的细节在高比例可再生能源系统中储能系统如果配置了的充放电行为本身也会改变火电机组的运行模式。具体来说储能充满电意味着“等效负荷增加”放电话意味着“等效负荷减少”这相当于在净负荷曲线上叠加了一个“可平移的负荷块”。如果储能充放电策略不够优化反而可能造成火电机组频繁调整出路增加调峰成本。所以在分摊模型中我建议把储能也列为分摊参与主体而不只是把它当作“调峰服务提供者”。我自己加了一版含储能的算例储能场站分摊到的调峰成本为负值——意味着储能确实缓解了系统调峰压力应该获得补偿而不是承担成本。这个结果对市场机制设计很有参考价值。写在最后调峰成本量化与分摊模型说到底是“用清晰的数据回答一个充满争议的分配问题”。模型只是一条技术路径真正有说服力的是你能否把每个参数、每项成本背后的物理过程讲清楚。我在实际使用中发现公式再漂亮不如把系统净负荷曲线和机组组合结果摆出来一对比、把各场站的出力特性数据一展示问题就清楚了一大半。如果你现在正准备复现这套模型我的建议是先别急着把模型做全先把5个机组的微型算例跑通确认各个成本项和分摊逻辑都正确再逐步扩大规模。跑通一个最小可行版本比追求大而全重要得多。后续如果大家在代码实现过程中遇到具体报错也欢迎带着求解日志来交流调试这种模型自己闷头搞容易走进死胡同。