# Graph-based Global Robot Simultaneous Localization and Mapping using
Architectural Plans

Muhammad Shaheer1, Jose Andres Millan-Romera1, Hriday Bavle1, Jose Luis
Sanchez-Lopez1,
Javier Civera2 and Holger Voos1 1Authors are with the Automation and Robotics
Research Group, Interdisciplinary Centre for Security, Reliability and Trust
(SnT), University of Luxembourg (UL). Holger Voos is also associated with the
Faculty of Science, Technology and Medicine, University of Luxembourg,
Luxembourg. {muhammad.shaheer, jose.millan, hriday.bavle,
joseluis.sanchezlopez<EMAIL_ADDRESS>is with I3A, Universidad de
Zaragoza, Spain<EMAIL_ADDRESS>work was partially funded by the Fonds
National de la Recherche of Luxembourg (FNR) under the project
17097684/RoboSAUR, by a partnership between the SnT-UL and Stugalux
Construction S.A., and by the Spanish and Aragón governments (projects
PID2021-127685NB-I00, TED2021-131150B-I00 and DGA_FSE-T45_20R). For the
purpose of Open Access, the author has applied a CC BY 4.0 public copyright
license to any Author Accepted Manuscript version arising from this
submission.

###### Abstract

In this paper, we propose a solution for graph-based global robot simultaneous
localization and mapping (SLAM) using architectural plans. Before the start of
the robot operation, the previously available architectural plan of the
building is converted into our proposed architectural graph (A-Graph). When
the robot starts its operation, it uses its onboard LIDAR and odometry to
carry out an online SLAM relying on our situational graph (S-Graph), which
includes both, a representation of the environment with multiple levels of
abstractions, such as walls or rooms and their relationships, as well as the
robot poses with their associated keyframes. Our novel graph-to-graph matching
method is used to relate the aforementioned S-Graph and A-Graph, which are
aligned and merged, resulting in our novel informed Situational Graph (iS-
Graph). Our iS-Graph not only provides graph-based global robot localization,
but it extends the graph-based SLAM capabilities of the S-Graph by
incorporating into it the prior knowledge of the environment existing in the
architectural plan.

##  I Introduction

The construction industry is increasingly using mobile robots, which offer
numerous potential benefits. These robots can significantly reduce costs by
regularly inspecting ongoing construction sites and monitoring progress.
However, most robots used in construction are either teleoperated or operate
semi-autonomously, primarily due to the perception challenges associated with
the ever-changing nature of construction sites. To enable fully autonomous
operation, it would be advantageous for these robots to possess comprehensive
prior knowledge of the construction site’s geometry. By combining this prior
knowledge with sensor readings during real-time operation, these robots could
achieve robust and accurate global localization within construction sites.

Digital architectural plans, such as Building Information Modelling (BIM) [1],
provide a means of capturing and communicating information about a
construction site and incorporating it as prior knowledge about the scene.
Works such as [2, 3, 4] have addressed the problem of extracting relevant
structural knowledge from BIM and using it for real-time robot localization.
However, these methods only extract geometric information from the BIM and do
not leverage the topological and relational information also available in it,
which limits the robustness and accuracy in complex and changing construction
sites.

To tackle this problem, we present a novel approach to localize robots
leveraging not only geometry but also higher-level hierarchical information
from architectural plans. We present in this paper how to model the BIM
information in the form of a graph that we denote as Architectural Graph
(A-Graph), and then match and merge with the online Situational Graph
(S-Graph) [5], [6] that the robot builds as it navigates the environment. As a
key aspect, translating low-level geometry into high-level features in both
graphs is what allows a robust matching between such different inputs.

Our method can be divided into three main stages. In the first one, an A-Graph
is created for a given environment with nodes representing the semantic
features available in a BIM model, specifically, wall-surfaces, doors, and
rooms as the graph nodes, and edges containing the relevant relational
information such as two wall surfaces comprising a wall, four wall surfaces
connecting room and rooms connected through doors. In the second stage,
running in real-time onboard the robot, a S-Graph is estimated using 3D LiDAR
measurements. The nodes of our S-Graph correspond to semantics such as wall
surfaces and rooms, and the edges correspond to constraints between these wall
surfaces and the relevant room nodes. Finally, to localize the robot within
this environment, a graph-matching algorithm is proposed utilizing
hierarchical information from both graphs to provide the best match candidates
finally resulting in informed (iS-Graphs) that fuses the information of both.
This last graph will be the one used for global localization.

Figure 1: Overview of our approach. We generate offline an Architectural Graph
(A-Graph) from a BIM model. A robot estimates online a Situational Graph
(S-Graph) from its sensors. We do graph matching between the two, align them,
and merge their information. This generates the final iS-Graphs, which is
utilized by the robot to be localized with respect to the BIM.

##  II System Overview

Fig. 1 shows an overview of the proposed approach. Firstly the BIM information
is extracted for a given environment to create the two-layered Architectural
Graphs (Section. III) in an offline manner, with all the elements extracted in
the BIM frame of reference $B$. As the robot navigates within the given
environment, an online Situational Graph (Section. II-A) is estimated by the
robot in the map frame of reference $M$. In parallel, we run our the Graph
Matching method (Section. IV) to provide match candidates between the existing
A-Graphs and the current S-Graph. Finally, after retrieving the best match
candidate, our Graph Merging (Section. V) provides the merged iS-Graph
utilized for the global localization with respect to the BIM frame of
reference $B$.

###  II-A Situational Graphs (S-Graphs)

S-Graphs are four-layered optimizable hierarchical graphs built online using
3D LiDAR measurements. The full details of the S-Graphs we use in this work
can be found in [5],[6]. In brief, their four layers can be summarized as:

Keyframes Layer. It consists of the robot poses factored as
$\leftidx{{}^{M}}{\mathbf{x}}_{R_{i}}\in SE(3)$ nodes in the map frame $M$
with pairwise odometry measurements constraining them.

Walls Layer. It consists of the planar wall-surfaces
$\leftidx{{}^{M}}{\boldsymbol{\pi}}_{i}\in\mathbb{R}^{3}$ extracted from the
3D LiDAR measurements and factored using minimal plane parameterization. The
planes observed by their respective keyframes are factored using pose-plane
constraints.

Rooms Layer: It consists of two-wall rooms
$\leftidx{{}^{M}}{\boldsymbol{\gamma}}_{i}\in\mathbb{R}^{2}$ or four-wall
rooms $\leftidx{{}^{M}}{\boldsymbol{\rho}}_{i}\in\mathbb{R}^{2}$, each
constraining either two or four detected wall-surfaces respectively.

Floors Layer: It consists of a floor node
$\leftidx{{}^{M}}{\boldsymbol{\xi}}_{i}\in\mathbb{R}^{2}$ positioned at the
center of the current floor level and constraining all the rooms present at
that floor level.

##  III Architectural Graphs (A-Graphs)

We extract relevant information from BIM models into two-layered optimizable
graphs denoted as A-Graphs. In the lowest-level layer, we will model the
geometry of the walls, and in the highest level the rooms. Room-to-wall
constraints connect the two layers and neighboring rooms are constrained by
doorways. The specific formulation is detailed below.

A-Walls Layer. This layer extracts all the information about the walls and
wall surfaces from the BIM and connects them with appropriate wall-to-wall-
surface edges.

Wall-Surfaces: Wall-surfaces are planar entities
$\leftidx{{}^{B}}{\boldsymbol{\pi}}$ extracted in the BIM frame of reference
$B$. All the wall-surfaces are converted to their Closest Point (CP)
representation, as in [6]. Wall-surface normals with their component
$\leftidx{{}^{B}}{{n}_{x}}$ greater than $\leftidx{{}^{B}}{{n}_{y}}$ are
classified as $x$-wall-surfaces, and wall surfaces whose normal component
$\leftidx{{}^{B}}{{n}_{y}}$ is greater than the normal component
$\leftidx{{}^{B}}{{n}_{x}}$ are defined as $y$-wall-surfaces. These are
initialized in the graph as
$\leftidx{{}^{B}}{\boldsymbol{\pi}}=[\leftidx{{}^{B}}\phi,\leftidx{{}^{B}}\theta,\leftidx{{}^{B}}d]$,
where $\leftidx{{}^{B}}\phi$ and $\leftidx{{}^{B}}\theta$ stand for the
azimuth and elevation of the plane in frame $B$ and $\leftidx{{}^{B}}d$ is the
perpendicular distance in $B$.

Walls: We introduce a novel semantic entity with respect to [6] in the form of
a Wall $\boldsymbol{\omega}\in\mathbb{R}^{3}$, consisting of two planar wall-
surfaces. Two opposed planar wall-surface entities either in $x$-direction or
$y$-direction with similar perpendicular distance $\leftidx{{}^{B}}d$ can be
classified as a part of a single wall entity. The wall center
$\leftidx{{}^{B}}{\boldsymbol{\omega}_{x_{i}}}$ for two opposed $x$-direction
wall-surfaces is computed as:

$\leftidx{{}^{B}}{\mathbf{w}_{x_{i}}}=\frac{1}{2}\big{[}\lvert{\leftidx{{}^{B}}{d_{x_{1}}}\rvert}\cdot\leftidx{{}^{B}}{\mathbf{n}_{x_{1}}}-{\lvert\leftidx{{}^{B}}{d_{x_{1}}}\rvert}\cdot\leftidx{{}^{B}}{\mathbf{n}_{x_{2}}}\big{]}+\lvert\leftidx{{}^{B}}{{d_{x_{2}}}\rvert}\cdot\leftidx{{}^{B}}{\mathbf{n}_{x_{2}}}$

$\displaystyle\leftidx{{}^{B}}{\boldsymbol{\omega}_{x_{i}}}=\leftidx{{}^{B}}{\mathbf{w}_{x_{i}}}+\big{[}\leftidx{{}^{B}}{\mathbf{\mathfrak{s}}_{i}}-[\
\leftidx{{}^{B}}{\mathbf{\mathfrak{s}}_{i}}\cdot\leftidx{{}^{B}}{\hat{\mathbf{w}}_{x_{i}}}]\
\cdot\hat{\leftidx{{}^{B}}{\mathbf{w}}_{x_{i}}}\big{]}$ (1)

where $\leftidx{{}^{B}}{\mathbf{\mathfrak{s}}_{i}}\in\mathbb{R}^{3}$ is the
starting point for a given BIM wall and $\mathbf{n}$ and $d$ are the plane
normals and distance. For Eq. III to hold true, all plane normals are
converted to point away from the BIM frame of reference as in [6]. The wall
center along with it wall-surfaces is factored in the graph as:

$c_{\boldsymbol{\omega}}(\leftidx{{}^{B}}{\boldsymbol{\omega}_{i}},\big{[}\leftidx{{}^{B}}{\boldsymbol{\pi}_{x_{1}}},\leftidx{{}^{B}}{\boldsymbol{\pi}_{x_{1}}},\leftidx{{}^{B}}{\mathfrak{s}}_{i}])\\\
=\sum_{i=1}^{K}\|\leftidx{{}^{B}}{\hat{\boldsymbol{\omega}}_{i}}-f(\leftidx{{}^{B}}{\tilde{\boldsymbol{\pi}}_{x_{1}}},\leftidx{{}^{B}}{\tilde{\boldsymbol{\pi}}_{x_{1}}},\leftidx{{}^{B}}{\mathfrak{s}_{i}})\|^{2}_{\mathbf{\Lambda}_{\boldsymbol{\tilde{\boldsymbol{\omega}}}_{i,t}}}$
(2)

Where
$f(\leftidx{{}^{B}}{\tilde{\boldsymbol{\pi}}_{x_{1}}},\leftidx{{}^{B}}{\tilde{\boldsymbol{\pi}}_{x_{1}}},\leftidx{{}^{B}}{\mathfrak{s}_{i}})$
is the function mapping the wall center using the wall-surfaces and its
starting point following Eq. III. Wall factors add an additional layer of
structural consistency to the graph. A Wall center for opposed planes in
$y$-direction is computed following Eq. III.

A-Rooms Layer. The second layer of the graph extracts all the information
about the rooms along with the door-ways interconnecting the rooms.

Rooms: We use the similar concept of a four-wall room
$\leftidx{{}^{B}}{\boldsymbol{\rho}}\in\mathbb{R}^{2}$ as presented in [6],
where each room comprises the four-wall surfaces extracted in the first layer
of the graph.

Door-Ways: We incorporate an additional entity in the graph called door-ways
interconnecting the room nodes, easily available from BIM. The position of a
door-way node $\leftidx{{}^{B}}{\boldsymbol{\mathcal{{D}}}}\in\mathbb{R}^{3}$
is directly extracted from BIM in the frame of reference $B$. Using the
semantic information from BIM of the rooms connected by a given door-way, the
door-way-to-rooms factor can be formulated as:

$c_{\boldsymbol{\mathcal{D}}}(\leftidx{{}^{B}}{\boldsymbol{\rho}_{1}},\leftidx{{}^{B}}{\boldsymbol{\rho}_{2}},\leftidx{{}^{B}}{\boldsymbol{\mathcal{D}}_{i}})=\\\
\|\
f(\leftidx{{}^{B}}{\hat{\boldsymbol{\rho}}_{1}},\leftidx{{}^{\rho_{1}}}{\hat{\boldsymbol{{\mathcal{D}}}}_{i}})-\
f(\leftidx{{}^{B}}{\hat{\boldsymbol{\rho}}_{2}},\leftidx{{}^{\rho_{2}}}{\hat{\boldsymbol{\mathcal{D}}}_{i}})\|$
(3)

Where $\leftidx{{}^{B}}{\boldsymbol{\rho}_{1}}$ and
$\leftidx{{}^{B}}{\boldsymbol{\rho}_{2}}$ are the four-wall rooms connected to
the door-way $\leftidx{{}^{B}}{\boldsymbol{\mathcal{D}}_{i}}$.
$\leftidx{{}^{\rho_{1}}}{\boldsymbol{\mathcal{D}}_{i}}$ and
$\leftidx{{}^{\rho_{2}}}{\boldsymbol{\mathcal{D}}_{i}}$ are the positions the
door-way nodes estimated with respect to rooms
$\leftidx{{}^{B}}{\boldsymbol{\rho}_{1}}$ and
$\leftidx{{}^{B}}{\boldsymbol{\rho}_{2}}$.

##  IV Graph Matching

Our second contribution to this paper is a novel approach to graph matching in
which we match an architectural A-graph, $\mathcal{G}_{a}$, and a S-Graph,
$\mathcal{G}_{s}$. We compare and match the room ($\rho$) entities and for
each room entity its corresponding wall-surfaces ($\pi$), extracted from both
graphs. The correspondences form a bipartite graph connecting the nodes of
some parts of $\mathcal{G}_{a}$ with all, or almost all nodes, in
$\mathcal{G}_{s}$, at the rooms and walls layers. $\mathcal{G}_{s}$ is built
incrementally as the robot navigates the environment, the graph matching is
run after every map update until a successful match is obtained. A schema of
the entire graph-matching process is described in Fig. 2.

Notation. Let $\boldsymbol{V}$ be any set of nodes, with $\boldsymbol{V}_{a}$
and $\boldsymbol{V}_{s}$ the sets of nodes in the A-Graph and the S-Graph
respectively. Let $\boldsymbol{m}={(v_{a},v_{s}):v_{a}\in V_{a},v_{s}\in
V_{s}}$, and $\boldsymbol{M}$ be any set of $m$ such as
$\boldsymbol{M}=\\{\boldsymbol{m}_{1},\boldsymbol{m}_{2},...,\boldsymbol{m}_{n}\\}$.
Let $\boldsymbol{\mathcal{M}}$ be any set of $\boldsymbol{M}$ such as
$\boldsymbol{\mathcal{M}}=\\{\boldsymbol{M}_{1},\boldsymbol{M}_{2},...,\boldsymbol{M}_{n}\\}$.
Local candidates are $\boldsymbol{M}$ including $\boldsymbol{m}$ referring to
a small part of the input graphs. Global candidates are $\boldsymbol{M}$
including $\boldsymbol{m}$ referring to the entire input graphs.

Figure 2: Graph matching schema. a) Downwards and at each level, different
combinations of matches are proposed and selected by their geometrical
affinity either at the same level or with the associated upper-level pair. b)
The match graph is traversed upwards while combined with same-level nodes to
define all-level match candidates. c) The lowest-level pairs of every
candidate are scored in global affinity. That score is clustered to find
symmetries in the best cluster.

