Nicolas Gast
email: nicolas.gast@epfl.chjean-yves.leboudec@epfl.chdan-cristian.tomozei@epfl.ch

Jean-Yves Le Boudec

Dan-Cristian Tomozei

Impact of Demand-Response on the Efficiency and Prices in Real-Time Electricity Markets *

Keywords: C.4 [Performance of Systems]: Modeling techniques., G.3 [Probability and Statistics]: Stochastic processes Electricity Pricing, Demand-response, Market Efficiency, Energy Economics

We study the effect of Demand-Response (DR) in dynamic real-time electricity markets. We use a two-stage market model that takes into account the dynamical aspects of generation, demand, and DR. We study the real-time market prices in two scenarios: in the former, consumers anticipate or delay their flexible loads in reaction to market prices; in the latter, the flexible loads are controlled by an independent aggregator. For both scenarios, we show that, when users are price-takers, any competitive equilibrium is efficient: the players' selfish responses to prices coincide with a socially optimal policy. Moreover, the price process is the same in all scenarios. For the numerical evaluation of the properties of the equilibrium, we develop a solution technique based on the Alternating Direction Method of Multipliers (ADMM) and trajectorial forecasts. The forecasts are computed using wind generation data from the UK. We challenge the assumption that all players have full information. If the assumption is verified, then, as expected, the social welfare increases with the amount of DR available, since DR relaxes the ramping constraints of generation. However, if the dayahead market cannot observe how elastic loads are affected by DR, a large quantity of DR can be detrimental and leads to a decrease in the welfare. Furthermore, the DR operator has an incentive to under-dimension the quantity of available DR. Finally, we compare DR with an actual energy storage system. We find that storage has a faster responsetime and thus performs better when only a limited amount is installed. However, storage suffers from charge-discharge inefficiency: with DR, prices do concentrate on marginal cost (for storage, they do not) and provide a better welfare.

INTRODUCTION

Electricity markets are developing worldwide, replacing tightly regulated systems by decentralized control mechanisms [START_REF] Wilson | Architecture of power markets[END_REF]. Demand-response (DR) mechanisms have emerged in order to take advantage of flexible demand in a system where demand is traditionally considered as inflexible [START_REF] Spees | Demand response and electricity market efficiency[END_REF]. DR can compensate mismatches between production and consumption, thus increasing the potential deployment of renewable generation by reducing the necessity for regulation reserve. It can be implemented by having consumers react to prices [START_REF] Mohsenian-Rad | Optimal residential load control with price prediction in real-time electricity pricing environments. Smart Grid[END_REF] or to congestion signals [START_REF] Meyn | Ancillary service to the grid from deferrable loads: the case for intelligent pool pumps in Florida[END_REF].

The authors of [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF] present a model of real-time electricity market that incorporates the dynamic constraints of generation and the uncertainties dues to forecast errors. The system is composed of two market actors (a supplier and a consumer) that exchange energy at a spot price. They show that there exists a competitive equilibrium, i.e., a price process such that market actors agree on the quantity of energy exchanged. This market is efficient: for any competitive equilibrium, the selfish decisions of both market actors coincide with the hypothetical decisions taken by a social planner that aims to maximize the sum of every actor's utility. However, the equilibrium price exhibits large fluctuations and is never equal to the marginal production cost.

We extend [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF] to model a two-stage electricity market that features generation constraints, inelastic loads, and elastic loads corresponding to a set of DR appliances. Our elastic load model captures several key-features of controllable thermostatic loads. The loads are elastic in the sense that it is possible to perform price arbitrage by delaying or anticipating the consumption of each load, but an applicance's consumption cannot be arbitrarily reduced by augmenting prices. Each appliance has an internal state representing, e.g., its temperature. A load that is delayed for too long reaches an undesirable state. For thermostatic loads, this corresponds to having the temperature outside a given deadband. We add a second dimension to the state space: a counter indicates the waiting time before the appliance is allowed to respond to DR signals, thus modeling mini-cycle avoidance (an appliance that has just been switched on cannot be switched off immediately). In our numerical evaluations we use a mean-field model to keep track of the empirical distribution of the appliances' states.

We consider two market scenarios. In Scenario 1, the flexible appliances are controlled by the consumer who reacts to real-time prices. She takes advantage of price difference by anticipating or delaying the consumption. In Scenario 2 a DR operator controls all the flexible appliances. She charges a fixed rate to users and makes a profit by price arbitrage: scheduling the appliances when the market price is low. We show that in both cases, the market is efficient, i.e., the selfish decisions of players coincide with the ones of a hypothetical social planner.

We develop a numerical methodology based on the Alternating Direction Method of Multipliers (ADMM) to compute the market equilibrium for a realistic forecast error model. The forecast errors are modeled using multiple trajectories [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF], generated from multivariate Gaussian random variables. The correlation matrices are derived from traces of data from the UK.

We study numerically the system as a function of the DR capacity, which is the agregate maximal power of flexible appliances. Our results show that, when the DR capacity is low, the gain obtained by using demand-response is almost linear in the capacity. For large capacity, if the state of the appliances can be observed by all players, the gain saturates to a value that does not depend on the appliances' nature. However, if the states cannot be observed by the day-ahead market, increasing the capacity of DR can be detrimental: after a certain capacity, it decreases the social welfare, as the consumption of DR appliances is harder to predict. Finally, we also show that DR operators have an incentive to undersize the DR capacity.

We also compare DR with classical energy storage [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF]. We find that, at low capacity, a storage system offers better performance than DR because it reacts faster. However, for large capacities of installed DR, DR behaves similarly to a perfectly efficient storage. As such, for large capacities, DR outperforms energy storage systems, since these have a charge-discharge efficiency of 70 -90%.

Roadmap. We discuss related work in Section 2. We present the model, the assumptions, and the main definitions in Section 3, and the social welfare theorem in Section 4. We describe our statistical model of error and our numerical methodology in Section 5. Finally, we give the numerical parameters and present the numerical evaluation of the model in Section 6. We conclude in Section 7.

RELATED WORK

The authors of [START_REF] Albadi | A summary of demand response in electricity markets[END_REF] present an overview of DR programs and their integration within the electricity market. They compare price-based and incentive-based programs and present their benefits (peak reduction, bill savings, reliability) and drawbacks (deployment cost, rescheduling, metering infrastructure). New market models tailored for DR have also been investigated. For example, the authors of [START_REF] Su | Quantifying the effect of demand response on electricity markets[END_REF] develop a new market-clearing mechanism for load-shifting. However, it assumes perfect forecast and is limited to a small number of players. Hence, a key research topic is the means of presenting DR as a service transparent to use via smart home controllers or by having a DR operator with direct control of appliances. Multiple models are proposed, using response to prices [START_REF] Mohsenian-Rad | Optimal residential load control with price prediction in real-time electricity pricing environments. Smart Grid[END_REF] or to congestion signals [START_REF] Meyn | Ancillary service to the grid from deferrable loads: the case for intelligent pool pumps in Florida[END_REF].

Demand response fatigue and rebound effect are a large concern for the DR operator and the grid regulator, in particular because it leads to observability problem. In [START_REF] Boudec | Stability of a Stochastic Model for Demand-response[END_REF], the authors use a simplified model and find that when the state of flexible appliances is not taken into account in the control, a large accumulated delayed load may manifest unexpectedly and randomly. In this paper we model the internal states of each appliance in detail and, in Section 6.3 study whether such an effect continues to occur.

Our appliance model is similar to [START_REF] Koch | Modeling and control of aggregated heterogeneous thermostatically controlled loads for ancillary services[END_REF][START_REF] Meyn | Ancillary service to the grid from deferrable loads: the case for intelligent pool pumps in Florida[END_REF], who also consider thermostatic controllable loads. A difference with our model is that we add a second dimension for modeling minicycle avoidance. We also take into account the undesirable states that could be reached in blackout conditions. Moreover, we consider that the elastic loads can be anticipated or delayed, but that in the long run they consume a fixed amount of energy, e.g., the average consumption of a heater or cooler. This is not the case in [START_REF] Meyn | Ancillary service to the grid from deferrable loads: the case for intelligent pool pumps in Florida[END_REF]. In many papers about demand-response in electricity markets, it is often considered that higher prices result in reduced consumption [START_REF] Albadi | A summary of demand response in electricity markets[END_REF].

Our work builds on various papers about efficiency in electricity markets (e.g., [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF][START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF] and the reference therein). One key contribution of the current paper is that we are able to handle a more realistic forecast model in our numerical methodology. The numerical results produced in [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF][START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF] were performed by assuming that the forecast error can be represented by a Brownian motion with stationary increments. We use a non-stationary model of error that is becoming a standard for wind forecast. It uses a probabilistic forecast developed by Pinson et al. [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF], where the errors are represented by a finite number of possible trajectories.

