Self-organized Collective Motion with a Simulated Real Robot Swarm

TL;DR

This paper demonstrates a fully autonomous swarm of simple robots achieving collective motion without orientation sensors, using a modified AES method optimized by TCACS, tested in simulation environments.

Contribution

It introduces a novel autonomous collective motion approach for simple robots without orientation sensors, utilizing a modified AES method and TCACS optimization.

Findings

Successful implementation of autonomous flocking without orientation sensors

Effective movement towards a desired direction and rotation capabilities

Validated results in Matlab and Webots simulations

Abstract

Collective motion is one of the most fascinating phenomena observed in the nature. In the last decade, it aroused so much attention in physics, control and robotics fields. In particular, many studies have been done in swarm robotics related to collective motion, also called flocking. In most of these studies, robots use orientation and proximity of their neighbors to achieve collective motion. In such an approach, one of the biggest problems is to measure orientation information using on-board sensors. In most of the studies, this information is either simulated or implemented using communication. In this paper, to the best of our knowledge, we implemented a fully autonomous coordinated motion without alignment using very simple Mona robots. We used an approach based on Active Elastic Sheet (AES) method. We modified the method and added the capability to enable the swarm to move toward a desired direction and rotate about an arbitrary point. The parameters of the modified method are optimized using TCACS optimization algorithm. We tested our approach in different settings using Matlab and Webots.

Tables

Table 1. Table 1: AES Simulation Settings

Parameter	Description	Value
$v_0$	Absolute linear velocity of each agent	0.05 [m/sec]
$t_f$	Final time of simulation	30 [sec]
$N$	Number of agents	10*10
$d_{init}$	Initial distance between two side-by-side agents	0.2 [m]
${\widehat{𝒗}}_d$	Desired linear velocity unit vector	${[-1.0,0.0]}^T$ [m/sec]
$w_l$	Weight for desired linear velocity	0.8
$\omega$	Desired rate of rotation	0.7 [rad/sec]
$w_r$	Weight for desired rotational velocity	1
$D_r$	Strength coefficient of measurement noise	0.5
$D_{\theta }$	Strength coefficient of actuation noise	0.02

Equations

$${\dot{\mathbf{x}}}_i(t)=\left\{v_0(t)+\alpha \left[\left({\mathbf{f}}_i(t)+D_r{\hat{\mathbf{\xi }}}_r(t)\right)\cdot {\hat{\mathbf{n}}}_i(t)\right]\right\}{\hat{\mathbf{n}}}_i(t)$$

$${\dot{\theta }}_i(t)=\beta \left\{\left[{\mathbf{f}}_i(t)+D_r{\hat{\mathbf{\xi }}}_r(t)\right]\cdot {\hat{\mathbf{n}}}_i^{\perp }(t)\right\}+D_{\theta }{\mathbf{\xi }}_{\theta }(t)$$

$${\hat{\mathbf{n}}}_i(t)=\left[\begin{array}{c}\cos ({\theta }_i(t))\\ \sin ({\theta }_i(t))\end{array}\right]$$

$${\mathbf{f}}_i(t)=\sum _{j\in S_i}-\frac{k}{l_{ij}}\left(\Vert {\mathbf{r}}_{ij}(t)\Vert -l_{ij}\right)\frac{{\mathbf{r}}_{ij}(t)}{\Vert {\mathbf{r}}_{ij}(t)\Vert }$$

$${\mathbf{r}}_{ij}(t)={\mathbf{x}}_j(t)-{\mathbf{x}}_i(t)$$

$${\mathbf{x}}_i^{t+1}={\mathbf{x}}_i^t+\left\{v_0{\hat{\mathbf{n}}}_i^t+\alpha \left[\left({\mathbf{f}}_i^t+\frac{D_r}{\sqrt{\Delta t}}{\hat{\mathbf{\xi }}}_r^t\right)\cdot {\hat{\mathbf{n}}}_i^t\right]{\hat{\mathbf{n}}}_i^t\right\}\Delta t$$

