Coordination of operational planning and real-time

optimization in microgrids
Jonathan Dumas∗ , Selmane Dakir∗ , Clément Liu† , Bertrand Cornélusse∗
∗ Department of electrical engineering and computer science
ULiège, Liège, Belgium
{jdumas, s.dakir, bertrand.cornelusse}@uliege.be
† clement.liu@polytechnique.edu

Abstract—Hierarchical microgrid control levels range from layer by considering power flow and voltage limits. A two-
distributed device level controllers that run at a high frequency stage dispatch strategy for grid connected systems is discussed
to centralized controllers optimizing market integration that run in [4], where the first stage deals with the day ahead schedule,
much less frequently. Centralized controllers are often subdivided
in operational planning controllers that optimize decisions over optimizing capital and operational cost, while the lower level
a time horizon of one or several days, and real-time optimization handles the rescheduling of the units for few hours ahead with
controllers that deal with actions in the current market period. a time resolution of 15 min. A two-stage control strategy for
The coordination of these levels is of paramount importance. a PV BESS-ICE (Internal Combustion Engine) microgrid in
In this paper we propose a value function based approach as implemented in [5]. Discrete Dynamic Programming is used
a way to propagate information from operational planning to
real-time optimization. We apply this method to an environment in the first layer, while the second layer problem is posed as
where operational planning, using day-ahead forecasts, optimizes a Boundary Value Problem. An approach with a high-level
at a market period resolution the decisions to minimize the total deterministic optimizer running at a slow timescale (15 min)
energy cost and revenues, the peak consumption and injection coupled to a low-level stochastic controller running at higher
related costs, and plans for reserve requirements. While real- frequency (1 min) is studied in [6]. A two-layer predictive
time optimization copes with the forecast errors and yields
implementable actions based on real-time measurements. The energy management system for microgrids with hybrid energy
approach is compared to a rule-based controller on three use storage systems consisting of batteries and supercapacitors is
cases, and its sensitivity to forecast error is assessed. considered in [7]. This approach incorporates the degrada-
Index Terms—Hierarchical control, microgrid, optimization tion costs of the hybrid energy storage systems. A practical
Energy Management System for isolated microgrid which
consider the operational constraints of Distributed Energy
I. I NTRODUCTION
Resources, active-reactive power balance, unbalanced system
The hierarchical microgrid control levels divide a global configuration and loading, and voltage dependent loads is
microgrid control problem in time and space [1]. Control levels studied in [8]. A two-layer mixed integer linear programming
range from distributed device level controllers that run at a predictive control strategy was implemented and tested in
high frequency to centralized controllers optimizing market simulation and experimentally in [9], and [10] implemented
integration that run much less frequently. For computation a two-layer predictive management strategy for an off-grid
time reasons, centralized controllers are often subdivided in hybrid microgrid featuring controllable and non-controllable
operational planning controllers that optimize decisions over generation units and a storage system.
a time horizon of one or several days but with a market In this paper we propose a two-layer approach with a value
period resolution (e.g. 15 minutes), and real-time optimization function to propagate information from operational planning
controllers that deal with actions within the current market to real-time optimization. The value function based approach
period. The coordination of these two levels is of paramount shares some similarities with the coordination scheme pro-
importance to achieve the safest and most profitable opera- posed in [11], which is based on stochastic dual dynamic
tional management of microgrids. Microgrid control and man- programming. This paper brings new contributions:
agement can be achieved in several ways. Control techniques
and the principles of energy-storage systems are summarized • The approach is tested by accounting for forecasting
in [1]. A classification of microgrid control strategies into errors and high resolution data monitored on site cor-
primary, secondary, and tertiary levels is done in [2]. The two responding to a "real life" case.
level approach has been intensively studied. A double-layer • The value function approach allows to deal with indeter-
coordinated control approach, consisting of the schedule layer minacy issues. When there are several optimal solutions
and the dispatch layer is adopted in [3]. The schedule layer to the upper level problem, this is accounted for in the
provides an economical operation scheme including state and lower level part, and a bias term can be added to favor
power of controllable units-based on the look-ahead multi-step one type of behavior over another (e.g. charge early).
optimization, while the dispatch layer follows the schedule • This methodology is fully compatible with the energy

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 markets as it is able to deal with imbalance, reserve and related to the current market position, such as for instance the
dynamic selling/purchasing prices. nominated net position of the microgrid over the next market
This paper reports results on an industrial microgrid capable periods. The cost function c gathers all the economical criteria
of on/off grid operation. Generation and consumption fore- considered. The transition function f describes the physical
casts are based on weather forecasts obtained with the MAR and net position evolution of the system. At time instants
model [12]. It is organized as follows. Section II formulates t ∈ {∆τ , 2∆τ , . . .}, with ∆τ the market period, some costs
the problem in an abstract manner. Section III introduces are incurred based on the value of some state variables, which
the novel two-level value function based approach and the are then reset for the next market period. This problem is very
assumptions made. Section IV describes the numerical tests. difficult to solve since the evolution of the system is uncertain,
Section V reports the results. Conclusions are drawn in Section actions have long-term consequences, and are both discrete
VI. Section VIII summarizes the notation. The methodology and continuous. Furthermore, although functions f and c are
used for forecasting is reported as an annex. assumed time invariant, they are generally non-convex and
parameterized with stochastic variables (ωt ).
II. P ROBLEM STATEMENT
A global microgrid control problem can be defined, for III. P ROPOSED METHOD
a given microgrid design and configuration, as operating a
microgrid safely and in an economically efficient manner, by In practice, solving the microgrid optimization problem
harvesting as much renewable energy as possible, operating above amounts, at every time t, to forecasting the stochastic
the grid efficiently, optimizing the service to the demand side, variables ω Tl (t) , then solving the problem1
and optimizing other side goals. We refine this definition below X
and start by making a few assumptions. a?Tl (t) = arg min c(at0 , st0 , ω̂t0 ) (2a)
t0 ∈Tl (t)
A. Assumptions s.t. ∀t ∈ Tl (t), st0 +∆t0 = f (at0 , st0 , ω̂t0 , ∆t0 ),
0
(2b)
In this paper, the control optimizes economical criteria, st0 ∈ St0 , (2c)
which are only related to active power. All devices are
supposed to be connected to the same electrical bus, which can and applying a?t (potentially changing am,?t at some specific
be connected or disconnected from a public grid permanently moments only). As forecasts are valid only for a relatively
or dynamically. Device level controllers offer an interface near future and optimizing over a long time horizon would
to communicate their operating point and constraints (e.g. anyway be incompatible with real-time operation, this prob-
maximum charge power as a function of the current state) lem is approximated by cropping the lookahead horizon to
and implement control decisions to reach power set-points. Ta (t) ⊂ Tl (t). However, market decisions must be refreshed
Fast load-frequency control, islanding management, as well much less frequently than set-points. We thus propose to
as reactive power control, are not in scope. The microgrid is further decompose the problem in an operational planning
a price taker in energy and reserve markets. problem (OPP) for Tam (t) that computes market decisions
B. Formulation X
am,?
T m (t) = arg min cm (am
t0 , st0 , ω̂t0 ) (3a)
In an abstract manner, a microgrid optimization problem a
t0 ∈Tam (t)
can be formulated as follows 0
X s.t. ∀t ∈ Tam (t), st0 +∆τ = f m (am
t0 , st0 , ω̂t0 , ∆τ ) (3b)
min c(at , st , ωt ) (1a) st0 ∈ St0 , (3c)
a
t∈Tl
s.t. ∀t ∈ Tl , st+∆t = f (at , st , ωt , ∆t), (1b) and a real-time problem (RTP) that computes set-points for
st ∈ St . (1c) time t

