# A combined Control by Interconnection-Model Predictive Control design for constrained Port-Hamiltonian systems

T.H. Pham, N.M.T. Vu, I. Prodan, L. Lefèvre

► To cite this version:

T.H. Pham, N.M.T. Vu, I. Prodan, L. Lefèvre. A combined Control by Interconnection-Model Predictive Control design for constrained Port-Hamiltonian systems. Systems and Control Letters, 2022, 167, pp.105336. 10.1016/j.sysconle.2022.105336. hal-04030283

HAL Id: hal-04030283

https://hal.science/hal-04030283

Submitted on 15 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# A combined *Control by Interconnection* - *Model Predictive Control* design for
Constrained Port-Hamiltonian Systems

T. H. Phamª, T. Vuᵇ, I. Prodanª, L. Lefèvreª

ªUniv. Grenoble Alpes, Grenoble INP, LCIS, 26000 Valence, France

ᵇÉcole Polytechnique Fédérale de Lausanne (EPFL), Swiss Plasma Center (SPC), CH-1015 Lausanne, Switzerland

## Abstract

This paper proposes a Control by Interconnection design, for a class of constrained Port-Hamiltonian systems, which is based on an associated Model Predictive Control optimization problem. This associated optimization problem allows to consider both state and input constraints simultaneously. Based on the first order Karush-Kuhn-Tucker optimality condition, the primal-dual gradient method is then used to build a passive feedback controller from the MPC-induced optimization problem. The resulting passive controller is coupled with the original Port-Hamiltonian system through a power-preserving interconnection, in order to guarantee the closed-loop stability. Comments on parameters tuning for the proposed control design, together with validations of the approach through simulations on a LC circuit, the simplified model of a DC-DC buck converter, and comparisons with a classical MPC design, are provided to discuss the effectiveness of the approach.

**Keywords:** Constrained port-Hamiltonian systems, Control by Interconnection, Model Predictive Control, primal-dual gradient method

## 1. Introduction

Port-Hamiltonian (PH) modeling is often a fruitful approach for the stability analysis and control design of nonlinear multiphysics systems [1, 2]. The approach is based on the modular power-preserving interconnection of passive subsystems (and external power supplies). Therefore, a PH system is intrinsically passive and the Hamiltonian function (energy, entropy, etc.) may be interpreted as a Lyapunov function to tackle the stability issue. Many control methods from the literature are developed based on this property [3], e.g. Control by Interconnection (CBI, [4]), Energy Shaping [5] or Interconnection and Damping Assignment Passivity-Based Control (IDA PBC, [6]).

Recently, various industrial applications which make use of this formalism have been shown to require constraints handling [7, 8, 9]. On the other hand, investigations on the connections between feedback and optimal control designs have a long history [10]. The *Inverse problem of optimal control* is investigated for dissipative affine nonlinear system in [11]. More recently, optimization-based control designs for PH systems without constraints are developed as linear quadratic (LQ) design in [12] or linear quadratic Gaussian (LQG) control design in [13]. In [7], an $H_\infty$ control law is proposed for a class of switched PH systems where the input saturation is considered. In [8, 9], the authors investigate the benefits of a passive

dynamical controller, designed by applying the primal-dual gradient method to finite-dimensional optimization problems. This construction guarantees the constraint satisfactions of the instantaneous input and the steady state. Furthermore, an off-line optimal controller for PH systems is designed in [14]. In [15], the optimization problem is solved by a numerical tool equivalent to a Model Predictive Control (MPC) solver which, however, does not take advantage of the PH formalism. Note that, in all these approaches except for [15], no prediction of the states is taken into account. Therefore, they can only deal with input constraints and not with state constraints which should be satisfied at all times.

To deal with this issue, a well-known method is the MPC [16]. Although the theory on linear MPC gained ground over the last decades, stability analysis and high computation effort of nonlinear MPC are still challenging. Furthermore, finding a Lyapunov function to analyze the stability of the closed loop system is one of those popular questions which is relatively simple to formulate but not trivial to solve. A possible solution for this issue is exploiting the passivity property of the closed-loop system as studied in [17, 18], where constraints on the supplied energy are added to the MPC formulation to facilitate the stability illustration. However, this technique reduces the feasibility region of the MPC optimization problem, and thus, the controller may have no solution. Moreover, MPC solves an optimization problem at each time instant, which requires a suitable optimizer and a considerable computational effort. The authors in [19] proposed an

Email address: trang.vu@epfl.ch (T. Vu)

Preprint submitted to Elsevier
instant-MPC to deal with this drawback by using the primal-dual gradient method to solve online the MPC optimization problem. As a result, the computation time can be drastically reduced, about hundred times faster. Nevertheless, in the aforementioned work, the supply rate determination for the dissipativity condition is not trivial, and the stability is not generally guaranteed.

This work aims at a control design methodology for PH systems with constraints using the advantages of MPC in combination with the PH formalism. Our work inspires from a result developed in [20] where the application of the primal-dual gradient method to a convex optimization problem leads to a passive dynamical controller. The main contribution of this work is to propose a Control by Interconnection (CBI) method combined with the MPC principles, leading to the following advantages:

* the system state constraints are taken into account. It is important to note that we do not try to find the exact MPC law with the same optimization problem, rather we are concentrating on enforcing state and input constraints satisfaction for the controlled systems.

* the proposed dynamical controller provides the instant control action without any iterative optimizer as used in MPC. This significantly reduces the computational effort.

* the stability analysis is facilitated, and the convergence of the closed-loop system is guaranteed thanks to the passivity property of the PH formulation.

The paper is organized as follows. In Sec. 2, we briefly remind the finite dimensional Port-Controlled Hamiltonian (PCH) systems definition, the primal-dual gradient method to solve optimization problems, and the problem formulation with MPC technique. In Sec. 3, we propose a dynamical feedback control design, discuss the closed-loop system stability and comment the control tuning parameters. Numerical demonstrations are shown in Sec. 4. Finally, we conclude the paper with some prospects for future work in Sec. 5.

## 2. Prerequisites

In this section, we briefly recall the definition of finite dimensional port-controlled Hamiltonian systems and the passivity with respect to the Hamiltonian function and the power conjugate input-output variables. Then the primal-dual gradient method for solving finite dimensional convex optimization problems and the MPC principle to deal with system constraints are shortly presented.

### 2.1. Finite dimensional port-controlled Hamiltonian system

In this work, we consider finite dimensional port-controlled Hamiltonian (PCH) systems described in the

following explicit input-state-output form:

$$ \left\{ \begin{aligned} \dot{\mathbf{x}}(t) &= [\mathbf{J}_x(\mathbf{x}) - \mathbf{R}_x(\mathbf{x})] \nabla H_x(\mathbf{x}) + \mathbf{G}_x(\mathbf{x}) \mathbf{u}(t), \\ \mathbf{y}(t) &= \mathbf{G}_x^\top(\mathbf{x}) \nabla H_x(\mathbf{x}), \end{aligned} \right. \quad (1) $$

where $\mathbf{x}(t) \in \mathbb{R}^n$ and $\mathbf{u}(t) \in \mathbb{R}^m$ are the state and input vectors, respectively, $\mathbf{J}_x(\mathbf{x}) = -\mathbf{J}_x^\top(\mathbf{x}) \in \mathbb{R}^{n \times n}$ is the skew-symmetric interconnection matrix, $\mathbf{R}_x(\mathbf{x}) = \mathbf{R}_x^\top(\mathbf{x}) \in \mathbb{R}^{n \times n}$ is the symmetric and non-negative dissipation matrix, $\mathbf{G}_x(\mathbf{x}) \in \mathbb{R}^{n \times m}$ is the input matrix and $H_x(\mathbf{x}) \in \mathbb{R}$ is the positive Hamiltonian, e.g. the system's energy. As one of the main properties of PH systems, the plant (1), with conjugate input $\mathbf{u}(t)$ and output $\mathbf{y}(t)$, is passive with respect to the storage function $H_x(\mathbf{x})$, since $dH_x(\mathbf{x})/dt \le \mathbf{u}^\top(\mathbf{x})\mathbf{y}(t)$. We will therefore take into account the following assumption.

**Assumption 1.** The Hamiltonian $H_x(\mathbf{x})$ is bounded from below, strictly convex, and minimized at the origin $\mathbf{x}^* = \mathbf{0}$, which is the equilibrium of the autonomous system corresponding to $\mathbf{u}(t) = \mathbf{0}$.

Also note that system (1) is completely integrable when $\mathbf{J}_x(\mathbf{x})$ satisfies the Jacobi identities [3].

### 2.2. Primal-dual gradient method

We recall hereafter the primal-dual gradient method [21] which is used to solve the following finite-dimensional optimization problem:

$$ \begin{array}{ll} \mathbf{z}^* = \underset{\mathbf{z}}{\operatorname{argmin}} & f(\mathbf{z}) \\ \text{s.t.} & \mathbf{A}_z \mathbf{z} + \mathbf{b}_z = \mathbf{0}, \\ & \mathbf{g}(\mathbf{z}) \le \mathbf{0}, \end{array} \quad (2) $$

where $\mathbf{z} \in \mathbb{R}^{n_z}$, $f(\mathbf{z}) \in \mathbb{R}$, $\mathbf{A}_z \in \mathbb{R}^{n_z \times n_z}$, $\mathbf{b}_z \in \mathbb{R}^{n_z}$, $\mathbf{g}(\mathbf{z}) \in \mathbb{R}^{n_\mu}$, and $n_z, n_\lambda, n_\mu \in \mathbb{N}$. The following assumption is necessary for a feasible optimization problem in (2).

**Assumption 2.** The cost function $f(\mathbf{z})$ is strictly convex and continuously differentiable, $\mathbf{g}(\mathbf{z})$ is convex, continuously differentiable and $\mathbf{g}(\mathbf{0}) < \mathbf{0}$.

Let $L(\mathbf{z}, \lambda, \mu) \in \mathbb{R}$ denote the Lagrangian function associated to problem (2), i.e.

$$ L(\mathbf{z}, \lambda, \mu) = f(\mathbf{z}) + \lambda^\top (\mathbf{A}_z \mathbf{z} + \mathbf{b}_z) + \mu^\top \mathbf{g}(\mathbf{z}), \quad (3) $$

with $\lambda \in \mathbb{R}^{n_\lambda}$ and $\mu \in [0, +\infty)^{n_\mu}$. For all optimal solutions $\mathbf{z}^*$ of (2), there exist $\lambda^*$ and $\mu^*$ satisfying the first-order Karush-Kuhn-Tucker (KKT) conditions [22]:

$$ \begin{array}{ll} \nabla L(\mathbf{z}, \lambda, \mu) = 0 & \\ \Leftrightarrow \begin{cases} \nabla f(\mathbf{z}^*) + \mathbf{A}_z^\top \lambda^* + \nabla \mathbf{g}^\top(\mathbf{z}^*) \mu^* = \mathbf{0}, \\ \mathbf{A}_z \mathbf{z}^* + \mathbf{b}_z = \mathbf{0}, \\ \mathbf{g}(\mathbf{z}^*) \le \mathbf{0}, \mu^* \ge \mathbf{0}, \mu^{*T} \mathbf{g}(\mathbf{z}^*) = \mathbf{0}. \end{cases} & (4) \end{array} $$

Based on the previous KKT conditions, the primal-dual gradient algorithm is described by the following dynamical
system (similar to the one proposed in [20]):

$$
\left\{
\begin{array}{l}
\tau_z \dot{\mathbf{z}}(t) = -\nabla f(\mathbf{z}) - \mathbf{A}_z^\top \lambda(t) - \nabla \mathbf{g}^\top(\mathbf{z}) \mu(t), \\
\tau_\lambda \dot{\lambda}(t) = \mathbf{A}_z \mathbf{z}(t) + \mathbf{b}_z, \\
\tau_\mu \dot{\mu}(t) = [\mathbf{g}(\mathbf{z})]_\mu^+,
\end{array}
\right.
\qquad (5)
$$

where the $i^{th}$ element ($i \in \{1, \dots, n_\mu\}$) of the vector $[\mathbf{g}(\mathbf{z})]_\mu^+ \in \mathbb{R}^{n_\mu}$ is defined as:

$$
[g_i(\mathbf{z})]_{\mu}^{+} = \begin{cases} g_i(\mathbf{z}), & \text{if } \mu_i > 0, \\ \max\{0, g_i(\mathbf{z})\}, & \text{if } \mu_i = 0, \end{cases} \quad (6)
$$

and where $\tau_z \in \mathbb{R}_+^{n_z \times n_z}$, $\tau_\lambda \in \mathbb{R}_+^{n_\lambda \times n_\lambda}$ and $\tau_\mu \in \mathbb{R}_+^{n_\mu \times n_\mu}$ are symmetric positive matrices, characterizing the different timescales appearing in the dynamics.

**Proposition 1.** *The states of the dynamics (5) converge to the set of equilibrium points.*

Proof. See Appendix A. □

Since the equilibrium points of the dynamics (5) are also the solutions of the KKT equations (4), any numerical integration method for (5) can be used to solve the optimization problem (2). Moreover, the autonomous system (5) may be cast as a closed loop PH system, which simplifies the demonstration of the convergence of the states to the equilibrium [20].

2.3. *Model predictive control*

In the following, we briefly recall the general
optimization problem formulation for constrained systems
using MPC technique. We also show how the optimization
problem is transformed to fit into the finite dimensional
framework studied in this work.

Let $\mathbf{U}(t) = \{\mathbf{u}(\cdot|t) : [t, t+h] \to \mathbb{R}^m : \tau \mapsto u(\tau|t)\}$ and $\mathbf{X}(t) = \{\mathbf{x}(\cdot|t) : [t, t+h] \to \mathbb{R}^n : \tau \mapsto x(\tau|t)\}$ denote respectively the sets of input and state functions for current time $\tau \in [t, t+h]$ (i.e. over a prediction horizon $h$), where $\mathbf{x}(\tau|t)$ and $\mathbf{u}(\tau|t)$ are respectively the values of the system states and inputs at the time instant $\tau \in [t, t+h]$ which are predicted at time $t$. Then consider the following constrained optimization problem:

$$
\{\mathbf{U}^*(t), \mathbf{X}^*(t)\} = \underset{\mathbf{U}(t), \mathbf{X}(t)}{\operatorname{argmin}} V_f(\mathbf{x}(t+h|t)) + \int_t^{t+h} l_{xu}(\mathbf{x}, \mathbf{u}) \, d\tau \quad (7a)
$$

$$
\text{s.t. } \dot{\mathbf{x}}(\tau|t) = [\mathbf{J}_x(\mathbf{x}) - \mathbf{R}_x(\mathbf{x})] \nabla H_x(\mathbf{x}) \\
\qquad + \mathbf{G}_x(\mathbf{x}) \mathbf{u}(\tau|t), \quad \forall \tau \in [t, t+h], \qquad (7b)
$$

$$
\mathbf{g}(\mathbf{x}, \mathbf{u}) \le \mathbf{0}, \quad \forall \tau \in [t, t+h], \qquad (7c)
$$

The MPC feedback control at time $t$, is then defined as $\mathbf{u}_{\text{MPC}}(t) = \mathbf{u}^*(t|t)$ where $\mathbf{u}^*(t|t)$ denotes the value of the optimal input trajectory $\mathbf{u}^*(\tau|t)$ for the current time value $\tau = t$. In (7), the stage and final cost functions $l_{xu}(\mathbf{x}, \mathbf{u})$ and $V_f(\mathbf{x}(t+h|t))$ penalize the state error and the control deviation.

*Problem linearization.* Note that (7) is an infinite-
dimensional optimization problem which is not the case
of the problem (2) solved by the primal-dual gradient

method described in section 2.2. Therefore it is
necessary to approximate (7) by a finite-dimensional
optimization problem. In this work, simple piecewise-
constant approximations are used for the state and control
time profiles on the prediction horizon [t, t + h]. Hence,
we will consider **x**(*τ*|*t*) = ∑<sub>*k* = 1</sub><sup>*N*</sup> **x**(*k*|*t*) *β*<sub>*k*</sub>(*τ*), **u**(*τ*|*t*) = ∑<sub>*k* = 1</sub><sup>*N*</sup> **u**(*k*|*t*) *β*<sub>*k*</sub>(*τ*), where *β*<sub>*k*</sub>(*τ*) are the window functions
described as:

$$
\beta_k(\tau) = \begin{cases} 1, & \text{if } t + (k - 1)\Delta t \le \tau < t + k\Delta t, \\ 0, & \text{else, } \forall k \in \{1, \dots, N\}, \end{cases} \quad (8)
$$

with time step $\Delta t$ and $N = \frac{h}{\Delta t} \in \mathbb{N}$. Existence of the solution of the MPC formulation defined in (7) also requires that the plant states $\mathbf{x}(t)$ are fully observable [23].

*Problem discretization.* On the other hand, regarding the linear equality constraint in (2), the plant (1) or (7b) also needs to be represented in a linearized discrete-time form:

$$
\mathbf{x}(k+1|t) = \mathbf{A}\mathbf{x}(k|t) + \mathbf{B}\mathbf{u}(k|t), \quad (9)
$$

where the matrices $\mathbf{A} \in \mathbb{R}^{n \times n}$ and $\mathbf{B} \in \mathbb{R}^{n \times m}$ are constants and where $\mathbf{x}(k|t)$ and $\mathbf{u}(k|t)$ (with $k \in \mathbb{N}$) denote respectively the predicted values of the state and input variables at instant $t+k\Delta t$. This linear discrete-time model is obtained through linearization and subsequent structure-preserving time discretization. The latter is a symplectic Runge-Kutta method defined in order to preserve the intrinsic geometric interconnection (Dirac) structure of the original PCH system [24]. In this approach, the local error of the stored energy is consistent with the numerical integration scheme [24, Theorem 2]. This discrete scheme is briefly recalled in Appendix B.

We consider hereafter the recursive construction of a discrete-time optimal open-loop state and control sequence $\mathbf{z}(t) \in \mathbb{R}^{(m+n)N}$:

$$
\mathbf{z}(t) = [\mathbf{u}^{\top}(0|t), \mathbf{u}^{\top}(1|t), \dots, \mathbf{u}^{\top}(N-1|t), \\
\qquad \mathbf{x}^{\top}(1|t), \dots, \mathbf{x}^{\top}(N|t)]^{\top} \quad (10)
$$

at each time instant $t$ over a finite prediction horizon
$[t, t + \Delta t, \dots, t + N \Delta t]$, $N \in \mathbb{Z}_+$. The feedback control law
of the plant is thus the first element of $\mathbf{z}(t)$:

$$
\mathbf{u}(t) = \mathbf{u}(0|t) = \mathbf{E}\mathbf{z}(t), \qquad (11)
$$

with $\mathbf{E} = [\mathbf{I}_m \mathbf{0}] \in \mathbb{R}^{m \times (m+n)N}$. Moreover, the equivalent MPC law is:

$$
\mathbf{u}_{\text{MPC}}(t) = \mathbf{E}\mathbf{z}^*(t), \qquad (12)
$$

where $\mathbf{z}^*(t)$ is the optimal solution of the following
optimization problem:

$$
\mathbf{z}^*(t) = \underset{\mathbf{z}(t)}{\operatorname{argmin}} f(\mathbf{z}) \qquad (13\text{a})
$$

$$
\text{s.t.} \quad \mathbf{A}_z \mathbf{z}(t) + \mathbf{B}_z \mathbf{x}(t) = \mathbf{0}, \qquad (13\text{b})
$$

$$
\mathbf{g}(\mathbf{z}) \le \mathbf{0}. \qquad (13c)
$$

The matrices $\mathbf{A}_z \in \mathbb{R}^{n_z \times (m+n)N}$ and $\mathbf{B}_z \in \mathbb{R}^{n_z \times n}$ are
defined as:

$$ \mathbf{A}_z = \begin{bmatrix} \mathbf{B} & 0 & \dots & 0 & -\mathbf{I}_n & 0 & \dots & 0 \\ 0 & \mathbf{B} & \dots & 0 & \mathbf{A} & -\mathbf{I}_n & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & \mathbf{B} & 0 & \dots & \mathbf{A} & -\mathbf{I}_n \end{bmatrix}, \quad (14a) $$

$$ \mathbf{B}_z = \begin{bmatrix} \mathbf{A} \\ \mathbf{0} \end{bmatrix}, \quad (14b) $$

where the matrices **A** and **B** are defined in (9). The cost function $f(\mathbf{z})$ now corresponds to the discrete-time form of the cost in (7a), i.e.:

$$ f(\mathbf{z}) = V_f(\mathbf{x}(N|t)) + \Delta t \sum_{k=0}^{N-1} l_{xu}(\mathbf{x}(k|t), \mathbf{u}(k|t)). \quad (15) $$

**Remark 1.** Usually, the cost functions

$l_{xu}(\mathbf{x}(k|t), \mathbf{u}(k|t))$ and $V_f(\mathbf{x}(N|t))$ are chosen

quadratic, i.e., $l_{xu} = \mathbf{x}^\top \mathbf{Q}_x \mathbf{x} + \mathbf{u}^\top \mathbf{Q}_u \mathbf{u}$ and $V_f = \mathbf{x}^\top \mathbf{Q}_f \mathbf{x}$, where the weight matrices $\mathbf{Q}_x \in \mathbb{R}^{n \times n}$, $\mathbf{Q}_u \in \mathbb{R}^{m \times m}$ and $\mathbf{Q}_f \in \mathbb{R}^{n \times n}$ are symmetric and positive. Hence, the cost function $f(\mathbf{z})$ in (13a) is also quadratic, i.e., $f(\mathbf{z}) = \mathbf{z}^\top \mathbf{Q}_z \mathbf{z}$, where the weight matrix $\mathbf{Q}_z \in \mathbb{R}^{(m+n)N \times (m+n)N}$ has the block-diagonal form:

$$ \mathbf{Q}_z = \mathrm{diag} \left\{ \mathbf{Q}_u, \dots, \mathbf{Q}_u, \mathbf{Q}_x, \dots, \mathbf{Q}_x, \frac{\mathbf{Q}_f}{\Delta t} \right\}. \quad (16) $$

**Remark 2.** More linear equality constraints can easily be taken into account in the optimization problem (13) by adding more rows in the matrices $\mathbf{A}_z$ and $\mathbf{B}_z$.

### 3. Main idea

#### 3.1. Controller design

This work focuses on the design of a dynamic feedback control law, named CbI-MPC, which on the one hand stabilizes the state vector $\mathbf{x}(t)$ of system (1) to the origin $\mathbf{x}^e = \mathbf{0}$ (using Control-by-Interconnection (CbI) technique [2]), and on the other hand respects inequality constraints $\mathbf{g}(\mathbf{x}, \mathbf{u}) \le \mathbf{0}$, both on the system state and input (using MPC technique).

The controller dynamics are derived from the primal-dual gradient method for the MPC optimization problem (13) (see also Fig. 1). Note that step reference tracking is a particular case of this work. However, time-varying reference tracking or economic MPC are excluded.

From (5) and (13), the controller dynamics are derived as:

$$ \left\{ \begin{array}{l} \tau_z \dot{\mathbf{z}}(t) = -\nabla f(\mathbf{z}) - \mathbf{A}_z^\top \lambda(t) - \nabla \mathbf{g}^\top(\mathbf{z}) \mu(t), \\ \tau_\lambda \dot{\lambda}(t) = \mathbf{A}_z \mathbf{z}(t) + \mathbf{B}_z \mathbf{x}(t), \\ \tau_\mu \dot{\mu}(t) = [\mathbf{g}(\mathbf{z})]_\mu^+, \end{array} \right. \quad (17) $$

Unlike the autonomous system (5), the controller system (17) has an input $\mathbf{u}_c(t)$ to get the plant information and consequently a corresponding output $\mathbf{y}_c(t)$ for the control action. To apply the CbI technique, the controller dynamics (17) must be a passive system where its input $\mathbf{u}_c(t)$ and the output $\mathbf{y}_c(t)$ are power-conjugate variables, i.e., their product is the supplied power to the controller system. The plant (1) and the controller

<figure><img src="image_1.png" /><figcaption>Figure 1: Dynamic controller coupled to the PH system using CbI</figcaption></figure>

(17) is then coupled together using a power-preserving interconnection, in order to form a passive closed loop system. A simple form of such interconnection is defined as:

$$ \left\{ \begin{aligned} \mathbf{u}_c(t) &= \mathbf{y}(t), \\ \mathbf{u}(t) &= -\mathbf{y}_c(t). \end{aligned} \right. \quad (18) $$

According to (11), (17) and (18), the input $\mathbf{u}_c(t)$ and the output $\mathbf{y}_c(t)$ should respect two following conditions:

$$ \mathbf{u}_c(t) = \mathbf{x}(t), \quad (19) $$

$$ \mathbf{y}_c(t) = -\mathbf{E}\mathbf{z}(t). \quad (20) $$

**Remark 3.** Condition (19) requires a direct construction of the plant state $\mathbf{x}(t)$ from the plant output $\mathbf{y}(t)$, which is, in general, not trivial, for instance in the case of under actuated systems (where the input dimension is smaller than the state one). However, this issue can be tackled using an additional state observer defined in such a way that the augmented system, including the plant and the observer, is also passive (see [25, 26, 27] and the references therein). As a result, the main principle of the presented CbI-MPC controller design will not be affected. However, some parameter tuning may need to be adapted according to the augmented system. This will be discussed with more details in Sec. 3.3.

In this work, for the sake of simplicity, such observer is not considered and thus the following assumption is admitted in order to derive the state $\mathbf{x}(t)$ from the plant output $\mathbf{y}(t)$ in (1).

**Assumption 3.** There exists an invertible constant matrix $\mathbf{M} \in \mathbb{R}^{n \times n}$ such that:

$$ \mathbf{y}(t) = \mathbf{G}_x^\top(\mathbf{x}) \nabla H_x(\mathbf{x}) = \mathbf{M}\mathbf{x}(t). \quad (21) $$

This assumption implies that the plant input, output and state have the same dimension, i.e., $m = n$.

Similar to Appendix A, the Hamiltonian function $H_r(\mathbf{r})$
of the controller dynamics (5) is simply chosen as:

$$H_r(\mathbf{r}) = \frac{1}{2}\mathbf{r}_z^\top(t)\tau_z^{-1}\mathbf{r}_z(t) + \frac{1}{2}\mathbf{r}_\lambda^\top(t)\tau_\lambda^{-1}\mathbf{r}_\lambda(t) + \frac{1}{2}\mathbf{r}_\mu^\top(t)\tau_\mu^{-1}\mathbf{r}_\mu(t), \quad (22)$$

with the transformed state vector $\mathbf{r}(t) \in \mathbb{R}^{3nN+n_\mu}$ defined by:

$$\mathbf{r}(t) = \begin{bmatrix} \mathbf{r}_z(t) \\ \mathbf{r}_\lambda(t) \\ \mathbf{r}_\mu(t) \end{bmatrix} = \begin{bmatrix} \tau_z \mathbf{z}(t) \\ \tau_\lambda \lambda(t) \\ \tau_\mu \mu(t) \end{bmatrix}. \quad (23)$$

Based on (17)-(18), and (21)-(23), the controller dynamics are rewritten as:

$$\left\{ \begin{aligned} \dot{\mathbf{r}}(t) &= \mathbf{f}_r(\mathbf{r}) + \begin{bmatrix} \mathbf{0} \\ \mathbf{B}_z \\ \mathbf{0} \end{bmatrix} \mathbf{M}^{-1} \mathbf{u}_c(t), \\ \mathbf{y}_c(t) &= \mathbf{M}^{-\top} [\mathbf{0} \quad \mathbf{B}_z^\top \quad \mathbf{0}] \nabla H_r(\mathbf{r}), \end{aligned} \right. \quad (24)$$

with

$$\mathbf{f}_r(\mathbf{r}) = \begin{bmatrix} -\mathbf{A}_z^\top \partial_{\mathbf{r}_\lambda} H_r(\mathbf{r}) - \nabla f(\mathbf{z}) - \nabla^\top \mathbf{g}(\mathbf{z}) \mu(t) \\ \mathbf{A}_z \partial_{\mathbf{r}_z} H_r(\mathbf{r}) \\ [\mathbf{g}(\mathbf{z})]_\mu^\top \end{bmatrix}. \quad (25)$$

It is important to note that the requirement (20) can not be respected according to equations (11), (18) and (24). As a result, $\mathbf{u}(t)$ does not satisfy the constraint in (7c) even though $\mathbf{z}(t)$ satisfies the constraint (13c). In order to tackle this issue, we propose in the following to add an extra term $\mathbf{G}_z(\mathbf{z}, \lambda)$ to the input matrix of the controller dynamics (24) such that:

$$\left\{ \begin{aligned} \dot{\mathbf{r}}(t) &= \mathbf{f}_r(\mathbf{r}) + \begin{bmatrix} \mathbf{G}_z(\mathbf{z}, \lambda) \\ \mathbf{B}_z \\ \mathbf{0} \end{bmatrix} \mathbf{M}^{-1} \mathbf{u}_c(t), \\ \mathbf{y}_c(t) &= \mathbf{M}^{-\top} [\mathbf{G}_z^\top(\mathbf{z}, \lambda) \quad \mathbf{B}_z^\top \quad \mathbf{0}] \nabla H_r(\mathbf{r}), \end{aligned} \right. \quad (26)$$

where $\mathbf{G}_z(\mathbf{z}, \lambda) \in \mathbb{R}^{2nN \times n}$ is non-linear and satisfies the following condition:

$$-\mathbf{M}^{\top} \mathbf{E} \mathbf{z}(t) = \mathbf{G}_z^{\top}(\mathbf{z}, \lambda) \mathbf{z}(t) + \mathbf{B}_z^{\top} \lambda(t). \quad (27)$$

The matrices $\mathbf{M}$, $\mathbf{E}$ and $\mathbf{B}_z$ are defined in (21), (11) and (14b), respectively. The condition (27) implies that the control law given in (17)-(18) is equal to the first element of $\mathbf{z}(t)$ at all time. Note that with given values of $\mathbf{z}$ and $\lambda$, (27) is actually a linear equation of $\mathbf{G}_z(\mathbf{z}, \lambda)$ (see discussion in Sec. 3.3).

**Proposition 2.** The controller system defined by (25), (26) and (27) is passive.

*Proof.* From(22), (23) and (26), we have:

$$\begin{aligned} \dot{H}_r(\mathbf{r}) &= \nabla^\top H_r(\mathbf{r}) \dot{\mathbf{r}}(t) \\ &= -\mathbf{z}^\top(t) \nabla f(\mathbf{z}) - \mathbf{z}^\top(t) \nabla^\top \mathbf{g}(\mathbf{z}) \mu(t) \\ &\quad + \mu^\top(t) [\mathbf{g}(\mathbf{z})]_\mu^\top + \mathbf{y}_c^\top(t) \mathbf{u}_c(t). \end{aligned} \quad (28)$$

With derivations similar to those in appendix Appendix A to obtain (A.11), we obtain:

$$\dot{H}_r(\mathbf{r}) \le \mathbf{y}_c^\top(t) \mathbf{u}_c(t), \quad (29)$$

and thus, the proposition is concluded. □

**Remark 4 (Convergence).** Assume there exists an equilibrium $\mathbf{r}^*(\mathbf{u}_c)$ of (26), which includes the predicted input and state vectors completely respecting the constraints. Despite the controller's passivity, the convergence of the controller state $\mathbf{r}(t)$ to $\mathbf{r}^*(\mathbf{u}_c)$ is not guaranteed. Indeed, using Proposition 1 and the corresponding proof in Appendix A, the shifted controller state is defined as $\tilde{\mathbf{r}}(t) = \mathbf{r}(t) - \mathbf{r}^*(\mathbf{u}_c)$, which leads to the shifted controller dynamics:

$$\dot{\tilde{\mathbf{r}}}(t) = \mathbf{f}_r(\mathbf{r}) - \mathbf{f}_r(\mathbf{r}^*(\mathbf{u}_c)) + \begin{bmatrix} \mathbf{G}_z(\mathbf{z}, \lambda) \\ \mathbf{0} \end{bmatrix} (\mathbf{z}^*, \lambda^*) \mathbf{u}_c.$$

Using the Hamiltonian $\tilde{H}_r(\tilde{\mathbf{r}})$ defined in (A.4) with the result proved in (A.11), we derive that:

$$\dot{\tilde{H}}_r(\tilde{\mathbf{r}}) \le \tilde{\mathbf{z}}^\top(t) [\mathbf{G}_z(\mathbf{z}, \lambda) - \mathbf{G}_z(\mathbf{z}^*, \lambda^*)] \mathbf{u}_c.$$

Since the right-hand side of the previous inequality is not generally negative, the shifted controller dynamics are not proved passive, and thus, the convergence of the state $\mathbf{r}(t)$ to the equilibrium $\mathbf{r}^*(\mathbf{u}_c)$ is not ensured. Nonetheless, in simulations we observe empirically that convergence holds.

**Remark 5 (Optimality).** Due to the presence of the nonlinear matrix $\mathbf{G}_z^\top(\mathbf{z}, \lambda)$, it is not easy to find the optimization problem corresponding to the controller dynamics (26) through the relation using the primal-dual gradient method presented in the section 2.2. However, since constraint consideration is the main objective in this work, finding such an equivalent optimization problem is not mandatory.

### 3.2. Closed-loop system

Based on the previously designed controller, the closed-loop system is defined by coupling the plant (1) and the controller dynamics (26) through the power-preserving interconnection (18). The resulting closed loop system reads:

$$\left\{ \begin{aligned} \begin{bmatrix} \dot{\mathbf{r}}(t) \\ \dot{\mathbf{r}}_z(t) \\ \dot{\mathbf{r}}_\lambda(t) \end{bmatrix} &= [\mathbf{J}(\mathbf{x}, \mathbf{z}, \lambda) - \mathbf{R}(\mathbf{x})] \begin{bmatrix} \partial_{\mathbf{x}} H(\mathbf{x}, \mathbf{r}) \\ \partial_{\mathbf{r}_z} H(\mathbf{x}, \mathbf{r}) \\ \partial_{\mathbf{r}_\lambda} H(\mathbf{x}, \mathbf{r}) \end{bmatrix} \\ &\quad + \begin{bmatrix} \mathbf{0} \\ -\nabla f(\mathbf{z}) - \nabla^\top \mathbf{g}(\mathbf{z}) \mu(t) \end{bmatrix}, \\ \dot{\mathbf{r}}_\mu(t) &= [\mathbf{g}(\mathbf{z})]_\mu^\top, \end{aligned} \right. \quad (30)$$

where $\mathbf{r}_z(t) \in \mathbb{R}^{2nN}$, $\mathbf{r}_\lambda(t) \in \mathbb{R}^{nN}$, $\mathbf{r}_\mu(t) \in \mathbb{R}^{n_\mu}$ and $\mathbf{r}(t) \in \mathbb{R}^{3nN+n_\mu}$ are defined in (23); $\mathbf{J}(\mathbf{x}, \mathbf{z}, \lambda)$, $\mathbf{R}(\mathbf{x}) \in \mathbb{R}^{(n+3nN) \times (n+3nN)}$ and the closed loop Hamiltonian $H(\mathbf{x}, \mathbf{r})$ are defined as follows:

$$\mathbf{J} = \begin{bmatrix} \mathbf{J}_x & -\mathbf{G}_x \mathbf{M}^{-T} \mathbf{G}_z^\top & -\mathbf{G}_x \mathbf{M}^{-T} \mathbf{B}_z^\top \\ \mathbf{G}_z \mathbf{M}^{-1} \mathbf{G}_x^\top & \mathbf{0} & -\mathbf{A}_z^\top \\ \mathbf{B}_z \mathbf{M}^{-1} \mathbf{G}_x^\top & \mathbf{A}_z & \mathbf{0} \end{bmatrix}, \quad (31a)$$

$$\mathbf{R} = \text{blockdiag}\{\mathbf{R}_x(\mathbf{x}), \mathbf{0}, \mathbf{0}\}, \quad (31b)$$

$$H = H_x(\mathbf{x}) + H_r(\mathbf{r}) \quad (31c)$$

with the Hamiltonians $H_x(\mathbf{x})$ and $H_r(\mathbf{r})$ given in (1) and (22). Note that the term $-\nabla f(\mathbf{z}) - \nabla^\top \mathbf{g}(\mathbf{z}) \mu(t)$ contributes to the dissipation of the closed-loop system. The stability and the convergence of the closed-loop system are proved in the following proposition.
**Proposition 3.** The closed-loop system (30)-(31):

i) is passive.

ii) converges to the origin if $ker (\mathbf{A}_z^\top) = \{0\}$.

*Proof.* i. Since the plant (1) and the controller system (26) are passive, and the interconnection (18) is power-preserving, the closed-loop system is also passive, i.e., $\dot{H}(\mathbf{x}, \mathbf{r}) \le 0$ [2].

ii. Consequently, according to the LaSalle's invariance principle, the states vector of the closed-loop system (30) converges to the largest invariant set $\mathbb{M}$ such that

$$ \mathbb{M} = \{ (\mathbf{x}, \mathbf{r}) | \dot{H}(\mathbf{x}, \mathbf{r}) = 0 \}. $$

In this largest invariant subset, we may conclude:

$$ (\mathbf{r}_z(t), \mathbf{r}_\mu(t)) = \mathbf{0}, \forall (\mathbf{x}, \mathbf{z}, \lambda, \mu) \in \mathbb{M} \quad (32) $$

$$ \Rightarrow \mathbf{r}_\lambda(t) = \mathbf{0}, \quad (33) $$

$$ \Rightarrow \nabla H_x(\mathbf{x}) = \mathbf{0} \Leftrightarrow \mathbf{x}(t) = \mathbf{0}, \quad (34) $$

(32) thanks to Assump. 2, (6), (28) and (A.10)

(33) thanks to (27), (30)-(31) and $ker (\mathbf{A}_z^\top) = \{0\}$

(34) thanks to Assumption 3, (14b) and (30)-(31)

Finally, we obtain $(\mathbf{x}(t), \mathbf{r}(t)) \underset{t \to \infty}{\longrightarrow} \mathbf{0}$ which concludes the proposition. $\square$

### 3.3. Parameter tuning and discussions

The efficiencies of the proposed controller depend on the discrete-time system model (9), the prediction step $\Delta t$, the prediction horizon $N$, the cost function $f(\mathbf{z})$ in (13a), the non-linear matrix $\mathbf{G}_z(\mathbf{z}, \lambda)$ in (27), the timescales matrices $(\tau_z, \tau_\lambda, \tau_\mu)$ and the initial controller states $(\mathbf{z}(0), \lambda(0), \mu(0))$.

* It is worth noting that we are dealing with **continuous** systems for both the plant and the controller. The discrete-time scheme (9) with respect to the time step $\Delta t$ is only used to define the finite-dimensional MPC optimization problem (13). Choosing the appropriate time-discretization scheme for the constrained optimal control is a hard question which will not be rigorously discussed in this work. However, different methods developed for PH systems should be used to preserve intrinsic system properties, e.g., the power-preserving structure and the energy conservation, as mentioned in Appendix B or discussed in [24].

* The choices of the prediction horizon $N$ and the cost function $f(\mathbf{z})$ in (13a) are not specific features of the proposed controller. They are key challenges for MPC designs. It is indeed not trivial to select these parameters in order to obtain a feasible optimization problem. In practice, the “trials and errors” approach is adopted the most frequently, combined with the extensive use of numerical simulations. When the MPC optimization problem is not feasible, no specific parameters tuning direction can be determined,

since the MPC solution does not exist. One of the advantages of the proposed CbI-MPC method is precisely that the closed-loop system behaviour may be obtained with an arbitrary parameter choice, thanks to the constraint relaxation (gradient method). This provides a guideline to adjust these tuning parameters which will be shown during the control implementation in the next section.

* The timescales matrices $(\tau_z, \tau_\lambda, \tau_\mu)$ are chosen with respect to the time constant of the controlled system. If the time scales are too high, the controller dynamics are much slower than the plant dynamics. Therefore, the constraints may be seriously violated. On the contrary, if they are small enough, the controller dynamics, in theory, rapidly converge to the instantaneous equilibrium corresponding to the input $\mathbf{u}_c(t)$. Hence, the constraints on the predicted plant dynamics, input and output are respected before the control action application. This implies that the constraints are better taken into account. Moreover, if the timescales are small enough, $\mathbf{G}_c(\mathbf{z}, \lambda) \simeq \mathbf{0}$, the control law is then directly defined in (11) and the controller states $(\mathbf{z}(t), \lambda(t), \mu(t))$ quickly converge to the optimum values $(\mathbf{z}^*(t), \lambda^*(t), \mu^*(t))$ given in (4). In that case, the control law will converge to the conventional MPC law given in (12). However, in practice, small time scales will increase the computational time which may exceed the time limit, e.g. in real time applications. The compromise between performance and rapidity thus depends on each application.

* The matrix $\mathbf{G}_z(\mathbf{z}, \lambda)$, which must be computed at each time step, is a solution of the $n$ linear equations (27). This matrix has $2nN \times n$ elements and therefore many degrees of freedom exist for its choice. The detailed analysis of the influences of these choices on the control performance, which is quite complicated due the nonlinearity, is beyond the scope of this paper and left for future research. To the best of our knowledge, in CbI technique, the input matrix is usually chosen constant due to the fact that no input constraints are considered so far. This work hence confirms the flexibility of the CbI method, which can be further developed for more applications in the future.

* The influence of the initial controller states $(\mathbf{z}(0), \lambda(0), \mu(0))$ on the system stability is less important than the previous tunable parameters. $\mathbf{z}(0)$ just needs to satisfy the constraints in the optimization problem (13), and $\mu(0)$ must not be negative. However, bad choices of these parameters may lead to an invalid $\mathbf{G}_z(\mathbf{z}, \lambda)$ in the condition (27). A possible solution is to choose the initial controller states $\mathbf{r}(0)$ at the equilibrium $\mathbf{r}^*(\mathbf{u}_c)$ of the controller dynamics (26) where $\mathbf{u}_c = \mathbf{x}(0)$.
Besides, regarding Remark. 3, the proposed controller design can also be extended to general systems where $m \neq n$. By adding an appropriate observer, e.g. PH structure-preserving observer, we can guarantee the passivity property of the plant-observer augmented system. Similar ideas of such observer design are presented in [26, 25]. However, the output of these augmented systems is the difference between the plant output $\mathbf{y}$ and the estimated output $\hat{\mathbf{y}}$ which can not be directly used by the proposed CbI-MPC controller. In an ongoing work, we define new observer conjugate input-output pairs so-that the estimated state $\hat{\mathbf{x}}$ can be easily extracted from the observer outputs while the augmented system remains passive. The proposed observer will facilitate state-feedback controller design. In particular, controller laws based on CbI technique will take charge of stabilizing the closed loop system, as well as ensuring the convergence of the observer.

In order to illustrate the effectiveness of the proposed CbI-MPC method, we will compare in the next section the performances of different control methods through a qualitative evaluation with four criteria: computational effort, input constraint consideration, state constraint consideration and stability illustration (see Table 1).

## 4. Numerical examples

In the following we validate the proposed method over an electrical system which is in the PH system class defined Subsection 2.1 (more precisely, it is linear, with $\mathbf{u}(t) \in \mathbb{R}^n$ and $\mathbf{G}(\mathbf{x}) = \mathbf{G}$ is constant and invertible).

### 4.1. LC circuit

<figure><img src="image_2.png" /><figcaption>Figure 2: Simple LC circuit with 2 control signal $v_i(t)$ and $i_o(t)$.</figcaption></figure>

A LC circuit with two control inputs is described in Fig. 2. Usual Kirchoff's balance equations may be written in the form of the following PH system:

$$ \begin{bmatrix} \dot{\phi}(t) \\ \dot{q}(t) \end{bmatrix} = \mathbf{JQ} \begin{bmatrix} \phi(t) \\ q(t) \end{bmatrix} + \begin{bmatrix} v_i(t) \\ i_o(t) \end{bmatrix}, \quad (35) $$

where $\phi(t) \in \mathbb{R}$ is the magnetic flux of the inductance $L$, $q(t) \in \mathbb{R}$ is the electric charge of the capacitance $C$, and the

matrices $\mathbf{J}, \mathbf{Q} \in \mathbb{R}^{2 \times 2}$ are given as: $\mathbf{J} = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$, $\mathbf{Q} = \text{diag} \left\{ \frac{1}{L}, \frac{1}{C} \right\}$. Since the CbI-MPC controller of Sec. 3.2

has been designed to stabilize the PH system state around the origin, a change of state variables is considered for (35), i.e. it shifts the desired (reference) equilibrium value of

the state to the origin. Therefore, the shifted state vector $\mathbf{x}(t) \in \mathbb{R}^2$, the corresponding input vector $\mathbf{u}(t) \in \mathbb{R}^2$ and Hamiltonian function $H_x(\mathbf{x})$ are given as:

$$ \begin{aligned} \mathbf{x}(t) &= [\phi(t) + L i_o^* \quad q(t) - C v_L^*]^\top, \\ \mathbf{u}(t) &= [v_i(t) - v_i^* \quad i_o(t) - i_o^*]^\top, \\ H_x(\mathbf{x}) &= \frac{1}{2} \mathbf{x}^\top(t) \mathbf{Q} \mathbf{x}(t). \end{aligned} $$

System dynamics (35) then read:

$$ \left\{ \begin{aligned} \dot{\mathbf{x}}(t) &= \mathbf{J} \nabla H_x(\mathbf{x}) + \mathbf{u}(t), \\ \mathbf{y}(t) &= \nabla H_x(\mathbf{x}), \end{aligned} \right. \quad (36) $$

The following constraints of the state and input will be considered: $\mathbf{x}_{min} \leq \mathbf{x}(t) \leq \mathbf{x}_{max}$, $\mathbf{u}_{min} \leq \mathbf{u}(t) \leq \mathbf{u}_{max}$.

**Remark 6.** The LC circuit can be considered as a simplified buck DC-DC converter described in Fig. 3 where $L, C, V_{dc}, r$ and $S_1, S_2, S_3$ denotes, respectively, the inductance, the capacitance, the input DC voltage, the resistance load and the ideal switches [28]. Usually, the switches are alternatively switched at high frequency by using the Pulse Width Modulation technique. For simplicity, we can consider a slower timescale where the input voltage is represented by the continuous average value $v_i(t)$. Moreover, according to the studied example Fig. 3, the passive load is replaced by an active current source $i_o(t)$.

<figure><img src="image_3.png" /><figcaption>Figure 3: DC-DC buck converter [28].</figcaption></figure>

### 4.2. Simulation results

In the following simulations, the results are obtained using both the MPC and the CbI-MPC laws. The simulations are implemented using MATLAB 2017b, and the MPC optimization problem is solved using the quadprog function. The values of the plant, controller and simulation parameters are given in Table 2. Furthermore, we simply use the mid-point discretization method to determine the constant matrices $\mathbf{A}$ and $\mathbf{B}$ in (9):

$$ \mathbf{A} = [2\mathbf{I}_2 - \Delta t \mathbf{JQ}]^{-1} [2\mathbf{I}_2 + \Delta t \mathbf{JQ}], $$

$$ \mathbf{B} = [2\mathbf{I}_2 - \Delta t \mathbf{JQ}]^{-1} 2\Delta t \mathbf{I}_2. $$

The cost function $f(\mathbf{z})$ defined in (13a) is chosen quadratic as presented in Remark 1, i.e.,

$$ f(\mathbf{z}) = \mathbf{z}^\top(t) \mathbf{Q}_z \mathbf{z}(t). $$

Three simulation scenarios are considered as presented in Table 3: small limits of inputs, small limits of inputs and states, and critical (even smaller) limits of inputs and states, respectively. In all cases, the controller equilibrium
<table><caption>Table 1: Qualitative comparison of different control methods.</caption><thead><tr><th>Criteria</th><th>MPC [16]</th><th>Instant MPC [19]</th><th>Optimal CbI [20]</th><th>CbI-MPC</th></tr></thead><tbody><tr><td>Computational effort</td><td>high</td><td>medium</td><td>low</td><td>medium</td></tr><tr><td>Input constraint consideration</td><td>yes</td><td>yes</td><td>yes</td><td>yes</td></tr><tr><td>State constraint consideration</td><td>yes</td><td>yes</td><td>no</td><td>yes</td></tr><tr><td>Stability illustration</td><td>hard</td><td>hard</td><td>easy</td><td>easy</td></tr></tbody></table>

<table><caption>Table 2: Parameter values.</caption><thead><tr><th>Description</th><th>Notation</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td colspan="4">System</td></tr><tr><td>Inductance</td><td>L</td><td>1</td><td>[H]</td></tr><tr><td>Capacitance</td><td>C</td><td>1</td><td>[F]</td></tr><tr><td>State and input dim.</td><td>n</td><td>2</td><td></td></tr><tr><td colspan="4">Controller</td></tr><tr><td>Prediction time step</td><td>Δt</td><td>0.5</td><td>[s]</td></tr><tr><td>Prediction horizon</td><td>N</td><td>10</td><td></td></tr><tr><td>Weight matrix</td><td>Q<sub>z</sub></td><td>I<sub>40</sub></td><td></td></tr><tr><td>Time scale matrices</td><td>τ<sub>z</sub></td><td>0.01 × I<sub>40</sub></td><td></td></tr><tr><td></td><td>τ<sub>λ</sub></td><td>0.01 × I<sub>20</sub></td><td></td></tr><tr><td></td><td>τ<sub>μ</sub></td><td>0.01 × I<sub>80</sub></td><td></td></tr><tr><td colspan="4">Simulation</td></tr><tr><td>Simulation duration</td><td></td><td>5</td><td>[s]</td></tr><tr><td>Initial state</td><td>x(0)</td><td>0.8 × I<sub>2</sub></td><td></td></tr></tbody></table>

<table><caption>Table 3: Simulation scenarios.</caption><thead><tr><th>Description</th><th>Scenario 1</th><th>Scenario 2</th><th>Scenario 3</th></tr></thead><tbody><tr><td><b>u</b><sup>max</sup><sub>max</sub></td><td>[1 1]</td><td>[1 1]</td><td>[0.35 0.35]</td></tr><tr><td><b>u</b><sup>min</sup><sub>min</sub></td><td>−[0.7 0.4]</td><td>−[0.7 0.4]</td><td>−[0.7 0.4]</td></tr><tr><td><b>x</b><sup>max</sup><sub>max</sub></td><td>[1 1]</td><td>[1 1]</td><td>[1 1]</td></tr><tr><td><b>x</b><sup>min</sup><sub>min</sub></td><td>−[1 1]</td><td>−[0.1 0.2]</td><td>−[0.1 0.2]</td></tr></tbody></table>

when $\mathbf{u}_c(t) = \mathbf{x}(0)$ is chosen as the initial conditions for the controller dynamics as mentioned in Sec. 3.3.

In Scenario 1 (Fig. 4), small limits of inputs are considered. Profiles of the input and output variables with the MPC and CbI-MPC laws are described by the green dashed and blue continuous lines, respectively. The results illustrate the input constraint consideration in the CbI-MPC controller as well as the stability and the convergence to the references. Note that, since a relaxation is used to deal with the constraints, the constraints are not always respected. To improve the constraint satisfaction, we can reduce the time scale $\tau_z, \tau_\lambda, \tau_\mu$ as discussed in Sec. 3.3.

In Scenario 2 (Fig. 5), small limits of both inputs and states are considered. Comparing to Fig. 4, we can see that, besides the input constraint which is satisfied, the state constraint is also taken into account by the proposed controller.

Scenario 3 (Fig. 6) shows a clear advantage of the proposed CbI-MPC method with respect to the MPC method and the optimal CbI method developed in [20] (see Appendix C for the definition of the corresponding controller). Here, with this critical $\mathbf{u}_{max}$ value, the MPC optimization problem is not feasible. The optimal CbI controller tries to keep the inputs between their limits with remarkable oscillations during the first two seconds. Note also that this CbI controller does not handle

<figure><img src="image_4.png" /><figcaption>Figure 4: The profiles of the input and state vectors in Scenario 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 5: The profiles of the input and state vectors in Scenario 2.</figcaption></figure>

state constraints, which are completely violated in this approach.

## 5. Conclusion

This paper presents a novel control design to deal with system constraints using a Port-Hamiltonian formulation based on Model Predictive Control (MPC). The state
<figure><img src="image_6.png" /><figcaption>Figure 6: The profiles of the input and state vectors in Scenario 3.</figcaption></figure>

and input constraints are firstly taken into account by
formulating a MPC-type optimization problem. Then,
an open dynamical controller system is constructed
based on the primal-dual gradient method with an
additional nonlinear input. The controlled system and the
controller are finally coupled together using the Control by
Interconnection technique. The proposed control method
deals with both state and input constraints while explicitly
admitting the Hamiltonian as a Lyapunov function for the
closed-loop system. Moreover, a guideline to tune different
controller parameters is presented. The effectiveness of
the control design is illustrated in simulation through a
qualitative comparison with different control methods. As
future work, we aim at extending the proposed Cbl-MPC
method to more general systems where the input matrix
is not necessarily invertible. This can be realized by
replacing the plant with the passive augmented system
which includes the plant and an appropriate observer.

Acknowledgment

This work was supported in part by the Swiss National
Science Foundation.

References

[1] B. Maschke, A. V. D. Schaft, Port-controlled hamiltonian systems: modelling origins and system theoretic properties, in: Nonlinear Control Systems Design 1992, Elsevier, 1993, pp. 359–365.

[2] A. V. D. Schaft, Port-Hamiltonian Modeling for Control, Annual Review of Control, Robotics, and Autonomous Systems 3 (2020) 393–416.

[3] A. V. D. Schaft, D. Jeltsema, Port-hamiltonian systems theory: An introductory overview, Foundations and Trends in Systems and Control 1 (2014) 173–378.

[4] R.Ortega, A. V. D. Schaft, F. Castanos, A. Astolfi, Control by interconnection and standard passivity-based control of port-hamiltonian systems, IEEE Transactions on Automatic control 53 (2008) 2527–2542.

[5] P. Borja, R. Cisneros, R. Ortega, A constructive procedure for energy shaping of port—hamiltonian systems, Automatica 72 (2016) 230–234.

[6] R. Ortega, E. Garcia-Canseco, Interconnection and damping assignment passivity-based control: A survey, European Journal of control 10 (2004) 432–450.

[7] Z. Wang, A. Wei, G. Zong, X. Zhao, H. Li, Finite-time stabilization and control for a class of switched nonlinear port-controlled Hamiltonian systems subject to actuator saturation, Journal of the Franklin Institute 357 (2020) 11807–11829.

[8] T. Stegink, C. D. Persis, A. V. D. Schaft, A unifying energy-based approach to stability of power grids with market dynamics, IEEE Trans. on Automatic Control 62 (2017) 2612–2622.

[9] E. Benedito, D. del Puerto-Flores, A. Doria-Cerezo, J. Scherpen, Port-hamiltonian based optimal power flow algorithm for multi-terminal dc networks, Control Engineering Practice 83 (2019) 141–150.

[10] R. E. Kalman, When is a linear control system optimal?, in: Joint Automatic Control Conference, 1, 1963, pp. 1–15.

[11] A. V. D. Schaft, L2-gain and passivity techniques in nonlinear control, volume 2, Springer, 2000.

[12] N. Vu, L. Lefèvre, A connection between optimal control and ida-pbc design, 6th IFAC Workshop on Lagrangian and Hamiltonian Methods for Non Linear Control, Valparaiso, Chile, May 1-4 (2018) 1017–1022.

[13] Y. Wu, B. Hamroun, Y. Le-Gorrec, B. Maschke, Reduced order LQG control design for infinite dimensional Port-Hamiltonian systems, IEEE Transactions on Automatic Control 66 (2021) 865–871.

[14] L. Kölsch, P. Jané-Soneira, F. Strehle, S. Hohmann, Optimal control of port-hamiltonian systems: A time-continuous learning approach, Automatica 130 (2021) 109725.

[15] L. Gao, B. Shi, M. Kleeberger, J. Fottner, Optimal control of the hydraulic actuated boom system based on port-hamiltonian formulation, 12th International Fluid Power Conference 1 (2020) 489–498.

[16] D. Q. Mayne, Model predictive control: Recent developments and future promise, Automatica 50 (2014) 2967–2986.

[17] P. Falugi, Model predictive control: a passive scheme, 19th IFAC World Congress 47 (2014) 1017–1022.

[18] H. C. Pangborn, J. P. Koeln, A. G. Alleyne, Passivity and decentralized MPC of switched graph-based power flow systems, Annual American Control Conference (2018) 198–203.

[19] K. Yoshida, M. Inoue, T. Hatanaka, Instant MPC for linear systems and dissipativity-based stability analysis, IEEE Control. Syst. Lett. 3 (2019) 811–816.

[20] T. Stegink, C. D. Persis, A. V. D. Schaft, Port-hamiltonian formulation of the gradient method applied to smart grids, IFAC-PapersOnLine 48 (2015) 13–18.

[21] K. Arrow, L. Hurwicz, H. Uzawa, Studies in Linear and Non-Linear Programming, Cambridge Univ. Press, 1958.

[22] S. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.

[23] J. Rawlings, D. Mayne, Model Predictive Control: Theory and Design, Nob Hill Publishing, 2009.

[24] P. Kotyczka, L. Lefèvre, Discrete-time port-Hamiltonian systems: A definition based on symplectic integration, Systems & Control Letters 133 (2019).

[25] A. Venkatraman, A. V. D. Schaft, Full-order observer design for a class of port-hamiltonian systems, Automatica 46 (2010) 555–561.

[26] B. Vincent, N. Hudon, L. Lefèvre, D. Dochain, Port-hamiltonian observer design for plasma profile estimation in
tokamaks, IFAC-PapersOnLine 49 (2016) 93–98.

[27] M. Pfeifer, Thesis: Automated model generation and observer design for interconnected systems: a port-Hamiltonian approach, Ph.D. thesis, Karlsruhe institute of technology, Karlsruhe, Germany, 2021.

[28] R. P. Aguilera, D. E. Quevedo, On stability and performance of finite control set MPC for power converters, 2011 Workshop on Predictive Control of Electrical Drives and Power Electronics (2011) 55–62.

## Appendix A. Primal-dual gradient convergence proof

This section proves the result presented in Proposition 1, inspired from the ideas of [20]. Let $\Omega$ denote the set of equilibrium points of the dynamics (5), and $\Omega_\mu$ denote the following set:

$$ \Omega_{\mu} = \left\{ (\mathbf{z}, \lambda, \mu) \mid \begin{array}{l} \{\mu \ge 0, \mathbf{g}(\mathbf{z}) = 0\} \\ \text{or } \{\mu = 0, \mathbf{g}(\mathbf{z}) < 0\} \end{array} \right\} \quad (A.1) $$

From (4), we can see that $\Omega \subset \Omega_\mu$. Consider an equilibrium point $(\mathbf{z}^*, \lambda^*, \mu^*) \in \Omega$. The state deviations $(\tilde{\mathbf{z}}(t), \tilde{\lambda}(t), \tilde{\mu}(t))$ are defined as:

$$ (\tilde{\mathbf{z}}(t), \tilde{\lambda}(t), \tilde{\mu}(t)) = (\mathbf{z}(t), \lambda(t), \mu(t)) - (\mathbf{z}^*, \lambda^*, \mu^*). \quad (A.2) $$

From (5) and (A.2), the deviation dynamics are derived as:

$$ \left\{ \begin{array}{l} \left[ \begin{array}{c} \tau_z \dot{\tilde{z}} \\ \tau_\lambda \dot{\tilde{\lambda}} \end{array} \right] = \begin{bmatrix} \mathbf{0} & -\mathbf{A}_z^\top \\ \mathbf{A}_z & \mathbf{0} \end{bmatrix} \left[ \begin{array}{c} \tilde{\mathbf{z}} \\ \tilde{\lambda} \end{array} \right] + \\ \qquad -\nabla f(\mathbf{z}) + \nabla f(\mathbf{z}^*) + \nabla \mathbf{g}^\top(\mathbf{z}) \mu - \nabla \mathbf{g}^\top(\mathbf{z}^*) \mu^* \\ \\ \tau_\mu \dot{\tilde{\mu}} = [\mathbf{g}(\mathbf{z})]_\mu^+ . \end{array} \right. \quad (A.3) $$

The corresponding shifted Hamiltonian with respect to the equilibrium point is chosen as:

$$ \tilde{H}_r(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}) = \frac{1}{2} \tilde{\mathbf{z}}^\top(t) \tau_z \tilde{\mathbf{z}}(t) + \frac{1}{2} \tilde{\lambda}^\top(t) \tau_\lambda \tilde{\lambda}(t) \\ + \frac{1}{2} \tilde{\mu}^\top(t) \tau_\mu \tilde{\mu}(t). \quad (A.4) $$

Firstly, we admit the following inequalities (see the proof in [20]):

$$ \tilde{\mu}^{\top} [\mathbf{g}(\mathbf{z})]_{\mu}^{+} \le \tilde{\mu}^{\top} \mathbf{g}(\mathbf{z}), \quad \text{from (6, A.1)}, \quad (A.5) $$

$$ \mathbf{g}(\mathbf{z}) \le \mathbf{g}(\mathbf{z}^*) + \tilde{\mathbf{z}}^\top \nabla^\top \mathbf{g}(\mathbf{z}), \quad \mathbf{g}(\mathbf{z}) \text{ is convex}, \quad (A.6) $$

$$ \mathbf{g}(\mathbf{z}) \ge \mathbf{g}(\mathbf{z}^*) + \tilde{\mathbf{z}}^\top \nabla^\top \mathbf{g}(\mathbf{z}^*), \quad \mathbf{g}(\mathbf{z}) \text{ is convex}, \quad (A.7) $$

$$ \tilde{\mu}^{\top} \mathbf{g}(\mathbf{z}^*) \le 0, \quad \text{from (4, 6).} \quad (A.8) $$

From (A.2), (A.3) and (A.4), we obtain:

$$ \begin{aligned} \dot{\tilde{H}}_r(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}) = & -\tilde{\mathbf{z}}^\top[\nabla f(\mathbf{z}) - \nabla f(\mathbf{z}^*)] \\ & -\tilde{\mathbf{z}}^\top[\nabla^\top \mathbf{g}(\mathbf{z}) - \nabla^\top \mathbf{g}(\mathbf{z}^*)] \mu^* \\ & -\tilde{\mathbf{z}}^\top \nabla_z^\top \mathbf{g}(\mathbf{z}) \tilde{\mu} + \tilde{\mu}^\top[\mathbf{g}(\mathbf{z})]_\mu^+. \end{aligned} \quad (A.9) $$

We also have the following inequalities:

$$ \left\{ \begin{array}{ll} -\tilde{\mathbf{z}}^{\top} [\nabla f(\mathbf{z}) - \nabla f(\mathbf{z}^*)] & \le 0, \quad \text{from Assump. 2,} \\ -\tilde{\mathbf{z}}^{\top} [\nabla^\top \mathbf{g}(\mathbf{z}) - \nabla^\top \mathbf{g}(\mathbf{z}^*)] \mu^* & \le 0, \quad \text{from (A.1),} \\ -\tilde{\mathbf{z}}^{\top} \nabla_z^\top \mathbf{g}(\mathbf{z}) \tilde{\mu} + \tilde{\mu}^\top [\mathbf{g}(\mathbf{z})]_\mu^+ & \le 0, \quad \text{from (A.5-A.8).} \end{array} \right. \quad (A.10) $$

From (A.9) and (A.10), we obtain:

$$ \dot{\tilde{H}}_r(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}) \le 0, \quad \forall (\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}). \quad (A.11) $$

Let $\mathbb{M} = \{(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu})\}$ denote the largest invariant set of the system (A.3) such that $\dot{\tilde{H}}_r(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}) = 0, \forall (\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}) \in \mathbb{M}$. From Assumption 2, (A.9) and (A.10), we derive that $\forall (\tilde{\mathbf{z}}(t), \tilde{\lambda}(t), \tilde{\mu}(t)) \in \mathbb{M}, \tilde{\mathbf{z}}(t) = \mathbf{0}, \text{ or } \mathbf{z}(t) = \mathbf{z}^*$. Let $\mathbb{M}_r$ denote the set of $(\mathbf{z}, \lambda, \mu)$ such that $(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu}) \in \mathbb{M}$. From (A.4) and (A.11), by LaSalle's invariance principle we may conclude that $(\tilde{\mathbf{z}}, \tilde{\lambda}, \tilde{\mu})$ converges to $\mathbb{M}$, i.e.:

$$ \begin{array}{l} (\tilde{\mathbf{z}}(t), \tilde{\lambda}(t), \tilde{\mu}(t)) \xrightarrow{t \to \infty} \mathbb{M}, \\ \text{or} \\ (\mathbf{z}(t), \lambda(t), \mu(t)) \xrightarrow{t \to \infty} \mathbb{M}_r. \end{array} \quad (A.12) $$

When $(\mathbf{z}, \lambda, \mu) \in \mathbb{M}_r$, we consider the dynamics of $\mu(t)$ in (5), that is $\dot{\mu}(t) = [\mathbf{g}(\mathbf{z}^*)]_\mu^+$. If $\mathbf{g}(\mathbf{z}^*) = \mathbf{0}$ and $\mu(t) \ge \mathbf{0}$, $\dot{\mu}(t) = \mathbf{0}$. If $\mathbf{g}(\mathbf{z}^*) < \mathbf{0}$ and $\mu(t) = \mathbf{0}$, $\dot{\mu}(t) = \mathbf{g}(\mathbf{z}^*) < \mathbf{0}$. Therefore, it is easy to see that when $\mathbf{g}(\mathbf{z}^*) < \mathbf{0}$, $\mu(t) \xrightarrow{t \to \infty} \mathbf{0}$, e.g.:

$$ (\mathbf{z}^*, \lambda(t), \mu(t)) \xrightarrow{t \to \infty} \Omega_\mu \cap \mathbb{M}. \quad (A.13) $$

When $(\mathbf{z}, \lambda, \mu) \in (\Omega_\mu \cap \mathbb{M}), \dot{\mathbf{z}}(t) = \mathbf{0}, \dot{\lambda}(t) = \mathbf{0}$ and $\dot{\mu}(t) = \mathbf{0}$, and thus:

$$ (\mathbf{z}(t), \lambda(t), \mu(t)) \in \Omega. \quad (A.14) $$

From (A.12), (A.13) and (A.14), we conclude that:

$$ (\mathbf{z}(t), \lambda(t), \mu(t)) \xrightarrow{t \to \infty} \Omega. $$

## Appendix B. Discrete-time Port-Hamiltonian system

This section briefly recalls a definition of discrete-time port-Hamiltonian systems based on the symplectic integration presented in [24]. Using the collocation method, we define $s \in \mathbb{N}$ collocation points $\{t_1^k, \dots, t_s^k\}$ over a time step $[k\Delta t, (k+1)\Delta t]$ such that:

$$ k\Delta t < t_i^k < t_{i+1}^k < (k+1)\Delta t, \text{ with } i \in \{1, \dots, s-1\}. $$

