A Unifying Passivity-Based Framework for Pressure and Volume Flow Rate Control in District Heating Networks

A fundamental precondition for the operation of district heating networks (DHNs) is a stable hydraulic behavior. However, the ongoing transition toward a sustainable heat supply, especially the rising integration of distributed heat sources and the increasingly meshed topologies, introduces complex and potentially destabilizing hydraulic dynamics. In this work, we propose a unifying, equilibrium-independent passivity (EIP)-based control framework, which guarantees asymptotic stability of any feasible, hydraulic DHN equilibrium for a wide range of DHN setups covering different DHN generations, meshed, time-varying topologies, and multiple, dynamically interacting distributed heat sources. The obtained results hold for the state of the art as well as future DHN generations featuring, for example, multiple distributed heat sources, asymmetric pipe networks, and multiple temperature layers.


I. INTRODUCTION
District heating networks (DHNs) are a key element for a holistic energy transition, particularly in densely populated areas [1], [2], [3], [4]. For their operation, well-defined and stable hydraulic conditions are a fundamental requirement as they form the basis for the actual thermal power flows [3], [5].
In the operation of traditional 2nd or 3rd generation DHNs, the hydraulics and thus the thermal power flows are wellunderstood (see e.g., [6, pp. 52-54]). However, emerging 4th generation DHNs bring about challenges that call for new strategies and methods of operating, controlling and analyzing DHNs [1], [2], [3], [7]. 1 Most prominently, we can observe a decentralization trend with several interacting subsystems and controllable components (pumps, control valves). Primarily, this is due to the integration of ever more renewable and michele.cucuzzella@unipv.it distributed heat generation units (DGUs) such as heat pumps, waste heat recycling, combined heat and power, waste/biomassto-energy, solar, and geothermal heat plants, and [1], [2]. Additionally, distributed variable-speed pumps (DVSPs) installed at every DGU and some (up to all) consumer substations, have shown considerable potential to reduce the overall electrical energy required to operate DHN pumps [4], [8], [9].
In summary, from a control point of view, the decentralization trend with the increasingly actuated subsystems, the naturally intermittent supply behavior of DGUs [2], [14], and the new operation strategies lead to more frequent topology changes and plug-in and plug-out operations of units. Additionally, complex pressure and volume flow dynamics and interactions on small time scales are introduced [3]. For example, it is possible to observe more frequent volume flow reversals in pipes [15], [16] and interactions between the pump and valve controllers, which may lead to severe hydraulic oscillations [4], [8], [17].
Upon closer examination, the discussion above reveals parallels to the trends and challenges in modern power systems such as microgrids (see for example [18] or the discussion in [7]). To provide stabilization in the face of intermittent generation and distributed units, the power system control community is looking for decentralized control solutions. Thus, the similar trends and challenges motivate to explore decentralized control designs for novel pressure and volume flow rate controllers. Such decentralized controllers provide scalability and allow for an easy addition or removal of subsystems in a plug-andplay fashion without requiring communication, without adapt-ing the other controllers, and without endangering pressure and volume flow rate stability. Recent works on microgrids (see e.g., [19], [20]) have shown that the compositional properties of passive systems present a promising, unifying framework for designing decentralized controllers and realizing modular energy systems in which topologies change frequently and different controllers interact dynamically via the network.
In the literature, the field of decentralized, passivity-based hydraulic control of DHNs has recently been explored in [21]. However, the considered DHN model exhibits a number of restrictions such as symmetric DHN topologies 2 , two temperature layers only, static pump models, no pressure holding units, and valves modelled as non-actuated components. Furthermore, pumps are assumed to be installed at every producer and every consumer, which excludes traditional DHNs in which consumers regulate their volume flow rates only via control valves. The same restrictions underlie the DHN models used in [13], [22], [23], [24], where, additionally, only single producer DHNs are considered. Further noteworthy works that contribute towards a passivity-based control design and analysis of DHNs by introducing port-Hamiltonian system (PHS) models are [25], [26]; see also [27], where water distribution networks are modeled as PHSs. However, [25] exhibits the same setup restrictions as [13], [22], [23], [24], while [26] overcomes most of these restrictions with the exception of modeling valves as non-actuated components and considering pumps installed at every producer and consumer.

A. Contributions
In this work, we present a unifying framework based on equilibrium-independent passivity (EIP) which allows for decentralized pressure and volume flow rate control in DHNs with meshed, time-varying topologies and different, dynamically interacting DGUs.
Firstly, we formalize the considered DHN setups as digraphs. Our considerations overcome the restrictions in the above-referenced literature and are based on a comprehensive dynamic hydraulic DHN model that encompasses traditional 2nd and 3rd generation DHNs, future 4th generation DHNs as well as intermediate stages. That is, we allow for an arbitrary number of DGUs and consumers to be connected in a meshed DHN topology with multiple temperature layers (temperature cascading), DVSP configurations, and optional booster pumps at consumers and pipes. Moreover, we incorporate dynamic models for the pumps, consider valves as actuators, and explicitly account for the presence of pressure holding units.
Secondly, we derive explicit input-state-output port-Hamiltonian system (ISO-PHS) models for each of the DHN subsystems. Although not strictly necessary for the subsequent control design and modular stability analysis, the systematic PHS modeling is nevertheless useful during these steps. In particular, it gives a clear perspective on which input-output ports are accessible for control and over which ports subsystems 2 In symmetric DHN topologies, supply and return pipes are laid in parallel. This excludes practically relevant cases with meshed supply pipe networks and tree-like return pipe networks or more complex structures arising in multilayer topologies. interact with each other. Furthermore, the passivity properties with respect to these ports and the Hamiltonian as storage function candidate are directly visible.
Thirdly, we propose pressure controllers for pumps, based on an algebraic interconnection and damping assignment (IDA) [28] extended with integral action on the non-passive output [29]; volume flow rate controllers for pumps, based on a state-feedback with integral action on the passive output; and volume flow rate controllers for valves, based on a proportional-integral (PI) control of a modified passive output inspired by [30].
Fourthly, by leveraging the EIP properties of the subsystems, the skew-symmetric nature of their interconnection structure, and LaSalle's invariance principle, we prove asymptotic stability of any forced hydraulic DHN equilibrium in a modular, bottom-up manner. This establishes a unifying framework for pressure and volume flow rate control, where multiple DHN subsystems, if EIP, can enter or leave the DHN in a plugand-play fashion without having any impact on the stability properties of the hydraulic equilibrium. Simulations based on realistic DHN data validate our theoretical results.

II. SYSTEM SETUP
In this section, we outline the general DHN setups that are covered by our approach and formally describe them as weakly connected digraphs. We allow for DHNs with asymmetric, meshed topologies, multiple temperature layers, multiple DGUs, pressure holding units, booster pumps at consumers and pipes, DVSP configurations as well as mixing connections used to increase the water temperature in topologies with multiple temperature layers. To the best of our knowledge, such a comprehensive collection of DHN setups covering 2nd, 3rd, 4th generation DHNs as well as intermediate development stages has not been considered before in any control-oriented literature.