$${\theta }_i^{t+1}={\theta }_i^t+\left\{\beta \left[\left({\mathbf{f}}_i^t+\frac{D_r}{\sqrt{\Delta t}}{\hat{\mathbf{\xi }}}_r^t\right)\cdot {\hat{\mathbf{n}}}_i^{\perp ,t}\right]+\frac{D_{\theta }}{\sqrt{\Delta t}}{\mathbf{\xi }}_{\theta }^t\right\}\Delta t$$

$${\mathbf{f}}_{aux,l}=w_l{\hat{\mathbf{v}}}_d$$

$${\mathbf{v}}_r=\omega \left[\begin{array}{c}{\mathbf{r}}_{\mathbf{ic},\mathbf{y}}\\ -{\mathbf{r}}_{\mathbf{ic},\mathbf{x}}\end{array}\right]$$

$${\mathbf{r}}_{\mathbf{ic}}={\mathbf{x}}_{\mathbf{i}}-{\mathbf{x}}_{\mathbf{c}}$$

$${\mathbf{f}}_{aux,r}=w_r{\mathbf{v}}_r$$

$$\mathbf{f}={\mathbf{f}}_{neighors}+{\mathbf{f}}_{aux,l}+{\mathbf{f}}_{aux,r}$$

$$\psi =\frac{1}{N{v}_0}\Vert \sum _{i=1}^N{\mathbf{v}}_i\Vert$$

$$\psi =\frac{1}{N}\Vert \sum _{i=1}^N{\hat{\mathbf{n}}}_i\Vert$$

$${\psi }_i=\left|cos({\gamma }_i)\right|=\left|\frac{{\mathbf{f}}_i\cdot {\hat{\mathbf{n}}}_i}{\Vert {\mathbf{f}}_i\Vert }\right|$$

$$\psi =\frac{1}{N}\Vert \sum _{i=1}^N{\psi }_i\Vert$$

$$J=\frac{1}{N_{MC}}\sum _{m=1}^{N_{MC}}\left\{\sum _{t=0}^{t_f}\left[w_1\left(\sum _{i=1}^N\Vert {\mathbf{f}}_i^t\Vert \right)+w_2{\psi }^t\right]\right\}$$

Full text

1 Introduction

Collective motion (CM) is an eye-catching demonstration of a more broad phenomenon called collective behaviour [22]. CM is observed in diverse domains such as: physics [11], chemistry [12], from micro-size creatures such as bacteria colonies [20], cells [1], macro-molecules [18] to handful macro-size examples of fish [24], birds [4, 5], even humans [15] and to non-living systems [21]. The principle feature of CM is that each individual agent behaves based on the interaction between its neighbors; resulting in polarized motion of the group. This interaction might be either simple (attraction/repulsion) or complex (combinations of simple interactions) [23].

In most of the works on CM, agents use both orientation and proximity information of their neighbors [8, 6]. In contrast to these works, some works introduced methods that do not rely on this information. An individual-based model is defined in [17], where the interaction between neighbor agents is defined only by escape-pursuit behavior. In [9], inelastic collisions between isotropic agents, and in [13] a pairwise repulsive force between deform-able, self-propelled particles are introduced as interaction mechanism regardless of the orientation of neighbors. These methods leave an indirect effect on individual agents, and as a result, enable them to establish aligning interaction in an effective, implicit manner. By so doing, not only the swarm achieves CM, but also the reduced dynamics reaches consensus in the heading direction of agents. An example of this implicit aligning is [7], in which, Ferrante et al. identified a new elasticity-based approach for achieving collective motion and illustrated the behavior by introducing an Active Elastic Sheet (AES) model. In their proposed model, the motion of an individual is determined by only attraction-repulsion forces. This feature is very useful in practical cases, such as a swarm of robots in which a specific robot have not any information of its neighbors’ orientation [3].

