F Leboeuf
email: leboeuf@lms.univ-poitiers.fr

G Bessonnet
email: bessonnet@lms.univ-poitiers.fr

P Seguin

P Lacouture

Energetic Versus Sthenic Optimality Criteria for Gymnastic Movement Synthesis

Keywords: dynamic synthesis, parametric optimization, minimum energy cost, minimum effort cost, handstand exercise

Dynamic synthesis of human movements raises the question of the selection of a suitable performance criterion able to generate proper dynamic behaviors. Two quite different criteria are likely to be appropriate candidates: the minimum effort cost (or sthenic criterion) and the minimum energy cost. The paper is aimed at clarifying the dynamic effects of these two fundamental criteria when considering movements executed with liveliness as they are in gymnastic. It is well known that the former cost generates movements with smooth dynamics. A special attention is devoted to the latter. The optimal control theory shows that minimizing the energy consumption results in actuating inputs of bang-off-bang type producing momentum impulses. When achieving dynamic synthesis, this criterion makes necessary to account for bounds set on driving torques. Moreover, when dealing with one-sided contacts, as in floor handstands, the unilaterality of contact forces must be explicitly accounted for since it tends to be infringed by impulsive efforts.

Numerical simulations of these formal properties are carried out using a parametric optimization technique, and considering the raising phase of floor handstands. It is shown that the energetic criterion tends to generate movements which exhibit similarities with their real counterpart performed by an expert gymnast. Conversely, the sthenic criterion produces movements quite different. But, a salient fact is that these ones proved to be easier to perform by young beginners. As a result, they could help to coach novice gymnasts.

Introduction

Dynamic analysis of human movement has been extensively developed during the last three decades. The central problem is the determination of internal forces, especially the driving 2 joint torques which control and give momentum to the motion. Two main approaches are at stake: inverse dynamics and optimal dynamic synthesis.

Due to the fact that internal efforts can not be directly assessed, inverse dynamics is an essential approach to be used to deal with human movement analysis. It needs complete information on motion kinematics which is the key point to be dealt with in this case. A video recording system provides, according to an acquisition frequency, a sequence of positions of anatomic markers. Using data filtering together with numerical differentiation techniques allows velocities and accelerations to be derived from recorded positions. It should be noted that both measurement and data processing hold uncertainties. The first suffers mainly from its lack of accuracy for identifying relative positions of joint rotation axes and centers of rotation. The second is very likely to produce discrepancies between real velocities and accelerations and their computed counterparts of which maximum values generally truncate the real ones. Consequently, the evaluated kinematic parameters could provide, through an inverse dynamics model, joint actuating torques with values fairly different from the real ones. Nevertheless, the inverse dynamics approach generally gives valuable results representing a basic assessment of internal efforts which were at work to generate the recorded movement. It can help to understand the joint motion coordination as during normal and pathological gait [START_REF] Alkjaer | Comparison of inverse dynamics calculated by two and three-dimensional models during walking[END_REF][START_REF] Selles | Lower leg inertial properties in transtibial amputees and control subjects and their influence on the swing phase during gate[END_REF][START_REF] Sanderson | Lower extremity kinematic and kinetic adaptations in unilateral below-knee amputees during walking[END_REF]. It is helpful also to analyze ergonomic conditions, especially while performing load lifting [START_REF] Plamondon | Validation of two 3D segments models to calculate the net reaction forces and moments at the L5S1 joint in lifting[END_REF][START_REF] Desjardins | Sensitivity analysis of segment models to estimate the net reaction force at the L5S1 joint in lifting[END_REF]. Another objective may consist in identifying some performance parameters when considering athletic exercises [START_REF] Arampatzis | A mathematical high bar-human body model for analysing and interpreting mechanical energetic processes on the high bar[END_REF][START_REF] Johnson | Muscle patterns in the mid-acceleration phase of sprinting[END_REF][START_REF] Purkiss | Methods for calculating internal mechanical work : comparison using elite runners[END_REF].

Unlike inverse dynamics, optimal dynamic synthesis needs either few experimental data, or not at all. However, it requires a great amount of computational effort. The movement to be dealt with is generated using an optimization technique which minimizes some performance criterion. The method consists in extracting an optimal solution from the equations governing the motion dynamics, while satisfying typical constraints the movement must comply with. This approach is not dependent on experimental uncertainties, and it correlates at best kinematics with forces that create the movement. Consequently, it could be an alternative to the inverse dynamics method. However, the minimizing criterion could generate, according to its nature, motions noticeably different from related human movements. The problem is that we do not know really what dynamic criterion governs human motricity. Thus, the choice of a performance criterion is simply a working basis. Its relevance will help to generate movements with kinematics and dynamics having insightful similarities with their real counterparts.

In this way, a variety of motion-synthesis attempts were carried out in order to get better knowledge of human movement dynamics. Two general methods have been used: optimal control theory and parametric optimization. The first method was implemented in [START_REF] Chow | Studies of Human Locomotion via Optimal Programming[END_REF][START_REF] Bessonnet | Optimal Motion Synthesis -Dynamic Modelling and Numerical Solving Aspects[END_REF][START_REF] Bessonnet | Optimal Gait Synthesis of a Seven-Link Planar Biped[END_REF] using the Pontryagin Maximum Principle (PMP) to study human locomotion [START_REF] Chow | Studies of Human Locomotion via Optimal Programming[END_REF] and for generating human-like gait of a planar biped [START_REF] Bessonnet | Optimal Motion Synthesis -Dynamic Modelling and Numerical Solving Aspects[END_REF][START_REF] Bessonnet | Optimal Gait Synthesis of a Seven-Link Planar Biped[END_REF]. The PMP is a powerful mathematical tool which 3 gives formally exact optimal solutions. However, its implementation may present some intricacy when dealing with one-sided contacts and considering multibody systems having a great number of degrees of freedom. This double situation is encountered when taking into consideration complex movements such as complete gait cycles and various gymnastic exercises. Parametric optimization makes it possible to cope with such dynamic and kinematic aspects of movements to be generated. This approach was developed on the basis of quite different intersegmental actuating models and related cost functions.

Authors like Pandy et al. [START_REF] Pandy | A Parameter Optimization Approach for the Optimal Control of Large-Scale Musculoskeletal Systems[END_REF][START_REF] Pandy | Optimal Control of Non-Ballistic Muscular Movements: A Constraint-Based Performance Criterion for Rising from a Chair[END_REF], Anderson and Pandy [START_REF] Anderson | Dynamic Optimization of Human Walking[END_REF], and Eberhard et al. [START_REF] Eberhard | Investigations for the Dynamical Analysis of Human Motion[END_REF] considered musculoskeletal models accounting for distinctive mechanical and physiological characteristics of each muscle, to generate human jumping [START_REF] Pandy | A Parameter Optimization Approach for the Optimal Control of Large-Scale Musculoskeletal Systems[END_REF][START_REF] Eberhard | Investigations for the Dynamical Analysis of Human Motion[END_REF], sit-to-stand movements [START_REF] Pandy | Optimal Control of Non-Ballistic Muscular Movements: A Constraint-Based Performance Criterion for Rising from a Chair[END_REF] and walking cycles [START_REF] Anderson | Dynamic Optimization of Human Walking[END_REF]. Various performance criteria were taken into account according to the movement considered: minimal time integral of square muscle forces [START_REF] Pandy | A Parameter Optimization Approach for the Optimal Control of Large-Scale Musculoskeletal Systems[END_REF][START_REF] Pandy | Optimal Control of Non-Ballistic Muscular Movements: A Constraint-Based Performance Criterion for Rising from a Chair[END_REF] and their time derivatives [START_REF] Pandy | Optimal Control of Non-Ballistic Muscular Movements: A Constraint-Based Performance Criterion for Rising from a Chair[END_REF], minimal task-time [START_REF] Pandy | A Parameter Optimization Approach for the Optimal Control of Large-Scale Musculoskeletal Systems[END_REF], minimal metabolic energy expended per unit distance covered [START_REF] Anderson | Dynamic Optimization of Human Walking[END_REF]. Movement synthesis based on musculoskeletal models is aimed at revealing neuromuscular coordination. Yet, this approach suffers from high complexity modeling together with uncertainties about the mechanical and physiological muscle behavior. Although this method deals with essential aspects of human movement organization, it is not very well adapted for generating wide-range movements with strong dynamic effects.

