Total energy-shaping control for mechanical systems via Control-by-Interconnection

Joel Ferguson¹ ¹Joel Ferguson is with the School of Engineering, The University of Newcastle, Australia Email: [email protected]

Abstract

Application of IDA-PBC to mechanical systems has received much attention in recent decades, but its application is still limited by the solvability of the so-called matching conditions. In this work, it is shown that total energy-shaping control of under-actuated mechanical systems has a control-by-interconnection interpretation. Using this interpretation, alternate matching conditions are formulated that defines constraints on the added energy, rather then the total closed-loop energy. It is additionally shown that, for systems that are under-actuated degree one with the mass matrix depending on a single coordinate, the kinetic energy matching conditions resolve to ODEs which can be evaluated numerically. Using this approach controllers are proposed for the benchmark cart-pole and acrobot systems.

I INTRODUCTION

Energy-based methods for controlling nonlinear physical systems have been shown to be effective in a variety of physical domains [1]. Such methods consider the energy and structure of the system to be controlled to derive control strategies that exploit the natural system behaviours. Interconnection and Damping Assignment, Passivity-Based Control (IDA-PBC) is one such control methodology where the control input is designed such that the closed-loop can be interpreted as an alternate physical system with a different energy, interconnection and damping structure [2].

While IDA-PBC has been applied to a broad range of systems, particular attention has been given to mechanical systems which exhibit a rich canonical structure [3, 4, 5]. In the case of fully-actuated systems, IDA-PBC allows a user to arbitrarily modify the potential and kinetic energy of the closed-loop system [6], a process known as total energy shaping [3]. For under-actuated systems, however, application of IDA-PBC is limited by solutions that satisfy a set of PDEs, the so-called matching conditions. Much research effort has been committed to solving these equations, with solutions posed in several special cases [4, 5, 6]. This design methodology has been applied to a number of benchmark examples such as the cart-pole, acrobot, spider crane, amongst others.

Control-by-Interconnection (CbI) describes a sub-class of energy-based control methods that falls under the umbrella of IDA-PBC [7, 8]. Under this scheme, the controller is assumed to be a passive system that is interconnected with a passive plant to be controlled via the passive input-output pair. Casimirs, conserved quantities between the control sub-system and the plant, can be constructed to help shape the energy of the closed-loop system. It is known that potential energy shaping of fully-actuated mechanical systems falls into the class of CbI [7, 9]. Control of underactuated mechanical systems has been explored in the context of CbI by applying nonlinear PID controllers to both the standard and alternate passive outputs [10, 11, 12]. The idea of using PID for stabilisation of passive systems was formalised in [13] and a general characterisation of all passive outputs from a given system was characterised.

In this work, the connection between IDA-PBC and CbI for under-actuated mechanical systems is explored. Using the bond graph formalism (see [14] for introduction), a control sub-system is proposed that allows for shaping of the kinetic and potential energies of the closed-loop system. By representing the controller as a passive interconnection, the requisite matching conditions are reformulated in terms of the added mass and added potential energy. Equivalence between the CbI and IDA-PBC is then established in the case of mechanical systems by identifying Casimirs relating the controller states to those of the plant. Finally, using the reformulated matching conditions, it is shown that in the case that the mass matrix depends on only one coordinate that the kinetic energy matching conditions can be formulated as an ODE that can be evaluated numerically for implementation.

Notation. Function arguments are declared upon definition and are omitted for subsequent use. $0_{n\times m}$ denotes a $n\times m$ zeros matrix whereas $I_{n}$ denotes a $n\times n$ identity matrix. For mappings $\mathcal{H}:\mathbb{R}^{n}\to\mathbb{R}$ , we denote the transposed gradient as $\nabla\mathcal{H}:=\left(\frac{\partial\mathcal{H}}{\partial x}\right)^{\top}$ . For $P=P^{\top}\in\mathbb{R}^{n\times n}$ , $\lambda_{min}\left[P\right],\lambda_{max}\left[P\right]$ denotes the minimum and maximum (real) eigenvalues of P, respectively.

II BACKGROUND AND PROBLEM FORMULATION

In this section a number of key concepts necessary for the subsequent developments are briefly revised.

II-A Control-by-interconnection

In this work we consider input-state-output port-Hamiltonian systems (ISO-PHS) of the form

\begin{split}\dot{x}_{p}&=F_{p}(x_{p})\nabla_{x_{p}}H_{p}(x_{p})+G_{p}(x_{p})u_{p}\\ y_{p}&=G_{p}^{\top}(x_{p})\nabla_{x_{p}}H_{p}(x_{p})\end{split}

(1)

where $x_{p}\in\mathbb{R}^{p}$ is the state of the plant, $F_{p}(x_{p})\in\mathbb{R}^{p\times p}$ is the combined interconnection and damping matrix satisfying $F_{p}(x_{p})+F_{p}^{\top}(x_{p})\leq 0$ , $H_{p}(x_{p})\in\mathbb{R}$ is the Hamiltonian, $u_{p}\in\mathbb{R}^{m}$ is the input, $G_{p}(x_{p})\in\mathbb{R}^{p\times m}$ is the input mapping matrix and $y_{p}\in\mathbb{R}^{m}$ is the natural passive output corresponding to the input $u_{p}$ .

CbI assumes that the controller is a passive system that is interconnected with the plant (1) via a passive interconnection. In this work, we consider a controller subsystem with two input-output ports, described by the ISO-PHS

\begin{bmatrix}\dot{x}_{c}\\ -y_{c1}\\ -y_{c2}\end{bmatrix}=\underbrace{\begin{bmatrix}K_{11}(x_{c})&K_{12}(x_{c})&K_{13}(x_{c})\\ K_{21}(x_{c})&K_{22}(x_{c})&K_{23}(x_{c})\\ K_{31}(x_{c})&K_{32}(x_{c})&K_{33}(x_{c})\end{bmatrix}}_{:=K(x_{c})}\begin{bmatrix}\nabla_{x_{c}}H_{c}\\ u_{c1}\\ u_{c2}\end{bmatrix}

(2)

where $x_{c}\in\mathbb{R}^{c}$ is the state of the controller, $H_{c}(x_{c})\in\mathbb{R}$ is the controller Hamiltonian, $u_{c1},y_{c1}\in\mathbb{R}^{m}$ and $u_{c2},y_{c2}\in\mathbb{R}^{r}$ are passive input-output pairs and $K(x_{c})\in\mathbb{R}^{(p+m+r)\times(p+m+r)}$ satisfies $K(x_{c})+K^{\top}(x_{c})\leq 0$ [15].

The controller system (2) can be interconnected with the plant (1) via the passive interconnection

\begin{split}u_{p}&=-y_{c1}\\ u_{c1}&=y_{p},\end{split}

(3)

resulting in the closed-loop dynamics

\begin{split}\begin{bmatrix}\dot{x}_{p}\\ \dot{x}_{c}\\ -y_{c2}\end{bmatrix}&=\underbrace{\begin{bmatrix}F_{p}+G_{p}K_{22}G_{p}^{\top}&G_{p}K_{21}&G_{p}K_{23}\\ K_{12}G_{p}^{\top}&K_{11}&K_{13}\\ K_{32}G_{p}^{\top}&K_{31}&K_{33}\\ \end{bmatrix}}_{F_{cl}}\begin{bmatrix}\nabla_{x_{p}}H_{p}\\ \nabla_{x_{c}}H_{c}\\ u_{c2}\end{bmatrix}\end{split}

(4)

where $u_{c2},y_{c2}$ is a passive input-output pair to the interconnected system. Noting that $K+K^{\top}\leq 0$ , the closed-loop interconnection and damping structure $F_{cl}$ satisfies $F_{cl}+F_{cl}^{\top}\leq 0$ also.

In the case of stabilisation, the objective is to construct the plant functions $H_{c}(x_{c})$ , $K(x_{c})$ to ensure the existence of Casimirs which statically relate the controller states to functions of the plant states

x_{c}=f_{c}(x_{p}),

(5)

for $f_{c}(x)\in\mathbb{R}^{c}$ . The Casimir functions and controller initial conditions are then designed to assign a desirable minimum to the total energy function

W(x_{p})=H(x_{p})+H_{c}(x_{c})|_{x_{c}=f_{c}(x_{p})}.

(6)

It is noted that the Lyapunov candidate $W(x_{p})$ can be generalised to a function of $H,H_{c}$ and the Casimirs $x_{c}-f_{c}(x_{p})$ [15]. Methods to ensure the existence of and constructing Casimirs have been reported in [7] and the references therein.

II-B Underactuated mechanical systems

The primary objective of this work is to apply CbI to the class of underactuated mechanical systems, described by the dynamics

\begin{split}\begin{bmatrix}\dot{q}\\ \dot{p}\end{bmatrix}&=\begin{bmatrix}0_{n\times n}&I_{n}\\ -I_{n}&0_{n\times n}\end{bmatrix}\begin{bmatrix}\nabla_{q}H\\ \nabla_{p}H\end{bmatrix}+\begin{bmatrix}0_{n\times m}\\ G\end{bmatrix}u\\ H(q,p)&=\underbrace{\frac{1}{2}p^{\top}M^{-1}(q)p}_{:=T(q,p)}+V(q)\\ y&=G^{\top}\nabla_{p}H,\end{split}

(7)

where $q\in\mathbb{R}^{n},p\in\mathbb{R}^{n}$ are the configuration and momentum vectors, respectively, $u\in\mathbb{R}^{m}$ in the input, $M(q)=M^{\top}(q)>0$ is the inertia matrix and $y$ is the natural passive output corresponding the the input $u$ . The input mapping matrix $G$ is assumed to be constant and have the structure

G=\begin{bmatrix}I_{m}\\ 0_{(n-m)\times m}\end{bmatrix},

(8)

where $n-m<n$ is the degree of underactuation of the system¹¹1The assumed structure of $G$ requires that the first $m$ configuration coordinates are chosen to be collocated with the actuators. This class of dynamics falls into the broader class of ISO-PHS (1). For a more general input mapping matrix $\bar{G}(q)\in\mathbb{R}^{n\times m}$ , there exists a change of coordinates recovering the structure (8) if the columns of $\bar{G}(q)$ are involute.. The Hamiltonian $H(q,p)$ is the sum of the kinetic energy $T(q,p)$ and the potential energy $V(q)$ , which allows the gradient of $H$ with respect to $q$ to be written as

\nabla_{q}H(q,p)=\nabla_{q}T(q,p)+\nabla_{q}V(q).

(9)

A full-rank left-annihilator for the input mapping matrix (8) is defined as

G^{\perp}=\begin{bmatrix}0_{(n-m)\times m}&I_{(n-m)}\end{bmatrix}

(10)

which satisfies $G^{\perp}G=0_{(n-m)\times m}$ .

In the subsequent development, we will require an alternate representation of the gradient of the kinetic energy with respect to configuration $\nabla_{q}T(q,p)$ . Noting that the kinetic energy is quadratic in $p$ , the gradient $\nabla_{q}T(q,p)$ can always be factored into the form

\nabla_{q}T(q,p)=E(q,p)M^{-1}(q)p,

(11)

for some matrix $E(q,p)\in\mathbb{R}^{n\times n}$ . This has been previously noted in [16] using the Christoffel symbols. Note, however, that the matrix $E(q,p)$ is non-unique and in this work we will use the representation given by

