Dorothée Normand-Cyrot
email: dorothee.normand-cyrot@centralesupelec.fr

Salvatore Monaco
email: salvatore.monaco@uniroma1.it

Mattia Mattioni
email: mattia.mattioni@uniroma1.it

Alessio Moreschini
email: a.moreschini@imperial.ac.uk

Passivity techniques and Hamiltonian structures in discrete time

Keywords: Discrete-time systems, Port-Hamiltonian structures, Passivity based control, Energy-based control

Passivity techniques and Hamiltonian structures in discrete time

Dorothée Normand-Cyrot 1 , Salvatore Monaco 2 , Mattia Mattioni 2 , and Alessio Moreschini 3

Introduction

Energy-based modeling and control are fundamental evergreen concepts that are extensively investigated because of their sustained impact towards the newest methodologies and technologies (see the textbooks or survey-oriented contributions [START_REF] Willems | Dissipative dynamical systems part i: General theory[END_REF]1,[START_REF] Sepulchre | 50 years of dissipativity theory, part ii [about this issue[END_REF][START_REF] Ortega | Putting energy back in control[END_REF][START_REF] Van Der | L2-gain and passivity techniques in nonlinear control[END_REF][START_REF] Sepulchre | Constructive nonlinear control[END_REF][START_REF] Brogliato | Dissipative Systems Analysis and Control: Theory and Applications[END_REF]). The underlying idea consists in deducing a system representation explicitly catching the energy features of the dynamics through dissipative or conservative components. Then, energy based control strategies can be designed. In this framework, the class of port-Hamiltonian systems is paradigmatic [START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF][START_REF] Ortega | PID Passivity-Based Control of Nonlinear Systems with Applications[END_REF][START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF]. From a theoretical point of view, most of the literature is devoted to the continuous-time setting in spite of a pervasive interest in computer-oriented applications and thus in digital structures at large. The general obstacles to perform an equivalent analysis in a discrete-time framework are well-known and related to the difficulty to describe the geometric structure behind the state and output discrete evolutions [START_REF] Monaco | Advanced tools for nonlinear sampleddata systems' analysis and control[END_REF]. In this sense, the definition of dissipation itself is an open question [START_REF] Byrnes | Discrete-time lossless systems, feedback equivalence and passivity[END_REF][START_REF] Byrnes | Losslessness, feedback equivalence, and the global stabilization of discrete-time nonlinear systems[END_REF][START_REF] Navarro-López | Several dissipativity and passivity implications in the linear discrete-time setting[END_REF][START_REF] Monaco | Nonlinear representations and passivity conditions in discrete time[END_REF]. As a consequence, a shared definition of discrete-time port-Hamiltonian structures has been missing in spite of a variety of approaches often inspired by discretization schemes aimed to preserve the energy properties and/or power-balance exchanges (to cite a few see [START_REF] Itoh | Hamiltonian-conserving discrete canonical equations based on variational difference quotients[END_REF][START_REF] Gonzalez | Time integration and discrete Hamiltonian systems[END_REF][START_REF] Stramigioli | Sampled data systems passivity and discrete port-Hamiltonian systems[END_REF][START_REF] Laila | Discrete-time IDA-PBC design for separable Hamiltonian systems[END_REF][START_REF] Laila | Construction of discrete-time models for portcontrolled Hamiltonian systems with applications[END_REF][START_REF] Talasila | Discrete port-Hamiltonian systems[END_REF][START_REF] Sümer | A direct discrete-time IDA-PBC design method for a class of underactuated Hamiltonian systems[END_REF][START_REF] Sümer | A direct discrete-time IDA-PBC design method for a class of underactuated Hamiltonian systems[END_REF][START_REF] Šešlija | Port-Hamiltonian systems on discrete manifolds[END_REF][START_REF] Seslija | Explicit simplicial discretization of distributed-parameter port-Hamiltonian systems[END_REF][START_REF] Castaños | Discrete-time models for implicit port-Hamiltonian systems[END_REF][START_REF] Yalçin | Discrete-time modeling of Hamiltonian systems[END_REF][START_REF] Kotyczka | Discrete port-Hamiltonian formulation and numerical approximation for systems of two conservation laws[END_REF][START_REF] Aoues | Hamiltonian systems discrete-time approximation: Losslessness, passivity and composability[END_REF][START_REF] Celledoni | Energy-preserving and passivity-consistent numerical discretization of port-Hamiltonian systems[END_REF]).

The present work presents recent advances in the definition of dissipative notions for discrete-time systems when adopting a Differential Difference Representation (DDR) and their impact for characterizing port-Hamiltonian structures. The idea behind such a representation consists in separating discrete-time dynamics into two components: a difference equation that describes the control-free evolution and a differential equation that models the variation of the dynamics with respect to the control variable. This hybrid structure is adapted to cope with the intrinsic nonlinearity (that becomes more and more complex through iterative composition along successive time steps) in both the state and input variables when representing discrete-time dynamics in the map form. Splitting the free evolution from the controlled part leads naturally to a description of a discrete-time dynamics through two coupled difference and differential equations. An exponential representation of the discrete-time flow results through the integration of the differential controlled part. This is useful for further composing the state evolutions along successive time steps; the composition of exponential flows replace the composition of nonlinear functions. Analogously, when considering an output function, one gets an exponential form representation of the Volterra series expansions characterizing the input-to-output evolutions. Recalling, if necessary, that typical numerical burdens have to be faced when considering both purely discrete systems or systems issued from sampling due to their intrinsic nonlinearities, these exponential forms provide efficient computational tools. Accordingly, the combinatorial properties of the series expansions that describe the solutions can be exploited to qualify suitable polynomial approximations of increasing order. Section 2 recalls these developments as useful prolegomena to the paper while details can be found in [START_REF] Monaco | A lie exponential formula for the nonlinear discrete time functional expansion[END_REF][START_REF] Monaco | Nonlinear systems in discrete time[END_REF][START_REF] Monaco | A unified representation for nonlinear discrete-time and sampled dynamics[END_REF][START_REF] Monaco | From chronological calculus to exponential representations of continuous and discrete-time dynamics: A lie-algebraic approach[END_REF].

In Section 3, the aforementioned representation is used to characterize passivity. This notion, and more generally the one of dissipativity, relies upon energy exchanges and the way the system interacts with the environment. Roughly speaking, a system is passive when the internally stored energy does not exceed the one externally supplied [START_REF] Willems | Dissipative dynamical systems part i: General theory[END_REF][START_REF] Willems | Dissipative dynamical systems part ii: Linear systems with quadratic supply rates[END_REF]. This property is caught by characterizing the variation of a particular function along the system's trajectories [START_REF] Van Der | L2-gain and passivity techniques in nonlinear control[END_REF][START_REF] Sepulchre | Constructive nonlinear control[END_REF][START_REF] Ortega | Putting energy back in control[END_REF]. Such a function, referred to as the storage function, generally represents the energy and is strictly linked to Lyapunov and/or Hamiltonian functions. With this in mind and for a fixed storage quantity, the DDR form we propose immediately allows us to characterize average passivity [START_REF] Monaco | Nonlinear representations and passivity conditions in discrete time[END_REF]; a novel notion of passivity in discrete time. This is done by isolating the control dependent part and thus defining in a very natural way the corresponding average passive output that is, in turn, a conjugate quantity whose product with the control variable is a power unit. The notion of average dissipation can be then exploited for control purposes by extending the usual concept of negative output damping to cope with stabilization at the origin. This approach can be further generalized to describe Passivity-Based Control (PBC) at large for nonlinear discrete-time systems towards the so-called second generation of passivity-based control, aimed at managing the system energy to satisfy control specifications. Those techniques endow Energy Balancing (EB) and Interconnection and Damping Assignment (IDA) to deal with more complex systems including networked or cascade dynamics. The foundations of average passivity-based techniques are recalled in Section 3, while more detailed studies are in [START_REF] Monaco | From passivity under sampling to a new discrete-time passivity concept[END_REF][START_REF] Monaco | Nonlinear average passivity and stabilizing controllers in discrete time[END_REF].

As already mentioned, port-Hamiltonian dynamics are of pervasive interest due to their mathematical structure and their foundations in physics (to cite a few [START_REF] Maschke | Port-controlled Hamiltonian systems: modelling origins and system theoretic properties[END_REF][START_REF] Maschke | Port-controlled Hamiltonian systems: modelling origins and systemtheoretic properties[END_REF][START_REF] Maschke | An intrinsic Hamiltonian formulation of network dynamics: Non-standard poisson structures and gyrators[END_REF][START_REF] Ortega | Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian systems theory: An introductory overview[END_REF][START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF][START_REF] Rashad | Twenty years of distributed port-Hamiltonian systems: a literature review[END_REF]). Hereinafter, a novel state-space representation of discrete-time port-Hamiltonian structures that represents a breakthrough in the literature is naturally deduced from the DDR form. The so defined forms are endowed with average passivity properties that validate the proposed choice of conjugate output. Further on, the fundamental characteristics of port-Hamiltonian structures as the qualifying closeness property under power-preserving interconnection [START_REF] Moreschini | Interconnection through u-average passivity in discrete time[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF], are validated. As a consequence, these forms are efficient for the design of average passivity based control strategies for complex and networked discrete-time dynamics as illustrated in [START_REF] Moreschini | Stabilization of discrete port-Hamiltonian dynamics via interconnection and damping assignment[END_REF][START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF][START_REF] Mattioni | Quaternion-based attitude stabilization via discrete-time IDA-PBC[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF]. Port-Hamiltonian structures are discussed in Section 4 while dedicated studies are in [START_REF] Moreschini | Discrete port-controlled Hamiltonian dynamics and average passivation[END_REF][START_REF] Moreschini | Stabilization of discrete port-Hamiltonian dynamics via interconnection and damping assignment[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF].

Notations

Throughout the paper all the functions and vector fields defining the dynamics are assumed smooth and complete over the respective definition spaces. The sets R, N and Z denote the set of real, natural numbers including 0 and integers respectively. For any vector v ∈ R n , |v| and v ⊤ define the norm and transpose of v respectively. I d denotes the identity function on the definition space while I denotes the identity operator and the identity matrix when related to a linear operator. The symbols " > 0" and " < 0" denote positive and negative definite functions (or matrices), respectively. The symbol • denotes the composition of two functions or operators, depending on the context. Given a realvalued function V : R n → R assumed differentiable, setting ∇ = col{ ∂ ∂xi } i=1,n , ∇V (x) represents the gradient column-vector. Given a vector-valued function F (x) = col(F 1 (x), . . . , F n (x)), the operator

J x [F ](x) = ∂F (x) ∂x = ∂F i ∂x j (x) i,j=1,n
denotes the Jacobian of the function F evaluated at x. Given a smooth vector field over R n , e f ) denotes the exponential Lie operator e f := I + i≥1

