#  Measures of Entropy and Complexity in altered states of consciousness

D. M. Mateos${}^{1}*$, R. Guevara Erra2, R. Wennberg3, J.L. Perez Velazquez1
1 Neuroscience and Mental Health Programme, Division of Neurology, Hospital
for Sick Children.
Institute of Medical Science and Department of Paediatrics, University of
Toronto, Toronto, Canada.
2 Laboratoire Psychologie de la Perception, CNRS and Université Paris
Descartes, Sorbonne Paris Cité, Paris, France.
3 Krembil Neuroscience Centre, Toronto Western Hospital, University of
Toronto, Toronto, Canada.
*<EMAIL_ADDRESS>
###### Abstract

Quantification of complexity in neurophysiological signals has been studied
using different methods, especially those from information or dynamical system
theory. These studies revealed the dependence on different states of
consciousness, particularly that wakefulness is characterized by larger
complexity of brain signals perhaps due to the necessity of the brain to
handle varied sensorimotor information. Thus these frameworks are very useful
in attempts at quantifying cognitive states. We set out to analyze different
types of signals including scalp and intracerebral electroencephalography
(EEG), and magnetoencephalography (MEG) in subjects during different states of
consciousness: awake, sleep stages and epileptic seizures. The signals were
analyzed using a statistical (Permutation Entropy) and a deterministic
(Permutation Lempel Ziv Complexity) analytical method. The results are
presented in a complexity vs entropy graph, showing that the values of entropy
and complexity of the signals tend to be greatest when the subjects are in
fully alert states, falling in states with loss of awareness or consciousness.
These results are robust for all three types of recordings. We propose that
the investigation of the structure of cognition using the frameworks of
complexity will reveal mechanistic aspects of brain dynamics associated not
only with altered states of consciousness but also with normal and
pathological conditions.

##  1 Introduction

Multitude of studies focus on the investigation of patterns of correlated
activity among brain cell ensembles based on magnitudes of a variety of
synchrony indices or similar measures. A prominent common aspect that is
emerging from those studies is that of the importance of the variability in
the brain coordination dynamics. In general, neurophysiological signals
associated with normal cognition demonstrate fluctuating patterns of activity
that represent interactions among cell networks distributed in the brain [1].
This variability allows for a wide range of configurations of connections
among those net-works exchanging information, and thus it supports the
flexibility needed to process sensory inputs. Therefore, it has been argued
that a certain degree of complexity in brain signals will be associated with
healthy cognition, whereas low complexity may be a sign of pathologies [2, 3,
4]. We sought to obtain evidence for the correlation between complexity in
brain signals and conscious states, using brain electrophysiological
recordings in conscious and unconscious states.

There exist a number of statistical measures to analyses electrophysiological
recording [5]. In our work we use two well knows measures, one statistical
–Shannon entropy, a measure of unpredictability of information content in a
message [6] and the other deterministic the Lempel-Ziv complexity based in the
minimum information required to recreate the original signal [7]. For both
measures, we use use the quantifiers introduced by Bandt and Pompe [8], called
permutation vectors, these are based on the relationship of the neighbour
values belonging a time series. The Shannon entropy applying to the
permutation vector is knowing as permutation entropy (HPE) [8]. In a similar
manner the Lempel-Ziv complexity applied to the permutation vectors is called
permutation Lempel-Ziv complexity (PLZC) [9]. We used these two method to
obtain information about the signals dynamics from two differents perspective,
probabilistic (HPE) and deterministic (PLZC). The permutation entropy and the
Lempel-Ziv complexity have been employed in previous studies analyzing
electrophysiological recording in epilepsy, coma or sleep stages [10, 11, 12,
13, 14]. Moreover , there is an interesting relation, under certain
restrictions, between the Shannon entropy and the Lempel-Ziv complexity that
naturally can extend to the HPE and PLZC [15, 9].

The result we obtain are shown in a complexity-entropy graph. This kind of
representation allows us to visualize better the results. In a recent study on
chaotic maps and random sequences, it was shown that the complexity-entropy
graph allows the distinction of different dynamics which are impossible to
discern using each analysis separately ( [16], unpublished results). In our
work we analyze brain signals recorded using scalp (EEG), intracranial
electroencephalogram (iEEG) and magnetoencephalogram (MEG), in fully alert
states and in two conditions where consciousness is impaired: seizures and
sleep. The hypothesis derived from the previous consideration on the
variability of activity is that the brain tends towards larger complexity and
entropy in wakefulness as compared with the altered states of consciousness.

##  2 Methods

### Electrophysiological recordings

Recordings were analysed from 9 subjects using magnetoencephalography (MEG),
scalp electroencephalography (EEG) or intracranial EEG (iEEG). Three epilepsy
patients were studied with MEG; one epilepsy patient was studied with iEEG; 3
epilepsy patients were studied with simultaneous iEEG and scalp EEG; and 2
nonepileptic subjects were studied with scalp EEG.

For the study of seizures versus alert states, the three subjects with MEG
recordings and the one with iEEG were used. Details of the patients’
epilepsies, seizure types and the recording specifics have been presented in
previous studies (MEG patients in [17], 2005; iEEG patients in [18]). For the
study of sleep versus alert states, the 3 patients with combined iEEG and
scalp EEG have been described previously (patients 1, 3, 4 in [19]); the 2
subjects studied with scalp EEG alone had been investigated because of a
suspected history of epilepsy, but both were ultimately diagnosed with
syncope, with no evidence of epilepsy found during prolonged EEG monitoring.
In brief, the MEG seizure recordings were obtained in one patient with primary
generalized absence epilepsy, in one patient with symptomatic generalized
epilepsy, and in one patient with frontal lobe epilepsy. The iEEG seizure
recordings were obtained from a patient with medically refractory temporal
lobe epilepsy as part of the patient’s routine clinical pre-surgical
investigation.

MEG recordings were obtained using a whole head CTF MEG system (Port
Coquitlam, BC, Canada) with sensors covering the entire cerebral cortex,
whereas iEEG electrodes were positioned in various locations including the
temporal lobe epilepsy patient, the amygdala and hippocampal structures of
both temporal lobes. EEG recordings were obtained using an XLTEK EEG system
(Oakville, ON, Canada). The details of the acquisitions varied from patient to
patient (e.g., acquisition rate varied from 200 to 625 Hz) and were taken into
consideration for the data analyses. The duration of the recordings varied as
well: for the seizure study, the MEG sample epochs were of 2 minutes duration
each, with total recording times of 30-40 minutes; the iEEG patient sample was
of 55 minutes duration. The sleep study data segments were each 2-4 minutes in
duration, selected from continuous 24-hour recordings.

### Data analysis

The data were analyzed throw the permutation Lempel Ziv complexity (PLZC) and
the permutation entropy (PE). Due the relationship existing between these
quantities the result were shown in a complexity-entropy graph, to extract
information from the signals either deterministic to statistic. In this
section we give a breve explanation of both method and the relationship
between then.

#### Permutation entropy

The permutation entropy (HPE) is a measure develop by Bandt and Pompe [8], for
time series based on comparing neighboring values. The continuous time series
is mapped onto a sequence of symbols which describe the relationship between
present values and a fixed number of equidistant values at a given past time.

To understand the idea let us consider a real-valued discrete-time series
$\\{X_{t}\\}_{t\geq 0}$ , and let $d\geq 2$ and $\tau\geq 1$ be two integers.
They will be called the embedding dimension and the time delay, respectively.
From the original time series, we introduce a $d$-dimensional vector
$Y^{(d,\tau)}_{t}$:

$\mathbf{Y}_{t}^{(d,\tau)}\rightarrow(X_{t-(d-1)\tau},...,X_{t-\tau},X_{t});t\geq(d-1)\tau$

There are conditions on $d$ and $\tau$ in order that the vector
$\mathbf{Y}^{(d,\tau)}_{t}$ preserves the dynamical properties of the full
dynamical system111For EEG signals values of $d=3,...,7$ have been recommended
[8]; For the time lag, it is adequate to use a value of $\tau=1$ [20], For all
signals in this work we used the parameter $d=3,..,6$ and $\tau=1$.. The
components of the phase space trajectory $\mathbf{Y}^{(d,\tau)}$ are sorted in
ascending order. Then, we can define a permutation vector, $\Pi^{d,\tau}_{t}$,
with components given by the position of the sorted values of the component of
$\mathbf{Y}^{(d,\tau)}_{t}$ Each one of these vectors represents a pattern (or
motif). There are $d!$ possible patterns. It is possible to calculate the
frequencies of occurrence of any of the $d!$ possible permutation vectors.
From these frequencies, we can estimate the Shannon entropy associated with
the probability distributions of permutation vector. If we denote the
probability of occurrence of the $i$-th pattern by $P(\Pi^{d,\tau})_{i}=P_{i}$
with $i\leq d!$ then the (normalized) permutation entropy associated with the
time series $\\{X_{t}\\}$ is (measured in bits):

$H_{PE}=\frac{-\sum^{d!}\nolimits_{i=1}P_{i}~{}log_{2}P_{i}}{log_{2}d!}$ (1)

The fundamental assumption behind the definition of HPE is that the $d!$
possible permutation vectors might not have the same probability of
occurrence, and thus, this probability might unveil knowledge about the
underlying system.

#### Permutation Lempel-Ziv complexity

Entropy is a statistical characterization of a random variable and/or
sequence. An alternative caracterization of time series is the deterministic
notion of complexity of sequences due to Kolomogorof. In this view, complexity
is defined as the size of the minimal (deterministic) program (or algorithm)
allowing to generate the observed sequence [15, Chap. 14]. Later on, Lempel
and Ziv proposed to define such a complexity for the class of “programs” based
on recursive copy-paste operators [7].

