G Rousseau
email: gregory.rousseau@ifp.fr

D Sinoquet

P Rouchon

Constrained Optimization of Energy Management for a Mild-Hybrid Vehicle

Keywords: Hybrid vehicle, Modelling, Optimal control, Optimization

Hybrid-electric vehicles appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions. The presented work focuses on a power train architecture for mild hybrid vehicles where the kinematic energy in the breaking phases is stored in a battery to be re-used later via the electric motor. This additional traction power allows to downsize the engine and still fulfill the power requirements. Moreover, the engine can be turned off in idle phases. The complete mild-hybrid vehicle is modelled in AMESim environment and the fuel consumption for given driving cycles is estimated. The control strategies for the energy management between the two power sources are optimized with respect to fuel consumption with a classical dynamic programming (DP) method. We propose an other method based on Pontryagin Minimum Principle which furnishes results very close to the DP results for a significantly reduced calculation time. These optimization results furnish the optimal control laws from which could be derived the control laws to be implemented on the vehicle. To illustrate the potential of optimization for component design, mild hybrid vehicles with varying battery and electric motor sizes, with different types of engine (gasoline / natural gas), are evaluated in terms of consumption gain with the presented methodology.

INTRODUCTION

Growing environmental concerns coupled with concerns about global crude oil supplies stimulate research on new vehicle technologies. Hybrid-electric vehicles appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions [START_REF] German | Hybrid powered vehicles[END_REF] :

• Fuel economy. It is well known that the kinetic energy gained by the classical vehicle after an acceleration can not be recuperated when braking. One of the advantage of the hybrid vehicle is the ability to recuperate a part (depending of efficiencies) of this energy, to store it in a battery, and to re-use it latter. This electric energy can therefore be used by an electric motor to power the vehicle and to save fuel. • Reduction of pollutant emissions. As using an electric motor to supply the traction power, it becomes possible to distribute the power in order to reduce the pollutant emissions. The engine can also be downsized, and still fulfill the maximum power requirements of the vehicle, which can decrease again the pollutant emissions. Last but not least, hybrid vehicles allow to eliminate the idle fuel consumption by turning off the engine when no power is required (system called "stop'n go").

The control of hybrid power trains is more complex than control of classic engine. Indeed, the control laws have to deal with the state of charge of the battery, which provides the level of energy remaining, and with the variable efficiency of each element of the power train. Optimization of energy management strategies on given driving cycles is often used to derive sub-optimal control laws to be implemented on the vehicle (see among others [START_REF] Antonio | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF], [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], [START_REF] Delprat | Evaluation de stratégies de commande pour véhicules hybrides parallèles[END_REF]).

IFP, in partnership with Gaz de France and the Ademe, has combined its downsizing technology with a natural gas engine in a small urban demonstrator vehicle (VEHGAN vehicle), equipped with a starter alternator and super capacity manufactured by Valeo ( [START_REF] Tilagone | Development of natural gas demonstrator based on an urban vehicle with a down-sized turbocharged engine[END_REF]). In this paper, we focus on this hybrid architecture which is modelled in AMESim environment. Two different optimization algorithms are presented and applied on this model:

a Dynamic Programming algorithm, classically used to handle this kind of optimal control problem ( [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF], [START_REF] Antonio | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF]),

an algorithm based on Pontryagin Minimum Principle.

In a third section, we will show some results on the architecture of the VEHGAN vehicle with different kinds of engine, different sizes of battery and of electric motor.

SYSTEM MODELLING AND OPTIMAL CONTROL PROBLEM

AMESim Model

AMESim environment combined with IFP libraries Drive, Engine and Exhaust allows to model anything from a whole vehicle to the fine behavior of the engine. Figure 1 presents the scheme of our AMESim model of the VEHGAN : a mild hybrid vehicle with a natural gas engine and a starter alternator ( [START_REF] Tilagone | Development of natural gas demonstrator based on an urban vehicle with a down-sized turbocharged engine[END_REF]). This AMESim model consists in a model of :

• the vehicle (mass, wheel inertia, resistance force) • a driver model (which can anticipate the driving cycle)

• a driving cycle to follow (imposing vehicle speed and gear shifts) • a manual gearbox coupled with a clutch • a 660cm 3 Natural gas engine, with fuel consumption map displaying in Figure 2 • a 3kW starter alternator • a battery (voltage source and an internal resistance) of 0.4Ah • and the control unit which manages the energy repartition between the engine and the electric motor. Control laws obtained by solving the optimal control problem presented in section (1.2) can be tested thanks to this control unit.

Optimal Control Problem

The optimal control problem under study consists in minimizing the fuel consumption of the vehicle along a given vehicle cycle, taking into account physical constraints from battery, engine and electric motor. The general optimization problem is the following :

8 > > > > > > > > > < > > > > > > > > > : min ( J(u) = Z T L(x(t), u(t), t)dt + g(x(T ), T )
) subject to :

