Manh Quan Nguyen

Olivier Sename

Luc Dugard

Actuator fault estimation based on a switched LPV extended state observer

Keywords: Actuator fault estimation, switched LPV observer, dwell time

Introduction

Fault estimation is a step in Fault Diagnosis and plays a key role in designing a fault tolerant control. Many different approaches have been developped to estimate a fault which can be either actuator or sensor malfunction. Let us mention some classical methods, based on the parity space theory [START_REF] Gertler | Fault detection and isolation using parity relations[END_REF] to generate the residues and approximate the fault, or the bank of observers approach [START_REF] Isermann | Fault-diagnosis applications: model-based condition monitoring: actuators, drives, machinery, plants, sensors, and faulttolerant systems[END_REF], or the sliding mode observers [START_REF] Edwards | Sliding mode observers for fault detection and isolation[END_REF]. Recently, [START_REF] Hamdi | Fault detection and isolation in linear parameter-varying descriptor systems via proportional integral observer[END_REF] used an Unknown Input Proportional Integral Observer for actuator fault detection and estimation but this method is more oriented for constant faults. To deal with the time-varying faults, [START_REF] Zhang | Adaptive observer-based fast fault estimation[END_REF] presented a method using a Fast Adaptive Fault Estimation (FAFE) methodology based on an adaptive observer. Therein, the authors solved the problem with a regular LTI system without considering the disturbances. Then, [START_REF] Rodrigues | Actuator fault estimation based adaptive polytopic observer for a class of lpv descriptor systems[END_REF] proposed an adaptive polytopic unknow input observer for time-varying fault estimation, for a class of descriptor LPV systems. Besides, several works have been done for fault estimation and fault tolerant control for LPV systems. [START_REF] Abdullah | Sensor-fault-tolerant control for a class of linear parameter varying systems with practical examples[END_REF] proposes the design of a low-order LPV observer to estimate the unmeasured states of the system and to estimate the sensor faults The material in this paper was partially presented at 1st IFAC Workshop on Linear Parameter Varying systems, Grenoble, France, October, 7-9 2015 Email addresses: manh-quan.nguyen@gipsa-lab.fr (M.Q. Nguyen), olivier.sename@gipsa-lab.fr (O.Sename), L.Dugard@gipsa-lab.fr (L.Dugard).

for a class of uncertain LPV systems. An LMI-based poleplacement robust LPV estimator is presented in [START_REF] Patton | An lpv poleplacement approach to friction compensation as an ftc problem[END_REF], an interval observer for LPV systems in [START_REF] Montes De Oca | Robust fault detection based on adaptive threshold generation using interval lpv observers[END_REF], an LPV sliding mode observer in [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF]. Moreover, a new approach in [START_REF] Shi | Active fault tolerant control of lpv descriptor systems based on extended state observers[END_REF] considers the fault element as a state of the augmented system and an LPV extended observer is designed to estimate, at the same time, the state and the fault of system where constant faults are considered. Regarding to the actuator fault estimation problem, along with the additive model, the fault can be also written in a multiplicative form by using fault effect factors which are assumed to be constant or slow-varying. Besides, inspired by the fact that the control input is known, the considered system can be rewritten as an extended LPV system while considering the control inputs as scheduling paramters. Moreover, it has to be noticed that when the control input is zero, the actuator fault information is not available. Therefore, to deal with this problem, the considered system is modeled as a switched LPV system . It is well known that for the switched system, the main challenge for a controller or observer design problem for this type of system is to deal with the stability analysis, especially for the continuous time switched system due to the discontinuities of the Lyapunov function at switching instants. A lot of studies have been investigated for such a problem during recent decades. Let us mention firstly the LTI case. The stability of continuous-time LTI switched systems has been addressed in [START_REF] Michael | Multiple lyapunov functions and other analysis tools for switched and hybrid systems[END_REF], [START_REF] Hespanha | Stability of switched systems with average dwell-time[END_REF], [START_REF] Liberzon | Switching in systems and control[END_REF], [START_REF] Jose | Stability and stabilization of continuous-time switched linear systems[END_REF]. Therein, the stability under dwell time, average dwell time constraints using multiple Lyapunov functions are studied. The main idea is to ensure that the Lyapunov function is non-increasing at switching instants or to allow some short-time increasing but the dwell time between two consecutive switching instants is sufficient large in order to compensate for possible increase of Lyapunov functions. In particular, [START_REF] Jose | Stability and stabilization of continuous-time switched linear systems[END_REF] provides an efficient way to get the stability for the switched LTI system under minimum dwell-time by mean of a family of quadratic Lyapunov function. However, such an approach is not easy to extend to the switched uncertain systems because of nonconvex dependence on the system matrices. To overcome this, recently, [START_REF] Briat | Affine characterizations of minimal and mode-dependent dwell-times for uncertain linear switched systems[END_REF] presented the Lyapunov looped-functionals approach, a new type of functionals leading to stability conditions that are affine in the system matrices. [START_REF] Briat | Convex conditions for robust stabilization of uncertain switched systems with guaranteed minimum and mode-dependent dwell-time[END_REF] proposes an alternative solution using lifted conditions which are convex in the system matrices and shown to be equivalent to the nonconvex conditions proposed in [START_REF] Jose | Stability and stabilization of continuous-time switched linear systems[END_REF]. Another possibility is presented in [START_REF] Allerhand | Robust stability and stabilization of linear switched systems with dwell time[END_REF], [START_REF] Liron | Robust control of linear systems via switching. Automatic Control[END_REF] where a piecewise linear in time, quadractic Lyapunov form function is used to derive convex conditions for both stability and stabilization problems with dwell time. For the class of switched LPV systems, most works have considered the control problem but very few the observation one. [START_REF] Lu | Switching lpv control designs using multiple parameter-dependent lyapunov functions[END_REF] proposed a switching LPV controller based on the multiple Lypunov functions under hysteresis switching and average dwell time constraint. [START_REF] Zhang | Model reduction for switched lpv systems with average dwell time. Automatic Control[END_REF] dealt with a model reduction problem for switched LPV system in which the stability was derived also by using the average dwell time technique. The main purpose of this paper is to propose a methodology to estimate multiplicative faults for the actuators. An actuator time-varying fault is considered in the form of actuator power loss. Then, effectiveness factors are used to model the efficiency of actuators. The paper contributions emphasize on the following aspects. Firstly, the system in the presence of actuator faults is modeled in the form of a switched LPV system by considering the control inputs as scheduling parameters. Secondly, a switched LPV extended observer is designed to estimate both the actuator faults and the system state. The stability of the switched LPV observer will be guaranteed using the dwell time constraint. To this aim, a non-increasing piece-wise linear in time Lyapunov function is assigned for each subsystem. The H ∞ performance is used to minimize the L 2 gain from the disturbance to the estimation error. The rest of this paper is organized as follows: the next section presents the problem formulation. Section 3 gives some preliminaries for this work. Section 4 gives a full description for multiplicative faults estimation based on the switched LPV observer approach. In Section 4, two numerical examples are presented to illustrate the effectiveness the proposed approach. Finally, some conclusions are drawn in the section 5.

Problem Formulation

System definition

Consider a continuous time linear invariant system:

ẋ(t) = Ax(t) + B 2 u(t) + B 1 w(t) y(t) = Cx(t) (1)
where x(t) ∈ R n , u(t) ∈ R m , w(t) ∈ R q and y(t) ∈ R p are the state, the control input, the input disturbance and the measured output vectors, respectively. Matrices A ∈ R n×n , B 2 ∈ R n×m , B 1 ∈ R n×q ,C ∈ R p×n are known matrices of appropriate dimensions. In this work, only actuator faults are considered. Assume that the system (1) is in the faulty actuator situation, e.g loss of actuator power (Fig. 1). These actuator faults are modeled in a multiplicative representation. In fact, denoting that ūi is the output of i th faulty actuator , then:

ūi (t) = λ i u i (t) (2)
where λ i stands for the efficiency coefficient of the i th actuator and λ i is assumed to be constant. λ i = 1 implies that the i th actuator is fault-free, λ i = 0 mean that the i th actuator is in total failure. And 0 < λ i < 1 represents the fact that the fault of the ith acutator is a partial loss of control effectiveness, e.g if λ i = 0.8, the i th actuator loses 20% of its effectiveness.

It is worth noting also that even if λ i is assumed to be constant, the corresponding additive fault magnitude on ith actuator given by f i (t) = (1-λ i )u i (t) is a time varying signal and depends on the value of the control input u i (t). Thanks to the multiplicative representation, the information about the actuator fault λ i is considered as constant or slow-varying and λi = 0 which will be used later for the extended system.

In the presence of the actuator faults, the input matrix B 2 becomes B 2 Λ where Λ ∈ R m×m is a diagonal matrix representing the impact factors of the actuator faults, i.e:

Λ = diag([λ 1 λ 2 ...λ m ])
Then, the following equation gives the representation of the faulty system subject to the actuator faults:

ẋ(t) = Ax(t) + B 2 Λu(t) + B 1 w(t) y(t) = Cx(t) (3)
Since Λ is a diagonal matrix, (3) can be rewritten as:

ẋ(t) = Ax(t) + B 2        u 1 u 2 . . . u m        U        λ 1 λ 2 . . . λ m        λ +B 1 w(t) (4)
The objective of this work is to estimate the vector λ = [λ 1 λ 2 ... λ m ] T where λ i represents the effectiveness factor of the actuator i th . The estimation, based on an extended switched observer, is presented in the sequel.

Switched LPV system

Thanks to the multiplicative fault representation, considering λi = 0, the system (4) can be augmented as follows:

             ẋ λ = A B 2 U 0 0 x λ + B 1 0 w y = [C 0] x λ (5)
where the dimension of the augmented system is n + m.

In order to estimate the vector λ , the system (5) must be observable or at least detectable. However, it can be seen that if the control input u i (t) = 0, then the fault information of the ith actuator λ i in (5) becomes unobservable. It makes the problem unfeasible in the observer design step. Thus, in this proposed method, the persistent excitation condition on the control input u is required. Moreover, an interesting remedy is to take into account the effect of the sign of u i (t) in the observer synthesis step, i.e the observer will be designed for different cases where the control inputs are positive and negative. To account for the change the sign of u i (t), the system will be rewritten as a switched system and a switched observer will be designed in the sequel.It is well known that in a real mechatronic system, the actuators always admit some physical constraints, i.e the control input u(t) satisfies the following saturation constraint:

u(t) ∈ U = u ∈ R m |u i ≤ u i ≤ u i (6)
where u i , u i are the lower, upper bounds of the ith actuator u i . Moreover, since the control input u(t) is known, the system (3) can be represented as an LPV system by choosing u(t) as the vector of scheduling parameters. Let us rewrite u i (t) = |u i (t)|sign(u i (t)), and denote ρ i (t) = |u i (t)| as a scheduling parameter. Then,

u i (t) = ρ i (t)sign(u i (t)) = ρ i (t) i f u(t)) ≥ 0 -ρ i (t) i f u(t) < 0 (7)
Thus, the scheduling parameter vector ρ = [ρ 1 ρ 2 ... ρ m ] T is assumed to read for the following constraint:

ρ(t) ∈ Ω = {ρ ∈ R m |ε ≤ ρ i ≤ ρ i } (8)
where ε > 0, ρ i are the bounds of the parameter ρ i .

The faulty system ( 4) is now rewritten as the following LPV system:

ẋ(t) = Ax(t) + B σ (ρ)λ + B 1 w(t) (9)
where

B σ (ρ) = B 2        ρ 1 sign(u 1 ) ρ 2 sign(u 2 ) . . . ρ m sign(u m )        and λ = [λ 1 λ 2 ... λ m ] T .
The system (9) can then be augmented into the following form:

                     ẋ λ = A B σ (ρ) 0 0 A e (ρ) x λ + B 1 0 w y = [C 0] C e x λ (10)
The system [START_REF] Gertler | Fault detection and isolation using parity relations[END_REF] is actually a switched system where σ (t) is the switching rule that depends on the value of the function sign(u i (t)) and takes values on the discrete set {1, 2, ..., M}. M = 2 m is the number of "modes" that compose the overall switched dynamics. Indeed, the matrix B σ (ρ) switches between different matrices and the switching moments depend on the sign of the actuators sign(u i (t)). For example, if we consider a system with 2 actuators (m = 2), the matrix B σ (ρ) belongs to the following set:

B 2 ρ 1 ρ 2 , B 2 -ρ 1 ρ 2 , B 2 ρ 1 -ρ 2 , B 2 -ρ 1 -ρ 2
Therefore, we can rewrite the system (10) in the following switched LPV system form:

             ẋ λ = A e,σ (ρ) x λ + B 1 0 w y = C e x λ (11)
where A e,σ (ρ) ∈ R (n+m)×(n+m) switches among the subsystems {A e,1 (ρ), A e,2 (ρ), ..., A e,M (ρ)}.

Problem Statement

The actuator fault estimation problem is recast into that of the parameter vector λ . To this aim, a switched LPV Extended State Observer (ESO) will be designed for the switched LPV system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF].

Then, the following detectability condition is assumed for the design of the ESO of extended LPV system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF]:

A1: rank sI -A e,σ (ρ)
C e = n + m, s > 0, ∀ρ ∈ Ω (12)
Remark 1 The observability for the LPV system is not a trivial problem because the variation of the parameter ρ(t) is considered in computing the observability matrix. Several definitions of observability for the LPV systems have been given in the literature such as Quadratic detectability [START_REF] Wu | Control of linear parameter varying systems[END_REF] in the sense of Lyapunov function and Structural Observability [START_REF] Toth | Modeling and identification of linear parameter-varying systems[END_REF] i.e the observability matrix is full rank in function sense of ρ(t) but one can lose the observability in some frozen point.

Under the Assumption A1, the problem consists now in designing a switched LPV observer in order to estimate the effectiveness factors of the actuator faults. Therefore, considering the switched LPV system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF], the following switched LPV extended observer is proposed to estimate the system's state and the effectiveness coefficient vector λ :

             ẋ λ = A e,σ (ρ) x λ + K σ (.)(y -ŷ) ŷ = C e x λ (13)
From ( 11) and ( 13), the estimation error e(t) is calculated by:

ė = ėx ėλ = A e,σ (ρ) e x e λ -K σ (.)(y -ŷ) + B 1 0 w (14)
Or equivalently:

ė = ėx ėλ = (A e,σ (ρ) -K σ (.)C e ) e x e λ + B 1e w (15)
where B 1e = B 1 0 and K σ (.) is the observer gain which has to be determined. Problem definition: Let consider the switched LPV system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF]. The system ( 13) is said to be an ESO for the system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF] if it satisfies the following conditions:

• when w(t) ≡ 0, the estimation errors ( 15) is asymptotically stable

• when w(t) = 0, the estimation error satisfies the following L 2 -induced gain performance criterion:

min γ s.t sup w =0,w∈L 2 z 2 w 2 < γ (16)
where . 2 stands for L 2 norm and z is given by: z = e if both state and fault estimation error are minimized z = [0 I m ]e if the fault estimation error only is to be minimized Therefore, the main problem now is to design the switched LPV observer [START_REF] Hespanha | Stability of switched systems with average dwell-time[END_REF]. As mentioned previously, regarding to a switched observer, the main challenge is to ensure the stability requirement. To deal with, the main result presented in the the next section is a switched LPV observer with dwell-time. In order to ensure the stability of the switched observer, a piecewise linear in time Lyapunov function in quadratic form is used. The applied Lyapunov function is non-increasing at switching instants and is assigned seperately to each subsystem. During the dwell time, this function varies piecewise linearly in time and after the dwell time, it becomes time invariant. Such a Lyapunov function allows to derive the stability conditions for the switched LPV system. Moreover, the minimization of L 2 induced gain is performed in order to minimize the effect of disturbance on the estimation error.

3 Preliminaries on the stability of switched LPV system 3.1 Recall for the LTI case This section is devoted to recall some results on the stability analysis for the continuous time, switched system by using the multiple Lyapunov function. The interested readers can refer to several works of [START_REF] Hespanha | Stability of switched systems with average dwell-time[END_REF], [START_REF] Liberzon | Switching in systems and control[END_REF], [START_REF] Jose | Stability and stabilization of continuous-time switched linear systems[END_REF], [START_REF] Allerhand | Robust stability and stabilization of linear switched systems with dwell time[END_REF]. Let us consider the following switched LTI system:

ẋ(t) = A σ (t) x(t), x(0) = x 0 , (17)
where σ (t) is the switching signal and

A σ (t) ∈ {A 1 , ..., A M } , A i ∈ R n×n , i = 1...M.
Obviously, this model is discontinuous w.r.t A σ (t) since this matrix jumps instanteously from A i to A j for i = j. It is shown in [START_REF] Jose | Stability and stabilization of continuous-time switched linear systems[END_REF] that the stability of this switched LTI system is ensured under a minimum dwell time i.e, if, for some T > 0, there exists a family of symmetric and positive Lyapunov matrices {P 1 , ...P M }, such that:

A i P i + P i A i < 0, ∀i = 1, ...M e A i T P j e A i T -P i < 0, ∀i = j = 1, ...M (18)
Then the system is globally asymptotically stable for a dwell time greater than or equal to T . However, since the condition ( 18) is non convex in A i , it is not easy to generalize to a system with uncertainties or to an LPV system. More recently, an alternative solution was proposed in [START_REF] Allerhand | Robust stability and stabilization of linear switched systems with dwell time[END_REF] to deal with both nominal and uncertain systems using a piecewise linear in time Lyapunov quadratic function. This Lyapunov function is non-increasing at the switching instants and is assigned to each subsystem. This method could be more conservative than the one presented in [START_REF] Jose | Stability and stabilization of continuous-time switched linear systems[END_REF], however it provides an efficient way to deal with uncertainties and LPV systems. Inspired by this approach, one will use this kind of Lyapunov function in order to guarantee the stability of the switched LPV system in the next section.

Stability condition for switched LPV system

Now, the LPV case is taken into account with the following switched LPV system:

ẋ(t) = A σ (t) (ρ)x(t), x(0) = x 0 , (19)
where

σ (t) is the switching signal, ρ = [ρ 1 ρ 2 . . ρ m ] T ∈ Ω ⊂ R m is the scheduling parameter vector and A σ (t) (ρ) ∈ {A 1 (ρ), ..., A M (ρ)} , A i ∈ R n×n , i = 1...M.
Let consider the a piecewise linear in time Lyapunov candidate function V (x(t)) = x (t)P σ (t) x(t) where P σ (t) is given as:

P σ (t) (t) =          P i,k + (P i,k+1 -P i,k ) t -τ s,k T /G := Pi,k , i f t ∈ [τ s,k , τ s,k+1 ) P i,G i f t ∈ [τ s,G , τ s+1,0 ) P i 0 ,G i f t ∈ [0, τ 1 ) (20)
where i = 1, ...M with M is number of subsystems, i 0 = σ (0). τ 1 , τ 2 , ... are the switching instants. T is the dwell time satisfying τ s+1τ s ≥ T , and τ s,k = τ s + k(T /G) for k = 0, ...G , τ s = τ s,0 , τ s,G = τ s + T ). P i,k are symmetric matrices of compatible dimensions, where G is an integer that may be chosen a priori. The switched LPV system [START_REF] Rodrigues | Actuator fault estimation based adaptive polytopic observer for a class of lpv descriptor systems[END_REF] is asymptotically stable if all subsystems are stable and the Lyapunov function V (x(t)) is non-increasing at the switching instants, i.e:

A i (ρ) P σ (t) + P σ (t)A i (ρ) + Ṗσ (t) < 0, ∀i = 1, ...M (21) V (x(τ k )) ≤ V (x(τ - k )), τ k is the switching instant (22)
In order to transform these conditions into convex conditions, the following lemma is used: Lemma 1 ( [START_REF] Boyarski | Time-convexity and time-gain-scheduling in finite-horizon robust hinf-control[END_REF]) Assume that for some interval t ∈ [t 0 t f ], and δ = t ft 0 , there exist two symmetric positive matrices P 1 , P 2 of appropriate dimensions that satisfy the following conditions:

P 2 -P 1 δ + P 1 A + A P 1 < 0, P 2 -P 1 δ + P 2 A + A P 2 < 0 (23)
Then, for the system ẋ = Ax, the Lyapunov function V (t) =

x (t)P(t)x(t), where P(t

) = P 1 + (P 2 -P 1 ) t -t 0 δ , is strictly decreasing over the interval [t 0 t f ].
Then, by using the formula of P σ (t) in [START_REF] Shi | Active fault tolerant control of lpv descriptor systems based on extended state observers[END_REF] and applying the polytopic approach for the LPV system, the following theorem gives stability conditions for the switched LPV system (19): Theorem 2 Consider the switched LPV system [START_REF] Rodrigues | Actuator fault estimation based adaptive polytopic observer for a class of lpv descriptor systems[END_REF], if there exists a collection of matrices P i,k > 0, k = 0, ...G, i = 1, ...M, of appropriate dimensions and G is a prescribed integer, such that for all i = 1, ..M and j = 1, ...N, (N = 2 m : number of the vertices of the polytope), the following LMIs hold:

(P i,k+1 -P i,k ) T /G + A ( j) i P i,h + P i,h A ( j) i < 0, (24)
k = 0, ...G -1, h = k, k + 1 A ( j) i P i,G + P i,K A ( j) i < 0 (25) P i,G -P l,0 ≥ 0 ∀l = 1, ..i -1, i + 1, ...M. ( 26
)
then the switched LPV system ( 19) is asymptotically stable.

Proof: The switched LPV system ( 19) is asymptotically stable if:

V (x) = x(t) (A i (ρ) P σ (t) + P σ (t)A i (ρ) + Ṗσ (t))x(t) < 0 V (x(τ k )) ≤ V (x(τ - k )), τ k is the switching instant (27)
It is assumed that at switching instant τ k , the system switches from A i (ρ) to A l (ρ), i.e we can write the Lyapunov function at instant τ k as:

V (x(τ - k )) = x(τ k )P i,G x(τ k ), V (x(τ k )) = x(τ k )P l,0 x(τ k ) (28
) Then the non-increasing Lyapunov function condition holds if P i,G -P l,0 ≥ 0 which is actually (27). Now, we apply Lemma 1 for each switching instant τ s . Indeed, during the dwell time, considering the time interval [τ s,k , τ s,k+1 ), the Lyapunov matrix P σ (t) changes linearly from P i,k to P i,k+1 . Then from Lemma 1 with δ = τ s,k+1τ s,k = T /G and the polytopic approach for LPV system, (27) is true if (24) holds. After the dwell time and before the next switching instant, i.e t ∈ [τ s,G , τ s+1,0 ), the Lyapunov matrix P σ (t) becomes time invariant P i,G , then (27) holds if A i (ρ) P i,G + P i,G A i (ρ) < 0 that is equivalent to (25) thanks to polytopic approach for LPV system. Remark 1 Such a choice of piecewise Lyapunov function as mentioned above gives an efficient way to deal with the stability analysis problem for the class of switched LPV system. Indeed, the Lyapunov matrix P depends on time but not on parameter ρ(t), then the derivative of P depends on the varition of the piecewise Lyapunov function during the dwell time and not on ρ(t). An alternative solution is to use the Lifted Conditions with Lyapunov dependant parameter function to convexify the stability condition as in [START_REF] Briat | Convex conditions for robust stabilization of uncertain switched systems with guaranteed minimum and mode-dependent dwell-time[END_REF]. Extension to the observation problem: Consider the system [START_REF] Rodrigues | Actuator fault estimation based adaptive polytopic observer for a class of lpv descriptor systems[END_REF] and y = C σ (t) (ρ)x ∈ R p is the output vector of the system. Then, the states of the system (19) are reconstructable if there exist a matrix P σ (t) (t) ∈ R n×n of form of (20) and a function K σ (t) (t) ∈ R n×p , such that:

P σ (t)[A σ (ρ) -K σ (t)C σ (ρ)] +[A σ (ρ) -K σ (t)C σ (ρ)] P σ (t) + Ṗσ (t) < 0 (29)
The following theorem allows to solve the problem in (29):

Theorem 3 Consider the switched LPV system [START_REF] Rodrigues | Actuator fault estimation based adaptive polytopic observer for a class of lpv descriptor systems[END_REF], if there exists a collection of matrices P i,k > 0,Y i,k , k = 0, ...G, i = 1, ...M, of appropriate dimensions and G is a prescribed integer, such that for all i = 1, ..M and j = 1, ...N, (N = 2 m : number of the vertices of the polytope), the following LMIs hold:

(P i,k+1 -P i,k ) T /G + A ( j) i P i,h -C i Y i,h + P i,h A ( j) i -Y i,h C i < 0,( 30
) k = 0, ...G -1, h = k, k + 1 A ( j) i P i,G -C i Y i,G + P i,K A ( j) i -Y i,G C i < 0(31) P i,K -P l,0 ≥ 0 ∀l = 1, ..i -1, i + 1, ...M. (32)
then the states of the system (19) are reconstructable. And K σ (t) (t) is given by:

K σ (t) (t) = P σ (t) (t) -1 Y σ (t) (t) =        P-1 i,k Ŷi,k , i f t ∈ [τ s,k , τ s,k+1 ) P -1 i,G Y i,G , i f t ∈ [τ s,G , τ s+1,0 ) P -1 i 0 ,G Y i 0 ,G , i f t ∈ [0, τ 1 ) (33) where: Y σ (t) (t) =          Y i,k + (Y i,k+1 -Y i,k ) t -τ s,k T /G := Ŷi,k , i f t ∈ [τ s,k , τ s,k+1 ) Y i,G i f t ∈ [τ s,G , τ s+1,0 ) Y i 0 ,G i f t ∈ [0, τ 1 ) (34) .
Proof: The proof can be infered easily from the theorem 2.

4 Switched LPV observer under a dwell-time constraint In this section, the design of the switched observer using the dwell time notion is presented. Considering the switched system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF], as presented in section 2. The following switched LPV extended observer is proposed to estimate the system's state and the effectiveness coefficient vector λ :

ẋ λ = A e,σ (ρ) x λ + K σ (t)(y -ŷ) (35) ŷ = C e x λ
and the estimation error is calculated by:

ė = ėx ėλ = (A e,σ (ρ) -K σ (t)C e ) e x e λ + B 1e w ( 36
)
where B 1e = B 1 0 and K σ (t) is the observer gain which has to be determined. Now, regarding to the Problem definition in section 2, and thank to the Bounded Real Lemma and the Polytopic approach, the following theorem allows to compute the observer gain. Theorem 4 Consider the switched system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF] and the switched extended observer (35). If there exists a collection of matrices P i,k > 0,Y i,k , k = 0, ...G, i = 1, ...M, of appropriate dimensions and G is a prescribed integer, such that for all i = 1, ..M and j = 1, ...N, (N = 2 m : number of the vertices of the polytope), the following LMIs hold:

     (P i,k+1 -P i,k ) T /G + He[P i,h A ( j) e,σ -Y i,h C e ] B 1e P i,h -γ 2 I I 0 -I      < 0 (37) for k = 0, ...G -1, h = k, k + 1,     A ( j) e,σ P i,G -C e Y i,G + P i,G A ( j) e,σ -Y i,G C e B 1e P i,G -γ 2 I I 0 -I     < 0 (38) P i,G -P l,0 ≥ 0 ∀l = 1, ..i -1, i + 1, ...M. ( 39
) then K σ (t) (t) = P σ (t) (t) -1 Y σ (t) (t) =        P-1 i,k Ŷi,k , i f t ∈ [τ s,k , τ s,k+1 ) P -1 i,G Y i,G , i f t ∈ [τ s,G , τ s+1,0 ) P -1 i 0 ,G Y i 0 ,G , i f t ∈ [0, τ 1 ) (40)
is the gain of the extended observer (35) and the error estimation asymptotically converges to zero for a dwell time of T , where:

Y σ (t) (t) =          Y i,k + (Y i,k+1 -Y i,k ) t -τ s,k T /G := Ŷi,k , i f t ∈ [τ s,k , τ s,k+1 ) Y i,G i f t ∈ [τ s,G , τ s+1,0 ) Y i 0 ,G i f t ∈ [0, τ 1 ) (41)
Proof: Let V (e(t)) = e (t)P σ (t)e(t) be the Lyapunov candidate function for the estimation error system (36) where P σ (t) is a piecewise linear in time Lyapunov matrix defined as in [START_REF] Shi | Active fault tolerant control of lpv descriptor systems based on extended state observers[END_REF]. From the Bounded Real Lemma, the condition ( 16) is satisfied if the following condition holds:

V + e eγ 2 w w < 0

i.e:

    He[P σ (t)(A e,σ (ρ) -K σ (t)C e )] + Ṗσ (t) P σ (t)B 1e I * -γ 2 I 0 * * -I     < 0 (43)
From the formula of P σ (t) in ( 20), ( 43) is satisfied if

     He[P i,h A e,σ (ρ) -Y i,h C e ] + (P i,k+1 -P i,k ) T /G B 1e P i,h -γ 2 I I 0 -I      < 0 (44) holds for h = k, k + 1, i = 1, 2...M, k = 0, ...G -1. and     A e,σ (ρ) P i,G -C e Y i,G + P i,G A e,σ (ρ) -Y i,G C e B 1e P i,G -γ 2 I I 0 -I     < 0 (45) hold for i = 1, 2...M.
The equation (44) guarantees that the Lyapunov function V σ (t) decreases and that (42) holds during the time intervals t ∈ [τ s,0 , τ s,G ). The LMIs (45) ensure that V σ (t) decreases and that (42) holds after the dwell time and before the next switching instant, i.e t ∈ [τ s,G , τ s+1,0 ). From the definition of P σ (t), consider that at instant τ k , that the system switches from the mode i to the mode l. To guarantee the Lyapunov function is non-increasing at the switching instants, we must ensure:

P i,G -P l,0 ≥ 0 ∀l = 1, ..i -1, i + 1, ...M. (46)
Now, in order to resolve the LMIs in ( 44), (45), we apply the polytopic solution for the LPV system where the polytope is given by Ω ρ = ρ ρ and we obtain the LMIs in (37), (38).

Extension with a decay rate on the convergence:

The next result extends the previous one, imposing a prefixed decay rate on the convergence of the estimation error.

Theorem 5 Consider the switched system [START_REF] Hamayun | An lpv fault tolerant control scheme using integral sliding modes[END_REF] and the switched observer (35). If there exists a collection of matrices P i,k > 0,Y i,k , k = 0, ...G, i = 1, ...M, of appropriate dimensions, G is a prescribed integer, and a positive scalar β such that for all i = 1, ..M and j = 1, ...N, (N = 2 m : number of the vertices of the polytope), the following LMIs hold:

     (P i,k+1 -P i,k ) T /G + He[P i,h A ( j) e,σ -Y i,h C e ] + 2β P i,h B 1e P i,h -γ 2 I I 0 -I      < 0 (47) for k = 0, ...G -1, h = k, k + 1,     A ( j) e,σ P i,G -C e Y i,G + P i,K A ( j) e,σ -Y i,G C e + 2β P i,G B 1e P i,G -γ 2 I I 0 -I     < 0 (48) P i,K -P l,0 ≥ 0 ∀l = 1, ..i -1, i + 1, ...M. ( 49
) then K σ (t) (t) = P σ (t) (t) -1 Y σ (t) (t)
is the gain of the extended observer (35) and the error estimation asymptotically converges to zero for a dwell time of T . Proof: The proof is similar to the last cases and is omitted here for the simplification. Proposition 1 The design of the switched observer (35) can be performed by solving the following optimization problem:

min P i,k ,Y i,k γ 2
subject to (47), ( 48), (49) and P i,k > 0 (50)

By solving this optimization problem, one can derive P σ (t),Y σ (t) and the observer gain is calculated by

K σ (t) = P σ (t) -1 Y σ (t).

Numerical Exemples

In this section, two different numerical examples will be presented in order to illustrate the effectiveness of the proposed switched LPV observer. Firstly, the present approach is applied for a semi-active suspension system that is actually a single input-multiple output (SIMO) system where the actuator is the semi-active damper. Then, an example with a multiple inputs multiple outputs (MIMO) system is used to emphasize the interesting of the presented apporach.

