P Malisani
email: paul.malisani@ifpen.fr

J Zhu

P Pognant-Gros
email: philippe.pognant-gros@ifpen.fr

Optimal charging scheduling of electric vehicles: the co-charging case

Keywords: Co-charging stations, electric vehicle, battery charging, optimal scheduling, Branch and Bound methods, mixed integer optimal control

A co-charging parking system is a charging system in which the chargers are automated, i.e., they can automatically connect to or disconnect from vehicles presented in the parking. This system makes it possible to have less chargers than the number of parking places. This paper studies an optimal scheduling problem aiming at charging all vehicles to required battery level while minimizing the total charging cost. We formulate this scheduling problem as a non-linear optimal control problem with mixed integer constraints and mixed discrete and continuous time. The problem is solved using a Branch and Bound Strategy coupled with Pontryagin Maximum Principle (PMP) and Interior Point Methods (IPMs) for optimal control. The presented method is tested on numerous configurations and shows promising results both in term of precision and in term of convergence rate.

I. INTRODUCTION

The increasing use of electric vehicles (EVs) is a major component of improving air quality in cities. To support this development actual car parks must be renewed with charging equipment. In many cases, this renewal consists in equipping a certain number of parking places with chargers. This solution can be expensive due to the large number of chargers to be deployed and lacks of flexibility in separating parking spaces for electric and thermal vehicles. In this context, so-called co-charging solutions can help to alleviate these drawbacks. In a co-charging system, chargers can move throughout the parking and plug-in and recharge chosen electric vehicles. This autonomous moving chargers allows to recharge a number of vehicles greater than the number of chargers (typically one charger for 3 to 5 parking places) by means of planning the sequence of plug-in/charging/unplug cycles of each charger. In the case we are interested in, the autonomous chargers can switch from a vehicle to another only at a regular time steps. In order to operate a co-charging system one has to schedule the charging of the EVs in order to ensure that they leave the parking lot with required state of charge while minimizing a user defined objective function, typically the charging cost. At present, this charging technology is still in the early planning stage and has not yet been tested. However, some manufacturers are starting to take an interest in it which motivates this contribution [START_REF]Initial contact[END_REF].

Optimal charging of EVs has been extensively studied for years. On the one hand, some of this literature focus on single EV optimal charging using different methods ranging from genetic algorithm to deep-learning-based methods [START_REF] Liu | Online energy management for multimode plug-in hybrid electric vehicles[END_REF]- [START_REF] Xiong | Reinforcement learningbased real-time power management for hybrid energy storage system in the plug-in hybrid electric vehicle[END_REF]. On the other hand, another part of the literature focus on EV fleet charging scheduling [START_REF] Morstyn | Conic optimization for electric vehicle station smart charging with battery voltage constraints[END_REF], [START_REF] Rui | A distributed charging strategy based on day ahead price model for pv-powered electric vehicle charging station[END_REF]. The problem we are interested in lies in between both cases and concerns the optimal charging of a few EVs, i.e. more than one vehicle but not enough to consider the statistical properties of the EVs (not enough to make a fleet). Indeed, the charging scheduling of a small group of EVs has also been addressed in recent works. In [START_REF] Zhang | Optimal management for parking-lot electric vehicle charging by two-stage approximate dynamic programming[END_REF], the authors studied a stochastic problem of optimal charging of plug-in electric vehicles at a parking-lot. In [START_REF] You | Optimal cooperative charging strategy for a smart charging station of electric vehicles[END_REF], [START_REF] Sabillon | A new methodology for the optimal charging coordination of electric vehicles considering vehicle-to-grid technology[END_REF] the optimal charging problem of electric vehicles that are able to exchange energy between them and with the electrical grid (V2G) has been addressed. However, these studies did not consider the case where the number of vehicles to be charged is larger than the number of chargers whereas it is a key feature regarding the co-charging case. In [START_REF] Zheng | Online distributed mpc-based optimal scheduling for ev charging stations in distribution systems[END_REF] the problem of the scheduling of several chargers distributed over a an electrical grid is addressed and the scheduling is handled using a heuristic prioritizing the vehicles with smaller departure times. This charging method does not allow to take advantage of a time varying price of electricity or does not allow the charging strategy to reduce energy losses when possible. In [START_REF] Salah | Simultaneous electric vehicles scheduling and optimal charging in the business context: case study[END_REF], the authors address the optimal assignment of tours to be processed of a mixed fleet of combustion and electric vehicles together with the optimal charging scheduling for the latter. The battery model used in the study is a discrete integrator model, which does not take into account the non-linearities in the battery dynamics. But, optimal control strategy of linear systems consists in bang-bang control [START_REF] Bryson | Applied optimal control[END_REF], that is to say always charging EVs at maximal power. This bangbang strategy does not allow to take advantage of the higher efficiency of a low power charging and hence yields an over consumption of energy.

The dependence of the charging efficiency on the state-ofcharge and on the charging power makes the total amount of energy required to charge a vehicle depends on the charging trajectory, i.e. the optimization variables. Indeed, the faster the vehicle is recharged, the greater the energy losses. Therefore, the problem we are interested in is not a combinatorial problem such as the bin-packing problem but must be considered as a mixed integer nonlinear optimal control problem. Several contributions in the literature focus on mixed integer optimal control problems referred as hybrid optimal control [START_REF] Sussman | A maximum principle for hybrid optimal control problem[END_REF], [START_REF] Buss | Nonlinear hybrid dynamical systems: Modeling, optimal control, and applications[END_REF] or mixed integer optimal control problems [START_REF] Sager | Reformulations and algorithms for the optimization of switching decisions in nonlinear optimal control[END_REF], [START_REF] Sager | The integer approximation error in mixed-integer optimal control[END_REF]. In [START_REF] Sussman | A maximum principle for hybrid optimal control problem[END_REF], [START_REF] Buss | Nonlinear hybrid dynamical systems: Modeling, optimal control, and applications[END_REF], the authors considered a pure continuous time problem, in which switching times of the integer control variable are free and must be computed. In [START_REF] Sager | Reformulations and algorithms for the optimization of switching decisions in nonlinear optimal control[END_REF], [START_REF] Sager | The integer approximation error in mixed-integer optimal control[END_REF], approximations methods such as sum-up rounding strategies are presented and successfully applied. These methods consist of a timediscretization of the mixed integer optimal control problem together with a continuous relaxation of the integer constraints followed by rounding them into integer solutions. These methods converge to the optimal solution as the switching time steps converge to 0. Unfortunately, in the co-charging case, these switching times are fixed and are much larger than the time scale of the battery dynamics. In this situation, there is no guarantee that a rounding strategy converges to an optimal solution or even a feasible one. Consequently, the co-charging problem is a Non-Linear Optimal Control Problem with Mixed Integer constraints and Mixed discrete and continuous Time (NLOCPMIMT). The integer part of the problem stems from the fact that only some of the vehicles can be charged at the same time (since there are more parking slots than chargers). The continuous time optimal control part stems from the fact that the charging of electrical batteries can be modeled as a non-linear dynamical system. Finally, the discrete time part stems from the fact that the chargers can move from a vehicle to another at certain time steps. To the knowledge of the authors, this problem has not been addressed in the literature.

