# Unified Graph Neural Network Force-field for the Periodic Table

Kamal Choudhary (0000-0001-9737-8074<EMAIL_ADDRESS>[  Brian
DeCost (0000-0002-3459-5888)  [  Lily Major (0000-0002-5783-8432)  [  Keith
Butler (0000-0001-5432-5597)  [  Jeyan Thiyagalingam (0000-0002-2167-1343)  [
Francesca Tavazza (0000-0002-5602-180X)  [

###### Abstract

Classical force fields (FF) based on machine learning (ML) methods show great
potential for large scale simulations of materials. MLFFs have hitherto
largely been designed and fitted for specific systems and are not usually
transferable to chemistries beyond the specific training set. We develop a
unified atomisitic line graph neural network-based FF (ALIGNN-FF) that can
model both structurally and chemically diverse materials with any combination
of 89 elements from the periodic table. To train the ALIGNN-FF model, we use
the JARVIS-DFT dataset which contains around 75000 materials and 4 million
energy-force entries, out of which 307113 are used in the training. We
demonstrate the applicability of this method for fast optimization of atomic
structures in the crystallography open database and by predicting accurate
crystal structures using genetic algorithm for alloys.

###### keywords:

American Chemical Society, LaTeX

National Institute of Standards and Technology] Material Measurement
Laboratory, National Institute of Standards and Technology, Gaithersburg,
20899, MD, USA. [Theiss Research] Theiss Research, La Jolla, 92037, CA, USA.
National Institute of Standards and Technology] Materials Measurement
Laboratory, National Institute of Standards and Technology, Gaithersburg,
20899, MD, USA. Aber] Department of Computer Science, Aberystwyth University,
SY23 3DB, UK [Rutherford Appleton Laboratory] Scientific Computing Department,
Rutherford Appleton Laboratory, Science and Technology Facilities Council,
Harwell Campus, Didcot, OX11 0QX, UK Rutherford Appleton Laboratory]
Scientific Computing Department, Rutherford Appleton Laboratory, Science and
Technology Facilities Council, Harwell Campus, Didcot, OX11 0QX, UK Rutherford
Appleton Laboratory] Scientific Computing Department, Rutherford Appleton
Laboratory, Science and Technology Facilities Council, Harwell Campus, Didcot,
OX11 0QX, UK National Institute of Standards and Technology] Materials
Measurement Laboratory, National Institute of Standards and Technology,
Gaithersburg, 20899, MD, USA. IR,NMR,UV

##  1 Introduction

Large scale atomistic simulation of multi-component systems is a difficult
task but they are highly valuable for industrial applications such as
designing alloys, designing electrical contacts, touch screens, transistors,
batteries, composites and catalysts 1, 2, 3, 4. Quantum chemistry methods,
such as density functional theory, are obvious approaches to simulate such
systems, however, they are computationally very expensive for large systems 5.
Classical force-fields, or interatomic potentials, such as embedded-atom
method (EAM), modified embedded-atom method (MEAM), reactive bond-order
(ReaxFF), charge-optimized many-body (COMB) etc., 6, 7, 8, 9, 10 can be used
for such simulations but they are usually parameterized for a very narrow
chemical phase-space limiting their applicability and transferability.
Moreover, it can be quite strenuous and time consuming to develop such
traditional classical FFs.

Recently machine learning based FFs 11, 12 are being used to systematically
improve accuracy of FFs and have successfully been used for multiple systems.
One of the pioneer MLFFs were developed by Behler-Parinello in 2007 using
neural network 13. It was initially used for molecular systems and now has
been extended to numerous other applications 14. Although neural network is
one of the most popular regressors, other methods such as Gaussian processes-
based Gaussian approximation potentials (GAP) 15, as well as linear regression
and basis functions-based spectral neighbor analysis potential (SNAP) 16 have
also been thoroughly used. Such FFs use two and three body descriptors to
describe local environment. Other popular MLFF formalisms include smooth
overlap of atomic positions (SOAP) 17, moment tensor potentials (MTP) 18, 19,
symbolic regression 20 and polynomial-based approaches 21. One of the critical
issues in developing and maintaining classical force-fields is that they are
hard to update with software and hardware changes. Luckily, MLFFs are more
transparently developed and maintained than other classical FFs. A review
article on this topic can be found elsewhere 14. Nevertheless, early-
generation MLFFS are also limited to a narrow chemical space and may require
hand-crafted descriptors which may take time to be identified.

Graph neural network (GNN) based methods have shown remarkable improvements
over descriptor based machine learning methods and can capture highly non-
Euclidean chemical space 22, 23, 24, 25, 26, 12, 27, 28, 29. GNN FFs have been
also recently proposed and are still in the development phase 30, 31. We
developed an atomistic line graph neural network (ALIGNN) in our previous
work32 which can capture many body interactions in graph and successfully
models more than 70 properties of materials, either scalar or vector
quantities, such as formation energy, bandgap, elastic modulus,
superconducting properties, adsorption isotherm, electron and density of
states etc 32, 33, 34, 35, 36, 37. The same automatic differentiation
capability that allows training these complex models allows for physically
consistent prediction of quantities that must follow conservation laws, such
as forces and energies; this enables GNNs to be used in molecular dynamic
simulations to quickly identify relaxed or equilibrium states of complex
systems. However, there is a need for a large and diverse amount of data to
train unified force-fields.

In this work, we present a dataset of energy and forces with 4 million entries
for around 75000 materials in the JARVIS-DFT dataset which have been developed
over past 5 years 38. We extend the ALIGNN model to also predict derivatives
that are necessary for FF formalism. There can be numerous applications of
such a unified FF, however, in this work we limit ourselves to pre-
optimization of structures, genetic algorithm based structure, and molecular
dynamics applications. The developed model will be publicly available on the
ALIGNN GitHub page (https://github.com/usnistgov/alignn) with several examples
and a brief documentation.

##  2 Methodology

A flow-chart for training ALIGNN-FF is shown in Fig. 1a. To train the FF, we
use a large DFT dataset, JARVIS-DFT, which contains about 75000 materials with
a variety of atomic structures and chemistry and has been generated over the
last 5 years. JARVIS-DFT is primarily based on Vienna Ab initio Simulation
Package (VASP) 39, 40 software and OptB88vdW 41 functional but also contains
data obtained using other functionals and methods. In this work, only
OptB88vdW-based data has been used. The OptB88vdW functional was shown to be
very well applicable to solids in Ref. 41 and, ever since, it has been used to
model rare-gas dimers and metallic, ionic, and covalent bonded solids,
polymers, and small molecular systems 42. Energies and forces are available
for each structure optimization and elastic constant calculation runs. The
total number of such entries is around 4 million. Although it would be
justified to train on the entire dataset, we choose to use only a subset of it
because of the computational budget and hardware requirements available to us.
Instead of the 4 million datapoints, we use 307113 points i.e. more than an
order magnitude less, by taking unique set of first, last, middle, maximum
energy and minimum energy structures only. If some snapshots for a run are
identical, say the last step and the minimum energy configuration are the
same, then we only include one of them. The dataset consists of perfect
structures only.

We convert the atomic structures to a graph representation using atomistic
line graph neural network (ALIGNN). Details on ALIGNN can be found in the
related paper 32. In brief, each node in the atomistic graph is assigned 9
input node features based on its atomic species: electronegativity, group
number, covalent radius, valence electrons, first ionization energy, electron
affinity, block and atomic volume. The inter-atomic bond distances are used as
edge features with radial basis function up to 8 Å cut-off. We use a periodic
12-nearest-neighbor graph construction. This atomistic graph is then used for
constructing the corresponding line graph using interatomic bond-distances as
nodes and bond-angles as edge features. ALIGNN uses edge-gated graph
convolution for updating nodes as well as edge features. One ALIGNN layer
composes an edge-gated graph convolution on the bond graph with an edge-gated
graph convolution on the line graph. The line graph convolution produces bond
messages that are propagated to the atomistic graph, which further updates the
bond features in combination with atom features.

In this work, we developed the functionality for atomwise and gradient
predictions. Quantities related to gradients of the predicted energy, such as
forces on each atom, are computed by applying the chain rule through the
automatic differentiation system used to train the GNN, thus preserving the
law of conservation of energy. In a closed system, the forces on each atom $i$
depend on its position with respect to every other particle $j$ through a
force-field as:

$m_{i}\frac{d^{2}r_{i}(t)}{dt^{2}}=\sum_{j}F_{ij}(t)=-\sum_{j}\nabla_{i}U(r_{ij}(t))$
(1)

where $r_{ij}$ is the distance between atom $i$ and $j$. For stress
predictions, we use the Virial stress 43 which is a measure of mechanical
stress on an atomic scale for homogeneous systems:

$\sigma_{\alpha\beta}=-\frac{1}{V}\sum_{i}\sum_{i\neq
j}F_{ij}^{\alpha}r_{ij}^{\beta}+m_{i}v_{i}^{\alpha}v_{i}^{\beta}$ (2)

The 307113 data points are split into 90:5:5 ratio for training, validation
and testing. We train the model for 250 epochs using the same hyper-parameters
as in the original ALIGNN model 32. ALIGNN is based on deep graph library
(DGL) 44, PyTorch 45 and JARVIS-Tools packages 38. We optimize a composite
loss function ($l$) with weighted mean absolute error terms for both forces
and energies:

$l=\left|E^{DFT}-E^{GNN}\right|+w\sum_{i}^{N_{atoms}}|F_{i}^{DFT}-F_{i}^{GNN}|$
(3)

where, $E^{DFT}$ and $E^{GNN}$ are energies per atom using DFT and ALIGNN,
$F_{i}^{DFT}$ and $F_{i}^{GNN}$ are forces acting on an atom using DFT and
ALIGNN, $w$ is a tunable weighting factor scaling the force contribution to
the loss relative to the energies; in this work we set $w=10$. Without such
weighting, it is difficult to learn forces as they vary within a wide range.

The ALIGNN-FF model has been integrated with atomic simulation environment
(ASE) 46 as an energy, force and stress calculator for structure optimization
and MD simulation. This calculator can be used for optimizing atomic
structures, using genetic algorithm 47, and running molecular dynamics
simulations, for example constant-temperature, constant-volume ensemble (NVT)
simulations. The structural relaxations are carried out with the fast inertial
relaxation engine (FIRE) 48, available in ASE. In order to predict equation of
state/ energy-volume-curve (EV) simulation, we apply volumetric strains
between ranges of -0.1 to 0.1 with an interval of 0.01.

The interfaces are generated using Zur algorithm 49, 50, 51, 52 as implemented
in the JARVIS-Tools. For a few cases, we compare the ALIGNN-FF data with
classical-forcefield embedded atom method (EAM) 53, 54, 55 and density
functional theory based GPAW packages 56.

##  3 Results and discussion

###  3.1 Performance on test set

Figure 1: Schematic showing flow chart for developing the force-field,
visualization of available data and model performance on test dataset. The
left panel shows how dataset was obtained from JARVIS-DFT and model was
trained. Fig. a and b shows energy and force-distribution in the dataset.
While energy dataset ranges from -9.98 to 9.97 eV and is well-dispersed over
this range, the force dataset varies from -591.12 to 591.40 eV/Å and is highly
localized around zero. The force-dataset contains the x,y,z force on each
atom. Fig. d and Fig. e shows energy and force predictions on the test set.

We find that the energy per atom dataset (with 307113 entries) varies from
-9.98 eV/atom to 9.97 eV/atom. The force dataset contains 9593385 entries
(resulting from x,y,z forces on each atom) is mostly centered around zero, but
ranges from -591.12 eV/Å to 591.40 eV/Å. The mean absolute deviation (MAD) for
energies and forces are 1.80 eV and 0.10 eV/Å respectively. The distributions
of the energy and force used in the training are shown in Fig. 1b and 1c,
respectively. The total dataset was split in 90:5:5 train-validation-test
sets. We show the performance on test set for energies and forces in Fig. 1d
and 1e respectively. We note that these points represent a variety of
chemistry and structures unlike usual MLFFs which are usually focused on a
specific chemistry. The mean absolute error (MAE) of the predicted energy per
atom and force component per atom are 0.086 eV and 0.047 eV/Å respectively. In
comparison to the previous ALIGNN model for energy only, which has an MAE of
0.03 eV/atom, 0.086 eV/atom might seem high, however, we note that the
previous energy model was trained on relaxed structures only while the current
model also captures several un-relaxed structures which can be much higher in
energy scale. The MAD:MAE ratios for energies and forces are 20.93 and 2.12
which are high.

While the above results are for force tunable weighting factor of 10, we also
train the models with other weigthing factors as shown in Table 1. We find
that as we increase the weighting factors, the MAEs increase for energies
increase but decrease for forces. As the MAD for forces is 0.1 eV/Å, we choose
to work with the model obtained for the lowest MAE for forces and to analyze
its applications in the rest of the paper. Nevertheless, we share model
parameters for other weighting factors for those interested in analyzing its
effect on property predictions.

Table 1: Effect of different weighting factors for energy and force predictions. Weight | MAE-Energies (eV/atom) | MAE-Forces (eV/$\AA$)
---|---|---
0.1 | 0.034 | 0.092
0.5 | 0.044 | 0.089
1.0 | 0.051 | 0.088
5.0 | 0.082 | 0.054
10.0 | 0.086 | 0.047

###  3.2 Energy-volume curves

The energy-volume (E-V) curves are crucial to understand the behavior of an
FF. We obtain EV-curves for a few test case materials by applying volumetric
strains between ranges of -0.1 to 0.1 with an interval of 0.01. Specifically,
we compare the energy volume curves for Ni3Al (JVASP-14971), Al2CoNi
(JVASP-108163), CrFeCoNi (4 atoms cell with spacegroup number 216 and lattice
parameter of 4.007 Å), NaCl (JVASP-23862), MgO (JVASP-116), and BaTiO3
(JVASP-8029) using EAM potential 54, 55, GPAW DFT, and ALIGNN-FF in Fig. 2.
The energy scales for these methods differ so we align them with respect to
the corresponding minimum energies, for comparison. We used less EV-curve
points for a structure in GPAW to save computational cost. We notice all the
EV-curves are parabolic in nature and smooth, indicating a smooth potential
energy surface.

We find that the EV curves from these methods coincide near the minimum for
all systems but for Al2CoNi and CrFeCoNi, the GPAW equilibrium volume is
slightly smaller than for EAM and ALIGNN-FF, suggesting that the lattice
constants for ALIGNN-FF and EAM might be overestimated compared to GPAW.
Nevertheless, EAM and ALIGNN-FF data agree well. Comparing EAM and ALIGN-FF,
it’s important to remember that GNNs have no fundamental limitation to number
of species they can model (i.e. high chemical diversity), and can in principle
even extrapolate to species not contained in the training set, which is
extremely powerful compared to conventional FFs like EAM.

Note that there are many other conventional FF repositories available (such as
Inter-atomic Potential Repository 57 and JARVIS-FF 58) which contain data for
a variety of systems. It is beyond the scope of the current work to compare
all of them with ALIGNN-FF, however, it would be an interesting effort for the
future work.

Figure 2: Energy-volume/expansion-contraction curves for a few example
systems: a) Ni3Al, b) Al2CoNi, c) CrFeCoNi, d) NaCl, e) MgO, f) BaTiO3 with
ALIGNN-FF, EAM force-fields and GPAW DFT methods. We apply volumetric strains
between ranges of -0.1 to 0.1 with an interval of 0.01. For GPAW, we use an
interval of 0.05.

