Alain Billionnet

Marie-Christine Costa

Pierre-Louis Poirion
email: plpoirion@gmail.com

Robust optimal sizing of an hybrid energy stand-alone system *

This paper deals with the optimal design of a stand-alone hybrid system composed of wind turbines, solar photovoltaic panels and batteries. To compensate for a possible lack of energy from these sources, an auxiliary fuel generator guarantees to meet the demand in every case but its use induces important costs. We have chosen a two-stage robust approach to take account of the stochastic behavior of the solar and wind energy production and also of the demand. We seek to determine the optimal system, i.e. the one that generates a minimum total cost when the worst case scenario relating to this system occurs. We use a constraint generation algorithm where each sub-problem (the recourse problem) can be reformulated by a mixed-integer linear program and hence solved by a standard solver. We also propose a polynomial time dynamic programming algorithm for the recourse problem and show that, in some cases, this algorithm is much more efficient than mixed-integer linear programming. Finally, we report computational experiments on instances constructed from real data, that show the efficiency of the proposed approach and we study the addition of constraints linking the uncertainty in consecutive time periods.

Introduction

The fast development of renewable energies brought new complex problems of combinatorial optimization in particular as regards the autonomous hybrid energy systems. These systems involve several energy sources as wind, biomass or sun and are not connected to the grid. They are particularly useful in islands as La Réunion [START_REF] Praene | Renewable energy: Progressing towards a net zero energy island, the case of reunion island[END_REF] and in remote areas as there are in Canada [START_REF] Hertz | Optimizing the design of a wind farm collection network[END_REF]. Mathematical techniques have been used to optimize either the operation of the system [START_REF] Ai | Computer-aided design of pv/wind hybrid system[END_REF][START_REF] Gupta | Modelling of hybrid energy system[END_REF][START_REF] Morais | Optimal scheduling of a renewable micro-grid in an isolated load area using mixedinteger linear programming[END_REF][START_REF] Zhang | Scheduling electric power production at a wind farm[END_REF] or the design of the park which this paper is dealing with. The most recent references on this topic include [START_REF] Diaf | A methodology for optimal sizing of autonomous hybrid pv/wind system[END_REF], [START_REF] Ferrer-Marti | A milp model to design hybrid wind-photovoltaic isolated rural electrification projects in developing countries[END_REF], [START_REF] Hertz | Optimizing the design of a wind farm collection network[END_REF], [START_REF] Mustakerov | Wind turbines type and number choice using combinatorial optimization[END_REF] and [START_REF] Prasad | Optimization of integrated photovoltaic-wind power generation systems with battery storage[END_REF]. In [START_REF] Diaf | A methodology for optimal sizing of autonomous hybrid pv/wind system[END_REF], a methodology is introduced to perform the optimal sizing of an autonomous photovoltaic/wind system. In [START_REF] Ferrer-Marti | A milp model to design hybrid wind-photovoltaic isolated rural electrification projects in developing countries[END_REF], the authors present a mathematical programming model to optimize the design of hybrid wind-PV systems that solves the location of the wind-PV generators and the design of the microgrids. In [START_REF] Hertz | Optimizing the design of a wind farm collection network[END_REF], integer linear programming is used to optimize the design of wind farm collection networks. In [START_REF] Mustakerov | Wind turbines type and number choice using combinatorial optimization[END_REF] the aim is to determine the types, numbers and placement of wind turbines to install, considering investment costs and power production criteria; a mixed-integer nonlinear model is proposed and tested. In [START_REF] Prasad | Optimization of integrated photovoltaic-wind power generation systems with battery storage[END_REF], an integrated photovoltaic/wind system with battery storage is considered: a heuristic approach is proposed to find the sizes of the wind farm, photovoltaic array and battery.

In this paper, we study a stand-alone hybrid system composed of wind turbines, solar photovoltaic (or PV) panels and batteries. To compensate for a lack of energy from these sources, an auxiliary fuel generator guarantees to meet the demand in every case but its use induces important costs. We consider here that the type of wind turbines, PV panels, and batteries have been predetermined by the users according to the lands on which the park will be settled. The aim is to determine the optimal number of photovoltaic panels, wind turbines and elements of battery to install in order to serve a given demand while minimizing the total cost of investment and use.

Moreover, the stochastic behavior of the solar and wind energy production on the one hand, and the demand on the other hand, needs to search for a robust solution, i.e. a solution which is good enough whatever the scenario that occurs (see [START_REF] Gabrel | Recent advances in robust optimization: An overview[END_REF]). We assume that there is no known probability distributions of the data and following the approach proposed in [START_REF] Bertsimas | The price of robustness[END_REF] and [START_REF] Thiele | Robust linear optimization with recourse[END_REF], we consider that the uncertain data can vary between given bounds and that there are limits to the total variation of each kind of data. We propose a mixed-integer program in two stages to model the problem: investment decision variables also called here-and-now variables must be fixed in the first stage while operating recourse variables also called wait-and-see variables will be determined once the uncertainty has been revealed. Then we follow the approach proposed in [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF] to solve the problem.

In Section 2 we describe the system operation and the notations. In Section 3 we give a mixed-integer program for the problem without uncertainty. In Section 4 we propose a model in two stages based on mixed-integer programming where the decision variables are integer and the recourse variables are real. In Section 5 we show that, in this case the recourse problem, i.e. the second stage problem, can be solved in polynomial time by using dynamic programming, contrary to the general case presented in [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF]. In Section 6 we propose an exact approach based on constraints generation to solve the robust problem. In Section 7 we test our method on real instances obtained in [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF]. In the last section, we take into account some dependencies between the uncertainty in consecutive time periods before concluding.

System operation

We have to design an energy system for a period spanning many years. The optimization model we propose focuses on one year which is decomposed in T time periods of one hour, where a time period t goes from time t -1 to time t. Considering that the main part of the data depends on the climate, this allows to take into account the variations of weather during the year. As explained before, the types of wind turbines and PV panels that will be installed are predetermined. They are defined by their expected nominal output power, respectively E w t and E p t (in Kw.h) for each time period t, which are functions of the characteristics of the equipments, the land where they will be installed and the mean meteorological data over the past few years for each time period. The costs of a wind turbine and a PV panel are denoted by C w and C p respectively: these costs include purchase and installation costs (reduced to one year according to the lifetime of equipments) and annual maintenance cost including the lease of land.

The purchase and maintenance of the diesel generator induce a fixed cost which is not involved in the optimization. However, its use induces a cost proportional to the energy it provides: let C g denote the unit cost (for 1 Kw.h).

The system is described in Figure 1. When weather conditions are favorable, the energy produced by wind turbines and PV panels is sufficient to serve the demand and the excess energy can be used to charge the battery. In case of unfavorable weather, the energy stored in the battery is used to serve the demand. It is only when battery is empty that the fuel generator is used: the cost-in-use of the generator is very high but it allows to meet the demand in any case.

The battery bank is composed of elements connected together. The cost of one element is denoted by C b and it has a minimal load denoted by K min and a (maximal) capacity denoted by K max : w.l.o.g., in the following we consider Figure 1: A stand-alone energy park that, for each element, the minimal load is equal to 0 and the capacity is equal to K = K max -K min . We assume that the initial load of each element is equal to 0. Each element has a maximum charge per hour, denoted by E in , and a maximum discharge per hour, denoted by E out , in Kw.h. The battery bank operating induces a loss of energy. The rate of return is denoted by γ < 1: for 1 Kw.h charged in the battery, only γ Kw.h can be effectively used. The physical constraints impose that either all the elements are charging (or inactive) or all are discharging; notice that, since the rate of return is less than 1, it is easy to prove that this battery bank operating is optimal. Then, if there are x b elements in the battery bank, we can consider that there is only one battery with capacity equal to x b K, with a maximum charge per hour equal to x b E in and a maximum discharge per hour equal to x b E out . As in [START_REF] Morais | Optimal scheduling of a renewable micro-grid in an isolated load area using mixedinteger linear programming[END_REF], we assume that the state of the battery does not change during a time period: either the battery is charging (or inactive) or it is discharging.

From the size of the lands, the characteristics of equipments and their placement on the land, the maximum number N w max of wind turbines, N p max of PV panels and N b max of elements in the battery bank which can be installed have been calculated. The expected demand in energy has also been evaluated from the data of previous years for each time period and is denoted by D t , t = 1, ..., T .

The problem without uncertainty

If there is no uncertainty, the objective is to determine the design of the park, that is the number of wind turbines, PV panels and elements in the battery bank, in order to meet the demand with a minimal global cost. We propose a mixed-integer model where variables x w , x p and x b are respectively the (integral) number of wind turbines, PV panels and elements in the battery to install in the park, and for t = 1, ..., T , variable e g t denotes the amount of energy produced by the generator during the time period t (from t -1 to t), e in t (resp. e out t ) denotes the amount of energy being charged (resp. discharged) in the battery during the time period t and for t = 0, ..., T , e b t denote the load of the battery at time t. The problem can be written as the following mixed-integer linear program:

LP min x,e C p x p + C w x w + C b x b + C g T t=1 e g t E p t x p + E w t x w -e in t + γe out t + e g t ≥ D t , t = 1, ..., T (1)
e in t ≤ x b E in , t = 1, ..., T (2)
e out t ≤ x b E out , t = 1, ..., T (3)
e b t ≤ x b K, t = 1, ..., T (4)
e b t = e b t-1 + e in t -e out t , t = 1, ..., T

x p ≤ N p max (6) x w ≤ N w max (7) x b ≤ N b max (8) x w , x p , x b ∈ N, (5)
e in t , e out t , e b t , e g t ∈ R + , t = 1...T ; e b 0 = 0 [START_REF] Hertz | Optimizing the design of a wind farm collection network[END_REF] In the objective function, C p x p + C w x w + C b x b corresponds to the investment cost as defined in Section 2, and C g T t=1 e g t represents the usage cost. For each time period, constraints [START_REF] Ai | Computer-aided design of pv/wind hybrid system[END_REF] impose that the demand is met: the amount of energy obtained from wind and sun is equal to E p t x p + E w t x w ; -e in t + γe out t + e g t is the energy supplied by the battery (see Remark 1) and the auxiliary generator. Constraints (2), (3) and (4) limit the charge, discharge and load of the battery. Constraints [START_REF] Ferrer-Marti | A milp model to design hybrid wind-photovoltaic isolated rural electrification projects in developing countries[END_REF] gives the load of the battery at time t: it is equal to the load at time t -1 (e b t-1 ), plus the energy charged in the battery (e in t ) minus the energy provided by the battery (e out t ); as explained in Section 2, we suppose w.l.o.g. that e b 0 = 0. Constraints (6), ( 7) and ( 8) bounds the number of installed equipments. Finally, constraints [START_REF] Gupta | Modelling of hybrid energy system[END_REF] and [START_REF] Hertz | Optimizing the design of a wind farm collection network[END_REF] impose that all variables are positive, variables x are integer and variables e are real. [START_REF] Ai | Computer-aided design of pv/wind hybrid system[END_REF] are verified for Ŝ. The value of the objective function and all the other variables are unchanged and Ŝ is an optimal solution of LP . So the model ensures that from any optimal solution we can deduce a solution such that the battery cannot be simultaneously in charge and discharge.

From any optimal solution of LP we can deduce another optimal solution such that e in t e out t = 0 for all t = 1, ..., T . That is why these constraints, implying that the battery is either in charge or in discharge during a time period, may be omitted in LP .

In the following, we denote by P x the polyhedron defined by constraints (6,7,8 9) (on x) and by P e the polyhedron defined by constraints [START_REF] Bertsimas | The price of robustness[END_REF][START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF][START_REF] Diaf | A methodology for optimal sizing of autonomous hybrid pv/wind system[END_REF][START_REF] Ferrer-Marti | A milp model to design hybrid wind-photovoltaic isolated rural electrification projects in developing countries[END_REF][START_REF] Hertz | Optimizing the design of a wind farm collection network[END_REF] (on e).

Here, the problem has only three integer variables and it can be easily solved with a mixed-integer linear programming software. More generally, the problem is N P -hard when there are n sources of energy (see Annex 2) but it can be solved in polynomial time when n is fixed [START_REF] Lenstra | Integer Programming with a Fixed Number of Variables[END_REF].

A robust model

Now, let us consider that a part of the data are uncertain. For sake of clarity, we assume for the moment that there is uncertainty only on the demand which can vary in a given domain D. Our robustness objective is to find a feasible solution (x, e) that minimizes the total cost involved by the worst possible scenario of D in connection with x. We can state the robust problem as the following mathematical program:

RP min x∈P x C p x p + C w x w + C b x b + max d∈D min e∈P e C g T t=1 e g t E p t x p + E w t x w -e in t + γe out t + e g t ≥ d t t = 1, ..., T
Variables x are the decision variables and variables e are the recourse variables. For any feasible x, the following program R(x) is called the "Recourse Program":

R(x) max d∈D min e∈P e C g T t=1 e g t E p t x p + E w t x w -e in t + γe out t + e g t ≥ d t t = 1, ..., T
As in [START_REF] Thiele | Robust linear optimization with recourse[END_REF] or [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF], we define the uncertainty set by:

D = {d ∈ R T + : d t = D t + δ t ∆ t , T t=1 δ t ≤ δ, 0 ≤ δ t ≤ 1, ∀t = 1, ..., T },
where D t , ∆ t and δ are data, δ being integer. Then, the demand d t will vary between D t and D t + ∆ t , for t = 1, ..., T . ∆ t is the maximum variation of d t , δ t represents the uncertainty on d t , and δ fixes a bound to the cumulated variations. Since we only consider the worst scenarios, δ t vary between 0 and 1 (and not between -1 and 1 as it could be expected). If δ ≥ T , then whatever the values of x, the worst-case scenario will be d t = D t + ∆ t for all t; if δ = 0 then there is only one scenario: d t = D t for all t; thus δ will be chosen between 0 and T : the choice of δ is discussed in Section 7. This definition of the uncertainty implies that the total variation of the demand relative to its reference value D is bounded.

We can generalize the model to the case where there is also uncertainty on the generation of solar and wind energies. Let e p t (resp. e w t ) be the uncertain energy produced by PV arrays (resp. wind turbines). Similarly to D, we define the uncertainty sets E p and E w associated to E p and E w as:

E p = {e p : e p t = E p t -φ t Φ t , T t=1 φ t ≤ φ, 0 ≤ φ t ≤ 1, ∀t = 1, ..., T }, and
E w = {e w : e w t = E w t -ω t Ω t , T t=1 ω t ≤ ω, 0 ≤ ω t ≤ 1, ∀t = 1, ..., T }.
As for the demand, since we only consider the worst scenarios, for given uncertainty budgets φ and ω, we only consider cases where e p t ≤ E p t and e w t ≤ E w t . From [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF] the method proposed hereafter to solve the problem for an uncertain demand (right-hand side of contraints) can be extended for uncertain data on the left-hand side of the constraints. For sake of clarity, we present the method only for uncertain demands, but the tests in Section 7 are carried out for uncertainty on e p and e w as well.

For fixed values of x, the recourse problem R(x) is generally a difficult problem [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF]; nevertheless, we show in the next section that it can be solved in polynomial time in our case.

Solving the recourse problem

From [START_REF] Thiele | Robust linear optimization with recourse[END_REF], we know that there is an optimal solution of R(x) such that δ t ∈ {0, 1} for all t = 1, ..., T and we can rewrite R(x) as follows:

R(x) max T t=1 δt≤ δ δt∈{0,1}, t=1,...,T min e∈P e C g T t=1 e g t E p t x p + E w t x w -e in t + γe out t + e g t ≥ D t + δ t ∆ t t = 1, ..., T
Clearly, for any x, the worst scenario is obtained with T t=1 δ t = δ and then it will be determined by setting δ t = 1 for δ periods and δ t = 0 for the others, or equivalently by setting d t = D t + ∆ t for δ periods and d t = D t for the others. The problem now is to select the δ periods for which δ t = 1. We propose a polynomial dynamic programming approach to answer this question and solve the recourse problem.

We have to solve the recourse problem for given values of x, thus

E p t x p +E w t x w is fixed and we denote Dt = D t -E p t x p -E w t x w , dt = d t -E p t x p -E w t x w = Dt + δ t ∆ t . In addition, since x b is fixed, we assumed w.l.o.g. that x b = 1, or equivalently we set K ← x b K, E in ← x b E in and E out ← x b E out .
Since the unit diesel cost does not depend on t, we can determine the battery load at time t (e b t ) and the amount of energy supplied by the generator during period t (e g t ), from the battery load at t -1 (e b t-1 ) and the amount of energy either in excess or in lack ( dt ). We have the following proposition : Proposition 1. For given values of d and β = e b t-1 (0 ≤ β ≤ K), for t = 1, ..., T, we have:

e b t = f ( dt , β) = max 0, β -E out , β - dt γ e g t = g( dt , β) = dt -γ min β, E out , dt γ            if d ≥ 0 e b t = f ( dt , β) = min K, β + E in , β -dt e g t = g( dt , β) = 0    if d < 0
Proof. During time period t,

• If dt ≥ 0, then the energy produced by wind turbines and PV panels is not sufficient to serve the demand. The battery discharges a quantity of energy equal to min β, E out , dt γ , and the generator must provide an amount of energy equal to g( dt , β) = dt -γ min β, E out , dt γ ≥ 0.

• If dt < 0, then the energy produced is greater than the demand and the system stores in the battery an amount of energy equal to min K -β, E in , -dt .

The generator is not used and g( dt , β) = 0.

The battery load at time t, e b t = f ( dt , β), is immediately deduced from the amount of energy charged or discharged from the battery during time period t.

