capbtabboxtable[][]

# Proactive slip control
by learned slip model and trajectory adaptation

Kiyanoush Nazari
School of Computer Science
University of Lincoln United Kingdom
<EMAIL_ADDRESS>
&Willow Mandil
School of Computer Science
University of Lincoln United Kingdom
<EMAIL_ADDRESS>
&Amir Ghalamzan E.
Lincoln Institute for Agri-Food Technology
University of Lincoln United Kingdom
<EMAIL_ADDRESS>

###### Abstract

This paper presents a novel control approach to dealing with object slip
during robotic manipulative movements. Slip is a major cause of failure in
many robotic grasping and manipulation tasks. Existing works increase grip
force to avoid/control slip. However, this may not be feasible when (i) the
robot cannot increase the gripping force– the max gripping force is already
applied or (ii) increased force damages the grasped object, such as soft
fruit. Moreover, the robot fixes the gripping force when it forms a stable
grasp on the surface of an object, and changing the gripping force during
real-time manipulation may not be an effective control policy. We propose a
novel control approach to slip avoidance including a learned action-
conditioned slip predictor and a constrained optimizer avoiding a predicted
slip given a desired robot action. We show the effectiveness of the proposed
trajectory adaptation method with receding horizon controller with a series of
real-robot test cases. Our experimental results show our proposed data-driven
predictive controller can control slip for objects unseen in training.

(Video and Github)

> Keywords: Slip avoidance, Robot Control, Manipulation

##  1 Introduction

Sense of touch in humans and primates is key to building sophisticated
dexterous manipulation skills [1]. It provides localized contact information
on human skin in particular as an effective perception in our daily lives [2].
A slip event occurs when the ratio of applied shear to the normal force on the
contact area exceeds the static friction coefficient. Recent studies show that
human reflexes to the feeling of slip (we use the term slip to refer to
slipping action of an object between fingers) during manipulative movements
not only include grip adjustment but also specific movements of the shoulder
and elbow [3]. Existing slip avoidance control methods in robotic manipulation
using artificial tactile sensing are limited to: (1) grip force adaptation in
case of slip [4, 5], (2) reactive algorithms that activate after slip onset
and detection [6, 7], and (3) gentle manipulation tasks, e.g., slow lifting
movements for 10-20 [cm] [5, 6]. However, increased grip force may not be an
effective action when (1) the gripper has already reached its force limits,
(2) the grasped object is delicate/fragile such as soft foods (i.e. strawberry
where the average grip force to damage the fruit is 1.5 N [8]). Moreover, many
grippers, such as Robotiq’s Hand-E Gripper or Franka Emika Gripper, do not
have enough force resolution, or (3) the gripper cannot be controlled in real-
time due to connection limitations i.e. blocking TCP/IP connections in the
Franka Emika gripper.

We propose two novel means of controlling slip by adapting the trajectory of
manipulative movements: (i) A Reactive Slip Control (RSC) policy, which uses
detected slip – slip classifier at time $t$– to reduce the robot velocity to
avoid slip; (ii) A Proactive Slip Control (PSC) policy for avoiding slip
during robotic manipulative movements (see schematic of the approach in Fig.
1). PSC includes an action-conditioned slip predictor, a slip classifier at
$t+h$ where $h$ is the prediction horizon using an LSTM (Long-Short Term
Memory) based neural network and a predictive controller. The controller
minimizes the likelihood of predicted slip in a prediction horizon by
adjusting the robot’s planned movements and deviating from the reference
trajectories. We use a quadratic term penalizing deviation from a reference
trajectory; hence, PSC avoids slip while closely following a given reference
trajectory. Our constrained optimization enables the robot to satisfy its
constraints in generating optimal future robot actions. We study the
effectiveness of our on-the-fly adaptation of a real robot’s pre-planned
trajectory (which may be minimum time or minimum jerk) in a series of real
robot experiments. We compare _RSC_ and _PSC_ in terms of their performance in
slip avoidance, computation complexity, and optimality. Our experimental
results indicate PSC outperforms RSC and can change the game in robotic
manipulation success where increased gripping force is not possible.

##  2 Related works

Slip detection algorithms are thoroughly reviewed in [9, 10]. [11] use Support
Vector Machine (SVM) and Random Forest (RF) on SynTouch BioTac’s raw tactile
readings and [12] use SVM on TacTip signals for slip classification. SVM and
RF have computational complexity problems in training for large data sets with
high feature dimensionality. Specifying a threshold on the rate of shear
forces for slip detection is a common approach for slip classification [13,
4]. However, it needs initial friction estimation, and finding the threshold
for large object sets can be expensive. Other techniques rely on analyzing
micro-vibration in the incipient slip phase such as spectral analysis [14,
15]. These approaches can fail when robot vibrations interfere with sensor-
object vibrations. Multi-modal approaches are also proposed to use proximity
[16] or visual sensing data [17] to improve slip detection. Recurrent Neural
Networks and Graph Neural Networks are also used for slip classification [18,
19]. While [18, 19] only use tactile readings for slip detection, we use an
LSTM-based model with combined tactile data and robot actions to predict slip
in a prediction horizon.

Grip force controllers based on tactile readings are usually designed to avoid
slip [11]. [20, 4, 7] increase the normal force until the ratio of normal to
shear force becomes smaller than the known/estimated static friction
coefficient. This method can lead to deformation or damage of delicate
objects, e.g., soft fruits. [4] propose a heuristic to limit deformations of a
deformable object but it needs to measure the initial grasp width and object
deformation. [21] suggest increasing the number of gripping fingers to avoid
slip by using a multi-finger hand. In bi-manual manipulation, [6] performs
hand pose correction to bring the resultant force inside the friction cone as
a slip avoidance controller. The above approaches have certain limitations
where increased grip force is not an option (e.g. in Franka Emika gripper, or
in delicate object manipulation).

Online trajectory re-planning is effectively used in areas such as agile
quadrotor control [22, 23], but to the authors’ knowledge, it is the first
time it is being used for slip avoidance in dynamic manipulation tasks. The
main difficulty relates to the lack of analytical models for tactile/slip
dynamics. By the proposed constrained optimization setting, we seek to benefit
from the data-driven models for slip prediction and model-based approaches for
online trajectory optimization.

Figure 1: A schematic of the proposed _proactive_ slip control (blue block)
that is running at 30 Hz and includes a slip prediction module (grey block),
and predictive controller for trajectory adaptation (orange block).

##  3 Methodology

Our proposed slip avoidance controllers are built in two stages: first,
offline slip classification phase where the LSTM-based networks (Section 3.1)
are trained on a collected dataset (Section 4) and second, online trajectory
adaptation phase where the trained models are used in a constraint
optimization setting (Section 3.2) to minimize the expected value of the slip
signal and deviation from a reference (desired) trajectory.

###  3.1 Slip detection and slip prediction models

Slip classification is explored by three main approaches in the literature:
(i) incipient slip detection [10], (ii) full slip detection [13, 4], and (iii)
slip prediction [5]. Incipient slip detection requires precise slip
calibration and sensor mechanical compliance [7, 24]. Our focus is using data-
driven approaches for full slip detection and slip prediction. Tactile shear
forces have proven to be a useful feature for slip classification [10]. This
grounds the higher rate of changes for shear force when the taxels are losing
contact with the objects, which is classified as an anomaly versus a stuck
phase.

[25, 26] show the benefits of using a history of tactile features and temporal
processing by modules such as LSTM cells. As such we use the shear forces of
the uSkin tactile sensor (introduced in Section 4) with two LSTM-based models
for slip _detection_ and _prediction_ problems. In _slip detection_ , a time
window of length $C$ of shear tactile forces
$\textbf{x}_{t-C:t}\in\mathbb{R}^{C\times 48}$ are used as input features to
classify the binary slip class $c_{t}=f(\textbf{x}_{t-C:t})$ where
$c_{t}\in\\{c_{slip},c_{non-slip}\\}$ at time step $t$. For _slip prediction_
, we use a multi-modal LSTM-based model (inspired by [11, 26]). This model
uses time windows of tactile features
$\textbf{x}_{t-C:t}\in\mathbb{R}^{C\times 48}$ and robot future actions
$\textbf{a}_{t+1:t+T}\in\mathbb{R}^{T\times 6}$ to classify the slip state
$c_{t+T}$ at time step $t+T$ (where $T$ is the prediction horizon). The
architecture of the _slip prediction_ model is shown in Fig. 1. The _slip
detection_ model has a similar architecture without the robot data input and
the corresponding dense layer. Fusing the latent vectors of two input
modalities seeks to exploit the future robot motion for slip prediction.

###  3.2 Closed-loop slip control by trajectory adaptation

The studies in closed-loop control design for slip avoidance are limited to
gentle short-range lifting motions [4, 20, 27]. Our proposed approach extends
slip control beyond simple slow lifting motions. First, we collected a data
set of dynamic moving tasks by pre-defined task space trajectories, i.e. a
linear motion with various velocity profiles and ranges of
acceleration/deceleration values. We learned slips are initiated by the
dynamic motions of the robot as opposed to being caused by object weight or a
non-stable grasp for slow lifting tasks in literature [4, 20, 27]. Therefore,
we propose two trajectory optimization/adaptation methods incorporating
tactile and robot proprioceptive data for slip control. We tested/compared PSC
and RSC approaches in real-time scenarios on real robot experimentation.

