Stratigraphy L Pan

C Quantin

REGIONAL GEOLOGICAL CONTEXT OF THE INSIGHT LANDING SITE FROM MINERALOGY AND

Introduction: Planned to launch in May 2018, the InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport) robotic lander will explore the interior of Mars with a set of geophysical instruments [START_REF] Banerdt | LPS XLV Abstract #[END_REF][2]. Orbital data for the InSight landing site may provide crucial information on the regional context, rock properties and potential structural discontinuities in the crust, which can be utilized to better interpret the future seismic data collected by SEIS (Seismic Experiment for Interior Structures) instrument on board InSight lander. Complementary to the previous studies on the surface properties in the landing ellipses [START_REF] Golombek | [END_REF][4], here we present a study on the regional geologic context of the landing site using a combination of orbital datasets and discuss the interpretation for the lithology and structure in the subsurface of the landing site.

Geologic context: Situated between the Elysium volcanic structure and the dichotomy boundary, the In-Sight landing site region has been shaped by the early formation and subsequent modification of the Martian dichotomy. The landing site is situated in an Early Hesperian transition unit [5][6], inferred to be volcanic in origin [START_REF] Golombek | [END_REF], overlain by the Elysium Mons volcanic unit [7] and younger lavas from Cerberus Fossae [5,7] to the north and the Medusae Fossae Formation to the east [8].

The surface of the landing site is covered by >3 m thick regolith, evidenced by the onset diameter of small rocky ejecta craters [4]. As a consequence, this investigation is focused on the limited exposures in the craters and local scarps in this region.

Method: We analyzed 79 CRISM (Compact Reconnaissance Imaging Spectrometer for Mars) targeted images (18-36 m/pixel, 438 channels) in the landing site region (125°E to 145°E, 2°S to 15°N), following the data processing method in the previous northern lowlands survey [9]. For locations where subsurface mineralogy is detected, we analyzed the accompanying HiRISE (High Resolution Imaging Science Experiment) and CTX (Context Camera) images to investigate the correlation of the mineral detections with morphological units. The annual median thermal inertia of relevant units is derived from measurements by the Thermal Emission Spectrometer [10].

Results: Mafic detections in the volcanic units. Olivine spectral signatures with weak absorptions are detected on small craters rims close to the landing site, on the ridged plains unit (Fig. 1). The mineral detection is well-correlated with the bedrock exposure on the crater wall. The olivine-rich composition; surface expression of wrinkle ridges and the relatively older age compared to nearby volcanic units [5][6] imply their relation to the Hesperian lava flow units widespread in the northern lowlands [11][12][13]. Both olivine and low-calcium pyroxene are found in a small crater within the Elysium Mons unit to the north of the landing ellipse, indicating either a different composition or heavier dust cover masking pyroxene absorptions in the ridged plains unit. More definitive olivine and clay detections are found in the transitional unit close to the dichotomy, which could be the remnants of highlands materials, possibly of similar origin as the adjacent stratigraphy in the Sharp-Knobel watershed [14].

Phyllosilicate detection and layered unit in a 50-km crater. Of particular interest is the Fe/Mg phyllosilicate identified in the central peak of a 50-km crater on the southwest flank of Elysium Mons. The phyllosilicate detection based on the 2.3 and 1.9 µm bands show the strongest absorptions in a topographic depression within the central peak. Similar spectral features with weaker absorptions are found in the adjacent central peak, which consists of a fragmented unit characterized by relatively light-toned features and heavily eroded morphology (Fig. 2C-D). The alternating resistant (cliffforming) and friable (slope-forming) layers in this unit form an irregular cliff-bench morphology, exhibiting extensive layering of <10 m thick. The fragments of 2-5 meters in size within the unit show evidence of little to no relative displacement (Fig. 2E).

Discussion: Possible origins of layered units. The clay-bearing unit in the central peak of a 50-km crater in the Elysium volcanic unit, is found with light-toned, cliff-bench morphology, consistent with the attributes of sedimentary rocks identified elsewhere on Mars [15]. Though the morphology resembles certain locations in the lower member of the Medusae Fossae Formation [16], the unit has distinctively higher thermal inertia (up to 500-1000 J m -2 K -1 s -½ ) than that of the Medusae Fossae Formation.

The feature of incipient fracturing (Fig. 2E) in the layers suggests they likely result from the structural uplift of a homogeneous rock mass at the impact event, instead of deposition at a later stage. The original depth of the central peak unit is >3-4 km, inferred from the size of the impact crater. Such layering units may form via subaerial (aeolian, volcaniclastic) or subaqueous processes before the emplacement of Elysium Mons lavas. The detection of Fe/Mg phyllosilicates associated with the central peak indicates aqueous alteration has occurred, which may be directly related to the presence of water during deposition. Alternatively, the layering structure and fragmentation in the unit could make it more susceptible to alteration after its deposition or excavation by the impact. Though the exposure is limited, it is possible that such unit is more extensive within the subsurface of Elysium Planita. Other large craters in the landing site vicinity also show possible detections of hydrated minerals, but the absorptions are obscured, likely due to a combination of noise and dust mixing.

Inferred stratigraphy and implications for InSight. The overall stratigraphy of the landing site shows a gradational sequence which becomes younger toward Elysium Mons. This analysis suggests a clay-bearing sedimentary unit underlies lava flows of several kilometers thick and a regolith layer close to the InSight landing site. The presence of these different geologic units in the subsurface of the landing site can result in variations in density or seismic anisotropy, which may be analyzed using data collected by the seismic instrument on board InSight lander.

Future work: To better understand the origin of the layered unit, further analysis will focus on its distribution and the thickness and orientation of the layers through continued survey of HiRISE images and digital elevation models. We will also extend the investigation to regions around Elysium Mons and the adjacent highland crust, to gain a wider view of crustal structure and geologic history of the InSight landing site vicinity.

Figure 1 :

 1 Figure 1: Geological map of the InSight landing site vicinity based on [6], with major geological units labeled and inferred age of each unit given in brackets (Am: Amazonian, Hr: Hesperian, N: Noachian). Stars show the mineral detections from CRISM data.

  Planetary Science Conference 2018 (LPI Contrib. No. 2083)

Acknowledgement: This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 751164.


Quang Grégory Rousseau

Huy Tran

Delphine Sinoquet

SCOP: a Sequential Constraint-free Optimal control Problem algorithm

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Optimal control problems involving dynamic systems and state boundary constraints can be solved following various approaches. These approaches can be mainly divided into two families: (i) direct methods, which assume that the optimal control problem has been discretized in time and space, and which lead to classical but often large-sized nonlinear optimization problems; (ii) indirect methods ( [START_REF] Bryson | Applied Optimal control[END_REF]), which take advantage of Pontryagin's maximum principle at the continuous level, using shooting algorithms.

While direct methods are well suited for non-analytic models, their main drawback is the computational time. A good example is the Dynamic Programming algorithm (DP) ( [START_REF] Bryson | Applied Optimal control[END_REF], [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]). Nevertheless DP can be easily extended to take state constraints into account.

Besides, Pontryagin's principle gives optimality conditions which are often solved by a shooting method. However, the presence of state constraints is a tremendous obstacle, and requires the use of a multiple shooting algorithm. Furthermore, the active state constraints should be known ( [START_REF] Bonnans | Well-Posedness of the Shooting Algorithm for State Constrained Optimal Control Problems with a Single Constraint and Control[END_REF]). This contribution presents optimal control problems with a bounded scalar state. We devise a numerical method that does not require any a priori knowledge on active constraints. Proof of convergence is given for two applications: a hybrid vehicle and an elastic rope. Finally, we carry out a numerical test in order to show that the method still works fine on a realistic case: this application is the optimization of torque split between the engine and the electric motor of a hybrid vehicle with respect to fuel consumption with constraints on the state of charge of the battery ([1], [START_REF] Guzzella | Vehicle Propulsion Systems: Introduction to Modeling and Optimization[END_REF], [START_REF] Sciarretta | A Real-Time Optimal Control Strategy for Parallel Hybrid Vehicles with on-board Estimation of the Control Parameters[END_REF]).