##  V Graph Merging

Our global state $\mathbf{s}$ at time $T$, before graph merging, contains all
the nodes of the A-Graph, generated offline, as well as the current nodes
estimated online by the S-Graph

$\displaystyle\mathbf{s}$ $\displaystyle=$
$\displaystyle[\leftidx{{}^{M}}{\mathbf{x}}_{R_{1}},\ \ldots,\
\leftidx{{}^{M}}{\mathbf{x}}_{R_{T}},$
$\displaystyle\leftidx{{}^{M}}{\boldsymbol{\pi}}_{1},\ \ldots,\
\leftidx{{}^{M}}{\boldsymbol{\pi}}_{P},\
\leftidx{{}^{B}}{\boldsymbol{\pi}}_{1},\ \ldots,\
\leftidx{{}^{B}}{\boldsymbol{\pi}}_{Q},$
$\displaystyle\leftidx{{}^{M}}{\boldsymbol{\rho}}_{1},\ \ldots,\
\leftidx{{}^{M}}{\boldsymbol{\rho}}_{S},\
\leftidx{{}^{B}}{\boldsymbol{\rho}}_{1},\ \ldots,\
\leftidx{{}^{B}}{\boldsymbol{\rho}}_{R}$
$\displaystyle\leftidx{{}^{M}}{\boldsymbol{\gamma}}_{1},\ \ldots,\
\leftidx{{}^{M}}{\boldsymbol{\gamma}}_{G},\leftidx{{}^{M}}{\boldsymbol{\xi}}_{1},\
\ldots,\ \leftidx{{}^{M}}{\boldsymbol{\xi}}_{E},$
$\displaystyle\leftidx{{}^{B}}{\boldsymbol{\omega}_{1}},\ \ldots,\
\leftidx{{}^{B}}{\boldsymbol{\omega}_{W}},\
\leftidx{{}^{B}}{\boldsymbol{\mathcal{D}}}_{1},\ \ldots,\
\leftidx{{}^{B}}{\boldsymbol{\mathcal{D}}}_{D},$
$\displaystyle\leftidx{{}^{B}}{\mathbf{x}}_{M}]^{\top},$

where $\leftidx{{}^{B}}{\mathbf{x}}_{M}$ is the estimated transformation
between the map frame $M$ of the S-Graph and the BIM frame B of the A-Graph,
which is set to identity before graph merging. The graph matching method from
Section. IV provides match candidates between the room nodes and the wall-
surface nodes of the S-Graph and the A-Graph. To efficiently merge the two
graphs to generate the iS-Graph, we introduce room-to-room constraints as well
as a wall-surface-to-wall-surface constraints between the matched candidates.
The room-to-room constraint is defined as

$c_{\boldsymbol{\rho}}(\leftidx{{}^{B}}{\boldsymbol{\rho}_{1}},\leftidx{{}^{M}}{\boldsymbol{\rho}_{2}})=\|\
\leftidx{{}^{B}}{\hat{\boldsymbol{\rho}}_{{1}}}-\leftidx{{}^{M}}{\hat{\boldsymbol{\rho}}_{{2}}}\|^{2}_{\mathbf{\Lambda}_{\boldsymbol{\tilde{\boldsymbol{\rho}}}_{1,2}}}~{},$
(5)

where $\leftidx{{}^{B}}{\boldsymbol{\rho}_{1}}$ is the room node in the
A-Graph and $\leftidx{{}^{M}}{\boldsymbol{\rho}_{2}}$ is the corresponding
room node in S-Graph. Similarly, for all correspondences between wall-surface
candidates, the constraint is formulated as

$c_{\boldsymbol{\pi}}(\leftidx{{}^{B}}{\boldsymbol{\pi}_{1}},\leftidx{{}^{M}}{\boldsymbol{\pi}_{2}})=\|\
\leftidx{{}^{B}}{\boldsymbol{\pi}_{1}}-\leftidx{{}^{M}}{\boldsymbol{\pi}_{2}}\|^{2}_{\mathbf{\Lambda}_{\boldsymbol{\tilde{\boldsymbol{\pi}}}_{1,2}}}~{},$
(6)

where $\leftidx{{}^{B}}{\boldsymbol{\pi}_{1}}$ and
$\leftidx{{}^{M}}{\boldsymbol{\pi}_{2}}$ are the wall-surfaces in the A-Graph
and the S-Graph.

With the constraints between the two graphs,
$\leftidx{{}^{B}}{\mathbf{x}}_{M}$ can be estimated and all the robot poses,
the wall-surfaces, rooms, and floors of the _S-Graph_ can be referred
accurately with respect to the BIM frame of reference $B$ of the _A-Graph_ ,
resulting in the final improved situational graph iS-Graph. In this manner the
robot is localized with respect to the global reference of the architectural
plan. Our approach thus can perform global localization exploiting the
hierarchical high-level information in the environment without the need for
appearance-based loop closure constraints at keyframe level, which are more
variable.

##  VI Experimental Evaluation

###  VI-A Experimental Setup

We validated our proposed approach in multiple simulated and real-world
construction environments. In all experiments, we generated A-Graphs from
various building models that were created in Autodesk Revit. The Boston
Dynamics Spot robot, equipped with a Velodyne VLP-16 3D LiDAR, was
teleoperated for data collection. We compared our approach to two 2D LiDAR-
based localization algorithms (AMCL [7] and Cartographer [8]) and one 3D
LiDAR-based localization (UKFL [9]) algorithm. In simulated datasets we
measured the Absolute Pose Error (APE) [10] with respect to the ground truth,
while in real-world due to the absence of ground truth pose information, we
compare the point cloud RMSE of the generated map with the available ground
truth map from BIM. The proposed methodology was implemented in C++, and the
experiments were validated on an Intel i9 16-core workstation.

TABLE I: Absolute Pose Error (APE) [m] for several LiDAR-based localization baselines and our iS-Graphs Datasets have been recorded in simulated environments. ‘$-$’ stands for localization failure. Method |  | APE [m] $\downarrow$ |
---|---|---|---
|  | Datasets |
| D1 | D2 | D3 | D4 | D5 | D6
AMCL [7] | 2.04 | 1.71 | 2.03 | 2.01 | - | -
UKFL [9] | 0.97 | 0.78 | - | 0.74 | 0.70 | 0.88
Cartographer [8] | 0.10 | 0.16 | 0.12 | - | 0.13 | -
iS-Graphs (ours) | 0.08 | 0.01 | 0.02 | 0.04 | 0.09 | 0.12
TABLE II: Point cloud RMSE [m] for our real-world dataset. Best results are boldfaced. ‘$-$’ stands for localization failure. Method | Alignment Error $\downarrow$ |  |
---|---|---|---
|  | Datasets |  |  |
| D1 | D2 | D3 |  |
AMCL [7] | - | 0.90 | 0.98 |  |
UKFL [9] | - | 0.86 | 0.69 |  |
Cartographer [8] | - | 0.58 | 0.64 |  |
iS-Graphs (ours) | 0.17 | 0.20 | 0.21 |  |

###  VI-B Results and Discussion

Simulated Datasets. Table. I presents the APE of our proposed iS-Graphs and
state-of-the LiDAR-based localization algorithms for six simulated datasets.
Each simulated dataset represents a single floor level of a given construction
environment with varying configurations of wall-surfaces and rooms. Given the
output from 2D LiDAR algorithms, the APE of all the algorithms is computed in
2D ($x$, $y$, $\theta$). Table. I shows that our iS-Graphs shows higher
robustness against localization failure and outperforms all the localization
baselines using both 2D and 3D LiDARs.

