Port-Hamiltonian Systems on Graphs

1. PORT-HAMILTONIAN SYSTEMS ON GRAPHS*

A. J. VAN DER SCHAFT† AND B. M. MASCHKE‡

Abstract. In this paper we present a unifying geometric and compositional framework for modeling complex physical network dynamics as port-Hamiltonian systems on open graphs. The basic idea is to associate with the incidence matrix of any directed graph a Dirac structure relating the flow and effort variables associated to the edges and vertices of the graph, and to formulate energy-storing or energy-dissipating relations between the flow and effort variables of the edges and the internal vertices. This allows for state variables associated to the edges and formalizes the interconnection of networks. Examples from different origins, such as consensus algorithms, that share the same structure are shown. It is shown how the identified Hamiltonian structure offers systematic tools for the analysis and control of the resulting dynamics.

Key words. physical systems, Hamiltonian dynamics, Dirac structures, network dynamics, stability, symmetry reduction

AMS subject classifications. 05C21, 37J99, 53D20, 70H05, 93A15, 93A30, 93C15, 93D20

DOI. 10.1137/110840091

1. Introduction. Discrete topological structures arise abundantly in physical systems modeling. The classical approach to the analysis of electrical circuits, dating back to Kirchhoff, is based on the circuit graph. Similar approaches apply to many other cases, including mass-spring-damper mechanical systems, multibody systems, hydraulic networks, chemical reaction networks, and power systems. A common feature is that the discrete structures, in particular graphs, are blended with dynamical relations, leading to various sorts of network dynamics.

During the last two decades the study of network dynamics has received ever-increasing attention, with input from, among others, the fields of graph theory, multiagent systems, dynamical systems, and statistical mechanics. In this paper we formulate a general geometric framework for defining physical dynamics on directed open graphs.¹ The generalized Hamiltonian nature of the resulting dynamical models is due to the assumption that the constitutive relations between the variables corresponding to storage at the vertices and/or edges are derivable from an energy (Hamiltonian) function, while the remaining variables are related by static energy-dissipating relations. This will imply that the total energy itself satisfies a conservation law: the increase of the total energy is equal to externally supplied power (through the boundary vertices of the graph) minus the power lost in the dissipative elements (associated to some of the edges or vertices of the graph). The resulting generalized Hamiltonian systems, allowing for energy-dissipation and interaction with the environment, fall within the class of port-Hamiltonian systems, as coined and explored in, e.g.,

*Received by the editors July 8, 2007; accepted for publication (in revised form) November 28, 2012; published electronically March 14, 2013.

http://www.siam.org/journals/sicon/51-2/84009.html

†Johann Bernoulli Institute for Mathematics and Computer Science, University of Groningen, 9700 AK, Groningen, The Netherlands (A.J.van.der.Schaft@rug.nl). The research of this author was supported by the European Union Seventh Framework Programme [FP7/2007–2013] under grant agreement n257462 HYCON2 Network of Excellence.

‡Laboratoire d'Automatique et de Génie des Procédés, Université Claude Bernard Lyon-1, F-69622 Villeurbanne, Cedex, France (maschke@lagep.univ-lyon1.fr).

¹Note that this does not include the (random) evolution of the graphs themselves, as studied in random graph theory and statistical mechanics.

[36, 10, 32, 37, 14].

From a geometric point of view the generalized Hamiltonian structure of the network dynamics is defined, apart from its Hamiltonian function and energy-dissipating relations, by a Dirac structure. This Dirac structure (generalizing the symplectic or Poisson structure from classical mechanics) is directly defined by the incidence matrix of the directed graph, and thus captures the conservation laws. In fact, we will show how a directed graph gives rise to three canonically defined Dirac structures on its vertex and edge spaces. The first two differ only in the different role of the boundary vertices, while the third, the Kirchhoff–Dirac structure, captures the special case where no storage or dissipation is associated with the vertices of the graph (corresponding to Kirchhoff’s current laws).

We will illustrate this framework on some of the physical examples mentioned above. Furthermore, we will show how the same port-Hamiltonian structure is shared by network dynamics with a different origin, such as consensus and clustering algorithms, and how the identification of the underlying port-Hamiltonian structure provides powerful tools for analysis and control, which unify and go beyond existing approaches.

While all examples given in the paper are simple, and could be approached from other angles as well, we believe that a major contribution of the paper resides in the identification of a common mathematical structure in all of these examples, which is, moreover, closely related to classical Hamiltonian systems. Furthermore, the approach formalizes network dynamics as open systems, while due to the compositionality properties of port-Hamiltonian systems, it is easily scalable and extends to heterogeneous and multiscale systems as well.

In a companion paper we will describe how the geometric framework as developed in this paper for graphs can be extended to arbitrary $k$ -complexes. Among others, this will allow for a structure-preserving spatial discretization of distributed-parameter physical systems, otherwise described by partial differential equations; see [38, 39].

Preliminary work regarding sections 3.4 and 3.5 can be found in [40, 38, 39].

2. From directed graphs to Dirac structures. As a guiding example let us consider a mass-spring-damper system, for example, the one depicted in Figure 1.

The underlying directed graph of such a system is defined by vertices corresponding to the masses and edges corresponding to the springs and dampers, leading to the graph in Figure 2.

How do we formalize such a system as a port-Hamiltonian system? A key ingredient in the definition of a port-Hamiltonian system is the geometric notion of a Dirac structure, generalizing the symplectic structure from classical Hamiltonian dynamics. In section 2.4 we will define two canonical Dirac structures on the combination of the vertex, edge, and boundary spaces of a directed graph, and their dual spaces. These two Dirac structures will differ only in the role of the boundary vertices, which for a mass-spring-damper system either will be associated to boundary masses (with inputs being the external forces on them) or will be massless (with inputs being their velocities).

graph LR
    m1[m1] ---|spring 1| m2[m2]
    m1 ---|damper 1| m2
    m2 ---|damper 2| m3[m3]
    m1 ---|spring 2| m3

Diagram of a mass-spring-damper system with three masses m1, m2, and m3. Mass m1 is connected to mass m2 by a spring labeled 'spring 1' and a damper labeled 'damper 1'. Mass m2 is connected to mass m3 by a damper labeled 'damper 2'. Mass m1 is also connected to mass m3 by a spring labeled 'spring 2'.

FIG. 1. Mass-spring-damper system.

A directed graph with three vertices arranged horizontally. The leftmost and middle vertices are connected by two curved edges: one above and one below. The middle and rightmost vertices are connected by a single straight edge. All three vertices are marked with a circle containing a cross, indicating they are boundary vertices.

FIG. 2. The corresponding graph.

We first recall some basic notions of graph theory (see, e.g., [4]) and Dirac structures (see, e.g., [9, 13, 10]).

2.1. Directed graphs and their vertex and edge spaces. A directed graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ consists of a finite set $\mathcal{V}$ of vertices (nodes) and a finite set $\mathcal{E}$ of directed edges (branches or links), together with a mapping from $\mathcal{E}$ to the set of ordered pairs of $\mathcal{V}$ , where no self-loops are allowed. Thus to any branch $e \in \mathcal{E}$ there corresponds an ordered pair $(v, w) \in \mathcal{V} \times \mathcal{V}$ (with $v \neq w$ ), representing the tail vertex $v$ and the head vertex $w$ of this edge.

A directed graph is completely specified by its incidence matrix $\hat{B}$ , which is an $N \times M$ matrix, where $N$ is the number of vertices and $M$ is the number of edges, with the $(i, j)$ th element equal to $-1$ if the $j$ th edge is an edge toward vertex $i$ , equal to $1$ if the $j$ th edge is an edge originating from vertex $i$ , and $0$ otherwise. Since we will consider only directed graphs in what follows “graph” will mean “directed graph.”

Given a graph, we define its vertex space $\Lambda_0$ as the vector space of all functions from $\mathcal{V}$ to some linear space $\mathcal{R}$ . In the examples, $\mathcal{R}$ will be mostly $\mathcal{R} = \mathbb{R}$ or $\mathcal{R} = \mathbb{R}^3$ . In the first case, $\Lambda_0$ can be identified with $\mathbb{R}^N$ . Furthermore, we define its edge space $\Lambda_1$ as the vector space of all functions from $\mathcal{E}$ to the same² linear space $\mathcal{R}$ . Again, if $\mathcal{R} = \mathbb{R}$ , then $\Lambda_1$ can be identified with $\mathbb{R}^M$ .

The dual spaces of $\Lambda_0$ and $\Lambda_1$ will be denoted by $\Lambda^0$ and $\Lambda^1$ , respectively. The duality pairing between $f \in \Lambda_0$ and $e \in \Lambda^0$ is given as

$\langle f | e \rangle = \sum_{v \in \mathcal{V}} \langle f(v) | e(v) \rangle,$

where $\langle \cdot | \cdot \rangle$ on the right-hand side denotes the duality pairing between $\mathcal{R}$ and $\mathcal{R}^*$ , and a similar expression holds for $f \in \Lambda_1$ and $e \in \Lambda^1$ (with summation over the edges).

The incidence matrix $\hat{B}$ of the graph induces a linear map $B$ from the edge space to the vertex space as follows. Define $B : \Lambda_1 \rightarrow \Lambda_0$ as the linear map with matrix representation $\hat{B} \otimes I$ , where $I : \mathcal{R} \rightarrow \mathcal{R}$ is the identity map and $\otimes$ denotes the Kronecker product. $B$ will be called the incidence operator. For $\mathcal{R} = \mathbb{R}$ the incidence operator reduces to the linear map given by the matrix $\hat{B}$ itself, in which case we will use $B$ for both the incidence matrix and the incidence operator. The adjoint map of

²In principle we could also associate with the edges a linear space $\mathcal{R}'$ which is different from the space $\mathcal{R}$ associated with the vertices. In that case the definition of the incidence operator needs an additional linear map from $\mathcal{R}'$ to $\mathcal{R}$ .

$B$ is denoted as

$B^* : \Lambda^0 \rightarrow \Lambda^1$

and is called the coincidence operator. For $\mathcal{R} = \mathbb{R}^3$ the coincidence operator is given by $\hat{B}^T \otimes I_3$ , while for $\mathcal{R} = \mathbb{R}$ the coincidence operator is simply given by the transposed matrix $B^T$ , and throughout we will use $B^T$ for both the coincidence matrix and the coincidence operator.

We will use the terminology³ flows for the elements of $\Lambda_0$ and $\Lambda_1$ (notation $f_0$ and $f_1$ ) and efforts for the elements of their dual spaces $\Lambda^0$ and $\Lambda^1$ (notation $e^0$ , respectively, $e^1$ ).

2.2. Open graphs. An open graph $\mathcal{G}$ is obtained from an ordinary graph with a set of vertices $\mathcal{V}$ by identifying a subset $\mathcal{V}_b \subset \mathcal{V}$ of $N_b$ boundary vertices. The interpretation of $\mathcal{V}_b$ is that these are the vertices that are open to interconnection (i.e., with other open graphs). The remaining subset $\mathcal{V}_i := \mathcal{V} - \mathcal{V}_b$ consists of the $N_i$ internal vertices of the open graph.

The splitting of the vertices into internal and boundary vertices induces a splitting of the vertex space and its dual, given as

$\begin{aligned}\Lambda_0 &= \Lambda_{0i} \oplus \Lambda_{0b}, \\ \Lambda^0 &= \Lambda^{0i} \oplus \Lambda^{0b},\end{aligned}$

where $\Lambda_{0i}$ is the vertex space corresponding to the internal vertices and $\Lambda_{0b}$ the vertex space corresponding to the boundary vertices. Consequently, the incidence operator $B : \Lambda_1 \rightarrow \Lambda_0$ splits into

$B = B_i \oplus B_b$

with $B_i : \Lambda_1 \rightarrow \Lambda_{0i}$ and $B_b : \Lambda_1 \rightarrow \Lambda_{0b}$ .

Furthermore, we will define the boundary space $\Lambda_b$ as the linear space of all functions from the set of boundary vertices $\mathcal{V}_b$ to the linear space $\mathcal{R}$ . Note that the boundary space $\Lambda_b$ is isomorphic to the linear space $\Lambda_{0b}$ , and that by using this isomorphism the linear mapping $B_b$ can also be regarded as a mapping

$B_b : \Lambda_1 \rightarrow \Lambda_b$

called the boundary incidence operator. Nevertheless, we will be careful in distinguishing the two isomorphic linear spaces $\Lambda_b$ and $\Lambda_{0b}$ because of their different interpretations in physical examples (e.g., for mass-spring-damper systems $\Lambda_b$ will denote the space of external forces as exerted on the boundary masses, and $\Lambda_{0b}$ will denote the space of momenta of the boundary masses). The dual space of $\Lambda_b$ will be denoted as $\Lambda^b$ . The elements $f_b \in \Lambda_b$ are called the boundary flows and the elements $e^b \in \Lambda^b$ the boundary efforts.

³This terminology stems from port-based and bond-graph modeling [27], where it has a slightly more specific connotation than in our case. The space $\Lambda_0$ is also called the space of 0-chains, while the elements of $\Lambda_1$ are called the 1-chains. Furthermore, the dual spaces $\Lambda^0$ and $\Lambda^1$ are called the space of 0-cochains, respectively, 1-cochains. In [23] this will be generalized to higher-order chains and cochains. In generalized circuit theory, $f_1 \in \Lambda_1$ are referred to as through variables and $e^1 \in \Lambda^1$ as across variables.

2.3. Dirac structures. Recall (see [36, 9, 32]) the definition of a (constant⁴) Dirac structure. Consider a vector space $\mathcal{F}$ with dual space $\mathcal{F}^*$ . As before, the variables $f \in \mathcal{F}$ are called the flow variables, while the conjugate variables $e \in \mathcal{F}^*$ are called the effort variables. On the total space $\mathcal{F} \times \mathcal{F}^*$ define the indefinite inner product $\langle\langle \cdot, \cdot \rangle\rangle$ as

$\langle\langle (f_a, e_a), (f_b, e_b) \rangle\rangle := \langle e_a | f_b \rangle + \langle e_b | f_a \rangle, \quad f_a, f_b \in \mathcal{F}, \quad e_a, e_b \in \mathcal{F}^*,$

where $\langle \cdot | \cdot \rangle$ denotes the duality product between $\mathcal{F}$ and $\mathcal{F}^*$ .

DEFINITION 2.1. A subspace $\mathcal{D} \subset \mathcal{F} \times \mathcal{F}^*$ is a Dirac structure if $\mathcal{D} = \mathcal{D}^\perp$ , where $^\perp$ denotes the orthogonal complement with respect to $\langle\langle \cdot, \cdot \rangle\rangle$ .

In the finite-dimensional case an equivalent, and often easier, characterization of Dirac structures is given as follows (see, e.g., [8, 14] for a proof).