A controller has to return a set of actions at = (am d

t , at ) at ad,?
t = arg min cd (adt , st , ω̂t ) + vτ (t) (sτ (t) ) (4a)
any time t over the life of the microgrid (Tl ). Actions should s.t. sτ (t) = f d (adt , st , ω̂t , τ (t) − t) (4b)
be taken as frequently as possible to cope with the economical
sτ (t) ∈ Sτ (t) , (4c)
impact of the variability of the demand and generation sides,
but not too often to let transients vanish (e.g. every few with c(at , st , wt ) = cm (am d d
t , st , wt ) + c (at , st , ωt ). The
seconds). The time delta between action at and the next action function vt is the cost-to-go as a function of the state of the
taken is denoted by ∆t , and is not necessarily constant. Some system at the end of the ongoing market period, it regularizes
of these actions are purely market related (am t ), while other decisions of RTP so as to account for the longer term evolution
actions are communicated as set-points to the devices of the of the system. We detail hereunder how we obtain vt . An
microgrid (adt ). The state st = (sm d
t , st ) of the microgrid at overview of the approach is depicted in Figure 1.
time t is thus also made of two parts: sdt represents the state
of the devices, such as the state of charge of a storage system 1 Which is here expressed as a deterministic problem for simplicity, but
or the state of a flexible load, while sm t gathers information should be treated as a stochastic problem in practice.

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 the grid, the costs of purchasing energy from the grid and the
costs for using storage
X  X
OP she she she
Ct0 = ∆τ πd,t0 Cd,t0 ad,t0

t0 ∈Tam (t) d∈D she

Operational Real-time
X
ste ste ste
Forecasting
planning optimization + ∆τ πd,t 0 Pd,t0 ad,t0

d∈D ste
Fig. 1: Hierarchical control procedure illustration. X
nst nst nst
+ ∆τ πd,t 0 Pd,t0 ad,t0

d∈D nst
 
A. Computing the cost-to-go function vτ (t)
X P
+ ∆τ γdsto P d ηdcha acha
d,t0 + disd adis 0
ηd d,t
d∈D sto
The function vt represents the optimal value of (3) as a 
function of the initial state sτ (t) of this problem. If we make − πte0 egri
t0 + π i gri
t t0 .
0 i (8)
the assumption that (3) is modeled as a linear program, the
function vτ (t) is thus convex and piecewise linear. Every DtOP is composed of the peak cost and symmetric reserve
0
evaluation of (3) with the additional constraint2 revenue

sτ (t) = s0 ⊥ µ (5) DtOP

0 = π p δpt0 − πOP
s
rtsym
0 , (9)