L i f i! with L f = n i=1 f i (x) ∂
∂xi and indicating by L i f the operator at power i with respect to the usual composition of vector fields (for a linear vector field, the exponential Lie operator recovers the exponential of the matrix representing the operator). For any smooth function h : R n → R, one verifies e f h(x) = h(e f (x)) = e f h| x where | x denotes the evaluation of the function at x. Given two vector fields f, g, their Lie bracket is again a vector field denoted by

ad f g = [f, g] = (L f •L g I d -L g •L f )I d and, iteratively for i ≥ 1, ad i f g = ad i-1 f • ad f g with ad 1 f g = ad f g. A function R(x, u) is said in O(u p ) for p ≥ 1 if,
whenever it is defined, it can be written as R(x, u) = u p-1 R(x, u) and there exist a function θ ∈ K ∞ and u

⋆ > 0 such that ∀u ≤ u ⋆ , | R(x, u)| ≤ θ(u).
2 Basic on discrete-time dynamics

Differential/Difference Representation

A single-input nonlinear discrete-time dynamics over R n is usually represented by a function F : R n × R → R n , smooth in both the state and input variables

x k+1 = F (x k , u k ) = x k + F 0 (x k ) + g(x k , u k )u k .
(1)

For any pair of state and input variables (x k , u k ) fixed at time instant k ∈ N, x k+1 denotes the state reached at time k + 1 from x k , under the action of the control u k . For convenience that will be clear later on and without loss of generality, the free evolution is decomposed as F (x, 0) = x + F 0 (x) while g(x, u)u (with g(x, 0) ̸ = 0) represents the control dependent part of the dynamics.

In [START_REF] Monaco | A unified representation for nonlinear discrete-time and sampled dynamics[END_REF], we proposed to represent discrete dynamics as two coupled difference and differential equations (DDR). More in detail, for all time step k and under mild conditions (e.g., submersivity of F (x, u)), (1) can be represented as

x + (0) = x + F 0 (x) (2a) dx + (u) du = G(x + (u), u) (2b)
when denoting (x, u) = (x k , u k ) ∈ R n × R, any pair of state and input variables at generic time instant k, x + (u) = x + (u k ) = x k+1 , the state reached at time k+1 starting from x k under the action of u k and x + (0), the state reached along free evolution u k = 0. In doing so, we underline that x + (u) is viewed as a curve in R n , parameterized by the control variable u ∈ R. G(•, u) is a properly defined vector field on R n , parameterized by u, satisfying G(F (x, u), u) = ∂F ∂u (x, u). With this in mind, the difference equation (2a) describes the jump of the state evolving in free evolution while the differential equation (2b) models the rate of change of the state dynamics with respect to control variation.

Remark 1. Provided F (x, 0) be invertible, the vector field G(x, u) can be defined for u sufficiently small as

G(x, u) := ∂F (x, u) ∂u x=F -1 (x,u) , (3)
and expanded, around u = 0, as

G(x, u) = G 1 (x) + i≥1 u i i! G i+1 (x); G 1 (x) = G(x, 0). ( 4
)
The vector fields G i in the power expansion (4) define a family of canonical vector fields associated to the dynamics [START_REF] Aoues | Hamiltonian systems discrete-time approximation: Losslessness, passivity and composability[END_REF]. The (G i ) ′ s have been proved useful in the study of analysis and design problems linked to the geometry of the evolutions in the state space (e.g. [START_REF] Monaco | A lie exponential formula for the nonlinear discrete time functional expansion[END_REF][START_REF] Monaco | Nonlinear systems in discrete time[END_REF][START_REF] Monaco | A unified representation for nonlinear discrete-time and sampled dynamics[END_REF][START_REF] Califano | Non-linear non-interacting control with stability in discrete-time: a geometric framework[END_REF][START_REF] Monaco | Normal forms and approximated feedback linearization in discrete time[END_REF]).

Remark 2. The DDR form can be generalized to multi-input dynamics by modeling the rate of change of the dynamics under the action of each control. Setting u = (u 1 , . . . , u m ) ⊤ , one replaces (2b) by the set of partial derivative equations ∂x + (u) ∂u j = G j (x + (u), u) for j = (1, . . . , m) with x + (0) = F (x, 0) [START_REF] Byrnes | Discrete-time lossless systems, feedback equivalence and passivity[END_REF] and G j (x, u) satisfying the condition G j (F (x, u), u) := ∂F (x,u) ∂u j .

Some useful manipulations

To better grasp the DDR forms, some useful concepts are specified below. Particular classes can be discussed besides the linear one that corresponds to set in (2), F 0 (x) = Ax with constant control vector G(x, u) = B and matrices (A, B) of suitable dimensions. A peculiar class is represented by the input-affine one. 2), the dynamics in map form (1) preserves a nonlinear dependence on u which admits the exponential representation