Reactive slip control by trajectory adaptation Inspired by episodic policy
search in [28], we represent the robot actions (e.g. task space velocity) by
linear Gaussian parametric models. As such, robot’s action at time $t$ can be
represented as $\dot{X}_{t}\;:=\sum_{j=1}^{N}w_{j}\phi_{j}(t)$, where
$w_{j}\in\mathbb{R}$ are the weight parameters and $\phi_{j}$s are the basis
functions with Gaussian forms
$\phi_{j}(t)=\exp{(-\frac{(t-\mu_{j})^{2}}{2\sigma})}$ where $j=1,2,..,N$ is
the number of basis functions.

A sequence of future robot actions with length $T$ in a matrix form will
become $\dot{X}_{\scriptstyle{t+1:t+T}}\;:=\Phi^{T}\times W$ where
$\Phi(t)=[\phi_{1}(t),\phi_{2}(t),...,\phi_{n}(t)]\in\mathbb{R}^{N\times T}$
and $W=\\{w_{1},w_{2},...,w_{n}\\}\in\mathbb{R}^{N}$ are the matrices of basis
functions and weight values respectively. The mean values of the Gaussians are
spread with equal distance in the horizon length $1:T$ and the $\sigma$
parameters are chosen such that Gaussians have reasonable widths and overlap.

The main policy in grip force adaptation for slip control is to incrementally
increase the normal force ($N(t)$) until the output of the slip detection
model ($S$) is zero [26]. This is usually formalized as
$N(t+1)=N(t)+S\times\delta N$, where $\delta N$ is the force increasing step.
Likewise, we implement RSC by adapting the robot movements trajectory by
decreasing task space velocity by fine steps (we learned slip happens during
the initial parts of movements with high acceleration) until the output of the
slip detection model is zero. If no slip is detected, the robot will follow
the desired reference trajectory. We empirically come up with a formulation in
Eq. 1 after testing different ones. This includes an affine combination of
slip cost and reference trajectory tracking cost. Although in RSC we do not
need optimizing over a receding horizon, this constrained optimization will be
later extended in PSC (in Eq. 2) allowing easy comparison between _PSC_ and
_RSC_.

$\displaystyle\underset{\textbf{w}}{\arg\min}\quad\;\lVert\dot{X}_{t+1:t+T}(\textbf{w},\phi)-\dot{X}^{ref}_{t+1:t+T}\rVert_{2}+a^{\frac{1}{S+\epsilon}}\lVert\dot{X}_{t+1:t+T}(\textbf{w},\phi)\lVert_{2}$
(1) $\displaystyle subject\;to:\quad
lb<\dot{X}_{t+1}-\dot{X}^{obs}_{t}<ub\qquad\qquad$

Where $\dot{X}$ is the generated trajectory by optimization, $\dot{X}^{ref}$
is the reference trajectory, $S\in\\{0,1\\}$ is the output of the slip
detection model introduced in 3.1, $a$ and $\epsilon\in\mathbb{R}^{+}$ are
hyper-parameters, and $\dot{X}^{obs}_{t}$ is robot’s observed velocity at time
$t$. The first term in the objective function aims to minimize the Euclidean
distance between the generated and the reference trajectory. The second term
in Eq. 1 results in generated smaller velocities (where the object is moved
with a large acceleration). As such, by choosing $a$ to be in $0<a<1$ and
$\epsilon$ as a small positive value i.e. 0.01, for non-slip mode when $S=0$,
the effect of the second term vanishes ($a^{\infty}\>,\>a\in(0,1)$) and if
$S=1$, it constantly reduces robot’s task space velocity until there is no
slip. The imposed constraint keeps the robot action at $t$ bounded to avoid
abrupt deceleration and ensures the admissibility condition of robot’s low-
level controller is met. We selected the exponential term for the slip signal
in that it resulted in more stable optimization and better slip avoidance
behavior than other forms, e.g., linear cost of slip. The hyper-parameters are
empirically chosen as $a=0.5$ and $\epsilon=0.01$.

Proactive slip control by trajectory adaptation Predictive controllers, e.g.
[29], showed success when dealing with nonlinear systems with the risk of
being stuck in a local minimum of the corresponding utility function. They can
achieve a solution close to the global optimum by finite-horizon optimization.
We use the same approach by solving the constrained optimization in a
prediction horizon and applying the first component of the generated action
sequence at each time step.

Parametric action representation has benefits in exploration complexity versus
direct search in action space [28]. The benefit is more effective when the
action is optimized in a future time horizon rather than a single time step.
In this regard, in the trajectory adaptation problem, we use linear Gaussian
parametric representation [30] for the robot action sequence. This leads to
fewer search parameters in the optimization problem relative to searching in
the action space. In this scenario, by the _slip prediction_ model introduced
in section 3.1, we close the loop for PSC by solving the constraint
optimization in Eq. 2. The predicted slip value at time $t$ is a function of
robot’s future trajectory, $S(\dot{X}_{t+1:t+T}(\textbf{w},\phi))$.

$\displaystyle\underset{\textbf{w}}{\arg\min}\;\lVert\dot{X}_{t+1:t+T}(\textbf{w},\phi)-\dot{X}^{ref}_{t+1:t+T}\rVert_{2}$
(2) $\displaystyle
subject\;to:\qquad\quad\;\operatorname{\mathbb{E}}[S(\dot{X}_{t+1:t+T}(\textbf{w},\phi))]\>=\>0\qquad$
$\displaystyle\qquad\qquad\quad lb<\dot{X}_{t+1}-\dot{X}^{obs}_{t}<ub\qquad$

Our PSC aims at minimizing the expected slip over the prediction horizon while
keeping close to the given reference trajectory as shown in Fig. 1. In
contrast to RSC in Eq. 1, we consider the slip term as a non-linear constraint
in Eq.2. This forms a more effective optimization for slip minimization since
the optimization library we use forms a Lagrangian function with optimizing
the Lagrange multipliers as opposed to the heuristic form of the multiplier in
the second term of Eq.1.

(a)

(b)

Figure 2: (a) Experimental setup including Franka Robot, 4x4 uSkin tactile
sensor on its fingers, box-shaped object, ArUco marker on the object, and
RealSense camera. A sample of trapezoidal task space velocity trajectory is
shown on the top left. (b) Data set features: tactile readings, robot
trajectories, marker trajectory, and slip signal.

##  4 Data set

###  4.1 Tactile sensor and manipulation task

Spatial resolution, reading frequency, and mechanical compliance are the main
features of a tactile sensor reviewed in [31]. We use uSkin, a magnetic-based
tactile sensor from XELA ROBOTICS that works with 100 Hz frequency. It
publishes 3-axis point-wise force measurements for 16 taxels arrayed in a 4x4
matrix (Fig. 2(a)). A pair of sensors are mounted on 3-D printed fingers for a
jaw gripper. We use a box-shaped object (4$\times$15$\times$25 [cm] and 400
[gr]) for training and to encourage more slip-cases, we reduce the friction
coefficient by covering the uSkin sensor with a _latex_ cover and increase the
object weight with large metal bolts fixed inside the box that do not move
during manipulation. A wrist camera reads the pose of the ArUco marker
attached on the top side of the object measuring the relative motion of the
object w.r.t the hand. This is used to create labels for ground truth slip
values. Maximum noise amplitude for marker readings is 0.38 [degrees] for
rotation and 0.5 [mm] for the position in this range of distance.

The robot grasps the object from a top grasp pose with the jaw gripper, lifts
it for 10 [cm] (in $+Z$ direction) and executes a linear motion in the
horizontal ($XY$) plane. The linear motion encourages slip events which are
direct outcome of dynamic manipulation forces rather than object static
weight. As such, for trajectories with small accelerations the object keeps it
upright orientation; While high accelerations cause object’s large rotation
(as a result of objects centre of mass moment of inertia around the grip axis)
and slip. We use task space velocity controller with several reference
trajectories with trapezoidal profile and various maximum velocity ($V_{MAX}$)
values. The velocity commands are converted to joint torques through the
robot’s low-level joint impedance controller. Fig. 2(a) shows a sample of the
reference velocity profile (start and end point of the motion are
$[0.4m,0.25m,0.3m]$ and $[0.1m,-0.25m,0.3m]$ respectively in robot’s base
frame); As $V_{MAX}$ changes, the acceleration and deceleration values change
accordingly such that the area under the curve which is equal to the overall
translation remains constant.

###  4.2 Data collection and pre-processing

We collected a data set of 660 linear motion tasks containing (1) tactile
readings from two fingers (at 100 [Hz]), (2) robot state data (at 1,000 [Hz])
including joints position, joints velocity, task space pose, and task space
velocity, (3) Marker pose data (at 60 [Hz]), and (4) Slip label values based
on marker’s rotation (see Fig. 2(b)). We use ROS ApproximateTime policy to
synchronize the data. We applied Kalman Filter on the raw tactile readings to
delete the measurement noise on the data (See Fig. 2(b)). Considering the
linear robot motion, slip can happen either as pure rotational slip, or
rotational slip followed by translational slip (which usually leads to object
dropping and task failure). As such, a threshold on the rotation of the ArUco
marker is used to create the slip labels. Based on marker rotations for slower
ranges of motions with no slip, rotations $>6^{\circ}$ are labeled as slip-
case. The data set along with preprocessing and models will be _publicly
available_ for future investigations of action-conditioned tactile or slip
prediction [32, 33] or slip detection models in a dynamic manipulation task.