While the above examples are for individual crystals, it is important to
distinguish different polymorphs of a composition system for materials
simulation (i.e. structural diversity). In Fig. 3, we analyze the energy-
volume (EV) curve of four systems and their polymorphs using ALIGNN-FF. We
choose such four example systems because representative of different stable
structures. In general, however, the EV-curve can be computed for any
arbitrary system and structure. In Fig. 3a, we show the EV-curve for 4 silicon
materials (JARVIS-IDs: JVASP-1002, JVASP-91933, JVASP-25369, JVASP-25368) with
diamond cubic correctly being the lowest in energy. Similarly, the EV-curve
for naturally prevalent SiO2 systems (JARVIS-IDs:JVASP-58349, JVASP-34674,
JVASP-34656, JVASP-58394), binary alloy Ni3Al (JARVIS-IDs:JVASP-14971,
JVASP-99749, JVASP-11979) and vdW bonded material MoS2 (JARVIS-
IDs:JVASP-28733, JVASP-28413, JVASP-58505) all have the correct structure
corresponding to the minimum energy. Therefore, while the MAE for our overall
energy model is high, such model is able to distinguish polymorphs of
compounds with meV level accuracy which is critical for atomistic
applications.

Figure 3: Energy-volume/expansion-contraction curves for a few example
systems: a) silicon, b) SiO2, c) Ni3Al, d) MoS2 polymorphs. We optimize the
structures and then apply volumetric strains between ranges of -0.05 to 0.05
with an interval of 0.01.

###  3.3 Lattice constants and formation energies

Figure 4: Comparison of DFT and FF data for a) lattice constant in
x-direction, b) lattice-constant in y direction, c) lattice constant in
z-direction, d) formation energies for stable binary solids in JARVIS-DFT.
Comparison of crystallography open database (COD) and FF lattice constants in
e) x-direction, f) y-direction, g) z-direction and h) formation energy
distribution.

In this section, we compare the lattice constants and formation energies for
elemental and multi-component solids from the JARVIS-DFT after optimizing them
using ALIGNN-FF as shown in Fig. 4a-4d. For ALIGNN-FF we choose solids with
less than 10 atoms in a cell from the JARVIS-DFT database. We optimize 23495
such materials, and present the results in Fig 4a-4d. The formation energies
require chemical potential of elemental systems. The energy per atom of
elemental solid systems with minimum energy are used as chemical potentials.
We optimize the lattices with the FIRE algorithm as implemented in ASE and
find reasonable agreement between the DFT and ALIGNN-FF lattice constants. The
MAE for a,b,c lattice constants are 0.11, 0.11 and 0.13 Å respectively.The
mean absolute error (MAE) for formation energies between ALIGNN-FF and JARVIS-
DFT is 0.08 eV/atom which is reasonable for pre-screening applications. Note,
the above validation is different from the performance measurement in Fig. 1
for the 5 % test dataset because we optimize the crystals rather than directly
predicting the formation energies on unrelaxed structures. Similarly, we apply
the ALIGNN-FF on the crystallography open database (COD). The COD database
contains more than 431778 atomic structures with different types of chemical
bondings and environments. We optimize 34615 structures in the COD database
with number of atoms in a cell less than 50 and show the results in Fig.
4e-4h. Here, we find the MAE for a, b, c lattice parameter as 0.20, 0.20 and
0.23 Å respectively. Most of the systems in COD have been derived from
experiments, and hence we see many of them have negative formation energies as
shown in Fig. 4h.

###  3.4 Genetic algorithm based structure search

Computational prediction of the ground-state structures of a chemical system
is a challenging task. Some of the common methods for this task include
genetic algorithm (GA), simulated annealing and basin or minima hopping 59. In
the following examples, we show using GA together with ALIGNN-FF to search for
crystal structures of Ni-Al and Cu-Al example systems. Genetic algorithms
mimic the biological evolution process to solve optimization problems. A GA
optimization consists of 1) inheritance, 2) mutation, 3) selection, and 4)
crossover operations to produce new structures and search for better survivors
from generation to generation based on ”survival of the fittest” idea which in
our case would be energetics criteria. While all GAs follow similar strategy,
the details of the individual operations can vary a lot from problem to
problem and can be critical to search efficiency. We start with face-centered
cubic (FCC) structures of the individual components (FCC Al (JVASP-816) and
FCC Ni (JVASP-943) for Ni-Al search and FCC Al and FCC Cu (JVASP-867) for Cu-
Al search) as the initial population. We perform relaxation of these systems
with ALIGNN-FF. We choose a population size of 10 individuals and create the
initial population by randomly selecting the elements. We use 40 generations
to evolve the system and store the entries which are also relaxed with ALIGNN-
FF.

