# Vacancy localization effects on MX2 transition metal dichalcogenides: a
systematic ab-initio study

Rafael L. H. Freire<EMAIL_ADDRESS>Felipe Crasto de Lima
<EMAIL_ADDRESS>Adalberto Fazzio<EMAIL_ADDRESS>Ilum
School of Science, CNPEM, 13083-970 Campinas, Brazil

###### Abstract

Two-dimensional transition metal dichalcogenides (MX2) vacancy formation
energetics is extensively investigated. Within an ab-initio approach we study
the MX2 systems, with M=Mo, W, Ni, Pd and Pt, and X=S, Se, and Te. Here we
classify that chalcogen vacancies are always energetic favorable over the
transition metal ones. However, for late transition metals Pd $4d$, and Pt
$5d$ the metal vacancy are experimentally achievable, bringing up localized
magnetic moments within the semiconducting matrix. By quantifying the
localization of the chalcogen vacancy states we evidentiate that it rules the
intra- and inter-vacancy interactions that establish both the number of
vacancy states neatly lying within the semiconducting gap, as well as its
electronic dispersion and SOC splitting. Combining different vacancies and
phase variability 1T and 1H of the explored systems allow us to construct a
guiding picture for the vacancy states localization.

##  I Introduction

Defects on materials, and particularly atomic vacancies, arises as a possible
route for an experimental control of properties [1; 2]. Although the plethora
of possible materials defects, ranging from localized to extended, even the
elementary structural defect (atomic single vacancy) allow for an complex
phenomena to arise, for instance ruling the molecular self-assemblies over
two-dimensional (2D) substrates [3]. Additionally, controlling such punctual
defects on the precise manipulation of vacancies have allowed the construction
and engineering of metamaterials [4; 5]. In order to further explore materials
modification through vacancies, a clear picture of the intrinsic vacancy
properties are needed.

Two-dimensional transition metal dichalcogenides (TMD’s) are an interesting
class of materials holding tremendous potential for applications, encompassing
catalysis [6; 7], (opto)electronics [8; 9; 10], spintronics [10; 11],
magnetism [12] and energy storage [13; 14]. Particularly, the presence of
defects in 2D-TMD materials can directly impact their properties, causing a
variety of phenomena that could be either detrimental acting like carrier
scattering and localization sites [15] or beneficial as is the case of active
catalytic sites, and magnetic orderings [16]. Furthermore, spin-orbit effects
on those materials combined with point defects may give rise to interesting
phenomena such as giant magnetoresistance [17], topological phases [18], all
of which could be employed on spintronic devices. Additionally, recent
experimental works have shown that the localization character of the vacancy
states induces electronic transport [19], and magnetism ordering on transition
metal vacancies [20].

In this study, we perform ab-initio calculations based on the density
functional theory to clarify and quantify the localization character of
vacancies states on TMD. We considered the pristine and defected MX2 2D
systems (for M $=$ Ni, Mo, Pd, W, and Pt, and X $=$ S, Se, and Te) in their
1H, 1T, and 1T’ phases (see Fig. 1). We highlight the most stable structural
phase and most stable vacancies (of metals or chalcogen). We present an
analysis of the localization nature of such vacancy states and their
consequences to the energetic, electronic, and magnetic properties.

Figure 1: 1H, 1T, and 1T’ MX2 2D structural phases.

##  II Computational Approach

Spin-polarized calculations based on density functional theory (DFT) [21; 22]
were performed within the semi-local exchange-correlation functional proposed
by Perdew–Burke–Ernzerhof (PBE) [23]. To treat the long range dispersion van
der Waals (vdW) interactions, the pairwise D3 correction framework as proposed
by Grimme were considered [24; 25].

Figure 2: (a) Relative cohesive energies for all pristine systems and
structural parameters, namely (b) in-plane lattice and (c) monolayer
thickness, for the most stable MX2 structures.

For the total energies, $E_{tot}^{\text{DFT+vdW}}$, the electron-ion core
interactions were considered within the the projector augmented wave (PAW)
method [26; 27], as implemented in the Vienna Ab-Initio Simulation Package
(VASP) [28; 29]. The spin-orbit coupling was considered for the most stable
configuration. For all calculations, cutoff energy for the plane-wave
expansion of the Kohn–Sham orbitals was set to $520$ eV, under an energy
convergence parameter of $10^{-6}$ eV, with all atoms relaxed until the atomic
forces on every atom were smaller than $10^{-3}$ eV Å-1. A uniform $4\times
4\times 1$ k-point mesh was considered for the Brillouin zone (BZ)
integration, while a thinner mesh grid, $8\times 8\times 1$, was used for
Density of States (DOS) calculations.

The 2D in-plane lattice parameters ($a$ and $b$) were optimized for a $1\times
1$ unit cell with a fixed vacuum distance ($c$) of at least $20$ Å. Those
lattice parameters were used to build up the pristine and single vacancy
$5\times 5$ and $3\times 3$ unit cells to the 1H and 1T (1T’) structures, with
the former being large enough to avoid spurious interaction between periodic
images.

##  III Results and discussion

###  III.1 Structural Stability

To analyze the structural stability of the pristine structures, we calculated
the cohesive energies – being the total energy difference between the material
and each isolated atom – for the 1H, 1T, and, 1T’ phases. To better present
the relative stability of one phase against each other, we evaluated their
relative cohesive energy, namely $\Delta E^{1H-1T}$, $\Delta
E^{1H-1T^{\prime}}$, and $\Delta E^{1T-1T^{\prime}}$, as shown in Fig. 2 (a).
That is, $\Delta E^{1H-1T}>0$ indicates the 1H phase is more stable than the
1T phase, and similarly to the other relative energies. As one can see, the
negative values for 1H-1T ($200$–$400$ meV atom-1) and 1H-1T´ ($0$–$200$ meV
atom-1) for MoX2 and WX2 indicate the 1H phase as the most stable, with this
stability decreasing with the increase of the chalcogen atomic radius. Despite
the slightly positive values observed for NiX2, PdX2, and PtX2, the cohesive
energy differences point to the 1T phase as the most stable, with the 1T-1T’
differences being close to zero. Although the calculation started in the 1T’
phase, after relaxation it converged barrierlessly into the 1T phase.

In general, the 1T phase has larger in-plane lattice parameters [Fig. 2(b)]
with a smaller monolayer thickness as compared to the 1H phase and its more
open structure (hexagonal holes). Consequently the monolayer thickness [Fig.
2(c)] decreases as the chalcogen planes move down to stabilize the equilibrium
covalent bond length. In the 1H phase, the in-plane lattice parameter and
monolayer thickness increase going from S$\rightarrow$Te, as a result of
accommodating larger chalcogen atoms within the structure. In the case of the
1T phase, only PtX2 follows a similar trend as the 1H systems. On the other
hand, 1T-NiX2 and 1T-PdX2 undergo larger relaxation effects.

###  III.2 Vacancy formation energy

After pristine systems optimization, we built up a $5\times 5$ defective
monolayers considering the most stable phases for each system, namely the 1H
phase for MoX2 and WX2 systems, and the 1T phase for NiX2, PdX2, and PtX2.
Both native point defects comprising chalcogen (VX) and transition metal (VM)
vacancies were considered with the vacancy-vacancy distances in the 1H (1T)
around $15.8$–$17.6$Å ($15.7$–$20.0$Å).

In Figure 3, we present the formation energies, $E_{f}$, for both M and X
vacancies for the selected systems, evaluated according to

$E_{f}=E_{def}^{MX_{2}}-(E_{pristine}^{MX_{2}}-E_{M,X}^{free-atom}),$ (1)

in which $E_{def}^{MX_{2}}$ is the energy of the defected system;
$E_{pristine}^{MX_{2}}$ is the energy of the pristine system; and
$E_{M,X}^{free-atom}$ is free-atom energy of the corresponding atom M or X
which generates the vacancy defect.

As indicated in Fig. 3 both formation energies are endothermic, with the
formation energies of transition metal vacancies [Fig. 3 (a)] are higher than
that for the chalcogen [Fig. 3 (b)]. Considering the different structure
phases, one can realize that the defect formation energies in the 1T phase are
always smaller than those for the 1H phase. Particularly, Ni and Pd systems
with the heavier chalcogen (Se and Te) present transition metal vacancy
formation energy close to the scale values of the chalcogen vacancies,
indicating a possible occurrence of such defects (although less favorable than
X-vacancies). Additionally, for the X-vacancies, we present the formation
energy in a $3\times 3$ supercell, i.e. increased vacancy density, where the
interaction between neighboring vacancies is greater. As we can see the NiSe2,
NiTe2, and PdTe2 structures present a larger difference in the chalcogen
vacancy formation energy in the $3\times 3$ cell when compared with the lower
vacancy density case. As we will show in the next sections, such behavior is a
consequence of the more delocalized nature of the vacancy states in those
systems.