A. DHN Setup and Digraph Representation
We describe a DHN as a weakly connected digraph G = (N , E) without self-loops as shown in Fig. 1. The edges E are partitioned into four sets: D = {1, . . . , D}, D ≥ 1, represents the DGUs, L = {D + 1, · · · , D + L}, L ≥ 1, the consumers (loads), P = {D + L + 1, · · · , D + L + P }, P ≥ 2, the pipes, and M = {D + L + P + 1, · · · , D + L + P + M }, M ≥ 0, the mixing connections. Conventionally, the nodes N correspond to ideal system junctions interconnecting DGUs, consumers, and mixing connections with the pipe network of the DHN. At an ideal junction, all volume flow rates sum up to zero. However, in this work, we also view pressure holding units and elasticity capacitors arising from equipment in the DGU and consumer circuits as nodes. Therefore, we assume that N = H ∪ C ∪ K, where H are pressure holding units, C are elasticity capacitors, and K are the remaining ideal junctions.
The orientation of the edges represents the arbitrary reference direction of positive flows. Moreover, for any i ∈ E, N − i and N + i denote its source and target node, respectively. Analogously, for a given node j ∈ N , E − j and E + j denote Illustration of a traditional 2nd or 3rd generation DHN with high temperature supply and medium temperature return (left) and a 4th generation DHN with medium temperature supply and low temperature return (right); between the temperature layers there may be an any number of D ≥ 1 DGU and L ≥ 1 consumer edges; the temperature layers coincide with the two hydraulic layers (dashed grey bubbles).
the sets of edges with j as source node and j as target node, respectively. 3

B. Temperature and Hydraulic Layers
As illustrated in Figs. 1-3, a DHN may comprise different temperature layers. We distinguish between three temperature layers, i.e., high temperature (80°C-120°C), medium temperature (40°C-70°C), and low temperature (20°C-40°C) [1], [11], [12], [3] [31, pp. 16-17] [6, p. 44]. The high and medium temperature layers form the supply and return of the dominating 2nd and 3rd generation DHNs, while the medium 3 In practice DGUs and consumers are always connected via pipes and never directly connect to the same node (see Fig. 1). Figure 3. Illustration of a DHN with the three temperature layers high (red), medium (orange), low (blue), and the two hydraulic layers (dashed grey bubbles); between the temperature layers there may be any number of D ≥ 1 DGU and L ≥ 1 consumer edges; between the high and medium layer may be any number of M ≥ 0 mixing connections. and low temperature layers form the supply and return of the emerging 4th generation DHNs (see Fig. 2) [1] [31, pp. 16-17].
In future DHNs, the medium temperature return of a 2nd or 3rd generation DHN may additionally serve as supply for (new) low-temperature DHN sections, yielding a three temperature layer topology as illustrated in Figs. 1 and 3. Such low-temperature DHN sections allow to efficiently use the heat energy in a DHN (temperature cascading) and integrate renewable heat sources (e.g., waste heat, solar thermal, heat pumps) and new consumers (e.g., low-energy buildings) into existing DHNs [6, p. 44] [10], [11], [12]. However, due to the ongoing trend of decreasing DHN temperatures, particularly in summer, the medium temperature might not be sufficiently high to cover the heat demand of some low-temperature consumers. Thus, low-temperature DHN sections in a three layer topology typically have at least one mixing connection, i.e., an edge i ∈ M, that allows to boost the temperature by mixing high with medium temperature water (see, e.g., node 7 in Fig. 1) [12], [10], [11].
Furthermore, in a three temperature layer topology, the low temperature water is typically fed directly into the medium temperature layer (see node 5 in Fig. 1 and the set of pipe edges P between low and medium temperature in Fig. 3) [10], [11], [12]. This implies that despite there possibly being three temperature layers, there are exactly two hydraulic layers (see Fig. 3). The number of hydraulic layers can be defined as follows: Definition 1: A DHN has n l ≥ 2 hydraulic layers, where n l is the number of weakly connected subgraphs G 1 , . . . , G n l obtained by removing all edges i ∈ D ∪L∪M, i.e., all DGUs, consumers, and mixing connections, from G.
Remark 1: We want to emphasize that our setup does not require to have three temperature layers. Traditional DHNs with high temperature supply and medium temperature return are covered as well as standalone 4th generation DHNs with medium temperature supply and low temperature return. In any case, there are always exactly two hydraulic layers containing all pipe edges.

III. HYDRAULIC MODELING
With the DHN setup formalized as a digraph, we now present the mathematical models describing the hydraulic dynamics of the edges and nodes. In Sections III-A1 and III-A2, we first model the main actuators responsible for pressure and volume flow rate control, i.e., pumps and valves. They serve as building blocks for the hydraulic models of the DGU, consumer, pipe, and mixing connection edges in Sections III-B2-III-B4, as well as the pressure holding nodes in Section III-C1. Lastly, the capacitive, and simple junction nodes are modeled in sections III-C2 and III-C3. For the modeling, we make the following assumptions which are valid under normal operating conditions (see also [13], [32]): Assumption 1: The compressibility of water is neglected. Any reference and nominal pressure values as well as all model parameters are strictly positive. Pressure losses inside pipes λ(q) : R → R and valves µ(q, s) : R × R ≥0 → R caused by volume flow rates q ∈ R are continuously differentiable Figure 4. Equivalent circuit of a linear, second-order approximation of pump dynamics (cf. [33]).
As a starting point for such an improved design and analysis, we follow [26], [33] and model each pump by a linear equivalent RLC circuit as shown in Fig. 4. The RLC circuit arises by approximating the complex arrangement of power electronics, speed-controlled AC motor, and pump hydraulics by a linear second-order system. By applying Kirchhoff's voltage law (KVL) and KCL to Fig. 4, we obtain the model d dt where x i is the state vector, u i the control input, and y i the measurable output vector. The additional input d i and output z i model the interaction or physical interconnection between the pump and other subsystems, e.g., a DGU. Furthermore, R P,i , J P,i , C P,i are the model parameters, p P,i is the pressure difference produced by the pump between its terminals, q i is the volume flow rate through the pump, and q P,i is an auxiliary variable without physical interpretation. The control input u P,i can be interpreted as a pressure source originating from the rotational speed of the pump produced by an AC motor [33]. 4 Remark 2: In (1), R P,i , J P,i , and C P,i are black box parameters without physical interpretation. This is due to the fact that in the RLC circuit representation, the speed control and AC motor dynamics are merged with the hydraulic pump dynamics comprising fluid mass inertia, pressure losses, and hydraulic capacities due to fluid compressibility and fluid volume. In practice, R P,i , J P,i , and C P,i can be identified from measurement data obtained by operating the respective pump in typical scenarios (see e.g., [33]). Alternatively, they can be fitted in simulations to match characteristic curves provided in data sheets.
2) Control valves: Besides pumps, valves are the main actuators in DHNs. Their main task is the regulation of volume flow rates [6, pp. 143-145,151] [31, p. 19,29][2] [37]. In order to establish a desired volume flow rate, valves adjust their pressure drop by varying their stem position between fully closed (s v,i = 0) and fully open (s v,i = 1). Thus, they behave as variable, nonlinear flow resistors. In order to avoid volume surge behavior around their closing point, valves are typically designed such that the stem has a lower limit just above zero in normal operation [6, pp. 145]. Consequently, the following assumption can be made: Assumption 2: Valves are designed appropriately such that in normal operation the stem position never reaches the zero value, i.e., The nonlinear characteristic pressure drop equation of any valve is given by [37,Eq. (18)], [9,Eq. (5)] where s v,i ∈ s min v,i , 1 , 0 < s min v,i ≤ ǫ v,i is the stem position, q i ∈ R the volume flow rate through the valve, C v,i > 0 the flow capacity of the valve, and in (2a), the pressure drop can be written as where u v,i (s i ) : s min v,i , 1 → R ≥0 is a bijective mapping of the actual stem position s v,i to the virtual control input u v,i , andμ i (q i ) a continuously differentiable, strictly monotonically increasing function satisfyingμ i (0) = 0. Note that (3) is affine in the virtual control input u v,i . Figure 5. Equivalent circuit of a hydraulic DGU model (i ∈ D) with pressure holding (black voltage source) and circulation circuit (red) [26, Fig. 2]; without loss of generality, the capacitance C j may be lumped with the pressure holding (see Section III-C1); j ∈ N − i and k ∈ N + i . Figure 6. Equivalent circuit of a hydraulic pipe model (i ∈ P) with optional booster pump; j ∈ N − i and k ∈ N + i .

