Abstract
The kinematic behavior of mobile robots can be represented as functions of time. During the operation of a multirobot system, the orbit of a special robot is recorded. The embedding dimension and the delay time are chosen based on the correlation integral method. A chaotic attractor equivalent to the original system is reconstructed in phase space. The multirobot system can be adequately described based on the phase space information, and the dynamic system states can be forecast based on this information. The eigenvalues of the attractor are calculated including the maximum Lyapunov exponent and correlation dimension. The robot collective behavior is described and analyzed quantitatively based on the eigenvalues. The critical factor that affects the interaction of robots is investigated based on quantified parameters. Our analysis results can be used to improve the understanding of robot interaction mechanisms.
1. Introduction
With the development of robotics, multirobot systems have received considerable attention [1]. A single robot has many limitations in the ability of the acquisition, processing, and control for the information [2]. The research of multirobot system began in the late 1970s. Multiagent system theory provides a new way in thought and application [3]. The distributed artificial intelligence, complex system theory and biochemical approaches are recently introduced in robotics [4, 5]. The scientists began studying the organization, interaction, and learning mechanism for multirobot system [6–8]. The current solutions of multirobot have made impressive progress. But the theory framework and implementation need further improvement.
A multirobot system is a set comprised of robots that interrelate and act reciprocally [9]. The system properties can be studied using pertinent state variables. The current state of the system is entirely determined by the previous states at a given time [10]. Research in kinetic systems investigates the regularity with respect to time of system state variables [11]. The state variables are an object set that can describe the system’s dynamic behavior completely. In a multirobot system, any physical object that can determine the system’s properties can serve as the status information such as speed, acceleration, gesture, position, behavior, and the mutual relationship among robots [12, 13].
Multirobot learning is the process of acquiring new cooperative behaviors for a particular task [14]. Typically, a multirobot learning method can be classified as simulation, an experiment, or mathematical modeling [15]. In a multirobot system, a large number of homogeneous mobile robots interact with each other based on a shared environment. Beckers et al. conducted initiative simulations and physical experiments to interpret the nesting behavior of termites with a stigmergy mechanism [16]. Seo et al. discussed the use of evolutionary psychology in order to select a set of personality traits that could evolve due to a learning process based on reinforcement learning [17]. Kobayashi et al. proposed an objective-based reinforcement learning system for multiple autonomous mobile robots in order to acquire cooperative behavior. The proposed system used profit sharing (PS) as a learning method [18]. Lee et al. presented an algorithm for behavior learning and online distributed evolution for cooperative behavior of a group of autonomous robots. Individual robots improved the state-action mapping through online evolution with the crossover operator based on -values and their update frequencies [19]. Generally, the aforementioned methods can be used to prove that a robot system is meeting specific tasks. However, there is no guarantee that the control strategies can be implemented effectively on different robot platforms or under different environmental conditions.
Many studies of robot cooperative behavior using mathematical methods can be found to determine the major factor for the system behaviors and have enhanced the understanding of emergent properties of robot collective behavior. Probabilistic robotics is concerned with perception and control in the face of uncertainty [20]. Fox et al. presented a statistical algorithm for collaborative mobile robot localization. The approach is capable of localizing mobile robots at any time by using a sample-based version of Markov localization [21]. Martinoli et al. proposed a time-discrete macroscopic model to capture the dynamics of a robotic swarm system engaged in a collaborative manipulation task [22]. The critical parameters that affect the robots’ collective behavior performance can be analyzed by using a mathematical model of multirobot collective behavior. The optimal parameters of the system can be chosen via mathematical analysis. This approach can provide the essential theory basis for the design and analysis of multirobot collective behavior.
The collective behavior of mobile robots emerges from the interaction among the individuals and also between the individuals and the environment. The evolution of such behavior is a highly complex dynamic process. The movement form of behavior is often chaotic [23, 24], so most existing modeling and design methods for robot behavior are insufficient to describe the complexity of robot collective behavior on a mechanism. The methods used to examine mobile robotic collective behavior include mathematical modeling and quantitative analysis of mobile robot behavior [25]. And this is an urgent problem of theory and technology in the practical applications of robot behavior learning.
We analyzed the interaction between robots and the environment based on a proposed dynamic system theory developed in this study. A system law in multidimensional phase space is investigated based on the evolution track of a special robot. Data points over time are first obtained for the special robot. The appropriate embedding dimension and delay time are chosen. The phase space equivalent to the original system is reconstructed. The multirobot system can be adequately described based on the information of the phase space. The dynamic system states can be forecast based on this information. Then, the property of the attractor in the phase space is analyzed. The eigenvalues of the attractor, including the Lyapunov exponent and correlation dimension, are calculated. A quantitative description and an analysis of multirobot collective behavior are described, based on the eigenvalues. The key that affects the interaction of robots is finally investigated based on the quantified parameters. Our analysis results help us to better understand robot interaction mechanisms.
2. Phase Space Reconstruction
In the application process, the univariate time series for a given time period can generally be measured from a real physical system [26, 27]. However, the kinetic equations of the system cannot be determined. Therefore, the phase space reconstruction is the first stage in the analysis and modeling of a chaotic time series observed from nonlinear systems. Typically, all possible state variables cannot be obtained through the observation of a real process. Due to the way in which the components of a physical system are coupled to each other, the phase space reconstructed from a single observation by time delay embedding can actually reflect the law of a chaotic system. Packard et al. suggested using the delay coordinate of a variable reconstructed in phase space for the original system [28]. The appropriate embedding dimension delay time can be estimated according to Taken’s embedding theorem. The reconstructed attractor ensures that the topological properties are essentially in agreement with the original dynamic system [29].
When reconstructing the attractor from the scalar time series, the most critical issue is the selection of the delay time and embedding dimension [20]. However, the finite real data and noise cause difficulty in estimating the parameters for phase space reconstruction. Normally, the selection of parameters for phase space reconstruction can be classified in two categories: one in which the parameters are estimated independently and one in which they are mutually dependent. The conventional C-C method can be practicable for the selection of a delay time and the delay time window. Then, the embedding dimension can be estimated based on the delay time and delay time window.
The time series , where , can be reconstructed in a phase space. The system model can be described as follows: where ; ; and and are the embedding dimension and delay time, respectively.
The correlation integral of the embedded time series is defined as follows:where , , for , and for .
The time series , where , is subdivided into t disjoint time series in order to investigate the dependence of the nonlinear system. The subdivided time series are represented as follows:where is the length of subseries, , where INT denotes reserving integer of the value, and is the delay time.
For a general , the statistic for every subsequence is defined as follows:where is the correlation dimension; is the radius of the region in multidimensional space.
As goes to infinity, the statistic can be represented as follows:
If the time series is an independent and identical distribution, is equal to zero constantly for a fixed value , and and . However, since real data sets are finite and correlated between the components of a series, is not equal to zero in common practice. Thus, the locally optimal times can be either the zero crossings of or the times at which shows the least variation with , since this indicates a nearly uniform distribution of points. The dispersion of the statistic corresponding to the maximum radius and minimum radius is defined as follows:where and are the maximum and the minimum radii, respectively.
According to statistics results, values of range from 2 to 5, values of range from 1 to 4, and . The averages of the statistics can be defined aswhere is the average of the statistics, is the dispersion of the statistics, and is the correlation statistics.
The delay time can be accurately estimated based on the first zero crossing of or the first local minimum of . However, the optimal embedding window will be relatively inaccurate corresponding to the global minimum of . Although chaotic systems oscillate without periodicity, low-dimensional chaotic systems show pseudoperiodicity. The mean orbital period could naturally be associated with the mean time between two consecutive visits to a Poincare section. For a time series with a mean orbital period , all points are in the same Poincare section in phase space for fixed values and , , and , where is an integer greater than zero. The local maximum of correlation integral is obtained since the distance is most proximate between points that are in the same section. Then, we define the mean correlation integral as follows:where is the mean orbital period of the time series, , , and k is an integer greater than zero.
The statistics of the mean correlation integral are defined as follows based on the disjoint time series given by (3):where and are the variables for statistics.
The value of has many fluctuations for the chaotic system, which has stronger inner randomness. The value of is consistent with that of apart from the periodic points in simulations. But the curve of has completely opposite trends compared to in the vicinity of periodic points. To filter the fluctuations in the neighborhood of the periodic points, is defined as follows:where is the statistics of the mean correlation integral and is the average of the statistics.
The curve of has local peaks without the fluctuations in the neighborhood of the periodic points. The optimal embedding window will thus be estimated accurately. The embedding dimension can be obtained according to the relation between the embedding window and delay time. The combined relation can be represented as follows:where is the embedding window and is delay time.
The phase space can be reconstructed as follows:where is a scalar time series, ; m is the embedding dimension; and is the delay time. The reconstructed system preserves the topological properties of the original dynamic system.
The periodic points of the mean correlation integral have significance when choosing the optimal embedding window. The quantity of has local peaks without the fluctuations, apart from the periodic points. The estimation of the optimal embedding window will be more accurate compared to the classical C-C method. The improved C-C method makes use of the inherent property of chaotic systems. The optimal embedding window will be estimated first when reconstructing the phase space, and subjectivity will be avoided when choosing embedding dimension and delay time.
3. Eigenvalues of Attractor for Multirobot Collective Behavior
The phase space reconstruction theory, first proposed by Lv et al. [24], is the origin of research on chaotic systems from a time series point of view. Any component of a system is determined by other components that interact with each other. The relevant information of the system is implicit in the evolution of any component. Therefore, the original law of the chaotic system can be returned based on a time series of a variable. The movement orbits are the important parameters of a multirobot system; these determine robot behavior. The time series of a robot position on the - or -axis can be used to investigate the robot’s behavior. The attractor of the original system can be reconstructed from the input and output time series in a multidimensional phase space. Robot behavior can be described quantitatively based on the eigenvalues of the attractor.
The interaction among robots can be analyzed based on dynamics system theory. The kinetic equations of a multirobot system are normally unknown. In the present study, the original law of the system was investigated based on the movement orbit of a special robot. During the operation of a multirobot system, the orbit of a special robot was recorded as shown in Figure 1. A rich supply of data, which included 1000 data points, was collected. The embedding dimension and the delay time were chosen based on the aforementioned method, wherein m = 10 and delay time . The phase space can be reconstructed based on the embedding dimension and delay time. The attractor of the robot movement orbit is shown as Figure 2. As a result of the corrupting noise, the real time series will contain some randomness. The evolution of robot behavior is chaotic, which means that every point in the reconstructed space is approached arbitrarily but is not coincident on the pseudoperiodic orbits, and the orbit is restricted in a special area of the phase space.
The first stage of our study of nonlinear systems is to detect and quantify chaos. The trademark of a chaotic system is its sensitivity to initial conditions. The Lyapunov exponents are the important indices with which to quantify the sensitivity of a dynamical system. If the nonlinear system is chaotic, the maximum Lyapunov exponent must be positive. Therefore, the maximum Lyapunov exponent must be calculated in order to judge a dynamic system. The maximum Lyapunov exponent is a quantity that characterizes the rate of separation of infinitesimally close trajectories in phase space. There are many algorithms available that can be used to estimate the maximum Lyapunov exponent of a dynamical system. In this study, the maximum Lyapunov exponent was calculated from small data sets using the following procedure:(1)A fast Fourier transform is applied to time series , and the delay time and mean period are calculated.(2)By calculating the correlation dimension , the embedding dimension m can be determined according to the formula .(3)The phase space is reconstructed as based on the delay time and embedding dimension .(4) is calculated for every point in phase space, which is the distance between this point and the position after the th discrete time step:(5)The average of can be calculated as follows:where is the number of and is not equal to zero. The model parameters are estimated using a least squares method. The slope of the straight line is the maximum Lyapunov exponent. The maximum Lyapunov exponent for the orbit of multirobot collective behavior is shown in Figure 3. The estimate of the largest Lyapunov exponent is about 0.12 bits/s.
The existence of the attractors in multidimensional phase space is the primary characteristic of a chaotic dynamical system. The dimension of the attractors reflects the complexity of the topology structure and the kinetic information of the attractors. In the evolution of a chaotic dynamical system, any state cannot return to the previous states accurately. The behavior of a chaotic system is nonperiodic. The correlation dimension reflects the proximity of a system state to the corresponding previous states at a pseudoperiod point.
The geometric texture of the chaotic system attractors is formed by the orbit cuddling and separating over and over again. The attractors have a self-similar structure with an infinite number of levels. The dimension of the attractors reflects the geometrical properties of a nonlinear dynamic system, and is associated with the pseudoperiodicity of chaotic systems. The stronger the nonperiodicity, the greater are the dimensions of the system. The G-P algorithm can applied to calculate the dimension of the attractors based on the observed data. An dimensional space is reconstructed by phase space reconstruction technology. The number of all associated vectors of a pair is calculated in the phase space. The distance of the associated vectors should be no greater than the given radius . The distance is calculated by adopting the -norm. The distance of the associated vectors is the maximum distance between the pair vectors. The correlation integration is defined as the ratio of the number of all associated vectors to all pair vectors in phase space, as follows: where is a unit function of Heaviside, which is defined as follows:
Supposing that the center of sphere is any point in the attractor, the radius is equal to . Then, the correlation dimension is the ratio of the number of points in the sphere to all pair vectors in the attractor. The correlation dimension depends on the value of the radius . If the radius exceeds the limit, then all the distances among the pair points are less than the radius, and the correlation dimension is equal to 1. If the radius is too small, then all of the distances among the pair points are greater than the radius, and the correlation dimension is equal to 0. Therefore, in practice, a suitable radius should be chosen with specific conditions. In the chaotic attractor, the correlation dimension increases with the value of radius . The correlation dimension of the attractor is equal to the slope of the curve of . The slope can be estimated as follows: where is roughly equivalent to , is the standard deviation of the data series, and .
The correlation dimension depends on the embedded dimension and the distance of the associated vectors. It is an uncertainty to calculate the embedded dimension and the associated distance simultaneously. The embedded dimension increases for some given radius in numerical work. A logarithmic relationship is established based on every embedded dimension, which is ~. A straight line is fitted by the method of least squares. The slope of the fitted straight line is the associated exponent. The associated exponent increases with the embedded dimension up to a certain threshold value. The threshold value is the associated dimension of the chaotic attractor. The embedded dimension is the minimum embedded dimension which ensures that the topological structures of chaotic attractors can be unfolded completely. During the operation of the multirobot system, 1000 data points on the x-axis of the position of a special robot were recorded. Figure 4 shows the associated dimension when the range of the embedded dimension is from 2 to 10. The slope of the red straight lines is the associated exponent. The associated exponent increases with the embedded dimension until reaching a maximum value. The maximum associated exponent is the associated dimension. The associated dimension tends to be stable when the range of the embedded dimension is from 2 to 10, as shown as Figure 5. The estimate of the associated dimension of robot cooperative hunting behavior is about 5.5. When topological structures of chaotic attractors can be unfolded completely, the associated dimension is the minimum embedded dimension. Therefore, the associated dimension can measure the complexity of the chaotic attractor in phase space.
4. Relationship between the Eigenvalues and Multirobot Collective Behavior
The attractor is a set of numerical values, which projects onto dynamic system. The dynamic system has a multidimensional space represented by the evolving variables [30]. And it is difficult to understand everything about the system. The attractor is a subspace which has main characters of the original data.
In a multirobot system, every individual is capable of exerting influences upon others. Therefore, the behavior of any individual comes as the result of the mutual influences of the individuals. The attractor is reconstructed based on the orbit of a special robot to reduce the dimensionality of the complex collective behavior. And the attractor is the projection of the physical system with the principal component.
The robot behavior is chaotic since the maximum Lyapunov exponent is positive. The motion of the robot is governed by a definable law. Therefore, the behavior of the robot is predictable in some time domain. The greater the maximum value of the Lyapunov exponent, the less time the robot behavior can be predicted in. The maximum forecast period is approximately equal to the inverse of the maximum Lyapunov exponent. The maximum forecast period of the robot behavior is 8.3 s.
The associated dimension is the maximum associated exponential as the embedded dimension increases. When topological structures of chaotic attractors can be unfolded completely, the associated dimension is the minimum embedded dimension. Therefore, the associated dimension can measure the complexity of the chaotic attractor in phase space. On the other hand, the minimal and maximum independent variables to describe the system can be determined based on the associated dimension. The number of the independent variables is estimated as follows:where INT represents rounding to integers; is the associated dimension. Therefore, the minimal and maximum independent variables to describe the system should be 7 and 13, respectively.
5. Conclusion
The kinematic behavior of mobile robots can be represented as functions of time. Therefore, a multirobot system is a nonlinear dynamic system. The system law in multidimensional phase space was investigated based on the evolution track of a special robot. During the operation of the multirobot system, 1000 data points on the -axis of the position of the special robot were recorded. The appropriate embedding dimension and delay time were chosen based on the correlation integral method. The chaotic attractor equivalent to the original system was reconstructed in phase space. The multirobot system can be adequately described based on the information of the phase space. The dynamic system states can be forecast based on this information. Then, the property of the attractor in the phase space was analyzed. The eigenvalues of the attractor were calculated including the maximum Lyapunov exponent and the correlation dimension. The robot collective behavior can be described and analyzed quantitatively based on the eigenvalues. The key that influences the interaction of robots was finally investigated based on the quantified parameters. This analysis results from this study will aid in a better understanding of robot interaction mechanisms.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (61573213, 61233014, 61473179, and 61473174); the Natural Science Foundation of Shandong Province (ZR2015PF009 and ZR2014FM007); the Shandong Province Science and Technology Development Program (2014GGX103038); the Special Technological Program of Transformation of Initiatively Innovative Achievements in Shandong Province (2014ZZCX04302); the China Postdoctoral Science Foundation (2014M551907); and the Independent Innovation Foundation of Shandong University (2013ZRQP002).