# Heat conduction in a three-dimensional momentum-conserving fluid

Rongxiang Luo<EMAIL_ADDRESS>Department of Physics, Fuzhou University,
Fuzhou 350108, Fujian, China  Fujian Science and Technology Innovation
Laboratory for Optoelectronic Information of China, Fuzhou 350108, Fujian,
China  Lisheng Huang  Department of Physics, Fuzhou University, Fuzhou 350108,
Fujian, China  Fujian Science and Technology Innovation Laboratory for
Optoelectronic Information of China, Fuzhou 350108, Fujian, China  Stefano
Lepri<EMAIL_ADDRESS>Consiglio Nazionale delle Ricerche, Istituto
dei Sistemi Complessi, via Madonna del Piano 10, I-50019 Sesto Fiorentino,
Italy  Istituto Nazionale di Fisica Nucleare, Sezione di Firenze, via G.
Sansone 1, I-50019 Sesto Fiorentino, Italy

###### Abstract

Size-dependence of energy transport and the effects of reduced dimensionality
on transport coefficients are of key importance for understanding
nonequilibrium properties of matter on the nanoscale. Here, we perform
nonequilibrium and equilibrium simulations of heat conduction in a 3D fluid
with the multiparticle collision dynamics, interacting with two thermal-walls.
We find that the bulk 3D momentum-conserving fluid has a finite non-diverging
thermal conductivity. However, for large aspect-ratios of the simulation box,
a crossover from 3D to one-dimensional (1D) abnormal behavior of the thermal
conductivity occurs. In this case, we demonstrate a transition from normal to
abnormal transport by a suitable decomposition of the energy current. These
results not only provide a direct verification of Fourier’s law but also
further confirm the validity of existing theories for 3D fluids. Moreover,
they indicate that abnormal heat transport persists also for almost 1D fluids
over a large range of sizes.

###### pacs:

Introduction. Understanding thermal transport is of key importance both for
basic science and for technological development S. Lepri2016 ; Benenti2020 ;
N. B. Li2012 ; Z. W. Zhang2020 ; MacDonald1978 . Macroscopic behavior of
thermal transport in the linear response regime is governed by Fourier’s law

$J=-\kappa\nabla T,$ (1)

where $J$ is the heat current, $\nabla T$ is the spatial temperature gradient,
and $\kappa$ is a finite constant termed as the thermal conductivity. The heat
conduction behavior is also known as “normal heat conduction” if it obeys the
Fourier law or “abnormal heat conduction” otherwise.

It is well established that for low-dimensional systems with a conserved total
momentum, heat conduction is generically abnormal S. Lepri2016 . Specifically,
for one-dimensional (1D) momentum-conserving systems, $\kappa$ usually
diverges with the system length $L$ as a power-law: $\kappa\sim L^{\alpha}$ S.
Lepri1997 . The exponent $\alpha=1/3$ is related to the fact that fluctuations
of the conserved field belong to the universality class of the famous Kardar-
Parisi-Zhang (KPZ) equation H. van Beijeren2012 ; Spohn2014 . Thus, the
phenomenon should be largely independent of the microscopic details. For
anharmonic chains this has been well numerically verified in several works
both for hard-core A. Dhar2001 ; G. Casati2003 ; Cipriani2005 ; S. Chen2014 ;
Hurtado2016 and smooth interaction potentials Das2014 ; Wang2016 . Actually, a
different $\alpha$ value signals a different universality class, as it may
occur at specific points of the phase diagram or for systems with a different
set of conserved quantities. A known instance is the case of nonlinear chains
with even potentials, where the value of $\alpha\neq 1/3$ G. R. Lee-
Dadswell2005 ; L. Delfini2006 ; Wang2011 although its precise value of
$\alpha$ is not yet fully assessed theoretically (see Ref.Benenti2020 and
references therein for a review).

For two-dimensional (2D) momentum-conserving systems, a $\kappa\sim\ln L$
divergence of the thermal conductivity is predicted from different theoretical
apporaches S. Lepri1998 ; O. Narayan2002 and also from an exactly solved
stochastic model G. Basile2006 . As for 1D case, such prediction has been
numerically verified in 2D lattices A. Lippi2000 ; Delfini2005 ; L. Yang2006 ;
D. Xiong2010 , and only recently, also in 2D fluids P. D. Cintio2017 ; R. X.
Luo2020 . However, there are some cases in which a power-law divergence have
been reported, see Wang2020 and references therein. Therefore, it is worth
noting that these theoretical predictions equally apply to both lattices and
fluids. With respect to the latter case, abnormal heat transport is related to
the famous long-time-tail problem, long known in statistical mechanics
Pomeau1975 .

For the three-dimensional (3D) momentum-conserving systems, which is the
subject of the present paper, all analytical predictions of existing theories
S. Lepri2016 ; G. Basile2006 support that the heat conduction behavior is
normal. In 2010, the first numerical evidence for the validity of Fourier’s
law was given by Saito and Dhar K. Saito2010 . They performed nonequilibrium
simulations of heat conduction in a 3D anharmonic lattice and found that the
system, under certain conditions, has a finite nondiverging $\kappa$ but can
exhibit a crossover from 3D to 1D behavior. In the same year, similar results
for the same class of lattice systems were obtained by Wang _et al_. by means
of equilibrium simulations L. Wang2010 . However, to the best of our
knowledge, at present only the 3D lattices have been examined. Thus, to verify
whether the above picture of low-dimensional heat conduction applies equally
well to both lattices and fluids, numerical evidences for the counterpart 3D
fluids are still in urgent need. This then raises the relevant question: Can
the finite thermal conductivity also be confirmed in 3D fluids? Also, can the
dimensional-crossover behavior of heat conduction be observed in 3D fluids?

In this work, we give a positive answer to the above two questions. To this
aim, we consider a 3D system of particles in a cuboid-shaped box, interacting
through the multi-particle collision (MPC) dynamics A. Malevanets1999 . The
MPC method consists in replacing the conventional deterministic molecular
dynamics with a stochastic rotation of velocities. It can correctly capture
the hydrodynamic equations with great advantages in the numerical simulations
J. T. Padding2006 . By using the MPC dynamics, researchers achieved a
considerable understanding of various aspects of transport G. Gompper2009 .
More recently, it has been adopted to study the coupled particle and heat
transport Benenti2014 ; R.X. Luo2018 . Importantly, the MPC dynamics conserves
the total momentum and energy of the system, and thus it can help us test the
theoretical conjecture for momentum-conserving systems P. D. Cintio2017 ; R.
X. Luo2020 .

Here, we study heat transport properties of the 3D MPC fluid by means of
extensive numerical simulations. We employ the nonequilibrium thermal-wall
method where the system interacts with two heat baths at different
temperatures and compare the results with those obtained by equilibrium Green-
Kubo method. We show that the 3D fluid has finite thermal conductivity. Upon
increasing the aspect ratio of the simulation box, the fluid system can also
exhibit the crossover from 3D to 1D behavior of the thermal conductivity.

The 3D fluid model. We consider a 3D system of $N$ interacting point particles
with equal mass $m$. All particles are confined in a cuboid-shaped volume of
length $L$, width $W$, and height $H$ in the $x$, $y$ and $z$ coordinate,
respectively. At $x=0$ and $x=L$ in longitudinal direction, the particles
exchange heat with two heat baths at temperatures $T_{h}$ and $T_{c}$, modeled
as thermal walls Lebowitz1978 . When a particle hits the $x=0$ ($x=L$)
boundary (of area $WH$), it is reflected back with new velocity components
(denoted by $v_{x}$, $v_{y}$ and $v_{z}$ of its $x$, $y$ and $z$ directions,
respectively) drawn from a distribution with probability densities
Lebowitz1978

$\displaystyle P\left(v_{x}\right)=$
$\displaystyle\frac{m|v_{x}|}{k_{B}T_{\iota}}\textrm{exp}\left(-\frac{mv^{2}_{x}}{2k_{B}T_{\iota}}\right),$
(2) $\displaystyle P(v_{y,z})=$ $\displaystyle\sqrt{\frac{m}{2\pi
k_{B}T_{\iota}}}\textrm{exp}\left(-\frac{mv^{2}_{y,z}}{2k_{B}T_{\iota}}\right),$

where $T_{\iota}$ ($\iota=h,c$) is the temperature of the respective heat bath
in dimensionless units and $k_{B}$ is the Boltzmann constant. In the $y$ and
$z$ directions the particles are subject to periodic boundary conditions (even
though, it has been verified numerically that the results also apply to fixed
boundary conditions, since in both cases $v_{y,z}>0$ and $v_{y,z}<0$ are with
the same probability.).

The evolution in the bulk occurs by MPC dynamics A. Malevanets1999 ; J. T.
Padding2006 ; G. Gompper2009 , which simplifies the numerical simulation of
interacting particles by coarse-graining the time and space at which
interaction occur. Specifically, by MPC the system evolves in discrete time
steps, consisting of non-interacting propagation during a time $\tau$ followed
by instantaneous collision events. During the propagation, a particle keeps
its velocity $\mathbf{v}_{i}$ unchanged and its position is updated as
$\mathbf{r}_{i}\rightarrow\mathbf{r}_{i}+\tau\mathbf{v}_{i}$. As for the
collisions, the system’s volume is partitioned into cubic cells of linear size
$a$ and for all particles in a cell, their velocities are rotated around a
randomly chosen axis, with respect to their center of mass velocity
$\mathbf{V}_{CM}$ by an angle, $\theta$ or $-\theta$, randomly chosen with
equal probability. The velocity of a particle in a cell is thus updated from
$\mathbf{v}_{i}$ to
$\mathbf{V}_{CM}+\hat{\mathcal{R}}^{\pm\theta}\left(\mathbf{v}_{i}-\mathbf{V}_{CM}\right)$,
where $\hat{\mathcal{R}}^{\theta}$ is the rotation operator by the angle
$\theta$. Such moves preserve the total momentum and energy of the fluid
system. The time interval between successive collisions $\tau$ and the
collision angle $\theta$ tune the strength of the interactions and
consequently affect the transport of the particles. Note that the angle
$\theta=\pi/2$ corresponds to the most efficient mixing of the particle
momenta.