ẋ = f (x(t), u(t), t), x(0) = x 0 x min x(t) xmax u min (t) u(t) umax(t) (1
Trq(t) = Te(t) + Tm(t) + T bk (t), (2)
• the engine torque is

Te(t) = u(t)Trq(t), (3)
• the electric motor torque is

Tm(t) = (1 -u(t))Trq(t) (4)
During braking, Trq(t) < 0, the braking torque completes the electric motor torque to slow down the vehicle.

T bk (t) = Trq(t) -T min m (t) (1) (5)

(1) We neglect any friction torque of the engine with T min m (t) is the minimal electric motor torque depending on the engine speed, thus depending on time.

A simple model is implemented for the battery, the state of charge of the battery evolves according to

ẋ(t) = - ω(t)Tm(t)K ′ U batt (t)ncapa (6)
with ω(t) : Electric motor and engine speed, -U batt (t) : Battery voltage, -K ′ : Scaling constants, -ncapa : Nominal capacity of the battery.

In the following, we assume that the battery voltage is constant along the cycle. We set K = K ′ U batt .ncapa and by replacing Tm(t) by its expression, we obtain

ẋ(t) = -Kω(t)(1 -u(t))Trq(t). (7)
In this optimal control problem, we make several hypothesis • the pollutant emissions are not taken into account in the optimization process, • the engine speed ωe and the electric motor speed ωm are equal, ωe = ωm = ω, • as the optimization problem is applied on a mild hybrid vehicle, the engine can not be stopped when the torque is provided only by the electric motor, except for the stop'n go mode at the idle speed. So, for a command u(t) that cancels the engine torque and for positive torque request, the fuel consumption does not vanish (Figure 2), • recharging the battery is only possible for negative torques (breaking request), we did not consider recharging by an additional engine torque beyond the driver request torque. Thus the command u(t) remains between 0 and 1, • we neglect any friction torque of the electric motor.

As in the general optimal control problem (1), the command variable u(t) but also the state variable x(t) are submitted to several constraints :

• the engine can only produce a positive torque, and is limited to a maximum torque which depends on engine speed ω(t) :

0 Te(t) T max e (ω(t)), (8)
• the electric motor has a maximum torque, and also a minimum torque during regeneration breaking :

T min m (ω(t)) Tm(t) T max m (ω(t)), (9)
• the storage capacity implies a minimum and a maximum state of charge of the battery (which are fixed to 0% and 100% in our example)

x min x(t) xmax, (10)
• the final state of charge is constrained to be equal to the initial state of charge to maintain a null electrical energy balance (to avoid to discharge totally the battery for minimizing the consumption)

x(0) = x(T ). (11)
Constraints on command u(t) are directly derived from (8), ( 9) and ( 10)

u(t)Trq(t) T max e (t), (12)
T min m (t) (1 -u(t))Trq(t) T max m (t). (13)
In the following, we will write U (t) in continuous time (respectively U k in discrete time) for the feasible domain for u(t) respecting the constraints (12) and (13).

DYNAMIC PROGRAMMING OPTIMIZATION

Classically, the Dynamic Programing method (DP) is used to solve the problem (1) ( [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF]) : it relies on a very simple idea, the principle of optimality, stated by Richard Bellman as follows : An optimal policy has the property that, whatever the initial state and initial decision are, the remaining decision must constitute an optimal policy with regard to the outcome resulting from the first decision. The principle of optimality (or Bellman principle) suggests that an optimal policy can be built step by step sequentially, first by building an optimal policy for the "tail subproblem" involving the last step, then by extending the optimal policy to the "tail subproblem" involving the last two stages, and by iterating until an optimal policy is built for the complete time interval. DP algorithm applied to our problem furnished an optimal trajectory of the state of charge of the supercapacitor. First, the optimal control problem (1) is discretized in time :

min u∈U J(u) := N -1 X k=0 L k (x k , u k ), (14)
where

L k (x k , u k ) is the cumulated fuel consumption over the time interval [k, k + 1],
x k is the state of charge of the super capacity at time k, and N is the final time of the driving cycle. The constraint on the final SOC (11) is introduced by a penalizing term g(x N ) = β.(x N -x 0 ) 2 (β, a constant to be chosen (2) ) in (14)

min u∈U J(u) := N -1 X k=0 L k (x k , u k ) + g(x N ). ( 15
)
(2) In the following results, the value of 10 has been implemented

IFP -E-COSM -(ScripT E X o570104 : 1 er specimen) -4 -
From Bellman principle, the minimum cost V k (x k ) at the time step k, 0 k N -1, is then expressed as follows :

V k (x k ) = min u k ∈U (L k (x k , u k ) + V k+1 (f k (x k , u k ))) (16)
f k being the function that modelizes the battery state of charge dynamic in the discrete form of [START_REF] Guzzella | Vehicle propulsion systems[END_REF]. At time N , the cost function is

V N (x N ) = g(x N ). ( 17
)
This optimization problem is solved backward from final time step to initial time step using a discretization of function V in the command space and in the state space. It allows a fast optimization : (N -1) constrained optimization problems (16) of one parameter, u k , to be solved with a fine time discretization of the controller. It furnishes the optimal power repartition u k at each time step regarding fuel consumption with constraints on the battery state of charge. (we refer to [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] for some theoretical results on the convergence of the method and error estimations).

A direct constrained optimization of (15) is an alternative method but leads to a large non linear optimization problem with a large number of inequality constraints. This is the reason why the DP algorithm is often preferred when the number of state variables is small (one or two state variables).

DP Optimization algorithm

As to get V k (x k ) for different values of x k , a 2-D grid (time and SOC nodes) is set. A standard time step used in our examples is 1s, and the SOC step is 0.5%. A classical DP algorithm consists in computing each feasible command (to go from a point x i k to an other point x j k+1 ), finally taking the best trajectory (3) . In such a method, the state of charge trajectory remains on the points of the defined grid in the state space. We propose, in our algorithm, to interpolate the function V (x k , k) in the state space, for each time step k. Therefore, it is possible to use a (state) continuous constrained optimization algorithm to solve each problem (16). The optimization algorithm is only used when Trq > 0 : when the request torque is negative or null, the optimal command u k is completely known, as the battery is recharged as much as possible, the command u k being constrained by the minimal electric motor torque from [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF] and by maximum SOC from (10).

Results obtained with DP method

DP algorithm described in previous section is applied on the VEHGAN vehicle problem : the optimal SOC trajectory, depending of the fuel consumption map and of the vehicle cycle is displayed for NEDC cycle in Figure 5.

(3) the trajectory which minimizes at each step k the sum Optimization results obtained with DP method are presented and commented in section 4.

(L k (x k , u k ) + V k+1 (f k (x k , u k )))

PONTRYAGIN MINIMUM PRINCIPLE OPTIMIZATION

In this section, we propose an alternative method to solve the optimal control problem (1). It relies on the Pontryagin Minimum Principle (PMP) and unlike the DP method does not require any discretization scheme.

Pontryagin Minimum Principle

The optimization problem (1) with the constraints defined in section (1.2) is rewritten as :

8 > > > > > > > < > > > > > > > : min u∈U ( J(u) = Z T 0 L(x(t), u(t), t)dt + g(x(T ), T ) ) subject to ẋ = f (x(t), u(t), t), x(0) = x 0
x min x(t) xmax (18) We introduce the Hamiltonian function H(u(t), x(t), p(t)) = L(x(t), u(t), t)+p(t) ẋ(t)+η(t)h(x(t), t), (19) p(t) is called the co-state or adjoint variable of our system. η(t) is the Lagrange multiplier associated with the state inequality constraints x min x(t) xmax. In this section, we simplify the problem by cancelling this constraint. We assume here that L is a smooth convex function of u and f is affine linear versus u.

The Pontryagin Minimum Principle states that a necessary condition of the problem (18) is the existence of p(t) such that

∂H ∂x = -ṗ. ( 20
)
If there is no constraint on command variable, we have also

∂H ∂u = 0. ( 21
)
We refer to [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF] and [START_REF] Bryson | Applied optimal control[END_REF] for further details about Pontryagin Principle.

Application

The fuel consumption L(t) to be minimized in (18), is defined by a discrete map L(ω, Te) that can be modelled by a 2-order polynomial [START_REF] Guzzella | Vehicle propulsion systems[END_REF], to be used for continuous optimization. Choosing a 2-order polynomial ensures an unique solution (under some assumptions on polynomial coefficients) and allows to model a large variety of engine maps as shown in section 4. The fuel map is therefore written as :

L(ω, Te) = 2 X i,j=0 K ij .ω i .T j e . (22)
The problem (18) is written as

8 > > > > > > > < > > > > > > > : min u∈U ( J(u) = Z T 0 L(ω(t), Te(t))dt + g(x(T ), T ) ) subject to ẋ = -Kω(t)(1 -u(t))Trq(t), x(0) = x 0
x min x(t) xmax (23) From (20) and (7) we obtain

ṗ = 0 ⇒ p = constant = p 0 . ( 24
)
Without any constraint on the state and on the command, the minimization problem can be easily solved. The minimum fuel consumption is then reached for u * so as

∂H ∂u = ∂L ∂u + p ∂f ∂u = 0. ( 25
)
The optimal command u * can be calculated easily by solving the equation (25), which depends linearly on u. u * finally depends on p(t), Trq(t) and ω(t)

u * (t) = - 2 X i=0 K i1 ω(t) i + p 0 .K.ω(t) 2 2 X i=0 K i2 ω(t) i .Trq(t) . ( 26
)
The expression of p 0 (27) is obtained by replacing u * (t) by its expression in the state equation [START_REF] Guzzella | Vehicle propulsion systems[END_REF], and by integrating this equation in time

p 0 (τ ) = x(0) -x(τ ) - Z τ 0 (K.ω(t)(Trq(t) + F (t))dt Z τ 0 ω(t) 2 K 2 2 2 X i=0 K i2 ω(t) i .Trq(t) dt , (27) with F (t) = 2 X i=0 K i1 ω(t) i 2 2 X i=0 K i2 ω(t) i .Trq(t) , ( 28
)
and τ is the considered final time which can be set to a different value than T as we will see in next section.

Handling constraints on control and state variables

The previous section presents the computation of the optimal control of the continuous problem in a restricted case where no constraint is introduced. While control constraints are generally easily taken into account, handling the state constraints is cumbersome in the continuous optimal control problem. To handle control constraints, several singular cases can be found in [START_REF] Bryson | Applied optimal control[END_REF]. In our application, we are not able to find an analytic solution of the optimal control problem with control constraints : indeed, these constraints depends on time (via their dependency on engine speed) and depends on p 0 which depends on final SOC (cf. previous section). By an iterative method (called algo1 in the following), we can compute the value of p 0 in order to reach the desired SOC at final time with the control, computed thanks to (26), projected on bound constraints (coming from the feasible domain U ).

Concerning the inequality state constraints, [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF], [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF], [START_REF] Evans | An introduction to mathematical optimal control theory[END_REF], [START_REF] Bryson | Applied optimal control[END_REF], [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] have studied the general problem (18). In our application, we can show that p(t) presents discontinuities at the time steps where the state inequality constraints are saturated. The time steps where such discontinuities occur are not a priori known : this prevents us to solve explicitly the continuous optimal control problem with these state constraints.

PMP Optimization algorithm

Considering the difficulties described in the previous section, we propose a heuristic iterative method that allows to find a sub-optimal trajectory from the constrained continous optimal control problem (18). The proposed algorithm consists in an initialization step and 3 steps : The proposed algorithm based on Pontryagin Minimum Principle. 0. algo1 is applied on the NEDC driving cycle [0, T ] (see Figure 6 Step 0). The obtained optimal trajectory violates the state constraints, the farthest SOC (ie the "most violated point") from the bounds being for instance on the figure at point (x(tv) = -37%, tv = 818s). The initial time is called t i , here set to 0.

1. The SOC at tv is projected on the nearest bound of the feasible state domain (for instance, SOC is fixed to x min = 0 at point tv). 2. algo1 is applied again on [t i , tv] (see Figure 6 Step 2). If the obtained trajectory still violates the state constraints on [t i , tv], steps 1 and 2 are applied again on the farthest SOC from the bounds (defining a new point tv). This procedure is repeated until the trajectory remains on the feasible domain. Then the last point tv becomes the new initial time t i in step 3. 3. algo1 is applied on [t i , T ] (see Figure 6 Step 3). If the obtained optimal trajectory still violates the state constraints, steps 1 and 2 are repeated. This sequence is repeated until we reach the final step T at the desired final SOC, without violating the state constraints (Figure 6 bottom right).

PMP compared to DP Results

Using the same fuel consumption map and vehicle cycle than with Dynamic Programming, we can compare the two optimization algorithms along the NEDC driving cycle on Figure 7.

The two curves are very similar, the commands u(t) being very similar too. We now compare the two algorithms on Urban Artemis driving cycle [START_REF] André | The artemis european driving cycles for measuring car pollutant emissions[END_REF], with the same fuel consumption map, on Figure 8.

We notice again that the 2 curves are almost superimposed. In the following will be presented some results obtained with Pontryagin Principle algorithm, as it is about 1000 times faster than DP (DP algorithm computation time is about 1 hour), and gives quite the same results than DP. The small differences between the 2 algorithms can be minimized by taking smaller time and space discretization steps in DP, but it leads to a too expensive computation time.

RESULTS

In this section we present results for different kinds of engine, different sizes of battery, and different sizes of electric motor. As our mild-hybrid vehicle is a urban vehicle, we will present results on the urban Artemis driving cycle. For all simulations, we will assume that initial SOC is 80%, final SOC of 80% being introduced as a constraint to be close to initial SOC.

Different kind of engine

The both following results, natural gas engine and gasoline engine, were obtained with a battery capacity of 0.4 Ah, and a 3 kW electric motor, as on the VEHGAN vehicle, on Figure 9 and 11.

Natural Gas engine

The first engine used is a natural gas engine, from the VEHGAN project. Urban Artemis driving cycle with natural gas engine : optimal state of charge trajectory (Top); electric motor torque (Bottom)

As the natural gas engine has an efficiency growing with the torque (the best operating points being between 2500 and 4500 rpm and over 60Nm, see Figure 2), the optimization decreases the fuel consumption using the electric motor in the worth efficiency points, i.e. when ω < 2500 rpm. As using electric motor on low torque operating points would lead to worth operating points as illustrated on Figure 10, only high torque operating points are modified by an active hybrid mode.

Small Gasoline Engine

We now use a gasoline fuel map to optimize the fuel consumption for a 600cm 3 gasoline engine .

As displayed on Figure 12, the optimized operating points are quite the same than ones of natural gas engine. Nevertheless, as the gasoline engine does not present a good efficiency at high torque operating points, the electric motor is also useful in high torque requests.

Urban Artemis driving cycle with gasoline engine : optimal state of charge trajectory (Top); electric motor torque (Bottom)

Different sizes of battery

In this section we optimize the fuel consumption for the natural gas engine, with different sizes of battery. As we can see on Figure 9, the state of charge crosses the point x min = 0 %, meaning that the whole range of the state of charge of the battery is used. Then, if we increase the battery capacity, the optimization algorithm should give better results in term of fuel consumption. With a battery capacity of 0.4 Ah, we obtain a fuel consumption of 5.90 l/100km, while a pure thermic mode consumption gives 6.08 l/100km on the whole cycle, so we decreased the consumption of 2.9 %. Optimizing the torque repartition with a battery capacity two times bigger, i.e 0.8 Ah, we obtain the same fuel consumption : 5.90 l/100km. We obtain the same consumptions with the two battery sizes because of the low power of the electric motor, which is often saturated at regenerative breaking.

Increasing the Electric motor torque power

From the conclusion of these last results, we multiply by two the minimum and maximum electric motor power, to take advantage of a higher energy when regenerative breaking. We also increase the size of the battery, to have coherent battery capacity regarding to regenerative breaking and maximum electric motor power. Figure 14 presents the SOC evolution and the electric motor torque. The fuel consumption with this new electric motor reaches 5.87 l/100km, i.e Urban Artemis driving cycle with natural gas engine, a 6 kW electric motor and a 1.2 Ah battery capacity : optimal state of charge trajectory (Top); electric motor torque (Bottom) Results with different sizes of battery and of electric motor on Artemis Urban Cycle.

Discussion

The previous results show that it is possible to decrease by 2 to 4% of fuel consumption along a urban driving cycle with this mild hybrid vehicle. Activating the stop'n go mode would save around 10 to 15 %, but it would also decrease the available electric energy from battery, as this energy must be used to turn on the engine after a stop.

CONCLUSIONS

In this study, we have presented two of the three most used methods in optimal control optimization. The heuristic method based on Pontryagin Minimum Principle, well known in the free state constraint case, has been applied successfully to our state constrained problem, with very similar results to Dynamic Programming methods and a computation time divided by 1000. Nevertheless, there is currently no theoretical proof to confirm the presented validation results. Moreover, there are some limitations to this approach, mainly the assumptions on the fuel consumption map, modelled by a smooth convex function of control u (2-order polynomial) ; this limitation could lead to a bad approximation of the real fuel consumption for some particular engines.

The case study, a mild hybrid vehicle, does not present as good results as could be obtained with a full hybrid, in terms of fuel consumption gain. So in future work will be studied optimization applied to full-hybrid vehicles, as the ability to turn off the engine while the electric motor only power the vehicle is a key advantage. Other degrees of freedom, as the gear-shifting sequence should also be taken into account in the optimization problem. Reduction of pollutant emissions will also be studied by considering a second state.

From optimization results could be derived a suboptimal feedback law [START_REF] Antonio | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF], [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF] based on state of charge and/or speed measurements.

Figure 1 AMESim

 1 Figure 1 AMESim Model of the VEHGAN vehicle

Figure 2 Fuel

 2 Figure 2 Fuel consumption map of natural gas engine of VEHGAN vehicle

Figure 3 EngineFigure 4

 34 Figure 3Engine and electric motor torque simulated with AMESim model of VEHGAN vehicle of Figure1

Figure 5 NEDC

 5 Figure 5 NEDC cycle(Top); Optimal state of charge trajectory of VEHGAN vehicle computed for NEDC cycle with DP algorithm (Bottom)

  Figure 6

Figure 7 NEDCFigure 8 Urban

 78 Figure 7 NEDC cycle(Top); Optimal state of charge trajectory of VEHGAN vehicle computed for NEDC cycle with PMP & DP algorithm (Bottom).

  Figure 9

Figure 13 Urban

 13 Figure 13Urban Artemis driving cycle with natural gas engine and a 0.8 Ah battery capacity : optimal state of charge trajectory (Top); electric motor torque (Bottom)

Figure 15

 15 Figure 15

IFP -E-COSM -(ScripT E X o570104 : 1 er specimen) -2 -

IFP -E-COSM -(ScripT E X o570104 : 1 er specimen) -11 -

ACKNOWLEDGMENTS

We would like to thank Gilles Corde, Philippe Moulin and Antonio Sciarretta for helpful discussions and advice at various stages of the elaboration of this work. We acknowledge Jean-Charles Dabadie and Luc Voise for their precious help on Amesim modelling, and Quang Huy TRAN for his advice on numerical methods. We would like to thank the referee Antonio Sciarretta and Gérard Gissinger for their constructive remarks on this paper.

Operating points of engine in hybrid mode obtained by PMP algorithm for the urban Artemis Driving Cycle with a gasoline engine.

IFP -E-COSM -(ScripT E X o570104 : 1 er specimen) -10 -


J Doridam

A Macron

C Vergari

A Verney

P.-Y Rohan

H Pillet

Feasibility of sub-dermal soft tissue deformation assessment using B-mode ultrasound for pressure ulcer prevention

Keywords:

Pressure Ulcer (PU) prevention remains a main public health issue. The physio-pathology of this injury is not fully understood, and a satisfactory therapy is currently not available. Recently, several works suggested that mechanical strains are responsible of deformation-induced damage involved in the initiation of Deep Tissue Injury (DTI). A better assessment of the internal behavior could allow to enhance the modeling of the transmission of loads into the different structures composing the buttock. A few studies focused on the experimental in vivo buttock deformation quantification using Magnetic Resonance Imaging (MRI), but its use has important drawbacks. In clinical practice, ultrasound imaging is an accessible, low cost, and real-time technic to study the soft tissue. The objective of the present work was to show the feasibility of using B-mode ultrasound imaging for the quantification of localised soft-tissue strains of buttock tissues during sitting. An original protocol was designed, and the intra-operator reliability of the method was assessed. Digital Image Correlation was used to compute the displacement field of the soft tissue of the buttock during a full realistic loading while sitting. Reference data of the strains in the frontal and sagittal planes under the ischium were reported for a population of 7 healthy subjects. The average of shear strains over the region of interest in the fat layer reached levels up to 117% higher than the damage thresholds previously quantified for the muscular tissue in rats. In addition, the observation of the muscles displacements seems to confirm previous results which already reported the absence of muscular tissue under the ischium in the seated position, questioning the assumption commonly made in Finite Element modeling that deep tissue injury initiates in the muscle underlying the bone.

Introduction

Pressure Ulcer (PU) prevention remains a main public health issue, by its human and financial cost with extended hospitalization, decreased quality of life and loss of autonomy [START_REF]Agence Nationale d'Accreditation et d'Evaluation en Sante. Prevention and treatment of pressure ulcers of adults and elderly persons[END_REF][START_REF] Van Nie-Visser | An International prevalence measurement of care problems: study protocol[END_REF][START_REF] Cogan | Systematic review of behavioral and educational interventions to prevent pressure ulcers in adults with spinal cord injury[END_REF]. Two types of PU can be identified: one superficial developing from the skin and another developing from the deep soft tissue which, in an advanced stage, progressively extends to the skin [START_REF] Ceelen | Validation of a numerical model of skeletal muscle compression with MR tagging: a contribution to pressure ulcer research[END_REF][START_REF] Ankrom | Pressure-related deep tissue injury under intact skin and the current pressure ulcer staging systems[END_REF][START_REF] Bouten | The etiology of pressure ulcers: skin deep or muscle bound?[END_REF][START_REF] Gefen | The biomechanics of sitting-acquired pressure ulcers in patients with spinal cord injury or lesions[END_REF]. These deep lesions, called Deep Tissue Injuries (DTI), are more severe because of the late time of diagnosis and sometimes may require surgical intervention. Spinal Cord Injury (SCI) patients are especially concerned. The development of PU due to reduced mobility and neurosensorial loss remains frequent in these patients, with a prevalence of 33% in the first year, and an incidence of 80% during the whole life [START_REF] Garber | Pressure ulcer risk in spinal cord injury: Predictors of ulcer status over 3 years[END_REF][START_REF] Krause | Patterns of recurrent pressure ulcers after spinal cord injury: Identification of risk and protective factors 5 or more years after onset[END_REF]. Due to prolonged sitting position and wheelchair use, the pelvic area is the most frequent location of PU in these patients, with about 64% in the ischiatic area [START_REF] Gefen | The biomechanics of sitting-acquired pressure ulcers in patients with spinal cord injury or lesions[END_REF][START_REF] Richardson | Prevalence and incidence of pressure sores in acute spinal cord injuries[END_REF]. Recommendations for PU prevention have been established in clinical practice based on pressure relief strategies with patient mobilization and specific pressure-relieving cushion [START_REF]Agence Nationale d'Accreditation et d'Evaluation en Sante. Prevention and treatment of pressure ulcers of adults and elderly persons[END_REF][START_REF] Brem | Protocol for the successful treatment of pressure ulcers[END_REF]. Despite these, PU prevalence still remains high [START_REF] Van Nie-Visser | An International prevalence measurement of care problems: study protocol[END_REF][START_REF] Cogan | Systematic review of behavioral and educational interventions to prevent pressure ulcers in adults with spinal cord injury[END_REF]. More efficient prevention necessitates a better understanding of the underlying physiopathological mechanism.

PU are provoked by compression of the soft tissues between an external surface and a bony prominence. Compression-induced ischemia is known to be a major factor of PU onset [START_REF] Stekelenburg | Role of ischemia and deformation in the onset of compressioninduced deep tissue injury: MRI-based studies in a rat model[END_REF][START_REF] Loerakker | The effects of deformation, ischemia, and reperfusion on the development of muscle damage during prolonged loading[END_REF]. More recently, several experimental works suggested that damage involved in the initiation of DTI are tightly related to localized soft tissue strains [START_REF] Loerakker | The effects of deformation, ischemia, and reperfusion on the development of muscle damage during prolonged loading[END_REF][START_REF] Ceelen | Compression-induced damage and internal tissue strains are related[END_REF][START_REF] Gawlitta | The relative contributions of compression and hypoxia to development of muscle tissue damage: an in vitro study[END_REF]. Indeed, Stekelenburg et al. showed for the left tibialis of rats that the combination of deformation with ischemia decreased the time for irreversible damage compared to ischemia alone [START_REF] Stekelenburg | Role of ischemia and deformation in the onset of compressioninduced deep tissue injury: MRI-based studies in a rat model[END_REF]. Moreover, the damages were evident when the shear strains (Cauchy Green formulation) endured by the muscle were above 75% [START_REF] Ceelen | Validation of a numerical model of skeletal muscle compression with MR tagging: a contribution to pressure ulcer research[END_REF]. Gefen et al.

T

further showed that at 77% of deformation, instantaneous cell death could be observed in bio-artificial muscles [START_REF] Gefen | Strain-time cell-death threshold for skeletal muscle in a tissue-engineered model system for deep tissue injury[END_REF].

Few works focused on in vivo buttock deformation quantification by MRI imaging. Sonenblum et al. used 3D seated Magnetic Resonance Imaging (MRI) to describe the anatomy and deformation of the buttocks for 4 able-bodied and 3 SCI subjects during sitting in loaded and unloaded conditions [START_REF] Sonenblum | 3D anatomy and deformation of the seated buttocks[END_REF]. They highlighted important variation between soft tissue behaviors under the ischial tuberosity (IT) in the different subjects. The authors showed that, in the loaded configuration, 5 subjects out of 7 only had fat tissue beneath the IT, the Gluteus Maximus (GM) being positioned more laterally and posteriorly to the IT. In the loaded configuration, only one of the subject had a GM muscle thickness of 8 mm beneath the IT. The authors highlighted that this observation was inconsistent with the assumptions dominating the literature, focusing on localized muscle strains either in experimental studies on rats or in numerical simulations in humans. Yet, the reported results suggest that the muscle, particularly the GM is often not directly loaded by the bone. Likewise, Call et al. observed, using MRI, that, in a population of 10 spinal cord injury subjects and one able-bodied person, no muscle were present beneath the IT for 40% of the SCI subjects during sitting [START_REF] Call | Off loading wheelchair cushion provides best case reduction in tissue deformation as indicated by MRI[END_REF]. They concluded that this lack of muscle in the SCI population compared to healthy people decreased the overall soft tissue thickness, increasing the risk of tissue damage over bony prominences. None of these studies quantified internal deformation. On the contrary, Solis and al. assessed internal deformation in vivo in 5 healthy and 2 SCI pigs using tagged MRI [START_REF] Solis | Distribution of internal strains around bony prominences in pigs[END_REF]. They demonstrated that relatively low levels of external loads are capable of generating high magnitudes of strain in the loaded muscles, particularly in atrophied muscles. They also observed that compressive strains were higher in the SCI pigs than in the healthy sample. However, this technic is based on the assumption of small deformations, which is not the case in biological soft tissues in normal loading. As the authors reported, when the deformation in the tissue becomes too large, the tag lines become indiscernible and strain vectors can no more be measured. It is important to stress that Solis et al. are the only research team to have quantified localized soft tissue strains in the buttock region in vivo but only on an animal model. In short, if MRI is a potential tool for the quantitative evaluation of subdermal soft tissue strains, it has important drawbacks including long acquisition time, examination cost and confined environment.

Because of the experimental difficulties associated with the estimation of the localised strains deformation in soft tissue during sitting, some research teams have developed finite element models of the buttock region [START_REF] Luboz | Personalized modeling for real-time pressure ulcer prevention in sitting posture[END_REF][START_REF] Linder-Ganz | Patient-specific modeling of deep tissue injury biomechanics in an unconscious patient who developed myonecrosis after prolonged lying[END_REF][START_REF] Reed | Development and validation of a high Anatomical Fidelity FE model for the buttock and thigh of a seated individual[END_REF]. FE is a powerful tool to improve the understanding of pressure ulcer formation near the ischium because it allows to study the impact of relative influence of individual parameters, soft tissue mechanical properties and testing in many conditions. Even if partial validation based on external envelope shape or external pressure measurement [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF] were reported, no direct validation of localized soft tissue strains has been performed. Yet, these localized soft tissue strains are paramount because they are currently the only mechanical parameter known to be responsible of the onset of DTI. Moreover, the assumptions on which these models rely often simplified both the anatomy and the relative interaction between the internal components such as the sliding between muscles. A better assessment of the internal behavior could allow to enhance the modeling of the transmission of loads into the different structures composing the buttock.

In clinical practice, ultrasound imaging is an accessible, low cost, and real-time technic to study the soft tissue [START_REF] Klauser | Sonoelastography: musculoskeletal applications[END_REF]. Moreover, ultrasound imaging has been shown to be promising for the quantification of localized soft tissue strains using Digital Image Correlation (DIC) techniques [START_REF] Gennisson | Ultrasound elastography: Principles and techniques[END_REF]. It has been used, for example, in vitro, in tissue-mimicking phantom [START_REF] Zhu | An efficient optimal method for a 2D strain estimation of ultrasound tissue-mimicking material phantom[END_REF] and, in vivo, in the human Achilles tendon [START_REF] Chimenti | Ultrasound strain mapping of Achilles tendon compressive strain patterns during dorsiflexion[END_REF] and in the quadriceps muscle [START_REF] Affagard | Measurement of the quadriceps muscle displacement and strain fields with ultrasound and Digital Image Correlation (DIC) techniques[END_REF]. The extension to buttock tissue should contribute to a better understanding of the behavior of subdermal tissues in this region during loading.

The objective of the present work was to determine the feasibility of B-mode ultrasound imaging for the quantification of localized soft tissue strains of buttock soft tissues during sitting. The intra-operator reliability of the method will be assessed and reference data on a population of 7 healthy subjects will be given.

Materials and method

Subjects and preparation

Seven subjects participated to the experiment, 1 woman and 6 men (Age: 28 ± 4 yrs, Weight: 75 ± 10 kg, BMI: 23.2 ± 2.0 kg/m 2 ). After approval by the ethics committee (Comité de protection des Personnes CPP NX06036) and written informed consent, the subjects were equipped with 4 reflective markers on the anterior and posterior iliac spines of the pelvis and one on the spinal process of the 4 t h lumbar vertebra.

Ultrasound acquisitions

An ultrasound acquisition of the subdermal tissue in the region beneath the ischial tuberosity was performed using a commercial device (Aixplorer, SuperSonic Imagine, France) with a linear ultrasound probe of 8 MHz central frequency (SuperLinear SL 15-4). Simultaneously, the 3D position of the markers were recorded using a optoelectronic system (Vicon, Vicon Motion Systems Ltd., Oxford, UK).

A custom-made stool was used for the experiment. A cross-shaped notch was made in the seat to fix the ultrasound probe in two orthogonal positions (Fig. 1). In order to ensure a good fit between the crossshaped notch and the probe, the space was filled with a custom made 3D printed part. Special care was paid to the design so as to ensure a flat surface in the vicinity of the probe. The stool was equipped with 3 optoelectronic markers.

The subject was asked to sit on the stool in front of the on-screen display of the ultrasound device. After covering of the probe with echography gel, two video clips were acquired: the first one with the probe in the frontal plane and the second one in the sagittal plane (Fig. 1). For each acquisition, the subject was instructed to align his left ischium above the probe, to adjust his position to see the bone at the center of the screen and then to slowly unload his weight with his arms while visually keeping the ischium as aligned as possible with the probe.

Data analysis 2.3.1. Image processing and strain calculation

Digital Image Correlation (DIC) was used to quantify localized soft tissue strains in the muscle and fat layers using a custom software written in Matlab (2014b) (MathWorks, Inc, Natick MA) [START_REF] Vergari | Lamellar and fibre bundle mechanics of the annulus fibrosus in bovine intervertebral disc[END_REF]. First, a rectangular Region Of Interest (ROI) was defined manually in the first image and discretized into a grid of 21 points (3 × 7). These points were then automatically tracked in each successive frame using a normalized 2-D cross-correlation algorithm. The position of the points could be manually corrected if necessary. From the positions of the points, a displacement field was calculated at each frame as described below.

The deformation gradient tensor F and the Green-Lagrange strain tensor E were computed from the displacement field for each point of the grid as following: (1)

= E F F I 1 2 ( ) T (2) = E E E 1 2 Shear I II (3)
Here, x i is the position of the grid points in the ith image and X in the initial image. The Green Lagrange tensor E were diagonalized to quantify the first and the second principal strains: respectively E I and E II . These were used to estimate the shear strains E Shear using equation [START_REF] Cogan | Systematic review of behavioral and educational interventions to prevent pressure ulcers in adults with spinal cord injury[END_REF].

The first and the second principal strains in the ultrasound frontal and sagittal planes are illustrated in Fig. 2.

Reproducibility

Intra-operator reproducibility of the shear strain was also assessed. For each subject and for each acquisitions (in respectively the frontal and sagittal planes), one operator performed 4 times the whole procedure. Each time, the initial grid was automatically positioned in the fat layer. For each point i of the grid, the variance (S i 2 ) among the trials was calculated as follows:

= x E E ik ik i (4) = = S x x N ( ) i k N ik i i 2 1 2 (5)
where E ik is the shear strain of the ith point for the kth trial, N the number of trials and E i the mean shear strain of the ith point among the trials.

Si was computed as an estimator of the error of reproducibility.

Calculation of the horizontal rotation of the pelvis

In order to quantify the change in orientation of the pelvis between the two ultrasound acquisitions, a local pelvis frame (R p ) was constructed from the anterior and posterior iliac spines markers. For each acquisition, the position of this frame was expressed in a local frame attached to the stool (R c ) and the orientation matrix of the relative position of R p between both acquisitions was calculated. Then, the angle around the vertical axis of the pelvis between both configurations was computed for each subject by decomposition of the rotation matrix.

Results

Calculation of the horizontal rotation of the pelvis

For all the subjects, the horizontal rotation of the pelvis between the frontal and sagittal acquisitions was below 8°. Strain measurement was possible for all subjects in the fat layer. Due to the limited acquisition field inherent to the use of ultrasound probes, the monitoring of the muscle tissue during the whole acquisition was not always possible. The distribution of the shear strain E Shear is given for the 7 subjects in Figs. 5 and6 (respectively in the frontal and sagittal planes) below.

Distribution of the internal deformation of the soft tissues

Except for subject 1 and 6, the shear strain was larger in the sagittal plane as compared to the frontal plane. On average, the shear strain was equivalent for subject 1 in both planes and larger in the frontal plane for subject 6. Results also show a high inter-individual variability in the mechanical response.

Table 1 below gives the average values of the strains in the first and second principal directions over the grid for each subject for both planes. For the first principal direction, it can be observed that strains in both planes are of the same intensity. It can also be observed that the strains in the first principal direction (compression) are lower than the strains in the second principal direction (tension in frontal and sagittal planes).

Intra operator reproducibility

The distribution of the error of reproducibility S i of the shear strain is reported in Fig. 3 below in the form of a box-plot in the frontal and sagittal planes for each subject. The average error was below 2% in the frontal plane and 3% in the sagittal plane.

Discussion

The present work aimed to show the feasibility of using B-mode ultrasound imaging for the quantification of internal soft-tissue strains of buttock tissues in two perpendicular planes during sitting. A methodology was proposed for measuring subdermal soft tissue displacements and mapping strains in the frontal and sagittal planes respectively using DIC. Our results show that soft tissue in the fat layer can be tracked. On the contrary, the soft tissue displacement in the muscle region beneath the ischial tuberosity could not always be assessed due to the limited lateral field of view of the device. The method included a manual correction phase which was performed by the same operator for all the subjects. To ensure the reliability of the strain measurements, the repeatability of this operator was also quantified and shown to be inferior to 3% for the shear strain in the fat tissue in both the frontal and sagittal planes.

If the quantification of the soft tissue strain field has gained increasing interest in static, dynamic and transient elastography for the clinical diagnosis of various diseases in the last years, this is, as far as the authors are aware of, the first time that B-mode ultrasound imaging is used for the quantification of subdermal soft tissue strains in the buttock region. Previous studies focused on the ultrasound strain mapping of the human Achilles tendon [START_REF] Chimenti | Ultrasound strain mapping of Achilles tendon compressive strain patterns during dorsiflexion[END_REF], quadriceps muscle [START_REF] Affagard | Measurement of the quadriceps muscle displacement and strain fields with ultrasound and Digital Image Correlation (DIC) techniques[END_REF], carotid artery [START_REF] Gahyun | Strain elastography using dobutamine-induced carotid artery pulsation in canine thyroid gland[END_REF] and breast [START_REF] Mutala | Comparison of qualitative and semiquantitative strain elastography in breast lesions for diagnostic accuracy[END_REF] for example. The assessment of the strains in the buttock region implies several specific difficulties due to the limited accessibility to the zone, the deep position of the osseous structures and the necessity to localize the ischial tuberosity. Other attempts have been made with MRI to assess the deformation of the tissue during sitting. Most of these studies did not report local strains [START_REF] Sonenblum | 3D anatomy and deformation of the seated buttocks[END_REF][START_REF] Call | Off loading wheelchair cushion provides best case reduction in tissue deformation as indicated by MRI[END_REF]. Only Solis et al. used tagged MRI to quantify internal deformation in the buttocks of pigs [START_REF] Solis | Distribution of internal strains around bony prominences in pigs[END_REF]. However, this technic cannot be easily transposed in humans while sitting because of several limitations. First, the sequential acquisition of slices requires relatively large acquisition time at each load step, which prevents obtaining a continuous assessment over the sitting motion. Second, as the technic relies on the assumption of small deformations, the assessment of internal deformation during the sitting would necessitates a large number of loading step given the large deformation experienced by buttock soft tissue during sitting. Third, the initial unloaded sitting position in the MRI can only be simulated [START_REF] Reed | Deformation of the gluteal soft tissues during sitting[END_REF]. On the contrary, the present work allowed to image subdermal soft tissue response during the full progressive loading of the buttock during sitting motion and, from these, to compute internal soft-tissue strains of buttock in two perpendicular planes in a population of 7 healthy subjects.

In the present contribution, the monitoring of the muscle tissue during the whole acquisition was not always possible due to the limited acquisition field inherent to the use of ultrasound probes. Moreover, in both the frontal and sagittal planes, the muscle tissue motion in the second principal direction (perpendicular to the pelvis motion) was important suggesting there was a non-negligible out-of-plane motion of the material particles. As a result, tracking muscle features using image registration techniques in each plane would introduce biases. Accordingly, computing muscle tissue strains of buttock was not reported. Yet, the analysis of the videos allowed to observe the complex

Table 1

Average and standard deviation of the principal strains in the frontal and sagittal planes respectively. behaviour of this tissue layer during the loading. In particular, it was observed that the gluteus maximus and hamstring muscles did not always remain underneath the IT, but instead, in certain subjects, slided out of plane, leaving the IT on the fat tissue. This observation is close to previous reported results. Sonenblum et al. observed that the gluteus maximus had a posterior and lateral displacement between simulated unloaded and loaded sittings in an MRI [START_REF] Sonenblum | 3D anatomy and deformation of the seated buttocks[END_REF][START_REF] Sonenblum | 3-dimensional buttocks response to sitting: a case report[END_REF]. This contributes to the decrease of the thickness of soft tissue muscle under the ischial tuberosity. These results should be further investigated because, if confirmed, it could question the assumptions, commonly made by teams developing finite element model of the buttock that deep tissue injury initiates in the muscle underlying the bone. Because the localised tissue strain is obviously directly linked to the morphology of the bones future work focusing on the relationship between localised soft tissue strains and bone morphology and surface orientation is also useful for better understanding the determinants of the onset of PU.

The values of shear strains in the fat layer reported in the present contribution in two perpendicular planes can reach levels higher than the damage thresholds determined in vitro on rat models [START_REF] Loerakker | Temporal effects of mechanical loading on deformation-induced damage in skeletal muscle tissue[END_REF]. It must be noticed that these thresholds concern muscle tissues. However, in a case study, Sonenblum et al. observed that a patient with a history of pelvic pressure ulcer had no muscle loaded while sitting, the gluteus maximus resting entirely posterior and lateral to the ischium [START_REF] Sonenblum | 3D anatomy and deformation of the seated buttocks[END_REF]. This suggest that pressure ulcers can develop in the fat tissue. Considering the level of strains quantified in our study, we could expect that the strain damage thresholds would be different for the fat than for the muscle. Further studies are necessary to evaluate these thresholds.

The proposed methodology also allowed to quantify the deformation in two orthogonal planes. The quantification of the rotation of the pelvis between the two acquisitions confirm that the deviation from the orthogonal assumption could be neglected. Although analysis of full 3D strain field is more appropriate to characterize soft tissue deformation, the proposed protocol is a good alternative easier to implement and process. The results obtained also underlined that the shear strain intensity is of comparable magnitude in both the frontal and sagittal planes, which disconfirmed the assumption of plane deformation, usually made in finite element simulations.

In addition to the improvement of the proposed methodology to address the biases discussed in this section, future work should also focus on the validation of the localised strain field. One promising avenue is the adaptation of the "truth cube" method proposed by Kerdok et al. [START_REF] Kerdok | Truth cube: Establishing physical standards for soft tissue simulation[END_REF] where a physical standard (cube of silicone rubber with a pattern of embedded Teflon spheres loaded in uniaxial compression) was developed to validate a FE model. This raises several issues related to the realistic representation of the complexity of the subdermal soft tissues in silico and which will have to be addressed.

To conclude, this study gives one of the first experimental evaluation of the internal deformation in the buttock region in humans during a real sitting motion. One of the future perspective to assess deformation in the muscle would be to combine this protocol with an acquisition of the load exerted on the buttock to be able to assess the deformation at intermediate strain ranges when muscles remain visible in the field of view of the ultrasound device.

Fig. 4

 4 Fig. 4 below shows the deformation of the subdermal soft tissues during sitting for one typical subject.Strain measurement was possible for all subjects in the fat layer. Due to the limited acquisition field inherent to the use of ultrasound probes, the monitoring of the muscle tissue during the whole acquisition was not always possible.

Fig. 1 .Fig. 2 .

 12 Fig. 1. a) Diagram showing the position of the subject during the ultrasound acquisition b) and c) sagittal and frontal positions of the ultrasound probe with respect to the pelvis.

