# An Event-based Algorithm for Simultaneous 6-DOF Camera Pose Tracking and
Mapping

Masoud Dayani Najafabadi, Mohammad Reza Ahmadzadeh The authors are with the
Department of Electrical and Computer Engineering, Isfahan University of
Technology, Isfahan, Iran (e-mail<EMAIL_ADDRESS>ahmadzadeh@iut.ac.ir).

###### Abstract

Compared to regular cameras, Dynamic Vision Sensors or Event Cameras can
output compact visual data based on a change in the intensity in each pixel
location asynchronously. In this paper, we study the application of current
image-based SLAM techniques to these novel sensors. To this end, the
information in adaptively selected event windows is processed to form motion-
compensated images. These images are then used to reconstruct the scene and
estimate the 6-DOF pose of the camera. We also propose an inertial version of
the event-only pipeline to assess its capabilities. We compare the results of
different configurations of the proposed algorithm against the ground truth
for sequences of two publicly available event datasets. We also compare the
results of the proposed event-inertial pipeline with the state-of-the-art and
show it can produce comparable or more accurate results provided the map
estimate is reliable.

###### Index Terms:

Event-based vision, event camera, KLT tracker, simultaneous localization and
mapping (SLAM), state estimation.

##  I Introduction

In applications such as robot navigation and control, augmented reality, and
3D reconstruction, we need to know a device’s location according to a map
representation of the environment. Sometimes we are interested in estimating
either the map or the robot’s pose. Localization or odometry is the problem of
estimating the robot’s pose when the map is known. On the other hand, the
general problem of estimating the device’s pose and the scene structure in
parallel is called simultaneous localization and mapping (SLAM).

The main challenge in an odometry/SLAM system is to achieve accurate and
robust performance in real-time, no matter how severe the conditions are.
Certain circumstances, including poor lighting, fast motion, dynamic scene,
and long outdoor distances, can adversely complicate the problem. This problem
is challenging because the estimation generally requires the optimization of
nonlinear equations with a large number of unknown parameters.

One way to address this issue is to use algorithmic techniques and constructs
such as parallelization of different parts of the system or the use of
efficient data structures [1], [2]. Another approach to achieve acceptable
accuracy and robustness against various conditions is to use a range of
complementary sensors.

A Dynamic Vision Sensor (DVS) such as [3] is a visual sensor in which each
pixel detects the log-intensity change at the pixel location and outputs an
_event_ independent of the neighboring pixels. Event cameras are ideal for
robotic applications because they allow low latency, high event generation
rates, high dynamic range, and low power consumption.

After surveying the related work in Section II, a short review of event data
and event generation mechanism is presented in Section III. We then explain
different components of the proposed algorithm in Section IV in detail.
Section V evaluates the performance of the proposed algorithm. Finally, we
conclude our discussion in Section VI.

##  II Related Work

Because our algorithm first represents events as intensity images and then
tracks detected features, we review both classical feature-based SLAM and
event-based algorithms. Our aim here is to consider only the most relevant
work rather than being comprehensive and complete.

PTAM [1] is a feature-based algorithm that separates concepts of localization
and mapping and parallelizes both tasks in concurrent threads. This technique
not only improves the overall performance but also provides more flexibility;
the estimated parameters from one module can be exploited in another without
waiting for the process to finish.

ORB-SLAM [2] gradually built on PTAM to address some issues. The last version
of this algorithm at the time of this writing, ORB-SLAM3 [4], supports a range
of sensors including regular cameras and the inertial measurement unit (IMU)
with different configurations (monocular, stereo, fisheye, distorted pinhole).
This algorithm uses the estimated map to localize pose, supports multi-session
mapping, can relocalize when tracking is lost, and reduce the accumulated
error with the _loop closing_ method. ORB-SLAM performance is limited due to
the specifications of regular cameras, despite its robustness and accuracy.

VINS-Mono [5] and its enhanced version, VINS-Fusion [6], are other feature-
based methods that fuse different sensors to estimate the pose and structure
using a window-based optimization scheme. Unlike ORB-SLAM, which associates
features with their descriptors, VINS-Fusion is a KLT-based pipeline. Although
it outperforms ORB-SLAM in some situations [4], it still suffers from the same
restrictions.

In [7], Gallego et al. surveyed event-based vision and its applications. In
the case of event-based odometry and SLAM, researchers have begun by
simplifying assumptions for the camera motion or scene structure. For
instance, Gallego et al. [8] proposed an algorithm for angular velocity
estimation in which only pure 3D rotations were considered. However, most
practical applications require more degrees of freedom, and a broader motion
model should be considered.

Mueggler et al. [9] proposed a continuous framework for pose estimation based
on event data only. The information of every individual event is involved in
state estimation. As noted in [10], each event contains little information, so
it is better to process a group of events. Additionally, some implementations
of the continuous algorithms demand estimating many parameters in a very
narrow time window.

In [11], Kim et al. reconstructed intensity images using event data and the
features were detected and tracked in these frames. Although the algorithm
performs well under fast motion and poor lighting conditions, it needs GPUs
for real-time performance.

Rebecq et al. [12] used IMU measurements and events to track 6-DOF camera
motion in a window-based optimization framework. They reconstructed motion-
compensated images from consecutive overlapped spatiotemporal event windows.
Then the features were extracted and tracked using the KLT [13] method.
Similarly, Vidal et al. [10] further utilized intensity images. Unlike [12],
event windows were synchronized with the timestamp of intensity images. This
synchronization can obscure the advantages of asynchronous low-latency event
data. Both schemes require IMU measurements to initialize motion parameters
and reconstruct frames. They also manually set the event window size for each
sequence, which restricts the flexibility of their algorithm.

Similar to our algorithm, the event-based visual-inertial odometry pipeline
(EVIO) in [14] adaptively selects the best spatiotemporal event window length
based on the events’ optical flow. However, while the proposed method depends
only on input events in its basic configuration and addresses a range of
conditions, EVIO relies on IMU readings for state estimation and does not
explicitly consider situations where, for example, the camera is motionless.

In this paper, we propose an event-based SLAM algorithm with the following
features:

  * •

A novel image reconstruction algorithm adaptively selects the event window
size based on camera motion and scene structure and converts them into a
motion-compensated image.

  * •

A higher-level KLT-based localization and mapping module exploits map data
derived from MC images to estimate the current pose of the camera.

  * •

An event-inertial version of the proposed event-only pipeline is presented to
show how additional sensors can improve the algorithm.

  * •

Distorted pinhole and Kannala Brandt [15] camera models are supported, the
camera motion has six degrees of freedom (6-DOF), and there is no assumption
about the type of scene.

  * •

The proposed algorithm has robust performance in different conditions as long
as the map estimate is reliable and there is relative motion between the
camera and the scene.

##  III Event Data

If there is a relative motion between the DVS and the scene, the illumination
$I$ at each pixel location changes. When there is enough change in the log-
intensity, $L(t_{k})=\log(I(t_{k}))$ at current time $t_{k}$, relative to a
reference time $t_{r}$, it outputs an event at each pixel coordinate
$\boldsymbol{x}_{k}=(x_{k},y_{k})$. Formally, if the change in $L$ is greater
than a threshold, $C$,

$|\Delta L(t_{k},t_{r})|=|L(t_{k})-L(t_{r})|>C$ (1)

there is an event

$\boldsymbol{e}_{k}:\\{t_{k},x_{k},y_{k},p_{k}\\}$

where $p_{k}$ shows the sign of change, i.e., if $L(t_{k})$ is bigger than
$L(t_{r})$, $p_{k}$ is positive, and it is negative otherwise.

Fig. 1 shows a slice of all events in the sequence shapes_6dof from the Public
Event Dataset [16]. Each 3D point in the spatiotemporal space in Fig. 1(1(a))
designates an event, color-coded based on event polarities. Each slice of this
space represents a spatiotemporal event window.

(a)

(b)

Figure 1: Representation of a spatiotemporal event window from shapes_6dof
[16]. (1(a)) Spatiotemporal window of 2000 events starting from a specific
timestamp. Blue and red points show events with positive and negative
polarity, respectively. (1(b)) 2D view of (1(a)) in the image ($x$-$y$) plane.

Since event generation depends on intensity change in each pixel, there is no
spatial correlation between adjacent pixels. Each pixel generates events
asynchronously as soon as a change in its intensity value is detected. On the
other hand, uniform areas in the image cannot contribute to event generation.
This redundancy reduction allows lower latencies and faster event generation
rates. Events compactly summarize information regarding the scene’s high-
contrast areas.

Note how the event generation rate depends on the relative motion between the
camera and the scene and structure in normal illumination conditions. For a
fixed threshold, $C$, if the camera traverses slowly through the environment,
the event generation rate is low, and event data is noisy. On the other hand,
in highly textured scenes, the event generation rate could be high at regular
speeds. For specific camera settings, if the camera moves fast in a textured
environment, the event generation rate eventually reaches the maximum
allowable output bandwidth of the event camera.

We can benefit from available image-based SLAM algorithms by representing
event data as a 2D image. To produce the event histogram $I_{r}$ for a
specific event window $W_{r}$ at reference time $t_{r}$, all events in each
pixel location are added

$I_{r}(\boldsymbol{x})=\sum_{\boldsymbol{e}_{k}\in
W_{r}}p_{k}\delta(\boldsymbol{x}-\boldsymbol{x}_{k}^{\prime}),$ (2)

where $\delta(.)$ is the continuous Dirac’s delta kernel. It is common to
replace the continuous delta operator with a discrete-time sampled Gaussian
kernel with an arbitrary standard deviation $\sigma_{I}$. The parameter
$p_{k}$ can be changed based on event polarity or set to constant 1 for all
events in the window.

Depending on the camera speed and the length of the event window, the event
histograms may suffer from motion distortion. In this case, event coordinates
can be corrected before reconstructing these images. The resulting image is
called a Motion-Compensated Image or MCI for short.

For a carefully selected event window, we can mitigate the motion distortion
using a two-dimensional similarity or rigid body transformation [Sim(2) or
SE(2)]. Assuming a constant rotational speed $\omega$ and translational speed
$\boldsymbol{v}=(v_{x},v_{y})$ for each event $\boldsymbol{e}_{k}$ in the
window, the warp

$\begin{bmatrix}x^{\prime}\\\
y^{\prime}\end{bmatrix}=s\begin{bmatrix}\cos\theta_{k}&-\sin\theta_{k}\\\
\sin\theta_{k}&\cos\theta_{k}\end{bmatrix}\begin{bmatrix}x\\\
y\end{bmatrix}+\begin{bmatrix}t_{k,x}\\\ t_{k,y}\end{bmatrix}$ (3)

maps the event location $\boldsymbol{x}$ at $t_{k}$ to
$\boldsymbol{x}^{\prime}$ at time $t_{r}$ where $\theta_{k}=\omega\Delta t$,
$\boldsymbol{t}_{k}=\boldsymbol{v}\Delta t$, $s$ is the arbitrary scale, and
$\Delta t=|t_{k}-t_{r}|$. Event coordinates should be normalized based on the
camera intrinsics beforehand.

