Geometric Modeling for Control of Thermodynamic Systems

Abstract

This paper discusses the way that energy and entropy can be regarded as storage functions with respect to supply rates corresponding to the power and thermal ports of the thermodynamic system. Then, this research demonstrates how the factorization of the irreversible entropy production leads to quasi-Hamiltonian formulations, and how this can be used for stability analysis. The Liouville geometry approach to contact geometry is summarized, and how this leads to the definition of port-thermodynamic systems is discussed. This notion is utilized for control by interconnection of thermodynamic systems.

1. Introduction

Since the early 1970s [1] contact geometry has been recognized as underlying macroscopic thermodynamics, starting from Gibbs’ fundamental thermodynamic relation. This has spurred a series of papers on the geometry of thermodynamics; including [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]; see [28] for an introduction and survey. Nevertheless, this literature points to major differences with, for example, the geometric theory of classical mechanics (using symplectic geometry), and hints at aspects which have not yet been addressed. First, the thermodynamic phase space (which is formulated as a contact manifold) comprises the extensive and intensive variables, and thus, its dimension is more than twice the minimal number of variables to describe the thermodynamic system at any moment of time. Second, most of the theory is about thermostatics, and the proper geometric formulation of the dynamics is much less clear. Third, the contact geometric approach to thermodynamics is usually based on the energy representation of thermodynamic systems and its corresponding Gibbs one-form. On the other hand, there is also the entropy formulation which corresponds to another (although conformally equivalent) one-form. This led [29] to the use of homogeneous coordinates for the intensive variables, and thus to extend the thermodynamic phase space by one more degree of freedom. This was followed up in [25,26] by emphasizing the formulation of the thermodynamic phase space as the projection of the cotangent bundle over the space of extensive variables. Thus contact geometry is approached from the vantage point of the geometry of cotangent bundles with their Liouville one-form. Fourth, until now, not much work has been performed regarding the geometry of irreversible thermodynamics, based on the factorization of the irreversible entropy production. Fifth, how to use these geometric frameworks for the control of thermodynamic systems has not yet been addressed.

The present paper continues the investigation of all of these aspects. In Section 2, a systems and control perspective on macroscopic thermodynamics is emphasized by primarily regarding thermodynamic systems as systems interacting with their surroundings via heat, mechanical work, exchange of chemical species, etc. A classical example is, of course, the heat engine. A summary of how dissipativity theory provides a natural framework for interpreting and formulating the first and second laws of thermodynamics, Clausius’ inequality, and eventually entropy is provided. Indeed, energy and entropy reveal themselves to be the storage functions corresponding to two supply rates involving the thermal and mechanical ports of the thermodynamic system. Finally, this leads to Gibbs’ fundamental relation and to the definition of the thermodynamic phase space. Section 3 focuses on geometric descriptions of irreversible thermodynamic systems. The way that the classical factorization of the irreversible entropy production suggests quasi-Hamiltonian formulations (somewhat resembling GENERIC [30]) based on energy conservation and the increase of entropy of the autonomous part of the dynamics are discussed. This paper also indicates how such formulations may be used for stability analysis. Section 4 starts with the geometry of the thermodynamic phase space from the point of view of the Liouville geometry of the cotangent bundle over the space of extensive variables. Identifying the constitutive relations (‘thermostatics’) of the thermodynamic system as a Liouville submanifold, and the dynamics as homogeneous Hamiltonian dynamics lead to the definition of a port-thermodynamic system. Such systems interact with their environment via power ports and/or entropy flow ports. In Section 5, this is used for ‘control by interconnection’ of port-thermodynamic systems, where the dynamics of the system are sought to be controlled by interconnection with a suitable controller port-thermodynamic system. Finally, Section 6 contains conclusions and a discussion of venues for further research.

The sections are illustrated by three running examples: the gas-piston-damper system, chemical reaction networks, and the heat exchanger. Overall, the paper heavily builds upon previous papers [25,26,31,32,33,34], in which further details and background can be found.

2. The First and Second Law from the Point of View of Dissipativity Theory

The first law of thermodynamics expresses two fundamental properties: (1) the different types of interaction of a thermodynamic system with its surroundings (e.g., heat flow, mechanical work, flow of chemical species, etc.) all result in an exchange of a common quantity called energy, (2) there exists a function of the macroscopic thermodynamic variables that represents the energy stored in the system, and the increase of this function during any time interval is equal to the sum of the energies supplied to the system during this time interval by its surroundings (conservation of energy). Thus, energy manifests itself in different physical forms, which are equivalent and to a certain extent exchangeable. ‘To a certain extent’ because, as expressed by the second law of thermodynamics, there are limitations to the conversion of heat to other forms of energy.

The first law can be mathematically formulated through the use of dissipativity theory as formulated in [35]; see also [31,36,37]. Consider a simple thermodynamic system such as a gas, described by three variables: volume V, pressure P, and temperature T. Then, mechanical power (rate of mechanical work) provided by the surroundings to the thermodynamic system is given by $-P{u}_V$, where $u_V:=\dot{V}$ is the rate of volume change. (In the physics convention for the pressure P, $P{u}_V$ is the rate of mechanical work exerted by the system on the surroundings). The second type of interaction with the surroundings comes from heat delivered to the system (for instance, from a heat source). Let us denote, using q, the heat flow (heat per second) from the heat source into the system. Then the first law is expressed by the existence of a function $E\left(x\right)$ of the thermodynamic state x (e.g., $(V,P,T)$ satisfying the equation of state), expressing the energy of the system and satisfying, at all times, t

$$\frac{d}{dt}E\left(x\left(t\right)\right)=q\left(t\right)-P\left(t\right)u_V\left(t\right)$$

(1)

That is, the increase of the total energy E of the thermodynamic system is equal to the incoming heat flow (through the thermal port) minus the mechanical work performed by the system on its surroundings (through the mechanical port). Equivalently, in the terminology of dissipativity theory, the first law amounts to the system being cyclo-lossless for the supply rate $q-P{u}_V$, with storage functionE. This is directly extended to more involved thermodynamic systems. For example, suppose that apart from mechanical and thermal interaction with the surroundings, there is additional mass inflow of chemical species. Then, the supply rate $q-P{u}_V$ is extended to $q-P{u}_V+{\sum }_k{\mu }_k{\nu }_k$. Here, ${\nu }_k=\frac{d{N}_k}{dt}$, with $N_k$ the mole number of the k-th chemical species, and ${\mu }_k$ its chemical potential.

The first law emphasizes the role of thermodynamic systems as devices for energy conversion; energy from one physical domain is converted into energy in another domain. ‘Optimal’ conversion of heat into mechanical work, motivated by the design of steam engines in the beginning of the 19th century, was one of the starting points of thermodynamic theory. Electro-chemical devices such as batteries, and electro-mechanical systems including electrical motors and generators, are among the many other classical examples of energy-converting devices [38]. On the other hand, almost from the very start of thermodynamic theory, it was realized that there are intrinsic limitations to energy conversion. In particular, heat cannot just be converted into mechanical work. This is the origin of the second law of thermodynamics. The second law also admits a dissipativity interpretation; however, more involved than that of the first law. Let us start with the formulation of the second law, as given by Lord Kelvin (see [39]):

A transformation of a thermodynamic system whose only final result is to transform into work heat extracted from a source which is at the same temperature throughout is impossible.

Since the work done during a time interval $[t_1,t_2]$ is equal to ${\int }_{t_1}^{t_2}-P\left(t\right)dV\left(t\right)={\int }_{t_1}^{t_2}-P\left(t\right)u_V\left(t\right)dt$, Kelvin’s formulation immediately implies that for each constant temperature T, any thermodynamic system is cyclo-passive with respect to the supply rate $-P{u}_V$. However, the second law is stronger than that. Namely, Kelvin’s formulation also forbids the conversion into work of heat from a source at constant temperature for all transformations in which the system interacts as well with a second heat source at another temperature, as long as the net heat taken from this second heat source is zero. As demonstrated by Carnot, the interaction with heat sources at different temperatures is crucial for the conversion of heat into mechanical energy. This led to the famous Carnot cycle which can be described as follows: consider a simple thermodynamic system, in particular, a fluid or gas in a confined space of a certain volume. Control of the system functions in two ways: (1) via isothermal transformations, where heat is supplied to, or taken from, the system at a constant temperature (classically described as the interconnection of the thermodynamic system with an infinite heat reservoir at the temperature of the isothermal process), (2) via adiabatic transformations, where the only interaction with the surroundings is via work supplied to, or taken from, the system (classically described by the movement of a piston that changes the volume of the system, with a pressure equal to the pressure of the gas). A cycle consists of two isothermal transformations and two adiabatic transformations: first, an isothermal transformation at temperature $T_h$ (‘hot’) takes the system from an initial state to another state, secondly, an adiabatic transformation lowers the temperature of the system to $T_c$ (‘cold’), thirdly, an isothermal transformation at temperature $T_c$ takes the system to a state from which, fourthly, an adiabatic transformation takes the system back to the original initial state; see Figure 1.

The cycle is called a Carnot cycle if it is reversible; i.e., can be traversed in the opposite direction as well.

Remark 1. 

In the exposition of the Carnot cycle, often terminologies such as ‘infinitesimally slow’, ‘quasi-reversible’, ‘quasi-static’, etc., are used. This is largely with regard to the interaction of a system with its surroundings as being implemented by actual physical devices. For example, an isothermal transformation is viewed as the result of the ‘real’ physical action of a force exerted by a piston on the gas (implying that the pressure delivered by the piston could be different from the pressure of the gas). Furthermore, the system is considered to be in ‘real’ physical contact with a heat reservoir at a certain temperature (which could differ from that of the gas). In contrast in, e.g., electrical network theory and control theory the concept of an ‘ideal’ control action is employed, where, for instance, the pressure and the temperature are directly controlled.

The heat delivered to the system during the first isothermal at temperature $T_h$ is denoted by $Q_h$, and during the second isothermal at $T_c$ by $Q_c$ (generally $Q_c$ is negative). Then, by the first law, since the final state is equal to the initial state, $Q_h+Q_c=W$, where $W=\int PdV$ is the mechanical work that is done by the thermodynamic system on its surroundings.

By an intricate reasoning from [39], see also [31], Kelvin’s formulation of the second law yields for any cycle the fundamental inequality

