Keywords: Mathematics Subject Classification (2010). Primary 33C75, 33C20; Secondary 65B10 Alternating harmonic number, Infinite series, Integral transform, Harmonic series. H 2n, g n

introduced in our article. Our integration-based method, when applied in conjunction with the main lemma from [Ramanujan J 46:2, 2018], also provides a new way of evaluating series involving even-indexed harmonic numbers.

Series containing squared central binomial coefficients and alternating harmonic numbers

John M. Campbell

Abstract. We introduce an integration method for evaluating infinite series involving alternating harmonic numbers that generalizes many of the results in [Integral Transforms Spec Funct 28 :7, 2017]. Using this method, we provide new evaluations for series containing factors of the form 2n n 2 H 2n that cannot be evaluated using known generating functions involving harmonic-type numbers. A closed-form evaluation is given for the series and we describe why the problem of symbolically computing this particular series is especially difficult. The integration technique given in our article may be applied to evaluate natural generalizations and variants of the above series, such as the binomial-harmonic series

Introduction

The symbolic computation of infinite series containing harmonic numbers is a fascinating subject that traces back to Euler's seminal work on the evaluation of series of the form ∞ n=1

Hn n m . Of course, the study of the closedform evaluation of infinite series with summands containing factors given by harmonic-type numbers remains an active area of research. In our present article, we introduce a useful technique for evaluating infinite series involving alternating harmonic numbers.

New techniques for computing series containing products of harmoniclike numbers and central binomial coefficients are given in [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF][START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF][START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF][START_REF] Campbell | On the interplay between hypergeometric functions, complete elliptic integrals and Fourier-Legendre series expansions[END_REF][START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF][START_REF] Cantarini | On the interplay between hypergeometric series, Fourier-Legendre expansions and Euler sums[END_REF][START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF][START_REF] Chu | Infinite series with harmonic numbers and central binomial coefficients[END_REF][START_REF] Guillera | More hypergeometric identities related to Ramanujan-type series[END_REF][START_REF] Liu | Generalized harmonic number summation formulae via hypergeometric series and digamma functions[END_REF][START_REF] Nicholson | Quadratic Transformations of Hypergeometric Function and Series with Harmonic Numbers[END_REF][START_REF] Wan | Random Walks, Elliptic Integrals and Related Constants[END_REF]. Letting G denote Catalan's constant, our article is inspired in part by the recent formula

4G -12 ln 2 + 6 π = ∞ n=1 2n n 2 16 n (2n -1) 2 H 2n (1.1)
that was proved in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] following two separate techniques, through an application of the operator T ln,arcsin introduced in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] and also through the use of a generating function given by Chen in 2016 in [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], in which a variety of generating functions involving harmonic-like numbers and binomial-type expressions are given. The interesting formula given in (1.1) was subsequently noted in [START_REF] Sofo | General order Euler sums with multiple argument[END_REF], in which the evaluation of series of the form

letting H m denote the alternating harmonic number of order m ∈ N 0 , with

H 2n = H 2n -H n .
The integral transform T ln,arcsin introduced in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] in 2017 is used to prove new evaluations for series involving factors of the form indicated in (1.2), such as the series

n∈N 2n n 4 n (n + 1)(2n -1)(2n -3) 2 H 2n
introduced in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF]. In our present article, we introduce an integration-based method that may be applied to summations of the form n∈N 2n n for sequences (g n : n ∈ N), with this new technique generalizing and building upon a number of the results from [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF]. Since the integration method from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] is used to compute series of the form

n∈N 2n n 2 H n • g n ,
we thus obtain a new approach toward the problem of evaluating series involving even-indexed harmonic numbers.

In [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF], it is noted that there appears to be an interesting connection between series involving It is noted in [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF] that the integration method from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] cannot be applied directly to evaluate the series

∞ n=1 1 32 n 2n n 2 H n 2n -3 , (1.5)
in the sense that it is unclear how the required integral expression given by this method may be evaluated. If the Mathematica function Simplify is applied to this expression, we obtain a complicated expression involving numerous parameter derivatives of 3 F2 (-1) functions. A re-indexing argument is used in [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF] to evaluate (1.5), and this is used, together with the expression for (1.5) in terms of parametric derivatives on 3 F2 (-1) mappings, to show that the alternating series

∞ n=1 - 1 16 n 2n n 2 (2n + 3)H 2n (2n -1)(n + 1) (1.6)
is equal to

- 2 3 + 5(π -4 ln(2))Γ 2 1 4 12 √ 2π 3/2 - √ 2π(-20 + 3π + 12 ln(2)) 3Γ 2 1 4 ,
which motivates the study of summations containing (1.4) as a factor. We introduce an interesting integration technique in our present article to evaluate series of this form, offering something of a partial solution to an open problem given in [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF]. As noted in [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF], it is surprising that the fundamental lemma from [3] can be used to evaluate (1.6), since this lemma is specifically "designed" for series containing H n for n ∈ N, as opposed to series involving evenindexed alternating harmonic numbers. The evaluation of generalizations and variations of (1.6) is left as an open problem in [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF], and serves as a basis for our present article. This article is also inspired by the formula

6 ln(2) -2 π = ∞ n=1 2n n 2 H 2n 16 n (2n -1) (1.7)
which was recently proved in [START_REF] Campbell | On the interplay between hypergeometric functions, complete elliptic integrals and Fourier-Legendre series expansions[END_REF] through the use of Fourier-Legendre theory.

In [START_REF] Campbell | On the interplay between hypergeometric functions, complete elliptic integrals and Fourier-Legendre series expansions[END_REF], it is noted that from the FL expansion

- 3 2 + ∞ n=1 (-1) n 2 1 n - 1 n + 1 - 1 n + 1 n + 1 P n (2x -1) (1.8)
for the expression ln(1 -√ x), together with the formula

1 0 x n ln(1 - √ x) dx = - H 2n+2 n + 1 ,
we may evaluate the series in (1.7), by using the formula

1 0 E √ x ln(1 - √ x) dx = - π 2 ∞ n=0 2n n 2 H 2n+2 16 n (n + 1)(1 -2n)
given in [START_REF] Campbell | On the interplay between hypergeometric functions, complete elliptic integrals and Fourier-Legendre series expansions[END_REF] and applying the FL series provided in (1.8). It is not clear how one would go about proving new FL expansion formulas in order to evaluate more general classes of harmonic-binomial sums of the form noted in (1.3). This inspires us to construct a simple integration technique that easily allows us to evaluate series of this form, by analogy with the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], which is given below. From something of a heuristic perspective, this lemma shows us that if we apply the binomial expansion for √ 1 -x2 and then multiply the Maclaurin series for this elementary function by the inverse of this mapping, we can manipulate the resultant equation so as to be able to express a given harmonic-binomial series in terms of a relatively manageable integral and the underlying non-harmonic summation of the initial series.

Lemma 1.1. If the series ∞ n=0 (-1) n 1 2 n f (n) x 2n ln(1 -x 2 ) √ 1 -x 2 is integrable on [0, 1] for a sequence (f n ) n∈N0 , and if we let g n = fn 16 n (2n-1) , then ∞ n=0 g(n) 2n n 2 H n equals the product of 2 π with the sum of 1 0 ∞ n=0 (-1) n 1 2 n f (n) x 2n ln(1 -x 2 ) √ 1 -x 2 dx and π ln(2) ∞ n=0 ( 2n n )
As discussed in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], the above lemma provides us with a useful way of producing new Ramanujan-like series for 1 π that cannot feasibly be proved using the generating functions from papers such as [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF].

Main results

The following lemma may be regarded as a direct analogue of the main technique from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF].

Lemma 2.1. If we let (f n : n ∈ N 0 ) denote a sequence such that ∞ n=0 (-1) n x 2n 1 -x 2 -1 2 n f n ln (x) (2.1)
is integrable on [0, 1], writing g n in place of fn n+1 , we have that the series

n∈N0 1 16 n H 2n 2n n 2 g n is equal to 4 π times the definite integral with respect to x of (2.1) over [0, 1] plus 1 2 ∞ n=0 1 16 n 2n n 2 (2 ln(2)(n + 1) + 1) (n + 1) 2 f n .
Proof. This follows immediately from the identity whereby

1 0 x 2n 1 -x 2 ln(x) dx evaluates as √ π H n-1 2 -H n+1 Γ n + 1 2 8Γ(n + 2) for (n) > -1 2 .

2.1.

A seemingly recalcitrant series containing alternating harmonic numbers As described below, the problem of evaluating the infinite series given in (2.2) below is very difficult, and gives rise to challenging new problems related to the subject matter from [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF][START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF][START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF]. The problem of determining the symbolic value of the infinite series in Question 2.2 below is inspired in part by the evaluation (1.1) from [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] and the equation (1.7) introduced in [START_REF] Campbell | On the interplay between hypergeometric functions, complete elliptic integrals and Fourier-Legendre series expansions[END_REF]. If we apply partial fraction decomposition in the summand in (1.6), we obtain the series given below.

Question 2.2. How can the series

∞ n=0 - 1 16 n 2n n 2 H 2n n + 1 (2.2)
be evaluated in closed form?

Since there are known algorithms for computing series of the forms noted below for z ∈ Z >0 [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF][START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], it seems natural to consider the evaluation of variants of the following series, further motivating the problem given in Question 2.2.

∞ n=1 2n n 2 H n 16 n (n + z) ∞ n=1 2n n 2 H n 16 n (2n -2z + 1) ∞ n=1 2n n 2 H n 32 n (n + z) ∞ n=1 2n n 2 H n 32 n (2n -2z + 1) ∞ n=1 2n n 2 H n 16 n (n + z) 2 ∞ n=1 2n n 2 H n 16 n (2n -2z + 1) 2 ∞ n=1 2n n 2 H n 32 n (n + z) 2 ∞ n=1 2n n 2 H n 32 n (2n -2z + 1) 2
Mathematica 11 is not able to evaluate the equivalent formulation