We consider a system with a large number of flexible appliances. In order to keep the model tracktable, we use a mean-field approximation to approximate their dynamics. More details about these models and their relations with stochastic optimal control are presented in [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Gast | A mean field approach for optimization in discrete time[END_REF][START_REF] Gast | Mean field for markov decision processes: from discrete to continuous optimization[END_REF].

SYSTEM MODEL

We consider a two-stage electricity market (day-ahead and real-time) with three actors: the consumer, the supplier, and the DR operator. The supplier has two roles: she generates the bulk of the electricity determined via the dayahead scheduling and she provides regulation electricity to compensate for mismatches in real-time. The consumer aggregates the consumption of a population of end-users. The controllable demand of end-users is managed by the DR operator who uses the flexibility of the controllable appliances as a means to perform arbitrage on the real-time electricity market. The role of the DR operator can be assumed by the consumer, or by an independent actor. In Section 6, we will explore two types of controllable loads: fridges and boilers.

Two-stage Electricity Markets

Two-stage electricity markets are used to determine the price of electricity supplied to the consumers. In the dayahead stage, the actors use forecasts of consumption and of renewable production to schedule the bulk production for the next day. In the real-time stage, last-minute decisions are taken to compensate mismatches.

We start from the day-ahead market of [START_REF] Cho | Efficiency and marginal cost pricing in dynamic competitive markets with friction[END_REF]. The forecast demand for the next day at time t is used to set the scheduled production g da (t). The scheduled production incorporates both renewable energy (volatile) and conventional energy sources. The forecast demand has two components: the non-controllable demand d da (t) and the controllable (flexible) demand f da (t). An additional fixed quantity r da is produced as a precaution against forecast errors. Thus, in the day-ahead market the actors agree to trade an amount g da (t)+f da (t) the next day at time t at a price p da (t), where g da (t):=d da (t)+r da .

The real-time market deals with the inevitable mismatches that arise from forecast errors. Thus, at time t the endusers express a total aggregated non-controllable demand D a (t) = d da (t) + D(t) and controllable (flexible) demand F a (t) = f da (t) + F (t). The quantities D(t) and F (t) are the real-time components of the demand, and they can be positive or negative. While D(t) is given by nature, F (t) can be controlled to some extent, depending on the state of the controllable appliances and on the control signal decided by the DR operator.

The supplier deploys real-time production G(t) in order to compensate for mismatches. This part of the generation comes from conventional sources, and is subject to ramping constraints. The total amount of produced electricity at time t is thus G a (t) = f da (t) + g da (t) + Γ(t) + G(t), where Γ(t) is the forecast error of renewable sources (i.e., the difference between actual and forecast production). The energy produced in real-time is traded at price P (t) at time t.

In this paper, we study the effect of demand-response on the real-time market. We assume that the market actors are price takers: they strategically define their actions based on prices, but they cannot influence these prices. Moreover, we assume that the market actors base their decisions at time t on the knowledge of past data and actions (up to time t) and the statistics of future data. This knowledge is shared by all actors, unless otherwise specified. In our mathematical model, this means that there exists a filtration (Ft) such that the decision processes are adapted to this filtration.

Controllable Appliance Model

We consider N appliances that can be controlled via demandresponse. Their maximum aggregated power consumption is denoted Pon. At any time step, each appliance is fully characterized by its state i = (s, x, y) ∈ M. The state space is defined as M = {on, off} × Z × {0, . . . , Ymax}, where

• s ∈ {on, off} indicates whether the appliance is on, in which case it consumes a power Pon/N , or whether it is off, in which case it consumes 0; if s = on then we denote by s = off and vice-versa; • x ∈ Z reflects the internal state of the appliance, e.g., the temperature in the case of a fridge; • y ∈ {0, . . . , Ymax} represents the amount of time that needs to pass before the device can react to demandresponse. The device can be switched from on to off (or vice-versa) only if y = 0. Moreover, as soon as the state is switched as a result of a demand-response action, the value y jumps from 0 to Ymax. By accounting for y, we avoid operation in mini-cycles which might damage the appliance. In Figure 1, we represent the Markov chain that gives the evolution of the state of an appliance. Each node on the diagram corresponds to a possible state, with

• s = on for the nodes in the upper rectangle and s = off for the ones in the lower rectangle, • x is given by the position on the x-axis, and • y is given by the corresponding y-axis. Any state i = (s, x, y) has a horizontal transition to a neighboring state (on, x + 1, y) if s = on and (off, x -1, y) if s = off. The probability of this horizontal transition is The x-axis represents the internal state (e.g., temperature). The y-axis is the time before the device can be switched on or off.

γs if y > 0, or (1 -ai)γs if y = 0. We consider that γon is typically larger than γ off . This allows to model the different speed at which the internal state of the appliance evolves when it is on or off (e.g., a fridge spends less time on than off). When y > 0, we add a vertical transition to a state (s, x, y -1), which decreases the y counter with a certain probability θs.

We partition the state space into "responsive" states M1

def. = {(s, x, y) : y = 0} and "non-responsive" states M2 = M \ M1. A demand-response action is defined as a vector of probabilities a = (ai)i∈M indexed on the states. For a state i = (s, x, y), ai is the probability that an appliance switches its s state:

• If the appliance is in the responsive state (s, x, 0), its state becomes the non-responsive state (s, x, Ymax). • If the state is non-responsive, the action has no effect. We consider that if the same action a is sent to two appliances, their transitions are independent. Moreover, the effect of the actions is applied after the horizontal transitions. The actions are represented as thick dotted transitions in Figure 1.

The appliance remains in the same state with the probability that none of the transitions previously described occur: 1 -γs if y = 0 and 1 -γs -θs if y > 0. These self-transitions are not represented in Figure 1 for the sake of readability. The transition probability of an appliance from state i to state j can be written:

π(i, j, α) = P (m(t+1) = j|m(t) = i, ai(t) = α) = (1 -α)π(i, j, 0) + απ(i, j, 1). ( 1
)
Some states of the appliance are undesirable, e.g., if the internal temperature of a fridge is too high or too low. We model this by associating a penalty with certain values of x. Specifically, as long as x ∈ {0 . . . Xmax}, this penalty is 0. When x < 0 or x > Xmax, the penalty is proportional to the gap between x and the desirable states. Explicitly, we define a penalty vector κ = (κi)i∈M indexed by the states

i = (s, x, y): κi =    0 if x ∈ {0, . . . , Xmax} -δ • x/Xmax if x < 0 δ • (x -Xmax)/Xmax if x > Xmax for some positive constant δ > 0.
When there is no DR, the appliance automatically switches on or off just before it reaches an undesirable state. In this specific scenario there is no penalty.

Mean-field Approximation

The role of each DR appliance is symmetrical. Hence, the state of the appliances is fully characterized at time t by an occupancy measure M N (t) = (M N i (t))i∈M, where M N i (t) is the fraction of appliances that are in state i at time t. We represent M N (t) as a row vector (of infinite dimension). Unless otherwise specified, all actors know M N (0) = M0.

We consider that the number of flexible appliances N is large. It is shown in [7, Section 3.1] that, as N grows, the M N (t) goes to M (t), which evolves as:

Mj(t + 1) = i∈M Mi(t)π(i, j, ai(t)). (2)
Moreover, as N goes to infinity, the optimal control for the mean field limit M is asymptotically optimal for the system with a finite number of appliances [START_REF] Gast | A mean field approach for optimization in discrete time[END_REF][START_REF] Gast | Mean field for markov decision processes: from discrete to continuous optimization[END_REF]. We denote by Π 0 and Π1 the two matrices:

Π 0 (i, j) = π(i, j, 0) and Π 1 (i, j) = π(i, j, 1) -π(i, j, 0).
Using Equation ( 1), the evolution of the occupancy measure is

Mj(t + 1) = i∈M Mi(t)Π 0 (i, j) + i∈M Mi(t)ai(t)Π 1 (i, j).
In what follows we assume that the occupancy measure M (t) can be fully observed at t. Thus, we replace a(t) by a more convenient control U (t), Ui(t) = Mi(t)ai(t). Formally, we define the set XU of processes (U (t)) such that M (0) = M0 and for all t:

M (t + 1) = M (t)Π 0 + U (t)Π 1 , (3)
Ui(t) ∈ [0; Mi(t)] ∀i ∈ M (4)
The two controls a and U are equivalent in the following sense: a process (M (t), M (0) = M0) describing the evolution of the occupancy measure can be obtained via (2) using a control process (a(t)) if and only if there exists U ∈ XU such that Ui(t) = Mi(t)ai(t) for all t, i. These equations remain valid if the same actions a(t) are sent to all appliances, or if each individual received an action. In that case, a(t) is the average over all actions. See the discussion in [7, Section 5.1].

Finally, the power consumption and penalty, defined in §3.2, can be expressed as scalar products 1 with M (t):

• The power consumption of flexible appliances is

F a (t) = f da (t) + F (t) = PonM (t) • ς , (5)
where the row vector ς = (ςi)i∈M is such that ςi=1 if i is an on state and ςi=0 otherwise. • The total penalty at time t is

M (t) • κ . (6)

Market Actors

We now describe the market actors and their objectives in two different scenarios. In the first scenario, the consumer also manages the flexible demand. Based on the current or anticipated market prices and on the blackout cost she decides whether to anticipate or delay a flexible load. The amount of flexible load close to undesirable states also influences this decision. In the second scenario, the DR operator is a stand-alone market actor that directly buys energy from the market.

Scenario 1: Consumer Controls DR Actions

The consumer controls the quantity of energy ED(t) purchased on the real-time market at time t. In this scenario, she also manages the the flexible demand F (t) by influencing the state of the flexible appliances. The payoff of the consumer at time t is a sum of three terms:

• We consider linear utilities of satisfied non-controlled consumption v min{D a (t), ED(t) + g da (t) -F (t)} and of satisfied controlled consumption v F a (t). • The disutility of a blackout is expressed as a sum of two terms: a linear term in the amount of unsatisfied non-controllable demand -c bo (D a (t) + F (t) -ED(t)g da (t)) + and a penalty of pushing controllable demand in undesirable states -M (t) • κ defined in [START_REF] Feinberg | Load forecasting[END_REF]. As the consumer aggregates many individuals, she may curtail only part of their demand. For example, this corresponds to interrupting the service for a neighborhood in case of insufficient resources, rather than interrupting an entire city. We treat controllable demand differently, since it is actively interrupted or encouraged to consume. • The amount of money spent on buying energy on the two markets at prices P (t) and p da (t) is captured as -P (t)ED(t) -p da (t)(g da (t) + f da (t)). To summarize, the total payoff of the consumer is

W tot D (t) := v min{D a (t), ED(t)+g da (t)-F (t)}+v F a (t) -c bo (D a (t)+F (t)-ED(t)-g da (t)) + -M (t) • κ -P (t)ED(t) -p da (t)(g da (t)+f da (t)).
We strip away all the terms that are not controllable in real-time. Taking into account ( 5) and ( 6), the real-time component of the payoff is:

WD(t) = M (t) • (v Ponς -κ) -P (t)ED(t) (7)
-(v+c bo )(ED(t)-D(t)-PonM (t) • ς +f da (t)+r da ) -.
where x -= max(-x, 0), for any real number x.

In the real-time market the consumer maximizes her expected payoff over the duration T of a day, given the scheduled day-ahead decisions. Her objective is arg max

E D ,U ∈X U (F t )-adapted WD, where WD := E T t=0 WD(t) .
The supplier sells a quantity ES(t) on the real-time market, and produces G(t) subject to the following constraints at each time t:

1. the sold quantity cannot exceed the produced quantity, i.e., ES(t) ≤ G(t) + Γ(t), 2. the real-time production has ramping limitations ζ -< 0 and

ζ + > 0: ζ -≤ G(t )-G(t) t -t ≤ ζ + , ∀t > t.
We denote the set of processes (ES(t), G(t)) that satisfy these constraints by XS. Thus, taking into account the marginal generation cost of the day-ahead production c da and that of the real-time production c, the payoff of the supplier is:

W tot S (t) = (p da (t) -c da )(f da (t) + g da (t)) + P (t)ES(t) -cG(t).
Removing the terms that are not controllable in real-time, this payoff becomes

WS(t) = P (t)ES(t) -cG(t). (8)
Finally, the supplier maximizes her expected payoff arg max

(E S ,G)∈X S (F t )-adapted WS, where WS := E T t=0 WS(t) .
Dynamic Competitive Equilibrium. We recall the notion of dynamic competitive equilibrium [START_REF] Cho | Efficiency and marginal cost pricing in dynamic competitive markets with friction[END_REF][START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF].

Definition 1. A dynamic competitive equilibrium (P e , E e D , E e S , G e , U e
) is a set of (Ft)-adapted price and control processes that satisfy:

(E e D , U e ) ∈ arg max E D ,U ∈X U WD, (9)
(E e S , G e ) ∈ arg max

(E S ,G)∈X S WS, (10)
E e D = E e S . (11)
In the above definition, [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF] means that (E e D , U e ) constitutes an optimal control from the selfish consumer's perspective. Similarly, [START_REF] Koch | Modeling and control of aggregated heterogeneous thermostatically controlled loads for ancillary services[END_REF] states that (E e S , G e ) is optimal from the supplier's perspective. Finally, [START_REF] Koutitas | A delay based optimization scheme for peak load reduction in the smart grid[END_REF] is the market constraint. Note that in (9), the consumer is not subject to the supplier's constraints and vice-versa for [START_REF] Koch | Modeling and control of aggregated heterogeneous thermostatically controlled loads for ancillary services[END_REF]. See [START_REF] Cho | Efficiency and marginal cost pricing in dynamic competitive markets with friction[END_REF] for a discussion.

Scenario 2: Standalone DR operator

In this scenario, we consider that end-users are charged at two different retail prices per energy unit: one for the non-controllable appliances, and another one for the controllable appliances. The bill for the controllable appliances is collected by the DR operator. The rationale of the DR operator is that controllable appliances (e.g., fridges, boilers) are always operational and that they consume roughly a constant (predictable) quantity of energy per day. Hence, from the DR operator's viewpoint, it is as if consumers pay a constant subscription fee per month for a provided service (e.g., to keep the fridge cool). The DR operator has the possibility of performing arbitrage on the real-time market thanks to the flexibility of the appliances.

The consumer controls only the quantity of purchased real-time energy ED(t) for serving the non-controllable demand. In this scenario, her payoff at each time t is:

W tot D (t) :=v min(D a (t), ED(t) + g da (t))+ -c bo (D a (t) -ED(t) -g da (t)) + -[P (t)ED(t) + p da (t)g da (t)].
The real-time component of the consumer's payoff is

WD(t)=-(v+c bo )(ED(t)-D(t)+r da ) --P (t)ED(t). ( 12
)
The DR operator is an independent actor that controls the real-time consumption of controllable demands F (t). On the real-time market, when F (t) < 0, she acts as a virtual energy supplier and gets payed for generating "nega-Watts". She generates revenue via arbitrage on the real-time market. Her payoff function is

W tot F (t)= -p da (t)f da (t)-P (t)F (t)+v F a (t)-M (t)•κ .
By Equation [START_REF] Cho | Efficiency and marginal cost pricing in dynamic competitive markets with friction[END_REF],

F (t) = PonM (t)κ T -f da (t)
. Hence, her real-time controllable payoff is

WF (t)=-P (t)PonM (t) • ς +M (t) • (v Ponς-κ) (13)
The DR operator maximizes her expected payoff subject to (Ft)-adapted feasible controls U ∈ XU , specifically:

arg max

(U )∈X U (F t )-adapted WF , where WF := E T t=0 WF (t) .
The supplier behaves the same way as in the previous scenario and has the same payoff [START_REF] Gast | Mean field for markov decision processes: from discrete to continuous optimization[END_REF].

A similar definition of a dynamic competitive equilibrium can be written in this scenario. A competitive equilibrium is now the tuple (P e , E e D , F e , E e S , G e , U e ). The difference with respect to the previous case is that the consumer who was adjusting the amount of bought energy via the real-time flexible demand F e now needs to pay for F e (t) = PonM e (t) • ς -f da (t) at the price of the real-time market as a form of virtual energy generation. The market constraint [START_REF] Koutitas | A delay based optimization scheme for peak load reduction in the smart grid[END_REF] becomes

E e S (t) = E e D (t) + F e (t) for all t. (14)

SOCIAL WELFARE THEOREM

In both presented scenarios we compute the social welfare as the sum of the actors' payoffs when the market clearing constraint is satisfied. In the case of a standalone DR operator (Scenario 2), the total payoff is the sum of Equations ( 8), [START_REF] Boudec | Stability of a Stochastic Model for Demand-response[END_REF] and [START_REF] Meyn | Ancillary service to the grid from deferrable loads: the case for intelligent pool pumps in Florida[END_REF]. It is equal to

Wtot(t) := -(v+c bo )(ED(t)-D(t)+r da ) --cG(t) (15)
+ M (t)•(v Ponς-κ) + P (t)(ES(t)-ED(t)-F (t)).
When the market constraints ( 14) are satisfied, the bought energy is equal to the sold energy at each time step. Hence, the price disappears from the total welfare. A socially optimal allocation maximizes the total expected payoff. As Wtot(t) is increasing in ES(t), we can further simplify [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF]. Taking into account the constraints on ES(t), we get that ES(t) = G(t) + Γ(t) maximizes Wtot(t). Using [START_REF] Mohsenian-Rad | Optimal residential load control with price prediction in real-time electricity pricing environments. Smart Grid[END_REF] to express ED and denoting the random (uncontrolled) component of the social welfare by Z(t) := Γ(t) -D(t) + f da (t) + r da , a socially optimal allocation maximizes the total expected payoff:

W tot (t) =-(v+c bo )(G(t) + Z(t)-PonM (t) • ς ) - + M (t) • (v Ponς-κ) -cG(t). (16)
We call (G * , U * ) a socially optimal allocation if it is a solution to the following problem:

arg max G∈X S , U ∈X U are (F t )-adapted E T t=0 Wtot(t) (17)
A direct computation shows that the total welfare is the same in Scenario 2. This implies that the socially optimal allocations are the same in both scenarios.

The following result holds. It is proved in Appendix A.

Theorem 1. If (P e , E e D , E e S , G e , U e ) is a dynamic competitive equilibrium in Scenario 1, then the allocation (G e , U e ) is socially optimal.

Similarly, if (P e , E e D , F e , E e S , G e , U e ) is a dynamic competitive equilibrium in Scenario 2, then the allocation (G e , U e ) is socially optimal.

NUMERICAL SOLUTION OF THE PROBLEM

In the rest of this paper we compute numerically the equilibrium prices and welfare given by Theorem 1. To this end we use a specific model for the random component Z(t) that is both tractable and realistic.

Branching Trajectory Forecasts

We use discrete time and evaluate our model over a period of 24h. Each time slot represents a five minute interval [START_REF]Pjm operational data[END_REF]. Hence, the real-time market lasts for T = 288 time slots.

Inspired by the methodology of Pinson et al. [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF], we assume that the forecast error process Z can be represented by a finite number of 2 τ trajectories. For a trajectory ω ∈ {1 . . . 2 τ }, the value of the forecast error at time t is denoted Z ω (t). To model the fact that the forecast for the near future (e.g., the next hour) is more accurate than the forecast for a distant future (e.g., 12 hours from the present), we consider that the possible trajectories for the forecast error coincide initially before separating in distinct branches at certain moments {t1 . . . tτ }. Thus, over time the number of branches increases exponentially: between times t k-1 and time t k , there are 2 k possible such branches of the forecast errors. We consider that initially, between times t0 = 0 and t1, there are two equally likely branches: the odd trajectories and the even trajectories.

At each time t k , each group of trajectories separates into two equiprobable groups of trajectories. Formally, between t k-1 and t k , all trajectories ω that have the same remainder modulo 2 k share the same values:

∀ω < 2 τ -2 k , ∀t ∈ [t k-1 ; t k ] : Z ω (t) = Z ω+2 k (t). ( 18
)
An example of trajectory is represented in Figure 2. For more clarity, we only represent eight trajectories. In the rest of the paper, the numerical evaluations use 512 trajectories. On this figure, we see that before t1 := 4.5h, there are two possible forecast errors (the "odd" and "even" trajectories) Between t1 and t2 = 12h, there are four branches. After t2, all trajectories are distinct.

These trajectories are generated using the covariance of real forecast errors in the UK. The methodology and the algorithm used will be described in §6.1.

Observability assumptions

We assume that at each of the branching instants t k , players can observe the forecast error and can know on which branch the system is evolving at t k . However, they cannot predict which branch the forecast errors will follow after t k . For example, a time t1 = 4.5h, players can observe if the forecast error is one of the four odd (Z 1 , Z To ease the presentation, this example has only eight trajectories. The numerical results of Section 6 use 512 trajectories. These trajectories are generated using data from the UK (see §6.1).

0 t_1 8h t_2 16h 24h -10 -5 0 5 10 G 1 =...=G 8 G 1 =G 3 =G 5 =G 7 G 2 =G 4 =G 6 =G 8 G 1 =G 5 G 3 =G 7 G 2 =G 6 G 4 =G 8
time (in hours) Generation G (in GW) of the four even trajectories. The only information about the future that they have is that each trajectory is equiprobable.

The set of control variables (G, U ,. . . ) is restricted to the set of admissible controls, which are the sequence of causal decisions: a decision taken at time t can use the information available at time t as well as the statistics about future trajectories but not the exact future realization of the forecast errors. Using our observation model, this means an admissible control for G and U is a sequence of values G ω (t) and U ω (t) that satisfies the generation and DR constraints [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] such that for all ω < 2 τ -2 k and all t ∈ [t k ; t k+1 ]:

G ω (t) = G ω+2 k (t); U ω (t) = U ω+2 k (t). ( 19
)
For example, if the forecast error is represented by the trajectories of Figure 2, an example of admissible control for G is depicted in Figure 3. This control is causal. At time t0 = 0 it knows that the forecast error may follow one of the two branches. The chosen control for the interval [0; t1] is a compromise that yields the best average performance. At the observation point t1, players observe whether the branch that was followed was that of an odd or of an even trajectory. Since each branch separates again in two branches, there are two causal sequences of controls for the interval [t1; t2] that can be applied, depending on the observations made at time t1. Similarly, at time t2 the four possible observations lead to four possible decisions.

Casting the Stochastic Problem into Coupled Deterministic Problems using ADMM

To compute the equilibrium, we use an optimization algorithm called the alternating direction method of multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. ADMM is an iterative process to solve deterministic constrained maximization problems with concave 2 objective and linear constraints:

max x∈R d ,z∈R e f (x) + g(z) subject to Ax + Bz = c, (20)
where f and g are two concave functions that take real values or -∞; A and B are two matrices of size d × dx and d × dz; and c is a vector of size d.

Let ρ > 0. For each (x, z, λ) ∈ R dx+dz +d , the augmented Lagrangian Lρ(x, z, λ) is defined as:

Lρ(x, z, λ) = f (x) + g(z) - ρ 2 Ax+Bz+λ-c 2 . ( 21
)
ADMM performs the following iterations at step k + 1:

x k+1 := arg max x Lρ(x, z k , λ k ) ( 22
)
z k+1 := arg max x Lρ(x k+1 , z, λ k ) ( 23
)
λ k+1 := λ k + Ax k+1 + Bz k+1 -c, (24)
where x 1 , z 1 and λ 1 have been initialized to any vector of the proper size.

It is shown in [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] that under very general assumptions, the ADMM iterations converge to a solution of the optimization problem [START_REF] Wilson | Architecture of power markets[END_REF], regardless of the initial conditions. ADMM also converges when f or g are not strictly concave and when they take infinite values.

The social maximization problem [START_REF] Spees | Demand response and electricity market efficiency[END_REF], is a stochastic optimization problem where the expectation is taken over all possible forecast processes. As we assume that the forecast errors can be represented by a finite number of trajectories, this problem can be written arg max

G,U ∈admissible ω Wtot(G ω , U ω ).proba(ω), ( 25
)
where the set admissible controls G, M denotes controls that satisfy generation constraints XS, DR constraints XU , as well as causality constraints [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF].

In particular, Equation [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF] implies that the decisions taken at a time t depend only on the information up to time t. As such, they are a compromise between the possible optimal decisions taken by a hypothetical omniscient controller that can observe the future.

Equilibrium Computation using ADMM

The formulation of ADMM, given by Equation ( 20), requires the optimization space to be partitioned into two groups of variables with no objective functions that depends on both groups, like h(x, z). To perform this separation, we enlarge our optimization space. In addition to the original variables G ω (t), M ω (t) and U ω (t), we add, for each t < T and ω, the variables M

ω (t), G ω (t), G ω (t), E ω D (t), E ω F (t), E ω S (t), U ω (t), U ω (t)
and the linear constraints given by Equation [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF] and

M (t + 1) = M (t)Π 0 + U (t)Π 1 ; M (t)=M (t); ( 26
)
2 ADMM is usually presented as a solution to solve convex minimization problems. In this paragraph, we keep the concave maximization formulation of Section 3.

U (t)=U (t); U (t)=M (t)-U (t); EF (t)=PonM (t)•ς (27) G(t) = G(t) = G(t); G(t) = ES(t). ( 28
)
The first group of variables (the "X variables") contains the variables M ω (t), U ω (t) and U ω (t), G ω (t). The corresponding objective function f (M, U , U , G) is:

t,ω M ω (t) • (v Ponς-κ) -cG ω (t)+Pos(U ω , U ω ), (29)
where Pos(x, x ) = 0 if all coordinates of the vectors x and x are non-negative and -∞ otherwise.

The second group of variables contains the variables M , U , G, G, ED, EF , ES. The corresponding objective function g(U, M , G, G, ED, EF , ES) is equal to:

t,ω -(v+c bo )ED(t) -+ Pos(ζ + -G ω (t+1) + G ω (t)) + Pos(G ω (t+1) -G ω (t) -ζ -) + Equal0(E ω S (t) + Z ω (t) -E ω D (t) -E ω F (t)) ( 30
)
where Equal0(x) = 0 if all coordinates of x equal 0 and -∞ otherwise.

Analysis of the ADMM Algorithm

When the constraints given by Equations (26-27) and the two terms Pos(•) of Equation ( 29) are satisfied, the evolution of M (t), given by Equation ( 3) is respected. Similarly, the constraints given by ( 28) and the two term Pos(•) of Equation ( 30) ensure that the ramping constraints of G are satisfied.

This implies that our algorithm converges to the set of socially optimal allocations. Moreover, it also computes competitive equilibria. The following theorem summarizes these results. It is proved in Appendix B.

Theorem 2. Let E k D , E k S , G k , U k
be the values after the kth iteration and let (λ k G ) ω (t) be the multipliers associated with the constraints E ω S (t) = G ω (t). Then: (i) There exists P such that lim k→∞ ρ(λ k G ) ω (t) = P ω (t).

(ii) Any subsequence of E k D , E k S , G k , U k has a subsequence that converges. Let E ∞ D , E ∞ S , G ∞ , U ∞ be its limit. Then • (P, E ∞ D + E ∞ F , E ∞ S , G ∞ , U ∞ ) is a competitive equi- librium for Scenario 1. • (P, E ∞ D , E ∞ F , E ∞ S , G ∞ , U ∞ ) is a competitive equilib- rium for Scenario 2.
As a consequence, G ∞ , U ∞ is socially optimal. (iii) If the maximization algorithm has a unique solution, the algorithm converges to this solution. If there are multiple optimal solutions, they form a convex set and the algorithm converges to this set (for the Euclidean distance).

Remark. Even if we consider piecewise-linear objective functions, ADMM only requires the objective to be concave. Hence, the same algorithm can be applied to more general cost functions and Theorem 2 also holds.

NUMERICAL EVALUATIONS

Parameters and Trajectory Generation

As in [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF], we use ζ + =1GW/h, ζ -=3ζ + . We set the ratio between the blackout cost and the generation cost to (v + c bo )/c=10. The values of v, c bo and c are normalized so the flexibility provided by fridges L parametrized below leads to a welfare gain of 100 when Pon is large.

We compare three cases of appliances [START_REF] Koutitas | A delay based optimization scheme for peak load reduction in the smart grid[END_REF]:

• Fridges L (large inertia fridges) -the on period is set to Xmax/γon = 60min and the off period is 120min, with a mini-cycle prevention period Ymax/θon=0 or 20min. It corresponds to Xmax=6, γon=θon=θ off =0.5, γ off =0.25. • Fridges S (small inertia fridges) -same as Fridges L but with Xmax = 3. The only difference is that the on period lasts 30min and the off period lasts 60min. • Boilers -the on period is 4h and the off period is 20h.

It corresponds to Xmax = 12, γon = 0.25, γ off = 0.05 and θon = θ off = 0.5. The random transitions of our model of appliance account for various random evolutions of the internal temperature, e.g., due to fridge doors opening. The second fridge model corresponds to a fridge with a smaller inertia. It consumes the same average power but has less flexibility, typically because the temperature deadband is narrower. The ratio between the duration of off and on periods is 2 for fridges and 5 for boilers. Our numerical evaluation shows that, qualitatively, the behaviors of the three models are the same.

Trajectory Generation using UK Data

To construct wind forecast errors, we use wind production and day-ahead wind production forecast in the time interval from June 2009 to April 2012, obtained from the BMRA data archive (elexonportal.co.uk). We normalize the values of production and forecast to maintain a constant wind capacity3 of 26GW all throughout the 1300 days. We obtain 1300 samples of forecast error trajectories {ε d (t)}. ε d (t) is the error for the day d ∈ {0 . . . 1300} at time t ∈ [0; 24h].

A generic method to generate forecast error trajectories has been introduced in [START_REF] Pinson | From probabilistic forecasts to statistical scenarios of short-term wind power production[END_REF]. To generate one trajectory, the authors compute the covariance Σ of the process ε d (a 288 × 288 matrix in our case) and then generate a multivariate normal vector of covariance Σ. We adapt their method to our case of branching forecast. Algorithm 1 generates a set of trajectories such that each one has covariance Σ, and such that two trajectories that share a branching point depend on each other only via their values before the branching point.

Input: Covariance matrix Σ, sequence t1, . . . , tτ-1 Output: Trajectories Z 1 , . . . , Z 2 τ 1 A ← lower Cholesky decomposition of Σ; t0 ← 0; 2 N1 ← vector of T i.i.d. normal random variables;

3 for k ← 0 to τ -1 do 4 for ω ← 1 to 2 k do 5 N ω+2 k ← [N ω (1 . . . t k ); N ω+2 k ], where N ω+2 k is a vector of T -t k i.i.d. normally distributed random variables; 6 end 7 end 8 for ω ← 1 to 2 τ do 9 Z ω ← A × N ω ; 10 end
Algorithm 1: Generation of 2 τ branching trajectories satisfying Theorem 3. The notation N ω (1 . . . t k ) denotes the first t k elements of the vector N ω . All vectors are column vectors. [A; B] is the concatenation of two column vectors to form one column vector.

These properties are summarized in Theorem 3.

Theorem 3. Assume that Σ has full rank. Then, Algorithm 1 generates a set of trajectories such that:

(i) individually, each trajectory Z ω is a multivariate normal vector of covariance Σ; (ii) If a subset of trajectories Z i 1 . . . Z i k is fixed, then the distribution of another trajectory Zi is the distribution of a multivariate normal vector of covariance Σ subject to respecting Equation (18).

This theorem is proved in Appendix C. The proof can be adapted if Σ is not full rank but it is more technical.

As a consequence, all trajectories Z 1 . . . Z 2 τ are equally likely. Moreover, as τ grows, the process obtained by following the branching trajectories "approaches" a normally distributed random process with covariance Σ.

Impact of DR Capactiy and of the Nature of Elastic Loads

For a given power flexibility Pon, we denote W * (Pon) the optimal value of social welfare in Equation [START_REF] Spees | Demand response and electricity market efficiency[END_REF]. W * (0) corresponds to the social welfare without DR. For Pon > 0, the optimal welfare W * (Pon) depends on the type of appliance, while W * (0) does not.

In Figure 4 we plot the "relative" social welfare, i.e., the difference W * (Pon) -W * (0) as a function of Pon. We compare the three types of appliances. For fridges L, we plot both the case without mini-cycle avoidance (Ymax = 0) and with mini-cycle avoidance (Ymax > 0). As expected, the welfare is then increasing and concave in Pon. We observe that, as Pon grows, the social welfare increases and tends to stabilize at the value 100, regardless of the appliance type or the value of Ymax, although for a fixed capacity of DR, increasing the mini-cycle avoidance parameter Ymax diminishes its benefits. Moreover, for a given appliance type, the gain of using demand-response is almost linear in Pon for small values of the installed capacity Pon. In view of the similarity of the appliance responses, for the rest of the numerical evaluation we only show the results for the fridges L.

Non-Observability of DR

The previous figure assumes that the day-ahead market is able to predict without error the consumption of demandresponse appliances. In practice, this can be done because each appliance evolves mainly independently of the others. As a results, the relative error of current load-prediction techniques is less than 1% [START_REF] Feinberg | Load forecasting[END_REF]. However, one of the concerns of demand-response is that the presence of a large player that controls many appliances will introduce synchronization between devices, which will make the behavior of the system less predictable. For example, a small congestion can cause lots of DR loads to be delayed, which can cause a bigger problem in the future because these delayed loads were not forecast at this date [START_REF] Boudec | Stability of a Stochastic Model for Demand-response[END_REF].

In this section, we show that when the day-ahead market cannot fully observe the demand-response state, adding too much demand-response can be harmful, even if the real-time market has full information about the states of appliances. More precisely, we assume that:

• the day-ahead market cannot observe M (0), the state of the demand-response at time 0. It assumes that f da (t) is equal to the overall average consumption of all demand-response appliances and plans the day-ahead generation g da (t) accordingly. • the that controls flexible loads can observe the states of all DR appliances in real time. We plot the relative social welfare W * (Pon) -W * (0) for this scenario in Figure 5. When the flexible power is small (Pon < 10GW), the benefit of demand-response is close to the one with full information. However, when the power flexibility is large, adding more demand response decreases the total welfare. The non-observability problem is probably hidden today as the power flexibility of DR is small but this is a threat for the future. The results for fridges S and boilers are very similar and are not shown here.

Comparison with Energy Storage Systems

Demand-response is often viewed as a virtual energy storage system, which is potentially cheaper than real storage systems such as batteries. Instead of charging and discharging a battery, DR allows to anticipate or delay the consumption of appliances. In this section, we compare energy storage and flexible loads. We highlight two differences: for low capactiy, storage provides more flexibility because it reacts faster. However, at high capactiy, storage is less efficient because of the energy losses at each charging/discharging cycle.

We consider a storage model like in [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF]: the storage has an energy capacity of Bmax and a maximum charging and discharging power capacity Cmax and Dmax. The cycle efficiency is η: only a fraction η of the stored energy can be retrieved. To perform a fair comparison of DR and storage, we consider a storage system that has the same flexible power: we set Cmax = Dmax = Pon/2. We chose Bmax to reflect the quantity of the energy stored in DR appliances. Recall that when an appliance is on, its internal state x increases with probability γon. When off, its internal state decreases with probability γ off . Thus, the average consumption of the DR appliances is Ponγ off /(γon + γ off ). The internal state of an appliance that avoids the undesirable states oscillates between 0 and Xmax. We set Bmax equal to the difference in consumption between the case where all DR appliances have an internal state of Xmax and the case where they are all in state 0. Hence, we take:

Bmax = PonXmax γon + γ off × 5min (31)
With the fridges L model, this implies Bmax= Pon 2 GWh. We adapt our ADMM algorithm of §5.3 for the case where a storage system replaces DR. The results are reported in Figure 6. We plot the welfare as a function of the flexible power Pon in four scenarios: presence of DR but no storage (two dashed curves); or presence of storage but no DR (two solid curves).

We observe that, at low power capacities, storage provide higher gain than DR. Two reasons explain this:

• the fatigue effect: because of internal state constraints, some appliances cannot be switched anymore. • the mini-cycle effect Ymax > 0: appliances that have just been switched cannot be switched instantaneously. For large capacities, the situation is reversed. The welfare of DR and that of idealized storage saturate at the same value. This value is strictly larger than the one of 70%-efficient storage. For large power capacity, DR outperforms realistic storage (typical efficiency 70 -90%).

Price Equilibrium

When there is no demand-response or storage, it has been shown in [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF] that the price process oscillates between 0 and the choke price v + c bo . In this case it is never equal to the marginal production cost. The situation is different in the presence of storage. In [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF], we exhibit two different situations, depending on the charge-discharge efficiency of storage η: a) When η = 100% (perfect storage), the prices do in fact concentrate on the marginal production cost c as the quantity of storage gets large. b) When η < 100% (realistic storage), the prices do not concentrate on c, but exhibit two modes around c √ η and c/ √ η. These results were obtained assuming that the forecast error is stationary and can be represented by a Brownian motion with stationary increments. We simulate the same scenarios using our more realistic non-stationary branching forecast model and our ADMM implementation. The price distributions are reported in Figure 7. They confirm that the results of [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF] are robust to the forecast model: for η = 1, the prices concentrate on c = 1. For η < 1, the price distribution exhibit two modes around √ η ≈ 0.89 and 1/ √ η ≈ 1.12. This price spread is explained by the fact that only a fraction η of the stored energy can be retrieved. As a consequence, a storage owner will store energy only if it can be sold at a 1/η times larger price. In the case of DR, the energy is virtually stored or retrieved by delaying or anticipating consumption. The corresponding prices are reported in Figure 8. It confirms that the system behaves like a 100% efficient storage system (like in §6.4): even when Ymax > 0, the prices concentrate on the marginal production cost c.

