Thomas H Élie

Elementary tools on Port-Hamiltonian Systems with applications to audio/acoustics

This document presents elementary tools on Port-Hamiltonian Systems (PHS) dedicated to the energy-consistent modelling and simulation of linear or nonlinear open physical systems. It focuses on the case of finite-dimensional input-state-output representations of systems constructed as an assembly of multi-physical lumped-components, which limits the prerequisites to differential calculus and matrix algebra. The presentation and exercises are chosen with the aim of providing a toolbox to address audioacoustic problems with real time application perspectives.

Introduction

Port-Hamiltonian Sytems (PHS) introduced by A. van der Schaft and B. Maschke in [START_REF] Bernhard | Port-controlled hamiltonian systems: modelling origins and systemtheoretic properties[END_REF] (see also [START_REF] Van | Port-hamiltonian systems: an introductory survey[END_REF][START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF][START_REF] Van | Port-hamiltonian systems theory: An introductory overview[END_REF]) describe open energetically-balanced physical systems in a structured form, which naturally encodes the power balance into conservative, dissipative and external parts. This form guarantees passivity and provides a modular composable framework, passivity being inherited through composition of PHS. This document presents elementary tools on linear and nonlinear Port-Hamiltonian Systems (PHS) for the modelling and simulation, that can be used in some audio/acoustic applications.

Why PHS for musical audio/acoustic applications ? Musical instruments and audio systems are multiphysical systems, which naturally fulfill a power balance. Most of them involve nonlinearities responsible for the timbre evolution versus the excitation amplitude and for self-oscillations (e.g. bowed instruments, wind instruments, voice). Power-balanced numerical methods naturally preserve passivity, avoiding numerical instability artifacts that are often a stumbling block in nonlinear systems. Port-Hamiltonian systems provide a power-balanced and component-based approach, with a modular framework similar to that experienced by instrument makers or electronics engineers in their workshop: choose, build, refine your components and assemble them. Another important issue is control (mentioned but not addressed in this document), in order to design correctors for audio-transducers, and augmented or hybrid musical instruments (with power-balanced reprogrammed physics). This document focuses on the modelling (linear and nonlinear) and simulation, introducing tools based on differential calculus and matrix algebra. Some of these audio/acoustic applications concern (non-exhaustive list):

• Nonlinear analog audio circuits: wah pedal [START_REF] Falaize | Simulation of an analog circuit of a wah pedal: a port-Hamiltonian approach[END_REF], ondes Martenot [START_REF] Najnudel | Analog circuits and Port-Hamiltonian realizability issues: a resolution method for simulations via equivalent components[END_REF][START_REF] Najnudel | Simulation of the Ondes Martenot Ribbon-Controlled Oscillator Using Energy-Balanced Modeling of Nonlinear Time-Varying Electronic Components[END_REF][START_REF] Najnudel | Simulation of an ondes Martenot circuit[END_REF], specific components such as operational amplifiers [START_REF] Müller | A minimal passive model of the operational amplifier: application to Sallen-Key analog filters[END_REF] and nonlinear coils [START_REF] Najnudel | A power-balanced dynamic model of ferromagnetic coils[END_REF][START_REF] Najnudel | Power-Balanced Modeling of Nonlinear Coils and Transformers for Audio Circuits[END_REF][START_REF] Najnudel | From statistical physics to macroscopic port-Hamiltonian Systems: A roadmap[END_REF], circuit graph generation tools for explicit [START_REF] Falaize | Passive Guaranteed Simulation of Analog Audio Circuits: A Port-Hamiltonian Approach[END_REF] or implicit [START_REF] Müller | Fully-implicit algebro-differential parametrization of circuits[END_REF] representations, and identification of nonlinear components in circuit [START_REF] Najnudel | Identification of nonlinear circuits as port-Hamiltonian systems[END_REF], (see also the PhDs [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF][START_REF] Müller | Time-continuous power-balanced simulation of nonlinear audio circuits: realtime processing framework and aliasing rejection[END_REF])

• Electro-acoustics: loudspeaker [START_REF] Falaize | Modélisation d'un haut parleur électrodynamique: approche dans le cadre des Systemes à Hamitoniens à Ports[END_REF][START_REF] Falaize | Passive simulation of electrodynamic loudspeakers for guitar amplifiers: a port-Hamiltonian approach[END_REF][START_REF] Falaize | Passive modelling of the electrodynamic loudspeaker: from the Thiele-Small model to nonlinear port-Hamiltonian systems[END_REF], electro-mechanic piano [START_REF] Falaize | Guaranteed-passive simulation of an electromechanical piano: A port-Hamiltonian approach[END_REF][START_REF] Falaize | Passive simulation of the nonlinear port-Hamiltonian modeling of a Rhodes Piano[END_REF], (see also the PhDs [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF][START_REF] Lebrun | Modélisation multi-physique passive, identification, simulation, correction et asservissement de haut-parleur sur des comportements cibles[END_REF])

• Vibration mechanics: Duffing oscillator [START_REF] Jossic | Energy shaping of a softening Duffing oscillator using the formalism of Port-Hamiltonian Systems[END_REF], nonlinear strings and membranes [START_REF] Hélie | Corde non linéaire amortie : formulation hamiltonienne à ports, réduction d'ordre exacte et simulation à passivité garantie[END_REF][START_REF] Roze | Simulation passive d'un modèle réduit exact de plaque de Berger en grandes déformations[END_REF] with, in the case of the violin, friction modelling [START_REF] Falaize | A generic passive-guaranteed structure for elastoplatic friction models[END_REF], tom-drum [START_REF] Wijnand | Active control of the axisymmetric vibration modes of a tom-tom drum[END_REF][START_REF] Wijnand | Active control of the axisymmetric vibration modes of a tom-tom drum using a modalbased observer-regulator[END_REF], and generalised (linear and nonlinear) damping models [START_REF] Matignon | On damping models preserving the eigenfunctions of conservative systems: a port-Hamiltonian perspective[END_REF][START_REF] Matignon | A class of damping models preserving eigenspaces for linear conservative port-Hamiltonian systems[END_REF][START_REF] Hélie | Nonlinear damping models for linear conservative mechanical systems with preserved eigenspaces: a port-Hamiltonian formulation[END_REF][START_REF] Hélie | Physically-based dynamic morphing of beam sounds[END_REF], (see also the PhD [START_REF] Wijnand | Contrôle en temps fini de systèmes vibratoires hybrides couplant équations aux dérivées partielles et équations aux dérivées ordinaires : les cas du tom et du câble pesant[END_REF])

• Fluid-Structure interaction and acoustics: Brass instruments [START_REF] Lopes | Modèle d'interaction Jet/Lèvre préservant le bilan de puissance pour les instruments de type cuivre[END_REF][START_REF] Lopes | Modélisation et simulation à passivité garantie d'un cuivre[END_REF][START_REF] Lopes | Energy Balanced Model of a Jet Interacting With a Brass Player's Lip[END_REF] and vocal apparatus (larynx and vocal tract) [START_REF] Hélie | Self-oscillations of a vocal apparatus: a port-Hamiltonian formulation[END_REF][START_REF] Hélie | Port-Hamiltonian approach to self-sustained oscillations in the vocal apparatus[END_REF][START_REF] Silva | CFA2018/165 Modélisation physique, simulation à bilan de puissance garanti et examen de régimes d'un appareil vocal simplifié[END_REF][START_REF] Silva | Port-Hamiltonian Representation of Dynamical Systems. Application to Self-Sustained Oscillations in the Vocal Apparatus[END_REF][START_REF] Wetzel | Power balanced time-varying lumped parameter model of a vocal tract: modelling and simulation[END_REF][START_REF] Wetzel | Power-balanced modelling of the vocal tract: a recast of the classical lumped-parameter model[END_REF][START_REF] Silva | Energy-consistent modelling of the fluid-structure interaction in the glottis[END_REF], (see also the PhDs [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF][START_REF] Wetzel | Lumped Power-Balanced Modelling and Simulation of the Vocal Apparatus: A Fluid-Structure-Interaction Approach[END_REF])

• Control applications: loudspeaker correction with differential-flatness derived from PHS [START_REF] Falaize | Compensation of loudspeaker's nonlinearities based on flatness and port-Hamiltonian approach[END_REF], passive finite-time controlled systems for a mechanical oscillator [START_REF] Wijnand | Contrôle des vibrations d'un oscillateur passif : stabilisation en temps fini et par remodelage d'énergie[END_REF], for acoustic absorbers [START_REF] Lebrun | Electroacoustic absorbers based on passive finite-time control of loudspeakers: a numerical investigation[END_REF] and for hybrid tom-drums [START_REF] Wijnand | Active control of the axisymmetric vibration modes of a tom-tom drum[END_REF][START_REF] Wijnand | Active control of the axisymmetric vibration modes of a tom-tom drum using a modalbased observer-regulator[END_REF], real-time passive digital controller robust to computational latency [78], (see also the PhDs [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF][START_REF] Lebrun | Modélisation multi-physique passive, identification, simulation, correction et asservissement de haut-parleur sur des comportements cibles[END_REF][START_REF] Wijnand | Contrôle en temps fini de systèmes vibratoires hybrides couplant équations aux dérivées partielles et équations aux dérivées ordinaires : les cas du tom et du câble pesant[END_REF]),

• Numerical methods and automatic code generation: power-balanced numerical scheme with non-iterative solver based on quadratisation method [START_REF] Lopes | Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems[END_REF], passive integration scheme of skew-gradient systems [START_REF] Muller | Power-Balanced Modelling Of Circuits As Skew Gradient Systems[END_REF], automatic code generation with PyPHS library [77] (see also [START_REF] Falaize | PyPHS: Un module Python pour la modélisation et la simulation à passivité garantie de systèmes multi-physiques[END_REF][START_REF] Falaize | Simulation of finite-dimensional multi-physical systems described by networks of components[END_REF]), anti-aliasing for discrete-gradient method [START_REF] Muller | Trajectory Anti-Aliasing on Guaranteed-Passive Simulation of Nonlinear Physical Systems[END_REF] and regular power-balanced method [30, chap. 5], (see also the PhDs [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF][START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF][START_REF] Müller | Time-continuous power-balanced simulation of nonlinear audio circuits: realtime processing framework and aliasing rejection[END_REF]).

It is worth mentioning that PHS and quadratisation methods have recently started to be used for musical applications in other laboratories (see [START_REF] Danish | Applications of port hamiltonian methods to non-iterative stable simulations of the korg35 and moog 4-pole vcf[END_REF][START_REF] Ducceschi | Simulation of the geometrically exact nonlinear string via energy quadratisation[END_REF]). See also [START_REF] Mora | Port-Hamiltonian modeling of fluid-structure interactions in a longitudinal domain[END_REF] for results on fluid-structure interactions and voice.

Throughout this document, suggestions for exercises are indicated by the symbol ▷ .

Differential-Algebraic formulation (semi-explicit DAF)

In this formulation, an open multi-physical system is represented by a network of (i) storage components of state x and energy E := H(x) ≥ 0, (ii) memoryless passive components described by an effort law z(w) for flow w, such that the dissipated power is P diss := z(w) ⊺ w ≥ 0 (P diss = 0 for conservative ones);

(iii) external components receiving power P ext := u ⊺ y through ports represented by 1system inputs u and outputs y, (ii) memoryless → algebraic eq.

(iii) ports → physical signals [START_REF] Bernhard | Port-controlled hamiltonian systems: modelling origins and systemtheoretic properties[END_REF] through a skew-symmetric interconnection matrix S possibly depending on (x, w, u).

Then, the powers received by (i-iii) are balanced:

P stored +P diss +P ext = 0, where P stored := Ė (= ∇H(x) ⊺ ẋ ) (2)
denotes the power received by energy-storing components. Indeed, P stored + P diss +

P ext = e ⊺ f (1)
= e ⊺ Se is zero since e ⊺ Se = (e ⊺ Se) ⊺ = -(e ⊺ Se) due to the skew-symmetry of S. This ensures passivity as Ė = -(P diss +P ext ) ≤ -P ext makes the internal energy non-increasing when external sources are switched off (P ext = 0).

