The paper presents a new stiffness modelling method for overconstrained parallel manipulators, which is applied to 3-d.o.f. translational mechanisms. It is based on a multidimensional lumped-parameter model that replaces the link flexibility by localized 6-d.o.f. virtual springs. In contrast to other works, the method includes a FEA-based link stiffness evaluation and employs a new solution strategy of the kinetostatic equations, which allows computing the stiffness matrix for the overconstrained architectures and for the singular manipulator postures. The advantages of the developed technique are confirmed by application examples, which deal with comparative stiffness analysis of two translational parallel manipulators.

I. INTRODUCTION / RELATED WORKS

elative to serial manipulators, parallel manipulators are claimed to offer an improved stiffness-to-mass ratio and better accuracy. This feature makes them attractive for innovative machine-tool structures for high speed machining [START_REF] Tlusty | Fundamental Comparison of the Use of Serial and Parallel Kinematics for Machine Tools[END_REF][START_REF] Wenger | A Comparative Study of Serial and Parallel[END_REF][START_REF] Majou | The design of Parallel Kinematic Machine Tools using Kinetostatic Performance Criteria[END_REF]. When a parallel manipulator is used as a Parallel Kinematic Machine (PKM), stiffness becomes a very important issue in its design [START_REF] Pritschow | Systematic Design of Hexapods and Other Parallel Link Systems[END_REF][START_REF] Company | Modelling and Design Issues of a 3-axis Parallel Machine-Tool[END_REF][START_REF] Brogardh | PKM Research -Important Issues, as seen from a Product Development Perspective at ABB Robotics[END_REF][START_REF] Paskhevich | Kinematic and stiffness analysis of the Orthoglide, a PKM with simple, regular workspace and homogeneous performances[END_REF]. This paper presents a general method to compute the stiffness analysis of 3-dof overconstrained translational parallel manipulators.

Generally, the stiffness analysis of parallel manipulators is based on a kinetostatic modeling [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF], which proposes a map of the stiffness by taking into account the compliance of the actuated joints. However, this method is not appropriate for PKM whose legs are subject to bending [START_REF] Kong | Kinematics and Singularity Analysis of a Novel Type of 3-CRR 3-DOF Translational Parallel Manipulator[END_REF].