If one wants to optimize efficiently the global dynamics of complex movements, a simple joint intersegmental actuating model is much easier to handle than its musculoskeletal counterpart. This approach was particularly used for generating gait cycles designed for bipedal walking [START_REF] Beletskii | Parametric optimization in the problem of biped locomotion, Izs[END_REF][START_REF] Channon | Derivation of optimal walking motions for a bipedal walking robot[END_REF][START_REF] Chevallereau | Optimal reference trajectories for walking and running of a biped robot[END_REF][START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF][START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF] as well as for quadrupedal locomotion [START_REF] Muraro | Optimal Trajectories for a Quadruped Robot with Trot, Amble and Curvet Gaits for Two Energetic Criteria[END_REF]. Kuželički et al. [START_REF] Kuželički | Synthesis of standing-up trajectories using dynamic optimization[END_REF] achieved the synthesis of sit-to-stand movements for trans-femoral amputees with legprosthesis and for normal persons using a dynamic model comprising 11 degrees of freedom. The cost functional to be minimized consisted of two main terms representing the integrals of both joint actuating torques and their time derivatives. As shown in [START_REF] Bessonnet | Optimal Motion Synthesis -Dynamic Modelling and Numerical Solving Aspects[END_REF], the latter term has the effect of smoothing the driving torque variations, which is a necessary precaution for persons fitted with prosthesis. The performance criterion generally used for generating gait cycles is the minimum-effort cost. Since during gait the body segments are essentially submitted to gravity, the minimization of joint driving torques tends to generate upright walking patterns which require moderate actuating efforts to counterbalance the gravity effects [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF][START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF]. However, Chevallereau and Aoustin [START_REF] Chevallereau | Optimal reference trajectories for walking and running of a biped robot[END_REF], and Muraro et al. [START_REF] Muraro | Optimal Trajectories for a Quadruped Robot with Trot, Amble and Curvet Gaits for Two Energetic Criteria[END_REF] made also use of an energy cost to generate walking and running cycles for simplified bipeds. In fact, in section 3, we will recall that such a criterion generates bang-off-bang actuating inputs, which raises some questions regarding its application to human movement synthesis performed through a parametric optimization technique. 4 Indeed, the present paper is specifically aimed at showing that the minimum-effort cost and the minimum-energy cost yield movements exhibiting quite different kinematic features and dynamic characteristics. Using an accurate state parameterization, we will also emphasize that dealing with an optimal control problem using a parametric optimization technique produces suboptimal solutions, owing to the fact that unknown functions defined along the motion time are characterized by a finite number of optimized discrete parameters. Moreover, approximating state functions using high order polynomials as in [START_REF] Beletskii | Parametric optimization in the problem of biped locomotion, Izs[END_REF][START_REF] Channon | Derivation of optimal walking motions for a bipedal walking robot[END_REF][START_REF] Chevallereau | Optimal reference trajectories for walking and running of a biped robot[END_REF][START_REF] Muraro | Optimal Trajectories for a Quadruped Robot with Trot, Amble and Curvet Gaits for Two Energetic Criteria[END_REF], or splines of [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF][START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF], imply the continuity of joint accelerations. As a result, optimized joint driving torques will be continuous as well, whatever the minimized criterion may be.

class n C , 2  n
Especially, instead of bang-off-bang optimal controls, the energy criterion will produce continuous actuating inputs, even though they are not subjected to numerical bounds, which has no mathematical sense from the optimal control theory viewpoint [START_REF] Lewis | Optimal Control[END_REF]. This remark makes necessary the introduction of explicit limits for the driving torques in order to find relevant suboptimal solutions. An outline of the next sections is as follows. In the section 2, the content of the optimaldynamics problem we intend to deal with is detailed. We explicitly show in section 3 how the energetic cost we want to implement generates bang-off-bang actuating torques. Next, in section 4, the original dynamic optimization problem is converted into a minimization problem of mathematical programming using a state-parameterization method. Numerical simulations are presented in section 5. They were carried out considering a basic gymnastic movement, the handstand, performed as a floor exercise.

Stating a dynamic optimization problem

Optimal motion synthesis is based on stating and solving a dynamic optimization problem. This can be dealt with as an optimal control problem, or converted into a parametric optimization problem. We will briefly consider the former problem in Section 3 in order to gain an insight into the formal nature of optimal control inputs which result from the criterion minimized, and the latter will be developed in Section 4 to compute suboptimal solutions using standard computing codes implementing quadratic sequential programming algorithms.

As regards its kinematics and dynamics, the human body is commonly described as a multibody system made of rigid segments and actuated rotational joints. In Section 5, a gymnast is depicted as a rooted planar system with an open tree-like topology (Fig. 3).

Dynamics equations and boundary conditions