\begin{split}\nabla_{q}T(q,p)&=\frac{1}{2}\frac{\partial^{\top}}{\partial q}\left(M^{-1}(q)p\right)p\\ &=\underbrace{\frac{1}{2}\frac{\partial^{\top}}{\partial q}\left(M^{-1}(q)p\right)M(q)}_{:=E(q,p)}M^{-1}(q)p.\end{split}

(12)

In constructing a CbI interpretation to total energy shaping it is useful to define a virtual input-output pair for the system (7) by defining the input

u_{v}=Gu,

(13)

which allows the system to written similarly to a fully-actuated system as

\begin{split}\begin{bmatrix}\dot{q}\\ \dot{p}\end{bmatrix}&=\begin{bmatrix}0_{n\times n}&I_{n}\\ -I_{n}&0_{n\times n}\end{bmatrix}\begin{bmatrix}\nabla_{q}H\\ \nabla_{p}H\end{bmatrix}+\begin{bmatrix}0_{n\times n}\\ I_{n}\end{bmatrix}u_{v}\\ y_{v}&=\nabla_{p}H,\end{split}

(14)

where $u_{v},y_{v}\in\mathbb{R}^{n}$ . From the definition (13), it is clear that any input $u_{v}$ must satisfy

G^{\perp}u_{v}=0_{m\times 1},

(15)

which will be ensured in subsequent control design. Assuming that (15) holds, the input $u$ can be described as a function of $u_{v}$ by

\begin{split}u&=G^{\top}u_{v}.\end{split}

(16)

The advantage of constructing the virtual input-output pair is that the virtual output now describes the full velocity vector

y_{v}=M^{-1}(q)p=\dot{q},

(17)

a property that will be exploited when shaping the potential energy.

II-C IDA-PBC for underactuated mechanical systems

IDA-PBC is a control design methodology whereby the control signal is designed such that the closed-loop dynamics have a port-Hamiltonian (pH) structure. When applied to underactuated mechanical systems, the target closed-loop dynamics have the structure

\begin{split}\begin{bmatrix}\dot{q}\\ \dot{p}\end{bmatrix}&=\begin{bmatrix}0_{n\times n}&M^{-1}(q)M_{d}(q)\\ -M_{d}(q)M^{-1}(q)&J_{2}(q,p)-GK_{d}G^{\top}\end{bmatrix}\begin{bmatrix}\nabla_{q}H_{d}\\ \nabla_{p}H_{d}\end{bmatrix}\\ H_{d}&(q,p)=\frac{1}{2}p^{\top}M_{d}^{-1}(q)p+V_{d}(q),\end{split}

(18)

where $M_{d}(q)=M_{d}^{\top}(q)>0,V_{d}(q)$ are the desired closed-loop inertia matrix and potential energies, respectively, $J_{2}(q,p)=-J_{2}^{\top}(q,p)$ is skew-symmetric and $K_{d}=K_{d}^{\top}\geq 0$ is a tuning parameter used for damping injection. If $V_{d}$ is minimised at the target configuration, $H_{d}$ qualifies as a Lyapunov function for the closed-loop system [4].

The complexity of applying IDA-PBC to underactuated system is satisfying the so-called matching conditions. This conditions requires that the dynamics of the open-loop and closed-loop systems must agree on the spaces perpendicular to the control signal. The structure chosen fo the closed-loop system in (18) ensures that the dynamics of $q$ agree with (7). Comparing the dynamics of $p$ results in the condition

\begin{split}G^{\perp}\left\{\nabla_{q}H-M_{d}(q)M^{-1}(q)\nabla_{q}H_{d}-J_{2}(q,p)\nabla_{p}H_{d}\right\}\\ =0_{(n-m)\times 1},\end{split}

(19)

which defines a PDE that should be solved for $M_{d}(q),V_{d}(q),J_{2}(q,p)$ . Noting the structure of the Hamiltonians, this PDE can be separated into the components involving $p$ , and those that do not by

\begin{split}G^{\perp}\left\{\nabla_{q}T-M_{d}(q)M^{-1}(q)\nabla_{q}T_{d}-J_{2}(q,p)M_{d}^{-1}(q)p\right\}\\ =0_{(n-m)\times 1}\\ G^{\perp}\left\{\nabla_{q}V-M_{d}(q)M^{-1}(q)\nabla_{q}V_{d}\right\}\\ =0_{(n-m)\times 1}.\end{split}

(20)

These expressions are known as the kinetic energy and potential energy matching equations, respectively.

II-D Contributions

The objective of this work is to construct a CbI interpretation of IDA-PBC when applied to underactuated mechanical systems. The contributions of this work are threefold:

C.1

ISO-PHS with Casimirs that statically relate states are considered and a closed-form solution to remove the Casimirs by reducing the dimension of the state vector is proposed. This solution can be applied to the closed-loop dynamics of CbI implementations of the form (4) to describe the resulting dynamics as a function of $x_{p}$ only.
C.2

A CbI controller for underactuated mechanical system of the form (2) is proposed and the resulting closed-loop is shown to be equivalent to the well-known dynamics (18). The CbI interpretation generates alternate matching conditions to the expressions (20), describing constraints on the added mass and added potential energy.
C.3

Using the alternate matching conditions, it is shown that the kinetic energy matching equations reduce to ODEs in the special case of underactuation degree one where the mass matrix is a function of only one configuration coordinate. It is demonstrated that numerical methods can be utilised in such cases to avoid solving these expressions analytically.

II-E Related works

Significant attention has been given to solving the matching equations (20) in recent decades. In [3], [17] it was shown that is the system is under-actuated degree one and the mass matrix depends on a single un-actuated coordinate, the kinetic energy matching condition can be simplified to an ODE. Using a novel parametrisation of $J_{2}(q,p)$ , a general solution for under-actuated degree one system was proposed in [4] under the assumption that the mass matrix depends only on the actuated coordinates. This approach was extended in [18] using a momentum transformation to simplify the matching equations. More recently, solutions to the potential energy matching equations were considered in [19] under the assumption that the mass matrix and potential energy functions were dependent on only one variable. Finally, a general solution to the special case of 2 degree-of-freedom system was proposed in [6]. The studies [20], [21] considered the effects of friction on the closed-loop stability using IDA-PBC.

In recent works, alternate approaches to constructing solutions to the matching equations have been explored. The existence of conservative forces that cannot be factorised into a skew-symmetric matrix $J_{2}(q,p)$ was investigated in [22], which resulted in alternate matching equations. Implicit system representations were used in [23] to construct solutions in an over-parameterised space where the closed-loop dynamics were subject to constraints. By working in the larger dimension, a solution to a under-actuated degree 2 crane system was proposed. Pfaffian differential equations were utilised in [24] which resulted in the kinetic energy PDEs being converted to an alternate form which admits simpler solutions.

Some authors have investigated the possibility of avoiding the matching equations altogether by considering the control signal to be a CbI. The work [10] relied on a Lagrangian structure and several technical assumptions to verify the existence of a second passive output corresponding to the input $u$ . Using this second output, a stabilising control was designed that ensured stability without requiring a solution to the matching PDEs. A similar approach was proposed in [12, 11] where a second passive output was utilised and the control assumed to have a PID structure.

III CASIMIR REDUCTION

In this section, a method for reducing the state dimension of ISO-PHS with Casimirs is derived. The reduction method applies to general ISO-PHS with Casimirs and can be directly applied to the resulting closed-loop dynamics of CbI schemes of the form (4) to describe the system as a function of $x_{p}$ only. In the sequel, the reduction method will be used to show equivalence between the CbI controller for underactuated mechanical systems and the IDA-PBC dynamics (18). Before introducing the state reduction solution, a useful lemma is required.

Lemma 1

Consider a square block matrix of arbitrary dimension

A=\begin{bmatrix}A_{11}&A_{12}\\ A_{21}&A_{22}\end{bmatrix}

(21)

and assume that $A_{22}$ is invertible. If the symmetric component of $A$ is negative semi-definite, $A+A^{\top}\leq 0$ , the symmetric component of the Schur complement

A_{11}-A_{12}A_{22}^{-1}A_{21}

(22)

is negative semi-definite also.

Proof:

First note that the Schur complement of $X$ can be computed by

\begin{split}\begin{bmatrix}I&-A_{21}^{\top}A_{22}^{-\top}\end{bmatrix}A\begin{bmatrix}I\\ -A_{22}^{-1}A_{21}\end{bmatrix}=A_{11}-A_{12}A_{22}^{-1}A_{21}.\end{split}

(23)

The symmetric component of this expression is negative semi-definite as $A+A^{\top}\leq 0$ . ∎

The solution for reducing the dimension of ISO-PHS which exhibit Casimirs is now introduced. This development applies to systems of the form

\begin{split}\begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\\ -y\end{bmatrix}&=\underbrace{\begin{bmatrix}F_{11}(x)&F_{12}(x)&F_{13}(x)\\ F_{21}(x)&F_{22}(x)&F_{23}(x)\\ F_{31}(x)&F_{32}(x)&F_{33}(x)\end{bmatrix}}_{F(x)}\begin{bmatrix}\nabla_{x_{1}}H\\ \nabla_{x_{2}}H\\ u\end{bmatrix}\end{split}

(24)

where $x\in\mathbb{R}^{p+c}$ is the state of the system which has been partitioned into $x_{1}\in\mathbb{R}^{p},x_{2}\in\mathbb{R}^{c}$ , $H(x_{1},x_{2})\in\mathbb{R}$ is the Hamiltonian, $F(x)\in\mathbb{R}^{(p+c)\times(p+c)}$ is the full-rank interconnection and damping matrix satisfying $F(x)+F^{\top}(x)\leq 0$ , $u\in\mathbb{R}^{m}$ is the input and $y\in\mathbb{R}^{m}$ is the corresponding passive output. It is assumed that the system contains a Casimir and the states have been partitioned such that the Casimir can be written as

x_{2}=f_{c}(x_{1}),

(25)

where $f_{c}(x_{1})\in\mathbb{R}^{c}$ is differentiable.

The first step in constructing a minimal system representation is defining a new set of coordinates given by

w=x_{2}-f_{c}(x_{1})=0_{c\times 1},

(26)

which is identically equal to zero by construction. The system (24) can be described in the coordinates $(x_{1},w)$ by

\begin{split}\begin{bmatrix}\dot{x}_{1}\\ -y\\ \dot{w}\end{bmatrix}=&\begin{bmatrix}I_{p}&0_{p\times c}&0_{p\times m}\\ 0_{m\times p}&0_{m\times c}&I_{m}\\ -\frac{\partial f_{c}}{\partial x_{1}}&I_{c}&0_{c\times m}\end{bmatrix}\begin{bmatrix}\dot{x}_{1}\\ \dot{x}_{2}\\ -y\end{bmatrix}\\ =&\begin{bmatrix}I_{p}&0_{p\times c}&0_{p\times m}\\ 0_{m\times p}&0_{m\times c}&I_{m}\\ -\frac{\partial f_{c}}{\partial x_{1}}&I_{c}&0_{c\times m}\end{bmatrix}F\\ &\times\begin{bmatrix}I_{p}&0_{p\times m}&-\frac{\partial^{\top}f_{c}}{\partial x_{1}}\\ 0_{c\times p}&0_{c\times m}&I_{c}\\ 0_{m\times p}&I_{m}&0_{m\times c}\end{bmatrix}\begin{bmatrix}\nabla_{x_{1}}H_{r}\\ u\\ \nabla_{w}H_{r}\end{bmatrix}\\ =&\underbrace{\begin{bmatrix}\bar{F}_{11}&\bar{F}_{12}&\bar{F}_{13}\\ \bar{F}_{21}&\bar{F}_{22}&\bar{F}_{23}\\ \bar{F}_{31}&\bar{F}_{32}&\bar{F}_{33}\end{bmatrix}}_{\bar{F}}\begin{bmatrix}\nabla_{x_{1}}H_{r}\\ u\\ \nabla_{w}H_{r}\end{bmatrix},\\ \end{split}