Several methods exist for the computation of the stiffness matrix: the Finite Element Analysis (FEA) [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF], the matrix structural analysis (SMA) [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot uncoupled kinematics[END_REF], and the virtual joint method (VJM) that is often called the lumped modeling [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF].

The FEA is proved to be the most accurate and reliable, however it is usually applied at the final design stage because of the high computational expenses required for the repeated remeshing of the complicated 3D structure over the whole workspace. The SMA also incorporates the main ideas of the FEA, but operates with rather large elements -3D flexible beams describing the manipulator structure. This leads obviously to the reduction of the computational expenses, but does not provide clear physical relations required for the parametric stiffness analysis. And finally, the VJM method is based on the expansion of the A. Pashkevich is with the IRCCyN (UMR CNRS 6597), Nantes, France and with the Department of Automatics and Production Systems, École des Mines de Nantes, France (anatol.pashkevich@emn.fr); D. Chablat is with the IRCCyN (UMR CNRS 6597), Nantes, France (Damien.Chablat@irccyn-nantes.fr); P. Wenger is with the IRCCyN (UMR CNRS 6597), Nantes, France (Philippe.Wenger@irccyn-nantes.fr).

traditional rigid model by adding the virtual joints (localized springs), which describe the elastic deformations of the links. The VJM technique is widely used at the pre-design stage.

Next section introduces a general methodology to derive the kinematic and stiffness model. Section 3 describes the manipulator compliant elements and the link stiffness evaluation methods. Finally in section 4, we apply our method on two application examples.

II. GENERAL METHODOLOGY

A. Manipulator Architecture

Let us consider a general 3 d.o.f. translational parallel manipulator, which consists of a mobile platform connected to a fixed base by three identical kinematics chains (Fig. 1). Each chain includes an actuated joint "Ac" (prismatic or rotational) followed by a "Foot" and a "Leg" with a number of passive joints "Ps" inside. Certain geometrical conditions are assumed to be satisfied with respect to the passive joints to eliminate the platform rotations and to achieve stability of its translational motions. ; where each leg consists of a rod ended by two U-joints (with parallel intermediate and exterior axes), and active joint is driven by linear actuator [START_REF] Li | Stiffness Analysis for a 3-PUU Parallel Kinematic Machine[END_REF];

(b) Delta parallel robot (Fig 2b) that is based on the 3-RRPaR architecture with parallelogram-type legs and rotational active joints [START_REF] Clavel | DELTA, a fast robot with parallel geometry[END_REF];

(c) Orthoglide parallel robot (Fig 2c) that implements the 3-PRPaR architecture with parallelogram-type legs and translational active joints [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF].

Here R, P, U and Pa denote the revolute, prismatic, universal and parallelogram joints, respectively.

It should be noted that examples (b) and (c) illustrate overconstrained mechanisms, where some kinematic constrains are redundant but do not affect the resulting degrees of freedom. However, most of the past works deal with non-overconstrained architectures, which motivates the subject of this paper [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF].

B. Basic Assumptions

To evaluate the manipulator stiffness, let us apply a modification of the virtual joint method (VJM), which is based on the lump modeling approach [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF][START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF]. According to this approach,

Stiffness Analysis of 3-d.o.f. Overconstrained Translational Parallel Manipulators

Anatoly Pashkevich, Damien Chablat, Philippe Wenger R the original rigid model should be extended by adding the virtual joints (localized springs), which describe elastic deformations of the links. Besides, virtual springs are included in the actuating joints to take into account stiffness of the control loop. To overcome difficulties with parallelogram modeling, let us first replace the manipulator legs (see Fig. 3) by rigid links with configuration-dependent stiffness. This transforms the general architecture into the extended 3-xUU case allowing treating all the considered manipulators in the similar manner. Under such assumptions, each kinematic chain of the manipulator can be described by a serial structure (Fig. 3), which includes sequentially: ( )

i i a q θ + V
where 0 i q is the actuated coordinate and 0 i θ is the virtual spring coordinate; (c) a rigid "Foot" linking the actuating joint and the leg, which is described by the constant homogenous transformation matrix foot T ;

(d) a 6-d.o.f. virtual joint defining three translational and three rotational foot-springs, which are described by the homogenous matrix function { , } i i q q , which is described by the homogenous matrix function 1 1 2 ( , )

i i u q q V ; (f) a rigid "Leg" linking the foot to the movable platform, which is described by the constant homogenous matrix transformation leg T ;

(g) a 6-d.o.f. virtual joint defining three translational and three rotational leg-springs, which are described by the homogenous matrix function ( , ) (h) a 2-d.o.f. passive U-joint at the end of the leg allowing two independent rotations with angles 3 4 { , } i i q q , which is described by the homogenous matrix function 2 3 4 ( , )

i i s θ θ V K ,
i i u q q V ; (i) a rigid link from the manipulator leg the end-effector (part of the movable platform) described by the constant homogenous matrix transformation i tool T ;

The expression defining the end-effector location subject to variations of all coordinates of a single kinematic chain may be written as follows

0 0 1 6 1 1 2 7 12 2 3 4 ( ) ( , ) ( , ) ( , ) ( , )
i i i i i i base a foot s i i i i i i i u leg s u tool q q q q q θ θ θ θ θ = ⋅ + ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ T T V T V V T V V T K K (1)
where matrix function (.) a V is either an elementary rotation or translation, matrix functions 1 (.) u V and 2 (.) u V are compositions of two successive rotations, and the spring matrix (.) s V is composed of six elementary transformations. In the rigid case, the virtual joint coordinates 0

1 2 , i i θ θ
K are equal to zero, while the remaining ones (both active 0 i q and passive 0 4 , i i q q K ) are obtained through the inverse kinematics, ensuring that all three matrices , 1, 2,3 i i = T are equal to the prescribed one that characterizes the spatial location of the moving platform (kinematic loop-closure equations). Particular expressions for all components of the product (1) may be easily derived using standard techniques for the homogenous transformation matrices. It should be noted that the kinematic model ( 1) includes 18 variables (1 for active joint, 4 for passive joints, and 13 for virtual springs). However, some of the virtual springs are redundant, since they are compensated by corresponding passive joints with aligning axes or by combination of passive joints. For computational convenience, nevertheless, it is not reasonable to detect and analytically eliminate redundant variables at this step, because the developed below technique allows easy and efficient computational elimination.

C. Differential Kinematic Model

To evaluate the manipulator ability to respond to the external forces and torques, let us first derive the differential equation describing relations between the end-effector location and small variations of the joint variables. For each ith kinematic chain, this equation can be generalized as follows , 1 ,2,3

i i i i q i i θ δ = ⋅δ + ⋅δ = t J θ J q , ( 2
)
where the vector δ

(δ , δ , δ , δ , δ , δ ) T i x i y i z i x i y i z i p p p ϕ ϕ ϕ = t describes the translation δ (δ , δ , δ ) T i x i y i z i p p p = p and the rotation δ (δ , δ , δ ) T i x i y i z i ϕ ϕ ϕ = ϕ
of the end-effector with respect to the Cartesian axes; vector

0 1 2 ( , ) i i T i θ θ δ = δ δ θ K collects all virtual joint coordinates, vector 1 4
( , )

i i T i q q δ = δ δ q K includes all passive
joint coordinates, symbol ' ' δ stands for the variation with respect to the rigid case values, and ι θ J , q ι J are the matrices of sizes 6×13 and 6×4 respectively. It should be noted that the derivative for the actuated coordinate 0 i q is not included in q ι J but it is represented in the first column of ι θ J through variable 0 i θ . The desired matrices ι θ J , q ι J , which are the only parameters of the differential model (2), may be computed from (1) analytically, using some software support tools, such as Maple, MathCAD or Mathematica. However, a straightforward differentiation usually yields very awkward expressions that are not convenient for further (a) 3-PUU translational PKM [START_REF] Li | Stiffness Analysis for a 3-PUU Parallel Kinematic Machine[END_REF] (b) Delta parallel robot [START_REF] Clavel | DELTA, a fast robot with parallel geometry[END_REF] (c) Orthoglide parallel robot [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF] Fig. 2. Typical 3 d.o.f. translational parallel mechanisms computations. On the other hand, the fractionized structure of [START_REF] Tlusty | Fundamental Comparison of the Use of Serial and Parallel Kinematics for Machine Tools[END_REF], where all variables are separated, allows applying an efficient semi-analytical method. To present this technique, let us assume that for the particular virtual joint variable 0 i θ the model ( 1) is rewritten as

1 2
( )

i i i j j j i j θ = ⋅ θ ⋅ T H V H , (3)
where the first and the third multipliers are the constant homogenous matrices, and the second multiplier is the elementary translation or rotation. Then the partial derivative of the homogenous matrix i T for the variable i j θ at point 0 i j θ = may be computed from a similar product where the internal term is replaced by (.)

j θ ′ V
that admits very simple analytical presentation. In particular, for the elementary translations and rotations about the X-axis, these derivatives are:

0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 x Tran ⎡ ⎤ ⎢ ⎥ ′ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ V ; 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 x Rot ⎡ ⎤ ⎢ ⎥ = - ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ V . ( 4
)
Furthermore, since the derivative of the homogenous matrix 1 2

( )

i i i j j j i j θ ′ ′ = ⋅ θ ⋅ T H V H may be presented as 0 0 0 0 0 0 0 iz iy ix iz ix iy i iy ix iz p p p ϕ ϕ ϕ ϕ ϕ ϕ ′ ′ ′ - ⎡ ⎤ ′ ′ ′ ⎢ ⎥ - ′ = ⎢ ⎥ ′ ′ ′ - ⎢ ⎥ ⎣ ⎦ T , (5)
then the desired jth column of ι θ J can be extracted from i ′ T (using the matrix elements 14

T ′ , 24 T ′ , 34 T ′ , 23 T ′ , 31 T ′ , 12
T ′ ). The Jacobians q ι J can be computed in a similar manner, but the derivatives are evaluated in the neighborhood of the "nominal" values of the passive joint coordinates i j nom q corresponding to the rigid case (these values are provided by the inverse kinematics). However, simple transformation i i i j j j nom q q q δ = + and corresponding factoring of the function ( ) ( ) ( ) i i i q j j q j j q j j nom q q q δ = V V V allow applying the above approach. It is also worth mentioning that this technique may be used in analytical computations, allowing one to avoid bulky transformations produced by the straightforward differentiating.

D. Kinetostatic and Stiffness Models

For the manipulator kinetostatic model, which describes the force-and-motion relation, it is necessary to introduce additional equations that define the virtual joint reactions to the corresponding spring deformations. In accordance with the adopted stiffness model, three types of virtual springs are included in each kinematic chain:

• 1-d.o.f. virtual spring describing the actuator compliance; • 6-d.o.f. virtual spring describing compliance of the foot; • 6-d.o.f. virtual spring describing compliance of the leg.

Assuming that the spring deformations are small enough, the required relations may be expressed by linear equations

0 i i act K τ θ 0 = θ ⎡ ⎤ ⎡ ⎤ ⎣ ⎦ ⎣ ⎦ ; 1 6 i i Foot i i τ τ θ1 θ 6 ⎡ ⎤ θ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ θ ⎣ ⎦ ⎣ ⎦ K M M ; 7 12 i i Leg i i τ τ θ 7 θ12 ⎡ ⎤ θ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ θ ⎣ ⎦ ⎣ ⎦ K M M , (6)
where i j τ θ is the generalized force for the jth virtual joint of the ith kinematic chain, act K is the actuator stiffness (scalar), and, Foot K , Leg K are 6×6 stiffness matrices for the foot and leg respectively. It should be stressed that, in contrast to other works, these matrices are assumed to be non-diagonal. This allows taking into account complicated coupling between rotational and translational deformations, while usual lump-based approach does consider this phenomena [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF]. For analytical convenience, expressions (6) may be collected in a single matrix equation θ , 1 ,2,3

i i i θ = ⋅δ = τ K θ (7)
where 0 1 2 ( , )

i i i T τ τ θ θ θ = τ K
is the aggregated vector of the virtual joint reactions, and ( , , )

act Foot Leg diag K = θ K K K
is the aggregated spring stiffness matrix of the size 13×13. Similarly, one can define the aggregated vector of the passive joint reactions ( , )

i i i T q q q τ τ = τ
K but all its components must be equal to zero: , 1 ,2,3

i q i = = τ 0 (8)
To find the static equations corresponding to the end-effector motion i δt , let us apply the principle of virtual work assuming that the joints are given small, arbitrary virtual displacements ( , )

i i
Δ Δ θ q in the equilibrium neighborhood. Then the virtual work of the external force i f applied to the end-effector along the corresponding displacement

i i i i q i θ Δ = Δ + Δ t J θ J q is equal to the sum ( ) ( ) T i T i i i i q i θ Δ + Δ f J θ f J q .
For the internal forces, the virtual work is

θ T i i
-⋅ Δ τ θ since the passive joints do not produce the force/torque reactions (the minus sign takes into account the adopted directions for the virtual spring forces/torques). Therefore, because in the static equilibrium the total virtual work is equal to zero for any virtual displacement, the equilibrium conditions may be written as

T i i i θ θ ⋅ = J f τ ; T i q i ⋅ = J f 0 . (9)
This gives additional expressions describing the force/torque propagation from the joints to the end-effector. Hence, the complete kinetostatic model consists of five matrix equations ( 2), ( 7)…( 9) where either i f or i δt are treated as known, and the remaining variables are considered as unknowns. Obviously, since separate kinematic chains posses some degrees-of-freedom, this system cannot be uniquely solved for given i f . However, vice versa, for given end-effector displacement i δt , it is possible to compute both the corresponding external force i f and the internal variables, i δθ , i θ τ , i δq (i.e. virtual spring reactions and displacements in passive joints, which may also provide useful information for the designer). Since matrix θ K is non-singular (it describes the stiffness of the virtual sprigs), the variable i δθ can be expressed via i f using equations

θ i i θ = ⋅δ τ K θ and T i i i θ θ ⋅ = J f τ . This yields substitution 1 θ ( ) T i i i - θ δ = ⋅ θ K J
f allowing reducing the kinetostatic model to system of two matrix equations 1 θ ( )

T i i i i q i i - θ θ ⋅ + ⋅δ =δ J K J f J q t ; T i q i ⋅ = J f 0 (10)
with unknowns i f and i Δq . This system can be also rewritten in a matrix form

θ i i q i i T i i q ⎡ ⎤ δ ⎡ ⎤ ⎡ ⎤ ⋅ = ⎢ ⎥ δ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦ S J f t q 0 J 0 (11)
where the sub-matrix

1 θ θ T i i i - θ θ = S J K J
describes the spring compliance relative to the end-effector, and the sub-matrix i q J takes into account the passive joint influence on the end-effector motions. Therefore, for a separate kinematic chain, the desired stiffness matrix i K defining the motion-to-force mapping

i i i = ⋅ δ f K t ,
can be computed by direct inversion of relevant 10×10 matrix in the left-hand side of [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot uncoupled kinematics[END_REF] and extracting from it the 6×6 sub-matrix with indices corresponding to θ i S . It is also worth mentioning that computing θ i S requires 6×6 inversions only, since ,  ) [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot uncoupled kinematics[END_REF] in general case, i.e. for any given i θ J and i q J , cannot be proved. Moreover, if the matrix i q J is singular, the passive joint coordinates i q can not be found uniquely. From a physical point of view, it means that if the kinematic chain is located in a singular posture, then certain displacements i δt can be generated by infinite combinations of the passive joints. But for the variable i f the corresponding solution is unique (since the matrix i θ J is obviously non-singular if at least one 6 d.o.f. spring is included in a serial kinematic chain). On the other hand, the singularity may produce an infinite number of stiffness matrices for the same spatial location of the end-effector and for different values i q provided by the inverse kinematics. A special technique to tackle this case, based on the singular value decomposition, has been also developed. After the stiffness matrices i K for all kinematic chains are computed, the stiffness of the entire manipulator can be found by simple addition

act Foot Leg diag K - - - - θ = K K K . Solvability of system
3 1 i m i = = ∑ K K .
This follows from the superposition principle, because the total external force corresponding to the end-effector displacement δt (the same for all kinematic chains) can be expressed as

3 1 i i = = ∑ f f where i i = ⋅δ f K t .
It should be stressed that the resulting matrix i K is not invertible, since some motions of the end-effector do not produce the virtual spring reactions (because of passive joints influence). However, for the entire manipulator, the stiffness matrix m K is s positive definite and invertible for all non-singular (for the rigid model) postures.

E. Comparison with Other Results

The main advantage of the proposed methodology is its applicability to overconstrained mechanisms. To describe it in more details, let us briefly review an alternative technique [START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF]. The latter is originated from the same principal equations but the solution strategy includes straightforward elimination of the passive joint variables i q using the differential kinematic equations (2) only. Obviously, the feasibility of this step depends on the solvability of the equivalent matrix system

1 1 1 2 2 1 2 3 2 3 3 3 q q q θ θ θ δ ⎡ ⎤ - ⎡ ⎤ ⎡ ⎤ δ ⎡ ⎤ δ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ - ⋅ = ⋅ δ ⎢ ⎥ ⎢ ⎥ δ ⎢ ⎥ ⎢ ⎥ δ ⎢ ⎥ ⎣ ⎦ - ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ δ ⎣ ⎦ ⎣ ⎦ t I J J θ q I J J θ q θ J I J q (12)
where δt and i δq are treated as unknowns. In the non- constrained case (for the 3-PUU architecture, for instance) the matrix in the left-hand side of ( 14) is square, of size 18×18, so it can be inverted usually. However, for overconstrained manipulators, this matrix is non-square, so the system cannot be solved uniquely. For example, for manipulators with the parallelogram-type legs (Orthoglide, Delta, etc.) the matrix size is 18×15. So, in [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF] three additional (virtual) passive joints were introduced to solve the problem. But, obviously, such a modification changes the manipulator architecture and its stiffness matrix, doubting validity of the corresponding model. Besides, the developed technique allows computing the stiffness matrix even for the singular manipulator postures and does not incorporate the least-square pseudo-inversions applied by other authors. This is achieved by applying another solution strategy, which considers simultaneously the kinematic and static-equilibrium equations for each kinematic chain separately. Some hidden conveniences are included in the modeling stage. In particular, the kinematic models of the chains may include several redundant springs that are totally compensated by relevant passive joints. However, there is no need to eliminate these springs from the model manually, since they do not increase the matrix sizes in system [START_REF] Bouzgarrou | Rigidity analysis of T3R1 parallel robot uncoupled kinematics[END_REF]. This allows including in the model 6-d.o.f. virtual springs of general type, without any modifications. Another advantage of the proposed technique is that it can be generalized easily. Within this paper, it is applied to the stiffness modeling of 3-d.o.f. translational manipulators with actuators located between the base and the foot. However, it can be easily modified to cover other actuator locations, which may be included in the foot or in the leg. A further generalization is related to a number of kinematic chains and their similarity. They are also not crucial assumptions and influence on the Jacobian computing only. But after the Jacobians are determined, the stiffness matrices for separate chains may be computed in the same manner and then aggregated.

III. PARAMETERS OF THE COMPLIANT ELEMENTS

The adopted stiffness model of each kinematic chain includes three compliant components, which are described by one 1-d.o.f. spring and two 6-d.o.f. springs corresponding to the actuator, and to the foot/leg links (see Fig. 3). Let us describe particular techniques for their evaluation.

A. Actuator Compliance

The actuator compliance, described by the scalar parameter

1 act act k K - =
, depends on both the servomechanism mechanics and the control algorithms. Since most of modern actuators implement the digital PID control, the main contribution to act k is done by the mechanical transmissions. The latter are usually located outside the feedback-control loop and consist of screws, gears, shafts, belts, etc., whose flexibility is comparable with the flexibility of the manipulator links. Because of the complicated mechanical structure of the servomechanisms, the parameter act k is usually evaluated from static load experiments, by applying the linear regression to the experimental data.

B. Link Compliance

Following a general methodology, the compliance of a manipulator link (foots and legs) is described by 6×6 symmetrical positive definite matrices [START_REF] Tlusty | Fundamental Comparison of the Use of Serial and Parallel Kinematics for Machine Tools[END_REF] 1

, leg foot --K K corresponding to 6-d.o.f. springs with relevant coupling between translational and rotational deformations. This distinguishes our approach from other lumped modeling techniques, where the coupling is neglected and only a subset of deformations is taken into account (presented by a set of 1-d.o.f. springs). The simplest way to obtain these matrices is to approximate the link by a beam element for which the non-zero elements of the compliance matrix may be expressed analytically: Here L is the link length, A is its cross-section area, I y , I z , and J are the quadratic and polar moments of inertia of the cross-section, and E and G are the Young's and Coulomb's modules respectively. However, for certain link geometries, the accuracy of a singlebeam approximation can be insufficient. In this case the link can be approximated by a serial chain of the beams, whose compliance is evaluated by applying the same method (i.e. considering the kinematic chain with 6-d.o.f. virtual springs, but without passive joints). This leads to the resulting compliance matrix

C. FEA-based evaluation of stiffness

For complex link geometries, the most reliable results can be obtained from the FEA modeling. To apply this approach, the CAD model of each link should be extended by introducing an auxiliary 3D object, a "pseudo-rigid" body, which is used as a reference for the compliance evaluation. Besides, the link origin must be fixed relative to the global coordinate system. Then, sequentially and separately applying forces , ,

x y z F F F and torques , ,

x y z M M M to the reference object, it is possible to evaluate corresponding linear and angular displacements, which allow computing the stiffness matrix columns. The main difficulty here is to obtain accurate displacement values by using proper FEA-discretization ("mesh size"). Besides, to increase accuracy, the displacements must be evaluated using redundant data set describing the reference body motion. For this reason, it is worth applying a dedicated SVDbased algorithm. As follows from our study, the single-beam approximation of the Orthoglide foot gives accuracy of about 50%, and the four-beam approximation improves it up to 30% only. While the FEA-based method is the most accurate one, it is also the most time consuming. However, in contrast to the straightforward FEA-modeling of the entire manipulator, which requires recomputing for each manipulator posture, the proposed technique involves a single evaluation of link stiffness.

IV. APPLICATION EXAMPLES

To demonstrate efficiency of the proposed methodology, let us apply it to the comparative stiffness analysis of two 3-d.o.f. translational mechanism, which employ Orthoglide architecture. CAD models of these mechanisms are presented in Fig. 4.

θ θ = - ⋅ +Δ ⋅ ⋅ ⋅ ⋅ -+ Δ ⋅ T T R T V R T K (15) 1 1 6 2 ( /2) ( ) ( ) ( , ) ( ) ( /2) dn dn z y x dn dn dn s y z d q q L q q d θ θ = ⋅ +Δ ⋅ ⋅ ⋅ ⋅ -+ Δ ⋅ - T T R T V R T K (16)
where L, d are the parallelogram geometrical parameters, 1 2 , , { , } i i q q i up dn Δ Δ ∈ are the variations of the passive joint coordinates and the sub/superscripts "up" and "dn" correspond to the upper and lower chain respectively. Hence, the parallelogram compliance matrix may be also derived using the proposed technique that yields an analytical expression

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 4 8 0 0 0 0 0 4 0 0 0 8 4 q q Plg q q q K K K d C K d S K K d C K d S K d S K K K ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ = ⋅ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ + ⎢ ⎥ ⎣ ⎦ K (17)
where cos( ); sin( ) q q C q S q = = . Using this model and applying the proposed technique, we computed the compliance matrices for three typical manipulator postures (see table Table 1). As follows from the comparison with the U-joint case, the parallelograms allow increasing the rotational stiffness roughly in 10 times. This justifies application of this architecture in the Orthoglide prototype design [START_REF] Chablat | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF].

V. CONCLUSIONS

The paper proposes a new systematic method for computing the stiffness matrix of overconstrained parallel manipulators. It is based on multidimensional lumped model of the flexible links, whose parameters are evaluated via the FEA modeling and describe both the translational/rotational compliances and the coupling between them. In contrast to previous works, the method employs a new solution strategy of the kinetostatic equations, which considers simultaneously the kinematic and static relations for each separate kinematic chain and then aggregates the partial solutions in a total one. This allows computing the stiffness matrices for overconstrained mechanisms for any given manipulator posture, including singular configurations and their neighborhood. Another advantage is computational simplicity that requires low-dimensional matrix inversion compared to other techniques. Besides, the method does not require manual elimination of the redundant spring corresponding to the passive joints, since this operation is inherently included in the numerical algorithm. The efficiency of the proposed method was demonstrated through application examples, which deal with comparative stiffness analysis of two parallel manipulators of the Orthoglide family (with U-joint based and parallelogram based links). Relevant simulation results have confirmed essential advantages of the parallelogram based architecture and validated adopted design of the Orthoglide prototype. Another contribution is the analytical stiffness model of the parallelogram, which was derived using the same methodology. While applied to the 3-d.o.f. translational mechanisms, the method can be extended to other parallel architectures composed of several kinematic chains with rotational/prismatic joints and limb-or parallelogram-based links. So, future work will focus on the stiffness modeling of more complicated parallel mechanism with another actuator location (such as the Verne machine [START_REF] Kanaan | Kinematics analysis of the parallel module of the VERNE machine[END_REF]) and also on the experimental verification of the stiffness models for the Orthoglide robot.

x z y B 1 P A 1 C 1 A 2 C 2 A 3 C 3 B 3 B 1 i 1 P x z y j 1 A 1 C 1 A 2 B 2 C 2 A 3 C 3 B 3 B 1 i 1 P x z y j 1 A 1 C 1 A 2 B 2 C 2 A 3 C 3 B 3 (A) U-

Fig. 1 .

 1 Fig. 1. Schematic diagram of a general 3-d.o.f. translational parallel manipulator (Ac -actuated joint, Ps -passive joints, F -foot, L -Leg) Typical examples of such architectures are: (a) 3-PUU translational PKM (Fig 2a); where each leg consists of a rod ended by two U-joints (with parallel intermediate and exterior axes), and active joint is driven by linear actuator [13]; (b) Delta parallel robot (Fig 2b) that is based on the 3-RRPaR architecture with parallelogram-type legs and rotational active joints [14]; (c) Orthoglide parallel robot (Fig 2c) that implements the 3-PRPaR architecture with parallelogram-type legs and translational active joints [10].Here R, P, U and Pa denote the revolute, prismatic, universal and parallelogram joints, respectively.It should be noted that examples (b) and (c) illustrate overconstrained mechanisms, where some kinematic constrains are redundant but do not affect the resulting degrees of freedom. However, most of the past works deal with non-overconstrained architectures, which motivates the subject of this paper[START_REF] Gosselin | Stiffness mapping for parallel manipulators[END_REF].

Fig. 3 .

 3 Fig. 3. Flexible model of a single kinematic chain (a) a rigid link between the manipulator base and the ith actuating joint (part of the base platform) described by the constant homogenous transformation matrix i base T ; (b) a 1-d.o.f. actuating joint with supplementary virtual spring, which is described by the homogenous matrix function 0 0

K

  incorporate the Jacobian and the compliance matrices for all virtual springs.

  JOINT BASED ARCHITECTURE ( B) PARALLELOGRAM BASED ARCHITECTURE (C) WORKSPACE AND CRITICAL POINTS Q 1 AND Q 2 FIG. 4. KINEMATICS OF TWO 3-DOF TRANSLATIONAL MECHANISMS EMPLOYING THE ORTHOGLIDE ARCHITECTURE

TABLE I :

 I TRANSLATIONAL AND ROTATIONAL STIFFNESS OF THE 3-PUU AND 3-PRPAR MANIPULATORS

	MANIPULATOR ARCHITECTURE	Point Q 0 , , 0.00 x y z =	mm	Point Q 1 73.65 x y z , , = -	mm	, , x y z	Point Q 2 126.35 = +	mm
		tran k [N/mm]	rot k [N⋅mm/rad]	tran k [N/mm]	rot k [N⋅mm/rad] tran k [N/mm]	rot k [N⋅mm/rad]
	3-PUU manipulator	2.78⋅10 -4	20.9⋅10 -7	10.9⋅10 -4	24.1⋅10 -7	71.3⋅10 -4		25.8⋅10 -7
	3-PRPaR manipulator	2.78⋅10 -4	1.94⋅10 -7	9.86⋅10 -4	2.06⋅10 -7	21.2⋅10 -4		2.65⋅10 -7

A. Stiffness of U-Joint Based Manipulator

First, let us derive the stiffness model for the simplified Orthoglide mechanics, where the legs are comprised of equivalent limbs with U-joints at the ends. Accordingly, to retain major compliance properties, the limb geometry corresponds to the parallelogram bars with doubled cross-section area. Let us assume that the world coordinate system is located at the end-effector reference point corresponding to the isotropic manipulator posture (when the legs are mutually perpendicular and parallel to relevant actuator axes). For this assumption, the geometrical models of separate kinematic chains can be described by the expression (1) Because for the rigid manipulator the end-effector moves with only translational motions, the nominal values of the passive joint coordinates are subject to the specific constraints 3 2 4 1

; q q q q = -= -, which are implicitly incorporated in the direct/inverse kinematics [START_REF] Majou | Parametric stiffness analysis of the Orthoglide[END_REF].

However, the flexible model allows variations for all passive joints. Using the link stiffness parameters obtained by the FEA-modeling and applying the proposed methodology, we computed the compliance matrices for three typical manipulator postures, the principal components of which are presented in Table 1. Below, they are compared with the compliance of the parallelogram-based manipulator.

B. Stiffness of Parallelogram Based Manipulator

Before evaluation the compliance of the entire manipulator, let us derive the stiffness matrix of the parallelogram. Using the adopted notations, the parallelogram equivalent model may be written as

where, compared to the above case, the third passive joint is eliminated (it is implicitly assumed that 3 2 q q = -). On the other hand, the original parallelogram may be split into two serial kinematic chains (the "upper" and "lower" ones)


Hedy Attouch
email: hedy.attouch@umontpellier.fr

• Szilárd

Csaba László

Szilárd Csaba László
email: szilard.laszlo@math.utcluj.ro

Continuous Newton-like Inertial Dynamics for Monotone Inclusions

Keywords: Damped inertial dynamics, Hessian damping, maximally monotone operators, Newton method, vanishing viscosity, Yosida regularization AMS subject classification 37N40, 46N10, 49M30, 65K05, 65K10, 90B50, 90C25

In a Hilbert framework H, we study the convergence properties of a Newton-like inertial dynamical system governed by a general maximally monotone operator A : H → 2 H . When A is equal to the subdifferential of a convex lower semicontinuous proper function, the dynamic corresponds to the introduction of the Hessian-driven damping in the continuous version of the accelerated gradient method of Nesterov. As a result, the oscillations are significantly attenuated. According to the technique introduced by Attouch-Peypouquet (Math. Prog. 2019), the maximally monotone operator is replaced by its Yosida approximation with an appropriate adjustment of the regularization parameter. The introduction into the dynamic of the Newton-like correction term (corresponding to the Hessian driven term in the case of convex minimization) provides a well-posed evolution system for which we will obtain the weak convergence of the generated trajectories towards the zeroes of A. We also obtain the fast convergence of the velocities towards zero. The results tolerate the presence of errors, perturbations. Then, we specialize our results to the case where the operator A is the subdifferential of a convex lower semicontinuous function, and obtain fast optimization results.

Introduction

Let H be a real Hilbert space endowed with the scalar product •, • and norm • . Given a general maximally monotone operator A : H → 2 H , we will study the asymptotic behavior, as t → +∞, of the second-order in time evolution equation (DIN-AVD) α,β,λ,e ẍ(t) + α t ẋ(t) + β d dt A λ(t) (x(t)) + A λ(t) (x(t)) = e(t).

The operators J λA : H → H and A λ : H → H which are defined by

J λA = (I + λA) -1 and A λ = 1 λ (I -J λA ) ,
are respectively the resolvent of A and the Yosida regularization of A of index λ > 0. The coefficients α, β are positive damping parameters. The tuning of the time dependent parameter λ(t) which enters the Yosida regularization of A will play a crucial role in the asymptotic analysis. The second member e takes account of perturbations, errors. Without ambiguity, we refer briefly to the dynamic as (DIN-AVD). The terminology reflects the link of this dynamic with the Dynamic Inertial Newton method and to the Asymptotic Vanishing Damping, as explained in the next paragraph. According to the Lipschitz continuity property of the Yosida approximation, (DIN-AVD) is relevant to the Cauchy-Lipschitz theorem, which provides existence and uniqueness of the corresponding Cauchy problem. On the basis of an appropriate adjustment of the parameters, we will obtain the weak convergence of the generated trajectories towards the zeroes of A, i.e. solutions of the monotone inclusion 0 ∈ Ax.

(

) 1
This dynamic is the support of a recent study by the authors concerning rapidly converging algorithms to solve [START_REF] Abbas | Newton-like dynamics and forward-backward methods for structured monotone inclusions in Hilbert spaces[END_REF] and which make use only of the resolvents of A, see [START_REF] Attouch | Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators[END_REF]. It is a difficult problem of fundamental importance in optimization, equilibrium theory, economics and game theory, partial differential equations, statistics, among other subjects. An in-depth study of (DIN-AVD) is a key to going further in the analysis of the associated algorithms. Our study is based on several recent advances in the study of damped inertial dynamics for solving optimization problems and monotone inclusions. We describe them briefly in the following paragraphs.

Inertial dynamics and algorithms with vanishing damping coefficient

The use of inertial dynamics to accelerate the gradient method in optimization was first considered by B. Polyak in [START_REF] Polyak | Some methods of speeding up the convergence of iterative methods[END_REF]. He considered the inertial system with a fixed viscous damping coefficient γ > 0 (HBF) ẍ(t) + γ ẋ(t) + ∇f (x(t)) = 0, which, because of its natural mechanical interpretation, is called the Heavy Ball with Friction method. This system was further developed by Attouch-Goudou-Redont [START_REF] Attouch | The heavy ball with friction method. The continuous dynamical system, global exploration of the local minima of a real-valued function by asymptotical analysis of a dissipative dynamical system[END_REF] as a tool to explore the local minima of f . For a strongly convex function f , and γ judiciously chosen, (HBF) provides convergence at exponential rate of f (x(t)) to min H f . For a general convex function f , the asymptotic convergence rate of (HBF) is O(1 t ) (in the worst case). This is however not better than the steepest descent. A decisive step to obtain a faster asymptotic convergence was taken by Su-Boyd-Candès [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF] who considered the case of an Asymptotic Vanishing Damping coefficient γ(t) = α t , that is (AVD) α ẍ(t) + α t ẋ(t) + ∇f (x(t)) = 0.

For a general convex differentiable function f , and α = 3, it provides a continuous version of the accelerated gradient method of Nesterov [START_REF] Nesterov | A method for solving the convex programming problem with convergence rate O(1/k 2 ). (Russian)[END_REF]. 1 For α ≥ 3, each trajectory x(•) of (AVD) α satisfies the asymptotic convergence rate of the values

f (x(t)) -inf H f = O 1/t 2 as t → +∞.
As a specific feature, the viscous damping coefficient α t vanishes (tends to zero) as time t goes to infinity, hence the terminology. The case α = 3, which corresponds to Nesterov's historical algorithm, is critical. In the case α = 3, the question of the convergence of the trajectories remains an open problem (except in one dimension where convergence holds [START_REF] Attouch | Rate of convergence of the Nesterov accelerated gradient method in the subcritical case α ≤ 3[END_REF]). For α > 3, it has been shown by Attouch-Chbani-Peypouquet-Redont [START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF] that each trajectory converges weakly to a minimizer. For α > 3, it is shown in [START_REF] Attouch | The rate of convergence of Nesterov's accelerated forwardbackward method is actually faster than 1/k 2[END_REF] and [START_REF] May | Asymptotic for a second-order evolution equation with convex potential and vanishing damping term[END_REF] that the asymptotic convergence rate of the values is actually o(1/t 2 ). The subcritical case α ≤ 3 has been examined by Apidopoulos-Aujol-Dossal [START_REF] Apidopoulos | Convergence rate of inertial Forward-Backward algorithm beyond Nesterov's rule[END_REF] and Attouch-Chbani-Riahi [START_REF] Attouch | Rate of convergence of the Nesterov accelerated gradient method in the subcritical case α ≤ 3[END_REF], with the convergence rate of the objective values O t -2α 3 . These rates are optimal, that is, they can be reached, or approached arbitrarily close. The corresponding inertial algorithms are in line with the Nesterov accelerated gradient method

     y k = x k + 1 - α k (x k -x k-1 ) x k+1 = y k -s∇f (y k ).
They enjoy similar properties to the continuous case. They were first obtained by Chambolle-Dossal [START_REF] Chambolle | On the convergence of the iterates of the Fast Iterative Shrinkage Thresholding Algorithm[END_REF], see [START_REF] Attouch | Convergence rates of inertial forward-backward algorithms[END_REF], [START_REF] Attouch | Fast convergence of inertial dynamics and algorithms with asymptotic vanishing viscosity[END_REF], and [START_REF] László | Convergence rates for an inertial algorithm of gradient type associated to a smooth nonconvex minimization[END_REF] for further results, and the extension to proximal-gradient algorithms for structured optimization.

Hessian damping

The following inertial system combines asymptotic vanishing damping with Hessiandriven damping

ẍ(t) + α t ẋ(t) + β∇ 2 f (x(t)) ẋ(t) + ∇f (x(t)) = 0.
It was considered by Attouch-Peypouquet-Redont in [START_REF] Attouch | Fast convex minimization via inertial dynamics with Hessian driven damping[END_REF]. At first glance, the presence of the Hessian may seem to entail numerical difficulties. However, this is not the case as the Hessian intervenes in the above ODE in the form ∇2 f (x(t)) ẋ(t), which is nothing but the derivative with respect to time of ∇f (x(t)). So, the temporal discretization of this dynamic provides first-order algorithms of the form

y k = x k + α k (x k -x k-1 ) -β k (∇f (x k ) -∇f (x k-1 )) x k+1 = y k -s∇f (y k ).
As a specific feature, and by comparison with the accelerated gradient method of Nesterov, these algorithms contain a correction term which is equal to the difference of the gradients at two consecutive steps. While preserving the convergence properties of the Nesterov accelerated method, they provide fast convergence to zero of the gradients, and reduce the oscillatory aspects. Several recent studies have been devoted to this subject, see Attouch-Chbani-Fadili-Riahi [START_REF] Attouch | First-order algorithms via inertial systems with Hessian driven damping[END_REF], Bot ¸-Csetnek-László [START_REF] Bot | Tikhonov regularization of a second order dynamical system with Hessian damping[END_REF], Kim [START_REF] Kim | Accelerated Proximal Point Method for Maximally Monotone Operators[END_REF], Lin-Jordan [START_REF] Lin | A Control-Theoretic Perspective on Optimal High-Order Optimization[END_REF], Shi-Du-Jordan-Su [START_REF] Shi | Understanding the acceleration phenomenon via high-resolution differential equations[END_REF].

Inertial dynamics and cocoercive operators

Let's come to the transposition of these techniques to the case of maximally monotone operators. Álvarez-Attouch [START_REF] Alvarez | An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping[END_REF] and Attouch-Maingé [START_REF] Attouch | Asymptotic behavior of second order dissipative evolution equations combining potential with non-potential effects[END_REF] studied the equation

ẍ(t) + γ ẋ(t) + A(x(t)) = 0, ( 2
)
when A is a cocoercive 2 (and hence maximally monotone) operator, (see also [START_REF] Bot | Second order forward-backward dynamical systems for monotone inclusion problems[END_REF]).

Cocoercivity plays an important role in the study of [START_REF] Alvarez | An inertial proximal method for maximal monotone operators via discretization of a nonlinear oscillator with damping[END_REF], not only to ensure the existence of solutions, but also to analyze their long-term behavior. Assuming that the cocoercivity parameter λ and the damping coefficient γ satisfy the inequality λγ 2 > 1, they showed that each trajectory of (2) converges weakly to a zero of A. Since for λ > 0, the operator A λ is λ-cocoercive and A -1 λ (0) = A -1 (0), we immediately deduce that, under the condition λγ 2 > 1, given a general maximally monotone operator A, each trajectory of ẍ(t) + γ ẋ(t) + A λ (x(t)) = 0 converges weakly to a zero of A. In the quest for faster convergence, the system ẍ(t) + α t ẋ(t) + A λ(t) (x(t)) = 0, t > t 0 > 0, involves a time-dependent regularizing parameter λ(•) satisfying the inequality

λ(t) × α 2 t 2 > 1, t > t 0 > 0,
see Attouch-Peypouquet [START_REF] Attouch | Convergence of inertial dynamics and proximal algorithms governed by maximal monotone operators[END_REF]. Time discretization of this dynamic gives the Relaxed Inertial Proximal Algorithm (RIPA)

y k = x k + α k (x k -x k-1 ) x k+1 = (1 -ρ k )y k + ρ k J µ k A (y k ),
whose convergence properties have been analyzed by Attouch-Cabot [START_REF] Attouch | Convergence of a relaxed inertial proximal algorithm for maximally monotone operators[END_REF], Attouch-Chbani-Riahi [START_REF] Attouch | Fast proximal methods via time scaling of damped inertial dynamics[END_REF], Attouch-Peypouquet [START_REF] Attouch | Convergence of inertial dynamics and proximal algorithms governed by maximal monotone operators[END_REF].

Link with Newton-like methods for solving monotone inclusions

Let us specify the link between our study and Newton's method for solving [START_REF] Abbas | Newton-like dynamics and forward-backward methods for structured monotone inclusions in Hilbert spaces[END_REF]. To overcome the ill-posed character of the continuous Newton method for a general maximally monotone operator A, the following first order evolution system was studied by Attouch-Svaiter [START_REF] Attouch | A continuous dynamical Newton-Like approach to solving monotone inclusions[END_REF],

v(t) ∈ A(x(t)) γ(t) ẋ(t) + β v(t) + v(t) = 0.
This system can be considered as a continuous version of the Levenberg-Marquardt method, which acts as a regularization of the Newton method. Remarkably, under a fairly general assumption on the regularization parameter γ(t), this system is well posed and generates trajectories that converge weakly to equilibria. Parallel results have been obtained for the associated proximal algorithms obtained by implicit temporal discretization, see [START_REF] Abbas | Newton-like dynamics and forward-backward methods for structured monotone inclusions in Hilbert spaces[END_REF], [START_REF] Attouch | A dynamic approach to a proximal-Newton method for monotone inclusions in Hilbert Spaces, with complexity O(1/n 2 )[END_REF], [START_REF] Attouch | Global convergence of a closed-loop regularized Newton method for solving monotone inclusions in Hilbert spaces[END_REF]. Formally, this system writes as

γ(t) ẋ(t) + β d dt (A(x(t))) + A(x(t)) = 0.
Thus, (DIN-AVD) α,β can be considered as an inertial and regularized version of this dynamical system.

Organization of the paper

In Section 2, we show the existence and the uniqueness of a strong global solution to the Cauchy problem associated with (DIN-AVD). In Section 3, based on an appropriate tuning of the parameters, we study the convergence properties as t → +∞ of the trajectories generated by (DIN-AVD). Our study takes into account the presence of perturbations, errors. Section 4 is devoted to numerical experiments. In Section 5, we specialize our study in the case where the operator A is the subdifferential of a convex lower semicontinuous function. In this case, we get fast minimization properties. Finally, we present some lines of research for the future.

2 Existence and uniqueness results for (DIN-AVD)

We rely on the first-order equivalent formulation of (DIN-AVD) which is valid when β > 0. It was first considered by Alvarez-Attouch-Bolte-Redont [START_REF] Alvarez | A second-order gradient-like dissipative dynamical system with Hessian-driven damping. Application to optimization and mechanics[END_REF] and Attouch-Peypouquet-Redont [START_REF] Attouch | Fast convex minimization via inertial dynamics with Hessian driven damping[END_REF] in the framework of convex minimization, that is A = ∂f with f convex. Specifically, in our context, we have the following equivalence, which follows from a simple differential and algebraic calculation.

Proposition 1 The following are equivalent:

(i) ⇐⇒ (ii) (i) ẍ(t) + α t ẋ(t) + β d dt A λ(t) (x(t))) + A λ(t) (x(t)) = e(t).
(ii)

     ẋ(t) + βA λ(t) (x(t)) -1 β -α t x(t) + 1 β y(t) = 0; ẏ(t) -1 β -α t + αβ t 2 x(t) + 1 β y(t) = -βe(t).