PROPOSITION 2.2. A subspace $\mathcal{D} \subset \mathcal{F} \times \mathcal{F}^*$ is a Dirac structure if and only if the following two conditions are satisfied:

(1) $\langle e | f \rangle = 0 \quad \text{for all } (f, e) \in \mathcal{D},$
(ii) $\dim \mathcal{D} = \dim \mathcal{F}.$

Note that the first equation in (1) can be regarded as a power-conservation property. The second equation states that a Dirac structure has maximal dimension with respect to this power-conserving property [10, 32].

While Dirac structures thus formalize power-conserving interconnections of maximal dimension, the following special type of Dirac structure can be seen to be a generalization of Tellegen's theorem in circuit theory (stating that the product $\langle V_a | I_b \rangle = 0$ for any two vectors of voltages $V_a$ and currents $I_b$ satisfying Kirchhoff's laws).

DEFINITION 2.3. A Dirac structure $\mathcal{D} \subset \mathcal{F} \times \mathcal{F}^*$ is separable if

(2) $\langle e_a | f_b \rangle = 0 \quad \text{for all } (f_a, e_a), (f_b, e_b) \in \mathcal{D}.$

Separable Dirac structures have the following simple geometric characterization, similar to Kirchhoff's current and voltage laws.

PROPOSITION 2.4. Consider a separable Dirac structure $\mathcal{D} \subset \mathcal{F} \times \mathcal{F}^*$ . Then

(3) $\mathcal{D} = \mathcal{K} \times \mathcal{K}^\perp$

for some subspace $\mathcal{K} \subset \mathcal{F}$ , where $\mathcal{K}^\perp = \{e \in \mathcal{F}^* \mid \langle e | f \rangle = 0 \text{ for all } f \in \mathcal{K}\}$ . Conversely, any subspace $\mathcal{D}$ as in (3) for some subspace $\mathcal{K} \subset \mathcal{F}$ is a separable Dirac structure.

Proof. It is immediately seen that any subspace $\mathcal{K} \times \mathcal{K}^\perp$ satisfies (2) and is a Dirac structure since it satisfies (1). Conversely, let the Dirac structure $\mathcal{D}$ satisfy (2). Define the following subspaces:

$\begin{aligned} \mathcal{F}_0 &= \{f \in \mathcal{F} \mid (f, 0) \in \mathcal{D}\}, & \mathcal{F}_1 &= \{f \in \mathcal{F} \mid \exists e \in \mathcal{F}^* \text{ such that } (f, e) \in \mathcal{D}\}, \\ \mathcal{E}_0 &= \{e \in \mathcal{F}^* \mid (0, e) \in \mathcal{D}\}, & \mathcal{E}_1 &= \{e \in \mathcal{F}^* \mid \exists f \in \mathcal{F} \text{ such that } (f, e) \in \mathcal{D}\}. \end{aligned}$

⁴This definition can be extended [13, 9] to (nonconstant) Dirac structures on manifolds: a Dirac structure $\mathcal{D}$ on a manifold $\mathcal{M}$ is defined as a vector subbundle of the Whitney sum $T\mathcal{M} \oplus T^*\mathcal{M}$ such that for each $m \in \mathcal{M}$ the linear space $D(m) \subset T_m\mathcal{M} \times T_m^*\mathcal{M}$ is a constant Dirac structure. This will be needed in the treatment of spatial mechanisms in section 3.3.

It is readily seen [10] that for any Dirac structure $\mathcal{E}_1 = (\mathcal{F}_0)^\perp, \mathcal{E}_0 = (\mathcal{F}_1)^\perp$ . We will now show that (2) implies that $\mathcal{F}_0 = \mathcal{F}_1 =: \mathcal{K}$ (and hence $\mathcal{E}_0 = \mathcal{E}_1 =: \mathcal{K}^\perp$ ). Clearly, $\mathcal{F}_0 \subset \mathcal{F}_1$ . Now let $(f_a, e_a) \in \mathcal{D}$ and thus $f_a \in \mathcal{F}_1$ . Then for all $(f_b, e_b) \in \mathcal{D}$ ,

$\langle (f_a, 0), (f_b, e_b) \rangle := \langle e_b | f_a \rangle + \langle 0 | f_b \rangle = \langle e_b | f_a \rangle = 0$

by (2). Hence, also $(f_a, 0) \in \mathcal{D}$ and thus $f_a \in \mathcal{F}_0$ . By definition $\mathcal{F}_0 \times \mathcal{E}_0 \subset \mathcal{D}$ , and hence $\mathcal{K} \times \mathcal{K}^\perp \subset \mathcal{D}$ . Finally, since the dimension of $\mathcal{K} \times \mathcal{K}^\perp$ equals the dimension of $\mathcal{F}$ , equality results. $\square$

A typical instance of a separable Dirac structure, which will be frequently used in the remainder, is the following.

PROPOSITION 2.5. Let $A : \mathcal{V} \rightarrow \mathcal{W}$ be a linear map between the linear spaces $\mathcal{V}$ and $\mathcal{W}$ with adjoint mapping $A^* : \mathcal{W}^* \rightarrow \mathcal{V}^*$ , that is

$(4) \quad \langle w^* | Av \rangle = \langle A^* w^* | v \rangle$

for all $v \in \mathcal{V}, w^* \in \mathcal{W}^*$ (where, as before, $\langle \cdot | \cdot \rangle$ denotes the duality product between the dual spaces $\mathcal{W}$ and $\mathcal{W}^*$ , respectively $\mathcal{V}$ and $\mathcal{V}^*$ ). Identify $(\mathcal{V} \times \mathcal{W})^* = \mathcal{V}^* \times \mathcal{W}^*$ . Then

$(5) \quad \begin{aligned} \mathcal{D} &:= \{(v, w, v^*, w^*) \in (\mathcal{V} \times \mathcal{W}) \times (\mathcal{V}^* \times \mathcal{W}^*) \mid \\ &Av = w, v^* = -A^* w^*\} \end{aligned}$

is a separable Dirac structure.

Proof. Define $\mathcal{K} := \{(v, w) \in \mathcal{V} \times \mathcal{W} \mid Av = w\}$ . Then $\mathcal{K}^\perp = \{(v^*, w^*) \in \mathcal{V}^* \times \mathcal{W}^* \mid v^* = -A^* w^*\}$ . $\square$

A key feature of Dirac structures is that their composition is again a Dirac structure (in contrast with symplectic or Poisson structures, where this is not generally the case). Let $\mathcal{D}_A \subset \mathcal{F}_A \times \mathcal{F}_c \times \mathcal{F}_A^* \times \mathcal{F}_c^*$ and $\mathcal{D}_B \subset \mathcal{F}_B \times \mathcal{F}_c \times \mathcal{F}_B^* \times \mathcal{F}_c^*$ be two Dirac structures with shared space of flow and effort variables $\mathcal{F}_c$ , respectively, $\mathcal{F}_c^*$ . Define their composition as

$(6) \quad \begin{aligned} \mathcal{D}_A \circ \mathcal{D}_B &= \{(f_A, e_A, f_B, e_B) \in \mathcal{F}_A \times \mathcal{F}_B \times \mathcal{F}_A^* \times \mathcal{F}_B^* \mid \exists (f, e) \in \mathcal{F}_c \times \mathcal{F}_c^* \text{ such that} \\ &(f_A, e_A, f, e) \in \mathcal{D}_A, (f_B, e_B, -f, e) \in \mathcal{D}_B\}. \end{aligned}$

It has been shown in [8, 31] that $\mathcal{D}_A \circ \mathcal{D}_B$ is again a Dirac structure. Separable Dirac structures turn out to have the following special compositional property.

PROPOSITION 2.6. Let $\mathcal{D}_A \subset \mathcal{F}_A \times \mathcal{F}_c \times \mathcal{F}_A^* \times \mathcal{F}_c^*$ and $\mathcal{D}_B \subset \mathcal{F}_B \times \mathcal{F}_c \times \mathcal{F}_B^* \times \mathcal{F}_c^*$ be two separable Dirac structures given as

$\mathcal{D}_i = \mathcal{K}_i \times \mathcal{K}_i^\perp \quad i = A, B,$

where $\mathcal{K}_i \subset \mathcal{F}_i \times \mathcal{F}_c, i = A, B$ . Define the composition

$\mathcal{K}_A \circ \mathcal{K}_B = \{(f_A, f_B) \in \mathcal{F}_A \times \mathcal{F}_B \mid \exists f \in \mathcal{F}_c \text{ such that } (f_A, f) \in \mathcal{K}_A, (f_B, -f) \in \mathcal{K}_B\}.$

Then the composition $\mathcal{D}_A \circ \mathcal{D}_B$ is the separable Dirac structure

$(7) \quad \mathcal{D}_A \circ \mathcal{D}_B = (\mathcal{K}_A \circ \mathcal{K}_B) \times (\mathcal{K}_A \circ \mathcal{K}_B)^\perp.$

For explicit equational representations of compositions of Dirac structures we refer the reader to [8].

The compositionality property of Dirac structures is a key ingredient of port-Hamiltonian systems theory, and implies that the standard interconnection of port-Hamiltonian systems results in another port-Hamiltonian system with Dirac structure being the composition of the Dirac structures of the component port-Hamiltonian systems, and Hamiltonian equal to the sum of the Hamiltonians of the component systems [31, 8].

2.4. The graph Dirac structures. We now have all ingredients to define Dirac structures corresponding to the incidence structure of a directed graph.

DEFINITION 2.7. Consider an open graph $\mathcal{G}$ with vertex, edge, and boundary spaces, incidence operator $B$ , and boundary incidence operator $B_b$ . The flow-continuous⁵ graph Dirac structure $\mathcal{D}_f(\mathcal{G})$ is defined as

$\begin{aligned} \mathcal{D}_f(\mathcal{G}) &:= \{(f_1, e^1, f_{0i}, e^{0i}, f_b, e^b) \\ (8) \quad &\in \Lambda_1 \times \Lambda^1 \times \Lambda_{0i} \times \Lambda^{0i} \times \Lambda_b \times \Lambda^b \mid \\ &B_i f_1 = f_{0i}, B_b f_1 = f_b, e^1 = -B_i^* e^{0i} - B_b^* e^b\}. \end{aligned}$

The effort-continuous graph Dirac structure $\mathcal{D}_e(\mathcal{G})$ is defined as

$\begin{aligned} \mathcal{D}_e(\mathcal{G}) &:= \{(f_1, e^1, f_0, e^0, f_b, e^b) \\ (9) \quad &\in \Lambda_1 \times \Lambda^1 \times \Lambda_0 \times \Lambda^0 \times \Lambda_b \times \Lambda^b \mid \\ &B_i f_1 = f_{0i}, B_b f_1 = f_{0b} + f_b, e^1 = -B^* e^0, e^b = e^{0b}\}. \end{aligned}$

By Proposition 2.5 both $\mathcal{D}_f(\mathcal{G})$ and $\mathcal{D}_e(\mathcal{G})$ are separable Dirac structures. Note that $\mathcal{D}_f(\mathcal{G})$ and $\mathcal{D}_e(\mathcal{G})$ differ only in the role of the boundary flows and efforts, and that $\mathcal{D}_f(\mathcal{G}) = \mathcal{D}_e(\mathcal{G})$ if there are no boundary vertices. For mass-spring-damper systems the flow-continuous Dirac structure will correspond to the case where the boundary vertices are massless, while the effort-continuous Dirac structure corresponds to boundary masses, with momenta in $\Lambda_{0b}$ .

2.5. Interconnection of open graphs and composition of graph Dirac structures. Interconnection of two open graphs $\mathcal{G}^\alpha$ and $\mathcal{G}^\beta$ is performed by identifying some of their boundary vertices, and equating (up to a minus sign) the boundary efforts and flows corresponding to these boundary vertices, resulting in a new graph. For simplicity of exposition consider the case where the open graphs have all their boundary vertices in common, resulting in a (closed) graph with a set of vertices $\mathcal{V}_i^\alpha \cup \mathcal{V}_i^\beta \cup \mathcal{V}$ , where $\mathcal{V} := \mathcal{V}_b^\alpha = \mathcal{V}_b^\beta$ denotes the set of boundary vertices of both graphs.

The incidence operator of the interconnected (closed) graph is obtained as follows. For simplicity of notation consider the case where $\mathcal{R} = \mathbb{R}$ . Let $\mathcal{G}^j$ have incidence operators

$B^j = \begin{bmatrix} B_i^j \\ B_b^j \end{bmatrix}, \quad j = \alpha, \beta.$

⁵The terminology flow-continuous and effort-continuous stems from the fact that in the first case the boundary flows $f_b$ are exclusively linked to the edge flows $f_1$ , while in the second case the boundary efforts $e^b$ are determined by a part of the internal vertex efforts $e^0$ . Note that the spaces of involved flow and effort variables of $\mathcal{D}_f(\mathcal{G})$ and $\mathcal{D}_e(\mathcal{G})$ are different.

The incidence operator $B$ of the interconnected graph is then given as

$(10) \quad B = \begin{bmatrix} B_i^\alpha & 0 \\ 0 & B_i^\beta \\ B_b^\alpha & B_b^\beta \end{bmatrix},$

corresponding to the interconnection constraints on the boundary potentials and currents given by

$(11) \quad e^{b\alpha} = e^{b\beta}, \quad f_b^\alpha + f_b^\beta = 0.$

Of course, several extensions are possible. For example, one may retain the set of shared boundary vertices $\mathcal{V}_b := \mathcal{V}_b^\alpha = \mathcal{V}_b^\beta$ as being boundary vertices (instead of internal vertices as above) by extending (11) to

$(12) \quad e^{b\alpha} = e^{b\beta} = e^b, \quad f_b^\alpha + f_b^\beta + f_b = 0,$

with $f_b, e^b$ the boundary flows and efforts of the interconnected graph.

Comparing the interconnection of open graphs with the composition of their graph Dirac structures (see, e.g., Proposition 2.6) it is readily seen that the flow/effort-continuous graph Dirac structure of an interconnected graph equals the composition of the flow/effort-continuous graph Dirac structures of $\mathcal{G}^\alpha$ and $\mathcal{G}^\beta$ ; we leave the straightforward proof to the reader.

2.6. Derived graph Dirac structures. Other Dirac structures can be derived from the flow/effort-continuous Dirac structure by constraining some of the flows and the efforts to zero. For example, the composition of the flow/effort-continuous Dirac structure with the trivial separable Dirac structure

$\{(f_{0i}, e^{0i}) \in \Lambda_{0i} \times \Lambda^{0i} \mid f_{0i} = 0\}$

will result by Proposition 2.6 in another separable Dirac structure called the Kirchhoff–Dirac structure, which will be discussed in detail in section 6.

