# Domain-aware Control-oriented Neural Models for Autonomous Underwater
Vehicles

Wenceslao Shaw Cortez  Soumya Vasisht  Aaron Tuor  Ján Drgoňa  Draguna Vrabie
Pacific Northwest National Laboratory, Richland, WA 99354 USA (e-mail:
{w.shawcortez, soumya.vasisht, aaron.tuor, jan.drgona,
draguna.vrabie}@pnnl.gov).

###### Abstract

Conventional physics-based modeling is a time-consuming bottleneck in control
design for complex nonlinear systems like autonomous underwater vehicles
(AUVs). In contrast, purely data-driven models, though convenient and quick to
obtain, require a large number of observations and lack operational guarantees
for safety-critical systems. Data-driven models leveraging available partially
characterized dynamics have potential to provide reliable systems models in a
typical data-limited scenario for high value complex systems, thereby avoiding
months of expensive expert modeling time. In this work we explore this middle-
ground between expert-modeled and pure data-driven modeling. We present
control-oriented parametric models with varying levels of domain-awareness
that exploit known system structure and prior physics knowledge to create
constrained deep neural dynamical system models. We employ universal
differential equations to construct data-driven blackbox and graybox
representations of the AUV dynamics. In addition, we explore a hybrid
formulation that explicitly models the residual error related to imperfect
graybox models. We compare the prediction performance of the learned models
for different distributions of initial conditions and control inputs to assess
their accuracy, generalization, and suitability for control.

###### keywords:

Deep learning, System identification, Nonlinear predictive control

††thanks: The research described in this work is part of the Mathematics for
Artificial Reasoning in Science Initiative at Pacific Northwest National
Laboratory. It was conducted under the Laboratory Directed Research and
Development Program at PNNL, a multiprogram national laboratory operated by
Battelle for the U.S. Department of Energy under contract DE-AC05-76RL01830.

##  1 Introduction

Underwater vehicle control is a vital aspect of maritime applications, such as
marine environmental monitoring, seabed mapping, geohazard assessment, and
remote exploration and sensing (Wynn et al., 2014; Yu et al., 2002). Control
design for autonomous underwater vehicles (AUVs) is a complex task due to the
underactuated nonlinear dynamics and the extreme operating environment. Most
applications are often restricted to simplistic control strategies that are
lacking in speed, versatility, and reliable autonomy in extreme environments
and situations. Advanced control methods such as predictive or adaptive
control require reliable system models for prediction and to evaluate system
properties such as stability and safety. Classical system identification
techniques that rely on transfer function learning or parameter estimation of
linear system models (Ljung, 2010; Van Overschee and De Moor, 1994) are
inefficient and non-trivial for complex systems.

Methods that discover governing equations underlying a dynamical system
directly from system measurement data (Brunton et al., 2016; Raissi et al.,
2017) have revolutionized the domain of scientific machine learning and data-
driven dynamics discovery. Physics-informed neural networks (PINNs) (Raissi et
al., 2019) incorporate physical domain knowledge in the form of soft
constraints to solve or obtain parameters of nonlinear partial differential
equations (PDEs) from observed noisy data using automatic differentiation for
training the deep learning models. However, PINNs have been shown to fail for
complex problems, making the loss function harder to optimize (Krishnapriyan
et al., 2021). Moreover, PINNs have been predominantly used for learning the
solutions to PDEs for a single initial condition, and thus are not directly
suitable for use in control that requires a generative model trained over a
distribution of initial conditions.

Deep learning-based alternatives to traditional system identification and
control have become a leading research avenue due to their rich expressive
power and their ability to capture complex dynamical interactions from data.
Neural state-space models (NSSMs) represent an example of such methods that
use discrete-time formulations (Rangapuram et al., 2018). Skomski et al.
(2021) showed that NSSMs can embed domain knowledge, including structural
priors and block nonlinearities. These models have been shown to work well in
predictive control and differentiable programming settings (Drgoňa et al.,
2021). A continuous-time and depth alternative to popular deep learning models
such as residual networks (ResNets) was introduced in Chen et al. (2018)
called neural ordinary differential equations (NODEs). NODEs have garnered
immense interest for their desirable properties such as invertibility (Deng et
al., 2020), suitability for control (Asikis et al., 2022), and parameter
efficiency and have been found more robust (Yan et al., 2019) than other deep
neural network model variants. Primarily used for developing black-box
representations, NODEs can now be augmented to include domain priors and
system structure via universal differential equations (UDEs) (Rackauckas et
al., 2020), making them an attractive proposition over other identification
methods. Hybrid neural networks form another class of flexible nonlinear
models that combine a given, partially correct model with an additional neural
network block to estimate the unmeasured, unmodeled dynamics that are
difficult to capture from first principles or whitebox models Psichogios and
Ungar (1992). This model is intuitive and easier to analyze, and has better
prediction performance than the standard single blackbox or graybox models and
has been applied in a number of domains (Meleshkova et al., 2021; Gutiérrez
and Hervás-Martínez, 2011).

In this work, we explore control-oriented NODE models for autonomous systems,
such as AUVs, and study the effects of increasing physical system knowledge
ranging from “blackbox” models with no embedded physics knowledge to “graybox”
models which rely on varying degrees of physical knowledge. We consider six
model types, including the vanilla NODE blackbox model, with and without
domain priors and boundary constraints, the NODE graybox model, and multiple
hybrid models with a fixed graybox block with varying parameter uncertainties
augmented with a blackbox block that learns the residual errors. The key
contributions include:

  * •

Development of a suite of domain-aware differentiable neural system models for
AUVs

  * •

Empirical analysis of the effect of applying domain-aware priors and system
knowledge during training on open-loop prediction performance.

We begin by introducing the ground truth AUV system model in Section 2, which
is used for data generation and model evaluation. Section 3 describes the
different system models and the training procedure. We present results and
discuss the model performance in Section 4.

##  2 Problem Setup

The problem addressed here is to develop a system identification approach for
AUVs using NODEs that accommodates varying degrees of a priori information
available for the system. Such system information considered here includes
dynamical properties, such as state constraints and stability-related
properties, and system model structure. In this section, we present the
dynamical model of an AUV as well as the development of the data set used for
training in the system identification approach.

###  2.1 AUV dynamic model

We consider a torpedo-style vehicle, such as the Iver family of AUVs from
OceanServer Technology (2022), to generate data for our modeling and analysis.
The nonlinear kinematic and hydrodynamic models we use are extended from
Stankiewicz et al. (2021) and summarized here.

$\displaystyle\dot{p}_{x}=u\cos{\psi}\cos{\theta}+w\cos{\psi}\sin{\theta},$
(1a)
$\displaystyle\dot{p}_{y}=u\sin{\psi}\cos{\theta}+w\sin{\psi}\sin{\theta},$
(1b) $\displaystyle\dot{p}_{z}=w\cos{\theta}-u\sin{\theta},$ (1c)
$\displaystyle\dot{\theta}=q,$ (1d) $\displaystyle\dot{\psi}=r/\cos{\theta},$
(1e) $\displaystyle\dot{u}=\overline{X}_{uu}u^{2}+\overline{k}\delta_{u_{c}},$
(1f) $\displaystyle\dot{w}=\overline{Z}_{w|w|}w|w|+\overline{WB}\cos{\theta},$
(1g)
$\displaystyle\dot{q}=\overline{M}_{uq}uq+\overline{M}_{q}q-\overline{Bz_{B}}\sin{\theta}+\bar{b}u^{2}\delta_{q_{c}},$
(1h) $\displaystyle\dot{r}=\overline{N}_{ur}ur+\bar{c}u^{2}\delta_{r_{c}}$
(1i)
$\displaystyle\dot{\delta}_{u_{c}}=\overline{K}_{\delta_{u}}(\delta_{u}-\delta_{u_{c}})$
(1j)
$\displaystyle\dot{\delta}_{q_{c}}=\overline{K}_{\delta_{q}}(\delta_{q}-\delta_{q_{c}})$
(1k)
$\displaystyle\dot{\delta}_{r_{c}}=\overline{K}_{\delta_{r}}(\delta_{r}-\delta_{r_{c}}).$
(1l)

The state is x $=$ $(p_{x},p_{y},p_{z},\theta,\psi,u,w,q,r,$ $\delta_{u_{c}},$
$\delta_{q_{c}},$ $\delta_{r_{c}})$, where
$(p_{x},p_{y},p_{z})\in\mathbb{R}^{3}$ are the respective north, east, down
inertial coordinates, $\theta\in[-\pi/2,\pi/2],\psi\in\textsf{S}^{1}$ are the
respective pitch and yaw orientation angles (Euler angles), $u,w\in\mathbb{R}$
are the respective surge and heave velocities, $q,r\in\mathbb{R}$ are the
pitch and yaw angular velocities, and
$\delta_{u_{c}}\in[0,1],\delta_{q_{c}},\delta_{r_{c}}\in[-1,1]$ are the
delayed control inputs for thrust, elevator deflection, and rudder deflection,
respectively. Note that $(u,w,q,r)\in\mathbb{R}^{4}$ are defined in the body
frame. Similar to Stankiewicz et al. (2021), we assume roll and sway
velocities are negligible. We also assume the vehicle has been trimmed to be
neutrally buoyant and have no heave velocities, i.e. $w\equiv 0$. Refer to
Stankiewicz et al. (2021) for the hydrodynamic coefficients
$(\overline{\cdot})$ of the Iver3 AUV. The control inputs
$\textbf{u}=(\delta_{u},\delta_{q},\delta_{r})\in\mathcal{U}\mathrel{\mathop{\mathchar
58\relax}}=[0,1]\times[-1,1]\times[-1,1]$ are the normalized thrust,
deflections in the diving plane and rudder fins, respectively. The system (1)
considers input delays that are common in practice as a low level controller
tracks the given input commands. Here the input delay is modeled as a linear
tracking control law as seen in the dynamics of
$\delta_{u_{c}},\delta_{q_{c}}$, and $\delta_{r_{c}}$ with
$\overline{K}_{\delta_{u}}=\overline{K}_{\delta_{q}}=\overline{K}_{\delta_{r}}=-10$.

A typical AUV hardware setup includes sensors like an IMU to measure angular
rates, a compass for heading information, an acoustic doppler current profiler
(ADCP) for measuring the vehicle’s velocity through water and a depth sensor.
We assume that the positions in the inertial coordinates are not directly
measurable. Finally, the actuators contain sensors to determine the true
actuator values in order to track the commanded input. Consistent with this
realistic scenario, we consider the vehicle output to be velocities, angular
rates and heading,
$\mathbf{y}=(\theta,\psi,u,q,r,\delta_{u_{c}},\delta_{q_{c}},\delta_{r_{c}})\in\mathcal{M}\mathrel{\mathop{\mathchar
58\relax}}=[-\pi/2,\pi/2]\times\mathbb{S}^{1}\times\mathbb{R}^{3}\times[0,1]\times[-1,1]\times[-1,1]$,
where $\mathbf{y}$ consists of $n_{y}=8$ elements.

###  2.2 Dataset

For the modeling process, we generate a dataset $\mathcal{D}$ of input-output
pairs obtained by applying appropriate excitation signals to explore the AUV
state and input space:

$\mathcal{D}=\\{[(\mathbf{u}_{0},\mathbf{y}_{0}),(\mathbf{u}_{\delta},\mathbf{y}_{\delta}),\ldots,(\mathbf{u}_{N\delta},\mathbf{y}_{N\delta})]^{(i)}\\},$
(2)

with sampling time $\delta\in\mathbb{R}$ and data batches $i=1,2,\ldots,n$, of
length $N$ each, and let $\mathbf{u}_{k}$ denote the vector $\mathbf{u}$ at
time instant $k$. We use the mathematical model in (1) representing the AUV
dynamics to generate the training and testing datasets and for ground truth
trajectory generation for model evaluation. We sample control input
trajectories and initial conditions from a distribution of pre-defined
trajectories and initial conditions based on domain knowledge, that generates
a combination of steady flight, diving, climbing, and port and starboard side
turning maneuvers. Each input trajectory is constructed by combining a
sequence of $50$ base trajectories for a length of $\frac{N}{50}$ seconds
total per batch. Each base trajectory is randomly selected from a step
response, periodic response, or spline constructed such that
$\textbf{u}_{\tau}\in\mathcal{U}$ for all $\tau\in[0,N\delta]$. The training
set was determined by integrating (1) over a horizon of $N\delta$ s for
$\delta=0.01$ using ‘odeint’ from Python’s scipy module.

##  3 Data-driven Models and Training

In this section we present the suite of models and how they are trained to
approximate the true system (1). The suite of models lie on a spectrum of
domain-aware representations to reflect real-world scenarios depending on what
degree of model knowledge is known a priori. The models are implemented as
NODE.

###  3.1 Blackbox model

