Michele Focchi, Francesco Roscia, Claudio Semini 11institutetext: Dynamic
Legged Systems (DLS), Istituto Italiano di Tecnologia (IIT), Genoa, Italy
22institutetext: Department of Information Engineering and Computer Science
(DISI), University of Trento, Trento, Italy
33institutetext: Department of Informatics, Bioengineering, Robotics and
Systems Engineering (DIBRIS), University of Genoa, Genova, Italy
33email: {michele.focchi}, {francesco.roscia}<EMAIL_ADDRESS>
${}^{*}\,$The authors equally contributed to this paper

# Locosim: an Open-Source Cross-Platform Robotics Framework

Michele Focchi 1122** Francesco Roscia 1133** Claudio Semini 11

###### Abstract

The architecture of a robotics software framework tremendously influences the
effort and time it takes for end users to test new concepts in a simulation
environment and to control real hardware. Many years of activity in the field
allowed us to sort out crucial requirements for a framework tailored for
robotics: modularity and extensibility, source code reusability, feature
richness, and user-friendliness. We implemented these requirements and
collected best practices in Locosim, a cross-platform framework for simulation
and real hardware. In this paper, we describe the architecture of Locosim and
illustrate some use cases that show its potential.

###### keywords:

Computer Architecture for Robotics; Software Tools for Robot Programming;
Software-Hardware Integration for Robot Systems

##  1 Introduction

Writing software for robotic platforms can be arduous, time-consuming, and
error-prone. In recent years, the number of research groups working in
robotics has grown exponentially, each group having platforms with peculiar
characteristics. The choice of morphology, actuation systems, and sensing
technology is virtually unlimited, and code reuse is fundamental to getting
new robots up and running in the shortest possible time. In addition, it is
pervasive for researchers willing to test new ideas in simulation without
wasting time in coding, for instance, using high-level languages for rapid
code prototyping. To pursue these goals, in the past years several robotics
frameworks have been designed for teaching or for controlling specific
platforms, e.g., OpenRAVE [1], Drake [2] and SL [3].

To avoid roboticists reinventing the wheel whenever they buy or build a new
robot, we present our framework Locosim111Locosim can be downloaded from
www.github.com/mfocchi/locosim.. Locosim is designed with the primary goal of
being platform-independent, dramatically simplifying the task of interfacing a
robot with planners and controllers. Locosim consists of a ROS control node
[4] (the low-level controller), written in C++, that interfaces a custom
Python ROS node (the high-level planner/controller) with either a Gazebo
simulator [5] or the real hardware. The planner relies on Pinocchio [6] for
computing the robot’s kinematics and dynamics and closes the control loop at a
user-defined frequency.

Figure 1: Examples of robots already included in Locosim (from left to right,
top to bottom): Aliengo [7], Go1 [8], HyQ [9], Starbot, CLIO [10], UR5 [11]
with Gripper, Solo with Flywheels [12], Tractor (images are not in scale).

###  1.1 Advantages of Locosim

The benefits of the proposed framework are multiple.

  * •

Locosim is platform-independent. It supports a list of robots with different
morphology (e.g., quadrupeds, arms, hybrid structures, see Fig. 1) and
provides features for fast designing and adding new robots.

  * •

Locosim implements functions needed for all robots. Once the robot description
is provided, no effort is spent on libraries for kinematics, dynamics,
logging, plotting, or visualization. These valuable tools ease the synthesis
of a new planner/controller.

  * •

Locosim is easy to learn. The end user invests little time in training and
gets an open-source framework with all the benefits of Python and ROS.

  * •

Locosim is modular. Because it heavily uses the inheritance paradigm, classes
of increasing complexity can provide different features depending on the
nature of the specific robotic application. For instance, the controller for
fixed-base robotic arms with a grasping tool can be reused for a four-legged
robot with flywheels since it is built upon the same base class.

  * •

Locosim is extensible. Our framework is modifiable and the end-user can add
any supplementary functionality.

  * •

Locosim is easy to install. It can be used either inside a virtual machine, a
docker container, or natively by manually installing dependencies.

###  1.2 Outline

The remainder of this paper is organized as follows. In Section 2 we highlight
the critical requirements of a cross-platform robotics framework. In Section 3
we detail structure and features of Locosim. In Section 4 we discuss use-case
examples of our framework, either with the real robot or with its simulated
digital twin. Eventually, we condense the results and present future works in
Section 5.

Figure 2: End-users, robotic platform and simulation environment make a triad
only if an effective robotics framework can join them.

##  2 Key aspects of a robotics framework

In the most general sense, a robotics framework is a software architecture of
programs and data that adhere to well-defined rules for operating robots. End-
users are people who will ultimately use the robotics framework and
potentially bring some modifications. The robotics framework is the center of
a triangle having at its vertices the end-user, the robotic platform and the
simulation environment (see Fig. 2). The simulation must replicate the
behaviour of the robot with a sufficient degree of accuracy. The end-user must
be able to test new features and ideas in the simulation environment before
running them on the robot platform. A robotics framework should provide the
link among these three. In this context, we identify a list of crucial
requirements that a robotics framework must possess: generality, modularity,
reusability, extensibility, rapid prototyping vs. performances, feature-rich,
and end-users development.

####  2.0.1 Generality.

It is essential to release software free from unnatural restrictions and
limitations. An end-user may require to design software for controlling a
single robot, for multi-robot cooperation, for swarm coordination, or for
reinforcement learning, which requires an abundant number of robots in the
training phase. It should be possible to model any kinematic structure
(floating base/fixed base robot, kinematic loops, etc.).

####  2.0.2 Modularity.

A robotics framework should provide separate building blocks of
functionalities according to the specific needs of the robot or the
application. These building blocks must be replaceable and separated by clear
interfaces. The end-user may want to only visualize a specific robot in a
particular joint configuration, move the joints interactively, or plan motions
and understand the effects of external forces. Each of these functions must be
equipped with tools for debugging and testing, e.g., unit tests. Replacing
each module with improved or updated versions with additional functionalities
should be easy.

####  2.0.3 Reusability.

Pieces of source code of a program should be reused by reassembling them in
more complex programs to provide a desired set of functionalities, with minor
or no modifications. From this perspective, parametrization is a simple but
effective method, e.g., the end-user should be able to use the same framework
with different robots changing only a single parameter. In the same way,
digital twins [13] can be realized by varying the status of a flag that
selects between the real hardware and the simulation environment. This avoids
writing different codes for the simulator and the real robot.

####  2.0.4 Extensibility.

A robotics framework must be designed with the foresight to support the
addition and evolution of hardware and software that may not exist at
implementation time. This property can be achieved by providing a general set
of application programming interfaces (APIs). Concepts typical of Object-
Oriented Programming, such as inheritance and polymorphism, play a crucial
role in extensibility.

####  2.0.5 Rapid Prototyping vs. Performances.

A framework should allow for fast code prototyping of researchers’ ideas. More
specialized controllers/planners are built from simpler ones in the form of
recursive inheritance (matryoshka principle). In this way end–users have unit
tests at different levels of complexity. With fast code prototyping, end-users
can quickly write software without syntax and logic errors. However, they do
not have any assurance about the performance: such code is just good enough
for the case of study in a simulation environment. Stringent requirements
appear when executing codes on real robots, e.g., short computation time and
limited memory usage. Thus, the framework must expose functionalities that can
deal with performance.

####  2.0.6 Feature-rich.

Most of the end-users need a sequence of functionalities when working with
robots. These include but are not limited to the computation of kinematics and
dynamics, logging, plotting, and visualization. A robotics framework should
provide them, and they must be easily accessible.

####  2.0.7 End-users Development.

Besides its implementation details, a robotics framework should provide
methods, techniques and tools that allow end-users to create, modify, or
extend the software [14] in an easy way. It should run on widely used
Operating Systems and employ renowned high-level programming languages that
facilitate software integration. Clear documentation for installation and
usage must be provided, and modules should have self-explanatory names and
attributes.

##  3 Locosim Description

Locosim was born as a didactic framework to simulate both fixed- and floating-
base robots. Quickly it evolved to be a framework for researchers that want to
program newly purchased robots in a short time. Locosim runs on machines with
Ubuntu as Operating System, and it employs ROS as middleware. Within Locosim,
end-users can write robot controllers/planners in Python.

###  3.1 Architecture

Locosim consists of four components: Robot Descriptions, Robot Hardware
Interfaces, ROS Impedance Controller, and Robot Control. We illustrate each
component in the following111Some of the functions in Locosim components,
which are quite established for the robotics community, are named after [3]..

####  3.1.1 Robot Descriptions.