= π π 2 /(2m) ẋ1 = π H ′ 1 (x 1 )= x 1 /m sp x 2 := ξ k ξ 2 /2 ẋ2 = ξ H ′ 2 (x 2 )= k
π = F m ξ = V sp V dp V ext      ẋ1 ẋ2 w y      f =       0 -1 -1 -1 1 0 0 0 1 0 0 0 1 0 0 0       S = -S ⊺        H ′ 1 (x 1 ) H ′ 2 (x 2 ) z(w) u        V m = π/m F sp = kξ F dp = rV dp -F ext e
This matches (1) with energy function H(x) = H 1 (x 1 ) + H 2 (x 2 ) and restores the usual ODE with2 z := ξ .

In this example, nonlinear passive oscillators can be obtained by replacing component laws (see H ′ n or z in the table above) satisfying (i-ii). For instance, H 2 (ℓ) =

k 1 ξ 2 2 + k 3 ξ 4 4
with k 1 , k 3 > 0 restores a stiffening Duffing oscillator 3 . Remark 1 (receiver convention). The flows f and efforts e in (1) denote those received by components, whatever their type (i: storing, ii: memoryless, iii: exterior). In physics, this convention is classical for (i-ii), but not systematic for (iii). Note that in example 1, the quantities (iii) are actually those (velocity V ext and force -F ext ) received by the external actuator.

Differential formulation (DF)

A largely-used formulation is the state-space representation

ẋ = f(x, u) := J xx -R xx ∇H(x) + Gu, y = g(x, u) := -G ⊺ ∇H(x), (3)
where matrices

J xx = -J ⊺ xx (skew-symmetric), R xx = R ⊺ xx ⪰ 0 (non-negative symmet- ric)
and G can depend on x, u. This is a particular case of the differential state-space representation ẋ y

=: f D = J xx J xu * J yu =: J D = -J ⊺ D - R xx R xu * R yu =: R D = R ⊺ D ⪰ 0 ∇H(x) u . =: e D (D for Differential formulation) (4)
Remark 2. The case (3) corresponds to matrices

J D = J xx J xu = G * J yu = 0 and R D = R xx R xu = 0 * R yu = 0 ,
involving no feedforward (J yu = R yu = 0) and a conservative routing between ports and energy-storing components (R xu = 0).

Only the flows and efforts of (i) and (iii) are represented, whereas laws of (ii) are converted into connections through matrices J (conservative part) and R (dissipative part), as

P stored + P ext = e ⊺ D f D = e ⊺ D J e D =0 -e ⊺ D R e D ≥0 . (5)
The dissipated power due to (ii) reformulates as

P diss,D := e ⊺ D R e D ≥ 0. It reduces to ∇H(x) ⊺ R xx ∇H(x) for equation (3).
Example 2. The oscillator corresponds to such a case with [START_REF] Duindam | Modeling and control of complex physical systems: the port-Hamiltonian approach[END_REF]. Note that contrary to J (or S in example 1), matrix R is not independent of all the component laws: it mixes the routing (V m → V dp → w) and the reconstruction of the damping force (z(w) = rw) to produce rV m in the force balance equation.