The algorithm operates from τ = T to τ = 1 and considers at each step a "truncated recourse problem" R x (τ, ζ, β) defined on the (T -τ + 1) last time periods. At each step, it looks for an optimal operating for the considered period, i.e. from time τ -1 to time T , in function of the two parameters β and ζ: β is the battery load at time τ -1, i.e. at the beginning of period τ, and ζ is the uncertainty budget for these periods, that is the number of time periods with δ t = 1 among the (T -τ + 1) time periods; then, δ represents the global "uncertainty budget" and the only value of ζ to consider for

τ = 1 is ζ = δ. Then R x (τ,
e in t ≤ E in , t = τ, ..., T (1)
e out t ≤ E out , t = τ, ..., T (2)
e b t ≤ K, t = τ, ..., T (3)
e b t = e b t-1 -e out t + e in t , t = τ, ..., T (4)
e b τ -1 = β (5 ) (5)
e in t , e out t , e b t , e g t ∈ R + , t = τ, ..., T

Notice that (τ, ζ, β) represents the "current state" of the decision process and that

R(x) = R x (1, δ, 0). Let us denote v(R x (τ, ζ, β)) by v(τ, ζ, β). We obtain the following recurrence relation verified by v(τ, ζ, β) for 0 ≤ β ≤ K : v(τ, 0, β) = C g g( Dτ , β) + v(τ + 1, 0, f ( Dτ , β)), for 1 ≤ τ < T , v(T, ζ, β) = C g g( DT + ζ∆ T , β), for 0 ≤ ζ ≤ 1, v(τ, ζ, β) = max(C g g( Dτ , β) + v(τ + 1, ζ, f ( Dτ , β)), C g g( Dτ + ∆ τ , β)+ v(τ + 1, ζ -1, f ( Dτ + ∆ τ , β))), for 1 ≤ τ < T and 0 < ζ ≤ T -τ + 1.
Now, let us study v(τ, ζ, β). We are going to prove that there is a constant B (0 ≤ B ≤ K) such that v(τ, ζ, β) is a function of β linearly decreasing on [0, B] and constant on [B, K] (see Figure 2).

Proposition 2. For any τ, 1 ≤ τ ≤ T , any ζ, 0 ≤ ζ ≤ T -τ + 1, and any β ∈ [0, K],
• there are [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF].

B ≥ 0 and C ≥ 0 s.t. v(τ, ζ, β) = C -γC g min(β, B) (11), • B = v(τ,ζ,0)-v(τ,ζ,K) γC g and C = v(τ, ζ, 0) is a solution of
Proof. The complete proof is given in Annex 1. We just give here the main ideas it uses.

First, the proof is given for ζ = 0: in this case, the problem reduces to a problem without uncertainty since δ t = 0 for all t = τ to T . The proof is by induction on t, from T to T -τ + 1; first we prove the existence of B and C and then we compute them. When ζ = T -τ + 1 we have δ t = 1 for all t = τ to T and Notice that since the computation of the initial states is O(T 2 ), the complexity of the algorithm is O( δ(T -δ+1)+T 2 ). Since δ ≤ T, this algorithm is polynomial.

O × β v(τ, ζ, β) v(τ, ζ, 0)× v(τ, ζ, β)× β × s l o p e : -γ C g v(τ, ζ, K) × K ×

Solving the robust problem

The robust problem can be rewritten as:

T T δ τ ζ • • • • • • • • • • • • • • • • • • β = 0 T T δ τ ζ × × • • • • • • • • • • • • • • • • • • β = K • : initial states × • : intermediate states • : final state

RP min

x C p x p + C w x w + C b x b + v(R(x)) x p ≤ N p max x w ≤ N w max x b ≤ N b max x b , x w , x p ∈ N
where v(R(x)) is the optimal value of R(x) which can be obtained in polynomial time (see Section 5 ). Since in addition, v(R(x)) is a decreasing function of x, we could use a branch and bound algorithm to solve RP . Nevertheless, we prefer the approach proposed in [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF] and [START_REF] Thiele | Robust linear optimization with recourse[END_REF] which proved efficient in practice. For any x ∈ P x and d ∈ D, we define the following linear program: R(x, d)

min e∈P e C g T t=1 e g t -e in t + γe out t + e g t ≥ d t -E p t x p -E w t x w , t = 1, ..., T (1)
Let us associate the dual variables λ to constraints (1) and α, β, µ, π respectively to constraints (2), ( 3), ( 5) and (4) of P e : we dualize R(x, d). Then R(x) becomes a max max instead of a max min problem and we can reformulate R(x) to obtain a quadratic maximization program DR(x) such that v(DR(x)) = v(R(x)):

DR(x) max α,β,δ,λ,µ,π T t=1 [(D t -E p t x p -E w t x w + δ t ∆ t )λ t -x b E in α t -x b E out β t -x b Kπ t ] s.c. λ t ≤ C g , t = 1, ..., T (13) -λ t -α t + µ t ≤ 0, t = 1, ..., T (14) γλ t -β t -µ t ≤ 0, t = 1, ..., T (15) µ t+1 -µ t -π t ≤ 0, t = 1, ..., T (16)
T t=1

δ t ≤ δ ( 18
)
δ t ∈ {0, 1}, t = 1, ..., T (19) α t , β t , λ t , µ t , π t ≥ 0 t = 1, ..., T (17)
As in [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF] we can linearize the quadratic terms in the objective function by substituting the variable ν t to the product δ t λ t . We obtain the following mixedinteger linear program

LDR(x) max α,β,δ,λ,µ,ν,π T t=1 [(D t -E p t x p -E w t x w )λ t + ∆ t ν t -x b E in α t -x b E out β t -x b Kπ t ] s.c. λ t ≤ C g , t = 1, ..., T (13) -λ t -α t + µ t ≤ 0, t = 1, ..., T (14) γλ t -β t -µ t ≤ 0, t = 1, ..., T (15) µ t+1 -µ t -π t ≤ 0, t = 1, ..., T (16)
T t=1 δ t ≤ δ ( 18
)
δ t ∈ {0, 1}, t = 1, ..., T (19) ν t ≤ C g δ t , t = 1, ..., T (20) ν t ≤ λ t , t = 1, ..., T (21) α, β, λ, µ, ν, π ≥ 0 ( 17
)
Let P Q be the polyhedron defined by the constraints (13), ..., (21) of LDR(x) where we replace (19) by 0 ≤ δ t ≤ 1, and let (P

Q ) I = conv(P Q ∩{δ ∈ {0, 1} T }),
be the convex hull of the feasible solution of LDR(x). Notice that this convex hull does not depend on x. (P Q ) I being a polyhedron we can rewrite the robust problem RP as the following linear program:

P ROB min x,z C p x p + C w x w + C b x b + z s.c. z ≥ T t=1 [(D t -E p t x p -E w t x w )λ s t + ∆ t ν s t -x b E in α s t -x b E out β s t -x b Kπ s t ], s = 1, ..., S. x b ≤ N b max , x w ≤ N w max , x p ≤ N p max , z ≥ 0, x b , x w , x p ∈ N
where S = |S| and S = {(α s , β s , λ s , µ s , ν s , π s ) 1≤s≤S } is the set of extreme points of (P Q ) I . However, due to the potentially tremendous number of constraints, we solve P ROB by a constraint generation algorithm as in [START_REF] Gabrel | Robust location transportation problems under uncertain demands[END_REF]. Initially, we consider a subset S 0 of S; at a step k, we consider a subset S k of S and we solve a relaxed program P ROB k of P ROB, called master problem, which consists in solving P ROB with the subset of constraints corresponding to S k . The obtained solution is denoted by (x k , z k ).

Then we solve LDR(x k ), called slave problem, to check if (x k , z k ) is optimal. If not, then a new constraint is added, i.e. an extreme point is added to S k (See Algorithm 1).

On the basis that the number of extreme points of (P Q ) I is finite, one can prove that this algorithm converges in a finite number of steps.

Two approaches can be applied to solve LDR(x): either we directly use Cplex to solve this mixed-integer program or we use the approach presented in Section 5 to solve it in polynomial time. Indeed, by using dynamic programming we obtain in polynomial time the optimal values of R(x), e and d; let d * be the optimal associated value of d. So we have the optimal value of R(x, d * ), which is, by strong duality, the optimal value of DR(x, d * ) : v(DR(x, d * )) = v(DR(x)) = v(LDR(x)). But we also need the values of the dual variables, α, ..., π : there are easy to obtain by applying the complementary slackness theorem to both dual programs R(x, d * ) and DR(x, d * ) which are in standard form. In our experiments, we directly get these values by Cplex.

Algorithm 1 Constraint generation algorithm

1: We denote: x = (x p , x w , x b ) and Cx = c p x p + C w x w + C b x b . 2: (α 0 , β 0 , λ 0 , µ 0 , ν 0 , π 0 ) = 0. Set L ← -∞, U ← +∞, k ← 1.
3: Solve the master problem :

P ROB k min x,z Cx + z s.c. z ≥ T t=1 [(Dt -E p t x p -E w t x w )λ s t + ∆tν s t -x b E in α s t -x b E out β s t -x b Kπ s t ], 0 ≤ s ≤ k -1. x b ≤ N b max , x p ≤ N p max , x w ≤ N w max , z ≥ 0, x b , x w , x p ∈ N. Let (x k , z k ) be the obtained solution. L ← Cx k + z k . 4: Solve LDR(x k ) either by Cplex or by dynamic programming. Let (δ k , λ k , ν k , α k , β k , π k ) be the optimal solution. U ← min{U, Cx k + v(DR(x k ))}. if U = L, then return (x k , z k ) else go to 5. 5: Add the constraint z ≥ T t=1 [(Dt -E p t x p -E w t x w )λ k t + ∆tν k t -x b E in α k t -x b E out β k t -x b Kπ k t ]
to the master problem P ROB k , k ← k + 1 and go to 2.

Results

The proposed wind-PV model was applied for designing off-grid networks from three very different sets of real instances given for the state of Montana (USA) which has a continental climate, an island in the Philippines which has a tropical climate and Dunkerque (France) which has an oceanic climate.

The characteristics of the system components are the same for the three regions and come from the Homer website [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF]. There are small wind turbines of type BWC XL.1 with a rated power of 1.24 kW, PV panels of size 1 kW and batteries Trojan L16P. The annual unit costs C p , C w and C b are computed in the following way: purchase cost lifetime + annual O&R cost (Operations and Regulatory); C g is the cost for producing 1 Kw.h by the generator. The number T of time periods is equal to 8760 which is the number of hours in one year. Table 1 gives the main parameter values. For the three regions, the demand values D t , and the solar and wind production values E p t and E w t , t = 1, ..., T = 8760, are the average values for each hour (or time period) over several passed years. For the Montana and the Philippines, the data come from the Homer website [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF]. The average values for Dunkerque were calculated from climatic data provided by Meteo France [START_REF] France | Méteo climat, histoclimat MCS, billions of data[END_REF] and demand data provided by RTE-France [START_REF]Réseau de transport d'électricité[END_REF]. We do not recall here the 8760 hourly mean values of these parameters. Finally, we bounded the number of PV panels and wind turbines to 120 and batteries to 700. We made several kinds of tests in order to verify the efficiency and the accuracy of our approach. First we study the influence of the demand uncertainty on the design of the park, as well for the parameter fixing the level of uncertainty ( δ) as for the parameter limiting the demand deviation from the mean demand value (∆ t ). Then we compare the dynamic programming approach to the approach based on the use of CPLEX from the computational point of view. Finally, we give some results where the uncertainty concerns both the demand and production by wind turbines and photovoltaic panels. Additional tests are presented in Section 8 to take into account some dependencies between demands in consecutive time periods.

C p C w C b C g K γ E in E out N p max , N

Results in function of the uncertainty level

In this section we give for the three regions the results for different values of δ (see Section 4). For each δ we give the optimal number of wind turbines, PV panels and elements in the battery. We also give the optimal value of the robust problem, that is the minimal cost including the investment cost for the park and the fuel consuming.

Table 2 gives the results for Montana. In the first part, we assume that the maximum deviation from the mean value D t of the demand at each time period is ∆ t = 10%D t , t = 1, ..., 8760. In the second part, we assume that ∆ t = 30%D t which corresponds to the actual deviation for the Montana. As expected, the optimal value of the robust problem increases as a function of δ corresponding to an increase of the demand but it reaches its maximum value for a certain threshold (here 900). Indeed, the optimal value of the robust problem cannot decrease when the global uncertainty level increases since any solution for δ = δ is admissible for δ < δ. Furthermore, once the park is sufficient to cover the demand for some critical set of periods, it is also sufficient to cover any other periods. Therefore the cost remains constant as soon as δ is large enough to cover the critical periods. In addition we note that there is no significant differences between the two cases ∆ t = 10%D t and ∆ t = 30%D t : evidently the cost and the number of elements in the park are slightly larger in the second case since the demand can be larger but there is also a stabilization of their values for a threshold (around 850).

In Tables 3 and4 we give the results for Dunkerque and the Philippines, respectively. We give the results only for ∆ t = 10%D t because the results obtained for ∆ t = 30%D t are very similar. We see that the stabilization is obtained for a lower threshold ( δ = 160) for Dunkerque. On the contrary, for the Philippines we must compute the park until δ = 6000 to obtain a stabilization but the increasing is very slow from δ = 2000 to 6000. The differences between the fuel costs and the thresholds of stabilization in the two regions are probably due to the fact that the climate is much more irregular in the Philippines than in Dunkerque.

Dynamic programming versus integer linear programming

In Table 5, we compare the CPU times required by Algorithm 1 when the recourse problem (LDR(x) at Step 4 of Algorithm 1) is solved either by using the polynomial time dynamic programming algorithm of Section 5 or by using CPLEX. We notice that for intermediate values of δ the dynamic programming approach is much faster than the approach using Cplex. In particular, for values of δ between 400 and 600 for the Montana, the Cplex approach can take several hours to solve the problem while the dynamic programming approach only needs a few minutes. Indeed, there are about 40000 variables including 8000 integer variables, and 50000 constraints. Moreover, solving LDR(x) is difficult because there is an important gap between the values of its continuous relaxation and its optimal integer solution: this gap is in part due to the linearization of the quadratic terms. Furthermore, with Cplex, the CPU time increases until δ = 900 and then decreases until δ = 8760 while it is almost constant for dynamic programming. We verified that the non efficiency of CPLEX in median cases is due to the very big number of nodes explored in the branch and bound. Indeed, median values correspond to difficult instances and extremal values correspond to cases where there is little uncertainty because the actual demands are all close either to the mean value or to the largest value.

Uncertainty in demand and energy production

Now we consider that the demand as well as the energy production are uncertain and we bound the uncertainty level respectively to δ for the demand, ω for the wind turbine energy production, and φ for the photo-voltaic production. From the results presented in [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF] we know that the our approach can be used to solve the problem when ω and φ correspond to the maximum number of periods where the unit productions E w t and E p t can reach their minimal values. But the size of the problems to solve is much more larger since we have to introduce three set of 0-1 variables, z t for the demand (as before), z ω t (for wind turbine production) and z φ t (for photo-voltaic production), t = 1, ..., T : the number of variables and linearization constraints is very large. Moreover, the dynamic programming approach becomes intractable to solve the recourse problem for large values of T , although the algorithm remains polynomial. Nevertheless, we could test the constraint generation algorithm for small values of ( δ, φ, ω) and some significant results are presented in Table 6. Solving the most difficult instances ((300,100,100), (300,300,300) or (500,500,500)) requires more than two hours and it was not possible to solve bigger instances. 6: Montana (∆ t = 10%D t ). Design and cost of the park when considering energy and demand uncertainties .

We do not present all the results but our tests show that, as previously the optimal value of the robust problem increases as a function of δ, φ and ω, until a certain threshold. Furthermore, we notice that the combined influence of uncertainties in the energy generated by wind turbines and PV-panels induces a larger cost augmentation than the sum of the ones induced by considering separate uncertainties in wind and solar energy generation. Indeed, the production of the park must cover the demand for the δ critical periods of demand as for the φ and ω critical periods of production.

Other definitions of demand uncertainty

In this section we only consider demand uncertainty and the uncertainty corresponding to a worst scenario is defined by the values taken by δ t , t = 1, ..., T in the recourse problem. Our tests were made on real world data sets which show some natural dependency between demands D t in consecutive time periods; but in the definition of uncertainty given in Section 4 we assume that the variations δ t on these demands are independent. Doing so, we probably oversize the park and obtain an upper bound of the worst cost: indeed, introducing new constraints linking δ t values will lead up to reduce the uncertainty domain and so to reduce the cost.

Several possible models allow to take into account some dependencies between variables δ t . For instance, we could consider additional constraints as (1 -ρ)δ t-1 ≤ δ t ≤ min((1 + ρ)δ t-1 , 1), ∀t = 2, ..., T , 0 < ρ < 1, where ρ is a fixed parameter which induces a small variation of uncertainty between consecutive time periods. But our approach is not compatible with such constraints because they imply non integer values of δ t in the optimal solutions of the recourse problem. Then R(x) cannot be easily solved: the dynamic programming approach is no more valid, and the dual DR(x) is a non convex quadratic program which cannot be linearized. So, since there is no specific real world justifications for such constraints, we try other ways to track this aspect.

Sliding time window constraints for uncertainty

A first idea to avoid erratic behaviour of the data variation is to forbid a sawtooth oscillation of the δ t values, for instance sequels as 1 0 1 0 1 0 1... . For this purpose, we add a new constraint limiting the total number of gaps 0 to 1 (or 1 to 0) between δ t-1 and δ t : T t=2 |δ t -δ t-1 | ≤ H where H < δ is a fixed integer parameter measuring the "oscillation level". The new definition of D is:

D = {d ∈ R T + : d t = D t + δ t ∆ t , δ t ∈ {0, 1}∀t, T t=1 δ t ≤ δ, T t=2 |δ t -δ t-1 | ≤ H}.
To introduce the constraint T t=2 |δ t -δ t-1 | ≤ H in the recourse problem R(x), we need to linearize it. We add new variables δt ≥ 0; t = 2, ..., T and the constraint becomes:

δt ≥ δ t -δ t-1 ∀t = 2, ..., T (a) δt ≥ δ t-1 -δ t ∀t = 2, ..., T (b) T t=2 δt ≤ H (c)
The recourse problem with these new constraints can still be solved in polynomial time by dynamic programming: we add a new parameter η = 1, ..., H and consider a new recurrence formulation of the recourse problem R x (τ, ζ, β, η), where η measures the oscillation level used in the time periods τ to T . The presentation of the algorithm would be tedious and is not given here. Moreover, the computing time increases greatly and is no more competitive with a linear programming approach using CPLEX, where the new constraints (a), (b) and (c) are added in R(x).

On another hand, too long sequels of time periods with δ t = 1 are unlikely and lead also up to oversize the park; to avoid such scenarios, we can add the constraints:

t 0 +T s t=t 0 δ t ≤ T , ∀t 0 = 1 to T -T s (d)
where T and T s are fixed parameters which limit to T the number of consecutive periods with δ t = 1 in any set of T s consecutive periods.

The costs are somewhat lower than the previous ones but there is no significant differences. In fact, for the main part of the results obtained in Section 7, the worst scenarios defined by δ t were obtained by sequels of 1 followed by sequels of 0, which can be explained as follows: firstly, a worst scenario corresponds to δ t = 1 when the demand at t is high (and the production is low), and an analysis of the data show that there are generally several consecutive such time periods; secondly, the recourse cost increases when the battery is empty and this occurs after several periods of high demands. Moreover, the computing time is much greater with these new constraints. In fact, when H increases, the problem is less constrained, which can explained that the computing time decreases while the cost increases.

We also tested the model with additional constraints (a, b, c, d) but in this case the computing time becomes very high (several hours) and there is no significant differences in the obtained solutions.

Conclusion

We studied a robust two-stage mixed integer optimization problem associated to the design of an hybrid energy system. This problem is characterized by the basic assumptions made in [START_REF] Billionnet | 2-stage robust milp with continuous recourse variables[END_REF] or [START_REF] Gabrel | Robust location transportation problems under uncertain demands[END_REF], i.e. a linear model of the problem without uncertainty, a smooth definition of uncertainty with no given distribution or mean values of the data, an objective which is to minimize the cost of a worst scenario, mixed integer first stage variables and continuous second stage variables. We proved that in this special case Dynamic Programming can be used to solve the recourse problem. Our results could be applied to other problems of the same type as inventory management problems.

The method allows to take into account sliding time window constraints between the variables defining the uncertainty. However, for our problem the dynamic programming approach fails to take into account these new constraints because of the great number of time periods (8760 ). Nevertheless, it could probably be used for smaller problems.

The next objective for solving robust two-stage mixed integer linear optimization problems will be to build an efficient method for integer recourse variables: that is a real challenge for operational research specialists. • there are B ≥ 0 and C ≥ 0 s.t. v(τ, ζ, β) = C -γC g min(β, B) [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF], [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF].

• B = v(τ,ζ,0)-v(τ,ζ,K) γC g and C = v(τ, ζ, 0) is a solution of
Proof. The problem reduces to a problem without uncertainty: if ζ = 0 then δ t = 0 for all t = τ, ..., T. We first prove that for each τ ∈ {1, ..., T } there are two constants B and C such that: v(τ, 0, β) = C -γC g min(β, B) (12), then we verify the expression of B and C. From the recurrence relation (12); we have, for τ < T : v(τ, 0, β) = C g g( Dτ , β) + v(τ + 1, 0, f ( Dτ , β)). We then proceed by induction from T to 1.

• Step 1: τ = T, v(T, 0, β) = C g g( DT , β). If DT < 0, then v(T, 0, β) = 0 (B 1 = 0, C 1 = 0). Else v(T, 0, β) = C g DT -γC g min β, E out , DT γ (B 1 = min E out , DT γ , C 1 = C g DT . ) • Step 2: τ = τ + 1. Assume there are B 2 ≥ 0 and C 2 ≥ 0 such that v(τ + 1, 0, β) = C 2 -γC g min(β, B 2 ). • Step 3: τ = τ , v(τ , 0, β) = C g g( Dτ , β) + C 2 -γC g min(f ( Dτ , β), B 2 ).
-If Dτ ≥ 0, we have: -If Dτ < 0, then:

v(τ , 0, β) = C g Dτ -γC g min β, E out , Dτ /γ + C 2 -γC g min max β -Dτ /γ, β -E out , 0 , B 2 = (C g Dτ + C 2 ) -γC g min β, E out ,
Q = 0 + C 2 -

Annex 2. Complexity of the general problem without uncertainty

Let us consider a general model where there are n sources of renewable energy. We will call LP gen the resulting model. In this case, one unit of source i costs C i and produces a quantity E i t of energy during period t, t = 1, ..., T . The objective function becomes n i=1 C i x i + C b x b + C g T t=1 e g t . Constraint (1) becomes n i=1 E i t x i -e in t + γe out t + e g t ≥ D t , t = 1...T and constraints ( 6)-( 7) are replaced by x i ≤ N i max , i = 1, ..., n. Constraints (2)-( 5), ( 8) and ( 10) are retained. We are going to show that this generalized problem is NP-hard. Thus, there is no polynomial time algorithm solving LP gen unless P = N P . Proof. We will show that the bounded knapsack problem BKP reduces to LP gen . Consider the decision problem DBKP associated with BKP : n i=1 a i x i + a n+1 x b ≥ V . Thus, the vector y defined by y i = x i for i = 1, ..., n, y n+1 = x b is a solution of DBKP , and the two problems are equivalent.

Since DBKP is NP-complete [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF], we have proved that the decision problem DLP gen associated with the generalized problem LP gen is also NP-complete.

Aknowledgment

We thank the referees for their very careful reading and for all their justified remarks which allowed us to improve the paper.

Figure 2 :

 2 Figure 2: Variations of v(τ, ζ, β) in function of β

Figure 3 :

 3 Figure 3: Pattern of calculation of v(τ, ζ, β).

x p x w x b

  Cost δ = φ = 0 64 45 467 49874 δ = 0, φ = 0, ω = 500 66 48 465 51348 δ = 8760, φ = 0, ω = 0 71 50 504 54884 δ = 0, φ = 8760, ω = 0 71 48 51852 δ = 0, φ = 0, ω = 8760 64 51 469 51438 δ = φ = ω = 8760 78 54 523 58672 δ = 100, φ = 100, ω = 100 67 47 468 53743 δ = 300, φ = 100, ω = 100 66 48 499 55566 δ = 100, φ = 300, ω = 100 70 46 489 54334 δ = 100, φ = 100, ω = 300 69 39 522 54297 δ = 300, φ = 300, ω = 300 70 49 504 56809 δ = 500, φ = 500, ω = 500 74 49 500 58183 Table

Annex 1 .Proposition 3 .

 13 Proof of Proposition 2First, we give the proof for ζ = 0: If ζ = 0 then for for any τ, 1 ≤ τ ≤ T and any β ∈ [0, K],

+ C 2 - 2 = (C g Dτ + C 2 )-= C 3 -

 2223 Dτ /γ + min max β -Dτ /γ, β -E out , 0 , B 2 By considering the three cases β ≤ min E out , Dτ γ , min E out , Dτ γ < β ≤ B 2 +min E out , Dτ γ , and B 2 +min E out , Dτ γ < β, it is easy to verify that:v(τ , 0, β) = (C g Dτ + C 2 ) -γC g min β, B 2 + min E out , Dτ /γ = C 3 -γC g min(β, B 3 )•Step 2: assume it is true for (τ + 1, ζ -1) and for (τ + 1, ζ), we have∀β ∈ [0, K], v(τ + 1, ζ -1, β) = C 2 -γC g min(β, B 2 ), v(τ + 1, ζ, β) = C 2 -γC g min(β, B 2 ),• Step 3: from the recurrence relation (1) we have:v(τ, ζ, β) = max(C g g( Dτ , β) + v(τ + 1, ζ, f ( Dτ , β)), C g g( Dτ + ∆ τ , β) + v(τ + 1, ζ -1, f ( Dτ + ∆ τ , β))). Let Q = C g g( Dτ , β) + v(τ + 1, ζ, f ( Dτ , β)).-If Dτ ≥ 0:Q = C g Dτ -γC g min β, E out , Dτ /γ γC g min max β -Dτ /γ, β -E out , 0 , B γC g min β, E out , Dτ /γ + min max β -Dτ /γ, β -E out , 0 , B 2As in the proof of Proposition 3, we getQ = (C g Dτ + C 2 ) -γC g min β, B 2 + min E out , Dτ /γ γC g min(β, B 3 ) with B 3 = min B 2 + min E out , Dτ γ , K and C 3 = (C g Dτ + C 2 ).

2 = C 2 -= C 2 -B 3

 2223 γC g min min(β -Dτ , β + E in , K), B γC g min(β -Dτ , β + E in , K, B 2 ) γC g min(β + min(-Dτ , E in ), min(K,B 2 )) = C 3 -γC g min(β, B 3 ) = min(K, B 2 ) -min(-Dτ , E in ) and C 3 = C 2 -γC g min(-Dτ , E in ) if min(K, B 2 ) ≥ min(-Dτ , E in ) B 3 = 0 and C 3 = C 2 -γC g min(-Dτ , E in ) -γC g (min(K, B 2 ) -min(-Dτ , E if min(K, B 2 ) < min(-Dτ , E in ) Let R = C g g( Dτ + ∆ τ , β) + v(τ + 1, ζ -1, f ( Dτ + ∆ τ , β)).We prove similarly that there exists C 3 " ≥ 0, and B 3 " ≥ 0, such thatR = C 3 " -γC g min(β, B 3 ").Therefore,v(τ, ζ, β) = max(C 3 -γC g min(β, B 3 ), C 3 " -γC g min(β, B 3 ")).By considering the four cases:C 3 < C" 3 and B 3 < B" 3 ; C 3 < C" 3 and B 3 ≥ B" 3 ; C 3 ≥ C" 3 and B 3 < B" 3 ; C 3 ≥ C" 3 and B 3 ≥ B" 3 , it is easy to verify that there is B 3 , B 3 ∈ [B 3 , B" 3 ], such that: v(τ, ζ, β) = C 3 -γC g min(β, B 3 ), with C 3 = max(C 3 , C" 3 ).we conclude the recurrence referring to the scheme shown on figure 2. Indeed, by proposition 3, the proposition is true for (τ, ζ) ∈ {(T, 0), (T, 1)}, and it is easy to verify that if the proposition is true for (τ, ζ) ∈ {(τ + 1, 0), (τ + 1, 1), ..., (τ + 1, T -(τ + 1) + 1)}, then it is true for (τ, ζ) ∈ {(τ , 0), (τ , 1), ..., (τ , T -τ )} and for (τ, ζ) = (τ , T -τ +1), by proposition 3.

Proposition 4 .

 4 LP gen belongs to the class of NP-hard problems.

  ζ, β) can be written as the following mathematical program where dt = Dt + δ t ∆ t :

				T
	max T t=τ δt≤ζ δt∈{0,1}, t=τ,...,T	min e	C g	t=τ	e g t
			-e in t + γe out t	+ e g t ≥ Dt + δ
	R x (τ, ζ, β)