Proposition 1. -Input-affine DDR -Assuming G(x, u) = G 1 (x) in (
F (x, u) = e uG1 (x) x+F0(x) (6)
where e uG1 denotes the flow associated with the solution to the differential equation (2b) when

G(x, u) = G 1 (x).
This result follows from the fact that integrating (2b) from 0 and u with initial condition x + (0) = x + F 0 (x), one recovers a dynamics in map form (1); i.e.

x

+ (u) = x + F 0 (x) + u 0 G(x + (s), s)ds = x + F 0 (x) + g(x, u)u. (7)
When assuming G(x, u) = G 1 (x), this integral form admits the simple exponential representation described in Proposition 1 that is generalized in the sequel to a control dependent vector field G(•, u).

Remark 3. We note for completeness that given an input-affine dynamics in map form, g(x, u) = g(x) in (1), the associated control vector field G(x, u) in ( 4) is no more linear in u and satisfies the algebraic constraint G(F (x, u), u) = g(x).

Easy manipulations show how DDR forms are transformed under coordinates change and feedback.

Lemma 1. Let the coordinates change z = T (x) defined by the diffeomorphism T : R n → R n , then the DDR dynamics ( 2) is transformed into

z + (0) = z + F0 (z) dz + (u) du = Ḡ(z + (u), u) (8)
with

z + F0 (z) = T (x + F 0 (x)) x=T -1 (z) Ḡ(z, u) = Ad T G(z, u)
where Ad T G(•, u) indicates the transport of the vector field G(•, u) along T ; i.e.

Ad T G(z, u) = ∂T ∂x G(•, u) x=T -1 (z) = L G(•,u) T (x) x=T -1 (z) . (9)
Lemma 2. Let the state feedback u(x, v) = α(x) + v with α : R n → R smooth and external control v ∈ R, then the DDR dynamics ( 2) is transformed into

x + α (0) = x + F α (x) dx + α (v) dv = G α (x + α (v), v) (10)
where

x + α (v) = F (x, α(x) + v) = x + (α(x) + v) x + α (0) = F (x, α(x)) = x + F α (x) = x + (α(x)) G α (x, v) = G(x, α(x) + v).
The proof works out showing that

dx + α (v) dv = ∂F (x, u) ∂u u=α(x)+v ∂(α(x) + v) ∂v = G(x + α (v), α(x) + v) = G α (x + α (v), v).

Input-to-state and input-to-output trajectories

In our context, a discrete-time system Σ d (h) is given by a state dynamics in the DDR form (2) and a smooth output function h : R n → R. The generalization of Proposition 1 to a non input-affine form is reported and further extended to describe input-output evolutions over one or several time-steps.

Theorem 1. Given Σ d (h) and a pair of state and input variables (x k , u k ) at generic time instant k, one gets at time k + 1

x k+1 = x + (u k ) = e u k G(•,u k ) I d | x k +F0(x k ) = F (x k , u k ) h(x k+1 ) = h(x + (u k )) = e u k G(•,u k ) h| x k +F0(x k ) = h(F (x k , u k ))
where the series exponent is a Lie series G(•, u) ∈ Lie{G 1 , . . . , G p , . . . } that can be described through its expansion in powers of u. For the first term one reports

G(•, u) = G 1 + u 2 G 2 + u 2 3! (G 3 + [G 1 , G 2 ]) + O(u 3 ).
The proof, developed in [START_REF] Monaco | A lie exponential formula for the nonlinear discrete time functional expansion[END_REF][START_REF] Monaco | From chronological calculus to exponential representations of continuous and discrete-time dynamics: A lie-algebraic approach[END_REF], follows from the Lie properties endowed by the flow characterizing the solutions to nonlinear u-dependent ordinary differential equations of the form (2b), when expressed through a chronological exponential.

A major property of these exponential representations of discrete flows is to be easily composed along successive time steps. For, it is instrumental to define, according to [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF], the transport of any vector field G i along the free evolution (Id + F 0 ) so getting for j ≥ 0, the family of transported vector fields (G j i ) ′ s with

G j+1 i (x) = Ad I d +F0 G j i (x) = L G j i (x + F 0 (x)) (I d +F0) -1 (x)
and iteratively

G j+1 i (x) = Ad j+1 I d +F0 G i (x) with G 0 i = G i .
The following result holds when denoting by G j (•, u); j ≥ 0, the series exponent

G j (•, u) ∈ Lie{G j 1 , . . . , G j p , . . . }
defined as in Theorem 1 with respect to the transported vector fields (G j i ) ′ s .

Theorem 2. Given Σ d (h), initial state value x 0 ∈ R n and input sequences {u 0 , . . . , u k }, k ∈ N, then at generic time instant k + 1 > 0, one gets

x k+1 = x + (u k , . . . , u 0 ) = e u0G k (•,u0) • • • • • e u k G 0 (•,u k ) | (I d +F0) k+1 (x0) y k+1 = y + (u k , . . . , u 0 ) = e u0G k (•,u0) • • • • • e u k G 0 (•,u k ) h| (I d +F0) k+1 (x0) .
According to these results it becomes clear that the vector fields (G i ) ′ s and their transport along the drift term enter in a differential geometric characterization of the structure of the accessible states. In fact, accessibility can be reported to the properties of the orbit of the associated Lie groups (see [START_REF] Fliess | A lie-theoretic approach to nonlinear discrete-time controllability via ritt's formal differential groups[END_REF][START_REF] Jakubczyk | Automatique théorique. orbites de pseudo-groupes de difféomorphismes et commandabilité des systèmes non linéaires en temps discret[END_REF] for further details). Further on, invariance can be characterized in terms of the Lie algebra Lie{G j 1 , . . . , G j p , . . . } as discussed in [START_REF] Monaco | Invariant distributions for nonlinear discrete-time systems[END_REF]. The same holds regarding the properties of controlled invariance, feeback linearization [START_REF] Monaco | Normal forms and approximated feedback linearization in discrete time[END_REF] or decoupling [START_REF] Califano | Non-linear non-interacting control with stability in discrete-time: a geometric framework[END_REF], up to characterize the corresponding control solutions.

Concluding comments

In this section, an alternative description to the usual map form is proposed for discrete-time dynamics. The free evolution that generates a purely discrete evolution defines the initial condition of the differential equation modeling the control action. This suggests a modeling approach that would separately identify the free evolution as a map F 0 and the variation of the dynamics with respect to u as a vector field G(•, u), that may depend on u. In doing so, a family of control vector fields, the (G i ) ′ s , that enters in the characterization of the structural and control properties of the dynamics, is defined. Accordingly, it becomes possible to combine a more visual geometric approach with a more computational algebraic one to provide intriguing relationships between the continuous-time and discrete-time settings that should be further investigated towards a unified vision. Further on, being the solution to an ordinary differential equation, the state dynamics can be rewritten in terms of the exponential operator characterizing the associated flow. It follows that the whole differential geometry apparatus behind this operator form can be used to describe the input-to-state and inputto-output behaviours along successive time steps. In particular, the composition of nonlinear maps is replaced with the composition of exponential forms that is more tractable in practice. This is well illustrated with the characterization of the Volterra series and the computation of its kernels in terms of the vector fields (G j i )

′ s and their Lie brackets. We refer to [START_REF] Monaco | On the realization of nonlinear discrete-time systems[END_REF][START_REF] Monaco | Finite volterra-series realizations and input-output approximations of non-linear discrete-time systems[END_REF] for further discussion along these lines in relation with realization problems.

Passivity techniques

The notion of average passivity has been introduced in [START_REF] Monaco | Nonlinear average passivity and stabilizing controllers in discrete time[END_REF] to weaken the necessary requirement of direct throughput when adopting the standard notion in a discrete-time framework. This definition is in fact directly inspired by the splitting of the state dynamics into free and control parts.

Average passivity

Denoting by Σ d (h) a discrete-time system described by the dynamics (1) (equivalently (2)) with output h : R n → R, the definition below is recalled from [START_REF] Monaco | Nonlinear average passivity and stabilizing controllers in discrete time[END_REF].

Definition 1 (Average passivity). Σ d (h) is said average passive if there exists a positive semi-definite function S : R n → R ≥0 (the storage function) such that for all (x, u) ∈ R n × R the following inequality holds

S(x + (u)) -S(x) ≤ u 0 h(x + (s))ds := uh av (x, u) (11)
with the average output defined as

h av (x, u) := 1 u u 0 h(x + (s))ds. ( 12
)
This definition is directly inspired by the DDR form of the dynamics that yields to rewrite the rate of change of the storage function S, between two successive states, in an integral form as

S(x + (u)) -S(x) = S(x + (0)) -S(x) + u 0 L G(•,s) S(x + (s))ds.
The so-defined average map h av (•, u) in [START_REF] Haier | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF], introduces a direct input-output link in such a way that average passivity with respect to h is in fact equivalent to usual passivity with respect to h av (•, u).

As an immediate consequence of average passivity, setting u = 0 in [START_REF] Gonzalez | Time integration and discrete Hamiltonian systems[END_REF], one verifies S(x + (0)) -S(x) ≤ 0, so concluding stability of any equilibrium x e ∈ R n when S qualifies as a Lyapunov function at x e (S(x e ) = 0 and S(x) > 0 for x ̸ = x e ) and asymptotic stability provided S(x + (0)) -S(x) < 0. More in general, average dissipativity can be defined making reference to the average output and a supply rate function s : R × R → R, so requiring that the dissipation inequality below be verified for all (x, u

) ∈ R n × R S(x + (u)) -S(x) ≤ s(u, h av (x, u)).
The notion of average passivity can be extended to average passivity from some nominal non-zero constant value ū ̸ = 0, [START_REF] Mattioni | Forwarding stabilization in discrete time[END_REF][START_REF] Moreschini | Interconnection through u-average passivity in discrete time[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF]. It is instrumental when discussing the action of a feedback law over passivity.

Definition 2 (Average passivity from ū). Σ d (h) is said average passive from a given ū with ū ∈ U ⊆ R, if there exists a positive semi-definite function S : R n → R ≥0 (the storage function) such that, for all

(x, u) ∈ R n × R S(x + ū (u)) -S(x) ≤ uh av ū (x, u) (13)
with u-average output from ū defined as

h av ū (x, u) := 1 u u 0 h(x + (ū + s))ds.
When ū = 0, one recovers average passivity since h av 0 (x, u) = h av (x, u).

Feedback stabilization and interconnection

On these bases, Passivity Based Control (PBC) techniques can be developed as discussed below. Two basic ingredients of control strategies exploiting passivity are revisited hereafter making reference to the average notion. Firstly, stabilization through negative output feedback and additional damping is specified. Secondly, closeness under power preserving input-output interconnection is discussed so enlarging the control design to interconnected dynamics. The next definition is instrumental.

Definition 3 (Zero-state detectability). Given Σ d (h) let x e ∈ R n be an equilibrium and Z ∈ R n be the largest invariant set contained in the set {x ∈

R n s.t. h(x + (0), 0) = 0}. Σ d (h) is said zero-state detectable (ZSD) if x = x e is
an asymptotically stable equilibrium conditionally to Z.

The following theorem characterizes the negative output damping feedback.

Theorem 3 (Negative average output feedback). Given Σ d (h) with equilibrium x e ∈ R n , (F 0 (x e ) = 0), assumed average passive with storage function S > 0, then the feedback u = α(x) solving the algebraic equality

u = -κh av (x, u) with κ > 0 ( 14
)
ensures asymptotic stability provided that Σ d (h) is Zero State Detectability (ZSD). Moreover, setting u(x, v) = α(x) + v with external control v ∈ R, then the closed loop dynamics is average passive again with respect to h

S(x + α (v)) -S(x) ≤ vh av α (x, v). (15)
To conclude average passivity of the closed-loop dynamics it is sufficient to note that under the feedback u(x, v) = α(x) + v, the average output associated with h along the closed-loop dynamics described in Lemma 2, recovers the average output from ū = α(x) defined in Definition 2. In fact, one has

h av α (x, v) := 1 v v 0 h(x + α (s))ds.
Then, by definition of α(x), one gets

α(x)+v 0 h(x + (s))ds = α(x) 0 h(x + (s))ds + α(x)+v α(x) h(x + (s))ds = α(x) 0 h(x + (s))ds + v 0 h(x + (α(x) + s))ds = -κ(h av (x, α(x))) 2 + v 0 h(x + α (s))ds
≤ vh av α (x, v) so concluding that average passivity from α(x) coincides with average passivity under preliminary feedback α(x) with respect to the same output map h.

The negative average output feedback solving [START_REF] Jakubczyk | Automatique théorique. orbites de pseudo-groupes de difféomorphismes et commandabilité des systèmes non linéaires en temps discret[END_REF] is the first step towards stabilizing strategy through additional damping or average PBC feedback. Remark 4. Computational aspects are beyond the scope of the paper. However, we underline that the control solution being expressed as the solution to the algebraic equality of the form [START_REF] Jakubczyk | Automatique théorique. orbites de pseudo-groupes de difféomorphismes et commandabilité des systèmes non linéaires en temps discret[END_REF], its computation may be a difficult. In practice, it can be performed according to suitable approximation methods as discussed in [START_REF] Mattioni | Forwarding stabilization in discrete time[END_REF][START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF].

The second fundamental property to verify regards the interconnection of two average passive systems through their respective input and output variables when setting u = ϕ(h av (x, u)). Given two average passive systems Σ i d (h i ) with storage function S i for i = 1, 2, a power-preserving input-output interconnection making reference to the average outputs can be naturally defined as in [START_REF] Moreschini | Interconnection through u-average passivity in discrete time[END_REF][START_REF] Moreschini | Dirac structures for a class of port-Hamiltonian systems in discrete time[END_REF]. Definition 4 (Power preserving interconnection). The input-output interconnection between Σ d (h 1 ) and Σ d (h 2 ) is said power preserving if it satisfies the integral equality

u 1 0 h 1 (x 1 + (s))ds + u 2 0 h 2 (x 2 + (s))dw = 0 ( 16
)
equivalently rewritten as

u 1 h 1 av (x 1 , u 1 ) + u 2 h 2 av (x 2 , u 2 ) = 0. ( 17
)
We easily note that the simplest way to solve ( 17) is to set

u 1 u 2 = ϕ h 1 av (x 1 , u 1 ) h 2 av (x 2 , u 2 ) = 0 -1 1 0 h 1 av (x 1 , u 1 ) h 2 av (x 2 , u 2 ) ( 18
)
so recovering the classical power preserving interconnection expressed with respect to the average outputs. The solution to the implicit equality (17) defines a preliminary power preserving state-feedback that we denote α(x) = (α 1 (x), α 2 (x)) ⊤ with x = (x 1 , x 2 ) ⊤ . The following Theorem can be stated.

Theorem 4 (Average passivity under power preserving interconnection). Let, for i = 1, 2, the systems Σ d (h i ) be average passive with respective storage functions S i . Let α(x) be the power-preserving interconnection satisfying [START_REF] Maschke | Port-controlled Hamiltonian systems: modelling origins and system theoretic properties[END_REF] and set u = α(x) + v with external control v = (v 1 , v 2 ) ⊤ . Then, the interconnected system

x 1+ α 1 (v 1 ) =F 1 α 1 (x, v 1 ) (19a) x 2+ α 2 (v 2 ) =F 2 α 2 (x, v 2 ) ( 19b
)
with output h = (h 1 , h 2 ) ⊤ is average passive with storage function S(x) := S 1 (x 1 ) + S 2 (x 2 ). Namely, the dissipation inequality holds; i.e.

S 1 (x 1+ α 1 (v 1 )) -S 1 (x 1 ) + S 2 (x 2+ α 2 (v 2 )) -S 2 (x 2 ) ≤ v ⊤ h av α (x, v) (20)
with the average output of the closed-loop dynamics defined as

h av α (x, v) = 1 v 1 v 1 0 h 1 (x 1+ α 1 (s))ds 1 v 2 v 2 0 h 2 (x 2+ α 2 (s))ds . ( 21
)
The proof, detailed in [START_REF] Moreschini | Interconnection through u-average passivity in discrete time[END_REF], works out in two steps. First, one concludes that under power-preserving input-output interconnection, the interconnected system is average passive from α(x) according to the Definition 2. Then, because of the feedback structure, average passivity from α recovers average passivity of the dynamics under preliminary feedback α(x), so concluding the claim.

Passivating output map

Theorem 3 shows how stabilization under feedback can be achieved by exploiting average passivity. However, stabilization to some target equilibria that are local extrema of suitably shaped energy functions can be requested. For, the following proposition is instrumental, it specifies a dummy output function that preserves average passivity and is zero at the local minima of the storage function. It is reminiscent of the continuous-time context (see [START_REF] Monaco | Nonlinear average passivity and stabilizing controllers in discrete time[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF] for details).

Proposition 2 (Average passivating output). Let Σ d (h) be average passive with storage function S, then it is also average passive with respect to the dummy output function

Y (x, u) = L G(•,u) S(x). (22)
Y (•, u), computed as the Lie derivative of S along G(•, u), is referred to as an average passivating output because it satisfies the Energy Balance Equality (EBE)

S(x + (u)) -S(x) stored energy = S(x + (0)) -S(x) dissipated energy + S(x + (u)) -S(x + (0)) supplied energy
.

Average passivity with respect to the output [START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF] holds since by definition

S(x + (u)) -S(x + (0)) = u L G(•,u) S av (x, u) = uY av (x, u)
with S(x + (0)) -S(x) ≤ 0 from the average passivity assumption of

Σ d (h). Moreover, Y (x, u) is zero at local extrema of S since Y (x, u) = ∂S ∂x G(x, u).
Remark 5. Proposition 2 generalizes to assuming Σ d (h) average passive from a given ū. In that case, average passivity from ū with respect to Y (•, u) follows

S(x + (ū + u)) -S(x) ≤ S(x + (ū + u)) -S(x + (ū)) = uY av ū (x, u)
because by assumption S(x + (ū)) -S(x) ≤ 0 and by definition

L Gū(•,u) S av (x, u) = 1 u u 0 L G(•,ū+s) S(x + (ū + s))ds = Y av ū (x, u) with G ū(•, s) = G(•, ū + s).
Remark 6. Specifying the result in Theorem 3 on such output function [START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF], the stabilizing feedback u = α(x) results to be the solution to the algebraic equality

u = -κ L G(•,u) S av (x); κ > 0
that can be solved in first approximation around u = 0 so getting

α ap (x) = -κλ(x) ∂S(x) ∂x G 1 (x) x=x+F0(x)
with a suitable gain λ(x) > 0 as discussed in [START_REF] Mattioni | Forwarding stabilization in discrete time[END_REF]. Setting u = α(x)+v, the closed loop dynamics with output L Gα(•,v) S(x) = Y α (x, v) remains average passive.

In the present paper, oriented to characterize Hamiltonian dynamics in discrete time, it is instrumental to relate the passivating output defined in [START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF] to a certain discrete gradient function. The definition below is recalled from [START_REF] Mclachlan | Geometric integration using discrete gradients[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF].

Definition 5 (Discrete gradient function). Given a smooth real-valued function S : R n → R, its discrete gradient is a function of two variables ∇S| z x :

R n × R n → R n satisfying for all x, z ∈ R n (z -x) ⊤ ∇S| z x = S(z) -S(x) with lim z→x ∇S| z x = ∇S(x). ( 23
)
Definition 5 properly states that the discrete gradient function satisfying [START_REF] Mattioni | Quaternion-based attitude stabilization via discrete-time IDA-PBC[END_REF] describes the rate of change of this function between two states. It is not uniquely defined and different methods to solve the equality can be worked out [START_REF] Gonzalez | Time integration and discrete Hamiltonian systems[END_REF][START_REF] Mclachlan | Geometric integration using discrete gradients[END_REF][START_REF] Haier | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]. Through component-wise integration, one gets the computable expression below