yields the value vτ (t) (s0 ) and a supporting inequality (a cut) δpt0 is the peak difference between the previous maximum
historic peak ph and the current peak within the market period
vτ (t) (s) ≥ vτ (t) (s0 ) + µT s. (6) t0 . rtsym
0 is the symmetric reserve provided to the grid within
the current market period t0 .
The algorithm to approximate vτ (t) (s0 ) works as follows:
C. OP constraints
1) estimate the domain of vτ (t) , i.e. the range of states
reachable at time τ (t) and the most probable state that The first set of constraints defines bounds on state and action
will be reached s?τ (t) . variables, ∀t0 ∈ Tam (t)
2) evaluate vτ (t) (s?τ (t) ) and the associated µ?
akd,t0 ≤ 1 ∀d ∈ Dk , ∀k ∈ {ste, she, nst} (10a)
3) repeat step 2 for other state values until all regions of
sto
vτ (t) are explored. acha
d,t0 ≤1 ∀d ∈ D (10b)
sto
adis
d,t0 ≤1 ∀d ∈ D (10c)
Note that if the state is of dimension one and (3) is a linear
sto
program, simplex basis validity information can be used to S d ≤ sd,t0 ≤ S d ∀d ∈ D . (10d)
determine for which part of the domain of vτ (t) the current
cut is tight, else a methodology such as proposed in [13] can The energy flows are constrained, ∀t0 ∈ Tam (t), by
be used. X X
(egri gri
t0 − it0 )/∆τ − (1 − anst )Pd,t
nst
0 + aste ste
d,t0 Pd,t0
d∈D nst d∈D ste
X X
nfl
B. OPP formulation + Cd,t 0 + (1 − ashe she
d,t0 )Cd,t0
d∈D nfl d∈D she
The objective function of the (OPP) implemented for the X
P d acha dis


case study is
+ d,t0 − P d ad,t0 = 0 (11a)
d∈D sto

(et0 − igri
gri
t0 )/∆τ ≤ Etcap (11b)
 
X 0
JTOP
m = CtOP
0 + DtOP
0 (7)
a (t)
t0 ∈Tam (t)
(it0 − egri
gri
t0 )/∆τ ≤ Itcap
0 . (11c)

The dynamics of the state of charge are, ∀d ∈ Dsto

with Operational Planner (OP) the name of this planer. Tam (t)
is composed of 96 values with ∆τ = 15 minutes and Ta = 24
 
P d dis
hours. CtOP models the immediate costs and DtOP the delayed sd,1 − ∆τ P d ηd ad,1 − dis ad,1 = Sdinit
cha cha
(12a)
0 0 η
0
costs at t . CtOP
0 takes into account different revenues and costs  d 
P d dis
related to energy flows: the costs of shed demand, steered and sd,t0 − sd,t0 −∆τ − ∆τ P d ηdcha acha
d,t0 − dis ad,t0 =0
ηd
non steered generation, the revenues from selling energy to
, ∀t0 ∈ Tam (t) (12b)
2 The ⊥ µ notation means that µ is the dual variable of the constraint. sd,τ (t+Ta ) = Sdend . (12c)

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
The set of constraints related to the peak power ∀t0 ∈ Tam (t) E. RTO constraints
The set of constraints that defines the bounds on state and
(igri gri
t0 − et0 )/∆τ ≤ pt0 (13a)
action variables and the energy flows are the same as the OP
−δpt0 ≤ 0 (13b) (10) and (11) by replacing t0 by t, ∆τ by ∆t and considering
−δpt0 ≤ −(pt0 − ph ). (13c) only one period of time t. The next constraint describes the
dynamics of the state of charge ∀d ∈ Dsto and ∀t ∈ Ti (t)
The last constraints define symmetric reserve ∀t0 ∈ Tam (t)  
P d dis
sd,τ (t) − ∆t P d ηdcha acha
d,t − a init
= Sd,t . (17)
s+ (sd,t0 − S d ) ηddis ηddis d,t
rd,t0 ≤ ∀d ∈ Dsto (14a)
∆τ
s+
The set of constraints related to the peak power ∀t ∈ Ti (t)
rd,t dis
0 ≤ P d (1 − ad,t0 ) ∀d ∈ Dsto (14b)
S d − sd,t0

(igri gri
t − et )/∆t ≤ pt,τ (t) (18a)
s−
rd,t0 ≤ ∀d ∈ Dsto (14c) −δpτ (t) ≤ 0 (18b)
ηdcha ∆τ
s− cha
∀d ∈ Dsto −δpτ (t) ≤ −(pτ (t−∆τ ),τ (t) − ph ) (18c)
rd,t0 ≤ P d (1 − ad,t0 ) (14d)
X
s+
X pτ (t−∆τ ),τ (t) = βpτ (t−∆τ ),t + (1 − β)pt,τ (t) (18d)
rsym, OP ≤ rd,t0 +
ste
Pd,t ste
0 (1 − ad,t0 )