Remark 3: It might not be intuitively clear that
However, when considering that valve characteristics are specified assuming a fixed pressure drop

Remark 4:
Since valves can only function as variable, nonlinear flow resistors, valve-based volume flow rate control requires sufficient differential pressure over the hydraulic DGU or consumer circuit the valve is part of. In future DHNs with frequently changing hydraulic conditions, this motivates the addition of booster pumps in some consumer circuits or pipes (see Sections III-B3, III-B1 and Assumption 3).
Remark 5: Assumption 2 implies that u v ∈ 1, u max v,i . However, for the sake of simplicity and in line with classical feedback control design, we do not consider the possibility of control input saturation explicitly during the control design stage. Instead, we discuss the well-behaved nature of input saturation by means of our numerical simulations in Section VII.
B. Edge Dynamics 1) Pipes: The hydraulic pipe model at an edge i ∈ P is illustrated in the equivalent circuit diagram in Fig. 6. Following the literature (see, e.g., [7], [13], [35], [23], [26]), we model the pipe friction by a nonlinear, volume flow-dependent resistance λ i (q i ) (see Assumption 1) and the volume inertia by the linear inductance J i . In contrast to prior works, we assume that some pipes might have booster pumps in series. Such pumps are represented by the voltage source in Fig. 6 and modeled by (1). By applying KVL and KCL to Fig. 6, we obtain the model for each i ∈ P as where x i is the state vector, u i the control input, y i the measurable output vector, (d i , z i ) the interaction (coupling) port pair, and (j, k) ∈ N − i × N + i are the source and target nodes of i with pressures p j and p k , respectively.
Remark 6: Any pipe i ∈ P without a booster pump can be modeled by (4) by fixing u P,i = p P,i = 0 and removing the part corresponding to the dynamics of q P,i and p P,i .
2) DGUs: From a hydraulic viewpoint, a DGU may comprise two main parts as illustrated in Fig. 5: a circulation circuit (red) (see [39], [40]) and an optional pressure holding unit (see [6, pp. 54-55]). We view the circulation circuit (in red) as the actual edge i ∈ D, which comprises a circulation pump and a control valve connected in series with pipes and a heat exchanger. All pipes are lumped into the nonlinear, volume flow-dependent resistance λ i (q i ) and the inductance J i , which represent the pipe friction and volume inertia, respectively. The control valve is modeled as a variable, nonlinear resistancê (2). The circulation pump is modeled by (1). By applying KVL and KCL to the red part in Fig. 5, we obtain the model for each where x i is the state vector, u i the control input vector, y i the measurable output vector, (d i , z i ) the interaction (coupling) port pair, and (j, k) ∈ N − i × N + i are the source and target nodes of i with pressures p j and p k , respectively.
Remark 7: The pressure holding unit of a given DGU is represented by the (black) voltage source p P,h shown in Fig. 5. The capacitances C j and C k model the hydraulic elasticity of all the components in the DGU circulation circuit, particularly of the heat exchanger (see [41], [38]). For simplicity, we assume that C j is lumped with C P,h of the pressure holding (see Fig. 4). Furthermore, to clearly describe the network interconnection among all the subsystems in the DHN, we view the pressure holding and the elasticity capacitances as nodes of the DHN graph; we elaborate on their models in Section III-C.

3) Consumers:
Most of nowadays consumers are connected indirectly to a DHN via heat exchangers in series with pipes and a control valve for volume flow rate control [6, pp. 87,143]. In future DHNs, however, additional pumps are expected to be included in some (up to all) consumer circuits: either for pressure boosting to ensure a proper functioning of the control valves under unclear and changing hydraulic conditions [1], [13] or in DVSP DHN configurations [4], [8], [9]. Consequently, the hydraulic consumer circuit at an edge i ∈ L is modeled similarly to the hydraulic DGU circulation circuit (red part in Fig. 7, [7, Fig. 2]). The only differences are on the working direction of the pump and the sign convention of the volume flow rate, which is reflected in the edge orientation in the DHN digraph (see Fig. 1 and [6, pp. 87,143]). Furthermore, pressure holding units are typically not installed at consumers. By applying KVL and KCL to Fig. 7, we obtain the model for each i ∈ L as in (5). Remark 8: Any consumer i ∈ L without a pump can be modeled by (5) by fixing u P,i = p P,i = 0 and removing the part corresponding to the dynamics of q P,i and p P,i . Such consumers regulate their flow rate through their respective control valve.

4)
Mixing connection: As outlined in Section II-B, future DHNs may have a topology with three temperature layers. In order to guarantee sufficient heat supply of the low-temperature sections, the medium temperature water is typically mixed with high temperature water via a mixing connection before it is fed into the low-temperature section (see, e.g., node 7 in Fig. 1). The hydraulic circuit of a mixing connection at an edge i ∈ M is illustrated in Fig. 8. It comprises a pipe in series with a control valve. By applying KVL to Fig. 8, we obtain the model for each i ∈ M as where x i is the state, u i the control input, y i the measurable output, (d i , z i ) the interaction (coupling) port pair, and (j, k) ∈ N − i × N + i are the source and target nodes of i with pressures p j and p k , respectively.

C. Node Models
As outlined in Sections II-A, III-B2 and III-B3, the set of nodes N of the DHN graph G = (N , E) is the union of three disjoint sets: where H is the set of pressure holding units, C is the set of elasticity capacitances, and K is the set of simple junctions. 1) Pressure holding: Pressure holding units are realized technically in two ways: dynamic pressure holding with a pressure dictation pump, and static pressure holding with a closed vessel [6, pp. 54-56]. Furthermore, pressure holding units are almost exclusively installed on the suction side of circulation pumps (pre-pressure control) (see Fig. 5) and are instrumental in preventing cavitation [6, pp. 54-55], [17].
Dynamic pressure holding is typically conducted in larger DGUs with powerful circulation pumps. It is realized by a pressure dictation pump located between a pressurized container and the DHN [6, pp. 54-55], [42, Fig. 1]. As outlined in Section III-A1, we approximate the dynamics of any pump by the linear second-order system (1). Thus, the case in which a dynamic pressure holding unit is installed at a DGU i ∈ D is equivalent to replacing the black voltage source in Fig. 5 by the RLC circuit in Fig. 4. Note that in contrast to the circulation pump, which is coupled with the circulation circuit (red part in Fig. 5), the black voltage source already represents the entire pressure holding unit, i.e., we assume that the dictation pump is lumped together with the pressurized container. Thus, the model for each j ∈ H is similar to (1) and given by d dt J P,j q P,j C P,j p P,j xj = −p P,j − R P,j q P,j q P,j fj (xj ) where x j is the state vector, u j the control input, y j the measurable output vector, (d j , z j ) the interaction (coupling) port pair, and I j ⊆ E the set of edges that are incident to j. Remark 9: A static pressure holding is used in smaller DGUs with compact circulation pumps. It is realized by directly adding a closed, pressurized vessel. In an equivalent circuit perspective, this can be understood as a preloaded capacitor. Thus, in case of static pressure holding, we simply replace the black voltage source in Fig. 5 with a capacitor C P,h that we consider to be lumped with C j .
Remark 10: The presented DHN model allows for pressure holding units to be installed at the suction side of each DGU i ∈ D. In practice, however, it is typically sufficient to operate only one of them at any given time (see, e.g., [6, p. 55][42, Fig. 1] and Fig. 1). This is due to the fact that in a connected and closed hydraulic network, the pressure values at each node j ∈ V are uniquely determined by fixing one node pressure value and expressing all edge subsystems by means of differential pressure dynamics as done in Section III-B (see, e.g., [43,Sec. 2

.3]).
2) Capacitive nodes: Invoking the volume balance, we obtain the model for each j ∈ C as where x j is the state, (d j , z j ) the interaction (coupling) port pair, and I j the set of edges that are incident to j.
3) Simple junctions: The model for each j ∈ K is analogous to (9) with C j fixed to zero and treating z j = p j as an algebraic variable. Note that there is no state variable to describe the behavior of simple junctions.

