Andres Ovalle
email: andres.ovalle-villamil@g2elab.grenoble-inp.fr

Seddik Bacha

Ahmad Hably

Kaustav Basu

On the most convenient Mixed Strategies in a Mixed Strategist Dynamics Approach for Load Management of Electric Vehicle Fleets

Keywords: Plug-In Hybrid Electric Vehicles, Smart charging, Evolutionary Game theory, Distributed Optimization, Mixed Strategist Dynamics

This manuscript explores the selection of appropriate mixed strategies (MSs) in a Mixed Strategist Dynamics (MSD) application for load management of Plug-in Electric Vehicle (PEV) fleets. This selection is based on the convenience of PEV owners, aiming to choose those MSs that privilege early high (or fast) charging rates when it is possible. The previously published MSD and Maximum Entropy principle (MSD-MEP) approach is revised and illustrated with several examples, specially in the context of selection of MSs. This revision allows a wider understanding of the method, and aims to inspire new contributions on domains where distributed optimization methods are pertinent. Results obtained without any management structure are compared to those obtained with the MSD-MEP approach under different scenarios, where full sets of MSs and reduced sets of convenient MSs are applied. The performance of the method, using conveniently reduced sets of MSs, is tested with real historical active power measurements, provided by the SOREA utility grid company in the region of Savoie, France.

I. INTRODUCTION

Maturing and emerging technologies have boosted the markets of plug-in electric vehicles (PEVs) and plug-in hybrid electric vehicles in the recent years, in terms of battery capacities and autonomy. These technological advancements have impacts not only in the development of sustainable mobility systems [START_REF] Helmbrecht | How electric vehicles affect driving behavioral patterns[END_REF], but also on the electrical infrastructure that provides power to these systems [START_REF] Jian | Regulated charging of plug-in hybrid electric vehicles for minimizing load variance in household smart microgrid[END_REF]- [START_REF] Nguyen | Optimal charging and discharging for multiple phevs with demand side management in vehicle-to-building[END_REF]. Connected PEVs represent distributed energy storage devices to the electricity distribution system, that can be employed to improve power quality and efficiency or reduce maintenance costs [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF]. Nevertheless, without any smart management strategy handling charging infrastructures, distribution systems are prone to have problems such as heavy load unbalances, high load peaks and fast transitions, or transformers lifespan reduction [START_REF] Dyke | The impact of transport electrification on electrical networks[END_REF]- [START_REF] Clement-Nyns | The impact of charging plug-in hybrid electric vehicles on a residential distribution grid[END_REF].

Several approaches have been proposed in order to properly handle PEV fleets load. In particular, decentralized approaches are preferable because of the benefits the provide for owners autonomy [START_REF] He | Optimal scheduling for charging and discharging of electric vehicles[END_REF]. In [START_REF] Rezaei | Packetized plug-in electric vehicle charge management[END_REF], a communication channel inspired approach is proposed, where each PEV load is divided in several packets. Each PEV asks for permissions in order to consume its programmed packets of load, according to the congestion of the system. Authors of [START_REF] Richardson | Local versus centralized charging strategies for electric vehicles in low voltage distribution systems[END_REF] propose an approach for low penetration of PEVs, where each one maximizes its charging rates taking into account voltage and current constraints based on offline parameter computations. Uncertain final states of charge, sudden variations in power consumption heavily affecting batteries and transformers, and inflexible charging rates, are some of the drawbacks of these approaches.

Other authors formulate the load scheduling problem as Nperson potential games (a concept from game theory), where each PEV is a player [START_REF] Nguyen | Optimal charging and discharging for multiple phevs with demand side management in vehicle-to-building[END_REF], [START_REF] Ma | Decentralized charging control of large populations of plug-in electric vehicles[END_REF]. These techniques optimize a certain criteria by defining proper load schedules in a decentralized way. Some weak points in these approaches arise when the number of PEVs is not large enough, or when parameters (total PEV load, charger power limits, connection times) are not homogeneous. Moreover, these methods unfairly allocate the available resources. In [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], an evolutionary game dynamics approach is proposed based on the application of the Mixed strategist Dynamics (MSD) and the Maximum Entropy principle (MEP) in a multi-population model. This method handles the distributed PEVs load allocation problem in a fair way, reducing the impact on the grid, and taking into account social constraints from owners side. It takes advantage of the concept of Mixed Strategies (Ms) and applies it in order to define local sets of feasible solutions as convex hulls, where the local dynamics evolve.

In this paper, the previously published MSD-MEP approach of [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], is explored in terms of the selection of appropriate MSs. This selection is based on the convenience of PEV owners, aiming to choose those MSs that privilege early high (or fast) charging rates when it is viable. Several details and examples of the MSD-MEP method and the proposed appropriate selection of MSs are provided throughout the manuscript.

The paper is organized as follows. Section II provides a brief introduction to the MSD, with some examples clarifying the definition and usage of MSs. Section III summarizes the MSD-MEP application, and introduces illustrative examples. Section IV describes the role of MSs in the MSD approach. Besides, it describes the selection of MSs privileging early high charging rates, among a full set of MSs. Section V provides results under several illustrative and realistic scenarios.

II. SUMMARY AND SOME DETAILS ON THE MSD

As in [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], [START_REF] Hofbauer | The Theory of Evolution and Dynamical systems[END_REF], [START_REF]Evolutionary Games and Population Dynamics[END_REF], the MSD describes the evolution of the distribution of a population whose individuals are able to choose among an amount M of MSs in a predefined set, depending on the payoff they get from each of them. These MSs are convex combinations of an original set of K pure strategies, and the payoffs associated to MSs depend on those associated to pure strategies. The set of MSs can be represented by the vertices of a standard simplex in R M ,

∆ M = y ∈ R M : y m ≥ 0, M m=1 y m = 1 , (1)
where every point y is a distribution. On the other hand, given a distribution y ∈ R M there is an associated distribution x ∈ R K which can be found as,

x = Cy, (2)
where the M columns of C represent all the MSs of the predefined set. Given that each MS is a convex combination of the original set of K pure strategies, each column

m = {1, • • • , M } of C has K elements c k,m such that, c k,m ≥ 0, ∀k = {1, • • • , K} and K k=1 c k,m = 1.
Furthermore, a distribution given the original pure strategies can also be represented by a point x inside the simplex in R K ,

∆ K = x ∈ R K : x k ≥ 0, K k=1 x k = 1 , (3)
where again each vertex represents a strategy, in this case, from the set of pure K strategies. As it can be noticed, MSs are points inside the simplex ∆ K , and distributions by MSs y ∈ ∆ M have representations in terms of pure strategies as points x ∈ ∆ K given by eq. ( 2). As it is described in [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], this is useful to define a subset Θ M ⊂ ∆ K of viable distributions x in terms of pure strategies, defined by,

Θ M = x = Cy ∈ R K : y m ≥ 0, M m=1 y m = 1 , (4)
which will be exemplified later in this section. The MSD is described in continuous time by,

ẏm = y m ((C T f (Cy)) m -ḡ(y)), (5)
where f (Cy) := f (x) is a K dimensional vector whose elements are the payoffs associated to pure strategies, given the population distribution y by MSs, and eq. ( 2). Correspondingly, vector g(y) := C T f (Cy) is an M dimensional vector whose elements are the payoffs associated to MSs. Given the nature of the columns of C, element m of g can be interpreted as the expected value of using the probability mix of pure strategies given by the elements of column m. Term ḡ(y) represents the weighted average value of payoffs. When payoff functions for pure strategies are defined as the elements of the Euclidean gradient vector of a given potential function, it can demonstrated that eq. ( 5) represents a different type of gradient vector called Shahshahani gradient. Thus, un der the MSD, y evolves and converges to a local maximum ŷ for initial distributions in its neighborhood [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], [START_REF] Hofbauer | The Theory of Evolution and Dynamical systems[END_REF]- [START_REF] Harper | Information geometry and evolutionary game theory[END_REF].