The blackbox model assumes minimal knowledge of the system dynamics. Here the
dynamics are modelled with the neural network mapping
$\textbf{f}_{W}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{M}\times\mathcal{U}\to\mathbb{R}^{n_{y}}$, where $W$
represents the weights and biases of a multi-layer perceptron (MLP) with four
hidden layers of size 128 each. The blackbox dynamical system is defined by:

$\dot{\textbf{z}}(t)=\textbf{f}_{W}(\textbf{z}(t),\textbf{u}(t)),\
\textbf{z}(0)=\textbf{z}_{0}\in\mathcal{M},$ (3)

where $\mathbf{z}\in\mathbb{R}^{n_{y}}$ is the output of the model.

###  3.2 Graybox model

The graybox model requires the most knowledge of the system dynamics as it
assumes the structure of (1) is known, but the model parameters are uncertain.
Let the vector of model parameters be:
$\boldsymbol{\mu}=(\overline{X}_{u}u,\overline{k},\overline{M}_{uq},\overline{M}_{q},\overline{B_{z_{B}}},\overline{b},\overline{N}_{ur},\overline{c})\in\mathbb{R}^{n_{\mu}}$,
and let
$\hat{\boldsymbol{\mu}}=(\hat{X}_{u}u,\hat{k},\hat{M}_{uq},\hat{M}_{q},\hat{B_{z_{B}}},\hat{b},\hat{N}_{ur},\hat{c})$
denote the approximate model parameters. Furthermore, let the approximate
model dynamics be:

$\textbf{f}_{\hat{\boldsymbol{\mu}}}(\textbf{y},\textbf{u})=\begin{pmatrix}q\\\
r/\cos{\theta}\\\ \hat{X}_{uu}u^{2}+\hat{k}\delta_{u}\\\
\hat{M}_{uq}uq+\hat{M}_{q}q-\hat{Bz_{B}}\sin{\theta}+\hat{b}u^{2}\delta_{q}\\\
\hat{N}_{ur}ur+\hat{c}u^{2}\delta_{r}\\\
\dot{\delta}_{u_{c}}=\hat{K}_{\delta_{u}}(\delta_{u}-\delta_{u_{c}})\\\
\dot{\delta}_{q_{c}}=\hat{K}_{\delta_{q}}(\delta_{q}-\delta_{q_{c}})\\\
\dot{\delta}_{r_{c}}=\hat{K}_{\delta_{r}}(\delta_{r}-\delta_{r_{c}})\end{pmatrix}.$
(4)

The graybox dynamical system is defined as:

$\dot{\textbf{z}}(t)=\textbf{f}_{\hat{\boldsymbol{\mu}}}(\textbf{z}(t),\textbf{u}(t)),\
\textbf{z}(0)=\textbf{z}_{0}\in\mathcal{M}.$ (5)

Assuming some degree of informed initial guess by the modeller, the graybox
model parameters, $\hat{\boldsymbol{\mu}}$, are initialized to within an
offset of the true model parameters, $\boldsymbol{\mu}$. Let $e_{\mu}\in[0,1]$
represent the error in the model parameters. For our experiments, parameters
are sampled uniformly from a set whose range is defined by $e_{\mu}$. For
example,
$\hat{X}_{uu}\sim\mathscr{U}([\min\\{X_{uu}(1-e_{\mu}),X_{uu}(1+e_{\mu})\\},\max\\{X_{uu}(1-e_{\mu}),X_{uu}(1+e_{\mu})\\}]$,
where $\mathscr{U}(S)$ denotes a uniform sampling of the set $S$. To state
this more formally, let
$\mathcal{S}^{\mu}_{e_{\mu}}=\mathcal{S}^{\mu^{1}}_{e_{\mu}}\times...\times\mathcal{S}^{\mu^{n_{\mu}}}_{e_{\mu}}$,
where
$\mathcal{S}^{\mu^{i}}_{e_{\mu}}=[\min\\{\mu_{i}(1-e_{\mu}),\mu_{i}(1+e_{\mu})\\},\max\\{\mu_{i}(1-e_{\mu}),\mu_{i}(1+e_{\mu})\\}]$
with $\mu_{i}$ as the $i$th element of $\boldsymbol{\mu}$. Thus for a given
${e_{\mu}}$,
$\hat{\boldsymbol{\mu}}\sim\mathscr{U}(\mathcal{S}^{\mu}_{e_{\mu}})$ implies
an estimated parameter uniformly sampled from a range $e_{\mu}$ of the true
model parameters. For the graybox model, the model parameters were initialized
as $\hat{\boldsymbol{\mu}}_{0}\sim\mathcal{S}^{\mu}_{e_{\mu}=1.0}$ to provide
the maximum offset of the true values and test if the estimated model
parameters can converge to the true model parameters.

###  3.3 Hybrid model

The hybrid model combines an uncertain approximation of the system dynamics,
similar to the graybox, with a blackbox model to compensate for residual
errors. This model considers the case when some, imperfect knowledge, of the
system is known a priori. Here the model parameters of the graybox model,
$\hat{\boldsymbol{\mu}}$, are initialized and fixed to a random offset of
their true parameters. The approximate model used is the same as the graybox
$\textbf{f}_{\hat{\boldsymbol{\mu}}}(\textbf{y},\textbf{u})$ from (4). The
blackbox model is the same as from the previous case with the neural network
mapping $\textbf{f}_{W}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{M}\times\mathcal{U}\to\mathbb{R}^{n_{y}}$, with $W$
representing the weights and biases in the neural network. In this hybrid
case, only the blackbox parameters, $W$, are trained and not the offset model
parameters $\hat{\boldsymbol{\mu}}$. The hybrid model dynamics are defined by:

$\dot{\textbf{z}}(t)=\textbf{f}_{\hat{\boldsymbol{\mu}}}(\textbf{z}(t),\textbf{u}(t))+\textbf{f}_{W}(\textbf{z}(t)),\textbf{u}(t)),\
\textbf{z}(0)=\textbf{z}_{0}\in\mathcal{M}.$ (6)

To compare various different errors in the graybox portion of the hybrid
model, we look at three implementations of the hybrid model which are
dependent on the percentage of error in the model parameters,
$\hat{\boldsymbol{\mu}}$, which is defined by $e_{\mu}$. Here we look at the
three cases ‘$\text{hybrid}_{e_{\mu}=1.0}$’, ‘$\text{hybrid}_{e_{\mu}=0.5}$’,
and ‘$\text{hybrid}_{e_{\mu}=0.3}$’.

###  3.4 Training

We perform batched training on the trajectories in $\mathcal{D}$. Each model
was trained for $n=5$ batches. The length, $N$, was increased by a factor of 2
for each batch yielding batches with $N=100$, $N=200$, $N=400$, $N=800$, and
$N=1600$. The purpose of using multiple batches of different lengths is to
prevent over-approximation over the same time length and it extends how long
the control inputs are applied to the system, which yields different system
responses. Training is performed using the Neuromancer differentiable
programming framework Tuor et al. (2022), wherein the following optimization
problem is solved to train the model parameters $\bf V\in\mathbb{R}^{n_{v}}$:

$\displaystyle\min_{{\bf V}}$
$\displaystyle\frac{1}{mN}\sum_{i=1}^{m}\sum_{k=0}^{N-1}\big{(}\ell_{\texttt{mse}}({\bf
z}_{k}^{i},{\bf y}_{k})+p_{z}(c({\bf z}_{k}^{i})\big{)}$ (7a) s.t.
$\displaystyle{\bf z}_{k+1}^{i}=ODESolve({\bf f}_{\bf V^{i}}({\bf
z}_{k}^{i},{\bf u}_{k})),\ k\in\mathbb{N}_{0}^{N-1}$ (7b) $\displaystyle{\bf
z}_{0}^{i}={\bf y}_{0}$ (7c) $\displaystyle{\bf
V}^{i}\in\Xi\subset\mathbb{R}^{n_{v}}$ (7d)

where

$\ell_{\texttt{mse}}({\bf z},{\bf y})=\|{\bf y}-{\bf z}\|_{2}^{2}$ (8)

is the objective loss function, $c\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{n_{y}}\to\mathbb{R}^{n_{c}}$ represents state
constraints with the penalty functions $p_{x}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{n_{c}}\to\mathbb{R}$, $\Xi$ is a random distribution
from which the model parameters ${\bf V}^{i}$ are sampled from, and
$\textbf{f}_{\textbf{V}}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{M}\times\mathcal{U}\to\mathbb{R}^{n_{y}}$ is a general NODE
model parameterized by V. Stochastic gradient descent is used to minimize the
loss function (7a) over the distribution of model parameters (7d), where $m$
represents the total number of model parameter samples, and $i$ denotes the
index of the sample. We note that ODESolve can represent any ODE solver such
as Runge Kutta or Euler integrators. In our implementation, an Euler
approximation is used to integrate the NODE models. The problem (7) is solved
for each batch, and the final trained model is output after the final batch.
The penalty functions are used to enforce boundary conditions, which will be
discussed in the next section.

We emphasize that the different model types will have different model
parameters that are trained. In the case of the blackbox model,
$\textbf{f}_{\textbf{V}}=\textbf{f}_{\textbf{W}}$, where
$\textbf{V}=\textbf{W}$ are the MLP parameters that will be trained. For the
graybox model, $\textbf{f}_{\textbf{V}}=\textbf{f}_{\hat{\boldsymbol{\mu}}}$,
where $\textbf{V}=\hat{\boldsymbol{\mu}}$. Finally, for the hybrid model
$\textbf{f}_{\textbf{V}}=\\{\textbf{f}_{\hat{\boldsymbol{\mu}}},\textbf{f}_{W}$},
where $\textbf{V}=\textbf{W}$.

###  3.5 Learning with Boundary Conditions

The Neuromancer framework allows for enforcing boundary conditions in the form
of state constraints $c({\bf z})$. In the case of (1), it is clear that if
$\theta\to\pm\frac{\pi}{2}$, then the system dynamics become undefined at the
singularity. We note that the randomized training sets are restricted to
$\mathcal{M}$ so such instabilities are never reached, and thus we aim to
enforce constraints on the models to also avoid the singularity. Similarly,
the states associated with input delays should remain inside the input set
$\mathcal{U}$.

To enforce the boundary conditions, we define
$c_{i}(\mathbf{z})\mathrel{\mathop{\mathchar
58\relax}}=\max\\{0,-\mathbf{z}_{i}+\underline{z}_{i}\\}+\max\\{0,\mathbf{z}_{i}-\overline{z}_{i}\\}$,
where $i\in\mathbb{N}_{0}^{n_{y}-1}$, $\mathbf{z}_{i}$ is the $i$th element of
$\mathbf{z}$, $\underline{z}$ $=$
$(-\frac{\pi}{2},\linebreak[1]-\infty,\linebreak[1]-\infty,\linebreak[1]-\infty,\linebreak[1]-\infty,\linebreak[1]0,\linebreak[1]-1,\linebreak[1]-1)$
is the lower bound on $\mathbf{z}$, and
$\overline{z}=(\frac{\pi}{2},\infty,\infty,\infty,\infty,1,1,1)$ is the upper
bound on $\mathbf{z}$. Furthermore, let
$c(\mathbf{z})=(c_{0}(\mathbf{z}),...,c_{n_{y}-1}(\mathbf{z}))$. The penalty
function used is: $p_{c}(\textbf{c})=\|\textbf{c}\|_{2}$.

###  3.6 Learning with Domain-aware Priors

In addition to state constraints, the Neuromancer framework enables the
training process to enforce system-type properties of the trained model. As
with many physical systems, the AUV system is dissipative, which means that
the system has no internal production of energy, and the energy stored by the
system at any time does not exceed the energy externally supplied at that time
(Willems, 2007). The trained models can be trained to enforce the this
property. For the blackbox model and the blackbox component of the hybrid
model, if the 2-norm of the neural network weights and activation matrices are
restricted to be less than $1$, then the NODE is dissipative (Drgona et al.,
2022). Here, the neural network weights are restricted such that the
eigenvalues of the weighting and activation matrices lie within
$[\sigma_{min},\sigma_{max}]$ for
$\sigma_{min},\sigma_{max}\in\mathbb{R}_{\geq 0}$. This restriction is applied
during training to always ensure the dissipativity property.

For the blackbox model, the eigenvalues are set to $\sigma_{min}=0.5$,
$\sigma_{max}=1.0$. The reason for setting $\sigma_{min}=0.5$ instead of zero
is to keep the model away from steady-state solutions since the training set
was designed to excite the system and avoid such solutions. The blackbox
component of the hybrid models had eigenvalues set to $\sigma_{min}=0.0$,
$\sigma_{max}=1.0$. Here we allow the eigenvalues to reach zero unlike in the
blackbox model because in the hybrid model the blackbox component is used to
compensate for the residual error of the graybox. Thus if the graybox is able
to significantly reduce the residual error, the blackbox should approximate a
zero map. To show the effects of the boundary conditions and the domain-aware
priors on the neural network NODE model, we compare in our experiments a NODE
model without boundary conditions/domain-aware priors denoted “blackbox”, and
the same architecture with boundary conditions/domain-aware priors denoted
“constr. blackbox”

##  4 Results and Discussion

(a) Pitch angle ($\theta$) residual MSE (b) constr. blackbox (CBB)

(c) Surge velocity ($u$) residual MSE

(d) graybox (GB)

Figure 1: Comparison of MSE for (a) Pitch angle, $\theta$, (c) Surge velocity,
$u$, for each model type. The open-loop response, showing the measured states
of the best model of type (b) constrained blackbox (CBB) and (d) graybox (GB),
evaluated on test data. The ground truth and the predicted trajectories are in
cyan and purple, respectively.

For our empirical analysis, 10 sets of randomly sampled initial parameters for
each of the six model types, blackbox, constr. blackbox ,
$\text{hybrid}_{e_{\mu}=1.0}$, $\text{hybrid}_{e_{\mu}=0.5}$,
$\text{hybrid}_{e_{\mu}=0.3}$, and graybox, were fit to the training dataset
using the ADAM-W (Loshchilov and Hutter, 2017) variant of stochastic gradient
descent. After training these 10 models instances per model type, they were
evaluated on a held out test set. The test set was developed by randomly
initializing five initial conditions from $\mathcal{M}$ each for five randomly
generated input trajectories, which yields a total of 25 trajectories. The
input trajectories were developed in the same manner as the training set with
$N=5000$. We note that the test set is the same for all model types.

For each model type, there are 10 models that were each evaluated on 25
trajectories. For each trajectory of each model, the mean squared error (MSE)
is:

$\frac{1}{N}\sum_{k=0}^{k=N}\ell_{\texttt{mse}}({\bf z}^{\prime}_{k},{\bf
y}^{\prime}_{k})$ (9)

where ${\bf z}^{\prime}_{k},{\bf y}^{\prime}_{k}$ are the normalized outputs
of the model and true system respectively. The outputs are normalized using
the training set data for a unit-less error metric. We note that the inter
quartile range approach was used to remove outlier models from the results.
For each model type, we consider the ‘best’ trained model to be the one in
which the mean of the mean square error (MSE) over all 25 trajectories is the
lowest.

Table 1: Comparison of blackbox plain, blackbox, hybrid, and graybox models  model | Mean MSE | Standard deviation
---|---|---
blackbox | 16.459841 | 61.742187
constr. blackbox | 7.985731 | 0.311266
$\text{hybrid}_{e_{\mu}=1.0}$ | 82.699305 | 11399.445734
$\text{hybrid}_{e_{\mu}=0.5}$ | 6.01626 | 21.648181
$\text{hybrid}_{e_{\mu}=0.3}$ | 2.626124 | 0.474631
graybox | 0.025351 | 0.000446

Table 1 shows the results of the blackbox, constr. blackbox ,
$\text{hybrid}_{e_{\mu}=1.0}$, $\text{hybrid}_{e_{\mu}=0.5}$,
$\text{hybrid}_{e_{\mu}=0.3}$, and graybox models. Figure 1 interprets Table 1
visually and focuses on the pitch and surge velocity residual error evolution
over time and the standard deviation over the 25 trajectories, for each
learned model type. The open loop response of all measured states for the
constr. blackbox and the graybox models are also presented.

Several conclusions can be drawn from Table 1. First, the graybox outperforms
all model types with respect to the mean of all trained models and associated
standard deviation. This shows that on average, the graybox is the best model
type to use and the trained models will consistently yield the same
performance in terms of mean squared error from the true model dynamics. This
fits with intuition in that the graybox requires the most domain information
to describe the system structure exactly, but with uncertain parameters.
However this may not be true in many scenarios especially when simplifications
are made in the graybox model.

In practice, there usually will be unmodeled dynamics that need to be
accounted for. This can be represented by the hybrid model, where imperfect
model parameters will cause the graybox system dynamics to deviate with the
true system and requires the blackbox component to compensate. The hybrid
model results from Table 1 also follow intuition in that as $e_{\mu}\to 0.0$,
i.e., the graybox component approaches the true system, the blackbox component
performs better with respect to reducing the MSE between the hybrid model and
the true system model. Furthermore, the trained models become more consistent
in their performance as $e_{\mu}\to 0.0$, which is represented by the decrease
in standard deviation of the model performances.

The blackbox/blackbox constrained models represent the case when no model
knowledge is assumed and they must try to learn the true system dynamics.
According to Table 1, the blackbox and blackbox constrained models actually
outperform the $\text{hybrid}_{e_{\mu}=1.0}$ in terms of mean MSE and standard
deviation of all trained models. This suggests that the blackbox type models
are preferred to the hybrid model if the graybox component poorly represents
the true system. In this experiment, a poor representation is defined by
$e_{\mu}=1.0$, wherein the model parameters are sufficiently inaccurate as to
impede the blackbox component’s ability to compensate for the error between
the graybox and true system dynamics. One possible reason for this is that the
model parameters may cause the graybox model to approach an unstable system
dynamics, which would mean that the blackbox component is attempting to
compensate for the near-instability while also eliminating the residual error
between the graybox and true system.

Finally, we can see the comparison between blackbox and blackbox constrained
models. It is clear that on average, the blackbox constrained models yield
lower MSE and are more consistent with respect to the standard deviation of
all models. This suggests that the use of a priori information of the system
model, specifically the state constraints and dissipativity of the system,
help improve the training of the blackbox models and can reduce the errors of
the blackbox model by a factor of 2. This result fits with intuition in that
more information leads to improved system identification.

##  5 Conclusion

In this paper we present a domain-aware framework for system identification of
autonomous underwater vehicles using neural ordinary differential equations
(NODEs). The framework allows for structured training of the NODE model by
enforcing state constraints and domain-aware priors to ensure dissipativity of
the closed-loop system. The proposed approach considers a suite of models that
lie on a spectrum of domain-aware representations to reflect real-world
scenarios depending on what degree of model knowledge is known a priori. The
models range from a blackbox model, a hybrid model, and a graybox model. The
blackbox, graybox, and several implementations of the hybrid model were
compared in an experiment using training sets and evaluation sets. Future work
will investigate the effect of noise on training and incorporate the domain-
aware models in a predictive control structure.

## References

  * Asikis et al. (2022) Asikis, T., Böttcher, L., and Antulov-Fantulin, N. (2022).  Neural ordinary differential equation control of dynamics on graphs.  _Physical Review Research_ , 4(1), 013221.  10.1103/physrevresearch.4.013221.
  * Brunton et al. (2016) Brunton, S.L., Proctor, J.L., and Kutz, J.N. (2016).  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 113(15), 3932–3937.  10.1073/pnas.1517384113.
  * Chen et al. (2018) Chen, R.T., Rubanova, Y., Bettencourt, J., and Duvenaud, D.K. (2018).  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31.
  * Deng et al. (2020) Deng, R., Chang, B., Brubaker, M.A., Mori, G., and Lehrmann, A.M. (2020).  Modeling continuous stochastic processes with dynamic normalizing flows.  _CoRR_ , abs/2002.10516.
  * Drgoňa et al. (2021) Drgoňa, J., Tuor, A., Skomski, E., Vasisht, S., and Vrabie, D. (2021).  Deep learning explicit differentiable predictive control laws for buildings.  _IFAC-PapersOnLine_ , 54(6), 14–19.
  * Drgona et al. (2022) Drgona, J., Vasisht, S., Tuor, A., and Vrabie, D. (2022).  Dissipative Deep Neural Dynamical Systems.  URL https://arxiv.org/abs/2011.13492.
  * Gutiérrez and Hervás-Martínez (2011) Gutiérrez, P.A. and Hervás-Martínez, C. (2011).  Hybrid artificial neural networks: Models, algorithms and data.  In J. Cabestany, I. Rojas, and G. Joya (eds.), _Advances in Computational Intelligence_ , 177–184. Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Krishnapriyan et al. (2021) Krishnapriyan, A., Gholami, A., Zhe, S., Kirby, R., and Mahoney, M.W. (2021).  Characterizing possible failure modes in physics-informed neural networks.  _Advances in Neural Information Processing Systems_ , 34.
  * Ljung (2010) Ljung, L. (2010).  Perspectives on system identification.  _Annual Reviews in Control_ , 34(1), 1–12.  https://doi.org/10.1016/j.arcontrol.2009.12.001.
  * Loshchilov and Hutter (2017) Loshchilov, I. and Hutter, F. (2017).  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_.
  * Meleshkova et al. (2021) Meleshkova, Z., Ivanov, S.E., and Ivanova, L. (2021).  Application of neural ode with embedded hybrid method for robotic manipulator control.  _Procedia Computer Science_ , 193, 314–324.
  * OceanServer Technology (2022) OceanServer Technology (2022).  Iver3 Platform.  https://auvac.org/293-2/.  Accessed: 2022-06-11.
  * Psichogios and Ungar (1992) Psichogios, D.C. and Ungar, L.H. (1992).  A hybrid neural network-first principles approach to process modeling.  _AIChE Journal_ , 38(10), 1499–1511.
  * Rackauckas et al. (2020) Rackauckas, C., Ma, Y., Martensen, J., Warner, C., Zubov, K., Supekar, R., Skinner, D., Ramadhan, A., and Edelman, A. (2020).  Universal differential equations for scientific machine learning.
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., and Karniadakis, G. (2019).  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  _Journal of Computational Physics_ , 378, 686–707.
  * Raissi et al. (2017) Raissi, M., Perdikaris, P., and Karniadakis, G.E. (2017).  Inferring solutions of differential equations using noisy multi-fidelity data.  _Journal of Computational Physics_ , 335, 736–746.
  * Rangapuram et al. (2018) Rangapuram, S.S., Seeger, M.W., Gasthaus, J., Stella, L., Wang, Y., and Januschowski, T. (2018).  Deep state space models for time series forecasting.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), _Advances in Neural Information Processing Systems 31_ , 7785–7794. Curran Associates, Inc.
  * Skomski et al. (2021) Skomski, E., Vasisht, S., Wight, C., Tuor, A., Drgoňa, J., and Vrabie, D. (2021).  Constrained block nonlinear neural dynamical models.  In _2021 American Control Conference (ACC)_ , 3993–4000.  10.23919/ACC50511.2021.9482930.
  * Stankiewicz et al. (2021) Stankiewicz, P., Tan, Y.T., and Kobilarov, M. (2021).  Adaptive sampling with an autonomous underwater vehicle in static marine environments.  _Journal of Field Robotics_ , 38(4), 572–597.
  * Tuor et al. (2022) Tuor, A., Drgona, J., Skomski, M., Koch, J., Chen, Z., and Vrabie, D. (2022).  NeuroMANCER: Neural Modules with Adaptive Nonlinear Constraints and Efficient Regularizations.  URL https://github.com/pnnl/neuromancer.
  * Van Overschee and De Moor (1994) Van Overschee, P. and De Moor, B. (1994).  N4SID: Subspace algorithms for the identification of combined deterministic-stochastic systems.  _Automatica_ , 30(1), 75–93.  Special issue on statistical signal processing and control.
  * Willems (2007) Willems, J.C. (2007).  Dissipative dynamical systems.  _European Journal of Control_ , 13(2), 134–151.
  * Wynn et al. (2014) Wynn, R.B., Huvenne, V.A., Le Bas, T.P., Murton, B.J., Connelly, D.P., Bett, B.J., Ruhl, H.A., Morris, K.J., Peakall, J., Parsons, D.R., Sumner, E.J., Darby, S.E., Dorrell, R.M., and Hunt, J.E. (2014).  Autonomous underwater vehicles (auvs): Their past, present and future contributions to the advancement of marine geoscience.  _Marine Geology_ , 352, 451–468.
  * Yan et al. (2019) Yan, H., Du, J., Tan, V.Y.F., and Feng, J. (2019).  On robustness of neural ordinary differential equations.  _CoRR_ , abs/1910.05513.
  * Yu et al. (2002) Yu, X., Dickey, T., Bellingham, J., Manov, D., and Streitlien, K. (2002).  The application of autonomous underwater vehicles for interdisciplinary measurements in massachusetts and cape cod bays.  _Continental Shelf Research_ , 22(15), 2225–2245.


#  Enhancing Robotic Navigation: An Evaluation of Single and Multi-Objective
Reinforcement Learning Strategies