(a)

(b)

(c)

(d)

(e)

Figure 3: Object rotation around grasp axis during following a reference
trajectory (a) without slip control and with (b) RSC, and (c) PSC, where
$>6^{\circ}$ are slip events. Trajectories generated by RSC (d) and PSC (e):
the blue areas show variations across 10 trails; the green bands show
approximate time slip is detected (predicted) for the first time.

##  5 Experimental results and discussion

We tested the proposed trajectory adaptation methods to control slip in real-
time in several test cases of manipulative trajectories with Franka Emika
robot connected to a PC with AMD® Ryzen 16-core processor and 64 GB RAM
running ROS Noetic. We use ONNX library to reduce the inference time of the
slip classification PyTorch models from 1 [ms] to 0.2 [ms]. Scipy optimization
library in python 3 (i.e. minimize() function with ’trust-constr’) is used to
solve the constrained optimization in Eq. 1 and 2. Each test trial is executed
10 times on the robot to achieve statistically meaningful results.

The metrics to compare the performance of slip controllers are (i) Resulted
optimality value (ROV) from the optimization library (i.e. gradient of the
Lagrangian function, and indicative of the convergence of the optimization),
(ii) Maximum amount of object rotation (MOR) in the trial, (iii) Execution
time (ET) in millisecond which shows the computation complexity, (iv)
instances of rotations larger than $6^{\circ}$, ”RTS” (Rotation $>6^{\circ}$
Time Steps) score which is the number of time steps that the object had
rotations larger than 6 degrees. Larger rotations and large number of cases of
RTS indicate higher probability of failure, and (v) distance to reference
trajectory (DRT) which is the sum of Euclidean distance between the reference
and adapted trajectories over the task execution time. In total, 300 testing
trials are executed with the train and test objects. Slip classification
accuracy and f-score on the test set (20% of the data) are 98% and 0.81 for
the detection and 94% and 0.74 for the prediction models, respectively. The
planning horizon $h$ is chosen by trading off between classification scores
and controller success rate. While longer horizons worsen the classification
performance, they give the proactive controller more reaction time and
increase its success rate by having a larger safety margin. By testing
horizons of length 5, 10, 15, and 20, we picked 10 as the best case. For the
context frames length $C$, considering the sensory data frequency (60 hz) and
task completion time (1.2-3 seconds), we chose $C=10$ as a reasonable value
including past 0.17 seconds tactile data.

Model # Basis Functions ROV MOR ET RTS DRT Reactive 2 0.80 $\pm$ 0.72 12.80
$\pm$ 7.58 9.79 $\pm$ 1.33 35 $\pm$ 8 1.31 $\pm$ 0.15 3 0.71 $\pm$ 0.56 9.06
$\pm$ 0.56 10.03 $\pm$ 1.57 31 $\pm$ 8 1.20 $\pm$ 0.15 4 0.67 $\pm$ 0.45 9.65
$\pm$ 1.50 10.49 $\pm$ 1.62 30 $\pm$ 5 1.16 $\pm$ 0.22 5 0.63 $\pm$ 0.31 8.82
$\pm$ 0.61 9.96 $\pm$ 1.49 28 $\pm$ 7 1.04 $\pm$ 0.23 6 0.47 $\pm$ 0.23 9.09
$\pm$ 1.19 10.15 $\pm$ 1.70 32 $\pm$ 10 1.05 $\pm$ 0.19 7 0.37 $\pm$ 0.22 9.66
$\pm$ 1.17 9.70 $\pm$ 2.66 33 $\pm$ 6 1.06 $\pm$ 0.17 8 0.28 $\pm$ 0.19 9.80
$\pm$ 0.75 9.99 $\pm$ 1.50 33 $\pm$ 10 0.99 $\pm$ 0.21 Proactive 2 0.75 $\pm$
0.60 3.99 $\pm$ 0.58 11.73 $\pm$ 5.35 0 $\pm$ 0 1.67 $\pm$ 0.05 3 0.72 $\pm$
0.67 6.14 $\pm$ 1.52 22.57 $\pm$ 1.79 4 $\pm$ 2 1.56 $\pm$ 0.24 4 0.51 $\pm$
0.41 7.16 $\pm$ 1.07 25.73 $\pm$ 2.14 5 $\pm$ 2 1.49 $\pm$ 0.27 5 0.40 $\pm$
0.34 9.06 $\pm$ 0.55 29.61 $\pm$ 3.26 8 $\pm$ 4 1.26 $\pm$ 0.20 6 0.33 $\pm$
0.22 9.75 $\pm$ 0.65 32.09 $\pm$ 5.05 13 $\pm$ 5 1.23 $\pm$ 0.18 7 0.30 $\pm$
0.22 9.41 $\pm$ 0.55 35.27 $\pm$ 3.65 11 $\pm$ 3 1.17 $\pm$ 0.36 8 0.25 $\pm$
0.20 9.14 $\pm$ 0.37 37.55 $\pm$ 3.84 10 $\pm$ 6 1.21 $\pm$ 0.21

Table 1: Trajectory optimization performance by different number of basis
functions. For all of the metrics, smaller values show better performance.
(all metrics are defined in Section 5)

Reactive versus proactive slip control with trajectory adaptation: As the
first test case, we compare RSC and PSC with the same reference trajectory and
design parameters. Resulted object rotation for RSC and PSC are presented in
Fig. 3(b) and Fig. 3(c) respectively; which can be compared to the case of
executing the reference trajectory without a slip control in Fig. 3(a)
(horizontal lines in Fig. 3(a) show the time instance the object is dropped
off the robot’s hand due to slip). RSC activates after slip is observed in
reality (i.e. rotation $>6^{\circ}$ around time step 70 in Fig. 3(b)), but it
keeps the rotation bounded in the rest of the trial. PSC (shown in Fig. 3(c))
outperforms RSC in the acceleration phase of the motion as the slip control
activates before slip actually initiates. Both controllers show 100% success
in avoiding task failure while in the mode without slip controller (Fig. 3(a))
19 trials out of 20 led to object dropping and failure.

The velocity profiles generated by RSC and PSC are shown in Figs. 3(d) and
3(e). Although these figures show RSC keeps closer than PSC to the reference
trajectory, it yields much more rotation of the object and none-smooth (with a
higher jerk that is undesirable) movements. PSC yields trajectories farther
from the reference trajectory but has smoother acceleration and deceleration
trends which can lead to smaller slip instances. This is mostly because the
weight of the slip avoidance relative to Euclidean distance term in Eq. 2 is
set by the Lagrange multipliers of the optimizer in PSC. This makes the
trajectories get more distant from the reference trajectory to satisfy the
slip constraint.

The effect of the number of basis functions on trajectory adaptation: We
tested _RSC_ and _PSC_ (_Reactive_ and _Proactive_ , respectively, in Table 1)
with number of basis functions ranging from 2 to 8. While PSC sees a
prediction of slip it has enough time to deviate enough from the reference to
avoid slip more effectively. This leads to having smaller RTS and larger DRT
for PSC. On the other hand, RSC has shorter reaction time relative to PSC,
which leads to larger RTS and smaller DRT. In terms of improving RTS (slip
prevention performance), maximum improvement in RSC by changing the number of
basis functions is $\frac{35-28}{35}\times 100=20\%$. While switching from RSC
to PSC with 5 basis functions for both, improves slip prevention by
$\frac{28-8}{28}\times 100=71\%$. This improvement can be higher for other
number of basis functions.

Depending on how many basis functions are used, the overlapping area of the
first Gaussian with the next ones can change (since we take only the first
element of the generated action sequence in the horizon, the change in the
beginning part of the horizon is most important). For RSC, this change has
positive effect on having smaller RTS up to 5 basis functions. For reference
tracking (DRT), increasing the number of basis functions improves tracking
performance. For PSC, 2 basis functions resulted zero RTS with largest
distance to the reference. Since the controller has fewer degrees of freedom
in this case, the cost it needs to pay for preventing slip is to get much
farther from the reference (the generated trajectories can be seen in Figure 2
(a) in the Appendix). This is a very conservative proactive controller and
based on how much weight we want to put on slip prevention w.r.t reference
tracking, we could prefer this controller. In fact, the number of basis
functions is a control parameter that we can use to specify how much
conservative we want to be about slip avoidance. The safety margin is achieved
by reducing the number of basis functions, but the cost is getting farther
from the reference trajectory. Changing the basis functions from 2 to 5 for
PSC, improves tracking performance by 6%, 4%, 15% respectively, but going
further to 6 basis functions has large negative effect on RTS
($\frac{13-8}{8}\times 100=62\%$).  Figure 4: Train and test objects. Object
ROV MOR RTS DRT theBest 0.51 $\pm$ 0.41 7.16 $\pm$ 1.07 4 $\pm$ 2 1.53 $\pm$
0.2 THINS 0.54 $\pm$ 0.43 7.88 $\pm$ 1.38 8 $\pm$ 3 1.71 $\pm$ 0.1 Carrs 0.52
$\pm$ 0.39 8.02 $\pm$ 2.12 7 $\pm$ 3 1.61 $\pm$ 0.3 RICE 0.61 $\pm$ 0.44 10.58
$\pm$ 3.56 10 $\pm$ 2 1.51 $\pm$ 0.2 Monster 0.66 $\pm$ 0.46 36.20 $\pm$ 13.93
28 $\pm$ 5 1.24 $\pm$ 0.3 Slug killer 0.52 $\pm$ 0.41 8.23 $\pm$ 1.86 6 $\pm$
2 1.62 $\pm$ 0.3 Table 2: Tests on novel objects. As such we choose 5 basis
functions to be the best value for PSC based on this quantitative analysis.
While ROV improves by an increased number of basis functions for both
controllers, MOR achieved by RSC does not show a clear improvement after 5
basis functions. In contrast, it leads to larger MOR for PSC. The computation
time increases proportionally in PSC, while it remains roughly constant in
RSC. The reason for larger computation time for PSC is that the Jacobians of
the slip constraint should be numerically estimated by the optimization
library which needs 30-50 calls of the slip detection/prediction networks in
each optimization iteration.