However, there are other possibilities which we will only indicate. One, somewhat dual to the Kirchhoff–Dirac structure, is to constrain (some of) the edge efforts in the flow/effort-continuous graph Dirac structure to zero. Another interesting option is to constrain some of the edge flows in the flow/effort-continuous graph Dirac structure to zero. Considering the description of the flow/effort-continuous graph Dirac structure, this effectively reduces (by disregarding the associated edge efforts) to the flow/effort-continuous graph Dirac structure of the reduced graph where the edges corresponding to the zero edge flows have been left out. Alternatively, one may constrain some of the internal vertex efforts to zero. Again considering the description of the flow/effort-continuous graph Dirac structure, this amounts to deleting the corresponding internal vertices, turning them into boundary vertices with prescribed zero efforts. Note that this yields a setting for dealing with dynamic graphs.

3. Port-Hamiltonian systems on graphs. First (section 3.1) we will describe how port-Hamiltonian systems can be defined with respect to the canonical graph Dirac structures defined above. In the subsequent subsections this will be illustrated on a number of typical examples, ranging from mass-spring-damper systems and spatial mechanisms to consensus and clustering algorithms.

3.1. Definition of port-Hamiltonian systems with regard to the graph Dirac structures. In this subsection we will apply the general definition of port-Hamiltonian systems with regard to an arbitrary Dirac structure (see, e.g., [36, 10, 32]) to the graph Dirac structures as defined above.

For clarity of exposition throughout we consider the effort-continuous graph Dirac structure $\mathcal{D}_e(\mathcal{G})$ involving the flow and effort variables

$(f_1, e^1, f_0, e^0, f_b, e^b) \in \Lambda_1 \times \Lambda^1 \times \Lambda_0 \times \Lambda^0 \times \Lambda_b \times \Lambda^b$

(the exposition is directly repeated for the flow-continuous graph Dirac structure $\mathcal{D}_f(\mathcal{G})$ ). A port-Hamiltonian system is specified by defining, between all the internal conjugate flow and effort variables $(f_1, e^1, f_0, e^0)$ , either an energy-storing relation or a purely dissipative relation. An energy-storing relation between a vector of flow variables $f$ and a conjugate vector of effort variables $e$ is of the form⁶

$\dot{x} = -f, \quad e = \frac{\partial H}{\partial x}(x),$

or dually,

$\dot{x} = e, \quad f = -\frac{\partial H}{\partial x}(x),$

where $x$ is a vector of energy variables (of the same dimension as $f$ and $e$ ), and $H(x)$ is any function, representing the energy stored in the system.

Furthermore, a dissipative relation between a vector of flow variables $f$ and a conjugate vector of effort variables $e$ is any static relation

$R(f, e) = 0$

with the property that $\langle e | -f \rangle \geq 0$ for all $(f, e)$ satisfying $R(f, e) = 0$ .

In the case of a mass-spring-damper system with boundary masses (see subsection 3.2) the vertex flow and effort variables $f_0, e^0$ will be related by energy-storing relations $\dot{p} = -f_0, e_0 = \frac{\partial K}{\partial p}(p)$ , with $p$ the momenta of the masses and $K(p)$ their kinetic energies; the flow and effort variables $f_{1s}, e^{1s}$ of the spring edges will correspond to energy-storing relations $\dot{q} = e^{1s}, f_{1s} = -\frac{\partial V}{\partial q}(q)$ , with $q$ the spring elongations and $V(q)$ the spring potential energies; while finally the flow and effort variables $f_{1d}, e^{1d}$ of the damper edges are connected by energy-dissipating relations $f_{1d} = -D(e^{1d})$ satisfying $(e^{1d})^T D(e^{1d}) \geq 0$ .

Thus a port-Hamiltonian system on a graph is defined by adding to the linear relations imposed by the graph Dirac structure constitutive relations between all the internal effort and flow variables, either of energy-storing or of energy-dissipating type.⁷ It is clear that this leaves many possibilities for defining port-Hamiltonian dynamics. In particular, energy-storage, respectively, energy-dissipation, can be associated to the vertices or the edges or both. The examples presented in the next subsections cover a number of these different possibilities.

⁶Throughout this paper $\frac{\partial H}{\partial x}(x)$ will denote the column vector of partial derivatives of $H$ , with $\frac{\partial^T H}{\partial x}(x)$ denoting the row vector of partial derivatives.

⁷Hence port-Hamiltonian dynamics generalizes both classical Hamiltonian dynamics (with no energy-dissipation), and gradient systems (where there is in general no oscillation between different energies, and energy-dissipation does take place); see [34] and the references therein.

The interpretation of the flow/effort-continuous graph Dirac structure as describing discrete conservation or balance laws becomes clearer from the above description of port-Hamiltonian dynamics. For example, consider for the effort-continuous graph Dirac structure the case of energy-storage associated to all the edges and vertices:

$\begin{aligned} \dot{x}^1 &= e^1, & f_1 &= -\frac{\partial H}{\partial x^1}(x^1, x_0), \\ \dot{x}_0 &= -f_0, & e^0 &= \frac{\partial H}{\partial x_0}(x^1, x_0) \end{aligned}$

for state variables $x^1 \in \Lambda^1$ and $x_0 \in \Lambda_0$ , and energy function $H$ . Then the relations imposed by the effort-continuous graph Dirac structure imply

$\dot{x}_0 + B_i f_1 = 0, \quad \dot{x}^1 + B^* e^0 = 0,$

expressing discrete conservation (or balance) laws between the storage of the quantities $x_0$ associated to the vertices and the flow $f_1$ through the edges, respectively, between the storage of the quantities $x^1$ associated to the edges and the effort $e^0$ at the vertices. The mass-spring system discussed in the next subsection will be of this type.

Furthermore, it is well known [36, 10, 32] that port-Hamiltonian systems may easily entail algebraic constraints on their state variables. Indeed, whenever some of the effort variables $e = \frac{\partial H}{\partial x}(x)$ or $f = -\frac{\partial H}{\partial x}(x)$ are constrained by the Dirac structure, this will generally lead (depending on the Hamiltonian $H$ ) to algebraic constraints on the state variables $x$ .

Finally, we note a fundamental property of any port-Hamiltonian dynamics. Let $H(x)$ denote the total energy of the port-Hamiltonian system. Then because of the power-conserving property of the Dirac structure, and denoting the flows and efforts of the dissipative elements by $f_R, e^R$ ,

$(13) \quad \frac{d}{dt}H(x) = \left\langle \frac{\partial^T H}{\partial x}(x) \mid \dot{x} \right\rangle = \langle e^R \mid f_R \rangle + \langle e^b \mid f_b \rangle \leq \langle e^b \mid f_b \rangle.$

Hence the total energy itself satisfies a conservation law: its increase is equal to the externally supplied power $\langle e^b \mid f_b \rangle$ minus the dissipated power $-\langle e^R \mid f_R \rangle$ .

Remark 3.1. One may directly extend the definition of port-Hamiltonian systems on graphs to the case where the graphs are dynamically changing in time, as briefly indicated in section 2.6. This leads to switching port-Hamiltonian systems on graphs; see [15, 35, 14].

3.2. Mass-spring-damper systems. The basic way of modeling a mass-spring-damper system as a port-Hamiltonian system on a graph is to associate the masses to the vertices, and the springs and dampers to the edges of the graph; cf. Figures 1 and 2. For clarity of exposition we will start with separate treatment of mass-spring (section 3.2.1) and mass-damper (section 3.2.2) systems, before their merging in section 3.2.3.

3.2.1. Mass-spring systems. Consider a graph $\mathcal{G}$ with $N$ vertices (masses) and $M$ edges (springs), specified by an incidence operator $B$ . First, consider the situation where the mass-spring system is located in one-dimensional space $\mathcal{R} = \mathbb{R}$ , and the springs are scalar. A vector in the vertex space $\Lambda_0$ then corresponds to the vector $p$ of the scalar momenta of all $N$ masses, i.e., $p \in \Lambda_0 = \mathbb{R}^N$ . Furthermore, a vector in the dual edge space $\Lambda^1$ will correspond to the total vector $q$ of elongations of all $M$ springs, i.e., $q \in \Lambda^1 = \mathbb{R}^M$ . The next ingredient is the definition of the Hamiltonian (stored energy) $H : \Lambda^1 \times \Lambda_0 \rightarrow \mathbb{R}$ (which normally splits into a sum of the kinetic and potential energies of each mass and spring). In the absence of boundary vertices the dynamics of the mass-spring system is then described as the port-Hamiltonian system

$(14) \quad \begin{bmatrix} \dot{q} \\ \dot{p} \end{bmatrix} = \begin{bmatrix} 0 & B^T \\ -B & 0 \end{bmatrix} \begin{bmatrix} \frac{\partial H}{\partial q}(q, p) \\ \frac{\partial H}{\partial p}(q, p) \end{bmatrix}$

defined with respect to the graph Dirac structure $\mathcal{D}_e(\mathcal{G}) = \mathcal{D}_f(\mathcal{G})$ . Note that in fact the skew-symmetric matrix

$(15) \quad J := \begin{bmatrix} 0 & B^T \\ -B & 0 \end{bmatrix}$

defines a Poisson structure on the state space $\Lambda^1 \times \Lambda_0$ .

The inclusion of boundary vertices, and thereby of external interaction, can be done in different ways. The first option is to associate boundary masses to the boundary vertices. Considering the effort-continuous graph Dirac structure $\mathcal{D}_e(\mathcal{G})$ , we are then led to the port-Hamiltonian system

$(16) \quad \begin{aligned} \dot{q} &= B^T \frac{\partial H}{\partial p}(q, p), \\ \dot{p} &= -B \frac{\partial H}{\partial q}(q, p) + E f_b, \\ e^b &= E^T \frac{\partial H}{\partial p}(q, p). \end{aligned}$

Here $E$ is a matrix with as many columns as there are boundary vertices; each column consists of zeros except for exactly one 1 in the row corresponding to the associated boundary vertex. $f_b \in \Lambda_b$ are the external forces exerted (by the environment) on the boundary masses, and $e_b \in \Lambda^b$ are the velocities of these boundary masses.

Another possibility is to start from the flow-continuous graph Dirac structure $\mathcal{D}_f(\mathcal{G})$ . In this case there are no masses associated to the boundary vertices, and we obtain the port-Hamiltonian system (with $p$ now denoting the vector of momenta of the masses associated to the internal vertices)

$(17) \quad \begin{aligned} \dot{q} &= B_i^T \frac{\partial H}{\partial p}(q, p) + B_b^T e^b, \\ \dot{p} &= -B_i \frac{\partial H}{\partial q}(q, p), \\ f_b &= B_b \frac{\partial H}{\partial q}(q, p), \end{aligned}$

with $e^b \in \Lambda^b$ the velocities of the massless boundary vertices and $f_b \in \Lambda_b$ the forces at the boundary vertices as experienced by the environment. Note that in this latter case the external velocities $e^b$ of the boundary vertices can be considered as inputs to the system and the forces $f_b$ as outputs, in contrast to the previously considered case (boundary vertices corresponding to boundary masses), where the forces $f_b$ are inputs and the velocities $e^b$ are the outputs of the system.⁸

⁸One can also consider the hybrid case where some of the boundary vertices are associated to masses while the remaining ones are massless.

The above formulation of mass-spring systems in $\mathcal{R} = \mathbb{R}$ directly extends to $\mathcal{R} = \mathbb{R}^3$ by using the incidence operator $B = \tilde{B} \otimes I_3$ as defined previously. Finally, we remark that in the above treatment we have considered springs with arbitrary elongation vectors $q \in \Lambda^1$ . For ordinary springs the vector $q$ of elongations is given as $q = B^T q_c$ , where $q_c \in \Lambda^0$ denotes the vector of positions of the vertices. Hence in this case $q \in \text{im } B^T \subset \Lambda^1$ . Note that the subspace $\text{im } B^T \times \Lambda_0 \subset \Lambda^1 \times \Lambda_0$ is an invariant subspace with regard to the dynamics (16) or (17). We will return to this in section 5.1.

3.2.2. Mass-damper systems. Replacing springs by dampers leads to mass-damper systems. In the case of the flow-continuous graph Dirac structure this yields the equations

$\begin{aligned} B_i f_1 &= -\dot{p}, \\ B_b f_1 &= f_b, \\ e^1 &= -B_i^T \frac{\partial H}{\partial p}(p) - B_b^T e^b, \end{aligned} \quad (18)$

where $f_1, e^1$ are the flows and efforts corresponding to the dampers (damping forces, respectively, velocities). For example, for linear dampers $f_1 = -R e^1$ , where $R$ is the positive diagonal matrix with the damping constants on its diagonal. Substitution into (18) then yields the port-Hamiltonian system

$\begin{aligned} \dot{p} &= -B_i R B_i^T \frac{\partial H}{\partial p}(p) - B_i R B_b^T e^b, \\ f_b &= B_b R B_i^T \frac{\partial H}{\partial p}(p) + B_b^T R B_b^T e^b, \end{aligned} \quad (19)$

where, as before, the inputs $e^b$ are the boundary velocities and $f_b$ are the forces as experienced at the massless boundary vertices.

3.2.3. Mass-spring-damper systems. For a mass-spring-damper system the edges will correspond partly to springs and partly to dampers. Thus a mass-spring-damper system is described by a graph $\mathcal{G}(\mathcal{V}, \mathcal{E}_s \cup \mathcal{E}_d)$ , where the vertices in $\mathcal{V}$ correspond to the masses, the edges in $\mathcal{E}_s$ to the springs, and the edges in $\mathcal{E}_d$ to the dampers of the system. This corresponds to an incidence matrix $B = [B_s \quad B_d]$ , where the columns of $B_s$ reflect the spring edges and the columns of $B_d$ reflect the damper edges. For the case without boundary vertices the dynamics of such a mass-spring-damper system with linear dampers takes the form

$\begin{bmatrix} \dot{q} \\ \dot{p} \end{bmatrix} = \begin{bmatrix} 0 & B_s^T \\ -B_s & -B_d R B_d^T \end{bmatrix} \begin{bmatrix} \frac{\partial H}{\partial q}(q, p) \\ \frac{\partial H}{\partial p}(q, p) \end{bmatrix}. \quad (20)$

In the presence of boundary vertices we may distinguish, as above, between massless boundary vertices, with inputs being the boundary velocities and outputs being the boundary (reaction) forces, and boundary masses, in which case the inputs are the external forces and the outputs are the velocities of the boundary masses. We leave the details to the reader.

3.3. Spatial mechanisms. In this section we briefly discuss the extension of mass-spring-damper systems in $\mathbb{R}$ or $\mathbb{R}^3$ to spatial mechanisms, that is, networks of rigid bodies in $\mathbb{R}^3$ related by joints. In this case, the linear space $\mathcal{R}$ is given by $\mathcal{R} := se^*(3)$ , the dual of the Lie algebra of the Lie group $SE(3)$ describing the position of a rigid body in $\mathbb{R}^3$ . A spatial mechanism (or multibody system) is a mechanical system consisting of rigid bodies related by joints (defined as kinematic pairs) restricting the relative motion between the rigid bodies. The reader may find numerous references about their definition and analysis in [2, 30], using different geometric representations of rigid body displacements. In this paper, however, we shall follow the exposition in, e.g., [24, 17], which is based on the Lie group of isometries in Euclidean space $\mathbb{R}^3$ .