F m V sp V ext    ẋ1 ẋ2 y    f D =      0 -1 -1 1 0 0 1 0 0    J = -J ⊺ -    r 0 0 0 0 0 0 0 0    R = R ⊺ ⪰ 0       H ′ 1 (x 1 ) H ′ 2 (x 2 ) u     V m F sp -F ext e D so that J xx = 0 -1 1 0 , R xx = r 0 0 0 and G = -1 0 in (

Relations between these formulations

A differential formulation (4) simply stems from formulation (1), under the following conditions (met in example 1):

z(w) = Γ(w) w with Γ + Γ ⊺ ⪰ 0, (6a) S ww = 0 and P := [-S ⊺ xw , S wu ] is independent of w. (6b)
In this case, the algebraic part of (1) admits a unique solution

w = W (e D ) := P e D . (7)
Expressing the efforts of memoryless components through the function composition4 z • W leads to (4) with

J D = S xx S xu * S yu -P ⊺ J Γ P with J Γ := 1 2 (Γ D -Γ ⊺ D ), (8a)
R D = P ⊺ R Γ P ⪰ 0 with R Γ := 1 2 (Γ D + Γ ⊺ D ), (8b)
and where

Γ D = Γ • W can depend on 5 e D .
Example 3. In example 1, Γ = r (J Γ = 0 and R Γ = r) and P = [1 0|0] yield the result stated in example 2.

▷You are now ready for exercice 1 (questions Q1 and Q2, see § 5).

Hybrid formulation (HF)

A reformulation can also be derived by eliminating a part A of components (w A → z A (w A , w B )) for which (6a-6b) are met, while the complementary part B (w

B → z B (w A , w B ))
is left untreated. Renaming w B and z B without label B for sake of conciseness, this case leads to the hybrid formulation

   ẋ w y    =: f h =      J xx J xw J xu * J ww J wu * * J yu    =: J h = -J ⊺ h -    R xx R xw R xu * R ww R wu * * R yu    =: R h = R ⊺ h ⪰ 0      ∇H(x) z(w) u    =: e h (h for hybrid) , (9)
that is, f h = M e h with M := J h -R h , for which the dissipated power corresponds to

P diss,h = z(w) ⊺ w ≥ 0 + e ⊺ h R h e h ≥ 0 ≥ 0. ( 10
)
Remark 3 (Link with formulation (1)). Formally, this hybrid formulation (9) can be viewed as [START_REF] Bernhard | Port-controlled hamiltonian systems: modelling origins and systemtheoretic properties[END_REF], in which S ≡ J h -R h is no longer skew-symmetric but includes a nonpositive symmetric part (-R h ) partly responsible for dissipation.

Interconnection of several PHS

Consider a collection of N independent systems S n=1,...,N , all described by [START_REF] Bernhard | Port-controlled hamiltonian systems: modelling origins and systemtheoretic properties[END_REF]. Denote their gathered variables

v ⊺ := [v ⊺ S 1 , . . . , v ⊺ S N ] ⊺ , for all labels v = x, w, u or y,
and the associated block-diagonal matrices

S v ′ v = diag(S v ′ S 1 v S 1 , . . . , S v ′ S N v S N
) for all labels v ′ and v, to form the description (1) of the global PHS S . Now, interconnect these systems according to a causal6 memoryless passive 7 law

u c = z c (y c ), (interconnection law) (11)
between some ports of S ≡ (S 1 , . . . , S N ), selected as

[y c , u c ] = C ⊺ [y, u], (connected ports) (12a) [y f , u f ] = F ⊺ [y, u], ( free ports) (12b)
where the concatenation of selection matrices [C T , F ⊺ ] defines a permutation matrix that rearranges the ports. Then, the interconnected global system is a PHS governed by

     x w y c y f      =      S xx S xw S xu C S xu F * S ww S wu C S wu F * * C ⊺ S yu C C ⊺ S yu F * * * F ⊺ S yu F      skew-symmetric      ∇H(x) z(w) z c (y c ) u f      . ( 13
)
Under condition (6a-6b) adapted to8 z c , the connection variable y c may be fully (otherwise, partially) reduced following the processes described in sections 2.3-2.4. In this case, z c defines a fully (otherwise, partially) causal interconnection. Moreover, in accordance with remark 3, formulation (13) can be adapted to the hybrid formulation.

Remark 4 (PyPHS). In [START_REF] Falaize | Passive Guaranteed Simulation of Analog Audio Circuits: A Port-Hamiltonian Approach[END_REF], a method which allows the automatic derivation of a complete model from components and interconnections is described for electronic circuits. A corresponding Python library PyPHS [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF]77,[START_REF] Falaize | Simulation of finite-dimensional multi-physical systems described by networks of components[END_REF] is available for multi-physical systems. It automates the derivation of the PHS, of the discrete gradient-based numerical scheme (presented below) with a Newton-Raphson solver, and it generates C++ code for simulation. It is aimed at numerical simulation, such as 20-sim based on bond graphs (see also the recent work [START_REF] Pfeifer | Automated generation of explicit port-hamiltonian models from multi-bond graphs[END_REF][START_REF] Pfeifer | Automated Model Generation and Observer Design for Interconnected Systems : a Port-Hamiltonian Approach[END_REF]), Modelica/Dymola. See also work on dedicated languages [START_REF] Lohmayer | Ephs: A port-hamiltonian modelling language[END_REF].

Elementary non-causal interconnections

Some systems involve non-causal interconnections (no relation f = Se is available).

Examples 9 are components that share the same effort (constraints of type [e] i = [e] j ), or, when not caused by the complementary part of the system, the same flow

([f] i = [f] j ).
The general case can be addressed by considering kernel representations of Dirac structures 10 (see e.g. [START_REF] Van | Port-hamiltonian systems: an introductory survey[END_REF]) or specific differential-algebraic formulations (see e.g. [START_REF] Van Der Schaft | Generalized port-hamiltonian dae systems[END_REF][START_REF] Van Der Schaft | Dirac and lagrange algebraic constraints in nonlinear port-hamiltonian systems[END_REF]). Two types of elementary (but common) interconnections of energy-storing components are detailed below, that lead to an equivalent component (single Hamiltonian).

Shared efforts

Hamiltonians H 1 and H 2 are supposed to be C 1 convex non-negative functions, so that their gradients are continuous monotone. The two components share the same effort

e := E 1 (x 1 ) = E 2 (x 2 ), with E i := ∇H i of inverse X i := E -1 i : e -→ x i . ( 14
)
The set of these components receives the total flow f = ẋ1 + ẋ2 . This extensivity property on the flow makes f interpretable as the time derivative of a total state

x := x 1 + x 2 = X(e), with X := X 1 + X 2 monotone of inverse E := X -1 : x → e. (15)
Finally, expressing the total energy H 1 (x 1 ) + H 2 (x 2 ) as a function of the total state x provides the Hamiltonian of the equivalent component

H : x ∈ R N -→ [H 1 • X 1 + H 2 • X 2 ] • E(x). (16)
Note that H is convex (sum of convex functions), continuous (functions in ( 16) are all continuous) and such that H 1 (0) = H 2 (0) = 0 ⇒ H(0) = 0. Its C 1 -regularity (expected to be inherited from H 1 , H 2 ) is also satisfied (a proof not detailed here, also valid for C k -cases, is based on the use of push-forward measures).

Shared flows

In the case of energy-storing components that share a flow, the two components with Hamiltonian H 1 (x 1 ) and H 2 (x 2 ) receive the same flow ẋ1 = ẋ2 , so that x 1 (t) -x 1 (0) = x 2 (t) -x 2 (0). An equivalent component can then be defined by a Hamiltonian

H X 0 : x ∈ R N -→ H 1 (x) + H 2 (x + X 0 ), with X 0 := x 2 (t = 0) -x 1 (t = 0), (17)
which depends on the initial condition through X 0 .

Note that hamiltonian H is chosen such that H X 0 (0) = 0. It admits a lower bound, possibly negative. Another choice for defining its non-negative version is H + X 0 (x) := H X 0 (x)min x H X 0 (x). For quadratic hamilonians, these two choices coincide.

Note also that, as mentioned above, when the flow f is caused by the complementary part of the system, this interconnection is not to be processed as a constraint (writing ẋ1 = f, ẋ2 = f), but [START_REF] Kotyczka | Discrete-time port-hamiltonian systems and control[END_REF] can still be used for order reduction.

PHS shifting

When the time-varying variables of a system have a bias 11 (non-centered signals), it can be beneficial, for e.g. analysis or numerical issues, to formulate the exact equations that govern their centered fluctuating parts. The shifting operation is the transformation that produces these equations from the formulation governing the original variables.

Consider a PHS described by the hybrid formulation [START_REF] Pfeifer | Automated generation of explicit port-hamiltonian models from multi-bond graphs[END_REF] for sake of generality 12 . Assume that M := J h -R h is constant. Decompose all the original variables (var(t)) into the sum of a bias value (var ⋆ ) and a fluctuation ( var(t)). Variables var ⋆ correspond to equilibrium values, as they coincides with the original variables when the system is at rest (no fluctuation). These variables are related as:

(PHS):    ẋ w y    f h (t) = M    ∇H(x) z(w) u    e h (t)
and (PHS) ⋆ :

   ẋ⋆ = 0 w ⋆ y ⋆    f ⋆ h = M    ∇H(x ⋆ ) z(w ⋆ ) u ⋆    e ⋆ h . (18)
The shifted system (PHS)

≡ (PHS)-(PHS) ⋆ reads f h = M e h ≡ f h (t)-f ⋆ h = M (e h (t)-e ⋆ h
). Writing efforts w.r.t. var ⋆ and var leads to ∇H( x + x ⋆ ) -∇H(x ⋆ ) for energy-storing components (i), z( w+w ⋆ ) -z(w ⋆ ) for memoryless components (ii), and u for ports (iii). (PHS) must involve functions x ⋆ → H x ⋆ ( x) and w ⋆ → z w ⋆ ( w), in which x ⋆ and w ⋆ are considered as parameters. Integrating expression (i) from 0 w.r.t. x defines H x ⋆ ( x), leading to (PHS):

   ẋ w y    f h (t) = M    ∇ H x ⋆ ( x) z w ⋆ ( w) u    e h (t) with H x ⋆ ( x) := H( x + x ⋆ ) -∇H(x ⋆ ) ⊺ x -H(x ⋆ ), z w ⋆ ( x) := z( w + w ⋆ ) -z(w ⋆ ). (19)
Remark 5 (Passivity). The shifted system (PHS) is passive w.r.t. input u and output y, if z w ⋆ ( w) ⊺ w ≥ 0.

Remark 6 (Bregman distance). As noticed in [14, Rk 6.5.2, p. 137], the shifted Hamiltonian can be interpreted as

H x ⋆ ( x) := D H (x ⋆ + x, x ⋆ ) with D H (x, x ⋆ ) := H(x) -H(x ⋆ ) -∇H(x ⋆ ) ⊺ (x -x ⋆ ), ( 20
)
where D H defines the Bregman distance for convex functions H (see also [START_REF] Jayawardhana | Passivity of nonlinear incremental systems: Application to pi stabilization of nonlinear rlc circuits[END_REF][START_REF] De | Bregman storage functions for microgrid control[END_REF]).

3 Discrete gradient method

Objective and principle

This section presents a numerical method to compute discrete-time trajectories of PHS, which preserves passivity, that is, for which (2) is satisfied 13 on a time grid t ∈ (δt)Z, and in which P stored still corresponds to a time-variation of energy E in discrete-time with consistent sign (including zeroness for conservative systems). Such an adequate and simple discrete-time approximation of the energy variation is

P stored (t) ≈ [ E(t + δt) -E(t) ] / δt, (21)
so that the energy still increases or decreases according to the sign of -(P diss +P ext ).

Introducing the discrete derivation operator D as

∀f ∈ C 1 (R, R N ), D f (s, δ ) := f (s+δ )-f (s) δ if δ ̸ = 0, f ′ (s) otherwise, (22)
this choice consists in considering the discrete variation D E(t, δt) instead of its continuous limit case D E(t, 0). The objective is to propose a numerical method to simulate the trajectories governed by (1, 4 or 9) and that satisfies this principle.

Case I: mono-variate Hamiltonian

For a mono-variate Hamiltonian H : x ∈ R → R + , the time-variation D E(t, δt) of energy E(t) := H • x(t) is given by

D[H • x](t, δt) = H( x(t + δt) ) -H( x(t) ) δt = D H x(t), x(t + δt) -x(t) x(t + δt) -x(t) δt = D H x(t), D x(t, δt) δt D x(t, δt). (23)
This formula generalises the chain rule for derivatives ( Ė = H ′ (x(t)) ẋ(t)) to operator D. Therefore, jointly replacing ẋ(t) and ∇H(x(t)) in ( 1) by their approximations D x(t, δt) and D H x(t), D x(t, δt) δt defines a family (continuously parameterised by δt > 0) of discrete-time PHS with the expected power balance. A lighter writing is obtained by denoting

E[k], x[k] (etc) the signals sampled at time t k := k δt (k ∈ Z or, given an initial condition x[0], k ∈ N), and δ x[k] the increment x[k + 1] -x[k] = D x(t k , δt k )δt.
The substitutions to operate in (1, 4 or 9) then rewrite

(I) ẋ(t) → δ x[k]/δt ∇H(x(t)) → D H(x[k], δ x[k]) with δ x[k] := x[k + 1] -x[k]. (24)

Case II: Hamiltonian with separated variables

Assuming that

H(x = [x 1 , . . . , x N ] ⊺ ) (II) := ∑ N n=1 H n (x n ), equation (23) becomes D[H • x](t, δt) (22)
=

N ∑ n=1 H n ( x n (t + δt) ) -H n ( x n (t) ) δt (23) = N ∑ n=1 D H n x n (t), D x n (t, δt) δt D x n (t, δt), (25)
rewritten as

D[H • x](t, δt) = ∇ D H x(t), D x(t, δt) δt ⊺ D x(t, δt), (26)
and the substitutions (24) become vectorial

(II or III) ẋ(t) → δ x[k]/δt ∇H(x(t)) → ∇ D H(x[k], δ x[k]) with δ x[k] := x[k + 1] -x[k], (27)
where operator ∇ D H can be defined from [START_REF] Yang | Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method[END_REF][START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF] by

∇ D H = ∇ II H with ∇ II H(x, δ x) :=     D H 1 (x 1 , δ x 1 ) . . . D H N (x N , δ x N )     under assumption (II). (28)
In general, any operator ∇ D satisfying ( 26) is called a discrete gradient.

▷ Exercise 1 (Q3a)

Case III: general case

In the case III for which H(x = [x 1 , . . . , x N ] ⊺ ) is any multivariate function (II is no longer assumed), equation ( 26) remains valid for the discrete gradient

∇ III H(x, δ x) :=     D[σ1H](x1, δ x 1 ) . . . D[σNH](xN, δ x N )     with σ n H : ξ → H        x 1 . . . x n-1 ξ x n+1 + δ x n+1 . . . x N + δ x N        . ( 29
) Indeed, since [σ n H](x n + δ x n ) ( * ) = [σ n-1 H](x n-1 ) for 2 ≤ n ≤ N, it follows that ∇ III H(x, δ x) ⊺ δ x = D[σ1H](x1, δ x 1 ) , . . . , D[σNH](xN, δ x N ) δ x (22) = N ∑ n=1 σ n H(x n + δ x n ) -σ n H(x n ) = σ 1 H(x 1 + δ x 1 ) + N ∑ n=2 σ n H(x n + δ x n ) - N-1 ∑ n=1 σ n H(x n ) -σ N H(x N ) ( * ) = σ 1 H(x 1 + δ x 1 ) -σ N H(x N ) = H(x + δ x) -H(x), (30)
from which (26) stems. Then, the substitution [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF] with ∇ D H = ∇ III H provides a powerbalanced numerical method.

Symmetric discrete gradient

For all permutations π ∈ P(N) (permutations on {1, 2, . . . , N}), consider the transformation that reorders the state variables and denote for simplicity

π(x = [x 1 , . . . , x N ] ⊺ ) :≡ [x π(1) , . . . , x π(N) ] ⊺ (without ambiguity), (31a)
and

H π := H • π -1 , (31b)
so that H π ( π(x) ) = H(x) restores the Hamiltonian for the reordered state variables π(x). From (31b), the n-th component of

∇H π ( π(x) ) is the π -1 (n)-th component of ∇H(x), that is, π • ∇H π • π = ∇H.
Technically, this means that the gradient ∇ is invariant under the transformation

T π : ∇ -→ (H → π • ∇H π • π), (32)
that is, T π (∇) = ∇. Physically, this means that reordering the state variables does not modify the effort laws, which consistently translates the physical principle of invariance to any arbitrary declaration order of the state variables.

To satisfy this physical principle, a discrete gradient ∇ D H(x, δ x) must satisfy the invariance w.r.t. to the joint reordering of (x, δ x), that is,

T π D (∇ D ) = ∇ D for the adapted transformation T π D : ∇ D -→ ( H → π • [∇ D H π ](π • • , π • •) ). (33)
As the sequence (σ π(n) ) n=1,...,N differs from (σ n ) n=1,...,N , ∇ III does not satisfy this invariance: for all π ∈ P(N), T π D (∇ III ) defines a discrete gradient 14 denoted

∇ π III below, but T π D (∇ III ) ̸ = ∇ III if π ̸ = I d .
For all discrete gradients ∇ D that apply to multivariate Hamiltonians, the symmetric discrete gradient ∇ D = ∇ D defined by

∇ D -→ ∇ D := 1 N! ∑ π∈P(N) T π D (∇ D ) (symmetric discrete gradient), (34)
restores the invariance of the original PHS to any arbitrary declaration order of the state variables.

Example 4 (Non-equivalent discrete gradients). Consider the (dimensionless) Hamiltonian H(x

:= [x 1 , x 2 ] ⊺ ) = 1 2 x 2 1 (1+εx 2
2 ) and denote π the swapping application (π(1) = 2, π(2) = 1). Deriving the discrete gradients yields 30) is satisfied replacing (H, x) by (H π , π • x) and leaving σ n unchanged.