∞ n=0 - 1 16 n 2n n 2 (H 2n -H n ) n + 1 (2.3)
of the series in Question 2.2 at all. We remark that it took Mathematica 11 about two straight hours of running time before providing any output after inputting the expression in (2.3), and only outputted the same symbolic form in (2.3), without providing any kind of simplification or evaluation. This illustrates that the problem of evaluating the series in Question 2.2 is highly nontrivial. Unsurprisingly, Maple 18 cannot evaluate the series in (2.3) at all. Even if we directly make use of the definition of the term alternating harmonic number in our attempts to evaluate (2.2) using computer algebra systems, we find that Mathematica 11 cannot evaluate

∞ n=0 - 1 16 n 2n n 2 2n i=1 (-1) i+1 i n + 1 , (2.4)
in that Mathematica 11 is only able to to compute (2.4) as

∞ n=0 - 1 16 n 2n n 2 (ln(2) -Φ(-1, 1, 2n + 1)) n + 1 ,
with this evaluation following trivially from the definition of the Lerch transcendent. Maple 18 is only able to provide the redundant evaluation

∞ n=0 - 1 16 n 2n n 2 1 2 Ψ n + 1 2 -1 2 Ψ(n + 1) + ln(2) n + 1
of (2.4), letting Ψ denote the digamma function. Now that we have discussed how current CAS software cannot evaluate (2.2) directly, we proceed to elaborate on the following claim concerning the difficulty of solving Question 2.2: It is not feasible to apply the generating functions from [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF][START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF] to evaluate the series in (2.2).

In [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], Chen noted that the power series

∞ n=1 2n n H 2n n + 1 x n (2.5)
may be evaluated as

1 2x 1 - √ 1 -4x + 1 + √ 1 -4x ln 1 + √ 1 -4x 2 , (2.6)
and we remark that Mathematica is able to compute (2.5) directly. However, there is no known symbolic evaluation for

∞ n=1 2n n 2 H 2n n + 1 x n .
In terms of power series that are known to have closed-form evaluations, it appears that (2.5) is by far the "closest" or most similar to (2.2), modulo scalar multiples of the variables involved. There are a number of binomialharmonic generating functions that are much less similar to (2.2) computed in [START_REF] Boyadzhiev | Series with central binomial coefficients, Catalan numbers, and harmonic numbers[END_REF][START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], including the generating functions noted below.

G 2n n H n ; x G 2n n n + 1 H n ; x G 2n n H 2n ; x G 2n n n (H 2n-1 -H n ) ; x
From the evaluation of the ordinary generating function

G 2n n H 2n n + 1 ; x given in (2.6), we find that ∞ n=1 - 1 16 n 2n n (H 2n -H n ) n + 1 x n is equal to 4 √ x + 4 - √ x + 4 + 2 ln 1 4 √ x + 4 + 2 -2 x (2.7)
From the above discussion, we see that

G - 1 16 n 2n n (H 2n -H n ) n +
π 3/2 - 4 π π 2 0 1 + cos 2 (t) ln 1 + 1 + cos 2 (t) 2 sec 2 (t) dt
which cannot be evaluated using Mathematica 11 or Maple 18. It is unclear as to how the above integral may be evaluated, especially since the two terms resulting from the expansion of the integrand in

π 2 0 1 + cos 2 (t) ln 1 + 1 + cos 2 (t) -ln(2) sec 2 (t) dt
cannot be separately integrated on 0, π 2 . So, we see that the problem of evaluating the series in Question 2.2 is equivalent to determining a closedfrom for the difficult integral

2 1 √ u ln 1+ √ u 2 √ 2 -u(u -1) 3/2 du, (2.8)
which cannot be evaluated by Mathematica 11 or Maple 18. The above integral cannot be evaluated by rewriting the expression ln 1+ √ u 2 as ln(1 + √ u) -ln(2) and expanding the resultant integrand, since the integral

2 1 √ u √
2 -u(u -1) 3/2 du does not converge. It is not at all obvious as to how to evaluate the definite integral displayed in (2.8). It seems as if it would not be feasible to expand a given combination of factors in the integrand in (2.8) as a Taylor series and integrate term-by-term. For example, we have that

1 √ 2 -u(u -1) 3/2 = ∞ n=0 (-1) n 2 n+2 u - 3 2 n 1 2 (n + ((n -1) mod 2)) 1 2
, but it seems that this expansion cannot be directly used to evaluate (2.8), especially since state-of-the-art computer algebra systems cannot evaluate

- 2 π 1 0 1-x x ln 1 + √ 1 + x √ 1 + x dx - 4 √ 2π ln(2) Γ 1 4 2 + Γ 1 4 2 ln(2) √ 2π 3/2 ,
with state-of-the-art computer algebra systems being unable to evaluate the above integral. The Maclaurin series for ln(1 + y) cannot be used to evaluate the above integral, as may be verified. Also, the Maclaurin series for 1-x 1+x cannot be used to compute this integral. Even if we try to make use of non-standard integral formulas for central binomial coefficients or Catalan-type sequences as a way of approaching Question 2.2, we obtain recalcitrant integrals. For example, if we try to apply the recently-discovered Catalan number formula

C n = 1 π ∞ 0 √ t t + 1 4 n+2 dt
given in [START_REF] Qi | An integral representation, complete monotonicity, and inequalities of the Catalan numbers[END_REF], we find that (2.2) is equal to

∞ 0 - 8 √ 2 √ x ln 1 2 1 4x+1 + 1 + 1 π √ 2x + 1(4x + 1) 3/2 dx,
with Mathematica 11 and Maple 18 unable to compute the above expression.

If we try to make use of the standard integral formula for the factor 1 n+1 in the summand in (2.2), we find that this would require the evaluation of

∞ n=0 - 1 16 n 2n n 2 (H 2n -H n ) x n , (2.10)
with Mathematica 11 and Maple 18 being unable to evaluate (2.10) at all. From the above discussion, it should be clear that we cannot use known generating functions for harmonic/binomial-type sequence to evaluate the seemingly recalcitrant series in (2.2). Moreover, it seems that it is not feasible to use known integral formulas for alternating harmonic numbers to evaluate the infinite series in (2.2). For example, if we try to make use of the known identity

1 0 ln(x + 1)x 2n dx = -1 + 2 ln 2 + 4n ln 2 (2n + 1) 2 - H 2n 2n + 1
noted in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF], we find that the problem of symbolically computing the aforementioned series is equivalent to the problem of evaluating the difficult integral

1 0 ln(x + 1) (E (ix) -K (ix)) x 2 dx,
which cannot be evaluated by Mathematica 11 or Maple 18. Question 2.2 is also interesting because it is unclear as to how the corresponding series

∞ n=0 - 1 16 n 2n n 2 H 2n n + 1 (2.11)
with non-alternating harmonic numbers of even index could be computed symbolically, as a corollary of the following theorem, using the main lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF].

Theorem 2.3. The series ∞ n=0 - 1 16 n 2n n 2 H 2n n + 1 (2.12) is equal to 2 + (4 ln(2) -π)Γ 2 1 4 4 √ 2π 3/2 - √ 2π(4 + π + 4 ln(2)) Γ 2 1 4 .
Proof. This follws immediately from Lemma 2.1 by letting f n = (-1) n .

It is truly remarkable that the integration technique given by Lemma 2.1 provides such a simple way of evaluating the infinite series in (2.2), since the problem of determining the symbolic value of the series in Theorem 2.2 may be regarded as being very difficult, as described above.

Generalizations and variants

Through the use of Lemma 2.1, we also obtain the following extensions of the above theorem.

∞ n=0 - 1 16 n 2n n 2 H 2n n + 2 = - 8 9 + √ 2π(4 + π + 4 ln(2)) 3Γ 2 1 4 + (4(4 + 3 ln(2)) -3π)Γ 2 1 4 108 √ 2π 3/2 ∞ n=0 - 1 16 n 2n n 2 H 2n n + 3 = Series containing 2n n 2 H 2n 11 128 225 - (32 + 39π -156 ln(2))Γ 2 1 4 540 √ 2π 3/2 - √ 2π(316 + 95π + 380 ln(2)) 375Γ 1 4 2 ∞ n=0 - 1 16 n 2n n 2 H 2n n + 4 = √ 2π(516 + 145π + 580 ln(2)) 875Γ 2 1 4 + (1360 -273π + 1092 ln(2))Γ 2 1 4 20580 √ 2π 3/2 - 512 1225
Applying partial fraction decomposition to the rational component of the series in (1.6), we obtain the expression

1 3 ∞ n=1 1 16 n (-1) n+1 2n n 2 H 2n n + 1 + 8 3 ∞ n=1 1 16 n (-1) n 2n n 2 H 2n 2n -1
and in view of Theorem 2.3 this strongly motivates the evaluation of the series given in the following theorem, since much of our present article is inspired by the open problem given in [START_REF] Campbell | New series involving harmonic numbers and squared central binomial coefficients[END_REF].

Theorem 2.4. The equality

∞ n=0 -1 16 n 2n n 2 H 2n 2n -1 = (π -4 ln(2))Γ 2 1 4 8 √ 2π 3/2 - π 2 (π + 4 ln(2) -4) Γ 2 1 4 holds.
Proof. This follows immediately from Lemma 2.1 by letting f n = (-1) n (n+1) 2n-1

.

The fundamental lemma introduced in our present article my also be used to evaluate

∞ n=0 - 1 16 n 2n n 2 H 2n 2n -2z -1 for z ∈ Z >0 , yielding the following evaluations. ∞ n=0 - 1 16 n 2n n 2 H 2n 2n -3 = √ 2π(-4 + π + 4 ln(2)) 6Γ 2 1 4 + Γ 2 1 4 (16 + 3π -12 ln(2)) 216 √ 2π 3/2 ∞ n=0 - 1 16 n 2n n 2 H 2n 2n -5 (39π -4(8 + 39 ln(2)))Γ 2 1 4 1080 √ 2π 3/2 - π 2 (-316 + 95π + 380 ln(2)) 375Γ 2 1 4 ∞ n=0 - 1 16 n H 2n 2n n 2 2n -7 = (1360 + 273π -1092 ln(2))Γ 2 1 4 41160 √ 2π 3/2 + π 2 (-516 + 145π + 580 ln(2)) 875Γ 2 1 4
2.3. On the evaluation of series involving even-indexed harmonic numbers and squared central binomial coefficients As a way of further demonstrating the utility of Lemma 2.1, we offer a simplified proof of the following result that had been introduced in [START_REF] Campbell | On the interplay between hypergeometric functions, complete elliptic integrals and Fourier-Legendre series expansions[END_REF].

