# Risk-Averse Self-Scheduling of Storage in Decentralized Markets ††thanks:
Mr. Yurdakul gratefully acknowledges the support of the German Federal
Ministry of Education and Research and the Software Campus program under Grant
01IS17052. Mr. Billimoria’s Doctor of Philosophy at the University of Oxford
is supported by the Commonwealth PhD Scholarship, awarded by the Commonwealth
Scholarship Commission of the UK (CSCUK).

Ogun Yurdakul12 and Farhad Billimoria3  1Energy Systems and Infrastructure
Analysis Division, Argonne National Laboratory, Lemont, IL 60439, USA.
2Department of Electrical Engineering and Computer Science, Technical
University of Berlin 10623 Berlin, Germany  3Department of Engineering
Science, University of Oxford, Oxford OX1 2JD, United Kingdom

###### Abstract

Storage is expected to be a critical source of firming in low-carbon grids. A
common concern raised from ex-post assessments is that storage resources can
fail to respond to strong price signals during times of scarcity. While
commonly attributed to forecast error or failures in operations, we posit that
this behavior can be explained from the perspective of risk-averse scheduling.
Using a stochastic self-scheduling framework and real-world data harvested
from the Australian National Electricity Market, we demonstrate that risk-
averse storage resources tend to have a myopic operational perspective, that
is, they typically engage in near-term price arbitrage and chase only few
extreme price spikes and troughs, thus remaining idle in several time periods
with markedly high and low prices. This has important policy implications
given the non-transparency of unit risk aversion and apparent risk in
intervention decision-making.

###### Index Terms:

electricity markets, electricity storage, self-scheduling, risk-aversion.

##  I Introduction

In renewable-dominated power systems, electricity storage is expected to play
an important role in maintaining system balance and resource adequacy. In
decentralized markets, system operators (SOs) do not issue commitment
instructions, instead relying upon full-strength price formation and iterative
pre-dispatch systems to signal the need for additional or reduced resource
commitment. During periods of scarcity, understanding the projected available
capacity of resources to service load at intraday system stress points is
critical to the SO’s decision-making on whether to intervene or apply
emergency protocols, such as load shedding.

Storage resources introduce new complexity to the problem because while the
resource may be operationally available, whether it can actually charge or
discharge at a point in time depends upon its state of charge (driven by prior
decisions) and technical constraints. As a case in point, the lack of
transparency on the availability of energy-limited plants was one of the key
drivers of the Australian Energy Market Operator’s (AEMO) decision to suspend
the National Electricity Market (NEM) in June 2022 during a period of extreme
scarcity. In certain cases, storage resources were observed to not discharge
during the highest prices of the day (while discharging at more muted prices
earlier in the day).111Our online companion [1] provides an empirical record
of many such observations of battery dispatch and energy prices for the NEM.
This seemingly non-intuitive behavior is commonly attributed to error or
failure (e.g., in forecasting or operational management) of storage
operations. Our work however proposes a second explanation: that of risk-
averse self-scheduling under uncertainty, a concomitant of which is that such
behavior could well be rational and utility-maximizing, and thus persist in
energy markets.

In the literature, risk-averse self-scheduling of thermal resources has been
extensively studied. Risk-constrained self-scheduling is explored in [2, 3]
given the internalization of the non-convex attributes of thermal resources.
In [4], the conditional Value-at-Risk measure is used to develop bidding
curves derived from an optimal self-schedule. For storage and in relation to
decentralized markets, inter-temporal tradeoffs are considered from the
perspective of look-ahead bidding to optimize state-of-charge and maximize
profits [5]. In [6] a risk-constrained problem is formulated for a pumped
storage plant to trade-off between expected profit and risks in energy and
ancillary service markets.

In this paper, we work out a novel risk-averse stochastic optimization
framework for the self-scheduling of storage in decentralized markets under
uncertain prices. In Section II, we develop the mathematical background and
market architecture, leading to the precise mathematical formulation of the
risk-averse scheduling problem in Section III. In Section IV, we conduct
several experiments using real-world actual data from the NEM and storage
resources bidding therein.

Using the insights we draw from the experiments, we lay out in Section IV two
key contributions of this work. First, we provide a novel explanation for the
seemingly non-intuitive behavior of storage. Specifically, we demonstrate how
the increasing uncertainty of prices with longer horizons can lead a risk-
averse decision-maker to adopt a rational but myopic approach to dispatch. We
illustrate how such a decision-maker can favor low-risk near-term dispatch
decisions at more moderate prices rather than higher-reward but more uncertain
peaks and troughs. Second, we present valuable insights into the sensitivity
of the expected profits to the duration and the degree of risk-aversion of a
storage resource. We observe that while increasing the capacity of a storage
resource can significantly boost profits for a risk-neutral decision-maker, it
barely makes a dent for relatively more risk-averse decision-makers. We set
out concluding remarks and policy implications in Section V.

##  II Mathematical Background

We consider a storage resource trading in a decentralized electricity market.
The adopted market setting does not include a centrally organized short-term
forward market for energy and relies solely on the real-time market (RTM) for
the spot trading of energy, settled under a marginal pricing scheme222In this
paper we are focused upon energy prices and do not consider, at this stage,
ancillary service or other spot markets for which storage is eligible. We
assume that the resource acts as a price-taker, that is, it has no capability
to exercise market power or alter market prices.

We denote by $k$ the index of each time period for which the RTM is cleared
and by $K$ the total number of time periods in the problem horizon. We
represent by $\kappa$ the duration of each time period $k$ in minutes, which
is the smallest indecomposable unit of time in our analysis, during which we
assume the system conditions hold constant. Using our notation, an RTM cleared
for a day at half-hourly intervals would correspond to $\kappa=30$ min and
$K=48$. Define the set $\mathscr{K}\coloneqq\\{k\colon k=1,\ldots,K\\}$.

We denote by $\tilde{\lambda}_{k},k\in\mathscr{K}$ as the uncertain energy
price in time period $k$ at the transmission system bus or zone into which the
storage resource is located. We construct the vector
$\bm{\tilde{\lambda}}\coloneqq[\tilde{\lambda}_{1}\cdots\tilde{\lambda}_{k}\cdots\tilde{\lambda}_{K}]^{\mathsf{T}}$.
We assume that the SO determines and publishes pre-dispatch energy prices over
the $K$ time periods in order to provide market participants and itself with
advance information necessary to plan the physical operation of the power
system. We write the relation
$\bm{\tilde{\lambda}}\coloneqq\bm{\bar{\lambda}}+\bm{\tilde{\epsilon}}$, where
$\bm{\bar{\lambda}}$ denotes the vector of pre-dispatch energy prices and
$\bm{\tilde{\epsilon}}$ is a random vector of forecast errors representing the
difference between the pre-dispatch and market-clearing prices.

Typically, the uncertain forecast error exhibits a greater degree of variance
with extending horizon [1]. While the forecast error between pre-dispatch and
market-clearing prices are likely small for the periods close to the dispatch
interval, the pre-dispatch prices that eventuate may significantly deviate
from the market price for forecast horizons well away from the actual dispatch
time. We use a set of scenarios to model $\bm{\tilde{\epsilon}}$, where
$\bm{{\epsilon}}^{\omega}$ denotes the vector of forecast errors in scenario
$\omega\in\Omega$ and ${\pi}^{\omega}$ denotes its associated probability of
realization. We define the vector of LMPs for each scenario by
$\bm{{\lambda}}^{\omega}\coloneqq\bm{\bar{\lambda}}+\bm{{\epsilon}}^{\omega}$
and write
$\bm{{\lambda}}^{\omega}\coloneqq[{\lambda}_{1}^{\omega}\cdots{\lambda}_{k}^{\omega}\cdots{\lambda}_{K}^{\omega}]$,
where ${\lambda}_{k}^{\omega}$ denotes the energy price in time period $k$ of
scenario $\omega$.

We next turn to the mathematical modeling of storage. We denote by $p^{c}_{k}$
and $p^{d}_{k}$ the charging and discharging power of the storage resource in
time period $k$ with maximum values $p^{c}_{M}$ and $p^{d}_{M}$, respectively.
To ensure the mutual exclusivity of the charging and discharging modes, and to
reduce the computational burden of the problem, we draw upon the single binary
variable storage formulation presented in [7] and enforce the following
constraints:

$\displaystyle p^{d}_{k}-p_{M}^{d}u_{k}\leq 0,$ $\displaystyle
p^{c}_{k}-p_{M}^{c}(1-u_{k})\leq 0$ $\displaystyle\forall k\in\mathscr{K}.$
(1) $\displaystyle u_{k}\in\\{0,1\\},$ $\displaystyle p^{c}_{k},p^{c}_{k}\geq
0$ $\displaystyle\forall k\in\mathscr{K}.$ (2)

We denote the efficiency of the storage unit by $\eta$ (assuming symmetry
between charging and discharging efficiencies). We represent by $E_{k}$ the
stored energy level at the end of the time period $k$ and write the
intertemporal operational constraints:

$\displaystyle E_{k}=E_{k-1}-\frac{1}{\eta}p_{k}^{d}\kappa\frac{1\text{
h}}{60\text{ min}}+{\eta}p_{k}^{c}\kappa\frac{1\text{ h}}{60\text{ min}}$
$\displaystyle\forall k\in\mathscr{K}\setminus 1,$ (3) $\displaystyle
E_{k}=E_{o}-\frac{1}{\eta}p_{k}^{d}\kappa\frac{1\text{ h}}{60\text{
min}}+{\eta}p_{k}^{c}\kappa\frac{1\text{ h}}{60\text{ min}}$
$\displaystyle\forall k\in\\{1\\},$ (4)

where $E_{o}$ denotes the initial stored energy level at the beginning of the
problem horizon. The stored energy level in each period $k$ is bounded from
above and below by $E_{M}$ and $E_{m}$, respectively:

$\displaystyle E_{m}\leq E_{k}\leq E_{M}\hskip 8.5359pt\forall
k\in\mathscr{K}.$ (5)

In the next section, we harness the mathematical models laid out in this
section to work out a framework for the risk-averse self-scheduling problem of
a storage resource.

##  III Risk-Averse Self-Scheduling

Absent a short-term forward market, the storage resource is directly exposed
to the uncertainty in the RTM prices. A great deal of studies in the
literature considers a risk-neutral market participant, for which the storage
resource may bring about large losses in some scenarios as long as those are
offset by even greater gains in other scenarios. Such studies however do not
cater to risk-averse decisions-makers, who constitute the focus of this work,
that prefer to ward off large losses independent of potential profits.

A widely used measure for incorporating risk into decision-making is the
Value-at-Risk (VaR). For a specified risk confidence level $\alpha$ in
$(0,1)$, $\alpha$-VaR is an upper estimate of losses that is exceeded with
probability $1-\alpha$. We denote the $\alpha$-VaR of the loss associated with
a decision $\bm{x}$ by $\zeta_{\alpha}(\bm{x})$. Despite presenting an
intuitive representation of losses, VaR exhibits several undesirable
properties, including not taking account of the losses suffered beyond
$\zeta_{\alpha}(\bm{x})$, non-coherence, and non-convexity when computed using
scenarios. Instead, we draw upon the CVaR measure to manage risk in our
framework.

For continuous distributions, the $\alpha$-CVaR of the loss for a decision
$\bm{x}$, which we denote by $\phi_{\alpha}(\bm{x})$, is the expected loss
given that the loss is greater than or equal to $\zeta_{\alpha}(\bm{x})$. The
definition of CVaR for discrete distributions is yet more subtle, which is the
case in our framework as we rely on scenarios to represent uncertain prices.
Rockafellar and Uryasev [8] define $\phi_{\alpha}(\bm{x})$ for general
distributions as the weighted average of $\zeta_{\alpha}(\bm{x})$ and the
expected loss strictly exceeding $\zeta_{\alpha}(\bm{x})$.

From a mathematical standpoint, CVaR presents several appealing features.
Pflug [9] shows that it is a coherent risk measure. Most notably,
$\phi_{\alpha}(\bm{x})$ can be efficiently computed by minimizing a piecewise
linear and convex function [8, Theorem 10], which can be cast as a linear
programming (LP) problem by introducing an additional variable.

A key thrust of our framework is to hedge the decision-maker against the risk
of incurring high charging costs. For an associated confidence level $\alpha$,
the $\alpha$-CVaR of the uncertain charging cost over the problem horizon can
be evaluated by solving the following LP problem:

$\displaystyle\underset{z^{\omega},\zeta}{\text{minimize}}$
$\displaystyle\mathcal{R}_{\alpha}(\bm{x})\coloneqq\zeta+\frac{1}{1-\alpha}\sum_{\omega\in\Omega}\pi^{\omega}z^{\omega},$
(6) subject to $\displaystyle
z^{\omega}\geq\sum_{k\in\mathscr{K}}{\lambda}^{\omega}_{k}p_{k}^{c}-\zeta,\;z^{\omega}\geq
0$ (7)

where the decision vector $\bm{x}$ succinctly represents the storage variables
$p_{k}^{c}$, $p_{k}^{d}$, and $E_{k}$ for all $k$ in $\mathscr{K}$, $\zeta$,
and the auxiliary variable $z^{\omega}$ introduced for each $\omega$ in
$\Omega$ and $\pi^{\omega}$ denotes the probability of the scenario
$\omega\in\Omega$.

In conjunction with minimizing the risk of incurring high charging costs,333We
would be remiss if we did not set forth that the resource faces a risk while
discharging as well, that of recording low revenues if the RTM prices
materialize at significantly lower levels than their pre-dispatch
counterparts. We choose to omit this risk in this paper and defer it to our
future work.  the decision-maker seeks to maximize the expected profits of the
resource over the $K$ time periods:

$\displaystyle\mathcal{P}(\bm{x})\coloneqq\sum_{\omega\in\Omega}\pi^{\omega}\Big{[}\sum_{k=1}^{K}\lambda^{\omega}_{k}\big{(}p^{d}_{k}-p^{c}_{k}\big{)}\Big{]}.$
(8)

We adjust the trade-off between these seemingly conflicting objectives by
minimizing the weighted combination of $\mathcal{R}(\bm{x})$ and
$\mathcal{P}(\bm{x})$ subject to the storage constraints laid out in Section
II and the CVaR constraints presented in this section. The risk-averse self-
scheduling (RASS) problem is expressed as:

$\displaystyle\text{RASS}:$ minimize $\displaystyle\hskip
14.22636pt-\mathcal{P}(\bm{x})+\beta\mathcal{R}(\bm{x}),$ subject to
$\displaystyle\hskip
14.22636pt\eqref{plim}\text{--}\eqref{enl},\eqref{cvarc1}.$

The weight parameter $\beta\in[0,\infty)$ tunes the decision-maker’s degree of
risk-aversion. While increasing values of $\beta$ underscores the decision-
maker’s desire to mitigate risk, driving down $\beta$ toward zero represents a
more risk-neutral decision-maker. Setting $\beta=0$ implies that the sole
objective of the decision-maker is to maximize her profits—independent of the
risk that her decisions entail.

The RASS problem is solved on a rolling-window basis with a shrinking horizon.
At the outset, it is solved for the time periods $k=1\cdots K$ before the
uncertain price for any of the time periods is revealed, yet the optimal
decisions for only the first time period are implemented. After the RTM price
for $k=1$ is observed, the RASS problem is solved again for $k=2\cdots K$,
this time implementing the optimal decision for only $k=2$, and so forth. The
process repeats until the RASS problem is solved for $k=K$ with a single-
period horizon. The binding decisions of each rolling window are passed along
to the subsequent window by explicitly enforcing that the stored energy level
at the end of the first, binding time period of each window be equal to the
initial stored energy level at the beginning of the ensuing window.

##  IV Case Study and Results

In this section, we conduct two case studies to gain insights into the self-
scheduling of a storage resource under different pre-dispatch price signals.
The price data of both case studies are from the NEM. As spelled out below, we
use the actual pre-dispatch prices observed in the NEM in two representative
days to form the vector of pre-dispatch prices in our experiments. To
construct the scenarios $\bm{{\lambda}}^{\omega},\,\omega\in\Omega$, we draw
upon the historical differences between the pre-dispatch and RTM prices that
eventuated in the NEM across 2019, yielding 11,680 observations, from which we
randomly select 100 observations to form the scenarios of each case study. The
price data have a temporal granularity of 30 min, that is, $\kappa=30$ min and
$K=48$. Commensurate with the price data, we harness the data of the storage
resources currently bidding in the NEM [10] in our experiments. Unless
otherwise stated, we pick the risk confidence level as $\alpha=0.95$ and the
storage efficiency as $\eta=0.85$ in all experiments. The data and the source
code of all experiments are furnished in the online companion [1].

###  IV-A Case Study I

The storage data of Case Study I are from that of the Victorian Big Battery,
which is a grid-connected storage resource in the Australian state of Victoria
with a charging/discharging power limit of 300.0 MW and an energy storage
capacity of 450.0 MWh. We use the pre-dispatch prices in Victoria for June 12,
2022, which constitutes the last day before the cumulative price threshold was
exceeded in Victoria, triggering the onset of a series of market interventions
that culminated in AEMO suspending the NEM on June 15, 2022. We start out by
solving the RASS problem for $\beta=0$ and $\beta=0.4$. We observe from Fig. 1
that the net discharging power (i.e., discharging power less charging power)
under $\beta=0$, by and large, closely follows the pre-dispatch prices,
effectively exploiting price differences across three time windows, with the
first being between $k=11$ and $k=19$, second between around $k=24$ and
$k=36$, and finally that around the price spike at $k=46$.

Figure 1: Case Study I optimal storage dispatch decisions for $\beta=0$. Each
time period $k$ has a duration of 30 minutes.

We note from Fig. 2 that whereas the risk-neutral decision-maker ($\beta=0$)
waits until the price reaches the early-hour minimum at $k=11$ to start
charging, the risk-averse decision-maker ($\beta=0.4$) charges at its maximum
limit right at the first two time periods, during which the pre-dispatch price
is higher than that at $k=11$. We attribute this seemingly counter-intuitive
behavior to the fact that when the RASS problem is solved at the beginning of
the horizon, charging at $k=1$ entails a lower risk compared to $k=11$, as the
market-clearing price for the initial time periods is expected hover closely
around the pre-dispatch price. Since the uncertainty in forecast error
increases with the length of the look-ahead period, the risk-averse decision-
maker is driven to store more energy at the beginning of the horizon vis-à-vis
the risk-neutral counterpart.