Generalization tests with novel objects: To evaluate generalization on novel
tactile features we test PSC with 5 basis functions on the novel objects shown
in Fig. 4 (from left to right the objects are called theBest, Slug killer,
Monster, RICE, THINS, and Carrs respectively). Three of the objects have a box
shape similar to the training object, and two are cylindrical for a more
difficult generalization cases with novel contact geometry. Table 2 shows the
resulted metrics for the train (first row) and the test objects. The first two
test objects have similar geometry and weight to the train object but with
different friction behavior where the PSC performance is close to the ones
with train object in terms of both RTS and MOR. Although the RICE box has
similar geometry, the PSC performance is relatively poor. This is partially
related to the lighter weight and different friction behavior of the box
cover. We observed our PSC can successfully avoid slip in manipulating objects
unseen in the training.

The Monster can with dynamic weight (liquid inside), deformable surface and
different geometry (cylindrical) and friction is the most challenging test
object for the controller and resulted in the worst performance metrics
values. The PSC can effectively control the Slug Killer can slip which has a
plastic and rigid material very different from the train object. The results
in Table 2 demonstrates PSC can generalize to a novel contact geometry for a
rigid object. In fact, the proposed constrained optimization setting is
dependent on how reliable the slip detection/prediction model performs. As
such, as far as the slip classification model can generalize to novel objects,
the proposed slip avoidance control is guaranteed to effectively avoid object
slip and task failure.

##  6 Conclusion

We presented a novel approach to slip control by online trajectory adaptation.
This is very important as the increased gripping force control policy is not
an effective solution in many settings. We presented two novel approaches for
slip control by robot trajectory adaptation (i) Reactive Slip Control and (ii)
Proactive Slip Control. We showed the success of the slip controllers in real-
time tests on a dynamic real-world manipulation task. We show Proactive Slip
Control outperforms the reactive one and it can generalize to objects unseen
in training. Our future works include extending our dataset to include
different tactile readings, e.g. camera-based tactile readings, and more real-
world manipulation experiments, i.e. different objects and manipulative
movements.

Limitations: Our approach is heavily relying on the slip data set to predict
the slip. We have only one class of trajectories (i.e. a path with different
velocity profiles) in our dataset which limits the prediction performance. We
will extend the dataset in our future works and collect data with different
classes of manipulative movement trajectories. The training dataset only
includes one object with varying weights. Nonetheless, the slip controller
could successfully generalize to unseen objects in our test cases. We will
improve our dataset by adding more objects (from multiple classes) to our
dataset in our future works. The controller is using one predicted slip signal
to adjust its next action. Using all the time steps in a prediction horizon
may improve the performance. In our future works, we will use a classifier for
different predicted tactile readings provided by an action-conditioned tactile
predictive model to improve the performance of slip control.

#### Acknowledgments

If a paper is accepted, the camera-ready version should include
acknowledgments. All acknowledgments go at the end of the paper, including
thanks to reviewers, to colleagues who contributed to the ideas, and to
funding agencies and corporate sponsors that provided financial support.

## References

  * Johansson and Flanagan [2008] R. Johansson and J. Flanagan.  6.05-tactile sensory control of object manipulation in humans.  _Senses Compr. Ref. Acad. Press NY NY_ , 67:86, 2008.
  * Dahiya et al. [2009] R. S. Dahiya, G. Metta, M. Valle, and G. Sandini.  Tactile sensing—from humans to humanoids.  _IEEE transactions on robotics_ , 26(1):1–20, 2009.
  * Hernandez-Castillo et al. [2020] C. R. Hernandez-Castillo, R. S. Maeda, J. A. Pruszynski, and J. Diedrichsen.  Sensory information from a slipping object elicits a rapid and automatic shoulder response.  _Journal of neurophysiology_ , 123(3):1103–1112, 2020.
  * Kaboli et al. [2016] M. Kaboli, K. Yao, and G. Cheng.  Tactile-based manipulation of deformable objects with dynamic center of mass.  In _2016 IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids)_ , pages 752–757. IEEE, 2016.
  * Veiga et al. [2020] F. Veiga, B. Edin, and J. Peters.  Grip stabilization through independent finger tactile feedback control.  _Sensors_ , 20(6):1748, 2020.
  * Hogan et al. [2020] F. R. Hogan, J. Ballester, S. Dong, and A. Rodriguez.  Tactile dexterity: Manipulation primitives with tactile feedback.  In _2020 IEEE international conference on robotics and automation (ICRA)_ , pages 8863–8869. IEEE, 2020.
  * Khamis et al. [2021] H. Khamis, B. Xia, and S. J. Redmond.  Real-time friction estimation for grip force control.  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 1608–1614. IEEE, 2021.
  * Hietaranta and Linna [1999] T. Hietaranta and M.-M. Linna.  Penetrometric measurement of strawberry fruit firmness: device testing.  _HortTechnology_ , 9(1):103–105, 1999.
  * Romeo and Zollo [2020] R. A. Romeo and L. Zollo.  Methods and sensors for slip detection in robotics: A survey.  _IEEE Access_ , 8:73027–73050, 2020.
  * Chen et al. [2018] W. Chen, H. Khamis, I. Birznieks, N. F. Lepora, and S. J. Redmond.  Tactile sensors for friction estimation and incipient slip detection—toward dexterous robotic manipulation: A review.  _IEEE Sensors Journal_ , 18(22):9049–9064, 2018\.
  * Veiga et al. [2015] F. Veiga, H. Van Hoof, J. Peters, and T. Hermans.  Stabilizing novel objects by learning to predict tactile slip.  In _2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pages 5065–5072. IEEE, 2015.
  * James et al. [2018] J. W. James, N. Pestell, and N. F. Lepora.  Slip detection with a biomimetic tactile sensor.  _IEEE Robotics and Automation Letters_ , 3(4):3340–3346, 2018.
  * Su et al. [2015] Z. Su, K. Hausman, Y. Chebotar, A. Molchanov, G. E. Loeb, G. S. Sukhatme, and S. Schaal.  Force estimation and slip detection/classification for grip control using a biomimetic tactile sensor.  In _2015 IEEE-RAS 15th International Conference on Humanoid Robots (Humanoids)_ , pages 297–303. IEEE, 2015.
  * Van Wyk and Falco [2018] K. Van Wyk and J. Falco.  Calibration and analysis of tactile sensors as slip detectors.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 2744–2751. IEEE, 2018.
  * Romano et al. [2011] J. M. Romano, K. Hsiao, G. Niemeyer, S. Chitta, and K. J. Kuchenbecker.  Human-inspired robotic grasp control with tactile sensing.  _IEEE Transactions on Robotics_ , 27(6):1067–1079, 2011.
  * Hasegawa et al. [2010] H. Hasegawa, Y. Mizoguchi, K. Tadakuma, A. Ming, M. Ishikawa, and M. Shimojo.  Development of intelligent robot hand using proximity, contact and slip sensing.  In _2010 IEEE International Conference on Robotics and Automation_ , pages 777–784. IEEE, 2010.
  * Li et al. [2018] J. Li, S. Dong, and E. Adelson.  Slip detection with combined tactile and visual information.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 7772–7777. IEEE, 2018.
  * Deng et al. [2020] Z. Deng, Y. Jonetzko, L. Zhang, and J. Zhang.  Grasping force control of multi-fingered robotic hands through tactile sensing for object stabilization.  _Sensors_ , 20(4):1050, 2020.
  * Garcia-Garcia et al. [2019] A. Garcia-Garcia, B. S. Zapata-Impata, S. Orts-Escolano, P. Gil, and J. Garcia-Rodriguez.  Tactilegcn: A graph convolutional network for predicting grasp stability with tactile sensors.  In _2019 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–8. IEEE, 2019.
  * Dong et al. [2019] S. Dong, D. Ma, E. Donlon, and A. Rodriguez.  Maintaining grasps within slipping bounds by monitoring incipient slip.  In _2019 International Conference on Robotics and Automation (ICRA)_ , pages 3818–3824. IEEE, 2019.
  * Kobayashi et al. [2014] F. Kobayashi, H. Kanno, H. Nakamoto, and F. Kojima.  Slip based pick-and-place by universal robot hand with force/torque sensors.  In _2014 IEEE Symposium on Robotic Intelligence in Informationally Structured Space (RiiSS)_ , pages 1–4. IEEE, 2014.
  * Romero et al. [2022] A. Romero, R. Penicka, and D. Scaramuzza.  Time-optimal online replanning for agile quadrotor flight.  _arXiv preprint arXiv:2203.09839_ , 2022.
  * Zhou et al. [2021] B. Zhou, J. Pan, F. Gao, and S. Shen.  Raptor: Robust and perception-aware trajectory replanning for quadrotor fast flight.  _IEEE Transactions on Robotics_ , 37(6):1992–2009, 2021.
  * Yuan et al. [2015] W. Yuan, R. Li, M. A. Srinivasan, and E. H. Adelson.  Measurement of shear and slip with a gelsight tactile sensor.  In _2015 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 304–311. IEEE, 2015.
  * Zapata-Impata et al. [2019] B. S. Zapata-Impata, P. Gil, and F. Torres.  Learning spatio temporal tactile features with a convlstm for the direction of slip detection.  _Sensors_ , 19(3):523, 2019.
  * Veiga et al. [2018] F. Veiga, J. Peters, and T. Hermans.  Grip stabilization of novel objects using slip prediction.  _IEEE transactions on haptics_ , 11(4):531–542, 2018.
  * Zhang et al. [2020] Y. Zhang, W. Yuan, Z. Kan, and M. Y. Wang.  Towards learning to detect and predict contact events on vision-based tactile sensors.  In _Conference on Robot Learning_ , pages 1395–1404. PMLR, 2020.
  * Plappert et al. [2017] M. Plappert, R. Houthooft, P. Dhariwal, S. Sidor, R. Y. Chen, X. Chen, T. Asfour, P. Abbeel, and M. Andrychowicz.  Parameter space noise for exploration.  _arXiv preprint arXiv:1706.01905_ , 2017.
  * Rawlings et al. [2017] J. B. Rawlings, D. Q. Mayne, and M. Diehl.  _Model predictive control: theory, computation, and design_ , volume 2.  Nob Hill Publishing Madison, 2017.
  * Deisenroth et al. [2013] M. P. Deisenroth, G. Neumann, J. Peters, et al.  A survey on policy search for robotics.  _Foundations and trends in Robotics_ , 2(1-2):388–403, 2013.
  * Yamaguchi and Atkeson [2019] A. Yamaguchi and C. G. Atkeson.  Recent progress in tactile sensing and sensors for robotic manipulation: can we turn tactile sensing into vision?  _Advanced Robotics_ , 33(14):661–673, 2019.
  * Mandil et al. [2022] W. Mandil, K. Nazari, et al.  Action conditioned tactile prediction: a case study on slip prediction.  _arXiv preprint arXiv:2205.09430_ , 2022.
  * Nazari et al. [2021] K. Nazari, W. Mandill, M. Hanheide, and A. G. Esfahani.  Tactile dynamic behaviour prediction based on robot action.  In _Annual Conference Towards Autonomous Robotic Systems_ , pages 284–293. Springer, 2021.


