Gregory Rousseau
email: gregory.rousseau@ifp.fr

Delphine Sinoquet
email: delphine.sinoquet@ifp.fr

Antonio Sciarretta
email: antonio.sciarretta@ifp.fr

Yohan Milhau
email: yohan.milhau@ifp.fr

Design Optimization and Optimal Control for Hybrid Vehicles

Keywords: Hybrid vehicles, optimal control, dynamic programming, Pontryagin's principle

Growing environmental and global crude oil supplies concerns are stimulating research on new vehicle technologies. Hybrid-electric vehicles appear to be one of the most promising technologies for reducing fuel consumption and pollutant emissions. Different types of hybrid-electric powertrains exist: from the mild-hybrid vehicle, equipped with a small electric motor, to the combined hybrid like the Toyota prius.

This paper presents a parametric study focused on variations of the size of the powertrain components, and optimization of the power split between the engine and electric motor with respect to fuel consumption.

To perform this optimization on a prescribed driving cycle (for instance, the New European Driving Cycle), a dynamic programming algorithm based on a reduced model is implemented. This simplified model allows a fast optimization with a fine parameterization of the controller: it furnishes the optimal power repartition at each time step regarding fuel consumption under constraints on the battery state of charge. The obtained results may be used to determine the best components of a given powertrain, for a prescribed vehicle cycle.

The optimal split obtained thanks to dynamic programming algorithm can not be used directly on a vehicle as a real time control law, as the future can not be known in advance in normal driving conditions. To overcome this difficulty, we implement, as a real-time strategy, the Equivalent Consumption Minimization Strategy (ECMS): the battery being considered as an auxiliary reversible fuel reservoir, an instantaneous minimization of ECMS is performed. This control law is inferred from Pontryagin's Minimum Principle, where the Lagrange multiplier can be deduced from previous optimization results on given driving cycles. Offline optimization results and real-time control laws are compared for a realistic hybrid vehicle application.

2.

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 [START_REF] German | Hybrid powered vehicles[END_REF]. Indeed, it is well known that the kinetic energy gained by the classical vehicle after an acceleration can not be recuperated when braking. One of the advantage of the hybrid vehicle is its ability to recuperate a part (depending of efficiencies) of this energy, to store it in a battery, and to re-use it latter. This electric energy can therefore be used by an electric motor to power the vehicle and to save fuel.

The control of hybrid powertrains is more complex than control of classic engine. Actually, the control laws have to deal with the state of charge of the battery, which provides the level of remaining electric energy, and with the variable efficiency of each element of the powertrain. Optimization of energy management strategies on given driving cycles is often used to derive sub-optimal control laws to be implemented on the vehicle (see among others [START_REF] Sciarretta | A Real-Time Optimal Control Strategy for Parallel Hybrid Vehicles with on-board Estimation of the Control Parameters[END_REF], [START_REF] Scordia | Systematic elaboration of online energy management laws for hybrid vehicles[END_REF], [START_REF] Wu | Optimization of Power Management Strategies for a Hydraulic Hybrid Medium Truck[END_REF], [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF]).

In this paper, we focus on a parallel hybrid powertrain. The vehicle used for this parametric study is a Citroën C1, with a gasoline engine capacity of 1000cc. From this vehicle, we consider an additional source of power given by an electric motor located near the wheels. This motor takes its energy from an additional battery. A model of this powertrain has been implemented in AMESim * environment for validation purposes, and also in Matlab, for optimization. The latter model is a simplified static model, that allows fast computations. Two different optimization algorithms are presented and applied on this model, depending on the needs:

• a classical Dynamic Programming algorithm ( [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 | Systematic elaboration of online energy management laws for hybrid vehicles[END_REF], [START_REF] Wu | Optimization of Power Management Strategies for a Hydraulic Hybrid Medium Truck[END_REF]), which provides the optimal control strategy on a prescribed cycle, and the potential of a powertrain.

• a sub-optimal control law, based on the Principle of Pontryagin, used as a real-time control law.

System Modelling and Optimal Control Problem

Optimal Control Problem

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

                 min u J(u) := T 0 L(x(t), u(t), t)dt + g(x(T ), T ) subject to: ẋ(t) = f (x(t), u(t), t), x(0) = x 0 x min ≤ x(t) ≤ x max u min (t) ≤ u(t) ≤ u max (t) (1)
where 0 and T are respectively the initial and the final times of the prescribed driving cycle, u(t) is the control, x(t) is the state, L(x(t), u(t), t) is the instantaneous fuel consumption, f (x(t), u(t), t) is the function that controls the variations of the dynamic system, and g(x(T ), T ) is a penalization term handling a constraint on the final state of charge of the battery.

In this problem, the state has two components:

x(t) =
x 1 (t) : the state of charge (SOC) of the battery x 2 (t) : the state of the engine (on/off).

(

) 2
While x 1 (t) is a continuous state governed by ẋ1 (t) = f 1 (u(t), t), the second state x 2 (t) takes only the values 0 or 1 corresponding to engine stop / engine start. The introduction of x 2 (t) permits to take into account an extra fuel consumption when the engine is being started.

The control vector has three components, to be optimized:

u =    u 1 :
the power split ratio between the engine and the electric motor u 2 : the gear u 3 : starting up / stopping the engine

where u 1 is the classical control for HEV, u 2 may be optimized too, especially for automatic transmissions or CVT transmissions (continuous variation of the gear ratio), and u 3 is linked to x 2 state.

In this optimization problem, the fuel consumption is given by

L(x(t), u(t), t) =    L 1 (u 1 (t), t) if x 2 = 1 (engine started) 0 if x 2 = 0 (engine stopped) L 2 if x 2 = 0 and u 3 = 1 (engine start-up) (4)
where L 1 (u 1 (t), t) is the fuel consumption given by a map depending on engine speed ω e (t) and provided torque T e (t), and L 2 is the extra fuel consumption at engine start-up.

This minimization problem must be solved on a prescribed cycle, for instance the NEDC cycle, Figure 1. Thanks to a vehicle model (mass, wheel diameter, aerodynamic coefficients, etc.) and the speed request V (t), the requested torque and speed at the wheels are computed: (T w (t), ω w (t))). In this optimization problem, we consider the cycle as a wheel torque T w (t) and wheel speed ω w (t) trajectory to be followed. In this application, we consider that the engine is located upstream of the gearbox, as on conventional vehicles. Besides, the electric motor is located downstream of the gearbox, beside the wheels. Let us introduce now some useful notations to describe the powertrain and to define the control u(t) :

• the requested torque at wheels is written as T w (t) = R f ar R gr T e (t) + T m (t) + T bk (t), and can be achieved thanks to the engine torque T e (t) multiplied by the front axle ratio R f ar and by the gear ratio R gr , and thanks to the motor torque T m (t). A brake torque T bk (t) can also be applied by braking system.

• the engine torque is T e (t) = u(t)T w (t)/(R f ar R gr ),

• the electric motor torque is

T m (t) = (1 -u(t))T w (t).
During braking phases, T rq (t) < 0, the braking torque can complete the electric motor torque to slow down the vehicle: T bk (t) = T w (t) -T min m (t) with T min m (t) the minimal electric motor torque depending on the engine speed, thus depending on time.

A simple model of the battery is implemented, the state of charge of the battery evolving as

ẋ1 (t) = - ω m (t)T m (t)K U batt (t)n capa (5)
with ω m (t) the electric motor speed, U batt (t) the battery voltage, K a scaling constants, and n capa the nominal capacity of the battery. In the following, we assume that the battery voltage is constant along the cycle. We set K = K U batt .ncapa and by replacing T m (t) by its expression, we obtain

