# Dynamic Price Incentivization for Carbon Emission Reduction using Quantum
Optimization

David Bucher1 , Jonas Nüßlein2, Corey O’Meara3, Ivan Angelov4, Benedikt
Wimmer1,
Kumar Ghosh3 , Giorgio Cortiana3 and Claudia Linnhoff-Popien2
0009-0002-0764-9606 0000-0002-4628-6951 1Aqarios GmbH
Prinzregentenstraße 120, Munich, Germany
Email<EMAIL_ADDRESS>2Mobile and Distributed Systems Chair
Ludwig-Maximilians Universität, Munich
Email<EMAIL_ADDRESS>3E.ON Digital Technology GmbH
Tresckowstraße 5, Hannover, Germany
Email<EMAIL_ADDRESS>4Comsysto Reply GmbH
Riesstraße 22, 80992 Munich, Germany
Email<EMAIL_ADDRESS>
###### Abstract

Demand Side Response (DSR) is a strategy that enables consumers to actively
participate in managing electricity demand. It aims to alleviate strain on the
grid during high demand and promote a more balanced and efficient use of
electricity resources. We implement DSR through discount scheduling, which
involves offering discrete price incentives to consumers to adjust their
electricity consumption patterns. Since we tailor the discounts to individual
customers’ consumption, the Discount Scheduling Problem (DSP) becomes a large
combinatorial optimization task. Consequently, we adopt a hybrid quantum
computing approach, using D-Wave’s Leap Hybrid Cloud. We observe an indication
that Leap performs better compared to Gurobi, a classical general-purpose
optimizer, in our test setup. Furthermore, we propose a specialized
decomposition algorithm for the DSP that significantly reduces the problem
size, while maintaining an exceptional solution quality. We use a mix of
synthetic data, generated based on real-world data, and real data to benchmark
the performance of the different approaches.

##  I Introduction

The increasing demand on energy resources and the growing adaption of
renewable electricity sources have prompted a search for innovative solutions
to optimize energy consumption in order to reduce grid congestion and carbon
emissions. Demand Side Response (DSR) [1] has emerged as a promising strategy
that focuses on actively managing and adjusting energy consumption patterns in
response to grid conditions. Multiple investigations on DSR exist in
literature, describing its impact on smart grid technology [2], load
scheduling [3], energy economics [4], optimal control and pricing schemes [5].

Price adjustment is a simple tool to steer consumer behavior. With the
emergence of smart devices and the electrification of heating and
transportation, the response to price incentives can progressively be
automated. Typically, DSR is achieved by handing out a dynamic price to all
customers simultaneously. Yet, the different usage patterns of the consumers
may favor alternate dynamic price policies. Therefore, we set out to find
individual price patterns on a per-customer basis such that an optimal load
shift can be achieved. We call the distribution of discounts or penalties to
specific customers the discount scheduling problem (DSP). The number of
customers that have to be considered in such a problem, i.e., an urban power
grid, can become prohibitively large to be solved by classical resources.

In recent years, Quantum Computing (QC) has garnered significant attention as
a potential game-changer in various domains, including optimization.
Leveraging the principles of quantum mechanics, quantum optimization
algorithms are hypothesized to solve complex optimization problems more
efficiently than their classical counterparts. Besides gate-based universal
quantum computing, Adiabatic Quantum Computing (AQC) has emerged, which can be
shown in general to be equivalent to gate-based approaches[6]. A subset of
AQC, Quantum Annealing (QA) [7, 8] has been widely adopted for solving
optimization problems [9, 10]. With D-Wave being the industry leader in
providing quantum annealing computing hardware, we use D-Wave’s quantum
annealer in this work to solve the DSP. The limited size of current quantum
computing hardware forces us to utilize hybrid quantum computing approaches,
like Leap, which is a Cloud service offered by D-Wave and is based on internal
problem size reduction [11]. In this work, we additionally develop a
customized hybrid approach that performs a problem-specific decomposition.

The paper is structured as follows: After giving a concise Literature Review
in Sec. II, we describe the problem formulation and mathematical modeling of
the DSP in Sec. III. Since the algorithm should be applicable to large amounts
of customers in realistic scenarios, Sec. IV motivates and develops a problem-
specific decomposition approach for problem size reduction. This decomposition
routine proves to be very effective, as the benchmarking of classical and
quantum-enhanced solvers, based on various criteria, in Sec. V shows.
Furthermore, we observe that Gurobi, as a mixed-integer classical solver,
reaches a limit for larger problem sizes, while D-Wave’s Leap Hybrid Quantum
solver still outputs workable results. Nevertheless, the decomposition routine
aided by a classical solver provides the overall best results.

##  II Literature Review

###  II-A Related Work

Recently, quantum computing applications in the power and energy sector [12,
13, 14, 15] are gaining attention for the development of smart grid
technology. Several important problems are addressed using quantum computing,
for example power flow calculations [16, 17] or energy grid classification
[18]. The traditional planning and scheduling tasks in power systems, such as
the minimization of generation cost or the maximization of revenue from
electricity generation, are generally formulated as combinatorial optimization
problems, which are often NP-hard. Using quantum-inspired optimization
algorithms is expected to outperform their classical counterparts [19, 13]. A
wide range of optimization problems can be converted into quadratic
unconstrained binary optimization (QUBO) problems [20], which can be
efficiently solved with the Quantum Approximate Optimization Algorithm (QAOA)
[21] using gate-based universal quantum computers or using D-Wave quantum
annealers. In the literature, there exist multiple quantum computing
approaches towards unit commitment [22, 23, 24, 25] and other mixed integer
problems [26], using quantum-inspired ADMM [27] or Benders’ decomposition
methods [28]. Quantum annealing approaches are also used for community
detection in electrical grids [29], peer-to-peer energy trading [30] or
coalition structure optimization [31, 32].

As one of this work’s main contributions is developing a problem-specific
decomposition method to solve large instances of the DSP on currently
available hardware, we give a brief overview of combinatorial problem
decomposition algorithms in the context of quantum optimization here. Divide-
and-conquer approaches have been used for various problem instances, such as
the MaxClique problem [33, 34, 35, 36], Minimum Vertex Cover [37, 36],
Community Detection [38] and MaxCut [38, 39]. They all combine the splitting
of the problem into sub-problems using problem-related methods. In special
cases, such as Ref. [39], quantum optimization is further utilized in
recombining the solution because of the special $\mathbb{Z}_{2}$ symmetry of
MaxCut solutions. Quantum Local Search (QLS) [40] takes local sub-problems of
a graph-based problem and iteratively improves a global solution vector.
Although applicable to any graph-based problem, QLS has been specifically
tested for the Maximum Independent Set problem. The recent emergence of
distributed quantum computing has led to the development of decomposition
algorithms that still allow for a limited amount of quantum information
exchange between the optimization of the sub-problems [41, 42], which was
successfully demonstrated for the Maximum Independent Set problem. Apart from
problem-specific methods, general QUBO decomposition methods have been
devised, like QBSolv [43]. Here, subsets of binary variables of the full QUBO
are selected as sub-problems, which are solved on a quantum annealer, while in
parallel, a classical heuristic optimizes the original problem. During the
process, solutions to the sub-problems will incrementally improve the current
solution state of the heuristic.

###  II-B Introduction to Quantum Annealing

Quantum annealing (QA) is a heuristic for solving combinatorial optimization
problems, first proposed in 1998 by Kadowaki and Nishimori [9]. QA utilizes
the adiabatic theorem to find the unknown ground state of an Ising Hamiltonian
$\mathcal{H}_{\text{Ising}}$, whose minimal energy state corresponds to the
solution of a target problem.

With $\mathcal{H}_{\text{Init}}$ being the initial Hamiltonian, the annealing
process can be described by the following dynamic Hamiltonian:

$\displaystyle\mathcal{H}(s)=A(s)\mathcal{H}_{\text{Init}}+B(s)\mathcal{H}_{\text{Ising}}$
(1) $\displaystyle\quad\mathcal{H}_{\text{Init}}=-\sum_{i}{\sigma}_{x}^{i}$
(2)
$\displaystyle\quad\mathcal{H}_{\text{Ising}}=-\sum_{i}h_{i}{\sigma}_{z}^{i}-\sum_{i>j}J_{ij}{\sigma}_{z}^{i}{\sigma}_{z}^{j},$
(3)

where ${\sigma}_{x,z}^{(i)}$ are Pauli matrices operating on qubit $i$, and
$h_{i}$ and $J_{i,j}$ are the qubit biases and coupling strengths, which
encode the specific problem. $A(s)$ and $B(s)$ are known as the annealing
schedule, with $s\in[0,1]$. At $s=0$, $A(s)\gg B(s)$, while $A(s)\ll B(s)$ for
$s=1$. As we increase $s$ from 0 to 1, the system undergoes a gradual change
from $\mathcal{H}_{\text{Init}}$ to $\mathcal{H}_{\text{Ising}}$. The
adiabatic theorem of quantum mechanics states that if that evolution happens
slowly enough and the system is initialized in the trivial ground state of
$\mathcal{H}_{\text{init}}$, then the state will remain in the ground state of
the momentary Hamiltonian [44]. Eventually, at $s=1$, the state will be in the
ground state of the $\mathcal{H}_{\text{Ising}}$. Finding the ground state of
the Ising model is isomorphic to QUBO [20], therefore, measuring the final
state will reveal the solution to an NP-hard optimization task.

