我自己接触综合能源优化调度这个方向也有几年了从最开始只做电热联供到后来把氢能加进去中间踩了不少坑也实实在在体会到氢能在整个系统里那种“四两拨千斤”的作用。这几年但凡涉及园区级低碳调度、新能源高比例消纳的课题几乎都绕不开氢能这个环节。很多同行朋友问我说看着论文里那些“电-热-氢”联合调度的结果很漂亮但自己拿到Matlab代码一跑要么约束写错要么求解器报错要么结果不合理根本不知道从哪里改起。这篇文章我就结合自己做过的实际项目把计及氢能的综合能源优化调度研究从头到尾捋一遍包括模型怎么建、约束怎么写、代码怎么组织、参数怎么调以及对新手最有用的排错经验一次性说清楚。1. 为什么要做“计及氢能”的综合能源调度先说个最直观的场景一个典型的风光互补微能源网白天光伏大发、晚上风电出力足但用户侧负荷根本吃不下这么多电尤其是夏季午间和冬季夜间弃风弃光率能到15%以上。传统的做法是配电池储能把多余的电存起来但电池的问题很明显——容量有限、循环寿命受充放电次数制约、长期运行还有容量衰减而且电池只能解决“电-电”的时空平移没法把富余电能转成别的能量形式去服务热负荷、气负荷更不能在更长的时间尺度上消纳波动。氢能的价值恰恰就在这里它通过电解水制氢把富余电能转化成化学能储存在氢气里之后再通过燃料电池发电、通过氢锅炉供热、甚至直接供应工业氢负荷。这就相当于给能源系统加了一个“跨介质缓存”既能缓解电力紧张又能参与热力平衡还能产出有市场价值的氢气副产品。我做过的几个实际案例里接入氢能之后的典型收益体现在三个维度弃风弃光率从原来的12%~18%压降到5%以内风光利用率显著提升。在实时电价较高的时段让燃料电池并网放电系统购电成本降低8%~12%左右。通过氢储能设备参与电热联供整个系统的运行碳排放强度下降了约20%。从调度建模的角度看引入氢能之后系统的耦合变量变多了。原来我们只需要处理电功率平衡、热功率平衡两个约束现在还要考虑氢气的生产、存储、消耗平衡以及电解槽、储氢罐、燃料电池、氢锅炉这四类核心设备的运行约束。变量类型也从连续变量扩展到了含0-1状态变量设备启停和整数变量储氢罐分级的混合整数优化问题。这正是很多初学者一上来就被绕晕的地方。2. 系统架构与建模思路先画清楚能量流再写方程2.1 物理架构怎么搭开始建模之前我习惯先把物理架构用能量流图的形式画出来。一套典型的“电-热-氢”综合能源系统通常由下面这些环节组成供能侧光伏PV、风电WT、上级电网购电Grid、天然气购气Gas。转换侧电解槽EL、燃料电池FC、燃气轮机GT、燃气锅炉GB、氢锅炉HB、电锅炉EB。储能侧蓄电池BESS、储氢罐HST、储热罐TST。负荷侧电负荷、热负荷、氢负荷。这个架构里最关键的设计决策是“氢能的去向”。我做过两种模式的对比纯燃料模式氢气只进燃料电池发电燃料电池余热再供热这种模式对电平衡帮助最大。综合供能模式氢气可以进燃料电池发电也可以进氢锅炉直接供热还可以直接卖给外部氢负荷。实际项目里第二种模式更贴近“综合能源系统”的设定因为氢负荷的存在相当于给系统增加了一个可调节的出力口——电价低时多制氢存起来热价高时氢气去供热电价高时氢气去发电灵活性非常强。2.2 目标函数怎么定计及氢能的调度问题目标函数一般有三种形态最小化运行成本、最小化碳排放、或者成本碳排的加权多目标。我平时最常用的是运行成本最小化因为数据好获取而且工程上最关心钱。运行成本包括购电成本分时电价下各时段从上级电网购电的电量乘以对应电价。购气成本燃气轮机和燃气锅炉消耗的天然气费用。设备运维成本电解槽、燃料电池、蓄电池等设备的运行维护费用通常按出力大小线性折算。弃风弃光惩罚成本为了保证模型更倾向于消纳可再生能源可以给弃电设一个惩罚项。氢负荷收益可选如果系统对外供氢可以在目标函数里减掉售氢收入。目标函数的数学形式大致是这样[ \min \quad C \sum_{t1}^{T} \left( c_{e,t} P_{grid,t} c_{g} V_{gas,t} \sum_{i} k_i P_{i,t} c_{curtail} P_{curtail,t} \right) ]其中 (c_{e,t}) 是t时段购电价(c_g) 是天然气单位热值价格(k_i) 是设备单位出力运维成本(c_{curtail}) 是弃电惩罚系数。2.3 约束条件怎么分类我把约束分成四类功率平衡约束、设备运行约束、储能动态约束和系统级耦合约束。下面逐一展开。功率平衡约束是模型中必须严格满足的等式约束电功率平衡光伏出力风电出力电网购电燃气轮机发电燃料电池发电蓄电池放电 电负荷电解槽耗电电锅炉耗电蓄电池充电弃电量。热功率平衡燃气轮机余热燃气锅炉产热氢锅炉产热电锅炉产热储热罐放热 热负荷储热罐吸热。氢气平衡电解槽产氢储氢罐放氢 燃料电池耗氢氢锅炉耗氢氢负荷储氢罐充氢。设备运行约束包括上下限约束和爬坡约束[ P_{i}^{min} \cdot u_{i,t} \le P_{i,t} \le P_{i}^{max} \cdot u_{i,t} ][ -\Delta P_i^{down} \le P_{i,t} - P_{i,t-1} \le \Delta P_i^{up} ]这里的 (u_{i,t}) 就是设备的启停0-1变量。之前有个读者跟我抱怨说他的模型总是无解后来一查发现是燃气轮机的状态变量只有0-1但没有给最小启停时间约束结果求解器给了一个“逻辑上不可能”的方案——设备在连续两个时段里频繁启停导致爬坡约束直接冲突。这个问题后面我还会细讲。储能动态约束主要针对蓄电池、储氢罐和储热罐[ SOC_{t1} SOC_t \eta_{ch} P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}} ][ SOC_{min} \le SOC_t \le SOC_{max} ][ SOC_{T} SOC_{0} ]最后这个SOC等式约束很关键——它规定了调度周期结束时储能状态必须回到初始值保证系统在一个调度周期内是“可持续”的。如果不加这个约束优化结果会把储能在最后一个时段全部放空虽然成本最优但完全不符合实际运行逻辑。2.4 氢能子系统的独特建模细节氢能这块和电池储能很不一样有几个容易忽略的点电解槽的效率不是常数而是和输入功率相关的曲线。很多论文为了简化用固定效率但实际做工程分析时我建议用分段线性化的效率曲线否则在部分负荷工况下制氢量的误差能到10%以上。简单处理可以用公式[ V_{H2,t} \eta_{EL}(P_{EL,t}) \cdot P_{EL,t} / HHV_{H2} ]其中HHV是氢气的高热值约等于39.4 kWh/kg。储氢罐的建模也需要特别注意。储氢压力与储氢量之间是近似线性的关系理想气体状态方程但储氢罐的工作区间受压力上限约束这要转化成储氢量的上下限[ m_{H2,t} \in [m_{H2}^{min}, m_{H2}^{max}] ]燃料电池发电效率则和输出功率正相关但存在一个最佳工作区间。负荷率低于20%的时候效率掉得很厉害所以我一般在约束里加一个最低出力限制避免求解器给出一个“发电1kW”这种工程上完全没有意义的解。3. 优化调度求解YalmipCplex/Gurobi的典型组合3.1 为什么选这个求解组合现在做混合整数线性规划MILP优化调度研究最主流的工具组合是Matlab Yalmip工具箱 Cplex或Gurobi求解器。这个组合的好处是Yalmip把建模和求解分离你写约束用的是人看的符号底层自动转换成求解器能认的标准形式。Cplex和Gurobi对MILP问题的求解效率远高于Matlab自带的intlinprog尤其当系统有数百个变量和约束时求解速度差距能达到几十倍。调试方便Yalmip会主动报错哪里维度不匹配、哪里定义了未声明的变量。如果你的Matlab版本不能用Cplex新版Matlab对Cplex兼容性确实一般Gurobi是很好的替代方案。我在2022b和2024a上都测过Gurobi 10.0以上的兼容性表现稳定。3.2 状态变量与连续变量的声明方式在Yalmip里变量的声明直接决定模型规模。我建议这样处理%% 定义时间尺度和参数 T 24; % 调度周期24小时 n_EL 1; % 电解槽数量 n_FC 1; % 燃料电池数量 %% 连续变量 P_PV sdpvar(1, T); % 光伏出力 P_WT sdpvar(1, T); % 风电出力 P_Grid sdpvar(1, T); % 电网购电 P_EL sdpvar(1, T); % 电解槽耗电 P_FC sdpvar(1, T); % 燃料电池发电 P_GB sdpvar(1, T); % 燃气锅炉产热 P_HB sdpvar(1, T); % 氢锅炉产热 SOC_BESS sdpvar(1, T1); % 蓄电池SOC m_HST sdpvar(1, T1); % 储氢量 %% 0-1变量 u_EL binvar(1, T); % 电解槽启停 u_FC binvar(1, T); % 燃料电池启停 u_GB binvar(1, T); % 燃气锅炉启停注意我习惯把储能类变量定义成 (T1) 维这样SOC初始值可以直接赋给第1个点然后从第t时刻到t1时刻的递推表达式写起来更方便不会出现索引越界的问题。3.3 核心约束的代码写法以电解槽为例它的运行约束在Yalmip里是这样组织的%% 电解槽约束 % 出力上下限 P_EL_min 0.2 * P_EL_max; % 最低运行负荷率20% Constraints [Constraints, P_EL_min.*u_EL P_EL P_EL_max.*u_EL]; % 爬坡约束 Constraints [Constraints, -delta_EL_down P_EL(2:T) - P_EL(1:T-1) delta_EL_up]; % 产氢量计算分段线性化效率 eta_EL 0.75 - 0.05 * (P_EL / P_EL_max); % 效率随负载率降低 V_H2_prod eta_EL .* P_EL / HHV_H2; % 单位kg/h这里有个小技巧分段线性化的效率表达式如果用变量直接算非线性项Yalmip会把它识别成非线性约束Cplex就没法直接求解。一种稳妥的做法是把效率设为常数比如0.7或者使用二进制展开法把非线性项线性化。对于硕士阶段的课题研究固定效率是可以接受的导师如果要严谨一点再上分段线性化。储氢罐的动态约束%% 储氢罐约束 m_HST sdpvar(1, T1); Constraints [Constraints, m_HST(1) m_HST_init]; Constraints [Constraints, m_HST_min m_HST m_HST_max]; % 氢气平衡单位kg % 产氢量 放氢量 燃料电池耗氢 氢锅炉耗氢 氢负荷 充氢量 % 注意这里的方向定义要统一 for t 1:T % HST充放氢用同一个变量hlf正为充、负为放 Constraints [Constraints, m_HST(t1) m_HST(t) V_H2_prod(t) - V_H2_fc(t) - V_H2_hb(t) - V_H2_load(t)]; end % 周期末SOC回到初值 Constraints [Constraints, m_HST(T1) m_HST_init];算例跑通之后我强烈建议你把功率平衡输出的各组成部分画在一张堆叠图上检查一下每个时段“供给需求”是否真的对上了。很多时候模型无解就是因为有个时段的等式约束差了一个常数项比如单位不一致MWh和kWh混了。4. 算例设计参数从哪来结果怎么解读4.1 基础数据怎么定做算例之前需要先确定系统的规模和数据来源。我的习惯是参考已有文献的测试系统或公开数据集负荷数据可以采用某园区典型日的电、热负荷曲线也可以自己构造带峰谷特性的曲线但一定要保证时序特征明显否则调度结果没有辨识度。新能源出力数据光伏可以按典型晴天的出力曲线构造早晚低、中午高风电用带有一定随机波动的时间序列。如果没有实测数据用德国Open Power System Data的公开数据是个不错的选择。分时电价按峰平谷三段设置比如峰时1.2元/kWh、平时0.8元/kWh、谷时0.4元/kWh具体数值根据研究场景调整。设备参数可以参考一些公开发表的综述文章中的典型参数。比如电解槽效率0.7~0.75燃料电池发电效率0.45~0.55蓄电池充放电效率0.95等。这里的核心逻辑是参数数值本身不需要极度精确但必须合理且自洽。审稿人更关心的是你模型结构对不对、机理是不是清楚而不是你的电解槽效率到底是0.71还是0.72。但如果你把电解槽最小运行负荷率设成0.8那模型基本就是摆设了——实际电解槽在20%~30%负荷率下仍然可以稳定运行。4.2 一个典型的算例结果长什么样我拿一个实际跑通的案例来说明。系统配置光伏200 kWp、风电150 kW、电解槽100 kW、燃料电池50 kW、储氢罐50 kg有效容量、燃气锅炉200 kW、蓄电池100 kWh/50 kW。在典型冬季日里调度结果呈现出几个有规律的现象凌晨0点到6点负荷低、风电出力高电价处于谷段。这个时段电解槽以接近满功率运行把多余风电转化为氢气储存蓄电池也在低电价时段充电。白天8点到16点光伏出力大但电价处于峰段。蓄电池开始放电燃料电池也在部分时段开启尽量替代高价购电。如果光伏有富余电解槽可以继续低功率运行制氢。傍晚17点到22点光伏出力下降热负荷升高进入晚高峰。此时储氢罐放氢燃料电池发电氢锅炉补充供热燃气锅炉作为支撑热源系统购电被抑制在较低水平。夜间23点后热负荷降低天然气价格相对低燃气锅炉和燃料电池的出力策略调整等待次日谷电时段重新制氢。如果你画出的调度图中电解槽的启停总是和谷电时段完美重合恭喜你说明模型求解逻辑是正常的。如果出现电解槽在峰电时段满负荷制氢的结果那就要检查电价序列是不是写反了或者惩罚项系数是否有问题。4.3 灵敏度分析与场景对比为了让研究更有深度我建议至少做两组对比有无氢能系统的对比Same系统分别在不含氢能只有电池和含氢能的情况下求解对比总成本、弃电率、碳排放。这个结果基本就是整个研究的核心卖点。储能容量灵敏度分析改变储氢罐容量和电解槽功率观察系统总成本的变化趋势。通常随着储氢容量增大成本先显著下降然后趋于饱和这个拐点就是最优配置点。这类分析用Matlab跑非常方便——写一个for循环遍历不同参数每次求解一次模型把结果保存下来就行。要注意的是每一次求解都要重新构建约束不要想着“只改一个参数复用之前的优化结果”那样会把前一组的变量和约束带进来导致结果错误。5. 工具选型与Matlab代码组织经验5.1 需要什么工具箱跑这套代码建议的环境是Matlab R2020a及以上版本新版本对Yalmip的支持更好Yalmip工具箱免费直接从GitHub下载添加到路径Cplex 12.10或Gurobi 10.0求解器学术版免费如果你的电脑装不了Cplex或者Gurobi退而求其次可以用Matlab内置的intlinprog。但要注意intlinprog对大规模MILP的求解速度确实不行如果模型的0-1变量超过100个求解时间可能会从几十秒变成几十分钟。5.2 代码结构怎么组织才不容易乱一个完整的优化调度项目我建议按下面的目录组织project/ ├── data/ % 存放所有输入数据 │ ├── load_curve.xlsx │ ├── pv_wt_data.mat │ └── price.m ├── scripts/ % 主脚本 │ ├── main_opt.m % 主程序入口 │ ├── build_model.m % 构建模型 │ ├── solve_model.m % 求解与结果处理 │ └── plot_results.m % 画图 ├── functions/ % 工具函数 │ ├── build_constraints.m │ ├── init_parameters.m │ └── calc_cost.m └── results/ % 输出结果 ├── figure/ └── data/这种结构最大的好处是如果参数变了你只需要改 init_parameters.m 里的值如果约束变了只需要改 build_constraints.m 里的对应块。Debug的时候也方便算错了知道去哪个文件找。5.3 求解效率优化的小技巧在实际调代码的过程中有几个提升求解效率的经验尽量少用for循环写约束能用向量化运算就用向量化。Yalmip的约束拼接对向量是友好的比如 P_min.*u P P_max.*u 直接一个表达式就完成整个时间段的约束。0-1变量越少越好。每个0-1变量都会增加分支定界的搜索空间。如果一个设备不需要启停优化比如燃气锅炉一直开着就没必要给它加0-1变量。给求解器设置合理的容差。在我的算例里MIP Gap设为0.5%通常足够再小只会拖慢速度而不会带来有意义的最优性提升。一个最实用的技巧是先在简化条件下调试比如把T从24改成4把设备数量减少确认模型逻辑正确后再跑完整算例。很多人上来直接跑24时段、20个设备的大模型一旦无解根本不知道是哪个约束的问题排查起来非常痛苦。6. 常见报错与调试记录6.1 模型无解Infeasible problem这是最常见的问题通常有以下几个原因某一时段的负荷值给得太高所有设备满出力也满足不了平衡方程。检查方式是看总供给上限是否大于最大负荷。0-1变量约束和连续变量约束冲突比如设备最小出力设得比最大出力还大或者启停变量被强制为0时连续变量又被要求大于某个正值。SOC等式约束写错了方向导致储能递推关系在时间轴上自相矛盾。单位混用比如电功率用了kW但储能容量用了MWh导致数量级差1000倍约束数值范围差距过大求解器数值稳定性变差。排查方法很简单把约束逐条注释掉从“所有设备自由出力”开始逐条加回约束看是哪一条加入后导致无解。Yalmip还会给出诊断信息提示“the constraints involving ... are infeasible”定位到具体那一块。6.2 求解器报错“License Error”或找不到求解器Cplex/Gurobi装好之后Yalmip偶尔会报找不到求解器。解决办法是检查Matlab路径里是否已经添加了求解器的安装目录然后在Matlab里运行yalmiptest这个命令会列出当前可用的所有求解器。如果Cplex/Gurobi没列出来就要检查环境变量或重新安装求解器接口组件。我之前遇到过一次Cplex在2024a里连不上的情况后来发现是需要以管理员身份重新运行安装脚本路径设置权限不够导致的。6.3 求解时间过长如果模型在半小时内出不来结果基本不是计算机性能问题而是模型规模或约束设置的问题。优先做两件事检查0-1变量数量。把不必要做启停优化的设备全部去掉0-1变量。检查是否有非线性的目标或约束混在模型里。Yalmip会把这类问题变成混合整数非线性规划MINLP求解难度骤然上升。如果真的需要非线性尝试用线性化手段比如分段线性化、大M法把问题重新表达为MILP。我碰到最夸张的一次一个看似很小的模型跑了近两个小时没出结果后来发现是储氢罐的压力—容量关系进了约束引入了非线性项。改成线性近似后求解时间直接降到40秒。7. 把代码灵活改造成自己的版本很多朋友拿到一套代码后第一反应是“我要换一个场景改哪里”这里我给出几个最常改的点7.1 怎么改负荷与新能源数据数据修改最简单直接刷新 data 目录下的输入文件就行。注意检查曲线的长度和采样周期如果原来是60个点10分钟间隔你改成24个点小时间隔那么所有T相关的变量都会自动匹配但需要确认爬坡约束、电价时段等参数是否还与原来的时间单位一致。7.2 怎么增加新设备增加新设备需要改四个地方参数区、变量声明、约束构建、目标函数。以增加一台微型燃气轮机为例% 参数区 P_MT_max 80; % kW P_MT_min 10; % kW eta_MT 0.35; % 发电效率 eta_MT_heat 0.45; % 余热回收效率 % 变量声明 P_MT sdpvar(1, T); u_MT binvar(1, T); % 约束构建 Constraints [Constraints, P_MT_min.*u_MT P_MT P_MT_max.*u_MT]; % 电平衡中加入P_MT % 热平衡中加入eta_MT_heat/(1-eta_MT)*P_MT*eta_waste_heat... 注意余热计算公式 % 目标函数中加入购气成本看起来简单但有一个容易错的地方燃气轮机的热电比是固定的电出力和热出力通过热电比耦合不能单独设定。所以增加燃气轮机时需要同时修改电平衡和热平衡两个等式约束并且在原方程中同步引入热电比参数。7.3 怎么改求解器相关配置Yalmip里改求解器很方便只需要在最后一行这样写optimize(Constraints, Objective, sdpsettings(solver, gurobi, verbose, 2));想要改MIPGap容差和最大运行时间可以在sdpsettings里加参数options sdpsettings(solver, gurobi, verbose, 2, ... gurobi.MIPGap, 0.005, ... gurobi.TimeLimit, 600); optimize(Constraints, Objective, options);每个求解器的参数名不太一样Gurobi用 gurobi.MIPGapCplex用 cplex.mip.tolerances.mipgap但Yalmip会自动把这些设置映射到对应求解器。8. 氢能调度研究还能怎么扩展这一节聊点我自己的体会。作为研究课题计及氢能的综合能源优化调度可以往几个方向继续做深一是多时间尺度优化。日前调度是小时级但氢能系统的响应秒级到分钟级都可以做到。把调度周期拉长、把时间粒度细化可以做“日前-日内滚动修正-实时调整”的多级协调框架这样就更有工程价值了。二是考虑不确定性。风光出力和负荷预测都有误差可以用场景法蒙特卡洛抽样场景削减或鲁棒优化来建模。氢能因为可以大规模长时间存储天然就是应对不确定性的优质资源。三是与碳交易机制结合。在目标函数里加入碳配额成本和碳交易价格让系统在做调度决策时主动倾向于低碳设备出力这个方向在“双碳”背景下非常热门。四是园区级或区域级的多主体博弈。多个微能源网之间可以共享氢能基础设施比如一个园区制氢另一个园区用氢主体之间的交互可以用纳什谈判或Stackelberg博弈建模。这些方向我大多在项目里试过matlab的Yalmip基本都能撑起来核心不会变——先把基础的电热氢耦合模型吃透再往外延伸就是水到渠成的事情。我自己刚上手时前后花了快两个月才把第一版模型完整跑通现在回过头来看其中一大半时间都花在了排错和调试参数上真正建模的思路梳理清楚之后代码反而是最容易的部分。希望这次分享的模型框架、约束写法、代码结构和调试经验能帮正在做或者准备做这个方向的同学少走一些弯路。特别是那些刚拿到一套代码不知道怎么下手的按照上面的思路先改数据跑通再逐个调整约束最后再换场景做对比分析基本一周内就能把一套代码吃透成自己的东西。