II. DESCRIPTION OF THE METHOD

A. Statement of the problem

Let us consider the abstract optimal control problem (P) min

u(•)∈U S 0 L(z(s), u(s), s) ds (1)
subject to

dz ds (s) = f (z(s), u(s), s) (2a) z(0) = ζ 0 (2b) z(S) = ζ S (2c) z(s) ≥ z min (2d) z(s) ≤ z max , (2e)
where z(•) ∈ R, u(•) ∈ U, and s > 0 stand for the state, the control and the evolution variable. We assume that

z min < min{ζ 0 , ζ S } ≤ max{ζ 0 , ζ S } < z max . (3)

B. The SCOP method

Our idea consists in solving problem (P) by a sequence of unconstrained subproblems (Q k ), defined as

(Q k )[s k , Z k ] min u(•)∈U s k 0 L(z k (s), u(s), s) ds (4) subject to dz k ds (s) = f (z k (s), u(s), s) (5a) z k (0) = ζ 0 (5b) z k (s k ) = Z k , ( 5c
)
where (s k , Z k ) are two parameters meant to evolve until s k reaches some "correct" value. The overall algorithm reads Algorithm 1. SCOP data: k = 0, s 0 = S, and

Z 0 = ζ S ; begin Solve (Q 0 )[s 0 , Z 0 ] while ∃s ∈ ]0, s k [ | z k (s) ≤ z min or z k (s) ≥ z max do Define ∆ k (s) = max{z min -z k (s), z k (s) -z max } Determine s k+1 = argmax s∈[0,s k ] ∆ k (s) Compute Z k+1 = Π [zmin,zmax] (z k (s k+1 )) Solve (Q k+1 )[s k+1 , Z k+1 ] end end
The symbol Π [zmin,zmax] (Z) denotes the truncation operator

Π [zmin,zmax] (Z) = z min if Z < z min z max if Z > z max . ( 6
)
Our claim is that if

∀s ∈ [0, s k ], z min ≤ z k (s) ≤ z max (7)
occurs at some intermediate iteration k or at the limit k → ∞, then z k (s) coincides with the optimal trajectory z(s) of problem (P) over the interval [0, s k ].

The benefit of this method lies in the fact that each unconstrained problem (Q k ) can be solved either analytically or by a single shooting algorithm via a constraint-free Pontryagin's principle.

III. PONTRYAGIN'S

) 9c
When there is no constraint on the control u, e.g. for U = R, condition (9c) can be expressed, at least formally, as

0 = H u (z(s), u(s), s, p(s)). ( 10
)

IV. ANALYSIS FOR A HYBRID VEHICLE MODEL

In this section, we set s ≡ t, z(s) ≡ x(t), S = T and Z k ≡ X k . Consider the optimal control problem (H ) min

u(•)∈R T 0 α 2 (t)u 2 (t) dt (11) subject to ẋ(t) = -α(t)(1 -u(t)) (12a) x(0) = χ 0 (12b) x(T ) = χ 0 (12c) x(t) ≤ x max . ( 12d
)
This problem corresponds to a simplified model of torque split for a hybrid vehicle powered by an engine and an electric motor (see section VI for a detailed description of a realistic hybrid vehicle). The integrand α 2 u 2 represents the fuel consumption to be minimized, -α(1u) is the battery law, and x(t) is the state of charge of the battery at time t.

Note that x(0) = x(T ) = χ 0 < x max . We assume that α(t) is positive, differentiable and strictly increasing with respect to t ∈ [0, T ]. For convenience, we introduce

β(t) = 1 t t 0 α(τ ) dτ, for t > 0, (13a)
β(0) = α(0). (13b)
It is straightforward to check that β is an increasing function of t ∈ [0, T ].

Lemma 4.1:

There exists a unique t ∈ ]0, T [ such that

α( t) = β(T ) (14)
Proof: This follows from the intermediate value theorem and from the strictly increasing property of α.

First, we write down the exact solution of (H ).

Proposition 4.1: Let x = χ 0 + t[β(T ) -β( t)] ( 15
)
where t is defined in Lemma 4.1. 1) If x max > x, then the optimal trajectory is given by

x(t) = χ 0 + t[β(T ) -β(t)]. (16)
In this case, the constraint x(t) ≤ x max is never active. 2) If x max ≤ x, then there exists a unique t * ∈ ]0, T ] such that, for t ∈ [0, t * ], the trajectory is given by

x(t) = χ 0 - p * 0 2 t - t 0 α(τ )dτ, (17)
where (t * , p * 0 ) is the unique solution to the system

x max -χ 0 = - p * 0 2 t * - t * 0 α(τ ) dτ, (18a) α(t * ) = - p * 0 2 . ( 18b
)
Proof: Case 1. Applying ( 9)- [START_REF] Bonnans | Commande et optimisation de systèmes dynamiques[END_REF] to the Hamiltonian

H(x, u, t, p) = α 2 u 2 -pα(1 -u), ( 19
)
we end up with

ṗ = 0, 2α 2 u + pα = 0, (20)
from which we infer that

p(t) = p 0 , u(t) = -1 2 p0 α(t) . ( 21
)
Integration with respect to time of the dynamic law

ẋ(t) = -(α + 1 2 p 0 ) (22) yields x(t) = χ 0 - p 0 2 t - t 0 α(τ ) dτ. ( 23
)
Invoking the final condition x(T ) = χ 0 , we obtain (16). Since

ẋ(t) = β(T ) -α(t), ẍ(t) = -α(t) < 0, (24)
we see that the only critical point t ∈ [0, T ] achieves a maximum for x(•), i.e.,

max t∈[0,T ] x(t) = x( t) = x, (25)
where x is defined in (15). Therefore, if x max > x, the optimal trajectory is (16). Case 2. If x max < x, let t * be the contact point where the trajectory hits the boundary. From the continuity of the Lagrange multipliers established in [START_REF] Bonnans | Well-Posedness of the Shooting Algorithm for State Constrained Optimal Control Problems with a Single Constraint and Control[END_REF] by Bonnans and Hermant, we derive the matching conditions

x(t * ) = x max (26a) ẋ(t * ) = 0. ( 26b
)
For t < t * , the optimal trajectory x(•) is of the form (23), in which we write the co-state as p * 0 instead of p 0 for clarity. Plugging (23) into (26), we obtain (18).

Let us eliminate p * 0 from (18a)-(18b), so as to have

t * α(t * ) - t * 0 α(τ ) dτ = x max -x 0 . ( 27
)
From ( 27), existence and uniqueness of t * can be obtained by arguing that the function

ψ(t) = tα(t) - t 0 α(τ ) dτ (28)
is strictly increasing and that

ψ(0) < x max -x 0 < ψ(T ). (29)
This completes the proof.

We now proceed to investigate the behavior of SCOP. It is obvious that if x max > x, then once (Q 0 )[t 0 = T, X 0 = χ 0 ] has been solved, the trajectory x 0 (•) is the optimal one (16) and the algorithm stops. Theorem 4.1: For x max ≤ x, the sequence {t k } k≥0 produced by SCOP is well-defined, strictly decreasing and converges to t * as k → ∞.

This sequence is governed by the implicit recursion

α(t k+1 ) = x max -x 0 t k + β(t k ) (30)
and can be interpreted as a fixed point approximation to the continuous relation

α(t * ) = x max -x 0 t * + β(t * ), (31)
which is a consequence of (27).

Proof: Implicit recursion. Except for X 0 = χ 0 , we have X k = x max for k ≥ 1, as long as the algorithm goes on. According to (23),

x k (t) = χ 0 - p k 2 t - t 0 α(τ ) dτ, ( 32
)
at the k-th step. By construction of Algorithm 1, we have