Note that (ii) is different from the Hamiltonian formulation of (i). On the one hand, the formulation (ii) is written as an evolution system in the product space H × H, which is governed by the sum of a time dependent maximally monotone operator and a time dependent continuous linear operator. From this, we will deduce in the following theorem the existence and uniqueness of a strong global solution for the associated Cauchy problem. This first order equivalent formulation offers many applications. It was used by Attouch-Maingé-Redont in [START_REF] Attouch | A second-order differential system with Hessiandriven damping; Application to non-elastic shock laws[END_REF] for the modeling of damped shocks in mechanics, and by Castera-Bolte-Févotte-Pauwels in deep learning [START_REF] Castera | An Inertial Newton Algorithm for Deep Learning[END_REF].

On the other hand, when one works with (ii), one loses the mechanical interpretation of the dynamic, and the intuition of the energy notions which are attached to it. So, in the following sections, to develop a Lyapunov analysis for (DIN-AVD), we will work with the initial formulation (i). The analysis of (DIN-AVD) is essentially based on the use of time-dependent Yosida parameters λ(t). We have gathered in the following lemma some technical results concerning the properties of the mapping (λ, x) → A λ (x), which will be useful later in the proofs.

Lemma 1 Let A : H → 2 H be a maximally monotone operator, let γ, ν > 0, and x, y ∈ H. Then, the following inequalities are satisfied