To be more precise, let us consider a finite-size sequence
$S_{1:T}=S_{1}...S_{T}$ of size $T$, of symbols $S_{i}$ that take their values
in an alphabet $\mathcal{A}$ of finite size $\alpha=|\mathcal{A}|$. The
definition of the Lempel–Ziv complexity lies in the two fundamental concepts
of reproduction and production:

  * •

Reproduction: it consists of extending a sequence $S_{1:T}$ by a sequence
$Q_{1:N}$ via recursive copy-paste operations, which leads to
$S_{1:_{T}+N}=S_{1:T}Q_{1:N}$, i.e., where the first letter $Q_{1}$ is in
$S_{1:T}$, let us say $Q_{1}=S_{i}$, the second one is the following one in
the extended sequence of size $T+1$, i.e., $Q_{1}=S_{i+1}$ , etc.: $Q_{1:N}$
is a subsequence of $S_{1:T+N-1}$. In a sense, all of the “information” of the
extended sequence $S_{1:T+N}$ is in $S_{1:T}$.

  * •

Production: the extended sequence $S_{1:T+N}$ is now such that $S_{1:T+N-1}$
can be reproduced by $S_{1:T}$, but the last symbol of the extension can
either follow the recursive copy-paste operation (thus we face to a
reproduction) or can be “new”. Note thus that a reproduction is a production,
but the converse is false. Let us denote a production by $S_{1:T}\Rightarrow
S_{1:N+T}$.

Any sequence can be viewed as constructed through a succession of productions,
called a history $\mathcal{H}$. For instance, a history of $S_{1:T}$ can be
$\mathcal{H}(S_{1:T}):\emptyset\Rightarrow S_{1}\Rightarrow
S_{1:2}\Rightarrow\cdots\Rightarrow S_{1:T}$. The number the productions used
for the generation $C_{\mathcal{H}(S_{1:T})}$ is in this case equals to the
size of the sequence. A given sequence does not have a unique history and in
the spirit of the Kolmogorov complexity, Lempel and Ziv were interested in the
optimal history, i.e., the minimal number of production necessary to generate
the sequence. The size of the shortest history is the so-called Lempel–Ziv
complexity, denoted as
$C[S_{1:T}]=\min_{\mathcal{H}(S_{1:T})}C_{\mathcal{H}(S_{1:T})}$ [7]. In a
sense, $C[S_{1:T}]$ describes the “minimal” information needed to generate the
sequence $S_{1:T}$ by recursive copy-paste operations.

As explained above, the Lempel–Zip complexity ($C_{LZ}$) needed a alphabet of
finite size to be used. In continuos time series as EEG or MEG it is necessary
to discretized the series before calculating the $C_{LZ}$. Using the same idea
that in permutation entropy can be taken the alphabet as the set of
permutation vectors $\mathcal{A}=\\{\Pi^{(d,\tau)}\\}$ and the alphabet large
$\alpha=|d!|$. This is called permutation Lempel–Ziv complexity (PLZC)222 From
now we call the permtation Lempel–Ziv complexity as $C_{LZ}$ [9]

The most interesting thing is although analyzing a sequence from a completely
deterministic point of view, it appears that $C_{LZ}[S_{1:T}]$ sometimes also
contains the concept of information in a statistical sense. Indeed, it was
shown in references [15, 7] that for a random stationary and ergodic process,
when correctly normalized, the Lempel-Ziv complexity of the sequence tends to
the entropy rate of the process; this result were extend to the permutation
Lempel-Ziv complexity and the permutation entropy [9]; i.e.,

$\lim_{T\to+\infty}C_{LZ}[S_{1:T}]\frac{\log(T)}{T}\>=\>\lim_{T\to+\infty}\frac{H_{PE}[S_{1:T}]}{T}$
(2)

where $H_{PE}[S_{0:T-1}]$ is the joint permutation entropy of the $T$ symbols,
and the righthand side is the permutation entropy rate (entropy per symbol) of
the process. Such a property gave rise to the use of the permutation Lempel-
Ziv complexity for permutation entropy estimation purposes.

##  3 Results

The results obtained with recordings acquired during conscious state are
compared with those acquired during unconscious states, with include sleep
(all stages) and epileptic seizures. We note that while we work at the signal
level we made the reasonable assumption that the MEG and scalp EEG sensors
record cortical activity underlying those sensors and thus throughout the text
we used the term brain signals. On the other hand, the iEEG, obviously,
records signals at the source level. For all the signals the permutation
vector parameter used were $d=3,...,6$ and $\tau=1$.

###  3.1 Entropy-complexity analysis from epileptic recordings

To visualize the dynamics of entropy and complexity in the time, we use a non
overlapping running window ($\Delta=625$) corresponding to $1s$ MEG recording
points. For each window the PLZC and HPE were calculated.Fig. 1 shows the
complexity (PLZC) and entropy (HPE) values correspond to a MEG recording from
a patient suffering primary generalized epilepsy (A), secondary generalized
epilepsy (B) and an frontal lobe epilepsy (C). For patients A and B the
entropy and complexity values represented were calculated the average over the
143 channels. For patient C the values in each plot correspond a particular
channel.

One MEG channel corresponding to patient A, are shown in the inset of Fig.
1A), where the seizure is visible as a high amplitude in signal. The
complexity-entropy graph depict clearly the dynamics of the ictal event.
During conscious states (baseline) – when patients remain conscious since they
are not having generalized seizures– the PLZC and HPE tent to maximum values,
but as the patients experience seizures both values decrease widely, returning
to the original baseline values after the event.

Similar result can be seen in patient B who had 7 seizures, the seizures are
visible in the inset of Fig. 1B. We can see in the graph that the baseline and
the interictal activity – the recording between to seizures – reach always the
highest values in entropy and complexity, declining to values well below in
the ictal state (seizure). This result is repeated for each of the seizures.

In Fig. 1C we show the analysis for 4 different MEG channels corresponding to:
left frontal (LF23), left temporal (LT5), left occipital (LO41) and right
occipital (RO43). The first two belong to the region where the seizure spread.
For all channels the values of HPE and PLZC are higher in baseline, however
the entropy and complexity decay in the most affected areas (LF23, LT5), while
for the other areas (LO41, RO43) the complexity doesn’t change, there being a
small decrease in entropy. Similar result were found in the signals of the
other epileptic patients, recorded with scalp EEG and iEEG.

A possible explanation for this decreas in complexity and entropy during
seizures, is that there is higher synchrony during ictal periods (seizures),
therefore this causes the recording signals become more stereotyped, the
number of permutation vectors used to quantized the signals are smaller and
more regular giving a lower entropy and complexity. This will be further
commented in the discussion.

Figure 1:  Represent the permutation Lempel Ziv complexity ($C_{LZ}$) vs
permutation entropy ($H_{PE}$) (with parameter $d=4$ and $\tau=1$) time
tracking values for a MEG signal in epileptic patients during conscious,
baseline (BL) and unconscious, seizure (Sz) states. A) Patient with primary
generalized epilepsy, the MEG signal for one channel is plotted in the inset
(the high amplitude represent the seizure). We observe that before the seizure
entropy and the complexity the values remains very high, decreasing in the
seizure epoch and return to the original values after the seizure. B) Patient
with secondary generalized epilepsy, who had 7 seizure during the recording
period, is in the inset. When the patient stay in the inter-ictal state (no
seizure period) entropy and complexity values are higher and decreasing in the
every attack. C) Patient suffering from frontal lobe epilepsy; 4 channels were
analysed separately, left frontal (LF23), left temporal (LT5), left occipital
(LO41), right occipital (RO43). For the two recording areas affecting by the
seizure (LF23 and LT5) entropy and complexity change in the ictal state, but
for the areas which are not affected (LO41 and RO43), the $C_{LZ}$ and
$H_{PE}$ values are the same that in baseline state. The same result we
obtained for the parameter $d=3,4,5,6$ and $\tau=1$.

###  3.2 Entropy-complexity analysis during sleep stages

Te recording in these cases were of 2-4 minutes duration during wakefulness
with eyes opened (’Awake’) or closed, and in sleep stages slow-wave 2 (Sws2),
slow-wave 3-4 (Sws3-4) and rapid eye movement (’REM’). Fig. 2A shows entropy
and complexity values applying to 4 whole recording (iEEG channels): left
frontal media (LFM1), right frontal media (RFM4), left temporal anterior
(LTA1), right temporal anterior (RTA4). The various stages of sleep are
remarkably differentiated in the graph. Note how during wakefulness entropy
and the complexity is in the higher region of the graph, whereas for the slow
wave stages, the values stay in the lower region. The deepest sleep stage,
slow wave 3-4 (sws 3-4), has consistently the lowest entropy and complexity.
Interestingly, entropy during REM sleep is very close, in most cases, to the
normal, alert state. This result may not be as surprising as it appeare, if we
consider the mental activity during REM episodes that are normally associated
with dreams. The results are in agreement with those reported in [12, 13].

The results for 4 scalp EEG channels are shown in Fig. 2B, where the same
result was obtained: higher complexity and entropy for awake state and lower
for deep sleep state. In this case, during REM, the values remains between
slow-wave period and wakefulness. For the other 2 subjects analyzed we
obtained similar results. This example demonstrates that the same qualitative
result is obtained with different recording techniques. The similarity of the
result indicate that these type of analysis is not influenced by the recording
methodology.

