Gregory Rousseau
email: gregory.rousseau@ifp.fr

Antonio Sciarretta

Delphine Sinoquet

REAL-TIME CONTROL STRATEGIES FOR HYBRID VEHICLES ISSUED FROM OPTIMIZATION ALGORITHM

Keywords: Hybrid vehicle, Optimal control, Dynamic programming

This paper focuses on a mild-hybrid city car (Smart), equipped with a starter-alternator, where the kinetic energy in the braking phases can be recovered to be stored in a supercapacitor, and re-used later via the electric motor. The 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 optimal control problem of the energy management between the two power sources is solved for given driving cycles by a classical dynamic programming method. From dynamic models of the electric motor and supercapacitor a quasistatic model of the whole system is derived and used in the optimization. The real time control law to be implemented on the vehicle is derived from the resulting optimal control strategies.

INTRODUCTION

Growing environmental concerns coupled with concerns about global crude oil supplies stimulate research on new vehicle technologies. Hybrid-Electric Vehicles (HEV) appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions. HEVs can save fuel thanks to: (i) the possibility of turning off the engine in idle phases (stop-and-go), (ii) the recuperation of braking energy to be stored in the battery and re-used later via the electric motor, (iii) the possibility of downsizing the engine, and (iv) the optimization of energy flows.

The energy management of hybrid powertrains requires some specific supervisory control laws. This controller relies on the estimation of the battery state of charge, and it must take into account the variable efficiency of each element of the powertrain. Optimization of energy-management strategies on prescribed driving cycles (offline optimization) is often used to derive sub-optimal control laws to be implemented on the vehicle [START_REF] Sciarretta | 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].

Optimization tools require a mathematical description of the system. In most offline optimization techniques, the quasistatic backward modelling approach is used. On the other hand, dynamic forward-facing models are necessary to accurately simulate the system and to identify the main system parameters.

In this paper, we present a dynamic model of the electric propulsion system of a mild hybrid vehicle. The quasistatic version of the model is used in a Dynamic Programming algorithm to optimize the energy management of the vehicle. Finally, we propose a real-time feedback control strategy derived from the offline optimization results, which is validated for typical driving conditions.

SYSTEM MODELLING

IFP, in partnership with Gaz de France and ADEME, has applied its downsizing technology to a natural gas engine, resulting in a small urban demonstrator vehicle (VEHGAN vehicle). The conventional powertrain is integrated with a belt-driven starter-alternator and a supercapacitor, both provided by Valeo [START_REF] Tilagone | Development of natural gas demonstrator based on an urban vehicle with a down-sized turbocharged engine[END_REF].

The resulting architecture is a mild-hybrid powertrain: the torque requested for traction cannot be provided by the electric motor only. Thus, the engine cannot be turned off, except during idling (stop-and-go).

In order to optimize the energy management of the VEHGAN, we use either a quasistatic backward model, or a more detailed AMESim model [START_REF] Dabadie | A new tool for advanced vehicle simulations[END_REF] to calculate the speed ω rq (t) and torque T rq (t) required at the engine output shaft along a prescribed drive cycle. This torque request can be split between the engine torque T e (t) and the electric motor torque T m (t). Finding the optimal splitting strategy is the objective of the supervisory controller. To evaluate the performance of the two components as a function of their respective torque, we build a simplified model of the system, which consists of:

• a model of the 0.660 l natural gas engine; the engine performance is characterized by the fuel consumption map displayed in Figure 1 and by its maximum torque T max e (ω rq ) that depends on the engine speed (see ( 18)),

• a model of the 3 kW starter alternator (described in section 2.1.1), • a model of the supercapacitor, which is made up of 12 modules of 1500 F; this model is described in section 2.1.2.

Submodels of the electric components

Electric motor

The starter-alternator used in the VEHGAN vehicle is a permanentmagnet motor with six magnets. A realistic model of torque evolution, as described in [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF], is given by

T m (t) = 3 2 pϕ m . . R st U q (t) -pω(t)R st ϕ m -pω(t)L st U d (t) R 2 st + p 2 ω 2 (t)L 2 st , (1)
where ω(t) is the rotor speed, R st the stator resistance, ϕ m (t) the mutual flux linkage, U q (t) and U d (t) are the d-q axis stator voltage components, L st is the stator inductance, and p is the number of pole pairs.

For simplification convenience, we choose to model this motor using DC motor equations by neglecting the non-linear contribution. We then replace (1) by

T m (t) = 1 R a (K a (t)U a (t) -K a (t)K m (t)ω(t)), (2)
where ω(t) is the motor speed, U a (t) plays the role of the DC armature voltage,

U a (t) = K m (t)ω(t) + R a I a (t), (3)
I a (t) plays the role of the armature current, R a of the armature resistance.

The quantities K a (t) and K m (t) are not constant as in DC motors, since they represent influence of the variable d,q voltages and currents performed by the motor controller. The equivalence between real AC quantities and fictitious DC counterparts U a (t), I a (t), K a (t), K m (t), and R a (t) is found by comparing ( 1) and (2-3).

Supercapacitor

