# Hierarchical Graph Modeling for
Multi-Scale Optimization of Power Systems

David L. Cole1, Harsha Gangammanavar2, Victor M. Zavala13  1 Department of
Chemical and Biological Engineering, University of Wisconsin-Madison, Madison,
WI, United States
<EMAIL_ADDRESS>2 Operations Research and Engineering Management, Southern
Methodist University, Dallas, TX, United States
<EMAIL_ADDRESS>3 Mathematics and Computer Science Division, Argonne National
Laboratory, Lemont, IL, United States
<EMAIL_ADDRESS>
###### Abstract

Hierarchical optimization architectures are used in power systems to manage
disturbances and phenomena that arise at multiple spatial and temporal scales.
We present a graph modeling abstraction for representing such architectures
and an implementation in the Julia package Plasmo.jl. We apply this framework
to a tri-level hierarchical framework arising in wholesale market operations
that involves day-ahead unit commitment, short-term unit commitment, and
economic dispatch. We show that graph abstractions facilitate the
construction, visualization, and solution of these complex problems.

###### Index Terms:

Graph Theory, Hirearchical Optimization, Multiscale, Power Systems

##  I Introduction

Hierarchical optimization architectures are used in power systems (and many
other industrial systems) for managing operations, disturbances, and phenomena
that arise at multiple spatial and temporal scales. These architectures
involve multiple decision-making layers where decisions of higher layers
influence or inform lower layers (and vice versa); for example, market
operations often involve the solution of a unit commitment (UC) problem whose
solution informs an economic dispatch (ED) problem [1]. Hierarchical
decomposition is often necessary for enabling scalable implementation (e.g.,
solving a combined UC/ED problem in real-time might be impossible) and for
providing intuitive decomposition of functionalities (which can aid
explainability). Capturing the unique characteristics of optimization problems
in different hierarchical layers (i.e., space/time resolution, data,
variables, objectives, constraints) and their hierarchical coupling is
essential for enabling decision-making consistency across scales. This has
motivated research in models and solution approaches that aim to identify how
to best design hierarchical architectures to manage diverse types of features
(e.g., identify the number of layers, resolutions, and decisions made by each
layer). For example, Atakan and co-workers [2] presented a stochastic
optimization framework that consists of a tri-level hierarchy of market
operations (day-ahead UC, short-term UC, and ED) that aims to handle high
renewable penetration. The authors demonstrate that the hierarchical framework
provides significant operational improvement over competing architectures. Guo
and co-workers [3] used a hierarchical architecture for decentralizing ED of a
large power networks; this was a tri-level architecture where local, clustered
agents (lowest layer) inform leader agents (middle layer), which in turn
inform a coordinating agent (top layer). Kong and co-workers [4] proposed a
hierarchical architecture for a network of electric vehicle charging stations
connected to the grid; the formulation considers the placement of stations,
the allocation of resources, and the operation policy of the stations on three
separate hierarchical layers. They found that the framework provided improved
system performance and quality of service.

As power systems become increasingly complex (e.g., they include new assets
and face new disturbances), it will become necessary to have modeling and
solution tools that enable the seamless construction, evaluation, and
benchmarking of different hierarchical architectures. In this work, we propose
a graph-based modeling framework for representing hierarchical optimization
structures arising in power system operations. The use of graphs to model
structured optimization problems has been recently explored [5, 6, 7, 8, 9]. A
variety of tools for exploiting graph and graph-like structures are also
available in open-source packages such as Plasmo.jl (in Julia) [5] and Pyomo
(in Python) [10, 11]. In this work, we focus on the use of Plasmo.jl; this
package uses an OptiGraph abstraction, where nodes of the graph contain
optimization subproblems (with their own objective functions, data, variables,
and constraints) and where edges capture connectivity (constraints) across
subproblems. The OptiGraph abstraction is flexible in that nodes can contain
subproblems of different granularity; moreover, the abstraction enables the
creation of hierarchical structures (a node can be a graph itself). The graph
abstraction provides the ability to build complex structures in a modular
manner (e.g., node by node) and the ability to visualize, decompose, and
aggregate the overall problem graph. We provide a case study show how the
graph representation can be used for expressing and solving complex
hierarchical problems arising in power systems.

##  II Graph-Based Modeling Overview

Plasmo.jl is a Julia package that models general optimization problems as
hypergraphs. This package has been described in detail by Jalving and co-
workers [5], but here we provide a short overview of how this can be used for
representing hierarchical problems. Plasmo.jl is built on an abstraction
called OptiGraphs, which are graphs containing OptiNodes ($\mathcal{N}$) and
OptiEdges ($\mathcal{E}$). OptiNodes contain subproblems (with their own
variables, constraints, data, and objective functions), and OptiEdges are
linking constraints that capture connectivity between Optinodes. We denote an
OptiGraph as $\mathcal{G}(\mathcal{N},\mathcal{E})$, where
$\mathcal{N}(\mathcal{G})$ is the set of OptiNodes in $\mathcal{G}$ and
$\mathcal{E}(\mathcal{G})$ is the set of OptiEdges in $\mathcal{G}$. A
visualization of an OptiGraph containing three OptiNodes is shown in Figure 1.
The optimization model associated with an OptiGraph can be represented as

$\displaystyle\begin{split}\min_{\\{x_{n}\\}_{n\in\mathcal{N}(\mathcal{G})}}&\;\sum_{n\in\mathcal{N}(\mathcal{G})}f_{n}(x_{n})\\\
\textrm{s.t.}&\;x_{n}\in\mathcal{X}_{n},\quad n\in\mathcal{N}(\mathcal{G})\\\
&\;g_{e}(\\{x_{n}\\}_{n\in\mathcal{N}(e)})\geq 0,\quad
e\in\mathcal{E}(\mathcal{G})\end{split}$ (1)

where $\mathcal{N}(e)$ is the set of OptiNodes that support OptiEdge $e$. The
notion of nodes and edges is highly flexible in this abstraction; for
instance, in a power system context, a node can represent a spatial location,
time instance, a specific asset, or an entire network. Moreover, each node can
have its own independent features (e.g., data, objective functions,
constraints). The edges (containing the constraints $g_{e}$) can be used to
link nodes across time, space, or hierarchical layers.