Fig. 3 .Fig. 4 .

 34 Fig. 3. a) box-plot of the error of reproducibility S i of the shear strain in the (a) frontal plane and (b) sagittal plane.

Fig. 5 .

 5 Fig. 5. Frontal plane. (a) initial grid superposed on the initial image (b) shear strain at each point of the grid for one typical subject (c) for each subject (abscissa 1 to 7 of the graph) a group of 7 bars gives the shear strain measured in each point of the grid.

Fig. 6 .

 6 Fig. 6. Sagittal plane. (a) initial grid superposed on the initial image (b) shear strain at each point of the grid for one typical subject (c) for each subject (abscissa 1 to 7 of the graph) a group of 7 bars gives the shear strain measured in each point of the grid.

Acknowledgments

This work was supported by the Fondation de l'Avenir (grant number AP-RM-2016-030), by la Fondation des Arts et Métiers and the Fond de dotation Clinatec. The authors are also grateful to the ParisTech BiomecAM chair program on subject-specific musculoskeletal modeling.

Conflicts of interest

The authors certify that no conflict of interest is raised by this work.


Gregory Rousseau

Delphine Sinoquet

Pierre Rouchon

HYBRID ELECTRIC VEHICLES : FROM OPTIMIZATION TOWARD REAL-TIME CONTROL STRATEGIES