The Robot Descriptions component contains dedicated packages for the
characterization of each robot. For instance, the package that contains files
to describe the robot myrobot, a generic mobile robot platform, is named
myrobot_description. With the main focus on fast prototyping and human
readability, the robot description is written in Xacro (XML-based scripting
language) that avoids code replication through macros, conditional statements
and parameters. It can import descriptions of some parts of the robot from
urdfs sub-folder or meshes describing the geometry of rigid bodies from the
meshes sub-folder. At run time, the Xacro file is converted into URDF,
allowing the end-user to change some parameters. The gazebo.urdf.xacro
launches ros_control package and the gazebo_ros_p3d plugin, which publishes
the pose of the robot trunk in the topic /ground_truth (needed only for
floating-base robots). The robot description directory must contain the files
upload.launch and rviz.launch. The former processes the Xacro generating the
URDF and loads into the parameter server, and the latter allows to visualize
the robot and interact with it by providing the conf.rviz file. This is the
configuration for the ROS visualizer RViz, which can be different for every
robot. Additionally, the file ros_impedance_controller.yaml must be provided
for each robot: it contains the sequence of joint names, joint PID gains and
the home configuration. The Python script go0 will be used during the
simulation startup to drive the robot to the home configuration.

####  3.1.2 Robot Hardware Interfaces.

This folder contains drivers for the real hardware platforms supported by
Locosim. They implement the interface that bridges the communication between
the controller and the real robot, abstracting the specificity of each robot
and exposing the same interface (e.g., EffortInterface). For instance, the UR5
robot through its driver provides three possible ROS hardware interfaces:
Effort, Position and Velocity, hiding the details of the respective underlying
low-level controllers.

Table 1: Main attributes and methods of the BaseControllerFixed (BCF) and of the BaseController (BC) classes. All vectors (unless specified) are expressed in world frame. For methods with the same name, the derived class loads the method of the parent class and adds additional elements specific to that class. | Name |  |  Meaning |  | Class
---|---|---|---|---|---
Attributes |  q, q_des |  |  actual / desired joint positions |  | BCF, BC
qd, qd_des |  |  actual /desired joint velocities |  | BCF, BC
tau, tau_ffwd |  |  actual / feed-forward joint torques |  | BCF, BC
x_ee |  |  position of the end-effector expressed in base frame |  | BCF
contactForceW |  |  contact force at the end-effector |  | BCF
contactMomentW |  |  contact moment at the end-effector |  | BCF
basePoseW |  |  base position and orientation in Euler angles |  | BC
baseTwistW |  |  base linear and angular velocity |  | BC
b_R_w |  |  orientation of the base link |  | BC
contactsW |  |  position of the contacts |  | BC
grForcesW |  |  ground reaction forces on contacts |  | BC
Methods | loadModelAndPublishers() |  |  creates the object robot (Pinocchio wrapper) and loads publishers for visual features (ros_pub), joint commands and declares subscriber to /ground_truth and /joint_states |  | BCF, BC
startFramework() |  |  launch ros_impedance_controller |  | BCF, BC
send_des_jstate() |  |  publishes /command topic with set-points for joint positions, velocities and feed-forward torques |  | BCF, BC
startupProcedure() |  |  initializes PID gains |  | BCF, BC
initVars() |  |  initializes class attributes |  | BCF, BC
logData() |  |  fill in the variables x_log with the content of x for plotting purposes, it needs to be called at every loop |  | BCF, BC
receive_jstate() |  |  callback associated to the subscriber /joint_states, fills in actual joint positions, velocities and torques |  | BCF, BC
receive_pose() |  |  callback associated to the subscriber /ground_truth, fills in the actual base pose and twist, and publishes the fixed transform between world and base |  | BC
updateKinematics() |  |  from q, qd, basePoseW, baseTwistW, computes position and velocity of the robot’s center of mass, the contacts position, the associated Jacobians, the ground reaction forces, the centroidal inertia, the mass matrix and the non-linear effects |  | BC
Figure 3: Schematic representation of a typical use-case of Locosim. The end-
user wants to simulate the UR5 robot arm. An instance of Ur5Generic, which is
a derived class of BaseControllerFixed, sends the command to the robot though
ros_impedance_controller and it receives back the actual state. Ur5Generic
implements features to manage the real robot and the gripper, perform a homing
procedure at startup and a class for inverse kinematics.

####  3.1.3 ROS Impedance Controller.

ROS Impedance Controller is a ROS package written in C++ that implements the
low-level joint controller, in the form of PID with feedforward effort (force
or torque). The /joint_state_publisher publishes the actual position, velocity
and effort for each robot joint. By default, the loop frequency is set to 1
kHz. This and other parameters can be regulated in the launch file called
ros_impedance_controller.launch. Robots with specific needs can be dealt with
by specifying a custom launch file. This is the case of the CLIO climbing
robot that requires the model of the mountain to which it is attached. The
robot_name parameter is used to load the correct robot description. If the
real_robot flag is set to true, the robot hardware interface is initialized;
otherwise, the Gazebo simulation starts running first the go0 script from the
robot description. In any case, Rviz will be opened with the robot’s specific
configuration file conf.rviz. The end-user can manually change the location
where the robot is spawned in Rviz with the spawn parameters. Physics
parameters for the simulator are stored in the sub-folder worlds.

####  3.1.4 Robot Control.

From the end-user perspective, this is the most crucial component. It embraces
classes and methods for computation of the robot’s kinematics and dynamics,
logging and plotting time series variables, and real-time visualization on
Rviz. Within this component, high-level planning/control strategies are
implemented. The codes of the component are entirely written in Python and
have a few dependencies: above all, NumPy [15] and Pinocchio [6]. The former
offers tools for manipulating multidimensional arrays; the latter implements
functions for the robot’s kinematics and dynamics. Pinocchio can be an
essential instrument for researchers because of its efficient computations.
Nevertheless, it can be time-consuming and cumbersome to understand for
newcomers. To facilitate the employment, we developed a custom_robot_wrapper
for building a robot object, computing robot mass, center of mass, Jacobians,
centroidal inertia, and so on with easy-to-understand interfaces.
For the end-user, the starting point for building up a robot planner is the
class BaseControllerFixed, suitable for fixed-base robots, and the derived
class BaseController, which handles floating-base ones. In Table 1, we report
a list of the main attributes and methods of BaseController and of its parent
BaseControllerFixed. For having a more complex and specific controller, the
end-user can create its own class, inheriting from one of the previous two and
adding additional functionalities. E.g., QuadrupedController inherits from
BaseController, and it is specific for quadruped robots. Ur5Generic adds to
BaseControllerFixed the features of controlling the real robot, a gripper and
a camera attached (or not) to the robotic arm UR5 (see Fig. 3). The controller
class is initialized with the string robot_name, e.g., we write
BaseController(myrobot) for the controller of myrobot. The end-user must pay
particular attention to this string because it creates the link with the robot
description and the robot hardware interface if needed. The robot_name is used
for accessing the dictionary robot_param too. Among the other parameters of
the dictionary, the flag real_robot permits using the same code for both the
real (if set to true) and simulated (false) robot, resulting in the digital
twin concept. The BaseControllerFixed class contains a ControllerManager to
seamlessly swap between the control modes and controller types if the real
hardware supports more than one. For instance, UR5 has two control modes
(point, trajectory) and two controller types (position, torque), whereas Go1
supports a single low-level torque controller. Additionally, the
GripperManager class manages the gripper in different ways for the simulation
(i.e., adding additional finger joints) and for the real robot (on-off
opening/closing service call to the UR5 driver as specified by the
manufacturer), hiding this complexity to the end-user. The method
startFramework() permits to launch the simulation or the driver, executing
ros_impedance_controller.launch. It takes as input the list additional_args to
propagate supplementary arguments, dealing with robot and task specificity. In
the components folder there are additional classes for inverse kinematics,
whole-body control, leg odometry, IMU utility, filters and more. Finally, the
folder lab_exercises contains an ample list of scripts, with didactic
exercises of incremental complexity, to learn Locosim and the main robotics
concepts.

###  3.2 Analysis of Design Choices

To fulfill the requirements stated in Section 2, we made a number of choices.
We want to focus on why we selected ROS as middleware, Python as (preferred)
programming language, and Pinocchio as computational tool for the robot’s
kinematics and dynamics.

####  3.2.1 Why ROS.

The ROS community is spread worldwide. Over the last decade, the community
produced hundreds or thousands of open-source tools and software: from device
drivers and interface packages of various sensors and actuators to tools for
debugging, visualization, and off-the-shelf algorithms for multiple purposes.
With the notations of nodes, publishers, and subscribers, ROS quickly solves
the arduous problem of resource handling between many processes. ROS is
reliable: the system can still work if one node crashes. On the other hand,
learning ROS is not an effortless task for newcomers. Moreover, modeling
robots with URDF can take lots of time, as well as starting simulations [16].
Locosim relieves end-users from these complications by adopting a common
skeleton infrastructure for the robot description and for the high-level
planner/controller.