After this example GA search for structures, we plot the convex hull diagram
of these systems in Fig. 5. We find that the GA predicts AB and A3B compounds
which are in fact observed experimentally in such binary alloys 47, 60.
Additionally, Ni3Al (spacegroup: Pm-3m) is known to be one of the best
performing super-alloys 47 which is reproduced in the above example. We also
found that the formation energy of this structure (-0.47 eV/atom) is similar
to Johannesson’s47 findings of -0.49 eV/atom.Although the above example is
carried out for binary systems,in principle the same methodology can be
applied for any other systems as well.

Figure 5: Convex hulls obtained from a genetic algorithm search for alloys
starting with elemental solids only. The x-axis values represent mole fraction
of corresponding element. We show the stable structures only for a) Ni-Al
system, b) Cu-Au system.

###  3.5 MD-NVT simulations

Figure 6: An example NVT simulation run for MoS2(001)-Al2O3(001) and
Cu(111)-TiN(111) interfaces. Fig. a and c show snapshots of the atomic
structure and Fig. b and d show the time-evolution of temperature for the
systems.

Next, we show example constant-temperature, constant-volume ensemble (NVT)
simulation runs for MoS2(001)-Al2O3(001) and Cu(111)-TiN(111) interface in
Fig.6 using NVT ensemble and Berendsen thermostat. Several ensembles such as
NPT, NVT and NVE with several structure optimizers are already available in
ALIGNN-FF using the ASE package. The ALIGNN-FF is currently implemented as an
ASE calculator which can be used for several tasks. In this particular
MoS2-Al2O3 example, we first generate the interfaces with JARVIS-Tools and
then convert them into ASE-atoms object which is then subjected to structure
optimizer and Berendsen thermostat at 300 K for 0.1 ps with 0.01 fs timestep.
We observe how the thermostat gradually increases the temperature of the
system to 300 K and then maintains the temperature, as expected. The structure
keeps its integrity, i.e., doesn’t collapse.

###  3.6 Timing study

Now, we compare the time for a single step potential energy calculation for
FCC Aluminum (JVASP-816) with varied supercell sizes with EAM, ALIGNN-FF and
GPAW DFT 56 methods. All of these calculations are performed on CPU on a
personal laptop; GPU performance scaling may differ, particularly for large
cell sizes. We start with a unit cell and start making supercells with
[2,2,2], [3,3,3]. [4,4,4] and [5,5,5] dimensions. For GPAW case we use 8x8x8
k-points and as we make supercells, we reduce the k-points inversely
proportional to the supercell size. In Fig. 7, we find that out of these
three, EAM is the fastest method. The ALIGNN-FF is an order of magnitude
slower than EAM, as expected. While EAM potentials are considerably faster, we
note that they are difficult to train for multi-component systems. The GPAW
single step calculations are almost 100 times slower than ALIGNN-FF. As we add
multiple electronic and ionic relaxation steps in DFT, the computational cost
drastically increases. Note that we chose a very generalized DFT set up of 330
eV plane wave cutoff and PBE 61 exchange correlation functional for GPAW
calculations. As we add more k-points and plane wave cut-off, the
computational cost will increase.

Figure 7: Timing comparison for FCC Al system for EAM, ALIGNN-FF and GPAW DFT
methods.

##  4 Conclusion

In summary, we have developed a unified atomistic line graph neural network
based force-field (ALIGNN-FF) that can model diverse set of materials with any
combination of 89 elements from the periodic table. Using several test cases,
we demonstrate the application of this method to determine the several
properties such as lattice constants, formation energies, EV-curve for
different materials including metallic, ionic and van der Waal bonded
materials.

Although, the above examples are based on a few test cases, ALIGNN-FF can in
principle be used for several applications such as investigating defective
systems, high-entropy alloy, metal-organic framework, catalyst, battery
designs etc. and its validity needs to be tested for other applications which
is beyond the scope of the present work. Also, as the methods for including
larger datasets improve (such as training on millions of data), and
integration of active learning and transfer learning strategies are achieved,
we believe we can train more accurate models. Moreover, such universal FF
model can be integrated with universal tight-binding model 62 so that for
classical and quantum properties can be predicted for large systems. In the
current state, such FFs can be very useful for structure optimization,
however, there is lot of room for improvement in terms of other physical
characteristics such as defects, magnetism, charges, electronic levels etc.
Nevertheless, we believe such multi component interface simulations will
greatly expedite high-throughput computational screening of advanced
industrial materials.

K.C. thanks National Institute of Standards and Technology for funding and
computational support. In addition, K.C thanks Ruth Pachter and Kiet Ngyuen
from the Air Force Research Laboratory for helpful discussions. This work was
also partially supported by Wave 1 of the UKRI Strategic Priorities Fund under
the EPSRC grant EP/T001569/1, particularly the ‘AI for Science’ theme within
that grant, by the Alan Turing Institute and UKRI AIMLAC CDT, grant no.
EP/S023992/1.

