# A Rolling Horizon Approach for a Bilevel Stochastic Pricing Problem for
Demand-Side Management

Luce Brotcorne
Inria Lille-Nord Europe
Avenue du Halley 40
F - 59650 Villeneuve d’Ascq
<EMAIL_ADDRESS>Sébastien Lepaul
EDF R&D OSIRIS
Campus EdF Paris-Saclay
13 boulevard Gaspart Monge
F - 91120 Palaiseau
<EMAIL_ADDRESS>The views and opinions expressed in this article are
those of the authors and do not necessarily reflect the official policy or
position of the EDF group. Assumptions made within the analysis of this
article are not reflective of the position of the EDF group. Léonard von
Niederhäusern
Inria Lille-Nord Europe
Avenue du Halley 40
F - 59650 Villeneuve d’Ascq
<EMAIL_ADDRESS>
###### Abstract

To guarantee the well-functioning of electricity distribution networks, it is
crucial to constantly ensure the demand-supply balance. To do this, one can
control the means of production, but also influence the demand: _demand-side
management_ becomes more and more popular as the demand keeps increasing and
getting more chaotic. In this work, we propose a bilevel model involving an
energy supplier and a smart grid operator (SGO): the supplier induces shifts
of the load controlled by the SGO by offering time-dependent prices. We assume
that the SGO has contracts with consumers and decides their consumption
schedule, guaranteeing that the inconvenience induced by the load shifts will
not overcome the related financial benefits. Furthermore, we assume that the
SGO manages a source of renewable energy (RE), which leads us to consider a
stochastic bilevel model, as the generation of RE is by nature highly
unpredictable. To cope with the issue of large problem sizes, we design a
rolling horizon algorithm that can be applied in a real context.

##  1 Introduction

The efficient production, distribution, and consumption of energy are some of
the most important challenges of our time. With the growing significance of
distributed generation (DG) in the global energy mix, the production levels
are less predictable than ever. To avoid both losses or blackouts, new
solutions must be found: adapting the demand to the production instead of the
opposite has become increasingly popular to improve the efficiency of the
electricity grid’s functioning. The notion of giving the desired shape to the
demand curve is known as _demand-side management_ (DSM) (Kreith and Goswami,
2016), which can be implemented through several techniques. Among those
techniques, _load shifting_ consists in shifting a part of the demand, either
by moving forward or by postponing the consumption of electricity (Wang et
al., 2016). In this work, load shifting is the only DSM technique that is
considered.

Naturally, applying DSM techniques can only be done in a _smart grid_ context,
where extensive means of communication ensure the transmission of data among
the various actors of the electricity distribution network, specifically
between energy suppliers and consumers (Farhangi, 2010; Kabalci, 2016).

In Alekseeva et al. (2018), a power supplier aims to maximize its profit,
knowing that its clients are going to optimize their consumption accordingly
to the prices that the supplier offers. A similar pattern is found in Afşar et
al. (2016), where the supplier’s objective consists in minimizing the peak
load. Finally, in Aussel et al. (2020), a model involving an energy supplier
in a best response situation and various types of clients is considered.
Common elements of those last three works are twofold. First, they all feature
an energy supplier that induces a load shifting from its clients through price
incentives. Second, they all rely on bilevel optimization to reach their
objectives.

  * •

Pricing can be a strong tool to incentivize the buyers to adopt a certain
behavior, as in general augmenting the prices will lead to a decrease in
demand, and vice versa. Numerous applications can be contemplated: from peer-
to-peer networks (Park and van der Schaar, 2010) to freight delivery (Holguín-
Veras and Aros-Vera, 2015) via communication networks (Ozdaglar and Srikant,
2007). In the electricity domain, a personalized real-time pricing mechanism
that aims to optimize the system’s functioning while preserving the users’
welfare is proposed in Tsaousoglou et al. (2019). Closer to this work, in Liu
et al. (2019), pricing-based demand response is implemented for a smart home
with various types of household appliances, together with energy storage units
and DG, and taking into account the consumer’s welfare, but without a bilevel
structure.

  * •

Bilevel optimization originates in the seminal work of Stackelberg (von
Stackelberg, 1934), and is used to model all kinds of hierarchical
interactions between a _leader_ making decisions first and a _follower_
reacting optimally to the leader’s decisions (general references to bilevel
programming include Bard (2010); Dempe (2002)). In the three works mentioned
above, the electricity supplier always acts as the leader: it offers time-
dependent prices, to which the clients and/or customers react in an optimal
way by changing their load demand. Although proven to be NP-hard even in the
simplest cases where everything is linear (Ben-Ayed and Blair, 1990; Labbé et
al., 1998), bilevel programs are appropriate tools, as they allow to
explicitly take into account the response of the follower to the leader’s
decisions.

Besides the three works mentioned above that combine demand response and
bilevel optimization, several researchers have investigated this specific
setting (Besançon et al., 2018; Kovács, 2019; Yuan et al., 2020; Shomalzadeh
et al., 2020). Next to those works, in Grimm et al. (2020), several pricing
schemes (time-of-use, critical-peak-pricing, real-time-pricing, and fixed-
price) are compared in a retailer-prosumer interaction. Two bilevel models are
presented in Alves et al. (2020) in the aim to optimize time-of-use prices. In
the first model, the periods where modified prices apply are predetermined,
whereas in the second model, they are variables as well. Finally, in Tang et
al. (2019), an interaction between a power grid and individual buildings is
shown, where the grid aims to optimize its profit and reduce the demand
fluctuation, while the buildings minimize their bill while modifying their
demand as less as possible. A more detailed overview of demand response
problems modeled through bilevel programs consists in Henggeler Antunes et al.
(2020).

In this article, we present the problem (SBPP), which is strongly inspired by
the problems studied in Afşar et al. (2016): an energy supplier aims to
maximize its profit (its sales minus its costs) by selling electricity to a
smart grid operator (SGO) that manages the consumption schedule of a set of
clients. The SGO aims to minimize a utility function that consists in the sum
of its clients’ electricity bills plus a so-called _inconvenience cost_. We
assume indeed that shifting the clients’ loads comes at a certain cost: if the
use of a client’s device has to be postponed, the client will undergo some
inconvenience. To satisfy its clients’ demand, we finally hypothesize that the
SGO manages a source of DG in the form of photovoltaic panels, and some
storage capacities - a configuration is considered in Xu et al. (2020) without
the bilevel setting. Specifically, managing the DG induces some difficulties,
as it is impossible to know with complete precision the energy that will be
produced. This is why (SBPP) involves stochasticity, under the form of
scenario trees. The models resulting from this method becoming exponentially
large as the size of the scenario trees grows, we present a _rolling horizon_
algorithm that does not solve a problem with a large scenario tree, but can be
applied in a real setting to decide the supplier’s prices on the go. Although
rolling horizons have been known and used for decades in the framework of
stochastic optimization (Sethi and Sorger, 1991), with various applications
e.g., in scheduling (Samà et al., 2013) or in vehicle routing (Crama and
Pironet, 2019), they have not been used yet in the framework of bilevel
optimization, to the notable exception of Kallabis et al. (2019). In the
latter, the leader is an investor who decides on investments in power
generation assets, and the follower is a market operator maximizing welfare
given consumer demand and installed generation assets.

This article is organized as follows: in Section 2, the problem (SBPP) is
first defined in its deterministic form, that is without scenario trees, for
pedagogical reasons. In Section 3, we introduce the stochasticity by exposing
the modifications implied on the deterministic form of (SBPP), and present a
one-level formulation of the problem in order to solve it computationally.
Section 4 is dedicated to the rolling horizon algorithm. In Section 5,
numerical results support our theoretical work: first, an analysis of the
deterministic form of (SBPP) justifies our approach, showing that pricing
incentives can have interesting outcomes for an energy supplier, then the
rolling horizon algorithm is tested and proves to be applicable in a real
setting. Finally, some conclusions are drawn in Section 6.

##  2 Deterministic Form of (SBPP)

###  2.1 Follower’s Problem

The follower is a smart grid operator (SGO) that has contracts with a set of
clients $N$. Each of the clients $n\in N$ has available a set of devices
$A_{n}$ that must be powered during an associated time window $T_{(n,a)}$. The
device $(n,a)$ ($a\in A_{n}$ for $n\in N$) must receive a quantity of energy
$E_{(n,a)}$ during its associated time window, but cannot receive more than
$\beta^{max}_{(n,a)}$ energy units. To power the devices, the SGO can take
energy from four sources:

  * •

$x^{h}_{(n,a)}$: energy purchased from the leader,

  * •

$\bar{x}^{h}_{(n,a)}$: energy purchased from the competitor,

  * •

$\lambda^{h}_{(n,a)}$: energy taken from the DG,

  * •

$s^{h}_{(n,a)}$: energy taken from the storage.

Therefore, for each device $(n,a)$, the SGO must satisfy the constraints:

$\displaystyle\displaystyle\sum_{h\in
T_{(n,a)}}\left(x^{h}_{(n,a)}+\bar{x}^{h}_{(n,a)}+\lambda^{h}_{(n,a)}+s^{h}_{(n,a)}\right)\geq
E_{(n,a)}$ (1) $\displaystyle
x^{h}_{(n,a)}+\bar{x}^{h}_{(n,a)}+\lambda^{h}_{(n,a)}+s^{h}_{(n,a)}\leq\beta^{max}_{(n,a)}$
$\displaystyle\forall h\in T_{(n,a)}.$ (2)

Next, the SGO manages a storage capacity in the form of a battery. The battery
state $S^{h}$ must be actualized at each hour, and naturally cannot exceed the
battery capacity, or be lower than the minimal required amount of stored
energy. The battery state being known at time $0$, the follower must ensure
that the following constraints are satisfied for all $h\in\mathcal{H}$:

$\displaystyle\displaystyle S^{0}=S^{start}$ (3) $\displaystyle\displaystyle
S^{h+1}=\rho^{d}S^{h}-\sum_{\begin{subarray}{c}n\in N\\\ a\in
A_{n}\end{subarray}}s^{h}_{(n,a)}+\rho^{c}\left(\lambda^{h}_{s}+x^{h}_{s}+\bar{x}^{h}_{s}\right)$
(4) $\displaystyle\displaystyle S^{\min}\leq S^{h}\leq S^{\max}.$ (5)

Here, $\rho^{c}$ and $\rho^{d}$ are respectively the charging and the
discharging coefficients of the battery, and $S^{\min}$ and $S^{\max}$ are the
lower and upper bound on the battery state. Furthermore, the variables
$\lambda^{h}_{s}$, $x^{h}_{s}$ and $\bar{x}^{h}_{s}$ represent the stored
energy respectively taken from the DG, purchased from the leader, and
purchased from the competitor.