In quantum annealing, this transition speed will typically be faster than
required for the adiabatic theorem, due to practical considerations.
Nevertheless, experimental evidence suggests that, depending on the spin glass
model, faster evolution times still output the solution with high probability
[45]. Thus, measuring the output repeatedly will eventually find the correct
solution.

##  III Problem Formulation

Since currently the distribution system operator (DSO) cannot automatically
influence the consumption of the customer at a certain time yet, we have to go
the leeway over price incentives. We assume a customer is strictly
economically motivated, i.e., he tries to minimize his electricity costs.
Thus, we aim to assign each consumer a custom price that dictates how much the
respective consumer varies its load. Of course, the convenience of having
access to electricity at all times is more important than saving on the cost,
such that, in reality, a customer cannot vary his consumption arbitrarily at a
certain time. However, with the emergence of electric vehicles (EV) with home
charging and heat pumps, automatically varying the load comes in the realm of
the possible. The given discounts (and penalties) serve as a protocol that
tells a smart home appliance on the customer side when to use electricity and
when not, e.g., start or stop charging the EV.

The goal on the DSO side is to reduce its CO2 emissions, which is linked to
production cost through carbon pricing. We, therefore, aim to provide
customers with individual dynamic tariffs in order to reduce the overall CO2
emissions. However, we strive to avoid change in the total consumed energy by
a single customer, since we only want to shift the time of consumption.
Beneficially, reducing the CO2 emissions leads to more consumption at times
with lots of local renewable energy production.

To formulate the problem, we discretize the optimization horizon in $N_{t}$
steps and assign each customer $c=1,\dotsc,N_{c}$ a discount (or penalty)
$z_{c,t}$ at each timestep. Furthermore, we have the forecasted consumption
data $d_{c,t}\geq 0\,[\text{kWh}]$ for each customer and predicted grid CO2
intensity $I_{t}\,[\text{g}/\text{kWh}]$ of the power generation in the
considered region. We start by introducing a suitable variable encoding for
the discounts and define the main optimization objective afterward. Finally,
we present grid and customer constraints, as well as lower-priority
optimization objectives.

###  III-A Discount Encoding

The given discounts $z_{c,t}$ are defined as discrete discount categories for
two reasons. First, discrete discounts are easily represented through integer
encoding as binary variables [20]. This makes it easy to translate the
following formulation into QUBO form, which is needed to employ quantum
optimization techniques. Secondly, discrete steps allow the user to change his
behavior more distinctly. For instance, if a thousand customers would receive
a tiny discount, we cannot expect all the thousand customers to increase their
consumption by just a tiny bit. Instead, suppose we offer a moderate discount
to only a few customers. In that case, we can expect those customers to adapt
their load such that the overall demanded effect of increasing consumption by
a small amount is achieved. Thus, limiting the discounts to a small set of
categories is useful.

In order to prepare the optimization formulation for quantum computation, we
need a full binary encoding of the variables. We have two options for that:
First, by assigning each discount category $Z$ ($z_{c,t}\in Z$) its own binary
variable. In a valid encoding, exactly one of the introduced binary variables
has to be set to 1, while the others must be 0. This so-called _one-hot
encoding_ has the advantage that we do not have to impose structure on $Z$,
and the effects of a given $z_{c,t}$ can also be non-linear. Nevertheless,
one-hot encoding requires an additional penalty term that enforces valid
bitstrings as solutions. Based on our initial experiments, we choose _integer
encoding_ as the preferred encoding since it provided better results as no
invalid solutions are possible. Here, we discretize a range
$[z_{\text{min}},z_{\text{max}}]$ into $N_{k}$ linearly spaced categories. We
choose a symmetric interval $z_{\text{max}}=-z_{\text{min}}=z_{m}>0$.
Therefore, $Z=\\{-z_{m}+i\Delta z\,|\,i=0,\dotsc\,N_{k}-1\\}$, with $\Delta
z=\frac{2z_{m}}{N_{k}-1}$. A $z_{k}<0$ refers to a discount, while $z_{k}>0$
is a penalty. This range can subsequently be expressed using
$Q=\lfloor\log_{2}N_{k}+1\rfloor$ binary variables $x_{c,t,k}$ for each
discount $z_{c,t}$

$\displaystyle z_{c,t}=\Delta{z}\sum_{k=0}^{Q-1}w_{k}x_{c,t,k}-z_{m},$ (4)
$\displaystyle\text{with }w_{k}=\begin{cases}2^{k}&\text{if }k<Q-1\\\
N_{k}-2^{Q-1}+1&\text{else}\end{cases}.$ (5)

Expressing the discounts this way, every bit combination $x_{c,t,k}$ results
in a valid encoding. Thus, there is no need for an additional penalty term in
the objective. Furthermore, this encoding is more space efficient, allowing
for an exponential number of categories to be represented with a linearly
growing number of qubits. Although typically an advantage, the impact on the
discussed problem is relatively small as the number of discount categories is
chosen small.

Given the customer initially receives a flat tariff $t_{0}$
$[\text{\texteuro}/\text{kWh}]$, the discount of penalty only affects the
consumption that deviates from the projected consumption $d_{c,t}$. Say the
new consumption is $\tilde{d}_{c,t}$, the customer pays

$\displaystyle
p_{c,t}=t_{0}d_{c,t}+(1+z_{c,t})t_{0}(\tilde{d}_{c,t}-d_{c,t}).$ (6)

A consumer will pay less if $z$ is negative and he increases his consumption
or when $z$ is positive and he decreases it.

###  III-B Objective

The main goal of the optimization problems is to lower the overall CO2
emissions of all customers by shifting the load to times when the grid CO2
intensity $I_{t}$ is small. We are not interested in decreasing the overall
consumption, but in shifting the load to an alternate timeslot. Furthermore,
let us introduce the average adjusted CO2 grid intensity $\Delta
I_{t}=I_{t}-\langle I_{t}\rangle_{t}$, which is negative whenever carbon
emissions are low and positive otherwise. $\langle I_{t}\rangle_{t}$ is the
average CO2 intensity.

The central assumption of the formulation is that a given discount (or
penalty) influences the customer’s behavior to alter their consumption
proportionally. This follows the inverse of the effect on the price change
from Eq. (6), such that the new consumption becomes

$\displaystyle\tilde{d}_{c,t}=(1-\chi_{c}z_{c,t})d_{c,t},$ (7)

when given a discount $z_{c,t}$. Here, $\chi_{c}>0$ is a scaling constant that
states how intensely a customer will respond to an incentive.

Finally, we can formulate the cost function for our optimization
$\operatorname*{arg\,min}_{z}R(z)$ as follows

$\displaystyle R(z)=\sum_{c,t}\Delta I_{t}[1-\chi_{c}z_{c,t}]d_{c,t},$ (8)

where $z(x)\in Z^{N_{c}\times N_{t}}$ is the discount matrix that can be
encoded through binary variables $x\in\mathbb{B}^{N_{c}\times N_{t}\times
N_{k}}$. As we later include soft constraints into the objective formulation,
we rescale and shift the objective (without changing the solution to the
problem) by employing the constants $R_{\text{min}}$ and $R_{\text{max}}$

$\displaystyle\text{obj}(z)=\frac{R(z)-R_{\text{min}}}{R_{\text{max}}-R_{\text{min}}}.$
(9)

Here, $R_{\text{min}}$ serves as a lower bound to the amount of CO2 that can
be reduced through our optimization, and the
$R_{\text{max}}=R(0)=\sum_{t}\Delta I_{t}D_{t}$ is the result of $R(z)$ prior
to optimization, i.e., $z_{c,t}=0$. Therefore, $\text{obj}(z)\in[0,1]$. Note
that valid solutions with $R(z)>R_{\text{max}}$ do exists, but since the
optimization should optimally take place $R(z)<R_{\text{max}}$, i.e., below
the naive baseline, we are satisfied with this scaling.

A naive way to compute the lower bound $R_{\text{min}}$ is to give maximal
discounts (penalty) when CO2 intensity is low (high).

$\displaystyle R_{\text{min}}=\sum_{c,t}\Delta
I_{t}[1-\chi_{c}\operatorname{sign}(\Delta I_{t})z_{m}]d_{c,t}.$ (10)

Of course, this lower-bound solution does not satisfy the constraints
introduced in Sec. III-C. Therefore, it is substantially smaller than the best
solution. Nevertheless, in Sec. IV-A1, we present a better method of finding a
more accurate lower bound.

###  III-C Constraints

The following constraints must be satisfied for $z$ being an applicable
discount matrix that can be forwarded to the customers. We decide between two
kinds of constraints: Hard constraints that must never be violated and soft
constraints, where small violations are tolerated but are aimed to be
minimized. This typically happens by adding the violation as a lower-priority
term to the objective function. An additional scaling constant then dials in
the strength of the violation penalty.

####  III-C1 Consumption deviation constraint

We do not want customers to change their total consumed energy over the
optimization horizon, i.e.

$\displaystyle\delta_{c}=\frac{1}{D_{c}}\sum_{t}z_{c,t}d_{c,t}=0\quad\forall
c\in\\{1,\dots,N_{c}\\},$ (11)