Figure 2:  A) Each windows the channel recording from one iEEG channels
analyzing by the permutation Lempel-Ziv complexity ($C_{LZ}$) vs permutation
entropy ($H_{PE}$) graph (with parameter $d=4$ and $\tau=1$), for a patient
recording during sleep. Data samples were of 2-4 minutes duration during
wakefulness with eyes open (’Aw Oe’) , and sleep stages slow-wave 2 (‘Sws2’),
slow-wave 3-4 (‘Sws3-4’) and rapid eye movement (’REM’). The electrode
localization are: left frontal media (LFM1), rigth frontal media (RFM4), left
temporal anterior (LTA1), right temporal anterior (RTA4), the yellow circle
show the position of the channel in the brain. When the patient go in deeper
sleep states, both PLZC and HPE decreases across all channels. B) The same
analysis as in A applied to another subject (scalp EEG channel recording), as
in the previous one, awake state has the higher entropy and complexity and the
values decrease for deeper state of sleep. For the REM stage the values are
remains between the sleep stages and the wakefulness stage. The same result we
obtained for all patient analyzed with the parameter $d=3,..,6$ and $\tau=1$.

##  4 Discusion

Our results indicate a pronounced loss of entropy and complexity in brain
signals during unconscious states or in states that do not represent full
alertness (eyes closed). This is consistent with what the signals represent:
the coordinated collective activity of cell ensembles, which, in alert states,
are responsible for optimal sensory processing. This optimality requires
certain variability in the interactions among those cell networks, which will
be conceivably represented in greater complexity.

Previous work has indicated less variability in the coordinated activity
patterns in altered states of consciousness, mainly derived from the analysis
of synchronization in patients in coma [21, 22], or during seizures [17, 23].
A common feature of several theories of consciousness is the notion of a broad
distribution of cellular interactions in the brain that results in conscious
awareness (reviewed in [24]). This requirement implies that a certain, high
degree of variability in the formation and dissolution of functional cell
ensembles should take place [25], and this variability will be reflected in
higher complexity of the brain signals during alert states. Moreover in
several computational studies have revealed as well the lower complexity
associated with epilepsy and abnormal cognitive states, like schizophrenia
[26]

In fully alert states, brain recordings exhibit higher frequencies of
relatively low amplitude, and are less regular than during other states where
alertness is perturbed, including closing the eyes (when a prominent periodic
alpha rhythm appears in parieto-occipital areas, for instance). Brain cell
ensembles that need to integrate and segregate sensorimotor transformations
while they receive rich sensory-motor inputs [27]; it is then conceivable that
these characteristics will be reflected in the high entropy and complexity
values we observe. As consciousness is gradually lost, during sleep, the
values of entropy and complexity decrease because brain networks do not need
the richness in states needed to process the sensorium. The lack of arrival of
multiple sensory inputs during unconscious states decreases the need for
neurons to display many different firing frequencies, since there is not much
integration/segregation being done at those stages and there is not much
sensory load. One consequence of this change in firing patterns during
unconscious states, particularly in sleep (for a comprehensive review of the
neurophysiological mechanisms leading to slow-wave sleep and other
thalamocortical phenomena see [28]) is that the high frequencies (gamma range)
become less prominent and there is higher synchrony at lower frequencies. As
well, the amplitude of the slow waves is now high since there are more
synchronized cells. Thus, all these events result in the recording becoming
more regular and exhibiting the typical slow wave frequencies, and therefore
our complexity measures decrease as compared to alert states. These results
are consistent with measures obtained from analysis of sleep EEG using
permutation entropy [12] and other nonlinear measures, such as approximate
entropy, correlation dimension, recurrence plots and Hurst exponent, amongst
others [29, 30, 31].

In the case of the epileptic recordings we have observed that the complexity
and entropy values are larger in the interictal stage (between seizures) and
decline sharply in the ictal stage (seizures). This may be due to the fact
that seizures are characterized by excessive synchronous neuronal activity,
which generates predominance of large amplitude waveforms, the frequencies
depending on the seizure type; e.g., the frequency is low in absence seizures
(3-4 Hz), but vary substantially in temporal lobe seizures. However, the
frequencies remain relatively constant for certain time periods (originating a
distribution of periodic epochs, or laminar phases), that have been used in
the characterization of dynamical regimes in epileptiform activity [32], and
therefore the complexity and entropy tend to decrease. During the sleep stages
we also found decreased entropy and complexity as compared with alert states,
a reflection of the aforementioned emergence of highly synchronous cell
activity during slow wave sleep. On the other hand, we found that complexity
during REM sleep is similar to that of the awake state. This is conceivable
since REM episodes are normally associated with dreaming, and there is certain
cognitive activity going on in dreams, when there is partial awareness.
Previous work has shown decreases in HPE and LZC in patients under anesthesia
effects [14, 10, 33], thus the decreased complexity of brain signals in
unconscious states may be a common phenomenon.

Hence in the final analysis what we measure, at the macro(meso)scopic level
(through the recording of collective cel activity in EEG or MEG), is a
reflection of that the brain handles more information during wakefulness. A
larger code is required to manipulate more information. The complexity/entropy
of the signals used in this work have been quantified through the Bandt and
Pompe method [8], which focuses on the relative values of neighbouring data
points in a time series. Every embedding vector (or motif $\Pi^{d,\tau}_{i}$)
gives an idea of how the waveform is, in a small section, of the original
signal. As the original signal carries more variable information, the waveform
tend to be more fluctuating, and the number of distinct motifs required to map
it increases. Because of that the probability distribution of motifs
$P(\Pi^{d,\tau})$ tends to be uniform, and this caused entropy increase.
Besides, due to of the waveform fluctuation, the PLZC increases too, since
much more information is required to reconstruct the signal. In contrast, for
monotonal repetitive signals which have little new information, just a limited
number of motifs are required, e.g. for a sinusoidal signal the PLZC and HPE
tend to be zero.

We note that our present resutls are complementary to those recently obtained
using measures of coordinated activity, namely the number of configurations of
connections derived from an index of phase synchronization [1]; we should
consider that the present analysis, done on the raw signals, represent too
correlated activity as each local field potential (in case of iEEG) or signals
recorded in scalp EEG or MEG represent the collective activity in cell
ensembles, thus these signals are themselves a measure of coordinated cell
activity, and therefore it is not surprising we obtain similar observations

It can be concluded that in the awake state, when information has to be
handled is larger, the complexity and entropy of the signals recorded from the
brain tend to be higher than in absence of consciousness, a result that stems
from the distinct waveforms recorded in these mental states.