In this paper, we tackle the NLOCPMIMT in its full complexity. More precisely, we consider a non-linear battery model calibrated on experimental data, arbitrary numbers of vehicles and chargers. Our main contribution is an algorithm that converges to the global optimum of the NLOCPMIMT. Our approach is as follows, we rely on a Branch and Bound approach, since the large switching time steps ensure a reasonable depth of the associated search-tree. This Branch and Bound method is combined with a continuous time interior point method for optimal control allowing one to compute the optimal solution of the Branch and Bound lower bounding procedure with high precision and small execution time. In addition, the presented methodology is still quite computationally expensive and additional work might be necessary for an online implementation. Nevertheless, the global convergence of the presented algorithm makes this methodology a meaningful benchmark against which other scheduling strategy can be compared. Such a benchmark is particularly useful to assess the performances of heuristics based methods.

This paper is organized as follows, in Section II the nomenclature is presented. In Section III, the optimal scheduling of vehicle charging problem is presented. In detail, in Section III-A the non-linear battery dynamics is presented, in Section III-B the NLOCPMIMT is presented. In Section IV, the Branch and Bound algorithm used to handle the integer variables is described. In Section V, the lower bounding part of the Branch and Bound algorithm is detailed. In Section VI the convergence of the Branch and Bound algorithm is investigated. A numerical study of an optimal vehicle charging scheduling is given in Section VII. Finally, conclusions and perspectives are discussed in Section IX.

II. NOMENCLATURE

• n x : number of vehicles to be charged • n c : maximum number of chargers • P c = {p 1 , . . . , p ncs }: set of maximal charging powers • n cs : number of groups of same maximum power chargers • G j : Group of chargers of same maximum power, j = 1, . . . , n

• k f = ⌊t f /∆T ⌋: maximum number of switching times • x 0,k : initial battery energy of vehicle k • x f,k : desired battery energy of vehicle k • p c,k ∈ C 2 (R, R): charging power of vehicle k • λ k ∈ L ∞ ([0, t f ], [0, 1]): charging rate of vehicle k • x k ∈ C 0 ([0, t f ], R): battery state of energy of vehicle k • p + c,k ∈ C 2 (R, R): maximal charging power of vehicle k • f c,k : R × R → R: dynamics of vehicle battery k • θ x ∈ R + :
M k = ⟨m 0 , . . . , m k ⟩, m i ∈ R nx×ncs , i = 0, . . . , k (1)
And we note the concatenation between a sequence M k and a matrix p ∈ R nx×ncs as follows

M k : p = ⟨m 0 , . . . , m k , p⟩ (2)
The concatenation of two sequences M k = ⟨m 0 , . . . , m k ⟩ and S l = ⟨s 0 , . . . , s l ⟩ writes as follows

M k : S l = ⟨m 0 , . . . , m k , s 0 , . . . , s l ⟩ (3)

III. PROBLEM PRESENTATION

The problem we are interested in consists in charging n x electric vehicles using n b autonomous charging stations that can automatically move from one parked vehicle to another with a regular time-step ∆T . The question of planing the order of vehicles charging arise naturally in this context. This paper provides an optimal control based algorithm to compute such a planing. This algorithms aims at minimizing a time-of-use cost of charging seen from the charging station owner/operator viewpoint. The battery model used for the optimal control algorithm is presented in Section III-A. Given the battery model, Section III-B describes the optimal control problem corresponding to the vehicles charging planing. That is to say, the cost to minimize, the equations of the dynamical systems and the optimization constraints of the problem.

A. Battery Model

In order to represent the battery dynamics, an equivalent circuit model comprising a voltage source u ocv , an electric resistance r c is used [START_REF] Guzzella | Vehicle Propulsion Systems[END_REF] and a maximal charging power p + c . These variables are dependant on the battery state of energy (SOE) x. In addition, these functions are approximated by polynomial functions of the SOE using experimental data. Using these functions and noting p c,k the charging power of vehicule k and x p c,k its corresponding SOE, one has the following dynamics

ẋk = u ocv (x k ) 2r c (x k ) × -u ocv (x k ) + u 2 ocv (x k ) + 4r c (x k )p c,k =f c (x k , p c,k ) (4)
Taking into account the SOE dependant charging power limitation one has

p c,k ≤ p + c,k (x k ) (5)

B. Optimal charging problem

In the problem we are interested in, chargers can be of different maximum charging power. Therefore, these chargers are gathered in n cs groups G j , j = 1, . . . , n cs , of same charging maximal charging power P + Gj . Therefore, the total amount of chargers is equal to the sum of elements of each groups G j :

j card(G j ) = n c (6)
The problem we are interested in consists in charging n x vehicles using n c chargers able to switch position to charge another vehicle every ∆T . Therefore, for all t ∈ [s∆T, (s + 1)∆T ), the connections between vehicles and chargers are represented by a matrix m s ∈ {0, 1} nx×ncs . Each element m s i,j of this matrix is equal to 1 if and only if vehicle i is connected to group charger j. Since each vehicle can be connected to at most one charger, one has

ncs j=1 m s i,j ≤ 1 ; ∀i = 1, . . . , n x , ∀s ∈ 0, . . . , k f (7)
In addition, the number of vehicles connected to the chargers of group j must be less or equal to the number of chargers within this group, which writes nx i=1 m s i,j ≤ card(G j ) ; ∀j = 1, . . . , n cs , ∀s ∈ 0, . . . , k f (8) Moreover, the total number of connected vehicles cannot exceed the minimum between the amount of vehicles or the amount of chargers:

nx i=1 ncs j=1 m s i,j ≤ min{n x , n c } ; ∀s = 0, . . . , k f (9)
Finally, the quantity to minimize is a trade-off between the co-charging operator cost and the charging speed. We assume a time-of-use charging rate to compute the charging cost and the charging speed is represented using a quadratic cost on difference between the current and the target state of charge. This trade-off cost writes

min pc,i,M k f k f s=0 (s+1)∆T s∆T price(t) nx i=1 ncs j=1 m s i,j p c,i (t) . . . + θ x (x i (t) -x i (t f,i )) 2 dt (10)
with k f = ⌈T /∆T ⌉. The weight θ x is a tuning parameter. For zero or small values of θ x emphasis is given to the charging cost while high values of θ x favor a faster charging. For the sake of simplicity, the following change of variable is performed