where $D_{c}=\sum_{t}d_{c,t}$ is the total energy usage of a certain customer.
Yet, this strict equality can generally not be exactly achieved through
discrete discounts $z_{c,t}$ on continuous data $d_{c,t}$, unless $z_{c,t}=0$.
Therefore, we impose a soft constraint

$\displaystyle\text{obj}(z)+\frac{\lambda}{N_{c}z_{m}^{2}}\sum_{c}\delta_{c}^{2}$
(12)

with a penalty factor $\lambda$ that steers the importance of the constraint.
Furthermore, the squared error is normalized to make its action problem size
independent.

####  III-C2 Global consumption deviation constraint

Even though the per-customer consumption deviation is soft-constrained, the
consumption deviation of all customers together may deviate. Globally, i.e.,
the combined view of all customers, we do not want any change in overall
consumption

$\displaystyle\frac{1}{D}\sum_{c,t}z_{c,t}d_{c,t}\chi_{c}=\frac{1}{D}\sum_{c}D_{c}\delta_{c}=0.$
(13)

Normally, this equality can be held up to numerical precision when many
customers are considered.

####  III-C3 Power restriction constraint

The momentary change in consumption of all customers should be bounded:

$\displaystyle p^{\text{low}}_{t}\leq\sum_{c}\chi_{c}z_{c,t}d_{c,t}\leq
p^{\text{high}}_{t}\quad\forall t\in\\{1,\dots,N_{t}\\}.$ (14)

This should be hard-constrained since we want to allow a range of consumption
deviation through scheduling. Physical grid limitations are the reason for the
limited flexibility here. The values for $p_{t}^{\text{high/low}}$ can be
determined using power flow computations. Of course, the presented global
restriction is a major simplification, but it suffices for an initial
investigation of the problem.

####  III-C4 Discount change penalty

We want to give the customer long periods of the same discount instead of
rapidly changing discounts in order to give them response time to the
discount. Therefore, we penalize the changes of discounts between two
consecutive timesteps:

$\displaystyle\text{obj}(z)+\frac{\tau}{4N_{c}(N_{t}-1)z_{m}^{2}}\sum_{c,t}(z_{c,t}-z_{c,t+1})^{2}\,.$
(15)

The normalized penalty can be dialed in using the factor $\tau$, which will be
small.

####  III-C5 Discount regularization

We do not want to give out discounts that only have little effect on the
objective. For example, suppose somebody does not consume any electricity at
some timestep. In that case, a discount won’t have an impact but may be given
anyways since it does not make a difference. An $L2$-regularization ensures
that discounts are only given if they benefit the overall goal to a certain
extent

$\displaystyle\text{obj}(z)+\frac{\rho}{N_{c}N_{t}z_{m}^{2}}\sum_{c,t}z_{c,t}^{2}\,.$
(16)

Again, the penalty is normalized and a penalty factor $\rho$ is used, which
has low priority.

###  III-D On customer savings

Discount scheduling incentivizes customers to reduce their CO2 through
discounts and penalties. Yet, the procedure should never happen to the
customer’s disadvantage, even if they do not change their behavior.

The customer’s cost change over the optimization horizon can be computed via
the sum of momentary price differences through Eq. (6) and Eq. (7):

$\displaystyle\Delta p_{c}$
$\displaystyle=\sum_{t}p_{c,t}-t_{0}\sum_{t}\tilde{d}_{c,t}$ (17)
$\displaystyle=t_{0}\sum_{t}z_{c,t}(\tilde{d}_{c,t}-d_{c,t}).$ (18)

Note that we have used the sum over the changed consumption as the baseline
for our comparison, since in any case
$\sum_{t}d_{c,t}\approx\sum_{t}\tilde{d}_{c,t}$ and we only want to compare
the cost for the same amount of purchased electricity.

Plugging in the altered consumption from Eq. (7), we get as a price change

$\displaystyle\Delta p_{c}=-t_{0}\chi_{c}\sum_{t}z^{2}_{c,t}d_{c,t}.$ (19)

As $z_{c,t}^{2}\geq 0$ and it is assumed that $\chi_{c}\geq 0$, the customer’s
price change is guaranteed to be $\Delta p_{c}\leq 0$, so a customer will save
money by complying to the incentives. The savings are exactly zero if the
customer does not change his behavior at all. As the absolute price change
quantity is dependent on the flat tariff and the total consumption of the
customer, we will look at the relative savings $s_{c}=-\Delta
p_{c}/\sum_{t}t_{0}\tilde{d}_{c,t}\geq 0$ in the evaluation section.

##  IV Problem Decomposition

The number of integer variables needed to construct the discount matrix is
$N_{c}\times N_{t}$. Given a one-day optimization horizon with 15-minute
timesteps, each customer requires 96 integer decision variables in the
problem. However, as the number of customers will grow quite
large111Typically, we want to consider more than 1000 customers., the number
of integers grows akin. Even worse, the number of qubits in the quantum
formulation is scarce, and every integer must be encoded with $\lfloor\log
N_{k}+1\rfloor$ qubits. Thus, the move to a hybrid quantum-classical
optimization scheme seems inevitable.

In this section, we propose a hybrid approach that is based on problem
decomposition. Despite the drawback that decomposition increases solution
bias, we find that we can manage the hard constraints of the DSP classically
in a pre-processing step. This eliminates the need for a costly reformulation
of inequality constraints with slack variables. Fig. 1 shows an overview of
the steps taken for the decomposition.

###  IV-A Motivation

####  IV-A1 Global Solution

Let us shift our perspective from the individual customer level to a global
view, where all customers are regarded as a unified entity. We consider the
overall consumption $D_{t}=\sum_{c}d_{c,t}$ and the mutable consumption, i.e.,
the consumption weighted by the individual customer susceptibilities
$\widetilde{D}_{t}=\sum_{c}\chi_{c}d_{c,t}$. Furthermore, we can express the
weighted average of all discounts given per customer—from now on called
effective discount—as follows

$\displaystyle\zeta_{t}=\langle
z_{c,t}\rangle_{c}=\frac{1}{\widetilde{D}_{t}}\sum_{c}\chi_{c}d_{c,t}z_{c,t}\in[-z_{m},z_{m}].$
(20)

Utilizing the formulation of the effective discount, we can transform the CO2
reduction from Eq. (8) into

$\displaystyle R(\zeta)=\sum_{t}\Delta
I_{t}\left(D_{t}-\widetilde{D}_{t}\zeta_{t}\right).$ (21)

The global consumption deviation constraint, Eq. (13), and the power
restriction constraint, Eq. (14) can be expressed solely in terms of the
effective discount. Therefore, we represent the global version of the DSP as a
linear program

$\displaystyle\begin{split}\text{minimize}\quad&R(\zeta)\\\ \text{subject
to}\quad&p_{t}^{\text{low}}\leq\widetilde{D}_{t}\zeta_{t}\leq
p_{t}^{\text{high}}\quad\forall t\in{1,\dots,N_{t}}\\\
&\sum_{t}\widetilde{D}_{t}\zeta_{t}=0.\end{split}$ (22)

This formulation disregards any per-customer constraints that are still part
of the DSP. Nevertheless, it is a useful tool to estimate how much CO2
reduction is maximally possible, with all the hard constraints in place. In
fact, the solution $\zeta_{t}^{*}$ is guaranteed to give an optimal lower
bound $R(\zeta^{*})$

$\displaystyle R(\zeta^{*})\leq R(z)\quad\forall z\in\mathcal{Z},$ (23)

where $\mathcal{Z}$ is the set of feasible discount matrix configurations. The
global DSP consists of only $N_{t}$ continuous variables. Thus, it can be
quickly and efficiently solved using standard procedures like the Simplex
method [46].

Given an optimal effective discount, $\zeta_{t}^{*}$, we can utilize Eq. (20)
to optimize the integers $z_{c,t}$ for the individual customers per timestep.
Additionally, we need to include the soft constraints from the DSP in the
subsequent optimization. However, doing so would yield an optimization problem
the same size as the original problem.

Nonetheless, the following section reveals that we can achieve a satisfactory
approximation of a continuous number by considering only a limited number of
customers. As a result, we can divide the customers into smaller groups or
chunks and optimize each chunk separately.

####  IV-A2 Representational Power

FullProblemSub-Problem $1$Sub-Problem $2$Sub-Problem $N$Sub-Solution $1$Sub-
Solution $2$Sub-Solution $N$$\vdots$$\vdots$Post-Processing$\cdots$ Figure 1:
Overview of the decomposition routine. The problem is split into sub-problems.
Solutions can influence the following sub-problems via sequential updating.
Finally, sub-solutions are gathered to a full solution and a post-processing
step is employed that improves the solution quality greedily while also making
the power restriction constraint is satisfied. Figure 2: The relative
approximation error for different values of $\zeta$ averaged over multiple
timesteps.

In this section, we motivate that Eq. (20) can be fulfilled for any arbitrary
$\zeta_{t}$ with sufficient accuracy given a small constant number of
customers. We will focus on a discount range $\zeta_{t}\in[-1/2,1/2]$ and five
discrete discounts $z_{c}\in\\{-1/2,-1/4,0,1/4,1/2\\}$. From the generated
consumption data, see Sec. V, we take a random set of customers and compute

$\displaystyle\min_{z_{c}}\frac{1}{\zeta}\left|\frac{1}{D}\sum_{c}d_{c,t}z_{c}-\zeta\right|$
(24)