Incentive to Install Demand Response

The storage/demand-response operator makes money via price arbitrage. An independent storage operator buys cheap energy and resells it when prices rise. Similarly, a demandresponse operator tries to shift the appliances' consumption to when the energy is inexpensive. The flexibility brought by demand-response or by storage results in more concentrated prices. In [START_REF] Gast | Impact of storage on the efficiency and prices in real-time electricity markets[END_REF], we have shown that, to maximize their revenue, independent storage operators will underdimension their storage system. In this section, we show that this also holds for DR operators. We simulate the same scenarios as in Figure 6 and we compute the expected revenue of the DR operator or of the storage operator as a function of the flexible power capacity. In all cases, this gain reaches a maximum and then decreases. This maximum is attained at Pon≈1GW for batteries and fridges L with Ymax=0. When Ymax=20min, this maximum is reached at Pon≈3GW. More problematically, the maximum is attained at a value that is far from optimal for the social welfare: This entails that even if we disregard installation costs, the operator seeking to maximize her revenue will not deploy more than 1GW (or 3GW) of DR capacity. In Figure 4, we observe that these values of Pon lead to a gain of only 40-60% of the potential welfare benefit of DR/storage. Hence, to maximize her revenue, a DR operator will deploy a suboptimal capacity from a social perspective.

