Formation Control Algorithm of Agents Based on Earth Mover’s Distance

Massive sport, such as unmanned aerial vehicle performance, often needs fast and efficient calculation of formation morphing and individual path planning.*is paper introduces a novel fast formation control method of a crowd. First, we get the agents’ location in a 2D polygon with centroidal Voronoi tessellation and L-BFGS techniques. *en, we transform crowd formation shapes with a global shortest motion path pair assignment using earth mover’s distance algorithm. Finally, the repulsing force between agents and obstacles is calculated based on the recursive velocity observer method control agents’ motion. Extensive experimental results show the effectiveness and usefulness of our algorithm in 2D group formation transformation.


Introduction
Crowd animation simulation, which is a fundamental problem in computer graphics, augmented, virtual, and mixed reality, has received much attention in recent years [1,2]. With the rapid development of computer hardware and software technology, especially the continuous improvement of GPU rendering technology, simulating animation of a massive number of virtual agents on a limited computing capability hardware is a problem that has recently become quite popular and crowd formation control has become one of the most important issues in the research field of Multirobot System (MRS). It plays an important role in a wide range of applications, including industry, military affairs, aviation, military drills, urban planning, art and sports, and video games. Especially, the manufacturing industry has an increasing demand for automatic processing of large objects or dangerous goods. In the military, the power and command efficiency can be improved by maintaining the military equipment or the army as reasonable and appropriate formation by the group simulation system. Agents usually interact with the user, and it not only improves operation efficiency but also ensures human safety. With the help of virtual agents, we can effectively complete the task and ensure the safety of humans in dangerous activities, such as exploring unknown environments, demining, patrolling, and reconnaissance. However, the crowd formation transform is still quite limited in terms of both agent number and individual intelligence. Designing and rehearsal of group formation are tedious and time-consuming; the crowd formation transformation simulation system is an effective solution for this problem. rough team formation control, formation design, and crowd simulation, the simulation system can combine with art and sports perfectly.
ere are three major classes of algorithms for crowd formation control, namely, formation generation, formation maintenance, and formation transformation. e main difficulties are collision avoidance and path planning. If the motion matching between individuals is not accurate, the collision avoidance between individuals is difficult to be realized. It is complex and unrealistic to use artificial methods to realize the matching. Leader-follower dynamic technique [3] uses a linearization method to design a velocity controller for transforming, but it depends on leader agents. When the leader follower fails, the formation will be out of control. Zheng et al. [4] proposed a group formation control method for formation based on geometric constraints, which can be applied to perform agent pairing, path planning, and collision avoidance for the crowd; the calculating process still needs a tedious iterative step and is timeconsuming in real-time simulation.
To form a uniform group formation distribution, efficient calculation of individual motion matching, and smooth formation transformation, we propose a novel group formation transform approach to automatically generate a homogeneously distributed crowd and visually pleasing formation transformation. e group is uniformly sampled by the centroid Voronoi structure centroidal Voronoi tessellation (CVT). Based on the earth mover's distance (EMD) algorithm, we precompute a global shortest motion path assignment in the offline process. During the real-time transformation simulation, our contributions are summarized as follows: (1) We introduce a novel group formation morphing approach that can automatically generate a global shortest motion path pair assignment based on geometric constraints (2) We introduce a new method to enable a swarm of insects to follow the motion of a moving shape by applying EMD method