Real-World Datasets. Table. II presents the point cloud RMSE for three
different construction environments, where our iS-Graphs is able to localize
the robot correctly while providing a more accurate 3D map of the environment
as compared to the ground truth. Because of noisy LiDAR measurements and the
clutter in a real construction environment traditional approaches to
localization fail to localize the robot (see D1 in Table. II). Since our
approach relies only on high-level entities like wall-surfaces and rooms and
their topological relationship instead of directly relying on low-level LiDAR
measurements, it is more robust to noise and clutter in the environment than
the baselines.

##  VII Conclusion

In this paper, we presented a novel method for global robot localization
utilizing prior information from architectural plans. We embeded the
architectural data from BIM models into optimizable graphs called
Architectural Graphs (A-Graphs). Using Situational Graphs (S-Graphs) estimated
by a robot as it navigates its environment, we present a novel graph matching
strategy to match the A-Graphs with the S-Graphs and we also present a graph
merging strategy to fuse both graphs. The result of the fusion is an informed
iS-Graph, which enables the robot to localize itself within the architectural
plan for the given environment.

## References

  * [1] S. Azhar, M. Hein, and B. Sketo, “Building information modeling (bim): benefits, risks and challenges,” in _Proceedings of the 44th ASC Annual Conference_ , 2008, pp. 2–5.
  * [2] R. W. M. Hendrikx, P. Pauwels, E. Torta, H. J. Bruyninckx, and M. J. G. van de Molengraft, “Connecting semantic building information models and robotics: An application to 2d lidar-based localization,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 11 654–11 660.
  * [3] M. S. Moura, C. Rizzo, and D. Serrano, “Bim-based localization and mapping for mobile robots in construction,” in _2021 IEEE International Conference on Autonomous Robot Systems and Competitions (ICARSC)_ , 2021, pp. 12–18.
  * [4] M. Shaheer, H. Bavle, J. L. Sanchez-Lopez, and H. Voos, “Robot localization using situational graphs and building architectural plans,” 2022. [Online]. Available: https://arxiv.org/abs/2209.11575
  * [5] H. Bavle, J. L. Sanchez-Lopez, M. Shaheer, J. Civera, and H. Voos, “Situational graphs for robot navigation in structured indoor environments,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 4, pp. 9107–9114, 2022.
  * [6] ——, “S-graphs+: Real-time localization and mapping leveraging hierarchical representations,” 2022. [Online]. Available: https://arxiv.org/abs/2212.11770
  * [7] D. Fox, “Kld-sampling: Adaptive particle filters,” _Advances in neural information processing systems_ , vol. 14, 2001.
  * [8] W. Hess, D. Kohler, H. Rapp, and D. Andor, “Real-time loop closure in 2d lidar slam,” in _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , 2016, pp. 1271–1278.
  * [9] K. Koide, J. Miura, and E. Menegatti, “A portable three-dimensional LIDAR-based system for long-term and wide-area people behavior measurement,” _International Journal of Advanced Robotic Systems_ , vol. 16, no. 2, Mar. 2019.
  * [10] M. Grupp, “evo: Python package for the evaluation of odometry and SLAM,” https://github.com/MichaelGrupp/evo, 2017.


# Virtual Linking Bids for Market Clearing with Non-Merchant Storage

Eléa Prat, Jonas Bodulv Broge, Richard Lusby  Technical University of Denmark,
Department of Technology, Management and Economics
{emapr<EMAIL_ADDRESS>
###### Abstract

In the context of energy market clearing, non-merchant assets are assets that
do not submit bids but whose operational constraints are included. Integrating
energy storage systems as non-merchant assets can maximize social welfare.
However, the disconnection between market intervals poses challenges for
market properties, that are not well-considered yet. We contribute to the
literature on market-clearing with non-merchant storage by proposing a market-
clearing procedure that preserves desirable market properties, even under
uncertainty. This approach is based on a novel representation of the storage
system in which the energy available is discretized to reflect the different
prices at which the storage system was charged. These prices are included as
virtual bids in the market clearing, establishing a link between different
market intervals. We show that market clearing with virtual linking bids
outperforms traditional methods in terms of cost recovery for the market
participants and discuss the impacts on social welfare.

###### Index Terms:

energy market design, non-merchant storage, passive storage

Thanks to Linde Frölke, Pierre Pinson and Jalal Kazempour for the first
discussions on the concept presented in this paper. Thanks to Pierre Pinson
for reviewing and providing feedback.

##  I Introduction

In order to enable the safe operation of future energy systems with a high
share of intermittent and stochastic renewable sources of energy production,
the share of large-scale energy storage is expected to increase significantly
in the coming years [1]. A major challenge that comes with this evolution is
how to best integrate storage in energy markets. This is emphasized by a
recent order by the Federal Energy Regulatory Commission in the United States,
urging system operators to implement changes in order to facilitate market
participation of electric storage systems [2]. As a step towards addressing
this challenge, [3] describes three different ways in which storage systems
can be included in the market clearing. In the first and second options,
storage systems participate similarly to conventional generators and loads,
submitting price and quantity bids. We refer to this setup as _merchant
storage_. In the third option, the storage systems’ operational constraints
are included but they do not need to submit bids, which we call _non-merchant
storage_. Market clearing with non-merchant storage can achieve the most
economically efficient outcomes and the highest social welfare [3], as opposed
to market clearing with merchant storage [4, 5]. However, this setup requires
more modifications to the current energy markets, which might explain why it
has not been considered in detail yet.

The main issue to address when clearing a market with non-merchant storage is
how to represent the time-linking effect of the storage system. Indeed, while
storage systems potentially connect an infinite number of time periods, the
market-clearing window is finite. In the literature on non-merchant storage, a
very common assumption is that the storage is initially empty and the state of
energy at the end of the market-clearing horizon is free111This is equivalent
to setting it to zero in the absence of negative prices. [6, 7, 8]. This leads
to myopic decision-making regarding the state of energy of the storage at the
end of the market interval. These assumptions are not always stated [9],
showing that this problem is disregarded.

This time-linking property also has to be considered in the pricing problem.
By default, the market does not remember the price at which the storage was
charged and thus considers any energy available at the beginning of a market
horizon as free. Therefore, the resulting market price might be lower than the
price the storage system paid to charge. In [10] we proposed a method to
reestablish the connection between the different market intervals and retrieve
prices that send the proper signal to the storage system. However, it is only
valid under the assumption of perfect information, which does not correspond
to a realistic setup. The problem of connection between current and future
time periods with uncertain realization has also been identified in the case
of remuneration of the storage system with financial storage rights [11]. This
situation has also been observed in the case of demand response in [12], which
could be seen as a virtual type of storage in the case of load shifting [13],
and is therefore also covered by this paper.

In real-time markets, which include ramping products that can span over
several time periods, a common approach is to use a multi-interval market
clearing, where only the decisions over the first time periods, termed the
_decision horizon_ , are implemented, while the decisions over the rest of the
horizon are advisory [14, 15, 16, 17]. It can help to make future-aware
decisions regarding the state of energy at the end of the decision horizon.
However, even in this framework, the final decision might have an impact if
the window is not long enough. Moreover, it also suffers from issues with
pricing, even in a deterministic setup, as shown in [16]. The work in [16]
proposes two methods to reintroduce a link to previous intervals but does not
prevent the need for non-transparent uplift payments. In [14], an argument for
non-uniform pricing is made, but the prices they obtain are highly dependent
on the forecasts used.

Because of these limitations, we focus on single-interval market clearing with
uniform pricing. Moreover, this is the usual approach for day-ahead auctions,
while multi-interval market-clearing models are used for real-time markets
[17]. Since single-interval markets are still widely in use, it is important
to study the integration of non-merchant storage into those. We move the
consideration of future intervals to the choice of end-of-horizon storage
parameters, which are to be decided in a separate problem222This other problem
is out-of-scope for this paper.. In this setup, a modeling approach to reflect
the transfer of energy between time periods has been introduced in [18].
However, the authors do not discuss the valuation of the charged energy in
subsequent market intervals. The problem of valuation has been studied in [12]
for demand response, in terms of deviation from a scenario with no
flexibility, which cannot always be applied in the case of a storage system.

We aim to help address the unanswered question: How to best design a market
clearing that includes non-merchant storage and fully exploit the potential of
storage systems to increase social welfare? Towards this, we propose a novel
market-clearing procedure with non-merchant storage that ensures cost recovery
for all the market participants, in particular for the storage system, and
even when considering uncertainty. The main idea is to remember the prices at
which the storage system charged and automatically create virtual linking bids
to reflect these values in the next market intervals.

First, in Section II, we present in more detail the challenges of market
clearing with non-merchant storage, with the help of an illustrative example.
We then introduce market clearing with virtual linking bids in Section III and
show promising results in terms of market properties in Section IV. We
conclude in Section V.

##  II Challenges of Market Clearing with Non-Merchant Storage

###  II-A Initial Model and Assumptions

In order to get a better understanding of the problem, we model a stylized
version of the storage system using a couple of assumptions. First, we
consider that there is only one, non-merchant, storage system in the market.
We also assume that there are no losses when charging and discharging, meaning
that the storage system is perfectly efficient, and has no leakage over time.
We consider a storage system without charging and discharging limits and with
no minimum on the state of energy. We thereby focus on the time-linking aspect
of the storage system and limit the subtleties that would be introduced by
considering each of these aspects. Indeed, the essence of the method presented
in this paper would be the same, but small modifications would have to be
introduced to deal with these different aspects, which would complicate the
understanding of the basics if included here. We furthermore assume that
minimum levels for loads and generators are 0, in order to avoid non-
linearities.

We refer to _split market clearing_ as the process of clearing the market
interval-by-interval, as opposed to an _ideal market clearing_ , where all
market intervals would be cleared at once. Here, we consider one market
interval. The set $\mathcal{T}$ gathers all the time periods $t$ of the given
market interval, where each time period has a duration of $\Delta t$ (in
hours). For a day-ahead market for instance, $\mathcal{T}$ would typically
correspond to one day and $\Delta t$ would be one hour.

Under these assumptions, the market clearing with non-merchant storage for a
given market interval can be modeled as follows:

$\displaystyle\max_{\mathbf{x}}\quad$ $\displaystyle\Delta
t\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}\right)$
(1a) s.t.
$\displaystyle\sum_{l\in\mathcal{L}}d_{lt}+p_{t}^{\mathrm{C}}-\sum_{g\in\mathcal{G}}p_{gt}=0,\hskip
112.73882pt\forall t\in\mathcal{T}$ (1b) $\displaystyle 0\leq
p_{gt}\leq\overline{P}_{gt},\hskip 151.76964pt\forall
g\in\mathcal{G},t\in\mathcal{T}$ (1c) $\displaystyle 0\leq
d_{lt}\leq\overline{D}_{lt},\hskip 156.10345pt\forall
l\in\mathcal{L},t\in\mathcal{T}$ (1d) $\displaystyle 0\leq
e_{t}\leq\overline{E},\hskip 208.13574pt\forall t\in\mathcal{T}$ (1e)
$\displaystyle e_{t}=e_{t-1}+p_{t}^{\mathrm{C}}\Delta t,\hskip
134.42113pt\forall t\in\mathcal{T}\setminus\\{1\\}$ (1f) $\displaystyle
e_{1}=E^{\mathrm{init}}+p_{1}^{\mathrm{C}}\Delta t.$ (1g)

Here, and in the following, $\mathbf{x}$ is a vector gathering all the
decision variables of the model at hand. The variables are the quantity
accepted for load $l\in\mathcal{L}$, $d_{lt}$, the quantity accepted for
generator $g\in\mathcal{G}$, $p_{gt}$, and for the storage system, the state
of energy $e_{t}$ and the amount charged $p_{t}^{\mathrm{C}}$. The latter can
be negative, indicating a discharge. The objective function (1a) is to
maximize the difference between load utilities $U_{lt}$ and generation costs
$C_{gt}$ for accepted offers. The maximum generation $\overline{P}_{gt}$ and
load $\overline{D}_{lt}$ are enforced by constraints (1c) and (1d)
respectively. Constraint (1e) sets the maximum state of energy $\overline{E}$.
Constraints (1f) and (1g) update the storage level, starting from the initial
level $E^{\mathrm{init}}$. The market price at $t$ is given by $\lambda_{t}$,
the dual variable of (1b).

