11institutetext: University of Camerino, Camerino MC 62031, Italy
11email: {flavio.corradini, michele.loreti<EMAIL_ADDRESS>
<EMAIL_ADDRESS>

# REPTILE: A Proactive Real-Time Deep Reinforcement Learning Self-adaptive
Framework

Flavio Corradini 11 0000-1111-2222-3333 Michele Loreti 11 1111-2222-3333-4444
Marco Piangerelli 11 2222–3333-4444-5555 Giacomo Rocchetti 11

###### Abstract

In this work a general framework is proposed to support the development of
software systems that are able to adapt their behaviour according to the
operating environment changes. The proposed approach, named REPTILE, works in
a complete proactive manner and relies on Deep Reinforcement Learning-based
agents to react to events, referred as novelties, that can affect the expected
behaviour of the system. In our framework, two types of novelties are taken
into account: those related to the context/environment and those related to
the physical architecture itself. The framework, predicting those novelties
before their occurrence, extracts time-changing models of the environment and
uses a suitable Markov Decision Process to deal with the real-time setting.
Moreover, the architecture of our RL agent evolves based on the possible
actions that can be taken.

###### Keywords:

Self-Adaptivity Deep Reinforcement Learning Proactivity

##  1 Intoduction

The last decades witnessed a significant evolution in software designing
enactment as the process of managing, developing and maintaining software is
becoming more and more complex. The always increasing expansion of Internet
together with the advent of a plethora of _smart devices_ , posed new
challenges to engineers and developers that have to design and maintain
software infrastructure that must be able to _adapt their behaviour_ to new
requirements and updates in the system architecture often unknown at design
time. These can have a great impact on the Quality of Service (QoS) offered by
the system, thus impacting the final user experience.

In this context, Kephart and Chess stated that “ _the only option remaining is
autonomic computing systems that can manage themselves given high-level
objectives from administrators_ ” [10]. Such systems are commonly called Self-
Adaptive Systems (SASs) and they have been extensively researched and
developed in the last twenty years. The main feature of SASs is to be able to
operate in an environment that is changing in an unpredictable way without
any, or limited, external control. For this reason, a big effort has made to
design methodologies and techniques supporting SASs adaptation.

Several approaches have been proposed to design and develop SASs [13]. In
particular we can mention here the _model-based_ , the _architecture-based_ ,
the _control theory-based_ , the _formal modeling and verification-approach_
and the Learning-based approaches. The _model-based_ approach uses model-
driven engineering, filling incomplete information at design time with the
usage of _run-time models_ (MUSIC framework [9]). The _architecture-based_
approach use the concept of the high level structure of software for different
activities. One of the most well-known frameworks that implements this
approach is the _Rainbow framework_ , which provides an architecture layer
whose components define adaptation plans [8]. In the _control theory-based_
approach the system continuously takes measurements (from sensors) and
performs adjustments to keep the measured variable in a defined range in a
process called _feedback loop_ [18]. Together with the implementation of SASs,
there is the need to prove the correctness of such systems using _formal
modeling and verification_ approaches. A common method is the _formal
reference model for self-adaptation_ (FORMS), which uses the Z notation in
order to define the whole system, from its interaction with the environment to
the adaptation [27]. The latter approach, the learning-based one, given the
increasing popularity of the field of _Machine Learning_ (ML), mainly in the
last decade, affected a lot the research of Self-Adaptive Systems using ML
algorithms for a variety of purposes: from the use of genetic algorithms for
the run-time adaptation of mobile applications in [20], or the use of
_Reinforcement Learning_ (RL) for an architecture-based adaptation [11], to
the application of _Deep Learning_ (DL) for SASs in network security [19, 15].

Moreover, challenges such as the analysis of streams of data in real-time and
the extraction of environmental models in order to aid a proactive
adaptability for unforeseen changes have to be adequately taken into account.
To this extent, we can mention a framework for self adaptive systems proposed
in [7]. The authors tackle the problem of unforeseen changes in the
environment by learning the impact of adaptation decisions on the system’s
goal applying a reactive, i.e. post-hoc or after-the-fact, model. A more
refined framework is the one proposed in [12]. It is based on the MAPE-K
framework, extended with a simplified RL component.

Finally, the IoTArchML framework refers to a proactive SAS applied to _iot_
systems [17]. In particular, it is designed for the control of system’s QoS
parameters, using a proactive approach to anticipate eventual changes before a
QoS deviation event.

In this work we propose a general _framework for software SASs dealing with
real-time streams of data_. Our framework works in a complete proactive
manner. It uses of _Deep Reinforcement Learning_ (DRL) algorithms, in order to
adapt itself to eventual _novelties_ , i.e. new events that can be managed in
order to work properly. In our framework, as already pointed out in [12], we
identify two types of novelties: those relating to the context/environment and
those relating to the physical architecture itself. The framework, predicting
those novelties before their occurrence, extract time-changing models of the
environment and uses a modified Markov Decision Process to deal with the real-
time setting as shown in [23]. Moreover the architecture of our RL agent
evolves based on the possible actions that can be taken.

The paper is structured as follows: in section 2 the self-adaptability
framework is described, section 3 presents the explanation of the adaptability
mechanisms, in section 4 we show the results obtained by applying the
framework to a real case study and, finally, in section 5 we discuss the
results.

##  2 REPTILE Architecture

In this Section we provide a general overview of our framework, REPTILE. It is
based on the well known IBM’s MAPE-K architecture enriched with a
_Reinforcement Deep Learning module_. The whole architecture of REPTILE is
depicted in Figure 1. It consists of six modules: the five from the standard
MAPE-K architecture (_Monitoring_ , _Analyzing_ , _Planning_ , _Executing_ ,
and _Knowledge_) plus a _Reinforcement Deep Learning_ one. Below, the role of
these modules is outlined together with their use in a simple running
scenario.

Figure 1: Overview of REPTILE.

Monitoring (yellow). It continuously collects data from the environment and
performs an extraction of useful features of the system, based on the domain
application. It then sends data to the Knowledge component.

######  Example 1 (Running Case Study (1/6))

To better illustrate the technicalities of our framework we present here a
running case study. We applied our framework to a Heating, Ventilation and Air
Conditioning (HVAC) system. The reason behind the choice of this particular
application domain is its simplicity. A simulator for the environment was
implemented in Python: it computes the next internal temperature of a single
room based on the actions performed by the system. The simulator is made up by
the following three devices:

  * •

Heater: it produces power (Watts) in order to heat the room.

  * •

Cooler: it produces power (Watts) in order to cool the room.

  * •

Window: managed by a binary signal (open or closed).

The gained/lost temperature inside a room is obtained by calculating the heat
loss between the room and the outside, through walls and windows, and
converting the power produced by the heater as a heat gain and the one
produced by the cooler as heat loss.

Both the description of the room and the available devices are customizable. A
device is any instrument that can modify the internal temperature with its
action. Table 1 shows the simulations parameters. External temperatures data
have been pre-processes from two existing datasets [1, 2].

__

Parameter | Value
---|---
Room Width | 5
Room Height | 3
Room Depth | 3
Wall R-Value | 3.6
Window R-Value | 0.176
Air Changes per Hour | 1.7
N. Heaters | 1
N. Coolers | 1
N. Windows | 1
Heater Power (W) | 0, 200, 400
Cooler Power (W) | 0, 400
Table 1: Room parameters used in simulations

Knowledge (yellow and blue). It receives data from the Monitoring module and
it builds two types of models of the environment: _time-varying models_ ,
$\mathcal{M}_{t,{\bar{\alpha}}}^{\mathcal{TV}}$. and _global models_ ,
$\mathcal{M}_{\xi}^{\mathcal{G}}$. The former is used for monitoring the model
in the short-period, while the latter represents the comprehensive knowledge
of the environment. In the previous notation $t$ represents a generic time
step, $\bar{\alpha}$ the vector containing the parameters of the model.
Superscript $\mathcal{TV}$ indicates that we are modeling the time-varying,
while $\mathcal{G}$ stands for Global, where $\xi=1,2,3...,k$ is the index
identifying each model . Moreover, this module contains some pre-trained
agents, $\mathcal{M}_{\xi}^{\mathcal{A}}$, where $\xi=1,2,3...,n$ is the index
identifying each model and the superscript $\mathcal{A}$ stands for agent.

######  Example 2 (Running Case Study (2/6))

In our running case study, the _knowledge_ consists in a set of _models_ (
$\mathcal{M}_{t,{\bar{\alpha}}}^{\mathcal{TV}}$,
$\mathcal{M}_{\xi}^{\mathcal{G}}$, $\mathcal{M}_{\xi}^{\mathcal{A}}$). Each
$\mathcal{M}_{t,{\bar{\alpha}}}^{\mathcal{TV}}$ is represented by using an
Auto-Regressive, Moving Average and Integrated (ARIMA) model [3].

Such models are described formally by the following equation:

$y^{\prime}_{t}=c+\epsilon_{t}+\sum_{i=1}^{p}\varphi_{i}y^{\prime}_{t-i}+\sum_{i=1}^{q}\theta_{i}\epsilon_{t-i}$
(1)

where $\theta_{i}$ and $\phi_{i}$ are the parameters of the model, $p$ and $q$
are the orders, respectively of the auto-regression and the moving average
contributions; The orders of an ARIMA model are defined as follows:

  * •

$p$: the number of lag observations included in the model, also called the lag
order.

  * •

$q$: the size of the moving average window, also called the order of moving
average.

$\epsilon_{t}\overset{iid}{\sim}N(0,\sigma^{2})$ and $c=\mu+\varphi_{0}$ with
$\mu$ the expected value of $y_{t}$ and $y_{t}^{\prime}=y_{t}-y_{t-1}$. To
obtain the model we used a number of data points $d=30$. Each time a new data
point arrives, we compute a new model discarding the oldest point end
incorporating the new one. The $\mathcal{M}_{\xi}^{\mathcal{G}}$ models can be
obtained using more sophisticated techniques. In this work, for simplicity, we
already had such models. Finally, $\mathcal{M}_{\xi}^{\mathcal{A}}$ is the
model obtained using some DRL-based model.

Analyzing (blue). It inspects the data stored in the Knowledge, that is a
_model_ of the collected data, and uses it to forecast possible evolution. At
each step this module is able to detect potential unseen behaviours, i.e.
novelties.

We can distinguish two kinds of novelties: Contextual novelty and
Architectural novelty [12]. The former indicates a _novelties_ are mainly
inferred by comparing the predictions done at the previous steps with the one
at the current step. For instance, the module can detect that an observable,
namely the _temperature_ , is going to run out of the expected range in the
next fifteen minutes, since the predicted value increased in the last
observations. The latter refer to a change in the interaction between the
system and the environment due to a new capability that is now available. For
instance, a new device is installed that can affect the observations.

Actually, this module is also able to detect novelties in a _reactive_ manner,
i.e. if the observable violates a certain ”soft” constraint for a given amount
of time. For instance, whenever _temperature_ is in the Intermediate zone (see
Section 2.1) for $t$ previous time steps.

######  Example 3 (Running Case Study (3/6))

In this case study, we identify novelties in the following way:

  * •

a new device has been inserted into the system (architectural novelty).

  * •

the internal temperature does not stay in the interval [$18$ °C, $22$ °C]
(contextual novelty).

First of all, REPTILE looks for new devices or changes in the devices
themselves, i.e. e new heater with different powers of use. If such a
component is detected, the system begins the adaptation procedure in order to
adapt the structure of the _learning_ module to manage new actions. If no new
device is found the system continues to behave as usual. Moreover, REPTILE
uses also the predictions of the next three points for deciding about the
onset of the adaptability process: if the third point (the farthest in time)
is not in the comfort zone then the adaptability process is triggered,
otherwise everything remains the same (_proactive_ process). In the meantime,
the $T=12$ previous time steps are monitored in order to establish if system
has been keeping violating constraints since 12 time steps: in that case the
adaptability process is triggered, too (_reactive_ process). The _executing_
module, i.e, the DRL-based agent performs an action every 15 minutes while the
ARIMA model is derived every 5 minutes meaning that the system had 3 time
steps to try to correct its behavior.

Learning (blue and grey). It is the crucial component of the framework and it
is activated when a novelty is foreseen. This module is responsible for the
choice of the optimal adaptive policy. It uses the Knowledge and a DRL-based
agent (the same used by the Planning module) for the adaptation process. Once
a novelty is detected, the Learning component is invoked, otherwise the Plan
component will decide the next action based on the current DRL-based agent.
Generally speaking, DRL is about maximizing a reward function, i.e. a cost
function, in order to make the DRL-based agent to learn how to behave in an
environment (refer to Subsection 2.1 for the reward function used in this
work).

######  Example 4 (Running Case Study (4/6))

We used the Deep Q-Learning algorithm for implementing the _Learning_ module:
in particular we implemented a Duelling Deep Q-network [26]. The Dueling Deep
Q-network splits the single stream of fully connected layers in two separate
streams, one for the value function, $V$, and one for the advantage function,
$A$, defined as:

$A^{\pi}(s,a)=Q^{\pi}(s,a)-V^{\pi}(s)$ (2)

In our implementation of the Dueling Deep-Q network we used the real time
markov decision process (RTMDP) framework [23], see Section 2.1 for more
details. Accordingly, states and actions evolve simultaneously. The state
vector is then denoted as a state-action pair
$\boldsymbol{x}_{t}=(s_{t},a_{t})$. The actions in input to the Duelling Deep
Q-network are encoded using a one-hot encoding of the combinations of the
basic actions of the three devices. For example, referring to Table 4, there
are $12$ actions (3 for the Heater, 2 for the Cooler and 2 for the Window).
Using a possible one-hot encoding, Table 2 shows all the possible actions:

Heater | Cooler | Window | $a_{t}$
---|---|---|---
0 | 0 | CLOSE | 1,0,0,0,0,0,0,0,0,0,0,0
0 | 0 | OPEN | 0,1,0,0,0,0,0,0,0,0,0,0
0 | 200 | CLOSE | 0,0,1,0,0,0,0,0,0,0,0,0
0 | 200 | OPEN | 0,0,0,1,0,0,0,0,0,0,0,0
200 | 0 | CLOSE | 0,0,0,0,1,0,0,0,0,0,0,0
200 | 0 | OPEN | 0,0,0,0,0,1,0,0,0,0,0,0
200 | 200 | CLOSE | 0,0,0,0,0,0,1,0,0,0,0,0
200 | 200 | OPEN | 0,0,0,0,0,0,0,1,0,0,0,0
400 | 0 | CLOSE | 0,0,0,0,0,0,0,0,1,0,0,0
400 | 0 | OPEN | 0,0,0,0,0,0,0,0,0,1,0,0
400 | 200 | CLOSE | 0,0,0,0,0,0,0,0,0,0,1,0
400 | 200 | OPEN | 0,0,0,0,0,0,0,0,0,0,0,1
Table 2: One-Hot encodind of the actions $a_{t}$

If the vector $a_{t}=(1,0,0,0,0,0,0,0,0,0,0,0)$ meaning Heater = 0 W, Cooler =
0 W, Window = CLOSE, and $s_{t}=(T_{ext}=17,T_{in}=19)$, then
$\textbf{x}_{t}=(17,19,1,0,0,0,0,0,0,0,0,0,0,0)$. As shown in Figure 2, our
Duelling Deep-Q Networks is made by 3 fully connected hidden layers each of
which contains 256 neurons. The Input layer contains 14 neurons and the Output
ones consist of 12 neurons.

Figure 2: The Dueling Deep Q-Network used in this work.

Planning (green). This component simply uses the DRL-based agent to choose the
next action to perform.

######  Example 5 (Running Case Study (5/6))

Given in input $\textbf{x}_{t}=(17,19,1,0,0,0,0,\\\ 0,0,0,0,0,0,0)$, then the
Deep Q-Network chooses the next action to perform as
$a_{t+1}=(0,0,0,1,0,0,0,0,0,0,0,0)$ (Heater = 0 W, Cooler = 200 W, Window =
OPEN).

Executing (red). Just as a usual MAPE-K, this component is in charge of
executing the action planned by the Planning component.

######  Example 6 (Running Case Study (6/6))

The Executor will closed the window and switch off the heater

###  2.1 Real Time Reinforcement Learning Agent

The classical idea of a RL-based agent is the following: the agent receives
the state from the environment, it selects an action according to a transition
distribution and it performs that action, changing the environment. This
scenario is modeled by using _Markov Decision Process_ (MDP) [25]. The
fundamental problem with this conception is that the environment is considered
paused during the action selection and that is why MDP is referred to _Turn-
based Markov Decision Process_ (TBMDP) in [23]. While it is suited for
problems such as board games, it is certainly not suited for real-time
applicationsIn order to derive the $\mathcal{M}_{\xi}^{\mathcal{A}}$, where
$\xi=1,2,3...,n$ is the index assign to each model and the superscript
$\mathcal{A}$ indicates that we are speaking about a model of the Agent, we
implemented a Dueling Deep Q-network (Figure 2), using the RTRL paradigm. In
RL, the goal is maximizing the reward function to let the agent learn how to
behave in a given environment. It has an experience replay and a target
network, performing a soft update, to optimize the training phase.

The exploration-exploitation trade-off is given by $\epsilon$: the chosen
action at time step $t$, will be a greedy action (exploit) with probability
(1-$\epsilon$) or may be a random action (explore) with probability of
$\epsilon$. Then it decrease the parameter $\epsilon$ by an $\epsilon$-decay
value. The soft update is regulated by the parameter $\tau$. The list of
parameters used in the study are shown in Table 3.

Parameter | Value
---|---
Memory Size | 1000000
Batch Size | 256
$\epsilon$-decay | 0.0002
$\tau$ | 0.005
Hidden Layers Neurons | 256
Training Episodes | 600
Table 3: Reinforcement Learning Agent parameters.

To enable the agent to calculate the reward of its actions and to manage the
training process, the Open AI Gym toolkit has been used [5]. The advantage of
using this tool is that it can be used by every type of algorithm and it does
not need any knowledge about the RL agent interacting with it. In fact, an
Open AI environment only needs the action chosen by the agent in order to
return the next state and the reward for that action. Every step of the
environment is considered to be scheduled at 15 minutes intervals.

####  2.1.1 Reward Function

A reward function balancing energy consumption and thermal comfort was chosen.
The two aspects are conflicting: in order to have an ideal thermal comfort,
consumption of energy is needed, while saving energy could result in thermal
discomfort. In this work, we modified the reward function defined in [4]. Let
us define the setpoint as the desired comfort temperature set by the user and
let us define two parameters $\epsilon,\delta>0$ such that $\epsilon<\delta$.

  * •

Comfort Zone: area in the interval $[setpoint-\epsilon,setpoint+\epsilon]$

  * •

Intermediate Zone: two areas enclosed in the intervals
$(setpoint-\delta,setpoint-\epsilon]$ and
$(setpoint+\epsilon,setpoint+\delta]$.

  * •

Danger Zone: the two remaining intervals $(-\infty,setpoint-\delta)$,
$(setpoint+\delta,+\infty)$.

The reward function is defined as follows:

$\begin{cases}0,&\text{if}\ setpoint-\epsilon\leq T_{in}\leq
setpoint+\epsilon\\\ -\beta\cdot energy-\rho\cdot(setpoint-
T_{in})^{2},&\text{if}\ setpoint-\delta\leq T_{in}<setpoint-\epsilon\\\ &\lor\
setpoint+\epsilon<T_{in}\leq setpoint+\delta\\\ -\beta\cdot
energy-\rho\cdot(|setpoint-T_{in}|)^{3},&\text{otherwise}\\\ \end{cases}$ (3)

where $T_{in}$ is the internal temperature, $\beta$ and $\rho$ are parameters
that permit to balance the trade-off between the energy consumption and the
thermal comfort. Table 4 shows the values of these parameters used in the
simulations.

Parameter | Value
---|---
Setpoint | 20
$\epsilon$ | 2
$\delta$ | 4
$\beta$ | 0.05
$\rho$ | 1
Table 4: Reward function parameters used.

##  3 Self-Adaptive strategy

To managing the two different kinds of novelties, we implement three different
self-adaptive strategies: the Model Switching Adaptation, the Architectural
Adaptation and the _DQN-Knowledge-based Adaptation_. The latter method always
starts only if the first two fail.

####  3.0.1 Model Switching Adaptation

REPTILE applies this strategy when a _context novelty_ is detected. All the
available pre-trained models $\mathcal{M}_{1}^{\mathcal{A}}$,
$\mathcal{M}_{2}^{\mathcal{A}}$,…,$\mathcal{M}_{k}^{\mathcal{A}}$ are
evaluated and the best one, in terms of some evaluation parameters, is
selected. Our selecting criterion is based on returning the greatest reward
using the current data.

####  3.0.2 Architectural Adaptation

This strategy of adaptation is based on the idea that the DRL-agent can learn
to perform new actions if it is aware of having instruments allowing them:
equipped with new devices there is a change in the number of potential actions
the agent can perform. The solution we adopted is pretty straightforward: we
modify the structure of our Dueling Deep-Q network, adding a number of input
neurons and output neuron that allow us to encode the new actions, too. Let us
consider a model $\mathcal{M}_{i}^{\mathcal{A}}$ that has been trained using
$N$ actions, if a new device is available we enrich the Dueling Deep-Q network
with a number $K$ of neurons such that the other $N+M$ actions can be taken
into account. Actually the re-modulation of the Dueling Deep-Q network
architecture does not imply automatically the re-training of the model
$\mathcal{M}_{i}^{\mathcal{A}}$: if the agent can perform well there in no
need to use the new device. That is the reason why the connections to and from
the new neurons are set to 0: in this way, the trained models can be used in a
context of $N+M$ actions, but maintaining the original behavior.

####  3.0.3 DQN-Knowledge-based Adaptation

This strategy answer the question: what if the other two adaptation methods
fail? If a novelty that cannot be managed using already available models, is
detected, REPTILE should use the data stored in the Knowledge for training a
new DRL-based agent both if there are new devices and if not.

###  3.1 Self-adaptation process

The complete self-adaptation process, shown in Figure 3, works in the
following way:

  1. 1.

The system continuously monitors the environment and updates the current
Knowledge.

  2. 2.

The Analyzing module checks for a possible novelty: if an architectural
novelty is detected, all current models are updated (neurons are added) and
the system goes back to the monitor phase (step 1); on the other hand, if it
is not an architectural novelty, the system selects the best model for the
current environment, which is the model allowing the RTRL agent to gain the
most cumulative reward; if no novelty is detected, the system goes back to the
monitor phase.

  3. 3.

Once the best model has been selected, the system waits $n$ steps in order to
observe the behavior of the agent using the selected model. The quantity $N$
can vary based on the application domain of the SAS.

  4. 4.

After $N$ steps, the system checks if the model fulfills all system’s goals
even in the presence of the novelty: if so, the system has successfully
adapted to the novelty and the monitor process starts again (step 1),
otherwise an alarm is sent to the maintenance: this is done because it means
that the system can not fulfill its goals with its current models,
communicating that it has to begin the training phase.

  5. 5.

Once the alarm has been fired, the system keeps selecting the best model at
its disposal in order to maintain a safe state while a new model is trained
using the new environment’s data stored in the Knowledge to generate episodes
for the RTRL agent.

  6. 6.

Once the training process is terminated, the system selects the new model and
observes the new behavior (step 3).

Figure 3: Model Switching Adaptation Diagram.

##  4 Results

In this paper we applied REPTILE to a HVAC system. Such a system was
implemented in Python and it is able to compute the internal temperature of a
room by simulating the actions of three devices: a heater, a cooler and a
window. The parameters used by the simulator are reported in Table 1. Since
REPTILE is based on DRL-basd agent we interfaced our simulator with the the
Open AI environment that invokes the temperature simulator to calculate the
Agents’s next state. REPTILE extends the MAPE-K framework with a _Learning_
module, based on DRL, and equips the _Knowledge_ module with a proactive
capability based on the possibility of creating on-the-fly time-varying models
of the environment and foreseen their next states.

Figure 4: Winter Agent. On the left, the monitoring of the temperature and the
actions taken by the agent for maintaining the temperature in the comfort
zone; on the right, the mazimization of the reward function is shown

The _Learning_ module comes into play whenever a novelty is detected and,
using the knowledge of the environment (or simulate the knowledge from the
collected data), it is able to trigger an adaptation strategy: the model
switching one, the system architectural one or the re-training one. The
_proactive_ Knowledge module is able to derive a model of the environment by
using the ARIMA formalism and, at the same time, contains some pre-trained
models of DRL-based agents, $\mathcal{M}_{1}^{\mathcal{A}}$,
$\mathcal{M}_{2}^{\mathcal{A}}$,…,$\mathcal{M}_{k}^{\mathcal{A}}$, and some
global models of the changing environment,
$\mathcal{M}_{1}^{\mathcal{G}}$,…,$\mathcal{M}_{j}^{\mathcal{G}}$, in order to
avoid to use some specific techniques to derive such models. In this work,
only the real data about external temperatures, $T_{ext}$, were used. They
consist of time series whose values were collected every hour [1, 2]. We pre-
processed those data by replicating each value 3 times in order to simulate a
sampling frequency of 15 minutes. We also have a model
$\mathcal{M}^{\mathcal{G}}$ and a model $\mathcal{M}^{\mathcal{A}}$. We also
have no available new devices. We show the results in two different operative
conditions: summer and winter. All the results described below are obtained
using a laptop with the Mint Linux OS, 8Gb RAM, a dual-core Intel<EMAIL_ADDRESS>CPU and a graphic card Nvidia GeForce GTX 940MX. Initially, we started by
having an agent that is trained to behave in winter conditions. We used 600
episodes, with a decreasing $\epsilon$-greedy exploration-exploitation trade-
off starting from 1 and an experience reply soft optimization. The learning
procedure took about 1 hour to be completed.

The behavior of such an agent is depicted in the left panel of Figure 4: the
upper block shows the internal temperature (blue line) remaining inside the
comfort zone, $18$°C - $22$°C (red horizontal lines), and the external
temperature (orange line).

(a) No adaptation

r107ptb25pt

(b) Adaptation

Figure 5: Test of Novelty Detection. The Agent can not go back in the Comfort
Zone, an alarm is sent and a new model is trained.

Moreover, the three blocks below display the actions performed by all the tree
devices. In particular only the heater keep on increasing and decreasing its
power while the cooler and the window do not act.

In order to test the _Model Switching Adaptation_ and the _DQN-Knowledge-based
Adaptation_ , we let the agent behaves in a summer scenario where external
temperatures range from 20 to 40 °C. This is a trick for being sure the agent
stumbles across a novelty. Figure 5(a) shows an example of triggering an
adaptation process. Horizontal red lines represent our comfort zone while the
blue line is the indoor temperature. After 15 time steps (red circle) the
_Analyzing_ modules detects, in a _reactive_ manner, a violation of our “soft”
constraint: not to be in the Intermediate zone for an hour. Now the system
starts the adaptation process finding no new devices (actually the search for
new devices is performed at each time step) and only one pre-trained model
available: $\mathcal{M}^{\mathcal{A}}$. Such a model has been prepared in a
way that the system can choose it (its reward is higher) but is not successful
in adapting. After 3 time steps (remember that the agent perform an action
every fifteen minutes), at time step 18, the indoor temperature is again out
of range; moreover, the third predicted point (the rightmost yellow point) is
in the Danger zone. At this point, given the impossibility of adapting with
the first two strategies, the system chooses to re-train a new DRL-based model
using the available $\mathcal{M}^{\mathcal{G}}$ notifying and alarm to the
system administrator. The new agent is trained with
$\mathcal{M}^{\mathcal{G}}$, containing summer environmental information, for
the same 600 episodes as the winter agent. The inset in Figure 5(b) shows the
maximization of the reward function. The training phase took 1 hour to be
completed in this case, too. Once the new agent has been trained, the system
switches form the current model to the new one and it is able to bring the
internal temperature back to the Comfort Zone, as depicted in Figure 5(b).

The above result underlines the ability of the system to adapt its behavior to
new operative conditions. In this case it is worth noticing that the knowledge
regarding the new environment (summer scenario) was already in the knowledge
module but REPTILE potentially allows to implement any technique to acquire
such a knowledge. Finally, in order to test the _Architectural Adaptation_ a
new heater with different power values (50 W, 200 W, 250W ,400W) was plugged
into the system.

Figure 6: Adaptation after a new device was plugged into the system. On the
left the dynamics of the monitoring of the system with the corresponding
actions taken by the DRL-based agent; on the right the maximization of the
reward function during the re-training process.

Once the system detected the new device, it immediately modified the
architecture of the Duelling Deep-Q network in order to be able to manage two
more actions and then it ran the adaptive procedure: in this the system was
running in summer, it was not able to manage the winter external condition so
it started the adaptive procedure. The procedure took again about one hour and
produce a behavior like the one shown in Figure 6. Comparing Figure 6 (left)
with Figure 4 it is evident the different actions the agent is performing in
order to keep the temperature in the comfort zone. It is also evident that the
agent continues to open and close the window and synchronously raises and
lowers the heating. The meaning of this behavior could be that the agent,
despite managing to keep the temperature in the comfort zone, would need more
episodes to optimize its behavior. The previous hypothesis can be confirmed if
we look at Figure 6 (right): even though the reward function is maximizing,
there are many large amplitude fluctuations. It is worth noticing that for
showing the effectiveness of _Architectural Adaptation_ the system was induced
to apply the _DQN-Knowledge-based Adaptation_ , too.

##  5 Conclusions

In this paper, we have presented REPTILE, a framework for the design of
proactive, real-time and self-adaptive systems based on _Deep Reinforcement
Learning_ (DRL). The proposed framework extends IBM’s MAPE-K loop with a _DRL
module_ that is used to see in advance potential dangerous and unwanted
situations, named _novelties_ , that may affect the expected behaviour of the
considered system. Two types of novelties are taken into account in our
framework: those related to the context/environment and those regarding the
physical architecture itself [12]. In the first case, a set of _models_ , that
are either built via a pre-existing knowledge or via the one acquired _on the
fly_ , are used to identify the best countermeasures/actions to execute in
order to guarantee the satisfaction of some given requirements. If no model is
available, or when an _architectural novelty_ is detected, the _DRL module_ is
updated to learn new behavioral policies and, consequently, to derive new
_models_ that will support the adaptation in the new configuration. When both
these strategies fail, a _DQN-Knowledge-based Adaptation_ is executed to re-
train the _DRL module_.

The proposed framework can be applied to both simulated and real/physical
systems. However, in the latter case, the appropriate software interfaces
should be developed to let REPTILE controls physical devices. To show an
application of REPTILE, in the paper, a _Heating, Ventilation and Air
Conditioning system_ (HVAC) has been used. The performed experiments confirm
the capability of our framework to guide the adaptation procedure. We can
observe that in the considered example, both the _Model Switching Adaptation_
and the _DQN-Knowledge-based Adaptation_ seem working well, due to the global
model $\mathcal{M}^{\mathcal{G}}$ ussed for training the new _DRL module_.

As far as we know, this is the first proposal taking into account both the
proactivity and the use of real-time deep reinforcement learning paradigm,
applied to a Deep Q-Network, for managing data streams.

Even if the _Architectural Adaptation_ showed good performances, we believe
that further studies are needed to evaluate the possibility of exploiting the
_transfer learning paradigm_ to avoid the completely re-training of the Deep-Q
Network [14].

Finally, one of the key point in REPTILE is $\mathcal{M}_{\xi}^{\mathcal{G}}$
that is the model used whenever we need to re-train the DRL-Agent. The better
is $\mathcal{M}_{\xi}^{\mathcal{G}}$, the more accurate is
$\mathcal{M}_{\xi}^{\mathcal{A}}$. These models can be derived directly from
data collected from existing systems by means of different data analysis
techniques [16, 24, 21, 6, 22]. In the case of the
$\mathcal{M}_{\xi}^{\mathcal{G}}$ models considered in this paper, existing
time series have been used. Nevertheless, further studies are needed to
identify the algorithms that are able to derive a model from a (possibly
small) set of data, or fast enough to deal with a vary large amount of data.

## References

  * [1] Beach weather stations - automated sensors (2020), https://catalog.data.gov/dataset/beach-weather-stations-automated-sensors
  * [2] Beniaguev, D.: Historical hourly weather data 2012-2017 (2017), https://www.kaggle.com/selfishgene/historical-hourly-weather-data
  * [3] Box, G.E., Jenkins, G.M., Reinsel, G.C., Ljung, G.M.: Time series analysis: forecasting and control. John Wiley & Sons (2015)
  * [4] Brandi, S., Piscitelli, M.S., Martellacci, M., Capozzoli, A.: Deep reinforcement learning to optimise indoor temperature control and heating energy consumption in buildings. Energy and Buildings 224, 110225 (2020)
  * [5] Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J., Zaremba, W.: Openai gym. arXiv preprint arXiv:1606.01540 (2016)
  * [6] De Simone, A., Piangerelli, M.: A bayesian approach for monitoring epidemics in presence of undetected cases. Chaos, Solitons & Fractals 140, 110167 (2020)
  * [7] Elkhodary, A., Esfahani, N., Malek, S.: Fusion: a framework for engineering self-tuning self-adaptive software systems. In: Proceedings of the eighteenth ACM SIGSOFT international symposium on Foundations of software engineering. pp. 7–16 (2010)
  * [8] Garlan, D., Cheng, S.W., Huang, A.C., Schmerl, B., Steenkiste, P.: Rainbow: Architecture-based self-adaptation with reusable infrastructure. Computer 37(10), 46–54 (2004)
  * [9] Hallsteinsen, S., Geihs, K., Paspallis, N., Eliassen, F., Horn, G., Lorenzo, J., Mamelli, A., Papadopoulos, G.A.: A development framework and methodology for self-adapting applications in ubiquitous computing environments. Journal of Systems and Software 85(12), 2840–2859 (2012)
  * [10] Kephart, J.O., Chess, D.M.: The vision of autonomic computing. Computer 36(1), 41–50 (2003)
  * [11] Kim, D., Park, S.: Reinforcement learning-based dynamic adaptation planning method for architecture-based self-managed software. In: 2009 ICSE Workshop on Software Engineering for Adaptive and Self-Managing Systems. pp. 76–85. IEEE (2009)
  * [12] Klös, V., Göthel, T., Glesner, S.: Adaptive knowledge bases in self-adaptive system design. In: 2015 41st Euromicro Conference on Software Engineering and Advanced Applications. pp. 472–478. IEEE (2015)
  * [13] Krupitzer, C., Roth, F.M., VanSyckel, S., Schiele, G., Becker, C.: A survey on engineering approaches for self-adaptive systems. Pervasive and Mobile Computing 17, 184–206 (2015)
  * [14] Kwiatkowski, R., Lipson, H.: Task-agnostic self-modeling machines. Science Robotics 4(26), eaau9354 (2019)
  * [15] Maimó, L.F., Gómez, Á.L.P., Clemente, F.J.G., Pérez, M.G., Pérez, G.M.: A self-adaptive deep learning-based system for anomaly detection in 5g networks. IEEE Access 6, 7700–7712 (2018)
  * [16] Merelli, E., Piangerelli, M.: Rnn-based model for self-adaptive systems (2014)
  * [17] Muccini, H., Vaidhyanathan, K.: A machine learning-driven approach for proactive decision making in adaptive architectures. In: 2019 IEEE International Conference on Software Architecture Companion (ICSA-C). pp. 242–245. IEEE (2019)
  * [18] Müller, H.A., Pezzè, M., Shaw, M.: Visibility of control in adaptive systems. In: ULSSIS@ ICSE. pp. 23–26 (2008)
  * [19] Papamartzivanos, D., Mármol, F.G., Kambourakis, G.: Introducing deep learning self-adaptive misuse network intrusion detection systems. IEEE Access 7, 13546–13560 (2019)
  * [20] Pascual, G.G., Pinto, M., Fuentes, L.: Run-time adaptation of mobile applications using genetic algorithms. In: 2013 8th International Symposium on Software Engineering for Adaptive and Self-Managing Systems (SEAMS). pp. 73–82. IEEE (2013)
  * [21] Piangerelli, M., Maestri, S., Merelli, E.: Visualising 2-simplex formation in metabolic reactions. Journal of Molecular Graphics and Modelling 97, 107576 (2020)
  * [22] Quadrini, M., Cavallin, M., Daberdaku, S., Ferrari, C.: Prosps: Protein sites prediction based on sequence fragments. In: International Conference on Machine Learning, Optimization, and Data Science. pp. 568–580. Springer (2021)
  * [23] Ramstedt, S., Pal, C.: Real-time reinforcement learning. In: Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 32. Curran Associates, Inc. (2019), https://proceedings.neurips.cc/paper/2019/file/54e36c5ff5f6a1802925ca009f3ebb68-Paper.pdf
  * [24] Rucco, M.: Survey of topdrim applications of topological data analysis
  * [25] Sutton, R.S., Barto, A.G.: Reinforcement learning: An introduction. MIT press (2018)
  * [26] Wang, Z., Schaul, T., Hessel, M., Hasselt, H., Lanctot, M., Freitas, N.: Dueling network architectures for deep reinforcement learning. In: International conference on machine learning. pp. 1995–2003. PMLR (2016)
  * [27] Weyns, D., Malek, S., Andersson, J.: Forms: a formal reference model for self-adaptation. In: Proceedings of the 7th international conference on Autonomic computing. pp. 205–214 (2010)


Also at ]Electronics Engineering Department, Universitat de València, Avinguda
de l’Universitat, 46100, Burjassot, Spain.