∇S| z x = ∇1 S| z1 x1 . . . ∇n S| zn xn ⊤ with ∇i S| zi xi = 1 z i -x i zi xi ∂S(x 1 , ..., x i-1 , s, z i+1 , ..., z n ) ∂s ds.
Remark 7. When S(x) = 1 2 x ⊤ P x with P ∈ R n×n , the discrete gradient is uniquely expressed as

∇S| z x = 1 2 P (x + z). ( 24
)
By definition of the discrete-gradient function, the EBE in Proposition 2 rewrites

S(x + (u)) -S(x) = S(x + (0)) -S(x) + S(x + (u)) -S(x + (0)) = S(x + (0)) -S(x) =F ⊤ 0 (x) ∇S| x + (0) x + u 0 L G(•,s) S(x + (s))ds =ug ⊤ (x,u) ∇S| x + (u) x + (0)
.

It is instrumental to describe the passivating output defined in [START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF] in terms of the discrete gradient function. One gets.

Lemma 3 (Average passivating output in discrete gradient form). Given the dynamics (1) and a real-valued smooth function map S : R n → R, the following equalities hold

S(x + (0)) -S(x) = S(x + F 0 (x)) -S(x) = F ⊤ 0 (x) ∇S| x + (0) x S(x + (u)) -S(x + (0)) = S(F (x, u)) -S(x + F 0 (x)) = ug ⊤ (x, u) ∇S| x + (u)
x + (0) .

with the relation

ug ⊤ (x, u) ∇S| x + (u) x + (0) = u 0 L G(•,s) S(x + (s))ds = u L G(•,u) S av (x, u). (25)
Example 1. The discrete integrator

x + (u) = x + u; y = h(x) = x
is the simplest storage element. Setting S(x) = 1 2 x 2 , as storage function, the system is average passive as

S(x + (u)) -S(x) = xu + 1 2 u 2 = u 0 x + (v)dv = u 0 (x + s)ds = uh av (x, u)
with h av (x, u) = x + 1 2 u. Accordingly, one conclude passivity of the input-stateoutput system

x + (u) = x + u; h av (x, u) = x + 1 2 u.
The associated negative average output feedback satisfies the algebraic equality

u = -κ(x + 1 2 u) so computing u = -κ 1+ κ 2
x; κ > 0 that recovers negative output feedback with suitably shaped gain.

Concluding comments

In this section, the notion of average passivity is introduced and is shown to be qualifying for the design of discrete-time average passivity based stabilizing strategies through damping (Theorem 3) or interconnection (Theorem 4). A passivating output and its related negative average output feedback are described in Proposition 2, so enlarging the control objectives to stabilization to target equilibria. This method of passivation can be further exploited to describe cascade stabilizing procedures for triangular state dynamics through backstepping or feedforward strategies in a discrete-time context, as in [START_REF] Mattioni | Forwarding stabilization in discrete time[END_REF]. Further on, the second generation of average passivity based control including an energy-shaping component to shape the energy of the system and fulfil required control specifications can be developed. Such extensions include Interconnection and Damping Assignment (IDA) techniques that modify the internal port-Hamiltonian structure to assign a new equilibrium, or Control by Interconnection (CbI) techniques that manage energy exchanges through an interconnection pattern. Preliminary works in this direction are in [START_REF] Moreschini | Discrete port-controlled Hamiltonian dynamics and average passivation[END_REF][START_REF] Moreschini | Interconnection through u-average passivity in discrete time[END_REF][START_REF] Moreschini | Stabilization of discrete port-Hamiltonian dynamics via interconnection and damping assignment[END_REF][START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF]. Finally, Lemma 3, that rewrites the average passivating output in its discrete gradient form, directly inspires the novel port-control Hamiltonian structure we propose in the next section.