for all available time steps. Fig. 2 shows the result with different numbers
of customers. The average over all timesteps is plotted, and the error bands
indicate a 95% confidence interval. It is evident that even with only ten
customers, the relative error remains consistently below 1%. As more customers
are added, the error decreases significantly, reaching a negligible level.
Therefore, we contend that by maintaining a small, constant number of
customers within a chunk (e.g., 20-50 customers), it is possible to obtain a
reliable approximation of an effective discount while still considering the
per-customer soft-constraints from the DSP.

###  IV-B The Full Decomposition Routine

Let us now assemble the pieces into a full hybrid routine for decomposition,
as seen in Fig. 1. The process begins with solving the global DSP (22),
followed by dividing customers into chunks. We sort the customers by total
consumption and split them into $M$-sized groups, s.t. the largest customers
are arranged in the first chunk, etc. We argue that it is better to have
customers with comparable consumption in one chunk because they can counteract
each other better than e.g. one industrial customer and 20 single households.
For each chunk, we can define sub-problems in which special effective
discounts per chunk are introduced in Sec. IV-B1. These sub-problems can be
transformed into QUBO format and solved on a QC.

Since we can solve the sub-problems sequentially, we can enhance the results
by incorporating the errors from prior optimizations into the subsequent sub-
problems in Sec. IV-B2. Eventually, all the chunks are collected, and a final
post-processing step shown in Sec. IV-B3 is applied to ensure that no
constraints are violated.

####  IV-B1 Chunk Problems

The customers are partitioned into $M=N_{c}/m$ mutually exclusive chunks
$C_{j}$, s.t. $\bigcup_{j}C_{j}=\\{1,\dots,N_{c}\\}$ and $C_{i}\cap
C_{j}=\emptyset\,\forall i\neq j$. Note, that we require and expect the chunk
size to be chosen, s.t. $N_{c}\mod m=0$.

Most likely the consumption deviation per chunk

$\displaystyle\sum_{c\in C_{j}}\sum_{t}\chi_{c}d_{c,t}\zeta^{*}_{t}\neq
0\quad\forall j$ (25)

is not zero, which, by default, violates the consumption deviation soft-
constraint (12). Thus, the first goal is to define chunk effective discounts
$\xi_{t}^{j}$ with the following properties:

$\displaystyle\sum_{t}\widetilde{D}_{t}^{j}\xi_{t}^{j}$
$\displaystyle=0\quad\forall j,$ (26)
$\displaystyle\frac{1}{\widetilde{D}_{t}}\sum_{j=1}^{M}\widetilde{D}^{j}_{t}\xi_{t}^{j}$
$\displaystyle=\zeta^{*}_{t}\quad\forall t,$ (27)

where we define an alterable consumption for one chunk
$\widetilde{D}^{j}_{t}=\sum_{c\in C_{j}}\chi_{c}d_{c,t}$, similar to the
definition of the total alterable consumption.

We define the chunk-effective discount as follows

$\displaystyle\xi_{t}^{j}=\zeta^{*}_{t}-\frac{\alpha_{t}}{\widetilde{D}_{t}^{j}}\sum_{t^{\prime}}\widetilde{D}_{t^{\prime}}^{j}\zeta_{t^{\prime}},$
(28)

where $\alpha_{t}$ are arbitrarily chosen constants, s.t.
$\sum_{t}\alpha_{t}=1$. The conditions (26) and (27) are satisfied with this
definition. The values $\alpha_{t}$ are chosen constant $\alpha_{t}=1/N_{t}$,
but we have to make sure that $\xi_{t}^{j}\in[-z_{m},z_{m}]\,\forall t,j$. If
this is not possible for one timestep $t$, we have to dial back the
$\alpha_{t}$ while equally increasing the remaining $\alpha$s.

The optimization objective is to approximate the following equality with the
chunk effective discount as well as possible

$\displaystyle\widetilde{D}^{j}_{t}\xi^{j}_{t}=\sum_{c\in
C_{j}}d_{c,t}\chi_{c}z_{c,t}\quad\forall t\in\\{1,\dots,N_{t}\\}.$ (29)

The objective can be reformulated as a least squares error problem

$\displaystyle
z^{*}_{c,t}=\operatorname*{arg\,min}_{z_{c,t}}\frac{1}{N_{t}z_{m}^{2}}\sum_{t}\left(\xi_{t}-\frac{1}{\widetilde{D}^{j}_{t}}\sum_{c\in
C_{j}}d_{c,t}\chi_{c}z_{c,t}\right)^{2}$ (30)

and is directly in QUBO form after the binary representation of the discounts
has been plugged in.

The previously discussed constraints and regularizations—consumption deviation
(12), discount regularization (16) and discount change penalty (15)—can be
carried over to this optimization problem.

####  IV-B2 Sequential updating

When the sub-problems are solved in sequence, the error between the true
achieved effective discount and the demanded one can be carried over into the
next optimization to be corrected. For optimizing $\xi_{t}^{j}$ can be adapted
as follows

$\displaystyle\xi^{j}_{t}\leftarrow\xi^{j}_{t}+\frac{1}{\widetilde{D}_{t}^{j}}\sum_{i=1}^{j-1}\left(\widetilde{D}_{t}^{i}\xi^{i}_{t}-\sum_{c\in
C_{i}}z^{*}_{c,t}d_{c,t}\right).$

Doing so will significantly improve the overall accuracy of the method. Of
course, one has to ensure that the altered $\xi$s do not exceed the bounds
$[-z_{m},z_{m}]$.

####  IV-B3 Post-processing

Finally, we describe a post-processing scheme that refines the result and
ensures that the power restriction constraint (14) is held. The Algorithm 1
describes the greedy improvement of the solution.

Data: $z_{c,t},r\in\mathbb{N}\,\,\,$

# $r$ is a parameter that dials the accuracy/runtime

Result: $z_{c,t}$

$\Delta z\leftarrow 2\,z_{m}/(N_{k}-1)$;

# Discount step

$\Delta_{c,t}\leftarrow\chi_{c}d_{c,t}\Delta z$;

# Possible deviations

$\delta_{c}\leftarrow\sum_{t}\chi_{c}d_{c,t}z_{c,t}$;

# Consumption deviation

for  _$t\in\\{1,\dots,N_{t}\\}$_ do

$p\leftarrow\sum\chi_{c}d_{c,t}z_{c,t}$;

# Power deviation

$\varepsilon\leftarrow\zeta_{t}^{*}\widetilde{D}_{t}-p$;

# Error from demanded

# Increasible customers

$C^{\uparrow}\leftarrow\\{c=1,\dots,N_{c}\,|\,z_{c,t}<z_{m},\delta_{c}<-\Delta_{c,t}/2\\}$;

# Decreasible customers

$C^{\downarrow}\leftarrow\\{c=1,\dots,N_{c}\,|\,z_{c,t}>-z_{m},\delta_{c}>\Delta_{c,t}/2\\}$;

$C^{\uparrow}\leftarrow\text{limit}(C^{\uparrow},r)$;

# reduce size $|C^{\uparrow}|=r$

$C^{\downarrow}\leftarrow\text{limit}(C^{\downarrow},r)$;

# Compute combinations of increasing and decreasing two customer discounts

$X_{c,c^{\prime}}\leftarrow\operatorname{sign}(\zeta^{*}_{t})(\varepsilon-[\Delta_{c,t}-\Delta_{c^{\prime},t}])$;

# Find positive (feasible) ones

$C_{2}\leftarrow\\{(c^{\uparrow},c^{\downarrow})\in C^{\uparrow}\times
C^{\downarrow}\,|\,X_{c^{\uparrow},c^{\downarrow}}>0\\}$;

# Get the best move

$c^{\uparrow},c^{\downarrow}\leftarrow\operatorname*{arg\,min}_{c^{\uparrow},c^{\downarrow}\in
C_{2}}X_{c^{\uparrow},c^{\downarrow}}$;

if  _$X_{c^{\uparrow},c^{\downarrow}}
>\operatorname{sign}(\zeta^{*}_{t})\varepsilon$_ then

continue __;

# No improvement

end if

# Update solution and consumption deviation

$z_{c^{\uparrow},t}\leftarrow z_{c^{\uparrow},t}+\Delta z$;

$z_{c^{\downarrow},t}\leftarrow z_{c^{\downarrow},t}-\Delta z$;

$\delta_{c^{\uparrow}}\leftarrow\delta_{c^{\uparrow}}+\Delta_{c^{\uparrow},t}$;

$\delta_{c^{\downarrow}}\leftarrow\delta_{c^{\downarrow}}-\Delta_{c^{\downarrow},t}$;

end for

Algorithm 1 The post-processing algorithm

Conceptually, it is quite simple: For each timestep, we extract those
customers whose discounts can be increased or decreased while also improving
the consumption deviation penalty (12). Then we try all combinations between
one increase and one decrease and investigate how the effective discount
behaves. If $\zeta^{*}_{t}$ is negative, we want the real effective discount
to be as close as possible but at least larger than $\zeta^{*}_{t}$. If it is
positive, the other way around. Doing so always satisfies constraint (14). We
find the combination that matches the requirements the best and update the
respective discounts if it achieves an improvement. Otherwise, the timestep is
skipped.