Vicki Young1, Jumman Hossain2, Nirmalya Roy2 2{jumman.hossain<EMAIL_ADDRESS>1Department of Computer Science, University of San Francisco, San Francisco,
USA  2Department of Information Systems, University of Maryland, Baltimore
County, USA<EMAIL_ADDRESS>
###### Abstract

This study presents a comparative analysis between single-objective and multi-
objective reinforcement learning methods for training a robot to navigate
effectively to an end goal while efficiently avoiding obstacles. Traditional
reinforcement learning techniques, namely Deep Q-Network (DQN), Deep
Deterministic Policy Gradient (DDPG), and Twin Delayed DDPG (TD3), have been
evaluated using the Gazebo simulation framework in a variety of environments
with parameters such as random goal and robot starting locations. These
methods provide a numerical reward to the robot, offering an indication of
action quality in relation to the goal. However, their limitations become
apparent in complex settings where multiple, potentially conflicting,
objectives are present. To address these limitations, we propose an approach
employing Multi-Objective Reinforcement Learning (MORL). By modifying the
reward function to return a vector of rewards, each pertaining to a distinct
objective, the robot learns a policy that effectively balances the different
goals, aiming to achieve a Pareto optimal solution. This comparative study
highlights the potential for MORL in complex, dynamic robotic navigation
tasks, setting the stage for future investigations into more adaptable and
robust robotic behaviors.

###### Index Terms:

Reinforcement Learning, Single-Objective Reinforcement Learning, Multi-
Objective Reinforcement Learning (MORL), Robotic Navigation, Dynamic
Environments.

##  I Introduction

Robotic navigation, particularly in dynamic and complex environments, is a
central problem in the field of robotics and artificial intelligence. With the
proliferation of robotics in a myriad of applications - from automated
vehicles and industrial automation to healthcare and service sectors, the
challenge of effective navigation has garnered significant attention. A
critical capability in these applications is the ability of the robot to
navigate reliably towards a specific target while avoiding obstacles in its
path, thereby ensuring efficient task completion and safety.

Figure 1: TurtleBot3 robot simulated in a simple environment with the Gazebo
framework.

Reinforcement Learning (RL) provides a robust framework for training agents to
learn from their interactions with the environment and make optimal decisions
[1]. With the advent of deep learning and the development of algorithms
capable of handling continuous and high-dimensional action spaces, RL has
demonstrated substantial success in various complex tasks, including robotic
navigation. Among these, Deep Q-Network (DQN) [2], Deep Deterministic Policy
Gradient (DDPG) [3], and Twin Delayed DDPG (TD3) [4], have been pivotal,
contributing significantly to advancements in the field.

These methods, however, primarily focus on optimizing a single objective
function. While this has proven effective in many scenarios, real-world
environments often present situations that require balancing multiple,
potentially conflicting objectives. For example, in the context of robotic
navigation, a robot might need to reach its destination quickly while also
minimizing energy consumption and collision risk. The need to address these
multifaceted objectives simultaneously has been recognized, leading to the
development of Multi-Objective Reinforcement Learning (MORL) [5].

MORL extends the traditional RL paradigm by optimizing a vector of rewards
instead of a single scalar reward, allowing an agent to balance multiple
objectives in a more sophisticated manner. While the body of work around MORL
has been growing, its applications to complex and dynamic robotic navigation
tasks remain relatively unexplored [6].

In this paper, we present a comparative analysis of single-objective and
multi-objective reinforcement learning strategies for training a robot to
navigate efficiently and effectively. Building upon the strengths of MORL, we
propose a novel application of this method in the context of robotic
navigation. We redefine the reward function to return a vector of rewards,
each reflecting a distinct objective, thereby enabling the robot to balance
the different goals and aim towards achieving a Pareto optimal solution.

The main contributions of our work are as follows:

  * •

A comprehensive comparative analysis between single-objective RL strategies
(DQN, DDPG, TD3) and a multi-objective reinforcement learning approach in the
context of robotic navigation tasks.

  * •

A novel application of MORL to complex and dynamic robotic navigation
scenarios that inherently involve multiple objectives, providing a method for
balancing and simultaneously optimizing these objectives.

  * •

An empirical demonstration of the effectiveness of our proposed MORL approach.
Through extensive simulations, we show that the proposed approach outperforms
traditional single-objective methods, offering more adaptable and robust
behaviors across varying levels of environmental complexity.

The remainder of this paper is organized as follows: Section 2 reviews related
work in the fields of reinforcement learning, single and multi-objective, with
an emphasis on their application to robotic navigation. Section 3 defines the
problem and presents the methodology employed in our study. The experimental
setup and results are detailed in Sections 4 and 5, respectively. The final
section, Section 6, concludes the paper and suggests future research
directions.

##  II Background Information

###  II-A Markov Decision Process (MDP)

In order to use RL to solve a problem, the problem must be structured as a
Markov Decision Process (MDP), which consists of five components: an agent,
environment, state, action, and reward. The agent is the learner that is to be
trained with RL, and it interacts with the environment by taking actions and
receiving rewards based on those actions. The environment is the real-world
environment with which the agent must operate in and execute actions. States
represents the situation of the environment at any point, with different
states capturing changes such as the agent’s position and direction as well as
the position of other objects around it. Actions represent what the agent can
do to interact with its environment, where actions may be a defined set of
possible action choices or an infinite set. The reward is the positive or
negative reinforcement the agent receives from executing a given action in a
given state of the environment, and it is the feedback that evaluates the
action quality.

In this study, the agent is the robot being trained with RL algorithms within
Gazebo environments which have various obstacles. The robot’s objective is to
navigate to the end goal as efficiently as possible while avoiding obstacles.
The robot executes actions by moving, with movements being represented through
a continuous, finite range of linear directions and velocities. The state can
be captured at any point, and it represents the current position of all
objects in the environment (the robot, the obstacles, and the end goal). A
penalty is given each time step the robot takes in order to encourage timely
navigation, and a large positive or negative reward is given based on the
episode outcome: the robot is rewarded for successfully reaching the goal and
penalized for failing to, such as by colliding with an obstacle or going
beyond the specified time limit.

An episode is every simulation run which ends with an outcome. Over a sequence
of time steps, the agent interacts with its environment through the MDP cycle,
doing so until the episode ends, or terminates. By repeatedly training the
agent over numerous episodes, the agent will by trial-and-error be able to
learn a policy, which is a strategy of deciding what action to take in a given
state. RL algorithms therefore aim to train the agent to learn the optimal
policy, since the agent will be able to determine what are the best actions to
take in any given situation, allowing it to efficiently achieve its objective
and maximize the total reward.

###  II-B Deep Q-Network (DQN)

While RL algorithms aim to learn the optimal policy, where the optimal action
for any given state is known, the optimal policy can be derived indirectly by
learning the values of all actions. The value of an action in a given state is
called the Q-value. By knowing the Q-values of all possible actions in a given
state, the optimal policy can be derived by taking the action with the highest
Q-value. Deep Q-Network (DQN) is an RL algorithm which maps a state to the
Q-values of all actions that can be taken from that state; it does so through
building a table with the Q-values of all actions that can be taken from every
state. DQN predicts the Q-values of actions by estimating the Q-function, a
process which occurs over training. As DQN is used to train the agent, the
feedback that is received is used to update the weights of the Q-function,
allowing it to predict the Q-values of actions with greater accuracy. DQN’s
architecture uses two neural networks, which are optimal function
approximators: a Q-network and the target network. The Q-network is the agent
that is trained to produce the optimal Q-values. The target network is
identical to the Q-network, in that it predicts target Q-values from the next
state. In addition, DQN has a component called Experience Replay, which
gathers training data by interacting with the environment, allowing this data
to be used to train the Q-network.

###  II-C Deep Deterministic Policy Gradient (DDPG)

Rather than learning the Q-values of actions, Deep Deterministic Policy
Gradient (DDPG) directly learns a policy. While DQN can only predict the
Q-values for discrete actions within a finite action setting, where actions
are discretely defined, DDPG is advantageous for being able to operate within
a continuous action setting [3]. For a given state, DDPG directly outputs the
best action to take from a continuous action space, hence the reason it is
called ”deterministic.” DDPG is an actor-critic algorithm, meaning it has two
neural networks, an actor and a critic. The actor is a policy network, where
it takes the state as input and outputs the exact action to take. The critic
is a Q-value network, which takes in the state and action and outputs the
Q-value. Essentially, the actor chooses what actions to take in a given state,
and the critic calculates the predicted Q-value of that action. In addition,
the critic and actor each have a target network, which help to stabilize
learning. The target networks are used to compute the target Q-values for the
next states. Using the feedback received from executing an action in the
state, the error between the predicted Q-value and the target Q-value is used
to update the actor and critic.

###  II-D Twin Delayed DDPG (TD3)

The successor to DDPG, Twin Delayed Deep Deterministic Policy Gradients (TD3)
is built to address an issue precedent RL algorithms have with continuously
overestimating the Q-values of the critic network [4]. A build up of these
estimation errors over time causes issues such as the agent falling into a
local optima, or extremes (minimum or maximum) of the objective function for a
given region of the input space. TD3 focuses instead on reducing
overestimation bias. This is done through the addition of three key features:
1) using two critic networks rather than one, 2) delaying updates of the actor
network, and 3) adding action noise regularization. For 1), TD3 uses two
critic networks to calculate the predicted Q-values and then takes the smaller
of the two Q-values, thus favoring underestimation of Q-values.
Underestimating values ensures overestimated values are not continuously
propagated through the algorithm. 2) deliberately ensures updates to the actor
(policy) network are less frequent than updates to the critic (value) network.
Training of the agent is badly affected by continuously overestimating poor
policies. Delaying updates to the actor allow the critic network to become
more stable and reduce errors in estimating Q-values before it is used to
update the actor. Lastly, 3) essentially adds noise to selected actions that
are fed to the target networks. By doing this, higher target Q-values will be
returned for actions that are more robust to noise and interference,
essentially giving higher value to more robust actions. This prevents having
high variance in the target Q-values that will be used to update the critic
networks.

##  III Related Work

Advancements in the field of reinforcement learning have led to various
strategies for optimal control in complex environments. The Q-learning
algorithm, as proposed by Watkins and Dayan [7], has become a cornerstone of
these strategies. Subsequent improvements led to the development of Deep
Q-Network (DQN) by Mnih et al. [2], where deep learning was combined with
reinforcement learning for improved performance on complex tasks. For handling
high-dimensional, continuous action spaces, Lillicrap et al. [3] introduced
the Deep Deterministic Policy Gradient (DDPG) algorithm, which uses deep
function approximators. Twin Delayed DDPG (TD3), an improvement over DDPG, was
proposed by Fujimoto et al. [4]. It mitigates the issue of overestimation bias
by using a pair of Q-functions and delaying policy updates.

However, in scenarios where tasks involve multiple, potentially conflicting
objectives, the limitations of the aforementioned single-objective
reinforcement learning techniques become evident. The concept of Multi-
Objective Reinforcement Learning (MORL) has been explored to tackle this
issue, as discussed by Roijers et al. [5]. The applications of MORL in robotic
navigation have been the subject of several studies. Van Moffaert and Nowé [6]
presented a scalarized MORL method optimizing energy consumption and travel
time, whereas Gomez et al. [8] employed MORL for optimizing safety and
efficiency in robot path planning. Xu et al. [9] proposed an efficient
evolutionary learning algorithm to find the Pareto set approximation for
continuous robot problems, demonstrating their approach produced denser and
higher-quality Pareto policies when balancing between robot movement and
energy efficiency. Lastly, Mannion et al. [10] discussed how making decisions
under outlined requirements for trustworthy AI is inherently a multi-objective
problem.

More recent works like those of Hwangbo et al. [11] and Abdolmaleki et al.
[12] demonstrate the potential of MORL for efficient robotic navigation and
control in real-world scenarios. Hwangbo et al. presented a learning framework
that combines model-free and model-based methods, demonstrating its
effectiveness in teaching a quadruped robot to run, climb, and jump over
obstacles. Meanwhile, Abdolmaleki et al. utilized MORL to optimize performance
in diverse simulated physical tasks.

Despite the aforementioned advancements, further exploration is required to
fully harness the potential of MORL in complex and dynamic robotic navigation
tasks, which is the focus of the current work.

##  IV Methods

###  IV-A System Architecture

In this study, the Gazebo simulation framework is used to simulate testing and
training the RL methods to control a TurtleBot3 robot simulation with Robot
Operating System (ROS). The TurtleBot3 robot reads data from LiDAR sensors
using ROS topics. Additionally, environment information about the end goal
location and obstacle locations are also published on ROS topics. ROS nodes
such as the DRL Environment Node and the DRL Agent Node subscribe to these
topics. ROS robot control features then allow the RL methods to use this
sensory and environment information in order to determine how the robot is to
navigate toward the end goal.

Figure 2: The system architecture describes how information about the
environment is input to the reward function, and its output helps the chosen
policy decide what action the robot should take.

###  IV-B Reward Function