$$\frac{Q_h}{T_h}+\frac{Q_c}{T_c}\le 0,$$

(2)

with equality in the case of a Carnot cycle. Furthermore, the reasoning can be extended to complex cycles, consisting of n isothermals at temperatures $T_i$ and absorbed heat quantities $Q_i$, $i=1,2,\cdots ,n$, interlaced by n adiabatics, leading to

$$\sum _{i=1}^n\frac{Q_i}{T_i}\le 0,$$

(3)

with equality in the case of reversibility. Finally, a slight extension (approximating continuous heat flow time-functions $q(·)$ by step functions with step values $Q_1,\cdots ,Q_n$) yields the celebrated Clausius inequality

$$\oint \frac{q\left(t\right)}{T\left(t\right)}dt\le 0$$

(4)

for all cyclic processes $\left(q(·),T(·)\right)$ (where q is the heat flow into the thermodynamic system, and T is the temperature of the system), with equality

$$\oint \frac{q\left(t\right)}{T\left(t\right)}dt=0$$

(5)

holding for all reversible cyclic processes (see [31] for details and refinements).

From the point of view of dissipativity theory [31,35] the Clausius inequality (4) is the same as cyclo-dissipativity of the thermodynamic system with respect to the supply rate $-\frac{q}{T}$. Thus, assuming reachability from and controllability from some ground state $x^{*}$ this means, see [40], that there exists a storage function F for the supply rate $-\frac{q}{T}$, that is $\frac{d}{dt}F\le -\frac{q}{T}$. Hence $S:=-F$ satisfies

$$\frac{d}{dt}S\ge \frac{q}{T}$$

(6)

The function S was called ‘entropy’ by Clausius, from the Greek word $\tau \rho o\pi \eta$ for ‘transformation’.

From the point of view of dissipativity theory, the storage function F need not be unique. In order to guarantee the uniqueness of F (modulo a constant), and therefore of the entropy S, we additionally assume [31,40] that, given some ground state, for every thermodynamic state there exists a reversible cyclic transformation through this state and the ground state satisfying

$$\oint \frac{q\left(t\right)}{T\left(t\right)}=0$$

(7)

This uniqueness of S is, explicitly or implicitly, always assumed in expositions of macroscopic thermodynamics, and also in this paper.

Remark 2. 

The dissipativity theory formulation of the second law already appears in [35], but under the additional assumption that F is nonnegative. In fact, in [35,37] there exists a nonnegative storage function for the supply rate $-\frac{q}{T}$ (and thus the system is dissipative instead of merely cyclo-dissipative) if and only if for all initial conditions x

$$F_a\left(x\right)=sup{\int }_0^{\tau }-\frac{q\left(t\right)}{T\left(t\right)}dt<\infty ,$$

(8)

where the supremum is taken over all $\tau \ge 0$ and all heat flow functions $q(·)$ on the time interval $[0,\tau ]$, and corresponding temperature profiles $T(·)$ resulting from $x\left(0\right)=x$. Furthermore, if (8) holds, then $F_a\ge 0$ is minimal among all nonnegative storage functions. It follows that $S_a=-F_a\le 0$ given by

$$S_a\left(x\right)=inf{\int }_0^{\tau }\frac{q\left(t\right)}{T\left(t\right)}dt>-\infty$$

(9)

is maximal among all nonpositive storage functions. Since an arbitrary constant may be added to S while still satisfying (6), the assumption that S is nonpositive is equivalent to S being bounded from above. However, in many thermodynamic systems the entropy is not bounded from above. Thus thermodynamic systems are generally only cyclo-dissipative with respect to the supply rate $-\frac{q}{T}$, and not dissipative.

The Thermodynamic Phase Space and Gibbs’ Relation

The next step is now to add the energy and entropy as extra extensive variables to the description of the thermodynamic system. In order to illustrate this, consider a simple thermodynamic system, with extensive variable V (volume) and intensive variables $P,T$ (pressure and temperature). The equation of state is an equation $f(V,P,T)=0$ for some scalar function f. (For example, for an ideal gas $PV=RT$ with R the universal gas constant.) Any $(V,P,T)$ satisfying $f(V,P,T)=0$ is called a state of the thermodynamic system. Hence, under regularity conditions the set of states of the thermodynamic system is a 2-dimensional submanifoldM of ${\mathbb{R}}^3$. Then, consider the functions $E:M\to \mathbb{R}$ (energy) and $S:M\to \mathbb{R}$ (entropy) as obtained from dissipativity theory. Then, we may equally well represent the set of states $M\subset {\mathbb{R}}^3$ by the 2-dimensional submanifold $L\subset {\mathbb{R}}^5$ comprising the total set of extensive and intensive variables ${\mathbb{R}}^5$:

$$L:=\left\{\right(E,S,V,T,P)\mid f(V,P,T)=0,E=E(V,P,T),S=E(V,P,T\left)\right\}$$

(10)

(With some abuse of notation, the extra variables $E,S$, are denoted by the same letters as used for the functions defined before.) The space ${\mathbb{R}}^5$ of all extensive and intensive variables is called the thermodynamic phase.

Furthermore, by the first law $\frac{d}{dt}E=-P\frac{d}{dt}V+q$, while for any state there exists a path through this state and the ground state such that $\frac{d}{dt}E=-P\frac{d}{dt}V+T\frac{d}{dt}S$. Taken together, this implies that the Gibbs one-form on the thermodynamic phase space ${\mathbb{R}}^5$ defined as

$$dE-TdS+PdV,Gibbsone-form,$$

(11)

is zero restricted to L. This is called Gibbs’ fundamental thermodynamic relation. The thermodynamic phase space, together with the Gibbs one-form, defines a contact manifold. Furthermore, a submanifold of the thermodynamic phase space ${\mathbb{R}}^5$ restricted to which the Gibbs one-form (11) is zero, and moreover has maximal dimension (in this case 2), is called a Legendre submanifold. Gibbs’ fundamental relation implies that any Legendre submanifold L is actually given as

$$L:=\{(E,S,V,T,P)\mid E=E(S,V),T=\frac{\partial E}{\partial S}(S,V),-P=\frac{\partial E}{\partial V}(S,V)\}$$

(12)

for some energy functions $E(S,V)$. Thus, L is completely described by expressing the energy E as a function $E(S,V)$ of the other two extensive variables $S,V$, hence the name energy representation. Instead of relying on such an energy function (or its partial Legendre transforms), there is still an alternative way of describing L. This is to start, not with $E(S,V)$, but instead with the expression of the entropy as a function $S(E,V)$. For a simple thermodynamic system this leads to the entropy representation of the submanifold $L\subset {\mathbb{R}}^5$ given as

$$L:=\{(E,S,V,T,P)\mid S=S(E,V),\frac{1}{T}=\frac{\partial S}{\partial E}(E,V),\frac{P}{T}=\frac{\partial S}{\partial V}(E,V)\}$$

(13)

3. Irreversible Thermodynamics

Clausius interpreted the term $\frac{q}{T}$ in the inequality (6) as the part of the infinitesimal transformation $\frac{d}{dt}S$ that is compensated by the opposite rate of change $-\frac{q}{T}$ of the entropy of the surroundings; that is, of the reservoir supplying the heat to the thermodynamic system. The remaining part

$$\sigma :=\frac{d}{dt}S-\frac{q}{T}\ge 0$$

(14)

was called the ‘uncompensated transformation’ by Clausius, and later the irreversible entropy production [38]. Irreversible thermodynamics is concerned with thermodynamics where $\sigma$ is different from zero, implying an autonomous (independent from external heat flow) increase of the entropy S. Sometimes it is also referred to as non-equilibrium thermodynamics, because the entropy increase is resulting from (internal) non-equilibrium conditions.

The standard postulate of irreversible thermodynamics (see e.g., [38]) is that $\sigma$ can be factorized as

$$\sigma =\sum _{k=1}^s{F}_k{J}_k\ge 0,$$

(15)

where $F_k$ are called the thermodynamic forces and $J_k$ are the thermodynamic flows (or fluxes), in such a way that

$$\sigma =0\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}{F}_k=0,k=1,\cdots ,s$$

(16)

In linear irreversible thermodynamics [38] it is furthermore assumed that the vectors F and J with components $F_k$ and $J_k,k=1,\cdots ,s,$ are related by a symmetric linear map

$$J=LF,L=L^{\top }$$

(17)

These are the celebrated Onsager reciprocity relations [38], corresponding to the symmetric factorization $\sigma =F^{\top }LF$.

Example 1 

(The heat exchanger). The perhaps simplest example of irreversible dynamics and irreversible entropy production is offered by the heat exchanger. Consider two heat compartments, having temperatures $T_h$ and $T_c$ (‘hot’ and ‘cold’), connected by a heat-conducting wall. In the absence of the conducting wall (and thus, without irreversible entropy production), these are two separate systems with entropies $S_h$ and $S_c$, each satisfying

$$\frac{d}{dt}{S}_h=\frac{q_h}{T_h},\frac{d}{dt}{S}_c=\frac{q_c}{T_c}$$

(18)

Due to the conducting wall, there is a heat flow q from the hot to the cold compartment, which is given by Fourier’s law for heat conduction as$q=\lambda (T_h-T_c)$for some positive constant$\lambda$. Furthermore, in view of the first law$q=-q_h=q_c$. Hence, the total entropy$S:=S_h+S_c$satisfies

$$\frac{d}{dt}S=-\frac{q}{T_h}+\frac{q}{T_c}=\left(\frac{1}{T_c}-\frac{1}{T_h}\right)q$$

(19)

This yields the following expression for the irreversible entropy production$\sigma =\frac{d}{dt}S$due to heat conduction (non-equilibrium conditions)

$$\sigma =\left(\frac{1}{T_c}-\frac{1}{T_h}\right)\lambda (T_h-T_c)=\lambda \frac{{\left(T_h-T_c\right)}^2}{T_h{T}_c}\ge 0$$

(20)

In this example the thermodynamic force is$F=\frac{1}{T_c}-\frac{1}{T_h}$, while the thermodynamic flow is$J=q=\lambda (T_h-T_c)$. Indeed$\sigma =0$if, and only if,$F=0$. Despite its simplicity, this is an example of nonlinear irreversible thermodynamics, since the thermodynamic flow$q=\lambda (T_h-T_c)$cannot be expressed as a linear function of the thermodynamic force$F=\frac{1}{T_c}-\frac{1}{T_h}$.