Figure 2: Case Study I optimal storage dispatch decisions for $\beta=0.4$

In this light, the risk-averse decision-maker incurs a greater lost
opportunity cost by storing a higher level of energy at $k=1$ and $k=2$ and
discharging very little before $k=11$. This is because, by doing so, it
reaches its maximum capacity at $k=12$ and foregoes its ability to store
energy when the prices are around their lowest ebb of the day at $k=13$ and
$k=14$. In contrast, the risk-neutral decision-maker manages to store 195 MWh
of energy at $k=13$ and $k=14$. Tellingly, in order to store energy when the
prices reach the minimum of the PM hours at $k=27$, the risk-averse decision
maker winds up discharging at $k=25$. Although it could have well discharged
between $k=15$ and $k=22$ at higher prices and could have recorded greater
revenues, it prefers to discharge at a time period closer to that during which
it charges again. Indeed, the dispatch decisions under $\beta=0.4$ is overall
marred by a myopic focus. The risk averse decision-maker exploits price
differences primarily at extremely short time windows, such as between $k=10$
and $k=11$ or between $k=25$ and $k=27$, unless the price significantly rises
as around $k=46$. The prevailing myopic focus under $\beta=0.4$ can be
ascribed to the increasing uncertainty in forecast error with longer horizons,
driving the risk-averse decision-maker to arbitrage primarily between shorter
windows. A ramification of this behavior is that the resource fails to respond
to pre-dispatch price signals at time periods during which the system is in
dire need. For instance, after the pre-dispatch price reaches A$481.1/MWh at
$k=33$ (more than a sevenfold increase from that at $k=27$) the risk-averse
decision-maker does not discharge any energy. Similarly, when the pre-dispatch
price precipitously falls around $k=41$, the risk-neutral decision-maker
stores 352.9 MWh, whereas the risk-averse decision-maker does not charge at
all, failing to respond to the available 65.5% price drop. As a result, after
the price spikes at $k=46$, the risk-neutral decision-maker discharges 176.5
MWh more energy compared to the risk-averse decision-maker, greatly aiding the
SO during a period of extreme scarcity.

Figure 3: Case Study I optimal net discharging levels (i.e., discharging less
charging power) under different values of the weight parameter $\beta$.

Fig. 3 makes evident that the so-called myopic focus of the decision-maker
becomes more conspicuous under increasing values of $\beta$. Around the price
spike at $k=19$, the storage resource discharges the highest level of energy
under $\beta=0$, followed by $\beta=0.1$ and $\beta=0.2$, whereas no energy is
discharged under $\beta=0.3$ and $\beta=0.4$. These decisions impart the
relatively less conservative decision-makers ($\beta=0$, $\beta=0.1$, and
$\beta=0.2$) the capability to store more energy during the price drop around
$k=24$ compared to those under $\beta=0.3$ and $\beta=0.4$. These observations
are echoed in the dispatch decisions throughout the rest of the day. For
instance, as $\beta$ is reduced, the storage resource manages to more closely
follow the pre-dispatch price signal before and around the sharp rise at
$k=46$. Under decreasing values of $\beta$, the storage resource discharges a
higher level of energy around $k=38$, gaining in turn the capability to store
more energy during the dip in prices after $k=41$, which it can discharge when
the pre-dispatch price abruptly climbs at $k=46$.

###  IV-B Case Study II

We next examine the battery dispatch decisions for a day in which the pre-
dispatch prices are highly volatile. For this purpose, we draw upon the pre-
dispatch prices in South Australia for January 16, 2019. The storage data are
taken from that of the ESCRI storage resource, which is connected to the grid
in South Australia and has a charging/discharging power limit of 30 MW and a
capacity of 8 MWh.

Figure 4: Case Study II optimal net discharging levels under different values
of the weight parameter $\beta$ for $\alpha=0.95$.

We begin by solving the RASS problem by increasing $\beta$ from $0$ to $0.5$
in $0.1$ increments. The results in Fig. 4 show that the relatively more risk-
neutral cases ($\beta=0$ and $\beta=0.1$) closely follow the pre-dispatch
price signals, charging/discharging at all of the seven price troughs/peaks of
the day. Under $\beta=0.2$ and $\beta=0.3$, however, the storage resource
prefers to arbitrage between only six and four time windows, respectively,
whereas under $\beta=0.4$, it capitalizes on only the widest price spread,
charging and discharging once throughout a highly volatile day. Most notably,
under $\beta=0.4$, the resource relinquishes the chance to take advantage of a
$56.7\%$ price difference between $k=22$ and $k=23$ when the price falls
precipitously, yet changes course three periods later and records a $126.7\%$
rise between $k=23$ and $k=26$.

Figure 5: Total expected profits under different values of $\beta$ and
$E_{M}$.

We observe that the storage fails to reach its maximum charging and
discharging power in the simulations, because if it were to sustain its
maximum charging power over 30 minutes, it would exceed its total energy
capacity. As such, we explore how the recorded profits would have evolved had
the storage resource had a different capacity. To this end, we repeat the
experiments by varying $E_{M}$ from $4.0$ MWh to $24.0$ MWh with $4.0$ MWh
increments. We note from Fig. 5 that, across most values of $\beta$, the total
expected profits rise under growing values of $E_{M}$, which can be ascribed
to an increased ability to leverage price differences under a larger storage
capacity. As $\beta$ increases, however, the profits seem to be plagued by
diminishing returns, and indeed plateau at a certain $E_{M}$ level as $\beta$
increases beyond $0.2$. These observations bring home that risk attitudes (via
increasing $\beta$ values) can create higher hurdles to notching greater
profits, which are not necessarily mitigated by a larger storage capacity.

Figure 6: Total expected profits under different values of $\alpha$ and
$\beta$.

In Fig. 6, we examine the influence of the risk confidence level $\alpha$ as
well as $\beta$ on the total expected profits over the $K$ time periods. Note
that, driving down $\alpha$ toward zero brings $\mathcal{R}_{\alpha}(\bm{x})$
toward the expected value of the charging costs, signifying a more risk-
neutral decision-maker. In contrast, increasing $\alpha$ toward one drives
$\mathcal{R}_{\alpha}(\bm{x})$ towards the highest value the costs can take,
representing a more conservative decision-maker. Fig. 6 bears out the
diminished capability of risk-averse decision-maker to respond to pre-dispatch
price signals and to arbitrage, manifesting itself through dwindling profits
under increasing values of $\beta$ and $\alpha$. These results reaffirm the
well-known relationship between risk and profit, whereby expected profits
climb under less conservative decisions, characterized by the values of
$\alpha$ and $\beta$ approaching zero.

##  V Conclusion and Policy Implications

We model the self-scheduling problem of a risk-averse storage resource in
decentralized markets. Four core results are gleaned from the analysis. First
that risk-aversion tends to lead to a myopic approach to dispatch most notably
evident in seeking to arbitrage between near-term intervals, rather than
taking charging positions in the expectation of more profitable but risky
revenues much later in the look-ahead period. Second that risk aversion tends
to lead to reduced profits because of the conservative operating stance
adopted by the resource. Third, that this conservatism can mean that the
resource forgoes opportunities to dispatch at peak spot prices, which are
times of highest system scarcity. Finally, and importantly, while higher
energy storage capacity can mitigate low profits, increasing the storage
capacity has virtually no impact on profitability for higher risk-aversion
levels.

There are important policy implications that flow from this study. First, SOs
need to be acutely aware of the role of risk aversion in storage resource
dispatch and bidding. In particular, higher levels of risk aversion may mean
that storage resources may not be available for dispatch at times where
scarcity price signals are most evident. As risk aversion itself is a non-
observable quantity, this introduces uncertainty into an SO’s short-term and
medium term reliability forecasts, and risk into its decision-making on market
intervention. This points to an increasing need for system transparency on key
storage parameters including state-of-charge so that SOs have a better view on
available storage capacity. Finally, it points to the need for the SO itself
to develop a set of tools that quantify such risks and allow for better risk-
adjusted decisions on intervention and other system operations during
scarcity.

## References

  * [1] O. Yurdakul and F. Billimoria. Online companion to risk-averse self-scheduling of storage in decentralized markets. [Online]. Available: https://github.com/oyurdakul/pesgm23
  * [2] A. Conejo, F. Nogales, J. Arroyo, and R. Garcia-Bertrand, “Risk-constrained self-scheduling of a thermal power producer,” _IEEE Transactions on Power Systems_ , vol. 19, no. 3, pp. 1569–1574, 2004.
  * [3] A. Papavasiliou, Y. He, and A. Svoboda, “Self-Commitment of Combined Cycle Units Under Electricity Price Uncertainty,” _IEEE Transactions on Power Systems_ , vol. 30, no. 4, pp. 1690–1701, 2015.
  * [4] R. A. Jabr, “Robust self-scheduling under price uncertainty using conditional value-at-risk,” _IEEE Transactions on Power Systems_ , vol. 20, no. 4, pp. 1852–1858, 2005.
  * [5] Y. Wang, Y. Dvorkin, R. Fernandez-Blanco, B. Xu, T. Qiu, and D. S. Kirschen, “Look-ahead bidding strategy for energy storage,” _IEEE Transactions on Sustainable Energy_ , vol. 8, no. 3, pp. 1106–1117, 2017.
  * [6] J. Kazempour, M. Moghaddam, M. Haghifam, and G. R. Yousefi, “Risk-constrained dynamic self-scheduling of a pumped-storage plant in the energy and ancillary service markets,” _Energy Conversion and Management_ , vol. 50, no. 5, pp. 1368–1375, 2009.
  * [7] Y. Chen and R. Baldick, “Battery storage formulation and impact on day ahead security constrained unit commitment,” _IEEE Transactions on Power Systems_ , 2022.
  * [8] R. T. Rockafellar and S. Uryasev, “Conditional value-at-risk for general loss distributions,” _Journal of banking & finance_, vol. 26, no. 7, pp. 1443–1471, 2002.
  * [9] G. C. Pflug, “Some remarks on the value-at-risk and the conditional value-at-risk,” in _Probabilistic constrained optimization_. Springer, 2000, pp. 272–281.
  * [10] AEMO. (2022) Market Data NEMWeb - Information Registry. [Online]. Available: https://https://aemo.com.au/en/energy-systems/electricity/national-electricity-market-nem/data-nem/market-data-nemweb