5.1 Actuator fault estimation for semi-active suspension system The simulation is performed in a small pilot SOBEN Car equipped with four electrorheological semi-active dampers (see Fig. where F sa is the semi-active damper force.

The physical parameters characterize the quarter SOBEN car is given in the Table 1. Let us consider now F sa as the control input u = F sa , and assume that a fault occured on the semi-active damper e.g an oil leakage which induces a lack of force modeled as:

F sa = λ F sa = λ u (52)
where F sa stands for the fault force expressed as a reduction of the nominal semi-active force and λ ∈ [0 1] is the oil leakage degree, e.g. λ = 0.8 means that the damping force will be of 80% of the nominal damper force F sa due to an efficiency loss of 20% .

Then the state space representation of the vertical dynamic using a quarter car model and taking into account a faulty semi-active damper, is given as follows: 3 is the output vector and λ stands for actuator fault.

ẋ = Ax + B 1 w + B 2 λ u (53) y = Cx where x = (z s , żs , z us , żs ) T ∈ R 4 is the state vector, w = z r is the input disturbance, u ∈ R is the control input, y = [z s - z us , żs -żus , z s ] = [z de f , żde f , z s ] ∈ R
A =        0 1 0 0 -k s m s 0 k s m s 0 0 0 0 1 k s m us 0 -k s +k t m us 0        , B 1 =        0 0 0 k t m us        ,B 2 =        0 -1 m s 0 1 m us       
, C = 1 0 -1 0; 0 1 0 -1; 1 0 0 0 . The semi-active damper force is assumed to be bounded by -10N ≤ F sa ≤ 10N. Thus, the scheduling parameter ρ(t) = |u(t)| = |F sa | is considered to belong to 0.001 ≤ ρ ≤ 10.

The quarter car model ( 53) is a SIMO system and can be rewritten in the form of the switched LPV system (10) where 2 subsystems A e,1 (ρ), A e,2 (ρ) are given by: Then, the switched LPV observer is designed using a dwell time constraint as presented in the last section. With dwell time T = 0.1s, a prescribed integer K = 1, and a decay rate β = 0.2, by solving the optimization problem (50) with LMIs (47,48, 49) (where the decay rate is taken into account) for the 2 modes corresponding to A e,1 (ρ), A e,2 (ρ) to obtain the matrices P i,k ,Y i,k , i = 1, 2 and k = 0, 1. Then the gain of the switched observer is calculated by K σ (t) = P σ (t) -1 Y σ (t) where P σ (t),Y σ (t) are given as in [START_REF] Shi | Active fault tolerant control of lpv descriptor systems based on extended state observers[END_REF] and (41).

A e,1 (ρ) = A B 2 ρ 0 1×4 0 ,A e,2 (
The road profile used in this test is a sinus profile as in Fig. 3. The control input and the scheduling parameter are given in Fig. 4. For the first scenario, it is assumed that an oil leakage occuring at t=5s causes a loss of 20% efficiency of damper force, the faulty damper remains 80% its healthy damper force (i.e λ = 0.8). It is noted that the switched observer is activated from t=0.5s. Fig. 5 shows the state estimation and estimation errors. It can be seen that the state of the system is well estimated. Fig. 6 demonstrates that the switched LPV observer allows to estimate the effectiveness factor λ of the actuator. Now, let us consider the second scenario where a gradual

Actuator faults estimation for a MIMO system

Let us consider the following MIMO system:

ẋ(t) = Ax(t) + B 2 Λu(t) + B 1 w(t) (54) y(t) = Cx(t) A = -1.25 1 0.1 -3 , B 2 = 5 1 3 10 , Λ = λ 1 0 0 λ 2 C = 1 1 1 -1 , B 1 = 0.5 1 , u = u 1 u 2 , x = x 1 x 2
The system is subject to the control inputs u 1 (t) = 20 sin(4πt), u 2 (t) = 30 sin(2πt) and the disturbance w(t).

Then, the scheduling parameters ρ 1 (t) = |u 1 (t)|, ρ 2 (t) = |u 2 (t)| are assumed to be bounded by: 0.001 ≤ ρ 1 ≤ 20, 0.001 ≤ ρ 2 ≤ 30. λ 1 , λ 2 are the effectiveness factors of the 2 control inputs u 1 , u 2 .

We have 2 actuators, so the system (54) is rewritten as a switched system with 4 subsystems, i.e A e,σ (ρ) is switched between four modes A e,1 (ρ), A e,2 (ρ), A e,3 (ρ), A e,4 (ρ) according to the switching signal σ (t) as follows:

               A e,σ (ρ) =   A B σ (ρ) 0 2×2 0 2×2   , B σ (ρ) ∈ B 2 ×      ρ 1 ρ 2   ,   -ρ 1 ρ 2   ,   ρ 1 -ρ 2   ,   -ρ 1 -ρ 2     
where ρ 1 , ρ 2 belong to the polytope:

Ω ρ = ρ 1 ρ1 ∨ ρ 2 ρ2 .
The Assumption A1 is also satisfied in this case. Then, in order to estimate the vector of effectiveness factors, the switched LPV observer is designed using the observer design procedure presented in Section 4. Let us choose: K = 1, T = 0.2, β = 1e -4. By solving the optimization problem (50) with LMIs (47,48, 49) for the 4 modes corresponding to A e,1 (ρ), A e,2 (ρ), A e,3 (ρ), A e,4 (ρ) to obtain the matrices P i,k ,Y i,k , i = 1, 2, 3, 4 and k = 0, 1, then the gain of the switched observer is calculated by K σ (t) = P σ (t) -1 Y σ (t). Fig. 8 shows disturbance w(t), the control inputs u 1 , u 2 , the correspondant scheduling parmeters ρ 1 , ρ 2 and the switching signal σ (t).

Conclusion

In this paper, actuator fault estimation problem has been proposed within the LPV approach. The actuator faults are modeled in a multiplicative way by using the effectiveness factors (λ i ∈ [0 1]). The fault estimation is based on a switched LPV observer. A non-increasing piece-wise linear time Lyapunov function is used to ensure the stability of the switched observer.The effectiveness of the proposed approach has been validated on a semi-active suspension system and an academic system. In the future work, some fault tolerant control strategies can be developped, based on the proposed fault diagnosis.

Fig. 1 .

 1 Fig. 1. System with actuator fault

  2) using the experimental data. Only the quarter car model is considered. In this model, the quarter vehicle body is represented by the sprung mass (m s ), the wheel and tire are represented by the unsprung mass (m us ). They are connected by a spring with the stiffness coefficient k s and a semi-active damper. The tire is modeled by a spring with the constant stiffness coefficient k t . As seen in the figure, z s (respectively z us ) is the vertical displacement around the equilibrium point of m s (respectively m us ) and z r stands for the road profile.

Fig. 2 .

 2 Fig. 2. Quarter-car vehicle model The dynamical equations of such a quarter vehicle model are given by:

  = [C 0 3×1 ] It can be easily obtained that the rank conditions in the Assumption A1 are satisfied: rank sI -A e,i (ρ) C e = 5 for i = 1, 2, all ρ ∈ Ω = [0.001 10].

Fig. 4 .

 4 Fig. 3. Road profile

4 Fig. 5 .Fig. 6 .

 456 Fig. 5. State estimation (above) and error estimation

Fig. 7 .

 7 Fig. 7. A gradual damper fault scenario fault is considered. The damper oil leaks slowly from t = 2-6s. The damper fault is well-estimated as plotted in Fig.7.5.2 Actuator faults estimation for a MIMO systemLet us consider the following MIMO system:

2 Fig. 8 .Fig. 9 .

 289 Fig.8. Disturbance, control inputs, varying parameters, switching signal Fig.9shows the estimation of the effectiveness factors λ 1 , λ 2 . Obviously, despite of the input disturbance w(t), the switched observer allows to have a good estimation of the coefficients λ 1 , λ 2 .

Table 1

 1 Parameters of SOBEN car model

	Parameters m s [kg]	m us [kg]	k s [N/m] k t [N/m]
	Value	2.28	0.25	1800	12269


Sébastien Marmin

David Ginsbourger

Efficient batch-sequential Bayesian optimization with moments of truncated Gaussian vectors

Keywords: Kriging, Expected Improvement, Parallel Optimization

We deal with the efficient parallelization of Bayesian global optimization algorithms, and more specifically of those based on the expected improvement criterion and its variants. A closed form formula relying on multivariate Gaussian cumulative distribution functions is established for a generalized version of the multipoint expected improvement criterion. In turn, the latter relies on intermediate results that could be of independent interest concerning moments of truncated Gaussian vectors. The obtained expansion of the criterion enables studying its differentiability with respect to point batches and calculating the corresponding gradient in closed form. Furthermore, we derive fast numerical approximations of this gradient and propose efficient batch optimization strategies. Numerical experiments illustrate that the proposed approaches enable computational savings of between one and two order of magnitudes, hence enabling derivativebased batch-sequential acquisition function maximization to become a practically implementable and efficient standard.

Introduction

Since their beginnings about half a century ago [START_REF] Kushner | A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise[END_REF][START_REF] Zhilinskas | Single-step bayesian search method for an extremum of functions of a single variable[END_REF][START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF], Bayesian optimization algorithms have been increasingly used for derivative-free global minimization of expensive to evaluate functions. Typically assuming a continuous objective function f : x ∈ D ⊂ R d -→ f (x) ∈ R, single-objective Bayesian optimization algorithms consist in sequentially evaluating f at promising points under the assumption that f is a sample realization (path or trajectory) of a random field (Y (x)) x∈D . Such algorithms are especially popular in the case where evaluating f (x) requires heavy high-fidelity numerical simulations (or computer experiments, see notably [START_REF] O'hagan | Curve fitting and optimal design for prediction[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Schonlau | Computer Experiments and global optimization[END_REF][START_REF] Jones | Efficient global optimization of expensive blackbox functions[END_REF]), where x stands for some design parameters to be optimized over. Such expensive simulations are classically encountered in the resolution of partial differential equations from physical sciences, engineering and beyond [START_REF] Forrester | Engineering design via surrogate modelling: a practical guide[END_REF]. In recent years, Bayesian optimization also has attracted a lot of interest from the machine learning community [START_REF] Lizotte | Practical Bayesian Optimization[END_REF][START_REF] Brochu | A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF][START_REF] Osborne | Bayesian Gaussian Processes for Sequential Prediction, Optimization and Quadrature[END_REF][START_REF] Srinivas | Gaussian process optimization in the bandit setting: No regret and experimental design[END_REF], be it to optimize simulation-based objective functions [START_REF] Lizotte | Automatic gait optimization with gaussian process regression[END_REF][START_REF] Tesch | Using response surfaces and expected improvement to optimize snake robot gait parameters[END_REF][START_REF] Romero | Navigating the protein fitness landscape with gaussian processes[END_REF] or even to estimate tuning parameters of machine learning algorithms themselves [START_REF] Bartz-Beielstein | Proc. of CEC-05[END_REF][START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF][START_REF] Snoek | Practical Bayesian optimization of machine learning algorithms[END_REF]. In both communities, a Gaussian random field (or Gaussian Process, GP) model is often used for Y , so that prior information on f is taken into account through a trend function m : D -→ R and a covariance kernel k : (x, x ′ ) : D × D -→ R. Once m and k are specified, possibly up to some parameters to be inferred based on data, the considered GP model can be used as an instrument to locate the next evaluation point(s) via so-called infill sampling criteria, also referred to as acquisition functions or simply as criteria. While a number of Bayesian Y ( x n + 1 )

T n
Improvement zone

(Y( x n + q )< T n ) Improv. zone (Y( x n + 1 )< T n ) Null improvement zone
Figure 1: Illustration of the principles underlying q-EI for d = 1, n = 4, q = 2. Left: Gaussian process prediction of a function f from observations A n (depicted by black crosses). The green horizontal line stands for T n , the smallest response value from A n . Three conditional simulation draws are plotted in orange and various point symbols represent their respective values at two unobserved locations x n+1 and x n+q . Right: distribution of the random vector (Y (x n+1 ), Y (x n+q )) ⊤ knowing A n (black contours). For each point symbol, the length of the purple segment represents the improvement realized by the corresponding sample path. The multipoint EI is the expectation of this length, or in other words, it is the integral of the improvement (grey-scale function) with respect to the conditional distribution of (Y (x n+1 ), Y (x n+q ))

⊤ knowing A n .

optimization criteria have been proposed in the literature (see, e.g., [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF][START_REF] Frazier | A knowledge-gradient policy for sequential information collection[END_REF][START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF][START_REF] Srinivas | Gaussian process optimization in the bandit setting: No regret and experimental design[END_REF][START_REF] Contal | Gaussian process optimization with mutual information[END_REF] and references therein), we concentrate here essentially on the Expected Improvement (EI) criterion [START_REF] Mockus | Bayesian Approach to Global Optimization[END_REF][START_REF] Jones | Efficient global optimization of expensive blackbox functions[END_REF] and on variations thereof, with a focus on its use in synchronous batch-sequential optimization. Denoting by x 1 , . . . , x n ∈ D points were f is assumed to have already been evaluated and by x n+1:n+q := (x n+1 , . . . , x n+q ) ∈ D q a batch of candidate points where to evaluate f next, the multipoint EI is defined as

EI n (x n+1:n+q ) = E n min i=1,...,n Y (x i ) - min j=n+1,...,n+q Y (x j ) + , (1)
where E n refers to the conditional expectation knowing the event A n := {Y (x 1 ) = f (x 1 ), . . . , Y (x n ) = f (x n )}. One way of calculating such criterion is to rely on Monte Carlo simulations. Figure 1 illustrates both what the criterion means and how to approach it by simulations, relying on three samples from the multivariate Gaussian distribution underlying Equation (1). Our main focus here, in contrast, is on deriving Equation (1) in closed form, studying the criterion's differentiability, and ultimately calculating and efficiently approximating its gradient in order to perform efficient batch optimization using derivative-based deterministic search. Now, for q = 1, it is well known that EI can be expressed in closed form as a function of m n (x) = E n (Y x ) and s n (x) = var n (Y (x)) as follows

EI n (x) = s n (x) (u n (x)Φ(u n (x)) + ϕ(u n (x))) if s n (x) = 0 and 0 else, (2)
where u n (x) = (min i=1,...,n f (x i )m n (x))/s n (x) (defined for s n (x) = 0) and Φ, ϕ are the cumulative distribution function and probability density function of the standard Gaussian distribution, respectively. When deriving Equation (2), Equation ( 1) happens (hence for q = 1) to involve a first order moment of the truncated univariate Gaussian distribution. As shown in [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF] and developed further here, it turns out that Equation (1) can be expanded in a similar way in the multipoint case (q ≥ 2) relying on moments of truncated Gaussian vectors. This is essential for the open challenges tackled here of efficiently calculating and optimizing the multipoint criterion of Equation [START_REF] Azimi | Batch bayesian optimization via simulation matching[END_REF].

The applied motivation for having batch-sequential EI algorithms is strong, as distributing evaluations of Bayesian optimization algorithms over several computing units allows significantly reducing wall-clock time and with the fast popularization of clouds, clusters and GPUs in recent years it is becoming always more commonplace to launch several calculations in parallel. Even at a slightly inflated price and scripting effort, reducing the total time off is often a primary goal in order to deliver conclusions involving heavy experiments, be they numerical or laboratory experiments, in studies subject to hard time limitations. Obviously, given its practical importance, the question of parallelizing EI algorithms and alike by selecting q > 1 points per iteration has been already tackled in a number of works from various disciplinary horizons (including notably [START_REF] Queipo | Assessing the value of another cycle in surrogate-based optimization[END_REF][START_REF] Azimi | Batch bayesian optimization via simulation matching[END_REF][START_REF] Desautels | Parallelizing exploration-exploitation trade-offs with Gaussian process bandit optimization[END_REF][START_REF] Contal | Parallel gaussian process optimization with upper confidence bound and pure exploration[END_REF][START_REF] González | Batch bayesian optimization via local penalization[END_REF]). Here we essentially focus on approaches relying on the maximization of Equation ( 1) and related multipoint criteria. The multipoint EI of Equation ( 1) has been defined in [START_REF] Mockus | Bayesian Approach to Global Optimization[END_REF][START_REF] Schonlau | Computer Experiments and global optimization[END_REF] and first calculated in closed form for the case q = 2 in [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF]. For the case q ≥ 3, a Monte Carlo scheme and some sub-optimal batch selection strategies were proposed. Further work on Monte Carlo simulations for multipoint EI estimation can be found in [START_REF] Janusevskis | Expected improvements for the asynchronous parallel global optimization of expensive functions : Potentials and challenges[END_REF][START_REF] Girdziusas | Parallel budgeted optimization applied to the design of an air duct[END_REF]; besides this, stochastic simulation ideas have been explored in [START_REF] Frazier | Parallel global optimization using an improved multi-points expected improvement criterion[END_REF] for maximizing this multipoint EI criterion via a stochastic gradient algorithm, an approach recently investigated in [START_REF] Wang | Parallel bayesian global optimization of expensive functions[END_REF]. Meanwhile, a closed-form formula for the multipoint EI relying on combinations of (q -1)and q-dimensional Gaussian cumulative distribution functions was obtained in [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF], a formula which applicability in reasonable time is however restricted to moderate q (say q ≤ 10) in the current situation. Building upon [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF], [START_REF] Marmin | Machine Learning, Optimization, and Big Data: First International Workshop, MOD 2015[END_REF] recently calculated the gradient of the multipoint EI criterion in closed form and obtained some first experimental results on (non-stochastic) gradient-based multipoint EI maximization.

Our aim in the present paper is to present a set of novel analytical and numerical results pertaining to the calculation, the computation, and the maximization of the multipoint EI criterion. As most of these novel results apply to a broader class of criteria, we first present in Section 2 a generalization of the multipoint EI that allows accounting for noise in conditioning observations and also exponentiating the improvement. This generalized criterion is calculated using moments of truncated Gaussian vectors in the flavour of [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF]. The obtained formula is then revisited in the standard case (noise-free with an exponent set to 1), leading to a numerical approximation of the multipoint EI with arbitrary precision and very significantly reduced computation time. Next, the (qd)-dimensional maximization of the multipoint EI criterion is discussed in Section 3, where the differentiability of the generalized criterion is studied, its analytical gradient is calculated, and further numerical approaches for fast gradient approximations with controllable accuracy are presented. Finally, Section 4 is dedicated to numerical experiments where, in particular, a multistart derivative-based multipoint EI maximization algorithm highlighting the benefits of the considered methodological principles and the proposed fast approximations is tested and compared to baseline strategies.

2 Criteria in parallel Bayesian optimization

General definition of Expected Improvement

Throughout this section the objective function f may be observed noise-free or in noise, meaning that at some arbitrary iteration i the observed value may be f (x i ) or f (x i ) + ε i where ε i is a realization of a zero mean Gaussian random variable with known (or estimated and plugged-in) variance. f is assumed to be one realization of a random field Y , where Y has a Gaussian random field (GRF) distribution conditionally to events of the form A n := {Y (x 1 ) = f (x 1 ), . . . , Y (x n ) = f (x 1 )} (with conditioning on Y (x i ) + ε i in the noisy case, see for instance [START_REF] Picheny | Quantile-based optimization of noisy computer experiments with tunable precision[END_REF]). This setup naturally includes the case where Y is a GRF, but also the so-called Universal Kriging settings where Y is the sum of a trend with an improper prior and a GRF [START_REF] O'hagan | Curve fitting and optimal design for prediction[END_REF][START_REF] Oakley | Bayesian inference for the uncertainty distribution of computer model outputs[END_REF]. Note that in noisy cases the ε i 's are generally assumed to be independent (although the case of ε i 's forming a Gaussian vector is tractable), but more essentially they are assumed independent of Y .

In batch-sequential Bayesian Optimization we are interested in computing sampling criteria J n depending on q ≥ 1 new points x n+1:n+q = (x n+1 , . . . , x n+q ) ∈ D q . At any step of corresponding synchronous parallel algorithms, the next batch of q points x ⋆ n+1:n+q is then defined by globally maximizing J n over all possible batches:

x ⋆ n+1:n+q ∈ arg max xn+1:n+q∈D q J n (x n+1:n+q ). (3)
Values of such criteria typically depend on x n+1:n+q through the conditional distribution Y (x 1:n+q )|A n , simplifying to Y (x n+1:n+q )|A n in the noiseless context. Conditional mean and covariance functions are analytically formulated via the so-called kriging equations, see e.g. [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by Kriging-Based Metamodelling and Optimization[END_REF]. Working out these criteria thus generally boils down to Gaussian vector calculus, which may become intricate and quite cumbersome to implement as q (or n + q, in noisy settings) increases. Our considered generalized version of the multipoint EI criterion, that allows accounting for a Gaussian noise in the conditioning observations and also for an exponentiation in the definition of the improvement, is defined as:

EI n (x n+1:n+q ) = E n min ℓ=1,...,n Y (x ℓ ) -min k=1,...,q Y (x n+k ) α + , (4)
where α ∈ N\{0}, E n (•) = E(•|A n ) and (•) + := max(0, •). This form gathers several sampling criteria notably including q-EI, both in noiseless and noisy settings, and also a multipoint version of the generalized EI of [START_REF] Schonlau | Computer Experiments and global optimization[END_REF]. In addition, the obtained results apply to batch-sequential versions of the Expected Quantile Improvement [START_REF] Picheny | Quantile-based optimization of noisy computer experiments with tunable precision[END_REF] (EQI) and variations thereof, by a simply change of process from Y to the quantile process. We will show in proposition 2 that such generalized multipoint EI criteria can be formulated as a sum of moments of truncated Gaussian vectors. In the next subsection, in order to get a closed form for the generalized EI we first define these moments and derive some first analytical formulas, that might also be of relevance in further contexts.

Preliminaries on moments of truncated Gaussian distribution

We fix α ∈ N\{0} and p = n + q in noisy settings or p = q in noiseless settings.

Definition 1. Let Z be a Gaussian vector with mean m ∈ R p and covariance matrix Σ ∈ S p ++ , where S p ++ is the cone of positive definite matrices of R p×p . For all positive integer k ≤ p, we define the function

M k,α on R p × S p ++ by M k,α : (m, Σ) → M k,α (m, Σ) = E n Z α k ½ {Z≤0} , (5)
where the inequality Z ≤ 0 is to be interpreted component-wise.

If Z is composed of values of a GRF at a batch of q locations x n+1:n+q , we use the notation M k,α (Z(x n+1:n+q )) := M k,α (m(x n+1:n+q ), Σ(x n+1:n+q )). We obtain the moments M k,α (m, Σ) of a truncated Gaussian distribution by an extension of Tallis' technique [START_REF] Tallis | The moment generating function of the truncated multi-normal distribution[END_REF] to any order, presented in the following proposition:

Proposition 1. The function G : R p × R p × S p ++ → R defined by G(t, m, Σ) = e 1 2 (t+Σ -1 m) ⊤ Σ(t+Σ -1 m)-m ⊤ Σ -1 m Φ p,Σ (-m -Σt) , (6)
where Φ p,Σ (•) is the cumulative distribution function of the centered p-variate normal distribution, is infinitely differentiable, and the moments M k,α are given by:

M k,α (m, Σ) = ∂ α G(•, m, Σ) ∂t α k t=0 . (7)
The proof of this Proposition is given in appendix B.1 and relies on calculating the moment generating function t → E exp t ⊤ Z ½ {Z≤0} . Even if an analytical formula can be obtained at any order of differentiation α, the complexity of derivatives in equation ( 7) increases rapidly. We give below the results for α equals 1 and 2.

Case α = 1

Differentiating G with respect to t yields:

∂G ∂t (t, m, Σ) = exp 1 2 t + Σ -1 m ⊤ Σ t + Σ -1 m -m ⊤ Σ -1 m × Σ t + Σ -1 m Φ p,Σ (-m -Σt) -Σ∇Φ p,Σ (-m -Σt)
where ∇Φ p,Σ is the gradient of Φ p,Σ (see appendix A.1 for an analytical derivation). Taking t = 0 in the previous equation gives

M k,1 (m, Σ) = m k Φ p,Σ (-m) -Σ ⊤ k ∇Φ p,Σ (-m) (8)
where Σ k is the k th column of Σ. It is shown in appendix A.1 that computing each of the p components of ∇Φ p,Σ requires to compute a multivariate CDF of the normal distribution in dimension p -1. The number of calls to this function for computing the first moment of the truncated Gaussian distribution is thus of O(p).

Case α = 2

Similarly, differentiating G twice with respect to t yields

M k,2 (m, Σ) = (Σ kk + m 2 k )Φ p,Σ (-m) + Σ ⊤ k ∇∇ ⊤ Φ p,Σ (-m)Σ k + 2m k M k,1 (m, Σ). ( 9
)
For readability, the detailed formula of ∇∇ ⊤ Φ p,Σ , the Hessian matrix of Φ p,Σ , is sent to Appendix A.2. The number of calls to the multivariate normal CDF is of O(p 2 ).

Analytic formulas for generalized q -EI

The previous results obtained for the moments of the truncated normal distribution turn out to be of interest for computing the generalized q-EI introduced in Equation ( 4), as shown by the following proposition.

Proposition 2. For x n+1:n+q ∈ D q , the criterion EI n defined by (4) exists for all α and can be written as a sum of moments of truncated normal distributions

EI n (x n+1:n+q ) = n ℓ=1 q k=1 M n+k-1,α Z (ℓ,k) (x n+1:n+q ) , (10)
with Z (ℓ,k) (x n+1:n+q ) a vector of size n + q -1 defined, noting

Y i := Y (x i ), by Z (ℓ,k) i =        Y ℓ -Y i if 1 ≤ i ≤ ℓ -1, Y ℓ -Y i+1 if ℓ ≤ i ≤ n -1, Y k -Y i+1 if n ≤ i ≤ n + q -1 and i = n + k -1, Y k -Y ℓ if i = n + k -1.
Moreover, in the noiseless case the random vector (Y (x 1 ), . . . , Y (x n )) becomes deterministic conditionally to A n . Denoting by ℓ 0 the (smallest) index of the minimal observation, i.e. Y ℓ0 = min ℓ=1,...,n Y ℓ , and writing Z (k) (x n+1:n+q ) the vector of the q last components of Z (ℓ0,k) (x n+1:n+q ), Equation (10) is simplified to:

EI n (x n+1:n+q ) = q k=1 M k,α Z (k) (x n+1:n+q ) . ( 11
)
Remark 1. In the rest of the article we also use the following compact notation for the (n + q -1)-dimensional vector Z (ℓ,k) (x n+1:n+q ):

Z (ℓ,k) (x n+1:n+q ) = A (ℓ,k) (Y (x 1 ), . . . , Y (x n+q )) ⊤ , (12)
where A (ℓ,k) is a matrix implicitly defined by

Z (l,k) i of proposition 2.
The proof of Proposition 2 is relegated to Appendix C for conciseness. Equation [START_REF] Cressie | The moment-generating function and negative integer moments[END_REF] highlights that the computation of the generalized q-EI in noisy settings is challenging since it involves computing nq different moments, each requiring (n + q) α calls to the multivariate normal CDF in a dimension close to n + q. Even for α = 1 and moderate q, the linear dependence in the number of observations n makes the use of this criterion challenging in application. Regarding the noiseless criterion, the computation of q moments is more affordable, at least for moderate q, but one has to keep in mind that the ultimate goal here is to perform global maximization of the considered criteria. It is thus important to bring further calculation speed-ups in order to perform this optimization in a reasonable time compared to the evaluation time of the objective function f , assumed expensive to evaluate. The next section discusses these matters and proposes faster formulas to compute both q-EI and its gradient.

3 Computing and optimizing the criteria

Generalities

Maximizing the EI n expressions given in Equation ( 10) (noisy settings) or (11) (noiseless settings) is difficult. These maximizations are performed with respect to a batch of q points x n+1:n+q ∈ (R d ) q , and are thus optimization problems in dimension dq. In this space, the objective function to be maximized, is not convex in general and has the interesting property that the q points in the batch can be permuted without changing the value of EI n ; i.e. EI n ((x n+1 , . . . , x n+q )) = EI n ((x n+σ(1) , . . . , x n+σ(q) )) for any permutation σ of {1, . . . , q}. With this property, one can reduce the measure of the search domain by q!, e.g. by imposing that the first coordinate of the q points in the batch are in ascending order. We will restrict our attention here to the use of multi-start gradient based local optimization algorithms acting on the whole input domain D q ⊂ R dq , that do not exploit the structure of the problem but do not seem to be affected by this, at least with the chosen settings regarding the starting designs. Our contribution here will be to propose a faster formula for computing the first moments M k,1 previously presented, as well as their derivatives. This will yield an easier computation of both the generalized EI and its dq-dimensional gradient. Besides, a second approximate but faster formula to further reduce the calculation time of the gradient will be introduced.

Gradient of the generalized q -EI

In this section, we extend the analytical gradient calculation of the q-EI performed in [START_REF] Marmin | Machine Learning, Optimization, and Big Data: First International Workshop, MOD 2015[END_REF] to the case of the generalized noisy and noise-free q-EI, and provide in turn a more concise formula. Again, the presented formulas rely on results on moments of truncated Gaussian distributions. Proposition 3. Let x n+1:n+q ∈ D q be a batch such that the conditional covariance matrix

(cov ( Y (x n+i ) , Y (x n+j )| A n )) 1≤i,j≤q is positive definite and the functions E ( Y (•)| A n ) and (cov ( Y (•), Y (x n+j )| A n )) j=1,.
..,q are differentiable at each point x n+i (1 ≤ i ≤ q). These derivatives are written m ′(i) ∈ R d and Σ ′(i) ∈ R q×d respectively. In this setup, the EI n function of Equation ( 10) is differentiable and its derivative with respect to the j th coordinate of the point x n+i is

∂EI ∂x ij (x n+1:n+q ) = n ℓ=1 q k=1 m ′(i) j A (l,k) i ⊤ ∂M n+k-1,1 ∂m Z (ℓ,k) + (13) tr A (l,k) Γ ′(i,j) A (l,k) ⊤ ∂M n+k-1,1 ∂Σ Z (ℓ,k) ,
where Γ ′(i,j) = Σ

′(i) u,j δ i,v + Σ ′(i) v,j δ i,u u,v
∈ R q×q , and δ is the Kronecker symbol. The derivatives

∂M n+k-1,1 ∂m
and

∂M n+k-1,1 ∂Σ are calculated in Appendix B.2.
This new expansion of the gradient of the generalized EI as a sum of derivatives of first order moments will prove to be very useful thanks to formulas presented next.

Fast numerical estimation of first order moments and their derivatives

Let us now focus on the practical implementation of the closed-form formula of Equation [START_REF] Cressie | The moment-generating function and negative integer moments[END_REF]. We take α = 1 and note p = n + q in noisy settings and p = q in noiseless settings. As mentioned before, the computation of the noisy or noiseless q-EI (see, Eqs. [START_REF] Cressie | The moment-generating function and negative integer moments[END_REF], [START_REF] Ginsbourger | DiceOptim: Kriging-Based Optimization for Computer Experiments[END_REF]) requires calls to the CDF of the p and (p -1)-variate normal distribution, Φ p and Φ p-1 . These CDFs are here computed using the Fortran algorithms of [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF] wrapped in the mnormt R package [START_REF] Azzalini | mnormt: The multivariate normal and t distributions[END_REF]. A quick look at Eqs. ( 8), [START_REF] Cressie | The moment-generating function and negative integer moments[END_REF] suggests that the noisy q-EI requires nq evaluations of Φ p and nq 2 evaluations of Φ p-1 . For the noiseless case (see, Equation [START_REF] Ginsbourger | DiceOptim: Kriging-Based Optimization for Computer Experiments[END_REF]), the number of calls are divided by n. In both cases, a slight improvement can be obtained by noticing a symmetry which reduces the number of Φ p-1 calls from nq 2 (resp. q 2 in the noiseless case) to nq(q + 1)/2 (resp. q(q + 1)/2). This symmetry is justified in Appendix E. Despite this improvement, and even in the classical noiseless case, the number of Φ p-1 calls is still proportional to q 2 . We now give new efficient and trustworthy expansion that enables a fast and reliable approximation of first order moments of truncated Gaussian vectors M k,1 by reducing this number of calls to O(q). Proposition 4. Let ε > 0, and let Z be a Gaussian random vector with mean vector and covariance matrix (m, Σ) ∈ R p × S p ++ . Then we have

M k,1 (m, Σ) = 1 ε (e m k ε Φ p,Σ (-εΣ k -m) -Φ p,Σ (-m)) + O(ε 2 ). ( 14
)
Proof. Let us consider the function

g k : t ∈ R → e m k t Φ p,Σ (-Σ k t -m).
This function g k is tangent at t = 0 with the function t ∈ R → G(te k ), where the function G is introduced in Proposition 1 and e k is the k th vector of the canonical basis. It follows from Proposition 1 that

M k,1 (m, Σ) = ∂G ∂t k (t, m, Σ) t=0 ,
and we obtain the announced result by Taylor expansion of g k .

The obtained formula simply uses the approximation of a moment with finite differences of the moment generating function. We showed here that instead of fully computing the moment generating function, we can expand the simpler tangent function g k . For conciseness, we name here the use of this formula "tangent moment method". This formula thus enables approximating the first order moment M k,1 at the cost of only two calls to Φ p . Hence, from Equation [START_REF] Ginsbourger | DiceOptim: Kriging-Based Optimization for Computer Experiments[END_REF], computing a noiseless q-EI can be performed at the cost of 2q calls to Φ q . Besides, a similar approach can be applied to approximate the gradient of q-EI through faster computations of

∂M k,1 ∂m =Φ p,Σ (-m)e k - 1 ε (e m k ε ∇Φ p,Σ (-Σ k ε -m) -∇Φ p,Σ (-m)) + O(ε 2 ) ( 15
)
∂M k,1 ∂Σ = - ∂Φ p,Σ ∂x v (-m) δ u,k + ∂Φ p,Σ ∂x u (-m) δ v,k u,v≤p (16)
+ 1 ε e m k ε ∇∇ ⊤ Φ p,Σ (-Σ k ε -m) -∇∇ ⊤ Φ p,Σ (-m) + O(ε 2 )
where ∇∇ ⊤ Φ p,Σ is the Hessian matrix of Φ p,Σ (see Appendix A.2 for details).

As before, these formulas enable reducing the number of calls to the multivariate CDF by an order q. For the computation of q-EI this number goes from O(q 2 ) to O(q). For computing its dq-dimensional gradient, it goes from O(q 4 ) to O(q 3 ). The latter complexity suggests restricting to moderate values of q in applications. In the next section we present further results that enable further reducing the complexity for numerically estimating the gradient.

A slightly biased but fast proxy of the gradient

The key idea to obtain further computational savings is summarized in this section. We first strategically decompose the gradient of moments as a sum of two terms. Proposition 6. Let us consider a Gaussian multivariate random field Z from R d to R p . For x ∈ R d , let us denote by m(x) and Σ(x) the mean and the covariance matrix of Z(x). Let x a ∈ R d and assume that Σ(x a ) is positive definite. Also, assume that the functions x → m(x), x → Σ(x) and x → (cov(Z i (x), Z j (x a ))) i,j≤p are differentiable at x = x a . Then the following decomposition holds for k = 1, . . . , p.

∇ x [M k,α (m(x), Σ(x))]| x=xa := ∇ x E Z α k (x)½ {Z(x)≤0} x=xa = ∇ x E Z α k (x)½ {Z(xa)≤0} x=xa + ∇ x E Z α k (x a )½ {Z(x)≤0} x=xa . (17)
Proof. Σ(•) is continuous at x a , so there exists a neightborhood V x a of x a such that for all x ∈ V x a , Σ(x) is positive definite. Let us define on V x a × V x a :

g(u, v) = E Z α k (u)½ {Z(v)≤0} .
Applying equation ( 22) of appendix B.1, for all u and v, g(u, v) is a moment generated by differentiation of the following function:

M u,v : t → e 1 2 (Σkk(u)t 2 +2tm k (u)) Φ p,Σ(v) -m(v) -t (cov(Z k (u), Z j (v))) ⊤ j≤p . (18)
The analytical form of equation ( 18) and the assumed differentiability at x a ensure existence of partial derivatives of g = (u, v) →

d α Mu,v dt α (0) at (x a , x a ). So to conclude, ∇ x [M k,α (m(x), Σ(x))]| x=xa = ∇ x [g(x, x)]| x=xa = ∂ ∂u [g(u, x a )] u=xa + ∂ ∂v [g(x a , v)] v=xa
The latter decomposition can be interpreted as follows: infinitesimal variations of (m(x), Σ(x)) around (m(x a ), Σ(x a )) modify the moments M k,α (m(x), Σ(x)) in two ways. First, it modifies the distribution of Z α k (x), second it changes the distribution of the truncation ½ {Z(x)≤0} . For the particular case of q-EI, we propose to neglect this second variation. Applying this approximation to [START_REF] Ginsbourger | DiceOptim: Kriging-Based Optimization for Computer Experiments[END_REF] gives for

X 0 ∈ D q , ∇ x n+j EI(x n+1:n+q ) xn+1:n+q=X0 = q k=1 ∇ xn+j E (T -Y (x n+k )) α ½ {A (k) Y (xn+1:n+q)≤0} x n+1:n+q =X0 ≈ q k=1 ∇ x n+j E (T -Y (x n+k )) α ½ {A (k) Y (X0)≤0} xn+1:n+q=X0 = -∇ x n+j E Y (x n+j ) α ½ {A (j) Y (X0)≤0} xn+1:n+q=X0 = -E ∇ xn+j Y (x n+j ) α | x n+1:n+q =X0 ½ {A (j) Y (X0)≤0} , (19)
where the last step is obtained by mean square differentiability of the process x → Y (x) α ½ {B} , with B an event constant with respect to x, see Appendix D. We can observe that this approximation makes a summation term disappear. The computation of this formula requires (d + 1) evaluations of q-variate Gaussian CDF. Indeed, Equation [START_REF] González | Batch bayesian optimization via local penalization[END_REF] indicates that each component of the gradient vector can be considered as a moment of a truncated Gaussian vector, so we can apply the results of section 2. In particular, when α = 1, applying Proposition 4, two Gaussian CDF calls are needed for each of the d components, leading to 2d evaluations. Besides, from Equation ( 14), the second CDF call does not depend on k, which implies that this term is common for every dimension. Thus the gradient of Equation ( 19) finally comes with d+1 CDF evaluations instead of 2d. For a full gradient with respect to all q points of the batch, we then need q(d+1) CDF evaluations -a substantial improvement compared to the O(q 4 ) obtained in [START_REF] Marmin | Machine Learning, Optimization, and Big Data: First International Workshop, MOD 2015[END_REF] and the O(q 3 ) obtained in the previous section. The complexities for computing moments, q-EI and its gradients, expressed in terms of number of calls to the Φ function, are summarized in Table 1. These new computational savings come at the price of a non-exact gradient calculation. A first numerical validation is represented in Figure 2. On this example, we observe small (1 × 10 -2 ) relative errors between the exact and approximate gradient of dimension q × d = 4 (the biggest difference vector has a norm of 0.13, compared to an exact gradient norm of 13.1). We also observe that the relative error appears to be typically smaller with higher q-EI, which is promising for q-EI maximisations. However, this apparently trustful but non-exact calculation naturally raises the question of the impact of such an approximation on the performances of gradient-based q-EI maximization algorithms. As we will see in the next section, this proxy gradient turned out to enable quite competitive q-EI maximization performances based on numerical experiments.

Application

The goal of this section is to illustrate the usability of the proposed gradient-based q-EI maximization schemes and in particular the improvements brought by the fast formulas detailed in the previous 19), with α = 1, q = 2, d = 2.

From left to right: 1) Norm of the q-EI gradient, with respect to the first batch point (the other point is fixed in the center of [0, 1] d ) ; 2) Norm of the difference vector between the analytical gradient and its approximation ; 3) Relative error (norm of the difference divided by the real norm) computed on 3000 random batches sampled uniformly in [0, 1] d×q , with respect to their q-EI.

