做微电网优化调度这几年电池储能几乎成了标配但真正把能量在时间维度上搬来搬去的这两年明显转向了氢储能。标题里这个基于氢储能的热电联供型微电网优化调度方法说白了就是在一个同时有电负荷和热负荷的小型园区或社区电网里把燃气轮机、电解槽、储氢罐、燃料电池、燃气锅炉这些设备放在同一个优化框架下用Matlab求解每个时段各个设备该发多少电、产多少热、制多少氢、存多少氢。这篇文章就围绕这套代码把建模思路、约束处理、YALMIP加CPLEX求解、调试经验和参数敏感性一次性讲清楚适合正在做综合能源系统、微电网、氢电耦合方向的研究生和相关工程师参考。1. 项目背景与核心问题拆解1.1 氢储能为什么会出现在热电联供型微电网里热电联供系统的最大特点是电热耦合。燃气轮机或内燃机在发电的同时必然产生大量余热典型的热电比在1.0到1.5左右也就是说发1千瓦电同时要面对1到1.5千瓦的热量释放。传统控制方式下系统一般以热定电先把热负荷满足再让电功率跟随热功率走。问题在于园区电负荷和热负荷在一天里往往不同步白天电负荷高但热负荷可能不高早晚热负荷高但电负荷又不高于是系统要么多发电导致浪费要么得从电网买电补缺。电和热需求峰的错位让以热定电这条路走得非常别扭整个系统的灵活性被热电比这个硬指标卡住了。电池储能只能缓解电力侧的不平衡储热罐只能缓解热力侧的不平衡两者之间没有能量形态的转换本质上还是在各自的能量域里做缓冲。氢储能的介入方式完全不同它把电解槽、储氢罐和燃料电池串起来形成一个电-氢-电/热的转换链电价低谷时段电解槽把多余的电变成氢存起来电价高峰或者电负荷紧张时段燃料电池再把氢变回电同时余热回收供热水。原本必须在同一时段满足的电负荷和热负荷通过氢这个东西实现了时间上的平移这就是氢储能对热电联供系统最大的价值。当然氢储能不是万能的。电解槽效率一般在60%到75%燃料电池电效率在40%到60%如果算上热回收综合效率能到80%以上但只看电-氢-电的往返效率30%到45%是常态。这意味着用氢储能搬电跟用锂电池搬电相比单位能量的电损失大得多。它的优势场景是长时储能、大容量存储以及热电联供的综合利用而不是短时高频的调峰。在优化调度模型里必须把这种效率特性写清楚否则算出来的策略会严重失真。1.2 优化调度问题到底是什么混合整数线性规划框架把上述物理系统转成数学问题本质上是一个带约束的最小化问题。每个时段系统要决定从电网买多少电或卖多少电CHP机组发多少电、产多少热电解槽输入多少电、制多少氢燃料电池发多少电、回收多少热燃气锅炉补多少热储氢罐向各时段分配多少氢。目标函数是总运行成本最低约束包括功率平衡、热平衡、设备功率上下限、爬坡速率、储氢罐容量、购售电上下限等。这里面有一类特殊变量——设备的启停标志只取0和1。比如燃气轮机早上8点要启动启动后出力连续可调那么是否启动是个二进制变量而出力多少是连续变量。数学模型包含二进制变量、连续变量和一系列线性等式不等式这就是标准的混合整数线性规划MILP问题。选择MILP而不是启发式算法粒子群、遗传算法主要考虑三点一是MILP能给出全局最优解成本比较时不会因为随机种子不同而打架二是CPLEX、Gurobi这些商业求解器对MILP的求解已经非常成熟几十到上百个变量的微电网模型往往几秒钟就能出结果三是YALMIP这类Matlab建模工具箱把约束和目标函数写得跟数学公式几乎一一对应代码可读性好后续加碳排放约束、需求响应约束都很方便。从问题规模上看一个典型的24时段调度模型设备6到8个连续变量大概100到200个二进制变量十几个约束条件300到500条。这种规模对现代求解器来说属于热身运动难点不在于算不动而在于把物理约束建得准确、把数值问题处理好。2. 数学模型把物理问题翻译成求解器能解的方程2.1 系统架构、能量流与设备参数先交代一下我这里用的系统结构后面所有代码都围绕这张拓扑展开。系统接入外部配电网园区内部有一个燃气轮机CHP机组发电加余热回收供热、一个电解槽、一个储氢罐、一个燃料电池、一个燃气锅炉以及电负荷和热负荷两类负荷。电网和天然气管道是外部输入其余能量在系统内部流动。设备之间的能量关系可以用一张表说明设备输入输出典型效率/热电比关键运行约束CHP机组天然气电热电效率35%热电比1.2出力区间、爬坡、启停电解槽电氢少量热制氢效率65%输入功率区间、爬坡、启停储氢罐氢氢存储效率98%SOC区间、充放速率燃料电池氢电余热电效率50%热回收效率40%输出功率区间、爬坡、启停燃气锅炉天然气热热效率90%热出力区间电网联络线电电效率100%购售电功率上限这里给一套我常用的基准参数CHP额定电功率100kW热功率120kW电解槽额定输入60kW燃料电池额定电功率40kW热回收16kW按电功率的40%储氢罐容量40kgSOC允许范围0.1到0.9燃气锅炉额定热功率80kW联络线购电和售电上限200kW。这套参数对应一个中小规模的园区放在论文和工程里都比较有代表性。值得注意CHP机组电功率和热功率不是相互独立的而是被限制在一个可行区间内。实际中这个区间通常不是简单的矩形而是由燃气轮机运行特性决定的多边形。建模时可以用一组线性不等式来描述这个多边形或者用它的几个顶点坐标做凸组合。后者在代码里实现起来比较方便第3.3节会给出具体写法。2.2 目标函数成本怎么一项一项加起来优化目标是最小化一个调度周期内的总运行成本我通常写成四个部分之和。购电和售电成本按照分时电价计算买电为正成本、卖电为负成本也就是把电网看作一个可正可负的电源和负荷节点。这部分是分时电价套利的主要驱动力谷时段多买电制氢、峰时段燃料电池多发点电策略的核心逻辑就在这里。燃料成本CHP机组和燃气锅炉消耗天然气成本按天然气单价乘消耗量计算。CHP的天然气消耗量可以由电功率和电效率推导。这里为了保持线性模型要么把效率当作常数要么对效率曲线做分段线性化。常数效率虽然有一点点误差但在这个规模的问题里完全可以接受。运维成本每个设备的运行维护成本按输出功率或输入功率的线性比例估算数值上通常很小但考虑了以后模型更贴近工程实际。比如电解槽的运维成本按输入功率乘0.01元/kWh计入成本项虽然不起眼却能让优化结果避免出现为了省几分钱电费让设备频繁启停的极端策略。启停成本每次机组启动或停机产生一个固定费用用二进制变量与连续功率变量的乘积来表示。这个乘积是非线性的处理办法是引入指示变量用Big-M约束把启停和出力范围绑在一起写法在3.3节展示。目标函数写成表达式就是这样min sum(购电价 * P_grid) 气价 * 天然气消耗 运维系数 * 各设备功率 启停成本 * 启停变量所有项对每个时段累加。没有加碳排放成本如果想做双碳方向的研究在目标函数里加一个CO2排放系数乘以天然气消耗量和购电量的线性项代码改起来非常快这里不展开。2.3 功率平衡与约束条件分组建模约束条件我习惯按五组来组织代码里也按这个结构写。第一组是电功率平衡。每个时段进入系统的电等于流出系统的电P_grid(t) P_chp(t) P_fc(t) P_load(t) P_el(t)这里P_el是电解槽消耗的电功率。如果有电池储能右边还要加上电池充电功率减去放电功率但这套方案里暂时没有电池所以公式保持最简形式。注意电网功率P_grid可以为负负数表示向电网卖电。第二组是热功率平衡。因为CHP机组、燃料电池都带余热回收再加上燃气锅炉作为补热设备Q_chp(t) Q_fc(t) Q_boiler(t) Q_load(t)Q_fc在这里指燃料电池的回收热功率。热网没有储热罐所以热平衡必须逐时段满足这也是为什么热负荷波动大的场景里必须靠CHP和锅炉协调。很多人在初版模型里把热平衡漏掉直接用CHP满足所有热负荷的假设代替这样算出来的策略在工程上是不成立的。第三组是CHP电热运行约束。电功率有上下限热功率与电功率满足一个线性关系启停状态决定出力是否能非零。用数学语言说就是z_chp(t) * P_chp_min P_chp(t) z_chp(t) * P_chp_max Q_chp(t) alpha * P_chp(t) beta * z_chp(t)第二行的式子描述热电耦合alpha是热电比beta用来调整空载热损失。如果直接采用线性可行域多边形也可以写成顶点凸组合的形式后文代码里给一个简化版。第四组是氢储能动态方程。这个方程是整个模型的核心我建议用氢质量的平衡来写而不是像电池SOC那样直接写能量平衡。储氢罐内氢质量的变化等于电解槽产氢量减去燃料电池耗氢量M_h2(t1) M_h2(t) eta_el * P_el(t) / HHV_h2 - P_fc(t) / (eta_fc * HHV_h2)eta_el是电解效率eta_fc是燃料电池电效率HHV_h2是氢的高热值。这样写的好处是储氢罐容量直接对应到实际氢气质量后面做储氢量约束、容量规划都不用再转换。为了数值稳定两个系数都可以事先算成一个常数。比如1kWh电能制出多少kg氢1kWh电能对应的氢耗是多少kg都提前算好写进参数表。储氢罐自身还有容量上下限和充放速率限制同时为了支持周期性调度通常还会加一个末时段储氢量等于初始储氢量的约束。这个约束看起来很合理但实务里是导致不可行的一大源头第4.1节专门讲。第五组是设备爬坡约束、联络线功率约束和设备出力边界。爬坡约束本质上是相邻时段功率差值的上下限常见写法是-ramp_chp P_chp(t1) - P_chp(t) ramp_chp对于启停频繁的电解槽和燃料电池爬坡值往往也是影响调度策略的关键参数。很多人在初版代码里会漏掉爬坡约束结果优化结果在相邻时段出现剧烈跳变这在物理上根本实现不了。3. Matlab代码实现从建模到求解的完整链路3.1 代码组织怎么搭一个能跑通、能改参数的工程拿到这套问题时我建议不要把所有逻辑挤在一个脚本里。我自己的代码通常按功能分文件大致结构如下文件作用main.m主程序按顺序调用各模块load_parameters.m集中定义所有设备参数、负荷曲线、分时电价build_problem.m用YALMIP构建决策变量、约束和目标函数返回优化问题solve_and_report.m调用求解器输出结果并生成图表plot_results.m绘制电/热功率平衡图、储氢SOC曲线、成本结构图load_parameters.m里把所有可调参数聚在一起做敏感性分析、改场景时只需要动这一个文件。负荷曲线和电价用数组存比如T24P_loadzeros(1,T)然后按照典型日曲线填数。要注意统一单位我习惯功率全用kW时间粒度1小时氢量用kg成本用元。如果时间粒度改成15分钟所有能量相关的系数都要相应缩放这是最常见的一个低级坑。main.m的核心流程就五六行加载参数、构建问题、配置求解器、求解、画图。这样分模块的好处是后面换场景、加需求响应、加碳排放约束都不需要改动主流程定位问题也快。3.2 用YALMIP定义决策变量YALMIP是Matlab环境下非常好用的建模工具箱它最大的优点就是变量定义和约束输入跟数学公式几乎一一对应。连续变量用sdpvar二进制变量用binvar约束直接用中括号拼起来目标函数直接写表达式。代码维护成本很低换求解器的时候也只需要改一行配置。以24时段模型为例关键决策变量的定义代码是这样T 24; P_grid sdpvar(1, T, full); % 联络线功率正为购电负为售电 P_chp sdpvar(1, T, full); % CHP电功率 Q_chp sdpvar(1, T, full); % CHP热功率 P_el sdpvar(1, T, full); % 电解槽输入电功率 P_fc sdpvar(1, T, full); % 燃料电池输出电功率 Q_fc sdpvar(1, T, full); % 燃料电池回收热功率 Q_boiler sdpvar(1, T, full); % 燃气锅炉热功率 M_h2 sdpvar(1, T 1, full); % 储氢罐氢质量 z_chp binvar(1, T, full); % CHP启停 z_el binvar(1, T, full); % 电解槽启停 z_fc binvar(1, T, full); % 燃料电池启停有些初学者会问为什么M_h2要多定义一位T1因为氢储能动态方程需要从t时刻过渡到t1时刻末时段的状态也得有地方放。T1个变量让初值和末值约束直接落在两个明确变量上代码写起来更清爽也方便后面做滚动优化的时候对接上一轮的状态。3.3 建模核心代码约束、目标与线性化约束构建的核心过程是先设一个空Cell数组Constraints {}然后把每一类约束往里加。以启停与大功率范围的绑定为例这是MILP建模里最经典的Big-M写法P_chp_min 0; P_chp_max 100; P_el_min 0; P_el_max 60; P_fc_min 0; P_fc_max 40; Constraints{end1} P_chp_min * z_chp P_chp P_chp_max * z_chp; Constraints{end1} P_el_min * z_el P_el P_el_max * z_el; Constraints{end1} P_fc_min * z_fc P_fc P_fc_max * z_fc;这一组约束的含义是如果z等于0则对应功率必须等于0如果z等于1则功率被限制在上下限之间。它把设备是否运行和出力范围两个信息绑定在一起比单独限制功率区间严谨得多。注意这里Big-M参数我直接用了各设备的最大功率这是最紧、最有利于求解的取值。CHP的电热耦合如果采用可行域顶点法可以先把电热运行域的顶点表示出来然后用一个0到1之间的辅助变量做凸组合。简化版可以直接写alpha_CHP 1.2; % 热电比 beta_CHP 2; % 空载热损失系数 Constraints{end1} Q_chp alpha_CHP * P_chp beta_CHP * z_chp;电解槽和燃料电池的耦合关系类似燃料电池的热回收功率取电功率的固定比例eta_fc_heat 0.4; Constraints{end1} Q_fc eta_fc_heat * P_fc;氢储能动态方程写成% 预先算好系数每kWh电能制氢质量、每kWh电能对应氢耗 k_el 0.65 * 0.033; % 电解效率65%每kWh电能对应约0.033kg氢 k_fc 0.033 / 0.5; % 燃料电池电效率50%每kWh电能对应约0.066kg氢 Constraints{end1} M_h2(2:T1) M_h2(1:T) k_el * P_el - k_fc * P_fc;这个写法的好处是如果时间颗粒度从1小时缩到15分钟只需要把k_el和k_fc乘以0.25其他公式完全不用动。目标函数部分先把成本各分项累加%% 分时电价假设为峰谷两段 price_buy [0.3 * ones(1, 8), 0.1 * ones(1, 8), 0.3 * ones(1, 8)]; % 元/kWh示例 price_sell price_buy * 0.8; Cost_purchase sum(price_buy .* max(P_grid, 0)); % 购电成本 Cost_sell -sum(price_sell .* min(P_grid, 0)); % 售电收益 Cost_gas gas_price * sum(P_chp / eta_chp Q_boiler / eta_boiler); Cost_om sum(om_chp * P_chp om_el * P_el om_fc * P_fc om_boiler * Q_boiler); Cost_start sum(start_cost_chp * z_chp start_cost_el * z_el start_cost_fc * z_fc); Objective Cost_purchase Cost_sell Cost_gas Cost_om Cost_start;这里没有采用更精细的启动时刻捕捉逻辑直接用启停变量乘以单位启停成本属于一种近似只要设备在运行就当它产生了成本。对容量规划或粗略策略评估来说足够了但要做精确的启停次数统计会有偏差。想做得更精细需要引入额外的启动动作二进制变量并加状态转移约束代码量会增加不少考虑到微电网设备启停次数本身不多这种近似在工程里被广泛接受。3.4 求解器配置与结果输出建模完成后调用求解器的代码非常简洁ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 0.001; ops.cplex.timelimit 600; sol optimize(Constraints, Objective, ops);这里我把MIP gap设置为0.1%也就是求解器只要找到的解与最优解差距小于0.1%就停。这样既保证精度又能避免在一些数值困难的模型上死磕。如果求解器换成Gurobi配置方式是ops.gurobi.MIPGap 0.001参数结构略有不同需要留意。代码跑完以后结果输出至少画三张图电功率平衡堆叠图每个时段各设备发用情况、热功率平衡图CHP/燃料电池/锅炉各承担多少热负荷以及储氢罐SOC曲线。这套代码跑完能很直观看到电解槽集中在电价低谷时段工作、燃料电池在电价高峰时段发电这就是分时电价驱动下的低买高用逻辑。4. 踩坑实录与调试经验4.1 不可行问题的快速定位方法优化调度代码90%的运行失败都是infeasible problem也就是约束之间自相矛盾。定位方法我一般按三步走。第一步把目标函数临时改成常数0让求解器只做可行性检查很多情况下能更快暴露矛盾。第二步在Constraints里按组注释掉氢储能动态方程、爬坡约束、末值约束逐一排除嫌疑。第三步如果还找不到就检查负荷数据比如热负荷曲线在某个时段为负数、电价数据有NaN、负荷峰值超过所有设备容量之和这类低级错误是最容易被忽略的。实务里最常见的不可行原因是末时段储氢量等于初始储氢量与储氢罐容量上限冲突。比如白天燃料电池把氢耗掉60%但储氢罐初始SOC设得过高导致无论怎么调度都无法在周期末回到初值。解决办法是检查初值设置或者把末值约束从严格相等放松为介于初始值的95%到105%之间。对于日调度模型连续规划运行建议采用滚动窗口不用死磕末值相等。4.2 数量级差异与Big-M的选择问题MILP求解器对数值尺度非常敏感。这个模型里功率是几十到几百kW氢质量是几kg到几十kg成本是几千元三者放在同一个模型里如果处理不好求解器内部数值误差会明显上升甚至出现明明有最优解但求解器判断不可行的情况。我的经验是尽量保持变量在同一数量级。比如储氢罐SOC不用氢质量而用归一化的0到1或者给氢储能动态方程除以储氢容量让状态量落在0到1之间。Big-M参数也要尽量取紧的值比如启停约束里的上限就用该设备的最大功率不要图省事统一用10000。松散的Big-M会让线性松弛质量变差分支定界树巨大求解时间倍增。这个问题在初学者代码里出现频率非常高。我做过一个对比实验同一个24时段模型Big-M从设备最大功率改成10000CPLEX求解时间从0.8秒变成23秒MIP gap收敛曲线也明显变差。所以数值尺度不是小事值得在建模阶段就注意。4.3 求解时间过长时怎么取舍如果模型规模加大、二进制变量变多求解时间会指数增长。一个典型的例子从24时段扩展到96时段设备数量不变但二进制变量变成原来的4倍CPLEX默认配置下可能从秒级变成分钟级甚至更久。这时候有几个实用手段。一是设置一个合理的MIP gap比如0.5%或1%牺牲一点点最优性换来时间大幅下降。二是给求解器传入初始可行解比如用启发式规则生成电解槽在谷时段固定运行、燃料电池在峰时段固定运行的简单策略作为MIP start。三是在模型层面削减二进制变量个数比如把燃料电池和电解槽的启停变量合并成同一个变量因为它们确实不会同时运行或者把设备的启停时间窗口固定化只优化连续变量。另外还要提一句Matlab环境下YALMIP建模本身有开销如果T8760做全年调度变量和约束的构建时间比求解时间还长。代码层面可以考虑把模型写到结构化稀疏矩阵中或者直接用YALMIP的optimizer函数做预测控制式滚动调用。这一块展开讲篇幅很大但知道有这个方向对做长期调度的人很有帮助。4.4 参数敏感性分析哪些参数决定了调度策略代码跑通以后一定要做参数敏感性分析这既是论文需要也更贴近工程实际。我测试下来影响最大的三个参数依次是电解槽效率、氢气价格、储氢罐容量。电解槽效率的影响不用多说效率每下降5个百分点调度结果里电解槽的利用率会明显下降系统倾向于更多依赖电网购电而不是自产氢。氢气价格的敏感性更微妙气价很低时系统甚至可能在谷时段制氢、峰时段不发电而是直接用氢气价很高时燃料电池几乎不开机氢储能退化成纯耗能设备。这两个场景从优化结果图上能看得很清楚。储氢罐容量的影响体现在SOC曲线的形态上小容量储氢罐的SOC会在一天内频繁触碰上下限调度策略被迫短平快大容量储氢罐则可以把便宜的谷时段电能大量存储SOC曲线呈现明显的充电-放电两段式。理解了这种规律再回头设计容量配置和运行策略心里就有底了。最后给还没开始动手的人一个工作流建议先用一组平庸参数把模型跑通不要贪多设备、不要加太多约束跑通以后再逐个加爬坡、加启停、加末值约束每加一组就保存一个可复现的case。这样出了问题永远知道是哪里引入的比一次性写完几百行代码再回头查要省力得多。我早期吃过这个亏一口气写完整个模型结果一个约束没设对排查花了两天。做这套模型给我最大的感受是氢储能本身的经济性目前还不太好单纯靠分时电价套利通常算不出回本但把它放到热电联供系统里燃料电池的余热回收、电解槽参与需求响应、储氢罐作为长时备用这些附加价值加在一起模型的调度策略才真正有意义。如果你正准备复现这套代码我的建议是从三台设备的小系统开始先把电功率平衡和氢储能动态方程吃透再逐步加上CHP和热平衡每一步都对一下物理直觉别急着堆复杂度。这样调出来的代码不管以后做鲁棒优化、多目标还是容量规划底座都是稳的。