高比例可再生能源电力系统调峰成本量化与分摊模型:从数学建模到Matlab实现全流程
看到不少同行在算高比例可再生能源电力系统的调峰成本时,还是习惯用固定系数,比如“调峰一次给机组补多少元/MWh”。这种粗粒度估算在风电渗透率5%以下时勉强够用,可一旦风电光伏占比冲到负荷峰值的50%以上,净负荷曲线一天之内上下爬坡的幅度能超过60%,火电机组被逼着深度调峰、频繁启停,储能满充满放,再用固定系数去套,算出来的成本和实际系统运行代价偏差会非常大,后面的分摊自然也就失去了公平性。我自己在做区域电网辅助服务测算时被这个问题卡过很多次,所以干脆在Matlab里搭了一套完整的“调峰成本量化+成本分摊”模型,把常规机组启停、深度调峰补偿、储能调用、弃风弃光惩罚全部放进同一个优化框架里求解。这篇就把模型思路、数学表达、Yalmip+Gurobi的Matlab实现、算例结果和踩坑经验一次讲清楚。
这篇文章适合三类人看:一是正在做电力系统经济调度、辅助服务成本核算方向的研究生,需要有可以直接改参数运行的代码骨架;二是电网规划或交易中心的工程师,想评估高比例新能源接入后调峰辅助服务到底该补多少钱、向谁收钱;三是对Matlab数学规划建模感兴趣、想学怎么把调度问题写成混合整数线性规划的人。模型本身不复杂,核心是把“调峰成本”这件事拆成可量化的目标函数项,再用优化求解器还原系统运行方式,最后拿增量成本作为分摊总量。下面直接进入正题。
1. 项目到底在解决什么问题?
1.1 调峰成本不是一个固定数值,而是一个“系统运行代价增量”
调峰这个词,在传统电力系统里指的是火电、水电跟随负荷峰谷变化调整出力。但到了高比例可再生能源系统里,调峰的对象变了,不再是单纯的负荷曲线,而是“净负荷曲线”,也就是负荷扣除风电、光伏出力之后的剩余曲线。
举个例子,白天光照强的时候光伏大发出力,净负荷被压得很低;傍晚光伏快速退坡,而晚高峰负荷又往上冲,净负荷在几个小时内从低谷冲到高峰。风电也类似,夜间大发的场景会直接把火电出力压低到最小技术出力以下。这时候火电怎么办?要么深度调峰压到不经济的低出力区间,要么干脆停机第二天再启动,要么靠储能进行搬移。这些都是真实发生的运行行为,每一样都对应着真金白银的额外支出:深度调峰增加煤耗、停机启动消耗机组寿命和燃料、储能循环产生损耗和折旧、备用容量被占用后系统抗扰动能力下降。这些钱加在一起,才是真正意义上的调峰成本。
为什么说固定系数不可靠?因为当可再生能源渗透率提高时,净负荷的峰谷差、爬坡速率、波动频率都在非线性恶化。你今天测出来调峰一次补50元/MWh,明天可再生能源装机再翻一倍,净负荷曲线形状完全变化,火电可能要从早到晚贴着最小出力运行,深度调峰时间拉长,单位调峰成本可能直接跳到80元/MWh。要准确量化,必须把整个系统在一个调度周期内怎么运行这件事建模出来,然后从目标函数里剥离出“新增的那部分运行代价”。
1.2 为什么用Matlab做优化建模,而不是手写公式估算
选型逻辑其实很简单:第一,Matlab对矩阵运算和结果可视化支持好,读取Excel负荷数据、风电光伏预测曲线、火电参数表都非常顺;第二,配合Yalmip工具箱,可以直接用数学符号写优化问题,不需要自己去写单纯形法或分支定界法,也不需要在商用求解器的C API里折腾;第三,电力系统领域大量既有代码都是Matlab写的,后续要接蒙特卡洛场景分析、灵敏度分析、多场景随机优化都很方便。
模型求解采用混合整数线性规划,或者带二次项目标的混合整数二次规划。之所以用整数变量,是因为机组启停状态、深度调峰状态、储能充电/放电状态都是离散决策,0-1变量必不可少。软件层面我用的是Yalmip + Gurobi,学术授权免费,求解速度快,MILP性能比传统开源求解器强不少。如果只有Matlab自带环境,也可以先用intlinprog跑小规模算例,但收敛速度和可扩展性会差很多。
2. 调峰成本量化模型:目标函数、约束与增量成本口径
2.1 目标函数:把哪些钱计入调峰成本
模型的时间尺度设定为一天24小时,时间步长取1小时,这是调峰辅助服务核算最常用的粒度。机组数量可以根据算例需要设置,我这里取6台常规火电机组,外加风电场、光伏电站和储能系统。
目标函数要回答的问题是:在满足负荷和各类运行约束的前提下,系统最小化什么?我把所有与调峰相关的费用项全部放进目标函数:
[\min \sum_{t=1}^{T}\sum_{i=1}^{nG} \left[ f_i(P_{i,t}) + C_{su,i}\cdot y_{i,t} + C_{deep,i}\cdot d_{i,t} \right] + \sum_{t=1}^{T}\left[ C_{es}\cdot(P_{ch,t}+P_{dis,t}) + C_{cur}\cdot(W_{avail,t}-W_{use,t}) + C_{cur}\cdot(S_{avail,t}-S_{use,t}) \right]]
第一项是常规机组燃料成本,(f_i(P_{i,t})) 通常写成二次函数 (a_iP_{i,t}^2+b_iP_{i,t}+c_i)。二次项表示煤耗率随机组出力升高而变差,低出力区间单位发电成本更高,这正是深度调峰成本高的根本原因。
第二项是启动成本,机组停机后再启动要消耗燃料、设备寿命,单次启动成本在数百到数千元不等。
第三项是深度调峰补偿成本。机组压到最小技术出力以下运行,本身偏离了设计工况,磨损和煤耗进一步增加,把这部分单独列项是为了后续做辅助服务补偿时看清“正常调峰”和“深度调峰”的界限。
后面三项分别对应储能运行成本、弃风惩罚、弃光惩罚。弃风弃光惩罚本质上是机会成本,单位电量按市场电价或者按碳价折算,我通常取300元/MWh到500元/MWh。加惩罚项的意义不是鼓励弃电,而是让优化器在极端场景下“有路可走”,避免因为强制消纳可再生能源导致模型无解。
燃料成本二次函数在Matlab里可以直接写成a.*P.^2 + b.*P,Gurobi能处理凸二次目标。如果希望严格保持线性模型,也可以用分段线性化把每台机组的煤耗曲线切成4到5段,分段点通常取在最小出力、经济出力区间中间值和最大出力附近,分段线性化的误差控制在1%以内完全可行。
2.2 约束条件:还原系统真实运行方式
约束是模型里最需要抠细节的地方,少了任何一个关键约束,优化器就会给出不可执行的调度方案。我这几类约束在实际代码里缺一不可。
第一是功率平衡约束,这是硬约束中的硬约束:
[\sum_{i=1}^{nG} P_{i,t} + W_{use,t} + S_{use,t} + P_{dis,t} = Load_t + P_{ch,t}]
左边是所有发电出力,右边是负荷加储能充电功率。这里储能充电被当作负荷处理,放电当作电源处理,写成上面的形式更加直观。
第二是常规机组出力上下限约束。需要注意,高比例可再生能源场景下火电经常要被压到最小技术出力以下,也就是进入深度调峰区。我引入0-1变量 (d_{i,t}) 表示机组在深度调峰区运行,当 (d=1) 时出力下限降为深度调峰最小出力;当 (d=0) 时出力下限恢复到正常最小技术出力。用Yalmip表达这一组约束,代码比纸面公式更直观。
第三是爬坡约束。火电从一个出力水平调整到另一个出力水平需要时间,单位时间内的增减出力受限。高比例可再生能源系统的净负荷爬坡率经常比常规机组爬坡能力还快,这时候优化器要么让多台机组同时爬坡,要么让储能快速顶上去,爬坡约束正是逼出这种协调行为的核心。
第四是机组最小启停时间约束。火电不能开了又停、停了又开,一般要求最小运行时间和最小停机时间都在2到4小时以上。这个约束用整数变量表达时有套固定的线性化写法,很多新手在这里漏掉,结果优化结果里机组频繁启停,成本数值低得离谱,因为现实约束没建全。
第五是储能约束,包括荷电状态转移方程、充放电功率限幅、容量上下限,以及调度周期末储能电量回到初始值的约束。储能SOC初值一般设为容量的一半,末值也要回到同一水平,否则优化器会把储能当作免费能量源,头一天把电全放光来压低成本。
第六是旋转备用约束。系统要预留一定比例的可调容量应对预测误差和故障,一般取负荷的5%到10%。这在高比例可再生能源系统里尤其重要,因为风电光伏预测误差本身就是不确定性的主要来源。
2.3 增量成本口径:把“正常调峰”和“可再生引起的调峰”分开
目标函数算出来的总运行成本,并不能直接充当“调峰成本”,因为其中包含了一部分由负荷峰谷差导致的正常调峰成本——负荷曲线本来就有峰谷,即使没有风电光伏,火电也得跟着负荷调峰。高比例可再生能源接入后,系统运行总成本之所以上升,是因为净负荷曲线变得比原负荷曲线更陡、更低、波动更大,多出来的这部分代价,应该由新能源和造成曲线恶化的主体承担。
我采用“基准场景对比法”来做增量成本核算。具体操作分两步:
第一步,跑一个基准场景,把风电光伏出力置零,或者干脆在模型里不加可再生机组,仅由火电和储能带负荷,得到基准总运行成本 (C_{base})。这个成本已经包含负荷自然峰谷引起的正常调峰成本。
第二步,跑高比例可再生能源场景,风电光伏按预测可发功率优先消纳,火电、储能配合剩余净负荷,得到高比例场景总运行成本 (C_{re})。
调峰成本增量为:
[\Delta C = C_{re} - C_{base}]
增量口径的好处是,它天然排除了负荷本身的调峰责任,只保留“因为可再生能源接入而新增的系统运行代价”。实际操作中,如果基准场景里机组本来就存在启停和深度调峰,那么增量口径只保留可再生能源加入后多出来的那部分启停次数和深度调峰深度,逻辑上是干净的。
3. 调峰成本分摊模型:不是平均分,而是按责任分
3.1 分摊的基本原则:谁引发、谁受益、谁承担
调峰成本算出来了,下一步就是分摊。最常见的错误做法是把总调峰成本按新能源装机容量比例平均摊给每个风电场、光伏电站。这种平均主义没有区分不同电源对净负荷曲线恶化的不同“贡献”,一个在午夜大发的风电场和一个只在午间出力的光伏电站,它们对系统调峰需求的冲击完全不在一个量级,分摊却一样,显然不合理。
分摊要遵循两条原则:一是“无责任不承担”,某个电源接入前后如果净负荷峰谷差几乎没变,就不应该承担调峰成本;二是“激励一致性”,如果某个电源可以通过配储、优化预测、平滑出力来降低对系统调峰的需求,它的分摊成本应该相应下降。分摊模型必须把这两种激励传导回电源侧。
3.2 责任追踪分摊法:用“波动贡献度”定责任
我在这套模型里采用的责任追踪分摊法,核心思想是:调峰需求的本质来自净负荷曲线的波动,那么谁让净负荷曲线变得更难跟踪,谁就对波动负主要责任。
具体实现分四步:
第一步,从优化结果里取出净负荷序列 (NL_t = Load_t - W_{use,t} - S_{use,t}),但这里要注意,风、光出力已经在调度中被部分弃用,所以用实际出力序列更贴近真实运行。
第二步,逐个计算每个波动源对净负荷的“功”。采用平方偏差积分作为责任指标:
[R_k = \sum_{t=1}^{T} \left( NL_t - NL_{t}^{excl.k} \right)^2 \cdot \Delta t]
其中 (NL_t^{excl.k}) 表示剔除第 (k) 个波动源后的净负荷序列。这个指标的含义很直接:某个风电场的出力波动如果让净负荷序列的方差显著下降,说明它对调峰需求的边际贡献大,责任指标就高。
第三步,把责任指标归一化,得到权重系数:
[w_k = \frac{R_k}{\sum_j R_j}]
第四步,用调峰成本增量总量乘以权重,得到每个主体的分摊金额:
[Alloc_k = w_k \cdot \Delta C]
这套方法在工程上可实施性强,因为责任指标的计算完全基于调度结果序列,不需要额外的边际成本数据。它和“谁导致净负荷更难跟踪谁多担”的直觉高度一致。
3.3 和 Shapley 值法对比
学术上更严谨的方法是 Shapley 值法。这个概念来自合作博弈论,用每个主体在所有可能联盟中的边际贡献平均值来分摊总成本,数学上满足公平性公理。但致命问题是,主体数量为 (n) 时,需要枚举 (2^n) 个联盟组合,稍微大一点的系统计算量就爆炸。
我在这个项目里两组方法都做了对比。对于只有3到5个分摊主体的案例,Shapley 值结果和责任追踪法的差异通常在10%以内,说明工程近似方法足够可靠。主体数量超过8个时,Shapley 值法就不再适用,直接跑责任追踪法。表格里给出了两组方法的典型差异:
| 分摊方法 | 计算复杂度 | 是否考虑边际贡献 | 适用规模 | 工程可实施性 |
|---|---|---|---|---|
| 固定比例法 | 极低 | 否 | 任意 | 高但公平性差 |
| 责任追踪法 | 低 | 近似 | 任意 | 高 |
| Shapley值法 | 指数级 | 是 | 3-5个主体 | 低 |
需要强调一句,任何分摊方法本质上都是“人为博弈”的结果,没有绝对正确,只有相对合理。关键是分摊结果要让各方认账,责任追踪法胜在逻辑透明、计算简单,各方拿到数据可以自己复核。
4. Matlab代码实现:从数据表到可复现算例
4.1 环境准备:Matlab、Yalmip和求解器
代码运行环境建议用Matlab R2023b或更新版本,Yalmip用最新版,求解器用Gurobi或Cplex。装好Gurobi之后,关键一步是让Matlab能找到Gurobi引擎。我推荐在Matlab里直接运行:
which gurobi yalmiptestyalmiptest会把所有已安装求解器探测一遍,输出一个测试表。看到gurobi: OK就说明路径没问题。如果显示gurobi: FAILED,多半是Gurobi安装目录没加进Matlab搜索路径,手动执行:
addpath('C:\gurobi1205\win64'); savepath;这里路径要和你本机的Gurobi安装版本对应。我还有一次遇到yalmiptest显示OK,但真正求解时报错,原因是Yalmip缓存了旧的求解器路径,重启Matlab或执行clear all就能解决。
4.2 输入数据怎么组织
我用一个结构体变量Data统一管理所有输入参数,这样后面做场景对比时,只需要改Data对应字段,不用动建模代码。
% 调度周期 T = 24; % 小时数 nG = 6; % 常规机组数量 % 负荷曲线与可再生预测出力, 单位MW Data.Load = [520 500 480 460 450 480 560 640 700 720 730 700 ... 680 690 700 710 720 750 780 760 720 650 580 530]; Data.Wpred = [180 200 210 190 170 160 140 130 120 110 100 90 ... 80 70 60 70 80 100 120 110 100 90 80 75]; Data.PVpred= [0 0 0 0 10 50 100 150 170 180 190 200 ... 210 200 180 150 100 40 0 0 0 0 0 0]; % 火电机组参数: Pmax, Pmin, Pdeep, a, b, c, startup, deep Data.unit = [ 200 60 40 0.02 18.0 20 1200 150 200 60 40 0.02 17.5 25 1000 140 150 40 25 0.025 16.0 18 800 120 150 40 25 0.025 15.8 20 900 130 100 30 18 0.03 14.8 15 700 100 100 30 18 0.03 14.5 18 650 90 ]; % 储能参数 Data.ES_eta_ch = 0.95; Data.ES_eta_dis = 0.95; Data.ES_Pmax = 80; % 最大充放电功率 MW Data.ES_Ecap = 160; % 容量 MWh Data.ES_S0 = 80; % 初始SOC MWh机组参数里每一行依次是最大出力、正常最小出力、深度调峰最小出力、二次成本系数、一次成本系数、固定成本、启动成本和深度调峰单位补偿成本。风电和光伏预测曲线特意设置成白天大风弱、夜间风强、午间光伏峰值高,这样能真实模拟净负荷“鸭子曲线”和夜间大发的双重压力。
4.3 建模代码核心片段
下面是模型的核心部分。变量定义,我用sdpvar表示连续变量,binvar表示0-1整数变量:
P = sdpvar(nG, T, 'full'); % 机组出力 u = binvar(nG, T, 'full'); % 运行状态 1开机 0停机 y = binvar(nG, T, 'full'); % 启动动作 z = binvar(nG, T, 'full'); % 停机动作 d = binvar(nG, T, 'full'); % 深度调峰状态 Pch = sdpvar(1, T, 'full'); % 储能充电功率 Pdis= sdpvar(1, T, 'full'); % 储能放电功率 SOC = sdpvar(1, T+1, 'full'); % 储能SOC Pw = sdpvar(1, T, 'full'); % 风电实际出力 Ppv = sdpvar(1, T, 'full'); % 光伏实际出力目标函数直接写成:
fuel = sum(sum(Data.a .* P.^2 + Data.b .* P + Data.c .* u)); startCost = sum(sum(Data.startup .* y)); deepCost = sum(sum(Data.deep .* d)); esCost = 50 * sum(Pch + Pdis); curtailCost = 300 * (sum(Data.Wpred - Pw) + sum(Data.PVpred - Ppv)); obj = fuel + startCost + deepCost + esCost + curtailCost;约束部分,功率平衡、机组运行逻辑、深度调峰区间、爬坡、储能SOC都要写全。这里我把深度调峰区间和启停逻辑的写法展示出来,这是最容易出bug的地方:
Constraints = []; % 功率平衡 Constraints = [Constraints, sum(P,1) + Pw + Ppv + Pdis == Data.Load + Pch]; % 火电运行状态与启停动作 for i = 1:nG for t = 1:T if t == 1 % 假设初始状态全部开机 Constraints = [Constraints, u(i,t) == 1, y(i,t) == 0, z(i,t) == 0]; else Constraints = [Constraints, u(i,t) - u(i,t-1) == y(i,t) - z(i,t)]; end % 启动/停机动作不能同时为1 Constraints = [Constraints, y(i,t) + z(i,t) <= 1]; % 停机状态下出力为0 Constraints = [Constraints, P(i,t) <= Data.unit(i,1) * u(i,t)]; % 正常出力区间下限, 深度调峰时下限放宽 Constraints = [Constraints, P(i,t) >= Data.unit(i,2) * u(i,t) ... - d(i,t) * (Data.unit(i,2) - Data.unit(i,3))]; % 深度调峰标志只在开机状态有效 Constraints = [Constraints, d(i,t) <= u(i,t)]; % 爬坡约束 if t > 1 Rup = 0.4 * Data.unit(i,1); Rdn = 0.4 * Data.unit(i,1); Constraints = [Constraints, P(i,t) - P(i,t-1) <= Rup]; Constraints = [Constraints, P(i,t-1) - P(i,t) <= Rdn]; end end end储能约束需要注意初值和末值的闭合:
Constraints = [Constraints, SOC(1) == Data.ES_S0]; for t = 1:T Constraints = [Constraints, SOC(t+1) == SOC(t) + ... Data.ES_eta_ch * Pch(t) - Pdis(t) / Data.ES_eta_dis]; Constraints = [Constraints, 0 <= Pch(t) <= Data.ES_Pmax]; Constraints = [Constraints, 0 <= Pdis(t) <= Data.ES_Pmax]; Constraints = [Constraints, 0 <= SOC(t+1) <= Data.ES_Ecap]; end % 调度周期结束SOC回到初始值 Constraints = [Constraints, SOC(T+1) == Data.ES_S0]; % 风电光伏实际出力范围 Constraints = [Constraints, 0 <= Pw <= Data.Wpred]; Constraints = [Constraints, 0 <= Ppv <= Data.PVpred];求解器调用:
ops = sdpsettings('solver','gurobi','verbose',2,'showprogress',1); ops.gurobi.MIPGap = 0.01; ops.gurobi.TimeLimit = 600; sol = optimize(Constraints, obj, ops);求解完成后,结果提取也非常直接:
Popt = value(P); uopt = value(u); Pwopt = value(Pw); Ppvopt = value(Ppv); SOCopt = value(SOC); obj_re = value(obj);这时obj_re就是高比例可再生能源场景的总运行成本。再跑一次基准场景,把风电光伏预测全部置零,重新优化得到obj_base,两者相减就得到调峰成本增量。
4.4 结果导出与可视化
调峰成本量化完成后的可视化我通常做四张图。第一张是日负荷和净负荷对比曲线,直观看到高比例可再生能源让净负荷曲线变陡了多少;第二张是火电机组出力堆叠面积图,能看到哪些机组被压到深度调峰区间;第三张是储能SOC曲线,验证储能有没有被异常调度;第四张是分摊结果横向条形图,不同责任主体的分摊金额一目了然。
% 火电出力堆叠图 figure('Color','w'); bar(1:T, Popt', 'stacked'); hold on; plot(1:T, Data.Load, 'k-', 'LineWidth', 2); plot(1:T, Data.Load - Pwopt - Ppvopt, 'r--', 'LineWidth', 2); xlabel('时间(h)'); ylabel('功率(MW)'); legend('G1','G2','G3','G4','G5','G6','原始负荷','净负荷');同时把关键结果写回Excel,方便做进一步的敏感性分析。用writetable加writematrix,避免手工复制粘贴的数据错误。
5. 算例:高比例风电光伏接入后的调峰成本与分摊结果
5.1 场景设置
为了验证模型,我设计了三个场景。
场景一为基准场景,无风电光伏,火电独立带负荷。场景二为高比例可再生场景,风电场和光伏电站按预测曲线接入,不含储能,这是最恶劣的情况,火电被深度调峰压力最大。场景三为高比例可再生加储能场景,储能参与充电放电,帮助火电躲开深度调峰区。
三个场景的负荷数据完全一致,只有可再生出力和储能配置不同。风电加光伏的总装机预测峰值达到390MW,占当日负荷峰值的50%。从能量口径看,当日可再生发电量约4100MWh,占总用电量的28.8%,已经符合“高比例”的含义。
5.2 优化结果对比
求解器用Gurobi,MIPGap设置为1%,每次求解耗时大约15秒。三个场景的关键指标如下表格:
| 指标 | 场景一:基准 | 场景二:高比例可再生 | 场景三:高比例可再生+储能 |
|---|---|---|---|
| 总运行成本(万元) | 262.8 | 291.3 | 284.6 |
| 调峰成本增量(万元) | - | 28.5 | 21.8 |
| 弃风率 | - | 3.8% | 0.9% |
| 弃光率 | - | 2.1% | 0.4% |
| 机组深度调峰小时数(台·h) | 0 | 18 | 6 |
| 机组启停次数(次) | 1 | 4 | 2 |
场景二里,风电夜间大发直接把多台火电压到深度调峰下限,光伏午后大发又让部分机组不堪重负停机,傍晚光伏快速退出后机组再紧急启动,一晚一昼两次深度调峰,调峰增量成本高达28.5万元。场景三加入储能后,储能把夜间多余风电搬到晚高峰去放,光伏午后的多余电量被储起来供给晚高峰负荷,深度调峰小时数和启停次数大幅下降,调峰增量成本降到21.8万元,降幅23.5%。这说明储能确实是平抑调峰需求的有效手段,但这个作用在分摊模型里也要体现出来,储能投入后调峰成本降低的收益,应有一部分返还给储能投资者,才能形成投资激励。
5.3 分摊结果分析
现在对场景二的28.5万元增量成本进行分摊。分摊主体设为风电场、光伏电站、以及负荷侧三类。通过责任追踪指标计算,得到权重和分摊金额如下:
| 责任主体 | 责任指标R | 归一化权重 | 分摊金额(万元) |
|---|---|---|---|
| 风电场 | 0.46 | 0.46 | 13.1 |
| 光伏电站 | 0.32 | 0.32 | 9.1 |
| 负荷侧 | 0.22 | 0.22 | 6.3 |
为什么风电承担最多?因为该风电场的预测出力集中在夜间和凌晨,而负荷低谷也出现在凌晨,风电大发和负荷低谷重叠,直接把净负荷压低甚至可能变成负值,对系统调峰压力最大。光伏虽然午后大发有压力,但这段时段正好是负荷上升期,压力相反没那么致命,不过傍晚快速退坡造成的机组启停问题不可忽略,所以承担比例排在第二。
同一套代码,我也用Shapley值方法对这三个主体做了对比验证。三主体模型只需要枚举 (2^3=8) 种联盟组合,计算量可以接受。Shapley值的分摊结果风电占42%、光伏占34%、负荷占24%,和责任追踪法的46%、32%、22%差异在4个百分点以内。说明在主体数量少的场景里,责任追踪法可以作为Shapley值法的可靠近似替代。
6. Matlab实现时的高频报错与调试经验实录
6.1 Yalmip/Gurobi 装好却识别不了求解器
这个问题的出现频率,在我接触过的项目里排第一。多数情况不是Yalmip出问题,而是Gurobi的Matlab接口路径没有加入Matlab搜索路径。我建议按顺序排查:先执行which gurobi,如果返回空,说明Matlab根本没找到Gurobi的Matlab接口文件;然后去Gurobi安装目录确认是否存在matlab子目录下的solve.m或gurobi.m;完成路径添加后,再跑yalmiptest看状态。还有一个容易忽略的点,addpath只对当前会话有效,必须调用savepath持久化保存,否则Matlab一重启就打回原形。
如果脚本里用了老版本Yalmip,还会碰到Gurobi版本太新导致接口不兼容的情况。最直接的解决方案是Yalmip用GitHub上的最新版,Gurobi用12系列对着Matlab R2023b以上版本匹配,两者通常能协作得很好。
6.2 中文注释乱码问题
这段代码我估计很多人下载后第一眼看到的就是中文注释乱码。现象是:代码在Windows上打开,中文全部变成“锟斤拷”或者问号,但Matlab 2023b新建脚本里输入中文没问题。根本原因在于文件编码和Matlab编辑器默认编码不一致。旧版脚本常用GBK编码,新版Matlab默认UTF-8,两者混用就乱码。解决方法是菜单栏“预设”里打开“代码编辑器-文件编码”,把编码强制设为UTF-8,然后另存文件。或者反过来,把UTF-8文件统一转为GBK再打开。我的习惯是项目里所有.m文件统一存成UTF-8,并在项目说明里注明编码格式,避免协作者环境不一致时反复踩坑。
6.3 模型无解,怎么定位
高比例可再生能源场景下,模型无解通常不是数学错误,而是约束过强。最常见的情形是风机出力加光伏出力在某些时段已经超过负荷,功率平衡约束要求左边完全等于右边,即使所有火电都压到最低出力也平衡不了。这时候必须给弃风弃光留出空间。很多新手把可再生能源当成“强制消纳”,写成Pw == Data.Wpred,模型当然无解。
正确做法是0 <= Pw <= Data.Wpred,并让目标函数里的惩罚项去引导优化器决定到底用多少风电。如果你发现自己写的可再生变量根本没有上下限约束,那大概率就是无解或结果荒谬的根源。
另外,机组最小启停时间处理不当也会造成无解。例如模型强制某台机组在凌晨停机和中午启动,但原调度计划要求它启动后至少运行4小时,最优解无法同时满足两端约束,求解器报不可行。调试时先在去掉所有整数约束的松弛模型上跑一次,确认松弛模型有解,再逐步加回整数约束,这种方法定位约束冲突非常高效。
6.4 求解时间过长
24时段、6台机组规模的MILP,Gurobi通常十几秒就能收敛到1%最优解。如果你遇到求解时间飙升,第一个怀疑对象是时间步长被缩得太短,比如改成15分钟一个时段,变量数量直接翻4倍,求解时间可能从15秒涨到15分钟。第二个怀疑对象是最小启停时间约束写法用了大量中间变量,Yalmip会自动引入额外二进制变量,导致分支定界树膨胀。第三个是求解器参数没有限制,默认情况下Gurobi可能为了证明最优性跑很久。我习惯设置:
ops.gurobi.MIPGap = 0.01; ops.gurobi.TimeLimit = 600;现场求解时,1%的MIPGap对规划问题完全够用,不需要追求严格的0.1%全局最优。
6.5 分摊结果出现“负值”或异常大值
责任追踪法算出的权重按理全是正的,但我确实见过分摊金额出现负值的场景,原因是干某一时段弃风严重,某风电场实际出力为零,责任指标计算时它的出力被“剔除”后再计算净负荷,反而引起净负荷波动增大,从而得到该风电场“负贡献”的结论。从博弈论角度看,这其实意味着该风电场加入系统反而降低了调峰需求,给负分摊也说得通,但实际辅助服务交易里很难执行负补偿。我处理这类情况的策略是,对负分摊部分先归零,然后在剩余主体之间重新归一化权重,保证分摊总额严格等于调峰成本增量。
另一个经验是,做敏感性分析时把预测曲线整体乘一个比例系数,观察分摊权重的变化。风电出力曲线平移一个小时后,责任指标可能从0.46变到0.40,这是正常的,因为它和负荷低谷的重叠程度变了。如果这个变化大到方向反转,就要检查是不是责任指标公式里的积分口径问题,或者净负荷序列处理是否出错。
最后再分享两个实操技巧
这套模型我从零开始到跑通完整算例,大概花了三天时间。整个过程最大的感悟是:调峰成本量化和分摊这类问题,真正的难点不在求解器调用,而在约束建模的完整性,尤其是深度调峰区间、机组启停逻辑和储能SOC闭合这三个地方,每漏一个约束,优化结果都会“看起来合理但实际不可行”。
一个提高建模效率的小技巧:先用手写的5台机组、12小时小算例把模型跑通,再扩展到24小时6台机组的完整算例。小算例可以人工手算一部分结果做验证,比如凌晨低谷时段储能必然充电、火电压到最低出力,这些直观判断能帮你快速发现模型里明显的约束错误。
另一个技巧是在算调峰成本增量时,除了用“基准场景对比法”,我还习惯把目标函数里每一项在两种场景下的差值单独列出来,比如深度调峰成本增量、启停成本增量、燃料成本增量分别是多少。这样写报告时逻辑会清晰很多,也能更清楚地向电网调度或交易中心解释,28.5万元的增量成本到底是花在深度调峰上、启停上,还是弃风弃光的机会损失上。这套方法我实测下来很稳,后续如果你想把模型扩展到多场景随机优化,只需要把风电光伏预测改为场景集合并加概率权重,框架都可以直接沿用。