# The shape of higher-dimensional state space: Bloch-ball analog for a qutrit

Christopher Eltschka  Institut für Theoretische Physik, Universität
Regensburg, D-93040 Regensburg, Germany  Marcus Huber  Atominstitut,
Technische Universität Wien, 1020 Vienna, Austria  Institute for Quantum
Optics and Quantum Information - IQOQI Vienna, Austrian Academy of Sciences,
Boltzmanngasse 3, 1090 Vienna, Austria  Simon Morelli  Institute for Quantum
Optics and Quantum Information - IQOQI Vienna, Austrian Academy of Sciences,
Boltzmanngasse 3, 1090 Vienna, Austria  Atominstitut, Technische Universität
Wien, 1020 Vienna, Austria  Jens Siewert<EMAIL_ADDRESS>Departamento de
Química Física, Universidad del País Vasco UPV/EHU, E-48080 Bilbao, Spain
IKERBASQUE Basque Foundation for Science, E-48009 Bilbao, Spain

(22 June, 2021)

###### Abstract

Geometric intuition is a crucial tool to obtain deeper insight into many
concepts of physics. A paradigmatic example of its power is the Bloch ball,
the geometrical representation for the state space of the simplest possible
quantum system, a two-level system (or qubit). However, already for a three-
level system (qutrit) the state space has eight dimensions, so that its
complexity exceeds the grasp of our three-dimensional space of experience.
This is unfortunate, given that the geometric object describing the state
space of a qutrit has a much richer structure and is in many ways more
representative for a general quantum system than a qubit. In this work we
demonstrate that, based on the Bloch representation of quantum states, it is
possible to construct a three dimensional model for the qutrit state space
that captures most of the essential geometric features of the latter. Besides
being of indisputable theoretical value, this opens the door to a new type of
representation, thus extending our geometric intuition beyond the simplest
quantum systems.

##  1 Introduction

Nowadays virtually every student of quantum mechanics learns about the Bloch
sphere and the Bloch ball as the geometrical representations of pure and mixed
states of qubits, respectively. These objects have become indispensable for
developing an intuition of elementary concepts such as basic quantum
operations and the action of decoherence [1], or more advanced topics like the
Majorana representation of symmetric multi-qubit states [2]. Nonetheless, the
more experienced practitioner in the field of quantum mechanics is aware of
various shortcomings of the Bloch ball if systems of higher dimension $d>2$
are to be discussed. For example, the entire surface of the Bloch ball is
covered by pure states, whereas most of the boundary of higher-dimensional
state spaces is formed by mixed states. Those parts of the boundary may either
be flat or curved, however, their curvature is different from that of the
pure-state surface parts, which actually represents a single unitary orbit.
Another important fact not shown by the Bloch ball is that not every
orthogonal transformation can be applied to any quantum state: Whenever a
state is not an element of the inscribed sphere of the state space of maximum
radius there exist rotations in Bloch space that take it outside the state
space and therefore are not allowed. Consequently, it is a long-standing
question whether it is possible to construct a consistent three-dimensional
model for higher-dimensional state spaces $\mathcal{Q}_{d}$, $d\geq 3$, that
captures at least a part of these important geometric features.

##  2 Bloch ball for qubits

Let us briefly recapitulate the properties of the Bloch ball for qubits,
$d=2$. According to Fano [3, 4] density operators of qubits are parametrized
by using the Pauli matrices $\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$ and the
identity $\mathds{1}_{2}$

$\displaystyle\rho\ $ $\displaystyle=\ \frac{1}{2}\left(s_{0}\mathds{1}_{2}+x\
\sigma_{x}+y\ \sigma_{y}+z\ \sigma_{z}\right)$ (1a) $\displaystyle s\ $
$\displaystyle=\ \operatorname{Tr}\left(\sigma_{s}\rho\right)\ \ ,$ (1b)

with real numbers $s\in\\{x,y,z\\}$, $|s|\leq 1$, and the normalization
$s_{0}=1$. The state space $\mathcal{Q}_{2}$ of all qubit density operators is
represented by a ball of radius 1 about the origin of $\mathbb{R}^{3}$, that
is, each point $(x,y,z)$ of this ball corresponds to exactly one state $\rho$.
Pure states lie on the surface (forming a connected set), whereas mixed states
are located inside the Bloch sphere, with the fully mixed state at the center.
Interestingly, the convex combinations of two states,
$\rho=\lambda\rho_{1}+(1-\lambda)\rho_{2}$ ($0\leq\lambda\leq 1$), are given
by the straight line connecting the points of $\rho_{1}$ and $\rho_{2}$. A
common choice of computational basis states is
$\\{\left|{0}\right\rangle,\left|{1}\right\rangle\\}$, which correspond to
north and south pole of the sphere and form the regular simplex $\Delta_{1}$,
the special case of $\Delta_{d-1}$ for $d=2$. This shows that the vectors in
$\mathbb{R}^{3}$ belonging to the basis states are not orthogonal. These
vectors are orthogonal in $\mathbb{R}^{4}$ including the direction of $s_{0}$,
but their projections into $\mathbb{R}^{3}$ lose this property.

Every density matrix can be obtained by unitarily rotating a diagonal density
matrix, that resembles a classical probability distribution. Since for $d=2$
the special unitary group SU(2) is the universal cover group of the rotation
group SO(3), every rotation of the Bloch ball $Q_{2}\subset\mathbb{R}^{3}$ has
a corresponding unitary rotation in the state space. From this perspective,
the Bloch ball is obtained by all possible rotations of the simplex
$\Delta_{1}$ in $\mathbb{R}^{3}$.

To the best of our knowledge, the idea that the parametrization in Eq. (1)
entails a useful visualization for the non-unitary dynamics of a spin
$\frac{1}{2}$, e.g., in a situation of radiation damping was put forward by
Feynman and co-workers [5].

##  3 A Bloch-ball analog for a qutrit

Over the years much work has been done to elucidate the geometry of the state
space of higher level systems, with special focus on the qutrit state space
$\mathcal{Q}_{3}$ [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25]. To develop an intuition for the full high-dimensional
geometry of the qutrit state space, subsets, cross sections and projections
onto two and three dimensions were extensively studied [6, 7, 10, 12, 15, 18,
17, 20, 21, 22, 23] and multi-parameter representations of qutrit states were
developed [14, 21, 26, 24, 25]. While these approaches can reproduce many
geometric properties correctly, they do not give a global view of the Bloch
body.

In this section we review known facts about the higher dimensional state
space, with focus on dimension $d=3$, collecting a list of requirements we
wish our model to reproduce. We then construct a three dimensional global
model of the state space, that reproduces astonishingly many properties of
$\mathcal{Q}_{3}$.

###  3.1 Qutrit geometry

In analogy to qubits, density matrices describing the state of a $d$-level
quantum system can be parametrized by the identity $\mathds{1}_{d}$ and
$d^{2}-1$ traceless Hermitian matrices [4, 6, 27]. So the quantum state space
$\mathcal{Q}_{d}$ can be represented by a subset in $\mathbb{R}^{d^{2}-1}$,
constrained by inequalities that arise from the positivity of the density
operators [10, 11]. While it is still true that the quantum state space
$\mathcal{Q}_{d}$ is obtained by rotating the simplex $\Delta_{d-1}$ in
$\mathbb{R}^{d^{2}-1}$, it is no longer a sphere as not all rotations are
allowed. In fact, for $d>2$ the special unitary group SU($d$) is a proper
subgroup of the rotation group SO($d^{2}$-1), meaning that the quantum state
space $\mathcal{Q}_{d}$ is a proper subset of the $d^{2}-1$ dimensional
(Hilbert-Schmidt) ball.

For qutrits the density matrices are parametrized by the identity
$\mathds{1}_{3}$ and the normalized Gell-Mann matrices
$\operatorname{\mathsf{X}}_{j}$, $\operatorname{\mathsf{Y}}_{k}$
($j,k=1,2,3$), and $\operatorname{\mathsf{Z}}_{1}$,
$\operatorname{\mathsf{Z}}_{2}$ (see Appendix),

$\displaystyle\rho\ =$ $\displaystyle\
\frac{1}{3}\left(\mathds{1}_{3}+\sum_{j=1}^{3}x_{j}\operatorname{\mathsf{X}}_{j}+\sum_{k=1}^{3}y_{k}\operatorname{\mathsf{Y}}_{k}+\sum_{l=1}^{2}z_{l}\operatorname{\mathsf{Z}}_{l}\right)$
(2a) $\displaystyle x_{j}=$
$\displaystyle\operatorname{Tr}\left(\operatorname{\mathsf{X}}_{j}\rho\right)\,,\
y_{k}=\operatorname{Tr}\left(\operatorname{\mathsf{Y}}_{k}\rho\right)\,,\
z_{l}=\operatorname{Tr}\left(\operatorname{\mathsf{Z}}_{l}\rho\right)\,,$ (2b)

with real numbers $x_{j}$, $y_{k}$, and $z_{l}$, hence their space has eight
dimensions. A Euclidean metric, corresponding to that of our everyday
geometric experience, is induced in this space by the Hilbert-Schmidt norm,
$\|\rho_{1}-\rho_{2}\|_{2}$$\equiv\sqrt{3\operatorname{Tr}\left[(\rho_{1}-\rho_{2})^{2}\right]}$.

The geometric properties of $\mathcal{Q}_{3}$ we wish for a global model to
reproduce were thoroughly discussed by Bengtsson et al. in Ref. [18].
i) Most importantly, $\mathcal{Q}_{3}$ is a convex set with the topology of a
ball, so the model should share this characteristics. There are no pieces of
lower dimension attached to it (“no hair” condition).
ii) The actual Bloch body is neither a polytope nor a smooth object.
iii) The Bloch body has an outer sphere of radius $\sqrt{2}$ and an inner
sphere of radius $1/\sqrt{2}$ [28].
iv) The pure states (rank 1) form a connected set on the surface at maximum
distance $\sqrt{2}$ from the completely mixed state
$\frac{1}{3}\mathds{1}_{3}$. Its measure is zero compared to that of
$\mathcal{Q}_{3}$. In particular we aim at prominently displaying the three
pure states corresponding to the preferred basis for the model.
v) Density matrices on the surface of $\mathcal{Q}_{3}$ are of rank 1 or 2,
whereas states inside $\mathcal{Q}_{3}$ are of full rank (rank 3).
Finally, there are some additional properties specifically related to the
nature of $\mathcal{Q}_{3}$ as a convex set:
vi) The set of quantum states is self-dual.
vii) All cross sections of $\mathcal{Q}_{3}$ do not have non-exposed faces.
All corners of two-dimensional projections of $\mathcal{Q}_{3}$ are
polyhedral.