x k (t k ) = x max (33a) ẋk (t k+1 ) = 0. (33b)
Plugging (32) into (33) leads to

x max -x 0 = - p k 2 t k - t k 0 α(τ ) dτ (34a) α(t k+1 ) = - p k 2 , ( 34b
)
from which we deduce (30) by eliminating p k . Existence and uniqueness. We are going to show that if t k > t * , then t k+1 exists and t k+1 > t * , so that by induction the whole sequence can be generated.

Equation ( 30) is guaranteed to have a unique solution as soon as

α(0) < x max -x 0 t k + β(t k ) < α(t k )
The left inequality is obvious, since α(0) = β(0) β(t k ) (and β is an increasing function). As for the right inequality, it is equivalent to

x max -x 0 < ψ(t k ), (35)
where ψ was introduced in (28). However, we have seen in ( 27) that

x max -x 0 = ψ(t * ). ( 36
)
Because ψ is strictly increasing,

ψ(t * ) < ψ(t k ) ⇔ t * < t k , ( 37
)
and this equivalence ensures existence and uniqueness of

t k+1 ∈ ]0, t k [ for t k > t * .
To prove that t k+1 > t * , we resort to the auxiliary function

Φ k (θ) = t k α(θ) - t k 0 α(τ )dτ. ( 38
)
Since α is strictly increasing, Φ k is strictly increasing too. Thus, we simply have to prove that

Φ k (t * ) < Φ k (t k+1 ). (39)
This amounts to

t k α(t * ) - t k 0 α(τ ) dτ < x max -x 0 (40a) = t * α(t * ) - t * 0 α(τ ) dτ, (40b)
where the last equality is due to (27). Now, (40b) can be re-written under the form

(t k -t * )α(t * ) < t k t * α(τ ) dτ, (41)
which holds true because α is strictly increasing. Convergence. Since {t k } is decreasing and bounded from below by t * , it has a limit t. Passing to the limit in the recursion (30), we get the desired result, namely t = t * .

V. ANALYSIS FOR AN ELASTIC ROPE

In this section, the notations are switched to s ≡ x, z(s) ≡ y(x), S ≡ 1 and Z k ≡ Y k . Consider the control problem

(E ) min u(•)∈R 1 0 [ 1 2 u 2 (x) + gy(x)] dx (42) subject to y (x) = u(x) (43a) y(0) = 0 (43b) y(1) = 0 (43c) y(x) ≥ -h. ( 43d
)
This problem corresponds to an elastic rope fixed at its endpoints and bending under a uniform gravity force g. The integrand 1 2 u 2 + gy represents the potential energy to be minimized in order to find the equilibrium vertical position y(x) as a function of abscissa x, under the constraint due to the level -h < 0 of the floor on which the rope can lie.

Proposition 5.1:

The exact solution of (E ) is given by the following rule.

1) If h > 1 8 g, then the equilibrium position is

y(x) = - 1 2 g x(1 -x). ( 44
)
In this case, the constraint y(x) ≥ -h is never active.

2) If h ≤ 1 8 g, let x * = 2h g . ( 45
)
Then, for x ∈ [0, x * ], the equilibrium position is

y(x) = -2hg x + 1 2 gx 2 . ( 46
)
Proof: The proof follows the same steps as in the previous section, first by considering the unconstrained case, then by expressing the matching condition at the contact point. We refer readers to [START_REF] Bonnans | Stability and Sensitivity Analysis for Optimal Control Problems with a First-order State Constraint[END_REF] for further details.

We now address the question of convergence for SCOP. It is obvious that if h > 1 8 g, then once (Q 0 )[x 0 = 1, Y 0 = 0] has been solved, the trajectory y 0 (•) is the optimal one (44) and the algorithm stops.

Theorem 5.1: For h ≤ 1 8 g, the sequence {x k } k≥0 produced by SCOP is well-defined, strictly decreasing and converges to x * as k → ∞.

This sequence is governed by the explicit recursion

x k+1 = 1 2 x k + 2h gx k (47)
and can be interpreted as a fixed point approximation to the continuous relation

x * = 1 2 x * + 2h gx * (48)
that is a consequence of (45).

Proof: Except for Y 0 = 0, we have Y k = -h for k ≥ 1, as long as the algorithm goes on. As a consequence of Pontryagin's principle ( 9)-( 10), the optimal trajectory is of the form

y k (x) = -p k x + 1 2 gx 2 (49)
at the k-th step. By construction of Algorithm 1, we have

y k (x k ) = -h (50a) y k (x k+1 ) = 0. ( 50b
)
Plugging ( 49) into (50) leads to

-h = -p k x k + 1 2 gx 2 k (51a) 0 = -p k + gx k+1 , (51b)
from which we deduce (47) by eliminating p k . It is now an easy algebra exercise to check that the sequence {x k } k≥0 is decreasing and converges to x * .

It is remarkable that the sequence (47) is the classical Babylonian algorithm that computes the square-root of the positive real number 2h/g. Fig. 1 and 2 display the SCOP results compared to those obtained with a multiple shooting algorithm which needs to know the active constraints. The two trajectories superimpose perfectly. Fig. 2 illustrates that SCOP retrieves the exact position of x * .

VI. NUMERICAL RESULTS FOR A REALISTIC HYBRID

VEHICLE

In this section, SCOP algorithm is applied to a realistic model of a hybrid vehicle. The engine is characterized by its static fuel consumption map, depending on engine speed ω and delivered torque T eng . For this application, a quadratic polynomial is used The electric motor is characterized by its torque T mot , and is fed by a battery. The state variable x, the state of charge of the battery, is defined by

L(ω, T eng ) = 2 i,j=0 K ij ω i T j eng .
ẋ(t) = - ω(t)T mot (t)K U batt (t)n capa = -Kω(t)T mot (t), (52)
with U batt the battery voltage, assumed to be constant, ω(t) the electric motor speed, K a scaling constant, and n capa the nominal capacity of the battery. Then the state x depends on the electric power.

In this problem, a prescribed vehicle cycle is imposed (the Artemis Urban cycle for instance, see [START_REF] André | The ARTEMIS European driving cycles for measuring car pollutant emissions[END_REF] for more details), and supplies the speed ω(t) and torque T rq (t) of the vehicle. We assume that both the engine and the electric motor have the same rotation speed ω, and that the requested torque T rq must be provided according to T rq = T eng + T mot .

Let us introduce the torque split control u(t), defined as

u(t)T rq (t) = T eng (t) (53a) (1 -u(t))T rq (t) = T mot (t). (53b)
Because of maximum and minimum electric motor torque (during battery regeneration), and of maximum engine torque, the control u(t) is constrained by values u min and u max depending on the speed ω(t).

The resulting optimal control problem is thus similar to problem (1-2) with

L(u(t), t) = 2 i,j=0 K ij ω i (t)T j rq (t)u j (t) (54)
and

f (u(t), t) = -Kω(t)T rq (t)(1 -u(t)). ( 55
)
We apply SCOP algorithm with the following expression of u * k , solution of (Q k ), deduced from ( 9):

u * k (t) = - 2 i=0 K i1 ω i (t) + p k Kω(t) 2 2 i=0 K i2 ω i (t)T rq (t) . ( 56
)
Fig. 3 shows for prescribed driving cycle (Artemis Urban cycle) (a) the requested speed of the vehicle, (b) the corresponding engine speed (depending on the vehicle characteristics), (c) the requested torque, that can be provided by the engine, or by the electric motor, or both. Fig. 4 shows a comparison between the state trajectory obtained with SCOP, and those obtained with a classical dynamic programming algorithm.

The trajectories are very close to one another. However, the computational times are not the same: the dynamic programming algorithm needs 135s for dx = 0.5, while SCOP needs only about 1s to obtain the optimal trajectory. One can also notice that the dynamic programming trajectory converges to SCOP trajectory, as the state step size decreases and tends to 0.