The basic topology of the mechanism is described by a directed graph, called the primary graph, whose vertices correspond to the rigid bodies and whose edges are associated with the kinematic pairs. This is similar to the mass-spring or mass-damper systems described in section 3.2, with the differences that the dynamical system associated with each vertex corresponds to rigid body dynamics instead of point-mass dynamics, and that the edges are, in the first instance, associated with kinematic constraints between the bodies instead of springs or dampers. We shall see how (spatial) springs may be included in the second instance.

3.3.1. The rigid body element. The configuration space of a rigid body is the Lie group of isometries in Euclidean space $\mathbb{R}^3$ , called the special Euclidean group and denoted by $SE(3) \ni Q$ (also called the space of rigid body displacements). Using the momentum map associated with the action of $SE(3)$ on its cotangent bundle $T^*SE(3)$ , following for instance [18, Chap. 4], one may define the state space of the rigid body as $SE(3) \times se^*(3) \ni (Q, P)$ by means of the left trivialization, where $P$ is called the momentum in body frame.

The kinetic energy of a rigid body is defined by

$(21) \quad K(P) = \frac{1}{2} \left\langle P, \left( I^b \right)^{-1} (P) \right\rangle,$

where $I^b : se(3) \rightarrow se^*(3)$ is a symmetric, positive-symmetric isomorphism, called the inertia operator of the rigid body in the body frame. The potential energy of the rigid body is defined by a function $U(Q)$ of the displacement $Q$ . The potential energy may be due to gravity or may be zero in the case of the Euler–Poincaré problem.

We assume that the rigid body is subject to an external force expressed as an element $W_e \in se^*(3)$ , called force in fixed frame [18] or wrench in fixed frame [17], which is obtained by the right trivialization of $T^*SE(3)$ . We shall associate a conjugate velocity to this external force, the velocity of the body $T_e$ in fixed frame [18] (also called twist in fixed frame [17]), and obtained by the right trivialization of $TSE(3)$ .

The dynamical equations of the rigid body elements may then be written as a port-Hamiltonian system [36], [20, eq. (1.37)]:

$(22) \quad \begin{aligned} \frac{d}{dt} \begin{pmatrix} Q \\ P \end{pmatrix} &= \begin{pmatrix} 0 & TL_Q \\ -T^*L_Q & -P^\times \end{pmatrix} \begin{pmatrix} dU(Q) \\ (I^b)^{-1}(P) \end{pmatrix} + \begin{pmatrix} 0 \\ Ad_Q^* \end{pmatrix} W_e, \\ T_e &= \begin{pmatrix} 0 & Ad_Q \end{pmatrix} \begin{pmatrix} dU(Q) \\ (I^b)^{-1}(P) \end{pmatrix}, \end{aligned}$

where $TL_Q$ denotes the tangent map to the left translation (mapping the velocities $T \in se(3)$ in body frame into the velocities $v \in T_Q SE(3)$ ), $T^*L_Q$ denotes its dual map (mapping forces $F \in T_Q^*SE(3)$ into forces in body frame $W \in se^*(3)$ ), $Ad_Q$ denotes the adjoint representation (mapping velocities in body frame into velocities in fixed frame), $Ad_Q^*$ denotes the adjoint map to $Ad_Q$ , while finally $\times$ is defined by the coadjoint representation of the Lie algebra $se(3)$ , that is, $W \times T = ad_T^*W$ , for any $(W, T) \in se^*(3) \times se(3)$ . The Dirac structure $\mathcal{D}_{RB}$ of this port-Hamiltonian system (22) is thus specified as⁹

$(23) \quad \begin{aligned} \mathcal{D}_{RB}(Q) = & \left\{ (v, W, T_e, F, T, W_e) \right. \\ & \left. \in T_QSE(3) \times se^*(3) \times se(3) \times T_Q^*SE(3) \times se(3) \times se^*(3) \mid \right. \\ & \left( \begin{array}{c} v \\ W \end{array} \right) = \left( \begin{array}{cc} 0 & TL_Q \\ -T^*L_Q & -P \times \end{array} \right) \left( \begin{array}{c} F \\ T \end{array} \right) + \left( \begin{array}{c} 0 \\ Ad_Q^* \end{array} \right) W_e, \\ & T_e = \left( \begin{array}{cc} 0 & Ad_Q \end{array} \right) \left( \begin{array}{c} F \\ T \end{array} \right) \left. \right\}. \end{aligned}$

In this way we have associated with every vertex of the primary graph of the spatial mechanism a dynamical system (22) with inputs and outputs $(W_e, T_e) \in se^*(3) \times se(3)$ .

3.3.2. The kinematic pair. Constraints between the rigid bodies of the mechanism will be specified by kinematic pairs corresponding to each edge of the primary graph. A kinematic pair is the idealization of a set of contacts that occur between two rigid bodies at some configuration of the bodies. It constrains the possible relative twists between the bodies as well as the possible transmitted wrenches. The wrench $W$ transmitted by a kinematic pair is constrained to a linear subspace of the space of wrenches $se^*(3)$ called the space of constraint wrenches and denoted by $CW$ . A relative twist between the two bodies is allowed by the kinematic pair when it produces no work with any transmissible wrench. The relative twist is thus constrained to a linear subspace $\mathcal{FT}$ of the space of twists $se(3)$ , called the space of freedom twists. Since an ideal kinematic pair is workless, the subspace $\mathcal{FT}$ is orthogonal (in the sense of the duality product) to the space of transmitted wrenches $CW$ , that is $\mathcal{FT} = CW^\perp$ .

We have defined the spaces of freedom twists and constraint wrenches as subspaces of the Lie algebra $se(3)$ and its dual. However, these spaces express constraints on the twists and wrenches of the rigid bodies related by the kinematic pairs, and hence are expressed in some common frame with configuration $Q_{KP}$ (in most cases equal to the configuration of one of the related bodies). Consequently, the constitutive relations of a kinematic pair are given in terms of its pair of twists and wrenches $(T_{KP}, W_{KP}) \in T_{Q_{KP}}SE(3) \times T_{Q_{KP}}^*SE(3)$ in the form

$(24) \quad Ad_{Q_{KP}}^* W_{KP} \in CW \quad \text{and} \quad Ad_{Q_{KP}^{-1}} T_{KP} \in \mathcal{FT}.$

Hence the constitutive equations of a kinematic pair may be expressed as the following nonconstant separable Dirac structure:

$(25) \quad \mathcal{D}_{CW}(Q_{KP}) = \left\{ (T_{KP}, W_{KP}) \in T_{Q_{KP}}SE(3) \times T_{Q_{KP}}^*SE(3) \mid \right. \\ \left. Ad_{Q_{KP}}^* W_{KP} \in CW \quad \text{and} \quad Ad_{Q_{KP}^{-1}} T_{KP} \in CW^\perp \right\}.$

⁹Note that this is a nonconstant Dirac structure on $SE(3)$ .

The kinematic pair introduced above represents ideal kinematic constraints; in general, however, mechanical work may be produced at the kinematic pair due to the presence of actuators or springs and dampers. Such an interaction is captured by considering the linear space $\mathcal{IW} := se^*(3)/\mathcal{CW}$ (which may be identified with a subspace of $se^*(3)$ complementary to the space of constraint wrenches $\mathcal{CW}$ ). The space of interaction twists is then defined as its dual space $\mathcal{IT} := \mathcal{IW}^* \simeq \mathcal{CW}^\perp$ . Using the canonical projection $\pi$ of $se^*(3)$ onto $\mathcal{IW}$ , together with its adjoint map $\pi^*$ , one may thus define an additional pair of port variables that are able to connect actuators, or damper or spring elements to the kinematic pairs. The resulting interacting kinematic pair is then defined as a 2-port element with constitutive relations defined by the following nonconstant separable Dirac structure:

$\begin{aligned}
\mathcal{D}_{\mathcal{CW}}^I(Q_{KP}) &= \{(T_{KP}, W_{KP}, T_I, W_I) \\
(26) \quad &\in T_{Q_{KP}}SE(3) \times T_{Q_{KP}}^*SE(3) \times \mathcal{CW}^\perp \times se^*(3)/\mathcal{CW} \mid \\
&W_I = \pi \circ Ad_{Q_{KP}}^*(W_{KP}), T_{KP} = -Ad_{Q_{KP}} \circ \pi^*(T_I)\}.
\end{aligned}$

It is easy to check that for $W_I = 0$ the interacting kinematic pair reduces to the kinematic pair as defined previously.

3.3.3. The kinestatic connection network. The primary graph of the mechanism together with the kinematic pairs is called the kinestatic model of the mechanical system. Its associated Dirac structure is the composition of the Dirac structures corresponding to the kinematic pairs with the flow-continuous¹⁰ graph Dirac structure of the primary graph.

Consider a mechanism defined by its primary graph $\mathcal{G}$ composed of $n_{RB}$ internal vertices (associated with the rigid bodies), $n_b$ boundary vertices corresponding to rigid bodies with zero inertia operator, and $n_{KP}$ edges (associated with the kinematic pairs). Define the vertex space $\Lambda_0 \ni T^{RB}$ and the edge space $\Lambda_1 \ni T^{KP}$ with respect to the Lie algebra $se(3)$ , which represent, respectively, the external twist of the rigid bodies and the kinematic pairs. The dual spaces $\Lambda^0 \ni W^{RB}$ , respectively, $\Lambda^1 \ni W^{KP}$ , then represent the external wrenches of the rigid bodies, respectively, the wrenches of the kinematic pairs. The twists and wrenches of the boundary vertices (the rigid bodies with zero inertia operator) are associated with the vertex space $\Lambda_b \ni T^b$ , respectively, its dual $\Lambda^b \ni W^b$ . Kirchhoff's laws on the twists and wrenches [11] amount to constraining these variables to belong to the flow-continuous graph Dirac structure, i.e.,

$(T^{KP}, W^{KP}, T^{RB}, W^{RB}, T^b, W^b) \in \mathcal{D}_f(\mathcal{G}).$

Composition of $\mathcal{D}_f(\mathcal{G})$ with the Dirac structures $\mathcal{D}_{\mathcal{CW}}(Q_{KP})$ corresponding to all the kinematic pairs then results in the Dirac structure $\mathcal{D}_{KS}$ of the kinestatic model:

$(27) \quad (T^I, W^I, T^{RB}, W^{RB}, T^b, W^b) \in \mathcal{D}_{KS}.$

3.3.4. Dynamics of spatial mechanisms. The state space $\mathcal{X}$ of the complete mechanism is the product space of the state spaces of all the rigid bodies, i.e., $\mathcal{X} = (SE(3) \times se^*(3))^{n_{RB}}$ , where $n_B$ denotes the number of rigid bodies (equal to the

¹⁰Or the effort-continuous graph Dirac structure in case the rigid bodies corresponding to the boundary vertices have nonzero inertia operator.

number of internal vertices of the primary graph). Recalling that the rigid body dynamics is defined as a port-Hamiltonian system with respect to the Dirac structure (23), one then obtains the overall Dirac structure $\mathcal{D}_M$ of the mechanism by composing the Dirac structure $\mathcal{D}_{KS}$ of the kinematic model with the Dirac structures $\mathcal{D}_{RB}$ of all the rigid bodies. Finally, defining the Hamiltonian $H_M(x)$ as the sum of the Hamiltonians of each body, one obtains the following port-Hamiltonian model of the mechanism:

$(28) \quad \left( -\frac{dx}{dt}, \frac{\partial H_M}{\partial x}(x), T^I, W^I, T^b, W^b \right) \in \mathcal{D}_M.$

3.4. Hydraulic networks. The interpretation of the flow/effort-continuous graph Dirac structures as capturing the basic conservation/balance laws of a network becomes especially tangible for hydraulic networks.

A hydraulic network can be modeled as a directed graph with edges corresponding to pipes; see, e.g., [29, 12]. The vertices may correspond to either connection points with fluid reservoirs (buffers) or merely connection points of the pipes; we concentrate on the first case (the second case corresponding to a Kirchhoff–Dirac structure; cf. section 6.1). Let $x_v$ be the stored fluid at vertex $v$ and let $\nu_e$ be the fluid flow through edge $e$ . Collecting all stored fluids $x_v$ into one vector $x$ , and all fluid flows $\nu_e$ into one vector $\nu$ , the mass-balance is summarized as

$(29) \quad \dot{x} = -B\nu$

with $B$ denoting the incidence matrix of the graph. In the absence of fluid reservoirs this simply reduces to Kirchhoff’s current laws $B\nu = 0$ .

For incompressible fluids a standard model of the fluid flow $\nu_e$ through pipe $e$ is

$(30) \quad J_e \dot{\nu}_e = P_i - P_j - \lambda_e(\nu_e),$

where $P_i$ and $P_j$ are the pressures at the tail, respectively, head, vertices of edge $e$ . Note that this in fact captures two effects: one corresponding to energy storage and one corresponding to energy dissipation. First, defining the energy variable $\varphi_e := J_e \nu_e$ , the stored energy in the pipe associated with edge $e$ is given as $\frac{1}{2J_e} \varphi_e^2 = \frac{1}{2} J_e \nu_e^2$ . Second, $\lambda_e(\nu_e)$ is a damping force corresponding to energy dissipation.

In the case of fluid reservoirs at the vertices the pressures $P_v$ at each vertex $v$ are functions of $x_v$ , and thus, being scalar functions, are always derivable from an energy function $P_v = \frac{\partial H_v}{\partial x_v}(x_v)$ , $v \in \mathcal{V}$ , for some Hamiltonian $H_v(x_v)$ (e.g., gravitational energy). The resulting dynamics (with state variables $x_v$ and $\varphi_e$ ) is port-Hamiltonian with respect to the graph Dirac structure $\mathcal{D}_f(\mathcal{G}) = \mathcal{D}_e(\mathcal{G})$ . The setup is immediately extended to boundary vertices (corresponding to either controlled fluid reservoirs or direct in/out-flows).

3.5. Port-Hamiltonian formulation of consensus algorithms. While all previous examples of port-Hamiltonian systems on graphs arise from physical modeling, the system treated in this subsection has a different origin. Nevertheless, it shares the same structure and, in fact, turns out to have dynamics equal to the mass-damper system treated previously.