IV. PROBLEM FORMULATION AND APPROACH
With the DHN formalized as a digraph and its hydraulic model established, we can now formulate the main pressure and volume flow rate control problems addressed in this work. As outlined in Section I, we focus on decentralized control schemes with plug-and-play capabilities. Furthermore, we keep to the terminology of electrical power systems and suppose that DGUs may operate either in grid-forming or grid- DGUs in grid-forming mode (i ∈ D form ) actively form the hydraulic conditions required to operate DHNs by regulating the differential pressure generated by their circulation pumps to desired setpoints p * P,i [6, p. 47-48] [44]. In this case, the control valves in their circulation circuits are fully open, i.e., u v,i = u v,i | sv,i=1 . 5 DGUs in grid-feeding mode (i ∈ D valve ⊆ D feed ) regulate the volume flow rate through their circulation circuits to desired setpoints q * i by means of their control valves. Under approximately constant water temperature, this is equivalent to controlling the thermal energy they feed into the DHN (see [45], [39][43, Sec. 2.3]). Note that for a proper functioning of the control valve, the circulation pump still introduces some differential pressure p * P,i , which is then throttled by the control valve such that the desired q * i is reached.
Consumers (i ∈ L) regulate the thermal energy they consume by controlling their volume flow rates to desired setpoints q * i [6, pp. 143-145,151] [31, p. 29]. Traditionally, this is done by control valves only. The set of consumers for which u v,i (s v,i ) is the main control input is thus denoted by L valve ⊆ L. However, as discussed in Section III-B3, booster pumps might be added to some consumer circuits. We identify these consumers by the set L boost ⊆ L. In each consumer i ∈ L boost , the pump pressure is controlled to some desired setpoint p * P,i , which is then throttled by the control valve such that the desired q * i is reached.
For pipes i ∈ P, we identify by P boost ⊆ P the subset of pipes that have a booster pump connected in series. These pumps are in charge of counteracting the differential pressure loss over the corresponding pipe by introducing a differential pressure p * P,i . In mixing connections i ∈ M, a desired volume flow rate q * i is stabilized such that a desired mixing ratio of high and medium temperature water is achieved (see Section III-B4).
Last but not least, we note that regardless of the operation mode of the DGUs, each pressure holding unit j ∈ H associated to a given DGU regulates the pressure p j at the suction side of the circulation pump j ∈ N − i (see Fig. 5) to a suitable setpoint p * j = p * P,h . This pressure also serves as the static pressure in a DHN (see [6, p. 55]).
Remark 11: In some DVSP setups, it is suggested to directly conduct the volume flow rate control (grid-feeding, consumer) by pumps without including any control valves in the circuits of DGUs i ∈ D VSP ⊆ D feed and consumers i ∈ L VSP ⊆ L [4], [8]. 6 In [9], a hybrid DVSP setup is suggested in which all DGUs and some consumers have only pumps, while some consumers have only control valves. Consequently, it is apparent that depending on the topology and producerconsumer configuration, different hydraulic designs of DGU and consumer circuits might be beneficial. In this work, we thus allow for all possible combinations (see Problem 1 below) and show that under the roof of the passivity-based framework presented in the sequel, asymptotic stability of the desired pressures and volume flow rates can always be guaranteed.
In summary, we note that the control tasks boil down to pressure and volume flow rate control of pumps and volume flow rate control of valves.
Problem 1: Consider a DHN as described in Section III. Design decentralized controllers for pumps and valves that asymptotically stabilize any forced, hydraulic DHN equilibrium with the following characteristics: (a) For each i ∈ D form , s v,i = 1 is fixed and u P,i is such thatp P,i = p * P,i . (b) For each i ∈ D valve , u P,i and u v,i are such thatp P,i = p * P,i andq i = q * i . (c) For each i ∈ D VSP , s v,i = 1 is fixed and u P,i is such i) For each j ∈ H, u P,h is such thatp P,j = p * P,j . Remark 12: Note that similar to the hierarchical control of power systems, the setpoints of the volume flow rates and pressures are assumed to be known and specified by a higher level control ensuring sensible DHN operation, i.e., ensuring that a feasible hydraulic equilibrium exists.

A. Approach
Inspired by [46,Theorem 3.1], [19], [7], we address Problem 1 in a modular, bottom-up manner by means of a passivitybased approach. The leading idea is that the hydraulic equilibrium is stable if two conditions are satisfied: C.1 each, possibly closed-loop subsystem at the edges and nodes of the DHN digraph is EIP with positive definite storage function; C.2 the interconnnection structure of the DHN subsystems is power-preserving.
Asymptotic stability of the hydraulic equilibrium can be investigated either by invoking LaSalle's Invariance Principle or checking for strict EIP of all subsystems.

V. PASSIVITY-BASED CONTROL DESIGN
In Section III, we have established that the dynamics of any edge or node i ∈ E ∪ N of a DHN, with the exception of simple junctions, can be written in the general forṁ Based on [26], it is possible to represent each of these subsystems as an input-state-output port-Hamiltonian system (ISO-PHS), which in general can be written as follows [36, p. 114]): is the nonlinear, vector-valued damping function. To improve readability all details are included in Appendix A.
The ISO-PHS representation gives a clear perspective on which input-output ports are accessible-or are convenient for control-and over which ports subsystems interact with each other: note that the output ζ i is not necessarily the measurable output. Furthermore, the passivity properties with respect to these ports and the Hamiltonian as storage function are directly visible [36]. Indeed, it is direct to see that along system solutions the Hamiltonian H i satisfies the following inequality:Ḣ Then, considering the approach described in Section IV-A, the main goal in this section is to design decentralized controllers of the forṁ for the pumps and valves such that each closed-loop subsystem i ∈ D ∪ L ∪ P boost ∪ M ∪ H fulfills its respective steadystate characteristic from Problem 1 (a)-(i) and is EIP with respect to its interaction port pair (d i , z i ) and some positive definite storage function. That is, it is to be shown that for any feasible equilibrium pair (x i ,d i ) withx i fulfilling the respective requirement (a)-(i) in Problem 1, there exists a scalar function x i →Ĥ i (x i ) that is positive definite with respect tox i and satisfieṡ Subsequently, in Section VI, it will be shown that the interconnection among the DHN subsystems is such that making i∈E∪NĤ i (x i ) a bona fide Lyapunov function. Even though the design that we propose in this section for u i follows passivity-based arguments, we do not restrict ourselves to use only the passive output ζ i . Moreover, the definition of the above-mentioned storage functionĤ i will be inspired by the Hamiltonian H i of the system, but will not be equal to it: note that H i is a Lyapunov function of any ISO-PHS only if u i and d i are zero at steady-state, which is not the situation we face in this paper.

