# Sizing Co-located Storage for Uncertain Renewable Energy Sold Through
Forward Contracts

Tomas Valencia Zuluaga and Shmuel S. Oren  Industrial Engineering & Operations
Research
University of California at Berkeley
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we propose a high-level Stochastic steady-state model to
analyze the value of co-located energy storage systems for wind power
producers that participate in an electricity market through Forward or Day
Ahead contracts. In particular, we try to find optimal sizing and contracting
and stationary operating policies for profit maximization in the long-run. We
obtain a stylized model calibrated to actual wind power production and
electricity wholesale price data that allows us to asses the value of storage
size and perform sensitivity analysis on key parameters such as contract
prices, storage cost and storage efficiency.

###### Index Terms:

Battery Storage, Forward Contracts, Optimal Sizing, Wind Generation

##  I Introduction

The impending consequences of climate change has driven massive proliferation
of renewable energy resources (RES) around the world. However, a key obstacle
to large-scale integration of RES in power systems is the short-term
uncertainty and variability of their output. This poses both technical
challenges for the reliable operation of the power system and financial
challenges for investors in RES, since it is hard to guarantee a reliable
income flow without RES implicit subsidies through policies such as feed-in
tariffs and mandatory contracting for retailers. Such policies however are not
sustainable and are being challenged due to their cost to consumers and many
systems are moving toward imposing scheduling requirements and forward
commitments on RES forcing them to compete on a level playing field with other
resources. The addition of Battery Energy Storage Systems (BESS) to renewable
power plants can help mitigate their uncertainty and thus can help towards
solving both the technical and financial issues. While it is intuitive to
understand how storage can mitigate uncertainty, sizing and managing BESS is
not obvious.

Optimizing the sizing and operation of energy storage for intermittent power
plants has been a very active research topic for several years, and there are
numerous references of relevant related work, both theoretical as well as in
applied cases. In [1], a theoretical analysis is performed for the case of a
wind farm that uses storage to optimize its bidding strategy on the day-ahead
market in order to minimize imbalance penalties. In this, and other works that
study the behavior of wind power producers (WPP) in day-ahead markets [2, 3],
it is observed that the optimal bidding strategy for WPP takes the form of an
optimal fractile, in the spirit of a Newsvendor problem. In fact, the bidding
problem of a WPP is a reverse newsvendor problem, where the uncertainty lies
not on the supply side.

There are also numerous references that address the issue of optimal sizing
and optimal management of co-located energy storage. In [4], the infinite-
horizon average cost of electricity purchases is minimized by finding an
optimal storage management policy and optimal storage size for a power plant
that serves a local demand and purchases any shortfall from the grid in
presence of dynamic pricing. The optimality of the balancing storage
management policy and the authors prove optimality of a dual threshold policy,
in the spirit of optimal $(s,S)$ policies for inventory. In [5], the optimal
size of a BESS is found for a grid-connected PV system that can purchase and
sell energy from the grid under time-of-use pricing, and the convexity of
profit in storage size is shown. In [6], the value of co-located storage is
analyzed for the case of the UK market as subsidies are phased out; a
comparison of participation in the wholesale market with co-located storage vs
a long-term power purchase agreement with contracts for differences is made
with help of a model using stochastic differential equations describing wind
power production and prices.

Despite the abundance of literature in the matter, the application of high-
scale, steady state models, such as fluid queues, to energy storage remains
relatively rare. Limiting distributions for fluid queue models have been most
commonly used in high-speed communication networks, but also for manufacturing
[7], and energy systems, with early applications in hydro dam management
models [8]. In [9], a two-state markov-modulated fluid queue (MMFQ) is used
for optimizing the size and management policy of the national strategic
petroleum reserve (SPR).

More recently, there have been application of fluid queue models to renewable
power plants with BESS. In [10], a model-predictive control (MPC) algorithm is
proposed for solving a joint storage configuration problem, i.e. sizing and
managing storage, and the MMFQ framework is used to analyze the performance of
the algorithm in terms of reliability of power availability. In [11], the
authors consider a case in which there is no control over the charge/discharge
rate (akin to a balancing policy), and use MMFQ to find a relation, in
particular an asymptotic relation, between the battery size and the loss of
load probability for a given grid configuration, so that the problem of
finding the optimal size of storage to hit a target loss-of-load probability
can be solved. Unlike in these cases, here we use the MMFQ framework in the
objective function of our optimization model directly.

While short term operation of co-located storage can be optimized using
dynamic programming techniques and accounting for current information on state
of charge, wind forecast and prices, the optimal sizing of the co-located
storage is based on long run average behavior of the production-storage system
under an optimized stationary policy. For this purpose we employ the spectral
method for characterizing limiting distributions of $n$-dimensional MMFQ
described in [7, 12]. Algorithms that are more numerically stable have been
developed [13, 14], but with a limitation of a single state with negative
drift. However, this limitation is too restrictive for our setting since it
would correspond to a single state where batteries can be discharged. In [15],
an algorithm is proposed for calculating limiting distributions of markov-
modulated fluid queues without solving an eigenvalue problem. The algorithm
has promising results in terms of stability to obtain results of a desired
precision and makes no assumptions regarding the structure of the fluid queue,
which is a promising option for general application models like ours.

In this paper, we use limiting distributions directly to obtain simple
expressions of long-run profits of the system, suitable for use in a
sensitivity analysis. In particular, we are interested in investigating the
effect of some key parameters such as cost of storage and available contract
prices on the financial outcome of the project, the optimal storage size and
the optimal contracting strategy of a wind power plant trying to maximize its
profit. This simple model would be of interest to a wind power plant operator
and project developer for early-stage project feasibility analysis. It could
also be of interest for policy designers, to evaluate the financial viability
of storage projects for wind power plants with and without special incentives,
e.g. by considering tentative subsidies in the form of reduced feed-in-tariffs
for surplus generation conditional on participation in long-term fixed
quantity contracts.

##  II Problem setting

We consider the case of a WPP that participates in the wholesale electricity
market and has access to a long-term forward market. The WPP is also
evaluating the construction of a collocated Battery Energy Storage System
(BESS), which will be used to mitigate imbalance costs. In our setting, the
wind farm operator needs to answer two main questions: how much storage should
be installed and how much energy should be sold in a long-term contract. To
answer these questions, it is also relevant to determine a stationary
management (charge/discharge) policy for the BESS. Further aspects of our
setting are described in more detail in the paragraphs below.

We consider that the WPP is a price-taker in the long-term forward market,
where it has access to three contracts:

  1. 1.

A long-term fixed quantity contract, where it can choose how much to sell for
a fixed price in every period of the wholesale electricity market.

  2. 2.

A pay-as-demanded forward contract, which can be used to cover any shortfall
in generation with respect to the pledge in contract 1.

  3. 3.

A pay-as-produced forward contract, which can be used to sell any excess
generation with respect to the pledge in contract 1.

In order to avoid arbitrage opportunities, contract 2 must have a higher price
than 1 and contract 3 must have a lower price than contract 1.

This setting can be seen as a special case of a WPP that participates in the
day-ahead market as described next. The WPP can commit to generate any
quantity up to its capacity. The producer has no control over the wind power
output and the commitment is made before the actual output is known. An
imbalance penalty is paid for the difference between the energy commitment and
the actual energy delivered. The case studied here is the special case where
the commitment is the same for all periods, the imbalance penalty for
shortfall is non-negative and constant, and the imbalance penalty for surplus
is non-positive (i.e. it is a reward, not a penalty) and constant.

To mitigate the imbalance cost, the wind power producer can charge/discharge
the BESS, so that the net output of the plant, i.e. that perceived by the
market is the combined production of the wind farm and the BESS. Having a co-
located BESS will have the effect of firming the WPP’s energy output, and thus
giving it access to other markets (e.g. a capacity market). However, we don’t
take any other sources of revenue into consideration in this setting.