When the relative transform $T_{t_{r},t_{k}}$ and the depth
$Z(\boldsymbol{x}_{k})$ of each event in the window are known, a SE(3) mapping
can be used to warp events:

$\boldsymbol{x}_{k}^{\prime}=\pi_{0}(T_{t_{r},t_{k}}[Z(\boldsymbol{x}_{k})\pi_{0}^{-1}(\boldsymbol{x}_{k})]).$
(4)

Here $\pi_{0}(.)$ projects the 3D point to the image plane viewed from the
current pose, and $\pi_{0}^{-1}(.)$ is its inverse operation. In this case, we
first retrieve the 3D feature point associated with each event, seen from the
current camera position, and map it to a reference location.

If we extract feature points from the reference reconstructed frame and track
them across the current image, it is possible to estimate the parameters of
(3) and (4) for all events between these frames. For this, we minimize the
reprojection error between each warped reference feature location
$\boldsymbol{x}_{i}$ and the corresponding match $\boldsymbol{x}_{j}$ in the
current frame using

$\boldsymbol{\theta}^{*}=\underset{\boldsymbol{\theta}^{*}}{\text{argmin.}}\sum_{i,j\in
J}\boldsymbol{e}_{ij}^{T}(\boldsymbol{x})\boldsymbol{\Omega}_{ij}\boldsymbol{e}_{ij}(\boldsymbol{x}),$
(5)

where $\boldsymbol{\theta}^{*}$ is the transformation parameters,
$\boldsymbol{e}_{ij}=\boldsymbol{x}_{i}-\boldsymbol{x}_{j}$,
$\boldsymbol{\Omega}_{ij}$ is the information matrix, and $J$ is the set
containing all matched feature pairs between the two frames.

To use the estimated parameters $\boldsymbol{\theta}^{*}$ to warp events in
the current window, we first convert the relative transform to speed,
$\boldsymbol{\psi}$, using

$\boldsymbol{\psi}=\frac{\text{Log}(\boldsymbol{\theta}^{*})}{\Delta T},$

where $\Delta T$ is the length of the window in units of time and
$\text{Log}(.)$ is the inverse algebra that maps an element of SE(2) or SE(3)
to a member of the corresponding tangent space. On the premise that the camera
moves slow enough to assume constant speed for all events, the unknown
transform between each event location and the reference timestamp,
$T_{t_{r},t_{k}}$, is

$T_{t_{r},t_{k}}=\text{Exp}(\boldsymbol{\psi}\Delta t),$

where $\Delta t=|t_{k}-t_{r}|$ and $\text{Exp}(.)$ is the exponential map for
the corresponding Lie group.

We also distinguish the notion of an image and a frame. While an image is
simply a two-dimensional array of numbers, a frame consists of other
information, including the timestamp, camera pose, and key points beside the
image.

Interested readers can refer to [7] for an in-depth review of the event
generation mechanism and event representation.

##  IV Algorithm

Figure 2: An overview of the proposed event-based SLAM pipeline.

Fig. 2 illustrates an overview of the proposed event-based SLAM algorithm. Our
pipeline consists of two main components, which run concurrently in separate
threads. The first part selects an appropriate spatiotemporal window of input
events and reconstructs an image for each event window. The second part of the
proposed algorithm extracts features from each input image and tracks them
using the KLT method across the following frames. The 6-DOF pose of the camera
is then estimated, and the local scene is reconstructed. Sections IV-A and
IV-B discuss each component in more detail.

###  IV-A Event Window Selection and Image Reconstruction

The basic idea of the algorithm is to process events in small chunks and
generate an MCI whenever there are enough events. The algorithm uses two sets
of event windows to select the appropriate window size and reconstruct the
MCI. It continually tracks and accumulates _tiny windows_ to determine whether
there are enough events. It reconstructs the MCI using the _reconstruction
window_ consisting of all previously collected events. Fig. 3 illustrates the
relationship between the event stream, tiny windows, and reconstruction
windows.

Figure 3: An illustration of the input event stream, tiny event windows, and
non-overlapped reconstruction windows relative to each other.

This algorithm first selects a spatiotemporal window of input events.
Initially, the length of the window, $N_{e}$, is selected arbitrarily and can
be adjusted later. At this stage, consecutive event windows have no overlap.
The proper event window size depends on the event generation rate. As
discussed in Section III, the maximum bandwidth of the event camera limits the
event generation rate. For example, in an event camera with a bandwidth of 1
million events per second and a resolution of $240\times 180$ pixels, we can
expect up to around 23 events per second per pixel. Thus the minimum event
generation rate can be chosen based on the bandwidth of the event camera,
e.g., one event per pixel per second. We select this threshold solely based on
the event camera specifications agnostic to the camera motion and scene.

Based on these observations, the event generation rate is calculated by

$r=\frac{N_{e}}{\Delta tWH}$ (6)

where $\Delta t=t_{e_{N_{e}}}-t_{e_{1}}$ is the difference between the first
event timestamp, $t_{e_{1}}$, and the last timestamp, $t_{e_{N_{e}}}$ in the
current window, and $W$ and $H$ are the width and height of the event image in
pixels, respectively.

If the camera does not move fast enough and the event generation rate is low
then the event data will be noisy. The algorithm rejects most tiny windows
with an event generation rate less than a threshold, $th_{e}$, and restarts
from the initial state. The only exception is when there are at least $N_{f}$
frames. In this case, the algorithm restarts after it generates and dispatches
the MCI.

If the event generation rate of the current tiny window is acceptable, a tiny
frame is reconstructed based on (2) without motion compensation. The following
sections review the main components of the image reconstruction algorithm.

####  IV-A1 KLT Optical Flow Initialization and Tracking

In each iteration, the algorithm detects FAST features [17] in the first tiny
frame. Since we use image features to determine the displacement in event
locations, fewer features are sufficient. If we cannot detect enough features,
this step is repeated with the subsequent tiny frames; otherwise, we proceed
to the next step.

The algorithm then tracks the reference features in the subsequent frames. We
can use the feature matches obtained from the feature tracks for the two-view
reconstruction. The map initialization algorithm in [2] is used to estimate
the relative pose between the current and the reference frames and the 3D map
points corresponding to the feature matches. Due to the small baseline between
the consecutive tiny frames, most scene reconstruction algorithms cannot yield
reliable results. So, we relieve some of the stringent conditions in [2],
i.e., the best model does not have to stand out by a large margin. Because
most implementations use the RANSAC method internally, feature matches are
still enhanced despite potential failures in scene reconstruction. If a
successful reconstruction of the 3D structure is possible, this information is
stored and used in the MCI reconstruction step.

We define the median feature displacement metric as the median distance of the
tracked image features in pixels. Whenever the median feature displacement in
the current tiny frame is higher than a threshold, e.g., five pixels, enough
events have been accumulated, and the upcoming MCI reconstruction is
commenced.

Note that depending on the type of the scene, the number of detected features
varies, but the median feature displacement is not affected. We select the
feature displacement threshold based on how much distortion we can tolerate in
event locations. This threshold is constant during the operation of the
algorithm for all scenes.

After the reconstruction and dispatching of the MCI, the algorithm restarts
from the next iteration. Before it begins the next iteration, it updates the
size of the tiny windows, $N_{e}$. If $N_{f}$ is the number of processed
frames and $N_{x}$ is the number of expected frames per iteration, the new
$N_{e}^{*}$ is

$N_{e}^{*}=\big{\lfloor}\frac{N_{f}N_{e}}{N_{x}}\big{\rfloor},$

where $\lfloor x\rfloor$ is the greatest integer less than or equal $x$.

In this case, $N_{e}$ is adjusted according to the camera speed automatically.
If the camera speed increases, fewer events satisfy the feature displacement
condition, and $N_{e}$ decreases. On the other hand, if the camera moves
slowly, we need more events for a defined feature displacement, and $N_{e}$
increases. Since we did not assume any scene structure, $N_{e}$ is independent
of the type of the scene.

####  IV-A2 Motion Compensated Images

We use different warps to rectify motion distortion in event coordinates and
use (2) to reconstruct the MCI. We consider the last event timestamp as the
reference and project preceding events in the forward direction.

If there is a reliable 3D reconstruction, we exploit this information in (4)
to warp the event locations. As in [10], the median depth of 3D map points is
used for all events in the window. Prior to this operation, motion and map
parameters can be enhanced using the Bundle Adjustment (BA) method.

On the other hand, if there exists no initial scene and motion estimate, but
we have tracked the features across several tiny frames, we can still
undistort events using a 2D optimization scheme. In this case, we first fit a
Sim(2) or SE(2) motion model to the feature tracks and then reproject events
to the reference frame according to (3). We also store this model as
initialization parameters for the subsequent iterations.

In some situations, only the relative pose estimate is known. When the camera
motion is adequately slow, we can either integrate the IMU measurements or use
the last pose estimate to infer the current camera location. In either case,
we can estimate the average depth of the scene using (5) and warp event
locations based on (4).

We run each reconstruction method concurrently in parallel threads. Besides
the MC images, we reconstruct the event histogram using (2) with no motion
compensation. Finally, we score each image based on its contrast or
_sharpness_ [18]. To measure image sharpness, we first divide each image into
non-overlapped patches and compute the local standard deviation (STD) for each
patch. We select and dispatch the image with the highest average local STD. We
choose the local STD because it yields more stable results than the global
STD.

The KLT algorithm has a poor performance when there is a large baseline
between consecutive frames [13]. Although there is not much motion between
tiny frames, the distance between the reconstructed frames can be much higher.
One way to overcome this issue is to send tiny frames in addition to MCIs to
the event-based tracking algorithm (Section IV-B).

Another approach to ensure smooth transition between consecutive MC frames is
to overlap the MC windows. We prefer this method to the non-overlapped scheme
because it results in more robust feature tracking in the second level.

###  IV-B Event-based Localization and Mapping

With the MCI from the previous step, the event-based localization module
follows an image-based KLT optical flow scheme to estimate the 6-DOF pose of
the camera and reconstruct the scene. The basic steps of this algorithm are
almost similar to the reconstruction algorithm discussed in Section IV-A.

####  IV-B1 Initialization

If there are not enough feature tracks, the algorithm detects new FAST
features in the input MCI. It then creates and merges new feature tracks to
manage feature locations across successive images. We use a bucketing grid
scheme to ensure a uniform distribution of detected features across images.

The last features are tracked in each subsequent image using the KLT method.
If the map is not initialized, we use the two-view reconstruction algorithm to
reconstruct the scene and recover the camera’s relative motion between the
current and reference frames. As in Section IV-A, we use the map
initialization algorithm in [2] for this step. After a successful
reconstruction, we initialize the map and perform a global bundle adjustment
optimization to enhance the estimate.

