# X-ray Scatter Estimation Using Deep Splines

Philipp Roser, Annette Birkhold, Alexander Preuhs, Christopher Syben, Lina
Felsner, Elisabeth Hoppe, Norbert Strobel, Markus Korwarschik, Rebecca Fahrig,
Andreas Maier P. Roser, A. Preuhs, C. Syben, L. Felsner, E. Hoppe, and A.
Maier are with the Pattern Recognition Lab, Department of Computer Science,
Friedrich-Alexander Universität Erlangen-Nürnberg, Erlangen, Germany. P. Roser
is funded by the Erlangen Graduate School in Advanced Optical Technologies
(SAOT), Friedrich-Alexander Universität Erlangen-Nürnberg, Erlangen, Germany.
A. Maier is principal investigator at the SAOT. A. Birkhold, M. Kowarschik,
and R. Fahrig are employees of Siemens Healthcare GmbH, 91301 Forchheim,
Germany. N. Strobel is with the Institute of Medical Engineering Schweinfurt,
University of Applied Sciences Würzburg-Schweinfurt, 97421 Schweinfurt,
Germany.

###### Abstract

Algorithmic X-ray scatter compensation is a desirable technique in flat-panel
X-ray imaging and cone-beam computed tomography. State-of-the-art U-net based
image translation approaches yielded promising results. As there are no
physics constraints applied to the output of the U-Net, it cannot be ruled out
that it yields spurious results. Unfortunately, those may be misleading in the
context of medical imaging. To overcome this problem, we propose to embed
B-splines as a known operator into neural networks. This inherently limits
their predictions to well-behaved and smooth functions. In a study using
synthetic head and thorax data as well as real thorax phantom data, we found
that our approach performed on par with U-net when comparing both algorithms
based on quantitative performance metrics. However, our approach not only
reduces runtime and parameter complexity, but we also found it much more
robust to unseen noise levels. While the U-net responded with visible
artifacts, our approach preserved the X-ray signal’s frequency
characteristics.

###### Index Terms:

Approximation, B-spline, Neural network, X-ray scatter.

##  I Introduction

Xray transmission imaging enables advanced diagnostic imaging or facilitates
instrument guidance during minimally-invasive interventions. Unfortunately,
primary photons scattered by the patient’s body impair image quality at the
detector. Especially cone-beam X-ray systems using flat-panel detectors suffer
from this phenomenon due to the large X-ray field and the resulting high
scatter-to-primary ratio. Besides contrast degradation in X-ray projections,
artifacts in cone-beam computed tomography (CBCT) deteriorate the image
quality in the reconstructed volume. Since the advent of CBCT, various methods
to compensate for scatter have been developed [1, 2]. These can be categorized
into physical and algorithmic scatter compensation methods.

###  I-A Physical Scatter Compensation

Physical scatter compensation refers to the direct manipulation of the X-ray
field to either suppress scattered photons from reaching the detector, or
modulate them to enable differentiation into scattered and primary photons. A
simple yet effective approach is to increase the patient to detector distance,
the so-called air gap [3, 4, 5, 6]. Since the number of scattered photons
depends on the irradiated volume, the scatter-to-primary ratio can be reduced
by shaping the X-ray field to a small volume of interest. In slit scanning,
the stitching of these smaller volumes of interest yields a normal-sized
reconstruction [7, 8]. Although slit scanning reduces patient dose and
scattered radiation, the prolonged acquisition time complicates motion
compensation and X-ray tube heat management [9]. Since relying on large air
gaps or strict collimation limits the flexibility of the imaging system,
clinical systems are equipped with anti-scatter grids that physically block
scattered photons [5, 10, 11, 12, 13, 14]. Although state-of-the-art for
clinical CBCT, anti-scatter grids have disadvantages. For example, they can
increase radiation exposure [15] and lead to Moiré and grid line artifacts
[16]. Overall, the efficiency and effectiveness of anti-scatter grids depends
on multiple factors, e.g., air gap, photon energy, source-to-detector
distance, and detector resolution [2, 17, 18]. In contrast to detector-side
grids, primary modulation follows a different operation principle [19, 20,
21]. In the shadow of this modulator grid, almost no primary radiation is
measured. Instead, information of the scatter properties of the irradiated
volume is encoded. Using demodulation algorithms, the X-ray projection and the
associated scatter distribution can be distinguished. Due to its complex
structure, primary modulation, especially with C-arm systems, has not made its
way into the clinical practice. In conclusion, hardware-based solutions induce
additional manufacturing costs and limit flexibility. Also, since most
approaches need algorithmic support or look-up tables, software-only solutions
are highly desirable [1, 2].

###  I-B Algorithmic Scatter Compensation

Algorithmic scatter compensation approaches typically try to estimate the
scatter signal from the acquired X-ray projections first. The estimated
scatter image is then subtracted to obtain the primary signal. For CBCT, using
Monte Carlo (MC) methods or Boltzmann transport solvers are common approaches
[22, 23, 24]. Both require an initial reconstruction. Their computational
complexity can partially be coped with using dedicated hardware and variance
reduction techniques [25] or by using a coarse simulation to derive the
scattering model online [26]. Trading off effectiveness for efficiency, model-
based approaches are typically preferred for clinical applications. Such
approaches rely on either simplified physical or analytical models [27, 28,
29] or convolution kernels [30, 31, 32] in combination with an iterative
correction scheme. Although being fast, their ability to generalize to
different imaging settings is limited.

Recently, learning-based methods found their way to modeling physical
processes [33, 34, 35]. Especially the application of the U-net [36] to the
task of scatter estimation, referred to as deep scatter estimation, yields
superior results compared to kernel-based or MC methods, respectively [37].
The U-net based method has the potential to become the new gold standard in
the domain of scatter correction. However, there are several drawbacks to
employing a deep U-net. First, deep convolutional neural networks are
challenging to comprehend in their operating principle. Given the high
parameter complexity of such a U-net, large amounts of training samples are
mandatory to arrive at a robust solution. Second, to maintain a fast
performance, the U-net requires a dedicated top tier graphics processing unit
(GPU), which might not be universally available. Third, the U-net is a
universal function approximator without including any known physical
characteristics of scattered radiation.

To come up with an alternative, it appears attractive to build on the rich
prior knowledge about X-ray scatter. Also, it has been shown that the
incorporation of prior knowledge into neural networks reduces error bounds
[38] and simplifies the analysis of the deployed network [39]. Although
subject to photon Poisson noise, X-ray scatter is mainly a low-frequency
signal in the diagnostic energy regime [40, 41]. In a previous proof-of-
concept study, we exploited this property by approximating X-ray scatter using
a smooth bivariate B-spline and directly inferred spline coefficients using a
lean convolutional encoder architecture [42]. This (a) allowed us to reduce
the number of parameters and computational complexity tremendously, (b)
ensured that no high-frequency details of diagnostic value can get
manipulated, while (c) still meeting the high accuracy of the U-net.

###  I-C Contributions

In this work, we extend the idea of approximating X-ray scatter as bivariate
B-spline by integrating its evaluation into the computational graph of a
neural network. We analyze our approach in depth and compare it to several
U-net realizations in a nested cross-validation study using synthetic head and
thorax X-ray image data. Besides benchmarking the parameter and computational
complexity, we investigate all networks in terms of their frequency response,
the power spectral density of the co-domain, and the response to unseen noise
levels. Finally, we demonstrate how our method performs when applied to real
data in an anthropomorphic thorax phantom study.

##  II Materials and Methods

###  II-A X-ray Projection Formation using B-splines

In general, in both techniques, X-ray fluoroscopy and CBCT it is assumed that
primary photons either reach the detector along a straight line or are
absorbed completely. Therefore, the observed X-ray projection is typically
described in terms of its primary component
$\boldsymbol{I}_{\text{p}}\in\mathbb{R}^{w\times h}$ with image width $w$ and
height $h$ in pixels. The intensity at pixel $(u,v)$ is given by the Beer-
Lambert law

${I}_{\text{p}}\left(u,v\right)=\int_{E}{I}_{0}\left(u,v,E\right)e^{-\int_{\lambda}\mu\left(\boldsymbol{s}+\lambda\cdot\boldsymbol{r}\left(u,v\right),E\right)d\lambda}dE\enspace.$
(1)

The polychromatic flat-field projection or X-ray spectrum $\boldsymbol{I}_{0}$
as well as the linear attenuation coefficient
$\mu\left(\boldsymbol{l},E\right)$ depend on the photon energy $E$. Note that,
due to the cone-beam geometry, $\boldsymbol{I}_{0}$ decreases towards the
borders and thus depends on the pixel position $(u,v)$. The linear photon
attenuation further depends on the media along the straight line
$\boldsymbol{l}:\mathbb{R}\mapsto\mathbb{R}^{3}$ from the X-ray source
$\boldsymbol{s}\in\mathbb{R}^{3}$ in direction
$\boldsymbol{r}\in\mathbb{R}^{3}$ to the pixel $(u,v)$. Unfortunately, in
reality, besides being photoelectrically absorbed, X-ray photons also undergo
Compton scattering, Rayleigh scattering, or multiple occurrences of both
effects. Therefore, a more realistic formation model $\boldsymbol{I}$ adds
scattered photons $\boldsymbol{I}_{\text{s}}$ leading to

${I}\left(u,v\right)={I}_{\text{p}}\left(u,v\right)+{I}_{\text{s}}\left(u,v\right)\enspace.$
(2)

Since, for CBCT, we are mostly interested in the low-frequency components of
the scatter, we neglect photon shot noise and detector noise below.

In a recent study [42], we experimentally showed that the main scatter
components, which are low-frequency for diagnostic X-rays, can be well
approximated by sparse bivariate B-splines with error rates below
$1\text{\,}\mathrm{\char 37\relax}$. The domain and co-domain of a B-spline of
degree $k$ and order $n=k+1$ are fully characterized by the knot vectors
$\boldsymbol{t}_{u}=(t_{u,1},t_{u,2},\dots t_{u,w_{c}-n})$ and
$\boldsymbol{t}_{v}=(t_{v,1},t_{v,2},\dots t_{v,h_{c}-n})$, as well as the
coefficient matrix $\boldsymbol{C}\in\mathbb{R}^{w_{\text{c}}\times
h_{\text{c}}}$ with width $w_{\text{c}}$ and height $h_{\text{c}}$,
respectively. Based on this, we can approximate the spatial scatter
distribution as a tensor product B-spline
$\tilde{\boldsymbol{I}}_{\text{s},n}$

$\tilde{I}_{\text{s},n}\left(u,v\right)=\sum_{i=1}^{w_{\text{c}}}\sum_{j=1}^{h_{\text{c}}}C\left(i,j\right)B_{i,n,\boldsymbol{t}_{u}}\left(u\right)B_{j,n,\boldsymbol{t}_{v}}\left(v\right)\,,$
(3)

with the basis splines $B$, which are zero for knots that do not affect the
spline at the pixel $(u,v)$. The basis splines are recursively defined by

$\displaystyle B_{i,n,\boldsymbol{t}}(x)$ $\displaystyle={}$
$\displaystyle\frac{x-t_{i}}{t_{i+n}-t_{i}}$ $\displaystyle\cdot
B_{i,n-1,\boldsymbol{t}}(x)$ (4) $\displaystyle+{}$
$\displaystyle\frac{t_{i+n+1}-x}{t_{i+n+1}-t_{i+1}}$ $\displaystyle\cdot
B_{i+1,n-1,\boldsymbol{t}}(x)\enspace,$

with

$B_{i,0,\boldsymbol{t}}(x)=\begin{cases}1,&\text{if }t_{i}\leq x<t_{i+1}\\\
0,&\text{otherwise}\end{cases}\enspace.$ (5)

Since the basis functions $B$ are equal to zero for coefficients that do not
contribute to a pixel $(u,v)$, only a $n\times n$ sub-grid of $\boldsymbol{C}$
needs to be considered for each pixel. Thus, using matrix notation, the tensor
product can be reformulated to

$\tilde{I}_{\text{s},n}\left(u,v\right)=\boldsymbol{u}^{\text{T}}\cdot\left[\boldsymbol{M}_{n,\boldsymbol{t}_{u}}\left(u\right)\right]\cdot\boldsymbol{C}_{\text{uv}}\cdot\left[\boldsymbol{M}_{n,\boldsymbol{t}_{v}}\left(v\right)\right]^{\text{T}}\cdot\boldsymbol{v}\enspace,$
(6)