∇ III H ∇ π III H := T π D ∇ III H ∇ III H   (x 1 + δ x 1 2 ) 1 + ε(x 2 + δ x 2 ) 2 ) x 2 1 2 ε (2x 2 + δ x 2 )   (x 1 + δ x 1 2 )(1 + εx 2 2 ) (x 1 +δ x 1 ) 2 2 ε (2x 2 + δ x 2 )   (x 1 + δ x 1 2 ) 1 + ε x 2 2 +(x 2 +δ x 2 ) 2 2 1 2 x 2 1 +(x 1 +δ x 1 ) 2 2 ε (2x 2 + δ x 2 )   14 (
Remark 7 (Midpoint discrete gradient). The operator given by

∇ mid H(x, δ x) :=    ∇H(x+ 1 2 δ x)+ H(x+δ x)-H(x)-∇H(x+ 1 2 δ x) ⊺ δ x δ x ⊺ δ x δ x, if δ x ̸ = 0, ∇H(x) otherwise, (35)
defines a discrete gradient which is naturally symmetric (see the lecture notes [17, (26-27)]).

Remark 8 (Equivalent discrete gradients). If a Hamiltonian H satisfies assumption II, then all its discrete gradients coincide. This property is also satisfied for quadratic Hamiltonians, even with non-separated variables (see remark 9 below). In summary:

∇ III H = ∇ π III H = ∇ III H = ∇ mid H, under assumption (II) or if H is quadratic. As a corollary, if H is a sum of K functions, namely, H(x = [x 1 , . . . , x N 1 +•••+N K ] ⊺ ) = ∑ K k=1 H k ([x n k +1 , . . . , x n k +N k ] ⊺ ) (with indexes n 1 = 0, n k+1 = n k + N k and N k ≥ 1)
, then the average in [START_REF] Lopes | Energy Balanced Model of a Jet Interacting With a Brass Player's Lip[END_REF] can be restricted to permutations operating separately inside the K subgroups of state variables. This advantageously involves N 1 ! . . . N K ! rather than

(N 1 + • • • + N K )! terms.
Remark 9 (Quadratic Hamiltonian and midpoint rule). For linear systems, the discrete gradient method corresponds to the midpoint rule. Indeed, for quadratic Hamiltonians, the continuous and discrete gradients are given by

H(x)= 1 2 x ⊺ Lx with L=L ⊺ ≻0 =⇒ ∇H(x)=Lx and ∇ D H(x, δ x)=L(x+ 1 2 δ x), (36)
so that ∇ D H(x, δ x) = ∇H(x mid ) with x mid = x + 1 2 δ x. This is not the case in general.

The definitions and properties of the various discrete gradients are summarised below.

Case on H non symmetric

∇ D symmetric ∇ D II ∇ II = ∇ III = ∇ π III = ∇ III = ∇ mid III-quadratic × ∇ III = ∇ π III = ∇ III = ∇ mid (≡ midpoint rule) III-general × ∇ III ̸ = ∇ π III ̸ = ∇ III ̸ = ∇ mid (all distinct) Definition (28) (29) (29,33) (29,34) (35)
3.6 Method (implicit numerical scheme)

Consider a PHS described by the hybrid formulation (9) for sake of generality 15 , where

M [k] := M h X[k], W[k], U[k] , with M := J h -R h (37a)
is evaluated at samples k, choosing the basic (causal non time-symmetric) estimates

X[k] := x[k], W[k] := w[k], U[k] := u[k], (37b)
that do not involve any time-increments for simplicity (see remarks 10-11 for comments on properties). Then, for a discrete gradient ∇ D , the discrete-time system described by, for all k ∈ Z (or, given an initial condition x[0], for all k ∈ N),

   δ x[k]/δt w[k] y[k]    = M [k]      ∇ D H x[k] , δ x[k] z w[k] u[k]      , (38a)
x

[k + 1] = x[k] + δ x[k], (38b)
defines a discrete-time PHS that approximates (9) and fulfills the discrete power balance

(E[k + 1] -E[k])/δt P stored [k] + z(w[k]) ⊺ w[k] + e h [k] ⊺ R h [k] e h [k] P diss [k] ≥ 0 + u[k] ⊺ y[k] P ext [k] = 0. (39)
Equations (38a-38b) define an implicit numerical scheme, since (38a) must be solved w.r.t. (δ x[k], w[k]) at each sample k. The next section is devoted to avoid the use of (fixed-point, Newton-Raphson, etc.) iterative solvers.

Remark 10 (Accuracy). The accuracy order of this method is shown to be 1 in general, and 2 if w is empty and M is constant. Methods based on the discrete gradient, of order higher than 1, can be built using more refined estimates than (37b) that involve increments (implicit estimates, see e.g. [17, § 3]) or that involve multiples stages (with possibly non implicit estimates, see [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF][START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF][START_REF] Lopes | Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems[END_REF]).

Remark 11 (About properties and features). In the context of PHS, important properties and features of numerical schemes are: (a) preserving passivity and the power balance (structured into conservative, dissipative/irreverible and external parts), (b) the accuracy order (error consistency, error stability, convergence), and possible additional features according to some objectives such as, for audio-acoustics, (c) the computation cost (with real time application perspectives), (d) the rejection of aliasing (artifact naturally produced when combining nonlinearities and sampling), or, some general properties such as (see the lecture notes [17, § 3]) (e) simplecticity (area/hypervolume conservation in the phase space) of Hamiltonian flow16 , (f) preserving timereversal symmetries (see J in [17, § 3.1]), etc. The method (38a-38b) achieves properties (a) (see [START_REF] Falaize | Passive modelling of the electrodynamic loudspeaker: from the Thiele-Small model to nonlinear port-Hamiltonian systems[END_REF]), (b) with a low but improvable order (remark 10) and prepares (c) ( § 4). Extensions to (c,d) with controllable orders are mentioned in the conclusion.

Concluding remarks

Work on structure-preserving numerical methods is extensive [START_REF] Haier | Geometric Numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF]. Energy-preserving methods based on discrete gradient [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] or average-vector field [START_REF] Quispel | A new class of energy-preserving numerical integration methods[END_REF] in order to derive discrete-time PHS have been used in several works (see e.g. [START_REF] Aoues | Schémas d'intégration dédiés à l'étude, l'analyse et la synthèse dans le formalisme Hamiltonien à ports[END_REF][START_REF] Aoues | Hamiltonian systems discrete-time approximation: Losslessness, passivity and composability[END_REF] or e.g. [24, see also other papers of this author] for high order accuracy). Complementary references can be found in the lecture notes [START_REF] Kotyczka | Discrete-time port-hamiltonian systems and control[END_REF]. Recent work introducing RPM (for Regular Power-Balanced Methods) [30, chap. 5], which deals with the properties (a-d,f) with time-continuous frames, should also be mentioned.

4 Quadratisation method and non-iterative solver

Consider the problem (38a-38b) restricted to the differential formulation (w is empty) as ( 4) for simplicity of presentation: omitting [k] for conciseness, (38a) is replaced by,

δ x /δt y = M xx M xu M yx M yu   ∇ D H x , δ x u   with M = J -R. (40)

Objective and principle

This section presents a numerical method with conditions under which the implicit numerical scheme [START_REF] Matignon | On damping models preserving the eigenfunctions of conservative systems: a port-Hamiltonian perspective[END_REF] admits an explicit solution.

The objective is to avoid the use of iterative solvers such as the Newton-Raphson or the fixed-point algorithms. Indeed, for non-sparse low-dimensional problems, this alternative can be useful to reduce the computational cost and guarantee real-time for some audio-acoustics applications.

The principle introduced in [START_REF] Lopes | Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems[END_REF] relies on ( § 4.2) the existence of such a solution for linear systems (quadratic hamiltonian) and ( § 4.3) the introduction of a change of variables that ( § 4.4) quadratises the hamiltonian while preserving the skew-symmetric (resp., non-negative) nature of matrix J (resp., R), for a large class of nonlinear systems.

Case of a quadratic hamiltonian

Assume that

H(x) = 1 2 x ⊺ Lx with L = L ⊺ ≻ 0. ( 41a
)
Then (see remarks 8-9), all its discrete gradients coincide and ∇ D H(x, δ x) = L(x +

is invertible, since ∆ is the sum of a positive matrix (I + δt 2 R xx ≻ 0) and a skewsymmetric one (-δt 2 J xx ). Equations (41b-41c) provide a non-iterative process to compute the increment δ x at each sample k, from which are computed the output (lower line in [START_REF] Matignon | On damping models preserving the eigenfunctions of conservative systems: a port-Hamiltonian perspective[END_REF]) and the state at sample k + 1 (updating equation (38b)).

Remark 12 (Effort equation). For the effort (e x = Lx), the equation on increment δ e x is the same as (41b) with matrices A e := M xx , B e := M xu and ∆ e := L -1 -δt 2 A e . This equation may be preferred as it potentially requires fewer flops, L being constant.

▷ Exercise 1 (Q3b -Q5).

Change of state (continuous time domain)

Consider a C 1 -regular bijective change of state 17 Q and the Hamiltonian H q Q : x -→ q of inverse denoted X := Q -1 : q -→ x, (42a)

H q : q -→ H( x = X(q) ), (42b)
and denote the (invertible) associated Jacobian matrix (from x to q)

J qx (x) :=∇ ⊺ x Q (x) =     ∂ x 1 Q 1 (x) . . . ∂ x N Q 1 (x) . . . . . . . . . ∂ x 1 Q N (x) . . . ∂ x N Q N (x)     . (42c)
Then, the following differential formulations represent the same PHS:

ẋ y = M (x, u) ∇H(x) u ⇐⇒ q y = M q (q, u) ∇H q (q) u with M q := QM , (43)
where the linear transformation Q defined by

[QM ](q, u) := J q (q) M (X(q), u) J q (q) ⊺ with J q (q) :=

J qx • X(q) 0 0 I dim(u) , (44)
preserves the skew-symmetry, the symmetry and the positivity of matrices, so that

M = J -R, J = -J ⊺ , R = R ⊺ ⪰ 0 Q -→ M q = J q -R q , J q = -J ⊺ q , R q = R ⊺ q ⪰ 0.
The proof stems from the mappings between the pairs of states, flows and efforts:

For all times, these pairs are related as

(states) q = Q(x) ( ⇔ x = X(q) ), (45) (flows)
q = J qx (x) ẋ, (46)
(efforts) ∇H(x) = ∇ x [H q • Q](x) = J qx (x) ⊺ [∇ q H q ](Q(x)). (47)
Expressing the structure of M as in (4) and its dependency w.r.t. (x, u), it comes

q (46,4) = J qx (x) M xx (x, u)∇H(x) + M xu (x, u)u (47) = J qx (x) M xx (x, u)J qx (x) ⊺ [∇ q H q ](Q(x)) + M xu (x, u)u (45) = J qx (X(q)) M xx (X(q), u)J qx (X(q)) ⊺ [∇ q H q ](q) + M xu (X(q), u)u , y (47,45)
= M yx (X(q), u)J qx (X(q)) ⊺ [∇ q H q ](q) + M yu (X(q), u)u, yielding [START_REF] Lopes | Modèle d'interaction Jet/Lèvre préservant le bilan de puissance pour les instruments de type cuivre[END_REF][START_REF] Falaize | Passive simulation of electrodynamic loudspeakers for guitar amplifiers: a port-Hamiltonian approach[END_REF]. The preservation of the symmetries and of the spectrum sign of a matrix K through Q stems from its form in J q K J ⊺ q with J q invertible.

Quadratisation method and solver

A port-Hamiltonian system (4) with a C 1 -regular hamiltonian H(x) is said to be quadratisable if it exists a C 1 -regular bijective change of state Q = X -1 such that

H q (q) := H • X(q) = 1 2 q T Lq, with L ≻ 0. (48)
Then, the q-representation of the PHS (right-hand side of equivalence ( 43)) has a quadratic hamiltonian so that the non-iterative solver proposed in section 4.2 applies to this q-representation

q y = M q (q, u) Lq u with M q := J q -R q , J ⊺ q = J q , R q = R ⊺ q ⪰ 0. ( 49
)
An important corollary is the existence and uniqueness of the discrete-time solution.

The quadratisation method can be summarised in the following sequence of steps.

Procedure:

Given a PHS under formulation (4) with hamiltonian H(x):

1. Find Q such that (48) is satisfied, 2. Build (=derive or implement) function X := Q -1 , 3. Build matrices M q := J q -R q , J ⊺ q = J q , R q = R ⊺ q as functions of (q, u), by applying transformation Q (see [START_REF] Falaize | Passive simulation of electrodynamic loudspeakers for guitar amplifiers: a port-Hamiltonian approach[END_REF]), 4. Build matrices A q , B q and ∆ q as functions of (q, u) for the q-representation (49) following § 4.2 (adapting labels q instead of x in (41c)), 5. Build the incremental state function q[k] from (41b), 6. Build the updating equation q[k+1] = q[k]+δ q[k] and the observation equation on y.

→ Your input-output power-balanced simulation is (almost) ready! → The state x can also be evaluated with x = X(q). Remark 13 (Units). In this definition, matrix L can be chosen such that q has the same physical units as x (or other properly-chosen physical units). It can also be chosen, without loss of generality, as the dimensionless identity matrix, meaning that the units of q are square-roots of Joules [ √ J = √ Kg.m.s -1 ].

5 Exercises Observe (a posteriori) that you could find this result directly from the circuit: the interest of the procedure (a,b,c) is its systematic property.

v G i G L 0 v L = φ i L R 0 v R i R C 0 v C q = i c
Question 3. Discrete gradient and numerical method. (a) Derive the discrete gradient associated with the Hamiltonian H using [START_REF] Lebrun | Modélisation multi-physique passive, identification, simulation, correction et asservissement de haut-parleur sur des comportements cibles[END_REF]. (b) Observe that it also satisfies the property [START_REF] Najnudel | Simulation of the Ondes Martenot Ribbon-Controlled Oscillator Using Energy-Balanced Modeling of Nonlinear Time-Varying Electronic Components[END_REF] in remark 9 (valid for all quadratic Hamiltonians), so that the method in § 4.2 applies.

Exercise 3 (nonlinear circuit): PHS shifting

(This exercise can be done after question 6)

Consider that the generator has a bias: v G (t) = v ⋆ G + v G (t) Question 9. Write the equations satisfied by x ⋆ . Is the solution unique ?

Question 10. Derive the shifted Hamiltonian H x ⋆ ( x) (see [START_REF] Itoh | Hamiltonian-conserving discrete canonical equations based on variational difference quotients[END_REF]).

Question 11. Can the shifted-PHS be quadratised ? Propose an implementation.

Exercise 4 (pyPHS)

The Python library pyPHS provides tools to automate the derivation of the PHS (in L A T E X), of the discrete gradient-based numerical scheme with a Newton-Raphson solver, and to generate C++ code for simulation. After this lecture, the reader is invited to follow the presentation, tutorials and examples (see link indications in [77]).

theS

  internal flows f(t) and efforts e(t) of which are coupled according to the following differential-algebraic equation xx S xw S xu * S ww S wu * storage → differential eq.

1 2

 1 δ x). The implicit equation to solve w.r.t. δ x (upper line in[START_REF] Matignon | On damping models preserving the eigenfunctions of conservative systems: a port-Hamiltonian perspective[END_REF]) rewrites δ x/δt = A(x+ 1 2 δ x) + Bu so thatδ x/δt =∆ -1 (Ax + Bu),(41b)with A :=M xx L, B := M xu , and where ∆ := Iδt 2 A

Figure 1 :Question 1 .

 11 Figure 1: Circuit and notations for exercises

  Differential formulation (DF) . . . . . . . . . . . . . . . . . . . . . . 2.3 Relations between these formulations . . . . . . . . . . . . . . . . . 2.4 Hybrid formulation (HF) . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Interconnection of several PHS . . . . . . . . . . . . . . . . . . . . . 2.6 Elementary non-causal interconnections . . . . . . . . . . . . . . . . 2.6.1 Shared efforts . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6.2 Shared flows . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 PHS shifting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	Contents 2 Input-state-output representations of port-Hamiltonian
	1 Introduction systems
	This section presents standard structured state-space representations of PHS in the 2 Input-state-output representations of port-Hamiltonian systems finite-dimensional case, according to several formulations: ( § 2.1) a component-based 2.1 Differential-Algebraic formulation (semi-explicit DAF) . . . . . . . . semi-explicit differential-algebraic formulation, ( § 2.2) a commonly used differential formulation. Then, relations between them ( § 2.3- § 2.4), interconnections of PHS ( § 2.5-2.2 References 2.6) and PHS shifting ( § 2.7).

3 Discrete gradient method 3.1 Objective and principle . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Case I: mono-variate Hamiltonian . . . . . . . . . . . . . . . . . . . 3.3 Case II: Hamiltonian with separated variables . . . . . . . . . . . . . 3.4 Case III: general case . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Symmetric discrete gradient . . . . . . . . . . . . . . . . . . . . . . 3.6 Method (implicit numerical scheme) . . . . . . . . . . . . . . . . . . 3.7 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Quadratisation method and non-iterative solver 4.1 Objective and principle . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Case of a quadratic hamiltonian . . . . . . . . . . . . . . . . . . . . 4.3 Change of state (continuous time domain) . . . . . . . . . . . . . . . 4.4 Quadratisation method and solver . . . . . . . . . . . . . . . . . . . 5 Exercises 5.1 Exercise 1 (linear circuit): formulations and power-balanced simulation 5.2 Exercise 2 (nonlinear circuit): quadratisation method . . . . . . . . . 5.3 Exercise 3 (nonlinear circuit): PHS shifting . . . . . . . . . . . . . . 5.4 Exercise 4 (pyPHS) . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Example 1 (damped mechanical oscillator). The oscillator governed by mz + r ż + kz = F ext can be formulated as

	• 4 separate components: (i) a mass m of momentum π, a spring sp of elonga-
	tion ξ ; (ii) a damper dp of velocity V dp ; (iii) an external actuator ext applying
	a force F ext :
	state	energy H n flow f	effort e
	m x 1 :

  +F sp +F dp +(-F ext )=0) and all velocities are equal (V m =V sp =V dp =V ext ), leading to

					x 2
	dp	blue : force	w := V dp	z(w) := r w
	ext	red : velocity	y := V ext u	:= -F ext
	• assembled with rigid connections, meaning that internal forces are balanced
	(F m

Variables (u, y) := ±(u, -y) are sometimes a preferred convention so that their product (= -P ext ) reports power supplied by the outside (see also remark 1).

Through the connections, the velocities coincide (V sp = V m = V dp = V ext = ż) and through the component laws, it comes z = ξ and mz = π.

See Ref.[START_REF] Jossic | Energy shaping of a softening Duffing oscillator using the formalism of Port-Hamiltonian Systems[END_REF] for a passive softening case.

As in example 2, this mixes a routing P and a component law z.

defining a PHS under an integrability condition on J D[4, p.48] or a so-called pseudo-PHS otherwise.

law feeding inputs of S from observed outputs of S .

as in § 2.1 (ii)

The condition is C ⊺ S yu C = 0 and that P = C ⊺ -S xu , -S wu , +S yu F does not depend on w c .

Shared flow: capacitors in series, coils in parallel, cascaded springs, etc. Shared effort: capacitors in parallel, coils in series, springs attached together at both ends, etc.

Efforts and flows are related as Kf + Le = 0 where K, L are N × N interconnection matrices, which satisfy KL ⊺ + LK ⊺ = 0 and rank[KL] = N. Matrix S = -K -1 L exists if K is invertible.

This occurs in electronics (operating point configuration through a power supply), electromagnetism (non-zero average charge, or magnets), in mechanics (pre-stress), in fluid mechanics, thermodynamics and acoustics (average density, velocity, temperature/entropy or atmospheric pressure), etc.

in the generic sense that w can be empty or R h can be zero.

and by extension (5) and[START_REF] Pfeifer | Automated Model Generation and Observer Design for Interconnected Systems : a Port-Hamiltonian Approach[END_REF]

in the generic sense that w can be empty or R h can be zero.

this property can be combined with (a), accepting to modify the original Hamiltonian.

We consider here diffeomorphisms Q : R N → R N .

Note that the receiver convention is used for all dipoles (currents and voltages are oppositely oriented, including for the voltage generator).

Acknowledgements

Note: some elements of section 2 are based on an article in preparation.

Acknowledgements

The author thanks (by alphabetical order) Ivan Cohen, Jean-Baptiste Dakeyo, Brigitte d'Andréa-Novel, Antoine Deschamps, Antoine Falaize, Thibault Geoffroy, Thomas Guennoc, Marguerite Jossic, Tristan Lebrun, Nicolas Lopes, Bernhard Maschke, Denis Matignon, Rémy Müller, Judy Najnudel, Nicolas Papazoglou, Mathis Raibaud, David Roze, Fabrice Silva, Tarik Usciati, Colette Voisembert, Victor Wetzel and Marc Wijnand for their collaborative work and publications on Port-Hamiltonian Systems within the S3AM team at STMS laboratory (IRCAM-CNRS-SU).

Software [77] Antoine Falaize and Thomas Hélie. PyPHS: Passive modeling and simulation in python, 2016. Python package on github: -https://github.com/pyphs/pyphs Complementary web pages: -https://pyphs.github.io/pyphs/ (more information on the algebraic-differential structure and the numerical method) -https://github.com/pyphs/pyphs/tree/master/pyphs/examples (many examples) -https://medias.ircam.fr/x152ab3 (presentation video in French): -https://medias.ircam.fr/x4d31af (tutorial video in French). Patent [78] Thomas Hélie and Tristan Lebrun. Dispositif de génération d'un signal de commande d'un système électrique. Patent (Déclaration d'invention -INPI -FR1913122), https://patentscope.wipo.int/search/fr/detail.jsf?docId=WO2021099509.

Step 1 is straightforward for the class of Hamiltonians such that (1 ≤ n ≤ N, N ≥ 1) (separated variables, H1) :

(C 1 -regularity, H2) :

(locally quadratic at 0, H3) :

(strict quasi-convexity, H4) :

Under hypotheses (H1-4), function

defines a C 1 -regular bijection that satisfies [START_REF] Hélie | Nonlinear damping models for linear conservative mechanical systems with preserved eigenspaces: a port-Hamiltonian formulation[END_REF] with L = I.

Proof:

• Function H n is continuous (H2), null at 0 (H3) and then, from (H4), non-negative strictly decreasing (resp. increasing) on R -(resp. R + ). It follows from ( 51) that q n is a real-valued strictly increasing continuous function.

• The C 1 -regularity of q n is obvious on R \ {0} (H2 and H(x) > 0), where

• Finally, Q is a collection of C 1 -regular bijections, that concludes the proof.

Remark 14 (non-separated variables). In the case of Hamiltonians with non-separated variables, sufficient conditions for the quadratisation can be found in [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF]Property 16,Example 19] and [27, Property 2.9].