In this paper, we improve the AES model presented in [7], and optimize some of the parameters of the proposed algorithm using a heuristic optimization algorithm. To illustrate the improved performance of the model, various simulations with different assumptions are performed. In addition, some simulations are conducted in Webots [14] using the model of Mona robot [2]. Finally, We implement the AES model on a swarm of simulated real Mona robots for the first time. To the best of our knowledge, this is the first implementation of fully autonomous CM using such a simple robotic platform without alignment information.

This paper is organized as follows: in the Sect. 2, we introduce the concept of AES mechanism to achieve CM, and the modifications are addressed there. Then, the optimization algorithm is discussed and the objective function is defined in Sect. 3. We mention the settings of simulations in the Sect. 4. In Sect. 5, the results are shown and we discuss them. Sect. 6 closes the paper with a conclusion about the current work.

2 Collective Motion Control

In this section, the dynamic model of AES along with some modifications is presented. The model is derived based on a simple two-dimensional active elastic sheet mechanism [7]. After defining the basic mechanism, the numerical dynamics of particles and modifications will be elaborated upon.

2.1 Active Elastic Sheet Model and Numerical Dynamics

Each agent $i$ out of $N$ agents, is affected by the attraction-repulsion force of its neighbors. This force will leave an effect on both the linear velocity $\boldsymbol{\dot{x}}_{i}$ and the rotational velocity $\dot{\theta}_{i}$ of the agent. This effect is formulated in continuous-time form as:

[TABLE]

[TABLE]

where, $v_{0}$ is the forward biasing speed, and parameters $\alpha$ and $\beta$ are the inverse transitional and rotational coefficients, respectively. The unit vector $\boldsymbol{\hat{n}}_{i}(t)$ is the heading vector which is parallel to the heading direction of agent $i$, and $\boldsymbol{\hat{n}}_{i}^{\perp}(t)$ is perpendicular to it.

[TABLE]

The term $\boldsymbol{f}_{i}\centerdot\boldsymbol{\hat{n}}_{i}^{\perp}$ of Eq. (2) indicates that the angle between force and heading vectors ($\gamma_{i}$) defines the rotation speed of the agent (Fig. 1). So, the rotation speed will be zero if two vectors $\boldsymbol{f}_{i}$ and $\boldsymbol{\hat{n}}_{i}$ are parallel (i.e. $\boldsymbol{f}_{i}$ and $\boldsymbol{\hat{n}}_{i}^{\perp}$ are perpendicular.)

Agents are subject to measurement and actuation noises ${D_{r}}\boldsymbol{\hat{\xi}}_{r}(t)$ and ${D_{\theta}}\boldsymbol{\xi}_{\theta}(t)$; in which, $\boldsymbol{\hat{\xi}}_{r}(t)$ is a randomly oriented unit vector and $\boldsymbol{\xi}_{\theta}(t)$ is a random variable with standard normal distribution. ${D_{r}}$ and ${D_{\theta}}$ are noise strength coefficients. Inducing total force vector $\boldsymbol{f}_{i}(t)$ over agent $i$ at time $t$ can be calculated by the summation of linear spring-like forces, which link agent $i$ to its interacting neighbor $j$, member of set $S_{i}$:

[TABLE]

[TABLE]

${l_{ij}}$ is equilibrium distance (or alternatively called natural length) of the spring that links agent $i$ and $j$. In order to simulate the model, we need to numerically integrate the dynamic equations of motion. Integrating Eq. (1) and (2) using Euler method at time step $t$, results in the governing discrete-time kinematic equations of each particle which is expressed by the following equations:

[TABLE]

[TABLE]

here, $\boldsymbol{x}_{i}^{t}$ and ${\theta}_{i}^{t}$ are the position and orientation of agent $i$ at step time $t$, respectively, and $\Delta t$ is the numerical time-step interval of integration. It is noteworthy that the velocity vector in Eq. (1) points parallel to the heading of the agent which means it is assumed that each agent can only turn, and move forward/backward. Modeling agents with 2 degrees of freedom provides many advantages to robotic applications in which the robots are not able to move omni-directionally, e.g. Mona [2] robot.

