Self-organized pattern formation in a swarm system as a transient phenomenon of non-linear dynamics

Abstract

This article presents a microscopic model (agent positions, directions and interactions are explicitly modelled) of mobile agents (or self-propelled particles) that is inspired by the ‘complex transport networks’ reported by Jones (2010; The emergence and dynamical evolution of complex transport networks from simple low-level behaviours, International Journal of Unconventional Computing 6, pp. 125–144). Here, the agents' positions are modelled continuously. This multi-agent system (or artificial swarm) shows a wide variety of self-organized pattern formations. The self-organization is based on the non-linearity of the agents' turns (discrete jumps in the agents' directions) and the indirect interactions of the agents via a potential field that determines their motion (high values are attractive) and which is changed by themselves (agents increase the value of the potential field at their positions). At least most of the irregular and complex patterns are transient. The patterns found during the transient are more complex than those the system converges to. Still, this transient behaviour is relevant. We empirically investigate the transient times in dependence of several system parameters and give examples.

Keywords:

• swarm system
• transient
• self-organization
• non-linear dynamics

1. Introduction

Natural and artificial systems consisting of many interacting parts, which are extended only to local interactions or local information, often show global, non-trivial structures and patterns. Such cases were, for example, reported for stone and soil formations [1], agitated rod-like particles [2], flow of granular substances [3,4], bacteria [5], ants [6] and cellular automata [7]. In this work, we focus on an abstract model of self-propelled [8,9], mobile particles or agents [10].

Our model is microscopic, that is, agent positions, directions and interactions are explicitly modelled. Its complement is a macroscopic model that does not describe individual properties of particles. An example from physics would be the Navier–Stokes equations representing a macro-model and the corresponding micro-model would be the classical mechanics approach describing every particle using Newton's laws.

Microscopic modelling of self-organizing multi-agent systems [11] and swarms goes at least back to Reynolds [12]. The field of physics called self-driven particle systems, which was started by Vicsek [13], is related (see also [14]). Generally, the field of multi-particle, multi-agent and swarm models is vast [15–21]. Concerning microscopic swarm models, focus is often on structure formation [6,21–23], influence of noise [13,24] or the measurement of physical quantities [25,26].

The patterns shown by our model are similar to those found in stone and soil formations [1]. Similar patterns are also found in density-dependent models of chemosensitive movement of bacteria [5]. These similarities indicate that our model shows natural and possibly generalizable behaviours despite its abstractness. Furthermore, it is a microscopic model that shows behaviours similar to those of macroscopic models [5], which, in turn, indicates the potential existence of a micro–macro link [15]. The investigated system is inspired by and similar to the ‘complex transport networks’ reported by Jones [27,28]. In the system reported here, however, the agents' positions are modelled continuously. In Jones' work originating in image processing, the agents are said to cover small squares (pixels) that can only hold one agent at a time. In this work, the agents have real-valued positions . As we are simulating this model on a computer, we end up with a discrete representation making our approach similar to that of Jones. However, the resolution of positions is much higher here and the observed patterns are different.

The investigated system shows a wide variety of patterns. However, it turns out that most of the irregular patterns are transient. The patterns found in the transient are more complex than those the system converges to. Many dynamical systems show transient behaviour before they converge to a fixed-point, periodic or chaotic behaviour (e.g. differential equations [29], coupled map lattices [30], cellular automata [7] or simple mappings [31]). Traditionally, transient phenomena are not investigated in fields such as mathematics (e.g. non-linear dynamics) or physics because they are a temporary phenomenon before a steady state is reached. Therefore, they were not considered to be interesting for research. However, it turned out that many phenomena are transient but should still be investigated as explained by Prigogine [32]:

After I had presented my own lecture on irreversible thermodynamics, the greatest expert in the field of thermodynamics made the following comment: ‘I am astonished that this young man is so interested in nonequilibrium physics. Irreversible processes are transient. Why not wait and study equilibrium as everyone else does?’ I was so amazed at this response that I did not have the presence of mind to answer: ‘But we are all transient. Is it not natural to be interested in our common human condition?’

We believe the transient behaviour observed in the model presented in this article is relevant because many complex self-organized phenomena seem to be based on transient behaviour [32–38]. Long transients were observed in high-dimensional dynamical systems, for example, coupled map lattices [37], they were observed in ecological models [33,34,36] and models of chemotaxis [38]. Transients become, for example, relevant, ‘if the transition time is comparable to the lifetime of the system’ [38]. This might be true for several multi-agent and swarm systems.

As we are more interested in complex patterns rather than regular or trivial patterns, the question arises whether there exist parameter settings such that permanent complex patterns develop. Or will all systems always degenerate to trivial patterns after a transient? If self-organized complex patterns were always transient we would like to control the transient length. This would be a critical issue for all evolved and engineered swarm-like systems.

2. The model of the multi-agent system

In this section, the investigated multi-agent or swarm system is explained and defined. The system consists of several autonomous, mobile, simple reflex agents [10]. Such agents are autonomous in their decisions (here, they decide on the direction of their motion), however, these decisions are purely reactive and rely only on the agents' current perception. The sensory perception is implemented by sensors that measure local values of a potential field. The potential field is an abstraction of the environment. In natural systems, this could, for example, be a temperature gradient (cf. [39,40]). The agents' behaviour is defined by condition–action rules (i.e. certain conditions trigger certain actions). These agents are modelled as self-propelled particles [8,9]. They have an unlimited energy reservoir that allows them to move autonomously. The system consists of a high number of such agents and is called multi-agent system [11,41]. In reference to similar natural systems, they are sometimes also called swarm system [42,43].

The agents move in two-dimensional space. The change of a single agent's position x over time is defined by

(1)

for constant velocity v. The change of the agent's direction φ over time is defined by

(2)

for a function that defines the direction of the agents' turns (counterclockwise, no turn or clockwise). α implements the agents' autonomous decisions based on the current sensor input and will be defined below. Δ defines the absolute value of the turn angles (see below). The sensor readings of the central sensor s _c, of the left sensor s _l, and of the right sensor s _r are defined by local values of the potential field as follows:

(3)

(4)

(5)

for a potential field P representing the environment, , sensor angle σ and sensor distance d. Figure 1 shows the sensor setting.

Figure 1. Sensor setting (cf. Equations (3)(3) (4) –(5)(5)).

The function α, that implements the condition–action rules of the agents, is defined closely following [27]:

(6)

Note that order matters here (i.e. Equation (6)(6) should be read as an algorithm): The first condition that is met starting from the top defines the agents' behaviour. The treatment of in the last condition is only for mathematical completeness because it is extremely improbable to occur in the simulation (for reasonable definitions of the potential field). The random turn has a probability of 50% for +1 and 50% for −1. The last two conditions define a gradient ascent. The first condition assures an exploitation of good directions and the second condition assures an exploration behaviour for situations when the agent is in a valley of the potential field.

Δ defines the absolute value of the turn angles:

(7)

for a constant rotation angle and a time interval τ at which the agents' turn and their directions are updated (in this work we set τ = 1). Thus, we obtain a synchronized system that is discrete in time (all agents decide on their turning behaviour at the same time). Note that the rotation angle φ_rot can effect qualitative differences in the agents' behaviour. Small rotation angles correspond to small angular speed of the turns. Related to the speed of the agents and the size of the world, the agents tend to move in lines for small rotation angles. For big rotation angles, they tend to cluster because they virtually turn on the spot. The system could, for example, be extended by defining Δ as a stochastic process. This would transform Equation (2)(2) into a stochastic differential equation.

The temporal evolution of the potential field P is basically defined by a diffusion process combined with additions at agent positions and subtractions (e.g. cooling, evaporation). The addition at agent positions introduces a self-enhancing process because peaks in the potential field attract other agents and aggregated agents create peaks in the potential field (see also [40]). In principle, it is defined by the following partial differential equation

(8)

for diffusion D, evaporation rate η, addition θ (the Dirac delta indicates that an agent only contributes to the potential field at its position), number of agents N and agents' positions x _i (t). However, for simplicity and to reduce the computational complexity, the diffusion and the evaporation (first two terms in Equation (8)(8)) are only executed at unlike the addition process (third term in Equation (8)(8)) that is executed at . Furthermore, in our implementation the diffusion is not normalized in correspondence to the resolution of the used grid that is used as spatial discretization. Note that we use standard numerical methods to compute the diffusion on a standard (discrete) computer. The grid is always quadratic and the number of grid points is given by L ². Periodic boundary conditions apply (i.e. space is a torus).

The observed behaviour is based on the non-linearity of the agents' turns (jumps in the direction) and the combined effects of averaging, blurring and evaporating in the potential field. The system is defined to be deterministic in most cases except for symmetry-breaking situations (random turns defined by Equation (6)(6)). Thus, the system's dynamics is mostly determined by the inhomogeneities in the initialization.

3. Pattern formation

Table 1 gives the typical parameters that were used, if not explicitly stated otherwise. In Figure 2, a typical example of a complex, most likely transient pattern is shown (N = 6750, L = 600, t = 3000). Note that the agents generate a potential field showing parallel nearby lines. This seems to be a qualitative difference to the patterns reported by Jones [27] and might, therefore, be induced by the higher resolution of agent positions (continuous positions approximated by floating point arithmetic).

Table 1. Typical parameters

Sensor angle σ	45°
Rotation angle φrot	45°
Rotation period τ	1 (time units)
Grid length L	400
Sensor distance d	0.0375 (units)
Velocity v	0.01 (units/τ)
Diffusion D	0.1 ((1/L) /(10τ))
Evaporation η	0.04 (1/10τ)
Addition θ	5 (1/τ)
Number of agents N	1000
Active cell threshold δactive	N/100
Simulated steps	1 × 10^5 (time units)

Figure 2. Typical example of a complex, most likely transient pattern for N = 6750, L = 600 and t = 3000.

Figure 3 shows typical examples of the patterns in the potential field obtained by different rotation angles φ_rot. These plots demonstrate the variety of patterns that the investigated systems produce.

Figure 3. Patterns for varied rotation angle φ_rot: (a) 15°, (b) 35°, (c) 45°, (d) 55°, (e) 80°, (f) 90°, (g) 105°, (h) 165°, (i) 255°, (j) 275°, (k) 285°, (l) 300°, (m) 330°, (n) 350° and (o) 355°

4. Investigation of the transients

In the following we want to investigate the transient behaviour of this system. We need a measure of complexity to analyse the dynamics of the transient. A subjective measurement of the pattern complexity is the simple estimation by looking at the plots of the potential field. However, we need an objective and automatable metric to measure the pattern complexity. We use a simple metric that just computes the percentage a of grid points of the potential field that are above a threshold δ_active (called active cells in the following). We obtained good results by setting δ_active = N/100. This simple metric is a valid measure of complexity for the investigated parameter sets because complex patterns are composed of many cells with values above the threshold, while simple patterns have fewer active cells. In general, a high active cell percentage does not assure a complex pattern. However, in the parameter region, that is investigated in the following, this is assured by simulation. Moreover, small changes in the temporal evolution of the active cell percentage correlate with a static pattern in the potential P. Thus, the temporal evolution of the active cell percentage indicates quasi-stationary and stationary states.

A typical example for L = 150 is shown in Figure 4. The percentage of active cells decreases over time until a steady state is reached. Figure 5 shows four plots of the potential field for the 4 times marked in Figure 4 (t ₁ = 450, t ₂ = 5 × 10³, t ₃ = 1 × 10⁴ and t ₄ = 2.5 × 10⁴).

Figure 4. Temporal evolution of the active cell percentage a for L = 150; crosses indicate times for which the potential field is given in Figure 5.

Figure 5. Four plots of the potential field of the example shown in Figure 4. (a) t ₁ = 450, a ₁ = 64.0%; (b) t ₂ = 5 × 10³, a ₂ = 31.2%; (c) t ₃ = 1 × 10⁴, a ₃ = 27.8%; and (d) t ₄ = 2.5 × 10⁴, a ₄ = 14.1%.

It is well known that the transient times scale exponentially with the dimensionality in many high-dimensional systems (e.g. see [44]). In the following, we want to estimate a lower bound of the average transient time depending on the system's size which is the grid length L here. We do this by measuring the time until a pattern without bifurcations is reached (cf. Figure 5(d)) for several samples. We believe that such bifurcation-free patterns are a steady state for systems of grid size L ≤ 200, which is supported by our experience with long runs of the simulation. However, it is, in principle, possible that fluctuations lead to the formation of a new (and possibly stable) bifurcation, albeit we have never observed that. In order to get an objective measure of when exactly a bifurcation-free pattern is reached and the simulation can be stopped, we use the active cell percentage a again. For each grid length L, a threshold is determined that characterizes the upper bound of the active cell percentage for bifurcation-free patterns. The bifurcation-free pattern of highest active cell percentage is a long line that twines diagonally around the torus. Also several unconnected lines occur (usually not more than two for the investigated grid lengths) that also have a high active cell percentage. For the five investigated grid lengths (L ∈ {100, 125, 150, 175, 200}) we got thresholds that decrease approximately proportional to increased grid lengths (a _thresh ∈ {31%, 25%, 20%, 18.5%, 17%}). Furthermore, we used an additional parameter to increase the detection rate of bifurcation-free patterns. We determine the order in the agents' directions by determining the length of the sum of the direction vectors of all agents:

(9)

If this value got bigger than 0.4 (i.e. a high percentage of the swarm moves in similar directions), we had the simulation stopped as well. For grid lengths of L ≥ 150, we found patterns with bifurcations for which we could not determine within t ≤ 1 × 10⁶ whether they are unstable and degenerate a bifurcation-free pattern. These samples were thrown out and were not included in the calculation of the lower bound of the average transient time. The result is shown in Figure 6 and it suggests an exponential increase of the actual average transient time.

Figure 6. Increase of a lower bound of the average transient time with increasing grid length L (not all outliers shown), n = 50 for each grid size.

In contrast to Figure 4, the behaviour of the system is more complex for bigger grid sizes and bigger agent numbers. In Figure 7 we give two examples of the active cell percentage for a bigger system. In case of Figure 7(a), the behaviour is dynamic most of the time. New patterns emerge and break down again. Even values of a < 2% are reached corresponding to patterns of low complexity. However, these states are not stable and more complex patterns emerge again. Also time intervals of little dynamics are found (e.g. t ≈ 2 × 10⁵, t ≈ 3 × 10⁶, t ≈ 5.2 × 10⁶). Whether this irregular behaviour is transient stays unanswered. In contrast, Figure 7(b) gives an example of a system that most likely converges to a stable pattern after a long transient of about 1.4 × 10⁶ time steps.

Figure 7. (a) Example with active cell percentage staying active; (b) example with active cell percentage converging.

In the following, we want to test the influence of other system parameters on the transient length. The influence of the diffusion process cannot be fully excluded from the above experiments when increasing the grid resolution. When varying the diffusion coefficient only, experiments showed no significant dependence concerning the transient length (data not shown). Therefore, we completely removed the diffusion process from the system by setting D = 0 and focused on the evaporation rate η. It turned out that there is a dependence of the transient length on the evaporation rate. High evaporation results in short transient lengths as expected because low-populated parts of patterns have low values in the potential field and are likely to evaporate to values close to 0. Low evaporation results in long transients. The transients increase exponentially with decreasing evaporation rate as shown in Figure 8. A macroscopic interpretation of the evaporation process is that it lowers the exploration rate. Fluctuations in the potential field (e.g. an agent that left a highly populated line of agents) will less frequently be reinforced. Thus, new bifurcations in the patterns are less likely to emerge and have a higher probability to vanish.

Figure 8. Influence of the evaporation rate η on the transient length t; n = 10 samples each.

5. Discussion

The vast variety of patterns found in this system by changing only one parameter is higher than expected. Even small changes of the rotation angle lead to qualitative changes in the observed patterns in certain parameter intervals. This is also the reason that we cannot tell whether all main classes of patterns are shown in this article. A parameter interval, within which a potentially yet unreported pattern might be found, could be small.

The exponential dependence of the transient lengths on the system size L is known and was reported, for example, in coupled map lattices before [37]. The exponential dependence of the transient lengths on the evaporation rate is a less typical phenomenon in dynamical systems. The system dynamics can be interpreted metaphorically as a cooling process, that is, the system has initially a high temperature that is cooled by evaporation. High temperatures correspond to complex patterns in the potential field and low temperatures correspond to simple, bifurcation-free patterns. The natural cooling process by evaporation is characterized by an exponential decay. Thus, the difference in the time when a certain value is reached exponentially scales with the evaporation rate.

Another way of explaining the role of evaporation in this scenario is to compare it with the averaging of velocity vectors as reported in [13]. The patterns in the potential field result in an implicit averaging process of velocity vectors of at least groups of agents (i.e. convergence to non-trivial distributions of velocities). High evaporation rates result in faster decays of the patterns and ultimately in a faster convergence of the velocity distributions. Lower evaporation rates slow this convergence down. At the current state of research, it stays unanswered whether there exist parameter settings that lead to permanent complex patterns. Due to the discrete simulation, all patterns will at least become cyclic eventually, because the number of possible states is finite. However, acyclic patterns might exist in principle.

Given the transient lengths can be arbitrarily increased it is possible that multi-agent systems as the one reported in this article show interesting behaviour only in their transients.

6. Summary and outlook

In this article, a self-organized system was reported that generates a variety of complex patterns. At least for small system sizes, these patterns are transient. A lower bound of the transient time was shown that suggests an exponential increase in the average transient time with increasing system size.

The presented microscopic model seems to be a good object of research to investigate the relevance and influence of transient behaviour in self-organized swarms. An investigation of the influence of noise to this system is pending. In first experiments, no qualitative changes for reasonable intensities of noise were found concerning the formation of patterns and the influence to transient lengths. This will be the focus of future work.