In our simulations, we set $T_{\iota}$ to be slightly biased from the nominal
temperature $T$, i.e., $T_{h,c}=T\pm\Delta T/2$, to measure the temperature
profile $T(x)$, where $x$ is the space variable, and the heat current $J$.
Numerically, $T(x)$ and $J$ are obtained in a same way as in the 2D case R. X.
Luo2020 . The thermal conductivity is finally computed by assuming Fourier’s
law, as $\kappa\approx JL/(T_{h}-T_{c})$, since we have checked that in the
linear response regime, the temperature jump between the heat bath and
simulated system’s boundary can be neglected (as shown in Fig. 3 below). We
set the main parameters as follows: $m=T=k_{B}=1$, $\Delta T=0.2$, $a=0.1$,
$\theta=\pi/2$, $\tau=0.1$, and the averaged particle number density
$\rho=N/(LWH)=88$. Note that the average distance,
$\lambda=\tau\sqrt{k_{B}T/m}$, that particles stream between rotations is
equal to $a$, meaning that our simulations satisfy Galiean-invariance of the
stochastic rotation dynamics T. Ihle2001 . We set $W=H$ henceforth, and focus
only on the the heat conduction in the $x$ direction for different $W$ and
$L$. So, the relevant parameter is the cell aspect-ratio $W/L$. Relatively
small values of it correspond to quasi-1D setups. In addition, long enough
relaxation times $(>10^{7})$ and integration times $(>10^{8})$ are utilized to
ensure the system reaches the steady state and the relative errors of all
statistical averages are smaller than 0.5$\%$.

Numerical results. First of all, we provide evidence that in the ballistic
(integrable) limit the thermal conductivity is a linear function of the system
length. If all particles do not interact, i.e., each particle keeps its
velocity unchanged from one heat bath to the other as it crosses the system,
then the system is integrable. Under this situation, $J$ can be computed by a
similar analysis as done in R. X. Luo2019 . Taking it into Eq. (1), we can
obtain an analytical expression for the thermal conductivity:

$\kappa=\rho L\sqrt{\frac{8k^{3}_{B}}{\pi
m}}/\left(\frac{1}{\sqrt{T_{h}}}+\frac{1}{\sqrt{T_{c}}}\right).$ (3)

In Fig. 1 this result is compared with our simulations (black circles) and the
agreement is perfect.

Figure 1: (Color online) The thermal conductivity $\kappa$ as a function of
the system length $L$ for the 3D fluid system with different $W$ values. The
symbols are for the numerical results, and for reference the black solid line
indicates the divergence with $L$ as $\sim L^{1/3}$. The red dot-dashed line
is a fit of the data with $W=1$ with the functional expression
$\kappa_{kin}=652.25/(1+12/L)$.

Now, let us turn to the interacting systems with collisions. It can be seen in
Fig. 1 that for $W=1$, $\kappa$ tends to saturate and becomes constant as $L$
is increased. This result indicates that heat conduction for the 3D momentum-
conserving fluid system is governed by the Fourier law, as expected S.
Lepri2016 ; G. Basile2006 . However, it can also be seen that as $W$
decreases, $\kappa$ is no longer constant but diverges with $L$. In
particular, for $W=0.1$, $\kappa$ eventually approaches the scaling
$\kappa\sim L^{1/3}$ predicted for 1D momentum-conserving fluids S. Chen2014 ;
O. Narayan2002 . Not surprisingly, the convergence is relatively slow and is
affected by finite-size corrections, since anomalous scaling is expected to
set in only at large scales. Note also that the scaling holds on a range of
about two decades. We also remark that the dimensional crossover from small to
large aspect ratio is similar to what observed in 3D lattices K. Saito2010 ;
L. Wang2010 .

Figure 2: (Color online) Scaling of the abnormal component of the heat
currents $J_{A}$ for three $W$ values. $J_{A}$ is computed subtracting the
normal component $J_{N}$ as detailed in the text. The dashed line is a best-
fit of the data, yielding a decay $J_{A}\propto L^{-0.61}$, while the solid
line is the expected decay $L^{-2/3}$.

To further support the scenario of quasi-1D transport, we perform a similar
analysis as done in S. Lepri2020 . The idea is to decompose the heat current
in a “normal” and “abnormal” components, i.e., $J=J_{N}+J_{A}$. To estimate
the abnormal part, we first numerically compute the heat current $J$ and
subtract from it the best estimate of the “normal” heat current, obtained as
$J_{N}=\kappa_{kin}\Delta T/L$, where $\kappa_{kin}=\kappa_{0}/(1+A/L)$. Such
formula is given by kinetic theory to account for the initial ballistic
regime. The parameter $A$ is related to boundary resistance between the fluid
and the reservoirs K. Aoki2001 ; S. Lepri2003 , and$\kappa_{0}$ is the bulk
value of the conductivity. As shown in Fig.1 the fit of the data with $W=1$ is
very accurate meaning that the abnormal contribution to the thermal
conductivity is negligible over the considered length range. The data in Fig.
2 show that for all small $W$ values, $J_{A}\sim L^{-0.61}$ when
$L\rightarrow\infty$, approaching the expected scaling $L^{-2/3}$, the value
predicted for 1D systems with the momentum conservation S. Lepri2016 . We thus
conclude that the measurements in 3D fluid system confirm the normal heat
conduction and the crossover from 3D to 1D behavior of the thermal
conductivity.

Figure 3: (Color online) Plot of temperature profiles for the 3D fluid system
with different $L$ values. Here our numerical results are compared with the
analytical Eq. (4). In (a) and (b) we fix $W=1.0$ and $W=0.1$, respectively.
Inset: Plot of the differences $\delta T$ between the data and the black line.

The difference between normal and abnormal heat conduction behaviors can be
further appreciated also in the steady-state kinetic temperature profiles
$T(x)$. For systems with normal heat conduction $T(x)$ is determined by
solving the stationary heat equation assuming that the thermal conductivity is
proportional to $\sqrt{T}$ as prescribed by standard kinetic theory, yielding
A. Dhar2001

$T(x)=\left[T_{h}^{3/2}\left(1-\frac{x}{L}\right)+T_{c}^{3/2}\frac{x}{L}\right]^{2/3}.$
(4)

In Fig. 3(a) this prediction is compared with our simulation results for
$W=1$. It is seen that there is a good agreement between the results of our
numerical simulations and Eq. (4). To better appreciate the deviations from
the prediction, we also plot the differences $\delta T$ between the data and
the black line. It is shown in the inset of Fig. 3(a) that $|\delta T|$
decreases with increasing $L$, as expected since $|\nabla T|=\Delta T/L$
decreases when $L$ increases, indicating that the linear response can
correctly describe the transport properties of the system for large enough
system length. On the other hand, for systems with abnormal heat conduction
the temperature profile is expected to be qualitatively different, being
solution of a fractional diffusion equation as demonstrated in several cases
S. Lepri2011 ; Kundu2019 . A typical feature is that the temperature profile
is concave upwards in part of the system and concave downwards elsewhere, and
this is true even for small temperature differences T. Mai2007 ; A. Lippi2000
; S. Lepri2009 . This is confirmed in Fig. 3(b) by our numerical simulations
for $W=0.1$. Note that the data for three different $L$ overlap with each
other, implying that the deviations from Fourier’s behavior are not finite-
size effects. Altogether, those numerical results again support our findings
based on the length-dependence of the thermal conductivity, that heat
conduction in 3D is normal while in lower dimensions it is abnormal.

We now turn to the comparison with linear-response results. Based on the
celebrated Green-Kubo formula, which relates transport coefficients to the
current time-correlation functions, the thermal conductivity can be expressed
as R. Kubo1991 ; S. Lepri2003 ; A. Dhar2008

$\kappa_{GK}=\frac{\rho}{k_{B}T^{2}}\lim_{\tau_{\mathrm{tr}}\rightarrow\infty}\lim_{N\rightarrow\infty}\frac{1}{N}\int_{0}^{\tau_{\mathrm{tr}}}\left<J(0)J(t)\right>dt.$
(5)

In this formula, $J\equiv\frac{1}{2}\Sigma_{i}^{N}\mathbf{v}_{i}^{2}v_{x,i}$
represents the total heat current along the $x$ coordinate and
$\left<J(0)J(t)\right>$ is its correlation function in the equilibrium state.
In the simulations, we consider an isolated fluid with periodic boundary
conditions also in the $x$ direction. The initial condition is randomly
assigned with the constraints that the total momentum is zero and the total
energy corresponds to $T=1$. The system is then evolved and after the
equilibrium state is attained, we compute $\left<J(0)J(t)\right>/N$ and the
integral in Eq. (5).

Figure 4: (Color online) (a) Correlation functions of the total heat current
for the 3D momentum-conserving fluid system. The black solid line and the blue
dashed line are the fit of the data with $W=1.0$ and $W=0.1$, respectively.
(b) the comparison of the thermal conductivity obtained by using the Green-
Kubo formula [Eq. (5)] and by using the nonequilibrium setting (symbols). The
integration is up a time $L/C$ with $C=12.0$ fixed empirically (see text).