Naturally, it is impossible to take more energy from the battery than the
quantity that is stored at the beginning of a time slot. Therefore, the
constraint

$\displaystyle\displaystyle\sum_{\begin{subarray}{c}n\in N\\\ a\in
A_{n}\end{subarray}}s^{h}_{(n,a)}\leq S^{h}$ (6)

must be satisfied for all $h\in\mathcal{H}$.

The last constraint of the follower concerns the DG: only the energy that is
indeed produced can be consumed or stored. Therefore, for all
$h\in\mathcal{H}$,

$\displaystyle\displaystyle\lambda^{h}_{s}+\sum_{\begin{subarray}{c}n\in N\\\
a\in A_{n}\end{subarray}}\lambda^{h}_{(n,a)}\leq\lambda^{h}_{max}$ (7)

has to hold.

Finally, the objective of the SGO is to minimize the generalized cost
consisting in the sum of the purchase costs and the inconvenience caused by
delaying the use of a device. We indeed assume that the clients prefer to use
their devices at the beginning of their respective time windows. To represent
this inconvenience, we introduce the coefficient $C^{h}_{(n,a)}$, that,
multiplied by the energy consumption of device $(n,a)$, gives the
inconvenience caused by the usage of $(n,a)$ at time $h\in T_{(n,a)}$. The
objective function of the follower is thus:

$\displaystyle\begin{array}[]{l}\displaystyle\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}C^{h}_{(n,a)}\left(x^{h}_{(n,a)}+\bar{x}^{h}_{(n,a)}+\lambda^{h}_{(n,a)}+s^{h}_{(n,a)}\right)\\\
\displaystyle+\sum_{\begin{subarray}{c}n\in N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}\left(p^{h}x^{h}_{(n,a)}+\bar{p}^{h}\bar{x}^{h}_{(n,a)}\right)+\sum_{h\in\mathcal{H}}\left(p^{h}x^{h}_{s}+\bar{p}^{h}\bar{x}^{h}_{s}\right).\end{array}$
(10)

In summary, the follower’s problem can be formulated as
$\left(P_{\operatorname{SGO}}\right)$:

$\begin{array}[]{c
c}\displaystyle\min_{\mathbf{x},\mathbf{\bar{x}},\mathbf{s},\boldsymbol{\lambda},\mathbf{S}\geq
0}&(\ref{fol_obj})\\\
\mathrm{s.t.}&(\ref{fol_con1})-(\ref{fol_con7}),\end{array}$

where bold variables are used to denote vectors (whenever indices and
exponents are clear from the context).

###  2.2 Leader’s Problem

The leader is an energy supplier who aims to maximize their revenue, which is
the difference between their sales and their costs. These costs can be
production costs, or, in our case, a linear function $K^{h}$ that represents
the purchase costs on the spot market. The prices of the spot market are here
supposed to be known in advance. The problem of the leader is thus
$(P_{\operatorname{S}})$:

$\begin{array}[]{r
l}\displaystyle\max_{\mathbf{p}}\max_{\mathbf{x},\mathbf{\bar{x}},\mathbf{s},\boldsymbol{\lambda},\mathbf{S}}&\displaystyle\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}p^{h}x^{h}_{(n,a)}+\sum_{h\in\mathcal{H}}p^{h}x^{h}_{s}-\sum_{h\in\mathcal{H}}K^{h}\left(x^{h}_{s}+\sum_{\begin{subarray}{c}n\in
N,a\in A_{n}\\\ \mathrm{s.t.\ }h\in
T_{(n,a)}\end{subarray}}x^{h}_{(n,a)}\right)\\\
\mathrm{s.t.}&\left(\mathbf{x},\mathbf{\bar{x}},\mathbf{s},\boldsymbol{\lambda},\mathbf{S}\right)\mathrm{\
solves\ }\left(P_{\operatorname{SGO}}\right).\end{array}$

Observe in particular that $(P_{\operatorname{S}})$ is the optimistic
formulation of the supplier’s problem, as the maximum over the follower’s
variables is taken. Remark as well that there are no constraints on the
leader’s variables, as upper and lower bounds on the prices are implicit.
Negative prices would indeed induce losses for the leader, whereas prices
higher than the competitor’s prices imply that the follower buys from the
competitor.

##  3 Scenario Tree Approach

To tackle the issue of having a stochastic bound in the follower’s problem, we
adopt a scenario tree approach. Scenario trees are a widely used method in
classical stochastic optimization (see e.g. Heitsch and Römisch (2009)),
although their application to bilevel stochastic problems remains
confidential.

Since the unknown information is part of the lower level, the lower level is
actually a multistage problem; thus the scenario trees are introduced into the
lower level problem. At the upper level, the prices of the leader are decided
at the beginning of the time horizon and do not change at a later stage, thus
do not vary depending on the scenario.

Formally, in our problem, a scenario $\sigma\in\Sigma$ takes the form of a
vector of bounds on the DG
$\left(\lambda^{1,\sigma}_{\max},\dots,\lambda^{H,\sigma}_{\max}\right)$,
where for all $h\in\mathcal{H}$,
$\lambda^{h,\sigma}_{\max}=\lambda^{h,\sigma_{i}}_{\max}$ for some
$i\in\\{1,\dots,n_{\Sigma}\\}$. The scenarios $\sigma_{i}$ for
$i\in\\{1,\dots,n_{\Sigma}\\}$ are called the _base_ scenarios: each scenario
$\sigma\in\Sigma$ switches from one base scenario to the other as time passes
by. The set of possible scenarios can thus be represented as the leaves of a
rooted tree, where each node has one ascendant and $n_{\Sigma}$ descendants.
Furthermore, the scenario $\sigma\in\Sigma$ has a probability $P[\sigma]$ to
occur.

Introducing scenario trees into the follower’s problem has multiple
repercussions. First, for each scenario $\sigma\in\Sigma$, we associate a set
of follower’s variables:

$\displaystyle\begin{array}[]{l c l l c
l}x^{h}_{(n,a)}&\mapsto&x^{h,\sigma}_{(n,a)}&x^{h}_{s}&\mapsto&x^{h,\sigma}_{s}\\\
\bar{x}^{h}_{(n,a)}&\mapsto&\bar{x}^{h,\sigma}_{(n,a)}&\bar{x}^{h}_{s}&\mapsto&\bar{x}^{h,\sigma}_{s}\\\
\lambda^{h}_{(n,a)}&\mapsto&\lambda^{h,\sigma}_{(n,a)}&\lambda^{h}_{s}&\mapsto&\lambda^{h,\sigma}_{s}\\\
s^{h}_{(n,a)}&\mapsto&s^{h,\sigma}_{(n,a)}&S^{h}&\mapsto&S^{h,\sigma}.\\\
\end{array}$ (15)

Naturally, constraints 1-6 must hold for all scenarios $\sigma\in\Sigma$, with
the variables replacement described in 15. We denote those constraints by
(1)($\sigma$)-(6)($\sigma$). The upper bound on the renewable energy
production in constraint (7) depends on the scenario; therefore the new
constraints (7)($\sigma$) are as follows:

$\displaystyle\lambda^{h,\sigma}_{s}+\sum_{\begin{subarray}{c}n\in N\\\ a\in
A_{n}\end{subarray}}\lambda^{h,\sigma}_{(n,a)}\leq\lambda^{h,\sigma}_{\max}\qquad\forall
h\in\mathcal{H}.$

Since scenario trees are considered, it is crucial that as long as two
scenarios are indistinguishable, the same decisions are made. Thus we define

$h(\sigma,\sigma^{\prime})=\max\left\\{h\in\mathcal{H}\mid\lambda^{h^{\prime},\sigma}_{\max}=\lambda^{h^{\prime},\sigma^{\prime}}_{\max}\quad\forall
h^{\prime}\leq h\right\\},$

which is the latest time slot for which the scenarios $\sigma$ and
$\sigma^{\prime}$ are indistinguishable. For all pairs of scenarios
$\sigma,\sigma^{\prime}\in\Sigma$, the so-called _nonanticipativity_
constraints have to be satisfied:

$\displaystyle\begin{array}[]{l c l l c
l}x^{h,\sigma}_{(n,a)}&=&x^{h,\sigma^{\prime}}_{(n,a)}&x^{h,\sigma}_{s}&=&x^{h,\sigma^{\prime}}_{s}\\\
\bar{x}^{h,\sigma}_{(n,a)}&=&\bar{x}^{h,\sigma^{\prime}}_{(n,a)}&\bar{x}^{h,\sigma}_{s}&=&\bar{x}^{h,\sigma^{\prime}}_{s}\\\
\lambda^{h,\sigma}_{(n,a)}&=&\lambda^{h,\sigma^{\prime}}_{(n,a)}&\lambda^{h,\sigma}_{s}&=&\lambda^{h,\sigma^{\prime}}_{s}\\\
s^{h,\sigma}_{(n,a)}&=&s^{h,\sigma^{\prime}}_{(n,a)}&S^{h,\sigma}&=&S^{h,\sigma^{\prime}}\\\
\end{array}$ $\displaystyle\forall h\leq h(\sigma,\sigma^{\prime}).$ (20)

The last necessary change induced by the introduction of scenario trees
concerns the objective functions of both the leader and follower. First, the
follower’s objective function is now defined as

$\displaystyle\begin{array}[]{l}\displaystyle\sum_{\sigma\in\Sigma}P[\sigma]\left(\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}C^{h}_{(n,a)}\left(x^{h,\sigma}_{(n,a)}+\bar{x}^{h,\sigma}_{(n,a)}+\lambda^{h,\sigma}_{(n,a)}+s^{h,\sigma}_{(n,a)}\right)\right.\\\
\displaystyle\hphantom{\sum_{\sigma\in\Sigma}P[\sigma]}\left.+\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}\left(p^{h}x^{h,\sigma}_{(n,a)}+\bar{p}^{h}\bar{x}^{h,\sigma}_{(n,a)}\right)+\sum_{h\in\mathcal{H}}\left(p^{h}x^{h,\sigma}_{s}+\bar{p}^{h}\bar{x}^{h,\sigma}_{s}\right)\right),\end{array}$
(23)

where $P[\sigma]$ is the occurring probability of $\sigma\in\Sigma$, so that
$\sum_{\sigma\in\Sigma}P[\sigma]=1$.

The follower’s stochastic problem is finally defined as
$\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right)$:

$\begin{array}[]{r
c}\displaystyle\min_{\mathbf{x},\mathbf{\bar{x}},\mathbf{s},\boldsymbol{\lambda},\mathbf{S}\geq
0}&(\ref{fol_obj_sto})\\\ \mathrm{s.t.}&\left\\{\begin{array}[]{c
l}(\ref{fol_con1})(\sigma)-(\ref{fol_con7})(\sigma)&\forall\sigma\in\Sigma\\\
(\ref{cons_nonant})&\forall\sigma,\sigma^{\prime}\in\Sigma.\end{array}\right.\end{array}$

At the upper level, the leader also optimizes their profit expectancy.
Therefore, their objective function becomes

$\displaystyle\sum_{\sigma\in\Sigma}P[\sigma]\left(\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}p^{h}x^{h,\sigma}_{(n,a)}+\sum_{h\in\mathcal{H}}p^{h}x^{h,\sigma}_{s}-\sum_{h\in\mathcal{H}}K^{h}\left(x^{h,\sigma}_{s}+\sum_{\begin{subarray}{c}n\in
N,a\in A_{n}\\\ \mathrm{s.t.\ }h\in
T_{(n,a)}\end{subarray}}x^{h,\sigma}_{(n,a)}\right)\right),$ (24)

and the leader’s stochastic problem is defined as
$\left(P^{\operatorname{sto}}_{\operatorname{S}}\right)$:

$\begin{array}[]{r
c}\displaystyle\max_{\mathbf{p}}\max_{\mathbf{x},\mathbf{\bar{x}},\mathbf{s},\boldsymbol{\lambda},\mathbf{S}}&(\ref{lead_obj_sto})\\\
\mathrm{s.t.}&\left(\mathbf{x},\mathbf{\bar{x}},\mathbf{s},\boldsymbol{\lambda},\mathbf{S}\right)\mathrm{\
solves\ }\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right).\end{array}$

###  3.1 One-level formulation

Now that the stochastic problem is defined, we intend to solve it. This can be
achieved by transforming it into an equivalent MIP that can be computationally
tackled by (commercial) solvers. Observe that for fixed decisions of the
leader, the follower’s problem is a linear one. Therefore, it can be replaced
by its optimality conditions in the leader’s problem, giving rise to a
mathematical problem with complementarity constraints (MPCC), which can be
linearized through the _big M_ method, yielding the expected MIP. This method
is widely used in bilevel programming, since MIPs can be solved efficiently by
commercial solvers. The optimality conditions consist in the primal
constraints, the dual constraints and the complementarity slackness
constraints.

Let us denote a tuple of primal variables by

$\varphi^{p}=\left(\mathbf{x},\mathbf{\bar{x}},\boldsymbol{\lambda},\mathbf{s},\mathbf{x_{s}},\mathbf{\bar{x}_{s}},\boldsymbol{\lambda_{s}},\mathbf{S}\right).$

This tuple $\varphi^{p}$ belongs to $Primal$ if it satisfies the constraints
(1)($\sigma$)-(7)($\sigma$) for all $\sigma\in\Sigma$ and (20) for all pairs
of scenarios $\sigma,\sigma^{\prime}\in\Sigma$. To each of the constraints, we
associate its dual variable: respectively $d1^{\sigma}_{(n,a)}$,
$d2^{\sigma}_{(n,a)}$, $d3^{\sigma}_{0}$, $d4^{h,\sigma}$,
$d5^{h,\sigma}_{min}$, $d5^{h,\sigma}_{max}$, $d6^{h,\sigma}_{max}$,
$d7^{h,\sigma}$, $dx^{\sigma,\sigma^{\prime},h}_{(n,a)}$,
$d\bar{x}^{\sigma,\sigma^{\prime},h}_{(n,a)}$,
$d\lambda^{\sigma,\sigma^{\prime},h}_{(n,a)}$,
$ds^{\sigma,\sigma^{\prime},h}_{(n,a)}$,
$d\lambda^{\sigma,\sigma^{\prime},h}_{s}$,
$dx^{\sigma,\sigma^{\prime},h}_{s}$,
$d\bar{x}^{\sigma,\sigma^{\prime},h}_{s}$, and
$dS^{\sigma,\sigma^{\prime},h}$. Thus, we denote a tuple of dual variables by

$\varphi^{d}=\left(\mathbf{d1},\mathbf{d2},\mathbf{d3},\mathbf{d4},\mathbf{d5},\mathbf{d6},\mathbf{d7},\mathbf{dx},\mathbf{d\bar{x}},\mathbf{d}\boldsymbol{\lambda},\mathbf{ds},\mathbf{d}\boldsymbol{\lambda_{s}},\mathbf{dx_{s}},\mathbf{d\bar{x}_{s}},\mathbf{dS}\right),$

and write $\varphi^{d}\in Dual$ if $\varphi^{d}$ satisfies the dual
constraints of $\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right)$,
which are described in detail in A.

Finally, let us denote by $\varphi$ a pair $(\varphi^{p},\varphi^{d})$ of
primal and dual variables. This tuple belongs to $Comp$ if it satisfies the
complementarity constraints detailed in B. If $(\varphi^{p},\varphi^{d})$
belongs to $\left(Primal\times Dual\right)\cap Comp$, then $\varphi^{p}$ is
primal optimal and $\varphi^{d}$ is dual optimal. Therefore, the leader’s
problem $\left(P^{\operatorname{sto}}_{\operatorname{S}}\right)$ can be
replaced by the single-level problem:

$\displaystyle\max_{\mathbf{p},\varphi}\displaystyle\sum_{\sigma\in\Sigma}P[\sigma]\cdot\left(\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}p^{h}x^{h,\sigma}_{(n,a)}+\sum_{h\in\mathcal{H}}p^{h}x^{h,\sigma}_{s}\right.$
$\displaystyle\left.\hphantom{\max_{\mathbf{p},\varphi}\displaystyle\sum_{\sigma\in\Sigma}P[\sigma]}-\sum_{h\in\mathcal{H}}K^{h}\left(x^{h,\sigma}_{s}+\sum_{\begin{subarray}{c}n\in
N,a\in A_{n}\\\ \mathrm{s.t.\ }h\in
T_{(n,a)}\end{subarray}}x^{h,\sigma}_{(n,a)}\right)\right)$
$\displaystyle\mathrm{s.t.}\ \varphi\in\left(Primal\times Dual\right)\cap
Comp.$

Nevertheless, this problem is nonlinear, due to the products of variables in
the complementarity constraints and in the objective function. To get a mixed
integer linear problem (MILP), we linearize those products. First, the dual
objective function of the follower’s problem is used. Thanks to strong
duality, an optimal pair $\left(\varphi^{p},\varphi^{d}\right)$ satisfies that

$\displaystyle(\ref{fol_obj_sto})=$
$\displaystyle\sum_{\sigma\in\Sigma}\left(\sum_{\begin{subarray}{c}n\in N\\\
a\in A_{n}\end{subarray}}\left(E_{(n,a)}d1^{\sigma}_{(n,a)}-\sum_{h\in
T_{(n,a)}}\beta^{max}_{(n,a)}d2^{\sigma}_{(n,a)}\right)\right.$
$\displaystyle\hphantom{\sum_{\sigma\in\Sigma}}\left.+\sum_{h\in\mathcal{H}}\left(-\lambda^{h}_{max}d7^{h,\sigma}+S^{\min}d5^{h,\sigma}_{min}-S^{\max}d5^{h,\sigma}_{max}\right)\vphantom{\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\end{subarray}}}\right).$

It follows that the leader’s objective function $(\ref{lead_obj_sto})$ can be
rewritten as

$F(\varphi)=\left\\{\begin{array}[]{l}\displaystyle\sum_{\sigma\in\Sigma}\left(\vphantom{\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in T_{(n,a)}\end{subarray}}}\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\end{subarray}}\left(E_{(n,a)}d1^{\sigma}_{(n,a)}-\sum_{h\in
T_{(n,a)}}\beta^{max}_{(n,a)}d2^{\sigma}_{(n,a)}\right)\right.\\\
\displaystyle+\sum_{h\in\mathcal{H}}\left(-\lambda^{h}_{max}d7^{h,\sigma}+S^{\min}d5^{h,\sigma}_{min}-S^{\max}d5^{h,\sigma}_{max}\right)\\\
\displaystyle-\left(\sum_{\begin{subarray}{c}n\in N\\\ a\in A_{n}\\\ h\in
T_{(n,a)}\end{subarray}}\left(\bar{p}^{h}\bar{x}^{h,\sigma}_{(n,a)}+C^{h}_{(n,a)}\left(x^{h,\sigma}_{(n,a)}+\bar{x}^{h,\sigma}_{(n,a)}+\lambda^{h,\sigma}_{(n,a)}+s^{h,\sigma}_{(n,a)}\right)\right)\right.\\\
\displaystyle\left.\left.+\sum_{h\in\mathcal{H}}\left(\bar{p}^{h}\bar{x}^{h,\sigma}_{s}+K^{h}\left(x^{h,\sigma}_{s}+\sum_{\begin{subarray}{c}n\in
N\\\ a\in
A_{n}\end{subarray}}x^{h,\sigma}_{(n,a)}\right)\right)\vphantom{\sum_{\begin{subarray}{c}n\in
N\\\ a\in A_{n}\\\ h\in T_{(n,a)}\end{subarray}}}\right)\right).\\\
\end{array}\right.$

Furthermore, the complementarity constraints are linearized with the big $M$
method, which consists in replacing the equation system

$\displaystyle\begin{array}[]{l}a\cdot b=0\\\ a,b\geq 0\\\ \end{array}$ (27)

by

$\displaystyle\begin{array}[]{l}a\leq M\cdot\delta\\\ b\leq
M\cdot\left(1-\delta\right)\\\ a,b\geq 0,\\\ \delta\in\\{0,1\\},\\\
\end{array}$ (32)

with $M$ sufficiently large. This process is not anodyne. Taking $M$ too small
might eliminate solutions, and even careful algorithms can be fooled, as shown
by Pineda and Morales (2019). On the other hand, choosing $M$ too large will
likely decrease the efficiency of the solver. Finally, choosing the right $M$
is NP-hard, as demonstrated by Kleinert et al. (2019). In our case, the big
$M$ is chosen empirically. If $\varphi$ satisfies the set of big $M$
constraints replacing the original complementarity constraints, we write
$\varphi\in Comp^{\prime}$. Hence, the leader’s problem
$\left(P^{\operatorname{sto}}_{\operatorname{S}}\right)$ can finally be solved
by finding an optimal solution of Problem
$\left(P^{\operatorname{sto},\operatorname{MIP}}_{\operatorname{S}}\right)$:

$\displaystyle\max_{\mathbf{p},\varphi}F\left(\varphi\right)$
$\displaystyle\mathrm{s.t.\ }\varphi\in\left(Primal\times Dual\right)\cap
Comp^{\prime}.$

Obviously, Problem
$\left(P^{\operatorname{sto},\operatorname{MIP}}_{\operatorname{S}}\right)$
becomes huge when the number of time periods increases. Therefore, solving
this MIP is not an option for instances comprising a large time horizon. This
fact motivates the following section, where we design a rolling horizon
algorithm.

##  4 Rolling Horizon Approach

Methods using rolling horizons are well-known in the topic of stochastic
optimization, see for example Pironet (2014) for applications in
transportation management. The idea behind rolling horizons is that it is
possible to build a solution for a large problem out of the solutions of a
sequence of smaller problems. More precisely, this sequence of smaller
problems follows a chronological logic: the first problem to be solved
considers a small time horizon from the beginning of the whole time horizon
($0$) to another close time slot ($l_{RH}$). Then, the second problem will
consider a time horizon going from $s_{RH}$ to $l_{RH}+s_{RH}$ and so on,
justifying the naming of _rolling horizon_.

We first define the main parameters:

  * •

The (whole) _time horizon_ is the set $\mathcal{H}=\\{0,1,\dots,H\\}$, where
$H$ can be infinite.

  * •

The _length_ of the rolling horizon $l_{RH}$ is the number of time slots that
will be considered at every iteration of the rolling horizon method. Clearly,
$l_{RH}\leq H$.

  * •

The rolling horizon _step_ $s_{RH}$ is the number of time slots between two
instances of the rolling horizon method. More precisely, if at the $n$th
iteration, the time horizon of the considered instance starts at $t$, then the
time horizon of the $n+1$st iteration will start at $t+s_{RH}$. To have
sensible parameters, it is necessary that $s_{RH}\leq l_{RH}$.

  * •

The length of the frozen horizon $l_{FH}$ determines how much of the results
computed at the $n-1$st iteration will be reused in the $n$th iteration. In
classical stochastic optimization, it means that at the $n$th iteration, on
the time horizon going from $t$ to $t+l_{RH}$, the decisions corresponding to
the time slots $t,\dots,t+l_{FH}$ were made during the $n-1$st iteration and
are not recomputed. In bilevel programming, the situation is slightly more
complex. Indeed, the lower level variables can be seen as recourse variables,
and should thus be modifiable at any iteration of the rolling horizon method.
On the other hand, the energy supplier offers prices to the SGO, and thus
should not be able to modify their prices at the last minute. Hence, during
the $n$th iteration of the rolling horizon method, the leader’s prices will be
fixed for the time slots $t,\dots,t+l_{FH}$, and only the prices corresponding
to later time slots might change, along with the SGO-decided consumption
schedule. It follows from the definition of $l_{FH}$ that $l_{RH}-l_{FH}\leq
s_{RH}$ must hold.

Figure 1 illustrates the situation.

Figure 1: The rolling horizon scheme

Let us denote by $P_{\mathcal{H}}$ the set of parameters defining an instance
of $(SBPP)$ over the full time horizon $\mathcal{H}$. $P_{\mathcal{H}}$ thus
comprises the information related to the devices, to the competitor’s prices,
to the storage capacities, and to the distributed generation scenarios. The
rolling horizon algorithm is described in Algorithm 1.

Algorithm 1 Rolling horizon algorithm

1:procedure RH($P_{\mathcal{H}}$, $l_{RH}$, $l_{FH}$, $s_{RH}$)

2: $t\leftarrow 0$

3: repeat

4: Generate $P^{t}_{s_{RH}}$

5: Solve $P^{t}_{s_{RH}}$

6: Select a scenario $\sigma$

7: Actualize data

8: $t\leftarrow t+s_{RH}$

9: until $t+l_{RH}\geq H$

10: return data

11:end procedure

In step 4, the subinstance $P^{t}_{s_{RH}}$ is generated. The time period
considered in $P^{t}_{s_{RH}}$ is the time period
$\mathcal{H}_{t}=\\{t,\dots,t+l_{RH}\\}$. Therefore, the only devices in
$P^{t}_{s_{RH}}$ are the devices $(n,a)$ such that
$T_{(n,a)}\cap\\{t,\dots,t+l_{RH}\\}\neq\emptyset$. Two cases may arise for a
given device $(n,a)$:

  * •

If $T_{(n,a)}^{\mathrm{first}}<t$, then the device was already considered in
the previous iterations of the algorithm. Therefore, it might have received
some energy during the previous time slots. It results that the required
energy has to be actualized: the energy demand of $(n,a)$ in $P^{t}_{s_{RH}}$
is denoted by $E_{(n,a)}^{t}$ and is worth
$E_{(n,a)}-\sum_{h<t}\left(x^{h,\sigma}_{(n,a)}+\bar{x}^{h,\sigma}_{(n,a)}+s^{h,\sigma}_{(n,a)}+\lambda^{h,\sigma}_{(n,a)}\right)$,
with $\sigma$ the scenario that has been selected at step 6 in the previous
iterations.

  * •

If $T_{(n,a)}^{\mathrm{last}}>t+l_{RH}$, then $T_{(n,a)}$ partly belongs to
the time period considered in $P^{t}_{s_{RH}}$. Several options are
conceivable. Here, we choose to set
$E_{(n,a)}^{t}=\min\left\\{E_{(n,a)},\left(t+l_{RH}-T_{(n,a)}^{\mathrm{first}}\right)\beta^{max}_{(n,a)}\right\\}$,
meaning that the device $(n,a)$ should consume as much energy as possible
during the time period $\mathcal{H}_{t}$. Indeed, delaying the demand leads to
inconvenience for the follower, thus powering the device as early as possible
minimizes the inconvenience.

The instances $P_{\mathcal{H}}$ and $P^{t}_{s_{RH}}$ for all $t\in\mathcal{H}$
all involve a scenario tree. Here, we assume that the scenario tree of
$P_{\mathcal{H}}$ is a complete tree with time periods of length $s_{RH}$.
Therefore, one of the base scenarios is selected at every iteration of step 6:
the scenario selected at time $t$ is the scenario that comes to reality for
the time slots $t,\dots,t+s_{RH}-1$. However, considering a complete tree on
$\mathcal{H}_{t}$ can be difficult, even if $\mathcal{H}_{t}$ is small, as
shown in Section 5. Therefore, the scenario tree of $P^{t}_{s_{RH}}$ is
assumed to only comprise the base scenarios. Furthermore, we assume that the
probability to choose one of the scenarios only depends on the scenario that
was selected at the previous iteration, as in a Markov process.

In step 5, the problem solved is a MIP, as defined in Section 3.1, to which we
add constraints ensuring that the prices belonging to the frozen horizon are
indeed the prices that have been determined in the previous iterations.

In step 7, the added information consists in the choice made for the scenario
in step 6, the follower’s decisions associated to the chosen scenario for the
time slots $t$ to $t+s_{RH}-1$ and the leader’s decision for the frozen
horizon (except for the last iteration of the algorithm, where the leader’s
and follower’s decisions are saved until the end of the horizon). At the end
of the rolling horizon method, the data thus contains one scenario and the
associated decisions for the whole time horizon $\mathcal{H}$.

##  5 Numerical Results

This section is divided into two parts. First, a sensitivity analysis is
conducted on a one-week instance with a single DG scenario to study the
influence of individual parameters on the results. Second, we evaluate the
efficiency of the rolling horizon method, focusing on how the length of the
frozen horizon affects the obtained results.

All the results presented here were obtained with CPLEX 12.6 on a Linux
virtual machine with 10 Go RAM working on a computer equipped with an Intel
i7-4600u processor at 2.10 GHz.

###  5.1 Test instances

All the numerical tests are conducted on an instance having the following
features:

  * •

Each time slot represents a 30-minute period,

  * •

The instance comprises 336 time slots, adding up to one full week,

  * •

The energy costs are based on the prices of electricity on the spot market
during a fall week (see Figure 2),

  * •

There are 120 devices. Their demand in the ideal case (i.e., without delay)
approximates the electricity consumption due to the heating of a dozen
households during a fall week (see Figure 3),

  * •

The DG scenario is based on the load factor profile for a fall week, with a
generation capacity large enough to cover the electricity consumption during
off-peak periods (see Figure 4).

The data used to create the test instance were supplied by the industrial
partner (EDF).

Figure 2: The spot prices retrieved from industrial data. Figure 3: The energy
consumption retrieved from industrial data, and implied by the devices if they
are used during the first time slots of their associated time window. Figure
4: The original DG scenario.

###  5.2 Sensitivity analysis

The parameters that are perturbed to conduct the sensitivity analysis are the
DG scenario, the inconvenience coefficients, the energy costs, the battery
size, and the time window sizes. In total, thirteen instances are solved. From
these thirteen instances, only one is solved to optimality within the time
limit of $1200$ seconds, namely the instance where the inconvenience
coefficients are all set to $0$. For the other twelve, the optimality gaps
range from $0.3\%$ to $1.9\%$. Usually, the solver can reach a relatively good
solution early in the process, but fails to close the optimality gap
subsequently.

The numbers of variables and constraints of the MIP not only depend on the
number of time slots and the number of devices, but also on the length of the
time windows and their position in the time horizon. A time window located in
the first time period will generate more nonanticipativity constraints than a
time window located near the end of the time horizon. As an indication,
solving the MIP relative to the original test instance requires the use of
35530 constraints, 33745 continuous variables and 1157 binary variables.

To perform the sensitivity analysis and quantify the impact of the bilevel
model, we compute the relevant values of the so-called _reference case_ , in
opposition with the _optimized case_ , which consists in the optimal solution
of (SBPP). The reference case is defined as follows: first, the leader sets
their prices at the same level as the competitor. Then the follower determines
the schedule of the devices in the most convenient way: each device $(n,a)$ is
powered at its maximum power $\beta^{max}_{(n,a)}$ during the first hours of
the related time slots until $E_{(n,a)}$ is reached, then the device is
stopped. To power the device, the follower makes DG their first choice, then
energy from the battery, and finally energy from the leader. If at a given
time slot, the DG exceeds the demand, the remaining energy is stored, unless
the battery is full. The energy demand of the follower in the reference case
is depicted in Figure 3. In particular, observe that the reference case
minimizes the inconvenience for the follower, but does not optimize the DG
use, which implies that in the reference case, the follower’s response is not
optimal with relation to the leader’s prices.