d∈D sto d∈D ste with β = 1 − ∆t /∆τ . The last set of constraints defining the
X
+ nst
Pd,t 0 (1 −a )nst
(14e) symmetric reserve are the same as the OP (14) by replacing
d∈D nst t0 by t, rsym, OP by rsym, RTO and adding ∀t ∈ Ti (t)
X
she she
+ Cd,t0 (1 − ad,t0 ) − ∆rsym ≤ 0 (19a)
d∈D she sym sym, OP sym, RTO
X X − ∆r ≤ −(r −r ). (19b)
s−
rsym, OP ≤ rd,t0 +
ste ste
Pd,t0 ad,t0

d∈D sto d∈D ste IV. T EST DESCRIPTION

X X
nst nst she she Our case study is based on the MiRIS microgrid located
+ Pd,t0a + Cd,t0 ad,t0 . (14f)
d∈D nst d∈D she at the John Cockerill Group’s international headquarters in
Seraing, Belgium3 . It is composed of PV, several energy
D. RTP formulation storage devices and a non-sheddable load. The load and PV
The objective function of the (RTP) implemented for the data we use come from on-site monitoring. All data, including
case study is the weather forecasts, are available on the Kaggle platform4 .
The case study consists in comparing RTO-OP to a Rule Based
JtRT O = CtRT O + DtRT O + vτ (t) (sτ (t) ) (15) Controller (RBC) for three configurations of the installed PV
capacity, cf. Table I. The RBC prioritizes the use of PV
with Real Time Optimizer (RTO) the name of this controller. production for the supply of the electrical demand. If the
CtRT O models the immediate costs, DtRT O the delayed costs microgrid is facing a long position, it charges the battery.
and vτ (t) (sτ (t) ) the cost-to-go function of the state of the And if this one is fully charged it exports to the main grid.
system at time t within a current market period. CtRT O is the If the microgrid is facing a short position it prioritizes the
same as CtOP 0 by replacing t0 by t, ∆τ by ∆t and considering use of battery to supply the demand. And if this one is fully
only one period of time t. DtRT O is composed of the peak discharged it imports from the main grid. This controller does
cost and symmetric reserve penalty costs not take into account any future information (PV, consumption
forecasts, energy prices, etc) or market information such as
DtRT O = π p δpτ (t−∆τ ),τ (t) + sTt SO πRT
s
O ∆r
sym
, (16) the peak of symmetric reserve. Case 3 is the result of a sizing
study that defined the optimal device sizes given the PV and
δpτ (t−∆τ ),τ (t) is the peak difference between the previous consumption data. The sizing methodology used is described
maximum historic peak ph and the current peak within the in [14].
market period computed by RTO. The difference with OP
Figure 2 shows the PV & consumption data over the
relies on its computation as at t the market period is not
simulation period: from M ay 20, 2019 to June 16, 2019. The
finished. Thus the peak within this market period is computed
selling price π e is constant, the purchasing price is composed
by adding the peak from the beginning of the market period
of a day πdi and night prices πni . Day prices are applicable
to t and the one resulting from the actions taken from t to the
from 5 a.m. to 8 p.m. (UTC) during week days and night prices
end of the market period. ∆rsym is the difference between the
are applicable from 8 p.m. to 5 a.m. during week days and
symmetric reserve computed by OP and the current reserve
during the entire weekend. The peak mechanism is taken into
within the market period computed by RTO. sTt SO is the
account with a constant peak price π p and an initial maximum
reserve activation signal to activate the tertiary symmetric
reserve. It is set by the TSO (0 if activated, else 1). The 3 https://johncockerill.com/fr/energy/stockage-denergie/

activation occurs at the beginning of the next market period. 4 https://www.kaggle.com/jonathandumas/liege-microgrid-open-data

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 TABLE I: Case studies parameters and data statistics. 1200
C
1000 PV 3
800 PV 2
Case PVp PV PVmax PVmin PVstd PV 1
600

kW
1 400 61 256 0 72 400
2 875 133 561 0 157 200
0
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
3 1750 267 1122 0 314 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
1200
Case Cp C Cmax Cmin Cstd C
1000 PV 3
800 PV 2
1-3 1000 153 390 68 72 PV 1
600

kW
cha dis init 400
Case Sp S, S P, P η ,η S
200
1-3 1350 1350, 0 1350, 1350 0.95, 0.95 100 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Case ph , π p
I cap
E cap
πdi , i
πn πe 00:0 01:0 02:0 03:0 04:0 05:0 06:0 07:0 08:0 09:0 10:0 11:0 12:0 13:0 14:0 15:0 16:0 17:0 18:0 19:0 20:0 21:0 22:0 23:0

1-3 150, 40 1500 1500 0.2, 0.12 0.035 Fig. 2: Top: PV & consumption simulation data. Bottom:
zoom on June 12, 2019.

TABLE II: Results without symmetric reserve.

historic peak ph . Storage systems are initially fully charged.
The PV and consumption data have a 1 second resolution, Case 1 cE cp ct ∆p Itot Etot

meaning the RTO could compute its optimization problem RBC 10.13 6.68 16.81 167 61 0