The results for different $W$ values are presented in Fig. 4. It can be seen
from Fig. 4(a) that for $W=1$, the correlation function undergoes a rapid
decay at short times, and eventually, when $t>20$, it begins to oscillate
around zero, (the negative values of $\left<J(0)J(t)\right>/N$ are not shown
in this log-log scale.). Fitting exhibits an exponential decay. However, it is
clear that for $W=0.1$, after a rapid decay at short times, the correlation
function eventually follows a power law decay of $\left<J(0)J(t)\right>/N\sim
t^{\gamma}$, fully compatible with the theoretical prediction $\gamma=-2/3$.

As a further check, in Fig. 4(b), we compare the nonequilibrium values of
$\kappa$ (symbols) with the finite-size Green-Kubo value $\kappa_{GK}(L)$
(lines). The latter is obtained truncating the integral in Eq.(5) up to
$\tau_{\mathrm{tr}}=L/C$ where $C$ is a constant that we fix empirically to
yield the best match between $\kappa$ and $\kappa_{GK}(L)$. This procedure is
justified by the fact that for finite systems the integral in Eq.(5) should be
extended up to times of order of the characteristic transit time of waves
$L/c_{s}$ ($c_{s}$ being the sound speed) S. Lepri2003 ; A. Dhar2008 . It can
be seen that $\kappa_{GK}(L)$ agrees with $\kappa$, especially the part for
large system length. Comparing the results for $W=1$ and $W=0.1$, we conclude
that for a fixed $L$ (see $L=4096$), a gradual transition from 3D to 1D
behavior of the thermal conductivity will take place when $W$ is decreased.
Note also that, both in the 3D and quasi-1D cases the conductivity should
approach the bulk value for $L$ larger than the mean-free-path which the same
in both cases. So it is reasonable that the two curves start to differ from a
size $L\approx 100$ much larger the mean free path for $\tau=0.1$. Once again,
the emerging picture supports the validity of the Fourier law for large $W$
and the dimensionality-crossover behavior for small $W$.

Summary and discussion. We have numerically studied heat conduction in a 3D
fluid by nonequilibrium thermal-wall and equilibrium Green-Kubo methods. Both
give sound evidence that in the 3D case i.e. for $W$ is large enough, the
system has a finite nondiverging thermal conductivity. As $W$ decreases, the
system exhibits a crossover from 3D to 1D behavior of the thermal
conductivity. These results are rooted in the condition for 3D fluid systems
described by MPC dynamics to keep the momentum conservation of the system. The
fact that the microscopic MPC dynamics is inherently stochastic does not alter
the scenario, as seen by comparing with anharmonic lattices with symplectic
dynamics K. Saito2010 . Also we provided a strong evidence that the Fourier
law is valid in 3D fluid systems, further confirming that the predictions of
existing theories S. Lepri2016 ; G. Basile2006 apply to fluids. In addition,
our results suggest that the dimensionality-crossover behavior of heat
conduction in 3D lattices K. Saito2010 ; L. Wang2010 can also take place in 3D
fluids. Remarkably, the quasi-1D fluid displays the same transition from a
normal to abnormal behavior for large enough longitudinal sizes,
characteristic of strictly 1D models (Zhao2018, ; Miron2019, ; S. Lepri2020, ;
Lepri2021, ). This stems from the fact that abnormal (hydrodynamic) effects
can be only observed on scales that can be much larger than the typical
kinetic mean-free-path (S. Lepri2020, ). Apart from general theoretical
implications in transport theory, our findings may have experimental relevance
as well, because in experimental setup researchers usually fix the cross-
sectional area of samples and study how the thermal conductivity changes with
the sample’s length for achieving the purpose of understanding and controlling
thermal transport.

Acknowledgments. We acknowledge support by the Education Department of Fujian
Province (Grant No. JAT200035).

## References

  * (1) S. Lepri (ed.), _Thermal Transport in Low Dimensions: From Statistical Physics to Nanoscale Heat Transfer_ (Springer, Berlin, 2016).
  * (2) G. Benenti, S. Lepri and R. Livi, Frontiers in Physics, Vol 8, p.292 (2020).
  * (3) N. B. Li, J. Ren, L. Wang, G. Zhang, P. Hanggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012).
  * (4) Z. W. Zhang, Y. L. Ouyang, Y. Cheng, J. Chen, N. B. Li, and G. Zhang, Phys. Rep. 860, 1 (2020).
  * (5) R. A. MacDonald and D. H. Tsai, Phys. Rep. 46, 1 (1978).
  * (6) S. Lepri, R. Livi, and A. Politi, Phys. Rev. Lett. 78, 1896 (1997).
  * (7) H. van Beijeren, Phys. Rev. Lett. 108, 180601 (2012).
  * (8) H. Spohn, J. Stat. Phys. 154, 1191 (2014).
  * (9) A. Dhar, Phys. Rev. Lett. 86, 3554 (2001).
  * (10) G. Casati and T. Prosen, Phys. Rev. E 67, 015203(R) (2003).
  * (11) P. Cipriani, S. Denisov, and A. Politi, Phys. Rev. Lett. 94, 244301 (2005).
  * (12) S. Chen, J. Wang, G. Casati, and G. Benenti, Phys. Rev. E 90, 032134 (2014).
  * (13) P. I. Hurtado and P. L. Garrido, Sci. Rep. 6, 38823 (2016).
  * (14) S. Das, A. Dhar, and O. Narayan, J. Stat. Phys. 154, 204 (2014).
  * (15) L. Wang, N. B. Li and P. Haenngi in S. Lepri (ed.), _Thermal Transport in Low Dimensions: From Statistical Physics to Nanoscale Heat Transfer_ (Springer, Berlin, 2016).
  * (16) G. R. Lee-Dadswell, B. G. Nickel, and C. G. Gray, Phys. Rev. E 72, 031202 (2005); J. Stat. Phys. 132, 1 (2008).
  * (17) L. Delfini, S. Lepri, R. Livi, and A. Politi, Phys. Rev. E 73, 060201(R) (2006); J. Stat. Mech. (2007), P02007.
  * (18) L. Wang and T. Wang, Europhys. Lett. 93, 54002 (2011).
  * (19) S. Lepri, R. Livi, and A. Politi, Europhys. Lett. 43, 271 (2002).
  * (20) O. Narayan and S. Ramaswamy, Phys. Rev. Lett. 89, 200601 (2002).
  * (21) G. Basile, C. Bernardin, and S. Olla, Phys. Rev. Lett. 96, 204303 (2006).
  * (22) A. Lippi and R. Livi, J. Stat. Phys. 100, 1147 (2000).
  * (23) L. Delfini, S. Lepri, R. Livi, J. Stat. Mech. 2005 (05), P05006.
  * (24) L. Yang, P. Grassberger, and B. Hu, Phys. Rev. E 74, 062101 (2006).
  * (25) D. Xiong, J. Wang, Y. Zhang, and H. Zhao, Phys. Rev. E 82, 030101(R) (2010).
  * (26) P. DiCintio, R. Livi, S. Lepri, and G. Ciraolo, Phys. Rev. E 95, 043203 (2017).
  * (27) R. X. Luo, Phys. Rev. E 102, 052104 (2020).
  * (28) J. Wang, T. X. Liu, X. Z. Luo, X. L. Xu, and N. B. Li, Phys. Rev. E 101, 012126 (2020).
  * (29) Y. Pomeau and P. Resibois, Phys. Rep. 19 (2), 63 (1975).
  * (30) K. Saito and A. Dhar, Phys. Rev. Lett. 104, 040601 (2010).
  * (31) L. Wang, D. He, and B. Hu, Phys. Rev. Lett. 105, 160601 (2010).
  * (32) A. Malevanets and R. Kapral, J. Chem. Phys. 110, 8605 (1999).
  * (33) J. T. Padding and A. A. Louis, Phys. Rev. E 74, 031402 (2006).
  * (34) G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler, _Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids_ (Springer, Berlin 2009).
  * (35) G. Benenti, G. Casati, C. Mejía-Monasterio, New J. Phys. 16, 015014 (2014).
  * (36) R. X. Luo, G. Benenti, G. Casati, and J. Wang, Phys. Rev. Lett. 121, 080602 (2018); Phys. Rev. Research 2, 022009(R) (2020).
  * (37) J. L. Lebowitz and H. Spohn, J. Stat. Phys. 19, 633 (1978); R. Tehver, F. Toigo, J. Koplik, and J. R. Banavar, Phys. Rev. E  57, R17(R) (1998).
  * (38) T. Ihle and D. M. Kroll, Phys. Rev. E 63, 020201(R) (2001).
  * (39) R. X. Luo, Phys. Rev. E 99, 032138 (2019).
  * (40) S. Lepri, R. Livi, and A. Politi, Phys. Rev. Lett. 125, 040604 (2020).
  * (41) K. Aoki and D. Kusnezov, Phys. Rev. Lett. 86, 4029 (2001).
  * (42) S. Lepri and A. Politi, Phys. Rev. E 83, 030107(R) (2011).
  * (43) A. Kundu, C. Bernardin, K. Saito, A. Kundu, A Dhar, J. Stat. Mech. 2019 (1), 013205.
  * (44) T. Mai, A. Dhar, and O. Narayan, Phys. Rev. Lett. 98, 184301 (2007).
  * (45) S. Lepri, C. Mejía-Monasterio, and A. Politi, J. Phys. A 42, 025001 (2009).
  * (46) R. Kubo, M. Toda, and N. Hashitsume, _Statistical Physics II: Nonequilibrium Statistical Mechanics_ , (Springer, New York, 1991).
  * (47) S. Lepri, R. Livi, and A. Politi, Phys. Rep. 377, 1 (2003).
  * (48) A. Dhar, Adv. Phys. 57, 457 (2008).
  * (49) H. Zhao and W. G. Wang, Phys. Rev. E  97, 010103(R) (2018).
  * (50) A. Miron, J. Cividini, A. Kundu, and D. Mukamel, Phys. Rev. E  99, 012124 (2019).
  * (51) S. Lepri, G. Ciraolo, P. Di. Cintio, J. Gunn, and R. Livi, Phys. Rev. Research 3, 013207 (2021).