with the coefficient patch $\boldsymbol{C}_{\text{uv}}\in\mathbb{R}^{n\times
n}$ impacting the pixel $(u,v)$, the general matrix representation of a
univariate B-spline
$\boldsymbol{M}_{n,\boldsymbol{t}}\left(\cdot\right)\in\mathbb{R}^{n\times n}$
[43], and the vectors
$\boldsymbol{u}=\begin{pmatrix}u^{0}&u^{1}&\cdots&u^{k}\end{pmatrix}^{\text{T}}$
and
$\boldsymbol{v}=\begin{pmatrix}v^{0}&v^{1}&\cdots&v^{k}\end{pmatrix}^{\text{T}}$.
The general matrix representation of B-splines allows for arbitrarily spaced
knot grids, i.e., in our case, it allows for endpoint interpolation, which
exposes a more convenient behavior towards the borders of the evaluation grid.
Besides the borders of the grid, we currently limit our approach to uniform
knot grids and cubic B-splines ($k=3$). Thus, for most pixels $(u,v)$, the
B-spline is evaluated using

$\boldsymbol{M}_{4}=\frac{1}{6}\begin{pmatrix}1&4&1&0\\\ -3&0&3&0\\\
3&-6&3&0\\\ 1&3&-3&1\end{pmatrix}\enspace.$ (7)

With using a fixed knot grid and evaluation grid, we can pre-calculate both
$\boldsymbol{u}^{\text{T}}\cdot\left[\boldsymbol{M}_{n,\boldsymbol{t}_{u}}\left(u\right)\right]$
and
$\boldsymbol{v}^{\text{T}}\cdot\left[\boldsymbol{M}_{n,\boldsymbol{t}_{v}}\left(v\right)\right]$
for each pixel $(u,v)$. Padding the resulting vectors with zeros to account
for all coefficients in $\boldsymbol{C}$ and stacking them row-wisely, we
obtain the evaluation matrices $\boldsymbol{U}_{n}\in\mathbb{R}^{w\times
w_{\text{c}}}$ and $\boldsymbol{V}_{n}\in\mathbb{R}^{h\times h_{\text{c}}}$,
respectively. Thus, we can calculate the cubic ($n=4$) scatter distribution
given the B-spline coefficients $\boldsymbol{C}$ via

$\tilde{\boldsymbol{I}}_{\text{s},4}=\boldsymbol{U}_{4}\cdot\boldsymbol{C}\cdot\boldsymbol{V}_{4}^{\text{T}}\enspace.$
(8)

Since neural networks can extract scatter distributions from a measured X-ray
projection either directly [37] or indirectly via B-spline coefficients [42]
we aim to combine deep learning with the matrix evaluation scheme. Using the
Kronecker product ‘$\otimes$’, the derivative
$\frac{\partial\tilde{\boldsymbol{I}}_{\text{s},4}}{\partial\boldsymbol{C}}$
is

$\frac{\partial\tilde{\boldsymbol{I}}_{\text{s},4}}{\partial\boldsymbol{C}}=\boldsymbol{V}_{4}\otimes\boldsymbol{U}_{4}\enspace.$
(9)

Thus, we can straightforwardly embed the B-spline evaluation into a
computational graph $f_{\boldsymbol{\theta}}:\mathbb{R}^{w\times
h}\mapsto\mathbb{R}^{w\times h}$ with parameters to train
$\boldsymbol{\theta}$ without breaking differentiability and thus back-
propagation.

###  II-B Network Architecture

Figure 1:  The general architecture of our proposed approach. First, a
(convolutional) encoder extracts the latent variables $\boldsymbol{Z}$ from
the input X-ray projection $\boldsymbol{I}$. Second, a bottleneck network maps
the latent space to bivariate spline coefficients $\boldsymbol{C}$. The
different choices for both networks are discussed later in more detail. Third,
we obtain the scatter estimate $\tilde{\boldsymbol{I}}_{\text{s},n}$ by
evaluating the spline coefficients $\boldsymbol{C}$ using the pre-calculated
evaluation grid defined by the pre-computed matrices $\boldsymbol{U}$ and
$\boldsymbol{V}$. The top (green) path refers to the forward pass, and the
bottom (orange) path to the backward pass, respectively.

In the following, we describe the devised neural network. The overall
architecture comprises a generic convolutional encoder followed by a
bottleneck network to infer spline coefficients from measured X-ray
projections, as proposed in our previous study [42]. This architecture is
depicted in Fig. 1. We employ a lean convolutional encoder
$f_{\boldsymbol{\theta}}:\mathbb{R}^{w\times
h}\mapsto\mathbb{R}^{w_{\text{c}}\times h_{\text{c}}}$ that consists of $d$
convolutional blocks. A block comprises two convolutional layers with $c$
feature channels and $3\times 3$ kernels. Each convolutional layer is followed
by a rectified linear unit (ReLU) [44] activation, and between two blocks,
$2\times 2$ average pooling is applied. As X-ray scatter is low-frequency and
to potentially limit the computational complexity of our model, we allow for
additional pooling layers before the first convolutional block. The encoder
$f_{\boldsymbol{\theta}}$ is completed by a $1\times 1$ convolution to build
the weighted sum of all channels $c$. Overall, based on convolutions only,
$f_{\boldsymbol{\theta}}$ only encodes a local scatter representation
$\boldsymbol{Z}$, since its receptive field does not necessarily cover the
whole input image. To establish a global context, we employ different types of
bottleneck networks $g_{\boldsymbol{\phi}}$ to map $\boldsymbol{Z}$ to spline
coefficients $\boldsymbol{C}$: (1) a constrained weighting matrix
$\boldsymbol{W}$ with $W_{i,j}>0$, (2) an unconstrained fully-connected layer
with ReLU activation (which merely relates to an unconstrained weighting
matrix with enforced non-negativity constraint), (3) two fully-connected
layers with ReLU activation, or (4) two additional convolutional blocks
followed by a fully-connected layer.

###  II-C Synthetic Dataset

Acquiring raw scatter-free X-ray projections and their scatter-contaminated
counterparts is tedious and time-consuming, especially at the scale to carry
out deep learning. As a solution, we leveraged the X-ray transport code MC-GPU
[25] to generate artificial pairs of scatter-contaminated and scatter-free
X-ray projections. We used openly available CT scans from The Cancer Imaging
Archive (TCIA) [45] as inputs to the simulation. In order to keep the
simulation time manageable, we selected 20 head scans from the HNSCC-3DCT-RT
dataset [46] and 15 thorax scans from the CT Lymph Nodes dataset [47],
respectively. To prepare the phantoms for MC simulation, we employed a basic
pre-processing pipeline [48] based on tissue and density estimation [49], and
connected component labeling [50]. For each CT scan, we simulated a stack of
260 X-ray projections ($w=1152$, $h=768$) over an angular range of
$200\text{\,}\mathrm{\SIUnitSymbolDegree}$. The source-to-isocenter and
source-to-detector distances are $785\text{\,}\mathrm{mm}$ and
$1300\text{\,}\mathrm{mm}$, respectively. Per X-ray projection, we simulated
$5\text{\times}{10}^{10}$ photons sampled from an $85\text{\,}\mathrm{kV}$
peak voltage tungsten spectrum. All projections were flat-field normalized.
Note that we used the data as is, without registering similar anatomies in a
common reference frame. To be more comparable to the related work [37], speed-
up training times, and suppress simulation noise, we down-sampled the
projections to $384\times 256$ pixels. Before the down-sampling, we applied
Gaussian filtering to the primary and scatter projections independently
($\sigma_{\text{p}}=2$, $\sigma_{\text{s}}=30$). Corresponding cross-sectional
slices were reconstructed on an isotropic ${256}^{3}$ grid with
$1\text{\,}{\mathrm{mm}}^{3}$ voxels.

###  II-D Real Dataset

To evaluate the proposed method on real data, we scanned the thorax of an
anthropomorphic phantom (PBU-60, Kyoto Kagaku Co., Ltd., Kyoto, Japan) using a
C-arm CBCT system (ARTIS icono floor, Siemens Healthineers AG, Forchheim,
Germany). In total, we acquired 12 datasets, three short scans for each full
view grid (referred to as full), and maximum collimation (referred to as
slit), both with and without an anti-scatter grid. Each short scan consists of
397 projections ($648\text{\times}472$ pixels,
$616\text{\,}\mathrm{\SIUnitSymbolMicro m}$ isotropic spacing) over an angular
range of $197.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ using
$85\text{\,}\mathrm{kV}$ peak tube voltage. The source-to-isocenter and
source-to-detector distances are $750\text{\,}\mathrm{mm}$ and
$1200\text{\,}\mathrm{mm}$, respectively. We reconstructed
$512\text{\times}512$ slices on an $484\text{\,}\mathrm{\SIUnitSymbolMicro m}$
isotropic voxel grid using an in-house reconstruction pipeline. We regard the
slit scan in conjunction with the anti-scatter grid as ground truth.

##  III Experiments

In the following, we describe the experiments carried out to evaluate our
proposed method. Since the U-net based approach outperformed other
computational scatter estimation methods by far [37], we considered different
configurations of the U-net as a baseline method. To account for our
relatively small training corpora, we evaluated our method and the baseline
using a nested cross-validation approach. For the head dataset (20 subjects),
we used a $4^{*}3$-fold cross-validation approach and, for the thorax dataset
(15 subjects), we used a $5^{*}4$-fold cross-validation approach. The real
dataset was only used for testing. To keep the training procedure manageable,
we divided the evaluation into meta-parameter search based on the scatter mean
absolute percentage error (MAPE), and further in-depth analysis.

###  III-A Meta-Parameter Search

Since clinic CBCT systems usually have fine-tuned acquisition protocols for
each anatomic region, we evaluated each network architecture for head and
thorax data separately. To fix the meta-parameters, we only used the synthetic
head dataset and the scatter MAPE as metric. We validated both our approach
and the U-net for different combinations of depths $d$ and feature channels
$c$ using Glorot initialization [51], which was also used by the baseline
[37]. Furthermore, we distinguished between a deep U-net (DU-net) and a
shallow U-net (SU-net), in which the number of feature maps is not doubled at
each level. For all experiments in this section, we performed a $4^{*}3$-fold
cross-validation and trained all networks using the adaptive moments optimizer
(Adam) [52] with an initial learning rate of ${10}^{-4}$ for 100 epochs. In
total, we trained 12 networks for each configuration and used the averaged
MAPE to assess their quality. To reduce the total computation time, we stopped
each training procedure when no significant performance increase was observed
for 20 epochs.

In a first step, we scanned different parametrizations of the network
architectures to find the most promising ones for in-depth comparison.
Overall, we tested both U-nets with $c=16$ and $d\in\\{4,5,6,7\\}$, and our
approach with $c=16$, $d\in\\{4,5,6\\}$, and additional pre-pooling
$p\in\\{0,1,2\\}$. In addition, to decouple the architecture of our
convolutional encoder $f_{\boldsymbol{\theta}}$ and the spline coefficient
dimensionality, we investigated the four bottleneck architectures
$g_{\boldsymbol{\phi}}$ as described in Sec. II-B.

###  III-B Qualitative and Quantitative Results

Based on the findings of the previous experiments, we adapted the learning
rate to ${10}^{-5}$ but kept the overall training routine. We separately
performed 4∗3-fold and 5∗4-fold cross-validations for the head and thorax
datasets, respectively. In addition to the scatter MAPE, we included the
structural similary index (SSIM) of scatter-compensated reconstructions with
respect to the ground truth in our evaluation.

###  III-C In-Depth Analysis

####  III-C1 Spectral Analysis

For clinical applications, data integrity is of utmost importance to ensure
that automated systems do not alter diagnostically relevant content. While the
predictions of neural networks may appear reasonable at first glance,
unrealistic perturbations can be unveiled by investigating the spectral
properties of (a) the predicted images [53] or (b) the neural network itself.
Therefore, we first investigate all networks’ performance concerning the
predicted scatter distributions’ power spectral density. Second, from scatter
estimation theory, we know that the scatter distribution can be recovered from
the measured signal by the convolution with a so-called scatter kernel [30].
This allows us to interpret a neural network for a specific pair of scatter
distribution and X-ray signal as a filtering operation and assess its
frequency response.

####  III-C2 Noise Analysis

Assessing a network’s robustness is inherently difficult given small training
corpora. Therefore, adversarial attacks are often used to expose weaknesses
deliberately. Since we trained all networks on noise-free data, testing them
on data with different unseen noise levels appeared appropriate. To this end,
we applied Poisson noise associated with different photon counts ranging from
${10}^{3}{10}^{5}$ to our head dataset and investigated the accuracy of the
networks’ predictions.

####  III-C3 Runtime Analysis

For interventional applications, the fast execution speed of computer programs
is essential. Therefore, we benchmarked the average execution time for all
networks for different batch sizes using a 12-core CPU (Intel(R) Xeon(R)
Silver 4116 CPU 2.10GHz).

