Sandrine Bernard

François Chaplais
email: francois.chaplais@mines-paristech.fr

Fabienne Chaplais
email: fabienne@mathedu.fr

François Chaplais

Averaging on simple windows in deterministic optimal control

A windowed averaged scheme is defined for general control systems. The same method is used to average costs in optimal control problems (OCPs). A numerical parameter α can be computed, which expresses the distance between the original system and the averaged system in a weak sense.

Then, if we use the optimal control of the averaged OCP in the original OCP, the suboptimality of the control is bounded by an expression of the form Cα 2 .

Introduction

Historically, the method of averaging was introduced to study the motion of celestial bodies by solving a simple two body equation which is perturbed by the influence of other bodies (see the section in [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF] about the history of averaging). As developed in [START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste[END_REF], the framework was that of the perturbation of an orbit by the small influence of a periodic input. Averaging was then generalized in a geometric framework, notably in [START_REF] Arnold | Mathematical Methods of Classical Mechanics[END_REF][START_REF] Arnold | Geometrical Methods in the Theory o f Ordinary Differential Equations[END_REF]. A comprehensive book on the subject is [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF].

All of the previous references deal with systems of ordinary differential equations, without any notion of control. In this article we apply averaging on control systems, specifically in optimal control problems (OCPs). In this context, the question is: given a control that is optimal for an averaged OCP, will it be almost optimal for the non averaged OCP? In this case, it is not efficient to average a dynamical system and then add a control. Indeed, the "optimal" control obtained this way would vary slowly. A very simple counter example on a LQ problems shows in [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF] that, in order to obtain a sub optimal control, the optimal control computed in the averaged problem must include fast varying components.

An early work [START_REF] Balachandra | An averaging theorem for two-point boundary value problems with applications to optimal control[END_REF] applies averaging to two point boundary value problems, but its application to optimal control is limited to what is essentially the LQ case. In [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF], the method of averaging is applied to optimal control, both in open loop (in the periodic case) and in closed loop (study of the Hamilton-Jacobi equation under an ergodicity assumption). The study of the HJB equation is improved in [START_REF] Barron | Averaging in lagrange and minimax problems of optimal control[END_REF]. Observe that, in these two references, the horizon is finite and the oscillatory input is fast. By contrast, the references [START_REF] Gaitsgory | Occupational measures formulation and linear programming solution of deterministic long run average problems of optimal control[END_REF][START_REF] Gaitsgory | Linear programming approach to deterministic long run average problems of optimal control[END_REF][START_REF] Finlay | Duality in linear programming problems related to deterministic long run average problems of optimal control[END_REF][START_REF] Finlay | Duality in linear programming problems related to deterministic long run average problems of optimal control with applications to periodic optimization[END_REF] consider an optimal control problem "in the long run" with averaging techniques. The convergence of the optimal cost is proved, but there is no study of the suboptimality of the optimal control of the averaged OCP when used in the original OCP. Optimal control of celestial objects (namely the optimal control of low thrust engines in space) have been studied from a geometric view point and in the periodic case in [START_REF] Caillau | Energy minimization of single input orbit transfer by averaging and continuation[END_REF][START_REF] Caillau | Averaging and optimal control of elliptic keplerian orbits with low propulsion[END_REF][START_REF] Bonnard | Riemannian metric of the averaged energy minimization problem in orbital transfer with low thrust[END_REF]. Averaging has also been used for similar problems [START_REF] Geffroy | Optimal low-thrust transfers with constraints-generalization of averaging techniques[END_REF][START_REF] Ch | Optimal control for engines with electro-ionic propulsion under constraint of eclipse[END_REF] in a spirit that is close to [START_REF] Chaplais | Averaging and deterministic optimal control[END_REF].

To be complete, averaging has been used is stochastic optimal control, notably of Markov chains (see for instance [START_REF] Hernández | Average cost optimal policies for markov control processes with borel state space and unbounded costs[END_REF][START_REF] Vargas | Average optimal stationary policies: convexity and convergence conditions in linear stochastic control systems[END_REF][START_REF] Tsai | Averaging, aggregation and optimal control of stochastic hybrid systems with singularly perturbed morkovian switching behavior[END_REF]). Indeed, when there exists a cycle in a discrete state Markov chain that has high transition probabilities, then this cycle is gone through very fast and averaging can be applied.

A common feature of the previous literature is that

• it relies on a periodicity or ergodicity assumption

• it only provides asymptotic results.

By contrast, here

• the averaged cost and dynamics are obtained by numerical averaging on contiguous windows over the horizon [0, T ]. These functions do not need to be periodic.

• we define a number α, which represents how close the original and averaged cost and dynamics are in a numerical "weak" sense, provided the functions are regular enough.The number α can by made small by using small windows. The number α only depends on the averaged problem solution and it can be defined for any smooth OCP.

• the number α 2 , multipied by a number that depends essentially on the regularity and the convexity of the original functions, provides an error estimate between the optimal cost of the original problem, and the original cost obtained by using the optimal control of the averaged problem. This error estimate holds for any OCP provided that α ≤ β 2k J1 . The number β is part of the convexity assumption and k J1 depends on the regularity of the original functions.

• for controls that are better than the optimal control of the averaged problem, error estimates on the trajectories and controls are exhibited and are proportionnal to α.

The paper is organized as follows. Section 2 presents the original OCP. It then presents windowed averaging for functions, differential equations and control systems. It then presents the averaged OCP that is studied in this paper. This is where α is defined. Section 3 makes formal expansions in α of the state and of the costate of the nominal problem around the state and costate of the averaged problem. Auxiliary variables are introduced there. Section 4 introduces an auxiliary problem of optimization as well as new auxiliary variables. The main assumptions (bounded derivatives and convexity) are given before we state the auxiliary problem. The main result is exposed in section 5. It is a result on the control cost, trajectories and the optimal control. Section 6 is devoted to the proof of the main theorem. A conclusion is presented is section 7.

Note that, for a fluent reading of the paper, the detailed computations are in the appendices.

2 Problem statement

Nominal Problem

We wish to minimize the following Optimal Control Problem (OCP) :

min u T 0 L(x, u, t)dt (1)
where x is is a finite dimensional state which satisfies the dynamics :

dx dt = f (x, u, t) , x(0) = χ 0 ( 2
)
where u is an integrable finite dimensional, unconstrained, control. The assumptions on f and L are given in section 4.1.

Averaged Problem

It is well known that, if the dependency of f or L with respect to time contains fast oscillations, the nominal problem may be difficult to solve. To avoid this, we define an averaged problem with the help of a very simple low pass filter. It is an averaging method for general functions, i.e. not necessarily periodic.

t k = k T N , k = 0..N (3)
The averages of g on these intervals define the low pass filter LP on g:

LP [g](t) = 1 t k+1 -t k t k+1 t k g(s)ds , t ∈ [t k , t k+1 ) , k = 0..N -1 (4)
The difference g -LP [g] defines the high pass filter HP on g:

HP = Id -LP (5)
We then denote I[g] the function:

I[g](t) = t 0 HP [g](s)ds , t ∈ [0, T ] (6)
The upper bound of the function I[g] is small when N is big and g is bounded.

Proposition 1. The following bound holds for any bounded function g and any t ∈ [0, T ]:

|I[g](t)| ≤ 2 g ∞ T N (7)
Proof. Let us first prove that for any k = 0..N

t k+1 t k HP [g](t)dt = 0 Indeed: t k+1 t k HP [g](t)dt = t k+1 t k g(t)dt - t k+1 t k 1 t k+1 -t k t k+1 t k g(s)ds] dt = 0 Hence, for t ∈ [t k t k+1 [: I[g](t) = t t k HP [g](s)ds But |LP [g](s)| ≤ g ∞ , so that |HP [g](s)| ≤ 2 g ∞ . This gives finally: |I[g](t)| ≤ 2 g ∞ (t -t n+1 ) ≤ 2 g ∞ T N
Example 1 Let us suppose that g is periodical with period T N . Then LP [g] is constant equal to the mean of g over a period. The function HP [g] is a periodic signal with 0 average. The function I[g] is the periodic antiderivative of HP (g) with value 0 at 0. Its upper bound is of order T N .

Example 2 Let us suppose that g is periodical with period a small divisor of T , and N = 1. Then LP [g] is constant equal to the mean of g over the many periods in [0, T ]. The function HP [g] is a periodic signal of small period with 0 average. The function I[g] is the periodic antiderivative of HP (g) with value 0 at 0. Its upper bound is of order . This example shows that the upper bound of I[g] can be small even with a small N .

Averaging errors for ordinary differential equation

Let us consider the ordinary differential equation:

dx dt = f (x, t) , x(0) = χ 0 (8)
with f Lipschitz with respect to x and integrable with respect to t. The averaged ODE is defined as:

dx 0 dt = LP [f ](x 0 , t) , x 0 (0) = χ 0 (9)
with the low-pass filter on a function g(x, t) defined as:

LP [g](ξ, t) = 1 t k+1 -t k t k+1 t k g(ξ, s)ds , t ∈ [t k , t k+1 ) , k = 0..N -1 (10)
As a summary, LP [g] averages g with respect to time on rectangular adjacent windows, leaving the state variable unchanged.

The state x 0 is well defined and bounded because f , and thus LP [f ], is Lipschitz.

Proposition 2. Let g(x, t) be bounded Lipschitz in x with a Lipschitz constant λ g , and let's suppose that f (x, t) is bounded. Let I[g, x 0 ] the function of t defined by:

I[g, x 0 ](t) = t 0 [g(x 0 (s), s) -LP [g](x 0 (s), s)] ds (11)
where x 0 is the solution of the averaged ODE (9) Then the following bound holds:

|I[g](t)| ≤ 2( g ∞ + λ g T f ∞ ) T N (12)
Proof. see appendix A.

Proposition 3. Let x defined by the ODE (8) and let x 0 defined by the averaged ODE [START_REF] Cohen | Wavelet bases on the interval and fast algorithms[END_REF]. Let λ the Lipschitz constant of f in x and α = sup t∈[0,T ] (|I[f ](t)|). Then the following bound holds:

x -x 0 ∞ ≤ α e λT -1 λ (13)
Proof. We have

|x(t) -x 0 (t)| = t 0 f (x(s), s) -f (x 0 (s), s) + f (x 0 (s), s) -LP [f ](x 0 (s) ds ≤ λ t 0 |x(s) -x 0 (s)|ds + α
which leads to the result thanks to the Gronwall lemma.

Low pass filtering of a controlled system

In the scope of this article, we define, for the integer N > 0, and any functions u(t) and g(x, u(t), t), an averaged function LP

[g, u] LP [g, u](x, t) = N T (k+1)T N kT N g(x, u(s), s)ds f or t ∈ kT N , (k + 1)T N and k ∈ [0, N -1] (14)
As a summary, LP [g] averages g with respect to time (this includes the open loop control) on rectangular adjacent windows, leaving the state (or costate) unchanged.

Statement of the averaged problem

We define the averaged OCP which minimizes the cost

J 0 (v) = T 0 LP [L, v](y(s), s)ds ( 15
)
where y is the state defined by

dy dt = LP [f, v](y, t) , y(0) = ξ 0 (16)
which is a well defined differential equation.

Assumption 1. The averaged OCP admits an optimal control u 0 with a corresponding trajectory x 0 .

The trajectory x 0 is defined by the ODE:

dx 0 dt = LP [f, u 0 ](x 0 , t) , x 0 (0) = ξ 0 (17)

Stationnarity condition for the averaged problem

Theorem 1. Let H be the Hamiltonian of the nominal problem:

H(x, u, p, t) = L(x, u, t) + p f (x, u, t) (18)
Let p 0 be the costate of the averaged problem, defined along the optimal trajectory x 0 by the ODE with final condition:

dp 0 dt = -LP ∂H ∂x , u 0 (x 0 , p 0 , t) , p 0 (T ) = 0 ( 19
)
Then the following stationarity condition holds:

∂H ∂u (x 0 (t), u 0 (t), p 0 (t), t) = 0 a.e. ( 20
)
where a.e. stands for almost everywhere in t ∈ [0, T ].

Proof. see appendix B

Introduction of a small α

For a function g(x, u, t), the difference g(x 0 (t), u 0 (t), t) -LP ]g, u 0 ](x 0 (t), t) is the result of high pass filtering HP ]g, u](x 0 (t), t). Let's define the antiderivative:

I[g, x 0 , u 0 ](t) = t 0 HP [g, u 0 ](x 0 (s), s)ds (21)
Then, as a consequence of Proposition 2, the following bound holds:

|I[g, x 0 , u 0 ](t)| ≤ (2 g ∞ + λ g T f ∞ ) T N (22)
In other words, I[g, x 0 , u 0 ] ∞ can be made small if N is large, f is bounded and g is bounded and Lipschitz.

Let's define similarly the backwards antiderivative

I T [g, x 0 , u 0 ](t) = T t HP [g, u 0 ](x 0 (s), s)ds (23)
Then, with a similar proof, the same inequality holds:

|I T [g, x 0 , u 0 ](t)| ≤ (2 g ∞ + λ g T f ∞ ) T N (24)
Consequently, I T [g, x 0 , u 0 ] ∞ can also be made small if N is large, f is bounded and g is bounded and Lipschitz.

Definition 2. For the rest of that document, we consider the small number α defined as:

α = sup I[f, x 0 , u 0 ] ∞ , I T ∂H ∂x , (x 0 , p 0 ), u 0 ∞ ( 25
)
This number is small because N is big and assumption 2 below holds (bounded functions and their derivatives).

A priori expansions and definitions of auxiliary variables

3.1 Notations Definition 3. We denote the following variables from the state x (x 0 ), the control u (u 0 ) and the costate p (p 0 ): σ = (x, u) , σ 0 = (x 0 , u 0 ) w = (x, u, p) , w 0 = (x 0 , u 0 , p 0 ) Note: p is the costate of the nominal problem defined by the ODE with final condition:

dp dt = - ∂H ∂x (x, u, p, t) , p(T ) = 0 (26)
Definition 4. We denote the derivatives up to second order of functions with respect to the variables x, u or σ with indexes, on the model:

f x = ∂f ∂x H uu = ∂ 2 f ∂u 2 H σσ = H xx H xu H ux H uu

Formal expansion in α

The variable α defined by the equation ( 25) is a small quantity. We thus develop the state x and the costate p at the first order in α.

Expansion in the state variable

The state variable x is the solution of the original dynamics equation [START_REF] Arnold | Geometrical Methods in the Theory o f Ordinary Differential Equations[END_REF]. It is developed on the first order in α, as well as the original control u:

x = x 0 + αx 1 (27) u = u 0 + αu 1 (28)
Note that the redundant definitions of x 0 and u 0 are consistent, as will be seen later.

We then distribute the coefficients in low (LP ) and high (HP ) frequencies signals:

x 0 = x0 + x0 , x 1 = x1 + x1 (29)
Because of the equations ( 27) and (29), and because the derivative of the high frequency signal x1 is in

1 α ( 1 )
, the derivative of x has the following expansion in α:

dx dt = d x0 dt + d x0 dt + α d x1 dt + α d x1 dt (30)
But, because x is the solution of the original dynamics equation (2), using the developments in α of x and u, and thanks to assumption 2, we have another development of dx dt at the first order in α:

dx dt = f (x 0 , u 0 , t) + α(f u (x 0 , u 0 , t)u 1 + f x (x 0 , u 0 , t)x 1 ) (31)
Consequently, identifying the zero order terms in equations (30) and (31), we have:

d x0 dt + d x0 dt + α d x1 dt = f (x 0 , u 0 , t)
But by definition of x0 as the low frequency part of x 0 , we have:

d x0 dt = LP [f, u 0 ](x 0 , t)
Consequently, by definition of HP , we have:

d x0 dt + α d x1 dt = HP [f, u 0 ](x 0 , t)
with initial value 0, that derives in:

x0 + α x1 = I[f, x 0 , u 0 ] But I[f, x0 , u 0 ]
is of order 1 in α as x0 is of order 0. Thus x0 = 0, that gives x 0 = x0 , and thus the definitions of x 0 and u 0 are consistent. Moreover, we have a definition of x1 :

α x1 = I[f, x 0 , u 0 ] (32) so that the derivative of x1 is 1 α HP [f, u 0 ](x 0 , t), that is in 1 α .
A consequence of the definition of α x1 in (32) and α in definition 2 is that:

x1 ∞ ≤ 1 (33)

Expansion in the costate variable

The costate variable p is the solution of the costate dynamics equation with ending condition (26). We develop it at the first oder in α:

p = p 0 + αp 1 (34)
We then distribute the coefficients in low (LP ) and high (HP ) frequencies signals:

p 0 = p0 + p0 , p 1 = p1 + p1 (35)
Because of the equations ( 34) and ( 35), and because the derivative of the high frequency signal p1 is in 1 α , the derivative of p has the following expansion in α:

dp dt = d p0 dt + d p0 dt + α d p1 dt + α d p1 dt (36)
On the other hand, p is the solution of the costate dynamics equation with ending condition (26). Moreover, we have defined the developments in α of x, p and u.

Thus we have another development of dp dt at the first order in α:

dp dt = -H x (x 0 , u 0 , p 0 , t) + α(-H xu (x 0 , u 0 , p 0 , t)u 1 -H xx (x 0 , u 0 , p 0 , t)x 1 -f x (x 0 , u 0 , t)p 1 ) (37)
Consequently, identifying the zero order terms in equations ( 36) and (37), we have:

d p0 dt + d p0 dt + α d p1 dt = -H x (x 0 , u 0 , p 0 , t)
But by definition of p0 as the low frequency part of p 0 , we have:

d p0 dt = -LP [H x , u 0 ](x 0 , p 0 , t)
Consequently, by definition of HP , we have:

d p0 dt + α d p1 dt = -HP [H x , u 0 ](x 0 , p 0 , t)
with final value 0, that derives in:

p0 + α p1 = I T [H x , (x 0 , p 0 ), u 0 ]
But I[-H x , x 0 , p 0 , u 0 ] is of order 1 in α as p0 is of order 0. Thus p0 = 0, that gives p 0 = p0 , and thus the definitions of p 0 are consistent. Moreover, we have a definition of p1 :

α p1 = I T [H x , (x 0 , p 0 ), u 0 ] (38) so that the derivative of p1 is -1 α HP [H x , u 0 ](x 0 , p 0 , t), that is in 1 α .
A consequence of the definition of α p1 in (38) and α in definition 2 is that:

p1 ∞ ≤ 1 (39)
4 Auxiliary Problem

Assumptions

Assumption 2 (smoothness). The derivatives, up to the third order, of f and L with respect to x and u are bounded by some k > 0.

A consequence of that assumption is the the value α defined by equation ( 25) is well defined and small if N is sufficiently big (or the problem is periodic of small period).

Another consequence is that:

Proposition 4. f σσ (X, U, t) Y V 2 ≤ k(|Y | + |V |) 2 (40)
for any (X, U, Y, V ).

Proof.

f σσ (X, U, t) Y V 2 = |f xx (X, U, t)Y 2 + 2f xu (X, U, t)Y V + f uu (X, U, t)V 2 | ≤ k(|Y | + |V |) 2
Moreover, as p 0 is the solution of the ODE with final condition [START_REF] Strang | A Fourier analysis of the finite element variational method[END_REF], it is differentiable and thus continuous of the bounded interval [0, T ], so that it is bounded. So that another consequence of the assumption 2 is so: Proposition 5. The hamiltonian H(x, u, p 0 , t) and its derivatives up to the third order in u and x are bounded by a constant K.

Proof. Take K = (1 + p 0 ∞ )k.
A consequence of that is:

Proposition 6. H σσ (X, U, p 0 t) Y V 2 ≤ K(|Y | + |V |) 2 (41)
for any (X, U, Y, V ).

Proof. Similar proof as for f σσ .

Assumption 3 (convexity).

There exists β > 0 so that for any (x, u), the following holds:

H uu (x, u, p 0 , t) ≥ βId (42)
and

H xx -H xu (H uu ) -1 H ux (x, u, p 0 , t) ≥ 0 (43)
The consequence of equation ( 42) is that H uu is invertible and H -1 uu ∞ ≤ 1 β on any (x, u, p 0 ). The consequence of equation ( 43) is :

Proposition 7. H σσ (x, u, p 0 , t) ≥ 0 (44)
Proof. We make a proof by contradiction. Let's suppose that (44) is not true. Then there exists a negative eigenvalue -γ of H σσ (x, u, p 0 , t), that is there exists an eigenvector y v so that:

H σσ (x, u, p 0 , t) y v = -γ y v
This implies that the two following equations hold:

H xx (x, u, p 0 , t)y + H xu (x, u, p 0 , t)v = -γy (45)
H ux (x, u, p 0 , t)y + H uu (x, u, p 0 , t)v = -γv (46)
But (γId + H uu) ≥ (γ + β)Id > 0, so that it is invertible and the equation ( 46) can be solved in v, giving:

v = (γId + H uu (x, u, p 0 , t)) -1 H ux y (47)
Then, replacing v by its value in the equation ( 45), we have:

(H xx (x, u, p 0 , t) -H xu (x, u, p 0 , t)(γId + H uu (x, u, p 0 , t)) -1 H ux (x, u, p 0 , t))y = -γy (48)
But as γ > 0, we have the succession of inequalities:

(γId + H uu (x, u, p 0 , t)) ≥ H uu (x, u, p 0 , t) then (γId + H uu (x, u, p 0 , t)) -1 ≤ H -1 uu (x, u, p 0 , t) and then (H xx (x, u, p 0 , t) -H xu (x, u, p 0 , t)(γId + H uu (x, u, p 0 , t)) -1 H 0ux (x, u, p 0 , t)) ≥ (H xx -H xu H -1 uu H ux ) (x, u, p 0 , t) ≥ 0
Thus the equation (48) can not hold, because -γ < 0 can not be an eigenvalue, and equation ( 44) is proved by contradiction.

Auxiliary Problem Statement

Definition 5. Let's define the following notations:

H 0σσ = H σσ (w 0 , t) = H 0xx H 0xu H 0ux H 0uu and f 0x = f x (σ 0 , t) , f 0u = f u (σ 0 , t)
Then we define the auxiliary problem as the linear quadratic OCP with state y and control v:

dy dt = f 0x (y + x1 ) + f 0u v , y(0) = 0 (49) J 1 (v) = T 0 1 2 y v H 0σσ y v + p1 (f 0x y + f 0u v) dt (50)
Proposition 8. There exists an optimal cost v 1 for the auxiliary problem.

Proof. It is a convex linear quadratic problem because H 0σσ is non-negative (equation ( 44)).

Let's denote v 1 an optimal control, y 1 the trajectory corresponding to v 1 and q 1 the corresponding costate. Then y 1 follows the dynamics of the auxiliary problem:

dy 1 dt = f 0x (y 1 + x1 ) + f 0u v 1 , y 1 (0) = 0 (51)
The hamiltonian of the auxiliary problem expands in:

H 1 (y, v, q) = 1 2 H 0xx y 2 + 2H 0xu yv + H 0uu v 2 + p1 (f 0x y + f 0u v) + q 1 [f 0x (y + x1 ) + f 0u v] (52)
Thus the stationarity condition ∂H 1 ∂u (y 1 , v 1 , p 1 ) = 0 may be written the following way:

H 0xu y 1 + H 0uu v 1 + ( p1 + q 1 )f 0u = 0 (53)
Moreover, the costate q 1 of the auxiliary problem follows the dynamics:

dq 1 dt = -H 0xx y 1 -H 0xu v 1 -( p1 + q 1 )f 0x , q 1 (T ) = 0 (54)
Moreover, we have:

Proposition 9. y 1 , v 1 and q 1 are bounded by a constant M .

Proof. The auxiliary problem is smooth and convex.

More auxiliary variables and their upper bounds

Definition 6. For any u ∈ L 2 [0,T ] ,
x is the trajectory of the nominal dynamics (2). Let's then define the following notations:

δx = x -x 0 x = δx -α x1 δu = u -u 0 ũ = δu δσ = (δx, δu) ρ(λ, µ) = σ 0 + λµδσ
Where (u 0 , x 0 ) is a solution of the averaged problem and α x1 is I[f, x 0 , u 0 ] (equation (32)).

Upper bounds for r and v Definition 7. We define the following data:

r = x -αy 1 , v = ũ -αv 1 Z[λ, µ](t) = v + [H -1 uu H ux ](ρ(λ, µ), p 0 , t)(r + α x1 ) z 2 = 1 0 1 0 λ Z[λ, µ] 2 2 dλdµ
Definition 8. We define the following constants:

k r1 = 4k 2 T e 2k 1+ K β , k r2 = K β + 1 2 (1 + 2M ) 2 α 2 1 + K β 2 e k 1+ K β -1 2 k v1 = 6 2 + K 2 T k r1 β 2 , k v2 = 6K 2 T (k r2 + 1) β 2
with:

-k is introduced in assumption 2 abound the bounded derivatives of f (x, u, t) and L(x, u, t).

-K = (1 + p 0 ∞ )k -M
is the upper bound of the optimal trajectory of the auxiliary problem introduced in section 4.2 (proposition 9).

-α is the small quantity defined in equation ( 25).

-β is the convexity constant of H uu introduced in equation ( 42).

Proposition 10. The following inequalities hold:

r 2 ∞ ≤ k r1 z 2 + k r2 α 2 (55)
and:

v 2 2 ≤ k v1 z 2 + k v2 α 2 (56) Proof. see appendix C.1
Upper bounds for r -r 1 and r 1 Definition 9. r 1 is defined by the following dynamics:

dr 1 dt = f (r 1 + x 0 + α( x1 + y 1 ), v + u 0 + αv 1 , t) -f (x 0 + α( x1 + y 1 ), u 0 + αv 1 , t) , r 1 (0) = 0 (57)
Definition 10. We define the following constants:

k r5 = T 2 (1 + 2M ) 2 (e kT -1) , k r3 = 2k r1 , k r4 = 2(k r2 + k 2 r5 α 2 )
with k, K and M as in definition 8.

Proposition 11. The following inequalities hold:

r -r 1 ∞ ≤ k r5 α 2 (58) and r 1 2 ∞ ≤ k r3 z 2 + k r4 α 2 (59)
Proof. see appendix C.2

Main theorem

Definition 11. We define the following constants:

k J1 = 4 √ 3KM T k r1 + 3k 2 + KM T k r3 + 3k 2 + (2 √ 3 + 1)KM k v1 k J2 = 2 √ 3KT k r5 + k(1 + 2M ) + 2KT M (M + 1) + 3k 2 T k r2 + kM + 4 √ 3KM T k r4 + 3k 2 + (2 √ 3 + 1)KM k v2 k J0 = e kT -1 + K 2k 2 T e kT -1 2 + k 2 α k J = k J2 + k J0 k x = 2 2k J k r1 β + k r2 + (M + 1) 2 k u = 2 2k J k v1 β + k v2 + T 2 M 2
where:

-k, K, M , α and β are as in definition 8.

-k r1 , k r2 , k v1 and k r2 are defined in definition 8 -k r3 , k r4 and k r5 are defined in definition 10.

Note that these constants depend only of k, K, M , α, β, and the horizon T .

Assumption 4. α ≤ β 2k J1
Theorem 2 (Main Theorem). Considering the nominal problem in section 2.1, let H(x, u, p, t) its Hamiltonian, and let J * = inf u J(u) be its infimum cost.

Let u 0 a solution of the averaged problem described in section 2.2.4 with its trajectory x 0 . Such a solution exists by assumption 1.

Let α be the small quantity defined in equation ( 25) and let β be the constant introduced in assumption 3.

Let the set of constants (k J1 , k J , k x , k u ) introduced in definition 11. Then, under the set of assumptions listed in section 4.1 and the assumption 4, the following inequalities hold:

-the suboptimality of the real system commanded by u 0 is limited to:

J * ≤ J(u 0 ) ≤ J * + k J α 2 (60)
-any trajectory x of the nominal problem for a u better than u 0 (J(u) ≤ J(u 0 )), is close to (x 0 , u 0 ), with:

x -x 0 ∞ ≤ k x α (61) u -u 0 2 ≤ k u α (62)
6 Proof of the main result

Proof Process

To prove the main theorem, we proceed the following way. The section 6.2 is devoted to the search of a lower bound of any real cost J(u) of the nominal problem. That lower bound contains two integral terms that do not depend on u, a term in z 2 that is the only one depending on u, and a term in α 2 , that dos not depend on u either.

The section 6.3 is devoted to the search of an upper bound of the real cost J(u 0 ) of the nominal problem controlled by u 0 . That upper bound contains the same two integral terms as in the lower bound of J(u) and a term in α 2 .

Then the section 6.4.1 uses the assumption α ≤ β k J1 to obtain a lower bound of J(u) independent of u, so that it is also a lower bound for J * . That lower bound is combined with the upper bound of J(u 0 ) to prove the suboptimality in α 2 of L(u 0 ) stated in the equation (60) of the first part of the main theorem 2.

Then the section 6.4.1 uses the stronger assumption α ≤ β 2k J1 to obtain, for any u better than u 0 , i.e. so that J(u) ≤ J(u 0 ), an upper bound of z 2 . With that bound of z 2 , upper bounds for x -x 0 ∞ and u -u 0 2 are found in the equations ( 61) and (62) of the second part of the main theorem, with the help of the definitions and bounds of r and v in section 4.3.

6.2 Lower bound on the real cost of the nominal problem 6.2.1 Expansion of the real cost J(u) Proposition 12. The cost J(u) = T 0 L(x, u, t)dt for any command u ∈ L 2 [0,T ] expands the following way:

J(u) = T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + T 0 HP [H x , u 0 ](x 0 , p 0 , t)xdt + T 0 1 0 1 0 λH σσ (ρ(λ, µ), p 0 , t)(δσ) 2 dλdµdt (63)
Proof. see appendix D.

Lower

Bound of the third term of the expansion of J(u) in equation ( 63) Proposition 13. The following inequality holds:

T 0 HP [H x , u 0 ](x 0 , p 0 , t)xdt ≥ α T 0 p1 [f 0x r 1 + f 0u v]dt - 3k 2 (T k r3 + k v1 )αz 2 (64)
-2KT k r5 + k(1 + 2M ) + k 2 (3(T k r4 + k v2 ) + 4T (1 + 2M ) 2 )α α 2
Proof. See appendix E.

Lower

Bound of the fourth term of the expansion of J(u) in equation ( 63) Proposition 14. The following inequality holds:

T 0 1 0 1 0 λH σσ (ρ(λ, µ), p 0 , t)(δσ) 2 dλdµdt ≥ -2KM T M + √ 3(2T k r2 + k v2 + T (M 2 + 2))α α 2 + β -2 √ 3KM (2T k r1 + k v1 )α z 2 +α T 0 p1 [y 1 , v 1 ] H 0σσ r + α x1 v dt (65)
Proof. see appendix F.

6.2.4

Bound in absolute value for the sum of the integral terms of the right hand sides of equations ( 64) and (65) Proposition 15. Let's define R as the sum of the integral terms of the right hand sides of equations ( 64) and (65):

R = α T 0 p1 f 0x r 1 + f 0u v + [y 1 , v 1 ] H 0σσ r + α x1 v dt
Then the following inequality holds:

|R| ≤ KM [T k r3 + k v1 ] + [2KT M + (2KT M k r5 + kT M k r4 + kM k v2 )α]α 2 (66)
Proof. see appendix G.

6.2.5 Lower Bound of the real cost J(u) Lemma 1. A lower bound of the cost J(u) of the nominal system for any control u is given by:

J(u) ≥ T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + (β -k J1 α)z 2 -k J2 α 2 (67)
where k J1 and k J2 are defined in definition 11.

Proof. This is a consequence of equations ( 64), ( 65) and (66)

6.3 Upper bound on the cost of the nominal system controlled by u 0 6.3.1 Expansion of the cost J(u 0 ) Definition 12. Let x 0 be the trajectory of the nominal problem controlled by u 0 . It is defined by the dynamics:

dx 0 dt = f (x 0 , u 0 , t) , x 0 (0) = χ 0 ( 68
)
For that trajectory, we set the notations:

δx 0 = x 0 -x 0 x0 = δx 0 -α x1 ρ 0 (λ, µ) = x 0 + λµ δx 0
Where (u 0 , x 0 ) is the solution of the averaged problem and α x1 is I[f, x 0 , u 0 ] (equation ( 32)).

Proposition 16. The cost J(u 0 ) = T 0 L(x 0 , u 0 , t)dt for the optimal command u 0 of the averaged problem expands the following way:

J(u 0 ) = T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + T 0 HP [H x , u 0 ](x 0 , p 0 , t)x 0 dt + T 0 1 0 1 0 λH xx (ρ 0 (λ, µ), u 0 , p 0 , t)(δx 0 ) 2 dλdµdt (69)
Proof. It is a consequence of proposition 12 with u = u 0 , so that δu = 0.

6.3.2 Comparison of x 0 and x 0 Proposition 17. The following inequality hold for δx 0 = x 0 -x 0 :

|δx 0 | ≤ α kT e kT -1 (70)
Proof. As δx 0 = x 0 -x 0 , x 0 follows the dynamics (68) and x 0 follows the averaged dynamics (17), we have the following integral equation:

δx 0 = t 0 [f (x 0 , u 0 , t) -LP [f, u 0 ](x 0 , t)]dt = t 0 [f (x 0 , u 0 , t) -f (x 0 , u 0 , t)]dt + t 0 [f (x 0 , u 0 , t) -LP [f, u 0 ](x 0 , t)]dt
Thus the following inequality holds:

|δx 0 | ≤ f x ∞ t 0 |δx 0 |ds + |I[f, x 0 , u 0 ](t)| ≤ kT t 0 |δx 0 |ds + α
Equation (70) follows from Gronwall lemma.

T 0 HP [H x , u 0 ](x 0 , p 0 , t)x 0 dt ≤ e kT -1 + k 2 α α 2 (71)
Proof. By definition of x0 , we have x0 = x 0 -x 0 -α x1 , so that:

dx 0 dt = f (x 0 + δx 0 , u 0 , t) -LP [f, u 0 ](x 0 , t) + HP [f, u 0 ](x 0 , t) = f (x 0 + δx 0 , u 0 , t) -f (x 0 , u 0 , t)
Thus, by Taylor expansion of f (x 0 , u 0 , t) with integral remainder, we have:

dx 0 dt = f x (x 0 , u 0 , t)δx 0 + 1 0 1 0 λf xx ((x 0 + λµ(δx 0 ), u 0 , t)(δx 0 ) 2 dλdµ
Thus, thanks to proposition 17, we have:

dx 0 dt ≤ (e kT -1)α + k 2 α 2 (72)
But an integration by part, together with the fact that HP [H x , u 0 ](x 0 , p 0 , t) = -α d p1 dt and that x0 (0) = p1 (T ) = 0 leads to:

T 0 HP [H x , u 0 ](x 0 , p 0 , t)x 0 dt = α T 0 p1 dx 0 dt dt (73)
λH xx (ρ 0 (λ, µ), u 0 , p 0 , t)(δx 0 ) 2 dλdµdt ≤ K 2k 2 T e kT -1 2 α 2 (74)
Proof. This is a consequence of the proposition 17.

Upper Bound of J(u 0 )

Inserting equations ( 71) and (74) into equation (69 leads to:

Lemma 2. An upper bound of the cost J(u 0 ) of the nominal system controlled by u 0 is given by:

J(u 0 ) ≤ T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + k J0 α 2 (75)
where k J0 is defined in definition 11.

Proof of the main theorem

The main theorem is proved in two steps. The first step compares the costs to estimate the suboptimality of the real system cotrolled by u 0 . The second step compares the trajectories and controls better than u 0 to estimate how close they are from the trajectory and control dealed by u 0 .

6.4.1 Comparison of the real cost controlled by u 0 and the infimum cost of the real system Let's now use the assumption α ≤ β 2k J1 of the first part of the Main Theorem 2 into the equation (67) of Lemma 1. The term in z 2 is then non negative, and we get the lower bound independent of u:

J(u) ≥ T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt -k J2 α 2
As that lower bound holds for any u, it is also a lower bound for the infimum cost J * = inf u J(u):

J * ≥ T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt -k J2 α 2 so that: T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt ≤ J * + k J2 α 2
Inserting this equation into the equation (75) of Lemma 2, together with the fact that J(u 0 ) ≥ J * , by definition of J * , proves the suboptimality equation (60) in he Main Theorem 2, since k J = K J0 + k J2 .

Comparison of the controls and trajectories with and without u = u 0

Let's consider a control u better than u 0 , i.e. such that J(u) ≤ J(u 0 ).

Let's now use in a stronger manner the assumption α ≤ β 2k J1 of the second part of the Main Theorem 2 into the equation (67) of Lemma 1. The coefficient z 2 is then lower than β 2 , and we get the lower bound dependent of z 2 , that depends on u:

J(u) ≥ T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + β 2 z 2 -k J2 α 2
Thus, together with the equation (75) of Lemma 2, we have the list of inequalities:

T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + β 2 z 2 -k J2 α 2 ≤ J(u) ≤ J(u 0 ≤ T 0 L(x 0 , u 0 , t)dt + α T 0 H 0x x1 dt + k J0 α 2
So that, with the definition of k J = K J0 + k J2 , we have

z 2 ≤ 2k J β α 2 (76)
On the other hand, by definition of r, we have:

x -x 0 = r + α( x1 + y 1 )
so that, together with the equation (55): 76) into that equation leads to equation (61) of the second part of the Main Theorem 2. Now let's consider the fact that, by definition of v: 76) into that equation leads to equation (62) of the second part of the Main Theorem 2.

x -x 0 2 ∞ ≤ 2 r 2 ∞ + α 2 (1 + M ) 2 ≤ 2 k r1 z 2 + (k r2 (1 + M ) 2 )α 2 Introducing equation (
u -u 0 = v + αv 1 so that: u -u 0 2 2 ≤ 2 v 2 2 + α 2 T 2 M 2 ≤ 2 k v1 z 2 + (k v2 T 2 M ) 2 α 2 Introducing equation (

Conclusion

We have shown that the method of averaging can be used very simply by performing averages of the dynamics on adjacent intervals. Its efficiency, notably in optimal control, is measured by α and by the convexity β of the cost function. Using α as a measure of the efficiency of averaging amounts to saying that, from this point of view, the simple integrator acts as reference for all state space models. An important point is that, provided that α ≤ β 2k J1 , the estimates (60,61,62) in the main theorem 2 hold for any system.

This method of averaging has been applied to the guidance of a low thrust satellite in the non keplerian case [START_REF] Goumri | Transfert orbital d'un engin à faible poussée[END_REF]. In this case, we cannot use periodic averaging on orbits because of the influence of the sun and of the moon, which have different periods.

A Proof of Proposition 2 about the upper bound of I[g, x 0 ](t) I[g, x 0 ](t) is defined by equation ( 11) and x 0 is the solution of the averaged ODE ( 9).

Let's fix t and let K be so that t ∈ [t K , t K+1 ). Then we have:

I[g, x 0 ](t) = K-1 k+0 I k [g, x 0 ](t k+1 ) + I K [g, x 0 ](t) ( 77
)
where I j [g, x 0 ](τ ) is defined for τ ∈ [t j , t i+1 ) as:

I j [g, x 0 ](τ ) = τ t j g(x 0 (s), s) - 1 t j+1 -t j t j+1 t j g(x 0 (s), σ)dσ ds
We can develop the terms I k [g, x 0 ](t k+1 ) and I K g, x 0 (t) in subtracting and adding the values at t k .

I k [g, x 0 ](t k+1 ) = t k+1 t k g(x 0 (t k , s) + 1 t k+1 -t k t k+1 t k g(x 0 (t k ), σ)dσ ds + t k+1 t k [g(x 0 (s), s) -g(x 0 (t k ), s)]ds + 1 t k+1 -t k t k+1 t k t k+1 t k [g(x 0 (s), σ) -g(x 0 (t k ), σ)]dσds
The first term is equal to 0.

The terms integrated once and twice in the second and third terms are both lower or equal to λ g x 0 (s) -x 0 (t k because g is Lipschitz in x, with Lipschitz constant λ g .

Moreover, because x 0 is the solution of the averaged ODE (9), we have:

x 0 (s) -x 0 (t k ≤ f ∞ (s -t k )
Thus, we have:

|I k [g, x 0 ](t k+1 )| ≤ 2λ g f ∞ t k+1 t k (s -t k )ds
That is:

|I k [g, x 0 ](t k+1 )| ≤ 2λ g f ∞ T N 2 (78)
Let's now develop I K [g, x 0 ](t) in an analog way starting from t K :

I K [g, x 0 ](t) = t t K g(x 0 (t K , s) + 1 t K+1 -t K t K+1 t K g(x 0 (t K ), σ)dσ ds + t t K [g(x 0 (s), s) -g(x 0 (t K ), s)]ds + 1 t K+1 -t K t t K t K+1 t K [g(x 0 (s), σ) -g(x 0 (t K ), σ)]dσds
The first term is not 0 here, but it is lower or equal to 2 g ∞ T N , and the other terms are bounded as the ones of I k [g, x 0 ](t k+1 ).

Thus we have:

|I K [g, x 0 ](t)| ≤ g ∞ T N + 2λ g f ∞ T N 2 (79)
Finally, because K + 1 ≤ N , the inequations ( 78) and ( 79) yield to the proposition 2.

B Proof of the stationarity condition of the averaged problem B.1 Averaged two boundaries problem

Let u 0 be the optimal control of the averaged problem ( 16) and let x 0 be the corresponding trajectory. x 0 is defined by the ODE with initial condition [START_REF] Poincaré | Les méthodes nouvelles de la mécanique céleste[END_REF]. Let p 0 be the costate of the optimal trajectory, defined by the ODE with final condition [START_REF] Strang | A Fourier analysis of the finite element variational method[END_REF], with H the hamiltonian [START_REF] Sanders | Averaging Methods in Nonlinear Dynamical Systems[END_REF].

The system constituted of the of equations ( 17) and ( 19) is a two boundaries problem. It is defined as the two boundaries problem corresponding to the averaged optimal control problem.

B.2 First variation in the direction of δu

Let δu ∈ L 2 [0,T ] a scalar square integrable function on [0, T ]. Let > 0 and let u = u 0 + δu the variation of u 0 in the direction of δu. As J 0 (u 0 ) = min u (J 0 (u), the following stationarity condition holds:

∀ δu ∈ L 2 [0,T ] , dJ 0 (u ) d =0 = 0 (80)
Let x be the trajectory corresponding to u , defined by the dynamics equation 16 with v = u :

dx dt = LP [f, u ](x , t) , x (0) = ξ 0 (81)
Let δx be the variation trajectory corresponding to the direction δu given by δx = dx d =0 . Lemma 3. δx respects the following dynamics function:

d(δx) dt = LP ∂f ∂x , u 0 (x 0 , t)δx + LP ∂f ∂u δu, u 0 (x 0 , t) (82) Proof. d(δx) dt = d dt dx d =0 = d d dx dt =0 = d d LP [f, u 0 + δu](x , t) =0
Let k be so that [t k , t k+1 )). Then:

d(δx) dt = d d 1 t k+1 -t k + t k+1 t k f (x (t), u 0 (s) + δu(s), s)ds =0 = 1 t k+1 -t k + t k+1 t k d d (f (x (t), u 0 (s) + δu(s), s)) =0 ds = 1 t k+1 -t k t k+1 t k ∂f ∂x
(x 0(t),u 0 (s),s δx(t) + ∂f ∂u (x 0(t),u 0 (s),s δu(s) ds

= LP ∂f ∂x , u 0 (x 0 , t)δx + LP ∂f ∂u δu, u 0 (x 0 , t)

B.3 Proof of the stationarity result

Let's make use of Equation (80) in developing dJ 0 (u )

d =0
for a given δu.

Lemma 4. The derivative at 0 of J 0 (u ) in is related to the hamiltonian by the following equation:

dJ 0 (u ) d =0 = T 0 LP ∂H ∂u δu, u 0 (x 0 , p 0 , t)dt (83)
Proof. Let's use the equation ( 15) and then commute the differentiation and integration:

dJ 0 (u ) d =0 = d d T 0 LP [L, u ](x , t)dt =0 = T 0 d d [LP [L, u ](x , t)] =0 dt
For any k ∈ [0, N -1] and for any t ∈ [t k , t k+1 ), we have the definition [START_REF] Gaitsgory | Occupational measures formulation and linear programming solution of deterministic long run average problems of optimal control[END_REF] of LP :

LP [L, u ](x , t) = 1 t k+1 -t k t k+1 t k L(x (t), u (s), s)ds
Thus if we commute again the integration and the differentiation:

d d [LP [L, u ](x , t)] =0 = 1 t k+1 -t k t k+1 t k dL(x (t), u (s), s) d =0 ds
But by definition of x , u , δx and δu, we have:

dL(x (t), u (s), s) d =0 = ∂L ∂x (x 0 (t), u 0 (s), s)δx(t) + ∂L ∂u (x 0 (t), u 0 (s), s)δu(s)
Thus by averaging on [t k , t k+1 ], the result is (LP is linear):

d d [LP [L, u 0 ](x , t)] =0 = LP ∂L ∂x , u 0 (x 0 (t), t)δx(t) + LP ∂L ∂u δu, u 0 (x 0 (t), t)
But because of the dynamics (19) of the averaged costate p 0 , we have, with the definition (18) of the Hamiltonian:

LP ∂L ∂x , u 0 (x 0 (t), t) = - dp 0 dt -p 0 LP ∂f ∂x , u 0 (x 0 (t), t) Thus integrating d d [LP [L, u 0 ](x , t)] =0
between 0 and T , we obtain:

dJ 0 (u ) d =0 = T 0 - dp 0 dt δx- T 0 p 0 LP ∂f ∂x , u 0 (x 0 (t), t)δx(t)dt+ T 0 LP ∂L ∂u δu, u 0 (x 0 (t), t)dt (84)
Let's make an integration by part for the first term of that equation:

T 0 - dp 0 dt δx = -[p 0 δx] T 0 + T 0 d(δx) dt p 0
The variation of p 0 δx between 0 and T is null because δx(0) = 0 and p 0 (T ) = 0. Thus, with the dynamics of δx given by the Lemma 3, the following holds:

T 0 - dp 0 dt δx = T 0 p 0 LP ∂f ∂x , u 0 (x 0 , t)δx(t)dt + T 0 p 0 LP ∂f ∂u δu, u 0 (x 0 , t)dt
Let's insert this equation in the first term of equation ( 84).It results in:

dJ 0 (u ) d =0 = T 0 LP ∂L ∂u δu, u 0 (x 0 , t) + p 0 LP ∂f ∂u δu, u 0 (x 0 , t) dt
This proves the equation ( 83) by definition of the Hamiltonian.

Let's now make use of equation ( 83). Let's first fix t and let k be so that t ∈ [T k , t k+1 ). Then we have:

LP ∂H ∂u δu, u 0 (x 0 , p 0 , t) = 1 t k+1 -t k t k+1 t k ∂H ∂u (x 0 (t), u 0 (s), p 0 (t), t)δu(s)ds
Let's specialize δu as a "needle variation":

δu = t k+1 -t k η 1 [t,t+η] δv
with η > 0 so that t + η < t k+1 and δv ∈ L 2 [0,T ] . Then we have:

LP ∂H ∂u δu, u 0 (x 0 , p 0 , t) = 1 η t+η t
∂H ∂u (x 0 (t), u 0 (s), p 0 (t), t)δv(s)ds a.e.

-→ η → 0 ∂H ∂u (x 0 (t), u 0 (t), p 0 (t), t)δv(t)

More precisely, the limit is the value of the function at t everywhere the function is continue, that is for any t possibly except for a countable number of "jumps". As any countable set is negligible, the limit holds almost everywhere. Thus, because of the equations ( 80) and (83), we have:

T 0 ∂H ∂u (x 0 (t), u 0 (t), p 0 (t), t)δv(t) = 0
and this is true for any δv ∈ L 2 [0,T ] . This proves the stationnarity result:

∂H ∂u (x 0 (t), u 0 (t), p 0 (t), t) = 0 a.e.
C Proof of the inequalities on r, v, r -r 1 and r 1

dr dt = f (r + x 0 + α( x1 + y 1 ), v + u 0 + αv 1 , t) -f (x 0 + α( x1 + y 1 ), u 0 + αv 1 , t) (85)
+α 2 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 , v + u 0 + αv 1 , t) x1 + y 1 v 1 2 dλdµ
Proof. By definition of r and v, we have:

r = x -x 0 -α( x1 + y 1 ), so that : x = r + x 0 + α( x1 + y 1 ) v = u -u 0 -v 1 , so that : u = v + u 0 + v 1
Moreover:

dr dt = dx dt - dx 0 dt -α d x1 dt -α dy 1 dt = f (r + x 0 + α( x1 + y 1 ), v + u 0 + v 1 , t) -LP [f, u 0 ](x 0 , t) -HP [f, u 0 ](x 0 , t) -α[f 0x (y1 + x1 ) + f 0u v 1 ] = f (r + x 0 + α( x1 + y 1 ), v + u 0 + v 1 , t) -f (x 0 , u 0 , t) (86)
-α[f 0x (y1 + x1 ) + f 0u v 1 ]
But a Taylor expansion of f (x 0 + α( x1 + y 1 ), u 0 + αv 1 , t) is so:

f (x 0 ( x1 + y 1 ), u 0 + αv 1 , t) = f (x 0 , u 0 , t) + α[f 0x (y1 + x1 ) + f 0u v 1 ] (87) +α 2 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 , v + u 0 + αv 1 , t) x1 + y 1 v 1 2 dλdµ
Introducing equation (87) in equation ( 86) proves equation ( 85).

Proposition 21. The following inequality holds:

r ∞ ≤ ke k 1+ K β √ T Z(λ, µ) 2 + T K β + 1 2 ((1 + 2M ) 2 α α (88)
Proof. Equation (40) about the upper bound of f σσ leads to:

f σσ (x 0 + λµα( x1 + y 1 , v + u 0 + αv 1 , t) x1 + y 1 v 1 2 ≤ k(| x1 + y 1 | + |v 1 |) 2 ≤ k(| x1 | + |y 1 | + |v 1 |) 2
Hence:

f σσ (x 0 + λµα( x1 + y 1 , v + u 0 + αv 1 , t) x1 + y 1 v 1 2 ≤ k(1 + 2M ) 2 (89)
Moreover, the function f is Lipschitz in x and u with Lipschitz constant the bound of the derivatives k, so that:

|f (r + x 0 + α( x1 + y 1 ), v + u 0 + αv 1 , t) -f (x 0 + α( x1 + y 1 ), u 0 + αv 1 , t)| ≤ k(|r| + |v|) (90)
The equations ( 89) and (90) in the equation (85) give, together with the fact that r(0) = 0 the following inequality:

|r(t)| ≤ t 0 (|r(s)| + |v(s)|)ds + kT 2 (2M + 1) 2 α 2 (91)
But by definition of Z[λ, µ], we have:

v = Z[λ, µ] -[H -1 uu H ux ](ρ(λ, µ), p 0 , t)(r + α x1 ) (92) so that: |v(s)| ≤ |Z[λ, µ](s)| + H -1 uu ∞ H ux ∞ (|r(s)| + α| x1 (s)|)
But the assumption 3 proves that H uu is invertible and that H -1 uu is bounded by 1 β , so that:

|v(s)| ≤ |Z[λ, µ](s)| + K β (|r(s)| + α) (93)
Including equation(93) in equation ( 91) leads to:

|r(t)| ≤ k 1 + K β t 0 (|r(s)|)ds + kKT β α + kT 2 (2M + 1) 2 α 2 + k T 0 |Z[λ, µ](t)||dt (94)
But Cauchy property leads to:

T 0 |Z[λ, µ](t)||dt ≤ √ T Z(λ, µ)] 2 (95)
Including equation (95) in equation ( 94) lead to:

|r(t)| ≤ k 1 + K β t 0 (|r(s)|)ds + k √ T Z(λ, µ) 2 + kT K β + 1 2 (2M + 1) 2 α α (96)
Equation ( 96), together with Gronwall lemma, proves equation (88).

Proof of Equation (55) Let's take the square of equation ( 88):

r 2 ∞ ≤ 2k 2 e 2k 1+ K β T Z(λ, µ) 2 2 + T 2 K β + 1 2 ((1 + 2M ) 2 α 2 α 2
Let's now multiply by λ and integrate relatively to λ and µ between 0 and 1:

1 2 r 2 ∞ ≤ 2k 2 e 2k 1+ K β T z 2 + 1 2 T 2 K β + 1 2 ((1 + 2M ) 2 α 2 α 2
Multiplying by 2 that equation proves Equation (55).

C.1.2 Upper bound for v

Let's apply the triangular inequality for the L 2 norm to the expression of v (92):

v 2 ≤ Z[λ, µ] 2 + [H -1 uu H ux ](ρ(λ, µ), p 0 , t)(r + α x1 ) 2 But: [H -1 uu H ux ](ρ(λ, µ), p 0 , t)(r + α x1 2 ≤ √ T [H -1 uu H ux ](ρ(λ, µ), p 0 , t)(r + α x1 ) ∞ ≤ K √ T β ( r ∞ + α x1 ∞ )
Thus, taking the squares:

v 2 2 ≤ 3 Z[λ, µ] 2 + K 2 T β 2 ( r 2 ∞ + α 2 x1 2 ∞ )
Let's multiply by λ and integrate relatively to λ and µ between 0 and 1:

1 2 v 2 2 ≤ 3 z 2 + K 2 T 2β 2 ( r 2 ∞ + α 2 ≤ 3 1 + K 2 T k r1 2β 2 z 2 + K 2 T 2β 2 (1 + k r2 )α 2
Multiplying by 2 that equation proves Equation (56).

C.2 Proof of proposition 11

r follows the dynamic (85) with r(0) = 0 and r 1 follows the dynamic (57). Thus r -r 1 follows the dynamics:

d(r -r 1 ) dt = f (r + x 0 + α( x1 + y 1 ), v + u 0 + αv 1 , t) -f (r 1 + x 0 + α( x1 + y 1 ), v + u 0 + αv 1 , t) +α 2 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 , v + u 0 + αv 1 , t) x1 + y 1 v 1 2 dλdµ
Thus, in a similar way than for the upper bound of r, the following inequality holds:

|r(t) -r 1 (t)| ≤ k t 0 |r(s) -r 1 (s)|ds + kT 2 (2M + 1) 2 α 2
The equation (58) follows from Gronwall lemma. The equation ( 59) is then the consequence of equations ( 55) and ( 58), together with:

r 1 2 ∞ ≤ 2( r 2 ∞ + r -r 1 2 ∞ )
D Proof of the expansion of the real cost (proposition 12)

Proposition 22. L(x, u, t) expands the following way:

L(x, u, t) = L(x 0 , u 0 , t) + L x (x 0 , u 0 , t)δx + L u (x 0 , u 0 , t)δu + 1 0 1 0 λL σσ (ρ(λ, µ), t)(δσ) 2 dλdµ (97)
Proof. It is a Taylor expansion of L(x, u, t) with integral remainder.

Proposition 23. The dynamics dx dt of x expands the following way:

dx dt = f x (x 0 , u 0 , t)δx + f u (x 0 , u 0 , t)δu + 1 0 1 0 λf σσ (ρ(λ, µ), t)(δσ) 2 dλdµ (98)
Proof. By definition of x, we have x = x -x 0 -α x1 , so that:

dx dt = f (x 0 + δx, u 0 + δu, t) -LP [f, u 0 ](x 0 , t) + HP [f, u 0 ](x 0 , t) = f (x 0 + δx, u 0 + δu, t) -f (x 0 , u 0 , t)
The equation (98) follows by Taylor expansion of f (x, u, t) with integral remainder.

Proposition 24. L(x, u, t)) rewrites the following way:

L(x, u, t) = L(x 0 , u 0 , t) + H 0x (x + α x1 ) -p 0 dx dt + 1 0 1 0 λH σσ (ρ(λ, µ), p 0 , t)(δσ) 2 dλdµ (99)
Proof. We have the following identities: δx = x + α x1 and:

L x (x 0 , u 0 , t) = H 0x -p 0 f x (x 0 , u 0 , t)
Moreover, because of the stationary condition of the averaged problem, we have:

L u (x 0 , u 0 , t) = -p 0 f u (x 0 , u 0 , t)
Finally we change the integral remainder of the expansion of L(x, u, t) in equation ( 97) with:

L σσ (ρ(λ, µ), t) = H σσ (ρ(λ, µ), p 0 , t) -p 0 f σσ (ρ(λ, µ), t)

Thus equation (97) leads to:

L(x, u, t) = L(x 0 , u 0 , t + H 0x (x + α x1 ) -p 0 [f x (x 0 , u 0 , t) + f u (x 0 , u 0 , t) + 1 0 1 0 λf σσ (ρ(λ, µ), t)](δσ) 2 dλdµ + 1 0 1 0 λH σσ (ρ(λ, µ), p 0 , t)(δσ) 2 dλdµ
Inserting equation (98) in this equation proves proposition 24.

Proposition 25. The following equality holds:

T 0 -p 0 dx dt dt = - T 0 LP [H x , u 0 ](x 0 , t)xdt (100)
Proof. Let's make an integration by part:

T 0 -p 0 dx dt dt = -[p 0 x] T 0 + T 0 dp 0 dt xdt
This leads to the equation (100) because x(0) = 0, p 0 (T ) = 0 and:

dp 0 dt = -LP [H x , u 0 ](x 0 , t)
Now inserting the equation (100) in the equation (99) integrated between 0 and T leads to equation (63), which ends the proof of proposition 12.

E Proof of the lower bound of the third term (Proposition 13)

Proposition 26. The following inequality holds:

T 0 HP [H x , u 0 ](x 0 , p 0 , t)xdt ≥ -(2KT k r5 + k(1 + 2M ))α 2 + α T 0 p1 dr 1 dt dt (101)
Proof. By definition of r, we have:

x = r + αy 1 = (r -r 1 ) + (r 1 + αy 1 ) So that T 0 HP [H x , u 0 ](x 0 , p 0 , t)xdt = T 0 HP [H x , u 0 ](x 0 , p 0 , t)(r -r 1 )dt + T 0 HP [H x , u 0 ](x 0 , p 0 , t)(r 1 + αy 1 )dt
But because of equation (58), we have

T 0 HP [H x , u 0 ](x 0 , p 0 , t)(r -r 1 )dt ≤ T HP [H x , u 0 ] ∞ k r5 α 2 ≤ 2KT k r5 α 2
Moreover, by definition of α p1 , we have:

HP [H x , u 0 ](x 0 , p 0 , t) = -α d p1 dt (t)
So that the following inequality holds:

T 0 HP [H x , u 0 ](x 0 , p 0 , t)xdt ≥ -2KT k r5 α 2 -α T 0 d p1 dt ( x1 + αy 1 )dt
If we make an integration by part and use the fact that p1 (T ) = x1 (0) = y 1 (0) = 0, we get:

T 0 HP [H x , u 0 ](x 0 , p 0 , t)xdt ≥ -2KT k r5 α 2 + α T 0 p1 d x1 dt + α dy 1 dt dt
y 1 follows the dynamics (51) and p1 ∞ ≤ 1, so that:

α 2 T 0 p1 dy 1 dt dt ≤ kα 2 ( y 1 + x1 ∞ + v 1 ∞ ≤ k(2M + 1)α 2
So that the inequality (101) is proved.

Proposition 27. The dynamics of r 1 expands the following way:

dr 1 dt = [f 0x r 1 + f 0u v] + 1 0 1 0 λf σσ (x 0 + λµ(r 1 + α( x1 + y 1 )), u 0 + λµ(v + αv 1 ), t) r 1 + α( x1 + y 1 ) v + αv 1 2 dλdµ -α 2 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 ), u 0 + λµαv 1 , t) x1 + y 1 αv 1 2 dλdµ (102)
Proof. The dynamics of r r 1 (equation ( 57)) is the difference between the quantities f (r 1 + x 0 + α(x 1 + y 1 ), v + u 0 + αv 1 , t) and f (x 0 + α(x 1 + y 1 , u 0 + αv 1 , t).

Let's make the Taylor expansions of these quantities at (x 0 , u 0 ):

f (r 1 + x 0 + α(x 1 + y 1 ), v + u 0 + αv 1 , t) = f (x 0 , u 0 , t) +f 0x (r 1 + α(x 1 + y 1 )) + f 0u (v + αv 1 ) + 1 0 1 0 λf σσ (x 0 + λµ(r 1 + α( x1 + y 1 )), u 0 + λµ(v + αv 1 ), t) r 1 + α( x1 + y 1 ) v + αv 1 2 dλdµ and f (x 0 + α(x 1 + y 1 ), u 0 + αv 1 , t) = f (x 0 , u 0 , t) +α[f 0x (x 1 + y 1 )) + f 0u v 1 ] +α 2 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 ), u 0 + λµαv 1 , t) x1 + y 1 v 1 2 dλdµ
The simplifications between the two expansions while we take their differences gives the dynamics of r 1 (102).

The consequence of that dynamics expansion is that:

α T 0 p1 dr 1 dt dt = α T 0 p1 [f 0x r 1 + f 0u v]dt +α T 0 p1 1 0 1 0 λf σσ (x 0 + λµ(r 1 + α( x1 + y 1 )), u 0 + λµ(v + αv 1 ), t) r 1 + α( x1 + y 1 ) v + αv 1 2 dλdµ dt -α 3 T 0 p1 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 ), u 0 + λµαv 1 , t) x1 + y 1 v 1 2 dλdµ dt (103)
Proposition 28. Let's define the terms:

R 1 = α T 0 p1 1 0 1 0 λf σσ (x 0 + λµ(r 1 + α( x1 + y 1 )), u 0 + λµ(v + αv 1 ), t) r 1 + α( x1 + y 1 ) v + αv 1 2 dλdµ dt and R 2 = α 3 T 0 p1 1 0 1 0 λf σσ (x 0 + λµα( x1 + y 1 ), u 0 + λµαv 1 , t) x1 + y 1 v 1 2 dλdµ dt
Then the following inequalities hold:

|R 1 | ≤ 3k 2 [(k r3 T + k v1 ]αz 2 + [(k r4 T + k v2 + T (1 + 2M ) 2 ]α 3 (104) and |R 2 | ≤ kT 2 (1 + 2M ) 2 α 3 (105)
Proof. Equation (40) about the upper bound of f σσ leads to:

f σσ (x 0 + λµα( x1 + y 1 ), u 0 + λµαv 1 , t) r 1 + α( x1 + y 1 v + αv 1 2 ≤ 3k[|r 1 + α( x1 + y 1 )| + |v + αv 1 |] 2 ≤ 3k[|r 1 | + |v| + α(| x1 + y 1 | + |v 1 |)] 2 ≤ 3k[r 2 1 + v 2 + α(| x1 + y 1 | + |v 1 |) 2 ]
Thus, together with the fact that p1 ≤ 1, we have:

|R 1 | ≤ 3kα 2 [ r 1 2 2 + v 2 2 + α 2 T (2M + 1) 2 ] |R 1 | ≤ 3kα 2 [ r 1 2 ∞ + v 2 2 + α 2 T (2M + 1) 2 ]
Equations ( 59) and (56) then lead to equation (104). Equation (40) about the upper bound of f σσ leads also to: Proof. Let's use the fact that:

δσ = r + α x1 v + α y 1 v 1
to expand H σσ (ρ(λ, µ), p 0 , t)(δσ) 2 :

H σσ (ρ(λ, µ), p 0 , t)(δσ)

≤ 2M K 2r 2 + 2α 2 + M 2 α 2 3(r 2 + v 2 + α 2 ) ≤ 2M K √ 3 (2r 2 + v 2 + (M 2 + 2)α 2 )
Thus, with a triple integration after multiplication by λ, we get:

|J 41 | ≤ 4 √ 3KM α T r 2 ∞ + T 2 v 2 2 + ((M 2 + 2)α 2
Then, using equations (55) and (56) lead to equation (113).

Applying propositions 29 and 30 lead to: That proves proposition 14 because H 0σσ = H σσ (w 0 , t) does not depend on λ and µ.

6. 3 . 3 Proposition 18 .

 3318 Upper Bound of the third term of the development of J(u 0 ) (69) The following inequality holds:

  Including equation (72) and the fact that p1 ∞ ≤ 1 into equation (72) proves equation (71) 6.3.4 Upper Bound of the fourth term of the development of J(u 0 ) (69) Proposition 19. The following inequality holds:

C. 1 Proof of proposition 10 C. 1 . 1

 11011 Upper bound for r Proposition 20. The dynamics of r is the following:

f 1 v + αv 1 2 ≤

 12 σσ (x 0 + λµα( x1 + y 1 ), u 0 + λµαv 1 , t) r 1 + α( x1 + y k(| x1 + y 1 | + |v 1 |) 2 ≤ k(1 + 2M ) 2That proves equation (105) by triple integration and multiplication by α 3 .The equations (104) and (105) included in equation (105) (103) lead to:0x r 1 + f 0u v]dt -3k 2 (k r2 T + k v1 )M az 2 -k 2 [3(k r4 T + k v2 ) + 4T (1 + 2M ) 2 ]α 3Introducing that equation into equation (101) leads to equation (64) and thus the proposition 13.

F0λ [y 1 , v 1 ]

 11 Proof of the lower bound of the fourth term (Proposition 14) Proposition 29. The following inequality holds:λH σσ (ρ(λ, µ), p 0 , t)(δσ) 2 dλdµdt ≥ -2KT M 2 α 2 + βz 2 H σσ (ρ(λ, µ), p 0 t) r + α x1 v dλdµdt

0 2 + β - 2 √λ

 22 λH σσ (ρ(λ, µ), p 0 , t)(δσ) 2 dλdµdt ≥ -2KM T M + √ 3(2T k r2 + k v2 + T (M 2 + 2))α α 3KM (2T k r1 + k v1 )α z 2 [y 1 , v 1 ] H σσ (w 0 , t) r + α x1 v dλdµdt

  2 = H σσ (ρ(λ, µ), p 0 , t) r + α x1 vThe second term is easily upper bounded in absolute value (equation (41) about the upper bound of H σσ ):≤ K(|y 1 | + |v 1 |) 2 ≤ 4KM 2 -1 uu H ux ](ρ(λ, µ), p 0 , t)(r + α x1 to expand the terms in v and v 2 : 2(r + α x1 )H xu (ρ(λ, µ), p 0 , t)v = 2(r + α x1 )H xu (ρ(λ, µ), p 0 , t)Z[λ, µ] -2[H ux H -1 uu H xu ](ρ(λ, µ), p 0 , t)(r + α x1 ) 2(110)andH uu (ρ(λ, µ), p 0 , t)v 2 = H uu (ρ(λ, µ), p 0 , t)Z[λ, µ] 2 +[H ux H -1 uu H xu ](ρ(λ, µ), p 0 , t)(r + α x1 ) 2 (111) -2(r + α x1 )H xu (ρ(λ, µ), p 0 , t)Z[λ, µ]Because of assumption 3 on the convexity, this proves that:H σσ (ρ(λ, µ), p 0 , t) r + α x1 v , v 1 ] [H σσ (ρ(λ, µ), p 0 t) -H σσ (w 0 , t)] Proof. Let's expand [y 1 , v 1 ] [H σσ (ρ(λ, µ), p 0 t)H σσ (w 0 , t)] r + α x1 v into its coordinates: [y 1 , v 1 ] [H σσ (ρ(λ, µ), p 0 t)H σσ (w 0 , t)] r + α x1 v = y 1 [H xx (ρ(λ, µ), p 0 t) -H xx (w 0 , t)](r + α x1 )But the second derivatives of H(x, u, p 0 , t) are Lipschitz in (x, u) of Lipschitz constant K because the third derivatives of H(x, u, p 0 , t) are bounded by K. Thus with the definition of ρ(λ, µ) and w 0 , equation (114) leads to:|[y 1 , v 1 ] [H σσ (ρ(λ, µ), p 0 t)H σσ (w 0 , t)] r + α x1 v ≤ λµK (r + α x1 ) 2 + α 2 v2 1 (|y 1 ||r + α x1 | + |y 1 ||v| + |v||r + α x1 | + |v 1 ||v| ≤ K (r + α) 2 + α 2 M 2 (2M (|r| + α) + 2M |v|)

	2 dλdµdt ≤ 2KT M 2 + α 2 H σσ (ρ(λ, µ), p 0 , t) +2α[y 1 , v 1 ]H σσ (ρ(λ, µ), p 0 , t) r + α x1 v H σσ (ρ(λ, µ), p 0 , t) y 1 v 1 2 T 0 1 0 1 0 H σσ (ρ(λ, µ), p 0 , t) y 1 v 1 2 To upper bound the first term, let's expand it in its components: so that: H σσ (ρ(λ, µ), p 0 , t) r + α x1 v ≥ βZ[λ, µ] 2 so that: T 0 1 0 1 0 H σσ (ρ(λ, µ), p 0 , t) r + α x1 v 2 dλdµdt ≥ βz 2 Including equations (113) and (108) into equation (107) leads to equation (106). Proposition 30. Let's denote: J 41 = 2α T 0 1 0 1 0 λ [y 1 r + α x1 v dλdµdt y 1 v 1 Then the following bound holds: Introducing 2 |J 41 | ≤ 4 √ 3KM α (T k r1 + 1 2 k v1 z 2 + T k r2 + 1 2 k v2 + T 2 (M 2 + 2) α 2	2	(107) (108) (112) (113)

2

= H xx (ρ(λ, µ), p 0 , t)(r

+ α x1 ) 2 + 2(r + α x1 )H xu (ρ(λ, µ), p 0 , t)v +H uu (ρ(λ, µ), p 0 , t)v 2 (109) Now let's use the definition of Z[λ, µ](t) = v + [H equations (

110

) and (111) in equation (

109

) leads to:

H σσ (ρ(λ, µ), p 0 , t) r + α x1 v 2 = [H xx -H ux H -1 uu H xu ](ρ(λ, µ), p 0 , t)(r + α x1 ) 2 +H uu Z[λ, µ] 2 +y 1 [H xu (ρ(λ, µ), p 0 t) -H xu (w 0 , t)]v (114) +v 1 [H ux (ρ(λ, µ), p 0 t) -H ux (w 0 , t)](r + α x1 )

+v 1 [H uu (ρ(λ, µ), p 0 t) -H uu (w 0 , t)]v

The integral of HP is in α, so we consider that the derivative of HP is formally in 1 α

G Proof of Bound in absolute value for the sum of the integral terms of the right hand sides of equations ( 64) and (65)

Thanks to the fact that r + α x1 = r 1 + ((r -r 1 ) + α x1 ), we have R = R 3 + R 4 , where R 3 and R 4 are defined as:

Proposition 31. The following inequality holds:

Now let's use the costate dynamics (54) of the auxiliary problem, together with its stationarity condition (53). Then we get:

An integration by parts, together with the fact that r 1 (0) = q 1 (T ) = 0 leads to:

But r 1 follows the dynamics (57), so that a Taylor expansion of f (r 1 + x 0 + α( x1 + y 1 ), v + u 0 + αv 1 , t) at (x 0 + α( x1 + y 1 ), u 0 + αv 1 ) leads to:

so that, thanks to equation (40) about the bound of f σσ :

An integration on [0, T ] and a multiplication by α leads to

That equation, together with the bounds on r 1 (59) and v (56) lead to equation (115).

Proposition 32. The following inequality holds:

Thus, thanks to the bound on r -r 1 (58), its absolute value can be bounded:

An integration between 0 and T and a multiplication by α lead to equation (116).

Proof of equation (66) This is a consequence of equations ( 115) and ( 116), together with the fact that R = R 3 + R 4 .


Réda Alami
email: reda1.alami@orange.com

Odalric Maillard
email: odalric.maillard@inria.fr

Raphael Féraud
email: raphael.feraud@orange.com

Memory Bandits: a Bayesian approach for the Switching Bandit Problem

. In this paper, we propose another extension of this approach based on a Bayesian aggregation framework. Experiments provide some evidences that in practice, the proposed algorithm compares favorably with the previous version of Thompson Sampling for the Switching Multi-Armed Bandit Problem, while it outperforms clearly other algorithms of the state-of-the-art.

Introduction

We consider the non-stationary multi-armed bandit problem with a set K of K independent arms. At each round t, the agent chooses an action k t K and observes a reward x kt . In a stationary environment, the agent has to explore to find the best arm and to exploit it to maximize his gain. Efficient algorithms [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF][START_REF] Cappé | Kullback-leibler upper confidence bounds for optimal sequential allocation[END_REF][START_REF] Kaufmann | Thompson sampling: An asymptotically optimal finite-time analysis[END_REF] have been proposed for handling the so-called exploration-exploitation dilemma. In an evolving environment, the best arm can change during time and hence the agent has to explore more. The adversarial bandits handles non-stationary environments by considering that a sequence of deterministic rewards is chosen in advance by an oblivious adversary. Rate optimal algorithms have been proposed to find the best arm or the best sequence of arms of the run in [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF][START_REF] Neu | Explore no more: Improved high-probability regret bounds for non-stochastic bandits[END_REF][START_REF] Allesiardo | The non-stationary stochastic multi-armed bandit problem[END_REF]. Another approach for handling non-stationary environment is to consider that the rewards are generated by an unknown stochastic process that evolves during time. In the switching bandit problem [START_REF] Hartland | Multiarmed bandit, dynamic environments and meta-bandits[END_REF][START_REF] Garivier | On upper-confidence bound policies for switching bandit problems[END_REF][START_REF] Allesiardo | Exp3 with drift detection for the switching bandit problem[END_REF], the mean rewards of arms change abruptly. In comparison to the adversarial approach, the advantage of non-stationary stochastic approach is that with some mild assumptions, stochastic algorithms, which are more efficient in practice than adversarial algorithms, can be used. For practical and theoretical reasons, the recent years have seen an increasingly interest for the oldest bandit algorithm, the Thompson Sampling [START_REF] William | On the likelihood that one unknown probability exceeds another in view of the evidence of two samples[END_REF]. It exhibits excellent results in practice [START_REF] Chapelle | An empirical evaluation of thompson sampling[END_REF], while it is asymptotically optimal [START_REF] Kaufmann | Thompson sampling: An asymptotically optimal finite-time analysis[END_REF]. In [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF], the authors propose an adaptation of the Thompson Sampling to the switching bandit problem. To be consistent with the Bayesian approach, rather than using a frequentist drift detector used in [START_REF] Hartland | Multiarmed bandit, dynamic environments and meta-bandits[END_REF][START_REF] Allesiardo | Exp3 with drift detection for the switching bandit problem[END_REF], the authors use a Bayesian online change point detection [START_REF] Prescott | Bayesian online changepoint detection[END_REF] combined with the Thompson Sampling algorithm. In this paper, we propose a similar approach of [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF] which is based on a Bayesian aggregation of a growing number of experts seen as learners. Our approach compares favorably with the one of [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF].

Problem formulation

Let us consider an agent facing a non-stationary multi-armed bandit problem with a set K t1, ..., Ku of K independent arms. At each round t rr1, T ss, the agent chooses to observe one of the K possible actions. When playing the arm k t at time t, a reward x kt is received, where x kt B pµ kt,t q is a random variable drawn from a Bernoulli distribution of expectation µ k,t . Let µ t max kK tµ k,t u denotes the best expected reward at round t, k t arg max kK tµ k,t u the best arm at round t, k t the action chosen by the decision-maker at time t and x kt the reward obtained at the same time.

31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA.

Changes in the Bernoulli distributions expectations It should be noted that µ k,t , i.e. the reward mean of arm k at time t changes over time according to a global abrupt switching model parametrized with an unknown hazard function hptq r0, 1s assumed to be a constant hptq ρ such that:

µ k,t 4 µ k,t¡1
with probability 1 ¡ ρ µ new Up0, 1q with probability ρ (1)

When the behavior's environment is modeled by equation ( 1) for all k K, the problem setting is called a Global Switching Multi-Armed Bandit (GS-MAB), i.e. when a switch happens all arms change their expected rewards. Where changes occur independently for each arm k (i.e. arms change points are independent from an arm to another), the problem setting is called a Per-arm Switching Multi-Armed Bandit. For the sake of clarity, in the following we will focus on GS-MAB.

Sequence of change points It should be noted that for each GS-MAB, it exists a non-decreasing change points sequence of length Υ T denoted by pτ κ q κrr1,Υ T 1ss N Υ T 1 where:

4 d κ rr1, Υ T ss, d t T κ rrτ κ 1, τ κ 1 ss, d k K, µ k,t µ k,rκs τ 1 1 τ 2 ... τ Υ T 1 T In this case, µ rκs max k 2 µ k,rκs
@ denotes the highest expected reward at epoch T κ . Pseudo Cumulative Regret for the switching environment In a switching environment, the pseudo cumulative regret RpT q up to time T is defined as the expected difference between the rewards obtained by our policy and those received by the oracle which always plays the best arm k rκs at each epoch T κ such as:

R pTq T ţ1 µ t ¡ E T ţ1 x kt ' Υ T κ1 £ |T κ | µ rκs ¡ E τκ 1 ţτκ 1 x kt '
3 Global Switching Thompson Sampling with Bayesian Aggregation

The Thompson Sampling algorithm

Unlike optimistic algorithms belonging to the UCB family [START_REF] Auer | Finite-time analysis of the multiarmed bandit problem[END_REF], which are often based on confidence intervals, the Thompson Sampling deals with Bayesian tools by assuming a Beta prior distribution π k,t1 Betapα 0 , β 0 q on each arm for some α 0 , β 0 ¡ 0. Based on the rewards observed x kt , the posterior distribution π k,t is updated such as: π k,t Beta pα k,t #preward 1q α 0 , β k,t #preward 0q β 0 q. At each time, the agent takes a sample θ k,t from each π k,t and then plays the arm k t arg max k θ k,t . Formally, by denoting D t¡1 t¡1 i1 x i the history of past rewards we write: θ t pθ 1,t , ..., θ K,t q Ppθ t |D t¡1 q ± K k1 π k,t . Recently, the Thompson Sampling has been shown to be asymptotically optimal [START_REF] Kaufmann | Thompson sampling: An asymptotically optimal finite-time analysis[END_REF], i.e. the expectation of the pseudo-cumulative regret reaches the Lai and Robbins lower bound on regret in the stochastic Bernoulli multi-armed bandit setting [START_REF] Leung | Asymptotically efficient adaptive allocation rules[END_REF].

Decision making based on Bayesian aggregation

Best achievable performance: the Thompson Sampling oracle Let T S denotes the oracle that knows exactly the change points τ κ . It simply restarts a Thompson Sampling at these change points.

Assume that Υ T is the overall number of change points observed until T , then T S runs successively Υ T Thompson Sampling processes starting at τ κ 1 and ending at τ κ 1 .

Notion of expert Let t N and i rr1, tss. An expert f i,t is a Thompson Sampling procedure which has started at time i. The expert f i,t observes exactly t ¡ i rewards from the environment.

Expert Aggregation Like [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF], to characterize the occurrence of changes, we use the expert f i,t as an index to access in the memory the parameters of the model created at time i. The computation of Ppθ t |D t¡1 q is done by taking into account the distribution w i,t P pf i,t |D t¡1 q of the expert f i,t such as:

Ppθ t |D t¡1 q t i1 Ppθ t |D t¡1 , f i,t qPpf i,t |D t¡1 q (2)
Then, unlike the sampling procedure used in [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF], we build the index θ k,t of arm k at time t by launching a Bayesian aggregation of a growing number of experts. The estimation of the expert distribution is done recursively according to the work of [START_REF] Prescott | Bayesian online changepoint detection[END_REF] where:

P pf i,t |D t¡1 q loooooomoooooon Expert distribution at t ∝ t¡1 i1 change point prior hkkkkkkikkkkkkj Ppf i,t |f i,t¡1 q Ppx kt |f i,t¡1 , D t¡2 q loooooooooomoooooooooon Instantaneous gain Expert distribution at t¡1 hkkkkkkkikkkkkkkj Ppf i,t¡1 |D t¡2 q (3) 2
The change point prior Ppf i,t |f i,t¡1 q is naturally computed following equation (1):

Ppf i,t |f i,t¡1 q p1 ¡ ρq 1 pi tq ρ1 pt 1q

Thus, the inference model takes the following form (Up to a normalization factor):

4

Growth probability:

Ppf i,t |D t¡1 q ∝ p1 ¡ ρq ¤ Ppx kt |f i,t¡1 , D t¡2 q ¤ Ppf i,t¡1 |D t¡2 q change point probability: Ppf t,t |D t¡1 q ∝ ρ °t¡1 i1 Ppx kt |f i,t¡1 , D t¡2 q ¤ Ppf i,t¡1 |D t¡2 q (5) It should be noted that Ppx kt |f i,t¡1 , D t¡2 q is a Bernoulli distribution of expectation α k t ,i,t¡1
α k t ,i,t¡1 β k t ,i,t¡1 , where α kt,i,t¡1 and β kt,i,t¡1 are the hyper-parameters of the arm k t learned by the expert f i,t¡1 . Let l i,t¡1 denotes the instantaneous logarithmic loss associated to the forecaster f i,t¡1 such as: Ppx kt |f i,t¡1 , D t¡2 q exp p¡l i,t¡1 q, then Algorithm 1 provides us an index prediction of each arm k at time t based on a Bayesian aggregation of the available experts.

Algorithm 1 Bayesian Aggregation with a growing number of experts for t 2, ... do -1-Update: di t1, ..., t ¡ 1u w i,t Ð p1 ¡ ρq w i,t¡1 exp p¡l i,t¡1 q -2-Create new expert starting at t: w t,t Ð ρ °t¡1 i1 w i,t¡1 exp p¡l i,t¡1 q -3-Predict arm index:

d k K θ k,t Ð °t¡1 i1 θ k,i,t wi,t °t¡1 i1 wi,t where: θ k,i,t Beta pα k,i,t , β k,i,t q
Finally, by plugging the Bayesian aggregation into the formalism of [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF], we get the Global Switching Thompson Sampling with Bayesian Aggregation (Global-STS-BA), described in Algorithm 2.

Algorithm 2 Global Switching Thompson Sampling with Bayesian Aggregation

1: procedure GLOBAL-STS-BA(K, T, α 0 , β 0 , ρ)

2: t Ð 1, w 1,t Ð 1, and dk K α k,1,t Ð α 0 , β k,1,t Ð β 0 Initializations 3:
for t ¤ T do Interaction with environment 4:

k t Ð CHOOSEARM ptwu t , tαu t , tβu t q 5:
x kt Ð PLAYARMpk t q Bernoulli trial 6:

twu t 1 Ð UPDATEEXPERTWEIGHT ptwu t , tαu kt,t , tβu kt,t , x kt , ρq

l i,t Ð ¡x kt log ¡ α k t ,i,t α k t ,i,t β k t ,i,t © ¡ p1 ¡ x kt q log ¡ β k t ,i,t α k t ,i,t β k t ,i,t © d i rr1,
w i,t 1 Ð p1 ¡ ρq w i,t exp p¡l i,t q d i rr1, tss
Increasing the size of expert f i,t

β kt,i,t 1 Ð β kt,i,t 1 px kt 0q d i rr1, tss 19: α k,t 1,t 1 Ð α 0 , β k,t 1,t 1 Ð β 0 d k K
Initializing new expert 20:

return tαu t 1 , tβu t 1
Global-STS-BA In the global switching setting, when a switch occurs, all arms change their expected pay-off at the same time. It should be noted that the data from all plays collaborate in the posterior of the expert distribution estimation (UPDATEEXPERTWEIGHT). The experts tell us how much previous observed data can be used in the arm indexes prediction, i.e. when a switch occurs all past data have not to be taken into account in the arm characterizations because of their obsolete information. The change point detection concept is based on a tracking of the optimal expert (see Appendix B.1). Indeed, the total mass of the expert distribution w i,t tends to focus around the optimal expert f τκ,t i.e. the expert starting at the most recent change point τ κ and corresponds to the most appropriate characterization of the environment. The Bayesian aggregation exploits this concentration to highlight the contribution of the most appropriate experts (starting around the most recent change point τ κ ) in the arm indexes prediction. The concentration around f τκ,t is possible thanks to the inference model of equation [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF]. In fact, when a change occurs the instantaneous gain P px kt |f i,t¡1 , D t¡2 q of all experts starting before the change point suddenly fall down because of their wrong estimation of the environment, giving the advantage to the experts newly created while annihilating the former ones (see Appendix B.1). At this point, we deviate from [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF] and instead of sampling the expert distribution and then sampling the arms, we index each arm k by launching an overall Bayesian aggregation of samples taken from the posterior distribution of arm k related to the hyper-parameters tαu k,t , tβu k,t (CHOOSEARM). Finally, the agent chooses to pull the arm with highest index. This process allows us to avoid the sampling noise induced by Global-STS and by this way the model chosen at time t tends to better fit the unknown environment (figure 1).

Experiments

In all the experiments, we consider a GS-MAB of three arms observing three change points occurring at each 1000 rounds. Experiments are run 100 times. The parameters of the state-of-the art algorithms are chosen to be experimentally optimal. Exp3, Exp3P and Exp3S [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF] are launched with an exploration rate (γ 5%). We run Exp3R [START_REF] Allesiardo | Exp3 with drift detection for the switching bandit problem[END_REF], Rexp3 [START_REF] Besbes | Stochastic multi-armed-bandit problem with nonstationary rewards[END_REF] and SW-UCB [START_REF] Garivier | On upper-confidence bound policies for switching bandit problems[END_REF] respectively with H 1000, ∆ T 1000 and τ 500. Global-STS [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF] and Global-STS-BA are outperforming the well parametrized state of art non stationary MAB algorithms (figure 1).

Replacing the expert distribution sampling used in [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF] with the Bayesian Aggregation allows us to obtain performances challenging those of the Thompson Sampling Oracle (figure 1).

A Extension to the per-arm switching setting

A.1 Per-Arm Switching Thompson Sampling with Bayesian Aggregation (Per-Arm-STS-BA)

In the per-arm switching setting, since the changes occur independently from an arm to another, each arm k will have its own expert distribution denoted by twu t k . To estimate twu t k at each time step t, we need to build some recursive message-passing algorithm. The main difference between the global and the per-arm version of the message-passing is that the instantaneous gain will be used in the update of the expert distribution associated to the pulled arm k t . For the other arms not pulled, the instantaneous gain won't intervene in the expert distribution [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF]. More formally, if we denote by f k,i,t the expert associated to the arm k at time t which has been introduced at time i, then the expert distributions are updated as follow:

Ppf k,i,t |D t¡1 q ∝ 5 °t¡1 i1 P pf k,i,t |f k,i,t¡1 q Ppx kt |f k,i,t¡1 , D t¡2 qPpf k,i,t¡1 |D t¡2 q if k is pulled °t¡1 i1 P pf k,i,t |f k,i,t¡1 q Ppf k,i,t¡1 |D t¡2 q
otherwise Then, we model the changes of the arm k with a constant switching rate γ k such that:

P pf k,i,t |f k,i,t¡1 q 4 1 ¡ ρ k if i t ρ k if i t
Then, following the inference model used in the global switch, we deduce the following expert distributions update rules (Up to a normalization factor):

Growth probability:

P pf k,i,t q ∝ 4 p1 ¡ ρ k qPpx kt |f k,i,t¡1 qPpf k,i,t¡1 q if k is pulled p1 ¡ ρ k qPpf k,i,t¡1 q otherwise change point probability: Ppf k,t,t q ∝ 5 ρ k °t¡1 i1 Ppx kt |f k,i,t¡1 qPpf k,i,t¡1 q if k is pulled ρ k °t¡1 i1 Ppf k,i,t¡1 q otherwise Notice that: Ppx kt |f k,i,t¡1 q still a Bernoulli distribution of expectation α k t ,i,t¡1
α k t ,i,t¡1 β k t ,i,t¡1 , where α kt,i,t¡1 , β kt,i,t¡1 are the hyper-parameters of arm k t at time t ¡ 1 learned by the expert f kt,i,t¡1

which is associated to the change point model of arm k t . We define:

Ppx kt |f kt,i,t¡1 q exp p¡l i,t¡1 q where: l i,t¡1 denotes the instantaneous logarithmic loss of the pulled arm k t associated to the forecaster f kt,i,t¡1 . For convenience, we write: w k,i,t Ppf k,i,t |D t¡1 q. Then, we naturally extend the Bayesian aggregation used in the global switching setting to the per-arm Bayesian Aggregation (Algorithm 3).

Algorithm 3 Per-arm Bayesian Aggregation

for t 2, ... do -1-Update: di t1, ..., t ¡ 1u d k K w k,i,t Ð p1 ¡ ρ k q w k,i,t¡1 exp p¡l i,t¡1 q 1pkktq -2-Create new expert starting at t: d k K w k,t,t Ð ρ k °t¡1 i1 w k,i,t¡1 exp p¡l i,t¡1 q 1pkktq -3-Predict arm index: d k K θ k,t Ð °t¡1 i0 θ k,i,t w k,i,t °t¡1 i0 w k,i,t where: θ k,i,t Beta pα k,i,t , β k,i,t q
Then, by plugging the per-arm Bayesian Aggregation into the formalism of the Global-STS-BA we get the Per-arm Switching Thompson Sampling with Bayesian Aggregation (Per-arm-STS-BA) described in Algorithm 4.

A.2 Experiments

In all the experiments, we consider a Per-Arm Switch MAB of three arms observing several change points. Experiments are run 100 times. Per-Arm-STS [START_REF] Mellor | Thompson sampling in switching environments with bayesian online change point detection[END_REF] and Per-Arm-STS-BA are outperforming the well parametrized state of art non stationary MAB algorithms (figure A.2). Exp3, Exp3P and Exp3S [START_REF] Auer | The nonstochastic multiarmed bandit problem[END_REF] are launched with an exploration rate (γ 5%). We run Exp3R [START_REF] Allesiardo | Exp3 with drift detection for the switching bandit problem[END_REF], Rexp3 [START_REF] Besbes | Stochastic multi-armed-bandit problem with nonstationary rewards[END_REF] and SW-UCB [START_REF] Garivier | On upper-confidence bound policies for switching bandit problems[END_REF] respectively with H 600, ∆ T 600 and τ 600. Optimal expert Let t ¡ 0 and let τ rts denotes the most recent change point before time t. At each time t, a set of exactly t experts is available (each expert has started at i rr1, tss). The optimal expert is the one which has started exactly at t τ rts . It is the expert which gives the best description of the environment.

Tracking the optimal expert Before the decision step, Global-STS-BA characterizes each arm k by aggregating all the contributions of the available experts. To do well, Global-STS-BA needs to highlight the contribution of the optimal expert because of its optimal description of the environment. This task is called the track of the optimal expert.

The estimation of the expert distribution (expert weight) is built according to the instantaneous gain of each expert. This gain takes into account the history of past rewards D t¡1 . When a switch occurs, the instantaneous gains of all experts fall down because of their obsolete estimation giving the advantage to the expert newly created which will become the optimal expert during the next rounds.

Let us consider a GS-MAB of three arms observing two change points at time t 25 and t 49.

Discussions

• Just before the change point, the optimal expert (index = 1 for the first switch and index = 25 for the second one) is given the highest weight. The lower the index, the higher the expert gain. This behavior is expected as much as the higher the number of the observations the higher the expert gain.

• When the switch occurs, all the experts created before the change point see their gains dropping sharply because of the abrupt change of the environment i.e. the instantaneous gain doesn't match anymore with the current environment, except the expert newly created which is given a gain based on the prior of a Thompson Sampling ¡ α0 α0 β0 © .

• Just after the change point, the new optimal expert see its weight starting to grow up while the weight of the previous optimal expert starts to fall down. • Largely after the change point, the optimal expert is given the highest weight, because of its well fit of the current environment. This is corresponding to the expected behavior of the Global-STS-BA: the optimal expert is well tracked.

B.3 Behavior of Global-STS-BA with respect to the switching rate ρ

We should also notice that the inference model of the experts needs the knowledge of the true switching rate (ρ true ). In real life, this value is unknow to the agent. [START_REF] Turner | Adaptive sequential Bayesian change point detection[END_REF] has proposed to learn the switching rate from the data via a gradient descent. This method appears to not perform particularly well if the switching rate has to be adapted at every time step. w kt,i,t 1 Ð p1 ¡ ρq w kt,i,t exp p¡l i,t q d i rr1, tss Increasing the size of expert f kt,i,t 16:

w kt,t 1,t 1 Ð ρ °i w kt,i,t exp p¡l i,t q Creating new expert for arm k t starting at t 1

  tss

Figure 1 :

 1 Figure 1: Overall comparison with the non-stationary state of art and the TS Oracle.

  Replacing the expert distribution sampling used in STS with the Bayesian Aggregation allows us to obtain performances challenging those of the Thompson Sampling Oracle (figure A.2).

Figure 2 :

 2 Figure 2: Overall comparison with the state of art and the oracle.

Figure B. 1

 1 Figure B.1 shows the behavior of the expert gain and the expert weight when a switch occurs. It should be noted that the experts are indexed by their starting time i.e. the most recent expert is the one with highest index and inversely.

  (a) Tracking the optimal expert during the first switch (b) Tracking the optimal expert during the second switch

B. 2

 2 Behavior of Global-STS-BA with respect to the memory size Regarding the implementation of the Global-STS-BA, we should notice that memory space and time requirements of the experts inference model grow linearly at each time step because of the creation of a new expert. This makes the size of the support set of the expert distribution twu t increasing by one.Thus, for computational reasons, we propose to restrict the number of experts by fixing a maximum memory size M . So, at each round, after launching the inference model, we delete the worst expert f min t arg min i w i,t . In all the experiments, we consider a GS-MAB of three arms. Changes occur at each 1000 rounds. We variate the memory size from M 15 to M 3000. Experiments are run 100 times.

Figure 3 :

 3 Figure 3: Behavior with respect to the switching rate ρ. For a not very poor memory size value (M ¥ 50), the performances of the Global-STS-BA remains stable.

  Figure B.3 shows the behavior of the Global-STS-BA for different values of the switching rate.Discussions The performances of the Global-STS-BA remains stable even if the switching rate is quite far from the true one.Algorithm 4 Per-Arm Switching Thompson Sampling with Bayesian Aggregation1: procedure PER-ARM-STS-BA(K, T, α 0 , β 0 , ρ)

  twu t 1 16: procedure UPDATEARMMODEL(tαu t , tβu t , x kt , k t ) 17:α kt,i,t 1 Ð α kt,i,t 1 px kt 1q d i rr1, tss

	18:

14:

w t 1,t 1 Ð ρ °i w i,t exp p¡l i,t q

Creating new expert starting at t 1

15:

return

  Ð 1, and dk K w k,1,t Ð 1, α k,1,t Ð α 0 , β k,1,t Ð β 0 Initializations t Ð CHOOSEARM ptwu t , tαu t , tβu t q twu t 1 Ð UPDATEEXPERTWEIGHT ptwu t , tαu kt,t , tβu kt,t , x kt , ρq 7:tαu t 1 , tβu t 1 Ð UPDATEARMMODELptαu t , tβu t , x kt , k t q k K d i rr1, tss θ k,i,t Ð Beta pα k,i,t , β k,i,t q 11:return arg max k °irr1,tss w k,i,t °jr r1,tss w k,j,t θ k,i,t procedure UPDATEEXPERTWEIGHT(twu t , tαu kt,t , tβu kt,t , x kt , ρ)

	2: 3:	for t ¤ T do					Interaction with environment
	4: 5:	x kt Ð PLAYARMpk t q				Bernoulli trial
	6:
								Bayesian aggregation
	12:
	13: 14:	l i,t Ð ¡x kt log	¡	α k t ,i,t α k t ,i,t β k t ,i,t	© ¡ p1 ¡ x kt q log	¡	β k t ,i,t α k t ,i,t β k t ,i,t	© d i rr1, tss
	15:

t k 8: 9: procedure CHOOSEARM(twu t , tαu t , tβu t ) 10: d

  k,i,t 1 Ð p1 ¡ ρq w k,i,t d i rr1, tss and d k $ k t Increasing the size of expert f k,i,t 18:w t 1,t 1,k Ð ρ °i w k,i,t d k $ k t Creatingnew expert for arm k starting at t 1 procedure UPDATEARMMODEL(tαu t , tβu t , x kt , k t ) 22: α kt,i,t 1 Ð α kt,i,t 1 px kt 1q d i rr1, tss 23: β kt,i,t 1 Ð β kt,i,t 1 px kt 0q d i rr1, tss 24: α k,t 1,t 1 Ð α 0 , β k,t 1,t 1 Ð β 0 d k K

	17:
	19:	return twu t 1
	20:
	21: Initializing new expert
	25:

w return tαu t 1 , tβu t 1

Acknowledgments

This work has been supported by the French Agence Nationale de la Recherche (ANR), under grant ANR-16-CE40-0002 (project BADASS).