11institutetext: Huawei Ireland Research Center
11email<EMAIL_ADDRESS>
Authors with equal contribution

# Continual Model-based Reinforcement Learning for Data Efficient Wireless
Network Optimisation

Cengis Hasan  Alexandros Agapitos  David Lynch  Alberto Castagna  Giorgio
Cruciata  Hao Wang  Aleksandar Milenovic

###### Abstract

We present a method that addresses the pain point of long lead-time required
to deploy cell-level parameter optimisation policies to new wireless network
sites. Given a sequence of action spaces represented by overlapping subsets of
cell-level configuration parameters provided by domain experts, we formulate
throughput optimisation as Continual Reinforcement Learning of control
policies. Simulation results suggest that the proposed system is able to
shorten the end-to-end deployment lead-time by two-fold compared to a
reinitialise-and-retrain baseline without any drop in optimisation gain.

##  1 Introduction

One of the major factors influencing the Quality of Experience (QoE) in
wireless networks is the parameter configuration of the cells in a base-
station. Incorrectly configured cells can interfere with neighbouring cells
and degrade quality of service through inadequate coverage or over-
utilisation. Traditionally, cell-level parameter configuration is realised
before deployment, at which time the engineers have to anticipate diverse
traffic conditions (i.e. user load), radio channel conditions, environment
conditions (i.e. physical location and surroundings of a cell), and the
complex relationship of QoE objective with other conflicting optimisation
objectives in a wireless network (i.e. coverage, utilisation, power
consumption).

Existing approaches for configuring cell-level parameters are based on a
discrete black-box optimisation problem formulation, which is often solved
with surrogate-driven optimisation methods. Due to the combinatorial parameter
space defined by the large number of cell-level configuration parameters (>
500 parameters), and the variable impact of the same set of parameters on the
QoE in different deployment sites, selecting and configuring parameters for a
new system deployment is a stage-wise process. In this process, priorities are
assigned to subsets of parameters based on domain knowledge, and these subsets
are then configured in order of priority in a sequence of independent
surrogate-driven optimisation runs. When QoS metrics reach target values, the
parameter selection and configuration process is terminated and the parameter
configuration is executed to each cell. Each independent optimisation run is
composed of 7 days of data collection, followed by 10 days of iterative
optimisation, amounting to a total of 17 days. Each consecutive parameter
subset optimisation will add this amount of additional lead-time before the
configured parameters are applied to each cell.

In this work, we developed a technology to reduce the overall lead-time
required to optimise cell-level parameter configurations for a new wireless
network site. We formulated the problem of stage-wise optimisation of
parameter subsets as a _Continual Reinforcement Learning_ problem. This
approach leverages forward transfer of knowledge between optimisation policies
with overlapping subsets of actions in order to learn the ultimate policy in a
data-efficient task-oriented fashion. Additionally it allows for a safe
rollback to a policy of a previous subset, if objective KPIs do not improve,
by avoiding catastrophic forgetting. Through a series of experiments, we
demonstrate a two-fold reduction in deployment lead-time compared to a
_Reinitialise-and-Retrain_ baseline. The main challenges addressed in this
paper, with corresponding technical solutions, are summarised in Table 4.2.

Table 1: Challenges and technical solutions.

o 1 | X[l] | X[l] | Challenge | Technical Solution
---|---
Data limitation in real wireless network trials, which are expensive and constrained to real-time. | Model-based RL.
Time constraints on policy deployment lead-time to address a new candidate configuration parameter subset. | Continual RL.
High levels of noise in objective KPIs. | Probabilistic reward model explicitly accounting for aleatoric and epistemic uncertainty.
Inference time constrained to under 5 minutes for 20,000 cells. | Learning-based solution as opposed to online planning at decision time.

The rest of the paper is organised as follows. Section 2 summarises previous
work on Reinforcement Learning (RL) for network parameter optimisation, and
outlines the main classes of methods for Continual RL (CRL). Section 3
presents the detailed description of the real-world dataset, the problem
formulation, and the solution methods. Section 4 describes the experiment
design, and Section 5 analyses the experiment results. Finally, we conclude in
Section 6 and propose future research.

##  2 Related Work

Previous research has addressed the problem of wireless network parameter
optimisation using rule-based methods [10, 11], mathematical models [19], or
RL [1, 3, 4, 5, 8, 17, 27, 30]. A common characteristic of the aforementioned
works is that the action space (adjustable network configuration parameters)
is defined a priori at design time and is kept fixed. Scenarios where one
wishes to dynamically extend the action space with additional configuration
parameters can be solved in a sample-efficient way via a Continual RL problem
formulation [13]. The main classes of methods for achieving positive forward
transfer while mitigating catastrophic forgetting in CRL tasks are summarised
as follows:

  * •

Parameter storage based methods require multiple independently trained models
to be stored for the different tasks. Catastrophic forgetting is overcome at
the cost of storing parameters for each model. The space needed to store
policies is linear in the number of tasks addressed. Unfortunately, this
technique does not support knowledge transfer across tasks regardless of their
potential similarities.

  * •

Distillation is the process of distilling a source model(s) to a target model
[22, 26, 29, 31]. Distillation mitigates the need to store multiple models by
compressing them into a single neural network. Hence, models trained on
several source tasks can be distilled into a single network which captures
shared experiences from all tasks. Despite reducing storage space,
distillation still requires task specific layers to extrapolate features or to
fit a task’s action policy.

  * •

Rehearsal consists of training using examples from both the current task and
old examples from previously encountered tasks. This requires either an
experience replay buffer to store old examples, or if storage space is
limited, pseudo-rehearsal whereby examples are synthesized using a generative
model [12, 15, 21]. When solving a novel task, models can retain performance
achieved in previous tasks by continuously revisiting examples from same.

  * •

Regularization based methods maintain a single model across multiple tasks.
Catastrophic forgetting is mitigated by preventing parameters that are
important for previous tasks from changing significantly when learning a new
task [14, 32]. When compared to Parameter Storage methods, this approach
reduces the storage space required since only one model is maintained. A
drawback is that the initial model must have sufficient capacity to
accommodate all future tasks, which are typically unknown upfront.
Furthermore, it is not clear how to handle extension of the action space.

  * •

Modular architectures facilitate CL by exploiting flexibility and
compositionability of neural networks. Model capacity can be adjusted
dynamically, making adaptation to unseen tasks possible without over-
parameterising the model initially [23, 26, 32]. It is sometimes helpful to
decompose complex problems into easier sub-problems. In this case, neural
modules can be combined and re-used across tasks [9, 16]. For instance, a
decision policy could be decomposed into one module to extract features and
another to select actions. Modules are combined as required to achieve
positive transfer, while catastrophic forgetting is prevented by storing
lightweight modules trained on prior tasks.

In our optimisation problem, the order and the type of overlap across the
sequence of cell-level configuration parameter subsets are defined by the
domain expert. We selected Progress-and-Compress (P&C) [26], a method that has
demonstrated good performance in CRL scenarios of evolving action space. P&C
is a hydrid between the methods of distillation, regularisation and modular
architectures, and it can be decomposed into two steps. Firstly, the progress
step, where the agent achieves positive transfer by expanding the architecture
horizontally via lateral connections from an existing knowledge base
(assimilates experience from all previous tasks) to a new active column
(learns skills on the new task). Secondly, in the compress step, model
distillation is used to distill the active column into the knowledge base
without disrupting performance on previous tasks. As such, P&C mitigates
unbounded growth in learned policies, while ensuring transfer across similar
tasks and preventing catastrophic forgetting.

##  3 Methods

In this section, the dataset is described and the problem is formulated as a
Markov Decision Process (MDP). Methods are introduced to achieve data
efficient continual learning of network control policies using model-based RL.

###  3.1 Description of the Dataset

A dataset $D_{init}$ was collected over 5 days in a real 5G network containing
966 cells, by executing random actions to random cells at hourly intervals. We
collected four different feature sets111The names of CPs, PCs, and EPs cannot
be disclosed.: cell-level configuration parameters (CPs) representing the
actions, performance counters (PCs), engineering parameters (EPs), and spatio-
temporal (ST) context. PCs include time-varying features (e.g. demand, channel
quality, etc.) which are reported by cells hourly. EPs are fixed
characteristics of a cell such as its tilt, antenna type, etc. The
optimisation objective KPI is cell-level throughput measured in Gigabits per
second (Gbps). Exploratory analysis revealed the following challenging
properties of the dataset.