Port-Hamiltonian structures in discrete time

Port-Hamiltonian structures have a pervasive impact in numerous applied domains enlarging the more traditional mechanical one. These structures are more essentially described in the continuous-time domain while in discrete time, a consensus on a specific structure is not reached in spite of a rich literature. In this section, a novel description of port-Hamiltonian structures is proposed exploiting the DDR form and the average passivating output map introduced in Proposition 2.

Control-free port-Hamiltonian dynamics

From [START_REF] Moreschini | Discrete port-controlled Hamiltonian dynamics and average passivation[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF][START_REF] Moreschini | Dirac structures for a class of port-Hamiltonian systems in discrete time[END_REF], we first recall the definition of a control-free port-Hamiltonian structure. A unified definition that mimics the continuous-time structure exists in the discrete-time framework, just replacing the gradient function with the discrete gradient function. Let H : R n → R ≥0 , be a smooth real-valued function that denotes the Hamiltonian function. Definition 6. A control-free discrete-time port-Hamiltonian dynamics over R n can be described by the first-order difference equation

x + -x = (J(x) -R(x)) ∇H| x + x ( 26
)
where J(x) = -J ⊤ (x), R(x) = R ⊤ (x) ⪰ 0, are matrices of functions representing the interconnection and resistive parts respectively.

By construction, one immediately verifies that:

• any local extremum of H(x) ( ∇H| xe xe = ∇H(x e ) = 0), is an equilibrium; • the rate of change of the Hamiltonian along the dynamics satisfies the equality

H(x + ) -H(x) = -∇H ⊤ | x + x J(x) ∇H| x + x =0 -∇H ⊤ | x + x R(x) ∇H| x + x ≤0
by skew symmetry of J(x) and semi-positiveness of R(x). Taking the sum of these increments, energy dissipation from time 0 to time k is described by the equality

H(x k ) -H(x 0 ) = - k-1 i=0 ∇H ⊤ | x + i xi R(x i ) ∇H| x + i xi dissipated energy≤0
.

• When J(x) = 0, the dynamics is dissipative. The simplest example is the gradient dynamics defined with R(x) = I and J(x) = 0 that is

x + -x = -∇H| x + x satisfying H(x + ) -H(x) = -∇H ⊤ | x + x ∇H| x + x = -|| ∇H| x + x || 2 ≤ 0.
• When R(x) = 0, the dynamics is conservative

H(x + ) = H(x)
so concluding that the Hamiltonian is a constant of motion for (26).

Canonical discrete Hamiltonian dynamics

Let us illustrate the proposed definition on a peculiar class. Canonical Hamiltonian dynamics are defined over R 2n , when setting as skew symmetric interconnection matrix

J c = 0 I d -I d 0 . Setting x = (x 1 , x 2 ) ⊤ , x i ∈ R n for i = 1, 2 and ∇H| x + x = col( ∇1 H| x 1+ x 1 , ∇2 H| x 2+
x 2 ), the canonical discrete Hamiltonian vector field associated to a given H, denoted by XH , satisfies

x 1+ -x 1 x 2+ -x 2 = J c ∇H| x + x = ∇2 H| x 2+ x 2 -∇1 H| x 1+ x 1 = XH (27)
For completeness, we note that for a given Hamiltonian function over R 2n , the canonical Hamiltonian dynamics is solution for all v of the condition

Ω(x + -x, v) = ∇H| x + x , v
where Ω(u, v) =< u, J c v > denotes the usual symplectic form. In fact, easy computations show that this equality rewritten as

< (x + -x), J c v >= ∇H(x), v
is solved by x + -x = XH defined in [START_REF] Monaco | A lie exponential formula for the nonlinear discrete time functional expansion[END_REF].

Further on, for any given real-valued smooth function C : R 2n → R, its rate of change along the Hamiltonian dynamics XH is given by

C(x + ) -C(x) = {C, H} D
where {C, H} D indicates the discrete Poisson bracket interestingly defined as the usual Poisson bracket but with respect to the discrete gradient; i.e.

{C, H}

D := n i=1 ∇1i C| x 1i+ x 1i ∇2i H| x 2i+ x 2i -∇2i C| x 2i+ x 2i ∇1i H| x 1i+ x 1i .
Canonical discrete Hamilton's equations can be alternately written as

C(x + ) -C(x) = {C, H} D , ∀C : R 2n → R.
In our formalism, any function C satisfying {C, H} D = 0 is referred to as a discrete integral or constant of the motion with respect to the discrete dynamics generated by XH .

Port-controlled Hamiltonian structures

A novel description of port-controlled Hamiltonian structures exploiting the DDR form and the passivating average output map defined in Proposition 2 can now be proposed. This form is further validated by the Energy Balance Equation that it satisfies and its relation with feedback strategies and power-preserving interconnection. The definition below is recalled from [START_REF] Moreschini | Discrete port-controlled Hamiltonian dynamics and average passivation[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF].

Definition 7 (Port-controlled Hamiltonian system (pH)). Given a smooth real-valued function H : R n → R ≥0 , a discrete-time port-Hamiltonian system Σ H d over R n can be described according to the input-state-output form below

x + (0) = x + (J(x) -R(x)) ∇H| x + (0) x (28a) dx + (u) du = G(x + (u), u) (28b) Y (x, u) = L G(•,u) H(x) (28c)
J(x) = -J ⊤ (x), R(x) = R ⊤ (x) ⪰ 0, are matrices of functions representing the interconnection and the resistive parts.

Accordingly, the following result holds.

Theorem 5. Given a discrete-time port-Hamiltonian system of the form (28), then the following holds:

• any local extremum of H(x) is an equilibrium;

• the rate of change of the Hamiltonian along the dynamics satisfies

H(x + (u)) -H(x) = -( ∇H| x + (0) x ) ⊤ R(x) ∇H| x + (0) x ≤0 + u 0 L G(•,s) H(x + (s)) =uY av (x,u) ds ( 29
)
so concluding average passivity with respect to the output map (28c).

Some comments are in order. We note that the so defined output (28c), that corresponds to the passivating output introduced in Proposition 2 when substituting the storage function S with the Hamiltonian function H, qualifies as conjugate output: the product uY av (x, u) describes the energy brought to the system through the external input and output variables between two successive time steps. Taking the sum of each increment (29) from time 0 to k, one gets the Energy Balance Equality in a form that perfectly splits in the total stored energy, the internally dissipated energy from the one supplied by the input/output variables:

H(x k ) -H(x 0 ) stored energy = - k-1 i=0 ∇H ⊤ | x + i (0) xi R(x i ) ∇H| x + i (0) xi dissipated energy + k-1 i=0 u i Y av (x i , u i ) supplied energy .
As an alternative to the DDR form of Σ H d in [START_REF] Monaco | Nonlinear systems in discrete time[END_REF], integration with respect to u transforms the port-Hamiltonian system into its map form. Adopting the discrete gradient form representation of the average output described in Lemma 3, one gets equivalently to (28) the port-Hamiltonian structure in map form. The following proposition specifies the equivalence between these two forms.

Proposition 3 (Port-Hamiltonian systems in map form). The pH structure (28) can be equivalently rewritten in map form as

x + (u) = x + (J(x) -R(x)) ∇H| x + (0) x + ug(x, u) (30a) y(x, u) = g ⊤ (x, u) ∇H| x + (u) x + (0) (30b)
where by definition and from Lemma 3

ug(x, u) := u 0 G(x + (s), s)ds; y(x, u) = Y av (x, u) = H(x + (u)) -H(x + (0)) u .
Remark 8. The average representation of the conjugate output is instrumental to describe its series expansion in power of u that gives for the first terms

y(•, u) = Y av (x, u) = L G1 H| x + + u 2 L 2 G1 + L G2 H| x + + O(u 2 ) (31)
where O(u 2 ) contains all the remaining terms of higher order in the control variable u. Further details regarding the complete series expansion and its iterative computation are in [START_REF] Monaco | Nonlinear average passivity and stabilizing controllers in discrete time[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF].

It is interesting to highlight a matrix representation of the port-Hamiltonian structure we propose as a preamble to describe the associated Dirac structure. Easy computations show that the equations (28) satisfy Lemma 4. Easy computations show that the equations (28) satisfy

  x + (0) -x dx + (u) -Y (x + (u), u)   =   J(x) -R(x) 0 0 0 0 G(x + (u), u) 0 -G ⊤ (x + (u), u) 0     ∇H| x + (0) x ∇H(x + (u)) du  
that underlines the hybrid representation that couples a one step ahead difference regarding the free evolution with a differential form with respect to u-variation.

Average PBC strategies for port-Hamiltonian systems

Port-controlled Hamiltonian systems represent a common class of average passive systems over which PBC strategies can be applied. Specifying the result in Theorem 3, one gets.

Theorem 6. Let the port-Hamiltonian system Σ H d described in (28) (equivalently in ( 30)) be ZSD with H having a minimum in x ⋆ . Then, the feedback α(x), solution of the implicit damping equality

α(x) = -κY av (x, α(x)) = - κ α(x) α(x) 0 L G(•,s) H(x + (s))ds (32)
equivalently rewritten in terms of the discrete gradient function as

α(x) = -κg ⊤ α (x) ∇H| x + (α(x)) x + . ( 33
)
with g α (x) := g(x, α(x)) achieves asymptotic stabilization of the equilibrium x ⋆ .

Accordingly the closed loop dynamics rewrites

x + (α(x)) = x + (J(x) -R(x)) ∇H| x + x -κg α (x)g ⊤ α (x) ∇H|

x + (α(x))