Keywords: Hybrid vehicle, Optimal control, Dynamic programming, Pontryagin, Control strategies

Hybrid-electric vehicles appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions. The presented work focuses on two types of architecture : a mild hybrid and a full hybrid where the kinetic energy in the breaking phases is stored in a battery to be re-used later via the electric motor. This additional traction power allows to downsize the engine and still fulfill the power requirements. Moreover, the engine can be turned off in idle phases for both architectures and for the parallel architecture, it may be turned off whereas the electric motor furnishes all the traction power. The optimal control problem of the energy management between the two power sources is solved for given driving cycles by a classical dynamic programming method and by an alternative method based on Pontryagin Minimum Principle. The real time control laws to be implemented on the vehicle are derived from the resulting optimal control strategies. These control laws are evaluated on another driving cycle which was not given a priori.

INTRODUCTION

Growing environmental concerns coupled with concerns about global crude oil supplies stimulate research on new vehicle technologies. Hybridelectric vehicles appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions [START_REF] German | Hybrid powered vehicles[END_REF] : mainly thanks to the system stop'n go that allows to turn off the engine in idle phases, to the recuperated braking energy to be stored in a battery and re-used later via the electric motor and to the possibility to downsize the engine.

The energy management of hybrid power trains requires then some specific control laws : they rely on the estimation of the battery state of charge which provides the remaining level of energy, and the variable efficiency of each element of the power train has to be taken into account. Optimization of energy management strategies on given driving cycles is often used to derive sub-optimal control laws to be implemented on the vehicle (see among others [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with onboard estimation of the control parameters[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF], [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], [START_REF] Delprat | Evaluation de stratégies de commande pour véhicules hybrides parallèles[END_REF]).

IFP, in partnership with Gaz de France and the Ademe, has combined its downsizing technology with a natural gas engine in a small urban demonstrator vehicle (VEHGAN vehicle), equipped with a starter alternator and supercapacitor manufactured by Valeo [START_REF] Tilagone | Development of natural gas demonstrator based on an urban vehicle with a down-sized turbocharged engine[END_REF].

In this paper, we present two different optimization algorithms and apply them to a simplified model of the VEHGAN vehicle and to a parallel architecture version of this vehicle: a classical Dynamic Programming algorithm ( [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF], [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with onboard estimation of the control parameters[END_REF]), and an original algorithm based on Pontryagin Minimum Principle that allows to handle constraints on the state and control variables. Finally, we propose two types of control strategies derived from the optimization results on given driving cycles and evaluate them as a real time strategy on a driving cycle which was not given a priori.

SYSTEM MODELLING AND OPTIMAL CONTROL PROBLEM

Characteristics of the considered hybrid vehicle

Two different architectures are modelled:

• a mild hybrid architecture : the engine can not be stopped when the requested torque is provided only by the electric motor, except for the stop'n go mode at the idle speed. So, for a control that cancels the engine torque and for positive torque request, the fuel consumption does not vanish (Figure 1), • a full parallel hybrid architecture : the engine can be stopped to let the electric motor power alone the vehicle. In that case, the fuel consumption vanishes.

In both cases, the battery is regenerated in braking phases accordingly to the available minimum electric torque at the considered engine speed.

In order to solve the optimal control problem of energy management, we build a simplified model which is composed of :

• a driving cycle to be followed (imposing vehicle speed and gear shifts), • a vehicle model defining its mass, wheel inertia, resistance force, • a manual gearbox with 5 gear ratios,

• a 660CC natural gas engine characterized by a fuel consumption map displayed in Figure 1 and a maximum torque depending on the engine speed (see ( 5)), • a starter alternator (3kW for mild-hybrid, 6kW for full-hybrid) characterized by a maximum torque and a minimum torque for regenerative braking phases, both depending on the engine speed (see ( 6)). Its efficiency is assumed to be 1 in the presented examples, • a battery characterized by a capacity of 0.4Ah for mild-hybrid architecture and 40Ah for full-hybrid one.

ẋ(t) = - ω(t)T m (t)K ′ U batt n capa (1)
with ω(t), the electric motor and engine speed (assumed to be equal), U batt , the battery voltage considered to be constant, K ′ , a scaling constant and n capa , the nominal capacity of the battery.

The driving cycle is converted in a (engine speed, torque) trajectory either thanks to a backward model based on the vehicle model, or thanks to a forward model as in AMESim Drive library which furnishes a more realistic trajectory taking into account a simulated behavior of a driver as the anticipation of the driving cycle.

Optimal Control Problem

The optimal control problem under study consists in minimizing the fuel consumption of the vehicle along a given driving vehicle cycle, taking into account physical constraints from battery, engine and electric motor. The control variable associated with this problem is called u(t). It represents the distribution of the requested torque T rq , between the engine torque T e and the electric motor torque T m , written as

   T rq (t) = T e (t) + T m (t) T e (t) = u(t)T rq (t) T m (t) = (1 -u(t))T rq (t).
(2)

The state variable is the battery state of charge x(t) and follows from ( 1)

ẋ(t) = -Kω(t)(1 -u(t))T rq (t) = f (u(t), t), (3)
where

K = K ′
U batt ncapa . The resulting optimization problem is then the following :

               min u    J(u) = T 0 L(u(t), t)dt + g(x(T ), T )    subject to : ẋ = f (u(t), t), x(0) = x 0 x min ≤ x(t) ≤ x max u min (t) ≤ u(t) ≤ u max (t) (4)
with 0 and T , respectively the initial and the final times of the given driving cycle, L(u(t), t), the instantaneous fuel consumption, computed from the map displayed in Figure 1, g(x(T ), T ), the penalization term that constrains the final state of charge to be close to the initial state of charge in order to maintain a null electrical energy balance (to avoid to discharge totally the battery for minimizing the consumption).

The bound constraints on the state and on the control in (4) are derived from the following constraints :

• the engine can only produce a positive torque, and is limited to a maximum torque which depends on engine speed ω(t), written as 0 ≤ T e (t) ≤ T max e (ω(t)), and leads to

0 ≤ u(t)T rq (t) ≤ T max e (ω(t)), (5)
• the electric motor torque is limited between a maximum torque and a minimum torque during regenerating breaking,

T min m (ω(t)) ≤ T m (t) ≤ T max m (ω(t))
, and leads to the control constraints

T min m (ω(t)) ≤ (1 -u(t))T rq (t) ≤ T max m (ω(t)),(6)
• the storage capacity implies a minimum and a maximum state of charge of the battery (which are fixed to 0% and 100% in our example)

x min ≤ x(t) ≤ x max . (7)
In this optimal control problem, we make several assumptions

• the pollutant emissions are not taken into account in the optimization process, • the engine speed and the electric motor speed are equal, • in the mild hybrid case, recharging the battery is only possible for negative torques (breaking request), we did not consider regeneration by an additional engine torque beyond the driver request torque. Thus the control u(t) remains between 0 and 1. In the full hybrid case, u(t) can take values larger than 1, allowing battery regeneration with additional engine torque.

In the following, we will call U (t) in continuous time (respectively U k in discrete time) the feasible domain for u(t) (respectively u k ) with respect to the constraints ( 5) and (6).

DYNAMIC PROGRAMMING OPTIMIZATION

The Dynamic Programming method (DP) is classically used to solve the problem (4) ( [START_REF] Wu | Optimization of power management strategies for a hydraulic hybrid medium truck[END_REF], [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF])) : it relies on the principle of optimality or Bellman principle. First, the optimal control problem (4) is discretized in time

         min u k ∈U k J(u) := N -1 k=0 L k (u k ) + g(x N )
subject to :

x k+1 = f k (x k , u k ), x(0) = x 0 x min ≤ x k ≤ x max (8)
where

L k (u k ) is the cumulated fuel consumption over the time interval [k, k + 1],
x k is the state of charge of the battery at time k, f k is the function that modelizes the battery state of charge evolution in the discrete form of (3) and g(x N ) = β.(x N -x 0 ) 2 is the penalization term for the constraint on final state of charge (β is a constant to be chosen1 ), N being the final time of the driving cycle.

From Bellman principle, the minimum cost

V k (x k ) at the time step k, 0 ≤ k ≤ N -1, is expressed as V k (x k ) = min u k ∈U k (L k (u k ) + V k+1 (f k (u k ))). (9) At time N , the cost function is V N (x N ) = g(x N ).
This optimization problem is solved backward from final time step to initial time step using a discretization of function V in the control space and in the state space.

DP Optimization algorithm

A standard time step used in our examples is 1s, and the step for state discretization is 0.5%. Two algorithms may be used to solve the DP problem :

• a classical DP algorithm, called Ford algorithm in the following [START_REF] Scordia | Approche systématique de l'optimisation du dimensionnement et de l'élaboration de lois de gestion d'énergie de véhicules hybrides[END_REF], consists in exploring all the feasible controls (to go from a point x i k to an other point x j k+1 ), finally taking the best trajectory (the trajectory which minimizes at each step k the sum

L k (u k ) + V k+1 (f k (u k ))
). In such a method, the state of charge trajectory remains on the points of the defined grid in the state space which may lead to inaccurate results.