Since all possible combinations of up and down moves have to be considered,
the complexity of the Algorithm scales at worst with
$\mathcal{O}(N_{t}N_{c}^{2}/4)$. Nevertheless, limiting the possible moves to
at most $r$ empirically provides sufficient accuracy. This then reduces the
complexity to $\mathcal{O}(N_{t}N_{c}+N_{t}r^{2})$.

##  V Experiments & Results

###  V-A Experimental Setup

To benchmark the performance of solving the DSP, we consider out-of-the-box
solvers, as well as our developed decomposition method, and evaluate the
results on a set of metrics that best represent the different goals described
in the DSP formulation.

####  V-A1 Investigated Solvers

An overview of the considered solvers and settings can be found in Table I. As
a state-of-the-art purely classical baseline, we used Gurobi222All experiments
with Gurobi were conducted on an M1 MacBook Pro (2020) with Gurobi Version
9.0. [47]. This was compared to D-Wave’s LeapHybridCQM solver [11] (called
just Leap in the following), which is a quantum-classical hybrid algorithm
that uses classical algorithms to optimize the problem whilst using quantum
computers to solve suitable sub-tasks. This has the benefit of solving larger
problems than possible directly on current quantum hardware while also
supporting more sophisticated optimization models that include hard
constraints. Leap is accessed through D-Wave’s Cloud service. These two out-
of-the-box solvers are compared against our own problem-specific decomposition
routine introduced in Sec. IV, subsequently called Decomp-Gurobi, Decomp-Leap
or Decomp-QPU, depending on the method considered for solving the chunk
problems (30). QPU refers to direct access to the D-Wave’s Quantum Annealing
processor Advantage 4.1 [48]. Whenever a decomposition solver is followed by
an integer, it refers to the chunk size $m$. The post-processing algorithm is
turned on and a cut-off value at $r=500$ is chosen.

To ensure a fair comparison, we gave each solver a time limit of
$0.1\,\mathrm{s}\times N_{c}$. The timeout needs scaling with the problem size
since the problem difficulty grows considerably with the problem size.
Nevertheless, we observed that Leap tends to overrun the set timeout, which is
the reason that we first run Leap with the linear growing timeout and then run
the remaining solvers with the timeout set to Leaps runtime. Since the
Decomposition solver consists of multiple sub-solver calls, we set the timeout
for each sub-solver as the whole timeout divided by the number of chunks,
i.e., a timeout of $0.1\,\mathrm{s}\times m$.

TABLE I: Overview of Solvers Name | Definition
---|---
Gurobi | Classical Solver
Leap | Quantum classical hybrid solver
Decomp-Gurobi | Decomposition with Gurobi as sub-solver
Decomp-Leap | Decomposition with Leap as sub-solver
Decomp-QPU | Decomposition with D-Waves Advantage 4.1
| quantum annealing processor as sub-solver
Decomp-SA | Decomposition with simulated annealing as
| sub-solver

####  V-A2 Metrics

Because we are considering an optimization task with multiple goals involved,
it does not suffice only to consider the objective value of our model as a
performance metric. Instead, we simultaneously investigate multiple metrics:

  * •

_CO 2 reduction_: First and foremost, the CO2 reduction is the central goal of
the DSP, hence it is also the main metric that is investigated. We compute the
relative CO2 reduction error by making use of the solution to the global DSP.
Therefore,

$\displaystyle\frac{R(\zeta^{*})-R(z)}{R(\zeta^{*})},$ (31)

is a positive quantity and tells us how good the optimization has performed,
in comparison to a theoretical maximal reduction.

  * •

_Energy_ : The energy, or objective, of the optimization problem consists of
the rescaled CO2 reduction with the penalties added. For easier comparison, we
again utilize the relative energy error for investigating a solver’s
performance. The baseline is taken from the best possible CO2 reduction; all
penalties are set to zero. This is a guaranteed lower bound to the energy.

  * •

_Consumption deviation standard deviation_ : We expect the consumption
deviations for each customer to be centered around zero since the problem is
constrained to have a zero total consumption deviation. Therefore, the spread
of $\delta_{c}$ around the zero may be a good measure to judge, whether a
result produces satisfactory little consumption deviations. That is, we
measure the standard deviation of the consumption deviations $\delta_{c}$.

  * •

_Average discount changes_ : Since we want to reduce the changes between two
discount categories as much as possible, we measure the average discount
changes.

$\displaystyle\frac{1}{N_{c}(N_{t}-1)}\sum_{c}\sum_{t=1}^{N_{t}-1}(1-\delta_{z_{c,t},z_{c,t+1}}).$
(32)

Here, $\delta$ refers to the Kronecker-Delta.

  * •

_Average relative cost savings_ : Not a quantity that is optimized for, but
very interesting for the DSO, is to measure the relative cost savings per
customer, as defined in Sec. III-D. To obtain a single indicator of the
performance, we evaluate the mean $\langle s_{c}\rangle_{c}$ of the relative
savings.

####  V-A3 Data Generation

To facilitate the benchmarking of the DSP, we first require realistic data for
the consumption of individual customers. For that, we take standard load
profiles of residential and industrial customers, to which we add noise.
Furthermore, we randomly shift the load profiles by small amounts in time.
Additionally, they get scaled according to various numbers of inhabitants of a
household. The number of inhabitants is taken from a residential area in
Munich. Moreover, we include photovoltaic (PV) electricity generation into the
mix, by estimating the potential based on roof data and simulating the
production from historic solar irradiance data. The PV production screens the
customer’s consumption. Grid infeed, i.e., if more PV is generated than
consumed, is not specially considered. The grid CO2 intensity is taken from
the real-world data in Munich333The data is provided by E.ON’s App for
monitoring local CO2 intensities: https://www.bayernwerk.de/de/fuer-
zuhause/oekoheld.html.

####  V-A4 Parameters

TABLE II: Parameter setting for the investigated problems Parameter | Value
---|---
$\chi_{c}$ | 1
$N_{k}$ | 5
$z_{m}$ | $1/2$
$p^{\mathrm{high}}_{t}$ | $0.1\times\frac{1}{N_{t}}\sum_{t}D_{t}$
$p^{\mathrm{low}}_{t}$ | $-p^{\mathrm{high}}_{t}$
$\lambda$ | $0.1$
$\rho$ | $10^{-5}$
$\tau$ | $10^{-4}$

If we want to solve the DSP for a given data frame, consisting of the
consumption of $N_{c}$ customers at $N_{t}$ time steps, we still need to fix a
set of open variables and parameters. In a real-world scenario, the customer
susceptibilities $\chi_{c}$ would be measured from the individual customer’s
behavior. However, as it only acts as a proportionality constant, we turn
their effect off and set them all to one. Next, we use five discount
categories, with a 50% discount maximally. That, in turn refers to the
following valid discounts $z_{c,t}\in\set{-1/2,-1/4,0,1/4,1/2}$. As a
consequence, a discount of, e.g., 50% would result in an increase in the
customer’s consumption by 50%.

The power deviation bounds $p_{t}^{\mathrm{high/low}}$ are set to a constant
10% of the average total consumption. Due to our goal not being an overly
accurate representation of the real world and more the analysis of the quantum
solver, we prefer the constant values here. In practice, however, those values
may be derived from real-world grid constraints that can be inferred through
power-flow calculations.

Finally, the remaining penalty parameters are fixed by analyzing a small-scale
example with Gurobi and dialing in the strengths of the penalties, such that
they have a reasonable effect for the Gurobi result. It is important to note
that a comprehensive examination of the solver’s response to parameter
settings is beyond the scope of the current investigation.

An overview of all parameter settings is given in Table II.

###  V-B Example with 100 Customers

Figure 3: The discount matrices $z_{c,t}$ found by the investigated solvers
for $N_{c}=100$. Blue indicates a discount and red corresponds to a penalty.
White means no discount given at all. Despite their effects on the overall
consumption (see Fig. 4) being the same, the discount matrices differ a lot
from each other. It is apparent that Gurobi hands out the discounts more
greedily than Leap, indicating a bigger impact of the regularization.
Nevertheless, a similar pattern is observable in the last three solutions.
Figure 4: The effect of the DSP solution for problem size $N_{c}=100$. The
plot shows the aggregated consumption with and without (Base) discounts in
place, as well as the grid CO2 intensity. The solutions of all solvers produce
a similar effective consumption change, as already predicted in Sec. IV. Times
with high CO2 emissions produce an effective decrease in consumption and vice
versa, just as expected.

Let us first take a look at how the optimization result of the different
solvers looks in detail before we only focus on the previously discussed
metrics. To this end, we take 100 random customers out of the full dataset
consisting of over 16,000 simulated customers. The 76 timesteps reach from 1
am to 7 pm for January 13, 2023. We analyze the solutions of four solvers,
Gurobi, Leap and two $m=50$ decomposition methods with the same solvers as the
sub-routine. The results for the discount matrices $z_{c,t}$ can be seen in
Fig. 3, while their overall effect on the consumption is displayed in Fig. 4.
Although the particular solutions differ quite a lot, the effective result
stays similar, regarding the CO2 reduction. Especially the difference between
Gurobi and the other solvers is notable.