## References

  * [1] R. M. Guevara Erra, D. M. Mateos, R. Wennberg, and J. L. Perez Velazquez.  Towards a statistical mechanics of consciousness: maximization of number of connections is associated with conscious awareness.  Physical Review E, (in press), 2016.
  * [2] Douglas D Garrett, Gregory R Samanez-Larkin, Stuart WS MacDonald, Ulman Lindenberger, Anthony R McIntosh, and Cheryl L Grady.  Moment-to-moment brain signal variability: a next frontier in human brain mapping?  Neuroscience & Biobehavioral Reviews, 37(4):610–624, 2013.
  * [3] Jose L Perez Velazquez, Miguel A Cortez, O Carter Snead, and Richard Wennberg.  Dynamical regimes underlying epileptiform events: role of instabilities and bifurcations in brain activity.  Physica D: Nonlinear Phenomena, 186(3):205–220, 2003.
  * [4] D. M. Mateos, J. M. Diaz, and P. W. Lamberti.  Permutation entropy applied to the characterization of the clinical evolution of epileptic patients under pharmacologicaltreatment.  Entropy, 16(11):5668–5676, 2014.
  * [5] Katerina Hlaváčková-Schindler, Milan Paluš, Martin Vejmelka, and Joydeep Bhattacharya.  Causality detection based on information-theoretic approaches in time series analysis.  Physics Reports, 441(1):1–46, 2007.
  * [6] C. E. Shannon.  A mathematical theory of communication.  The Bell System Technical Journal, 27(4):623–656, October 1948\.
  * [7] A. Lempel and J. Ziv.  On the complexity of finite sequences.  IEEE Transactions on Information Theory, 22(1):75–81, January 1976\.
  * [8] C. Bandt and B. Pompe.  Permutation entropy: A natural complexity measure for time series.  Physical Review Letters, 88(17):174102, April 2002.
  * [9] S. Zozor, D. Mateos, and P. W. Lamberti.  Mixing Bandt–Pompe and Lempel–Ziv approaches: another way to analyze the complexity of continuous-states sequences.  The European Physical Journal B, 87(5):107, kli 2014.
  * [10] E. Olofsen, J. W. Sleigh, and A. Dahan.  Permutation entropy of the electroencephalogram: a measure of anaesthetic drug effect.  British journal of anaesthesia, 101(6):810–821, 2008.
  * [11] E. Ferlazzo, N. Mammone, V. Cianci, S. Gasparini, A. Gambardella, A. Labate, M. A. Latella, V. Sofia, M. Elia, F. C. Morabito, et al.  Permutation entropy of scalp eeg: A tool to investigate epilepsies: Suggestions from absence epilepsies.  Clinical Neurophysiology, 125(1):13–20, 2014.
  * [12] N. Nicolaou and J. Georgiou.  The use of permutation entropy to characterize sleep electroencephalograms.  Clinical EEG and Neuroscience, 42(1):24–28, 2011.
  * [13] A. G. Casali, O. Gosseries, M. Rosanova, M. Boly, S. Sarasso, K. R. Casali, S. Casarotto, M. Bruno, S. Laureys, G. Tononi, et al.  A theoretically based index of consciousness independent of sensory processing and behavior.  Science translational medicine, 5(198):198ra105–198ra105, 2013\.
  * [14] X. S. Zhang, R. J. Roy, and E. W. Jensen.  Eeg complexity as a measure of depth of anesthesia for patients.  Biomedical Engineering, IEEE Transactions on, 48(12):1424–1433, 2001.
  * [15] T. M. Cover and J. A. Thomas.  Elements of Information Theory.  John Wiley & Sons, Hoboken, New Jersey, 2nd edition, 2006.
  * [16] D. Mateos and S. Zozor.  On the analysis of signals in a permutation lempel–ziv complexity - permutation shannon entropy plane.  Manuscript submitted for publication., 2016.
  * [17] L. Garcia Dominguez, R. A Wennberg, W. Gaetz, D. Cheyne, O. C. Snead, and J. L. Perez Velazquez.  Enhanced synchrony in epileptiform activity? local versus distant phase synchronization in generalized seizures.  The Journal of neuroscience, 25(35):8077–8084, 2005.
  * [18] J. L. Perez Velazquez, L. G. Dominguez, V. Nenadovic, and R. A. Wennberg.  Experimental observation of increased fluctuations in an order parameter before epochs of extended brain synchronization.  Journal of biological physics, 37(1):141–152, 2011.
  * [19] R. Wennberg.  Intracranial cortical localization of the human k-complex.  Clinical Neurophysiology, 121(8):1176–1186, 2010.
  * [20] A. A Bruzzo, B. Gesierich, M. Santi, C. A. Tassinari, N. Birbaumer, and G. Rubboli.  Permutation entropy to detect vigilance changes and preictal states from scalp eeg in epileptic patients. a preliminary study.  Neurological sciences, 29(1):3–9, 2008.
  * [21] Vera Nenadovic, James S Hutchison, Luis Garcia Dominguez, Hiroshi Otsubo, Martin P Gray, Rohit Sharma, Jason Belkas, and Jose Luis Perez Velazquez.  Fluctuations in cortical synchronization in pediatric traumatic brain injury.  Journal of neurotrauma, 25(6):615–627, 2008.
  * [22] V Nenadovic, J. L. Perez Velazquez, and J. S. Hutchison.  Phase synchronization in electroencephalographic recordings prognosticates outcome in paediatric coma.  PloS one, 9(4):e94942, 2014.
  * [23] J. L. Perez Velazquez, L. Garcia Dominguez, and R. Wennberg.  Complex phase synchronization in epileptic seizures: evidence for a devil’s staircase.  Physical Review E, 75(1):011922, 2007.
  * [24] P. C. Klink, M. W. Self, V. A.F. Lamme, P. R Roelfsema, and S.M. Miller.  Theories and methods in the scientific study of consciousness.  The Constitution of Phenomenal Consciousness: Toward a Science and Theory, 92, 2015.
  * [25] H. Flohr.  Sensations and brain processes.  Behavioural brain research, 71(1):157–161, 1995.
  * [26] G Karl Steinke and Roberto F Galán.  Brain rhythms reveal a hierarchical network organization.  PLoS Comput Biol, 7(10):e1002207, 2011.
  * [27] G. Tononi.  An information integration theory of consciousness.  BMC neuroscience, 5(1):42, 2004.
  * [28] Alain Destexhe and Terrence J Sejnowski.  Thalamocortical assemblies: How ion channels, single neurons and large-scale networks organize sleep oscillations.  2001\.
  * [29] J. Röschke and J.B. Aldenhoff.  A nonlinear approach to brain function: deterministic chaos and sleep eeg.  Sleep: Journal of Sleep Research & Sleep Medicine, 1992.
  * [30] R. Acharya, O. Faust, N. Kannathal, T. Chua, and S. Laxminarayan.  Non-linear analysis of eeg signals at various sleep stages.  Computer methods and programs in biomedicine, 80(1):37–45, 2005\.
  * [31] N. Burioka, M. Miyata, G. Cornélissen, F. Halberg, T. Takeshima, D. T. Kaplan, H. Suyama, M. Endo, Y. Maegaki, T. Nomura, et al.  Approximate entropy in the electroencephalogram during wake and sleep.  Clinical EEG and neuroscience, 36(1):21–24, 2005.
  * [32] J.L. Perez Velazquez, H. Khosravani, A. Lozano, B. Bardakjian, P. L. Carlen, R. Wennberg, et al.  Type iii intermittency in human partial epilepsy.  European Journal of Neuroscience, 11(7):2571–2576, 1999.
  * [33] D. Li, X. Li, Z. Liang, L. J. Voss, and J. W. Sleigh.  Multiscale permutation entropy analysis of eeg recordings during sevoflurane anesthesia.  Journal of neural engineering, 7(4):046010, 2010.


# Decentralized Voltage Control with Peer-to-peer Energy Trading in a
Distribution Network

Chen Feng
School of Industrial Engineering
Purdue University
<EMAIL_ADDRESS>
&Andrew L. Lu
School of Industrial Engineering
Purdue University
<EMAIL_ADDRESS>
&Yihsu Chen
Electrical and Computer Engineering
University of California Santa Cruz
<EMAIL_ADDRESS>

###### Abstract

Utilizing distributed renewable and energy storage resources via peer-to-peer
(P2P) energy trading has long been touted as a solution to improve energy
system’s resilience and sustainability. Consumers and prosumers (those who
have energy generation resources), however, do not have expertise to engage in
repeated P2P trading, and the zero-marginal costs of renewables present
challenges in determining fair market prices. To address these issues, we
propose a multi-agent reinforcement learning (MARL) framework to help automate
consumers’ bidding and management of their solar PV and energy storage
resources, under a specific P2P clearing mechanism that utilizes the so-called
supply-demand ratio. In addition, we show how the MARL framework can integrate
physical network constraints to realize decentralized voltage control, hence
ensuring physical feasibility of the P2P energy trading and paving ways for
real-world implementations.

#### Keywords:

peer-to-peer energy trading, distributed voltage control, multi-agent
reinforcement learning

##  1 Introduction

Distributed energy resources (DERs), such as photovoltaic (PV) panels and
energy storage systems, have grown rapidly and have become a key element of a
smart energy system. Many potential benefits of DERs include increased use of
renewable energy, making the energy system cleaner and more sustainable, and
improved resilience to natural disasters. From the perspectives of consumers,
such as households and commercial buildings, when equipped with DERs, they
become prosumers; that is, they can either buy from or sell energy to the
power grid. When compensated properly, consumers will have more incentives to
invest in DERs, which, in turn, will accelerate the many benefits of DERs.

From our perspective, however, there are at least two major challenges to the
broad adoption of DERs. How to properly compensate DER owners (aka prosumers)
represents one such challenge; the other is how to efficiently utilize the
distributed resources to maintain system reliability and feasibility, while
not incurring inconvenience or intruding privacy. To address the first
challenge, feed-in-tariff (FIT) has been a widely adopted approach to entice
DER investment, which is a non-market-based approach that just pays DER owners
a fixed rate for the energy they produce. However, government agencies or
utility companies usually do not know what a fair FIT rate is, and the tariffs
are often set too high. As an alternative, peer-to-peer (P2P) market
mechanisms have been proposed to integrate and incentivize DERs. Such
mechanisms include the bilateral matching approach [Liu et al., 2019], double
auctions [Zhao et al., 2022], and the supply-demand-ratio mechanism [Liu et
al., 2017b], among others. Extensive reviews of such work can be found in
[Tushar et al., 2020, Guerrero et al., 2020].

As characterized in [Tushar et al., 2020], a P2P energy market consists of two
layers: the virtual layer (aka the financial layer) and the physical layer.
Many works in the literature focus on one layer only. For example, both [Zhao
et al., 2022] and [Qiu et al., 2021] study decentralized agents in repeated
double-auctions in a P2P market, with the former using a multiagent, multi-
armed bandit approach, while the later using a multi-agent reinforcement
learning (MARL) framework. A similar MARL approach is applied in [Ye et al.,
2021] with the focus on scalability and privacy-preserving of the algorithm.
In terms of specific trading and pricing mechanisms in a P2P market, bilateral
contract approaches are studied in [Morstyn et al., 2018] and [Sorin et al.,
2018]. The former employs a trading mechanism from the economics literature
[Hatfield et al., 2013]; while the later uses a consensus-based algorithm.
Other market mechanisms, including bill sharing, mid-market rate and supply-
demand ratio (the same mechanism used in this paper), are compared in [Zhou et
al., 2018] using a multi-agent simulation approach. None of the mentioned
works consider physical distribution network constraints or voltage control.
Such physical-layer considerations are the main focus of [Liu et al., 2017a]
and [Liu and Wu, 2021], among others. The former paper uses alternating
direction method of multipliers (ADMM) to realize distributed control; while
the later one employs a MARL approach.

Much fewer workers consider both trading mechanisms and physical constraints
in a P2P energy market. Some notable works along this direction include the
following. [Guerrero et al., 2018] employs a matching mechanism similar to
stock trading to facilitate bilateral trading among prosumers, which are
modeled as zero-intelligence-plus agents. Physical constraints, including
voltage, distribution lines and energy losses, are considered through the so-
called sensitivity factors; that is, each bilateral transaction is analyzed
through their impacts on the network, and infeasible transactions are
rejected. In [Lu et al., 2019], network capacity constraints are considered,
but the pricing of DER-generated energy is somewhat arbitrary. Two most
relevant papers to our work are [Samende et al., 2022] and [Biagioni et al.,
2021]. The former paper uses the supply-demand ratio approach for market
clearing. Distribution network constraints are represented by an ACOPF
formulation. Prosumers bidding policies, along with the network constraints,
are trained through a centralized training, decentralized implementation
framework, the so-called multi-agent deep deterministic policy gradient
(MADDPG) algorithm. Such an algorithm requires the pooling of all agents’
state and action spaces for training, and can be a major hurdle for real-world
implementation. The later paper, [Biagioni et al., 2021], uses a completely
decentralized training-learning framework and imposes system voltage violation
on each agent to realize decentralized voltage control. However, their work
does not have a specific P2P trading mechanism.