Theorem 2.5.

∞ n=1 1 16 n ( 2n n ) 2 H2n
2n-1

= 6 ln(2)-2 π [4].
Proof. Letting f n = n+1 2n-1 , by Lemma 2.1 we have that

∞ n=1 1 16 n 2n n 2 H 2n 2n -1 = 2 -2 ln(2) π
and from the evaluation of

∞ n=1 1 16 n 2n n 2 H n 2n -1
given in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] we obtain the desired result.

Remark 2.6. It is interesting to note that Mathematica is able to evaluate the series in Theorem 2.5 as

1 2 2 F 1 (0,0,1,0) -1 2 , 1 2 1 1 + 2 ln(2) π and since ∞ n=1 2n n 2 H n 16 n (2n -1) = 2 F 1 (0,0,1,0) -1 2 , 1 2 1 1 ,
as noted in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] we thus obtain another proof of the above theorem.

Mimicking the proof of Theorem 2. The series in (1.1) from [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] is a natural extension of Chen's series for 1 π from [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], and it is natural to explore new methods of deriving the formula in (1.1), in the hope of arriving at similar results. is equal to 6+4G-12 ln(2) π [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF].

Proof. This follows from Lemma 2.1 by letting f n = n+1 (2n-1) 2 , using the symbolic form for

∞ n=1 1 16 n 2n n 2 H n (2n -1) 2
given in [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF].

Remark 2.8. By rewriting the expression H 2n in the summand in (2.13)

H 2n = 1 0 1 -x 2n 1 -x dx
and reversing the order of summation and integration, it can be shown that the above theorem also follows from Theorem 2.5.

Using the technique given in the proof of Theorem 2.7, we obtain the following results. Again through an application of Lemma 2.1 together with the fundamental lemma from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] as in the proof of Theorem 2.7, we obtain the following evaluations, none of which can be symbolically computed by Mathematica The integral transform introduced in our present article, together with the evaluations presented above, nicely generalize some of the concepts presented in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] introduced in [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] follows from some of the main results put forth in the present section of our article together with some of the main results from [START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF], by rewriting the above summand as 2n n 2

16 n (n + 1) 2 (2n -3) 2 (2n -1) 2 H 2n and applying partial fraction decomposition to the rational component 1 (n + 1) 2 (2n -3) 2 (2n -1) 2 in the above summand.

Conclusion

The main lemma introduced in our article can also be applied to series that do not involved squared central binomial coefficients. For example, a direct application of Lemma 2.1 shows that , and this can also be determined using Chen's generating function for

2n n n + 1 H 2n : n ∈ N 0
given in [START_REF] Chen | Interesting series associated with central binomial coefficients, Catalan numbers and harmonic numbers[END_REF], together with a Wallis-type integral. The exploration of further applications of Lemma 2.1 seems like a worthwhile area to pursue.

2 H2 H 2n and 2n n 2 H

 22 From the elegant evaluation provided in (1.1), together with the technique introduced in[START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] for evaluating series containing 2n n n as a factor for indices n ∈ N, we are inspired to explore the symbolic computation of series with summands involving factors such as 2n n 2n ,

  [START_REF] Campbell | An integral transform related to series involving alternating harmonic numbers[END_REF] given above, by letting f n = n+1 2n-(2z-1) for natural numbers z, we obtain the evaluations indicated below.

16G

  As discussed below, it appears that this would not be feasible, which is interesting since many similar series can be computed directly through the use of Wallis-type integrals.A standard version of Wallis's integral formula is given by the identity

	whereby				2n n	=	0	π 2	2(2 cos(t)) 2n π	dt
	for n ∈ N 0 . So, we find that the problem of evaluating the series in Question
	2.2 is equivalent to the problem of evaluating
	2 -	4	√	4 2Γ 3	2

1 ;

x is the "closest" possible generating function relative to (2.2) that can be computed in terms of elementary functions. So, given the relatively simple closed-form evaluation given in (2.7), it would seem as if the summation in Question 2.2 could be evaluated by making use of an integral formula for central binomial coefficients.

  It is interesting to note that the main lemma from[START_REF] Campbell | Ramanujan-like series for 1 π involving harmonic numbers[END_REF] cannot be applied directly to evaluate series of the form Based on the evaluation provided in the above theorem, through the use of iterative re-indexing, we obtain the following new evaluations.

	Series containing 2n n	2 H 2n	15
	However, if we make use of the canonical integral formula
	∞ n=0 2i i i + 1 for the Catalan sequence, we observe that the problem of finding a sym-1 16 n 2n n 2 H 2n (2.14) = 4 0 x i 4-x x dx 2π (n + z) 2
	bolc evaluation for the hypergeometric expression in (2.16) is equivalent to for z ∈ Z. A general method for evaluating series of this form is given in [2], evaluating the seemingly recalcitrant integral through the use of a clever re-indexing argument. Theorem 2.9. The infinite series ∞ n=1 2n n 4 π 4 0 x -4 + 2 √ 4 -x 1 -2 ln(2) + ln 2 + √ 4 -x dx x 3/2 2 H 2n which Mathematica, interestingly, is able to directly evaluate, as (2.15) 16 n (n + 1) 2 is equal to 16G+24-48 ln(2) π 16(-2G + 3 + π(ln(2) -1)) , π + 4 -8 ln(2). thus providing us with the desired result.
	Proof. Through a direct application of Lemma 2.1 we find that the series
	∞	n (H 2n -H n ) 2n	2
	n=1 -27 8 1 2 , 1 2 , 1, 1 16 n (n + 1) 2 32 ln(2) -9 2, 2, 2 1 -4 + H 2n 64G + 144 -176 ln(2) 9π is equal to 1 2 4 F 3 ∞ 2 2n = n=1 n 4 n (n + 2)	16 ln(2) π	.
	So, from the evaluation ∞ n=1 2n n 76800G + 200528 -203280 ln(2) 16875π 2 H n 16 n (n + 1) 2 = 16 + ∞ 2n 2 = n=1 n 4 n (n + 3) H 2n	-32G -64 ln(2) 2752 3375 512 ln(2) -225 π -16 ln(2)
	introduced in [2], we find that the problem of evaluating binomial-harmonic series in the above theorem reduces to the evaluation of a 4 F 3 (1) series with 45158400G + 128732144 -116933040 ln(2) 13505625π -115456 128625 2048 ln(2) -1225 parameters in 1 ∞ 2n 2 Z. However, the expression 4 F 3 1 2 , 1 2 , 1, 1 2, 2, 2 1 n = n=1 4 n (n + 4)
	cannot be evaluated by Mathematica, not even through the use of	11
	or Maple 18.
	92 + 24π -180 ln(2) 27π FunctionExpand = ∞ 1 16 n 2n n or through the use of the equivalent hypergeometric form in 2 H 2n n + 2 n=1 9628 + 1920π -16020 ln(2) 3375π = ∞ n=1 1 16 n 2n n ∞ 2n 2 2 H 2n n + 3 n=0 n 16 n (n + 1) 3 ,	(2.16)
	105788 + 17920π -160020 ln(2) 42875π which Maple is only able to evaluate in terms of the Barnes G-function as = ∞ n=1 1 16 n 2n 2 H 2n n n + 4 1 π • G 4,1 4,4 -1 1,2,2,2 1,1, 1 2 2 , 1 .

2 H 2n

  . Letting C n denote the Catalan number of order n, we see that

	the formula
	1500G -1946 + 2048 ln 2 3375π	-	244 3375	+	8 ln 2 225	=	∞ n=1	C 2 n 16 n (2n -3) 2 (2n -1) 2 H 2n

f (n)

H 2n


Wim Van Ackooij
email: wim.van-ackooij@edf.fr

Jérôme De Boeck
email: jdeboeck@ulb.ac.be

Boris Detienne

Stefania Pan

Michael Poss
email: michael.poss@lirmm.fr

Optimizing power generation in the presence of micro-grids

Keywords: integer programming, power generation, bilevel optimization, stochastic optimization, single leader multiple follower games 2010 MSC: 90C10, 90C27, Micro-grids: smaller nearly isolated sub-grids

come   L'archive ouverte pluridisciplinaire

Introduction

In the upcoming future the energy landscape will significantly change from the current picture by incorporating more and more decentralized elements. Of particular importance is the advent of so-called micro-grids (MG). These are subparts of the system with an advanced energy management system interacting with programmable elements in the grid including good monitoring and control functions, a pervasive communication system and specific items such as smart meters, programmable loads, switchable storage systems and a variety of controllable energy sources including solar, wind and wave generators. Some of the major changes introduced by smart grids are the following: load/offer mismatch occurs. Most importantly these sub-grids can be managed to follow a local economical target (which may be different and contrary to a system-wide interest).

• Partial storage, perhaps through electrical vehicles or powerful batteries. These storage devices can partially mitigate the intermittency of local decentralized production such as wind / solar generation.

• Demand management tools: use advanced information technology to pilot electricity use. For instance, shut down electrical heating, reprogram hot water tank recharging etc...

In this new setting it becomes of great interest to examine the interaction with more traditional elements composing the power system. For instance what will be the new role for large centralized generation assets such as nuclear, thermal or hydro generation? It also becomes of interest to examine how classical energy management questions, such as unit-commitment (e.g., [START_REF] Tahanan | Large-scale unit commitment under uncertainty: a literature survey[END_REF][START_REF] Van Ackooij | Large-scale unit commitment under uncertainty: an updated literature survey[END_REF]), scheduling maintenance of large power plants or cascaded reservoir management (e.g., [START_REF] Taktak | An overview on mathematical programming approaches for the deterministic unit commitment problem in hydro valleys[END_REF][START_REF] De Matos | Assessing policy quality in a multistage stochastic program for long-term hydrothermal scheduling[END_REF]) should evolve to account for this new context. Needless to say, these classical management questions are already difficult on their own without considering a potential interaction with micro-grids. On top of this comes also the need to consider and account for uncertainty which are of great importance for obtaining meaningful management solutions.