x + . ( 34
)
Setting now u(x, v) = α(x) + v and x + α (v) = x + (α(x) + v), the closed loop port-Hamiltonian structure can be represented in matrix form over R 3n+1 as Theorem 6 is the first step towards a variety of stabilizing techniques. Specifying the result in Theorem 4 to port-Hamiltonian structures, one gets preservation of the port-Hamiltonian structure under power-preserving interconnection; a qualifying property to discuss energy management based control schemes. Examples in this direction are developed in [START_REF] Moreschini | Stabilization of discrete port-Hamiltonian dynamics via interconnection and damping assignment[END_REF][START_REF] Moreschini | Dirac structures for a class of port-Hamiltonian systems in discrete time[END_REF][START_REF] Moreschini | A gradient descent algorithm built on approximate discrete gradients[END_REF][START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF][START_REF] Mattioni | Quaternion-based attitude stabilization via discrete-time IDA-PBC[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF].

    x + -x x + (α(x)) -x + dx + α (v) -Y α (x + α (v), v)
Remark 9. For completeness, we report the port-controlled Hamiltonian structure usually proposed in the literature [START_REF] Laila | Discrete-time IDA-PBC design for separable Hamiltonian systems[END_REF][START_REF] Aoues | Hamiltonian systems discrete-time approximation: Losslessness, passivity and composability[END_REF][START_REF] Yalçin | Discrete-time modeling of Hamiltonian systems[END_REF]. One sets

x + (u) = x + (J(x) -R(x)) ∇H| x + (u)

x + ug lit (x, u) (35a)

h lit (x, u) = g ⊤ lit (x, u) ∇H| x + (u) x (35b)
where, with respect to the form we propose, the discrete gradient of H from x to x + (0) is substituted with the discrete gradient of H from x to x + (u) plus an additive controlled part ug lit (x, u). It results that the rate of change of H between two successive time instants rewrites as

H(x + (u)) -H(x) = -( ∇H| x + (u) x ) ⊤ R(x) ∇H| x + (u)
x ≤0 +ug ⊤ lit (x, u) ∇H| x + (u)

x [START_REF] Monaco | From passivity under sampling to a new discrete-time passivity concept[END_REF] so naturally concluding passivity with respect to the output map h lit (x, u). However, the fact that the resistive part depends on the input variable in an unpredictable way through the term ∇H|

x + (u) x may represent an obstacle to managing damping or energy exchanges under feedback. As a consequence in [START_REF] Monaco | From passivity under sampling to a new discrete-time passivity concept[END_REF], the inner product ug ⊤ lit (x, u) ∇H|

x + (u) x

, does not contain the total power supplied to the system because the resistive part has a control dependent element too.

Example 2. Specifying the state equations (28) (equivalently [START_REF] Monaco | A unified representation for nonlinear discrete-time and sampled dynamics[END_REF]) to a linear dynamics with quadratic Hamiltonian function H(x) = 1 2 x ⊤ P x and symmetric positive matrix P , a linear Port Hamiltonian structure can be be described as

x + (0) = x + (J -R) P 2 (x + x + (0); dx + (u) du = B; Y (x) = B ⊤ P x [START_REF] Monaco | Nonlinear hamiltonian systems under sampling[END_REF] or equivalently in map form as

x + (u) = x + (J -R) P 2 (x + x + (0)) + Bu; Y av (x, u) = B ⊤ P 2 (x + (0) + x + (u)) [START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF] with matrices of appropriate dimensions and constant elements. Because the discrete gradient function can be explicitly expressed as a function of x and x + (0), the state equations (37) (equivalently [START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF]) can be rewritten in their explicit form so getting a linear dynamics with an output map admitting a feedthrough term as it is required to encompass a passivity property. One gets

x + (u) = A H x + Bu; y(x, u) = B ⊤ P A H x + B ⊤ P B 2 u (39)

v v 0 L

 0 α (x)g ⊤ α (x) ⪰ 0 expressing the modified closed-loop dissipation matrix. The conjugate output Y av α (x, v), computed as the average from α(x) of the output Y (x, u) = L G(•,u) H(x) can be rewritten in discrete gradient form asY av α (x, v) = 1 Gα(•,s) H(x + α (s))ds = g⊤ (x, v) ∇H| + (α(x) + s), α(x) + s)ds = v 0 G α (x + α (s), s)ds= α(x)g(x, α(x) + v) -α(x)g(x, α(x)) + vg(x, α(x) + v).

Acknowledgments

Dorothée Normand-Cyrot expresses her gratitude to the Committee Members for the invitation to give a plenary lecture at the ICDEA 2022.

with

Easy computations show that the output y(x, u) is exactly the average output associated with Y (x) = B ⊤ P x since by definition

It is also possible to raise the question: when and how a passive system satisfies a pH structure? In the linear case, the question can be answered.

Proposition 4. Consider the average passive linear system

with positive definite storage 1 2 x ⊤ P x, then it can be rewritten in the port-Hamiltonian form (37) (equivalently ( 38)) according to the decomposition in skew-symmetric and symmetric part as follows

Concluding comments

In the proposed differential algebraic framework, discrete-time Port-controlled Hamiltonian structures that validate the usually required energy balance properties are described in the proposed differential algebraic framework. Accordingly, the basic stabilizing techniques behind energy-based control strategies are revisited to confirm the open perspectives regarding control design through energy management along the lines developed in [START_REF] Moreschini | Stabilization of discrete port-Hamiltonian dynamics via interconnection and damping assignment[END_REF][START_REF] Mattioni | Discretetime energy-balance passivity-based control[END_REF][START_REF] Mattioni | Quaternion-based attitude stabilization via discrete-time IDA-PBC[END_REF][START_REF] Moreschini | Modeling and control of discrete-time and sampleddata port-Hamiltonian systems[END_REF]. All the material discussed in this paper regards a purely discrete time setting but it can be specified to the sampled-data context. How port-Hamiltonian structures are transformed under sampling when assuming both the measurements and control variables available at discrete time instants is thus a natural and challenging question addressed in [START_REF] Monaco | A unified representation for nonlinear discrete-time and sampled dynamics[END_REF][START_REF] Monaco | From chronological calculus to exponential representations of continuous and discrete-time dynamics: A lie-algebraic approach[END_REF][START_REF] Monaco | Nonlinear hamiltonian systems under sampling[END_REF]. In that digital framework, it comes out that the sampleddata dynamics are necessarily parameterized by the sampling period as well as the control solutions that are described around the continuous-time ones by infinite series expansions. As a result, the solutions can be computed through an iterative procedure and approximated at any order so rendering the approach constructive in a digital environment. Finally, we stress that in our opinion the proposed approach sets an unifying framework to investigate controlled finitedimensional dynamical systems in discrete-time as well as under sampling.


Morten Hovd

Giorgio Valmorbida

Solving LP-MPC problems using ramp functions*

A new implementation of the Simplex method for solving linear programming problems is developed, and its application for solving MPC problems with linear objective functions is described. A detailed comparison with existing implementations of the Simplex method is beyond the scope of this paper. However, initial numerical results indicate that gains in calculation times may be achieved for problems of small to moderate size.

I. INTRODUCTION

Model Predictive Control (MPC) has been a great success in industry [START_REF] Qin | A survey of industrial model predictive control technology[END_REF], and since its initial development in the 1970's (see, e.g. [START_REF] Richalet | Model predictive heuristic control: Application to industrial processes[END_REF]) it has found application in a wide range of industrial processes. The main feature of MPC distinguishing it from classical control design methods such as Linear Quadratic Regulator is the ability to take into account constraints in both inputs, state, and outputs. However, in the conventional MPC formulation, an optimization problem has to be solved online, thus incurring a significant computational cost. Many approaches have been studied in order to reduce the computational requirements of MPC, including input blocking [START_REF] Cagienard | Move blocking strategies in receding horizon control[END_REF] and utilizing structure in the optimization problem (e.g. [START_REF] Rao | Application of interior-point methods to model predictive control[END_REF]). This paper will make no attempt at covering all such approaches. Instead, the rampbased approach to solving quadratic programming (QP) problems, proposed recently in [START_REF] Valmorbida | Quadratic programming with ramp functions and fast online QP-MPC solutions[END_REF], is extended to linear programming (LP) problems.

The next section describes the LP-MPC formulation used. Thereafter, the main result of the paper, the rampbased solution to LP problems in standard form, is presented, followed by a description of how the Simplex method can be used with the ramp-based solution. Numerical experiments are performed to compare the calculation performance of the new approach to that of of the well known commercial solver CPLEX [START_REF] Ibm | IBM ILOG CPLEX optimizer[END_REF] from IBM. Initial results indicate computational advantages for small to moderate-size problems.

II. LP-MPC FORMULATION

A linear programming problem in standard for can be expressed as *This work was supported by the French-Norwegian exchange programme AURORA under contract 294676. 1 Morten Hovd is with the Department of Engineering Cybernetics, Norwegian University of Science & Technology, 7491 Trondheim, Norway. morten.hovd@itk.ntnu.no Saclay, 3 rue Joliot Curie, Gif-sur-Yvette, 91192 France. giorgio.valmorbida@centralesupelec.fr

minimize z c ⊤ z (1a) subject to Gz = b (1b) -z ≤ 0 (1c)
where

G ∈ R ne×nz , c ∈ R nz×1 , z ∈ R nz×1 , and b ∈ R ne×1 .
Naturally it is assumed that n e < n z , as otherwise there would be no degrees of freedom for optimization.

The system dynamics are given by the linear discretetime state space model

x(k + 1) = Ax(k) + Bu(k) (2a) y(k) = Cx(k) (2b)
The constraints on the inputs u and outputs y are given by polytopic sets 1 of the form

y min ≤ y ≤ y max (3a) u min ≤ u ≤ u max . (3b)
We wish to compute the input sequence that minimizes the cost

J(y, u, x(k + N )) =∥P x(k + N )∥ p + (4) N -1 i=1 ∥Qy(k + i)∥ p + ∥Ru k+i-1 ∥ p
where y and u denote the concatenated vectors y =