each five to ten seconds in operational mode5 . However, to RTO-OPRN N 10.37 3.62 13.99 91 64 1
GBR
maintain a reasonable simulation time RTO is called every RTO-OP 10.25 5.27 15.53 132 63 1
minute6 . The OP computes a planning on a quarterly basis RTO-OP? 10.24 0.99 11.23 25 64 1
corresponding to the belgian market period. The computation Case 2 cE cp ct ∆p Itot Etot
time of the RTO on a regular computer is around a few seconds RBC 3.19 4.85 8.04 121 22 7
and the OP around twenty seconds. In total the simulation RTO-OPRN N 4.78 2.87 7.65 72 31 15
computation time is up to a few hours. The OP computes RTO-OPGBR 4.30 4.90 9.2 123 28 13
quarterly a planning based on PV and consumption twenty- RTO-OP ?
4.06 0 4.06 0 26 10
four ahead forecasts. The weather based forecast methodology
Case 3 cE cp ct ∆p Itot Etot
is described in details in Annex IX. Two "classic" deterministic
RBC -2.13 4.12 1.99 105 3 77
techniques are implemented, a Recurrent Neural Network
RTO-OPRN N -1.66 4.12 2.46 105 7 80
(RNN) with the keras python library [15] and a Gradient
RTO-OPGBR -1.67 4.23 2.56 106 7 81
Boosting Regression (GBR) with the scikit-learn python li-
RTO-OP? -1.90 0 0 0 5 79
brary [16]. These models use as input the weather forecasts
provided by the Laboratory of Climatology of the university
of Liège, based on the MAR regional climate model [12]. It is
an atmosphere model designed for meteorological and climatic capacity the higher the peak and energy costs. The RTO-OP?
research, used for a wide range of applications, from km-scale provides the minimal peak cost whereas the RBC provides the
process studies to continental-scale multi-decade simulations. minimal energy cost on all cases. However, RTO-OP? achieves
To estimate the impact of the PV and consumption forecast the minimal total cost, composed of the energy and peak
errors on the controllers, the simulation is performed with costs. This simulation illustrates the impact of the forecasts on
the OP having access to the PV and consumption future the RTO-OP behavior. The RNN forecaster provides the best
values (RTO-OP? ). Then, the simulation is performed with results but the RTO-OPRN N is still a long way to manage
the symmetric reserve mechanisms to cope with the forecasts the peak as RTO-OP? due to the forecasting errors. The peak
s
errors. A constant symmetric reserve price πOP for the OP cost strongly penalizes the benefits as it applies on all the year
s
and a penalty reserve πRT O for the RTO are set to 20 (e/ ahead once it has been reached.
kW). In case 3, all the controllers except RTO-OP? reached
the maximum peak on June 12, 2019 around 10:30 a.m. as
V. N UMERICAL RESULTS
shown on Figure 4. Figure 2 shows a sudden drop in the PV
A. No symmetric reserve production around 10 a.m. that is not accurately forecasted
Table II provides the simulation results without taking into
account the symmetric reserve. The smaller the PV installed 1000
900 PV
800 RNN
5 CPLEX 12.9 is used to solve all the optimization problems, on an Intel 700 GBR
600
500
kW

core i7-8700 3.20 GHz based computer with 12 threads and 32 GB of RAM. 400
The average computation time per optimization problem composed of the OP 300
200
and RTO is a of few seconds. 100
6 The dataset is composed of 28 days with an average computation time 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
05:0 06:0 07:0 08:0 09:0 10:0 11:0 12:0 13:0 14:0 15:0 16:0 17:0 18:0 19:0 20:0 21:0 22:0 23:0 00:0 01:0 02:0 03:0 04:0 05:0 06:0 07:0
of two hours to solve 1440 optimization problems per day, with one minute
resolution, leading to a total of two days for the entire dataset. Fig. 3: Case 3 PV forecast on June 12, 2019, 06h00 UTC.

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 7
RTO OP
TABLE III: Results with symmetric reserve.
6 RTO OP RNN
5 RTO OP GBR
RBC Case 1 cE cp ct ∆p Itot Etot
4
3 RTO-OPRN N
k

10.50 2.12 12.62 53 65 3

2
1 RTO-OP? 10.47 2.75 13.22 69 65 2
0
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- Case 2 cE cp ct ∆p Itot Etot
7
RTO OP RTO-OPRN N 5.33 0.04 5.37 1 41 27
6 RTO OP RNN
5 RTO OP GBR RTO-OP? 4.78 0.99 5.77 25 35 20
4 RBC
3
k

Case 3 cE cp ct ∆p Itot Etot

2
1 RTO-OPRN N -0.04 0 -0.04 0 24 99
00 1 2 3 4 5 6 7 8 9 0 1 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6
2 2 2 2 2 2 2 2 2 2 3 3 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
05 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
- RTO-OP? -0.15 0 -0.15 0 23.2 98
7
RTO OP
6 RTO OP RNN
5 RTO OP GBR
RBC 100
4 90
3 80
k

2 70
1 60

SOC %
0
50
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 40
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 30
20 RTO OP RNN s = 0
Fig. 4: Case 1 (top), 2 (middle), 3 (bottom) cumulative peak 10 RTO OP RNN s = 20
00 1 2 3 4 5 6 7 8 9 0 1 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6
costs. 2 2 2 2 2 2 2 2 2 2 3 3 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
RN N
Fig. 6: Case 3 SOC comparison for RTO-OP with and
without symmetric reserve.
100
90
80
70
60
SOC %