CONCLUSION

In this paper we studied the effects of demand-response in a real-time electricity market. We modeled a population of flexible appliances that can anticipate or delay their consumption. This model accounts for fatigue via undesirable states and enforces mini-cycle avoidance via an additional dimension in the state space. We showed that when demand-response is traded by price-taking actors, the dynamic competitive equilibrium is socially efficient. We considered a realistic and tractable trajectorial forecast error model where trajectories separate into branches over time. We parametrized the model using real wind forecasts from the UK and solved numerically the model by using ADMM. We showed that if the appliance states are unknown, too much demand-response can be detrimental. We compared demand-response to an energy storage system and concluded that a sufficiently large amount of demand-response is as beneficial as a perfect storage system. Demand-response relaxes the ramping constraints of generation and smooths prices, that concentrate around the marginal generation cost.

APPENDIX A. PROOF OF THEOREM 1

We give the proof for Scenario 1. The proof for Scenario 2 is similar, with three actors instead of two. It is similar to the proof of [9, Theorem 2] and [START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF].

Let W(G, U ) be the social welfare. For a given price P , we denote WD(ED, U, P ) the welfare of the consumer when she takes the decisions ED, U and by WS(ES, G, P ), the welfare of the supplier for the decisions ES, G. We consider the social optimization problem [START_REF] Spees | Demand response and electricity market efficiency[END_REF] The last equality holds because a competitive equilibrium implies E e D = E e S . The one before last equality holds because each player maximizes his social welfare.