Group Formation Morphing.
Group formation control is a vital problem of many crowds and multiagent systems. Based on the Velocity Obstacle concept [5], Van den Berg et al. [6] proposed the concept of Reciprocal Velocity Obstacles for real-time multiagent navigation among autonomous decision-making entities, and it is widely used in the crowd simulation system. Anderson et al. [7] introduced a flocks animation simulation method with shape constraints to improve the motion concerning the underlying behavioral model. Ji et al. [8] proposed a method to simulate the movement modeling, path planning, and collision avoidance, which provided an auxiliary tool for the group calisthenics to improve the design of the formation and pattern. Zheng et al. [9] adopt geometric constraint idea for smooth formation animation of regulated crowds. Chen et al. [10] presented a robust multiagent model for large-scale controllable shape constrained simulation of flying insects. Wang et al. [11] proposed a linear formation control method under directed and switching topologies for distributed robot systems. Zheng et al. [4] proposed a geometric constraint group formation control method for group calisthenics, which is combined morphing, Centroidal Voronoi Tessellation, and Lloyd techniques and it is an efficient solution for individual distribution layout, individual matching, path planning, and collision avoidance. Furthermore, Zheng et al. [12] introduced a regulated heterogeneous flock formation control method based on a capacity-constrained centroidal power diagram to model the internal layouts of groups in the constraint shapes. Some researchers convert the formation problem into tracking problem, use the pseudorigid body method to calculate the path of the expected point of the agent, and use the reverse step method to design the formation control of the nonlinear control algorithm [13]. Based on backstepping control and bioheuristic neural network, Ding et al. [14] used a trajectory tracking method for threedimensional formation control and obstacle avoidance. Henry et al. [15,16] proposed a crowd control algorithm using deformed grids for single traversal search, and it can be used to control crowd interactive with the environment and collision avoidance. Xu et al. [1] proposed a natural group transformation and control method between source shape and target shape based on Kuhn-Munkres and mutual information method, which can be used for the interaction and collision avoidance of between individuals and environments.

Voronoi Diagram.
In the centroidal Voronoi diagram, the distribution of seeds is very uniform. Due to these properties, centroidal Voronoi graphs can be used in a wide range of applications, including image compression, numerical integration, finite difference, location problem, and meshing [17]. e Voronoi diagram is suitable for not only uniform sampling but also sampling with constraints; it can be used to calculate a distribution as uniform as possible and compute region partition with capacity constraint. Based on the quasi-Newton method, Liu et al. [18] proposed a superlinear convergence algorithm for centroidal Voronoi diagram and successfully applied it to the important problem of remeshing. Yan et al. [19] proposed a method to compute the centroidal Voronoi diagram limited to the curved surface, which improved the efficiency of remeshing. After that, the centroidal Voronoi diagram has been widely studied in digital geometry. Lévy and Liu [20] introduced a centroidal Voronoi diagram method based on p-norm and demonstrated its application in quadrilateral meshing. Pan et al. [21] used the centroidal Voronoi diagram to construct grid surfaces with constant mean curvature. Wang et al. [22] proposed a fast algorithm for calculating the Riemannian center and the center of mass for any geodesic Voronoi diagram. Recently, Liu et al. [23] further improved the calculating speed of the geodesic centroid Voronoi diagram based on the evolutionary algorithm. To compute the Voronoi diagram with capacity constraint, Sorkine et al. [24] and Zheng et al. [25,26] proposed some new optimization strategies, respectively, and their algorithms have obvious improvement in computing efficiency. e distribution problem with a capacity constraint can be interpreted by the optimal transmission theory, Desbrun et al. [27] provided a first-order convergent optimization algorithm and found that the blue noise generated by the capacity constraint has better spectral characteristics. Xin et al. [28] named this problem as the centroid power diagram problem and proposed an optimization algorithm with superlinear convergence. e group distribution sampling phase computes the sampling position based on the centroid Voronoi diagram. Pair assignment precalculation phase formulates the globally shortest path matching as a linear assignment optimization problem. e movement control phase combines target-oriented motion and interaction force between nearest agents to ensure collision-free movements for all the individual agents in real time.

Sampling Number Calculation.
For the source and target formations, we first convert the formation shapes to Delaunay triangulation (DT) representations [29] in the interior of the shape. en, we calculate the area s k of each triangulated subshape k for a formation with multiple shapes. Finally, we use the followed formula to calculate the number n k of sampled agents in each triangulated subshape: where N is the total number of sampled agents and m is the number of subshapes.