ẋ1 (t) = -Kω m (t)(1 -u(t))T w (t). (6)
In this optimal control problem, we make several assumptions about the powertrain: (i) the wheel speed ω w is equal to the electric motor speed ω m , (ii) as the optimization problem is applied on a full hybrid vehicle, the engine can be stopped when the torque is provided only by the electric motor. As mentioned above in the definition of L(x(t), u(t), t), a fuel consumption penalty, the necessary energy to start the engine, must be considered, (iii) recharging the battery is possible for negative torques (braking request) or by an additional engine torque beyond the driver request torque. In the last case, the command u 1 (t) is larger than 1, (iv) we consider that the average efficiency of the electric motor is 0.85.

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

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

• the electric motor torque is bounded:

T min m (ω(t)) ≤ T m (t) ≤ T max m (ω(t)),
• the storage capacity implies a minimal and a maximal state of charge of the battery (which are fixed to 50% and 70% in our study):

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

x 1 (0) = x 1 (T ).
Constraints on command u(t) are directly derived from the constraints on maximal and minimal engine and electric motor torque

u(t)T rq (t) ≤ T max e (t), (7)
T min m (t) ≤ (1 -u(t))T rq (t) ≤ T max m (t). (8)
In the following, we will note U (t) (respectively U k in discrete time) the feasible domain of u(t) (constraints ( 7) and ( 8)).

Dynamic Programming

The Dynamic Programming method (DP) is classically applied to solve the problem (1) ( [START_REF] Scordia | Systematic elaboration of online energy management laws for hybrid vehicles[END_REF], [START_REF] Wu | Optimization of Power Management Strategies for a Hydraulic Hybrid Medium Truck[END_REF]) : it relies on a very simple idea, the principle of optimality, stated by Richard Bellman. It states that an optimal policy can be built step by step, first by building an optimal control for the "tail subproblem" involving the last step, then by extending the optimal policy to the "tail subproblem" involving the last two stages, and by iterating until an optimal control has been built for the complete time interval. Thus, the optimization problem is solved backward from final time to initial time. The optimal control problem (1) is thus discretized in time:

min u∈U J(u) := N -1 k=0 L k (x k , u k )∆t + g(x N ), (9)
where

L k (x k , u k ) is the cumulated fuel consumption over the time interval [k, k + 1],
x k is the state vector at time k, and N is the final time of the driving cycle. The constraint on the final state of charge x 0 = x N has been introduced via a penalization term g(x N ) = β(x N -x 0 ) 2 (β, a constant to be chosen).

From Bellman principle, the minimum cost-to-go V k (x k ) for the time interval [k, N ], 0 ≤ k ≤ N -1, is then expressed as follows :

V k (x k ) = min u k ∈U (L(x k , u k )∆t + V k+1 (f k (x k , u k ))), (10)
where f k corresponds to the system dynamic, given by:

x k+1 = f k (x k , u k ) = x (1)k + f (1)k (u k , k) u (3)k (11)
f k being the function that modelizes the battery state of charge dynamic in the discrete form of (6), as well as the engine state. At time N , the cost function is

V N (x N ) = g(x N ). ( 12
)
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, that gives a 3-D grid (time t, state of charge x 1 , and engine state x 2 nodes). A standard time step used in our examples is ∆t = 0.5s, the discretization step ∆x 1 for the state x 1 is 0.25 %, and x 2 takes only 0 or 1 values. Thus, the number of state nodes is M = 2((x (1)max -x (1)min )/∆x 1 ).

A standard DP algorithm consists in computing, at each step k and at each state value, all feasible controls, and calculating the optimal control that minimizes the sum (

L k (x k , u k )+V k+1 (f k (x k , u k ))).
This method allows a fast optimization: (N -1)M constrained optimization problems [START_REF] Bonnans | Parametric optimization of hybrid car engines[END_REF] of only one parameter u k are solved with a fine time discretization of the controller. This provides the optimal power repartition u * k at each time step (we refer to [START_REF] Guilbaud | Méthodes numériques pour la commande optimale[END_REF] for some theoretical results on the convergence of the method and error estimations).

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

6. Variations of the size of some powertrain elements On conventional vehicles, the maximal engine power is mainly defined by its capacity. Economic vehicles have a small engine, with a small fuel consumption, but with a low maximal torque. When engine capacity is higher, fuel consumption is higher, but maximal engine torque is also higher.

In hybrid vehicles, engine capacity has still almost the same impact on fuel consumption. However, the presence of a second energy source, which implies at least an electric motor and a battery, gives other degrees of freedom to decrease the fuel consumption.

The following results have been obtained by testing different sizes of elements of the powertrain. This approach has been preferred to parametric optimization, which is a quite complicated problem. The reader can nevertheless see some examples in [START_REF] Bonnans | Parametric optimization of hybrid car engines[END_REF].

Reference case

The reference case is a Citroën C1 with a 1000 cc gasoline engine. We consider on this vehicle that an electric motor of 10 kW is located downstream of the gearbox, at the same level as the wheels, and that a battery of 1.3 kW/h provides the electric energy to the motor. The potential of this vehicle is assessed on the NEDC cycle, and gives the following results: The battery capacity affects the constraints on the state of charge: x min and x max . As the system dynamic does not depend on the state x, the fuel consumption will not evolve with the battery capacity as soon as the state of charge remains inside its bounds x min , x max without touching them. Table 2 gives optimal fuel consumption on NEDC cycle for varying capacities.

Variations of battery capacity and battery weight

By considering the variations of battery capacity and the corresponding variations of its weight, not only the available electric power is changed, but also the global weight of the vehicle. In this context, the necessary power to match the requested vehicle speed of the cycle becomes higher as the weight of the battery increases. Thus, the variation of battery capacity causes the engine to change its operating points: same engine speeds, but different requested torques. Also, the vehicle will potentially be able to recover much more electric energy during braking, if control constraints are not activated. Consequently, the analysis of these results is not trivial.

We consider in the following that the battery weight varies according to 1kg for 50Wh. The Table 3 shows obtained results. The car manufacturers are interested in estimating the best size of an electric motor for a hybrid vehicle, according to its use, as there is no need to have a very powerful motor if half of its maximum power is never reached.

The variations of maximal and minimal electric motor power may change the fuel consumption, but it also changes the state of charge trajectory (see results in Table 4). For this vehicle used on the NEDC cycle, the fuel consumption does not evolve anymore if the electric motor power is larger than P = 15 kW. Smaller powers of electric motor do not affect much the fuel consumption, as the gain for an electric motor is still 20.15% for a power of 5 kW. The global fuel consumption obtained with AMESim model is 310.2 g (3.81 l/100km) while the one with matlab model is 310.34 g (3.75 l/100km), but the vehicle linear displacement is not the same: Matlab model covers 11.01km, but AMESim model only covers 10.84km, because of inertia that cause a bad follow-up of the vehicle cycle at low speeds. The resulting error between the two models is about 1.5%.

Real-time control strategy

In this section, suboptimal control law is derived from the optimization results. The chosen realtime control law is the Equivalent Consumption Minimization Strategy: it is based on instantaneous equivalent fuel consumption, the battery being considered as a second energy source. ECMS aims to choose the cheapest energy source, between fuel tank and battery storage.