## References

  * Ogale 2006 Ogale, S. B. _Thin films and heterostructures for oxide electronics_ ; Springer Science & Business Media, 2006
  * Andersson et al. 2006 Andersson, M. P.; Bligaard, T.; Kustov, A.; Larsen, K. E.; Greeley, J.; Johannessen, T.; Christensen, C. H.; Nørskov, J. K. Toward computational screening in heterogeneous catalysis: Pareto-optimal methanation catalysts. _Journal of Catalysis_ 2006, _239_ , 501–506
  * Liang et al. 2013 Liang, T.; Shan, T.-R.; Cheng, Y.-T.; Devine, B. D.; Noordhoek, M.; Li, Y.; Lu, Z.; Phillpot, S. R.; Sinnott, S. B. Classical atomistic simulations of surfaces and heterogeneous interfaces with the charge-optimized many body (COMB) potentials. _Materials Science and Engineering: R: Reports_ 2013, _74_ , 255–279
  * Li et al. 2009 Li, X.; Zhu, Y.; Cai, W.; Borysiak, M.; Han, B.; Chen, D.; Piner, R. D.; Colombo, L.; Ruoff, R. S. Transfer of large-area graphene films for high-performance transparent conductive electrodes. _Nano letters_ 2009, _9_ , 4359–4363
  * Srolovitz and Vitek 2012 Srolovitz, D. J.; Vitek, V. _Atomistic Simulation of Materials: Beyond Pair Potentials_ ; Springer Science & Business Media, 2012
  * Daw and Baskes 1984 Daw, M. S.; Baskes, M. I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals. _Physical Review B_ 1984, _29_ , 6443
  * Liang et al. 2012 Liang, T.; Devine, B.; Phillpot, S. R.; Sinnott, S. B. Variable charge reactive potential for hydrocarbons to simulate organic-copper interactions. _The Journal of Physical Chemistry A_ 2012, _116_ , 7976–7991
  * Brenner 1990 Brenner, D. W. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films. _Physical review B_ 1990, _42_ , 9458
  * Van Duin et al. 2001 Van Duin, A. C.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: a reactive force field for hydrocarbons. _The Journal of Physical Chemistry A_ 2001, _105_ , 9396–9409
  * Brenner 2000 Brenner, D. W. The art and science of an analytic potential. _physica status solidi (b)_ 2000, _217_ , 23–40
  * Poltavsky and Tkatchenko 2021 Poltavsky, I.; Tkatchenko, A. Machine learning force fields: Recent advances and remaining challenges. _The Journal of Physical Chemistry Letters_ 2021, _12_ , 6551–6564
  * Choudhary et al. 2022 Choudhary, K.; DeCost, B.; Chen, C.; Jain, A.; Tavazza, F.; Cohn, R.; Park, C. W.; Choudhary, A.; Agrawal, A.; Billinge, S. J., et al. Recent advances and applications of deep learning methods in materials science. _npj Computational Materials_ 2022, _8_ , 1–26
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. _Physical review letters_ 2007, _98_ , 146401
  * Unke et al. 2021 Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müllerr, K.-R. Machine learning force fields. _Chemical Reviews_ 2021, _121_ , 10142–10186
  * Bartók et al. 2010 Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. _Physical review letters_ 2010, _104_ , 136403
  * Wood and Thompson 2018 Wood, M. A.; Thompson, A. P. Extending the accuracy of the SNAP interatomic potential form. _The Journal of chemical physics_ 2018, _148_ , 241721
  * Bartók et al. 2013 Bartók, A. P.; Kondor, R.; Csányi, G. On representing chemical environments. _Physical Review B_ 2013, _87_ , 184115
  * Shapeev 2016 Shapeev, A. V. Moment tensor potentials: A class of systematically improvable interatomic potentials. _Multiscale Modeling & Simulation_ 2016, _14_ , 1153–1173
  * Novikov et al. 2020 Novikov, I. S.; Gubaev, K.; Podryabinkin, E. V.; Shapeev, A. V. The MLIP package: moment tensor potentials with MPI and active learning. _Machine Learning: Science and Technology_ 2020, _2_ , 025002
  * Hernandez et al. 2019 Hernandez, A.; Balasubramanian, A.; Yuan, F.; Mason, S. A.; Mueller, T. Fast, accurate, and transferable many-body interatomic potentials by symbolic regression. _npj Computational Materials_ 2019, _5_ , 1–11
  * Drautz 2019 Drautz, R. Atomic cluster expansion for accurate and transferable interatomic potentials. _Physical Review B_ 2019, _99_ , 014104
  * Schütt et al. 2018 Schütt, K. T.; Sauceda, H. E.; Kindermans, P.-J.; Tkatchenko, A.; Müller, K.-R. Schnet–a deep learning architecture for molecules and materials. _The Journal of Chemical Physics_ 2018, _148_ , 241722
  * Xie and Grossman 2018 Xie, T.; Grossman, J. C. Crystal graph convolutional neural networks for an accurate and interpretable prediction of material properties. _Physical review letters_ 2018, _120_ , 145301
  * Chen et al. 2019 Chen, C.; Ye, W.; Zuo, Y.; Zheng, C.; Ong, S. P. Graph networks as a universal machine learning framework for molecules and crystals. _Chemistry of Materials_ 2019, _31_ , 3564–3572
  * Chen and Ong 2022 Chen, C.; Ong, S. P. A universal graph deep learning interatomic potential for the periodic table. _arXiv preprint arXiv:2202.02450_ 2022,
  * Kearnes et al. 2016 Kearnes, S.; McCloskey, K.; Berndl, M.; Pande, V.; Riley, P. Molecular graph convolutions: moving beyond fingerprints. _Journal of computer-aided molecular design_ 2016, _30_ , 595–608
  * Gilmer et al. 2017 Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; Dahl, G. E. Neural message passing for quantum chemistry. International conference on machine learning. 2017; pp 1263–1272
  * Klicpera et al. 2020 Klicpera, J.; Giri, S.; Margraf, J. T.; Günnemann, S. Fast and uncertainty-aware directional message passing for non-equilibrium molecules. _arXiv preprint arXiv:2011.14115_ 2020,
  * Batzner et al. 2022 Batzner, S.; Musaelian, A.; Sun, L.; Geiger, M.; Mailoa, J. P.; Kornbluth, M.; Molinari, N.; Smidt, T. E.; Kozinsky, B. E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials. _Nature communications_ 2022, _13_ , 1–11
  * Park et al. 2021 Park, C. W.; Kornbluth, M.; Vandermause, J.; Wolverton, C.; Kozinsky, B.; Mailoa, J. P. Accurate and scalable graph neural network force field and molecular dynamics with direct force architecture. _npj Computational Materials_ 2021, _7_ , 1–9
  * Chmiela et al. 2018 Chmiela, S.; Sauceda, H. E.; Müller, K.-R.; Tkatchenko, A. Towards exact molecular dynamics simulations with machine-learned force fields. _Nature communications_ 2018, _9_ , 1–10
  * Choudhary and DeCost 2021 Choudhary, K.; DeCost, B. Atomistic Line Graph Neural Network for improved materials property predictions. _npj Computational Materials_ 2021, _7_ , 1–8
  * Choudhary et al. 2022 Choudhary, K.; Yildirim, T.; Siderius, D. W.; Kusne, A. G.; McDannald, A.; Ortiz-Montalvo, D. L. Graph neural network predictions of metal organic framework CO2 adsorption properties. _Computational Materials Science_ 2022, _210_ , 111388
  * Choudhary and Garrity 2022 Choudhary, K.; Garrity, K. Designing High-Tc Superconductors with BCS-inspired Screening, Density Functional Theory and Deep-learning. _arXiv preprint arXiv:2205.00060_ 2022,
  * Choudhary and Sumpter 2022 Choudhary, K.; Sumpter, B. G. A Deep-learning Model for Fast Prediction of Vacancy Formation in Diverse Materials. _arXiv preprint arXiv:2205.08366_ 2022,
  * Kaundinya et al. 2022 Kaundinya, P. R.; Choudhary, K.; Kalidindi, S. R. Prediction of the electron density of states for crystalline compounds with Atomistic Line Graph Neural Networks (ALIGNN). _JOM_ 2022, _74_ , 1395–1405
  * Gurunathan et al. 2022 Gurunathan, R.; Choudhary, K.; Tavazza, F. Rapid Prediction of Phonon Structure and Properties using the Atomistic Line Graph Neural Network (ALIGNN). _arXiv preprint arXiv:2207.12510_ 2022,
  * Choudhary et al. 2020 Choudhary, K.; Garrity, K. F.; Reid, A. C.; DeCost, B.; Biacchi, A. J.; Hight Walker, A. R.; Trautt, Z.; Hattrick-Simpers, J.; Kusne, A. G.; Centrone, A., et al. The joint automated repository for various integrated simulations (JARVIS) for data-driven materials design. _npj Computational Materials_ 2020, _6_ , 1–13
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Physical review B_ 1996, _54_ , 11169
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Computational materials science_ 1996, _6_ , 15–50
  * Klimeš et al. 2009 Klimeš, J.; Bowler, D. R.; Michaelides, A. Chemical accuracy for the van der Waals density functional. _Journal of Physics: Condensed Matter_ 2009, _22_ , 022201
  * Choudhary et al. 2018 Choudhary, K.; Cheon, G.; Reed, E.; Tavazza, F. Elastic properties of bulk and low-dimensional materials using van der Waals density functional. _Physical Review B_ 2018, _98_ , 014107
  * Subramaniyan and Sun 2008 Subramaniyan, A. K.; Sun, C. Continuum interpretation of virial stress in molecular simulations. _International Journal of Solids and Structures_ 2008, _45_ , 4340–4346
  * Wang et al. 2019 Wang, M.; Zheng, D.; Ye, Z.; Gan, Q.; Li, M.; Song, X.; Zhou, J.; Ma, C.; Yu, L.; Gai, Y., et al. Deep graph library: A graph-centric, highly-performant package for graph neural networks. _arXiv preprint arXiv:1909.01315_ 2019,
  * Paszke et al. 2017 Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.; DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; Lerer, A. Automatic differentiation in pytorch. 2017,
  * Larsen et al. 2017 Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C., et al. The atomic simulation environment—a Python library for working with atoms. _Journal of Physics: Condensed Matter_ 2017, _29_ , 273002
  * Johannesson et al. 2002 Johannesson, G. H.; Bligaard, T.; Ruban, A. V.; Skriver, H. L.; Jacobsen, K. W.; Nørskov, J. K. Combined electronic structure and evolutionary search approach to materials design. _Physical Review Letters_ 2002, _88_ , 255506
  * Bitzek et al. 2006 Bitzek, E.; Koskinen, P.; Gähler, F.; Moseler, M.; Gumbsch, P. Structural relaxation made simple. _Physical review letters_ 2006, _97_ , 170201
  * Zur and McGill 1984 Zur, A.; McGill, T. Lattice match: An application to heteroepitaxy. _Journal of applied physics_ 1984, _55_ , 378–386
  * Mathew et al. 2016 Mathew, K.; Singh, A. K.; Gabriel, J. J.; Choudhary, K.; Sinnott, S. B.; Davydov, A. V.; Tavazza, F.; Hennig, R. G. MPInterfaces: A Materials Project based Python tool for high-throughput computational screening of interfacial systems. _Computational Materials Science_ 2016, _122_ , 183–190
  * Ding et al. 2016 Ding, H.; Dwaraknath, S. S.; Garten, L.; Ndione, P.; Ginley, D.; Persson, K. A. Computational approach for epitaxial polymorph stabilization through substrate selection. _ACS applied materials & interfaces_ 2016, _8_ , 13086–13093
  * Choudhary et al. 2020 Choudhary, K.; Garrity, K. F.; Pilania, G.; Tavazza, F. Efficient computational design of 2D van der Waals Heterostructures: band-alignment, lattice-mismatch, web-app generation and machine-learning. _arXiv preprint arXiv:2004.03025_ 2020,
  * Daw and Baskes 1983 Daw, M. S.; Baskes, M. I. Semiempirical, quantum mechanical calculation of hydrogen embrittlement in metals. _Physical review letters_ 1983, _50_ , 1285
  * Pun et al. 2015 Pun, G. P.; Yamakov, V.; Mishin, Y. Interatomic potential for the ternary Ni–Al–Co system and application to atomistic modeling of the B2–L10 martensitic transformation. _Modelling and Simulation in Materials Science and Engineering_ 2015, _23_ , 065006
  * Farkas and Caro 2020 Farkas, D.; Caro, A. Model interatomic potentials for Fe–Ni–Cr–Co–Al high-entropy alloys. _Journal of Materials Research_ 2020, _35_ , 3031–3040
  * Enkovaara et al. 2010 Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dułak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H., et al. Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. _Journal of physics: Condensed matter_ 2010, _22_ , 253202
  * Becker et al. 2013 Becker, C. A.; Tavazza, F.; Trautt, Z. T.; de Macedo, R. A. B. Considerations for choosing and using force fields and interatomic potentials in materials science and engineering. _Current Opinion in Solid State and Materials Science_ 2013, _17_ , 277–283
  * Choudhary et al. 2017 Choudhary, K.; Congo, F. Y. P.; Liang, T.; Becker, C.; Hennig, R. G.; Tavazza, F. Evaluation and comparison of classical interatomic potentials through a user-friendly interactive web-interface. _Scientific data_ 2017, _4_ , 1–12
  * Ji et al. 2010 Ji, M.; Wang, C.-Z.; Ho, K.-M. Comparing efficiencies of genetic and minima hopping algorithms for crystal structure prediction. _Physical Chemistry Chemical Physics_ 2010, _12_ , 11617–11623
  * van de Walle and Ceder 2002 van de Walle, A.; Ceder, G. Automating first-principles phase diagram calculations. _Journal of Phase Equilibria_ 2002, _23_ , 348–359
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Physical review letters_ 1996, _77_ , 3865
  * Garrity and Choudhary 2021 Garrity, K. F.; Choudhary, K. Fast and Accurate Prediction of Material Properties with Three-Body Tight-Binding Model for the Periodic Table. _arXiv preprint arXiv:2112.11585_ 2021,


# Ai4EComponentLib.jl: A Component-base Model Library in Julia

Yuebao Yang1, Jingyi Yang2, Mingtao Li111Corresponding author
Xi’an Jiaotong University

###### Abstract

Ai4EComponentLib.jl222https://ai4energy.github.io/Ai4EComponentLib.jl/dev/(Ai4EComponentLib)
is a component-base model library based on Julia language, which relies on the
differential equation solver DifferentialEquations.jl[1] and the symbolic
modeling tool Modelingtoolkit.jl[2]. To handle problems in different physical
domains, Ai4EComponentLib tries to build them with component-base model.
Supported by a new generation of symbolic modeling tools, models built with
Ai4EComponentLib are more flexible and scalable than models built with
traditional tools like Modelica. This paper will introduce the instance and
general modeling methods of Ai4EComponentLib model library.

##  1 Motivation