In order to illustrate these properties, let us consider a potential function 2 , for x ∈ R 3 . Now, let us define the payoffs vector for pure strategies as the euclidean gradient vector of this potential function. The function U (x) has a global maximum in x+ = [26/30, 11/30, 8/30], which is a point in the plane Both the unfeasible maxima and the feasible maximum are shown in Fig. 1(a). This example illustrates a comparison between the Euclidean gradient and the Shahshahani gradient vectors. For several points x in the feasible set Θ 4 , Fig. 1(a) shows the vector fields of both gradients. As it can be observed, the euclidean gradient vector ∇U (x) points to the plane {x 1 + x 2 + x 3 = 1.5} which contains the unconstrained maximum x+ . This happens because the euclidean gradient points in the direction of maximal increase of the potential function U (x). On the other hand, the Shahshahani gradient vector ∇ sh U (Cy) = C ẏ lies on the tangent space T x ∆ 3 = {x 1 + x 2 + x 3 = 0} of the simplex ∆ 3 , which in words means that its tail (i.e. x) and tip (i.e. x + C ẏ) both lie on the plane [START_REF] Harper | Information geometry and evolutionary game theory[END_REF], [START_REF] Lebanon | Metric learning for text documents[END_REF].

U (x) = -(x 1 -26/30) 2 -(x 2 -11/30) 2 - (x 3 -8/30)
{x 1 + x 2 + x 3 = 1.5} ∈ R 3 . However, if the feasible set is constrained to points x ∈ ∆ 3 = {x 1 ≥ 0, x 2 ≥ 0, x 3 ≥ 0, x 1 + x 2 + x 3 = 1}, then the feasible maximum is x * = [7/10, 2/10, 1/10]. If the feasible set is even more constrained to x ∈ Θ 4 = {0 ≤ x 1 ≤ 0.55, x 2 ≥ 0, x 3 ≥ 0, x 1 + x 2 + x 3 = 1.5},
{x 1 + x 2 + x 3 = 1} ⊃ ∆ 3 [15],
The sample points x ∈ ∆ 3 on Fig. 1(a) are obtained from uniformly distributed sample points y ∈ ∆ 4 , by applying eq. ( 2). On Fig. 1(b), it is possible to check the original sample points y ∈ ∆ 4 and the corresponding vector field of the Shahshahani gradient ẏ for the potential function of the example. Fig. 1(c) shows a profile view of the simplex ∆ 3 and the subset Θ 4 . It also illustrates the vector field of ẏ within the convex hull given by the MSs of C, i.e. ẋ = C ẏ.

In Fig. 1(b) it is possible to check 5 different trajectories, for the evolution of the population, starting from different initial distributions y within ∆ 4 . The equivalent trajectories in ∆ 3 obtained with eq. ( 2), are plotted in Fig. 1(c) as well. It is possible to check that trajectories converge to the feasible optimal.

These characteristics are exploited in [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF] in the design of a methodology for the distributed load management of electrical vehicle fleets. In the next section, this methodology and the proposed analogies are summarized.

III. SUMMARY OF THE MSD APPROACH FOR A PEV FLEET LOAD MANAGEMENT

The method proposed [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF] employs multiple analogies that will be summarized an illustrated in this section. In a first place, pure strategies are considered to be all the possible slots of time (hour steps) where PEVs can distribute their energy demand in order to fully charge their batteries. An example is shown on on Fig. 2(a), where a given population is distributed among three environments: 11.1% on the first, 55.5% on the second, and 33.4% on the third. Moreover, slots of time have previously allocated load corresponding to the forecasted demand of the transformer. In this scenario, PEVs load and forecasted load are both interpreted as populations, willing and not willing to migrate (nomad and sedentary populations) among environments (slots of time or pure strategies). An example illustrates these analogies on Fig. 2(b), where a sedentary population and two nomad populations are distributed in six environments (in the fourth environment there is only a portion of the sedentary population). The sedentary population represents the forecasted load profile for K = 6 time slots (hours for instance), while the two nomad populations represent the load of two PEVs. These two nomad population distributions over the six environments, will evolve in order to maximize their weighted average payoffs.

The MSD-MEP approach proposed in [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF] is a multipopulation model, where each single nomad population distribution evolves taking into account the portability of the environments and the payoffs they provide. These payoffs are defined by the trade-off between two measurements, a local load distribution entropy measurement and a total load distribution entropy measurement, given by,

f i k (x i k ) = -α ln   l k µ + J j=1 x j k µ   -(1-α) ln x i k Γ i -1, (6)
where, α i is the trade-off factor (0 ≤ α i ≤ 1) defined by the owner of PEV i, α is the mean of all the trade-off factors, l k is the forecasted load for time slot k, x i k is the load from PEV i allocated to time slot k. On the other hand, the total load µ If λ has changed From inputs define µ, Γ i , and C i End MSD Routine for vehicle i Inputs: λ, K i , soc i des , soc i ini , x i up y i [0] = y i * , from previously saved distribution. Uniformly distributed otherwise.

While V ar y i {h i } > tol or limit of iterations

f i [n] = -1 -α log((C i y i [n] + λ -x i * )/µ) -(1 -α) log(C i y i [n]/Γ i ) g i [n] = C iT f i [n] ḡi [n] = y i [n] T g i [n] y i [n + 1] = y i [n] + Nsy i [n] • (g i [n] -ḡi [n]) x i [n + 1] = C i y i [n + 1]
End Save x i * = x i , y i * = y i , and λ. including both forecast and PEVs, and the total local load Γ i for PEV i, are given by,

Send x i * i λ x i *
µ = K k=1 λ k , Γ i = soc i des -soc i ini τ , (7)
where

λ k = l k + J j=1
x j k is the total load at time slot k, soc i ini and soc i des are the initial and desired states of charge for PEV i, and τ is the length of time slots (hours). It should be noticed that total and individual PEV load distributions lie in simplices like the one of eq. ( 3), given by,

∆ K µ = λ ∈ R K : λ k ≥ 0, K k=1 λ k = µ , (8)
∆ K i Γ i =    x i ∈ R K i : x i k ≥ 0, K i k=1 x i k = Γ i    . (9)
Each individual population distribution evolves such that the trade-off is maximized. The trade-off parameter is defined by owners In the extreme case where α = 1, each local evolution is such that the total load reaches a distributtion as uniform as possible among the environments. In the other extreme case, where α = 0, each local evolution is such that the local load of the PEV is distributed as uniformly as possible among the environments without taking into account the sedentary population, or other PEVs populations. The evolution of each population distribution is locally controlled by each PEV which has an associated MSD routine like the one presented in Fig. 3. Further details on the MSD routine are provided in [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF].

IV. THE MOST CONVENIENT MSS FOR THE APPROACH

In order to illustrate the role of MSs in this approach, let us consider the example of Fig. 2(c). In this example the population of a PEV is going to be split among three predefined MSs. These MSs are possible ways of distributing the PEV load in the available environments, the first MS (m = 1) pre-defines 5.5% of the population allocated in the first environment, 5.5% in the second one, and 89% on the third; the second MS (m = 2) pre-defines 55.5%, 5.5%, and 39% for each environment; and the last MS (m = 3) pre-defines uniformly 33.3% for each environment. These MSs define a matrix C as in eq. (2). In this example, the PEV population is distributed among these three MSs by 11.1% for m = 1, 55.5% for m = 2, and 33.4% for m = 3. This distribution per MSs can be summarized in y = [0.111, 0.555, 0.334] T , and applying eq. ( 2), the total load is divided by 42.55% for k = 1, 14.79% for k = 2, and 42.66% for k = 3, or x = [0.4255, 0.1479, 0.4266] T . In this approach, MSs are defined based on the limitation of the population that can be assigned to an environment, i.e. the maximum instantaneous power allowed by the chargers. Based on the constraints 0 ≤ x i k ≤ x i up , the total number of MSs was defined by,