####  III-C4 Real Data Analysis

To confirm the generalizability to real data, we tested the networks, which
performed best on our synthetic thorax dataset, on the real thorax dataset.
For reference, we considered various different scatter suppression techniques,
namely using an anti-scatter grid or slit scanning. As an almost scatter free
baseline, we considered the configuration to use both, an anti-scatter grid
and slit scanning. Note that we neither performed an intensity or a geometry
calibration in between the scans to account for the missing anti-scatter grid.

##  IV Results

###  IV-A Meta-Parameter Search

####  IV-A1 Network Parametrization

$\displaystyle 10^{5}$$\displaystyle 10^{6}$$\displaystyle
10^{7}$$\displaystyle 10^{8}$No. parameters$\displaystyle 6$$\displaystyle
7$$\displaystyle 8$$\displaystyle 9$MAPE [$\displaystyle\%$]DU-netSU-netOurs
Figure 2: Distribution of different network configurations for our approach
and the U-net in terms of absolute percentage errors averaged over all folds
and patients. Note that, due to visualization purposes, the standard deviation
encoded by the circular margins does not correspond to absolute values but
rather indicates the relative spread between the networks. The standard
deviation is in the range of $1\text{\,}\mathrm{\char
37\relax}5\text{\,}\mathrm{\char 37\relax}$. The x-axis refers to the number
of parameters in the convolutional layers.

Figure 2 establishes the relationship between the number of convolutional
parameters to train for each network to the averaged error rates of all folds
and patients. All networks’ error rates range between $6\text{\,}\mathrm{\char
37\relax}9\text{\,}\mathrm{\char 37\relax}$, and we observed that overall the
more compact networks outperform the DU-nets. Based on these findings, we
selected two spline networks ($c=16$, $(d,p)\in\\{(4,1),(5,0)\\}$) and four
U-nets (deep and shallow, $c=16$, $d\in\\{6,7\\}$) for further investigations.

####  IV-A2 Bottleneck

1234Fold05101520MAPE [$\displaystyle\%$]constrainedunconstrainedfc-netconv-net
Figure 3: Boxplots of the mean absolute percentage error (MAPE) of our
proposed method for different bottleneck architectures averaged over the
validation folds for each training fold. The horizontal lines indicate the
mean value over all test and validation folds.

$\displaystyle 0$$\displaystyle 200$$\displaystyle 0$$\displaystyle
50$constrained$\displaystyle 0$$\displaystyle 200$$\displaystyle
0$$\displaystyle 50$unconstrained

Figure 4: Normalized constrained (left) and unconstrained (right) weighting
matrices used to map the latent variables $\boldsymbol{Z}$ to spline
coefficients $\boldsymbol{C}$.

Figure 3 shows the results for the different bottleneck architectures. Our
proposed constrained weighting matrix, homogeneously initialized, achieves the
best results on average, even surpassing the convolutional bottleneck followed
by a fully-connected layer. The unconstrained fully-connected layer
architectures yield considerably worse results. Figure 4, which shows the
weighting matrix for the constrained and unconstrained case, substantiates
this finding. While the constrained matrix converges to a block circulant
matrix, which merely relates to a convolution, the unconstrained one hardly
resembles a sensible operation and is overall noisy.

###  IV-B Qualitative and Quantitative Results

Head data | Thorax data
---|---
1234Fold05101520MAPE [$\displaystyle\%$] |  12345Fold05101520MAPE [$\displaystyle\%$]
1234Fold0.950.960.970.991.00SSIM |  12345Fold0.950.960.970.991.00SSIM

Figure 5: Quantitative boxplots for our synthetic datasets. The top row shows
the mean absolute percentage errors (MAPE) with respect to the scatter ground
truth for each test fold. The bottom row shows structural similarity indices
(SSIM) between reconstructed volumes from the simulated ideal primary signal
and the scatter-corrected ones using neural networks. The horizontal lines
indicate the overall average across all folds. For the sake of clarity, we
only depict the best performing networks in each category. For the head
dataset, we depict the DU-net with $d=6$, SU-net with $d=7$, and our method
with $d=4$ and $p=1$. For the thorax dataset, we depict the DU-net with $d=7$,
SU-net with $d=6$, and our method with $d=5$ and $p=0$.

Figure 6: Selected subjects of both datasets. The results of the learning-
based approaches are given in terms of error maps using the absolute
percentage error (APE) for scatter and the absolute error (AE) for the
reconstructed slices. For convenience, the associated mean APE (MAPE) and mean
AE (MAE) are also provided below each output.

Figure IV-B shows fold-wise boxplots for the two synthetic datasets. We
observe similar error rates of approximately $5\text{\,}\mathrm{\char
37\relax}$ across all folds and network configurations regarding the predicted
scatter distributions for the head dataset. In general, all networks achieve
high SSIM values above 0.99 when comparing the reconstructed volumes from the
simulated ideal primary signal to their counterparts obtained using the
scatter-corrected projections. Note that for the head data, our approach
performed equally well for all folds, whereas the results obtained with the
U-nets varied more widely. Processing the thorax datasets, on the other hand,
was more challenging. Again, all networks performed comparably well with error
rates of approximately $7.5\text{\,}\mathrm{\char 37\relax}$. In comparison to
the head dataset, larger error margins and outliers were found. This can also
be seen with the SSIM, which is, on average, just below 0.98. Again, we find
that our proposed method is on par with U-net-like structures from a
quantitative point of view.

Looking at selected subjects of both datasets in Fig. 6 reveals the potential
advantages of the proposed approach. Modeling the predicted scatter as a
B-spline intrinsically limits the network output to smooth surfaces. In
contrast, both U-nets preserve some details of the input, especially the
shallow architecture. However, the reconstructed slices show no systematic
trend, and all methods can adequately recover the desired signal.

###  IV-C In-Depth Analysis

####  IV-C1 Spectral Analysis

$\displaystyle 10^{-2}$$\displaystyle 10^{-1}$$\displaystyle 10^{0}$Normalized
frequency$\displaystyle 10^{-2}$$\displaystyle 10^{-1}$$\displaystyle
10^{0}$Normalized powerGTDU-netSU-netOurs Figure 7: Power spectral density
plots for the X-ray scatter distributions of the ground truth (black solid),
ours (light orange dashed), and the U-net (light purple dotted) averaged over
all test folds. Note that the graphs associated with both U-nets are almost
identical.

We calculated the power spectral density by azimuthally averaging the
magnitudes of the 2D Fourier transform of X-ray scatter distributions. Also,
we averaged all power spectral densities of all projection, patients, and test
folds. The resulting densities plotted in Fig. 7 support our previous
findings. While the U-nets yield numerically sound predictions, they
systematically boost the high frequencies in the scatter distributions. Our
approach, however, preserves the real power spectral density over the whole
frequency spectrum.

AmplitudeGTDU-netSU-
netOurs$\displaystyle-3$$\displaystyle-2$$\displaystyle-1$$\displaystyle
0$Phase$\displaystyle 0$$\displaystyle 1$$\displaystyle 2$ Figure 8: Frequency
responses of the systems for one test patient averaged over all validation
folds. From top to bottom: Normalized $\log$-amplitude, phase in radians. From
left to right: ideal system (GT), our spline-net (Ours), deep U-net (DU-net),
shallow U-net (SU-net).

As mentioned above, conventional convolutional neural networks can be
interpreted for a single input image in terms of a filtering operation. Thus,
we divided the Fourier transform of the output by the Fourier transform of the
input to obtain the respective frequency responses. Figure 8 shows the
amplitude and phase of an ideal system’s frequency responses, for our method
and both U-nets, averaged over all projections for one patient. Our spline-
net’s frequency response was closer to the ideal frequency response in both,
amplitude and phase. However, we observed a noticeable intensity shift in our
method. The U-nets, in contrast, both exhibited larger deviations in the
patterns of amplitude and phase, indicating that their represented operation
is less predictable.

####  IV-C2 Noise Analysis

$\displaystyle 10^{3}$$\displaystyle 10^{4}$$\displaystyle 10^{5}$Photon
count$\displaystyle 5.0$$\displaystyle 7.5$$\displaystyle 10.0$$\displaystyle
12.5$$\displaystyle 15.0$$\displaystyle 17.5$MAPE [$\displaystyle\%$]DU-netSU-
netOurs Figure 9: Error rates of predicted scatter distributions for different
noise levels averaged over all validation and test fold networks. All networks
have been trained with noise-free data. Note that a lower photon count relates
to a higher noise level.

Figure 9 shows plots of the networks’ error rates when confronted with noise.
We could confirm that the U-net is very sensitive to unseen noise levels in
both configurations, whereas our approach performs more robustly.

####  IV-C3 Runtime Analysis

12481632Batch size$\displaystyle 0$$\displaystyle 25$$\displaystyle
50$$\displaystyle 75$$\displaystyle 100$$\displaystyle 125$$\displaystyle
150$Runtime [ms]DU-net 6DU-net 7SU-net 6SU-net 7Ours 4\1Ours 5\0 Figure 10:
Runtimes per projection in $\mathrm{ms}$ for different architectures
(specified by the depth and optional pre-pooling, $d$ \ $p$) and batch sizes.

Figure 10 compiles the inference speed of all considered networks. As expected
from the parameter complexity, our approach is the fastest with
$4\text{\,}\mathrm{ms}30\text{\,}\mathrm{ms}$ and therefore $1.78.5$ times
faster than the SU-net with $34\text{\,}\mathrm{ms}50\text{\,}\mathrm{ms}$.
The DU-net is the slowest with
$89\text{\,}\mathrm{ms}144\text{\,}\mathrm{ms}$.

####  IV-C4 Real Data Analysis

ReconstructionGrid+SlitGrid+FullSlitFullDU-netSU-
netOursDifference$\displaystyle 39.31\pm 0.21$$\displaystyle 58.46\pm
0.09$$\displaystyle 123.84\pm 0.32$$\displaystyle 62.86\pm 0.25$$\displaystyle
64.52\pm 0.23$$\displaystyle 63.97\pm 0.19$ Figure 11: Central slices of
reconstructed volumes for different scatter compensation strategies and errors
taken with respect to the reconstruction result resulting from the grid + slit
data acquisition. Grid refers to employing a conventional anti-scatter grid.
Slit refers to the most narrow collimator setting available on the X-ray
imaging system. Full refers to using no collimation at all. Both, the
reconstructed slices as well as the difference images are shown using a gray
level window $[-1000,1000]$ HU. For convenience the absolute average HU errors
and the associated standard deviations with respect to three consecutive
measurements are provided.

As shown in Fig. 11, all methods were able to compensate for scatter artifacts
well. However, employing an anti-scatter grid still yielded the lowest overall
error in Hounsfield units (HU) as compared to using slit collimation in
addition. The learning-based approaches performed about as well as the slit
scanning without an anti-scatter grid, which is in accordance to previous
findings [37]. Overall, the networks achieved similar error rates, and no
systematic trend was observable.

##  V Discussion

X-ray scatter is a major source of artifacts in interventional CBCT. Deep-
learning-based approaches have shown the potential to outperform conventional
physical or algorithmic approaches to scatter compensation. Without
incorporating prior information, scatter distributions can be inferred from
the measured X-ray signal [37]. However, data integrity and robustness are
critical aspects of clinical imaging, which can be violated by deep neural
networks [54].

To ensure sound scatter estimates, we proposed to embed bivariate B-splines in
neural networks to constrain their co-domain to smooth results [42]. By
reformulating the spline evaluation in terms of matrix multiplications, we
were able to integrate B-splines with neural networks without further ado such
that end-to-end training was feasible. In an extensive cross-validation using
synthetic data, we showed that our proposed lean convolutional encoder using
B-spline evaluation performs on par with several U-net based architectures. We
substantiated this finding in a first phantom study. There, our approach
performed basically as well as the U-net and the slit scanning technique
(without anti-scatter grid), which was used as a baseline before [37].

Note, however, that the proposed method offers several advantages not present
in the U-net architectures. First, we considerably lowered the parameter and
runtime complexity, rendering our method suitable for a variety of hardware.
More importantly, we verified that our spline-based approach ensures data
integrity concerning the power spectral density of scatter estimates and the
overall network’s frequency response. This property ensures that no high-
frequency details, which relate to anatomic structures or pathologies, are
altered. In comparison, the U-net considerably changes the concerning spectral
contents, which was already shown for neural networks containing up-
convolutions [53]. As our spline network corresponds to a low-pass filtering
operation, it is robust towards noise even when trained on noise-free data.

While we implemented the U-net baselines to the best of our knowledge, our
error rates are higher than previously reported [37]. Potential reasons
include but are not limited to (a) different simulation codes, (b) our smaller
training corpora, and (c) the heterogeneity of our data. Our simulation setup
currently assumes an ideal detector, and we do not consider the domain shift
between synthetic and real data.