Example 2 

(The gas-piston-damper system). Another simple example is the gas-piston-damper system. Consider a cylinder containing a gas whose volume can be controlled by a piston actuated by an external force u, and is subject to linear damping. The total energy E of the system can be expressed as a function of the other extensive variables as

$$E(S,V,\pi )=U(S,V)+\frac{{\pi }^2}{2m},$$

(21)

with S representing entropy, V volume, $\pi$ momentum of the piston with mass m, and $U(S,V)$ representing the internal energy of the gas. Assuming that the heat as produced by the damping of the piston is fully absorbed by the gas in the cylinder, the dynamics are given as

$$\begin{array}{ccc}\hfill \dot{V}& =& A\frac{\pi }{m}\hfill \\ \hfill \dot{\pi }& =& -A\frac{\partial U}{\partial V}-dv+u\hfill \\ \hfill \dot{S}& =& \frac{d{v}^2}{T}\hfill \end{array}$$

(22)

with $v=\frac{\pi }{m}$ the velocity of the piston, A its area, d the damping constant, $T=\frac{\partial U}{\partial S}$ the temperature, and u the external force on the piston. The thermodynamic force F is identified as $\frac{dv}{T}$ and the thermodynamic flow as $J=v$. Clearly Onsager relations $J=LF$ are satisfied with $L=\frac{T}{d}$.

Example 3 

(Chemical reaction network). A third, more involved, example of irreversible thermodynamics are the dynamics of chemical reaction networks [34,41]. Consider an isolated (no incoming or outgoing chemical species, and no external heat flow) reaction network, with m chemical species and r reactions. Disregarding volume and pressure, consider the vector $x\in {\mathbb{R}}^m$ of concentrations of the chemical species. The dynamics take the form

$$\dot{x}=Nv\left(x\right),$$

(23)

where $v\in {\mathbb{R}}^r$ is the vector of reaction fluxes. The $m\times r$ stoichiometric matrixN, which consists of positive and negative integer elements, captures the structural balance laws of the reactions. Chemical reaction network theory, as originating from [42,43,44], identifies the edges of the underlying directed graph with the r reactions, and the nodes with the ccomplexes of the chemical reactions, i.e., the different left- and right-hand sides of the reactions in the network. This means that the stoichiometric matrix N is factorized as $N=ZB$, with B denoting the $c\times r$ incidence matrix of the graph of complexes, and Z denoting the $m\times c$ complex composition matrix (a matrix of nonnegative integers), whose $\rho$-th column captures the expression of the $\rho$-th complex in the m chemical species. It is shown in [45] that the dynamics $\dot{x}=Nv\left(x\right)$ of a large class of chemical reaction networks (including detailed-balanced mass action kinetics networks) can be written into the compact form

$$\dot{x}=Nv\left(x\right)=-Z\mathbb{L}\mathrm{Exp}\left(\frac{1}{RT}{Z}^{\top }\mu \left(x\right)\right),$$

(24)

where $\mathrm{Exp}$ is the vector exponential mapping $\mathrm{Exp}\left(z\right)={(exp{z}_1,\cdots ,exp{z}_c)}^{\top }$, R is the gas constant, T is the temperature, and $\mu$ is the m-dimensional vector of chemical potentials of the chemical species (for which e.g., in the case of detailed-balanced mass action kinetics explicit expressions are available). Furthermore, the matrix $\mathbb{L}:=B\mathcal{K}{B}^{\top }$ in (24) defines a weighted Laplacian matrix for the graph of complexes, with the diagonal elements ${\kappa }_1,\cdots ,{\kappa }_r$ of the diagonal matrix $\mathcal{K}$, depending on the temperature T and the reference state. We have the following fundamental property [45]

$${\gamma }^{\top }\mathbb{L}\mathrm{Exp}\gamma \ge 0\mathrm{for}\mathrm{all}\gamma \in {\mathbb{R}}^c,{\gamma }^{\top }\mathbb{L}\mathrm{Exp}\gamma =0\mathrm{iff}{B}^{\top }\gamma =0$$

(25)

Expressing the entropy S as a function of x and the total energy E, Gibbs’ fundamental relation yields$\frac{\partial S}{\partial x}(x,E)=-\frac{\mu }{T},\frac{\partial S}{\partial E}(x,E)=\frac{1}{T}$. This implies

$$\frac{d}{dt}S=\frac{1}{T}{\mu }^{\top }Z\mathbb{L}\mathrm{Exp}\left(\frac{Z^{\top }\mu }{RT}\right)=:\sigma \ge 0,$$

(26)

with equality if, and only if,$B^{\top }{Z}^{\top }\mu =N^{\top }\mu =0$, i.e., if and only if the chemical affinities$N^{\top }\mu$of the reactions are all zero. Hence the equilibria of the system correspond to states of minimal (i.e., zero) entropy production$\sigma$, in accordance with the theory of irreversible thermodynamics [38].

The vectors of thermodynamic forces F and thermodynamic flows J are given as

$$F=\frac{1}{T}{N}^{\top }\mu ,J=\mathcal{K}{B}^{\top }\mathrm{Exp}\frac{Z^{\top }\mu }{RT},$$

(27)

and indeed by (25)$\sigma =0$if and only if$F=0$. Note that Jcannot be expressed as a linear function of F and thus, in general, chemical reaction networks define nonlinear irreversible thermodynamics.

3.1. Quasi-Hamiltonian Formulation of Irreversible Thermodynamic Systems

Conservative mechanical systems are well-known to admit a Hamiltonian formulation. The same holds for many other physical systems. The Hamiltonian formulation of the dynamics of thermodynamic systems is, however, much more elusive. This has already studied and elaborated upon in, e.g., [16,24,41,46]. The present formulation emphasizes the factorization (15) of the irreversible entropy production.

Consider an isolated thermodynamic system with entropy S and energy E. Collect all other extensive variables in a vector denoted by z. The energy E can be expressed as a function $E=E(z,S)$ of z and S. Now consider the irreversible entropy production $\dot{S}=\sigma =J^{\top }F$, with J the vector of thermodynamic flows and F the vector of thermodynamic forces. Often (as illustrated by the examples to be discussed), the thermodynamic force F can be expressed as $C^{\top }\frac{\partial E}{\partial z}$ for some matrix C, whose elements are possibly depending on $\frac{\partial E}{\partial z},\frac{\partial E}{\partial S}$, as well as on $z,S$. Note that $\frac{\partial E}{\partial z}$ equals the vector of intensive variables associated with the extensive variables z, while the intensive variable $\frac{\partial E}{\partial S}$ equals the temperature T.

Energy conservation $\frac{d}{dt}E(z,S)=0$ together with $\frac{d}{dt}S(z,E)=J^{\top }F$ suggests writing the dynamics of z and S into the form

$$\left[\begin{array}{c}\dot{z}\\ \dot{S}\end{array}\right]=\underset{{\mathcal{J}}_e}{\underbrace{\left[\begin{array}{cc}\mathcal{J}& -CJ\\ J^{\top }{C}^{\top }& 0\end{array}\right]}}\left[\begin{array}{c}\frac{\partial E}{\partial z}\\ \frac{\partial E}{\partial S}\end{array}\right],$$

(28)

for some skew-symmetric matrix $\mathcal{J}$, possibly depending on $\frac{\partial E}{\partial z},\frac{\partial E}{\partial S}$ and $z,S$. This implies that the extended matrix ${\mathcal{J}}_e$ is also skew-symmetric, and thus indeed $\frac{d}{dt}E(z,S)=0$. Note however, that since the matrix ${\mathcal{J}}_e$ may depend on the intensive variables $\frac{\partial E}{\partial z},\frac{\partial E}{\partial S}$, it does not define a Poisson bracket on the state space with coordinates $z,S$. Therefore (28) will be called a quasi-Hamiltonian formulation.

This is illustrated by the previously discussed examples of the gas-piston-damper system, chemical reaction network, and heat exchanger as follows.

Example 4 

(Gas-piston-damper system continued). The dynamics of the gas-piston-damper system (22) can be written into the quasi-Hamiltonian form as (see also [41])

$$\left[\begin{array}{c}\dot{V}\\ \dot{\pi }\\ \dot{S}\end{array}\right]=\underset{{\mathcal{J}}_e}{\underbrace{\left[\begin{array}{ccc}0& A& 0\\ -A& 0& -d\frac{v}{T}\\ 0& d\frac{v}{T}& 0\end{array}\right]}}\left[\begin{array}{c}\frac{\partial E}{\partial V}\\ \frac{\partial E}{\partial \pi }\\ \frac{\partial E}{\partial S}\end{array}\right]+\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]u$$

(29)

with $v=\frac{\partial E}{\partial \pi }=\frac{\pi }{m}$ as the velocity of the piston and $T=\frac{\partial E}{\partial S}$ the temperature. The thermodynamic flow and force are $J=v$ and $F=\frac{d}{T}v$, respectively. Hence, ${\mathcal{J}}_e$ is of the form as given in (28) with $C^{\top }=\left[\begin{array}{cc}0& \frac{d}{T}\end{array}\right]$. ${\mathcal{J}}_e$ depends on the intensive variables T and v, therefore, it does not define a Poisson bracket.

Example 5 

(Chemical reaction network continued). In the case of chemical reaction networks, the vector of thermodynamic forces is given as $F=\frac{1}{T}{N}^{\top }\mu =C^{\top }\frac{\partial E}{\partial x}$ with $C=\frac{1}{T}N$, and $\frac{\partial E}{\partial x}=\mu$ the vector of chemical potentials. Furthermore, according to (27) the vector of thermodynamic flows is given as $J=\mathcal{K}{B}^{\top }\mathrm{Exp}\frac{Z^{\top }\mu }{RT}$. This leads to the quasi-Hamiltonian representation

$$\left[\begin{array}{c}\dot{x}\\ \dot{S}\end{array}\right]=\left[\begin{array}{cc}0& -\frac{1}{T}NJ\\ \frac{1}{T}{J}^{\top }{N}^{\top }& 0\end{array}\right]\left[\begin{array}{c}\frac{\partial E}{\partial x}\\ \frac{\partial E}{\partial S}\end{array}\right]$$