p c,i (t) = λ i (t)p + c,i (x i (t)) (11)
where

λ i (t) ∈ [0, 1]
. This change of variable allows the control variables to be all valued in the same set [0, 1]. Moreover, when connected to a station from group G j , the charging power of vehicle i is bounded as follows:

λ i (t)p + c,i (x i (t)) ≤ min{p + c,i (x i (t)), P + Gj } (12)
Finally, the NLOCPMIMT writes

min λi,M k f k f s=0 (s+1)∆T s∆T price(t) nx i=1 ncs j=1 m s i,j × . . . λ i (t)p + c,i (x i (t)) . . . + θ x (x i (t) -x i (t f,i )) 2 dt ( 13
)
under the following constraints

ẋi (t) = f c,i (x i (t), λ i (t)p + c,i (x i (t))) (14)
x

i (t) ∈ [0, E + i ] (15) λ i (t) ∈ [0, 1] (16) M k f = ⟨m 0 , . . . , m k f ⟩ (17) m s ∈ {0, 1} nx×ncs (18) ncs j=1 m s i,j ≤ 1 ( 19
) nx i=1 m s i,j ≤ card(G j ) ( 20
) nx i=1 ncs j=1 m s i,j ≤ min{n x , n c } (21)
P + Gj ≥ λ i (t)p + c,i (x i (t))m s i,j (22)
P s ≥ nx i=1 λ i (t) ncs j=1 m s i,j p + c,i (x i (t)) ( 23
)
where P s is the maximum power capacity of the charging station.

IV. A BRANCH AND BOUND APPROACH TO INTEGER OPTIMAL CONTROL PROBLEMS

A. Branch and Bound method overview

The problem presented in Section III-B is not a pure mixed integer optimal control problem such as in [START_REF] Sager | The integer approximation error in mixed-integer optimal control[END_REF]. Indeed, M k f is a discrete time variable which can switch from a discrete value to another only at regular time steps ∆T which can be large compared to the battery dynamics. Now, the method presented in [START_REF] Sager | The integer approximation error in mixed-integer optimal control[END_REF] converges to the optimal solution only if the time between consecutive switching times converges to zero. This convergence condition is not met for the presented application.

The mixed integer optimal control problem consists in finding both the optimal sequence of discrete valued control variable M k f and the continuous optimal controls λ i . To do so, we use a Branch and Bound approach [START_REF] Papadimitriou | Combinatorial Optimization, algorithms and complexity[END_REF]. Generally speaking, these techniques rely on a systematic exploration of a search-tree representing all possible values of the discrete variables. Throughout this exploration, at each node, the Branch and Bound algorithm completes two sequential tasks:

• Bounding : Computing the optimal solution of a relaxed problem for a node of depth k is an optimization problem where the discrete variables have been set over k stages, i.e. from the root to the current node of the search-tree. This relaxed optimization problem provides the algorithm a lower-bound of the actual optimal cost. • Branching : Once the bounding part of the node has been computed, the branching part can be performed. The branching of a node consists in two mutually exclusive actions :

-Exploring : If the branch rooted at the current node of depth k is worth exploring, the node's children are computed.

-Pruning : If the branch rooted at the current node is not worth exploring, i.e. if the optimal solution of the relaxed problem of the current node is higher than the cost of an admissible trajectory for the considered problem, the branch is then pruned and will not be explored. The efficiency of a Branch and Bound approach requires to have a relaxed problem which approximates the integer value problem and an exploring method, which explores primarily the nodes located on the optimal path of the tree. Now, let us see how to apply this general method to Problem ( 13)- [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF].

B. Lower bounding

The lower bounding part of the algorithm consists in solving a relaxed problem, i.e. a continuous problem deriving from the original problem where the integer constraints have been relaxed. In order to introduce the relaxed problem let us define the following useful functions

g 1 (λ i (t), x i (t), j, M k ) = . . . λ i (t)p + c,i (x i ) -max j P + Gj if t < k∆T λ i (t)p + c,i (x i )m s i,j -P + Gj , if t ≥ k∆T (24) g 2 (λ(t), x(t), M k ) = . . . i λ i (t)p + c,i (x i ) -P s , if t < k∆T i λ i (t)p + c,i (x i ) j m s i,j -P s , if t ≥ k∆T (25)
Function g 1 (resp. g2) corresponds to the relaxation of constraint (22) (resp. ( 23)) over the time interval [0, k∆T )

min λi J(λ, M k ) = . . . nx i=1 k s=0 (s+1)∆T s∆T price(t)λ i (t)p + c,i (x i (t))dt . . . + k f s=k+1 (s+1)∆T s∆T price(t) ncs j=1 m s i,j λ i (t)p + c,i (x i (t))dt . . . + t f 0 θ x (x i (t) -x i (t f,i )) 2 dt (26)
under constraints ( 14)-( 16) and where constraints ( 22)-( 23) are relaxed as follows:

g 1 (λ i (t), x i (t), j, M k ) ≤ 0 ∀t, i, j (27)
g 2 (λ(t), x(t), M k ) ≤ 0 ∀t, i (28)
M k = ⟨m 0 , . . . , m k ⟩ (29) m s ∈ {0, 1} nx×ncs (30) ncs j=1 m s i,j ≤ 1 (31) nx i=1 m s i,j ≤ card(G j ) (32) nx i=1 ncs j=1 m s i,j ≤ min{n x , n c } (33)
At this point, several remarks can be made regarding this optimal control problem

• Any node of depth k from the search-tree is uniquely associated to a particular sequence M k . • Notice from equations ( 24)-( 25) that the problem at depth k respects the integer constraints from stage k f -k to stage k f and is a continuous problem from stage 0 to stage k f -k. At depth k f the solution of the relaxed problem satisfies the integer constraints at all stages. Therefore, in this algorithm the setting of the integer matrices m k is a time-backward procedure. Although other methods of setting the integer variables are possible (forward time, by random,...), this method proved to be very efficient in practice. Each node generation consists in solving a problem of the form ( 26)-(33). Therefore, the computational tractability of the Branch and Bound method strongly relies on the ability to efficiently solve this problem and the solving algorithm must be designed with great attention. The chosen algorithm used to solve this problem is presented in details in Section V.

C. Upper bounding

A leaf node of the tree, i.e. a node of depth k f is a node for which all discrete variables have been set. Therefore this node satisfies the discrete constraints ( 17)-( 23), i.e. any optimal solution of a leaf node is admissible for problem ( 13)- [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF]. Thus, the optimal cost corresponding to the relaxed problem of a leaf node is an upper bound on the optimal cost of the original problem and at least one the leaf nodes is the optimal solution for Problem ( 13)- [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF].