In this study, the single-objective reward function takes in the robot’s
current action in the form of linear velocity (`action_linear`) and direction
of movement (`action_angular`), the robot’s distance and angle from the end
goal (`goal_dist` and `goal_angle`), its distance from the nearest obstacle
(`min_obstacle_dist`), and the status of the robot in terms of reaching the
end goal (`succeed`). The reward value is calculated based on the combination
of these terms, increasing or decreasing the cumulative result based on how
desirable or undesirable the input values are–thereby rewarding or penalizing
the behavior of the robot.

1\. Yaw Reward:

$r_{\text{yaw}}=-1\cdot\left|\text{goal\\_angle}\right|$

2\. Angular Velocity:

$r_{\text{vangular}}=-1\cdot\left(\text{action\\_angular}^{2}\right)$

3\. Distance Reward:

$r_{\text{distance}}=\frac{1}{2}\cdot\left(\frac{2\cdot\text{goal\\_dist\\_initial}}{\text{goal\\_dist\\_initial}+\text{goal\\_dist}}-1\right)$

4\. Obstacle Reward:

$r_{\text{obstacle}}=\begin{cases}-20&\text{if min\\_obstacle\\_dist}<0.22\\\
0&\text{otherwise}\end{cases}$

5\. Linear Velocity Reward:

$r_{\text{vlinear}}=-2\cdot\left(\left(0.22-\text{action\\_linear}\right)\times
10\right)^{2}$

6\. Total Reward without Success Condition:

$\text{reward}=r_{\text{yaw}}+r_{\text{distance}}+r_{\text{obstacle}}+r_{\text{vlinear}}+r_{\text{vangular}}-1$

7\. Final Reward with Success Condition:

$\text{reward}=\begin{cases}\text{reward}+2500,&\text{if succeed = SUCCESS}\\\
\hbox{\multirowsetup\text{reward} - 2000,}&\text{if succeed = }\\\
&\text{COLLISION\\_OBSTACLE or }\\\ &\text{COLLISION\\_WALL}\\\
\text{reward},&\text{otherwise}\end{cases}$

In step 6, after combining all the individually calculated values into the
final reward value, a penalty of `-1` is added to penalize the robot for each
time step that the reward function is called, indicating that the robot must
reach the end goal quickly in order to avoid the continuous decrease in
reward.

In step 7, a final large reward or penalty may be added to the cumulative
reward, depending on the robot’s success status: `+2500` is added if the robot
successfully reaches the goal (`succeed = SUCCESS`), and `-2000` is added if
the robot fails to, namely when it collides with an obstacle (`succeed =
COLLISION_OBSTACLE`) or with a wall (`succeed = COLLISION_WALL`). If neither
success condition applies, meaning the robot is still navigating to the end
goal without success or failure, then no additional value is added.

###  IV-C Environments

The RL methods are evaluated using the Gazebo simulation framework within two
different environments: stage A which features three non-moving obstacles and
four moving obstacles, and stage B, which features six moving obstacles. The
red circle indicates the end goal location which the robot aims to navigate
to.

Figure 3: Stage A has an environment with three non-moving obstacles (pink)
and four moving obstacles (white). Figure 4: Stage B has an environment with
six moving obstacles (white).

##  V Results

For training in both stages A and B, DDPG performed the best in terms of
having episodes yielding the highest rewards, followed by TD3, then DQN. All
three methods produced overall higher reward scores when training in stage B
than in stage A, as stage A is more difficult to navigate than stage B due to
having non-moving obstacles blocking key pathways of the environment. However,
all three algorithms still demonstrated consistent collisions with moving and
non-moving obstacles in both stages A and B.

For stage A, DDPG’s best performing episodes were demonstrated in later
episodes of training, with top four episodes `[3200, 3300, 3500, 3000]`
returning reward scores `[-2322, -2349, -2484, -2551]` respectively.

In contrast, TD3’s performance in stage A demonstrated better performing
episodes in earlier episodes of training, with top four episodes `[1200, 900,
600, 1000]` returning reward scores `[-2517, -2892, -2938, -2943]`
respectively. All of TD3’s best performing episodes had reward scores that
were lower than that of DDPG’s.

Lastly, the stage A performance of DQN showed better performing episodes in
the middle of training, with top four episodes `[2700, 2800, 3200, 2000]`
returning reward scores of `[-3579, -3766, -3782, -3823]` respectively. All of
DQN’s best performing episodes had the lowest reward scores of the three
algorithms in stage A.

Figure 5: The reward graph shows the performance of DQN, DDPG, and TD3 in
stage A over the course of 4000 episodes.

For stage B, DDPG again outperformed DQN and TD3, with best performing
episodes witnessed near the later part of training: the top four best
performing episodes `[2600, 2800, 2900, 3200]` returned reward scores `[-1126,
-1509, -1512, -1522]` respectively.

TD3’s performance in stage B also demonstrated the best performing episodes in
the later episodes of training, with top four episodes `[3600, 3500, 3300,
3400]` returning reward scores `[-1637, -1766, 2064, -2103]` respectively. The
top reward scores from TD3 in stage B were lower than the top reward scores
from DDPG.

Lastly, DQN’s best performing episodes from stage B were spread out from the
middle to the end of the training duration, with top four episodes `[2000,
3500, 3100, 3600]` returning reward scores `[-2401, -2532, -2578, -2629]`
respectively. DQN’s top performing episodes in stage B ultimately returned the
lowest reward scores of the three algorithms.

Figure 6: The reward graph shows the performance of DQN, DDPG, and TD3 in
stage B over the course of 4000 episodes.

Fig. 5 and Fig. 6 show the reward graphs for the performance of DQN, DDPG, and
TD3, in stage A and B respectively. These reward graphs demonstrate that DDPG
generally returned the highest rewards over the course of 4000 episodes. DQN
is noted to consistently return the lowest reward scores in both stages.
However, while DDPG generally outperformed TD3 in training, TD3 had a more
stable learning curve that may demonstrate greater improvement if testing
beyond 4000 episodes. This is likely due to TD3’s structure which reduces
overestimation of action values. Because DDPG’s structure does not prevent
against this overestimation bias, this may account for the volatility
witnessed in DDPG’s returned reward scores. For future implementation of MORL
algorithms, TD3 could be modified to accept a vector of reward values from an
MORL reward function.

While the individual performance of all algorithms were recorded, only DDPG’s
performance in stage A is shown as an example in Fig 7. For DDPG’s individual
performance in stage A, it is observed that the average actor loss gradually
decreased, meaning the actor generally learned the optimal robot actions,
though not at a consistent rate, as indicated by the spikes in the graph. On
the other hand, critic loss remained relatively high, indicating that the
critic failed to estimate the Q-values of these actions. Additionally, there
was a total of more collisions with walls and moving obstacles than there was
of successful episodes reaching the goal. The average reward also did not
increase significantly over episodes. Although only DDPG’s performance for
stage A is shown, the individual performances of the other algorithms in both
stages A and B demonstrated similar issues with actor and critic loss not
minimizing over the course of training, relatively low average rewards, and
consistently high collisions with walls and moving obstacles.

Figure 7: The performance of DDPG in stage A over 4000 episodes is shown.

Overall, DQN, DDPG, and TD3 did not perform optimally with consistent success
to reach the goal, particularly so when the end goal was located in more
difficult locations, such as the corners of the environment, where there were
more obstacles in the way. Frequent collisions with obstacles continued to be
observed even after constant training, and the robot did not reach the goal at
optimal speeds. While the algorithms gradually learned to reach the end goal
and therefore optimize the robot’s single objective with relative success,
they did so at the cost of neglecting other priorities, such as avoiding
obstacles and finding the shortest and quickest path to the end goal. These
algorithms consequently demonstrated overall poor performance, however these
results are to be expected, and in fact demonstrate the main issue: single-
objective RL strategies are unsuitable for providing a solution to the problem
of autonomous robotic navigation.

Improvements could be made with the implementation of an MORL algorithm and
reward function in order to optimize multiple objectives. DQN, DDPG, and TD3
demonstrated a failure to minimize collisions and reach the goal in a quicker
manner, and the single-objective reward function returned a reward score that
was composed of a large reward for successfully reaching the goal or a large
penalty for failing to. This large reward/penalty encouraged solely
prioritizing the objective at the cost of neglecting other priorities, as the
smaller rewards/penalties the robot gained for its movements, distance to the
goal, and distance to other obstacles were a minuscule component of the
cumulative reward. If the robot was able to reach the goal location despite
taking the longest path or navigating narrowly close to obstacles, it was
encouraged to do so by the large reward that it gained from successfully doing
so.

Consequently, the suitable solution is an MORL algorithm which uses a modified
reward function that returns a vector of rewards rather than a single reward.
This reward function allows a way to simultaneously balance and optimize
having multiple objectives, as each objective is represented by its own reward
score, therefore enabling the rewarding or penalizing of particular behaviors
in regards to each objective. This also allows modifying the reward function
and the modeled behavior of the robot according to preferences for particular
objectives. For example, a robot with safer navigation behavior can be
preferred by choosing behaviors which earn higher reward values for collision
avoidance.

##  VI Conclusion and Future Directions

We presented a comprehensive analysis of the performances between single-
objective RL strategies (DQN, DDPG, TD3) and highlighted the issues with using
such strategies for the problem of robot navigation, and examined how such
issues could be resolved through using multi-objective RL strategies instead.
We ran extensive training simulations with each single-objective algorithm and
examined the individual performances of each one in two different
environments, before comparing their performances through generated reward
graphs. We demonstrated how the performance of single-objective RL strategies
suffer as a result of being unable to provide an optimal solution which
balances and simultaneously optimizes multiple objectives.

In the future, we would like to test the trained models in order to fully
demonstrate the performance of each algorithm, in addition to finding ways to
increase the success rate of performance. Moreover, we plan to implement and
test multi-objective RL (MORL) strategies. Recent advancements in MORL have
shown promising results in resolving issues related to conflicting objectives,
a limitation inherent in single-objective RL strategies. By applying MORL, we
aim to provide a more holistic solution to the complex problem of robotic
navigation [13].

Additionally, we intend to extend our research to physical robots for further
testing and training. This step is crucial for understanding how these
algorithms perform in real-world scenarios, as opposed to controlled,
simulated environments [14]. By incorporating SynchroSim’s insights [15], we
plan to focus on synchronization and communication aspects of single-objective
and multi-objective RL strategies in multi-robot systems. The comparison of
single-objective and multi-objective RL strategies on different physical
robots will also provide valuable insights into the adaptability and
practicality of these approaches in varied robotic platforms.

We also plan to load the models onto a physical robot for further testing and
training, in addition to comparing the performance of single-objective and
multi-objective RL strategies on different physical robots as well.

## References

  * [1] R. S. Sutton and A. G. Barto.  Reinforcement learning an introduction– second edition.  MIT Press, 2020.
  * [2] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Andrei A Rusu, Joel Veness, Marc G Bellemare, Alex Graves, Martin Riedmiller, Andreas K Fidjeland, Georg Ostrovski, et al.  Human-level control through deep reinforcement learning.  Nature, 518(7540):529–533, 2015.
  * [3] Timothy P Lillicrap, Jonathan J Hunt, Alexander Pritzel, Nicolas Heess, Tom Erez, Yuval Tassa, David Silver, and Daan Wierstra.  Continuous control with deep reinforcement learning.  arXiv preprint arXiv:1509.02971, 2015.
  * [4] Scott Fujimoto, Herke Hoof, and David Meger.  Addressing function approximation error in actor-critic methods.  In International conference on machine learning, pages 1587–1596. PMLR, 2018.
  * [5] Diederik M Roijers, Peter Vamplew, Shimon Whiteson, and Richard Dazeley.  A survey of multi-objective sequential decision-making.  Journal of Artificial Intelligence Research, 48:67–113, 2013.
  * [6] Kristof Van Moffaert and Ann Nowé.  Multi-objective reinforcement learning using sets of pareto dominating policies.  Journal of Machine Learning Research, 15(1):3483–3512, 2014.
  * [7] Christopher JCH Watkins and Peter Dayan.  Q-learning.  Machine learning, 8(3-4):279–292, 1992.
  * [8] Faustino Gomez, Jürgen Schmidhuber, and Risto Miikkulainen.  Efficient non-linear control through neuroevolution.  European Journal of Operational Research, 185(2):593–603, 2008\.
  * [9] Jie Xu, Yunsheng Tian, Pingchuan Ma, Daniela Rus, Shinjiro Sueda, and Wojciech Matusik.  Prediction-guided multi-objective reinforcement learning for continuous robot control.  In Proceedings of the 37th International Conference on Machine Learning, 2020.
  * [10] Patrick Mannion, Fredrik Heintz, Thommen George Karimpanal, and Peter Vamplew.  Multi-objective decision making for trustworthy ai.  In Proc. of the 1st Multi-Objective Decision Making Workshop, 2021\.
  * [11] Jemin Hwangbo, Joonho Lee, Alexey Dosovitskiy, Dario Bellicoso, Vassilios Tsounis, Vladlen Koltun, and Marco Hutter.  Learning agile and dynamic motor skills for legged robots.  Science Robotics, 4(26), 2019.
  * [12] Abbas Abdolmaleki, Jost Tobias Springenberg, Yuval Tassa, Remi Munos, Nicolas Heess, and Martin Riedmiller.  Maximum a posteriori policy optimisation.  In International Conference on Learning Representations, 2018.
  * [13] Jumman Hossain, Abu-Zaher Faridee, Nirmalya Roy, Anjan Basak, and Derrik E. Asher.  Covernav: Cover following navigation planning in unstructured outdoor environment with deep reinforcement learning.  In 2023 IEEE International Conference on Autonomic Computing and Self-Organizing Systems (ACSOS), pages 127–132, 2023.
  * [14] Avi Spector, Wanying Zhu, Jumman Hossain, and Nirmalya Roy.  Simulated forest environment and robot control framework for integration with cover detection algorithms.  In 2022 IEEE/ACM International Conference on Big Data Computing, Applications and Technologies (BDCAT), pages 277–283. IEEE, 2022.
  * [15] Emon Dey, Jumman Hossain, Nirmalya Roy, and Carl Busart.  Synchrosim: An integrated co-simulation middleware for heterogeneous multi-robot system.  In 2022 18th International Conference on Distributed Computing in Sensor Systems (DCOSS), pages 334–341. IEEE, 2022.