(30)

Example 6 

(Heat exchanger continued). The quasi-Hamiltonian formulation of the heat exchanger is slightly different. This caused by the fact that, in this example, we have two entropies, $S_1$ and $S_2$, corresponding to the two compartments (and not a total entropy as in the previous two examples). In fact, the quasi-Hamiltonian formulation of the heat exchanger is given as (see [41])

$$\left[\begin{array}{c}{\dot{S}}_1\\ {\dot{S}}_2\end{array}\right]=\left[\begin{array}{c}-\frac{q}{T_1}\\ \frac{q}{T_2}\end{array}\right]=\left[\begin{array}{cc}0& \lambda \left(\frac{1}{T_1}-\frac{1}{T_2}\right)\\ -\lambda \left(\frac{1}{T_1}-\frac{1}{T_2}\right)& 0\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]$$

(31)

since $\frac{\partial E}{\partial S_i}=T_i,i=1,2$, and the heat flow from compartment 1 to 2 is given by $q=\lambda (T_1-T_2)$. Here, we recognize $\frac{1}{T_1}-\frac{1}{T_2}$ as the thermodynamic force.

A further structured form of quasi-Hamiltonian modeling of irreversible thermodynamic systems, called irreversible port-Hamiltonian systems, was introduced in [24]; see [41,46,47] for more developments and references.

A special case occurs if the total energy $E(z,S)$ splits as

$$E(z,S)=H\left(z\right)+U\left(S\right),$$

(32)

for some thermal energy $U\left(S\right)$ and remaining energy $H\left(z\right)$. In this case, one obtains the equations

$$\dot{z}=\mathcal{J}\frac{\partial H}{\partial z}\left(z\right)-TCJ,\dot{S}=J^{\top }{C}^{\top }\frac{\partial H}{\partial z}\left(z\right),F=C^{\top }\frac{\partial H}{\partial z}\left(z\right)$$

(33)

If, furthermore, $J=LF$ with $L=L^{\top }$ (Onsager’s reciprocity relations) then the dynamical equations for the extensive variables z can be combined into

$$\dot{z}=\mathcal{J}\frac{\partial H}{\partial z}\left(z\right)-TCL{C}^{\top }\frac{\partial H}{\partial z}\left(z\right)=\left(\mathcal{J}-TCL{C}^{\top }\right)\frac{\partial H}{\partial z}\left(z\right)$$

(34)

This is the standard internal dynamics of a port-Hamiltonian system with state vector z; see e.g., [37,48,49]. In this case, irreversibility means that, even though the total energy $E(z,S)=H\left(z\right)+U\left(S\right)$ is preserved, the part of the energy given by $H\left(z\right)$ is continuously transformed (by the resistive power flow $T\dot{S}$) into the thermal energy $U\left(S\right)$. Conversely, one can show [31] that any port-Hamiltonian system can be embedded into an energy-conserving thermodynamic system.

Example 7 

(Mass-spring-damper system). A simple example is the ubiquitous mass-spring-damper system. Its dynamics are very similar to that of the gas-piston-damper system, the difference being that the internal energy $U(V,S)$ of the gas is replaced by the sum $\frac{1}{2}k{x}^2+U\left(S\right)$, where $\frac{1}{2}k{x}^2$ is the potential energy of the spring (with x denoting the elongation of the spring), and $U\left(S\right)$ is the thermal energy of the system. This leads to the dynamics (compare with (29))

$$\left[\begin{array}{c}\dot{x}\\ \dot{\pi }\\ \dot{S}\end{array}\right]=\left[\begin{array}{ccc}0& 1& 0\\ -1& 0& -d\frac{\pi }{mT}\\ 0& d\frac{\pi }{mT}& 0\end{array}\right]\left[\begin{array}{c}kx\\ \frac{\pi }{m}\\ T\end{array}\right]+\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]u,$$

(35)

as well as the following port-Hamiltonian formulation of the mass-spring-damper system

$$\left[\begin{array}{c}\dot{x}\\ \dot{\pi }\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -1& -d\end{array}\right]\left[\begin{array}{c}kx\\ \frac{\pi }{m}\end{array}\right]+\left[\begin{array}{c}0\\ 1\end{array}\right]u$$

(36)

3.2. Stability Analysis

The quasi-Hamiltonian formulation can be readily used for stability analysis. Note, however, that the conditions $\frac{\partial E}{\partial z}=0$, $\frac{\partial E}{\partial S}=0$ for E having a minimum often do not correspond to equilibria of interest. This is illustrated by the gas-piston-damper system, where these conditions correspond to pressure, velocity, and temperature all being equal to zero. Instead, in such cases it is of much more interest to consider the stability of steady states$(\overline{V},\overline{\pi }=0,\overline{S})$ corresponding to a non-zero force $\overline{u}$ delivered by the piston. In view of (29) and the energy expression $E(V,\pi ,S)=\frac{{\pi }^2}{2m}+U(V,S)$, this corresponds to the steady state condition

$$A\frac{\partial U}{\partial V}(\overline{V},\overline{S})=\overline{u}$$

(37)

(Note that $\dot{S}=0$ is ensured by $\overline{\pi }=0$ implying that $\overline{v}=\frac{\overline{\pi }}{m}=0$, and thus corresponds to a singularity in the skew-symmetric matrix ${\mathcal{J}}_e$, instead of a vanishing of all the partial derivatives of E. In particular, the temperature $T=\frac{\partial E}{\partial S}$ at steady state will not be zero.) Instead of using $E(z,S)$ as a candidate Lyapunov function which leads to the consideration of the availability function [50] (also called Bregman divergence or shifted Hamiltonian [37])

$$\widehat{E}(V,\pi ,S):=\frac{{\pi }^2}{2m}+U(V,S)+\overline{P}(V-\overline{V})-\overline{T}(S-\overline{S})-E(\overline{V},\overline{\pi },\overline{S}),$$

(38)

where $P=-\frac{\partial U}{\partial V}(\overline{S},\overline{V})$ and $\overline{T}=\frac{\partial U}{\partial S}(\overline{S},\overline{V})$ are the pressure and temperature at steady state, for some arbitrary value $\overline{S}$. Indeed, using the steady state condition (37), a direct computation yields

$$\frac{d}{dt}\widehat{E}=-\frac{\overline{T}}{T}d{v}^2\le 0$$

(39)

for all values of the temperature $T>0$ and the steady state temperature $\overline{T}>0$. Furthermore, given that for thermodynamic systems the internal energy $U(V,S)$ (and therefore $E(V,\pi ,S)$) is a convex function, $\widehat{E}(V,\pi ,S)$ is also convex with minimum at $\overline{V},\overline{\pi }=0,\overline{S}$. Hence, if $\widehat{E}(V,\pi ,S)$ is strictly convex (which is often the case), then this proves the asymptotic stability of the steady state. The use of the availability function for stability analysis and stabilization was already advocated for in [51]; see also, e.g., [47,52] for related work using the availability function in the context of passivity-based control of irreversible port-Hamiltonian systems.

This is extended to general quasi-Hamiltonian systems

$$\left[\begin{array}{c}\dot{z}\\ \dot{S}\end{array}\right]={\mathcal{J}}_e\left[\begin{array}{c}\frac{\partial E}{\partial z}(z,S)\\ \frac{\partial E}{\partial S}(z,S)\end{array}\right]+Gu,$$

(40)

where the skew–symmetric matrix ${\mathcal{J}}_e$ and the input matrix G may both depend on the extensive variables $z,S$ and the intensive variables $\frac{\partial E}{\partial z}(z,S),\frac{\partial E}{\partial S}(z,S)=T$. The steady state condition for $u=\overline{u}$ is given as

$${\overline{\mathcal{J}}}_e\left[\begin{array}{c}\frac{\partial E}{\partial z}(\overline{z},\overline{S})\\ \frac{\partial E}{\partial S}(\overline{z},\overline{S})\end{array}\right]+\overline{G}\overline{u}=0,$$

(41)

where ${\overline{\mathcal{J}}}_e$ and $\overline{G}$ denote the values of ${\mathcal{J}}_e$ and G at steady state, i.e.,