D. Generating children of depth k + 1 from node of depth k

At this point, the relaxed problem of depth k has been described. One still has to define the branching part of the algorithm, i.e. how to generate the children of a node of depth k. Since Branch and Bound is an exhaustive exploration method, generating the children of depth k+1, of a parent node n(M k ) actually consists in generating all possible matrices m k+1 ∈ {0, 1} nx×ncs satisfying ( 19)-( 21) and solving the problem for each matrix. At first glance generating these matrices is not straightforward. But, there is a one to one correspondence between this set and the set of all possible vehicle to charger connection configurations and one can directly derive the connection matrix m s corresponding to any of these configurations.

E. Exploring the tree

The branching part consisting in lower-bounding the optimal cost and generating the children of a node has been covered. To fully describe the Branch and Bound algorithm, one has to specify the exploration part of this algorithm, that is to say in what order the nodes are going to be explored.

1) Depth-First-Search like strategy: In order to have an efficient pruning of the tree, it is very important to get a good upper bound on the optimal cost as fast as possible. Indeed, as long as no upper bound has been computed no pruning can be performed. In the same time, the convergence speed of the Branch and Bound method strongly relies on the ability to prune branches in early stages of the algorithm execution. But, as described in Section IV-C, upper bounds on the optimal cost are computed only in the leaf nodes of the tree. Therefore, it is very important to reach a leaf node as fast as possible. To do so, a Depth-First-Search like technique of exploration is used. Depth-First-Search techniques consist in putting generated children of a node in a stack and then exploring the node on top of the stack. Proceeding that way ensures that the current explored node is always one step deeper than the previous one, and that a leaf node is reached in k f steps. But, in order to get a good upper bound on the optimal cost, the generated nodes are stacked in decreasing order of their relaxed optimal cost. Doing this way ensures that the Depth-First-Search method follows a path of best children per generation. Indeed, a node is explored, then its best child is explored, then its best child's best child and so on. This heuristic allows the algorithm to find a good upper bound in k f steps. Moreover, using a backward strategy for setting the integer variable (see equation ( 24)-( 25)) often allows this strategy to find the optimal solution when reaching a leaf node.

2) Best first search strategy: Once the Depth-First-Search strategy has reached a leaf node, an upper bound is available.

To complete the optimization it is now time to prune as many branch as possible to ensure the optimal solution is reached as fast as possible. To do so, the exploration method is now a Best-First-Search method, the node with the lower relaxed optimal cost is explored first. The data structure behind this Best-First-Search is a min-priority-queue [START_REF] Cormen | Introduction to Algorithms[END_REF] sorting nodes in increasing order of their relaxed optimal cost.

F. Branch and Bound algorithm

The complete Branch and Bound algorithm is presented in Algorithm 1.

Remark 1. At the end of the Depth-First-Search part of the Branch and Bound algorithm a good solution (often optimal in practice) is reached and can be used as an approximate solution of the problem. This property is particularly useful if the execution time must be controlled.

V. SOLVING THE LOWER BOUNDING PROBLEM : AN INTERIOR POINT METHOD APPROACH

As stated in Section IV-B, the solving algorithm of the lower-bound problem must be chosen carefully. First, the algorithm must be fast enough to ensure that the node generation will not slow down the tree exploration. Interestingly, at each node generation, one can notice that the solution of the lowerbounding procedure of the parent node is a good initialization for any of its child nodes. Indeed, the only difference between the optimization problems of node k and its children of depth k + 1 are functions g 1 and g 2 from equations ( 24)-( 25) on the time interval [(k -1)∆T, k∆T ). In order to compute the lower bound of each node, Pontryagin Maximum Principle based methods are suitable [START_REF] Bryson | Applied optimal control[END_REF] because of their convergence speed when provided with a good initialization. However, in the framework of optimal control, handling pure state and control constraints such as in equations ( 15) and ( 16) is Algorithm 1: Branch and Bound algorithm end for 40: end while 41: return best node not straightforward [START_REF] Hartl | A survey of the maximum principles for optimal control problems with state constraints[END_REF] and classic off the shelves TPBVP solvers fail to compute the solution of such problems 1 . To overcome this difficulty, one can use Interior Point Methods (IPMs) adapted for state and input constraints optimal control problems.

