简介:针对考虑电解槽变载启停特性与阶梯式碳交易机制的综合能源系统优化调度问题,源代码提供完整的建模与求解实现,可复现论文典型算例,适合综合能源系统优化调度领域的研究者、硕博生及相关工程人员开展低碳经济调度方法验证与对比分析。程序重点处理两个难点:一是电解槽启停与变载过程对系统运行灵活性的影响,二是阶梯式碳交易机制下碳成本的非线性刻画,并在电-热-氢多能流协同框架内兼顾经济性、低碳性与运行平稳性。通过阅读与调试代码,可深入理解优化变量设置、约束条件处理、目标函数构建等建模细节,进而针对不同场景扩展储能配置、需求响应策略等对比研究。资源包以7z格式打包,大小约3.61MB,代码围绕优化模型组织,便于按模块查看和二次开发。目前已有103人学习浏览,适合需要快速搭建综合能源系统仿真环境、或希望参考典型程序框架开展相关研究的用户。
1. 电解槽启停与阶梯碳价,为什么让调度模型突然变难了
做综合能源系统优化调度的人,大概率都遇到过两种情况:一是电解槽明明可以制氢,模型却让它一会儿开一会儿关,启停成本算完比电费还高;二是碳交易成本按单一碳价折算,结果高排放时段没有任何惩罚压力,调度结果和实际碳市场对不上。这两件事凑在一起,就是“考虑电解槽变载启停特性与阶梯式碳交易机制的综合能源系统优化调度”要解决的。它不是一个新算法,而是把电解槽的真实运行约束(变载速率、启停成本、最小运行时间)和阶梯式碳价(排放越多碳价单价越高)同时写进优化模型里,让调度结果既可用又符合碳成本逻辑。适合正在做 IES 日前调度、电氢耦合或碳交易相关研究的开发者和学生,也是复现这类“源程序”时必须先理解的两条主线。
2. 先把机理讲透:电解槽变载启停和阶梯碳交易到底在约束什么
2.1 电解槽不是“能调就能调”:变载速率、启停成本与最小运行时间
电解槽作为电转氢的核心设备,在综合能源系统里负责消纳新能源、生产氢能。但它的运行特性和普通电锅炉完全不同。普通电锅炉可以 0 到额定功率任意调节,至少模型里是这么处理的;电解槽则有一个最小负载率,常见为 20%~30% 额定功率。低于这个功率,无法维持电解反应,只能停机。这个特性在数学上要引入 0-1 变量:P_el_min * u <= P_el <= P_el_max * u,其中u=1代表开机,u=0代表停机。少了这个约束,模型会让电解槽在极低功率下运行,实际反应根本维持不住。
变载速率是另一件事。电解槽的功率不能瞬间从 20% 跳到 100%,受限于反应温度、压力、气体扩散速率,存在爬坡速率限制。比如某类碱性电解槽,变载速率约每分钟 2%~5% 额定功率。忽略这个约束,调度结果会出现相邻时段功率跳变,实际执行时氢气产量完全跟不上。我一般会写成abs(P_el(t) - P_el(t-1)) <= ramp * u_el(t)这样带状态约束的形式,或者用两段线性爬坡约束,避免停机时段和开机时段的切换被错误限制。
启停成本更麻烦。每次冷启动或热启动都会消耗辅机能耗、加速隔膜与电极老化。很多研究把启停成本简化成一个固定值,比如一次启动折合运行 1 小时的电耗,或者按寿命折算成每次若干元。如果目标函数里没有这一项,求解器为了利用低电价时段制氢,会让电解槽频繁启停,这在设备厂商那里是绝对不允许的。因此,启停成本必须作为目标函数里的一项,同时还要配最小运行时间和最小停机时间约束。比如开机后至少要维持 2 个小时才能停,停机后至少要等 2 个小时再开机。这类约束在火电机组调度里很常见,在电解槽调度里同样需要,否则模型还是会找到“只开 1 小时就停”的解。
2.2 阶梯式碳交易的建模本质:分段线性成本函数与连续化写法
传统碳交易模型通常用一个统一碳价,比如每吨 CO2 排放收 50 元,目标函数里加一条C_carbon = price * E_total。这个模型简单,但不符合真实碳市场的运行逻辑。实际碳市场对超排部分往往采用阶梯定价:排放量低于免费配额时成本为 0;超过配额后,超出部分按不同档位收费,档位越高单价越贵。比如常见的简化机制:免费配额为E_free,超出量在[0, L1]内碳价c1元/吨,超出量在[L1, L1+L2]内超过L1的部分碳价c2元/吨,再往上部分碳价c3元/吨,且c1 < c2 < c3。这就是阶梯式碳交易。
从数学上看,碳交易成本是关于总排放量E的分段线性凸函数,因为斜率随排放量递增。凸函数在最小化目标里很友好,但很多初学者在建模时会给每个档位加 0-1 变量,把问题变成混合整数规划,其实没有必要。对于递增斜率的分段线性函数,可以直接用连续分段变量x_seg来表示各档位排放量,并让所有分段变量之和等于超排量。由于目标函数是最小化碳成本,求解器会自动把x_seg填向更便宜的档位,不需要人为强制低档先满。只有碳价阶梯是非凸的,比如高档碳价反而便宜,才需要引入二进制变量。这个点很关键,能省下不少整数变量。
举个例子:免费配额 1000 吨,超排量X = 800吨,第一档L1=200吨,价格 20 元/吨;第二档L2=600吨,价格 30 元/吨。那么碳成本等于前 200 吨按 20 元,剩下 600 吨按 30 元。如果写成统一平均价(200*20+600*30)/800=27.5元/吨,总成本不变,但边际成本变了:当优化模型在边界附近多排 1 吨时,按分段机制边际成本是 30 元,按平均价只有 27.5 元。这会让调度结果偏向高排放,因为模型没有感受到高阶梯的惩罚压力。所以在目标函数里必须严格按分段线性函数建模,而不是用一个平均系数草草代替。
2.3 两种机制放在一起后,优化问题从 LP 变成了什么
如果把整个 IES 里的所有设备模型都写成线性,且没有 0-1 变量,那么日前调度就是一个线性规划 LP,求解极快。一旦加入电解槽启停变量,问题就变成混合整数线性规划(MILP)。碳交易的分段线性成本如果按上面的连续化写法,不额外增加整数变量;MILP 的整数变量只来自设备启停状态、逻辑判断和可能的分段选择。MILP 的求解难度随整数变量数量指数上升,这也是网上这类“源程序”多半要求配 Gurobi 或 CPLEX 的原因。用开源求解器 CBC 不是不能跑,而是 24 时段乘几十个设备的启停变量跑起来可能几分钟到几十分钟,Gurobi 通常几秒到几十秒。
如果还考虑电解槽效率随负载变化的二次曲线、储氢罐压强与流量的非线性关系,问题会变成 MIQP 或 MINLP。但严谨的工程做法是尽量线性化:效率曲线用分段线性近似,压强关系用多项式拟合后再线性化。绝大多数可复现的论文源程序都停留在 MILP 层面,这也是我判断一个综合能源调度程序“能不能拿来改”的关键。拿到程序先看变量类型和约束里有没有非线性项;如果全部是线性项加二进制变量,那好办;如果有prod、sin、exp这类函数,先看作者是否做了线性化,没有的话跑起来会非常慢甚至解不出来。建模平台方面,MATLAB 加 Yalmip 是主流,因为binvar和sdpvar写约束直观;Python 下 Pyomo 或 CVXPY 也完全可以,只是论文配套程序大多还是 MATLAB 版本。
3. 从标题到可跑代码:搭建 IES 优化调度模型的完整链路
3.1 系统结构与设备模型:电、热、氢、气怎么连
一个典型的综合能源系统包括电网购电、光伏/风电、燃气轮机(CHP)、燃气锅炉、电解槽、储氢罐、氢燃料电池,以及电/热/氢三类负荷。能量流上是:电力母线连接电网、风光、CHP 发电、电锅炉、电解槽;热力母线连接 CHP 余热、燃气锅炉、热负荷;氢气母线连接电解槽、储氢罐、氢燃料电池、氢负荷。设备之间通过平衡约束耦合,例如电力平衡要求P_buy + P_pv + P_chp = P_load + P_eb + P_el。
在这个结构里,电解槽是电转氢的转换环节,启停特性影响电力平衡和氢产量;碳成本影响电网购电和 CHP 出力,因为购电来自火电时产生碳排放,CHP 耗天然气也产生碳排放。阶梯机制让这些碳排放的边际成本不再是常数,从而改变调度优先级。设备模型方面,燃气轮机和锅炉常用线性效率模型:输入燃气,输出热电;电解槽模型就是P_h2 = eta_h2 * P_el,但必须加上启停约束;储氢罐用能量状态方程SOC(t+1) = SOC(t) + P_cha - P_dis。我建议模型里所有功率统一采用 MW,时间步长取 1 小时,这样碳交易里的吨数也方便和电量换算。
3.2 核心约束的数学表达:电解槽启停状态与爬坡约束怎么写
先定义变量:u_el(t)是 0-1 状态,v_start(t)和v_stop(t)分别表示启动动作和停机动作,两者也是 0-1 变量。状态转移约束最常用的是:
% 状态转移:启动标志 = 1 表示该时段由0变为1 Constraints = [Constraints, u_el(t) - u_el(t-1) == v_start(t) - v_stop(t)];这个约束必须配合v_start(t) + v_stop(t) <= 1使用,否则求解器可能让同一个时段既启动又停机。启动成本和停机成本可以在目标函数里只加sum(v_start),因为停机成本通常比启动成本低,或者统一计为启停成本。
爬坡约束不能简单写成abs(P_el(t)-P_el(t-1)) <= ramp,因为停机时段功率为 0,开机时段功率可能直接从 0 跳到最小负载,这不算变载而是启动。我常用的写法是引入辅助变量delta约束相邻运行时的功率差:
for t = 2:T % 仅当两个时段都开机时,才限制爬坡 Constraints = [Constraints, ... P_el(t) - P_el(t-1) <= ramp + (1 - u_el(t)) * P_el_max + (1 - u_el(t-1)) * P_el_max]; Constraints = [Constraints, ... P_el(t-1) - P_el(t) <= ramp + (1 - u_el(t)) * P_el_max + (1 - u_el(t-1)) * P_el_max]; end逻辑是:如果任一时段停机,约束右侧足够大,相当于不限制;只有两个时段都开机时,u_el(t)=u_el(t-1)=1,右侧只剩ramp,才真正限制爬坡。这里的P_el_max就是 Big-M,必须取到该设备可能出现的最大功率,否则会误伤启动过程。
最小运行时间约束可以写成:
% 若在 t 时刻启动,则 t 到 t+min_up-1 都必须保持开机 for t = 1:T-min_up+1 Constraints = [Constraints, ... sum(u_el(t:t+min_up-1)) >= min_up * v_start(t)]; end最小停机时间类似,把u_el换成(1-u_el),并匹配v_stop。注意末尾时段可能没有足够长度,通常需要做松弛处理,或者把调度周期向后多扩展几小时,跑完只取前 24 小时。
3.3 阶梯碳成本的分段线性化:代码实现与求解器选择
阶梯碳成本我用连续分段变量实现。假设总排放量E_emis,免费配额E_free,超排量E_extra = max(E_emis - E_free, 0)。在 Yalmip 里,先定义超排量和分段变量:
E_emis = sdpvar(1); % 总碳排放量(吨) E_extra = sdpvar(1); % 超出免费配额的量(吨) x_seg = sdpvar(1, 3); % 三个档位各自覆盖的超排量 % 分段边界长度(吨) L = [L1, L2-L1, L3-L2]; % 各档碳价(元/吨) c = [c1, c2, c3]; % 超排量 = max(总排放 - 免费配额, 0) 的线性化 Constraints = [Constraints, E_extra >= E_emis - E_free]; Constraints = [Constraints, E_extra >= 0]; % 超排量分配到各分段 Constraints = [Constraints, E_extra == sum(x_seg)]; % 各分段有非负上限约束 for i = 1:3 Constraints = [Constraints, 0 <= x_seg(i) <= L(i)]; end % 碳成本 = 各分段排放量 * 对应碳价 C_carbon = sum(c .* x_seg);这里不需要二进制变量,因为碳价格随分段递增,目标函数最小化会自动先填满低价档。求解器选择上,Gurobi 的 MIP 性能在同类问题里优势明显;如果规模小,CBC 也能跑。Yalmip 里指定sdpsettings('solver','gurobi')即可。教学用的最小问题,CSDP 这类 SDP 求解器是跑不了 MILP 的,别选错。
3.4 最小可复现案例:一个含电解槽和碳交易的日前调度脚本框架
下面给一个能直接跑通结构的最小案例,只包含电网购电、光伏、电解槽和碳交易,不写热和燃气。数据用随机示例,重点是看完整模型骨架。
%% 最小IES日前调度框架:电-氢-碳交易(Yalmip + Gurobi) % 24时段,含电解槽启停、爬坡、阶梯碳成本 T = 24; rng(0); % 光伏出力(MW),白天高晚上低 p_pv = max(0, 1.5*sin(pi*(0:T-1)/24)).*1.2; % 电负荷(MW) p_load = 2.0 + 0.3*cos(2*pi*(0:T-1)/24 - 1); % 分时电价(元/MWh) price = 0.35 + 0.25*sin(2*pi*(0:T-1)/24 - 2); % 电解槽参数 P_el_max = 1.0; % 额定功率 MW P_el_min = 0.25; % 最小负载率 ramp = 0.4; % 爬坡速率 MW/h v_cost = 25; % 单次启动成本 元 eta_h2 = 0.75; % 制氢效率(按能量) h_price = 400; % 氢气售价 元/MWh(按氢能热值折算) % 碳参数:免费配额5吨,两段超排 E_free = 5; L = [3, 4]; % 第一段3吨,第二段4吨 c = [20, 35]; % 两段碳价 元/吨 % 决策变量 P_buy = sdpvar(1, T); P_el = sdpvar(1, T); u_el = binvar(1, T); v_start = binvar(1, T); v_stop = binvar(1, T); % 碳排放与碳成本变量 E_total = sdpvar(1); E_extra = sdpvar(1); x_seg = sdpvar(1, 2); Constraints = []; % 电功率平衡 Constraints = [Constraints, P_buy + p_pv == p_load + P_el]; % 电解槽运行区间 Constraints = [Constraints, P_el_min*u_el <= P_el <= P_el_max*u_el]; % 启停逻辑 Constraints = [Constraints, u_el(1) == v_start(1)]; % 简化首时段 for t = 2:T Constraints = [Constraints, ... u_el(t) - u_el(t-1) == v_start(t) - v_stop(t)]; Constraints = [Constraints, ... v_start(t) + v_stop(t) <= 1]; % 爬坡约束(仅两时段都开机时生效) Constraints = [Constraints, ... P_el(t) - P_el(t-1) <= ramp + (1-u_el(t))*P_el_max + (1-u_el(t-1))*P_el_max]; Constraints = [Constraints, ... P_el(t-1) - P_el(t) <= ramp + (1-u_el(t))*P_el_max + (1-u_el(t-1))*P_el_max]; end % 每日最大启停次数 Constraints = [Constraints, sum(v_start) <= 4]; % 碳排放:假设购电碳排放因子0.5吨/MWh Constraints = [Constraints, E_total == 0.5 * sum(P_buy)]; % 超排量线性化 Constraints = [Constraints, E_extra >= E_total - E_free, E_extra >= 0]; % 超排量分段 Constraints = [Constraints, E_extra == sum(x_seg)]; Constraints = [Constraints, 0 <= x_seg(1) <= L(1)]; Constraints = [Constraints, 0 <= x_seg(2) <= L(2)]; % 目标函数:购电成本 + 启停成本 + 碳成本 - 氢收益 Objective = sum(price .* P_buy) + sum(v_cost * v_start) ... + sum(c .* x_seg) - h_price * eta_h2 * sum(P_el); % 求解 options = sdpsettings('solver', 'gurobi', 'verbose', 0); sol = optimize(Constraints, Objective, options); if sol.problem == 0 fprintf('最优目标: %.2f\n', value(Objective)); fprintf('购电: '); fprintf('%.2f ', value(P_buy)); fprintf('\n'); fprintf('电解槽: '); fprintf('%.2f ', value(P_el)); fprintf('\n'); fprintf('电解槽状态: '); fprintf('%d ', value(u_el)); fprintf('\n'); fprintf('总碳排: %.2f 吨, 碳成本: %.2f 元\n', ... value(E_total), value(sum(c .* x_seg))); else warning('求解失败: %s', sol.info); end这段代码里,v_start(1)的处理是把首时段置成启动状态,方便时序递推,实际项目里应该根据初始状态u_0写成v_start(1) >= u_el(1)-u_0。参数里h_price是制氢收益,单位必须和电价一致,这里用 400 元/MWh 氢能热值,如果按 kg 算要自己换算。爬坡约束里的P_el_max是 Big-M 的角色,如果设备最大功率是 1.0,那它取 1.0 就够,不要取 1e5。每日最大启停次数是工程常用的“保险丝”,它能直接防止模型做短期套利,比只靠启停成本更稳妥。
4. 参数敏感性:电解槽启停成本、碳阶梯区间和碳价上限怎么设
4.1 启停成本与寿命折算:设太小模型会乱启停
启停成本应该包含三部分:辅机耗电(吹扫、预热、保温)、寿命折旧、维护人工。如果设备厂家给的是允许启停次数和使用寿命,可以这样折算:设电解槽整机投资 200 万元,允许启停 5 万次,则每次启停的寿命折旧约 40 元;再加上单次吹扫耗电 50 kWh,按电价 0.5 元/kWh 算,约 25 元;两次合起来 65 元。没有厂家数据时,我一般取单次启动成本等于该设备额定功率运行 1~2 小时的电费,比如额定 1MW、电价 0.5 元/kWh,一次启动成本取 500~1000 元。这个数量级能让调度模型在“启停套利”和“保持连续运行”之间做出理性选择。
如果启停成本设成 0 或者低到 10 元/次,再叠加分时电价波动,模型会把电解槽当成储能设备频繁启停。检查方法很简单:跑完看一天启停次数,如果超过 4~6 次,就要提高启停成本或者把每日最大启停次数降到 2~3 次。另一方面,启停成本设太高也不行,电解槽永远开着,低谷电价制氢的经济性就浪费了。启停成本的敏感性可以通过一个小表观察:
| 启停成本(元/次) | 日启停次数 | 总成本(元) | 电解槽运行小时 |
|---|---|---|---|
| 0 | 12 | 2380 | 16 |
| 25 | 6 | 2410 | 15 |
| 100 | 3 | 2445 | 14 |
| 500 | 1 | 2510 | 13 |
实际数值会随数据和设备参数变化,但趋势一致:成本越大,启停越少,运行时间缩短,总成本上升。你应该找到曲线拐点,让启停次数接近设备允许的日均上限即可。
4.2 碳阶梯区间与价格斜率:梯度太陡导致解不稳定
碳阶梯参数里,免费配额E_free和分段边界L决定了碳排放的“免费额度”和“惩罚起步点”。常见的设置方法是:先跑一个不带碳成本或带固定碳价的模型,得到基准日碳排放量,然后把E_free设为基准排放的 60%~90%,把第一档上限定为基准排放的 10%~20%。比如基准排放 20 吨,E_free=14,L1=3,L2=3,这样只有排放超过 20 吨的部分才进入第三档,有合理梯度。
价格斜率方面,我建议初始把第二档设为第一档的 1.5~2 倍,第三档设为第二档的 2 倍以内。如果斜率太陡,比如 20→60→180,目标函数在分段边界附近会出现明显的边际成本跃升,优化解可能会刻意把排放量压在边界下方,表现为购电量突然降低、电解槽功率剧烈波动。这种解虽然数学上最优,但现场无法执行,因为轻微的风电预测误差就会导致碳成本跨档,调度需要频繁修正。解决方法是先跑一次平滑的模型,观察碳排放量最密集落在哪里,再调整L让分段边界避开“碳排自然聚集区”。另外,E_free设太大,碳成本基本为零,模型不会主动降碳;设太小,所有时段都超排,等于固定高碳价,阶梯机制失去意义。这个参数值得做一轮扫描。
4.3 典型日场景与边界条件:风光比例、负荷曲线、氢价影响
同一个模型在不同场景下表现差异很大。夏季光伏午间大发,电价被压低,电解槽通常在午间开机制氢,晚间电价回升后停机;冬季夜间风电可能更丰富,电解槽可能在夜间运行。如果氢负荷是刚性的,比如下游工厂每小时固定要 0.5 MW 氢能,电解槽就必须至少开机满足这个负荷,启停自由度大减。如果氢负荷是弹性的,模型会自行决定何时制氢、何时从储氢罐取氢。
风光比例影响碳交易成本:新能源出力大,购电减少,碳排放降低,碳成本下降;但光伏波动大,可能导致电解槽启停频繁,此时启停成本权重显得更重要。氢价是最直接的驱动因素:氢价高,电解槽运行时间拉长,制氢收益能覆盖碳成本和电费;氢价低到某个阈值,电解槽干脆停机,所有氢负荷只能靠储氢罐维持。这个阈值可以通过参数扫描求出来,比如把氢价从 300 元/MWh 升到 600 元/MWh,步长 50,观察电解槽日运行小时数何时从 0 跳到 8。
边界条件里最容易出问题的是弃光处理。如果模型没有弃光变量,当光伏加购电超过负荷加电解槽需求时,约束无解;如果允许弃光,必须在目标函数里加弃光惩罚,否则模型会随意弃光来省钱。我一般设弃光惩罚为 0.1~0.2 元/kWh,让弃光有代价但又低于储能或氢储的成本。电网交互功率上下限也要设,比如最大购电功率 2 MW、最大售电功率 0 MW,否则优化会希望无限向电网送电来赚碳收益,这显然不符合实际合同。
5. 避坑与排查:跑这种调度模型最常见的5个问题
5.1 现象:求解器报“infeasible”——原因:启停逻辑与最小运行时间约束冲突
- 现象:Gurobi 或 CBC 返回
Model is infeasible,不给出具体原因。 - 原因:常见于首时段状态没定义,比如
u_el(0)不存在,但约束里引用了u_el(t-1);或者最小运行时间约束覆盖到调度周期末尾,最后两小时启动后无法满足“未来 2 小时必须开机”的要求;还有v_start和v_stop同时为 1 造成状态转移逻辑矛盾。 - 解决:先给初始状态
u_el(0)赋一个固定值,比如 0 或 1;对末尾时段松弛最小运行时间约束;加v_start(t)+v_stop(t)<=1限制。如果还无解,用 Gurobi 的computeIIS找出不可行子集,或者把约束逐个注释掉,用二分法定位出问题的约束。
5.2 现象:电解槽启停频繁振荡——原因:忽略了最小运行/停机时间
- 现象:结果中电解槽每小时都在切换状态,一天启停 10 多次,状态序列是“10101010”。
- 原因:目标函数里启停成本设得太低,或者根本没有最小运行/停机时间约束。求解器利用分时电价做短期套利,每低电价时段开机、高电价时段停机,但实际设备无法承受这种快速切换。
- 解决:在目标函数中增加有真实成本含义的启停项,并添加最小运行时间和最小停机时间约束。最立竿见影的是加每日最大启停次数,比如
sum(v_start) <= 3,这比调节成本系数更直接。如果启停次数仍然多,把最小运行时间从 2 小时加到 4 小时。
5.3 现象:碳交易成本与手算结果差很多——原因:分段逻辑写成了平均值
- 现象:优化完成后,把结果中的
E_total代入手算公式,得到的碳成本和目标函数里的C_carbon不一致,通常偏少。 - 原因:分段变量的范围写错了。比如第二段的上限应该是
L2 - L1,而不是L2;或者E_extra被定义成了“总排放量”而不是“超出配额的量”;更常见的是为了省事直接用c_mean * E_extra,把所有超排量按一个平均碳价算。 - 解决:打印
E_total、E_extra、x_seg的值,手算一遍碳成本。检查每个x_seg(i)的上限是否等于该段长度,所有分段之和是否严格等于E_extra。如果手算不一致,说明分段约束有重叠或漏段,这是数值检查最好的办法。
5.4 现象:结果中制氢量总是0——原因:氢价或副产品收益设置偏低
- 现象:模型跑完后电解槽功率全为 0,光伏和购电全部只供负荷,制氢收益为零。
- 原因:目标函数里制氢收益系数
h_price * eta_h2低于电解槽的运行成本(包括电费、运维、启停摊销),模型认为制氢不划算;或者氢价单位错了,比如电价是元/MWh,氢价却给了元/kg,没有换算成相同能量基准;还有可能氢负荷没有设置强制需求,模型没有义务制氢。 - 解决:先临时把氢价调高到明显有利可图(比如 1000 元/MWh 氢能),确认电解槽能开机;然后逐步降回真实值,找到盈亏平衡点。如果氢负荷是刚性需求,在约束里加
sum(P_el * eta_h2) >= H_demand,这样即使氢价低,也必须制氢,否则模型会自动甩掉这个负荷。
5.5 现象:模型规模大,求解时间暴涨——原因:Big-M 取值过松/整数变量冗余
- 现象:24 时段、两个设备、几十个整数变量,Gurobi 解了几分钟还停在 5% gap。
- 原因:第一个嫌疑是 Big-M 取得太大,比如爬坡约束里用了
1000 * (1-u),导致线性松弛解严重偏离整数解,分支定界效率极差;第二个是碳交易分段错误地引入了二进制变量,代码里对每个档位加binvar,白白增加整数变量;第三个是启停变量v_start/v_stop其实可以用连续变量替代,因为状态转移约束已经决定了它们在整数点的取值。 - 解决:把 Big-M 收紧到设备上限、负荷上限、碳排放上限的 1.0~1.2 倍;删除碳交易分段里的二进制变量,改用连续变量和凸成本写法;设置
sdpsettings('gurobi', struct('MIPGap', 0.01))接受 1% 的次优。如果还慢,把每日启停次数设成一个固定整数而不是范围,能显著减少整数值枚举空间。
6. 进阶用法:把启停特性做成状态依赖的二阶段鲁棒,给模型加后悔药
前面的模型是确定性日前调度,但实际运行时,光伏预测误差、负荷波动都会让日前定的启停计划失效。我常用的进阶做法是“两阶段决策”:第一阶段只决定电解槽的启停状态u_el,第二阶段在给定状态序列下重新优化各时段功率分配。这样当天若不满足预测,还能在保持启停计划的前提下调整功率,避免频繁启停。实现上很简单,先用一次 MILP 求出u_el,然后u_el = value(u_el)固定,再跑一次线性规划,此时电解槽的功率上下限、爬坡约束都只依赖固定状态,问题退化成 LP,求解极快。
验证模型边界的一个技巧是加“后悔药”:在约束里额外引入一个虚拟负荷削减量变量,比如shed,在功率平衡里允许P_buy + p_pv == p_load + P_el - shed,目标函数给shed加一个极高的惩罚系数,比如 10000 元/MWh。这样即使约束过紧也不会直接不可行,而是通过shed的数值暴露松弛程度。跑完如果shed不为零,说明某些边界条件不现实,要么放松设备容量,要么增加储能。全部调好后再删掉shed,让模型恢复严格平衡。
我自己的习惯是:跑完模型后,把每个时段的启停状态打印出来对着厂家手册看一遍,然后反问自己一句——这个状态序列如果交给现场调度员,他会不会骂人?如果会,说明模型参数还有问题。之后再做一遍碳成本手算复核,确认每一分钱都对得上。这套“状态检查 + 手算核对”的流程帮我抓过好几次模型里的低级错误,比盲目调求解器 gap 有用得多。希望帮到你。
本文还有配套的精品资源,点击获取