50
40 RTO OP
B. Results with symmetric reserve
30 RTO OP RNN
20
10
RTO OP GBR
RBC
Table III provides the simulation results by taking into
00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22:0 23:0
: 0 : 0 :0 : 0 :0 : 0 : 0 : 0 : 0 :0 : 0 : 0 : 0 : 0 : 0 : 0 : 0 : 0 : 0 : 0 : 0 : 0 account the symmetric reserve. Figure 6 depicts on case 3
1400 RTO OP
the behavior differences between RTO-OPRN N without and
1200 RTO OP RNN with symmetric reserve. Figures 7 and 8 show the SOC and
1000 RTO OP GBR
800 RBC peaks costs evolution of case 2 & 1. The controller tends to
600 maintain a storage level that allows RTO-OPRN N to better
kW

400
200 cope with forecast error. Indeed for case 3 there is no more
0 peak reached by RTO-OPRN N , only 1 kW for case 2 and
200
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 it has been almost divided by two for case 1. However, this
00:0 01:0 02:0 03:0 04:0 05:0 06:0 07:0 08:0 09:0 10:0 11:0 12:0 13:0 14:0 15:0 16:0 17:0 18:0 19:0 20:0 21:0 22:0 23:0
behavior tends to increase the energy cost if the PV production
Fig. 5: Case 3 SOC (top) and net export power (bottom) on
is important in comparison with the consumption, such as
June 12, 2019.
for case 3. Indeed, the controller will tend to store more
energy into the battery instead of exporting it. RTO-OP? did
not perform better with the symmetric reserve. In fact the
symmetric reserve is in competition with the peak management
by the RNN and GBR forecasters as shown in Figure 3. This and the RTO-OP? tends to not discharge completely the battery
prediction leads to a non accurate planning of OP. Thus, the even if it is required to avoid a peak. In case 2, the peak is
RTO cannot anticipate this drop and has to import at the reached on June 12, 2019 around 08:00. The controller could
last minute energy to balance the microgrid. Figure 5 shows have avoided it by totally discharging the battery but did not to
the controllers behavior on June 12, 2019 where the peak is maintain the reserve level. This is the same behavior on case
reached. In case 2 all controllers reached the same peak as 1 where the peak could have been limited if all the battery
in case 3 except RTO-OPRN N that reached a smaller one on was discharged. There is an economic trade-off to reach to
June 5, 2019. The forecasts accuracy explains this behavior as manage the peak and the reserve simultaneously depending
in case 3. Finally in case 1, each controller reached a different on the valorization or not on the market of the symmetric
peak. The smallest one is achieved by the RTO-OP? , followed reserve. The reserve can also be valorized internally to cope
by the RTO-OPRN N . These cases show that the RTO-OP with non or difficult forecastable events such as a sudden drop
controller optimizes PV-storage usage, and thus requires less of the export or import limits due to loss of equipment or grid
installed PV capacity for a given demand level. This result congestion.
was expected as the peak management is not achieved by the
RBC and becomes critical when the PV production is smaller VI. C ONCLUSION
than the consumption. This simulation also demonstrates the A two-level value function based approach was introduced
forecast accuracy impact on the RTO-OP behavior. as a solution method for a multi-resolution microgrid optimiza-

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 3
RTO OP RNN R EFERENCES
RTO OP
2 [1] O. Palizban, K. Kauhaniemi, and J. M. Guerrero, “Microgrids in active
network management part i: Hierarchical control, energy storage, virtual
k

1
power plants, and market participation,” Renewable and Sustainable
Energy Reviews, vol. 36, pp. 428–439, 2014.
[2] D. E. Olivares, A. Mehrizi-Sani, A. H. Etemadi, C. A. Cañizares,
0
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 R. Iravani, M. Kazerani, A. H. Hajimiragha, O. Gomis-Bellmunt,
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
M. Saeedifard, R. Palma-Behnke et al., “Trends in microgrid control,”
3
RTO OP RNN IEEE Transactions on smart grid, vol. 5, no. 4, pp. 1905–1919, 2014.
RTO OP [3] Q. Jiang, M. Xue, and G. Geng, “Energy management of microgrid in
2 grid-connected and stand-alone modes,” IEEE transactions on power
systems, vol. 28, no. 3, pp. 3380–3389, 2013.
k

1 [4] X. Wu, X. Wang, and C. Qu, “A hierarchical framework for genera-

tion scheduling of microgrids,” IEEE Transactions on Power Delivery,
0
vol. 29, no. 6, pp. 2448–2457, 2014.
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- [5] J. Sachs and O. Sawodny, “A two-stage model predictive control strategy
for economic diesel-pv-battery island microgrid operation in rural areas,”
Fig. 7: Case 1 (top) and 2 (bottom) cumulative peak costs. IEEE Transactions on Sustainable Energy, vol. 7, no. 3, pp. 903–913,
2016.
100 [6] S. R. Cominesi, M. Farina, L. Giulioni, B. Picasso, and R. Scattolini,
90 RTO OP RNN “A two-layer stochastic model predictive control scheme for microgrids,”
80 RTO OP
70
IEEE Transactions on Control Systems Technology, vol. 26, no. 1, pp.
60 1–13, 2017.
SOC %

50 [7] C. Ju, P. Wang, L. Goel, and Y. Xu, “A two-layer energy management