# Novel reaction force for ultra-relativistic dynamics of a classical point
charge

P. Martín-Luna  D. González-Iglesias  B. Gimeno  D. Esperante  [  C. Blanch
N. Fuster-Martínez  P. Martinez-Reviriego  J. Fuster  Instituto de Física
Corpuscular (IFIC), Universitat de València - Consejo Superior de
Investigaciones Científicas, Calle Catedrático José Beltrán, 2, 46980,
Paterna, València, Spain

###### Abstract

The problem of the electromagnetic radiation of an accelerated charged
particle is one of the most controversial issues in Physics since the
beginning of the last century, representing one of the most popular unsolved
problems of the Modern Physics. Different equations of motion have been
proposed throughout history for a point charge including the electromagnetic
radiation emitted, but all these expressions show some limitations. An
equation based on the principle of conservation of energy is proposed in this
work for the ultra-relativistic motion. Different examples are analyzed
showing that the energy lost by the charge agrees with the Liénard formula.
This proposed equation has been compared with the Landau-Lifshitz equation
obtaining a good agreement in the range of application of the Landau-Lifshitz
formula.

##  I Introduction

When a charged particle is accelerated it radiates electromagnetic energy
according to the well-known Liénard formula [1]. The energy lost comes from
the particle’s kinetic energy and consequently a reaction force must act on
the radiating particle for accounting such loss of kinetic energy. The
equation of motion of a charged particle including the reaction force has been
a subject of interest during the last 100 years and it has become one of the
most popular unsolved problems of Modern Physics. The first formula of the
reaction force was the Abraham-Lorentz formula [2, 3], which was extended to a
Lorentz-covariant form by Dirac in 1938 [4]. This expression is called the
Lorentz-Abraham-Dirac (LAD) equation and suffers from some problems. Firstly,
the LAD equation involves the temporal derivative of the acceleration and
hence the initial position and velocity of the particle does not determine the
solution uniquely. For this reason, one extra initial condition is needed
compared to the Newtonian mechanics. Secondly, if the external force is zero
the LAD equation admits solutions that accelerate exponentially in addition to
the trivial solution of uniform motion. These solutions are called runaway
solutions and violate the principle of conservation of energy. The runaway
solutions can be avoided by artificially introducing a preacceleration, but
then the causality is violated. Furthermore, the LAD equation does not agree
with the theory of Compton scattering except in the weak field case [5].

As a consequence of the inherent difficulties of the LAD equation, several
alternative formulations have been proposed. Mo and Papas proposed in 1971 a
new radiation reaction force formula that avoids the runaway solutions as well
as the preacceleration phenomenon, assuming a reaction force proportional to
the acceleration of the particle rather than the velocity [6]. Ford and
O’Connell proposed in 1991 an alternative equation based on the use of the
generalized quantum Langevin equation for an electron with a finite size [7].
Hartemann and Luhman [8], Yaghjian [9], and Hammond [10] have proposed other
equations, but the most famous reaction force is the Landau-Lifshitz (LL)
equation [11] that can be obtained assuming that the radiation reaction force
is much smaller than the external electromagnetic force. But the LL equation
can also be obtained as a limit of the self-force modeling the electron as a
sphere [12, 13]. Thus, the LL theory has been proposed as the correct equation
of motion of a classical point charge [14] and it has even been recently
verified experimentally measuring the emission spectra of electrons and
positrons penetrating into aligned single crystals [15]. However, the LL
equation is not totally correct since it predicts no reaction force when the
particle has a linear acceleration although in this case the radiation losses
are practically negligible for typical electric fields [1]. Moreover, the LL
equation fails for electromagnetic fields with abrupt changes, i.e. for high
frequencies and/or high intensities [16].

The purpose of this work is the derivation of a new equation that overcomes
these difficulties. The article is organized as follows. In Sec. II the basic
principles and the approximations to obtain the new proposed reaction force
are discussed. In Sec. III the proposed reaction force is solved for some
simple physical examples comparing the results with the LL equation. We
conclude the paper with final remarks in Sec. IV.

##  II The proposed reaction force

Radiation from charged particles becomes important for ultra-relativistic
particles and its emission is concentrated in the direction of motion of the
particle with a symmetric angular distribution around this direction [1].
Consequently, we will assume that the reaction force is antiparallel to the
velocity $\bf{v}$, i.e.

$\mathbf{F}_{\mathrm{R}}=-F_{\mathrm{R}}\hat{\mathbf{v},}$ (1)

with $\hat{\mathbf{v}}=\mathbf{v}/\|\mathbf{v}\|$, in order to satisfy in
average the principle of conservation of energy and momentum when the emission
is produced. Note that $\|\|$ indicates the euclidean norm of a vector on the
3-dimensional space $\mathbb{R}^{3}$. On the other hand, in the average
distance, called mean free path [17]

$\lambda=\frac{2\sqrt{3}\rho}{5\alpha\gamma},$ (2)

is emitted a photon of mean energy

$\langle E\rangle=\frac{8}{15\sqrt{3}}E_{c},$ (3)

where $E_{\mathrm{c}}=\frac{3}{2}\hbar c\gamma^{3}/\rho$ is the critical
energy, $\rho$ is the instantaneous radius of curvature of the trajectory,
$\gamma$ is the Lorentz relativistic factor, $\hbar$ is the reduced Planck
constant, $c$ is the speed of light in vacuum and $\alpha\approx\frac{1}{137}$
is the fine-structure constant. Thus, the principle of conservation of energy
during a small displacement $\Delta x$ gives the relation
$F_{\mathrm{R}}\Delta x=\frac{\Delta x}{\lambda}\langle E\rangle$, i.e.

$F_{\mathrm{R}}=\frac{\langle E\rangle}{\lambda},$ (4)

which depends only on $\gamma$ and $\rho$. The radiation emitted by an
extremely relativistic charge is approximately the same as that of a particle
moving instantaneously in a circular arc of radius of curvature [1]:

$\rho=\frac{\|\mathbf{v}\|^{2}}{\|\dot{\mathbf{v}}_{\perp}\|}\approx\frac{c^{2}}{\|\dot{\mathbf{v}}_{\perp}\|},$
(5)

where $\|\dot{\mathbf{v}}_{\perp}\|$ is the magnitude of the perpendicular
component with respect to $\mathbf{v}$ of the acceleration. Note that the dot
indicates a time derivative. The equation of motion of the particle will be

$\dot{\mathbf{p}}=m\dot{\gamma}\mathbf{v}+m\gamma\dot{\mathbf{v}}=\mathbf{F}_{\mathrm{ext}}+\mathbf{F}_{\mathrm{R}},$
(6)

where $\mathbf{p}$ is the relativistic linear momentum, $m$ is the rest mass
of the particle and $\mathbf{F}_{\mathrm{ext}}$ is the external force applied
to the charged particle, e.g. the Lorentz force. Therefore, projecting the
equation of motion in the perpendicular direction, we obtain

$\dot{\mathbf{v}}_{\perp}=\frac{\mathbf{F}_{\mathrm{ext},\perp}}{m\gamma}.$
(7)

Hence, substituting in Eq. (4), the following reaction force is obtained

$\mathbf{F}_{\mathrm{R}}=-\frac{\tau_{m}}{mc}\beta^{-4}\gamma^{2}\|\mathbf{F}_{\text{ext
},\perp}\|^{2}\hat{\mathbf{v}}\approx-\frac{\tau_{m}}{mc}\gamma^{2}\|\mathbf{F}_{\text{ext
},\perp}\|^{2}\hat{\mathbf{v}},$ (8)

where the characteristic time
$\tau_{m}=\frac{2}{3}\frac{1}{4\pi\varepsilon_{0}}\frac{e^{2}}{mc^{3}}$ has
been defined; $e$ is the elementary charge, $\varepsilon_{0}$ is the vacuum
electric permittivity and $\beta=v/c$. This proposed reaction force will be a
great approximation for ultra-relativistic particles. Note that if the
acceleration is parallel to the velocity, the Eq. (8) erroneously predicts a
null reaction force as e.g. the Landau-Lifshitz force. Nevertheless, for a
given magnitude of applied force, the radiation emitted due to a perpendicular
acceleration is a factor $\gamma^{2}$ larger than with a parallel acceleration
[1]. Consequently, (8) suggests the addition of an identical term for the
parallel component of the force without the factor $\gamma^{2}$. Thus, the
following final expression is obtained

$\mathbf{F}_{\mathrm{R}}=-\frac{\tau_{m}}{mc}\left(\gamma^{2}\|\mathbf{F}_{\mathrm{ext},\perp}\|^{2}+\|\mathbf{F}_{\mathrm{ext},\|}\|^{2}\right)\hat{\mathbf{v}}.$
(9)

As a first check, we see that in absence of the external force the reaction
force is zero. An alternative demonstration using four-vectors is included in
the Appendix.

##  III SPECIAL CASES AND COMPARISON WITH THE LANDAU-LIFSHITZ EQUATION

Now we will solve the equation of motion (6) assuming the proposed reaction
force (9) for an electron (with charge $q=-e$) in different basic physical
situations comparing the results with those of the Landau-Lifshitz reaction
force equation [11]

$\displaystyle F_{\mathrm{LL}}^{\mu}=\tau_{m}\left(\frac{q}{c}\frac{\partial
F^{\mu\nu}}{\partial
x^{\gamma}}u_{\nu}u^{\gamma}+\frac{q^{2}}{mc}F^{\mu\gamma}F_{\nu\gamma}u^{\nu}+\right.$
(10)
$\displaystyle\left.\frac{q^{2}}{mc^{3}}\left(F_{\nu\gamma}u^{\gamma}\right)\left(F^{\nu\alpha}u_{\alpha}\right)u^{\mu}\right),$

where $F^{\mu\nu}$ is the electromagnetic field tensor [1] and $u^{\mu}$ is
the four-velocity; the metric tensor is
$g^{\mu\nu}=\mathrm{diag}(1,-1,-1,-1)$. In the three-dimensional form it can
be written for an electron as [18]

$\displaystyle\mathbf{F}_{\mathbf{LL}}$
$\displaystyle=\tau_{m}\left\\{e\gamma\left(\left[\frac{\partial}{\partial
t}+(\mathbf{v}\cdot\nabla)\right]\mathbf{E}+\left\\{\mathbf{v}\times\left[\frac{\partial}{\partial
t}+(\mathbf{v}\cdot\nabla)\right]\mathbf{B}\right\\}\right)\right.$ (11)
$\displaystyle+\frac{e^{2}}{mc}\left(\mathbf{E}\times
c\mathbf{B}+c[\mathbf{B}\times(\mathbf{B}\times\mathbf{v})]+\frac{1}{c}\mathbf{E}(\mathbf{v}\cdot\mathbf{E})\right)$
$\displaystyle\left.-\frac{e^{2}\gamma^{2}}{mc^{2}}\mathbf{v}\left[(\mathbf{E}+\mathbf{v}\times\mathbf{B})^{2}-\frac{1}{c^{2}}(\mathbf{v}\cdot\mathbf{E})^{2}\right]\right\\}$

where $\mathbf{E}$ and $\mathbf{B}$ are the external electric and magnetic
field, respectively.

###  III.1 Motion perpendicular to uniform magnetic field

In this first example, we are going to assume a uniform magnetostatic field
$\mathbf{B}=B_{0}\hat{\mathbf{z}}$ and motion in the $xy$-plane. Then, the
proposed reaction force (9) and the LL force (11) are simply

$\mathbf{F}_{\mathbf{R}}=\mathbf{F}_{\mathbf{LL}}=-\frac{\tau_{m}e^{2}}{m}\gamma^{2}B_{0}^{2}\mathbf{v},$
(12)

where $v\approx c$ has been assumed. The solution of the motion for the
general initial conditions $(x_{0},y_{0}),(v_{x0},v_{y0})$ is

$\begin{gathered}v_{x}(t)=e^{-\sigma\tau}\left[v_{x0}\cos(\omega_{c}\tau)+v_{y0}\sin(\omega_{c}\tau)\right],\\\
v_{y}(t)=e^{-\sigma\tau}\left[-v_{x0}\sin(\omega_{c}\tau)+v_{y0}\cos(\omega_{c}\tau)\right],\\\
x(t)=x_{0}+v_{x0}I_{1}(\tau)+v_{y0}I_{2}(\tau),\\\
y(t)=y_{0}-v_{x0}I_{2}(\tau)+v_{y0}I_{1}(\tau),\end{gathered}$ (13)

where we have defined the auxiliary integrals

$\begin{gathered}I_{1}(\tau)=\int_{0}^{\tau}\frac{e^{-\sigma\tau}}{\sqrt{1-\beta_{0}^{2}e^{-2\sigma\tau}}}\cos(\omega_{c}\tau)d\tau,\\\
I_{2}(\tau)=\int_{0}^{\tau}\frac{e^{-\sigma\tau}}{\sqrt{1-\beta_{0}^{2}e^{-2\sigma\tau}}}\sin(\omega_{c}\tau)d\tau,\end{gathered}$
(14)

the following quantities $\omega_{c}=eB_{0}/m,\sigma=\tau_{m}\omega_{c}^{2}$
and the proper time $\tau$ that can be expressed in the laboratory frame time
$t$ as

$\tau=\frac{1}{\sigma}\ln\left(\frac{\delta-e^{2\sigma t}}{\delta-1}\right)-t$
(15)

with $\delta=\frac{1-\gamma_{0}}{1+\gamma_{0}}$. Note that
$\beta_{0}=v_{0}/c$, $v_{0}=(v_{x0}^{2}+v_{y0}^{2})^{1/2}$ and
$\gamma_{0}=(1-\beta_{0}^{2})^{-{1/2}}$.

Furthermore, we can calculate the energy lost per unit time
$-\frac{dE}{dt}=-mc^{2}\frac{d\gamma}{dt}$, obtaining

$-\frac{dE}{dt}=\frac{\tau_{m}\gamma^{2}}{mc}e^{2}B_{0}^{2}v^{2}$ (16)

that agrees with the Liénard formula [1]

$P=\frac{\tau_{m}\gamma^{2}}{mc}\left(\dot{\mathbf{p}}^{2}-\beta^{2}\dot{p}^{2}\right),$
(17)

if it is particularized for the considered uniform magnetic field.

Figure 1 shows the trajectory (13) in the $xy$-plane of an ultra-relativistic
electron. As it is expected, the motion is a spiral due to the radiated energy
by the accelerated electron.

Figure 1: Motion of an electron with an initial total energy of 100 GeV in a
perpendicular magnetic field of $B=1$ T for the first 50 $\mu$s (the vectors
indicate the initial and final velocity). The dashed line shows the circular
trajectory with radius $x_{0}=m\gamma_{0}v_{0}/(eB)$ obtained if the radiation
is neglected.

###  III.2 Motion parallel to uniform electric field

In this case we will assume a homogeneous static electric field
$\mathbf{E}=E_{0}\hat{\mathbf{z}}$ and motion in the $z$-axis. Hence, the
proposed reaction force (9) is

$\mathbf{F}_{\mathbf{R}}=-\frac{\tau_{m}e^{2}}{mc}E_{0}^{2}\hat{\mathbf{v}}.$
(18)

and the Landau-Lifshitz force (11) is $\mathbf{F}_{\mathbf{LL}}=\mathbf{0}$.
The general solution of this motion is

$u_{z}(t)=\gamma(t)v_{z}(t)=u_{z0}+\frac{qE_{0}}{m}(1-\mathrm{sgn}(v_{z})\epsilon)t,$
(19)

$z(t)=z_{0}+\frac{mc^{2}}{qE_{0}(1-\mathrm{sgn}(v_{z})\epsilon)}(\gamma(t)-\gamma_{0}),$
(20)

where the dimensionless parameter $\epsilon=\frac{\tau_{m}eE_{0}}{mc}\approx
3.68\times 10^{-21}E_{0}(\mathrm{V/m})$ has been defined and $\mathrm{sgn}(x)$
is the sign function. The solution of the Landau-Lifshitz force, i.e.
neglecting the losses due to the radiation, is obtained by substituting
$\epsilon=0$. Thus, it can be seen that for typical values of the electric
fields in conventional particle accelerators, i.e, $E_{0}\sim 100$ MV/m, this
parameter includes a small perturbation to the motion obtained neglecting the
effects of radiation. The cumulative energy lost as a function of time is
given by

$\Delta
E(t)=mc^{2}(\gamma(t,\epsilon=0)-\gamma(t,\epsilon))=qE_{0}v_{z}t\epsilon+O\left(\epsilon^{2}\right).$
(21)

If we make the approximation $v_{z}\approx c$ this equation is analogous to
the Liénard equation (17), as it is expected. Hence, our proposed reaction
force will be valid as long as $\epsilon\ll 1$, i.e. up to electric fields of
the order of $10^{19}$ V/m. It is important to remark that the effects of the
radiation are in general totally negligible in this case. For example, an
electron initially at rest accelerated with an electric field of 1 GV/m during
100 m has finally an energy of $\sim$100 GeV, but the predicted losses due to
the radiation are $\sim$0.36 eV in such distance. However, Eq. (21) shows that
the energy lost as a function of time is proportional to the square of the
electric field. Consequently, if the electric field is increased four or five
orders of magnitude (i.e. for electric fields $\sim$10-100 TV/m) the radiated
energy will become very important. These ultrahigh electric fields have been
predicted using X-ray wakefield acceleration in metallic crystals [19, 20] and
hence the LL equation cannot be used in their simulations, whereas the novel
reaction force presented in this manuscript will be a good approximation.

###  III.3 Motion in an electromagnetic plane wave

In this section we are going to study the motion of an electron due to a
linear polarized electromagnetic plane wave with fields
$\mathbf{E}=E_{0}\cos(\omega t-kz)\hat{\mathbf{x}}$,
$\mathbf{B}=B_{0}\cos(\omega t-kz)\hat{\mathbf{y}}$ with $\omega=ck$ and
$E_{0}=cB_{0}$. The fourth order Runge-Kutta method has been used to
numerically solve the electron relativistic dynamics. The radiated power can
be calculated as

$P_{\mathrm{rad}}=\mathbf{F}_{\mathbf{ext}}\cdot{\mathbf{v}}-\frac{dE}{dt},$
(22)

i.e., the difference between the power obtained from the external force (the
Lorentz force) and the variation of the total energy per unit of time.

Figure 2: Comparison between the trajectory obtained using the proposed
reaction force (R) and the LL equation for $E_{0}=1$ TV/m and (a)
$\omega=10^{16}$ rad/s, (b) $\omega=10^{17}$ rad/s, (c) $\omega=10^{18}$
rad/s. The initial conditions are $\mathbf{r_{0}}=\mathbf{0}$,
$\mathbf{v_{0}}=v_{0}\hat{\mathbf{x}}$ with an initial electron total energy
of 10 GeV.  Figure 3: Comparison between the Liénard formula and the radiated
power associated with the motion obtained using the proposed reaction force
(a) and the LL equation (b) for $E_{0}=1$ TV/m, $\omega=10^{18}$ rad/s, and
the initial conditions $\mathbf{r_{0}}=\mathbf{0}$,
$\mathbf{v_{0}}=v_{0}\hat{\mathbf{x}}$ and an initial electron total energy of
10 GeV. The time is normalized to the period of the plane wave
$T=2\pi/\omega$.

Figure 2 shows the comparison between the trajectories obtained using the
proposed reaction force (9) and the LL equation (11) for different frequencies
$\omega$, showing that the disagreement between the trajectories becomes very
important at higher frequencies. Furthermore, Fig. 3 shows the comparison
between the Liénard formula (17) and the radiated power (22) when the electron
relativistic dynamics is numerically solved using the proposed reaction force
and the LL equation. It can be seen that the new reaction force agrees with
the Liénard formula, while the LL equation does not have a good accordance. It
is interesting to remark that the radiated power associated with the proposed
reaction and the LL equation is similar, but the power predicted by the
Liénard formula is different in each case since the motion is different (see
Figure 2).

Figure 4: Discrepancies between the Liénard formula and the radiated power
associated with the motion obtained using the proposed reaction force (R) and
the LL equation for different frequencies $\omega$: (a) $\omega=10^{15}$
rad/s, (b) $\omega=10^{16}$ rad/s, (c) $\omega=10^{17}$ rad/s and (d)
$\omega=10^{18}$ rad/s. Same $E_{0}$ and initial conditions as in Fig. 3. The
discrepancies are normalized to the maximum value of the corresponding Liénard
formula and expressed in %, and the time is normalized to the period of the
plane wave $T=2\pi/\omega$.

Figure 4 shows the difference between the Liénard formula and the radiated
power for different frequencies $\omega$ in order to study the differences
between the two forces when the electromagnetic fields have very rapid
changes. Firstly, we can see that the discrepancy between the radiated power
obtained using the LL equation and the Liénard equation increases linearly
with the frequency. Nevertheless, the discrepancy is negligible for the
proposed reaction force. Therefore, it has been verified that the emitted
radiation calculated with the proposed reaction force agrees much better to
the Liénard formula for high frequencies than the LL equation which
definitively fails.

###  III.4 Motion in an undulator

Figure 5: Scheme of electron motion in an undulator. Adapted from [21].

Finally, we are going to consider an ultra-relativistic electron moving in a
planar undulator as it is depicted in Figure 5. In the planar undulator, near
the $z$-axis, the magnetic field can be approximated by [21]

$B_{y}(z)=B_{0}\sin\left(\frac{2\pi}{\lambda_{u}}z\right),$ (23)

where $\lambda_{u}$ is the undulator period length. If we consider a
trajectory in the $xz$-plane and assume that the energy is constant (i.e.
$\gamma$ is constant), the horizontal component is given by

$x(z)=-\frac{K}{k_{u}\gamma\beta_{z}}\sin\left(k_{u}z\right),$ (24)

where we have imposed an initial velocity
$v_{x}(z=0)=-Kc/\gamma,\|\mathbf{v}\|=\beta c$, $\mathbf{x_{0}}=\mathbf{0}$;
$k_{u}=2\pi/\lambda_{u}$, $K=\frac{eB_{0}}{m_{e}ck_{u}}$ is the undulator
parameter and

$\beta_{z}=\beta\left(1-\frac{K^{2}}{4\gamma^{2}}-\frac{K^{2}}{4\gamma^{2}}\cos\left(2k_{u}z\right)\right)$
(25)

is the velocity (in units of $c$) of the electron in the $z$-direction.

In this example, the proposed reaction force is given by

$\mathbf{F}_{\mathbf{R}}=-\frac{\tau_{m}e^{2}}{m}\gamma^{2}B_{0}^{2}\sin^{2}\left(k_{u}z\right)\mathbf{v},$
(26)

and the LL equation is this Eq. (26) plus a term due to the spatial derivative
of the magnetic field. Nevertheless, for typical values of the undulator
period (about tens of mm) this additional term is negligible since the
magnetic field has not a sufficient strong spatial variation. Consequently,
the discrepancies between the proposed reaction force and the LL equation when
they are numerically solved are totally negligible.

Figure 6 displays the trajectory obtained solving the proposed reaction force
compared to the Eq. (24) showing that the difference is negligible for typical
values even when the electron energy is 500 GeV. However, this approximation
assumes no radiation and therefore cannot be used in order to take into
account the effects of the electron electromagnetic energy loss in the
radiation spectrum.

Figure 6: (a) Comparison between the trajectory obtained using the proposed
reaction force (R) and the analytical approximation (24) of an electron with
an initial total energy of 17.5 GeV and in a planar undulator with $K=4$ and
$\lambda_{u}=40$ mm after 1000 periods. These numerical values correspond to
the European XFEL parameters [22]. (b) Same as (a) but with an initial energy
of 500 GeV.

The energy lost per one undulator period can be calculated using the Liénard
expression obtaining [22]

$\Delta E(\text{period})=\frac{2}{3}\pi\alpha
hc\frac{\gamma^{2}K^{2}}{\lambda_{u}},$ (27)

and the relative energy loss $\delta$ is given by

$\delta=\frac{\Delta E(\text{period})}{m\gamma
c^{2}}=\frac{(2\pi)^{2}}{3}\gamma K^{2}\frac{r_{e}}{\lambda_{u}},$ (28)

where $r_{e}=\frac{\alpha\hbar}{m_{e}c}$ is the classical electron radius.

Figure 7: Comparison between the energy lost per one undulator period using
the proposed reaction force (R) and the Liénard expression (27) as a function
of the number of period $n$. The electron has an initial total energy of 17.5
GeV and the planar undulator parameters are $K=4$ and $\lambda_{u}=40$ mm.

Figure 7 shows the energy lost per one period as a function of the number of
period of the undulator. It can be seen that the result obtained with the
proposed reaction force agrees with the Liénard expression. As a consequence
of the electron energy decrease along the trajectory, the electromagnetic
radiation from different undulators drops out of synchronism. Thus, the
proposed reaction force can be easily used to calculate the radiation spectrum
taking into account this effect. The interfering factor associated with the
emission of electromagnetic radiation in each period of the undulator is given
by [22]

$I_{\mathrm{int}}(\lambda)=\sum_{n=1}^{N}\exp\left(i\Phi\left(z_{n}\right)\right),$
(29)

where

$\Phi\left(z_{n+1}\right)-\Phi\left(z_{n}\right)=2\pi\frac{\lambda_{1,n}}{\lambda}$
(30)

is the difference of phase between the emission in the $(n+1)$-th and the
$n$-th period of the undulator, $N$ is the total number of periods, $\lambda$
is the wavelength of the emitted radiation and

$\lambda_{1,n}=\left(\frac{\widetilde{AC}}{\beta}-\overline{AC}\right)_{n},$
(31)

$\widetilde{AC}$ being the length of the electron trajectory and
$\overline{AC}=\lambda_{u}$ the length of the photon trajectory (see Figure
8).

Figure 8: Scheme of the electron trajectory in order to calculate the
resonance condition of the emitted radiation.

When we assume that $\gamma$ is constant, the interfering factor is given by

$I_{\text{int
}}(\lambda)=\frac{\sin\left(N\pi\frac{\lambda_{1}}{\lambda}\right)}{\sin\left(\pi\frac{\lambda_{1}}{\lambda}\right)}$
(32)

with
$\lambda_{1}=\frac{\lambda_{u}}{2\gamma_{0}^{2}}\left(1+\frac{K^{2}}{2}\right)$.
We can deduce from this expression that the resonance occurs at the
wavelengths $\lambda=\lambda_{1}/k$ being $k$ an integer (the $k$-th harmonic)
and the width of the maxima is inversely proportional to $N$. Otherwise, if we
assume that the relative energy $\delta$ is constant, the phase is given by
[22]

$\Phi\left(z_{n}\right)=\frac{\pi\lambda_{u}}{\lambda\gamma_{0}^{2}}\left(1+\frac{K^{2}}{2}\right)(n-1)(1+\delta(n-2)),$
(33)

where $\gamma_{0}$ is the initial relativistic factor. On the other hand,
using the proposed reaction force to numerically solve the electron dynamics,
we can calculate the value of $\lambda_{1,n}$ as a discretization of the
integral of the path described by the electron:

$\lambda_{1,n}=\left(\frac{\widetilde{AC}}{\beta}-\overline{AC}\right)_{n}=\left(\sum_{j}^{A\rightarrow
C}\frac{|\Delta\vec{r}|_{j}}{\beta_{j}}-\left|z_{C}-z_{A}\right|\right)_{n}.$
(34)

Thus, we are going to compare the interfering factor using the three different
approximations: Eq. (32)-(34). Figure 9 shows the normalized intensity as a
function of the energy of the emitted photons for the first and second
harmonic. It can be seen that the approximation of $\delta$ constant (33) and
using the proposed reaction force (34) agree very well. As a consequence of
the energy lost by the electron along the trajectory, the peak maximum
decreases and is redshifted. Furthermore, the spectrum is also wider compared
to the case neglecting the energy loss. These three effects can be clearly
seen in Figure 10 where the number of periods $N$ and the initial energy have
been increased. In this case, the first harmonic has a very wide spectrum and
the normalized intensity is lower than the 2 $\%$ of the peak maximum
predicted neglecting the energy loss. Moreover, it can be seen that the
approximation of $\delta$ constant agrees well with the case using the
proposed reaction force, although in this extreme case there is a small
discrepancy.

Thus, in this last example we have checked that the proposed reaction force is
consistent with the typical approximations that are made in the undulators
theory. Furthermore, the way that we have proposed in order to calculate the
spectrum using the reaction force can be evidently used for any magnetic
field. For example, the effects of the magnetic field nonuniformity along the
horizontal X-axis due to finite poles could be easily studied.

Figure 9: Spectra of the electromagnetic radiation as a function of the photon
energy $hc/\lambda$. Normalized intensity $|I_{\text{int }}/N|^{2}$ associated
with an electron with an initial total energy of 17.5 GeV and a planar
undulator with parameters: $K=4$ and $\lambda_{u}=40$ mm, $N=1000$. (a) First
harmonic and (b) second harmonic. Figure 10: Spectra of the electromagnetic
radiation as a function of the photon energy $hc/\lambda$. (a) Normalized
intensity $|I_{\text{int }}/N|^{2}$ associated with an electron with an
initial total energy of 175 GeV and a planar undulator with parameters: $K=4$
and $\lambda_{u}=40$ mm, $N=3000$. (b) Zoom of the spectra calculated using
Eq. (33-34).

##  IV Conclusions

A new reaction force has been proposed to take into account the
electromagnetic radiation emitted from an accelerated charged particle. The
proposed reaction force works equal or better than the Landau-Lifshitz force
in the analyzed cases. Firstly, this new reaction force allows to introduce
the radiation when the motion is parallel to an electric field, while the
Landau-Lifshitz equation predicts no radiation in such case. Secondly, the
proposed reaction force agrees much better than the Landau-Lifshitz force with
the Liénard formula for electromagnetic fields with rapid changes. Moreover,
this new reaction force is much easier to numerically solve and implement
since it does not depend on the derivatives of the electromagnetic fields.

###### Acknowledgements.

This work has been supported by Ministerio de Universidades (Gobierno de
España) under grant number FPU20/04958.

## Appendix: Alternative demonstration of the proposed reaction force

The Liénard formula for the radiated power in terms of the four-velocity is
given by [1]

$P_{\mathrm{rad}}=-\tau_{m}m\left(\frac{du^{\mu}}{d\tau}\right)\left(\frac{du_{\mu}}{d\tau}\right),$
(35)

where $\tau$ is the proper time. If we assume that the external force is the
Lorentz force $F_{\mathrm{L}}^{\mu}=qF^{\mu\nu}u_{\nu}$, the equation of
motion taking into account the radiation is

$m\frac{du^{\mu}}{d\tau}=F_{\mathrm{L}}^{\mu}+F_{\mathrm{rad}}^{\mu}.$ (36)

Thus, the temporal component is the variation of the total energy per unit of
time

$\frac{dE}{dt}=e\mathbf{E}\cdot{\mathbf{v}}+\frac{cF_{\mathrm{rad}}^{0}}{\gamma},$
(37)

where the second term is the contribution due to the radiated power
$P_{\mathrm{rad}}$. Consequently, the temporal component of the reaction force
is

$F_{\mathrm{rad}}^{0}=-\frac{\gamma
P_{\mathrm{rad}}}{c}=-\frac{P_{\mathrm{rad}}}{c^{2}}u^{0},$ (38)

where the minus sign indicates that the particle is radiating energy. This
equation suggests that the reaction force is proportional to the four-
velocity, obtaining the following four-force

$F_{\mathrm{rad}}^{\mu}=-\frac{P_{\mathrm{rad}}}{c^{2}}u^{\mu}=\frac{\tau_{m}m}{c^{2}}\left(\frac{du^{\nu}}{d\tau}\right)\left(\frac{du_{\nu}}{d\tau}\right)u^{\mu}.$
(39)

If we assume that the reaction force is much less important than the Lorentz
force, i.e. $m\frac{du^{\mu}}{d\tau}=F_{\mathrm{L}}^{\mu}$, the vector
component of (39) is

$\mathbf{F}_{\mathrm{rad}}=-\frac{\tau_{m}v}{mc^{2}}\left(\gamma^{2}\|\mathbf{F}_{\mathrm{L},\perp}\|^{2}+\|\mathbf{F}_{\mathrm{L},\|}\|^{2}\right)\hat{\mathbf{v}},$
(40)

that coincides with the proposed reaction force (9) if $v\approx c$ is
assumed.

## References

  * Jackson [1999] J. D. Jackson,  _Classical electrodynamics, 3rd ed._ (Wiley, New York, NY, 1999).
  * Abraham [1905] M. Abraham,  _Theorie der elektrizität: Bd. Elektromagnetische theorie der strahlung_ , Vol. 2 (BG Teubner, 1905).
  * Lorentz [1909] H. A. Lorentz,  _The theory of electrons and its applications to the phenomena of light and radiant heat_ , Vol. 29 (BG Teubner, 1909).
  * Dirac [1938] P. A. M. Dirac, Classical theory of radiating electrons, Proc. R. Soc. A 167, 148 (1938).
  * Hartemann  _et al._ [2005] F. V. Hartemann, D. J. Gibson, and A. K. Kerman, Classical theory of Compton scattering: Assessing the validity of the Dirac-Lorentz equation, Phys. Rev. E 72, 026502 (2005).
  * Mo and Papas [1971] T. C. Mo and C. H. Papas, New equation of motion for classical charged particles, Phys. Rev. D 4, 3566 (1971).
  * Ford and O’Connell [1991] G. Ford and R. O’Connell, Radiation reaction in electrodynamics and the elimination of runaway solutions, Phys. Lett. A 157, 217 (1991).
  * Hartemann and Luhmann Jr [1995] F. Hartemann and N. Luhmann Jr, Classical electrodynamical derivation of the radiation damping force, Phys. Rev. Lett. 74, 1107 (1995).
  * Yaghjian [2010] A. Yaghjian,  _Relativistic dynamics of a charged sphere: updating the Lorentz-Abraham model_ , Vol. 11 (Springer Science & Business Media, 2010).
  * Hammond [2010] R. T. Hammond, Relativistic particle motion and radiation reaction in electrodynamics., Electron. J. Theor. Phys. 7 (2010).
  * Landau and Lifshitz [1975] L. D. Landau and E. M. Lifshitz,  _The Classical Theory of Fields_ , Course of Theoretical Physics (Pergamon Press, Oxford, 1975).
  * Medina [2006] R. Medina, Radiation reaction of a classical quasi-rigid extended particle, J. Phys. A 39, 3801 (2006).
  * Griffiths  _et al._ [2010] D. J. Griffiths, T. C. Proctor, and D. F. Schroeter, Abraham–Lorentz versus Landau–Lifshitz, Am. J. Phys. 78, 391 (2010).
  * Rohrlich [2001] F. Rohrlich, The correct equation of motion of a classical point charge, Phys. Lett. A 283, 276 (2001).
  * Nielsen  _et al._ [2021] C. Nielsen, J. Justesen, A. Sørensen, U. Uggerhøj, and R. Holtzapple, Experimental verification of the Landau–Lifshitz equation, New J. Phys. 23, 085001 (2021).
  * Hammond [2008] R. T. Hammond, Radiation reaction force without runaway solutions, Nuovo Cimento B 123, 567 (2008).
  * Burkhardt [2007] H. Burkhardt,  _Monte Carlo generation of the energy spectrum of synchrotron radiation_, Tech. Rep. (CERN, Geneva, 2007).
  * Bulanov  _et al._ [2011] S. V. Bulanov, T. Z. Esirkepov, M. Kando, J. K. Koga, and S. S. Bulanov, Lorentz-abraham-dirac versus landau-lifshitz radiation friction force in the ultrarelativistic electron interaction with electromagnetic wave (exact solutions), Phys. Rev. E 84, 056605 (2011).
  * Tajima and Cavenago [1987] T. Tajima and M. Cavenago, Crystal x-ray accelerator, Phys. Rev. Lett. 59, 1440 (1987).
  * Zhang  _et al._ [2016] X. Zhang, T. Tajima, D. Farinella, Y. Shin, G. Mourou, J. Wheeler, P. Taborek, P. Chen, F. Dollar, and B. Shen, Particle-in-cell simulation of x-ray wakefield acceleration and betatron radiation in nanotubes, Phys. Rev. Accel. Beams 19, 101004 (2016).
  * Huang and Kim [2007] Z. Huang and K.-J. Kim, Review of x-ray free-electron laser theory, Phys. Rev. ST Accel. Beams 10, 034801 (2007).
  * Smolyakov [2014] N. Smolyakov, Undulator Radiation Spectral Broadening Due To Radiation Energy Loss, in  _5th International Particle Accelerator Conference_ (2014) p. WEPRO041.