(27)

where

\begin{split}H_{r}(x_{1},w):=&H(x_{1},w+f_{c}(x_{1}))\\ \bar{F}_{11}(x_{1})=&F_{11}(x)|_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{12}(x_{1})=&F_{13}(x)|_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{13}(x_{1})=&F_{12}(x)-F_{11}(x)\frac{\partial^{\top}f_{c}}{\partial x_{1}}\bigg{|}_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{21}(x_{1})=&F_{31}(x)|_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{22}(x_{1})=&F_{33}(x)|_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{23}(x_{1})=&F_{32}(x)-F_{31}(x)\frac{\partial^{\top}f_{c}}{\partial x_{1}}\bigg{|}_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{31}(x_{1})=&F_{21}(x)-\frac{\partial f_{c}}{\partial x_{1}}F_{11}(x)\bigg{|}_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{32}(x_{1})=&F_{23}(x)-\frac{\partial f_{c}}{\partial x_{1}}F_{13}(x)\bigg{|}_{x_{2}=f_{c}(x_{1})}\\ \bar{F}_{33}(x_{1})=&F_{22}(x)-F_{21}(x)\frac{\partial^{\top}f_{c}}{\partial x_{1}}-\frac{\partial f_{c}}{\partial x_{1}}F_{12}(x)\\ &+\frac{\partial f_{c}}{\partial x_{1}}F_{11}(x)\frac{\partial^{\top}f_{c}}{\partial x_{1}}\bigg{|}_{x_{2}=f_{c}(x_{1})}.\end{split}

(28)

Recalling that $w$ is identically equal to zero, $\dot{w}$ is also equal to zero. Consequently, the final row of (27) is a constraint that needs to be resolved to construct a minimal system representation. There are two methods of reduction that will be considered. Firstly, in the transformed coordinates (27) it can occur that one or more columns of $\bar{F}_{\star 3}$ is identically zero. Without loss of generality, it is assumed that the first $d$ columns of $\bar{F}_{\star 3}$ are equal to zero. To remove the zero rows, the full-rank matrix $B$ is defined as

B=\begin{bmatrix}0_{d\times(c-d)}\\ I_{(c-d)}\end{bmatrix},

(29)

which acts to select the non-zero columns of $\bar{F}_{\star 3}$ . The zero columns and corresponding rows are removed from (27) by

\begin{split}\begin{bmatrix}\dot{x}_{1}\\ -y\\ B^{\top}\dot{w}\end{bmatrix}=&\underbrace{\begin{bmatrix}\bar{F}_{11}&\bar{F}_{12}&\bar{F}_{13}B\\ \bar{F}_{21}&\bar{F}_{22}&\bar{F}_{23}B\\ B^{\top}\bar{F}_{31}&B^{\top}\bar{F}_{32}&B^{\top}\bar{F}_{33}B\end{bmatrix}}_{:=\bar{F}_{B}}\begin{bmatrix}\nabla_{x_{1}}H_{r}\\ u\\ B^{\top}\nabla_{w}H_{r}\end{bmatrix},\end{split}

(30)

which does not modify the system dynamics. Note that as $\bar{F}+\bar{F}^{\top}\leq 0$ , $\bar{F}_{B}+\bar{F}_{B}^{\top}\leq 0$ also. Using this representation, a method for resolving the remaining constraint equations is presented under the assumption that $B^{\top}\bar{F}_{33}B$ is full rank.

Proposition 1

Consider the pH system (24) with Casimir (25). If the matrix $B^{\top}\bar{F}_{33}B$ is full-rank for all $x_{1}$ the system can be described by a reduced-order model

\begin{split}\begin{bmatrix}\dot{x}_{1}\\ -y\end{bmatrix}&=F_{r}(x_{1})\begin{bmatrix}\nabla_{x_{1}}H_{r}\\ u\end{bmatrix},\end{split}

(31)

where

\begin{split}H_{r}(x_{1})=&H\left(x_{1},f_{c}(x_{1})\right)\\ F_{r}(x_{1})=&\begin{bmatrix}\bar{F}_{11}&\bar{F}_{12}\\ \bar{F}_{21}&\bar{F}_{22}\\ \end{bmatrix}\\ &-\begin{bmatrix}\bar{F}_{13}B\\ \bar{F}_{23}B\\ \end{bmatrix}(B^{\top}\bar{F}_{33}B)^{-1}\begin{bmatrix}B^{\top}\bar{F}_{31}&B^{\top}\bar{F}_{32}\end{bmatrix}\end{split}

(32)

and $\bar{F}(x_{1})$ satisfies $F_{r}(x_{1})+F_{r}^{\top}(x_{1})\leq 0$ .

Proof:

The expression $B^{\top}\dot{w}$ , defined in (30), is identically equal to $0_{(c-d)\times 1}$ by construction. Note, however, that the gradient $B^{\top}\nabla_{w}H_{r}$ is not necessarily equal to zero. Assuming that $B^{\top}\bar{F}_{33}B$ is full rank, the expression $B^{\top}\nabla_{w}H_{r}$ can be described as

\begin{split}B^{\top}\nabla_{w}H_{r}=&-(B^{\top}\bar{F}_{33}B)^{-1}\begin{bmatrix}B^{\top}\bar{F}_{31}&B^{\top}\bar{F}_{32}\end{bmatrix}\begin{bmatrix}\nabla_{x_{1}}H_{r}\\ u\end{bmatrix}\end{split}

(33)

Substituting this expression into the dynamics (30) resolves to the reduced dynamics (31).

To verify that $F_{r}+F_{r}^{\top}\leq 0$ , note that $F_{r}$ is the Schur complement of $\bar{F}_{B}$ which satisfies $\bar{F}_{B}+\bar{F}_{B}^{\top}\leq 0$ . It follows that $F_{r}+F_{r}^{\top}\leq 0$ by application of Lemma 1. ∎

Proposition 1 showed that an ISO-PHS that exhibits a Casimir function can be described in a reduced state-space. The class of dynamics that are derived from application of CbI (4) that result in a Casimir of the form (5) falls into the class of systems (24). The following Corollary tailors the Casimir reduction for this important sub-class of dynamics.

Corollary 1

If the closed-loop dynamics of a CbI scheme (4) exhibit a Casimir of the form (5), the system can be equivalently expressed in the form (31) where

\begin{split}x_{1}=&x_{p}\\ H_{r}(x_{p})=&H_{p}(x_{p})+H_{c}\left(f_{c}(x_{p})\right)\\ \bar{F}_{11}=&F_{p}+G_{p}K_{22}G_{p}^{\top}\\ \bar{F}_{12}=&G_{p}K_{23}\\ \bar{F}_{13}=&G_{p}K_{21}-\left[F_{p}+G_{p}K_{22}G_{p}^{\top}\right]\frac{\partial^{\top}f_{c}}{\partial x_{p}}\\ \bar{F}_{21}=&K_{32}G_{p}^{\top}\\ \bar{F}_{22}=&K_{33}\\ \bar{F}_{23}=&K_{31}-K_{32}G_{p}^{\top}\frac{\partial^{\top}f_{c}}{\partial x_{p}}\\ \bar{F}_{31}=&K_{12}G_{p}^{\top}-\frac{\partial f_{c}}{\partial x_{p}}\left[F_{p}+G_{p}K_{22}G_{p}^{\top}\right]\\ \bar{F}_{32}=&K_{13}-\frac{\partial f_{c}}{\partial x_{p}}G_{p}K_{23}\\ \bar{F}_{33}=&K_{11}-K_{12}G_{p}^{\top}\frac{\partial^{\top}f_{c}}{\partial x_{p}}-\frac{\partial f_{c}}{\partial x_{p}}G_{p}K_{21}\\ &+\frac{\partial f_{c}}{\partial x_{p}}\left[F_{p}+G_{p}K_{22}G_{p}^{\top}\right]\frac{\partial^{\top}f_{c}}{\partial x_{p}}\end{split}

(34)

and $B$ is suitably chosen as per (29) using the expressions $\bar{F}_{\star 3}$ . The arguments have been dropped from the definitions of $\bar{F}_{\star\star}(x_{p})$ for the sake of readability.

Proof:

The result follows from direct application of Proposition 1 to the dynamics (4). ∎

IV CONTROL-BY-INTERCONNECTION FOR MECHANICAL SYSTEMS

In this section, a control-by-interconnection scheme for under-actuated mechanical systems is presented. A dynamic 2-port control system is introduced with the intention that it will be interconnected to the plant (14) via one of the ports. The controller states are constructed to be statically related to the plant states after interconnection, resulting in Casimirs. By applying Proposition 1, the closed-loop dynamic are defined in a reduced space in which the dynamics coincide with standard total energy-shaping control (18).

The proposed CbI scheme is shown in Figure 1. The intention of this control subsystem is to interconnect with the plant (14) via the $u_{c1},y_{c1}$ power port.

Refer to caption — Figure 1: Total energy shaping as a CbI for under-actuated mechanical systems.

The second input $u_{c2}$ is available for subsequent control design, such as damping injection. The terms $M(q_{a_{2}}),E(q_{a2},p_{a})$ are the plant mass matrix (7) and factorisation of the kinetic energy gradient (12) evaluated at the controller states whereas $J(q_{a2},p_{a})=-J(q_{a2},p_{a})^{\top}\in\mathbb{R}^{n\times n}$ is a skew-symmetric matrix to be chosen. The three-port storage element $H_{a}(q_{a1},q_{a2},p_{a})$ has states $q_{a1},q_{a2},p_{a}\in\mathbb{R}^{n}$ and energy function similar to mechanical systems,

\begin{split}H_{a}(q_{a1},q_{a2},p_{a})&=\underbrace{\frac{1}{2}p_{a}^{\top}M_{a}^{-1}(q_{a2})p_{a}}_{:=T_{a}(q_{a2},p_{a})}+\underbrace{V_{d}(q_{a2})-V(q_{a1})}_{:=V_{a}(q_{a1},q_{a2})},\end{split}

(35)

where $M_{a}^{-1}(q_{a2})$ is the inverse added mass, $T_{a}(q_{a2},p_{a})$ is the added kinetic energy, $V_{d}(q_{a2})$ is the desired closed-loop potential energy, $V(q_{a1})$ is the plant potential energy function (7) evaluated at the plant state $q_{a1}$ and $V_{a}(q_{a1},q_{a2})$ is the total added potential energy. It is important to note that, although $M_{a}^{-1}(q_{a2})$ is represented as a matrix inverse, it need not be invertible nor positive. Indeed, it will be shown in subsequent developments that the key requirement is that

M_{d}^{-1}(q):=M^{-1}(q)+M_{a}^{-1}(q)