####  3.2.2 Why Python.

Among the general-purpose programming languages, Python is one of the most
used. It is object-oriented and relatively straightforward to learn and use
compared to other languages. The availability of open-source libraries and
packages is virtually endless. These reasons make Python perfect for fast
prototyping software. Being an interpreted language, Python may suffer in
terms of computation time and memory allocation, colliding with the real-time
requirements of the real hardware. In these cases, when testing the code in
simulation, the end-user may consider using profiling tools to look for the
most demanding parts of the code. Before executing the software on the real
robot, these parts can be optimized within Python or translated into C++ code,
providing Python bindings. For the same performance reasons, the most critical
part of the framework, the low-level controller, is directly implemented in
C++.

####  3.2.3 Why Pinocchio.

Pinocchio [6] is one of the most efficient tools for computing poly-
articulated bodies’ kinematics and dynamics. Differently from other libraries
in which a meta-program generates some source code for a single specific robot
model given in input, as in the case of RobCoGen [17], Pinocchio is a dynamic
library that loads at runtime any robot model. This characteristic makes it
suitable for a cross-platform framework. It implements the state-of-the-art
Rigid Body Algorithms based on revisited Roy Featherstone’s algorithms [18].
Pinocchio is open-source, mostly written in C++ with Python bindings and
several projects rely on it. Pinocchio also provides derivatives of the
kinematics and dynamics, which are essential to in gradient-based
optimization.

##  4 Use Cases

We want to emphasize the valuable features of Locosim with practical use
cases222A video showing the above-mentioned and other use cases can be found
here: https://youtu.be/ZwV1LEqK-LU .

###  4.1 Visualize a Robot: kinematics check

As first use case, we illustrate the procedure to visualize the quadruped
robot Aliengo (see Fig. 1) in RViz and how to manually interact with its
joints. This is a debugging tool which is crucial during the design process of
a new robot, because it allows to test the kinematics without added
complexity. In the Robot Description package, we create a folder named
aliengo_description. In the robots folder we add the XML file for describing
the robot’s kinematic and dynamic structure. We make use of the flexibility of
the open source Xacro language to simplify writing process: we include files
that describe a leg, the transmission model, and the meshes for each of the
bodies. We create the launch folder containing the files upload.launch and
rviz.launch. Launching rviz.launch from command line, the end-user can
visualize the robot in RViz and manually move the joints by dragging the
sliders of the joint_state_publisher_gui. The conf.rviz file helps the end-
user to restore previous sessions in RViz. With this simple use case we can
effectively understand the importance of the key aspects formalized so far.
Being extensible, Locosim allows for the introduction of any new robots,
reusing parts of codes already present.

Figure 4: Execution of the pick-and-place task with the anthropomorphic arm
UR5. The end-user can drive the real hardware (setting real_robot to true) or
perform a simulation (real_robot to false).

###  4.2 Simulation and Real Robot

As a second example, we present a pick-and-place task with the anthropomorphic
arm UR5 (6 degrees of freedom, see Fig. 1). The pipeline of planning and
control starts with the launch file ros_impedance_controller.launch. This is
common for all the robots: it loads the robot_description parameter calling
the upload.launch and it starts the Gazebo simulator or the robot driver
according to the status of the real_robot flag. Additionally, it loads the
controller parameters (e.g., PID gains), which are located in the
robot_description package of the robot. In the simulation case, the launch
file spawns the robot at the desired location. In the real robot case, the
robot driver is running (in a ROS node) on an onboard computer while the
planner runs on an external computer. In both cases, the RViz GUI shows the
robot configuration. Another node running on the same computer reads from a
fixed frame ZED2 camera and publishes a point cloud message of the environment
on a ROS topic. We extract the coordinates of the plastic bricks that are
present in the workspace. With Ur5Generic, we plan trajectories for the end-
effector position and orientation to grasp and relocate the bricks. We set an
inverse kinematic problem to find a joint reference trajectory. It is
published in the /ur5/command topic, together with feed-forward torques for
gravity compensation. The low-level ros_impedance_controller provides feedback
for tracking the joint references, based on the actual state in /joint_state.
On the real robot there is no torque control and only position set-points are
provided to the ur5/joint_group_pos_controller/command topic as requested by
the robot driver provided by the manufacturer. All this is dealt with by the
controller_manager class, transparent to the end users. The power of Locosim
lies on the fact that it is possible to use the same robot control code both
to simulate the task and to execute it on the real robot, as reported in Fig.
4.

##  5 Conclusions

Locosim is a platform-independent framework for working with robots, either in
a simulation environment or with real hardware. Integrating features for
computation of robots’ kinematics and dynamics, logging, plotting, and
visualization, Locosim is a considerable help for roboticists that needs a
starting point for rapid code prototyping. If needed, performances can be
ensured by implementing the critical parts of the software in C++, providing
Python bindings. We proved the usefulness and versatility of our framework
with use cases. Future works include the following objectives: Locosim will be
able to handle multiple platforms simultaneously to be used in the fields of
swarm robotics and collaborative robotics. We want to provide support for ROS2
since ROS lacks relevant qualifications such as real-time, safety,
certification, and security. Additionally, other simulator like Coppelia Sim
or PyBullet will be added to Locosim.

## References

  * [1] “OpenRAVE.” [Online]. Available: https://github.com/rdiankov/openrave
  * [2] “Drake.” [Online]. Available: https://github.com/RobotLocomotion/drake
  * [3] S. Schaal, “The sl simulation and real-time control software package,” Los Angeles, CA, Tech. Rep., 2009, clmc. [Online]. Available: http://www-clmc.usc.edu/publications/S/schaal-TRSL.pdf
  * [4] M. Quigley, K. Conley, B. Gerkey, J. Faust, T. Foote, J. Leibs, R. Wheeler, A. Y. Ng _et al._ , “Ros: an open-source robot operating system,” in _ICRA workshop on open source software_ , vol. 3, no. 3.2. Kobe, Japan, 2009, p. 5.
  * [5] N. Koenig and A. Howard, “Design and use paradigms for gazebo, an open-source multi-robot simulator,” in _2004 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)(IEEE Cat. No. 04CH37566)_ , vol. 3. IEEE, 2004, pp. 2149–2154.
  * [6] J. Carpentier, G. Saurel, G. Buondonno, J. Mirabel, F. Lamiraux, O. Stasse, and N. Mansard, “The Pinocchio C++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives,” in _2019 IEEE/SICE International Symposium on System Integration (SII)_. IEEE, 2019, pp. 614–619.
  * [7] “Aliengo.” [Online]. Available: https://www.unitree.com/en/aliengo
  * [8] “Go1.” [Online]. Available: https://www.unitree.com/en/go1
  * [9] C. Semini, N. G. Tsagarakis, E. Guglielmino, M. Focchi, F. Cannella, and D. G. Caldwell, “Design of hyq–a hydraulically and electrically actuated quadruped robot,” _Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering_ , vol. 225, no. 6, pp. 831–849, 2011.
  * [10] M. Focchi, M. Bensaadallah, M. Frego, A. Peer, D. Fontanelli, A. Del Prete, and L. Palopoli, “Clio: a novel robotic solution for exploration and rescue missions in hostile mountain environments,” _arXiv preprint arXiv:2209.09693_ , 2022.
  * [11] “UR5.” [Online]. Available: https://www.universal-robots.com/products/ur5-robot/
  * [12] F. Roscia, A. Cumerlotti, A. Del Prete, C. Semini, and M. Focchi, “Orientation control system: Enhancing aerial maneuvers for quadruped robots,” _Sensors_ , vol. 23, no. 3, p. 1234, 2023.
  * [13] A. K. Ramasubramanian, R. Mathew, M. Kelly, V. Hargaden, and N. Papakostas, “Digital twin for human–robot collaboration in manufacturing: Review and outlook,” _Applied Sciences_ , vol. 12, no. 10, p. 4811, 2022.
  * [14] H. Lieberman, F. Paternò, and V. Wulf, _End user development_. Springer, 2006, vol. 9.
  * [15] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with NumPy,” _Nature_ , vol. 585, no. 7825, pp. 357–362, Sep. 2020. [Online]. Available: https://doi.org/10.1038/s41586-020-2649-2
  * [16] L. Joseph and J. Cacace, _Mastering ROS for Robotics Programming: Design, build, and simulate complex robots using the Robot Operating System_. Packt Publishing Ltd, 2018.
  * [17] M. Frigerio, J. Buchli, D. G. Caldwell, and C. Semini, “Robcogen: a code generator for efficient kinematics and dynamics of articulated robots, based on domain specific languages,” _Journal of Software Engineering for Robotics (JOSER)_ , vol. 7, no. 1, pp. 36–54, 2016.
  * [18] R. Featherstone, _Rigid body dynamics algorithms_. Springer, 2014.