1: stack ← Stack() 2: ub ← +∞ 3: lb ← -∞ 4: E -1 ← ∅ 5: root ← node(E -

A. Presentation of the interior penalized optimal control problem

In order to solve the lower bounding problem ( 26)-(33) of the Branch and Bound algorithm, an interior point approach 1 The adjoint variables of the two-point boundary value problems can be discontinuous when the state constraint reaches 0.

[21] is used. This method consists in relaxing constraints (15)-( 16) and ( 27)-( 28) through an augmentation of the cost (26) using interior penalization functions. This relaxation depends on a parameter sequence (ε n ) converging to 0. As this sequence converges to 0 the solution of the relaxed optimal solution converges to the solution of the original problem. The penalized problem writes as follows

min λ∈Λ M (λ, M k , ε) = . . . t f t0 price i p + c,i (x i )λ i + θ x (x i (t) -x i (t f,i )) 2 dt + εp int (x, λ, M k ) (34)
under constraints [START_REF] Buss | Nonlinear hybrid dynamical systems: Modeling, optimal control, and applications[END_REF] and where

p int (x, λ, M k ) = t f 0 i γ x (-x i ) + γ x (x i -E + i ) + γ x • g 2 (λ, x, M k ) + i γ λ (λ i ) + i j γ x • g 1 (λ i , x i , j, M k ) dt (35) Λ =L ∞ ([0, T ]; [0, 1] nx ) (36) γ x (x) = (-x) -1.1 if x < 0 0 otherwise (37) γ λ (x) = -log(-x) -log(1 -x) (38)
This particular choice of penalty is given in [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of non linear systems[END_REF] and ensures that the optimal solution of the penalized problem strictly satisfies the constraints.

Theorem 1 (From [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of non linear systems[END_REF]). If penalty functions γ x and γ λ are chosen according to equations (37)-(38) then any optimal solution λ * of problem (34) with ε > 0 is such that

x * i (t) ∈ (0, E + i ) (39) λ * i (t) ∈ (0, 1) (40) g 1 (λ * i (t), x * i (t), j, M k ) < 0 (41) g 2 (λ * (t), x * (t), M k ) < 0 (42)
B. Solving algorithm 1) Saturation functions for a fully unconstrained problem: In order to get a fully unconstrained problem, we use saturation functions on the control:

λ = ϕ(ν) = 1 2 (1 + tanh(ν)) (43)
Using this change of variables the penalized optimal control problem (34) becomes

min ν M (ϕ(ν), M k , ε) = . . . t f 0 price i p + c,i (x i )ϕ(ν i ) + θ x (x i (t) -x i (t f,i )) 2 dt . . . + ε p int (x, ϕ(ν), M k ) (44)
under the following dynamic constraint adapted from ( 14):

ẋi = f c,i (x i , ϕ(ν i )) (45)
Thanks to Theorem 1, this change of variables is well posed since the optimal control λ * for problem (34) + ( 14) is in L ∞ ([0, t f ]; (0, 1)), i.e. is strictly interior.

2) Pontryagin Maximum Principle approach for solving the unconstrained optimal control problem: Solving the lower bounding problem consists in solving a sequence of unconstrained optimal control problems (UOCPs) (44)-(45). To solve this sequence of UOCPs, a Pontryagin Maximum Principle (PMP) based approach is used. To do so, let us first define the Hamiltonian [START_REF] Bryson | Applied optimal control[END_REF] of problem (44)-( 45)

H(t, x, ϕ(ν), p, M k , ε) = . . . price i p + c,i (x i )ϕ(ν i ) + θ x (x i (t) -x i (t f,i )) 2 + ε p int (x, ϕ(ν), M k ) + i p T i f c,i (x i , ϕ(ν i )) (46)
The PMP states that any optimal solution from problem (44)-( 45) is solution of the following TPBVP

ẋi = f c,i (x i , ϕ(ν i )) (47) ṗi = - ∂H(t, x, ϕ(ν), p, M k , ε) ∂x i (48) 0 = ∂H(t, x, ϕ(ν), p, M k , ε) ∂ν i (49) x i (0) = x 0,i (50)
x i (t f ) = x t f ,i (51)
To solve this TPBVP numerous algorithms such as shooting or collocation algorithms exist and an exhaustive presentation is given in [START_REF] Ascher | Numerical solution of boundary value problems for ordinary differential equations[END_REF]. For the problem under consideration, a collocation based solving algorithm such as [START_REF] Shampine | Solving boundary value problems for ordinary differential equations in matlab with bvp4c[END_REF] is a relevant choice. Indeed, these methods allow one to achieve a good compromise between computational speed and numerical sensitivity.

3) Presentation of the lower bounding algorithm: In the following we present the resolution algorithm for the lower bounding part of the Branch and Bound Algorithm 1 Algorithm 2: Lower bounding algorithm end if 11: end while 12: return

1: ε = ε 0 > 0, α > 1, tol > 0 2: Initialization of (ν, x) such that x i ∈ (0, E + i ); ϕ(ν i )p + c,i (x i ) ≤ max j P b,j and i ϕ(ν i )p + c,i (x i ) ≤ P s 3: C int ← ⊥ 4: while ¬C int do 5: (ν, x, p) is solution of (47)-(51) 6: if εp int (x, ϕ(ν), M k )/J • ϕ(ν) < tol
t f 0 price i p + ci (x i )ϕ(ν i )dt
The exit condition from Algorithm 2 allows to stop the algorithm when the perturbation on the optimal cost provided by the penalty function εp int (x, ϕ(ν), M k ) is negligible with respect to the original cost. The working justification of this exit condition will be detailed in Section VI-A.

VI. CONVERGENCE ANALYSIS

Before proving convergence of the presented method, let us introduce the following Definition Definition 1 (Lower bound solution). The lower bound solution of a node n of depth k is a 2-uple (M k n , λ n ), where M k n is the sequence of integer matrices corresponding to node n defined in equations (29)-( 33) and where λ n is an optimal solution of problem ( 26)-( 33) with M k = M k n . Now, the proof of convergence is presented in three parts. First, the convergence of the interior point method from Section V is investigated. Then, proof of path-increasing cost is given. Finally the proof of convergence of the method directly stems from the aforementioned proofs.

A. Convergence of the interior point method

In order to prove the convergence of the interior point method, one first needs to prove that the solution of the Two Point Boundary Value Problem (TPBVP) (47)-( 51) is indeed the optimal solution of the considered problem. This is the object of the following proposition Proposition 1. There exists a unique optimal solution λ * to Problem (34)-( 38) and λ * is solution of TPBVP (47)-(51)