40
30 system for microgrids with hybrid energy storage considering degrada-
20 tion costs,” IEEE Transactions on Smart Grid, vol. 9, no. 6, pp. 6047–
10
0 6057, 2017.
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- [8] B. V. Solanki, C. A. Cañizares, and K. Bhattacharya, “Practical energy
100 management systems for isolated microgrids,” IEEE Transactions on
90 RTO OP RNN
80 RTO OP Smart Grid, 2018.
70 [9] S. Polimeni, L. Moretti, G. Manzolini, S. Leva, L. Meraldi, and
60 P. Raboni, “Numerical and experimental testing of predictive EMS
SOC %

50
40 algorithms for PV-BESS residential microgrid,” IEEE-PES Powertech
30 (accepted), p. 6, 2019.
20
10 [10] L. Moretti, S. Polimeni, L. Meraldi, P. Raboni, S. Leva, and G. Man-
0 zolini, “Assessing the impact of a two-layer predictive dispatch algorithm
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
0 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
5 -
on design and operation of off-grid hybrid microgrids,” Renewable
Fig. 8: Case 1 (top) and 2 (bottom) SOC. Energy, vol. 143, pp. 1439–1453, 2019.
[11] R. Kumar, M. J. Wenzel, M. J. Ellis, M. N. ElBsat, K. H. Drees, and
V. M. Zavala, “A stochastic dual dynamic programming framework for
multiscale mpc,” IFAC-PapersOnLine, vol. 51, no. 20, pp. 493–498,
tion problem. The value function computed by the operational 2018.
[12] X. Fettweis, J. Box, C. Agosta, C. Amory, C. Kittel, C. Lang, D. van
planner based on PV and consumption forecasts allows to cope As, H. Machguth, and H. Gallée, “Reconstructions of the 1900–2015
with the forecasting uncertainties. The real time controller greenland ice sheet surface mass balance using the regional climate
solves an entire optimization problem including the future in- MAR model,” Cryosphere (The), vol. 11, pp. 1015–1033, 2017.
[13] A. Bemporad, A. Garulli, S. Paoletti, and A. Vicino, “A greedy approach
formation propagated by the value function. This approach has to identification of piecewise affine models,” in International Workshop
been tested on the MiRIS microgrid case study with PV and on Hybrid Systems: Computation and Control. Springer, 2003, pp.
consumption data monitored on site. The results demonstrate 97–112.
[14] S. Dakir and B. Cornélusse, “Combining optimization and simulation
the efficiency of this method to manage the peak in comparison for microgrid sizing,” Submitted to the 21st Power Systems Computation
with a Rule Based Controller. This test case is completely Conference (PSCC 2020), 2019.
reproducible as all the data used are open (PV, consumption [15] F. Chollet et al., “Keras,” https://keras.io, 2015.
[16] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion,
monitored and forecasted including the weather forecasts). O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vander-
The proposed approach can be extended in several ways. The plas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duch-
deterministic formulation of the operational planning problem esnay, “Scikit-learn: Machine learning in Python,” Journal of Machine
Learning Research, vol. 12, pp. 2825–2830, 2011.
could be extended to a stochastic formulation, in order to cope [17] S. B. Taieb, G. Bontempi, A. F. Atiya, and A. Sorjamaa, “A review and
with probabilistic forecasts. Balancing market mechanisms comparison of strategies for multi-step ahead time series forecasting
could be introduced. Finally, the approach could be extended based on the nn5 forecasting competition,” Expert systems with appli-
cations, vol. 39, no. 8, pp. 7067–7083, 2012.
to a community by considering several entities inside the
microgrid.