For future work, we indicate several research directions for either method or
data and experiments. First, since our preferred bottleneck weighting matrix
converges to a block circulant matrix, a fixed representation is desirable to
further reduce the number of trainable parameters and increase the
plausibility of our approach. For instance, training both the encoder and the
bottleneck separately or introducing additional constraints are promising
approaches to do so. Second, replacing the bottleneck fully-connected layer
with a small U-net reduces the number of parameters while still covering the
entire latent space with its receptive field. Third, since our network is end-
to-end trainable, it appears reasonable to include the reconstruction into the
computational graph to calculate the loss function in the CBCT domain [55,
56]. Last but not least, we believe our approach applies to low-frequency
signal estimation and correction in general, e.g., bias field correction in
magnetic resonance imaging, ultrasound imaging, or microscopy techniques.

##  VI Conclusion

Embedding B-splines in neural networks ensures data integrity for low-
frequency signals. This reduces the number of network parameters needed to
arrive at physically sensible results, and thanks to the reduced parameter
set, network inference can be made faster.

## Disclaimer

The concepts and information presented in this article are based on research
and are not commercially available.

## References

  * [1] E.-P. Rührnschopf and K. Klingenbeck, “A general framework and review of scatter correction methods in x-ray cone-beam computerized tomography. Part 1: Scatter compensation approaches,” _Med Phys_ , vol. 38, no. 7, pp. 4296–4311, 2011.
  * [2] ——, “A general framework and review of scatter correction methods in cone-beam CT. Part 2: Scatter estimation approaches,” _Med Phys_ , vol. 38, no. 9, pp. 5186–5199, 2011.
  * [3] F. Groedel and R. Wachter, “Bedeutung der Röhren-Fern- und Platten-Abstandsaufnahmen,” _Verh d Dt Röntgengesellschaft_ , vol. 17, pp. 134–135, 1926.
  * [4] F. Janus, “Die Bedeutung der gestreuten ungerichteten Röntgenstrahlen fur die Bildwirkung und ihre Beseitigung,” _Verh d Dt Röntgengesellschaft_ , vol. 17, pp. 136–139, 1926.
  * [5] U. Neitzel, “Grids or air gaps for scatter reduction in digital radiography: A model calculation,” _Med Phys_ , vol. 19, no. 2, pp. 475–481, 1992.
  * [6] J. Persliden and G. A. Carlsson, “Scatter rejection by air gaps in diagnostic radiology. calculations using a Monte Carlo collision density method and consideration of molecular interference in coherent scattering,” _Phys Med Biol_ , vol. 42, no. 1, pp. 155–175, 1997.
  * [7] G. T. Barnes, “Contrast and scatter in x-ray imaging.” _RadioGraphics_ , vol. 11, no. 2, pp. 307–323, 1991.
  * [8] R. Bhagtani and T. G. Schmidt, “Simulated scatter performance of an inverse-geometry dedicated breast CT system,” _Med Phys_ , vol. 36, no. 3, pp. 788–796, 2009.
  * [9] T. G. Schmidt, R. Fahrig, N. J. Pelc, and E. G. Solomon, “An inverse-geometry volumetric CT system with a large-area scanned source: A feasibility study,” _Med Phys_ , vol. 31, no. 9, pp. 2623–2627, 2004.
  * [10] G. Bucky, “Über die Ausschaltung der im Objekt entstehenden Sekundärstrahlen bei Röntgenstrahlen,” _Verh d Dt Röntgengesellschaft_ , vol. 9, pp. 30–32, 1913.
  * [11] W. A. Kalender, “Calculation of x-ray grid characteristics by Monte Carlo methods,” _Phys Med Biol_ , vol. 27, no. 3, pp. 353–361, 1982.
  * [12] H.-P. Chan and K. Doi, “Investigation of the performance of antiscatter grids: Monte Carlo simulation studies,” _Phys Med Biol_ , vol. 27, no. 6, pp. 785–803, 1982.
  * [13] H.-P. Chan, Y. Higashida, and K. Doi, “Performance of antiscatter grids in diagnostic radiology: Experimental measurements and Monte Carlo simulation studies,” _Med Phys_ , vol. 12, no. 4, pp. 449–454, 1985.
  * [14] H.-P. Chan, K. L. Lam, and Y. Wu, “Studies of performance of antiscatter grids in digital radiography: Effect on signal-to-noise ratio,” _Med Phys_ , vol. 17, no. 4, pp. 655–664, 1990.
  * [15] H. Aichinger, J. Dierker, S. Joite-Barfuss, and M. Saebel, _Radiation exposure and image quality in X-ray diagnostic radiology_. Berlin, New York: Springer, 2004.
  * [16] D. M. Gauntt and G. T. Barnes, “Grid line artifact formation: A comprehensive theory,” _Med Phys_ , vol. 33, no. 6, pp. 1668–1677, 2006.
  * [17] V. Singh, A. Jain, D. R. Bednarek, and S. Rudin, “Limitations of anti-scatter grids when used with high resolution image detectors,” in _Proc SPIE Int Soc Opt Eng_ , B. R. Whiting and C. Hoeschen, Eds., vol. 9033, International Society for Optics and Photonics. SPIE, 2014, pp. 1618–1626.
  * [18] R. Rana, A. Jain, A. Shankar, D. R. Bednarek, and S. Rudin, “Scatter estimation and removal of anti-scatter grid-line artifacts from anthropomorphic head phantom images taken with a high resolution image detector,” in _Proc SPIE Int Soc Opt Eng_ , D. Kontos and T. G. Flohr, Eds., vol. 9783, International Society for Optics and Photonics. SPIE, 2016, pp. 1619–1628.
  * [19] A. Bani-Hashemi, E. Blanz, J. Maltz, D. Hristov, and M. Svatos, “Tu-d-i-611-08: Cone beam x-ray scatter removal via image frequency modulation and filtering,” _Med Phys_ , vol. 32, no. 6, pp. 2093–2093, 2005\.
  * [20] L. Zhu, N. R. Bennett, and R. Fahrig, “Scatter correction method for x-ray CT using primary modulation: Theory and preliminary results,” _IEEE Trans Med Imaging_ , vol. 25, no. 12, pp. 1573–1587, 2006.
  * [21] B. Bier, M. Berger, A. Maier, M. Kachelrieß, L. Ritschl, K. Müller, J.-H. Choi, and R. Fahrig, “Scatter correction using a primary modulator on a clinical angiography c-arm CT system,” _Med Phys_ , vol. 44, no. 9, pp. 125–137, 2017.
  * [22] W. Zbijewski and F. J. Beekman, “Fast scatter estimation for cone-beam x-ray CT by combined Monte Carlo tracking and richardson-lucy fitting,” in _IEEE Nucl Sci Symp Conf Rec_ , vol. 5, 2004, pp. 2774–2777.
  * [23] G. Poludniowski, P. M. Evans, V. N. Hansen, and S. Webb, “An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT,” _Phys Med Biol_ , vol. 54, no. 12, pp. 3847–3864, 2009.
  * [24] A. Wang, A. Maslowski, P. Messmer, M. Lehmann, A. Strzelecki, E. Yu, P. Paysan, M. Brehm, P. Munro, J. Star-Lack, and D. Seghers, “Acuros cts: A fast, linear boltzmann transport equation solver for computed tomography scatter – part ii: System modeling, scatter correction, and optimization,” _Med Phys_ , vol. 45, no. 5, pp. 1914–1925, 2018.
  * [25] A. Badal and A. Badano, “Accelerating Monte Carlo simulations of photon transport in a voxelized geometry using a massively parallel graphics processing unit,” _Med Phys_ , vol. 36, no. 11, pp. 4878–4880, 2009.
  * [26] M. Baer and M. Kachelrieß, “Hybrid scatter correction for CT imaging,” _Phys Med Biol_ , vol. 57, no. 21, pp. 6849–6867, 2012.
  * [27] W. Swindell and P. M. Evans, “Scattered radiation in portal images: A Monte Carlo simulation and a simple physical model,” _Med Phys_ , vol. 23, no. 1, pp. 63–73, 1996.
  * [28] W. Yao and K. W. Leszczynski, “An analytical approach to estimating the first order scatter in heterogeneous medium. ii. a practical application,” _Med Phys_ , vol. 36, no. 7, pp. 3157–3167, 2009.
  * [29] M. Meyer, W. A. Kalender, and Y. Kyriakou, “A fast and pragmatic approach for scatter correction in flat-detector CT using elliptic modeling and iterative optimization,” _Phys Med Biol_ , vol. 55, no. 1, pp. 99–120, 2009\.
  * [30] B. Ohnesorge, T. Flohr, and K. Klingenbeck-Regn, “Efficient object scatter correction algorithm for third and fourth generation CT scanners,” _Eur Radiol_ , vol. 9, no. 3, pp. 563–569, Mar 1999.
  * [31] H. Li, R. Mohan, and X. R. Zhu, “Scatter kernel estimation with an edge-spread function method for cone-beam computed tomography imaging,” _Phys Med Biol_ , vol. 53, no. 23, pp. 6729–6748, 2008.
  * [32] M. Sun and J. M. Star-Lack, “Improved scatter correction using adaptive scatter kernel superposition,” _Phys Med Biol_ , vol. 55, no. 22, pp. 6695–6720, 2010.
  * [33] T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” in _Adv Neural Inf Process Syst_ , S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, Eds. Curran Associates, Inc., 2018, pp. 6571–6583.
  * [34] P. Roser, X. Zhong, A. Birkhold, N. Strobel, M. Kowarschik, R. Fahrig, and A. Maier, “Physics-driven learning of x-ray skin dose distribution in interventional procedures,” _Med Phys_ , vol. 46, no. 10, pp. 4654–4665, 2019.
  * [35] P. Roser, X. Zhong, A. Birkhold, A. Preuhs, C. Syben, E. Hoppe, N. Strobel, M. Kowarschik, R. Fahrig, and A. Maier, “Simultaneous estimation of x-ray back-scatter and forward-scatter using multi-task learning,” in _International Conference on Medical Image Computing and Computer-Assisted Intervention_. Springer, 2020, pp. 199–208.
  * [36] O. Ronneberger, P. Fischer, and T. Brox, “U-Net: Convolutional networks for biomedical image segmentation,” in _Med Image Comput Comput Assist Interv_ , ser. Lecture notes in computer science, vol. 9351. Springer, 2015, pp. 234–241.
  * [37] J. Maier, E. Eulig, T. Vöth, M. Knaup, J. Kuntz, S. Sawall, and M. Kachelrieß, “Real-time scatter estimation for medical CT using the deep scatter estimation: Method and robustness analysis with respect to different anatomies, dose levels, tube voltages, and data truncation,” _Med Phys_ , vol. 46, no. 1, pp. 238–249, 2019.
  * [38] A. Maier, C. Syben, B. Stimpel, T. Würfl, M. Hoffmann, F. Schebesch, W. Fu, L. Mill, L. Kling, and S. Christiansen, “Learning with known operators reduces maximum error bounds,” _Nat Mach Intell_ , vol. 1, pp. 373–380, 2019\.
  * [39] B. Stimpel, C. Syben, F. Schirrmacher, P. Hoelter, A. Dörfler, and A. Maier, “Multi-modal deep guided filtering for comprehensible medical image processing,” _IEEE Trans Med Imaging_ , vol. 39, no. 5, pp. 1703–1711, 2019.
  * [40] G. H. Glover, “Compton scatter effects in CT reconstructions,” _Med Phys_ , vol. 9, no. 6, pp. 860–867, 1982.
  * [41] J. M. Boone and J. A. Seibert, “An analytical model of the scattered radiation distribution in diagnostic radiology,” _Med Phys_ , vol. 15, no. 5, pp. 721–725, 1988.
  * [42] P. Roser, A. Birkhold, A. Preuhs, C. Syben, N. Strobel, M. Korwarschik, R. Fahrig, and A. Maier, “Deep scatter splines: Learning-based medical x-ray scatter estimation using b-splines,” in _The 6th International Conference on Image Formation in X-Ray Computed Tomography (CT-Meeting)_ , 2020\.
  * [43] K. Qin, “General matrix representations for b-splines,” _The Visual Computer_ , vol. 16, no. 3-4, pp. 177–186, 2000.
  * [44] X. Glorot, A. Bordes, and Y. Bengio, “Deep sparse rectifier neural networks,” in _Proceedings of the fourteenth international conference on artificial intelligence and statistics_ , 2011, pp. 315–323.
  * [45] K. Clark, B. Vendt, K. Smith, J. Freymann, J. Kirby, P. Koppel, S. Moore, S. Phillips, D. Maffitt, M. Pringle, L. Tarbox, and F. Prior, “The Cancer Imaging Archive (TCIA): Maintaining and operating a public information repository,” _J Digit Imaging_ , vol. 26, pp. 1045–1057, 07 2013.
  * [46] T. Bejarano, M. De Ornelas Couto, and I. Mihaylov, “Head-and-neck squamous cell carcinoma patients with CT taken during pre-treatment, mid-treatment, and post-treatment dataset. the cancer imaging archive,” 2018. [Online]. Available: http://doi.org/10.7937/K9/TCIA.2018.13upr2xf
  * [47] H. Roth, L. Le, S. Ari, K. Cherry, J. Hoffman, S. Wang, and R. Summers, “A new 2.5 d representation for lymph node detection in ct. the cancer imaging archive,” 2018. [Online]. Available: http://doi.org/10.7937/K9/TCIA.2015.AQIIDCNM
  * [48] P. Roser, A. Birkhold, A. Preuhs, B. Stimpel, C. Syben, N. Strobel, M. Kowarschik, R. Fahrig, and A. Maier, “Fully-automatic CT data preparation for interventional x-ray skin dose simulation,” in _Bildverarbeitung für die Medizin 2020_. Springer, 2020, pp. 125–130.
  * [49] W. Schneider, T. Bortfeld, and W. Schlegel, “Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions,” _Phys Med Biol_ , vol. 45, no. 2, pp. 459–478, 2000.
  * [50] L. He, X. Ren, Q. Gao, X. Zhao, B. Yao, and Y. Chao, “The connected-component labeling problem: A review of state-of-the-art algorithms,” _Pattern Recognit_ , vol. 70, pp. 25–43, 2017.
  * [51] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , 2010, pp. 249–256.
  * [52] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _3rd Int. Conf. on Learning Representations_ , Y. Bengio and Y. LeCun, Eds., 2015.
  * [53] R. Durall, M. Keuper, and J. Keuper, “Watch your up-convolution: Cnn based generative deep neural networks are failing to reproduce spectral distributions,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 7890–7899.
  * [54] Y. Huang, T. Würfl, K. Breininger, L. Liu, G. Lauritsch, and A. Maier, “Some investigations on robustness of deep learning in limited angle tomography,” in _Med Image Comput Comput Assist Interv_. Springer, 2018, pp. 145–153.
  * [55] C. Syben, M. Michen, B. Stimpel, S. Seitz, S. Ploner, and A. K. Maier, “Pyro-NN: Python reconstruction operators in neural networks,” _Med Phys_ , vol. 46, no. 11, pp. 5110–5115, 2019.
  * [56] C. Syben, B. Stimpel, P. Roser, A. Dörfler, and A. Maier, “Known operator learning enables constrained projection geometry conversion: Parallel to cone-beam for hybrid MR/x-ray imaging,” _IEEE Trans Med Imaging_ , 2020\.