M i = K i ! (γ i -1)!(K i -γ i )! , (10)
where K i is the total number of pure strategies for PEV i, and γ i is the minimum amount of slots of time needed to allocate the PEV load taking into account the charger limitation of power. As it was shown in [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], this number can be very large depending on the problem. However, in this manuscript the idea of only using the most interesting MSs is further explored.

In order to illustrate the concept of the most convenient MSs in this approach, let us consider a PEV with an energy demand of 8kWh. This PEV has to charge its battery in K i = 6 hours with a 3kW charger. The minimum amount of environments needed is then γ i = 3, and the total amount of MSs that represent completely the feasible set is M i = 60. Nevertheless, some of the most interesting MSs can be selected as those which are more convenient for the PEV owners. For instance, in Fig. 4, for the same example, 7 out of the 60 MSs are plotted. These MSs are those that charge the battery as fast as the charger allows it, and those that apply full charging power during the first hour. The idea behind selecting these MSs consists in diving the total load of each PEV among only the convenient MSs and neglect other MSs that do not benefit the PEV owners. At the end, as it was explained with the example on Fig. 2(c), the final distributions are linear convex combinations of the chosen MSs. Consequently selecting only those MSs having a positive effect for owners, will result in final load distributions more beneficial for them, while grid services are still provided by PEVs as much as possible.

Reduced set of MSs (M i = 4) V. RESULTS

k = 1 k = 2 k = 3 k = 4 k = 5 k = 6 m = 1 m = 4 Reduced set of MSs (M i = 4) m = 3 m = 2
Firstly, 5 PEVs are considered to be charged in a residential grid between 11h and 17h. These PEVs are connected to chargers limited in nominal power to 3kW. In these conditions, they have 6 hours to consume 8kWh in order to reach a desirable state of charge (soc) (initially they have 40% and they have to reach 80% of their storage capacities of 20kWh). Thus the total amount of MSs for each PEV is again M i = 60 like in the example of the previous section. In Fig. 5, it is possible to observe the final load distributions obtained by the MSD approach in three cases: (Case 1) all the 60 MSs are used, (Case 2) only the 7 convenient MSs of Fig. 4 are used, and (Case 3) only the 4 most convenient MSs are used (the first two and the last two MSs in Fig. 4). The mean trade-off factor α is assumed to reach α = 0.99, given certain effective incentive policies from grid managers for PEV owners.

As it can be observed on the upper diagrams on Fig. 5, on each one of the three cases, all the PEVs converge to common final local load distributions. As it can be expected, using reduced sets of MSs still maintains the resource allocation fairness property of the MSD approach. Observing the lower diagrams on Fig. 5 the final total load distribution is almost flat, while for the cases with reduced sets of MSs, the final reached total load is more uneven. This happens because the local load distributions evolve such that they try to reach a maximum total load distribution entropy (a total load profile as flat as possible) with only the available MSs. On the other hand, it can be observed that reduced sets of MSs privilege early high charging rates, specially on the third case where the first three hours show increased charging rates compared to the other cases. Figs. 6 and 7 illustrate the evolution of the convex combinations for the three cases. They show the mapping of the sets of MSs which are convex combinations as well. In these figures, each vertex represents a pure strategy or an environment (a slot of time k), and red points inside represent the corresponding MSs. It is also possible to check how the 5 load distributions evolve and converge to the same distribution. Given that the third case has a set of only four MSs, the 5 load distributions evolution can be more easily illustrated in the simplex ∆ 4 (a tetrahedron in R 3 ) as it is show in Fig. 7.

On the other hand, a realistic scenario is presented, using real active power measurements of a distribution transformer, from the SOREA utility grid company, in the region of Savoie, France. These measurements correspond to mostly office buildings. In this case, the arrival and departure of several PEVs is studied in a random way, following a poisson model with variable rate of arrivals (changing with time) and variable connection times. The highest rate of arrivals is 5 PEVs/h around 05h and it decays up to 0.5PEVs/h at 04h the next day. Vehicle battery capacities are defined to be 20kWh. Chargers are randomly chosen to have limits of power of 3.3kW with a probability of 80%, and 7.5kW with a probability of 20%. For reducing impact on batteries lifespan, states of charge are constrained between 30%-80%. The mean trade-off factor α is assumed to reach α = 0.9, given incentive policies from grid managers for PEV owners.

Fig. 8 shows the total load distribution using full sets of MSs, reduced sets of convenient MSs like those of Fig. 4, and reduced sets of only the most convenient MSs privileging early high charging rates. It is worth noting that the peak of connected vehicles is coincident with the peak of the base load demand, which makes a management indispensable. Despite the degradation observed in the final total load distributions with the reduction of the amount of used MSs, the proposed MSD approach with only convenient MSs, is still able to handle PEVs' load properly. Fairness, PEV owners' convenience, and impact reduction are still achieved. Taking into account that entropy of the total load distribution is given by S(λ) = -K k=1 λ k ln (λ k /µ) [START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF], entropy measurements were taken for each hour step of the first three days of the scenario on Fig. 8. Using horizons of K = 24 hours, the results are presented on Fig. 9 for the ideal case (uniform total load distribution), the unmanaged case, and the managed case using MSD-MEP with: full sets of MSs, Best MSs, and Reduced Best MSs. It is possible to observe that in terms of total load distribution entropy maximization, MSD-MEP with full sets of MSs achieves the closest entropy measurements to the ideal case (uniform distribution). It can be observed as well that the worst measurements correspond to the unmanaged case, as it is expected. The most important thing to highlight is the closeness of the entropy measurements for the MSD-MEP using full and reduced sets of MSs, proving that the the proposed MSD approach with only convenient MSs, is still able to handle PEVs' load properly.

State of charge profiles on Fig. 8 show how owners benefit from the most convenient MSs. In particular, profiles in bold lines (marked with arrows) are useful to check how the use of the most convenient MSs affect charging rates, allowing high consumption as early as possible. Reduced sets of the most convenient MSs allow to take advantage from fast charging infrastructures while the impact on the transformers is still reduced, as it is observed with the entropy maximization measurements on Fig. 9. Additionally, Fig. 10 shows the frequency distributions of slots of time required to overpass 50% of the required energy, for all the PEVs during the four days. It can be observed that these frequency distributions are similar for full sets of MSs and sets of Best MSs, with mean values of 4.96 and 4.92 slots of time, respectively. However, with sets of Reduced Best MSs early fast charging rates are privileged, showing a reduction of the mean value to 3.28 slots of time, and a concentration of more than 70 vehicles (over a total of 178) using only 2 slots of time to overpass the 50%.

VI. CONCLUSION

This paper explores the selection of appropriate MSs in a MSD application for load management of PEV fleets. The MSD approach, using reduced sets of the most convenient MSs, privileging early high charging rates, is still able to handle the PEVs' load properly. Besides, fairness, PEV owners' convenience, and impact reduction are still achieved by the approach. As it has been studied, the most convenient MSs allow to take advantage from fast charging infrastructures while the impact on the transformers is still reduced, compared to the case without management. Furthermore, this wider revision on the MSD approach, the definition and selection of MSs, and the illustration of Shahshahani gradient concepts, motivates the exploration and conception of new applications.