###  3.2 A three dimensional model for a qutrit

Surprisingly, it is indeed possible to find an object in $\mathbb{R}^{3}$
representing $\mathcal{Q}_{3}$ that obeys most of the requirements in this
list. In fact, there are (at least) two solutions with slightly different
advantages. First we construct an object that fulfills properties i–iv) and
also partially v). It represents a valid model for $\mathcal{Q}_{3}$ that we
call $Q_{3}^{(1)}$. Then we show that this first solution can be extended to a
model $Q_{3}^{(2)}$ that possesses in particular also the property vi). Yet
that model does not obey vii).

It is well known that the computational basis states in $d$ dimensions form
the corners of a regular simplex $\Delta_{d-1}$ in $\mathcal{Q}_{d}$ [6, 18],
that is, for a qutrit a basis will be represented by $\Delta_{2}$, an
equilateral triangle. We insist that our model for $\mathcal{Q}_{3}$ displays
and emphasizes one particular basis
$\\{\left|{0}\right\rangle,\left|{1}\right\rangle,\left|{2}\right\rangle\\}$,
because the mapping between physical states and points in the model will
depend on this choice. Hence, we use the coordinates [cf. Eq. (2b)]

$\displaystyle z_{1}\ =\
\operatorname{Tr}\left(\operatorname{\mathsf{Z}}_{1}\rho\right)\ \ ,\ \ \ \
z_{2}\ =\ \operatorname{Tr}\left(\operatorname{\mathsf{Z}}_{2}\rho\right)\ \
,$

to faithfully represent the diagonal matrices as a simplex $\Delta_{2}$ in the
horizontal coordinate plane. In particular, the basis state
$\left|{2}\right\rangle$ corresponds to $(z_{1},z_{2})=(0,-\sqrt{2})$, whereas
the states $\left|{0}\right\rangle$ and $\left|{1}\right\rangle$ are located
in $(\pm\sqrt{\frac{3}{2}},\sqrt{\frac{1}{2}})$, respectively. The completely
mixed state lies at the origin $(0,0)$.

For the remaining six coordinates we are left with only one direction in
$\mathbb{R}^{3}$. Here, we propose to use the coordinate

$\displaystyle w\ =$ $\displaystyle\ \sqrt{\sum_{j=1}^{3}x_{j}^{2}+y_{j}^{2}}\
\ ,$ (3)

which assumes only non-negative values. Hence we have our first model for
$\mathcal{Q}_{3}$,

$\displaystyle Q_{3}^{(1)}=\\{(z_{1},z_{2},w)$
$\displaystyle\in\mathbb{R}^{3}\ \mathrm{s.t.}\ $
$\displaystyle\mathrm{Eqs.}~{}\eqref{eq:bloch3b},\eqref{eq:vert}\
\forall\rho\in\mathcal{Q}_{3}\\}.$ (4)

The interpretation for the coordinates of a point $P=(z_{1},z_{2},w)$ is
simple. Imagine the state $\rho$ corresponding to $P$ written as a sum of its
diagonal and offdiagonal parts,

$\rho\ =\ \rho_{\mathrm{diag}}\ +\ \rho_{\mathrm{offdiag}}\ \ .$

While the distance of $P$ from the origin in the plane equals the Hilbert-
Schmidt length $\sqrt{3\operatorname{Tr}(\rho_{\mathrm{diag}}^{2})-1}$ of the
diagonal part of the Bloch vector for $\rho$, the vertical distance of $P$
from the plane is the Hilbert-Schmidt length of this Bloch vector’s
offdiagonal part, $\sqrt{3\operatorname{Tr}(\rho_{\mathrm{offdiag}}^{2})}$.
The result is shown in Fig. 1.

Figure 1: The model $Q_{3}^{(1)}$ of the qutrit Bloch body $\mathcal{Q}_{3}$
according to Eq. (4). The location of the basis states
$\left|{j}\right\rangle$ is specified. The semicircular surface in the
foreground is the image of the Bloch ball spanned by the states
$\left|{0}\right\rangle$ and $\left|{1}\right\rangle$ and, in this
representation, has the “north pole”
$\left|{+_{2}}\right\rangle=\frac{1}{\sqrt{2}}(\left|{0}\right\rangle+\left|{1}\right\rangle)$.
The point with the largest $w$ coordinate is the image of the maximally
coherent superposition
$\left|{+_{3}}\right\rangle=\frac{1}{\sqrt{3}}(\left|{0}\right\rangle+\left|{1}\right\rangle+\left|{2}\right\rangle)$.
Note that all the pure states with $|\langle j|\psi\rangle|^{2}=\frac{1}{3}$
get mapped to this point, in particular all the bases which are mutually
unbiased with
$\\{\left|{0}\right\rangle,\left|{1}\right\rangle,\left|{2}\right\rangle\\}$.

###  3.3 Global geometric properties

Let us analyze the properties of this first model for $\mathcal{Q}_{3}$ with
respect to our list of requirements. The surface of this object corresponds to
a hemisphere of radius $\sqrt{2}$, where the three spherical segments beyond
the triangle in the plane are cut off. Evidently, it is both convex and simply
connected, without anything else attached to it. The flat surfaces of the cuts
are connected with the smooth upper boundary by a sharp corner, so the object
is neither smooth nor a polytope. Since our model preserves the length of the
Bloch vector, it is still circumscribed by an outer sphere of radius
$\sqrt{2}$ and center at the origin. The radius of the inner (hemi-)sphere
coincides with that of the in-circle of the simplex $\Delta_{2}$ in the
horizontal plane and equals $1/\sqrt{2}$.

The spherical part of the surface above the ground plane corresponds to the
set of pure states, hence to density matrices of rank 1. They form a simply
connected surface of measure zero at maximal distance $\sqrt{2}$ from the
origin. The cuts are half circles, in fact, these flat surfaces are the images
of the three two-dimensional Bloch balls corresponding to the pairs of states
$\\{\left|{0}\right\rangle,\left|{1}\right\rangle\\}$,
$\\{\left|{0}\right\rangle,\left|{2}\right\rangle\\}$, and
$\\{\left|{1}\right\rangle,\left|{2}\right\rangle\\}$ (cf. also Fig. 1). This
means, these surfaces correspond to states of at most rank 2. It is understood
that the base of $Q_{3}^{(1)}$ represents an artificial cut – similar to the
base of the sculpture of a bust – that does not represent a boundary of
$\mathcal{Q}_{3}$. The states on the other surfaces are of lower rank, all
states of rank 3 reside in the interior of the model. However rank-2 states
that have all of the basis vectors $\left|{0}\right\rangle$,
$\left|{1}\right\rangle$, $\left|{2}\right\rangle$ in their span are located
inside $Q_{3}^{(1)}$, although in $\mathcal{Q}_{3}$ they are part of the
boundary. But the rank-2 states do not cover the entire interior of
$Q_{3}^{(1)}$. Evidently, all points inside the inner sphere are of rank 3,
since rank-2 states have purity of at least $\frac{1}{2}$. However the set of
points corresponding only to rank-3 states is larger than that; see Sec. 3.4
for details. Conversely, the interior points outside that set correspond both
to rank-2 and rank-3 states. In particular, in agreement with property v), the
images of rank-3 states cover the complete interior. Thus, we have established
that our construction $Q_{3}^{(1)}$ has all the desired properties i–iv) of
our list, property v) is at least partially satisfied.

We may ask how the coordinates in $Q_{3}^{(1)}$ are related to those of the
actual state space $\mathcal{Q}_{3}$. Clearly our model is neither a
projection nor a cross section. Instead, it resembles a representation in
cylindrical coordinates: Our “diagonal coordinates” $z_{1}$, $z_{2}$
correspond to two longitudinal axes, whereas each
$r_{j}=\sqrt{x_{j}^{2}+y_{j}^{2}}$ (where $j=1,2,3$) may be viewed as a radial
coordinate belonging to mutually orthogonal directions. Our model does not
display the polar angles and shows only the total radial distance.

###  3.4 Local algebraic properties

We have mentioned that the simplex $\Delta_{2}$ represents the diagonal states
faithfully, it is isomorphic to that set just as the Bloch ball is isomorphic
to $\mathcal{Q}_{2}$. But our three-dimensional model $Q_{3}^{(1)}$ cannot
represent all states of the eight-dimensional state space faithfully. In fact,
a general point of $Q_{3}^{(1)}$ corresponds to infinitely many states and the
simplex representing the basis states
$\left|{j}\right\rangle\\!\left\langle{j}\right|$, $j=0,1,2$ is the only set
of states with a unique preimage. States mapped onto the same point in the
model belong to the equivalence class of states with the same diagonal entries
and purity, they form a five-dimensional manifold in the original state space
$\mathcal{Q}_{3}$. These subspaces are closed under the action of unitary
operators that commute with $\operatorname{\mathsf{Z}}_{1}$ and
$\operatorname{\mathsf{Z}}_{2}$, in particular diagonal unitary operators. As
a consequence the model is invariant under these transformations.

The action of general unitary transformations, however, is displayed by our
model. Since the purity remains unchanged, all points in the unitary orbit
have the same distance from the origin. This orbit forms a subset of a sphere,
its shape depends on the eigenvalues of the transformed state. As previously
elucidated, a point in the interior represents a whole subspace of states with
possibly different eigenvalues, so it makes little sense to talk about the
orbit of a point in our model. However, diagonal states in the ground triangle
are depicted faithfully and therefore uniquely identify a unitary orbit. One
only needs to use the eigenbasis of the state to investigate which points can
be reached by unitary transformations. The permutations of the diagonal
entries correspond to SU(3) rotations, therefore the unitary orbit includes
six points on the triangle (for non-degenerate eigenvalues). Since the
eigenvalues of a matrix majorize the diagonal entries, the orbit remains on a
sphere above this hexagon. Vice versa, every point on the sphere above the
hexagon can be reached through a unitary transformation, see Fig. 2. If the
eigenvalues are degenerate, the hexagon becomes a triangle. A special case are
the pure states, the unitary orbit of a pure state coincides with the upper
surface of $Q_{3}^{(1)}$.