$${\overline{\mathcal{J}}}_e={\mathcal{J}}_e\left(\overline{z},\overline{S},\frac{\partial E}{\partial z}(\overline{z},\overline{S}),\frac{\partial E}{\partial S}(\overline{z},\overline{S})\right),\overline{G}=G\left(\overline{z},\overline{S},\frac{\partial E}{\partial z}(\overline{z},\overline{S},\frac{\partial E}{\partial S}(\overline{z},\overline{S})\right)$$

(42)

Assuming the energy function $E(z,S)$ to be convex (which is normally the case in thermodynamic systems), then the availability function is given as the convex function

$$\widehat{E}(z,S):=E(z,S)-\frac{\partial E}{\partial z^{\top }}(\overline{z},\overline{S})(z-\overline{z})-\frac{\partial E}{\partial S}(\overline{z},\overline{S})(S-\overline{S})-E(\overline{z},\overline{S}),$$

(43)

having a minimum at $(\overline{z},\overline{S})$. A key property of the availability function $\widehat{E}$ is that

$$\nabla \widehat{E}(z,S)=\nabla E(z,S)-\nabla E(\overline{z},\overline{S})$$

(44)

where $\nabla \widehat{E}(z,S)$ denotes the gradient vector of $\widehat{E}$ (written as a column vector). The computation of $\frac{d}{dt}\widehat{E}(z,S)$ yields, exploiting the steady state condition (41),

$$\frac{d}{dt}\widehat{E}(z,S)={(\nabla \widehat{E}(z,S))}^{\top }\left[{\mathcal{J}}_e\nabla E(z,S)-{\overline{\mathcal{J}}}_e\nabla E(\overline{z},\overline{S})\right]$$

(45)

It follows that $\frac{d}{dt}\widehat{E}(z,S)\le 0$ if, and only if

$${\left(\nabla E(z,S)-\nabla E(\overline{z},\overline{S})\right)}^{\top }\left[{\mathcal{J}}_e\nabla E(z,S)-{\overline{\mathcal{J}}}_e\nabla E(\overline{z},\overline{S})\right]\le 0$$

(46)

(This condition is similar to the condition for asymptotic stability of steady states of port-Hamiltonian systems as derived in [53]; see also [37].) Hence if (46) is satisfied and $\widehat{E}(z,S)$ is not only convex but even strictly convex, then $\widehat{E}(z,S)$ serves as a Lyapunov function for assessing the (asymptotic) stability of the steady state $(\overline{z},\overline{S})$.

Instead of expressing the energy E as a function $E(z,S)$ of the remaining extensive variables $z,S$ and writing the dynamics as a quasi-Hamiltonian system with Hamiltonian given by E, one may also write the entropy S as a function $S(z,E)$ and try to express the dynamics as being generated by the gradient of this entropy function. However, since (in the isolated case) $\frac{d}{dt}S\ge 0$, this constitutes quite a different scenario. An example where it is possible is a chemical reaction network as mentioned before. Instead of the quasi-Hamiltonian formulation (30), one rewrites the dynamics as (with z replaced by the vector of concentrations x)

$$\left[\begin{array}{c}\dot{x}\\ \dot{E}\end{array}\right]=\underset{{\mathcal{F}}_e}{\underbrace{T\left[\begin{array}{cc}0& -NJ\\ J^{\top }{N}^{\top }& T{J}^{\top }F\end{array}\right]}}\left[\begin{array}{c}\frac{\partial S}{\partial x}(x,E)\\ \frac{\partial S}{\partial E}(x,E)\end{array}\right],$$

(47)

where $\frac{\partial S}{\partial x}(x,E)=-\frac{\mu }{T}$ and $\frac{\partial S}{\partial E}(x,E)=\frac{1}{T}$. Consequently

$$\dot{E}=-J^{\top }{N}^{\top }\mu +T{J}^{\top }F=0,$$

(48)

because of $F=\frac{1}{T}{N}^{\top }\mu$. Since $\frac{d}{dt}S\ge 0$, the availability function corresponding to $V\left(z\right):=-S(x,\overline{E})$, with $\overline{E}$ the constant total energy of the system, can be used as a Lyapunov function for stability analysis; cf. [34] for details. Note that the matrix ${\mathcal{F}}_e$, in the right-hand of (47), is not a skew-symmetric matrix anymore. In fact, the formulation (47) resembles the formulation of thermodynamic systems as used in the GENERIC formalism; see, e.g., [30].

Another interesting case are isothermal chemical reaction networks. In this case [45] one considers the Gibbs free energy (Legendre transform of $E(z,S)$ with respect to S)

$$G(z,T)=E(z,S)-TS,T=\frac{\partial E}{\partial S}$$

(49)

for constant T. By the properties of the Legendre transform $\frac{\partial G}{\partial z}=\frac{\partial E}{\partial z}=\mu$ (the vector of chemical potentials). Hence, in view of (47) one obtains for constant T

$$\frac{d}{dt}G={\mu }^{\top }\dot{z}=-{\mu }^{\top }NJ=T{F}^{\top }J=-T\sigma$$

(50)

with $\sigma$ the irreversible entropy production. Alternatively expressed, whenever the temperature T is constant, $\frac{d}{dt}G=\frac{d}{dt}E-T\frac{d}{dt}S$, while $\frac{d}{dt}E=q$ (with q the heat flow needed to keep the temperature constant) and $\frac{d}{dt}S=q-\sigma$. Taken together this indeed yields $\frac{d}{dt}G=-T\sigma$.

4. Thermodynamic Phase Space and Liouville Geometry

A typical feature of thermodynamic systems modeling is the use of many more variables than the minimal number of variables to describe the ‘state’ of the system. For example, a simple thermodynamic system is described by a 2-dimensional submanifold L of the 3-dimensional space of macroscopic quantities $V,P,T$; one extensive, and two intensive. Then, based on the first and second laws of thermodynamics, two extra extensive variables $E,S$ are introduced. As a result, the system is described as a 2-dimensional submanifold L of ${\mathbb{R}}^5$; the thermodynamic phase space is generated by the three extensive variables $E,S,V$, and the two intensive variables $T,P$. This is immediately extended to the general thermodynamic case, where L is an n-dimensional submanifold of the $(2n+1)$-dimensional thermodynamic phase space (comprising $n+1$ extensive variables and n intensive variables).

4.1. Constitutive Relations and Liouville Submanifolds

The Legendre submanifold L only defines the constitutive relations of the thermodynamic system, i.e., its thermostatics. The first and second laws impose constraints on any possible dynamics of the thermodynamic system, but do not define it. On the other hand, two requirements for any dynamics on the full thermodynamic phase space are natural: (1) the dynamics should respect the ‘structure’ of the thermodynamic phase space, (2) it should respect the constitutive relations; i.e., should leave the submanifold L invariant. The proper geometry to address this is contact geometry. However, in order to unify the energy and entropy representation we will take one more abstraction step; from contact geometry to Liouville geometry. This will have the additional benefit of making a clear separation between extensive and intensive variables, and of being computationally more easy; see [25,26,54,55] for further details and ramifications.

For a simple thermodynamic system with extensive variables $E,S,V$ and intensive variables $T,-P$, the step from contact to Liouville geometry amounts to replacing the intensive variables $T,-P$ (in the energy representation) with their homogeneous coordinates$p_E,p_S,p_V$ with $p_E\ne 0$, i.e.,

$$T=\frac{p_S}{-p_E},-P=\frac{p_V}{-p_E},$$

(51)

and thereby to express the intensive variables $\frac{1}{T},\frac{P}{T}$ in the entropy representation as

$$\frac{1}{T}=\frac{p_E}{-p_S},\frac{P}{T}=\frac{p_V}{-p_S}$$

(52)

In this way, the two Gibbs one-forms $dE-TdS+PdV$ and $dS-\frac{1}{T}dE-\frac{P}{T}dV$ are replaced by a single symmetric expression, namely by the Liouville one-form

$$p_{E}dE+p_{S}dS+p_{V}dV,$$

(53)

on the cotangent bundle $T^{*}{\mathbb{R}}^3$, with ${\mathbb{R}}^3$ the space of extensive variables $E,S,V$. By definition of homogeneous coordinates, the vector $(p_E,p_S,p_V)$ is always different from the 0-vector. Hence, the space $\left\{(E,S,V,p_E,p_S,p_V)\right\}$ is actually the cotangent bundle $T^{*}{\mathbb{R}}^3$ minus its zero section. Using homogeneous coordinates, the 2-dimensional Legendre submanifold L is now replaced by the 3-dimensional submanifold $\mathcal{L}\subset T^{*}{\mathbb{R}}^3$, given as

$$\mathcal{L}=\{(E,S,V,p_E,p_S,p_V)\mid (E,S,V,\frac{p_S}{-p_E},\frac{p_V}{-p_E})\in L,(p_E,p_S,p_V)\ne 0\}$$

(54)

It turns out that $\mathcal{L}$ is a Lagrangian submanifold [56,57,58], which is moreover homogeneous, in the sense that whenever $(E,S,V,p_E,p_S,p_V)\in \mathcal{L}$ then also $(E,S,V,\lambda p_E,\lambda p_S,\lambda p_V)\in \mathcal{L}$, for any non-zero $\lambda \in \mathbb{R}$. Such submanifolds are fully characterized as maximal manifolds restricted to which the Liouville form $p_{E}dE+p_{S}dS+p_{V}dV$ is zero, and are therefore called homogeneous Lagrangian submanifolds or Liouville submanifolds [25].

In general, one considers the $(n+1)$-dimensional manifold $\mathcal{Q}$ of all the extensive variables (including E and S), and its $(2n+2)$-dimensional cotangent bundle without zero section, denoted by ${\mathcal{T}}^{*}Q$. The constitutive properties of the thermodynamic system are defined by a $(n+1)$-dimensional Liouville submanifold $\mathcal{L}$. Conversely, starting from ${\mathcal{T}}^{*}\mathcal{Q}$ we may define a contact manifold in the following way [57]. For each $q\in \mathcal{Q}$ and cotangent space $T_q^{*}\mathcal{Q}$ consider the projective space$\mathbb{P}\left(T_q^{*}\mathcal{Q}\right)$, given as the set of rays in $T_q^{*}\mathcal{Q}$, that is, all the non-zero multiples of a non-zero cotangent vector. Thus, the projective space $\mathbb{P}\left(T_q^{*}\mathcal{Q}\right)$ has dimension n, and there is the canonical projection ${\pi }_q:{\mathcal{T}}_q^{*}\mathcal{Q}\to \mathbb{P}\left(T_q^{*}\mathcal{Q}\right)$, where ${\mathcal{T}}_q^{*}\mathcal{Q}$ denotes the cotangent space without its zero vector. The fiber bundle of the projective spaces $\mathbb{P}\left(T_q^{*}\mathcal{Q}\right)$, $q\in \mathcal{Q},$ over the base manifold $\mathcal{Q}$ will be denoted by $\mathbb{P}\left(T^{*}\mathcal{Q}\right)$. Furthermore, the bundle projection obtained by considering ${\pi }_q:{\mathcal{T}}_q^{*}\mathcal{Q}\to \mathbb{P}\left(T_q^{*}\mathcal{Q}\right)$ for every $q\in \mathcal{Q}$ is denoted by $\pi :{\mathcal{T}}^{*}\mathcal{Q}\to \mathbb{P}\left(T^{*}\mathcal{Q}\right)$. As detailed in [26,57,58], $\mathbb{P}\left(T^{*}\mathcal{Q}\right)$ defines a contact manifold of dimension $2n+1$, and will serve as the canonical thermodynamic phase space for the thermodynamic system with space of external variables $\mathcal{Q}$. In the case of a simple thermodynamic system the bundle projection $\pi$ is given in coordinates as

$$(p_E,p_S,p_V)\mapsto (\frac{p_S}{-p_E}=T,\frac{p_V}{-p_E}=-P)$$

(55)

whenever $p_E\ne 0$ (energy representation), or as

$$(p_E,p_S,p_V)\mapsto (\frac{p_E}{-p_S}=-\frac{1}{T},\frac{p_V}{-p_S}=\frac{P}{T})$$

(56)

whenever $p_S\ne 0$ (entropy representation). The cotangent bundle $T^{*}\mathcal{Q}$, and therefore also ${\mathcal{T}}^{*}\mathcal{Q}$, are endowed with the natural one-form $\alpha$, called the Liouville form, which in natural cotangent bundle coordinates $(q,p)=(q_1,\cdots ,q_{n+1},p_1,\cdots ,p_{n+1})$ is given as $\alpha ={\sum }_{i=1}^{n+1}{p}_{i}d{q}_i$. A submanifold $\mathcal{L}\subset {\mathcal{T}}^{*}\mathcal{Q}$ is a Liouville submanifold if $\alpha$ restricted to $\mathcal{L}$ is zero, and furthermore $\mathcal{L}$ is maximal with respect to this property. It turns out that maximality of $\mathcal{L}$ is equivalent to $dim\mathcal{L}=n+1$.

For any point in a Liouville submanifold $\mathcal{L}$ there exists a neighborhood of this point, a splitting of the index set $\{1,\cdots ,n+1\}=I\cup J$ and a function $F(q_I,p_J)$ that is homogeneous of degree 1 in $p_J$ (in particular $J\ne \varnothing$), where $q_I$ denotes the vector of coordinates $q_i$ with $i\in I$ and $p_J$ denotes the vector of coordinates $p_j$ with $j\in J$, such that on this neighborhood

$$\mathcal{L}=\{(q,p)\in {\mathcal{T}}^{*}\mathcal{Q}\mid p_i=\frac{\partial F}{\partial q_i},i\in I,q_j=-\frac{\partial F}{\partial p_j},j\in J\}$$

(57)

By homogeneity of F in $p_J$ it follows that for any $j\in J$ we can write, whenever $p_j\ne 0$,

$$F(q_I,p_J)=-p_j\widehat{F}(q_I,\frac{p_{\ell }}{-p_j}),\ell \in J,\ell \ne j)$$

