Modeling and Control Design of Port-Hamiltonian Systems in Discrete-Time

1. Introduction

Port-Hamiltonian systems were introduced about thirty years ago to provide a unified framework for modeling lumped-parameter continuous-time physical systems, [1]. Over the past decades, several control synthesis methodologies have been developed, [2,3,4]. More recently, their discrete-time extension has become a popular framework for the geometric integration of ordinary differential equations (ODEs), paving the way for digital control system design. Time-integration schemes are designed to preserve either the “structure” or the “energy” of continuous-time systems, and the literature distinguishes between two main categories: symplectic and energy-preserving integrators. Discrete-time Hamiltonian and port-Hamiltonian systems based on symplectic integrators have been extensively employed, for example, in [5,6,7,8,9]. This work instead focuses on a class of port-Hamiltonian systems belonging to the family of energy-preserving integrators, [10,11], with a particular emphasis on discrete-time control design.

The discrete-time port-Hamiltonian systems considered in this paper constitute a simple extension of those proposed in [12,13,14,15,16,17,18] and are characterised by dynamics that depend on the discrete gradient of the Hamiltonian, ([10], Definition 3.1). Further details can be found in [19,20], while the extension to distributed-parameter systems is presented in [21,22]. This particular dynamic structure is useful for control synthesis because it directly yields an energy-balance relation for the discrete-time system, thereby simplifying both the design procedure and Lyapunov-based stability analysis. Its main drawback is that the state equation is implicit.

The control design procedure consists of three steps. The first is devoted to developing a discrete-time port-Hamiltonian approximation of the plant that belongs to the class illustrated before. Based on this discrete-time model, the second design step addresses the definition of the control action. One option is to rely on an extension of the energy-shaping plus damping-injection paradigm, [2], for which several contributions have already been presented in the literature, [9,12,19,20]. Here, we propose an extension of the Interconnection and Damping Assignment Passivity-based Control (IDA-PBC) framework, [3]. It is worth mentioning that notable contributions in this research area can also be found in [18,23,24]. The idea is to determine a state-feedback law capable of modifying the internal and dissipative structures of the open-loop system, as well as its Hamiltonian function. In particular, the “new” Hamiltonian is characterised by a (global) minimum at the desired configuration (energy shaping), while the dissipative structure is modified to introduce sufficient dissipation to make this minimum attractive (damping injection). As in continuous time, mapping the open-loop system into the desired closed-loop system requires solving a matching equation. In continuous time, this equation is a nonlinear partial differential equation (PDE), whereas in the digital case it is algebraic in the discrete gradient.

The feedback action acts on the system by modifying its energy function and “structure”. These are two effective degrees of freedom for shaping steady-state and transient behaviour. Since the design and tuning of the regulator are performed in discrete time, the remaining problem is to develop a strategy to couple the digital controller with the continuous-time plant. This constitutes the third and final step of the synthesis procedure. Within energy- or passivity-based frameworks, the problem of interfacing digital and physical systems is usually addressed by preserving the passivity of the interconnection, [25,26]. The idea, instead, is to design an interface that guarantees performance similar to that achieved when the (digital) controller is interconnected with the discrete-time approximation of the plant (i.e., in the reference design scenario), while introducing minimal variations in the nominal control action, [27]. One way to “measure” performance is to consider the variation of the closed-loop Hamiltonian (energy) function along system trajectories. If, at each sampling instant, it decreases at least as much as in discrete time, then no performance degradation occurs. These considerations lead to an optimization problem in which the cost function is the amplitude of the perturbation applied to the nominal control action, while the constraints are given by the continuous-time plant dynamics and an upper bound on the “desired” energy level at the end of the sampling interval. The design procedure is illustrated through the example of a two-degree-of-freedom planar manipulator.

The paper is organized as follows. In Section 2, a summary on discrete-time port-Hamiltonian systems is presented, while an overview of the energy-based control design and the IDA-PBC framework in discrete-time are discussed in Section 3. Section 4 contains the digital-to-analog interfaces arising from the solution of different optimization problems. The example is discussed in Section 5, where the complete design procedure is illustrated in detail. Conclusions are in Section 6.

2. Discrete-Time Port-Hamiltonian Systems

In this paper, we refer to the following class of continuous-time port-Hamiltonian systems, [1]:

$$\left\{\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\left[J\left(x\right(t\left)\right)-R\left(x\right(t\left)\right)\right]\nabla H\left(x\left(t\right)\right)+G\left(x\left(t\right)\right)u\left(t\right)\hfill \\ \hfill y\left(t\right)& =G\left(x\right(t\left)\right)\nabla H\left(x\right(t\left)\right)\hfill \\ \hfill x\left(0\right)& =x_0.\hfill \end{array}\right.$$

(1)

In (1), $t{\in }_{\ge 0}$ and $x\left(t\right)\in X{\subseteq }^n$ are the time and state variable, $u\left(t\right)\in U$ and $y\left(t\right)\in Y$, with $U,Y{\subseteq }^m$, are the input and output, $J\left(x\right),R\left(x\right){\in }^{n\times n}$ are such that $J\left(x\right)=-J\left(x\right)$ and $R\left(x\right)=R\left(x\right)\ge 0$ for all $x\in X$, $G\left(x\right){\in }^{n\times m}$ is full-rank, $H:X\to$ is the Hamiltonian (energy) function, and $x_0\in X$ the initial condition. The matrix-valued functions J and R represent the interconnection and the dissipative structures of the system, respectively. Along system’s trajectories we have that

$$\begin{array}{cc}\hfill \dot{H}\left(x\left(t\right)\right)& =-\nabla H\left(x\right(t\left)\right)R\left(x\right(t\left)\right)\nabla H\left(x\right(t\left)\right)+y\left(t\right)u\left(t\right)\hfill \\ \hfill & \le y\left(t\right)u\left(t\right),\hfill \end{array}$$

(2)

so if H is lower-bounded, (1) is passive and H is the storage function.

The design of a digital control law for (1) relies on the discrete-time approximation proposed in [19,20]. In brief, such a discrete-time model is obtained at first by replacing the time derivative by the finite difference

$$\dot{x}\left(t_k\right)\simeq \frac{1}{\tau }\left(x_{k+1}-x_k\right),$$

(3)

being $t_k=k\tau$ the time samples and $x_k=x\left(t_k\right)$, with $\tau >0$ and $k\in$. Secondly, a discrete approximation of the gradient operator is adopted, see e.g. ([10], Definition 3.1).

Definition 1. 

Let $\mathsf{\Phi }{:}^p\to$ be a continuously differentiable function. A discrete gradient$\nabla \mathsf{\Phi }{:}^p{\times }^p{\to }^p$ is a continuous map such that for all $\phi ,{\phi }_{+}{\in }^p$ we have

$$\begin{array}{cc}\hfill \left({\phi }_{+}-\phi \right)\nabla \mathsf{\Phi }(\phi ,{\phi }_{+})& =\mathsf{\Phi }\left({\phi }_{+}\right)-\mathsf{\Phi }\left(\phi \right),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \underset{{\phi }_{+}\to \phi }{lim}\nabla \mathsf{\Phi }(\phi ,{\phi }_{+})& =\nabla \mathsf{\Phi }\left(\phi \right).\hfill \end{array}$$

(5)

Typical examples are the mean value discrete gradient [28] defined as

$$\nabla \mathsf{\Phi }(\phi ,{\phi }_{+})={\int }_0^1\nabla \mathsf{\Phi }\left((1-\sigma )\phi +\sigma {\phi }_{+}\right)\sigma ,$$

(6)

or the Gonzalez discrete gradient, [10]. For a quadratic function $\mathsf{\Phi }\left(\phi \right)=\phi Q\phi$ with $Q=Q$, it is easy to check that $\nabla \mathsf{\Phi }(\phi ,{\phi }_{+})=Q(\phi +{\phi }_{+})$.

If we combine (1), (3) and Definition 1, for any $k\in$, the discrete-time formulation of the port-Hamiltonian dynamics for control design is, [19,20]:

$$\left\{\begin{array}{cc}\hfill x_{k+1}& =x_k+\tau \left[\overline{J}(x_k,x_{k+1})-\overline{R}(x_k,x_{k+1})\right]H(x_k,x_{k+1})+\tau \overline{G}(x_k,x_{k+1})u_k\hfill \\ \hfill y_k& =\overline{G}(x_k,x_{k+1})\nabla H(x_k,x_{k+1}),\hfill \end{array}\right.$$

(7)

with $u_k=u\left(t_k\right)$ and initial condition $x_0$. This class of discrete-time systems can be seen as an extension of [12,15,18]. The functions $\overline{J},\overline{R}:X\times X{\to }^{n\times n}$ and $\overline{G}:X\times X{\to }^{n\times m}$ are discrete approximations of J, R and G. More precisely, for all $x,x_{+}\in X$, we require that $\overline{J}(x,x_{+})=-\overline{J}(x,x_{+})$, $\overline{J}(x,x)=J\left(x\right)$, $\overline{R}(x,x_{+})=\overline{R}(x,x_{+})\ge 0$, $\overline{R}(x,x)=R\left(x\right)$, and $\overline{G}(x,x)=G\left(x\right)$. Admissible choices for $\overline{J}$ are $J\left(\frac{1}{2}\left(x+x_{+}\right)\right)$, $\frac{1}{2}J\left(x\right)+\frac{1}{2}J\left(x_{+}\right)$ or simply $J\left(x\right)$, and similarly for $\overline{R}$ and $\overline{G}$. From (4) and (7), we get the discrete-time version of (2):

$$\begin{array}{cc}\hfill H\left(x_{k+1}\right)-H\left(x_k\right)& ={\nabla }^{T}H(x_k,x_{k+1})(x_{k+1}-x_k)\hfill \\ \hfill & =-\tau {\nabla }^{T}H(x_k,x_{k+1})\overline{R}(x_k,x_{k+1})\nabla H(x_k,x_{k+1})+\tau y_k{u}_k\hfill \\ \hfill & \le \tau y_k{u}_k.\hfill \end{array}$$

(8)

Example 1. 

In the linear case, the port-Hamiltonian system (1) takes the following form:

$$\left\{\begin{array}{cc}\hfill \dot{x}\left(t\right)& =\left(J-R\right)Qx\left(t\right)+Gu\left(t\right)\hfill \\ \hfill y\left(t\right)& =GQx\left(t\right)\hfill \\ \hfill x\left(0\right)& =x_0,\hfill \end{array}\right.$$

(9)

where $J,R{\in }^{n\times n}$ are constant and such that $J=-J$ and $R=R\ge 0$, while $G{\in }^{n\times m}$ is constant and full rank. In (9), the Hamiltonian function is quadratic, more precisely $H\left(x\right)=\frac{1}{2}xQx$, with $Q=Q>0$, which implies that $\nabla H\left(x\right)=Qx$. It is immediate to see that the discrete-time version of (9) is given by

$$\left\{\begin{array}{cc}\hfill x_{k+1}& =x_k+\frac{\tau }{2}\left(J-R\right)Q(x_k+x_{k+1})+\tau G{u}_k\hfill \\ \hfill y_k& =\frac{1}{2}GQ(x_k+x_{k+1}).\hfill \end{array}\right.$$

(10)

Now, since the matrix $I-\frac{\tau }{2}(J-R)Q$ is invertible because the eigenvalues of the matrix $(J-R)Q$ cannot have positive real part, the implicit state equation can be made explicit. Consequently, the state equation in (10) can be rewritten as

$$x_{k+1}=\left[I-\frac{\tau }{2}(J-R)Q\right]\left[I+\frac{\tau }{2}(J-R)Q\right]x_k+\tau \left[I-\frac{\tau }{2}(J-R)Q\right]G{u}_k$$

since $\nabla H(x,x_{+})=\frac{1}{2}Q(x+x_{+})$, thus the explicit dynamics is described by a linear system with feedthrough.

Example 2. 

Let us consider a mechanical system locally described by n configuration variables $q=(q_1,\dots ,q_n){\in }^n$. The kinetic energy is $\frac{1}{2}\dot{q}M\left(q\right)\dot{q}$, being $M\left(q\right)=M\left(q\right)>0$ for all q the generalized mass matrix, while we let $V\left(q\right)$ to denote the potential energy. As illustrated in [1], if no constraints are present, the dynamics can be described by a port-Hamiltonian system in the form (1), more precisely:

$$\left\{\begin{array}{cc}\hfill \left(\begin{array}{c}\dot{q}\left(t\right)\\ \dot{p}\left(t\right)\end{array}\right)& =\left(\begin{array}{cc}0& I\\ -I& -D\end{array}\right)\left(\begin{array}{c}\frac{\partial H}{\partial q}\left(q\right(t),p(t\left)\right)\\ \frac{\alpha H}{\alpha p}\left(q\right(t),p(t\left)\right)\end{array}\right)+\left(\begin{array}{c}0\\ B\left(q\right)\end{array}\right)u\left(t\right)\hfill \\ \hfill y\left(t\right)& =B^T\left(q\right)\frac{\partial H}{\partial p}\left(q\right(t),p(t\left)\right),\hfill \end{array}\right.$$

(11)

where $p=M\left(q\right)\dot{q}$ are the generalized momenta, the Hamiltonian function is the total energy

$$H(q,p)={\displaystyle \frac{1}{2}}{p}^{T}M\left(q\right)p+V\left(q\right),$$

(12)

$u,y{\in }^m$ are the generalized forces and velocities, respectively. Besides, in (11), we assume that $B\left(q\right)$ is full rank, while $D=D^T\ge 0$ takes into account the viscous friction forces, if present. It is easy to see that the discrete-time approximation of (11) is

$$\left\{\begin{array}{cc}\hfill q_{k+1}& =q_k+\tau {\nabla }_{p}H(q_k,p_k,q_{k+1},p_{k+1})\hfill \\ \hfill p_{k+1}& =p_k-\tau \left[{\nabla }_{q}H(q_k,p_k,q_{k+1},p_{k+1})-D{\nabla }_{p}H(q_k,p_k,q_{k+1},p_{k+1})\right]+\hfill \\ & \tau \overline{B}(q_k,q_{k+1})u_k\hfill \\ \hfill y_k& =B^T(q_k,q_{k+1}){\nabla }_{p}H(q_k,p_k,q_{k+1},p_{k+1}),\hfill \end{array}\right.$$

(13)

where B is a discrete approximation of B. For simplicity, we can set $B(q_k,q_{k+1})=B\left(q_k\right)$. Since the Hamiltonian function in (12) is non-quadratic, it is not always possible to employ the mean value discrete gradient, as obtaining an explicit expression for the integral in (6) may be difficult or even impossible. In such cases, one may either resort to an approximation of the integral or adopt a different discrete gradient, such as the Gonzalez discrete gradient, which is constructed by sampling the continuous one. Finally, note that the state equation in (13) is given in implicit form. Consequently, unlike the linear case discussed in Example 1, an explicit state-update equation cannot, in general, be derived.

Remark 1. 

System (7) is not based on a particular discrete gradient, which turns out to be a degree-of-freedom at disposal to the designer, together with the expression of $J$, $R$ and $G$. Such a model for control design is a compromise between accuracy in the response and (numerical) complexity of the final control law. The drawback is that the state equation is implicit: this is the price to pay to have “for free” the energy-balance relation (8), which is the starting point for control design and stability analysis.

Remark 2. 

As mentioned in Remark 1 and reported in Examples 1 and 2, the state equation in (7) is in implicit form. It can be made explicit, i.e. solved for $x_{k+1}$ in the linear case discussed in Example 1 because H is quadratic and J, R and G in (1) are constants. This approach has been followed in [12,18]. In the nonlinear case, however, it is necessary to keep such an implicit formulation. Based on passivity arguments and inspired by [29], under mild conditions on the structure of the system and on the discrete gradient, in [20] has been shown that trajectories exist independently from the sampling time. This is the same rationale pursued in [15] where a similar result has been proved under the hypothesis that the sampling time is sufficiently small. In this paper we assume that the dynamical equations are well-posed, so the next state can be always (at least numerically) computed.

3. Energy-Based Control Design in Discrete-Time

This section aims at generalising the Interconnection and Damping Assignment Passivity-Based Control (IDA-PBC) design methodology [3], originally developed for continuous-time port-Hamiltonian systems, to the discrete-time setting. Remarkable contributions in this direction include [18,23] and [24], the latter being specifically focused on mechanical systems. These methodologies aim at assigning an equilibrium $x_{★}{\in }^n$ to (7) by means of state feedback, while preserving the port-Hamiltonian structure. If the closed-loop system retains the port-Hamiltonian form, asymptotic convergence to $x_{★}$ can be established by relying on a balance relation similar to (8).

We start from the balance relation (8) for the open-loop dynamics (7), which can be compactly rewritten as

$${\tau }^{-1}\left[H\left(x_{k+1}\right)-H\left(x_k\right)\right]=y_k^T{u}_k-d(x_k,x_{k+1}),$$

where $d(x_k,x_{k+1})\ge 0$ takes into account the dissipated power during the k-th sampling interval. Similarly to the continuous-time case [2,3], the idea behind the design of an energy-based state-feedback controller is to find

$$u_k=\beta (x_k,x_{k+1})+v_k$$

(14)

so that the closed-loop dynamics obey the “new” balance relation

$${\tau }^{-1}\left[H_d\left(x_{k+1}\right)-H_d\left(x_k\right)\right]=z_k^T{v}_k-d_d(x_k,x_{k+1}),$$

(15)

in which $H_d\left(x_k\right)$ is a “desired” energy with a (local) minimum at $x_{★}$, $d_d(x_k,x_{k+1})\ge 0$ replaces the natural dissipation, increased to improve the convergence rate, and $z_k$ is the new passive output.

For simplicity, we assume that the matrix-valued functions J, R and G in (1) are constant, and consequently that their discrete-time counterparts $J$, $R$ and $\overline{G}$ in (7) are also constant. The general case, in which these matrices depend on the state, is more involved. Nevertheless, most of the results and considerations developed for the simplified setting can be extended to the state-dependent case. The most natural extension of the IDA-PBC approach to the discrete-time setting requires to determine the feedback law $\beta (x_k,x_{k+1})$ in (14) so that the state equation in (7) is mapped to

$$x_{k+1}=x_k+\tau \left(J_d-R_d\right)\nabla H_d(x_k,x_{k+1})+\tau G{v}_k,$$

(16)

which represents the “target” or “desired” dynamics. In (16), we have that

$$\begin{array}{cc}\hfill H_d\left(x_k\right)& =H\left(x_k\right)+H_a\left(x_k\right)\hfill \\ \hfill J_d& =J+J_a,\mathrm{with}{J}_a=-J_a^T\hfill \\ \hfill R_d& =R+R_a,\mathrm{with}{R}_a=R_a^T\ge 0,\hfill \end{array}$$

(17)

are the desired Hamiltonian function, interconnection matrix and dissipative structure, respectively. It is immediate to see that, if we let

$$z_k=G^T\nabla H_d(x_k,x_{k+1})=y_k+G^T\nabla H_a(x_k,x_{k+1})$$

(18)

the target system (16) equipped with the output $z_k$ satisfies the balance relation (15).

The comparison between the open-loop dynamics (7) and the target one (16) shows that the control action $\beta (x_k,x_{k+1})$ modifies the energy and the interconnection and dissipative structures of the original system. As far as the so-called energy-shaping action is concerned, the function $H_a\left(x_k\right)$ selected so that the closed-loop energy $H_d\left(x_k\right)$ in (17) has a minimum in $x_{★}$. Besides, $R_a$ is chosen to improve the convergence rate towards the minimum of the closed-loop energy function. A similar effect is achievable via output feedback (damping injection) by imposing that $v_k=-\phi \left(z_k\right)$ in (14), with the function $\phi$ such that $z^T\phi \left(z\right)\ge 0$ for all $z{\in }^m$, and $z_k$ defined in (18). The simplest damping injection law is

$$v_k=-K_z{z}_k=-K_z{G}^T\nabla H_d(x_k,x_{k+1})=-K_z\left[y_k+G^T\nabla H_a(x_k,x_{k+1})\right],$$

(19)

with $K_z{\in }^{m\times m}$ and $K_z=K_z^T\ge 0$. Note that such a control action can be also applied to the open-loop system (7) by setting $u_k=-\phi \left(y_k\right)$. If $x_{★}$ is the only invariant solution of the steady state, the “new” dissipative structure in the system forces the asymptotic stability of the equilibrium. With simple calculations, we get that $H_a\left(x_k\right)$, $J_a$, $R_a$ and $\beta (x_k,x_{k+1})$ have to satisfy the matching equation:

$$G\beta (x_k,x_{k+1})=(J_a-R_a)\overline{\nabla }H(x_k,x_{k+1})+(J_d-R_d)\overline{\nabla }{H}_a(x_k,x_{k+1})$$

(20)

In the next proposition, sufficient conditions for the solution of (20) so that the state-feedback action (14) asymptotically stabilises a (forced) equilibrium $x_{★}{\in }^n$ of (7) are presented. This result is an immediate extension of [20].

Proposition 1. 

Denote by $x_{★}{\in }^n$ a (forced) equilibrium for (7), and let ${\overline{K}}_a{:}^n{\times }^n{\to }^n$ be a function such that

$$G^{\perp }\left[(J_a-R_a){\nabla }^{¯}H(x_k,x_{k+1})+(J_d-R_d){\overline{K}}_a(x_k,x_{k+1})\right]=0$$

(21)

being G^⊥ the full-rank left annihilator of G, i.e. $G^{\perp }G=0$. If

$$\frac{\partial K_a}{\partial x}\left(x\right)={\displaystyle \frac{{\partial }^T{K}_a}{\partial x}}\left(x\right),$$

with $K_a\left(x\right){\overline{K}}_a(x,x)$, then there exists $H_a{:}^n\to$ for which $K_a\left(x\right)=\frac{\partial H_a}{\partial q}\left(x\right)$. Besides, if for all $x,\xi {\in }^n$, we have that

$$(\xi -x){\overline{K}}_a(x,\xi )=H_a\left(\xi \right)-H_a\left(x\right),$$

then ${\overline{K}}_a$ is the discrete gradient of $H_a$, i.e. $\overline{\nabla }{H}_a(x,\xi )={\overline{K}}_a(x,\xi )$, and the control action (14) with

$$\begin{array}{cc}\hfill \beta (x_k,x_{k+1})& =G^{+}\left[(J_a-R_a)\nabla H(x_k,x_{k+1})+(J_d-R_d){\overline{K}}_a(x_k,x_{k+1})\right]\hfill \\ \hfill & =G^{+}\left[(J_a-R_a)\overline{\nabla }H(x_k,x_{k+1})+(J_d-R_d)\overline{\nabla }{H}_a(x_k,x_{k+1})\right],\hfill \end{array}$$

(22)

being G⁺ the pseudo-inverse of G, maps (7) to the discrete-time system (16) with Hamiltonian function, interconnection and dissipative structures defined in (17), and output (18). Finally, if

$$\begin{array}{cc}\hfill \nabla H(x_{★},x_{★})+{\overline{K}}_a(x_{★},x_{★})& =0\hfill \\ \hfill (x-x_{★})\left[\nabla H(x_{★},x)+{\overline{K}}_a(x_{★},x)\right]& \ge 0\hfill \end{array}$$

(23)

for all $x\in B_{\delta }\left(x_{★}\right)=\{\xi {\in }^n:x-\xi \le \delta \}$ and for some $\delta >0$, then $x_{★}$ is a locally stable equilibrium for the closed-loop system with $v_k=0$. Furthermore, if the largest invariant set contained in

$$\left\{x{\in }^n\mid \exists \xi {\in }^n\mathrm{s}.\mathrm{t}.{\nabla }^T{H}_d(x,\xi )R_d\nabla H_d(x,\xi )=0\right\}\cap B_{\delta }\left(x_{★}\right)$$

(24)

is $\left\{x_{★}\right\}$, then $x_{★}$ is a locally asymptotically stable equilibrium for the closed-loop system with $v_k=0$.

Proof. 

The result can be proved in the same way as ([20], Proposition 3.1). Note that asymptotic stability is a consequence of the La Salle’s invariance principle. In fact, for the closed-loop system we have $H_d\left(x_{k+1}\right)-H_d\left(x_k\right)=-{\nabla }^T{H}_d(x_k,x_{k+1})R_d\nabla H_d(x_k,x_{k+1})\le 0$, thus the corresponding energy function decrease until the trajectories evolve inside the set (24). Since the largest invariant set contained in (24) equals $\left\{x_{★}\right\}$, the result immediately follows. □

A crucial point in the methodology of Proposition 1 is the solvability of the matching equation (21), which can be re-written in terms of the discrete gradient of $H_a$ as

$$G^{\perp }\left[(J_a-R_a)\nabla H(x_k,x_{k+1})+(J_d-R_d)\nabla H_a(x_k,x_{k+1})\right]=0.$$

(25)

The simplest approach is to rely on the PDE that follows from (25) once $x_{k+1}\to x_k\equiv x$, i.e.:

$$G^{\perp }\left[(J_a-R_a)\frac{\partial H}{\partial x}\left(x\right)+(J_d-R_d)\frac{\partial H_a}{\partial x}\left(x\right)\right]=0.$$

(26)

When the matrix-valued functions J, R, and G in (1) depend on the state variable x, it is not straightforward to derive a solution to (25) from its continuous-time counterpart (26). However, in the constant case considered in this paper, the function $H_a$ and the matrices $J_a$ and $R_a$ that satisfy (26) also satisfy (25). This result is summarized in the following corollary.

Corollary 1. 

Let us assume that there exists $H_a:X\to$, $J_a,R_a{\in }^{n\times n}$ and $G{\in }^{m\times n}$ such that the continuous-time matching equation (26) holds for the port-Hamiltonian system (1). Then, $H_a$, $J_a$ and $R_a$ are solution also of the discrete-time reformulation (25).

Proof. 

If we consider, for example, the mean value discrete gradient (6), since (26) holds for all $x\in X$, the result follows by replacing x with $(1-\sigma )x_k+\sigma x_{k+1}$, where $\sigma \in$, and then integrating with respect to $\sigma$. □

Remark 3. 

Since the control action (14) with $\beta (x_k,x_{k+1})$ defined in (22) and $v_k$ in (19) depends on $x_{k+1}$, we have to investigate how to compute $u_k$ at step k. The first important assumption is that the closed-loop dynamics obtained from (7), (14) and (19) is well-posed, see ([20], Proposition 2.1). The corresponding state equation is

$$x_{k+1}=x_k+\tau \left[J_d-\left(R_d+G{K}_z{G}^T\right)\right]\nabla H_d(x_k,x_{k+1}).$$

(27)

For (27) at each step k, we let

$${\mathcal{F}}_{d,k}\left(\xi \right)x_k+\tau \left[J_d-\left(R_d+G{K}_z{G}^T\right)\right]\nabla H_d(x_k,\xi ),$$

where $\xi {\in }^n$. So, ξ is the solution of $\xi ={\mathcal{F}}_{d,k}\left(\xi \right)$, with $\xi \equiv x_{k+1}$ only in the nominal case, i.e. when the plant exactly behaves as (27), while (14) and (19) provide $u_k$. This approach is conceptually similar to [9], where the control design problem is tackled for discrete-time models of mechanical systems in port-Hamiltonian form arising from symplectic integration schemes (i.e., the implicit midpoint rule). On the other hand, it differs from [12,18] where the implicit state equation is approximated with an explicit one.

Remark 4. 

The key step in the procedure of Proposition 1 is to solve (21) for suitable matrices ${\overline{K}}_a$, $J_a$, and $R_a$, where ${\overline{K}}_a$ must be the discrete gradient of a function $H_a$. If the structure matrices of the port-Hamiltonian system (1) are not constant, Corollary 1 cannot be applied, and solving the matching equation may become considerably more involved. A different approach has been discussed in [20]. The idea is to exploit the algebraic nature of (21) to derive a broader class of control laws of the form (22), capable of performing an (approximate) energy-shaping action on the open-loop system (7) without requiring ${\overline{K}}_a$ to be the discrete gradient of an energy function. This approach is the discrete-time extension of the technique proposed in [30]. This topic is beyond the scope of the present paper.

4. “Performance-Preserving” Coupling Strategies

The aim of this section is to present a strategy for coupling the discrete-time regulator designed according to the procedure illustrated in Section 3 with the continuous-time plant (1). Since the controller is designed on the basis of the discrete-time model (7), its performance is typically assessed through numerical simulations using the same model as a reference. The proposed coupling strategy therefore aims at applying the “smallest” perturbation to the nominal stabilising control law that guarantees performance as close as possible to that of the reference case. In particular, the transient and steady-state behaviours are preserved if the variation of the closed-loop Hamiltonian function $H_d$ defined in (17) along the system trajectories remains unchanged. These considerations lead to an optimization problem that must be solved at each step $k\in$, resulting in a novel digital-to-analog interface.

For any energy-based control synthesis methodology such as the ones that fit into the framework illustrated in Section 3, let us consider a coupling algorithm based on the solution of the following nonlinear optimization problem at each step $k\in \mathbb{N}$,27]:

$$\begin{array}{cc}\hfill & \underset{u_{\delta },x\left(t\right)}{min}{u_{\delta }}^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \dot{x}\left(t\right)=\left[J\left(x\right(t\left)\right)-R\left(x\right(t\left)\right)\right]\nabla H\left(x\left(t\right)\right)+G\left(x\left(t\right)\right)\left[\beta (x_k,x_{k+1})+u_{\delta }\right]\hfill \\ \hfill & x\left(t_k\right)=x_k\hfill \\ \hfill & \beta (x_k,x_{k+1})+u_{\delta }\in U\hfill \\ \hfill & {u_{\delta }}^2\le \Delta \hfill \\ \hfill & H_d\left(x\left(t_{k+1}\right)\right)\le H_d\left(x_{k+1}\right).\hfill \end{array}$$

(28)

The idea is simple. At each $t_k$, the state variable of (1) is acquired, and the control action (14) is computed by relying on the discrete-time model (7) in which $x_k=x\left(t_k\right)$. For simplicity and without loss of generality, in (14) we assume that $v_k=0$, i.e. the required dissipation is introduced via the state feedback law $\beta$. As discussed in Remark 3, the design procedure provides the “next” state $x_{k+1}$ and $u_k=\beta (x_k,x_{k+1})$ so, in the ideal case, the value of the closed-loop energy in $t_{k+1}$ is known and equal to $H_d\left(x_{k+1}\right)$. The digital-to-analog interface requires to compute the “smallest” perturbation $u_{\delta }\in {\mathbb{R}}^m$ to the nominal control input $\beta (x_k,x_{k+1})$ so that the control action $u\left(t\right)=\beta (x_k,x_{k+1})+u_{\delta }$ applied to (1) during $[t_k,t_{k+1})$ is admissible, i.e. it belongs to the set U, and the value of $H_d$ in $x\left(t_{k+1}\right)$ is not greater than $H_d\left(x_{k+1}\right)$. The result is that the closed-loop energy function decreases at least as in the reference discrete-time case. The parameter $\Delta >0$ is an upper bound for $u_{\delta }$.

Remark 5. 

The feasibility of (28) is not always guaranteed. This is the major issue of the approach since at every $t_k$ the perturbation $u_{\delta }$ has to be computed to obtain the control action $u_k$. As in ([31] Theorem 2.5]), thanks to the Gronwall-Bellman inequality, it is possible to show that there could be sampling instants for which $H_d$ does not decrease as expected, i.e. that $x\left(t_{k+1}\right)$ does not belong to the “desired” level set, as specified by the last inequality constraint in (28). A similar problem has been tackled in [25] by dissipating the difference between achieved and desired energy at $t_{k+1}$ in the next sampling interval(s). This strategy is applicable here, but this topic is not addressed in this paper. Soft constraints are a simple way to tackle the feasibility issue of (28) in a simple way. The corresponding optimization problem is

$$\begin{array}{cc}\hfill & \underset{u_{\delta },\epsilon ,x\left(t\right)}{min}{u_{\delta }}^2+K_{\epsilon }{\epsilon }^2\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& \dot{x}\left(t\right)=\left[J\left(x\right(t\left)\right)-R\left(x\right(t\left)\right)\right]\nabla H\left(x\left(t\right)\right)+G\left(x\left(t\right)\right)\left[\beta (x_k,x_{k+1})+u_{\delta }\right]\hfill \\ \hfill & x\left(t_k\right)=x_k\hfill \\ \hfill & H_d\left(x\left(t_{k+1}\right)\right)\le H_d\left(x_{k+1}\right)+{\epsilon }^2.\hfill \end{array}$$

(29)

where we have removed the amplitude constraint on $u_{\delta }$ and assumed that $U\equiv {\mathbb{R}}^m$. In (29), we let $K_{\epsilon }\gg 1$ and $U\equiv {\mathbb{R}}^m$, so the variable $\epsilon \in \mathbb{R}$ allows to dynamically relax the constraint on the closed-loop energy function $H_d$ at $t_{k+1}$. The idea is to compute the “smallest” correction $u_{\delta }$ to the nominal control action $\beta (x_k,x_{k+1})$ so that the $H_d\left(x\left(t_{k+1}\right)\right)$ is as close as possible to the target value $H_d\left(x_{k+1}\right)$. The mismatch is reduced by selecting $K_{\epsilon }$ sufficiently large. It is easy to see that (29) is always feasible, and that existence of solution follows from the cost function being lower-bounded.

Noticeable limitations of (28) and (29) are that the constraint on $x\left(t\right)$ depends on a nonlinear ODE and that the closed-loop Hamiltonian function $H_d\left(x\right)$ is, in general, non-quadratic. Consequently, the real-time implementation of the coupling mechanism is potentially computationally demanding. Under a few conditions, an explicit solution is possible when the plant is linear and the closed-loop Hamiltonian is quadratic. The result is based on ([32], Proposition 1), reported below.

Proposition 2. 

Let $a\in {\mathbb{R}}^n$, $P\in {\mathbb{R}}^{n\times n}$ and such that $P=P^T>0$, and $r\in \mathbb{R}$, with $r>0$. The solution $v_{★}$ of the convex optimization problem

$$\begin{array}{cc}\hfill & \underset{v}{min}{(v-a)}^{T}P(v-a)\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {(v-c)}^{T}P(v-c)\le r^2\hfill \end{array}$$

(30)

is

$$v_{★}=\left\{\begin{array}{cc}a\hfill & \mathrm{if}{(a-c)}^{T}P(a-c)\le r^2\hfill \\ {\displaystyle \frac{r(a-c)}{\sqrt{(a-c)P(a-c)}}}+c\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(31)

Let us consider the linear port-Hamiltonian system (9), which can be rewritten in the $(A,B,C)$ form with $A=(J-R)Q$, $B=G$, and $C=G^{T}Q$. Let us also assume that the discrete-time control design procedure assigns target dynamics in port-Hamiltonian form with quadratic Hamiltonian

$$H_d\left(x\right)={\displaystyle \frac{1}{2}}{(x-x_{★})}^T{Q}_d(x-x_{★}),Q_d=Q_d^T>0.$$

(32)

With simple computations, it is easy to see that, at the end of the k-th sampling interval, starting from $x\left(t_k\right)=x_k$ and under the control input $u_k=\beta (x_k,x_{k+1})+u_{\delta }$, (9) reaches

$$x\left(t_{k+1}\right)={\mathcal{A}}_k+{\mathcal{B}}_{\delta }{u}_{\delta },$$

(33)

where

$$\begin{array}{cc}\hfill {\mathcal{A}}_k& =e^{A\tau }\left[x_k+\left({\int }_0^{\tau }{e}^{-A\sigma }\sigma \right)B\beta (x_k,x_{k+1})\right]\hfill \\ \hfill {\mathcal{B}}_{\delta }& =e^{A\tau }\left({\int }_0^{\tau }{e}^{-A\sigma }\sigma \right)B.\hfill \end{array}$$

Consequently, the energy constraint $H_d\left(x\left(t_{k+1}\right)\right)\le H_d\left(x_{k+1}\right)$ in (28) takes the form of the quadratic constraint in (30), with $P={\mathcal{B}}_{\delta }^T{Q}_d\delta$, $c={\mathcal{B}}_{\delta }^{+}\left(x_{★}{-}_k\right)$, and $r^2=2H_d\left(x_{k+1}\right)$. Thus, the coupling strategy is recast as the optimization problem of Proposition 2, in which $a=0$. Consequently,

$$u_{\delta }=\left\{\begin{array}{cc}0\hfill & \mathrm{if}\frac{1}{2}{\Vert x_{★}-{\mathcal{A}}_k\Vert }_{Q_d}^2\le H_d\left(x_{k+1}\right)\hfill \\ \left[1-{\displaystyle \frac{\sqrt{2H_d\left(x_{k+1}\right)}}{{\Vert x_{★}-{\mathcal{A}}_k\Vert }_{Q_d}}}\right]\left(x_{★}{-}_k\right)\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

where ${\Vert x_{★}-{\mathcal{A}}_k\Vert }_{Q_d}^2={(x_{★}-{\mathcal{A}}_k)}^T{Q}_d\left(x_{★}{-}_k\right)$. Note that this result can be used to address the general nonlinear case approximately, by relying, at each sampling interval $k\in \mathbb{N}$, on a linear approximation of the plant dynamics around $x_k$.

5. Numerical Example: The 2-Dof Planar Robotic Manipulator

Let us consider the two degrees-of-freedom planar manipulator of Figure 1, which belongs to the class illustrated in Example 2.

The generalised coordinates are $q=(q_1,q_2)\in {\mathbb{R}}^2$, and the total energy is (12) in which the inertia matrix is

$$M\left(q\right)=\left(\begin{array}{cc}{j}_1+j_2+m_1{a}_1^2+m_2\left(l_2^2+a_2^2+2a_2{l}_{1}cos{q}_2\right)& j_2+m_2\left(a_2^2+a_2{l}_{1}cos{q}_2\right)\\ j_2+m_2\left(a_2^2+a_2{l}_{1}cos{q}_2\right)& j_2+m_2{a}_2^2\end{array}\right)$$

and the potential energy is

$$V\left(q\right)=m_{1}g{a}_{1}sin{q}_1+m_{2}g\left[l_{1}sin{q}_1+a_{2}sin(q_1+q_2)\right],$$

(34)

with g the gravity acceleration. The continuous-time model is (11), with $B\left(q\right)=I$ and $D=dI$, where $d\ge 0$. The inputs $u\left(t\right)\in {\mathbb{R}}^2$ are the torques at each joint. The parameters are reported in Table 1.

The discrete-time approximation of the plant model is (13), and in the numerical examples, the mean value discrete gradient (6) is employed. Note that from (34), for any $q,q_{+}{\in }^2$, we have that

$$V(q,q_{+})=\left(\begin{array}{c}{\overline{\nabla }}_{q_1}V(q,q_{+})\\ {\overline{\nabla }}_{q_2}V(q,q_{+})\end{array}\right),$$

(35)

where $q_{+}=(q_{1+},q_{2+})$ and

$$\begin{array}{cc}\hfill {\overline{\nabla }}_{q_1}V(q,q_{+})& =\left(m_{1}g{a}_1+m_{2}g{l}_1\right){\displaystyle \frac{sin{q}_{1+}-sin{q}_1}{q_{1+}-q_1}}+m_{2}g{a}_2{\displaystyle \frac{sin(q_{1+}+q_{2+})-sin(q_1+q_2)}{q_{1+}+q_{2+}-q_1-q_2}}\hfill \\ \hfill {\overline{\nabla }}_{q_2}V(q,q_{+})& =m_{2}g{a}_2{\displaystyle \frac{sin(q_{1+}+q_{2+})-sin(q_1+q_2)}{q_{1+}+q_{2+}-q_1-q_2}}.\hfill \end{array}$$

Based on the general framework summarized in Proposition 1, it is easy to check that the matching equation (21) with $J_a=R_a=0$ is satisfied if $H_a$ depends on the displacements $q_k$ only. As a consequence, the control action $\beta$ only shapes the potential energy. The simplest choice is to have $H_a\left(q\right)=-V\left(q\right)+\frac{1}{2}(q-q_{★})K_p(q-q_{★})$, with $K_p=K_p>0$, so that the desired Hamiltonian function $H_d(q,p)=H(q,p)+H_a\left(q\right)$ has a global minimum in $x_{★}$, which means that the conditions (23) are satisfied.. It is interesting to note that an explicit expression for $\beta$ can be computed based on (35) and on the fact that the further term in $H_a$ is quadratic. This leads to the PD plus gravity compensation controller. Asymptotic convergence to the global minimum $x_{★}$ of $H_d(q,p)$ is obtained via the damping injection law (19), where $z_k=y_k={\overline{\nabla }}_{p}H(x_k,x_{k+1})$ is the output of either the plant and of the closed-loop system. As far as the controller parameters are concerned, we have that $K_p=60I$ and $K_z=20I$.

The block diagram of the complete system resulting from the coupling of the discrete-time controller and the continuous-time plant is reported in Figure 2 to show how the control action is implemented in practice.

Since the stabilising law is based on state-feedback, the starting point is to sample the state $x\left(t\right)$ of the plant in $t=t_k$. The sampled variable $x_k$ is used to determine the “desired” next state ${\xi }_k$ according to the target dynamics. As discussed in Remark 3, the discrete variable ${\xi }_k$ has been introduced to distinguish from $x_{k+1}$ that denotes the sampled value of $x\left(t\right)$ in $t=t_{k+1}$. Getting ${\xi }_k$ requires the solution of an implicit equation. This step is represented by the grey box. Based on $x_k$ and ${\xi }_k$, the control action $\beta$ is computed, while the “corrective” term $u_{\delta }$ follows from the strategies discussed in Section 4. In both the cases, we let $x_{k+1}\equiv {\xi }_k$.

At first, the behaviour obtained in the discrete-time case, i.e. when the control law is applied to the discrete-time formulation (13) of the plant, is reported in Figure 3.

As a matter of fact, this represent the desired response that the designer aims at achieving when the digital controller is coupled with the “real” plant or, in this case, to its continuous-time model (11). The synthesis of the control law and the tuning of the regulator parameters is carried out in the discrete-time setting, and the idea is to get analogous performances in the real-world scenario. Starting from the top graph of Figure 3, the evolution of the discrete displacements $q_k$, of the control input $u_k$ and of the closed-loop Hamiltonian function $H_d\left(x_k\right)$ are reported. This latter response summarizes the transient behaviour of the controlled system and it is the main information required by the coupling techniques discussed in Section 4.

We analyse the coupling strategy based on the optimization problem (29) in which $K_{\epsilon }={10}^4$. This sub-optimal formulation of (28) has been preferred because it does not suffer from any feasibility issue. It is worth mentioning, however, that for the specific case discussed here (28) is always feasible. The closed-loop system behaviour is reported in Figure 4 and Figure 5. A comparison between Figure 3 and Figure 4 shows that the performances obtained when the discrete-time controller is coupled with the continuous-time plant are practically the same as in the reference case. The role of the interconnection mechanism is summarized in Figure 5. The correction term $u_{\delta }$ is different from 0 only in the initial part of the transitory and its amplitude is noticeably lower than the nominal control action $\beta$. The design procedure carried out in discrete-time by relying on the class of implicit port-Hamiltonian systems (7) or (13), in this specific case, ensures very good performances by its own. Consequently, the constraint on the amplitude of $\Vert u_{\delta }\Vert$, i.e. the parameter $\Delta$ that appears in the optimization problem (28), is not strictly necessary from a “practical” point of view in a real-world scenario. In some sense, this motivates the choice carried out to define the optimization problem (29) with the goal of not having any feasibility issues.

In the graph at the bottom of Figure 5, the evolution of $\Delta H_d\left(t_k\right)=H_d\left({\xi }_{k-1}\right)-H_d\left(x\left(t_k\right)\right)$ is reported. Such a quantity represents the difference between target and achieved closed-loop Hamiltonian at each sample instant. In other words, at each sampling interval k, the ${\xi }_k$ variable is computed. It represents the next state according to the desired evolution of the controlled system in discrete-time. This means that $H_d\left({\xi }_k\right)$ is the target value for the closed-loop energy to be reached at the end of the sampling interval, i.e. in $t_{k+1}$. The coupling strategy should assure that $H_d\left(x\left(t_{k+1}\right)\right)$ is lower or equal than $H_d\left({\xi }_k\right)$. Note that $\Delta H_d\left(t_k\right)$ in Figure 5 is coherent with such a desired behaviour, which means that the reference performances of the discrete-time scenario are obtained in the mixed discrete-continous-time case.

6. Conclusions and Future Activities

This paper presented a complete synthesis procedure for the design and implementation of energy-based digital controllers for port-Hamiltonian systems. The procedure consists of three main steps. First, a discrete-time port-Hamiltonian approximation of the plant is obtained by “replacing” the Hamiltonian gradient with a discrete gradient. This construction preserves the energy-balance relation of the continuous-time system, although it generally leads to implicit state equations. Second, an extension of the IDA-PBC methodology to the resulting class of discrete-time systems is presented. The proposed design reshapes the Hamiltonian function and modifies the interconnection and dissipative structures in order to assign a desired equilibrium and guarantee its asymptotic stability under suitable conditions. The third step addresses the interconnection between the digital controller and the continuous-time plant. A performance-preserving digital-to-analog interface is proposed, based on the solution of an optimization problem that computes the smallest correction to the nominal control action required to enforce a desired decrease of the closed-loop Hamiltonian over each sampling interval. A soft-constrained formulation is also considered to address possible feasibility issues. For linear systems with quadratic desired Hamiltonians, the optimization problem admits an explicit solution. The complete procedure was illustrated using a two-degree-of-freedom planar manipulator. The numerical results show that the proposed interface allows the mixed discrete-continuous-time closed-loop system to reproduce closely the response obtained in the reference discrete-time design scenario, while requiring only a small correction to the nominal control action, mainly during the initial transient.

Future work will focus on extending the proposed discrete-time IDA-PBC framework to systems with state-dependent structure matrices, for which solving the matching equation is more challenging. Further research will also address computationally efficient implementations of the coupling strategy for nonlinear systems, the treatment of feasibility issues over successive sampling intervals, and experimental validation on physical systems.