[y(k + 1) ⊤ y(k + 2) ⊤ • • • y(k + N ) ⊤ ] ⊤ and u = [u(k) ⊤ u(k + 1) ⊤ • • • u(k + N -1) ⊤ ] ⊤ , respectively, and Q ∈ R nc×ny , R ∈ R nu×nu .
To have a sensible LP-MPC formulation, it is naturally assumed that (A, B) is stabilizable and (Q

1 2 C, A
) is detectable and R full rank. The purpose of this paper is not the LP-MPC solution per se, but rather to present a numerical solution strategy to compute the input values u. Although terminal sets and costs are important for recursive feasibility and stability, and systematic approaches to calculating these do exist also for LP-MPC [START_REF] Danielson | Terminal-cost design for model predictive control with linear stage-costs: A set-theoretic method[END_REF], they will for simplicity be ignored in the following, as they are not important for presenting the contributions of the paper. Furthermore, p = 1 will be assumed throughout the paper, although LP-MPC using p = ∞ can also be formulated in a straight forward way.

In the following, it will be described how to translate the objective function (4), dynamics (2a,2b) and constraints (3a,3b) into an LP problem formulation on the form (1a)-(1c). This development is rather standard, but the presented level of detail is necessary to present the LP-MPC solution procedure that is the main contribution of this paper. Repeated use of the system dynamics dynamics (2a,2b) allow us to express the future measurements as

y = C A x(k) + C B u (5)
where

A 0 =      A 0 . . . 0      B = diag{B} C = diag{C} I A =         I 0 • • • 0 0 -A I 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . -A I 0 0 • • • 0 -A I         and C A = CI -1 A A 0 and C B = CI -1 A B.
The output constraints trivially translate to

y ≤ 1 N ⊗ y max (6a) -y ≤ -1 N ⊗ y min , (6b)
where 1 N denotes a length N column vector of ones, and ⊗ denotes the Kronecker product. Similarly, the input constraints are expressed as

u ≤ 1 N ⊗ u max (7a) -u ≤ -1 N ⊗ u min . (7b)
In order to translate the objective function in (4) with p = 1 into the form used in (4), we split y and u into positive and negative parts

y = y p -y m u = u p -u m .
Introducing auxiliary variables s p , s m , s up , s um , the equality constraint (5) together with inequality constraints (6a-7b) are converted to the equality constraint (1b), with

G =       I -I -C B C B 0 0 0 0 I 0 0 0 I 0 0 0 0 -I 0 0 0 -I 0 0 0 0 I 0 0 0 I 0 0 0 0 -I 0 0 0 I       (8) b =       C A x(k) 1 N ⊗ y max 1 N ⊗ y min 1 N ⊗ u max 1 N ⊗ u min       (9) where z = [y p ⊤ y m ⊤ u p ⊤ u m ⊤ s p ⊤ s m ⊤ s up ⊤ s um ⊤ ] ⊤
in (1b) and (1c).

The objective function now follows from (1a) with

c ⊤ = c ⊤ y c ⊤ y c ⊤ u c ⊤ u 0 0 0 0 where c y = 1 ⊤ N ⊗ σ(Q) c u = 1 ⊤ N ⊗ σ(R)
and σ(Q) and σ(R) are row vectors of column sums of Q and R, respectively, such that, for p = 1, we obtain

J(y, u, x(k + N )) = c ⊤ z.

III. Main results

This section first provides some background on ramp functions, as presented in [START_REF] Valmorbida | Quadratic programming with ramp functions and fast online QP-MPC solutions[END_REF]. Next, it is shown how a Linear Complementarity Problem (LCP) can be associated to an LP to obtain a PWA representation implicitly defined in terms of ramp functions. Section IV details how this implicit PWA representation can be exploited to quickly calculate solutions to the LP-MPC problem.

A. Ramp functions

Definition 1: The ramp function r(y) is given by

r(y) = 0 if y < 0 y if y ≥ 0 (10) Lemma 1: The ramp function is the only function satisfying, ∀y ∈ R (r(y) -y)r(y) = 0 (11a) r(y) ≥ 0 (11b) (r(y) -y) ≥ 0. ( 11c
) Proof: First note that the ramp function can be expressed as the unique solution to the convex optimization problem parameterized in y (thus depending on y) with strictly convex objective function as follows minimize

r 1 2 (r -y) 2 subject to r ≥ 0. ( 12
)
With the Lagrangian associated to the optimization problem, L(r, λ) = 1 2 (r -y) 2 -λr, we obtain the Karush-Kuhn-Tucker (KKT) conditions (r -y) -λ = 0; λr = 0; r ≥ 0; λ ≥ 0 which are necessary for optimality and also sufficient since the problem is strictly convex. To obtain a description in the variables (y, r) one can use λ = (r -y) above to obtain r ≥ 0, (r -y) ≥ 0, r(r -y) = 0. ■ For y ∈ R m let us define the function r : R m → R m , the vector-valued ramp function, as

r(y) = r(y 1 ) r(y 2 ) • • • r(y m )
⊤ Remark 1: Note that, due to the piecewise definition of the ramp function in (10), we have that vector r can also be expressed as the product r(y) = I A y where I A ∈ D m , with its diagonal elements verifying I A(i,i) ∈ {0, 1}. Clearly, the set A ⊆ {1, 2, . . . , N } depends on y. ⌟

B. Expressing the solution to an LP problem using ramp functions

Theorem 1: The solution of the LP-MPC is a PWA function given by

u(x) = M 1 x + M 2 r(y) + m 0 y = M 3 x + M 4 r(y) + m 5 (14) Proof: From the Lagrangian function L(z, λ) = c ⊤ z + µ ⊤ (Gz -b) -λ ⊤ z, we get the KKT conditions c + G ⊤ µ -λ = 0 (15a) Gz = b (15b) -z ≤ 0 (15c) λ ≥ 0 (15d) λ i z i = 0 ∀i (15e)
With the matrices and vectors partitioned (and reordered as necessary) as

G = G a G b , λ = λ a λ b , z = z a z b
with square and invertible G b . Set

λ a λ b = r(y a ) -y a r(y b ) , (16)
we can reformulate the KKT conditions as

c p + G ⊤ a -I a G ⊤ b 0 µ λ a - 0 I b λ b = 0 (17a) G a G b z a z b = b (17b) z a z b ≥ 0 (17c) λ a λ b ≥ 0 (17d) λ a λ b × z a z b = 0 (17e)
where the subscript p on c p indicates that c has to be reordered after partitioning z. In contrast, it is only columns of G that are reordered, and hence there is no need for reordering vector b.

With the ramp functions assigned to λ a , λ b as in ( 16), we associate the complementarity conditions (15e) with the complementarity (11a) of the ramp function. Therefore, we obtain for z a and for z b λ ai z ai = (r(y a ) -y a ) i r(y a ) i = 0 ⇒ r(y a ) = z a (18a)

λ bi z bi = r(y b ) i (r(y b ) -y b ) i = 0 ⇒ r(y b ) -y b = z b . ( 18b
) Define Ξ = G ⊤ a -I a G ⊤ b 0 -1 which is partitioned as Ξ 1 Ξ 2 = 0 G -T b -I a G ⊤ a G -T b .
We obtain from (17a)

µ r(y a ) -y a = Ξ -c p + 0 I b r(y b )
from which we obtain

µ = -Ξ 1 c p + Ξ 1 0 I b r(y b ) (19a)
y a = Ξ 2 c p -Ξ 2 0 I b r(y b ) + r(y a ) (19b)
From (17b) we have

G a r(y a ) + G b (r(y b ) -y b ) = b
and since G b is invertible this can be expressed as

y b = r(y b ) + G -1 b G a r(y a ) -G -1 b b (20)
In summary, we obtain

z = z a z b = r(y a ) r(y b ) -y b (21a) λ = λ a λ b = r(y a ) -y a r(y b ) (21b) µ = -Ξ 1 c p + Ξ 1 0 I b r(y b ) (21c)
y a y b = -I a G ⊤ a G -⊤ b 0 0 c p - 0 G -1 b b + I -G ⊤ a G -⊤ b G -1 b G a I r(y a ) r(y b ) = M c c p + M b b + M r r(y a ) r(y b ) (21d)
Noting that for the MPC formulation above we have

b =       C A x(k) y max y min u max u min       =       C A 0 0 0 0       x(k) +       0 y max y min u max u min       = b x x + b 0 therefore, (21d) becomes y a y b = M 3 x + M 4 r(y a ) r(y b ) + m 5 (22a) M 3 = - 0 G -1 b (:, 1 : n • n y ) C A (22b) M 4 = I -G ⊤ a G -⊤ b G -1 b G a I (22c) m 5 = -I a G ⊤ a G -⊤ b 0 0 c p - 0 G -1 b       0 y max y min u max u min       (22d)
where G -1 b (:, 1 : n • ny) indicates that we only retain the n • n y leading columns of G -1 b . From (18a) and (18b) we obtain

z = I 0 -G -1 b G a 0 r(y a ) r(y b ) + G -1 b (b x x + b 0 )
and noting that u = 0 0 I -I 0 0 0 0 z = P z z we get

M 1 = -P z M 3 = P z 0 G -1 b (:, 1 : n • n y ) C A (23a)
M 2 = P z I 0 -G -1 b G a 0 (23b) m 0 = P z G -1 b b 0 (23c)

■

IV. Implementation In this section, we will first present a brief recap of the Simplex method for solving LPs, and next detail how the Simplex method can be implemented via the ramp-based description developed above.

A. The Simplex method