a) γAγ (x) -νAν (y) ≤ 2 x -y + |γ -ν| Aγ (x) . b) Aγ (x) -Aν (y) ≤ 2 γ x -y + |γ-ν| γ ( Aγ (x) + Aν (y) ). c) Consider x : [t 0 , +∞[→ H a differentiable function. Assume further λ : [t 0 , +∞[→ ]0, +∞[ is a derivable function. Then, for every t ∈ [t 0 , +∞[ and every z ∈ A -1 (0), ( c 1 )
d dt λ(t)A λ(t) (x(t)) ≤ 2 ẋ(t) + λ (t) A λ(t) (x(t)) . (c 2 ) d dt A λ(t) (x(t)) ≤ 2 λ(t) ẋ(t) + 2 |λ (t)| λ(t) A λ(t) (x(t)) . (c 3 ) d dt λ(t)A λ(t) (x(t)) ≤ 2 ẋ(t) + |λ (t)| λ(t) x(t) -z .
Proof a) Note that according to [21, Proposition 23.28 (iii)]

J γA (x) -J νA (x) ≤ |γ -ν| Aγ (x) .
Hence,

γAγ (x) -νAν (x) = (x -y) -(J γA (x) -J νA (x) ≤ x -y + |γ -ν| Aγ (x) .
Moreover, by the 1 ν -Lipschitz continuity of the Yosida approximation Aν , we have

νAν (x) -νAν (y) ≤ x -y .
Finally,

γAγ (x) -νAν (y) = (γAγ(x) -νAν (x)) + (νAν (x) -νAν (y)) ≤ x -y + |γ -ν| Aγ (x) + x -y .
b) By using a) we have

Aγ (x) -Aν (y) = 1 γ (γAγ(x) -νAν (y)) + (ν -γ)Aν (y) ≤ 1 γ γAγ (x) -νAν (y) + |ν -γ| γ Aν (y) ≤ 2 γ x -y + |γ -ν| γ ( Aγ (x) + Aν (y) ).
(c 1 ) Let h > 0. From a) we have

λ(t + h)A λ(t+h) (x(t + h)) -λ(t)A λ(t) (x(t)) ≤ 2 x(t + h) -x(t) + |λ(t + h) -λ(t)| A λ(t+h) (x(t + h)) .
Dividing by h, and letting h → 0, we get the claim.

(c 2 ) Let h > 0. From b) we have A λ(t+h) (x(t + h)) -A λ(t) (x(t)) ≤ 2 λ(t + h) x(t + h) -x(t) + |λ(t + h) -λ(t)| λ(t + h) ( A λ(t+h) (x(t + h)) + A λ(t) (x(t)) ).
Dividing by h and letting h → 0, we get the claim.

(c 3 ) Let z ∈ A -1 (0). According to A λ(t) (z) = 0, and the 1 λ(t) -Lipschitz continuity of A λ(t) we have A λ(t) (x(t)) ≤ 1 λ(t) x(t) -z .
Combining this inequality with (c 1 ) gives the claim.

Theorem 1 Suppose that β ≥ 0. Suppose that λ : [t 0 , +∞[→]0, +∞[ is a measurable function, and that there exists λ > 0 such that λ(t) ≥ λ for all t ≥ t 0 . Suppose that e ∈ L 1 (t 0 , T ; H) for all T > t 0 . Then, for any (x 0 , x 1 ) ∈ H × H, there exists a unique strong global solution x : [t 0 , +∞[→ H of the continuous dynamic (DIN-AVD) which satisfies the Cauchy data

x(t 0 ) = x 0 , ẋ(t 0 ) = x 1 .
Proof a) Case β > 0. According to Proposition 1, it is equivalent to solve the first-order system (ii) with the Cauchy data x(t 0 ) = x 0 and y(t 0 ) = y 0 :=

-β x 1 + βA λ(t0) (x 0 ) -1 β -α t0 x 0 . This system writes Ż(t) + F (t, Z(t)) = 0, Z(t 0 ) = (x 0 , y 0 )
where

Z(t) = (x(t), y(t)) ∈ H × H and F (t, (x, y)) = βA λ(t) (x) - 1 β - α t x + 1 β y, - 1
β - α t + αβ t 2 x + 1 β y + βe(t) .
According to the Lipschitz continuity property of the Yosida approximation of a maximally monotone operator, the existence and uniqueness of a strong global solution (x, y) : [t 0 , +∞[→ H × H is relevant of the Cauchy-Lipschitz theorem. Precisely, we can apply the non-autonomous version of this theorem given in [28, Proposition 6.2.1]. Thus, we obtain a strong solution, that is, t → ẋ(t) is locally absolutely continuous. If, moreover, we suppose that the functions λ(•) and e(•) are continuous, then the solution is a classical solution of class C 2 . b) Case β = 0. We then get the system

ẍ(t) + α t ẋ(t) + A λ(t) (x(t)) = e(t).
In this case, we use the equivalent Hamiltonian formulation

   ẋ(t) -y(t) = 0; ẏ(t) + α t y(t) + A λ(t) (x(t)
) -e(t) = 0, which, thanks to the Lipschitz continuity of A λ , is relevant of the classical Cauchy-Lipschitz theorem.

3 Asymptotic convergence properties of (DIN-AVD)

Given α > 1, β ≥ 0 and e ∈ L 1 loc (t 0 , +∞; H), we consider the evolution system

(DIN-AVD) ẍ(t) + α t ẋ(t) + β d dt A λ(t) (x(t)) + A λ(t) (x(t)) = e(t), t > t 0 > 0.
The existence of strong global solutions to this system is guaranteed by Theorem 1. The convergence properties as t → +∞ of the trajectories generated by this system are summarized in the following theorem.

Theorem 2 Let A : H → 2 H be a maximally monotone operator such that S = A -1 (0) = ∅. Consider the evolution equation (DIN-AVD) where the parameters satisfy the following conditions Then, for any trajectory x : [t 0 , +∞[→ H of (DIN-AVD) the following properties are satisfied:

α > 1, β ≥ 0, and λ(t) = λt 2 with λ > 1 (α -1) 2 .

Assume further that

(i) (convergence) x(t) is bounded, and x(t) converges weakly, as t → +∞, to an element of S.

(ii) (integral estimates) +∞ t0 t ẋ(t) 2 dt < +∞, +∞ t0 t 3 ẍ(t) 2 dt < +∞ and +∞ t0 t 3 A λ(t) (x(t)) 2 dt < +∞. (iii) (pointwise estimates) lim t→+∞ ẋ(t) = 0, ẋ(t) = o 1 t as t → +∞ and, for all 0 < η < 1, one has A λ(t) (x(t)) = o 1 t 2 , d dt A λ(t) (x(t)) = o 1 t 2+η
as t → +∞.

Proof Lyapunov analysis. Take z ∈ S. For 0 < b < α -1 consider the energy functional

E b (t) := 1 2 b(x(t) -z) + t( ẋ(t) + βA λ(t) (x(t))) 2 + b(α -1 -b) 2 x(t) -z 2 . (3)
Using the classical derivation chain rule and (DIN-AVD), we get

Ėb (t) = (b + 1 -α) ẋ(t) + (β -t)A λ(t) (x(t)) + te(t), b(x(t) -z) + t( ẋ(t) + βA λ(t) (x(t))) +b(α -1 -b) ẋ(t), x(t) -z . (4)
After reduction, we obtain

Ėb (t) = b(β -t) A λ(t) (x(t)), x(t) -z + (-t 2 + β(b + 2 -α)t) A λ(t) (x(t)), ẋ(t) + (b + 1 -α)t ẋ(t) 2 + β(β -t)t A λ(t) (x(t)) 2 + bt e(t), x(t) -z + t 2 e(t), ẋ(t) + βt 2 e(t), A λ(t) (x(t)) . (5)
We have for all

p 1 > 0 bt e(t), x(t) -z ≤ bt e(t) x(t) -z , (6)
t 2 e(t), ẋ(t) ≤ p 1 t 3 e(t) 2 + t 4p 1 ẋ(t) 2 , (7)
βt 2 e(t), A λ(t) (x(t)) ≤ β 2 t 3 e(t) 2 + t A λ(t) (x(t)) 2 . (8)
Combining ( 5), ( 6), ( 7) and ( 8), we obtain

Ėb (t) ≤ b(β -t) A λ(t) (x(t)), x(t) -z + b + 1 -α + 1 4p 1 t ẋ(t) 2 (9)
+ (-t 2 + β(b + 2 -α)t) A λ(t) (x(t)), ẋ(t) + β(β -t)t + β 2 t A λ(t) (x(t)) 2 + bt e(t) x(t) -z + p 1 + β 2 t 3 e(t) 2 .
Now, using the fact that A λ(t) is λ(t) cocoercive, and that z ∈ S, we get, for all

t ≥ t 1 = max(t 0 , β) b(β -t) A λ(t) (x(t)), x(t) -z ≤ b(β -t)λ(t) A λ(t) (x(t)) 2 . ( 10
) Let us choose b = α -1 2 > 0. According to the assumption λ > 1 (α-1) 2 , we can find such that 0 < < α -1 - 1 √ λ < 2(α -1 -b), (11)
where the last inequality comes from the choice of b. Further, take p 1 = 1 . Then, ( 9) and ( 10) lead to

Ėb (t) + 4 t ẋ(t) 2 + 2 tλ(t)) A λ(t) (x(t)) 2 ≤ b + 1 -α + 2 t ẋ(t) 2 (12)
+ (-t 2 + β(b + 2 -α)t) A λ(t) (x(t)), ẋ(t) + b(β -t) + 2 t λ(t) + β(β -t)t + β 2 t A λ(t) (x(t)) 2
+ bt e(t) x(t) -z + 1 +

β 2 t 3 e(t) 2 ,
for all t ≥ t 1 . By [START_REF] Attouch | The heavy ball with friction method. The continuous dynamical system, global exploration of the local minima of a real-valued function by asymptotical analysis of a dissipative dynamical system[END_REF] we have b + 1 -α + 2 < 0. Moreover, still by [START_REF] Attouch | The heavy ball with friction method. The continuous dynamical system, global exploration of the local minima of a real-valued function by asymptotical analysis of a dissipative dynamical system[END_REF] we have

-b + 2 = -α-1 2 + 2 < 0. Since λ(t) = λt 2
with λ > 0, we deduce that there exists

t 2 ≥ t 1 such that, for all t ≥ t 2 b(β -t) + 2 t λ(t) + β(β -t)t + β 2 t < 0.
According to Lemma 2, we deduce that the sum

S(t) = b + 1 -α + 2 t ẋ(t) 2 + (-t 2 + β(b + 2 -α)t) A λ(t) (x(t)), ẋ(t) + b(β -t) + 2 t λ(t) + β(β -t)t + β 2 t A λ(t) (x(t)) 2
in the right hand side of ( 12) is nonpositive whenever

R(t) :=(-t 2 + β(b + 2 -α)t) 2 -4 b + 1 -α + 2 t b(β -t) + 2 t λ(t) + β(β -t)t + β 2 t ≤ 0.
Taking into account the fact that λ(t) = λt 2 we obtain

R(t) = 1 + 4 b + 1 -α + 2 b -2 λ t 4 + O(t 3 ). Since < α -1 -1 √ λ , we get that λ > 1 (α-1-) 2 . From b = α-1 2 we deduce that 1 + 4 b + 1 -α + 2 b -2 λ = 1 -(α -1 -) 2 λ < 0.
Consequently, there exists t 3 ≥ t 2 such that R(t) < 0, for all t ≥ t 3 . Hence, [START_REF] Attouch | Newton-like inertial dynamics and proximal algorithms governed by maximally monotone operators[END_REF] leads to, for all t ≥ t 3

Ėb (t) + 4 t ẋ(t) 2 + 2 tλ(t)) A λ(t) (x(t)) 2 ≤ bt e(t) x(t) -z + 1 + β 2 t 3 e(t) 2 . ( 13
)
Estimates. By integrating (13) on an interval [t 3 , t], and by denoting

C 0 = 1 + β 2 +∞ t3 t 3 e(t) 2 dt + E b (t 3 ) < +∞
we obtain that for all t ≥ t 3

E b (t)+ 4 t t3 s ẋ(s) 2 ds+ 2 t t3 sλ(s) A λ(s) (x(s)) 2 ds ≤ C 0 +b t t3
s e(s) x(s)-z ds.

(14) Taking into account the form of the energy functional E b (t) and the fact that b = α-1 2 , the latter relation leads to

(α -1) 2 8 x(t) -z 2 ≤ C 0 + α -1 2 t t3 s e(s) x(s) -z ds.
More precisely, we have

1 2 x(t) -z 2 ≤ C 1 + 2 α -1 t t3 s e(s) x(s) -z ds, (15)
where C 1 = 4C0 (α-1) 2 . Now, applying the Gronwall lemma (see [START_REF] Brézis | Opérateurs maximaux monotones dans les espaces de Hilbert et équations d'évolution[END_REF]Lemma A.5]) to (15) we obtain

x(t) -z ≤ 2C 1 + 2 α -1 t t3 s e(s) ds.
Therefore, x(t) -z and consequently x(t) are bounded. Further, the boundedness of x(t) -z and the assumption on e leads to +∞ t3 s e(s) x(s) -z ds < +∞.

Therefore, ( 14) becomes

E b (t) + 4 t t3 s ẋ(s) 2 ds + 2 t t3 sλ(s) A λ(s) (x(s)) 2 ds ≤ C, (16)
where

C = C 0 + b +∞ t3 s e(s) x(s) -z ds.
From this we immediately deduce that

∞ t0 t ẋ(t) 2 dt < +∞, ( 17
)
+∞ t0 t 3 A λ(t) (x(t)) 2 dt < +∞, (18)
sup t≥t0 b(x(t) -z) + t( ẋ(t) + βA λ(t) (x(t))) 2 < +∞. (19)
Moreover, the 1 λ(t) Lipschitz continuity of A λ(t) , and z ∈ S leads to

A λ(t) (x(t)) = A λ(t) (x(t)) -A λ(t) (z) ≤ 1 λ(t) x(t) -z .
Taking into account that λ(t) = λt 2 and x(t) -z is bounded, we deduce that

A λ(t) (x(t)) = O 1 t 2 , as t → +∞. ( 20
)
Further, from the boundedness of the trajectory x(•) and ( 19) we deduce that

ẋ(t) = O 1 t , as t → +∞.
In particular, we have lim t→+∞ ẋ(t) = 0. According to Lemma 1 we have

d dt λ(t)A λ(t) (x(t)) ≤ 2 ẋ(t) + 2 |λ (t)| λ(t) x(t) -z = O 1 t , as t → +∞. ( 21
)
Combining the latter relation with [START_REF] Attouch | A continuous dynamical Newton-Like approach to solving monotone inclusions[END_REF], we deduce that

λ(t) d dt A λ(t) (x(t)) = O 1 t , as t → +∞.
Consequently, we have

d dt A λ(t) (x(t)) = O 1 t 3 , as t → +∞, (22)
which implies

d dt A λ(t) (x(t)) = o 1 t 2+η , as t → +∞, for all 0 < η < 1.
Let us improve the estimate [START_REF] Attouch | A continuous dynamical Newton-Like approach to solving monotone inclusions[END_REF], and show that

A λ(t) (x(t)) = o 1 t 2 , as t → +∞.
We have

d dt λ(t)A λ(t) (x(t)) 4 = 4 λ(t)A λ(t) (x(t)), d dt (λ(t)A λ(t) (x(t))) λ(t)A λ(t) (x(t)) 2 . ( 23
)
According to [START_REF] Attouch | A continuous dynamical Newton-Like approach to solving monotone inclusions[END_REF] and ( 21) there exists K > 0 such that

λ(t)A λ(t) (x(t)), d dt (λ(t)A λ(t) (x(t))) ≤ λ(t)A λ(t) (x(t)) d dt (λ(t)A λ(t) (x(t))) ≤ K t .
Hence, (23) leads to

d dt λ(t)A λ(t) (x(t)) 4 ≤ 4K t λ(t)A λ(t) (x(t)) 2 . ( 24
)
According to [START_REF] Attouch | Fast convex minimization via inertial dynamics with Hessian driven damping[END_REF] the right hand side of (24) belongs to L 1 (t 0 , +∞), which implies

d dt λ(t)A λ(t) (x(t)) 4 ∈ L 1 (t 0 , +∞).
Therefore lim t→+∞ λ(t)A λ(t) (x(t)) 4 exists.

But then, L := lim t→+∞ λ(t)A λ(t) (x(t)) 2 also exists. Using (18) again, i.e.

+∞ t0 1 t λ(t)A λ(t) (x(t)) 2 dt = λ +∞ t0 tλ(t) A λ(t) (x(t)) 2 dt < +∞,
we deduce that L = 0. Therefore, lim t→+∞ λ(t)A λ(t) (x(t)) 2 = 0, which gives

A λ(t) (x(t)) = o 1 t 2 , as t → +∞. ( 25
)
Finally, by using (DIN-AVD) we have

ẍ(t) 2 = e(t) - α t ẋ(t) -β d dt A λ(t) (x(t)) -A λ(t) (x(t)) 2 .
According to an elementary convexity inequality, we get

t 3 ẍ(t) 2 ≤ 4t 3 e(t) 2 + 4α 2 t ẋ(t) 2 + 4β 2 t 3 d dt A λ(t) (x(t)) 2 + 4t 3 A λ(t) (x(t)) 2 .
According to

+∞ t0
t 3 e(t) 2 dt < +∞, ( 17), ( 22) and ( 18) we obtain that

+∞ t0 t 3 ẍ(t) 2 dt < +∞. ( 26
)
Let us now prove that ẋ(t) = o 1 t , as t → +∞. We have

d dt t 2 ẋ(t) 2 = 2t ẋ(t) 2 + 2t 2 ẍ(t), ẋ(t) 2t 2 ẍ(t), ẋ(t) ≤ t 3 ẍ(t) 2 + t ẋ(t) 2 .
Hence,

d dt t 2 ẋ(t) 2 ≤ t 3 ẍ(t) 2 + 3t ẋ(t) 2 .
According to [START_REF] Castera | An Inertial Newton Algorithm for Deep Learning[END_REF] and [START_REF] Attouch | The rate of convergence of Nesterov's accelerated forwardbackward method is actually faster than 1/k 2[END_REF] we have t 3 ẍ(t) 2 + 3t ẋ(t) 2 ∈ L 1 (t 0 , +∞). Therefore, from [1, Lemma 5.1] there exists lim t→+∞ t 2 ẋ(t) 2 ∈ R. Using (17) again, we have

∞ t0 1 t (t 2 ẋ(t) 2 )dt = ∞ t0 t ẋ(t) 2 dt < +∞.
Therefore, lim t→+∞ t 2 ẋ(t) 2 = 0, and ẋ(t) = o 1 t as t → +∞.

The limit. To prove the existence of the weak limit of x(t), we use Opial lemma. Let us introduce the anchor function hz(t) = 1 2 x(t) -z 2 . The classical derivation chain rule gives

ḧz(t) + α t ḣz(t) = ẍ(t) + α t ẋ(t), x(t) -z + ẋ(t) 2 .

By using (DIN-AVD) we get

β d dt A λ(t) (x(t)), x(t) -z = e(t) -ẍ(t) - α t ẋ(t) -A λ(t) (x(t)), x(t) -z = e(t) -A λ(t) (x(t)), x(t) -z -ḧz(t) + α t ḣz(t) + ẋ(t) 2 .
Therefore

t ḧz(t) + α ḣz(t) + t A λ(t) (x(t)), x(t) -z (27) = t ẋ(t) 2 + t e(t), x(t) -z -βt d dt A λ(t) (x(t)), x(t) -z .
By Cauchy-Schwarz inequality

-βt d dt A λ(t) (x(t)), x(t) -z ≤ βt d dt A λ(t) (x(t)) x(t) -z , t e(t), x(t) -z ≤ t e(t) x(t) -z .
According to [START_REF] Chambolle | On the convergence of the iterates of the Fast Iterative Shrinkage Thresholding Algorithm[END_REF], we deduce that

t ḧz(t) + α ḣz(t) + t A λ(t) (x(t)), x(t) -z ≤ t ẋ(t) 2 + βt d dt A λ(t) (x(t)) x(t) -z + t e(t) x(t) -z .
According to x(t) -z is bounded, t e(t) ∈ L 1 (t 0 , +∞), ( 17) and ( 22) we have

t ẋ(t) 2 + βt d dt A λ(t) (x(t)) x(t) -z + t e(t) x(t) -z ∈ L 1 (t 0 , +∞).
Moreover t → t A λ(t) (x(t)), x(t) -z is a non-negative function. So, we can apply Lemma A.6 from [START_REF] Attouch | Convergence of inertial dynamics and proximal algorithms governed by maximal monotone operators[END_REF], and obtain that lim t→+∞ hz(t) exists. In other words, lim t→+∞ x(t) -z exists for all z ∈ S.

To complete the proof via the Opial's lemma, we need to prove that every weak sequential cluster point of x(t) belongs to S. To this end, we use the following property of the Yosida approximation of a maximally monotone operator:

A λ (x) ∈ A(x -λA λ (x)
), for all x ∈ H and λ > 0.

Let tn → +∞ such that x(tn) x * , n → +∞. Since the graph of A is demi-closed, by using [START_REF] Bot | Tikhonov regularization of a second order dynamical system with Hessian damping[END_REF] we have

0 = lim n→+∞ A λ(tn) (x(tn)) ∈ A( lim n→+∞ (x(tn) -λ(tn)A λ(tn) (x(tn))) = A(x * ).
Consequently, x(t) converges weakly to an element of S.

Some numerical experiments

Take A : R 2 → R 2 , A(x, y) = (-y, x), which is a linear skew symmetric operator. Clearly, A is a maximally monotone whose single zero is x * = (0, 0). Further, A can be identified with the matrix 0 -1 1 0

, and an easy computation shows that the Yosida regularization of A can be identified with the matrix

A λ = λ 1+λ 2 -1 1+λ 2 1 1+λ 2 λ 1+λ 2 .
As a standing assumption, throughout the following numerical experiments, we take the coefficients λ(t) of the form λ(t) = λt 2 , and the perturbation errors of the form e(t) = e 1 t p , e 2 t q , p, q > 2. Further, we consider the starting points u 0 = (1, 1), v 0 = (1, 1). Note that this choice of the initial speed vector v 0 means that, at the start, the trajectory tends to move away from the origin. Obviously, a trajectory of (DIN-AVD) in this case is of the form x(t) = (x 1 (t), x 2 (t)). In order to solve the dynamical system (DIN-AVD) with starting points x(t 0 ) = u 0 , ẋ(t 0 ) = v 0 on an interval [t 0 , T ], we use the MATLAB function ode45. To this purpose we rewrite (DIN-AVD) as the first order system used in Proposition 1, that is,

                                   ẋ1 (t) = 1 β -α t -βλt 2 1+λ 2 t 4 x 1 (t) + β 1+λ 2 t 4 x 2 (t) -1 β y 1 (t) ẋ2 (t) = -β 1+λ 2 t 4 x 1 (t) + 1 β -α t -βλt 2 1+λ 2 t 4 x 2 (t) -1 β y 2 (t) ẏ1 (t) = 1 β -α t + αβ t 2 x 1 (t) -1 β y 1 (t) -βe1 t p ẏ2 (t) = 1 β -α t + αβ t 2 x 2 (t) -1 β y 2 (t) -βe2 t q (x 1 (t 0 ), x 2 (t 0 ), y 1 (t 0 ), y 2 (t 0 )) = 1, 1, -β 2 λt 2 0 -1 1+λ 2 t 4 0 -αβ t0 -β + 1, -β 2 λt 2 0 +1 1+λ 2 t 4 0 -αβ t0 -β + 1 . (29)
In the case β = 0, we rewrite (DIN-AVD) as in Theorem 1 (b), that is

                         ẋ1 (t) = y 1 (t) ẋ2 (t) = y 2 (t) ẏ1 (t) = -λt 2 1+λ 2 t 4 x 1 (t) + 1 1+λ 2 t 4 x 2 (t) -α t y 1 (t) + e1 t p ẏ2 (t) = -1 1+λ 2 t 4 x 1 (t) -λt 2 1+λ 2 t 4 x 2 (t) -α t y 2 (t) + e2 t q (x 1 (t 0 ), x 2 (t 0 ), y 1 (t 0 ), y 2 (t 0 ) = (1, 1, 1 , 1). (30)
I. In our first experiment, we are interested in the asymptotic behavior of the components x 1 and x 2 , obtained by taking different values of the parameters α, β, λ, e 1 , e 2 , p, q. According to Theorem 2 (i), lim t→+∞ (x 1 (t), x 2 (t)) = (0, 0), under the conditions

α > 1, β ≥ 0, λ > 1 (α -1) 2 , +∞ t0 t 3 e(t) 2 dt < +∞, +∞ t0 t e(t) dt < +∞.
We are interested in the gain that the term β d dt A λ(t) (x(t)) brings in (DIN-AVD) (considered with error term or without error term). We consider the following cases which fit the assumptions of Theorem 2 (recalled just above). The trajectories obtained by solving ( 29) and ( 30) with the ode45 function in Matlab on the interval [0.1, 50] are depicted at Figure 1, (a)-(d), where we represent the component x 1 (t) with red and x 2 (t) with black. In our second experiment, we examine the rapid convergence of the speed towards 0. According to Theorem 2, we have ẋ(t) = o 1 t as t → +∞. So, next to ẋ(t) , we also represent the entity t ẋ(t) . We solve ( 29) and ( 30) with the ode45 function in Matlab on the interval [0.1, 50] by considering the following instances. Note that for these values also, the hypotheses of Theorem 2 are verified, and, therefore, its conclusions are valid. The results obtained are depicted at Figure 2, (a)-(d), where we represent t ẋ(t) with red and ẋ(t) with black. A conclusion that these experiment give, is that indeed the term β d dt A λ(t) (x(t)) in (DIN-AVD) has an accelerating effect, even in the presence of an error term e. III. In our third experiment, we are interested in the behavior of the error x(t) -x * . Since the operator A has a single zero at x * = (0, 0), according to Theorem 4, x(t) must converge to 0. Our experiment shows the importance of the correction term β d dt A λ(t) (x(t)) in (DIN-AVD) as well as the sensitivity of the trajectories with respect to the error e.

Hence, we solve ( 29) and ( 30) with the ode45 function in Matlab on the interval [0. [START_REF] Abbas | Newton-like dynamics and forward-backward methods for structured monotone inclusions in Hilbert spaces[END_REF]50] where we take at first the error e(t) ≡ 0 that is e 1 = e 2 = 0. For the values β ∈ {0, 0.25, 0.5, 0.75, 1} we consider (α, λ) = (2.5, 0.5) depicted at Figure 3,(a

= q = 3, that is e(t) = 1 t 3 , 1 t 3 .
For the values β ∈ {0, 0.25, 0.5, 0.75, 1} we consider (α, λ) = (2.5, 0.5) depicted at Figure 4,(a Besides the importance of the correction term (β > 0), these numerical experiments show that the trajectories of (DIN-AVD) are quite sensitive to the parameters α, β, λ as well as to the error e(t).

The convex case

Let us specialize the previous results to the case of convex minimization, and show the rapid convergence of values. Given a lower semicontinuous convex and proper function f : H → R ∪ {+∞} such that argmin f = ∅, we consider the minimization problem (P) inf

x∈H f (x).
Fermat's rule states that x is a global minimum of f if and only if

0 ∈ ∂f (x). ( 31
)
Therefore, (P) is equivalent to the monotone inclusion problem (1) with A = ∂f . Moreover, argmin f = (∂f ) -1 (0). The Yosida approximation of ∂f is equal to the gradient of the Moreau envelope of f : for any λ > 0

(∂f ) λ = ∇f λ . ( 32
)
Recall that f λ : H → R is a C 1,1 function, which is defined by: for any x ∈ H

f λ (x) = inf ξ∈H f (ξ) + 1 2λ x -ξ 2 .
We deduce at once that min x∈H f λ (x) = min x∈H f (x). In this context, (DIN-AVD) reads as, for t ≥ t 0 > 0

(DIN-AVD)-convex ẍ(t) + α t ẋ(t) + β d dt ∇f λ(t) (x(t)) + ∇f λ(t) (x(t)) = e(t),
where α > 1, β ≥ 0 and e ∈ L 1 loc (t 0 , +∞; H). As a direct consequence of Theorem 2 we have the following result. Then, for any trajectory x : [t 0 , +∞[→ H of (DIN-AVD)-convex the following properties are satisfied: (i) (convergence) x(t) is bounded, and x(t) converges weakly, as t → +∞, to an element of S.

(ii) (integral estimates)

+∞ t0 t ẋ(t) 2 dt < +∞, +∞ t0 t 3 ẍ(t) 2 dt < +∞ and +∞ t0 t 3 ∇f λ(t) (x(t)) 2 dt < +∞. (iii) (pointwise estimates) lim t→+∞ ẋ(t) = 0, ẋ(t) = o 1 t as t → +∞
and, for all 0 < η < 1, as t → +∞ we have

∇f λ(t) (x(t)) = o 1 t 2 , d dt ∇f λ(t) (x(t)) = o 1 t 2+η .
(iv) (fast convergence of the values) As t → +∞

f λ(t) (x(t)) -min H f = o 1 t 2 and f (prox λ(t)f (x(t))) -min H f = o 1 t 2 .
In addition, prox λ(t)f (x(t)) -x(t) → 0 as t → +∞.

Proof (i)-(iii) follow directly from Theorem 2 applied to ∂f and using [START_REF] May | Asymptotic for a second-order evolution equation with convex potential and vanishing damping term[END_REF]. (iv) Take x * ∈ argmin f. From the gradient inequality, and x(t) bounded, we have

f λ(t) (x(t)) -min H f = f λ(t) (x(t)) -f λ(t) (x * ) ≤ ∇f λ(t) (x(t)), x(t) -x * ≤ ∇f λ(t) (x(t)) x(t) -x * ≤ M ∇f λ(t) (x(t)) ,
where M := sup t≥t0 x(t) -x * . Combining the above relation with the estimate obtained in ( 25)

), ∇f λ(t) (x(t)) = o 1 t 2
as t → +∞, we obtain

f λ(t) (x(t)) -min H f = o 1 t 2 as t → +∞. ( 33
)
By definition of f λ(t) and of the proximal mapping, we have

f λ(t) (x(t)) -min H f = f (prox λ(t)f (x(t))) -min H f + 1 2λ(t) x(t) -prox λ(t)f (x(t)) 2 .
(34) Combining ( 33) with [START_REF] Polyak | Some methods of speeding up the convergence of iterative methods[END_REF], we obtain, as t → +∞

f (prox λ(t)f (x(t))) -min H f = o 1 t 2 , lim t→+∞ t 2 1 2λ(t) x(t) -prox λ(t)f (x(t)) 2 = 0.
Therefore, lim t→+∞ x(t) -prox λ(t)f (x(t)) = 0, which completes the proof.

Remark 1 When A = ∂f , f convex, we have additional tools, such as the gradient inequality. We will show in the following theorem that, in this case, some assumptions can be weakened.

Theorem 4 Let f : H → R∪{+∞} be a lower semicontinuous convex, proper function such that S = argmin H f = ∅. Assume that the parameters of (DIN-AVD) satisfy α > 3, β ≥ 0 and λ(t) = λt r , with λ > 0, r ≥ 0.

Suppose that the error terms satisfy the integrability properties

+∞ t0 t 3 e(t) 2 dt < +∞ and +∞ t0
t e(t) dt < +∞.

Then, for any trajectory x : [t 0 , +∞[→ H of (DIN-AVD)-convex the following properties are satisfied.

a) In the general case of β ≥ 0, r ≥ 0 we have the following properties:

(pointwise estimates) The trajectory x(t) is bounded, and as t → +∞

f λ(t) (x(t)) -min H f = O 1 t 2 . f (prox λ(t)f (x(t)))-min H f = O 1 t 2 , and
x(t)-prox λ(t)f (x(t)) = O √ λ(t) t . Further we have ∇f λ(t) (x(t)) = O 1 t √ λ(t)
, and ẋ(t) = O 1 t .

Whenever r = 2, one has ẋ(t) = o 1 t as t → +∞.

(integral estimates)

+∞ t0 t(f λ(t) (x(t)) -min H f )dt < +∞, ∞ t0 t ẋ(t) 2 dt < +∞, ∞ t0 tλ(t) ∇f λ(t) (x(t)) 2 dt < +∞ and β ∞ t0 t 2 ∇f λ(t) (x(t)) 2 dt < +∞. Further, if r ≤ 2 +∞ t0 tλ(t) ẍ(t) 2 dt < +∞.
b) In the case β = 0 or β > 0 and r > 1 we have the following properties:

f λ(t) (x(t))-min H f = o 1 t 2 , as t → +∞, f (prox λ(t)f (x(t)))-min H f = o 1 t 2 , x(t)-prox λ(t)f (x(t)) = o √ λ(t) t as t → +∞. Further we have ∇f λ(t) (x(t)) = o 1 t √ λ(t) as t → +∞. Moreover, ẋ(t) = o 1 t as t → +∞.
c) In the case β = 0 and r ∈ [0, 2] or β > 0 and r ∈]1, 2] the trajectory x(t) converges weakly, as t → +∞, to an element of S.

Proof Lyapunov analysis. Take z ∈ S. Since α > 3 we can take 2 < b < α -1. In this particular case, the energy functional (3) becomes

E b (t) = 1 2 b(x(t) -z) + t( ẋ(t) + β∇f λ(t) (x(t))) 2 + b(α -1 -b) 2 x(t) -z 2 . ( 35
)
Using the same arguments as in the proof of Theorem 2, we obtain ( 9)), which in this particular case reads as

Ėb (t) ≤ b(β -t) ∇f λ(t) (x(t)), x(t) -z + b + 1 -α + 1 4p 1 t ẋ(t) 2 (36)
+ (-t 2 + β(b + 2 -α)t) ∇f λ(t) (x(t)), ẋ(t) + β(β -t)t + β 2 t ∇f λ(t) (x(t)) 2
+ bt e(t) x(t) -z + p 1 +

β 2 t 3 e(t) 2 .
Further we have

(-t 2 + β(b + 2 -α)t) ∇f λ(t) (x(t)), ẋ(t) = d dt (-t 2 + β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) + (2t -β(b + 2 -α))(f λ(t) (x(t)) -min H f ). ( 37
)
(-t 2 + β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) +(b(β -t) + t + 2t -β(b + 2 -α))(f λ(t) (x(t)) -min H f ) -t ∇f λ(t) (x(t)), x(t) -z . (39)
Combining ( 36) and (39) we get, for all t ≥ t 1

Ėb (t) + d dt (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) + t ∇f λ(t) (x(t)), x(t) -z ≤ ((-b + + 2)t -β(2 -α))(f λ(t) (x(t)) -min H f ) + b + 1 -α + 1 4p 1 t ẋ(t) 2 + β(β -t)t + β 2 t ∇f λ(t) (x(t)) 2 + bt e(t) x(t) -z + p 1 + β 2 t 3 e(t) 2 .
(40)

By using the cocoerciveness of ∇f λ(t) and the gradient inequality, we obtain

t ∇f λ(t) (x(t)), x(t) -z ≥ 2 tλ(t) ∇f λ(t) (x(t)) 2 + 2 t(f λ(t) (x(t)) -min H f ).
Further, take

1 < 2(α -1 -b), p 1 = 1 1 , 0 < 2 < 1.
Then, (40) leads to

Ėb (t) + d dt (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) (41) + 2 tλ(t) ∇f λ(t) (x(t)) 2 + 2 t(f λ(t) (x(t)) -min H f ) + 1 4 t ẋ(t) 2 + β 2 t 2 ∇f λ(t) (x(t)) 2 ≤ ((-b + + 2)t -β(2 -α))(f λ(t) (x(t)) -min H f ) + b + 1 -α + 1 2 t ẋ(t) 2 + β(β -t)t + β 2 t 2 + β 2 t ∇f λ(t) (x(t)) 2
+ bt e(t) x(t) -z + 1

1 + β 2 t 3 e(t) 2 , for all t ≥ t 1 .
Now, obviously there exists t 2 ≥ t 1 such that for all t ≥ t 2 one has

t 2 -β(b + 2 -α)t > 0, (-b + + 2)t -β(2 -α) ≤ 0 β(β -t)t + β 2 t 2 + β 2 t ≤ 0.
Hence, for all t ≥ t 2 , it holds

Ėb (t) + d dt (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) + 2 tλ(t) ∇f λ(t) (x(t)) 2 + 2 t(f λ(t) (x(t)) -min H f ) + 1 4 t ẋ(t) 2 + β 2 t 2 ∇f λ(t) (x(t)) 2 ≤
+ bt e(t) x(t) -z + 1

1 + β 2 t 3 e(t) 2 .
(42)

Estimates. By integrating (42) on an interval [t 2 , t], and by denoting

C 0 = 1 + β 2 +∞ t2 t 3 e(t) 2 dt+E b (t 2 )+(t 2 2 -β(b+2-α)t 2 )(f λ(t2) (x(t 2 ))-min H f )
we obtain that, for all t ≥ t 2

E b (t) + (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) + 2 t t2 sλ(s) ∇f λ(s) (x(s)) 2 ds + 2 t t2 s(f λ(s) (x(s)) -min H f )ds + 1 4 t t2 s ẋ(s) 2 ds + β 2 t t2 s 2 ∇f λ(s) (x(s)) 2 ds ≤ C 0 + b t t2 s e(s) x(s) -z ds. (43)
Taking into account the form of the energy functional E b (t), (43) leads to

b(α -1 -b) 2 x(t) -z 2 ≤ C 0 + b t t2
s e(s) x(s) -z ds.

More precisely, we have

1 2 x(t) -z 2 ≤ C 1 + 1 α -1 -b t t2 s e(s) x(s) -z ds, (44)
where C 1 = C0 b(α-1-b) . Now, applying the Gronwall lemma (see [START_REF] Brézis | Opérateurs maximaux monotones dans les espaces de Hilbert et équations d'évolution[END_REF]Lemma A.5]) to (44), we obtain

x(t) -z ≤ 2C 1 + 1 α -1 -b t t2 s e(s) ds.
Therefore, x(t)-z and consequently x(t) are bounded. Further, the boundedness of x(t) -z leads to +∞ t2 s e(s) x(s) -z ds < +∞.

Consequently, (43) becomes

E b (t) + (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) + 2 t t2 sλ(s) ∇f λ(s) (x(s)) 2 ds + 2 t t2 s(f λ(s) (x(s)) -min H f )ds + 1 4 t t2 s ẋ(s) 2 ds + β 2 t t2 s 2 ∇f λ(s) (x(s)) 2 ds ≤ C := C 0 + b +∞ t2 s e(s) x(s) -z ds. (45)
From (45) we get

sup t≥t0 b(x(t) -z) + t( ẋ(t) + β∇f λ(t) (x(t))) < +∞ (46) f λ(t) (x(t)) -min H f = O 1 t 2 , as t → +∞. ( 47
) +∞ t0 t(f λ(t) (x(t)) -min H f )dt < +∞. ( 48
) ∞ t0 t ẋ(t) 2 dt < +∞. ( 49
) ∞ t0 tλ(t) ∇f λ(t) (x(t)) 2 dt < +∞. ( 50
) β ∞ t0 t 2 ∇f λ(t) (x(t)) 2 dt < +∞. (51)
Note that (51) provides information only in the case β > 0. Now, (46) gives

ẋ(t) + β∇f λ(t) (x(t)) = O 1 t , as t → +∞. (52)
Further, from (47) and

f λ(t) (x(t)) -min H f = f (prox λ(t)f (x(t))) -min H f + 1 2λ(t) x(t) -prox λ(t)f (x(t)) 2 ,
we deduce that, as t → +∞

f (prox λ(t)f (x(t)))-min H f = O 1 t 2 , x(t)-prox λ(t)f (x(t)) = O λ(t) t (53)
Further we have ∇f λ Finally, by using (DIN-AVD)-convex, we have

(t) = (∂f ) λ(t) = 1 λ(t) (I -prox λ(t)f ), hence ∇f λ(t) (x(t)) = O 1 t λ(t) as t → +∞. ( 54
ẍ(t) 2 = e(t) - α t ẋ(t) -β d dt ∇f λ(t) (x(t)) -∇f λ(t) (x(t)) 2 .
Therefore,

tλ(t) ẍ(t) 2 ≤ 4tλ(t) e(t) 2 + 4α 2 λ(t) t + 32β 2 t λ(t) ẋ(t) 2 + 32β 2 λ 2 (t)t λ(t) + 4tλ(t) A λ(t) (x(t)) 2 .
Recall that λ(t) = λt r and assume that r ≤ 2. Then, according to the fact that

t0 1 t (t 1+ r 2 ẋ(t) 2 )dt = ∞ t0 t r 2 ẋ(t) 2 dt < +∞, hence, lim t→+∞ t 1+ r 2 ẋ(t) 2 = 0. Consequently, ẋ(t) = o 1 t 1 2 + r 4
, as t → +∞.

The limit. To prove the existence of the weak limit of x(t), we use Opial lemma. Let us introduce the anchor function

hz(t) = 1 2 x(t) -z 2 .
The classical derivation chain rule gives

ḧz(t) + α t ḣz(t) = ẍ(t) + α t ẋ(t), x(t) -z + ẋ(t) 2 .
By using (DIN-AVD)-convex, we get

β d dt ∇f λ(t) (x(t)), x(t) -z = e(t) -ẍ(t) - α t ẋ(t) -∇f λ(t) (x(t)), x(t) -z = e(t) -∇f λ(t) (x(t)), x(t) -z -ḧz(t) + α t ḣz(t) + ẋ(t) 2 .
Therefore

t ḧz(t) + α ḣz(t) + t ∇f λ(t) (x(t)), x(t) -z = t ẋ(t) 2 + t e(t), x(t) -z -βt d dt ∇f λ(t) (x(t)), x(t) -z . (57)
According to Cauchy-Schwarz inequality It remains to show that every weak sequential cluster point of the trajectory x(t) belongs to argmin f. Let x * be a weak sequential cluster point of x(t). Then, there exists a sequence tn → +∞, n → +∞ such that x(tn) x * , n → +∞. According to (59), if r ≤ 2, we have lim n→+∞ x(tn) -prox λ(tn)f (x(tn)) = 0. Therefore, prox λ(tn)f (x(tn))

x * , n → +∞.

Since f is lower semicontinuous and convex, it is weakly lower semicontinuous.

Combined with lim n→+∞ (f (prox λ(tn)f (x(tn))) -min H f ) = 0, it yields 0 = lim inf x(tn) x * (f (prox λ(tn)f (x(tn))) -min

H f ) ≥ f (x * ) -min H f.
The latter relation shows that x * ∈ argmin f . Consequently, according to Opial lemma, x(t) converges weakly to an element x ∈ argmin f as t → +∞.

Conclusion, perspective

Recent developments in convex optimization show the importance of the introduction of the Hessian driven damping in the continuous versions of the Nesterov accelerated gradient method. It allows to control and attenuate the oscillations, resulting in faster methods. The extension of these results to general monotone inclusions is an important and non-trivial question. Our study of the continuous dynamic (DIN-AVD) gives a solid mathematical basis to the algorithmic results obtained by the authors for these questions, and confirm them. Dealing with these issues in the context of general maximally monotone operators offers a wide range of applications. It is a natural idea for further study to specialize this study in the case of convex-concave saddle value problems. The convergence results are valid in the presence of perturbations or errors. This is an important step to study other instances. Among them, the introduction of a Tikhonov regularization term with vanishing coefficient, in order to asymptotically obtain the solution of minimum norm. It also suggests developing stochastic versions of (DIN-AVD) and corresponding algorithms in the context of general maximally monotone operators. To deal with concrete examples, it would be very interesting to develop corresponding splitting methods to resolve structured monotone inclusions. Finally, it would also be interesting to consider the closed loop version of these dynamics and algorithms where the coefficient λ of the Yosida regularization is taken as a feedback control of the state or the velocity of the system.

+∞ t0 t 3

 3 e(t) 2 dt < +∞ and +∞ t0 t e(t) dt < +∞.

  Figure 2.5 0 0.5 0 0 --Fig.1 (a) 2.5 0.5 0.5 0 0 --Fig.1 (b) 2.5 0 0.5 1 1 3.1 3.1 Fig.1 (c) 2.5 0.5 0.5 1 1 3.1 3.1 Fig.1 (d)

Fig. 1 :

 1 Fig. 1: Trajectories of (DIN-AVD) for different instances of the parameters α, β, λ, e

Figure 2

 2 Figure 2.1 0 1 0 0 --Fig.2 (a) 2.1 0 1 1 1 3 3 Fig.2 (b) 2.1 0.25 1 0 0 --Fig.2 (c) 2.1 0.75 1 1 1 3 3 Fig.2 (d)

Fig. 2 :

 2 Fig. 2: Fast convergence of the velocity for different instances of α, β, λ, e

Fig. 3 :

 3 Fig. 3: Comparison of the iteration error x(t) -x * for different instances of (DIN-AVD) without error

  ) and λ) = (3.1, 1) depicted at Figure 4, (b), respectively.

Fig. 4 :

 4 Fig. 4: Comparison of the iteration error x(t) -x * for different instances of (DIN-AVD) with error

Theorem 3

 3 Let f : H → R∪{+∞} be a lower semicontinuous convex proper function such that S = argmin f = ∅. Assume further that λ(t) = λt 2 with λ > 1 (α -1) 2 and +∞ t0 t 3 e(t) 2 dt < +∞ and +∞ t0 t e(t) dt < +∞.

∞

  Take 0 < < b -2. Then there exists t 1 ≥ t 0 such that b(β -t) + t < 0 for all t ≥ t 1 . So, by the gradient inequality we get, for all t ≥ t 1(b(β -t) + t) ∇f λ(t) (x(t)), x(t) -z ≤ (b(β -t) + t)(f λ(t) (x(t)) -min

		H	f ). (38)
	By adding (37) with (38), we obtain, for all t ≥ t 1
	b(β -t) ∇f λ(t) (x(t)), x(t) -z + (-t 2 + β(b + 2 -α)t) ∇f λ(t) (x(t)), ẋ(t)
	≤	d dt

  ≤ t 1+r ẍ(t) 2 + t ẋ(t) 2 .

				+∞
				tλ(t) ẍ(t) 2 dt < +∞.	(56)
				t0
	Let us now prove that ẋ(t) = o	t	1 2 + r 1 4	, as t → +∞. We have
	d dt	t 1+ r 2	ẋ(t) 2 = 1 +	r 2	t	r 2	ẋ(t) 2 + 2t 1+ r 2 ẍ(t), ẋ(t)
	and 2 ẍ(t), ẋ(t) Hence, 2t 1+ r
	d dt	t 1+ r 2	ẋ(t) 2 ≤ t 1+r ẍ(t) 2 + 1 +	r 2	t	r 2 + t ẋ(t) 2 .
	According to (56)-(49), we have t 1+r ẍ(t) 2 + 1 + r 2 t	r 2 + t ẋ(t) 2 ∈ L 1 [t 0 , +∞[.
	Consequently, from [1, Lemma 5.1] we obtain that there exists
				lim t→+∞	t 1+ r 2	ẋ(t) 2 ∈ R.
	Now, using (49) again, we have

+∞ t0

t 3 e(t) 2 dt < +∞, (49) and (

50

) we obtain that

  Now, if β = 0 or β > 0 and r > 1 then, according to the fact that x(t) -z is bounded, t e(t) ∈ L 1 (t 0 , +∞), (49), (55) and (54) we have ∈ L 1 (t 0 , +∞). Moreover t ∇f λ(t) (x(t)), x(t) -z is a non-negative function. So, we can apply Lemma A.6 from[START_REF] Attouch | Convergence of inertial dynamics and proximal algorithms governed by maximal monotone operators[END_REF], and obtain that lim t→+∞ hz(t) exists. In other words, lim t→+∞ x(t) -z exists for all z ∈ S.Let us return to the fact that, according to (42), for allt ≥ t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -minThe right hand side of the latter inequality belongs to L 1 (t 0 , +∞). Therefore, according to [1, Lemma 5.1] we get that there exists lim prox λ(t)f (x(t))2 , we deduce that, as t → +∞f (prox λ(t)f (x(t))) -min H f = o 1 t 2 , x(t) -prox λ(t)f (x(t)) = oFurther we have ∇f λ(t) = (∂f ) λ(t) = 1 λ(t) (Iprox λ(t)f ), hence

	-βt t e(t), x(t) -z ≤ t e(t) x(t) -z . d dt ∇f λ(t) (x(t)), x(t) -z ≤ βt So, we deduce from (57) that Now, according to Lemma 1 d) we have d dt ∇f λ(t) (x(t)) ≤ 2 λ(t) ẋ(t) + 2 dt d |λ (t)| ∇f λ(t) (x(t)) x(t) -z , λ(t) ∇f λ(t) (x(t)) . Therefore, t ḧz(t) + α ḣz(t) + t ∇f λ(t) (x(t)), x(t) -z ≤ t ẋ(t) 2 + t e(t) x(t) -z + 2βt λ(t) ẋ(t) + 2β |λ (t)|t λ(t) ∇f λ(t) (x(t)) t ẋ(t) 2 + t e(t) + 2βt λ(t) ẋ(t) + 2β |λ (t)|t λ(t) ∇f λ(t) (x(t)) x(t)-z Ėb (t) + d dt (t 2 H f ) ≤ bt e(t) x(t) -z + 1 1 + β 2 t 3 e(t) 2 . (58) x(t) -z . Consequently, lim t→+∞ 1 2 H f = o 1 t 2 , and ẋ(t) = o 1 t . From f λ(t) (x(t)) -min H f = f (prox λ(t)f (x(t))) -min H f + 1 2λ(t) t (59) x(t) λ(t) ∇f λ(t) (x(t)) = o 1 t λ(t) as t → +∞. (60)

t ḧz(t) + α ḣz(t) + t ∇f λ(t) (x(t)), x(t) -z

≤ t ẋ(t) 2 + βt d dt ∇f λ(t) (x(t)) x(t) -z + t e(t) x(t) -z . t→+∞ E b (t) + (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) ∈ R.

Since lim t→+∞ x(t) -z exists for all z ∈ S, we obtain that there exists

lim t→+∞ 1 2 t( ẋ(t) + β∇f λ(t) (x(t))) 2 + (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) ∈ R.

On the other hand, from (48), (

49

) and (51) we get

+∞ t0 1 t 1 2 t( ẋ(t) + β∇f λ(t) (x(t))) 2 + (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) dt ≤ +∞ t0 t ẋ(t) 2 dt + +∞ t0 β 2 t ∇f λ(t) (x(t))) 2 dt + +∞ t0 (t -β(b + 2 -α))(f λ(t) (x(t)) -min H f )dt < +∞. t( ẋ(t) + β∇f λ(t) (x(t))) 2 + (t 2 -β(b + 2 -α)t)(f λ(t) (x(t)) -min H f ) = 0,

which shows that, as t → +∞ f λ(t) (x(t)) -min

When f is not convex, the convergence of the trajectories generated by (AVD) α is a largely open question. Recent progress has been made in[START_REF] Bot | A second-order dynamical approach with variable damping to nonconvex smooth minimization[END_REF] where the convergence of the trajectories of a system, which can be considered as a perturbation of (AVD) α , has been obtained in a non-convex setting.

A : H → H is λ-cocoercive (λ > 0) if for all x, y ∈ H Ay -Ax, y -x ≥ λ Ay -Ax 2 .

A Auxiliary results

In the proof of Theorem 2, we use the following straightforward result.

is satisfied for all X, Y ∈ H, if and only if C 2 -AB ≤ 0 and A, B ≥ 0.