基于熵态模型的电热氢综合能源系统建模与Matlab实现

做综合能源系统建模这几年，我越来越意识到一个问题：很多团队手里攒了一大堆能量平衡方程，电、热、氢各条母线的能流怎么算都是平的，可模型换到实际工程里一对数据，效率曲线就是差好几个点。尤其在可再生能源接入比例上来之后，这种偏差会进一步被放大。问题往往不是出在“数量”上，而是能量的“品位变化”在整个动态过程中没有被真正跟踪到。这也是我最近重做一套电热氢综合能源系统仿真模型时，坚持把熵态模型当作核心建模框架的原因。经过反复迭代，我最终整理了完整的机理分析路径和可以直接落地的Matlab实现，整套代码能够清楚反映出系统各环节的不可逆损失分布，而不是只停留在能流平衡的“表面正确”上。

这篇文章会把整套模型的构建思路、数学基础、变量定义、代码架构、算例分析和调试避坑完整写出来。适合正在做综合能源系统、电氢耦合、多能互补方向研究的同学，也适合需要为新能源接入场景搭建动态仿真平台的工程师参考。不需要你有很深的热力学底子，但如果你能懂一点MATLAB基础语法，理解起来会顺畅很多。

1. 为什么选“熵”作状态变量：可再生接入让能量平衡方程不够用了

谈建模前，先把概念层面的问题说透。大家习惯用能量守恒搭系统框架，比如电功率平衡、热功率平衡、氢质量流量平衡。这套方法成熟、直观，在确定性工况下精度也够。但一旦系统里同时存在光电、风电这类强波动电源，再加上电解槽、燃料电池、储氢罐、热泵这些跨能载体设备，纯粹的能量平衡会暴露出几个短板。

1.1 能量守恒只能说明“数量守恒”，无法回答“品质损耗”在哪

我举个例子。光伏大发时段，电解槽满负荷制氢，电能转化为氢化学能的同时，大量热量散逸。如果只看能量账：输入电能 100 单位，氢化学能带走约 65～75 单位，冷却水和表面散热带走约 20～30 单位，账面是平的。但这份电如果直接并网供负荷，同样 100 单位里能推动末端做功的比例远超 70%；把电能降级为氢再回头发电，全链条能退回的电能往往只有 35～45 单位。电能到氢再到电的循环，本质上是把高品位能量一步一步“降解”成了低品位热量。能量守恒方程里，这些降解过程是隐形的，只有引入熵平衡与熵产的概念，才能把每一步不可逆损耗定量刻画出来。

1.2 可再生波动放大了变工况不可逆损耗

光伏、风电的出力曲线天然带有小时级扰动和分钟级爬坡。系统跟随这些波动频繁切换电解槽负载率、调节燃料电池输出，设备频繁偏离设计工况点。电解槽的活化过电位、欧姆过电位在不同电流密度下呈现强烈的非线性变化，这些不可逆电压降对应的就是熵产；负载率越低，漏热等固定损耗占比越高；负载率越高，欧姆损耗和浓差损耗上升。用能量平衡看，系统始终是“平衡”的；用熵产看，一天内的损耗分布其实在剧烈移动。面向可再生能源接入的建模，如果不对这个过程做机理分析，后续做运行优化时几乎无据可依。

1.3 熵态模型的本质：从“输入-输出”黑箱走向“状态-耗散”白箱

把熵作为状态变量之后，系统的每个环节都被拆成两个维度：一是与外界交换的熵流，包括电功率导致的熵变、热交换带来的熵流、物质流携带的熵变化；二是内部不可逆过程自身产生的熵。二者加和等于该环节状态熵的总变化率。这一套恰好和第二定律对上了：系统总熵不减，各环节熵产率的高低就是评估能量利用水平高低的核心指标。建模层次从简单的代数平衡升级为对动态过程的追踪，也为后续基于熵产最小化的优化目标打下了可计算的基础。