Figure 3: Native point defect formation energies. (a) is the transition metal
vacancy defect, and (b) the chalcogen defect. In the latter, we show formation
energies for both $5\times 5$ (dashed lines) and $3\times 3$ (solid lines)
supercell sizes.

###  III.3 Band structure analysis

In Fig. 4, we show a defected $5\times 5$ supercell model for both 1H and 1T
structures alongside their characteristic band structures with and without
spin-orbit coupling. The introduction of the chalcogen vacancy leads to three
defect states, corresponding to the three transition metal dangling bonds over
the vacancy environment. The triangular $C_{3v}$ local environment split the
dangling bond states into $E$ and $A_{1}$ irreducible representations.
Additionally, the interaction of those three intra-vacancy dangling bonds in
the 1H structure is larger than the 1T, given the remaining chalcogen coupling
those TM. That is, this stronger coupling leads to a large energy separation
between the $E$ and $A_{1}$ irreducible representation. Such an interaction
picture leads to three vacancy states on the 1T phase neatly lying on the
system energy gap ($E$ and $A$ states), while for the 1H phase the two $E$
states remain in the energy gap, with the $A_{1}$ state being within the
valence bands. Additionally to this picture, the SOC effect has a role in
splitting the $E$ states into two sets. Such value of splitting is directly
related to the SOC contribution in each system. In Fig. 4(c) we show the SOC
splitting value at $\Gamma$ for the lower vacancy density system ($5\times 5$
cell). Here, the TM atom mostly rules the SOC contribution, with a lower
variance of the value with the chalcogen atom. Such variance is higher for WX2
systems, 50 meV, while lower for MoX2 systems, 25 meV. For the 5d TM, the 1H
tungsten phases present the highest splitting, with the WS2 reaching close to
0.2 eV, followed by the 1T Pt phases. For the 4d TM similar trend is observed
with the 1H Mo phases presenting higher SOC splitting than Pd 1T phases.
Interestingly, going towards a higher vacancy density, where the interaction
between adjacent vacancy states becomes important, the SOC trend is altered.
First note that the vacancy states become dispersive, Fig. 4(b) right panels.
The combined effect of the dispersion of the bands and SOC, induce a higher
gap opening at the $\Gamma$ point for the 1T phases compared with the 1H ones.
For instance, such a gap goes over 0.3 eV for PtTe2, with a dependence of the
chalcogen atom being more accentuated. The dependence arises given the
indirect coupling between adjacent vacancy states being ruled by the TMD
matrix [18]. To further explore such indirect coupling we quantify the
localization of such vacancy states.

Figure 4: (a) $5\times 5$ defect supercell for 1H and 1T phases as indicated;
(b), the universal schematic of band structure of each semiconducting phase,
respectively, with and without spin-orbit corrections. The green arrows
indicated spin-orbit splitting of the corresponding $k$-point; (c) Spin-orbit
splittings calculated for both $3\times 3$ and $5\times 5$ defect supercells.
$E^{\Gamma}_{g}$ is the spin-orbit splitting at the $\Gamma$ point.

###  III.4 (De)Localization

Figure 5: Inverse participation ratio (IPR) as a function of the Bloch
wavefunction ($|\psi_{n,k}\rangle$) energy. (a1)-(a3) The IPR for 1H-MX2
phases with M$=$Mo, and W and (b1)-(b3) The IPR for 1T-MX2 phases, with
M$=$Ni, Pd, and Pt from S to Te, respectively

To take into account (de)localization effects due to the vacancy formation, we
characterized the defective systems through the inverse participation ratio
(IPR) [30] given by

$\text{IPR}_{n,k}=\frac{\sum\limits_{i=1}^{N}|\langle
i|\psi_{n,k}\rangle|^{4}}{\left(\sum\limits_{i=1}^{N}|\langle
i|\psi_{n,k}\rangle|^{2}\right)^{2}},$ (2)

in which $\langle i|\psi_{n,k}\rangle$ was taken as the sum of the orbital
projected KS eigenstate for each site/atom $i$, such that $N$ is the total
number of atoms in the cell. Thus, for a fully localized state the IPR should
be one, while a delocalized one corresponds to the limit $\text{IPR}=1/N$. In
Fig. 5, we show the IPR values for each $|\psi_{n,k}\rangle$ state as a
function of its eigenvalue, for the most stable phases for chalcogen
vacancies. It is worth pointing out that those defective systems are in
general semiconducting, with exception of NiSe2, NiTe2, and PdTe2, which are
metallic, as indicated by the energy gap at Fig. 5 $x$-axis.

In the semiconducting cases, the localized states that appear within the gap
corresponding to the vacancy states as depicted in Fig. 4. The 1H phases
present a slightly larger IPR when compared to the 1T phases. Taking the 1H
phase as an example [Fig. 5(a1)-(a3)], after the chalcogen vacancy formation,
if (in the ideal case) only the three transition metal orbitals
($|M_{i}\rangle$) neighboring the vacancy contributes equally (with $\langle
M_{i}|\psi_{vac-state}\rangle=a$) to the localization, we have
$\text{IPR}=3a^{4}/(3a^{2})^{2}=1/3$. That is, the localization limit is $1/3$
at the chalcogen vacancy surroundings. As observed the IPR values are around
$0.20$ which is close to the ideal limit, being reduced due to the spread of
the vacancy states to hybridized neighboring orbitals. The same IPR limit
(IPR$=1/3$) is valid for the 1T structures [Fig. 5(b1)-(b3)]. However, despite
the vacancy states being within the bandgap, the values observed are slightly
small, lying around IPR$=0.15$. In this case, an enhanced environment
interaction is observed. To better visualize it, in Fig. 6 we show the vacancy
states squared wave function (partial charge density) of MoS2 and PtSe2 with S
and Se vacancy, which respectively present an IPR of $0.2$ and $0.1$. The 1T
phase forms a pyramidal-like configuration [see Fig. 6 (b)], in which the
localized state spreads to opposite surface chalcogen atoms neighboring the
dangling bond Pt atom. On the other hand, the 1H phase LDOS [see Fig. 6(a)]
are mostly localized in the M atoms close to the vacancy. Thus, the IPR for
these localized states decreases for the 1T phase as compared to the 1H one.
The spatial distribution of those vacancy states, although with different
spreads follow the same three-fold symmetry of the vacancy structure.

Figure 6: LDOS for the chalcogen-defected systems (a) 1H-MoS2 and (b)
1T-PtSe2, with the same isosurface of 0.007 e/Å3

.

In the light of the quantification of the localization through the IPR, the
chalcogen vacancy formation energy and the band structure dependence of the
vacancy density can be readily explained. For instance, given the more
localized nature of Mo, W, and Pt vacancy states all above IPR$=0.1$, dictates
that the adjacent vacancies will have lower interaction not changing the
formation energy significantly. However Ni and Pd, with their IPR$<0.1$
indicate a stronger interaction between adjacent vacancies reducing the
formation energy for higher vacancy densities. Particularly, for the 1T phases
with Te [Fig. 5(b3)] which presents the lower IPR values for the vacancy
states, presents also the greater variability on the chalcogen vacancy
formation energy, see Fig. 3(b). A similar analysis can be done following the
band structure. Here, the 1T phases allowing a longer range of interaction
between adjacent vacancies leads to higher dispersive states when compared
with the 1H structure [Fig. 4(b)]. Here the localization of the chalcogen
vacancy states leads to the interpretation of their interaction, additionally,
TM vacancies can also lead to localized magnetic effects.

###  III.5 Magnetism