This visualization establishes a direct connection to the action of doubly
stochastic matrices on classical probability distributions. A classical
probability distribution $\mathsf{p}$ majorizes exactly the probability
distributions $M.\mathsf{p}$, where $M$ is a doubly stochastic matrix. The set
of distributions majorized by $\mathsf{p}$ is called the Birkhoff polytope
[6]. The action of doubly stochastic matrices on classical probability
distributions corresponds to the action of unitary transformations on
normalized diagonal matrices. That is, for every doubly stochastic matrix $M$
there exists a unitary matrix $U$, such that
$M.\mathsf{p}=\operatorname{diag}(U.\mathsf{D}.U^{\dagger})$, where
$\mathsf{D}$ is a diagonal matrix with
$\mathsf{p}=\operatorname{diag}(\mathsf{D})$. In the model this is visualized
by unitarily transforming a diagonal state and then projecting it onto the
ground triangle. This exactly reproduces the Birkhoff polytope, but it also
generalizes it in the sense that the norm of the offdiagonal part remains
visible.

Figure 2: The SU(3) orbit of the state
$\rho_{0}=\frac{6}{10}\left|{0}\right\rangle\\!\left\langle{0}\right|+\frac{3}{10}\left|{1}\right\rangle\\!\left\langle{1}\right|+\frac{1}{10}\left|{2}\right\rangle\\!\left\langle{2}\right|$
is represented by the blue spherical surface. The image of $\rho_{0}$ together
with the other five points in the $z_{1}$-$z_{2}$ plane forms the Birkhoff
polytope (see text). The eigenvalue vectors of the diagonal states inside the
polytope are majorized by that of $\rho_{0}$,
$(\frac{6}{10},\frac{3}{10},\frac{1}{10})$.

The images of SU(3) orbits also allow us to identify the set of images of
rank-2 states. The rank-2 diagonal states are exactly the sides of the base
triangle. Therefore the rank-2 states are given by all the orbits of those
points. It is readily seen that the boundary arcs connecting two triangle
sides form half-circular cones that each have one of the basis states as apex,
and are tangential to the insphere at their base. The points inside those
cones, as well as the points inside the inner sphere, are not covered by those
orbits, and therefore correspond only to states of rank 3. This rank-3 only
region is convex; indeed it is the convex hull of the interior of the insphere
and the interior of the base triangle.

Note that the mirror image of this region and its boundary forms the lower
part of our second model $Q_{3}^{(2)}$ which will be described in Sec. 4.

The map that takes the qutrit state space to our model is decidedly nonlinear,
but amazingly retains some linear properties of the state space for the image.
Any convex combination of two diagonal qutrit states
$\delta=\lambda\delta_{1}+(1-\lambda)\delta_{2}$ with $\lambda\in(0,1)$ is
again a diagonal state and its image lies on the straight line connecting the
image points of $\delta_{1}$ and $\delta_{2}$. But even the convex
combinations $\rho=\lambda\sigma+(1-\lambda)\delta$ of an arbitrary qutrit
state $\sigma$ with a diagonal state $\delta$ are located on the straight line
connecting $\sigma$ and $\delta$. The mixture of the diagonal parts is
faithfully represented, but in this case, this holds also for the offdiagonal
parts, because $\delta_{\mathrm{offdiag}}$ simply vanishes. Consequently, for
any problem involving two arbitrary states we can find a visualization of
their mixture by means of a straight line: It is enough to consider the
problem in the eigenbasis of one of the states.

Another noteworthy aspect of our model $Q_{3}^{(1)}$ is the separate treatment
of diagonal and offdiagonal parts of the state. The diagonal matrices may be
viewed as states in a classical probability space [18]. They become
nonclassical by adding coherences, i.e., an offdiagonal part. The $w$
coordinate in our model is the 2-norm for the offdiagonal part of the state.
As the 1-norm of the offdiagonal part is established as a coherence measure
[29] and the 1-norm is lower-bounded by the 2-norm, our model explicitly shows
a lower bound to the coherence of a given quantum state. Correspondingly, the
maximally coherent state
$\left|{+_{3}}\right\rangle=\frac{1}{\sqrt{3}}(\left|{0}\right\rangle+\left|{1}\right\rangle+\left|{2}\right\rangle)$
is located at the north pole of the hemisphere.

##  4 Self-dual qutrit geometry

The model we have discussed so far does not meet the requirements v) and vi)
of our list. In particular it is not self-dual, that is, our three-dimensional
Bloch vectors do not fulfill the following condition. Let $\vec{\xi}$ be the
Bloch vector belonging to a state on the boundary of the state space. Then for
the set of Bloch vectors $\vec{\eta}$ defining the dual hyperplanes enveloping
the state space on the side opposite to $\vec{\xi}$ we have [18]

$\displaystyle\vec{\xi}\cdot\vec{\eta}\ =\ -1\ \ .$ (5)

Our model so far is not self-dual simply because for vectors $\vec{\xi}$ on
the surface with $w>0$ the dual planes do not touch the lower boundary of our
Bloch body model. We will now demonstrate that the model can be extended by
adding a “lower part” with coordinates $w<0$ so that the entire object becomes
self-dual.

###  4.1 Self-dual extension of $Q_{3}^{(1)}$

To achieve self-duality, consider first the pure states $\Pi$ with
$\operatorname{Tr}(\Pi^{2})=\operatorname{Tr}(\Pi)=1$. If we denote the vector
of Gell-Mann matrices by $\vec{\mathsf{h}}$ and the Bloch vector of the state
$\Pi$ by $\vec{\pi}$, we can write
$\Pi=\frac{1}{3}(\mathds{1}_{3}+\vec{\pi}\cdot\vec{\mathsf{h}})$. Recall that
reversing the sign of all coordinates in the qubit Bloch sphere amounts to a
point reflection operation at the maximally mixed state. This fact suggests
that we might try reversing the sign of $w$ in a process of reversing all of
the coordinate signs, for example in an operation of the kind $\rho\
\longrightarrow\ \alpha\mathds{1}_{3}-\rho$ (where $\alpha>0$). We have to
make sure that the result is again positive. In Ref. [30] the so-called
universal state inversion map $\mathcal{S}(\rho)=\nu(\mathds{1}-\rho)$ was
introduced; it guarantees that the result actually is a state. We choose the
prefactor $\nu=\frac{1}{2}$ so as to normalize the resulting qutrit state and
write

$\displaystyle\mathcal{S}(\rho)\ \equiv\ \tilde{\rho}\ =\
\frac{1}{2}(\mathds{1}_{3}-\rho)\ \ .$ (6)

The inverted state $\tilde{\rho}$ is positive and for pure states we have
$\operatorname{Tr}(\Pi\tilde{\Pi})=0$. Therefore,

$\displaystyle 0\ =\ 1\ +\ \vec{\pi}\cdot\vec{\tilde{\pi}}\ \ ,$ (7)

that is, the Bloch vector $\vec{\tilde{\pi}}$ defines the dual plane for
$\vec{\pi}$ (and vice versa). As the image of $\vec{\pi}$ determines a point
in the spherical part of the boundary of our model, also the image of
$\vec{\tilde{\pi}}$ lies on a spherical surface, however, with half the radius
because of the prefactor $\nu=\frac{1}{2}$. The result is that the spherical
part of our model $Q^{(1)}_{3}$ gets replicated in the region $w<0$, whereat
it is point-reflected at the origin and scaled down by a factor $\frac{1}{2}$.

Figure 3: Constructing the dual of the surface of the model $Q_{3}^{(1)}$.
Red: Cross section for determining the dual states $\tilde{\rho}(p)$ of the
mixtures $\rho(p)$ in Eq. (8). Gray: We have included the part that has been
cut from the full Bloch body (see Fig. 4, which shows the complete object
$Q_{3}^{(2)}$ from the same point of view). By rotating the section
highlighted in red by the angle $\beta$ (see text) one obtains the lower
conical part of $Q^{(2)}_{3}$.

This reasoning cannot be applied for the mixed boundary states of the model.
Rather, we apply state inversion $x_{j}\rightarrow-x_{j}$, $y_{k}\rightarrow-
y_{k}$, $z_{l}\rightarrow-z_{l}$ and determine the boundary states by explicit
calculation. Consider, for example, the mixtures (cf. Fig. 3)

$\displaystyle\rho(p)\ =\
p\left|{+_{2}}\right\rangle\\!\left\langle{+_{2}}\right|\ +\
\frac{1-p}{2}(\left|{0}\right\rangle\\!\left\langle{0}\right|+\left|{1}\right\rangle\\!\left\langle{1}\right|)$
(8)

with $0\leq p\leq 1$. It is straightforward to show that the dual states are
given by

$\displaystyle\tilde{\rho}(p)\ =\ $
$\displaystyle\frac{1}{3}\bigg{(}\mathds{1}_{3}-\frac{\sqrt{2}}{1+3p}\operatorname{\mathsf{Z}}_{2}-\frac{\sqrt{6}\
p}{1+3p}\operatorname{\mathsf{X}}_{1}\bigg{)}$ $\displaystyle=\
\frac{1}{3}\bigg{(}\mathds{1}_{3}$
$\displaystyle-(1-q)\sqrt{2}\operatorname{\mathsf{Z}}_{2}-q\left[\frac{\sqrt{2}}{4}\operatorname{\mathsf{Z}}_{2}+\frac{\sqrt{6}\
p}{4}\operatorname{\mathsf{X}}_{1}\right]\bigg{)}$ (9)

with $q=\frac{4p}{1+3p}$. That is, also $\tilde{\rho}(p)$ lies on a straight
line connecting $\left|{2}\right\rangle\\!\left\langle{2}\right|$ and the dual
of $\left|{+_{2}}\right\rangle\\!\left\langle{+_{2}}\right|$. An analogous
calculation applies to all mixtures of
$\frac{1}{2}(\left|{0}\right\rangle\\!\left\langle{0}\right|+\left|{1}\right\rangle\\!\left\langle{1}\right|)$
with pure states on the circumference of the
$\\{\left|{0}\right\rangle,\left|{1}\right\rangle\\}$ “Bloch sphere”. Instead
of considering the cross section in the $z_{2}$-$w$ plane for the states
$\tilde{\rho}(p)$ we can rotate this plane about the $z_{2}$ axis by an angle
$\beta$ (see Fig. 3), so that the pure state on the circumference becomes
$\cos{\left(\frac{\pi-2\beta}{4}\right)}\left|{1}\right\rangle+\sin{\left(\frac{\pi-2\beta}{4}\right)}\left|{0}\right\rangle$.

The shape of the corresponding parts of the cross section remains the same as
in the one shown in Fig. 3. Consequently, the dual states of the qubit Bloch
ball $\\{\left|{0}\right\rangle,\left|{1}\right\rangle\\}$ are located on a
half circular cone with apex in
$\left|{2}\right\rangle\\!\left\langle{2}\right|$ (cf. Fig. 4). Due to the
three-fold symmetry of $Q_{3}^{(1)}$ the regions near the other pure states
$\left|{0}\right\rangle\\!\left\langle{0}\right|$,
$\left|{1}\right\rangle\\!\left\langle{1}\right|$ can be constructed
analogously.