#  DSOL: A Fast Direct Sparse Odometry Scheme

Chao Qu, Shreyas S. Shivakumar, Ian D. Miller and Camillo J. Taylor We
acknowledge the support of Distributed and Collaborative Intelligent Systems
and Technology Collaborative Research Alliance. Ian Miller acknowledges the
support of a NASA Space Technology Research Fellowship.C. Qu, S. S.
Shivakumar, I. Miller and C. J. Taylor are with the GRASP Laboratory, School
of Engineering and Applied Sciences, University of Pennsylvania {quchao,
sshreyas, iandm, cjtaylor} @seas.upenn.edu

###### Abstract

In this paper, we describe Direct Sparse Odometry Lite (DSOL), an improved
version of Direct Sparse Odometry (DSO) [1]. We propose several algorithmic
and implementation enhancements which speed up computation by a significant
factor (on average 5x) even on resource constrained platforms. The increase in
speed allows us to process images at higher frame rates, which in turn
provides better results on rapid motions. Our open-source implementation is
available at https://github.com/versatran01/dsol.

##  I Introduction

Localization and mapping are key components of many robotic systems. In this
work we are motivated by the requirements of micro aerial vehicles where
payload and power limit the choice of sensors and the computational capacity
of the system. For these platforms, vision is an attractive sensing modality
since cameras are relatively inexpensive in terms of mass and power
consumption.

We choose to build our work on DSO and Stereo-DSO (SDSO) as described in [1,
2]. The direct approach to visual odometry minimizes the photometric error
instead of a geometric one. In their work Engel et al. provide extensive
experiments which show that the direct and sparse combination offers unique
advantages over state-of-the-art feature-based methods like ORB-SLAM [3] and
semi-dense direct methods such as LSD-SLAM [4] in terms of accuracy and speed.

The computational advantage is particularly important in the aerial context
that motivates this work, where agile robots often travel up to 10m/s [5]. We
show that the proposed approach can run at rates of up to 500Hz on common
computational hardware. For direct methods, the improvement in runtime also
provides concomitant improvements in accuracy and robustness [6]: faster
processing times allows for higher frame rates which can then estimate
moderate motions with greater fidelity and aggressive motions with higher
robustness.

For many robotics application, DSO is not suitable since it is a monocular
method which does not have a consistent, metric scale and has delayed
initialization and re-initialization. This poses great safety concerns,
especially for flying robots. For this reason, we design our system to be
directly initializable from stereo and/or depth images. In this work, we focus
on the stereo version and evaluate our system against SDSO, which is itself an
improvement upon DSO.

The remaining sections describe the proposed system in more details,
highlighting differences compared to DSO and SDSO, and explaining relevant
implementation choices that lead to improvements. We summarize our changes as
follow: (1) we adopt the inverse compositional alignment approach in frame
tracking, (2) we track the new image w.r.t the entire window instead of the
last keyframe, (3) we propose a better stereo photometric bundle adjustment
formulation compared to SDSO, (4) we greatly simplify the keyframe creation
and removal criteria from DSO, and (5) we parallelize our system to
effectively utilize all available computational resources. Together, these
changes lead to a simple and lightning-fast direct sparse odometry system,
which we call Direct Sparse Odometry Lite (DSOL).

Figure 1: Results of DSOL on selected sequences from the TartanAir dataset
[7]. Green lines are groundtruth trajectories, red lines are estimated
trajectories.

##  II Related Work

Visual odometry (VO) algorithms can be broadly categorized along the following
two axes: direct vs. indirect and dense vs. sparse. Direct methods recover
model parameters directly from images by minimizing photometric error based on
the brightness constancy assumption [8]. This is in stark contrast to indirect
methods, often called feature-based methods, where correspondences are first
established based on some intermediate representations, and the model
parameters are optimized by minimizing reprojection errors. Dense methods aim
to use all information from the image for better accuracy and robustness at
the cost of increased computation. Sparse methods, on the other hand,
recognize that image data is highly redundant and choose to only process a
selected yet informative subset of the image.

Early VO/V-SLAM systems were mostly sparse and indirect [9, 10, 11, 12]. This
was partly due to the computation limit at the time, but also dictated by the
needs of loop closure schemes in full-fledged SLAM systems, which often relied
on feature descriptors [13]. Among these systems ORB-SLAM [3] stood out as a
reference implementation of indirect approaches for its superior accuracy and
versatility.

Dense and direct methods started gaining traction with the advent of more
powerful compute systems. DTAM [14] was the first system to demonstrate real-
time camera tracking and dense scene reconstruction running on a GPU. Several
semi-dense systems [15, 4] further reduced the computational demands by only
working with high-gradient pixels.

SVO [16] adopts a hybrid approach starting with a direct image alignment
process to acquire an initial motion estimate before moving on to optimizing
reprojection errors of converged map points, hence the name semi-direct.

The direct and sparse combination that we are interested in was originally
introduced in [17] and later popularized by DSO [1]. It has then been extended
to work with fish-eye lenses [18], stereo [2] and rolling shutter cameras
[19], loop closure systems [20] and integration with deep networks [21].

This work does not aim to introduce an entirely new system, but rather to show
that with a number of optimizations and refinements, as well as thoughtful
implementation choices, we can dramatically improve the performance of a
direct sparse odometry system.

##  III Direct Methods for Visual Odometry

In this section, we outline the direct methods used in our system. Namely,
direct image alignment for frame tracking and photometric bundle adjustment
for joint optimization.

###  III-A Preliminary

Our system maintains a sliding window of $N$ keyframes
$\mathcal{K}=\\{K_{1},\cdots,K_{N}\\}$. Each keyframe $K_{i}$ is associated
with a Gaussian pyramid of $P$ images
$\mathcal{I}_{i}=\\{I^{0}_{i},\cdots,I^{P}_{i}\\}$, a set of affine brightness
parameters $\mathbf{a}_{i}=(a_{i},b_{i})^{\top}$, a camera pose
$\mathbf{T}_{i}^{\mathcal{W}}\in\mathbb{SE}(3)$ w.r.t. the world frame
$\mathcal{W}$, and a set of $m_{k}$ points parameterized by inverse depths
$\mathcal{P}_{i}=\\{\rho^{1}_{i},\cdots,\rho^{m}_{i}\\}$ hosted in this
keyframe. The total number of points in the sliding window is
$M=\sum_{\mathcal{K}}m_{k}$. Note that we need an additional image pyramid and
affine brightness parameters when stereo images are available, but only keep a
single set of poses and points per keyframe for the left camera.

We assume that the camera system has been calibrated before use and the images
undistorted such that each camera model can be described by an intrinsic
matrix $\mathbf{K}$. If a stereo rig is used we assume that the images are
rectified and that the extrinsic parameters are known. Any 3D point
$\mathbf{p}=(X,Y,Z)^{\top}$ in the camera frame can be mapped to pixel
coordinate $\mathbf{u}=(u,v)^{\top}$ via the projection (forward) function
$\Pi_{\mathbf{K}}:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}$, where

$\mathbf{u}=\Pi_{\mathbf{K}}(\mathbf{p})=\left(f_{x}\frac{X}{Z}+c_{x},\
f_{y}\frac{Y}{Z}+c_{y}\right)^{\top}.$ (1)

Similarly, given pixel coordinate $\mathbf{u}$ and its inverse depth $\rho$,
we can recover the 3D point $\mathbf{p}$ via the inverse projection (backward)
function $\Pi_{\mathbf{K}}^{-1}$, where

$\mathbf{p}=\Pi_{\mathbf{K}}^{-1}(\mathbf{u},\rho)=\left(\frac{u-c_{x}}{f_{x}\rho},\
\frac{v-c_{y}}{f_{y}\rho},\ \frac{1}{\rho}\right)^{\top}.$ (2)

When the context is clear, we simply refer to the inverse depth and its pixel
coordinate as a point $\mathbf{p}:=(\mathbf{u},\rho)$.

For photometric camera calibration, we adopt the image formation model from
[22], which consists of a non-linear response function and a lens attenuation
(vignetting) function. These parameters are acquired via the software provided
in [23]. Finally, we assume a well-calibrated camera and do not optimize any
calibration parameters online.

###  III-B Frame Tracking