High dimensional state and action spaces: The composite action space consists
of 19 CPs. Three different groups of CPs were adjusted during the data
collection experiment. The three groups include: 4 CPs for power control on
cells, 6 CPs for modulation code scheme selection (method used by a cell to
encode digital data), and 9 CPs for rank selection (number of spatial streams
used to transmit data from a cell to users). CP are categorical variables that
take between 2 and 13 different values – the number of unique action
combinations is $\approx 4.6\times 10^{18}$. Raw network state is the union of
PCs, EPs, and ST context. A total of 410 features constitute the network
state.

Throughput (TP) time-series exhibits high levels of noise: A seasonality-
trend-noise decomposition of the time-series was performed for each cell. It
was observed that $49.41\%\pm 10.04\%$ of the variance in TP is explained by
seasonal and trend components. The remaining variance can be attributed to the
noise component of the time-series.

The challenges and technical solutions are summarised in Table 4.2. A tabular
dataset was formed by aligning the four feature sets with respect to cell ID
and timestamp. The aligned dataset contained
$966\,\,[\mathrm{cells}]*5\,\,[\mathrm{days}]*24\,\,[\mathrm{hours}]$
examples. Rows with missing values due to energy saving enabled on some cells
in early morning hours were removed, resulting in a total of 103,648 examples.
CPs and EPs were transformed using min-max scaling. PCs were transformed using
the Yeo-Johnson method due to extreme values. Lastly, ST features were
cyclically encoded.

###  3.2 Problem Formulation

The problem is formulated as a MDP defined by the tuple
$\langle\mathcal{S},\mathcal{A},T,R,\gamma,D\rangle$, where $\mathcal{S}$ is
the state space, $\mathcal{A}$ is the action space,
${T}:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\rightarrow\mathbb{R}^{|\mathcal{S}|}$
is the transition function from a state $s$ to $s^{\prime}$ after taking
action $a$, $R:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is the
reward function, $\gamma\in(0,1]$ is the discount factor, and $D$ is the
initial state distribution. The goal is to find a policy
$\pi_{\theta}:S\xrightarrow{}\mathcal{A}$ that maximises expected cumulative
reward, defined as
$\mathbb{E}\left[\sum_{t=0}^{H}\gamma^{t}r(s_{t},a_{t})|s_{0}\sim
D,a_{t}\sim\pi_{\theta}(s_{t})\right]$, where $H$ is the planning horizon.

States: A raw state $s^{raw}_{t}\in\mathcal{S}$ at hour $t$ is composed of PCs
observed between $t-1$ and $t$, EPs, and ST features. $s^{raw}_{t}$ is
compressed into $s_{t}$ using an auto-encoder $g_{\psi}:S\rightarrow Z$,
$S\in\mathbb{R}^{410}$, $Z\in\mathbb{R}^{50}$.

Actions: The action space $\mathcal{A}$ is combinatorial, where each CP $i$
may take between 2 and 13 different discrete values represented by
$\mathit{CP}_{i}\equiv\\{\mathit{CP}_{i,1},\mathit{CP}_{i,2},$
$\ldots,\mathit{CP}_{i,n_{\mathit{CP}_{i}}}\\},$ giving a total of
$|\mathcal{A}|=\prod_{i=1}^{N}n_{\mathit{CP}_{i}}\approx 4.6\times 10^{18}$
distinct actions for all 19 CPs. We train the policy to generate an action
$a_{t}=(a_{1,t},a_{2,t},\ldots,a_{N,t})\in\mathcal{A}$ at the beginning of
hour $t$ given network state $s_{t}$ as input. Note that
$a_{i,t}\in\mathit{CP}_{i}$ is a discrete value from the set
$\mathit{CP}_{i}$.

Reward: Defined as the value of throughput KPI observed on a cell between $t$
and $t+1$.

Planning Horizon: Single-step episodes are considered, in which an action is
executed at time $t$, the reward is given at time $t+1$ and the episode
terminates.

###  3.3 Model-based Reinforcement Learning

Since we are optimising over a single-step horizon, the learned dynamics
function is limited to a one-step reward prediction model. This is defined as
a probabilistic model
$f_{\phi}(r_{t+1}|s_{t},a_{t})=Pr(r_{t+1}|s_{t},a_{t};\phi)$ that outputs the
conditional distribution of the reward given the current state and action.
Learning a reward model is a task of fitting an approximation model to the
true reward function given the training dataset
$D_{train}=\\{(s_{i}=g_{\psi}(s^{raw}_{i}),a_{i}=CP_{i},r_{i}=TP_{i})\\}_{i=1}^{M}$
collected from the real network, where $g_{\psi}:S\rightarrow Z$,
$S\in\mathbb{R}^{410}$, $Z\in\mathbb{R}^{50}$ is a feature extractor in the
form of an auto-encoder neural network trained to reconstruct the raw state
using dataset $\\{(s_{i}^{raw})\\}_{i=1}^{M}$. The motivation of pre-
processing of raw state using an auto-encoder is to produce compact control
polices for TP optimisation. Algorithm 1 presents the pseudo-code for model-
based RL.

Figure 1: Architecture of the policy neural network $\pi_{\theta}(s|a)$, which
is composed of $N$ sub-policies $\pi_{\theta}^{i}(a_{i}|s)$ each represented
by a separate output head. Output heads share an embedding of the state
computed by two hidden layers.

The architecture of the policy network $\pi_{\theta}(a|s)$ is shown in Figure
1. In our setting, the action components corresponding to CPs are chosen
independently from each other. Given $N$ CPs, we decompose policy
$\pi_{\theta}$ into $N$ sub-policies $\pi_{\theta}^{i}(a_{i}|s)$ such that
$\pi_{\theta}(a|s)=\prod_{i=1}^{N}\pi_{\theta}(a_{i}=CP_{i}|s)$. Decomposition
of the policy network is achieved by branching $N$ output heads from a shared
embedding of the state. The embedding is computed by two layers, where the
first layer takes the compressed state $s_{t}$ as input. The $i^{th}$ output
head returns a probability over the $i^{th}$ action component corresponding to
$CP_{i}$.

Proximal policy optimisation (PPO) [25] is employed to train
$\pi_{\theta}(a|s)$. In step $k$, the policy update is given by:

$\centering\theta_{k+1}={\arg\max}_{\theta}\,\,\mathbb{E}_{s,a\sim\pi_{\theta_{k}}}[{\cal
L}_{PPO}(s,a,\theta_{k},\theta)].\@add@centering$

where the loss function is defined as:

$\centering\mathcal{L}_{PPO}(s,a,\theta_{k},\theta)=\min\left(\frac{\pi_{\theta}(a|s)}{\pi_{\theta_{k}}(a|s)}A^{\pi_{\theta_{k}}}(s,a),\,\,\,\,g(\epsilon,A^{\pi_{\theta_{k}}}(s,a))\right).\@add@centering$
(1)

The advantage $A^{\pi_{\theta_{k}}}(s,a)$ is clipped using function $g(\cdot)$
to prevent large changes in the policy parameters:
$g(\epsilon,A)=(1+\epsilon)A$, if $A\geq 0$; otherwise
$g(\epsilon,A)=(1-\epsilon)A$. Advantages are estimated by subtracting the
mean reward in a mini-batch $\mathcal{B}$ from the current reward:

$A^{\pi_{\theta_{k}}}(s,a)=r(s,a)-\frac{1}{|\mathcal{B}|}\sum_{i\in\mathcal{B}}r(s_{i},a_{i}).$
(2)

Algorithm 1 Model-based Reinforcement Learning

1:Collect dataset
$D_{init}=\\{(s_{i}^{raw},a_{i}=\mathit{CP}_{i},r_{i}=\mathit{TP}_{i})\\}_{i=1}^{M}$
from real network.

2:Train feature extractor $g_{\psi}$ on dataset
$\\{(s_{i}^{raw})\\}_{i=1}^{M}$

3:Process $D_{init}$ into dataset
$D_{train}=\\{(s_{i}=g_{\psi}(s_{i}^{raw}),a_{i}=\mathit{CP}_{i},r_{i}=\mathit{TP}_{i})\\}_{i=1}^{M}$

4:Train reward model $f_{\phi}$ on $D_{train}$

5:Initialise control policy $\pi_{\theta}$

6:for E epochs do

7: shuffle $D_{train}$ and partition it into mini-batches.

8: for number of mini-batches do

9: initialise training set $D_{RL}$

10: for each example index i of a mini-batch do

11: $\hat{a}_{i}$ $\sim$ $\pi_{\theta}(s_{i})$

12: $\hat{r}_{i}$ = $f_{\theta}(s_{i},\hat{a}_{i})$

13: compute advantage $\hat{A}_{i}$ via Equation 2

14: add $(s_{i},\hat{a}_{i},\hat{A}_{i})$ to $D_{RL}$

15: update $\pi_{\theta}$ based on $D_{RL}$ using the Adam optimiser, given
loss in Equation 1

###  3.4 Compressing High Dimensional Network State

A total of 383 performance counters (PCs) are combined with 9 engineering
parameters (EPs), and 18 spatio-temporal features (ST) to form the raw network
state $s_{t}^{raw}\in\mathbb{R}^{410}$. An under-complete auto-encoder neural
network $g_{\psi}:S\rightarrow Z$, $S\in\mathbb{R}^{410}$,
$Z\in\mathbb{Z}^{50}$ is trained to compress the raw network state.

Compressed states produced by the encoder
$g_{\psi}^{\mathit{encoder}}(s_{t}^{\mathit{raw}})\rightarrow s_{t}$ are
passed to the decoder, which reconstructs the input:
$g_{\psi}^{\mathit{decoder}}(s_{t})\rightarrow\hat{s}_{t}^{\mathit{raw}}$. A
regulariser $g_{\psi}^{\mathit{reg}}$ predicts the throughput at $t+1$ given
$s_{t}$. Incorporating the regulariser was found to improve predictive
accuracy of the reward model and policy gain. Loss is defined as the mean
square error between the input state and reconstructed state:

$\mathcal{L}(s^{\mathit{raw}},r,\psi)=(\lambda_{\mathit{AE}})\left(\frac{1}{|S|}\sum_{i\in
S}\left(s^{\mathit{raw}}[i]-\hat{s}^{\mathit{raw}}[i]\right)^{2}\right)+(1-\lambda_{\mathit{AE}})(r-\hat{r})^{2},$

where $\hat{s}^{\mathit{raw}}$ is the reconstructed state produced by the
decoder at time $t$, $\hat{r}$ is predicted reward (throughput) at $t+1$ given
by the regulariser, and $\lambda_{\mathit{AE}}\in\left[0,1\right]$ weights
contribution of the decoder and regulariser. Subscripts $t$ are dropped for
clarity.

(a) Network State Compressor.

(b) Reward Model.

Figure 2: High dimensional raw network state $s^{raw}_{t}\in\mathbb{R}^{410}$
at hour $t$ is compressed using a modular auto-encoder to
$s_{t}\in\mathbb{R}^{50}$. The reward model is an ensemble of probabilistic
neural networks. A single component model $m$ is displayed.

###  3.5 Reward Model

The reward model $f_{\phi}:S\times A\rightarrow\mathbb{R}$ is based on the
dynamics model described by Chua et al. [6] termed as Probabilistic Ensemble
with Trajectory Sampling. Predicted reward is given by a probabilistic
ensemble consisting of 10 component models, each trained with different seeds
and bootstrap samples of $D_{train}$ collected from the real network. Figure
1(b) displays the architecture of a single component model. Each component
model $f_{\phi}^{i}$ parameterises a Gaussian distribution with diagonal
covariance:
$P_{\phi}(r_{t+1}|s_{t},\mathit{ST}_{t},a_{t})=\mathcal{N}(\mu_{\phi}(s_{t},\mathit{ST}_{t},a_{t}),\sigma_{\phi}(s_{t},\mathit{ST}_{t},a_{t})^{2})$,
from which the predicted reward (throughput) at $t+1$ can be sampled. Inputs
at hour $t$ include: the compressed network state $s_{t}$, current time and
location of the cell $\mathit{ST}_{t}$, and actions produced by the policy
$a_{t}=\mathit{CP}_{t}$.

We employ a causal inference-based regulariser to improve the generalisation
of the reward model similar to the work of [24]. Generalised Propensity Score
(GPS) [28] predicts the action $a_{t}=CP_{t}$ from a representation $\xi$ of
$a_{t}$, compressed state $s_{t}$, and $ST_{t}$. The composite loss function
contains terms for the negative log likelihood and GPS regulariser:

$\mathcal{L}(s,a,\phi)=\lambda_{\mathit{RM}}\left(\frac{\left({\mathit{r}}-\mu_{\phi}(s,a)\right)^{2}}{\sigma_{\phi}^{2}(s,a)}+\log{\sigma_{\phi}^{2}(s,a)}\right)+\left(1-\lambda_{\mathit{RM}}\right)\mathrm{KL}\left(\mathit{a}\|\widehat{\mathit{a}}\right),$

where $r$ is the ground truth reward (throughput), $\mu_{\phi}(s,a)$ and
$\sigma_{\phi}^{2}(s,a)$ are the predicted mean and variance, and
$\lambda_{\mathit{RM}}\in\left[0,1\right]$ balances the contribution of the
negative log likelihood and GPS regularisation terms. Subscripts $t$ are
dropped for clarify. The overall reward model architecture is shown in Figure
1(b).

Importantly, the probabilistic ensemble, consisting of component models
$m\in\mathcal{M}$, quantifies the uncertainty associated with the predicted
reward. Total uncertainty is decomposed into an epistemic term reflecting
model uncertainty due to limited data, and an aleatoric term which represents
irreducible noise:

$\sigma^{2}_{\mathit{epistemic}}(s,a)=\frac{1}{|{\mathcal{M}}|}\sum_{m\in\mathcal{M}}\left(\sigma^{2}_{\phi,m}(s,a)+\mu^{2}_{\phi,m}(s,a)\right)-\left(\frac{1}{|{\mathcal{M}}|}\sum_{m\in\mathcal{M}}\mu_{\phi,m}(s,a)\right)^{2},$
$\sigma^{2}_{\mathit{aleatoric}}(x)=\frac{1}{|{\mathcal{M}}|}\sum_{m\in\mathcal{M}}\sigma^{2}_{\phi,m}(s,a).$

Figure 3: Analysis of the reward model.

####  3.5.1 Reward Model Analysis:

Experiments were carried out to confirm that the reward model is an accurate
simulator of the network environment dynamics. The neural network ensemble was
bench-marked against a random forest containing 500 trees. Results from 10
fold cross-validation are displayed in Figure 3. Two conclusions are drawn
from the plots. Firstly, the lowest error is achieved by the reward model in
all folds. Secondly, error residuals are symmetric around zero, confirming
that throughput is neither systematically overestimated nor underestimated.

###  3.6 Continual Reinforcement Learning

Continual RL is formulated as a sequence of MDPs $MDP(1),\ldots,MDP(L)$ where
$\mathrm{MDP}(l)=\langle\mathcal{S}(l),\mathcal{A}(l),T(l),$
$R(l),\gamma,D(l)\rangle$ [13], $l\in\\{1,2,\ldots,L\\}$. In our case, each
action space $\mathcal{A}(l)$ is a subset of the next:
$\mathcal{A}(1)\subset\mathcal{A}(2)\subset\ldots\subset\mathcal{A}(L)$.
Therefore each consecutive MDP in the sequence represents a new task with an
expanding action space.

We employ the Progress-&-Compress (P&C) framework [26] for continual RL. P&C
enables the reuse of past information through layer-wise adaptors to a
knowledge base. The policy consists of two columns: the active column acquires
experience on a new task, and the knowledge base accumulates experience
acquired over all previous tasks. The active column is trained during the
progress phase. The knowledge base is then updated during the compress phase,
in order to subsume new knowledge acquired by the active column. Catastrophic
forgetting of the knowledge base is prevented through a regularisation
technique described below.

Figure 4: Training policies using progress and compress framework. Two hidden
layers, called as shared and branch, are dense layers with $128$ and $64$
neurons, respectively, with $\mathrm{tanh}$ non-linearity.

Progress Phase: When a new task is presented, the parameters of the knowledge
base are frozen – only those parameters in the active column are optimised.
Layer-wise connections between the knowledge base and the active column allow
reuse of representations produced by the knowledge base, thus enabling
positive transfer from previously learned tasks. Lateral adaptors are
implemented as multi-layer perceptrons. The $j$-th layer of the active column
is defined as follows:

$h_{j}=\sigma\left(W_{j}h_{j-1}+\alpha_{j}\odot
U_{j}\sigma\left(V_{j}h_{j-1}^{KB}+b_{j}^{KB}\right)+b_{j}\right),$

where $\alpha_{j}$ is a trainable vector initialised by sampling from ${\cal
U}(0,0.1)$ with size equal to number of units in layer $j$, and
$W_{j},U_{j},V_{j}$ are weight matrices. The active column is trained using
Algorithm 1.

Compress Phase: The compress phase distills the learned active column into the
knowledge base. At this stage, the parameters of the active column are frozen.
The compression optimisation objective is a distillation loss with the Elastic
Weight Consolidation (EWC) penalty. EWC protects the knowledge base against
catastrophic forgetting, such that all previously learned skills are
maintained:

${\cal
L}_{KB}\left(\theta_{T_{t}}^{KB}\right)=\mathbb{E}\left[\mathrm{KL}(\pi_{\theta_{T_{t}}}(\cdot|s)\|\pi_{\theta_{T_{t}}^{KB}}(\cdot|s))\right]+\mathrm{EWC}\left(\theta_{T_{t}}^{KB},\theta_{T_{s}}^{KB}\right)$

where $\pi_{\theta_{T_{t}}}$ is the policy of active column after learning
task $T_{t}$ while $\pi_{\theta_{T_{t}}^{KB}}$ is the knowledge base into
which the active column is distilled.

The authors in [26] proposed online EWC to reduce the computational cost of
calculating the Fisher information matrix:

$\mathrm{EWC}(\theta_{T_{t}}^{KB},\theta_{T_{s}}^{KB})=\frac{1}{2}\sum_{i\in{\cal
I}^{KB}}\gamma_{\mathit{Fisher}}\left(F^{*}_{T_{s},i}(\theta^{KB}_{T_{t},i}-\theta^{KB}_{T_{s},i})^{2}\right)$

where
$F^{*}_{T_{t}}=\gamma_{\mathit{Fisher}}\left(F^{*}_{T_{s}}\right)+\frac{1}{|\mathcal{B}|}\sum_{(s,a)\in\mathcal{S}\times\mathcal{A}}\nabla\log\pi_{\theta_{T_{t}}^{KB}}(a|s)\left(\nabla\log\pi_{\theta_{T_{t}}^{KB}}(a|s)^{\intercal}\right)$,
${\cal I}^{KB}$ are layers in the knowledge base excluding the output layers,
$\mathcal{B}$ is the set of examples in the mini-batch, and
$\gamma_{\mathit{Fisher}}\in\mathbb{R}$ is a hyper-parameter.

##  4 Experimental Evaluation