Once we have found the shape [i.e., the boundary coordinates
$w_{\mathrm{max}}(z_{1},z_{2})$ and $w_{\mathrm{min}}(z_{1},z_{2})$] of the
self-dual object, we also need to provide a rule to decide whether a given
state $\rho$ from inside $\mathcal{Q}_{3}$ is mapped to the “upper” ($w>0$) or
the “lower” ($w<0$) part of the object. In other words, from $\rho$ we can
determine the coordinates $z_{1}$, $z_{2}$, $|w|$ as well as
$w_{\mathrm{min}}(z_{1},z_{2})$, but how can we decide about the sign of $w$?

In the following we describe a procedure to define this mapping $\rho\
\longrightarrow\ P(z_{1},z_{2},w)$. It is both straightforward and consistent
(however, there might be other choices). Let $\rho$ be an arbitrary qutrit
state with smallest eigenvalue $\lambda_{\mathrm{min}}$. First we subtract as
much of the fully mixed state so that one eigenvalue vanishes, i.e.,

$\rho^{\prime}\ =\
\frac{1}{1-3\lambda_{\mathrm{min}}}\left(\rho-\lambda_{\mathrm{min}}\mathds{1}_{3}\right)$
(10)

is a rank-2 state with Bloch coordinates
$(z_{1}^{\prime},z_{2}^{\prime},|w^{\prime}|)$. Then

$\
\left.\begin{array}[]{l}|w^{\prime}|>|w_{\mathrm{min}}(z_{1}^{\prime},z_{2}^{\prime})|\\\\[5.69054pt]
|w^{\prime}|\leq|w_{\mathrm{min}}(z_{1}^{\prime},z_{2}^{\prime})|\end{array}\right\\}\
\Longrightarrow\ \left\\{\begin{array}[]{l}w>0\\\\[5.69054pt]
w<0\end{array}\right.\ \ .$ (11)

This implies that in the case of equality the state $\rho^{\prime}$ is located
at the lower boundary in $Q_{3}^{(2)}$. Moreover, together with
$\rho^{\prime}$ also all its mixtures with $\frac{1}{3}\mathds{1}_{3}$ belong
to the lower part of $Q_{3}^{(2)}$, ensuring its compactness. The meaning of
this definition is clear: In comparison to the structure of $Q_{3}^{(1)}$, in
$Q_{3}^{(2)}$ we redistribute the location of the rank-2 states
$\rho^{\prime}$. If the purity of $\rho^{\prime}$ is large enough for its
diagonal coordinates $(z_{1}^{\prime},z_{2}^{\prime})$ it needs to remain in
the upper part, $w>0$, otherwise it goes to the lower part, $w<0$.

This concludes the construction of the self-dual model $Q^{(2)}_{3}$ for the
qutrit Bloch body, see Fig. 4:

$\displaystyle Q_{3}^{(2)}=\big{\\{}(z_{1},$ $\displaystyle
z_{2},w)\in\mathbb{R}^{3}\ \mathrm{s.t.}\ $
$\displaystyle\mathrm{Eqs.}~{}\eqref{eq:bloch3b},\eqref{eq:vert},\eqref{eq:sign-w}\forall\rho\in\mathcal{Q}_{3}\big{\\}}\
.$ (12)

We note that the extension to a self-dual object can be applied as well to the
half-circle model $Q_{2}^{(1)}$ of the $d=2$ Boch ball. Then, however, the
rule given in Eq. (11) to determine the sign of $w$ (without any addendum)
does not improve the model. This is because in $d=2$ all offdiagonal elements
can be made real and positive by applying a single appropriate diagonal
unitary to the state. In contrast, for $d>2$ at most $2(d-1)$ offdiagonal
matrix elements of a generic state can simultaneously be made positive by
applying a diagonal unitary.

Figure 4: The self-dual model $Q_{3}^{(2)}$ of the qutrit state space
$\mathcal{Q}_{3}$.

###  4.2 More convexity properties and relation to previous work

Convexity properties of the state space in the context of the Bloch-ball
representation for $d=2$ (and even for $d=3$) were discussed early on, e.g.,
by Bloore [7] who analyzed the stratification of the state space with respect
to the matrix rank. For a qutrit, all states of reduced rank $r<3$ are part of
the boundary. The four-dimensional set of rank-1 states
$\left|{\psi}\right\rangle\\!\left\langle{\psi}\right|$ forms the extremal
states of the convex state space, as they do not have a convex decomposition
into other states $\rho_{1}$ and $\rho_{2}$,
$\left|{\psi}\right\rangle\\!\left\langle{\psi}\right|\neq\lambda\rho_{1}+(1-\lambda)\rho_{2}$
with $0\leq\lambda\leq 1$. Geometrically they are part of the outer sphere,
because they have distance $\sqrt{2}$ from the origin.

As described in the previous section, the lower part of our self-dual model
$Q_{3}^{(2)}$ consists of three half-circular cones, each with a basis state
as their apex, which end at a half-sized point mirror image of the pure-state
spherical section, that is, at the insphere. This description matches exactly
the one of the internal boundary separating the points corresponding to rank-2
states from those corresponding only to rank-3 states in the model
$Q_{3}^{(1)}$, except that it is below the base triangle rather than above.
Also in the previous section, it was established that rank-2 states are mapped
to the lower half exactly if doing so does not locate them outside the model.

This allows us to determine the image of the set of rank-2 states in
$Q_{3}^{(2)}$. The only rank-2 states that fit into the lower part are exactly
those which in $Q_{3}^{(1)}$ lie on the internal boundary between rank-2
states and rank-3 only states. Therefore all rank-2 states that get mapped
into the lower half of the self-dual model are mapped onto the model’s
surface, while the interior of the lower part consists only of rank-3 states,
just as criterion v) requires.

However, the same is not true in the upper part. Almost all the rank-2 states
in the interior of ${Q}_{3}^{(1)}$ are not moved to the lower part and
therefore are also in the interior of ${Q}_{3}^{(2)}$. Only the points on the
internal rank boundary in $Q_{3}^{(1)}$ are moved to the boundary of
$Q_{3}^{(2)}$, leaving only rank-3 states on the internal boundary. That is,
criterion v) is still not perfectly satisfied. However, the new model
$Q_{3}^{(2)}$ is self-dual by construction and hence obeys the criterion vi).

The last criterion vii) is critical because also the self-dual model
$Q_{3}^{(2)}$ cannot reproduce these properties. This can be observed in the
highlighted cross section in Fig. 3. The point corresponding to
$\frac{1}{2}(\left|{0}\right\rangle\\!\left\langle{0}\right|+\left|{1}\right\rangle\\!\left\langle{1}\right|)$
is non-exposed and counts as a non-exposed face of the two-dimensional cross
section [18]. As these non-exposed faces occur at points where the self-dual
extension is attached to $Q_{3}^{(1)}$ one might speculate that their
occurrence is a price to pay for having such an extension. Moreover, the
highlighted cross section in Fig. 3 coincides with the projection of
$Q_{3}^{(2)}$ in the direction of the $z_{1}$ axis. Neither the corner
corresponding to $\left|{2}\right\rangle\\!\left\langle{2}\right|$ nor that of
$\left|{+_{2}}\right\rangle\\!\left\langle{+_{2}}\right|$ is polyhedral.

The self-duality of $Q_{3}^{(2)}$ is particularly interesting, because it
enables us to analyze how findings from previous work are featured in our
model. As a first example we consider the duality of boundary states on outer
and inner spheres. As was elucidated in Refs. [11, 20], boundary states on the
outer sphere have their dual counterparts on the inner sphere, and vice versa.
This property can be observed in Figs. 3 and 4 and is guaranteed by
construction via the universal state inversion, Eq. (6), (7).

It is noteworthy that the universal state inversion (or reduction [31]) map
finds an explicit geometric application here. Until now it was mainly
associated with entanglement properties of multi-party states, although the
relevance of geometrical concepts was implicit also in that context (cf.,
e.g., [32]).

One of the most notable qutrit three-sections is the so-called obese
tetrahedron (or three-dimensional elliptope). It was already noted by Bloore
[7] and thoroughly studied by Goyal and co-workers [20] (see Sec. 5.4 and Fig.
5 in Ref. [20]). Its barycenter is the completely mixed state and the corners
are given by the four non-orthogonal states
$\frac{1}{\sqrt{3}}(\left|{0}\right\rangle\pm\left|{1}\right\rangle\pm\left|{2}\right\rangle)$.
The faces are bulgy, because the tetrahedron contains with its corner points
also their duals on the opposite side. The image of the tetrahedron in
$Q_{3}^{(2)}$ is the part of the vertical axis belonging to the model: The
corner points are all mapped to the north pole, the barycenter to the origin
and the centers of the faces to the point with $z_{1}=z_{2}=0$ and
$w=-\frac{1}{\sqrt{2}}$.

The last example we mention here is the conical three-section, Sec. 5.1 and
Fig. 2 in Ref. [20]. Instead of the [128] section we consider the (unitarily
equivalent) [138] section which indeed has a three-dimensional image in
$Q_{3}^{(2)}$. The states of the cone are formed by all convex combinations of
the real states of the $\\{\left|{0}\right\rangle,\left|{1}\right\rangle\\}$
Bloch ball (a circular disk) and the state
$\left|{2}\right\rangle\\!\left\langle{2}\right|$. Therefore, one might expect
the image to include for $w>0$ the corresponding half-cone with the
semicircular surface at $z_{2}=\frac{1}{\sqrt{2}}$ and the apex at
$z_{2}=-\sqrt{2}$. However, this is only partially correct, because a part of
the half-cone gets mapped to values $w<0$ according to our rule, Eq. (11).
This includes the complete half-cone attached to
$\left|{2}\right\rangle\\!\left\langle{2}\right|$ for values $w<0$ and
$z_{2}\leq-\frac{\sqrt{2}}{4}$. Moreover, also the half-cone with apex in the
origin and base at $z_{2}=-\frac{\sqrt{2}}{4}$ with base radius
$\frac{\sqrt{6}}{4}$ gets mapped to $w<0$. The union of both parts at $w<0$
and $w\geq 0$ corresponds to the half-cone mentioned before, that is, to one
half of the actual conic three-section.

##  5 Quantum channels

In order to give yet another demonstration that our method directly connects
to the well-known properties of the Bloch sphere/ball visualization we show
the action of three standard decohering channels (cf. Ref. [1]) on a qutrit.
The depolarizing channel corresponds to driving a state towards the completely
mixed state $\frac{1}{3}\mathds{1}_{3}$, whereas phase-damping amounts to
mixing a state with its diagonal part. Amplitude-damping describes the
relaxation to one of the basis states. Following the discussion of
$Q^{(1)}_{3}$ the action of these channels on the qutrit Bloch body is evident
and completely analogous to what is known for the $d=2$ Bloch ball (cf. Fig.
5).