Figure 5: Histogram of the relative consumption deviation. One can see that
both Gurobi solvers have relatively little spread. And are well centered
around zero. The Leap solvers, on the other hand, possess a large spread and
are additionally shifted away from zero. The shift away from zero reduces in
larger problem instances. Figure 6: A cumulative distribution plot of the
relative savings of the customers. In the Decomp solvers, the two chunks can
be well distinguished. Gurobi, because it is greedy with handing discounts,
only distributes savings to relatively few customers. On the other hand, Leap
distributes similar discounts to all customers. Remember: We do not optimize
for this metric. This is just an observation of the different strategies and
can be interpreted as a measure of the fairness of the optimization algorithms
with respect to different customers.

Apart from the global action of optimization, we are also interested in how
the optimization performs per customer. In Fig. 5, one can see how the final
relative consumption changes are distributed. Furthermore, Fig. 6 visualizes
the distribution of cost savings to the customers.

Lastly, it remains important to note that the results for the Leap solvers
vary throughout multiple runs. Here only a single run has been picked that is
characteristic of the behavior of these solvers. Furthermore, no investigation
towards direct QPU access has been made, since the space requirements for a
single customer are already 76 integer variables, i.e., 228 binary variables.
The problem after gathering multiple customers in a chunk is, hence, not
embeddable in the QPU, since we are facing quite dense connectivity in the
QUBO. For a reduced problem size, we perform investigations in Sec. V-F.

###  V-C Scaling Analysis

Figure 7: The investigated metrics for different problem sizes and different
solvers. The runtime of all solvers has been set to be equal for a certain
problem size, but grows with $N_{c}$. The top row shows the global metrics,
which tell the most about how the solver performed. The relative energy error
is the central objective that we try to minimize, while the energy error gives
an overview of the performance with regards to all optimization targets. The
second row shows the per-customer metrics, which are optimized for in the
problem formulation. That is the standard deviation of all customer
consumption deviations and the average discount changes, which we both want to
be small. The last row displays the average relative savings and the deviation
from the $0.1\,\mathrm{s}\times N_{c}$ runtime. The error bands indicate the
maximum and minimum of the three runs. The spike in relative runtime at 25
customers arises from the minimum Leap timeout of 5 s. As one can see, the
other solvers follow the time Leap took well.

To test the performance of different solvers, we created test instances using
the generated data described above with $N_{c}$ ranging from 25 to 3200
customers and considering the full 76 timesteps. Our problem instances,
therefore, consist of 1,900 to 243,200 integer variables444All solvers can
handle integers directly, so we do not perform the encoding explicitly, but
supply the solver with the full information.. To account for the stochasticity
of the results from the quantum solvers, we run the quantum solvers multiple
times. However, we need to limit ourselves to three execution of the quantum
routines because of cost considerations.

The results for the discussed metrics are visualized in Fig. 7. Each plot
shows a singular metric against the problem size for the considered solvers.
Focusing on the relative CO2 reduction, it is evident that a crossover in
performance between Gurobi and Leap happens between 100 and 200 customers.
After that size, Gurobi is not able to find converged results in the given
time limit. Although not a directly fair comparison since Gurobi is run on a
local machine whilst the Leap hybrid solver is run on an proprietary D-Wave
cloud architecture. However, we argue that the pattern generalizes, i.e., the
point where Gurobi doesn’t reach satisfactory results anymore shifts to the
right but eventually happens. Leap starts off with a relatively weak
performance in the small problem sizes but decreases its energy almost
continuously. Yet the decomposition routines greatly outperform the general-
purpose solvers. We can see steadily decreasing relative CO2 reduction errors,
which can indicate that the total error stays fixed, but gets diminishingly
small in comparison to the total quantity.

The relative energy error loosely shows the same picture, but there we can see
that the solvers utilizing Leap approach a similar level, which can be
explained through the dominating penalty terms of the other objectives, where
the decomposition does not interfere anymore. Decomp-Gurobi returns a very
good and constant optimization objective.

Moving on to the per-customer constraints, we notice that the Gurobi based
solvers outperform the quantum-enhanced routines (at least where they
converge). This is likely due to Gurobi being better at handling smaller
energy changes in the optimization objective. However, it is also important to
note that, as apparent from the discount matrices in Fig. 3, Gurobi gives many
customers not even a single discount. Hence they do not receive any discount
changes or consumption deviations, which reduces the average measure.
Furthermore, by investigating customer savings, it becomes clear that, on
average, Gurobi does provide discounts of well below 5% before becoming
extremely generous. The latter arises from the discount matrix being almost
completely filled with $z_{c,t}=\pm z_{m}$, and the customer savings are
dependent on a weighted average over
$z_{c,t}^{2}=\sum_{t}z_{m}^{2}=\mathrm{const.}$, remember (19). Therefore,
$s_{c}=z_{m}^{2}=1/4$ in the investigated case. Furthermore, which is also
visible in Fig. 3, the more generous customer savings of Decomp-Leap can be
explained the same way.

To conclude this analysis, we remark that Gurobi struggles with its solving
strategy for large problem sizes, which indicates a potential advantage of the
quantum-enhanced solver here. Yet, the biased domain-specific decomposition
routine provides even better results, especially with the classical solver
underneath. We argue, that since the decomposition-based solvers work so well,
the possible space of good solutions is rather big, which makes this problem a
fitting choice for heuristic-based solvers more than mathematical solvers,
like Gurobi.

###  V-D Chunk Size Effect

Figure 8: Per-customer metrics evaluated with different chunk sizes in the
decomposition. As expected, the metrics improve as the chunk sizes get larger
since more flexibility remains in the chunk. The global metric, i.e., the CO2
reduction, performs equally well for all chunk sizes.

After we have seen that the decomposition solver provides satisfying results
both with Gurobi and Leap employed as sub-solver, we are interested in what
impact the chunk size has on the result. For that, we only inspect Decomp-
Gurobi with different chunk sizes $m=5,10,25,50$ and focus on a reduced
problem size frame up until $N_{c}=800$. We have seen that the problem
complexity does not grow linearly with the problem size. Thus, we give a more
generous timeout of $0.5\,\mathrm{s}\times m$ in this investigation in order
to isolate the effects of the decomposition routine from the solver
performance555We observed a performance increase for the larger chunk sizes
when increasing the timeout. This increase came to a slowdown at around
$0.5\,\mathrm{s}\times m$.. The global effect, i.e., how much CO2 was reduced,
did not differ between the chunk sizes. They all performed equally well. The
constant sequential updating of the objective also helps a lot with finding
the best CO2 reduction, even with five customer chunks. Fig.8 shows the per-
customer metrics that are optimized for. Here, we can observe a clear tendency
that larger chunks result in better per-customer metrics, i.e., less
consumption deviation per customer and fewer overall discount changes.

###  V-E Fairness Analysis

Figure 9: A cumulative distribution plot of the relative savings of the
customers at $N_{c}=800$. As discussed earlier, Gurobi does not converge
anymore, which causes savings of around 25%. Leap produces fair discounts,
similar to Fig. 6. The other two solvers produce more complex, unfair savings
distributions.

The goal of this section is to investigate how the solvers strategically
distribute the discounts to the target customers. This is done by
investigating how the relative savings $s_{c}$ are distributed between
individual customers. Fig. 6 and Fig. 9 show two cumulative distribution plots
of the results from 100 and 800 customer problem sizes. Except for Gurobi,
which does not converge, the observable patterns of the solvers are similar.
Leap produces a fair savings distribution, which means that all customers
experience the same savings. In the cumulative plots, that means a straight
vertical line. The more that line is skewed, the more different cost savings
are registered between the customers.

In Fig. 6 the splitting in half of the decomposition can be observed quite
remarkably. The resolution of the 16 individual chunks in Fig. 9 is no more
possible. Although, a kink in Decomp-Leap can be observed, which means that
about 70% of the customers save a similar and relatively large amount, while
fewer savings are distributed to a smaller group. Decomp-Gurboi reveals a
straight but shallow curve, which means that customers will receive savings
between 0% and 20% almost equally likely.

###  V-F Direct QPU-Access with Decomposition

Figure 10: Embeddable sub-problem size for the D-Wave Advantage 4.1 QPU. The
left-hand matrix shows how many physical qubits are needed when a sub-problem
with $N_{c}$ customers and $N_{t}$ timesteps are embedded. A white field
indicates that no embedding has been found. The right-hand plot shows the
maximal chain length for the found embedding, i.e., how many qubits are
maximally connected to form one logical qubit. All embeddings were found using
D-Wave’s MinorMiner package.

A quantum annealing processor, such as D-Waves Advantage 4.1, suffers from
limited connectivity between the physical qubits. However, for our QUBO sub-
problems (30), we can analytically compute the number of couplings for a
single qubit as follows:

$\displaystyle N_{k}\left(N_{t}-1\right)+N_{k}\left(N_{c}-1\right)+N_{k}-1.$
(33)

This term is derived by inspecting the terms in the QUBO formula and observing
that we either have couplings within all customers at a single timestep or
couplings within all timesteps of a single customer. For the first case, one
qubit is connected to all $N_{k}$ qubits of the other $N_{c}-1$ customers and
to $N_{k}-1$ qubits of the same customer. Analogue for the second case, but
the $N_{k}-1$ connections within the timestep have already been covered in the
first case.

The derived quantity grows with the problem size, but the couplings per qubits
of the D-Waves Pegasus graph is a constant 15[48]. Thus, physical qubits have
to be chained together to logical qubits in order to allow for higher
connectivity. Finding the best, so-called embedding, is itself an NP-hard
optimization problem, for which we utilize D-Wave’s heuristic MinorMiner.

