简介:这份资源面向从事能源系统优化的研究人员与具备MATLAB基础的技术人员,聚焦含阶梯型碳交易机制与电制氢模块的综合能源系统优化调度问题。系统架构涵盖传统机组、风电、电解槽、储氢罐与燃料电池五个单元,重点解决阶梯碳成本分段线性化、氢能系统状态方程构建以及经济性与环保性兼顾的目标函数设计等难点。资源包为1个PDF文件,约129KB,内容以代码讲解与模型推导为主,便于对照阅读与复现。已有161人学习。读者可从中获得完整的YALMIP建模思路、big-M法处理分段碳价、储氢罐时间耦合约束的索引对齐技巧,以及风电大发时段优先制氢、碳排放接近阈值时火电出力调整等案例分析结论,并附有可视化结果与文献参考,适合作为综合能源调度课题的入门与进阶参考。
1. 阶梯型碳交易叠加电制氢:这套综合能源调度到底在算什么
园区级综合能源系统这两年有个明显变化:碳排不再是“算总账”,而是分段计价,超排越多单价越贵;同时电制氢设备从“可选附件”变成“可调负荷”。把这两件事放进同一个优化调度模型里,就是标题说的方向。它解决的是:给定电、热、氢三类负荷曲线和风光出力预测,怎么安排燃气轮机、电解槽、储氢罐、电池的逐时出力,让日运行成本最低,同时碳排落在阶梯区间内。适合做微网调度、园区能源管理、碳资产核算的工程师,也适合正在找“双碳+氢能”交叉课题落地点的研究生。MATLAB 在这里的角色不是仿真动画,而是把混合整数非线性问题拆成可求解的 MILP 或 MINLP,跑出 24 小时调度表。
2. 阶梯碳交易模型怎么建:从线性碳价到分段线性
2.1 为什么不能再用一个固定碳价乘总排放
固定碳价模型里,碳成本等于“总排放×单价”,目标函数是线性的,求解器几秒出结果。但阶梯碳交易把排放量切成若干区间,每个区间单价不同,总碳成本变成分段线性函数。常见做法是设三个区间:基准排放以内单价低,超出 10% 到 30% 单价翻倍,超出 30% 以上单价再翻。这个函数在数学上不是凸的,直接写进目标会让求解器报“非凸”或收敛到局部最优。我一般会引入非负辅助变量把每段排放量拆开,用二进制变量控制区间是否激活,转成 MILP。代价是变量数增加,但 CPLEX 或 Gurobi 对这类结构处理得很稳。
2.2 分段线性化的变量定义与约束写法
假设总碳排 E_total 由燃煤/燃气机组排放和从外部购电的间接排放组成。设三个区间上限为 E1、E2、E3,对应单价 c1、c2、c3,且 c1 < c2 < c3。引入 e1、e2、e3 表示落在各段的排放量,以及二进制变量 y1、y2、y3 表示该段是否被使用。约束如下:
% 阶梯碳交易分段线性化约束(MATLAB + YALMIP 写法) % e_seg: 1x3 各段实际排放量;y_seg: 1x3 二进制,段是否激活 % E_lim: 1x3 各段容量上限;E_total: 总排放量 constraints = []; constraints = [constraints, sum(e_seg) == E_total]; for k = 1:3 constraints = [constraints, e_seg(k) <= E_lim(k) * y_seg(k)]; constraints = [constraints, e_seg(k) >= 0]; end % 段激活顺序约束:前一段用满才允许后一段激活 constraints = [constraints, y_seg(2) <= y_seg(1), y_seg(3) <= y_seg(2)]; % 碳成本表达式 carbon_cost = c1*e_seg(1) + c2*e_seg(2) + c3*e_seg(3);逻辑说明:sum(e_seg) == E_total保证拆分后总量守恒;e_seg(k) <= E_lim(k)*y_seg(k)是 big-M 的紧凑写法,当 y_seg(k)=0 时该段排放被强制为 0;顺序约束y_seg(2) <= y_seg(1)防止出现“第二段有量但第一段没满”的物理不合理情况。参数说明:E_lim 的取值直接决定阶梯位置,常见做法是以系统历史平均排放为基准,第一段取基准值,第二段取基准的 1.1 到 1.3 倍,第三段取 1.3 倍以上。c1、c2、c3 的比值建议 1:2:4 或 1:3:6,具体看地方碳市场指导价,但不要设成 1:1.1:1.2,那样阶梯效果几乎消失,优化结果和固定碳价没区别。
2.3 碳排核算边界:直接排放和间接排放要分开列
很多初学者翻车在核算边界上:把电解槽耗电对应的间接排放和燃气轮机直接排放混在一起算,结果阶梯区间划分失去意义。正确做法是分别建立排放向量。燃气轮机排放 = 出力 × 排放因子;购电间接排放 = 购电量 × 电网排放因子;电解槽如果用的是系统内风光余电,这部分可以算零碳,但要在约束里明确“电解槽优先消纳风光弃电”。我一般会在目标函数里加一个很小的惩罚项,让电解槽在风光富余时优先启动,而不是单纯为了产氢去购电。
3. 电制氢设备建模:效率曲线、启停约束和储氢罐
3.1 电解槽效率不是常数,用分段线性拟合实测曲线
电解槽在低负载率下效率明显下降,如果模型里写“产氢量 = 输入电功率 × 固定效率”,优化结果会倾向于让电解槽频繁在低负载运行,实际产氢量对不上。常见做法是取厂家提供的效率-负载曲线,用 3 到 5 个断点做分段线性拟合。MATLAB 里可以用interp1做插值,但在优化模型里不能直接调插值函数,要手动写分段线性约束。下面是一个三断点拟合的写法:
% 电解槽分段线性效率拟合 % P_el: 电解槽输入功率;H_prod: 产氢量 % 断点:P_brk = [0.2, 0.5, 1.0]*P_rated % 对应效率:eta_brk = [0.55, 0.68, 0.72] P_brk = [0.2, 0.5, 1.0] * P_rated; eta_brk = [0.55, 0.68, 0.72]; % 引入权重变量 lambda(1:3),sum(lambda)=1 constraints = [constraints, sum(lambda) == 1, lambda >= 0]; constraints = [constraints, P_el == sum(lambda .* P_brk)]; constraints = [constraints, H_prod == sum(lambda .* (P_brk .* eta_brk)) / HHV]; % 注意:lambda 中最多两个非零,SOS2 约束 constraints = [constraints, sos2(lambda)];逻辑说明:sos2(lambda)是 YALMIP 里的特殊有序集约束,保证 lambda 中最多两个相邻元素非零,这样 P_el 和 H_prod 就落在同一条线性段上。参数说明:P_rated 是电解槽额定功率,HHV 是氢气高热值,取 39.4 kWh/kg 左右。断点数量和位置根据实际设备手册调整,但不要少于 3 个,否则低负载区误差太大。如果求解器不支持 SOS2,可以用二进制变量手动实现,代价是变量更多。
3.2 启停约束和最小运行时间不能省
电解槽频繁启停会缩短寿命,模型里必须加最小运行时间和最小停机时间约束。常见做法是引入启停状态二进制变量 u_on(t) 和 u_off(t),以及启动指示变量 y_start(t)。约束写成:u_on(t) - u_on(t-1) = y_start(t) - y_stop(t),并且sum(y_start(t:t+T_on-1)) <= 1之类。这部分和火电机组组合问题里的启停约束同构,可以直接套用。我一般会把最小运行时间设成 2 小时,最小停机时间设成 1 小时,具体看设备厂家建议。如果不加这个约束,优化结果会出现电解槽在 24 小时里启停十几次,调度表根本没法执行。
3.3 储氢罐模型:状态转移和容量约束
储氢罐的状态方程很简单:SOH(t) = SOH(t-1) + H_prod(t) - H_use(t) - H_sell(t),加上容量上下限和首末状态相等约束。容易忽略的是储氢罐的充放氢速率限制,以及压力等级对有效容量的影响。如果模型里只写容量约束不写速率约束,优化结果会出现“一小时内把罐子从空充满”这种不现实的操作。我一般会把充氢速率限制在额定产氢量的 1.2 倍以内,放氢速率限制在 1.5 倍以内,留一点裕度。
4. 综合能源系统优化调度:目标函数与求解器选择
4.1 目标函数里到底放哪些成本项
目标函数一般包含:燃料成本、购电成本、碳交易成本、设备运维成本、弃风弃光惩罚。电制氢带来的收益有两块:卖氢收入,以及电解槽作为可调负荷消纳弃电减少的惩罚。我一般会把卖氢收入写成负成本,但单价不要设得太高,否则优化结果会拼命购电制氢,反而增加间接碳排。常见做法是氢气售价按当地工业氢市场价的 0.8 到 0.9 倍取值,留出运输和储存成本空间。碳交易成本用第 2 章的分段线性表达式直接代入。
4.2 用 YALMIP + Gurobi 求解的完整框架
MATLAB 里做这类混合整数线性规划,我习惯用 YALMIP 建模,底层调 Gurobi 或 CPLEX。下面是一个最小可运行框架的骨架:
% 综合能源系统优化调度骨架 clear; clc; % 读取 24 小时风光预测、电热氢负荷 load('forecast_data.mat'); % 包含 P_wind, P_pv, L_ele, L_heat, L_h2 T = 24; % 定义决策变量 P_gt = sdpvar(1,T); % 燃气轮机出力 P_el = sdpvar(1,T); % 电解槽输入功率 P_bat_ch = sdpvar(1,T); % 电池充电 P_bat_dis = sdpvar(1,T); % 电池放电 P_grid = sdpvar(1,T); % 购电量 H_prod = sdpvar(1,T); % 产氢量 H_use = sdpvar(1,T); % 用氢量 SOH = sdpvar(1,T); % 储氢罐状态 u_on = binvar(1,T); % 电解槽运行状态 % 目标函数:燃料+购电+碳+运维+弃电惩罚 fuel_cost = sum(a_gt*P_gt.^2 + b_gt*P_gt + c_gt); grid_cost = sum(price_grid .* P_grid); carbon_cost = c1*e1 + c2*e2 + c3*e3; om_cost = sum(om_gt*P_gt + om_el*P_el + om_bat*(P_bat_ch+P_bat_dis)); curtail_penalty = sum(pen_curtail .* P_curtail); obj = fuel_cost + grid_cost + carbon_cost + om_cost + curtail_penalty; % 约束 constraints = []; constraints = [constraints, P_gt >= 0, P_gt <= P_gt_max]; constraints = [constraints, P_el >= 0.2*P_el_rated*u_on, P_el <= P_el_rated*u_on]; % ... 功率平衡、储能状态、碳排分段等约束 % 求解 ops = sdpsettings('solver','gurobi','verbose',1,'gurobi.TimeLimit',300); sol = optimize(constraints, obj, ops); if sol.problem == 0 disp('求解成功'); P_gt_val = value(P_gt); else disp('求解失败,检查约束可行性'); end逻辑说明:binvar定义二进制变量,sdpvar定义连续变量;目标函数里P_gt.^2是二次项,如果求解器不支持二次规划,可以分段线性化或直接忽略二次项用线性燃料成本。参数说明:gurobi.TimeLimit设 300 秒是防止大规模问题卡死,实际园区级 24 小时问题一般 30 秒内收敛。verbose设 1 可以看到求解器日志,排查不可行时很有用。如果报 infeasible,优先检查功率平衡约束和储能首末状态约束。
4.3 求解结果怎么验证:三个必查项
跑出结果后不要直接信。第一,查功率平衡:每一时刻电负荷+电解槽+电池充电 = 燃气轮机+风光+电池放电+购电,误差应在 1e-4 以内。第二,查碳排分段:e1、e2、e3 是否满足顺序约束,有没有出现 e2>0 但 e1 没满的情况。第三,查电解槽启停次数:如果 24 小时内启停超过 5 次,说明最小运行时间约束没生效或者设得太短。我一般会把这三项写成后处理脚本,每次求解后自动跑一遍。
5. 避坑与排查:阶梯碳交易+电制氢调度里最容易翻车的五件事
5.1 求解器报非凸,结果随机波动
现象:换一组风光数据,目标函数值跳变很大,或者求解器直接报“non-convex”。原因:碳交易分段线性化时用了max或if逻辑,或者电解槽效率曲线直接用了非线性函数。解决:把所有分段函数转成 MILP 形式,用二进制变量和 big-M 约束替代逻辑判断。检查 YALMIP 里有没有max、min、abs直接出现在约束里。
5.2 电解槽全天低负载运行,产氢量对不上
现象:优化结果显示电解槽 24 小时都在 20% 负载附近运行,总产氢量比预期低。原因:效率曲线没做分段线性化,模型以为低负载和高负载效率一样。解决:按 3.1 节的方法加 SOS2 约束,并且把最小负载率设成 30% 以上,低于这个值直接停机。
5.3 碳成本几乎为零,阶梯形同虚设
现象:优化结果里碳排始终落在第一段,碳成本占比不到 1%。原因:E1 设得太大,或者 c1 太低。解决:E1 按系统历史平均排放的 0.8 倍取值,c1 不要低于当地碳市场均价的 0.5 倍。如果还是没效果,检查间接排放有没有算进去。
5.4 储氢罐一天内充放几十次
现象:SOH 曲线锯齿状剧烈波动。原因:没加充放氢速率约束,或者速率约束太松。解决:把充氢速率限制在额定产氢量的 1.2 倍以内,放氢速率限制在 1.5 倍以内。另外可以在目标函数里加一个很小的充放氢惩罚项,抑制无意义的高频动作。
5.5 购电量在碳价高的时段反而增加
现象:阶梯碳价生效后,某些时段购电比燃气轮机发电更便宜,导致间接碳排上升。原因:购电成本里没考虑电网排放因子对应的碳成本,或者碳成本只算了直接排放。解决:把购电间接排放纳入碳排总量,并且在购电成本里叠加对应的碳成本。如果电网排放因子是分时的,按时段分别计算。
6. 进阶技巧:用滚动优化和场景削减让调度表更抗预测误差
前面跑的是日前确定性优化,风光预测误差一大,调度表就偏了。我一般会加一层滚动优化:每 4 小时用最新预测更新未来 4 小时的调度,只执行第一个小时的结果。MATLAB 里实现滚动优化不难,把第 4 章的框架包成一个函数,输入是预测数据和当前储能状态,输出是调度指令。关键点是每次滚动都要把储氢罐和电池的当前状态作为初始条件传进去,否则状态会跳变。
场景削减用来处理风光不确定性。常见做法是生成 100 到 200 个风光场景,用 k-means 或快速前代法削减到 5 到 10 个代表性场景,然后在目标函数里按场景概率加权。下面是一个简单的场景削减代码片段:
% 场景削减:快速前代法简化版 % scenarios: N x T 矩阵,N 个场景,T 个时段 % n_keep: 保留场景数 function kept = fast_forward_reduction(scenarios, n_keep) [N, T] = size(scenarios); prob = ones(N,1) / N; kept_idx = []; % 第一步:选概率距离最大的场景 dist = pdist2(scenarios, scenarios); [~, idx] = max(sum(dist .* prob', 2)); kept_idx = [kept_idx, idx]; % 迭代削减 while length(kept_idx) < n_keep remaining = setdiff(1:N, kept_idx); % 计算每个剩余场景到已选场景的最小距离 min_dist = min(dist(remaining, kept_idx), [], 2); [~, best] = max(min_dist .* prob(remaining)); kept_idx = [kept_idx, remaining(best)]; end kept = scenarios(kept_idx, :); end逻辑说明:pdist2算场景间欧氏距离,min_dist是每个场景到已选场景集的最小距离,选距离最大的加入保留集。参数说明:n_keep一般取 5 到 10,太少会丢失极端场景,太多计算量下不来。概率prob初始均匀,如果场景来自历史数据可以按频次加权。
验证滚动优化效果的方法:用后 7 天的实际风光数据回测,对比日前调度和滚动调度的日运行成本。我自己的经验是,滚动优化能把预测误差带来的成本偏差降低 15% 到 30%,但前提是预测模型本身不能太差。如果预测误差超过 20%,滚动优化也救不回来,得先换预测算法。
最后说一个习惯:每次改完模型,先跑一个 4 时段的简化算例,确认约束和变量没问题,再上 24 小时全量。直接上全量,一旦 infeasible,排查起来就是黑匣子。希望帮到你。
本文还有配套的精品资源,点击获取