The channels are characterized by their Kraus operators $K_{j}$ that obey the
relation $\sum_{j}K_{j}^{\dagger}K_{j}=\mathds{1}_{3}$. For the depolarizing
channel we have $K_{0}=\sqrt{1-8\gamma/9}\ \mathds{1}_{3}$,
$K_{j}=\sqrt{\gamma/9}\ \mathsf{h}_{j}$, where $\ \mathsf{h}_{j}$ $(j=1\ldots
8)$ stands as a shortcut for all the Gell-Mann matrices, so that

$\displaystyle\rho^{\prime}(\gamma)\ =\ (1-\gamma)\ \rho\ +\ \gamma\
\mathds{1}_{3}\ \ .$ (13)

The parameter $\gamma\in[0,1]$ describes the strength of the channel action.

The Kraus operators for phase damping (or pure dephasing) in its simplest
version are $K_{0}=\sqrt{1-2\gamma/3}\ \mathds{1}_{3}$,
$K_{1}=\sqrt{\gamma/3}\ \operatorname{\mathsf{Z}}_{1}$,
$K_{2}=\sqrt{\gamma/3}\ \operatorname{\mathsf{Z}}_{2}$.

Under the action of the depolarizing and the phase-damping channel, the
shrinking of the Bloch body model $Q_{3}^{(1)}$ occurs proportionally along
straight lines, and the images of states that are represented by the same
point will again be represented by the same point in the shrunken model. The
same is, however, not true for the amplitude-damping channel.

The amplitude damping channel considered is of the form given by Grassl et al.
[33], with $\gamma_{01}=\gamma_{02}=:\gamma$ and $\gamma_{12}=0$, so that the
states $\lvert 1\rangle$ and $\lvert 2\rangle$ are damped equally. That is,

$\displaystyle K_{0}\ =\ $
$\displaystyle\left|{0}\right\rangle\\!\left\langle{0}\right|+\sqrt{1-\gamma}\left|{1}\right\rangle\\!\left\langle{1}\right|+\sqrt{1-\gamma}\left|{2}\right\rangle\\!\left\langle{2}\right|$
$\displaystyle K_{1}\ =\ $
$\displaystyle\sqrt{\gamma}\left|{0}\right\rangle\\!\left\langle{1}\right|$
(14) $\displaystyle K_{2}\ =\ $
$\displaystyle\sqrt{\gamma}\left|{0}\right\rangle\\!\left\langle{2}\right|\ \
.$

Figure 5: Action of quantum channels on the qutrit state space in
$Q^{(1)}_{3}$ representation: a) depolarizing, b) phase-damping, and c)
amplitude-damping channel.

This results in non-equal scaling of the non-diagonal elements. Whereas
$\rho_{01}$ and $\rho_{02}$ scale with a prefactor of $\sqrt{1-\gamma}$, the
element $\rho_{12}$ gets a prefactor of $(1-\gamma)$.

This does not matter for diagonal states ($w=0$), which are mapped into the
model uniquely, and which are therefore again scaled linearly, by

$\displaystyle z_{1}$ $\displaystyle\mapsto
z_{1}^{\prime}=\sqrt{\frac{3}{2}}+\left(z_{1}-\sqrt{\frac{3}{2}}\right)(1-\gamma)$
(15) $\displaystyle z_{2}$ $\displaystyle\mapsto
z_{2}^{\prime}=\frac{1}{\sqrt{2}}+\left(z_{2}-\frac{1}{\sqrt{2}}\right)(1-\gamma)$
(16)

Note that these two coordinates are scaled the same way for all the states.

For pure states, the absolute value of the offdiagonal elements is uniquely
determined by the diagonal elements, and since only the absolute value enters
into the $w$ coordinate, again each point of the original model gets mapped to
only one point of the shrunken model. In particular, the $w$ coordinate of
pure states gets mapped according to

$\displaystyle w^{\prime 2}=(1-\gamma)^{2}$
$\displaystyle\left[\frac{2}{3}\left(1-z_{2}^{2}-\frac{z_{2}}{\sqrt{2}}\right)\
+\ z_{1}\left(\frac{2}{\sqrt{3}}z_{2}-\sqrt{\frac{2}{3}}\right)\right]$
$\displaystyle\ +\ (1-\gamma)\ \left[\frac{4}{3}\ -\
z_{1}^{2}+\left(\sqrt{\frac{2}{3}}-\frac{2}{\sqrt{3}}z_{2}\right)z_{1}-\frac{z_{2}^{2}}{3}+\frac{\sqrt{2}z_{2}}{3}\right]\
\ .$ (17)

For non-diagonal mixed states, states with different distribution of absolute
values in the offdiagonal matrix elements get mapped to the same point in the
model. Therefore states that are represented by the same point will get mapped
by the channel to states represented by different points. It is obvious that
the value of $w$ will always be mapped somewhere into the range
$[(1-\gamma)^{2}w,(1-\gamma)w)$, but that range will not be exhausted for all
source points.

##  6 Conclusion

We have presented two ways of visualizing the eight-dimensional state space of
a qutrit in $\mathbb{R}^{3}$, thereby preserving a number of essential
geometric properties of $\mathcal{Q}_{3}$, cf. Ref. [18]. Depending on the
context one or the other representation may appear more useful. Our findings
are relevant, because the properties of $Q^{(1)}_{3}$ and $Q^{(2)}_{3}$ are
much closer to those of arbitrary higher-dimensional state space than the
Bloch ball for qubits.

Thus we hope our results are helpful to develop a more precise intuition for
objects in higher dimensions, which are commonly encountered in quantum
information science.

It is conceivable that there are more possibilities for such visualizations
along the lines of this work. The essential features of this visualization
persist even if we use it to represent state spaces of dimensions $d>3$
(splitting the generalized Gell-Mann into a diagonal and off-diagonal part in
an identical fashion). Finding the most useful geometrical representation of
those and categorizing them is a direction that we believe to yield further
fruitful insight in the future.

This work was funded by the German Research Foundation Project EL710/2-1
(C.E., J.S.), by grant PGC2018-101355-B-100 (MCIU/FEDER,UE) and Basque
Government grant IT986-16 (J.S.). C.E. and J.S. acknowledge Klaus Richter’s
support for this project. M.H. and S.M. acknowledge funding from the Austrian
Science Fund (FWF) through the START project Y879-N27 and the stand-alone
project P 31339-N27.

##  7 Appendix

For completeness, we provide the definition for the Gell-Mann matrices that we
use in Eqs. (2a) and (9). Please note that we use a normalization to the
dimension $d=3$ (as opposed to the normalization to 2 conventionally used in
high-energy physics). Moreover, also the enumeration differs from the usual
one and is adapted to our coordinate description:

$\displaystyle\operatorname{\mathsf{X}}_{1}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}0&1&0\\\ 1&0&0\\\
0&0&0\end{array}\right]\ \ $ $\displaystyle\operatorname{\mathsf{Y}}_{1}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}0&-\mathrm{i}&0\\\
\mathrm{i}&0&0\\\ 0&0&0\end{array}\right]$
$\displaystyle\operatorname{\mathsf{X}}_{2}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}0&0&1\\\ 0&0&0\\\
1&0&0\end{array}\right]\ \ $ $\displaystyle\operatorname{\mathsf{Y}}_{2}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}0&0&-\mathrm{i}\\\
0&0&0\\\ \mathrm{i}&0&0\end{array}\right]$
$\displaystyle\operatorname{\mathsf{X}}_{3}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}0&0&0\\\ 0&0&1\\\
0&1&0\end{array}\right]\ \ $ $\displaystyle\operatorname{\mathsf{Y}}_{3}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}0&0&0\\\
0&0&-\mathrm{i}\\\ 0&\mathrm{i}&0\end{array}\right]$
$\displaystyle\operatorname{\mathsf{Z}}_{1}\ =\ $
$\displaystyle\sqrt{\frac{3}{2}}\left[\begin{array}[]{ccc}1&0&0\\\ 0&-1&0\\\
0&0&0\end{array}\right]\ \ $ $\displaystyle\operatorname{\mathsf{Z}}_{2}\ =\ $
$\displaystyle\frac{1}{\sqrt{2}}\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&-2\end{array}\right]\ \ .$