A. Operating Modes of the DHN Subsystems
In a first preliminary step, we investigate what kind of degrees of freedom are available in choosing the operating modes of the DHN subsystems. In particular, we are interested in determining (i) if DGUs are subject to any restrictions on their choice of operating modes, i.e., grid-forming or grid-feeding, and (ii) whether all consumers and mixing connections can in principle, i.e., disregarding any technical constraints, independently control their volume flow rates to desired steady-state values.
Proposition 1: Consider a DHN as modeled in Section III with D ≥ 1 DGUs, L ≥ 1 consumers, P ≥ 2 pipes, M ≥ 0 mixing connections, and two hydraulic layers G 1 , G 2 . Then, the following holds: (i) there must be at least one grid-forming DGU i ∈ D form connecting the two hydraulic layers G 1 , G 2 ; (ii) each steady-state volume flow rateq i , i ∈ D feed ∪ L ∪ M ∪ P loop , is an independent variable, where each P loop ⊆ P forms an independent loop within G 1 , or G 2 .
Proof: Following [47, pp. 477-482] (see also [13], [35]), we make use of the fundamental loop analysis of circuit theory. First we note that since the digraph G representing a DHN is weakly connected, it admits a generally non-unique spanning tree T (see, e.g., [13], [37]). The spanning tree T is a weakly connected subgraph of G that contains all nodes of G and no loops [47, p. 477]. Any edge of G not in T is referred to as chord and creates a (fundamental) loop when added to T . For G = (N , E), there are |E| − |N | + 1 chords and loops, respectively. By applying KCL, it can be shown that each steady-state flow through an edge in T is the superposition of one or more of the steady-state loop flows [47, p. 482]. By setting the loop flows equal to the chord flows, we thus find that the |E| − |N | + 1 chord flows form a complete set of independent variables.
To prove (i), we note that per Definition 1, only DGU, consumer, or mixing edges may connect the two hydraulic layers G 1 and G 2 (see also Figs. 2 and 3). Additionally, the hydraulic layers may have a meshed structure. Thus, in order to create a spanning tree T of G, we have to connect G 1 and G 2 via exactly one edge in the union D ∪ L ∪ M and might have to remove some pipe edges in P loop ⊆ P that form loops within G 1 and G 2 , respectively. From the fundamental loop analysis we know that the steady-state flow through the edge connecting G 1 and G 2 is not independent. Thus, the edge must be a pressure-controlled subsystem in D ∪ L ∪ M, which can only be fulfilled by a grid-forming DGU i ∈ D form .
In turn, we can choose the set of chords such that it contains all DGU edges except one, all consumer edges, possibly some pipe edges that form loops within a meshed, hydraulic layer G 1 or G 2 , respectively, and all mixing edges. Consequently, the steady-state flows through these edges are independent variables, which proves (ii).
Remark 13: The independent steady-state pipe flowsq i , i ∈ P loop , provide additional degrees of freedom that might be used, e.g., to minimize pumping costs (see [37,Sec. 2.2]).

B. Pressure and Volume Flow Rate Control of Pumps
Instrumental to solving Problem 1 is the ability to regulate the pressure p P,i or the volume flow rate q i of a given pump in the subsystems i ∈ D ∪ L boost ∪ L VSP ∪ P boost ∪ H towards desired constant setpoints. For the control design, we initially consider the isolated pump model (1) and address pressure and volume flow rate regulation. Then, the closed-loop pump dynamics are interconnected to the respective subsystems D ∪ L boost ∪ L VSP ∪ P boost ∪ H.
1) Pressure Control of Pumps: For the pump pressure p P,i , we propose a controller for u P,i that is based on algebraic IDA [28] passivity-based control extended with integral action on the non-passive output of the pump model [29].
Proposition 2: Consider the pump model (1). Assign the control input u P,i as with pressure setpoint p * P,i > 0 and control parameters R p i , Q I,i > 0. Then, the closed-loop system is given by Moreover, (13) is EIP with supply rate (z i −z i )(d i −d i ) and positive definite storage function for any (feasible) equilibrium pair (d i ,z i ) and associated equilibrium statē Proof: See Appendix B. Remark 14: The controller (12) is decentralized as it only requires knowledge of local variables and parameters, such as q P,i , p P,i and J P,i . Observe in particular that if we choose R P i = R P,i , (12) becomes independent of q P,i . This is desired as q P,i does not represent a physical quantity. Moreover, note that (12d) can be written as which clearly illustrates its composition as the combination of a state feedback term and a PI term.
2) Volume Flow Rate Control via Pumps: Next, we address the task of achieving volume flow rate control of pipes or heat exchangers via pumps. For this, we focus on the model of a pump in series with a pipe element, which is equivalent to (4) but with d i treated as an arbitrary external input.
Proposition 3: Consider the model of any pump in series with a pipe element (see (4)). Assign u P,i as with volume flow rate setpoint q * i and control parameters Q I,i , K P,i satisfying Q I,i > 0, 0 < Q I,i (K P,i + 1) − C P,i =: κ f i .
Then, the closed-loop system is given by Moreover, (19) is EIP with supply rate (z i −z i )(d i −d i ) and positive definite storage function for any (feasible) equilibrium pair (d i ,z i ) and associated equilibrium statē Proof: See Appendix C. Remark 15: The controller (17)-which is inspired by [48, Th. 2]-will be used in the sequel to regulate the flow through some DGUs and some consumers. Note that the dynamics (5) of DGUs and consumers are equivalent to the dynamics (4) of a pipe element in series with a pump, excluding the control valve.

C. Volume Flow Rate Control of Valves
Problem 1 also considers the regulation of volume flow rates q i through DGUs, consumers, and mixing connections via control valves. Note that the model of a control valve in series with a pipe element is equivalent to (6), but with d i treated as an arbitrary external input.
Proposition 4: Consider the model of any control valve in series with a pipe element (see (6)). Let and assign the control input u v,i as with volume flow rate setpoint q * i and control parameters Q I,i , K P,i > 0. Then, the closed-loop system is given by d dt Moreover, (24) is EIP with supply rate (z i −z i )(d i −d i ) and positive definite storage function for any (feasible) equilibrium pair (d i ,z i ) and associated equilibrium statē Proof: See Appendix D Remark 16: The design of the PI controller (23) is based on the observation that the dynamics of a control valve in series with a pipe are EIP with respect to the (control) inputoutput pair (u v,i ,ŷ i ). The definition of the outputŷ i originates from the fact that the input matrix G i (x i ) as in (6) is statedependent. This fact complicates the use of a standard PI controller around the shifted passive output (see also [36, p. 137] [30] [46, p. 26]). Following [30,Eq. (8)], we instead proposeŷ i as a new passive output, which is obtained from a suitable, shifted representation of the dynamics.
In steady-state,ŷ i = 0 generally allows for eitherq i = q * i or q i = 0 (see (22) and (24)). However,q i = 0 implies λ i (0) = 0, µ i (0) = 0 (see (2c) and Assumption 1) and thusd i = 0, where d i is the pressure difference over the serial connection of valve and pipe element (cf. (6)). This makes sense from a practical perspective, as a control valve with zero differential pressure available is not functional. In practice, such a situation is thus ruled out by a proper assignment of the pressure setpoints via a higher level control (see Remark 12), which motivates the following assumption.
Assumption 3: Any control valve in series with a pipe element (see (6)) is assumed to have a positive differential pressure d i > 0 for all time. This implies for DGUs and consumer circuits i ∈ D∪L that p j +p P,i −p k > 0 (see Figs. 5 and 7) and for mixing connections i ∈ M that p j − p k > 0 (see Fig. 8).