The scope of our work is a high-level analysis that could be of interest for
early-stage project developers or policy makers. We are interested in looking
at infinite-horizon average profits by considering the steady-state behavior
of the model proposed and performing sensitivity analyses on a number of key
parameters. To be consistent with this approach, we do not consider the
possibility of using storage for arbitrage, since arbitrage opportunities
cannot exist in steady state in an efficient market. This approach suggests
the assumption of existence of a long-run equilibrium in the forward market,
which depends on many factors that are not considered in our model, such as
long-run demand behavior and a stabilization of the cost of new capacity
(which in the case of storage, for instance, is actually expected to continue
declining for some years [16]). We do not make any such strong assumptions.
While such an equilibrium may be far from being reached, this model is a high-
level analysis for which this coarse approximation of reality is sufficient.

We do not specify any particular battery technology. We include in our
analysis the issues of charge/discharge conversion efficiency and energy
dissipation. However, we do not include the degradation of the battery because
of usage and aging in our model. It has been found that this can be an
important characteristic to take into account in sizing studies [6], so this
could be an interesting feature to add to future versions of our model, but is
considered out of the scope of this paper.

##  III Model

###  III-A General description

As mentioned before, the case we study here is a special case of a day-ahead
model. To emphasize the flexibility of our model and lay the groundwork for a
future extension considering variable prices, we describe first the model of
the day-ahead market allowing for the possibility of variable prices, and then
focus on the special case addressed in this paper.

####  III-A1 The market

We consider a day-ahead market in which the WPP is a price taker, so it
commits to produce an amount $q_{t}$ (in MWh) during market period $t$ at a
price $p_{t}$ (in $/MWh), which is known at the time the commitment is made.
During period $t$, the actual wind power output of the farm is $w_{t}$ (in
MWh), while $r_{t}$ is the amount energy injected into the BESS ($r_{t}<0$ if
the energy is extracted), so that the net output of the plant is
$w_{t}-r_{t}$. The imbalance is $y_{t}=q_{t}-(w_{t}-r_{t})$ and the imbalance
penalty charged to the WPP is $\Xi_{t}=\Xi(y_{t},p_{t})$. We assume
$\Xi(\cdot)$ to be a known, deterministic, time-invariant function of the
imbalance and the energy price.

For the case with constant prices, i.e. the case of forward contracts, $p_{t}$
is constant and the imbalance is defined by (1). Note that the negative sign
in front of $\kappa^{\prime}$ implies that the WPP is not penalized for excess
injection, but, on the contrary, sells it on the forward market. To avoid
arbitrage opportunities, we must have $\kappa^{\prime}<1$ and $\kappa>1$.

$\Xi(y_{t},p_{t})=\begin{cases}\Xi^{+}(y_{t},p_{t})=-\kappa^{\prime}p_{t}y_{t}&\text{
if }y_{t}\geq 0\\\ \Xi^{-}(y_{t},p_{t})=\kappa p_{t}y_{t}&\text{ if
}y_{t}<0\end{cases}$ (1)

, with $\kappa^{\prime},\kappa\geq 0$.

####  III-A2 The BESS

Let $\rho_{c}$ and $\rho_{d}$ be the conversion efficiencies of charge and
discharge respectively. Thus, the round-trip efficiency is
$\rho=\rho_{c}\rho_{d}$. The quantity $r_{t}$ is measured from the exterior of
the BESS, so the energy effectively stored in the battery is $\rho_{c}r_{t}$
(for $r_{t}>0$) and the energy effectively extracted from the battery is
$-r_{t}/\rho_{d}$ (for $r_{t}<0$).

The capacity (size) of the battery is $b$ (in MWh). In order to determine an
optimal size, we need to model the cost of installing and operating the
battery. We express this cost as an amortized cost $c$ in $/(MWh$\cdot$h). By
doing this, we assume that the battery is replaced at the end of its lifetime
with the same capacity and at the same cost. It is common to break down the
cost of storage into a cost for the energy storage capacity and a cost for the
inversion capacity. In this model, we do not consider inversion capacity
restrictions; this could be a feature to include in extensions of our model.

An important factor in BESS is energy dissipation, i.e. the proportion of
energy stored in the BESS that is spontaneously lost without any charging or
discharging performed. This is usually expressed as a fraction $\eta$ of
stored energy per unit time. This is inconvenient for our fluid queue model,
as described below. So instead, we model dissipation as a constant leakage.

####  III-A3 Objective

The profit during time period $t$ is thus
$\Pi_{t}=p_{t}q_{t}-\Xi_{t}\big{(}q_{t}-(w_{t}-r_{t})\big{)}$. In the previous
paragraphs, we have referred to $t$ as a period for ease of exposition, given
its similarity with standard electricity markets. However, we propose here a
continuous-time model, so $t$ actually, and in all instances in the remainder
of this document, refers to an instant, and, consequently, the quantities
$q_{t},w_{t},r_{t}$ are powers (in MW). Note that the price $p_{t}$ is indeed
in $/MWh, so that our profit $\Pi$ is an instantaneous profit rate, in $/h.

We are interested in the long-run average profit, as defined in (2). The
expectation is taken with respect to the stochastic process of interest here,
($w_{t},p_{t}$), as described next.

$\Pi=\lim_{T\to\infty}\mathbb{E}\left[\frac{1}{T}\int_{0}^{T}\Pi_{t}dt\right]-cb$
(2)

####  III-A4 Sources of uncertainty

We consider two sources of uncertainty in our model: wind power output and
energy prices. This is done by considering $(w_{t},p_{t})$ as a joint
continuous-time stochastic process. In particular, we model it as a
continuous-time Markov chain (CTMC) with a discrete state space
$\mathcal{S}=\mathcal{W}\times\mathcal{P}$, with $\mathcal{W}$ and
$\mathcal{P}$ being the discrete state spaces of $w_{t}$ and $p_{t}$
respectively. This stylized model allows for a convenient formulation of
limiting distributions, as described in the next subsection.

There are several examples of modeling wind speed and wind power output as a
Markov chain in the literature, with applications in simulation of wind data
series [17, 18, 19] but also in long-run analysis [11, 10]. It has been found
that simply performing a max likelihood parameter estimation can result in a
very good approximation of the limiting distribution, but that a key metric to
obtain a more accurate model is the autocorrelation [20]. In [20], in the
context of co-located storage sizing for robust operation of microgrids, it is
found that Markovian models with an autocorrelation that poorly reflects that
of the original data series can lead to underestimation of necessary storage
by as much as 50%. Autocorrelation performance can be improved by increasing
the order of the Markov chain, but this increases the size of the state space
exponentially, so that only chains of second or third order are of practical
relevance. In [21], a rolling-average method is proposed to obtain higher
autocorrelation performance in the lower range (0-6h) without increasing the
size of the state space. More recently, authors have proposed non-homogeneous
Markov chains [22, 19] to more accurately reproduce the autocorrelation of
real data series, with much better results replicating the daily behavior of
wind (i.e. the autocorrelation around integer multiples of 24h). These results
can also be leveraged to capture seasonal changes of longer duration in the
wind distribution. The expected range of optimal storage sizes should also
inform the decision of how much autocorrelation needs to be captured by the
model. If these sizes lie in the 0-6 hour range, capturing daily
autocorrelation is less crucial than if storage is expected to be in the 24h+
range.

The discussion above refers only to the Markovian modelling of wind speed or
wind power output and provides enough justification for our modelling of wind
power output as a CTMC in the constant prices case that is addressed in this
paper. We could find no references in the literature of modelling wind power
production and electricity prices as a joint CTMC as our general framework
proposes. This might be regarded as an acceptable approximation for the high-
level steady-state analysis considered here, but some additional
justification, including some empirical support, should be provided when
developing that stage of the model. We leave a more thorough discussion of
this for a future work in which the variable-price version of this model is
developed. For the case considered in this paper, we discuss some empirical
performance of the model in section IV-A3.

####  III-A5 Optimization model

Putting the previous pieces together, we are interested in solving the
infinite-horizon average profit optimization problem in (III-A5).

—s— q_t,r_t,bΠ b≥0 q_t∈[0,W] ∀t q_t,r_t,b∈R ∀t

$W$ is the plant capacity (in MW). For simplicity, we are neglecting any
ramping charge/discharge limits in this model. Furthermore, we intend here to
model an electricity market, so we will restrict our analysis to stationary
policies where the bid is a function of the price, i.e. $q_{t}=q(p_{t})$. We
can also write this as $q_{t}=q^{s}$ when $p_{t}=p^{s}$, for
$s\in\mathcal{S}$.