OptiGraphs enable hierarchical representations and modular model building via
the use of subgraphs. Specifically, within Plasmo.jl, an OptiGraph can be
embedded in another OptiGraph as a node. For example, consider the OptiGraphs
$\mathcal{G}_{i}$ and $\mathcal{G}_{j}$, each with an independent set of
OptiNodes and OptiEdges. We can consider these OptiGraphs as low-level graphs
(also referred to as subgraphs) that can be used to build a higher-level
OptiGraph which we denote by
$\mathcal{G}(\\{\mathcal{G}_{i},\mathcal{G}_{j}\\},\mathcal{N}_{g},\mathcal{E}_{g})$.
The set of nodes $\mathcal{N}_{g}$ are contained on $\mathcal{G}$ and are
separate from $\mathcal{N}(\mathcal{G}_{i})$ and
$\mathcal{N}(\mathcal{G}_{j})$, such that
$\mathcal{N}_{g}=\mathcal{N}(\mathcal{G})/\\{\mathcal{N}(\mathcal{G}_{i})\cup\mathcal{N}(\mathcal{G}_{j})\\}$.
Similarly,
$\mathcal{E}_{g}=\mathcal{E}(\mathcal{G})/\\{\mathcal{E}(\mathcal{G}_{i})\cup\mathcal{E}(\mathcal{G}_{j})\\}$,
meaning $\mathcal{E}_{g}$ may connect nodes across $\mathcal{N}_{g}$,
$\mathcal{N}(\mathcal{G}_{i})$, and/or $\mathcal{N}(\mathcal{G}_{j}$). The
OptiGraph $\mathcal{G}$ may also be placed in another higher-level OptiGraph
$\mathcal{G}^{\prime}(\\{\mathcal{G}\\},\mathcal{N}_{g}^{\prime},\mathcal{E}_{g}^{\prime})$.
Any subgraph can also be collapsed into a single OptiNode (containing the
entire problem of the subgraph); this is useful for visualizing hierarchies.
For instance, the hierarchical setting allows us to capture how assets can be
aggregated in space (e.g., assets can be embedded at a network location), or
how multiple time points can be embedded in another time point (e.g., multiple
5-min time periods can be embedded in an hour). This feature is key for
representing hierarchical structures that span multiple scales.

Figure 1: OptiGraph abstraction used in Plasmo.jl. Figure 2: OptiGraph
abstraction of bi-level problem.

The general approach for representing hierarchical problems as graphs is
illustrated in Figure 2. This is a bi-level hierarchical problem; the top-
level OptiGraph is given by
$\mathcal{G}_{1}(\\{\mathcal{G}_{1,1},\mathcal{G}_{1,2}\\},\emptyset,\mathcal{E}_{1})$,
the lower-level OptiGraph by
$\mathcal{G}_{2}(\\{\mathcal{G}_{2,1},\mathcal{G}_{2,2},\mathcal{G}_{2,3},\mathcal{G}_{2,4}\\},\emptyset,\mathcal{E}_{2})$,
and the overall OptiGraph by
$\mathcal{G}_{0}(\\{\mathcal{G}_{1},\mathcal{G}_{2}\\},\emptyset,\mathcal{E}_{0})$,
where $\mathcal{E}_{0}$ are the constraints linking the solutions of the upper
and lower layers.

OptiGraphs enable flexible partitioning of hierarchical structures, and this
can be used to implement different solution approaches. For instance, it has
been recently shown that graph structures facilitate the development of
decomposition algorithms [5, 6, 9, 12, 13, 14]. To be specific, any OptiNode
or subgraph can be treated as an individual optimization problem; for example,
the OptiGraph presented in Figure 2 can be solved in at least three different
ways (each likely resulting in different solutions): i) $\mathcal{G}_{0}$
could be solved as a single monolithic problem; ii) $\mathcal{G}_{1}$ and
$\mathcal{G}_{2}$ can be solved sequentially with the solution of
$\mathcal{G}_{1}$ passed via $\mathcal{E}_{0}$; iii) $\mathcal{G}_{1,1}$ and
$\mathcal{G}_{1,2}$ can be solved sequentially with the solution of
$\mathcal{G}_{1,1}$ passed via $\mathcal{E}_{1}$. Then those solutions to
$\mathcal{G}_{1}$ can be passed via $\mathcal{E}_{0}$ to $\mathcal{G}_{2}$,
and $\mathcal{G}_{2,1}$, $\mathcal{G}_{2,2}$, $\mathcal{G}_{2,3}$, and
$\mathcal{G}_{2,4}$ can be solved sequentially, with solutions passed via
$\mathcal{E}_{2}$. We can thus see that the OptiGraph abstraction offers
significant flexibility in modeling and solving hierarchical problems.

##  III Case Study

###  III-A Problem Overview

We consider the tri-level problem proposed in [2] for capturing coupling in
market operations (see Figure 3). Each layer is composed of subproblems at
different timescales and these are linked to subproblems in other layers. The
top layer is a day-ahead unit commitment (DA-UC) problem that schedules a
subset of conventional (non-renewable) generators (denoted as
$\Gamma^{d}_{c}$). The DA-UC layer has a 1-hour resolution and a 24-hour
horizon, and an entire time horizon is partitioned into periods of 24 hours
(DA-UC is solved every 24 hours). The second layer includes a short-term unit
commitment (ST-UC) problem; this schedules a subset of conventional generators
(denoted as $\Gamma_{c}^{s}$), such that
$\Gamma_{c}^{d}\cap\Gamma_{c}^{s}=\emptyset$, while also incorporating the
commitment decisions of the DA-UC subproblems. The ST-UC layer has a 15-min
resolution with subproblems containing a 4-hour horizons and solved every 3
hours (there is overlap). The bottom layer is an hour-ahead economic dispatch
(HA-ED) layer which determines the generation levels for units committed in
the DA-UC and ST-UC layers. The HA-ED subproblems have a 15-minute resolution
and a 75-min time horizon, and are solved every 15 minutes (there is overlap).
Thus, for a given day, there are: 1 DA-UC subproblem, 8 ST-UC subproblems, and
96 HA-ED subproblems (12 for each ST-UC subproblem). We highlight that this
architecture is just one design (of many possible ones). In other words, one
could design diverse hierarchical architectures (e.g., experimenting with the
types of variables, resolutions, and time horizons that each layer uses).