(58)

for some function $\widehat{F}$. The choice of $j\in J$ corresponds to a choice of the coordinates for the contact manifold $\mathbb{P}\left(T^{*}\mathcal{Q}\right)$ (for example, corresponding to the energy or the entropy representation).

The Liouville submanifold $\mathcal{L}$ projects under $\pi$ to a Legendre submanifold $L\subset \mathbb{P}\left(T^{*}\mathcal{Q}\right)$, and conversely any Legendre submanifold $L\subset \mathbb{P}\left(T^{*}\mathcal{Q}\right)$ is the projection of a Liouville submanifold $\mathcal{L}\subset {\mathcal{T}}^{*}Q$. Furthermore, the function $\widehat{F}$ serves as a generating function for the Legendre submanifold L. Although the close relation of contact geometry with the Liouville geometry of cotangent bundles is well-known in differential geometry [57,58], the use of homogeneous coordinates for thermodynamics was first advocated for in [29].

4.2. Homogeneous Hamiltonian Dynamics and Port-Thermodynamic Systems

The dynamics of the thermodynamic system should now satisfy the following two basic requirements. First, it should respect the structure of ${\mathcal{T}}^{*}Q$, and therefore of the contact manifold $\mathbb{P}\left(T^{*}\mathcal{Q}\right)$. Second, it should leave invariant the Liouville submanifold $\mathcal{L}$ specifying the constitutive relations. The first requirement amounts to a requirement that the dynamics are Hamiltonian on ${\mathcal{T}}^{*}Q$, with the additional property that the Liouville forms on ${\mathcal{T}}^{*}Q$ is preserved. This can be seen to correspond to Hamiltonian dynamics with a Hamiltonian K that is homogeneous of degree 1 in the p-variables. Thus, if q are coordinates for Q and $(q,p)$ are corresponding cotangent bundle coordinates for ${\mathcal{T}}^{*}Q$ (such that the Liouville form is $\alpha ={\sum }_{i=1}^{n+1}{p}_{i}d{q}_i$), then we consider Hamiltonians $K(q,p)$ satisfying $K(q,\lambda p)=\lambda K(q,p)$ for all $\lambda \in \mathbb{R}$ different from zero. Equivalently, by Euler’s theorem, the Hamiltonian $K(q,p)$ should satisfy

$$K(q,p)=\sum _{i=1}^{n+1}{p}_i\frac{\partial K}{\partial p_i}(q,p),$$

(59)

with the functions $\frac{\partial K}{\partial p_i}(q,p)$ homogeneous of degree 0 in p, $i=1,\cdots ,n+1$. The second requirement is equivalent to K being such that K restricted to $\mathcal{L}$ is zero. Finally, we will split K into two parts, i.e.,

$$K^a+K^{c}u,u\in {\mathbb{R}}^m$$

(60)

Here $K^a:{\mathcal{T}}^{*}\mathcal{Z}\to \mathbb{R}$ is the Hamiltonian corresponding to the autonomous dynamics due to internal non-equilibrium conditions, while $K^c=(K_1^c,\cdots ,K_m^c)$ is a row vector of Hamiltonians corresponding to dynamics arising from interaction with the surroundings of the system, parameterized by a vector u of control or input variables. (However, as we will notice in the context of the damper system (81), there are cases where the dependence on u is not affine.) Thus all these Hamiltonians $K^a,K_1^c,\cdots ,K_m^c$ are zero on $\mathcal{L}$. This implies that the dynamics of the extensive variables are given as

$${\dot{q}}_i=\frac{\partial K^a}{\partial p_i}(q,p)+\sum _{j=1}^m\frac{\partial K^{c_j}}{\partial p_i}(q,p),i=1,\cdots ,n+1,$$

(61)

where the right-hand side is homogeneous of degree 0 in p. (Note that this does not mean that the right-hand side is necessarily independent of p; it may depend on degree 0 variables $\frac{p_i}{-p_E}$; i.e., on the intensive variables!)

Finally, there are two important extra constraints on $K^a$ (the Hamiltonian governing the autonomous dynamics corresponding to $u=0$) which are directly imposed by the first and second laws. By the first law $0=\dot{E}=\frac{\partial K^a}{\partial p_E}$ on $\mathcal{L}$. Furthermore, by the second law necessarily $\frac{\partial K^a}{\partial p_S}{|}_{\mathcal{L}}\ge 0$. In fact, using the postulate of factorization of the irreversible entropy production as discussed in Section 3 one has

$$\frac{\partial K^a}{\partial p_S}{|}_{\mathcal{L}}=\sigma =\sum _{k=1}^s{F}_k{J}_k\ge 0$$

(62)

where $\sigma =0$ if and only if $F_k=0,k=1,\cdots ,s$.

Such constraints do not hold for the control (interaction) Hamiltonians $K^c$. In fact, the corresponding terms of the control Hamiltonians define natural outputs conjugated to the inputs u. The first option is to define the m-dimensional row vector

$$y_p=\frac{\partial K^c}{\partial p_E},$$

(63)

with the subscript p in $y_p$ standing for power. Then, it follows that $\frac{d}{dt}E=y_{p}u$, and thus, $y_p$ is the vector of power-conjugate (passive) outputs corresponding to the input vector u. Similarly, by defining the m-dimensional row vector

$$y_e=\frac{\partial K^c}{\partial p_S}$$

(64)

it follows that $\frac{d}{dt}S\ge y_{e}u$. Hence $y_e$ is the output vector which is conjugate to u in terms of entropy flow. This is summarized in the following definition of a port-thermodynamic system as given in [25,26].

Definition 1. 

Consider a manifold $\mathcal{Q}$ of extensive variables. A port-thermodynamic system on $\mathcal{Q}$ is defined by a Liouville submanifold $\mathcal{L}\subset {\mathcal{T}}^{*}\mathcal{Q}$ specifying the constitutive relations of the system, together with a Hamiltonian $K^a+K^{c}u,u\in {\mathbb{R}}^m$, homogeneous of degree 1 in p, which is zero on $\mathcal{L}$ for every u and satisfying $\frac{\partial K^a}{\partial p_E}=0$ on $\mathcal{L}$, and $\frac{\partial K^a}{\partial p_S}\ge 0$ on $\mathcal{L}$. Its power port is defined by u together with the output $y_p=\frac{\partial K^c}{\partial p_E}$, and its entropy flow port by u and $y_e=\frac{\partial K^c}{\partial p_S}$.

Remark 3. 

The Hamiltonian K generates the dynamics, but does not have an interpretation in terms of energy. In fact, K is dimensionless; see [25,26] for further information.

Remark 4. 

Through the use of entropy flow ports, one could express the irreversible entropy production $\sigma ={\sum }_{k=1}^s{F}_k{J}_k\ge 0$ as being the result of the interconnection of the system with a pure entropy producing element. This is especially clear if the vector of thermodynamic flows J can be expressed as a function of the vector of thermodynamic forces F, like in Onsager’s relations $J=LF,L=L^{\top }$. Namely, in this case one may consider $u=J$ and entropy conjugate outputs $y_e=F$, and then ‘close’ the loop by setting $u=L{y}_e$.

Because of the homogeneity of the Liouville submanifold $\mathcal{L}$ and of the Hamiltonian K, the port-thermodynamic system defined on ${\mathcal{T}}^{*}\mathcal{Q}$, including its power and entropy flow ports, projects to a system on the thermodynamic phase space $\mathbb{P}\left(T^{*}\mathcal{Q}\right)$ with Legendre submanifold L; see [25] for details. In fact, the resulting class of systems on $\mathbb{P}\left(T^{*}\mathcal{Q}\right)$ is very close to the classes of input–output contact systems and conservative control contact systems on contact manifolds as introduced and studied in [7,8,9,59,60].

The definition of port-thermodynamic systems (Definition 1) is illustrated by the examples of the gas-piston-damper system, chemical reaction network, and heat exchanger as follows.

Example 8 

(Gas-piston-damper system continued [26]). As discussed before, the extensive variables are E (energy), S (entropy), V (volume), and $\pi$ (momentum of the piston). For simplicity we will take $A=1$. The constitutive properties are given by the Liouville submanifold

$$\begin{array}{ccc}\hfill \mathcal{L}& =& \{(E,S,V,\pi ,p_E,p_S,p_V,p_{\pi })\mid E(S,V,\pi )=U(S,V)+\frac{{\pi }^2}{2m},\hfill \\ & & p_S=-p_E\frac{\partial U}{\partial S},p_V=-p_E\frac{\partial U}{\partial V},p_{\pi }=-p_E\frac{\pi }{m}\}\hfill \end{array}$$

(65)