(36)

should be positive definite.

In subsequent analysis it will be shown that the interconnection of the control system with the plant (14) via the interconnection

\begin{split}u_{v}&=-y_{c1}\\ u_{c1}&=y_{v},\end{split}

(37)

yields Casimirs

\begin{bmatrix}q_{a1}\\ q_{a2}\\ p_{a}\end{bmatrix}=\underbrace{\begin{bmatrix}I_{n}&0_{n\times n}\\ I_{n}&0_{n\times n}\\ 0_{n\times n}&I_{n}\end{bmatrix}\begin{bmatrix}q\\ p\end{bmatrix}}_{f_{c}(x_{p})}.

(38)

Assuming the Casimirs exist, some intuition regarding the construction of the control system in Figure 1 can be provided. Both $q_{a1},q_{a2}$ were constructed to be equal to $q$ . Firstly $\dot{q}_{a1}$ is equal to $\dot{q}$ by interconnection with the plant virtual output $y_{v}$ via a 1-junction. To verify a similar relation for $q_{a2}$ , assume that $q_{a2}=q,p_{a}=p$ holds which results in $\nabla_{p_{a}}H_{a}=M_{a}^{-1}(q)p$ and $u_{c1}+\nabla_{p_{a}}H_{a}=M_{d}^{-1}p$ . With this in mind, the transformer can be seen to reconstruct the velocity $\dot{q}=M^{-1}(q)p$ for the bottom 1-junction, resulting in $\dot{q}_{a1}=\dot{q}$ .

To construct a Casimir $p_{a}=p$ , first note from (7) and (12) that the plant momentum dynamics can be expressed as

\dot{p}=-\nabla_{q}V-E(q,p)M^{-1}(q)p+u_{v}.

(39)

The control structure acts to remove these forces from the plant via the right side of the control structure and re-introduce them via the top 0-junction where they are shared with the dynamics of $p_{a}$ . The $\dot{q}_{a1}$ bond acts to cancel the gravity term from the plant $-\nabla_{q}V$ . Recalling that the bottom 1-junction has flow equal to $M^{-1}(q)p$ , the right-side gyrator cancels the term $-E(q,p)M^{-1}(q)p$ from the plant. The left-side gyrator then re-introduces the force $-E(q,p)M^{-1}(q)p$ via the top 0-junction where it is shared between $\dot{p}$ and $\dot{p}_{a}$ , establishing the desired Casimir.

The claimed Casimir (38) is now formalised in the following Proposition. For this development, note that the gradients of the added energy $H_{a}(\cdot)$ satisfy

\begin{split}\nabla_{q_{a1}}H_{a}&=-\nabla_{q_{a1}}V\\ \nabla_{q_{a2}}H_{a}&=\nabla_{q_{a2}}T_{a}+\nabla_{q_{a2}}V_{d}\\ \nabla_{q_{a2}}T_{a}&=\frac{1}{2}\frac{\partial^{\top}}{\partial q_{a2}}\left(M_{a}^{-1}(q_{a2})p_{a}\right)p_{a}\\ \nabla_{p_{a}}H_{a}&=M_{a}^{-1}(q_{a2})p_{a}\end{split}

(40)

and the expressions $A(\cdot),B(\cdot),C(\cdot)$ in Figure 1 can be evaluated as

\begin{split}A(q_{a2},p_{a},u_{c1})=&M^{-1}(q_{a2})M_{d}(q_{a2})\left[u_{c1}+\nabla_{p_{a}}H_{a}\right]\\ B(q_{a2},p_{a},u_{c1})=&\nabla_{q_{a2}}H_{a}-E^{\top}(q_{a2},p_{a})\nabla_{p_{a}}H_{a}\\ &-M(q_{a2})J(q_{a2},p_{a})M_{d}(q_{a2})\\ &\times\left[u_{c1}+\nabla_{p_{a}}H_{a}\right]\\ C(q_{a2},p_{a},u_{c1})=&M_{d}(q_{a2})M^{-1}(q_{a2})\nabla_{q_{a2}}H_{a}\\ &-M_{d}(q_{a2})M^{-1}(q_{a2})E^{\top}(q_{a2},p_{a})\nabla_{p_{a}}H_{a}\\ &-M_{d}(q_{a2})J(q_{a2},p_{a})M_{d}(q_{a2})\\ &\times\left[u_{c1}+\nabla_{p_{a}}H_{a}\right],\end{split}

(41)

with $M_{d}(\cdot)$ defined in (36). To ensure that the Casimir exists, a number of requirements are imposed on the selection of the added inverse mass $M_{a}^{-1}(q)$ and closed-loop potential energy $V_{d}(q)$ which are equivalent of the standard matching conditions used in IDA-PBC (20).

Proposition 2

Consider the control system in Figure 1 and assume that it is interconnected to the plant (14) via the interconnection (37). If $M_{a}^{-1}(q_{a}),V_{d}(q_{a})$ are chosen such that

G^{\perp}C(q_{a2},p_{a},u_{c1})|_{q_{a2}=q,p_{a}=p}=G^{\perp}\nabla_{q}V

(42)

and the controller states are initialised as $q_{a1}(0)=q_{a2}(0)=q(0)$ , $p_{a}(0)=p(0)$ , the Casimir (38) holds for all time.

Proof:

Consider that at some time instant $T$ (38) holds, implying that

\begin{split}q_{a1}(T)=q_{a2}(T)=q(T),\ p_{a}(T)=p(T).\end{split}

(43)

It is shown that if (42) is satisfied, then the derivatives of the states also agree

\begin{split}\dot{q}_{a1}(T)=\dot{q}_{a2}(T)=\dot{q}(T),\ \dot{p}_{a}(T)=\dot{p}(T),\end{split}

(44)

establishing the existence of a Casimir for all future time.

We proceed by first establishing the relationship for the configuration vector. From (37) and (14) $u_{c}=M^{-1}(q)p$ which establishes $\dot{q}_{a1}(T)=\dot{q}(T)$ . The input $u_{c}=M^{-1}(q)p$ is substituted into $A(\cdot)$ (41) to find

\begin{split}A|_{t=T}&=M^{-1}(q_{a2})M_{d}(q_{a2})\left[M^{-1}(q)p+M_{a}^{-1}(q_{a2})p_{a}\right]|_{t=T}\\ &=M^{-1}(q)p|_{t=T}\\ &=\dot{q}|_{t=T},\end{split}

(45)

confirming that $\dot{q}_{a2}(T)=\dot{q}(T)$ .

Next we consider the behaviour of the momentum states. First note that, from the bond graph in Figure 1 and the definition (16), the plant input $u_{v}$ is given by

\begin{split}u&=G^{\top}u_{v}\\ &=G^{\top}\left[Gu_{c2}-C(q_{a2},p_{a},u_{c1})+\nabla_{q_{a1}}V(q_{a1})\right]\\ &=u_{c2}-G^{\top}\left[C(q_{a2},p_{a},u_{c1})-\nabla_{q_{a1}}V(q_{a1})\right].\end{split}

(46)

Using the control definition (46) and the condition (42), the plant dynamics (7) can be expanded as

\begin{split}\dot{p}=&-\nabla_{q}T(q,p)-\begin{bmatrix}G^{\top}\nabla_{q}V(q)\\ G^{\perp}\nabla_{q}V(q)\end{bmatrix}+Gu\\ =&-\nabla_{q}T(q,p)-\begin{bmatrix}G^{\top}\nabla_{q}V(q)\\ G^{\perp}\nabla_{q}V(q)\end{bmatrix}\\ &+G\left\{u_{c2}-G^{\top}\left[C(q_{a2},p_{a},u_{c1})-\nabla_{q_{a1}}V(q_{a1})\right]\right\}\\ =&-\nabla_{q}T(q,p)+Gu_{c2}\\ &-\begin{bmatrix}G^{\top}\left\{C(q_{a2},p_{a},u_{c1})+\nabla_{q}V(q)-\nabla_{q_{a1}}V(q_{a1})\right\}\\ G^{\perp}\nabla_{q}V(q)\end{bmatrix}\\ \end{split}

(47)

Note that at time $T$ , $\nabla_{q}V(q)|_{t=T}=\nabla_{q_{a1}}V(q_{a1})|_{t=T}$ . Additionally recall the assumption (42) which allows the simplification

\begin{split}\dot{p}|_{t=T}=&-\nabla_{q}T(q,p)-C(q_{a2},p_{a},u_{c1})+Gu_{c2}.\\ \end{split}

(48)

Recalling the identity (45), the dynamics of $p_{a}$ at time $T$ can be expanded to

\begin{split}\dot{p}_{a}&=Gv-E(q_{a2},p_{a})A(\cdot)|_{t=T}-C(q_{a2},p_{a},u_{c1})|_{t=T}\\ &=Gv-\nabla_{q}T(q,p)|_{t=T}-C(q_{a2},p_{a},u_{c1})|_{t=T},\end{split}

(49)

which agrees with (48). As (48) and (49) agree at time $T$ , (44) is verified for the momentum states. If at the initial time $t=0$ we have $q_{a}(0)=q(0)$ , $p_{a}(0)=p(0)$ , it follows that $q_{a}(t)=q(t)$ and $p_{a}(t)=p(t)$ for all time via integration, completing the proof. ∎

Proposition 2 has established that the Casimir (38) holds under some technical assumptions that will be verified in subsequent design. Before proceeding, we note that the control subsystem in Figure 1 can be written in the form (2) with

\begin{split}x_{c}&=\begin{bmatrix}q_{a1}^{\top}&q_{a2}^{\top}&p_{a}^{\top}\end{bmatrix}^{\top}\\ H_{c}(q_{a1},q_{a2},p_{a})&=H_{a}(q_{a1},q_{a2},p_{a})\\ K_{11}(q_{a2},p_{a})&=\begin{bmatrix}0_{n\times n}&0_{n\times n}&0_{n\times n}\\ 0_{n\times n}&0_{n\times n}&M^{-1}M_{d}\\ 0_{n\times n}&-M_{d}M^{-1}&D-D^{\top}+M_{d}JM_{d}\end{bmatrix}\\ K_{12}(q_{a2},p_{a})&=\begin{bmatrix}I_{n}\\ M^{-1}M_{d}\\ M_{d}JM_{d}-D^{\top}\end{bmatrix}\\ K_{13}&=\begin{bmatrix}0_{n\times m}\\ 0_{n\times m}\\ G\end{bmatrix}\\ K_{21}(q_{a2},p_{a})&=\begin{bmatrix}-I_{n}&-M_{d}M^{-1}&D+M_{d}JM_{d}\\ \end{bmatrix}\\ K_{31}&=\begin{bmatrix}0_{m\times n}&0_{m\times n}&-G^{\top}\end{bmatrix}\\ K_{22}(q_{a2},p_{a})&=M_{d}JM_{d}\\ K_{23}&=G\\ K_{32}&=-G^{\top}\\ K_{33}&=0_{m\times m}\\ D(q_{a2},p_{a})&=M_{d}M^{-1}E^{\top}.\end{split}

(50)

In the subsequent developments it is assumed that the requisite (42) of Proposition 2 holds, implying $q_{a1}(t)=q_{a2}(t)=q(t)$ , $p_{a}(t)=p(t)$ . Condition (42) will be verified by choice of $M_{a}^{-1}$ and $V_{d}$ . Assuming the Casimir holds, the expressions for $A(\cdot),B(\cdot),C(\cdot)$ in (41) can be simplified to