Group Distribution Sampling.
e Voronoi diagram is a partitioning of a given domain Ω ⊂ E d into regions based on distance to a set of points p i ⊂ E d n i�1 in a specific subset of the domain. at set of seeds is specified beforehand, and the region containing point p i is called the Voronoi cells; for each seed, there is a Voronoi cell consisting of all points closer to that seed than to any other. P � p 1 , . . . , p n represents the discrete point set on the domain Ω; the Voronoi cell of the site p i can be denoted as follows: Let the domain Ω be endowed with a density function ρ(x) > 0 which is bounded; then, the centroid c i of each Voronoi cell Ω i is given by that is, the points c i that serve as generators for the Voronoi regions Ω i are themselves the mass centroids of those regions. We call such a region a centroidal Voronoi diagram. en, the CVT energy function E of the Voronoi tessellation Ω i is defined as [17] In this case, minimizing the energy function E, we can get all the centroids c i of the tessellation which directly ensures that each position p i is a representative of c i . Liu et al. [18] and Xin et al. [28] use gradient descent method to calculate the center of CVT rapidly; the gradient of the energy with regard to x i is defined as where c i is the centroid and Ω i ρ(x)dx is the mass of Voronoi regions Ω i .
It is worth noting that all point sampling is located inside of shape using the above CVT method. However, groups who are strictly distributed according to user-designed shapes need to evenly distribute a part of the sampling on the shape boundary to achieve a more orderly crowd distribution in crowd formation transform simulation. erefore, the sampling points of the group formation need to be calculated in two parts. In the first part, you need to sample m points on the boundary; namely, it needs to generate m points on the boundary zΩ of the shape and the Voronoi region of zΩ denoted as V i . e remaining k of the sampling points need to be distributed in the internal region zΩ of shape and W j represents the Voronoi region of zΩ, as Figure 2 shows. e whole energy function of CVT is defined as whose gradient with regard to x i is where K � N − sign(⌊s/l⌋) × ( �� N √ − 1), m � N − k, s is the surface's area of shape, and l is the boundary's circumference of shape. We are able to adopt the L-BFGS solver [30] to minimize the objective function, leading to empirically superlinear convergence rate. e pseudocode of uniform sampling in shape with CVT is shown in Algorithm 1

Position Assignment.
e goal of group formation transformation is to change all the agent from source group formation to target group formation. e formation transformation not only realizes collision avoidance during crowd movement but also achieves the global shortest moving distance. In group formation transformation simulation, agents need to observe the movement and position of other agents, maintain the coordination of the team, and avoid the free and scattered individual behavior. erefore, it is necessary to realize the optimal movement control of the group and keep the number of the crowds unchanged. We apply the classical earth mover's distance (EMD) algorithm [31] to solve the global position assignment problem. e group formation transformation objective function is defined as where s i is the i-th position of from source group, t j is the j-th position of target group, δ i is the distance from s i to its nearest neighboring point, and δ j is the distance from t j to its nearest neighboring point. e limiting condition of objective optimization function is subject to where n i�1 X i,j � 1 denotes that each agent in the source group can only move to one target position and m j�1 X i,j � 1 means that a position of the target group can only be occupied by one agent, as shown in Figure 3.

Movement Control.
In movement control phase, we can also mix our method with other methods for agent-based planning [6,32]. Recursive Velocity Observer (RVO) method is used to calculate the actual speed of the agent. In addition to collision avoidance, the agent may change the target exit according to the updated environmental information or change the path through the area and then need to replan the path. For each time step, we first advance the individually controlled agents and then treat them as obstacles in the continuum model so that the crowd motion will flow around them.