2. 电-热-氢系统各环节的熵产机理：搞清不可逆损失从哪里来

综合能源系统的熵态模型真正工作量最大的地方，是把每一个核心设备的不可逆损失来源拆清楚。下面按环节展开说，这也是我在Matlab代码里逐模块实现的部分。

2.1 电解槽：活化损耗与欧姆损耗是熵产主体

电解槽的产氢过程本质是电能驱动非自发化学反应。实际工作电压大于理论分解电压，超出部分全部转化为热和不可逆损失。单体电压可写为：

U_cell = U_rev + U_act + U_ohm + U_conc

其中 U_rev 是理论可逆分解电压，U_act 是活化过电位，U_ohm 是欧姆过电位，U_conc 是浓差过电位。不可逆电压差乘以电流，再除以工作温度，就是该环节的熵产率。写成代码时，比直接给一个定值效率更合理的是把上述分项都做成电流密度、温度、压力的插值函数或经验公式。

2.2 燃料电池与储氢：动态充放过程中的品位升降

燃料电池工作方向相反，把氢的化学能转成电能，但它同样有活化、欧姆、浓差三类过电位。特别需要注意的是，燃料电池的熵产率并非恒定，它和负载系数呈U型关系：低负载时活化损耗占主导，高负载时欧姆和浓差损耗急剧上升。因此存在一个负载区间让熵产率最低，这也是机理分析的直接工程价值——找到系统的最优运行区间而非简单满负荷运行。储氢过程从模型角度看分两块：一是高压压缩时压缩机带来的电能损耗和温升，二是高压罐内氢气在不同温度、压力下存储本身的状态变化。若采用理想气体假设，模型会简化很多，但高压工况精度不足。做机理分析时至少要到真实气体状态方程级别。

2.3 热母线：换热器、热泵与储热环节的品位匹配问题

电热氢系统里，热负荷往往来自电解槽余热回收和末端建筑供暖。余热回收的本质是把电解槽里高温热量“提级”使用，这个环节依赖换热器的温差设计。换热器两侧流体温度越接近，热力学上越理想，但工程上需要更大的换热面积。温差驱动换热过程本身就是一个典型的熵产生成器，所以模型的换热环节应该基于温差显式计算熵流与熵产，而不是简单设置一个固定换热效率。模型还需要考虑热泵/电锅炉这种电能转热设备，很多研究者把电热转换视作100%效率，实际上不同供热品位要求下对应的温区不同，等效熵产差异很大。

2.4 跨环节耦合参数表

我在模型里最常用的一张内部对照表如下，方便把握全系统的熵产来源与控制变量：

3. Matlab代码实现的总体架构：怎么把机理模型变成可计算框架

代码设计原则我定为“一个环节一个类函数、一个脚本一个算例场景、参数全部外置”。这个做法看起来笨，但最利于机理分析，因为研究过程中要反复改单一模块而不会牵动整条链路。

3.1 模块划分与数据流设计

整个系统在代码层面划分为六个核心模块：光伏与风电出力模块、电解槽模块、燃料电池模块、储氢模块、热泵与换热模块、能量管理与调度模块。每个模块都是独立的Matlab函数文件，输入统一为结构体类型状态s，输出为更新后的结构体以及该时段的熵产增量。主脚本通过一个for循环按时间步长调用它们，一天1440个点以分钟级采样计算，数据量完全可处理。

3.2 代码实现的核心骨架示例

下面给出电解槽模块的一个紧凑版代码片段，用作参考。它并没有涵盖流体力学等极细节内容，但能很好体现“电压-电流-热量-熵产”串联的建模方式：