####  IV-B2 Tracking The Local Map

We can use feature tracks to identify feature and map point associations
across different frames. Therefore, the geometrical distance between the
projected map points in the current frame and the corresponding key point
observations is minimized to estimate the current pose of the camera.

Due to the camera motion, we might lose track of some features. In this case,
the algorithm fails to track features that are not visible in the subsequent
frames.

####  IV-B3 Key Frame Insertion and Local Mapping

The algorithm inserts a new keyframe whenever there is enough baseline between
the current tracked frame and the last keyframe. We use several heuristics to
decide when to spawn a new keyframe. Whenever the number of tracked map points
falls below a threshold or the median pixel displacement between key points in
the current frame and the preceding keyframe is greater than a certain limit,
a new keyframe is created.

The local mapping module is very similar to the method used in ORB-SLAM. While
our algorithm uses the KLT method to track and associate features, ORB-SLAM
detects new matches with descriptors.

In local mapping, we first cull outlier map points based on their rate of
observations. Additional map points are then triangulated for the new feature
matches between the current and previous keyframes. Finally, a local bundle
adjustment optimization enhances the map point estimates.

###  IV-C IMU Measurements

In this section, the necessary changes to incorporate IMU measurements are
introduced. The preintegration theory and the IMU initialization module of
[19] and [2] are used to accomplish this goal. These changes are as follows:

  * •

IMU measurements are integrated between the tiny frames, the reconstruction
frames, and the keyframes.

  * •

In the MCI reconstruction module, using the initial pose from the IMU
integrations and feature matches, we first find the average depth of the scene
in the 3D warp (4) through the optimization framework of (5). We also consider
IMU biases and the gravity direction as the state variables in this
optimization. The initial gravity direction is evaluated using the first
accelerometer reading for the current iteration. These parameters are stored
for the subsequent iterations, and the MCI is reconstructed.

  * •

We use an inertial bundle adjustment optimization framework in the MCI
generation step with a successfully reconstructed scene. Besides IMU biases
and the gravity direction, we also consider the relative scale.

  * •

We adapt the inertial local bundle adjustment and the mechanism of IMU
initialization discussed in [2] in the local mapping algorithm. In summary,
the IMU biases, the direction of the gravity, and the relative scale are
estimated using an adequate number of keyframes. Such optimization and scale
refinement are repeated for several predefined periods.

##  V Experiments

To evaluate the performance of the proposed pipeline, we use two publicly
available event datasets, the Public Event Dataset [16] and the Multi-vehicle
Stereo Event Dataset [20]. The former includes short-duration sequences
recorded in different scenes and challenging conditions by a monocular
240$\times$180 pixel DAVIS device. An embedded IMU measures acceleration and
angular velocity along three axes of freedom. Most sequences include events,
intensity images, IMU measurements, and the ground truth. We only consider the
ones for which both the ground truth and IMU measurements are available; hence
we exclude sequences from the “depth” group.

Sequences of [20] are recorded with two similar DAVIS 346$\times$260 pixel
cameras attached in the stereo configuration. It contains events and intensity
images for both cameras, IMU measurements, the ground truth trajectory, and
depth maps for each frame. We only use the left camera and several
representative outdoor sequences, noting that the proposed pipeline is not
designed to tackle the stereo case or track longer outdoor distances.

We use several objective criteria to verify the results of the proposed
algorithm quantitatively. Since our algorithm generally results in an atlas of
multiple disconnected pose graphs, we extend the Relative Pose Error (RPE) in
[21] to average the normalized relative errors as

$\overline{RPE}(\boldsymbol{\theta})=\frac{1}{|\mathcal{A}|}\sum_{\mathcal{G}\in\mathcal{A}}\frac{1}{D|\mathcal{G}|}\sum_{i,j\in\mathcal{G}}\boldsymbol{\theta}(\boldsymbol{s}(\Delta
T^{*}_{i,j})\ominus\Delta T_{i,j})$ (7)

where $\mathcal{A}$ is the set of all pose graphs, $\mathcal{G}$, $\Delta
T_{i,j}=T_{i}\ominus T_{j}$ is the relative transform between pairs of SE(3)
poses in $\mathcal{G}$, $T_{i}$ and $T_{j}$, $\ominus$ is the inverse of SE(3)
Lie algebra, $D$ is the total traversed distance for the current graph,
$\boldsymbol{s}()$ is the scaling operation for monocular event-only tracking,
$\boldsymbol{\theta}()$ returns either rotational or translational component
of error, and ∗ indicates the estimate. We normalize RPE by the total
traversed distance to reflect the effects of variable length pieces. For the
monocular event-only configuration, we calculate the unknown scale by
comparing the associated estimate and ground truth pairs and scale the
estimate before computing the average error.

Additionally, we consider measures to assess the stability of the pipeline.
Total traversed distance or time is the sum of all distances or delta times
between consecutive frames in each pose graph.

Before the algorithm processes input events, lens distortions in event pixel
locations are rectified using available calibration parameters for each
sequence. We choose $N_{x}=3$, and the initial value of $N_{e}$ is 2000 events
for sequences of [16] and 6000 for [20]. Since the event frames in [20] have a
higher resolution, a higher value for $N_{e}$ helps the convergence speed,
though starting from 2000 events should eventually converge to the optimal
value. Motion-compensated windows have 50% overlap, and Event frames are
reconstructed using a Gaussian kernel with $\sigma_{I}=1$ in both tracking
levels (tiny frames and MC images). The FAST feature detector threshold is set
to a small amount (around zero) because this setting yields the most features
and allows us to filter them by their response. We set the KLT tracker with
two pyramid levels, a block size of 23$\times$23 pixels, and a maximum
bidirectional error of one. These settings are chosen arbitrarily based on our
experiments.

TABLE I: Performance evaluation of the proposed algorithm for two event-based configurations compared against the ground truth for some sequences of [16] and [20]. |  | E-Only |  | E-I-C1
---|---|---|---|---
|  | $\overline{\text{{{RPE}}}}$ | $\overline{\text{{{RPE}}}}$ |  |  | $\overline{\text{{{RPE}}}}$ | $\overline{\text{{{RPE}}}}$ |
Dataset | Sequence | Position | Rotation | Stability |  | Position | Rotation | Stability
|  | (-) | (deg/m) | ($\times 10^{3}\text{m.s}$) |  | (-) | (deg/m) | ($\times 10^{3}\text{m.s}$)
| shapes_6dof | 0.048 | 0.209 | 2.31 |  | 0.050 | 0.114 | 2.50
| shapes_translation | 0.084 | 0.186 | 2.29 |  | 0.092 | 0.169 | 2.62
Public | poster_6dof | 0.017 | 0.062 | 3.18 |  | 0.018 | 0.057 | 3.14
Event | poster_translation | 0.005 | 0.006 | 2.29 |  | 0.005 | 0.006 | 2.33
[16] | hdr_poster | 0.011 | 0.025 | 2.67 |  | 0.002 | 0.007 | 2.73
| boxes_6dof | 0.002 | 0.004 | 3.91 |  | 0.006 | 0.008 | 3.95
| boxes_translation | 0.004 | 0.004 | 3.40 |  | 0.004 | 0.004 | 3.46
| hdr_boxes | 0.004 | 0.007 | 3.08 |  | 0.004 | 0.006 | 3.07
| dynamic_6dof | 0.019 | 0.041 | 2.09 |  | 0.032 | 0.117 | 2.11
| dynamic_translation | 0.012 | 0.018 | 1.24 |  | 0.034 | 0.032 | 1.38
| indoor_flying1 | 0.044 | 0.041 | 1.12 |  | 0.080 | 0.046 | 1.34
Multi- | indoor_flying2 | 0.041 | 0.032 | 1.86 |  | 0.029 | 0.019 | 2.11
Vehicle | indoor_flying3 | 0.018 | 0.006 | 3.23 |  | 0.029 | 0.008 | 3.63
Event | indoor_flying4 | 0.038 | 0.013 | 0.11 |  | 0.038 | 0.016 | 0.12
Stereo | outdoor_day1 | 0.889 | 0.033 | 3.14 |  | 1.181 | 0.041 | 6.38
[20] | outdoor_night1 | 0.678 | 0.003 | 3.29 |  | 1.131 | 0.008 | 3.94
| outdoor_night3 | 1.208 | 0.009 | 5.04 |  | 1.406 | 0.017 | 8.95

Table I summarizes the evaluation results for event-only (E-Only) and event-
inertial (E-I-C1) sensor configurations. The stability column of this table is
the multiplication of the total traversed distance in meters and time in
seconds. Although event-only tracking shows similar or superior results in
most cases, the inertial method is more stable. Since more challenging periods
can be tracked in the inertial case, the errors grow accordingly. Furthermore,
because we do not scale inertial pose graphs, the reported results also
include scaling errors.

Besides the challenging conditions of each sequence that can affect the
quality of the event image reconstruction and the accuracy of pose estimates,
we also note two degrading effects. For lower resolution images of [16], the
rays connecting the map points with their respective observations in each
frame are so close that the uncertainty in the depth of the subsequent map
point estimates grows [22]. Since our pipeline estimates the following poses
according to the map, the distance between consecutive poses diminishes, which
in turn results in less accurate map estimates. This cycle continues until the
tracking is lost. As discussed in [23], this situation can also happen for
larger images of [20] recorded by the fisheye camera due to inaccurate
translation estimates caused by far map points, especially for outdoor
sequences.

Although in the inertial case, the scale information embedded in IMU readings
prevents the _shrinking_ effect, inaccuracies due to IMU biases can affect map
point estimates. If there are not enough map points, the following pose and
map point estimates are affected by the growing bias of inertial measurements.
In this case, subsequent map point estimates continue to expand after the
insertion of each keyframe until the tracking is lost.

The inertial configuration discussed in Table I (E-I-C1) still assumes that
the tracking can be lost due to severe conditions, so it tries to reinitialize
the map as soon as possible. To further assess the limitations of the proposed
algorithm and compare it against the inertial state estimation in [10], we
consider two other configurations (E-I-C2 and E-I-C3). We configure our
pipeline to enforce continuous tracking without spawning and initializing new
maps in unfavorable conditions. In both cases, we disallow the
reinitialization of the map and continue to estimate 3D points using inertial
readings and feature matches between the nearby frames. In the last
configuration (E-I-C3), we also fix the tiny event window size to a predefined
value and restrict the withdrawal of noisy frames based on the event
generation rate. Table II summarizes our preferences for the length of fixed
tiny windows for each group of sequences in [16]. We select these values based
on our experience of E-I-C1 with variable window sizes.