## References

  * Nielsen and Chuang [2000] M.A. Nielsen and I.L. Chuang,  _Quantum Computation and Quantum Information_ (Cambridge University Press, 2000).
  * Martin  _et al._ [2010] J. Martin, O. Giraud, P.A. Braun, D. Braun, and T. Bastin,  _Multiqubit symmetric states with high geometric entanglement,_Physical Review A 81, 062347 (2010).
  * Fano [1954] U. Fano,  _A Stokes-Parameter Technique for the Treatment of Polarization in Quantum Mechanics,_Phys. Rev. 93, 121 (1954).
  * Fano [1957] U. Fano,  _Description of States in Quantum Mechanics by Density Matrix and Operator Techniques,_Rev. Mod. Phys. 29, 74 (1957).
  * Feynman  _et al._ [1957] R.P. Feynman, F.L. Vernon, and R.W. Hellwarth,  _Geometrical Representation of the Schrödinger Equation for Solving Maser Problems,_J. Appl. Phys. 28, 49 (1957).
  * Bengtsson and Życzkowski [2017] I. Bengtsson and K. Życzkowski,  _Geometry of Quantum States: An Introduction to Quantum Entanglement, 2nd edition_ (Cambridge University Press, 2017).
  * Bloore [1976] F.J. Bloore,  _Geometrical description of the convex sets of states for systems with spin-1/2 and spin-1,_J. Phys. A 9, 2059 (1976).
  * Ramachandran and Murthy [1979] G. Ramachandran and M.V.N. Murthy,  _A new representation for the density matrix,_Nucl. Phys. A 323, 403 (1979).
  * Byrd and Khaneja [2003] M.S. Byrd and N. Khaneja,  _Characterization of the positivity of the density matrix in terms of the coherence vector representation,_Phys. Rev.A 68, 062322 (2003).
  * Kimura [2003] G. Kimura,  _The Bloch vector for N-level systems,_Phys. Lett. A 314, 339 (2003).
  * Kimura and Kossakowski [2005] G. Kimura and A. Kossakowski,  _The Bloch-vector space for $N$-level systems: the spherical-coordinate point of view, _Open Syst. Inf. Dyn. 12, 207 (2005).
  * Mendaš [2006] I.P. Mendaš,  _The classification of three-parameter density matrices for a qutrit,_J. Phys. A: Math. .Gen. 39, 11313 (2006).
  * Boya and Dixit [2008] L.J. Boya and K. Dixit,  _Geometry of density matrix states,_Phys. Rev. A 78, 042108 (2008).
  * Vinjanampathy and Rau [2009] S. Vinjanampathy and A.R.P. Rau,  _Bloch sphere like construction of SU(3) Hamiltonians using Unitary Integration,_J. Phys. A: Math. Theor. 42, 425303 (2009).
  * Dunkl  _et al._ [2011] C.F. Dunkl, P. Gawron, J.A. Holbrook, J.A. Miszczak, Z. Puchała, and K. Życzkowski,  _Numerical shadow and geometry of quantum states,_J. Phys. A: Math. Theor. 44, 335301 (2011).
  * Kurzyński [2011] P. Kurzyński,  _Multi-Bloch vector representation of the qutrit,_Quantum Inf. Comp. 11, 361 (2011).
  * Sarbicki and Bengtsson [2012] G. Sarbicki and I. Bengtsson,  _Dissecting the qutrit,_J. Phys. A: Math. Theor. 46, 035306 (2012).
  * Bengtsson  _et al._ [2012] I. Bengtsson, S. Weis, and K. Życzkowski,  _Geometry of the Set of Mixed Quantum States: An Apophatic Approach_ (Springer Basel, 2012) p. 175.
  * Tabia and Appleby [2013] G.N.M. Tabia and D.M. Appleby,  _Exploring the geometry of qutrit state space using symmetric informationally complete probabilities,_Phys. Rev. A 88, 012131 (2013).
  * Goyal  _et al._ [2016] S. Goyal, B.N. Simon, R. Singh, and S. Simon,  _Geometry of the generalized Bloch sphere for qutrits,_J. Phys. A: Math. Theor. 49, 165203 (2016).
  * Kurzyński  _et al._ [2016] P. Kurzyński, A. Kołodziejski, W. Laskowski, and M. Markiewicz,  _Three-dimensional visualization of a qutrit,_Phys. Rev. A 93, 062126 (2016).
  * Szymański  _et al._ [2018] K. Szymański, S. Weis, and K. Życzkowski,  _Classification of joint numerical ranges of three hermitian matrices of size three,_Lin. Alg. Appl. 545, 148 (2018).
  * Xie  _et al._ [2020] J. Xie, A. Zhang, N. Cao, H. Xu, K. Zheng, Y.T. Poon, N.S. Sze, P. Xu, B. Zeng, and L. Zhang,  _Observing Geometry of Quantum States in a Three-Level System,_Phys. Rev. Lett. 125, 150401 (2020).
  * Dogra  _et al._ [2020] S. Dogra, A. Vepsäläinen, and G.S. Paraoanu,  _Majorana representation of adiabatic and superadiabatic processes in three-level systems,_Phys. Rev. Research 2, 043079 (2020).
  * Sharma and Ghosh [2021] G. Sharma and S. Ghosh,  _Four-dimensional Bloch sphere representation of qutrits using Heisenberg-Weyl Operators,_arXiv:2101.06408 [quant-ph] (2021).
  * Serrano-Ensástiga and Braun [2020] E. Serrano-Ensástiga and D. Braun,  _Majorana representation for mixed states,_Phys. Rev. A 101, 022332 (2020).
  * Bertlmann and Krammer [2008] R.A. Bertlmann and P. Krammer,  _Bloch vectors for qudits,_Journal of Physics A: Mathematical and Theoretical 41, 235303 (2008).
  * [28] Note that we use a different normalization compared to that in Ref. [18] because of a different choice of prefactor for the Hilbert-Schmidt norm.
  * Baumgratz  _et al._ [2014] T. Baumgratz, M. Cramer, and M.B. Plenio,  _Quantifying Coherence,_Phys. Rev. Lett. 113, 140401 (2014).
  * Rungta  _et al._ [2001] P. Rungta, V. Bužek, C.M. Caves, M. Hillery, and G.J. Milburn,  _Universal state inversion and concurrence in arbitrary dimensions,_Phys. Rev. A 64, 042315 (2001).
  * Horodecki and Horodecki [1999] M. Horodecki and P. Horodecki,  _Reduction criterion of separability and limits for a class of distillation protocols,_Phys. Rev. A 59, 4206 (1999).
  * Eltschka and Siewert [2018] C. Eltschka and J. Siewert,  _Distribution of entanglement and correlations in all finite dimensions,_Quantum 2, 64 (2018).
  * Grassl  _et al._ [2018] M. Grassl, L. Kong, Z. Wei, Z. Yin, and B. Zeng,  _Quantum Error-Correcting Codes for Qudit Amplitude Damping,_IEEE Transactions on Information Theory 64, 4674 (2018).


# Should we use model-free or model-based control? A case study of battery
management systems
††thanks: Mohamad Fares El Hajj Chehade and Young-ho Cho contribute equally to
this paper. (Corresponding author: Mohamad Fares El Hajj Chehade.)