Remark 15 (higher accuracy order). Following remark 10, due to the q-dependency of J q and R q the discrete gradient method is accurate only at order 1 (also with the quadratisation method). Higher-order methods that benefit from the quadratisation are available in [START_REF] Lopes | Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems[END_REF][START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF][START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios : approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF]: they involve multiple stages to refine the estimates and increments in (37a-38b), inspired from Runge-Kutta methods. These methods achieve properties (a-c) of remark 11 for orders 1, 2 or 3 without iterative solver (see [26, p.115-122]).

Note that after this method was introduced in [START_REF] Lopes | Explicit second-order accurate method for the passive guaranteed simulation of port-Hamiltonian systems[END_REF], an alternative method (also based on quadratisation but with auxiliary variables) was proposed in [START_REF] Yang | Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method[END_REF]. These methods are used in [START_REF] Danish | Applications of port hamiltonian methods to non-iterative stable simulations of the korg35 and moog 4-pole vcf[END_REF][START_REF] Ducceschi | Simulation of the geometrically exact nonlinear string via energy quadratisation[END_REF] for in audio-acoustic applications. Remark: observe that the conservative case (R 0 = 0) makes the representations derived in questions 1-2 degenerate in a singular way. Addressing this case requires to include constraints in the algebraic-differential formulations or to consider kernel representations of Dirac structures (see footnote 10 in § 2.6).

Exercise 2 (nonlinear circuit): quadratisation method

Consider the same circuit with a coil characterised by a saturating effort law of Hamiltonian H(φ ) = E 0 ln(cosh(φ /φ 0 )). Question 6. PHS. (a) How to modify the formulations (1) and ( 4) ? (b) What is the effort law of the coil ? (c) What is the equivalent inductance L 0 for small amplitudes ? Question 7. Quadratised PHS. Derive the quadratised PHS by applying the procedure of § 4.4 (steps 1-4). Question 8. Simulation and numerical experiments. Implement a code (by adapting that of exercise 1). Suggestion: to appreciate the nonlinear effect and simplify comparisons with exercise 1, choose φ 0 smaller than the maximal amplitude of your previous simulations, then fix E 0 to have the same L 0 as your previous code.


Philippe Lagacherie

Maider Arregui

David Fages

Malithi Weerasekara

Evaluating the quality of soil legacy data used as input of Digital Soil Mapping models

• Geo-referencing errors = differences of relative positions from remarkable points given by the geographical coordinates of soil profiles vs by the localisation plan • Measured on 43 randomly selected sites

Objective

Most of the Digital Soil Mapping products now available across the globe have been developed from legacy measured soil profiles inherited from several decades of soil survey activity. To our knowledge, no reference on the quality of these soil legacy data have been established yet. Our objectives was then i) to evaluate the differences between values of soil properties provided by legacy measured soil profiles and actual certified soil laboratory measurements , ii) to depict the different causes of these differences and iii) to investigate how to reduce some of them.

Philippe Lagacherie 1 , Maider Arregui • High relative Mean errors and Mean Square errors revealing great differences between legacy and actual soil property values

• Important contributions of analytical biases that could be corrected

• Variable and substantial contributions of the geo-referencing errors related with the relative importance of the local variability

• For some properties, a large amount of the difference remained undefined ( time evolution, analytical errors, recording errors,…)

Take home messages

• The legacy measured soil profiles usually used as input of our DSM models may have substantial biases and uncertainties that could be a major limitation of DSM prediction performances.