# Non-Stationary Policy Learning for Multi-Timescale Multi-Agent Reinforcement
Learning

Patrick Emami, Xiangyu Zhang, David Biagioni and Ahmed S. Zamzam This work was
authored in part by the National Renewable Energy Laboratory (NREL), operated
by Alliance for Sustainable Energy, LLC, for the U.S. Department of Energy
(DOE) under Contract No. DE-AC36-08GO28308. This work was supported by the
Laboratory Directed Research and Development (LDRD) Program at NREL. The views
expressed in the article do not necessarily represent the views of the DOE or
the U.S. Government. The U.S. Government retains and the publisher, by
accepting the article for publication, acknowledges that the U.S. Government
retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or
reproduce the published form of this work, or allow others to do so, for U.S.
Government purposes. A portion of this research was performed using
computational resources sponsored by the Department of Energy’s Office of
Energy Efficiency and Renewable Energy and located at the National Renewable
Energy Laboratory.Patrick Emami, Xiangyu Zhang, and Ahmed Zamzam are with the
National Renewable Energy Laboratory, Golden, CO 80401 USA (e-mail:
{patrick.emami, xiangyu.zhang<EMAIL_ADDRESS>Biagioni is with
Maplewell Energy, Broomfield, CO 80021, USA (email: dave@maplewelleng.com).

###### Abstract

In multi-timescale multi-agent reinforcement learning (MARL), agents interact
across different timescales. In general, policies for time-dependent
behaviors, such as those induced by multiple timescales, are non-stationary.
Learning non-stationary policies is challenging and typically requires
sophisticated or inefficient algorithms. Motivated by the prevalence of this
control problem in real-world complex systems, we introduce a simple framework
for learning non-stationary policies for multi-timescale MARL. Our approach
uses available information about agent timescales to define a periodic time
encoding. In detail, we theoretically demonstrate that the effects of non-
stationarity introduced by multiple timescales can be learned by a periodic
multi-agent policy. To learn such policies, we propose a policy gradient
algorithm that parameterizes the actor and critic with phase-functioned neural
networks, which provide an inductive bias for periodicity. The framework’s
ability to effectively learn multi-timescale policies is validated on a
gridworld and building energy management environment.

##  I INTRODUCTION

The ability to control multiple interacting components is essential to
efficiently manage complex systems. For instance, in power systems
applications, flexible loads and distributed generation operate within this
complex system with coupling in the dynamical models, constraints, and
objectives. Similar multi-component control challenges appear in
transportation systems, robotics, etc. Thus, interest in data-driven
approaches that try to learn distributed control policies from experience has
grown recently. Multi-agent reinforcement learning (MARL) is a promising
agent-based sequential decision-making framework for learning complex
coordination strategies. Crucially, it does not depend on having access to a
model of the (often stochastic and nonlinear) environment dynamics.

However, applying MARL to real-world problems is challenging because agents
typically only receive noisy, partial observations of the environment state
and have limited communication with each other. Moreover, from the perspective
of each agent, the environment dynamics appear to shift over time as other
agents learn to adapt their behaviors. All of these factors contribute to the
extreme _non-stationarity_ present in MARL, which makes learning good agent
policies notoriously challenging [1, 2]. Despite the dedication of significant
effort to discover practical MARL algorithms capable of learning policies in
the face of non-stationarity [3, 4], there remains a gap between the synthetic
environments used for algorithm development and the real world.

This work studies an under-explored MARL setting inspired by real world
applications where agents need to coordinate _time-dependent_ actions across
different timescales. This type of time-dependent coordination, arises, for
example, in:

  * Ex 1:

Power systems applications where we wish to learn a coordination strategy
between electrical devices that can be controlled at different timescales
(e.g., energy storage units, solar panels, and thermostatically controlled
loads) connected to the same micro-grid (Fig. 1),

  * Ex 2:

Robotic control tasks where multiple heterogeneous robots try to collaborate
to execute a task (e.g., moving an object) while their actuators are
controlled at different frequencies.

Introducing a time dependency via multiple timescales to MARL adds
_additional_ sources of non-stationarity, which makes learning policies
particularly challenging.

Figure 1: In Multi-timescale MARL, agents that act on different timescales try
to learn to coordinate to achieve a goal, e.g., controlling a heterogeneous
set of electronic devices for building energy optimization. Learning a _non-
stationary_ multi-agent policy allows multi-timescale agents to perform
complex time-dependent behaviors.

In this work, we formally define the multi-timescale decentralized partially
observable Markov decision process (Dec-POMDP) setting and propose a framework
for learning multi-agent time-dependent (i.e., non-stationary) policies to
solve such environments. We show that multiple timescales induces the optimal
multi-timescale policy to be periodic in nature. A practical policy gradient
method for learning periodic multi-agent policies based on phase-functioned
neural networks [5] is provided. Our framework’s ability to learn effective
multi-timescale policies with fewer environment interactions than key
baselines is validated on a gridworld and a building energy management
environment.

##  II Multi-timescale Dec-POMDPs

First, we define a Dec-POMDP [6], which is a multi-agent stochastic game
defined by a tuple:

$\texttt{DEC-
POMDP}:=(N,S,\mathbf{A},\mathbf{O},r,T,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}dt},\gamma,\rho,p,P).$
(1)

The number of agents is $N$, $S$ is the state space of the environment,
$\mathbf{A}:=\times_{i=1}^{N}A^{i}$ is a joint action space,
$\mathbf{O}:=\times_{i=1}^{N}O^{i}$ is a joint observation space,
$r(s,a):S\times\mathbf{A}\rightarrow\mathbb{R}$ is a global reward function,
$T$ is the (possibly infinite) horizon, $dt$ is the time discretization of the
environment, $\gamma\in[0,1]$ is a discount factor, and
$\rho:S\rightarrow[0,1]$ is the initial state distribution. The joint action
$\mathbf{a}\in\mathbf{A}$ induces a transition from state $s$ to state
$s^{\prime}$ according to the transition function
$p(s,\mathbf{a}):S\times\mathbf{A}\rightarrow S$. We assume that each agent is
only given a noisy, partial observation of the state governed by an
observation function $P(s,i):S\times\mathbb{N}\rightarrow O^{i}$, where $i$ is
the agent index, and that there is limited or no communication between agents.
Each agent learns a stochastic policy function $\pi^{i}(a^{i}\mid o^{i})$ that
maps an observation to a distribution over $A^{i}$, and the joint policy is
$\bm{\pi}=\\{\pi^{i},\dots,\pi^{N}\\}$.

In this work, we introduce multi-timescale Dec-POMDPs as an extension of Dec-
POMDPs:

$\texttt{MT-DEC-POMDP}:=[\texttt{DEC-POMDP};k,C],$ (2)

where $[;]$ is concatenation. The agents are defined with action frequencies
$k:=\\{k_{1},\dots,k_{N}\\}$ and the environment is defined to be periodic
with period $C\geq 1$, i.e., the reward function $r_{t}(s,\mathbf{a})$ and
transition functions $p_{t}(s,\mathbf{a})$ are $C$-periodic. When $C=1$, the
environment is aperiodic, as is typically assumed in a DEC-POMDP. Each agent’s
timescale is defined by the action frequency $k_{i}$ times the base timescale
discretization $dt$. For example, agent 1 with $k_{1}=2$ acts every two steps
(i.e., when $t\bmod 2dt=0$) and agent 2 with $k_{2}=3$ acts every three steps
(i.e., when $t\bmod 3dt=0$). We assume environments are defined with $dt=1$
for simplicity in the remainder of the paper. Between actions, “slow” agents
take a null action $a^{i}_{\texttt{null}}$ or repeat the most recently taken
action (e.g., when an action represents the setpoint for a device). Our
setting differs from the multi-timescale MARL setting considered in Wu et al.
[7] as they assume agents can communicate locally with their neighbors and
that the environment is aperiodic.

One can observe that the sequencing of actions across agent timescales repeats
periodically every $\tilde{K}=\text{LCM}(k_{1},\dots,k_{N})$ time steps, where
LCM is the least common multiple. Taking into account the periodicity of the
environment $C$, we see that the pattern of action and state transition
sequencing repeats every $K=\text{LCM}(\tilde{K},C)$ steps. When $K>1$,
solving a MT-DEC-POMDP is generally more difficult than solving a DEC-POMDP,
as the periodicity introduces a time dependency that compounds on top of the
non-stationarity caused by partial observability and limited communication
between agents [8, 3].

To help motivate our framework, we now briefly describe a class of MARL
environments where the optimal stationary joint policy is sub-optimal (i.e.,
when the agents ignore time dependencies). In detail, certain MARL
environments cause agents to suffer from what we call the _observation
aliasing problem_. In Sec. V-A, we provide an example multi-timescale
environment exhibiting observation aliasing. This aliasing problem, which is
closely related to the state aliasing problem that can arise in single agent
RL [9], occurs when agent $i$ receives a specific partial observation $o^{i}$
that has a different action under non-stationary policy $\pi^{i}$ depending on
the time step:

$\exists t,t^{\prime}\text{ with }t<t^{\prime}\text{ s.t.
}\pi_{t}^{i}(o^{i}_{t})\neq\pi_{t^{\prime}}^{i}(o^{i}_{t^{\prime}}).$

A _time-unaware_ agent is not able to learn a stationary policy that
distinguishes between $o^{i}_{t}$ and $o^{i}_{t^{\prime}}$ due to the
impossibility of proper credit assignment. Intuitively, time-unaware agents
perceive time-dependent rewards or dynamics as stochasticity in the
environment (e.g., from an exogenous source). Consequentially, these agents
learn suboptimal policies that try to explain the aliased observation as
perceived randomness.