t ∆ t , t = τ, ..., T

Table 1 :

 1 Main data values.

	w max	N b max

Table 2 :

 2 Montana. Park design and total associated robust cost as a function of the global level of uncertainty δ and the maximal variation on the demand ∆ t .

	δ

Table 3 :

 3 Dunkerque. Results as a function of δ.

	δ

Table 4 :

 4 Philippines. Results as a function of δ.

Table 5 :

 5 Montana and Philippines. Comparison of the two approaches from the CPU time point of view (CPU times in seconds).

	δ

* This work was partially supported by the PGMO Programme Gaspard Monge for Optimization and OR of Fondation Mathématique

Results

First we tested the complemented model obtained by adding constraints (a, b, c) to the definition of D and using CPLEX, on the data provided for Montana tested in Section 7. Some of the results are given in Table 7. with

-If Dτ < 0, then:

Now, let us verify the expression of B and C.

The proof if ζ = T -τ + 1 is similar and is not given here.

The proof for ζ = T -τ + 1 is very similar and is not given here. Just notice that δ t = 1 for all t. We now prove the proposition for any value of (ζ).

• there are [START_REF]Homer energy, energy modeling software for hybrid renewable energy systems[END_REF].

Proof. We prove the proposition by recurrence on (τ, ζ). Now we define the decision problem DLP gen associated with LP gen :