A number of experiments are designed to assess the sample-efficiency of
continual RL in the wireless optimisation scenarios, where the action space
$CP_{source}$ in a source task $T_{s}$ is expanded to a set $CP_{target}$ in
the target task $T_{t}$, where $CP_{source}\subset CP_{target}$. We are
contrasting P&C against a baseline _Reinitialise-and-Retrain_ (R&R) method,
which initialises $\pi_{\theta}(s|a)$ and trains it solely on the dataset
$D_{init}$ specific to the CP set of a task in Algorithm 1.

###  4.1 Wireless Network Optimisation Scenarios

Figure 5: Policies trained via Reinitialise-&-Retrain (R&R) use four days of
data from the target task $(T_{t})$. Policies trained via Progress-&-Compress
(P&C) use two days of data from $T_{t}$ (50% less), and transfer knowledge
from a source task $T_{s}$.

We defined three scenarios that demonstrate how domain knowledge of a network
operator is used to optimise sets of CPs for a new network site in a stage-
wise manner.

Scenario 1: Use of prior knowledge about the environment in which the cell is
deployed (e.g. height of surrounding buildings, density of buildings) and the
profile of services that are expected in this geographical area (e.g. ultra
reliable low-latency communication). In this scenario, based on domain
knowledge, the operator chooses initially a small subset of CPs known to have
high impact to the optimisation objective KPIs. The small subset is then
extended to the full set of CPs to leverage additional CP interactions.

Scenario 2: Use of prior knowledge about the potential of certain CPs to cause
network disruptions or network performance degradation (i.e. exploring CPs
related to energy saving sleep-modes of a base-station, or CPs that determine
antenna configurations directly impacting coverage). In this scenario, network
optimisation starts with a set of CPs known to have low impact on the
optimisation objective KPIs in order to reduce said risks. The initial set is
then gradually extended to the full set of CPs.

Scenario 3: A combination of considerations reflected in Scenarios 1 and 2.

###  4.2 Policy Training Setup

Five experiments were designed to assess policy deployment lead-time reduction
in the three scenarios described above. The experimental setups are outlined
in Table 4.2. In each experimental setup, policies trained using P&C were
initialised on a source task $T_{s}$ and updated on the target task $T_{t}$.
Policies trained using R&R were trained on the $T_{t}$ using twice the
training data volume as P&C (4 days as opposed to 2 days).

Table 2: Experimental setups to assess policy deployment lead-time reduction
after expanding the set of CPs from a source task ($T_{s}$) to target task
($T_{t}$).

o 1 | X[l] | X[l] | X[l] | X[l] | Experiment | Scenario | CPs in $T_{s}$ | CPs in $T_{t}$
---|---|---|---
1 | 1 | 4 high-impact | All 19
2 | 2 | 15 low-impact | All 19
3 | 1 | 2 high-impact | 4 high-impact
4 | 3 | 2 high-impact & 7 low-impact | All 19
5 | 2 | 7 low-impact | 15 low-impact

To ensure an unbiased evaluation, the dataset $D_{init}$ was partitioned into
disjoint subsets as follows:

  * •

30% of the examples were reserved to train a test system model with similar
structure to the reward model.

  * •

Of the 70% remaining; 70% of the examples were used to train the reward model,
state compressor, and policies, and 30% was used to assess TP gain as
described in the following section.

The datasets were further structured as outlined in Figure 5 in order to
evaluate P&C under a 50% reduction in data volume compared to R&R.

Table 3: Hyperparameters used in the experimental evaluation.

o 1 | X[l] | X[l] | X[l] | X[l] | X[l] | Hyperparameter | State Compressor | Reward Model | R&R | P&C
---|---|---|---|---
Training Iterations | $\mathit{epochs}=100$ | $\mathit{epochs}=100$ | $\mathit{epochs}=100$ | $\mathit{steps}=1000$
Learning Rate | 0.00015 | 0.00015 | 0.0001 | 0.0001
Batch Size | 128 | 128 | 64 | 64
$\gamma_{\mathit{Fisher}}$ | NA | NA | NA | 0.9
$\lambda$ | $\lambda_{\mathit{AE}}=0.9$ | $\lambda_{\mathit{RM}}=0.9$ | NA | NA

###  4.3 Performance Metrics

Policy deployment is constrained by real-time, and its lead-time is
approximately equivalent to the time required for the collection of dataset
$D_{init}$ in Algorithm 1. Our target is to demonstrate (i) 50% reduction in
the amount of data required by P&C relative to R&R, and (ii) positive TP gain.
TP gain in hour $t$ is defined as:

$\mathit{TP\,Gain}(t)\,\,[\%]=\frac{\sum_{c\in\mathcal{C}}\widehat{TP}_{c,t}^{P\&C}-\sum_{c\in\mathcal{C}}\widehat{TP}_{c,t}^{R\&R}}{\sum_{c\in\mathcal{C}}\widehat{TP}_{c,t}^{R\&R}}\times
100\,\%,$ (3)

where $\mathcal{C}$ is the set of all cells, $\widehat{TP}_{c,t}^{P\&C}$ is
the predicted TP when cell $c\in\mathcal{C}$ executes CPs given by the P&C
policy and $\widehat{TP}_{c,t}^{R\&R}$ is the predicted TP when cell $c$
executes CPs given by the R&R policy. Predictions are given by the test system
model. TP gain is computed in the test set as described in Section 4.2.

##  5 Results

Results are reported based on 20 independent runs of each experiment in Table
4.2.

####  5.0.1 Training performance:

Figure 6 shows the learning curves of P&C and R&R. First, we observe that the
runs of experiments for Scenario 1 that start with high-impact CPs in the
source task converge faster than those runs of Scenario 2, which start with
low-impact CPs is source task. A context change occurs between September 12th
and 13th when training R&R policies incrementally. The context change is
associated with switching from the weekend to a week day.

Figure 6: Learning curves. Blue shading shows standard deviation over 20
independent runs.

The jump in mean reward per epoch on September 13th is caused by a change in
the underlying user demand. This change is due to more users occupying
hotspots such as downtown business districts, transport hubs, etc. during
weekdays. These areas are served by high capacity cells, resulting in a higher
overall throughput.

Figure 7: Assessment of TP gain.  Figure 8: Training resource footprint.

####  5.0.2 Assessment of TP gain:

Figure 7 illustrates the distribution of median TP gain over 20 independent
runs calculated at hourly intervals as defined in Equation 3. The median of TP
over 20 independent runs for the P&C and R&R methods is given by
$TP_{c,t}^{X}=\mathrm{median}(TP_{c,t,1}^{X},TP_{c,t,2}^{X},\ldots,TP_{c,t,20}^{X}),\quad
X\in\\{P\&C,R\&R\\}$, where $TP_{c,t,run}^{X}$ is the throughput for cell $c$
at time $t$ for $run=1,2,\ldots,20$.

Results displayed in Figure 7 confirm that P&C achieves a positive median TP
gain along with a data reduction of 50% relative to the baseline R&R policy in
experiments 1, 2, 4, and 5. The highest median TP gain is attained in
experiment 4 due to the presence of both high-impact and low-impact CPs in the
source task. P&C is most effective in experiments 1, 2, and 4 with all 19 CPs
in the target task. We observe smaller positive median TP gain in experiment 5
because of the low impact CPs in both source and target tasks, while a
negative median TP gain is observed in experiment 3 due to few low-impact CPs
in the source task. This result suggests that the effectiveness of continual
RL as a sequence of MDPs with an expanding combinatorial action space is
mostly affected by the strength of the causal effect of each consecutive
action subset.

####  5.0.3 Training resource cost:

The amount of training RAM and convergence wall time are illustrated in Figure
8. $P\&C$ results in a reduction of approx. 31% in memory footprint compared
to $R\&R$ baseline, amounting to a total of $6.7$ GB. This is crucial for
training in resource-constrained hardware, which is typical for on-premise
deployments in telecommunication networks. Average training wall time for
$P\&C$ is measured at 80 minutes compared to 140 minutes for $R\&R$. Reduced
training time allows us to allocate more resources for hyper-parameter tuning
and model selection, and enables more frequent model updates in highly-dynamic
wireless networks.

####  5.0.4 Inference time:

Inference for setting configuration parameters reduces to a forward pass of
the compact neural net policy, and it was measured at 80 milliseconds of wall-
time for 20,000 cells.

##  6 Conclusions and Future Work

Learning optimisation policies from scratch in order to assess the performance
of different subsets of configuration parameters on the objective KPI is
costly is terms of time required for data collection, particularly for unknown
wireless network environments with exploration challenges (i.e. safety
constraints, limited budget of interactions with the wireless network, high
dimensional state/action spaces, high noise levels). The contribution of our
work is the application of Continual Reinforcement Learning to a wireless
network optimisation system where the number of decision variables gradually
increase throughout its lifetime. By refraining from reinitialising-and-
retraining the policy, we achieve a two-fold reduction in end-to-end
deployment time for each consecutive decision variable subset, which has
direct implications in reducing operational cost and overall data collection
requirement.

We defined three scenarios reflecting how domain knowledge of a network
operator can be leveraged to handcraft a series of CP subsets to be optimised
stage-wise. Five experiments based on these scenarios were conducted, and
empirical results demonstrated a positive throughput gain of up to 4% of P&C
continual RL method against the R&R baseline that used double the amount of
training data.