This strategy is derived from the Pontryagin's Principle, which set optimality conditions for the continuous optimization problem [START_REF] German | Hybrid powered vehicles[END_REF]. Some important theoretical results can be found in [START_REF] Pontryagin | The Mathematical Theory of Optimal Processes[END_REF]. To apply this strategy, we form the Hamiltonian function H given by H(x(t), u(t), p(t), t) = L(x(t), u(t), t)∆t + p(t) ẋ1 (t)∆t,

where L(x(t), u(t), t)∆t represents the variation of fuel consumption over a time step ∆t, given by

L(x(t), u(t), t) = L 1 (u(t), t) if x 2 = 0 0 if x 2 = 0 , (14)
ẋ1 ∆t corresponds to the variation of state of charge over a time step, and p is a Lagrange multiplier associated with the constraint ẋ1 (t) = f 1 (u(t), t). The state constraints have not been considered in the definition of (13). In this strategy, the optimal control is the admissible one which minimizes H, taking into account all the control constraints. The same models as the ones used in dynamic programming algorithm are used in ECMS, which is the advantage of this model-based strategy.

Meaning of p

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.

The variation of p depends on the time evolution of the dynamic state variable x 1 , according to the Euler-Lagrange equations, namely

ṗ(t) = - ∂H ∂x 1 (x(t), u * (t), p(t), t). (15)
As either L or ẋ1 do not depend on x 1 , Equation (15) gives

ṗ(t) = 0. ( 16
)
As a results, when no state constraint is active, p is constant for the considered system [START_REF] Guzzella | Vehicle Propulsion Systems[END_REF]. When a state constraint becomes active, the sign of ṗ is completely known, depending on the active constraint: x = x min or x = x max .

Control of p

In ECMS, the value of p is often controlled thanks to rules [START_REF] Delprat | Control of a parallel hybrid powertrain: optimal control[END_REF], [START_REF] Liu | Control Optimization for a Power-Split Hybrid Vehicle[END_REF], mainly depending of the state of charge. Indeed, as the future is not known, it can be better to keep the state of charge x 1 far from its bounds, and to control p as a function of x 1 with a feedback controller. However, in this paper, we consider that p remains constant, and we just try to show the potential of this strategy with an adequate initial value of p.

An estimation of p can be obtained by calculating the equivalent Lagrange multiplier with the dynamic programming results. From [START_REF] Bryson | Applied Optimal Control[END_REF], equivalence between Pontryagin's principle and dynamic programming can be found

p(t) = ∂V ∂x (t, x), (17)
where V (t, x) corresponds to the cost-to-go function in continuous time, already defined in [START_REF] Bonnans | Parametric optimization of hybrid car engines[END_REF]. Using theses results allows to estimate that p should be into the range [-3400; -3000]. Because of some differences between the model used with dynamic programming and the AMESim model, a trial & error method has been applied to end up with p = -3075 † . This value allows to reach a final state x(T ) close to the initial state x(0).

Results of the ECMS Real-time controller

The Figure 5 shows some results obtained with this real-time control strategy, with a correct value of p. shows that the hybrid vehicle still follows the speed request from the cycle. However, some differences appear on the state of charge, Figure 5(b). The fuel consumption obtained with this sub-optimal strategy is 326.8 g, corresponding to 3.97 l/100km. Compared to conventional vehicle, this controller -with an ideal value of p in this example -allows to save 16.6% of fuel on the NEDC cycle. Some rules to control p with the state of charge should be implemented to make this controller more robust.

Conclusions

In this paper, has been presented an optimization study applied on a full hybrid vehicle, with classical offline optimization tool, and an example of a real-time control strategy.

A simple model in Matlab has been used with a dynamic programming algorithm. This allows to get the optimal sizes and powers of elements of the hybrid powertrain. Nevertheless, these results are obtained on the NEDC cycle, so the optimal sizes can be a little different on other cycles.

After a validation of the optimal trajectory and torques applied on a more realistic model in AMESim, a model-based sub-optimal controller has been implemented, using Pontryagin's Principle and optimality conditions. Even if an average value of p can be estimated from offline optimization results, the choice of an adequate initial value remains a delicate task, and p should be controlled regarding to the evolution of the state of charge of the battery.

Figure 1 :

 1 Figure 1: NEDC driving cycle: vehicle speed at left / gear ratio at right

  capacity with constant weight variation of battery capacity and of its weight

Figure 2 :

 2 Figure 2: Variation of gain of fuel consumption with battery capacity 6.4. Variations of maximal and minimal electric motor torques

7 .

 7 Real-time control strategy 7.1. Validation of DP results on a more realistic model Results obtained with the dynamic programming applied on the matlab model are tested on a more realistic model implemented in AMESim, corresponding to the reference case, see Figure 3.

Figure 3 :

 3 Figure 3: Validation AMESim model

Figure 4 :

 4 Figure 4: Results of optimal control applied on AMESim model

  real-time controller Optimal SOC trajectory with offline optimization

Figure 5 :

 5 Figure 5: Results obtained with a sub-optimal strategy applied on AMESim model

Table 1 :

 1 Fuel consumption for the reference case

		Fuel consumption (l/100km) Gain (%)
	Conventional vehicle	4.76
	HEV	3.76	20.99 %
	6.2. Variations of battery capacity with constant battery weight

Table 2 :

 2 Study of the impact of the battery capacity

		Fuel consumption (l/100km) Gain (%)
	Conventional vehicle	4.76
	HEV with battery 0.325 kWh	3.78	20.40 %
	0.65 kWh	3.76	20.99 %
	1.3 kWh	3.76	21 %
	2.6 kWh	3.76	21 %
	5.2 kWh	3.76	21 %

Table 3 :

 3 Study of the impact of the battery capacity with weight variationsIt is noticeable that those results show a different trend than the results of Table2. An optimal value is reached for n capa = 1.04 kWh, as seen on Figure2.

			Fuel consumption (l/100km) Gain (%)
	Conventional vehicle		4.76
	HEV with battery	0.325 kWh	3.78	20.40 %
		0.65 kWh	3.76	20.84 %
		1.04 kWh	3.76	21.01 %
		1.3 kWh	3.76	20.99 %
		1.625 kWh	3.76	20.90 %
		1.95 kWh	3.76	20.79 %

Table 4 :

 4 Study of the impact of the maximal motor power

		Fuel consumption (l/100km) Gain (%)
	Conventional vehicle	4.76
	HEV with motor P= 5 kW	3.80	20.15 %
	HEV with motor P= 10 kW	3.76	20.99 %
	HEV with motor P= 15 kW	3.74	21.25 %
	HEV with motor P= 20 kW	3.74	21.25 %

† if the variations of fuel and chemical energies are used instead of the variations L(u(t), t)∆t and ẋ1 (t)∆t, the coefficient p can be scaled to a smaller value.


Fabien Chaudoye

Michel Castagne

Delphine Sinoquet

Francois Wahl

Michel Castagné

François Wahl

Modelling engine operating space for DoE calibration methods

come   L'archive ouverte pluridisciplinaire

Introduction

Nowadays, number of engine control parameters drastically increases. Consequently, DoE methods are now commonly used to reduce cost, delay and improve calibration process. Advanced methods [START_REF] Nelles | Non linear system identification[END_REF] based on modelling require some specific developments for several steps of the engine calibration: test planning, modelling, optimizing, etc. This paper will mainly discuss the first step of the calibration workflow, i.e. engine operating space (OS) determination. This step is crucial, since it allows to get the domain into which DoE could be built and control parameters could be optimized, with the insurance that every parameter combination will be performed correctly by the engine. Instead of an iterative approach such as Rapid Hull Determination [START_REF] Renninger | A new method to determine the design space for model based approaches[END_REF], we proposed to develop a physical and pragmatic approach to find the limits of OS. This will lead to basic OS modelling, easy to use during other calibration steps as modelling or optimization. This paper presents both local OS (with operating points defined by speed and load) and global OS (including speed and load as parameters) determination. Test plan problematic is also addressed. Purposes are illustrated with results from a EURO4 production 4 cylinders diesel engine operated at constant speed-load (stationary) in standard hot conditions on a fully instrumented test bed using Morphée as automation system.