Question: Is there a vector (x 1 , ..., x n ) ∈ N n , an integer x b , and nonnegative reals e in t , e out t , e b t , e g t , t = 1, ..., T, satisfying the constraints of LP gen and such that the value of the objective function is less than or equal to C?

From an instance of DBKP , let us construct the following instance of DLP gen : T = 1; there are n sources of energy with

Notice that we can set e in 1 = 0 and e out 1 = a n+1 , which results in e b 1 = 0 from constraint [START_REF] Ferrer-Marti | A milp model to design hybrid wind-photovoltaic isolated rural electrification projects in developing countries[END_REF]. Then constraints (2), (3), ( 4) and (5) are verified and the remaining problem is to determine if there are (x 1 , ..., x n ) ∈ N n , x b ∈ N, e g 1 ∈ R + such that the following set of constraints is satisfied :

Let y 1 , y 2 , ..., y n+1 be a solution of DBKP . It is clear that x i = y i for i = 1, ..., n, x b = y n+1 and e g 1 = 0 is a solution of DLP gen . Conversely, let ((x i ) 1≤i≤n , x b , e g 1 ) be a solution of DLP gen , i.e. of the above set of constraints. According to constraint (*), n i=1 a i x i + a n+1 x b ≥ V -e g 1 and, according to constraint (**), e g 1 ≤ 1/2; then we have n i=1 a i x i + a n+1 x b ≥ V -1/2. Since the quantities n i=1 a i x i +a n+1 x b and V are integers, we finally obtain:


Julien Lavauzelle
email: lavauzelle@lix.polytechnique.fr

¸oise Levy-Dit-Vehel

New Proofs of Retrievability using Locally Decodable Codes

Keywords: retrievability, locally decodable codes, lifted codes, cryptographic protocols, data storage

Proofs of retrievability (PoR) are probabilistic protocols which ensure that a client can recover a file he previously stored on a server. Good PoRs aim at reaching an efficient tradeoff between communication complexity and storage overhead, and should be usable an unlimited number of times.

We present a new unbounded-use PoR construction based on a class of locally decodable codes, namely the lifted codes of Guo et. al.. Our protocols feature sublinear communication complexity and very low storage overhead. Moreover, the various parameters can be tuned so as to minimize the communication complexity (resp. the storage overhead) according to the setting of concern.

I. INTRODUCTION

Proofs of Retrievability (PoR) are cryptographic protocols aiming at addressing the following issue: given a client who stored a file on a server and erased its local copy, how can the client check if he is able to retrieve his file from the server, without any loss? PoRs thus involve two parties: a client (or verifier) and a server (or prover). Three main phases can be highlighted: first, the client encodes and/or encrypts his file and sends it to the server; then, client and server run challengeresponse rounds to check the file's extractibility; finally, if he is convinced, the client runs an extraction procedure to recover his file. Given a security bound on the extraction reliability, PoRs are targeting to minimize the communication cost between both parties, their computation complexity and their storage overhead.

PoRs were formally introduced by Juels and Kaliski (JK) [START_REF] Juels | PORs: Proofs of Retrievability for Large Files[END_REF] in 2007. Their construction is based on randomly chosen sentinel symbols the client secretly keeps before uploading his file. Then the verification process consists in checking the integrity of some sentinels. Since the sentinel positions are unknown to the server, if the client file is not retrievable, then w.h.p. some sentinels should be corrupted, and the client can notice it. Therefore, this scheme has low communication but remains bounded-use, as the number of verifications depends on how many sentinels the client actually possess.

Aiming at unboundedness, Shacham and Waters (SW) [START_REF] Shacham | Compact proofs of retrievability[END_REF] proposed to append to the file some authenticator symbols, and then to check random linear combinations of file symbols and authenticators. Verifying the file just consists in comparing parity symbols and the decryption of the encrypted ones. File extraction is then possible by the additional use of an efficient erasure code.

Part of BJO scheme can be seen as a double-layer encoding: an inner code is used to recover information symbols from paritycheck ones; an outer code helps with correcting erasures the inner code could have left. Dodis et. al. [START_REF] Dodis | Proofs of Retrievability via Hardness Amplification[END_REF] follow this idea to formalize the verification process as a request to a PoRcode, i.e. an error correcting code which models the space of possible answers to a challenge. Varying the underlying codes, they find good theoretical PoRs (using list-decodability bounds). However, most of their constructions are boundeduse, and the unbounded-use one does not really improve on SW scheme, despite lighter cryptographic assumptions. Following Dodis et. al., we carry on with using codes in PoRs. We investigate code constructions which are naturally suited for these protocols, namely locally decodable codes (LDC) and more specifically the high rate lifted codes introduced by Guo, Kopparty and Sudan in [START_REF] Guo | New Affine-Invariant Codes from Lifting[END_REF]. In the sequel, we introduce PoRs whose verification process is based on the local structure of these codes (section III), and whose security is analysed in section IV. We then give theoretical parameters and practical settings for our PoRs in section V. We concede that our practical results are not especially better than previous ones. Nevertheless, using LDCs in PoRs is new and opens a new line of research.