# PAO: A general particle swarm algorithm with exact dynamics and closed-form
transition densities

Max D. Champneys
Dynamics Research Group
University of Sheffield
Mappin St, Sheffield, S1 4ET
<EMAIL_ADDRESS>
& Timothy J. Rogers
Dynamics Research Group
University of Sheffield
Mappin St, Sheffield, S1 4ET
<EMAIL_ADDRESS>
Corresponding author

(April, 2023)

###### Abstract

A great deal of research has been conducted in the consideration of meta-
heuristic optimisation methods that are able to find global optima in settings
that gradient based optimisers have traditionally struggled. Of these, so-
called particle swarm optimisation (PSO) approaches have proven to be highly
effective in a number of application areas. Given the maturity of the PSO
field, it is likely that novel variants of the PSO algorithm stand to offer
only marginal gains in terms of performance—there is, after all, no free
lunch. Instead of only chasing performance on suites of benchmark optimisation
functions, it is argued herein that research effort is better placed in the
pursuit of algorithms that also have other useful properties. In this work, a
highly-general, interpretable variant of the PSO algorithm —particle attractor
algorithm (PAO)— is proposed. Furthermore, the algorithm is designed such that
the transition densities (describing the motions of the particles from one
generation to the next) can be computed exactly in closed form for each step.
Access to closed-form transition densities has important ramifications for the
closely-related field of Sequential Monte Carlo (SMC). In order to demonstrate
that the useful properties do not come at the cost of performance, PAO is
compared to several other state-of-the art heuristic optimisation algorithms
in a benchmark comparison study.

_Keywords_ Meta-heuristic optimisation $\cdot$ Particle swarm optimisation
$\cdot$ Stochastic differential equations

##  1 Introduction

A great number of challenges in science, engineering and beyond can be cast as
optimisation problems. Indeed, a great deal of research effort has been
expended in the specification of methods that are able to provide optimal
parameters for a given objective function. Of particular interest are
approaches that can perform the optimisation without access to gradient
information. So called _heuristic methods_ have attracted a great deal of
research interest over many decades. Already, countless methods have been
proposed, ranging from early genetic algorithms [1] to more recent
evolutionarily inspired approaches. Recent comprehensive reviews of the field
can be found in [2, 3]. Despite strong empirical performance, the no-free-
lunch theorem for optimisation [4] precludes the existence of a single ‘best’
approach. For any two optimisation algorithms it is always possible to find an
optimisation task upon which one algorithm outperforms the other (and _vice
versa_).

Although no algorithm can be proven to be better than any other, this does not
mean that new approaches should not be sought out. Instead of placing the
focus solely on benchmark performance, it is useful to consider other factors
that make some heuristic methods more useful in practice. In this work, a
novel methodology is proposed that incorporates both _interpretable
hyperparameters_ , _exact dynamics_ and _closed-form transition densities_.

One issue (amongst many [5]) with several heuristic and meta-heuristic
optimisation strategies is that the performance of the approach comes to
depends strongly on the choice of several hyperparameters. In many
optimisation settings there exists a trade-off between so-called exploration
(sampling the search space for promising minima) and exploitation (greedily
searching within a minima for an optima) behaviours. Hyperparameter choice can
play a key role in biasing the search to one mode of optimisation or the
other. A motivation of the current work is that this trade-off is made as
interpretable as possible.

A further motivation for the approach presented in this work is direct closed-
form access to the transition densities of the particles. The position of any
particle in the next generation can be expressed as by a Gaussian distribution
that depends only on the previous value. This is desirable as it enables the
optimisation process to be used as a proposal step within a sequential Monte
Carlo (SMC) scheme [6] in order to provide robust uncertainty quantification
in the presence of noisy objective functions.

###  1.1 Background and related work

Optimisation algorithms inspired by nature and evolution have been a a topic
of considerable research since the early twentieth century. Although many
algorithms have been proposed, there has been recent concern that many of
these offer at best incremental novelty and at worst a deliberate deception
due to inherent bias towards optimisation benchmark functions where the optima
lie at the centre of the search space [5]. Although many new algorithms are
published every year, many application papers rely on a core number of highly-
cited methods and their variants including ant colony optimisation [7],
particle swarm optimisation (PSO) [8] and differential evolution (DE) [9].

Among the most intensively researched heuristics are those belonging to the
particle swarm class. First proposed by Kennedy and Eberhart in [8], the
approach has been extensively studied and extended. Applications of particle
swarm methods can be found in diverse fields from power systems [10] to image
segmentation [11]. Since the original proposition, a great deal of variants to
the method have been proposed, including additive noise [12], quantum leaps
[13], Gaussian attractors [14], ordinary differential equation attractors in
[15] and multiple constraints in [16] as well as many others.

In this study, the PSO algorithm is treated as a stochastic differential
equation (SDE). Casting the particle swarm as an SDE is not in itself a novel
approach. In a series of papers Grassi et al. [17, 18] cast several variants
of the particle swarm algorithm as SDEs in order to derive convergence results
from mean-field limit theory. However, the derivations in [17, 18] result in
nonlinear SDEs that cannot be solved exactly in order to recover closed-form
transition densities.

###  1.2 Contribution

Although particle swarms have been separately cast as SDEs [17, 18] and
employed with additive stochastic terms [12], the authors are not presently
aware of any methods that have employed both of these approaches with the
specific aim of providing a particle swarm algorithm with a closed-form
transition density.

In this paper, the authors propose a variant PSO approach – the particle
attractor optimisation (PAO)111Pronounced ‘pa’. algorithm. In the proposed
approach, the motions of each element of each particle are given by the
solution of a scalar, stochastic differential equation with additive noise.
Although the formulation bears a close resemblance to existing approaches, the
linear stochastic formulation here comes with a number of distinct advantages
over existing particle swarm methods:

  * •

The dynamics of the particles are computed exactly and no difference-equation
approximation is required.

  * •

The method is highly flexible through the choice of the attractors.

  * •

The hyperparameters of the approach are interpretable as the parameters of a
dynamic system (damping ratio, natural frequency etc.).

  * •

The forward transition densities of the particles are available in closed
form.

##  2 Methodology

In this work, a generalised version of a particle swarm optimisation algorithm
is proposed that enables the particle dynamics from one iteration to the next
to be represented as a linear, time-invariant SDE with additive stochasticity.
As in other PSO algorithms, the motions of the particles are driven by
restoring forces that act in the direction of a number of attractors in the
search space. The positions of the attractors are user-defined and updated
each iteration. Choices for the attraction points might include:

  * •