Finally, we will also be interested in stationary charge-discharge policies
$r_{t}=r(w_{t},q_{t},p_{t})$, which we write
$r_{t}=r(w_{t},q_{t},p_{t})=r^{s}$ when
$(w_{t},q_{t},p_{t})=(w^{s},q^{s},p^{s})$. We must note, however, that this
definition must be overriden if the storage is empty or full as summarized
below.

$\displaystyle r^{s}>0\text{ \& storage full }\Rightarrow r_{t}=0$
$\displaystyle r^{s}<0\text{ \& storage empty }\Rightarrow r_{t}=0$

In this paper, we focus on cases where strategic storage of wind power
generation is not attractive because prices and imbalance penalties are
constant in time. For this case, it is known that the optimal (cost-
minimizing) policy is a balancing policy, i.e. $r_{t}=q_{t}-w_{t}$ [4, 2].

###  III-B Steady-state approach

By the ergodicity of CTMC, and since we are restricting our analysis to
stationary policies, we can express the long-run average profit of (2) in
terms of limiting distributions as in (3). In words, the long-run average
profit is the sum over all states of the income minus the imbalance penalty,
for which there are two cases: if storage is available, the policy can be
followed; if storage is not available (empty or full), the policy must be
overriden.

$\displaystyle\Pi=$
$\displaystyle\sum_{s\in\mathcal{S}}\Big{(}p^{s}q^{s}\pi^{s}-\psi^{s}\Xi^{s}(q^{s}-w^{s})$
$\displaystyle-(\pi^{s}-\psi^{s})\Xi^{s}(q^{s}+r^{s}-w^{s})\Big{)}-cb$
$\displaystyle=$
$\displaystyle\sum_{s\in\mathcal{S}}\Big{(}p^{s}q^{s}\pi^{s}-\psi^{s}\Xi^{s}(q^{s}-w^{s})\Big{)}-cb$
$\displaystyle=$
$\displaystyle\sum_{s\in\mathcal{S}}\Big{(}p^{s}q^{s}\pi^{s}-\psi^{s}\kappa
p^{s}(q^{s}-w^{s})^{+}$
$\displaystyle+\psi^{s}\kappa^{\prime}p^{s}(w^{s}-q^{s})^{+}\Big{)}-cb$ (3)

, where $\boldsymbol{\pi}=(\pi^{s})_{s\in\mathcal{S}}$ is the limiting
distribution of the CTMC $(w_{t},p_{t})$ and
$\boldsymbol{\psi}=(\psi^{s})_{s\in\mathcal{S}}$ is the long-run probability
of storage unavailability (empty or full). The optimal charge/discharge policy
is balancing the output, so $q^{s}+r^{s}-w^{s}=0$ whenever storage is
available, and the imbalance penalty in that case is zero, which justifies the
second equality. The third equality is obtained by plugging in the definition
of the imbalance in (1). $\boldsymbol{\pi}$ can be easily determined from the
generator of the CTMC $(w_{t},p_{t})$. The long-run probability of unavailable
storage $\boldsymbol{\psi}$ can be determined from some results of fluid queue
theory, as shown next.

####  III-B1 Limiting distribution

Our model corresponds to the model of a Markov-modulated fluid queue with
finite buffer. A characterization of the long-run distribution of this process
can be obtained through a spectral analysis, which we overview next. This is
based on the presentation in [12], with more details of the proofs available
in [7]. For ease of exposition, in the following overview we omit the
efficiency factor $\rho$.

Define $\mathbf{r}=[r^{s}]_{s\in\mathcal{S}}$, a vector with the discrete
values taken by $r_{t}$, $\mathbf{D}=\mathrm{diag}(\mathbf{r})$, a diagonal
matrix with $\mathbf{r}$ in its diagonal. In the field of fluid queues,
$\mathbf{r}$ is called the drift vector. We assume for now that $r^{s}\neq
0\forall s\in\mathcal{S}$. The special case with $r^{s}=0$ is considered at
the end. Let $\mathbf{Q}$ be the infinitessimal generator matrix of the CTMC
and $\mathbf{F}$ the limiting distribution of the level state of the battery,
i.e.: $F(x,s)=\lim_{t\to\infty}\mathbb{P}\Big{(}X_{t}\leq
x,(w_{t},p_{t})=s\Big{)}$, $\mathbf{F}(x)=[F(x,s)]_{s\in\mathcal{S}}$. Then,
it can be shown [12] that $\mathbf{F}$ satisfies the differential equation

$\frac{d\mathbf{F}}{dx}\mathbf{D}=\mathbf{F}\mathbf{Q}$ (4)

with boundary conditions

$\displaystyle F(0,s)=0$ $\displaystyle\text{ if }r^{s}>0$ $\displaystyle
F(b,s)=\pi^{s}$ $\displaystyle\text{ if }r^{s}<0$

A spectral solution to these equations can be obtained introducing generalized
eigenvalues $\lambda$ and eigenvectors $\mathbf{u}$, so that
$\lambda\mathbf{u}\mathbf{D}=\mathbf{u}\mathbf{Q}$. The general solution to
(4) takes then the form

$\mathbf{F}(x)=\sum_{i=1}^{|\mathcal{S}|}a_{i}\exp(\lambda_{i}x)\mathbf{u}^{i}$

, where the values of coefficients $a_{i}$ can be found by solving a linear
system from the boundary conditions. The long-run probability of unavailable
storage can then be found as

$\boldsymbol{\psi}=\mathbf{F}(0)+\boldsymbol{\pi}-\mathbf{F}(b)$ (5)

In this manner, the long-run probability of unavailable storage can be found
for given $\mathbf{r}$ and $b$. Note, however, that the method requires
solving a generalized eigenvalue problem, which has three implications of
importance for our work. First, we cannot obtain a closed-form expression of
(5) in terms of $b$ and $\mathbf{r}$, so that numerical calculating approaches
are necessary. Secondly, this function is not convex in general, which makes
our optimization problem non-convex as well. Finally, the linear system posed
by the boundary conditions can be very ill-conditioned because of the presence
of both very large and very small eigenvalues, which is a major challenge for
the method. We address this in more detail in the description of the
algorithm.

####  III-B2 Cases with zero drift

If $r^{s}=0$ for some state $s\in\mathcal{S}$, then $F(x,s)$ can be expressed
as a linear combination of $F(x,j)$ for states $j:r^{j}\neq 0$. Thus, for
these cases, the method described above is performed on a reduced system that
includes only states $\\{j:r^{j}\neq 0\\}$. Then those values are used to find
the distribution for the null states. Details are omitted here, but can be
found in the Appendix of [14].

####  III-B3 Dissipation

We next briefly address the issue of modelling dissipation under our model. As
mentioned before, it is standard to consider that energy dissipates from the
battery at a rate that is a fix multiple of the current storage level. This
would mean that an additional term should be included in the drift vector:
$r_{t}=r^{s}-\eta x_{t}$, where $x_{t}$ is the current storage level in the
battery (in MWh). The addition of this term makes the drift dependent on the
storage level, so that the spectral analysis performed earlier would be no
longer valid. To avoid this, we model dissipation as a state-independent
constant leakage. So that $r_{t}=r^{s}-\eta b$.

This affects the optimality of the balance policy. Indeed, for sufficiently
large $\eta$, it could be more profitable to sell surplus energy immediately
than to store it and have a large proportion of that lost to dissipation. To
simplify our analysis, we make the following restriction: we make a
sensitivity analysis on $\eta$ with $\kappa^{\prime}=0$ and a sensitivity
analysis on $\kappa^{\prime}$ with $\eta=0$, so that the balance policy is
optimal for all these situations.

###  III-C Units

In the previous paragraphs, we defined all quantities in their appropriate
physical units. However, it is more convenient and illustrative for the
purposes of this work to express all quantities in per unit of power plant
capacity and storage capacity, so we introduce:

$\displaystyle\mathbf{w}$ $\displaystyle=\tilde{\mathbf{w}}W$
$\displaystyle\quad q$ $\displaystyle=\tilde{q}W$ $\displaystyle\quad b$
$\displaystyle=\tilde{b}W$ $\displaystyle c$ $\displaystyle=\tilde{c}p^{max}$
$\displaystyle\quad p$ $\displaystyle=\tilde{p}p^{max}$
$\displaystyle\quad\mathbf{r}$ $\displaystyle=\tilde{\mathbf{r}}b$

, so that $\tilde{\mathbf{w}},\tilde{q},\tilde{p}\in[0,1]$, $\tilde{b}$ is in
hours of storage of full plant capacity, and $\tilde{\mathbf{r}}$ is in units
of [p.u. of $b$]/h, which is interpreted as the number of times that the total
storage would be charged starting from an empty state in one hour at full
plant capacity. Equation (3) becomes:

$\displaystyle\frac{\Pi}{p^{max}W}=$
$\displaystyle\sum_{s\in\mathcal{S}}\Big{(}\tilde{p}^{s}\tilde{q}^{s}\pi^{s}-\psi^{s}\kappa\tilde{p}^{s}(\tilde{q}^{s}-\tilde{w}^{s})^{+}$
$\displaystyle+\psi^{s}\kappa^{\prime}\tilde{p}^{s}(\tilde{w}^{s}-\tilde{q}^{s})^{+}\Big{)}-\tilde{c}\tilde{b}$
(6)

It is not hard to check that this change of units does not affect the spectral
decomposition and hence the values of $\boldsymbol{\psi}$. In this paper we
are interested in the case where the price is constant, i.e.
$|\mathcal{P}|=1$, $p^{s}=p^{max}\;\forall s\in\mathcal{S}$ and hence
$\tilde{p}^{s}=1$ for all $s$. Finally, for readability and ease of notation,
the tildes will be omitted in the remainder of this paper, but we will always
refer to values in per unit.

With these unit changes and simplifications, the optimization model (III-A5)
is written as

—s— q,b ∑_s∈S( qπ^s - ψ^sκ(q-w^s)^+ +ψ^sκ’ (w^s-q)^+ ) \- cb b≥0 q∈[0,1]

###  III-D Solution algorithm

There are two challenges in finding the optimal solution to our problem.
First, the spectral method used for computing the objective function is
numerically unstable for certain values of $\mathbf{r}$, for which the
solution involves a very ill-conditioned matrix. Using Matlab’s symbolic
toolbox allows working with a degree of precision that reduces this problem.
However, this significantly increases the computation time of the objective
function. This difficulty can be circumvented with heuristic approaches as
described next. First, a computation in normal precision is tried. If boundary
conditions fail to be respected within a predefined tolerance, the computation
is repeated using the symbolic toolbox. We impose some precision limitations
to reduce execution time.

It was observed that the drift vectors that elicit this behavior are those
that include entries with very small absolute value. In the algorithm proposed
in [13], to guarantee obtaining results within a desired precision, the author
imposes the restriction that entries in the drift vector should have a known
infimum, which is in line with the behavior observed here. We use the
following heuristic: if numeric instability persists after the previous
procedure, a lower threshold is applied to the $\mathbf{r}$ vector.
Specifically, all entries with an absolute value relative to the maximum
absolute entry less than a certain threshold $r_{inf}$ are forced either to 0
or to the threshold, whichever is closer. If the problem persists, the
threshold is increased until a stable solution is found.

This restriction does not make our solution too limited or unrealistic, since
real charge/discharge equipment often does have a minimum threshold for
operation, because of inverter and transformer’s limitations, so that working
with $r_{inf}=0.1$ is reasonable. A more transparent way of introducing this
heuristic rule would be to include it in the constraints of the optimization
problem. However, this would make our problem a non-linear mixed-integer
program, further complicating it without much perceived benefit. The more
informal heuristic described above is thus preferred. In a strict sense,
enforcing this limitation implies not using the optimal balancing policy.
However, in the cases where this happens, we disregard the effect that using
the suboptimal policy has on the optimality of the result.

An interesting alternative to address this numerical difficulty is using the
algorithm proposed by Akar in [15]. This is discussed in more detail in the
following subsection.

The second difficulty is that we have a non-linear, non-convex function, so
that using a gradient descent algorithm does not provide a guarantee of global
optimality. For the scenario with constant prices, our search space only has a
dimension of 2; it is thus reasonable to start exploring the search space by
means of a grid, to then use the best candidate as starting point for a
finite-difference descent algorithm, which is the algorithm that was
implemented in our tests.

####  III-D1 Eigenvalue-free algorithm

In [15], Akar and Sohraby propose an algorithm that allows obtaining the
limiting distribution of a Markov-modulated fluid queue without solving an
eigenvalue problem, to get around the problem of ill-conditioned matrices
present in the spectral solution. They exploit a decomposition of the limiting
distribution in what they call a constant, stable and antistable subsystem.
This decomposition allows posing boundary conditions as a linear system
without exponentially growing terms, which is thus better suited for numerical
solution. Although they propose it in the context of communication networks,
their algorithm does not assume any special structure for the driving Markov
process or the drift vector, and can thus be readily used for our framework as
well.

We implemented this algorithm but obtained mixed results, which are presented
in Section IV. In consequence, although this method seems promising for a
continuation of this project, in the results shown in this paper we use the
heuristic approach described earlier.

##  IV Results and discussion

###  IV-A Model data estimation

####  IV-A1 Wind model

To obtain the values of the generator matrix that defines the CTMC, the
methodology outlined in [21] is followed. The methodology consists in first
passing the wind power output through an averaging window of one hour, then
discretizing the output with $N=15$ levels, and finally performing a max-
likelihood estimation on the resulting data sequence, i.e. finding the
transition probabilities from counting transitions in the sequence. As in
[10], the wind data available in [23] was used.

####  IV-A2 Reference parameters

We perform sensitivity analyses on the other model parameters. Reference
values are given in Table I. We comment on the choice of the reference cost of
storage in the following paragraphs.

TABLE I: Reference values for sensitivity analyses Parameter | Value | Parameter | Value
---|---|---|---
$\kappa$ | 1.35 | $\rho_{0}$ | 0.95
$\kappa^{\prime}$ | 0 | $\rho_{1}$ | 0.95
$c_{s}$ | 0.005 | $\eta$ | 0

####  IV-A3 Empirical validation

To provide an empirical reference for comparison with our model, we use the
same data series to solve the ex-post sizing and contracting optimization
problem IV-A3.

—s— q,b g(q,b) b≥0 q∈[0,1]

, where $g(q,b)$ is a function that computes the ex-post average profit over
the time available in the data series for contract quantity $q$ and storage
size $b$. The value of function $g$ is computed by the algorithm below, where
we use the fact that the the balancing policy is known to be optimal.

0: $q,b,\mathbf{w},\kappa,\kappa^{\prime},\rho,\Delta$

$\mathbf{w}${Wind power output data series}

$\Delta${Wind power output data series sampling period}

$T\leftarrow\mathtt{length}(\mathbf{w})$

$x_{0}\leftarrow 0$ {Initial state of charge}

$x,x_{prev}$ {Current, previous state of charge}

$\Pi$ {Average profit}

$\Pi\leftarrow q\times T\times\Delta$

for $t=1$ to $T$ do

if $x_{prev}+\Delta t(w_{t}-q)<0$ then

$\mathtt{empty}\leftarrow$ true{empty storage}

else {Not empty}

$\mathtt{empty}\leftarrow$ false

end if

if $x_{prev}+\rho\Delta(w_{t}-q)>b$ then

$\mathtt{full}\leftarrow$ true{surplus}

else {No surplus}

$\mathtt{full}\leftarrow$ false

end if

if $\mathtt{full}$ then

$\Pi\leftarrow\Pi+\kappa^{\prime}\max\\{\Delta(w_{t}-q)-(b-x_{prev})/\rho,0\\}$

end if

if $\mathtt{empty}$ then

$\Pi\leftarrow\Pi-\kappa\max\\{0,\Delta(q-w_{t})-(x_{prev})\\}$

end if

if $w_{t}-q<0$ then

$x\leftarrow\max\\{\min\\{x_{prev}+\Delta(w_{t}-q),b\\},0\\}$

else