# Ferromagnetic Heisenberg model with the Dzyaloshinskii-Moriya interaction

E. Albayrak  $|$ Corresponding author<EMAIL_ADDRESS>
(Received May 30, 2022, in final form July 25, 2022)

###### Abstract

The spin-1/2 Heisenberg model is formulated in terms of a mean-field
approximation (MFA) by using the matrix forms of spin operators
$\hat{S}_{x},\hat{S}_{y}$ and $\hat{S}_{z}$ in three-dimensions. The
considered Hamiltonian consists of bilinear exchange interaction parameters
$(J_{x},J_{y},J_{z})$, Dzyaloshinskii-Moriya interactions
$(\Delta_{x},\Delta_{y},\Delta_{z})$ and external magnetic field components
$(H_{x},H_{y},H_{z})$. The magnetization and its components are obtained in
the MFA with the general anisotropic case with $J_{x}\neq J_{y}\neq J_{z}$ for
various values of coordination numbers $q$. Then, the thermal variations of
magnetizations are investigated in detail to obtain the phase diagrams of the
model for the isotropic case with $J_{x}=J_{y}=J_{z}>0$. It is found that the
model exhibits ferromagnetic, paramagnetic, random phase regions and an extra
ferromagnetic phase at which the components of magnetizations present
branching.

Key words: spin 1/2, ferromagnetic, Dzyaloshinskii-Moriya, XYZ model,
magnetization, phase diagrams

###### Abstract

Ñïí-1/2 ìîäåëü Ãàéçåíáåðãà ñôîðìóëüîâàíî â íàáëèæåíí ñåðåäíüîãî ïîëÿ ç
âèêîðèñòàííÿì ìàòðèíî¿ ôîðìè ñïíîâèõ îïåðàòîðâ $\hat{S}_{x},\hat{S}_{y}$
$\hat{S}_{z}$ ó òðèâèìðíîìó ïðîñòîð. Ó ãàìëüòîíàí, ùî ðîçãëÿäàòüñÿ, âõîäÿòü
ïàðàìåòðè áëéíéíî¿ îáìííî¿ âçàìîä¿ $(J_{x},J_{y},J_{z})$, âçàìîä¿
Äçÿëîøèíñüêîãî-Ìîðÿ $(\Delta_{x},\Delta_{y},\Delta_{z})$, à òàêîæ êîìïîíåíòè
çîâíøíüîãî ìàãíòíîãî ïîëÿ $(H_{x},H_{y},H_{z})$. Êîìïîíåíòè íàìàãíåíîñò
îòðèìàí ó íàáëèæåíí ñåðåäíüîãî ïîëÿ äëÿ çàãàëüíîãî àíçîòðîïíîãî âèïàäêó, êîëè
$J_{x}\neq J_{y}\neq J_{z}$, äëÿ ðçíèõ êîîðäèíàöéíèõ èñåë $q$. Òåìïåðàòóðí
çàëåæíîñò íàìàãíåíîñòåé äîñëäæåíî äåòàëüíî ç ìåòîþ ïîáóäîâè ôàçîâèõ äàãðàì
ìîäåë äëÿ çîòðîïíîãî âèïàäêó $J_{x}=J_{y}=J_{z}>0$. Âèÿâëåíî, ùî äëÿ ìîäåë
ñíóþòü ôåðîìàãíòíà, ïàðàìàãíòíà, âèïàäêîâà ôàçîâ îáëàñò, à òàêîæ äîäàòêîâà
ôåðîìàãíòíà ôàçà, â ÿêé ó êîìïîíåíò íàìàãíåíîñòåé ñïîñòåðãàòüñÿ ãàëóæåííÿ.

Ключов слова: ñïí 1/2, ôåðîìàãíåòèêè, âçàìîäÿ Äçÿëîøèíñüêîãî-Ìîðÿ, $XYZ$
ìîäåëü, íàìàãíåíñòü, ôàçîâ äàãðàìè

##  1 Introduction

The Dzyaloshinskii-Moriya (DM) interaction [1, 2] is an antisymmetric exchange
interaction with a contribution to total magnetic exchange interaction between
two neighboring magnetic spins as a source of weak ferromagnetic (FM) behavior
in an antiferromagnet [3]. The helical structure is induced by the DM
interaction which breaks the chiral symmetry, and thus, the two same helical
structures with different winding directions do not degenerate [4]. It is
revealed theoretically that a chiral magnetic soliton lattice is formed with a
finite magnetic field perpendicular to the axis of the helical structure, and
a continuous phase transition to a forced ferromagnetic phase occurs with
increasing magnetic field [5]. These topologically protected magnetic
structures are called skyrmions and their existence was proven experimentally
[6, 7, 8, 9, 10, 11, 12, 13]. Due to its potential applications in spintronics
devices, they have drawn much attention. Note also that the competition
between DM and exchange interactions is important for skyrmion-lattice phases.
Even though the DM interaction has been well studied in crystalline magnets,
there are also numerous experimental [14, 15] and theoretical [16, 17] studies
demonstrating its importance in spin glasses.

The DM interaction may be rather difficult to deal with analytically and thus
it may require numerical investigations. Because of the nature of DM
interaction, it is necessary to study it with the spins at least in two-
dimensions in contrast to the Ising spins. The quantum phase transitions (QPT)
in a bond-alternative antiferromagnetic (AFM) Ising chain were considered
including the DM interaction [18] by using the fidelity based on the infinite
matrix product states algorithm. It was found that antiferromagnetic and
disordered phases exist in the ground state and the transition between them is
continuous. The low-energy excitation spectrum and the ground-state magnetic
phase diagram of the spin-1/2 ferromagnetic Ising chain with the DM
interaction were considered where a first-order metamagnetic phase transition
between FM and a spiral phases were reported [19].

The spin-1/2 models with the inclusion of the DM interaction (DMI) were also
established in some works. The Ising–Heisenberg model on the triangulated
kagome lattice was exactly solved by establishing a precise mapping
correspondence to the simple spin-1/2 Ising model [20]. The quantum
anisotropic Heisenberg antiferromagnet model in the presence of a DMI and a
uniform longitudinal magnetic field was considered by using the effective-
field theory with a finite cluster $N=2$ spin [21]. The QPTs and ground-state
phase diagram of the HeisenbergIsing alternating chain with uniform DM
interaction were investigated by a matrix-product-state method [22]. The pair
XYZ Heisenberg interaction and quartic Ising interactions were exactly solved
by establishing a precise mapping relationship with the corresponding zero-
field eight-vertex model [23]. The anisotropic Heisenberg model with
antiferromagnetic exchange interactions in the presence of a longitudinal
external magnetic field and a DMI was studied by employing the usual mean-
field approximation [24]. The ground-state magnetic phase diagram of an
antiferromagnetic two-leg ladder with period two lattice units modulated DMI
along the legs was considered [25]. The results of a combined analytical and
density matrix renormalized group study of the AFM XXZ Heisenberg chain
subject to a uniform DMI and a transverse magnetic field were reported [26].
The linked cluster expansion techniques was applied to study the polarized
high-field phase of an antiferromagnet on the kagome lattice with Heisenberg
and DMI [27]. The phase diagram in the $H$-$T$ plane of the potassium jarosite
compound KFe3(OH)6(SO4)2 for the antiferromagnetic XY model with DMI was
investigated by using the mean-field theory for different values of DM [28].
On the two-dimensional non-linear $\Sigma$-model describing a ferromagnet with
DMI, three families of exact static solutions depending on a single Cartesian
variable were obtained [29]. The roles of spatial anisotropy, DM interactions
and quantum fluctuations on the magnetization process of a triangular
antiferromagnet were considered [30]. The ground-state phase diagram of a one-
dimensional XXZ chain with a spatially modulated DMI in the presence of an
alternating magnetic field was calculated [31]. The steady-state phase diagram
of the dissipative XYZ model on a two-dimensional triangular lattice by means
of cluster mean-field approximation was constructed [32]. The last work that
we can mention is the study of the quantum spin liquid material
herbertsmithite described by an AFM Heisenberg model on the kagome lattice
[33].

In addition to the above theoretical works, we can also mention some
experimental works such as: magnetization of triangular-lattice
antiferromagnets Ba3CoSb2O9 and CsCuCl3 with a $120^{\rm o}$ spin structure in
the $ab$ plane [34], the nature of possible magnetic phases in the frustrated
hyperkagome iridate, Na4Ir3O8, based on the Kitaev-Heisenberg model with DM
interactions [35]; the results of magnetic measurements performed on
geometrically frustrated Ni3V2O8 and Ni3(1-x)Co3xV2O8 single crystals with
$x=0.03$ were presented [36].

In this work, the magnetization components $(M_{x},M_{y},M_{z})$ and thus
magnetization ($M_{T}$) are obtained in terms of the MFA for the general
values of the $(J_{x},J_{y},J_{z})$, $(\Delta_{x},\Delta_{y},\Delta_{z})$ and
$(H_{x},H_{y},H_{z})$ and $q$. Then, thermal variations of magnetizations are
studied for ferromagnetic interactions between the nearest-neighbor spins,
i.e., $J=J_{x}=J_{y}=J_{z}>0$, with equal DM interactions and external
magnetic field components. In addition, the phase diagrams are obtained on the
planes of ($\Delta_{m}/J,T/J$) for $q=3,4$ and 6 when $H$ is set equal to
zero. Even though our equations are obtained for the anisotropic case, the
results are only presented for the isotropic case for simplicity. The reason
of this study is that our literature search did not even reveal any results
for the three-dimensional isotropic case.

The rest of the work is set up as follows: the formulation for the MFA in
terms of the spin operators is presented in section 2, thermal variations of
magnetizations are demonstrated in section 3, phase diagrams are illustrated
in section 4 and the last section includes a brief summary and conclusions.

##  2 The formulation

The anisotropic spin-1/2 XYZ Heisenberg Hamiltonian in terms of the bilinear
exchange interaction parameters and the DM interactions between the nearest-
neighbor (NN) spins and the external magnetic field components acting on each
spin site is given as

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle-J_{x}\sum_{\langle
i,j\rangle}S_{i}^{x}S_{j}^{x}-J_{y}\sum_{\langle
i,j\rangle}S_{i}^{y}S_{j}^{y}-J_{z}\sum_{\langle
i,j\rangle}S_{i}^{z}S_{j}^{z}-\Delta_{x}\sum_{\langle
i,j\rangle}(S_{i}^{y}S_{j}^{z}-S_{i}^{z}S_{j}^{y})-\Delta_{y}\sum_{\langle
i,j\rangle}(S_{i}^{z}S_{j}^{x}-S_{i}^{x}S_{j}^{z})$ (2.1) $\displaystyle-$
$\displaystyle\Delta_{z}\sum_{\langle
i,j\rangle}(S_{i}^{x}S_{j}^{y}-S_{i}^{y}S_{j}^{x})-H_{x}\sum_{i}S_{i}^{x}-H_{y}\sum_{i}S_{i}^{y}-H_{z}\sum_{i}S_{i}^{z},$

where $\langle i,j\rangle$ refers to the summation over the NN spins.
$\hat{S}_{i}^{x},\hat{S}_{i}^{y}$ and $\hat{S}_{i}^{z}$ are the components of
spin-1/2 operator at site $i$ which are given in the matrix forms as