We make further assumptions regarding the market operation, namely that there
is perfect competition and that market participants bid their true costs.

To evaluate the efficiency of the market, one of the results to consider is
the social welfare, $SW$, which is calculated as the sum of surpluses of loads
and generators and the payments to the storage system. With (1b), it reduces
to:

$\small SW=\Delta
t\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}^{*}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}^{*}\right),$
(2)

where the superscript ∗ indicates that the optimal value of the considered
variables is used.

###  II-B Approaches to Avoid Myopic Decisions

As mentioned in [3] and [10], the model in (1) is not sufficient to avoid
myopic decisions regarding the state of energy at the end of the market
interval. Indeed, with this model the storage level will most likely return to
zero, and there will not be stored energy available in the first time period
of the next market interval. In [3], several options to avoid myopic decisions
are listed. The first one is to impose a final state of energy
$E^{\mathrm{end}}$ at the end of the market interval, which is determined by
considering information about future market intervals. This is done by adding
to (1) the constraint

$\small e_{\text{T}}=E^{\mathrm{end}},$ (3)

where $t=\text{T}$ corresponds to the last time period of the market interval.
Another option is to steer the level to the desired value by adding a penalty
term in the objective function, with a cost $S^{\mathrm{end}}$. In (1), the
objective function (1a) becomes

$\small\max_{\mathbf{x}}\quad\Delta
t\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}\right)-S^{\mathrm{end}}e_{\text{T}}.$
(4)

###  II-C Limits on an Illustrative Example

We consider a storage system that has a capacity $\overline{E}=2.5$ MWh and
which is initially empty. We clear the market for two market intervals (MI) of
one time period of one hour each. One load and two generators participate in
the market. The related parameters are listed in Table I. The code for this
example, as well as all the other examples presented in this paper, is
available online at https://github.com/eleaprat/MC_non_merchant_stg.

TABLE I: Load and generators parameters for the illustrative example. Prices are in €/MWh and quantities in MW. MI | $t$ | $U_{1t}$ | $\overline{D}_{1t}$ | $C_{1t}$ | $\overline{P}_{1t}$ | $C_{2t}$ | $\overline{P}_{2t}$
---|---|---|---|---|---|---|---
1 | 1 | 12 | 0 | 5 | 2 | 10 | 2
2 | 1 | 12 | 3 | 2 | 2 | 9 | 2

We first solve (1) with (3), where $E^{\mathrm{end}}=1$ MWh for the first MI
and $E^{\mathrm{end}}=0$ MWh for the second MI 333These are found to maximize
the total social welfare when clearing the two MIs together. We refer the
interested reader to [10] for more information.. The resulting dispatch and
market prices are shown in Table II.

TABLE II: Results for the illustrative example. The state of energy is in MWh, the market price in €/MWh and the rest of the variables in MW. | MI | $t$ | $e_{t}$ | $d_{1t}$ | $p_{1t}$ | $p_{2t}$ | Market price
---|---|---|---|---|---|---|---
With (3) | 1 | 1 | 1 | 0 | 1 | 0 | 5
2 | 1 | 0 | 3 | 2 | 0 | [2,9]
With (4) | 1 | 1 | 0 | 0 | 0 | 0 | 5
2 | 1 | 0 | 3 | 2 | 1 | 9

We can see that on MI 2, there is a price multiplicity, where any price
between 2 and 9€/MWh is valid. It can be a problem if the final price chosen
is below 5€/MWh because the storage system would then not recover its charging
cost. On MI 2, the information that the storage system charged at 5€/MWh is
not accessible. As a consequence, the price is chosen without accounting for
it. The total social welfare is 27€.

On the other hand, we can use the penalty term by solving (1) with (4).
Knowing about the potential price multiplicity on MI 2, we set the penalty
$S^{\mathrm{end}}=2$€/MWh on the first MI, to make sure the storage system
will then recover its costs if this lower price gets selected on MI 2, and
$S^{\mathrm{end}}=0$ on MI 2. The resulting dispatch and market prices are
shown in Table II. In this case, the storage system does not charge at all,
and as a consequence, the social welfare is reduced to 23€. Due to this
limitation, in the rest of the paper ”split market clearing” refers to solving
(1) with (3).

##  III Market Clearing with Virtual Linking Bids

In this section, we modify the market-clearing model with final storage level,
to ensure that the cost at which the storage system charged in one market
interval is accounted for in the following market intervals. To do so, we
introduce a new representation of the non-merchant storage system.

###  III-A Inter- and Intra-Storage

In order to save the value at which the storage charges in view of subsequent
intervals, we introduce the concepts of net charge and net discharge over a
market interval. Those are indicated by the difference between the initial and
the final state of energy. If it is positive, it corresponds to a net charge
and if negative, it corresponds to a net discharge. We can then consider
separately the exchanges of energy within the market interval, which we
associate to an _intra-storage_ , and the exchanges of energy with past or
future intervals, which correspond to _inter-storage_. We can conceptually
split the storage, where the quantity charged previously is equivalent to a
generator, bidding with the saved price, and the capacity for intra-storage
corresponds to the available capacity at the beginning of the market interval.
This is illustrated in Figure 1.

Figure 1: Storage system at the beginning of a market interval, with the
previously charged quantities in red and the associated price in parenthesis.

In the case of net charge, this quantity is saved along with the corresponding
charging price, which is then used as a virtual bid in future market
intervals. In this way, we ensure that the storage will get paid at least what
it paid for charging. Examples of net charge and net discharge are shown in
Figure 2. In the example of net charge from Figure 2a, we see that over the
market interval, the quantity in the inter-storage is not used. At the end of
the market interval, the quantity net charged is added to the intra-storage.
In the example of net discharge from Figure 2b, we can see how intra- and
inter-storages can be used independently. Similarly, the level at the end of
the market interval determines the intra-storage for the next interval.

(a) Example of net charge

(b) Example of net discharge

Figure 2: Examples of net charge and net discharge, for a fictional market
interval of 3 time periods of one hour each. Quantities are in MW.

###  III-B Model for Market Clearing with Virtual Linking Bids

We now introduce a model for market clearing with virtual linking bids (VLB),
based on the model in (1) and including the representation of the storage
system with intra- and inter-storage components.

$\displaystyle\max_{\mathbf{x}}\quad$ $\displaystyle\Delta
t\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}-\sum_{v\in\mathcal{V}}S_{v}p_{vt}^{\mathrm{D,e}}\right)$
(5a) s.t.
$\displaystyle\sum_{l\in\mathcal{L}}d_{lt}+p_{t}^{\mathrm{C,a}}-\sum_{g\in\mathcal{G}}p_{gt}-\sum_{v\in\mathcal{V}}p_{vt}^{\mathrm{D,e}}=0,\hskip
26.01613pt\forall t\in\mathcal{T}$ (5b)
$\displaystyle\eqref{eq:gen_bound}-\eqref{eq:load_bound}$ $\displaystyle
e^{\mathrm{a}}_{t}=e^{\mathrm{a}}_{t-1}+p_{t}^{\mathrm{C,a}}\Delta t,\hskip
121.41306pt\forall t\in\mathcal{T}\setminus\\{1\\}$ (5c) $\displaystyle
e^{\mathrm{a}}_{1}=p_{1}^{\mathrm{C,a}}\Delta t,$ (5d) $\displaystyle
e^{\mathrm{a}}_{\text{T}}\geq 0,$ (5e) $\displaystyle
e^{\mathrm{e}}_{vt}=e^{\mathrm{e}}_{v,t-1}-p_{vt}^{\mathrm{D,e}}\Delta
t,\hskip 56.37271pt\forall v\in\mathcal{V},\,t\in\mathcal{T}\setminus\\{1\\}$
(5f) $\displaystyle
e^{\mathrm{e}}_{v,1}=E^{\mathrm{init}}_{v}-p_{v,1}^{\mathrm{D,e}}\Delta
t,\hskip 143.09538pt\forall v\in\mathcal{V}$ (5g) $\displaystyle 0\leq
e^{\mathrm{a}}_{t}+\sum_{v\in\mathcal{V}}e^{\mathrm{e}}_{vt}\leq\overline{E},\hskip
147.4292pt\forall t\in\mathcal{T}$ (5h) $\displaystyle
e^{\mathrm{e}}_{vt},\,p_{vt}^{\mathrm{D,e}}\geq 0,\hskip 151.76964pt\forall
v\in\mathcal{V},\,\forall t\in\mathcal{T}$ (5i) $\displaystyle
e^{\mathrm{a}}_{\text{T}}+\sum_{v\in\mathcal{V}}e^{\mathrm{e}}_{v\text{T}}\geq
E^{\mathrm{end}}.$ (5j)

Here again, $\mathbf{x}$ is a vector gathering all the decision variables of
the model. There are no changes regarding loads and generators and their
decision variables. We have new variables for tracking the storage system. For
the intra-storage system, $e^{\mathrm{a}}_{t}$ gives the state of energy and
$p_{t}^{\mathrm{C,a}}$ the quantity charged (negative for a discharge). For
the inter-storage system, we introduce $v\in\mathcal{V}$, which corresponds to
the different values saved in the storage system, similarly to what was shown
in Figure 1. The corresponding variables are $e^{\mathrm{e}}_{vt}$, the state
of energy for value $v$ and $p_{vt}^{\mathrm{D,e}}$, the quantity discharged
from inter-storage with value $v$. The objective function (5a) is modified to
include the artificial bids from the inter-storage, with prices $S_{v}$. Here
and in the balance constraint (5b), the inter-storage appears similarly to
conventional generators. Constraints (5c) and (5d) give the update of the
state of energy for the intra-storage. Note that the intra-storage is by
definition empty at the beginning of the market interval. However, it can
occasionally take negative values, which corresponds to temporarily using some
of the capacity that is reserved by the inter-storage. This is necessary to
ensure that the arbitrage opportunities within a given market interval are not
limited due to the intra-storage being initially empty, while there is, in
practice, some energy to discharge. To limit these operations to arbitrage
within this market interval only, (5e) specifies that the final state of
energy of the intra-storage cannot be negative.

Constraints (5f) and (5g) give the update of the state of energy for each
inter-storage. The initial energy available at each value is given by
$E^{\mathrm{init}}_{v}$. Constraint (5h) enforces the storage capacity. The
states of energy and the quantity discharged from each inter-storage are
positive with (5i).

Constraint (5j) is the updated version of constraint (3). With the greater or
equal sign, the artificial bid in the objective function guides the final
level of the storage. This formulation thus combines some of the aspects of
(3) and (4).

There is no charge of the inter-storage during the market clearing. In the
case of net charge, the inter-storage is modified outside of the market
clearing, as shown in Section III-D.

###  III-C Simultaneous Charge and Discharge

Note that in the model in (5), $p_{t}^{\mathrm{C,a}}$ and
$p_{vt}^{\mathrm{D,e}}$ can be positive at the same time. However, this
separation of the storage system is only a modeling artifice and has no
physical meaning. It is the sum of the two that corresponds to the instruction
of charge or discharge for the storage system. Allowing the intra-storage
system to take negative values generates a multiplicity of optimal solutions,
for which the sum of intra- and inter-storage is the same. It is possible to
select among all those solutions one for which there is no simultaneous charge
and discharge, by minimizing the total quantity exchanged or by introducing
binary variables to choose between charge and discharge. We formalize this
result with a proposition.

######  Proposition 1.

For the market clearing in (5), and under the assumptions that the value of
stored energy is strictly positive, i.e., $S_{v}>0$, $\forall
v\in\mathcal{V}$, and that a feasible solution exists, there will always be an
optimal solution where the intra-storage does not charge when the inter-
storage discharges.

###### Proof.