TABLE II: Fixed tiny event window size specification for each group of sequences in [16]. Sequence | Number of Events
---|---
shapes | 2000
dynamic | 8000
poster | 10000
boxes | 12000
TABLE III: Comparison of the performance of different inertial configurations of the proposed algorithm for sequences of [16]. | E-I-C1 |  | E-I-C2 |  | E-I-C3
---|---|---|---|---|---
| $\overline{\text{{{RPE}}}}$ | $\overline{\text{{{RPE}}}}$ |  | $\overline{\text{{{RPE}}}}$ | $\overline{\text{{{RPE}}}}$ |  | $\overline{\text{{{RPE}}}}$ | $\overline{\text{{{RPE}}}}$
Sequence | Position | Rotation |  | Position | Rotation |  | Position | Rotation
| ($\times 10^{-3}$) | ($\times 10^{-3}$ deg/m) |  | ($\times 10^{-3}$) | ($\times 10^{-3}$ deg/m) |  | ($\times 10^{-3}$) | ($\times 10^{-3}$ deg/m)
shapes_6dof | 50 | 114 |  | 1 | 6 |  | 2 | 6
shapes_translation | 92 | 169 |  | 1 | 3 |  | 1 | 2
poster_6dof | 18 | 57 |  | 1 | 6 |  | 2 | 7
poster_translation | 5 | 6 |  | 4 | 5 |  | 4 | 5
hdr_poster | 2 | 7 |  | 3 | 7 |  | 3 | 7
boxes_6dof | 6 | 8 |  | 1 | 3 |  | 1 | 2
boxes_translation | 4 | 4 |  | 2 | 4 |  | 2 | 4
hdr_boxes | 4 | 6 |  | 2 | 6 |  | 2 | 6
dynamic_6dof | 32 | 117 |  | 2 | 5 |  | 2 | 5
dynamic_translation | 34 | 32 |  | 4 | 5 |  | 4 | 5

Table III compares different inertial configurations using the same RPE
metrics of Table I. Both E-I-C2 and E-I-C3 configurations show a significant
performance boost to the E-I-C1 configuration. In this case, the prolonged
tracking period of E-I-C2 and E-I-C3 configurations provokes more IMU
initialization and refinement steps and improves the accuracy of the IMU
parameters and the whole system.

The last two columns of Table III also investigate the efficacy of the
adaptive selection of the tiny window size. Although the results are similar
in most sequences, noting that the window length in E-I-C3 is fixed to the
optimal value, a wrongly selected size can impact the performance of the
system. Therefore, we report the E-I-C2 configuration as our preferred method.

Next, we compare our pipeline against the most relevant algorithm that
provides a stable open-source implementation. For this reason, we exclude the
methods that utilize the intensity images with the events.

(a)

(b)

Figure 4: The comparison between the proposed algorithm (triangles, solid
lines) and the event-inertial method of [10] (circles, dashed lines) in terms
of the average relative pose error over a range of distances for the color-
coded sequences of [16]: (4(a)) average relative translation error, (4(b))
average relative yaw error.

Fig. 4 contrasts the E-I-C2 configuration of the proposed algorithm against
the event-inertial method in [10]. We run the latest open-source release of
their algorithm111https://github.com/uzh-rpg/rpg_ultimate_slam_open for each
sequence in [16] using the default configurations (the sequence poster_6dof is
excluded because this implementation of [10] fails to produce reliable
results). Similar to [10], we measure and report the average relative
translation and yaw errors over a range of distances. We calculate the
relative pose error between the first one hundred pairs that meet a specific
distance range.

Although the window-based state estimator of [10] estimates the pose of the
camera and map points simultaneously, our algorithm relies on an accurate map
to find the next pose. As a result, any condition that degrades map estimate
accuracy will affect our results. For the sequence shapes_6dof, when the
camera moves fast in front of a low-textured poster, the tracking fails due to
the lack of detected features, and errors grow exponentially. Even though the
scene in sequence dynamic_6dof is textured, the change in 3D feature locations
affects the performance of the proposed method. Despite the poor lighting
condition in hdr_boxes, there are enough reliable map points, and the proposed
algorithm outperforms the state-of-the-art based on the translation errors.

Based on our experiments, while the state estimator of [10] can produce more
accurate results locally, the proposed algorithm can outperform over longer
distances as long as the map estimates are reliable. As mentioned before, when
the map estimates are valid, it will affect the subsequent pose estimates even
when some map points are no longer visible. Consequently, although the results
are not explicitly shown here, we claim that the proposed pipeline can reduce
the absolute pose error (APE) when these conditions are met.

Figure 5: The absolute estimated positions of the proposed algorithm and the
event-inertial method of [10] are compared against the ground truth in terms
of the APE for the first 15 seconds of the sequence boxes_6dof.

Fig. 5 compares the absolute estimated position of the results of Fig. 4 for
the sequence boxes_6dof. We use the APE tools in [2] to align the first 15
seconds of the estimates with the ground truth. In this case, the proposed
method can track the ground truth more accurately.

TABLE IV: Average tracking time statistics in milliseconds for three main modules of the proposed algorithm, calculated for representative sequences of [16] and [20]. | Event-Only |  | Event-Inertial-C1
---|---|---|---
| MTH | MTH | L1 | L2 |  | MTH | MTH | L1 | L2
Sequence | /MCF | /TF | /TF | /TF |  | /MCF | /TF | /TF | /TF
shapes_6dof | 72 | 48 | 12 | 6 |  | 82 | 56 | 13 | 7
poster_6dof | 158 | 106 | 26 | 7 |  | 173 | 119 | 27 | 9
boxes_6dof | 158 | 106 | 26 | 7 |  | 174 | 118 | 28 | 9
dynamic_6dof | 138 | 92 | 23 | 7 |  | 157 | 107 | 25 | 9
indoor_flying1 | 139 | 93 | 23 | 8 |  | 147 | 99 | 24 | 15
outdoor_day1 | 171 | 115 | 28 | 5 |  | 185 | 127 | 29 | 10
outdoor_night1 | 187 | 125 | 31 | 5 |  | 202 | 138 | 32 | 14

Although we have not tried to optimize the performance of the proposed
algorithm, we present the average tracking time statistics to contrast the
computation cost of each module. In particular, we measure and report the
performance of the main tracking thread (MTH), the first-level event image
reconstruction module (L1), and the second-level event frame tracking thread
(L2). Since it almost takes several tiny frames in L1 and one iteration of L2
for the main tracking thread to process input events, we include the timing as
per motion-compensated frames (MCF) and as per tiny frame (TF). Assuming there
are three TF for each MCF on average, we subtract double L1 values from the
MCF results and report the per TF results. The timing statistics for the local
mapping thread are not reported here. We run our pipeline on a system with an
Intel Core-i7 9700K CPU, 64 GB RAM, and Ubuntu 20.04 LTS operating system.

Table IV summarizes the timing results of different components of our pipeline
for two sensor configurations and some representative sequences of both
datasets [16] and [20]. The overall performance of our algorithm is several
times lower than the real-time. Note that more textured sequences demand more
processing power. The addition of the IMU slightly increases the performance
cost of our method. Based on these results, the reconstruction of each MCF is
generally the most costly operation. It takes about the same amount of time to
reconstruct an MCF, around 28 ms on average, consistent across all
configurations and datasets. The timing results of L2 are only comparable
across the Public Event dataset. For sequences of [20], the cost of the
inertial L2 module is almost twice the event-only case. One reason for this
difference could be the high depth variation inherent in the fisheye camera
and the difficulty of scale estimation and refinement involved in inertial
tracking in this case (although the local mapping is performed in a separate
thread, the inertial L2 must wait because of the map change).

##  VI Conclusion

We proposed an algorithm that reconstructs MC images from the adaptively
selected event windows and uses them to estimate the structure and trajectory.
The first module tracks FAST features across multiple tiny frames to choose
the best event window size and resolve the parameters needed for MCI
generation. After reconstruction, it sends the best MCI representation to an
image-based SLAM to initialize the map and track the 6-DOF pose of the camera
in it. Furthermore, we showed how to utilize inertial measurements to improve
the performance of the event-only algorithm.

We compared the estimated trajectory of different configurations of the
proposed pipeline with the ground truth for sequences of two publicly
available event datasets. Although the algorithm produces accurate results in
most cases, its performance depends on the reliability of the map estimate. We
also showed that the proposed algorithm outperforms the state-of-the-art as
long as a valid map is available.

Although the event window selection module is costly, it can still be used in
most event-based algorithms to select the best event window without
reconstructing the MCI. It is possible to use different techniques such as
relocalization or loop-closing to improve the accuracy of the map, the overall
performance, and the algorithm’s stability. These techniques generally extract
and use descriptors to link different segments of a pose graph. One can
consider event-based descriptors for that purpose or extract descriptors from
regular intensity images, which are processed along with the event stream.