The selection of CP subsets for optimisation policies deployed to new wireless
network sites is traditionally an iterative trial-and-error process that
predominantly relies on domain knowledge. Going forward, there are two
approaches that will be investigated to further improve sample efficiency. The
first approach is based on causal structure learning to automatically identify
the CPs that exert the strongest causal effect on the optimisation objective
KPI in a new environment, and hence complement and fine-tune operator’s domain
knowledge during selection of CP subsets. Of applicability here are methods
that admit latent confounders when learning causal graphs, for example the
works of [7, 18, 20] and also the auto-causal tuning method of [2]. The second
approach is based on policy warm-starting through the technology of _sim-to-
real transfer_ in RL [33], some examples of which include domain adaptation,
imitation learning and policy distillation.

## References

  * [1] Balevi, E., Andrews, J.G.: Online antenna tuning in heterogeneous cellular networks with deep reinforcement learning. IEEE Trans. Cogn. Commun. Netw. 5(4), 1113–1124 (2019). https://doi.org/10.1109/TCCN.2019.2933420, https://doi.org/10.1109/TCCN.2019.2933420
  * [2] Biza, K., Tsamardinos, I., Triantafillou, S.: Out-of-sample tuning for causal discovery. IEEE Transactions on Neural Networks and Learning Systems pp. 1–11 (2022). https://doi.org/10.1109/TNNLS.2022.3185842
  * [3] Bothe, S., Masood, U., Farooq, H., Imran, A.: Neuromorphic AI empowered root cause analysis of faults in emerging networks. In: IEEE International Black Sea Conference on Communications and Networking, BlackSeaCom 2020, Odessa, Ukraine, May 26-29, 2020. pp. 1–6. IEEE (2020). https://doi.org/10.1109/BlackSeaCom48709.2020.9235002, https://doi.org/10.1109/BlackSeaCom48709.2020.9235002
  * [4] Bouton, M., Farooq, H., Forgeat, J., Bothe, S., Shirazipour, M., Karlsson, P.: Coordinated reinforcement learning for optimizing mobile networks (09 2021)
  * [5] Calabrese, F.D., Wang, L., Ghadimi, E., Peters, G., Hanzo, L., Soldati, P.: Learning radio resource management in rans: Framework, opportunities, and challenges. IEEE Communications Magazine 56, 138–145 (2018)
  * [6] Chua, K., Calandra, R., McAllister, R., Levine, S.: Deep reinforcement learning in a handful of trials using probabilistic dynamics models. Advances in neural information processing systems 31 (2018)
  * [7] Colombo, D., Maathuis, M.H., Kalisch, M., Richardson, T.S.: Learning high-dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics 40(1), 294–321 (2012)
  * [8] Dandanov, N., Al-Shatri, H., Klein, A., Poulkov, V.: Dynamic self-optimization of the antenna tilt for best trade-off between coverage and capacity in mobile networks. Wirel. Pers. Commun. 92(1), 251–278 (jan 2017). https://doi.org/10.1007/s11277-016-3849-9, https://doi.org/10.1007/s11277-016-3849-9
  * [9] Devin, C., Gupta, A., Darrell, T., Abbeel, P., Levine, S.: Learning modular neural network policies for multi-task and multi-robot transfer. In: 2017 IEEE international conference on robotics and automation (ICRA). pp. 2169–2176. IEEE (2017)
  * [10] Eckhardt, H., Klein, S., Gruber, M.: Vertical antenna tilt optimization for lte base stations. In: VTC Spring. pp. 1–5. IEEE (2011), http://dblp.uni-trier.de/db/conf/vtc/vtc2011s.html#EckhardtKG11
  * [11] Eisenblätter, A., Geerdes, H.: Capacity optimization for UMTS: bounds and benchmarks for interference reduction. In: Proceedings of the IEEE 19th International Symposium on Personal, Indoor and Mobile Radio Communications, PIMRC 2008, 15-18 September 2008, Cannes, French Riviera, France. pp. 1–6. IEEE (2008). https://doi.org/10.1109/PIMRC.2008.4699919, https://doi.org/10.1109/PIMRC.2008.4699919
  * [12] Isele, D., Cosgun, A.: Selective experience replay for lifelong learning. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 32 (2018)
  * [13] Khetarpal, K., Riemer, M., Rish, I., Precup, D.: Towards continual reinforcement learning: A review and perspectives. J. Artif. Intell. Res. 75, 1401–1476 (2020)
  * [14] Kirkpatrick, J., Pascanu, R., Rabinowitz, N., Veness, J., Desjardins, G., Rusu, A.A., Milan, K., Quan, J., Ramalho, T., Grabska-Barwinska, A., et al.: Overcoming catastrophic forgetting in neural networks. Proceedings of the national academy of sciences 114(13), 3521–3526 (2017)
  * [15] Li, Z., Hoiem, D.: Learning without forgetting. IEEE transactions on pattern analysis and machine intelligence 40(12), 2935–2947 (2017)
  * [16] Mendez, J.A., van Seijen, H., Eaton, E.: Modular lifelong reinforcement learning via neural composition. arXiv preprint arXiv:2207.00429 (2022)
  * [17] Nasir, Y.S., Guo, D.: Multi-agent deep reinforcement learning for dynamic power allocation in wireless networks. IEEE Journal on Selected Areas in Communications 37(10), 2239–2250 (2019). https://doi.org/10.1109/JSAC.2019.2933973
  * [18] Ogarrio, J.M., Spirtes, P., Ramsey, J.: A hybrid causal search algorithm for latent variable models. In: Proceedings of the Eighth International Conference on Probabilistic Graphical Models. vol. 52, pp. 368–379. PMLR (06–09 Sep 2016)
  * [19] Partov, B., Leith, D.J., Razavi, R.: Utility fair optimization of antenna tilt angles in lte networks. IEEE/ACM Trans. Netw. 23(1), 175–185 (feb 2015). https://doi.org/10.1109/TNET.2013.2294965, https://doi.org/10.1109/TNET.2013.2294965
  * [20] Raghu, V.K., Ramsey, J., Morris, A., Manatakis, D., Spirtes, P., Chrysanthis, P.K., Glymour, C., Benos, P.: Comparison of strategies for scalable causal discovery of latent variable models from mixed data. International Journal of Data Science and Analytics 6, 33 – 45 (2018)
  * [21] Rolnick, D., Ahuja, A., Schwarz, J., Lillicrap, T., Wayne, G.: Experience replay for continual learning. Advances in Neural Information Processing Systems 32 (2019)
  * [22] Rusu, A.A., Colmenarejo, S.G., Gulcehre, C., Desjardins, G., Kirkpatrick, J., Pascanu, R., Mnih, V., Kavukcuoglu, K., Hadsell, R.: Policy distillation. arXiv preprint arXiv:1511.06295 (2015)
  * [23] Rusu, A.A., Rabinowitz, N.C., Desjardins, G., Soyer, H., Kirkpatrick, J., Kavukcuoglu, K., Pascanu, R., Hadsell, R.: Progressive neural networks. arXiv preprint arXiv:1606.04671 (2016)
  * [24] Saini, S.K., Dhamnani, S., Ibrahim, A.A., Chavan, P.: Multiple Treatment Effect Estimation using Deep Generative Model with Task Embedding. In: The World Wide Web Conference. pp. 1601–1611 (2019)
  * [25] Schulman, J., Wolski, F., Dhariwal, P., Radford, A., Klimov, O.: Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347 (2017)
  * [26] Schwarz, J., Czarnecki, W., Luketina, J., Grabska-Barwinska, A., Teh, Y.W., Pascanu, R., Hadsell, R.: Progress & Compress: A scalable framework for continual learning. In: International Conference on Machine Learning. pp. 4528–4537. PMLR (2018)
  * [27] Shafin, R.S.B., Chen, H., Nam, Y., Hur, S., Park, J., Zhang, J., Reed, J.H., Liu, L.: Self-tuning sectorization: Deep reinforcement learning meets broadcast beam optimization. IEEE Trans. Wirel. Commun. 19(6), 4038–4053 (2020). https://doi.org/10.1109/TWC.2020.2979446, https://doi.org/10.1109/TWC.2020.2979446
  * [28] Shi, C., Blei, D., Veitch, V.: Adapting Neural Networks for the Estimation of Treatment Effects. Advances in neural information processing systems 32 (2019)
  * [29] Traoré, R., Caselles-Dupré, H., Lesort, T., Sun, T., Cai, G., Díaz-Rodríguez, N., Filliat, D.: Discorl: Continual reinforcement learning via policy distillation. arXiv preprint arXiv:1907.05855 (2019)
  * [30] Vannella, F., Jeong, J., Proutière, A.: Off-policy learning for remote electrical tilt optimization. In: 92nd IEEE Vehicular Technology Conference, VTC Fall 2020, Victoria, BC, Canada, November 18 - December 16, 2020\. pp. 1–5. IEEE (2020). https://doi.org/10.1109/VTC2020-Fall49728.2020.9348456, https://doi.org/10.1109/VTC2020-Fall49728.2020.9348456
  * [31] Yin, H., Pan, S.J.: Knowledge transfer for deep reinforcement learning with hierarchical experience replay. In: Thirty-First AAAI conference on artificial intelligence (2017)
  * [32] Yoon, J., Yang, E., Lee, J., Hwang, S.J.: Lifelong learning with dynamically expandable networks. arXiv preprint arXiv:1708.01547 (2017)
  * [33] Zhao, W., Queralta, J.P., Westerlund, T.: Sim-to-real transfer in deep reinforcement learning for robotics: a survey. In: 2020 IEEE Symposium Series on Computational Intelligence (SSCI). pp. 737–744 (2020). https://doi.org/10.1109/SSCI47803.2020.9308468