Figure 2: Outlier rejection. In this visualization solid squares denote inlier
feature points, while hollow ones are outliers. The top image shows outliers
on moving objects. The bottom image shows outliers due to occlusion. Better
viewed in color and zoomed in.  Figure 3: Tracking a new image $I$ w.r.t. the
latest keyframe ($K_{1}$ only) vs. the entire window ($K_{0}$ and $K_{1}$).
The former can only track points in area 1 and 2, while the later can also
utilize points from area 3.

####  III-B1 Direct Image Alignment

Each new image $I_{j}$ is tracked w.r.t. a local map
$\mathcal{M}=\\{\mathcal{P}_{1},\cdots,\mathcal{P}_{N}\\}$ via direct image
alignment [24]. The parameters we wish to optimize are
$\mathbf{x}_{j}=(\mathbf{T}^{\mathcal{W}}_{j},\mathbf{a}_{j})$ while all
keyframe poses stay fixed. A simple constant velocity model provides an
initial guess to $\mathbf{T}^{\mathcal{W}}_{j}$ in the absence of an Inertial
Measurement Unit (IMU).

Following standard practice [1], we define the photometric error
$E_{\mathbf{p}}$ of a point $\mathbf{p}=(\mathbf{u},\rho)\in\mathcal{P}_{i}$
in a reference frame $I_{i}$ (a keyframe), observed in the target frame
$I_{j}$ (the new image), as the weighted sum of squared differences (SSD) over
a small patch of pixels $\mathcal{N}_{\mathbf{u}}$ . Let

$\displaystyle E_{\mathbf{p}}$
$\displaystyle=\sum_{\mathbf{u}_{k}\in\mathcal{N}_{\mathbf{u}}}w_{k}r_{k}^{2}(\mathbf{x}_{j})$
(3)
$\displaystyle=\sum_{\mathbf{u}_{k}\in\mathcal{N}_{\mathbf{u}}}w_{k}\left((I_{i}[\mathbf{u}_{k}]-b_{i})-\frac{e^{a_{i}}}{e^{a_{j}}}(I_{j}[\mathbf{u}^{\prime}_{k}]-b_{j})\right)^{2}$
(4)

where $r_{k}$ denotes the photometric residual associated with pixel $u_{k}$,
$I[\cdot]$ is the pixel index operator and $\mathbf{u}^{\prime}_{k}$ the
warped pixel of $\mathbf{u}_{k}$ from $I_{i}$ to $I_{j}$, given by the warping
function $W(\cdot,\cdot)$

$\mathbf{u}^{\prime}_{k}=W(\mathbf{p}_{k},\mathbf{T}_{i}^{j})=\Pi_{\mathbf{K}}(\mathbf{p}^{\prime}_{k})=\Pi_{\mathbf{K}}(\mathbf{T}_{i}^{j}\Pi_{\mathbf{K}}^{-1}(\mathbf{u}_{k},\rho_{k})).$
(5)

$w_{k}=w_{k}^{g}\cdot w_{k}^{r}$ is a combined weighting factor where

$\displaystyle w_{k}^{g}=\frac{c^{2}}{c^{2}+\|\nabla
I_{i}[\mathbf{u}_{k}]\|^{2}_{2}},\quad
w_{k}^{r}=\frac{\nu+1}{\nu+(\frac{r_{k}}{\sigma_{k}})^{2}}.$ (6)

$w_{k}^{g}$ penalizes pixels with high gradients [15] and $w_{k}^{r}$ is a
robust weight assuming residuals of the sparse model follow a t-distribution
with degree of freedom (DoF) $\nu$ [25, 26]. In addition, we employ a simple
gradient-based outlier rejection scheme to handle large errors due to
occlusions or moving objects, as shown in Fig. 2. An observation is considered
bad if $r_{k}^{2}>\|\nabla I_{i}[\mathbf{u}_{k}]\|^{2}_{2}$. If more than one
pixel within a patch is bad, the entire patch is discarded from the
optimization.

The total error for direct image alignment is

$E_{\mathrm{align}}=\sum_{I_{i}\in\mathcal{K}}\sum_{\mathbf{p}\in\mathcal{P}_{i}^{j}}E_{\mathbf{p}},$
(7)

where $\mathcal{P}_{i}^{j}$ is the set of points in $K_{i}$ that are visible
in $I_{j}$.

In contrast to DSO we adopt the inverse compositional alignment [24] method,
which swaps the roles of $I_{i}$ and $I_{j}$, allowing computationally
demanding steps to be performed ahead of time [27]. Our system also supports
tracking using both stereo images. This provides improved accuracy at the cost
of extra computation. It also improves robustness in situations where one of
the cameras is temporarily blocked: map points that are still visible in the
other image can create enough constraints to estimate model parameters.

####  III-B2 Frame-to-Window Alignment

In DSO, the new image is tracked w.r.t. the last keyframe. All active points
are projected into this keyframe. Frame parameters are then optimized using
conventional two-frame direct image alignment [28]. Although this worked well
in practice, we found it to be sub-optimal for the following reasons: (1)
Tracking w.r.t. a single image limits the number of points that can be used,
since points that are visible in the new image but not in the last keyframe
are discarded. (2) Every point will be projected twice: first into the last
keyframe, then into the new image. This incurs computational overhead and
prevents the reuse of previously computed gradients and Jacobians.

Instead, we propose tracking the new image w.r.t. the entire window. By using
all available keyframes, we essentially create a larger virtual field of view
(FOV), allowing more points to be used during alignment. This is illustrated
in Fig. 3, where tracking image $I$ w.r.t. both $K_{0}$ and $K_{1}$
incorporates more points into the optimization compared to using only $K_{1}$.
This improves robustness and reduces the number of keyframes created in cases
where the camera undergoes repetitive motion. Note that since all keyframe
poses are in the same reference frame, the inverse compositional formulation
still applies.

###  III-C Bundle Adjustment

####  III-C1 Photometric Bundle Adjustment

Upon adding a new keyframe, we perform windowed photometric bundle adjustment
(PBA) to update all model parameters.

$\displaystyle\mathbf{x}$
$\displaystyle=(\mathbf{x}_{1},\cdots,\mathbf{x}_{N})^{\top}$ (8)
$\displaystyle=(\mathbf{T}_{1},\mathbf{a}_{1},\mathcal{P}_{1},\cdots,\mathbf{T}_{N},\mathbf{a}_{N},\mathcal{P}_{N})^{\top}.$
(9)

We use the same cost function in Eqn. 3 with total error

$E_{\mathrm{pba}}=\sum_{I_{i}\in\mathcal{K}}\sum_{I_{j}\in\overline{\mathcal{K}}_{i}}\sum_{\mathbf{p}\in\mathcal{P}_{i}^{j}}E_{\mathbf{p}},$
(10)

where $\overline{\mathcal{K}}_{i}$ is the set of keyframes excluding $K_{i}$.
By forgoing explicit data association between images, PBA tries to project all
points from each keyframe to the others, forming a large number of implicit
photometric correspondences.

The factor graph representation of PBA is shown in Fig. 4(a). A directed edge
indicates the source of projection and an undirected one combines projections
from both keyframes. A window of size $N$ requires $N(N-1)$ pairwise
projections between keyframes. For example, if we have 5 keyframes, each
hosting 500 points with 5 pixels per patch, optimizing over 4 pyramid levels
with 4 iterations per level, there will be $5\times 4\times 500\times 5\times
4\times 4=8\times 10^{5}$ total residuals. This is several orders of magnitude
more compared to indirect methods [29], and poses significant challenges to
real-time operation. For this reason, DSO limits the total number of points to
800 in its fast setting. However, we show that it is possible to achieve even
a 100fps keyframe rate without limiting the number of points by leveraging
cache-friendly data structures and a parallel implementation, which we
describe in Sec. V.

(a) DSO [1]

(b) SDSO [2]

(c) DSOL (ours)

Figure 4:  Factor graphs of PBA with 3 keyframes. A directed edge in the graph
highlights the source of projection. For example, the edge
$\mathbf{x}_{0}^{L}\rightarrow\mathbf{x}_{0}^{R}$ represents the set of
factors formed by projecting $\mathcal{P}_{0}$ into $I_{0}^{R}$. An undirected
edge combines the pairwise projections between two keyframes.

####  III-C2 Stereo PBA

Using stereo images with a known baseline in PBA renders the scale of the
system directly observable and significantly reduces scale-drift compared to
its monocular counterpart. In SDSO [2], static stereo costs were added
alongside temporal ones from DSO, which constrain a subset of the total
parameters. The corresponding factor graph is shown in Fig. 4(b), where static
stereo factors are depicted by direct edges and temporal ones by undirected
ones. A carefully tuned stereo coupling factor was introduced to balance the
weight between these two types of cost.

In contrast, we choose to project every point into every other image (both
left and right) and from an even denser network of constraints, as shown in
Fig. 4(c). This eliminates the need for a coupling factor as all costs are now
temporal, with half of them depending on the fixed baseline. Note that this
also doubles the runtime of the algorithm. However, we believe that the time
penalty is more than justified by the improved robustness, direct
initialization, and scale recovery, which are indispensable to real-world
applications.

###  III-D Optimization and Marginalization

We optimize both $E_{\mathrm{align}}$ and $E_{\mathrm{pba}}$ using a few
iterations of the Gauss-Newton algorithm at each pyramid level. We stop early
if the optimization plateaus. The linearized system (at linearization point
$\mathbf{x}_{0}$) is given by