2.2 Modifications to AES

In this section, we present some modifications to AES model, which result in collective linear and rotational motion of swarm. Firstly, we add a stimulating force to compel the swarm to move toward a specific direction. This modification leads to faster convergence of collective motion. Furthermore, the desired direction of swarm movement is also achievable by such modification. To model this auxiliary inducing force, at a certain time, we simply added the auxiliary force $\boldsymbol{f}_{aux,\ l}$ to the induced force of Eq. (4). This force is parallel to the desired velocity unit vector $\boldsymbol{\hat{v}}_{d}$ and its magnitude is defined by a weighting parameter $w_{l}$ as follows:

[TABLE]

The second modification enables the swarm to rotate collectively by adding another auxiliary force. This force is also related to another desired linear velocity which itself is proportional to the rate of rotation. According to the kinematics of 2D rotation, the linear speed of each agent is proportional to its distance with respect to the center of rotation $\boldsymbol{x}_{c}$. Thus, pointed in the following equation, the linear velocity of rotation $\boldsymbol{v}_{r}$ is related to $\boldsymbol{r_{ic}}$, the position vector of agent $i$ wrt the center of rotation, and $\omega$ is the rate of rotation. $\boldsymbol{r_{ic,x}}$ and $\boldsymbol{r_{ic,y}}$ are the x- and y-component of $\boldsymbol{r_{ic}}$.

[TABLE]

[TABLE]

Rotation center $\boldsymbol{x_{c}}$ is an arbitrary point, yet it is common to consider it on the center of swarm. Hence, same as the previous approach, we calculated the modified force for rotation by the following equation:

[TABLE]

where, $w_{r}$ is the weighting parameter of induced force for swarm rotation. It is worth mentioning that $\boldsymbol{f}_{aux,\ r}$ should be calculated for each agent, but $\boldsymbol{f}_{aux,\ l}$ is a universal force and is same for all agents. The linearity feature of the aforementioned phenomena enables us to add both terms together and form the total auxiliary force, which still satisfies both motions. So, the final modified inducing force is:

[TABLE]

$\boldsymbol{f}_{neighors}$ is exactly the same as Eq. (4), but we added the subscript to show that this is the force which is defined by the influence of interacting neighbors.

2.3 Degree of Alignment

Since the symmetry is a principal feature for achieving flocking behavior, a criteria is introduced to measure it. This metric is defined by a specific system variable, called the degree of alignment. One of the most common, and widely used definition is represented in [23]. In this case the variable $\psi$ is the average normalized velocity vectors:

[TABLE]

where $N$ is the total number of the agents, $v_{0}$ is the average absolute velocity, and $\boldsymbol{v}_{i}$ is the velocity vector of the agent $i$. However, in methods with 2-DoF agents, like AES, as previously mentioned, the velocity vector is parallel to the heading direction. So, the degree of alignment of the system is defined as [7]:

[TABLE]

Because our purpose is to simulate the AES method on robots, we added another criterion which is defined from a control viewpoint. In this respect, the alignment of each particle $\psi_{i}$ is defined by its desired direction of movement, i.e. an agent is called aligned when its orientation is parallel to its desired direction of motion. Referred to Fig. 1, the alignment is defined by the angle $\gamma_{i}$ as follows:

[TABLE]

In this definition, $\psi_{i}$ is between zero (agent’s orientation is perpendicular to the force, $\gamma_{i}=\pi/2$) and one (agent is directed toward the force, $\gamma_{i}=0$ or $\gamma_{i}=\pi$). Note that, as the norm of the unit vector $\boldsymbol{\hat{n}}_{i}$ is equal to one, it is eliminated from the denominator. Therefore, the degree of system alignment is calculated by the following equation:

[TABLE]

Mathematically, the value of this variable is between zero (asymmetric, disordered, random orientation) and one (fully symmetric orientation).

3 Optimization