with generating function $-p_E\left(U(S,V)+\frac{{\pi }^2}{2m}\right)$. The dynamics are given by the Hamiltonian (homogeneous of degree 1 in p)

$$K=p_V\frac{\pi }{m}+p_{\pi }\left(-\frac{\partial U}{\partial V}-d\frac{\pi }{m}\right)+p_S\frac{d{\left(\frac{\pi }{m}\right)}^2}{\frac{\partial U}{\partial S}}+\left(p_{\pi }+p_E\frac{\pi }{m}\right)u,$$

(66)

which obviously is zero on $\mathcal{L}$. The power-conjugate output $y_p=\frac{\pi }{m}$ is the velocity of the piston. One could also add an extra control Hamiltonian $\left(\frac{p_S}{T}+p_E\right)v$, where $T=\frac{\partial U}{\partial S}$ is the temperature, and v is the heat flow from an external heat source into the cylinder. The corresponding entropy conjugate output is $y_e=\frac{1}{T}$.

Example 9 

(Chemical reaction network continued [33]). Consider a chemical reaction network in entropy representation, with the entropy S represented as a function $S=S(E,x)$ of the vector of chemical concentrations x and energy E. Then the Liouville submanifold describing the state properties of the reaction network is given as

$$\mathcal{L}=\{(x,S,E,p_x,p_S,p_E)\mid S=S(E,x),p_x=-p_S\frac{\partial S}{\partial x}(E,x),p_E=-p_S\frac{\partial S}{\partial E}(E,x)\}$$

(67)

with generating function $-p_{S}S(x,E)$ and $\frac{\partial S}{\partial x}(E,x)=-\frac{\mu }{T},\frac{\partial S}{\partial E}(E,x)=\frac{1}{T}$. The internal dynamics of the chemical reaction network are generated by the Hamiltonian

$$K^a=-p_x^{\top }Z\mathbb{L}\mathrm{Exp}\frac{-Z^{\top }}{R}\frac{\partial S}{\partial x}(E,x)-p_S\frac{\partial S}{\partial x^{\top }}(E,x)ZL\mathrm{Exp}\frac{-Z^{\top }}{R}\frac{\partial S}{\partial x}(E,x)$$

(68)

Furthermore, the control Hamiltonian

$$K^c=p_S\frac{\partial S}{\partial E}(E,x)+p_E,$$

(69)

corresponds to a heat flow input, and an entropy flow conjugate output$y_e=\frac{\partial S}{\partial E}{(x,E)}_{|\mathcal{L}}$equal to the reciprocal temperature. Another possible choice is

$$K^c=p_S\frac{\partial S}{\partial x_i}(E,x)+p_{xi},$$

(70)

corresponding to material in/outflow of the i-th chemical species, with entropy flow conjugate output$y_e=\frac{\partial S}{\partial x_i}{(E,x)}_{|\mathcal{L}}$equal to the chemical potential${\mu }_i$of the i-th chemical species divided by the temperature T.

Example 10 

(Heat exchanger continued [26]). The extensive variables are $S_1,S_2$ (entropies of the two compartments), and E (total internal energy). The state properties are described by the Liouville submanifold

$$\mathcal{L}=\{(S_1,S_2,E,p_{S_1},p_{S_2},p_E)\mid E=E_1\left(S_1\right)+E_2\left(S_2\right),p_{S_1}=-p_E{E}_1^{\prime }\left(S_1\right),p_{S_2}=-p_E{E}_2^{\prime }\left(S_2\right)\},$$

(71)

corresponding to the generating function $-p_E\left(E_1\left(S_1\right)+E_2\left(S_2\right)\right)$, with $E_1,E_2$ as the internal energies of the two compartments. Denoting the temperatures $T_1=E_1^{\prime }\left(S_1\right),$$T_2=E_2^{\prime }\left(S_2\right)$, the internal dynamics corresponding to Fourier’s law is given by the Hamiltonian

$$K^a=\lambda (\frac{1}{T_1}-\frac{1}{T_2})(p_{S_1}{T}_2-p_{S_2}{T}_1),$$

(72)

with $\lambda$ Fourier’s conduction coefficient.

5. Control by Interconnection

Control by interconnection is the paradigm of controlling a system by interconnecting it (through its inputs and outputs) to an additional controller system. The aim is to influence the dynamics of the original system by shaping the dynamics of the interconnected system by a proper choice of the controller system. Applied to port-thermodynamic systems, this means that given a plant thermodynamic system we interconnect it to a controller port-thermodynamic system such that in the closed-loop port-thermodynamic system the plant states converge to the desired set-point values. Port-thermodynamic systems can be interconnected, either by their power ports or by their entropy flow ports; cf. [26] for details. For example, the power port interconnection of two systems with variables

$$(E_i,S_i,q_i,p_{E_i},p_{S_i},p_i)\in {\mathcal{T}}^{*}{\mathcal{Q}}_i,i=1,2,$$

(73)

is defined as follows. With the homogeneity assumption in p in mind, impose the following constraint

$$p_{E_1}=p_{E_2}=:p_E$$

(74)

This leads to the summation of the Liouville one-forms ${\alpha }_1$ and ${\alpha }_2$ given by

$${\alpha }_{\mathrm{sum}}:=p_{E}d(E_1+E_2)+p_{S_1}d{S}_1+p_{S_2}d{S}_2+p_{1}d{q}_1+p_{2}d{q}_2$$

(75)

on the composed space ${\mathcal{T}}^{*}{\mathcal{Q}}_1\circ {\mathcal{T}}^{*}{\mathcal{Q}}_2$ defined as

$${\mathcal{T}}^{*}{\mathcal{Q}}_1\circ {\mathcal{T}}^{*}{\mathcal{Q}}_2:=\left\{(E,S_1,S_2,q_1,q_2,p_E,p_{S_1},p_{S_2},p_1,p_2)\right\}$$

(76)

Let the constitutive relations of the two port-thermodynamic systems be defined by the Liouville submanifolds ${\mathcal{L}}_i\subset {\mathcal{T}}^{*}{\mathcal{Q}}_i,i=1,2$. Then, the constitutive relations of the interconnected system are defined by the composition

$$\begin{array}{ccc}\hfill {\mathcal{L}}_1\circ {\mathcal{L}}_2& :=& \{(E,S_1,S_2,q_1,q_2,p_E,p_{S_1},p_{S_2},p_1,p_2)\mid \hfill \\ & & E=E_1+E_2,(E_i,S_i,q_i,p_{E_i},p_{S_i},p_i)\in {\mathcal{L}}_i,i=1,2\}\hfill \end{array}$$

(77)

Furthermore, consider the dynamics on ${\mathcal{L}}_i$ defined by Hamiltonians $K_i=K_i^a+K_i^c{u}_i,i=1,2$, where $K_i^c$ is the row vector of control Hamiltonians of system $i,i=1,2$. Additionally assume that the functions $K_i$ do not depend on the energy variables $E_i,i=1,2$. Then $K_1+K_2$ is well-defined on ${\mathcal{L}}_1\circ {\mathcal{L}}_2$ for all $u_1,u_2$. Next, consider the power conjugate outputs $y_{p1},y_{p2}$. By imposing interconnection constraints on the power port variables $u_1,u_2,y_{p1},y_{p2}$ satisfying the power preservation property

$$y_{p1}^{\top }{u}_1+y_{p2}^{\top }{u}_2=0$$

(78)

one obtains an interconnected port-thermodynamic system with constitutive relations described by ${\mathcal{L}}_1\circ {\mathcal{L}}_2$. Similarly, interconnecting the inputs $u_1,u_2$ to the entropy flow outputs $y_{e1},y_{e2}$ in such a way that

$$y_{e1}^{\top }{u}_1+y_{e2}^{\top }{u}_2\ge 0,$$

(79)

leads again to a port-thermodynamic system.

A basic control problem concerns the stabilization of a system to a desired set-point value (regulation). How can we use control by interconnection to address this problem? Suppose we want to stabilize the system at some set-point value $(z^{*},S^{*})$. If $E(z,S)$ already has a strict minimum at $(z^{*},S^{*})$ then one may asymptotically stabilize $(z^{*},S^{*})$ by the interconnection with a damper system [32]. In fact, assume for simplicity of exposition that $m=1$ (scalar output $y_p$). Then, consider an additional linear damper system (cf. [26]), whose Liouville submanifold is given as

$${\mathcal{L}}_d=\{(U_d,S_d)\mid U_d=U_d\left(S_d\right),p_{S_d}=-p_{U_d}{U}_d^{\prime }\left(S_d\right)\},$$

(80)

with entropy as $S_d$, internal energy as $U_d\left(S_d\right)$, and $U_d^{\prime }\left(S_d\right)$ its temperature. The dynamics of this damper system are generated by the Hamiltonian (see [26])

$$K=(p_{U_d}+p_{S_d}\frac{1}{U_d^{\prime }\left(S_d\right)})d{u}_d^2$$

(81)

(note the quadratic dependence on the input $u_d$), with power conjugate output $y_d=d{u}_d$ (damping force). Then interconnect the plant port-thermodynamic system $(\mathcal{L},K=K^a+K^{c}u)$ to this damper system by setting

$$u=-y_d,u_d=y$$

(82)

This results (after setting $p_{U_d}=p_E$) in the interconnected port-thermodynamic system with total Hamiltonian given as

$$K^a(E,S,z,p_E,p_S,p_z)-K^c(E,S,z,p_E,p_S,p_z)dy+(p_{U_d}+p_{S_d}\frac{1}{U^{\prime }\left(S_d\right)})d{y}^2$$

(83)

with total energy $E(S,z)+U_d\left(S_d\right)$. This implies

$$\frac{d}{dt}E(S,q)=-\frac{d}{dt}{U}_d\left(S_d\right)=-U_d^{\prime }\left(S_d\right)\dot{S_d}=-d{y}^2\le 0$$

(84)

Hence, by an application of LaSalle’s invariance principle, the system converges to the largest invariant set within the set where the power conjugate output $y_p$ is zero. Note that $y=0$ corresponds to zero entropy production $\dot{S_d}=0$; in accordance with irreversible thermodynamics. If the largest invariant set where y is zero equals the singleton $(E^{*},S^{*},q^{*})$ then asymptotic stability of $(E^{*},S^{*},q^{*})$ results; for some limiting value $S_d^{*}$ of the entropy $S_d$ of the damper system; see also [32].