D. Properties of the Closed-Loop Systems
In this section, we deploy the controlled pumps and valves from Propositions 2-4 in the corresponding actuated edges and nodes i ∈ D ∪ L ∪ P boost ∪ M∪ ∈ H. In particular, we show that these closed-loop systems are also EIP, a fact which is central for the subsequent stability analysis of the overall, interconnected DHN model.
Lemma 1: Assign the pump controllers (12), (17), and the valve controller (23) to the respective edge and node models i ∈ D ∪ L ∪ P boost ∪ M ∪ H according to the control tasks in Problem 1. Then, the resulting closed-loop subsystems can be written asẋ with appropriate vectors and matrices. Moreover, for any i ∈ D ∪ L ∪ P boost ∪ M ∪ H, (27) is EIP with supply rate (z i − z i )(d i −d i ) and positive definite storage function whereQ i is a suitable positive definite diagonal matrix, and x i is any (feasible) equilibrium value ofx i with associated d i . In addition, under Assumption 3, the equilibrium valuex i of any i ∈ D ∪ L ∪ P boost ∪ M ∪ H is such thatq i = q * i and p P,i = p * P,i in accordance with Problem 1. Proof: See Appendix E Remark 17: Note that in Lemma 1 we have intentionally kept the original names for the interaction input and output variables d i and z i . This has been done to emphasize that our control designs and the resulting closed-loop systems have not altered the description of the physical interconnection among the DHN subsystems.

VI. MODULAR STABILITY ANALYSIS
In this section, we prove that the overall closed-loop dynamics of the DHN admit an asymptotically stable hydraulic equilibrium.

A. EIP of the Unactuated Edges and Nodes
Firstly, we show that the unactuated pipe edges i ∈ P \P boost and the capacitive nodes i ∈ C ∪ K are EIP as well, although with non-assignable steady-states.
Lemma 2: The model of any i ∈ P \ P boost ∪ C, i.e., any pipe edge without booster pump and any capacitive node, can be written as (27) with an appropriate choice of the defining vectors and matrices. Furthermore, these models are EIP with some positive definite storage function in the form of (28).
Proof: See Appendix F.

B. Interconnection Structure
Next, we illustrate the power-preserving nature of the interconnection structure that describes the interaction between edge and node subsystems of the DHN. In order to facilitate subsequent elaborations and provide a compact representation, we assemble the overall system dynamics in vector form. Recall that the set of edges E is the ordered union of D, L, P and M. Let ∆ be the ordered union of the sets of nodes H and C. For andĤ i are as in Lemmas 1 and 2. Then, the overall, closed-loop DHN can be written aṡ Let us denote by B the incidence matrix of the DHN digraph G (see Section II-A). With the considered ordering of the edges and nodes of G, the incidence matrix can be written as With (30), the interaction inputs of each of the subsystems in (29) can be represented as Note, e.g., that for any i ∈ E, (5). That is, d i is expressed as the product of an appropriate column of B and the vector z ⊤ Consider the following lemma, which will be useful later to analyze the stability of (29).
Proof: As (31) is a skew-symmetric interconnection structure, it holds for all time that Hence, the interconnection among the subsystems in (29) is power-preserving. 7 7 The equivalence (32) also holds if the respective (interaction) input-output pairs are shifted with respect to any feasible equilibrium values.
The system conformed by (29) and (31) is an index-2 differential algebraic system in Hessenberg form (see, e.g., [49]). Based on [49], the time derivative of the algebraic constraint (29c) can be computed once to explicitly obtain z K in terms ofx E andx ∆ as follows: where we have used the fact thatK K is an identity matrix andĈ EKE = diag( 1 Ji ) i∈E . Then, (29), (31) is equivalent to the ODĖ defined on the invariant manifold The invertibility-and also positive-definiteness-of (33) follows from the fact that B KE B ⊤ KE is a principal submatrix of the DHN's graph's Laplacian L = BB ⊤ . This submatrix is invariant under removal of any row of B not associated to K. Removal of any such row can be understood as a grounding of G (see [50] for a definition). Laplacians of grounded connected graphs, or grounded Laplacians, are positive definite [50]. Note directly that B KE ΘB ⊤ KE is positive definite for any symmetric, positive definite matrix Θ.

C. Asymptotic Stability of the Hydraulic DHN Equilibrium
As a last step, we combine the EIP properties of the DHN subsystems analyzed in Lemmas 1 and 2, the powerpreserving property of their interconnection structure shown in Lemma 3, and LaSalle's Invariance principle to prove asymptotic stability of any feasible, forced hydraulic DHN equilibrium in a modular, bottom-up manner.
Theorem 1: Consider a DHN as described in Section III with pumps and valves controlled as proposed in Lemma 1. Then, under Assumption 3, such a DHN admits an asymptotically stable hydraulic equilibrium (x E ,x ∆ ) in which the pump pressures as well as the pump and valve volume flow rates take on the desired steady states specified in Problem 1. That is, for any i ∈ D ∪ L ∪ P boost ∪ M ∪ H, the respectivē x i in (x E ,x ∆ ) are such thatq i =q P,i = q * i andp P,i = p * P,i in accordance with Problem 1.
Proof: See Appendix G.

VII. SIMULATION
In this section, we demonstrate the stabilizing properties, plug-and-play capabilities, and robustness of the proposed pressure and volume flow rate controllers via simulations in MATLAB/SIMULINK using SIMSCAPE components. In Section VII-A, we present a scenario with plug-and-play operations and varying reference values. In Section VII-B, the first scenario is repeated, albeit with parameter uncertainties and a saturation to the valve input.
The simulations are conducted by means of the DHN depicted in Fig. 1 which shows all structural features discussed in Sections II and III. Furthermore, we cover all control problems outlined in Problem 1 by assigning appropriate DGU, consumer, and pressure holding configurations, i.e., The model and controller parameters used in the simulations are reported in Tables I and II. The pump parameters follow from [33,Eq. (42)] by considering that RP,i JP,i = 7.2878, 1 JP,iCP,i = 341.4283, and setting R P,i = 1 · 10 6 Pa s/m 3 for pressure-controlled pumps, i ∈ D form ∪ D valve ∪ L boost ∪ P boost ∪ H and R P,i = 1 · 10 10 Pa s/m 3 for volume flow rate-controlled pumps, i ∈ D VSP ∪ L VSP . 8 The flow capacity of the control valves is obtained by considering a maximum volume flow rate of k vs = 90 m 3 /h = 0.025 m 3 /s at full valve opening s v,i = 1 → f v,i (s v,i = 1) = 1 and µ i (s v,i = 1, k vs = 0.025 m 3 /s) = 1 · 10 5 Pa valve pressure (see also [6, p. 144]). The value in Table I then follows from Any pipe resistance λ i (q i ) and fluid inertia J i , i ∈ E, are modelled by using the hydraulic pipe resistance and hydraulic fluid inertia SIMSCAPE components with standard values and no elevation. The diameters, roughness, and lengths given in Tables I and II are in line with typical values (see [51], [35]) and correspond to DN 32 and DN 80 pipes, respectively. The elasticities lumped into the capacitive nodes j ∈ C are chosen according to exemplary values for 1-family installations given in [38], [41]. Note that for the pressure-controlled pump, we follow Remark 14 and do not assign any damping, i.e., R P i = R P,i .