Hence, the allocation (G e , U e ) is socially optimal.

B. PROOF OF THEOREM 2 B.1 Convergence of ADMM

The problem formulated in §5.3 is a maximization problem with a concave objective and linear constraints. Moreover, there exists a allocation such that all constraints are satisfied (for example G ω (t) = G(0), and M ω (t+1) = M ω (t)Π 1 ). Hence, the Slater's conditions are satisfied and the unaugmented Lagrangian has a saddle point [START_REF] Boyd | Convex optimization[END_REF]Chapter 5].

We use the generic notation of §5.2 and denote the variables after iteration k by x k , z k , λ k . By [3, §3.2.1], the existence of the saddle point implies:

(I) x k , z k approach feasibility: lim k→∞ Ax k + Bz k = c. (II) λ k converges to a dual optimal point λ ∞ . (III) f (x k ) + g(z k ) converges to the optimal value of the social optimal problem: lim k→∞ f (x k ) + g(z k ) = sup The augmented Lagrangian Lρ(x, z, λ) goes to infinity as x or z goes to infinity. Hence, the iterations x k , z k are bounded and any subsequence x k , z k has converging subsequence. Let x ∞ , z ∞ be its limit. By (I), (x ∞ , z ∞ ) is a feasible allocation. Moreover, by (III), it is socially optimal allocation. This shows (i) and (iii). Since λ ∞ is a dual optimal point, (x ∞ , z ∞ , λ ∞ ) is a saddle point. In particular, (x ∞ , z ∞ ) is a maximizer of Lρ(x, z, λ ∞ ).

B.2 Definition of the Price as a Multiplier

In our ADMM formulation of §5.3, the augmented Lagrangian Lρ(x, z, λ), is equal to: Let (λ ∞ G ) ω (t) be the limit of (λ k G ) ω (t) as k grows. We define the price process P ω (t) as P ω (t) = ρ(λ ∞ G ) ω (t).

B.3 This Price Leads to an Equilibrium

We now show that this price leads to a competitive equilibrium for Scenario 2. The proof is similar for Scenario 1. We denote the limiting values of λ k , x k , z k and by λ ∞ ,

x ∞ = (M ∞ , U ∞ , U ∞ , G ∞ ), z ∞ =(M ∞ , U ∞ , G ∞ , G ∞ , E ∞ D , E ∞ F , E ∞ S
). Supplier's problem -We first show that G ∞ is an optimal schedule for the supplier's -The vector (x ∞ , z ∞ ) maximizes Lρ(x, z, λ ∞ ). Hence, (G ∞ , G ∞ , G ∞ ) maximizes The quantity ρ 2 E ω S -G ω +λ ω G 2 rewrites as

P ω •(E ω S -G ω ) + ρ 2 ( λ ω G 2 + E ω S -G ω 2
). Since G ∞ maximizes (33) and satisfies G ∞ = G ∞ = G ∞ , it also maximizes fρ(G), where

fρ(G) = t,ω (P ω (t) -c)G ω (t) - ρ 2 E ω S (t)-G ω (t) 2 .
Let GS be an optimal schedule for the supplier and let us show that G ∞ also optimal for the supplier.

Figure 1 :

 1 Figure1: Markov model of an appliance. Each state is characterized by a tuple (s, x, y), with s ∈ {on, off}. The x-axis represents the internal state (e.g., temperature). The y-axis is the time before the device can be switched on or off.

Figure 2 :

 2 Figure 2: Possible trajectories of the forecast errors.To ease the presentation, this example has only eight trajectories. The numerical results of Section 6 use 512 trajectories. These trajectories are generated using data from the UK (see §6.1).

Figure 3 :

 3 Figure 3: Example of an admissible control for G corresponding to the forecast errors of Figure 2.

Figure 4 :

 4 Figure 4: Relative social welfare W * (Pon) -W * (0) as a function of Pon -the available power capacity of DR -for the three types of appliances.

Figure 5 :

 5 Figure 5: Relative social welfare W * (Pon) -W * (0) as a function of available power flexibility Pon when the day-ahead market cannot observe DR's initial state, compared to the case when the day-ahead can observe this initial state. The x-axis is between 0 and 100GW.

Figure 6 :

 6 Figure 6: Social welfare as a function of available capacity of demand-response Pon or of storage. The storage flexibility is similar to the one of DR: Cmax = Dmax=Pon/2 and Bmax=(Pon/2)GWh.

Figure 7 :

 7 Figure 7: Distribution of prices for a system equipped with storage of capacity Cmax = 5GW, Bmax = 5GWh and charge-discharge efficiency η = 1 or 0.8. For readability, we scale the values c, v and c bo to c = 1.