A common heuristic is to make each agent time aware by appending the current
time step $t$ to the observation. Alternatively, recurrent neural networks can
be used to infer a _belief state_ for each agent based on its history [10]. In
this work, we explore how information about agent timescales and environment
periodicity can be used to more effectively learn non-stationary policies in
MT-DEC-POMDPs.

Recent work has introduced custom multi-timescale solutions for power systems
problems that simply attempt to learn a stationary policy [11] or use
recurrent neural networks to model temporal information for a fast agent and a
slow agent [12]. The single agent RL setting with multiple action frequencies
is also a closely related setting. Multi-timescale MDPs (MMDPs) [13, 14]
involve an agent formulated as a hierarchical MDP with fast and slow actions.
MMDPs aim to learn a (top-down) hierarchy over action timescales, in the sense
that actions taken on slower timescales influence the actions taken on faster
timescales, but not vice versa. Also related is a setting where a factored-
action MDP agent can choose actions that persist for various lengths of time
[15].

##  III $K$-periodic Non-stationary Joint Policies

We now demonstrate that, under simplifying assumptions, policy iteration in
the space of _$K$ -periodic_ non-stationary joint policies converges to the
optimal multi-timescale policy and value function.

###### Definition 1

For any $K\geq 1$, a $K$-periodic non-stationary joint policy satisfies

$\bm{\bar{\pi}}(\mathbf{o}_{t})=\bm{\bar{\pi}}(\mathbf{o}_{t+{K}}).$ (3)

Let
$\bm{\bar{\Pi}}_{K}=\big{\\{}\bm{\bar{\pi}}:\bm{\bar{\pi}}(\mathbf{o}_{t})=\bm{\bar{\pi}}(\mathbf{o}_{t+K})\big{\\}}$
be the set of all such policies. We use this to define a _multi-timescale non-
stationary joint policy_ $\bm{\pi}_{t}\in\bm{\Pi}$ with timescale action
frequencies $k$ as the policy _induced_ by a $K$-periodic non-stationary
policy:

$\displaystyle\forall t,\bm{\pi}_{t}=\\{\pi_{t}^{1},\dots,\pi_{t}^{N}\\}\text{
s.t. }$ (4) $\displaystyle\pi_{t}^{i}:=\begin{cases}\bar{\pi}_{t}^{i}\hskip
5.0pt&\text{if }t\bmod k_{i}=0\\\ \delta_{t-(t\bmod k_{i})}(a^{i})\text{ or
}\delta_{t}(a^{i}_{\texttt{null}})\hskip 5.0pt&\text{otherwise.}\end{cases}$

Here, $\delta_{t}$ is the Dirac delta function. In general, the inducing
policy for a multi-timescale policy $\bm{\pi}_{t}$ does not need to be
$K$-periodic. However, later in this section we prove that policy iteration in
the space of $K$-periodic policies converges to the optimal multi-timescale
policy assuming cooperative agents and full observability.

###### Definition 2

The projection operator for joint actions $\Gamma_{t}^{k}(\mathbf{a}_{t})$ is
defined as

$\displaystyle\Gamma_{t}^{k}(\mathbf{a}_{t})=\\{\bar{a}^{1}_{t},\dots,\bar{a}^{N}_{t}\\}\text{,
where }$ (5) $\displaystyle\bar{a}^{i}_{t}:=\begin{cases}a^{i}_{t}\hskip
5.0pt&\text{if }t\bmod k_{i}=0\\\ \bar{a}^{i}_{t-1}\text{ or
}a^{i}_{\texttt{null}}\hskip 5.0pt&\text{otherwise.}\end{cases}$ (6)

Notice that $\Gamma_{t}^{k}$ is $\tilde{K}$-periodic. That is,
$\Gamma_{t}^{k}(\mathbf{a})=\Gamma_{t+\tilde{K}}^{k}(\mathbf{a})$ for any $t$
and joint action $\mathbf{a}$.

###### Definition 3

The projection operator for the reward and transition function is defined as

$\Theta_{t}^{C}(f_{t})=\begin{cases}f_{t}^{0}\hskip 5.0pt&\text{if }t\bmod
C=0\\\ \vdots\\\ f_{t}^{C-1}\hskip 5.0pt&\text{if }t\bmod C=C-1.\end{cases}$
(7)

Notice that $\Theta_{t}^{C}$ is $C$-periodic, i.e.,
$\Theta_{t}^{C}(r(s,\mathbf{a}))=\Theta_{t+C}^{C}(r(s,\mathbf{a}))$ for any
$t$ and $(s,\mathbf{a})$.

The agents aim to learn an optimal multi-timescale policy

$\bm{\pi}^{*}=\operatorname*{arg\,max}_{\bm{\pi}\in\Pi}\mathbb{E}_{\bm{\pi},\Theta_{t}^{C}(p_{t}(s,a))}\Biggl{[}\sum_{t=0}^{T}\gamma^{t}\Theta_{t}^{C}\bigl{(}r_{t}(s_{t},\Gamma_{t}^{k}(\mathbf{a}_{t}))\bigr{)}\Biggr{]},$
(8)

and the optimal multi-timescale value function is

$Q_{t}^{\bm{\pi}^{*}}(a,s)=\mathbb{E}_{\mathbf{a_{t+1}}:\mathbf{a_{T}},s_{t+1}:s_{T}}\bigl{[}R_{t}|\Gamma_{t}^{k}(\mathbf{a_{t}}),s_{t},\bm{\pi}^{*}\bigr{]},$
(9)

where
$R_{t}=\sum_{\tau=t}^{T}\gamma^{\tau-t}\Theta^{C}_{\tau}r_{\tau}(s_{\tau},\mathbf{a}_{\tau})$.
The following theorem establishes convergence to the optimal multi-timescale
value function via _policy iteration in the space of $K$-periodic non-
stationary policies_ assuming agent cooperation and full observability.

###### Assumption 1

The agents are cooperative and have full observability of the environment.

###### Theorem 1

Under assumption 1, for any $\bm{\pi}^{0}\in\Pi$ induced by $K$-periodic non-
stationary policy $\bar{\bm{\pi}}^{0}\in\bar{\Pi}_{K}$,
$K=\text{LCM}(\tilde{K},C)$, the sequence of value functions
$Q^{\bm{\pi}^{n}}$ and improved policies $\bm{\pi}^{n+1}$ due to policy
iteration converges to the optimal multi-timescale value function and optimal
multi-timescale policy $\bm{\pi}^{*}$, i.e.,
$Q^{\bm{\pi}^{*}}(s,a)=\lim_{n\rightarrow\infty}Q^{\bm{\pi}^{n}}(s,a)\geq
Q^{\bm{\pi}}(s,a)$. Moreover, the optimal multi-timescale value function
always exists and is induced by a $K$-periodic non-stationary policy.

###### Lemma 1

Under assumption 1, MT-DEC-POMDP reduces to a multi-timescale multi-agent MDP
[16]. Furthermore, a multi-timescale multi-agent MDP is equivalent to an
action-persistent factored action (FA) MDP [15].

Proof  A multi-agent MDP can be described by the tuple
$\\{S,N,\mathbf{A},p,r\\}$ with elements defined as in Sec. II [16]. A multi-
timescale multi-agent MDP is defined similarly to MT-DEC-POMDP (Eq. (2)),
i.e., by expanding a multi-agent MDP with $k$ and $C$. Assuming agent
cooperation, MT-DEC-POMDP has a single reward shared by all agents. Assuming
full observability, the observation function $P(s,i)$ in a MT-DEC-POMDP is the
identity mapping. Therefore, each agent’s policy in MT-DEC-POMDP becomes
$\pi^{i}(a^{i}|s)$, which completes the reduction. A factored action (FA) MDP
with a fully factorized policy over $N$-dimensional actions
$\pi(a^{1},\dots,a^{N}|s)=\prod_{n=1}^{N}\pi^{n}(a^{n}|s)$ is a single agent
MDP that can be thought of as equivalently having $N$ agents, i.e., a multi-
agent MDP. A $k$-persistent FA MDP assumes that action $a^{i}$ is decided
every $k_{i}$ steps and repeated otherwise [15], and can similarly be extended
to periodic environments with period $C$. By modifying the persistence
property to allow action $a^{i}$ to be repeated $k_{i}$ times _or_ for a null
action $a^{i}_{\texttt{null}}$ to be subsequently taken $k_{i}-1$ times, the
$k$-persistent FA MDP is equivalent to the multi-timescale multi-agent MDP
setting, completing the proof.

Given Lemma 1, our proof for Theorem 1 follows the same proof technique used
to prove Theorem 3 in Lee et al. [15]. Differently, our proof handles
$C$-periodic environments. To that end, we define here the one-step multi-
timescale Bellman _optimality_ operator $\bar{\mathcal{T}}_{t}^{*}$ induced by
$\bm{\bar{\pi}}$ for $t\in\\{0,\dots,K-1\\}$:

$\displaystyle(\bar{\mathcal{T}}_{t}^{*}Q)(s,\mathbf{a}):=$ (10)
$\displaystyle\Theta_{t}^{C}(r_{t}(s_{t},\mathbf{a}_{t}))+\gamma\mathbb{E}_{s_{t+1}\sim\Theta_{t}^{C}(p_{t})}\bigl{[}\max_{\mathbf{a}_{t+1}}Q(s_{t+1},\Gamma_{t}^{k}(\mathbf{a}_{t+1})\bigr{]}.$

Notice that $\bar{\mathcal{T}}_{t}^{*}$ is $K$-periodic due to the
$\tilde{K}$-periodic action projection $\Gamma_{t}^{k}$ and the $C$-periodic
projection operator $\Theta_{t}^{C}$. Thus,
$\bar{\mathcal{T}}_{t}^{*}Q=\bar{\mathcal{T}}_{t+K}^{*}Q$ for any $t$ and $Q$.
Now, we define the $K$-step multi-timescale Bellman optimality operator
$\bar{H}_{t}^{*}$ by composing one-step Bellman optimality operators as
follows:

$\displaystyle(\bar{H}_{0}^{*}Q)(s,a)$
$\displaystyle:=(\bar{\mathcal{T}}_{0}^{*}\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-2}^{*}\bar{\mathcal{T}}_{K-1}^{*}Q)(s,a)$
(11) $\displaystyle(\bar{H}_{1}^{*}Q)(s,a)$
$\displaystyle:=(\bar{\mathcal{T}}_{1}^{*}\bar{\mathcal{T}}_{2}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}\bar{\mathcal{T}}_{K}^{*}Q)(s,a)$
$\displaystyle\phantom{:}\vdots$ $\displaystyle(\bar{H}_{K-1}^{*}Q)(s,a)$
$\displaystyle:=(\bar{\mathcal{T}}_{K-1}^{*}\bar{\mathcal{T}}_{K}^{*}\cdots\bar{\mathcal{T}}_{K-3}^{*}\bar{\mathcal{T}}_{K-2}^{*}Q)(s,a).$

The next lemma establishes that $K$-step multi-timescale Bellman optimality
operators are a contraction mapping.

###### Lemma 2

For all $t\in\\{0,\dots,K-1\\}$, the $K$-step multi-timescale Bellman
optimality operator $\bar{H}_{t}^{*}$ is a $\gamma^{K}$-contraction with
respect to infinity norm with $\bar{H}_{t}^{*}Q_{t}^{*}=Q_{t}^{*}$ as the
unique fixed point solution. That is, for any $Q_{t}^{0}$, define
$Q_{t}^{n+1}=\bar{H}_{t}^{*}Q_{t}^{n}$. Then, the sequence $Q_{t}^{n}$
converges to the $t\textsuperscript{th}$ multi-timescale optimal value
function as $n\rightarrow\infty$.

Proof  Without loss of generality, it is sufficient to prove the case $t=0$.

For any $Q_{1},Q_{2}$ and $s_{0}\in\mathcal{S},\mathbf{a}_{0}\in\mathbf{A}$,