Figure 11: All metrics for the small problem sizes, also considering QUBO sub-
problems that can be solved using D-Wave’s QPU. Comparing both Quantum
routines, one can observe that Decomp-QPU returns slightly lower energy
results than Leap. Yet, the CO2 reduction is vastly better in the decomposed
case. Notably, the difference of the Decomp solvers in the relative energy
error is mostly due to the higher consumption deviation from the QPU results.
To match the runtime of the Decomp-SA to Decomp-QPU, we had to manually adjust
the number of reads from the simulated annealing routine to 30. We also
observe that Leap outruns the set timeout here, which we do not adjust since
the baseline is Decomp-QPU.

Fig. 10 shows the computed embeddings for the sub-problem QUBOs with different
problem sizes. It is apparent that we are very limited to small problem sizes.
Since we do not want too few customers in a chunk to preserve flexibility, we
settle at a reasonable middle ground of chunk size six and $12$ timesteps. We
interpolate the original data to 12 timesteps and use various (multiples of 6)
customer sizes to compare the performance of Decomp-QPU against the other
solvers. For each sub-problem, we take 100 readings from the QPU. We cannot
directly steer the timeout in this case. Thus, we chose to first run Decomp-
QPU and then set the timeout of the remaining solvers to exactly that time.
However, Leap has a minimum runtime of 5 s, which is the reason why we only
include Leap in the cases where the Decomp-QPU time is more than 5 s, being
the case from $N_{c}=480$ onwards.

Again, we perform the analysis for different problem sizes, reaching from 60
to 1920 customers, or 720 to 23,040 integer variables. The sub-problems
comprise 72 integer variables, resulting in 216 binary variables in the QUBO
formulation. In contrast to the previous analysis, we additionally investigate
Simulated Annealing (SA) as a sub-problem QUBO solver in this instance. Due to
the larger problem sizes of the previous sub-problems, the SA routine could
not return results within the runtime boundaries we had set. However, with the
smaller problem sizes, this was not a problem. Fig. 11 displays the results of
these experiments. The previously discussed solvers (Gurobi, Decomp-Gurobi,
Decomp-Leap) exhibited similar patterns to the investigation done for the
larger problem-sizes (Fig. 7). Therefore, we only focus on the QPU and SA-
based decomposition routines and Leap.

The two Decomp solvers exhibit similar good performance judging by the CO2
reduction. Interestingly, the per-customer metrics, and therefore also the
dominating factor in the energy, measure at a very constant level between the
problem sizes. Curiously, SA as the sub-solver performs better in regard to
per-customer metrics than the QPU does, leading to a gap in the energy. Leap
exhibits similar performance as in our previous experiments. Most notably,
although, is that Decomp-QPU seems to perform slightly better than Leap
regarding the optimization objective (much better in regard to the CO2
reduction). That means that our developed hybrid quantum routine does seem to
outperform the general-purpose Leap.

##  VI Conclusion

Our investigation into the feasibility of current quantum computing techniques
for DSR has begun by developing a mathematical formulation utilizing discount
scheduling to shift grid load to more appropriate times. Our formulation
involves providing discretized discounts to multiple customers at different
times to incentivize a change in consumption while ensuring all customers
receive the same amount of electricity. Our central goal is to reduce overall
CO2 emissions while maintaining grid stability and customer well-being,
leading to a constraint quadratic integer optimization problem.

Upon close inspection of the problem, we developed a custom decomposition
algorithm that compartmentalizes the problem into customer chunks. These sub-
problems involve unconstrained integer optimization and can be effectively
addressed on quantum computers if encoded correctly. Moreover, since the
problems are solved sequentially, we incorporated the accumulated errors into
the subsequent optimization problems. Lastly, we developed a post-processing
algorithm that further refines the solution.

In the end, we benchmarked the performance of a classical general-purpose
solver against D-Wave’s Leap hybrid quantum-classical solver and our
customized decomposition method with various (quantum or classical) sub-
solvers employed. We found that, with a linearly increasing timeout for the
problem size, the classical solver fails to produce acceptable results after a
certain problem size, while the quantum-enhanced Leap continues to provide
adequate results. This indicates a potential advantage of solving this
particular problem using Leap over the purely classical counterpart, Gurobi.
Nonetheless, the decomposition method with the classical solver as sub-solver
developed the best-achieved results over the range of problem sizes we
investigated. Furthermore, using quantum or simulated annealing for the QUBO
problems has resulted in good performance. We found that decomposition paired
with quantum annealing returned slightly better energies than Leap.

We remark that the pairing of the decomposition method with Leap with large
chunk sizes might be a promising pathway for utilizing the quantum-enhanced
method for huge instances of this problem. This statement requires further
experiments to be conducted, but we reason that large sub-problems may have
issues to be solved within the time constraints using Gurobi whereas Leap
could find results. Furthermore, the response of the solvers to different
problem parameter settings is a topic for future work. Lastly, making the grid
constraints more realistically sound instead of our reasonably but arbitrarily
chosen constant band is a topic for future research.

## Acknowledgements

The authors acknowledge funding from the German Fed- eral Ministry of
Education and Research under the funding program ”Förderprogramm
Quantentechnologien – von den Grundlagen zum Markt” (funding program quantum
technologies—from basic research to market), project Q-Grid, 13N16177.