The detailed model can be found in [2]; here, we provide a high-level
perspective to illustrate how complex models are embedded in the different
layers and how coupling arises between layers. We use the sets $\Gamma$ for
the set of all generators, $\Gamma_{r}$ for the set of renewable generators,
and define $\Gamma_{c}^{h}=\Gamma_{c}^{d}\cup\Gamma_{c}^{s}$. We also use $*$
and $**$ to define variables, sets, or functions corresponding to a layer;
here, the symbols $*$ or $**$ are exchanged for $d$, $s$, or $h$ to denote the
DA-UC, ST-UC, or HA-ED layers, respectively. As each subproblem considers
different sets of times, we use $\mathcal{T}^{*}_{i}$ for the set of times (in
hours) of the $i$th subproblem of the $*$ layer. We also define the sets
$\bar{\mathcal{T}}^{*}_{i}$ as the set of times without the first time point
of the subproblem. We define $\Delta^{*}$ as the time step for the subproblem
in hours($\Delta^{d}=1$, $\Delta^{s}=0.25$, and $\Delta^{h}=0.25$). Note that
we set $\Delta^{s}=\Delta^{h}$, and this has a small influence on the
formulation of the problem presented below. The decision variables are given
by:

$\displaystyle\begin{split}\boldsymbol{x}^{d}_{i}=&\;\big{(}x_{g,t},s_{g,t},z_{g,t}\big{)}_{\forall
g\in\Gamma^{d}_{c},t\in\mathcal{T}^{d}_{i}}\\\
\boldsymbol{y}^{d}_{i}=&\;\Big{(}(G^{+,d}_{g,t},G^{-,d}_{g,t})_{\forall{g\in\Gamma^{d}_{c}\cup\Gamma_{r}}},(F^{d}_{j,k,t})_{\forall(j,k)\in\mathcal{L}},\\\
&\;(D^{d}_{j,t},\theta^{d}_{j,t})_{\forall j\in\mathcal{B}}\Big{)}_{\forall
t\in\mathcal{T}^{d}_{i}}\\\
\color[rgb]{1,0,0}\boldsymbol{x}^{s}_{i}=&\;\color[rgb]{1,0,0}\big{(}x_{g,t},s_{g,t},z_{g,t}\big{)}_{\forall
g\in\Gamma^{s}_{c},t\in\mathcal{T}^{s}_{i}}\\\
\color[rgb]{1,0,0}\boldsymbol{y}^{s}_{i}=&\;\color[rgb]{1,0,0}\Big{(}(G^{+,s}_{g,t},G^{-,s}_{g,t})_{\forall{g\in\Gamma_{c}^{h}\cup\Gamma_{r}}},(F^{s}_{j,k,t})_{\forall(j,k)\in\mathcal{L}},\\\
&\;\color[rgb]{1,0,0}\qquad(D^{s}_{j,t},\theta^{s}_{j,t})_{\forall
j\in\mathcal{B}}\Big{)}_{\forall t\in\mathcal{T}^{s}_{i}}\\\
\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{i}=&\;\color[rgb]{0,0,1}\Big{(}(G^{+,h}_{g,t},G^{-,h}_{g,t})_{\forall{g\in\Gamma_{c}^{h}\cup\Gamma_{r}}},(F^{h}_{j,k,t})_{\forall(j,k)\in\mathcal{L}},\\\
&\;\color[rgb]{0,0,1}\qquad(D^{h}_{j,t},\theta^{h}_{j,t})_{\forall
j\in\mathcal{B}}\Big{)}_{\forall t\in\mathcal{T}^{s}_{i}}\end{split}$

Symbols $\boldsymbol{x}^{d}_{i}$ and $\boldsymbol{y}^{d}_{i}$ are binary and
continuous decision variables for the $i$th subproblem of DA-UC,
$\color[rgb]{1,0,0}\boldsymbol{x}^{s}_{i}$ and
$\color[rgb]{1,0,0}\boldsymbol{y}^{s}_{i}$ are decision variables for the
$i$th subproblem of ST-UC, and $\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{i}$ are
decision variables for the $i$th subproblem o HA-ED. Symbols $x_{g,t}$,
$s_{g,t}$, and $z_{g,t}$ are binary variables indicating for time $t$ whether
generator $g$ is on/off, was turned on, or was turned off. $G^{+,*}_{g,t}$ is
the power of generator $g$ consumed by the grid at time $t$ and
$G^{-,*}_{g,t}$ is the power overgenerated (for conventional generators) or
curtailed (for renewable generators) from $g$ at time $t$. $F^{*}_{j,k,t}$ is
the power flow of transmission line $(j,k)$ from bus $j$ to bus $k$ during
time $t$. $D^{*}_{j,t}$ is the amount of load shed at bus $j$ for time $t$,
and $\theta^{*}_{j,t}$ is the bus angle for bus $j$ at time $t$. We also use
black, red, and blue color to denote variables for DA-UC, ST-UC, and HA-ED,
respectively.

Figure 3: The tri-level hierarchical architecture of Atakan et al. [2]
(Reproduced with permission from Elsevier).

The objective functions in the different layers are comprised of a UC part,
$f_{u}$, and an ED part, $f_{e}$:

$\displaystyle
f_{u}(\boldsymbol{x}^{*}_{i}):=\sum_{t\in\mathcal{T}^{*}_{i}}\sum_{g\in\Gamma^{*}_{c}}\phi_{g}^{s}s_{g,t}+\phi_{g}^{f}x_{g,t}\Delta^{*}$
(2)
$\displaystyle\begin{split}f_{e}(\boldsymbol{y}^{*}_{i}):=&\;\sum_{t\in\mathcal{T}^{*}_{i}}\Big{(}\sum_{g\in\Gamma^{*}_{c}}\phi^{v}_{g}G^{+,*}_{g,t}+\sum_{j\in\mathcal{B}}\Big{(}\phi^{u}_{j}D^{*}_{j,t}+\\\
&\;\sum_{g\in\Gamma_{j}\cap\Gamma^{*}_{c}}\phi^{o}_{g}G^{-,*}_{g,t}+\sum_{g\in\Gamma_{j}\cap\Gamma_{r}}\phi^{c}_{g}G^{-,*}_{g,t}\Big{)}\Big{)}\end{split}$
(3)

The function $f_{u}$ accounts for the startup cost $\phi^{s}_{g}$ and the no-
load cost $\phi^{f}_{g}$ for generator $g$ at time $t$. The no-load cost is
multiplied by $\Delta^{*}$ since the DA-UC and ST-UC levels have different
time resolutions. The function $f_{e}$ accounts for the variable cost,
$\phi^{v}_{g}$, of the energy consumed by the grid, the cost of
overgeneration, $\phi^{o}_{g}$, and the cost of curtailment, $\phi^{c}_{g}$,
for generator $g$ at time $t$. It also accounts for the cost of unmet demand
$\phi^{u}_{j}$. We next define constraints for the layers:

$\displaystyle\begin{split}c_{d}(\boldsymbol{y}_{i}^{*}):=&\;\Big{(}\sum_{j\in\mathcal{B}:(j,k)\in\mathcal{L}}F^{*}_{j,k,t}-\sum_{j\in\mathcal{B}:(k,j)\in\mathcal{L}}F^{*}_{k,j,t}+\\\
&\;\sum_{g\in\Gamma_{k}}G^{+,*}_{g,t}+D^{*}_{k,t}-\hat{D}^{*}_{k,t}-\hat{R}^{*}_{k,t}\Big{)}_{k\in\mathcal{B},t\in\mathcal{T}_{i}^{*}}\end{split}$
(4)

This requires that, at each time point, the flows coming into the bus, the
power consumed by the grid, and the amount of unmet demand is equal to the
demand, $\hat{D}^{*}_{k,t}$, and the reserve requirements, $\hat{R}^{*}_{k,t}$
for $k\in\mathcal{B}$. The power flow equations use a DC approximation:

$\displaystyle
c_{f}(\boldsymbol{y}^{*}_{i}):=\left(F^{*}_{j,k,t}-B_{j,k}(\theta^{*}_{j,t}-\theta^{*}_{k,t})\right)_{(j,k)\in\mathcal{L},t\in\mathcal{T}^{*}_{i}}$
(5)

The renewable resources are also restricted to a specific value; this is
enforced by constraint $c_{r}$, where $\hat{G}^{*}_{g,t}$ is the amount of
power produced by renewable generator $g$ at time $t$:

$\displaystyle
c_{r}(\boldsymbol{y}^{*}_{i}):=\left(G^{+,*}_{g,t}+G^{-,*}_{g,t}-\hat{G}^{*}_{g,t}\right)_{g\in\Gamma_{r},t\in\mathcal{T}^{*}_{i}}$
(6)

The ramp-up and ramp-down constraints ($c_{ru}$ and $c_{rd}$) have different
forms because DA-UC and ST-UC/HA-ED have different time resolutions. Ramping
constraints containing start-up and shut-down constraints are:

$\displaystyle\begin{split}c_{su}(\boldsymbol{y}^{*}_{i},\boldsymbol{x}^{**}_{j}):=&\;\Big{(}G^{+,*}_{g,t}+G^{-,*}_{g,t}-G^{+,*}_{g,t-\Delta^{*}}-G^{-,*}_{g,t-\Delta^{*}}-\\\
&\;(\overline{S}_{g}-\overline{R}_{g}\Delta^{**}-\underline{C}_{g})s^{**}_{g,t+\Delta^{**}}-\\\
&\;(\overline{R}_{g}\Delta^{**}+\underline{C}_{g})x^{**}_{g,t+\Delta^{**}}+\underline{C}_{g}x^{**}_{g,t}\Big{)}_{g\in\Gamma_{c}^{**},t\in\Theta}\end{split}$
(7)
$\displaystyle\begin{split}c_{sd}(\boldsymbol{y}^{*}_{i},\boldsymbol{x}^{**}_{j}):=&\;\Big{(}G^{+,*}_{g,t-\Delta^{*}}+G^{-,*}_{g,t-\Delta^{*}}-G^{+,*}_{g,t}-G^{-,*}_{g,t}-\\\
&\;(\underline{S}_{g}-\underline{R}_{g}\Delta^{*}-\underline{C}_{g})z^{**}_{g,t+\Delta^{**}}-\\\
&\;(\underline{R}_{g}\Delta^{**}+\underline{C}_{g})x^{**}_{g,t}+\underline{C}_{g}x^{**}_{g,t+\Delta^{**}}\Big{)}_{g\in\Gamma_{c}^{**},t\in\Theta}\end{split}$
(8)

Here, $\overline{S}_{g}$ and $\underline{S}_{g}$ are the startup/shutdown
limits for $g$, $\overline{R}_{g}$ and $\underline{R}_{g}$ are the ramp-up and
ramp-down limits for $g$ as a function of time, and $\underline{C}_{g}$ is the
minimum capacity of $g$. The set $\Theta$ is defined in this context as
$\Theta=\\{t_{i}:t_{i}-\Delta^{*}=t_{j}-\Delta^{**},\forall
t_{i}\in\bar{\mathcal{T}}^{*}_{i},t_{j}\in\mathcal{T}_{j}^{**}\\}$. The model
uses constraints
$c_{su,0}(\boldsymbol{y}^{*}_{i},\boldsymbol{y}^{h}_{k},\boldsymbol{x}_{j}^{**},\boldsymbol{x}_{j-1}^{**})$
and
$c_{sd,0}(\boldsymbol{y}^{*}_{i},\boldsymbol{y}^{h}_{k},\boldsymbol{x}_{j}^{**},\boldsymbol{x}_{j-1}^{**})$
to link the first time point of the $i$th subproblem with the solutions of the
previoussubproblems (this introduces complex time coupling). We also define
operating regions:

$\displaystyle\mathcal{X}_{i}^{*}$
$\displaystyle:=\left\\{\boldsymbol{x}^{*}_{i}|(x_{g,t},s_{g,t},z_{g,t})\in\\{0,1\\}^{3},\quad\begin{subarray}{c}\forall
g\in\Gamma_{c}^{*},\\\ t\in\mathcal{T}_{i}^{*}\end{subarray}\right\\}$ (9)
$\displaystyle\mathcal{Y}^{*}$
$\displaystyle:=\left\\{\begin{array}[]{c|l}&\underline{\theta}_{j}\leq\theta_{j,t}^{d}\leq\overline{\theta}_{j},\quad\forall
j\in\mathcal{B},t\in\mathcal{T}_{i}^{d}\\\ &\underline{F}_{j,k}\leq
F^{*}_{j,k,t}\leq\overline{F}_{j,k},\quad\forall(j,k)\in\mathcal{L},t\in\mathcal{T}_{i}^{d}\\\
\boldsymbol{y}^{*}_{i}&(G^{+,*}_{g,t},G^{-,*}_{g,t})\in\mathbb{R}_{+}^{2},\quad\forall
g\in\Gamma_{c}^{d}\cup\Gamma_{r},t\in\mathcal{T}_{i}^{d}\\\
&F^{*}_{j,k,t}\in\mathbb{R},\quad\forall(j,k)\in\mathcal{L},t\in\mathcal{T}^{d}_{i}\\\
&\theta^{*}_{j,t}\in\mathbb{R},D^{*}_{j,t}\in\mathbb{R}_{+},\quad\forall
j\in\mathcal{B},t\in\mathcal{T}^{*}_{i}\end{array}\right\\}$ (15)

The $i$th DA-UC subproblem is given by (16)