_Global best_ : The best particle location found across all particles.

  * •

_Local best_ : The best location in the location history of a single particle.

  * •

_Average local best_ : The local best positions averaged across all
particles222Care must be taken to avoid introducing centre-bias to the
optimisation..

  * •

_Average particle_ : The average particle location.

  * •

_Weighted average particle_ : The average particle location weighted by the
fitness scores.

  * •

Attractors based on other meta-heuristic methods i.e ‘rand/1/bin’ from
differential evolution [9].

  * •

Stochastic attractors such as those in [14].

However, the formulation of the approach here does not depend on the choice of
attractors. In order to ensure that the underlying SDE has an exact solution,
stochasticity is introduced to the dynamics by the application of independent
forcing terms that are applied to each element of each particle independently
and additively. The noise process is modelled as a zero-mean Gaussian with
variance given by,

$\sigma^{2}=q\nu(\alpha)$ (1)

where $\nu$ is a user-defined function that depends only on the positions of
the attractors, scaled by the hyperparameter $q_{0}$. In this work,
$\nu(\alpha)$ is defined as the sum of the squared distances between the
average and global best particle locations. The reasoning here is that the
amount of noise should decrease as these attraction centres converge onto each
other. The forces acting on each particle are depicted in Figure 1.

$\bm{x}$$\epsilon\sim\mathcal{N}(0,\sigma^{2})$$m\ddot{\bm{x}}+c\dot{\bm{x}}$$\alpha_{1}$$\alpha_{2}$$\alpha^{\prime}$$k_{1}$$k_{2}$$k^{\prime}\bm{x}^{\prime}$$\Delta\bm{x}$
Figure 1: Visualisation of the motions of a single particle in the proposed
algorithm. The particle moves as if attached to linear elastic springs at each
attraction point. The movement of the particle is subject to inertial and
viscosity terms that resist motion. Stochasticity is introduced as a Gaussian
white noise excitation term. The red arrow represents the restoring force from
equivalent single spring that acts towards the overall attractor weighted by
the stiffness terms.

Mathematically, the motions of the particles can be described by a second-
order differential equation. In order to simplify the forthcoming notation,
let, $x=x_{ij}$ be the position of the $j$th element of the $i$th particle in
the search space. Similarly, let, $\alpha_{r}=\alpha_{rij}$ be the position of
the $j$th element of the $i$th particle for the $r$th attractor. Because the
dynamics of the motions of the elements of the particles are all independent
from each other, they can be individually represented by the second-order
scalar differential equation,

$m\ddot{x}+c\dot{x}+\sum_{r}k_{r}(x-\alpha_{r})+w(t)$ (2)

where $m$, $c$ and $k_{i}$ are hyperparameters relating to the inertia,
viscous drag and relative weights of the attractors $\alpha_{r}$ and each
over-dot represents a time derivative. $w(t)$ is a continuous time white-noise
process with zero mean and spectral density given by $\sigma^{2}$.

The equation of motion above can be trivially recast into a first-order state-
space form as,

$\dot{\bm{x}}=\begin{bmatrix}\dot{x}^{\prime}\\\
\dot{z}\end{bmatrix}=\begin{bmatrix}z\\\
\frac{1}{m}(-cz-k^{\prime}x^{\prime}-\epsilon)\end{bmatrix}=\bm{0}$ (3)

where the following coordinate transformations have been made,

$k^{\prime}=\sum_{r}k_{r}$ (4)

$x^{\prime}=x-\frac{1}{k^{\prime}}\sum_{r}k_{r}\alpha_{r}$ (5)

$z=\dot{x}^{\prime}$ (6)

To aid with the interpretability of the hyperparameters, it is convenient to
express the viscous damping term in the above in terms of the damping ratio
[19],

$\zeta=\frac{c}{2\sqrt{k^{\prime}m}}$ (7)

The above differential equation admits an Itô process representation of the
form,

$\mathop{}\mathrm{d}\bm{x}=F\bm{x}\mathop{}\mathrm{d}t+\bm{L}\mathop{}\mathrm{d}\beta$
(8)

whereby

$F=\begin{bmatrix}0&1\\\
-\frac{k^{\prime}}{m}&-2\sqrt{\frac{k^{\prime}}{m}}\zeta\end{bmatrix}$ (9)

$\bm{L}=\begin{bmatrix}0\\\ 1\end{bmatrix}$ (10)

and $\beta$ is the formal Brownian motion process in one dimension with
diffusion coefficient given by $Q=\sigma^{2}$ (the spectral density of the
driving white noise process). If one assumes that all $m$, $c$, $k^{\prime}$,
$\alpha_{r}$ and $q$ are constant in the interval $[t_{0},t]$,333This is an
implicit assumption in the discretisation of the dynamics in practically every
particle swarm optimisation method. then the solution to (8) is a standard
result in the SDE literature [20]. In order to keep this paper self contained,
the salient results are reproduced here. The general solution to (8) is given
by,

$\mathbf{x}(t)=\exp\left(F\left(t-t_{0}\right)\right)\bm{x}\left(t_{0}\right)+\int_{t_{0}}^{t}\exp(F(t-\tau))\bm{L}\mathrm{d}\beta(\tau)$
(11)

where $\exp(\cdot)$ represents the matrix exponential operation. Denoting now,

$A=\exp\left(F\Delta t\right)=\exp\left(F\left(t-t_{0}\right)\right)$ (12)

then the time evolution of each element of each particle is given by,

$\bm{x}_{t}=A\bm{x}_{0}+\int_{0}^{\Delta t}A\bm{L}\mathrm{d}\beta$ (13)

where $\bm{x}_{0}$ is a given (deterministic) initial condition. The mean and
covariance of the above can now be computed,

$\mathbb{E}[\bm{x}_{t}]=A\bm{x}_{0}$ (14)

$\mathbb{E}[\bm{x}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{x}_{t}]=\int_{0}^{\Delta
t}\exp\left(F\left(\Delta
t-\tau\right)\right)\bm{L}Q\bm{L}^{\mkern-1.5mu\mathsf{T}}\exp\left(F\left(\Delta
t-\tau\right)\right)^{\mkern-1.5mu\mathsf{T}}\mathrm{d}\beta$ (15)

$\mathbb{E}[\bm{x}_{t}^{\mkern-1.5mu\mathsf{T}}\bm{x}_{t}]=\Sigma=\bm{L}Q\bm{L}^{\mkern-1.5mu\mathsf{T}}-A\bm{L}Q\bm{L}^{\mkern-1.5mu\mathsf{T}}A^{\mkern-1.5mu\mathsf{T}}$
(16)

Which implies that the dynamics for each element of each particle can be
updated according to,

$\bm{x}_{t+1}=A\bm{x}_{t}+\bm{d},\quad\bm{d}\sim\mathcal{N}(\bm{0},\Sigma)$
(17)

The transition probability from $\bm{x}_{t}$ to $\bm{x}_{t+1}$ can therefore
be written in the form a a Gaussian distribution,

$p(\bm{x}_{t+1}|\bm{x}_{t})=\mathcal{N}(A\bm{x}_{t},\Sigma)$ (18)

A naïve implementation of the proposed algorithm is sketched in 1.

Algorithm 1 Particle attractor optimisation (PAO)

Number of particles $N$, hyperparameters $\\{m,\zeta,k_{i},q,\Delta t\\}$

Initialise an initial swarm of $N$ particle vectors
$\bm{x}_{j}\in\mathbb{R}^{D}$

Initialise the initial velocities of the $N$ particle vectors.

Compute the value of the objective function for each particle.

for each generation do

for each particle do

Compute the positions attraction centres $\alpha_{r}$.

Compute the position of the particle in the transformed coordinates as
$\bm{x}^{\prime}=\bm{x}-\frac{1}{k^{\prime}}\sum_{r}k_{r}\alpha_{r}$.

Compute the mean and the covariances for the transition density as above.

Sample a new particle from the $\mathcal{N}(A\bm{x}_{t},\Sigma)$

Recover the positions of the new particle in the original coordinates as
$\bm{x}=\bm{x}^{\prime}-\frac{1}{k^{\prime}}\sum_{r}k_{r}\alpha_{r}$.

Compute the value of the objective function for the new particle.

end for

end for

### Notes on implementation

Since the positions of the global and local best particles change in each
iteration (and therefore the value of $Q$), it is necessary to recompute
$\Sigma$ each time the swarm moves (Note that $A$, and $Q$ only come to depend
on fixed hyperparameters.). Unfortunately, the relation,

$\Sigma=\bm{L}Q\bm{L}^{\mkern-1.5mu\mathsf{T}}-A\bm{L}Q\bm{L}^{\mkern-1.5mu\mathsf{T}}A^{\mkern-1.5mu\mathsf{T}}$
(19)

is prone to numerical instability as both terms in the difference are similar
in magnitude. Fortunately, the method of _matrix fraction decomposition_ can
be employed to jointly calculate $A$ and $\Sigma$ employing only a single
matrix exponential. The interested reader is directed to [20] for a full
treatment, but the salient result is,