What can be done if $E(S,z)$ does not have a strict minimum at $(S^{*},z^{*})$ ? This can be approached via the (generalized) Energy-Casimir method; similar to the theory of control by interconnection for port-Hamiltonian systems, see e.g., [37,61]. Consider a port-thermodynamic system with Liouville submanifold $\mathcal{L}\subset {\mathcal{T}}^{*}\mathcal{Q}$ with the generating function (in the energy representation) $-p_{E}E(S,z)$. A classical tool in the stability analysis of ordinary Hamiltonian dynamics is to consider additional conserved quantities; see e.g., [56,57,58]. In order to extend this idea to the present case, let us strengthen our assumption on $K^a$ by requiring that $\frac{\partial K^a}{\partial p_E}=0$everywhere on ${\mathcal{T}}^{*}\mathcal{Q}$; i.e., not just on $\mathcal{L}$. Next, consider additional conserved quantities for the dynamics $X_{K^a}$ that are only depending on the extensive variables $S,z$; i.e., functions $C(S,z)$ such that $\{K^a,C\}=0$, where $\{·,·\}$ is the standard Poisson bracket on ${\mathcal{T}}^{*}\mathcal{Q}$. Hence, also $\{K^a,E+C\}=0$. Subsequently, note that [32]

$$\alpha =p_{E}dE+p_{S}dS+p_{z}dz=p_{E}d(E+C(S,z))+(p_S-p_E\frac{\partial C}{\partial S})dS+(p_z-p_E\frac{\partial C}{\partial z})dz$$

(85)

Hence, the transformation

$$(E,S,z,p_E,p_S,p_z)\mapsto (E+C,S,z,p_E,p_S-p_E\frac{\partial C}{\partial S},p_z-p_E\frac{\partial C}{\partial z})=:(\tilde{E},\tilde{S},\tilde{z},{\tilde{p}}_E,{\tilde{p}}_S,{\tilde{p}}_z)$$

(86)

is a point transformation (that is, leaving the Liouville form invariant). Note that in the new coordinates the intensive variables $\frac{p_S}{-p_E},\frac{p_z}{-p_E}$ are transformed into the new intensive variables

$$\begin{array}{c}\frac{{\tilde{p}}_S}{-p_E}=\frac{p_S-p_E\frac{\partial C}{\partial S}}{-p_E}=\frac{p_S}{-p_E}+\frac{\partial C}{\partial S}\hfill \\ \frac{{\tilde{p}}_z}{-p_E}=\frac{p_z-p_E\frac{\partial C}{\partial z}}{-p_E}=\frac{p}{-p_E}+\frac{\partial C}{\partial z}\hfill \end{array}$$

(87)

In these new coordinates the generating function for $\mathcal{L}$ in entropy representation is given by $\tilde{E}(S,z)=E(S,z)+C(S,z)$. Furthermore, since $\{K^a,E+C\}=0$, the transformed Hamiltonian

$$\widetilde{K^a}(\tilde{E},\tilde{S},\tilde{q},{\tilde{p}}_E,{\tilde{p}}_S,\tilde{p}):=K^a(E,S,q,p_E,p_S,p)$$

(88)

satisfies $\{\widetilde{K^a},\tilde{E}\}=0$. Hence, in the new coordinates we are back to the situation considered before: if $\tilde{E}(S,z)$ has a strict minimum at $(S^{*},z^{*})$, then $\tilde{E}$ is a Lyapunov function for the dynamics restricted to $\mathcal{L}$, and the equilibrium $(E^{*},S^{*},z^{*})$ with $E^{*}=E(S^{*},z^{*})$, is stable with respect to the dynamics on $\mathcal{L}$.

Finally, note that the row vector $K^c$ in the new coordinates transforms to $\widetilde{K^c}(\tilde{E},\tilde{S},\tilde{z},{\tilde{p}}_E,{\tilde{p}}_S,{\tilde{p}}_z)$, leading to the transformed power conjugate outputs

$${\tilde{y}}_p:=\frac{\partial {\tilde{K}}^c}{\partial {\tilde{p}}_E},$$

(89)

and the transformed entropy flow conjugate outputs

$${\tilde{y}}_e:=\frac{\partial {\tilde{K}}^c}{\partial {\tilde{p}}_S}$$

(90)

All this is illustrated by the stabilization of the gas-piston system in the following example.

Example 11 

(Regulation of gas-piston system). Consider the gas-piston system (without a damper) with extensive variables $(E,S,V,\pi )$, as before. The constitutive properties of the system are given by the Liouville submanifold as in (65) with energy expression $E_p(S,V,\pi ):=U(S,V)+\frac{{\pi }^2}{2m}$ (‘p’ for plant). Without damping ($d=0)$ the dynamics are generated by the Hamiltonian

$$K=p_V\frac{\pi }{m}-p_{\pi }\frac{\partial U}{\partial V}+\left(p_{\pi }+p_E\frac{\pi }{m}\right)u_p,$$

(91)

with power conjugate output $y_p=\frac{\pi }{m}$ (velocity of the piston). A scalar controller system with extensive variables $(E_c,z_c)$ is given by the port-thermodynamic system $({\mathcal{L}}_c,K_c)$, with energy $E_c=E_c\left(z_c\right)$, and dynamics

$$K_c=\left(p_c+p_{E_c}{E}_c^{\prime }\left(z_c\right)\right)u_c$$

(92)

with output $y_c=E_c^{\prime }\left(z_c\right)$. The function $E_c\left(z_c\right)$ is a design parameter, specifying the controller system.

The closed-loop system is obtained by the negative feedback (with v a new external input)

$$u_p=-y_c+v=-E_c^{\prime }\left(z_c\right)+v,u_c=y_p=\frac{\pi }{m},$$

(93)

together with

$$E:=E_p+E_c,p_{E_p}=p_{E_c}=:p_E$$

(94)

This leads to the closed-loop Hamiltonian

$$K=p_V\frac{\pi }{m}-p_{\pi }\frac{\partial U}{\partial V}+\left(p_{\pi }+p_{E_p}\frac{\pi }{m}\right)\left(-E_c^{\prime }\left(z_c\right)+v\right)+\left(p_{z_c}+p_{E_c}{E}_c^{\prime }\left(z_c\right)\right)\frac{\pi }{m}$$

(95)

It is immediately seen that $C(V,\pi ,z_c)=\Phi (V-z_c)$ for any function $\Phi :\mathbb{R}\to \mathbb{R}$ is a conserved quantity. This motivates a consideration of new canonical coordinates $(\tilde{E},\tilde{V},\tilde{\pi },{\tilde{q}}_c,{\tilde{p}}_E,{\tilde{p}}_V,{\tilde{p}}_{\pi },{\tilde{p}}_c)$, where

$$\tilde{E}=E+\Phi (w=V-z_c),{\tilde{p}}_V=p_V-p_E\frac{\partial \Phi }{\partial w},{\tilde{p}}_{z_c}+p_E\frac{\partial \Phi }{\partial w},$$

(96)

while $\tilde{V}=V,\tilde{\pi }=\pi ,{\tilde{z}}_c=z_c,{\tilde{p}}_E=p_E,{\tilde{p}}_{\pi }=p_{\pi }$. In the new coordinates we compute$\tilde{K}$as

$$\begin{array}{c}\tilde{K}=\left({\tilde{p}}_V+{\tilde{p}}_E\frac{\partial \Phi }{\partial w}\right)\frac{\tilde{\pi }}{m}-{\tilde{p}}_{\pi }\frac{\partial U}{\partial V}-\left({\tilde{p}}_{\pi }+{\tilde{p}}_E\frac{\tilde{\pi }}{m}\right){\overline{E}}_c^{\prime }\left({\tilde{z}}_c\right)+\hfill \\ \left({\tilde{p}}_{z_c}-{\tilde{p}}_E\frac{\partial \Phi }{\partial w}+{\tilde{p}}_E{\overline{E}}_c^{\prime }\left({\tilde{z}}_c\right)\right)\frac{\tilde{\pi }}{m}+\left({\tilde{p}}_{\pi }+{\tilde{p}}_E\frac{\tilde{\pi }}{m}\right)v,\hfill \end{array}$$

(97)

leading to the same power conjugate output ${\tilde{y}}_p=y_p=\frac{\pi }{m}$ (velocity of the piston). For any set-point $V^{*}$ the functions $\Phi$ and $E_c$ should be chosen in such a way that the function $E\left(S,V,\pi ,\Phi (V-z_c)\right)$ has a strict minimum at $(S^{*},V^{*},{\pi }^{*}=0,z_c^{*})$ for some value of $S^{*}$ and state $z_c^{*}$ of the controller system. As discussed before, this can be turned into asymptotic stabilization by additionally interconnecting the obtained closed-loop system with a damper system through the power port $(v,{\tilde{y}}_p)$.

6. Conclusions and Outlook

Ever since the fundamental contributions of Gibbs and Maxwell to thermodynamics, geometry has played an essential role. Nevertheless, the development of the geometric theory of macroscopic thermodynamics still poses fundamental questions, especially when it comes to thermodynamics instead of thermostatics. In this paper, the focus has been on quasi-Hamiltonian formulations based on the factorization of the irreversible entropy production, and on the contact-geometric approach using Liouville geometry. The first topic is intimately related to the GENERIC framework, as well as to the theory of irreversible port-Hamiltonian systems. With respect to the second topic, Liouville geometry offers a versatile framework for dealing with the general thermodynamic phase space, in particular by providing a unification of the energy and entropy representation. Although in this approach the necessary conditions for the dynamics on the thermodynamic phase space are clear, natural specifications of the dynamics are still somewhat lacking. In this regard, a combination of the quasi-Hamiltonian and GENERIC formulations with contact and Liouville geometry should be promising.

Finally, it should not be forgotten that thermodynamics started as an engineering subject (dealing with the efficiency of the steam engine). The interaction of thermodynamic systems with their surroundings is key to the theory. This has been demonstrated in this paper through a discussion of the definitions of the energy and entropy as storage functions with respect to supply rates corresponding to the thermal and the (mechanical) power port, and through the definition of port-thermodynamic systems. Furthermore, it naturally leads to control of thermodynamic systems, including the theory of ‘control by interconnection’.