## References

  * [1] G. Klein and D. Murray, “Parallel tracking and mapping for small ar workspaces,” in 2007 6th IEEE and ACM International Symposium on Mixed and Augmented Reality, pp. 225–234, 2007.
  * [2] R. Mur-Artal, J. M. M. Montiel, and J. D. Tardos, “ORB-SLAM: A versatile and accurate monocular SLAM system,” IEEE Transactions on Robotics, vol. 31, no. 5, p. 1147–1163, 2015.
  * [3] P. Lichtsteiner, C. Posch, and T. Delbruck, “A 128$\times$128 120 db 15 $\mu$s latency asynchronous temporal contrast vision sensor,” IEEE Journal of Solid-State Circuits, vol. 43, no. 2, pp. 566–576, 2008.
  * [4] C. Campos, R. Elvira, J. J. G. Rodríguez, J. M. M. Montiel, and J. D. Tardós, “ORB-SLAM3: An accurate open-source library for visual, visual–inertial, and multimap SLAM,” IEEE Transactions on Robotics, vol. 37, no. 6, pp. 1874–1890, 2021.
  * [5] T. Qin, P. Li, and S. Shen, “VINS-Mono: A robust and versatile monocular visual-inertial state estimator,” IEEE Transactions on Robotics, vol. 34, no. 4, pp. 1004–1020, 2018.
  * [6] T. Qin, J. Pan, S. Cao, and S. Shen, “A general optimization-based framework for local odometry estimation with multiple sensors,” 2019.
  * [7] G. Gallego, T. Delbrück, G. Orchard, C. Bartolozzi, B. Taba, A. Censi, S. Leutenegger, A. J. Davison, J. Conradt, K. Daniilidis, and D. Scaramuzza, “Event-based vision: A survey,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, no. 1, pp. 154–180, 2022.
  * [8] G. Gallego and D. Scaramuzza, “Accurate angular velocity estimation with an event camera,” IEEE Robotics and Automation Letters, vol. 2, no. 2, pp. 632–639, 2017.
  * [9] E. Mueggler, G. Gallego, and D. Scaramuzza, “Continuous-time trajectory estimation for event-based vision sensors,” in Robotics: Science and Systems, 2015.
  * [10] A. R. Vidal, H. Rebecq, T. Horstschaefer, and D. Scaramuzza, “Ultimate SLAM? combining events, images, and IMU for robust visual SLAM in HDR and high-speed scenarios,” IEEE Robotics and Automation Letters, vol. 3, no. 2, pp. 994–1001, 2018.
  * [11] H. Kim, A. Handa, R. Benosman, S.-H. Ieng, and A. Davison, “Simultaneous mosaicing and tracking with an event camera,” in Proceedings of the British Machine Vision Conference, BMVA Press, 2014.
  * [12] T. H. Henri Rebecq and D. Scaramuzza, “Real-time visual-inertial odometry for event cameras using keyframe-based nonlinear optimization,” in Proceedings of the British Machine Vision Conference (BMVC) (G. B. Tae-Kyun Kim, Stefanos Zafeiriou and K. Mikolajczyk, eds.), pp. 16.1–16.12, BMVA Press, September 2017.
  * [13] B. D. Lucas and T. Kanade, “An iterative image registration technique with an application to stereo vision,” in Proceedings of the 7th International Joint Conference on Artificial Intelligence, vol. 2 of IJCAI’81, p. 674–679, Morgan Kaufmann Publishers Inc., 1981.
  * [14] A. Z. Zhu, N. Atanasov, and K. Daniilidis, “Event-based visual inertial odometry,” in 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 5816–5824, 2017.
  * [15] J. Kannala and S. Brandt, “A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 28, no. 8, pp. 1335–1340, 2006\.
  * [16] E. Mueggler, H. Rebecq, G. Gallego, T. Delbruck, and D. Scaramuzza, “The event-camera dataset and simulator: Event-based data for pose estimation, visual odometry, and SLAM,” The International Journal of Robotics Research, vol. 36, no. 2, p. 142–149, 2017.
  * [17] P. Wyatt and H. Nakai, “Fast feature extraction using approximations to derivatives with summed-area images,” in Proceedings of the 7th Asian Conference on Computer Vision - Volume Part I, ACCV’06, (Berlin, Heidelberg), p. 776–786, Springer-Verlag, 2006.
  * [18] G. Gallego, H. Rebecq, and D. Scaramuzza, “A unifying contrast maximization framework for event cameras, with applications to motion, depth, and optical flow estimation,” 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, Jun 2018.
  * [19] C. Forster, L. Carlone, F. Dellaert, and D. Scaramuzza, “On-manifold preintegration for real-time visual–inertial odometry,” IEEE Transactions on Robotics, vol. 33, no. 1, p. 1–21, 2017.
  * [20] A. Z. Zhu, D. Thakur, T. Ozaslan, B. Pfrommer, V. Kumar, and K. Daniilidis, “The multivehicle stereo event camera dataset: An event camera dataset for 3D perception,” IEEE Robotics and Automation Letters, vol. 3, no. 3, p. 2032–2039, 2018.
  * [21] J. Sturm, N. Engelhard, F. Endres, W. Burgard, and D. Cremers, “A benchmark for the evaluation of RGB-D SLAM Systems,” in Proc. of the International Conference on Intelligent Robot Systems (IROS), Oct. 2012.
  * [22] R. H. A. Zisserman, Multiple View Geometry in Computer Vision, 2nd Edition.  2 ed., 2004.
  * [23] R. Mur-Artal and J. D. Tardós, “ORB-SLAM2: an open-source SLAM system for monocular, stereo and RGB-D cameras,” IEEE Transactions on Robotics, vol. 33, no. 5, pp. 1255–1262, 2017.
  * [24] A. Geiger, P. Lenz, and R. Urtasun, “Are we ready for autonomous driving? the KITTI vision benchmark suite,” in Conference on Computer Vision and Pattern Recognition (CVPR), 2012.


# Graph Representation Learning for Energy Demand Data:
Application to Joint Energy System Planning under Emissions Constraints

Aron Brenner 1, Rahman Khorramfar2, Dharik Mallapragada3, Saurabh Amin4

###### Abstract

A rapid transformation of current electric power and natural gas (NG)
infrastructure is imperative to meet the mid-century goal of CO2 emissions
reduction requires. This necessitates a long-term planning of the joint power-
NG system under representative demand and supply patterns, operational
constraints, and policy considerations. Our work is motivated by the
computational and practical challenges associated with solving the generation
and transmission expansion problem (GTEP) for joint planning of power-NG
systems. Specifically, we focus on efficiently extracting a set of
representative days from power and NG data in respective networks and using
this set to reduce the computational burden required to solve the GTEP. We
propose a Graph Autoencoder for Multiple time resolution Energy Systems
(GAMES) to capture the spatio-temporal demand patterns in interdependent
networks and account for differences in the temporal resolution of available
data. The resulting embeddings are used in a clustering algorithm to select
representative days. We evaluate the effectiveness of our approach in solving
a GTEP formulation calibrated for the joint power-NG system in New England.
This formulation accounts for the physical interdependencies between power and
NG systems, including the joint emissions constraint. Our results show that
the set of representative days obtained from GAMES not only allows us to
tractably solve the GTEP formulation, but also achieves a lower cost of
implementing the joint planning decisions.

## Introduction

One of the most significant societal challenges that we currently face is to
transition to a reliable, low-carbon, and sustainable energy system as soon as
possible, and to meet the mid-century goal of limiting global warming below
2∘C (UN-FCCC 2015; Gielen et al. 2019). This requires a significant use of
renewable energy resources and well-planned integration of various energy
vectors, including emerging clean energy sources such as hydrogen and other
renewable energy sources. Our work is motivated by the enormous potential of
machine learning (ML) models in promoting sustainable energy systems. In
particular, we focus on ML modeling for extracting a set of representative
days from heterogeneous demand data associated with real-world electric power
and natural gas (NG) systems, and using this set for joint power-NG network
planning under emissions constraints. In doing so, we leverage ML-extracted
representative days to tractably solve an optimization problem that determines
a capacity and network expansion plan for regional-scale energy systems such
as that of New England.

Broadly speaking, our work addresses several practical and computational
challenges associated with capacity expansion models (CEMs) for
decarbonization of interdependent power-NG infrastructures. Classical examples
of such models include the generation expansion problem (GEP) and generation
and transmission expansion problem (GTEP), both of which are well-studied in
the context of power systems (Li et al. 2022a; He et al. 2018). Our
optimization model is a GTEP that determines the optimal location and timing
of generation units, transmission lines, and pipelines to meet future energy
demands under a range of operational and policy constraints such as joint
emission constraints. In our work, we extend the model to include two main
interdependencies between power and NG systems. The first interdependency
captures the increasing role of gas-fired power plants in the generation mix
of electricity production ((2022) EIA; He et al. 2018). The second
interdependency reflects the joint emission of CO2 in both systems.

The key _computational challenge_ in solving the GTEP arises from the fact
that it links long-term investment decisions (e.g. capacity and network
expansion) to short-term operational ones (e.g. unit commitment, power
production, and energy storage). The former decisions have a planning horizon
of 10-30 years with yearly granularity, while the latter usually require
hourly or sub-hourly resolution. Under reasonable assumptions, we can express
the GTEP as a large-scale mixed-integer linear program (MILP), but current
literature has limited success in tractably solving these problems to an
adequate level of spatial and temporal resolution. In our case, the
computational difficulty in solving the GTEP increases further because we
model both power and NG networks. Thus, taking into account (projected) demand
information on a day-to-day basis becomes prohibitively expensive from a
computational viewpoint. In the classical GTEP problems for power systems, the
computational challenge is addressed by aggregating power system nodes (buses)
within a geographical neighborhood (power zone) to a single node (Li et al.
2022a) and by solving the GTEP for a set of representative days (Hoffmann et
al. 2020). Crucially, the set of representative days needs to capture demand
and supply patterns. To the best of our knowledge, the notion of
representative days has not been clearly defined and developed in the context
of joint power-NG planning problem – this is where we leverage our graph
representation learning approach.

Figure 1: Adjacent nodes in the power network demonstrate similar variations
in demand over the course of the day. These spatial dependencies are modeled
explicitly by graph convolutional layers in the GAMES architecture.

Our work also addresses the _practical issues_ arising from coarse data
availability from the NG network. Firstly, we do not have access to the
detailed connectivity and transmission information in the NG network while
this information is readily available for the power network. Secondly, power
systems typically collect demand and generation data at a fine temporal
resolution (hourly or less), but this data is usually not publicly accessible
for NG systems. These issues thus require us to (a) formulate network
constraints based on loosely specified information on power and NG node
connectivity and (b) develop an approach to leverage demand and supply data
from the power system with demand data of NG system despite their different
temporal resolutions.

We address the aforementioned challenges by developing a _graph representation
learning approach_ that captures the physical interdependencies between power
and NG networks, and also handles the different granularity of data at each
network. We consider demand data for both systems, and consider capacity
factor (CF) data for solar and wind plants to reflect the supply pattern in
the renewable-dominated future grid. We utilize graph convolutions to capture
the network interactions both within and across power and NG networks, and
adopt an autoencoder architecture with tuneable reconstruction losses for the
respective demand and CF data. We demonstrate that the resulting Graph
Autoencoder for Multiple time resolution Energy Systems (GAMES) model is
ideally suited to handle embedding the spatio-temporal patterns in power and
NG demand as well as wind and solar CF data into a lower-dimensional
representation, which can be readily clustered to extract the set of
representative days. Furthermore, our approach to computing the set of
representative days can also enable an accurate estimation of the trade-off
between costs (both investment and operational) and joint emissions from power
and NG systems.111We believe this capability can have a significant societal
impact by lowering the barriers to investment in renewable energy resources
and alleviating reliability concerns in a low-carbon energy system.

Previous studies for selecting representative days propose variants of k-means
(Mallapragada et al. 2018; Li et al. 2022b; Teichgraeber and Brandt 2019;
Barbar and Mallapragada 2022), k-medoids (Scott et al. 2019; Teichgraeber and
Brandt 2019), and hierarchical clustering (Liu, Sioshansi, and Conejo 2017;
Teichgraeber and Brandt 2019). The distance matrices used in clustering
algorithms for most previous works are constructed based on a set of time
series inputs such as load data and variable renewable energies (VRE) capacity
factors (Li et al. 2022a; Hoffmann et al. 2020). Notably, these approaches
neither account for demand data with multiple time resolutions nor account for
network interdependencies. Hence, they cannot be readily extended to address
the task of extracting representative days for joint power-NG systems – an
aspect that is crucial for realism and tractability in joint planning
optimization models for decarbonizing these systems. We believe that our GAMES
model addresses these challenges and provides a promising path to better
extract representative days in interdependent power and NG systems.

## Graph Convolutional Autoencoder Approach

In this section, we describe the Graph Autoencoder for Multiple time
resolution Energy Systems (GAMES) model, a simple graph autoencoder with
linear graph convolutions. We argue that this architecture efficiently
captures spatio-temporal demand patterns in power and NG systems.

### Autoencoders