The results for the objective values of the leader are given in Table 1,
whereas Table 2 gathers information about the follower’s objective. In the
first column of both tables, "bat." stands for battery, "inc." for
inconvenience, "Spot" for the instance with higher market prices and "TW" for
the instances where the sizes of the time windows vary. "BC ref" represents
the billing cost of the follower in the reference case, "IC ref" the
inconvenience cost in the reference case, "BC opt" the billing cost when the
prices of the leader are optimized, and "IC opt" the inconvenience cost when
the prices of the leader are optimized. Finally, % BC, % IC and % GC show the
difference (in percentage) between the reference case and the optimized case
of the billing cost, the inconvenience cost and the generalized cost
respectively (e.g., % BC $=100\cdot$BC opt$/$BC ref).

Instance | Ref. obj. | Opt. obj. | % diff.
---|---|---|---
Base | $34172.68$ | $34676.53$ | $1.47$
Zero bat. | $34261.4$ | $34372.9$ | $0.33$
Small bat. | $34201.32$ | $34521.54$ | $0.94$
Large bat. | $34182.14$ | $34714.16$ | $1.56$
Zero inc. | $34172.68$ | $35753.54$ | $4.63$
Low inc. | $34172.68$ | $35385.64$ | $3.55$
High inc. | $34172.68$ | $34389.21$ | $0.63$
Zero DG | $46597.9$ | $47101.38$ | $1.08$
Low DG | $40390.54$ | $40721.06$ | $0.82$
High DG | $29084.13$ | $28611.95$ | $-1.62$
Spot | $33240.14$ | $34323.98$ | $3.26$
Small TW | $34172.68$ | $34668.36$ | $1.45$
Large TW | $34172.68$ | $34531.8$ | $1.05$
Table 1: The leader’s objective values in the reference case, in the optimized case, and the difference in percentages. Instance | BC ref | IC ref | BC opt | IC opt | % BC | % IC | % GC
---|---|---|---|---|---|---|---
Base | $46573.2$ | $525.13$ | $46025.7$ | $877.31$ | $98.82$ | $167.07$ | $99.59$
Zero bat. | $46659.2$ | $525.13$ | $45964.7$ | $856.52$ | $98.51$ | $163.11$ | $99.23$
Small bat. | $46589.2$ | $525.13$ | $46246.8$ | $737.52$ | $99.27$ | $140.45$ | $99.72$
Large bat. | $46573.2$ | $525.13$ | $46379.4$ | $646.92$ | $99.58$ | $123.19$ | $99.84$
Zero inc. | $46573.2$ | $0$ | $46573.2$ | $0$ | $100$ | $-$ | $100$
Low inc. | $46573.2$ | $262.56$ | $46343.7$ | $540.29$ | $99.51$ | $205.77$ | $100.1$
High inc. | $46573.2$ | $787.69$ | $46052.1$ | $1047.97$ | $98.88$ | $133.04$ | $99.45$
Zero DG | $63882$ | $525.13$ | $63737.8$ | $655.5$ | $99.77$ | $124.83$ | $99.98$
Low DG | $55227.6$ | $525.13$ | $54784.7$ | $712.29$ | $99.2$ | $135.64$ | $99.54$
High DG | $39317.2$ | $525.13$ | $37755.5$ | $786.96$ | $96.03$ | $149.86$ | $96.74$
Spot | $46573.2$ | $525.13$ | $46002.4$ | $905.36$ | $98.77$ | $172.41$ | $99.6$
Small TW | $46573.2$ | $525.13$ | $46383.4$ | $681.32$ | $99.59$ | $129.75$ | $99.93$
Large TW | $46573.2$ | $525.13$ | $45543.8$ | $1031.83$ | $97.79$ | $196.49$ | $98.89$
Table 2: The follower’s detailed objective values.

####  5.2.1 DG scenarios

Besides the original test instance, three more instances are solved with DG
scenarios equivalent to respectively $0$, $0.5$, and $1.5$ times the original
DG scenario which is represented in Figure 4.

As much as possible of the energy coming from DG is used (that is, directly
consumed or stored for later consumption), as this energy is considered free
by the SGO. This behavior is illustrated in Figure 5. It follows directly that
if the DG increases, the amount of energy bought from the grid will decrease
accordingly. This naturally widely affects the leader’s profits, as Table 1
shows: those profits range from $28611$ when the DG is high to $47101$ when
the DG is nonexistent.

Figure 5: The DG consumption for the base instance, for time slots $54$ to
$96$.

As for the demand, the fact that the follower’s problem is linear induces an
"all or nothing" situation: if for time slots $h_{1}<h_{2}$, the generalized
cost of powering a device $(n,a)$ during $h_{2}$ is smaller than during
$h_{1}$ (i.e.,
$\min\\{p^{h_{1}},\bar{p}^{h_{1}}\\}+C^{h_{1}}_{(n,a)}>\min\\{p^{h_{2}},\bar{p}^{h_{2}}\\}+C^{h_{2}}_{(n,a)}$),
then all the energy that can possibly be shifted from $h_{1}$ to $h_{2}$ will
be. This explains the huge variations in the follower’s demand (e.g., from $0$
during time slot $23$ to $30$ during the next time slot) that can be observed
in Figure 6. Furthermore, it can be observed in the same figure that power
coming from the DG is entirely used: smooth curves following the DG curve can
be observed between time $16$ and $29$, especially on the curve illustrating
the demand with a high DG scenario.

Figure 6: The follower’s demands for the instances with various DG scenarios
for the time slots $15$ to $40$ (bottom).

####  5.2.2 Inconvenience coefficients

To study the influence of the inconvenience coefficients on the results, we
consider four instances. Recall beforehand that the inconvenience coefficients
increase linearly with the time for each device, therefore the difference
between the inconvenience coefficients of two consecutive hours is constant,
and assumed to be the same for all devices: in the original case (normal
inconvenience), this value is $0.0625$. In the three other cases, it is
respectively $0$, $0.5$, and $1.5$ times the original value, that is to say
$0$ in the zero inconvenience case, $0.03125$ in the low inconvenience case,
and $0.09375$ in the high inconvenience case. Unsurprisingly, the leader’s
profit decreases as the inconvenience coefficients grow (see Table 1): if the
consumers are more willing to shift their loads, less incentive is necessary
to induce the same shift, and thus more profit is achievable for the leader.

The latter fact can be easily observed in Figure 7. The optimal leader price
profile when there is no inconvenience is constant at the level of the
competitor’s prices, because no incentive is required to induce a load shift.
The optimistic assumption implies that the follower’s answer is precisely what
is best for the leader. At the contrary, inducing load shifts for high
inconvenience values is difficult. Prices resulting in a load shift despite a
high inconvenience are clearly lower than prices achieving a similar load
shift with a low inconvenience. This is conspicuous for example on time slot
$33$. For other time slots, the inconvenience is too high to induce shift
(during time slot $7$ for instance) which is why the producer lowers their
prices during these periods only for the instances with low inconvenience
values.

Finally, let us observe that the prices tend to decrease linearly (with the
slope depending on the inconvenience) as the supplier induces delays. This can
be observed in Figure 7: for example, the optimal prices for the instance with
high inconvenience linearly decrease between time slots $11$ and $19$, with
the noticeable exceptions constituted by time slots $13$ and $18$. The
explanation lies in the optimistic assumption of the problem: when the
follower gets their energy at the same generalized cost at two time slots,
they consumes during the time slot that is most advantageous for the leader.
When the prices decrease linearly over a given period, they compensate the
inconvenience costs that grow linearly with the delay, so that the generalized
cost remains the same over the period.

Figure 7: The leader price profiles for the instances with various
inconvenience coefficients for the $40$ first time slots.

####  5.2.3 Energy costs

In this paper, the energy supplied by the leader is assumed to be purchased on
the spot market at prices that are known in advance. The spot market prices
thus have a large influence, as they determine which time slots are the most
profitable for the leader, or at the contrary which time slots are the least
valuable. To show this influence, we consider an instance where the spot
market prices are $20\%$ higher during peak periods, i.e., periods when the
prices are particularly high (more than $4.5$). In that case, the leader has
an increased motivation to induce a load shift from peak to off-peak periods,
as such a shift induces a larger profit.

Like for the two parameters previously studied, it is not surprising that the
leader’s optimal value is lower when the spot prices are higher: with the
original prices, the leader’s objective reaches $34677$, whereas it only
attains $34324$ with the higher peak prices. Observe, however, that the
difference is not as significant as when the DG scenario varies: whereas
differences in the leader’s optimal value can reach several thousand, it is
limited to $353$ in the case of higher market prices. This value is in the
same order of magnitude as the differences between the instance with various
inconvenience coefficients ($939$ between the instance with low inconvenience
and the instance with high inconvenience).

Figure 8: The follower’s demands and the leader price profiles for the
instances with various spot market prices.

In Figure 8, observe that the demand curves differ for the high market prices
instance and the base instance: the differences occur almost always during the
peak periods, i.e., when the market prices in both considered instances are
different. The higher need to shift the load from peak to off-peak periods
becomes clear while having a look at the optimal price profiles. The optimal
prices for the high market prices instance are indeed generally lower than the
optimal prices of the base instance, e.g., during time slots $25$ and $170$.

####  5.2.4 Battery Size

Besides the original instance, three instances with various battery sizes are
considered: a zero storage capacity, a small battery, and a large battery,
which are equal to respectively $0$, $0.5$, and $1.5$ times the size of the
original battery. Clearly, storage capacities are filled when electricity
prices are low, and the stored energy is consumed when prices are higher. One
might think that a greater storage capacity represents an advantage for the
follower. However, it turns out that a greater capacity is especially
advantageous for the leader, as the leader’s objective values indicate in
Table 1. The reason is that a larger battery allows for more freedom.
Specifically, the follower can purchase more energy during cheap time slots,
when the leader makes the most profit, and store this energy for expensive
time slots.

In Figure 9, the battery states are shown for the four considered storage
sizes. The origin of the stored energy is indicated for the instance with the
largest storage capacity. Observe that most of the stored energy is bought
from the grid, instead of being taken from the DG. As for the demand curves,
the "all or nothing" nature of linear programming is clearly visible, since
the battery often oscillates between being full and being empty. These
oscillations are particularly visible for the small battery instance during
time slots $97$ to $109$ (Figure 9, bottom graph).