II. DEFINITIONS AND SECURITY MODEL

Proofs of retrievability involve two parties: a client who owns a (usually huge) file F , and a server on which F is stored. We essentially follow and lighten BJO [START_REF] Bowers | Proofs of Retrievability: Theory and Implementation[END_REF] definitions. A PoR system is composed of three main procedures:

• An initialisation phase. The client encodes his file F with an initialisation function Init(F ) = ( F , data). He keeps data (e.g. keys, authenticator symbols, etc.) for himself, then he sends F to the server and erases F .

• A verification phase. The client produces a challenge c with a randomized Chall function, then he sends c to the server. The latter creates a response r = Resp( F , c) and sends it back to the client. The client checks if r is correct by running Verif (c, r), which also access data, and outputs true if r is considered as a correct answer, false otherwise. • An extraction phase. If the client has been convinced by the verification phase, he can use his Extract algorithm to recover his whole file with high probability.

The client wants to be (almost) sure that he retrieves his file undamaged when using the Extract algorithm. But in practice, the server can always make the file unretrievable, for example by erasing it entirely. The point is that we want to model the fact that, if the server's answers to client's challenges make him look to own the file, then the client must be able to recover it entirely. Thus, an adversarial model (representing the server's behaviour) must be defined first.

Once again, we follow BJO [START_REF] Bowers | Proofs of Retrievability: Theory and Implementation[END_REF] context.

Definition II.1 ( -adversary). Let P be a PoR system and X be the space of challenges generated by Chall. An -adversary A for P is an algorithm such that, for all files F :

P x∈X Verif (x, A(x)) = false ≤ .
The key idea is the following: the client models the server as an -adversary, and his verification process is used to maintain an approximation of . Depending on this estimate, he can decide whether his file is retrievable or not. We thus define a way to measure PoRs' security:

Definition II.2 (PoR security). Let , ρ ∈ [0, 1]. A PoR system is said to be ( , ρ)-sound if, for all -adversaries A and for all files F , we have:

P Extract A = F ≥ ρ ,
where the probability is taken over the internal randomness of Extract A .

Parameter ρ represents the soundness of the scheme: we want ρ > 1 -2 -λ , where λ is a desired security parameter.

III. OUR POR CONSTRUCTION

In [START_REF] Dodis | Proofs of Retrievability via Hardness Amplification[END_REF], error correcting codes are used as a framework to model the verification process. We will show that some codes which naturally possess a local decodability property can have a greater role in PoRs.

A. Locally decodable codes and lifted codes

Error correcting codes are usually used to ensure reliable communication/storage in a noisy environment. To fully decode a message, best decoding algorithms reach almost-linear complexity in the length of the code. However, with the emergence of massive data storage, the question of correcting errors in sublinear time appeared. Hadamard code and Reed-Muller codes are known for a long time to provide this local decodability property, but their rate is stuck below 1/2. Major breakthroughs were rencently made by Kopparty et. al. [START_REF] Kopparty | High-rate codes with sublineartime decoding[END_REF], Guo et. al. [START_REF] Guo | New Affine-Invariant Codes from Lifting[END_REF], and Hemenway et. al. [START_REF] Hemenway | Local correctability of expander codes[END_REF]. They all give high-rate locally decodable codes (LDC) constructions, in a sense that, asymptotically in the code length, the code rate can be made close to 1. For more information on LDCs, see Yekhanin's survey [START_REF] Yekhanin | Locally decodable codes[END_REF]. It has to be pointed out that high rates can practically be reached with reasonable lengths, as we show in section V.

Here we especially focus on Guo et. al. construction, named lifted codes [START_REF] Guo | New Affine-Invariant Codes from Lifting[END_REF]. Let F q denote the finite field with q elements and S be any finite set. From now on we identify a function f : S → F q with its vector of evaluations (f (x)) x∈S . Thus, some error correcting codes can be represented as subsets of functions. For example, {f : F q → F q , deg f < k} is the q-ary Reed-Solomon code of length n = q and minimum distance

d = q -k + 1. Definition III.1 (Affine-invariant code). Let F Q be an ex- tension of F q . A code C ⊆ {f : F t Q → F q } is said to be affine-invariant if, for all affine permutations T : F t Q → F t Q , we have f • T ∈ C.
GKS [START_REF] Guo | New Affine-Invariant Codes from Lifting[END_REF] builds a LDC by picking an affine-invariant base code C ⊆ {f : F Q → F q }, and then lifting it to a subset L of {g : F m Q → F q }, such that each affine one-dimensional restriction of functions in L belongs to C. More formally: Definition III.2 (GKS affine lifting [START_REF] Guo | New Affine-Invariant Codes from Lifting[END_REF]). Let C ⊆ {f : F Q → F q } be an affine-invariant code. The m-th affine lift of C is the affine-invariant code:

Lift m (C) = {g : F m Q → F q | ∀ affine injections T : F Q → F m Q , g • T ∈ C} .
It is clear that these lifted codes have a local decodability property: to locally decode a codeword c ∈ L = Lift m (C), one just restricts c to an affine line and uses the decodability of the base code to recover some symbols.

B. Our construction

We now build our PoR based on lifted codes. Recall that in a PoR scheme, we want to know with high probability and sublinear communication whether a large file is still retrievable. Our main idea is to use the local structure of a lifted code to check the file integrity. Thus, the client encodes his file F as a lifted codeword F ∈ Lift m (C), and the verification procedure essentially consists in checking if symbols aligned on a random line form a codeword of C. However, as the client wants to prevent from malicious strategies, he will also need to encrypt his encoded file.

We thus assume the client has access to a family of symmetric stream ciphers φ κ , where κ belongs to some key space K chosen according to the context, with the following property: one can get the i-th key stream symbol in constant time. For instance, block ciphers in counter (CTR) mode are suitable for such stream ciphers (see example 9.9 in [START_REF] Schneier | Protocols, Algorithms, and Source Code in C[END_REF]).

We sum up the initialisation procedure in the algorithm of Fig. 1.

Input: a file F . Output: an encryption key κ and the ready-to-send file F .

1: Pick uniformly at random an encryption key κ ∈ K.

2: Encode F to F using an encoding procedure of a lifted code Lift m (C). 3: Encrypt each symbol c i of F with a block-cipher After his file is sent, the client can check its integrity: he simply needs to ask for a random 1-dimensional restriction of F , and check that the associated decrypted word lies in C. We summarize the verification phase in Fig. 2.

φ κ : ci = φ κ (c i , i). 4: return (data = κ, F = {c i } i ).
Client Server κ F Pick (P, v) ∈ F m q × (F m q ) × Compute points on line = {P + tv, t ∈ Fq} . Read symbols c = { F [x]} x∈ . Compute c = {φ -1 κ (cx, x)} x∈ . Return true if c ∈ F, false otherwise. (P, v) c Fig. 2. Verification protocol
Finally, if the verification protocol succeeds a sufficient number of times, the extraction algorithm of Fig. 3 ensures that with very high probability, the file will be recovered.

IV. SECURITY ANALYSIS

In this section, we prove that our PoR succeeds with very high probability. We denote by (C, m, γ)-LiftPoR our PoR construction which uses a base code C ⊆ {F Q → F q }, a lifting parameter m and voting parameter γ.

Informally, our PoR security is based on the fact that it is implausible for the server to create answers that pass the verification test but do not correspond to correct symbols. To represent these answers, we need to introduce what we call spurious codewords: Definition IV.1. Let P = (C, m, γ)-LiftPoR be our PoR scheme, and F the encrypted encoded file held by the server. A spurious codeword for (P, F ) related to a line ⊂ F m Q is a word w ∈ F Q Q such that w = F| and Verif ( , w) = true.

Input: γ ∈ [0, 1] a parameter for voting (γ = 0.5 means majority voting). Output: the file F .

---1. Initialisation ---1: for all x ∈ F m q do 2:

v x ← [0] t∈Fq 3: end for ---2. Client-server interaction ---4: for all lines ⊂ F m q do 5:

Run the verification procedure with some random line . Let c be the decrypted server's answer.

6:

if Verif data ( , c) = true then 7:

for all x ∈ do 8:

Increment v x [c x ] 9:
end for 10:

end if 11: end for ---3. Vote ---12: for all x ∈ F m q do 13:

if ∃s ∈ F q , v x [s] > γL then 14: We set L = q m -1 q-1 the number of lines passing through any point of F m q .

Tab[x] = argmax s v x [s]
Let us try to estimate the best chance for any adversary A to provide a spurious codeword on a challenge . Upon receiving , the adversary may know all the encrypted symbols y = (φ κ (c x , x)) x∈ asked by the client. If he wants to fool the client, the adversary needs to modify y on some positions in order to construct y = (φ κ (c x , x)) x∈ with c ∈ C \ {c}. Now, if we assume that the cipher φ κ looks to A like a random permutation, changing y x = φ κ (c x , x) into y x = φ κ (c x , x) is the same as uniformly picking a random c x ∈ F q \ {c x }.

Thus, with this assumption on φ κ , the only strategy for A is:

1) Pick a random subset I ⊆ ;