# Rotational symmetry breaking potential for two-dimensional magnets

Cláudio J. DaSilva<EMAIL_ADDRESS>Instituto Federal de Goiás, Rua
76, Centro, Goiânia - GO, Brazil  L. S. Ferreira  Instituto de Física,
Universidade Federal de Goiás, Av. Esperança s/n, 74.690-900, Goiânia, GO,
Brazil  A. A. Caparica  Instituto de Física, Universidade Federal de Goiás,
Av. Esperança s/n, 74.690-900, Goiânia, GO, Brazil

###### Abstract

Here we present a new perspective to the breakdown of ferromagnetic order in
two-dimensional spin-lattice models employing the rotation of the underlying
lattice. Using an Ising spin system on a square lattice as a prototype, we
demonstrate that an additional low-symmetry interaction may lead to the
absence of the truly long-range order and forms aperiodic structure, such as
magnetic stripes. Employing annealing and entropic Monte Carlo simulations, we
show that our model allows tuning between different phases, magnetically
ordered as well as more exotic nonmagnetic phases such as Ising-nematic by
changing only one control parameter, which is responsible for the arising of
magnetic frustration. In addition, our methodology of considering the coupling
between the magnetic structure and the host material can be extended to the
study of any type of spin-exchange model in two dimensions and has many
potential interesting ramifications and applications.

The last several years have witnessed a renewed interest in frustrated
magnetism[1, 2, 3, 4, 5], which was propelled by the discovery of actual
compounds presenting new magnetic exotic phases which can be described by two-
dimensional spin-exchange models[6]. These novel excitations inherent to
frustrated systems include the possible quantum spin-liquid state[7] and
chiral ordered states[8].

Experimentally observable, materials containing layers of transition metal
ions on a square lattice are considered quasi-two-dimensional (quasi-2D) for
magnetism, as the interactions between layers are much smaller than those
within the layers. At sufficiently low energies their magnetic degrees of
freedom can be mapped to effective spin-exchange models into 2D Bravais
lattices. Commonly, the expected long-range magnetic order is suppressed by
magnetic frustration, which arises when there is no unique ground state which
simultaneously minimizes the energy of all exchange bonds. Generally, spin-
exchange frustration can be due to either the geometry of the lattice (as in
pyrochlore[9, 10], spinel, or Kagome systems) or where the competing
interactions between magnetic moments cannot be mutually satisfied (as the
competing nearest-neighbor and next-nearest-neighbor interactions). In its
turn, frustration enhances spin fluctuations, and novel complex phase
diagrams, new emergent critical phenomena, and technical applications can
emerge.

Nonetheless, mechanical vibrations of the crystal lattice (phonons) modulate
the crystal electric field of the magnetic ion, thus inducing a direct
relaxation between two different spin states, and unconventional magnetic
properties could be linked to phonon renormalization. It has been known for a
while that spin-phonon coupling induces complex collinear order in triangular
antiferromagnets[11] and noncollinear order in pyrochlore magnets [10, 9]. For
strong spin-lattice coupling, striped collinear states can be observed in
compounds such as $\alpha$NaFeO2 and MnBr2[11]. Yet, despite the vast number
of models describing spatial anisotropy of spin fluctuations[12, 5, 13], a
general lattice spin model considering the effects of spin-host lattice
excitation coupling is still missing.

Therefore, here we present a new perspective on the breakdown of ferromagnetic
order in two-dimensional lattice spin models through the rotation of the
underlying lattice. Our model adds a rotational symmetry-breaking potential
native to the host material structure, which is arguably necessary to orient
magnetic phases in macroscopic samples. We will show that such a model
possesses ground states with various patterns such as stripes, checkerboards,
and more complicated morphologies such as Ising-nematic. All these phases are
tuned by changing only a single controlling parameter.

Altogether, one can write a Hamiltonian for a model whose interaction depends
on the symmetry of the layer as follows

$\displaystyle\mathcal{H}=-\sum_{ij}J_{ij}^{\text{eff}}\vec{S}_{i}\cdot\vec{S}_{j}$
(1)

with $J_{ij}^{\text{eff}}=J+J_{a}\cos{n\phi_{ij}}$ and $J$ is the regular
magnetic exchange coupling. The constant $J_{a}$, determined by the strength
of spin coupling to acoustic phonons, will be considered here as a
phenomenological parameter. $\phi_{ij}$ is the angle between one of the plane
axes and the line joining any two spins $i$ and $j$ as illustrated in Fig. 1.
$n$ is an integer number related to the degree of symmetry of the underlying
crystal structure, namely $n=2$ for a tetragonal symmetry and $n=4$ for a
square symmetry. Generally, the dependence of exchange on displacements may
involve changes in the bond angles as well as distances[14, 15, 11]. Such a
theoretical procedure has already been successfully applied to strongly
correlated system in Refs. [16, 17]. Besides, Ref. [18] considered magnon-
phonon interactions in Iron and observed that the effect of phonons on magnons
depends strongly on the symmetry of the underling lattice. Also, Ref. [19]
proposed a Heisenberg model that couples the spin and lattice degrees of
freedom for BCC iron. They have demonstrated by using Monte Carlo simulations
that a higher transition temperature is obtained when a background
contribution from phonons is considered.

In the following, we will consider a single square lattice with only Ising
spin states $\sigma=\pm 1$ with ferromagnetic and the spin-phonon interactions
given by Eq. 1. Hence, the dimensionless interaction energy will be
$-\sum_{\langle ij\rangle}(1+\kappa\cos{2\phi_{ij}})\sigma_{i}\sigma_{j}$ for
the tetragonal symmetry and $-\sum_{\langle
ij\rangle}(1+\kappa\cos{4\phi_{ij}})\sigma_{i}\sigma_{j}$ for the square
symmetry, where $\kappa=J_{a}/J$. $\langle ij\rangle$ indicates that the sum
runs over all pairs of nearest-neighbor spins. Here, the parameter $\kappa$
quantifies the level of magnetic frustration. Although the geometry of the
square lattice does not lead to frustration, the relative strengths of $J$ and
$J_{a}$ determine the magnetic ground state, and a frustrated state can arise
from their competition. The case of a square lattice deserves our initial
attention since recent inelastic neutron scattering measurements confirmed an
effective one-dimensional coupling in the highly frustrated itinerant magnet
CaCo2-yAs2[6] where the authors have identified Néel-(or checkerboard)type
antiferromagnetic (AFM), stripe-type AFM, and ferromagnetic (FM)-type ordering
occurring and been suppressed near the region of maximum frustration.
Remarkably, within our approach we were also able to identify all those phases
including a spin-nematic state at the maximum frustration which occurs when
$\kappa=\pm 1$ for $n=2$ and $\kappa=-1$ for $n=4$.

Figure 1:  A schematic illustration of the interaction symmetry.

First, we focus on the zero-temperature phase diagram as a function of the
single parameter $\kappa$ obtained by two different numerical simulations
techniques, namely Monte Carlo annealing and entropic sampling. Simulated
annealing (SA) on each size was done using a uniform slow annealing schedule
as follows: The simulation consists of $M$ segments of $N$ Monte Carlo steps
each so that the total number of Monte Carlo steps is $NM$. The inverse
temperature is increased after each Monte Carlo step as
$\beta=1/k_{B}T=\beta_{0}+i\Delta\beta_{M}\quad(i=1,...,NM)$, where
$\Delta\beta_{M}=(\beta_{M}-\beta_{0})/NM$, with
$\beta_{0}=0.0,\beta_{M}=8.0,M=300$ and $N=10^{6}$ steps of single spin flips
updates with Metropolis dynamics. We executed ten independent annealing runs
to ensure the stability of the results.