Figure 9: The battery states for the instances with various storage sizes and
the origin of the stored energy for the large storage instance, for time slots
$97$ to $149$.

####  5.2.5 Time Window Sizes

Finally, the last parameters that we considered are the sizes of the time
windows. Whereas the length of the time windows varies between $1.8$ and $2.5$
times the number of time slots needed to power each appliance in the base
instance, these factors are in a range from $1.4$ to $2$ in the instance with
narrower time windows, and between $2.14$ and $3$ in the instance with longer
time windows. Observe, however, that the difference of profit for the leader
is rather small ($34677$ with normal time windows, $34668$ with narrow time
windows and $34532$ for longer time windows), and thus that this parameter
does not much influence the leader’s profit, nor the follower’s results.

####  5.2.6 A Word on the Objective Values

The solutions of (SBPP) and the reference case are compared on the thirteen
instances considered in this section. The objective values of the leader are
given in Table 1. Except for the instance with high DG, the leader’s profit is
always higher in the optimized case than in the reference case. Although the
difference in profit is rather small (between $0.33\%$ and $4.63\%$), it is
not negligible. Unsurprisingly, the largest difference between the profits in
the reference case and in the optimized case lies in the instance where the
follower does not face inconvenience. As argued in Section 5.2.2, no incentive
is required to make the follower react in the most favorable way for the
leader. Therefore, $4.63\%$ is the best possible increase for the leader’s
profit, in the case where the inconvenience factors vary. Concerning the
instance with higher peak energy costs, the difference in profit is
significant as well, with a difference equal to $3.26\%$. The energy cost
differences indeed make load shifts more advantageous for the leader than in
the instance with normal inconvenience. The example with high DG resulting in
a lower leader profit can easily be explained by the fact that the DG use in
the reference case is not optimal, which leads to significant differences as
the DG quantities become important.

On the follower’s side, even though (SBPP) is not a zero-sum game, there is a
clear relation between the follower’s and the leader’s objectives. The billing
cost of the follower constitutes the revenues of the leader, as no energy is
purchased from the competitor due to the optimistic assumption and the fact
that the competitor prices are greater than the energy costs. Therefore,
intuitively, the leader’s optimized prices are going to induce a follower’s
reaction that will bring approximately the same objective value as in the
reference case. Of course, the follower’s optimal objective value in the
reference case constitutes an upper bound for the follower’s optimal value
associated with any leader price profile, since the consumption schedule of
the reference case is always feasible. The results for the follower’s optimal
value are illustrated in Table 2.

All follower’s generalized costs are smaller in the optimized case than in the
reference case, except for the low inconvenience instance, which is probably
due to the commercial solver finding a local optimum instead of a global one.
Furthermore, the percentages are all very close to $100\%$, the furthest being
the instance with large time windows ($98.89\%$), which confirms the above-
mentioned intuition. Furthermore, the billing cost is always smaller in the
optimized case, but the inconvenience cost increases. This is due to the fact
that in the reference case, the follower’s inconvenience is actually
minimized.

Note that in the reference case, the parameters influencing the billing cost
are the size of the battery and the DG, whereas only changes in the
inconvenience coefficients induce a change in the inconvenience cost.

###  5.3 Test of the RH Method

To evaluate the efficiency of the rolling horizon methods, a set of tests has
been determined. The parameter we choose to analyze is the length of the
frozen horizon $l_{FH}$. The considered instance is the one described in
Section 5.1 with the three nonzero scenarios used in Section 5.2.1 as base
scenarios. For the various parameters, it has been decided that the length of
the rolling horizon would be $12$, whereas the iteration step would be $1$.
The values tested for the length of the frozen horizon are all even numbers
between (and including) $0$ and $10$. The number of possible scenarios over
the whole time horizon being obviously gigantic, a set of five random
scenarios have been generated as follows: during the run with $l_{FH}=0$, at
each iteration of the rolling horizon algorithm, the scenario that gets real
between times $t$ and $t+1$ is chosen following a Markov process, the
probability to keep the same scenario as between $t-1$ and $t$ being equal to
$0.4$, and the probability to switch to another scenario being equal to $0.3$.
In the runs for longer frozen horizons, the choice of the scenario at each
iteration is forced as the scenario that should be followed has been
previously determined.

To evaluate the efficiency of the rolling horizon method, two comparison
values are computed. Ideally, these values should represent bounds, but in
fact, they are not.

  * •

The first one is the _reference case_ , as presented in Section 5.2.6,

  * •

The second one is the _perfect case_ , that is the case where the scenario is
known in advance, and both the leader and the follower completely optimize
their decisions.

Observe that the second comparison value is strongly inspired by the well-
known expected value of perfect information (see e.g. Birge and Louveaux
(2011)): the value computed here would be the _wait-and-see_ solution.
However, the bilevel structure of the problem prevents this value to be a
bound, even without considering the effects of the rolling horizon method: the
information is shared both at the leader’s and the follower’s level, which
might benefit more to the follower than to the leader (see von Niederhäusern
(2019) for more details).

In Figure 10, the leader’s objective values are presented. In the above graph,
observe that no pattern is followed as the length of the frozen horizon grows:
for the first scenario, the leader’s objective value is higher with $l_{FH}=8$
than with $l_{FH}=6$, whereas for the third run, the opposite occurs. We first
expected the leader’s profit to be higher with shorter frozen horizons, as the
leader has more freedom, but it seems this is not the case in general. A
potential reason for this is that many of the problems that are considered
during the process are not solved to optimality and only feasible solutions
are provided, as the solver only gets 150 seconds to find a solution at each
iteration.

Figure 10: Leader’s objective values
The above graph indicates only the leader values, whereas the chart below
shows the comparison with the two computed comparison values.

Observe, furthermore, that first, the leader’s objective value with the
rolling horizon method is always higher than the objective value in the
reference case, which proves that it is worthwhile trying to optimize the
prices through DSM. Second, the objective value in the perfect case is also
surpassed by the leader’s objective value obtained with the rolling horizon
method, except for the second scenario with $l_{FH}=6$. Obviously, the
leader’s price profile obtained with the rolling horizon method is feasible
when computing the perfect case; thus one might think that the perfect case
would bring a better profit to the leader. However, the stochastic nature of
the follower’s problem combined with the functioning of the rolling horizon
method implies that the follower’s response to the leader’s prices might not
be optimal. In our experiments, this happens regularly: often, the follower’s
response with the rolling horizon method even performs worse than the
follower’s response in the reference case, as can be seen in Figure 11.

Figure 11: Follower’s objective values
Comparison for the follower’s objective values for all runs of the rolling
horizon methods and for the reference case.

Let us recall that the follower’s objective consists in the purchase costs
plus the inconvenience costs. The latter only amount to a small fraction of
the generalized cost: they range from $525.13$ in the reference case (in that
case, the inconvenience is minimal) to $878.13$ among generalized costs
ranging from $45830$ to $50597$, which means less than two percent.
Furthermore, it is more than logical to have higher inconvenience costs with
the rolling horizon method, since the follower’s response in the reference
case actually minimizes the inconvenience.

Besides the nature of the rolling horizon method that avoids any optimality
guarantee of the follower’s response, the fact that each step of the rolling
horizon method has to be solved in a time limit of $150$ seconds even hinders
CPLEX to find optimal solutions for all the subproblems encountered during the
algorithm’s run. Several artifacts confirm this. First, sometimes, energy is
bought from the competitor, which should not happen: at the price offered by
the competitor, the leader makes profit. Therefore, buying energy from the
competitor contradicts the optimistic assumption. This case, however, does not
often happen: among the $30$ runs of the algorithm, $17$ have strictly
positive values for $\sum_{n\in N,a\in A_{n},h\in
T_{(n,a)}}\bar{x}^{h}_{(n,a)}$, but only seven have values exceeding $1$, for
a total demand of $4563$, with a maximum at $27.25$ (see Figure 12 for their
repartition), which represents less than $381.5$ in the generalized cost.
Although not negligible, the importance of these flaws is relatively small.
Observe, furthermore, that problems occur mostly in the runs with smaller
lengths of frozen horizons. The justification is easy: the smaller the frozen
horizon, the more leader variables, and thus the higher is the difficulty to
obtain an optimal solution to each of the bilevel subproblems.

Figure 12: Energy quantity bought from the competitor
The bars represent the sum over all hours of the energy quantities bought from
the competitor.

The other obviously nonoptimal behavior concerns the usage of renewable
energy. In some of the algorithm’s runs, the renewable energy, that is
considered as free by the follower, is not used completely. Besides, some
cases are more problematic, as the follower consumed more renewable energy
than what is available (see Figure 13). This means that the equation that
bounds the usage of renewable energy is not always satisfied. In fact, this
happens on average between $15$ and $16$ times per run. Observe, however, that
when $l_{FH}=0$, the above-mentioned equation is always satisfied. Moreover,
the over- and underconsumption of renewable energy represent small percentages
of the total disposable energy: the values are bounded by $50$ in most cases,
which means less than $4$ percent. To justify those computing mistakes, we
assume that again, the time left to the solver at each iteration is too small.

Figure 13: Unused renewable energy
The bars represent the difference between the renewable energy that is
produced and the one that is either consumed or stored.

Despite the previously mentioned flaws, using the rolling horizon method is
worth it. Besides the obvious fact that time goes on, implying a constant need
for reoptimization, the method allows for leader profits that exceed the
losses (with relation to the reference case) of the follower. This situation
is illustrated in Figure 14. In a realistic setting, the possibility of
refunding the follower for the losses should be contemplated, in order to
motivate the clients to play along. No rational client would take the risk to
pay more (which is often the case in our setting) to help their supplier to
make more profit. However, since the leader’s profit easily covers the
follower’s losses, the deal would be harmless for the leader.

Figure 14: The advantages of the rolling horizon method
Above, the difference of profit for the leader with the rolling horizon method
and in the reference case.
Below, the difference in the objective function of the follower between the
rolling horizon method and the reference case, which can be seen as a loss.

Let us finally mention the computing time to run the rolling horizon method.
Unfortunately, it takes a very long time: between $5666$ and $11540$ seconds.
Observe that with a time limit at each iteration of $150$ seconds, the
theoretical maximum time for a run of the method is around $50000$.
Furthermore, as a rule of thumb, the needed time decreases as the length of
the frozen horizon increases: a greater $l_{FH}$ means a smaller number of
leader variables, and thus easier problems to solve at each iteration.
Although the running times may seem large, they actually allow considering the
horizon of a whole week, with a tremendous number of possible scenarios that
would be totally impossible to handle otherwise.

Figure 15: The running times for the various runs of the rolling horizon
method.