$\exp(\Phi\Delta
t)=\exp\left(\begin{bmatrix}F&\bm{L}Q\bm{L}^{\mkern-1.5mu\mathsf{T}}\\\
\bm{0}&-F^{\mkern-1.5mu\mathsf{T}}\end{bmatrix}\Delta
t\right)=\begin{bmatrix}A&\Sigma(A^{-1})^{\mkern-1.5mu\mathsf{T}}\\\
\bm{0}&(A^{-1})^{\mkern-1.5mu\mathsf{T}}\end{bmatrix}$ (20)

Thus, $\Sigma$ can be read off as the upper right block of $\exp(\Phi\Delta
t)$ post-multiplied by the upper left block.

Computing the matrix exponential every iteration still adds undesirable
computational complexity to the method. In order to significantly alleviate
this, one notices that the covariance structure is independent for every
element of every particle and can therefore be extracted as a constant in the
computation. If one pulls out Q as a factor from the above, the evolution of
the dynamics can thus be written,

$\bm{x}_{t+1}=A\bm{x}_{t}+q_{0}\nu(\alpha)\bm{d},\quad\bm{d}\sim\mathcal{N}(\bm{0},\Sigma)$
(21)

where $A$ and $\Sigma$ can be precomputed. Further numerical advantage can be
achieved by pre-computing the Cholesky decomposition of $\Sigma$,

$\Sigma=HH^{\mkern-1.5mu\mathsf{T}}$ (22)

so that the above can be written,

$\bm{x}_{t+1}=A\bm{x}_{t}+\sqrt{q_{0}\nu(\alpha)}H\bm{d},\quad\bm{d}\sim\mathcal{N}(0,I_{2})$
(23)

The overall computation (for every element of every particle) can be
efficiently implemented on modern architecture using tensor operations. The
proposed algorithm is sketched in algorithm 2. In the notation below, $A\odot
B$ refers to a broadcasted matrix-vector multiplication between the last two
indices of the tensor $A$ and the last axis of the Tensor $B$. $A\otimes B$
represents a tensor outer product.

Algorithm 2 An efficient implementation of PAO

Number of particles $N$, hyperparameters $\\{m,\zeta,k_{i},q,\Delta t\\}$

Precompute transition matrix $A\in\mathbb{R}^{2\times 2}$ and covariance
$\Sigma\in\mathbb{R}^{2\times 2}$ for a single element of a single particle.

Take the Cholesky decomposition of $\Sigma$ to compute
$H\in\mathbb{R}^{2\times 2}$.

Generate $N$ such particle vectors of position $\bm{x}_{j}\in\mathbb{R}^{D}$
and velocity $\dot{\bm{x}}_{j}\in\mathbb{R}^{D}$ according to some
initialisation scheme. Assemble these vectors into the tensor
$X_{ijk}\in\mathbb{R}^{N\times D\times 2}$.

Compute the value of the objective function for each particle in the tensor
$X$.

for each generation do

Compute the $r$ attraction centres $\bm{\alpha}\in\mathbb{R}^{r\times N\times
D}$ as required.

Transform the positions of the particles by
$\bm{x}_{ij1}^{\prime}=\bm{x}_{ij1}-\frac{1}{k^{\prime}}\sum_{r}k_{r}\alpha_{rij}$.

Move the swarm by $X^{\prime}_{t+1}=A\odot
X^{\prime}_{t}+(q\nu(\bm{\alpha})\otimes H)\odot D,\quad
D\in\mathbb{R}^{N\times D\times 2},\quad d_{ijk}\sim\mathcal{N}(0,1)$.

Recover the positions of the particles in the original coordinates as
$\bm{x}_{ij1}=\bm{x}^{\prime}_{ij1}+\frac{1}{k^{\prime}}\sum_{r}k_{r}\alpha_{rij}$.

Compute the value of the objective function for each particle in the tensor
$X$.

end for

##  3 Benchmark comparison study

In order to demonstrate that the proposed approach is competitive with other
approaches a simple benchmark study is presented. Although it is not the
primary objective of this algorithm to produce a heuristic method that out-
performs state-of-the art methods, it is important to establish that the
proposed method remains effective at the task of optimisation and that the
benefits of interpretable hyperparameters, exact dynamics and closed form
transition densities do not come at the cost of optimisation performance.

In order to assess the performance of PAO, a suite of nine standard
optimisation benchmark problems are considered. The problems are selected in
order to demonstrate the effectiveness of the model in a number of challenging
scenarios. The considered problems are collected in Table 1 in approximate
order of difficulty.

Table 1: Benchmark optimisation functions considered in this study. Problem | Domain | $f(\bm{x})$ | $\bm{x}^{*}$
---|---|---|---
De Jong’s function | $x_{i}\in[-5.12,5.12]$ | $\sum_{i}^{N}x_{i}^{2}$ | $x_{i}=0$
Hyper-ellipsoid function | $x_{i}\in[-5.12,5.12]$ | $\sum_{i}^{N}ix_{i}^{2}$ | $x_{i}=0$
Rotated Hyper-ellipsoid | $x_{i}\in[-65.54,65.54]$ | $\sum_{i}^{n}\sum_{j}^{i}x_{j}^{2}$ | $x_{i}=0$
Power-sum function | $x_{i}\in[-1,1]$ | $\sum_{i}^{n}|x_{i}|^{i+1}$ | $x_{i}=0$
Rosenbrock function | $x_{i}\in[-2.048,2.048]$ | $\sum_{i}^{n}[100(x_{i+1}-x_{i}^{2})^{2}+(1-x_{i})^{2}]$ | $x_{i}=1$
Griewangk’s function | $x_{i}\in[-600,600]$ | $\frac{1}{400}\sum_{i}^{n}x_{i}^{2}-\prod_{i}^{n}\cos(\frac{x_{i}}{\sqrt{i}})+1$ | $x_{i}=0$
Rastrigin’s function | $x_{i}\in[-5.12,5.12]$ | $10n+\sum_{i}^{n}[x_{i}^{2}-10*\cos(2\pi x_{i})]$ | $x_{i}=0$
Ackley’s function | $x_{i}\in[-32.77,32.77]$ | $-20e^{-0.2\sqrt{\frac{1}{n}\sum_{i}^{n}x_{i}^{2}}}-e^{\frac{1}{n}\sum_{i}^{n}cos(2\pi x_{i})}+20+e$ | $x_{i}=0$
Schwefel’s function | $x_{i}\in[-500,500]$ | $\sum_{i}^{n}[-x_{i}\sin(\sqrt{|x_{i}|})]$ | $x_{i}\approx 420.97$

For each benchmark function, both two and eight dimensional problems are
considered. The performance of the proposed algorithm (PAO) is compared to
several very popular optimisation approaches that have proven effective in a
wide range of fields:

  * •

Particle swarm optimisation (PSO) [8].

  * •

Quantum particle swarm optimisation (QPSO) [13].

  * •

Differential Evolution (DE) [9].

  * •

Self-adapting differential evolution (SADE) [21].

All of the above algorithms (including PAO) are implemented with default
hyperparameters from the freelunch python package developed by the
authors.444Implementations of all algorithms used in this study are freely
available as part of the open-source optimisation library ‘FreeLunch’,
available at https://pypi.org/project/freelunch/. To best highlight the
performance of the PAO approach, the attractors used in this study are simply
the local best and global best particle locations as in the standard PSO
algorithm. For each of the nine benchmark functions, each optimiser was
randomly initialised 100 times and run for 100 generations with a population
size of 100. In all cases this corresponds to 100 runs of each optimiser with
a budget of $1.01\times 10^{3}$ function evaluations. Figures 2 and 3 depict
convergence histories averaged over the 100 runs for each algorithm on each
optimisation problem for the 2D and 8D suites respectively. To aid comparison,
the convergence histories in the figures are shifted such that value of the
objective function at the true optimum lies at zero, i.e. $f(\bm{x}^{*})=0$.
For the convenience of the reader, all parameters relating to the benchmark
study and PAO are collected in Table 2

Table 2: Parameters relating to the benchmark comparison study. Parameter | Description | Value
---|---|---
$m$ | PAO inertia coefficient | 1.0
$\zeta$ | PAO damping ratio | 0.2
$k_{r}$ | PAO stiffness parameters | 1.0
$q_{0}$ | PAO stochastic parameter | 1.0
$\Delta t$ | PAO integration interval | 1.0
$N$ | Population size (all optimisers) | 100
$G$ | Iterations per run (all optimisers) | 100

(a) DeJong function.

(b) Hyper Ellipsoid function.

(c) Rotated Hyper Ellipsoid.

(d) Power sum function.

(e) Rosenbrock function.

(f) Griewangk function.

(g) Rastrigin function.

(h) Ackley function.

(i) Schwefel function.

Figure 2: Convergence histories of the best particle (averaged over 100 runs
of each optimiser) for the Molka benchmark suite in 2D.

(a) DeJong function.

(b) Hyper Ellipsoid function.

(c) Rotated Hyper Ellipsoid.

(d) Power sum function.

(e) Rosenbrock function.

(f) Griewangk function.

(g) Rastrigin function.

(h) Ackley function.

(i) Schwefel function.

Figure 3: Convergence histories of the best particle (averaged over 100 runs
of each optimiser) for the Molka benchmark suite in 8D.

##  4 Discussion