In this work, we propose a framework that combines the strengths of the two
papers mentioned above. On the algorithm side, we adopt a completely
decentralized MARL approach that does not require centralized learning, and
integrate both the virtual and physical layers into the algorithmic framework.
At the physical layer, we explicitly include smart inverters’ capability to
generate reactive power to aid in decentralized Volt-VAR control. At the
virtual layer, we represent each prosumer’s bidding problem as a Markov
decision process (MDP) and employ the supply-demand ratio mechanism to clear
the P2P market. As a result, complete decentralization is realized in both the
virtual and physical layers. While having a centralized and independent
distribution system operator (DSO) to operate a P2P energy market has been
suggested for a long time, there are many issues facing a potential DSO
[Baker, 2020], and there are no functional DSOs in the US market. Hence, we
aspire that our fully decentralized framework will pave the way for easier
participation and better utilization of the vast DER resources.

The rest of the paper is organized as follows. The P2P market and individual
prosumer’s decision-making problem are presented in Section 2. The multi-agent
learning framework and the corresponding algorithm are described in Section 3.
The numerical simulations and results are presented and analyzed in Section 4,
followed by the conclusion and future work in Section 5.

##  2 Problem Formulation

In this section, we first describe a clearing mechanism, referred to as the
supply-demand-ratio (SDR) mechanism, that can address the two potential issues
faced by a renewable energy-dominated P2P market: (i) breakdown of marginal-
cost-based pricing due to the all-zero-marginal-cost resources, and (ii) the
potential bang-bang outcomes in a double-auction-based mechanism as the
reserve prices for buyers (the UR) and the sellers (the FIT) are publicly
known.

Based on the SDR clearing mechanism, we present the formulation of each agent
(aka, a prosumer)’s problem as a Markov decision process, which serves as a
building block for the MARL framework.

###  2.1 Market Clearing Mechanism

Consider a P2P energy market in a distribution network in which consumers and
prosumers (all referred to as agents hereinafter) can buy or sell energy. Let
$\mathcal{I}=\\{1,2,...,I\\}$ be the set of all agents in the distribution
network. Market clearing is organized in fixed time steps, such as every hour
or every 15 minutes, and trade is generally carried out at some time ahead,
such as day-ahead or hour-ahead. At each time step $t$, every agent
$i\in\mathcal{I}$ can submit bids, denoted by $b_{i,t}$, to buy or sell
energy. If $b_{i,t}$ is positive, it represents a bid to sell; while a
negative quantity represents a bid to buy. In addition, we use
$\mathcal{B}_{t}$ and $\mathcal{S}_{t}$ to denote the set of buyers and
sellers at a time step $t$, respectively. By the notation above, we have that
$\mathcal{B}_{t}=\\{i:b_{i,t}<0\\}$ and $\mathcal{S}_{t}=\\{i:b_{i,t}\geq
0\\}$. Note that the set of buyers and sellers can certainly change over time,
reflecting real-world situations where prosumers may buy or sell energy in a
particular time period, depending on their own demand and energy generation.

To clear the market at time $t$, we adopt the supply-demand ratio (SDR)
mechanism, as proposed in [Liu et al., 2017b]. The mechanism is
straightforward to describe and to implement. First, the SDR at $t$ is defined
as the ratio between total sell and total buy bids; that is,

$\displaystyle
SDR_{t}:=\frac{\displaystyle\sum_{i\in\mathcal{S}_{t}}b_{i,t}}{\displaystyle-\sum_{i\in\mathcal{B}_{t}}b_{i,t}}.$
(1)

To ensure that $SDR_{t}$ is well-defined, we assume that
$\mathcal{B}_{t}\neq\emptyset$ for any $t$. This is a reasonable assumption as
it is highly unlikely in a distribution network that no customer needs energy.
(In the extreme case where every customer in the market is a prosumer and they
all have excess energy to sell at a time $t$, then the P2P market can simply
be suspended and all excess energy sold to a utility or an DER aggregator at a
pre-defined price (or the real-time wholesale market price). This will not
affect our modeling and algorithm framework in any way.) On the other hand,
$\mathcal{S}_{t}$ can be $\emptyset$ at a given $t$ (meaning that no prosumer
has excess energy to sell in $t$). In this case, $SDR_{t}=0$.

With the definition of $SDR_{t}$, in any particular round $t$, if $SDR_{t}\geq
1$, which means there are more sell bids than buy bids, all excess energy in
the P2P market is assumed to be sold to a utility company, a DSO or an
aggregator at a pre-defined rate, which we just use $FIT$ to denote in
general.111As mentioned in Section 1, the $FIT$ rate is an artificial rate set
by utilities/policy makers. For all intents and purposes, it can be set as
zero and will not affect the pricing mechanism or our algorithmic framework in
any way. If $0\leq SDR_{t}<1$, which means that the sell bids are less than
buy bids in round $t$, then the unmet demand bids will purchase energy at a
pre-defined utility rate, denoted as $UR$. Without loss of generality, we
assume that $FIT<UR$.222If $FIT>UR$, it means that energy consumers pay more
than the utility rate to purchase energy from the prosumers, which is
equivalent to direct subsidy from energy consumers (including low-income
consumers who do not have the means to install solar panels) to prosumers.
This cannot be justified from equity perspective.

Based on the above rules, the market price in a SDR mechanism can be defined
as follows:

$\displaystyle P_{t}:=P(SDR_{t})$ $\displaystyle:=\begin{cases}(FIT-UR)\cdot
SDR_{t}+UR,&0\leq SDR_{t}\leq 1\\\ FIT,&SDR_{t}>1.\end{cases}$ (2)

$P_{t}$ is the market clearing price in round $t$, which is a piece-wise
linear function that takes the value of $SDR_{t}$ as input. as illustrated in
Figure 1.

Figure 1: SDR market price function

The interpretation of the $SDR$-based market clearing price is
straightforward. When the total supply bids are less than the total demand
bids (that is, $0\leq SDR_{t}<1$), the market price is simply a downward-
sloped linear function connecting the de facto price ceiling ($UR$) and price
floor ($FIT$). When the total supply bids exceed total demand bids (that is,
$SDR_{t}>1$), the market price is simply the price floor $FIT$. There may be
critisms of the $SDR$-approach as the pricing mechanism may seem to be
arbitrary (for example, instead of having a linear function, the pricing
function could be replaced by any downward-shaped nonlinear functions, which
would lead to different market prices with the same $SDR$) and is not based on
any economic theories. However, from our perspectives, the simplicity (and
transparency) of the SDR mechansim is well-suited for a P2P energy market with
all-zero marginal cost resources. In [Zhao et al., 2022], it is identified
that a local P2P energy market using traditional market clearing mechanisms
may not work well due to the special characteristics of such a market. Since
uncleared uncleared electricity demand bids have to be met through buying from
a utility company at the rate of UR, and uncleared energy supply bids have to
be sold to a utility at FIT. Both the UR and the FIT become the de facto price
ceiling and price floor (or reservation prices) for all bids, respectively.
However, both the UR and the FIT are public information, coupled with the fact
that all the supply sources are of zero-marginal costs, making the P2P market
outcomes likely of the bang-bang type; that is, when the total demand bids are
more than supply bids, the market price would be UR; reversely, the price
would be FIT. We view the SDR approach as a fair and reasonable approach to
avoid such bang-bang outcomes, especially for a new market in which the supply
resources are likely much less than demand.

Based on the above descriptions, we can also easily write out a prosumer’s
cost of energy purchase (as a buyer) or profit of energy sales (as a seller)
from the P2P market in $t$:

$\displaystyle
r^{m}_{i,t}:=\begin{cases}\mathbb{I}_{i\in\mathcal{B}_{t}}\times\Big{[}SDR_{t}\cdot
P_{t}\cdot b_{i,t}\\\ +\ (1-SDR_{t})\cdot UR\cdot b_{i,t}\Big{]}\\\
+\mathbb{I}_{i\in\mathcal{S}_{t}}\times\Big{(}P_{t}\cdot b_{i,t}\Big{)},&0\leq
SDR_{t}\leq 1\\\ FIT\cdot b_{i,t},&SDR_{t}>1,\end{cases}$ (3)

where $\mathbb{I}_{i\in\mathcal{B}_{t}}$ is an indicator function such that
$\mathbb{I}_{i\in\mathcal{B}_{t}}=1$ if $i\in\mathcal{B}_{t}$, and 0
otherwise. Similarly, $\mathbb{I}_{i\in\mathcal{S}_{t}}=1$ if
$\in\mathcal{S}_{t}$ and 0 otherwise.

###  2.2 Markov Decision Process for a Prosumer

In a local P2P energy market, we assume that each prosumer (such as a
household or a commercial building) locates at one of the buses of the
distribution network, and has three components under their control: an HVAC
system, a photovoltaic (PV) system (either on rooftop or ground-mounted) with
a smart inverter, and an energy storage system. Each prosumer also has fixed
base load of energy consumption. In a particular time period $t$, whether a
prosumer is a net energy seller or a buyer depends on their energy consumption
from base load and HVAC versus the level of PV generation and energy storage.

Due to the presence of energy storage, each agent’s bidding and charging
decisions are linked over time, which naturally calls for a dynamic
programming framework to model the decision process. As such, we formulate an
agent’s energy trading and control process as a discrete-time MDP, and the key
modeling components are described in the following.

$\bullet$