$\displaystyle\min$
$\displaystyle\;f_{u}(\boldsymbol{x}^{d}_{i})+f_{e}(\boldsymbol{y}_{i}^{d})\Delta^{d}$
(16a) s.t. $\displaystyle\;c_{d}(\boldsymbol{y}^{d}_{i})=0,\quad
c_{f}(\boldsymbol{y}^{d}_{i})=0,\quad c_{r}(\boldsymbol{y}^{d}_{i})=0$ (16b)
$\displaystyle\;c_{su}(\boldsymbol{y}^{d}_{i},\boldsymbol{x}^{d}_{i})\leq
0,\quad
c_{su,0}(\boldsymbol{y}^{d}_{i},{\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{k}},\boldsymbol{x}^{d}_{i},\boldsymbol{x}^{d}_{i-1})\leq
0$ (16c)
$\displaystyle\;c_{sd}(\boldsymbol{y}^{d}_{i},\boldsymbol{x}^{d}_{i})\leq
0,\quad
c_{sd,0}(\boldsymbol{y}^{d}_{i},{\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{k}},\boldsymbol{x}^{d}_{i},\boldsymbol{x}^{d}_{i-1})\leq
0$ (16d) $\displaystyle\;x_{g,t}-x_{g,t-1}=s_{g,t}-z_{g,t}\quad\forall
g\in\Gamma_{c}^{d},t\in\mathcal{T}^{d}_{i}$ (16e)
$\displaystyle\;\sum_{j\in{\mathcal{U}^{d}_{g,t}}}s_{g,j}\leq
x_{g,t}\quad\forall g\in\Gamma_{c}^{d},t\in\mathcal{T}^{d}_{i}$ (16f)
$\displaystyle\;\sum_{j\in\mathcal{D}^{d}_{g,t}}^{t}z_{g,j}\leq
1-x_{g,t}\quad\forall g\in\Gamma_{c}^{d},t\in\mathcal{T}_{i}^{d}$ (16g)
$\displaystyle\;\underline{C}_{g}x_{g,t}\leq
G^{+,d}_{g,t}+G^{-,d}_{g,t}\leq\overline{C}_{g}x_{g,t},\quad\begin{subarray}{c}\forall
g\in\Gamma_{c}^{d},\\\ t\in\mathcal{T}_{i}^{d}\end{subarray}$ (16h)

The $i$th subproblem of ST-UC is given by (17). This has a similar structure
as DA-UC; however, there are now DA-UC variables that are incorporated into
this lower layer solution through (17c), (17d), (17k), and (17m). The last
constraint ensures that the generation amounts for the DA-UC generators in ST-
UC are within a certain bound ($\epsilon^{s}_{g}$) of the DA-UC subproblem
solutions. This helps avoid myopic solutions, since the ST-UC time horizon is
shorter than that of DA-UC.

$\displaystyle\color[rgb]{1,0,0}\min$
$\displaystyle\;\color[rgb]{1,0,0}f_{u}(\boldsymbol{x}^{s}_{i})+f_{e}(\boldsymbol{y}_{i}^{s})\Delta^{s}$
(17a) s.t.
$\displaystyle\;\color[rgb]{1,0,0}c_{d}(\boldsymbol{y}^{s}_{i})=0,\quad\color[rgb]{1,0,0}c_{f}(\boldsymbol{y}^{s}_{i})=0,\quad\color[rgb]{1,0,0}c_{r}(\boldsymbol{y}^{s}_{i})=0$
(17b)
$\displaystyle\;\color[rgb]{1,0,0}c_{su}(\boldsymbol{y}^{s}_{i},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}})\leq
0,\quad
c_{su,0}(\boldsymbol{y}^{s}_{i},{\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{k}},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j-1}})\leq
0$ (17c)
$\displaystyle\;\color[rgb]{1,0,0}c_{sd}(\boldsymbol{y}^{s}_{i},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}})\leq
0,\quad
c_{sd,0}(\boldsymbol{y}^{s}_{i},{\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{k}},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j-1}})\leq
0$ (17d)
$\displaystyle\;\color[rgb]{1,0,0}c_{su}(\boldsymbol{y}^{s}_{i},\boldsymbol{x}^{s}_{i})\leq
0,\quad
c_{su,0}(\boldsymbol{y}^{s}_{i},{\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{k}},\boldsymbol{x}^{s}_{i},\boldsymbol{x}^{s}_{i-1})\leq
0$ (17e)
$\displaystyle\;\color[rgb]{1,0,0}c_{sd}(\boldsymbol{y}^{s}_{i},\boldsymbol{x}^{s}_{i})\leq
0,\quad
c_{sd,0}(\boldsymbol{y}^{s}_{i},{\color[rgb]{0,0,1}\boldsymbol{y}^{h}_{k}},\boldsymbol{x}^{s}_{i},\boldsymbol{x}^{s}_{i-1})\leq
0$ (17f) $\displaystyle\;\color[rgb]{1,0,0}c_{ru}(\boldsymbol{y}^{s}_{i})\leq
0,\quad\quad\color[rgb]{1,0,0}c_{rd}(\boldsymbol{y}^{s}_{i})\leq 0$ (17g)
$\displaystyle\;\color[rgb]{1,0,0}x_{g,t}-x_{g,t-1}=s_{g,t}-z_{g,t}\quad\forall
g\in\Gamma_{c}^{s},t\in\mathcal{T}^{s}_{i}$ (17h)
$\displaystyle\;\color[rgb]{1,0,0}\sum_{j\in\mathcal{U}^{s}_{g,t}}s_{g,j}\leq
x_{g,t}\quad\forall g\in\Gamma_{c}^{s},t\in\mathcal{T}^{s}_{i}$ (17i)
$\displaystyle\;\color[rgb]{1,0,0}\sum_{j\in\mathcal{D}^{s}_{g,t}}z_{g,j}\leq
1-x_{g,t}\quad\forall g\in\Gamma_{c}^{s},t\in\mathcal{T}_{i}^{s}$ (17j)
$\displaystyle\;\color[rgb]{1,0,0}\underline{C}_{g}{\color[rgb]{0,0,0}x_{g,\lceil
t\rceil}}\leq
G^{+,s}_{g,t}+G^{-,s}_{g,t}\leq\overline{C}_{g}{\color[rgb]{0,0,0}x_{g,\lceil
t\rceil}}\quad\begin{subarray}{c}\forall
g\in{\color[rgb]{0,0,0}\Gamma_{c}^{d}},\\\
t\in\mathcal{T}_{i}^{s}\end{subarray}$ (17k)
$\displaystyle\;\color[rgb]{1,0,0}\underline{C}_{g}x_{g,t}\leq
G^{+,s}_{g,t}+G^{-,s}_{g,t}\leq\overline{C}_{g}x_{g,t}\quad\begin{subarray}{c}\forall
g\in\Gamma_{c}^{s},\\\ t\in\mathcal{T}_{i}^{s}\end{subarray}$ (17l)
$\displaystyle\;\color[rgb]{1,0,0}|G^{+,s}_{g,t}+G^{-,s}_{g,t}-{\color[rgb]{0,0,0}G^{+,d}_{g,\lceil
t\rceil}}-{\color[rgb]{0,0,0}G^{-,d}_{g,\lceil
t\rceil}}|\leq\epsilon^{s}_{g},\quad\begin{subarray}{c}\forall
g\in{\color[rgb]{0,0,0}\Gamma_{c}^{d}},\\\
t\in\mathcal{T}_{i}^{s}\end{subarray}$ (17m)