Number of CDF evaluations Φ

q-3 Φ q-2 Φ q-1 Φ q Total M k,1 analytic q 1 O(q) tangent moment 2 2 EI analytic q+1 2 q O(q 2 ) tangent moment 2q O(q) ∇M k,1 analytic 3 q 3 3 q 2 2q 1 O(q 3 ) tangent moment 2 q 2 2q 2 O(q 2 ) proxy d + 1 O(d) ∇EI analytic 6 q+1 4 3 q+1 3 (3q 2 + q)/2 q O(q 4 ) tangent moment q 2 (q -1) 2q 2 2q O(q 3 ) proxy q(d + 1) O(qd)
Table 1: In noiseless settings, total number of calls to the CDF of the multivariate Gaussian distribution for computing M k,1 , q-EI, their gradients and their approximations, depending on q and d. For q-EI in noisy setting, replace q by p = n + q and multiply each number of calls by n. Figure 3: Computation times for q-EI or its gradient as a function of the batch size q (logarithmic scale). We take an averaged computation time over 1000 batches (except for points marked with a * , averaged over 150 batches).

sections. The relevance of using sequential sampling strategies based on the q-EI maximization has already been investigated before (see, [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF][START_REF] Wang | Parallel bayesian global optimization of expensive functions[END_REF][START_REF] Marmin | Machine Learning, Optimization, and Big Data: First International Workshop, MOD 2015[END_REF]) and all these articles pointed out the importance of calculation speed which often limits the use of q-EI based to moderate q. We do not aim again at proving the performance of q-EI based sequential strategies. Instead we aim at illustrating the gain, in computation time, brought by the fast formulas and show that using the approximate gradient obtained in Equation [START_REF] González | Batch bayesian optimization via local penalization[END_REF] does not impair the ability to find batches with (close to) maximal q-EI.