A. Scenario A: Plug-and-Play and Setpoint Changes
In this scenario, the simulation starts off with DGU 3 disconnected. The pressure and volume flow rate setpoints for the controlled pumps and valves are assigned as in Table III. At the indicated times, the following events occur: • t = 5 s: Consumers {6, 7, 8} increase their their volume flow rates by 100 % until t = 10 s. • t = 20 s: To help cover the increased demand, DGU 3 connects and the mixing connection 23 increases its volume flow rate to 3 · 10 −3 m 3 /s. • t = 30 s: DGU 2 increases its input volume flow rate to 4.5 · 10 −3 m 3 /s. • t = 40 s: Consumer 4 disconnects. The pressure and volume flow rate trajectories shown in Figs. 9 and 10 confirm the theoretical stability statements. Despite plug-and-play operations and changing operating conditions, 8 Note that these hydraulic resistances values are in the order of magnitude of DN 80 and DN 20 pipes, respectively, for a length of 1 m and an external pressure of 1 · 10 5 Pa. J P,i = 1.37 · 10 5 Pa s 2 /m 3 C P,i = 2.13 · 10 −8 m 3 /Pa Flow-controlled R P,i = 1 · 10 10 Pa s/m 3 pumps (1) J P,i = 1.37 · 10 9 Pa s 2 /m 3 C P,i = 2.13 · 10 −12 m 3 /Pa Control valves (2) C v,i = 7.9 · 10 −     Table III For the volume flow rates, larger deviations can be observed. In particular at t = 20 s and t = 30 s during the connection of DGU 3 and the setpoint changes the error plots in Fig. 10 shows large outliers. However, from a practical perspective, this is natural as abrupt setpoint changes cannot be realized instantly by the respective volume flow rate controllers. More importantly, except for the load ramps at t ∈ [5 s, 10 s], the volume flow rates settle to within a 1.5 % band with respect to the setpoints after at most 5 s. During the load ramps, the errors are slightly higher, but remain below 8 %. This shows that the volume flow rate controllers for both pumps and valves, although not specifically designed for it, are sufficiently fast to adquately track setpoints that vary on a time scale of seconds.

B. Scenario B: Parameter Uncertainty and Valve Saturation
Scenario B is similar to Scenario A except for two modification: firstly, a 10 % uncertainty is added to the pump parameters R P,i , J P,i , C P,i and the valve parameters C v,i ; secondly, the virtual valve control input of all valves is saturated to u v ∈ 1, u max v,i (see Assumption 2). 9 9 In line with classical feedback control design, we did not consider the possibility of control input saturation explicitly during the control design stage in Section V. Instead, we analyze its impact by means of the numerical simulation in this section.  Table III The resulting pressure and volume flow rate trajectories shown in Figs. 11 and 12 are similar to those shown in Figs. 9. The main difference introduced by the valve saturation is an impaired convergence performance of the volume flow rate control via valves, particularly at DGU 2 and Consumer 6 (see the orange and black lines in Fig. 12). In practice, an appropriate redesign of the valves or by increasing the available pressure, e.g., via the booster pump in Pipe 15 or a separate booster pump in the respective consumer, control performance can be improved.
Overall, the results of the two scenarios illustrate that the passivity-based pressure and volume flow rate controllers indeed asymptotically stabilize the hydraulic variables while allowing for plug-and-play operations of the different DHN subsystems. Furthermore, the integral parts of the proposed controllers provides robustness against parameter uncertainties and changing hydraulic conditions naturally occuring during the operation of DHNs.
VIII. CONCLUSION In this work, we have proposed a unifying control framework that guarantees pressure and volume flow rate stability based on the EIP properties of the DHN subsystem models. We provided a comprehensive hydraulic model covering stateof-the-art as well as future DHN generations and formalized the hydraulic control problems arising in such systems. We demonstrated how each of the subsystem models can be represented in port-Hamiltonian form to facilitate the subsequent passivity-based control design and stability analysis. Subsequently, we showed how the hydraulic control tasks can be simplified to designing decentralized, passivity-based controllers for the pumps and valves in a DHN and proposed three controllers: two for pressure and volume flow rate control via pumps and one for volume flow rate control via valves. Based on the EIP properties of the (actuated and unactuated) subsystem models, the skew-symmetry of their interconnection structure, and LaSalle's Invariance principle, we then proved asymptotic stability of any forced, hydraulic DHN equilibrium in a modular manner. In conclusion, we want to highlight that the presented EIP-based control framework is not restricted to the models and controllers proposed in this work. In fact, the modular approach of the EIP-based stability analysis allows to incorporate more detailed models (e.g., for the pumps and valves) and presents general guidelines for developing other decentralized pressure and volume flow rate controllers. Proof: To see that (i) holds, consider first the pump dynamics in (1). By direct substitution, it can be verified that this model can be written as an ISO-PHS by defining Moving on, the model of any pipe in series with a valve is equivalent to (6). Then, the ISO-PHS representation is attained by defining with (j, k) ∈ N − i × N + i . Now, we show that (ii) holds. Any DGU or consumer i ∈ D∪C is modeled by (5). Then, ISO-PHS representation is attained by defining where (j, k) ∈ N − i × N + i . The model (4) of each pipe i ∈ P can be written as a ISO-PHS by defining where (j, k) ∈ N − i × N + i . Each mixing valve i ∈ M modeled by (6) can be represented as an ISO-PHS by defining where (j, k) ∈ N − i × N + i . Any pressure holding unit j ∈ H modeled by (8) can be represented as an ISO-PHS by defining x j = J P,j q P,j C P,j p P,j ⊤ , u j = u P,j , ζ i = q P,j , d j = i∈Ij q i , z j = p P,j , where I j ⊆ E is the set of edges that are incident to j.
Finally, each capacitive node j ∈ C modeled by (9) can be written as an ISO-PHS by defining where I j is defined as before.