Context

Most of industrial calibration methods based on DoE use simple domain definition, frequently a simple hypercube defined with minimal and maximal fixed bounds for every control parameter. The main advantage is a simple and fast determination, as well as an easy modelling of both engine domain and engine responses because of the limited variations of the control parameters. However, an important drawback of such domain is that it is significantly reduced in comparison with physical engine limits. Thus, depending on knowledge about engine behaviour, optimal parameter settings can be out of the defined hypercubic domain. Practically, the modelling space is limited to the domain used for DoE (validity domain of the model) in order to avoid extrapolation. Sometimes, the optimizer will find optimal settings at the limit of the domain. In such case, the full process (research of bounds, design of experiments, tests, engine responses modelling and optimization) has to be done again, until optimal parameters setting are into the limits of the domain ( [START_REF] Castagné | Engine calibration : Toward an integrated approach[END_REF]). Figure 1 shows this iterative process, which can be cumbersome and time-consuming.

Figure 1: progressive optimization

The aim of our study for local approach is to determine the complete OS with reduced test duration and to model the limits of this OS with simple mathematical expressions, easy to use in further operations as design of experiments or optimization. Such wide OS determination is very useful for several applications. First if we use the entire OS for designing the experiment, we can ensure that optimal settings can not be out of this domain and we avoid the iterative process described above. If a restricted domain is preferred for design of experiments, complete OS knowledge allows to define a new domain for iterative process without redoing preliminary tests. When using advanced calibration methods, large OS can be more easily used for merging local models to build a global model ([3]).

Local operating space limits

Limits of OS can either be hard limits as maximal air mass flow, minimal/maximal boost pressure, or limits defined with subjective criteria (AFR, opacity, stability, etc.) as minimal air mass flow, or minimal main injection advance. Figure 2 shows the limits of operating space in the 2D space (air mass flow and boost pressure), for given injection parameter setting in comparison with a usual hypercubic reduced domain designed into the OS. Experience shows that physical engine air loop limits could be modelled with straight lines for fixed injection parameter settings. The maximal boost pressure directly depends on turbocharger characteristics and evolves linearly with air mass flow. The upper air mass flow limit depends on global engine permeability. Thus it mainly evolves linearly with boost pressure, and also varies slightly with engine speed and load. This limit is not very sensitive to injection parameters. The lower boost pressure limit is related to turbocharger permeability, and also depends on atmospheric pressure. The last air loop limit is the minimal air mass flow. It is a subjective limit unlike previous described limits which are physical limits. Thus, it mainly depends on criteria chosen to define it. AFR limit is often use for this purpose. Linear models also present good results for this limit. However, if OS limits seem to be easily determined with a simple shape in a 2D space, it becomes more complex with an increasing number of parameters. Figure 3 shows OS shape in a 3D space air mass flow, minimal/maximal boost pressure/main injection advance.

Figure 3 : example of enlarged OS in a 3D space

With one more parameter, OS general shape is significantly modified. Indeed, every limit is impacted by other parameters. For example, when main injection advance is decreasing, maximal admissible boost pressure is increasing since there is more energy contained in gas for turbocharger. In another way, main injection advance can be reduced when air mass flow increases (i.e. EGR rate decreases) since combustion is more stable. Thereby, with many parameters, OS shape is much more complex since there are various interactions between parameters.

The pragmatic approach presented in this paper only considers relevant physical limits. For example, it is assumed that there is no need to widely increase main injection advance for low load operating points, although it is not limited with too high cylinder pressure, but due to high NOx emissions and noise value. Thus, this approach tries to define engine limits as wide as possible, by introducing engineering knowledge to keep most relevant part of the OS.

Test facilities and procedure of OS determination

All tests presented in this paper were performed with a 1.6L diesel engine equipped with a common rail injection, a variable geometry turbocharger and an EGR cooling system. The engine control is mainly based on maps. All tests were carried out with a simple injection strategy, i.e. a single pilote injection besides the main injection. Thus, six engine control parameters were used to define OS limits: air loop parameters, including boost pressure and air mass flow, and injection parameters including injection pressure, main injection advance and pilote injection advance and quantity. The test cell used during tests is equipped with Morphee [4] as cell supervisor, including Automapping component to design tests. This configuration makes it possible to perform tests in an automated way without the presence of operator, which drastically improves the productivity of the global process. The cell is also equipped with several measurement systems: engine and fluids temperatures, cylinder pressures, cell data, pollutants analyser, smokemeter and fuel consumption measurement. The OS determination is experimentally performed in two steps.

1. The first step is executed manually and consists of the determination of injection parameters bounds. Most of these bounds are considered as fixed. These are determined to allow a wide and relevant (considering pollutants emissions, noise, consumption, etc.) range of settings. It also consists of the determination of a criterion to define minimal air mass flow limit. 2. The second step consists of looking for air loop limits for given injection parameters.

As explained before, air loop limits are straight lines: thus, two points are enough to model limit for given injection parameters.

Specific developments were done to determine minimal air mass flow limit. This limit is gradually determined with continuous measurements (opacity, AFR or engine stability) until a bound defined at the first step is reached. Figure 4 shows principle of automatic test and iterative determination of minimal air mass flow limit. In order to shorten time spent in OS determination, no stabilisation is expected during test, since pollutants, temperature and other cell data are not analysed. In this way, about one minute is required to find one point at the maximal air mass flow limit, and about two minutes to find one point at the minimal air mass flow limit, due to iterative process described previously. Thereby, it is possible to perform OS determination for each operating point in a short duration. Typically six minutes are needed to define air loop limits for a given injection parameter settings, and a standard of fifteen combinations of injection parameters are usually performed. Thus, around three hours (including transitory states between combinations) are necessary to define a complete OS.

Figure 4: Principle of automatic test to determine air loop limit with iterative method

The minimal air mass flow limit is not the only one to be determined with an iterative process. For example, lower limit of main injection advance is also determined with such a process. When air mass flow is increasing, the combustion is more stable, and CA50 (crank angle where 50 percent of burnt fuel) which is a relevant combustion criterion increases. Thus, it is possible to reduce main injection advance. Criterion to define lower main injection limit is indeed the CA50 equivalent of the value obtained for minimal air mass flow.

Local operating space modelling