Objective function and pure calculation speed

The objective function is the so-called Borehole function [START_REF] Harper | Sensitivity/uncertainty analysis of a borehole scenario comparing latin hypercube sampling and deterministic sensitivity approaches[END_REF]. It has been previously used for testing methods using a surrogate model [START_REF] Worley | Deterministic uncertainty analysis[END_REF][START_REF] Gramacy | Gaussian process single-index models as emulators for computer experiments[END_REF]

φ = 2πT u (H u -H l ) ln r rw 1 + 2LTu (ln( r rw )r 2 w Kw) + Tu T l . ( 20
)
Here, the objective function f is obtained by rescaling φ on the input domain D = [0, 1] 8 . An analytical study of variations shows that there is a unique global minimum at x * = (0, 1, 0, 0, 0, 1, 1, 0) ⊤ , with f (x * ) ≈ 1.1918. Before using sequential strategies to minimize f , we look at empirical computation times for evaluating q-EI and its gradient as a function of the batch size q. For the computations, the so-called "analytic" method relies on the state of the art formulas of [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF][START_REF] Marmin | Machine Learning, Optimization, and Big Data: First International Workshop, MOD 2015[END_REF] with a number of calls to the multivariate normal cdf of respectively O(q 2 ) and O(q 4 ). The "tangent moment" method uses our formula for moment calculation to yield q-EI and its gradient (see, Equations ( 14) and ( 15),( 16)). Finally, for computing the gradient only, the "proxy" method relies on Equation [START_REF] González | Batch bayesian optimization via local penalization[END_REF].