Fig. 1 .

 1 Fig. 1. a) Euclidean and Shahshahani gradient vector fields (scaled) for an example of potential function with domain x ∈ R 3 . The feasible maximum x lies in the subset Θ 4 of the simplex ∆ 3 , while unfeasible maxima x * , and x+ lie in planes {x 1 + x 2 + x 3 = 1}, and {x 1 + x 2 + x 3 = 1.5} respectively, in R 3 . b) Representation of the simplex ∆ 4 ∈ R 4 as a tetrahedron in R 3 , Shahshahani gradient vector field (scaled) in ∆ 4 , and 5 different trajectories converging to the maximum. c) Equivalent Shahshahani gradient vector field, and equivalent 5 trajectories, as seen in ∆ 3 . (+: initial distributions, : maxima)

Fig. 2 .

 2 Fig. 2. a) Example of a distribution of a single population (PEV load) in three environments (pure strategies, or slots of time). b) Example of a multipopulation model, with populations willing to migrate (PEV load) through available environments, and sedentary populations (Base forecast load). c) Example of distribution of a single population (PEV load) among three predefined MSs, and the resulting distribution among pure strategies.

Fig. 3 .

 3 Fig.3. Description of a local MSD routine[START_REF] Ovalle | An electric vehicle load management application of the mixed strategist dynamics and the maximum entropy principle[END_REF]. Operator a • b represents the entry-wise product among vectors a and b. In this manuscript, the trade-off factor α is common to all single populations.

Fig. 4 .

 4 Fig. 4. Reduced set of the convenient MSs, for a PEV with an energy demand of 8kWh, K i = 6 hours, and a 3kW charger. The full set has M i = 60 MSs. ( : Most convenient MSs privileging early high charging rates.)

Fig. 5 .Fig. 6 .

 56 Fig.5. Example of 5 PEVs, all of them with energy demands of 8kWh, K i = 6 hours to charge their batteries, and 3kW chargers. Left: All of them using full sets of 60 MSs. Middle: All of them using reduced sets of only 7 convenient MSs. Right: All of them using reduced sets of 4 of the most convenient MSs, privileging early high charging rates. Upper: final PEV load distributions and state of charge profiles of the 5 PEVs, for the three cases. Lower: final total load distributions, for the three cases, compared to the non-managed case.

Fig. 7 .

 7 Fig. 7. Similar representation to that of Fig. 6, for the evolution of load distributions of the 5 PEVs using the 4 most convenient MSs privileging early high charging rates. (+: initial distributions, : final distributions)

Fig. 9 .

 9 Fig.9. Entropy measurements of total load distributions λ (Forecast + PEVs), for the unmmanaged case, the ideal case (uniform distribution), and for the MSD approach with full sets of MSs, and convenient MSs. Measurements are taken each hour on three of the four study days (with 24 hours horizons).

  , when full sets of MSs are employed, Fig.8. From top to bottom: 1) Profiles of arrivals and connected PEVs throughout the evaluation time. 2) Load forecast (blue dots), PEV aggregated load without management (green +), PEV aggregated load with MSD and full sets of MSs (purple ), and corresponding state of charge profiles (red for 3.3kW, and blue • for 7.5kW chargers). 3) Similar profiles, for MSD with reduced sets of convenient MSs. 4) Similar profiles, for MSD with reduced sets of the most convenient MSs, for early high charging rates. Bold line state of charge profiles: examples of 2 PEVs and the effect of convenient MSs.

		10																	30
	PEV Arrivals	2 4 6 8																Connected PEVs	5 10 15 20 25
		0																	0
		17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h
										Time																	Time
	Active Power [kW]	10 20 30 40 50 60 70 80																SoC [kWh]	6 16 8 10 12 14
		0
		17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h
										Time																	Time
	Active Power [kW]	10 20 30 40 50 60 70 80																SoC [kWh]	6 16 8 10 12 14
		0
		17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h
										Time																	Time
	Active Power [kW]	10 20 30 40 50 60 70 80																SoC [kWh]	6 16 8 10 12 14
		0
		17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h	23h	05h	11h	17h
										Time																	Time
		25
					Uniform		All MSs	Best MSs		Reduced Best MSs		Unmanaged
		24
		23
		22
		21
		20
		17h	23h		05h	11h	17h	23h	05h	11h	17h	23h	05h	11h		17h	23h	05h
										Time

VII. ACKNOWLEDGMENT

This work was supported by the ARC Energies Région Rhône-Alpes, and the Project ANR PARADISE, France.


Kaustav Basu

Ahmad Hably

Vincent Debusschere

Geert Seddik Bacha

Jan Geert

Andres Dirven

Ovalle

Jan Dirven

Andres Ovalle

A comparative study of low sampling non intrusive load dis-aggregation

Keywords: intrusive load monitoring, load disaggregation, Signal Processing, Smart Meters, Energy Management, Smart Grids

I. INTRODUCTION

Smart meters are getting deployed worldwide on a large scale. This large scale generation of digital energy data mandates a deeper look inside the consumption patterns of different appliances present inside the house. The insight into the residential load consumption has advantages for both the energy utility provider as well as for the consumer. In the context of smart grids, insights on the energy consumption patterns through appliance usage helps to manage the energy distribution [START_REF] Strasser | Review of trends and challenges in smart grids: An automation point of view[END_REF], [START_REF] Palensky | Demand side management: Demand response, intelligent energy systems, and smart loads[END_REF], especially for the integration of more fluctuating energy sources (i.e. renewable). The strategies commonly employed are demand response to reduce peak demand by eliminating electricity use, or by shifting it to non-peak times [START_REF] He | Incorporating nonintrusive load monitoring into building level demand response[END_REF], [START_REF] Siano | Demand response and smart grids -a survey[END_REF]. The time of use pricing has also been successfully employed to this end. From the point of view of energy providers, load identification can also play an important part in future prediction of usages of particular appliances where the process of historical data collection is made as less intrusive as possible [START_REF] Basu | A prediction system for home appliance usage[END_REF], [START_REF] Basu | Predicting home service demands from appliance usage data[END_REF]. The current smart meters are capable to report the total electrical load, to further develop energy management strategies, appliance usage insight needs to be provided. The appliances in the residence could be directly monitored but the associated cost and inconvenience to be used are the primary reasons that smart meters are unable to succeed in the market. Non-intrusive methods propose an attractive alternative with reduced cost and manual overheads.

To achieve this goal new data analysis mechanisms have to be proposed to inhabitants for their satisfaction and energy costs reduction. Just a transfer from an analogue to a digital system is not good enough for the customers. A comprehensive and qualitative data analysis mechanism has to be proposed K. Basu, A. Hably, V. Debusschere, S. Bacha and Andres Ovalle are with the Grenoble Alpes University -G2ELab, France. e-mail: ahmad.hably @gipsalab.grenoble-inp.fr GJ. Dirven is the CEO of Greeniant B.V coupled with the subsequent load management strategies. The residential building's sector is mainly considered in this work but the tools and methodology are fully applicable to any other kind of buildings. Next to the efforts necessary on the dataacquisition side, it is equally important to develop necessary insights in consumer needs. The NILM technologies heavily rely on the granularity of data. This work primarily digs deep into the two commonly found sampling rates and highlights the use cases that can be generated from it. The methodology used is data driven and focuses on the fundamental signal parameters that differentiates between appliances.

II. PAST WORK

The NILM can broadly segmented into six data flow phases [START_REF] Birt | Disaggregating categories of electrical energy end-use from whole-house hourly data[END_REF], namely, data acquisition followed by data processing, event detection, feature extraction, event classification and finally energy computation. In the early 90's one of the pioneering work in load disaggregation was published [START_REF] Hart | Nonintrusive appliance load monitoring[END_REF]. It proposed to identify appliances by their respective ON/OFF transitions. It tried to consider power consumption changes (i.e. events) both in the active and reactive power of the signal. The current output of a NILM algorthm, indentfying the various high consuming appliances in the residence is shown in Fig. 1. Methods were then proposed to identify individual appliances from their On/Off transitions. From that time, most of the approaches were event-based and at a high sampling rate, typically less than one second. Even though is considered one field, there are two fundamentally different directions. The main difference is in the granularity of the energy data that is considered, which drives the use of tools and analysis techniques. The two granularities are approximately:

• High sampling rate ≥50 Hz (≤ 0.02 s) (up to several kHz). • Low sampling rate ≤1 Hz (≥ 1 s). Signals are primarily analysed using event and non-event based techniques. All high sampling rate data are generally analysed using event based techniques but in the low sampling case it is analysed by both event based and non-event based techniques. These techniques are discussed in details in the following sections.

1) High sampling rate NILM (event based): Highfrequency measurements allow the analysis of not just the energy used, but also the analysis of the structure of the current-voltage waveform itself [START_REF] Hassan | An empirical investigation of vi trajectory based load signatures for non-intrusive load monitoring[END_REF]. Running different appliances introduce different signatures of disturbances of the pure sine waveform, which can be used to disaggregate them. Getting this high-frequency data requires specializing measuring equipment such as current clamp. In [START_REF] Dong | Non-intrusive signature extraction for major residential loads[END_REF], a method has been proposed to generate signatures based on the active power range, reactive power range, harmonic content range, with or without spike, single phase or double-phase and searching time. Approaches typically consist of identifying the steady state or in some cases transient state features [START_REF] Zeifman | Nonintrusive appliance load monitoring: Review and outlook[END_REF], [START_REF] Li | Building recognition using local oriented features[END_REF]. Subsequently, these signatures are matched with earlier learned models using a pattern recognition algorithm [START_REF] Berges | Enhancing electricity audits in residential buildings with nonintrusive load monitoring[END_REF], [START_REF] Bier | Smart meter systems measurements for the verification of the detection & classification algorithms[END_REF]. The drawbacks of these approaches are mainly hardware requirement due to high sampling rates and the impracticality of the process being totally non-intrusive [START_REF] Fernandes | Load profile identification interface for consumer online monitoring purposes in smart grids[END_REF], [START_REF] Norford | Non-intrusive electrical load monitoring in commercial buildings based on steady-state and transient load-detection algorithms[END_REF]. These methods do not fit well into the smart meter sampling rate, so separate device has to be installed for training, visualisation and communication to the grid. This is a major drawback for these methods, commercially and practically speaking. The load separation at a high sampling rate of all the appliances also raise privacy concerns as user activity can be easily detected, interpreted and monitored [START_REF] Birt | Disaggregating categories of electrical energy end-use from whole-house hourly data[END_REF]. The techniques rely on first identifying the power level transition of different appliances which in the NILM literature is called events. In Fig. 2 the event based technique is shown. Fig. 2a shows the generic appliance On/Off locations which are identified and subsequently matched to detect the appliance.

2) Low sampling rate NILM (event and non-event based): Low-frequency measurements obviously complicate the disaggregation process, but have one major advantage: They are easily available from smart meters from all over the world. Typically, these meters supply whole-home electricity usage every 1s, 10s or more. Since all the information of the waveform is lost, and typically only active power is available (and not reactive power), different techniques are required for disaggregation. The major issue at low sampling rate is that low energy consuming devices are difficult to be detected as the switching events are not very prominent. However, high energy consuming appliances, such as water heater or washing machine can still be identified with reasonable precision even at sampling rate of 15 minutes for example [START_REF] Kalogridis | Privacy for smart meters: Towards undetectable appliance load signatures[END_REF], [START_REF] Prudenzi | A neuron nets based procedure for identifying domestic appliances pattern-of-use from energy recordings at meter panel[END_REF]. In [START_REF] Basu | Nonintrusive load monitoring: A temporal multilabel classification approach[END_REF], [START_REF] Basu | Time series distance-based methods for non-intrusive load monitoring in residential buildings[END_REF] the states of various residential appliances were identified for consumption reading every 10-minutes in France.

The results indicated that high consuming appliances specially water heater could be identified even at that rate. The above methods were based on energy consumption reading above 10-minutes. Considering the constraint of low sampling rate, the differentiation of the methods is directly dependent on the choice of algorithms. Algorithms have further been implemented and tested in the field of load monitoring at 1 or 10 seconds. [START_REF] Ahmadi | Load decomposition at smart meters level using eigenloads approach[END_REF] has proposed a method generating appliance feature using the eigen vectors and during testing time using a pattern recognition to match the signatures. A method partially disaggregating total household electricity usage into five load categories has been proposed at a low sampling rate in [START_REF] Kolter | Energy disaggregation via disciminative sparce coding[END_REF] where different sparse coding algorithms are compared and a Discriminative Disaggregation Sparse Coding algorithm is tested. A featurebased Support Vector Machine classifier accuracy is also mentioned but is not presented. The method of [START_REF] Kolter | Energy disaggregation via disciminative sparce coding[END_REF] is an implementation of the blind source separation problem, which aims at disaggregating mixture of sources into its individual sources. A classic example for this would be the problem of identifying individual speakers in a room having multiple mikes placed at different locations. In the NILM context, the problem is undermined as there is only one mixture and a large number of sources. Another issue using blind source separation is the assumption of no prior information about the sources. On the contrary, in the NILM context, the sources (appliances) do have separate usage patterns which could be used. Nevertheless, blind source separation still remains a promising direction of research in this domain. Temporal graphical models such as Hidden Markov Models (HMM) also have been promisingly used in this domain as they are a classical method for sequence learning [START_REF] Parson | Using hidden Markov models for iterative non-intrusive appliance monitoring[END_REF]. They have been successfully used in many domains, especially speech recognition. In the NILM context, the problem is to learn the model parameters given the set of observations as input sequence and appliances states as output. HMM also considers sequential patterns in consumption but in the NILM context, at a very low sampling rate it seems to have a sensibility to training noise.

[24] proposed a technique using subtractive clustering and the maximum likelihood classifier. The features used to identify appliances are the power levels and the ON/OFF durations. The power levels are computed using a normal distribution and the ON/OFF duration using a Weibull distribution. On average for the six appliances under consideration, the average accuracy is 86%. In the proposed work the technique is enhanced to take into account the various power levels or states at which the appliance operate and consider temporal correlation among them. It is tested for a commercially viable solution considering all the practical challenges encountered during implementation. The remainder of the paper concerns the identification of electrical appliances usages from the smart meter monitoring. The main objective of this work is to compare NILM techniques for both 10 second and 15 minutes data.The proposed algorithm at 10-second sampling is a novel event based algorithm which improves the existing literature by modelling temporal relations between appliances electrical components. This work implements the one proposed in [START_REF] Basu | Nonintrusive load monitoring: A temporal multilabel classification approach[END_REF], (b) The On/Off events generated from raw data. Fig. 2: Event based NILM, both generic and practical example. [START_REF] Basu | Load identification from power recordings at meter panel in residential households[END_REF] for consumption readings at 15-minutes for houses in the Netherlands and training the model on benchmark Dutch houses using plugs. The results obtained from the data at these commonly observed sapling rates are then compared. In Fig. 3 the non-event based method is shown and the electricalsignature of a washing machine at 15-minutes.

III. DATABASE ARCHITECTURE

The database is provided from Greeniant B.V, which is a smart meter analytics company based in the Netherlands. The database consists of more than three thousand residential customers with the corresponding meter reading hosted in cloud based server. The training of the proposed model is built on a subset of the Greeniant B.V database (approx. 20 houses) with plugs for various appliance reading. For the purpose of validation the technology, four new houses were selected with additional plugs for the appliance power reading were installed. The database corresponded to varying houses in terms of consumption level, number of inhabitants and area of the house. The results are the average for appliances across all the houses.

IV. METHODOLOGY

The methodology can be described as a part based data driven approach. The meter reading is a time series from which daily overlapping sub-sequences are generated. The time series, sub-sequence and sliding time windows are described as follows:

Time Series: Ordered set of n real-valued variables T = t 1 , . . . , t n . Sub-sequence: For a given time series T of length n, a sub-sequence C k of T is a sampling of length w ≤ n of contiguous position from T , that is, C k = t k , . . . , t k+(w-1) for 1 ≤ k ≤ n -(k + 1). Sliding Time Window: Given a time series T of length n and a subsequence length of w, a matrix M of all possible subsequence can be built by "sliding time window" across T and placing subsequence C k in the k th row of M . The size of matrix M is (n -k + 1)/n × w. An appliance is identified to be used in a particular day if it starts before the end of the day (24:00). Four hours of additional time overlap at the end of the day is taken for the cross-over cases. The raw training data consists of 10-seconds active power reading together with the energy consumption during the period. The 10-second active power reading is used for the 10-second analytics whereas the 15-minutes consumption reading is used for the 15-minutes analytics.

A. 10-second analytics

The 10-second smart meter analytics consists of the following stages:

1) Preprocessing of the data.

2) Element generation from the ON/OFF events.

3) Appliance-wise filter parameter generation from training database. (Both smart meter and appliance data) 4) Element detection based on trained filter parameters 5) Calculation of start time and energy after appliance detection. 1) Pre-processing: At the preprocessing stage the signal noise such as data spikes and outages are processed and removed.

It is important to remove the spikes as they interfere with element detection. The spikes are removed using a median filter. In the data received from the server if the outage in the data is above a pre-defined threshold then appliance identification during the period is not considered.

2) Element generation from the ON/OFF events: The event is registered when the change in power level is above 35 W. This minimum threshold is chosen to reduce event complexity and to reduce the error of minor fluctuations. In Fig. 4, low power fluctuation seen in some smart meter data is shown.

The delta calculation can be expressed by the following equation:

∇(t i ) = P (t i ) -P (t i -1) ∇(t i ) > 35W t i ∈ [0, w] (1)
where 'w' is the window size. In Fig. 5 the events generated for a typical 10 second smart meter reading is shown. It can be seen clearly that the generated events correspond to various appliance states. The blue color corresponds to 'no appliance indentification' and the other colors are for 'detected appliances'. Element generation: Within the time-window (experimentally set to four hours based on maximum appliance duration), a list of consecutive positive events are sorted, and a list of consecutive negative events are sorted that come after the last positive. One negative event at a time is taken and subsequently matched with a positive event in the list of positive event. The events matched are removed and the process followed iteratively.

3) Appliance signature parameter calculation: Once the elements are generated the subsequent step is to generate its features or signature. The features are evaluated to consider the power levels, duration, number of occurrence and gaps within appliance states.

• Mean and Standard deviation of the ON and OFF duration in the window. • Number of occurrences within the window.

• Mean and standard deviation of the power level.

• Time difference between the various appliance electrical elements. As previously stated the power level based clustering roughly corresponds to various states or elements of a multi-state appliance and helps us to consider the temporal relation between the various states.

4) Training of part based model parameters:

In Fig. 6 it can be seen how the grouped elements roughly can be interpreted as appliance states. Once the features are generated for each time-window, it becomes training instance for the classifier and all the training instances are generated using the sliding time window (typically window jump of 30-minutes). The output or appliance state is taken to be ON if the appliance is 100 % present in the time-window. Once, the training instances are generated, the task is to learn a classifier function which can best map the input to the output (appliance state).

B. 15 minutes NILM

The 15-minutes load dis-aggregations algorithm is trained from the consumption reading at 15-minutes. The data is only pre-processed for outlier removal but no thresholding is performed.

1) The aggregate energy readings are extracted for the 20 houses at a sampling rate of 15 minutes. 2) Multiple sub-sequences and one class output are generated from each data-set using temporal sliding window with a window size of 6 units.

3) The classifier is trained with the generated features as input attributes and high energy appliances as output classes [START_REF] Basu | Nonintrusive load monitoring: A temporal multilabel classification approach[END_REF].

4) The model is trained from the 20 houses for a period of 4-weeks. A simplified graphical work-flow of this approach is illustrated in Fig. 7 from step one to four. The features generated step is implemented from the [START_REF] Basu | Nonintrusive load monitoring: A temporal multilabel classification approach[END_REF] article but the sampling rate are different.

V. RESULT

In this section, firstly, the measure used to evaluate the performance is discussed followed by the observed performance metrics using two different sampling rates.

A. Measures

Given a dataset of labelled instances, supervised machine learning algorithms seek a hypothesis that will correctly predict the class of future unlabelled instances. In order to compare structures of predictors, we need indicators that will give a quantitative way of assessing the classifier performances. While comparing these indicators values, the best predictor can be found for a given appliance.

In order to properly define the performance indices of the classification algorithms used for prediction, we introduce the confusion matrix [START_REF] Kohavi | Glossary of terms[END_REF].

A confusion matrix contains information about the actual and the predicted results obtained by a classification system. The performances of such systems are commonly evaluated using the data contained in this confusion matrix. Table I shows the confusion matrix for a two-class classifier. The classes that can be predicted are "positive" or "negative" instances, which in this case signifies the appliance consumes or does not consume energy. In the context of this study, the entries defined in the confusion matrix reported in Table I have the following meaning: a : is the number of correct predictions where an instance is negative, b :

is the number of incorrect predictions where an instance is positive, c :

is the number of incorrect of predictions where an instance negative, d :

is the number of correct predictions where an instance is positive. Several standard terms have been defined for this 2 class matrix:

The true positive rate (recall): is the proportion of positive cases that were correctly identified, as calculated using the equation T P = d/(c + d).

Represents the ratio between the predicted positives states of the appliances (ON) and the total number of correct positives states of the appliances. The precision: is the proportion of the predicted positive cases that were correct, as calculated using the equation P = d/(b + d). Represents the fraction of the positives states (ON) of the appliances correctly predicted. A correct appliance identification is defined as being one where the start time is +/-thirty minutes of the actual start time and for the correct appliance. The false-positive impact for use-cases is much more than the impact of false-negative thereby the precision scores are presented.

B. Dis-aggregation performance

In Table II the appliance-wise dis-aggregation performance is shown for the major white goods appliances present in a typical Dutch residence. The washing machine is present in all the houses and the electric oven in just one as in most residence it is run on gas. It can be observed from the results that the performance is better for the washing machine and dish washer, it reduces significantly for the clothes dryer and particularly electric oven. In case of 15 minutes sampling rate, the washing machine and dish-washer are considerably reduced and are mostly detected during the non-peak times. The major cases of error are also the clothes drier vis-a-vis Oven error due to mis-labelling. This can be improved by taking context into account, e.g. a Drier is generally used after the washing machine.

C. Discussion

The proposed method is an industrial application at 10second and 15-minutes sampling showing the pros and cons of using load dis-aggregation in residences. The existing literature is based on acquired dataset but not actual customers, which is a challenge in itself, addressed in this work. This method is an implementation of existing methods enhanced by the data-processing and a large training database. This work also emphasizes the gap between the trained models in a benchmark dataset and the intricacies of an application to actual people who have different appliances combinations and usages.

VI. CONCLUSION

NILM has acquired considerable amount of interest in the smart-energy domain and this work presents a comparison of results for two sampling rate for various appliances present in a typical house. It is observed that the performance is correlated to the number of appliances being present in the house as the signal of appliances other than the one which is being identified is considered noise.

The results also highlight the fact that performance of load dis-aggregation decrease as we move from the laboratory setting to the real residences. This may happen due to error at various stages; data acquisition, inter-appliance similarity, presence of previously unobserved appliances etcetera. The recording of household meta-data (appliances present in the house) could be used to reduce mis-labelling among appliances and thereby increase the performance. Finally, the results indicate that the proposed 10-second NILM algorithm is deployable for majority of houses (> 80 %) with a use case for the customer. In case of 15 minutes consumption data, the appliances detection are limited to the non-peak periods of the day.