• the chosen algorithm interpolates the function V (x k , k) in the state space, for each time step k thanks to an upwind scheme [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] :

V k (x i k ) = min u k ∈U k [∆tL k (u k ) + V k+1 (x i k+1 ) +f k (u k ) V k+1 (x i k+1 ) -V k+1 (x i-1 k+1 ) ∆x ∆t], (10)
where ∆x and ∆t are respectively the state and the time discretization step size. We refer to [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] for some theoretical results on the convergence of this method and error estimations. Therefore, it is possible to use a (state) continuous constrained optimization algorithm to solve each problem (9) which should furnish more accurate results than Ford algorithm. Nevertheless, this algorithm is generally more expensive in terms of computing time.

These two optimization algorithms are only used when T rq > 0 : when the requested torque is negative, the optimal control u k is completely known, as the battery is regenerated as much as possible, the control u k being constrained by the minimal electric motor torque from (6) and by maximum SOC from (7).

Optimization results obtained with DP method are displayed on Figure 2.

PONTRYAGIN MINIMUM PRINCIPLE OPTIMIZATION

In this section, we propose an alternative method to solve the optimal control problem (4). It relies on the Pontryagin Minimum Principle (PMP) and unlike the DP method does not require any discretization scheme.

Pontryagin Minimum Principle