\begin{split}A(q,p)=&M^{-1}(q)p\\ B(q,p)=&\nabla_{q}T_{a}(q,p)+\nabla_{q}V_{d}(q,p)-E^{\top}(q,p)M_{a}^{-1}(q)p\\ &-M(q)J(q,p)p\\ C(q,p)=&M_{d}(q)\left\{M^{-1}(q)\nabla_{q}T_{a}(q,p)+M^{-1}(q)\nabla_{q}V_{d}(q,p)\right.\\ &\left.-M^{-1}(q)E^{\top}(q,p)M_{a}^{-1}(q)p-M^{-1}(q)J(q,p)\right\}\\ \end{split}

(51)

Recalling the definition of $\nabla_{q_{a}}T_{a}$ in (40), it is noted that $C(\cdot)$ contains some terms which are quadratic in $p$ and some that are functions of $q$ only. The function $C(\cdot)$ is divided into

\begin{split}C(q,p)=&C_{KE}(q,p)+C_{PE}(q)\\ C_{KE}(q,p)=&\underbrace{\left[M_{a}^{-1}(q)+M^{-1}(q)\right]^{-1}}_{M_{d}(q)}\left\{Y(q,p)-J(q,p)\right\}p\\ C_{PE}(q)&=\left[M_{a}^{-1}(q)+M^{-1}(q)\right]^{-1}M^{-1}(q)\nabla_{q}V_{d},\end{split}

(52)

where $KE$ represents kinetic energy, $PE$ represents potential energy and $Y$ is defined as

\begin{split}Y(q,p)=&\frac{1}{2}M^{-1}(q)\frac{\partial^{\top}}{\partial q}\left(M_{a}^{-1}(q)p\right)\\ &-\frac{1}{2}\frac{\partial}{\partial q}\left(M^{-1}(q)p\right)M_{a}^{-1}(q).\end{split}

(53)

As $Y$ is linear in $p$ it can be written as

\begin{split}Y(q,p)=\sum_{i=1}^{n}p_{i}Y^{i}(q),\end{split}

(54)

where

\begin{split}Y^{i}(q)=&\frac{1}{2}M^{-1}(q)\frac{\partial^{\top}}{\partial q}\left(M_{a}^{-1}(q)e_{i}\right)\\ &-\frac{1}{2}\frac{\partial}{\partial q}\left(M^{-1}(q)e_{i}\right)M_{a}^{-1}(q).\end{split}

(55)

The key constraint for control design is choosing $M_{a}^{-1}(q),V_{d}(q)$ satisfying the matching condition (42). From the definition of $C(\cdot)$ in (51), the constraint equation is a function of both $M_{a}^{-1}(q)$ and $\left[M_{a}^{-1}(q)+M^{-1}(q)\right]^{-1}$ , making direct design of this matrix difficult. To simplify the design process, an alternate characterisation of (42) is introduced.

In the following proposition, the inverse mass matrix, inverse added mass matrix, interconnection matrix and $Y(\cdot)$ are partitioned as

\begin{split}\begin{bmatrix}m_{11}(q)&m_{21}^{\top}(q)\\ m_{21}(q)&m_{22}(q)\end{bmatrix}&=\begin{bmatrix}G^{\top}\\ G^{\perp}\end{bmatrix}M^{-1}(q)\begin{bmatrix}G&G^{\perp\top}\end{bmatrix}\\ \begin{bmatrix}m_{a11}(q)&m_{a21}^{\top}(q)\\ m_{a21}(q)&m_{a22}(q)\end{bmatrix}&=\begin{bmatrix}G^{\top}\\ G^{\perp}\end{bmatrix}M_{a}^{-1}(q)\begin{bmatrix}G&G^{\perp\top}\end{bmatrix}\\ \begin{bmatrix}J_{11}(q,p)&-J_{21}^{\top}(q,p)\\ J_{21}(q,p)&J_{22}(q,p)\end{bmatrix}&=\begin{bmatrix}G^{\top}\\ G^{\perp}\end{bmatrix}J(q,p)\begin{bmatrix}G&G^{\perp\top}\end{bmatrix}\\ \begin{bmatrix}Y_{11}(q,p)&Y_{12}(q,p)\\ Y_{21}(q,p)&Y_{22}(q,p)\end{bmatrix}&=\begin{bmatrix}G^{\top}\\ G^{\perp}\end{bmatrix}Y(q,p)\begin{bmatrix}G&G^{\perp\top}\end{bmatrix}.\end{split}

(56)

Using the above definitions, an alternate characterisation of (42) is presented.

Proposition 3

The matching condition (42) is satisfied if:

•

The added mass matrix $M_{a}^{-1}(q)$ is chosen such that

$\begin{split}D(q)\left[Y^{i}(q)+Y^{i\top}(q)\right]D^{\top}(q)=0_{(n-m)\times(n-m)}\end{split}$ (57)

for all $i\in\left\{1,\dots,n\right\}$ where

$D(q)=\begin{bmatrix}(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}&-I_{n-m}\end{bmatrix}.$ (58)

•

The desired potential energy $V_{d}(q)$ satisfies

\begin{split}s_{1}(q)G^{\perp}\nabla_{q}V&=-s_{2}(q)G^{\top}\nabla_{q}V_{d}-s_{3}(q)G^{\perp}\nabla_{q}V_{d}\\ &=-D(q)M^{-1}(q)\nabla_{q}V_{d},\end{split}

(59)

where

\begin{split}s_{1}&(q)=(m_{22}+m_{a22})\\ &-(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}(m_{21}^{\top}+m_{a21}^{\top})\end{split}

(60)

is the Schur complement of $M^{-1}(q)+M_{a}^{-1}(q)$ and

\begin{split}s_{2}(q)&=(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}m_{11}-m_{21}\\ s_{3}(q)&=(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}m_{21}^{\top}-m_{22}.\end{split}

(61)

Proof:

From the graph in Figure 1 and the interconnection (37), the virtual input $u_{v}$ is given by

\begin{split}u_{v}=Gu=&Gv-C+\nabla_{q}V.\end{split}

(62)

Recalling the definition of $G$ , (62) is equivalent to (42). Collecting the terms $u,v$ and left multiplying by $M_{a}^{-1}+M^{-1}$ results in

\begin{split}&\left[M_{a}^{-1}+M^{-1}\right]G(u-v)\\ &\phantom{--}=\left[M_{a}^{-1}+M^{-1}\right]\left\{-C_{KE}-C_{PE}+\nabla_{q}V\right\}\end{split}

(63)

Due to the structure of $G$ in (8), (63) has the left annihilator $D(\cdot)$ , defined in (58).

Left multiplying (63) by $D(\cdot)$ and separating the components into those relating to the kinetic and potential energies result in

	$\displaystyle 0_{(n-m)\times 1}$	$\displaystyle=-D\left[M_{a}^{-1}+M^{-1}\right]C_{KE}$		(64)
	$\displaystyle 0_{(n-m)\times 1}$	$\displaystyle=-D\left[M_{a}^{-1}+M^{-1}\right]\left\{C_{PE}-\nabla_{q}V\right\}.$		(65)

Using (52), (65) is expanded to

\begin{split}0_{(n-m)\times 1}=&D\left[M_{a}^{-1}+M^{-1}\right]\nabla_{q}V\\ &-DM^{-1}\nabla_{q}V_{a},\end{split}

(66)

which can be seen to agree with (59) after expanding.

Now considering the constraint on the kinetic energy expression (64), the definition (52) is substituted to find

\begin{split}0_{(n-m)\times n}=&D\left\{-Y+J\right\}.\\ \end{split}

(67)

Using the relevant definitions, the first component of (67) can be solved for $J_{21}(q,p)$ as

\begin{split}J_{21}(q,p)=&(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}(-Y_{11}+J_{11})\\ &+Y_{21}.\end{split}

(68)

Substituting this expression back into the second component of (67) reveals the constraint

\begin{split}0&{}_{(n-m)\times(n-m)}\\ =&(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}(-Y_{12}-J_{21}^{\top})\\ &\phantom{=}-(-Y_{22}+J_{22})\\ =&(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}\left[-Y_{12}-Y_{21}^{\top}\right]\\ &-(m_{21}+m_{a21})(m_{11}+m_{a11})^{-1}(-Y_{11}+J_{11})^{\top}\\ &\times(m_{11}+m_{a11})^{-1}(m_{21}^{\top}+m_{a21}^{\top})-(-Y_{22}+J_{22})\\ =&-D\begin{bmatrix}-Y_{11}+J_{11}^{\top}&-Y_{12}-Y_{21}^{\top}\\ 0_{(n-m)\times m}&-Y_{22}+J_{22}\end{bmatrix}D^{\top}.\end{split}

(69)

The term $J_{22}$ is taken as below to solve the skew-symmetric part of this expression,

\begin{split}&J_{22}\\ &=-\frac{1}{2}D\begin{bmatrix}-Y_{11}+Y_{11}^{\top}+J_{11}^{\top}-J_{11}&-Y_{12}-Y_{21}^{\top}\\ Y_{12}^{\top}+Y_{21}&-Y_{22}+Y_{22}^{\top}\end{bmatrix}D^{\top},\end{split}

(70)

where $J_{11}\in\mathbb{R}^{m\times m}$ is a free skew-symmetric term. The symmetric part of (69) must also be equal to zero, implying that

\begin{split}D\left[Y+Y^{\top}\right]D^{\top}=0_{(n-m)\times(n-m)}.\end{split}

(71)

Finally, noting that this must be true for each $p_{i}$ , the condition (57) follows. ∎

Remark 1

The expression (57) implicitly defines a set of PDEs that must be satisfied by any choice of $M_{a}^{-1}(q)$ . From the definition of $Y^{i}$ in (55), the first $m$ equations are describe partial differential equations involving the partial derivatives of $m_{a11},m_{a21}$ . The remaining $n-m$ equations describe partial differential equations involving the partial derivatives of $m_{a21},m_{a22}$ . This structure can be useful for resolving the equations into a standard representation for solving.

Corollary 2

In the special case of under-actuation degree 1, if $M^{-1}$ , $M_{a}^{-1}$ is a function of only 1 configuration variable $q_{i}$ , the kinetic energy matching equations (57) can be reduced to a set of ODEs.

\begin{split}\frac{d}{dq_{i}}\begin{bmatrix}m_{a21}^{\top}\\ m_{a22}\end{bmatrix}=g\left(m_{a11},\frac{d}{dq_{i}}m_{a11},M^{-1},\frac{d}{dq_{i}}M^{-1}\right),\end{split}

(72)

where $g(\cdot)\in\mathbb{R}^{n}$ is a function implicitly defined by the matching conditions (57) and $m_{a11}(q_{i})$ can be chosen freely.

Proof:

Assuming that the mass matrix $M^{-1}$ is function only of a single configuration variable $q_{i}$ , we will also impose that the added mass $M_{a}^{-1}$ is a function only of the same variable. As a consequence, the matching expression (57) is now only a function in the single variable $q_{i}$ . Notably, all partial derivatives of $M_{a}^{-1}$ with respect to $q_{k}$ , where $k\neq i$ , are equal to zero.