$x\leftarrow\max\\{\min\\{x_{prev}+\rho\Delta(w_{t}-q),b\\},0\\}$

end if

$x_{prev}\leftarrow x$

end for

$\Pi\leftarrow\frac{\Pi}{T\Delta}$

return $g(q,b)=\Pi$

The results of the model (III-C) and (IV-A3) are shown in Figure 1a and 1b.
The results themselves are discussed in the next paragraphs; in terms of model
validation, it is worth mentioning that there is a good level of agreement
between the curves given by the model and the empirical ex-post best. An
interesting continuation of this work would be to perform a similar analysis
using a non-homogeneous Markov chain to model the wind power output as in [22,
19].

###  IV-B Value of storage

As a first step before doing sensitivity analyses, we are interested in
observing the value of storage for the WPP in the setting described. To do
this, we fix $b$ at different values and solve (III-C) and (IV-A3) for $q$
only to find the optimal commitment and corresponding profit. The profit is
compared to the profit that would be obtained from a feed-in-tariff contract
with the same price. In the setting that we consider here, there is no access
to markets other than energy, so this profit is an upper-bound reference for
comparison.

We can make four important observations from the results in Figure 1. First,
as expected, even with unlimited storage, recovering feed-in-tariff profits is
not possible, which is explained by efficiency losses in storage/discharge.
Secondly, it is also worth mentioning that in our model, we do not find an
optimal bidding function having the form of an optimal fractile of wind power
production. Observe that since the support $\mathcal{W}$ is discrete, such a
function would be step-shaped, while our optimal curve is not.

Third, in line with the results of other authors under different structural
assumptions [4, 2], the value of storage in the reference case studied has
diminishing marginal profit. This allows identifying two quantities of
interest. First, the marginal profit at $b=0$ gives the critical amortized
cost of storage: if storage costs more than this quantity, it is not worth
having co-located storage. The second quantity is the optimal size; for an
amortized cost of storage $c$, the optimal size of storage is the quantity for
which the marginal profit is equal to $c$. These results justify our choice of
a descent algorithm to find the optimal storage size.

Finally, as expected, expected imbalances are reduced as storage sizes grows
without limit. It is worth observing, however, that positive imbalances
increase for the lowest range of storage sizes. This can be explained by the
fact that in our model for the given reference values, the optimal contracting
strategy in absence of storage is to commit the entire capacity of the plant,
which eliminates positive imbalances. As storage drives the optimal commitment
down, positive imbalances appear. With enough storage, these imbalances are
then brought down too.

For the reference case shown here, the critical cost is found at
$c_{o}=0.0193$, which for a reference price of energy of 60$/MWh, corresponds
to an amortized cost of 10.14$/kWh-yr. Note that this is a whole order of
magnitude below the reference cost of 100$/kWh-yr for 2020 [24], meaning that
under this setting and the reference values used, installing co-located
storage is not profitable without access to additional revenue streams (such
as a capacity market or arbitrage). Even taking into account that the cost of
storage is declining rapidly, a tenfold reduction is beyond what we should
expect to see in the near future. So additional revenue streams or a different
setting need to be in place for co-located storage to be attractive.

In fact, in our setting, the WPP is not sufficiently exposed to uncertainty to
justify paying for expensive storage. Indeed, the WPP is already covered
against prices higher than $\kappa p$ by the pay-as-demanded forward contract,
so that further risk hedging via storage is only modestly attractive. In
Figure 2, we see that as the coverage cap $\kappa$ is increased and the WPP is
exposed to higher penalties, the critical cost $c_{o}$ at which it is willing
to invest in storage also rises. However, the pay-as-demanded forward contract
would need to be available at a premium of several times the WPP’s normal
selling price to get anywhere near the critical cost.

A setting in which the WPP is exposed to much higher penalties is thus
necessary for the co-located investment to be attractive, i.e. a setting with
variable prices, which is the natural continuation of the work presented in
this paper.

(a)

(b)

(c)

Figure 1: (1a): Long-run average profit in excess of no-storage profit for
different storage sizes. The feed-in-tariff benchmark is equal to the average
output of the power plant (price is normalized to 1). (1b): Optimal energy
pledge for different storage sizes. (1a) and (1b) correspond to the reference
values of all parameters. (1c): Expected positive and negative imbalance for
different storage sizes. Figure 2: Critical storage cost for different values
of shortfall penalty factor. Critical storage cost is the slope of the profit
curve vs. storage size at $b=0$.

Finally, we obtain the identical curves in Figure 1 using the spectral method
described in Section III-D and the alternative algorithm of [15], with the
latter having a much quicker execution. However, this agreement of results did
not hold for some other values of the parameters, such that using that method
for sensitivity analyses was not possible. Indeed, it seems that for some
values of the parameters that were varied, the alternative algorithm can
return invalid limiting distributions (outside $[0,1]$). This behavior has not
been investigated any further but could be of interest for future extensions
of the model.

###  IV-C Sensitivity analyses

As we mentioned, the main scope of the model proposed here is performing
sensitivity analyses at a high level, such as during a feasibility study of a
project, or for policy evaluation. In this section, we show the results of
some of these analyses. Despite the discussion of the previous section,
sensitivity analyses to changes in some key parameters are of interest, since
they provide some understanding regarding the interplay between some key
quantities which a project developer could be interested in evaluating in
early stages of an investment.

We are interested in observing the behavior of the model as some key
parameters are modified. The current cost of storage, $c\approx 0.2$ is above
the critical cost for the reference case, so the optimal size would be 0, and
the behavior of the results to changes in other parameters would be hidden. To
avoid this, we take a reference storage cost of $c=0.005$. Although
unrealistically low, it allows observing the behavior of the optimal size as
some parameters of interest change. This could also be seen as the case of co-
located storage with access to multiple revenue streams (e.g. capacity
market), and only a portion of the full cost of the BESS needs to be recovered
via risk mitigation in forward markets.

####  IV-C1 Forward contract prices

The behavior is exactly as intuition would dictate: larger storage becomes
more attractive as the penalty for energy shortfall becomes higher, i.e. as
the exposure to high penalties increases (Figure 3). The change with respect
to the price of the pay-as-generated forward contract (Figure 4) shows a more
interesting dynamic. As $\kappa^{\prime}$ becomes positive and the value of
surplus energy increases, it becomes less undesirable to have excess energy,
so that the optimal commitment $q$ decreases. This in turn leaves room for
increasing the optimal storage size, as more energy is available for
accumulation. After a certain point, however, the price of the pay-as-
generated contract is so close to the price of the long-term fixed contract
that the latter becomes less and less attractive, so that the commitment
decreases significantly. With a very low commitment in the long-term contract,
the motivation for investing in storage also disappears progressively, leading
to a smaller optimal storage size.

####  IV-C2 Efficiency and dissipation

The behavior when the round-trip efficiency and dissipation vary is in line
with our expectation. As the system becomes more inefficient, storage becomes
less attractive. Figure 6 shows this behavior as dissipation varies. For the
round-trip efficiency parameter $\rho$, the relation between optimal size and
efficiency appears to be close to linear (Figure 5).

(a)

(b)

(c)

Figure 3:  Sensitivity to changes in shortfall penalty $\kappa$. (3a): Profit
gain with respect to no-storage profit (3b): Optimal storage size (3c):
Optimal quantity that should be sold in long-term fixed contract.

(a)

(b)

(c)

Figure 4:  Sensitivity to changes in discount factor $\kappa^{\prime}$ of
surplus energy sales. (4a): Profit gain with respect to no-storage profit
(4b): Optimal storage size (4c): Optimal quantity that should be sold in long-
term fixed contract.

(a)

(b)

(c)

Figure 5:  Sensitivity to changes in round-trip efficiency $\rho$. (5a):
Profit gain with respect to no-storage profit (5b): Optimal storage size (5c):
Optimal quantity that should be sold in long-term fixed contract.

(a)

(b)

(c)

Figure 6:  Sensitivity to changes in storage dissipation factor $\eta$. (6a):
Profit gain with respect to no-storage profit (6b): Optimal storage size (6c):
Optimal quantity that should be sold in long-term fixed contract.

##  V Conclusion and next steps