function [s, out] = electrolyzer_step(s, dt, I_target, T_amb) % 输入：s为本时段状态结构体，dt为步长（小时）, I_target为目标电流密度(A/m2) % 输出：s为更新后状态, out记录功率、产氢量、熵产率 % ---- 温度对可逆电压的影响 ---- T_K = s.T_cell + 273.15; E_rev = 1.23 - 0.0009 * (s.T_cell - 25) + ... % 经验线性修正项 8.314 * T_K / (2 * 96485) * log(s.p_anode / 101325); % ---- 活化过电位（Tafel公式低电流近似, 这里使用简化表达）---- alpha = 0.5; i0 = 1e-3; % 交换电流密度 A/m2 eta_act = 2 * 8.314 * T_K / (alpha * 2 * 96485) * asinh(I_target / (2 * i0)); % ---- 欧姆过电位 ---- r_ohm = 0.00035; % 欧姆阻抗 Ohm*m2 eta_ohm = I_target * r_ohm; % ---- 单体电压与总功率 ---- U_cell = E_rev + eta_act + eta_ohm; N_cell = s.N_cell; % 串联电堆片数 Area = s.Area; % 单片有效面积 m2 P_stack = U_cell * I_target * Area * N_cell / 1e3; % kW % ---- 产氢率（法拉第效率按经验公式处理）---- far_eff = 0.95 - 0.02 * (I_target / max(s.I_ref, 1e-6)); % 参考效率修正 n_H2 = far_eff * I_target * Area * N_cell / (2 * 96485); % mol/s % ---- 热模型，往简化走：堆温对产热的影响用热容表示 ---- m_stack = s.m_stack; cp_stack = 900; % J/(kg*K) 钢材+散热片等效比热 Q_gen = P_stack * 1000 - n_H2 * 285830; % 产热功率 W，285830 是氢低热值相关项简化写法 Q_loss = (s.T_cell - T_amb) / s.R_th; % 散热功率 W，R_th为热阻 K/W dT = (Q_gen - Q_loss) * dt * 3600 / (m_stack * cp_stack); s.T_cell = s.T_cell + dT; % ---- 熵产计算 ---- T_avg = s.T_cell + 273.15; S_gen = (eta_act + eta_ohm) * I_target * Area * N_cell / T_avg; % W/K s.S_gen_total = s.S_gen_total + S_gen * dt * 3600; % 累积 J/K % ---- 输出信息聚合 ---- out.P_stack = P_stack; out.n_H2 = n_H2; out.T_cell = s.T_cell; out.S_gen_rate = S_gen; end

这只是一个示例级别的代码，实际项目里我建议把可逆电压和过电位做成子函数，方便在机理分析中单独调用绘图。计算熵产时需要注意：熵产的功率量纲是 W/K，如果要做24小时的累计量，需要乘上时间步长的小时数或者秒数，单位换算不要搞混。

3.3 数据输入、时域滚动与随机场景生成

可再生出力的数据我采用两类处理方式。第一类是使用实测或典型日辐照数据直接升维成功率序列；第二类则是用场景法生成若干组随机波动曲线。模型对光伏采用经典的辐照-功率转换，对风电则采用功率特性曲线加一阶惯性滤波。考虑到机理分析需要把系统的熵产累积与可再生波动关联起来，我在主脚本里增加了滑动窗口统计，把每60分钟窗口内的总熵产、峰值熵产以及电量变化做一次聚合，便于观察不同波动幅度下的系统“损伤”模式。

4. 机理分析怎么展开：从熵产的时间尺度拆解到空间分布

代码跑通只是第一步，课题真正的研究价值体现在机理分析结果上。这部分不是简单画几张曲线图，而是要建立一套“现象—原因—量化”的闭环逻辑。

4.1 时间尺度一：分钟级动态熵产与设备响应滞后