OS model has to meet several requirements. It has to be simple enough to be easily used during the next step of the workflow: design of experiments. And it also has to be accurate enough, to ensure that every set of parameters defined for the test is physically realizable. After processing the experiments, we can model the constraint area in all directions. The main idea is to model the maximal and minimal limits for each variable: Results are fairly good for the three physical limits: lower and upper boost pressure limit and maximal air mass flow. Minimal air mass flow results present a worse accuracy, since this boost pressure air mass flow limit is subjective, and directly depends on the way it is defined. However, accuracy obtained with these models is sufficient to ensure feasibility of all the parameter settings. Figure 7 presents several constraints projected on 2D planes for a given OP. On the left, air loop limits are presented for several main injection advance and rail pressure values. As described in section 2.2, OS limits depend on these parameters: trends are well reproduced by the models. On the right, influence of air mass flow on main injection advance is represented: when air mass flow is increasing, combustion is more stable, thus, minimal main injection advance limit is decreasing. Otherwise, when main injection advance is decreasing, there is more energy available for VGT, and for a given boost pressure, minimal air mass flow is increasing.

) , , ( ) , , ( 1 max

Figure 7: example of 2D projection of local OS model

The advantage of this modelling clearly appears in comparison with hypercubic approach, as presented in Figure 8. For a given OP, standard hypercubic domain and expanded OS are compared. The OS enlargement avoids hypercubic OS drawbacks, i.e. the possibility that the optimum lies outside the OS. However, this expansion requires modifying test plan methods. This point will be discussed in next section.

Test plan development

The aim of the test plan is to build response surfaces approaching as much as possible the real behaviour of the engine in the chosen local area. In the remainder of this section, we start from a simple design in a hypercube that we make more complex step by step to finally designing the experiment in the entire OS. If domain was hypercubic, strategies for designing experiments would be well established. Usually in this kind of setting, the area to explore is relatively narrow in comparison to the maximal possible variation of each parameter, as it has been noted a couple of times. Thus, a polynomial representation of the studied engine response (CO2 emission for example) performs locally well. In this case, for six engine control parameters, and a quadratic polynomial, 28 coefficients have to be determined, meaning that at least 28 experimental points have to be gathered. However, it is wise and good practice to add some points within the domain and to keep an

AVI=-4°V RP=800bars AVI=+4°V RP=1200b

other set of experiments for validation. Indeed, we need 40 points as a minimal. To choose the experimental trials, D-optimal designs [START_REF] Atkinson | Optimum Experimental Designs[END_REF] can be put into practice. But even in a hypercubic domain, the real behaviour of the responses may be not so comfortable. In certain cases (especially for particles representation), we need a cubic polynomial to take into account all possible variations. Neglecting all the order three interactions, 34 coefficients are now to be determined, and the number of needed experimental points increases accordingly.

If we still want to refine our model, it can be noticed that the shape of the answers will depart more from the expected model as we approach the edges and the corners of the domain. So it sounds like a good strategy to include all the corners in the design to potentially avoid important residues for these areas. In this case, still within a six dimension parameter space, there are 64 corners to be tested. Furthermore, added points inside domain are still needed. Now, it appears that the real OS is not hypercubic but is limited with linear constraints, which are not parallel to the axes. Thus, basic hypercube becomes a complex polyhedra, which can present sharp and thin edges and the number of possible corners increases. For example, we had to treat an example with six additional constraints, and 89 corners were counted. A way to valid engine responses model consists of comparison between model prediction and simple experimental trends. For instance, experts build profile variations of the answers, meaning that all the parameters are fixed except one, which is varying between the allowed limits. Practically, a possible idea is to include in the test design a few points dedicated to represent the central profiles. Then, starting from the central point (P) of the domain, two points are taken on both sides of P along each direction. Thus 24 experimental trials are added. At the end, we have around 110 imposed points for a six dimensions case, counting all the corners and the experiments for the profiles. Besides, some few points are defined inside the OS by some space filling algorithm, and some security checks to be sure that all the experimental apparatus does not show any deviation during our experimental campaign, we eventually get 140 experimental trials to run. Figure 9 gathered all experimental points planed for a given speed and load operating point, in the 2D plane boost pressure and air mass flow.

Operating space definition: global approach

This section describes global OS modelling method, i.e. validity domain determination including speed and load as parameters.

Context

The time and cost constraints imposed by the automotive industry, lead to ever shorter product development. Thus, advanced calibration methods answer to these constraints, and aim to significantly reduce the time development of vehicles. Global approach for calibration, considering load and speed as parameters, reduces time spent for whole calibration workflow, since a unique DoE is necessary to model engine responses for entire domain (or several DoE if this domain is divided in several zones with common parameters). However, problematic is more complex than for local approach, and specific development were realised. Once again, OS determination is an important step for this calibration approach, and directly influences quality of global DoE and potential of optimization. For this last step, it is necessary to ensure that bounds for every parameter are wide enough to complete objectives of calibration campaign. Good OS models are also suitable for test realisation. Indeed, for mid and high load and speed zone (considering normalised cycle load and speed zone), many criteria are closed to engine limits, like intake manifold temperature, temperature before VGT, etc. If OS models are not precise enough, many settings of parameters planed by DoE will be out of OS leading to test bed stops, significantly lengthen test development. Figure 10 presents OS evolution (2D plane Air mass flow / Boost pressure) with load and speed. The aim of the approach presented in this paper is to get model of the OS boundaries. Such a model has to keep simple mathematical expression to be compatible with test plan design and optimisation methods.

Global operating space modelling

The OS modelling has to fulfil several requirements. It has to be compatible with simple DOE, model has to be simple enough to be used for DOE tools and it also has to be precise enough to ensure good development of test. Pragmatic approach was favoured and different types of modelling were defined for several parameters. Local approach and its very satisfying results, significantly inspire global approach method. Considering load and speed limits, simple linear models delimiting normalised cycle working area were used. For injection parameters, lower and upper limits are defined with maps. Amplitude of variations is defined thanks to engineering experience, and has to be compatible with calibration objectives. This type of model was chosen from local approach experience where fixed bounds for injection parameters led to satisfying results in term of potential for optimization. These models are also mathematically simple, since limits come down to simple bounds for fixed speed and load. Figure 11 presents example of limit maps for injection rail pressure and maximal main advance injection.

Figure 11 : example of limit modelled with cartography

As in local approach, the most critical point for global OS is air loop limits. Two methods were tested. The first method consists in using local data gathered for every operating point, and in modelling limits with low degree polynomials. Some interaction terms of the polynomial models were removed in order to obtain linear models for fixed speed and load, to be similar with local approach. This first type of model led to satisfying results, however, a lack of data within critical zone in mid and high speed and load drove to a poor precision of the model in this zone. Consequently some problems appeared during DoE test. The second method is based on a dedicated automatic test covering entire normalised cycle area. This second method presents the advantage of getting data well distributed all over cycle area, and not only for a limited number of operating points. This is particularly important for precision, and it also reduces overall test duration. The difficulty consists in defining a criterion to find minimal air mass flow limit for several speed and load combinations. A model of AFR limit, defined from local approach experience, has been used as criterion. In a limited time, it is possible to gather enough information to model global OS. Figure 12 and Figure 13 present modelling experimental correlation diagrams for the two limits maximal and minimal air mass flow. Crosses represent point used for modelling, and circles represent validation points. Results of maximal air mass flow limit are very good, while minimal limit presents worse precision, but still acceptable. This requires to define correctly global AFR criteria in order to ensure that planned parameter settings are physically realisable. The correlation diagram presents some high residues for air mass flow value around value of two hundred milligrams per cycle. This is due to permeability limit of EGR system, which does not allow to reduce air mass flow for very low speed and load OP. Model does not well represent this trend, with no consequences about engine security since air mass flow is physically limited in this area and no unfeasible point could be performed.. The next figure presents comparison for air loop limit between local and global OS models for a given operating point. As expected with good results presented previously, global OS is similar with local one, which validates global OS model. Furthermore, duration dedicated to this global approach is very encouraging and meets industry requirement. Test planning methods were modified to take into account these more complex models; this point is discussed in next section.

Global test plan

Design of experiments in this global setting rests on the same considerations as in local one. At the beginning, the space to be covered is a hypercube, which is limited by a set of constraints as it has been described in the previous sections. But due to the increase of dimensions, eight instead of six (including speed and load), new difficulties arise. Some constraints are not linear, due to the failure of simple models to respond to our objective: much more experimental points are requested to adequately represent the domain. When all the constraints are non linear, the basic strategy for designing the test plan described above has to be modified. Indeed, it is no longer possible to define the corners, as the edges of the domain can be curved. In our case, experts had to resort to quadratic polynomial to model the limits of the OS. The proposed approach consists in proceeding in two hierarchical steps, considering local approach experience. Indeed, for a given speed and load, the constraints delimiting OS have to be linear. Thus in a first step, a set of constraints is built over the speed and load domain. In this two dimensions space, we proceed as usual: the basic square (in two dimensions) is limited by constraints which have to remain linear in order to ensure that corners can be defined. All other constraints which have to be applied on the other parameters can depend nonlinearly on speed and load. However, once the speed and the load are fixed in a second step, these expressions must degrade to a linear convenient form for the six other parameters. The way an experimental point is chosen is quite simple in this way: firstly, a couple of values for speed and load is defined, and secondly, six other parameters are fixed as in the local test plan method. In our case, with the set of constraints we had to deal with, we found around 400 corners. We added three set of mono parametric variations, which count for about 100 more points. So 500 points are imposed. As we stand in a eight dimensional space, it is not enough to describe it fully: indeed, no specific parametric form of model seems adequate and one has to resort to kriging techniques and space filling design [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. If two points per dimension are chosen, the whole test design would be composed of 256 experiences. With three points, an overall of 6561 experiences is requested, which can not be afforded. At the end, we compromised for 2000 points (including corners), from which 200 were left out the calibration set and kept for validation only. To fix these 1500 imposed points, the two steps described above were run iteratively, by random sampling. Clearly, the overall design is satisfactory in the sense that all the prerequisite are fulfilled. But as can be seen in Figure 15, this design is not optimal, some holes and some highly concentrated regions are present. This figure shows also that the limits are non linear. The final step, once all the requested trials are run, is to build the model, which is done by kriging on 1800 calibration points, the validation is done on the set of 200 points.

Boost pressure

Air mass flow

Enhancements and outlooks

Method to determine local OS could be improved with a reduced duration dedicated to this step. Time could be saved with optimization of experiments performed during tests. Indeed, about 50 experiments are carried out to determine all linear constraints, including some validation points. This number can probably be reduced to 40 with experience and optimization of design of experiments. Then, test plan method can also be improved. If a wide OS leads to an increasing number of experiments to be done, experience will allow to reduce this number to improve productivity of calibration process and especially test bed experiences.

Besides, methods have to be tested with more parameters, since their number will sensibly increase in the next few years. However, approaches presented seem to be compatible with more injection parameters, and air loop limit determination method will still be valid with a low increasing of time and cost spent for this task.

An important way of enhancement is coupling experimental tests with simulation. Indeed, physical 0D or 1D simulation of the engine could help for determining both local and global OS boundaries. For example upper limits of air loop can be easily determinated by the use of a simulator. This would significantly save time at test bench and cost spent for this step of calibration workflow.

Conclusion

This paper presents simple and practical methods for determining the limits of engine operating space in order to use DoE calibration methods. Local approach and global approach are both addressed. Experimental test procedures, as well as specific tools have been developed to model these limits. These procedures are run automatically in order to save time and use a test bench 24h a day. Simple modelling of the operating space limits has been set up from the experimental data. Results are satisfying for both approaches in terms of prediction and time spent for tests. These models are used to design the experiments but also as constraints during optimization process [START_REF] Sinoquet | Multiobjective constrained optimization of engine maps[END_REF].

Test planning methods were successfully developed into the modelled operating spaces, in accordance with the type of engine response model anticipated. These methods have now to be adapted for an increasing number of parameters.

Figure 2 :

 2 Figure 2: limits of operating space and example of hypercubic domain in 2D space (air mass flow / boost pressure)

1 ==

 1 y is a given variable and ( ) n i i p are the other variables. As shown in Figure2, OS shape seems simple, thus, linear constraints were used first, and expressed as followed: are estimated from experiments of the specific space design algorithm implemented on the test bench. In the case of hypercubic limits, we have: Figure5and Figure6present modelling experimental correlation diagrams for the four air loop limits. Dots represent points used for modelling, and circles represent validation points which are not used for modelling, and thus are used to judge model prediction capability.

Figure 5 :

 5 Figure 5 : correlation diagram for three physical limits

Figure 6 :

 6 Figure 6: correlation diagram for minimal air mass flow limit

Figure 8 :

 8 Figure 8: comparison of local hypercubic domain (left) and enlarged OS (right)

Figure 9 :

 9 Figure 9: example of local DoE in 2D space projection

Figure 10 :

 10 Figure 10: evolution of local OS with load and speed

Figure 12 :

 12 Figure 12: correlation diagram for maximal air mass flow limit (global OS)

Figure 14 :

 14 Figure 14: comparison of local and global OS model for given speed and load

Figure 15 :

 15 Figure 15 : 2000 experimental points chosen in the Prail (rail pressure) and air flow rate (Qair) space

Acknowledgments:

The authors would like to acknowledge Y.Bentolila, S.Gindre, J.M. Hascoet, F. Nicolas and R.Dauphin for their contribution to this work.

The authors:

Fabien Chaudoye, research engineer in engine development Dr. François Wahl, project leader in applied mathematics Michel Castagné, project leader in engine and vehicle calibration Dr. Delphine Sinoquet, project leader in applied mathematics All the authors work at IFP, 1 et 4 avenue de Bois Préau, 92582 Rueil-Malmaison Cedex -France.


Multi-objective constrained optimization of engine maps

Delphine Sinoquet, Hoël Langouët, Fabien Chaudoye and Michel Castagné

A bstract

Nowadays, automotive manufacturers are submitted to strong constraints in engine calibration: lowest fuel consumption, emission-control legislation and driver requests for driving comfort and performances. These constraints lead to an increasingly complexity of the engines and thus an increasingly 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 parameter maps that minimize the cumulated fuel consumption and pollutant emissions, under combustion noise constraints, on a driving cycle. The usual way to get this result is to select specific operating points representing this cycle in the engine working range and to define upper bounds applied on the different engine responses (allocations) for each of them, in order to obtain a weighted sum of these local responses respecting the global targets. The underlying problem is a multi-objective optimization problem: different compromises between fuel consumption, noise and pollutant emissions on each operating point are possible. We propose an adapted optimization method based on the MO-CMA-ES method (Multi-objective Covariance Adaptation Evolution Strategy) which takes into account the non trivial limits of the engine parameter variation domains and some robustness constraints. An other point addressed in this paper is the map optimization which consists in a global optimization of engine responses cumulated on the driving cycle. This method avoids the cumbersome choice of allocations for each considered operating point and includes directly the map regularity constraints in map parameterizations. Finally, application on real dataset obtained at automated test-bench for a diesel engine are presented.

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 ( [START_REF] Rö Pke | DoE Design of Experiments, Methods and applications in engine development[END_REF], [START_REF] Isermann | The use of stationary and dynamic emission models for an improved engine performance in legal test cycles[END_REF], [START_REF] Seabrook | A Calibration Methodology to Ensure Low Variation in Emissions between Engines in Production[END_REF], [START_REF] Castagné | Comparison of Engine Calibration Methods B ased on Design of Experiments[END_REF]). 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 adapt the Multi-Objective Covariance-Adaptation Evolution Strategy method, proposed by [START_REF] Ostermeier | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF], 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.

E ngine calibration

Sketch of the engine calibration process

The emission calibration workfl ow is classically divided into three steps ( [START_REF] Rö Pke | DoE Design of Experiments, Methods and applications in engine development[END_REF]):

1. a preliminary phase consists in choosing a sample of operating points (referred to as OP in the following) to be studied and emissions targets associated with each OP, these targets being called allocations hereafter. [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[END_REF]. 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). 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 (New European Driving Cycle) simulation and a selection of 16 OP. The cumulated levels of pollutant emissions along the cycle are computed as weighted sums of the pollutant emissions for each chosen OP. Note that the transient effects during the cycle as well as the effect of the warming up of the engine and of the after-treatment on the emissions are not taken into account in these weighted sums. They must be introduced when defining the objective of optimisation (the allocation) for each OP.

Phase 2 consists of five steps: i) defining the domain of variations of the engine control parameters ( [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[END_REF]: 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 are usually sufficient to accurately model engine responses. H owever, choosing too small domains leads to difficulty in coherently fulfilling sub-optimal engine maps ( [START_REF] Castagné | Comparison of Engine Calibration Methods B ased on Design of Experiments[END_REF]).

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 ([10]). iii) running this test matrix on the test bench: as the tests are predefined, the experiments 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. 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 experiments. For some response as particulates or CO emissions, this type of models may be too limited, then non-parametric response surfaces may be chosen as RB F or kriging ( [START_REF] Seabrook | A Calibration Methodology to Ensure Low Variation in Emissions between Engines in Production[END_REF]). 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é | Comparison of Engine Calibration Methods B ased on Design of Experiments[END_REF] or in [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[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 on the whole engine operating domain (see Table 2 for examples of engine maps). 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 remain as close as possible to the local optima while preserving a smooth shape of the map, 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 4.

T he O P optim iz ation problem form u lation

The OP optimization problem consists in minimizing some engine responses under constraints on other 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 concerns the NOx/particulates trade-off. It could be formulated as the minimization of the particulate emissions and the NOx emissions under constraints on fuel consumption (CO 2 emissions), on CO and H C emissions and engine noise level. G enerally, 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 Sequential Q uadratic Programming method run several times with different initial points to search for the global optimum.

S olving th e O P multi-objective optimization problem

In this section, we present a method to solve the constrained multi-objective problem (1). It is based on the algorithm proposed by [START_REF] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF] namely the Multi-objective Covariance Matrix Adaptation Evolution Strategy method (MO-CMA-ES). The MO-CMA-ES method belongs to the class of evolution strategy methods which allow to solve global optimization problems. These methods are popular for their ability to find a global minimum but their slow convergence rate is often criticized. In [START_REF] Ostermeier | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF], H ansen 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.

.1 C M A -E S m ethod for single-objectiv e optim iz ation

In the CMA-ES evolutionary algorithm, a new population is created from the current population by generating realizations of a G aussian 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 to be minimized. The next population is thus directed towards the best points.

An ideal choice for the covariance matrix would be the inverse of the H essian matrix of the objective function, the idea proposed by H ansen is then to approximate this matrix as for a Q uasi-Newton method but without any computation of the derivatives. The update formula for the covariance matrix is then

C (g+1) = (1 -c co v )C (g) + c co v µ 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 H ansen 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] Ostermeier | 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 toward the global minimum.