The $i$th subproblem of HA-ED is given by (18). This formulation is similar to
ST-UC but without the binary variables and their accompanying constraints and
objective function. In addition, there are now links between both the ST-UC
and DA-UC layers.

$\displaystyle\color[rgb]{0,0,1}\min$
$\displaystyle\;\color[rgb]{0,0,1}f_{e}(\boldsymbol{y}_{i}^{h})\Delta^{h}$
(18a) s.t.
$\displaystyle\;\color[rgb]{0,0,1}c_{d}(\boldsymbol{y}^{h}_{i})=0,\quad
c_{f}(\boldsymbol{y}^{h}_{i})=0,\quad c_{r}(\boldsymbol{y}^{h}_{i})=0$ (18b)
$\displaystyle\;\color[rgb]{0,0,1}c_{su}(\boldsymbol{y}^{h}_{i},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}})\leq
0,\quad
c_{su,0}(\boldsymbol{y}^{h}_{i},\boldsymbol{y}^{h}_{i-1},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j-1}})$
(18c)
$\displaystyle\;\color[rgb]{0,0,1}c_{sd}(\boldsymbol{y}^{h}_{i},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}})\leq
0,\quad
c_{sd,0}(\boldsymbol{y}^{h}_{i},\boldsymbol{y}^{h}_{i-1},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j}},{\color[rgb]{0,0,0}\boldsymbol{x}^{d}_{j-1}})$
(18d)
$\displaystyle\;\color[rgb]{0,0,1}c_{su}(\boldsymbol{y}^{h}_{i},{\color[rgb]{1,0,0}\boldsymbol{x}^{s}_{i}})\leq
0,\quad
c_{su,0}(\boldsymbol{y}^{h}_{i},\boldsymbol{y}^{h}_{i-1},{\color[rgb]{1,0,0}\boldsymbol{x}^{d}_{j}},{\color[rgb]{1,0,0}\boldsymbol{x}^{d}_{j-1}})$
(18e)
$\displaystyle\;\color[rgb]{0,0,1}c_{sd}(\boldsymbol{y}^{h}_{i},{\color[rgb]{1,0,0}\boldsymbol{x}^{s}_{i}})\leq
0,\quad
c_{sd,0}(\boldsymbol{y}^{h}_{i},\boldsymbol{y}^{h}_{i-1},{\color[rgb]{1,0,0}\boldsymbol{x}^{d}_{j}},{\color[rgb]{1,0,0}\boldsymbol{x}^{d}_{j-1}})$
(18f) $\displaystyle\;\color[rgb]{0,0,1}c_{ru}(\boldsymbol{y}^{h}_{i})\leq
0,\quad\quad c_{rd}(\boldsymbol{y}^{h}_{i})\leq 0$ (18g)
$\displaystyle\;\color[rgb]{0,0,1}\underline{C}_{g}{\color[rgb]{0,0,0}x_{g,\lceil
t\rceil}}\leq
G^{+,h}_{g,t}+G^{-,h}_{g,t}\leq\overline{C}_{g}{\color[rgb]{0,0,0}x_{g,\lceil
t\rceil}}\quad\begin{subarray}{c}\forall
g\in{\color[rgb]{0,0,0}\Gamma_{c}^{d}},\\\
t\in\mathcal{T}_{i}^{h}\end{subarray}$ (18h)
$\displaystyle\;\color[rgb]{0,0,1}\underline{C}_{g}{\color[rgb]{1,0,0}x_{g,t}}\leq
G^{+,h}_{g,t}+G^{-,h}_{g,t}\leq\overline{C}_{g}{\color[rgb]{1,0,0}x_{g,t}}\quad\begin{subarray}{c}\forall
g\in{\color[rgb]{1,0,0}\Gamma_{c}^{s}},\\\
t\in\mathcal{T}_{i}^{h}\end{subarray}$ (18i)
$\displaystyle\;\color[rgb]{0,0,1}|G^{+,h}_{g,t}+G^{-,h}_{g,t}-{\color[rgb]{0,0,0}G^{+,d}_{g,\lceil
t\rceil}}-{\color[rgb]{0,0,0}G^{-,d}_{g,\lceil
t\rceil}}|\leq\epsilon_{g},\quad\begin{subarray}{c}\forall
g\in{\color[rgb]{0,0,0}\Gamma_{c}^{d}},\\\
t\in\mathcal{T}_{i}^{h}\end{subarray}$ (18j)
$\displaystyle\;\color[rgb]{0,0,1}|G^{+,h}_{g,t}+G^{-,h}_{g,t}-{\color[rgb]{1,0,0}G^{+,s}_{g,\lceil
t\rceil}}-{\color[rgb]{1,0,0}G^{-,s}_{g,t}}|\leq\epsilon_{g},\quad\forall\begin{subarray}{c}g\in{\color[rgb]{1,0,0}\Gamma_{c}^{s}},\\\
t\in\mathcal{T}_{i}^{h}\end{subarray}.$ (18k)

###  III-B Graph Representation