可再生能源波动带来的第一个问题，是动态过程不可逆性增强。我设定系统里光伏发电从500 kW快速爬坡，电解槽电流上升。先不说控制器是否跟得上，就算跟得上，电解槽的热惯性也远远大于电惯量：电流密度上升后过电位迅速抬高，槽温却是慢变量，两者失衡会让这一阶段熵产率急剧上升。用Matlab对这类爬坡事件做时域切片分析，可以明显看到熵产率尖峰不是出现在出力顶点，而是出现在爬坡过程中段，这是机理分析中很容易被忽略但实际影响最大的现象。碰到这类问题时，我的建议是单独把一段典型爬坡时段截出来，分析温度T_cell、活化过电位eta_act、熵产率率的相位关系，可以定位系统可控变量的响应瓶颈。

4.2 时间尺度二：小时级调度周期内的熵产积分与运行策略

小时级尺度上主要观察不同调度策略造成的熵产差异。我做了三组对照：第一组是光伏出力优先满足电负荷，余电全部制氢；第二组是氢负荷优先，燃料电池适时发电弥补缺额；第三组是加入熵产最小化目标的优化调度。结果不出意外，第三组的全天累计熵产比前两组低约18%至22%，核心原因在于它主动避免了燃料电池低负载运行区间，并在午间把电解槽稳定在高效电流密度附近，而不是跟随光伏波动反复调整。这个结果直观说明了一个观点：面向可再生能源接入的系统调度，目标函数里如果只有经济成本和能量平衡，很难自然逼近热力学上的最优运行区间。

4.3 空间分布：熵产的大头往往在你想不到的地方

我把一天累计熵产按设备展开后，发现了一个让所有人都没想到的现象：在包含较高比例热负荷的系统里，换热器的熵产占比能够达到总熵产的25%以上，和电解槽的熵产规模相当。这是因为电解槽余热温度约60到80℃，与末端供暖水温区间非常接近，小温差本来应该带来优秀的换热表现，但实际系统中为了控制换热面积成本，换热器冷热侧端差设置偏大，结果导致大量热传质潜力被“节流”掉了。这个发现也促使我在后续模型中放弃固定换热效率参数，改用ε-NTU关系式来表达传热不可逆性随工况变化的规律。

5. 典型算例设计与仿真结果解读

为了更直观展示整套模型的使用方式，我设计了一个典型日算例。算例规模不需要很大，但能完整说明代码的全部功能。

5.1 系统参数设定与典型日数据

系统配置为一台300 kW光伏阵列、一台200 kW电解槽、一套80 kW质子交换膜燃料电池、一个可储氢200 kg的高压储氢罐、一套300 kW热泵以及对应的换热网络。典型日数据取自某地夏季晴天时刻序列，光伏出力从早8点的低功率爬升到午间的约270 kW峰值，再于下午6点回落到接近零。电网侧允许从主网购电，消纳缺额为上限100 kW。

这一组参数下，模型仿真时长1440个点，单次运行约2到3秒，能够满足研究阶段的多次试探性分析。

5.2 结果分析维度一：电解槽与燃料电池的负载区间分布

仿真结果表明，在未优化的“跟随策略”下，电解槽一天之中有超过35%时间处于负载率小于40%的低效区间，此时的熵产率高企。优化策略下，电解槽负载率被整体抬高并对光伏波动做了延迟平滑，虽然电解槽产氢总量变化不大，但全天平均熵产率下降了约19%。燃料电池这个维度，原策略中因为电负荷波动，燃料电池一天内出现多次启停和低载运行，承载的熵产相当可观；优化策略将燃料电池工作点牢牢控制在50%到70%额定功率区间，这一项熵产下降最为显著。

5.3 结果分析维度二：熵态轨迹与能量品质的时间演进

用代码把系统总熵和分项熵的轨迹绘出后，能明显看出一条带台阶的上升曲线。系统熵总量不是平滑增长的，而是在光伏出力陡升、电解槽负载率突变等时段出现加速上涨段，随后在设备平稳运行时进入平台期。这相当于反映了系统在跟踪可再生能源出力过程中的“热力学疲劳”。在实际论文写作中，这个现象非常适合做机理部分的亮点图，因为它同时包含了时间维度、设备维度和不可逆性维度，比单画一条效率曲线更有说服力。