First, we give an equivalent formulation of the model where the state of
charge is replaced, using that
$e^{\mathrm{a}}_{t}=\sum_{i=1}^{t}p_{i}^{\mathrm{C,a}}$ and
$e^{\mathrm{e}}_{vt}=E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e}}$.
We also notice that in (5i), $e^{\mathrm{e}}_{vt}\geq 0$ can be equivalently
replaced by $e^{\mathrm{e}}_{vT}\geq 0$, since the update quantities
$p_{vt}^{\mathrm{D,e}}$ are non-negative. Lastly, we assume that the time
periods of the market have a duration of one hour in order to lighten the
notations with $\Delta t=1$. The same proof can be made for any value of
$\Delta t$. We obtain the following:

$\displaystyle\max_{\mathbf{x}}\quad$
$\displaystyle\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}-\sum_{v\in\mathcal{V}}S_{v}p_{vt}^{\mathrm{D,e}}\right)$
(6a) s.t.
$\displaystyle\sum_{l\in\mathcal{L}}d_{lt}+p_{t}^{\mathrm{C,a}}-\sum_{g\in\mathcal{G}}p_{gt}-\sum_{v\in\mathcal{V}}p_{vt}^{\mathrm{D,e}}=0,\hskip
26.01613pt\forall t\in\mathcal{T}$ (6b) $\displaystyle 0\leq
p_{gt}\leq\overline{P}_{gt},\hskip 151.76964pt\forall
g\in\mathcal{G},t\in\mathcal{T}$ (6c) $\displaystyle 0\leq
d_{lt}\leq\overline{D}_{lt},\hskip 156.10345pt\forall
l\in\mathcal{L},t\in\mathcal{T}$ (6d)
$\displaystyle\sum_{i=1}^{T}p_{i}^{\mathrm{C,a}}\geq 0,$ (6e) $\displaystyle
0\leq\sum_{i=1}^{t}p_{i}^{\mathrm{C,a}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e}}\right)\leq\overline{E},\hskip
0.43668pt\forall t\in\mathcal{T}$ (6f) $\displaystyle
p_{vt}^{\mathrm{D,e}}\geq 0,\hskip 177.78578pt\forall
v\in\mathcal{V},\,\forall t\in\mathcal{T}$ (6g) $\displaystyle
E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e}}\geq 0,\hskip
156.10345pt\forall v\in\mathcal{V}$ (6h)
$\displaystyle\sum_{i=1}^{T}p_{i}^{\mathrm{C,a}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e}}\right)\geq
E^{\mathrm{end}}.$ (6i)

Let’s consider an optimal solution to (6). We denote it with the superscript
∗, for example, $p_{t}^{\mathrm{C,a*}}$. Under the assumption that the
feasible set is not empty, such a solution exists. If it is such that the
intra-storage never charges when the inter-storage discharges, we are done. We
consider the case where there is at least one time period $\tau\in\mathcal{T}$
such that the intra-storage charges when the inter-storage discharges, meaning
that $p_{\tau}^{\mathrm{C,a*}}>0$ and
$\sum_{v\in\mathcal{V}}p_{v\tau}^{\mathrm{D,e*}}>0$.

Let’s call $q_{v}^{\tau}$ the quantity discharged from the inter storage with
value $v$ at $\tau$, i.e. $q_{v}^{\tau}=p_{v\tau}^{\mathrm{D,e*}}$. We also
call $q^{\tau}$ the total quantity discharged from inter storage at $\tau$,
$q^{\tau}=\sum_{v\in\mathcal{V}}q_{v}^{\tau}$.

We now identify another time period $\kappa\in\mathcal{T}$ such that
$p_{\kappa}^{\mathrm{C,a*}}<0$.

#### Existence of $\kappa$

We first prove that $\kappa$ exists. Let’s suppose that it does not, i.e. we
suppose that $p_{t}^{\mathrm{C,a*}}\geq 0$, $\forall t\in\mathcal{T}$. We have
$\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}>0$ since $p_{\tau}^{\mathrm{C,a*}}>0$.

We build a new solution that we identify with the superscript ′. It is
identical to the previous solution, except for
$p_{\tau}^{\mathrm{C,a^{\prime}}}=p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}$
and
$p_{v\tau}^{\mathrm{D,e^{\prime}}}=p_{v\tau}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}}$,
$\forall v\in\mathcal{V}$, with
$q^{\tau^{\prime}}=\min\\{q^{\tau},\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}\\}$ and
$q_{v}^{\tau^{\prime}}$ are such that
$q^{\tau^{\prime}}=\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$ and $0\leq
q_{v}^{\tau^{\prime}}\leq q_{v}^{\tau}$, $\forall v\in\mathcal{V}$. It is
possible to find such $q_{v}^{\tau^{\prime}}$, since
$\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}>0$ and $q^{\tau}>0$, so
$0<q^{\tau^{\prime}}\leq q^{\tau}$, meaning that
$0<\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}\leq\sum_{v\in\mathcal{V}}q_{v}^{\tau}$.
We check that this new solution is feasible.

For (6b), at $t=\tau$, we now have

$\displaystyle\sum_{l\in\mathcal{L}}d_{l\tau}^{\prime}+p_{\tau}^{\mathrm{C,a^{\prime}}}-\sum_{g\in\mathcal{G}}p_{g\tau}^{\prime}-\sum_{v\in\mathcal{V}}p_{v\tau}^{\mathrm{D,e^{\prime}}}=$
(7a) $\displaystyle=$
$\displaystyle\sum_{l\in\mathcal{L}}d_{l\tau}^{*}+p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}-\sum_{g\in\mathcal{G}}p_{g\tau}^{*}-\sum_{v\in\mathcal{V}}(p_{v{\tau}}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}})$
(7b) $\displaystyle=$
$\displaystyle\sum_{l\in\mathcal{L}}d_{l\tau}^{*}+p_{\tau}^{\mathrm{C,a*}}-\sum_{g\in\mathcal{G}}p_{g\tau}^{*}-\sum_{v\in\mathcal{V}}p_{v\tau}^{\mathrm{D,e*}}-q^{\tau^{\prime}}+\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$
(7c) $\displaystyle=$
$\displaystyle\sum_{l\in\mathcal{L}}d_{l\tau}^{*}+p_{\tau}^{\mathrm{C,a*}}-\sum_{g\in\mathcal{G}}p_{g\tau}^{*}-\sum_{v\in\mathcal{V}}p_{v\tau}^{\mathrm{D,e*}},$
(7d)

so the constraint (6b) is satisfied.

Constraints (6c) and (6d) still hold, since the solution is not modified for
these variables.

For constraint (6e), we have

$\displaystyle\sum_{i=1}^{T}p_{i}^{\mathrm{C,a^{\prime}}}$
$\displaystyle=\sum_{i=1}^{\tau-1}p_{i}^{\mathrm{C,a^{\prime}}}+p_{\tau}^{\mathrm{C,a^{\prime}}}+\sum_{i=\tau+1}^{T}p_{i}^{\mathrm{C,a^{\prime}}}$
(8a)
$\displaystyle=\sum_{i=1}^{\tau-1}p_{i}^{\mathrm{C,a*}}+p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}+\sum_{i=\tau+1}^{T}p_{i}^{\mathrm{C,a*}}$
(8b) $\displaystyle=\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}-q^{\tau^{\prime}},$
(8c)

which we know is positive, since
$q^{\tau^{\prime}}=\min\\{q^{\tau},\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}\\}$.

For constraint (6f), there are no changes for $t<\tau$. For $t\geq\tau$, we
have

$\displaystyle\sum_{i=1}^{t}p_{i}^{\mathrm{C,a^{\prime}}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e^{\prime}}}\right)=$
(9a) $\displaystyle=$
$\displaystyle\sum_{i=1}^{\tau-1}p_{i}^{\mathrm{C,a^{\prime}}}+p_{\tau}^{\mathrm{C,a^{\prime}}}+\sum_{i=\tau+1}^{t}p_{i}^{\mathrm{C,a^{\prime}}}$
$\displaystyle+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{\tau-1}p_{vi}^{\mathrm{D,e^{\prime}}}-p_{v\tau}^{\mathrm{D,e^{\prime}}}-\sum_{i=\tau+1}^{t}p_{vi}^{\mathrm{D,e^{\prime}}}\right)$
(9b) $\displaystyle=$
$\displaystyle\sum_{i=1}^{\tau-1}p_{i}^{\mathrm{C,a*}}+p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}+\sum_{i=\tau+1}^{t}p_{i}^{\mathrm{C,a*}}$
$\displaystyle+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{\tau-1}p_{vi}^{\mathrm{D,e*}}-p_{v\tau}^{\mathrm{D,e*}}+q_{v}^{\tau^{\prime}}-\sum_{i=\tau+1}^{t}p_{vi}^{\mathrm{D,e*}}\right)$
(9c) $\displaystyle=$
$\displaystyle\sum_{i=1}^{t}p_{i}^{\mathrm{C,a*}}-q^{\tau^{\prime}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e*}}\right)+\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$
(9d) $\displaystyle=$
$\displaystyle\sum_{i=1}^{t}p_{i}^{\mathrm{C,a*}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e*}}\right),$
(9e)

so the constraint still stands.

For (6g), we only have to check at $t=\tau$. We defined
$q_{v}^{\tau^{\prime}}\leq q_{v}^{\tau}$, to make sure that this constraint is
satisfied.

Constraint (6h) still holds, since we discharge less. We have

$\displaystyle
E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e^{\prime}}}$
$\displaystyle=E^{\mathrm{init}}_{v}-\sum_{i=1}^{\tau-1}p_{vi}^{\mathrm{D,e^{\prime}}}-p_{v\tau}^{\mathrm{D,e^{\prime}}}-\sum_{i=\tau+1}^{T}p_{vi}^{\mathrm{D,e^{\prime}}}$
(10a)
$\displaystyle=E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e*}}+q_{v}^{\tau^{\prime}})$
(10b) $\displaystyle\geq
E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e*}}.$ (10c)

For (6i), using previous calculation, we have that

$\small\sum_{i=1}^{T}p_{i}^{\mathrm{C,a^{\prime}}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e^{\prime}}}\right)=\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e*}}\right),$
(11)

so this constraint is also satisfied. We can conclude that our new solution is
feasible.

The value of the objective function for this solution is

$\displaystyle\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}^{{}^{\prime}}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}^{{}^{\prime}}-\sum_{v\in\mathcal{V}}S_{v}p_{vt}^{\mathrm{D,e^{\prime}}}\right)$
(12a) $\displaystyle=$
$\displaystyle\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}^{*}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}^{*}\right)-\sum_{v\in\mathcal{V}}S_{v}(\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}})$
(12b) $\displaystyle=$
$\displaystyle\sum_{t\in\mathcal{T}}\left(\sum_{l\in\mathcal{L}}U_{lt}d_{lt}^{*}-\sum_{g\in\mathcal{G}}C_{gt}p_{gt}^{*}-\sum_{v\in\mathcal{V}}S_{v}p_{vi}^{\mathrm{D,e*}}\right)+\sum_{v\in\mathcal{V}}S_{v}q_{v}^{\tau^{\prime}}.$
(12c)

We have $\sum_{v\in\mathcal{V}}S_{v}q_{v}^{\tau^{\prime}}>0$, under the
assumption that $S_{v}>0$, $\forall v\in\mathcal{V}$ and since
$q_{v}^{\tau^{\prime}}\geq 0$ , $\forall v\in\mathcal{V}$ and
$\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}>0$, meaning that there is at
least one $v$ for which $q_{v}^{\tau^{\prime}}>0$.

It is greater than the value of the objective function for the ∗ solution,
which is optimal, indicating a contradiction.

Note 1: This also excludes the case where $|\mathcal{T}|=1$, meaning that it
cannot be optimal to have $p_{\tau}^{\mathrm{C,a*}}>0$ and
$p_{\tau}^{\mathrm{D,e*}}>0$ in this case.

Note 2: $\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}>0$ also in case of net charge, so
there will not be simultaneous charge of intra-storage and discharge of inter-
storage then.

#### Building a new solution