Proof: See Appendix A The proof of convergence of the interior point method is a direct use of [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of non linear systems[END_REF].

Theorem 2 (From [START_REF] Malisani | An interior penalty method for optimal control problems with state and input constraints of non linear systems[END_REF]). If the penalty functions are chosen according to equations (37)-(38), any optimal solution

ν * ε (M k ) of problem (44)-(45) is such that, ∀k ≤ k f lim ε→0 M (ϕ(ν * ε (M k )), M k , ε) = min λ J(λ, M k ) (52) lim ε→0 J(ϕ(ν * ε (M k )), M k ) = min λ J(λ, M k ) (53) lim ε→0 εp int (x, ϕ(ν * ε (M k )), M k ) = 0 (54)
under constraints ( 14)-( 16)+( 27)-( 33).

This result yields lim

ε→0 εpint(x,ϕ(ν * ε (M k )),M k ) J(ϕ(ν * ε (M k )),M k )
= 0, thus the exit condition from Algorithm 2 B. Path increasing lower bound procedure Proposition 2. Let n 1 and n 2 be two nodes on the same simple path from the root to a leaf of the tree. Without loss of generality, assume that n 1 (resp. n 2 ) has depth k (resp. s, s > k) and let (M k n1 , λ n1 ) (resp. (M s n2 , λ n2 )) be its lower bounding solution from Definition 1 . One has J(λ n1 , M k n1 ) ≤ J(λ n2 , M s n2 ). In other words, the lower bounding cost is increasing along any simple path of the search-tree.

Proof: Since n 1 and n 2 are on the same path connecting the root and a leaf of the search tree, one has M s n2 = M k n1 : ⟨m k+1 , . . . , m s ⟩. Hence, using equation ( 24) one has

g 1 (λ n1 i , x n1 i , j, M k n1 ) ≥ . . . (55) g 1 (λ n2 i , x n2 i , j, M k n2 ) ∀t ∈ [0, t f ], . . . (56) i = 1, . . . , n x ; j = 1, . . . , n b (57)
where x n1 = (x n1 i , . . . , x n1 nx ) T is solution of equations ( 14), with control variable λ n1 = (λ n1 i , . . . , λ n1 nx ) T and using equation ( 25) one also has

g 2 (λ n1 , x n1 , M k n1 ) ≥ . . . ( 58
)
g 2 (λ n2 , x n2 , M k n2 ) ∀t ∈ [0, t f ] (59)
Therefore, the search space of problem ( 26)-( 33) associated to M k n1 contains the search space associated to M s n2 . This yields J(λ n1 , M k n1 ) ≤ J(λ n2 , M s n2 ).

C. Convergence of the Branch and Bound method

We are now ready to state our main convergence result

lim ε↓0 J(ϕ(ν * ϵ (M k f o )), M k f o ) = J(λ o , M k f o ) = min M k f min λ J(λ, M k f )(60)
under constraints ( 14)-( 23)

Proof: From Theorem 2 and Definition 1 one has

lim ε↓0 J(ϕ(ν * ϵ (M k f o )), M k f o ) = min λ J(λ, M k f o ) = J(λ o , M k f o ) (61
) Now, let us prove that node o is the leaf node with lowest lower bounding cost among all leaf nodes. The proof is by contradiction. Assume that there exists a leaf node l with lower bounding solution (M

k f l , λ l ) such that J(λ l , M k f l ) < J(λ o , M k f o ).
This configuration is possible only if during the Branch and Bound algorithm execution a branch containing l rooted at a node m of depth s < k f with lower bounding solution (M s m , λ m ) is pruned. Let c ub be the current upper bound when this pruning happened, according to Proposition 2 one has c ub ≤ J(λ m , M s m ) ≤ J(λ l , M k f l ) otherwise no pruning would have occurred. When the execution is completed, the final upper bound is the lower bounding cost of the leaf node o. From the updating method of the Branch and Bound algorithm upper bound one has

J(λ o , M k f o ) ≤ c ub ≤ J(λ m , M s m ) ≤ J(λ l , M k f l )
for any leaf node l, which contradicts the initial assumption and yields

J(λ o , M k f o ) = min n∈L J(λ n , M k f n ) ( 62
)
where L is the set of all leaf nodes. By construction of the Branch and Bound search-tree, the set {M k f n } n∈L is exactly the set of all discrete sequences of length k f which yields

J(λ o , M k f o ) = min n∈L J(λ n , M k f n ) = min M k f min λ J(λ, M k f )
(63) and concludes the proof.

VII. NUMERICAL PERFORMANCE ANALYSIS

A. Example presentation

In this section the following n x -dependent problem is considered n c = 2 ; P c = {50, 22} ; P s = 60 (64) ∆T = 30 minutes (65)

θ x = 0 (66) t 0,k = 0, ; t f,k ∈ U (14, 18) (67) E + k ∈ U (30, 90) (68)
x 0,k ∈ U E + k 10 , E + k 5 , ; x f,k ∈ U 7 10 E + k , 9 10 E + k (69) price = 0, 095 if t ∈ [6, 9) ∪ [12, 14) 0.07 otherwise (70)
where U (x, y) is the uniform distribution on [x, y]. From this setting with n x vehicles and 2 different chargers, the size of the Branch and Bound search-tree is exactly the number of permutations of 2-elements among n x to the power of the time steps number, that is to say nx! (nx-2)! 36 .

B. Reference charging strategy

In order to benchmark the proposed method, we consider the following reference strategy. The idea is to charge two vehicles which need the larger amount of energy at the switching moment. More precisely, at switching time t s , we note i 1 (t s ) (resp. i 2 (t s )) the vehicle such that x f,i1(ts) -x i1(ts) (t s ) is maximum (resp. second to maximum) and let p c,k be the charging power of vehicle k. The charging strategy is as follows

p c,i1(ts) (t) = min{50, p c,i1(ts) (x i1(ts) (t))} (71) p c,i2(ts) (t) = min{22, p c,i2(ts) (x i2(ts) (t)); 60 -p c,i1(ts) (t)} (72) p c,k (t) =0 ∀k ∈ {1, . . . , n x } \ {i 1 (t s ), i 2 (t s )} (73)
where t ∈ [t s , t s+1 ). Note that once the state of energy of vehicle k has reached x f,k , the charging power of this vehicle is always 0.

C. Performance analysis indicators

For each value of n x we solve 50 instances of problem (64)-(70) both with reference and presented method. Let e i opt (nx) (resp. e i ref (nx)) be the optimal solution of instance i = 1, . . . , 50, of problem (64)-(70). In the same way, we note c i opt (nx) (resp. c i ref (nx)) the optimal (resp. reference) charging cost corresponding to the solving of instance

(resp. δ i c (n x ) = c i ref (nx)-c i opt (nx) c i ref (nx)
) the relative energy (resp. cost) savings. In order to quantify the performances of the proposed method, the following indicators are used

• Relative Mean Energy Savings, RMES(n x ):

RMES(n x ) = 100E i δ i e (nx) (74)
• Energy Relative Standard Deviation, ERSD(n x ):

ERSD(n x ) = 100 E i {δ i e (nx) 2 } -E i {δ i e (nx)} 2 (75)
• Relative Mean Cost Savings, RMCS(n x ):

CRMS(n x ) = 100E i δ i c (nx) (76)
• Cost Relative Standard Deviation, CRSD(n x ):

CRSD(n x ) = 100 E i {δ i c (nx) 2 } -E i {δ i c (nx)} 2 (77)
where E i is the mathematical expectation with respect to the instance number i.

D. Compared performance analysis 1) Energy consumption comparison:

On Table I, one can see that the relative energy savings decreases with the number of vehicle to be charged. This phenomenon is a consequence of the fact that fewer vehicles to charge over a given time period allows to use lower power charging limiting the energy losses in the process. In contrast, the reference method aims at charging EVs as fast as possible to the expense of energy losses.

2) Cost savings comparison: On Table I, the behavior of the RMCS as a function of n x is more complex, it decreases from n x = 3 to n x = 7 and increases from n x = 7 to n x = 9. First, for n x = 3, . . . , 7, the cost savings is a decreasing function of n x . On Table II, one can see that the percentage of energy consumed during peak hours is zero 2 for reference and optimal strategies when n x < 7. Since the whole energy is consumed outside peak hours, the percentage of cost decrease is equal to the percentage of energy decrease. From Table I, the CRMS is now an increasing function of n x , for n x ≥ 7. This increase phenomenon is due to the fact that the main 2 In the optimal case, the peak-hours energy is not exactly zero. This stems from using IPMs to solve the optimal control which leaves a residual non zero consumption during these periods. source of cost savings does not come from limiting energy losses anymore but from shifting energy consumption from the peak periods. On Table II, one can see that the percentage of energy consumed during peak hours increases when n x ≥ 7 in the reference case but is still 0 with the optimal strategy. Therefore, the more energy is consumed during peak hours for the reference strategy the greater the cost savings using the proposed method.