By its turn, the entropic Monte Carlo simulations, based on a modified Wang-
Landau (WL) scheme[20], are mainly applied to study the thermodynamics of the
system at the critical values of $|\kappa|=1$. In the WL method, a random walk
is performed in the energy space, which allows estimating the density of
states $g(E)$, while a flat histogram for the energy distribution is
generated. With the density of states in hands, one can calculate the
canonical averages of any thermodynamic quantities. In our simulations, some
easily implementable changes to the method are included, which lead to
improved accuracy and to noticeable savings in CPU time. In particular, ($i$)
the density of states is updated only after every Monte Carlo sweep, such that
we discard very correlated configurations, ($ii$) the microcanonical averages
are accumulated beginning from the eighth WL level ($f_{7}$), so that we avoid
considering the initial configurations that do not match with those of maximum
entropy[21], ($iii$) a checking parameter $\varepsilon$ is used for halting
the simulation [22] (the integral of the specific heat over a range of
temperatures is calculated using the current density of states during the
simulations and the simulations are halted if this quantity varies less then
$10^{-4}$ during a whole WL level), and ($iv$) a single run is performed for
all lattice sizes up to the Wang-Landau level $f_{6}$ and then, the further
simulations begin from these outputs, since up to this point the current
density of states is not biased yet and the results we reach are similar to
those that would be obtained beginning from the first WL level $f_{0}$[23], a
measure that saves about 60% in computational time.

Figure 2: Zero-temperature phase diagram for the symmetries $n=2$ and $4$
obtained using both WL and SA procedures. For $n=2$, on increasing the
frustration parameter $\kappa$, an stripe XS state is followed first by a
fully frustrated state (shaded area), second by a ferromagnetic FM state,
third by another fully frustrated region, then an stripe YS state. For $n=4$
an antiferromagnetic/checkerboard (AF) state is followed first by a fully
frustrated state and then a ferromagnetic (FM) state.

The phase diagrams shown in Fig. 2 were obtained using lattices with periodic
boundary condition and size $42\times 42$, and choosing the state with the
lowest energy per site $E_{0}/JN$. With our numerical procedures, we have
identified four (two) distinct ordering phases of the ground states for $n=2$
($n=4$). For $n=2$, the stripes morphology – i.e., magnetic domains displayed
in stripes of alternating spins – becomes energetically favorable compared to
the uniform ferromagnetic when $|\kappa|>1$. As the value of the frustration
parameter increases, the stripe XS phase is followed first by a typical Ising
ferromagnetic FM state (with $E_{0}/JN=-2.0$), and then a stripe YS state. For
$n=4$, a checkerboard pattern associated with an antiferromagnetic AF ground
state is identified in the region where $\kappa<-1$ followed by a
ferromagnetic FM state for $\kappa>-1$. More precisely, SA results reveal a
frustration region (shaded area) when $0.75<|\kappa|<1.25$ with a maximum
frustration at $|\kappa|=1$. In this region the stripe width $h$ changes with
$\kappa$. In this range of $\kappa$, striped configurations of width $h>1$
coexist with the ferromagnetic state and a first-order phase transition seems
to take place, as we shall see later. While the SA algorithm does not
guarantee convergence to the ground state (as one can see the system get stuck
in a value of $E_{0}/JN$ close to zero for $|\kappa|=1$.), the WL procedure
always finds an ordered phase even in the frustrated region where highly
degenerated ground states emerge. We nevertheless believe that an accurate
picture emerges from the SA procedure. A frustration region is also identified
for $n=4$ case at $-1.25<\kappa<-0.75$. However, the maximum frustration
condition produces a system of non-interacting spins since $E_{0}/JN=0.0$.
Hence, a thermodynamical analysis in this case will become complete when
considering external fields coupled to all the spins, which is out of the
scope of the present paper.

Figure 3: Typical equilibrium spin configurations of the ground state (top)
for both symmetries and close to an order-order transition (bottom left) for
$n=2$. The presence of high density of topological defects, mainly
dislocations in the directions of the underlying striped structure (bottom
middle) is in agreement with the qualitative description of the Ising nematic
phase.

Fig. 3 shows typical equilibrium spin configurations of the non FM ground
state patterns, namely stripes and checkerboard, for both symmetries obtained
after a WL simulation[24]. Also, it is shown some configurations in the region
of maximum frustration for $n=2$. At the maximum frustration and low
temperatures, stripes and ferromagnetic phases coexist. As the temperature
increases, the striped structure acquires several dislocations in the
directions of its underlying structure. The presence of such defects reduces
the orientational order, which is in agreement with the qualitative
description of the Ising nematic phase[25]. With a further increase of the
temperature, domains with predominantly square corners pop out before the
system reaches the paramagnetic state. Such a phase is known as a tetragonal
liquid[26].

To better characterize the different ordering phases and possibly the order-
disorder phase transition proper order parameters were assigned to different
values of $\kappa$. The XS and YS stripe phases observed for $n=2$ are
characterized by a two-component order parameter
$m_{stp}^{2}=m_{i}^{2}+m_{j}^{2}$ with
$m_{i}=N^{-1}\sum_{i=1}^{L}\sum_{j=1}^{L}\sigma_{i,j}(-1)^{i}$,
$m_{j}=N^{-1}\sum_{i=1}^{L}\sum_{j=1}^{L}\sigma_{i,j}(-1)^{j}$, and $(i,j)$ is
the coordinate of a site on the lattice[5]. Hence, the stripe susceptibility
is defined as $\chi_{stp}=N(\langle
m_{stp}^{2}\rangle-\langle|m_{stp}|\rangle^{2})/T$. For $n=4$ case and
$\kappa<-1$, a staggered magnetization
$m_{stg}=N^{-1}\sum_{i=1}^{L}\sum_{j=1}^{L}\sigma_{i,j}(-1)^{i+j}$ and
consequently the staggered susceptibility $\chi_{stg}=N(\langle
m_{stg}^{2}\rangle-\langle|m_{stg}|\rangle^{2})/T$. For both symmetries in the
region of maximum frustration, we defined an order parameter as following
$m_{\kappa}=\sum_{j=1}^{L}\lvert\sum_{i=1}^{L}\sigma_{i,j}\rvert$ with the
corresponding susceptibility $\chi_{\kappa}=N(\langle
m_{\kappa}^{2}\rangle-\langle|m_{\kappa}|\rangle^{2})/T$. This order parameter
is able to capture emergent order induced by the second interaction term in
Eq. 1 and is alone the proper way to reveal an order-disorder phase transition
in a region where several ordered phases coexist.

Figure 4: Scaling behaviour of the susceptibility for $n=2$ symmetry at the
maximum frustration ($\kappa=-1$). The inset shows a FSS analysis to obtain
the transition temperature.

In this paper, we are not going to concentrate on the thermodynamics of the
phase transitions in the range of the $\kappa$ out of the frustration region.
We aim at the determination of the critical behavior only at the maximum
frustration. For the other values of $\kappa$ we calculate the above order
parameters only to obtain the full $T_{c}\times\kappa$ phase diagram for both
symmetries. To perform the finite-size scaling (FSS) analysis at the maximum
frustration, we used the outcome of the entropic simulations based on the WL
algorithm. We ran sampling for lattice sizes of $L=28$, 32, 36, 40, 44, 48,
52, 56, 64 e 72 with $N=24$, 24, 24, 20, 20, 20, 16, 16, 16, and 12
independent runs, respectively. The flatness criterion used was 80%. Every run
was stopped at the twentieth level of the WL, $f_{21}$, where the parameter
$\varepsilon\sim 10^{-4}$. The order parameter $m_{\kappa}$ was certain in
characterizing an order-disorder phase transition at the maximum frustration,
as can be seen by the behavior of the susceptibility $\chi_{\kappa}$ shown in
Fig. 4. To obtain the critical temperature, we first estimate the value of the
critical exponent $\nu=0.6101(16)$ using the quantities defined from the
derivatives of the logarithm of the order parameter $m_{\kappa}$[27]. Then, we
used the scaling law $T_{c}(L)=T_{c}+aL^{-\frac{1}{\nu}}$, where $T_{c}$ is
the critical temperature of the infinite system. In the inset of Fig. 4 we
show the expected FSS behavior for a second-order phase transition and the
linear regression yields $T_{c}=0.6441(25)$ for $n=2$. For $n=4$ at the
maximum frustration, since the system behaves as an ensemble of non-
interacting spins, the expected critical temperature is exactly $T_{c}=0.0$.

Figure 5: Phase diagram in the coupling-temperature plane for both symmetries.
Dashed and solid curves indicate first-order and continuous phase transitions,
respectively. The square symbols indicate the critical points after a FSS
analysis.