Consider a network of $N$ agents moving in linear space $\mathcal{R}$ , whose interaction topology is described by an undirected graph $\mathcal{G}$ (symmetric interaction). Denote by $E(\mathcal{G})$ the edges of this undirected graph, consisting of unordered pairs $(v, w)$ of vertices $v, w$ . Hence $(v, w) \in E(\mathcal{G})$ if and only if $(w, v) \in E(\mathcal{G})$ . Thus the vertices of the graph correspond to the agents, and the edges to the symmetric interactions between them. Distinguish between leader and follower agents (see, e.g., [28]), and associate the leader agents to the boundary vertices and the follower agents to the internal vertices.

Associated to each agent $v$ is a vector $x_v \in \mathcal{R}$ describing the motion in the linear space $\mathcal{R}$ . In the standard consensus algorithm (see, e.g., [25]), the vector $x_v$ of each follower agent $v$ satisfies the dynamics

$(31) \quad \dot{x}_v(t) = - \sum_{(v,w) \in E(\mathcal{G})} g_{(v,w)}(x_v(t) - x_w(t)),$

where $g_{(v,w)} > 0$ denotes a certain positive-definite weight matrix associated to each edge. For simplicity of exposition let us take the linear space $\mathcal{R}$ to be equal to $\mathbb{R}$ in the rest of this section, implying that $g_{(v,w)} > 0$ are just positive constants. Collecting all follower variables $x_v$ into one vector $x \in \mathbb{R}^{N_i}$ , and all leader variables $x_v$ into one vector $u \in \mathbb{R}^{N_b}$ , it is readily checked that the dynamics can be written as

$(32) \quad \dot{x} = -B_i G B_i^T x - B_i G B_b^T u,$

with $B$ the incidence matrix of the graph endowed with an arbitrary orientation,¹¹ and $G$ the diagonal matrix with elements $g_{(v,w)}$ corresponding to each edge $(v,w)$ . This defines a port-Hamiltonian system with respect to the flow-continuous graph Dirac structure $\mathcal{D}_f(\mathcal{G})$ and the Hamiltonian $H(x) := \frac{1}{2} \|x\|^2$ . Indeed, (32) is equal to

$(33) \quad \dot{x} = -B_i G B_i^T \frac{\partial H}{\partial x}(x) - B_i G B_b^T u,$

which are the same equations as those for the mass-damper system (19), with $u = e_b \in \Lambda^b$ . Note that the corresponding artificial output vector $y = f_b \in \Lambda_b$ given as

$y := B_b G B_i^T \frac{\partial H}{\partial x}(x) + B_b G B_b^T u$

equals minus the rate of the leader variables if the leader variables were supposed to obey the consensus algorithm with regard to the follower agents (which is not the case). Hence this artificial output measures the discrepancy between the leaders and the followers.

3.5.1. Network clustering dynamical models. A dynamical model for network clustering, where the network splits into subnetworks which separately reach consensus, was recently proposed and discussed in [6]. Consider again a multiagent system with $N$ agents and state variables $x_i \in \mathbb{R}$ , $i = 1, \dots, N$ , whose dynamics is described as

$(34) \quad \dot{x}_i = -\frac{dJ_i}{dx_i}(x_i) + u_i, \quad i = 1, \dots, N,$

where the functions $J_i(x_i)$ are certain objective functions. Let the vector $u$ with components $u_i$ be determined as

$(35) \quad u = B \frac{\partial V}{\partial z}(z), \quad z = -B^T x,$

¹¹It is easily seen [4] that the Laplacian matrix $BGB^T$ is independent of the chosen orientation.

where $V(z) = V_1(z_1) + \dots + V_M(z_M)$ for certain functions $V_1, \dots, V_M$ . This is readily seen to result in a port-Hamiltonian system with total Hamiltonian $H(x, z) = \frac{1}{2}\|x\|^2 + V(z)$ , and a nonlinear resistive characteristic associated to each $i$ th vertex defined by the functions $J_i(x_i)$ , interpreted as Rayleigh dissipation functions.¹² Clustering may occur once the energy functions $V_i$ define bounded constitutive relations $e_{1i} = \frac{dV_i}{dz_i}(z_i)$ for the edge efforts. Depending on the strength of the objective functions $J_i$ , this will imply that consensus among the $x_i$ -variables will only be reached for subnetworks.

Many other models of network dynamics, of a “nonphysical” background, can be formulated as port-Hamiltonian systems on graphs. Examples include coordination control [1] and edge agreement [47].

4. Dynamical analysis. In this section we will investigate the dynamical properties of a paradigmatic example of a port-Hamiltonian system on a graph, namely the mass-spring-damper system as discussed in section 3.2.3. As we have seen, many other examples share the same mathematical structure, and the dynamical analysis will follow along the same lines.

Thus we will consider a mass-spring-damper system as described by a graph $\mathcal{G}(\mathcal{V}, \mathcal{E}_s \cup \mathcal{E}_d)$ , where the vertices in $\mathcal{V}$ correspond to the masses, the edges in $\mathcal{E}_s$ to the springs, and the edges in $\mathcal{E}_d$ to the dampers of the system, with incidence matrix $B = [B_s \ B_d]$ , where the columns of $B_s$ reflect the spring edges and the columns of $B_d$ reflect the damper edges. Without boundary vertices the dynamics takes the form (see (20) in section 3.2.3)

$(36) \quad \begin{bmatrix} \dot{q} \\ \dot{p} \end{bmatrix} = \begin{bmatrix} 0 & B_s^T \\ -B_s & -B_d R B_d^T \end{bmatrix} \begin{bmatrix} \frac{\partial H}{\partial q}(q, p) \\ \frac{\partial H}{\partial p}(q, p) \end{bmatrix}.$

Throughout this section we make the following simplifying assumption.¹³

ASSUMPTION 4.1. The graph $\mathcal{G}(\mathcal{V}, \mathcal{E}_s \cup \mathcal{E}_d)$ is connected or, equivalently, $\ker B_s^T \cap \ker B_d^T = \text{span } \mathbb{1}$ .

1.1. 4.1. Equilibria and Casimirs.

PROPOSITION 4.1. The set of equilibria $\mathcal{E}$ of (36) is given as

$(37) \quad \mathcal{E} = \left\{ (q, p) \in \Lambda^1 \times \Lambda_0 \mid \frac{\partial H}{\partial q}(q, p) \in \ker B_s, \frac{\partial H}{\partial p}(q, p) \in \text{span } \mathbb{1} \right\}.$

Proof. $(q, p)$ is an equilibrium whenever

$B_s^T \frac{\partial H}{\partial p}(q, p) = 0, \quad B_s \frac{\partial H}{\partial q}(q, p) + B_d R B_d^T \frac{\partial H}{\partial p}(q, p) = 0.$

Premultiplication of the second equation by the row vector $\frac{\partial^T H}{\partial p}(q, p)$ , making use of the first equation, yields $\frac{\partial^T H}{\partial p}(q, p) B_d R B_d^T \frac{\partial H}{\partial p}(q, p) = 0$ or, equivalently, $B_d^T \frac{\partial H}{\partial p}(q, p) = 0$ . This in turn implies $B_s \frac{\partial H}{\partial q}(q, p) = 0$ . $\square$

In other words, for $(q, p)$ to be an equilibrium, $\frac{\partial H}{\partial p}(q, p)$ should satisfy the consensus conditions corresponding to the spring-damper graph $\mathcal{G}(\mathcal{V}, \mathcal{E}_s \cup \mathcal{E}_d)$ , whereas

¹²The condition of convexity imposed in [6] on $J_i$ thus corresponds to incremental passivity.

¹³This assumption can be made without loss of generality, since otherwise the same analysis can be repeated for each connected component.

$\frac{\partial H}{\partial q}(q, p)$ should be in the space of cycles of the spring graph $\mathcal{G}(\mathcal{V}, \mathcal{E}_s)$ (corresponding to zero net spring forces applied to the masses at the vertices).

Similarly, the Casimirs (conserved quantities independent of the Hamiltonian $H$ ) are computed as follows.

PROPOSITION 4.2. The Casimir functions are all functions $C(q, p)$ satisfying

$(38) \quad \frac{\partial C}{\partial p}(q, p) \in \text{span } \mathbb{1}, \quad \frac{\partial C}{\partial q}(q, p) \in \ker B_s.$

Proof. $C(q, p)$ is a Casimir if

$\begin{bmatrix} \frac{\partial C}{\partial q}(q, p) & \frac{\partial C}{\partial p}(q, p) \end{bmatrix} \begin{bmatrix} 0 & B_s^T \\ -B_s & -B_d R B_d^T \end{bmatrix} = 0$

or, equivalently,

$\frac{\partial^T C}{\partial p}(q, p) B_s = 0, \quad \frac{\partial^T C}{\partial q}(q, p) B_s^T + \frac{\partial^T C}{\partial p}(q, p) B_d R B_d^T = 0.$

Postmultiplication of the second equation by $\frac{\partial C}{\partial p}(q, p)$ , making use of the first equation, gives the result. $\square$

Therefore all Casimir functions can be expressed as functions of the linear Casimir functions

$(39) \quad C(q, p) = \mathbb{1}^T p, \quad C(q, p) = k^T q, \quad k \in \ker B_s.$

This implies that starting from an arbitrary initial position $(q_0, p_0) \in \Lambda^1 \times \Lambda_0$ , the solution of the mass-spring-damper system (36) will be contained in the affine space

$(40) \quad \mathcal{A}_{(q_0, p_0)} := \begin{bmatrix} q_0 \\ p_0 \end{bmatrix} + \begin{bmatrix} 0 \\ \ker \mathbb{1}^T \end{bmatrix} + \begin{bmatrix} \text{im } B_s^T \\ 0 \end{bmatrix};$

i.e., for all $t$ the difference $q(t) - q_0$ remains in the space of cocycles of the spring graph, while $\mathbb{1}^T p(t) = \mathbb{1}^T p_0$ .

4.2. Stability analysis. Under generic conditions on the Hamiltonian $H(q, p)$ , each affine space $\mathcal{A}_{(q_0, p_0)}$ will intersect the set of equilibria $\mathcal{E}$ in a single point $(q_\infty, p_\infty)$ , which will qualify as the point of asymptotic convergence starting from $(q_0, p_0)$ (provided there is enough damping present). In order to simplify the statement of the results, throughout this subsection we will consider linear mass-spring systems, corresponding to a quadratic and decoupled Hamiltonian function

$(41) \quad H(q, p) = \frac{1}{2} q^T K q + \frac{1}{2} p^T G p,$

where $K$ is the positive diagonal matrix of spring constants, and $G$ is the positive diagonal matrix of reciprocals of the masses. It follows that the set of equilibria is given as $\mathcal{E} = \{(q, p) \in \Lambda^1 \times \Lambda_0 \mid Kq \in \ker B_s, Gp \in \text{span } \mathbb{1}\}$ , while for each $(q_0, p_0)$ there exists a unique point $(q_\infty, p_\infty) \in \mathcal{E} \cap \mathcal{A}_{(q_0, p_0)}$ . In fact, $q_\infty$ is given by the spring graph cocycle/cycle decomposition

$(42) \quad q_0 = v_0 + q_\infty, \quad v_0 \in \text{im } B_s^T \subset \Lambda^1, \quad Kq_\infty \in \ker B_s \subset \Lambda_1,$

while $p_\infty$ is uniquely determined by¹⁴

$(43) \quad Gp_\infty \in \text{span } \mathbb{1}, \quad \mathbb{1}^T p_\infty = \mathbb{1}^T p_0.$

This leads to the following asymptotic stability theorem. First, note that the energy $H(q, p) = \frac{1}{2}q^T Kq + \frac{1}{2}p^T Gp$ (which obviously is radially unbounded) satisfies

$(44) \quad \frac{d}{dt}H(q, p) = -\frac{\partial^T H}{\partial p}(q, p)B_d R B_d^T \frac{\partial H}{\partial p}(q, p) = -p^T B_d G R B_d^T G p \leq 0,$

and thus qualifies as a Lyapunov function, showing at least stability.

THEOREM 4.3. Consider a linear mass-spring-damper system with $H(q, p) = \frac{1}{2}q^T Kq + \frac{1}{2}p^T Gp$ , where $K$ and $G$ are diagonal positive matrices. Then for every $(q_0, p_0)$ there exists a unique equilibrium point $(q_\infty, p_\infty) \in \mathcal{E} \cap \mathcal{A}_{(q_0, p_0)}$ , determined by (42), (43). Define the spring Laplacian matrix $L_s := B_s K B_s^T$ . Then for every $(q_0, p_0)$ the following holds: the trajectory starting from $(q_0, p_0)$ converges asymptotically to $(q_\infty, p_\infty)$ if and only if the largest $GL_s$ -invariant subspace contained in $\ker B_d^T$ is equal to $\text{span } \mathbb{1}$ .

Proof. By Lasalle's invariance principle and (44) the trajectory converges to the largest invariant subspace contained in $\{(q, p) \mid B_d^T G p = 0\}$ . Differentiation of $B_d^T G p = 0$ yields

$0 = \frac{d}{dt}B_d^T G p = B_d^T G(-B_s K q - B_d R B_d^T G p) = -B_d^T G B_s K q,$

while further differentiation gives

$0 = B_d^T G B_s K B_s^T G p = B_d^T G L_s G p.$

By repeated differentiation one thus concludes that $Gp(t)$ for $t \rightarrow \infty$ will converge to the largest $GL_s$ -invariant subspace $\mathcal{V}$ contained in $\ker B_d^T$ , while $q(t)$ will converge to the subspace $\{q \mid GB_s K q \in \mathcal{V}\}$ . Thus if $\mathcal{V} = \text{span } \mathbb{1}$ , then $Gp(t) \rightarrow Gp_\infty$ with $Gp_\infty \in \text{span } \mathbb{1}$ . Furthermore $q(t) \rightarrow q_\infty$ with $GB_s K q_\infty \in \text{span } \mathbb{1} \subset \ker B_s^T$ , and thus $B_s^T GB_s K q_\infty = 0$ or, equivalently, $B_s K q_\infty = 0$ . $\square$

The condition that the largest $GL_s$ -invariant subspace contained in $\ker B_d^T$ is equal to $\text{span } \mathbb{1}$ amounts to pervasive damping: the influence of the dampers spreads through the whole system.

Remark 4.2. Theorem 4.3 is closely related to recent results on partial consensus for double-integrator multiagent systems [16, 7], as will become clear from the discussion in section 5.1.

Another feature of the dynamics of the mass-spring-damper system (36) is its robustness with regard to constant external (disturbance) forces. Indeed, consider a mass-spring-damper system with boundary masses (see section 3.2) and general Hamiltonian $H(q, p)$ , subject to constant forces $\bar{f}_b$ ,

$(45) \quad \begin{bmatrix} \dot{q} \\ \dot{p} \end{bmatrix} = \begin{bmatrix} 0 & B_s^T \\ -B_s & -B_d R B_d^T \end{bmatrix} \begin{bmatrix} \frac{\partial H}{\partial q}(q, p) \\ \frac{\partial H}{\partial p}(q, p) \end{bmatrix} + \begin{bmatrix} 0 \\ E \end{bmatrix} \bar{f}_b,$