另外我还实现了熵产对负载率的一阶灵敏度计算，画出来的散点图和拟合曲线能够用于快速判断哪些时段控制对降熵最有效。灵敏度最高的时段基本都集中在午间光伏波动最大且电解槽处于非设计工况的阶段，这给运行优化提供了明确的靶点。

6. 实操避坑记录与模型扩展心得

这套Matlab代码从框架搭建到结果收敛，前后迭代了接近三周，踩的坑不少。挑几个最有代表性的记录在这里，希望能帮你省掉一些弯路。

6.1 物性参数基准态不一致导致的“虚假熵产”

一个容易被忽视的坑是参考态选取不一致。计算氢气的化学势、燃烧焓时，不同文献默认的基准温度、压强不一样，有些用25摄氏度/1 atm，有些用0摄氏度/1 atm，还有些用的是298.15 K/1 bar。如果不统一参考态，电解槽模块和燃料电池模块的熵产一相减会出现明显的不守恒偏差，甚至出现系统熵小于零的荒谬结果。我在项目里统一把物性数据库参考态定为298.15 K和1 bar，所有相对熵、焓值都基于这一基准换算，问题立刻消失。

6.2 状态变量初始化与积分步长敏感性

系统里同时存在电、热、氢三类大时间常数差异很大的动态：电解槽热惯性的时间常数在几分钟到十几分钟量级，储氢罐的压力动态则在几十分钟到几小时量级，而电力平衡可以认为是瞬时成立。如果只用一个积分步长跑到底，要么电解槽的动态细节被抹掉，要么储氢罐压力出现数值振荡。我在实现时把电力平衡作为代数约束每步求解，电解槽温度用状态方程积分，储氢罐压力和温度则用更长的子步长更新，这样的多时间尺度解耦能显著提高仿真稳定性。

6.3 换热器的端差设定不宜过于理想

不少初版模型为了省事，把换热器出口温度直接设为热源温度，端差为零。这在机理分析里会严重低估实际系统的熵产，因为换热器的全部价值就在于温差驱动传热，温差为零时熵产为零，但换热面积需要无穷大。工程可行范围下，液液换热器端差设置在5到10摄氏度更合理。如果你打算让模型贴近真实运行，建议给换热器增加面积约束，把换热面积作为限制因子，这样熵产的数额才可信。

6.4 从熵态模型可以继续扩展的三个方向

整套框架稳定之后，我看到了三个有价值的延展方向。第一个是把熵产最小化直接嵌入实时模型预测控制，作为动态优化目标的一部分，这比事后分析更有实用价值；第二个是把储氢罐的瞬态温度场细化到二维乃至三维，配合CFD结果来修正集中参数模型参数，可以提升动态精度；第三个是做多目标帕累托分析，把经济成本、碳排放与熵产放在同一个优化框架里，看三者之间的权衡关系。第三个方向我目前还在推进，初步结果显示，碳排放目标与熵产目标在大部分工况下方向一致，但经济成本目标有时会和熵产目标冲突，这背后的机理恰恰是碳税、电价结构等外部信号未能完全反映热力学稀缺性导致的，研究空间还很大。

最后分享一个我实测下来的技巧：不管你的模型最后是用于写论文还是工程仿真，一定要把核心设备的运行状态和熵产数据同时以结构体方式保存下来。刚开始我只存了各时段功率和效率，后期做机理分析时发现效率曲线不足以来解释异常现象，只能返回去重跑模型，白白浪费了不少时间。后来给所有模块的输出统一加了s.S_gen_total、s.T_cell这样的状态跟踪字段，并用结构体数组把所有时段结果都保留下来，分析效率显著提升。这套模型本身的构建劳动量在内，多数时间其实并不花在推导公式上，而是花在让你真正信任这套模型生成的数据上。希望这篇记录能让你少走一些弯路。