As can be seen from the figures, the PAO algorithm has performed well on the
benchmark optimisation problems. In the 2D set of problems, the performance of
PAO on the continuous, unimodal problems (from DeJong to Power sum inclusive)
is excellent, surpassed only by the SADE method which is more computationally
demanding. On the nonconvex and multimodal problems, the performance of PAO is
very much in line with the other approaches, finding exact optima (within
numerical tolerance) on both the Rastrigin and Ackley problems.

In the 8D suite of problems, the performance of PAO is excellent on the
continuous multimodal benchmarks, clearly outstripping the other methods. In
the other problems PAO performs similarly to other methods and manages to
avoid becoming trapped in poor minima on the Rastrigin and Griewangk functions
as is observed with the PSO and DE methods. In the challenging Schwefel and
Ackley functions, almost all methods (PAO included) were unable to locate the
global minima in the higher dimensional problems with default hyperparameters.

The benchmark results are especially impressive given that the attractors in
PAO implementation were simply set to the be local and global best particles
as in the standard PSO algorithm. Indeed, PAO has outperformed the PSO
approach in every problem. It might be expected that the the selection of more
specialised attractors may enable the practitioner to input domain knowledge
into the optimisation procedure resulting in greater performance. For example,
one could imagine that local minima might be avoided by using a stochastic
attractor (in which the position of the attractor is drawn from some
distribution for each particle) when the problem is known to be highly multi-
modal.

Although the benchmark performance is encouraging, the authors would make the
argument that it is not the only feature that should be considered in the
assessment of a novel heuristic. As well as optimisation power, PAO has been
shown to have a number of desirable features including interpretable
hyperparameters, exact dynamics and closed-form transition densities.

The hyperparameters of the PAO approach are familiar to engineers as the
inertial, damping and stiffness terms of a linear dynamic system. This choice
is deliberate on the part of the authors in order to permit an interpretable
way to select hyperparameters. Large choices for the inertial parameter slow
converge and encourage exploration. Whereas large choices for the values of
the stiffness parameters accelerate convergence and promote exploitation.
Particularly interesting is the choice of the damping ratio parameter $\zeta$.
As is the case for linear dynamics, the selection of values $\zeta<1$ give
rise to underdamped dynamics and oscillation, values $\zeta>1$ lead to
overdamped dynamics and slow convergence. The choice of the stochastic scaling
parameter $q$ can be related to the excitation level of the dynamics and
controls the extent to which the PAO algorithm is dominated by random search.

An advantage of PAO compared to other PSO methods is that the motions of the
particles can be expressed exactly in terms of a time-step $\Delta t$. Unlike
methods that rely on discretisation to move the particles in successive
iterations, the exact nature of the dynamics means that no error is accrued
even as the value of $\Delta t$ becomes very large. This means that the
$\Delta t$ could conceivably be altered during the run in either a scheduled
or an adaptive fashion similar to the energy function in simulated annealing
[22]. A scheduling approach could used to promote large movements during the
initial stages of the optimisation procedure (thus quickly moving to promising
regions of the search space) and then more fine-grained motions close to
convergence. In highly multi-modal problems, the opposite approach might be
taken in order to attempt to escape minima as the particles converge.

A major advantage of the proposed approach is that the transition densities
(parametrised by the integration interval $\Delta t$) are available in closed
form. One particular advantage is that at the end of the optimisation run, the
final transition density (averaged over every particle, weighted by objective
score) might be viewed as a kind of uncertainty on the position of the optima.
This could provide the practitioner with insight in the case of noisy
objective functions such as those that are typically encountered when
performing parameter identification tasks from data.

A more formal treatment of uncertainty can also be envisaged. Access to the
transitions densities would allow the PAO method to be sued as the proposal
step within an SMC scheme. This would enable a Bayesian viewpoint to be taken
and would give access to valid posterior distributions over the positions of
optima. The application of PAO within an SMC framework has the potential to
drastically improve the convergence of SMC algorithms in line with other
heuristic proposal schemes such as the DE-based approach in [23].

In this paper, a novel variant of the particle swarm algorithm has been
proposed with a number of distinct advantages including interpretable
hyperparameters, exact dynamics and closed-form access to the transition
density. The performance of the proposed approach has been shown to be in line
with other common choices for heuristic optimisation, indicating the
advantages of the proposed approach do not come at the cost of reduced
performance.

## Author contributions

The authors contributed equally to the conception and development of the
method. MDC collected the data and wrote the manuscript.

## Acknowledgements

The authors would like to thank the Isaac Newton Institute for Mathematical
Sciences, Cambridge, for support and hospitality during the programme Data
Driven Engineering (DDE) where work on this paper was undertaken. This work
was supported by EPSRC grant no EP/W002140/1.

## References

  * Koza [1994] J. Koza.  Genetic programming as a means for programming computers by natural selection.  _Statistics and Computing_ , 1994.
  * Abdel-Basset _et al_. [2018] M. Abdel-Basset, L. Abdel-Fatah, and A. Sangaiah.  _Metaheuristic Algorithms: A Comprehensive Review_.  Elsevier, 2018.
  * Hussain _et al_. [2019] K. Hussain, M. Mohd Salleh, S. Cheng, and Y. Shi.  Metaheuristic research: a comprehensive survey.  _Artificial Intelligence Review_ , 2019.
  * Joyce and Herrmann [2018] T. Joyce and J. Herrmann.  _A Review of No Free Lunch Theorems, and Their Implications for Metaheuristic Optimisation_.  Springer International Publishing, 2018.
  * Kudela [2023] J. Kudela.  A critical problem in benchmarking and analysis of evolutionary computation methods.  _Nature Machine Intelligence_ , 2023.
  * Chopin and Papaspiliopoulos [2020] N. Chopin and O. Papaspiliopoulos.  _An Introduction to Sequential Monte Carlo_.  Springer International Publishing, New York, 2020.
  * Dorigo and Gambardella [1997] M. Dorigo and L. Gambardella.  Ant colony system: a cooperative learning approach to the traveling salesman problem.  _IEEE Transactions on Evolutionary Computation_ , 1997.
  * Kennedy and Eberhart [1995] J. Kennedy and R. Eberhart.  Particle swarm optimization.  In _Encyclopedia of Machine Learning and Data Mining_ , volume 4. Springer US, 1995.
  * Storn and Price [1997] R. Storn and K. Price.  On the usage of differential evolution for function optimization.  In _Proceedings of North American Fuzzy Information Processing_ , volume 11. IEEE, 1997.
  * Del Valle _et al_. [2008] Y. Del Valle, G. Venayagamoorthy, S. Mohagheghi, J.-C. Hernandez, and R. Harley.  Particle swarm optimization: Basic concepts, variants and applications in power systems.  _IEEE Transactions on Evolutionary Computation_ , 2008.
  * Omran _et al_. [2006] M. Omran, A. Salman, and A. Engelbrecht.  Dynamic clustering using particle swarm optimization with application in image segmentation.  _Pattern Analysis and Applications_ , 2006.
  * Kennedy [2003] J. Kennedy.  Bare bones particle swarms.  In _Proceedings of the 2003 IEEE Swarm Intelligence Symposium. SIS'03 (Cat. No.03EX706)_. IEEE, 2003.
  * Sun _et al_. [2004] J. Sun, B. Feng, and W. Xu.  Particle swarm optimization with particles having quantum behavior.  In _Proceedings of the 2004 Congress on Evolutionary Computation (IEEE Cat. No.04TH8753)_ , volume 1. IEEE, 2004.
  * Sun _et al_. [2011] J. Sun, W. Fang, V. Palade, X. Wu, and W. Xu.  Quantum-behaved particle swarm optimization with gaussian distributed local attractor point.  _Applied Mathematics and Computation_ , 2011.
  * Yang _et al_. [2014] W. Yang, W. Wu, Y. Fan, and Z. Li.  Particle swarm optimization based on local attractors of ordinary differential equation system.  _Discrete Dynamics in Nature and Society_ , 2014.
  * Duan _et al_. [2022] B. Duan, C. Guo, and H. Liu.  A hybrid genetic-particle swarm optimization algorithm for multi-constraint optimization problems.  _Soft Computing_ , 2022.
  * Grassi and Pareschi [2020] S. Grassi and L. Pareschi.  From particle swarm optimization to consensus based optimization: stochastic modeling and mean-field limit.  2020\.
  * Grassi _et al_. [2021] S. Grassi, H. Huang, L. Pareschi, and J. Qiu.  Mean-field particle swarm optimization.  _arXiv preprint arXiv:2108.00393_ , 2021.
  * Thomson [1996] W. Thomson.  _Theory of vibration with applications_.  CrC Press, 1996.
  * Särkkä and Solin [2019] S. Särkkä and A. Solin.  _Applied stochastic differential equations_ , volume 10.  Cambridge University Press, 2019.
  * Qin and Suganthan [2005] A. Qin and P. Suganthan.  _Self-adaptive Differential Evolution Algorithm for Numerical Optimization_ , volume 2.  IEEE, 2005.
  * Kirkpatrick _et al_. [1983] S. Kirkpatrick, C. Gelatt, and M. Vecchi.  _Optimization by Simulated Annealing_ , volume 220.  WORLD SCIENTIFIC, 1983.
  * Turner and Sederberg [2012] B. Turner and P. Sederberg.  Approximate bayesian computation with differential evolution.  _Journal of Mathematical Psychology_ , 2012.