Although less energetically favorable, TM vacancies can also be found on TMD,
where the introduction of such vacancies can induce a local net magnetic
moment. In Table 1 we summarize the systems presenting a net magnetic moment
after a TM vacancy is formed, while the ones not present on the table did not
present any magnetic properties. After the TM vacancy is formed, localized
magnetic moments arise on the neighboring chalcogen atoms. Here a
ferromagnetic (FM) phase can be stabilized for some explored systems, being
the only observed phase for the 1H structure. Interestingly, some of the 1T
systems can present an antiferromagnetic (AFM) arrange of such chalcogen
magnetic moments, which can be more stable than the ferromagnetic, as already
observed for 1T-PtSe2 [12]. Aware of this behavior, we probed this
antiferromagnetic configuration for our 1T-based based systems. The systems
which presented a possible antiferromagnetic phase were 1T-NiS2, 1T-PdS2,
1T-PdSe2, 1T-PtS2, and 1T-PtSe2. However, as shown in Table 1, for Ni- and Pd-
based systems the antiferromagnetic phase ($\Delta E^{AFM-FM}<0$) were not the
most stable. On the other hand, for Pt-based systems we found more stable
antiferromagnetic phases for 1T-PtS2 and 1T-PtSe2 systems with an energy
difference about $-54$ meV/vacancy, and $-33$ meV/vacancy, respectively. That
is, such values for the AFM phase together with the FM of 1T-PdSe2 [with
$\Delta E^{AFM-FM}=43$ meV/vacancy], dictates such magnetic configurations to
be robust close to the ambient temperature.

Table 1: Magnetic moments, $m$ ($\mu_{B}$), induced in the most stable phases with the introduction of vacancies. The last column, $\Delta E^{AFM-FM}$ (meV/vacancy) , is the energy difference between the antiferromagnetic (AFM) and ferromagnetic (FM) phases, with $\Delta E^{AFM-FM}<0$ indicating the AFM phase being more stable. MX2 | vacancy | $m$ | $\Delta E^{AFM-FM}$
---|---|---|---
1H-MoSe2 | Mo | 4.00 | –
1H-MoTe2 | Mo | 2.00 | –
1H-WTe2 | W | 1.96 | –
1T-NiS2 | Ni | 4.00 | 19
1T-PdS2 | Pd | 4.00 | 17
1T-PdSe2 | Pd | 4.00 | 43
1T-PtS2 | Pt | 4.00 | -54
1T-PtSe2 | Pt | 4.00 | -33

##  IV Conclusions

We systematically investigated the energetic and electronic properties of a
series of two-dimensional transition metal dichalcogenides (MX2, with M$=$Ni,
Mo, Pd, W, and Pt; and X$=$S, Se, and Te) presenting native point defects,
namely chalcogen and transition metal vacancies, in different structural
phases. Here, we found the chalcogen vacancy as the most stable for all
systems, with lower formation energy in the 1T phase (Ni, Pd, and Pt systems)
than those in the 1H phase (Mo and W systems). However, transition metal
vacancies can still be found under experimental conditions. In this sense, our
results show the appearance of localized magnetic moments induced by metal
vacancies for 1T-PtS2 and 1T-PtSe2 that could be applied to design new 2D
magnets. Furthermore, we have explored the localization effects of the
chalcogen vacancy states. Such, localized states give rise to three energy
levels that can be neatly lying on the TMD matrix energy gap. Here the
localization strength, quantified by the defect states inverse participation
rate, is shown to be greater in the 1H phases. This leads to both (i) a
stronger repulsion between the three defect states (with $C_{3v}$, $E$ and
$A_{1}$ irreducible representations) increasing the gap between $E$ and
$A_{1}$ states, and (ii) giving rise to lower dispersion for higher vacancy
densities, that is, present a lower vacancy-vacancy interaction. For the 1T
phases, the more delocalized nature of the vacancy states gives rise to a
stronger hopping-like interaction between adjacent vacancies. Additionally, we
have shown that vacancy-vacancy interactions (given by their localization)
ruled not only the materials band dispersion but also the SOC splittings. This
investigation brings insightful discussions on the nature of energetic and
electronic effects of vacancy defects within the different 2D-TMD material
phases and vacancy concentration.

## Acknowledgments

The authors acknowledge financial support from the Brazilian agencies FAPESP
(grant 20/14067-3 and 17/02317-2), INCT-Nanomateriais de Carbono, and
Laboratório Nacional de Computação Científica for computer time.