A simple equivalent circuit of a supercapacitor consists of a capacitor and a resistor in series [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with on-board estimation of the control parameters[END_REF]. The Kirchhoff's voltage law yields

U s (t) = Q(t) C -R s I s (t), I s (t) = - dQ(t) dt , (4)
where U s (t) is the terminal voltage, I s (t) is the terminal current, R s is the equivalent resistance, C is the capacitance, and Q(t) is the charge.

2.1.3. DC/AC Link The electric power link between electric motor power P a and supercapacitor power P s is achieved with an inverter. We assume here that the motor torque is controlled by the voltage ratio λ(t) of the inverter:

U s (t) = 1 λ(t) U a (t), I s (t) = λ(t)I a (t). (5)
Because of internal physical constraints, it is not possible to control the electric motor with a continuous control variable λ(t) when positive torque is requested. Instead, five discrete values of λ are available at each time t, as a function of the system operating conditions. This discrete controller permits to obtain five values of torque, depending on the supercapacitor voltage and the motor speed. The difference between minimum and maximum available torque does not exceed 2 Nm.

Dynamic model of the whole electric system

The physical causality representation of the whole electric system is sketched in Figure 2. The speed is an input variable, together with the controller λ. The torque is the output variable. To obtain the dynamic dependency between torque and speed, we differentiate each term of (4), which leads to

Us (t) = - I s (t) C -R s İs (t). (6)
From ( 3) and ( 5) we obtain

I s (t) = λ(t) λ(t)U s (t) -K a (t)ω(t) R a . (7)
From ( 6) and ( 7) we finally obtain

Us (t) = λ(t) 1 + λ(t) 2 Rs Ra - λ(t)U s (t) -K a (t)ω(t) CR a + R s R a K a (t) ω(t) . (8)
The model given by ( 8) is validated against measurements taken on a dedicated test bench at Valeo. Figure 3 shows the transient behavior of the supercapacitor voltage resulting from a motor speed step occurring in the first seconds of test.

The experimental data are post-processed with a filter having a time constant of about 100 ms. Four speed levels are shown in the figure, ranging from 3000 rpm to 6000 rpm. Starting from an initial open-circuit value of 24 V, the supercapacitor voltage exhibits an abrupt drop of about 4 V, followed by a slower decrease. Consequently, the supercapacitor current (Figure 4) exhibits an initial step appearing as a smooth peak due to the low-pass filtering, followed by a slower decrease that would virtually end when the supercapacitor is completely discharged.

The comparison between experimental data and simulations of Figure 3 and 4 shows that the simplified approach of ( 8) is indeed able to represent the main dynamics of the system. The control parameters K a and K m have been parameterized as

K a = K a0 ω(t) , K m = K m0 U s (t) ω(t) . ( 9
)
The fitted values of the model parameters are listed in Table 1. The tests of Figures 3 and4 were performed for a given setting of the discrete motor controller.

In the simulations, that corresponds to a specific value λ 4 of the voltage ratio of the inverter. The torque curves calculated for all of the controller settings are shown in Figure 5. For the specific value λ 4 , the simulation results can be compared to experimental data.

Quasistatic model

A quasistatic counterpart of the model presented above is developed for the offline optimization.

The causality representation of the electric system is sketched in figure 6. From the motor speed ω(t) and torque T m (t) required, provided they are compatible with the physical limits of the motor, the armature voltage and current are computed with (2) and (3). The requested electric power is computed as

P (t) = U a (t)I a (t), (10)
and, assuming no conversion losses as in (5), is also

P (t) = U s (t)I s (t).
From ( 4) and ( 10) we obtain a differential equation for U s (t) as a function of P (t),

(1 -

R s P (t) U 2 s (t) ) dU 2 s (t) dt = -2 P (t) C -2R s Ṗ (t). ( 11
)
Using numerical integration methods, this equation is integrated at any time using the value of U s at the previous time step

U 2 k+1 = - 2P k U 2 k ∆t C(U 2 k -R s P k ) + U 2 k , (12)
with:

• ∆t: constant time step,

• U k ,U k+1 : supercapacitor voltage at time step k and k + 1, • U k : average voltage during time interval [k, k + 1],

• P k : electric power during time interval [k, k + 1].

In (12), we have neglected the term proportional to Ṗ of (11) since R s C ≪ ∆t.

The consequent value of λ to obtain the requested torque T m is calculated using

I sk = P k /U k , (13)
then, with (5):

λ k = I sk /I ak . (14)

DYNAMIC PROGRAMMING OPTIMIZATION

Optimal Control Problem

The optimal control problem under study consists in minimizing the fuel consumption of the vehicle along a prescribed driving cycle, while providing charge sustenance to the supercapacitor, and taking into account physical limits of the supercapacitor, the engine, and the electric motor.

We define the control variable u(t) which represents the split factor of the requested torque T rq , between the engine torque T e and the motor torque T m :

T rq (t) = T e (t) + ρT m (t) T e (t) = u(t)T rq (t) , (15)
where ρ corresponds to the ratio between engine speed and electric motor speed.

The only relevant dynamics considered is that of supercapacitor voltage, called U (t) in the following, which follows from ( 11)

U (t) = f (U (t), u(t), t). (16)
The resulting optimization problem is then the following :

       min u∈W (t)    J(u) = T 0 L(u(t), t)dt + g(U (T ), T )    subject to : U (t) = f (U (t), u(t), t), U (0) = U 0 (17)
where T is the final time of the drive cycle, L(u(t), t) is the fuel consumption rate, computed from the data displayed in Figure 1,g(U (T ),T ) is a term that penalizes any deviation of the final voltage from the initial voltage.

In addition to the latter "soft" constraint over the final state, the problem ( 17) is also subjected to the following inequality constraints, leading to the feasible control space W (t) :

• the engine can only produce a positive torque, which is limited to a maximum torque that depends on engine speed ω rq (t), leading to the control constraints 0 ≤ u(t)T rq (t) ≤ T max e (ω rq (t)), (18)

• the electric motor torque is limited between a maximum torque and a minimum torque during generating operation, leading to the control constraints

ρT min m (ω(t)) ≤ (1 -u(t))T rq (t) ≤ ρT max m (ω(t)), ( 19
)
• the motor controller can only be set to one over five discrete values, thus

λ(t) ∈ (λ i ) i∈ 0,4 , (20)
as already discussed; that also leads to corresponding control constraints.

The state variable is also subjected to inequality constraints, deriving from physical limits of the supercapacitor and of the motor,

U min ≤ U (t) ≤ U max . ( 21
)
The values of U min and U max are set to 16 V and 24 V respectively.

DP Optimization algorithm

The Dynamic Programming method (DP) is classically used to solve the problem ( 17) [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 ( 17) is discretized in time

         min u k ∈W k J(u) := N -1 k=0 L k (u k ) + g(U N ) subject to : U k+1 = f k (U k , u k ), U (0) = U 0 U min ≤ U k ≤ U max (22)
where L k (u k ) is the fuel consumption over the time interval [k, k + 1], U k is the voltage of the supercapacitor at time k, f k is the function ( 16) at time k, and g(U N ) = β(U N -U 0 ) 2 is the penalization term (β is a constant to be chosen), N is the final time, W k the feasible domain for u k with respect to constraints ( 18), ( 19), and (20).

From Bellman principle, the cost-to-go function

V k (x k ) at the time step k, 0 ≤ k ≤ N -1, is expressed as V k (U k ) = min u k ∈W k (L k (u k ) + V k+1 (f k (U k , u k ))).(23) At time N , the cost-to-go 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.

Model implementation

As to control the electric system, both quasistatic and dynamic models are used. First, from a requested torque T rq , the quasistatic model ( 10)-( 14) is used to calculate the corresponding necessary value of λ, namely λ rq . Then, only the control set u adm that corresponds to the set λ adm as

λ adm = {λ |λ ∈ (λ i ) i∈ 0,4 , λ ≤ λ rq } (24) is tested.
A standard time step used in our examples is 1 s. The solving algorithm used is that of [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] to which we refer for some theoretical results on the convergence and error estimations.

Optimization Results

The optimal system trajectories are calculated for the Urban Artemis cycle, shown in Figure 7.

The figure also shows the optimal trajectory of motor torque, which exhibits a somehow intuitive behavior. The motor is used to assist the engine during accelerations, and as a generator during The optimal voltage trajectory of the supercapacitor is presented in Figure 8. The figure clearly shows that, since the efficiency of the electric powertrain is higher at high voltage, the electric motor is used mostly between 20 V and 24 V.

The fuel consumption gain between pure engine mode and optimal hybrid mode is about 11%. Clearly, the system is not designed to achieve substantial improvements in fuel economy with respect to the conventional architecture, but rather to improve its performance, e.g., in terms of acceleration responses at low speed. In this section, suboptimal control law is derived from the optimization results.

Figure 9 shows the optimal values of the electric motor torque as a function of requested torque and speed. When negative torques are required, i.e., during regenerative braking, the whole torque available is used to recharge the supercapacitor.

The factor "2" between negative electric motor torque and negative requested torque comes from the speed ratio between engine and electric motor (ρ = 1.96). For positive requested torque, the following behavior can be observed:

• the electric machine is only used as a motor above 2200 rpm (engine speed). Below this speed, the electric machine can be used to recharge the supercapacitor, • the electric machine is only used as a motor when the requested torque is higher than 20 Nm, while it is used in the recharging mode in all the requested torque range, with a predominance of low requested torques (0 to 20 Nm).

An intuitive controller can be derived from these considerations, based on operating domains as a function of speed and torque requests, as shown in Figure 10.

In the boosting operating domain (positive motor torque area), the motor torque is subjected to the constraint (20). Figure 11 shows the optimal values of the discrete controller setting λ(t) along the Artemis Cycle. Clearly, the optimal choice consists in selecting always the minimum admissible value λ 0 , which corresponds to the minimum admissible motor torque. It can be proven that this choice maximizes the efficiency of the electric subsystem.

In the regenerating operating domain (negative motor torque area) of Figure 10, the controller setting λ can be varied continuously, meaning that the recharging torque can have every value below the physical limit -T min m . Searching for a simple heuristic rule to represent the optimal motor behavior, it seems rather intuitive to compare the optimal regenerating torque of Figure 9,T • m with the torque values that maximize the global efficiency of the system, Tm . When the engine recharges the supercapacitor through the electric machine, the global efficiency is calculated as

η charge = 1 2 CU s ∆U s H i ∆L∆t , (25)
where ∆U s the variation of the supercapacitor voltage for the time step ∆t and ∆L is the additional fuel consumption rate that is used to recharge the supercapacitor.

A comparison between T • m and Tm is shown in Figure 12 for the Artemis cycle. Clearly, the two data sets are very different. Consequently, the optimal motor behavior during regeneration phases cannot be simply described by a rule based on (25). The reason is that the optimal motor torque when recharging the supercapacitor does not depend on local conditions only, but it is a function of the optimal motor torque during boosting phases, as well as of the state boundaries and final state constraint. Noncausal solutions of the optimal control problem (17) may lead to approximate, suboptimal controllers that perform better than heuristic feedback controllers, like the one described in this section, in representing the optimal behavior of the system.

CONCLUSIONS

In this paper, we have proposed a dynamic model of the electric system used in a mild hybrid vehicle. The model is based on physical equations of a DC electric motor, a supercapacitor, and an inverter. The model parameters have been successfully fitted to experimental data in order to simulate accurately the behavior of the electric system.

This dynamic model has been used in a Dynamic Programming algorithm to minimize the fuel consumption, taking into account the main physical constraints of the system (especially the discrete control settings). The optimal trajectory leads to an interesting gain in fuel consumption (9%) despite the reduced size of the electric motor and capacity of the supercapacitor.

From the optimization results on the Artemis cycle operating domains for the different electric motor modes (assist or recharge) were derived. For each domain a heuristic control law adapted to the electric motor mode was defined. Although the heuristics is able to represent a good approximation of the optimal switching between operating modes, it appeared to be a too simplified approach to calculate the optimal motor torque trajectories.

Future work will therefore include the synthesis of a noncausal controller based on an appropriate, although approximate, solution of the optimal control problem.

Fig. 1 .

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

Fig. 2 .

 2 Fig. 2. Dynamic representation of the electric system

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 6. Quasistatic representation of the electric system

Fig. 7 .

 7 Fig. 7. Vehicle speed on Artemis Cycle / Optimal trajectory of motor torque decelerations and during low-speed, low-torque phases.

Fig. 8 .

 8 Fig. 8. Optimal voltage trajectory on Artemis cycle (minimal and maximal bounds on voltage: 16V and 24V)

Fig. 9 .

 9 Fig. 9. Optimal electric motor torque T • m as a function of torque and speed requests, calculated with DP for the Artemis cycle.

Fig. 11 .

 11 Fig. 10. Operating domains of the electric motor as a function of torque and speed requests

  Fig. 12. Relation between calculated electric torque

Table 1 .

 1 Numerical values of the fitted model parameters

	Parameter	Value
	C	196 C
	Ra	30.7 Ω
	Rs	0.6 mΩ
	Ka 0	8
	Km 0	125

ACKNOWLEDGMENTS

We would like to thank Paolo Tona and Pierre Rouchon 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. We thank Sylvain Delion from Valeo for providing data and for his advice about the model.


Gregory Rousseau
email: gregory.rousseau@ifp.fr

Antonio Sciarretta

Delphine Sinoquet

OPTIMAL ENERGY MANAGEMENT OF A MILD-HYBRID VEHICLE

The paper presents the development of a supervisory controller for a mild-hybrid vehicle, a hybrid natural gas SMART, equipped with a starter alternator and supercapacitor manufactured by Valeo. This electric additional power can be used to stop and start quickly the engine and also to power the vehicle alongside with the engine. The electric motor can also be used to recharge the supercapacitor. After a description of the models developed for the electric motor dynamics, a dynamic programming algorithm is applied for the optimization of power split, based on these models. The resulting optimal power split is compared to a real-time control law. Among the available control laws, the choice of the Equivalent Consumption Minimization Strategy (ECMS) allows to keep the same models that have been used for dynamic programming algorithm. Moreover, some road tests show the resulting behavior of the powertrain, in terms of supercapacitor voltage, motor and engine torque and speed.

INTRODUCTION

Growing environmental and global crude oil supplies concerns are stimulating research on new vehicle technologies. Hybridelectric vehicles (HEV) appear to be one of the most promising technologies for reducing fuel consumption. In this context, mild hybridization gives numerous potential advantage : (i) (stop-and-start) improves driving agreement by reducing noise and vibrations, (ii) the recuperation of braking energy to be stored in the battery and re-used later via the electric motor, (iii) the possibility of downsizing the engine, and (iv) the ability to choose high-efficiency engine running modes to reduce energy consumption.

The energy management of hybrid powertrains requires dedicated supervisory control laws. The controller must actually rely on the estimation of the state of charge of the energy storage system, and it must take into account the variable efficiency of each element of the powertrain. Optimization of energy-management strategies on prescribed driving cycles (offline optimization) is often used to derive sub-optimal control laws to be implemented on the vehicle [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].

Optimization tools require a mathematical description of the components of the system. In this paper, a model-based control of the propulsion system of a mild hybrid vehicle is presented. This model is used in dynamic programming algorithm, and also with a real-time ECMS-type control strategy. With respect to other implementations of this concept [START_REF] Delprat | Evaluation de strat égies de commande pour v éhicules hybrides parall èles[END_REF], [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with onboard estimation of the control parameters[END_REF]) usually adopted for battery-based parallel hybrids, here the studied powertrain includes a supercapacitor as an energy storage system. Due to its small capacity, in the normal situation, the supercapacitor is often close to its physical voltage limits. That is a major difficulty when using real-time control, which deserves particular techniques detailed in the following.

SYSTEM MODELLING

IFP has developed a mild hybrid vehicle, equipped with a natural gas engine, a reversible starter-alternator, and a supercapacitor (Valeo Starter Alternator Reversible System, StARS). [START_REF] Tilagone | Development of natural gas demonstrator based on an urban vehicle with a downsized turbocharged engine[END_REF]. The partners of this project are Gaz de France, INRETS and Valeo, with financial support by ADEME.

The resulting architecture is a mild-hybrid powertrain: the torque requested for traction cannot be provided by the electric motor only. Thus, the engine cannot be turned off, except during idling (stop-and-start). The electric motor can be controlled by several modes:

(1) the boost mode, where the electric motor provides a traction power alongside with the engine; (2) the regenerative mode, where a the motor is used as a generator to decrease vehicle speed and recharge the supercapacitor (regenerative braking), or b the generator drains an extra power from the engine (torque defined by energy management); (3) the alternator mode, where a constant supercapacitor voltage setpoint is imposed (typically 16.2 V); this mode is only used when the voltage drops below a certain threshold.

In the following, the electric power required by the auxiliaries (supply of electric power for lights, vehicle control unit, etc.), and drained from a 14 V auxiliary battery is taken into account as well. However, for simplicity, we assume that it is the supercapacitor that provides the electric power (via the DC/DC module) to maintain the battery voltage at 14 V. We also assume that the requested power is constant at 140 W.

In order to optimize the energy management of the VEHGAN, either a quasistatic backward model is used, or a more detailed AMESim model [START_REF] Dabadie | A new tool for advanced vehicle simulations[END_REF] to calculate the speed ω rq (t) and torque T rq (t) required at the engine output shaft along a prescribed drive cycle. This torque request can be split between the engine torque T e (t) and the electric motor torque T m (t). Finding the optimal splitting strategy is the objective of the supervisory controller. To evaluate the performance of the two components as a function of their respective torque, a simplified model of the system is built, which consists of:

• a model of the 0.660 l natural gas engine; the engine performance is characterized by the fuel consumption map displayed in Figure 1 and by its maximum torque

T max e (ω rq ) that depends on the engine speed, see ( 13), • a model of the 3 kW starter-alternator (described in section 2.1.1),

• a model of the supercapacitor, which is made up of 12 modules of 1500 F; this model is described in section 2.1.2.

S ubmodels of th e electric components

All the models described in the following section are taken from [START_REF] Rousseau | Realtime control srategies for hybrid vehicles issued from optimization algorithm[END_REF]. To model the electric motor in a realistic way, two different models depending on the sign of the power are considered. If the delivered power is positive, then the electric motor is used to provide some additional power to the engine, and the torque is positive. Otherwise, when the electric motor is used as a generator, maps provided by the manufacturer are used, which give the delivered current as a function of voltage, speed, and torque.

Electric motor

The starter-alternator used in the VEHGAN vehicle is a permanent-magnet motor with six magnets.

To limit the complexity of the model, the motor is modelized using DC motor equations by neglecting the non-linear contribution:

T m (t) = 1 R a (K m (t)U a (t) -K a (t)K m (t)ω(t)), (1)
where ω(t) is the motor speed, U a (t) plays the role of the DC armature voltage,

U a (t) = K a (t)ω(t) + R a I a (t), (2)
I a (t) plays the role of the armature current, R a of the armature resistance.

The quantities K a (t) and K m (t) are not constant as in DC motors, since they represent influence of the variable d,q voltages and currents performed by the motor controller.

S upercapacitor

A simple equivalent circuit of a supercapacitor consists of a capacitor and a resistor in series [START_REF] Sciarretta | A real-time optimal control strategy for parallel hybrid vehicles with onboard estimation of the control parameters[END_REF]. The K irchhoff's voltage law yields

U s (t) = Q(t) C -R s I s (t), I s (t) = - dQ(t) dt , (3)
where U s (t) is the terminal voltage, I s (t) is the terminal current, R s is the equivalent resistance, C is the capacitance, and Q(t) is the charge.

. D C/AC L ink

The electric power link between electric motor power P a and supercapacitor power P s is achieved with an inverter. We assume here that the motor torque is controlled by the voltage ratio λ(t) of the inverter:

U s (t) = 1 λ(t) U a (t), I s (t) = λ(t)I a (t). (4)
B ecause of internal physical constraints, it is not possible to control the electric motor with a continuous control variable λ(t) when positive torque is requested (boost mode). Instead, five discrete values of λ are available at each time t, as a function of the system operating conditions. This discrete controller permits to obtain five values of boost torque, depending on the supercapacitor voltage and the motor speed. The difference between minimum and maximum available torques does not exceed 2 Nm. In regenerative mode, however, it is possible to control the electric motor with a continuous λ, meaning that all admissible regenerative torques are available.

D y namic model of th e w h ole electric sy stem

The physical causality representation of the whole electric system is sketched in Figure 2. The speed is an input variable, together with the controller λ. The torque is the output variable.

ω

T m I a U a motor DC/AC conv'r λ supercap I s U s Fig. 2

. Dynamic representation of the electric system

To obtain the dynamic dependency between torque and speed, we differentiate each term of (3), which leads to

Us (t) = - I s (t) C -R s İs (t).
(5)

From ( 2) and ( 4) we obtain

I s (t) = λ(t) λ(t)U s (t) -K a (t)ω(t) R a . (6)
From ( 5) and ( 6) we finally obtain

Us (t) = λ(t) 1 + λ(t) 2 R s R a - λ(t)U s (t) -K a (t)ω(t) CR a + R s R a K a (t) ω(t) . (7)
The model given by ( 7) is validated against boost measurements (see [START_REF] Rousseau | Realtime control srategies for hybrid vehicles issued from optimization algorithm[END_REF]) taken on a dedicated test bench at Valeo. Figure 3 shows the transient behavior of the supercapacitor voltage resulting from a motor speed step occurring in the first seconds of test. The experimental data are post-processed with a filter having a time constant of about 100 ms. Four speed levels are shown in the figure, ranging from 3000 rpm to 6000 rpm. Starting from an initial open-circuit value of 24 V, the supercapacitor voltage exhibits an abrupt drop of about 4 V, followed by a slower decrease. Consequently, the supercapacitor current (not represented) exhibits an initial rise, appearing as a smooth peak due to the low-pass filtering, followed by a slower decrease that would virtually end when the supercapacitor is completely discharged.

The comparison between experimental data and simulations of Figure 3 shows that the simplified approach of ( 7) is indeed able to represent the main dynamics of the system. The control parameters K a and K m have been parameterized as

K a = K a0 ω(t) , K m = K m0 U s (t) ω(t) . (8 )
The tests of Figure 3 were performed for a given setting of the discrete motor controller. In the simulations, that corresponds to a specific value λ 4 of the voltage ratio of the inverter. For the specific value λ 4 , the simulation results can be compared to experimental data.

To model the regenerative mode behavior, instead of using ( 7), we replace the DC/AC model and the electric motor model by a map provided by Valeo. Equation ( 6) is then replaced by the following

I s (t) = f (U s (t), ω(t), T m (t), (9)
where f (•) corresponds to the map depending on speed ω, voltage U s , and torque T m . The motor torque is calculated as matching exactly the torque setpoint T • m , which thus substitutes λ as an input variable to the motor/converter submodel, see

Figure 2.

The overall model has been implemented in Matlab Simulink, and used within a dynamic programming algorithm.

DYNAMIC PROGRAMMING OPTIMIZ ATION

.1 O ptimal Control Problem

The optimal control problem under study consists in minimizing the fuel consumption of the vehicle along a prescribed driving cycle, while providing charge sustenance to the supercapacitor, and taking into account physical limits of the supercapacitor, the engine, and the electric motor.

We define the control variable u(t) which represents the split factor of the requested torque T rq , between the engine torque T e and the motor torque T m :

T rq (t) = T e (t) + ρT m (t) T e (t) = u(t)T rq (t) , (10)
where ρ corresponds to the ratio between electric motor speed and engine speed, ρ = ω(t)/ω rq . This ratio is equal to 1.9 6

on the studied vehicle.

The only relevant dynamics considered is that of supercapacitor voltage, called U(t) in the following

U(t) = f (U(t), u(t),t). ( 11
)
The resulting optimization problem is then the following : In addition to the latter "soft" constraint over the final state, the problem ( 12) is also subject to the following inequality constraints, leading to the feasible control space W (t):

       m in u∈W (t)    J(u) = T 0 L(u(t),t)dt + g(U(T ), T )    s u b je c t to : U(t) = f (U(t), u(t),t), U(0) = U 0 (12)
• the engine can only produce a positive torque, which is limited to a maximum torque that depends on engine speed ω rq (t), leading to the control constraints

0 ≤ u(t)T rq (t) ≤ T max e (ω rq (t)), (13)
• the electric motor torque is limited between a maximum torque and a minimum torque during regenerative mode, leading to the control constraints

ρT min m (ω(t)) ≤ (1 -u(t))T rq (t) ≤ ρT max m (ω(t)), (14)
• the motor controller, during boost mode, can only be set to one over five discrete values, thus

λ(t) ∈ (λ i ) i∈ 0,4 , (15)
as already discussed; that also leads to corresponding control constraints on the allowed positive torques.

The state variable is also subject to inequality constraints, deriving from physical limits of the supercapacitor and of the motor,

U m in ≤ U(t) ≤ U m a x . (16)
The values of U m in and U m a x are set to 16 V and 24 V respectively.

.2 D P O ptimiz ation alg orith m

The Dynamic Programming method (DP) is usually used to solve the problem (12) ( [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 B ellman principle. First, the optimal control problem ( 12) is discretized in time

       m in u k ∈W k J(u) := N-1 ∑ k=0 L k (u k ) + g(U N ) s u b je c t to : U k+1 = f k (U k , u k ), U(0) = U 0 U m in ≤ U k ≤ U m a x (17)
where L k (u k ) is the fuel consumption over the time interval [k, k + 1], U k is the voltage of the supercapacitor at time k, f k is the function (11) at time k, and g(U N ) = β(U N -U 0 ) 2 is the penalization term (β is a constant to be chosen), N is the final time, W k the feasible domain for u k with respect to constraints (13), ( 14), and (15).

From B ellman principle, the cost-to-go function

V k (x k ) at the time step k, 0 ≤ k ≤ N -1, is expressed as V k (U k ) = m in u k ∈W k (L k (u k ) +V k+1 ( f k (U k , u k ))). ( 18
) At time N, the cost-to-go is V N (U N ) = g(U 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. A description of the classical DP algorithm can be found in (B ertsekas, 2001), we refer to [START_REF] Guilbaud | M éthodes num ériques pour la commande optimale[END_REF] for some theoretical results on the convergence and error estimations.

In the following, the control u is specialized as the feasible λ values in boost mode, and by the torque setpoint T • m in regenerative mode.

.2.1. Model implementation

A tool developed by IFP, based on dynamic programming, has been used to solve the optimization problem. The general outline of the optimization tool is described in Figure 4. It allows to plug any vehicle and powertrain model with the dynamic programming algorithm.

Each control u k of the electric motor is tested for each point of the discretization grid (t k ,U k ). As for a real-time controllers, operating conditions will forbid some or all the values of the control. The tested controls are then really applied only if these are feasible controls. Otherwise, the electric motor can not be used. Figure 5 represents the tested controls on a particular operating point. Eleven control values are tested, with five positive torques (corresponding to the feasible values of λ), five negative torques (arbitrarily chosen among feasible regenerative torque values) and one last value where the electric motor is not used. The model, executed at time t k , provides the instantaneous fuel consumption L k and the variation of the state, U k+1 , for each control u k . Finally, the DP algorithm selects the optimal control u * according to (17). model can be smaller than the dynamic programming one. Considering a smaller time step for the dynamic model permits to limit approximation errors, and to calculate a more realistic variation of voltage over a time step. In our model, the chosen time step is 0.1 s.

.O ptimiz ation Results

To calculate fuel consumption over a whole driven cycle, we use a additional matlab model, which gives exact wheel speed and torque corresponding to the vehicle parameters.

The optimal system trajectories are calculated for the NEDC cycle, and shown in Figures 6 and7. B ecause of the low global efficiency (≈ 0.8 0) of the electric system and of the large number of constraints, notably the only five allowed torques when boost is required, boost is never selected as an option. This particularity is also a consequence of the fact that the vehicle is able to follow the NEDC cycle without any need of a secondary power source to help the engine.

Figure 7 shows the transient behavior of supercapacitor voltage. Each time a regeneration is applied, the supercapacitor voltage exhibits a rise followed by a constant-rate increase. A voltage drop then ends the regeneration phase. Q uantitative results of DP vary according to the level of complexity used to represent the auxiliary power needs and torque demand from the driver. Improvements in fuel consumption obtained with hybridization range around 1% of the fuel consumption for a baseline case. In the latter, the engine only drives the vehicle while the motor uses regenerative braking to drive the auxiliaries.

Thus, accordingly to the DP results, energy management seems not much helpful in reducing the fuel consumption. However, these results must be put into perspective. First, the stop-and-start has not be taken into account in the offline optimization. This functionality should lead to interesting fuel consumption reduction. Moreover, the most stop-and-start is used, the most the energy management will be helpful to recharge and maintain a correct supercapacitor voltage, and to optimize fuel consumption. Second, some improvements on the electric models must be carried on. B asically, because of internal temperature constraints inside the electric motor, the real boost torque is smaller than the boost torque given by the model. As a consequence, during simulations, the boost mode could be applied only on few operating points, when done on NEDC cycle. Finally, the initial aim of adding an electric motor to this vehicle is to maintain the same functionalities of a conventional Smart, by using a downsized engine coupled with a starter-alternator. In this perspective, the energy management become an interesting add-on to decrease much more the fuel consumption.

REAL-TIME FEEDB ACK CONTROL

ECMS

In this section, suboptimal control law is derived from the optimization results. The chosen real-time control law is the Equivalent Consumption Minimization Strategy (ECMS): it is based on instantaneous equivalent fuel consumption. The battery is considered as a second energy source, and ECMS aims to choose the cheapest energy source, between fuel tank and battery storage. To do so, we form the Hamiltonian function H given by:

H = ṁ∆t + p Ėelec ∆t, (19)
where ṁ∆t represents the variation of fuel consumption over a time step, and Ėelec ∆t corresponds to the variation of electric energy of the supercapacitor over a time step. In this strategy, the optimal control is the admissible one which minimizes H , taking into account all constraints. The same models as the ones used in dynamic programming algorithm are here used in ECMS.

Meaning of p (L ag rang e multiplier)

The coefficient p corresponds to the price of electric energy. The larger is p, the more expensive is the electrical energy, and the more interesting to be recovered. The smaller is p, the cheaper is the electrical energy, and the more interesting to be used to power the vehicle conjointly with the engine. Some values of p corresponds to situations where it is better to just use the engine without any help of electric motor.

The variation of p depends on the time evolution of the dynamic state variable, according to the Euler-Lagrange equations. While for batteries the state variable is usually the charge, the supercapacitor state is rather in terms of voltage U. However, using simply U as the state, as done in section 3.2, would imply that U varies with U for a given electric power. This would result in an accordingly varying p.

The energy stored in the supercapacitor, E elec = 1 2 CU 2 is chosen as a state instead, with a dynamic state equation

Ėelec = CU U. (20)
With this choice, which also justifies the form (19) of the Hamiltonian, p is ideally a constant. However, this is no longer true in the presence of active state constraints. Thus, p must be determined along the cycle.

D etermination of p

The determination of p is an open question, and the initial value of p, p 0 , at the beginning of the cycle, is much more difficult to know. Many authors [START_REF] Liu | Control optimization for a powersplit hybrid vehicle[END_REF], [START_REF] Delprat | Evaluation de strat égies de commande pour v éhicules hybrides parall èles[END_REF] control p as a linear function of state of charge or considerer dependance with other measurements. Even if this 0 2 00 4 00 6 00 8 00 1 000 1 2 00 1 6

1 7 1 8 1 9 2 0 2 1 2 2 2 3 2 4 tim e v o lta g e (V )
Fig. 9. Sub-optimal voltage with ECMS choice is quite well-founded to be used in a real-time strategy, it is not optimal. However, it permits to enhance the regeneration phases when state of charge is low.

An adequate initial value can be obtained by calculating the equivalent Lagrange multiplier when using the dynamic programming results. From (B ryson and Ho, 1975) equivalence between Pontryagin's maximum principle and dynamic programming can be showed if

p(t) = 1 U ∂V ∂U (t,U), (21)
where V (t,U) corresponds to the cost-to-go function in continuous time, already defined in (18 ). The term 1 U appears here because of the difference of variation taken into account: ∆U in dynamic programming, but U∆U in ECMS.

In the following, we consider that p remains constant as long as the supercapacitor voltage trajectory stays between allowed limits, 16 V et 24 V. Previous work of the authors [START_REF] Rousseau | Hybrid electric vehicles : From optimization toward real-time control strategies[END_REF] and [START_REF] Rousseau | Scop : a sequential constraint-free optimal control problem algorithm[END_REF] has shown the behavior of p when the state touches a boundary, in simpler optimization problems. However, the optimal trajectory of p can be very complex, as it depends on the dynamics ( 11)

and also on the future active state constraints. Therefore, p is controlled by an integrator as long as the voltage U goes beyond given voltage limits.

S imulation results

The control law described above has been implemented in a complete vehicle model. To do so, we use the co-simulation between an AMESim model, and the Simulink controller. The Simulink controller includes the control of all the engine actuators, and corresponds to the one that is used on vehicle. The model, which can be seen on Figure 8 , reproduces all the main dynamics of the vehicle and of the powertrain. This model, including a description of the vehicle, has already been described in [START_REF] Tona | Amt control for a mild-hybrid urban vehicle with a downsized turbocharged cng engine[END_REF].

Figure 9 clearly shows that the electric motor is only used to recharge the supercapacitor, and that the boost is never used.

The main constraints on positive torques (only five possible controls) does not allow the electric motor to be used as a boost. B asically, the torque provided by the electric motor, as it has been considered in the model, would decrease the engine torque of about 25 Nm to 30 Nm at low engine speed, which is a high value compared to the potential of the vehicle. This consideration, coupled with global efficiencies of electric motor of about 8 0% are enough to explain the results obtained.

B etween time t = 1000 s and t = 104 0 s, the alternator mode, corresponding to a simple voltage sustaining at U = 16 .2 V, is active. The fuel consumption results are indicated in Table 1. The thermic mode result corresponds to the fuel consumption obtained with the conventional powertrain, without regenerative braking. In this configuration, the engine is used to yield the requested power for the auxiliaries, via the starter-alternator.

The second result corresponds to fuel consumption obtained while considering the previous, and the recovered energy from braking situations (the baseline in (3.3)). In normal driving conditions, the voltage decreases until the alternator gives the requested power as to maintain the voltage. The third results gives fuel consumption obtained with energy management (E.M.). In this strategy, the controller chooses the best moment to regenerate the supercapacitor thanks to the engine, instead of using the alternator mode all the time. Although the fuel consumption gain is not significant, this strategy allows to store some electric energy that can be reused later for the boost.

Veh icle results

The implemented control law has been successfully tested on the vehicle, and the correct behavior of supercapacitor voltage has been observed. Figure 10 shows the results obtained by testing the vehicle on a short portion of road. The decrease of the voltage due to the auxiliaries, which is quite high in this test because of the prototype system to be powered, is clearly observed. The variation of p remains constant while the voltage stays within 16.5 V and 23.5 V, which are its target bounds. When the voltage crosses these limits, the energy management controls p thanks to an integrator, as to decrease p and to change the price of electric energy. It should be noticed that, due again to the small energy stored in the supercapacitor, the system is very sensitive even to very small variations of p, as those shown in the figure. The electric motor torque setpoints are also represented in Figure 10 with some regenerative braking (here, the torques are limited to -4 Nm) during vehicle decelerations, which corresponds to null engine torque (fuel cut-off). The last braking, between t = 5 0 s and t = 5 1 s, does not allow to regenerate, as the supercapacitor voltage is still above its limit of 23.5 V. The energy management also controls both engine and electric The last plot of Figure 10 shows that the engine speed is not affected by the starts and the stops of the electric motor.

CONCLUSIONS

In this paper, has been presented a model-based of a mild hybrid vehicle. The model is based on physical equations of a DC electric motor, a supercapacitor, and an inverter. The model parameters have been successfully fitted to experimental data in order to simulate accurately the behavior of the electric system.

This dynamic model has been used in a dynamic programming algorithm to minimize the fuel consumption, taking into account the main physical constraints of the system (especially the discrete control settings).

The same models have been used to build an energy management controller, based on ECMS. The key parameter, the Lagrange multiplier p has been pre-established from dynamic programming results.

The first results, experimental and simulation, show that the room left for energy management in reducing the fuel consumption is rather limited. That is manly due to two reasons: (i) the supercapacitor, while sufficient for improving acceleration performance through boosting or allowing some regenerative braking due to its rather good power has a limited energy to act properly as a buffer; (ii) the stop-and-start capability has not been taken into account in the optimization, in order to focus only on the energy flows during traction phases.

The stop and start functionality, which allows to turn off the engine at idles, will also have a great impact on fuel consumption, and the availability of the electric power necessary to turn on the engine can also be optimized thanks to energy management.

Future work will therefore include some improvements on the models and some additional tests on the vehicle to explore the real potential of energy management.

  Fig. 1. Fuel consumption map of natural gas engine of VE-HGAN vehicle

TimeFig. 3 .

 3 Fig. 3. Experimental data compared to model: evolution of the voltage U s (t) where T is the final time of the drive cycle, L(u(t),t) is the fuel consumption rate, computed from the data displayed in Figure 1, g(U(T ), T ) is a term that penalizes any deviation of the final voltage from the initial voltage.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Diagram of the optimization tool A standard time step used in our examples, for dynamic algorithm, is 1 s. However, the time step of the Simulink

Fig. 6 .

 6 Fig. 6. Optimal trajectory of the supercapacitor voltage Figure 6 clearly shows that the final voltage constraint is approximately verified. The electric energy delivered by the supercapacitor to feed the auxiliaries (constant-rate voltage decreases) is compensated by the energy recuperated either from regenerative braking or via the engine (generator mode). B ecause of the low global efficiency (≈ 0.8 0) of the electric system and of the large number of constraints, notably the only

Fig. 7 .

 7 Fig. 7. Optimal trajectory of the supercapacitor voltage (zoom)

Fig. 8 .

 8 Fig. 8 . AMESim model used in co-simulation motor to recover some electric energy, see time t = 3 0 s and t = 4 2 s, by increasing the engine torque, with respect to the driver torque, to increase engine efficiency. In this plot, the driver torque corresponds to the sum of the engine torque and the electric motor torque, T e = T rq -ρT m .

Fig. 10 .

 10 Fig. 10. Evolution of voltage, p, torques, and speed during a road test

Table 1 .

 1 Sub-optimal fuel consumption values

	Th. mode	Th. mode & b ra k e rec u p . (b a s elin e)	E .M .
	4 5 8 g	4 5 4 .1 g	4 5 2 .7 g

ACK NOWLEDGMENTS

We would like to thank Paolo Tona and Rachid Amari for their help and for the time spent during the tests on vehicle. We also acknowledge Dominique Soleri to let us quite free to test our strategy on the VEHGAN vehicle. We would like to thank Pierre Rouchon for helpful discussions and advice at various stages of the elaboration of this work.


Frédéric Delbos

Kaoutar Sghiouer

Uncertainty Analysis in Prestack Stratigraphic Inversion: a first insight Authors

Stratigraphic inversion of prestack seismic data allows the determination of subsurface elastic parameters (density, P and S-impedances). Based on a Bayesian approach, the problem is formulated as a non-linear least-squares local optimization problem. The objective function to be minimized is composed of two terms, the first one measures the mismatch between the synthetic seismic data (computed via a forward operator) and the observed seismic data, the second one models geological a priori information on the subsurface model. It is crucial to estimate the a posteriori uncertainties because the solution model of the inversion is only one solution among the range of admissible models that fit the data and the a priori information . The goal of this paper is to propose an optimized deterministic method to estimate a posteriori uncertainties in stratigraphic inversion. The proposed method is based on the hypothesis that the covariance matrices describing the uncertainties on the data and on the model are laterally uncorrelated (no cross correlation among parameters of different traces). Moreover, the covariance matrix on the data is also supposed laterally stationary. Application on 2D synthetic PP data illustrates the performances of the method. Extension and limitations of the method are discussed.

Introduction

In a Bayesian approach, which is a natural formulation for many geophysical problems, solving an inverse problem consist not only in estimating the best probable set of subsurface parameters but to estimate the posterior probability density function (PDF) on the model space. This is particularly true in seismic facies analyses of 4D data where the noise level is in the magnitude of interpreted parameters : from on vintage to the other very small differences in subsurface elastic parameters are correlated a posteriori to production effects (see [START_REF] Nivlet | Uncertainties in Seismic Facies Analysis for Reservoir Characterisation or Monitoring : Causes and Consequences[END_REF]. Presently, numerical methods for estimation of posterior PDF can be gathered into at least two classes: the class of stochastic approaches which is mainly formed of Markov Chain Monte Carlo method, and the class of deterministic approaches which includes linearized a posteriori analysis. Often actual methods combine elements of the two classes, but their main features make them belonging to one of them. In stochastic approaches, no assumption on the posterior PDF function is needed. The Monte Carlo method consists in sampling the prior PDF function in the model space and then uses a Metropolis type algorithm with a random walk to map the posterior PDF function (see Mosegaard 1998 for example). Considering the large size of the model and data space (up to 10 8 unknowns and data samples in stratigraphic inversion) and the computational cost of the forward modelling, this class of method is often rejected. Deterministic approaches are based on a gaussian assumption of the posterior PDF and also of all the 'input' probabilities densities. This assumption provides an explicit analytical form of the posterior distribution which can be computed at a lower computational cost (see [START_REF] Tarantola | Inverse Problem Theory[END_REF]. The computation of the posterior distribution has been the subject of many algorithm studies, we mention the following techniques: frequency methods using Fast Fourier Transform (FFT), and direct methods using factorisation techniques (LU factorisation for example).

In frequency methods, one solves a sequence of small decoupled linear equations obtained by transforming a linearized inversion problem into the Fourier domain (see for instance [START_REF] Buland | Rapid spatially coupled AVO inversion in the Fourier domain[END_REF]. Thanks to FFT, this method has the advantage to be very cheap in computational cost. Its drawback comes from the restrictive assumptions to be postulated: the covariance function of the prior model has to be stationary and homogeneous (which not allows applications to complex media with high lateral variations or with more than one geological unit), the ratio between the P-wave velocity and the S-wave velocity is assumed to be constant (which is not the case in general) and finally this method works on a linearized AVO inversion problem (with the loss of the complexity of wave propagation). Direct methods do not need to work on a linearized version of the inversion problem. Then in a non linear context, uncertainty analysis consists in assuming a gaussian approximation in the vicinity of the optimal subsurface model which has been previously estimated via a non linear optimization method (see [START_REF] Gouveia | Bayesian seismic waveform inversion: Parameter estimation and uncertainty analysis[END_REF]. Then this approach is less restrictive than the frequency one. Nevertheless, the estimation of the posterior distribution is a more difficult task since it requires the inversion of a very large matrix (its size is n×n where n is the number of unknowns).

In the following, we first briefly describe an optimized method of posterior uncertainty analysis, based on the deterministic approach coupled with a direct. Then the method is applied on a 2D synthetic seismic data set corresponding to 2D realistic reservoir zones where optimal model and its error bars are plotted.

Uncertainty Analysis for Prestack Stratigraphic Inversion

For uncertainty analysis we use the prestack stratigraphic inversion algorithm described by [START_REF] Tonellot | Joint stratigraphic inversion of angle-limited stacks[END_REF] which minimizes the objective function

(m) J (m) J m m d w * (m) R J(m) G S 2 C prior θ 2 C obs θ θ θ 1 m 1 d + = - + - = - - ∑ P S ρ P S ρ
composed of two terms. J S (m) is the seismic term where R Ө (m) is the Aki-Richards reflectivity series corresponding to the current model m and to the angle Ө, w Ө is the wavelet, and d obs Ө is the observed seismic trace. J G (m) is the geologic term where m prior is the elastic a priori model. C d and C m are the prior covariance matrices describing uncertainties on the data and on the model, they are assumed here laterally uncorrelated, C d is also assumed laterally stationary. The maximum a posteriori (MAP) model is the model m post which minimizes J(m). Then, by a linearization of R Ө (m) operator, a gaussian approximation of the posterior PDF centered on m post is obtained. This distribution is then quantified by the posterior covariance matrix

1 1 m S 1 1 m 1 d T post ] C H [ C G W C W G [ C - - - - Τ - + = + = ∑ θ θ θ
, where G is the jacobian matrix of R Ө (m) in m post and W Ө is the convolution matrix of the wavelet at angle Ө. It is straightforward to show that the first term H S is a Gauss-Newton approximation of the hessian of J S in m post . The second term, which is the inverse of the covariance prior matrix on the model is exactly the hessian of J G . As the prior covariance matrices are diagonal, C post is a n T ×n T block diagonal matrix (the parameters are laterally uncorrelated) where n T is the number of traces . Then computing C post requires two steps:

(1) the computation of the matrix H S , and

(2) the inversion of the n T diagonal blocks of C post . The first step can be very high in computational cost if not made properly. We have optimized it using two improvements. The first one is to store in a matrix Z Ө the result of the matrix products W Ө C -1 d W T Ө . The second one is to use the following explicit analytical formula

( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 1 , , 2 , 2 , 2 , , 1 , 1 , 1 1 2 1 1 1 , , 2 , 2 , 1 2 , , 1 , 1 , 1 1 ) ( ) ( ) ( ) ( - - - - - - - - ⊗ ∂ ∂ ∂ ∂ + ⊗ ∂ ∂ ∂ ∂ t t
1 1 2 1 1 , , ( ) ( ) ( ) ( ) ( ) 2 1 1 , , 2 , 2 , 2 , , 1 , 1 ,
m m R m m R m m R m m R - + - - ⊗ ∂ ∂ ∂ ∂ + ⊗ ∂ θ θ θ ω ω ω ω
which gives the value of the element ((t1,p1,x,y),(t1,p1,x,y)) of the matrix H S (Ө). The 4 uplet (t,p,x,y) is a parameter coordinates in the subsurface grid (t is the time index, p the parameter type index, x and y the crossline and inline indices). have been previously stored, the computational cost of H S is then optimized (one no need to recompute these values from one trace to the other). The second step (the inversion of the matrix C post ) is more conventional: it can be processed very fast using direct method such as a LU factorisation (see [START_REF] Anderson | Lapack user Guide, Third Edition[END_REF]. In fact, the size of each diagonal block is small (around 1000×1000) and this step is naturally parallelisable (each matrix inversion being send to one processor).

A typical diagonal block of the matrix H S is shown in Figure 1. This 3×3 block matrix shows the autocorrelation of the P-wave impedance, S-wave impedance, and density on its main diagonal and the cross-correlation among different parameters in the offdiagonal blocks.

2D synthetic study

The exact model (see raw 1 of Fig. 3 for the elastic parameters) is a 2D crossline of a 3D realistic model representative of a turbiditic complex channel in a deep offshore environment (see [START_REF] Bourgeois | How to build a geological and geophysical realistic synthetic seismic data set for benchmarking reservoir studies?[END_REF]. The inversion window is 3km long and 960m thick (832 ms) from 1.5km to 2.46km in depth (time range is 600-1432ms). This area is discretized on regular grid of 151 traces (one each 20m) and of 417 time samples (one each 2ms). Using the exact model, the P-P Aki-Richards reflectivity formula and a 1D convolution operator, a seismic data set is generated. Four angle-limited stacks are generated: 05°-15° (see Fig. 2), 15°-25°, 25°- The prior covariance matrices C d and C m are choosen diagonal and respectively filled with a 1% noise-to-signal ratio and a 5% standard deviation with respect to the average values of exact elastic parameters: σ P prior =465, σ S prior =270 and σ ρ prior =0.12.

Then, the MAP model m post (see raw 2 of Fig. 3) and its associated residuals (see Fig. 2) are computed through the inversion process. The estimated elastic parameters are well recovered : m post is very close from m exact and the residuals are almost null (noiseless seismic data). Using the method described previously an uncertainty analysis is processed : absolute (σ) and relative (σ rel : σ normalysed by σ prior ) standard deviations for m post are respectively displayed in raw 3 and 4 of Figure 2. That provides insight into the different degree of data resolution for the 3 elastic parameters. The prior standard deviations are reduced to the ranges of : 40-100 (8%-20%) for the P-impedance parameters, 44-80 (16%-30%) for the S-impedance and 0.066-0.08 (54%-67%) for the density. It confirms the already known fact that for P-P data the P-impedance is the best resolved parameter, followed by the S-impedance and then the density : left and right parts of Figure 4 respectively shows that the P-impedance is 1.2 to 1.3 times better resolved than S-impedance and 3 to 5.5 times better resolved than the density. Note that, this analysis is neglecting the off-diagonal terms of the posterior covariance matrix.

A more detailed analysis would consist in generating pseudo-random realizations of posterior subsurface models. This analysis is straightforward once C post has been evaluated.

Conclusions

The proposed method, optimized for laterally uncorrelated prior matrices and a laterally stationnary prior matrix on the data, gives a first insight towards a complete uncertainty analysis. Generalization to multicomponent (see [START_REF] Devault | Prestack 9-C joint inversion for stratigraphic prediction in the Williston Basin[END_REF] or 4D prestack joint inversion is in principle straightforward. The stated assumptions may be restrictive. Nethertheless, for some realistic case, where for example the horizontal correlation length is small with respect to the horizontal dimensions of the model, choosing diagonal prior matrices may be appropriate. Moreover, a patent is being written on the extension of this method using exponential covariance operator on the prior model.

  t

Ө

  and ⊗ is the cross-correlation product. As the cross-correlation values

Figure 1 :

 1 Figure 1: one diagonal block of the posterior covariance matrix

Figure 2 :

 2 Figure 2: (near section) synthetic seismic (left); residuals (right)

  (

	t			x	y	t			x	y	t	t			x	y	t	x	y	t
											t	t					t	t
	t	p	x	y		t	p	x	y			t	p	x	y	t	p	x	y

Acknowledgments

We would like to thank Patrick Lailly and Thierry Tonellot for helpful discussions. We acknowledge Karine Labat, Karine Broto and Vincent Clochard for their precious help for the application.


Hoël Langouët
email: hoel.langouet@ifp.frludovic.metivier@ifp.frdelphine.sinoquet@ifp.frq-huy.tran@ifp.fr

Ludovic Métivier

Delphine Sinoquet

Quang-Huy Tran

Optimization for engine calibration

Keywords: Engine calibration, LoLiMoT, Multi-objective optimization, Evolutionary algorithm

Nowadays, automotive manufacturers are submitted to strong constraints in engine calibration such as: low fuel consumption, emission-control legislation and driver requests for driving comfort and performances. These constraints lead to an increasing complexity of the engines and thus an increasing number of parameters to be tuned, making the empirical engine calibration by a scan of parameter values impossible at engine test-bench. New methodologies in automated engine calibration based on statistics and optimization have emerged in order to limit the number of experimental tests to be run.

The optimization problem of engine calibration consists in the determination of engine tuning parameters that minimize the cumulated fuel consumption and pollutant emissions on a driving cycle generally associated with legislation norms. This cycle is decomposed in a set of stationary operating points of the engine characterized by its speed and its torque (the transient behaviors of the engine are not taken into account in the stabilized calibration). Then, the optimal tuning parameters of the engine should be defined for each operating points, the functions defining these parameters on the whole engine operating domain are called the engine maps. These two-dimensional optimal engine maps are then integrated in the engine control unit in the vehicle.

We illustrate the difficulties associated with this application and propose adapted optimization methodologies: LoLiMoT models for engine map parameterization in order to handle intrinsic constraints on the map regularity, multi-objective optimization method based on CMA-ES approach. Finally, application on real dataset obtained at IFP automated test-bench for a diesel engine are presented.

2.

Introduction

Engine calibration consists in fulfilling the engine tuning maps that are used in engine controls of the vehicle, i.e. in defining the optimal tuning of parameters used by engine control strategies. Due to the highly increased number of these parameters (especially for diesel engines but spark ignition engines are following the same trend) and the reduction of the development schedule available for the calibration process, manual tuning of engine parameters is now replaced by mathematically assisted calibration process. Such a process is based on the design of experiments with associated modeling methods, in order to reduce the number of tests used to build engine response models depending on engine control parameters, and optimization techniques to determine the optimal settings within the model definition domain. In order to perform the tests in a more productive way, these mathematical techniques are generally associated with test automation, requiring well controlled measurement methods and reliable test equipments.

This paper describes the optimization methods developed for this application and illustrates their effectiveness on a real case of a common rail diesel Engine. The first section introduces the classical steps of the calibration process and discusses the associated difficulties. In the second section, we propose the Multi-Objective Covariance-Adaptation Evolutionary Strategy method for solving the optimization problem associated with a given engine operating point defined by the engine speed and the engine load. In the third part, an integrated approach is proposed in order to directly optimize the engine maps on the whole driving cycle (associated with legislation norms) instead of the individual optimization of each engine operating point.

Engine calibration 4.1. Sketch of the engine calibration process

The emission calibration workflow is classically divided into four steps:

1. a preliminary phase consisting in choosing a sample of operating points (referred to as OP in the Figure 2: Projection of the domain of parameter variations in the (boost pressure, air mass) plane. The limits of the domain are generally modeled by linear constraints on the engine control parameters. The ellipse indicates the parameter dispersions (here projected on the 2D space) modeled by a multidimensional Gaussian probability law. The ellipsoïd should remain in the engine physical domain in order to avoid that any control parameter takes some value outside this domain after the implementation in the control unit. following) to be studied and emissions targets associated with each OP, these targets being called allocations hereafter.

2. the optimization of engine responses on each OP according to these targets, 3. the fulfilling of the maps with a smoothing step between these optimal settings. The preliminary step deduces from a simple model of the vehicle the trajectory of the driving cycle in the engine speed-load operating domain (see Figure 1): the transient effects on the accelerations of the cycle are neglected and the cycle is thus considered as a sum of stabilized points. A limited number of specific Operating Points are chosen to represent the cycle in the engine working range. Figure 1 gives an example of NEDC simulation and a selection of 17 OP. The cumulated levels of pollutant emissions along the cycle are computed as weighted sums of the pollutant emissions for each chosen OP. i) defining the domain of variations of the engine control parameters: this is an essential step of the process as it defines the validity domain of the models. The complexity of the models to be used for engine response depends on the size of this domain: for tiny domains low order polynomials (second order) are usually sufficient to accurately model engine responses. However, choosing too small domains leads to difficulty in coherently fulfilling sub-optimal engine maps.

ii) building the test matrix: various types of experimental designs can be used to build a test matrix: D-Optimal, space filling... The choice of the type of design as well as the number of tests to be done are directly correlated with the assumed complexity of the model and thus with the size of the considered domain. D-Optimal test designs are often used with hyper-cubic tiny domains.

iii) running this test matrix on the test bench: as the tests are predefined, the experiment can be performed in an automated way, which drastically improves the productivity of the global process. In this case, special attention must be paid to the validation of the experimental data.

iv) modeling the engine responses: the type of mathematical models depends on the complexity of the engine responses and on the size of the domain of parameter variations. In practice, the engineer often prefers to limit the parameter variations in order to avoid critical engine tunings that may lead to difficulty to realize the experiment at the test bench. However, as already mentioned, reducing the size of the parameter variations may lead to difficulties to build coherent optimal engine maps from the obtained optimal settings in those tiny domains. Often, low order polynomial functions are used coupled with D-optimal design of experiment. For some response as HC or CO emissions, this type of models may be too limited, then non-parametric response surfaces may be chosen as RBF or kriging. We will not detail this step of the calibration process in this paper, the reader may find more details on this topic in [START_REF] Castagné | Engine calibration: towards an integrated approach[END_REF].

v) optimizing the engine control parameters to meet the allocations. The problem may be formulated as a classical mathematical problem of optimization under constraints or as a multi-objective optimization (searching for compromises between antagonist objectives) as it is detailed in the following section. For this classical approach, the optimization is performed one OP after the other, considering the allocations of each OP as constraints. We refer in the following to this optimization problem as the OP optimization problem.

When the optimal settings are found, the last step consists in integrating them in the reference engine maps (if available) or building maps from these settings. For the drivability of the target vehicle and because sharp evolution of air loop parameters are not easily feasible during transient, it is necessary to provide smooth engine maps. Thus, the settings are often moved away from their optimal values in order to build smooth engine maps, especially for air loop parameters. This smoothing process can be performed with respect to local constraints (such as maximum gradients), as well as to keep some predefined shapes. The difficulty of this step is thus to stay as close as possible to the local optima while preserving a smooth map shape, in order to keep all the benefits of the optimization work and satisfy the targets. This delicate step is usually time consuming and is likely to deteriorate the work performed during the optimization phase by providing sub-optimal engine tunings. This is the motivation for our proposal of an integrated method of map optimization described in section 6.

The OP optimization problem formulation

The OP optimization problem consists in minimizing some engine responses under constraints on others engine responses whereas the engine control parameters are kept inside the domain usually modeled by linear constraints coupling the parameters (see Figure 2). It leads to the following multi-objective constrained optimization problem:

                         min p∈R Np (f 1 (p), f 2 (p), . . . , f m (p)) subject to l ≤ Ap ≤ u f 1 (p) ≤ s 1 f 2 (p) ≤ s 2 . . . f d (p) ≤ s d (1)
In turbo-diesel engine, a classical bi-objective optimization problem is the minimization of the particulate emissions and the NOx emissions under constraints on fuel consumption (CO 2 emissions), on CO and HC emissions and engine noise level. Generally, this multi-objective optimization problem ( 1) is replaced by a constrained optimization problem where only one of the objective is effectively minimized, the other being constrained to remain below some thresholds:

             min p∈R Np (f i (p)) subject to l ≤ Ap ≤ u f j (p) ≤ s j j ∈ [1, d] (2)
This constrained single-objective optimization problem ( 2) is then solved by a classical SQP method run several times with different initial points to search for the global optimum.

Solving the OP multi-objective optimization problem

In this section, we present a method to solve the constrained multi-objective problem [START_REF] Nelles | Non linear system identification[END_REF]. It is based on the algorithm proposed by [START_REF] Hansen | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF] namely the Multi-objective Covariance Matrix Adaptation Evolutionary Strategy method (MO-CMA-ES). The MO-CMA-ES method belongs to the class of evolutionary strategy methods which allow to solve global optimization problem. These methods are popular for their ability to find a global minimum but their slow convergence rate is often criticized. In [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF], Hansen has proposed the CMA-ES method for optimization of real-space functions, based on an original mutation operator that increases the convergence rate thanks to the adaptation of the mutation distribution at each generation.

CMA-ES method for single-objective optimization

During the CMA-ES evolutionary algorithm, a new population is created from the current population by generating realizations of a Gaussian probability law

x (g+1) k ∼ N m (g) , σ (g)2 k • C (g) k (3)
where m (g) is the barycenter of the best µ individuals with respect to the objective function to be minimized,

m (g) = µ i=1 w i x (g) i:λ , with µ i=1 w i = 1 and w µ ≤ w µ-1 ≤ . . . ≤ w 1 ≤ 1, (4)
x g i:λ being the i th best individual of the current population x (g) with respect to the function. The next population is thus directed towards the best points. An ideal choice for the covariance matrix would be the inverse of the Hessian matrix of the objective function, the idea proposed by Hansen is then to approximate this matrix as for a Quasi-Newton method but without any computation of the derivatives. The update formula for the covariance matrix is then

C (g+1) = (1 -c cov )C (g) + c cov µ i=1 w i x (g+1) i:λ -m (g) σ (g) x (g+1) i:λ -m (g) σ (g) T (5)
which amounts to adding to the current covariance matrix a µ-rank term that takes into account the repartition of the best µ individuals of the new population. An other update formula is proposed by Hansen which does not take into account only the last iteration as in ( 5) but all the previous iterations from the beginning. In the same way, an update formula of the standard deviation is built: the main idea is here to avoid that two consecutive iterations lead to a displacement of the barycenter of the population m (g) in the same direction, replacing these two iterations by one with a larger standard deviation would have been a better choice. Likewise, a small standard deviation is preferred to a large one which leads to two iterations with opposite displacements of the barycenter of the population. The reader is referred to [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF] for more details. Figure 3 illustrates on a 2D single-objective optimization example the adaptation of the covariance matrix and of the standard deviation leading to an efficient evolution of the population towards the global minimum. The non-dominated solutions of A are given by: ndom(A) = {a ∈ A | a ∈ A : a ≺ a}. We consider that these non-dominated solutions have a non-dominance level of 1. Recursively, we define the points of non-dominance level equal to 2 as the non-dominated points of the population A without the non-dominated solutions. We thus obtain the recursive definition:

ndom i (A) = ndom(dom i-1 (A)) with dom i (A) = dom i-1 (A)\ndom i (A) and dom 0 (A) = A. (7)
2. The second criterion is a diversity criterion of the solutions to obtain an uniform repartition of the solutions along the Pareto front. A measure of the dominance hyper-volume at a non-dominated solution is defined by considering the surface of the rectangles defined in Figure 4 (one rectangle for each non-dominated point). Combining these two criteria allows to sort the individuals of the current population with respect to the multiple objectives and then to apply the CMA-ES update formula of the covariance matrix adaptation for the mutation. For more details on the MO-CMA-ES method, the reader is referred to [START_REF] Hansen | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF].

Application of MO-CMA-ES method to a real case of engine calibration

In this section we present results obtained from a real dataset associated with a turbo-diesel engine. We consider the optimization problem (1) at a given engine operating point (1750 rd/mn -1 bar). Models of the engine responses have been built from experimental data, the tests being designed thanks to a classical D-optimal criterion: the models are cubic polynomials.

The optimization is performed on six engine control parameters: the main injection timing, the pilot fuel injection quantity, the pilot injection timing, the fuel injection pressure, the mass air flow and the boost pressure.

The optimization objectives are the particulate emissions and the NOx emissions, CO2, HC and CO emissions and the engine noise being constrained to remain smaller than a given threshold. Linear constraints defining physical limits of the domain of variations of the parameters are introduced. These linear constraints may take into account the parameter dispersions, these new constraints are denoted by robustness constraints hereafter.

Thus, the multi-objective optimization algorithm has been modified to deal also with linear and non linear inequality constraints: they have been introduced via a l 1 penalization term added to the objectives (if the i th constraint is not satisfied, the misfit with the bound is added to the objective functions, the penalization weight varying during the iterations).

The results are presented in Figure 5 : two optimizations are compared, one without taking into account the parameter dispersions and one with the robustness constraint which consists in taking into account the possible dispersion of the engine control parameters when applied on the vehicle. This dispersion is modeled by an uncorrelated Gaussian probability law, thus by one standard deviation by parameter. As shown in Figure 2, a vicinity of the engine setting is defined and this vicinity should remain inside the physical limits. The impact of this constraint is noticeable in Figure 5 by an important shift of the Pareto front and a modification of its shape. Engine experts were very interested in the different compromises between the objectives that have been obtained by the multi-objective optimization method. The points indicated in red, green and blue are associated with successive runs of a local single-objective optimization with different initial points and with a weighted sum of the two objectives (different weight values have been tested). We notice that these optimizations lead to a very few points on the Pareto front: the other points are some local minima. Moreover, as expected, the points obtained on the Pareto front are located on the convex parts only.

The MO-CMA-ES method has performed on this example 3 times more evaluations of the objective functions than the local optimization with different initial points and different weights. Even if this performance is already very encouraging, some work on the stopping criteria may reduce the number of evaluations of MO-CMA-ES.

Solving the map optimization problem

In the classical calibration process described in section 4, the engine map building in the smoothing step is delicate and time consuming. In this section, we propose an alternative method that consists in directly optimizing the cumulated engine responses over the cycle via deformations of the engine maps themselves instead of optimizing individually the selected OP and building afterward the engine maps by the smoothing step. The map optimization problem is formulated as

                   min m p ∈R Np T 0 F i (r(t), c(t), m p1 (r(t), c(t)), m p2 (r(t), c(t)), . . . , , m pN p (r(t), c(t))) dt subject to l(r, c) ≤ Am p (r, c) ≤ u(r, c) T 0 F j (r(t), c(t), m p1 (r(t), c(t)), m p2 (r(t), c(t)), . . . , , m pN p (r(t), c(t))) dt ≤ S j j ∈ [1, d] (8)
where (r(t), c(t)) indicate the trajectory of the cycle within the engine speed-load domain, F i is the model of the engine response i depending on the engine control parameters but also on the speed and load, m pi are the 2D engine maps of the control parameters. The objectives to be minimized (or constrained) are the engine responses cumulated on the considered driving cycle: the cumulated responses are here (see [START_REF] Castagné | like to thank Yohan Entails[END_REF]) the weighted sums of the local models defined at chosen representative OP but if global models of the engine responses are available (models including engine speed and load dependency), a finer optimization may be performed with a fine sample of the integrals in the original optimization problem (8). Additional smoothing constraints such as global smoothing constraints (to preserve the regularity of the original maps) or more local constraints (for example limits on the gradients of the maps) can also be introduced. min

m p ∈R Np ×R N OP NOP l=1 ω l F (l)
i (m p1 (r l , c l ), m p2 (r l , c l ), . . . , , m pN p (r l , c l )) (9)

Modeling the engine maps

The formulations (8) and ( 9) require an adapted parameterization of the engine maps that must be flexible enough to model the very different shapes of engine map surfaces (Figure 6) and should not require too many parameters to limit the number of unknowns in the optimization process. LoLiMoT models ([1], [START_REF] Hafner | The use of stationary and dynamic emission models for an improved engine performance in legal test cycles[END_REF]) seem to be a good compromise between flexibility, accuracy and complexity: some very simple local models (linear or bilinear) are combined by a weighted sum

m(r, c) = M i=1 mi (r, c)Φ i (r, c) with mi (r, c) = ω 0i + ω ri r + ω ci c (10)
where the weights Φ i (r, c), normalized Gaussian functions, control the degree of smoothness of the global surface: This representation allows an adaptive refinement of the surface: the patching associated with the definition domains of the local models may be refined during the optimization process. A finer patching allows a finer optimization (the number of degrees of freedom being increased) but may lead to a cumbersome optimization of a large number of parameters. The parameterization, namely the patch definition for LoLiMoT description of the maps should reflect the degree of smoothness the user expects for the maps: for some parameters like boost pressure, the map should remain smooth, for others like main injection timing, the smoothing constraint is not as strong. Thus, from some reference engine maps or some a priori information, a LoLiMoT parameterization of each map of engine tunings is defined (Figure 6). The unknowns of the optimization are the LoLiMoT parameters (coefficients of local linear models ω .i ). 6.2. Application of map optimization on a real case The map optimization has been applied on the diesel engine application presented in section 5.3. The optimization problem is defined as a constrained single-objective optimization problem where the cumulated particulate emissions are minimized and the cumulated CO2, CO, HC and NOx emissions are constrained to remain under thresholds (given in table 1). Linear constraints are added to take into account the limits of the physical domain. The obtained results are listed in table 1) and the optimized engine maps are displayed in Figure 7. These engine maps have been run at the test bench and particulate emissions have been measured, those measures are displayed in Figure 8.

Φ i (r, c) = µ i (r, c) M j=1 µ j (r, c) with µ i (r, c) = exp - 1 2α (r -r i 0 ) 2 σ r i 2 + (c -c i 0 ) 2 σ c i 2 (11)

Conclusions

In this paper, the optimization problems encountered in engine calibration are addressed. The constrained multi-objective optimization has been successfully solved by the MO-CMA-ES method modified to handle inequality constraints. The performance of this method on the considered real case application is very encouraging. The resulting Pareto front gives to the engine experts a worth information on the possible compromises between antagonist engine responses as particulate emissions and NOx emissions of a diesel engine. In a second part, a new formulation of the engine map optimization is proposed based on a cycle optimization instead of an optimization at given engine operating points. This formulation allows to handle the smoothing constraints on the maps directly in the optimization process, thus to avoid the generally cumbersome step of map smoothing in classical engine calibration process. Results of a constrained single-objective map optimization have shown the efficiency of the method. A next study will be to apply a multi-objective approach as for the OP optimization, the difficulty will be the large size of the parameter space ( 500 parameters to describe all the engine maps). An other improvement in this approach will be the use of global models of engine responses (defined on the whole engine operating domain instead of the local models at a limited number of OP): this optimization is expected to yield finer results.

Figure 1 :

 1 Figure 1: Simulation of the cycle from European legislation (NEDC: New European Driving Cycle) and selection of operating points

Figure 3 : 1 2Definition 1 .Definition 2 .

 3112 Figure 3: An application of CMA-ES method on the analytical function f (x, y) = x 2 +y 2 10 -

Figure 4 :

 4 Figure 4: Contributing dominance hyper-volumes for the evaluation of the diversity criterion

Figure 5 :

 5 Figure 5: Pareto fronts (in dark blue) obtained with the MO-CMA-ES method for the OP 1750 rd/mn -1bar. The front at the right side of the Figure is associated with an optimization taking into account the parameter dispersion (robustness constraints that ensure the effective engine control parameter values to remain inside the engine physical domain). The other front results from an optimization without these constraints. The pink points are solutions obtained by single-objective optimizations (weighted sum of the two objectives) associated with different initial points and different penalization weights.

Figure 6 :

 6 Figure 6: Examples of LoLiMoT parameterization of engine map. Left: patching of the domain. A local linear model is defined on each patch. Right: the resulting surfaces modeled by LoLiMoT for two different engine control parameters from given discrete reference maps.

Figure 7 :

 7 Figure 7: Optimal engine maps obtained by the map optimization on the NEDC driving cycle. 378 LoLiMoT parameters are necessary to model the 6 maps with the expected regularity. The green points indicate the engine operating points for which models of engine responses are available.

Figure 8 :

 8 Figure 8: Comparison of measures of particulate emission on the test bench for the initial engine maps and the optimal ones obtained by the proposed map optimization method.

Table 1 :

 1 Results of map optimization on a diesel engine application.

		Initial maps Optimized maps Upper bounds
	Particulate (g/km)	0.128	0.090	0.090
	NOx (g/km)	0.219	0.240	0.241
	HC (g/km)	0.137	0.150	0.150
	CO (g/km)	0.815	0.847	0.896
	CO2 (g/km)	168.89	164.40	168.89