Mohamad Fares El Hajj Chehade*, Young-ho Cho*, Sandeep Chinchali, and Hao Zhu
Chandra Department of Electrical and Computer Engineering
The University of Texas at Austin
Austin, TX, USA
{chehade, jacobcho, sandeepc<EMAIL_ADDRESS>
###### Abstract

Reinforcement learning (RL) and model predictive control (MPC) each offer
distinct advantages and limitations when applied to control problems in power
and energy systems. Despite various studies on these methods, benchmarks
remain lacking and the preference for RL over traditional controls is not well
understood. In this work, we put forth a comparative analysis using RL- and
MPC-based controllers for optimizing a battery management system (BMS). The
BMS problem aims to minimize costs while adhering to operational limits. by
adjusting the battery (dis)charging in response to fluctuating electricity
prices over a time horizon. The MPC controller uses a learning-based forecast
of future demand and price changes to formulate a multi-period linear program,
that can be solved using off-the-shelf solvers. Meanwhile, the RL controller
requires no time-series modeling but instead is trained from the sample
trajectories using the proximal policy optimization (PPO) algorithm. Numerical
tests compare these controllers across optimality, training time, testing
time, and robustness, providing a comprehensive evaluation of their efficacy.
RL not only yields optimal solutions quickly but also ensures robustness to
shifts in customer behavior, such as changes in demand distribution. However,
as expected, training the RL agent is more time-consuming than MPC.

###### Index Terms:

Battery management system, reinforcement learning, model predictive control,
time-series forecast.

##  I Introduction

In the domain of control problems, reinforcement learning (RL) and model-
predictive control (MPC) have their merits and drawbacks. RL deals with agents
in uncertain settings and optimizes the actions over infinite horizons without
modeling the dynamics of the environment. It excels in environments with large
modeling errors but suffers from data inefficiency and safety due to possibly
undesirable actions. On the other hand, MPC optimizes actions over a finite
horizon, with heavy reliance on accurate environmental models. It is less
concerned with safety and data efficiency, due to the deterministic nature of
the optimization problem and its use of hard constraints. However, it could
perform poorly due to modeling inaccuracies and the limitations of a finite
horizon[1, 2, 3].

In power and energy systems, both RL and MPC have been explored for the
battery control problem, because of its sequential decision-making nature. In
the problem of battery management system (BMS), the charging and discharging
schedule of an energy storage system is controlled to minimize the cost of
electricity purchase over a certain period. For example, one work based on MPC
uses a long short-term memory (LSTM) forecast to model the system; then, the
multi-period optimization problem is solved [4]. Despite showing its
superiority over conventional methods like dynamic programming and fuzzy
logic, this work does not check the impact of forecasting errors on the
optimality of MPC and other methods. On the other hand, RL in [5],
particularly deep-Q networks (DQN), has improved the battery’s operation and
minimized its degradation. However, there is a lack of benchmark in general
for justifying the choice of RL over currently used control methods.

In this work, a comparative analysis is conducted between the model-based MPC
and the model-free RL approaches to the battery control problem. The closest
to our work is [6], which compares the performance of DQN with one-step MPC.
Nonetheless, the comparisons therein are limited in three aspects. First, DQNs
have been outperformed by several policy-based algorithms, such as proximal
policy optimization (PPO) and deep deterministic policy gradient (DDPG), which
are currently the popular RL algorithms in power and energy systems [7, 8].
Second, one-step MPC is essentially a greedy algorithm that does not benefit
from planning over multiple time steps, thereby compromising the optimality of
the MPC results. Third, the comparison is solely based on the control cost,
and does not include some important criteria, such as robustness and
computational efficiency. To the best of our knowledge, our work is the first
to address these three aspects and conduct a full comparative analysis between
RL and MPC for BMS. The main contributions of this work are summarized below:

  1. 1.

We study the optimality of the decisions made by RL and MPC for the battery
control problem. We adopt the well-known, widely used algorithms for each
approach.

  2. 2.

We extend the basis of the comparison from control cost to include multiple
other criteria: data efficiency, online computation time, and robustness.

  3. 3.

Our results show that RL performs better in terms of optimality and online
computation time. It is also more robust to shifts in the demand distribution.
However, it is much more data-inefficient than MPC, whose optimality is
drastically affected by forecasting errors.

  4. 4.

Our methodology and analysis go beyond the BMS problem and could be applicable
to general control problems in power and energy systems.

The rest of the paper is organized as follows: Section II defines the problem,
Sections III and IV respectively introduce the MPC and RL approaches, Section
V discusses the main results, and Section VI concludes the paper.

##  II Preliminaries

The battery control problem is formulated for a system of three components: a
load with a time-varying demand for electricity, denoted by $d_{t}$ (in kW)
for a given time $t$, an electric utility that sells electricity at a time-
varying price $\rho_{t}$ (in $/kWh), and a battery storage system of a
capacity $E$ (in kWh) and state of charge $\text{SOC}_{t}$. The latter is
capable of charging electricity from the utility and discharging to the load
at a rate $a_{t}=a_{\text{max}}$ (in kW), as long as $\text{SOC}_{t}$ remains
within operational limits $\text{SOC}_{\text{min}}$ and
$\text{SOC}_{\text{max}}$. The electricity tariff $g_{t}$ at a given time $t$
is then linear in $a_{t}$: $g_{t}=\rho_{t}(d_{t}+a_{t})$. The interaction
between the three components of the system is best illustrated in Fig. 1.

To minimize the cost of purchasing electricity, the battery must be discharged
during periods of high tariffs, and charged otherwise. However, due to the
battery’s limited capacity, the controller must plan over a certain time
horizon, to see when charging and discharging are most economical. This
results in a multi-step optimization problem with cost function $f$, where $f$
is the sum of individual tariffs at the considered timesteps:
$f=\sum_{t}g_{t}$. Despite the linearity of the cost function in terms of the
decision variables $a_{t}$, the problem is non-trivial due to the
uncertainties in the demand $d_{t}$ and prices $\rho_{t}$. As a result, two
approaches are considered for solving the battery control problem: model-based
and model-free.

Figure 1: The components of the system. At a given time $t$, the load has a
demand $d_{t}$ and the utility sells electricity at a price $\rho_{t}$. To
minimize the total cost of purchasing electricity over a certain period, the
battery controller chooses between charging from the grid ($a_{t}>0$) and
discharging to the load ($a_{t}<0$).

##  III Model-Predictive Control (MPC)

The MPC problem aims at minimizing the operating cost of the energy system
over a finite time horizon $t\in\mathcal{T}=\\{1,2,\ldots,T\\}$, while
maintaining the battery within its operational limits. At each time step
$t\in\mathcal{T}$, the first element $b_{t|t}$ in the vector of solutions
obtained from each optimization problem is stored. Subsequently, another
optimization problem is solved starting at the time step $t+1$. The demand and
price in the problem are sampled from two time-series distributions:
$\rho_{{t+k}|t}\sim f_{\rho}(t+k,\theta)$ and $d_{{t+k}|t}\sim f_{d}(t+k,\mu)$
for all $k\in\\{1,2,3,\ldots,T\\}$, where $\theta$ and $\mu$ are the
parameters of the price and demand distributions, respectively. We obtain
estimated values for price and demand from forecasters, i.e.,
$\hat{\rho}_{t+k}|t\sim f_{\hat{\rho}}(t+k,\hat{\theta})$ and
$\hat{d}_{t+k}|t\sim f_{\hat{d}}(t+k,\hat{\mu})$ for all $t\in\mathcal{T}$,
where $\mathcal{T}$ is the full horizon of the problem.

To tackle the inherent unpredictability of future demand and electricity
prices, we employ a predictive framework relying on an LSTM network, well-
suited for sophisticated time series forecasting. This neural network acts as
a forecaster, projecting demand and cost factors from the current time step
$t+1$ to a specified horizon $\mathcal{T}$. The MPC framework leverages the
forecasted data to make informed decisions about battery charging and
discharging, aiming to minimize electricity expenses and stabilize battery
energy levels. To enhance the reliability of these forecasts, we engage in
iterative refinement of the predictive models. This critical step involves
continuously aligning the models with observed data, thereby reducing
forecasting errors and improving decision-making accuracy in the MPC
framework.

Given the predicted electricity cost $\hat{\rho}_{t+k|t}$ at time $t+k$, which
is predicted at time $t$, the BMS optimization problem seeks to minimize the
total cost while satisfying the operational limits, as given by:

$\displaystyle\underset{a_{t|t},\ldots,a_{t+T|t}}{\text{min}}$
$\displaystyle\sum_{k=0}^{T}\hat{\rho}_{t+k|t}(E\cdot
a_{t+k|t}+\hat{d}_{t+k|t})$ (1a) subject to
$\displaystyle\text{SOC}_{t+k+1|t}=\text{SOC}_{t+k|t}+a_{t+k|t},$ (1b)
$\displaystyle\text{SOC}_{\min}\leq\text{SOC}_{t+k+1|t}\leq\text{SOC}_{\max},$
(1c) $\displaystyle\|a_{t+k|t}\|\leq a_{\text{max}},$ (1d)
$\displaystyle\forall k\in\\{0,\ldots,T\\},~{}\forall t\in\mathcal{T}.$ (1e)

The constraints in (1b) represent the remaining battery level,
$\text{SOC}_{t+k+1|t}$, at $t+k+1$ after making a charging or discharging
decision at $t+k$. The battery’s operational limits in (1c) range from
$\text{SOC}_{\min}$ to $\text{SOC}_{\max}$. Moreover, the maximum allowable
charge/discharge rate is limited by the constant $a_{\text{max}}$ (in kW), as
denoted in (1d). The problem (1) is a multi-step linear optimization that can
be efficiently solved by off-the-shelf solvers.

##  IV Reinforcement Learning (RL)

Due to its sequential decision-making nature and Markovian state space, the
problem can be modeled as a Markov decision process (MDP) [9]. An MDP is
defined as the tuple $(\mathcal{S},\mathcal{A},p,r,\gamma)$, where
$\mathcal{S}$ and $\mathcal{A}$ are sets of states and actions, respectively.
For an action $a\in\mathcal{A}$ that an agent takes in a given state
$s\in\mathcal{S}$, $p(.|s,a)$ models the distribution over the next states.
For the transition $s\xrightarrow{a}s^{\prime}$, the agent receives a scalar
reward $r(s,a,s^{\prime})$. Finally, the discount factor $\gamma\in[0,1]$ can
be used to give smaller weights for future rewards. The elements of the MDP
tuple for the BMS problem are defined below:

  * •

State space ($\mathcal{S}$): The state is characterized by the tuple
$(\text{SOC}_{t},\rho_{t},d_{t},h_{t},D_{t})$, where:

    * –

$\text{SOC}_{t}$ represents the state of charge of the battery at time $t$.

    * –

$\rho_{t}$ is the electricity price at time $t$.

    * –

$d_{t}$ indicates the energy demand at time $t$.

    * –

$h_{t}$ is the hour $t$ of the day.

    * –

$D_{t}$ distinguishes between weekdays and weekends.

  * •

Action space ($\mathcal{A}$): The action space is one-dimensional and
represents the control variable $a_{t}$, defined as the amount of
charge/discharge at a given time $t$. Since the optimal solution is always one
of three values: charge with the maximum amount ($a_{t}$ = $a_{\text{max}}$ =
0.1E), discharge with the maximum amount ($a_{t}$ = - $a_{\text{max}}$ = -
0.1E), or idle ($a_{t}$ = 0), the action space is discrete
$A=\\{-0.1E,0,0.1E\\}$.

  * •

Transition function ($p$): Time elements in the state space incur
deterministic transitions. The state of charge progression is described below:

$\displaystyle\text{SOC}_{t+1}=\begin{cases}\text{SOC}_{\min}&\text{if
}\text{SOC}_{t}+a_{t}<\text{SOC}_{\min}\\\ \text{SOC}_{\max}&\text{if
}\text{SOC}_{t}+a_{t}>\text{SOC}_{\max}\\\
\text{SOC}_{t}+a_{t}&\text{otherwise}\end{cases}$ (2)

This ensures that SOC stays within the operational limits $\text{SOC}_{\min}$
and $\text{SOC}_{\max}$. On the other hand, the electricity price and energy
demand are assumed to be drawn from two time series distributions:
$\rho_{t}\sim f_{\rho}$ and $d_{t}\sim f_{d}$. Modeling $f_{\rho}$ and $f_{d}$
is challenging and may incur significant errors. Fortunately, model-free RL
allows for solving the problem without any knowledge of the time-series
distributions.

  * •

Reward function ($r$): the reward is defined as the negative of the
instantaneous electricity cost of the energy consumed. The power consumed is
used to satisfy the load demand and to charge the battery:
$r_{t}=-\rho_{t}(d_{t}+a_{t}E)$.

  * •

Discount factor ($\gamma$): Since electricity costs at future time steps
matter as much as costs in the present, no discounting is applied, i.e.
$\gamma=1$.

In an MDP framework, an RL agent takes its actions based on a policy
distribution $\pi$, i.e. $a_{t}\sim\pi(\cdot|s_{t})$, where $\pi(a|s)$ is the
probability of choosing action $a$ from state $s$ under policy $\pi$. Given
that the agent is at state $s_{t}$ at a given time $t$, the return $G_{t}$ is
defined as the total reward the agent receives starting at time $t$ and until
the end of the trajectory, when taking actions based on policy $\pi$. The
expected value of the return is known as the value function $V^{\pi}(s_{t})$.
The goal of the RL problem is to find the optimal policy $\pi^{*}$, which
maximizes the expected total reward the agent receives over its trajectory,
i.e. $\pi^{*}=\max_{\pi}V^{\pi}(s)\>\forall s\in\mathcal{S}$.

###  IV-A Methodology

In modern RL, policy-based deep RL methods have enjoyed more success than
classical value-based methods [10, 11, 12]. In the discrete action-space
setting, the policy distribution is represented by a deep neural network with
parameters $\theta$, as shown in Fig. 2. The input vector to the network is
the state $s$, and the output vector consists of the probability masses
$\pi_{\theta}(a|s)$. For a given state $s_{t}$, the action is sampled from the
output vector distribution: $a_{t}\sim\pi_{\theta}(\cdot|s_{t})$. The optimal
policy $\pi^{*}=\pi_{\theta^{*}}$ is obtained by maximizing the expected
return received along its experienced trajectories. In other words, if the
agent experiences a trajectory $\tau$, which produces a large return
$G^{\tau}_{0}$, the actions that yield such a trajectory are enforced by the
policy network, i.e. their probabilities increase. The objective function is
formally represented below:

$\displaystyle\max_{\theta}J(\theta)$
$\displaystyle=\max_{\theta}\mathbb{E}_{\tau}[G_{0}\mid s_{0}]$ (3)
$\displaystyle=\max_{\theta}\sum_{\tau}P(\tau)G_{0}$ (4)

where $P(\tau)$ is the probability of trajectory $\tau$, and is a function of
the actions chosen by the agent, i.e. $\pi_{\theta}$ and the environment
transitions $p$.

Figure 2: The policy network $\pi_{\theta}$. The state vector $s_{t}$
propagates through the neural network layers to output the probability masses
$\pi_{\theta}(a|s_{t})\>\forall a\in\mathcal{A}$. $\sigma$ is a non-linear
activation function to allow the network to represent non-linear input-output
relationships. A softmax function is applied at the output to convert the raw
numbers (logits) to probabilities.

This can be achieved with gradient ascent in the $\theta$ space:

$\displaystyle\theta$
$\displaystyle\leftarrow\theta+\alpha\sum_{t=0}^{T}\nabla_{\theta}\log\pi_{\theta}(a_{t}\mid
s_{t})G_{t}$ (5)

where $\alpha$ is the step size, and $T$ is the number of steps taken in the
trajectory.

To reduce the variance in the return during learning, an action-independent
baseline is subtracted from the return in (5). In most settings, this baseline
is the value function $V(s_{t})$, and the objective function becomes:

$\displaystyle\theta$
$\displaystyle\leftarrow\theta+\alpha\sum_{t=0}^{T}\nabla_{\theta}\log\pi_{\theta}(a_{t}\mid
s_{t})\left(G_{t}-V(s_{t})\right)$ (6)

The value function itself can be represented as a separate neural network with
parameter $\phi$, and is optimized by minimizing the temporal difference (TD)
error:

$\displaystyle\min_{\phi}\left(r+\gamma
V_{\phi}(s^{\prime})-V_{\phi}(s)\right)^{2}$ (7)

In some situations, large parameter updates result in unstable learning and
potential divergence. As a remedy, these updates can be ensured to be within a
trust region, by enforcing the KL-divergence between the parameters before and
after the update, $D_{KL}(\theta_{new}||\theta_{old})$, to be within a
threshold $\delta$. Therefore, the parameters take the largest step (towards
the steepest ascent) that keeps them inside trust region. This translates to
the constrained optimization problem below:

$\displaystyle\theta_{\text{new}}=\arg\max_{\theta}$ $\displaystyle\quad
J(\theta)$ (8) s.t. $\displaystyle\quad
D_{\text{KL}}(\theta\parallel\theta_{\text{old}})\leq\delta$

An approximate form of this problem is solved by the proximal policy
optimization (PPO) algorithm [13], whose pseudocode is shown in Algorithm 1.

!t

Algorithm 1 Proximal Policy Optimization (PPO)

1:Initialize parameters $\theta$ for policy $\pi_{\theta}$ and $\phi$ for
value function $V_{\phi}$

2:repeat

3: Collect trajectories $D_{k}=\\{s_{t},a_{t},r_{t},s^{\prime}_{t}\\}$ using
$\pi_{\theta}$

4: Compute $G_{t}$, $A_{t}=G_{t}-V_{\phi}(s_{t})$

5: for each epoch over $D_{k}$ do

6: Update $\theta$ by (approximately) solving (8)

7: Update $\phi$ minimizing $\left(r_{t}+\gamma
V_{\phi}(s^{\prime}_{t})-V_{\phi}(s_{t})\right)^{2}$

8: end for

9:until convergence

##  V Numerical Comparisons

In this section, we compare the performance of reinforcement learning (RL) and
model-predictive control (MPC) on a residential building. We consider two
stages in the evaluation: training and testing. In the training phase, the RL
agent learns a policy network, as outlined in Algorithm 1, while a forecaster
for price and demand is trained for MPC. In the testing phase, RL and MPC are
run on a new data set. The performance of each method is judged based on four
criteria:

  1. 1.

control performance (optimality): the total cost incurred on the test dataset.

  2. 2.

data efficiency: the data needed to train the RL agent or forecaster

  3. 3.

testing time: the time needed to make control decisions for the test dataset

  4. 4.

robustness: the degree of optimality when the distribution of the model
changes from the training dataset to the testing dataset

For this purpose, two datasets (each with training and testing sets) are
considered. In the first dataset, the model of the system remains the same
over the training and testing sets, and the first three criteria are tested.
In the second dataset, the distribution of the demand changes from the
training to the testing set, therefore allowing us to examine the robustness
of each controller.

As a baseline, we consider a pre-defined policy that discharges the battery
when the price is above average (computed based on the training dataset) and
charges the battery otherwise - as long as the state of charge is within
operational limits.

In both experiments, the data is hourly, and the horizon for MPC is 24 hours.
An LSTM network is trained in Tensorflow as the forecaster, while
StableBaselines3 [14] with default parameters is used to train the RL agent.
Five runs are performed for each method, for which a 95% confidence interval
of the cost is computed. The simulations took place on a machine equipped with
an Intel® CPU @ 2.10 GHz and 32 GB RAM.

TABLE I: The performance metrics for the different controllers. The ground
truth is computed by considering the actual values of demand and prices and
optimizing over the entire test set. RL performs best in terms of optimality
and testing time, while MPC is more data-efficient. The optimality of MPC
(with true values of demand and price) indicates that errors in the regular
MPC model are due to the forecaster.

(a) Optimality (control cost) Controller | Cost ($) | Optimality Gap
---|---|---
RL | $\textbf{85,100}\mathbf{\pm}\textbf{200}$ | 7.36%
MPC | $89,650\pm 665$ | $13.10\%$
MPC (exact model) | $79,268$ | $0\%$
Baseline | $90,989$ | $14.77\%$
Ground truth | $79,268$ | -
No BMS | $114,065$ | $43.89\%$

(b) Data Used and Testing Time Controller | Data Used | Testing Time (s)
---|---|---
RL | $3\times 10^{6}$ | $5$
MPC | $3\times 10^{3}$ | $30$

###  V-A Dataset 1: Consistent Model through Training and Testing Sets

The load corresponds to a residential building, and the real-time (RT) price
are from the Pennsylvania, New Jersey, and Maryland (PJM) region. The training
and testing sets are for July 2017 and 2018 respectively.

Table I shows the performance of each controller based on the first three
criteria. While the table shows that any control method that employs a battery
storage system is economic (as compared with no BMS), RL achieves the highest
degree of optimality, as compared with MPC and the baseline. This reinforces
the capabilities of trust-region policy-based methods in reaching the optimal
solution without any model of the environment [13, 10]. On the other hand, MPC
suffers from modeling errors incurred by the forecaster and performs just
better than the baseline pre-defined policy. The errors in the forecast are
verified by the optimality of the MPC solution when considering the actual
values for demand and price, rather than the forecasted ones.

In the meantime, RL requires about 1000 times more data than MPC, as the RL
agent requires a lot of interactions with the environment during training [2].
Nonetheless, training is performed only once, after which RL is 6 times faster
during the testing phase. This is because decision-making for RL involves a
forward propagation over the policy network (Fig. 2), while MPC has to solve a
multi-step optimization problem every time a decision is made.

###  V-B Dataset 2: Distributional Shifts in the Demand from the Training to
the Testing Sets

The training and testing data are for July and August of 2017, respectively.
The price data is again retrieved from PJM, while the demand is for a
commercial building in South Korea. Fig. 3 illustrates the differences in the
distributions of demand between training and testing.

Table II displays the total cost incurred by each controller, where RL is
shown to maintain its high level of optimality. This not only emphasizes the
advantages of model-free learning utilized by RL, but also RL’s ability to
learn general rules that transcend changes in the model. On the other hand,
the errors in the forecaster are magnified by this distributional shift, which
resulted in a relatively high control cost for MPC.

Figure 3: The shift in the distribution of the demand from the training to the testing datasets. TABLE II: The optimality of each controller when considering distributional shifts in the demand from training to testing. RL is robust to such shifts, as it maintains a high level of optimality due to its generalization capabilities. The performance of MPC is again drastically affected by forecasting errors. Controller | Cost ($) | Optimality Gap
---|---|---
RL | $\textbf{158,977}\mathbf{\pm}\textbf{389}$ | 2.58%
MPC | $161,806\pm 740$ | $4.40\%$
MPC (exact model) | $154,981$ | $0\%$
Baseline | $166,490$ | $7.43\%$
Ground Truth | $154,981$ | -
No BMS | $179,958$ | $16.12\%$

##  VI Conclusion

To sum up, we have developed and compared the model-based MPC and model-free
RL controllers for the battery management system. Although RL requires longer
training time compared to MPC, it can achieve optimal control and offers
faster testing time than MPC. We summarize the general comparison of RL and
MPC for optimal control in Table III.

TABLE III: Comparison of RL and MPC for the optimal control of the battery storage system. RL yields more optimal solutions in less time. Furthermore, its solutions are robust to changes in customer behavior, represented by shifts in the demand distribution. On the other hand, training the RL agent requires much more time than MPC. Comparison criteria | RL | MPC
---|---|---
Optimality | ✓ | ✗
Data Efficiency | ✗ | ✓
Online Computation | ✓ | ✗
Robustness | ✓ | ✗

Our future research directions include considering more generalized energy
management systems with photovoltaic panels and multiple loads, as well as
implementing different types of forecasters to handle uncertainties arising
from those components.

## References

  * [1] H. Zhang, S. Seal, D. Wu, F. Bouffard, and B. Boulet, “Building energy management with reinforcement learning and model predictive control: A survey,” _IEEE Access_ , vol. 10, pp. 27 853–27 862, 2022.
  * [2] S. Kamthe and M. Deisenroth, “Data-efficient reinforcement learning with probabilistic model predictive control,” in _International conference on artificial intelligence and statistics_. PMLR, 2018, pp. 1701–1710.
  * [3] Y. Lin, J. McPhee, and N. L. Azad, “Comparison of deep reinforcement learning and model predictive control for adaptive cruise control,” _IEEE Transactions on Intelligent Vehicles_ , vol. 6, no. 2, pp. 221–231, 2020.
  * [4] H. Chen, R. Xiong, C. Lin, and W. Shen, “Model predictive control based real-time energy management for hybrid energy storage system,” _CSEE Journal of Power and Energy Systems_ , vol. 7, no. 4, pp. 862–874, 2020.
  * [5] K.-b. Kwon and H. Zhu, “Reinforcement learning-based optimal battery control under cycle-based degradation cost,” _IEEE Transactions on Smart Grid_ , vol. 13, no. 6, pp. 4909–4917, 2022.
  * [6] A. S. Zamzam, B. Yang, and N. D. Sidiropoulos, “Energy storage management via deep Q-networks,” in _2019 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2019, pp. 1–5.
  * [7] J. Feng, Y. Shi, G. Qu, S. H. Low, A. Anandkumar, and A. Wierman, “Stability constrained reinforcement learning for decentralized real-time voltage control,” _IEEE Transactions on Control of Network Systems_ , 2023.
  * [8] Y. Chen, Y. Shi, D. Arnold, and S. Peisert, “SAVER: Safe learning-based controller for real-time voltage regulation,” in _2022 IEEE Power & Energy Society General Meeting (PESGM)_. IEEE, 2022, pp. 1–5.
  * [9] M. L. Puterman, _Markov decision processes: discrete stochastic dynamic programming_. John Wiley & Sons, 2014\.
  * [10] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust region policy optimization,” in _International conference on machine learning_. PMLR, 2015, pp. 1889–1897.
  * [11] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu, “Asynchronous methods for deep reinforcement learning,” in _International conference on machine learning_. PMLR, 2016, pp. 1928–1937.
  * [12] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International conference on machine learning_. PMLR, 2018, pp. 1861–1870.
  * [13] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017\.
  * [14] Stable-Baselines3, “Stable-Baselines3 Docs - Reliable Reinforcement Learning Implementations — Stable Baselines3 2.2.1 documentation,” Online, Accessed 2023. [Online]. Available: https://stable-baselines3.readthedocs.io/en/master/