In the field of modeling and simulation, there are already many modeling
tools, such as Modelica[3], gProms[4], Simulink[5] and so on. They play an
important role in different fields, and these tools have become essential
tools for scientific researchers. But that doesn’t mean they are perfect.
Taking the open source Modelica as an example, Modelica has its own compiler
that can compile the modeling language into C language, which is a complete
solution chain for simulation problems. Its disadvantage is that it is not
easy to expand and it cannot generate solutions of optimization problems
(parameter identification, optimal control problems), etc. The ModelingToolkit
in the Julia ecosystem solves the problem of insufficient scalability. Based
on the work of ModelingToolkit, Ai4EComponentLib has built some model
libraries, hoping to take advantage of its high scalability in the future.

##  2 Ai4EComponentLib System Design

The method of Ai4EComponentLib design system is the same with the traditional
component modeling idea. But due to the high scalability of ModelingToolkit,
we will look at this problem with a higher level of abstraction. The core of
component-base model is the process, and the changes in energy and matter in
the process. Two key points when designing the system are: internal process
and external connection. To design a system is to design their processes and
connections.

In a system, there are 3 components A, B and C in the system. As shown in the
Figure1, each component has its own input and output nodes (such as a1, a2,
etc.), and its inputs and outputs are shown by arrows.

Figure 1: System with 3 components A, B and C.

###  2.1 Internal Process

Each component has its own characteristics, which is reflected in the process
of the reaction between the input and output in the node or the exchange of
material and energy. And it varies from component to component. These
components are the mapping of real world to models and are abstractions of
different physical phenomena.

For example, a1, a2, a3 of Component A in the above Figure1. Inside the
components, they have a specific mathematical relationship.

$f(a1,a2,a3)=0$

Here, $f(a1,a2,a3)$ is a general functional form, which can be a differential
equation or an algebraic equation. Take a resistor component as an example,
the same voltage difference is applied across a different resistor, and the
current flowing through it is different. If the resistance changes with time,
then the current of the system also changes with time. "Different resistance"
and "time-varying resistance" are both characteristics of components, which
are described by different equations $f$ when designing components. These
equations determine what happens to "matter and energy" as it flows through
the component. Once the components in the system are designed, the function of
the whole system can be determined by these components.

###  2.2 External connection

When building system connections, we focus on the matter and energy of the
connection points.

  * •

If it is a circuit system. The a2 in Figure1 of the circuit system is used as
the outlet of the component, and the current and voltage values are necessary
attributes. And currents and voltage exist at the inlet and outlet of each
component. So current and voltage are the connetor of the circuit system.

  * •

If it is a pipe system(without considering potential energy). The a2 in
Figure1 of the pipe system is used as the outlet of the component, and the
pressure and velocity of the water flow are necessary attributes. Water
pressure and velocity properties exist at the inlet and outlet of each
component. So pressure and velocity are the connetor of the pipe system.

The specific components needs to be build according to the physical model of
the system. For example, Kirchhoff’s laws in circuits and Bernoulli’s equation
for pipe system. The governing equations are the basic theory for designing
the system.

In the external connection, the conservation law must generally be followed -
the conservation of energy and the conservation of mass:

$\left\\{\begin{array}[]{l}|a2|=|b1|\\\ |b3|=|c2|\\\
|a3|=|c1|\end{array}\right.$

For state variables, such as voltage, water pressure, gas pressure. Generally
there are:

$\left\\{\begin{array}[]{l}a2=b1\\\ b3=c2\\\ a3=c1\end{array}\right.$

For process variables such as current, water flow, gas flow. Generally there
are (specify that inflow is positive and outflow is negative):

$\left\\{\begin{array}[]{l}a2+b1=0\\\ b3+c2=0\\\ a3+c1=0\end{array}\right.$

##  3 Case Study

###  3.1 Incompressible System

When the liquid in the pipe is an incompressible fluid, the governing equation
of the system is Bernoulli’s equation[6].

$\frac{p}{\rho g}+\frac{v^{2}}{2g}+h=\text{ constant }$

Where, $p$ is pressure,$v$ is velocity, $h$ is height. Bernoulli’s equation
describes the law of conservation of energy for the flow of liquid inside a
pipe. The internal process of the component is the increase or decrease of
energy.

####  3.1.1 Pipe Components

For a straight pipeline, there are generally resistance losses along the route
and local resistance losses inside the component. As the liquid flows from one
port of the pipe to the other, the total energy of the liquid decreases. Along
the way resistance loss and local resistance loss are internal processes of
the pipeline assembly. Its equation is:

$\frac{p_{in}}{\rho
g}+\frac{8q_{in}^{2}}{2\pi^{2}D^{4}g}+z_{in}=\frac{p_{out}}{\rho
g}+\frac{8q_{out}^{2}}{2\pi^{2}D^{4}g}+z_{out}+h_{f}+h_{m}$

Where the subscript in represents the pipeline inlet, out represents the
pipeline outlet, ${h_{f}}$ represents the resistance loss along the way, and
${h_{m}}$ represents the local resistance loss.

$h_{f}=f\frac{L}{D}\frac{8q^{2}}{\pi^{2}D^{4}g}$

Where $f$ is the friction factor, $L$ is the pipe length, and $D$ is the pipe
diameter.

$h_{m}=K\frac{8q^{2}}{\pi^{2}D^{4}g}$

Where $K$ is the local resistance coefficient.

The external connection of the pipeline is mainly pressure and flow (the
height is not considered). Pressure is a state variable, and flow is a process
variable. Between the two components a, b, there are:

$p_{a}=p_{b}\\\ q_{a}+q_{b}=0$

####  3.1.2 Centrifugal Pump Components

The internal process of a centrifugal pump is to add energy to the fluid.

$\frac{p_{in}}{\rho
g}+\frac{8q_{in}^{2}}{2\pi^{2}D^{4}g}+z_{in}=\frac{p_{out}}{\rho
g}+\frac{8q_{out}^{2}}{2\pi^{2}D^{4}g}+z_{out}+H_{t}$

where $H_{t}$ represents the added energy when passing through the centrifugal
pump.

$H_{t}=\frac{(r\omega)^{2}}{g}-\frac{\omega\cot\beta}{2\pi
bg}Q=c_{0}\omega^{2}-c_{1}\omega Q=a_{0}-a_{1}Q$

Where $a_{0},a_{1}$ represent the characteristic parameters of the pump, and
Equation 10 is actually the theoretical head- flow curve of the pump.

####  3.1.3 Components and Systems

Before building the main components of the system such as pumps and pipes, we
should build the Connector in ModelingToolkit (the variables in the Connector
are pressure and flow in the pipeline system). Detailed code can be viewed in
the code repository.

⬇

# Component: SimplePipe(pipe with fixed friction factor `f`)

function SimplePipe(; name, L=10.0, D=25E-3, f=0.01, $\rho$=1.0E3, zin=0.0,
zout=0.0, K_inside=0.0)

@named in = PipeNode(z=zin)

@named out = PipeNode(z=zout)

ps = @parameters D = D L = L f = f K_inside = K_inside

eqs = [

_NodeEnergy(in, D, $\rho$) ~ _NodeEnergy(out, D, $\rho$) + _h_f(in, f, L, D) +
_h_m(in, K_inside, D)

0 ~ in.q + out.q

]

compose(ODESystem(eqs, t, [], ps, name=name), in, out)

end

# Component: CentrifugalPump

function CentrifugalPump(; name, D=25E-3, $\omega$=2500, c_0=4.4e-4,
c_1=5.622, $\rho$=1.0E3)

@named in = PipeNode()

@named out = PipeNode()

a_0 = c_0 * abs2($\omega$ * 2$\pi$ / 60)

a_1 = c_1 * $\omega$ * 2$\pi$ / 60

ps = @parameters D = D

eqs = [

_NodeEnergy(in, D, $\rho$) + a_0 - a_1 * abs(in.q) ~ _NodeEnergy(out, D,
$\rho$)

0 ~ in.q + out.q

]

compose(ODESystem(eqs, t, [], ps, name=name), in, out)

end

# Component: Sink_P

function Sink_P(; name, p=101325)

@named port = PipeNode(z=0.0)

eqs = [

port.p ~ p

]

compose(ODESystem(eqs, t, [], [], name=name), port)

end

It becomes simple and fast to build the following system from the above
components, as shown in Figure2.

Figure 2: An Example of Pipe System

⬇

@named A = Sink_P()

@named B = Sink_P()

system = [A, B, Pump]

@named Pipe1 = SimplePipe(L=2.0);

push!(system, Pipe1);

@named Pipe2 = SimplePipe(L=3.0);

push!(system, Pipe2);

@named Pipe3 = SimplePipe(L=7.0);

push!(system, Pipe3);

@named Pipe4 = SimplePipe(L=9.0);

push!(system, Pipe4);

@named Pipe5 = SimplePipe(L=5.0);

push!(system, Pipe5);

@named Pipe6 = SimplePipe(L=4.0);

push!(system, Pipe6);

@named Pipe7 = SimplePipe(L=5.0);

push!(system, Pipe7);

@named Pipe8 = SimplePipe(L=1.0);

push!(system, Pipe8);

@named Pipe9 = SimplePipe(L=10.0);

push!(system, Pipe9);

@named Pipe10 = SimplePipe(L=2.0);

push!(system, Pipe10);

@named Pipe11 = SimplePipe(L=2.0);

push!(system, Pipe11);

@named Pipe12 = SimplePipe(L=3.0);

push!(system, Pipe12);

@named Pipe13 = SimplePipe(L=12.0);

push!(system, Pipe13);

@named Pipe14 = SimplePipe(L=1.0);

push!(system, Pipe14);

@named Pipe15 = SimplePipe(L=2.0);

push!(system, Pipe15);

@named Pipe16 = SimplePipe(L=3.0);

push!(system, Pipe16);

@named Pipe17 = SimplePipe(L=6.0);

push!(system, Pipe17);

@named Pipe18 = SimplePipe(L=6.0);

push!(system, Pipe18);

@named Pipe19 = SimplePipe(L=6.0);

push!(system, Pipe19);

@named Pipe20 = SimplePipe(L=1.0);

push!(system, Pipe20);

@named Pipe21 = SimplePipe(L=1.0);

push!(system, Pipe21);

@named Pipe22 = SimplePipe(L=7.0);

push!(system, Pipe22);

@named Pipe23 = SimplePipe(L=3.0);

push!(system, Pipe23);

@named Pipe24 = SimplePipe(L=3.0);

push!(system, Pipe24);

@named Pipe25 = SimplePipe(L=2.0);

push!(system, Pipe25);

eqs = [

connect(A.port, Pump.in)

connect(Pump.out, Pipe1.in)

connect(Pipe1.out, Pipe2.in, Pipe5.in)

connect(Pipe2.out, Pipe3.in, Pipe6.in)

connect(Pipe3.out, Pipe4.in, Pipe7.in)

connect(Pipe4.out, Pipe10.out, Pipe14.in)

connect(Pipe5.out, Pipe11.in, Pipe12.in)

connect(Pipe6.out, Pipe8.in, Pipe9.in)

connect(Pipe7.out, Pipe9.out, Pipe10.in)

connect(Pipe12.out, Pipe8.out, Pipe13.in)

connect(Pipe13.out, Pipe14.out, Pipe15.in)

connect(Pipe11.out, Pipe19.in, Pipe16.in)

connect(Pipe16.out, Pipe17.in, Pipe20.in)

connect(Pipe17.out, Pipe18.in, Pipe21.in)

connect(Pipe18.out, Pipe15.out, Pipe22.in)

connect(Pipe19.out, Pipe20.out, Pipe23.in)

connect(Pipe21.out, Pipe22.out, Pipe24.in)

connect(Pipe23.out, Pipe24.out, Pipe25.in)

connect(B.port, Pipe25.out)

]

###  3.2 Thermodynamic System

####  3.2.1 Math work

In a thermodynamic cycle system, there are several typical internal
processes[7].

  * •

Isothermal Process: $\frac{T}{P}=constant$

  * •

Isobaric Process: $\frac{T}{v}=constant$

  * •

Isentropic Process: $pv^{k}=constant$

  * •

Isometric Process: $pv=constant$

The energy changes in these processes are ultimately reflected in the
pressure, density, enthalpy, entropy, and temperature of the fluid. In the
thermodynamic cycle, the external connection contains these 5 state variables.

When build a thermodynamic system, although governing equations for different
internal processes can be used, there will always be some discrepancies
between the theoretical equations and the real world. Therefore, when building
a thermal system, we can take advantage of the scalability feature of
ModelingToolkit—calling an external property library. Calling the property
library can reduce the error of some theoretical calculations.

Taking an isothermal process as an example, the state of the starting point of
the process is known. The governing equation of its internal process is:

$T_{1}=T_{2}$

Then determine another state variable at the end of the process, the other
states at the end can be obtained from two known states. Then there are the
following 4 combinations:

$\begin{matrix}T_{2},p_{2}\Rightarrow s_{2},D_{2},h_{2}\\\
T_{2},s_{2}\Rightarrow p_{2},D_{2},h_{2}\\\ T_{2},D_{2}\Rightarrow
s_{2},p_{2},h_{2}\\\ T_{2},h_{2}\Rightarrow s_{2},p_{2},D_{2}\end{matrix}$

Using the CoolProp property library, enter the two parameters on the left side
of the equation, and we can get any parameter on the right side of the
equation. Therefore, the internal equations can be replaced by the external
CoolProp library.

####  3.2.2 Components and Systems

Taking the isenthalpic process as an example, the internal equation is:

$h_{1}=h_{2}$

If a parameter $p$ is given, the function choose_equations will automatically
generate code that calls CoolProp through $p,h$ to get the values of other
parameters.

⬇

function IsoenthalpyProcess(; name, inter_state="Q_0", fluid="Water")

@assert inter_state != "H" "IsoenthalpyProcess can't accept H. Please chose
another state."

@named oneport = StreamPort()

@unpack $\Delta$h, out = oneport

eqs = [

$\Delta$h ~ 0

]

push!(eqs, chose_equations(out, inter_state, "H", fluid)...)

return extend(ODESystem(eqs, t, [], []; name=name), oneport)

end

Consider the reheat cycle depicted in Figure3:

Figure 3: Reheat Rankine Cycle

The code is as follows. For more detailed, please check the code repository
and documentation.

⬇

@named pump = IsentropicProcess(inter_state="P")

@named pump_P = DThermalStates(state="P", value=-1.0E5, u0=18.0E6)

@named boiler = IsobaricProcess(inter_state="T")

@named boiler_T = ThermalStates(state="T", value=550+ 273.15)

@named turbine = IsentropicProcess(inter_state="P")

@named turbine_P = ThermalStates(state="P", value=3.0e6)

@named reboiler = IsobaricProcess(inter_state="T")

@named reboiler_T = ThermalStates(state="T", value=450 + 273.15)

@named returbine = IsentropicProcess(inter_state="P")

@named returbine_P = ThermalStates(state="P", value=4.0e3)

@named condenser = IsothermalProcess(inter_state="Q_0")

eqs = [

connect(pump.out, boiler.in, pump_P.node)

connect(boiler.out, turbine.in, boiler_T.node)

connect(turbine.out, reboiler.in, turbine_P.node)

connect(reboiler.out, returbine.in, reboiler_T.node)

connect(returbine.out, condenser.in, returbine_P.node)

connect(condenser.out, pump.in)

]

##  4 Advantages of Component-base Models

To summarize, components are essentially equations that describe the real
world. The governing equations inside the component are the most important
part, which lead the change of the system. External connections describe the
relationship between components by connection equations. The connection
equation does not affect the real system evolution, but only assists building
the system.

Components-base model are designed to disassemble large systems into small
components. The advantage is that it is convenient for users to build the
system, and the disadvantage is that these connection variables will consume
computing resources to be simplified after the system is built. In fact, the
benefits of a componentized model far outweigh the cost of consuming a little
computing resources.

For example, if the system composed of components A,B,C are regarded as a
whole X. Then A,B,C become the internal structure of X, and its mathematical
form is also a differential algebraic equation, which is not different from
the mathematical form of components A,B,C. The only difference is the number
of equations. System X can be treated as a component which has three external
interfaces: input interface x1, x2 and output interface x3, as the Figure4
shows.

Figure 4: System X composed of components A,B,C

Now, the system and subsystem have a complete concept of closure. If A,B,C are
not a minimum system at this time, it can be further subdivided into smaller
components. A,B,C can be called a subsystem of system X, and system A,B,C also
has its own subsystem.

If A,B,C have only a few equations in total, then A,B,C are not required to be
constructed as a subsystem of X. We directly flatten the equation inside X
manually. But if A,B,C have thousands of equations, and they also have their
own subsystems. Is flattening them manually still a good way?

Therefore, the problems that are divided into components is of great
significance when build a giant system. Once the most critical internal
structure are determined, people can concentrate on implementing the system
(subsystem) at a high level of abstraction.

The rest of the simplification work is left to the computer! If people
constantly flatten the equation and build the system from scratch, wouldn’t it
be a waste of human resources? The value of human thinking should be reflected
in the design of the system, not in the flattening of the equation.

##  5 Conclusion

The component-base model has its own advantages and has a broader application
prospect. The callback system of ModelingToolkit and the solver of
DifferentialEquations can handle simulation problems more flexibly.
OptControl.jl[8] based on ModelingToolkit can transform the component model
system into an optimal control problem solved by JuMP. In addition, there is
SciML’s333https://sciml.ai/ parameter identification toolkit and so on. With
the support of the Julia community, a component-base model is not only
equations for simulation, but it can also be a state space model in control
problems, optimizing constraints in problems. Because the essence of
components - equations can be described by the most basic Julia element -
functions.

## References

  * [1] Christopher Rackauckas and Qing Nie.  Differentialequations.jl--a performant and feature-rich ecosystem for solving differential equations in julia.  Journal of Open Research Software, 5(1), 2017.
  * [2] Yingbo Ma, Shashi Gowda, Ranjan Anantharaman, Chris Laughman, Viral Shah, and Chris Rackauckas.  Modelingtoolkit: A composable graph transformation system for equation-based modeling.  arXiv preprint arXiv:2103.05244, 2021.
  * [3] Michael Tiller.  Introduction to physical modeling with Modelica.  Springer Science & Business Media, 2001.
  * [4] Mariano Asteasuain, Stella Maris Tonelli, Adriana Brandolin, and Jose Alberto Bandoni.  Dynamic simulation and optimisation of tubular polymerisation reactors in gproms.  Computers & Chemical Engineering, 25(4-6):509--515, 2001.
  * [5] Devendra K Chaturvedi.  Modeling and simulation of systems using MATLAB® and Simulink®.  CRC press, 2017.
  * [6] Bruce Roy Munson, Theodore Hisao Okiishi, Wade W Huebsch, and Alric P Rothmayer.  Fluid mechanics.  Wiley Singapore, 2013.
  * [7] Michael J Moran and George Tsatsaronis.  Engineering thermodynamics.  In CRC Handbook of Thermal Engineering Second Edition, pages 1--112. CRC Press, 2017.
  * [8] Jingyi Yang, Yuebao Yang, and Mingtao Li.  Optcontrol. jl: An interpreter for optimal control problem.  arXiv preprint arXiv:2207.13229, 2022.


# Capacity Cost of Fulfilling the URLLC Performance in Industrial 5G New Radio
Deployments

Ali A. Esswie
Cellullar Standards Technical Lead, Advanced Air Interface, Future Wireless,
InterDigital Communications<EMAIL_ADDRESS>
###### Abstract

The development of the 5G new radio specifications has been derived by the the
deterministic low latency use-cases such as the ultra-reliable and low-latency
communications (URLLC). A URLLC application requires a stringent radio latency
and reliability performance, e.g., one-way radio latency of 1 ms with 99.999%
success probability. Furthermore, there is a concurrent progressive demand for
broadband capacity cellular applications, e.g., enhanced mobile broadband
(eMBB) use-cases. The coexistence among the URLLC and eMBB service classes
over a single radio spectrum is a challenging task since achieving the tight
URLLC radio targets typically results in a capacity loss. Hence, it is vital
for telecom operators to understand the capacity cost of fulling the various
URLLC requirements in order to sufficiently plan the corresponding pricing
models. Hence, in this work, a comprehensive analysis of the system capacity
loss is presented to achieve the various requirements of the different URLLC
use-cases. An extensive set of realistic system level simulations is performed
and introduced where valuable insights and system design recommendations on
the URLLC-eMBB quality of service coexistence are presented.

Index Terms— 5G new radio; Indoor factory automation (InF); URLLC; eMBB; 3GPP.

##  I Introduction

Unlike the former radio network generations, the fifth generation (5G) new
radio (NR) specifications support revolutionary multi quality-of-service (QoS)
classes such as the ultra-reliable low-latency communications (URLLC), and the
enhanced mobile broadband (eMBB) [1]. The URLLC applications demand a
stringent set of the radio latency and reliability targets while the eMBB
services demand broadband communication data rates [2]. This multi-QoS-class
coexistence enables novel cellular applications and use-cases, e.g., tactile
internet, true smart infrastructures, and virtual reality communications [3].

However, the deterministic low-latency industrial applications, i.e., industry
4.0 deployments [4, 5], are the lead drivers of the 5G and 5G-advanced
(release-18 and beyond) developments. Within those deployments, a wide range
of; though, conflicting, QoS classes, radio performance requirements, and
device types should be supported. That is, a deterministic ultra-low radio
latency with a certain reliability level is always guaranteed while a pre-
determined capacity/data rate target(s) should be preserved for simultaneous
capacity demanding services [6]. Fundamentally, achieving an ultra-low radio
latency leads to an insufficient radio capacity of the same communication
bandwidth. Thus, it is of a significant importance for telecom operators and
cellular service providers to identify the system capacity cost, which is
paid/lost in order to fulfill a certain ultra-low latency service performance
target, and hence, adapt their pricing models accordingly.

In the state-of-the-art open literature, there are many contributions of MAC
and PHY schemes that trade-off the achievable system capacity with the lowest
guaranteed radio latency. Examples include reinforcement learning based link
adaptation for a faster URLLC packet transmission [7, 8], QoS-aware scheduling
schemes [9], and on-the-fly resource preemption techniques [1, 10]. Those
proposals seek to achieve a guaranteed maximum radio latency for the URLLC
critical services while maximizing the overall system capacity, i.e., by
incurring the minimum possible capacity loss of the eMBB traffic.

In this paper, a comprehensive analysis of the capacity loss, due to
fulfilling a certain URLLC latency and reliability performance target, is
presented. Various URLLC latency and reliability targets are considered to
match the different URLLC use-cases and deployments. The reference case
considered in this work is the best effort eMBB, where the corresponding radio
latency performance is relaxed. Thus, the paper answers the following
question: ”what is the system capacity lost, compared to the best effort case,
to satisfy a guaranteed maximum URLLC radio latency with a certain radio
reliability level?”. An extensive set of realistic system level simulations is
performed to obtain a set of statistically reliable results and the respective
conclusions.

This paper is organized as follows. Section II presents the system model and
the main performance indicator considered in this work. Section III introduces
the simulation methdology and the performance results. Section IV presents the
acknowledgments and Section IV concludes the paper.

##  II System Model

###  II-A Setting the scene

In this work, we consider an industrial factory automation deployment with $C$
indoor cells, horizontally inter-distanced by $d=20$ meters, as shown by Fig.
1. Each cell serves an average number of $K$ uniformly distributed user-
equipment’s (UEs), where even UE dropping for all cells is adopted. In this
study, we assume only downlink (DL) traffic towards active UEs. To emulate the
URLLC use cases, the traffic is characterized by the FTP3 traffic model, where
data packets of a finite payload size B-Bytes are considered. The DL packets
arrive at the cell following a Poisson Arrival Process with a predefined mean
arrival rate $\lambda$ (packets/sec). Therefore, the total offered load
$\Omega$ in bits/sec for the entire factory deployment is calculated as:

Figure 1: Industrial factory deployment.

$\Omega=C\times\left(K\times B\times 8\times\lambda\right).$ (1)

We follow the 3GPP assumptions and guidelines for URLLC simulations. URLLC UEs
are dynamically multiplexed using the orthogonal frequency division
multiplexing (OFDMA). The 30 kHz numerology (i.e., sub-carrier spacing) is
used in this work, in line with [1], since it offers a sufficiently short OFDM
symbol duration in order to fulfill the stringent radio latency requirements
of the URLLC services. The minimum schedulable resource unit is the physical
resource block (PRB), consisting of 12 successive sub-carries. Moreover, we
consider a short transmission time interval (TTI) duration of 4 OFDM symbols
for faster URLLC transmission and scheduling.

Thus, an arriving DL packet is processed as follows: the packet is first
received and prepared for transmission by the serving cell RF stack processor.
The delay consumed before the DL transport block is ready for transmission is
taken explicitly into account in this study, in line with [2]. Then, the
dynamic scheduler of the serving cell multiplexes all pending DL packets for
transmission, using the proportional fair (PF) scheduling criterion. Hence,
packets of the highest PF metric, are transmitted during the current TTI
interval while other packets are further buffered towards the upcoming
scheduling instants, i.e., queuing delay of the dynamic scheduling. We
consider an adaptive modulation and coding selection (MCS) corresponding to a
first-transmission block error rate (BLER) of 1%. Therefore, the MCS level is
dynamically selected based on the latest reported channel quality indication
(CQI) from each UE.

Accordingly, the serving cell sends a scheduling grant, using the lower-layer
downlink control information (DCI) signaling, to notify the corresponding UEs
of the DL resource allocations. The resource overhead of the physical-layer
control signaling is explicitly taken into account; however, in this work, we
assume the delay of transmitting and processing the scheduling grants is
negligible. At the UE side, the delay consumed to process and decode received
DL packets is explicitly considered into account of the total radio latency.
If the transmitted DL packet is not successfully decoded at the UE, the UE
sends a hybrid automatic repeat request (HARQ) negative acknowledgment (NACK).
The delay for transmitting the HARQ feedback for each packet is also
explicitly considered. In such case, the serving cell re-transmits the
respective packet for soft-combining at the UE side. At the cell side, the
HARQ re-transmissions are always prioritized over new packet transmission in
order to reduce the overall radio latency for pending packets.

###  II-B Main performance metric indicators (KPIs)

The major KPIs of this work are the achievable one-way URLLC outage latency
and the corresponding network throughput/capacity, respectively. Thus, the
total radio latency of each DL packet $\varphi$ is measured and tracked, for
various outage probabilities $\rho$. In particular, the considered latency
metric denotes the delay from the moment a DL packet is generated and arrives
at the packet data convergence protocol (PDCP) layer of the serving cell until
it is successfully decoded at the intended UE. This sums up the cell and UE
processing delays, the transmission delay, the dynamic scheduling queuing
delay, and HARQ re-transmission delay, respectively. Therefore, for a certain
target outage probability, we present the maximum supported offered
load/capacity of the network such as to fulfill a maximum guaranteed URLLC
radio latency.

Furthermore, in this work, we adopt the throughput cost metric $\mathit{\Psi}$
of the URLLC and best effort (BE) cases, respectively. The BE case denotes the
transmissions of an infinite packet size and without a target outage latency
requirement such that the network capacity is maximized. Hence, the cost
metric implies how much network throughput is lost (paid) in order to fulfill
a certain URLLC outage latency target, compared to the achievable network
throughput of the BE case, and is calculated as follows:

$\mathit{\Psi}\left(B,\varphi,\rho\right)=\left(1-\frac{\mu_{urllc}\left(B,\varphi,\rho\right)}{\mu_{BE}}\right)\times
100\%,$ (2)

where $\mathit{\Psi}\left(B,\varphi,\rho\right)$ denotes the inflicted network
throughput cost of the URLLC use cases, with a packet size of B-Bytes in order
to fulfill a maximum guaranteed radio latency of $\varphi$ ms for the target
outage probability $\rho$. $\mu_{urllc}$ and $\mu_{BE}$ are the achievable
mean network throughput metrics of the URLLC and BE cases, respectively.

##  III Performance Evaluation

###  III-A Simulation Methodology

We comprehensively evaluate the performance of the URLLC service class using
extensive system level simulations. The major system settings are presented in
Table I. The conducted simulations follow the system modeling assumptions,
presented in Section II, and the general 3GPP simulation methodology. The
industrial factory deployment layout is considered alongside with employing
the state of the art industrial factory channel model. Dynamic FTP3 traffic
arrivals are considered at each UE, where the arrival rate follows a Poisson
Point Arrival Process. When a UE is created, it connects to the surrounding
cell with the highest received reference signal power (RSRP). The simulations
explicitly include the major functionalities of the PHY and MAC stack layers.
For each transmitted packet, the sub-carrier signal to interference noise
ratio (SINR) is calculated using the linear minimum mean squared error
interference rejection and combining (LMMSE-IRC) receiver. The effective SINR
is calculated by combining the estimated sub-carrier SINRs using the mean
mutual information per coded bit (MMIB) mapping. Based on the effective SINR,
the packet error probability (PEP) is calculated from predefined look-up
tables, obtained from extensive link level simulations. Therefore, based on
PEP, the packet is determined as successfully received or not. If the packet
is not successfully decoded from first transmission, the corresponding HARQ
re-transmissions are triggered. Every DL TTI, UEs are dynamically scheduled
based on the proportional fair criterion. The corresponding modulation and
coding scheme (MCS) is selected based on the latest available channel quality
indicator (CQI) reports.

Finally, the simulator is periodically calibrated by reporting and comparing
baseline performance statistics among the 3GPP partners. Moreover, in line
with [2], we run simulations for a sufficiently long period of time in order
to ensure high statistical confidence of the achievable results before drawing
solid conclusions. That is, the default simulation time is at least 5 million
successfully-decoded URLLC packets.

Table I: Simulation parameters. Parameter | Value
---|---
Environment | 3GPP-InF, one cluster, 12 cells
DL channel bandwidth | 40 MHz, SCS = 30 KHz, FDD
Channel model | InF-DH (dense clutter and high BS) [11]
BS transmit power | BS: 25 dBm
Carrier frequency | 4 GHz
BS heights | BS: 10m
Antenna setup | $2x2$
Average UEs per cell | 10
TTI configuration | 4-OFDM symbols
URLLC Traffic model | FTP3, paket size = 50/1500 Bytes
eMBB Traffic model | Best effort with infinite payload size
DL scheduling | proportional fair
Processing time | $\begin{array}[]{c}\textnormal{PDSCH prep. delay: 2.5-OFDM symbols}\\\ \textnormal{PDSCH decoding : 4.5-OFDM symbols}\\\ \end{array}$
DL receiver | L-MMSE-IRC

###  III-B Performance Results

In this section, we mainly evaluate the inflicted network capacity loss in
order to fulfill certain URLLC outage latency targets. In particular, the
maximum supported offered capacity is investigated for URLLC and best effort
cases, respectively. The URLLC services are evaluated to achieve various
outage latency targets. That is, for different maximum guaranteed radio
latency budgets for several outage probability levels, and using two payload
sizes (50 and 1500 Bytes). The best effort reference denotes the case where
the requirement on the radio latency is relaxed, and is evaluated under two
different UE schedulers: the proportional fair (PF) and equal throughput (ET)
schedulers, respectively. The former enables a fair scheduling criterion
without an inter-UE throughout regularization while the latter seeks to
achieve a guaranteed equal throughput per UE.

Fig. 2 depicts the maximum supported offered load of the URLLC traffic, with a
large 1500-byte payload size, to achieve a maximum guaranteed outage latency
of 1, 3 and 10 ms, respectively, and under various outage probabilities as
$10^{-5}$ and $10^{-2}$. As can be noticed, achieving the 1 ms latency
deadline with $10^{-5}$ outage probability significantly degrades the
achievable capacity compared to both cases of the best effort class, with PF
and ET schedulers, respectively. That is, to achieve a maximum guaranteed 1 ms
of the radio latency, the maximum supported offered loads are 2.03 Mbps and
11.65 Mbps for the $10^{-5}$ and $10^{-2}$ outage probabilities, compared to
71.06 Mbps and 93 Mbps with the best effort case of the ET and PF packet
schedulers. Furthermore, the best effort case with ET scheduler clearly
exhibits ~ 23.6% reduction of the achievable capacity compared to the case
with the PF scheduler, due to sacrificing part of the system capacity in order
to guarantee an equal-UE perceived throughput.

Figure 2: URLLC outage latency performance, large payload case.

Looking at the case of the URLLC small payload, Fig. 3 shows the achievable
load performance of the 50-byte URLLC case, where similar conclusions as of
Fig. 1 are observed. Relaxing the stringent requirement of the outage
probability offers a clear capacity improvement over the low radio latency
region. For instance, to achieve a maximum guaranteed 1-ms target, the
supported offered capacity is improved by ~ 76.4% with the $10^{-2}$ outage
probability, compared to that is of the $10^{-5}$ outage level.

Figure 3: URLLC outage latency performance, small payload case.

The packet size is deemed to have a vital impact on the achievable joint
outage latency and capacity performance. Therefore, Fig. 4 depicts a capacity
comparison of the URLLC class with small and large payload sizes, i.e., 50 and
1500 bytes, respectively, at the $10^{-5}$ outage probability. As clearly
seen, the small packet transmissions are observed to be more efficient in
supporting more offered capacity to achieve the same outage latency target,
compared to that is of the large payload. For instance, to achieve the
stringent 1-ms latency target, the small-payload transmissions supports ~
56.4% more offered load than the large-payload case.

Figure 4: URLLC outage latency performance, with payload size.

This is mainly attributed to the required multiple resource allocations over
multiple TTIs for the large-size packets to get fully transmitted. Fig. 5
depicts the empirical cumulative distribution function (ECDF) of the total
number of scheduled PRBs per each FTP3 packet, for different supported offered
load levels, which correspond to a certain maximum guaranteed outage latency,
i.e., 1, 1.5, and 3 ms radio latency targets, respectively, at the $10^{-5}$
outage probability. At the 50%ile level, a single 1500-byte URLLC packet
requires 63 and 78 PRBs (out of total 100 PRBs) of a single TTI, respectively,
for the three latency targets under evaluation to be fulfilled. This implies a
maximum of a single packet transmission per TTI while the other concurrent
packets from other active UE are being buffered towards the upcoming
transmission opportunities. However, at the tail distribution, i.e., 95%ile
level, for the offered loads of 70 and 220 Mbps, respectively, a single FTP3
packet requires 138 and 237 PRBs, respectively. This denotes that a single
1500-byte packet is scheduled over multiple TTIs, consuming most of the
available bandwidth. This is particularity relevant to the packets from the
cell-edge UEs, where their decoding ability is highly deteriorated due to the
strong inter-cell interference. Thus, to achieve the target 1% BLER, the
serving BS relies on a conservative MCS selection for the respective packet
transmissions, leading to a significant degradation of the network spectral
efficiency.

Figure 5: Resource utilization performance of the URLLC large packets, for
different outage latency targets.

On the opposite side, and as shown by Fig. 6, the small-payload URLLC
transmissions, i.e., 50-byte payload, are more efficient than the large
1500-byte payload from the minimum-delay allocation perspective. To fulfill
the stringent 1-ms latency target at the $10^{-5}$ outage probability, the
small-payload transmissions require 3 PRBs on average compared to 63 PRBs with
the large payload. This allows for: (a) scheduling the entire small-payload
URLLC packets within a single TTI duration without segmentation, and (b) co-
scheduling multiple packets from multiple active UEs at the same TTI, hence,
reducing the average packet queuing delay accordingly.

Figure 6: URLLC resource utilization performance, for small and large payload
packets.

Finally, Fig. 7 and 8 show the overall throughput cost of the URLLC service
class, compared to the best effort case. We consider various outage levels and
latency budgets of the large and small payload URLLC cases, respectively. In
particular, such cost denotes how much network throughput is lost (paid),
relative to the best effort case, in order to fulfill the corresponding
stringent URLLC outage latency target. For example, with the large-payload
URLLC case in Fig. 7, the throughput cost spans the range from 20% up to 97%.
To fulfill a maximum guarantee latency of 0.5 ms while employing the URLLC
large payload, the supported network offered load is 97% less than that is of
the best effort case. This is majorly because, with the large payload URLLC
transmissions, and to fulfill such stringent latency deadline, URLLC packets
must be scheduled within a single TTI duration, and without further packet
queuing or re-transmissions. Therefore, a single large-payload URLLC packet
scheduling consumes the entire bandwidth, i.e., the majority of the available
PRBs. Hence, the maximum supported offered capacity is drastically reduced
compared to the best effort case. Similar conclusions can be drawn from Fig.
8, for the case of the URLLC small payload size, i.e., 50 Bytes.

Figure 7: Spectral efficiency cost of achieving the URLLC outage targets,
large payload size. Figure 8: Spectral efficiency cost of achieving the URLLC
outage targets, small payload size.

##  IV Acknowledgments

This work has been partially completed while the main author was with Nokia
Bell Labs, Standardization and Research Access Lab, Aalborg, Denmark.

##  V Concluding Remarks

In this paper, a comprehensive analysis of the system capacity loss, due to
satisfying the various URLLC targets, is presented by means of extensive
system level simulations. The paper answers the following research question:
”what is the network capacity loss, compared to best effort case, due to
fulfilling a guaranteed maximum URLLC latency and reliability target?”, and
hence, the paper offers valuable insights for telecom operators to adapt their
pricing models within their multi-QoS 5G deployments. Various URLLC
configurations and settings are considered to match the different realistic
URLLC use cases, and the achievable system capacity is compared to that is of
the best-effort case, where the radio latency performance is relaxed.

## References

  * [1] A. A. Esswie and K. I. Pedersen, "Opportunistic Spatial Preemptive Scheduling for URLLC and eMBB Coexistence in Multi-User 5G Networks," in IEEE Access, vol. 6, pp. 38451-38463, 2018\.
  * [2] A. A. Esswie and K. I. Pedersen, "On the Ultra-Reliable and Low-Latency Communications in Flexible TDD/FDD 5G Networks," _in Proc. IEEE CCNC_ , 2020, pp. 1-6.
  * [3] A. Ghosh, A. Maeder, M. Baker and D. Chandramouli, "5G Evolution: A View on 5G Cellular Technology Beyond 3GPP Release 15," in IEEE Access, vol. 7, pp. 127639-127651, 2019\.
  * [4] M. Gundall et al., "Introduction of a 5G-Enabled Architecture for the Realization of Industry 4.0 Use Cases," in IEEE Access, vol. 9, pp. 25508-25521, 2021.
  * [5] W. Trneberg et al., "Towards Intelligent Industry 4.0 5G Networks: A First Throughput and QoE Measurement Campaign," _in Proc. IEEE SoftCOM_ , 2020, pp. 1-6.
  * [6] A. A. Esswie and K. I. Pedersen, "Analysis of Outage Latency and Throughput Performance in Industrial Factory 5G TDD Deployments," _in Proc. IEEE VTC-Spring_ , 2021, pp. 1-6.
  * [7] A. A. Esswie, K. I. Pedersen and P. E. Mogensen, "Online Radio Pattern Optimization Based on Dual Reinforcement-Learning Approach for 5G URLLC Networks," in IEEE Access, vol. 8, pp. 132922-132936, 2020.
  * [8] C. She et al., "A Tutorial on Ultrareliable and Low-Latency Communications in 6G: Integrating Domain Knowledge Into Deep Learning," in Proceedings of the IEEE, vol. 109, no. 3, pp. 204-246, March 2021.
  * [9] L. Wang and H. Zhang, "Analysis of joint scheduling and power control for predictable URLLC in industrial wireless networks," _in Proc. IEEE ICII_ , Orlando, FL, USA, 2019, pp. 160-169
  * [10] M. Morcos, M. Mhedhbi, A. Galindo-Serrano and S. Eddine Elayoubi, "Optimal resource preemption for aperiodic URLLC traffic in 5G Networks," _in Proc. IEEE PIMRC_ , 2020, pp. 1-6.
  * [11] _Study on channel model for frequencies from 0.5 to 100 GHz_ ; Release 16, 3GPP, TR 38.901, V16.1.0, Dec. 2019.