VII. ACKNOWLEDGEMENT

This work has been a joint collaboration by the project PARADISE which is sponsored by the National Research Agency (ANR), France and the smart meter analytics enterprise Greeniant B.V operating from the Netherlands.

Fig. 1 :

 1 Fig. 1: Non-Intrusive Load Monitoring (blue: no appliances detected).

  (a) NILM event based (generic).

  (a) NILM non-event based (generic). (b) The 15-minutes smart-meter and washing machine consumption.

Fig. 3 :

 3 Fig. 3: Non-event based NILM, both generic and practical example.

Fig. 4 :

 4 Fig. 4: 32 W fluctuation in the meter-data.

Fig. 5 :

 5 Fig. 5: Events generated from the three consecutive use of washing machine (In red).

Fig. 6 :

 6 Fig. 6: The generated elements relation to appliance states.

Fig. 7 :

 7 Fig. 7: Methodology: 15-minutes load disaggregation.

TABLE I :

 I Confusion matrix

			Predicted
			Negative Positive
	Actual	Negative Positive	a c	b d

TABLE II :

 II Dis-aggregation performance, Precision (measure).

	Sampling Rate		Appliances
		Washing	Clothes	Dish	Electric
		Machine	Drier	Washer	Oven
	10 second	.89	.64	.79	.50
	15 minutes	.5	.6	.5	.32


Khaled Hajar
email: khaled.hajar@gipsa-lab.grenoble-inp.fr

Ahmad Hably

Seddik Bacha

Ahmad Elrafhi

Ziad Obeid

An application of a centralized model predictive control on microgrids

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

A small-scale power supply network where small communities are provided by their own power is named microgrid (MG). A MG contains loads, distributed generation (DG), renewable energy sources (RES) and energy storage systems (ESS) like electrical vehicles, hydraulic storage, batteries, etc. Microgrid applications have different objectives and characteristics but they all share some common technical problems due to the use of RES. These RES have an intermittent, random, non programmable nature, which makes them difficult to provide a continuous power supply to isolated loads and to assure energy dispatch of RES in distribution grids [START_REF] Hajar | Optimal centralized control application on microgrids[END_REF]. MGs can be considered as intelligent distribution systems with two different modes of operation: the isolated mode and the grid-connected mode [START_REF] Al-Nasseri | Protecting micro-grid systems containing solid-state converter generation[END_REF]. However, it is also expected to provide sufficient generation capacity, controls and operational strategies to supply at least a portion of the load, after being disconnected from the distribution system and remains operational as an autonomous (isolated) entity [START_REF] Shah | Decentralized power flow control for a smart micro-grid[END_REF]. The general objective of energy management of a microgrid is to minimize operating costs such as fuel, maintenance, and the purchased cost of exchange power from the main grid [START_REF] Yue Xiang | Robust energy management of microgrid with uncertain renewable generation and load[END_REF]. As for the distribution system operators (DSO), the presence of MGs provides many advantages. The DSO and MGs owners can benefit from the lower operation costs. The customers can benefit from a more reliable and economical power supply. Therefore, it is necessary to consider the DSO and networked MGs altogether [START_REF] Ahmed Ouammi | Coordinated model predictive-based power flows control in a cooperative network of smart microgrids[END_REF]. There are several works that address MPC and optimization of a network of MGs. In [START_REF] Venkat | Distributed mpc strategies with application to power system automatic generation control[END_REF], an MPC has been proposed where controllers are supposed to operate separately but their design is based on a global model of the overall system.

In [START_REF] Falahi | Decentralized energy management system for networked microgrids in grid-connected and islanded modes[END_REF], the authors have proposed an MPC for the optimal power exchanges in a smart network of MGs. Falahi et al. [START_REF] Zhaoyu Wang | Dynamic model predictive-based energy management of dg integrated distribution systems[END_REF] proposed an MPC based energy management system to regulate the active and reactive power in a microgrid. An investigation on decentralized linear quadratic Gaussian control (LQG) approach has been made in [START_REF] De Brabandere | A voltage and frequency droop control method for parallel inverters[END_REF].

In this present paper, we apply a centralized MPC of the power flow among a network of MGs and with DSO. This study is based on the economic profit of all the actors of the Grid (DSO and all the MGs) while minimizing the exchange between MGs and DSO from one side, and maximizing the use of renewable energy in each MG and the exchange of power between neighboring MGs from the other side. The paper is organized as follows. In Section II, MG units are presented. MG modeling as well as MPC problem are developed in Section III. The results of the application of the methodology are presented and commented in Section IV. Section V provides some conclusions and perspectives.

II. MICROGRID UNITS

Several models for MG have been proposed in previous studies. We have used the modeling procedure proposed by [START_REF] Del Real | An integrated framework for distributed model predictive control of large-scale power networks[END_REF] to create our model of MG. Energy storage systems (ESS), photovoltaic (PV), wind power, local loads and distributed generation are integrated in each MG. The variation of wind speed and solar irradiance affect the wind power and PV efficiency respectively. Depending on the loads and the production ESS has two status: charging and discharging. In our previous work [START_REF] Hajar | Optimization of a microgrid with renewable energy and distributed generation: A case study[END_REF], a control unit (CU) is used in order to maximize the power stored in ESS to assure the local consumption. Global Management Controller (GMC) is notified of the state of CU which optimize the control of each MG. It verifies if the MG has excess of power or viseversa. For each MG the goal can be achieved by controlling the charge/discharge of ESS, starting or stopping DG, and receiving/sending power from/to neighboring MGs. The GMC is dedicated to manage the power in the grid, the relation between MGs as well as with the DSO with a predicted control. It takes also the decision wherever the excess power is to be sent to other MGs or to DSO. The different elements of MG are presented in Fig. 1 from our previous work [START_REF] Hajar | Optimal centralized control application on microgrids[END_REF]. This figure also shows the units of producing or consuming energy. A microgrid can also be connected or disconnected from DSO or/and can be connected to other MGs. To summarize, an MG integrates:

• DSO providing/purchasing power for/from the MGs.

• Loads representing consumers demand.

• RES providing renewable energy.

• ESS improving the quality, and reliability of supply.

• DG improving the reliability of supply.

• Control Unit (CU) transferring the data to GMC.

III. MICROGRID MODELING AND OPTIMIZATION

A. MPC algorithm steps

The MPC is a control strategy that is suitable for optimizing the performance of constrained systems. Constraints are present in all control systems. Through handling of constraints, MPC can improve the performance of a system by allowing it to safely operate near constraint boundaries. To apply MPC, some steps must be followed. At the first step of time t, the system actual state is taken as the initial point of the problem as the RES, ESS, Loads, power generation. Then, for the control horizon N c we compute an optimal control sequence for the next prediction period N p based on RES, ESS, Loads, power generation in addition to the prediction of energy price. After that the first control step of all MGs will be applied. Finally, an update is made for all the information for next time step while moving to the next sampling time to re-apply the same optimization.

B. Modeling

Figure 2 shows the interaction among MGs on one hand and between DSO and MGs on the other hand. The power flow is as the following:

• P S,D : Predicted power sold by the DSO (kW).

• P S,M : Predicted power sold by the MG (kW).

• P B,D : Predicted power bought by the DSO (kW).

• P B,M : Predicted power bought by the MG (kW).

• P eq : Power equilibrium. The stored energy equation in each MG:

x(m, t + k) = x(m, t + k -1) - (1) δ dis,m E dis (m, t + k) + δ char,m E char (m, t + k)

C. Optimization Formulation