An interesting result is that the analytic expression of the optimal trajectory (56) can be used to build a realtime control strategy parameterized by p (see [START_REF] Rousseau | Constrained Optimization of Energy Management for a Mild-Hybrid Vehicle[END_REF], [START_REF] Guzzella | Vehicle Propulsion Systems: Introduction to Modeling and Optimization[END_REF], [START_REF] Sciarretta | A Real-Time Optimal Control Strategy for Parallel Hybrid Vehicles with on-board Estimation of the Control Parameters[END_REF]). For instance, an efficient control is ECMS (Equivalent fuel Consumption Minimization Strategy), where a parameter needs to be defined to determine the equivalence between fuel consumption and battery energy. This parameter is proportional to p.

VII. CONCLUSIONS AND FUTURE WORKS

In this paper, we have proposed an original method to solve optimal control problems with state constraints. A proof of convergence, in some particular cases, has been established. Although the proof of convergence needs some strong assumptions, SCOP has been successfully applied to more general cases.

The first application is the energy minimization of a rope fixed at its endpoints. This problem can also be solved by using a multiple shooting algorithm, however for this algorithm the shape of the trajectory (with number of junction points) needs to be known.

The second presented application is the optimization of fuel consumption of a hybrid vehicle. In this problem, it is often impossible to know the number of active constraints and when a state constraint becomes active. Classically a dynamic programming (DP) algorithm is used to solve this kind of problem. SCOP allows to retrieve the optimal trajectory of state-constrained optimal problems in a very small computational time, about 100 times faster than DP. Thus the presented applications of SCOP illustrate its potentiality compared to classical methods. However, the range of applicability of the method remains to be clarified.

Fig. 1 .

 1 Fig. 1. Elastic line trajectory constrained by the floor level

Fig. 2 .

 2 Fig. 2. Elastic line trajectory constrained by the floor level (zoom)

Fig. 3 .

 3 Fig. 3. Vehicle speed requested by Artemis Urban cycle (top) / requested engine speed (middle) / requested torque (bottom)

Fig. 4 .

 4 Fig. 4. Optimal state trajectories obtained with DP and SCOP (empty battery : x min = 0, full battery : xmax = 100).

VIII. ACKNOWLEDGMENTS

The authors wish to thank Pierre Rouchon and Audrey Hermant for several stimulating discussions.


Alexis Dehais Underdown
email: vsq.alexis.dehais-underdown@sorbonne-nouvelle.fr

Lise Crevier Buchman

Didier Demolin

ACOUSTICO-PHYSIOLOGICAL COORDINATION IN THE HUMAN BEATBOX: A PILOT STUDY ON THE BEATBOXED CLASSIC KICK DRUM

Keywords: Human Beatbox, phonetics, aerodynamics, laryngoscopy, acoustics

This paper presents a pilot study of physiological, acoustic and aerodynamic characteristics of the beatboxed classic kick drum by a single artist excerpted from our database of his beatboxing repertoire. We recorded 25 repetitions with synchronized aerodynamic (intraoral pressure, oral airflow) acoustic, electroglottographic and laryngoscopic data. Results show that the classic kick drum is produced either as a bilabial voiceless glottalic egressive plosive or as a bilabial voiceless glottalic egressive affricate. Laryngoscopic data show a glottal adduction, a slight supraglottal constriction and a laryngeal raising. We also observed an occlusion, a high amplitude burst and low frequencies in the spectrum. We discuss the coordination of the gestures involved in the classic kick production based on the multi-instrumental recording. Future analysis of the beatboxed database will help us understand beatboxing physiology and the extent of the vocal tract capacities.

INTRODUCTION

The Human Beatbox (HBB) make use of the vocal tract in order to imitate musical sonorities. Few studies mention physiological aspects of this peculiar technique. MRI studies of beatboxing [START_REF] Blaylock | Sounds of the Human Vocal Tract[END_REF][START_REF] Patil | Comparison of Basic Beatboxing Articulations between Expert and Novice Artists using Real-Time Magnetic Resonance Imaging[END_REF][START_REF] Proctor | Paralinguistic mechanisms of production in human "beatboxing": A real-time magnetic resonance imaging study[END_REF] give a description of basic and complex articulations of HBB production. Laryngoscopic studies of laryngeal behaviour during beatboxing [START_REF] De Torcy | A videofibroscopy study of laryngo-pharyngeal behaviour in the human beatbox[END_REF][START_REF] Sapthavee | Functional endoscopic analysis of beatbox performers[END_REF] show that beatboxers manage to recruit laryngopharyngeal articulators (i.e. pharynx constrictor, glottis, ventricular folds, aryepiglottic folds) separately and optimally. Understanding this technique would provide valuable knowledge about the extent of the vocal tract capacities' and its limits aside from its use for linguistic purpose.

This paper aims to describe physiological aspects of the beatboxed classic kick drum by a French beatboxer using multi-instrumental simultaneous recordings with aerodynamic, acoustic and laryngoscopic data. To our knowledge there is no existing study of the beatboxing vocal technique using multidimensional synchronous data. This multiparametric study was found useful to understand the physiology of the beatboxed classic kick drum, in particular the coordination of supralaryngeal and laryngeal gestures. The beatboxed classic kick drum is described [START_REF] Blaylock | Sounds of the Human Vocal Tract[END_REF][START_REF] De Torcy | A videofibroscopy study of laryngo-pharyngeal behaviour in the human beatbox[END_REF][START_REF] Patil | Comparison of Basic Beatboxing Articulations between Expert and Novice Artists using Real-Time Magnetic Resonance Imaging[END_REF][START_REF] Proctor | Paralinguistic mechanisms of production in human "beatboxing": A real-time magnetic resonance imaging study[END_REF][START_REF] Sapthavee | Functional endoscopic analysis of beatbox performers[END_REF] as a voiceless glottalic egressive bilabial produced as a plosive [p'] or an affricate [p'ɸ].

METHODS

The subject was a 35 y.o. male who volunteered for the experiments after signing an informed consent in accordance with the ethical committee of our hospital (N° 1922081). The subject was asked to perform his repertoire of beatboxed sounds. First he gave a description of the sounds, then produced them in isolation and finally, he combined the sounds in a beatboxed sequence of various sounds from his repertoire. For analytical purposes and to understand how the synchronisation of simultaneous parameters (articulation, aerodynamics, acoustics) is involved, we chose to focus only on the classic kick drum production.

Aerodynamic and electroglottographic data

Aerodynamic data was collected with the EVA2 workstation (SQLab-LPL, Aix en Provence, France) which allows synchronized aerodynamic, acoustic and electroglottographic (EGG) data collection. Intraoral pressure (Po), expressed in hectopascal (1 hPa = 1,02 cm H2O), was obtained inserting a small tube into the mouth; we took care to place the tube at the corner of the lip so it would not interfere with Oral airflow. Oral Airflow (Oaf), expressed in cubic decimeter per second (dm 3 /s), was collected using a flexible silicone mask pressed on the subject's mouth.