In sum, we are proposing here a model to perform a high-level steady-state
analysis of the value of co-located storage for a wind power producer that
participates in the electricity market through long-term forward contracts. In
particular, we assess the optimal size of storage and optimal quantity to sell
in forward contracts under different values of key parameters regarding
contract prices and storage efficiency. According to the model presented here,
a setting such as the one considered in this paper does not expose the
producer to enough uncertainty in income for storage to be attractive at
current prices. This is likely to be different under a setting with variable
prices, which is the natural extension of the model presented here, and what
we intend to develop in future research.

## References

  * [1] J. H. Kim and W. B. Powell, “Optimal Energy Commitments with Storage and Intermittent Supply,” _Operations Research_ , vol. 59, no. 6, pp. 1347–1360, Dec. 2011, publisher: INFORMS. [Online]. Available: https://pubsonline.informs.org/doi/abs/10.1287/opre.1110.0971
  * [2] E. Y. Bitar, R. Rajagopal, P. P. Khargonekar, K. Poolla, and P. Varaiya, “Bringing Wind Energy to Market,” _IEEE Transactions on Power Systems_ , vol. 27, no. 3, pp. 1225–1235, Aug. 2012, conference Name: IEEE Transactions on Power Systems.
  * [3] C. J. Dent, J. W. Bialek, and B. F. Hobbs, “Opportunity Cost Bidding by Wind Generators in Forward Markets: Analytical Results,” _IEEE Transactions on Power Systems_ , vol. 26, no. 3, pp. 1600–1608, Aug. 2011, conference Name: IEEE Transactions on Power Systems.
  * [4] P. Harsha and M. Dahleh, “Optimal Management and Sizing of Energy Storage Under Dynamic Pricing for the Efficient Integration of Renewable Energy,” _IEEE Transactions on Power Systems_ , vol. 30, no. 3, pp. 1164–1181, May 2015, conference Name: IEEE Transactions on Power Systems.
  * [5] Y. Ru, J. Kleissl, and S. Martinez, “Storage Size Determination for Grid-Connected Photovoltaic Systems,” _IEEE Transactions on Sustainable Energy_ , vol. 4, no. 1, pp. 68–81, Jan. 2013, conference Name: IEEE Transactions on Sustainable Energy.
  * [6] A. Loukatou, P. Johnson, S. Howell, and P. Duck, “Optimal valuation of wind energy projects co-located with battery storage,” _Applied Energy_ , vol. 283, p. 116247, Feb. 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0306261920316391
  * [7] D. Mitra, “Stochastic Theory of a Fluid Model of Producers and Consumers Coupled by a Buffer,” _Advances in Applied Probability_ , vol. 20, no. 3, pp. 646–676, 1988, publisher: Applied Probability Trust. [Online]. Available: http://www.jstor.org/stable/1427040
  * [8] D. Gaver and R. Miller, “Limiting distributions for some storage problems,” in _Studies in applied probability and management science_ , K. Arrow, S. Karlin, and H. Scarf, Eds. Stanford University Press, 1962, pp. 110–126.
  * [9] S. S. Oren and S. H. Wan, “Optimal Strategic Petroleum Reserve Policies: A Steady State Analysis,” _Management Science_ , vol. 32, no. 1, pp. 14–29, Jan. 1986, publisher: INFORMS. [Online]. Available: https://pubsonline.informs.org/doi/abs/10.1287/mnsc.32.1.14
  * [10] X. Chen, K.-C. Leung, and A. Y. S. Lam, “Power Output Smoothing for Renewable Energy System: Planning, Algorithms, and Analysis,” _IEEE Systems Journal_ , vol. 14, no. 1, pp. 1034–1045, Mar. 2020, conference Name: IEEE Systems Journal.
  * [11] V. Deulkar, J. Nair, and A. A. Kulkarni, “Sizing Storage for Reliable Renewable Integration,” in _2019 IEEE Milan PowerTech_ , Jun. 2019, pp. 1–6.
  * [12] V. Kulkarni, “Fluid Models for Single Buffer Systems,” in _Frontiers in Queueing: Models and Applications in Science and Engineering_ , J. H. Dshalalow, Ed. CRC Press, Jan. 1997, pp. 321–338, google-Books-ID: t4S4i8HkJXkC.
  * [13] B. Sericola, “A Finite Buffer Fluid Queue Driven by a Markovian Queue,” _Queueing Systems_ , vol. 38, no. 2, pp. 213–220, Jun. 2001. [Online]. Available: https://doi.org/10.1023/A:1010962516045
  * [14] B. Sericola and B. Tuffin, “A fluid queue driven by a Markovian queue,” _Queueing Systems_ , vol. 31, no. 3, pp. 253–264, Jul. 1999. [Online]. Available: https://doi.org/10.1023/A:1019114415595
  * [15] N. Akar and K. Sohraby, “Infinite- and finite-buffer Markov fluid queues: a unified analysis,” _Journal of Applied Probability_ , vol. 41, no. 2, pp. 557–569, Jun. 2004, publisher: Cambridge University Press.
  * [16] M. Bolinger, J. Seel, C. Warner, and D. Robson, “Utility-Scale Solar, 2021 Edition,” Lawrence Berkeley National Laboratory, Technical Report, Oct. 2021. [Online]. Available: https://emp.lbl.gov/sites/default/files/utility_scale_solar_2021_edition_slides.pdf
  * [17] V. V. Lopes, T. Scholz, A. Estanqueiro, and A. Q. Novais, “On the use of Markov chain models for the analysis of wind power time-series,” in _2012 11th International Conference on Environment and Electrical Engineering_ , May 2012, pp. 770–775.
  * [18] J. Tang, A. Brouste, and K. L. Tsui, “Some improvements of wind speed Markov chain modeling,” _Renewable Energy_ , vol. 81, pp. 52–56, Sep. 2015\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0960148115001871
  * [19] K. Xie, Q. Liao, H.-M. Tai, and B. Hu, “Non-Homogeneous Markov Wind Speed Time Series Model Considering Daily and Seasonal Variation Characteristics,” _IEEE Transactions on Sustainable Energy_ , vol. 8, no. 3, pp. 1281–1290, Jul. 2017, conference Name: IEEE Transactions on Sustainable Energy.
  * [20] K. Brokish and J. Kirtley, “Pitfalls of modeling wind power using Markov chains,” in _2009 IEEE/PES Power Systems Conference and Exposition_ , Mar. 2009, pp. 1–6.
  * [21] G. Papaefthymiou and B. Klockl, “MCMC for Wind Power Simulation,” _IEEE Transactions on Energy Conversion_ , vol. 23, no. 1, pp. 234–240, Mar. 2008, conference Name: IEEE Transactions on Energy Conversion.
  * [22] T. Scholz, V. V. Lopes, and A. Estanqueiro, “A cyclic time-dependent Markov process to model daily patterns in wind turbine power production,” _Energy_ , vol. 67, pp. 557–568, Apr. 2014. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0360544214000048
  * [23] NREL, “Wind generation profile in 2012,” National Renewable Energy Laboratory, Tech. Rep., Mar. 2016.
  * [24] K. Mongird, V. Viswanathan, J. Alam, C. Vartanian, V. Sprenkle, and R. Baxter, “2020 Grid Energy Storage Technology Cost and Performance Assessment,” Pacific Northwest National Laboratory, Tech. Rep. DOE/PA-0204, Dec. 2020.


# Particle track reconstruction with noisy intermediate-scale quantum
computers

T Schwägerl1,2 C Issever1,2 K Jansen2 T J Khoo1 S Kühn2 C Tüysüz1,2 and H
Weber1 1 Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15,
12489 Berlin, Germany  2 Deutsches Elektronen-Synchrotron DESY, Platanenallee
6, 15738 Zeuthen, Germany<EMAIL_ADDRESS>
###### Abstract

The reconstruction of trajectories of charged particles is a key computational
challenge for current and future collider experiments. Considering the rapid
progress in quantum computing, it is crucial to explore its potential for this
and other problems in high-energy physics. The problem can be formulated as a
quadratic unconstrained binary optimization (QUBO) and solved using the
variational quantum eigensolver (VQE) algorithm. In this work the effects of
dividing the QUBO into smaller sub-QUBOs that fit on the hardware available
currently or in the near term are assessed. Then, the performance of the VQE
on small sub-QUBOs is studied in an ideal simulation, using a noise model
mimicking a quantum device and on IBM quantum computers. This work serves as a
proof of principle that the VQE could be used for particle tracking and
investigates modifications of the VQE to make it more suitable for
combinatorial optimization.