In this section, the parameter optimization of the algorithm will be introduced. For so doing, we utilized a heuristic optimization algorithm called Tabu Continuous Ant Colony System (TCACS) [10]. TCACS combines two algorithms: Continuous Ant Colony System (CACS) [16], and Tabu Search (TS) [19]. We will mention a brief description of the algorithm in the next subsection. Our aim is to find the optimal decision parameters for AES model with a specific setup configuration. The decision parameters to be tuned are $\alpha$, $\beta$, and $k$.

3.1 TCACS Optimization Algorithm

The structure of the algorithm is very similar to CACS; however, the advantages of tabu concept are also added. Borrowing the concept of tabu balls (meta-spherical shapes in the decision-space), prevents the ants choosing their destination within the tabu regions. The flowchart of TCACS is illustrated in Fig. 2.

In this paper, we considered the parameters of TCACS as $N_{ants}=15$, $\gamma=0.5$, $m=2$ (which are number of ants, weighting factor, and PCA factor, respectively) and used Roulette as the weighting strategy.

3.2 Objective Function

One of the most important steps in optimization is the definition of objective function, which can leave a considerable effect on the performance of the algorithm to solve the problem. Various objective functions can stand as candidates, depending on what the aim of the problem is and which conditions needs to be satisfied. In case of collective motion using the modified AES method, we have two main objectives:

• •

To minimize the total induced force, which, generally speaking, means that agents are at their ideal position according to the elastic force of the sheet. Satisfying this objective will result in maintaining the shape of swarm same as its natural shape.

• •

To maximize the rate of convergence, which will be reached when the orientation of agents are about the same as the desired ones. To fulfill this objective, the optimization algorithm will try to maximize the degree of system alignment.

Although the defined problem is multi-objective, we combine both previous functions, multiplied by their weights, into a single objective function. This will help us to simplify and solve the parameter tuning problem in a more straight-forward manner.

It is notable that the system is subject to noise; consequently, the objective function gets stochastic, and it is necessary to reduce the effect of noise on it, otherwise, the performance of optimization will be reduced. Therefore, we calculate the average of multiple Monte Carlo simulations as the output of the objective function. With this in mind, the objective function is as follows:

[TABLE]

In which, the first term of the RHS is calculated by Eq. (12), and the second term by Eq. (16). Additionally, $N_{MC}$ is the number of Monte Carlo simulations which is set to 10, and $t_{f}$ is the final time of simulation. Considering the fact that the first objective is a minimization problem, and the second objective is a maximization problem, we determined $w_{1}$ (the weighting parameter of the first function) positive and $w_{2}$ negative. As the result, the optimization can be alternatively called minimization, and the objective function is a cost function. Determining the value of weighting parameters depends mainly on the purpose of the optimization, however, a common initial setting is normalizing the order of each term; so, we considered them as $w_{1}={1}/{N}$, and $w_{2}=-1$.

4 Experimental Setup

To model Mona robot [2] in Webots, we assembled different parts including upper and lower boards, motors, and wheels which is shown in Fig. 3-b. A population of 100 Monas, where presented in the simulations and each one has its own specific ID in order for the supervisor to recognize them in the simulation world-model, and calculate the adequate displacement of linear and rotational motion. Besides the simulations in Webots, we also conducted some simulations in Matlab to check the modifications and to run optimization programs. In this paper, we represented three different setups for simulations in both Webots and Matlab, which will be described in the following subsections.

4.1 Setup 1: Linear Motion

The first simulation is executed to prove the linear motion of swarm is achievable thanks to the first modification. The configuration of this setup is similar to Table 1, but the rate of rotation is set to zero, $\omega=0$. The result of this Matlab simulation is illustrated in the first row of Fig. 5.

4.2 Setup 2: Rotational Motion

Similarly, the second modification is verified by another simulation, in which the setting makes the swarm rotate about its center. The setting of Table 1 is analogous to the current setting, yet the desired linear speed is set to zero, $\lVert{\boldsymbol{\hat{v}}_{d}=0}\rVert$. The second row of Fig. 5 represents a time lapse of this simulation in Matlab.