Noting Remark 1, the first $n-1$ expressions of (57) produce differential equations involving the partial derivatives of $m_{a11},m_{a21}$ . The dimension of $m_{a21}$ is $1\times(n-1)$ , so the first $n-1$ equations can be solved simultaneously to find an expression for $\frac{d}{dq_{i}}m_{a21}^{\top}$ . The $n^{th}$ expressions of (57) can then be resolved for an expression for $\frac{d}{dq_{i}}m_{a22}$ , which has dimension 1. Combining these expressions, the matching equations (57) can be resolved into an ODE of the form (72). ∎

Remark 2

Corollary 2 describes situations in which the kinetic energy matching equations can be reduced to an ODE. The solution, however, will depend on the choice of $m_{a11}(q_{i})$ and may not be globally defined. This poses the question of how should the function $m_{a11}(q_{i})$ be chosen to ensure an appropriate solution $M_{a}^{-1}(q_{i})$ —a nonlinear control problem!

The results of Corollary 1 describe the degrees of freedom that exist when constructing a solution to the added inverse mass matrix. Similar degrees of freedom exist in the definition of the closed-loop potential energy that can be exploited to ensure positivity of the chosen function. The following Corollary defines a free function that can be utilised to this effect.

Corollary 3

Suppose that there exists a full rank matrix-valued function $K(q)\in\mathbb{R}^{m\times m}$ such that the integral

\begin{split}\Gamma(q)&=\int K(q)G^{\top}\left[M_{a}^{-1}(q)+M^{-1}(q)\right]M(q)\ dq,\end{split}

(73)

exists. The desired closed-loop potential energy can be chosen as

V_{d}(q)=V_{m}(q)+V_{f}(\Gamma(q)),

(74)

where $V_{m}(\cdot)$ must be chosen to satisfy the potential energy matching conditions (59) and $V_{f}(\cdot)$ is a free function that does not impact the matching equations. Consequently, the matching equation (59) can be equivalently written as

\begin{split}s_{1}(q)G^{\perp}\nabla_{q}V&=-s_{2}(q)G^{\top}\nabla_{q}V_{m}-s_{3}(q)G^{\perp}\nabla_{q}V_{m}\\ &=-D(q)M^{-1}(q)\nabla_{q}V_{m}.\end{split}

(75)

Proof:

Computing the gradient of $V_{d}$ results in

\begin{split}\nabla_{q}V_{d}&=\nabla_{q}V_{m}+\frac{\partial^{\top}\Gamma}{\partial q}\nabla_{\Gamma}V_{f}\\ &=\nabla_{q}V_{m}+M\left[M_{a}^{-1}+M^{-1}\right]GK^{\top}\nabla_{\Gamma}V_{f}.\end{split}

(76)

From the definition of $D(q)$ in (58), we have the identity

\begin{split}D(q)M^{-1}M(q)\left[M_{a}^{-1}(q)+M^{-1}(q)\right]G&=\begin{bmatrix}I_{m}&\star\end{bmatrix}G\\ &=0_{(n-m)\times m}\end{split}

(77)

Substituting the expression (76) into (59) and noting the above expression results in the simplified matching equation (75). ∎

Remark 3

$V_{f}(\cdot)$ is a free function precisely because $\Gamma$ is an integral of the passive output $y_{c2}$ . The potential energy $V_{f}$ could be alternatively constructed as a capacitor element added to the input $u_{c2}$ in Figure 1.

Now we arrive at one of the key results of this work, the equivalence of the proposed CbI scheme and total energy-shaping control of underactuated mechanical systems. Assuming that the CbI scheme has been constructed to satisfy the required matching conditions to ensure the existence of a Casimir of the form $q_{a}=q,p_{a}=p$ , Proposition 1 is applied to reconstruct the reduced closed-loop structure (18).

Proposition 4

Consider the underactuated mechanical system with virtual input (14) and assume that $M_{a}(q),V_{d}(q)$ are chosen such that the conditions of Proposition 3 are satisfied in some neighbourhood of a point $(q,p)=(q^{\star},0_{n\times 1})$ . If the control signal is chosen as

\begin{split}u(q,p)=&v-G^{\top}\left\{M_{d}(q)M^{-1}(q)\left[-E^{\top}(q,p)M_{a}^{-1}(q)p\right.\right.\\ &\left.\left.+\nabla_{q_{a}}H_{a}(q_{a},p_{a})-M(q)J(q,p)p\right]-\nabla_{q}V\right\}\end{split}

(78)

where

M_{d}(q)=\left[M_{a}^{-1}(q)+M^{-1}(q)\right]^{-1},

(79)

the following hold:

•

The closed-loop dynamics have the form

\begin{split}\begin{bmatrix}\dot{q}\\ \dot{p}\end{bmatrix}=&\begin{bmatrix}0_{n\times n}&M^{-1}(q)M_{d}(q)\\ -M_{d}(q)M^{-1}(q)&J_{2}(q,p)\end{bmatrix}\begin{bmatrix}\nabla_{q}H_{d}\\ \nabla_{p}H_{d}\end{bmatrix}\\ &+\begin{bmatrix}0_{n\times m}\\ G\end{bmatrix}v\\ H_{d}(q,p&)=\frac{1}{2}p^{\top}M_{d}^{-1}(q)p+V_{d}(q)\\ &y=G^{\top}\nabla_{p}H_{d},\end{split}

(80)

where

\begin{split}J_{2}(q,p)=&M_{d}(q)\left\{J(q,p)+M^{-1}(q)\left[E(q,p)-E^{\top}(q,p)\right]\right.\\ &\left.\times M^{-1}(q)\right\}M_{d}(q)+M_{d}(q)M^{-1}(q)E^{\top}(q,p)\\ &-E(q,p)M^{-1}(q)M_{d}(q)\end{split}

(81)

•

If $M_{d}(q),V_{d}(q)$ satisfy

$\begin{split}M_{d}(q)>0,\ \ V_{d}(q)>0\end{split}$ (82)

in some neighbourhood of $(q,p)=(q^{\star},0_{n\times 1})$ , $(q^{\star},0_{n\times 1})$ is a stable equilibrium of the closed-loop system for $v=0_{m\times 1}$ .
•

If the input signal $v$ is used for damping injection

$v=-K_{d}G^{\top}y$ (83)

for some positive $K_{d}\in\mathbb{R}^{m\times m}$ and the equilibrium $(q,p)=(q^{\star},0_{n\times 1})$ , $(q^{\star},0_{n\times 1})$ is locally detectable from the output $y$ , the point $(q^{\star},0_{n\times 1})$ is asymptotically stable.

Proof:

Interconnection of the mechanical system with the control subsystem results in a closed-loop of the form (4), where $x_{c}$ , $H_{c}$ and $K_{\star\star}$ are defined in (50) and

\begin{split}x_{p}=&\begin{bmatrix}q\\ p\end{bmatrix}\\ F_{p}=&\begin{bmatrix}0_{n\times n}&I_{n}\\ -I_{n}&0_{n\times n}\end{bmatrix}\\ G_{p}=&\begin{bmatrix}0_{n\times n}\\ I_{n}\end{bmatrix}.\end{split}

(84)

From (38), we have that

\frac{\partial f_{c}}{\partial x_{p}}=\begin{bmatrix}I_{n}&0_{n\times n}\\ I_{n}&0_{n\times n}\\ 0_{n\times n}&I_{n}\end{bmatrix}.

(85)

To verify the claim, Corollary 1 is applied which requires a suitable definition of $B$ . Expanding the definitions of $\bar{F}_{\star 3}$ from (34) reveals

\begin{split}\bar{F}_{13}&=\begin{bmatrix}0_{n\times n}&0_{n\times n}&-I_{n}\\ 0_{n\times n}&I_{n}-M_{d}M^{-1}&D\end{bmatrix}\\ \bar{F}_{23}&=\begin{bmatrix}0_{n\times n}&0_{n\times n}&0_{n\times n}\end{bmatrix}\\ \bar{F}_{33}&=\begin{bmatrix}0_{n\times n}&0_{n\times n}&0_{n\times n}\\ 0_{n\times n}&0_{n\times n}&I_{n}\\ 0_{n\times n}&-I_{n}&0_{n\times n}\\ \end{bmatrix},\end{split}

(86)

resulting in the choice

B=\begin{bmatrix}0_{n\times n}&0_{n\times n}\\ I_{n}&0_{n\times n}\\ 0_{n\times n}&I_{n}\\ \end{bmatrix}.

(87)

Expanding the expression $B^{\top}\bar{F}_{33}B$ results in

B^{\top}\bar{F}_{33}B=\begin{bmatrix}0_{n\times n}&-I_{n}\\ I_{n}&0_{n\times n}\end{bmatrix}

(88)

which is invertible, ensuring that Corollary 1 can be applied. Expanding the definitions of $F_{r}$ in (32) results in the reduced dynamics

\begin{split}\begin{bmatrix}\dot{q}\\ \dot{p}\\ -y\end{bmatrix}&=\underbrace{\begin{bmatrix}0_{n\times n}&M^{-1}M_{d}&0_{n\times n}\\ -M_{d}M^{-1}&\bar{J}_{2}&G\\ 0_{n\times n}&-G^{\top}&0_{n\times n}\end{bmatrix}}_{F_{r}}\begin{bmatrix}\nabla_{q}H_{d}\\ \nabla_{p}H_{d}\\ v\end{bmatrix}\\ \bar{J}_{2}&=M_{d}JM_{d}+D-D^{\top}+M_{d}M^{-1}D^{\top}-DM^{-1}M_{d},\end{split}

(89)

which agrees with (80) when substituting in the definition for $D$ in (50). Stability and asymptotic stability of the point $(q^{\star},0_{n\times 1})$ follows from Proposition 1 of [17]. ∎

Remark 4

From Proposition 4 it is clear that $M_{a}^{-1}(q)$ does not need to be a positive matrix. Rather, the closed-loop mass $M_{d}^{-1}$ must be positive to ensure stability fo the system. In cases that $M_{a}^{-1}$ is positive, the control sub-system in Figure 1 is passive.

V Example applications

In this section the matching conditions derived in Proposition 3 are used to construct stabilising control laws for the cart-pole and acrobot systems. In both cases, the mass matrix depends on only one configuration variable, so the kinetic energy matching conditions can be reduced to ODEs as detailed in Corollary 2. This enables the solutions to be constructed numerically, removing the need to analytically solve the equations.

Both examples were prepared in Matlab 2022a and the source code is available via https://github.com/JoelFerguson/Underactuated_Mechanical_CbI.

V-A Cart-pole example

The cart-pole system, shown in Figure 2, attempts to balance the pole of length $\ell$ and mass $m_{p}$ in the upright position by applying a force $F$ to the cart with mass $m_{c}$ . The state $q_{1}$ describes the horizontal displacement of the cart whereas $q_{2}$ describes the angle of the pole from vertical in the clockwise direction.

The cart-pole system can be written as a pH system of the form (7) with

\begin{split}q&=\begin{bmatrix}q_{1}\\ q_{2}\end{bmatrix}\\ M(q)&=\begin{bmatrix}m_{c}+m_{p}&m_{p}l\cos q_{2}\\ m_{p}l\cos q_{2}&m_{p}l^{2}\end{bmatrix}\\ V(q)&=m_{p}gl\cos q_{2}\\ G&=\begin{bmatrix}1\\ 0\end{bmatrix}.\end{split}