The dynamics of movements we consider over an interval of time ] , 0 [ T is governed by a set of equations summarized as the single n-vector relationship

) ( )) ( , ) ( , ) ( ( ], , 0 [ t τ t q t q t q B T t      , (1)
where 5 -

T n q q q ) ,..., ( 1 
is the vector of joint coordinates i q , n i  , n being the number of degrees of freedom.

-q  and q  are the first and second order time-derivatives of q.

- T n ) ,..., ( 1    
is the vector of joint driving torques, assuming that all degrees of freedom are actuated.

The phase trajectory

)) ( ), ( ( t q t q t   must comply with some initial and final conditions we formally write under the form

             ) 2 , ( 0 )) ( ), ( ( ) 2 , ( 0
)) 0 ( ), 0 ( ( 0 0 0 n n T q T q n n q q T n T n T     . ( 2
)

State constraints

Generally, joint motions have limited range and must avoid hyperextension, which results in formulating inequality constraints such as

              0 ) ( : )) ( ( 0 ) ( : )) ( ( , ], , 0 [ max min k k k n k k k q t q t q g t q q t q g n k T t . (3)
Let us notice that the joint velocities could be submitted to similar limitations. We formally combine inequalities (3) into the vector function

T n g g g g ) ,..., ( 1  , 0 )) ( ( ], , 0 [   t q g T t , (4)
where g n stands for the total number of constraints. It should be noted that min k q and max k q represent given data (see Table 1 in section 5.2).

Sthenic constraints

In addition, joint driving torques must be bounded too, especially when minimizing an energetic cost (see section 3). Consequently, there is the need for taking into account limiting constraints as

              0 ) ( : )) ( ( 0 ) ( : )) ( ( , ], , 0 [ max min k k k n k k k t t h t t h n k T t       , ( 5
)
where the lower bounds min k  are assumed to be negative. Generally, such limits are not symmetric:

max min k k     . 6
Moreover, interaction forces such as ground reactions due to one-sided contacts must fulfill unilaterality requirement and non sliding conditions (see section 5). Constraints of this type depend explicitly on phase variables and actuating torques. Thus, they can be written as 0

)) ( ), ( ), ( ( , ], , 0 [ 2     t t q t q h n j T t j n r   , (6)
r n being the number of constraints.

We will term sthenic constraints the conditions ( 5) and ( 6) that we summarize using the vector function

T n h h h h ) ,..., ( 1  ; 0 )) ( ), ( ), ( ( ], , 0 [   t t q t q h T t   . (7)
The subscript h n represents the number of scalar constraints.

Performance criteria

We consider the following two criteria -A sthenic cost defined as the time integral of quadratic driving torques

  T s T s dt t D t J 0 2 1 ) ( ) (   , ) ,..., ( diag 1 n s D    , (8)
-An energetic cost which is the time integral of absolute joint actuating powers

    T n i i i i e dt t t q J 0 1 ) ( ) (    , ( 9
)
where i  and i  are weighting coefficients.

In fact, we will take into consideration a combination of these two costs, that is to say the

mixed criterion ] 1 , 0 [ , ) 1 (        e s J J J ( 10
)
which gives s J and e J when  takes the values 0 and 1 respectively.

Finally, the dynamic optimization problem can be summarized as follows: a value of α being given, find a joint trajectory ) (t q t  and a vector-function

) (t t  
solution of (1), minimizing the criterion (10) while satisfying the boundary conditions ( 2), together with the double set of constraints ( 4) and ( 7). This problem is dealt with in the following two sections.

Optimal actuating inputs

The analytic nature of optimal driving torques minimizing the criterion (10) can be revealed using the Pontryagin Maximum Principle [START_REF] Lewis | Optimal Control[END_REF] without having to solve the optimal control problem. We refer the reader to [START_REF] Bessonnet | Optimal Motion Synthesis -Dynamic Modelling and Numerical Solving Aspects[END_REF][START_REF] Bessonnet | Optimal dynamics of actuated kinematic chains -Part 2: Problem statements and computational aspects[END_REF] for a general approach developed to deal with dynamic optimization problems related to multibody systems. It is shown that the necessary conditions for optimality are easier to derive from a Hamiltonian dynamic model than from its Lagrangian counterpart. Thus, let us briefly introduce:

7 -The Lagrangian: ) ( ) , ( ) , ( q V q q T q q L    
, T being the kinetic energy and V the gravity potential of the mechanical system.

-The conjugate momentum:

q L p     / .
-The Hamiltonian vector phase variable:

T T T p q x ) , (  .
Then, it can be easily shown [START_REF] Bessonnet | Optimal Motion Synthesis -Dynamic Modelling and Numerical Solving Aspects[END_REF] that the Hamilton equations equivalent to (1) take the simple form of the 2n-order vector-state equation

 A x F x   ) (  , (11)
where

           n n n n I A 0 .
Moreover, through ( 8) and ( 9), the criterion J can be written under the standard form

  T dt t t x l J 0 )) ( ), ( (  . ( 12
)
Then, defining the Pontryagin function [START_REF] Lewis | Optimal Control[END_REF] )

, ( ) ) ( ( ) , , ( , 2    x l A x F y y x H y T n      , (13)
and putting aside momentarily the state constraints (4), the Pontryagin Maximum Principle states that for any solution ) , (  x of ( 11) satisfying the constraints ( 2) and ( 7), and minimizing [START_REF] Pandy | A Parameter Optimization Approach for the Optimal Control of Large-Scale Musculoskeletal Systems[END_REF], there exists a 2n-vector function ) (t y t  satisfying the so-called adjoint system

T x H y T t ) / ( ], , 0 [       , (14)
and the maximality condition for H

) ), ( ), ( ( )) ( ), ( ), ( ( ], , 0 [ u t y t x H t t y t x H T t Max U u    , ( 15
)
where U stands for the set of feasible control variables s i  defined by the sthenic constraints ( 5) and ( 6). If only the constraints ( 5) are taken into account, the set U appears as a parallelepiped in the n-dimensional control space. Constraints (6) eliminate half-spaces (see [START_REF] Bessonnet | Optimal dynamics of actuated kinematic chains -Part 2: Problem statements and computational aspects[END_REF]) which truncate U and reduce it to a polyhedron. We need to examine the maximality condition [START_REF] Eberhard | Investigations for the Dynamical Analysis of Human Motion[END_REF] when the only constraints ( 5) are accounted for. Knowing the expression of the Lagrangian l in ( 12) and ( 13) through ( 8) and ( 9), the condition (15) allows the optimal inputs s i  to be explicitly formulated as functions of x and y. When J is reduced to s J , it is well known that optimal inputs appear as saturating functions which are continuous functions truncated by upper and lower bounds of the s i  (see [START_REF] Lewis | Optimal Control[END_REF]). At this point, we want to be more specific about the optimal inputs which result from the purely energetic cost we intend to implement as the limiting case 0   in [START_REF] Bessonnet | Optimal Motion Synthesis -Dynamic Modelling and Numerical Solving Aspects[END_REF].

First, from [START_REF] Bessonnet | Optimal Gait Synthesis of a Seven-Link Planar Biped[END_REF], we have )

(x F x q i i i     , n i  .
Then, the Lagrangian of e J in (9) can be written as

8    n i i i i e x F x l 1 ) ( ) , (    .
Accordingly, the Pontryagin function takes the form

) ( ) , ( ) , , ( 1 i n i i i i i n u F u y y x C u y x H        , in which ) , ( y x C
stands for terms independent of u.

This relationship shows that maximizing H amounts to maximizing the sum in the right hand member. This maximization problem can be solved explicitly when the only active sthenic constraints are reduced to (5), i.e. when the set U is a parallelepiped. In that case, it is possible to write

) ( ) ( 1 1 max min i n i i i i i n u i n i i i i i n U u u F u y Max u F u y Max i i i                 , ( 16
)
because the s u i vary in U independently of each other.

Next, the problem to be solved consists of finding a real variable i u solution of the elementary maximization problem

) ( 0 , max min i i i i u u u Max i i i i i             . ( 17
)
As shown in the appendix, the solution can be expressed using a dead zone function and the sign function. Since i u is determined for every fixed time t, we have

i i u t  ) (  .
Consequently, from the relationships (A2) and (A3), we get

                     0 )) ( ( if , )) ( ( sign 0 )) ( ( )) ( ( if , , , )) ( ( ) ( dez ) ( max max min t x F t y t q t x F t x F t y t i i n i i i i i i i i n i       . ( 18
)
The figure 1 shows typical variations of an optimal actuating torque resulting from the energy cost minimization. It must be emphasized that whatever the limiting values min i  and max i  may be, the corresponding optimal input will be of bang-off-bang type as shown in Figure 1.

If limiting values are not explicitly taken into account, singular solutions will result, taking the zero value over the entire time interval except for discrete points where the s i  will generate infinite impulses. However, optimization carried out using state-parameterization techniques has the effect of smoothing the suboptimal solutions which result. Thus, suboptimal actuating inputs are smoothed as well. Nevertheless, if they are not subjected to predefined bounds, they will be likely to reach disproportionate values over brief timesubintervals. This situation must be avoided. 9 Fig. 1. Two examples of bang-off-bang optimal input i  minimizing the energy cost, according to the values of lower and upper bounds to be respected.

To finish, let us add that when the constraints ( 6) are accounted for, optimal inputs are no longer of purely bang-off-bang type since the set U is a truncated parallelepiped. However, the problem is not fundamentally changed, because as soon as constraints in [START_REF] Arampatzis | A mathematical high bar-human body model for analysing and interpreting mechanical energetic processes on the high bar[END_REF] will not be active, the bang-off-bang nature of optimal inputs will reappear.

Parametric optimization method

Parameterizing an optimization problem consists first in approximating some variables to be optimized using a finite set of discrete parameters and, second, in recasting the criterion to be minimized into a simple function of these parameters. The reader is referred to [START_REF] Hull | Conversion of Optimal Control Problems into Parameter Optimization Problems[END_REF] for a short presentation of general ideas commonly used to convert optimal control problems into parameter optimization problems. Practical methods for solving such problems could be found in [START_REF] Betts | Practical Methods for Optimal Control Using Nonlinear Programming[END_REF]. Considering a dynamic optimization problem, the parameterization may be carried out according to quite different approaches depending on whether state variables or control variables are parameterized using approximation functions. The latter method was used notably in [START_REF] Pandy | A Parameter Optimization Approach for the Optimal Control of Large-Scale Musculoskeletal Systems[END_REF][START_REF] Pandy | Optimal Control of Non-Ballistic Muscular Movements: A Constraint-Based Performance Criterion for Rising from a Chair[END_REF][START_REF] Anderson | Dynamic Optimization of Human Walking[END_REF] where the values of control variables defined at a finite number of points along the motion time are taken into account as optimization variables. The control history is then reconstructed by linear interpolation between nodal points. Then, an initial state being given, a forward integration of motion equations yields a generalized final state together with the value of the cost function with respect to the optimization variables. The parameterization of the original problem is thus completed. This technique is appealing because it eliminates untimely oscillations of state variables between knots. However, its accuracy strongly depends on a given initial state together with the precision of the forward integration technique which could be very time consuming. It does not seem suitable for achieving the optimal synthesis of movements of which end states need to be optimized as required, for example, for a gait cycle [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF][START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF].

) (t i  ) , , ) ( ( dez ) ( max min i i i i t t      t max i  min i  0 min 2 i  max 2 i  ) 2 , 2 , ) ( ( dez ) ( max min i i i i t t      10
Parameterizing the state variables needs to use smooth approximation functions, at least twice differentiable in order to ensure the existence of accelerations. This preliminary operation being done, actuating inputs can be easily computed by performing inverse dynamics. The conversion of the original dynamic optimization problem into parametric optimization is straightforwardly achieved by the computation of the cost function with respect to the new optimization variables.

A number of approaches have been used to approximate the state variables which can be chosen in the task space to describe, for instance, swing foot and hip trajectories for generating gait steps [START_REF] Beletskii | Parametric optimization in the problem of biped locomotion, Izs[END_REF][START_REF] Channon | Derivation of optimal walking motions for a bipedal walking robot[END_REF] or, more generally, in the joint space in order to be in command of the joint motion coordination. In either cases, the most frequently used approximation functions are polynomials [START_REF] Beletskii | Parametric optimization in the problem of biped locomotion, Izs[END_REF][START_REF] Channon | Derivation of optimal walking motions for a bipedal walking robot[END_REF][START_REF] Chevallereau | Optimal reference trajectories for walking and running of a biped robot[END_REF][START_REF] Muraro | Optimal Trajectories for a Quadruped Robot with Trot, Amble and Curvet Gaits for Two Energetic Criteria[END_REF] and splines [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF][START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF]. The former are easy to use but are prone to undesirable oscillations, especially if they have high order. This could result both in a limited minimizing effect of the cost function, and in sizeable constraint infringements. Spline functions represent a more accurate choice, either under the form of B-splines used as base functions [START_REF] Martin | Minimum effort motions for open chain manipulators with task-dependent end-effector constraints[END_REF][START_REF] Lo | Human motion planning based on recursive dynamics and optimal control techniques[END_REF], or low order polynomials successively linked to each other at connecting points (or knots) defined along the motion time [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF][START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF]. In the latter case, polynomial functions of order three are sufficient to ensure the existence of second order derivatives (i.e. accelerations) at connecting points. This way, one gets spline functions of class C 2 over the traveling time. However, the non differentiability of accelerations at knots results in jerky variations which can be transferred to actuating torques. A more satisfactory approach developed in [START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF] consists of using splines of class C 3 obtained by connecting four order polynomials up to their third derivatives. Rough variations of actuating torques at knots are thus avoided. We have implemented this more accurate technique in order to generate smooth movements on the overall motion time.

Defining a set of approximation functions

We adapt and summarize the method developed in [START_REF] Seguin | Development of a parametric optimization technique for dynamic synthesis of smooth movements; application to gait synthesis[END_REF]. First, a set of connecting times must be defined over the interval [0, T] as

  ) ( ,..., ),..., 0 ( 1 1 T t t t N j    , N T t t j j / 1    ,
where N is the number of subintervals of equal length.

Second, over each subinterval ] , [ :

1   j j j t t I
, and for every joint coordinate i q , a four order polynomial ij P is defined over the

normalized interval ] 1 , 0 [ such that                    4 4 3 3 2 2 1 0 1 ) , ( ) ( , , ,       ij ij ij ij ij ij ij i j j j j c c c c c C P t q t t t t I t N j , ( 19
)
where ij C is the 5-order vector of coefficients:

T ij ij ij ij ij ij c c c c c C ) , , , , ( 4 3 2 1 0  . 11
Next, for every subscript i, the polynomials ij P are successively linked to each other, up to their third derivatives, as follows

                   ) 0 ( ) 1 ( ) 0 ( ) 1 ( ) 0 ( ) 1 ( ) 0 ( ) 1 ( , 2 1 , 1 , 1 , 1 , ij j i ij j i ij j i ij j i P P P P P P P P N j             . ( 20
)
Furthermore, we set

           ) ( ) 1 ( ) ( ) 1 ( ) 0 ( ) 0 ( ) 0 ( ) 0 ( 1 1 T q P T q P q P q P i iN i iN i i i i     , N j   2 , ) ( ) 0 ( j i ij t q P  . ( 21
)
Relationships in [START_REF] Muraro | Optimal Trajectories for a Quadruped Robot with Trot, Amble and Curvet Gaits for Two Energetic Criteria[END_REF] mean that the end states )) 0 ( ), 0 ( (

i i q q  and )) ( ), ( ( T q T q i i



, together with inner values of the s q i at connecting points will be considered as free parameters to be optimized, which we formally assemble in the set of vectors

T i i j i i i i N i i T q T q t q q q x x X n i )) ( ), ( ),..., ( ),..., 0 ( ), 0 ( ( : ) ,..., ( , 3 1 
     . (22)
Now, using the relationships ( 20) and ( 21), the coefficients in ( 19) can be computed as linear functions of the above parameters. Actually, for any subscript i, there are 5N coefficients introduced in [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF]. They are subjected, through [START_REF] Bessonnet | Parametric Optimization Approach to Walking Pattern Synthesis[END_REF] and ( 21), to 5N-1 linear equations.

Reducing to three the order of one of the polynomials ij P , for instance the last, iN P , by setting

0 4  iN c
, the set of equations ( 20) and ( 21) can be solved in the 5N-1 remaining coefficients that we bring together in the vector

T T iN T i i C C C ) ,..., ( 1 
. Thus, this one appears as a function of the vector i X defined in [START_REF] Kuželički | Synthesis of standing-up trajectories using dynamic optimization[END_REF], that is

) ( i i X C   . ( 23
)
It should be noted that the function  is formally the same for every index i. Accounting for [START_REF] Lewis | Optimal Control[END_REF] in [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF], each polynomial ij P can be recast as a function ij

 such as ) , ( ) , (    i ij ij ij X C P  .
Next, defining a function i  over the time interval [0, T] by setting

) , ( ) , ( , ,    i ij i i j X t X I t N j     , the generalized coordinate i q is approximated by i  : ] , 0 [ T t  , ) , ( ) ( t X t q i i i   . ( 24
)
Moreover, if we set

12 T T n T X X X ) ) ( ,..., ) (( 1  , T n n t X t X t X )) , ( ),..., , ( ( ) , ( 1 1
    , the configuration vector q is approximated by the function  as

) , ( ) ( ], , 0 [ t X t q T t    . (25)

Stating a parametric optimization problem

Through (25), the vector of actuating torques  in (1) can be written as a vector   depending on X and t as follows

) , ( ) ( )) , ( ), , ( ), , ( ( : ) , ( 2 2 t X t t X t t X t t X B t X               .
The criteria s J and e J in ( 8) and ( 9) are then approximated by the functions of

X      T s T s s dt t X D t X X F J 0 2 1 ) , ( ) , ( ) (   , ( 26
)       T e T e e dt t t X D t X X F J 0 ) , ( ) , ( ) (   . ( 27
)
The mixed criterion (10) becomes

) ( ) 1 ( ) ( ) ( X F X F X F J e s       . (28)
The boundary constraints (2) can be converted as

0 ) ) , ( ), , ( (
0 ) ) 0 , ( ), 0 , ( ( 0       t T X T X t X X T       , ( 29
)
where it should be seen that, through [START_REF] Hull | Conversion of Optimal Control Problems into Parameter Optimization Problems[END_REF] and ( 22)

T n x x t X ) ,..., ( ) 0 , ( 1 1 1    , T n x x X ) ,..., ( ) 0 , ( 2 1 2   , T n N N x x T X ) ,..., ( ) , ( 2 1 2     , T n N N x x t T X ) ,..., ( ) , ( 3 1 3      .
Now, the constraint functions g and h in ( 4) and ( 7) can be recast as the following functions of X and t

)) , ( ), , ( ( ) , ( , t X t X g t X g t     , )) , ( ), , ( ), , ( ( ) , ( , t X t X t X h t X h t       .
Finally, the distributed constraints ( 4) and ( 7) will be accounted for as the discrete conditions set at connecting times 13

           0 ) , ( 0 ) , ( , 1 k k t X h t X g N k . ( 30
)
To conclude, the original dynamic optimization problem defined in ( 10), ( 2), ( 4) and ( 7) is converted into a standard minimization problem which consists of finding a discrete set X of optimization variables minimizing the cost function ( 28) (( 26) or [START_REF] Martin | Minimum effort motions for open chain manipulators with task-dependent end-effector constraints[END_REF] as well) while satisfying equality constraints [START_REF] Seguin | Development of a parametric optimization technique for dynamic synthesis of smooth movements; application to gait synthesis[END_REF] and inequality constraints [START_REF] Clement | Contribution of eye positioning to control of the upside down standing posture[END_REF]. This constrained non linear problem of mathematical programming can be solved using existing computing codes implementing sequential quadratic programming algorithms which have proved to be quite efficient.

Dynamic synthesis of a gymnastic movement: the handstand

The handstand may be executed at all apparatus such as rings or high bar. As a gymnastic floor exercise (Fig. 2), it can be divided into three phases. First, the gymnast puts his hands in floor support. Second, after feet takeoff, he reaches the inverse balance. Third, the gymnast must perform a stationary handstand. This movement is considered as a fundamental gymnastic sequence which is difficult to perform by young beginners as it requires a specific development of muscular coordination. For this reason, it held researchers' attention. For example, in [START_REF] Clement | Contribution of eye positioning to control of the upside down standing posture[END_REF][START_REF] Pozzo | Head and trunk movements in frontal plane during complex dynamic equilibrium tasks in humans[END_REF], the authors focused their analyses on kinematic parameters and reaction forces measured using a video recording system and a force plate, respectively. In [START_REF] Kerwin | Strategies for maintaining a handstand in the anteriorposterior direction[END_REF], the computation of joint actuating torques was carried out to complete the analysis of the motionless final posture. The authors used a simple dynamic model with three degrees of freedom, considering both legs and arms together. Through such analyses, one cannot know how the initial momentum of the raising phase and its segmental coordination influence the dynamic control of the handstand. In this section, our purpose is to achieve the dynamic synthesis of the raising phase of floor handstands in order to provide a new insight into the way the joint motions can be coordinated. An essential further objective of the paper is to show that minimizing a sthenic cost or an energetic cost as defined in section 3, generates movements with distinctive features that highlight the effects of each criterion and which have proved to be useful in coaching young gymnasts.

Kinematic modeling and state constraints

Through the raising phase of the handstand, the gymnast is modeled in his sagittal plane as a 9-link multibody system (Fig. 3). His hands are flat on the ground (link L 0 ), and both arms are assumed to perform the same movement (links L 1 and L 2 ). Trunk and head are considered as an only rigid link (L 3 ). As schematized in Fig. 3, the model has nine active joints.

  1 O 6 A 0 x  0 y   6 O
15 Fig. 3. Planar segmental model of a gymnast performing the raising phase of a handstand on horizontal ground.

The kinematic model is described by a set of nine generalized coordinates which are the joint angles of relative rotations between adjacent links. Thus, we introduce the configuration vector T n i q q q q ) ,..., ,..., ( 1  with 9  n and

i i q  
as defined in Fig. 3.

At initial time of the raising phase, the foot of the lower leg takes off. Only three equality constraints in [START_REF] Pozzo | Head and trunk movements in frontal plane during complex dynamic equilibrium tasks in humans[END_REF] specify the initial state of the gymnast: the first indicates that the tip 6 A of the foot is still in contact with the supporting ground (Fig. 2), while the next two ones mean that its initial velocity is zero.

                0 )) 0 ( ( ( )) ( ), 0 ( ( 0 )) 0 ( ( ( )) ( ), 0 ( ( 0 )) 0 ( ( )) ( ), 0 ( ( 0 6 3 , 0 0 6 2 , 0 0 6 1 1 , 0 y q A V T q q x q A V T q q y q A O T q q          . ( 31
)
As joint positions and velocities are not specified, we emphasize that the initial posture and the initial momentum of the simulated gymnast will be optimized. At final time, the handstand is achieved as prescribed by the complete set of constraints [START_REF] Dempster | Space Requirements of the Seated Operator[END_REF] where the s q iT stand for given values.

           0 ) ( )) ( ), ( ( ) ( )) ( ), ( ( , , , T q T q T q q T q T q T q n i i i n T iT i i T      . ( 32
) 7 4 O O  5 O 6 O 8 O 9 O 4  5  6  7  8  9  3 x  4 x  5 x  6 x  7 x  8 x  9 x  ) ( 4 L ) ( 5 L ) ( 6 L ) ( 7 L ) ( 8 L ) ( 9 L  6 A 2 O 3 O 1  2  3  0 x  1 x  2 x  ) ( 0 L ) ( 1 L ) ( 3 L 1 O 0 y  0 A          1 C 2 C 3 C 4 C 5 C 8 C
MUBO 05026

Final version [START_REF] Beletskii | Parametric optimization in the problem of biped locomotion, Izs[END_REF] The state constraints (4) are reduced to the limitation of joint motion ranges at knees and ankles, that is q q q q q q q q . (33)

          .) (deg ) 60 , 20 ( ) , ( ) , ( .) (deg

Dynamic modeling and sthenic constraints

We used the Lagrange equations which give to (1) the formal structure

    ) ( ) , ( ) ( q G q q C q q M    , ( 34
)
where M is the ) ( n n  -inertia-matrix, C is the n-vector of Coriolis and centrifugal terms, and G stands for the gravity terms.

Bounds set on actuating torques in ( 5) are given in Table 1. Indicated values are simply based on preliminary optimization tests carried out using a small number of nodal points, exactly five. Table 1. Limiting values of actuating torques.

Contact forces exerted by the ground on the hands can be represented by a system of three forces, two of them being vertical forces applied to the hand extremities, the third representing the tangential friction component acting along the axis ) ; ( 0 1 x O  (Fig. 4).

Unilaterality of contact and non sliding condition can be expressed as constraints of type [START_REF] Arampatzis | A mathematical high bar-human body model for analysing and interpreting mechanical energetic processes on the high bar[END_REF], that is:

0 ) ( ) , , ( 0 ) , , ( 0 ) , , ( 0 1 0
1 3 2 2 2 1 2               A O n A n O n N N f T q q h N q q h N q q h       , ( 35
)
where f is a dry friction coefficient. These interaction forces were computed with respect to ) , , (  q q  using Newton-Euler equations formulated for the whole system.

 0 A 1 O 0 1 y N O  0 0 y N A  0 x T   ) ( 0 L ) ( 1 L
Fig. 4. System of forces exerted by the ground on the "hands" ) ( 0 L .

Numerical results

We used the biometric data, given in table A1 of the appendix, of an expert gymnast.

A first group of simulations was carried out without accounting for the state constraints [START_REF] Athans | Optimal Control -An Introduction to the Theory and Its Applications[END_REF]. Therefore, the legs are free to flex without limitations. Moreover, the relative position of the lower foot tip is not specified. Consequently, the initial hand-foot distance will be optimized. Results of nine simulations are shown in Table 2. They were obtained by varying two main parameters:

-The number N of control points in [START_REF] Saidouni | Generating globally optimised sagittal gait cycles of a biped robot[END_REF] which is a fundamental data for the optimization technique employed. The actual values are 5, 10 and 20 which generate 72, 117 and 207 optimization parameters, respectively.

-The value of the weighting factor  in the mixed criterion [START_REF] Lo | Human motion planning based on recursive dynamics and optimal control techniques[END_REF] which is an essential characteristic of the minimization problem to be solved.

Four values of  were taken into account. The first, equal to 1, defines a purely sthenic criterion. The second value, and particularly the third, represents a nearly energy cost. The fourth which defines a purely energetic criterion was used only for the second and third values of N which provide the better accuracy. Table 2. Minimal costs J, J s , J e , and optimal initial foot-hand distance with respect to the number N of control points and the values of the weighting factor .

On each of the first three lines in  , the fluctuations of results versus the number N of control points reflects the suboptimal nature of the minimization technique. On the columns, the variations are quite significant and express the change from the sthenic criterion to the energetic one. Results reveal that there are large increases of the total actuating effort, and at the same time the actuating work is greatly decreased. Fig. 5 shows three optimal movements computed for N equal to 20. The corresponding time charts of actuating torques and contact forces are displayed in Fig. 6. The change from the minimum effort criterion to the nearly minimum energy cost shows two different strategies for achieving the final balance of the handstand (Fig. 5): in the first . In the latter case, joint torques do not cease oscillating between the zero value and their maximum authorized values. This denotes a tendency for bang-off-bang dynamic behavior, as shown in section 3. Accordingly, contact forces are also oscillating. It should be noted that the unilaterality requirement is satisfied over the whole time interval. However, the bound constraint set on the lower leg hip-torque is slightly infringed twice between control points at beginning of the movement. 2). Unfortunately, the unilaterality condition is frequently infringed both at wrist 1 O and at tip 0 A of the hand (Fig. 7). Furthermore, the hip-torque limit of the lower leg is exceeded between knots. In brief, the minimum energy cost generates fast varying torques which are incompatible in the present case with one-sided contact conditions. Through the above examples, we can conclude that the energetic cost creates movements with little internal gesticulation and fast varying actuating torques. Quite the reverse, the sthenic cost generates more gesticulation while driving torques show small variations.

In two further simulations we introduce state constraints in order for the optimized motions to mimic the gymnast movement. Since the gymnast performs slight leg flexions, we limit the knee flexions to 30 degrees as given in [START_REF] Athans | Optimal Control -An Introduction to the Theory and Its Applications[END_REF]. Moreover, the initial hand-foot distance is no longer free but is fixed to 0.95m, the value observed in the real movement (Fig. 2).

Figure 8 shows that the kinematic differences between both optimized motions are attenuated by the new constraints. However, there is still more leg flexion generated by the sthenic criterion than its nearly energetic counterpart at beginning and end of movements. In the same way, the main differences shown previously in Fig. 6 reappear in Fig. 9: there are more variations for joint torques and contact forces in the energetic case than in the sthenic one. Furthermore, in the minimum energy case, the lower leg hip-torque limit is not respected between knots and the unilaterality condition is slightly infringed at the beginning of the motion.

To finish, we make a comparison between the dynamics of both the simulated movement and the movement performed by the gymnast. The torques generated by the gymnast were computed using inverse dynamic analysis. The kinematics was determined by means of a video recording system, data filtering and numerical differentiations techniques to compute joint velocities and accelerations. Fig. 10 shows that torques at hips and knees exhibit moderate extremal values and few variations. Arm driving torques are more fluctuating than the previous ones and reach sizeable values showing that the gymnast uses his arms very much to control his movement. The comparison between these results with their counterparts shown in Fig. 9 is not obvious.

Nevertheless, in the nearly energetic case ) 01 . 0 (   , one can see that torques at wrists and shoulders take important fast-varying values, especially at the beginning, which is fairly similar to the situation observed in Fig. 10. More striking similarities appear in the global results shown in Table 3: the actuating work and the initial kinetic energy are roughly of the same order in the simulated energetic case and in the real movement. But, in this respect, one can wonder why the mechanical energy expended in the latter case is less than in the former which represents an idealized situation. This is probably due to the fact that real actuating torques are undervalued through the numerical treatment undergone by the kinematic data. Indeed, data filtering and numerical differentiations yield velocities and even more accelerations of which peak values are attenuated. This certainly affects the computed dynamics which is found weaker than it is really.

Although the above comparisons are not quite clearly definite, they show that the dynamics of the physical movement has real similarities with the minimum-energy optimal motion. These results show also the need for improving both the accuracy of experimental and numerical techniques required to carry out dynamic analysis of human movements, and the precision of optimization methods which provide only local suboptimal solutions. Table 3. Minimal costs obtained with two additional constraints: foot-hand distance is fixed at 0.95m, and knee flexion is limited to 30°.

Concluding remarks

Optimal dynamic synthesis is a significant means to gain an insight into the way kinematic organization and dynamic coordination of human movements are linked together. Quite different results may be found according to the dynamic cost minimized. Indeed, the analysis of numerical simulations showed that the choice of a performance criterion is an ambiguous question because human movements are in the heart of a compromise between the least energy expended and the least effort put in. We paid a special attention to dealing with the energy cost. Especially, the effect of such a criterion was emphasized, and the problem of its implementation was explicitly clarified. As the minimum energy cost generates bang-off-bang actuating torques, we underline that it requires an accurate parametric optimization method to reveal the fast and complex variations 24 of suboptimal inputs which result. Moreover, it appeared that a purely energy cost could be not satisfactory for generating movements with one-sided contact conditions which might be infringed as during the raising phase of floor handstands.

A better approach is based on considering a mixed criterion associating actuating effort with driving work. This way of tackling dynamic movement synthesis offers a helpful adaptability of the optimization process. In numerical simulations, the progressive change from one cost to the other reveals that the movement is controlled in quite different ways, which results in noticeable modifications of its kinematics. Moreover, while the criterion with predominant energy generates movements with kinematic aspects similar to the real handstand performed by a gymnast, the movement produced by the sthenic criterion exhibits quite dissimilar features. Nevertheless, this movement characterized by smooth actuating torques proved to be easier to execute by novice gymnasts. Conversely, movements governed by a minimum energy cost are harder to perform probably for the reason that they need to control accurately successive brief impulses of actuating torques. -Case 2:

      Then 0 2  a , 0 1  a
, and the maximum of φ is reached for 0  u .

-Case 3:

   Then 0 2 1   a a
, and φ reaches its maximum for a u  . These results can be summarized using the dead zone function (see [START_REF] Lewis | Optimal Control[END_REF]) generalized here as the following asymmetrical function of the three variables

 Biometric data

Reference [START_REF] Dempster | Space Requirements of the Seated Operator[END_REF] was used to compute the biometric data of the expert gymnast.

Fig. 2 .

 2 Fig. 2. Stick diagram of a floor-handstand raising phase performed by an expert gymnast.

Fig. 5 .

 5 Fig. 5. Stick diagrams of optimal motions computed on the basis of three different criteria defined by the values of the weighting factor  .

Fig. 6 .

 6 Fig. 6. Time charts of actuating torques and contact forces for three values of the weighting factor ; vertical hatched lines are drawn at control points.

Fig. 7 .

 7 Fig. 7. Stick diagram of movement, and time charts of driving torques and contact forces generated by a purely energetic criterion ) 0 (  

Fig. 8 .

 8 Fig. 8. Stick diagram of constrained motions mimicking the gymnast movement.

Fig. 9 .

 9 Fig. 9. Time charts of driving torques and contact forces generated by the sthenic ) 1 (  

Fig. 10 .

 10 Fig. 10. Time charts of actuating torques exerted by the gymnast during the raising phase of the floor handstand.

27 -

 27

  Fig. A1. Variations of the function φ over the interval ] , [ b a

28 Fig. A2 .

 28A2 Fig. A2. Graph of an asymmetrical dead zone function dez(x,a,b)

Table A1 .

 A1 Gymnast biometric data.

	Body segment	Mass (kg) Length (m)	Local abscissa of CoG C i (m)	Moment of inertia at C i (kg.m 2 )
	Hand	0.408	0.093	0.046	0.005
	Fore-arm	3.808	0.201	0.114	0.016
	Arm	2.176	0.263	0.150	0.014
	Trunk-head	39.30	0.568	0.193	3.213
	Thigh	6.800	0.434	0.188	0.134
	Shin	3.162	0.372	0.161	0.040
	Foot	0.986	0.126	0.063	0.035

Appendix  Solving an elementary maximization problem

The presentation which follows generalizes to the cases b a   and 0   (see below) a result originally established in a somewhat different form in [START_REF] Athans | Optimal Control -An Introduction to the Theory and Its Applications[END_REF] and in [START_REF] Lewis | Optimal Control[END_REF].

In [START_REF] Channon | Derivation of optimal walking motions for a bipedal walking robot[END_REF], setting -

, solution of the maximization problem:

First, assuming that 0   , we define the function


Bastien Sérée

Loïg Jezequel

Didier Lime

An algorithm for single-source shortest paths enumeration in parameterized weighted graphs

Keywords: Shortest paths, weighted graphs, parameterized graphs

We consider weighted graphs with parameterized weights and we propose an algorithm that, given such a graph and a source node, builds a collection of trees, each one describing the shortest paths from the source to all the other nodes of the graph for a particular zone of the parameter space. Moreover, the union of these zones covers the full parameter space: given any valuation of the parameters, one of the trees gives the shortest paths from the source to all the other nodes of the graph when the weights are computed using this valuation.

Introduction

For many real-world systems and problems there are natural discrete-event abstractions, which can be modelled by a graph or a derived formalism. In many cases, we can also identify resources that need to be optimized (distance, memory, energy, time, etc.) and then (extensions of) weighted graphs are a formalism of choice.

One of the most basic problems is to find optimal paths, for which the accumulated weight (also often called cost) is minimal.

When addressing systems that are not well-known, maybe because we are in the early phases of a design process, one way to cope with this uncertainty is to use parameters for the weights. The interesting problems are parameter synthesis problems, in which one tries to find the values of parameters such that some path is optimal, or such that a target vertex can be reached within a given bound on the accumulated weight, etc.

Surprisingly, those problems have not been studied in detail for the setting for which cost themselves are parameters. Parametric timed automata (PTA) [START_REF] Alur | Parametric real-time reasoning[END_REF] allow clocks to be tested against parameters, which may allow simulation of parametric weights to some extent. The optimal time reachability problem has recently been studied for PTA [START_REF] André | Minimaltime synthesis for parametric timed automata[END_REF]. Similarly, the bounded-cost reachability problem has been addressed for a formalism related to PTA called parametric time Petri nets in [START_REF] Lime | Parameter synthesis for bounded cost reachability in time Petri nets[END_REF]. In [START_REF] André | Parametric analyses of attack-fault trees[END_REF], the authors do extend PTA with parametric weights on edges, but the topology of the systems considered is always that of a tree.

Much differently, the parametric one-counter machines of [START_REF] Bundala | On parametric timed automata and one-counter machines[END_REF] should allow one to model a parametric cost, however, the parameter synthesis problems are not addressed in that article.

Finally, in [START_REF] Karp | Parametric shortest path algorithms with an application to cyclic staffing[END_REF] the authors consider graphs in which weights are expressed as a function of a single parameter. They partition the real numbers in a finite way, such that for two parameter values in a given partition, the optimal paths from a given source vertex to all other vertices are the same, and exhibit the corresponding trees those paths form. Such parameteric graphs, with a single parameter, were later studied in [START_REF] Young | Faster parametric shortest path and minimum balance algorithms[END_REF] (directly improving [START_REF] Karp | Parametric shortest path algorithms with an application to cyclic staffing[END_REF]), and in [START_REF] Chakraborty | Two-phase algorithms for the parametric shortest path problem[END_REF] for example. Extending the results of [START_REF] Karp | Parametric shortest path algorithms with an application to cyclic staffing[END_REF] to parameterized graphs with multiple parameters makes the problem more complex since one needs to partition the n-dimensional real space (where n is the number of parameters). This is the subject of our work. To the best of our knowledge, this has not been done prior to this paper. The algorithm that we propose is exponential in the number of vertices and is polynomial in the number of edges, and in the number of parameters of the considered graph.

This article is organised as follows: in Section 2 we introduce the basic notations and definitions; in Section 3 we informally present our algorithm on a comprehensive example; in Section 4, we give the algorithm together with the associated proofs of correctness, completeness, and termination, as well as the complexity; in Section 5 we conclude.

Definitions

For all (a, b, i) ∈ N, we denote by i ∈ a, b any integer i such that a ⩽ i ⩽ b.

Parametric graphs Definition 1 (Parametric graph). A parametric graph with n parameters is

a tuple G = (V, E, f, (λ i ) n 1 , (Λ i ) n 1 ) where V is a set of vertices, E ⊆ V × V is a set of edges, f : E → R
is a non-parametric cost function, and for every i ∈ 1, n , λ i is a parameter, and Λ i : E → {0, 1} is a parametric cost function.

The parameters can take any value in R. With a slight abuse of notations, we denote by λ i not only the parameters but also their values. We also write -→ λ = (λ 1 . . . λ n ). All the following definitions are written for a parametric graph

G = (V, E, f, (λ i ) n 1 , (Λ i ) n 1 ).
Definition 2 (Edge cost). The cost of an edge e ∈ E is:

c(e) = f (e) - n i=1 λ i Λ i (e).
A Path p in G is a sequence of edges p = e 0 e 1 e 2 . . . e k such that ∀i ∈ 0, k , e i = (v i , v i+1 ). In such a path v 0 is the initial vertex and v k+1 is the terminal vertex. The integer k is called the length of p. Definition 3 (Path cost). Let p = e 0 e 1 . . . e k be a path of G, the cost of p is:

c(p, - → λ ) = k j=0 c(e j ) = k j=0 f (e j ) - n i=1 λ i Λ i (e j ) = k j=0 f (e j ) - n i=1 λ i k j=0 Λ i (e j ).
We call path of minimal cost for -→ λ ∈ R n a path of G such that there is no other path of G with the same initial and terminal vertices, and with a strictly lower path cost for -→ λ . Finally, we call a cycle a path where the initial vertex and the terminal vertex are the same. And we call a negative cycle for -→ λ ∈ R n a cycle with a negative cost for -→ λ .

Trees over parametric graphs

In the following we suppose that there is a distinguished vertex s in the graph G, from which we will search for paths of minimal cost toward all other vertices. As we are interested in paths from s to each vertex, in the following we assume that the graphs we consider are such that these paths exist.

Definition 4 (Tree). We define a tree of G rooted at s (or with source s) as any tuple

T = (V T , E T , f T , (λ T,i ) n 1 , (Λ T,i ) n 1 ) such that V T = V , E T ⊆ E is such that for all v ∈ V \ {s}, |{(u, v) : (u, v) ∈ E T }| = 1, {(u, s) : (u, s) ∈ E T } = ∅ and there is no cycle in T , f T = f |E T , ∀i, λ T,i = λ i , and ∀i, Λ T,i = Λ T,i|E T .
Notice that, in such a tree, there are always |V | -1 edges. Moreover, there is exactly one path from s to each vertex. In Figure 1, we have represented an example of a graph with one parameter, λ 1 (a) and a tree T of G rooted at s (b).

s u1 u2 u3 u4 1 -λ1 2 -λ1 3 -λ1 2 1 3 2 -λ1 3 (a) -A graph G s u1 u2 u3 u4 1 -λ1 2 -λ1 2 3 (b) -A tree T of G
In all the following definitions T is a tree rooted at s.

Definition 5 (Distance).

Let v be a vertex of G, the distance d(T, v, -→ λ ) is the cost of the unique path from s to v in T . Definition 6 (Partial distances). Let v be a vertex of G and let e 0 e 1 • • • e k be the unique path from s to v in T . The partial non-parametric distance between s and v in T is

d f (T, v) = k j=0 f (e j ).
The partial parametric distances between s and v are the

d Λi (T, v) = k j=0 Λ i (e j ),
for all i ∈ 1, n . Notice that d(T, v, - → λ ) = d f (T, v) - n i=1 λ i d Λi (T, v).
Definition 7 (Tree of minimal distances). We say that T is a tree of minimal distances for

- → λ ∈ R n if for all v ∈ V , the unique path from s to v in T is a path of minimal cost - → λ from s to v in G.
Moreover, if S ⊆ R n , T is a tree of minimal distances for all -→ λ ∈ S, we say that T is a tree of minimal distances on S.

Definition 8 (Neighbour)

. Let e = (u, v) be an edge in E, the neighbour of T generated by e is the tuple

N (T, e) = (V N , E N , f N , (λ N,i ) n 1 , (Λ N,i ) n 1 )
where:

-V N = V -E N = (E T \ {(u ′ , v) : (u ′ , v) ∈ E T }) ∪ {e}, f N = f |E N -∀i, λ N,i = λ i -∀i, Λ N,i = Λ i|E N
In other words, N (T, e) is obtained from T by deleting the only edge e ′ = (u ′ , v ′ ) such that v ′ = v and adding e.

Notice that for all e ∈ E T , N (T, e) = T and that, for (u, v) ∈ E \ E T , an edge (u ′ , v) does not necessarily exist in T , since T is rooted at s. Indeed an edge e such that e = (u, s) will not be in T (as T is a tree rooted at s) and can generate neighbours as any edge. In particular, this means that the neighbour of a tree is not necessarily a tree.

As an example, consider G from Figure 1 and let e = (s, u 3 ), e 1 = (u 1 , u 3 ) and e 2 = (u 4 , u 3 ). Figure 1 (b) represents a tree T rooted at s, Figure 2 (a) represents N (T, e 1 ), which is also a tree, and Figure 2 (b) represents N (T, e 2 ), illustrating the fact that a neighbour of a tree is not necessarily a tree itself.

The following proposition specifies in which case a neighbour of a tree is actually a tree. Its proof is omitted due to space constraints Proposition 1. Let e = (u, v) ∈ E be an edge, N (T, e) is a tree if and only if v is not on the unique path from s to u in T .

Before giving other properties of trees and their neighbours, we have to take a few notations. For an edge e = (u, v) ∈ E we note:

∆ f (T, e) = d f (T, u) + f (e) -d f (T, v), ∀i ∈ 1, n , ∆ Λi (T, e) = d Λi (T, u) + Λ i (e) -d Λi (T, v).
These deltas represent the differences in distance from s to v between T and its neighbour generated by e. ∆ f is the difference in the non-parametric part of the distance. Each ∆ Λi represents the difference of the number of occurrences of the corresponding parameter λ i . Notice that for any e ∈ E T , one always has ∆ Λi (T, e) = 0.

For example in Figure 1, with the same notations as before, ∆ f (T, e 1 ) =

d f (T, u 1 ) + f (e 1 ) -d f (T, u 3 ) = 1 + 3 -2 = 2, and ∆ Λ1 (T, e 1 ) = d Λ1 (T, u 1 ) + Λ 1 (e 1 ) -d Λ1 (T, u 3 ) = 1 + 1 -1.
- → λ is such that ∆ f (T, e) - n i=1 λ i ∆ Λi (T, e) = 0.
The proofs of these propositions are omitted due to space constraints.

Constraints and zones associated to trees

Definition 9 (Constraint). Let e ∈ E \ E T be an edge, the constraint associated with e is

C T,e = ∆ f (T, e) - n i=1 λ i ∆ Λi (T, e).
Definition 10 (Zone). Let E c ⊆ E \ E T be a set of edges such that ∀e c ∈ E c , ∃i ∈ 1, n , ∆ Λi (T, e c ) > 0. The zone defined by the constraints associated to the edges in E c is the set S such that:

- → λ ∈ S if and only if ∀e c ∈ E c , ∆ f (T, e c ) - n i=1 λ i ∆ Λi (T, e c ) ≥ 0.
We can notice that zones are convex by construction.

Definition 11 (Active constraint).

Let E c ⊆ E \ E T be a set of edges such that ∀e c ∈ E c , ∃i ∈ 1, n , ∆ Λi (T, e c ) > 0.
Let e c ∈ E c be an edge. Let S be the zone defined by the constraints associated to the edges in E c . Let S /ec be the zone defined by the constraints associated to the edges in E c \ {e c }. The constraint C T,ec is said to be active if and only if S /ec ̸ = S.

3 Presentation of our algorithm for minimal distances

In section 3, we propose an algorithm which, given a graph G returns a list of trees and a list of disjoint zones. Each tree T is associated with one zone S, such that T is a tree of minimal distances on S. Moreover, the union of the returned zones is the zone of R n for which there is no negative cycle in G. In this section, we demonstrate how this algorithm works on the example of G ex , a parametric graph with two parameters presented in Figure 3. The first step of the algorithm consists in finding the zone Z noCycle for which the concept of minimal distances makes sense, that is the values of λ 1 and λ 2 for which there is no negative cycle. Here, the only cycle is e 5 e 5 e 5 . . . . So, the only possible negative cycle is when the cost of e 5 is negative, therefore, when 10 -λ 1 -λ 2 < 0. Hence, we have

Z noCycle = {(λ 1 , λ 2 ) ∈ R 2 : λ 1 + λ 2 ⩽ 10}.
The goal will be to cover Z noCycle with zones associated to trees of minimal distances. For that the algorithm enumerates zones, associated with trees, going from one zone to another by considering the neighbours of the associated tree. The algorithm begins by computing a first tree T 0 . This tree is a tree of minimal distances for a particular pair (λ 1,init , λ 2,init ), where -λ 1,init = -λ 2,init = 1 + e∈E |f (e)|. Here we have (λ 1,init , λ 2,init ) = (-20, -20) and T 0 is represented in Now look at the constraints associated with the neighbours of T 0 to characterize the associated zone S 0 . The active constraints also tell which trees will be considered next. Here, the edges that generate neighbours are e 0 , e 3 , e 4 and e 5 . So, one has to look at the constraints C 0 , C 3 , C 4 and C 5 associated respectively with N (T 0 , e 0 ), N (T 0 , e 3 ), N (T 0 , e 4 ) and N (T 0 , e 5 ).

We have C 0 : 1 + λ 1 -λ 2 = 0, C 3 : 6 -λ 1 = 0, C 4 : 7 -λ 1 = 0 and C 5 : 10 -λ 1 -λ 2 = 0. Among these constraints, C 0 , C 3 and C 5 are active. Thus, we take S 0 = {(λ 1 , λ 2 ) ∈ R 2 : 1 + λ 1 -λ 2 ⩾ 0, 6 -λ 1 ⩾ 0, 10 -λ 1 -λ 2 ⩾ 0}, as represented in Figure 5.

As N (T 0 , e 0 ) and N (T 0 , e 3 ) are associated with active constraints and have not been considered yet, they are added in a list of trees to be considered later, called listT oDo. N (T 0 , e 5 ) is not added to listT oDo because it is not a tree. The pair (T 0 , S 0 ) is added to a list called listExplored. This list will be returned by the algorithm at the end of its computation.

From now on, the algorithm iteratively considers the trees of listT oDo. Assume, for example, that it begins with N (T 0 , e 0 ) = T 1 , represented in Figure 6. The edges that generate neighbours are e 1 , e 3 , e 4 and e 5 . The associated constraints are C 1 : -1 -λ 1 + λ 2 = 0, C 3 : 7 -λ 2 = 0, C 4 : 8 -λ 2 = 0 and C 5 : 10 -λ 1 -λ 2 = 0. Among them, C 1 , C 3 , and C 5 are active. From that we can define S 1 , represented in Figure 7.. As N (T 1 , e 1 ) = T 0 has already been considered, only N (T 1 , e 3 ) is added to listT oDo (recall that N (T 1 , e 5 ) is not a tree). (T 1 , S 1 ) is added to listExplored.

At the moment, listT oDo = {N (T 0 , e 3 ), N (T 1 , e 3 )} and listExplored = {(T 0 , S 0 ), (T 1 , S 1 )}. Assume that the algorithm considers N (T 0 , e 3 ) = T 2 next. This tree is represented in Figure 8. The constraints used to define S 2 are the ones associated with e 0 , e 2 and e 5 . Then, it remains to consider N (T 1 , e 3 ) = T 3 , represented in Figure 10. Three constraints are considered, this time associated with e 1 , e 2 and e 5 . All these constraints are active. The obtained S 3 is represented in Figure 11. No tree is added to listT oDo because all neighbours either have already been considered or are not trees. (T 3 , S 3 ) is added to listExplored.

At that point, the algorithm terminates because S 0 ∪ S 1 ∪ S 2 ∪ S 3 = Z noCycle . It returns listExplored. Note that listT oDo is empty.

Formal presentation of our algorithm

In this section, we formalize the algorithm that has been presented on the example in the previous section. We then prove that this algorithm is correct. 4.1 The algorithm Algorithm 1 is an algorithm that, given a parametric graph G and a vertex s, returns a list of pairs (T, S) such that every T is a tree of minimal distances on S and that the union of all S in the list is equal to the zone such that there is no negative cycle in G. In this algorithm, -→ λ init is the n-components vector so that each component is equal to

-1 + e∈E |f (e)| . Algorithm 1 trees of minimal distances of G = (V, E, f, (λ i ) n 1 , (Λ i ) n 1 )
1: let listExplored = ∅.

2: let T0 be a tree of minimal distances for

-→ λ init 3: let listT oDo = {T0} 4: while listT oDo ̸ = ∅ do 5:
choose a tree T in listT oDo (and delete it from listT oDo)

6:
let Neighbours be the set of all possible edges e such that ∃i ∈ 1, n , ∆Λ i (T, e) > 0.

7:

let S be the zone defined by the constraints associated to the edges of Neighbours In the following, we will refer to the zone where there are no negative cycles as Z noCycle and the union of the zones in listExplored as Z explored .

Algorithm 1 starts by computing T 0 (the tree from which the space will be explored) as a tree of minimal distances for -→ λ init . This initial tree is chosen to ensure that the algorithm will cover Z noCycle by only exploring toward greater λ i , i ∈ 1, n as proven in the next section. T 0 is added to listT oDo which is the list of the trees the algorithm needs to consider.

For each tree T , the algorithm begins by enumerating all the edges e that can generate a neighbour. From the associated constraints it characterizes S, a zone where T is a tree of minimal distances. Then the algorithm does two things:

(1) it adds to listT oDo all the neighbour trees that have not been considered yet (those that are not already in listT oDo or in listExplored) and ( 2) it adds the new result to the returned list by adding T and S to listExplored. We can notice that the zones in listExplored are disjoint by construction.

When listT oDo is empty the algorithm returns listExplored.

Notice that when there is only one parameter this algorithm is equivalent to the algorithm presented in [START_REF] Karp | Parametric shortest path algorithms with an application to cyclic staffing[END_REF] as constraints have a better form that make computing active constraints equivalent to look for a maximum. So it is possible to do a more efficient exploration because of the fact that there is only one dimension.

Proof of the algorithm

A first Lemma expresses that -with the value that has been chosen for --→ λ initit is not necessary to consider all the neighbours of each tree in the main loop of the algorithm. It is sufficient to consider neighbours with increasing number of occurrences of (at least) one parameter, as enforced by line 6 of Algorithm 1.

Lemma 1. Let S 0 be the first zone computed by the algorithm. ∀

- → λ ∈ Z noCycle , ∃ - → λ ′ ∈ R n and -→ λ S0 ∈ S 0 such that (1) ∃i ∈ 1, n , λ ′ i ⩾ λ i,S0 and (2
) for all trees T of minimal distances for -→ λ ′ , T is also a tree of minimal distances for -→ λ .

Proof is omitted due to space constraints A second lemma exhibits a loop invariant that will be instrumental in showing the correctness of the results of the algorithm.

Lemma 2 (loop invariant). Let Z notExplored = Z noCycle \ Z explored , at each loop of the while loop we have:

1: Z explored ∪Z notExplored = Z noCycle (in particular there is no -→ λ ∈ Z explored such that there is a negative cycle in G for -→ λ ) 2: for all (T, S) ∈ ListExplored, T is a tree of minimal distances for all -→ λ ∈ S.

Proof is omitted due to space constraints Building on the two above lemmas, we give our main theorem, that states that the proposed algorithm terminates and returns correct results. Proof. If the algorithm does not terminate it means that there is an infinite loop.

As we consider a new tree in each loop and there is a finite number of (possible) trees it is impossible to have an infinite loop so the algorithm does terminate. We also need that Z explored = Z noCycle at the end of the algorithm, i.e. when listT oDo is empty. If it is not the case it means that the algorithm has missed one or more zones and this is only possible if there is some zones such that -→ λ ∈ R n , ∄ -→ λ ′ ∈ R n and -→ λ S0 ∈ S 0 such that ∃i ∈ 1, n , λ ′ i ⩾ λ i,S0 . Which is impossible by Lemma 1 so the algorithm terminates and we have For each (T, S) ∈ ListReturned we also have that T is a tree of minimal length for all -→ λ ∈ S and

(T,S)∈ListReturned S = Z noCycle by the loop invariant.

Complexity

We conclude the presentation of Algorithm 1 by giving its worst-case complexity.

Conclusion

We have proposed an algorithm to find the optimal paths from a single source to all other vertices in a weighted graph in which weights involve an arbitrary number of real-valued parameters. Since those paths change with the values of the parameters, the result of our algorithm is a finite set of trees, each with a zone of the parameter space on which it is optimal. Those zones cover the parameter space for which there are no negative cost cycles in the graph. This algorithm generalizes a previous work by Karp and Orlin in which only one parameter was considered [START_REF] Karp | Parametric shortest path algorithms with an application to cyclic staffing[END_REF].

Further work includes implementing the algorithm, and evaluating its efficiency on real-world case-studies.

Fig. 1 .

 1 Fig. 1. A graph G and an example of a tree T of G.

Fig. 2 .

 2 Fig. 2. Some neighbours of the tree T of Figure 1.

Proposition 2 .Proposition 3 .

 23 Let e = (u, v) be an edge not in T . If N (T, e) is a tree then ∀w ∈ V, d(T, w, -→ λ ) = d(N (T, e), w, -→ λ ) if and only if ∆ f (T, e) -n i=1 λ i ∆ Λi (T, e) = 0. Let e = (u, v) be an edge. If N (T, e) is not a tree then N (T, e) has a cycle of cost 0 if and only if

Fig. 3 .

 3 Fig. 3. Graph Gex

Figure 4 .

 4 Figure 4.

Fig. 5 .

 5 Fig. 4. T0

Fig. 7 .Fig. 9 .

 79 Fig. 6. T1

Fig. 11 .

 11 Fig. 10. T3

8 :let

 8 Active be the subset of Neighbours containing the edges giving active constraints 9:for each edge e in Active such that N (T, e) is a tree do 10: let TN = N (T, e)

Theorem 1 .

 1 The algorithm terminates and returns ListReturned such that for all (T, S) ∈ ListReturned, T is a tree of minimal length for all -→ λ ∈ S, (T,S)∈ListReturned S = Z noCycle and for all (T, S) and (T ′ , S ′ ) in ListReturned such that S ̸ = S ′ , S ∩ S ′ = ∅.

  (T,S)∈ListReturned S = Z noCycle .

Theorem 2 .

 2 The worst case complexity of Algorithm 1 is exponential in the number |V | of vertices and is polynomial in the number |E| of edges and in the number n of parameters. Moreover it is logarithmic in the largest constant value M appearing in the constraints. The proof of this theorem is omitted due to space constraints and shows that this complexity is O(|V | |V |-1 (|E| -|V |) 3 n 3 log(M )).