##  6 Conclusion

In this paper, we have studied a bilevel pricing problem for demand-side
management with a strong stochastic component under the form of scenario
trees. Besides proving the validity of the pricing approach through numerical
analyses, we designed a rolling horizon method that numerically proved to be
applicable in real situations, for which the original problem cannot be solved
due to its exponential size.

Being a first attempt at applying rolling horizons in a bilevel framework,
prospects remain numerous. Specifically, using efficient heuristics at every
iteration of the rolling horizon method could help reach better solutions, and
possibly avoid feasibility problems when the time allowed at each iteration is
not enough for the MILP solver to find a feasible point. From a problem design
point of view, it could also be interesting to consider robust optimization to
solve (SBPP), in order to maximize the worst-case revenue of the leader.

## Acknowledgment

This research benefited from the support of the FMJH Program Gaspard Monge in
Optimization and Operations Research, and from the support to this program
from EDF.

##  Appendix A Dual constraints of
$\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right)$

The dual constraints of
$\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right)$ are as follows:

$\displaystyle\displaystyle
d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)dx^{\sigma,\sigma^{\prime},h}_{(n,a)}\leq
P[\sigma]\left(p^{h}+C^{h}_{(n,a)}\right)$ $\displaystyle\forall n\in N,a\in
A_{n},h\in T_{(n,a)},\sigma\in\Sigma$ $\displaystyle\displaystyle
d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)d\bar{x}^{\sigma,\sigma^{\prime},h}_{(n,a)}\leq
P[\sigma]\left(\bar{p}^{h}+C^{h}_{(n,a)}\right)$ $\displaystyle\forall n\in
N,a\in A_{n},h\in T_{(n,a)},\sigma\in\Sigma$ $\displaystyle\displaystyle
d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}-d7^{h,\sigma}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)d\lambda^{\sigma,\sigma^{\prime},h}_{(n,a)}\leq
P[\sigma]C^{h}_{(n,a)}$ $\displaystyle\forall n\in N,a\in A_{n},h\in
T_{(n,a)},\sigma\in\Sigma$ $\displaystyle\displaystyle
d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}+dS-ds^{h,\sigma}_{(n,a)}$
$\displaystyle\displaystyle\qquad+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)ds^{\sigma,\sigma^{\prime},h}_{(n,a)}\leq
P[\sigma]C^{h}_{(n,a)}$ $\displaystyle\forall n\in N,a\in A_{n},h\in
T_{(n,a)},\sigma\in\Sigma$
$\displaystyle\displaystyle-\rho^{c}d4^{h,\sigma}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)dx^{\sigma,\sigma^{\prime},h}_{s}\leq
P[\sigma]p^{h}$ $\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$
$\displaystyle\displaystyle-\rho^{c}d4^{h,\sigma}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)d\bar{x}^{\sigma,\sigma^{\prime},h}_{s}\leq
P[\sigma]\bar{p}^{h}$ $\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$
$\displaystyle\displaystyle-d7^{h,\sigma}-\rho^{c}d4^{h,\sigma}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)d\lambda^{\sigma,\sigma^{\prime},h}_{s}\leq
0$ $\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$
$\displaystyle\displaystyle
d3^{\sigma}_{0}-\rho^{d}d4^{h,\sigma}+ds^{h,\sigma}_{(n,a)}+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)dS^{\sigma,\sigma^{\prime},h}\leq
0$ $\displaystyle\forall\sigma\in\Sigma,h=0$ $\displaystyle\displaystyle
d4^{h-1,\sigma}-\rho^{d}d4^{h,\sigma}+ds^{h,\sigma}_{(n,a)}+d5^{h,\sigma}_{min}-d5^{h,\sigma}_{max}$
$\displaystyle\displaystyle\qquad+\sum_{\sigma^{\prime}\in\mathrm{Com_{h}\left(\sigma\right)}}sgn\left(\sigma,\sigma^{\prime}\right)dS^{\sigma,\sigma^{\prime},h}\leq
0$ $\displaystyle\forall h\in\mathcal{H}\setminus\\{0\\},\sigma\in\Sigma$
$\displaystyle\displaystyle
d4^{h-1,\sigma}+d5^{h,\sigma}_{min}-d5^{h,\sigma}_{max}\leq 0$
$\displaystyle\forall\sigma\in\Sigma,h=|\mathcal{H}|+1$
$\displaystyle\displaystyle
d1^{\sigma}_{(n,a)},d2^{\sigma}_{(n,a)},d7^{h,\sigma},ds^{h,\sigma}_{(n,a)},d5^{h,\sigma}_{min},d5^{h,\sigma}_{max}\geq
0$ $\displaystyle\forall\sigma\in\Sigma,h\in\mathcal{H},n\in N,a\in A_{n},$

where

$\mathrm{Com_{h}\left(\sigma\right)}=\left\\{\sigma^{\prime}\in\Sigma\mid\left(\sigma,\sigma^{\prime}\right)\in\mathrm{Com}_{h}\mathrm{\
or\ }\left(\sigma^{\prime},\sigma\right)\in\mathrm{Com}_{h}\right\\},$

and the sign function is defined as follows:

$sgn\left(\sigma,\sigma^{\prime}\right)=\left\\{\begin{array}[]{l
l}1&\mathrm{if\ }\sigma\prec\sigma^{\prime}\\\ 0&\mathrm{if\
}\sigma=\sigma^{\prime}\\\ -1&\mathrm{if\ }\sigma\succ\sigma^{\prime}.\\\
\end{array}\right.$

##  Appendix B Complementarity constraints of
$\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right)$

The complementarity constraints of
$\left(P^{\operatorname{sto}}_{\operatorname{SGO}}\right)$ are as follows:

$\displaystyle
x^{h}_{(n,a)}\left(d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}-P[\sigma]\left(p^{h}+C^{h}_{(n,a)}\right)\right)=0$
$\displaystyle\forall n\in N,a\in A_{n},h\in T_{(n,a)},\sigma\in\Sigma$
$\displaystyle\bar{x}^{h}_{(n,a)}\left(d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}-P[\sigma]\left(\bar{p}^{h}+C^{h}_{(n,a)}\right)\right)=0$
$\displaystyle\forall n\in N,a\in A_{n},h\in T_{(n,a)},\sigma\in\Sigma$
$\displaystyle\lambda^{h}_{(n,a)}\left(d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}-d7^{h,\sigma}-P[\sigma]C^{h}_{(n,a)}\right)=0$
$\displaystyle\forall n\in N,a\in A_{n},h\in T_{(n,a)},\sigma\in\Sigma$
$\displaystyle
s^{h}_{(n,a)}\left(d1^{\sigma}_{(n,a)}-d2^{\sigma}_{(n,a)}+dS-d6^{h,\sigma}_{max}-P[\sigma]C^{h}_{(n,a)}\right)=0$
$\displaystyle\forall n\in N,a\in A_{n},h\in T_{(n,a)},\sigma\in\Sigma$
$\displaystyle x^{h}_{s}\left(-\rho^{c}d4^{h,\sigma}-P[\sigma]p^{h}\right)=0$
$\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$
$\displaystyle\bar{x}^{h}_{s}\left(-\rho^{c}d4^{h,\sigma}-P[\sigma]\bar{p}^{h}\right)=0$
$\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$
$\displaystyle\lambda^{h}_{s}\left(-d7^{h,\sigma}-\rho^{c}d4^{h,\sigma}\right)=0$
$\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$ $\displaystyle
S^{h}\left(d3^{\sigma}_{0}-\rho^{d}d4^{h,\sigma}+d6^{h,\sigma}_{max}\right)=0$
$\displaystyle\forall\sigma\in\Sigma,h=0$ $\displaystyle
S^{h}\left(d4^{h-1,\sigma}-\rho^{d}d4^{h,\sigma}+d6^{h,\sigma}_{max}+d5^{h,\sigma}_{min}-d5^{h,\sigma}_{max}\right)=0$
$\displaystyle\forall h\in\mathcal{H}\setminus\\{0\\},\sigma\in\Sigma$
$\displaystyle
S^{h}\left(d4^{h-1,\sigma}+d5^{h,\sigma}_{min}-d5^{h,\sigma}_{max}\right)=0$
$\displaystyle\forall\sigma\in\Sigma,h=|\mathcal{H}|+1$ $\displaystyle
d1^{\sigma}_{(n,a)}\left(\sum_{h\in
T_{(n,a)}}\left(x^{h,\sigma}_{(n,a)}+\bar{x}^{h,\sigma}_{(n,a)}+\lambda^{h,\sigma}_{(n,a)}+s^{h,\sigma}_{(n,a)}\right)-E_{(n,a)}\right)=0$
$\displaystyle\forall n\in N,a\in A_{n},\sigma\in\Sigma$ $\displaystyle
d2^{\sigma}_{(n,a)}\left(x^{h,\sigma}_{(n,a)}+\bar{x}^{h,\sigma}_{(n,a)}+\lambda^{h,\sigma}_{(n,a)}+s^{h,\sigma}_{(n,a)}-\beta^{max}_{(n,a)}\right)=0$
$\displaystyle\forall n\in N,a\in A_{n},h\in T_{(n,a)},\sigma\in\Sigma$
$\displaystyle d7^{h,\sigma}\left(\lambda^{h,\sigma}_{s}+\sum_{n\in
N}\sum_{a\in
A_{n}}\lambda^{h,\sigma}_{(n,a)}-\lambda^{h,\sigma}_{\max}\right)=0$
$\displaystyle\forall h\in\mathcal{H},\sigma\in\Sigma$ $\displaystyle
d6^{h,\sigma}_{max}\left(\sum_{n\in N}\sum_{a\in
A_{n}}s^{h,\sigma}_{(n,a)}-S^{h,\sigma}\right)=0$ $\displaystyle\forall
h\in\mathcal{H},\sigma\in\Sigma$ $\displaystyle
d5^{h,\sigma}_{min}\left(S^{\min}-S^{h,\sigma}\right)=0$ $\displaystyle\forall
h\in\mathcal{H},\sigma\in\Sigma$ $\displaystyle
d5^{h,\sigma}_{max}\left(S^{h,\sigma}-S^{\max}\right)=0$ $\displaystyle\forall
h\in\mathcal{H},\sigma\in\Sigma.$