First we consider the optimization problem (4) and introduce the Hamiltonian function, without considering state and control constraints H(u(t), x(t), p(t)) = L(u(t), t) + p(t) ẋ(t). ( 11) p(t) is called the co-state of our system. We assume here that L is a smooth convex function of u.

The Pontryagin Minimum Principle states the following conditions for the unconstrained optimal control problem :

∂H ∂x = -ṗ and ∂H ∂u = 0. ( 12
)
We refer to [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF] and [START_REF] Bryson | Applied Optimal Control[END_REF] for further details about Pontryagin Principle.

Application

The fuel consumption L(u(t), t) to be minimized in (4), is defined by a discrete map L(ω, T e ), mod-elled by a 2-order polynomial, which is represented as

L(ω, T e ) = 2 i,j=0 K ij ω i T j e , (13)
which allows to model a large variety of engine maps [START_REF] Rousseau | Constrained optimization of energy management for a mild-hybrid vehicle[END_REF].

Mild-Hybrid case

In the mild-hybrid vehicle case, the fuel consumption can not be cancelled. We do not consider the stop and start, as well as the possibility to power the vehicle only with the electric motor.

From ( 12) and ( 3) we obtain

ṗ = 0 ⇒ p = constant = p 0 . ( 14
)
Without any constraint on the state and on the control, the problem of minimizing H can be easily solved. The minimum fuel consumption is then reached for u * so as