## References

  * Liang  _et al._ [2021] Qijie Liang, Qian Zhang, Xiaoxu Zhao, Meizhuang Liu, and Andrew T. S. Wee, “Defect Engineering of Two-Dimensional Transition-Metal Dichalcogenides: Applications, Challenges, and Opportunities,” ACS Nano 15, 2165–2181 (2021).
  * Cavallini and Gentili [2022] Massimiliano Cavallini and Denis Gentili, “Atomic Vacancies in Transition Metal Dichalcogenides: Properties, Fabrication, and Limits,” ChemPlusChem 87, e202100562 (2022).
  * Lin  _et al._ [2017] X. Lin, J. C. Lu, Y. Shao, Y. Y. Zhang, X. Wu, J. B. Pan, L. Gao, S. Y. Zhu, K. Qian, Y. F. Zhang, D. L. Bao, L. F. Li, Y. Q. Wang, Z. L. Liu, J. T. Sun, T. Lei, C. Liu, J. O. Wang, K. Ibrahim, D. N. Leonard, W. Zhou, H. M. Guo, Y. L. Wang, S. X. Du, S. T. Pantelides, and H.-J. Gao, “Intrinsically patterned two-dimensional materials for selective adsorption of molecules and nanoclusters,” Nature Materials 16, 717–721 (2017).
  * Nguyen  _et al._ [2018] Giang D. Nguyen, Liangbo Liang, Qiang Zou, Mingming Fu, Akinola D. Oyedele, Bobby G. Sumpter, Zheng Liu, Zheng Gai, Kai Xiao, and An-Ping Li, “3D Imaging and Manipulation of Subsurface Selenium Vacancies in ${\mathrm{PdSe}}_{2}$,” Phys. Rev. Lett. 121, 086101 (2018).
  * Li  _et al._ [2015] Zhe Li, Hsin-Yi Tiffany Chen, Koen Schouteden, Koen Lauwaet, Ewald Janssens, Chris Van Haesendonck, Gianfranco Pacchioni, and Peter Lievens, “Lateral Manipulation of Atomic Vacancies in Ultrathin Insulating Films,” ACS Nano 9, 5318–5325 (2015).
  * Voiry  _et al._ [2016] Damien Voiry, Jieun Yang, and Manish Chhowalla, “Recent Strategies for Improving the Catalytic Activity of 2D TMD Nanosheets Toward the Hydrogen Evolution Reaction,” Advanced Materials 28, 6197–6206 (2016).
  * Tsai  _et al._ [2017] Charlie Tsai, Hong Li, Sangwook Park, Joonsuk Park, Hyun Soo Han, Jens K. Nørskov, Xiaolin Zheng, and Frank Abild-Pedersen, “Electrochemical generation of sulfur vacancies in the basal plane of MoS2 for hydrogen evolution,” Nature Communications 8, 15113 (2017).
  * Cheng  _et al._ [2014] Rui Cheng, Dehui Li, Hailong Zhou, Chen Wang, Anxiang Yin, Shan Jiang, Yuan Liu, Yu Chen, Yu Huang, and Xiangfeng Duan, “Electroluminescence and Photocurrent Generation from Atomically Sharp WSe2/MoS2 Heterojunction p–n Diodes,” Nano Letters 14, 5590–5597 (2014).
  * Cui  _et al._ [2015] Xu Cui, Gwan-Hyoung Lee, Young Duck Kim, Ghidewon Arefe, Pinshane Y. Huang, Chul-Ho Lee, Daniel A. Chenet, Xian Zhang, Lei Wang, Fan Ye, Filippo Pizzocchero, Bjarke S. Jessen, Kenji Watanabe, Takashi Taniguchi, David A. Muller, Tony Low, Philip Kim, and James Hone, “Multi-terminal transport measurements of MoS2 using a van der Waals heterostructure device platform,” Nature Nanotechnology 10, 534–540 (2015).
  * Wang  _et al._ [2012] Qing Hua Wang, Kourosh Kalantar-Zadeh, Andras Kis, Jonathan N. Coleman, and Michael S. Strano, “Electronics and optoelectronics of two-dimensional transition metal dichalcogenides,” Nature Nanotechnology 7, 699–712 (2012).
  * Avsar  _et al._ [2020a] A. Avsar, H. Ochoa, F. Guinea, B. Özyilmaz, B. J. van Wees, and I. J. Vera-Marun, “Colloquium: Spintronics in graphene and other two-dimensional materials,” Rev. Mod. Phys. 92, 021003 (2020a).
  * Avsar  _et al._ [2020b] Ahmet Avsar, Cheol-Yeon Cheon, Michele Pizzochero, Mukesh Tripathi, Alberto Ciarrocchi, Oleg V. Yazyev, and Andras Kis, “Probing magnetism in atomically thin semiconducting PtSe2,” Nature Communications 11, 4806 (2020b).
  * Du  _et al._ [2010] Guodong Du, Zaiping Guo, Shiquan Wang, Rong Zeng, Zhixin Chen, and Huakun Liu, “Superior stability and high capacity of restacked molybdenum disulfide as anode material for lithium ion batteries,” Chem. Commun. 46, 1106–1108 (2010).
  * Chang and Chen [2011] Kun Chang and Weixiang Chen, “l-Cysteine-Assisted Synthesis of Layered MoS2/Graphene Composites with Excellent Electrochemical Performances for Lithium Ion Batteries,” ACS Nano 5, 4720–4728 (2011).
  * Edelberg  _et al._ [2019] Drew Edelberg, Daniel Rhodes, Alexander Kerelsky, Bumho Kim, Jue Wang, Amirali Zangiabadi, Chanul Kim, Antony Abhinandan, Jenny Ardelean, Micheal Scully, Declan Scullion, Lior Embon, Rui Zu, Elton J. G. Santos, Luis Balicas, Chris Marianetti, Katayun Barmak, Xiaoyang Zhu, James Hone, and Abhay N. Pasupathy, “Approaching the Intrinsic Limit in Transition Metal Diselenides via Point Defect Control,” Nano Letters 19, 4371–4379 (2019).
  * Cai  _et al._ [2015] Liang Cai, Jingfu He, Qinghua Liu, Tao Yao, Lin Chen, Wensheng Yan, Fengchun Hu, Yong Jiang, Yidong Zhao, Tiandou Hu, Zhihu Sun, and Shiqiang Wei, “Vacancy-Induced Ferromagnetism of MoS2 Nanosheets,” Journal of the American Chemical Society 137, 2622–2627 (2015).
  * Bai  _et al._ [2018] Chunxu Bai, Yanling Yang, and Lin Bai, “Giant magnetoresistance control and nontrivial metallic state manipulation in a transition-metal dichalcogenide spin-valve using a gate voltage,” Journal of Physics: Condensed Matter 30, 495801 (2018).
  * Crasto de Lima and Fazzio [2021] Felipe Crasto de Lima and Adalberto Fazzio, “At the Verge of Topology: Vacancy-Driven Quantum Spin Hall in Trivial Insulators,” Nano Letters 21, 9398–9402 (2021).
  * Qiu  _et al._ [2013] Hao Qiu, Tao Xu, Zilu Wang, Wei Ren, Haiyan Nan, Zhenhua Ni, Qian Chen, Shijun Yuan, Feng Miao, Fengqi Song, Gen Long, Yi Shi, Litao Sun, Jinlan Wang, and Xinran Wang, “Hopping transport through defect-induced localized states in molybdenum disulphide,” Nature Communications 4, 2642 (2013).
  * Absor  _et al._ [2017] Moh. Adhib Ulil Absor, Iman Santoso, Harsojo, Kamsul Abraha, Fumiyuki Ishii, and Mineo Saito, “Defect-induced large spin-orbit splitting in monolayer ${\mathrm{PtSe}}_{2}$,” Phys. Rev. B 96, 115128 (2017).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, “Inhomogeneous Electron Gas,” Phys. Rev. 136, B864–B871 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, “Self-Consistent Equations Including Exchange and Correlation Effects,” Phys. Rev. 140, A1133–A1138 (1965).
  * Perdew  _et al._ [1996] John P. Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett. 77, 3865–3868 (1996).
  * Grimme  _et al._ [2010] Stefan Grimme, Jens Antony, Stephan Ehrlich, and Helge Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu,” The Journal of Chemical Physics 132, 154104 (2010).
  * Grimme  _et al._ [2011] Stefan Grimme, Stephan Ehrlich, and Lars Goerigk, “Effect of the damping function in dispersion corrected density functional theory,” Journal of Computational Chemistry 32, 1456–1465 (2011).
  * Blöchl [1994] P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953–17979 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, “From ultrasoft pseudopotentials to the projector augmented-wave method,” Phys. Rev. B 59, 1758–1775 (1999).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, “Ab initio molecular dynamics for open-shell transition metals,” Phys. Rev. B 48, 13115–13118 (1993).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, “Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set,” Phys. Rev. B 54, 11169–11186 (1996).
  * Focassio  _et al._ [2021] Bruno Focassio, Gabriel R. Schleder, Felipe Crasto de Lima, Caio Lewenkopf, and Adalberto Fazzio, “Amorphous ${\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$ structural, electronic, and topological nature from first principles,” Phys. Rev. B 104, 214206 (2021).


# Model-Based Reinforcement Learning Framework of Online Network Resource
Allocation

Bahador Bakhshi, Josep Mangues-Bafalluy This work has been partially funded by
the EC H2020 5Growth Project (grant no. 856709), and Generalitat de Catalunya
grant 2017 SGR 1195. Centre Tecnologic de Telecomunicacions de Catalunya
(CTTC/CERCA), Spain
{bbakhshi<EMAIL_ADDRESS>
###### Abstract

Online Network Resource Allocation (ONRA) for service provisioning is a
fundamental problem in communication networks. As a sequential decision-making
under uncertainty problem, it is promising to approach ONRA via Reinforcement
Learning (RL). But, RL solutions suffer from the _sample complexity_ issue;
i.e., a large number of interactions with the environment needed to find an
efficient policy. This is a barrier to utilize RL for ONRA as on one hand, it
is not practical to train the RL agent offline due to lack of information
about future requests, and on the other hand, online training in the real
network leads to significant performance loss because of the sub-optimal
policy during the prolonged learning time. This performance degradation is
even higher in non-stationary ONRA where the agent should continually adapt
the policy with the changes in service requests. To deal with this issue, we
develop a general resource allocation framework, named RADAR, using model-
based RL for a class of ONRA problems with the known immediate reward of each
action. RADAR improves sample efficiency via exploring the state space in the
background and exploiting the policy in the decision-time using synthetic
samples by the _model_ of the environment, which is trained by real
interactions. Applying RADAR on the multi-domain service federation problem,
to maximize profit via selecting proper domains for service requests
deployment, shows its continual learning capability and up to 44% performance
improvement w.r.t. the standard model-free RL solution.

###### Index Terms:

Online Resource Allocation, Model-based Reinforcement Learning, Service
Federation, Continual Learning

##  I Introduction

Resource allocation is a fundamental problem in communication networks [1, 2],
wherein network resources are allocated for service requests to optimize an
objective function, e.g., the service provider’s profit, subject to satisfying
Quality of Service (QoS) requirements. Resource allocation problems are
categorized into two categories, namely online and offline. In the former,
service requests arrive one-by-one over time, and the Resource Allocator (RA)
allocates resources for each request without knowledge of future ones. But in
offline problems, it is assumed that all requests are given/known in advance.
In many practical use cases, where the exact information of future requests is
not available, online RA solutions are the de facto approach.

Due to the central role of the resource allocation problem, various solutions
have been developed over the last decades including problem-specific heuristic
RA algorithms, e.g., QoS routing [3], game theory and meta-heuristic based
solutions [4], and various optimization theory techniques [5]. Recently, by
the emergence of Machine Learning (ML) and its successful applications, ML-
based solutions for network resource allocation and management have also been
proposed [6]. Reinforcement Learning (RL) is a class of ML solutions that
contrary to supervised learning methods, learns the suitable action in each
state, known as policy, via the rewards it receives from interactions with the
environment [7]. RL is an efficient tool to approach a class of problems known
as _sequential decision making under uncertainty_ (SDMU), where the decision-
maker needs to take actions without knowing the uncertain future. In applying
RL to SDMU, the decision-maker is an RL agent that in a trial-and-error
process explores the state space, and eventually finds an appropriate policy.

Online Network Resource Allocation (ONRA) problems are indeed instances of
SDMU as the RA should allocate resources while the future requests are
uncertain. So, RL is a promising approach for ONRA problems [8, 9, 10].
Although the main advantage of RL is that it does not need explicitly labeled
training data, to find an efficient policy, the agent requires a huge number
of interactions with the environment, which is known as the _sample
complexity_ issue [11]; for example, in the well-known deep RL solution for
playing Atari games [12], the agent is trained using $5\times 10^{7}$ frames
of the game, or in our proposed solution for multi-domain service federation
(MDFS) [8], it takes about 9 $\times 10^{5}$ interactions to find the near-
optimal policy. In applying RL to ONRA, where the agent learns the policy in
an online manner in the real network, sample complexity can cause significant
performance degradation because before finding the (near) optimal policy,
agent’s decisions are sub-optimal for a high number of real requests. Dealing
with the RL sample complexity in ONRA is the research gap that we aim to
address in this paper.

_Model-Based_ RL (MBRL) is an approach to improve sample efficiency of RL by
integrating a model of the environment in the RL [11]. The core idea is that
in each interaction with the real environment, besides improving the policy,
the agent also learns a model of the environment. Then it generates synthetic
interactions by this model which are used to improve the policy. MBRL has been
used in contexts like playing Atari games [13], and robotic [14]; recently, it
has also got attention in communication networks. In [15], an MBRL approach
was developed for microservice-based applications resource allocation. Dynamic
computational resource allocation in the cloud was studied in [16], where the
problem is formulated as a Markov Decision Process (MDP) and by the value
iteration algorithm, an MBRL method was developed. In [17], again resource
allocation in the cloud was investigated, the authors developed an MBRL
approach to satisfy the application constraints on the rate of allocation
changes. These solutions are not applicable to ONRA as they are specifically
designed for cloud environment to model computational resources and workloads.
In this paper, we develop a _general framework_ for ONRA using MBRL, and make
the following contributions:

  * •

a general framework, called RADAR, with the capability of continual learning
is developed for ONRA problems;

  * •

by exploiting the flexibility of the framework, four algorithms are
implemented to be used in different use cases;

  * •

as a proof of concept, the RADAR-based algorithms are applied and evaluated in
the MDSF problem.

The remainder of this paper is organized as follows. In Section II, the system
model and the sample complexity issue are elaborated in more detail. The RADAR
framework is presented in Section III. The MDSF problem is solved by RADAR in
Section IV, which is numerically evaluated in Section V. Finally, Section VI
concludes this paper.

##  II System Model and Problem Statement

In this paper, we study a class of ONRA problems, named Known-Reward Online
Resource Allocation (KRORA); which is depicted in Fig. 1, and identified by
the following characteristics. First, the problem is online; i.e., requests
for the services defined in the service catalog arrive one-by-one. For each
request, without knowledge of future requests, the RA algorithm runs, and
after time $t_{1}$, takes a resource allocation action. Second, each RA action
has a _known_ immediate reward, as in practice the service provider either
knows or can determine the immediate consequences of each RA action. This is
exemplified in the following paragraph. Third, RA’s objective is to optimize a
cumulative long-term reward. Fourth, the next request arrives after time
$t_{2}$. Meanwhile, during this period, some alive services terminate and
depart the network. If an arrival or departure occurs in the time interval
$t_{1}$, it will be taken into account as soon as the RA algorithm ends.

A large number of practical ONRA problems are indeed in the KRORA class as
they have the mentioned characteristics. For example, in the online flow
routing problem, each request is a demand for a flow from a source node to a
destination node with known predefined QoS requirements. The objective is to
maximize the total number of accepted flows; so, the immediate reward is 1.
Another example is the MDSF problem, which is elaborated in Section IV.

As a subclass of ONRA, KRORA problems can be modeled as instances of SDMU. MDP
is a framework to formulate SDMU as a tuple
$(\mathcal{S},\mathcal{A},\mathfrak{P},\mathfrak{R})$, where $\mathcal{S}$ is
the environment states set, $\mathcal{A}=\\{\mathcal{A}(s)\ \forall
s\in\mathcal{S}\\}$, where $\mathcal{A}(s)$ is the set of the actions in state
$s$. Let $\mathcal{S}^{\prime}_{s,a}$ be the set of possible next states in
the case of taking action $a$ in state $s$,
$\mathfrak{P}(s,a,s^{\prime})\\!\\!:\mathcal{S}\times\mathcal{A}(s)\times\mathcal{S}^{\prime}_{s,a}\rightarrow[0,1]$
is the probability of the transition from $s$ to
$s^{\prime}\in\mathcal{S}^{\prime}_{s,a}$, and
$\mathfrak{R}(s,a)\\!\\!:\mathcal{S}\times\mathcal{A}(s)\rightarrow\mathbb{R}$
is the reward function. In KRORA, the service provider offers
$\mathcal{I}=\\{1,\ldots,I\\}$ types of services, which are defined in the
service catalog. Each state $s$, as shown in Fig. 1, corresponds to arrival of
a service request wherein the RA takes an action. State is defined as
$s=(\bm{c},\bm{l},\bm{d})$ where vector $\bm{c}$ represents the available
network resources, $\bm{l}$ is the vector representing the active requests of
each service type $i\in\mathcal{I}$, and $\bm{d}$ is a 0/1 vector identifying
the service type of the current request. $\mathcal{A}(s)$ is the set of
feasible different resource allocation patterns, e.g., different feasible
paths for the request. $\mathfrak{R}(s,a)$ is the reward of each RA pattern,
which is already known or can be determined. However,
$\mathfrak{P}(s,a,s^{\prime})$, as it is shown in Fig. 1, is determined by the
arrival and departure rates of the request during period $t_{2}$, which are
not known in advance. An example of such formulation is presented in Section
IV.

Figure 1: The timeline of the KRORA problems

RL is a common approach to solve MDP. In each state $s$, the agent takes
action $a$ and, consequently, the environment determines the next state
$s^{\prime}$ and reward $\mathfrak{R}(s,a)$. The fundamental concept in RL is
Temporal Difference (TD), i.e., the difference between the current value of
the state-action pair $Q[s,a]$ and the new estimation obtained according to
this interaction. This difference is used in a TD-update rule to improve the
policy, e.g., in the Q-Learning, the TD-update rule is

$Q[s,a]=(1-\alpha)Q[s,a]+\alpha\Big{(}\mathfrak{R}(s,a)+\gamma\max_{a^{\prime}}Q[s^{\prime},a^{\prime}]\Big{)},$

where $\alpha$ is the learning rate and $\gamma$ is the discount factor [7].

Despite successful applications of LR in contexts like games and robotics,
there are concerns about the practicality of applying RL to ONRA due to the
sample complexity issue. The RL agent can be trained either _offline_ or
_online_. In the offline approach, learning and evaluating are separated
stages, where in the first stage, the agent learns the policy by exploring
state space via a large number of trial-and-error interactions with the
(simulated) environment but it is not evaluated; in the second stage, the
learned policy is evaluated without learning. In the online approach, the
agent should learn the policy in the real environment, so, it is evaluated
while it is learning.

The offline approach is applicable only if there is an isolated training
environment that is not practical in KRORA due to the lack of information
about future requests. Moreover, in the offline approach, the agent should be
re-trained if the environment changes, which is quite often in KRORA due to
changes in request rates111Please note that the approach composed of these
steps: (1) learning the stochastic behavior of requests, (2) developing a
simulator using the stochastic models, (3) training the agent offline in the
simulated environment, (4) using the trained agent/policy in the real
environment; is not a practical solution either, as it again needs
(considerable) data for learning the stochastic behavior and also re-training
the agent in non-stationary environments.. So, the online training approach is
more suitable in KRORA but due to the sample complexity issue, the agent’s
decisions are sub-optimal for a high number of requests that decrease the
cumulative reward significantly. When request rates are non-stationary, the
agent should _continually learn_ the policy, which makes the problem worse
because, if dynamics of the environment change faster than the time it takes
to find the optimal policy, then the policy is always sub-optimal.

In summary, online agent training is the practical approach to deal with KRORA
via RL, but a solution is needed to alleviate the prolonged learning period.
In the following sections, we propose the RADAR (Resource Allocation via moDel
leARning) framework built on model-based RL to improve sample efficiency of
RL.

##  III MBRL Framework for KRORA

###  III-A Design Concepts

In traditional model-free reinforcement learning, like Q-Learning and DQN, it
is assumed that no information about the dynamics of the environment is known
in advance. Thus, the agent has to learn the policy by a TD-update rule in a
high number of real interactions with the environment, which leads to the
sample complexity problem. While this is a reasonable assumption in some
problems, we can do much better in KRORA, as some information about the
environment is known in advance, and some information can also be learned via
the real interactions. This is the design idea of the RADAR framework built on
MBRL—in each real interaction with the environment, the agent updates the
policy by a given TD-update rule exactly in the same way that model-free RL
does; but moreover, it also learns a _model_ of the environment, and then
exploits it to update the policy furthermore using samples generated by the
model instead of relying only on the real interactions with the environment.
So, it can efficiently deal with the sample complexity issue by generating
(guided) synthetic samples. However, to utilize this approach, a few
fundamental questions should be answered: ($i$) what is the model of the
environment? ($ii$) how to learn the model? ($iii$) how to integrate the model
into the RL architecture? In the following, we answer these questions in
RADAR.

In general, developing an analytical model for a (stochastic) environment
that, for given $s$ and $a$, determines $\mathfrak{P}(s,a,s^{\prime})$
$\forall s^{\prime}\in\mathcal{S}^{\prime}_{s,a}$ is not an easy task.
Therefore, in RADAR, we resort to a _sample_ model that in each invocation,
for given $s$ and $a$, returns just a
$s^{\prime}\in\mathcal{S}^{\prime}_{s,a}$ and $\mathfrak{R}(s,a)$. Referring
back to Fig. 1, the steps of this model are as follows: ($i$) action $a$ is
applied in state $s=(\bm{c},\bm{l},\bm{d})$ that changes $\bm{c}$ and $\bm{l}$
to $\tilde{\bm{c}}$ and $\tilde{\bm{l}}$, ($ii$) the reward
$\mathfrak{R}(s,a)$ is determined, ($iii$) sequentially, samples are taken
from the arrival and departure stochastic processes. If the sample is
departure, it is applied in the network, it changes $\tilde{\bm{c}}$ and
$\tilde{\bm{l}}$ to $\bm{c}^{\prime}$ and $\bm{l}^{\prime}$, and another
sample is taken; however, if the sample is arrival, it determines the next
state $s^{\prime}=(\bm{c}^{\prime},\bm{l}^{\prime},\bm{d}^{\prime})$.

In KRORA, as the service provider knows the available network resources
$\bm{c}$ and alive services $\bm{l}$, so can it determine the effect of each
RA action for the given request for service type $i$; i.e., it can perform
step ($i$). Moreover in KRORA, it is assumed that the immediate reward is
known; so, step ($ii$) is straightforward. The only remaining part is the
requests arrival/departure _stochastic processes_ , which are not known in
advance. In the RADAR framework, these processes are _learned_ from the
interactions with the real environment. There are various approaches to learn
a stochastic process from data [18, 19]. In RADAR, it is a black-box which on
one hand, takes information about the arrival/departure times of real
requests, and internally learns a stochastic process for it; on the other
hand, it can generate samples of the stochastic process. So, until now, we
answered the first two questions—the model is a sample model and model
learning is indeed learning the requests arrival/departure stochastic
processes.

The last issue is the integration of the sample model in the RL architecture
to improve sample efficiency by updating the policy using the synthetic
interactions by the model. In RADAR, this can be conducted in both time
intervals $t_{1}$ and/or $t_{2}$ depicted in Fig. 1. In $t_{1}$, the request
is known, so we seek the best action via _decision-time planning_ composed of
exploration and exploitation steps, but in $t_{2}$ the action has already been
taken, so we aim to improve the policy for the next unknown requests by
_background exploration_ of the MDP.

To explore state space, either in background or decision-time, we start from a
given state and update the policy by the RL TD-update rule in a number of
trajectories of synthetic samples where, to emphasize exploration, random
actions are taken. In exploitation, we update the value of each action in the
given state by generating a number of synthetic samples, but to exploit the
policy, we greedily select the action with the highest value. Details are
elaborated in the following.

###  III-B RADAR Framework Architecture

The architecture of RADAR built on the explained solutions, is shown in Fig.
2. The user’s request in combination with the available resources of the
network and alive services are given as a new state to the framework. The
arrival time is used to update the model of the arrival stochastic process.
The state itself is processed by the decision-time planning module. The
components of the module utilize the sample model and update the policy. Then,
the state is given to the RL algorithm, which selects the best action
according to its exploration strategy. The action is given to the environment
and, after a while, when the new state is given, the policy will be updated by
the RL TD-update rule. During this time, i.e., $t_{2}$ in Fig. 1, the state
and the selected action are processed by the background exploration module,
which also uses the sample model to update the policy. Finally, when a service
departs the network, its departure time is used to update the departure
stochastic model. The sample model module, besides the arrival and departure
models, uses the service catalog to determine the $\mathfrak{R}(s,a)$ and
$s^{\prime}$ for each given $s$ and $a$.

Figure 2: The architecture of the RADAR framework

RADAR is quite flexible; according to the constraints on $t_{1}$ and $t_{2}$,
not only the number and length of the exploration trajectories can be adjusted
per module but also each module can be enabled/disabled independently. Table I
shows four RL algorithms (out of 8 possible options) as representative
solutions based on RADAR to be used in different practical situations. MFRL is
the standard model-free RL where all modules are disabled. MB-BGEX is the
algorithm that can be used in real-time KRORA problems where $t_{1}=0$. If the
state instantly changes after taking action, i.e., $t_{2}=0$, the MB-DTP can
be used. Finally, MB-Full can be used when $t_{1},t_{2}>0$.

Table I: Algorithms based on the RADAR Framework

Algorithm | BG | Decision-Time |  Comments
---|---|---|---
| Explore | Explore | Exploit |
MFRL | $\times$ | $\times$ | $\times$ | Standard model-free RL
MB-BGEX | $\checkmark$ | $\times$ | $\times$ |  Real-time RA ($t_{1}=0$)
MB-DTP | $\times$ | $\checkmark$ | $\checkmark$ |  Instance next state ($t_{2}=0$)
MB-Full | $\checkmark$ | $\checkmark$ | $\checkmark$ | $t_{1},t_{2}>0$

###  III-C RADAR Framework Implementation

The main function of the sample model module is the step procedure, depicted
in Algorithm 1, that returns $s^{\prime}$ and $\mathfrak{R}(s,a)$ for a given
$s$ and $a$. It first looks up the new request in the service catalog and
finds the associated reward and required resources. Then, it applies the given
action that updates $\bm{c}$ and $\bm{l}$. In lines 5 and 6, the next request
is determined by sampling the arrival stochastic processes. Finally, according
to the departure stochastic process samples, if an active service terminates
before the next arrival, its resources are released.

Exploration, either in background time $t_{2}$ or decision-time $t_{1}$, is
conducted by the Explore procedure shown in Algorithm 2. It explores the state
space in $\theta$ trajectories starting from state $s_{0}$ with length
$\kappa$. In the beginning, the sample model is initialized by the state
$s_{0}$. Then, in each trajectory, if the action in state $s_{0}$ is already
known, the policy is fixed in this state in line 3. As the objective is to
explore the state space, a random action is selected in each state. The next
state and the corresponding reward are given by the sample model; and finally,
in line 8, the policy is updated. The background exploration in state $s$,
where the action $a$ has already been taken, is implemented by invoking
Explore$(s,a,\theta,\kappa)$; and the decision-time exploration in state $s$
is indeed Explore$(s$, NULL, $\theta,\kappa)$ as the action in this state is
not known in $t_{1}$.

Algorithm 1 SampleModelStep$(s=(\bm{c},\bm{l},\bm{d}),a)$

1:$\mathfrak{R}(s,a)\leftarrow$ lookup reward of new request $\bm{d}$ in
service catalog

2:$\bm{u}\leftarrow$ lookup required resources of request $\bm{d}$ in service
catalog

3:$\tilde{\bm{c}}\leftarrow$ apply action $a$ and resources $\bm{u}$ on
$\bm{c}$

4:$\tilde{\bm{l}}\leftarrow$ apply action $a$ on $\bm{l}$

5:$\mathcal{T}=\\{t_{a,1},\ldots t_{a,I}\\}\leftarrow$ sample of arrival
process $\forall i\in\mathcal{I}$

6:$\bm{d}^{\prime}\leftarrow\operatorname*{argmin}\mathcal{T}$,
$t_{a}\leftarrow\min\mathcal{T}$

7:for $j\in\mathcal{I}$ s.t. $d_{j}>0$ do

8: $t_{d,j}\leftarrow$ sample of departure stochastic process $j$

9: Release request $j$, update $\tilde{\bm{l}}$ and $\tilde{\bm{c}}$ if
$t_{d,j}<t_{a}$

10:$s^{\prime}=(\tilde{\bm{c}},\tilde{\bm{l}},\bm{d}^{\prime})$

11:return $s^{\prime}$, $\mathfrak{R}(s,a)$

Algorithm 2 Explore$(s_{0},a_{0},\theta,\kappa)$

1:for $\theta$ times do

2: Initialize the sample model by the given state $s_{0}$

3: $\mathcal{A}(s_{0})\leftarrow\\{a_{0}\\}$ if $a_{0}\neq$ NULL

4: $s\leftarrow s_{0}$

5: for $\kappa$ times do

6: $a\leftarrow$ random action of $\mathcal{A}(s)$

7: $s^{\prime},\mathfrak{R}(s,a)\leftarrow$ SampleModelStep$(s,a)$

8: TD-update$(s,a,s^{\prime},\mathfrak{R}(s,a))$

9: $s^{\prime}\leftarrow s$

Algorithm 3 Exploit$(s_{0},\theta,\kappa)$

1:for $\bar{a}\in\mathcal{A}(s_{0})$ do

2: for $\theta$ times do

3: Initialize the sample model by the given state $s_{0}$

4: $\mathcal{A}(s_{0})\leftarrow\\{\bar{a}\\}$, $s\leftarrow s_{0}$

5: for $\kappa$ times do

6: $a\leftarrow$ the highest value action in $\mathcal{A}(s)$

7: $s^{\prime},\mathfrak{R}(s,a)\leftarrow$ SampleModelStep$(s,a)$

8: stack.push$(s,a,s^{\prime},\mathfrak{R}(s,a))$

9: $s^{\prime}\leftarrow s$

10: while stack is not empty do

11: $s,a,s^{\prime},\mathfrak{R}(s,a)\leftarrow$ stack.pop()

12: TD-update$(s,a,s^{\prime},\mathfrak{R}(s,a))$

In the decision-time exploitation, we aim to evaluate the quality of each
action, so we generate a number of trajectories per action in the state to
update the value of the action according to the value of the next states.
Consider a trajectory $(s_{0},a_{0})$ $\rightarrow$ $(s_{1},a_{1})$
$\rightarrow$ $\ldots$ $\rightarrow$ $(s_{\kappa},a_{\kappa})$. In the
standard operation of RL algorithms, value $(s_{0},a_{0})$ is updated
according to value $s_{1}$ and _then_ the value of $(s_{1},a_{1})$ is updated.
So, action value updates in states $s_{1}$ …$s_{\kappa}$ do not have any
effect on the value $(s_{0},a_{0})$ in this trajectory. To take these updates
into account, in the decision-time exploration, we update the values
_backward_ , i.e., from $(s_{\kappa},a_{\kappa})$ to $(s_{0},a_{0})$. The
details of the implementation is shown in Algorithm 3. To exploit the current
policy, the action with the highest value is selected in line 6; and the
backward updates are implemented using a stack.

##  IV Multi-Domain Service Federation

In this section, as a proof of concept, we apply RADAR on MDSF [8] that
enables the service provider to collaborate with other providers in service
provisioning where the federation contract provides extra resources for the
consumer domain at the cost of federation. In MDSF, the RA either determines
the domain to deploy the given request or rejects it while it is not aware of
future requests; i.e., $\mathcal{A}(s)$={reject, local deploy, federate}.
Here, we consider MDSF between a consumer domain and a provider domain. Each
service type $i\in\mathcal{I}$ is specified by a tuple $(w_{i},r_{i})$ in the
service catalog, where $w_{i}$ is the total amount of required resources e.g.,
CPU cores, and $r_{i}$ is the revenue if a request of type $i$ is accepted;
moreover, according to the federation contract, $\varphi_{i}$ is the cost of
deploying request $\delta_{i}$ of type $i$ in the provider domain. So,
$\mathfrak{R}(s,a)$ of the reject, local deploy, and federate actions are
respectively known as 0, $r_{i}$, $r_{i}-\varphi_{i}$. The objective is to
maximizes the service provider’s average profit

$\frac{1}{|\mathcal{D}|}\sum_{i\in\mathcal{I}}\Big{(}\sum_{\delta_{i}\in\mathcal{L}}r_{i}+\sum_{\delta_{i}\in\mathcal{F}}(r_{i}-\varphi_{i})\Big{)};$
(1)

where, $\mathcal{D}$ is the set of requests; and $\mathcal{L}$ ($\mathcal{F}$)
is the set of the services deployed in the consumer (provider) domain. It is
easy to show that MDSF is a KRORA problem.

In MDSF, the state is defined as
$s=(\bm{c}_{c},\bm{c}_{p},\bm{l}_{c},\bm{l}_{p},\bm{d})$, where $\bm{c}_{c}$
and $\bm{c}_{p}$ are respectively the current available capacity of the
consumer and provider domains. $\bm{l}_{c}$ and $\bm{l}_{p}$ are vectors
wherein the $i$-th element is the number of alive services of type $i$ in the
consumer domain and in the provider domain, respectively.

To solve MDSF through RADAR, we need to specify the RL algorithm and the
arrival/departure stochastic process models; the remaining components of the
framework are problem-independent. In [8], we showed that average-reward
reinforcement learning, i.e., the R-Learning algorithm [20], is an efficient
solution for the problem. In this algorithm, in addition to the values of the
actions $Q[s,a]$, a parameter $\rho$, which is the average reward of the MDP,
is also learned. So the TD-update rule of the RL algorithm is:

$Q[s,a]\leftarrow(1-\alpha)Q[s,a]+\alpha\Big{(}\mathfrak{R}(s,a)-\rho+\max_{a^{\prime}}Q[s^{\prime},a^{\prime}]\Big{)},$

and

$\rho\leftarrow(1-\beta)\rho+\beta\Big{(}\mathfrak{R}(s,a)-\max_{a}Q[s,a]+\max_{a^{\prime}}Q[s^{\prime},a^{\prime}]\Big{)}.$

For the arrival and departure stochastic processes, we assumed that they are
Poisson processes, as usually done in the literature. Under this assumption,
we only need to learn the expected value of the distribution. Here, we use
exponential moving average to estimate it. It must be noted that any other
method can also be used to learn stochastic processes. In section V-D, we
evaluate the effect of this assumption.

##  V Simulation Results

In this section, we evaluate the RADAR-based RA algorithms, shown in Table I,
in the MDSF problem. As the _theoretical_ performance bound, the results of
the offline R-Learning algorithm, shown is near-optimal [8], are also included
by legend RL-Offline. In RL-Offline, the agent is trained offline using
$\theta$ = $10^{3}$ trajectories with $\kappa$ = $10^{4}$ requests per
trajectory. The performance metric is the average profit defined in (1), and
the capacity of the consumer and provider domains are respectively 100 and 50
CPU cores. The default settings of the simulated service types are shown in
Table II. The following results are the average of 20 experiments.

###  V-A Exploration Breadth vs. Depth

In the RADAR’s modules, the parameters $\theta$ and $\kappa$ determine the
breadth and depth of exploration trajectories. The effects of these parameters
on the framework modules’ performance are shown in Fig. 3 (the result of the
decision-time exploration is similar; omitted due to space limitations). While
evaluating each module, the other modules are disabled. It is seen that
$\theta$ and $\kappa$ should be larger than a threshold to effectively explore
the state space. The background exploration needs larger $\theta$ as due to
the unknown next state, it needs more trajectories to visit possible next
states. Increasing the values of the parameters, which increases the solution
complexity, beyond certain values does not yield better performance. Based on
these results, in the following simulations, we set $\theta$ = 5 and $\kappa$
= 3 in background exploration, $\theta$ = 3 and $\kappa$ = 2 in decision-time
exploration, and $\theta$ = 1, $\kappa$ = 3 in decision-time exploitation.

###  V-B Learning Capability

In this section, we evaluate how fast and efficient the RADAR-based algorithms
find the policy utilizing the sample model. For this assessment, learning of
the sample model and also policy update are disabled after a given time period
from simulation start time. Fig. 6 shows the performance of the algorithms
with respect to the length of the learning period (in terms of percentage of
simulation time). For the learning period less than 30% of the simulation
time, all algorithms are almost the same, and model-based RL does not perform
better than model-free RL because, in this case, the algorithms do not have
sufficient opportunity to explore the huge state space of the MDP, i.e.,
$|\mathcal{S}|=O(10^{9})$. However, by providing more opportunities for
learning, the MBRL algorithms learn a better policy faster, which leads to a
higher average profit than MFRL. These results also show the sample efficiency
of RADAR; MBRL-Full, which is learning during the whole simulation time using
about 16900 real samples, achieves an 11% optimality gap compared to RL-
Offline, which was trained with $10^{7}$ samples; i.e., 590 times fewer
samples.

Table II: Simulated Service Type Parameters

$i$ |  $w_{i}$ (core) |  $r_{i}$ ($) |  $\varphi_{i}$ ($) |  $\lambda_{i}$ (req $/$ h) |  $\mu_{i}$ (req $/$ h)
---|---|---|---|---|---
1 | 2 | 100 | 30 | 10 | 0.4
2 | 1 | 20 | 5 | 5 | 0.05
3 | 3 | 50 | 45 | 2 | 0.2

(a) Background Exploration (b) Decision-Time Exploitation

Figure 3: The effect of $\theta$ and $\kappa$ on the performance of each
modules

Figure 4: Effect of length of learning interval

Figure 5: Effect of variable traffic

Figure 6: Effect of inaccurate traffic model

###  V-C Adaptiveness Capability

As discussed, one of the advantages of MBRL is its continual learning
capability to adapt the policy to the changes in the environment. In this
section, we evaluate the performance of the RADAR-based algorithms in the case
that arrival rates of requests are changing over time. For this purpose, the
simulation time is divided into 5 equal time slots, and if a service type is
variable, then $\lambda_{1}=[6,8,10,8,6]$, $\lambda_{2}=[3,1,3,1,3]$, and
$\lambda_{3}=[4,2,0,2,4]$ req/h during the time slots. Fig. 6 shows the
average profit of the algorithms w.r.t. the number of variable service types.
As it is seen, MB-Full considerably outperforms MFRL that means RADAR can
improve the performance up to 44% in comparison to the standard model-free RL.
Moreover, the gap between MB-Full, using about $4\times 10^{5}$ real samples,
and RL-Offline, trained offline with $10^{7}$ samples, is independent of the
number of the variable service types, hence showing the excellent sample
efficiency of RADAR and its capability to continually learn the sample model
and use it to adapt the policy with the changes of the environment.

###  V-D Model Inaccuracy Tolerance

As mentioned, to apply RADAR on MDSF, we assumed that the arrival and
departure of requests are Poisson processes. In this section, we evaluate the
capability of the RADAR-based solutions in tolerating the inaccuracy of the
model. The average profits of the algorithms w.r.t. different stochastic
processes are shown in Fig. 6. In the Poisson process, the Inter-Arrival Time
(IAT) and Holding Time (HT) are respectively exponential distributions with
mean $1/\lambda_{i}$ and $1/\mu_{i}$, shown in Table II. In the case of
uniform distribution, IAT and HT are uniformly distributed in intervals
$[0,2/\lambda_{i}]$ and $[0,2/\mu_{i}]$, and for the normal distribution, IAT
= $\max(0,\mathcal{N}(10/\lambda_{i},(5/\lambda_{i})^{2})$ and HT =
$\max(0,\mathcal{N}(10/\lambda_{i},(5/\lambda_{i})^{2})$. These results show
that modeling the normal and uniform distributions as a Poisson process does
not significantly affect the performance; i.e., the sample model learned by
RADAR and the corresponding synthetic samples can remarkably improve the
performance.

##  VI Conclusion

In this paper, for KRORA problems, where the immediate rewards of RA’s actions
are known, we developed the RADAR framework built on MBRL, which is capable of
continual learning. In combination with the information about the network
resources, which the service provider has, the known rewards are utilized to
develop a sample model of the environment. In this model, the arrival and
departure stochastic processes of service requests are learned via
interactions with the environment. The model is integrated with an RL
algorithm to generate synthetic samples used in both background and decision-
time for updating the policy. Evaluation of RADAR in the MDSF problem showed
that it is not only applicable to KRORA problems, but also it can improve the
average profit up to 44% compared to the standard model-free RL.

## References

  * [1] Y. Xu, G. Gui, H. Gacanin, and F. Adachi, “A survey on resource allocation for 5g heterogeneous networks: current research, future trends and challenges,” IEEE Communications Surveys & Tutorials, 2021.
  * [2] S. Yang, F. Li, S. Trajanovski, R. Yahyapour, and X. Fu, “Recent advances of resource allocation in network function virtualization,” IEEE Transactions on Parallel and Distributed Systems, vol. 32, no. 2, pp. 295–314, 2020.
  * [3] J. W. Guck, A. Van Bemten, M. Reisslein, and W. Kellerer, “Unicast qos routing algorithms for sdn: A comprehensive survey and performance evaluation,” IEEE Communications Surveys & Tutorials, vol. 20, no. 1, pp. 388–415, 2017\.
  * [4] F. Rahmani, J. H. Joloudari, S. Shamshirband, and S. Mostafavi, “Game theory and evolutionary-optimization methods applied to resource allocation problems in emerging computing environments: A survey,” arXiv preprint arXiv:2012.11355, 2020.
  * [5] W. Ejaz, S. K. Sharma, M. Naeem, and N. A. Chughtai, “A comprehensive survey on resource allocation for cran in 5g and beyond networks,” Journal of Network and Computer Applications, 2020.
  * [6] M. E. Morocho-Cayamcela, H. Lee, and W. Lim, “Machine learning for 5g/b5g mobile and wireless communications: Potential, limitations, and future directions,” IEEE Access, vol. 7, pp. 137184–137206, 2019.
  * [7] R. S. Sutton and A. G. Barto, Reinforcement learning: An introduction.  MIT press, 2018.
  * [8] B. Bakhshi, J. Mangues-Bafalluy, and J. Baranda, “R-learning-based admission control for service federation in multi-domain 5g networks,” in IEEE GLOBECOM, 2021.
  * [9] X. Xiong, K. Zheng, and L. Hou, “Resource allocation based on deep reinforcement learning in iot edge computing,” IEEE Journal on Selected Areas in Communications, vol. 38, no. 6, pp. 1133–1146, 2020.
  * [10] H. Ye, G. Y. Li, and B.-H. F. Juang, “Deep reinforcement learning based resource allocation for v2v communications,” IEEE Transactions on Vehicular Technology, vol. 68, no. 4, pp. 3163–3173, 2019.
  * [11] T. M. Moerland, J. Broekens, and C. M. Jonker, “Model-based reinforcement learning: A survey,” arXiv preprint arXiv:2006.16712, 2020.
  * [12] V. Mnih and et al., “Human-level control through deep reinforcement learning,” nature, vol. 518, no. 7540, pp. 529–533, 2015.
  * [13] L. Kaiser, M. Babaeizadeh, P. Milos, B. Osinski, R. H. Campbell, K. Czechowski, D. Erhan, C. Finn, P. Kozakowski, S. Levine, et al., “Model-based reinforcement learning for atari,” arXiv preprint arXiv:1903.00374, 2019\.
  * [14] A. S. Polydoros and L. Nalpantidis, “Survey of model-based reinforcement learning: Applications on robotics,” Journal of Intelligent & Robotic Systems, vol. 86, no. 2, pp. 153–173, 2017.
  * [15] Z. Yang, P. Nguyen, H. Jin, and K. Nahrstedt, “Miras: Model-based reinforcement learning for microservice resource allocation over scientific workflows,” in IEEE International Conference on Distributed Computing Systems (ICDCS), pp. 122–132, IEEE, 2019.
  * [16] B. Chen, Y. Zhang, G. Iosifidis, and M. Liu, “Reinforcement learning on computational resource allocation of cloud-based wireless networks,” in IEEE WF-IoT, pp. 1–6, IEEE, 2020.
  * [17] Y. Jiang, M. Kodialam, S. Mukherjee, and L. Tassiulas, “Resource allocation in data centers using fast reinforcement learning algorithms,” IEEE Transactions on Network and Service Management, 2021.
  * [18] D. Koller and R. Fratkina, “Using learning for approximation in stochastic processes.,” in ICML, pp. 287–295, 1998.
  * [19] A. Y. Foong, W. P. Bruinsma, Y. Dubois, J. Requeima, and R. E. Turner, “Meta-learning stationary stochastic process prediction with convolutional neural processes,” arXiv preprint arXiv:2007.01332, 2020.
  * [20] A. Schwartz, “A reinforcement learning method for maximizing undiscounted rewards,” in International conference on machine learning, pp. 298–305, 1993.