The EGG signal (Glottal Enterprise) was obtained by placing two electrodes on both sides of the thyroid cartilage. The EVA2 station is equipped with a microphone allowing to retrieve the acoustic signal. We analyzed manually the data with PHONEDIT Signaix [START_REF]PHONEDIT Signaix (Version 4.2.0.8[END_REF] to get (1) the duration based on the audio waveform, (2) the peak of Po for each token, (3) the peak of Oaf and the air volume in cm 3 for each token. We analyzed the global shape of the EGG signal to confirm that the beatboxed classic kick drum uses a glottalic egressive airstream.

Acoustic data

We used Praat [START_REF] Boersma | Praat: doing phonetics by computer[END_REF] to run an acoustic analysis based on the audio signal together with the EGG signal that were both recorded during the aerodynamic experiment. First we segmented the audio on a textGrid and for each token we labelled the acoustic phase of the classic kick drum: the occlusion, the burst (i.e. impulsional noise produced when the lips come apart), the release noise (i.e. noise generated by the air passing through the lips after the lip aperture) and the frication noise (i.e. the noise generated during a fricative consonant, only for [p'ɸ]). Finally, for each phase we were able to extract the duration in milliseconds (ms), the intensity in decibels (dB) and the FFT spectrum (25ms window) in order to get the centre of gravity (CoG in Hz), the skewness (i.e. spectral asymmetry distribution) and kurtosis (i.e. peakedness of the distribution). Since a silicone mask was used during the experiment, we were able to keep the same distance between the participant mouth and the microphone, thus allowing us to analyze the intensity values.

Laryngoscopic data

The laryngeal examination was performed by the second author using a nasal flexible fiberscope (Kay-Pentax®FNL10RP3) with a DigitalStrobe®, RLS91000 (Kay Elemetrics, Lincoln Park, NJ, USA) station. The video rate was 25 frames per second (fps). Laryngoscopy is an invasive procedure that allows visualizing laryngeal structures without interfering with phonation or articulation. The subject did not ask for any anesthetic. We used VirtualDub software to visualize the video (13,72s duration) and to extract each image: 13,72s x 25 fps (frames per second) = 343 images. We focused on describing glottal adduction and supraglottal contraction.

RESULTS

We found two types of classic kick drum, one is a bilabial voiceless glottalic egressive [p']; the second was produced as a bilabial voiceless affricate [p'ɸ].

The second sound appears only in a beatboxed sequence of kick-drums, hi-hats (i.e. [t ͡ s']) and snare drum (i.e. [k']); for now, we do not know if it is a stylistic variant or an effect of respiratory control.

Aerodynamic and EGG data

We analyzed 24 productions of beatboxed classic kick drum, 18 were unaffricated; 6 were affricated. On Fig. 1, Po increases up to 20 hPa indicating that an occlusion occurs. The EGG signal shows that the glottis is closed until the end of [p']. Shortly after the occlusion occurs, we note that Oaf is positive, there is an egressive airstream. Taken together, the data prove that the beatboxer is using a glottalic egressive airstream to produce the classic kick drum effect. As shown in Fig. 1, the affricated classic kick drum adds a pulmonic fricative [ɸ]. After the release of this type of kick drum, Po starts dropping but as soon as the glottis opens on the EGG signal the Po drops more slowly and Oaf exhibits another peak. This second peak of Oaf indicates a second constriction is created in the vocal tract. It results in an increased velocity of the volume passing through the constriction point, made by the lips. The coordination of the occlusion and release of the beatboxed classic kick drum is very peculiar. Indeed, while Po still builds up, Oaf shows that air is coming out of the mouth. There seem to be a leakage during the lip occlusion. We hypothesize that while maintaining the occlusion, the beatboxer's lips become slightly opened by a little slit, thus allowing the pressure to keep raising and the air to come out. This type of coordination (i.e. closure and leakage) is normally not found in normal and casual speech. This appears to be a specificity of beatboxing or at least a specificity of the beatboxer's production since all tokens showed this type of coordination.

Table 1 shows the mean and standard deviation of the duration (ms), based on the acoustic waveform, Po (hPa), Oaf (dm 3 /s) and the air volume for [p'] and [p'ɸ]. Although the Po saturated at 20 hPa, it is consistent with an ejective production [START_REF] Demolin | Acoustic and aerodynamic characteristics of ejectives in Amharic[END_REF]. Concerning Oral Airflow, the values are stable and similar for both [p'] and [ɸ]. The volume of air for the fricative is higher and more variable than for the ejective. The use of different airstream mechanisms, one glottalic egressive and the other pulmonic, explains the discrepancy between volume values. Indeed, the pulmonic fricative uses a larger air volume in the lungs whereas the ejective uses only the volume contained between the glottis and the lips. We hypothesize that this alternation between pulmonic and glottalic airstreams may help to keep a fluid rhythm while beatboxing.

) 48( 11) 43( 6) 102(38)

Acoustic data

We extracted the audio data collected during the aerodynamic experiments and analyzed (1) the duration of the occlusion, burst and frication, (2) the intensity, (3) the resonance frequency. Concerning the duration of the occlusion (Table 2), due to the particular coordination of occlusion/release, not all of the classic kick drum showed an occlusion preceding the release. Indeed, only 18 productions (16 unaffricated and 2 affricated classic kick drums) showed evidence of closure. The mean duration and standard deviation of the closure is 7(1) ms which indicates a quick occlusion and little variation. The mean duration and standard deviation of the burst is 14(15)ms which indicates there is variation across the duration of the release. Finally, the duration of the frication of [p'ɸ] is longer, the mean and standard deviation are 137(49)ms. This difference of duration may be explained by the fact that to sustain a continuous frication, articulatory adjustments are necessary. Indeed, as pointed out by Signorello et al. [START_REF] Signorello | Toward an aerodynamic model of fricative consonants[END_REF], the size of the constriction is to be controlled and both intraoral differential pressure (i.e. ΔPintraoral = Pintraoral -Patmospheric) and subglottic differential pressure (ΔPsubglottic = Psubglottic -Pintraoral) must reach a specific threshold to produce and sustain audible frication.

Regarding intensity values, the burst shows the highest values with little variation across tokens. It can be explained by the fact that the classic kick drum uses a glottalic egressive airstream in order to increase Po resulting in a strong intensity burst. The noise release intensity is lower and indicates that intensity is dropping after the lips come apart. Finally, the frication noise is the less intense but the more variable. Signorello et al. [START_REF] Signorello | Toward an aerodynamic model of fricative consonants[END_REF] reported values between 52 dB and 64 dB for [f]; they used the same experimental methods. We hypothesize that the frication noise is in fact a bilabial fricative rather than a labiodental one. Indeed, the turbulence may be created by soft tissues (i.e. lips) and not by hard tissues (i.e. teeth).

Mean (standard deviation) in the spectral analysis revealed low frequencies and a peaked spectrum. The CoG is 209(35) Hz and 142(72) Hz at the offset. The burst and noise of release show similar values concerning the CoG and the skewness. The classic kick seems to be characterized by low frequency components. Concerning the kurtosis values we values are different, the burst is a bit flatter than the noise. Conversely, the frication noise shows a different pattern. The fricative component has a higher CoG, a lower skewness and a very low kurtosis when compared to the burst and the noise. It indicates a flat low frequency spectrum.

Laryngoscopic data

The laryngoscopic experiment took place during a second session on the same day. We analyzed 25 repetitions (22 unaffricated and 3 affricated) of the beatboxed classic kick drum in order to describe the laryngeal behavior. Fig. 2 presents an unaffricated classic kick drum together with the acoustic signal extracted from the video. During the initial position (Fig. 2 frame 1) the vocal folds are abducted and the arytenoid cartilages are apart. The pharyngeal cavity is 'open', indeed the hypopharynx shows no sign of contraction. Then, a full vocal fold adduction and an incomplete adduction of the ventricular folds is achieved (Fig. 2 frame 1 to 3). We noted an incomplete anteroposterior closure (i.e. arytenoids are getting closer to the epiglottis) of the arytenoid cartilages and the aryepiglottic folds (Fig. 2 frame 4 to 6). This diminishes both the vocal folds' and the epilaryngeal tube's length. Light contrast (i.e. increased luminosity) on frame 7, 8, 9 may suggest that the larynx has raised. The acoustic burst occurs at the end of frame 6 and the acoustic release corresponds to frame 7 where the glottis opens. During the unaffricated classic kick drum (figure 2) the vocal folds remain adducted while the larynx is raising, whereas during the affricated classic kick drum (not shown) the larynx raising starts with an adducted glottis that opens before the larynx reaches its highest position. This difference may be a sign of coordination between a glottalic egressive airstream and a pulmonic egressive airstream. The raising ends with a wide opened glottis compared to the neutral position. This wide aperture results in a larger volume of air coming from the lung to produce the turbulent airflow of the bilabial fricative component.

DISCUSSION

With a multiparametric methodology we were able to describe the physiology the beatboxed classic kick drum produced by a single beatboxer. Our findings concur with other studies about beatboxing. Concerning the laryngeal behavior, [START_REF] De Torcy | A videofibroscopy study of laryngo-pharyngeal behaviour in the human beatbox[END_REF][START_REF] Sapthavee | Functional endoscopic analysis of beatbox performers[END_REF] also found a glottal adduction and the absence of major supraglottal compression. We found a retraction of the epiglottis, possibly due to a tongue posteriorization. This might be a specific gesture of this present study's subject. This laryngeal behavior is not specific to beatboxing since the same gesture is com-monly used to achieve glottal stops [START_REF] Edmondson | The valves of the throat and their functioning in tone, vocal register and stress: laryngoscopic case studies[END_REF][START_REF] Esling | Glottal stop, glottalized resonants, and pharyngeals: A reinterpretation with evidence from a laryngoscopic study of Nuuchahnulth (Nootka)[END_REF]. Beatboxing studies [START_REF] Blaylock | Sounds of the Human Vocal Tract[END_REF][START_REF] De Torcy | A videofibroscopy study of laryngo-pharyngeal behaviour in the human beatbox[END_REF][START_REF] Patil | Comparison of Basic Beatboxing Articulations between Expert and Novice Artists using Real-Time Magnetic Resonance Imaging[END_REF][START_REF] Proctor | Paralinguistic mechanisms of production in human "beatboxing": A real-time magnetic resonance imaging study[END_REF][START_REF] Sapthavee | Functional endoscopic analysis of beatbox performers[END_REF] concluded that the classic kick drum is produced with a glottalic egressive airstream, our data suggest the same conclusion.

Aerodynamic data brought light on the coordination of gestures used to achieve the classic kick drum. Normally the production of a plosive is defined by a closure in the vocal tract increasing the pressure behind the constriction point and a decrease in Po after the constrictions' release. In our case the bilabial closure is simultaneous with a leakage suggesting an atypical coordination of common gestures (i.e. closure and aperture the Vocal Tract).

The discussion about the coordination of the gestures implicated in the classic kick production brings us to the heart of the problem: what are the primitives of the system? Is Beatboxing relying on gestures or features? On which mental representations does the Beatboxer rely on? All those issues will be resolved by future studies on this technique with multi-instrumental experiments. Indeed, our multiparametric study was useful because we could not only come to the same general conclusion as in [START_REF] Blaylock | Sounds of the Human Vocal Tract[END_REF][START_REF] De Torcy | A videofibroscopy study of laryngo-pharyngeal behaviour in the human beatbox[END_REF][START_REF] Patil | Comparison of Basic Beatboxing Articulations between Expert and Novice Artists using Real-Time Magnetic Resonance Imaging[END_REF][START_REF] Proctor | Paralinguistic mechanisms of production in human "beatboxing": A real-time magnetic resonance imaging study[END_REF][START_REF] Sapthavee | Functional endoscopic analysis of beatbox performers[END_REF] about the classic kick drum (i.e. glottalic egressive airstream) but also provide finer phonetic details such as the gestures involved and their coordination. In future studies we will use the same methods in order to gather more data from more subjects and assess the reproducibility of the present results. We hope to be able to answer the previous questions about the primitives of beatboxing.

CONCLUSION

This pilot study based on aerodynamics, EGG, acoustics and laryngoscopy showed the valuable information extracted from the data to discuss fine phonetic details in beatboxed sounds and to understand the gestures and their coordination.

Figure 1 :

 1 Figure 1: Acoustic waveform, EGG signal, Po (hPa) and Oaf (dm 3 /s) of an affricated classic kick [p'ɸ]. The red dotted line indicates the acoustic release.

Figure 2 :

 2 Figure 2: Acoustic waveform and its synchronized laryngoscopic images of an unaffricated classic kick-drum. Frame 1 to 3 show glottal adduction involving vocal folds and ventricular folds, frame 4 to 6 show supraglottal compression involving arytenoid cartilages and the aryepiglottic folds, frame 7 to 9 indicate glottal abduction involving vocal folds and ventricular folds (Frame 1: Vf = Vocal folds, VeF = Ventricular Folds, Ary = Arytenoids, Ary_F = aeryepiglottic folds, Epi = Epiglottis).

Table 1 :

 1 Mean (standard deviation) duration (ms), Po (hPa), Oaf (dm 3 /s) and air volume (cm 3 ) for [p'] and [p' ɸ].

		unaffricated	affricated
		[p']	[p']	[ɸ]
	Duration (ms)	65(9)	50(7)	133(40)
	Po (hPa)	19,7(0,7)	19,7(0,4) 4,1(1,4)
	Oaf (dm 3 /s)	1,1(

Table 2 :

 2 Mean

		(standard deviation) of duration
	(ms), intensity (dB), centre of gravity (Hz),
	Skewness and Kurtosis.
		Closure	Burst	Noise	Frication
		(n=18)	(n=25)	(n=25)	(n=8)
	Duration	6(1)	14(15)	41(12)	137(49)
	Intensity		68(4)	60(5)	49(6)
	CoG		243(44) 176(91)	531(268)
	Skewnes		9(1,8)	11(4,6)	4,5(2)
	s
	Kurtosis		133(41) 239(121)	39(32)

AKNOWLEDGEMENT

We thank the French beatboxer JOOS for participating to this experiment.


Hoël Langouët

Ludovic Métivier

Delphine Sinoquet
email: delphine.sinoquet@ifpen.fr

Quang-Huy Tran

Q.-H Tran

Engine calibration: multi-objective constrained optimization of engine maps

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

We present two new approaches to address the optimization problem associated with engine calibration. In this area, the tuning parameters are traditionally determined in a local way, i.e., at each engine operating point, via a single-objective minimization problem. To overcome these restrictions, the first method we propose is able to cope with several objective functions simultaneously in the local formulation. The second method we put forward relies on a global formulation, which allows the whole driving cycle to be taken into account while remaining single-objective. At the practical level, the two methods are implemented by combining various existing techniques such as the LoLiMoT (Local Linear Model Tree) parameterization and the MO-CMA-ES (Multi-Objective Covariance Matrix Adaptation Evolution Strategy) algorithm. A better compromise appears to be achieved on real case applications.

 (r, c)normalized Gaussian function N λ (r, c) elementary shape function

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, which allows us 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, methoden und anwendungen in der Motorenentwicklung[END_REF][START_REF] Hafner | 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 based on design of experiments[END_REF][START_REF] Sampson | Recent developments in techniques for model-based calibration[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.

The classical approach is based on a local statistical modeling of the engine responses at a limited set of operating points (given engine speed and engine torque or load) from experimental measurements at engine test bench. Then, for each operating point, a local optimization of these models is performed. This optimization problem is often a bi-objective problem which is transformed into a mono-objective problem with a weighted sum of the objectives. The choice of the weights being far from obvious, several trials are often necessary to find an acceptable solution [START_REF] Das | Normal-boundary intersection: a new method for generating the Pareto surface in nonlinear multicriteria optimization problems[END_REF][START_REF] Leyffer | A complementarity constraint formulation of convex multiobjective optimization problems[END_REF]. In this paper, we propose to adapt the Multi-Objective Covariance-Adaptation Evolution Strategy method (MO-CMA-ES) proposed by [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation[END_REF] to the problems arising in this application: linear and nonlinear constraints are introduced via a suitable penalty procedure.

As pointed out by many previous works [START_REF] Castagné | Comparison of engine calibration methods based on design of experiments[END_REF]Buzy et al. 2005;[START_REF] Schmied | A global constrained optimization algorithm for engine calibration[END_REF]), the local approach leads to several difficulties: the definition of the objectives of the local optimization at each operating point is cumbersome, the synthesis of the resulting local optima in engine maps (maps that are implemented in the engine control unit) leads to a difficult smoothing phase that usually ends up with a deterioration of the optimized operating points. Moreover, this method is time-consuming due to the modeling and the optimization phases for each operating point [START_REF] Castagné | Comparison of engine calibration methods based on design of experiments[END_REF]. In this respect, [START_REF] Schmied | A global constrained optimization algorithm for engine calibration[END_REF] proposed a global optimization approach which directly optimizes the engine maps on the whole driving cycle instead of the individual optimization of each engine operating point. In one single step, his method optimizes all the chosen operating points and directly delivers the engine maps that are implemented in the unit control of the vehicle. Nevertheless, the implementation of this approach is difficult due to the large number of unknowns of the optimization problem and the number of smoothing constraints to ensure feasibility and drivability [START_REF] Neßler | Global evaluation of the drivability of calibrated diesel engine maps[END_REF][START_REF] Schloßer | Model-based calibration of SI-engines map optimization[END_REF]). To circumvent these difficulties, we advocate the use of an adapted parameterization of the maps based on LoLiMoT models [START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF][START_REF] Hafner | The use of stationary and dynamic emission models for an improved engine performance in legal test cycles[END_REF], which allows a priori information on the engine maps to be incorporated naturally and the number of parameters to be reduced significantly.

Section 2 recalls the classical steps of the calibration process, then discusses the difficulties that have motivated our study. In Sect. 3, we outline the Multi-Objective Covariance-Adaptation Evolution Strategy method. Finally, in Sect. 4, the map optimization of the whole driving cycle is described. The potential of the two methods is illustrated by an application on a real calibration case for a turbo-charged commonrail diesel engine.

Engine calibration

Sketch of the engine calibration process

The emission calibration workflow is classically divided into three steps [START_REF] Röpke | DoE-design of experiments, methoden und anwendungen in der Motorenentwicklung[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.

(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 Fig. 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 (New European Driving Cycle) simulation and a selection of 17 OPs. The cumulated levels of pollutant emissions along the cycle are computed as weighted sums of the pollutant emissions for each chosen OP.

Phase 2 consists of five steps:

(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 [START_REF] Castagné | Comparison of engine calibration methods based 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 [START_REF] Röpke | DoE-design of experiments, methoden und anwendungen in der Motorenentwicklung[END_REF]). (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. 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 experiments. 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 [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 based on design of experiments[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 Fig. 7 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 Sect. 4.

The OP optimization problem formulation

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 Fig. 2). It leads to the multi-objective optimization problem

min p∈R N {f 1 (p), f 2 (p), . . . , f n (p)} (1a) subject to ≤ Ap ≤ u, (1b) f j (p) ≤ s j , j ∈ [1, d], (1c)
and d ≤ n. A classical example in turbo-diesel engine is the two-objective problem in which one attempts to minimize the particulate emissions and the NO x emissions under constraints on fuel consumption (CO 2 emissions), on CO and HC emissions and engine noise level.

Fig. 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

Traditionally, the multi-objective problem (1) is replaced by a constrained optimization problem where only one of the objectives is actually minimized, the other being constrained to remain below the prescribed thresholds. In other words, if the top priority is given to the i 0 -th objective function, then what we perform is min

p∈R N f i 0 (p) (2a) subject to ≤ Ap ≤ u, ( 2b
)
f j (p) ≤ s j , j ∈ [1, d]. ( 2c
)
This constrained single-objective 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 (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] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation[END_REF] was 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

In the CMA-ES evolutionary algorithm, a new population x (g+1) is created from the current population x (g) by generating realizations of a Gaussian probability law

x (g+1) k ∼ N (e (g) , [σ (g) ] 2 C (g) ).
(3)

Each subscript k ∈ [1, K] corresponds to an individual of the population. In the distribution law N , the expected value e (g) is chosen to be

e (g) = μ k=1 w k x (g) k:K , ( 4a
)
where the notation x (g) k:K stands for the k-th best individual of the current population x (g) with respect to the function to be minimized, and

μ k=1 w k = 1 and 0 ≤ w μ ≤ w μ-1 ≤ • • • ≤ w 1 ≤ 1. ( 4b
)
In doing so, we ensure that the next population is 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 μ k=1 w k x (g+1) k:K -e (g) σ (g) x (g+1) k:K -e (g) σ (g) T , (5)
which amounts to adding to the current covariance matrix a μ-rank term from the best μ individuals of the new population. An alternative update formula is proposed by Hansen, which takes into account not only the last iteration as in ( 5) but also all the previous iterations from the beginning. In a similar fashion, 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 e (g) of the population 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 toward the global minimum.

MO-CMA-ES method for multi-objective optimization

Adapting the CMA-ES method to multi-objective optimization requires the definition of new criteria to define what is a good and a bad individual in the 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:

1. The first criterion is the Pareto dominance. Let us recall this notion.

Definition 1 An individual x is said to dominate another individual x , which is symbolized by x ≺ x , if and only if

• ∀i ∈ {1, . . . , n}, f i (x) ≤ f i (x ); (6a) • ∃i ∈ {1, . . . , n}, f i (x) < f i (x ). ( 6b
)
Definition 2 The non-dominated subset ndom(X) of a given set X is defined as

ndom(X) = x ∈ X | x ∈ X, x ≺ x . ( 7
)
The elements of this subset are said to be Pareto-optimal points.

Additionally, the Pareto front is composed of the values taken by the objectives for all the Pareto optimal solutions. Examples of Pareto fronts are given in Fig. 6. The individuals of a given population X are sorted by their non-dominance level. The notion of level is defined by induction as follows. First, we consider that the elements of ndom(X) have a non-dominance level L = 1. Then, we introduce the points of level L = 2 as the non-dominated points of the population X without the non-dominated solutions. For L ≥ 1, we thus obtain the recursive definition

dom 0 (X) = X; (8a) ndom L (X) = ndom(dom L-1 (X)); (8b) dom L (X) = dom L-1 (X) \ ndom L (X). ( 8c
)
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 Fig. 4 (one rectangle for each non-dominated point).

Combining these two criteria allows to be sorted 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]. In this section, we present some 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. Figure 5 illustrates the quality of the models at some validation points chosen randomly in the admissible domain (those points were not used for model building). An average value of the measurement error has been estimated from some repeated measurements. The model of NO x emissions matches very well the data within the error bars whereas the model of particulate emissions is less accurate: this is due to large error measurements for the engine response and a bad approximation of the measurement error when considering this error to be constant. Nevertheless, engine experts have validated this model considering that the trends were well modeled. 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 flow and the boost pressure.

The optimization objectives are the particulate emissions and the NO x emissions, CO 2 , 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 nonlinear inequality constraints: they have been introduced via a l 1 penalty term added to the objectives (if the j -th constraint is not satisfied, the misfit with the bound is added to the objective functions, the penalty weight varying during the iterations).

The results are presented in Fig. 6: two optimizations are compared, one without taking into account the parameter dispersions and one with the robustness constraint Fig. 5 Comparison of experimental data and model predictions for given parameter sets (those points were not used for model building) at OP 1750 rd/mn-1 bar for NO x (a) and particulate emissions (b) in g/h. Error bars correspond to measurement error estimates from some repeated experiments 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 Fig. 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 Fig. 6 by an important shift of the Pareto front and a modification of its shape. The points indicated in pink 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.

Finding a NO x /particulate emissions compromise is a classical calibration problem for diesel engine. Engine experts were thus very interested in the different compromises between the two objectives that have been obtained by the multi-objective optimization method. We clearly observe 3 interesting valleys on the Pareto front with a significant decrease of NO x emissions for a small increase of particulate emissions. This method allows then a finer choice of the calibration control parameter set than a classical method based on a mono-objective formulation.

Solving the map optimization problem

In the classical calibration process described in Sect. 2, 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 smooth-ing step. The optimization problem is formulated as

min m:R 2 →R M T 0 F i (r(t), c(t), m(r(t), c(t))) dt 1≤i≤n (9a) subject to (r, c) ≤ Am(r, c) ≤ u(r, c), (r, c) ∈ R 2 , ( 9b
) T 0 F j (r(t), c(t), m(r(t), c(t))) dt ≤ S j , j ∈ [1, d]. (9c)
Here,

• t ∈ R → (r(t), c(t)) ∈ R 2 denote the (known) trajectory of the cycle within the engine speed-load domain; r,c) ∈ R M are the 2D engine maps; the components of m are called engine control parameters; • F i is the model of the engine response i, as a function of the engine control parameters but also of the speed and load.

• (r, c) ∈ R 2 → m(
The objectives to be minimized (or constrained) are the engine responses cumulated on the considered driving cycle: the cumulated responses are here (see ( 10)) 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 (9). 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.

Using a quadrature rule involving the OPs (r q , c q ), q ∈ [1, N OP ], to approximate the continuous integral (9a), we end up with the discrete version min m:R 2 →R M N OP q=1 ξ q F (q) i (m(r q , c q )) 1≤i≤n (10a) subject to (r q , c q ) ≤ Am(r q , c q ) ≤ u(r q , c q ), q ∈ [1, N OP ],

(10b)

N OP q=1 ξ q F (q) j (m(r q , c q )) ≤ S j , j ∈ [1, d]. (10c)
4.1 Modeling the engine maps

In the formulations ( 9) and ( 10), minimization has to be carried out with respect to a function m(., .) : R 2 → R M , which is a continuous object. In order to fully discretize the problem, we have to consider a parameterization of the engine maps that must be flexible enough to model the very different shapes of engine map surfaces (Fig. 7) and that should not require too many parameters to limit the number of unknowns in the optimization process. LoLiMoT models [START_REF] Nelles | Nonlinear system identification: from classical approaches to neural networks and fuzzy models[END_REF][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) = λ=1 m λ (r, c)Φ λ (r, c), ( 11a
)
with

m λ (r, c) = θ 0,λ + θ r,λ r + θ c,λ c. ( 11b
)
The weights Φ λ (r, c) are normalized Gaussian functions defined as

Φ λ (r, c) = N λ (r, c) κ=1 N κ (r, c) , ( 12a
)
with

N λ c) = exp ⎨ ⎩ - 1 2α ⎡ ⎣ r -r 0 λ σ r λ 2 + c -c 0 λ σ c λ 2 ⎤ ⎦ ⎫ ⎬ ⎭ . ( 12b
)
Their purpose is to control the degree of smoothness of the global surface. α is a global smoothing parameter. As far as the coefficients θ 0,λ , θ r,λ and θ c,λ are concerned, each of them is a vector of size M. Therefore, the continuous function m(., .) has been replaced by a vector

p = {θ 0,λ , θ r,λ , θ c,λ } 1≤λ≤ ( 13
)
whose size is N = 3M . Substituting ( 11) into (10) brings us back to problem (1).

The advantage of the LoLiMoT representation (11) lies in its ability to adaptively refine 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. An example of the LoLiMoT parameterization in the context of engine tunings is displayed in Fig. 7.

Application of map optimization on a real case

The map optimization has been applied on the turbo-diesel engine application described in Sect. 3.3. The optimization problem is here defined as a constrained singleobjective optimization problem on the whole driving cycle whereas a two-objective formulation was used for the local optimization problem at a fixed OP. The cumulated particulate emissions on the cycle are minimized and the cumulated CO 2 , CO, HC and NO x emissions are constrained to remain under thresholds (given in Table 1). The unknowns of the optimization problem are the 378 LoLiMoT parameters describing the 6 initial engine maps. Linear constraints are added to take into account the limits of the physical domain for each point of the engine maps. The obtained results are listed in Table 1 and the optimized engine maps are displayed in Fig. 8. These engine maps have been run at the test bench and particulate emissions have been measured, those measures are displayed in Fig. 9. Table 1 shows that the objectives of optimization are reached: the particulate emissions are minimized whereas all the other cumulated engine responses remain under the constraint thresholds. Moreover the chosen map parametrization allows us to match the constraints, which is crucial for the feasibility of the engine tunings. The comparison of the particulate emissions obtained at test bench with the initial maps and the optimized on Fig. 9 illustrates the obtained gain: in the area covered by the optimized driving cycle (under 3000 rd/mn and under 10 bar), the emissions are obviously minimized. In the other parts of the operating domain, the differences are the result of the smoothing constraints applied on the maps introduced in the optimization formulation by means of second derivatives minimization. In contrast to successive local OP optimizations, the global optimization offers the flexibility to increase slightly the particulate emissions locally if it is necessary to match the smoothing constraints. It is the trend we observe on Fig. 9 in the area located around 1750 rd/mn-7 bar.

Conclusions

Throughout this contribution, we have demonstrated the relevance of the two approaches proposed for engine calibration. By extending the standard local singleobjective formulation in two directions, namely, the local multi-objective formulation and the global single-objective formulation, we seek to provide engine experts with a more powerful tool that enables them to achieve better results.

The technical challenges we had to face were solved by existing dedicated algorithms. For the local multi-objective approach, the constrained multi-objective optimization was successfully carried out by the MO-CMA-ES method, which we slightly modified to handle inequality constraints. The performance of this method on the considered real case application is very encouraging, in the sense that the Pareto front obtained gives worthy information on the possible compromises between antagonist engine responses, such as particulate emissions and NO x emissions of a diesel engine. As far as the global single-objective approach is concerned, the fact that it rests upon a cycle optimization of the engine maps instead of a pointwise optimization allows us to avoid the cumbersome step of map smoothing in classical engine calibration process. The numerical results testify to a notable gain in efficiency thanks to intrinsic smoothing constraints introduced via the LoLiMoT parameterization.

An upcoming study for the map optimization will be to apply a multi-objective approach as for the OP optimization. Here, the difficulty lies in the large size of the parameter space (500 parameters to describe all the engine maps). Another 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.

Fig. 1

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

Fig. 3

 3 Fig. 3 An application of CMA-ES method on the analytical function f (x, y) = x 2 +y 2 10 -1 2 (cos( 2π 5 x) + cos( 2π 5 y)) + 3. From top left to bottom right, the population at iterations 1, 5, 8 and 20 are represented. The ellipse represents the 95%-isovalue of the Gaussian probability density of the mutation operator with the adaptive covariance matrix

Fig. 4

 4 Fig. 4 Contributing dominance hyper-volumes for the evaluation of the diversity criterion

Fig. 6

 6 Fig. 6 Pareto fronts (crosses) of bi-objective optimization of particulate emissions and NO x emissions (in g/h) 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 indicate solutions obtained by single-objective optimizations (weighted sum of the two objectives) associated with different initial points and different penalty weights

Fig. 7

 7 Fig. 7 of LoLiMoT parameterization 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

Fig. 8

 8 Fig. 8 Optimal engine maps by the map optimization on the NEDC driving cycle. 378 LoLiMoT 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

Table 1

 1 Results of map optimization on a diesel engine application

		Initial maps	Optimized maps	Upper bounds	Error bar
	Particulate (g/km)	0.128	0.090	-	0.003
	NO x (g/km)	0.219	0.240	0.241	0.025
	HC (g/km)	0.137	0.150	0.150	0.014
	CO (g/km)	0.815	0.847	0.896	0.039
	CO 2 (g/km)	168.89	164.40	168.89	1.055

Acknowledgements

The authors are grateful to thank Yohan Bentolila, Michel Castagné, Fabien Chaudoye, Roland Dauphin, Adrien Hallé and Frédéric Nicolas for introducing us to the engine calibration problem and for supplying us with the real dataset used in this study. We also wish to thank the two anonymous referees for their constructive comments, which helped us to improve this paper.