¹⁴ $Gp_\infty = c\mathbb{1}$ , where the constant $c$ is determined by the initial value vector $p_0$ via the formula $c \sum_{i=1}^N \prod_{j \neq i} g_j = (\Pi_{i=1}^N g_i) (\sum_{i=1}^N p_{0i})$ .

where we assume¹⁵ the existence of a $\bar{q}$ such that

$(46) \quad B_s \frac{\partial H}{\partial q}(\bar{q}, 0) = E \bar{f}_b.$

Then the availability function

$(47) \quad \bar{H}(q, p) := H(q, p) - (q - \bar{q})^T \frac{\partial H}{\partial q}(\bar{q}, 0) - H(\bar{q}, 0)$

satisfies

$(48) \quad \frac{d}{dt} \bar{H}(q, p) = -\frac{\partial^T H}{\partial p}(q, p) B_d R B_d^T \frac{\partial H}{\partial p}(q, p) \leq 0.$

Specializing to $H(q, p) = \frac{1}{2}q^T K q + \frac{1}{2}p^T G p$ , in which case $\bar{H}(q, p) = \frac{1}{2}(q - \bar{q})^T K (q - \bar{q}) + \frac{1}{2}p^T G p$ , we obtain the following analogue of Theorem 4.3.

PROPOSITION 4.4. Consider a linear mass-spring-damper system (45) with constant external disturbance $\bar{f}_b$ and Hamiltonian $H(q, p) = \frac{1}{2}q^T K q + \frac{1}{2}p^T G p$ , where $K$ and $G$ are diagonal positive matrices, and with $\text{im } E \subset \text{im } B_s$ . The set of controlled equilibria is given by $\mathcal{E} = \{(q, p) \in \Lambda^1 \times \Lambda_0 \mid B_s K q = E \bar{f}_b, G p \in \text{span } \mathbb{1}\}$ . For every $(q_0, p_0)$ there exists a unique equilibrium point $(\bar{q}_\infty, p_\infty) \in \mathcal{E} \cap \mathcal{A}_{(q_0, p_0)}$ . Here $p_\infty$ is determined by (43), while $\bar{q}_\infty = \bar{q} + q_\infty$ , with $\bar{q}$ such that $B_s K \bar{q} = E \bar{f}_b$ and $q_\infty$ the unique solution of (42) with $q_0$ replaced by $q_0 - \bar{q}$ . Furthermore, for each $(q_0, p_0)$ the trajectory starting from $(q_0, p_0)$ converges asymptotically to $(\bar{q}_\infty, p_\infty)$ if and only if the largest $GL_s$ -invariant subspace contained in $\ker B_d^T$ is equal to $\text{span } \mathbb{1}$ .

Note that the above proposition has a classical interpretation in terms of the robustness of integral control with regard to constant disturbances: the springs act as integral controllers which counteract the influence of the unknown external force $\bar{f}_b$ so that the vector of momenta $p$ will still converge to consensus.¹⁶

Remark 4.3. The analysis for a mass-damper systems with constant external velocities (19) and, equivalently, for the leader-follower network (33) is somewhat different. Assuming the graph to be connected, it is well known [4] that for each vector $\bar{e}^b$ there exists a unique equilibrium vector $\bar{p}$ such that

$(49) \quad 0 = -B_i R B_i^T \bar{p} - B_i R B_b^T \bar{e}^b.$

Asymptotic stability of $\bar{p}$ can then be proved by defining the availability function $\bar{H}(p) := H(p) - (p - \bar{p})^T \frac{\partial H}{\partial p}(\bar{p}) - H(\bar{p})$ , satisfying

$(50) \quad \begin{aligned} \frac{d}{dt} \bar{H}(p) &= -[(p - \bar{p})^T \quad (e^b - \bar{e}^b)^T] B G B^T \begin{bmatrix} p - \bar{p} \\ e^b - \bar{e}^b \end{bmatrix} \\ &\quad + (f_b - \bar{f}_b)^T (e^b - \bar{e}^b), \end{aligned}$

where $\bar{f}_b = B_b G B_b^T \bar{p} + B_b G B_b^T \bar{e}^b$ is the output equilibrium value. Since the set $\{p \mid (p - \bar{p})^T B_i G B_i^T (p - \bar{p}) = 0\}$ is equal to the single point $\bar{p}$ (because $\ker B G B^T = \text{span } \mathbb{1}$ ), this shows asymptotic stability of the controlled equilibrium $\bar{p}$ for $e^b = \bar{e}^b$ .

¹⁵If the mapping $q \rightarrow \frac{\partial H}{\partial q}(q, 0)$ is surjective, then there exists for every $\bar{f}_b$ such a $\bar{q}$ if and only if $\text{im } E \subset \text{im } B_s$ .

¹⁶The above proposition can be also applied to leader-follower networks (section 3.5.1), implying that constant leader inputs can be counteracted by integral control (on top of the standard consensus algorithm).

5. Port-Hamiltonian systems on graphs obtained by symmetry reduction. In this section we will show how port-Hamiltonian systems on graphs, such as the mass-spring-damper systems, can be alternatively obtained by symmetry reduction from a symplectic formulation, exploiting the invariance of the Hamiltonian function (in particular, of the spring potential energies).

5.1. Symmetry reduction from the symplectic formulation. Let us return to the formulation of a mass-spring system in section 3.2, where the vertices correspond to the masses, and the edges to the springs in between them. An alternative is to consider the configuration¹⁷ vector $q_c \in \Lambda^0 =: Q_c$ , describing the positions of all the masses. In fact, this is the classical starting point for Lagrangian mechanics, where we do not start with the energy variables $q$ and $p$ , but instead we start with the configuration vector $q_c$ and the corresponding velocity vector $\dot{q}_c$ . The classical Hamiltonian formulation is then obtained by defining the vector of momenta $p \in \Lambda_0 = Q_c^*$ as $p = M\dot{q}_c$ (with $M$ the diagonal mass matrix), resulting in the symplectic phase space $Q_c \times Q_c^* = \Lambda^0 \times \Lambda_0$ . For ordinary springs the relation between $q_c \in \Lambda^0$ and the vector $q \in \Lambda^1$ describing the elongations of the springs is given as $q = B^T q_c$ . Hence in this case the Hamiltonian can also be expressed as a function $H_c$ of $(q_c, p)$ by defining

$(51) \quad H_c(q_c, p) := H(B^T q_c, p).$

It follows that the equations of motion of the mass-spring system (with boundary masses) are given by the canonical Hamiltonian equations

$(52) \quad \begin{aligned} \dot{q}_c &= \frac{\partial H_c}{\partial p}(q_c, p), \\ \dot{p} &= -\frac{\partial H_c}{\partial q_c}(q_c, p) + E f_b, \\ e^b &= E^T \frac{\partial H_c}{\partial p}(q_c, p), \end{aligned}$

where, as before, $f_b$ are the external forces exerted on the boundary masses and $e_b$ are their velocities.

What is the relation of (52) with the port-Hamiltonian formulation given in section 3.2? It turns out that this relation is precisely given by the standard procedure of symmetry reduction of a Hamiltonian system.¹⁸ Indeed, since $\mathbb{1}^T B = 0$ , the Hamiltonian function $H_c(q_c, p)$ given in (51) is invariant under the action of the group $\mathfrak{G} = \mathbb{R}$ acting on the phase space $\Lambda^0 \times \Lambda_0 \simeq \mathbb{R}^{2N}$ by the symplectic group action

$(53) \quad (q_c, p) \mapsto (q_c + \alpha \mathbb{1}, p), \quad \alpha \in \mathfrak{G} = \mathbb{R}.$

From standard reduction theory (see, e.g., [19, 18] and the references therein), it follows that we may factor out the configuration space $Q_c := \Lambda^0$ to the reduced configuration space

$(54) \quad Q := \Lambda^0 / \mathfrak{G}.$

¹⁷Note that $q_c \in \Lambda^0$ is defined to be a function $q_c : \mathcal{V} \rightarrow \mathbb{R}$ , assigning to each vertex its position in $\mathcal{R}$ .

¹⁸This relation can be regarded as the discrete, graph-theoretic version of the correspondence between the port-Hamiltonian formulation of the Maxwell equations (using the Stokes-Dirac structure) and its symplectic formulation using the vector potential of the magnetic field; cf. [19, 44].

Let us first assume that the graph is connected, in which case (see, e.g., [4]) $\ker B^T = \text{span } \mathbf{1}$ . Then we have the following identification:

$(55) \quad Q := \Lambda^0 / \mathfrak{G} \simeq B^T \Lambda^0 \subset \Lambda^1.$

Hence the reduced state space of the mass-spring system is given by $\text{im } B^T \times \Lambda_0$ , where $\text{im } B^T \subset \Lambda^1$ . Furthermore, under the symmetry action, the canonical Hamiltonian equations (52) on the symplectic space $\Lambda^0 \times \Lambda_0$ reduce to the port-Hamiltonian equations (16) on $\text{im } B^T \times \Lambda_0 \subset \Lambda^1 \times \Lambda_0$ obtained previously:

$(56) \quad \begin{aligned} \dot{q} &= B^T \dot{q}_c = B^T \frac{\partial H_c}{\partial p}(q_c, p) = B^T \frac{\partial H}{\partial p}(q, p), \\ \dot{p} &= -\frac{\partial H_c}{\partial q_c}(q_c, p) + E f_b = -B \frac{\partial H}{\partial q}(q, p) + E f_b, \\ e^b &= E^T \frac{\partial H}{\partial p}(q, p). \end{aligned}$

In the case when the graph is not connected, the above symmetry reduction can be performed for each component of the graph (i.e., the symmetry group is $\mathbb{R}^{c_g}$ , with $c_g$ denoting the number of components of the graph $\mathcal{G}$ ), yielding again the reduced state space¹⁹ $\text{im } B^T \times \Lambda_0$ .

For a mass-spring-damper system, although not considered as a Hamiltonian system in the standard symmetry reduction framework, the same reduction procedure can still be applied. A mass-spring-damper system in coordinates $(q_c, p)$ takes the form

$(57) \quad \begin{aligned} \dot{q}_c &= \frac{\partial H_c}{\partial p}(q_c, p), \\ \dot{p} &= -\frac{\partial H_c}{\partial q_c}(q_c, p) - B_d R B_d^T \frac{\partial H_c}{\partial p}(q_c, p) + E f_b, \\ e^b &= E^T \frac{\partial H_c}{\partial p}(q_c, p), \end{aligned}$

where $H_c(q_c, p) = H(B_s^T q, p)$ , with $q = B_s^T q_c$ the spring elongations. Here $B_s$ and $B_d$ denote, as before, the incidence matrices of the spring, respectively, damper graph. Under the same symmetry action as above this reduces to (36) on the reduced state space $\text{im } B_s^T \times \Lambda_0$ .

The precise relation between Theorem 4.3 and the results obtained in [7, 16] now becomes clear. Indeed, the double-integrator networks studied in [16, 7] correspond to linear mass-spring-damper systems with unit masses, unit spring constants, and unit damping coefficients, expressed in the position variables $q_c$ and the velocities $\dot{q}_c$ , which are equal to the momenta $p$ . Thus Theorem 4.3 can be seen as a direct extension of the velocity consensus result expressed in Corollary 10 of [7]. Note, on the other hand, that thanks to the systematic use of the port-Hamiltonian structure, the stability treatment given in section 4 is directly extendable to the nonlinear case. Furthermore we obtain the following corollary to Theorem 4.3 regarding “second-order consensus.”

¹⁹Note that in fact the subspace $\text{im } B^T \subset \Lambda^1$ is determined by the Casimirs $k^T q, Bk = 0$ in the sense that $\text{im } B^T = \{q \in \Lambda^0 \mid k^T q = 0 \text{ for all } k \in \ker B\}$ . Furthermore, $\text{im } B^T = \Lambda^1$ if and only if the graph does not contain cycles.

COROLLARY 5.1. Consider the mass-spring-damper system (57) in coordinates $(q_c, p)$ , where we assume the spring graph to be connected. Then for all initial conditions $q_c(t) \rightarrow \text{span } \mathbb{1}, p(t) \rightarrow \text{span } \mathbb{1}$ if and only the largest $GL_s$ -invariant subspace contained in $\ker B_s^T$ is equal to $\text{span } \mathbb{1}$ , and moreover $\ker B_s = 0$ .

Proof. Compared to Theorem 4.3 on “velocity consensus” only the condition $\ker B_s = 0$ (no cycles in the spring graph) is new. However, it follows from the proof of Theorem 4.3 that $q(t) = B_s^T q_c(t) \rightarrow q_\infty$ with $B_s K q_\infty = 0$ . Thus $q(t) = B_s^T q_c(t) \rightarrow 0$ if and only if $\ker B_s = 0$ . Finally, by connectedness of the spring graph, $q(t) = B_s^T q_c(t) \rightarrow 0$ if and only if $q(t) \rightarrow \text{span } \mathbb{1}$ . $\square$

5.2. Further reduction. It is well known that symmetry reduction of a Hamiltonian system entails two steps [19]. Roughly speaking, the first step, as discussed above, deals with factoring out the state space by the action of the symmetry group, leading to a Hamiltonian system defined with respect to a Poisson structure possessing Casimirs. The second step deals with restricting the obtained Hamiltonian dynamics to the level sets of these Casimirs, thereby obtaining a reduced symplectic Hamiltonian system.

In the case of a mass-spring system (with connected graph) the second step is performed as follows. Note that the dual space $(\text{im } B^T)^*$ can be identified with

$(58) \quad (\text{im } B^T)^* = \Lambda_1 / (\text{im } B^T)^\perp = \Lambda_1 / \ker B \simeq \text{im } B.$

Thus the reduced state space can be identified with

$(59) \quad \text{im } B^T \times \text{im } B \subset \Lambda^1 \times \Lambda_0,$

which is again a symplectic space. We leave the extension to the nonconnected case and the presence of dampers to the readers.

6. The Kirchhoff–Dirac structure on graphs and its port-Hamiltonian dynamics. In this section we consider a third canonical graph Dirac structure, which results from constraining the flows at the internal vertices to zero (and thus there is no energy-storage or dissipation associated with the vertices for the corresponding port-Hamiltonian system).

6.1. The Kirchhoff–Dirac structure. As already alluded to at the end of section 2.6 the Kirchhoff–Dirac structure is defined as

$(60) \quad \begin{aligned} \mathcal{D}_K(\mathcal{G}) &:= \{(f_1, e^1, f_b, e^b) \in \Lambda_1 \times \Lambda^1 \times \Lambda_b \times \Lambda^b \mid \\ B_i f_1 &= 0, B_b f_1 = f_b, \exists e^{0i} \in \Lambda^{0i} \text{ such that } e^1 = -B_i^* e^{0i} - B_b^* e^b\}. \end{aligned}$