Observation: agent $i$ at time step $t$ can observe
$s_{i,t}:=(h,temp^{o}_{t},p_{i,t},d^{p}_{i,t},d^{q}_{i,t},v_{i,t},e_{i,t},P_{t-24})$
from an observation space $\mathcal{S}_{i}$, where $h$ is the time of a day
(such as a specific hour in a day), $temp^{o}_{t}$ is the outside temperature,
$p_{i,t}$ is the PV (real) power generation, $d^{p}_{i,t}$ and $d^{q}_{i,t}$
are the inflexible demand of active and reactive power of agent $i$,
respectively, $v_{i,t}$ is the magnitude of the voltage, $e_{i,t}$ is the
state of charge of agent $i$’s energy storage system, and $P_{t-24}$ is the
past P2P energy market price at the same time on the most recent day (or a
random price between $FIT$ and $UR$ if $t-24\leq 0$; namely, the first day)
observed at time step $t$. We assume that each agent controls only one PV
system and one energy storage system. Then each element of the observation is
a scalar, and the dimension of the observation space is 8; that is,
$\mathcal{S}_{i}\subset\Re^{7}$. Note that among the observed states, only
$e_{i,t}$ is the state variable affected by agent $i$’s action and state at
the previous time step, whose transition function is given below. The other
observations, $temp^{o}_{t}$ and $p_{i,t}$, are exogenous to agent $i$. It is
reasonable to assume that the outside temperature $temp^{o}_{t}$ is bounded by
the historical highest and lowest temperature, and PV generation $p_{i,t}$ is
subject to the PV capacity. Therefore, we can assume that the observation
space is bounded. In our simulation, we use simple random samples to create
variations of the temperature and PV generation for each agent. The details
are given in Section 4.1.

$\bullet$

Action/Control variables: agent $i$ in the time step $t$ takes action
$a_{i,t}:=(a^{flow}_{i,t},a^{temp}_{i,t},a^{q}_{i,t},a^{s}_{i,t},b_{i,t})$
from action space $\mathcal{A}_{i}$, with $a^{flow}_{i,t},a^{temp}_{i,t}$
being $Z_{i}$-dimension vectors where $Z_{i}$ is the total number of zones in
agent $i$’s HVAC system, and the other elements being scalars. Hence, the
dimension of the action space is $2Z_{i}+3$; that is,
$\mathcal{A}_{i}\subset\Re^{2Z_{i}+3}$. The variable $a^{flow}_{i,t}$ is the
vector of HVAC flow rate set for each zone of agent $i$’s HVAC system,
$a^{temp}_{i,t}$ is the vector of HVAC discharge temperature for each zone of
agent $i$’s HVAC system, $a^{q}_{i,t}$ is the reactive power generation by
agent $i$’s smart inverter, $a^{s}_{i,t}$ is energy charged (if
$a^{s}_{i,t}>0$) or discharged (if $a^{s}_{i,t}<0$) of the battery, and
$b_{i,t}$ is agent $i$’s energy bid to the P2P market in time $t$ (either buy
or sell). The reactive power generation and energy charge or discharge are
subject to PV and battery capacity, respectively. It is also reasonable to
assume that there are lower and upper bounds for the control of HVAC flow rate
and discharge temperature. Therefore, we can assume that the action space
$\mathcal{A}_{i}$ is bounded.

$\bullet$

State transition: The state transition function for agent $i$ is:

$\displaystyle\begin{split}e_{i,t+1}=&\max\Big{\\{}\min\Big{(}e_{i,t}+\eta^{c}_{i}\max(a^{s}_{i,t},0)\\\
&+\frac{1}{\eta^{d}_{i}}\min(a^{s}_{i,t},0),\overline{e}_{i}\Big{)},0\Big{\\}},\end{split}\vspace*{-5pt}$
(4)

where $\eta^{c}_{i}$ is the charging efficiency for agent $i$ defined by the
ratio of energy charged to the battery divided by the energy consumed in
charging, $\eta^{d}_{i}$ is the discharging efficiency for agent $i$ defined
by the ratio of energy output from the battery divided by the energy consumed
in discharging, and $\overline{e}_{i}$ is the battery capacity of agent $i$.

$\bullet$

Reward: Each agent $i$’s reward function has three components as follows:

$\displaystyle r_{i,t}:=r^{c}_{i,t}+r^{m}_{i,t}+r^{v}_{t}/I.$ (5)

The first term, $r^{c}_{i,t}$, is the reward for the thermal comfort of agent
$i$’ from HVAC system, which follows the reward model in the NREL project
[Biagioni et al., 2021]:

$\displaystyle\begin{split}r^{c}_{i,t}:=&-\sum_{j=1}^{Z_{i}}\Big{[}\max(0,temp^{j}_{i,t}-\overline{temp})^{2}\\\
&+\max(0,\underline{temp}-temp^{j}_{i,t})^{2}\Big{]},\end{split}$ (6)

where $temp^{j}_{i,t}$ is the temperature of zone $j$ of agent $i$ at time
step $t$, which is determined by the outside temperature, the flow rate of
zone $j$, and the discharge temperature of zone $j$ of agent $i$’s HVAC
system.333Details of the HVAC system model can be found at
https://github.com/NREL/PowerGridworld The two parameters, $\overline{temp}$
and $\underline{temp}$, represent the upper and lower limits of the comfort
temperature range, respectively.

The second term $r^{m}_{i,t}$ is the energy purchase cost or sales profit of
agent $i$ from the P2P energy market, as described in Equation (3) in Section
2.1.

The third term, $r^{v}_{t}$, is the system voltage violation penalty defined
as in [Biagioni et al., 2021]:

$\displaystyle\begin{split}r^{v}_{t}:=&-\lambda\sum_{j:Bus}\Big{[}\max(0,v^{j}_{t+1}-\bar{v})\\\
&+\max(0,\underline{v}-v^{j}_{t+1})\Big{]},\end{split}$ (7)

where $v^{j}_{t+1}$ is the voltage magnitude of Bus $j$ after each agent makes
the decision, and $\bar{v},\underline{v}$ is the upper and lower limit of the
bus voltage. $I$ is the total number of agents in the system, and $\lambda$ is
a sufficiently large number. We assume that the violation penalty is shared
equally by all agents in the system. To calculate $v^{j}_{t}$ at each bus,
after all supply and demand bids have been cleared in $t$, we solve the
standard bus-injection model as follows:

$\displaystyle\begin{split}p_{k}=&\sum_{j=1}^{N}|V_{k}||V_{j}|(G_{kj}\cos(\theta_{k}-\theta_{j})\\\
&+B_{kj}\sin(\theta_{k}-\theta_{j})),\\\
q_{k}=&\sum_{j=1}^{N}|V_{k}||V_{j}|(G_{kj}\sin(\theta_{k}-\theta_{j})\\\
&-B_{kj}\cos(\theta_{k}-\theta_{j}))\\\ &\text{ for
}k=1,2,\dots,N,\end{split}$ (8)

where $N$ is the number of buses in the network, $j$ and $k$ are bus
subscripts, $p_{k}$ is net real power injection at bus $k$, $q_{k}$ is the net
reactive power flow at bus $k$, $G_{kj}$ and $B_{kj}$ are the real and
imaginary parts of the admittance of branch $kj$ respectively, $V_{k}$ is the
voltage of bus $k$, and $\theta_{k}$ is the phase angle at bus $k$.

After market clearing and solving for voltage magnitudes, the total reward for
the agent $i$ at time step $t$ is then realized. Note that while we do not
explicitly assume to have a DSO in the distribution network, there still has
to be an entity to solve Equation (8). A utility company can certainly assume
such a role. The solution process may even be carried out by a distributed
ledger system, such as on a blockchain. We want to emphasize that since the
market clearing happens before real-time delivery, the voltage regulation
penalty is not real financial penalty; it is only some information to be sent
to each agent’s reinforcement learning algorithm for training and learning
purposes.

$\bullet$

Objective function: At each time step $t$, agent $i$ observes the environment
(temperature and PV generation) and its state of charge in the battery
$e_{i,t}$, decides on an action $a_{i,t}\in\mathcal{A}_{i}$, and then receives
a reward $r_{i,t}$. A policy $\pi_{i}$ in this case is a conditional joint
probability density function of actions in $\mathcal{A}_{i}$, given an
observation in $\mathcal{S}_{i}$; that is, agent $i$’s choice of $a_{i,t}$ at
time $t$ given an observation $s_{i,t}\in\mathcal{S}_{i}$ follows the joint
density function $\pi_{i}(\cdot|s_{i,t}):\mathcal{A}_{i}\to[0,\infty)$. The
expected long-term payoff of agent $i$ under policy $\pi_{i}$ is defined as:

$\displaystyle
J_{i}(\pi_{i})=\mathbb{E}_{\pi_{i}}[\sum_{t=0}^{\infty}\gamma_{i}^{t}r_{i,t}],$
(9)

where $\gamma_{i}\in(0,1)$ is agent i’s discount factor.

Note that the temperature is assumed to be finite, the voltage is based on a
per-unit system and hence bounded, and the reward function on the action and
observation space is bounded. Therefore, $J_{i}$ is finite and there exists an
optimal policy $\pi_{i}^{*}$ such that
$\pi_{i}^{*}=\sup_{\pi_{i}}J_{i}(\pi_{i})$. Each agent employs certain
algorithms to try to find such an optimal policy.

##  3 Multi-agent Reinforcement Learning

While the previous section describes the individual agent’s decision-making
problem as an MDP, the agents’ rewards (and long-term payoff) depend on the
collective actions of all the agents. Before introducing the multi-agent
framework, we first discuss algorithms for a single agent to solve their MDP
problem.