Let $\mathbf{x}_i^k \in \mathbb{R}^n$, $\mathbf{u}_i^k \in \mathbb{R}^m$ represent the state and input vectors at the collocation point $t_i^k$ with $i \in \{1, \dots, s\}$. According to [24], the discrete-time state $\mathbf{x}(k+1|t)$ of the port-Hamiltonian system (1) is determined as:

$$ \left\{ \begin{array}{l} \mathbf{x}_i^k = \mathbf{x}(k|t) - \Delta t \sum_{j=1}^{s} a_{ij} [\mathbf{J}(\mathbf{x}_j^k) - \mathbf{R}(\mathbf{x}_j^k)] \nabla H_x(\mathbf{x}_j^k) \\ \qquad + a_{ij} \mathbf{G}(\mathbf{x}_j^k) \mathbf{u}_j^k, \\ \\ \mathbf{x}(k+1|t) = \mathbf{x}(k|t) - \Delta t \sum_{j=1}^{s} b_j [\mathbf{J}(\mathbf{x}_j^k) - \mathbf{R}(\mathbf{x}_j^k)] \nabla H_x(\mathbf{x}_j^k) \\ \qquad + b_j \mathbf{G}(\mathbf{x}_j^k) \mathbf{u}_j^k, \end{array} \right. \quad (B.1) $$

where $a_{ij}, b_j \in \mathbb{R}$ are constants computed from the collocation functions and points, and $i, j \in \{1, \dots, s\}$.
The system (B.1) is linear if there exist constant matrices $\mathbf{A} \in \mathbb{R}^{n \times n}$ and $\mathbf{B} \in \mathbb{R}^{n \times m}$ such that:

$$ 
\left\{ 
\begin{aligned}
\mathbf{A}\mathbf{x}(k|t) &= \mathbf{x}(k|t) - \Delta t \sum_{j=1}^{s} b_j [\mathbf{J}(\mathbf{x}_j^k) - \mathbf{R}(\mathbf{x}_j^k)] \nabla H_x(\mathbf{x}_j^k) \\
\mathbf{B}\mathbf{u}(k|t) &= \Delta t \sum_{j=1}^{s} b_j \mathbf{G}(\mathbf{x}_j^k) \mathbf{u}_j^k
\end{aligned}
\right.
\qquad (B.2)
$$

## Appendix C. Optimal Control by Interconnection

This subsection recalls the optimal CbI method previously presented in [20]. The control design is based on the following optimization problem:

$$ \min_{\mathbf{u}} f(\mathbf{u}) \qquad (C.1a) $$

$$ \text{s.t.} \quad [\mathbf{J}_x - \mathbf{R}_x] \nabla H_x + \mathbf{G}_x \mathbf{u} = \mathbf{0}, \qquad (C.1b) $$

$$ \mathbf{g}(\mathbf{u}) \le \mathbf{0}, \qquad (C.1c) $$

where $\mathbf{u}$ is the input of the controlled system. The cost function $f(\mathbf{u}) \in \mathbb{R}$ is a convex function, derived from $f(\mathbf{z})$ in the optimization problem (13) with $\mathbf{Q}_x = \mathbf{0}$ and $\mathbf{Q}_f = \mathbf{0}$. The equality constraint (C.1b) is the equilibrium condition, and the inequality constraint (C.1c) considers the same limits of the input $\mathbf{u}$ as in (13c), and $\mathbf{g}(\mathbf{0}) \le \mathbf{0}$, as presented in Assumption 2.

Using primal-dual gradient method presented in Subsection 2.2 with additional conjugate input and output, the controller is given as:

$$ 
\left\{ 
\begin{aligned}
\tau_u \dot{\mathbf{u}}(t) &= -\nabla f(\mathbf{u}) - \mathbf{G}_x^\top \lambda(t) - \nabla \mathbf{g}^\top(\mathbf{u}) \mu(t), \\
\tau_\lambda \dot{\lambda}(t) &= [\mathbf{J}_x - \mathbf{R}_x] \nabla H_x + \mathbf{G}_x \mathbf{u}(t), \\
\tau_\mu \dot{\mu}(t) &= [\mathbf{g}(\mathbf{u})]_\mu^+, \\
\mathbf{y}_c(t) &= \mathbf{u}(t).
\end{aligned}
\right.
\qquad (C.2)
$$

Note that this control considers the input constraint but can not deal with the state constraint.