To begin with, we note that direct use of clustering algorithms to identify
representative days for any large-scale energy system is prone to the “curse
of dimensionality” due to the high dimensionality of time series data. In such
settings, it is desirable to first extract low-dimensional and denoised
representations of the data before clustering (Parsons, Haque, and Liu 2004).
To identify a set of representative days, we choose to utilize a state-of-the-
art autoencoder architecture for learning low-dimensional embeddings for
power-NG systems (that have different time resolutions) prior to clustering.

Given a high-dimensional input such as a time series of graph signals,
$X\in\mathbb{R}^{p}$, an autencoder can be trained to jointly learn an
encoder, $g:\mathbb{R}^{p}\to\mathbb{R}^{k}$, and a decoder,
$f:\mathbb{R}^{k}\to\mathbb{R}^{p}$ that minimize the reconstruction loss
function $\|X-\hat{X}\|_{2}^{2}$, where $\hat{X}=f(g(X))$ is the reconstructed
signal. Here, $k\ll p$ denotes the dimension of the learned latent space.

Variable | Interpretation | Granularity | Nodes
---|---|---|---
$X_{\mathrm{E}}$ | Electricity | Hourly | 188
$X_{\mathrm{W}}$ | Wind | Hourly | 188
$X_{\mathrm{S}}$ | Solar | Hourly | 188
$X_{\mathrm{G}}$ | Natural Gas | Daily | 18
Table 1: Notation for input variables.

We denote by $X_{\mathrm{E}}\in\mathbb{R}^{d\times n_{\mathrm{E}}\times
t_{\mathrm{E}}}$ the data tensor of electricity demands over all days $d$,
nodes $n_{\mathrm{E}}$, and times $t_{\mathrm{E}}$. Similarly, we denote the
natural gas data tensor by $X_{\mathrm{G}}\in\mathbb{R}^{d\times
n_{\mathrm{G}}\times t_{\mathrm{G}}}$, the wind capacity factor tensor by
$X_{\mathrm{W}}\in\mathbb{R}^{d\times n_{\mathrm{W}}\times t_{\mathrm{W}}}$,
and the solar capacity factor data tensor by
$X_{\mathrm{S}}\in\mathbb{R}^{d\times n_{\mathrm{S}}\times t_{\mathrm{S}}}$
(see Table 1). Because the GTEP considers different associated costs for
investment and operational decisions related to power, NG, wind, and solar, we
introduce hyperparameters $\alpha_{G},\alpha_{W},\alpha_{S}$ in the
autoencoder objective function to tune the trade-off between the multiple
reconstruction losses. This parameter reflects the contribution of each system
towards the total cost. For example, if the NG system cost is twice the power
system cost, then higher values of $\alpha_{G}$ ensure that the reconstruction
cost is penalized more when deviating from the data of the NG system. This
gives us the following loss function:

$\displaystyle\sum_{i=1}^{d}\Big{(}\frac{1}{dn_{\mathrm{E}}t_{\mathrm{E}}}\|X_{\mathrm{E}}^{(i)}-\hat{X}_{\mathrm{E}}^{(i)}\|_{F}^{2}+\frac{\alpha_{G}}{dn_{\mathrm{G}}t_{\mathrm{G}}}\|X_{\mathrm{G}}^{(i)}-\hat{X}_{\mathrm{G}}^{(i)}\|_{F}^{2}$
$\displaystyle+\frac{\alpha_{W}}{dn_{\mathrm{W}}t_{\mathrm{W}}}\|X_{\mathrm{W}}^{(i)}-\hat{X}_{\mathrm{W}}^{(i)}\|_{F}^{2}+\frac{\alpha_{S}}{dn_{\mathrm{S}}t_{\mathrm{S}}}\|X_{\mathrm{S}}^{(i)}-\hat{X}_{\mathrm{S}}^{(i)}\|_{F}^{2}\Big{)},$

where $\|\cdot\|_{F}$ denotes the Frobenius norm.

In our case study, we set $\alpha_{G}=2,\alpha_{S}=0.5,\alpha_{W}=0.5$.
However, we note that it is possible to choose the hyperparameters by
evaluating the downstream GTEP objective for different values. Specifically,
this can be performed using a grid search in which the quality of a
combination of hyperparameters $\\{\alpha_{G},\alpha_{W},\alpha_{S}\\}$ is
measured by GTEP objective costs given by solving the optimization model
rather than the GAMES validation loss directly.

### Graph Representation Learning

Next, we provide a brief introduction to modeling with graph convolutional
networks (GCNs).

#### Preliminaries

We encode the network topology with the binary adjacency222Ideally, one should
construct an affinity matrix $A$ with a Gaussian kernel such that
$A_{ij}=\exp{\left(-\frac{\mathrm{dist}(i,j)^{2}}{\sigma^{2}}\right)}$ for all
edges $(i,j)$, where $\mathrm{dist}(i,j)$ denotes the distance of edge $(i,j)$
and $\sigma$ denotes the standard deviation of distances in the network
(Shuman et al. 2012). Since we do not have access to edge distance data in our
case study, we proceed with the binary adjacency matrix. matrix $A$, which we
construct such that

$\displaystyle A_{ij}=\begin{cases}0&(i,j)\notin\mathcal{E}\\\
1&(i,j)\in\mathcal{E}.\end{cases}$

We also construct the diagonal degree matrix $D$ such that
$D_{ii}=\sum_{j}A_{ij}$.

#### Graph Convolutions

Our graph autoencoder approach follows (Kipf and Welling 2017) in utilizing
Chebyshev convolutional filters, which approximate spectral convolutions to
learn node embeddings as weighted local averages of embeddings of adjacent
nodes. This is ideal for learning low-dimensional embeddings of energy
networks as neighborhoods of nodes typically exhibit similar energy demands
patterns and can thus be represented jointly. Chebyshev filters operate on the
“renormalized” graph Laplacian
$\tilde{L}=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}$, where
$\tilde{D}=I+D$ and $\tilde{A}=I+A$, and perform a form of Laplacian smoothing
(Li, Han, and Wu 2018; Taubin 1995). We initialize $H^{(0)}=X$ and apply
convolutional filters to learn subsequent node embeddings as follows:

$\displaystyle H^{(l+1)}=\sigma(\tilde{L}H^{(l)}\Theta^{(l)}),$

where $\Theta^{(l)}$ is a trainable weight matrix and $H^{(l)}$ is a matrix of
node embeddings in layer $l$. $\sigma(\cdot)$ is typically a nonlinear
activation function, such as $\mathrm{ReLU}$ or $\mathrm{tanh}$.

In each layer, GCNs aggregate features from the immediate neighborhood of each
node. Deep GCNs stack multiple layers with nonlinear activations to learn node
embeddings as nonlinear functions of both local and global node features. In
contrast, (Salha, Hennequin, and Vazirgiannis 2019) propose a simpler graph
autoencoder model, which they demonstrate to have competitive performances
with multilayer GCNs on standard benchmark datasets despite being limited to
linear first-order interactions. Shallow neural architectures are also better
suited for settings where data is scarce. This is particularly significant in
modeling energy systems whose data may only be available for a few historical
years. Indeed, we find this simpler GCN approach to perform well for our case
study. We now introduce GAMES, an augmented version of the linear GCN
autoencoder for energy systems with multiple time resolutions.

### GAMES

Figure 2: The GAMES Architecture. The electric power, wind CF, solar CF, and
NG time series are combined into the block matrix $X$ with $n$ rows and
$t_{\mathrm{E}}+t_{\mathrm{W}}+t_{\mathrm{S}}+t_{\mathrm{G}}$ channels. A
single linear graph convolutional layer constructs matrix $Z$ by embedding
each row of $X$ into $k$ dimensions. Another graph convolutional layer scales
each row of $Z$ back to
$t_{\mathrm{E}}+t_{\mathrm{W}}+t_{\mathrm{S}}+t_{\mathrm{G}}$ dimensions,
which are then separated and fed through fully connected layers to reconstruct
the two time series. After the model is trained, the embeddings are extracted
by feeding the daily time series inputs through the encoder, at which point
clustering is applied.

Our proposed GAMES architecture is designed as follows and illustrated in Fig.
2.

#### Encoder

Consider the power, wind CF, solar CF, and NG time series corresponding to day
$i$, $X_{\mathrm{E}}^{(i)}$, $X_{\mathrm{W}}^{(i)}$, $X_{\mathrm{S}}^{(i)}$,
$X_{\mathrm{G}}^{(i)}$. We begin by constructing the data matrix $X^{(i)}$ as

$\displaystyle
X^{(i)}=\begin{pmatrix}X_{\mathrm{E}}^{(i)}&X_{\mathrm{W}}^{(i)}&X_{\mathrm{W}}^{(i)}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{0}&X_{\mathrm{G}}^{(i)}\end{pmatrix}.$

Note that $X^{(i)}\in\mathbb{R}^{n\times t}$, where $n\coloneqq
n_{\mathrm{E}}+n_{\mathrm{G}}$ and $t\coloneqq
t_{\mathrm{E}}+t_{\mathrm{W}}+t_{\mathrm{W}}+t_{\mathrm{G}}$. This is because
capacity factor data exists for all nodes in the power network and utilizes
the same network topology. $X^{(i)}$ is then passed through a single
convolutional layer to produce the low-dimensional embedding
$Z^{(i)}\in\mathbb{R}^{n\times k}$. The hyperparameter $k$ defines the
bottleneck of the autoencoder architecture (i.e. the dimension of each node
embedding) and consequently the tradeoff between compression and
reconstruction loss. In our case study, we find $k=3$ to show a sufficient
performance for our application of identifying representative days.

#### Decoder

$Z^{(i)}$ is passed through a convolutional layer to produce the embedding
$H^{(i)}\in\mathbb{R}^{n\times t}$. This reconstructed matrix is then split
along the second dimension into two blocks:
$H_{\mathrm{E,W,S}}^{(i)}\in\mathbb{R}^{(n_{\mathrm{E}}+n_{\mathrm{W}}+n_{\mathrm{S}})\times
t}$ and $H_{\mathrm{G}}^{(i)}\in\mathbb{R}^{n_{\mathrm{G}}\times t}$. Each
block is then passed to a separate series of fully connected layers with
$\mathrm{tanh}$ activations that map the node embeddings in
$H_{\mathrm{E,W,S}}^{(i)}$ and $H_{\mathrm{G}}^{(i)}$ respectively to the
reconstructions $\hat{X}_{\mathrm{E,W,S}}^{(i)}$ and
$\hat{X}_{\mathrm{G}}^{(i)}$. Finally, the tensor
$\hat{X}_{\mathrm{E,W,S}}^{(i)}$ is split into the reconstructions
$\hat{X}_{\mathrm{E}}^{(i)},\hat{X}_{\mathrm{W}}^{(i)},\hat{X}_{\mathrm{S}}^{(i)}$.

### Clustering