We now outline how we represent the tri-level hierarchical architecture in
Plasmo.jl. We represent each time point $t$ as a subgraph, and nodes are
placed on this subgraph for each bus and each transmission line. The nodes
corresponding to buses contain the variables $D^{*}_{k,t}$,
$\theta^{*}_{k,t}$, $G^{+,*}_{g,t}$, and $G^{-,*}_{g,t},\forall
g\in\Gamma_{k}\cap\\{\Gamma_{c}^{*}\cup\Gamma_{r}\\}$. In the case of DA-UC
and ST-UC problems, the bus nodes also contain $x_{g,t},s_{g,t}$, and
$z_{g,t}\>\forall g\in\Gamma_{k}\cap\Gamma_{c}^{*}$ and any constraints for
these variables. The nodes corresponding to transmission lines contain
variables $F^{*}_{i,j,t},\forall(i,j)\in\mathcal{L}$. For each node
representing line $(i,j)\in\mathcal{L}$, edges (linking constraints) are also
placed connecting to bus $i$ and to bus $j$. The resulting subgraph is shown
in Figure 4. The DA-UC, ST-UC, and HA-ED subproblems were constructed from
these time point subgraphs. The DA-UC subgraph has 24 time point subgraphs
(i.e., 24 replicates of the network shown in Figure 4) each representing one
hour, the ST-UC subgraph had 16 time point subgraphs with each representing 15
minutes, and the HA-ED subgraph had 5 time point subgraphs with each
representing 15 minutes. Linking constraints were also placed between time
point subgraphs where applicable, such as for $c_{ru}$, $c_{rd}$, (16e) -
(16g), or (17h) - (17j). After subproblem subgraphs were created, the
subproblems were combined onto another OptiGraph corresponding to one day of
operation. A single day graph contains one DA-UC subgraph, 8 ST-UC subgraphs,
and 96 HA-ED subgraphs (105 subgraphs in total). Figure 5 shows the complexity
of the resulting graph.

Figure 4: A graph visualization of the 118-bus system, where nodes correspond
to buses and to transmission lines. Figure 5: Representation of a single
subproblem subgraph for the DA-UC, ST-UC and HA-ED subproblems. Figure 6: An
example of the hierarchical linking between subproblems. Five time points from
an ST-UC subproblem are shown linked to a full HA-ED subproblem, with linking
constraints highlighted in black.

With the single-day graph formed, subproblems are linked together according to
the formulations given in (16), (17), and (18). Figure 6 shows an example of
this linking for part of an ST-UC subproblem and one HA-ED subproblem. The
linking constraints are highlighted in black; these constraints correspond to
$c_{su}$ and $c_{sd}$ in (18e) and (18f) and to the linking constraints in
(18i) and (18k). The full problem graph is shown in Figure 7a, with
accompanying representations highlighting the hierarchical structure. The full
graph contains 192,432 OptiNodes and 292,587 OptiEdges. Subgraphs can be
collapsed or aggregated into OptiNodes without changing the problem
formulation and this facilitates visualization. Figure 7b shows the graph with
all time subgraphs aggregated into nodes. Figure 7c shows all subproblem
subgraphs aggregated into nodes; this reveals the hierarchical structure and
the linking between layers, where the central black node corresponds to the
single-day DA-UC subproblem (top layer), the red nodes correspond to the 8 ST-
UC subproblems (middle layer), and the blue nodes correspond to the 96 HA-ED
subproblems (bottom layer).

((a))

((b))

((c))

((d))

Figure 7: Representation of (a) the 1-day monolithic graph highlighted by the
subproblems (DA-UC, ST-UC, HA-ED), and equivalent representations where (b)
individual time point subgraphs are aggregated into single nodes and where (c)
subproblem subgraphs are aggregated into single nodes.

###  III-C Decomposition Approaches

Graph representations facilitate the implementation of different decomposition
approaches. For example, [2] decomposed the hierarchy by solving the
subproblems in each layer in series (in a receding horizon approach). This
sequential decomposition approach can be easily visualized using graphs. For
Figure 7c, this is equivalent to solving the central DA-UC node, then solving
the first ST-UC node, and then solving the 12 connected HA-ED nodes. The next
ST-UC node is then solved followed by its 12 connected HA-ED nodes (again in
order) and so forth until all 8 ST-UC subproblems and their corresponding HA-
ED nodes are solved. The solutions of these problems are then passed to the
next 1-day graph and the process is repeated. However, there are other
decomposition approaches that could be used; for example, instead of solving
each subproblem in a receding horizon approach, we could solve each time-point
subgraph in a receding horizon approach. This would result in much smaller
optimization problems, but likely worse economic performance. In contrast, we
could instead solve the entire 1-day monolithic problem as a single
optimization problem rather than solving each supbroblem one at a time. The
implementation of these strategies can help study trade-offs between
tractability and performance.

###  III-D Results

In this section, we present the results for two different solution approaches.
The first approach is to solve the subproblems in a receding horizon approach
as done in [2] (”Receding Horizon”). The second approach is to solve each
1-day graph as a single, monolithic optimization problem (”Monolithic”).
Because the monlithic problem is a very large mixed-integer problem (MIP), it
took hours to solve, so we used a MIP gap termination criteria of 5%. In
contrast, we used a MIP gap of 0.5% or less for the receding horizon MIPs as
they were smaller and faster to solve. The 1-day monolithic graph contained
641,709 variables (41,727 binary) and 1,103,654 constraints. The code for
reproducing these results can be found at
https://github.com/zavalab/JuliaBox/tree/master/hierarchical_graphs. We used
the data provided by the 118-bus case study [15] used in [2]. This included a
day-ahead (forecasted) load demand, and a real time realized load demand. We
used the day-ahead demand for the DA-UC subproblems and the real time demand
for the HA-ED subproblems. Because there is no intermediate “short term”
demand data, we used the average of the day ahead and real time demands for
the ST-UC subproblems. The three load demands are shown in Figure 8. This data
was on an hour resolution, so we interpolated the data for higher resolutions.
In addition, reserve requirements can vary by system operator, but we chose to
use 10% of the demand for the reserve requirement for UC subproblems and 2.5%
of the demand for ED subproblems which corresponds to the “low reserve
requirements” scenario in [2].

Figure 8: Load demand used in the DA-UC, ST-UC, and HA-ED subproblems for the
32-day horizon.

The results of the receding horizon and monolithic approaches are shown in
Figure 9 which includes the number of committed DA-UC generators (a), the
number of committed ST-UC generators (b), the overgenerated or curtailed power
(c), and the amount of load shed (d). The overgenerated/curtailed power and
the load shed shown are from the first time point of each HA-ED subproblem. As
each HA-ED subproblem had significant overlap with the next problem, we only
consider the first HA-ED time point (this is the realized operation). The
overall cost of economic dispatch (based on the first time point of each HA-ED
subproblem) was $ 219.4 million and $ 199.7 million for the receding horizon
and monolithic approaches, respectively.

Figure 9: Results of the tri-level problem using a ”receding horizon” type
approach and a ”monolothic” approach. a) number of committed DA-UC generators
over time, b) number of committed ST-UC generators over time, d) overgenerated
or curtailed power, and d) load shedding.

###  III-E Discussion