The interaction between such micro-grids and centralized assets may be mathematically cast into the setting of a "game". We will make the specific choice of considering hierarchical games wherein a (single) large centralized operator (the leader or upper level) interacts with one or several micro-grids (the follower(s) or lower level) in a specific way. Such problems also go by the name of single leader multifollower or in the case of a single follower by the naming of Stackelberg game or principal agent problems. We follow the terminology of [START_REF] Dempe | Bilevel Programming Problems: Theory, Algorithms and Applications to Energy Networks[END_REF] speaking of bilevel optimization problems, i.e., an optimization problem containing constraints that certain variables belong to the optimal set of other parametric optimization problems. Generally speaking these problems fall into the class of mathematical problems with equilibrium constraints, where [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity problems: Volume I[END_REF][START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity problems[END_REF] are key references. The situation with multiple leaders is significantly harder to analyse and falls into the class of Equilibrium problems with Equilibrium constraints (EPECs). We refer to, e.g., [START_REF] Henrion | Analysis of m-stationary points to an epec modeling oligopolistic competition in an electricity spot market[END_REF][START_REF] Aussel | Nash equilibrium in a pay-as-bid electricity market: Part 1-existence and characterization[END_REF][START_REF] Aussel | Nash equilibrium in a pay-as-bid electricity market: Part 2 -best response of a producer[END_REF][START_REF] Heymann | Mechanism design and auctions for electricity network[END_REF] for some applications in energy and to e.g., [START_REF] Pang | Quasi-variational inequalities, generalized nash equilibria, and multi-leader-follower games[END_REF][START_REF] Pang | Quasi-variational inequalities, generalized nash equilibria, and multi-leader-follower games[END_REF] and references therein for structural studies. Under some appropriate structure, a bilevel optimization problem involving only continuous variables can be cast as a DC (difference of convex) -problem. We refer to [START_REF] Pang | Computing b-stationary points of nonsmooth DC programs[END_REF][START_REF] Van Ackooij | DC programming techniques with inexact subproblems' solution for general DC programs[END_REF] for more on such methodology.

In [START_REF] Hobbs | A nonlinear bilevel model for analysis of electric utility demand-side planning issues[END_REF] the authors claim to investigate bilevel programming in the electric utility industry for the first time. The lower level (or follower), representing customers deals with maximizing benefit while investing in energy conservation. From a technical viewpoint, the authors substitute the Karush-Kuhn-Tucker (KKT) conditions for the lower level problem and process the subsequent mathematical program with equilibrium constraints (MPEC), by employing an exhaustive search on the complementarity conditions (one classic problem is solved for each state of the complementarity conditions). This procedure means that for an n dimensional complementarity constraint, 2 n problems need to be solved. Substituting the lower level KKT conditions in order to obtain a one-level more classic problem is commonplace in the energy literature (e.g., [START_REF] Kardakos | Optimal offering strategy of a virtual power plant: A stochastic bi-level approach[END_REF][START_REF] Cervilla | Bi-level optimization of electricity tariffs and pv distributed generation investments[END_REF][START_REF] Asimakopoulou | Leader-follower strategies for energy management of multimicrogrids[END_REF] just to name a few). However it should be noted that, following [START_REF] Dempe | Is bilevel programming a special case of a mathematical program with complementarity constraints?[END_REF], this does not necessarily lead to an equivalent formulation of the original problem. In other words, the solution resulting from the KKT reformulation need not be a solution to the original bilevel program. This phenomenon can occur even in the simplest linear setting. We also refer to the recent [START_REF] Adam | On M-stationarity conditions in MPECs and the associated qualification conditions[END_REF] illustrating a similar phenomenon. The essential observation is that the original problem is augmented with a given set of Lagrange multipliers. Consequently any local solution to the original problem must be a local solution to the augmented problem for all such Lagrange multipliers. However, whenever the augmented problem is solved, one ends up with just one specific Lagrange multiplier. In principle it should thus be verified if the thus obtained solution remains a local solution for all other Lagrange multipliers. Yet these are in general not readily available. For a simple 2 variable example of a situation wherein the KKT reformulation provides local solutions that are not local solutions of the original problem, see [START_REF] Dempe | Bilevel Programming Problems: Theory, Algorithms and Applications to Energy Networks[END_REF]Example 3.1]. The authors of [START_REF] Haurie | A two-player game model of power cogeneration in new england[END_REF] also consider a leader-follower framework in the electrical industry but rather in the objective of benchmarking a given investment strategy in cogeneration against a centrally planned strategy. In order to do so, they develop an iterative algorithm, which numerically seems to converge to a Betrand-Cournot equilibrium. The comparison with the centrally planned solution reveals inefficiencies of the given investment strategy.

Recently, demand side management (or demand response) was recognized to be an important aspect of the upcoming electrical system. As such, the authors of [START_REF] Fernández-Blanco | Incorporating price-responsive demand in energy scheduling based on consumer payment minimization[END_REF] claim to be the first to consider such demand response within a bilevel framework. The lower level (or follower) problem is an optimal power flow problem under a direct current assumption. The lower level power balance equation and its dual multiplier in turn help fixing the global price of energy in the upper level (or leader) problem (e.g., [START_REF] Hao | New models for integrated short-term forward electricity markets[END_REF][START_REF] Luh | Payment cost minimization with demand bids and partial capacity cost compensations for day-ahead electricity auctions[END_REF] suggest a similar model but without solving it). From a methodological viewpoint, bilinear terms are linearized using McCormick envelopes, the lower level KKT conditions are substituted for the lower level problem and the special model structure can be exploited through several substitutions.

In view of those claims our work is one of the first to consider the effect of demand side management on classical energy management questions. Moreover we will highlight that demand side management is responsible for ensuring that the optimistic and pessimistic solutions of the bilevel program are not identical. Although in our work we will focus on the interaction between a utility and several end-users grouped within one or several smart-grids, we will not consider the interaction between users as done in [START_REF] Mohsenian-Rad | Autonomous demand-side management based on game-theoretic energy consumption scheduling for the future smart grid[END_REF] where the focus is on reducing peak-to-average load ratios while accounting for potential games between users. As usual in these works, we too will adopt a market-like setting wherein a predetermined set of contracts is offered to the customers.

In this paper we thus contribute to this difficult question by considering a stylized interaction between a generation company (GenCo) owning a (large) set of centralized assets and a set of micro-grids. We also account for potential competition at the centralized level, but only in a simplified way. Indeed we will make the assumption that any competitors to the GenCo have a fixed predetermined interaction with the micro-grids. In this work the GenCo can offer contracts to the micro-grid that detail the price of buying/selling electricity to the network. From a mathematical viewpoint we suggest a bilevel stochastic mixed integer program (BMIP) which in principle can account for a variety of instantiations thus covering any of the above management questions. In our numerical experiments however, we will focus on unit-commitment. The first general methods for obtaining "optimistic" optimal solutions BMIPs can be traced back to [START_REF] Moore | The mixed integer linear bilevel programming problem[END_REF] and [START_REF] Bard | An algorithm for the discrete bilevel programming problem[END_REF], which could solve only small instances, with up to 10 general integer variables and 35 binary variables for the upper level problem. While more recent works have shown significant improvements (e.g., [START_REF] Fischetti | Intersection cuts for bilevel optimization[END_REF] based on intersection cuts), they are unable to handle the complex MILPs involved in unitcommitments or more complex power flow optimization problems. Our technical contribution therefore lies in suggesting a tractable reformulation of the bilevel problem that leads to a solution that is neither "optimistic" nor "pessimistic". This "heuristic" reformulation is tractable thanks to a preprocessing step that is not harder than solving the original management problem. We also compare the solution obtained from this reformulation with the "optimistic" solution that can only be obtained in reasonable times for small systems. Extensive numerical experiments confirm the interest of the here suggested reformulation. Essentially our reformulation takes inspiration, and to some degree is equivalent with, the value function approach originally designed by Outrata [START_REF] Outrata | On the numerical solution of a class of stackelberg games[END_REF]. An advantage of this approach over the usual substitution of the KKT conditions is that we do obtain a solution of the original problem, albeit, not necessarily the "optimistic" solution.

The outline of this work is as follows. The assumptions of our model are described in Section 2, where we present the ingredients of our model in a deterministic setting. Section 3 deals with the suggested reformulation technique for the bilevel problem. In Section 4, we show how our models can be adapted to handle uncertainty. We also show how our suggested reformulation can be adapted to this new setting. We assess our models and methodology in Section 5 on a case study built from realistic data. Some of the details of the model underlying the case study are provided in Appendix A. Finally, we conclude the paper in Section 6.

The problem

Description

In this paper we propose a stylized model wherein the current electrical system has transmuted into a likely potential future. In this future, several nearly independent subsystems, -micro-grids -, exist and interact with classical generation companies. We will make the assumption that micro-grids typically dispose of a set of generation assets mostly comprised of renewable intermittent sources (wind, solar).

Remaining energetically independent of the rest of the system is then possible, up to a certain extent, by also disposing of a set of batteries. The remaining time, i.e., in case of production surplus or lack of generation, an interaction with the classical generators will provide useful back up to meet the total demand in energy. In our model we will focus on this last type of interaction from an abstract level and suggest several convenient reformulations so as to make the model tractable. This is non-trivial, since the formulated model will be a bilevel stochastic mixed integer optimization problem. Needless to say, the price to pay for tractability is a heuristic means of solving the problem, but by no means an unrealistic one. Indeed, our suggestion solution will be "feasible"1 for the original bilevel problem, but neither be the "optimistic" nor "pessimistic" solution. From the viewpoint of practice, this is actually not unrealistic since both extremes (the "optimistic" and "pessimistic" views) are stylized anyway. Indeed there is no reason to believe that an independent actor pursuing his own goals (e.g., in the case of micro-grids, this could be to remain independent energetically speaking most of the time), will purposely select the most advantageous response for the leader, nor any reason to believe it would select the most adverse option.

Having something in between can thus be assumed reasonably realistic. Before providing a more mathematical description of our model, let us begin by introducing the key elements of it.

Parties involved. We consider an extension of electricity production problems that involve two types of parties:

• GenCos are big producers of electricity in the network (they own large assets such as nuclear, thermal, hydro and other renewable energies). In order to avoid having to face an even tougher class of problems (an equilibrium problem with equilibrium constraints, or a coupled set of MPECs, e.g., [START_REF] Surowiec | Explicit stationarity conditions and solution characterization for equilibrium problems with equilib-rium constraints[END_REF]), we focus on the decision of a single operator that faces competing companies that have a fixed, perfectly known, policy. Hence, in what follows, GenCo refers to that single operator.

• Each micro-grid q ∈ Q consists of a small subnetwork that has highly volatile generation capacities (solar, wind), and two types of demands to be attended, which we denote as hard and elastic demands. On the one hand, hard demands must be met strictly at all times. On the other hand, elastic demands (heating up water, recharging electric cars, ...) can be shifted within a certain time window. This models the fact that environmentally aware users may be willing to postpone / anticipate an essential electrical consumption in order to reduce overall environmental burden (here translated through the system cost). Since this use is essential it will take place somewhere else within a given time window. Micro-grids are assumed to be relatively autonomous in terms of energy. However, due to the uncertain nature of their production, they need to buy or sell electricity from the GenCo.

Contract. Micro-grids and GenCos interact through contracts that specify the costs of buying/selling electricity from/to the GenCo for each period of the time horizon T . Specifically, each contract k ∈ K is specified by (i) the price of contract c k paid by the micro-grid to the GenCo that proposes it, (ii) a linear cost function x → f kt x for buying the amount x of electricity during time period t, and (iii) a linear function y → g kt y for selling the amount y of electricity during time period t. We denote by K 0 the subset of contracts possibly proposed by the GenCo whose decisions are being optimized, while K also contains contracts of competing companies.

Objective. On the one hand, the objective of the GenCo is to propose the least cost production schedule based on (i) the classical costs of unit-commitment and related problems (hydro optimization and nuclear maintenance scheduling, for instance) and (ii) the cost/benefit of buying/selling electricity to the microgrids. On the other hand, the objective of each micro-grid is to minimize its total cost of buying/selling electricity to the GenCos, by choosing a contract offered by the GenCo or one of the competing companies (in the latter case, the GenCo does not produce, buy or earn anything for/from the micro-grid). These two conflicting objectives can be naturally modeled as a bilevel optimization problem, introduced in the next section.

Coupling among the different time periods. We are given a time horizon T = {1, . . . , T } that is further partitioned into days:

T = D 1 ∪ • • • ∪ D J .
Elastic loads couple all time periods of a given day. However, they do not couple time periods of different days as we may reasonably assume that the required load will be met during one full day. Nevertheless, time periods of different days may be coupled together in the presence of batteries with significant storage capacity. In this setting, the resulting multi-stage stochastic optimization problem is significantly harder, preventing us from decomposing the problem in time (by day). This is an immediate result of the temporal coupling of adjacent days through the large storage capacity batteries and the "strongly increasing complexity" of multi-stage stochastic programs with the number of stages (e.g., [START_REF] Shapiro | On complexity of multistage stochastic programs[END_REF]).

A bilevel formulation

We describe below our mathematical model, represented as a bilevel multi-stage stochastic program. Let x ∈ R Q×T represent the electricity production of the GenCo. We are mainly interested here in the interaction between the GenCo and the micro-grids, so the value x qt represents the electricity produced by the GenCo and fed into micro-grid q during period t. Similarly, we can define y qt as the amount of electricity bought by the GenCo from micro-grid q during time period t. We denote by F : R Q×2T → R ∪ {∞} the cost of producing x -y. Hence, one can think of the problem min

x,y≥0

F (x, y) (1)
as a compact abstract representation for the combination of unit-commitment, nuclear power plant maintenance planning, hydro power generation and other related problems, feeding the connected microgrids with the power described by x -y. The mapping F can take the value ∞ for given vectors x, y if certain constraints cannot be met. For instance, whenever x -y exceeds the generation capacity of the system. With this convention any constraints on generation (or other constraints) can be readily incorporated in the framework by adding the characteristic function of these constraints to F . Let us note that computing F (0, 0) amounts to solving the power generation problem without considering microgrids. In particular, computing F is therefore as difficult as computing the least cost power generation schedule. We also remark that in view of the above discussion, assuming that the GenCo problem is feasible without micro-grids means that F (0, 0) < ∞. We note that from a computational viewpoint it is better to incorporate constraints directly in the model so as to avoid having to deal with very large numerical values. This is what is done in the numerical experiments. Handling constraints by allowing F to take on the value ∞ is just convenient for the mathematical presentation.

In our model the vectors x and y indicate the power flow between the micro-grids and the (GenCo). Within the model it is not optimal to buy and sell energy at the same time. Hence at any optimal solution, x and y will not both be positive. Moreover, these variables are related to the internal functioning of the micro-grids and in particular their internal constraints. To this end we introduce an abstract constraint set M q for each q ∈ Q, modelling in particular the demand-power balance. These power flows induce a certain cost governed by a contract. In our model, the GenCo has to decide which contracts it proposes to the micro-grids. We introduce an additional binary variable Z qk that is equal to 1 if contract k is offered to micro-grid q and zero otherwise. The index set K will denote the set of all contracts and K 0 ⊆ K, the set of contracts suggested by the GenCo, i.e., K\K 0 is the set of contracts suggested by the competitors of the GenCo. We denote the set of contracts that the GenCo can offer as Z ⊆ {0, 1} |Q|×|K0| . For instance, Z could contain all binary vectors, where the GenCo offers a fixed number of contracts N q to each micro-grid q ∈ Q. In that case, we would have Z = {Z ∈ {0, 1} |Q|×|K0| : k∈K0 Z qk ≥ N q , ∀q ∈ Q}.

For each micro-grid q ∈ Q, the binary variable z qk indicates if micro-grid q subscribes to contract k. Notice that micro-grid q interacts with the GenCo only if it subscribes to one of its contracts, that is, if k∈K0 z qk ≥ 1. We will introduce local versions of the variables x, y, called x and ỹ that describe locally the power status. In particular, x = x and y = ỹ only if a contract with the GenCo is subscribed. Summarizing, the deterministic bilevel problem below considers the following optimization variables:

• Z qk : 1 if contract k is offered to micro-grid q (leader )

• z qk : 1 if contract k is subscribed by micro-grid q (follower )

• xqt : power consumed by micro-grid q during period t (follower )

• ỹqt : power produced by micro-grid q during period t (follower )

• x qt : power consumed by micro-grid q during period t from the GenCo (follower )

• y qt : power produced by micro-grid q during period t for the GenCo (follower )

We are now in measure to formulate our bilevel program:

min F (x, y) - q∈Q k∈K0 z qk (c k + t∈T (f kt x qt -g kt y qt )) (2a) s.t. Z ∈ Z (2b) (x q , y q , z q ) ∈ arg min zq.∈{0,1} K ,xq,ỹq k∈K z qk (c k + t∈T (f kt xqt -g kt ỹqt )), ∀q ∈ Q (3a) s.t. (x q , ỹq ) ∈ M q , (3b)
x q = xq k∈K0 z qk , (3c)
y q = ỹq k∈K0 z qk , (3d)
z qk ≤ Z qk , k ∈ K (3e) k∈K z qk = 1, (3f)
where the constraints (2b), (3b)-(3f) are as explained below. The objective function of the leader, (2a), minimizes the production cost already mentioned in (1) combined with transaction costs related to the micro-grids. Similarly, the objective functions of the micro-grids, (3a), amounts to minimizing their total transaction costs with respect to all GenCos, not only the one represented by the leader. Hence, (3a) involves all contracts in K while (2a) considers only the contracts in K 0 . Constraints (2b) and (3b) restrain, respectively, the set of contracts offered by the GenCo and the feasible power consumption/production for the micro-grids. Constraints (3c) and (3d) model that x q = xq and y q = ỹq if and only if micro-grid q subscribes to one of the contracts offered by the GenCo. Constraints (3e) impose that only the contract offered by the GenCo can be subscribed by the micro-grid, while constraints (3f) forces each micro-grid to select exactly one contract.

MILP reformulations

In this section we will provide two reformulations of the bilevel problem (2)- [START_REF] Taktak | An overview on mathematical programming approaches for the deterministic unit commitment problem in hydro valleys[END_REF]. In order to understand these, it is to be observed that problem (2)-( 3) is not "unequivocally" defined in the sense that it is not clear what is meant with min in (2a), unless (x q , y q , z q ) in (3a) is always unique. A classic way to remove this ambiguity is to consider a so called optimistic formulation wherein among all (x q , y q , z q ) in the lower level problem (3) optimal set, the most advantageous for the leader is chosen. The pessimistic formulation then corresponds to picking the most disadvantageous case. Obviously both notions coincide whenever the lower level problem has a unique solution. We refer to [START_REF] Dempe | Bilevel Programming Problems: Theory, Algorithms and Applications to Energy Networks[END_REF] for further information on this topic. It is generally recognized that the pessimistic formulation is significantly harder to analyze than the optimistic formulation. In section 3.1 we will suggest a reformulation that is neither optimistic nor pessimistic, but somewhere in between. We also suggest a formulation leading to optimistic solutions in section 3.2 for comparative reasons. In particular the numerical experiments carried out in section 5 will show that results from both reformulations differ, hence implying that the arg min in (3a) is indeed not unique and justifying that care should be taken in the interpretation.

Heuristic reformulation

The bilinear optimization problem described in the previous section involves binary variables and nonlinear constraints, all that being built on the top of the already difficult optimization problem [START_REF] Tahanan | Large-scale unit commitment under uncertainty: a literature survey[END_REF]. Hence, we address the problem heuristically rather than exactly and we show below how it is possible to exploit the somewhat simple linking constraints (3e) to provide a heuristic reformulation for the bilinear problem.

We propose below a one-level reformulation of the bilevel problem that may lead to solutions that are not solutions to the optimistic problem. The key aspect of our reformulation relies on pre-processing. Specifically, for each k ∈ K and q ∈ Q, we solve the restricted follower problem where z qk is equal to 1, namely:

min xq,ỹq c k + t∈T (f kt xqkt -g kt ỹqkt ) s.t. (x qk , ỹqk ) ∈ M q .
We assume that the set of contracts and micro-grids is of reasonable cardinality so that this preprocessing step remains tractable. Let (x qk , y qk ) be an optimal solution of the above problem and C qk be its solution cost. Then, a heuristic solution to the bilevel problem ( 2)-( 3), which ensures the optimality of the lower level, can be found by solving:

min z,Z F (x, y) - q∈Q k∈K0 C qk z qk (4a) s.t. x qt = k∈K0 x qkt z qk , ∀q ∈ Q, t ∈ T (4b)
y qt = k∈K0 y qkt z qk , ∀q ∈ Q, t ∈ T (4c) z qk ≤ Z qk , ∀q ∈ Q, k ∈ K (4d) k∈K z qk = 1, ∀q ∈ Q (4e) z qk ≤ 1 -Z q , ∀k ∈ K, l ∈ K, q : C qk > C ql (4f) Z ∈ Z (4g) z ∈ {0, 1} Q×K
The objective (4a) is obtained from (2a) by replacing the micro-grid cost for the GenCo by C qk z qk . Similarly, (4b) and (4c) are obtained from (3c) and (3d) by replacing the variables x and ỹ with their fixed values x and y computed in the pre-processing phase. It is exactly in these constraints that lies the heuristic aspect of our reformulation since the followers no longer choose the optimal solution that most benefits the leader, but instead take the one computed in the pre-processing phase, when (3a) admits several solutions. Finally, constraint (4f) ensures that the micro-grids choose the contracts that lead to the cheapest solutions. Indeed if Z ql = 1 and the cheap solution Cql is available, any costlier solution is ruled out by equation (4f).

Comparative exact formulation

To assess the quality of our heuristic approach, we have also devised an exact one-level reformulation for the problem where (4b) and (4c) are replaced by non-linear constraints involving the aforementioned decision variables x and ỹ that are restricted by the sets M q for each q ∈ Q. In addition, the formulation contains restrictions enforcing x and ỹ to be optimal for the subproblems.

x qt = k∈K0 xqkt z qk , ∀q ∈ Q, t ∈ T (5b)
y qt = k∈K0 ỹqkt z qk , ∀q ∈ Q, t ∈ T (5c)
C qk = c k + t∈T (f kt xqkt + g kt ỹqkt ), ∀q ∈ Q, ∀k ∈ K (5d) (x qk , ỹqk ) ∈ M q (5e) z ∈ {0, 1} Q×K
Constraints (5b) and (5c) contain bilinear terms that can be linearized using classical techniques. Specifically, let us introduce variables X qkt and Y qkt , respectively equal to products xqkt z qk and ỹqkt z qk . These variables can be substituded in constraints (5b) and (5c), adding also the constraints

X qkt ≤ xqkt X qkt ≤ M z qk X qkt ≥ xqkt -M (1 -z qk ) Y qkt ≤ ỹqkt Y qkt ≤ M z qk Y qkt ≥ ỹqkt -M (1 -z qk ),
that models the products through linear constraints involving big-M coefficients. The difficulty of the resulting MILP is that big-M coefficients often lead to numerical instability and weak LP relaxations.

In our numerical experiments we have observed no specific difficulties related to the big-M coefficients.

Stochastic extension of the model

Power generation problems that involve renewable intermittent energy like wind and solar are subject to uncertainty, since the output of the renewable power plants depends on the weather conditions. This power output is therefore only partially known and should in principle be considered uncertain. We will set up a model in this section that accounts for uncertainty in generation. The resulting bilevel stochastic optimization problem will again become tractable due to a reformulation akin to the one used in the deterministic case. Obviously the main impact of considering uncertainty is that the preprocessing step will become much harder.

Model

The most important assumption for our stochastic model below is how we set up the structure of information. The main assumption, making a reformulation as in the deterministic case feasible, is that we will restrict the choice of contracts to the first stage. This means that all contracts have to be decided prior to observing any uncertainty. The second choice involves the structure of information within the remaining stochastic program. Therein we will assume that decent weather forecasts are available for an upcoming day and hence that all information for a given day is known at the beginning of that day.

The stochastic program therefore has as many stages as there are days. In what follows we will denote with ξ j the uncertainty related to day j ∈ {1, . . . , J}. Let us summarize these assumptions:

• First stage -The GenCo chooses a set of contracts that are compatible with its objective.

-Each micro-grid chooses a contract among the contracts offered by the GenCos.

• Subsequent stages The exact productions and demands are known for all entities for all time periods that belong to the current day. Hence, the GenCo can produce the electricity and the micro-grids manage their elastic loads and batteries according to the chosen contracts and generated power.

A second important assumption that we will make is that uncertainty ξ has finite support. It can therefore be represented by a finite set of scenarios Ξ and uncertainty can be assumed to be represented by a scenario tree. For any given scenario ξ ∈ Ξ, we define x(ξ) and y(ξ) as the scenario dependent power flows between the GenCo and micro-grids. Following e.g., [START_REF] Rockafellar | Solving stochastic programming problems with risk measures by progressive hedging, Set Valued and Variational Analysis[END_REF], we choose to not make explicit the non-anticipativity constraints which can be appended to the model through an additional abstract subspace to which the decisions should belong. Obviously on the scenario tree these conditions are readily enforced. As stated earlier, the optimization problem has as many stages as days and all information within a day is assumed to be known. The introduction of uncertainty affects both the objective function and the constraints of the follower problem, as well as the objective function of the leader. Thanks to our specific measurability assumption concerning the days, uncertainty in constraints has an easy interpretation. In particular the abstract set M q (ξ) (covering for instance power balance equations) is such that M q (ξ) = J j=1 M q (ξ j ). Consequently, we can write (x q (ξ), ỹq (ξ)) ∈ M q (ξ) a.s., which for a random realization at day j, simply means that (x q (ξ j ), ỹq (ξ j )) ∈ M q (ξ j ) as a plain deterministic constraint. As uncertainty is also present in the objective function, this calls for a consideration of a tradeoff between risk and average performance (see e.g., [START_REF] Rockafellar | The fundamental risk quadrangle in risk management, optimization and statistical estimation[END_REF] for a thorough discussion on such matters). As is common in practice, we will adopt a Markowitz type of idea so that we will weigh a measure of risk and average performance. To this end we will endow ourselves with a convex weight λ, that will weigh the average performance and the risk measure for the leader (the micro-grids are assumed risk-neutral for now, we discuss how to relax this assumption later). We choose to use the CVaR as a measure of risk as it preserves convexity and moreover allows for a linear reformulation due to the following characterization (e.g., [START_REF] Rockafellar | Optimization of conditional value-at-risk[END_REF]):

CVaR [X(ξ)] = min v v + 1 1 - E[(X(ξ) -v) + ] (7)
for any random variable X and probability level ∈ (0, 1).

We can now present the following risk-averse version of the bilevel energy management problem:

min λE   F (x(ξ), y(ξ)) - q∈Q k∈K0 (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk   (8a) + (1 -λ)CVaR   F (x(ξ), y(ξ)) - q∈Q k∈K0 (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk   (8b) s.t. Z ∈ Z (8c) (x q (ξ), y q (ξ), z q ) ∈ arg min zq.∈{0,1} K ,xq,ỹq E k∈K (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk , q ∈ Q (9a) s.t. (x q (ξ), ỹq (ξ)) ∈ M q (ξ)a.s., ξ ∈ Ξ (9b) x q (ξ) = xq (ξ) k∈K0 z qk , ξ ∈ Ξ (9c) y q (ξ) = ỹq (ξ) k∈K0 z qk , ξ ∈ Ξ (9d) z qk ≤ Z qk , k ∈ K 0 , ξ ∈ Ξ (9e) k∈K z qk = 1, ξ ∈ Ξ (9f)

Heuristic reformulation

As in the deterministic case, the problem ( 8)-( 9) can be solved heuristically by solving follower problems in a pre-processing phase. For each k ∈ K and q ∈ Q, we solve the restricted follower problem where z qk is equal to 1 min xq(ξ),ỹq(ξ)

c k + E t∈T (f kt xqt (ξ) -g kt ỹqt (ξ)) (10a) s.t. (x q (ξ), ỹq (ξ)) ∈ M q (ξ)a.s., ∀ξ ∈ Ξ, (10b)
under non-anticipativity constraints. Let (x qk , y qk ) be an optimal solution and C qk be its optimal value, namely,

C qk = c k + E t∈T (f kt x qt (ξ) + g kt y qt (ξ)) .
Then, a heuristic solution to the stochastic bilevel problem ( 8)-( 9), which ensures the optimality of the lower level, can be found by solving:

min z,Z λE   F (x(ξ), y(ξ)) - q∈Q k∈K0 (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk   + (1 -λ)CVaR   F (x(ξ), y(ξ)) - q∈Q k∈K0 (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk   (11a) s.t. z qk ≤ Z qk , ∀q ∈ Q, k ∈ K 0 (11b) k∈K z qk = 1, ∀q ∈ Q (11c) x qt (ξ) = k∈K0 x qkt (ξ)z qk , ∀q ∈ Q, t ∈ T , ξ ∈ Ξ (11d) y qt (ξ) = k∈K0 y qkt (ξ)z qk , ∀q ∈ Q, t ∈ T , ξ ∈ Ξ (11e) z qk ≤ 1 -Z q , ∀k, q : C qk > C ql (11f) Z ∈ Z (11g) z ∈ {0, 1} Q×K . ( 11h
)
The above program can be further reformulated by replacing CVaR with its characterization [START_REF] Facchinei | Finite-Dimensional Variational Inequalities and Complementarity problems[END_REF] and by introducing the optimization variables w(ξ) to linearize the function [•] + . This yields the regular stochastic program:

min z,Z,w(ξ) E   λF (x(ξ), y(ξ)) -λ q∈Q k∈K0 (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk + (1 -λ) v + w(ξ) 1 -   s.t. (11b) -(11h) w(ξ) ≥ F (x(ξ), y(ξ)) - q∈Q k∈K0 (c k + t∈T (f kt x qt (ξ) -g kt y qt (ξ)))z qk -v, ∀ξ ∈ Ξ w(ξ) ≥ 0, ∀ξ ∈ Ξ.
We mention that a similar MILP reformulation can be obtained whenever the micro-grids are also assumed to be risk-averse. For instance, we could consider the same Markowitz type model and endow ourselves with a set of convex multipliers {λ q } q∈Q weighting the expectation and CVaR for each microgrid. The resulting counterpart of ( 10) is more demanding computationally due to the difficulty of decomposing the problem per scenario. This difficulty is well-known in risk-averse stochastic optimization, see for instance the recent [START_REF] Rockafellar | Solving stochastic programming problems with risk measures by progressive hedging, Set Valued and Variational Analysis[END_REF]. As this goes beyond the scope of this paper, we disregard this extension in what follows.

Comparative reformulation

As in the deterministic case, the above heuristic reformulation can be made exact by replacing constraints (11d) and (11e) with the bilinear constraints

x qt (ξ) = k∈K0 xqkt (ξ)z qk , ∀q ∈ Q, t ∈ T , ξ ∈ Ξ ( 13
)
y qt (ξ) = k∈K0 ỹqkt (ξ)z qk ∀q ∈ Q, t ∈ T , ξ ∈ Ξ ( 14
)
and adding constraints

C qk = c k + E t∈T (f kt xqt (ξ) + g kt ỹqt (ξ)) , ∀q ∈ Q, k ∈ K (15) (x qk (ξ), ỹqk (ξ)) ∈ M q (ξ), ∀q ∈ Q, k ∈ K, ξ ∈ Ξ. ( 16
)
Moreover optimization is now to be carried out over (z, Z, x(ξ), ỹ(ξ)). Evidently the resulting problem is significantly harder!

Case study based on thermal power unit-commitment

We assess below our heuristic reformulation on a power generation problem where F represents a thermal unit-commitment problem, following closely [START_REF] Carrión | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF]. The latter minimizes the total operation cost, which is defined as the sum of the production cost, the startup cost, and the shutdown cost. The production cost is piece-wise linear, while the startup cost is piece-wise linear approximation of an exponential function of the off-line time prior to the startup. Since the MILP model is rather technical, and because computing F (x, y) is not the focus of this work, we refer the interested reader to, e.g., [START_REF] Carrión | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF] for the details about the mixed-integer linear program behind F . We also refer to [START_REF] Tahanan | Large-scale unit commitment under uncertainty: a literature survey[END_REF][START_REF] Van Ackooij | Large-scale unit commitment under uncertainty: an updated literature survey[END_REF] and the many references therein for further alternatives to the model considered here. In addition, we detail the constraints used in our micro-grid model in Appendix A.

Data

General system

The time horizon is one week starting on Monday at 12 am, each day is split in periods of one hour, i.e., |T | = 168. The data used for the unit-commitment part is provided by [START_REF] Carrión | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF], 50 generators are used giving a daily fixed load of 135,5 GWh.

Considered contracts

Several types of contracts are proposed combining the possibility of lower prices during nighttime and/or the possibility of lower prices during the weekend. This leads to a total of four types of contracts. The competitors offer one contract of each type to each micro-grid. The GenCo generates four possible contracts of each type, thus |K 0 | = 16 and |K| = 20. For all contracts, the average cost of 1kW is about 0,10 EUR

Description of the micro-grids

In what follows we will define devices as an abstract group of objects fulfilling a purpose of consumption, production or storage under time constraints within a micro-grid. The devices considered, their consumption and the time periods where they are used fit numbers reported in a recent study of electricity demand/offer equilibrium in France 2 .

• Consumption devices are devices appearing in a common household, which we split into three further sub-classes:

-Constant: these devices consume power constantly. They represent 30% of the daily consumption.

-Comfort: devices that consume power during a predetermined period each day, representing 30% of the daily consumption.

-Elastic: devices for which predetermined windows of usage are defined for each day as well as a total daily load. They represent 40% of the daily consumption.

The periods where comfort devices are used and elastic devices can be used as well as their consumptions are randomly generated such that the graph of the total consumption of a micro-grid fits a classical duck curve for each day. The average consumption of 200 devices over one week is 1GWh, which corresponds to the consumption of 5000 common (European) households. Consumption of elastic devices can be reorganized to reduce costs by taking advantage of low hourly prices, production devices and storage devices. Figure 1 illustrates the hourly demand of the micro-grid composed of 10000 devices.

• Production devices are equally split into solar panels and wind turbines. The production capacity of each production device is generated randomly. For each micro-grid, their production capacity is sufficient to satisfy on average 50% of the consumption of micro-grids.

• Storage devices are of a single type, their capacity is sufficient to store on average 30% of the daily electricity load of a micro-grid. We consider that during the storage process there is a total loss of 10% of energy due to physical constraints which corresponds to the storage efficiency of state-of-the-art domestic batteries such as Tesla's Powerwall3 .

In the largest instance, 90,000 devices are considered, corresponding to over 2, 000, 000 households. The total consumption of devices in this instance is on average 64 GWh per day. As on average 50% of the daily load is supplied by production devices, about 32 GWh must be provided from the global network each day.

The penetration of the MGs in the total load of the GenCo is defined as penetration = load of the M Gs total load f or the GenCo If all MG considered choose a contract from the GenCo, the penetration rate is about 19%.

Uncertainty

Four simple independent scenarios are considered based on the possibility of good sunshine and good wind that influence production of micro-grids. Weather conditions are considered as identical each day of the week in a scenario. In bad conditions (no sunshine or no wind), a production device produces only 50% of its maximum production capacity. Notice that this implies that the multi-stage structure of the problem is simplified to a two-stage problem for which no non-anticipativity constraints are required. The value of used in CVaR is equal to 0.1, representing a high aversion to risk. The full instance details are available from the authors upon request.

Numerical results

All the methods were implemented using Julia 0.4.6 and ILOG CPLEX 12.7 as MIP solver. Tests were made on a 4-core i7 2.30 GHz processor with 16Go of RAM memory. Maximum computation time is set to 3600s.

Table 1 presents the average computation times for solving the follower problems. Times and objectives are averaged over the 20 possible contracts for the 5 instances of micro-grids of each size. The objectives are the average consumption costs for micro-grids of a given size which is also the average income for a GenCo supplying those micro-grids. They are provided in euros. The standard deviation (σ) of the objective value is also provided. These results show that solving the restricted follower problem can be done fairly quickly. The objective value for all contracts is positive and for all instances considered, there exists at least one contract of the GenCo that is cheaper and one that is more expensive than those of the competitors. Each micro-grid is financially interesting for the GenCo if the additional production costs are not too high. In terms of size, the formulation of the follower problem on a micro-grid with 100 devices and a single contract available contains 59136 variables and 35390 constraints. Figure 2 illustrates the optimal demand strategies over time of two micro-grid of 5000 devices for two contracts. The left graph represents a micro-grid composed of households and the right one a micro-grid composed of offices. One contract considers constant buying and selling prices over time (full line) and the other considers smaller prices during nighttime and weekends (dashed line). The different demand strategies will influence the production of the GenCo if the two contracts are proposed. Notice that the curves have a different aspect than a classical duck curve illustrated in Figure 1. This is due to production devices as solar panel that can be active only during daytime. As a result, the demand is negative almost all daytime and reaches its biggest values during nighttime. The right hand side graph is mostly flat during the weekend as the micro-grid considered is composed of offices with very low consumption during this period. Table 2 presents results solving small instances with the exact (Exact optimistic) and the heuristic (Heuristic) bilevel reformulations. Here, the value of λ is set to 1 in order to optimize the expectation. The first two columns report the number of devices in the considered micro-grids and the number of micro-grids of each size. Hence, the first line considers one micro-grid with 100 devices and the last line considers five micro-grids with 100 devices plus five micro-grids with 200 devices. It has been observed in the results that for each instance, both approaches (Heuristic and Exact) select the same contracts for each micro-grid. The difference between optimal values lies in the extra production cost there is for the GenCo to satisfy the demand of the micro-grids who subscribed a contract. Specifically, in Exact, the planning of elastic devices of those micro-grids can be reorganized to reduce the cost of the UC problem, preserving optimal cost for the micro-grids. This observation shows that the bilevel problem ( 8)-(9) (or (2)-( 3)) is indeed ambiguous in so much that the argmin of problem ( 9) is not unique. Consequently the pessimistic and optimistic versions of the model differ and our heuristic reformulation offers something in between. However, this gap, as reported in the column ∆(F -U C) is not very large. Indeed, let us denote by (x opt , y opt ) and (x heur , y heur ) the power generation returned by Exact and Heuristic, respectively. Columns F (x opt , y opt ) -F (0, 0) and F (x heur , y heur ) -F (0, 0) correspond to GenCo's extra production costs due to the consumption of the micro-grids in the solutions of Exact and Heuristic respectively. The value

∆(F -U C) = F (x heur , y heur ) -F (x opt , y opt ) F (x opt , y opt ) -F (0, 0)
is the relative gap between the costs of extra production in Heuristic and Exact. The reported computation times of Heuristic do not include the preprocessing times. Computation times of Heuristic (without considering preprocessing) are quite stable while those of Exact tend to grow very quickly and the last instance considering 10 micro-grids cannot be solved in one hour. Concerning the instance with a single micro-grid with 200 devices, the extra production cost is 0 because the micro-grid chose a contract of another company. The gap on extra production cost is very low in all tests. The Heuristic/Exact formulation of the instance with five micro-grids of 100 devices contains 200/300 binary variables, 1680/37255 continuous variables and 1975/279725 constraints. As no variables are related to devices in the Heuristic formulation, the number of devices does not impact its the size. Table 3 reports the results of Heuristic on larger instances. We consider three values for λ: 1, 0.5 and 0. When set to 1 (resp. 0), the formulation optimizes the expectation (resp. CVaR ). We have chosen a value of 0.1, representing a high aversion to risk. In each instance, five micro-grids of each size are used. Columns E report the expectation of the optimal solutions, columns CVaR report their conditional value-at-risk, both are given in thousands of e. Columns sold report the number of contracts of the GenCo chosen by the MG in the optimal solution found. These values are reported only for λ = 1 as they were identical for other values of λ after rounding. The first line considers only the original UC problem while the last line considers 20 micro-grids (5 of each size) with a maximum possible penetration of 19%. Notice the negative solution cost of the last line, which is obtained from a revenue from the micro-grids that is higher than the production cost. Computation times tend to grow slowly when adding microgrids and stay close to the computation time of the original UC problem when λ = 1. For λ = 0.5, the expectation is similar, the CVaR decreases significantly for some instances and the computation time increases of about 20%. With λ = 0, results for the expectation and CVaR are similar than for λ = 0.5 and computation times are more than doubled for most instances in comparaison of λ = 1. Our formulation does not increase much the difficulty of the UC problem on the tested instances when minimizing the expectation. Integrating CVaR in the objective function with λ = 0.5 seems interesting as it is not very time consuming and reduces CVaR almost to its best attainable value.

The number of contracts sold per instance reflects the quadratic aspect of production costs considered in the formulation for the UC [START_REF] Carrión | A computationally efficient mixed-integer linear formulation for the thermal unit commitment problem[END_REF]. Summing the numbers of contracts sold in instances where all microgrids have the same size (lines 2 to 5 of Table 3), a total of 19 contracts are sold to the 20 micro-grids.

When solving the problem with all 20 micro-grids together (last line of the table), only 10 contracts are sold. In this solution, respectively 2, 1, 4 and 3 contracts are sold to the MGs of size 1000, 2000, 5000 and 10,000. The results of this last instance let us suppose the GenCo cannot afford having a penetration rate higher than about 13% with the contracts proposed by the competitors. Even the small MG with 1000 devices go to competitors illustrating that the GenCo can only propose contracts more expensive than those of competitors. If many micro-grids must be supplied with electricity, the average production costs increases for the GenCo in the UC problem. At some point, the GenCo cannot afford proposing low price contracts to additional micro-grids, otherwise it would produce at a loss. Our approach prevents the GenCo from falling in this situation by proposing contracts that are more expensive than those of the competitors to micro-grids that are not financially interesting.

For each instance, the same MGs choose a contract of the GenCo for each value of λ considered but the contracts sold are not all identical. Table 4 illustrates the solution of the instance with 5 micro-grids of size 5000. The contracts are represented in the columns and are grouped by type. Constant price, daytime-nighttime price, week-weekend prices and daytime-nighttime plus week-weekend prices. A blank dot appears when a contract is proposed to a micro-grid by the GenCo, a black dot appears when a micro-grid selects a contract. The five micro-grids select a contract from the GenCo for all values of λ but the contracts proposed and chosen for the fourth and fifth micro-grid vary. Integrating CVaR can discredit contracts that present a high risk in some scenarios.

Conclusions and perspectives

We have addressed a variant of power generation optimization problems where the GenCo interacts with micro-grids by buying or selling power to them. We have modeled the problem as a bilevel stochastic optimization problem, built on the top of the already difficult power generation optimization problem that the GenCo must solve to produce its energy. For realistic size data, solving the bilevel problem exactly is out of reach so that we have focused on a heuristic reformulation that produces a solution to the bilevel problem that is neither "optimistic" nor "pessimistic". Although our numerical results exhibit that the resulting solution is indeed different from the "optimistic" solution, they also indicated that the gap is small (roughly 0.05 -0.1 % variation on the upper level objective function). Interestingly enough, we could link this discrepancy between the optimistic bilevel solution and our solution to the existence of load-shifting devices within the micro-grid that are piloted differently in the "optimistic" situation.

We have assessed our reformulation on a case-study based on a thermal unit-commitment problem using realistic data. Our numerical results have confirmed the tractability of the heuristic approach, since the solution time of our reformulation is at most 50% higher than the one required for solving the thermal unit-commitment problem. In contrast, the solution of an exact reformulation of the bilevel problem increases the solution times by more than 10 000% in the presence of only 8 micro-grids. Our results also confirm the quality of the approximation provided by our heuristic, which provides a solution very close to the solution of the exact "optimistic" model. An interesting venue for future research would seek to extend the kind of reformulations proposed in this paper to other bilevel optimization problems with integer variables, which is a class of optimization problems that are notoriously difficult to solve exactly.

C

  qk z qk (5a) s.t. (4d) -(4g)

Figure 1 :

 1 Figure 1: Duck curve example for a micro-grid with 10000 devices

Figure 2 :

 2 Figure 2: Micro-grid demand management

Table 4 :

 4 Contracts proposed and selected to five MG with |D| = 5000.

Table 1 :

 1 Average computation time and objective value for the micro-grid power planning preprocessing step.

						# Devices Avg. Obj	σ		time(s)
							100		34.99	2.25	0.15
							200		37.41	3.67	0.32
							1000		207.51	14.97	1.91
							2000		332.93	21.32	3.95
							5000		653.61	37.36	14.57
							10000		957.65	40.56	24.01
										100
		50
										50
	MW	0								MW	0
		-50								-50
			Constant							-100	Cst
		-100	Varying									Varying
		0	24	48	72	96	120	144	168		0	24	48	72	96	120	144	168
					Time periods								Time periods

Table 2 :

 2 Comparing Exact and Heuristic on small instances.

	Test set # Devices	nb	Exact F (xopt, yopt) -F (0, 0)	time(s)	Heuristic F (x heur , y heur ) -F (0, 0)	time(s)	∆(F -U C) (%)
	100	1	8.134968	76.26	8.146005	26.9	0.136
	100	2	29.413995	202.3	29.438447	28.63	0.083
	100	3	67.440871	312.46	67.471015	28.07	0.045
	100	4	52.290345	678.9	52.334735	26.67	0.085
	100	5	24.543927	690.63	24.588317	27.22	0.181
	200	1	0	135.01	0	21.16	-
	200	2	71.09899	242.34	71.1142	25.88	0.021
	200	3	128.69007	635.02	128.717851	28.05	0.022
	200	4	188.672761	1684	188.715209	25.85	0.022
	200	5	171.455653	1919.43	171.498219	25.95	0.025
	100-200	1	60.079001	223.74	60.090038	25.7	0.018
	100-200	2	100.510032	893.82	100.550933	25.93	0.041
	100-200	3	196.12615	1697.23	196.184324	27.19	0.030
	100-200	4	240.957646	3520.43	241.045643	27.38	0.037
	100-200	5	-	> 3600	196.081488	28.7	-

Table 3 :

 3 Solution of Heuristic for larger instances.

	Test set	E	CVaR	λ = 1 time(s)	sold	penetration	E	λ = 0.5 CVaR	time(s)	E	λ = 0 CVaR	time(s)
	UC	3599	-	20.98	0	0	3599	-	20.98	3599	-	20.98
	1000	3221	3580	26.31	5	1.2%	3221	3250	30.09	3221	3250	47.41
	2000	2973	3304	22.98	4	2.1%	2973	3018	26.37	2973	3018	46.71
	5000	1373	1526	24.38	5	6.0%	1374	1510	27.33	1375	1509	46.54
	10000	-421	768	24.98	5	10.9%	-421	-144	29.87	-420	-145	56.54
	1000-2000	2595	2884	25.00	9	3.3%	2595	2670	28.09	2595	2670	50.53
	1000-2000-5000	378	960	25.61	14	9.1%	379	589	30.24	379	589	50.63
	1000-2000-5000-10000	-2055	145	32.42	10	12.6%	-2055	-1768	44.19	-2053	-1769	82.29

Feasible here means that the given solution satisfies both upper and lower level constraints. Moreover at the given fixed upper level variables, the lower level variables are optimal for the lower level problem.

http://www.rte-france.com/fr/document/bilan-previsionnel-de-l-equilibre-offredemande-2016

https://www.tesla.com/fr BE/powerwall?redirect=no

Acknowledgments. The authors are grateful for the comments of two anonymous referees.

Appendix A. Micro-grids

For simplicity, this specific study considers the bus model for power flows, that is, the power network is not taken into consideration. We define two types of unit components in the micro-grid, which are called devices in the sequel. Storage devices typically represent batteries, whose status may switch between online (connected to the grid) and offline during the time horizon. During its offline periods, a storage device can be unloaded: for example, an electric vehicle is loaded during the night. During this time period, it can be used to store and serve power, but it must be fully loaded at the end of the night. During the day, the vehicle is used and its battery is emptied so that its storage level is low when it is back online. In our setting, we associate with each storage device a set of online time intervals. At the beginning of such a time interval, the charge level is a stochastic input parameter (for example, it depends on how much the vehicle has been used during the day). The required charge level at the end of the online period is modeled through more general parameters giving minimum and maximum acceptable charge levels for each online time step. We also assume that each storage device has limited capacity, charging and discharging speed and a power loss factor that is the proportion of power stored to the power consumed during the charge. The regular devices come with stochastic consumptions and productions of power during each time step. Some of their consumption can be partially delayed (elastic demand). We model this feature by defining a set of time intervals for each device (for example, a water heater must heat the water during the night). During each of them, the required total power consumption is known, as a stochastic input data. The maximum power consumption of devices is limited during each time step. The decisions to be taken in the micro-grid problem are, first, to choose a contract among those proposed by the GenCos. Then, for each regular device and each time step, the amount of elastic power consumption must be determined. For each storage device, the amounts of power consumed (to charge) and released must be fixed. We provide below a model for the problems faced by each one of the micro-grids.

Parameters. Scenario-independent parameters:

• K: set of available contracts, defined by c k , f k and g k . We assume that ∀t ∈ T , f kt ≥ g kt . -∀k ∈ K: z k = 1 if contract k is chosen by the micro-grid, 0 otherwise

x, y, w, r, s, l, u ≥ 0 (Am)

The objective of the problem (Aa) is to minimize the total cost for the micro-grid, which is composed of the fixed cost of the contract, the cost of buying power from the GenCos minus the income obtained from selling the over-production. Constraint (Ab) states that exactly one contract must be chosen. Constraints (Ac) ensure that the power flow into/out of the micro-grid is equal to its production/consumption during each time step. In the right-hand-side, the summation over D (resp. S) represents the total consumption/production of regular (resp. storage) devices. Constraints (Ad) and (Ae) define the level of power stock for each device and time step. Constraints (Af) fix the correct total amount of power that must be consumed by a device during an elastic consumption interval. The instantaneous power consumed by a device is limited by constraints (Ag). The acceptable stock levels are bound by constraints (Ah). Constraints (Ai) and (Aj) define maximum charging and discharging speeds for the storage devices, while constraints (Ak) are just a way to state that an offline device cannot be charged or discharged (the corresponding variables may as well be omitted in the model). The domains of the variables are given in constraints (Al) and (Am).