Note that in contrast to the flow/effort-continuous graph Dirac structures the Kirchhoff–Dirac structure only involves the flow and effort variables of the edge and boundary vertex spaces (not of the internal vertex spaces).

PROPOSITION 6.1. $\mathcal{D}_K(\mathcal{G})$ is a separable Dirac structure.

Proof. The Kirchhoff–Dirac structure is equal to the composition of the flow-continuous²⁰ graph Dirac structure $\mathcal{D}_f(\mathcal{G})$ with the trivial separable Dirac structure defined as

$\{(f_{0i}, e^{0i}) \in \Lambda_{0i} \times \Lambda^{0i} \mid f_{0i} = 0\}.$

²⁰Or the composition of the effort-continuous graph Dirac structure with $\{(f_0, e^0) \in \Lambda_0 \times \Lambda^0 \mid f_0 = 0\}$ .

The result then follows from Proposition 2.6. $\square$

Port-Hamiltonian systems with respect to the Kirchhoff–Dirac structure are defined completely similar to the case of the flow/effort-continuous graph Dirac structure, with the difference being that energy-storing or -dissipative relations are now only defined for the flow and effort variables corresponding to the edges.

6.2. Electrical circuits. The prime example of a port-Hamiltonian system²¹ with respect to a Kirchhoff–Dirac structure is an electrical RLC-circuit, with circuit graph $\mathcal{G}$ . In this case the elements of $\Lambda_1$ and $\Lambda^1$ denote the vectors of currents through, respectively, the voltages across, the edges, and the Kirchhoff–Dirac structure amounts to Kirchhoff’s current and voltage laws (whence its name). Furthermore, the effort variables $e^0$ are the potentials at the vertices, while the boundary flows and efforts $f_b, e^b$ are the boundary currents, respectively, boundary potentials at the boundary vertices (the terminals of the electrical circuit).

On top of Kirchhoff’s laws, the dynamics is defined by the energy-storage relations corresponding to either capacitors or inductors, and dissipative relations corresponding to resistors. The energy-storing relations for a capacitor at edge $e$ are given by

$(61) \quad \dot{Q}_e = -I_e, \quad V_e = \frac{dH_{Ce}}{dQ_e}(Q_e)$

with $Q_e$ the charge, and $H_{Ce}(Q_e)$ denoting the electric energy stored in the capacitor. Alternatively, in the case of an inductor one specifies the magnetic energy $H_{Le}(\Phi_e)$ , where $\Phi_e$ is the magnetic flux linkage, together with the dynamic relations

$(62) \quad \dot{\Phi}_e = V_e, \quad -I_e = \frac{dH_{Le}}{d\Phi_e}(\Phi_e).$

Finally, a resistor at edge $e$ corresponds to a static relation between the current $I_e$ through and the voltage $V_e$ across this edge, such that $V_e I_e \leq 0$ . In particular, a linear (ohmic) resistor at edge $e$ is specified by a relation $V_e = -R_e I_e$ , with $R_e \geq 0$ .

Alternatively, we can decompose the circuit graph $\mathcal{G}$ as the interconnection of a graph corresponding to the capacitors, a graph corresponding to the inductors, and a graph corresponding to the resistors. For simplicity let us restrict ourselves to the case of an $LC$ -circuit without boundary vertices. Define $\hat{\mathcal{V}}$ as the set of all vertices that are adjacent to at least one capacitor as well as to at least one inductor. Then split the circuit graph into an open circuit graph $\mathcal{G}^C$ corresponding to the capacitors and an open circuit graph $\mathcal{G}^L$ corresponding to the inductors, both with a set of boundary vertices $\hat{\mathcal{V}}$ . Denote the incidence matrices of these two circuit graphs by

$B^C := \begin{bmatrix} B^C \\ B_b^C \end{bmatrix}, \quad B^L := \begin{bmatrix} B^L \\ B_b^L \end{bmatrix}.$

Assuming for simplicity that all capacitors and inductors are linear, we arrive at the following equations for the $C$ -circuit:

$\begin{aligned} B_b^C \dot{Q} &= I_b^C, & B_i^C \dot{Q} &= 0, \\ B_b^{CT} \psi_b^C &= C^{-1} Q - B_i^{CT} \psi_i^C \end{aligned}$

²¹The terminology “port-Hamiltonian” may be confusing in this context, because “ports” in electrical circuits are usually defined by pairs of terminals, that is, pairs of boundary vertices with external variables being the currents through and the voltages across an edge corresponding to each such port. See also the discussions in [45, 46, 39].

with $Q$ the vector of charges of the capacitors and $C$ the diagonal matrix with diagonal elements given by the capacitances of the capacitors. Similarly, for the $L$ -circuit we obtain the equations

$\begin{aligned}\dot{\Phi} &= B_b^{LT} \psi_b^L + B_i^{LT} \psi_i^L, \\ 0 &= B_i^L L^{-1} \Phi, \\ I_b^L &= -B_b^L L^{-1} \Phi\end{aligned}$

with $\Phi$ the vector of fluxes and $L$ the diagonal matrix of inductances of all the inductors.

The equations of the $LC$ -circuit are obtained by imposing the interconnection constraints $\psi_b^C = \psi_b^L =: \psi_i$ and $I_b^C + I_b^L = 0$ . By eliminating the boundary currents $I_b^C, I_b^L$ one thus arrives at the differential-algebraic port-Hamiltonian equations

$\begin{bmatrix} B_i^C & 0, \\ 0 & B_i^L \\ B_b^C & B_b^L \end{bmatrix} \begin{bmatrix} -\dot{Q} \\ L^{-1}\Phi \end{bmatrix} = 0,$ $\begin{bmatrix} C^{-1}Q \\ -\Phi \end{bmatrix} = \begin{bmatrix} B_i^{CT} & 0 & B_b^{CT} \\ 0 & B_i^{LT} & B_b^{LT} \end{bmatrix} \begin{bmatrix} \psi_i^C \\ \psi_i^L \\ \psi_i \end{bmatrix}.$

For a formulation of pure $R$ , $L$ , or $C$ circuits and their weighted Laplacian matrices, we refer the reader to [33].

6.3. Mass-spring systems with regard to a Lagrangian tree. An alternative port-Hamiltonian formulation of mass-spring(-damper) systems, in terms of the Kirchhoff-Dirac structure, can be given as follows. Recall the port-Hamiltonian formulation on $\Lambda^1 \times \Lambda_0$ with respect to the effort-continuous graph Dirac structure $\mathcal{D}_e(\mathcal{G})$ , in which case the masses correspond to the vertices, and the springs to the edges of the graph $\mathcal{G}$ , which we assume to be connected.²² This graph can be extended to an augmented graph $\mathcal{G}_{\text{aug}}$ by adding a ground vertex $g$ and adding edges from every vertex $v$ of $\mathcal{G}$ toward this ground vertex. (The augmented graph is called a Lagrangian tree.) Furthermore, by constraining the effort $e_g$ at the ground vertex to be zero we can equate the efforts $e_v$ at every vertex $v$ of $\mathcal{G}$ with the effort $e_{vg}$ at the edge from $v$ to $g$ of the augmented graph $\mathcal{G}_{\text{aug}}$ . In this way we can identify the effort-continuous graph Dirac structure $\mathcal{D}_e(\mathcal{G})$ with the Kirchhoff-Dirac structure $\mathcal{D}_K(\mathcal{G}_{\text{aug}})$ with the additional constraint $e_g = 0$ . (Note that this is again a separable Dirac structure since it equals the composition of the Kirchhoff-Dirac structure $\mathcal{D}_K(\mathcal{G}_{\text{aug}})$ with the trivial Dirac structure $\{(f_g, e_g) \mid e_g = 0\}$ .)

In this way, the masses become associated with the edges $e_{vg}$ from every vertex $v$ to the ground vertex $g$ . The interpretation of the ground vertex $g$ is that it represents the reference point (with velocity $e_g$ being zero). The flow $f_g$ at the ground vertex $g$ equals the total force exerted on a mass located at this reference point.

6.4. Properties of the boundary flows and efforts of the Kirchhoff-Dirac structure. The fact that the internal vertex flows in the definition of the Kirchhoff-Dirac structure are all zero (and consequently no storage or dissipation at the vertices takes place) has a number of specific consequences for the behavior of the boundary flows and efforts (see [46] for closely related considerations).

²²For nonconnected graphs $\mathcal{G}$ the same construction can be done for every connected component.

Assume (for simplicity of exposition) that $\mathcal{R} = \mathbb{R}$ . From the definition of the Kirchhoff–Dirac structure and $\mathbb{1}^T B = 0$ it follows that

$(63) \quad 0 = \mathbb{1}^T B f_1 = \mathbb{1}_b^T B_b f_1 = -\mathbb{1}_b^T f_b$

with $\mathbb{1}_b$ denoting the vector with all ones of dimension equal to the number of boundary vertices. Hence the boundary part of the Kirchhoff–Dirac structure of an open graph is constrained by the fact that the boundary flows add up to zero. Dually, we may always add to the vector of vertex efforts $e^0$ the vector $\mathbb{1}$ leaving invariant the edge efforts $e^1 = B^T e^0$ . Hence, to the vector of boundary efforts $e^b$ we may always add the vector $\mathbb{1}_b$ .

PROPOSITION 6.2. Consider an open graph $\mathcal{G}$ with Kirchhoff–Dirac structure $\mathcal{D}_K(\mathcal{G})$ . Then for each $(f_1, e^1, f_b, e^b) \in \mathcal{D}_K(\mathcal{G})$ it holds that

$\mathbb{1}_b^T f_b = 0,$

while for any constant $c \in \mathbb{R}$ ,

$(f_1, e^1, f_b, e^b + c\mathbb{1}_b) \in \mathcal{D}_K(\mathcal{G}).$

This proposition implies that we may restrict the dimension of the space of boundary flows and efforts $\Lambda_b \times \Lambda^b$ of a connected graph by two. Indeed, we may define

$\Lambda_{\text{bred}} := \{f_b \in \Lambda_b \mid f_b \in \ker \mathbb{1}_b^T\}$

and its dual space

$\Lambda_{\text{red}}^b := (\Lambda_{\text{bred}})^* = \Lambda^b / \text{im } \mathbb{1}_b.$

It is straightforward to show that the Kirchhoff–Dirac structure $\mathcal{D}_K(\mathcal{G})$ reduces to a linear subspace of the reduced space $\Lambda_1 \times \Lambda^1 \times \Lambda_{\text{bred}} \times \Lambda_{\text{red}}^b$ , which is also a Dirac structure. An interpretation of this reduction is that we may consider one of the boundary vertices, say the first one, to be the reference vertex, and that we may reduce the vector of boundary efforts $e^b = (e^{b1}, \dots, e^{bb})$ to a vector of voltages $(e^{b2} - e^{b1}, \dots, e^{bb} - e^{b1})$ . A graph-theoretic interpretation is that instead of the incidence matrix $B$ we consider the restricted incidence matrix [3].

For a graph $\mathcal{G}$ with more than one connected component, the above holds for each connected component.²³ It follows that there are as many independent constraints on the boundary flows $f_b$ as the number of the connected components of the open graph $\mathcal{G}$ . Dually, the space of allowed boundary efforts $e^b$ is invariant under translation by as many independent vectors $\mathbb{1}_b$ as the number of connected components.

A complementary view on Proposition 6.2 is the fact that we may close an open graph $\mathcal{G}$ so that it becomes a closed graph $\tilde{\mathcal{G}}$ as follows. Consider first the case that $\mathcal{G}$ is connected. Then we may add one virtual ground vertex $v_g$ , and virtual edges from this virtual vertex to every boundary vertex $v_b \in \mathcal{V}_e$ , in such a manner that the Kirchhoff–Dirac structure of this graph $\tilde{\mathcal{G}}$ extends the Kirchhoff–Dirac structure of the open graph $\mathcal{G}$ . In fact, to the virtual vertex $v_g$ we may associate an arbitrary

²³The rank of the incidence matrix is equal to the number of vertices minus the number of connected components [4]. In fact, each connected component of the graph satisfies the property $\ker B^T = \text{span } \mathbb{1}$ with $B$ the (restricted) incidence matrix of this component and the dimension of $\mathbb{1}$ equal to its number of vertices.

potential $e^0(v_g)$ (the ground potential), and we may rewrite the externally supplied power $\langle e^b | f_b \rangle$ as (since by (63), $\sum_{v_b} f_b(v_b) = 0$ )

$(64) \quad \langle e^b | f_b \rangle = \sum_{v_b} (e^b(v_b) - e^0(v_g)) f_b(v_b) = \sum_{v_b} e^{1b}(v_b) f_{1b}(v_b),$

where $e^{1b}(v_b) := e^b(v_b) - e^0(v_g)$ and $f_{1b}(v_b) := f_b(v_b)$ denote the effort across and the flow through the virtual edge toward the boundary vertex $v_b$ . It is clear that for every element $(f_1, e^1, f_b, e^b) \in \mathcal{D}_K(\mathcal{G})$ corresponding to the open graph $\mathcal{G}$ there exists $e^{1b}$ such that $(f_1, e^1, f_{1b}, e^{1b}) \in \tilde{\mathcal{G}}$ for the closed graph $\tilde{\mathcal{G}}$ , and conversely for every $(f_1, e^1, f_{1b}, e^{1b}) \in \tilde{\mathcal{G}}$ there exists $e^b$ such that $(f_1, e^1, f_b, e^b) \in \mathcal{D}_K(\mathcal{G})$ . This construction is extended to nonconnected graphs by adding a ground vertex to each component containing boundary vertices.

6.5. Physical analogies. From the above formulation of an RLC-circuit in section 6.2 we conclude that the structure of the dynamical equations of an inductor is different from that of a capacitor. In order to elucidate this basic difference we zoom in on the description of an inductor and a capacitor as two-terminal elements. To this end consider the elementary open graph consisting of one edge with two boundary vertices $\alpha, \beta$ , described by the incidence matrix $b = [1 \ -1]^T$ . It follows that an inductor with magnetic energy $H(\Phi)$ is described by the equations

$(65) \quad \begin{aligned} \dot{\Phi} &= b^T \begin{bmatrix} \psi_\alpha \\ \psi_\beta \end{bmatrix}, \\ \begin{bmatrix} I_\alpha \\ I_\beta \end{bmatrix} &= b \frac{dH}{d\Phi}(\Phi), \end{aligned}$

whereas a capacitor with electric energy $H(Q)$ is described as

$(66) \quad \begin{aligned} b\dot{Q} &= \begin{bmatrix} I_\alpha \\ I_\beta \end{bmatrix}, \\ \frac{dH}{dQ}(Q) &= b^T \begin{bmatrix} \psi_\alpha \\ \psi_\beta \end{bmatrix}. \end{aligned}$