To construct the $T_{c}\times\kappa$ phase diagrams we estimate the joint
density of states $g(E_{1},E_{2})$, where $E_{1}$ and $E_{2}$ are,
respectively $E_{1}=\sum_{\langle ij\rangle}\sigma_{i}\sigma_{j}$ and
$E_{2}=\sum_{\langle ij\rangle}\sigma_{i}\sigma_{j}\cos n\phi_{ij}$[28].
Hence, the thermodynamic averages can be obtained for any value of the
parameters $T$, $J$, and $J_{a}$[24, 29]. For these cases we used the lattice
sizes $L=16$, 20, 24, 28, 32, 36 e 40. Due to the high computational cost, the
WL runs were done up to $f=f_{15}$. Fig. 5 shows phase diagrams for both $n=2$
(top) and $n=4$ (bottom). The solid line represents the temperature of the
maximum of the susceptibility for $L=32$. The symbols with the corresponding
error bars are the FSS analysis for different values of $\kappa$. For $n=2$,
the discontinuity of the solid line indicates the limit of validity of the
order parameters $m_{stp}$ as $\kappa\rightarrow(-1)^{-}$ and $m$ as
$\kappa\rightarrow(-1)^{+}$. The point at the maximum frustration $\kappa=-1$
was obtained by an FSS analysis of the $m_{\kappa}$ order parameter. In the
region of phase coexistence (dashed line in Fig. 5), the energies of the
ordered phases are too close and there is a huge number of configurations with
even closer energies. Hence, small thermal fluctuations is enough to allow low
energy configurations to be filled. That is the reason for critical
temperatures in the frustrated region be close to zero. For $n=4$, as expected
the lines of phase transition follow the evolution of the ground state energy
$E_{0}(\kappa)$ as a function of $\kappa$ and the transition temperature at
the maximum frustration occurs only at $T=0$. While the lines of phase
transitions shown are related to an order-disorder transition, it is also
expected an Ising-nematic phase between the stripes and disordered phases[13,
26]. However, an analysis of the nature of such transition is left to be dealt
with in future work.

In summary, we have proposed a simple spin lattice model that couples the spin
degree of freedom to the underlying lattice symmetry. Our model allows tuning
between different phases, magnetically ordered as well as more exotic
nonmagnetic phases by changing only one control parameter. Also, it allows the
manifestation of a spin-nematic state at the maximum frustration coexisting
with two other ordered phases. In the present form, our model may offer a
generic explanation as to why exotic spin patterns are so frequently observed
in various natural magnetic systems. Besides, our results motivate further
systematic studies including spin systems with 2D (XY) or 3D (Heisenberg)
degrees of freedom, other lattice geometries such as honeycomb, triangular, or
Kagome as well as systems with long-range interactions.

We thank the computer support from LaMCAD/UFG.

## References

  * Moessner and Ramirez [2006] R. Moessner and A. P. Ramirez, Geometrical frustration, Physics Today 59, 24 (2006).
  * Lacroix  _et al._ [2011] C. Lacroix, F. Mendels, and F. Mila, eds.,  _Introduction to Frustrated Magnetism_ (Springer-Verlag, 2011).
  * Kohama  _et al._ [2019] Y. Kohama, H. Ishikawa, A. Matsuo, K. Kindo, N. Shannon, and Z. Hiroi, Possible observation of quantum spin-nematic phase in a frustrated magnet, Proc Natl Acad Sci USA 116, 10686 (2019).
  * Klein  _et al._ [2019] D. R. Klein, D. MacNeill, Q. Song, D. T. Larson, S. Fang, M. Xu, R. A. Ribeiro, P. C. Canfield, E. Kaxiras, R. Comin, and P. Jarillo-Herrero, Enhancement of interlayer exchange in an ultrathin two-dimensional magnet, Nature Physics 15, 1255 (2019).
  * Jin  _et al._ [2012] S. Jin, A. Sen, and A. W. Sandvik, Ashkin-teller criticality and pseudo-first-order behavior in a frustrated ising model on the square lattice, Phys. Rev. Lett. 108, 045702 (2012).
  * Sapkota  _et al._ [2017] A. Sapkota, B. G. Ueland, V. K. Anand, N. S. Sangeetha, D. L. Abernathy, M. B. Stone, J. L. Niedziela, D. C. Johnston, A. Kreyssig, A. I. Goldman, and R. J. McQueeney, Effective one-dimensional coupling in the highly frustrated square-lattice itinerant magnet caco2as2, Phys. Rev. Lett. 119, 147201 (2017).
  * Watanabe  _et al._ [2016] D. Watanabe, K. Sugii, M. Shimozawa, Y. Suzuki, T. Yajima, H. Ishikawa, Z. Hiroi, T. Shibauchi, Y. Matsuda, and M. Yamashita, Emergence of nontrivial magnetic excitations in a spin-liquid state of kagomé volborthite, Proceedings of the National Academy of Sciences 113, 8653 (2016).
  * Fert  _et al._ [2017] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nat. Rev. Mater. 2, 17031 (2017).
  * Sohn  _et al._ [2017] C. H. Sohn, C. H. Kim, L. J. Sandilands, N. T. M. Hien, S. Y. Kim, H. J. Park, K. W. Kim, S. J. Moon, J. Yamaura, Z. Hiroi, and T. W. Noh, Strong spin-phonon coupling mediated by single ion anisotropy in the all-in$--$all-out pyrochlore magnet cd2os2o7, Phys. Rev. Lett. 118, 117201 (2017).
  * Shinaoka  _et al._ [2012] H. Shinaoka, T. Miyake, and S. Ishibashi, Noncollinear magnetism and spin-orbit coupling in $5d$ pyrochlore oxide cd2os2o7, Phys. Rev. Lett. 108, 247204 (2012).
  * Wang and Vishwanath [2008] F. Wang and A. Vishwanath, Spin phonon induced collinear order and magnetization plateaus in triangular and kagome antiferromagnets: Applications to cufeo2, Phys. Rev. Lett. 100, 077201 (2008).
  * Shannon  _et al._ [2006] N. Shannon, T. Momoi, and P. Sindzingre, Nematic order in square lattice frustrated ferromagnets, Phys. Rev. Lett. 96, 027213 (2006).
  * Guerrero  _et al._ [2015] A. I. Guerrero, D. A. Stariolo, and N. G. Almarza, Nematic phase in the $j_{1}-j_{2}$ square-lattice ising model in an external field, Phys. Rev. E 91, 052123 (2015).
  * Bergman  _et al._ [2006] D. L. Bergman, R. Shindou, G. A. Fiete, and L. Balents, Models of degeneracy breaking in pyrochlore antiferromagnets, Phys. Rev. B 74, 134409 (2006).
  * Edlund and Jacobi [2010] E. Edlund and M. N. Jacobi, Universality of striped morphologies, Phys. Rev. Lett. 105, 137203 (2010).
  * Sherman [1995] E. Y. Sherman, Angular pinning and domain structure of a two-dimensional wigner crystal in a iii-v semiconductor, Phys. Rev. B 52, 1512 (1995).
  * DaSilva  _et al._ [2010] C. J. DaSilva, J. P. Rino, P. F. Damasceno, and E. Y. Sherman, Two-dimensional coulomb solid with interaction anisotropy, Phys. Rev. B 81, 153307 (2010).
  * Sabiryanov and Jaswal [1999] R. F. Sabiryanov and S. S. Jaswal, Magnons and magnon-phonon interactions in iron, Phys. Rev. Lett. 83, 2062 (1999).
  * Yin  _et al._ [2013] J. Yin, M. Eisenbach, D. M. Nicholson, and A. Rusanu, Effect of longitudinal degree of freedom of magnetic moment in body-centered-cubic iron, Journal of Applied Physics 113, 17E112 (2013), https://doi.org/10.1063/1.4794136 .
  * Wang and Landau [2001] F. Wang and D. P. Landau, Determining the density of states for classical statistical models: A random walk algorithm to produce a flat histogram, Phys. Rev. E 64, 056101 (2001).
  * Caparica and Cunha-Netto [2012] A. A. Caparica and A. G. Cunha-Netto, Wang-landau sampling: Improving accuracy, Phys. Rev. E 85, 046702 (2012).
  * Caparica [2014] A. A. Caparica, Wang-landau sampling: A criterion for halting the simulations, Phys. Rev. E 89, 043301 (2014).
  * Ferreira  _et al._ [2018] L. Ferreira, L. Jorge, S. L. ao, and A. Caparica, Wang-landau sampling: Saving cpu time, Journal of Computational Physics 358, 130 (2018).
  * Ferreira  _et al._ [2021] L. Ferreira, L. Jorge, C. J. DaSilva, and A. A. Caparica, Thorough analysis of the phase diagram for the bell–lavis model: An entropic simulational study, Braz. J. Phys. 51, 532–543 (2021).
  * Cannas  _et al._ [2006] S. A. Cannas, M. F. Michelon, D. A. Stariolo, and F. A. Tamarit, Ising nematic phase in ultrathin magnetic films: A monte carlo study, Phys. Rev. B 73, 184425 (2006).
  * Cannas  _et al._ [2004] S. A. Cannas, D. A. Stariolo, and F. A. Tamarit, Stripe-tetragonal first-order phase transition in ultrafilm magnetic films, Phys. Rev. B 69, 092409 (2004).
  * Caparica  _et al._ [2000] A. A. Caparica, A. Bunker, and D. P. Landau, Classical ferromagnet with double-exchange interaction: High-resolution monte carlo simulations, Phys. Rev. B 62, 9458 (2000).
  * Ferreira  _et al._ [2019] L. S. Ferreira, Álvaro A. Caparica, L. N. Jorge, and M. A. Neto, Thermodynamic properties of interacting like-rod chains: Entropic sampling simulations, Chemical Physics 517, 119 (2019).
  * Jorge  _et al._ [2021] L. Jorge, P. Martins, C. DaSilva, L. Ferreira, and A. Caparica, An entropic simulational study of the spin-1 baxter–wu model in a crystal field, Physica A: Statistical Mechanics and its Applications 576, 126071 (2021).