##  1 Introduction

At the Large Hadron Collider (LHC) thousands of particles are produced in
collisions of two opposing proton beams. The collision points are enclosed by
detectors whose inner-most components are multiple silicon layers with
cylindrical symmetry around the beam-pipe. One key computational challenge is
particle tracking, i.e. the reconstruction of trajectories of charged
particles from small energy deposits in the detector, referred to as hits. In
Fig. 1 hits and reconstructed trajectories of particles in the transverse
$x-y$-plane of the detector are shown, perpendicular to the beam-direction.

Figure 1: Hits and reconstructed trajectories of particles in the transverse
plane of the detector. The detector layers are indicated by the grey
concentric circles. Real tracks (indicated in green) are successfully
reconstructed. Tracking algorithms can miss or make up tracks, which are
indicated in blue and red. Generated using the Python framework hepqpr-qallse
[1].

Future colliders, such as the high-luminosity LHC (HL-LHC) [2], are in general
designed to produce many particles, rendering the reconstruction more
difficult and asking for improvements, especially in terms of CPU resources.
Quantum computing is an aspiring technology enabling exponential speedups for
certain problems [3]. Currently available quantum devices can have up to
several hundred qubits, but still suffer from a considerable level of noise,
and are commonly referred to as noisy intermediate-scale quantum (NISQ)
hardware [4]. One algorithm for NISQ computers is the hybrid quantum-classical
variational quantum eigensolver (VQE) that was originally proposed to estimate
the ground state energy of molecules [5]. In our study, the potential of VQE
to solve the problem of particle track reconstruction is investigated. There
are several works on quantum algorithms for particle tracking that investigate
different approaches [6, 7, 8, 9, 10]. For a review of quantum computing for
data analysis in high-energy physics, we refer to Ref. [11]. In Ref. [6], the
problem was formulated as a quadratic unconstrained binary optimization (QUBO)
problem, tailored for D-Wave quantum annealers. In this work, we use the same
formulation and investigate the performance of the VQE on gate-based quantum
computers to solve the QUBO.

The paper is organized as follows. In Sec. 2, we explain the QUBO formulation
of particle tracking, our approach for dividing the full QUBO into smaller
sub-QUBOs, and the modifications of the VQE for combinatorial problems we use.
In Sec. 3, we study the impact of sub-QUBO size on the efficiency and the
purity, and show results for the performance of VQE on sub-QUBO level in an
ideal simulation, using a noise model mimicking a quantum device and on real
quantum hardware. Finally, we conclude in Sec. 4.

##  2 Methodology

This work builds upon the approach introduced in Ref. [6] using the data set
of the TrackML challenge [12]. We use its implementation in the Python
framework hepqpr-qallse [1] to construct the QUBO and to compute efficiency
and purity. The main idea is to reconstruct the trajectories from smaller
track segments. The smallest constituents of tracks are single hits, but it is
difficult to reconstruct tracks from hits directly. Two hits $a$ and $b$ in
different layers of the detector form a doublet $(a,b)$. Doublets $(a,b)$ and
$(b,c)$ are combined to triplets $(a,b,c)$. Now, the algorithm aims at
identifying true $(t)$ triplets that are part of trajectories of charged
particles in contrast to false $(f)$ triplets that are random combinations of
three hits. This is done by rating the quality of individual triplets and by
comparing triplets to each other. One key component is to identify
combinations of triplets $(a,b,c)$ and $(b,c,d$) that overlap by two hits and
form a quadruplet $(a,b,c,d)$.

###  2.1 Quadratic unconstrained binary optimization

QUBO problems are combinatorial optimization problems arising in a wide range
of applications. In particular, the problem of particle tracking can be
expressed as a QUBO with a cost function $Q(T)$, where $T$ is a binary vector
of length $N$ whose entries correspond to the triplet candidates $T_{i}$.
Triplets identified by the algorithm to be true (false) triplets are called
positives $(p)$ (negatives $(n)$) and are denoted $T_{i}=1$ ($T_{i}=0$).
Triplet parameters, such as the curvature, are used to compute the
coefficients $a_{i}$ that rate the quality of individual triplets, and the
coefficients $b_{ij}$ that express the compatibility of two triplets. The
details of the computations can be found in Ref. [6]. The coefficients are
constructed such that the minimum of the function $Q(T)$ holds preferably all
true and no false triplets. Using these coefficients, particle tracking is
reduced to the problem of finding the minimum of the cost function

$\displaystyle
Q(T)=\sum_{i}^{N}a_{i}T_{i}+\sum_{i}^{N}\sum_{j<i}^{N}b_{ij}T_{i}T_{j}.$ (1)

Due to the limited size of current quantum devices, we use a geometric
approach to divide the full QUBO into smaller sub-QUBOs that fit on NISQ
computers. The triplets are sorted according to their angle in the
$x-y$-plane, and slices of $k$ consecutive triplets are selected. The slices
overlap by $k/2$ triplets. Subsequently, sub-QUBOs are created for each slice.
We chose this approach to give a perspective on how large the sub-QUBOs have
to be to reasonably solve the tracking problem. $Q(T)$ can be transformed to a
Hamiltonian by mapping $T_{i}\to(1-Z_{i})/2$, where $Z_{i}$ denotes the third
Pauli matrix. Every triplet candidate is represented by one qubit. Now, VQE
can be used to find the ground state of the Hamiltonian and thus the minimum
of $Q(T)$.

###  2.2 Layer variational quantum eigensolver

In Ref. [13] the Layer VQE (L-VQE) approach was introduced for combinatorial
optimization problems. The circuit structure for the L-VQE approach is shown
in Fig. 2.2.

Figure 2: The circuit of the L-VQE algorithm (a). The circuit consists of
single-qubit $R_{Y}$ rotation gates and entangling two-qubit CNOT gates. The
ansatz is grown by adding layers (dashed box on the right hand side) to the
initial rotation layer (dashed box on left hand side). Layer with fewer
$R_{Y}$ gates used for the runs on real hardware (b). Figure adapted from Ref.
[13].

First, a simple rotation layer of $R_{Y}$ gates is initialized with random
parameters uniformly drawn from the interval $[0,2\pi]$. Then, the VQE is run
for this single layer. The optimization process is stopped after a fixed
number of iterations, and not carried out until full convergence. A second
layer including $R_{Y}$ gates and CNOT gates is added to the circuit. Its
initial parameters are set to zero, not affecting the state reached by the
previous optimization. The optimization is then carried out for all
parameters, again for a fixed number of iterations. This procedure is repeated
for every additional layer, and in the last step optimization is carried out
until convergence. In principle, the initial rotation layer is sufficient to
express the solution, since it is a computational basis state because the
Hamiltonian is diagonal. However, the optimization for such an ansatz gets
easily stuck in local minima [14]. The idea behind the additional layers is,
that they might allow the algorithm to escape from these minima.

###  2.3 Conditional value at risk cost function

Originally, the VQE was designed to approximate the ground state of a
Hamiltonian by minimizing the energy expectation value. For diagonal
Hamiltonians it is enough to reach a state with sufficiently large ground
state component that can be revealed by taking enough measurements. Assuming
the final VQE solution has 1% ground state component and we take 1024
measurements, the probability for obtaining the optimal solution in one of the
measurements is greater than 99.996%. Thus, the authors of Ref. [15] proposed
to use a conditional value at risk (CVaR) as a cost function, which uses the
fraction $\alpha$ of measurements with the lowest energy instead of the full
sample mean.

##  3 Results

The geometric sub-QUBO approach allows for rating the full tracking
algorithm’s performance and the performance of VQE on sub-QUBO level
separately. This is an important step because current quantum computers are
still of limited size, thus rendering it unpractical to run the full tracking
algorithm on hardware devices. First, we introduce the relevant metrics to
rate the tracking performance and the performance of VQE on sub-QUBO level.
Subsequently, we present results for the tracking performance for different
sub-QUBO sizes and for different numbers of particles per event, ranging from
current conditions at the LHC (around 2k particles per event) to conditions
expected at the HL-LHC. Finally, we show results of VQE for small sub-QUBOs
and assess how the previously introduced modifications of VQE affect its
performance.