A nonlinear centralized optimization is used to formulate the energy management in MG. Our objective is to optimize the grid in a predictive way that the power generation, power balance, ESS, energy exchange between MGs are satisfied under certain constraints that will be presented later. The cost function F to be minimized is given by the following equation: Fig. 3: MGs power equilibrium.

F = Nc k=1 M m=1 M n,n =m [ψ M (k)P S,M (m, n, t + k)T S,M (n, t + k) (2) -ψM (k)P B,M (m, n, t + k)T B,M (n, t + k)]

D. Constraints

First, we must consider that the power sold by the m-th MG is equal the power bought by the n-th MG: An upper and a lower bound (x m,min , x m,max ) constrained the stored energy in each ESS. In the m-th MG:

P B,M (m, n, t + k) = P S,M (n, m, t + k)
x m,min ≤ x(m, t + k) ≤ x m,max
The charging power must be lower than certain limited power value P charmax,m :

0 ≤ P char (m, t + k) ≤ P charmax,m MG1
0 ≤ P dis (m, t + k) ≤ P dismax,m
At last an MG cannot buy and sell energy at the same time. buying can be done when the MG is in need of power otherwise MG can sell power when its production is higher than its consumption:

• If ∆P M G,n (t + k) > 0 then P B,D (m, t + k) = 0 and P B,M (n, m, t + k) = 0 • If ∆P M G,n (t + k) < 0 then P S,D (m, t + k) = 0 and P S,M (n, m, t + k) = 0 with ∆P M G,n (t + k) = (P EP,n (t + k) -P EL,n (t + k)), Where • ∆P M G,n : Power equilibrium of n-th MG at (t + k)
instant. • P EP,n : Power of the expected production of n-th MG at (t + k) instant. • P EL,n : Power of the expected load of n-th MG at (t + k) instant.

E. Numerical values

The values used in simulation are given in Table I. They give global values for some power constraints and some characteristics of ESS. The maximal capacity of the ESS available in each MG is comprised between 20 and 100 kWh. For each MG, the maximal charge and discharge powers are limited. As example, MG5 can discharge in 2 hours (δ dis,m = 0.5) and need more than 3 hours to be full charged (δ cha,m = 0.3).

IV. SIMULATION RESULTS

The above described optimization problem has been implemented in Matlab. We have used a modified version of the study case presented in [START_REF] Tsikalakis | Centralized control for optimizing microgrids operation[END_REF]. The simulation will take in consideration five interconnected MGs with N c = 3 and N p = 24. When the supply in a MG is less the demand, the solution is to discharge the ESS or purchase power from DSO or many neighbors MG. Each MG generates power and feeds a group of consumer in an interval time (t, t + 1). Also, it is connected to all neighboring MGs. As for the tariff vector of selling and buying energy, (T B,M , T S,M , T S,D , T B,D ) are represented for each MG in a random value between 0.07e/kWh and 0.2e/kWh at day time, and 0.075 e/kW from midnight to 6am (Fig. 4). In Fig. 3, we show the predicted value of P eq for each MG. The value is calculated based on the load and production curves of [START_REF] Tsikalakis | Centralized control for optimizing microgrids operation[END_REF] where we made some random calculation to create more MGs. The ESS can go from 25% to 90% of power storage. The power of ESS in each MG is different than the other. The sold and bought power from/to DSO is limited to kWh for each MG. optimization algorithm searches to use the surplus power available in each MG, under constraints, and to be exchanged. If the constraints are not satisfied, the DSO interferes while respecting its constraints. The optimal results shown in Fig. 5 gives us an overview about the interaction of MG1 with the other MGs. As result of selling and buying, MG1 sells 396kWh in total especially 102kWh to MG2 and 94kWh to MG3, 35kWh for MG5 and finaly to MG4 160kWh. In addition, MG1 buys 112kWh from MG2, 99kWh from MG4,50kWh from MG3 and 64kWh MG5. So in total, MG1 buys 391kWh.

Fig. 6 shows us the optimal energy exchanges of MG2. The total energy bought is 405kWh, the main bought energy is from MG5 with 150kWh, the rest of bought energy comes from MG1 and MG3 with 105kWh and 107kWh respectively.

Optimal energy strategy for MG3 is shown Fig. 7. It is obvious that MG3 buys energy more than selling. In total, it buys 563Kwh from other MGs where the big amount comes from MG4 with 190kWh, from MG5 141KWh, 109KWh from MG2, and 114kWh from MG1. It sells an amount of 250KWh in total, distributed between MG2 and MG4 with 90kWh and 119kWh respectively. For MG4 (Fig. 8), the sold energy represents 422Kwh where MG3 buys 192Kwh, MG1 buys 107kWh and finally MG5 and MG2 buy 58kWh for each one. Total bought energy is purchased from MG1 with 168kWh, MG2 68kWh, MG3 141kWh and MG5 188kWh. MG5 represents the maximum benefit of sold/bought energy to other MGs with a total 668/57kWh (Fig. 9). A. Discussion on the results

ESS charge and discharge of MG1 on Fig. 12 shows between the hour 20 and 23 a brutal discharge and charge respectively. The discharge is due to the need of power at this time in MG1 and when the demand is covered by MG1 production unit or by purchasing power from neighboring MGs, the CU gives a command to charge the ESS of MG1. The time of charge and discharge is the same as δ dis,1 = δ cha,1 (see Table I). On the same graph, it i s clear that the ESS value couldn't exceed the 90% of its total value, nor going beyond 25%. A similar analysis can be done for other microgrids.

V. CONCLUSION

The paper presents a model predictive control in a centralized manner to optimize the interaction between a network of MGs and DSO. Results have shown that the control performance of a network of MGs is achieved by maximizing the benefits of all the constituted elements of network. It is also shown the control of the ESS discharging/charging, and the predicted interaction of power exchange between MGs. The cooperation between MGs demonstrates its importance due to its benefit according the selling of the excess of power produced or its operation in standalone mode. A decentralized MPC version of this algorithm is under study. This control strategy can be used in case of failure of communication between microgrids.

VI. ACKNOWLEDGEMENT

Fig. 1 :

 1 Fig. 1: Microgird elements.

Fig. 2 :

 2 Fig. 2: Power flow in the grid (MG+DSO).

  where • x: Stored energy (kWh). • δ dis,m : Discharging efficiency of the ESS in the m-th MG. • δ char,m : Charging efficiency of the ESS in the m-th MG. • E dis : Discharged energy (kWh). • E cha : Charged energy (kWh).

  (k)P S,D (m, t + k)T S,D (m, t + k) -ψD (k)P B,D (m, t + k)T B,D (m, t + k)] with • m, n: MG index. • M : MGs number. • t : Control step. • N c : Control horizon. • k: Prediction step. • ψM (k): Scale factor of buying power from the MGs. • ψ M (k): Scale factor of selling power to the MGs. • ψD (k): Scale factor of buying power from DSO. • ψ D (k): Scale factor of selling power to DSO. • T S,D : Expected tariff value of selling energy by the DSO (e/kWh). • T S,M : Expected tariff value of selling energy by the MG (e/kWh). • T B,D : Expected tariff value of buying energy by the DSO (e/kWh). • T B,M : Expected tariff value of buying energy by the MG (e/kWh).As described in Fig.2, both terms in the objective function are related to energy sold or bought. The first term concern neighboring MGs and the second one concerns the DSO. In the following section we introduce the constraints for this optimization problem.

Fig. 4 :

 4 Fig. 4: The tariff vector.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Predictive control results for MG1.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Predictive control results for MG3.

Fig. 9 :Fig. 10 :

 910 Fig. 9: Predictive control results for

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Total predicted power exchange between MGs and DSO.

TABLE I :

 I Numerical value used for simulation Also the discharging power must be lower than certain limited power value P dismax,m :

This work has been a joint collaboration by the project PARADISE which is sponsored by the National Research Agency (ANR).