$\displaystyle\hat{S}_{i}^{x}=\frac{1}{2}\left(\begin{array}[]{cc}0&1\\\
1&0\\\ \end{array}\right),\hskip
1.42271pt\hat{S}_{i}^{y}=\frac{1}{2}\left(\begin{array}[]{cc}0&-i\\\ i&0\\\
\end{array}\right)\hskip 1.42271pt,\hskip
1.42271pt\hat{S}_{i}^{z}=\frac{1}{2}\left(\begin{array}[]{cc}1&0\\\ 0&-1\\\
\end{array}\right).$ (2.8)

In the MFA, the Hamiltonian $\mathcal{H}$ in equation (2.1) can be written in
the MF form as

$\displaystyle-\beta\mathcal{H}_{\rm{MFA}}=-\beta\sum_{i}\mathcal{H}_{\rm{MFA}}^{(i)},$
(2.9)

in which

$\displaystyle-\beta\mathcal{H}_{\rm{MFA}}^{(i)}$ $\displaystyle=$
$\displaystyle\beta
q(J_{x}M_{x}S_{i}^{x}+J_{y}M_{y}S_{i}^{y}+J_{z}M_{z}S_{i}^{z})+\beta(H_{x}S_{i}^{x}+H_{y}S_{i}^{y}+H_{z}S_{i}^{z})$
(2.10) $\displaystyle+$
$\displaystyle\Delta_{x}(S_{i}^{y}M_{z}-S_{i}^{z}M_{y})+\beta\Delta_{y}(S_{i}^{z}M_{x}-S_{i}^{x}M_{z})+\beta\Delta_{z}(S_{i}^{x}M_{y}-S_{i}^{y}M_{x}),$

where $M_{\mu}=\langle S_{j}^{\mu}\rangle$ are the magnetization components
with ${\mu}=x,y,z$ and $\beta=1/({\rm k}T)$ with ${\rm k}$ being the Boltzmann
constant set equal to 1 for convenience.

The matrix representation of $-\beta\mathcal{H}_{\rm{MFA}}^{(i)}$ is obtained
by using the spin operators, i.e., equation (2.8), and is found as

$\displaystyle-\beta\mathcal{H}_{\rm{MFA}}^{(i)}=\left(\begin{array}[]{cc}H_{11}&H_{12}\\\
H_{21}&H_{22}\\\ \end{array}\right),$ (2.13)

where the matix elements are given as