• Biases and a part of these uncertainties can however be corrected by using control sampling. Such control step should become systematic in DSM applications using legacy data

• Geopositionning is an important source of uncertainty for old soil datasets (< 2010) collected with manual geo-referencing

  of soil legacy data used as input of Digital Soil Mapping models

•

  2 , David Fages 1 & M. Weerasekara 1 Legacy measured soil profiles Control soil sampling • 8200 soil profiles with laboratory analysis collected between 1955 and 1963 by BRL soil surveyors • Integrated in the French Soil database using text recognition (Chemchem et al, this workshop) 137 sites over the BRL soil survey area (6 636 km 2 ) • Four samples per sites (local variations of soil properties) • Analysis of topsoil (0-20 cm) properties Differences between legacy and actual soil property values rME: relative Mean error : rME = mean(Zleg -Zref )/ Mean (Zref) rMSE: relative Mean Square Error = mean(Z 2 leg -Z 2 ref )/ variance (Zref) with Zleg: legacy soil property values, Zref: actual soil property values Correcting analytical biases Z corrected = a + b Z leg + c CaCO3 leg soil properties (% of total variance for 5m, 25m and 100m spacing) and impacts of a 31 m geo-referencing error on the differences of soil property values.


S E Sebastian

N Harrison

F Lévy

I Sheikin

C Berthier

M Horvatić

M Takigawa

Keywords: 75, 10, Jm -Quantized spin models 75, 30, Kz -Magnetic phase boundaries

In their comment, Sebastian and Harrison (SH) suggest that our torque and force measurements in SrCu 2 (BO 3 ) 2 [1] are incorrect due to a non-linear regime in the cantilever operation and deny the importance of Dzyaloshinski-Moriya (DM) interaction in torque measurements. In this reply, we show that their arguments are incorrect and why neither the torque measurements [1,2] nor the magnetization measurements in pulsed magnetic field [3,4] can give access to the field dependence of the longitudinal magnetization at thermodynamic equilibrium.

First of all, any non-linearity is perfectly excluded in all our measurements. Our largest signal reached only 10 % of the whole linear response range, and the variation of the capacitance has always been smaller than 1 % of the zero field capacitance. This means that the deflection of the cantilever was smaller that 10 -4 rad, which is comparable to the numbers given in the comment (0.5×10 -4 rad). Next, the angle θ between the c-axis of the sample and the applied magnetic field H in our experiment was smaller than 0.4 • , which is much smaller than the 2 • mentioned in [2] and comparable to the ≪ 1 • declared for the experiments reported in the comment.

SH base their interpretation of the torque results on the assumption that DM effects are negligible, as their data follow closely those obtained by inductive method in pulsed magnetic field, which couples only to the total magnetization M z . This statement is obviously wrong since in the gapped phase below 16 T the pulsed field measurements correctly record zero magnetization, while the torque measurements do record a non-zero signal -that is the signal due to DM interaction. Indeed, Ref. [2] explicitly mentions that at low fields the data have been corrected (put to zero) by hand, and in the Fig. 1 of the comment it is obvious that below 16 T there is an important variation of the torque signal. Furthermore, in the same figure the steep jump towards the 1/8 plateau Fig. 1: Magnetization from NMR and force measurements The field dependance of magnetization, and in particular the jump towards the 1/8 plateau, recorded by our force measurement [1] is in excellent agreement with values determined from the NMR data of Ref. [5] as explained in the text.

is obviously not the same in the torque and in the pulse field data. So the main assumption of SH is contradicted by their own data. Finally, Fig. 3 in Ref. [3] shows that below 30 T the magnetization measurements do not scale with the g-tensor anisotropy, at variance with their claims.

SH criticize the large jump in the magnetization curve observed in our force measurements [1], in which the deviation of the cantilever is dominated by the force F = dH dz M z . The variation of M z during this "jump" indeed corresponds to the coexistence of the uniform magnetization phase and the 1/8 plateau phase as already shown by copper NMR [5]. Kodama et al. have shown that up to p-1 Reply 26.5 T the system is in a uniform phase in which all copper electronic spins bear identical magnetization gµ B S z . In this phase, S z values have been accurately determined by NMR and found to increase up to the value of 0.034 at 26.5 T (Fig. 4A of Ref. [5]). Between 26.5 T and 27 T, the increase of the volume fraction x(H) of the plateau phase as a function of H has also been determined by NMR (Fig. 4C of Ref. [5]). One can thus easily deduce the field dependence of the average bulk magnetization S z = 0.034 (1x) + 0.0625 x, which is reported in Figure 1. Clearly, the agreement between the force measurement and the NMR data is excellent, meaning that the rapid variation of M z (H) in the field range 26.5 -27 T does correspond to the physics of the system. This clearly indicates that all previously reported measurements, whatever using torque technique in static field or flux integration in pulsed field, were unable to give the correct equilibrium M z (H) dependence.

In the following, we explain the possible origin of discrepancies between different techniques. In pulsed magnetic field measurements, during the fast increase of the field the energy levels of the lower triplet states are lowered towards the singlet level, but their populations have not enough time to relax to its thermal equilibrium value (adiabatic process). The effective temperature describing this non-equilibrium level populations is thus much lower than that of the experiment. This explains the observation of 1/8 magnetization plateau in pulsed magnetic field experiments performed at 1.4 K [3], 1.6 K and 0.6 K [4], although this phase is not stable above 0.55 K. This also implies that the values of M z in these experiments are incorrect.

Let us come to the torque measurements, and why the results of different groups can differ. There are three contributions to the torque applied on the cantilever. The first one is well known and related to the anisotropy of the diagonal part of the g-tensor. The second one is due to the field gradient (dH/dz = 0) if the sample is not rigorously at the center of the field. The third one is related to the presence of the DM interaction and the staggered g-tensor in the system. They induce a staggered transverse magnetization, but also a uniform transverse magnetization contributing to the torque [6]. This contribution vanishes at zero θ, but is otherwise always present in experimental results [1]. The interplay between the three contributions to the torque is very dependent on the position z of the sample with respect to the center of the field and the angle θ, which can be varied independently. We first focus on the first two terms. For small values of θ, the torque divided by the magnetic field can be written as: τ /H = [(g c -g a )/g c sin 2θ+ l/H dH/dz]M z , where l is the length of the cantilever, g c = 2.28 and g a = 2.05 [7]. Close to the center the field profile is parabolic, so that dH/dz = z d 2 H/dz 2 . Clearly one can find a position z 0 where the two terms cancel:

z 0 = -2θ/l (g c -g a )/g c [d 2 H/Hdz 2 ] -1 . For typical values of d 2 H
Hdz 2 = -50 ppm/mm 2 and l = 5 mm, this leads to z 0 (mm) = 14 θ( • ), e.g. z 0 = + 6 mm for θ = 0.4 • . That is, depending on the position of the sample, one can not simply predict the relative magnitude and sign of the contribution proportional to M z and that due to the transverse magnetization induced by the DM interaction. This can well explain why our torque measurements exhibit a negative slope within the plateau, while those of SH exhibit a positive slope. Unfortunately, there is no calculation yet of the transverse uniform magnetization generated by the DM and the staggered g tensor on a Shastry-Sutherland lattice, in presence of a superstructure of the magnetization corresponding to a plateau at M z /M saturation = 1/8 [5] or 1/9 [8].

One puzzling feature in our results, however, is that the signal amplitudes obtained in the force and in the torque measurements are comparable, although in principle one would expect the contribution of the DM interaction to be much smaller [6]. Since we have not done a systematic study as a function of the angle θ and the position z, we cannot disentangle the relative contributions of the longitudinal magnetization and the transverse one. Furthermore, it has been shown recently [9] that in frustrated spin ladders the DM interaction can give a huge contribution to the torque when the angle between the magnetic field and the DM vector is small. In SrCu 2 (BO 3 ) 2 the interdimer DM vector is nearly parallel to the c-axis [10], and may thus be a possible source of a strong torque signal. As regards our force measurement, considering their agreement with the NMR data, they appear to give a correct field dependence of the magnetization, at least up to the 1/8 plateau. The theoretical prediction for the absolute value of the magnetization in the plateaus in SrCu 2 (BO 3 ) 2 is currently controversial [2,8,11], which makes their experimental determination highly desirable. This requires a special experimental setup to minimize as much as possible any torque contribution with respect to the effect of force due to the field gradient applied on the sample [12].


Wichien Sriwichai

Myriam Collin

Timothy J Tranbarger

Jacques Berger

Sylvie Avallone

Improvement of the content in bioaccessible lipophilic micronutrients in raw and processed drumstick leaves (Moringa oleifera Lam.)

Keywords: cooking, vitamin retention, in vitro digestion, food destructuration, diffusion

Our objective was to evaluate the influence of cooking (boiling, steaming, frying) and processing (freezing, freeze-drying, sun drying, oven drying, sterilization or powdering) on the micronutrient profile of drumstick leaves (β-carotene, lutein, α-tocopherol) and the bioaccessibility of the latter after in vitro digestion. Fresh drumstick leaves contain a high level of carotenoids and α-tocopherol with a very low bioaccessibility (<1%). Freezing, steaming and sterilization preserve the more the micronutrients and most of the processing methods improve the contents in bioaccessible β-carotene, lutein and α-tocopherol. Among the cooking practices, only steaming was efficient to increase the compound release. The structural changes of the leaf tissue during processing partially explain the modulation of the bioactive compound release during the course of in vitro digestion. After boiling and steaming, the cells of epidermis and mesophyll were slightly disorganized while frying induced a deep destructuration of the overall tissue with cell wall rupture.

Introduction