## References

  * Kreith and Goswami [2016] Frank Kreith and D. Yogi Goswami.  _Energy Management and Conservation Handbook_.  CRC Press, October 2016.  ISBN 978-1-4665-8517-1.
  * Wang et al. [2016] Yunpeng Wang, Walid Saad, Narayan B. Mandayam, and H. Vincent Poor.  Load Shifting in the Smart Grid: To Participate or Not?  _IEEE Transactions on Smart Grid_ , 7(6):2604–2614, November 2016.  ISSN 1949-3053, 1949-3061.  doi: 10.1109/TSG.2015.2483522.
  * Farhangi [2010] H. Farhangi.  The path of the smart grid.  _IEEE Power and Energy Magazine_ , 8(1):18–28, January 2010.  ISSN 1540-7977.  doi: 10.1109/MPE.2009.934876.
  * Kabalci [2016] Yasin Kabalci.  A survey on smart metering and smart grid communication.  _Renewable and Sustainable Energy Reviews_ , 57:302–318, May 2016.  ISSN 1364-0321.  doi: 10.1016/j.rser.2015.12.114.
  * Alekseeva et al. [2018] Ekaterina Alekseeva, Luce Brotcorne, Sébastien Lepaul, and Alain Montmeat.  A bilevel approach to optimize electricity prices.  _Yugoslav Journal of Operations Research_ , (00):2–2, 2018.
  * Afşar et al. [2016] Sezin Afşar, Luce Brotcorne, Patrice Marcotte, and Gilles Savard.  Achieving an optimal trade-off between revenue and energy peak within a smart grid environment.  _Renewable Energy_ , 91(Supplement C):293–301, June 2016.  ISSN 0960-1481.  doi: 10.1016/j.renene.2016.01.055.
  * Aussel et al. [2020] Didier Aussel, Luce Brotcorne, Sébastien Lepaul, and Léonard von Niederhäusern.  A trilevel model for best response in energy demand-side management.  _European Journal of Operational Research_ , 281(2):299–315, 2020.  ISSN 0377-2217.  doi: 10.1016/j.ejor.2019.03.005.
  * Park and van der Schaar [2010] Jaeok Park and Mihaela van der Schaar.  Pricing and Incentives in Peer-to-Peer Networks.  In _2010 Proceedings IEEE INFOCOM_ , pages 1–9, March 2010.  doi: 10.1109/INFCOM.2010.5461952.
  * Holguín-Veras and Aros-Vera [2015] José Holguín-Veras and Felipe Aros-Vera.  Self-supported freight demand management: Pricing and incentives.  _EURO Journal on Transportation and Logistics_ , 4(2):237–260, June 2015.  ISSN 2192-4384.  doi: 10.1007/s13676-013-0041-1.
  * Ozdaglar and Srikant [2007] Asuman Ozdaglar and Rayadurgam Srikant.  Incentives and pricing in communications networks.  _Algorithmic Game Theory_ , pages 571–592, January 2007.  doi: 10.1017/CBO9780511800481.024.
  * Tsaousoglou et al. [2019] Georgios Tsaousoglou, Nikolaos Efthymiopoulos, Prodrommos Makris, and Emmanouel Varv Arigos.  Personalized real time pricing for efficient and fair demand response in energy cooperatives and highly competitive flexibility markets.  _Journal of Modern Power Systems and Clean Energy_ , 7(1):151–162, January 2019.  ISSN 2196-5420.  doi: 10.1007/s40565-018-0426-0.
  * Liu et al. [2019] Yi Liu, Liye Xiao, Guodong Yao, and Siqi Bu.  Pricing-Based Demand Response for a Smart Home With Various Types of Household Appliances Considering Customer Satisfaction.  _IEEE Access_ , 7:86463–86472, 2019.  ISSN 2169-3536.  doi: 10.1109/ACCESS.2019.2924110.
  * von Stackelberg [1934] Heinrich von Stackelberg.  _Marktform Und Gleichgewicht_.  J. Springer, 1934.
  * Bard [2010] Jonathan F. Bard.  _Practical Bilevel Optimization: Algorithms and Applications_.  Springer Publishing Company, Incorporated, 1st edition, 2010.  ISBN 978-1-4419-4807-6.
  * Dempe [2002] S. Dempe.  _Foundations of Bilevel Programming_ , volume 61 of _Nonconvex Optimization and Its Applications_.  Kluwer Academic Publishers, Boston, 2002.  ISBN 978-1-4020-0631-9.  doi: 10.1007/b101970.
  * Ben-Ayed and Blair [1990] Omar Ben-Ayed and Charles E. Blair.  Computational Difficulties of Bilevel Linear Programming.  _Operations Research_ , 38(3):556–560, 1990.  ISSN 0030-364X.  doi: 10.1287/opre.38.3.556.
  * Labbé et al. [1998] Martine Labbé, Patrice Marcotte, and Gilles Savard.  A Bilevel Model of Taxation and Its Application to Optimal Highway Pricing.  _Management Science_ , 44(12-part-1):1608–1622, December 1998.  ISSN 0025-1909.  doi: 10.1287/mnsc.44.12.1608.
  * Besançon et al. [2018] Mathieu Besançon, Miguel F. Anjos, Luce Brotcorne, and Juan A. Gomez-Herrera.  A Bilevel Framework for Optimal Price-Setting of Time-and-Level-of-Use Tariffs.  _arXiv:1809.00512 [math]_ , September 2018.
  * Kovács [2019] A. Kovács.  Bilevel programming approach to demand response management with day-ahead tariff.  _Journal of Modern Power Systems and Clean Energy_ , 7(6):1632–1643, November 2019.  ISSN 2196-5420.  doi: 10.1007/s40565-019-0569-7.
  * Yuan et al. [2020] Guanxiu Yuan, Yan Gao, Bei Ye, and Ripeng Huang.  Real-time pricing for smart grid with multi-energy microgrids and uncertain loads: A bilevel programming method.  _International Journal of Electrical Power & Energy Systems_, 123:106206, December 2020.  ISSN 0142-0615.  doi: 10.1016/j.ijepes.2020.106206.
  * Shomalzadeh et al. [2020] Koorosh Shomalzadeh, Jacquelien M. A. Scherpen, and M. Kanat Camlibel.  Energy Market Optimization: From Bilevel to Convex.  _arXiv:2004.08612 [cs, eess, math]_ , November 2020.
  * Grimm et al. [2020] Veronika Grimm, Galina Orlinskaya, Lars Schewe, Martin Schmidt, and Gregor Zöttl.  Optimal Design of Retailer-Prosumer Electricity Tariffs Using Bilevel Optimization.  _Omega_ , page 102327, August 2020.  ISSN 0305-0483.  doi: 10.1016/j.omega.2020.102327.
  * Alves et al. [2020] Maria João Alves, Carlos Henggeler Antunes, and Inês Soares.  Optimizing Prices and Periods in Time-of-use Electricity Tariff Design Using Bilevel Programming.  In Luís Paquete and Christine Zarges, editors, _Evolutionary Computation in Combinatorial Optimization_ , Lecture Notes in Computer Science, pages 1–17, Cham, 2020. Springer International Publishing.  ISBN 978-3-030-43680-3.  doi: 10.1007/978-3-030-43680-3_1.
  * Tang et al. [2019] Rui Tang, Shengwei Wang, and Hangxin Li.  Game theory based interactive demand side management responding to dynamic pricing in price-based demand response of smart grids.  _Applied Energy_ , 250:118–130, September 2019.  ISSN 0306-2619.  doi: 10.1016/j.apenergy.2019.04.177.
  * Henggeler Antunes et al. [2020] Carlos Henggeler Antunes, Maria João Alves, and Billur Ecer.  Bilevel optimization to deal with demand response in power grids: Models, methods and challenges.  _TOP_ , July 2020.  ISSN 1863-8279.  doi: 10.1007/s11750-020-00573-y.
  * Xu et al. [2020] Zhihong Xu, Yan Gao, Muhammad Hussain, and Panhong Cheng.  Demand Side Management for Smart Grid Based on Smart Home Appliances with Renewable Energy Sources and an Energy Storage System.  https://www.hindawi.com/journals/mpe/2020/9545439/, April 2020.  ISSN 1024-123X.
  * Sethi and Sorger [1991] Suresh Sethi and Gerhard Sorger.  A theory of rolling horizon decision making.  _Annals of Operations Research_ , 29(1):387–415, December 1991.  ISSN 0254-5330, 1572-9338.  doi: 10.1007/BF02283607.
  * Samà et al. [2013] Marcella Samà, Andrea D’Ariano, and Dario Pacciarelli.  Rolling Horizon Approach for Aircraft Scheduling in the Terminal Control Area of Busy Airports.  _Procedia - Social and Behavioral Sciences_ , 80:531–552, June 2013.  ISSN 1877-0428.  doi: 10.1016/j.sbspro.2013.05.029.
  * Crama and Pironet [2019] Yves Crama and Thierry Léon Augustin Pironet.  Vehicle allocation problem with uncertain transportation requests over a multi-period rolling horizon.  2019\.  doi: 10.23773/2019_1.
  * Kallabis et al. [2019] Thomas Kallabis, Steven A. Gabriel, and Christoph Weber.  Strategic generation investment using a stochastic rolling-horizon MPEC approach.  July 2019.
  * Heitsch and Römisch [2009] Holger Heitsch and Werner Römisch.  Scenario tree modeling for multistage stochastic programs.  _Mathematical Programming_ , 118(2):371–406, May 2009.  ISSN 1436-4646.  doi: 10.1007/s10107-007-0197-2.
  * Pineda and Morales [2019] Salvador Pineda and Juan Miguel Morales.  Solving Linear Bilevel Problems Using Big-Ms: Not All That Glitters Is Gold.  _IEEE Transactions on Power Systems_ , 34(3):2469–2471, May 2019.  ISSN 1558-0679.  doi: 10.1109/TPWRS.2019.2892607.
  * Kleinert et al. [2019] Thomas Kleinert, Martine Labbé, Fränk Plein, and Martin Schmidt.  There’s No Free Lunch: On the Hardness of Choosing a Correct Big-M in Bilevel Optimization.  April 2019.
  * Pironet [2014] Thierry Pironet.  _Multi-Period Stochastic Optimization Problems in Tranportation Management_.  PhD thesis, HEC Liège, Liège, 2014.
  * Birge and Louveaux [2011] John R. Birge and François Louveaux.  _Introduction to Stochastic Programming_.  Springer Series in Operations Research and Financial Engineering. Springer-Verlag, New York, second edition, 2011.  ISBN 978-1-4614-0236-7.
  * von Niederhäusern [2019] Léonard von Niederhäusern.  _Design and Pricing of New Energy Services in a Competitive Environment_.  PhD thesis, Ecole Centrale, Lille, April 2019.