Figure 3 exhibits computation times averaged over 1000 batches drawn uniformly. The Gaussian process model is based on an initial design of n 0 = 10d = 80 points drawn from a optimum-LHS procedure [START_REF] Kenny | Algorithmic construction of optimal symmetric latin hypercube designs[END_REF]. We use the Matérn (ν = 3/2) tensor product covariance function and estimate the hyperparameters by maximum likelihood using the DiceKriging R package [START_REF] Ginsbourger | DiceOptim: Kriging-Based Optimization for Computer Experiments[END_REF]. Figure 3 shows significant computational savings. For instance with q = 8, one gradient computation takes respectively 0.04s, 0.33s and 1.33s using respectively the proxy, the tangent moment and analytic methods. Since the complexity for computing a gradient with the proxy is of O(qd) against O(q 3 ) and O(q 4 ) for the two other methods, the savings of the proxy tend to increase with q. It should also be noted that these savings will be larger with decreasing domain dimension d. If we look at q-EI computations, the tangent moment method is 3.3 times faster than the analytic one when q = 8 and 6.5 times faster when q = 20; thanks to an O(q) complexity against O(q 2 ).

Experimental setup: sequential minimization strategies

We now perform a total of 50 minimizations of f , each using an initial design of experiments of n 0 = 80 points drawn from an optimum-LHS procedure with a different seed. Three different batch-sequential strategies are investigated.

The first one -serving as a benchmark -is a variation of the "Constant Liar Mix" heuristic [START_REF] Chevalier | Learning and Intelligent Optimization -7th International Conference[END_REF][START_REF] Wang | Parallel bayesian global optimization of expensive functions[END_REF] where, at each iteration, the batch of size q is chosen among several batches obtained from the Constant Liar heuristic [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF] with different lie levels. We use 7 lie levels fixed to the current maximum observation the current minimum observation, and the 2.5%, 10%, 50%, 90%, 97.5% quantiles of the conditional distribution of the point selected in the batch. A total of 7 batches are proposed at each iteration and the CL-mix heuristic picks the one with maximum q-EI.

The two other strategies considered here rely on pure q-EI maximization using a multistart BFGS algorithm with a stopping criterion of precision 2.2×10 -7 (parameter control$factr of the R function optim [START_REF]R: A Language and Environment for Statistical Computing[END_REF]). The gradients involved in the optimization are computed either with the tangent moment formula or the proxy. For the gradient-based q-EI maximization, we use a total of 10 starting batches obtained, again, using a Constant Liar heuristic with random lies sampled from the conditional distribution at the selected point. Finally we use two different batch sizes. When q = 8 we run a total of 10 iterations and when q = 4 we run 20 iterations. The hyperparameters of the GP model are re-estimated at each iteration after having incorporated the new observations.

First q -EI maximization

We first compare the performances, in terms of q-EI, of the multistart BFGS algorithm when the proxy gradient and the tangent moment methods are used. Table 2 compares the results at iteration 1 for these two methods and the CL-mix strategy. The results are averaged over the 50 initial designs. q = 4 q = 8 tangent moment 12. [START_REF] Tesch | Using response surfaces and expected improvement to optimize snake robot gait parameters[END_REF] Table 2: Average q-EI value of the optimal batches found for each of the 50 initial designs. The numbers between brackets are the average computation times.

As expected, the CL-mix heuristic yields batches with lower q-EI than the strategies directly maximizing q-EI. Also, for both q = 4 and q = 8, the two q-EI based methods have the same performance, which stresses out the relevance of the proxy method since the latter is about 1.6 times faster when q = 4 and 5.5 faster when q = 8.

Several q -EI maximization steps

We now compare the performances of the different q-EI maximization approaches after multiple batch evaluations. Figure 4 displays the average regret as a function of the iteration number (first row) and the total computation time (i.e. the time to evaluate f and find the next batch to evaluate) assuming respectively that the computation time of f is 0 seconds (i.e. instantaneous), two minutes and one hour (rows 2, 3, 4 respectively). Looking at the performances as a function of the iteration number (first row on Figure 4), the CL-mix heuristic, which samples a batch with lower q-EI at each step, leads in average to a slower convergence than the two other methods, for both q = 4 and q = 8. In contrast, the two strategies based on q-EI maximization have similar performances.

However, these conclusions do not hold when the regret is plotted as a function of the total computation time (rows 2, 3, 4 on Figure 4). First, when the computation time t eval of f is null (row 2) it is clear that q-EI-based sequential strategies are not adapted since they are too expensive. In this case, the CL-mix heuristic performs better and some other optimization strategies which are not metamodel-based would probably be more relevant. Second, when f is moderately expensive (i.e. t eval = 2 minutes), the proxy method and CL-mix have comparable performances when q = 8, but the proxy outperforms when q = 4. Besides, the proxy shows a much faster convergence than the tangent moment method when q = 8. The use of q-EI based strategies thus becomes relevant when t eval is larger than a few minutes, if the proxy is used. Finally, when t eval is equal to one hour, the use of q-EI based strategies is particularly recommended. In that case the relative improvement of the proxy compared to the tangent moment method tends to naturally vanish because of the long computation time of f . When f extremely expensive to compute, using the proxy is thus not essential. However, since it does not impair the ability to find a batch with large q-EI we still recommend to use it, especially when q is large.

Conclusion

In this article we provide a closed-form expression of generalized q-points Expected Improvement criterion for batch-sequential Bayesian global optimization. An interpretation based on moments of truncated Gaussian vectors yields fast q-EI formulas with arbitrary precision. Furthermore an new approximation for the gradient is shown to be even faster while preserving ability to find batches close to maximal q-EI. As the use of these strategies was previously considered cumbersome from a dozen of batch points, these formulas happen to be of particular interest to run q-EI based batch-sequential strategies for larger batch sizes. We show that these methods are implementable and efficient on a classic 8-dimensional test case. Additionally, some of the intermediate results established here might be of interest for other research questions involving moments of truncated Gaussian vectors and their gradients. Perspectives include deriving second order derivatives of q-EI and fast numerical estimates thereof. Also, we aim at improving the sampling of initial batches in multistart derivative-based q-EI maximization.

Acknowledgements: Part of this work has been conducted within the frame of the ReDice consortium, gathering industrial (CEA, EDF, IFPEN, IRSN, Renault) and academic ( École des Mines de Saint-Étienne, INRIA, Universität Bern) partners around advanced methods for computer experiments.

A Differentiating multivariate Gaussian CDF

We consider the CDF dimension p ≥ 2. We use the convention Φ 0 = 1. Figure 4: Log-scaled average regret of the three considered optimization strategies as a function of the iteration number (row 1) and the total computation time (rows 2, 3, 4) assuming that the computation times of f , t eval , are respectively 0 seconds, 2 minutes and 1 hour. Experiments are performed with q = 4 (left column) and q = 8 (right column).

A.1 Gradient

Using the following identity, derived from conditional distributions of a Gaussian vector, ∀i = 1, . . . , p, ϕ p,Σ (

x) = ϕ 1,Σii (x i ) ϕ p-1,Σ | i x -i -m |i,xi , with m |i,u = u Σii Σ -i,i and Σ |i = Σ -i,-i -1 Σii Σ -i,i Σ ⊤ -i,i
, we reformulate the integral of the Gaussian CDF:

∀i = 1, . . . , p, Φ p,Σ (x) = xi -∞ ϕ 1,Σii (u i ) Φ p-1,Σ |i x -i -m |i,ui du i .
Here indexed minus symbols, e.g. in Σ -i,i , refer to exclusions of a line or a column.

Finally we have

∇Φ p,Σ (x) = ϕ 1,Σii (x i ) Φ p-1,Σ|i x -i -m |i,xi i=1,...,p . (21)

A.2 Hessian

As for the computation of the gradient, we write ∀i, j = 1, . . . , p, i = j,

Φ p,Σ (x) = xi -∞ xj -∞ ϕ 2,Σij,ij ui u j Φ p-2,Σ |ij x -{i,j} -m |(i,j),(ui,uj ) du j du i , with m |(i,j),(u,u ′ ) = Σ -{ij},ij Σ -1 ij,ij u u ′ and Σ |ij = Σ -{ij},-{ij} -Σ -{ij},ij Σ -1 ij,ij Σ ⊤ -{ij},ij . So ∀i, j = 1, . . . , p, i = j, ∂ 2 Φ q ∂x i ∂x j (x) = ϕ 2,Σij,ij xi x j Φ p-2,Σ |ij (x -{i,j} -m |(i,j),(xi,xj) ).
When i = j, the differentiation of equation ( 21) gives,

∂ 2 Φ q ∂x 2 i (x) = - 1 Σ ii     x i ∂Φ p,Σ ∂x i (x) + p j=1 j =i Σ ij ∂ 2 Φ p,Σ ∂x i ∂x j (x)     .

B Moments of truncated multivariate Gaussian distribution B.1 Analytical formula (propositions 1 and 6)

We see here why we can derive an analytical formula of M k,α (m, Σ), with k ≤ s ∈ N\{0}, m ∈ R s and Σ ∈ S s ++ , by differentiating G, defined in equation [START_REF] Brochu | A tutorial on bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[END_REF]. It is known, see e.g. [START_REF] Cressie | The moment-generating function and negative integer moments[END_REF], that moments can be obtained differentiating the moment generating function G m,Σ,s :

M k,α (m, Σ) = ∂ α G m,Σ,s ∂t α k (0),
with, for r ∈ {1, . . . , s}, G m,Σ,r : t → E exp t ⊤ Z ½ {(Z1,...,Zr ) ⊤ ≤0} , Z ∼ N (m, Σ). We derive now an analytical formula for G m,Σ,r . As needed in Proposition 6, we derive an analytical formula for any r, and not only for r = s.