Results and Performance
We implemented and experimented with our algorithm on a computer with a 64-bit version of Win7 system, a 3.07 GHz Intel(R) Core(TM) i7 CPU, and 6-GB memory. e coding language is C++. We shall use extensive experimental results to exhibit and demonstrate the algorithm's efficiency, as well as its insensitiveness to scene resolution. To demonstrate the effectiveness and practicability of our method, we carry out the evacuation simulation in different scenarios.
In our method, the most computationally time-consuming step is the position assignment, which has the time complexity of O(N 3 ); we precalculate the result of the position assignment in the offline process. We recorded the time cost for the different simulation scenarios (Table 1); the numbers of agents are 200, 300, 400, 600, and 1000, respectively. As shown in Table 1, the cost time difference is very small on a different shape which contains the same number of agents. For example, the solving time of the objective optimization function is about 3 seconds for pair assignment of a group with 200 agents, while the cost time is 5 seconds and 45 seconds for groups which contain 300 and 1000 agents, respectively. It is suitable for offline formation preordering operations. Input: collective crowd shape; an error tolerance ε � 10 − 6 . Output: group sampling points c i n i . Generate the initial points p i n i as the seed; Compute the Voronoi region Ω i for each seed p i Compute the centroid c i of each Voronoi region Ω i ; While ∇E > ε do / * L-BFGS * / Compute the energy function E(x i , Ω i ) in equation (6) Compute the gradient w.r.t. x; see equation (7) Update x by step search end return the final centroid points c i n i ALGORITHM 1: CVT-based method for group sampling. 4 Mathematical Problems in Engineering 4.1. Performance. Figure 4 illustrates the process of motion transformation of different groups in a single shape. Figure 4(a) shows that the "A" word shape formation transformed into "H" word formation and then into "H" word formation with 200 agents. Figure 4(b) shows a 300-person square formation transformed into a five-star formation and then into a circular formation. Figures 4(c) and 4(d) illustrate the group with 400person formation from pig shape to rabbit shape and rabbit shape to skiing logo shape, respectively. Figure 5 shows the process of motion transformation of crowds in multiple shapes. Figure 5(a) shows that the group has 200-agent formation from "2018" word with 4 subshapes to the Olympic rings logo which contains 5 subshapes. Figure 5(b) shows the transform process from the shapes of English word "welcome" to the shapes of Chinese word "welcome." From Figure 5, we can observe that our method has good adaptability in nonconvex and irregular formations with multiple shapes, and the transformation process is smooth and without any confusion. Figure 6 shows the process of formation transformation with 400 agents in diamond shape under nonuniform density. We set the shape with an e density function which is defined as [4] ρ(x, y) � | (x, y)(x + y)|. en, the diamond can also be decomposed into four symmetrical shapes. It can be seen that our algorithm can get different formation transformation effects.

Evaluations and Comparisons.
To quantitatively analyze the effect of a crowd formation transformation, we introduce two different objective measures, as described as follows.

e Rangeability of Sampling Points.
For agent s i , the neighboring points set of s i within its the nuclear radius r i is denoted as s i np � s p | ‖s i − s p ‖ < r i . Let min dis(s i ) be its the minimum distance from position s i to its neighbors s i np , and t i is the current position of s i in the target group G(T, P T ). e source position s p ∈ s i np will move to its target position t p . e standard deviation σ s and the average value μ s can be used to quantify the stability rangeability of source group sampling during the transformation:   Clearly, the lower value of σ s indicates that the agents have more similar distances from their neighbors.

Transform Distance from Source to Target Group.
We also employ the standard deviation σ dt and the average value μ dt of transform distance to quantify the motion effect of agents. For an agent p i , its transform distance is defined as follows: Here, s i and t i denote its source position and its current position, respectively. e average value μ dt of all the agents' transform distance is defined as e standard deviation σ dt of transform distance is calculated as follows: Clearly, the lower value of σ dt and μ dt indicates that the agents have less time consumption during the crowd transformation process.
We compare our method with the other three crowd formation transform simulation algorithms including RVO [6], geometry-constrained mechanism method [4], and mutual information-based method [1]. As shown in Figure 7, the four different methods are applied to the same source and target formations. e visual results show that our approach can achieve the best visual effect among the four approaches. Table 2 shows the comparison of results from our method on crowd simulation and other algorithms using simulation time, μ s , σ s , μ dt , and σ dt . From this table, we can observe that our simulation method on crowd formation transform outperforms other approaches no matter what kind of evaluation metric is used, except the Zheng method.  Figure 4: Group formation between single formations. Group formation from (a) "A" word to "H" word and to "Z" word, (b) rectangle to five-star and to circular, (c) pig to rabbit, and (d) rabbit to skiing logo.  (b) group formation from English word "welcome" to Chinese word "welcome."

Mathematical Problems in Engineering 7
When comparing with the Zheng method, our algorithm takes less simulation time.

Conclusions
In this paper, we introduce a novel crowd formation transform method based on Earth mover's distance. e main precompute step of our approach is to sample shapes with CVT technique and use EMD to calculate the global shortest position assignment. In the real-time simulation, the interactions between individuals and obstacles are calculated based on RVO which realize the collision avoidance effectively. Large experimental results show that our algorithm has an obvious speed advantage in the realtime simulation process, and our method can also be used to realize the formation transformation in three dimensions.
For large-scale group, formation transformation simulation, the precomputation process, and the collision avoidance calculation process are still time-consuming. In the future, we will keep improving the performance by advanced parallelization techniques. We will also explore more applications in visual computing.  Data Availability e simulation data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.