By constructing the hierarchical architecture as a graph, different solution
schemes were enabled which provide different insights into the problem.
Despite the higher MIP gap used for the monolithic approach, it still
performed better than the receding horizon problem and had a lower cost in the
economic dispatch by more than $19 million. The monolithic approach was
expected to perform better as the lower layers and upper layers are solved in
the same problem, allowing the performance of the lower layer to inform the
upper layers. This is also likely why the monolithic approach has less load
shedding compared with the receding horizon approach (83.6 MWhr compared with
3113.5 MWhr). The monolithic approach did have a very large peak of
overgenerated/curtailed power, but the cost of load shedding (using the costs
from [2]) was 200 times more than the cost of overgenerated/curtailed power.
In addition, the monolithic approach had less fluctuation in the number of
generators turned on or off.

The results on load shedding were dependent on the reserve requirements used.
In this case, the higher reserve requirements on the UC layers compared with
the ED layer (10% vs. 2.5 %) reduced some of the apparent differences between
the demand in the HA-ED layer and the DA-UC layer (e.g., the gap between
demand in the DA-UC and the HA-ED layers in day 20 in Figure 8 would be
reduced). If we adjust the reserve requirements and use the “very low reserve
requirements” scenario from [2] (5% of load for UC layer and 1.25% of load for
ED layer), the load shed in the serial problem increases by more than 10
times. While not tested, it is possible that further increasing the reserve
requirements could reduce load shedding and/or overgeneration.

As expected, the monolithic approach took much longer to solve than the
receding horizon decomposition approach. In addition, the monolithic approach
experienced complications with memory management in the MIP solver. These
computational issues, combined with the performance comparisons between the
receding horizon and monolithic approaches, highlight the need for
decomposition schemes. The receding horizon approach results in a suboptimal
solution, but it could be possible to use a decomposition scheme that gives
results closer to the monolithic approach but with the computational
performance closer to that of the receding horizon problem. Constructing these
problems as graphs provides a framework under which a decomposition scheme
could be optimized. Overall, this work highlights the utility of representing
hierarchical optimization problems using graphs. These graph representations
provide a modular way to construct complex (but structured) problems. Each
time point can be constructed in a modular manner, and then each time point
can be embedded to a modular representation of each subproblem. Graphs are
also intuitive to visualize, potentially leading to insights into the problem
structure. They provide a framework for manipulating problem structure, such
as partitioning/aggregating subgraphs. Graphs also provide a structure that
could be exploited via decomposition schemes such as Benders decomposition and
Lagrangian relaxation.

##  IV Conclusions and Future Work

We discussed how hierarchical optimization problems can be represented with
graphs. We used the package Plasmo.jl to build a tri-level hierarchical
optimization problem arising in market operations and presented different
approaches to solve the problem. We presented visualizations of these graph
representations in Plasmo.jl, and we presented the results of the two solution
approaches. As part of future work, we are interested in using the graph
representation for applying and combining decomposition schemes (e.g.,
Lagrangian decomposition, Benders decomposition, dual dynamic integer
programming) to solve large-scale problem instances.

##  V Acknowledgements

This work was supported by the U.S. Department of Energy under grant
DE-0002722.

## References

  * [1] A. J. Conejo and L. Baringo, _Unit Commitment and Economic Dispatch_. Cham: Springer International Publishing, 2018, pp. 197–232.
  * [2] S. Atakan, H. Gangammanavar, and S. Sen, “Towards a sustainable power grid: Stochastic hierarchical planning for high renewable integration,” _European Journal of Operational Research_ , vol. 302, no. 1, pp. 381–391, 2022.
  * [3] F. Guo, C. Wen, J. Mao, J. Chen, and Y.-D. Song, “Hierarchical decentralized optimization architecture for economic dispatch: A new approach for large-scale power system,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 2, pp. 523–534, 2017.
  * [4] C. Kong, R. Jovanovic, I. S. Bayram, and M. Devetsikiotis, “A hierarchical optimization model for a network of electric vehicle charging stations,” _Energies_ , vol. 10, no. 5, p. 675, 2017.
  * [5] J. Jalving, S. Shin, and V. M. Zavala, “A graph-based modeling abstraction for optimization: Concepts and implementation in plasmo. jl,” _Mathematical Programming Computation_ , vol. 14, no. 4, pp. 699–747, 2022.
  * [6] D. L. Cole, S. Shin, and V. M. Zavala, “A julia framework for graph-structured nonlinear optimization,” _Industrial & Engineering Chemistry Research_, vol. 61, no. 26, pp. 9366–9380, 2022.
  * [7] M. Berger, A. Bolland, B. Miftari, H. Djelassi, and D. Ernst, “Graph-based optimization modeling language: A tutorial,” _ORBi_ , 2021. [Online]. Available: https://hdl.handle.net/2268/256705
  * [8] M. Berger, D. Radu, G. Detienne, T. Deschuyteneer, A. Richel, and D. Ernst, “Remote renewable hubs for carbon-neutral synthetic fuel production,” _Frontiers in Energy Research_ , vol. 9, p. 671279, 2021.
  * [9] A. Allman, W. Tang, and P. Daoutidis, “Decode: a community-based algorithm for generating high-quality decompositions of optimization problems,” _Optimization and Engineering_ , vol. 20, no. 4, pp. 1067–1084, 2019.
  * [10] W. E. Hart, J.-P. Watson, and D. L. Woodruff, “Pyomo: modeling and solving mathematical programs in python,” _Mathematical Programming Computation_ , vol. 3, no. 3, pp. 219–260, 2011.
  * [11] W. E. Hart, C. D. Laird, J.-P. Watson, D. L. Woodruff, G. A. Hackebeil, B. L. Nicholson, and J. D. Siirola, _Pyomo-optimization modeling in python_ , ser. Springer Optimization and Its Applications. Cham: Springer, 2017, vol. 67.
  * [12] P. Daoutidis, W. Tang, and A. Allman, “Decomposition of control and optimization problems by network structure: Concepts, methods, and inspirations from biology,” _AIChE Journal_ , vol. 65, no. 10, p. e16708, 2019.
  * [13] S. Shin, C. Coffrin, K. Sundar, and V. M. Zavala, “Graph-based modeling and decomposition of energy infrastructures,” _IFAC-PapersOnLine_ , vol. 54, no. 3, pp. 693–698, 2021.
  * [14] S. Shin, V. M. Zavala, and M. Anitescu, “Decentralized schemes with overlap for solving graph-structured optimization problems,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 3, pp. 1225–1236, 2020.
  * [15] I. Pena, C. B. Martinez-Anido, and B.-M. Hodge, “An extended ieee 118-bus test system with high renewable penetration,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 281–289, 2017.