As described in Section 2.2, each agent has only partial observation of the
system (that is, an agent does not observe other agents’ state variables). In
addition, each agent’s reward is determined by a complex dynamic HVAC model
which may not be known in real life, and it depends on the solution to a non-
linear power flow equation system which is hard to solve for large scale.
Therefore, it is unrealistic to figure out the value function and optimal
policy from the state-transition and reward function in our case by using a
model-based reinforcement learning (RL) approach, including classic MDP or
general dynamic programming algorithms. In contrast, model-free algorithms,
such as actor-critic reinforcement learning algorithms, do not need the
information about state-transition and reward function when solving the MDP
problem. In general, such algorithms learn the value function and optimal
policy by trying different policies many times and measuring their rewards. It
has been shown to perform well when dealing with problems where state-
transition and reward function are complex or inaccessible [Konda and
Tsitsiklis, 1999]. Therefore, the actor-critic reinforcement learning
algorithm can be a good approach to our problem. To start with, first define
the state-action value function, or Q-function under policy $\pi$ for agent
$i$, $Q^{\pi}(s_{i},a_{i})$ as the expected long-term payoff of taking action
of $a_{i}$ in state $s_{i}$, and following policy $\pi$ afterward:

$\displaystyle
Q_{i}^{\pi}(s_{i},a_{i})=\mathbb{E}_{\pi}\Big{[}\sum_{t=0}^{\infty}\gamma^{t}r_{i,t}|s_{i,0}=s_{i},a_{i,0}=a_{i}\Big{]}.$
(10)

Then define the state value function of a state $s$ under policy $\pi_{i}$ for
agent $i$ as:

$\displaystyle
V_{i}^{\pi_{i}}(s_{i})=\mathbb{E}_{\pi}[\sum_{t=0}^{\infty}\gamma^{t}r_{i,t}|s_{i,0}=s_{i}].$
(11)

In an actor-critic algorithm, for agent $i$, “actor” is a function with
parameter $\theta_{i}$:
$\pi_{\theta_{i}}:\mathcal{S}_{i}\times\mathcal{A}_{i}\to[0,\infty)$ that
characterizes the agent’s policy $\pi_{i}$, and “critic” is a function with
parameter $\omega_{i}$: $V_{\omega_{i}}:\mathcal{S}_{i}\to\mathcal{R}$ that
approximates agent $i$’s state value function under policy $\pi_{i}$. In most
of the state-of-the-art actor-critic algorithms such as proximal policy
optimization (PPO)[Schulman et al., 2017], both the “actor” and “critic” are
represented by neural networks. Algorithm 1 shows the pseudocode for the
training of actor and critic networks using PPO approach [Schulman et al.,
2017].

In this work, the traditional single-agent RL framework is extended to a
decentralized multi-agent reinforcement learning framework, in which agent $i$
still learns an individual actor network $\pi_{\theta_{i}}$ and critic network
$V_{\omega_{i}}^{\pi_{\theta_{i}}}$ based on the agent’s own observations and
actions. This is illustrated in Figure 2.

Figure 2: Architecture of the MARL framework

Compared to MARL algorithms for centralized training and decentralized
execution, such as the MADDPG algorithm [Lowe et al., 2017], our fully
decentralized approach has the advantage of less computational burden and
greater protection of agents’ privacy. The drawback of fully decentralized
MARL is that it may have poor performance due to the non-stationary
environment from the perspective of each single agent. However, in our case,
agents are coupled through the P2P market and system voltage in their reward
functions. Each agent is able to partially observe the market and voltage
information via the historical prices and system voltage violation penalty;
hence the impact of non-stationarity of the state variables on each agent’s
decision-making can be reduced. Therefore, the fully decentralized MARL,
rather than centralized training and decentralized execution MARL approach, is
chosen for our case. Specifically, each agent trains an individual PPO policy
separately, as in Algorithm 1 (the agent index $i$ is omitted).

Initialize training dataset $\mathcal{D}$, critic and actor parameters
$\theta$ and $\omega$ ($\theta$ and $\omega$ are initialized by Glorot uniform
initializer as in [Glorot and Bengio, 2010]);

for  _each time step t = 0,1,…,T_ do

$\text{1: observe}\ s_{t}$ and adopt
$a_{,t}\sim\pi_{\theta_{i}}(\cdot|s_{t});$

2: The power flow equation is solved and $r_{t}$ is computed according to Eq.
(7) ; $(s_{t},a_{t},r_{t},s_{t+1})$ is stored in $\mathcal{D}$;

3: Update policy parameters $\theta$ by solving the problem:

$\max_{\theta}\frac{1}{|\mathcal{D}|}\sum_{\tau\in\mathcal{D}}\min\Big{(}M,g(\epsilon,A^{\pi_{\theta}}(s^{\tau},a^{\tau}))\Big{)}$
via stochastic gradient ascent, where  $M=\frac{\pi_{\theta}(a^{\tau}\mid
s^{\tau})}{\pi_{\theta}(a^{\tau}\mid
s^{\tau})}A^{\pi_{\theta}}(s^{\tau},a^{\tau}),$
$g(\epsilon,A)=\begin{cases}(1+\epsilon)A&A\geq 0\\\
(1-\epsilon)A&A<0\end{cases},$ $\epsilon$ is a hyperparameter, say,
$\epsilon=0.2$, and $A^{\pi}(s,a):=Q^{\pi}(s,a)-V^{\pi}(s)$;

4: Fit value function parameters by regression using mean-squared error:
$\omega^{*}=\arg\min_{\omega}\frac{1}{|\mathcal{D}|}\sum_{\tau\in\mathcal{D}}(V_{\omega}(s^{\tau})-r^{\tau})^{2}$;

end for

Algorithm 1 PPO algorithm with distributed voltage control in a P2P market for
each agent

##  4 Simulation and Numerical Results

###  4.1 Test Case

We test our problem using the IEEE 13-bus test case. The topology of the test
network is shown in Figure 3.

Figure 3: IEEE-13 test feeder

Detail information on the test feeder, such as the base load of each bus and
line data, can be found in [Saha and Johnson, 2016]. In our case, each bus has
one prosumer, except for the substation bus; hence, there are 12 prosumers
(aka agents) in the test case.

The time step is set to be one hour. The baseload ($d_{i,t}$) and PV
generation ($p_{i,t}$) of each prosumer, and the outside temperature
($temp^{o}_{t}$) are all random variables. We assume that the expected values
of these random variables follow fixed daily shapes, which are the same for
all prosumers. More specifically, the shape of the average daily temperature
is shown in Figure 4.

Figure 4: Average outside temperature over one day

The average daily base load and PV generation shapes are represented in the
form of percentage to their peak values, instead of absolute values, as shown
in Figure 6 and Figure 5, respectively.