We now know that there exists another time period $\kappa\in\mathcal{T}$ and
$\kappa\neq\tau$ such that $p_{\kappa}^{\mathrm{C,a*}}<0$. We build a new
solution to our problem, which we identify with the superscript ′. It is
identical to the previous solution, except for $t=\tau$ and $t=\kappa$. We
want to discharge less the inter-storage at $t=\tau$ and discharge it more at
$t=\kappa$, and in turn charge less the intra-storage at $t=\tau$ and
discharge it less at $t=\kappa$, in order to keep the same storage level when
both are summed.

For $t=\tau$,
$p_{\tau}^{\mathrm{C,a^{\prime}}}=p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}$
and
$p_{v\tau}^{\mathrm{D,e^{\prime}}}=p_{v\tau}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}}$,
$\forall v\in\mathcal{V}$, with
$q^{\tau^{\prime}}=\min\\{q^{\tau},-p_{\kappa}^{\mathrm{C,a*}}\\}$ and
$q_{v}^{\tau^{\prime}}$ are such that
$q^{\tau^{\prime}}=\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$ and $0\leq
q_{v}^{\tau^{\prime}}\leq q_{v}^{\tau}$, $\forall v\in\mathcal{V}$. It is
possible to find such $q_{v}^{\tau^{\prime}}$, since
$-p_{\kappa}^{\mathrm{C,a*}}>0$ and $q^{\tau}>0$, so $0<q^{\tau^{\prime}}\leq
q^{\tau}$, meaning that
$0<\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}\leq\sum_{v\in\mathcal{V}}q_{v}^{\tau}$.
For $t=\kappa$,
$p_{\kappa}^{\mathrm{C,a^{\prime}}}=p_{\kappa}^{\mathrm{C,a*}}+q^{\tau^{\prime}}$
and
$p_{v\kappa}^{\mathrm{D,e^{\prime}}}=p_{v\kappa}^{\mathrm{D,e*}}+q_{v}^{\tau^{\prime}}$,
$\forall v\in\mathcal{V}$.

Let’s check that this solution is feasible.

For (6b), at $t=\tau$, we have the same as in (7), which is feasible.

At $t=\kappa$,

$\displaystyle\sum_{l\in\mathcal{L}}d_{l\kappa}^{\prime}+p_{\kappa}^{\mathrm{C,a^{\prime}}}-\sum_{g\in\mathcal{G}}p_{g\kappa}^{\prime}-\sum_{v\in\mathcal{V}}p_{v\kappa}^{\mathrm{D,e^{\prime}}}$
(13a) $\displaystyle=$
$\displaystyle\sum_{l\in\mathcal{L}}d_{l\kappa}^{*}+p_{\kappa}^{\mathrm{C,a*}}+q^{\tau^{\prime}}-\sum_{g\in\mathcal{G}}p_{g\kappa}^{*}-\sum_{v\in\mathcal{V}}(p_{v{\kappa}}^{\mathrm{D,e*}}+q_{v}^{\tau^{\prime}})$
(13b) $\displaystyle=$
$\displaystyle\sum_{l\in\mathcal{L}}d_{l\kappa}^{*}+p_{\kappa}^{\mathrm{C,a*}}-\sum_{g\in\mathcal{G}}p_{g\kappa}^{*}-\sum_{v\in\mathcal{V}}p_{v\kappa}^{\mathrm{D,e*}}+q^{\tau^{\prime}}-\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$
(13c) $\displaystyle=$
$\displaystyle\sum_{l\in\mathcal{L}}d_{l\kappa}^{*}+p_{\kappa}^{\mathrm{C,a*}}-\sum_{g\in\mathcal{G}}p_{g\kappa}^{*}-\sum_{v\in\mathcal{V}}p_{v\kappa}^{\mathrm{D,e*}},$
(13d)

so the constraint (6b) is satisfied.

Constraints (6c) and (6d) still hold, since the solution is not modified for
these variables.

Let’s call $t_{1}=\min\\{\tau,\kappa\\}$ and $t_{2}=\max\\{\tau,\kappa\\}$.
The total charged quantity in the intra-storage is unchanged:

$\displaystyle\sum_{i=1}^{T}p_{i}^{\mathrm{C,a^{\prime}}}$
$\displaystyle=\sum_{i=1}^{t_{1}-1}p_{i}^{\mathrm{C,a^{\prime}}}+\sum_{i=t_{1}+1}^{t_{2}-1}p_{i}^{\mathrm{C,a^{\prime}}}+\sum_{i=t_{2}+1}^{T}p_{i}^{\mathrm{C,a^{\prime}}}$
$\displaystyle+p_{\tau}^{\mathrm{C,a^{\prime}}}+p_{\kappa}^{\mathrm{C,a^{\prime}}}$
(14a)
$\displaystyle=\sum_{i=1}^{t_{1}-1}p_{i}^{\mathrm{C,a*}}+\sum_{i=t_{1}+1}^{t_{2}-1}p_{i}^{\mathrm{C,a*}}+\sum_{i=t_{2}+1}^{T}p_{i}^{\mathrm{C,a*}}$
$\displaystyle+p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}+p_{\kappa}^{\mathrm{C,a*}}+q^{\tau^{\prime}}$
(14b) $\displaystyle=\sum_{i=1}^{T}p_{i}^{\mathrm{C,a*}}.$ (14c)

It is also the case for the inter-storage:

$\displaystyle\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e^{\prime}}}$
$\displaystyle=\sum_{i=1}^{t_{1}-1}p_{vi}^{\mathrm{D,e^{\prime}}}+\sum_{i=t_{1}+1}^{t_{2}-1}p_{vi}^{\mathrm{D,e^{\prime}}}+\sum_{i=t_{2}+1}^{T}p_{vi}^{\mathrm{D,e^{\prime}}}$
$\displaystyle+p_{v\tau}^{\mathrm{D,e^{\prime}}}+p_{v\kappa}^{\mathrm{D,e^{\prime}}}$
(15a) $\displaystyle=$
$\displaystyle\sum_{i=1}^{t_{1}-1}p_{vi}^{\mathrm{D,e*}}+\sum_{i=t_{1}+1}^{t_{2}-1}p_{vi}^{\mathrm{D,e*}}+\sum_{i=t_{2}+1}^{T}p_{vi}^{\mathrm{D,e*}}$
$\displaystyle+p_{v\tau}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}}+p_{v\kappa}^{\mathrm{D,e*}}+q_{v}^{\tau^{\prime}}$
(15b) $\displaystyle=$ $\displaystyle\sum_{i=1}^{T}p_{vi}^{\mathrm{D,e*}}.$
(15c)

As a consequence, constraints (6e), (6h) and (6i) are satisfied by the new
solution.

For constraint (6f), there are no changes for $t<t_{1}$. For $t_{1}\leq
t<t_{2}$. If $t_{1}=\tau$, we have the same as in (9). If $t_{1}=\kappa$,

$\displaystyle\sum_{i=1}^{t}p_{i}^{\mathrm{C,a^{\prime}}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e^{\prime}}}\right)$
(16a) $\displaystyle=$
$\displaystyle\sum_{i=1}^{\kappa-1}p_{i}^{\mathrm{C,a^{\prime}}}+p_{\kappa}^{\mathrm{C,a^{\prime}}}+\sum_{i=\kappa+1}^{t}p_{i}^{\mathrm{C,a^{\prime}}}$
$\displaystyle+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{\kappa-1}p_{vi}^{\mathrm{D,e^{\prime}}}-p_{v\kappa}^{\mathrm{D,e^{\prime}}}-\sum_{i=\kappa+1}^{t}p_{vi}^{\mathrm{D,e^{\prime}}}\right)$
(16b) $\displaystyle=$
$\displaystyle\sum_{i=1}^{\kappa-1}p_{i}^{\mathrm{C,a*}}+p_{\kappa}^{\mathrm{C,a*}}+q^{\tau^{\prime}}+\sum_{i=\kappa+1}^{t}p_{i}^{\mathrm{C,a*}}$
$\displaystyle+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{\kappa-1}p_{vi}^{\mathrm{D,e*}}-p_{v\kappa}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}}-\sum_{i=\kappa+1}^{t}p_{vi}^{\mathrm{D,e*}}\right)$
(16c)
$\displaystyle=\sum_{i=1}^{t}p_{i}^{\mathrm{C,a*}}+q^{\tau^{\prime}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e*}}\right)-\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$
(16d)
$\displaystyle=\sum_{i=1}^{t}p_{i}^{\mathrm{C,a*}}+\sum_{v\in\mathcal{V}}\left(E^{\mathrm{init}}_{v}-\sum_{i=1}^{t}p_{vi}^{\mathrm{D,e*}}\right).$
(16e)

Then, it is easy to see that the constraint also stands for $t\geq t_{2}$.

For (6g), at $t=\tau$, the constraint is satisfied since
$q_{v}^{\tau^{\prime}}\leq q_{v}^{\tau}$. At $t=\kappa$ it is also satisfied
since $q_{v}^{\tau^{\prime}}\geq 0$.

All constraints are satisfied, so this new solution is feasible. The value of
the objective function for this solution is the same since $d_{lt}$ and
$p_{gt}$ are the same and the total energy discharged for the inter-storage is
the same.

We have thus found another solution in which the quantity discharged in $\tau$
is smaller. We can repeat the same procedure until there is no quantity
discharged at $\tau$.

#### Recursion

We can repeat the following until there is no simultaneous charge of the
intra-storage and discharge of the inter-storage:

  1. 1.

We consider $\tau$, for which $p_{\tau}^{\mathrm{C,a*}}>0$ and
$\sum_{v\in\mathcal{V}}p_{v\tau}^{\mathrm{D,e*}}>0$. We have
$q_{v}^{\tau}=p_{v\tau}^{\mathrm{D,e*}}$, $\forall v\in\mathcal{V}$, and
$q^{\tau}=\sum_{v\in\mathcal{V}}q_{v}^{\tau}$.

  2. 2.

We identify $\kappa$ such that $p_{\kappa}^{\mathrm{C,a*}}<0$.

  3. 3.

We modify the solution. We identify
$q^{\tau^{\prime}}=\min\\{q^{\tau},-p_{\kappa}^{\mathrm{C,a*}}\\}$ and
$q_{v}^{\tau^{\prime}}$, such that
$q^{\tau^{\prime}}=\sum_{v\in\mathcal{V}}q_{v}^{\tau^{\prime}}$ and $0\leq
q_{v}^{\tau^{\prime}}\leq q_{v}^{\tau}$, $\forall v\in\mathcal{V}$. For
$t=\tau$,
$p_{\tau}^{\mathrm{C,a^{\prime}}}=p_{\tau}^{\mathrm{C,a*}}-q^{\tau^{\prime}}$
and
$p_{v\tau}^{\mathrm{D,e^{\prime}}}=p_{v\tau}^{\mathrm{D,e*}}-q_{v}^{\tau^{\prime}}$.
For $t=\kappa$,
$p_{\kappa}^{\mathrm{C,a^{\prime}}}=p_{\kappa}^{\mathrm{C,a*}}+q^{\tau^{\prime}}$
and
$p_{v\kappa}^{\mathrm{D,e^{\prime}}}=p_{v\kappa}^{\mathrm{D,e*}}+q_{v}^{\tau^{\prime}}$,
$\forall v\in\mathcal{V}$.

  4. 4.

We update $q^{\mathrm{\tau,new}}=q^{\mathrm{\tau,old}}-q^{\tau^{\prime}}$.

  5. 5.

If $q^{\mathrm{\tau,new}}>0$, we go back to step 2 and identify a new
$\kappa$. If $q^{\mathrm{\tau,new}}=0$, we go to step 6.

  6. 6.

Check if there is another $t$ for which $p_{t}^{\mathrm{C,a*}}>0$ and
$\sum_{v\in\mathcal{V}}p_{vt}^{\mathrm{D,e*}}>0$. If so, it is the new $\tau$
and we go back to 1.

We have thus shown that there is always an optimal solution in which there is
no simultaneous charge and discharge, under the assumption that $S_{v}>0$,
$\forall v\in\mathcal{V}$.

∎