$\displaystyle\mathbf{H}\delta\mathbf{x}=\mathbf{b}\quad\mathrm{or}\quad\mathbf{J}^{\top}\mathbf{W}\mathbf{J}\delta\mathbf{x}=-\mathbf{J}^{\top}\mathbf{W}\mathbf{r},$
(11)

where $\mathbf{H}$ is the Hessian, $\mathbf{J}$ the Jacobian, $\mathbf{W}$ the
weight matrix, and $\mathbf{r}$ the residual vector. For direct image
alignment, the linear system is typically small because the variables are
parameters of the current frame $\mathbf{x}_{j}$. For PBA the dimension of
$\mathbf{H}$ scales linearly with the size of the sliding window and the
number of pixel features and is typically on the order of $10^{3}$. Although
large, this Hessian exhibits certain exploitable sparsity structures.
Specifically, we can divide the linear system into the following block form

$\displaystyle\begin{bmatrix}\mathbf{H}_{\mathrm{pp}}&\mathbf{H}_{\mathrm{pm}}\\\
\mathbf{H}_{\mathrm{mp}}&\mathbf{H}_{\mathrm{mm}}\end{bmatrix}\begin{bmatrix}\delta\mathbf{x}_{\mathrm{p}}\\\
\delta\mathbf{x}_{\mathrm{m}}\end{bmatrix}=\begin{bmatrix}\mathbf{b}_{\mathrm{p}}\\\
\mathbf{b}_{\mathrm{m}}\end{bmatrix},$ (12)

where subscript $\mathrm{p}$ denotes the set of pose and affine variables, and
$\mathrm{m}$ denotes the set of inverse depths. The bottom right block
$\mathbf{H}_{\mathrm{mm}}$ is a diagonal matrix because points are
parameterized by inverse depths. To solve the full system, we first
marginalize out all depth variables by applying the Schur complement and
arrive at the reduced system
$\mathbf{H}_{\mathrm{pp}}^{*}\delta\mathbf{x}_{\mathrm{p}}=\mathbf{b}_{\mathrm{p}}^{*}$,

$\displaystyle\mathbf{H}_{\mathrm{pp}}^{*}$
$\displaystyle=\mathbf{H}_{\mathrm{pp}}-\mathbf{H}_{\mathrm{pm}}\mathbf{H}_{\mathrm{mm}}^{-1}\mathbf{H}_{\mathrm{mp}}$
(13) $\displaystyle\mathbf{b}_{\mathrm{p}}^{*}$
$\displaystyle=\mathbf{b}_{\mathrm{p}}-\mathbf{H}_{\mathrm{pm}}\mathbf{H}_{\mathrm{mm}}^{-1}\mathbf{b}_{\mathrm{m}}$
(14)

$\delta\mathbf{x}_{\mathrm{p}}$ can be solved for efficiently and back-
substituted to yield

$\displaystyle\delta\mathbf{x}_{\mathrm{m}}=\mathbf{H}_{\mathrm{mm}}^{-1}(\mathbf{b}_{\mathrm{m}}-\mathbf{H}_{\mathrm{mp}}\delta\mathbf{x}_{\mathrm{p}})$
(15)

The new state is then updated via
$\mathbf{x}\leftarrow\mathbf{x}\boxplus\delta\mathbf{x}$, where $\boxplus$ is
the box-plus operator defined in [30].

When we need to remove variables from the active set (removing a keyframe
$K_{i}$ from the window), we compute a linear system $\mathbf{H}_{i}$, similar
to Eqn. 13, containing residuals that depend on $\mathbf{x}_{i}$ (variables to
be removed). We first marginalize all points from $\mathbf{H}_{i}$ and then
the frame variables. The result is a quadratic form on the rest of the
keyframe parameters $\hat{\mathbf{x}}_{i}$, which can be converted to a set of
linear constraints
$\hat{\mathbf{H}}_{i}\delta\hat{\mathbf{x}}_{i}=\hat{\mathbf{b}}_{i}$. This
can then be directly added to the normal equation of subsequent optimization
as a linear prior. Note that for all keyframes that have a prior term, we need
to fix their first estimates at $\overline{\mathbf{x}}$ and only update the
quadratic marginalization prior $\Delta\mathbf{x}$ to avoid inconsistencies
introduced by repeated linearizations [31]. The evaluation point can then be
recovered by
$\mathbf{x}=\overline{\mathbf{x}}\boxplus(\Delta\mathbf{x}+\delta\mathbf{x})$.
We refer the reader to [1, 32] for more details on marginalization.

##  IV Visual Odometry Front-end

Figure 5:  Keyframe initialization. Selected points with large local gradients
are in cyan. Colored square around each pixel reflects its initial inverse
depth. Depths are initialized from stereo matching (thick border) and local
map (thin border). Better viewed in color and zoomed.

In this section, we describe our VO front-end, which manages keyframes and
points within the window, including their creation, initialization and
removal.

###  IV-A Keyframe Creation

DSO employs a complicated set of rules for keyframe creation, which involves
computing the mean squared optical flow, the mean flow without rotation and
the change in affine brightness factor, all from the current image to the last
keyframe [1]. The weighted sum of these quantities is then compared to another
threshold to determine whether a new keyframe is needed, resulting in 4 extra
parameters to tune. In addition, DSO follows the philosophy as ORB-SLAM [3,
29], which is to take many keyframes upfront and sparsify them with early
marginalization. However, this strategy is ineffective when computational
resources are limited.

Keyframe creation is arguably the most time-consuming part of the entire
pipeline (often 5-10x slower than tracking, see Table I). Upon adding a new
keyframe, one needs to select candidate pixels, precompute patch gradients,
initialize point depths, and perform PBA and other necessary bookkeeping.
Creating too many keyframes (even in a separate thread) will cause the system
to eventually lag behind the image frame rate.

We argue that a keyframe is only needed when the current image cannot be
reliably tracked w.r.t. the sliding window. If enough points from the local
map can be successfully projected into the image, then we can simply keep
using these keyframes. This prevents the addition of new ones that do not
contribute very much to frame tracking. Quantitatively, we define the tracking
ratio $Q$ as the quotient between the number of tracked points and the number
of selected points from all keyframes in the window. We create a new keyframe
if $Q$ falls below $Q_{\min}$. A higher $Q_{\min}$ therefore results in more
keyframes taken.

###  IV-B Keyframe Initialization

When a new keyframe is created, we select points with high local-gradients,
extract a small patch around it and initialize its inverse depth if possible,
followed a final PBA (Sec. III-D). Once a keyframe is added it is immediately
used for frame tracking (Sec. III-B).

####  IV-B1 Point Selection

We divide the image into a grid of cells with size 16$\times$16 pixels and
select at most one point from each cell. We start by finding one candidate
pixel per cell with a large image gradient. We then choose candidates with
gradient greater than $g_{\min}$. If this first round of selection results in
too few points, we do a second round with $g_{\min}/2$. However, $g_{\min}$ is
a running quantity, so we also decrease it by $\delta_{g}$ for the next image.
On the other hand, if too many are chosen, we simply increase $g_{\min}$ by
$\delta_{g}$. In this way $g_{\min}$ is constantly adjusted to provide the
desired number of points. Hence $g_{\min}$ is continuously adapted to select a
sufficient number of points with a good spatial distribution and is upper-
bounded by the total cells in the grid.

Unlike DSO, we avoid selecting points from image areas where previous map
points were found. This is to prevent creating duplicate points that
correspond to the same scene surface area, which negatively impact the
linearized system by introducing unwanted dependencies [33]. To this end, we
project all map points into the current image and slightly dilate them to
create a binary mask. Only pixels outside of this mask will participate in the
selection process. Fig. 5 shows selected pixels for sample images.

####  IV-B2 Patch Extraction

We extract a small patch around each selected point. We adopt a simple cross-
shaped patch with a total of 5 pixels (see pattern 1 in [1]). This is done at
every pyramid level by down-scaling the pixel coordinate accordingly.
Intensity values and gradients are computed for each pixel within the patch
using bilinear interpolation. These quantities are re-used during tracking to
reduce computation.

###  IV-C Point Depth Initialization

Points in a new keyframe have their inverse depths initialized in the
following order if available: depth image, stereo images, map points. Any
point depth set by a previous method will not be considered by latter ones.
Note that we do not perform small baseline stereo matching to initialize
depth, as was done in DSO [1]. We also do not distinguish active points from
candidate ones: all points with valid depths are active and can be used in
tracking and PBA.

####  IV-C1 From depth image

If an aligned depth image is available (either from a depth camera, a lidar
projection, or even a depth prediction network), we directly set the inverse
depth of any pixel with a valid measurement.

####  IV-C2 From stereo images

If a pair of rectified stereo images is available, a simple sparse stereo
matching scheme is employed on pixels that are yet to be initialized. Since
the scale pyramid of both images are already constructed, we carry out a
correspondence search in a coarse-to-fine manner to reduce the number of
comparisons. In particular, we only search the max disparity range at the
lowest pyramid level, and gradually refine the disparity estimate while going
up. Similar to SDSO [2], we use zero-normalized cross correlation (ZNCC) as
the similarity metric on a small patch of pixels ($5\times 7$). Note that our
goal here is to initialize as many points as possible, therefore we do not
enforce any strict outlier rejection scheme, but rely on PBA to remove points
with large errors.