3) Execution time analysis: All instances of Problem ( 64)-( 70) are solved on a Linux Ubuntu operating system endowed with 35 Intel®Xeon®Gold 5120 CPU @ 2.2 GHz. The solver is based on homemade PYTHON solvers both for the Branch and Bound algorithm and the Differential Algebraic Equations. On Table I one can see that the growth in execution time suddenly increases for n x ≥ 7. This phenomenon stems from the fact that for n x ≥ 7 there are more children by generation than CPUs available. Hence the children generation part of the Branch and Bound algorithm cannot be fully parallelized. That said, one can see that the mean execution time is always less than 1% of the time horizon. Finally, with this hardware setting the computation time is compatible with an online implementation. Nevertheless, the presented method requires lots of computing power. Specifically, this methods relies on the availability of using many CPUs at the same time which can be a limit for online implementation.

VIII. MPC IMPLEMENTATION SIMULATION

A. MPC implementation algorithm

In this Section, we are interested in comparing the performances of the reference method presented in Section VII-B and the presented algorithm in an online implementation simulation. For this simulation, we implement the algorithm in a straightforward Model Predictive Control manner [START_REF] Morari | Model predictive control: past, present and future[END_REF] consisting in the following step:

• At each vehicle arrival solve the optimal charging problem corresponding to the charging of the cars parked in the co-charging station. • Execute the computed planning until new vehicle arrival. This implementation does not require any prediction on the arrival time of a vehicle. Each vehicle is taken into account in the optimization problem as soon as it arrives at the charging station. Therefore, if the charging speed is too low handling a large flow of vehicle arrival might be impossible. In order to gain some robustness with respect to the flow of vehicles the weight on the quadratic cost on the state of charge θ x can be tuned to achieve a good trade-off between minimization of the time-of-use charging cost and the charging speed.

B. Simulation Parameters

The MPC simulation is carried out with the following parameters n c = 2 ; P c = {20, 10} ; P s = 25

(78) ∆T = 30 minutes, ; θ x = 10 -3 (79)

E + k ∈ U (5, 15) (80)
x

0,k ∈ U E + k 10 , E + k 5 , ; x f,k ∈ U 7E + k 10 , 9E + k 10 (81) price =    0, 095 if t ∈ [6, 9) 0, 150 if t ∈ [18, 21) 0.07 otherwise (82)
In addition the probability of a car arrival of a car between two plug-times (30 min) of the charger is 40%. The simulation is carried out on 48 hours.

C. Simulation Results

This comparative simulation required the solving of 41 problems corresponding to 41 vehicle arrival. At max, there were 6 cars charging at the same time. Table IV shows that the optimal strategy is both energy and time-of-use cost efficient. The overall energy savings stems from a slower but more efficient charging of the vehicles. The lower charging rate is particularly visible by comparing Figures 1 and2 where it is clear that the required final state-of-energy is reached sooner with reference than optimal strategy. Moreover, Table III shows that the optimal strategy achieves to reduce by 75% the energy consumption during peak periods. This energy consumption is simply transferred to off-peak periods. This transfer is the most important factor of charging-cost reduction. In addition, one can see on Figures 3 and4 that the optimal charging strategy takes advantage of the possibility of regularly changing the vehicles connected to the charging stations. This increase rate of plug-in plug-out of vehicles allows to reduce the charging speed while satisfying the required state of energy target.

IX. CONCLUSION

This paper presents a methodology to compute an optimal charging scheduling of electric vehicles at a co-charging station. The presented methodology relies on a Branch and Bound strategy allowing to compute the global optimum. This methodology strongly relies on Maximum Principle based optimal control solver allowing to compute quickly a tight lower bound at each node of the Branch and Bound tree. It also relies on the slower time scale of the charging stations compared to the time scale of vehicles batteries to produce an efficient Branch and Bound algorithm. However, in order to be working at full potential this algorithm requires to run on a cluster so that the children generation can be fully parallelized. This need of computing power might hold back online implementations of this method. Nevertheless, thanks to the global convergence of the method, it can be used offline to assess the performances of heuristic-based methods in terms of distance the optimal solution. In future work, the authors will to take into account uncertainties on the departure time and arrival time of the electric vehicles. A progressive hedging approach [START_REF] Rockafellar | Scenarios and policy aggregation in optimization under uncertainty[END_REF] might be considered to use the presented algorithm in a stochastic framework.

APPENDIX

A. Proof of Proposition 1

First, let us recall the Inverse Function Theorem which will be used throughout the proof.

Theorem 4 (Inverse Function Theorem [26]). Let f : R n → R n be continuously differentiable on some open set containing a, and suppose det(∇f (a)) ̸ = 0. Then there is some open set V containing a and an open set W containing f (a) such that f : V → W has a continuous inverse f -1 : W → V which is differentiable for all y ∈ W and one has

∇ f -1 (y) = [∇f ] -1 (f -1 (y))
In addition, the Ky Fan Dominance Theorem recalled here after will also be used Theorem 5 (K. Fan Dominance Theorem [START_REF] Fan | Some metric inequalities in the space of matrices[END_REF]). Let M n be the space of n × n real matrices. For all A ∈ M n , we note (s i (A)) i=1,...,n the sequence of singular values in decreasing order. For 1 ≤ k ≤ n, the Ky Fan k-norm on M n is defined as

∥ A ∥ (k) = k i=1 s i (A)
These norms are unitarily invariant, i.e., ∥ U AV ∥ (k) =∥ A ∥ (k) and for given matrices A, B ∈ M n ∥ A ∥ (k) ≤∥ B ∥ (k) ⇔∥ A ∥≤∥ B ∥ ∀k = 1, . . . , n and for all unitarily invariant norms ∥ . ∥, such as the Euclidian ∥ . ∥ 2 norm.

The proof of Proposition 1 follows the following outline

• First, we prove that any optimal trajectory from Problem (34)-(38) does not contain any singular arc in Section A1. • Then, we prove that any optimal solution of Problem (34)-( 38) is characterized by the stationarity conditions (47)-(51) in Section A2. • Finally, we prove that there exists a unique solution from Problem (47)-(51) in Section A3.

1) Absence of singular arcs: Since the control variable λ appears explicitely in ∂ ∂λ H(t, x(t), λ, p(t), M k , ε) for all values of t, x(t), p(t), M k , ε, there are not singular arcs along the optimal trajectory, i.e. it is always possible to compute a minimizer of H(t, x(t), λ, p(t), M k , ε) with respect to λ ∈ [0, 1] nx .

2) Stationarity conditions of Pontryagin's Maximum Principle: From Theorem 1, any optimal solution of Problem (34)-( 38) is strictly interior to the constraints. Therefore, any optimal solution λ * of this problem satisfies the following stationarity condition: ∂H(t, x(t), λ * (t), p(t), M k , ε) ∂λ = 0