The potential of indigenous trees to combat desertification and contribute to the well-being and food security of local population was studied since decades [START_REF] Nesbitt | Linking biodiversity, food and nutrition: The importance of plant identification and nomenclature[END_REF]. The drumstick tree (Moringa Oleifera Lam.) is interesting because it provides environmental services and could bring other added values for human well-being. The integration of Moringa tree in farming and pastoral systems could improve the sustainability of local food systems by providing livelihoods and diet diversification [START_REF] Morton | The horseradish tree, Moringa pterygosperma (Moringaceae). A boon to arid lands?[END_REF]. In traditional medicine, Drumstick is known for its effect on the digestive tract (bactericide, antidiarrheal, laxative), kidney function (diuretic) and anti-inflammatory effect (Arun [START_REF] Giridhari | Antidiabetic property of drumstick (Moringa oleifera) leaf tablets[END_REF][START_REF] Paliwal | Elucidation of free radical scavenging and antioxidant activity of aqueous and hydro-ethanolic extracts of Moringa oleifera pods[END_REF][START_REF] Sharma | Renoprotective effects of Moringa oleifera pods in 7,12-dimethylbenz [a] anthracene-exposed mice[END_REF][START_REF] Coppin | Determination of flavonoids by LC/MS and anti-inflammatory activity in Moringa oleifera[END_REF]. Several human studies have highlighted the anti-hyperglycemic, anti-dyslipidemic and antioxidant activities properties of drumstick [START_REF] William | Effect of some Indian vegetables on the glucose and insulin response in diabetic subjects[END_REF][START_REF] Kumari | Hypoglycemic effect of Moringa oleifera and Azadirachta indica in type-2 diabetes[END_REF][START_REF] Nambiar | Impact of antioxidants from drumstick leaves on the lipid profile of hyperlipidemics[END_REF][START_REF] Giridhari | Antidiabetic property of drumstick (Moringa oleifera) leaf tablets[END_REF][START_REF] Kushwaha | Effect of supplementation of drumstick (Moringa oleifera) and amaranth (Amaranthus tricolor) leaves powder on antioxidant profile and oxidative status among postmenopausal women[END_REF]. In northern countries, encapsulated drumstick powder is already commercialized by the pharmaceutical industry for its antioxidant capacity, anti-hyperglycemic as well as lipid and cholesterol lowering effect [START_REF] William | Effect of some Indian vegetables on the glucose and insulin response in diabetic subjects[END_REF][START_REF] Kumari | Hypoglycemic effect of Moringa oleifera and Azadirachta indica in type-2 diabetes[END_REF][START_REF] Nambiar | Impact of antioxidants from drumstick leaves on the lipid profile of hyperlipidemics[END_REF][START_REF] Giridhari | Antidiabetic property of drumstick (Moringa oleifera) leaf tablets[END_REF][START_REF] Kushwaha | Effect of supplementation of drumstick (Moringa oleifera) and amaranth (Amaranthus tricolor) leaves powder on antioxidant profile and oxidative status among postmenopausal women[END_REF]. In Africa and Asia, distinct parts of the tree are consumed such as flowers, leaves, fruits, bark, grains or roots. The leaves can be integrated in salad, traditional sauces, or soups in raw or processed form; young leaves can also be steam cooked [START_REF] Sallau | Effect of boiling, simmering and blanching on the antinutritional content of Moringa oleifera leaves[END_REF]. The drumstick leaves have a quantifiable amounts of proteins, essential fatty acids, minerals, vitamins (B, C, E, provitamin A) and a wide variety of polyphenols (flavonoids, chlorogenic acid) [START_REF] Ramachandran | Drumstick (Moringa oleifera): a multipurpose Indian vegetable[END_REF][START_REF] Makkar | Nutrional value and antinutritional components of whole and ethanol extracted Moringa oleifera leaves[END_REF][START_REF] Aslam | Mineral composition of Moringa oleifera leaves and pods from different regions of Punjab, Pakistan[END_REF][START_REF] Amaglo | Profiling selected phytochemicals and nutrients in different tissues of the multipurpose tree Moringa oleifera L., grown in Ghana[END_REF][START_REF] Coppin | Determination of flavonoids by LC/MS and anti-inflammatory activity in Moringa oleifera[END_REF]. The consumption of the leaves could fulfill the recommended daily requirements of several macro and micronutrients. In fresh or dried forms, the addition of drumstick leaves increases the content in iron, proteins and vitamin A of traditional dishes [START_REF] Ejoh | Small quantities of carotenoid-rich tropical green leafy vegetables indigenous to Africa maintain vitamin A status in Mongolian gerbils (Meriones unguiculatus)[END_REF]. In animal models (gerbil, mouse), the carotene bioavailability was high with raw and dehydrated drumstick leaf consumption [START_REF] Nambiar | Bioavailability trials of β-carotene from fresh and dehydrated drumstick leaves (Moringa oleifera) in a rat model[END_REF][START_REF] Ejoh | Small quantities of carotenoid-rich tropical green leafy vegetables indigenous to Africa maintain vitamin A status in Mongolian gerbils (Meriones unguiculatus)[END_REF]. After few weeks of depletion, the addition of tropical leafy vegetables to the gerbil's diet increased the vitamin A liver stores while serum retinol was not affected by the treatment. In the latter study, the high carotenoid bioavailability was attributed to the small particle size of powdered leaves [START_REF] Ejoh | Small quantities of carotenoid-rich tropical green leafy vegetables indigenous to Africa maintain vitamin A status in Mongolian gerbils (Meriones unguiculatus)[END_REF]. Vitamin deficiencies are public health issues in Africa and South-East Asia and still concern an estimated 190 million preschool-age children and 19.1 million pregnant women [START_REF] Who | Global prevalence of vitamin A deficiency in populations at risk 1995-2005[END_REF]. Populations with low incomes have insufficient micronutrient intakes especially in vitamin A. Some carotenoids with provitamin A activity (β-carotene) and vitamin E (α-tocopherol) have antioxidant properties and contribute to several functions such as cell growth and integrity, tissue differentiation, immunological defense and vision. The drumstick leaves are rich in these bioactive compounds but their biological effects are still controversial. Leafy vegetables could have a low nutritional efficacy because of their structure and biochemical composition (rich in fibers and phenolic compounds). The effect of processing on the nutritional value of drumstick leaves has been investigated for several years but the capacity of the leaves to deliver their micronutrient during digestion is not documented. The aim of the present work was to assess the influence of several treatments on the content in total and bioaccessible β-carotene, lutein and α-tocopherol of drumstick leaves.

Materials and methods

Materials and chemicals

Solvents, HPLC standards (β-carotene, α-carotene, lutein, α-tocopherol, δ-tocopherol) and enzymes (pepsin 800-2,500 U/mg, pancreatin B1750, bile extract porcin B8631) were obtained from Sigma-Aldrich (Saint Quentin Fallavier, France). The polytetrafluoroethylene (PTFE) membranes were obtained from Sartorius (Palaiseau, France). Capsule in hydroxypropyl methylcellulose (HPMC) were bought in a local drugstore. Drumstick leaves (1 kg) were purchased from Gopal and Co supermarket (Paris) and leaves with same color and size were chosen to limit experimental variations. Vegetables were cleaned with distilled water and the stem and other inedible parts were removed. Raw leaves (5 × 1 g) were directly submitted to in vitro digestion to assess the content in bioaccessible micronutrients. Then, treatments were applied to 11 batches of leaves (10 g) for further analysis.

Dry matter assessment

The dry matter (DM) contents of the raw and processed leaves were determined by drying the leaves in oven (24 h at 105 °C). Dry matter was then weighed precisely on a Mettler AE166 balance (Viroflay, France).

Processing methods and culinary treatments

Ten treatments were applied to samples of fresh drumstick leaves (2 g) with 5 repetitions for each one (Fig. 1). Some treatments correspond to Thai cooking and the others are preservation methods. Boiling was realized in milli-Q water, stir-frying in colza oil, and the equipment used for the treatments were as follows:

• freeze-drying Alpha 1-2 LD plus, Christ (Osterode, Germany) • oven dryer (Radnor, Pennsylvania, USA) • grinder from IKA (Staufen, Germany). Grinding for 1 and 3 min provides particle size of 1 mm and 200 µm respectively For encapsulation, 0.5 g of fine powdered leaves was introduced in hydroxypropyl methylcellulose (HPMC) capsules. After treatments, the samples were weighed to determine the cooking yields.

In vitro digestion protocol (IVD)

The IVD protocol was adapted from a published procedure with the addition of 80 µL of colza oil (diluted in acetone 32%) at the beginning of digestion to enhance carotenoid micellarisation [START_REF] Hedrén | In Vitro accessibility of carotenes from green leafy vegetables cooked with sunflower oil or red palm oil[END_REF][START_REF] Reboul | Bioaccessibility of carotenoids and vitamin E from their main dietary sources[END_REF]. Samples (2 g) were mixed in an amber screw-capped bottle with 15 mL of oral solution (pyrogallol 1.3% in NaCl 0.9%). Samples were incubated for 10 min at 37 °C with gentle stirring (IKAMAG® RT 15 power, IKA, Staufen, Germany). The gastric phase started then by adding 500 µL of 1 M HCl and adjusting the pH to 4.0 ± 0.1, and then 2 mL of porcine pepsin (4% in 0.1 M HCL). The mixture was incubated at 37 °C for 30 min with gentle stirring. Thereafter, the pH was raised to 6.0 ± 0.1 by adding 800 µL of 0.9 M of sodium bicarbonate buffer before the addition of 9 mL of the intestinal solution: porcine pancreatin (0.2%) and porcine bile extract (1.2%) dissolved in 0.1 M of trisodium citrate, pH 6. Then, 4 mL of porcine bile extract at 0.1 g/mL were added. The samples were incubated for 30 min in a stirring water bath at 37 °C. Digesta were centrifuged at 4000 g for 1 h at 10 °C (Heraeus Multifuge X1R, Thermo Fisher Scientific, Villebon sur Yvette, France). Finally, 10 mL was collected with a needle fitted to a syringe and filtered through 0.2 µm cellulose acetate membrane to recover the micellar fraction. Lipophilic compounds were then extracted.

Determination of total and bioaccessible carotenoids and tocopherols

Carotenoids and tocopherols were extracted from 50 mg of raw or processed leaves. Samples were ground (Turrax, Prolabo, Coueron, France) for 1 min with 15 mL of ethanol/hexane (4:3, v/v) and then centrifuged for 30 min at 10 °C and 4000 g (Heraeus Multifuge, Thermo Fisher Scientific, Villebon sur Yvette, France) [START_REF] Taungbodhitham | Evaluation of extraction method for the analysis of carotenoids in fruits and vegetables[END_REF]. The supernatant was transferred into a glass funnel and 10 mL of NaCl 10% was added (p/v). After 10 min, the lower phase was discarded and the organic phase was washed with 10 mL of distilled water. Organic extracts were recovered, dried under nitrogen and residues were dissolved in 2 mL of acetone. After filtration on 0.2 µm PTFE membranes (minisart SRP4, Sartorius), sample (20 µL) was injected into the HPLC. After in vitro digestion, the micronutrients were extracted on 10 mL of digesta by adding 20 mL of ethanol/hexane (4:3, v/v). Samples were homogenized and 20 mL of NaCl 25% were added (p/v). After 10 min, the organic phase was re-extracted with 20 mL of ethanol/hexane (4:3, v/v) and washed with 15 mL of distilled water. The organic phase was dried under nitrogen and residues were dissolved in 400 µL of acetone. The extracts were filtered through 0.2 µm PTFE membranes and analyzed with HPLC. Elution was performed with a Dionex P680 pump. The column was a polymeric C30 (4.6 mm i.d × 250 mm, 5 μm particle size, YMC Europe, Schermbeck, Germany). The mobile phase comprised of two mixes at a flow rate of 1 mL min-1. Mix A constituted of methanol and milli-Q water, 60:40 (v/v), and B mix constituted of methanol, methyl tert-butyl-ether and milli-Q water, 28.5:67.5:4 (v/v/v). The gradient applied progressed from 100% to 0% (A/B) over a period of 31 min and the temperature was approximately 25 °C. A UV-visible photodiode array detector (Dionex UVD 340U) was used and chromatograms were analyzed at the wavelength of maximum absorption of the carotenoids and α-tocopherol in the mobile phase (450 and 292 nm respectively). External calibration was used in the range of 0.5 to 25 mg/L. Retinol equivalents (RE) were calculated using a conversion factor of 12 for β-carotene and of 24 for the other provitamins A (α-carotene) for 1 retinol and expressed in mg RE/100 g [START_REF] Murphy | Dietary reference intakes for the U.S. and Canada: update on implications for nutrient databases[END_REF]. The content in bioaccessible β-carotene and α-tocopherol corresponds to the amounts of compounds transferred to the aqueous micellar phase at the end of IVD. The bioaccessibility (in %) corresponds to the amount of compounds transferred to the micellar phase compared to the initial contents in the fresh leaves.

Light microscopy

Samples (raw or processed leaves) were fixed for 2 weeks in 0.2 M phosphate buffer at pH 7.2, supplemented with 2% (v/v) paraformaldehyde, 1% (w/v) caffeine, and 1.5% (v/v) glutaraldehyde at 4°C (Buffard-Morel et al., 1992). Progressive dehydration of samples was realized with ethanol from 50 to 100%, transfer to butanol/ethanol 1:1 (v/v), and then 100% butanol. Each sample was impregnated in Technovit 7100 resin (Heraeus Kulzer GmbH, Wehrhrim/Ts, Germany). Semi-thin sections of 4 µm were cut using a Historange microtome (LKB Broma, Uppsala, Sweden). Sections were stained with specific reagents as follow: i) periodic acid-Schiff's reagent (Sigma-Aldrich, Lyon, France) and with naphtol blue black (NBB) to specifically reveal polysaccharides (in pink) and proteins (in dark blue) (Buffard-Morel et al., 1992); ii) ruthenium red to stain pectins (in pink) [START_REF] Sterling | Crystal-Structure of Ruthenium Red and Stereochemistry of its Pectic Stain[END_REF]; iii) toluidine blue to reveal tannins in greenish blue (Brillouet et al., 2012). Observations and microphotographs were performed with a Leica camera (Wetzlar, Germany) DMRB 301-370 light microscope (Evolution MP 5.0 color media Cybernetics Cooled RTV 32-0040C-140).