# Real-time Drift Detection on Time-series Data

Nandini Ramanan Rasool Tahmasbi Marjorie Sayer
Deokwoo Jung Shalini Hemachandran
Claudionor Nunes Coelho Jr

Advanced Applied AI Research
Palo Alto Networks
https://www.paloaltonetworks.com

###### Abstract

Practical machine learning applications involving time series data, such as
firewall log analysis to proactively detect anomalous behavior, are concerned
with real time analysis of streaming data. Consequently, we need to update the
ML models as the statistical characteristics of such data may shift frequently
with time. One alternative explored in the literature is to retrain models
with updated data whenever the model’s accuracy is observed to degrade.
However, these methods rely on near real-time availability of ground truth,
which is rarely fulfilled. Further, in applications with seasonal data,
temporal concept drift is confounded by seasonal variation. In this work we
propose an approach called Unsupervised Temporal Drift Detector or UTDD to
flexibly account for seasonal variation, efficiently detect temporal concept
drift in time series data in the absence of ground truth, and subsequently
adapt our ML models to concept drift for better generalization.

##  1 Introduction

In statistical machine learning, drift detection is defined as the detection
of the deterioration of model performance, requiring users to retrain the
model based on new inputs. As we move from a model-centric to a data-centric
approach, users are more and more using pre-defined pipelines with large
amounts of data. In this scenario, we expect model evaluation to play an
important role in determining when model performance deteriorates, as this is
ultimately related to training cost in the cloud. Usually, this deterioration
is detected at the input level, as a deviation from input parameters used at
training time is the first indication that previous model training may no
longer be valid.

We are particularly interested in detecting model drift in temporal series
(one where a portion of the inputs refer to previous values of a variable we
want to forecast or to detect anomalies), where the machine learning model is
used to estimate a possibly non-linear function of previous $p$ instances,
such as given by Formula 1, where $y_{i}$ and $x_{j}$ are referred usually as
an endogenous variable and an exogenous variable respectively.

$y_{t}=f(y_{t-p:t-1},x_{t-p:t})$ (1)

When doing predictions or detecting anomalies in time series, one needs to
consider seasonality and cycles, which may present themselves in a non-linear
relation, and where their relation may not be known in a general case, such as
in the example presented in Figure 1. In this case, in order to detect drift,
we need to consider the unknown effects of the seasonality and cycles, as it
can be easily shown in this figure that the model contains two behaviors - one
during the week, and another one during the weekends.

This paper proposes an algorithm Unsupervised Temporal Drift Detector or
interchangeably referred as UTDD to address the problem of modeling
seasonalities and cycles in a generic way so that we can detect the drift in
time series variables by detecting anomalies in their behavior, after we
remove the seasonalities and cycles. We consider a generic framework that
accommodates unknown relations and non-linear relations as well.

Figure 1: Synthetic time-series data with seasonality/cycles between the month
of Aug to Oct 2020

This paper is organized as follows. In the next section, we present related
work to drift detection, including common techniques used in drift detection.
Then, we present our work to detect model deterioration based on detecting
anomalies in previous variables once we remove seasonalities and cycles in a
black-box way. We do this using artifacts of Boosted Embeddings. Then, we
present results and conclusions.

##  2 Related Work

Near-limitless volumes of data in temporal order, which are referred to as
data streams, are generally nonstationary as the characteristics of data
evolve over time. This phenomenon is called concept drift, and is an issue of
great importance in the literature, because it makes models obsolete by
decreasing their predictive performance. In the presence of concept drift, it
is necessary to adapt to changes in data to build more robust and effective
classifiers.

Gemaque [3] presented a comprehensive overview of approaches that tackle
concept drift in classification problems in an unsupervised manner. Brockhoff
et.al [1] used Earth Mover’s Distance, a measure of the distance between two
probability distributions, to detect the drift. In 2020, Liu [7] proposed a
cluster-based histogram, called equal intensity k-means space partitioning
(EI-kMeans). Pinagé [8] proposed a semi-supervised drift detector that uses an
ensemble of classifiers based on self-training online learning and dynamic
classifier selection. Xuan [9] presented a Bayesian non-parametric
unsupervised concept drift detection method based on the Poly tree hypothesis
test. The basic idea is to decompose the underlying data distribution into a
multi-resolution representation that transforms the whole distribution
hypothesis test into recursive and simple binomial tests.

On the one hand, a high accuracy detection approach usually requires labeled
data, possibly involving high cost for labeling. On the other hand, a variety
of methods have been devoted to the topic of concept drift detection with
unlabeled data, but these approaches often are most suited for only a subset
of the concept drift types. Hu et. al. [5] did a survey to present these
methods, categorize them and give recommendations of usage based on their
behaviors under different types of concept drift. Recently, Gözüaçık and Can
[4] presented an implicit (unsupervised) algorithm called One-Class Drift
Detector (OCDD), which uses a one-class learner with a sliding window to
detect concept drift.

##  3 Temporal Drift Detection

### Concept Drift in Time-series Data

Concept drift is a phenomenon in which the statistical properties of a domain
change over time in an arbitrary way. Given a set of $k$ samples from a joint
distribution over a time interval $[t_{1},t_{k}]$,

$\displaystyle\mathcal{S}=\\{\mathcal{D}_{t_{1}},\ldots,\mathcal{D}_{t_{k}}\\}\sim
F_{[t_{1},t_{k}]}(X,y)$ (2)

where dataset $\mathcal{D}_{t_{i}}$ is comprised of both a feature vector
$X_{t_{i}}$ and label $y_{t_{i}}$ such that
$\mathcal{D}_{t_{i}}=(X_{t_{i}},y_{t_{i}})$ and the notation
$F_{[t_{1},t_{k}]}$ making the connection between a joint distribution $F$ and
the time interval $[t_{1},t_{k}]$ explicit. Concept drift can then be defined
as,

$\displaystyle F_{[t_{1},t_{k}]}(X,y)\neq F_{t_{[k+1,\infty]}}(X,y),\text{
\,\,or equivalently,}$ (3) $\displaystyle\exists t_{j}:Pr_{t_{j}}(X,y)\neq
Pr_{t_{j}+1}(X,y)$ (4)

In this paper we study unsupervised model adaptation under concept drift. To
be able to detect drift in the underlying concept based on the variables X
alone, we must observe some drift or changes in the feature distribution
$F_{[t_{1},t_{k}]}{[X]}$ with time.

### Boosted Embeddings

In this section we leverage embeddings to learn seasonality (e.g., daily,
weekly, monthly) or unknown cycles by multiple categorical features. Let’s
assume $\theta^{T}$ to capture time-categorical features (e.g., months of the
year, days of the week, and hours of the day) and $\theta^{I}$ to represent
other independent categorical features. Let assume that a model trained is
trained at at $t_{1},t_{2},...t_{m}$ where $t_{m-1}<t_{m}$. Then, the model
trained at time $t_{m}$ can be formulated as
$f^{m}:=f^{m}_{emb}(\mathbf{x};\theta_{m}^{T},\theta_{m}^{I})+f^{m}_{res}(\mathbf{x})$
where $f_{emb}$ is embedding model and $f_{res}$ is a residual model. If the
model drift occurs at time $t_{m+1}$, it will incur a large
$f^{m}_{res}(\mathbf{x})$ that will be captured by our UTDD algorithm. We
leveraged the DeepGB algorithm in [6], where gradient boosting trains several
simple models sequentially. The key idea of boosting is that each subsequent
model trains only on the difference of the output and previous model, to
leverage each model’s strengths and minimize regression error. Their approach
is conducting gradient boosting to fit weak learners on residuals to improve
the previous models. They propose a loop wherein, at each iteration, they
freeze the previous embedding, and add the new embedding to the sequence of
models. Such that, $f^{m}=[e^{m}_{1},\ldots,e^{m}_{L},r^{m}].$ where
$e_{i},i=1,\cdots,L$ are embedding models to capture categorical data, and $r$
represents the residual model. The summary of their boosted embedding approach
is presented in Algorithm 1.