After the model is trained, the power and NG time series from each day $i$,
i.e.
$(X_{\mathrm{E}}^{(1)},X_{\mathrm{W}}^{(1)},X_{\mathrm{S}}^{(1)},X_{\mathrm{G}}^{(1)}),\dots,(X_{\mathrm{E}}^{(N)},X_{\mathrm{W}}^{(N)},X_{\mathrm{S}}^{(N)},X_{\mathrm{G}}^{(N)})$,
is passed through the encoder to generate the embedding matrices
$Z^{(1)},\dots,Z^{(N)}$. Then, k-medoids clustering is applied to select a set
of $K$ cluster medians, denoted by $\mathcal{S}\subset\\{1,\dots,N\\}$, and
assign each day $i$ to a corresponding cluster $j\in\mathcal{S}$. We denote
the set of days assigned to the cluster defined by day $j$ as
$\mathcal{C}_{j}$. Given the number of clusters $K$, the k-medoids algorithm
aims to minimize the objective function

$\displaystyle\min\sum_{j\in\mathcal{S}}\sum_{i\in\mathcal{C}_{j}}\|Z^{(i)}-Z^{(j)}\|_{F}^{2}$
(1)

(Hastie, Tibshirani, and Friedman 2001). Note that every day in the dataset
must be assigned to exactly one cluster. Semantically, (1) can be understood
as aiming to ensure that the set of representative days $\mathcal{S}$
proportionately partitions the full set of days in the dataset by minimizing
squared Euclidean distances in the latent space as constructed by the
autoencoder.

## Capacity Expansion Model

The result of the clustering algorithm is used to solve the CEM for joint
power and NG planning, which is formulated as a GTEP. The problem determines
the minimum investment cost and operational decisions for the year 2050 under
various investment, operational, and policy constraints. The investment
decisions include establishing new power plants, transmission lines, and
pipelines as well as decommissioning existing plants. The operational
constraints include minimum production, ramping, energy balance, transmission,
and storage. We consider emission limits and minimum share of VREs as policy
constraints. Importantly, in our formulation, the emissions constraint limits
CO2 emissions incurred by the consumption of NG in _both_ networks.

We introduce our model with simplified notation in this section and provide a
detailed formulation in the supplementary material (SI 2022). Let
$\boldsymbol{z}^{e}=(\boldsymbol{x}^{e},\boldsymbol{y}^{e},\boldsymbol{p})$
represent the set of variables for the power system. The integer variable
$\boldsymbol{x}^{e}$ is the variable establishing plants, decommissioning
plants, and establishing new transmission lines. The continuous variable
$\boldsymbol{p}$ captures the power generation in NG-fired plants while
$\boldsymbol{y}^{e}$ is a continuous variable that captures all the remaining
variables including power generation from non NG-fired plants and power flow
between nodes, storage, and load shedding variables. We use
$\boldsymbol{z}^{g}=(\boldsymbol{x}^{g},\boldsymbol{y}^{g},\boldsymbol{f})$ to
denote the set of variables associated with the NG system. The mixed-integer
variable $\boldsymbol{x}^{g}$ is the set of all investment, storage, and load
shedding decisions. The continuous variable $\boldsymbol{y}^{g}$ represents
the intra-network flow, i.e. the flow between NG nodes or the flow between NG
nodes and NG storage facilities. The flow between NG and electricity systems
is denoted by $\boldsymbol{f}$. We formulate the joint power-NG system as
follows:

$\displaystyle\min\ $
$\displaystyle(\boldsymbol{c}^{e}_{1}\boldsymbol{x}^{e}+\boldsymbol{c}^{e}_{2}\boldsymbol{y}^{e}+\boldsymbol{c}^{e}_{3}\boldsymbol{p})+(\boldsymbol{c}^{g}_{1}\boldsymbol{x}^{g}+\boldsymbol{c}^{g}_{2}\boldsymbol{y}^{g}+\boldsymbol{c}^{g}_{3}\boldsymbol{f})$
(2a) s.t.
$\displaystyle\boldsymbol{A}^{e}\boldsymbol{x}^{e}+\boldsymbol{B}^{e}\boldsymbol{y}^{e}+\boldsymbol{D}^{e}\boldsymbol{p}\leq\boldsymbol{b}^{e}_{1}$
(2b)
$\displaystyle\boldsymbol{H}^{e}\boldsymbol{y}^{e}\geq\boldsymbol{b}^{e}_{2}$
(2c)
$\displaystyle\boldsymbol{A}^{g}\boldsymbol{x}^{g}+\boldsymbol{B}^{g}\boldsymbol{y}^{g}+\boldsymbol{D}^{g}\boldsymbol{f}\leq\boldsymbol{b}^{g}_{1}$
(2d) $\displaystyle\boldsymbol{f}=\boldsymbol{E}_{1}\boldsymbol{p}$ (2e)
$\displaystyle\boldsymbol{G}_{2}\boldsymbol{y}^{g}+\boldsymbol{E}_{2}\boldsymbol{p}\leq\eta$
(2f)
$\displaystyle\boldsymbol{x}^{e}\in\mathbb{Z}^{+},\boldsymbol{y}^{e},\boldsymbol{x}^{g}\in\mathbb{Z}^{+}\times\mathbb{R}^{+},\boldsymbol{p},\boldsymbol{y}^{g},\boldsymbol{f}\in\mathbb{R}^{+}$
(2g)

The objective function (2a) minimizes the investment and operational costs for
the power system (first term) and NG system (second term). The constraint (2b)
represents all investment, commitment, and operational constraints for the
power system including the production limit, ramping, storage, and energy
balance constraints. The constraint (2c) enforces policy considerations such
as the minimum requirement for renewable portfolio standard (RPS). The NG
constraints are reflected in constraint (2d), which includes technological and
operational constraints such as the supply limit at each node, flow between NG
nodes, and storage.

The coupling constraint (2e) ensures that NG-fired plants operate based on the
gas flow they receive from the NG network. The second coupling constraint (2f)
is the decarbonization constraint that limits emissions resulting from NG
consumption to serve both electricity (via NG power plants) and non-power
related NG loads to $\eta$. The coefficient matrices $\boldsymbol{E}_{1}$,
$\boldsymbol{G}_{2}$, and $\boldsymbol{E}_{2}$ represent the heat rate,
emission factors for NG usage, and emission factor for NG-fired plants,
respectively. Indeed, emissions from coal-fired plants is a major driver for
decarbonization efforts and NG remains as primary fuel for which emissions
need to be regulated. Therefore, given the declining role of coal in the US
energy system, the constraint (2f) reflects a futuristic setting where such
plants are already decommissioned.

## Input Data

Using publicly available data, we consider the New England region and
construct its corresponding power and NG network. We then calibrate the
resulting networks using historical data. The power network consists of 188
nodes with 338 existing and candidate transmission lines. The NG network
consists of 18 NG nodes and 7 storage nodes. We assume that each NG node is
connected to two other storage nodes. We also assume that each power node is
connected to three of its closest NG nodes. The Supplementary Information
provides the details of the input data for the joint power-NG planning model
(SI 2022).

## Computational Experiments

### GAMES Performance

We train GAMES on a dataset of 292 days using the Adam optimizer with a
learning rate of $0.001$. We use the full batch of 292 data points for each
update step and perform early stopping to end training when the validation
loss no longer decreases. We report the validation reconstruction loss on a
set of 73 days for various node embedding dimensions $k$ in Table 2.

Embed. Dim. | $k=1$ | $k=2$ | $k=3$ | $k=4$
---|---|---|---|---
MSE Loss | 0.727 | 0.398 | 0.244 | 0.160
Table 2: The reconstruction loss shows diminishing returns for $k>3$ node
embedding dimensions.

We observe slightly diminishing returns for the validation reconstruction loss
for $k>3$. Consequently, we proceed with our representative day selection
using embeddings generated by the model corresponding to $k=3$.

### Representative Days Comparison

#### Setup

(a) GAMES vs. raw data clustering comparison under an 80% carbon reduction
goal.

(b) GAMES vs. raw data clustering comparison under a 95% carbon reduction
goal.

Figure 3: Various costs and power emission for different number of
representative days under different decarbonization goals.

We use the k-medoids clustering algorithm to obtain different sets of
representative days. We apply the clustering algorithm to both raw data and
the embeddings obtained from the GAMES model to compare the results of the
proposed model. Accordingly, two different sets are obtained for each number
of representative days. The optimization model over the full power network is
prohibitively challenging even for a very small number of days. Therefore, we
aggregate all buses in each state of the New England region to obtain a 6-node
power network. This aggregation allows us to run the formulation for up to 42
representative days.

We obtain a feasible solution in two steps for each set of representative
days: (1) The optimization model is aggregated to the set of representative
days for tractability and then solved. (2) Next, we consider the full planning
horizon (the entire year of 2050) and set the integer decision variables (i.e.
investment decisions) to the values determined in the first step. We note that
the investment decision variables in our formulation are (a) the only integer-
valued decision variables and (b) independent of planning periods. Therefore,
fixing them reduces the remaining operational problem to a linear program
(LP), which can be solved considerably faster. The resulting solution from the
second step is a feasible solution to the full-year problem, with which we can
analyze resulting costs and decisions.

Table 3: Average percentage change when using GAMES approach for for various costs and power emissions. Reduction Goal | Total | Power | NG | Inv-FOM (Power) | Shedding cost (Power) | Emission from Power Sys
---|---|---|---|---|---|---
80$\%$ | -5.14 | -7.03 | 0.24 | -4.64 | -24.13 | -9.87
95$\%$ | -7.27 | -10.50 | 1.5 | -8.51 | -27.80 | -3.31

In our computational experiments, we consider two decarbonization goals of
80$\%$ and 95$\%$ where the former is the projected target for New England
states (Weiss and Hagerty 2019), and the latter aims reflects a radical
decarbonization goal. Figures 3(a) and 3(b) show the results under 80$\%$ and
95$\%$ emission reduction goals respectively. Both figures evaluate the
following quantities for the clusters obtained from GAMES and raw data: i)
“Total Cost” which is the objective function of model 2; ii) “Power System
Cost” which is the first term in the objective function (2a); iii) “NG System
Cost” which is the second term in the objective function (2a); iv) “Investment
and FOM for Generation and Storage (Pow. Sysm)” (investment-FOM) which is part
of the power system cost and captures the capital investment and fixed
operating and maintenance (FOM) costs of installing new power plants and
storage systems; v) “Power System Load Shedding Cost” which is part of the
power system cost and reflects the cost of unsatisfied electricity demand; and
v) “Emission from Power System” which is the tonnage of emission as a result
of operating NG-fired power plants in the power system. We use “GAMES” to
denote the feasible solution for the set of days obtained by GAMES. We do not
report the wall-clock times, but all instances are solved under 5 hours. As
expected, run-times vary significantly depending on the number of
representative days utilized; instances with 2 representative days typically
run in fewer than 350 seconds, whereas 30-day instances may need to run for
2800 seconds. All instances are implemented in Python using Gurobi 9.5 and are
run on the MIT Supercloud system with an Intel Xeon Platinum 8260 processor
containing up to 96 cores and 192 GB of RAM (Reuther et al. 2018).

