A Port-Hamiltonian Modeling Approach for Integrated Hydrogen Systems

A Port-Hamiltonian Modeling Approach for Integrated Hydrogen Systems Abdullah Shahin1, Hannes Gernandt1,2, Anton Plietzsch1, and Johannes Schiffer1,3 Abstract— Hydrogen’s growing role in the transition towards climate-neutral energy systems necessitates structured modeling frameworks. Existing gas network models, largely developed for natural gas, fail to capture hydrogen systems distinct properties, particularly the coupling of hydrogen pipes with electrolyzers, fuel cells, and electrically driven compressors. In this work, we present a unified systematic port-Hamiltonian (pH) framework for modeling hydrogen systems, which inherently provides a passive input-output map of the overall interconnected system and, thus, a promising foundation for structured analysis, control and optimization of this type of newly emerging energy systems. I. INTRODUCTION The shift to climate-neutral energy systems requires integrating various energy carriers, with hydrogen playing a crucial role in decarbonization [1]. Yet, for hydrogen to indeed become a key energy carrier, a dedicated generation and transportation infrastructure along with suitable operation and control structures must be established [2]. Such development is, e.g., already underway in Germany, where a “hydrogen core network” is currently being implemented with planned commissioning in 2032 [3]. For climate-neutrality, the hydrogen production needs to be “green”, i.e., based on renewable energy [4]. Likewise, hydrogen is also foreseen as an important carrier for providing long-term energy storage to electric power systems [2]. This results in a bidirectional energy conversion between electrical and hydrogen domains, so that these systems exhibit complex interdependencies that require structured, modular modeling approaches. Hence, hydrogen network components, which include production, storage, and distribution, must interact seamlessly with electrical grids to facilitate efficient and climate-neutral energy provision. Clearly, many new control and operational challenges arise for these type of networks. A first step towards addressing these challenges in a structured manner is the development of suitable modeling procedures for systematically assembling This work was supported by the European Union’s Horizon Europe Framework Programme (HORIZON) under the GA n. 101120278 - DENSE. HG and JS acknowledge funding from the German Federal Government, the Federal Ministry of Education and Research, and the State of Brandenburg within the framework of the joint project EIZ: Energy Innovation Center (project numbers 85056897 and 03SF0693A) 1Fraunhofer IEG, Fraunhofer Research Institution for Energy Infrastructures and Geotechnologies IEG, 03046 Cottbus, Germany, {abdullah.shahin, anton.plietzsch}@ieg.fraunhofer.de, 2Institute of Mathematical Modelling, Analysis and Computational Mathematics IMACM, University of Wuppertal, 42119 Wuppertal, Germany, [email protected] 3Brandenburg University of Technology Cottbus-Senftenberg, 03046 Cottbus, Germany, [email protected] mathematical models of hydrogen systems that reliably capture their fundamental dynamics. This is the main objective pursued in the present paper. In this context, port-Hamiltonian (pH) systems offer a powerful framework to represent energy-based multi-domain systems, while preserving structural and systemic properties, such as interconnection patterns and passivity [5]. The pH framework has been used successfully to model electric power [6]–[8], heat [9]–[11] and gas networks [12]. In the gas network domain, recently, in [13] a pH representation of lumped gas pipe models and their interconnection has been developed. But therein the focus lies on pipeline dynamics for natural gas systems, omitting hydrogen storage, generation, and load integration, such as electrolyzers and fuel cells. Building on [13], in [14] compressors were introduced into the gas grid model. In [14] the control theoretic properties of a lumped compressor-gas network were analyzed through the establishment of equilibrium-independent passivity (EIP). Additionally, a pH PDE formulation of the gas network including compressors and based on the Euler equations was derived in [15]. Their model explicitly accounts for the enthalpy added to the gas in the compression process, providing a detailed representation of the involved thermodynamics. Only a few recent studies have considered hydrogenelectric systems within the pH framework. The existing works [16] and [17] primarily focus on the component-level modeling of PEM fuel cells and electrochemical processes, offering insight into the internal dynamics and energy conversion behavior of these individual units. However, these models are not readily extendable to network-scale applications, where the interaction between multiple components and energy carriers must be coherently represented. In light of this, the main contribution of our work is to adopt a system-level perspective and systematically develop a unified pH model that captures the dynamics of interconnected hydrogen infrastructures. To this end, we first derive individual pH representations for the components of the core hydrogen system, such as pipes, electrolyzers, fuel cells, and electrically driven compressors. Next, by using algebraic graph theory and Kirchhoff’s law to describe the network interconnections, we provide an overall pH model, which ensures structure-preserving interconnections among system components. In this way, the presented approach reveals that hydrogen systems admit a compositional pH structure, which inherently provides a passive input-output map of the overall interconnected system and, thus, a promising foundation for structured analysis, control, and optimization of this type of arXiv:2512.02717v1 [math.OC] 2 Dec 2025 newly emerging energy systems. The remainder of this paper is structured as follows: In Section II, we recall the class of pH systems. The pH models of the hydrogen grid components are presented in Section III. In Section IV the interconnection and integration of hydrogen components in an interconnected system is described. Notation. We employ the following short-hands. For ℓ > 1 scalars vi, i = 1, . . . , ℓ, v = (vi)ℓ i=1 ∈Rℓdenotes a column vector, whose entries are vi.For v∈Rℓand w∈Rℓ′, where ℓ′>1, we denote (v, w)∈Rℓ+ℓ′as a column vector containing the entries of vand w. Given ℓ>1matrices Mi∈Rni×mi,i= 1, . . . , ℓ, the notation diag(Mi)ℓ i=1 ∈ R(Pℓ i=1 ni)×(Pℓ i=1 mi)denotes the block matrix, whose block elements are Mi. The identity matrix of size k×kis denoted by Ik. Furthermore, the zero matrix of size k×lis denoted by 0k×land in the case l= 1 we simply write 0k∈Rk. II. PORT-HAMILTONIAN SYSTEMS In the following, we recall from [18] the form of pH systems employed throughout the paper ˙x= (J(x)−R(x))∇H(x)+Bu +d, y=B⊤∇H(x) + Du, (1) where H:Rn→[0,∞)is the Hamiltonian, which is assumed to be continuously differentiable with gradient ∇H(x),B∈Rn×mis the input matrix, D∈Rm×mis the feedthrough matrix and d∈Rnis a constant disturbance. Furthermore, the interconnection and damping matrices J: Rn→Rn×nand R:Rn→Rn×n, respectively, satisfy the following structural conditions for all x∈Rn, J(x)=−J(x)⊤R(x) = R(x)⊤≥0, D +D⊤≥0. This directly implies that in the unperturbed case, i.e. d≡0, the system (1) is passive [18], i.e. the time derivative of the Hamiltonian Halong the solutions of the system (1) satisfies ˙ H=−∇H⊤R(x)∇H+∇H⊤Bu =−∇H⊤R(x)∇H+y⊤u−1 2u⊤(D+D⊤)u≤y⊤u. III. PORT-HAMILTONIAN MODELING OF HYDROGEN SYSTEMS In this section, we introduce pH models for pipes, compressors and storage units and also develop pH representations for electrolyzers and fuel cells. A. Network topology The hydrogen network topology is described by an undirected, connected graph G= (N,E)with set of nodes N={n1,...,nN}and set of edges E={e1,...,eM}. To each node ni∈ N, we associate a node pressure pi∈R. Moreover, to each edge eℓ={ni,nk}∈E, we associate a volumetric flow rate qℓ∈Rand assign an arbitrary direction to each edge. Then, if qℓis directed from nito nk, we call ni the source of eℓand nkthe sink of eℓand we write eℓ=−−→ nink. It is convenient to introduce the node-edge incidence matrix BG= (bij)∈RN×M, which is defined by bij =     1if node niis the sink of edge ej, −1if node niis the source of edge ej, 0otherwise. (2) We consider a hydrogen system that contains S≥1 storage units, F≥1junctions, and C≥1compressors. The set of nodes Nis thus subdivided into the set NS={n1,...,nS}that represents storage units, the set NF={nS+1,...,nS+F}that represents junctions within the hydrogen grid and the set of compressors NC= {nS+F+1,...,nS+F+C},with nS+F+C= nNand N= NS∪NF∪NC. The hydrogen system is coupled to the electric domain via E≥0electrolyzers and L≥0fuel cells. The interconnection topology between the electric power system and the hydrogen system is modeled by a second undirected graph GP= (NS∪NP,EP).As commonly done in practice, we assume that each power conversion device is always connected to a dedicated storage unit on the hydrogen side, i.e., to a node ni∈ NS.Hence, E+L≤S. Consequently, the corresponding buses in the electrical power system are represented by the set NP={nP,1,...,nP,E+L}. To each nP,i ∈ NP, we associate a voltage vi∈R. Furthermore, each pair of nodes {ni,nP,i}is interconnected by an edge eP,k ={ni,nP,i} ∈ EPwith EP={eP,1,...,eP,E+L}, where the edge eP,k ∈ EPrepresents the volumetric flow rate between a fuel cell or an electrolyzer and the hydrogen grid, respectively. In the sequel, we focus on developing the main components of a sector coupled hydrogen system. The vector of pressures of the hydrogen system is denoted by p= (pi)N i=1 and that of the volumetric flow rates by q= (qi)M i=1.Moreover, we introduce the vector of exogenous volumetric flow rates qex = (qex,i)S+F i=1 ∈RS+F, that is used to represent exogenous hydrogen injections or demands stemming, e.g., from onshore terminals. n4 n1 nP,1 n5 n2 nP,2 n6 n7n3 e1q1 eP,1qsc,1 e3q3 eP,2qsc,2 e2 q2 e4 q4 e5q5 e6 q6 qex,3 qex,4qex,5qex,6 Fig. 1. Example hydrogen network topology. The symbols denote the following elements: storage unit NS, junction NF, compressor NC, electrolyzer NP, fuel cell NP. Arrows indicate the volumetric flow rates in pipelines q1,q2,q3,q4, compressors q5and q6, electrolyzer outlet qsc,1, fuel cells inlet qsc,2and exogenous volumetric flow rates qex,i. B. Pipe model The hydrogen pipe associated with the ith edge ei= {nl,nr}∈E,nl∈ NS∪NF,nr∈ NS∪NF,is considered to be a spatially one-dimensional object and the temporal change of pressure and volumetric flow rate of the hydrogen along the pipe can be described in terms of the isothermal compressible Euler equation, see e.g. [12]. A simplified lumped model which models the volumetric flow rate qi at standard conditions and with a constant compressibility factor was described in [13, Proposition 2] and is given by ρ˙ qi Ai =pl−pr Li−ˆ λi|qi|qi−gsin(θi) c2pM,i,(3) with pressures pland prat the left and right ends of the pipeline segment, respectively, the cross-sectorial area Ai> 0, the diameter Di>0, the pipe-length Li>0, the speed of sound in hydrogen c > 0, hydrogen density at standard conditions ρ>0, a friction coefficient ˆ λi=λc2ρ2 2DiA2 ipM,i , with Darcy friction factor λ>0, pipe inclination angle θi∈ [−π/2, π/2] and the gravitational acceleration g > 0. Moreover, pM,i >0is the mean pressure across the ith pipeline segment given by the Weymouth mean pressure pM,i =2 3 p3 l−p3 r p2 l−p2 r =2 3pl+pr−plpr pl+pr.(4) To model the hydrogen flow rate within a pipe, we consider the following standard modeling assumptions, see also [13]. Assumption 1: In the model (3), the inclination angle θi and the mean pressure pM,i given by (4) are constant. With Assumption 1 and by introducing the states, Hamiltonian and co-states xe,i =ρLi Ai qi, He,i(xe,i) = Ai 2Liρ∥xe,i∥2,∇He,i(xe,i)=qi, the pipe model (3) can be cast in the following pH form, see also [13, Theorem 3], ˙xe,i = (Je,i −Re,i(xe,i))∇He,i(xe,i)+Be,iue,i +de,i, ye,i =B⊤ e,i∇He,i(xe,i)=qi, ue,i =pl−pr,(5) de,i =−gLisin(θi) c2pM,i, Be,i = 1 Je,i = 0, Re,i(xe,i) = ˆ λiLi Ai Liρxe,i. As shown in [13], the lumped model (5) achieves an accuracy comparable to advanced discretized pipeline models [19]. Also, with Assumption 1, Re,i(xe,i)≥0holds for all xe,i ∈R. In comparison with [13], we do not incorporate the pressure dynamics of the boundary nodes of the pipe directly into the model (3). Instead, we incorporate these pressure dynamics in the junction and storage unit models, since this simplifies the modular interconnection of the individual network components. C. Hydrogen storage unit model In this section, we describe a general model of a hydrogen storage unit connected at node ni∈ NSwith pressure piand which is interconnected via two edges el∈ E and ek∈ E. Then, the dynamics is given by [20] MH2Vs,i ρRTs,i ˙ pi=−ρ−1rs,ipi+qin,i −qout,i +qex,i,(6) where Vs,i >0is the storage unit volume, MH2>0is the molar mass of the stored hydrogen, qin,i is the ingoing volumetric flow rate from the grid to the node ni,qout,i is the outgoing volumetric flow rate from the node niinto the grid, qex,i is a volumetric flow rate that represents exogenous hydrogen injections or demands, Ts,i >0is the temperature of the stored hydrogen, which is assumed to be constant, and R= 8.314 J/(mol ·K) is the universal gas constant. Moreover, to account for potential storage unit losses, such as leakage, we introduce a dissipation constant rs,i >0. To write the dynamics (6) in pH form, we introduce xn,i =MH2Vs,i ρRTs,i pi, Hn,i(xn,i) = 1 2 ρRTs,i MH2Vs,i x2 n,i, Rn,i =rs,i ≥0, Jn,i = 0, Bn,i =1 1, and define the output as well as the input yn,i =B⊤ n,i∇Hn,i(xn,i) = pi pi, un,i =qin,i −qout,i qex,i . With this, we obtain from (6) the following pH model ˙xn,i =−Rn,i∇Hn,i(xn,i)+Bn,iun,i, yn,i =B⊤ n,i∇Hn,i(xn,i).(7) D. Junction model Based on the considerations in [13] we describe a junction between different pipes at a node ni∈ NFby the following equation for the node pressure Ci˙ pi=qin,i −qout,i +qex,i, Ci= M X l=1,ni∈el LlAl 2ρc2, where Ci>0models a lumped storage capacity at the junction node niand we sum in the definition of Ciover all edges that are incident with the node ni. Our junction model can be viewed as a lossless storage with a comparably small storage capacity Ciand accordingly, the pH formulation of junctions is given similar to (7) with xn,i =Cipi, Hn,i(xn,i) = 1 2C−1 ix2 n,i, Rn,i = 0, Jn,i = 0, Bn,i =1 1. (8) E. Compressor unit model We consider a centrifugal compressor at node ni∈ NC, which is modeled according to [21]. That is, we describe the pressure dynamics within the plenum piand using the input pressure plat the inlet grid node nl∈ NS∪NFand the pressure prat the outlet node nr∈ NS∪NFthat receives the volumetric flow rate from the plenum pi. This means that the corresponding volumetric flow rates qftowards the plenum and from the plenum to the grid, qm, are represented via the volumetric flow rates over the edges ef={nl,ni} ∈ E and em={ni,nr}∈E. Thus, the dynamics are Vp,i ρa2 01,i ˙ pi=−ρ−1rpl,ipi+qf−qm, ρLc,f A1,f ˙ qf=p2,f −pi, ρLo,m Ao,m ˙ qm=pi−pr, (9) where a01,i >0is the inlet stagnation sonic velocity, Vp,i > 0is the volume of the plenum, Lc,f >0is the length of compressor and duct, Lo,m >0the length of the compressor outlet, A1,f >0,Ao,m >0are the reference areas, p2,f is the output pressure of the compressor and rpl,i >0is a coefficient accounting for pressure losses in the plenum. Adherent to the methodology proposed in [15], the differential pressure across the compressor ∆pjis used as a control input, i.e., p2,f =p2,f −pl+pl,∆pi−N+C:= p2,f −pl,(10) where the node index ifulfills i∈ {N−C+ 1, . . . , N}and therefore i−N+C∈ {1, . . . , C}. In the following, we split the compressor model into three pH submodels: the pressure node dynamics in the plenum which can be modeled as a storage node (7) with xn,i =Vp,i ρa2 01,i pi, Hn,i(xn,i) = 1 2 ρa2 01,i Vp,i x2 n,i, dn,i = 0, Rn,i =rpl,i ρ, Jn,i = 0, Bn,i = 1, un,i =qf−qm. (11) Secondly, the volumetric flow rate through the throttle can be modeled analogously to the pH pipe model (5) with xe,m =ρLo,m Ao,m qm, He,m(xt,m) = 1 2 Ao,m ρLo,m x2 e,m, de,m = 0, Re,m =Je,m = 0, Be,m = 1, ue,m =pi−pr.(12) Similarly, the volumetric flow rate dynamics through the compressor towards the plenum can be modeled analogously to the pH pipe model (5) with xe,f =ρLc,f A1,f qf, He,f (xe,f ) = A1,f 2ρLc,f x2 e,f , de,f = 0, (13) Re,f =Je,f = 0, Be,f =1 1, ue,f =pl−pi ∆pf. Remark 1: Following [21] one can also impose closed coupled valve control to replace the pressure difference ∆pf in (10) by a desired dissipating term. Additionally, the model could be augmented by incorporating the compressor shaft’s rotational speed ωas a state variable, with compressor torque as the control input. F. Sector coupling components In this section, we present lumped pH electrolyzer and fuel cell models that facilitate the integration of the hydrogen system with the electrical domain that can be used for control from a macroscopic grid perspective of an integrated hydrogen network. Recall from Section III-A that the interconnection topology between the hydrogen system and the electric power system is modeled by the graph GP= (NF∪NP,EP). First, we describe the electrolyzer model. Due to the increased integration and the increase in the dimensioning of the electrolyzers, it becomes crucial to model the electric dynamic behavior of the electrolyzer [22]. We focus on polymer exchange membrane (PEM) electrolyzers and follow the modeling in [23] to obtain a model that provides a dynamic description of the voltage and current dynamics as well as the outgoing hydrogen volumetric flow rates. The electrical dynamics of the electrolyzer represented by the node nP,i ∈ NPis essentially described by the activation overpotential va,i, which is a key component of the electrolyzer voltage that is mainly influenced by the reaction kinetics of the electrochemical process. It reflects the additional energy required to overcome the activation energy barrier for the hydrogen and oxygen evolution reactions. The activation overpotential va,i can be modeled by [23] Ca,i ˙ va,i =ii−R−1 a,i va,i,(14) where ii∈Ris the total electrolyzer current, Ra,i >0is the activation resistance, Ca,i =CDL,cell,iAi nc,i >0is the total double-layer capacitance, based on the cell surface area Ai and CDL,cell,i is the cell capacitance and nc,i is the number of cells. The total electrolyzer voltage viat the electrical node nP,i ∈ NPis given as [23] vi=voc,i +va,i +voh,i,(15) with voc,i being the open-circuit voltage of the electrolyzer, as described by the Nernst equation, i.e., voc,i =nc,i vst −βT(Ti−Tst) + RTi 2Fln pH2 √pO2pH2O, where βT= 0.0009 V/K is the temperature coefficient accounting for the variation of the open-circuit voltage with temperature Ti>0and F≈9.6485 ·104C/mol is the Faraday constant. The cell voltage under standard conditions vst is typically 1.23 V. Standard conditions are defined as a temperature of Tst = 298.15 K and a pressure of Pst = 101325 Pa, and pH2, pO2, pH2Oare the partial pressures of hydrogen, oxygen, and water [24]. The Ohmic overpotential voh,i accounts for resistive losses in the electrolyte and electrodes and is given by voh,i =nc,i δm,i σm,iAi ii,(16) with membrane thickness and conductivity δm,i, σm,i >0. The fuel cell can be thought of as the counterpart to an electrolyzer, i.e. hydrogen is consumed to produce electrical power. Thus, the dynamics of the activation overpotential va,i at the fuel cell node nP,i ∈ NPin terms of the fuel cell current iiis identical to the electrolyzer behavior (14) and the total fuel cell voltage viis given by [25] −vi=−voc,i +va,i +voh,i,(17) with its output power depending on the hydrogen input flow rate. The volumetric flow rate of hydrogen can be expressed in terms of the electrolyzer and fuel cell currents iias [26] qsc,i =MH2 zρF ii,(18) where MH2= 2.016 g/mol is the molar mass of hydrogen, z= 2 is the number of electrons transferred per hydrogen molecule. Assumption 2: Consider the electrolyzer and fuel cell dynamics (14), (15), and (17). (i) The open circuit voltage voc,i is constant for all i= 1, . . . , E +L; (ii) The activation resistance Ra,i is constant for all i= 1, . . . , E +L. Assumption 2 is valid for moderate temperature changes within the electrolyzer and the fuel cell stacks, i.e., if the electrolyzers and fuel cells are operated relatively close to some desired operating state for the hydrogen volumetric flow rates and the temperature dynamics. This is reasonable in many settings, since often additional temperature control is applied to achieve a specific temperature set-point [23], [27]. For larger temperature changes, the irreversible pH framework [28] could be employed to derive a model that captures the associated thermodynamic more accurately. Furthermore, a suitable constant value for the activation resistance Ra,i can be identified in experimental setups from measurement data [29]. With Assumption 2, we can rewrite the electrolyzer dynamics (14) and (15) and the fuel cell dynamics (17), (18) as the following pH system with feed-through ˙xsc,i = (Jsc,i −Rsc,i)∇Hsc,i(xsc,i)+Bsc,iusc,i,(19) ysc,i =B⊤ sc,i∇Hsc,i(xsc,i)+Dsc,iusc,i +dsc,i, Hsc,i(xsc,i) = 1 2C−1 a,i x2 sc,i, xsc,i =Ca,iva,i, usc,i =qsc,i, dsc,i =(zρF MH2 voc,i if 1≤i≤E, −zρF MH2 voc,i if E+ 1 ≤i≤E+L, Rsc,i =R−1 a,i , Jsc,i = 0, Bsc,i =zρF MH2 , Dsc,i =nc,i δm,i σm,iAi z2ρ2F2 M2 H2 . Note that the product of the electrolyzer input and output is equal to the electrical power ysc,iusc,i =zρF MH2 viqsc,i =iivi, for i= 1,...,E, whereas the output of the fuel cell model is ysc,i =−zρF MH2 vifor all i=E+ 1, . . . , E +Lreflecting the inverse relationship. IV. PORT-HAMILTONIAN MODELING OF INTEGRATED HYDROGEN SYSTEMS In this section, we interconnect the pH component models from Section III to obtain a pH model of the overall hydrogen grid. An example of such a grid is shown in Fig. 1. The interconnected hydrogen grid without the sector coupling components is given by a combination of the pH models for the volumetric flow rate dynamics for qthat are given by (5), (12), (13) and the pH models for the pressure dynamics for pthat is given by (7), (8), (11). Thus, we define the state vector assigned to nodes and representing (weighted) pressures at storage units, junctions and compressors as xn= ((xn,i)S i=1,(xn,i)S+F i=S+1,(xn,i)N i=S+F+1). Likewise, the state vector collecting the volumetric flow rates at the different pipes and compressors is defined as xe= ((xe,i)M−2C i=1 ,(xe,i)M i=M−2C+1), where without loss of generality we group the edges ef∈ E and em∈ E corresponding to any compressor ni∈ NC, such that m=f+ 1.Then, the overall system state vector is given by xg= (xn, xe). The pressure dynamics of all components of the hydrogen grid can be written as a diagonal combination of all subsystems (7), (8) and (11), i.e., ˙xn= (Jn−Rn)∇Hn(xn)+Bnun,(20) yn=B⊤ n∇Hn(xn) = y1 n y2 n=p (pi)S+F i=1 , Jn= diag(Jn,i)N i=1 = 0N×N, Rn= diag(Rn,i)N i=1, Bn=B1 nB2 n=INIS+F 0C×(S+F), un= (u1 n,qex) := ((u1 n,i)N i=1,(qex,i)S+F i=1 ). Furthermore, we rearrange the entries of the input vector such that u1 ncollects all first entries of un,i for all i= 1, . . . , N, which are precisely the volumetric flow differences and the Hamiltonian is given by Hn=PN i=1 Hn,i. Likewise, the volumetric flow rate dynamics of all components of the hydrogen grid can be written as a diagonal combination of all subsystems (5), (13) and (12), i.e., ˙xe= (Je−Re(xe))∇He(xe)+Beue+de,(21) ye=B⊤ e∇He(xe) = y1 e y2 e=q (qM+2(i−C)−1)C i=1, Je=diag(Je,i)M i=1 =0M×M, Re(xe)=diag(Re,i(xe,i))M i=1, Be=B1 eB2 e=IM0(M−2C)×C diag(1 0⊤)C i=1, ue= (u1 e,∆p) := ((u1 e,i)M i=1,(∆pi)C i=1), and the Hamiltonian is given by He=PM i=1 He,i. The interconnection of the node dynamics (20) and edge dynamics (21) is established via u1 n=BGq=BGy1 e, u1 e=−B⊤ Gp=−B⊤ Gy1 n, which leads to the following pH model for the interconnected hydrogen grid ˙xg= (Jg−Rg(xg))∇Hg(xg)+Bgug+dg,(22) yg=B⊤ g∇Hg(xg) = (pi)S+F i=1 (qM+2(i−C)−1)C i=1, ug=qex ∆p, Hg(xg) = N X i=1 Hn,i(xn) + M X i=1 He,i(xe),∇Hg(xg)=p q, Rg(xg) = Rn0N×M 0M×NRe(xe), Jg=0N×NBG −B⊤ G0M×M, Bg=B2 n0N×C 0M×(S+F)B2 e, dg=0N (de,j)M j=1. In the following, we add E≥0electrolyzers which are connected to the storage units at n1,...,nE∈ NS and L≥0fuel cells connected to the storage units at nE+1,...,nE+L∈ NSby setting qsc = (qsc,i)E+L i=1 = (qex,i)E+L i=1 . For the resulting system, we keep the hydrogen volumetric flow rates qsc as input variables. This leads to the following sector-coupled hydrogen system with state vector xcg = (xg,(xsc,i)E+L i=1 ) ˙xcg = (Jcg −Rcg(xcg))∇Hcg(xcg)+Bcgucg +dcg,1, ycg =B⊤ cg∇Hcg(xcg)+Dcgucg +dcg,2,(23) Rcg(xcg) = Rg(xg) 0(N+M)×(E+L) 0(E+L)×(N+M)Rsc , Jcg =Jg0(N+M)×(E+L) 0(E+L)×(N+M)0(E+L)×(E+L), ucg =  qsc (qex,i)S+F i=E+L+1 (∆pi)C i=1  , Bcg ="Bg hzρF MH2 IE+L0(E+L)×(N−E−L)i#, dcg,1= (dg,0E+L), dcg,2= ((dsc,i)E+L i=1 ,0N−E−L), Dcg =diag(Dsc,i)E+L i=1 0(E+L)×(N−E−L) 0(N−E−L)×(E+L)0(N−E−L)×(N−E−L), and with the overall Hamiltonian Hcg(xcg) = Hg(xg) + E+L X i=1 Hsc,i(xsc,i). V. CONCLUSION In this work, we have presented a unified systematic port-Hamiltonian (pH) framework for modeling hydrogen systems, including storage units, electrolyzers, fuel cells, compressors, and pipelines. The inherent passivity properties of pH systems open the door to a structured analysis and controller synthesis for this emerging class of new energy systems. In future work, we will focus on further refining the presented modeling approach, especially by relaxing some of the assumptions made and complementing the hydrogen system with an explicit pH representation of the electrical power system. REFERENCES [1] F. Neumann, E. Zeyen, M. Victoria, and T. Brown, “The potential role of a hydrogen network in Europe,” Joule, vol. 7, no. 8, pp. 1793–1817, 2023. [2] M. Yue, H. Lambert, E. Pahon, R. Roche, S. Jemei, and D. Hissel, “Hydrogen energy systems: A critical review of technologies, applications, trends and challenges,” Renewable and Sustainable Energy Reviews, vol. 146, p. 111180, 2021. [3] Vereinigung der Fernleitungsnetzbetreiber Gas e.V., “Joint application for the hydrogen core network,” July 2024. [4] F. Dawood, M. Anda, and G. Shafiullah, “Hydrogen production for energy: An overview,” International Journal of Hydrogen Energy, vol. 45, no. 7, pp. 3847–3869, 2020. [5] A. Schaft and D. Jeltsema, “Port-Hamiltonian systems theory: An introductory overview,” Found. Trends Syst. Control., 2014. [6] S. Fiaz, D. Zonetti, R. Ortega, J. M. A. Scherpen, and A. J. van der Schaft, “A port-Hamiltonian approach to power network modeling and analysis,” Eur. J. Control, vol. 19, no. 6, pp. 477–485, 2013. [7] J. Schiffer, R. Ortega, A. Astolfi, J. Raisch, and T. Sezi, “Conditions for stability of droop-controlled inverter-based microgrids,” Automatica, vol. 50, no. 10, pp. 2457–2469, 2014. [8] H. Gernandt, B. Severino, X. Zhang, V. Mehrmann, and K. Strunz, “Port-Hamiltonian modeling and control of electric vehicle charging stations,” IEEE Transactions on Transportation Electrification, vol. 11, no. 1, pp. 2897–2907, 2025. [9] S.-A. Hauschild, N. Marheineke, V. Mehrmann, J. Mohring, A. M. Badlyan, M. Rein, and M. Schmidt, “Port-Hamiltonian modeling of district heating networks,” in Progress in differential-algebraic equations II, pp. 333–355, Springer, 2020. [10] F. Strehle, J. E. Machado, M. Cucuzzella, A. J. Malan, J. M. Scherpen, and S. Hohmann, “Port-Hamiltonian modeling of hydraulics in 4th generation district heating networks,” in 2022 IEEE 61st Conference on Decision and Control (CDC), pp. 1182–1189, IEEE, 2022. [11] A. Krishna and J. Schiffer, “A port-Hamiltonian approach to modeling and control of an electro-thermal microgrid,” IFAC-PapersOnLine, vol. 54, no. 19, pp. 287–293, 2021. [12] P. Domschke, J. Giesselmann, J. Lang, T. Breiten, V. Mehrmann, R. Morandin, B. Hiller, and C. Tischendorf, “Gas network modeling: An overview (extended english version),” February 2023. Preprint. [13] A. J. Malan, L. Rausche, F. Strehle, and S. Hohmann, “PortHamiltonian modelling for analysis and control of gas networks,” IFAC-PapersOnLine, vol. 56, no. 2, pp. 5431–5437, 2023. [14] A. J. Malan, A. Gießler, F. Strehle, and S. Hohmann, “Passivity-based pressure control for grid-forming compressors in gas networks,” in 2024 European Control Conference (ECC), pp. 1097–1104, 2024. [15] T. Bendokat, P. Benner, S. Grundel, and A. S. Nayak, “Modelling gas networks with compressors: A port-Hamiltonian approach,” PAMM, vol. 24, no. 4, p. e202400164, 2024. [16] L. Kumar, J. Chen, C. Wu, Y. Chen, and A. van der Schaft, “A segmented model based fuel delivery control of PEM fuel cells: A port-Hamiltonian approach,” Automatica, vol. 168, p. 111814, 2024. [17] D. Sbarbaro, “On the port-Hamiltonian models of some electrochemical processes,” IFAC-PapersOnLine, vol. 51, no. 3, pp. 38–43, 2018. [18] A. van der Schaft, L2-Gain and Passivity Techniques in Nonlinear Control. Springer, 2017. [19] K. Pambour, R. Bolado-Lavin, and G. Dijkema, “An integrated transient model for simulating the operation of natural gas transport systems,” Journal of Natural Gas Science and Engineering, vol. 28, pp. 672–690, 2016. [20] N. Klopˇ ciˇ c, K. Esser, J. F. Rauh, M. Sartory, and A. Trattner, “Modelling hydrogen storage and filling systems: A dynamic and customizable toolkit,” International Journal of Hydrogen Energy, vol. 49, pp. 1180–1195, 2024. [21] J. T. Gravdahl and O. Egeland, “Passivity based compressor surge control using a close-coupled valve,” Proceedings of the 1997 COSY Workshop on Control of Nonlinear and Uncertain Systems. Zurich, Switzerland, pp. 139–143, 1997. [22] M. Espinosa-L´ opez, C. Darras, P. Poggi, R. Glises, P. Baucour, A. Rakotondrainibe, S. Besse, and P. Serre-Combe, “Modelling and experimental validation of a 46 kW PEM high pressure water electrolyzer,” Renewable Energy, vol. 119, pp. 160–173, 2018. [23] G. Lichtenberg, G. Pangalos, C. C. Y´ a˜ nez, A. Luxa, N. J¨ ores, L. Schnelle, and C. Kaufmann, “Implicit multilinear modeling: An introduction with application to energy systems,” at - Automatisierungstechnik, vol. 70, no. 1, pp. 13–30, 2022. [24] M. Espinosa-L´ opez, C. Darras, P. Poggi, R. Glises, P. Baucour, A. Rakotondrainibe, S. Besse, and P. Serre-Combe, “Modelling and experimental validation of a 46 kW PEM high pressure water electrolyzer,” Renewable Energy, vol. 119, pp. 160–173, 2018. [25] J. Pukrushpan, H. Peng, and A. Stefanopoulou, “Simulation and analysis of transient fuel cell system performance based on a dynamic reactant flow model,” ASME Int. Mech. Eng. Congress Expo., vol. 2, 01 2002. [26] M. Carmo and D. Stolten, “Chapter 4 - energy storage using hydrogen produced from excess renewable electricity: Power to hydrogen,” in Science and Engineering of Hydrogen-Based Energy Technologies (P. E. V. de Miranda, ed.), pp. 165–199, Academic Press, 2019. [27] M. Pfennig, B. Schiffer, and T. Clees, “Thermodynamical and electrochemical model of a PEM electrolyzer plant in the megawatt range with a literature analysis of the fitting parameters,” International Journal of Hydrogen Energy, vol. 104, pp. 567–583, 2025. [28] H. Ramirez and Y. Le Gorrec, “An overview on irreversible portHamiltonian systems,” Entropy, vol. 24, no. 10, 2022. [29] B. Yodwong, D. Guilbert, M. Hinaje, M. Phattanasak, W. Kaewmanee, and G. Vitale, “Proton exchange membrane electrolyzer emulator for power electronics testing applications,” Processes, vol. 9, no. 3, 2021.