APPENDIX B PROOF OF PROPOSITION 2
Following [29], we begin by introducing a change of coordinates from q P,i to χ i as (see (12)) Then, by computingχ i , we get where ν i is given in (12). Following the IDA-PBC design methodology [28], we assign u P,i as in (12) to obtain By propagating the coordinate transformation (44) to the p P,idynamics, we can write the closed-loop dynamics as in (13). In order to show that (13) is EIP with supply rate (d i −d i )(z i − z i ) and positive definite storage function H p i in (14), let d i be fixed to an arbitrary equilibrium valued i with associated equilibriaz i for the output andx p i (see (15)) for the state vector. Sincex p i satisfies f p i (x p i ) + K p id i = 0, we can write (13) equivalently aṡ For the time derivative of H p i in (14), it holds thaṫ where we have used the identity By combining (4) with the controller (17), the closed-loop system (19) follows directly. In order to show that system (19) is EIP with supply rate (d i −d i )(z i −z i ) and positive definite storage function H f i in (20), let d i be fixed to an arbitrary equilibrium valued i with associated equilibriaz i for the output andx f i (see (21)) for the state vector. Sincex f i satisfies f i (x f i )+ K f id i = 0, we can write (19) equivalently aṡ For the time-derivative of H f i in (20), it holds thaṫ where we have used the identity Following the same reasoning as in the proof of Proposition 3 and usingū i =r i , ±μ i (q i )r i , we can write the closedloop system (24) equivalently as Note that (51) is equivalent tȯ x for any arbitrary equilibrium pair (d i ,z i ) and associated equilibrium state vectorx v i (see (26)). For the time-derivative of H v i in (25), it holds thaṫ where we have used the identity Since both λ i (q i ) andμ i (q i ) are strictly increasing (see Assumption 1 and (2c)),ū i =r i > 0 per definition (see Section III-A2), and K P,i > 0, it follows that In the following, we successively consider the closed-loop systems according to their order in Problem 1 (a)-(i). Thus, we begin with i ∈ D form . By combining the open-loop DGU model (5) with u P,i as in (12) and fixing u v,i =ū v,i > 0 (due to s v,i = 1), we can write the closed-loop as in (27) d dt with d i as in (5). To show that (55) is EIP, we follow the same reasoning as in the proofs of Proposition 2 and 3. The storage functionĤ i is as in (28) witĥ The time derivative ofĤ i along solutions of (55) satisfieṡ where we have used the identity As λ i andμ i are strictly increasing andū v,i , is EIP. Finally, observe that due to the integral action (12b), it holds for any feasible equilibrium value ofx i that p P,i = p * P,i .
Next, we consider i ∈ D valve . By combining the openloop DGU model (5) with u P,i and u v,i as in (12) and (23), respectively, we can write the closed-loop as in (27) with d i as in (5). Note that we have used the indices α and β to distinguish between the integral actions of (12) and (23), respectively. To show that (59) is EIP, we proceed as before.
The storage functionĤ i is as in (28) witĥ The time derivative ofĤ i along solutions of (59) satisfieṡ where we have used the identity With the same reasoning as for (54) and R P i > 0, it follows thatψ i (x i ) ≥ 0 andḢ i ≤ (z i −z i )(d i −d i ). Hence, EIP is established. Finally, observe that due to the integral actions (12b) and (23a) and Assumption 3, it holds for any feasible equilibrium value ofx i thatq i = q * i andp P,i = p * P,i .
For any i ∈ D VSP , we combine the open-loop DGU model (5) with u P,i as in (3) and fix u v,i =ū v,i > 0 (due to s v,i = 1) to write the closed-loop as in (27) d dt     J i q i J P,i q P,i C P,i p P,i Q I,i r i with d i as in (5). To show that (63) is EIP, we note thatĤ i is given by (28) withQ i = Q f i from (20). The time derivative ofĤ i along solutions of (63) satisfieṡ where we have used the identity As λ i andμ i are strictly increasing,ū v,i > 0, and (18) holds, it follows thatψ i ( is EIP. Finally, observe that due to the integral action (17a), it holds for any feasible equilibrium value ofx i that q i =q P,i = q * i . For describing the consumers L in closed-loop, we recall that the open-loop model of any i ∈ L is identical to that of any DGU i ∈ D, i.e., to (5). For any i ∈ L boost , we assign u P,i and u v,i as in (12) and (23), respectively. Then, the resulting closed-loop is equivalent to (59) with the same EIP and equilibrium properties. Next, we consider i ∈ L valve . In this case, the associated pump can either be turned off or is not present. Thus, the open-loop of any i ∈ L valve is equivalent to that of a control valve in series with a pipe element, i.e., to (6). By closing the loop with u v,i as in (23), we obtain a closed-loop system equivalent to (24), which is clearly of the form (27). Moreover, EIP andq i = q * i directly follow from Proposition 4 and Assumption 3. Lastly, for any i ∈ L VSP , we fix u v,i =ū v,i > 0 and assign u P,i as in (3). Then, the resulting closed-loop is equivalent to (63) with the same EIP and equilibrium properties.
Next, we consider i ∈ P boost , i.e., an arbitrary pipe in series with a booster pump. By combining the open-loop pipe model (4) with u P,i as in (12), we can write the closed-loop dynamics as in (27) d dt     J i q i J P,i χ i C P,i p P,i Q I,i r i with d i as in (4). To show that (66) is EIP, we note thatĤ i is given by (28) witĥ Q i = diag −1 (J i , J P,i , C P,i , Q I,i ).
The time derivative ofĤ i along solutions of (66) satisfieṡ where we have used the identity As λ i is strictly increasing and R P i > 0, it follows that Hence, (66) is EIP. Finally, observe that due to the integral action (12b), it holds for any feasible equilibrium valuex i thatp P,i = p * P,i . The final type of actuated edges corresponds to mixing valves i ∈ M. By combining the open-loop model (6) with u v,i as in (23), we obtain a closed-loop system equivalent to (24), which is clearly of the form (27). Moreover, EIP and q i = q * i directly follow from Proposition 4 and Assumption 3. Finally, we consider pressure holding units i ∈ H, which are the only actuated nodes. By combining the open-loop model (8) with u P,i as in (12), we obtain a closed-loop system equivalent to (13). Clearly, this system can be written as in (27). Moreover, EIP andp P,i = p * P,i follow directly from Proposition 2.

APPENDIX F PROOF OF LEMMA 2
For any pipe i ∈ P \ P boost without booster pump, we use the open-loop model (4) and set u P,i = 0 to obtain a model of form (27).
with d i as in (4). EIP of (70) follows directly by usingĤ i in (28) as storage function withQ i = 1 Ji and noting thaṫ Writing the model (9) of any capacitive node i ∈ C as in (27) is trivial. Furthermore, EIP of system (9) follows directly by usingĤ i in (28) as a storage function withQ i = 1 Ci and noting thatḢ APPENDIX G PROOF OF THEOREM 1 Let (x E ,x ∆ ) denote a feasible equilibrium of the closedloop system (34). Due to Lemma 1, this equilibrium is such that each bullet point in Problem 1 is fulfilled. Consider the following Lyapunov function candidate: where eachĤ i is defined in Lemmas 1 or 2. Note that V (x E ,x ∆ ) > 0 for all (x E ,x ∆ ) = (x E ,x ∆ ) and V (x E ,x ∆ ) = 0. The time derivative of V along solutions of (34) on the invariant set M in (35) is given bẏ To get this equality, we have used three facts about any i ∈ E ∪ ∆. (I) That for any any feasible equilibrium input d i , Stability of (x E ,x ∆ ) follows directly from (74), as its righthand side is upper bounded by zero according to the EIP properties established in Lemmas 1 and 2.
We move on to show that (x E ,x ∆ ) is in fact asymptotically stable. Considering LaSalle's invariance principle, at this point it is sufficient to show that (x E ,x ∆ ) is the largest invariant set of (34) whereV is zero. Next we characterize the conditions under whichV (x) = 0. For that we split the summands in the right-hand side of (74) as we explain next.
Firstly, from the proof of Lemma 1 we have that