Figure 8 :

 8 Figure 8: Distribution of prices for a system equipped with DR appliances Fridges L with Ymax = 20min. For readability, we scale the values c, v and c bo to c = 1.

Figure 9 :

 9 Figure 9: Revenue for the DR or storage operator as a function of the installed capacity. The studied scenarios are the same as in Figure 6.

  and we relax the constraint ES = ED. Denoting the corresponding Lagrange multiplier by P , the Lagrangian L(ED, ES, U, G, P ) is equal to:W(G, U ) + E t t=0 P (t)(ED(t) -ES(t)) = WD(ED, U, P ) + WS(ES, G, P ).When the constraint ED = ES is respected, the social welfare W(G, U ) equals WD(ED, U, P ) + WS(ES, G, P ).If (E e D , E e S , G e , U e , P e ) is a competitive equilibrium, then supG,U W(G, U ) = sup G,U,E S ,E D s.t. E S = E D L(ED, ES, G, U, P e ) ≤ sup G,U,E S ,E D L(ED, ES, G, U, P e ) = sup E D ,U WD(ED, U, P e ) + sup E S ,G WS(ES, G, P e ) = W(E e D , U e , P e ) + WS(E e S , G e , P e ) = W(G e , U e ).

x

  ,z s.t.Ax+Bz=c f (x) + g(z)

S and E k+1 F after the kth iteration maximize - ρ 2 (G 2 +F 2 )

 k+1222 ES-G k+1 +λ k EF -F k+1 + λ k -(v+c bo )ED(t) -, (32) subject to ES + Z -EF -ED = 0.A necessary condition for that is thatE k+1 S -G k+1 +λ k G = -(E k+1 F -F k+1 + λ k F ). By Equation (24), the multipliers λ k+1 G and λ k+1 F are λ k+1 G := λ k G + E k+1 S -G k+1 and λ k+1 F := λ k F + E k+1 F -F k+1 .Hence, this implies that after this iteration, the two multipliers λ k+1 G and λ k+1 G are opposite (for all t, ω).

  For θ ∈ [0; 1], we define an allocationG θ := θGS + (1θ)G ∞ and g0(θ) = f0(G θ ) and gρ(θ) = fρ(G θ ). As G ∞ = ES we have g0(θ) = t,ω (P ω (t) -c)G ω θ (t);(34)gρ(θ) = t,ω (P ω (t) -c)G ω θ (t) -θ 2 ρ 2 GS-G ∞ 2 .(35)

  3 , Z 5 , Z 7 ) or one

	forecast error (in GW)	-10 -5 0 5 10	Z 2 =Z 4 =Z 6 =Z 8 Z 1 3 =Z =Z 5 7 =Z	Z 4 =Z 8 Z 2 =Z 6 Z 1 =Z 5 Z 3 =Z 7	Z 6 Z 1	Z 3	Z 2 Z 7
		0	t_1	8h	t_2 16h		24h
						time (in hours)

  -(v+c bo )E ω D (t) --cG ω (t) + Pos(U ω )+Pos(U ω ) + Pos(ζ + -G ω (t+1)+G ω (t))+Pos(G ω (t+1)-G ω (t)-ζ -) Bz -c + λ 2 ,where the linear equation Ax + Bz = c represents all linear constraints (26,27,28) and[START_REF] Wang | Dynamic competitive equilibria in electricity markets[END_REF]. Let (λ k G ) ω (t) and (λ k F ) ω (t) be the multipliers after the kth iteration of the ADMM, corresponding to the constraints E ω S (t) = G ω (t) and EF (t) = F (t) respectively. The values of E k+1

	t,ω
	+ Equal0(E ω S (t) + Z ω (t) -E ω F (t) -E ω D (t))
	-	ρ 2	Ax +

For two row vectors x and y, y denotes the transpose of y and x • y is the scalar product between x and y.

This scenario is envisioned for the UK in 2020, where 20% of the total electricity consumption is covered by wind.

* With financial support by Swiss Nano-Tera project "Smart Grid" and EU project QUANTICOL, 600708.

As G θ is a convex combination of two feasible allocation. Hence, it also satisfies the generation constraints. As GS is optimal for the supplier, the function g0(θ), attains its maximum in θ = 1. Moreover, this function is concave in θ and therefore, it has a right derivative in θ = 0 that satisfies g 0 (0) ≥ 0. Similarly, the function gρ is concave and attains its maximum in θ = 0. Hence, it has a right derivative in θ = 0 that satisfies g ρ (0) ≤ 0.

It should be clear from the expressions of g0 and gρ, given by Equations (34) and (35), that g 0 (0) = g ρ (0). This implies that g 0 (0) = 0. Hence, θ = 0 is also an optimum for g0, which implies that G ∞ is optimal for the supplier. DR operator's problem -We showed in §B.2 that the multipliers for EF = F and for ES = G are opposite. Hence, the optimality of E ∞ F , U ∞ for the DR operator under the price P follows from the same approach as for the supplier (by comparing the situations ρ = 0 and ρ > 0).

Demand's problem -Replacing ES by ED + EF -Z, maximizing Equation (32) on ED subject to ES + Z -ED -EF = 0 is equivalent to maximizing without constraint:

As the second term does not depend on ED, it maximizes

Hence, the optimality of ED for the Demand under the price P follows from the same approach as for the supplier (by comparing ρ = 0 and ρ > 0).

C. PROOF OF THEOREM 3

The matrix Σ is the covariance of the forecast error process ε. As such, Σ is a positive-definite symmetric matrix of size T × T and has real entries. Thus, there exists a real matrix A of size T × T such that Σ = A × A and A is a lower-triangular matrix, i.e., for all t ∈ {1, . . . , T } and s ∈ {t + 1, . . . , T }, At,s = 0. A is called the lower Cholesky decomposition of Σ. It also exists when Σ is not full rank.

Let Z ω be a trajectory. By construction, the vector N ω is a vector of T i.i.d. normally distributed random variables. Moreover, we have Z ω (t) = s At,sN ω (s). Hence, the covariance E [Z ω (t)Z ω (t )] is equal to AA = Σ, which implies that property (i) holds.

Let us fix a sub-set of trajectories Z ω 1 . . . Z ω k and let Z ω be another trajectory. Let τ be the largest branching time of Z ω shared with one of these trajectories. By construction, Z ω = AN ω . The matrix A and the vector N ω can be decomposed as:

where the matrices B, C and D have size τ × τ , (T -τ ) × τ and (T -τ ) × (T -τ ) and N1 has size τ . By construction, the distribution of Z ω given trajectories Z ω 1 , . . . , Z ω k is the distribution of Z ω given N ω (1), . . . , N ω (τ ). Let z = B.[N ω (1), . . . , N ω (τ )]. Because A has full rank, B has also full rank and is invertible. Hence, the distribution of Z ω conditioned on Z ω (t) = zt for t ≤ τ is the same as the distribution of Z ω given N ω (1) . . . N ω (τ ). The former is the conditional distribution of a multivariate Gaussian random variable of covariance Σ subject to Equation [START_REF] Su | Quantifying the effect of demand response on electricity markets[END_REF].


Christine Fricker
email: christine.fricker@inria.fr.

Nicolas Gast
email: nicolas.gast@epfl.ch.

Hanene Mohamed
email: hanene.mohamed@gmail.com

Hanene Mohamed Mean

Mean field analysis for inhomogeneous bike sharing systems

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Motivation. Bike sharing systems have generated research interest recently as such programs have been launched in many cities around the world. The concept is quite simple: the user takes a bike to a station, uses it and returns it at a station of his choice.

Despite their success, it is commonly observed that they do not behave very well. Users have to address the lack of resources, finding an empty station when in need of a bike, or a full station when returning a bike. The intuition is that the system is highly inhomogeneous, with very different arrival rates and destinations, between housing and working areas for instance, or up-hill and down-hill stations. These parameters may also vary with the time of day (see Borgnat et al. [START_REF] Borgnat | Shared bicycles in a city: A signal processing and data analysis perspective[END_REF]), but this issue is not considered here.

In the literature on this topic, few analytical models have been proposed, especially stochastic models, mainly because this leads to complex model that are are difficult to analyze. Nevertheless, such mathematical models could give insight on the behavior of the system or study different incentives algorithms for taking or returning bikes, designed to improve overall performance.

Goal. The paper deals with a stochastic model taking into account the finite number of bike posts at the stations. To our knowledge, it is the first attempt to analyze the problem of returning the bikes. Our aim is to study the impact of system inhomogeneity. We use an asymptotic approach, studying a system in which the number of stations and the fleet get large together. The proportion of bikes per station and the number of spots per station, which are two key parameters of the system, remain fixed. This asymptotic study is quite adequate for bike sharing systems, that have hundreds or thousands of stations.

Main result. The main result is the limiting stationary distribution of the number of bikes in the stations as the system gets large. The proof is based on mean field analysis. It amounts to considering the empirical measure process. In the special case of homogeneous model, it is simply the process Y (t) = (Y 0 (t) . . . Y K (t)) where Y i (t) is the proportion of stations with i bikes at time t. It converges when the system gets large to a deterministic dynamical system, solution of an ODE. Using a Lyapunov function for the dynamical system, we show that it has a unique attractor. This yields the convergence of the stationary distribution to this fixed point. Notice that it is often difficult to exhibit such a function.

The challenge is to deal with inhomogeneity of the system that makes the analysis of the system of differential equations difficult. While other works use numerical simulations of the system [START_REF] Gast | A mean field model of work stealing in large-scale systems[END_REF], we obtain analytically the limiting stationary distribution, by characterizing the equilibrium point and providing a Lyapunov function. Despite the cluster structure, our method uses a compact way to write the differential equation.

Mean field limit is an indirect way to obtain the limiting stationary distribution: from a system of fixed size, we take the limit as the size gets large, then the limit as time tends to infinity, identify the limits, and of course prove convergences. In our model, say I, the stationary distribution for a fixed size has a closed form expression. A more direct approach would be to obtain an asymptotic from this expression by an analytical method. The benefit of mean field analysis is that it works for models where the stationary probability measure is out of reach. As an example, we study an model where users take bikes from nonempty stations and return them directly to non-saturated stations. For this new model, say model II, the same results are only partially proved. However its homogeneous version is completely understood.

Impact on performance. These results allow us to study the performance of large bike sharing systems. The limiting probability that a station is empty or full is taken as the performance measure. In the homogeneous case, the optimal proportion of bikes per station is slightly more than half the capacity of a station. Moreover performance is very close for both models I and II for a large range of values around this optimum. When the system is inhomogeneous, there is a proportion of bikes per station s, depending on each cluster, which minimizes the ratio of problematic stations inside a cluster. However, even if, for a given proportion s, performance is not bad in one cluster, the ratio of problematic stations in the other clusters can be very high, and overall performance is then very poor. The main conclusion of the study is that performance collapses due to heterogeneity.