Statistical analysis

The micronutrient content of the moringa leaves was compared after 3 culinary treatments and 7 preservation methods applied with 5 repetitions. Data were assessed by F-test in one-way analysis of variance (ANOVA) using Statgraphics plus 5.1 (Virginia, USA). Significance was accepted at probability P < 0.05.

3.

Results and discussion

Effects of cooking and processing methods on the dry matter content and micronutrient profile

In the raw leaves, the major lipophilic compounds were β-carotene, lutein and ɑ-tocopherol and their respective quantities were 63, 32, and 57 mg for 100 g DM (Table 1). These results are in accordance with previously published data [START_REF] Subadra | Retention and storage stability of beta-carotene in dehydrated drumstick leaves (Moringa oleifera)[END_REF][START_REF] Nambiar | A study on β-carotene content of some green leafy vegetables of western India by high performance liquid chromatography[END_REF][START_REF] Saini | Effect of dehydration methods on retention of carotenoids, tocopherols, ascorbic acid and antioxidant activity in Moringa oleifera leaves and preparation of a RTE product[END_REF]. The DM was rather high 24 % as well as the RE content (5.3 mg for 100 g DM) as previously observed (2.5 mg for 100 g DM) [START_REF] Kidmose | Content of carotenoids in commonly consumed Asian vegetables and stability and extractability during frying[END_REF].

The DM in cooked leaves ranged from 8.5-22.2 % and drying leaded to dehydrated products with low residual water contents (4 to 6 %). The steam cooking, freezing and sterilization were the treatments which preserve the best β-carotene and retinol equivalent (RE), with respective loss ranging from 19-27 % and 16-37 %. In contrast, drastic losses of β-carotene and RE were observed after frying and sun drying.

Thermal and mechanical treatments (oven drying-fine grinding and encapsulation) also induced compound loss (>45%). During the leaf grinding and encapsulation, the surface area in contact with atmosphere and light was increased which caused the micronutrients loss by oxidation reactions. As observed in our work, sun drying (with light exposures) can induce drastic loss reaching 80 % [START_REF] Park | Effect of freezing, thawing, drying, and cooking on carotene retention in carrots, broccoli and spinach[END_REF]Craft, Wise, & Soares, 2002;[START_REF] Saini | Effect of dehydration methods on retention of carotenoids, tocopherols, ascorbic acid and antioxidant activity in Moringa oleifera leaves and preparation of a RTE product[END_REF][START_REF] Knecht | Tocopherol and tocotrienol analysis in raw and cooked vegetables: a validated method with emphasis on sample preparation[END_REF].

Steaming preserved more carotenoids than cooking in water where losses ranged from 24 to 83% [START_REF] Jayarajan | Effect of dietary fat on absorption of beta carotene from green leafy vegetables in children[END_REF][START_REF] Dikshit | Separation of carotenoids and estimation of beta-carotene content of selected Indian food and food preparations by HPLC[END_REF][START_REF] Grimme | Function of chlorophylls and carotenoids in thylakoid membranes: chlorophylls between pigment-protein complexes might function by stabilizing the membrane structure[END_REF][START_REF] Dietz | Reversed phase HPLC analysis of ɑ and β-carotene from selected raw and cooked vegetables[END_REF][START_REF] Mazzeo | Effect of two cooking procedures on phytochemical compounds, total antioxidant capacity and colour of selected frozen vegetables[END_REF]. Frying was the treatment that degraded the more the biochemical compounds.

Influence of cooking and processing on the contents in bioaccessible micronutrients

Once the IVD was completed, the lipophilic bioaccessible compounds are released from the food matrix and transferred to the micelles of the digestive medium. In raw leaves, the contents in bioaccessible β-carotene, lutein and ɑ-tocopherol were low (< 2 mg/100 g DM) as well as their relative bioaccessibility (< 3 %) (Fig. 2, 3 and 4).

Xanthophylls were more bioaccessible than carotenes as observed in numerous studies [START_REF] Pérez-Gálvez | Esterification of xanthophylls and its effect on chemical behavior and bioavailability of carotenoids in the human[END_REF][START_REF] Reboul | Bioaccessibility of carotenoids and vitamin E from their main dietary sources[END_REF][START_REF] Granado-Lorencio | In vitro bioaccessibility of carotenoids and tocopherols from fruits and vegetables[END_REF][START_REF] Palmero | Role of carotenoid type on the effect of thermal processing on bioaccessibility[END_REF]. The content in bioaccessible ɑ-tocopherol was always higher than carotenoids whatever the treatment. Most of the preservation methods (freeze drying, sterilization, sun drying, oven drying and encapsulation) increased significantly the contents in bioaccessible micronutrients and the effect of the preservation step was always higher for lutein (12-fold increase) than β-carotene (4-fold increase).

When drying was combined with grinding, the contents in bioaccessible micronutrients were improved. Grinding improved the diffusion of the compounds from the cells to the digestive medium; it opened the cells by breaking their walls and reducing the particle size to less than 1 mm or 200 µm according to the grinding time. It was previously demonstrated in wheat, that fine grinding allows disruption of plant cell walls and improves phenolic acid bioaccessibility (Hemery et al., 2010).

After steaming and sterilization, the leaves became softer than the raw ones (data not shown) probably due to the long duration of the thermal treatment (15 min). Improvement of carotenoid and ɑ-tocopherol micellarisation was already observed after boiling in water [START_REF] Nambiar | Bioavailability trials of β-carotene from fresh and dehydrated drumstick leaves (Moringa oleifera) in a rat model[END_REF][START_REF] Lemmens | Particle size reduction leading to cell wall rupture is more important for the β-Carotene bioaccessibility of raw compared to thermally processed carrots[END_REF][START_REF] Moelants | Relation between particle size and carotenoid bioaccessibility in carrot and tomato derived suspensions[END_REF]. Frying did not improve the content in bioaccessible compounds probably because the treatment induces important losses of carotenoids by chemical reactions [START_REF] Achir | Kinetic study of β-carotene and lutein degradation in oils during heat treatment[END_REF]. Encapsulation of dried and powdered leaves was an efficient way to deliver bioaccessible compounds. The HPMC capsule was not a physico-chemical barrier for the micronutrient release; experimentally, it was dissolved in 5 min in the digestive medium.

Structural changes after culinary treatments

To improve the micronutrient bioaccessibility, the cells have to be opened by destructuration.

The light micrographs clearly show that the leaf was organized with an upper and lower epidermis, and a mesophyll tissue comprised of palisade and spongy cells. The vascular bundles were also visible at the middle of the cross section (Fig. 5A-B). The palisade mesophyll was composed of elongated cells perpendicular to the upper epidermis; the contents of palisade cells stain blue with periodic acid-Schiff's reagent/NBB suggesting a high amount of protein (Fig. 5A). Tannin containing cells appeared greenish blue with toluidine blue, green with ruthenium red, but were also discernible as distinct cells with periodic acid-Schiff's reagent/NBB within the spongy mesophyll (Fig. 5A andG). According to the scale, we clearly see that the cells have a size (length and width) lower than 100 µm. Grinding the leaves to 1 mm or 200 µm particle only opened a fraction of the total leaf cells. Morphological and qualitative changes were observed after all the treatments. In boiled and steamed leaves, epidermal and palisade mesophyll cells appeared intact and were similar to untreated leaves (Fig. 5C, D, F, H and I). Some spongy cells appeared slightly disorganized after boiling and steaming. Frying promoted a loss of integrity of both epidermal and mesophyll cells, and a loss of protein and tannin. Upper and lower epidermal layers were no longer visible while the palisade and spongy cells appeared ruptured and disorganized. In contrast, both boiling and steaming did not appear to rupture the cells, while the blue staining pattern showed the protein content was more diffuse (Fig. 5C, D and F). No apparent structural difference was observed between the boiled and steamed that could explain the higher bioaccessibility measured in the latter. However, a comparison of sections of boiled and steamed samples stained with Ruthenium Red and NBB indicated a loss of tannin from cells in the steamed samples, while the tannins were still observed in the boiled samples (Fig. 5C,D,F).

In raw leaves, cell wall and plastid substructure are a physicochemical barrier for bioaccessibility [START_REF] Vishnevetsky | Carotenoid sequestration in plants: the role of carotenoid-associated proteins[END_REF][START_REF] Palmero | Novel targeted approach to better understand how natural structural barriers govern carotenoid in vitro bioaccessibility in vegetable-based systems[END_REF]. In general, cooking such as boiling improves the carotene bioaccessibility by changing cell morphology and softness [START_REF] De Belie | A kinetic model for turgor loss in red cabbage cells during mild heat treatment[END_REF][START_REF] De Belie | A kinetic model for turgor loss in red cabbage cells during mild heat treatment[END_REF]Miglio, Chiavaro, Visconti, Fogliano, & Pellegrini, 2008;[START_REF] Palmero | Role of carotenoid type on the effect of thermal processing on bioaccessibility[END_REF]. This was not observed in our case maybe because drumstick leaves are tree leaves with a higher mechanical resistance; their cell wall thickness is probably higher than other soft leafy vegetables. Cooking (boiling, steaming) with moist heat (with water or steam) involves gentle heat transfer from a liquid or vapor to food; the temperatures used are rather moderate (<110°C) and the vitamins are rather well retained by the leaves. On the contrary, during frying, the heat transfer is more intense because the oil temperature was maintained at 180°C; in these conditions, the compounds are highly degraded by chemical reactions (isomerisation and cleavage) and the food matrix is destructured. As a consequence, the vitamin bioaccessibility in % is increased because some cells were opened by frying but the overall content in bioaccessible vitamins is very low because most of them were degraded by chemical reactions. The methods of preservation tested are based on water activity decrease (drying, lyophilization), microorganism destruction (sterilization) or conservation in negative temperatures (freezing, lyophilization). Except for sterilization, the temperatures applied during the operation units were moderate and the vitamins were well retained. The combination of thermal and mechanical treatments by grinding is the best way to increase bioaccessible carotenoid and tocopherol contents because it induces some cell ruptures. The sterilization at 121°C under pressure also favors the leaf destructuration; they are easily dispersed in small pieces during the in vitro digestion (data not shown).

Conclusion

As other leafy vegetables, drumstick leaves contain a high level of carotenoids and vitamin E (ɑ-tocopherol) but, the contents in bioaccessible micronutrients are very low in the fresh product. The operation units which preserve the more the micronutrient and RE contents were freezing, steaming and sterilization while frying and sun-drying cause drastic loss. Drying leaves at moderate temperature (60°C) and grinding them was the best way to increase the content in bioaccessible β-carotene, lutein and tocopherol. Encapsulation of dried-powdered leaves was also an efficient way to deliver bioaccessible compounds during in vitro digestion. Even with the best combination of operation units (drying plus grinding), a great proportion of the bioactive compounds are still not released from the moringa leaves; grinding in smaller particle sizes (< 200 µm) would allow to better functionalize drumstick leaf powder.

Figure 1 .

 1 Figure 1. The treatments applied to drumstick leaves before in vitro digestion.

Figure 5 .

 5 Figure 5. Light micrographs of drumstick leaves (cross section, x20): raw leaves (A -C and G) and cooked leaves (D and I: boiled, E and J: fried, F and H: steamed). Periodic acid-Schiff's and naphtol blue black staining (A, C, D, E and F), toluidine blue staining (B), ruthenium red staining (G, H, I, J). E: epidermis. PC: palisade cells; SC: spongy cells; VB: vascular bundle, CT: condensed tannins, NCT: none contain condensed tannins (orange arrow pointing). Scale bars = 100 µm.

Acknowledgements

This work was supported by Thai Ministry of Sciences and Technology and Franco-Thai scholarship program. The authors report no conflicts of interest. The authors are responsible for the content and writing of the paper. The authors thank Marie-Jeanne Billet for her technical assistance. *Corresponding author. Tel.: 33(0)467 41 61 91; fax: 33(0)467 41 63 30. E-mail address: sylvie.avallone@supagro.fr