In the Simplex method, the decision variables z are partitioned into n e basis variables z B and n z -n e variables z 0 that are set to zero. To the variables z 0 (corresponding to active inequality constraints) we can associate Lagrange multipliers λ 0 . We have from (1b) that z B = G -1 B b, where G B consists of the columns of G corresponding to z B . A solution for which z B ≥ 0 is termed a basic feasible solution. If a Lagrange variable λ 0i < 0, this indicates that the objective function can be reduced by allowing z 0i to increase from 0, i.e., by allowing z 0i to enter the set of basis variables. If there are more than one negative Lagrange multiplier, we will assume that the first to enter the set of basis variables is the one with the most negative value. The corresponding variable to leave the active set is the first variable to become zero as variable entering the active set is increased. We find that the first basis variable to become zero is the variable z Bj corresponding to the minimum value of

δ j = z Bj G -1 B(j,:) G (:,i) = z Bj α j
The division indicated above need only be performed if α j > 0. The Simplex method assumes that a basic feasible solution is available. In each iteration, new basis variables and new Lagrange multipliers are calculated. An optimal solution is found if all Lagrange multipliers are non-negative.

B. A ramp-based implementation of the Simplex method

To recap, we are looking for solutions to (35):

y a y b = M c c p + M b b + M r r(y a ) r(y b )
We will (with some abuse of concept) call the set of variables y for which y = r(y) as the active set, denoted A, and the set of inactive variables, for which r(y) = 0, denoted A C . A solution to the LP is found if y i > 0 ∀i ∈ A and y i ≤ 0∀i ∈ A C . We also have the set of variables A for the y's corresponding to z a as well as B for the y's corresponding to z b . Thus for any consistent candidate active set, we may solve for the corresponding y by solving

(I -M r I A )y = M c c p + M b b where I A is a diagonal matrix with I A(i.i) = 1 for i ∈ A and I A(i.i) = 0 for i ∈ A C
The set of basis variables of the Simplex method correspond to

B = {i|i ∈ A ∧ i ∈ A} ∪ {i|i ∈ B ∧ i ∈ A C }
A consistent candidate active set is therefore one for which the number of basis variables equals n e , and in addition Q = (I -M r I A ) is invertible. For a consistent candidate active set we can calculate the corresponding

y from y = Q -1 (M c c p + M b b)
We then have

z i = y i , λ i = 0 if i ∈ A ∧ i ∈ A z i = 0, λ i = -y i if i ∈ A ∧ i ∈ A C z i = 0, λ i = y i if i ∈ B ∧ i ∈ A z i = -y i λ i = 0 if i ∈ B ∧ i ∈ A C
Denote by G B (as above) the matrix consisting of the columns of the matrix G for the variables in B. The inverse of this matrix is required in the Simplex method. The fact that we can extract z directly from the calculated y means that G -1 B can be found from

G -1 B = I B Q -1 B M b (24)
where I B is a diagonal matrix with +1 on the diagonal for row i if i ∈ A ∧ i ∈ A, and -1 on the diagonal of row i if i ∈ B ∧ i ∈ A C , and Q -1 B consists of the rows i of Q -1 for which i ∈ B. Furthermore, since the leading n a rows of M b are identically zero, we can reduce the computation of the

G -1 B to G -1 B = -I B Q -1 B,na+1:na+n b G -1 b
i.e., involving only the last n b columns of Q -1 B . Since y contains the Lagrange multipliers, these do not need to be calculated separately, and the identification of the variable to enter the basis is trivial. However, we still need G -1 B to calculate the α j 's required for identifying the variable to leave the active set.

Note, however, that α = G -1 B G (:,i) . Instead of calculating the whole of G B as explained above, we may instead pre-compute

G I = G -1 b G = G -1
b G a I , as storing G I (at least when storing as a full matrix) does not require more memory than storing G. At runtime, we therefore first calculate

Q I B = -I B Q -1
B,na+1:na+n b (i.e., extracting a sub-matrix from Q -1 , and flipping the sign of some rows) and then calculate α directly from

α = Q I B G I (:,i) (25)
where i is the index of the variable entering the basis.

Next, having described how to identify the variable to enter and the variable to leave the basis -and hence what variables enter and/or leave the active set, we need an efficient method for updating Q -1 .

Let k be an iteration counter, and at each iteration one variable enters and one variable leaves the basis. We then have

Q k+1 = Q k + M 4,i M 4,j s i 0 0 s j e ⊤ i e ⊤ j ( 26
)
where e i and e j are the ith and jth unit vectors, s i = 1 if variable i leaves the active set, and s i = -1 if variable i enters the active set (and similarly for s j ). Clearly, Q k+1 is a rank-2 update to Q k . Knowing Q -1 k we may thus use the Matrix inversion lemma to efficiently calculate Q -1 k+1 :

Q -1 k+1 = Q -1 k -Q -1 k U (S -1 + V Q -1 k U ) -1 V Q -1 k U = M 4,i M 4,j S = s i 0 0 s j V = e ⊤ i e ⊤ j
When calculating Q k+1 using the matrix inversion lemma, we make use of the fact that the product V Q -1 k is obtained by simply extracting rows i and j from Q -1 k . Also, if constraint i (or j) is already in the active set,

then column i of Q k is Q k(:,i) = e i -M 4(:,i) and hence Q -1 k M 4(:,i) = Q -1 k(:,i) -e i and Q -1 k(i,:) M 4(:,i) = Q -1 k(i,i) -1
Further, noting that S -1 = S, we calculate Q -1 k+1 by 1) Extract rows i and j from Q

-1 k to obtain V Q -1 k . 2) Multiply V Q -1
k and U (accounting for whether column i or j is already in the active set, as described above) 3) Invert (S -1 + V Q -1 k U ). This is a 2 × 2 matrix, so the inversion is trivial. 4) Calculate Q -1 k U (accounting for whether column i or j is already in the active set, as described above)

5) Calculate Q -1 k U (S -1 + V Q -1 k U ) -1 V Q -1 k . 6) Add Q -1 k and Q -1 k U (S -1 + V Q -1 k U ) -1 V Q -1 k .

C. Initial feasible solution

The Simplex method iteratively improves on a basic feasible solution, and hence an initial basic feasible solution must be available. The standard approach to obtaining a basic feasible solution is known as the Phase 1 of the Simplex method. In Phase 1, the problem description is augmented by slack variables, in such a way that it is trivial to define a basic feasible solution in terms of these slack variables. The objective function for Phase 1 weighs only these slack variables, and a basic feasible solution for the original problem is found if the value of objective function for Phase 1 is zero. However, the Phase 1 procedure can introduce a significant number of slack variables, and can easily find a basic feasible solution far away from the optimal point of the original problem.

Instead, we start with the basic feasible solution from the previous timestep, and follow the approach described by Maros [START_REF] Maros | Computational Techniques of the Simplex Method[END_REF] for Phase 1, where the objective function only weighs the negative basis variables. This means that the parameters of the objective function may change for each iteration of the modified Phase 1 procedure. Observe that once the updated Q -1 is available, one can extract the rows involved in calculating the basis variables. Having calculated the values of the basis variables, one can then modify the vector of weights for the objective function (based on the value of the basis variables), before updating the Lagrange multipliers. Once a basic feasible solution is found, with no negative basis variables, the objective function weights are fixed to the values in the original LP formulation.

V. NUMERICAL RESULTS

In this section, computational times for solving LP-MPC problems are compared, for the proposed rampbased LP solver as well as the commercial solver CPLEX [START_REF] Ibm | IBM ILOG CPLEX optimizer[END_REF]. It should be noted that the ramp-based solution is programmed in Matlab m-code, whereas CPLEX is called as a compiled .mex file. The MPC is simulated for 100 timesteps, starting from the same initial state for both solvers. For CPLEX the primal simplex solver is chosen2 . Both solvers are warm started from the solution at the previous timestep, for the ramp based solution this is done by specifying the initial basis, whereas for CPLEX the optimal solution from the previous timestep is specified. The two solvers produce the same system responses (as they should), and system responses are therefore not shown. The resulting computation times obtained on a Windows PC are given in Table I. To account for the fact that Windows is not a real time operating system, the average computation time over multiple runs is listed.

Example 1: Consider the double integrator example from [START_REF] Valmorbida | Quadratic programming with ramp functions and fast online QP-MPC solutions[END_REF], [START_REF] Hovd | On the design of exact penalty functions for MPC using mixed integer programming[END_REF]. The system dynamics are given by

x k+1 = 1 1 0 1 x k + 1 0.3 u k
with the input constrained between -1 and 1, and each state constrained between -5 and 5. The weight on the input is 1 and the weight on each output is 10. Prediction horizons N = {1, 5, 10} are used, to compare the growth in computation times with problem size for the two solvers. The initial state is x 0 = [5 -2] ⊤ . ⌟ Example 2: We now study the four-state system from [START_REF] Valmorbida | Quadratic programming with ramp functions and fast online QP-MPC solutions[END_REF], [START_REF] Hovd | On the design of exact penalty functions for MPC using mixed integer programming[END_REF]. The discrete-time dynamics are given by

  and each measurement constrained to lie between -1 and 1. The weight on each input is 1, and the weight on each output is 10.The prediction horizon is N = 30, and the initial state x 0 = [25.5724 25.3546 9.7892 0.2448] ⊤ . ⌟

Giorgio Valmorbida is with L2S, Laboratoire des Signaux et Systèmes, CentraleSupelec, CNRS, Université Paris-

One may have constraints on both x and y, but often it is more natural to specify the constraints in terms of y. If there are state constraints that do not directly correspond to a constraint on an output, the output vector can simply be augmented with one or more 'dummy outputs'.

It should be noted that this is not necessarily the fastest LP solver available with CPLEX.

VI. CONCLUSIONS A new implementation for the Simplex method for solving LPs, based on ramp functions, is proposed. Initial results indicate computational advantages for problems of small to moderate size. Although significantly reduced computation times may be obtained by programming the ramp-based solution method in a high-performance computer language, results on large problems indicate that significant further improvements in the solution method are required to make it competitive for large problem.