The assumption that $S_{v}>0$, $\forall v\in\mathcal{V}$ is reasonable since
there is no interest to charge the storage system when prices are non-
positive. Moreover, charging at a value of 0 is equivalent to not valuing the
stored energy, in which case we can equivalently use the original formulation,
in which there is no separation between intra- and inter-storage and thus no
simultaneous charge and discharge of intra- and inter-storage.

###  III-D Update of the Inter-Storage

In the case of net discharge, the inter-storage can be easily updated by
subtracting the sum of $p_{vt}^{\mathrm{D,e}}$ to $E^{\mathrm{init}}_{v}$. If
doing so, the inter-storage for value $v$ is empty, the corresponding index is
dropped in $\mathcal{V}$.

In the case of net charge, the quantity net charged is added to the inter-
storage, with the corresponding charging price. However, since the storage
might have been charging at different prices during the market interval, the
question of which price to save arises. To address this, we introduce the
following optimization model

$\displaystyle\min_{\mathbf{p^{\mathrm{C,loc}}},\mathbf{p^{\mathrm{C,mov}}}}\quad$
$\displaystyle-\Delta
t\sum_{t\in\mathcal{T}}\lambda^{*}_{t}p^{\mathrm{C,loc}}$ (17a) s.t.
$\displaystyle-\sum_{t\in\mathcal{T}}\lambda^{*}_{t}p^{\mathrm{C,loc}}\geq 0$
(17b) $\displaystyle\sum_{t\in\mathcal{T}}p_{t}^{\mathrm{C,loc}}=0$ (17c)
$\displaystyle
p_{t}^{\mathrm{C,loc}}+p_{t}^{\mathrm{C,mov}}=p_{t}^{\mathrm{C,a*}},\,\hskip
69.38078pt\forall t\in\mathcal{T}$ (17d) $\displaystyle
p_{t}^{\mathrm{C,loc}}=p_{t}^{\mathrm{C,a*}},\,\hskip 60.70653pt\forall
t\in\mathcal{T},\,p_{t}^{\mathrm{C,a*}}\leq 0$ (17e) $\displaystyle
p_{t}^{\mathrm{C,loc}}\geq 0,\,\hskip 91.0631pt\forall
t\in\mathcal{T},\,p_{t}^{\mathrm{C,a*}}>0$ (17f) $\displaystyle
p_{t}^{\mathrm{C,mov}}\geq 0,\,\hskip 147.4292pt\forall t\in\mathcal{T}.$
(17g)

The parameters in this model are obtained from the solution of the market
clearing, and indicated with ∗. The idea is to split the quantity charged in
the intra-storage into a local and a moved quantity, $p_{t}^{\mathrm{C,loc}}$
and $p_{t}^{\mathrm{C,mov}}$. Then, $p_{t}^{\mathrm{C,mov}}$ is added to the
inter-storage, with the market price at $t$. The total moved quantity must
correspond to the net charge, which is equivalent to setting the total local
quantity to zero (sum of total local charge and discharge), in (17c). The
inter-storage is only updated by quantities charged, which is ensured by
(17g): the quantity moved can only be positive. The sum of local and moved
quantities has to correspond to the quantity charged in the inter-storage,
with (17d). In case of discharge, constraint (17e) applies, ensuring that the
moved quantity is equal to zero. In case of charge, constraint (17f) prevents
discharge for the local quantity.

The objective (17a) is to minimize the storage profit over this interval,
while still ensuring that it is positive with (17b). We are minimizing in
order to avoid situations where the storage system would be setting prices
higher than necessary.

##  IV Study of Market Properties

We study cost recovery and social welfare for the market clearing with VLB,
comparing to the split market clearing, (1) with (3). We discuss the impact of
imperfect foresight on these properties. Finally, we discuss a limitation of
both models, opening directions for future research.

###  IV-A Cost Recovery

There is cost recovery for a market participant, if their surplus is always
non-negative. We can easily show that this is the case for the generators and
loads444We prove cost recovery for loads and generators by formulating their
individual surplus maximization problem and its dual problem and using the
strong duality theorem..

For a storage system, it is not relevant to ensure that the surplus is always
non-negative. For example, the storage system could be only paying to charge
in one market interval, to later discharge at a higher price and make a profit
in a subsequent market interval. In this example, the surplus of the storage
in the first market interval would be negative. Hence, we redefine cost
recovery for a non-merchant storage system.

######  Definition 1 (Cycle and cost recovery for a non-merchant storage).

We define a cycle as a group of consecutive market intervals for which the
storage system is initially empty and finally returns to this same state. We
say that there is cost recovery for a non-merchant storage system if its
surplus over a cycle is non-negative.

For a non-merchant storage system, and for a fair comparison, it only makes
sense to look at the surplus over a cycle. Otherwise, we need to know the
value of the stored energy to include it, which is actually the complex
problem that we are trying to solve here.

We argue that cost recovery stands by design of the market clearing and of the
update of inter-storage from (17), under the assumption that looking far
enough into the future, the inter-storage will eventually be empty. In this
model, constraint (17b) ensures that the profit of the storage system is non-
negative for the quantity exchanged over the time interval. Regarding the
quantity moved to the inter-storage, the charging price is saved and later
used as a bid. Since discharge is never imposed, this ensures that for this
quantity the price received will be at least equal to the charging price. Note
that this result is not based on the assumption of perfect foresight and is
thus valid in uncertain settings. The assumption that the storage will
eventually be empty is mild. It would be challenged in the case that a
decision was made to charge the storage in a day with very high prices that
never occur again. A complete day with very high prices is very unlikely to
happen without being foreseen. And if foreseen, this situation would
illustrate very poor decision-making on the final level of the storage system.

We saw in the example of Section II-C that cost recovery is not ensured for
the storage system in the split market clearing. In that example, the storage
system starts empty and finishes empty, so the two market intervals considered
are a cycle. However, the surplus of the storage in this cycle can be as low
as -3€. We run that same example with the market clearing with VLB, to show
that this situation does not arise anymore. The results are shown in Table
III. At the end of the first MI, the inter-storage is charged with a price of
5€/MWh and it is discharged in the next hour. The storage is now marginal and
gets paid at least 5€/MWh, so that the minimum surplus of the storage system
over this cycle is 0€, and cost recovery stands.

TABLE III: Results for the example of Section II-C, for the market clearing with VLB MI | $t$ | $d_{1t}$ | $p_{1t}$ | $p_{2t}$ | $p_{t}^{\mathrm{C}}$ | $p_{t}^{\mathrm{D,e}}$ | $e^{\mathrm{a}}_{t}$ | $e^{\mathrm{e}}_{t}$ | $\lambda_{t}$
---|---|---|---|---|---|---|---|---|---
1 | 1 | 0 | 1 | 0 | 1 | 0 | 1 | 0 | 5
1 | 2 | 3 | 2 | 0 | 0 | 1 | 0 | 0 | [5,9]

###  IV-B Social Welfare

If the level of the storage system at the end of each interval of (5) is set
to the value that is optimal for the ideal market clearing, the storage system
will follow the same trajectory, which will result in the same social welfare.
Indeed, since this level is obtained in a way that the storage system recovers
its costs in the ideal clearing, the inter-storage will discharge the same
quantity.

We now discuss different scenarios that can occur in case of imperfect
information and error in setting the level of the storage system at the end of
each market interval, and illustrate them with simple examples. We look at the
social welfare for the market clearing with VLB, in comparison to the one
obtained for the split market clearing using the same values for the final
levels $E^{\mathrm{end}}$. Since social welfare is the sum of the surpluses of
all the market participants, its calculation has to be carried out over a
cycle, as defined in Definition 1. We compare both methods on a cycle for the
market-clearing with VLB 555A cycle for the split market-clearing is not
necessarily a cycle for the market-clearing with VLB but the opposite is true.
If $E^{\mathrm{end}}=0$, the final level with VLB might be higher because of
(5j), while if the final level with VLB is equal to zero, it means that
$E^{\mathrm{end}}=0$ and the final level for the split market is also equal to
zero.. We also compare to the outcome of the ideal market clearing. In both
illustrative examples, the storage has the same capacity $\overline{E}=2.5$
MWh and it is initially empty, each MI consists of one time period of one
hour, and there is one load and one generator participating in the market. The
rest of the data and results are given in Tables IV and V, including the
number of MIs. The results in column “$e_{t}$” for the split market correspond
to the values set for $E^{\mathrm{end}}$, which are also used for VLB.

The first scenario corresponds to the case where discharge is imposed in a
period of low prices, also limiting the availability of storage for future
periods with low generation. In this case, the social welfare can be higher
for (5) than for (1) with (3). We show this in the example of Table IV. The
total social welfare for these three MIs is -1€ for the split market clearing
and 16€ for the market clearing with VLB, compared to 21€ for the ideal
clearing. Due to the error in predicting the ideal final levels for the
storage, the social welfare is lower than in the ideal clearing. The
difference in social welfare between split and VLB is due to the fact that in
the split market clearing, the storage discharges on the second MI, while it
does not in the market clearing with VLB. Indeed, the prices are too low for
the storage to recover the costs from charging on the first MI. Rather, it
discharges on the third MI, thereby allowing for the load to be completely
supplied, which increases social welfare. We also observe that the formulation
with VLB reestablishes cost recovery for the storage system.

TABLE IV: Data and results for the first example of Section IV-B. Here
$p^{\mathrm{D,e}}_{t}=\sum_{v\in\mathcal{V}}p^{\mathrm{D,e}}_{vt}$ and
$e^{\mathrm{e}}_{t}=\sum_{v\in\mathcal{V}}e^{\mathrm{e}}_{vt}$. The state of
energy is in MWh, prices in €/MWh, and the rest is in MW.

|  | Data | Results ideal | Results split | Results VLB
---|---|---|---|---|---
MI | $t$ | $D_{t}$ | $U_{t}$ | $P_{t}$ | $C_{t}$ | $d_{t}$ | $p_{t}$ | $p^{\mathrm{C}}_{t}$ | $e_{t}$ | $\lambda_{t}$ | $d_{t}$ | $p_{t}$ | $p^{\mathrm{C}}_{t}$ | $e_{t}$ | $\lambda_{t}$ | $d_{t}$ | $p_{t}$ | $p^{\mathrm{C,a}}_{t}$ | $p^{\mathrm{D,e}}_{t}$ | $e^{\mathrm{a}}_{t}$ | $e^{\mathrm{e}}_{t}$ | $\lambda_{t}$
1 | 1 | 1 | 5 | 4 | 5 | 0 | 0 | 0 | 0 | 5 | 0 | 2.5 | 2.5 | 2.5 | 5 | 0 | 2.5 | 2.5 | 0 | 2.5 | 0 | 5
2 | 1 | 3 | 4 | 6 | 3 | 3 | 5.5 | 2.5 | 2.5 | 3 | 3 | 0.5 | -2.5 | 0 | 3 | 3 | 3 | 0 | 0 | 0 | 2.5 | 3
3 | 1 | 3 | 10 | 1 | 9 | 3 | 0.5 | -2.5 | 0 | 9 | 1 | 1 | 0 | 0 | 9 | 3 | 0.5 | 0 | 2.5 | 0 | 0 | 9

The second scenario corresponds to the case where due to an error when
choosing the final level, charge is imposed in a period of high prices, which
do not occur again soon enough, preventing the complete use of the storage
capacity in the following market intervals. Then, the social welfare can be
lower for (5) than for (1) with (3). We show this in the example of Table V.
The total social welfare for these three MIs is 842.5€ for the split market
clearing, and 772.5€ for the market clearing with VLB, compared to 855€ for
the ideal clearing. The difference in social welfare between split and VLB is
due to the fact that the inter-storage does not discharge on the second and
fourth MIs. Indeed, the prices are too low for the storage to recover the
costs from charging on the first MI. More expensive generators are used to
supply the load on these MIs, thereby decreasing social welfare.