$\displaystyle\phantom{=}\lvert(\bar{H}_{0}^{*}Q_{1})(s_{0},\mathbf{a}_{0})-(\bar{H}_{0}^{*}Q_{2}(s_{0},\mathbf{a}_{0})\rvert$
$\displaystyle=\lvert(\bar{\mathcal{T}}_{0}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s_{0},\mathbf{a}_{0})-(\bar{\mathcal{T}}_{0}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s_{0},\mathbf{a}_{0})\rvert$
$\displaystyle=\Biggl{\lvert}\mathbb{E}_{s_{1}\sim\Theta_{0}^{C}(p_{0}(s_{0},\mathbf{a}_{0}))}\Bigl{[}\Theta_{0}^{C}(r_{0}(s_{0},\mathbf{a}_{0}))$
$\displaystyle\phantom{=}+\gamma\max_{\mathbf{a}_{1}}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s_{1},\Gamma_{1,\mathbf{a}_{0}}^{k}(\mathbf{a}_{1}))\Bigr{]}$
$\displaystyle\phantom{=}-\mathbb{E}_{s_{1}\sim\Theta_{0}^{C}(p_{0}(s_{0},\mathbf{a}_{0}))}\Bigl{[}\Theta_{0}^{C}(r_{0}(s_{0},\mathbf{a}_{0}))$
$\displaystyle\phantom{=}+\gamma\max_{\mathbf{a}_{1}}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s_{1},\Gamma_{1,\mathbf{a}_{0}}^{k}(\mathbf{a}_{1}))\Bigr{]}\Biggr{\rvert}$
$\displaystyle=\gamma\Biggl{\lvert}\mathbb{E}_{\Theta_{0}^{C}(p_{0})}\Bigl{[}\max_{\mathbf{a}_{1}}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s_{1},\Gamma_{1,\mathbf{a}_{0}}^{k}(\mathbf{a}_{1}))$
$\displaystyle\phantom{=}-\max_{\mathbf{a}_{1}}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s_{1},\Gamma_{1,\mathbf{a}_{0}}^{k}(\mathbf{a}_{1}))\Bigr{]}\Biggr{\rvert}$
$\displaystyle\leq\gamma\Biggl{\lvert}\mathbb{E}_{\Theta_{0}^{C}(p_{0})}\Bigl{[}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s_{1},\mathbf{a}_{1}^{*})$
$\displaystyle\phantom{\leq}-(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s_{1},\mathbf{a}_{1}^{*})\Bigr{]}\Biggr{\rvert}$
$\displaystyle\phantom{\leq}\text{ where
}\mathbf{a}_{1}^{*}=\operatorname*{arg\,max}_{\mathbf{a}}\Bigl{[}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s_{1},\Gamma_{1,\mathbf{a}_{0}}^{k}(\mathbf{a}_{1}))$
$\displaystyle\phantom{\leq}-(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s_{1},\Gamma_{1,\mathbf{a}_{0}}^{k}(\mathbf{a}_{1}))\Bigr{]}$
$\displaystyle\leq\gamma\max_{s,\mathbf{a}}\Biggl{\lvert}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s,\mathbf{a})$
$\displaystyle\phantom{\leq}-(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s,\mathbf{a})\Biggr{\rvert}.$

We can continue to expand the inequality in a similar manner:

$\displaystyle\forall s_{0},\mathbf{a}_{0},$
$\displaystyle\phantom{\leq}\lvert(\bar{H}_{0}^{*}Q_{1})(s_{0},\mathbf{a}_{0})-(\bar{H}_{0}^{*}Q_{2})(s_{0},\mathbf{a}_{0})\rvert$
$\displaystyle\leq\gamma\max_{s,\mathbf{a}}\bigl{\lvert}(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s,\mathbf{a})-(\bar{\mathcal{T}}_{1}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s,\mathbf{a})\bigr{\rvert}$
$\displaystyle\leq\gamma^{2}\max_{s,\mathbf{a}}\bigl{\lvert}(\bar{\mathcal{T}}_{2}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{1})(s,\mathbf{a})-(\bar{\mathcal{T}}_{2}^{*}\cdots\bar{\mathcal{T}}_{K-1}^{*}Q_{2})(s,\mathbf{a})\bigr{\rvert}$
$\displaystyle\phantom{\leq}\vdots$
$\displaystyle\leq\gamma^{K}\max_{s,\mathbf{a}}\bigl{\lvert}Q_{1}(s,\mathbf{a})-Q_{2}(s,\mathbf{a})\bigr{\rvert},$

which implies
$||\bar{H}_{0}^{*}Q_{1}-\bar{H}_{0}^{*}Q_{2}||_{\infty}\leq\gamma^{K}||Q_{1}-Q_{2}||_{\infty}$.
Therefore $\bar{H}_{t}^{*}$ is a $\gamma^{K}$-contraction with respect to
infinity norm, and by the Banach fixed-point theorem,
$\bar{H}_{t}^{*}Q_{t}^{*}=Q_{t}^{*}$ is the unique fixed point solution for
all $t$.
It follows that the fixed points of $\bar{H}_{0}^{*},\dots,\bar{H}_{K-1}^{*}$
together make up the optimal multi-timescale value function, which is
represented by the $K$ values $Q_{0}^{\pi^{*}},\dots,Q_{K-1}^{\pi^{*}}$. Next,
we show that these fixed points have the largest value compared to any other
multi-timescale value function for any history-dependent policy
$\bar{\pi}\in\Pi$.

###### Lemma 3

Let $\bar{H}_{t\bmod K}^{*}=\bar{\mathcal{T}}_{t\bmod
K}^{*}\dots\bar{\mathcal{T}}_{(t+K-1)\bmod K}^{*}$ be the $K$-step multi-
timescale Bellman optimality operator and $Q^{\pi^{*}}_{t\bmod K}$ be its
fixed point. Then, for any history-dependent policy $\bar{\pi}\in\bar{\Pi}$,
$Q^{\pi^{*}}_{t\bmod K}(s,\mathbf{a})\geq Q_{t}^{\pi}(s,\mathbf{a})$.

Proof  For any $\bar{\pi}\in\bar{\Pi}$, $t,s,\mathbf{a}$ and $Q$, the
following inequality holds:

$\displaystyle\phantom{\leq}(\bar{\mathcal{T}}_{t}^{\bar{\pi}}Q)(s_{t},\mathbf{a}_{t})$
$\displaystyle:=\Theta^{C}_{t}(r_{t}(s_{t},\mathbf{a}_{t}))$
$\displaystyle\phantom{\leq}+\gamma\mathbb{E}_{s_{t+1}\sim\Theta^{C}_{t}(p_{t}),\mathbf{a}_{t+1}=\bar{\pi}}\bigl{[}Q(s_{t+1},\Gamma^{k}_{{t+1},\mathbf{a}_{t}}(\mathbf{a}_{t+1})\bigr{]}$
$\displaystyle\leq\Theta^{C}_{t}(r_{t}(s_{t},\mathbf{a}_{t}))$
$\displaystyle\phantom{\leq}+\gamma\max_{\mathbf{a}_{t+1}}\mathbb{E}_{s_{t+1}\sim\Theta^{C}_{t}(p_{t})}\bigl{[}Q(s_{t+1},\Gamma^{k}_{t+1}(\mathbf{a}_{t+1})\bigr{]}$
$\displaystyle=(\bar{\mathcal{T}}_{t\bmod K}^{*}Q)(s_{t},\mathbf{a}_{t}).$

This implies

$\displaystyle\phantom{\leq}(\bar{\mathcal{T}}_{t}^{\bar{\pi}}\bar{\mathcal{T}}_{t+1}^{\bar{\pi}}\dots\bar{\mathcal{T}}_{t+K-1}^{\bar{\pi}}Q)(s,\mathbf{a})$
$\displaystyle\leq(\bar{\mathcal{T}}_{t\bmod
K}^{*}\bar{\mathcal{T}}_{(t+1)\bmod K}^{*}\dots\bar{\mathcal{T}}_{(t+K-1)\bmod
K}^{*}Q)(s,\mathbf{a})$ $\displaystyle=(\bar{H}_{t\bmod
K}^{*}Q)(s,\mathbf{a}).$

Therefore,

$\displaystyle\phantom{=}Q^{\pi}_{t}(s,\mathbf{a})$
$\displaystyle=\lim_{n\rightarrow\infty}(\bar{\mathcal{T}}_{t}^{\bar{\pi}}\bar{\mathcal{T}}_{t+1}^{\bar{\pi}}\dots\bar{\mathcal{T}}_{t+Kn-1}^{\bar{\pi}}Q)(s,\mathbf{a})$
$\displaystyle\leq\lim_{n\rightarrow\infty}((\bar{H}_{t\bmod
K}^{*})^{n}Q)(s,\mathbf{a})=Q_{t}^{\pi^{*}}(s,\mathbf{a})$

holds, which concludes the proof.

Finally, to prove the main claim in Theorem 1, by Lemma 3 it is sufficient to
show $\lim_{n\rightarrow\infty}Q_{t}^{\pi^{n}}=Q_{t}^{\pi^{*}}$ for all
$t\in\\{0,\dots,K-1\\}$. To prove this, first, we must establish monotonic
improvement of multi-timescale policies induced by $\bar{\pi}^{n}$ during
policy iteration, i.e., $Q_{t}^{\pi^{n+1}}(s,\mathbf{a})\geq
Q_{t}^{\pi^{n}}(s,\mathbf{a})$ always holds for all $t,s,\mathbf{a},n$. This
result follows by proving a similar result to Theorem 2 in Lee et al. [15] in
the same manner as the proofs for Lemmas 2 and 3, i.e., by projecting the
reward and transition functions with $\Theta^{C}_{t}$. We refer to [15] for
the details. When the policy is no longer improving, i.e., $\pi^{n+1}=\pi^{n}$
and $Q_{t}^{\pi^{n+1}}=Q_{t}^{\pi^{n}}$, then by definition, $Q_{t}^{\pi^{n}}$
satisfies the multi-timescale Bellman optimality equation. By Lemma 2, the
multi-timescale Bellman optimality equation has a unique solution $Q_{t}^{*}$,
implying $Q_{t}^{*}=Q_{t}^{\pi^{n}}$. Thus, $\pi^{n}$ is the optimal multi-
timescale policy. Moreover, at every step of policy iteration, we have a
policy $\bar{\pi}^{n}$ which is in $\bar{\Pi}$ and which induces a multi-
timescale policy $\pi^{n}$ that is guaranteed to eventually converge to the
optimal multi-timescale policy. Therefore, the optimal multi-timescale policy
always exists.

The known inequality $Q^{*}_{\text{Dec-POMDP}}\leq Q^{*}_{\text{ MDP }}$
establishes that the optimal value function of a multi-timescale multi-agent
MDP is an upper bound of the optimal value function for MT-DEC-POMDP. This
implies that the optimal value function obtained by policy iteration under
Assumption 1 may overestimate the true optimal value function. See Theorem 5.1
in Oliehoek et al. [8] for proof details.

##  IV Phasic Policy Gradient Method

The theory from the previous section suggests encoding $K$-periodicity into
learning algorithms for MT-DEC-POMDP to encourage learning the optimal policy
when $k$ and $C$ are known. Let the integer
$\triangle_{t}\in\\{0,\dots,K-1\\}$ indicate the current _phase_ , i.e.,
$t\bmod K$. A straightforward approach is to encode this phase as a one-hot
vector of size $K$, $\texttt{O.H.}(\triangle)$, which can be concatenated to
each agent’s observation $[o^{i}_{t};\texttt{O.H.}(\triangle_{t})]$. However,
when the mapping between the phase-augmented observation and the optimal
action is complex, encoding the phase as a one-hot vector may not be
sufficient to provide a good inductive bias for $K$-periodicity.

Alternatively, we can parameterize each agent with phase-functioned neural
networks [5] (PFNNs). PFNNs are spline-based neural architectures whose
weights smoothly vary as a function of the current phase. This provides an
inductive bias of using similar weights for adjacent phases and reusing the
same network weights at time steps separated by a specified period. Favorably,
the number of parameters in PFNNs scales proportionally with the number of
spline control points (a constant) and _not_ with the period $K$. The use of
PFNNs in RL is under-explored, with only one known previous use for training
single agents in cyclic environments [17]. Each layer $l$ of a PFNN has a
weight matrix $\alpha$ computed by a _phase function_
$\alpha_{l}=\Theta(\beta_{l};{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2\pi\triangle_{t}/K)}$
conditioned on learnable weight matrices $\beta_{l}$ and phase
$2\pi\triangle_{t}/K\in[0,2\pi]$. Following [5], we use a Catmull-Rom spline
for $\Theta$, which is a cubic spline with four learnable spline control
points $\beta_{l}=[\beta_{l}^{0},\beta_{l}^{1},\beta_{l}^{2},\beta_{l}^{3}]$.
The weight for layer $l$ is

$\displaystyle\alpha_{l}$
$\displaystyle=\beta_{l}^{x_{1}}+w(\frac{1}{2}\beta_{l}^{x_{2}}-\frac{1}{2}\beta_{l}^{x_{0}})$
$\displaystyle+w^{2}(\beta_{l}^{x_{0}}-\frac{5}{2}\beta_{l}^{x_{1}}+2\beta_{l}^{x_{2}}-\frac{1}{2}\beta_{l}^{x_{3}})$
$\displaystyle+w^{3}(\frac{3}{2}\beta_{l}^{x_{1}}-\frac{3}{2}\beta_{l}^{x_{2}}+\frac{1}{2}\beta_{l}^{x_{3}}-\frac{1}{2}\beta_{l}^{x_{0}}),$

where $w=4\triangle_{t}/K\hskip 5.0pt(\texttt{mod }1)$ and $x_{n}=\lfloor
4\triangle_{t}/K\rfloor+n-1\hskip 5.0pt(\texttt{mod }4)$. The bias for layer
$l$ is computed in a similar fashion. The start and end control points for
each layer are the same, making each PFNN layer cyclic. In this work, we adapt
the actor-critic policy gradients method COMA [4] by using PFNNs with period
$K=\text{LCM}(\tilde{K},C)$ for the actor and critic networks.

##  V Experiments

###  V-A The Move Box Problem

Setup: We adapted a gridworld environment called Move Box [18] to create a toy
multi-timescale environment with a time-dependent optimal policy. That is, the
two agents need to coordinate their actions to push a green box to a goal
location (“G”) within a maximum of 20 steps. In the easy version, one agent
(red) is a “fast” agent that uses $k_{1}=1$ and one agent (blue) is a “slow”
agent that acts every two steps $k_{2}=2$. The period provided to time-aware
agents is therefore $K:=\texttt{LCM}(1,2,C=1)=2$. For the hard version of this
task, the fast agent uses $k_{1}=2$ and slow agent uses $k_{2}=3$, thus
$K:=\texttt{LCM}(2,3,C=1)=6$. Each agent receives a 4D partial observation
consisting of its own position and the position of the box. To avoid the need
for a complex exploration strategy, we restrict the action space to 3 discrete
actions: move up, move down, or null (do nothing). To push the green box up or
down, the agents have to be on either side of the box and move in the same
direction at the same time. Agents receive a reward of +1 whenever the box
moves towards the goal and a reward of +20 once the box reaches the goal. The
easy version can be solved with 7 actions while the hard version requires 19
actions; there is a significant increase in exploration difficulty between the
easy and hard versions. To implement multiple timescales, the only legal
action available to the slow agent between steps is the null action. _Move Box
is designed so that a time-unaware fast agent suffers from observation
aliasing_ (Sec. II). A time-unaware fast agent’s limited information means it
cannot precisely determine whether to move up or do nothing to synchronize its
actions with the slow agent.

Agents: We train four COMA-based agents:

  * •

Basic COMA, a time-unaware agent that uses a feedforward neural network for
the actor and critic networks.

  * •

Recurrent COMA, a time-aware agent that uses an RNN for the actor and critic
networks to condition on the full history up to the current time step [10].

  * •

One-Hot (O.H.) phase-aware COMA, an agent whose observations are augmented
with a one-hot encoding of the current phase $\triangle_{t}$.

  * •

Phasic COMA, an agent whose actor and critic networks are PFNNs with weights
indexed by $2\pi\triangle_{t}/K$.

We also run a variant of PFNN with period 4 instead of $K$, Phasic COMA (4),
to explore performance sensitivity to this hyperparameter. In the easy
environment $4>K$ and in the hard environment $4<K$. All actor and critic
networks share parameters. We take the standard approach of providing a one-
hot agent ID as an auxiliary input to distinguish between agents.

Figure 2: Move Box qualitative analysis. In the easy setting, the red agent
acts every step ($k_{1}=1$) and the blue agent acts every two steps
($k_{2}=2$). The joint action is shown in white on each agent. The time-
unaware Basic COMA red agent tries to move the box up at $t=2$, which causes
it to drop the box at $t=3$. The time-aware Phasic COMA red agent learns to
take a null action at $t=2$. The hard setting requires more sophisticated
coordination between agents. Best viewed in color.

Figure 3: Move Box return vs. steps. Mean return over 8 random seeds of the
learned greedy policy at various steps during training (shaded region is the
95% confidence interval (CI)). Best possible return shown by black dotted
line. The PFNN-based agent Phasic COMA with correct periods $K=2,6$ achieves
the highest reward in the fewest environment steps. TABLE I: Move Box results.
Success (%) is the fraction out of 8 random seeds that the agents learned to
move the box to the goal. The optimal return is 23.0.

Move Box [Easy] Move Box [Hard] COMA Success (%) Avg. return Success (%) Avg.
return Basic 0 1.0$\pm 0.0$ 0 1.0$\pm 0.0$ Recurrent 75 17.5$\pm 10.2$ 0
3.75$\pm 7.8$ O.H. phase-aware 100 23.0$\pm 0.0$ 0 3.45$\pm 7.9$ Phasic (4)
100 23.0$\pm 0.0$ 50 12.0$\pm 11.8$ Phasic (15) 100 23.0$\pm 0.0$ 50 12.0$\pm
11.8$

Results: Table I shows quantitative results and Fig. 3 compares test return as
a function of environment steps. Both O.H. phase-aware COMA and Phasic COMA
learn to reliably solve the easy Move Box environment across all random seeds,
with Phasic COMA demonstrating a small advantage in terms of efficiency.
Recurrent COMA needs more steps to achieve a good test return yet ultimately
performs less reliably. Basic COMA fails on this environment as expected due
to observation aliasing (Fig. 2). In the hard version (Fig. 3), only Phasic
COMA learns to solve the environment on just 50% of the training runs. The
PFNN-based actor and critic is robust to a slightly smaller or larger period
than $K$, although it appears to require more environment steps.

###  V-B Building Energy Management

Figure 4: BEM results. Mean return over 8 random seeds of the learned greedy
policy at various steps during training (shaded region is the 95% CI). Best
possible is the black dotted line. The Phasic COMA agent with correct period
$K=15$ outperforms non-phasic variants by a wide margin. Figure 5: BEM
qualitative results. Total power vs. power reference signal from one of the
Phasic COMA training runs. Figure 6: BEM qualitative results. Agent actions
from one of the Phasic COMA training runs.

Setup: In this environment, agents attempt to coordinate the control of HVAC
and energy storage (ES) for a five-zone small office building. See Biagioni et
al. [19] for details about the reduced order model used to simulate the
building. There are 7 agents: an HVAC agent per zone able to change the mass
flow rate (kg/s) every 5 mins, an HVAC chiller agent that can change the
discharge air temperature (°C) every 5 mins, and an ES agent that can change
its charging or discharging power (kW) every 1 minute. The goal is for the
agents to coordinate their total power consumption to track a reference power
signal that changes every 3 mins while minimizing discomfort to building
occupants. The control horizon is set to 30 mins ($dt=1$min). The global
reward function at time step $t$ is defined as

$r_{t}=-\sum_{\text{zone}_{i}}\big{(}(T_{t}^{i}-\overline{T})^{+}+(\underline{T}-T_{t}^{i})^{+}\big{)}^{2}-\alpha(p_{t}-p_{t}^{\text{ref}})^{2},$

where $\alpha=0.01$, $\overline{T}:=26$ °C and $\underline{T}:=24$ °C is the
thermal comfort band, $T_{t}^{i}$ is zone $i$’s temperature, and $p_{t}$ is
the total power. To implement multiple timescales, agents repeat their
previous action between steps. The period used for learning periodic non-
stationary policies is $K:=\texttt{LCM}(1,5,C=3)=15$, where $C$ encodes the
cyclic power reference signal.

Results: Out of all agents, only Phasic COMA is able to reliably learn a near-
optimal joint policy (Fig. 4). The variant with arbitrary PFNN period $4\ll K$
is the second best agent. We visualize the control actions selected by the
joint policy from one of the Phasic COMA runs in Fig. 6. The slow HVAC agents
have successfully learned to coordinate with the faster ES agent to track the
reference signal (Fig. 5) without violating the thermal comfort band; for
example, by increasing their power consumption between 10-15 mins while the ES
agent is already maximally discharging.

##  VI CONCLUSIONS

In this work, we proposed a multi-timescale MARL framework for learning
policies that can represent complex time-dependent behaviors. We introduced
the multi-timescale non-stationary joint policy as the policy induced by a
$K$-periodic non-stationary joint policy, where the period $K$ is given by
knowledge about agent timescales and cyclic environment components, both of
which are typically known _a priori_. We use phase-functioned neural networks
to introduce an inductive bias for learning a periodic non-stationary joint
policy. Our results on grid world and building energy management environments
establish the effectiveness of our framework, suggesting that follow-up work
could explore using it to solve more advanced power systems problems.

## References

  * [1] G. Papoudakis, F. Christianos, A. Rahman, and S. V. Albrecht, “Dealing with non-stationarity in multi-agent deep reinforcement learning,” _arXiv preprint arXiv:1906.04737_ , 2019.
  * [2] P. Hernandez-Leal, M. Kaisers, T. Baarslag, and E. M. de Cote, “A survey of learning in multiagent environments: Dealing with non-stationarity,” _ArXiv preprint_ , vol. abs/1707.09183, 2017.
  * [3] R. Lowe, Y. Wu, A. Tamar, J. Harb, P. Abbeel, and I. Mordatch, “Multi-agent actor-critic for mixed cooperative-competitive environments,” in _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , 2017, pp. 6379–6390.
  * [4] J. N. Foerster, G. Farquhar, T. Afouras, N. Nardelli, and S. Whiteson, “Counterfactual multi-agent policy gradients,” in _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), New Orleans, Louisiana, USA, February 2-7, 2018_ , 2018, pp. 2974–2982.
  * [5] D. Holden, T. Komura, and J. Saito, “Phase-functioned neural networks for character control,” _ACM Trans. Graph._ , vol. 36, no. 4, pp. 1–13, 2017\.
  * [6] D. S. Bernstein, R. Givan, N. Immerman, and S. Zilberstein, “The complexity of decentralized control of markov decision processes,” _Mathematics of operations research_ , vol. 27, no. 4, pp. 819–840, 2002.
  * [7] J. Wu, K. Li, and Q.-S. Jia, “Decentralized multi-agent reinforcement learning with multi-time scale of decision epochs,” in _2020 59th IEEE Conference on Decision and Control (CDC)_ , 2020, pp. 578–584, ISSN: 2576-2370.
  * [8] F. A. Oliehoek, M. T. J. Spaan, and N. Vlassis, “Optimal and approximate q-value functions for decentralized POMDPs,” _JAIR_ , vol. 32, pp. 289–353, 2008.
  * [9] F. Pardo, A. Tavakoli, V. Levdik, and P. Kormushev, “Time limits in reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 4045–4054.
  * [10] R. E. Wang, M. Everett, and J. P. How, “R-maddpg for partially observable environments and limited communication,” _arXiv preprint arXiv:2002.06684_ , 2020.
  * [11] Z. Wu, Y. Li, W. Gu, Z. Dong, J. Zhao, W. Liu, X.-P. Zhang, P. Liu, and Q. Sun, “Multi-timescale voltage control for distribution system based on multi-agent deep reinforcement learning,” _International Journal of Electrical Power & Energy Systems_, vol. 147, p. 108830, 2023.
  * [12] T. Ochoa, E. Gil, A. Angulo, and C. Valle, “Multi-agent deep reinforcement learning for efficient multi-timescale bidding of a hybrid power plant in day-ahead and real-time markets,” _Applied Energy_ , vol. 317, p. 119067, 2022.
  * [13] Hyeong Soo Chang, P. Fard, S. Marcus, and M. Shayman, “Multitime scale markov decision processes,” _IEEE Trans. Automat. Contr._ , vol. 48, no. 6, pp. 976–987, 2003.
  * [14] M. He, S. Murugesan, and J. Zhang, “A markov decision process approach to multi-timescale scheduling and pricing in smart grids with integrated wind generation,” in _2011 4th IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , 2011, pp. 125–128.
  * [15] J. Lee, B.-J. Lee, and K.-E. Kim, “Reinforcement learning for control with multiple frequencies,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 3254–3264, 2020.
  * [16] C. Boutilier, “Planning, learning and coordination in multiagent decision processes,” in _TARK_ , vol. 96. Citeseer, 1996, pp. 195–210.
  * [17] A. Sharma and K. M. Kitani, “Phase-parametric policies for reinforcement learning in cyclic environments,” in _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, (AAAI-18), the 30th innovative Applications of Artificial Intelligence (IAAI-18), and the 8th AAAI Symposium on Educational Advances in Artificial Intelligence (EAAI-18), New Orleans, Louisiana, USA, February 2-7, 2018_ , 2018, pp. 6540–6547.
  * [18] S. Jiang and C. Amato, “Multi-agent reinforcement learning with directed exploration and selective memory reuse,” in _Proceedings of the 36th Annual ACM Symposium on Applied Computing_ , 2021, pp. 777–784.
  * [19] D. Biagioni, X. Zhang, C. Adcock, M. Sinner, P. Graf, and J. King, “From model-based to model-free: Learning building control for demand response,” _arXiv preprint arXiv:2210.10203_ , 2022.