####  IV-C3 From map points

Finally, if there are points remain uninitialized, we attempt to estimate
their depths from the local map. The direct image alignment step generates a
set of valid projections from all map points onto the current image, which can
be used for this purpose. We average the inverse depths of all points
(observed from the new image) that fall into a single cell and use it to
initialize the new point selected from this cell. Again, accuracy is not
important here, since they will be refined or removed by PBA. Fig. 5 also
shows the depth initialization process for some sample images.

###  IV-D Keyframe Removal

To bound computational complexity, the window size is fixed at $N$. DSO
proposed several heuristics for keyframe removal, which could remove multiple
keyframes at the same time. Specifically, they remove keyframes where less
than $5\%$ of their points are visible in the most recent keyframe. This was
partly due to their tracking strategy, which requires projection of active
points into the latest keyframe (Sec. III-B). Thus it is natural to remove
keyframes that do not have enough overlap with the last keyframe.

We argue that this is sub-optimal, despite having a slight improvement in
speed (due to fewer keyframes in the window). Since points anchored in a
particular keyframe are projected into every other keyframe in PBA, removing a
keyframe weakens the links among the network of constraints formed by all
pairwise projections. Although proper marginalization does retain a linear
prior connecting the remaining keyframes, it is weaker compared to the
original nonlinear constraints. Therefore we choose to remove only one
keyframe at a time with the fewest points visible in the current image. By
always having a full window, we ensure strong connections among keyframes.
Finally, before actually discarding the keyframe, we perform keyframe
marginalization according to Sec. III-D.

##  V Implementation Details

Figure 6:  Failure cases of our method on the TartanAir dataset [7]. Causes of
failure from left to right are: repetitive texture, blocked camera and low
texture.

We adopt the data-oriented design philosophy used in [34], where we leverage
array-based data structures for better cache-locality, ease of parallelism and
zero dynamic allocation at runtime. Each major stage in our pipeline can
parallelize its computation, including point selection, stereo matching, frame
tracking, and bundle adjustment. We apply a task-based, fine-grained
parallelism, where each task contains a single row of cells defined by the
point selection process. We found that this granularity struck a good balance
between context-switching overhead and resource utilization.

For direct image alignment, each task accumulates its own Hessian and the
final result is reduced by a fork-join model [35]. For PBA, due to the
$O(N^{2})$ complexity induced by pairwise projections, we add an extra
keyframe-level parallelization. Specifically, we spawn tasks to accumulate
block Hessians as long as they do not depend on the same set of points. The
final Hessian is updated by $N(N-1)$ block Hessians, which requires locking to
prevent data races. However, since $N$ is typically small, this mutex has a
low contention rate and the performance impact is negligible.

Our implementation follows modern C++ design principles [36], with no manual
resource or thread management. Under default settings, our system allocate
less than 10MB of memory, which fits into the L3 cache of most modern
processors. The small memory footprint is crucial to deployment on resource
constrained platforms.

We use a window of $N=4$ keyframes and a scale pyramid of $P=5$ levels. We
found that further increasing $N$ has diminishing returns in accuracy since we
already form a set of strongly-connected keyframes. We benchmark against an
open-source implementation of SDSO 111https://github.com/RonaldSun/stereo_DSO,
which is based on the official DSO with stereo integration.
Code222https://github.com/ShreyasSkandanS/stereo_DSO and
data333https://tinyurl.com/2z3nzyre to reproduce our results are also publicly
available.

##  VI Evaluation Results

###  VI-A Datasets

####  VI-A1 KITTI Odometry

The KITTI Odometry dataset [37] is an urban driving dataset with stereo
sequences and associated ground truth trajectories. We evaluate on all 11
(00-10) training sequences.

####  VI-A2 Virtual KITTI-2

The VKITTI2 dataset [38] is a synthetic video dataset intended to appear in
the same environment as KITTI with simulated lighting and weather conditions.
We use all 5 sequences, each with 10 different variations (including
challenging ones like fog and rain).

####  VI-A3 Tartan Air

TartanAir [7] is a synthetic dataset for robot navigation tasks collected in
photo-realistic simulation environments with various styles and scenes. It
contains diverse 6DoF motion patterns, making it more challenging compared to
driving datasets. We hand-picked 6 scenes (carwelding, gascola, oldtown,
office, office2, hospital) from the dataset, comprised of a total of 73
sequences covering both indoor, outdoor, static, and dynamic environments. We
only use the Easy variation of each sequence.

Figure 7:  Cumulative error plots of SDSO vs. DSOL on all runs (268 total).
y-axis shows percentage of runs. For an error value $x$, the corresponding $y$
value is the percentage of runs that has smaller error.  Figure 8:  Error
heatmaps of Trans- and Rot-APE for all 268 runs. Each pixel represents one
metric of a single run by one method. Trans- and Rot-APE are capped at $10\%$
and $25^{\circ}$ respectively, which typically indicate gross tracking
failure. The order of runs from top left to bottom right is VKITTI, KITTI,
TartanAir. This figure is best viewed in color

###  VI-B Accuracy and Robustness

(a) KITTI 00

(b) KITTI 05

(c) gascola 004

(d) office 004

Figure 9:  Qualitative results of SDSO and DSOL on selected sequences from
KITTI [37] and TartanAir [7].

To evaluate odometry accuracy, we run SDSO and DSOL on VKITTI2, KITTI and
TartanAir, both forward and backward, resulting in a total of
$(50+11+73)\times 2=268$ runs with a wide variety of scenarios. We report both
translational and rotational Absolute Pose Error (APE) and Relative Pose Error
(RPE), computed by the evo444https://github.com/MichaelGrupp/evo package. APE
measures the absolute pose error between the reference and the estimated
trajectory, whereas RPE measures the relative pose error between two
consecutive frames. We use Umeyama alignment [39] as a pre-processing step to
align the estimation to groundtruth. Note that due to the vast scale
differences between different sequences, we normalize the translational
metrics by the corresponding sequence length to obtain percentage errors. Both
methods are allowed to re-initialize when tracking is lost. However, if a
method cannot track the full sequence (due to the program crashing), we
consider it a failure only if it tracks no more than $80\%$ of the run and
assign it a fixed large RMSE.

Fig. 7 shows the cumulative error plots of SDSO and DSOL on all runs. A point
$(x,y)$ on a curve indicates that $y\%$ of runs has error less than $x$.
Therefore a curve closer to the top-left corner has better overall accuracy.
We see that DSOL outperforms SDSO by a large margin, as evidenced by the area
between the two curves.

We also show error heatmaps of all runs in Fig. 8, where each pixel represents
the corresponding metric from one run. The runs are in the order of VKITTI,
KITTI, and TartanAir from top left to bottom right within each plot. We see
that both methods perform relatively well on KITTI and VKITTI, where camera
motions are smooth and planar. However, SDSO struggles when facing the diverse
6DoF motion patterns from TartanAir, resulting in frequent tracking failures.
Our method was able to track most of the sequences from TartanAir with good
accuracy and without re-initialization, indicating a higher level of
robustness. However, there are still several failure cases that are generally
impossible for a vision-based odometry system, which we visualize in Fig. 6.
Qualitative results on selected sequences are shown in Fig. 9.

###  VI-C Speed

TABLE I: Runtime (in ms) of SDSO and DSOL on different processors under single and multi-threaded settings. Results are averaged among several testing sequences from TartanAir (with stereo images at $640\times 480$ resolution). S denotes single-thread mode, M-$n$ denotes multi-thread mode with at most $n$ threads.  Runtime [ms] | E5-2683 | i7-11800H | ARMv8.2
---|---|---|---
name-stage | S | M-32 | S | M-16 | S | M-6
SDSO-track | 13.95 | 15.02 | 7.18 | 7.29 | 34.39 | 32.73
DSOL-track | 4.81 | 2.27 | 2.39 | 1.23 | 9.35 | 5.08
Speed up | 2.90x | 6.67x | 3.00x | 5.93x | 3.68x | 6.44x
SDSO-kf | 79.07 | 52.64 | 71.59 | 40.56 | 282.6 | 195.5
DSOL-kf | 46.77 | 7.82 | 25.45 | 9.96 | 113.4 | 41.66
Speed up | 1.69x | 6.73x | 2.81x | 4.07x | 2.49x | 4.69x

We benchmark SDSO and DSOL on various compute units with single- and multi-
threaded settings. Processor types includes an Intel Xeon E5-2683
(16$\times$3.0GHz), an Intel i7-11800H (8$\times$4.6GHz) and an Nvidia Xavier
NX ARMv8.2 (6$\times$1.4GHz), which covers a wide range of core counts and
clock rates. We run both systems with stereo inputs on the gascola scene from
the TartanAir dataset and average timing results over several runs. We use the
default settings for SDSO (with 7 keyframes and 2000 points max) and do not
enforce real-time execution (no skipping frames). Note that it is difficult to
ensure a completely fair comparison, as each system uses slightly different
window sizes, pyramid levels, number of iterations and other hyper-parameters
that may affect its performance.