TABLE V: Data and results for the second example of Section IV-B. Here
$p^{\mathrm{D,e}}_{t}=\sum_{v\in\mathcal{V}}p^{\mathrm{D,e}}_{vt}$ and
$e^{\mathrm{e}}_{t}=\sum_{v\in\mathcal{V}}e^{\mathrm{e}}_{vt}$. The state of
energy is in MWh, prices in €/MWh, and the rest is in MW. Load values are
$d_{t}=D_{t}$ and market prices are $\lambda_{t}=C_{t}$, the same in all the
cases.

|  | Data | Results ideal | Results split | Results VLB | Results VLB - Discount
---|---|---|---|---|---|---
MI | $t$ | $D_{t}$ | $U_{t}$ | $P_{t}$ | $C_{t}$ | $p_{t}$ | $p^{\mathrm{C}}_{t}$ | $e_{t}$ | $p_{t}$ | $p^{\mathrm{C}}_{t}$ | $e_{t}$ | $p_{t}$ | $p^{\mathrm{C,a}}_{t}$ | $p^{\mathrm{D,e}}_{t}$ | $e^{\mathrm{a}}_{t}$ | $e^{\mathrm{e}}_{t}$ | $p_{t}$ | $p^{\mathrm{C,a}}_{t}$ | $p^{\mathrm{D,e}}_{t}$ | $e^{\mathrm{a}}_{t}$ | $e^{\mathrm{e}}_{t}$
1 | 1 | 10 | 25 | 15 | 20 | 10 | 0 | 0 | 12.5 | 2.5 | 2.5 | 20 | 2.5 | 0 | 2.5 | 0 | 12.5 | 2.5 | 0 | 2.5 | 0
2 | 1 | 10 | 25 | 15 | 15 | 10 | 0 | 0 | 7.5 | -2.5 | 0 | 10 | 0 | 0 | 0 | 2.5 | 10 | 0 | 0 | 0 | 2.5
3 | 1 | 10 | 25 | 15 | 1 | 12.5 | 2.5 | 2.5 | 12.5 | 2.5 | 2.5 | 10 | 0 | 0 | 0 | 2.5 | 10 | 0 | 0 | 0 | 2.5
4 | 1 | 10 | 25 | 15 | 15 | 7.5 | -2.5 | 0 | 7.5 | -2.5 | 0 | 10 | 0 | 0 | 0 | 2.5 | 10 | 0 | 2.5 | 0 | 0
5 | 1 | 10 | 25 | 15 | 1 | 12.5 | 2.5 | 2.5 | 12.5 | 2.5 | 2.5 | 10 | 0 | 0 | 0 | 2.5 | 10 | 2.5 | 0 | 2.5 | 0
6 | 1 | 10 | 25 | 15 | 21 | 7.5 | -2.5 | 0 | 7.5 | -2.5 | 0 | 7.5 | 0 | 2.5 | 0 | 0 | 7.5 | 0 | 2.5 | 0 | 0

However, this can be limited by introducing a discount on stored value, which
is shown in the last column of Table V. In this case, we decreased by 25% the
value of the energy stored in the inter-storage after each MI, except after
the MI in which it is stored. The social welfare in this case is 807.5€, as
the inter-storage can now be used on the fourth MI.

###  IV-C Limitation of Split Models

We have seen that our new method performs better than the traditional split
market clearing in terms of cost recovery for the market participants.
However, not all the limitations of the split market clearing are overcome,
which we see here in a last illustrative example.

We use again a storage system with capacity $\overline{E}=2.5$ MWh and
initially empty. The rest of the data and the results are shown in Table VI.
There is one load and one generator. The market is cleared for two MIs of
three hours each and the final storage level is set to its optimal value,
obtained when clearing for the two MIs together. Note that the final level at
the end of the second MI is equal to 2.5 MWh because of subsequent MIs, which
are not represented here.

We are in the perfect foresight set-up and we can see that all quantities
agree. As a consequence, social welfare is the same for all these types of
market clearing. However, we see that the prices in the last hour of the first
MI can potentially be higher than in the ideal clearing, for both other
approaches. This difference is due to future market intervals not being
considered for pricing in these cases. In the ideal clearing, the price for
that hour is at most 4€/MWh, which corresponds to the utility of the load in
the first hour of the next MI, which the other approaches do not take into
consideration. We argue for a limited impact of future uncertain information
on the formation of current prices, to ensure transparency. However, we see
here that it comes with a cost that somebody will ultimately have to pay. The
study of this trade-off is a topic for future research.

TABLE VI: Data and results for the example of Section IV-C. Here $p^{\mathrm{D,e}}_{t}=\sum_{v\in\mathcal{V}}p^{\mathrm{D,e}}_{vt}$ and $e^{\mathrm{e}}_{t}=\sum_{v\in\mathcal{V}}e^{\mathrm{e}}_{vt}$. The state of energy is in MWh, prices in €/MWh, and the rest is in MW. Load values $d_{t}$ and generation values $p_{t}$ are the same in all the cases and are not repeated |  | Data | Results ideal | Results split | Results VLB
---|---|---|---|---|---
MI | $t$ | $D_{t}$ | $U_{t}$ | $P_{t}$ | $C_{t}$ | $d_{t}$ | $p_{t}$ | $p^{\mathrm{C}}_{t}$ | $e_{t}$ | $\lambda_{t}$ | $p^{\mathrm{C}}_{t}$ | $e_{t}$ | $\lambda_{t}$ | $p^{\mathrm{C,a}}_{t}$ | $p^{\mathrm{D,e}}_{t}$ | $e^{\mathrm{a}}_{t}$ | $e^{\mathrm{e}}_{t}$ | $\lambda_{t}$
1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 1 | 2 | 1 | 1 | 2 | 1 | 0 | 1 | 0 | 2
1 | 2 | 1 | 5 | 2.5 | 2 | 1 | 2.5 | 1.5 | 2.5 | 2 | 1.5 | 2.5 | 2 | 1.5 | 0 | 2.5 | 0 | 2
1 | 3 | 2 | 6 | 0 | 0 | 2 | 0 | -2 | 0.5 | [3,4] | -2 | 0.5 | [3,6] | -2 | 0 | 0.5 | 0 | [3,6]
2 | 1 | 2.5 | 4 | 2 | 3 | 2.5 | 2 | -0.5 | 0 | [3,4] | -0.5 | 0 | [3,4] | -0.5 | 0 | -0.5 | 0.5 | [3,4]
2 | 2 | 0 | 0 | 1 | 2 | 0 | 1 | 1 | 1 | [3,7] | 1 | 1 | [3,7] | 1 | 0 | 0.5 | 0.5 | [3,7]
2 | 3 | 4 | 7 | 5.5 | 3 | 4 | 5.5 | 1.5 | 2.5 | [3,7] | 1.5 | 2.5 | [3,7] | 1.5 | 0 | 2 | 0.5 | [3,7]

##  V Discussion and Conclusion

We introduced a novel procedure for clearing an energy market with non-
merchant storage, using virtual linking bids. This is based on an artificial
representation of the storage system, dividing it into a component for local
arbitrage, within this market interval, and a component for arbitrage between
market intervals. We showed that it outperforms traditional approaches when it
comes to cost recovery. Indeed, it ensures cost recovery for the storage
system, even over multiple market intervals, which common split market
clearing does not. More importantly, we showed that this property also stands
when forecast errors are made when calculating the final state of energy of
the storage system, which corresponds to a realistic setup. It still remains
to study how this final state of energy should be determined. This should also
come with a study of the impact of the storage level on pricing. In
particular, a critical next step would be to investigate the potential impact
of a strategic choice of this level on prices and social welfare, and how it
compares to having merchant storage. We also discuss the impacts of
uncertainty on social welfare compared to traditional approaches. In the case
of forcing the discharge of the storage system at a disadvantageous price, we
show that our approach comes closer to closing the gap with an ideal oracle
market clearing. We also showed that this method does not solve the problem of
accounting for prices in future market intervals, thereby leading to higher
prices compared to an ideal market clearing. This problem needs to be further
studied. Another limitation of the method is that it might happen that a
stored quantity is kept in store for too long, due to a very high value. We
showed that a discount factor could be applied to the value of stored energy
over time to avoid this. We used illustrative small-scale examples to give
better intuition on how the method introduced behaves compared to a
traditional split-horizon market clearing. Since this method does not
introduce non-linearities, the computational complexity is similar to
traditional approaches. It only involves solving one more linear program for
the update of storage values, which should also scale well. This procedure was
introduced on an idealized representation of a storage system, and the
promising results are a good motivation for extending it to more general
storage system models and to multiple storage systems.

## References

  * [1] IEA, “World Energy Outlook,” 2022, accessed 13-09-2023. [Online]. Available: https://www.iea.org/reports/world-energy-outlook-2022
  * [2] FERC, “Docket Nos. RM16-23-000 and AD16-20-000, Order No. 841: Electric Storage Participation in Markets Operated by Regional Transmission Organizations and Independent System Operators,” 2108, accessed 13-09-2023. [Online]. Available: www.ferc.gov/media/order-no-841
  * [3] N. G. Singhal and E. G. Ela, “Pricing impacts of state of charge management options for electric storage resources,” in _2020 IEEE Power & Energy Society General Meeting (PESGM)_, 2020, pp. 1–6.
  * [4] R. Sioshansi, “When energy storage reduces social welfare,” _Energy Economics_ , vol. 41, pp. 106–116, 2014.
  * [5] K. Hartwig and I. Kockar, “Impact of strategic behavior and ownership of energy storage on provision of flexibility,” _IEEE Transactions on Sustainable Energy_ , vol. 7, no. 2, pp. 744–754, 2016.
  * [6] J. A. Taylor, “Financial storage rights,” _IEEE Transactions on Power Systems_ , vol. 30, no. 2, pp. 997–1005, 2014.
  * [7] Y. Jiang and R. Sioshansi, “What duality theory tells us about giving market operators the authority to dispatch energy storage,” _The Energy Journal_ , vol. 44, no. 3, pp. 1–19, 2023.
  * [8] D. Muñoz-Álvarez and E. Bitar, “Financial storage rights in electric power networks,” _Journal of Regulatory Economics_ , vol. 52, no. 1, pp. 1–23, 2017.
  * [9] M. Weibelzahl and A. Märtz, “On the effects of storage facilities on optimal zonal pricing in electricity markets,” _Energy Policy_ , vol. 113, pp. 778–794, 2018.
  * [10] L. Frölke, E. Prat, P. Pinson, R. M. Lusby, and J. Kazempour, “On the efficiency of energy markets with non-merchant storage,” 2022. [Online]. Available: https://doi.org/10.48550/arXiv.2210.05305
  * [11] S. Bose and S. H. Low, _Some Emerging Challenges in Electricity Markets_. Springer International Publishing, 2019, pp. 29–45.
  * [12] L. Werner, A. Wierman, and S. H. Low, “Pricing flexibility of shiftable demand in electricity markets,” in _Proceedings of the Twelfth ACM International Conference on Future Energy Systems_ , 2021, pp. 1–14.
  * [13] G. Morales-España, R. Martínez-Gordón, and J. Sijm, “Classifying and modelling demand response in power systems,” _Energy_ , vol. 242, p. 122544, 2022.
  * [14] C. Chen and L. Tong, “Pricing real-time stochastic storage operations,” _Electric Power Systems Research_ , vol. 212, p. 108606, 2022.
  * [15] J. Zhao, T. Zheng, and E. Litvinov, “A multi-period market design for markets with intertemporal constraints,” _IEEE Transactions on Power Systems_ , vol. 35, no. 4, pp. 3015–3025, 2020.
  * [16] B. Hua, D. A. Schiro, T. Zheng, R. Baldick, and E. Litvinov, “Pricing in multi-interval real-time markets,” _IEEE Transactions on Power Systems_ , vol. 34, no. 4, pp. 2696–2705, 2019.
  * [17] J. Cho and A. Papavasiliou, “Pricing under uncertainty in multi-interval real-time markets,” _Operations Research_ , 2022.
  * [18] W. Zhang, P. A. Tominac, and V. M. Zavala, “Pricing and remunerating electricity storage flexibility using virtual links,” _Computers & Chemical Engineering_, vol. 168, p. 108019, 2022.