Related work. As noted above, studies of such systems are rare. Bike sharing systems can be modelled as closed stochastic networks since the fleet is constant. If the capacity of stations is infinite, they are simply closed Jackson networks where the nodes are a mixture of one-server queues for the stations and infinite-server queues for the routes. The work by Malyshev and Yakovlev [START_REF] Vadim | Condensation in large closed Jackson networks[END_REF] on closed Jackson networks with one-server queues can be extended in this case. The limiting behavior of the network, namely condensation in some stations while the others exhibit a queue length with geometric distribution, can be similarly obtained. The technical arguments developed in [START_REF] Vadim | Condensation in large closed Jackson networks[END_REF] are analytical: from the product form of the stationary invariant measure for fixed size, a saddle point method is used (see [START_REF] Flajolet | Analytic Combinatorics[END_REF] for an introduction) to obtain the asymptotic of the partition function. Note that George and Xia [START_REF] George | Asymptotic analysis of closed queueing networks and its implications to achievable service levels[END_REF], in the context of bikesharing systems, provide a simpler proof in a more elementary framework: a case where the number of stations is fixed while the total fleet becomes large.

Mean field limits are applied to study a large class of models in statistical physics or applied probability. It has been also introduced in analysis of algorithms (see e.g. [START_REF] Mitzenmacher | The Power of Two Choices in Randomized Load Balancing[END_REF][START_REF] Mitzenmacher | Studying balanced allocations with differential equations[END_REF]). An issue raised in [START_REF] Tibi | Metastability in communications networks[END_REF] (see also [START_REF] Dupuis | On the construction of lyapunov functions for nonlinear markov processes via relative entropy[END_REF]) is to give an interpretation of Lyapunov functions related to relative entropy in some communication networks. It is pointed out in [START_REF] Tibi | Metastability in communications networks[END_REF] that the special feature of the model is a strong requirement for obtaining the Lyapunov function. Our model has a similar behavior and provides new insights on that topic.

Outline of the paper. Section 2 gives model I description and the main result. Section 3 is devoted to its proof, based on a mean field approach. Model II is described and investigated in Section 4 and practical results about performance of both models are presented in Section 5.

Model I: description and main result

This section gives a description of model I and the main result. This model has reversible Markovian dynamics which allows to obtain the invariant distribution.

We consider a bike sharing system with N stations and a fleet of M bikes. Each station i has a capacity K i,N . The dynamics of the system is as follows. Users arrive at the stations according to independent Poisson processes with rate λ i,N at station i. When a user arrives at a station, if there is no available bike, then he leaves the system. Otherwise, he takes a bike and chooses station j with probability p j,N . The travel time has an exponential distribution with parameter µ N , whatever the station he comes from. When he arrives at station j, if there are less than K j,N bikes in this station, he returns his bike and leaves the system. If there are K j,N bikes (i.e. the station is full), the user chooses again a station say k with probability p k,N and goes to this station. As before, it takes a time with exponential distribution with parameter µ N . He rides like this again until he can return his bike.

This model is a generalization of the homogeneous case where the arrival rates do not depend on the stations and where the user returns the bike at a station chosen at random. Even if the different probabilities to go from a station to another do not appear exactly, as a routing matrix, the popularity of a station is taken into account.

Markovian reversible dynamics

Let N be fixed and X i,N (t) be the number of bikes at station i, 1 ≤ i ≤ N . The process (X N (t)) = (X i,N (t)) 1≤i≤N is an irreducible Markov jump process on the finite state space

S N = {x = (x 1 , . . . , x N ) ∈ N N , 1 ≤ x i ≤ K i,N , N i=1 x i ≤ M }
thus it has a unique invariant probability measure π N . The process is reversible, so it is easy to obtain π N . Let R i,N = µ N p i,N /λ i,N be the utilization at station i. For x in S N and |x|

= N i=1 x i , π N (x) = 1 Z N N i=1 R xi i,N (M -|x|)! (1)
where the normalizing constant is

Z N = x∈S N ( N i=1 R xi i,N )/(M -|x|)!.

Main result

It gives an asymptotic of the invariant measure π N for N and M large, with M/N tending to s > 0. Let r i,N = R i,N /max i R i,N be the relative utilization at station i. The following assumption is crucial:

Assumption W . There exists I a probability measure on ]0, 1] × N and Λ > 0 such that, as N tends to infinity,

1 N N i=1 δ r i,N ,K i,N (w)
-→ I and N R max,N → Λ -1

where w means the weak convergence of measures and R max,N = max i R i,N . Note that this assumption is fundamental for the convergence to the solution of the differential equations and is related to the system topology (see for illustration the examples above).

Theorem 1 Under assumption (W ), if I has a finite support, as N and M tend to infinity with M/N tending to some s, the number of bikes at a station with parameters r and K has a geometric distribution ν ρr,K on {0, . . . , K} with parameter ρr where ρ is the unique solution of

s = ρΛ + ]0,1]×N m(ν ρr,K ) dI(r, K) (2)
where the mean of a random variable with distribution µ is denoted by m(µ).

Theorem 1 is proved in Section 3. Let us first give two examples.

Examples

To simplify, λ i,N , K i,N , µ N and N p i,N do not depend on N . It can be straightforwardly generalized.

The homogeneous case. For each i, K i = K, λ i = λ and p i,N = 1/N thus r i,N = 1. Assumption (W ) is obvious with I = δ 1,K and Λ = λ/µ. Equation (2) rewrites

s = λ µ ρ + m(ν ρ,K ). (3)
The 2) can be rewritten

s = ρΛ + C c=1 α c m(ν ρrc,Kc ) (4)
which generalizes equation (3) obtained for the homogeneous case.

A mean field analysis

Subsection 3.1 is devoted to the proof of Theorem 1, divided into several steps. Some connections to an analytical method are given in 3.2.

Mean field limit

A simplification. For fixed N , the invariant measure π N given by equation ( 1) depends on the different parameters only by R i,N . Thus, to obtain the limiting invariant measure, one can assume without loss of generality that µ N = 1 and a uniform routing p i,N = 1/N . In what follows, we assume that K i,N = K N to avoid cumbersome notations but it remains valid replacing r by (r, K) and proves Theorem 1.

The empirical measure process. Due to heterogeneity with respect to the N stations, the quantity of interest is the following empirical distribution of the stations as function of time

Y N (t) = 1 N N i=1 δ (r i,N ,X i,N (t)) .
The empirical measure process (Y N (t)) is a Markov process on the finite state space

Y N = {y = 1 N N i=1 δ (r i,N ,ni) , (n i ) 1≤i≤N ∈ N, n i ≤ K i,N , N i=1 n i ≤ M },
included in Y = P(]0, 1] × N), the set of probability measures on ]0, 1] × N, with transitions

y → y + 1 N (δ (r,n-1) -δ (r,n) ) 1 rR max,N y(r, n)1 n>0 y → y + 1 N (δ (r,n+1) -δ (r,n) ) N y(r, n)(M - (r ,n ) n y(r , n ))1 n<K
The first transition is an arrival at a station with utilisation r and with n bikes, n > 0. The second one is a return of a bike in a station with utilisation r and with n bikes, n < K.

Remark. The first marginal of Y

N (t) is Y N 1 (t) = 1 N N i=1 δ r i,N
which is deterministic. Thus, (W ) is necessary for the convergence of process Y N .

Convergence to the dynamical system. Under assumption (W ), as I has a finite support, by the same arguments as for an homogeneous model (see for example Darling and Norris [START_REF] Darling | Structure of large random hypergraphs[END_REF]), Y N converges in distribution to a dynamical system, unique solution of the following differential equation

ẏ = ]0,1]×{0,...,K} Λ r (δ (r,n-1) -δ (r,n) )1 n>0 + 1 n<K (δ (r,n+1) -δ (r,n) )(s -m(y 2 )) dy(r, n) (5)
where m(y 2 ) = ndy(r, n). Equilibrium point for the dynamical system. If ȳ in Y is an equilibrium point for the dynamical system, applying the right-hand side of equation ( 5) to A × {p}, A ⊂]0, 1], p ∈ {0, . . . , K},

0 = (1 p>0 ȳ(A × {p -1}) -1 p<K ȳ(A × {p})) (s -m(ȳ 2 )) + 1 p<K A×{p+1} Λ r dȳ(r, n) -1 p>0 A×{p} Λ r dȳ(r, n).
It means that, for r > 0, ȳ(r, .) is the invariant measure of a M/M/1/K queue with arrival to service rate ratio r(s -m(ȳ 2 ))Λ -1 . Thus, given r, ȳ(r, .) is a geometrical distribution on {0, . . . , K} with parameter rρ(ȳ) where ρ(y) is defined on Y by

ρ(y) = (s -m(y 2 ))Λ -1 . ( 6
)
Recall that the geometric distribution on {0, . . . , K} with parameter a is denoted by ν a,K . One gets that, for the equilibrium point ȳ,

m(ȳ 2 ) = m(ν rρ(ȳ),K )dI(r)
which, plugging it in equation ( 6), gives

ρ(ȳ) = (s -m(ν rρ(ȳ),K )dI(r))Λ -1 .
Thus, to prove the existence and uniqueness of ȳ, it is sufficient that the following equation

s = Λρ + m(ν rρ,K )dI(r)
which is equation ( 2), has a unique root ρ > 0. It is straightforward that the right-hand side of this equation is an increasing function of ρ on R + * , from 0 to +∞. This guarantees the existence and the uniqueness of the solution.

Limiting invariant measure. The purpose then is to prove that this equilibrium point is the concentration point as N tends to infinity of the invariant measure of process Y N . A Lyapunov function g of the dynamical system (y(t)) is a function which decreases along its trajectories, more precisely, for each t ≥ 0, d dt g(y(t)) ≤ 0 with equality only if y is the unique equilibrium point of (y(t)). Finding a Lyapunov function is a way to prove that the equilibrium point is the concentration point of the invariant measure of process Y N (see details for example in [START_REF] Tibi | Metastability in communications networks[END_REF]). The following proposition provides a Lyapunov function. Let us recall some definition. For µ and µ probability measures on {0, . . . , K}, H(µ/µ ) = K n=0 µ n log(µ n /µ n ) is the relative entropy of µ with respect to µ . Theorem 2 Function g defined on Y by g(y) = h(y) -S(ρ(y)), where ρ(y) is defined by equation [START_REF] Gast | A mean field model of work stealing in large-scale systems[END_REF],

S(u) = uΛ -s log u + ]0,1]×N log 1 -(ur) K+1 1 -ur dI(r, K), (7)
h(y) = ]0,1]×N
H(y(r, K, .)/ν rρ(y),K )dI(r, K), and ν a,K is defined in Theorem 1, is a Lyapunov function for the dynamical system (y(t)). Scketch of the proof. The main argument is the decreasing of the relative entropy along the trajectories of a Markov process. Here, the dynamical process can be interpreted as a Markov process but with transitions depending on time. This adds a supplementary term, deriving the relative entropy, which is the second term of the expression of g. Remark. Theorem 1 could be extended when I has no more a finite support. The previous approach would work, with further technical arguments for proving the convergence to the dynamical system, and also that g is a Lyapunov function.