3) Uniqueness of solutions of the stationarity conditions: In the following, for any mapping φ : X × Y → Z, we note φ[x](.), the mapping φ(x, .) : Y → Z where x is fixed and considered as a parameter. Using that notation we introduce the following functions From Theorem 4 one has ∥ r[t, x(t), p(t), M k , ε](λ 1 ) -r[t, x(t), p(t),

F [t
M k , ε](λ 2 ) ∥ 2 ≤ 1 8 ∥ λ 1 -λ 2 ∥ 2
Mapping r[t, x(t), p(t), M k , ε](.) is 1 8 -Lipschitz, hence it is a contraction. Using Banach fixed point Theorem one has that for all t, x(t), p(t), ε, M k , there exists a unique solution λ * to TPBVP (47)-(51).

  then

Theorem 3 .

 3 Let node o be the solution of the Branch and Bound algorithm, let (λ o , M k f o ) its lower bounding solution and let ν * ϵ (M k f o ) be the optimal solution of the lower bounding problem (47)-(51) for node o. The Branch and Bound algorithm converges to the global solution of problem (13)-(23) in the following sense

  nx

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .

 1234 Fig. 1. Optimal energy time series for MPC simulation

1 (

 1 2 E[x(t), M k ](λ)) ≥ s n (∇ 2 D[x(t), M k ](λ)), where (s n (.)) n is the sequence of eigenvalues in decreasing order. For any Ky-Fan norm[START_REF] Bhatia | Matrix Analysis[END_REF] ∥ . ∥ (k) = k i=1 s i (.), one has ∥ ∇ 2 D[x(t), M k ] + ∇ 2 E[x(t), M k ] λ) ∥ (k) = . . .

nx i=nx-k+1 1 s 1 s

 11 i ∇ 2 D[x(t), M k ](λ) + ∇ 2 E[x(t), M k ](λ) i ∇ 2 D[x(t), M k ](λ) =∥ ∇ 2 D[x(t), M k ]

  cs • P + Gj : Maximum power of group charger j = 1, . . . , n cs • m ∈ {0, 1} nx×ncs : Vehicle-charger connection matrix. Vehicle i is connected to group chargerj if and only if the element m i,j of m is 1. • P s ∈ R: maximal power available for the whole charging station • t 0,k : earliest plugging time of the vehicle k • t f,k : latest unplugging time of the vehicle k

• t f : optimization horizon • ∆T : minimum duration of a vehicle connection to a charger

ref (n x ) = nx k=1 p c,k dt: Sequence of total energy consumption with reference charging powers p c . • c ref (n x ) = price nx k=1 p c,k dt: total charging cost with reference charging powers p c .

  weight of state of charge penalization • e opt (n x ) =

	nx k=1 λ k p + c,k (x λ k )dt: total energy con-
	sumption with optimal charging rates λ.
	• c opt (n x ) = price	nx k=1 λ k p + c,k (x λ k )dt: total charging
	cost with optimal charging rates λ.
	• ⊤, ⊥: True and False boolean variables respectively
	• ¬: logical not
	• RMES : Relative Mean Energy Savings
	• ERSD : Energy Relative Standard Deviation
	• RMCS : Relative Mean Cost Savings
	• CRSD : Cost Relative Standard Deviation
	In the rest of the paper, bold lowercase variables such as m
	indicate matrices of R nx×ncs . Bold capital variables such as
	M k indicates a sequence of such matrices also noted

• e

  , x(t), p(t)](λ) = price(t) + ∂f c (x(t), λp + M k ](λ) is the gradient of the sum of convex functions over a convex set. Hence, the corresponding Hessian is strictly definite positive and invertible. Moreover, one haslim λ→1 -Γ i [x(t), M k ](λ) = +∞, i = 1, . . . , n x Therefore, Γ[x(t), M k ] is a continuous mapping from (0, 1) nx to R nx . Since x, p are bounded and since F [t, x(t), p(t)] is continuous F [t,x(t),p(t)]ε is also bounded for all ε > 0. Therefore, one can use the Inverse Function Theorem which givesλ = Γ[x(t), M k ] [t, x(t), p(t), M k , ε](λ) ∥= . . . D[x(t), M k ] + ∇ 2 E[x(t), M k ] Since ∇ 2 D[x(t), M k ](λ) and ∇ 2 E[x(t), M k ](λ)are positive symmetric definite matrices, their singular are the eigenvalues and s n (∇ 2 D[x(t), M k ](λ) + ∇

	sup λ	∥ ∇ Γ[x(t), p(t), M k ] -1 -	F [t, x(t), p(t)](λ) ε	∥
	From the inverse function theorem one has
	sup λ	∥ ∇ Γ[x(t), p(t), M k ] -1 -	F [t, x(t), p(t)](λ) ε	∥= . . .
	sup	∥ ∇Γ[x(t), p(t), M k ]	-1	. . .
	λ
	Γ[x(t), p(t), M k ]	-1	-	F [t, x(t), p(t)](λ) ε	∥
	= sup	∥ ∇ 2 -1	. . .
		λ
	Γ[x(t), p(t), M k ]	-1	-	F [t, x(t), p(t)](λ) ε	∥
						∂λ	c (x(t))	p(t)
						D[x(t), M k ](λ) =	∂ ∂λ
						lim λ→0 +	Γ -1	-	F [t, x(t), p(t)](λ) ε
						= r[t, x(t), p(t), M k , ε](λ)	(83)

i γ λ (λ i ) E[x(t), M k ](λ) = ∂g 2 ∂λ (λ, x(t), M k )γ ′ x • g 2 (λ, x(t), M k ) + i,j ∂g 1 ∂λ (λ i , x i (t), j, M k )γ ′ x • g 1 (λ i , x i (t), j, M k ) Γ[x(t), M k ](λ) = D[x(t), M k ](λ) + E[x(t), M k ](λ)

and one has:

∂H(x(t), λ, p(t), M k , ε, t) ∂λ = . . . F [t, x(t), p(t)](λ) + εΓ[x(t), M k ](λ) Γ[x(t), i [x(t), M k ](λ) = -∞, i = 1, . . . , n x

In order to prove existence and uniqueness of solutions of (83), one needs to prove that the mapping r[t, x(t), p(t),

M k , ε](.) is a contraction. ∥ r[t, x(t), p(t), M k , ε](λ 1 ) -r[t, x(t), p(t), M k , ε](λ 2 ) ∥ . . . ≤ sup λ ∥ ∇r[t, x(t), M k , p(t), ε](λ) ∥∥ λ 1 -λ 2 ∥

Plus, one has sup λ ∥ ∇r