∂H ∂u = ∂L ∂u + p ∂f ∂u = 0. ( 15
)
The optimal control u * can be calculated easily by solving the equation ( 15), which depends linearly on u (thanks to (3) and ( 13)) . u * finally depends on p(t), T rq (t) and ω(t)

u * (t) = - 2 i=0 K i1 ω(t) i + p 0 .K.ω(t) 2 2 i=0 K i2 ω(t) i .T rq (t) . ( 16
)
The expression of p 0 is obtained by replacing u * (t) by its expression in the state equation ( 3), and by integrating this equation in time, between T init and τ , T init and τ being respectively the considered initial and final times.

Full-Hybrid case

With the full-hybrid case, we have to consider the possibility to power the vehicle only with the electric motor. The previous expression of Hamiltonian becomes unadapted, as the fuel consumption can be completely cancelled. The fuel consumption function is then discontinuous

L f h (ω(t), T e(t)) = 0 if u(t) = 0 L(ω(t), T e(t)) if u(t) = 0. ( 17
)
The Hamiltonian, in the only electric motor case (u(t) = 0), is then written

H m (x(t), p(t)) = p(t) ẋ(t). ( 18
)
The optimal control u * must then be written as The previous section presents the computation of the optimal control of the continuous problem in a restricted case where no constraint is introduced. While control constraints are generally easily taken into account, handling the state constraints in the continuous optimal control problem is cumbersome: several singular cases can be found in [START_REF] Bryson | Applied Optimal Control[END_REF]. In our application, we are not able to find an analytic solution of the optimal control problem with control constraints : indeed, these constraints depends on time and depends on p 0 which depends on final SOC (cf. previous section). By an iterative method (called algo1 in the following), we can compute the value of p 0 in order to reach the desired SOC at final time with the control, expression ( 16), projected on its bound constraints. [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF], [START_REF] Pontryagin | Théorie mathématique des processus optimaux[END_REF], [START_REF] Evans | An Introduction To Mathematical Optimal Control Theory[END_REF], [START_REF] Bryson | Applied Optimal Control[END_REF], [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] have studied the general problem (4) with the state constraints. In our application, we can show that p(t) presents discontinuities at the time steps where the state inequality constraints are saturated. These time steps are not a priori known : this prevents us to solve explicitly the continuous optimal control problem with these state constraints.