parameter rρ, where r = pµ/λ and ρ is a root of equation ( 9). Thus, using the closed-form expressions for ν a,K (0) and ν a,K (K), ρ satisfies:

s = m(ν rρ,K )dI(r) + λ 1 -(rρ) K 1 -(rρ) K+1 I(dx) / λr 1 -(rρ) K 1 -(rρ) K+1 I(dx) . (11)
For the homogeneous case, it is then straightforward that equation ( 11) has a unique root ρ because, in this case, its right-hand side has a first term which reduces to a constant, 1/r, thus is increasing from 0 to K when ρ goes from 0 to +∞. It is conjectured in the cluster case.

Limiting invariant measure. For the homogeneous model, we refer to [START_REF] Fricker | Incentives and regulations in bike-sharing systems with stations of finite capacity[END_REF] for details. The key argument is that m(ȳ) = s -λ/µ is constant, which is no more valid even in the cluster case. 2

Performance

Performance of the system is measured by the stationary probability that a station is problematic, i.e. has no bikes, which means the loss of arriving users, or saturated, i.e. with no lockers available, so users cannot return their bike. For model I, by Theorem 1, when the system gets large, if the station has relative utilization r and capacity K, this probability tends to ȳ(r, K, 0) + ȳ(r, K, K). Recall that ȳ(r, K, .) = ν ρr,K where ρ is the root of equation ( 2). Model II is also detailed.

Model I in the one cluster case

Proposition 1 Consider model I in the one cluster case.

(i) The limiting proportion of problematic stations as a function of s is given by the parametric curve ρ →(s(ρ), ν ρ,K (0

) + ν ρ,K (K)) = λ µ ρ + m(ν ρ,K ), 1 -ρ 1 -ρ K+1 (1 + ρ K ) .
(ii) It is minimal when s = K/2 + λ/µ and the minimum is equal to 2/(K + 1). (iii) As K grows, the performance around s = K/2 + λ/µ becomes flatter and insensitive to λ/µ.

Proof: Assertion (i) is a simple consequence of Theorem 1. For assertion (ii), the intuitive argument is that the behavior of this geometric distribution ν ρ,K for the number of bikes per station depends on the position of ρ with respect to 1: more concentrated on small values if ρ is less than one, on large values, up to K, if ρ is greater than one. Indeed, it can be proved that ϕ(ρ) = ν ρ,K (0) + ν ρ,K (K) has a minimum at 1 which is 2/(K + 1). And thus s(1) = K/2 + λ/µ is the optimal proportion of bikes per station. For (iii), it can be proved that, for

ψ = ϕ • s -1 , ψ (s(1)) = K(K -1) 3(K + 1)(λ/µ + K 2 /12 + K/6) 2 ≥ 48 K(K + 2) 2 .
It gives that the performance curve can never be sharp at the optimal value, for the range of values considered here, i.e. K ≥ 20. See [START_REF] Fricker | Incentives and regulations in bike-sharing systems with stations of finite capacity[END_REF] for details.

2

In practical situations, a system with 30 lockers per station would lead to a proportion of problematic stations of at least 2/31 ≈ 6.5%.

Model I in the cluster case

When we have multiple clusters, the proportion of problematic stations in a cluster c can be expressed as a function of ρ and is equal to

ν ρrc,Kc (0) + ν ρrc,Kc (K c ) = 1 -ρr c 1 -(ρr c ) Kc+1 (1 + (ρr c ) Kc ).
The study of the one cluster case gives that this quantity has a minimum 2/(K c +1) for ρr c equal to 1 and, in this case, the proportion of problematic stations in cluster c is exactly ν r c /rc,K c (0)+ν r c /rc,K c (K c ).

As for the one cluster case, the proportion of problematic stations as a function of s can be plotted using the parametric curve given by, for ρ > 0,

ρ → Λρ + C c=1 α c m(ν ρrc,Kc ), C c=1 α c 1 -ρr c 1 -(ρr c ) Kc+1 (1 + (ρr c ) Kc ) . (12)

Model II

When studying the equilibrium point of the system, the only difference with model I studied in Sections 2 and 3 is that ρ is given by equation ( 11), conjectured to have a unique solution. Therefore, we can plot easily the performance of the system by a parametric curve similar to (12) but where s(ρ) in ( 12) is given by equation [START_REF] Mitzenmacher | Studying balanced allocations with differential equations[END_REF]. For the one cluster case, it means simply to replace λ/µρ by λ/µ in s(ρ) and, moreover, when s ∈ [λ/µ, K + λ/µ], the proportion of problematic stations is one. It shows that the proportion of problematic stations has a similar shape for both models. Indeed, as before, the minimal proportion is 2/(K + 1) and is obtained for s = K/2 + λ/µ. When s is not equal to K/2 + λ/µ, this proportion will be higher than in model I.

This shows that although forcing people to go to a non-saturated or non-empty station reduces the unhappy users since everyone can take or leave a bike at anytime, it makes the system more congested and does not improve overall performance. A user that is not aware of the state of the system will have more problems in this system than in the first one.

An example.

This is illustrated for a two cluster case where the performance for both models is compared. The two clusters have the same size N/2, the travel times have the same exponential distribution with parameter 1, the probability to join each station p 1,N = p 2,N = 1/N and the number of slots K 1 = K 2 = 30 is the same but the arrival rates are λ 1 = 1.3, respectively 2, while λ 2 = 1. Figure 1 is obtained. In Figure 2, the proportion of empty and full stations is detailed. For example, it can be seen, from Figure 2, that even when the system is close to symmetry (λ 1 /λ 2 = 1.3), if the number of bikes is optimal, near 13, a user going from an underloaded to an overloaded station has more than 25% chance of not finding a bike and more than 10% of finding a saturated station. This example shows that even slight asymmetry implies poor system performance.

Conclusion and future work

By extending the mean field methods for a model with clusters, we obtained performance of an inhomogeneous bike sharing system. Our argument would apply to other large networks. To make the model analytically tractable, the way to return the bike then is simplified and we do not take into account geometry. For example, in a real system, a user finding a saturated station would search for an available spot in a local neighborhood. Our simulations on more realistic models, not presented here, suggest that the performance is very close. Future work would include evaluation of the influence of geometry. Moreover, due to intrinsically poor performance, there is actually a pressing need to find efficient algorithms based on incentives to the users or relying on redistribution of bikes by trucks. This question has been investigated in [START_REF] Fricker | Incentives and regulations in bike-sharing systems with stations of finite capacity[END_REF] for an homogeneous system. Designing algorithms for inhomogeneous systems and studying then the system behavior is a challenging problem.

  cluster case. Assume that there are C clusters and each cluster c has N c stations with N c /N which tends to α c . Each station in cluster c has the same parameters λ c , K c and p c,N = β c /N . Thus each station in cluster c has the same utilization R c,N with N R c,N = µβ c /λ c . Assumption (W ) holds with I = C c=1 α c δ rc,Kc , with r c = Λµβ c /λ c where Λ = 1/ max c (µβ c /λ c ). Equation (

2 Fig. 1 :

 21 Fig.1: Performance of the system as a function of the number of bikes in the city with two areas that have different arrival rates for model I (solid) and model II (dash).

2 Fig. 2 :

 22 Fig. 2: Proportion of empty and full stations as a function of the number of bikes in a city with two areas that have different arrival rates for model I (solid) and model II (dash).

The asymptotic of the partition function

As the invariant measure has an explicit form, an analytical method can be considered. It consists in a saddle point method to obtain the asymptotic for Z N when N gets large with M/N which tends to s. If the capacity of the stations is infinite, the asymptotic was given in the paper by Kogan [8] , mentioned also in [9, Section 6] where a result similar to Theorem 1 appears.

Let Z N be indexed also by M . Substituting equation for Z M,N ,

which can be rewritten Θ N (u)/u M = e N S N (u) , where

By assumption (W ), S N (u) converges to S(u) given by equation [START_REF] George | Asymptotic analysis of closed queueing networks and its implications to achievable service levels[END_REF]. By derivation and some algebra,

Thus, S (u) = 0 if and only if u is a root of equation [START_REF] Darling | Structure of large random hypergraphs[END_REF]. A saddle point method will give that, as

It is worth noting the presence of S both in the partition function asymptotic and the Lyapunov function.

Model II

In this section, a second model is proposed. In model II, the total arrival process to the system is exactly the same as in Model I: a Poisson process with rate N i=1 λ i,N , but users arrive only at non-empty stations, with a rate at a non empty station i proportional to λ i,N . Thus he always find a bike and takes it for an exponential time with parameter µ N , then returns it at a non saturated station j, with probability proportional to p j,N . It takes into account the case where the users have information about the problematic stations and avoid them.

In model II, the Markovian dynamics are no more reversible. Moreover, the invariant distribution is out of reach and does not in general depend only on the utilization of the stations. But a result, similar to Theorem 1, is only partially obtained for model II. The idea is , to introduce in the empirical measure a deterministic coordinate more general than the utilization, to obtain the behavior of the system as it gets large by a mean field approach. What can be proved is given by the following proposition.

Theorem 3 Assume that it exists I probability measure on R 3 × N such that

If I has a finite support, as N and M tend to infinity with M/N tending to some s, the empirical measure process converges to a dynamical process. There exists an equilibrium point ȳ and for x = (λ, p, µ, K), ȳ(x, .) has a geometric distribution with parameter rρ, where r = pµ/λ and ρ is a root of

In the homogeneous case, the equilibrium point is unique and is the concentration point of the sequence of the invariant measures of the number of bikes at a station.

Remark. The uniqueness of the equilibrium point and the concentration of the invariant measures are conjectured for the cluster case (i.e. when I has a finite support).

Proof:

Non reversible dynamics. The empirical measure process (Y N (t)) to consider is defined by

where x i,N = (λ i,N , N p i,N , µ N , K i,N ) is called the parameter of station i. It is also a Markov process, on a finite state space Y N included in Y = P(R 3 × N) with transitions, if x = (λ, p, µ, K),

The first transition is an arrival at a station with parameter x and with n bikes, n > 0. The second one is a return of a bike in a station with parameter x and with n bikes, n < K.

Under assumption (W ), for I with finite support, standard arguments give that the empirical measure process converges to the unique solution of Equilibrium point. Let ȳ be an equilibrium point of the dynamical system given by equation [START_REF] Mitzenmacher | The Power of Two Choices in Randomized Load Balancing[END_REF]. As in Section 3.1, it can be proved that, for x = (λ, p, µ, K), ȳ(x, .) has a geometric distribution with