.2 M O -C M A -E S m ethod for m u lti-objectiv e optim iz ation

Adapting the CMA-ES method to multi-objective optimization requires the definition of new criteria to define what is a g ood and a bad individual for selection step of the evolutionary algorithm (see [START_REF] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF]). Two criteria are implemented in MO-CMA-ES:

n do m(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:

n do m i (A) = n do m(do m i-1 (A)) with do m i (A) = do m i-1 (A)\n do m i (A) and do m 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] Igel | Covariance Matrix Adaptation for Multi-objective Optimization[END_REF].

.A pplication of M O -C M A -E S m ethod 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 rpm -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 6 engine control parameters: the main injection timing, the pilot fuel injection quantity, the pilot injection timing, the fuel injection pressure, the mass air fl ow and the boost pressure.

The optimization objectives are the particulate emissions and the NOx emissions, CO2, H C 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 G aussian 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 two objectives that have been obtained by the multi-objective optimization method. The points indicated with dots 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 points not located on the Pareto front: these 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.

S olving th e map optimization problem

In the classical calibration process described in section 2.1, 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 p 1 (r(t), c(t)), m p 2 (r(t), c(t)), . . . , , m p Np (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 p 1 (r(t), c(t)), m p 2 (r(t), c(t)), . . . , , m p Np (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 p i are the 2D engine maps of the control parameters in the engine (speed,load) operating domain.

The objectives to be minimized (or constrained) are the engine responses cumulated on the considered driving cycle: the cumulated responses may be either the weighted sums of the local models defined at chosen representative OP or the integral of global models defined on the whole operating domain (depending on engine speed and engine load or some models defined on zones depending on the parameterizations: for instance number of injections). The latter formulation allows a finer optimization with a fine sample of the integrals in [START_REF] Matheron | Principles of geostatistics[END_REF]. 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. Note that, as with local optimization problem, the transient effects are not taken into account in engine response modeling and must be introduced when defining the objective of global optimization.

.1 M od eling the engine m aps

The formulation (8) requires an adapted parameterizations of the engine maps that must be fl exible 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 ( [START_REF] Nelles | and other modern development methods[END_REF], [START_REF] Isermann | 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 fl exibility, 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 (9)
where the weights Φ i (r, c), normalized G aussian 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 parametrization, namely the patch definition for LoLiMoT description of the maps should refl ect 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.

Φ i (r, c) = µ i (r, c) M j= 1 µ j (r, c) with µ i (r, c) = e x p - 1 2α (r -r i 0 ) 2 σ r i 2 + (c -c i 0 ) 2 σ c i 2 (10)
Thus, from some reference engine maps or some a priori information, a LoLiMoT parametrization 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 in equation ( 9)).

.2 A pplication of m ap optim iz ation on a real case

The map optimization has been applied on a turbo-diesel engine application. The optimization problem is defined as a constrained single-objective optimization problem where the cumulated CO2 emissions are minimized and the cumulated particulates, CO, H C and NOx emissions are constrained to remain under thresholds. Constraints are added to take into account the limits of the physical domain: linear and non linear inequality constraints as described in [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[END_REF]. Different methods have been compared as described in Table 1: • the loc a l m eth od is the most simple one but requires a smoothing step to build engine Table 1: Engine map optimization methods to be compared

m p ∈R Np N OP l= 1 ω l F (l) i m p 1 (r l , c l ), m p 2 (r l , c l ), . . . , , m p Np (r l , c l )
Mix ed g lob a l m eth od g lob al/z on al models weig h ted su m of g lob al resp on ses

(k rig in g ) m in m p ∈R Np N OP l= 1 ω l F i m p 1 (r l , c l ), m p 2 (r l , c l ), . . . , , m p Np (r l , c l )
G lob a l m eth od g lob al/z on al models in teg ral over th e cy cle of g lob al resp on ses

(k rig in g ) m in m p ∈R Np T 0 F i r(t), c(t), m p 1 (r(t), c(t)), . . . , , m p Np (r(t), c(t)) dt
maps. This step may destroy the optimization work by modifying the obtained local optima at a set of 16 selected OP as explained in Section 2. [START_REF] Castagné | Comparison of Engine Calibration Methods B ased on Design of Experiments[END_REF]. The engine maps are sub-optimal maps because of the smoothing step but also because of the limited number of selected OP.

• the m ix ed loc a l m eth od allows to build optimized engine maps by a global map optimization on cumulated responses on selected OP, responses being modeled only locally on each selected OP. The models are cubic polynomials built from data obtained thanks to a D-optimal design. This method may suffer from a lack of accuracy because of the limited number of OP as for the local method but, on the contrary to the latter, this method avoids the delicate smoothing step.

• the m ix ed g lob a l m eth od is very similar to the mixed local method with responses modeled globally on the operating domain. We use kriging models ( [START_REF]A statistical approach to some mine valuation and allied problems on the Witwatersrand[END_REF], [START_REF] Matheron | Principles of geostatistics[END_REF]) that allow to model complex responses and to take into account an estimate of measurement errors (see [START_REF] Chaudoye | Modelling engine operating space for DoE calibration methods[END_REF] for more details on the design of experiments). This method is not very interesting in practice: since we have at our disposal global models, an optimization of a finer sampling of the driving cycle is possible and should give more accurate results.

• the g lob a l m eth od is the most complex method but the most accurate: it consists in the map optimization of the cumulated responses on a fine sample of the cycle (in our example, 1H z sampling i.e. 1019 OP, see Figures 12 and13). The models are the same as in the mixed global models.

As shown by Figures 7 and8, the four methods match the constraints on the cumulated engine responses. Some of the engine maps obtained with these methods are displayed in Table 2: we observe that the maps obtained with mixed methods and global method are smoother than the maps obtained with local method. This is not surprising: the smoothing step of the local method is cumbersome, a compromise between smoothing and keeping the maps close to the local optima has to be found. A first advantage of the engine map optimization is thus the additional degree of freedom introduced by optimizing directly the cumulated responses instead of successive optimizations of local responses at each OP as in the local method. A response at an OP may be slightly increased in order to keep the regularity of maps. We observe the gain on Figure 7: the mixed local method gives smoother engine maps (Table 2) and allows to obtain smaller CO2 emissions.

Finally, the estimate of cumulated responses on the driving cycle should be more accurate: as shown on Figure 9, the approximation of the integrals by a weighted sum on 16 OP gives very different results from the 1H z discretisation. The choice of the 16 OP and the associated weights is cumbersome and may explain partly the differences. The accuracy of these approximations depends also on the accuracy of the models of the engine responses. The global models of the engine responses are not as accurate as the local models when considered locally. A special care should be taken in the validation of the global model.

Figure 8 illustrates the results obtained with mixed global method and global method: the method differs only by the discretisation of the driving cycle for computing the cumulated responses to be optimized (discretised with 16 OP or with 1019 OP for a 1H z sampling). The results are very similar: for this example, the mixed global method does not suffer from the rough sampling of the cycle. Nevertheless, the cumulated CO2 emissions are smaller for the maps obtained by the global method. Figures 10 and11 show the gain on the maps of the responses: the global method minimized more globally the CO2 emissions than the mixed global method.

O utlook

The methods proposed in this paper are based on experimental tests in stabilized conditions. The transient effects during the cycle as well as the effect of the warming up of the engine and of the post-treatment of the emissions are not taken into account during modeling process. They are only introduced in the definition of the objectives of optimization (the allocation for each OP in local approach or the objectives of map optimization in mixed or global approaches), with empirical methods based on the know-how of engine experts. These methods consist generally of minimizing each emission target by transfer coefficient between hot and cold behavior and between stabilized and transient behavior. Such coefficient are often determined by experience on similar engines or preliminary tests.

New challenge is to develop methods allowing to take into account in a more robust way the behavior of the engine during transient operations (load, speed and/or thermal variations). Experimental methods are in development, based on the realization of transient operations at test bench, for example by using DoE. Such method could be proposed in addition to methods based on stabilized tests (map corrections to take into account transient effect) or be used directly to obtain optimal maps. Some others methods proposed in the literature are based on generating APRB S signals (Amplitude modulated Pseudo Random B inary Signals) at each operating point for the identification of the non-linearity of the engine ( [START_REF] Martini | Dynamic Tuning: A new Approach for Fast Engine Emission Calibration[END_REF]). We also consider that coupling between experimental (statistical) models and physical modeling of the engine is a very promising way. The idea is to simulate the behavior of vehicle and engine in order to get the conditions in the combustion chamber before combustion(gas composition, injection behavior,...) and to use statistical models to obtain engine-out emissions. Such a simulation tool could be used for optimizing control parameters on the cycle, taking into account transient behavior.

C onclusions

In this paper, the optimization problems encountered in stationary 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, an other 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 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 compared to local methods. A comparison of the different approaches on this example shows the efficiency of the map optimization. A next study will be to apply a multi-objective approach as for the OP optimization, the difficulty being the large size of the parameter space (500 parameters to describe all the engine maps). Taking into account transient behavior of the engine during optimization process by the use of dedicated tests or by the use of simulation tool is an other challenge for calibration engineers.

Figure 1 :

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

Figure 2 :

 2 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 G aussian probability law. The ellipse 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.

  1. A classical criterion for multi-objective optimization is the Pareto dominance: Definition 1. A n ind iv id u al x is said to d om inate anoth er ind iv id u al x , w h ich is sym boliz ed by x ≺ x , if and only if ∀m ∈ {1, ..., M} : f m (x) ≤ f m (x ) ∃m ∈ {1, ..., M} : f m (x) < f m (x ).

( 6 )

 6 Definition 2. T h e set of Pareto op tim al p oints of th e set X is th u s d efi ned by th e set {x | x ∈ X : x ≺ x}.

Figure 3 : 10 - 1 2

 3101 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 (stars) obtained with the MO-CMA-ES method for the OP 1750 rd/mn -1 bar. 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 dots 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 parameterizations 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.

Models O p tim

  iz ed r esp on ses L oc a l m eth od local models local resp on ses (3 rd deg ree p oly n omials) for each O P , m in p∈R Np (F 1 (p), F 2 (p), . . . , Fm(p)) Mix ed loc a l m eth od local models weig h ted su m of local resp on ses (3 rd deg ree p oly n omials) m in

Table 2 :

 2 Examples of engine maps obtained with the different optimization methods. The green dots are the OP used to sample the integral of the cumulated responses (16 OP pr 1019 OP depending on the optimization method).

  al me th o d M ix e d lo c al me th o d M ix e d g lo b al me th o d G lo b al me th o d

Figure 7 :

 7 Figure 7: Comparison of results obtained by local method and mixed local method. The cumulated responses are computed by a weighted sum of 16 OP with local models. Upper bounds of engine responses introduced in optimization are displayed in dotted red lines.

Figure 8 :

 8 Figure 8: Comparison of results obtained with mixed global method and global method. The cumulated responses are computed by a 1H z sampling of the integral of global models (1019 OP). Upper bounds of engine responses introduced in optimization are displayed in dotted red lines.

Figure 9 :

 9 Figure 9: Comparison of estimation of cumulated responses on the driving cycle with a weighted sum of 16 OP and the integral with a discretisation of 1H z (1019 OP).

  Figures 12 and 13 display the optimized engine settings and the gain on engine responses obtained with global map optimization.

Figure 10 :Figure 11 :

 1011 Figure 10: Maps of differences of engine responses (modeled by kriging global models) obtained with initial maps and maps obtained with mixed global methods. The white dots indicate the driving cycle OP.

from top to bottom, the injected air quantity, the fuel injection pressure, the pilot fuel injection quantity, the main injection timing, the pilot injection timing, the position of variable geometry turbine, the boost pressure.