(90)

In the subsequent control design, the parameters $m_{c}=m_{p}=l=1,g=9.8$ have been used.

The mass matrix of the cart-pole system depends only on $q_{2}$ , the unactuated coordinate. The added inverse mass is assumed to also be a function of $q_{2}$ also, allowing it to be written as

M_{a}^{-1}(q_{2})=\begin{bmatrix}m_{a11}(q_{2})&m_{a21}^{\top}(q_{2})\\ m_{a21}(q_{2})&m_{a22}(q_{2})\end{bmatrix}.

(91)

As noted in Corollary 2, the kinetic energy matching equations (57) can be reduced to an ODE as both $M^{-1},M_{a}^{-1}$ are a function of only one variable. The associate ODE is of the form (72) for $q_{i}=q_{2}$ where $m_{a11}(q_{2})$ is a free function to be chosen. The ODE can be evaluated using numerical solvers.

Before solving the ODE associated with the kinetic energy matching equations, consideration should be given to how the resulting mass matrix impacts the closed-loop potential energy $V_{d}$ . Recalling (74), the closed-loop potential energy is composed of a free term $\Gamma(\cdot)$ and a term $V_{m}(q)$ which must satisfy (75), where $s_{1},s_{2},s_{3}$ are defined in (60), (61). As the potential $V,M^{-1},M_{a}^{-1}$ are all functions of only $q_{2}$ , $V_{m}$ is also assumed to be a function of $q_{2}$ only, reducing (75) to the ODE

\begin{split}\nabla_{q_{2}}V_{m}&=-\frac{s_{1}(q_{2})}{s_{3}(q_{2})}\nabla_{q_{2}}V,\end{split}

(92)

which can be evaluated numerically once a solution for $M_{a}^{-1}(q_{2})$ , and hence $s_{1}(\cdot),s_{3}(\cdot)$ , are found. noting that the vector field $\nabla_{q_{2}}V$ is divergent from the point $q_{2}=0$ , the closed-loop vector field $\nabla_{q_{2}}V_{m}$ should reverse the direction locally. This is ensured if the ratio $\frac{s_{1}(q)}{s_{3}(q)}$ is positive in some neighbourhood of the origin. Recalling that $s_{1}(q)$ is the Schur complement of $M^{-1}+M_{a}^{-1}$ , which is necessarily positive, it is required that $s_{3}(q)$ be positive in some neighbourhood of $q_{2}=0$ . The values

\begin{split}m_{a11}(0)&=0\\ m_{a21}(0)&=-2\\ m_{a22}(0)&=8,\end{split}

(93)

where chosen which result in $s_{1}(0)=1$ , $s_{3}(0)=1$ and $\lambda_{min}\left[M^{-1}(0)+M_{a}^{-1}(0)\right]=0.917>0$ .

The added inverse mass matrix can now be found by numerically evaluating the ODE (72). The term $m_{a11}(q_{2})$ is a free function that was chosen to be constant $m_{a11}(q_{2})=0,\frac{\partial}{\partial q_{2}}m_{a11}=0$ for this example. The resulting functions for $m_{a21}(q_{2}),m_{a22}(q_{2})$ were found to exist on the interval $q_{2}\in\left[-0.48,0.48\right]$ and are shown in Figure (3). From Proposition 4, $M^{-1}(q_{2})+M_{a}^{-1}(q_{2})$ should be positive to ensure stability, so the minimum eigenvalue of this expression is shown in the same figure.

The closed-loop potential energy $V_{m}(q_{2})$ can now be obtained by numerically by evaluating the ODE (92). The terms $s_{1}(\cdot),s_{3}(\cdot)$ are evaluated using the solutions to $M_{a}^{-1}$ shown in Figure (3). The resulting function $V_{m}(q_{2})$ is shown in Figure (4). As expected, the function is positive in some neighbourhood of $q_{2}=0$ due to the choice of the added mass at $q_{2}=0$ in (93).

The proposed functions of $M_{a}^{-1}$ , $V_{m}$ can be used to construct a controller to stabilise the pendulum in the upright position. To ensure stability of $q_{1}=0$ also, the free term $V_{f}(\Gamma(q))$ , defined in (74), is constructed. The function $\Gamma(\cdot)$ defined by the integral (73), where $K(q)$ is a free function chosen to ensure solvability. Noting that $M^{-1},M_{a}^{-1}$ are functions of $q_{2}$ only, the parametrisation

\begin{split}\begin{bmatrix}\beta_{1}(q_{2})&\beta_{2}(q_{2})\end{bmatrix}=G^{\top}\left[M_{a}^{-1}(q_{2})+M^{-1}(q_{2})\right]M(q_{2})\end{split}

(94)

is introduced. The free function is chosen as $K(q_{2})=\frac{1}{\beta_{1}(q_{2})}$ , resulting in

\begin{split}\Gamma(q)&=\int\begin{bmatrix}1&\frac{\beta_{2}(q_{2})}{\beta_{1}(q_{2})}\end{bmatrix}\ dq\\ &=q_{1}+\int\frac{\beta_{2}(q_{2})}{\beta_{1}(q_{2})}\ dq_{2},\end{split}

(95)

which can be solved numerically from the initial condition $\Gamma(0_{2\times 1})=0$ . The function $V_{f}(\cdot)$ was taken as $V_{f}(\Gamma(q))=\frac{1}{2}\kappa\Gamma(q)^{2}$ with $\kappa=5$ for simulation. A contour plot of the resulting closed-loop potential energy is shown in Figure (5). Note that a minimum has been assigned to $q=0_{2\times 1}$ .

As a final control design stage, damping is injected via the new passive input/output pair with

v=-5G^{\top}(M_{a}^{-1}+M^{-1})p.

(96)

The complete control signal is defined by the expression (46).

The cart-pole system was simulated for 5 seconds from initial conditions $q(0)=(0,0.3)$ , $p(0)=(0,0)$ . The resulting state evolution and closed-loop energy $H_{d}$ is shown in Figure 6. As expected, the proposed controller stabilises the origin and the closed-loop energy $H_{d}$ decreases monotonically.

V-B Acrobot example

The acrobot system, shown in Figure 7, consists of 2 links with an actuator supplying a input torque $\tau$ fixed between the base and second links. The base link has displacement of $q_{2}$ , measured from vertical, length $\ell_{2}$ , mass $m_{2}$ , moment of inertia $J_{\ell 1}$ and centre of mass $\ell_{c2}$ from the base pivot point. The actuated link has displacement of $q_{1}$ measured relative to the base link, length $\ell_{1}$ , mass $m_{1}$ , moment of inertia $J_{\ell 1}$ and centre of mass $\ell_{c1}$ from the actuated pivot point. The control objective of this system is to stabilise the upright equilibrium position $(q_{1},q_{2})=(0,0)$ .

The acrobot system can be written as a pH system of the form (7) with

\begin{split}M(q)&=\begin{bmatrix}c_{2}&c_{2}+c_{3}\cos q_{1}\\ c_{2}+c_{3}\cos q_{1}&c_{1}+c_{2}+2c_{3}\cos q_{1}\end{bmatrix}\\ V(q)&=c_{4}g\cos q_{2}+c_{5}g\cos(q_{1}+q_{2})\\ G&=\begin{bmatrix}1\\ 0\end{bmatrix},\end{split}

(97)

where

\begin{split}c_{1}&=m_{2}\ell_{c2}^{2}+m_{1}\ell_{2}^{2}+J_{\ell 2}\\ c_{2}&=m_{1}\ell_{c1}^{2}+J_{\ell 1}\\ c_{3}&=m_{1}\ell_{2}\ell_{c1}\\ c_{4}&=m_{2}\ell_{c2}+m_{1}\ell_{1}\\ c_{5}&=m_{1}\ell_{c1}.\end{split}

(98)

For the purposes of simulation, we take the values $g=9.8$ , $c_{1}=2.3333,c_{2}=5.3333,c_{3}=2,c_{4}=3,c_{5}=2$ which were previously used in [5], [25].

In this example, the total energy-shaping controller proposed in [5] is reconstructed as a CbI control scheme by solving the matching conditions of Proposition 3. In that work, the closed-loop mass matrix was chosen to be the constant matrix

\begin{split}M_{d}^{-1}&=\begin{bmatrix}0.3385&-0.9997\\ -0.9997&5.9058\\ \end{bmatrix}\end{split}

(99)

which will be recovered in subsequent computations.

The mass matrix of the acrobot system depends only on $q_{1}$ , the actuated coordinate. The added inverse mass matrix is assumed to be a function of only $q_{1}$ also, resulting in the structure

M_{a}^{-1}(q_{1})=\begin{bmatrix}m_{a11}(q_{1})&m_{a21}^{\top}(q_{1})\\ m_{a21}(q_{1})&m_{a22}(q_{1})\end{bmatrix}.

(100)

As the system is underactuated degree 1 and the mass matrix is a function of only one variable, the kinetic energy matching equations can be reduced to an ODE as per Corollary 2. The resulting ODE has the form (72) with $q_{i}=q_{1}$ and where $m_{a11}(q_{1})$ is a free function.

In order to recover the result (99), this free function $m_{a11}(q_{1})$ is chosen as

\begin{split}m_{a11}(q_{1})&=G^{\top}\left[M_{d}^{-1}-M^{-1}(q_{1})\right]G\\ &=0.3385-\frac{c_{1}+c_{2}+2c_{3}\cos q_{1}}{c_{1}c_{2}-c_{3}^{2}\cos^{2}(q_{1})}.\end{split}

(101)

The initial conditions $m_{a12}(0),m_{a22}(0)$ are similarly defined as

\begin{split}m_{a12}(0)&=G^{\perp}\left[M_{d}^{-1}-M^{-1}(0)\right]G=-0.1313\\ m_{a12}(0)&=G^{\perp}\left[M_{d}^{-1}-M^{-1}(0)\right]G^{\perp\top}=5.2743.\end{split}

(102)

The added inverse mass was evaluated numerically and the results are shown in Figure 8. As the previously reported solution (99) is globally defined, it is unsurprising that the inverse added mas is also globally defined. As expected, the minimum eigenvalue of $M^{-1}+M_{a}^{-1}$ is constant also.

Solving the potential energy PDE (75) is difficult due to the open-loop potential energy being a function of both $q_{1}$ and $q_{2}$ . This dependence implies that $V_{m}$ cannot be resolved directly using an ODE solver. Considering the structure of $V$ in (97), it is proposed that the closed-loop energy $V_{m}$ has the structure

V_{m}(q)=f_{1}(q_{1})\sin(q_{2})+f_{2}(q_{1})\cos(q_{2}),

(103)

which has derivatives

\begin{split}\nabla_{q_{1}}V_{a}=&\frac{\partial f_{1}}{\partial q_{1}}\sin(q_{2})+\frac{\partial f_{2}}{\partial q_{1}}\cos(q_{2})\\ \nabla_{q_{2}}V_{a}=&f_{1}(q_{1})\cos(q_{2})-f_{2}(q_{1})\sin(q_{2}).\end{split}

(104)

The open-loop potential energy has gradients

\begin{split}\nabla_{q_{1}}V=&-c_{5}g\sin(q_{1})\cos(q_{2})-c_{5}g\cos(q_{1})\sin(q_{2})\\ \nabla_{q_{2}}V=&-c_{4}g\sin(q_{2})-c_{5}g\sin(q_{1})\cos(q_{2})\\ &-c_{5}g\cos(q_{1})\sin(q_{2}).\end{split}