$\displaystyle H_{11}$ $\displaystyle=$
$\displaystyle{\beta}/{2}[q(\Delta_{y}M_{x}-\Delta_{x}M_{y}+J_{z}M_{z})+H_{z}],$
$\displaystyle H_{12}$ $\displaystyle=$
$\displaystyle{\beta}/{2}[q(J_{x}M_{x}+\Delta_{z}(\mathrm{i}M_{x}+M_{y})-\mathrm{i}J_{y}M_{y}-\mathrm{i}\Delta_{x}M_{z}-\Delta_{y}M_{z})+(H_{x}-\mathrm{i}H_{y})],$
$\displaystyle H_{21}$ $\displaystyle=$
$\displaystyle{\beta}/{2}[q(J_{x}M_{x}+\Delta_{z}(-\mathrm{i}M_{x}+M_{y})+\mathrm{i}J_{y}M_{y}+\mathrm{i}\Delta_{x}M_{z}-\Delta_{y}M_{z})+(H_{x}+\mathrm{i}H_{y}),$
$\displaystyle H_{22}$ $\displaystyle=$
$\displaystyle{\beta}/{2}[q(-\Delta_{y}M_{x}+\Delta_{x}M_{y}-J_{z}M_{z})-H_{z}].$

It is clear that $H_{12}=H_{21}^{*}$, which ensures that the eigenvalues of
this matrix are real. Thus, the eigenvalues are given as

$\displaystyle\Lambda_{1,2}$ $\displaystyle=$
$\displaystyle\pm{\beta}/{2}\left[\left(H_{x}^{2}+H_{y}^{2}+H_{z}^{2}\right)+2q\left(-\Delta_{z}H_{y}M_{x}+\Delta_{y}H_{z}M_{x}+H_{x}J_{x}M_{x}+\Delta_{z}H_{x}M_{y}\right.\right.$
$\displaystyle-$
$\displaystyle\left.\Delta_{x}H_{z}M_{y}+H_{y}J_{y}M_{y}-\Delta_{y}H_{x}M_{z}+\Delta_{x}H_{y}M_{z}+H_{z}J_{z}M_{z}\right)+q^{2}\left(\Delta_{y}^{2}M_{x}^{2}+\Delta_{z}^{2}M_{x}^{2}\right.$
$\displaystyle+$ $\displaystyle
J_{x}^{2}M_{x}^{2}-2\Delta_{x}\Delta_{y}M_{x}M_{y}+2\Delta_{z}J_{x}M_{x}M_{y}-2\Delta_{z}J_{y}M_{x}M_{y}+\Delta_{x}^{2}M_{y}^{2}+\Delta_{z}^{2}M_{y}^{2}+J_{y}^{2}M_{y}^{2}$
$\displaystyle-$ $\displaystyle
2\Delta_{x}\Delta_{z}M_{x}M_{z}-2\Delta_{y}J_{x}M_{x}M_{z}+2\Delta_{y}J_{z}M_{x}M_{z}-2\Delta_{y}\Delta_{z}M_{y}M_{z}+2\Delta_{x}J_{y}M_{y}M_{z}$
$\displaystyle-$
$\displaystyle\left.\left.2\Delta_{x}J_{z}M_{y}M_{z}+\Delta_{x}^{2}M_{z}^{2}+\Delta_{y}^{2}M_{z}^{2}+J_{z}^{2}M_{z}^{2}\right)\right]^{1/2}.$

Having obtained the eigenvalues ($\Lambda_{1},\Lambda_{2}$), we are now ready
to obtain the partition function which is given as

$\displaystyle Z_{i}$ $\displaystyle=$
$\displaystyle\mathrm{(}Tr)_{(i)}\exp[-\beta\mathcal{H}_{\rm{MFA}}^{(i)}]$
(2.14) $\displaystyle=$
$\displaystyle\sum_{n=1}^{2}\mathrm{e}^{\Lambda_{n}}=\mathrm{e}^{\Lambda_{1}}+\mathrm{e}^{\Lambda_{2}}=2\cosh(\Lambda).$

since $\Lambda=\Lambda_{1}=-\Lambda_{2}$. The free energy of the model is
found from the well-known definition by using the partition function as

$\displaystyle f=-\frac{1}{\beta}\ln Z_{i}$ (2.15)

which will be utilized to find the formulations of the order-parameters. The
dipolar moments or magnetization components, $M_{\mu}=\langle
S_{j}^{\mu}\rangle$ with ${\mu}=x,y,z$ are found from

$\displaystyle M_{\mu}$ $\displaystyle=$ $\displaystyle\langle
S_{i}^{\mu}\rangle=-\frac{\partial f}{\partial
H_{\mu}}=\frac{1}{\beta}\frac{\partial\ln Z_{i}}{\partial
H_{\mu}}=\frac{\mathrm{(}Tr)_{(i)}\left\\{S_{i}^{\mu}\exp\left[-\beta\mathcal{H}_{\rm{MFA}}^{(i)}\right]\right\\}}{Z_{i}}$
(2.16) $\displaystyle=$
$\displaystyle\frac{1}{\beta}\left[\frac{({\partial\Lambda_{1}}/{\partial
H_{\mu}})\mathrm{e}^{\Lambda_{1}}+({\partial\Lambda_{2}}/{\partial
H_{\mu}})\mathrm{e}^{\Lambda_{2}}}{\mathrm{e}^{\Lambda_{1}}+\mathrm{e}^{\Lambda_{2}}}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\beta^{\prime}}\left[\frac{({\partial\Lambda_{1}}/{\partial
h_{\mu}})\mathrm{e}^{\Lambda_{1}}+({\partial\Lambda_{2}}/{\partial
h_{\mu}})\mathrm{e}^{\Lambda_{2}}}{\mathrm{e}^{\Lambda_{1}}+\mathrm{e}^{\Lambda_{2}}}\right],$

where $\beta^{\prime}=\beta J_{z}$ and $h_{\mu}=H_{\mu}/J_{z}$ which are the
reduced temperature and external magnetic field components. The explicit forms
of the magnetization components are too long to be given explicitly. The
magnitude of the magnetization vector, $\vec{M}_{T}$=
$M_{x}\hat{i}+M_{y}\hat{j}+M_{z}\hat{k}$, is obtained by using the
magnetization components given in equation (2.16) as

$\displaystyle
M_{T}=\sqrt{\sum_{\mu=x,y,z}M_{\mu}^{2}}=\sqrt{M_{x}^{2}+M_{y}^{2}+M_{z}^{2}}.$
(2.17)

Having obtained the formulations for the magnetization components and
magnetization in the MFA for the general case, we are now ready to study their
thermal variations for the given values of $\Delta_{\mu}$, $J_{\mu}$,
$H_{\mu}$ and the coordination number $q=3,4$ and 6 corresponding to the
honeycomb, square and simple cubic lattices. An iterative procedure is
followed for the calculation of our order-parameters with the given values of
the system parameters under temperature variations. The obtained results are
only presented for the FM case with $J=J_{x}=J_{y}=J_{z}>0$ corresponding to
the isotropic case, as seen in the next section.

Figure 1: (Colour online) The thermal variations of magnetization components
and magnetization when $H=0$ and $q=6$ for the given values of $\Delta_{m}$ as
(a) $\Delta_{x}=\Delta_{y}=\Delta_{z}=0$ and $J_{x}=J_{y}=J_{z}$ case with
$H=0$ and in the inset $H=1$, (b) 1.3, (c) 1.2, (d) 1, (e) 0.94, (f) 0.9 and
(g) 0.6.

##  3 Thermal variations of magnetizations

In this section, the characteristic thermal variations of magnetizations are
illustrated when the external magnetic field $H$ is turned off and on. They
are presented for $q=6$ only, because of the qualitative similarities with the
$q=3$ and 4. There is only quantitative difference between them, i.e., the
critical temperatures are observed at higher values for higher $q$ which is
expected.

First, we perform a numerical calculation test for the simplified case of
$\Delta_{x}=\Delta_{y}=\Delta_{z}=0$ and $J_{x}=J_{y}=J_{z}$ when the
Hamiltonian (2.1) reduces to the isotropic Heisenberg model in which the MFA
method is well-established and can be used as reference. Figure 1 a shows the
case with zero external magnetic field with all the magnetizations tending to
zero at the same $T_{c}$ and the inset is obtained for $H=1$ shows that
magnetizations do not tend to zero. These are the well-known results of the
isotropic Heisenberg model. It should be noted that the behavior of
magnetizations changes for the values of
$\Delta_{m}=\Delta_{x}=\Delta_{y}=\Delta_{z}$ being large (I), greater than 1
but close to it (II), equal to one (III) and less than 1 (IV) in some
characteristic forms for $H=0$ as we shall see later. Figure 1 b is obtained
for the case (I) with $\Delta_{m}=1.3$ and shows that magnetizations are less
than 0.5 at zero temperature, and they decrease as the temperature decreases
and eventually terminate at the second-order phase transition temperature
($T_{c}$) separating the FM and PM phases. It is also clear that
$M_{T}>M_{x}=M_{z}=-M_{y}$. Case (II) is presented in figure 1 c for
$\Delta_{m}=1.2$ and shows that magnetizations first exhibit some random
behaviors moving up and down which terminates at the critical temperature
called $T_{R}$, then their behaviors become similar to the case (I) with all
the lines terminating at the $T_{c}$. When $\Delta_{m}$ is set equal to 1 for
the case (III) as shown in figure 1 d, the magnetizations show only random
behaviors from the beginning to the end. Now, $T_{R}$ separates the random
phase region from the PM phase. Figure 1 e shows some branching of
magnetizations calculated for $\Delta_{m}=0.94$. The lines start with random
behaviors in the interval of some values corresponding to
$M_{T}>M_{x}=M_{z}>M_{y}$ with branching. The random behavior terminates at
the $T_{R}$, then the lines become regular curves but branching continues
which terminates at a temperature called $T_{B}$. At $T_{B}$, we see that
$M_{x}=M_{y}=M_{z}$. Then, they follow each other terminating at the $T_{c}$
as $M_{T}$. The branching is still seen for $\Delta_{m}=0.9$ as shown in
figure 1 f with no random behaviors anymore. Again, the branching terminates
at $T_{B}$, then the rest is as in figure 1 e. Finally, for $\Delta_{m}=0.6$,
all the components of magnetizations follow the same curve, i.e.,
$M_{T}>M_{x}=M_{y}=M_{z}$ and terminates at the $T_{c}$. It is interesting to
note that when the random behavior or branching appears, $M_{T}$ may be a
little larger than 0.5, otherwise it is either smaller or equal to 0.5. It may
be interesting to further search the reason of these fluctuations, and it is
obvious that it has quantum mechanical origin. It should also be mentioned
that these random behaviors in magnetizations must be caused by the existence
of the DM interactions leading to the skyrmions which will not be examined in
this work.

Thermal behaviors of magnetizations when $H$ is turned on shows two
characteristic forms. When $\Delta_{m}>1$ and $H=1$, we see that branching
continues as $M_{x}=M_{z}>M_{y}$ which terminates at $T_{B}$ where
$M_{x}=M_{y}=M_{z}$ (see figure 2 a). Afterwards, they follow each other,
never tending to the zero which is expected when $H$ is on. It is also clear
that $M_{T}$ is a little larger than 0.5. When $\Delta_{m}\leqslant 1$ and
$H=1$, we always see $M_{T}>M_{x}=M_{y}=M_{z}$ as shown in figure 2 b. Again,
they do not tend to zero with increasing temperature.

In the next section, we combine all these critical temperatures to construct
the phase diagrams on the ($\Delta_{m}/J,T/J$) planes when $H$ is turned off.
The combinations of these points make up the possible phase lines and their
combinations lead to some critical points.

Figure 2: (Colour online) Thermal variations of magnetization components and
magnetization when $H=1$ and $q=6$ for the given values of $\Delta_{m}$ as (a)
1.125 and (b) 0.75.

##  4 The phase diagrams

Now, we construct the phase diagrams on the $(\Delta_{m}/J,T/J)$ planes for
the given values of the coordination numbers $q=3,4$ and 6. In the phase
diagrams we identify five different phase regions:

  * •

The FM phase with constant $T_{c}$-line corresponding to the values of
$\Delta_{m}\leqslant 1$ where $M_{x}=M_{y}=M_{z}$ as seen in figure 1 f (see
the part after $T_{B}$) and in figure 1 g.

  * •

The FM phase with linearly increasing $T_{c}$-line corresponding to the values
of $\Delta_{m}>1$ where $M_{x}=M_{z}=-M_{y}$ as seen in figure 1 b and c after
$T_{R}$ with the slopes of about 1.5, 2.0 and 3.0 for $q=3,4$ and 6,
respectively. It should be noted that these values are also the same values in
the Ising model where the FM phase finishes and the PM phase starts, i.e.,
order-disorder phase transition temperatures, for $q=3,4$ and 6, respectively.

  * •

The random phase regions indicated with (R) where magnetizations behave
randomly going up and down as seen in figure 1 c–e before $T_{R}$. It is also
interseting to see the branching of magnetization components in this phase.

  * •

The FM phase region with branching indicated with (B) corresponding to
$M_{x}=M_{z}\neq M_{y}$ as indicated in figure 1 e and f.

  * •

The PM phase region with zero magnetizations caused by the thermal agitations.
The border lines between these phase regions, i.e., phase transition lines,
are indicated with solid, dotted-dashed and dashed lines indicating the
$T_{c},T_{R}$ and $T_{B}$-lines, respectively.

As seen in figure 3, the phase diagrams are similar for all values of the
coordination numbers. The $T_{c}$-lines consist of a straight part and a
linearly increasing part. The straight parts are found at temperatures
$T_{c}=0.7675,1.02$ and 1.52 for $q=3,4$ and 6, respectively. The $T_{R}$\-
and $T_{B}$-lines originate from 1 and are located around it. The two portions
of the $T_{R}$-lines starting from 1 terminate at zero temperature for about
$\Delta_{m}=0.97$ and $\Delta_{m}=1.27$ making a closed loop enclosing the R
phase for all $q$, but with higher temperatures for higher $q$. The
$T_{B}$-lines terminate at the same value of $\Delta_{m}$ being about 0.83 for
all $q$. The $T_{B}$-lines separate FM and R phases when $\Delta_{m}<1$. The
$T_{R}$-line when $\Delta_{m}>1$ separates the FM phase from the R phase. It
is clear that the model does not produce any first-order phase transition
lines as in the well-known Ising model. It should also be noted that the
$\Delta_{m}=1$ is a special critical point from where two $T_{c}$-, two
$T_{R}$\- and one $T_{B}$-lines emerge.

Figure 3: The phase diagrams on the $(\Delta_{m}/J,T/J)$ planes for the
isotropic FM phase when $H$ is turned off for (a) $q=3.0$, (b) $q=4.0$ and (c)
$q=6.0$.

##  5 Summary and conclusions

The spin-1/2 Heisenberg Hamiltonian for the FM case with
$J=J_{x}=J_{y}=J_{z}>0$ is considered to study the effects of DM interactions
with $\Delta_{m}=\Delta_{x}=\Delta_{y}=\Delta_{z}$ on the thermal variations
of the magnetization components and magnetization when the external magnetic
field components $H=H_{x}=H_{y}=H_{z}$ are turned off and on. The phase
diagrams of the model are calculated on the ($\Delta_{m}/J,T/J$) planes for
$q=3,4$ and 6. Three different FM phase regions are observed with
$M_{x}=M_{y}=M_{z}$, $M_{x}=M_{z}=-M_{y}$ and the one with $M_{x}=M_{z}\neq
M_{y}$ exhibits branching in addition to the PM and random phase regions. It
should be noted that this model has not been studied before, so the comparison
is not possible. As a last word, the anisotropic case of this model is going
to be considered as a continuation of this work.

## Acknowledgements

This work was supported by the Research Fund of Erciyes University with
Project Identification Number: FBA-2021-11571.

## References

  * [1] Dzyaloshinsky I., J. Phys. Chem. Solids, 1958, 4, 241, doi:10.1016/0022-3697(58)90076-3.
  * [2] Moriya T., Phys. Rev. Lett., 1960, 4, 228, doi:10.1103/PhysRevLett.4.228.
  * [3] Komatsu H., Nonomura Y., Nishino M., Phys. Rev. B, 2021, 103, 214404, doi:10.1103/PhysRevB.103.214404.
  * [4] Nishikawa Y., Hukushima K., Phys. Rev. B, 2016, 94, 064428, doi:10.1103/PhysRevB.94.064428.
  * [5] Shinozaki M., Hoshino S., Masaki Y., Kishine J., Kato Y., J. Phys. Soc. Jpn., 2016, 85, 074710,
doi:10.7566/JPSJ.85.074710.

  * [6] Mühlbauer S., Binz B., Jonietz F., Pfleiderer C., Rosch A., Neubauer A., Georgii R., Böni P., Science, 2009, 323, 915, doi:10.1126/science.1166767.
  * [7] Yu X. Z., Onose Y., Kanazawa N., Park J. H., Han J. H., Matsui Y., Nagaosa N., Tokura Y., Nature, 2010, 465, 901, doi:10.1038/nature09124.
  * [8] Jonietz F., Mühlbauer S., Pfleiderer C., Neubauer A., Münzer W., Bauer A., Adams T., Georgii R., Böni P., Duine R. A., Everschor K., Garst M., Rosch A., Science, 2010, 330, 1648, doi:10.1126/science.1195709.
  * [9] Yu X. Z., Kanazawa N., Onose Y., Kimoto K., Zhang W. Z., Ishiwata S., Matsui Y., Tokura Y., Nat. Mater., 2011, 10, 106, doi:10.1038/nmat2916.
  * [10] Schulz T., Ritz R., Bauer A., Halder M., Wagner M., Franz C., Pfleiderer C., Everschor K., Garst M., Rosch A., Nat. Phys., 2012, 8, 301, doi:10.1038/nphys2231.
  * [11] Sampaio J., Cros V., Rohart S., Thiaville A., Fert A., Nat. Nanotechnol., 2013, 8, 839,
doi:10.1038/nnano.2013.210.

  * [12] Jiang W., Upadhyaya P., Zhang W., Yu G., Jungfleisch M. B., Fradin F. Y., Pearson J. E., Tserkovnyak Y., Wang K. L., Heinonen O., te Velthuis S. G. E., Hoffmann A., Science, 2015, 349, 283,
doi:10.1126/science.aaa1442.

  * [13] Boulle O., Vogel J., Yang H., Pizzini S., de Souza Chaves D., Locatelli A., Mentes T. O., Sala A., Buda-Prejbeanu L. D., Klein O., Belmeguenai M., Roussigné Y., Stashkevich A., Chérif S. M., Aballe L., Foerster M., Chshiev M., Auffret S., Miron I. M., Gaudin G., Nat. Nanotechnol., 2016, 11, 449,
doi:10.1038/nnano.2015.315.

  * [14] Préjean J. J., Joliclerk M. J., Monod P., J. Phys. (Paris), 1980, 41, 427, doi:10.1051/jphys:01980004105042700.
  * [15] Hippert F., Alloul H., Préjean J. J., Physica B, 1981, 107, 427.
  * [16] Levi P. M., Fert A., Phys. Rev. B, 1981, 23, 4667, doi:10.1103/PhysRevB.23.4667.
  * [17] Levi P. M., Morgan-Pond C., Fert A., J. Appl. Phys., 1982, 53, 2168, doi:10.1063/1.330770.
  * [18] Li B., Cho S. Y., Wang H. L., Hu B. Q., J. Phys. A: Math. Theor., 2011, 44, 392002, doi:10.1088/1751-8113/44/39/392002.
  * [19] Soltani M. R., Mahdavifar S., Akbari A., Masoudi A. A., J. Supercond. Novel Magn., 2010, 23, 1369,
doi:10.1007/s10948-010-0785-x.

  * [20] Strečka J., Čanová L., J. Phys.: Conf. Ser., 2009, 145, 012012, doi:10.1088/1742-6596/145/1/012012.
  * [21] Parente W. E. F., Pacobahyba J. T. M., Arújo I. G., Neto M. A., de Sousa J. R., Akinci Ü., J. Magn. Magn. Mater., 2014, 355, 235, doi:10.1016/j.jmmm.2013.12.041.
  * [22] Liu G. H., You W. L., Li W., Su G., J. Phys.: Condens. Matter, 2015, 27, 165602, doi:10.1088/0953-8984/27/16/165602.
  * [23] Strečka J., Čanová L., Minami K., Phys. Rev. E, 2009, 79, 051103, doi:10.1103/PhysRevE.79.051103.
  * [24] Parente W. E. F., Pacobahyba J. T. M., Neto M. A., Araújo I. G., Plascak J. A., J. Magn. Magn. Mater., 2018, 462, 8, doi:10.1016/j.jmmm.2018.04.054.
  * [25] Avalishvili N., Beradze B., Japaridze G. I., Eur. Phys. J. B, 2019, 92, 262, doi:10.1140/epjb/e2019-100323-1.
  * [26] Chan Y. H., Jin W., Jiang H. C., Starykh O. A., Phys. Rev. B, 2017, 96, 214441,
doi:10.1103/PhysRevB.96.214441.

  * [27] Flynn M. O., Singh R. R. P., Phys. Rev. B, 2019, 100, 121108(R), doi:10.1103/PhysRevB.100.121108.
  * [28] Freitas A. S., de Albuquerque D. F., Phys. Rev. E, 2015, 91, 012117, doi:10.1103/PhysRevE.91.012117.
  * [29] Grandi N., Lagos M., Oliva J., Vera A., Eur. Phys. J. B, 2019, 92, 244, doi:10.1140/epjb/e2019-100395-3.
  * [30] Griset C., Head S., Alicea J., Starykh O. A., Phys. Rev. B, 2011, 84, 245108, doi:10.1103/PhysRevB.84.245108.
  * [31] Japaridze G. I., Cheraghi H., Mahdavifar S., Phys. Rev. E, 2021, 104, 014134, doi:10.1103/PhysRevE.104.014134.
  * [32] Li X., Jin J., Phys. Rev. B, 2021, 103, 035127, doi:10.1103/PhysRevB.103.035127.
  * [33] Messio L., Bieri S., Lhuillier C., Bernu B., Phys. Rev. Lett., 2017, 118, 267201,
doi:10.1103/PhysRevLett.118.267201.

  * [34] Sera A., Kousaka Y., Akimitsu J., Sera M., Kawamata T., Koike Y., Inoue K., Phys. Rev. B, 2016, 94, 214408,
doi:10.1103/PhysRevB.94.214408.

  * [35] Shindou R., Phys. Rev. B, 2016, 93, 094419, doi:10.1103/PhysRevB.93.094419.
  * [36] Szymczak H., Baran M., Szymczak R., Barilo S. N., Bychkov G. L., Shiryaev S. V., Acta Phys. Pol., A, 2007, 111, 71, doi:10.12693/APhysPolA.111.71

Ukrainian 0

Ôåðîìàãíòíà ìîäåëü Ãàéçåíáåðãà ç âçàìîäþ Äçÿëîøèíñüêîãî-Ìîðÿ E. Àëáàéðàê

Ôçèíèé ôàêóëüòåò, Óíâåðñèòåò Åðèñ, 38039 Êàéñåð, Òóðåèíà