u * =

PMP Optimization algorithm

Considering the difficulties described in previous section, we propose a heuristic iterative method that allows to find a sub-optimal trajectory from the constrained continuous optimal control problem (4). The proposed algorithm consists in an initialization step and 3 steps : (0) algo1 is applied on the driving cycle [0, T ] (see Figure 3 Step 0). The obtained optimal trajectory violates the state constraints, the farthest SOC (ie the "most violated point") from the bounds being for instance at point (x(t v ) = -37%, t v = 818s). The initial time is called t i , here set to 0.

(1) The SOC at t v is projected on the nearest bound of the feasible state domain (for instance, SOC is fixed to x min = 0 at point t v ).

(2) algo1 is applied again on [t i , t v ] (see Figure 3 Step 2). If the obtained trajectory still violates the state constraints on [t i , t v ], steps 1 and 2 are applied again on the farthest SOC from the bounds (defining a new point t v ). This procedure is repeated until the trajectory remains on the feasible domain. Then the last point t v becomes the new initial time t i in step 3. (3) algo1 is applied on [t i , T ] (see Figure 3 Step 3). If the obtained optimal trajectory still violates the state constraints, steps 1 and 2 are repeated. This sequence is repeated until we reach the final step T at the desired final SOC, without violating the state constraints (Figure 3 bottom right).

Some optimization results

Mild Hybrid case

We can compare the two optimization algorithms (DP and PMP) on the Urban Artemis driving cycle [START_REF] André | The artemis european driving cycles for measuring car pollutant emissions[END_REF], in the mild Hybrid case, on Figure 2. The curves are very similar; we can notice that smaller is the state step size, nearer to the PMP curve are the DP curves.

Figure 4 presents the operating points (OP) of the engine obtained with PMP algorithm.

In this vehicle configuration, the state constraints are active 5 times, giving 6 different values of the Lagrange multiplier p(t). We display the six curves (green lines) ∂H ∂Te (p) = 0, which give optimal engine torque, function of engine speed. The engine OP are thus moved toward the green optimal curves when it is possible: the OP located below the curves remain unchanged (no battery regeneration being possible for positive torque requests for mild hybrid) whereas the OP located above are moved toward the curves by decreasing the engine torque as much as possible (saturating electric motor torque constraints). 4.3.2. Full Hybrid case Figure 5 gives optimized operating points for the engine and the electric motor (PMP algorithm is used). In addition to kinetic energy, we assume that it is possible to recharge the battery by using the engine at better OP, with an ideal efficiency of 1.

As for mild-hybrid case, the optimal trajectory (continuous green line) gives the optimal operating points of the engine by finding the solution of ∂H ∂Te = 0. Thus, many of low torque OP are moved to the optimal trajectory, recharging the battery by imposing a negative electric motor torque. As the full-hybrid configuration allows to turn off the engine for non-zero vehicle speed (pure electric mode), most of OP associated with engine speed below 3000 rpm and requested torque below 20Nm, lead to turn off the engine (points where engine torque is zero) : turning off the engine is more efficient than the optimal engine torque (green curve : ∂H ∂Te = 0).

REAL-TIME CONTROL

From optimization results on Urban Artemis cycle, we derive suboptimal control laws that will be tested on an other cycle. In this section, the FTP72 cycle has been chosen, for its realism of urban driving. For the full hybrid architecture, the two control laws give degraded results compared to optimal results. Many reasons can explain these differences. First, even if Urban Artemis cycle and FTP72 cycle are both realistic of an urban driving, operating points are very different. While requested operating points of Artemis cycle are almost uniformly located in the whole engine speed and torque space, all requested operating points of FTP72 are below ω = 3200 rpm, with a majority below ω = 2000 rpm. The consequence is a unadapted electric motor map for the second control law. Concerning the first control law, the optimal p (obtained with PMP algorithm on FTP72) is quite different from the optimal p obtained for Artemis cycle, leading to degraded results. Nevertheless, the consumption gain remains high : -14.76%.

These results illustrate that several driving cycles are needed to develop efficient suboptimal control laws based on p-control or electric motor map.

The vehicle speed (related to engine speed by gear ratios) could also be taken into account to improve fuel consumption gains.

CONCLUSIONS

In this study, we have presented two methods for optimal control optimization. The heuristic method based on Pontryagin Minimum Principle, well known in the free state constraint case, has been applied successfully to our state constrained problem, with very similar results to Dynamic Programming methods and a computation time divided by 100. Nevertheless, there is currently no theoretical proof to confirm the presented validation results. Moreover, there are some limitations to this approach, mainly the assumptions on the fuel consumption map, modelled by a smooth convex function of control u (2-order polynomial) ; this limitation could lead to a bad approximation of the real fuel consumption for some particular engines.

Other degrees of freedom, as the gear-shifting sequence should also be taken into account in the optimization problem to improve the fuel consumption gain. Reduction of pollutant emissions will also be studied by considering a second state based on exhaust temperature. From optimization results are derived two types of suboptimal feedback laws based on state of charge measurements. These laws give encouraging results even if it needs to be improved in the full hybrid case.

Fig. 1 .

 1 Fig. 1. Fuel consumption map of natural gas engine of VEHGAN vehicle

  Fig. 2. Urban Artemis cycle (Top); Optimal state of charge trajectory of VEHGAN vehicle computed with PMP & DP algorithm (Bottom).

Fig. 3 .

 3 Fig. 3. The proposed algorithm based on Pontryagin Minimum Principle.

Fig. 5 .

 5 Fig. 5. Operating points of engine in Full-hybrid mode obtained by PMP algorithm for the urban Artemis Driving Cycle.

Table 1 .

 1 Two different control laws will be tested : the first one, based on Optimization results from Pontryagin principle, consists of varying the value of p regarding to the state of charge, to control u(t), then the electric motor. The reference Lagrange multiplier value p is the mean of optimal values of p, obtained on Artemis Urban cycle with off-line optimization using PMP algorithm. The second one uses a map of electric motor torque created by the optimization results on Urban Artemis cycle. The electric motor torque from the map is then weighted by the state of charge of the battery : reduced if the SOC is low, increased if the SOC is high. The obtained results are displayed in Table1. For the mild hybrid configuration, the suboptimal laws give fuel consumptions which are close to the optimal one. Fuel Consumption

	Consump.	Th.	Optimal	p-control	Elec. mot.
	(l/100km)	veh. control	based	torq. map
	Mild-H.	3.32	3.22	3.23	3.23
			(-3,01%)	(-2,71%)	(-2,71%)
	Mild-H with		2.86	2.87	2.88
	Stop'n go.		(-13,62%)	(-13,49%)	(-13,33%)
	Full-H.		2.70	2.83	2.86
			(-18.67%)	(-14,76%)	(-13,85%)

In the following results, a value depending of battery capacity has been implemented

ACKNOWLEDGMENTS

We would like to thank Gilles Corde, Philippe Moulin and Antonio Sciarretta for helpful discussions and advice at various stages of the elaboration of this work. We acknowledge Quang Huy Tran for his advice on numerical methods.