VII. ACKNOWLEDGMENT
The authors would like to thank John Cockerill and Nethys
for their financial support, and Xavier Fettweis of the Labo-
ratory of Climatology of ULiège who produced the weather
forecasts based on the MAR regional climate model.

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 VIII. N OTATION • st microgrid state at time t
Set and indices • smt information related to the current market position
• sdt state of the devices
• d index of a device • vt the cost-to-go function
• t, t0 indexes of a RTO and OP time periods • ω̂ forecast of a random vector ω
• τ (t) beginning of the next market period at time t • X average of a variable X (kW)
• Ti (t) = {t, t + ∆t, ..., t + Ti } set of RTO time periods • Xmax , Xmin maximum and minimum of X (kW)
• Tam (t) = {τ (t), τ (t) + ∆τ, ..., τ (t + Ta )} set of OP time • Xstd standard deviation of X (kW)
periods • cE , cp , ct energy, peak and total costs (ke)
• Ta , Tl time durations, with Ta  Tl • ∆p peak increment (kW)
• Dk set of non-flexible loads (k = nfl), sheddable loads • Itot , Etot total import and export (MWh)
(k = she), steerable generators (k = ste), non-steerable
generators (k = nst), storage devices (k = sto) IX. A NNEX : FORECASTING METHODOLOGY
The inputs of the forecasting method are historical and
Parameters
external data, a forecasting horizon HT , a resolution, and a
• ∆t time delta between t and the market period (minutes) forecast frequency. The outputs are the PV production and
• ∆τ market period (minutes) the consumption. In this study, the input data are weather
• HT forecasting horizon (hours) forecasts and past PV production and consumption series.
• ω̂ forecast of a random vector ω The horizon is the time range of the forecasts from a few
• η cha , η dis charge and discharge efficiencies (%) hours to several hours or days. The resolution is the time
• P , P maximum charging and discharging powers (kW) discretization of the forecast from a few minutes to several
nfl
• Cd,t non-flexible power consumption (kW) hours. The forecast frequency indicates the periodicity at
she
• Cd,t flexible power consumption (kW) which the forecasts are computed. For instance, a forecasting
init
• Sd,t initial state of charge of battery d (kWh) module with HT = 24 hours, a resolution and periodicity of
• ph maximum peak over the last twelve months (kW) 15 minutes, computes each quarter, a quarterly forecast for
• π p yearly peak power cost (e/kW) the twenty-four hours ahead. This paper focuses on the real
s
• πOP unitary revenue for providing reserve (e/kW) time control of microgrids based on a planning that requires
s
• πRT O unitary RTO symmetric reserve penalty (e/kW) a forecast horizon of a few hours up to a few days.
k
• πd,t cost of load shedding (k = she), generating energy Two "classic" deterministic techniques are implemented,
(k = ste), curtailing generation (k = nst) (e/kWh) a Recurrent Neural Network (RNN) with the keras python
sto
• γd,t fee to use the battery d (e/kWh) library [15] and a Gradient Boosting Regression (GBR) with
• πte , πti energy prices of export and import (e/kWh) the scikit-learn python library [16]. The RNN is a Long Short
• πdi , πni energy prices of day and night imports (e/kWh) Term Memory (LSTM) with one hidden layer composed of
• I cap , E cap maximum import and export limits (kW) 2 × n + 1 neurons with n the number of input features.
• PVp , Cp PV and consumption capacities (kW) Both techniques are implemented with a Multi-Input Multi-
• Sp storage capacity (kWh) Output (MIMO) approach [17]. The MIMO strategy consists
• S, S maximum and minimum battery capacities (kWh) in learning only one model, m, b as follows
 
Forecasted or computed variables τ1 τH T τ0 τ−4 τ1 τH T
[wb , ..., w
b ]=m
b w , ..., w , w bi , ..., w
bi . (20)
• at action at t
• amt purely market related actions With w b the variable to forecast (PV, consumption, etc), w bi
• adt set-points to the devices of the microgrid the forecast of the ith weather variable such as direct solar
• akd,t fraction of load shed (k = she), generation activated irradiance, wind speed, air ambient temperature, etc. The
(k = ste), generation curtailed (k = nst) ([0, 1]) forecast is computed each quarter and composed of HT /∆τ
• acha dis
d,t , ad,t fraction of the maximum charging and discharg- y τ1 , ..., ybτHT ]. In our case study, HT /∆τ = 96 with
values [b
ing powers used for battery d ([0, 1]) HT = 24 h and ∆τ = 15 minutes. The forecasting process is
• egri gri
t , it energy export and import (kWh) implemented as a rolling forecast methodology. The Learning
• δpt0 OP peak difference between peak at t0 and ph (kW) Set (LS) is refreshed every six hours. The LS is limited to
• δpτ (t−∆τ ),τ (t) RTO peak difference between peak at τ (t) the week preceding the forecasts, to maintain a reasonable
and ph (kW) computation time.
• sTt SO TSO symmetric reserve signal (0; 1) The forecasts are evaluated using three deterministic met-
• rsym symmetric reserve (kW) rics: the Normalized Mean Absolute Error (NMAE), the
• ∆rsym reserve difference between OP and RTO (kW) Normalized Root Mean Squared Error (NRMSE) and the
s+ s−
• rd,t 0 , rd,t0 upward and downward reserves of power Normalized Energy Measurement Error (NEME). The NEME
available and provided by storage device d (kW) is an NMAE of the energy summed over the entire forecasting
• sd,t state of charge of battery d (kWh) horizon. The mean scores N M AEHT , N RM SEHT and

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020
 200
180 NMAE
NRMSE
NMAE mean
NRMSE mean
N EM EHT for a forecasting horizon HT are computed over
160
140 NEME NEME mean the entire simulation data set. The normalizing coefficient
120
100
for computing the NMAE and the NRMSE is the mean of
%

80
60
40 the absolute value of the PV and consumption over all the
20
0 simulation data set. Figures 9 and 10 provide the scores for
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
200 both GBR and RNN techniques computed for each quarter of
180 NMAE NMAE mean
160 NRMSE NRMSE mean
140 NEME NEME mean the simulation data set.
120
100
%

80
60
40
20
0
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-

Fig. 9: PV forecast scores for GBR (top) and RNN (bottom).

80
NMAE NMAE mean
NRMSE NRMSE mean
60 NEME NEME mean
40
%

20

0
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-
80
NMAE NMAE mean
NRMSE NRMSE mean
60 NEME NEME mean
40
%

20

0
20 21 22 23 24 25 26 27 28 29 30 31 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16
05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 05- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06- 06-

Fig. 10: Consumption forecast scores for GBR (top) and

RNN (bottom).

21st Power Systems Computation Conference Porto, Portugal — June 29 – July 3, 2020
PSCC 2020