(105)

Substituting the expressions (104) and (105) into (59) and matching coefficients results in the system of equations

\begin{split}\begin{bmatrix}\frac{\partial f_{1}}{\partial q_{1}}\\ \frac{\partial f_{2}}{\partial q_{1}}\end{bmatrix}&=\frac{1}{s_{2}(q_{1})}\\ &\times\begin{bmatrix}c_{4}gs_{1}(q_{1})+c_{5}gs_{1}(q_{1})\cos(q_{1})+s_{3}(q_{1})f_{2}(q_{1})\\ c_{5}gs_{1}(q_{1})\sin(q_{1})-s_{3}(q_{1})f_{1}(q_{1})\end{bmatrix},\end{split}

(106)

which can be evaluated numerically. The values of $f_{1},f_{2}$ at the origin should be chosen to ensure that the origin is an equilibrium point and $V_{m}$ is positive in $q_{2}$ . Considering the expressions (105), (106), the origin is an equilibrium for $f_{1}(0)=0$ . The energy function (103) is locally positive with respect to $q_{1}$ for $f_{2}(0)$ negative. For the purpose of simulation, $f_{2}(0)=-50$ was used. The resulting function $V_{m}$ is shown in Figure 9.

Considering Figure 9, it is clear that $V_{m}$ is not positive definite with respect to the origin. Note, however, that $q_{2}=0$ has been stabilised. To ensure stability of $q_{1}=0$ also, the free term $V_{f}(\Gamma(q))$ , defined in (74), is constructed. The function $\Gamma(\cdot)$ defined by the integral (73), where $K(q)$ is a free function chosen to ensure solvability. Noting that $M^{-1},M_{a}^{-1}$ are functions of $q_{1}$ only, the parametrisation

\begin{split}\begin{bmatrix}\beta_{1}(q_{1})&\beta_{2}(q_{1})\end{bmatrix}=G^{\top}\left[M_{a}^{-1}(q_{1})+M^{-1}(q_{1})\right]M(q_{1})\end{split}

(107)

is introduced. The free function is chosen as $K(q_{1})=\frac{1}{\beta_{2}(q_{1})}$ , resulting in

\begin{split}\Gamma(q)&=\int\begin{bmatrix}\frac{\beta_{1}(q_{1})}{\beta_{2}(q_{1})}&1\end{bmatrix}\ dq\\ &=\int\frac{\beta_{1}(q_{1})}{\beta_{2}(q_{1})}\ dq_{1}+q_{2},\end{split}

(108)

which can be solved numerically from the initial condition $\Gamma(0_{2\times 1})=0$ . The function $V_{f}(\cdot)$ was taken as $V_{f}(\Gamma(q))=\frac{1}{2}\kappa\Gamma(q)^{2}$ with $\kappa=250$ for simulation. A contour plot of the resulting closed-loop potential energy on a log scale is shown in Figure 10. Note that a minimum has been assigned to $q=0_{2\times 1}$ . As a final control design stage, damping is injected via the new passive input/output pair with

v=-5G^{\top}(M_{a}^{-1}+M^{-1})p.

(109)

The complete control signal is defined by the expression (46).

The acrobot system was simulated for 20 seconds from initial conditions $q(0)=(0,0.5)$ , $p(0)=(0,0)$ . The resulting state evolution and closed-loop energy $H_{d}$ is shown in Figure 11. As expected, the proposed controller stabilises the origin and the closed-loop energy $H_{d}$ decreases monotonically.

VI CONCLUSIONS AND FUTURE WORKS

In this work total energy shaping has been shown to have a CbI interpretation which results in alternate matching equations related to the added inverse mass. These equations were utilised to construct controllers for the cart-pole and acrobot systems, both of which have the property that the mass matrix depends on only one variable, using numerical methods. While the proposed approach is effective, a number of technical aspects of this approach require further investigation. In particular:

•

As detailed in Corollary 2, The kinetic energy matching equations can be posed as ODEs in the special case that the mass matrix depends on only one configuration variable. This property allows the matching equations to be evaluated numerical using ODE solvers. Further investigation into solving the matching equations in the case that the mass matrix is a function of multiple configuration variables is required. In some cases it may be possible to decouple the dependence on each coordinate, recovering equivalent ODEs. Alternatively, the numerical evaluation of the matching PDEs should be investigated.
•

When evaluating the kinetic energy matching equations in (72), the term $m_{a11}(q_{i})$ is a free function that can be used to control the resulting added inverse mass. As seen in the cart-pole example of Section V-A, poor choice of this function results in the solution only being defined on a small domain. Conversely, in the acrobot example of Section V-A this term was chosen to ensure a global solution to the matching equations. Choice of this function defines a nonlinear control problem that should be investigated to ensure desirable behaviour of the result.
•

In both examples of Section V the controllers were designed to stabilise the origin of the respective systems. While this was achieved and verified numerically, asymptotic stability was not established. Asymptotic stability requires that the passive output of the closed-loop system is zero-state detectable, a task that is non-trivial for underactuated systems. Further investigation into methods for injecting damping into the unactuated momentum channels of the closed-loop system is required. It is hoped that the CbI interpretation of the controller shown in Figure 1 may provide new insight into how this might be achieved.

References

[1] R. Ortega and E. Garcia-Canseco, “Interconnection and damping assignment passivity-based control: A survey,” European Journal of control, vol. 10, no. 5, pp. 432–450, 2004.
[2] R. Ortega, A. Van der Schaft, B. Maschke, and G. Escobar, “Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems,” Automatica, vol. 38, no. 4, pp. 585–596, 2002.
[3] F. Gómez-Estern, R. Ortega, F. R. Rubio, and J. Aracil, “Stabilization of a class of underactuated mechanical systems via total energy shaping,” Proceedings of the IEEE Conference on Decision and Control, vol. 2, no. December, pp. 1137–1143, 2001.
[4] J. Acosta, R. Ortega, and A. Astolfi, “Interconnection and damping assignment passivity-based control of mechanical systems with underactuation degree one,” IEEE Transactions on Automatic Control, vol. 50, no. 12, pp. 1936–1955, 2005.
[5] A. D. Mahindrakar, A. Astolf, R. Ortega, and G. Viola, “Further constructive results on interconnection and damping assignment control of mechanical systems: The Acrobot example,” International Journal of Robust and Nonlinear Control, vol. 18, no. July, pp. 557–569, 2010.
[6] P. Arpenti, F. Ruggiero, and V. Lippiello, “A Constructive Methodology for the IDA-PBC of Underactuated 2-DoF Mechanical Systems with Explicit Solution of PDEs,” International Journal of Control, Automation and Systems, vol. 20, no. 1, pp. 283–297, 2022.
[7] R. Ortega, A. van der Schaft, F. Castaños, and A. Astolfi, “Control by interconnection and standard passivity-based control of port-Hamiltonian systems,” IEEE Transactions on Automatic Control, vol. 53, no. 11, pp. 2527–2542, 2008.
[8] R. Ortega and L. P. Borja, “New results on Control by Interconnection and Energy-Balancing Passivity-Based Control of port-hamiltonian systems,” Proceedings of the IEEE Conference on Decision and Control, vol. 2015-Febru, no. February, pp. 2346–2351, 2014.
[9] V. Duindam, A. Macchelli, S. Stramigioli, and H. Bruyninckx, Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach. Berlin Heidelberg: Springer-Verlag, 2009.
[10] A. Donaire, R. Mehra, R. Ortega, S. Satpute, J. G. Romero, F. Kazi, and N. M. Singh, “Shaping the Energy of Mechanical Systems Without Solving Partial Differential Equations,” IEEE Transactions on Automatic Control, vol. 61, no. 4, pp. 1051–1056, 2016.
[11] J. G. Romero, A. Donaire, and R. Ortega, “Global Stabilisation of Underactuated Mechanical Systems via PID Passivity-Based Control,” pp. 1–27, 2016.
[12] J. G. Romero, A. Donaire, R. Ortega, and P. Borja, “Global Stabilisation of Underactuated Mechanical Systems via PID Passivity-Based Control,” IFAC-PapersOnLine, vol. 50, no. 1, pp. 9577–9582, 2017.
[13] M. Zhang, P. Borja, R. Ortega, Z. Liu, and H. Su, “PID Passivity-Based Control of Port-Hamiltonian Systems,” IEEE Transactions on Automatic Control, vol. PP, no. 99, pp. 1–1, 2017.
[14] P. J. Gawthrop and G. P. Bevan, “Bond-Graph Modeling: A tutorial introduction for control engineers,” IEEE Control Systems Magazine, vol. 27, no. 2, pp. 24–45, 2007.
[15] A. van der Schaft, L2-Gain and Passivity Techniques in Nonlinear Control, 3rd ed. Springer, 2017.
[16] R. Reyes-Báez, A. van der Schaft, and B. Jayawardhana, “Tracking Control of Fully-actuated port-Hamiltonian Mechanical Systems via Sliding Manifolds and Contraction Analysis,” in Proc. IFAC World Congress. Toulouse, France: Elsevier, 2017, pp. 8256–8261.
[17] R. Ortega, M. W. Spong, F. Gómez-Estern, and G. Blankenstein, “Stabilization of a class of underactuated mechanical systems via interconnection and damping assignment,” IEEE Transactions on Automatic Control, vol. 47, no. 8, pp. 1218–1233, 2002.
[18] G. Viola, R. Ortega, and R. Banavar, “Total energy shaping control of mechanical systems: simplifying the matching equations via coordinate changes,” IEEE Transactions on Automatic Control, vol. 52, no. 6, pp. 1093–1099, 2007.
[19] M. Ryalat and D. S. Laila, “A simplified IDA-PBC design for underactuated mechanical systems with applications,” European Journal of Control, vol. 27, pp. 1–16, 2016.
[20] F. Gómez-Estern and A. van der Schaft, “Physical damping in IDA-PBC controlled underactuated mechanical Systems,” European Journal of Control, vol. 10, no. 5, pp. 451–468, 2004.
[21] J. Sandoval, R. Kelly, and V. Santibanez, “Interconnection and damping assignment passivity-based control of a class of underactuated mechanical systems with dynamic friction,” International Journal of Robust and Nonlinear Control, vol. 21, no. 7, pp. 738–751, 2010.
[22] A. Donaire, R. Ortega, and J. G. Romero, “Simultaneous interconnection and damping assignment passivity-based control of mechanical systems using dissipative forces,” Systems & Control Letters, vol. 94, pp. 118–126, 2016.
[23] O. B. Cieza and J. Reger, “IDA-PBC for underactuated mechanical systems in implicit port-hamiltonian representation,” 2019 18th European Control Conference, ECC 2019, pp. 614–619, 2019.
[24] M. R. J. Harandi and H. D. Taghirad, “Solution of matching equations of IDA-PBC by Pfaffian differential equations,” International Journal of Control, pp. 1–11, 2021.
[25] A. Donaire, J. G. Romero, R. Ortega, B. Siciliano, and M. Crespo, “Robust IDA-PBC for underactuated mechanical systems subject to matched disturbances,” International Journal of Robust and Nonlinear Control, vol. 27, no. 6, pp. 1000–1016, 2017.