∀t ∈ R s , G m,Σ,r (t) = r times 0 -∞ ... 0 -∞ ∞ -∞ ... ∞ -∞ exp (t ⊤ z)ϕ Σ (z -m) dz 1 ... dz s = ϕ Σ (0) 0 -∞ ... 0 -∞ ∞ -∞ ... ∞ -∞ exp - 1 2 (z -m) ⊤ Σ -1 (z -m) -2t ⊤ z dz = e -1 2 -(t+Σ -1 m) ⊤ Σ(t+Σ -1 m)+m ⊤ Σ -1 m ϕ Σ (0) 0 -∞ ... 0 -∞ ∞ -∞ ... ∞ -∞ e -1 2 ((z-m-Σt) ⊤ Σ -1 (z-m-Σt) ⊤ ) dz = e 1 2 (t+Σ -1 m) ⊤ Σ(t+Σ -1 m)-m ⊤ Σ -1 m Φ r,(Σij ) i,j≤r -m -(Σ ij ) i≤r,j≤s t . (22)
In the frame of the proof of Proposition 6,

• if Σ k (u, v), the covariance matrix of (Z(v) ⊤ , Z k (u)) ⊤ , is positive definite, we take M u,v = t → G (m(v),m k (u)),Σ k (u,v),p ((0, . . . , 0, t) ⊤ ),
• else, as Σ(v) is definite positive, there exists only one index k 0 such as Z k (u) = Z k0 (v) almost surely (for example k = k 0 when u = v), and we have M u,v = t → G (m(v)),Σ(v),p ((0, . . . , t, . . . , ↑ k th 0 position 0) ⊤ ).

In both cases, equation ( 22) leads to equation [START_REF] Girdziusas | Parallel budgeted optimization applied to the design of an air duct[END_REF].

B.2 Differentiation with respect to mean and covariance

We differentiate here the equation ( 8) with respect to m and Σ.

With respect to the mean m +Σ ki ϕ Σii (-m i ) ∂ ∂Σ Φ p-1,Σ|i (-m |i ). [START_REF] Jones | Efficient global optimization of expensive blackbox functions[END_REF] with m| i = m -i -mi Σii Σ -i,i and Σ| i = Σ -i,-i -1 Σii Σ -i,i Σ ⊤ -i,i . Writting d Σ [m| i ] d Σ [Σ| i ] the differential of the functions Σ → m| i and Σ → Σ| i , we have: with:

d Σ [m| i ] (H) = m i Σ ii H -i,i (25)
d Σ [Σ| i ] (H) = H -i,-i + H ii Σ 2 ii Σ -i,i Σ ⊤ -i,i - 2 Σ ii H -i,i Σ ⊤ -i,i (26)
• E (r,s) = (δ ij ) i,j=1,...,p ,

• ∂ ∂Γ Φ p-1,Σ |i (-m |i ) the derivative of Γ → Φ p-1,Γ (-m |i ) evaluated at Σ |i . We use the Plackett's differential equation, extended by [START_REF] Berman | An extension of plackett's differential equation for the multivariate normal density[END_REF], to find

∂ ∂Γ Φ p-1,Σ |i (-m |i ) = ∇∇ ⊤ Φ p-1,Σ |i (-m |i ),
∇∇ ⊤ Φ is given in appendix A.2.

C Generalized q -EI as a sum of moments E n (Y (x ℓ ) -Y (x n+k )) α ½ {E ℓ,k } .

For each term (ℓ, k) of the sum, the conditioning event can be rewritten E ℓ,k = {Z (ℓ,k) (x n+1:n+q ) ≤ 0}, with Z (ℓ,k) a random vector of size n + q -1, defined by the following linear transformation of Y = (Y (x 1 ), . . . , Y (x n+q )) ⊤ : ∀i = 1, . . . , n + q -1, Z

(ℓ,k) i =        Y ℓ -Y i if 1 ≤ i ≤ ℓ -1 Y ℓ -Y i+1 if ℓ ≤ i ≤ n -1 Y n+k -Y i+1 if n ≤ i ≤ n + q -1, i = n + k -1 Y n+k -Y ℓ if i = n + k -1
Indeed, the first n -1 components of Z (l,k) ≤ 0 reflect {∀i ≤ n, i = ℓ; Y (x ℓ ) ≤ Y (x i )}, and the last components reflect {∀j ≤ q, j = k; Y (x n+k ) ≤ Y (x n+j )} and {Y (x n+k ) ≤ Y (x ℓ )}.

D Mean square differentiability of Y (x) α ½ {B} Let B be an event, Y be a mean-squared differentiable Gaussian process and α ∈ N. Then we have:

E Y (x + h) α -Y (x) α h ½ {B} - dY α dx (x)½ {B} 2 ≤ E Y (x + h) α -Y (x) α h - dY α dx (x)

E Symmetry argument

The term q(q+1) 2 comes from a symmetry occurring when summing terms with different index but actually equal. At fixed summation index ℓ in [START_REF] Cressie | The moment-generating function and negative integer moments[END_REF], we denote ω ki the i th term in the scalar product in [START_REF] Contal | Parallel gaussian process optimization with upper confidence bound and pure exploration[END_REF] for each M m+k-1,1 required for q-EI: ∀i, k = 1, . . . , q, ω ki = Σ (ℓ,k) ki ∇Φ p,Σ (ℓ,k) (-m (ℓ,k) ) i .

Then the following symmetry between indices i and k occurs:

∀i, k = 1, . . . , q, ω ki Σ (ℓ,k) ki ϕ 1,Σ (ℓ,k) ii (-m (ℓ,k) i ) = ω ik Σ (ℓ,i) ik ϕ 1,Σ (ℓ,i) kk (-m (ℓ,i) k
))

Indeed, using the formula of the derivative of CDF, (appendix A.1), leads to:

Proposition 5 .

 5 The following equations hold:

Figure 2 :

 2 Figure 2: Numerical validation of the approximation from equation (19), with α = 1, q = 2, d = 2. From left to right: 1) Norm of the q-EI gradient, with respect to the first batch point (the other point is fixed in the center of [0, 1] d ) ; 2) Norm of the difference vector between the analytical gradient and its approximation ; 3) Relative error (norm of the difference divided by the real norm) computed on 3000 random batches sampled uniformly in [0, 1] d×q , with respect to their q-EI.

  0

∂M k, 1 23 )ϕ

 123 ∂m (m, Σ) = Φ p,Σ (-m)e km k ∇Φ p,Σ (-m) + ∇∇ ⊤ Φ p,Σ (-m)Σ k . (With respect to the covariance Σ ∂M k,1 ∂Σ (m, Σ) = m k ∂ ∂Σ Φ p,Σ (-m) -p i=1 Σii (-m i ) Φ p-1,Σ|i (-m |i )E (k,i) +Σ ki ∂ ∂Σ ii ϕ Σii (-m i ) Φ p-1,Σ|i (-m |i )E (i,i)

  [m| i ] (E (r,s) ).∇Φ p-1,Σ|i (-m| i ) +tr ∂ ∂Γ Φ p-1,Σ|i -m |i .d Σ Σ |i E (r,s) E (r,s)

Proof.

  For given (ℓ, k) in {1, . . . , n} × {1, . . . , q}, we consider E ℓ,k the event that the random variable inside the expectation term of equation (4) equals (Y (x ℓ ) -Y (x n+k )) α . We haveE ℓ,k = {Y (x n+k ) ≤ Y (x ℓ )} ∩ {∀i ≤ n, i = ℓ; Y (x ℓ ) ≤ Y (x i )} ∩ {∀j ≤ q, j = k; Y (x n+k ) ≤ Y (x n+j )}Considering all pairs (ℓ, k), we have:EI n (x n+1:n+q ) =

  squared differentiability of Y α .

Y

  (x n+k ) ≤ Y (x ℓ ), Y (x n+j ) ≤ Y (x n+k ), ∀j = 1 . . . q, j = k, j = i Y (x n+i ) = Y (x n+k ), which is clearly symmetrical between i and k.


I Shili-Touzi

S De Tourdonnet

M Launay

T Doré

Objectives

To study four management scenarios of the system using STICS-crop model, to analyze the impact of the timing of the fescue cycle on the agronomic and environmental performance of the system. To find the best compromise between competitive effects and facilitative effect in intercropping winter wheat and red-fescue as cover crop (biomass accumulation and radiation interception).

Results

Conclusions

Advancing the sowing date of the fescue increases both competition effects (decreasing wheat biomass and yield) and facilitation effects (increasing total biomass and soil cover, decreasing the amount of solar radiation reaching the soil). It also increases the efficiency of radiation and nitrogen use. In most climatic scenarios the emergence of the fescue crop during the late tillering phase of the wheat crop gave the best results.

Materials and methods

Model used: STICS-CROP MODEL (Brisson et al., 2009) STICS was parameterized and evaluated using data obtained in two field experiments at the INRA Grignon experimental station (Picard et al., 2009) involving winter wheat and red-fescue grown alone and winter wheat intercropped with red fescue. The results of model evaluation are detailed in Shili-Touzi et al. (2010).

The scenarios simulated:

The wheat/red-fescue model was used to compare four management scenarios:

Scenario 0 : wheat grown as sole crop (reference scenario).

Scenario 1 : fescue emergence one week after the harvest of the wheat crop, corresponding to a double-cropping system.

Scenario 2 : fescue emergence in the spring (18 March), corresponding to relay intercropping:, wheat grown as a sole crop from October to March and then intercropped with fescue.

Scenario 3 : simultaneous emergence of wheat and fescue (10 October), corresponding to full intercropping.

For scenarios 2 and 3, fescue remainis alive after wheat harvest. These simulations were run over 35 years of climatic data from Versailles (48°48'N, 2°04'E), to enable us to take climate variability into account in the scenarios assessment.

Cumulative probability

Wheat + fescue dry matter (t.ha -1 )

EXPLORING MANAGEMENT SCENARIOS FOR INTERCROPPING WINTER WHEAT AND RED FESCUE AS COVER CROP USING STICS-MODEL Introduction

The use of a living cover crops in cropping system may help to preserve biodiversity, increase soil organic matter content, carbon sequestration and provide other ecosystem services such as natural pest regulation or nutrient recycling by maintaining actively growing vegetation during the whole year. But competition for light, water or nutrients between cover crop and commercial crop may decrease cash crop yield and returns Very few studies have been carried out to analyse the impact of various ways for managing this intercropping system.

Simulation results showed that wheat yield was not overly affected by fescue development in the intercropping system (Sc 2 and 3) in comparison with sole wheat yield (Sc 0 and 1).

Yield variability over time was similar for all four scenarios. the system was highly sensitive to the position of the fescue growth cycle, particularly in terms of dry matter production : overall dry matter production (wheat plus fescue) was greater than for sole-crop by 1.26, 1.84 and 2.19 t ha -1 on average for scenarios 1, 2 and 3, respectively. The difference between the scenarios is linked to differences in the timing of fescue dry matter production.

The simulated PAR reaching ground level over the entire simulated period showed that intercropping system ( Sc 2 and 3) decreased radiation transmission in comparison with sole wheat ( Sc 0 and 1) It has a direct effect on the reduction of radiation available for weed growth corresponding of 7%, 18% and 31%for scenarios 1, 2 and 3 , respectively.

Over the 35-year period, an increase in N acquisition by both crops was possible only in scenarios 1 and 2 , for which intercropping increased total N acquisition.

At the end of the simulated period (31 December), regardless of the scenario considered, intercropping with fescue resulted in lower soil inorganic N content than leaving the soil bare after the wheat harvest .

  , Launay M., Mary B., Beaudoin N., 2009. Conceptual basis, formalisations and parameterization of the STICS crop model. Quae, Paris, 304 p. Picard D., Ghiloufi M., Saulas P., de Tourdonnet S 2009. Does undersowing winter wheat with a cover crop increase competition for resources and is it compatible with high yield? Field Crops Res. Shili-Touzi I., De Tourdonnet S., Launay M., Doré T. 2010. Does intercropping winter wheat (Triticum aestivum) with red fescue (Festuca rubra) as a cover crop improve agronomic and environmental performance? A modeling approach. . Field Crops Res.

0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 5 6 7 8 9 10 11 Cumulative probability Wheat yield (t.ha -1 ) 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% 100% 15 17 19 21 23 25