1:function BoostedEmbeddings($\mathbf{X}$,$N$)

2: where $\mathbf{X}$ = $((t_{n},\mathbf{y}_{n}))^{N}_{n=1}$

3: $f^{m}=[]$

4: $F_{0}:=y$

5: for $1\leq l\leq L+1$: do $\triangleright$ iteration over the embedding
models

6: $e_{l}.fit(t,F_{l-1})$ $\triangleright$ fitting the selected embedding
model

7: $F_{l}=F_{l-1}-e_{l}.predict(t)$ $\triangleright$ residual computation

8: if $|F_{l}-F_{l-1}|<\epsilon:$ then $\triangleright$ Check termination
condition

9: break

10: end if

11: $f^{m}.append(e_{l})$

12: end for

13: return $f^{m}$

14:end function

Algorithm 1 Boosted Embeddings Algorithm

### Drift Detection with Boosted Embeddings

The proposed UTDD algorithm 2 is an unsupervised temporal data drift detection
algorithm where we followed a sequence of steps to deseasonalize the time-
series data. In line[3], we estimate the number of differences required to
make the given time series stationary using Augmented Dickey–Fuller Functions.
We observe from our experiments that the $max.diff=4$ is accurate for our
domain. In line[4], we apply the estimated differencing.

1:function UnsupervisedTemporalDriftDetector($\mathbf{X}$,$N$)

2: where $\mathbf{X}$ = $[X_{1},X_{2},..,X_{n=N}]$ such that
$X_{n}\in\mathbb{R}$, $X_{n}$ is the datapoint at time $t$

3: $k=\textsc{ndiffs}(x,max.diff=4,test=``adf")$ $\triangleright$ get min # of
differences estimates required to make a time series stationary

4: $X$ = $X.diff(k)$

5: $P$ = $\textsc{BoostedEmbeddings}(\mathbf{X},N)$ $\triangleright$ Fit BE to
seasonality + additional categorical variables

6: $Error=X-P$ $\triangleright$ Compute residual

7: $ZScoreCurr=\textsc{ComputeZScore}(Error)$

8: return $ZScoreCurr$

9:end function

Algorithm 2 Unsupervised Temporal Drift Detector

In line[5], we employ BE to fit to the seasonality along with any other
additional categorical (that can be from domain expertise) in the resultant
data. The BE model enables our approach to learn from the distinct time-series
signature at once by encoding the categorical features (linear time, weekday
or weekend effect, holiday effect, user-defined exogenous categorical
variables etc) in a meaningful way, in a lower dimensional space to extract
valuable information. We compute the residual, $Error$, by removing the
seasonal component from the raw data $X$. Finally, we compute the z_score of
this residual, otherwise referred to as resultant stationary series, and
evaluate if it has changed from the $z\\_score$ data on which we had trained
our model previously. Consequently, our approach is the first to do temporal
drift detection in the context of unsupervised time series data, enabling
automated cost-efficient retraining when the data has drifted.

##  4 Experimental Evaluation

In order to simulate time series data with multiple seasonal components with
different periodicities, let’s denote the time series by $x_{t}$, and the
seasonal component $i$ by $\gamma_{t}^{(i)}$. For simplicity assume we just
have two seasonal components:

$x_{t}=\mu_{t}+\gamma_{t}^{(1)}+\gamma_{t}^{(2)}+\epsilon_{t}$

where $\mu_{t}$ represents the trend or level, $\gamma_{t}^{(1)}$ represents a
seasonal component with a relatively short period, and $\gamma_{t}^{(2)}$
represents another seasonal component of longer period and $\epsilon_{t}$ is
white noise. Following equations (3.7) and (3.8) in Durbin and Koopman [2], we
can simulate the seasonality component $i$ with seasonal length of $s$ by

$\displaystyle\gamma^{(i)}_{t}$
$\displaystyle=\sum_{j=1}^{p}\gamma^{(i)}_{j,t}$ (5)
$\displaystyle\gamma^{(i)}_{j,t+1}$
$\displaystyle=\gamma^{(i)}_{j,t}\cos(\lambda_{j})+\gamma^{*,(i)}_{j,t}\sin(\lambda_{j})+\omega^{(i)}_{j,t},$
(6) $\displaystyle\gamma^{*,(i)}_{j,t+1}$
$\displaystyle=-\gamma^{(i)}_{j,t}\sin(\lambda_{j})+\gamma^{*,(i)}_{j,t}\cos(\lambda_{j})+\omega^{*,(i)}_{j,t},$
(7)

where $p=\lfloor s/2\rfloor$ is the periodicity, $1\leq j\leq p$,
$\lambda_{j}=\frac{2\pi j}{s}$, $\omega^{(i)}_{j,t}\sim
N(0,\sigma^{2}_{\omega})$ and $\omega^{*,(i)}_{j,t}\sim
N(0,\sigma^{2}_{\omega^{*}})$.

We generated artificial data with seasonal patterns that differ in their peaks
and shapes during weekdays and weekends as well as holidays as shown in Fig 1
from August through October. We estimate the drift between the data in August
to September, and September to October respectively. To this effect, Fig 5 is
the result of fitting BE to the seasonality in data between Aug to Sep. Fig 5
is the residual as a result of applying line[6] of algorithm 2. We repeat the
same with the data in the range between Sep to Oct, such that, Fig 5 after
fitting BE to the seasonality and Fig 5 is the residual. Such that the
resultant $z\\_score$ of $0.53$ and $0.65$ for Fig 5 & Fig 5 respectively show
that the data has drifted, as evident from the plots.

Figure 2: Fitted BE to the seasonality to the data between Aug to Sep.
Seasonality in Red. Residual in Blue.

Figure 3: Residual in Yellow for Aug to Sep. Residual Series will be used to
compute z score.

Figure 4: Fitted BE to the seasonality to the data between Sep to Oct.
Seasonality in Red. Residual in Blue.

Figure 5: Residual in Yellow for Sep to Oct. Residual Series will be used to
compute z score.

##  5 Conclusion

The literature for concept drift detection mostly relies on discovering
significant drop in model generalization. In this paper, we follow the
hypothesis that it is too unrealistic to assume that the labels are readily
available as with the case with many real time problems. Therefore, we tackle
the concept drift detection problem in an unsupervised manner. Unsupervised
temporal drift detection methods are widely used in a variety of research
areas as well as practical application domains. In this paper, we propose a
novel label-free drift detection algorithm, UTDD, for time series data. In our
knowledge, this would be the first work to do temporal drift detection in the
context of unsupervised time series data by subtracting the effect of
seasonality and modeling for the residual. Future work will focus on extending
this preliminary study to evaluate and compare highlighted methods in related
section by carrying out extensive experiments.

## References

  * [1] T. Brockhoff, M. S. Uysal, and W. M. van der Aalst.  Time-aware concept drift detection using the earth mover’s distance.  In 2020 2nd International Conference on Process Mining (ICPM), pages 33–40. IEEE, 2020.
  * [2] J. Durbin and S. J. Koopman.  Time series analysis by state space methods.  Oxford university press, 2012.
  * [3] R. N. Gemaque, A. F. J. Costa, R. Giusti, and E. M. dos Santos.  An overview of unsupervised drift detection methods.  WIREs Data Mining and Knowledge Discovery, 10(6):e1381, 2020.
  * [4] Ö. Gözüaçık and F. Can.  Concept learning using one-class classifiers for implicit drift detection in evolving data streams.  Artificial Intelligence Review, 54(5):3725–3747, 2021.
  * [5] H. Hu, M. Kantardzic, and T. S. Sethi.  No free lunch theorem for concept drift detection in streaming data classification: A review.  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery, 10(2):e1327, 2020.
  * [6] S. R. Karingula, N. Ramanan, R. Tahsambi, M. Amjadi, D. Jung, R. Si, C. Thimmisetty, and C. N. Coelho Jr.  Boosted embeddings for time series forecasting.  arXiv preprint arXiv:2104.04781, 2021.
  * [7] A. Liu, J. Lu, and G. Zhang.  Concept drift detection via equal intensity k-means space partitioning.  IEEE transactions on cybernetics, 51(6):3198–3211, 2020.
  * [8] F. Pinagé, E. M. dos Santos, and J. Gama.  A drift detection method based on dynamic classifier selection.  Data Mining and Knowledge Discovery, 34(1):50–74, 2020.
  * [9] J. Xuan, J. Lu, and G. Zhang.  Bayesian nonparametric unsupervised concept drift detection for data stream mining.  ACM Transactions on Intelligent Systems and Technology (TIST), 12(1):1–22, 2020.