4.3 Setup 3: Combination of Linear & Rotational Motion

This simulation is conducted in Webots. We considered the robots as independent particles, each of which is determined by its orientation and position in the xy-plane. In this paper, a supervisor code is programmed to calculate the position and orientation of robots and dictate them to the corresponding agent. In other words, it is assumed that the controller of robots are ideal and the desired states of robots which supervisor calculates, is the same as the actual state of robots.

The settings of this simulation is placed at Table 1. By defining such settings, we expect the swarm to rotate about its center while moving collectively toward the west direction. To depict the result of the simulation, we saved a log file from Webots, and then plot the position and rotation of each robots in Matlab, results of which are placed in the third row of Fig. 5.

5 Results & Discussion

Applying TCACS optimization algorithm on the previously defined problem provided us the tuned parameter of the algorithm for a specific setting of simulation. The setting was considered as Table 1. In the optimization, the maximum number of function evaluations was limited to 300. Finally, the optimized parameters achieved as $\alpha^{*}=0.066$, $\beta^{*}=0.97$, and $k^{*}=1.28$. The shape of swarm was considered as a square with 10 rows and 10 columns of agents.

In order to show the effect of optimized parameters on the performance of flocking, we investigated those objectives with both the optimized and non-optimized set of parameters. Fig. 4-b shows the significant improvements, in comparison with the default parameters Fig. 4-a, which is achieved by optimized parameters. As it is shown, the rate of convergence is increased, so the CM achieved faster. Besides, the total force of the network is decreased, meaning that the shape of the swarm is maintained better.

We also studied the effect of objective weights $w_{1}$, and $w_{2}$ on the optimized parameters. Since the behavior of swarm varies for each obtained parameters, we defined two more different weight sets to show the effect of different weights on swarm behavior:

• •

More emphasis on the first objective, $w_{1}=10/N$, $w_{2}=-1$ (Fig. 4-c)

• •

More emphasis on the second objective, $w_{1}=1/N$, $w_{2}=-10$ (Fig. 4-d)

As we expected, the former setting resulted in the more reduction of total force, and the latter one, increased the rate of convergence of alignment, yet made the system more sensitive to noise (the effect of which is shown on the fluctuating behavior of force on Fig. 4-d). The first row of Fig. 5 proves that the initially perturbed formation of the swarm is eliminated, then the collective motion of swarm is achieved. The inducing linear force caused the swarm to move toward the west direction, as it was expected by $\boldsymbol{\hat{v}}_{d}$. Due to the robustness of AES to noises, noises were not effective enough to prevent the swarm from flocking. The time history of the second row, for rotational motion, proves that the second modification is verified, similarly. In this simulation, the swarm rotates about its center and maintain its shape, meanwhile.

The third row shows that the combination of linear and rotational motions made the swarm to move toward the desired direction, and in the mean time, rotate about its center. The same optimized parameters were used for this simulation, in Webots, proving that similar settings results in similar behavior in Matlab and Webots.

6 Conclusion

We have addressed the parameter optimization and Webots simulations of AES model, which is an elasticity-based approach for achieving CM. Some of its advantages such as not relying on exchanging orientation information, and assuming 2-DoF agents attracted our attention. So, considering our robotic case, we added some modifications helping the swarm of robots to move toward the desired direction and turn about an arbitrary rotation center point. In addition, a new criterion for measuring the alignment of the swarm is defined based on the desired direction of motion.

utilizing TCACS optimization algorithm, we tuned the parameter of AES model for a specific setup containing combination of linear and rotational motion of swarm. The effect of optimization on the behavior of CM, such as increasing the rate of convergence and decreasing the total force, has been shown. Various weighting strategies for optimization has been studied, too.

Finally, three different scenarios have been conducted in Matlab and Webots. For this end, we have designed a CAD model of Mona robot in Webots and have simulated a supervised program to control agents. The results proved that our modifications are applicable and the optimized parameters improved the performance of the AES mechanism.