###  3.1 Performance metrics

To rate the tracking algorithm’s performance and to investigate the impact of
dividing the full QUBO into smaller sub-QUBOs, we compute the
$\text{efficiency}=\frac{tp}{tp+fn}$ and the $\text{purity}=\frac{tp}{tp+fp}$
defined by the true positives $(tp)$, the false positives $(fp)$ and the false
negatives $(fn)$ in the final set of doublets. For details on which doublets
that are included in the final set of triplets are considered in the rating of
the algorithm’s performance, we refer to Ref. [6]. To rate the performance of
VQE on sub-QUBO level, we calculate the fraction of instances with at least 1%
ground state component.

###  3.2 Impact of the sub-QUBO size

To assess the impact of the sub-QUBO size on the algorithm’s performance,
various sub-QUBO sizes are studied and solved using neal [16], an
implementation of a simulated annealing sampler. All triplets that are
identified as positives in at least one of the overlapping slices are marked
as positives. Then, efficiency and purity are computed. This estimates the
performance of VQE for particle tracking under the assumption that VQE finds
the ground state of every sub-QUBO. Efficiency and purity for sub-QUBO sizes
ranging from 16 to 512 triplets as well as for the full QUBO are shown in Fig.
3.

(a)(b)

Figure 3: Efficiency (a) and purity (b) as a function of track density for
slices of increasing size in the $x-y$-plane and for the full QUBO. For slices
of size 512 and small track densities, the performance is comparable to
solving the full QUBO.

For sub-QUBO size 16 the efficiency reaches 0.7 for the lowest track densities
and then drops rapidly. Purity is above 0.8 for almost the full range of track
densities. Doubling the size to 32 improves efficiency and purity. The
performance of the VQE for slices of size 16 to 32 is examined in Sec. 3.3.
For size 128, the efficiency is above 0.9 for up to 5k particles per event and
purity is comparable to solving the full QUBO. A sub-QUBO size of 512 is
comparable to solving the full QUBO for up to 5k particles per event in terms
of efficiency and purity.

###  3.3 Performance of the VQE on sub-QUBO level

The performance of the VQE for small sub-QUBOs of size 16 to 32 and for track
densities of 1k particles per event is determined for three scenarios: i) an
ideal simulation assuming a perfect, noise free quantum computer, ii) using a
noise model mimicking the quantum computer ibmq_kolkata, and iii) on the real
device. We use Qiskit [17] for the simulations and hardware runs. Moreover, we
adapt the L-VQE approach with up to 2 additional layers and the CVaR cost
function with $\alpha=0.1$ and $\alpha=1$ (the latter corresponding to running
a conventional VQE). The number of reps indicates the additional layers that
are added to the initial rotation layer. For the initial rotation layer and
after adding an additional layer, 128 iterations of the VQE are performed.
Then, the circuits are further optimized until a total of 1024 iterations is
reached. In every iteration we use 1024 measurements to evaluate the cost
function in all three scenarios. To update the circuit parameters we use
Constrained Optimization BY Linear Approximation (COBYLA) [18]. In Fig. 4 the
mean and 95% confidence interval of fractions of instances with at least 1%
ground state component for 50 slices (5 slices on real hardware) and 10 random
initial points are shown.

Figure 4: Fraction of instances with at least 1% ground state component at the
end of the VQE as a function of sub-QUBO size in an ideal simulation, using
the full noise model of the quantum computer ibmq_kolkata and the real device.
Results for two values of $\alpha$ in the CVaR cost function and for up to two
additional layers in L-VQE are shown. Sizes 28 and 32 show only the ideal
scenario.

Using the CVaR as a cost function with $\alpha=0.1$ yields a significant
improvement for the ideal and noisy simulations. In the ideal setting,
additional layers are beneficial for all system sizes. In the noisy simulation
adding one additional layer yields an improvement for up to 20 qubits. For 32
qubits, around 70% of instances reach at least 1% ground-state component in an
ideal simulation. The data obtained on the real quantum device is compatible
with our classical simulation.

##  4 Discussion and outlook

This work serves as a proof of principle that the VQE can be used for particle
tracking. It was shown that sub-QUBO sizes compatible with NISQ hardware yield
reasonable efficiency and purity. In an ideal simulation without noise, the
VQE finds the solution of small sub-QUBOs for the majority of instances.
Noise, either coming from a noise model mimicking a real quantum device or the
quantum hardware, significantly affects the results. Thus, in future studies
we aim at investigating different techniques for error mitigation. Two
modifications of VQE, the CVaR cost function and L-VQE, were tested. CVaR
shows a clear advantage for combinatorial optimization. While the additional
layers in L-VQE also show a small advantage, the role of entanglement for
classical optimization problems remains an open question asking for further
research.

This work is supported by: i) the Einstein Research Unit -“Perspectives of a
quantum digital transformation: Near-term quantum computational devices and
quantum processors”, ii) the European Union’s Horizon Europe research and
innovation funding programme under the ERA Chair scheme with grant agreement
No. 101087126 and under the European Research Council (ERC) with grant
agreement No. 787331, iii) the Helmholtz Association -“Innopool Project
Variational Quantum Computer Simulations (VQCS)”, and iv) the Ministry of
Science, Research and Culture of the State of Brandenburg within the Centre
for Quantum Technologies and Applications (CQTA).

## References

  * [1] Linder L 2019 hepqpr-qallse URL https://github.com/derlin/hepqpr-qallse
  * [2] Apollinari G, Bruening O, Nakamoto T and Rossi L 2015 CERN Yellow Report 5 1–19
  * [3] Shor P W 1997 SIAM Journal on Computing 26 1484–1509
  * [4] Preskill J 2018 Quantum 2 79
  * [5] Peruzzo A, McClean J, Shadbolt P, Yung M H, Zhou X Q, Love P J, Aspuru-Guzik A and O’Brien J L 2014 Nature Communications 5 4213
  * [6] Bapst F, Bhimji W, Calafiura P, Gray H, Lavrijsen W, Linder L and Smith A 2019 Computing and Software for Big Science 4 1
  * [7] Tüysüz C, Rieger C, Novotny K, Demirköz B, Dobos D, Potamianos K, Vallecorsa S, Vlimant J R and Forster R 2021 Quantum Machine Intelligence 3 29
  * [8] Duckett P, Facini G, Jastrzebski M, Malik S, Rettie S and Scanlon T 2022 Reconstructing charged particle track segments with a quantum-enhanced support vector machine arXiv:2212.07279 [hep-ex]
  * [9] Funcke L, Hartung T, Heinemann B, Jansen K, Kropf A, Kühn S, Meloni F, Spataro D, Tüysüz C and Yap Y C 2023 J. Phys. Conf. Ser. 2438 012127
  * [10] Zlokapa A, Anand A, Vlimant J R, Duarte J M, Job J, Lidar D and Spiropulu M 2021 Quantum Machine Intelligence 3 27
  * [11] Delgado A et al. 2022 Quantum Computing for Data Analysis in High-Energy Physics arXiv:2203.08805 [physics.data-an]
  * [12] Kaggle 2018 Trackml particle tracking challenge URL https://www.kaggle.com/c/trackml-particle-identification
  * [13] Liu X, Angone A, Shaydulin R, Safro I, Alexeev Y and Cincio L 2022 IEEE Transactions on Quantum Engineering 3 1–20
  * [14] Anschuetz E R and Kiani B T 2022 Nature Communications 13 7760
  * [15] Barkoutsos P K, Nannicini G, Robert A, Tavernelli I and Woerner S 2020 Quantum 4 256
  * [16] D-Wave Systems Inc 2019 dwave-neal URL https://github.com/dwavesystems/dwave-neal
  * [17] tA v A et al. 2018 Qiskit: An open-source framework for quantum computing URL https://github.com/Qiskit
  * [18] Powell M J D 1998 Acta Numerica 7 287–336