2) For all i ∈ I, change y x into random y x ∈ F q \ {y x }.

Hence, the only degree of freedom for A to optimize his strategy remains in his choice of the subset I.

In order to make the security analysis easier, we specify our PoR construction by instantiating the base code C as a Reed-Solomon code. Besides their MDS property, they are practically the best base codes w.r.t. the rate of their lifts.

Proposition IV.2. Let A be an adversary for a (C, m, γ)-LiftPoR where C is the q-ary Reed-Solomon code of minimum distance d. Then, for any line ⊂ F m q , the probability (over A's randomness) that A produces a spurious codeword on challenge is upper bounded by (q -1) -d+1 .

Proof. The proof is quite long and technical so we only sketch it here. For I ⊆ F q , let W I denote words in F q q with support I. Following our previous remarks on possible adversarial strategies, the aimed probability is upper bounded by: max

I⊆Fq |W I ∩ C| |W I | .
Now remark that W I ∩ C corresponds to maximum-weight codewords in the code C punctured on F q \ I. Studying the weight distribution of this punctured code (it is well-known when C is a Reed-Solomon code, see example 4.6 in [START_REF] Roth | Introduction to coding theory[END_REF]) and a careful computation provide the expected result.

Informally, Proposition IV.2 shows that it is unlikely that adversary A produces a spurious codeword. The voting phase in the extraction algorithm amplifies this hardness.

Proposition IV.3. Let L = q m -1 q-1 .
Then in the extraction algorithm of Fig. 3, we have for all x ∈ F m q , P[Tab[x] is wrong after step 3] ≤ 2 -L(γ(d-1) log(q-1)-1) .

Proof. Thanks to proposition IV.2, the aimed probability is bounded by the one defined by the binomial distribution with parameters n = L and p = (q -1) -d+1 . Adapting Lemma 4.13 of [START_REF] Roth | Introduction to coding theory[END_REF] on error distributions on codewords, we obtain the following bound:

P[Tab[x] is wrong after step 3] ≤ 2 -L•D(γ||p) ,
where D(γ||p) is the Kullback-Leibler divergence. In our case, a computation gives D(γ||p) > γ(d -1) log(q -1) -1 and achieves the proof.

Union bound implies that the probability there exists a position x ∈ F m q where Tab[x] is wrong after step 3 is bounded by 2 -L(γ(d-1) log(q-1)-1)+m log q . Thus, we just bounded the probability of having an error after the third step of the extraction algorithm. Let us now focus on potential erasures.

Lemma IV.4. Let γ ∈ [0, 1] and S ⊆ F m q with |S| = k. We set L = q m -1 q-1 . Let D be a set of lines of F m q such that ∀s ∈ S, |{ ∈ D, s ∈ }| ≥ (1 -γ)L. Then, |D| ≥ (1 -γ)kL - k(k -1) 2 .
Proof. By induction on k. Case k = 1 follows from the assumption on D. Assume now the proposition is true for k.

For S ⊆ F m q , let ∆ S = {lines ⊆ F m q , | ∩ S| ≥ (1 -γ)L}. Let S ⊆ F m q , |S| = k + 1
, and y ∈ S. We then have:

|∆ S | = |∆ S\{y} | + |∆ {y} | -|∆ S\{y} ∩ ∆ {y} | ≥ (1 -γ)kL - k(k -1) 2 + (1 -γ)L -k ≥ (1 -γ)(k + 1)L - (k + 1)k 2 .
Note that L is exactly the number of lines passing through a point of F m q . Hence, there remains an unfilled position (represented by an erasure ⊥) in s ∈ F m q if and only if at least (1 -γ)L challenges passing through s have been rejected by the Verif procedure.

Let S denote the set of remaining erasures and D the set of challenges (lines) rejected by Verif . We then have ∀s ∈ S, |{ ∈ D, s ∈ }| ≥ (1 -γ)L. Thus, by Lemma IV.4,

(1 -) q m-1 L # lines ⊂F m q ≥ |D| ≥ (1 -γ)kL - k(k -1) 2 , with k = |S|. It implies, if < (1 -γ)/2, that k < L(d -1)
. Now, let s ∈ S. Then, there exists a line s ⊂ F m q passing through s such that s has no more than d -1 erasures (otherwise, we refute k < L(d -1)), and we can decode the symbol on position s with C's decoding algorithm which can decode up to d -1 erasures.

Bringing together all our results, we have the following theorem:

Theorem IV.5. Let P be a (C, m, γ)-LiftPoR with C a qary Reed-Solomon code of minimum distance d. Then, for all < (1 -γ)/2, P is ( , 1 -2 -λ )-sound with λ = L((d -1)γ log(q -1) -1) -m log q.

Proof. The discussion after Proposition IV.3 implies that with the claimed probability, for any adversary A, there is no error in the temporary file Tab after step 3 of the extraction algorithm. Lemma IV.4 and following discussion show how to deal with potential erasures.

V. LIFTED-CODE-POR: THEORY AND PRACTICE

A. Theoretical parameters Fig. 4 gives the exact and asymptotic values of client and server storage overhead as well as the communication complexity of our PoRs.

We also emphasize that our PoR is unbounded-use, as a consequence of its "information-theoretical security" feature. Indeed, challenges are randomly chosen among lines of F m q and uniformly cover the file symbols; thus, the adversary cannot learn anything from past challenges. Fig. 5 asymptotically compares our scheme to Bowers et. al. [START_REF] Bowers | Proofs of Retrievability: Theory and Implementation[END_REF] and Shacham-Waters [START_REF] Shacham | Compact proofs of retrievability[END_REF] ones.

Server overhead |κ|N + R * |F | |F | 1-β R + R * |F | R * |F | Communication |κ| |F | β + |κ| R * |F | R 1/m Remark β ∈ ]0, 1[ m ≥ 2

B. Practical settings

Lifted codes are known to reach rates close to 1 when their length tends to infinity [START_REF] Guo | New Affine-Invariant Codes from Lifting[END_REF]. However, we need to highlight some good practical codes (not too long but with good rate and known encoding/decoding algorithms) in order to use them in real PoRs.

Recall that we instantiate our base codes with Reed-Solomon codes defined over the whole field F q . Thus, we have three parameters to deal with: field size q (which is also the base code length), minimum distance d < q and lifting parameter m ≥ 2. Fig. 6 presents some efficient choices for these parameters. On the one hand, the top grey row shows a very storage-efficient setting with only 3% server storage overhead; on the other hand, the bottom grey row highlights a good lowcommunication PoR with challenges of size less than 1/10000 fraction of the file size.

VI. CONCLUSION

To sum up, we have defined new unbounded-use Proofs of Retrievability with sublinear communication complexity, constant client storage and low server storage overhead. A first implementation also provides a proof of concept of our PoR; more efficient implementations of lifted codes underlying algorithms could then lead to competitive results. Besides, since the idea to use locally decodable codes in PoRs is quite generic, it would be interesting to design PoRs based on other such codes.

VII. ACKNOWLEDGEMENTS

The authors would like to thank Daniel Augot for his valuable and helpful comments. Fig. 6. Some practical settings for our PoR construction.

Fig. 1 .

 1 Fig. 1. Initialisation algorithm. Code C, integer m and key space K are system parameters.

4 .

 4 Decoding ---19: for all x ∈ F m q such that Tab[x] = ⊥ do 20: Find a line x passing through x such that #{y ∈ x , Tab[y] = ⊥} < d 21: Decode the word Tab | x with C's decoding algorithm, and update Tab on x with the codeword. 22: end for 23: return Tab.

Fig. 3 .

 3 Fig.3. Extraction algorithm. The parameters q, m and the code C are public. We set L = q m -1 q-1 the number of lines passing through any point of F m q .

Fig. 5 .

 5 Fig.5. Asymptotic comparison of some existing PoRs. log terms have been removed for clarity. Notations are identical to those in Fig.4, except that R * = 1 R -1 > 0 is the redundancy of the code.

  Fig. 4. Our PoR parameters. |F | denotes the file size in bits, q the field size, m ≥ 2 the lifting parameter, R the lifted code rate and |κ| the key size. We also have Rq m log q = |F |.

	C. storage overhead S. storage overhead C. → S. comm. S. → C. comm.	Exact value Asymptotics (|F | → ∞) O(1) |κ| 1 R -1 |F | O(|F |) 2m log q O(log |F |) q log q O(|F | 1/m )
	Construction	BJO [5]	SW [2]	Our work
	Unbounded-use No (N uses)	Yes	Yes
	Client storage	|κ|	|F | β + |κ|	|κ|