## References

  * [1] Jacopo Torriti.  Peak energy demand and demand side response.  Routledge, 2015.
  * [2] Pierluigi Siano.  Demand response and smart grids—a survey.  Renewable and sustainable energy reviews, 30:461–478, 2014.
  * [3] Haider Tarish Haider, Ong Hang See, and Wilfried Elmenreich.  A review of residential demand response of smart grid.  Renewable and Sustainable Energy Reviews, 59:166–178, 2016.
  * [4] Pierre Pinson, Henrik Madsen, et al.  Benefits and challenges of electrical demand response: A critical review.  Renewable and Sustainable Energy Reviews, 39:686–699, 2014.
  * [5] John S Vardakas, Nizar Zorba, and Christos V Verikoukis.  A survey on demand response programs in smart grids: Pricing methods and optimization algorithms.  IEEE Communications Surveys & Tutorials, 17(1):152–178, 2014.
  * [6] Ari Mizel, Daniel A. Lidar, and Morgan Mitchell.  Simple proof of equivalence between adiabatic quantum computation and the circuit model.  Phys. Rev. Lett., 99:070502, Aug 2007.
  * [7] Arnab Das and Bikas K. Chakrabarti.  Colloquium: Quantum annealing and analog quantum computation.  Rev. Mod. Phys., 80:1061–1081, Sep 2008.
  * [8] Tameem Albash and Daniel A. Lidar.  Adiabatic quantum computation.  Rev. Mod. Phys., 90:015002, Jan 2018.
  * [9] Tadashi Kadowaki and Hidetoshi Nishimori.  Quantum annealing in the transverse ising model.  Phys. Rev. E, 58:5355–5363, Nov 1998.
  * [10] A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll.  Quantum annealing: A new method for minimizing multidimensional functions.  Chemical Physics Letters, 219(5):343–348, March 1994.
  * [11] Catherine McGeoch, Pau Farre, and William Bernoudy.  D-Wave Hybrid Solver Service + Advantage: Technology Update.
  * [12] R. Eskandarpour, Kumar Ghosh, A. Khodaei, A. Paaso, and L. Zhang.  Quantum-enhanced grid of the future: A primer.  IEEE Access, pages 188993 – 189002, 2020.
  * [13] Y. Zhou, Z. Tang, N. Nikmehr, P. Babahajiyani, F. Feng, T. C. Wei, H. Zheng and P. Zhang.  Quantum computing in power systems.  iEnergy, 1, 2022.
  * [14] D.J. Tylavsky and G.T. Heydt.  Quantum computing in power system simulation.  2003 IEEE Power Engineering Society General Meeting, 2003.
  * [15] R. Eskandarpour, A. Khodaei, L. Zhang, E. A. Paaso, and S. Bahramirad.  Quantum computing applications in power systems.  CIGRE US National Committee, 2019 Grid of the Future Symposium, 2019\.
  * [16] Rozhin Eskandarpour, Kumar Ghosh, Amin Khodaei, and Aleksi Paaso.  Experimental Quantum Computing to Solve Network DC Power Flow Problem, 2021.
  * [17] Rozhin Eskandarpour, Kumar Ghosh, Amin Khodaei, Liuxi Zhang, Aleksi Paaso, and Shay Bahramirad.  Quantum Computing Solution of DC Power Flow, 2020.
  * [18] Stephen DiAdamo, Corey O’Meara, Giorgio Cortiana, and Juan Bernabé-Moreno.  Practical quantum k-means clustering: Performance analysis and applications in energy grid classification.  IEEE Transactions on Quantum Engineering, 3:1–16, 2022.
  * [19] Vasil S Denchev, Sergio Boixo, Sergei V Isakov, Nan Ding, Ryan Babbush, Vadim Smelyanskiy, John Martinis, and Hartmut Neven.  What is the computational value of finite-range tunneling?  Physical Review X, 6(3):031015, 2016.
  * [20] Andrew Lucas.  Ising formulations of many NP problems.  Frontiers in Physics, 2, 2014.
  * [21] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A Quantum Approximate Optimization Algorithm, November 2014.  arXiv:1411.4028 [quant-ph].
  * [22] N. Nikmehr, P. Zhang, and M. A. Bragin.  Quantum distributed unit commitment: An application in microgrids.  IEEE Transactions on Power Systems, 37:3592 – 3603, 2022.
  * [23] N. Nikmehr, P. Zhang, and M. A. Bragin.  Quantum-enabled distributed unit commitment.  2022 IEEE Power and Energy Society General Meeting (PESGM), 2022\.
  * [24] Samantha Koretsky, Pranav Gokhale, Jonathan M. Baker, Joshua Viszlai, Honghao Zheng, Niroj Gurung, Ryan Burg, Esa Aleksi Paaso, Amin Khodaei, Rozhin Eskandarpour, and Frederic T. Chong.  Adapting Quantum Approximation Optimization Algorithm (QAOA) for Unit Commitment.  In 2021 IEEE International Conference on Quantum Computing and Engineering (QCE), pages 181–187, October 2021.
  * [25] F. Feng, P. Zhang, M. A. Bragin, and Y. Zhou.  Novel resolution of unit commitment problems through quantum surrogate lagrangian relaxation.  IEEE Transactions on Power Systems (Early Access), 2022.
  * [26] A. Ajagekar and F. You.  Quantum computing for energy systems optimization: Challenges and opportunities.  Energy, 179:76–89, 2019.
  * [27] Claudio Gambella and Andrea Simonetto.  Multiblock ADMM heuristics for mixed-binary optimization on classical and quantum computers.  IEEE Transactions on Quantum Engineering, 1:1–22, 2020.
  * [28] Zhongqi Zhao, Lei Fan, and Zhu Han.  Hybrid Quantum Benders’ Decomposition For Mixed-integer Linear Programming, December 2021.  arXiv:2112.07109 [math].
  * [29] Corey O’Meara, Marina Fernández-Campoamor, Giorgio Cortiana, and Juan Bernabé-Moreno.  Quantum software architecture blueprints for the cloud: Overview and application to peer-2-peer energy trading.  In 2023 IEEE Conference on Technologies for Sustainability (SusTech), pages 191–198, 2023.
  * [30] Marina Fernández-Campoamor, Corey O’Meara, Giorgio Cortiana, Vedran Peric, and Juan Bernabé-Moreno.  Community Detection in Electrical Grids Using Quantum Annealing, December 2021.  arXiv:2112.08300 [quant-ph].
  * [31] Supreeth Mysore Venkatesh, Antonio Macaluso, and Matthias Klusch.  BILP-Q: quantum coalition structure generation.  In Proceedings of the 19th ACM International Conference on Computing Frontiers, CF ’22, pages 189–192, New York, NY, USA, May 2022\. Association for Computing Machinery.
  * [32] Supreeth Mysore Venkatesh, Antonio Macaluso, and Matthias Klusch.  GCS-Q: Quantum Graph Coalition Structure Generation.  In Computational Science – ICCS 2023: 23rd International Conference, Prague, Czech Republic, July 3–5, 2023, Proceedings, Part V, pages 138–152, Berlin, Heidelberg, July 2023. Springer-Verlag.
  * [33] Guillaume Chapuis, Hristo Djidjev, Georg Hahn, and Guillaume Rizk.  Finding Maximum Cliques on a Quantum Annealer.  In Proceedings of the Computing Frontiers Conference, CF’17, pages 63–70, New York, NY, USA, May 2017. Association for Computing Machinery.
  * [34] Georg Hahn and Hristo N. Djidjev.  Reducing Binary Quadratic Forms for More Scalable Quantum Annealing.  Technical Report arXiv:1801.08652, arXiv, January 2018.  arXiv:1801.08652 [quant-ph] type: article.
  * [35] Elijah Pelofske, Georg Hahn, and Hristo Djidjev.  Solving Large Maximum Clique Problems on a Quantum Annealer.  In Sebastian Feld and Claudia Linnhoff-Popien, editors, Quantum Technology and Optimization Problems, Lecture Notes in Computer Science, pages 123–135, Cham, 2019. Springer International Publishing.
  * [36] Elijah Pelofske, Georg Hahn, and Hristo Djidjev.  Decomposition Algorithms for Solving NP-hard Problems on a Quantum Annealer.  Journal of Signal Processing Systems, 93(4):405–420, April 2021\.
  * [37] Elijah Pelofske, Georg Hahn, and Hristo Djidjev.  Solving large minimum vertex cover problems on a quantum annealer.  In Proceedings of the 16th ACM International Conference on Computing Frontiers, CF ’19, pages 76–84, New York, NY, USA, April 2019\. Association for Computing Machinery.
  * [38] Gian Giacomo Guerreschi.  Solving Quadratic Unconstrained Binary Optimization with divide-and-conquer and quantum algorithms, January 2021.  arXiv:2101.07813 [quant-ph].
  * [39] Zeqiao Zhou, Yuxuan Du, Xinmei Tian, and Dacheng Tao.  QAOA-in-QAOA: Solving Large-Scale MaxCut Problems on Small Quantum Machines.  Physical Review Applied, 19(2):024027, February 2023.  Publisher: American Physical Society.
  * [40] Teague Tomesh, Zain H. Saleem, and Martin Suchara.  Quantum Local Search with the Quantum Alternating Operator Ansatz.  Quantum, 6:781, August 2022.  arXiv:2107.04109 [quant-ph].
  * [41] Zain H. Saleem, Teague Tomesh, Michael A. Perlin, Pranav Gokhale, and Martin Suchara.  Divide and Conquer for Combinatorial Optimization and Distributed Quantum Computation, July 2022.  arXiv:2107.07532 [quant-ph].
  * [42] Keisuke Fujii, Kaoru Mizuta, Hiroshi Ueda, Kosuke Mitarai, Wataru Mizukami, and Yuya O. Nakagawa.  Deep Variational Quantum Eigensolver: A Divide-And-Conquer Method for Solving a Larger Problem with Smaller Size Quantum Computers.  PRX Quantum, 3(1):010346, March 2022.  Publisher: American Physical Society.
  * [43] Michael Booth, Steven P Reinhardt, and Aidan Roy.  Partitioning Optimization Problems for Hybrid Classical/Quantum Execution.  D-Wave Technical Report Series, October 2017.
  * [44] M. Born and V. Fock.  Beweis des Adiabatensatzes.  Zeitschrift für Physik, 51(3):165–180, March 1928.
  * [45] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose.  Quantum annealing with manufactured spins.  Nature, 473(7346):194–198, May 2011.  Number: 7346 Publisher: Nature Publishing Group.
  * [46] J. A. Nelder and R. Mead.  A Simplex Method for Function Minimization.  The Computer Journal, 7(4):308–313, January 1965.
  * [47] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2023.
  * [48] Catherine McGeoch and Pau Farré.  The Advantage System: Performance Update.

##  Appendix A Glossary

For a better overview of the used symbols and parameters in the formulation of
the DSP, we provide an overview in Table III.

TABLE III: The symbols and parameters in the DSP formulation Label | Range | Description | Units
---|---|---|---
$c$ | $\set{1,\dots,N_{c}}$ | Customer index |
$t$ | $\set{1,\dots,N_{t}}$ | Time step index |
$k$ | $\set{1,\dots,N_{k}}$ | Discount type index |
$j$ | $\set{1,\dots,M}$ | Chunk index |
$d_{c,t}$ | $\mathbb{R}^{+}$ | Customer power draw | kWh
$\chi_{c}$ | $[0,1]$ | Customer susceptibility |
$I_{t}$ | $\mathbb{R}^{+}$ | Grid CO2 intensity | $\frac{\mathrm{g}}{\mathrm{kWh}}$
$\Delta I_{t}$ | $\mathbb{R}$ | Average adjusted CO2 intensity | $\frac{\mathrm{g}}{\mathrm{kWh}}$
$x_{c,t,k}$ | $\\{0,1\\}$ | Binary decision variables |
$z_{c,t}$ | $\\{-z_{\mathrm{m}},\dots,z_{\mathrm{m}}\\}$ | Discounts |
$D_{c}$ | $\mathbb{R}^{+}$ | Total consumption per customer | kWh
$D_{t}$ | $\mathbb{R}^{+}$ | Total consumption per timestep | kWh
$\widetilde{D}_{t}$ | $\mathbb{R}^{+}$ | Mutable consumption | kWh
$\zeta_{t}$ | $[-z_{\mathrm{m}},z_{\mathrm{m}}]$ | Effective discount |
$\xi_{t}^{j}$ | $[-z_{\mathrm{m}},z_{\mathrm{m}}]$ | Effective chunk discount |
$N_{k}$ | $\mathbb{N}$ | Number of discount types |
$z_{m}$ | $[0,1]$ | Maximal discount |
$\lambda$ | $\mathbb{R}^{+}$ | Consumption deviation penalty |
$p_{t}^{\mathrm{high}}$ | $\mathbb{R}^{+}$ | Upper momentary consumption | kWh
|  | deviation bound (like lower) |
$\rho$ | $\mathbb{R}^{+}$ | Discount $L2$-Regularization |
$\tau$ | $\mathbb{R}^{+}$ | Discount change penalty |
$m$ | $\mathbb{N}$ | Chunk size |
$M$ | $\mathbb{N}$ | Number of chunks |
$r$ | $\mathbb{N}$ | Post-processing cut-off |
$Q$ | $\mathbb{N}$ | Number of bits in encoding |