Runtime results are shown in Table I. We divide each system into two stages:
track is performed at the arrival of each image, while kf time is only
recorded when adding a new keyframe. DSOL runs faster than SDSO in all
benchmarks we performed. It is on average around 4.8x faster in frame tracking
and 3.8x faster in keyframe creation. Our system can track frames at more than
500 frames per second (FPS) and even handle keyframes at 100 FPS on laptop-
grade hardware. It is also worth noting that SDSO has relatively poor parallel
scalability (by dividing numbers under M by S in Table I), whereas ours can
often achieve further speedup where more computation resources are available.

###  VI-D High FPS Tracking

Figure 10: Qualitative results of DSOL on our Realsense D455 dataset with 60Hz
stereo images at 640$\times$480 resolution.

As a final experiment, we collected datasets to showcase tracking at high
frame rates with the Intel Realsense D455 sensor. We disabled the IR emitter
on the D455 and collected stereo images from the two grayscale cameras at
60FPS with an image resolution of 640$\times$480\. Qualitative results
obtained on two hand-held sequences are show in Fig. 10. Our system is able to
process the 60Hz stream in real time without dropping any frames.

##  VII Conclusion

In this paper, we have presented a fast, accurate and robust direct sparse
odometry (DSOL) that is suitable for robotics applications on resource
constrained platforms. Our system builds on ideas from DSO and SDSO with
several algorithmic and implementation improvements that enable it to handle
high frame rates. The significant speed boost was also made possible by
designing with parallelization in mind, allowing efficient usage of multi-core
processors. Our experiments have shown that our stereo version out-performed
the baseline SDSO on challenging datasets in terms of accuracy and speed. We
open source our implementation as a contribution to the community.

## References

  * [1] J. Engel, V. Koltun, and D. Cremers, “Direct sparse odometry,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 40, pp. 611–625, 2018.
  * [2] R. Wang, M. Schwörer, and D. Cremers, “Stereo dso: Large-scale direct sparse visual odometry with stereo cameras,” _2017 IEEE International Conference on Computer Vision (ICCV)_ , pp. 3923–3931, 2017.
  * [3] R. Mur-Artal, J. M. M. Montiel, and J. D. Tardós, “Orb-slam: A versatile and accurate monocular slam system,” _IEEE Transactions on Robotics_ , vol. 31, pp. 1147–1163, 2015.
  * [4] J. Engel, T. Schöps, and D. Cremers, “Lsd-slam: Large-scale direct monocular slam,” in _ECCV_ , 2014.
  * [5] K. Mohta, M. Watterson, Y. Mulgaonkar, S. Liu, C. Qu, A. Makineni, K. Saulnier, K. Sun, A. Z. Zhu, J. A. Delmerico, K. Karydis, N. A. Atanasov, G. Loianno, D. Scaramuzza, K. Daniilidis, C. J. Taylor, and V. R. Kumar, “Fast, autonomous flight in gps-denied and cluttered environments,” _ArXiv_ , vol. abs/1712.02052, 2018.
  * [6] A. Handa, R. A. Newcombe, A. Angeli, and A. J. Davison, “Real-time camera tracking: When is high frame-rate best?” in _ECCV_ , 2012.
  * [7] W. Wang, D. Zhu, X. Wang, Y. Hu, Y. Qiu, C. Wang, Y. Hu, A. Kapoor, and S. A. Scherer, “Tartanair: A dataset to push the limits of visual slam,” _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 4909–4916, 2020.
  * [8] M. Irani and P. Anandan, “All about direct methods,” 1999.
  * [9] A. J. Davison, I. D. Reid, N. Molton, and O. Stasse, “Monoslam: Real-time single camera slam,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 29, pp. 1052–1067, 2007.
  * [10] G. S. W. Klein and D. W. Murray, “Parallel tracking and mapping for small ar workspaces,” _2007 6th IEEE and ACM International Symposium on Mixed and Augmented Reality_ , pp. 225–234, 2007.
  * [11] H. M. Strasdat, J. M. M. Montiel, and A. J. Davison, “Scale drift-aware large scale monocular slam,” in _Robotics: Science and Systems_ , 2010.
  * [12] H. M. Strasdat, A. J. Davison, J. M. M. Montiel, and K. Konolige, “Double window optimisation for constant time visual slam,” _2011 International Conference on Computer Vision_ , pp. 2352–2359, 2011.
  * [13] M. J. Cummins and P. Newman, “Fab-map: Probabilistic localization and mapping in the space of appearance,” _The International Journal of Robotics Research_ , vol. 27, pp. 647 – 665, 2008.
  * [14] R. A. Newcombe, S. Lovegrove, and A. J. Davison, “Dtam: Dense tracking and mapping in real-time,” _2011 International Conference on Computer Vision_ , pp. 2320–2327, 2011.
  * [15] J. Engel, J. Sturm, and D. Cremers, “Semi-dense visual odometry for a monocular camera,” _2013 IEEE International Conference on Computer Vision_ , pp. 1449–1456, 2013.
  * [16] C. Forster, M. Pizzoli, and D. Scaramuzza, “Svo: Fast semi-direct monocular visual odometry,” _2014 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 15–22, 2014.
  * [17] H. Jin, P. Favaro, and S. Soatto, “A semi-direct approach to structure from motion,” _The Visual Computer_ , vol. 19, pp. 377–394, 2003.
  * [18] H. Matsuki, L. von Stumberg, V. C. Usenko, J. Stückler, and D. Cremers, “Omnidirectional dso: Direct sparse odometry with fisheye cameras,” _IEEE Robotics and Automation Letters_ , vol. 3, pp. 3693–3700, 2018.
  * [19] D. Schubert, N. Demmel, V. C. Usenko, J. Stückler, and D. Cremers, “Direct sparse odometry with rolling shutter,” in _ECCV_ , 2018.
  * [20] X. Gao, R. Wang, N. Demmel, and D. Cremers, “Ldso: Direct sparse odometry with loop closure,” _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 2198–2204, 2018.
  * [21] N. Yang, L. von Stumberg, R. Wang, and D. Cremers, “D3vo: Deep depth, deep pose and deep uncertainty for monocular visual odometry,” _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 1278–1289, 2020.
  * [22] J. Engel, V. C. Usenko, and D. Cremers, “A photometrically calibrated benchmark for monocular visual odometry,” _ArXiv_ , vol. abs/1607.02555, 2016\.
  * [23] P. Bergmann, R. Wang, and D. Cremers, “Online photometric calibration of auto exposure video for realtime visual odometry and slam,” _IEEE Robotics and Automation Letters_ , vol. 3, pp. 627–634, 2018.
  * [24] S. Baker and I. Matthews, “Lucas-kanade 20 years on: A unifying framework,” _International Journal of Computer Vision_ , vol. 56, pp. 221–255, 2004.
  * [25] C. Kerl, J. Sturm, and D. Cremers, “Robust odometry estimation for rgb-d cameras,” _2013 IEEE International Conference on Robotics and Automation_ , pp. 3748–3754, 2013.
  * [26] J. A. Zubizarreta, I. Aguinaga, and J. M. M. Montiel, “Direct sparse mapping,” _IEEE Transactions on Robotics_ , vol. 36, pp. 1363–1370, 2020\.
  * [27] S. Klose, P. Heise, and A. Knoll, “Efficient compositional approaches for real-time robust direct visual odometry from rgb-d data,” _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pp. 1100–1106, 2013.
  * [28] B. D. Lucas and T. Kanade, “An iterative image registration technique with an application to stereo vision,” in _IJCAI_ , 1981.
  * [29] 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, pp. 1255–1262, 2017.
  * [30] C. Hertzberg and G. Piazzi, “A framework for sparse, non-linear least squares problems on manifolds,” 2008.
  * [31] M. Li and A. I. Mourikis, “Consistency of ekf-based visual-inertial odometry,” 2012.
  * [32] S. Leutenegger, S. Lynen, M. Bosse, R. Y. Siegwart, and P. T. Furgale, “Keyframe-based visual–inertial odometry using nonlinear optimization,” _The International Journal of Robotics Research_ , vol. 34, pp. 314 – 334, 2015.
  * [33] H. Alismail, B. Browning, and S. Lucey, “Photometric bundle adjustment for vision-based slam,” _ArXiv_ , vol. abs/1608.02026, 2016.
  * [34] C. Qu, S. S. Shivakumar, W. Liu, and C. J. Taylor, “Llol: Low-latency odometry for spinning lidars,” _ArXiv_ , vol. abs/2110.01725, 2021.
  * [35] M. Voss, R. Asenjo, and J. Reinders, _Pro TBB: C++ Parallel Programming with Threading Building Blocks_ , 1st ed. USA: Apress, 2019.
  * [36] B. Stroustrup, _A Tour of C++_ , ser. Always learning. Addison-Wesley, 2014. [Online]. Available: https://books.google.com/books?id=-VNsAQAAQBAJ
  * [37] A. Geiger, P. Lenz, C. Stiller, and R. Urtasun, “Vision meets robotics: The kitti dataset,” _The International Journal of Robotics Research_ , vol. 32, pp. 1231 – 1237, 2013.
  * [38] Y. Cabon, N. Murray, and M. Humenberger, “Virtual kitti 2,” _ArXiv_ , vol. abs/2001.10773, 2020.
  * [39] S. Umeyama, “Least-squares estimation of transformation parameters between two point patterns,” _IEEE Transactions on Pattern Analysis & Machine Intelligence_, vol. 13, no. 04, pp. 376–380, 1991.