This difference stems from the fact that the energy variable $Q$ of a capacitor, as well as the current $I$ , takes values in the linear space $\Lambda_1$ , while the state variable $\Phi$ of an inductor, as well as the voltage $V$ , takes values in the dual space $\Lambda^1$ . Recalling from section 3.2.1 the description of a spring system

$(67) \quad \begin{aligned} \dot{q} &= b^T \begin{bmatrix} v_\alpha \\ v_\beta \end{bmatrix}, \\ \begin{bmatrix} F_\alpha \\ F_\beta \end{bmatrix} &= b \frac{dH}{dq}(q), \end{aligned}$

with $q$ the elongation of the spring and $H(q)$ its potential energy, we conclude that there is a strict analogy between a spring and an inductor.²⁴ On the other hand, a moving mass is not a strict analogue of a capacitor. Instead, it can be considered as

²⁴Thus we favor the so-called force-current analogy instead of the force-voltage analogy.

the analogue of a grounded capacitor, while the strict analogue of a capacitor (66) is the so-called inertor [43]

$b\dot{p} = \begin{bmatrix} F_\alpha \\ F_\beta \end{bmatrix}, \quad \frac{dH}{dp}(p) = b^T \begin{bmatrix} v_\alpha \\ v_\beta \end{bmatrix},$

where $p$ is the momentum of the inertor and $H(p)$ its kinetic energy, while $F_\alpha, F_\beta$ and $v_1, v_2$ denote the forces, respectively, velocities, at the two terminals of the inertor.

7. Conclusions. We have laid down a general geometric framework for the description of physical network dynamics on (nonrandom) graphs. Starting points are the conservation laws corresponding to the incidence matrix of the graph. These define three canonical Dirac structures on the combined vertex, edge, and boundary spaces and their duals, where the last one (the Kirchhoff–Dirac structure) corresponds to the absence of energy storage or energy dissipation at the vertices. Relating the internal flows and efforts by either energy-storing or energy-dissipating relations yields various forms of port-Hamiltonian systems. We have illustrated the approach on a number of typical physical examples. Examples that have not been discussed include supply-chain models and compartmental systems. We have shown how examples with a different origin, such as consensus algorithms, can be formulated and analyzed within the same framework. Furthermore we have shown how classical techniques from Hamiltonian dynamical systems can be exploited for the analysis of the resulting port-Hamiltonian systems.

For clarity of exposition we have considered only the basic building blocks of port-Hamiltonian systems on graphs. Indeed, because the interconnection of port-Hamiltonian systems again defines a port-Hamiltonian system [10, 31, 8], the framework also covers heterogeneous and multiscale situations, where several of the constructs considered in the present paper are connected to each other. Moreover, as already indicated in section 2.6 and Remark 3.1, various interesting extensions to dynamical graphs and switching port-Hamiltonian systems on graphs can be made.

All of the models treated in this paper correspond to conservation/balance laws within a particular physical domain. Furthermore, the energy-balance of the system components can be seen to result from the underlying conservation laws and the assumption of integrable constitutive relations for energy-storage. On the other hand, port-based (bond-graph) modeling as originating in the work of Paynter [27] is aimed at providing a unifying modeling framework for multiphysics systems by directly starting from energy-flows between system components from different physical domains. This also results in port-Hamiltonian models, as has been amply demonstrated in, e.g., [21, 22, 36, 32, 14]. It is well known that bond-graph modeling involves an additional abstraction step (e.g., different electrical circuits may lead to the same bond-graph, and, conversely, different bond-graphs may correspond to the same electrical circuit). Furthermore, in the case of electrical circuits, port-based modeling starts with a port description (pairs of terminals) instead of the more basic starting point of terminals corresponding to conservation laws. Although in most situations the resulting port-Hamiltonian systems are the same, this leaves some questions to be answered; see also [46, 39]. Another interesting venue for further research [34] is the precise relation between port-Hamiltonian systems (on graphs) and gradient dynamical systems; see especially [5, 42] for the gradient formulation of RLC circuits.

The identification of the port-Hamiltonian structure, as already crucially used in section 4 for (stability) analysis, offers important tools for simulation and control. Port-Hamiltonian systems theory has been successful in exploiting the physical structure for control and design purposes (see, e.g., [32, 26]), using various forms of passivity-based control, control by interconnection, and tools originating in network synthesis theory. The applications of this control methodology to port-Hamiltonian systems on graphs is an important area for further research. The combination with graph theory, and the inclusion of constraints on the flow and storage variables, are very promising; see [41] for preliminary work in this direction.

In a companion paper we will extend the framework from directed graphs to general $k$ -complexes. This allows us to give a spatially discretized model of the two-dimensional Maxwell equations and of general diffusive systems; see [38, 39]. Pertinent questions include identifying the relation of these methods to structure-preserving spatial discretization methods of partial differential equation models.

2. REFERENCES

[1] M. ARCAK, Passivity as a design tool for group coordination, IEEE Trans. Automat. Control, 52 (2007), pp. 1380–1390.
[2] R. S. BALL, A Treatise on the Theory of Screws, Cambridge University Press, Cambridge, UK, 1998. Reprint of the 1900 edition.
[3] P. BAMBERG AND S. STERNBERG, A Course in Mathematics for Students of Physics 2, Cambridge University Press, Cambridge, UK, 1990.
[4] B. BOLLOBAS, Modern Graph Theory, Grad. Texts Math. 184, Springer, New York, 1998.
[5] R. BRAYTON AND J. MOSER, A theory of nonlinear networks. I, II, Quart. Appl. Math., 22 (1964), pp. 1–33, 81–104.
[6] M. BÜRGER, D. ZELAZO, AND F. ALLGÖWER, Network clustering: A dynamical systems and saddle-point perspective, in Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference (CDC-ECC), Orlando, FL, 2011, pp. 7825–7830.
[7] M. K. CAMLIBEL AND S. ZHANG, Partial consensus for heterogeneous multi-agent systems with double integrator dynamics, in Proceedings of the 4th IFAC Conference on Analysis and Design of Hybrid Systems, Eindhoven, The Netherlands, 2012, pp. 121–126.
[8] J. CERVERA, A. J. VAN DER SCHAFT, AND A. BAÑOS, Interconnection of port-Hamiltonian systems and composition of Dirac structures, Automatica J. IFAC, 43 (2007), pp. 212–225.
[9] T. J. COURANT, Dirac manifolds, Trans. Amer. Math. Soc., 319 (1990), pp. 631–661.
[10] M. DALSMO AND A. VAN DER SCHAFT, On representations and integrability of mathematical structures in energy-conserving physical systems, SIAM J. Control Optim., 37 (1998), pp. 54–91.
[11] T. H. DAVIES, Mechanical networks—I: Passivity and redundancy. II: Formulae for the degrees of mobility and redundancy. III: Wrenches on circuit screws, Mech. Mach. Theory, 18 (1983), pp. 95–101, 103–106, 107–112.
[12] C. DE PERSIS AND C. S. KALLESOE, Pressure regulation in nonlinear hydraulic networks by positive and quantized controls, IEEE Trans. Control Systems Technology, 19 (2011), pp. 1371–1383.
[13] I. YA. DORFMAN, Dirac Structures and Integrability of Nonlinear Evolution Equations, John Wiley, New York, 1993.
[14] V. DUINDAM, A. MACCHELLI, S. STRAMIGIOLI, AND H. BRUYNINCKX, EDs., Modeling and Control of Complex Physical Systems—The Port-Hamiltonian Approach, Springer, New York, 2009.
[15] K. GERRITSEN, A. J. VAN DER SCHAFT, AND W. P. M. H. HEEMELS, On switched Hamiltonian systems, In Proceedings of the 15th International Symposium on Mathematical Theory of Networks and Systems (MTNS 2002), D. S. Gilliam and J. Rosenthal, eds., South Bend, IN, 2002.
[16] D. GOLDIN, S. A. ATTIA, AND J. RAISCH, Consensus for double integrator dynamics in heterogeneous networks, in Proceedings of the 49th IEEE Conference on Decision and Control (CDC), Atlanta, GA, 2010, pp. 4504–4510.
[17] A. KARGER AND J. NOVAK, Space Kinematics and Lie Groups, Gordon and Breach Science Publishers, New York, 1985.
[18] P. LIBERMANN AND C.-M. MARLE, Symplectic Geometry and Analytical Mechanics, D. Reidel Publishing Company, Dordrecht, The Netherlands, 1987.

[19] J. E. MARSDEN AND T. S. RATIU, Introduction to Mechanics and Symmetry: A Basic Exposition to Classical Mechanics, 2nd ed., London Math. Soc. Lecture Notes Ser., Springer, New York, 1999.
[20] B. M. MASCHKE AND A. J. VAN DER SCHAFT, Interconnected mechanical systems. Parts 1 and 2, in Modelling and Control of Mechanical Systems, Imperial College Press, London, 1997, pp. 1–30.
[21] B. M. MASCHKE, A. J. VAN DER SCHAFT, AND P. C. BREEDVELD, An intrinsic Hamiltonian formulation of network dynamics: Non-standard Poisson structures and gyrators, J. Franklin Institute, 329 (1992), pp. 923–966.
[22] B. M. MASCHKE, A. J. VAN DER SCHAFT, AND P. C. BREEDVELD, An intrinsic Hamiltonian formulation of the dynamics of LC-circuits, IEEE Trans. Circuits Systems I Fund. Theory Appl., 42 (1995), pp. 73–82.
[23] B. MASCHKE AND A. J. VAN DER SCHAFT, Conservation Laws and Port-Hamiltonian Systems on Complexes, in preparation.
[24] R. M. MURRAY, X. LI, AND S. S. SASTRY, A Mathematical Introduction to Robotic Manipulation, CRC Press, Boca Raton, FL, 1994.
[25] R. OLFATI-SABER, J. A. FAX, AND R. M. MURRAY, Consensus and cooperation in networked multi-agent systems, Proc. IEEE, 95 (2007), pp. 215–233.
[26] R. ORTEGA, A. J. VAN DER SCHAFT, I. MAREELS, AND B. MASCHKE, Putting energy back in control, IEEE Control Systems Mag., 21 (2001), pp. 18–32.
[27] H. M. PAYNTER, Analysis and Design of Engineering Systems, MIT Press, Cambridge, MA, 1961.
[28] A. RAHMANI, M. JI, M. MESBAHI, AND M. EGERSTEDT, Controllability of multi-agent systems from a graph-theoretic perspective, SIAM J. Control Optim., 48 (2009), pp. 162–186.
[29] J. A. ROBERSON AND C. T. CROWE, Engineering Fluid Mechanics, Houghton Mifflin Company, Boston, MA, 1993.
[30] J. M. SELIG, Geometric Methods in Robotics, Monogr. Comput. Sci., Springer-Verlag, New York, 1996.
[31] A. J. VAN DER SCHAFT, Interconnection and geometry, in The Mathematics of Systems and Control: From Intelligent Control to Behavioral Systems, J. W. Polderman and H. L. Trentelman, eds., University of Groningen, Groningen, The Netherlands, 1999, pp. 203–218.
[32] A. J. VAN DER SCHAFT, L₂-Gain and Passivity Techniques in Nonlinear Control, 2nd revised and enlarged ed., Comm. Control Engrg. Ser., Springer-Verlag, London, 2000.
[33] A. J. VAN DER SCHAFT, Characterization and partial synthesis of the behavior of resistive circuits at their terminals, Systems Control Lett., 59 (2010), pp. 423–428.
[34] A. J. VAN DER SCHAFT, On the relation between port-Hamiltonian and gradient systems, in Proceedings of the 18th IFAC World Congress, Milan, Italy, 2011, pp. 1321–1326.
[35] A. J. VAN DER SCHAFT AND M. K. CAMLIBEL, A state transfer principle for switching port-Hamiltonian systems, in Proceedings of the 48th IEEE Conference on Decision and Control, Shanghai, China, 2009, pp. 45–50.
[36] A. J. VAN DER SCHAFT AND B. M. MASCHKE, The Hamiltonian formulation of energy conserving physical systems with external ports, Arch. Elek. Übertr., 49 (1995), pp. 362–371.
[37] A. J. VAN DER SCHAFT AND B. M. MASCHKE, Hamiltonian formulation of distributed parameter systems with boundary energy flow, J. Geom. Phys., 42 (2002), pp. 166–174.
[38] A. J. VAN DER SCHAFT AND B. M. MASCHKE, Conservation laws and open systems on higher-dimensional networks, in Proceedings of the 47th IEEE Conference on Decision and Control, Cancun, Mexico, 2008, pp. 799–804.
[39] A. J. VAN DER SCHAFT AND B. M. MASCHKE, Conservation laws and lumped system dynamics, in Model-Based Control: Bridging Rigorous Theory and Advanced Technology, P. M. J. Van den Hof, C. Scherer, and P. S. C. Heuberger, eds., Springer, Berlin, Heidelberg, 2009, pp. 31–48.
[40] A. J. VAN DER SCHAFT AND B. MASCHKE, Port-Hamiltonian dynamics on graphs: Consensus and coordination control algorithms, in Proceedings of the 2nd IFAC Workshop on Distributed Estimation and Control in Networked Systems, Centre de Congrès de L'Impérial Palace, Annecy, France, 2010, pp. 175–178.
[41] A. J. VAN DER SCHAFT AND J. WEI, A Hamiltonian perspective on the control of dynamical distribution networks, in Proceedings of the 4th IFAC Workshop on Lagrangian and Hamiltonian Methods for Non Linear Control (LHMNLC 2012), Bertinoro, Italy, 2012, pp. 24–29.
[42] S. SMALE, On the mathematical foundations of electrical circuit theory, J. Differential Geom., 7 (1972), pp. 193–210.

[43] M. C. SMITH, Synthesis of mechanical networks: The inerter, IEEE Trans. Automat. Control, 47 (2002), pp. 1648–1662.
[44] J. VANKERSCHAUER, H. YOSHIMURA, H. LEOK, AND J. E. MARSDEN, Stokes-Dirac structures through reduction of infinite-dimensional Dirac structures, in Proceedings of the 49th IEEE Conference on Decision and Control, Atlanta, GA, 2010, pp. 6265–6270.
[45] J. C. WILLEMS, The behavioral approach to open and interconnected systems: Modeling by tearing, zooming, and linking, IEEE Control Syst. Mag., 27 (2007), pp. 46–99.
[46] J. C. WILLEMS, Terminals and ports, IEEE Circuits Systems Mag., 10 (2010), pp. 8–16.
[47] D. ZELAZO AND M. MESBAHI, Edge agreement: Graph-theoretic performance bounds and passivity analysis, IEEE Trans. Automat. Control, 56 (2011), pp. 544–555.