#### Results

Table 3 presents the percentage change in various quantities yielded by the
GAMES representative days solution as compared to the solution using
representative days selected from clustering the raw data. The cost
comparisons are also plotted in Figures 3(a) and 3(b). We observe, on average,
a 5.14$\%$ and 7.27$\%$ improvement (decrease) in the total cost when using
GAMES under 80$\%$ and 95$\%$ decarbonization goals, respectively. This
improvement may be attributed to GAMES’ ability to model dependencies between
power and NG system data. Under more stringent decarbonization targets, the
share of VRE increases and the role of dispatchable power plants, such as NG-
fired plants, diminishes. As a result, modeling the influence of capacity
factors and their interactions with power and gas demands becomes more
essential. This phenomenon may underlie our observation for the 22-day
instance in which, while both approaches provide similar results under the
80$\%$ decarbonization goal, GAMES significantly outperforms the raw data
clustering as measured by total cost for the higher decarbonization goal. As
shown in Figures 3(a) and 3(b), the total cost from GAMES outperforms or
matches the performance of the raw data clustering in all instances (except
the 30-day instance under an 80$\%$ reduction goal). Interestingly, this
disparity in performance is most drastic when 15 or fewer representative days
are utilized under both decarbonization goals. This is worth noting as the
optimization model instantiated on the full network topology (i.e. without
aggregating nodes by state) is only tractable over a small set of
representative days (i.e. after applying a very coarse temporal aggregation).
It is especially important when the a model-year model only affords to
consider a handful of representative days for each year.

The power system cost largely drives variation in the total cost under both
decarbonization goals – the total cost is lower for all solutions with a lower
power system cost. Note that the difference in performance is more pronounced
in the power system cost compared to the total cost as indicated by the
7.03$\%$ and 10.50$\%$ power system cost improvement for GAMES under the
80$\%$ and 95$\%$ decarbonization goals. In Figure 3(a), this trend aligns
with load shedding costs except for the 14-day instance. However, as the
24.13$\%$ decrease shows, the GAMES approach results in significantly lower
load shedding on average. The 27.80$\%$ improvement in the load shedding cost
for GAMES under the 95$\%$ goal is plotted in detail in Figure 3(b); GAMES
outperforms the raw data clustering for all instances. Moreover, the GAMES
approaches converges after 14 days with load shedding cost significantly lower
than those instances utilizing fewer representative days.

In both figures the trends of investment-FOM cost and power system cost are
the same, indicating that the power system cost is largely driven by
investment-FOM cost, and to a lesser extent, by load shedding cost. This is
expected as future energy systems will rely significantly on VREs such as
solar and wind power, which only incur investment and FOM costs. Another
interesting observation pertains to the quantity of emissions in the power
system caused by operating NG-fired plants. Emissions for the power system are
on average 9.87$\%$ and 3.31$\%$ lower for GAMES under the two decarbonization
goals. This indicates a greater share of VREs in the GAMES approach, and
correspondingly, a higher share of gas-fired plants in the raw data clustering
approach. This is an interesting observation that may have significant
implications for energy policy-making. In particular, it suggests that the
results from the raw data clustering approach may be misleading as they do not
sufficiently convey the radical changes required to transform the system from
the current gas-dominant generation portfolio to a renewable-dominant power
grid.

NG system cost is another essential component of the total costs. Although NG
costs are similar for GAMES and raw data clustering for each instance, the NG
cost increases with the number of representative days. A possible explanation
might be that neither GAMES nor raw data clustering aim to capture extreme
days with separate clusters. Therefore, days with loads similar to extreme
days are more likely to be selected as a cluster’s medoid as the number of
representative days increases, which inevitably raises the NG system cost.
This consideration is also consistent with the observed load shedding cost for
the power system, which is significantly higher for instances with fewer than
15 representative days, indicating that both approaches fail to account for
extreme days in cluster medoids.

## Conclusion

In this work, we propose GAMES, a graph convolutional autoencoder for modeling
energy demand in interdependent electric power and natural gas systems with
heterogeneous nodes and different time resolutions. GAMES is able to exploit
spatio-temporal demand patterns to learn efficient embeddings of
interdependent power and NG networks. We apply the k-medoids clustering
algorithm to these embeddings to identify a set of representative days with
which we are able to tractably solve an energy system infrastructure planning
problem calibrated for the joint power-NG system in New England. Our
computational results show that the proposed framework outperforms clustering
methods applied to the raw data and is effective in selecting a small number
of representative days to provide high-quality feasible solutions for the
optimization problem.

The current work can be extended in multiple directions. The immediate
extension of the GCN architecture is to explore alternative approaches to
graph representation learning such as Laplacian sharpening (Park et al. 2019).
The extraction and inclusion of extreme days, or low-frequency days with
unusually low or high demand is another potential next step which could
prevent high load shedding costs and better represents the NG system’s load
patterns.

## References

  * Barbar and Mallapragada (2022) Barbar, M.; and Mallapragada, D. S. 2022.  Representative period selection for power system planning using autoencoder-based dimensionality reduction.  _arXiv preprint arXiv:2204.13608_.
  * (2) (EIA), E. I. A. 2022.  EIA Website.  Website.  Accessed: 2022-2-18.
  * Gielen et al. (2019) Gielen, D.; Gorini, R.; Wagner, N.; Leme, R.; Gutierrez, L.; Prakash, G.; Asmelash, E.; Janeiro, L.; Gallina, G.; Vale, G.; et al. 2019.  Global energy transformation: a roadmap to 2050.
  * Hastie, Tibshirani, and Friedman (2001) Hastie, T.; Tibshirani, R.; and Friedman, J. 2001.  _The Elements of Statistical Learning_.  Springer Series in Statistics. New York, NY, USA: Springer New York Inc.
  * He et al. (2018) He, C.; Zhang, X.; Liu, T.; Wu, L.; and Shahidehpour, M. 2018.  Coordination of interdependent electricity grid and natural gas network—a review.  _Current Sustainable/Renewable Energy Reports_ , 5(1): 23–36.
  * Hoffmann et al. (2020) Hoffmann, M.; Kotzur, L.; Stolten, D.; and Robinius, M. 2020.  A review on time series aggregation methods for energy system models.  _Energies_ , 13(3): 641.
  * Kipf and Welling (2017) Kipf, T. N.; and Welling, M. 2017.  Semi-Supervised Classification with Graph Convolutional Networks.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net.
  * Li et al. (2022a) Li, C.; Conejo, A. J.; Liu, P.; Omell, B. P.; Siirola, J. D.; and Grossmann, I. E. 2022a.  Mixed-integer linear programming models and algorithms for generation and transmission expansion planning of power systems.  _European Journal of Operational Research_ , 297(3): 1071–1082.
  * Li et al. (2022b) Li, C.; Conejo, A. J.; Siirola, J. D.; and Grossmann, I. E. 2022b.  On representative day selection for capacity expansion planning of power systems under extreme operating conditions.  _International Journal of Electrical Power & Energy Systems_, 137: 107697.
  * Li, Han, and Wu (2018) Li, Q.; Han, Z.; and Wu, X.-M. 2018.  Deeper Insights into Graph Convolutional Networks for Semi-Supervised Learning.  In _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence and Thirtieth Innovative Applications of Artificial Intelligence Conference and Eighth AAAI Symposium on Educational Advances in Artificial Intelligence_ , AAAI’18/IAAI’18/EAAI’18. AAAI Press.  ISBN 978-1-57735-800-8.
  * Liu, Sioshansi, and Conejo (2017) Liu, Y.; Sioshansi, R.; and Conejo, A. J. 2017.  Hierarchical clustering to find representative operating periods for capacity-expansion modeling.  _IEEE Transactions on Power Systems_ , 33(3): 3029–3039.
  * Mallapragada et al. (2018) Mallapragada, D. S.; Papageorgiou, D. J.; Venkatesh, A.; Lara, C. L.; and Grossmann, I. E. 2018.  Impact of model resolution on scenario outcomes for electricity sector system expansion.  _Energy_ , 163: 1231–1244.
  * Park et al. (2019) Park, J.; Lee, M.; Chang, H. J.; Lee, K.; and Choi, J. Y. 2019.  Symmetric Graph Convolutional Autoencoder for Unsupervised Graph Representation Learning.  _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , 6518–6527.
  * Parsons, Haque, and Liu (2004) Parsons, L.; Haque, E.; and Liu, H. 2004.  Subspace Clustering for High Dimensional Data: A Review.  _SIGKDD Explor. Newsl._ , 6(1): 90–105.
  * Reuther et al. (2018) Reuther, A.; Kepner, J.; Byun, C.; Samsi, S.; Arcand, W.; Bestor, D.; Bergeron, B.; Gadepally, V.; Houle, M.; Hubbell, M.; et al. 2018.  Interactive supercomputing on 40,000 cores for machine learning and data analysis.  In _2018 IEEE High Performance extreme Computing Conference (HPEC)_ , 1–6. IEEE.
  * Salha, Hennequin, and Vazirgiannis (2019) Salha, G.; Hennequin, R.; and Vazirgiannis, M. 2019.  Keep it simple: Graph autoencoders without graph convolutional networks.  _arXiv preprint arXiv:1910.00942_.
  * Scott et al. (2019) Scott, I. J.; Carvalho, P. M.; Botterud, A.; and Silva, C. A. 2019.  Clustering representative days for power systems generation expansion planning: Capturing the effects of variable renewables and energy storage.  _Applied Energy_ , 253: 113603.
  * Shuman et al. (2012) Shuman, D.; Narang, S. K.; Frossard, P.; Ortega, A.; and Vandergheynst, P. 2012\.  The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains.  _IEEE Signal Processing Magazine_ , 30.
  * SI (2022) SI. 2022.  Supplementary material available at: https://shorturl.at/bkHOU.
  * Taubin (1995) Taubin, G. 1995.  A Signal Processing Approach to Fair Surface Design.  In _Proceedings of the 22nd Annual Conference on Computer Graphics and Interactive Techniques_ , SIGGRAPH ’95, 351–358. New York, NY, USA: Association for Computing Machinery.  ISBN 0897917014.
  * Teichgraeber and Brandt (2019) Teichgraeber, H.; and Brandt, A. R. 2019.  Clustering methods to find representative periods for the optimization of energy systems: An initial framework and comparison.  _Applied energy_ , 239: 1283–1293.
  * UN-FCCC (2015) UN-FCCC. 2015.  Decision 1/CP. 21, Adoption of the Paris Agreement.  In _Report of the Conference of the Parties on Its Twenty-First Session, Held in Paris from_ , volume 30.
  * Weiss and Hagerty (2019) Weiss, J.; and Hagerty, J. M. 2019.  Achieving 80$\%$ GHG Reduction in New England by 2050.