Figure 5: Average base load shape over one day (the shape is the same for all
buses). The load coefficient of a bus in hour $t$ is the ratio of hour $t$’s
base load to the base load of the corresponding bus in IEEE 13
(https://github.com/tshort/OpenDSS/tree/master/Distrib/IEEETestCases/13Bus).
Figure 6: Average PV generation shape over 1 day. The generation coefficient
of Hour $h$ is the ratio of PV real power generation of Hour $h$ to the PV
real power capacity, which is assumed to be $30kW$ for each agent.

To generate specific observations by agent $i$ at time $t$, a random sample
based on a uniform distribution is drawn at each $t$, with the mean being the
values indicated in Figures 4, 5 and 6, and the range is 0.95 to 1.05 times of
the mean values.

The energy storage setting is assumed to be the same for all prosumers, with
the storage capacity set at $50kWh$, and the charging and discharging
efficiency being 0.95 and 0.9, respectively. The output of the smart inverter
and PV panel is contrained by a apparent power of 50 VA. The initial state of
storage can be arbitrary; we just set it to zero for each agent at the very
beginning. Assume that each agent has the same 5-zone HVAC system. The
distribution network’s voltage limit is set to $[0.96\ pu,1.04\ pu]$. The
lower and upper bound of the comfort temperature range for each prosumer is
set at 22∘C and 28∘C, respectively. $\lambda$ is set to $10^{4}$ and $\gamma$
is set to 0.999 for all agents. $UR$ and $FIT$ are set at 14 cents / kWh and 5
cents / kWh for all time, respectively.

Policy training is implemented in RLLib, 444Details of RLLib can be found at
https://docs.ray.io/en/latest/rllib/index.html which is a platform for
reinforcement learning algorithm training built on Python. It allows users to
create custom “environment” that can describe their MDP problems on the
platform. Furthermore, in the custom “environment”, users can train existing
RL algorithms in RLLib or self-developed algorithms in TensorFlow or PyTorch.
To create an environment for our own problem, we add a P2P market and VVC
features to the power system environment built by NREL [Biagioni et al.,
2021]555Details of the package can be found at
https://github.com/NREL/PowerGridworld on RLLib. The power flow equations are
solved by the open distribution system simulator (OpenDSS) in
Python.666Details of OpenDSS package in Python can be found at
https://pypi.org/project/OpenDSSDirect.py/ The algorithm used to train the
policy is the PPO built in RLLib.

###  4.2 Result Analysis

Figure 7: 30-episode moving average of episodic total cost Figure 8: Hourly
prices over the last three days Figure 9: System voltage deviation. The total
system voltage deviation (p.u.) in the $n$-th episode is calculated as
$\sum_{t=24(n-1)}^{24n-1}\sum_{j:Bus}[\max(0,v^{j}_{t}-\bar{v})+\max(0,\underline{v}-v^{j}_{t})]$.

To present the numerical results, we first introduce the term episodic total
cost. The episodic total cost of the $n$-th episode is defined as the total
discounted cost of all agents over the $n$-th episode (24 hours):
$\sum_{t=24(n-1)}^{24n-1}\sum_{i=1}^{12}\gamma^{t}r_{i,t}$. The convergence
curves of the 30-episode moving average of the episodic total cost in the
distribution network with and without a P2P energy market are shown in Figure
7. In the case without P2P markets, agents can only buy energy from the grid
at the utility rate and sell energy to the grid at FIT. Agents still use the
distributed PPO algorithm to determine their buying/selling and charging
strategies. For the purpose of comparison, we also present the results for the
case without a P2P market. In such a case, the sellers receive payments at
$FIT$ and the buyers pay the utility rate $UR$. It can be seen that the mean
episodic total costs gradually converge to a low level at the end of the
training in both cases, which shows that the MARL decentralized training
approach performs well in terms of convergence in our case. Furthermore,
agents in the distribution network with a P2P energy market reach a lower
stationary cost level than those without a P2P energy market. This is expected
as the introduction of the P2P market allows agents to purchase energy at a
price lower than $UR$ and sell at a price higher than $FIT$.

Figure 8 presents the market price of the P2P market in chronological order
over the last three days after thousands of training episodes. The market
price remains at $UR$ one third of the day and drops from it only around noon
on each day. This is reasonable because the peak hours of PV generation are
around 12 pm, as shown in Figure 6. The market price remains at $UR$ during
the rest of the day, as energy storage and PV generation cannot meet the
demand for each prosumer and hence, no supply bids into the P2P market.

Figure 9 shows the gradual reduction of the system voltage deviation after
training. It can be seen that the voltage violation of the system is high at
the beginning. However, it drops quickly and converges to almost zero after
half of the training process. On the one hand, the results are remarkable in
the sense that voltage violation reduction is realized through purely
decentralized learning. As specified in Eq. (7), each agent will receive a
penalty if the collective bids submitted by all agents lead to a voltage
violation. The agents then simply try to submit a different bid. This is
entirely a trial-and-error process (with the underlying assumption that a
feasible collection of bids exists, and hence can be eventually reached). This
process is possible since a bid round always takes place some time before the
physical delivery of energy. On the other hand, this trial-and-error approach
to reduce voltage violation may be viewed as inefficient, and it is hard to
establish any theoretical results to guarantee the convergence (to a steady
state) of the decentralized learning process. An immediate next step of this
research is to investigate if there is any (privacy conserving) communication
mechanism among agents that would help accelerate the voltage violation
reduction, as well as leading to provable convergence of the overall MARL
framework. Two possible general approaches along this line have been proposed
in [Zhang et al., 2018, Qu et al., 2019]. Nevertheless, we hope that our
current work can serve as a benchmark to be compared with other decentralized
MARL algorithms that integrate market clearing and physical constraints.

##  5 Conclusion and Future Work

In this work, we introduced a decentralized framework to integrate DERs,
including energy storage, in a market-based approach while maintaining the
feasibility of a distribution system. Within such a framework, each agent’s
decision problem in a repeated P2P market can be automated through an RL-based
algorithm, as the learning and training of the algorithms are completely
decentralized. In addition, the privacy of the agents can be preserved, as no
exchange of private information is needed. The initial numerical simulation
shows encouraging results, as this framework is able to achieve lower energy
costs for all prosumers and realize decentralized voltage regulation for the
distribution network upon convergence.

Future work will focus on two directions. First, as mentioned earlier, whether
there exists a provably convergent MARL algorithm that is applicable to a
distribution-level energy market is an open question. We are actively
investigating in this direction. Second, the scalability of our framework (or
any MARL-based approach) to the size of real-world systems still needs to be
tested. One potential approach to improve the performance and scalability is
to study the possibility of increasing communication between neighbors without
disclosure of privacy to improve the performance or our MARL framework.

## References

  * [Baker, 2020] Baker, P. (2020).  Challenges facing distribution system operators in a decarbonised power system.  Regulatory Assistance Project.
  * [Biagioni et al., 2021] Biagioni, D., Zhang, X., Wald, D., Vaidhynathan, D., Chintala, R., King, J., and Zamzam, A. S. (2021).  Powergridworld: A framework for multi-agent reinforcement learning in power systems.  arXiv preprint arXiv:2111.05969.
  * [Glorot and Bengio, 2010] Glorot, X. and Bengio, Y. (2010).  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings.
  * [Guerrero et al., 2018] Guerrero, J., Chapman, A. C., and Verbič, G. (2018).  Decentralized p2p energy trading under network constraints in a low-voltage network.  IEEE Transactions on Smart Grid, 10(5):5163–5173.
  * [Guerrero et al., 2020] Guerrero, J., Gebbran, D., Mhanna, S., Chapman, A. C., and Verbič, G. (2020).  Towards a transactive energy system for integration of distributed energy resources: Home energy management, distributed optimal power flow, and peer-to-peer energy trading.  Renewable and Sustainable Energy Reviews, 132:110000.
  * [Hatfield et al., 2013] Hatfield, J. W., Kominers, S. D., Nichifor, A., Ostrovsky, M., and Westkamp, A. (2013).  Stability and competitive equilibrium in trading networks.  Journal of Political Economy, 121(5):966–1005.
  * [Konda and Tsitsiklis, 1999] Konda, V. and Tsitsiklis, J. (1999).  Actor-critic algorithms.  Advances in neural information processing systems, 12.
  * [Liu and Wu, 2021] Liu, H. and Wu, W. (2021).  Online multi-agent reinforcement learning for decentralized inverter-based volt-var control.  IEEE Transactions on Smart Grid, 12(4):2980–2990.
  * [Liu et al., 2017a] Liu, H. J., Shi, W., and Zhu, H. (2017a).  Distributed voltage control in distribution networks: Online and robust implementations.  IEEE Transactions on Smart Grid, 9(6):6106–6117.
  * [Liu et al., 2017b] Liu, N., Yu, X., Wang, C., Li, C., Ma, L., and Lei, J. (2017b).  Energy-sharing model with price-based demand response for microgrids of peer-to-peer prosumers.  IEEE Transactions on Power Systems, 32(5):3569–3583.
  * [Liu et al., 2019] Liu, Y., Wu, L., and Li, J. (2019).  Peer-to-peer (P2P) electricity trading in distribution systems of the future.  The Electricity Journal, 32(4):2–6.
  * [Lowe et al., 2017] Lowe, R., Wu, Y. I., Tamar, A., Harb, J., Pieter Abbeel, O., and Mordatch, I. (2017).  Multi-agent actor-critic for mixed cooperative-competitive environments.  Advances in neural information processing systems, 30.
  * [Lu et al., 2019] Lu, X., Xiao, X., Xiao, L., Dai, C., Peng, M., and Poor, H. V. (2019).  Reinforcement learning-based microgrid energy trading with a reduced power plant schedule.  IEEE Internet of Things Journal, 6(6):10728–10737.
  * [Morstyn et al., 2018] Morstyn, T., Teytelboym, A., and McCulloch, M. D. (2018).  Bilateral contract networks for peer-to-peer energy trading.  IEEE Transactions on Smart Grid, 10(2):2026–2035.
  * [Qiu et al., 2021] Qiu, D., Wang, J., Wang, J., and Strbac, G. (2021).  Multi-agent reinforcement learning for automated peer-to-peer energy trading in double-side auction market.  In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI, pages 2913–2920.
  * [Qu et al., 2019] Qu, C., Mannor, S., Xu, H., Qi, Y., Song, L., and Xiong, J. (2019).  Value propagation for decentralized networked deep multi-agent reinforcement learning.  Advances in Neural Information Processing Systems, 32.
  * [Saha and Johnson, 2016] Saha, S. and Johnson, N. (2016).  Modeling and simulation in xendee ieee 13 node test feeder.  XENDEE, IRA A. Fulton School of Engineering, Arizona State University.
  * [Samende et al., 2022] Samende, C., Cao, J., and Fan, Z. (2022).  Multi-agent deep deterministic policy gradient algorithm for peer-to-peer energy trading considering distribution network constraints.  Applied Energy, 317:119123.
  * [Schulman et al., 2017] Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O. (2017).  Proximal policy optimization algorithms.  arXiv preprint arXiv:1707.06347.
  * [Sorin et al., 2018] Sorin, E., Bobo, L., and Pinson, P. (2018).  Consensus-based approach to peer-to-peer electricity markets with product differentiation.  IEEE Transactions on Power Systems, 34(2):994–1004.
  * [Tushar et al., 2020] Tushar, W., Saha, T. K., Yuen, C., Smith, D., and Poor, H. V. (2020).  Peer-to-peer trading in electricity networks: An overview.  IEEE Transactions on Smart Grid, 11(4):3185 – 3200.
  * [Ye et al., 2021] Ye, Y., Tang, Y., Wang, H., Zhang, X.-P., and Strbac, G. (2021).  A scalable privacy-preserving multi-agent deep reinforcement learning approach for large-scale peer-to-peer transactive energy trading.  IEEE Transactions on Smart Grid, 12(6):5185–5200.
  * [Zhang et al., 2018] Zhang, K., Yang, Z., Liu, H., Zhang, T., and Basar, T. (2018).  Fully decentralized multi-agent reinforcement learning with networked agents.  In International Conference on Machine Learning, pages 5872–5881. PMLR.
  * [Zhao et al., 2022] Zhao, Z., Feng, C., and Liu, A. L. (2022).  Comparisons of auction designs through multi-agent learning in peer-to-peer energy trading.  IEEE Transactions on Smart Grid.  Online first.
  * [Zhou et al., 2018] Zhou, Y., Wu, J., and Long, C. (2018).  Evaluation of peer-to-peer energy sharing mechanisms based on a multiagent simulation framework.  Applied energy, 222:993–1022.