A High-Order Kalman Filter Method for Fusion Estimation of Motion Trajectories of Multi-Robot Formation

Multi-robot motion and observation generally have nonlinear characteristics; in response to the problem that the existing extended Kalman filter (EKF) algorithm used in robot position estimation only considers first-order expansion and ignores the higher-order information, this paper proposes a multi-robot formation trajectory based on the high-order Kalman filter method. The joint estimation method uses Taylor expansion of the state equation and observation equation and introduces remainder variables on this basis, which effectively improves the estimation accuracy. In addition, the truncation error and rounding error of the filtering algorithm before and after the introduction of remainder variables, respectively, are compared. Our analysis shows that the rounding error is much smaller than the truncation error, and the nonlinear estimation performance is greatly improved.


Introduction
As a typical complex system, swarm dynamics systems have been widely studied since first being was proposed. A swarm is composed of many simple intelligent individuals. The interaction between individuals based on basic mechanical rules can stimulate highly coordinated swarm behavior [1][2][3][4]. The dynamic behavior of swarming in nature shows amazing charm. The system composed of interconnected and constantly moving individuals emerges as colorful and highly coordinated swarming behavior, which provides a rich source of ideas for the cognition and optimization of industrial and social groups. It has considerable prospects in the fields of multi-agent collaboration [5], UAV formation [6], multi-robot collaboration [7], and intelligent grids [8].
In a multi-robot formation system, because the robot can detect the external environment it is a control system with pattern recognition, complex task execution and allocation, autonomous behavior decision-making, and other functions necessary to operate in a hostile environment and complete a variety of complex tasks. Therefore, a multi-robot formation system exhibits the characteristics of high consistency and collaboration with the outside world. In this relatively complex control system, it is necessary to obtain its precise state quantity first in order for control to be applied and to achieve the purpose of completing complex tasks. Therefore, it is necessary to perform more accurate state estimation for multi-robot motion.
Among the state estimation methods, in 1960 R. E. Kalman et al. [9,10] proposed a recursive state estimation method suitable for a system in which the model is linear and the noise obeys a Gaussian distribution. Kalman filtering (KF)is a time domain filtering method; when the system process noise and measurement noise are mutually independent Gaussian white noise with zero mean, a Kalman filter is the optimal unbiased estimator in the sense of minimum variance [11].
In general, the motion equation of a multi-robot formation and the measurement equation of its sensor are both nonlinear equations. The motion equation has weak nonlinear characteristics because only simple operations such as turning and going straight are available. The sensor observes the robot target, which shows strong nonlinear characteristics. For nonlinear state estimation problems, the traditional filtering methods are mainly extended Kalmanl filter (EKF), unscented Kalman filter (UKF) and cubature Kalman filter (CKF). EKF uses a Taylor series first-order approximation expansion of the nonlinear equations, which introduces linearization errors and has low estimation accuracy for estimating systems with strong nonlinearity [12]. Over time, Julier et al. [13] proposed another nonlinear filtering method, the unscented Kalman filter. The core of the UKF algorithm is the UT transformation, which uses an ensemble of sigma sampling points to approximate the nonlinear system model [14][15][16]. The UKF algorithm is more computationally intensive than the EKF algorithm, and its filtering accuracy is better, at least to the second-order accuracy of the Taylor expansion. However, the shortcomings of the UKF algorithm are obvious in that it can only be applied to Gaussian systems, and a problem with the nonpositive definite error covariance matrix is caused by the emergence of negative weight in the implementation of the algorithm [17]. The frequency of this problem becomes greater as the order of the system increases and the nonlinearity grows stronger. In order to better solve the problem of UKF being non-positive definite when dealing with high-dimensional covariance matrices, resulting in divergent filtering results, Ienkaram proposed a weightselective UKF filtering method called Cubature Kalman Filter (CKF) [18,19]. CKF optimizes the sigma point sampling method and weight distribution in UKF using a spherical integral and radial integral. Compared with traditional EKF, it has improved estimation accuracy and improved stability of filtering methods [20].
In recent years, in the treatment of nonlinear estimation problems, the structure of the extended Kalman filter (EKF) is used to design recursive estimators as a new idea to solve the problem.The main idea is to transform the linear error of the nonlinear equation after doing the first-order Taylor expansion into the form of a product of the scaling matrix associated with the problem solved and some time-varying positive definite matrix, then set the upper bound of the covariance matrix [21][22][23], and get the upper bound of the minimum covariance matrix by numerically solving the Ricatti equation, and finally use recursive filtering, which cannot be accurately resolved to get this scaling matrix, and then the process of solving the inverse matrix in the covariance matrix obtained by the Ricatti equation is more complicated, which greatly increases the computational burden in some high-dimensional occasions and cannot meet the real-time requirements [24]. The method cannot accurately analyze the scaling matrix, and the process of solving the inverse matrix in the covariance matrix obtained by the Ricatti equation is relatively complicated. The other is to borrow the Taylor expansion idea, expand the higher-order term with the original equation, design the new state and measurement equation, and use recursive filtering to get the required state estimate [25,26]. It can choose the expansion order reasonably according to the nonlinear strength and weakness of the actual problem, which is a good balance between the requirements of arithmetic power and estimation accuracy, and has a good real-time filtering feature [27,28].At the same time, this method has also been extended to correlation entropy filtering [29] , Kronecker product [30] and lithium battery life prediction [31] , which all show good estimation performance.
In traditional state estimation for multi-robot operation, the nonlinear state estimation is performed separately for each robot using methods such as EKF, UKF, or CKF. In addition to the problems already mentioned in this paper, performing individual state estimation for each robot reduces the estimation accuracy [32][33][34]. In recent years, the use of trigger probability in the multi-robot formation problem as the connection relationship between robots has gradually become a mainstream idea [24]. The reason for real-time estimation is that at a certain moment the connection relationship between multiple robots is deterministic, and there is only a connection or no connection, which is equivalent to a binary connection.
On this basis, this paper addresses the multi-robot formation state estimation problem by proposing a centerless multi-robot state joint estimation method based on high-order extended Kalman filter (HEKF). This improves the original multi-robot connection relationship based on trigger probability and improves the multi-robot state estimation accuracy. The main algorithm of this paper is as follows. In view of the nonlinear characteristics of multi-robot motion and observation, Taylor's higher-order expansion of the robot's state equation and measurement equation is performed to solve the problem of different nonlinear degrees of the two nonlinear equations. The high-order Kalman filtering algorithm is improved by expanding the truncation error with the state for dimensional estimation and fusing the state information of neighboring node robots to improve the multi-robot state estimation accuracy. A numerical simulation experiment proves the effectiveness of the algorithm. Therefore, we choose the second method based on high-order Kalman filtering. The main work of this paper includes: The new high-order linear equation is used to obtain the state estimation value using the recursive filtering algorithm, and the effectiveness of the algorithm is analyzed.
The rest of this paper is organized as follows. Section 2 proposes state equations and observation equations for the multi-robot formation operating system and designs a recursive estimator based on the EKF structure, while Section 3 presents the recursive estimator form derived from the second section, introduces the remainder variables into the state equation and the observation equation, and presents a new high-order linear filter derived through dynamic modeling, then analyzes its performance. In Section 4, an indoor multi-robot simulation numerical experiment is used to verify the effectiveness of the algorithm proposed in this paper. Finally, conclusions are provided in Section 5. In addition, the important mathematical symbols in Appendix A.

Problem Description
The movement of the robot is generally carried out by operations such as going straight and turning, and the observation system is generally composed of radar or visual sensors. Therefore, both the robot's motion and the observation system have nonlinear characteristics. In a multi-robot formation environment, there is a certain degree of motion consistency between the robots, and thus a certain coupling relationship. The dynamic equations can be described by the following N node coupling equations: where x i (k) ∈ R n and y i (k + 1) ∈ R q denote the state vector and the measurement vector, respectively, i and k denote the node index and the time instant, respectively, f (·) and h(·) are known nonlinear functions that are assumed to be continuously differentiable, c > 0 is the overall coupling strength, and Γ is the inner-coupling matrix. The process noise, w i (k),and the measurement noise, v i (k + 1), are assumed to be mutually uncorrelated zero-mean white Gaussian noise with the covariances Q i (k) and R i (k + 1), respectively.
The EKF structure is used to design the recursive estimation of the multi-robot formation trajectories (1) and (2): wherex i (k + 1|k) andx i (k + 1|k + 1) denote the predicted and updated estimates at time instant k + 1, respectively., and K i,k+1 is the gain matrix to be determined As in the EKF, we define the updated estimation error and the corresponding covariance as follows:

Algorithm Description
In the second section, we use the EKF structure to design the recursive filter. This idea generally expands the nonlinear function to obtain the first-order terms and truncation errors. The traditional EKF algorithm directly discards the truncation errors and keeps only the first-order terms. In the transfer function with strong non-linearity a large amount of information will be lost, which makes the estimation result inaccurate. In recent years, in the study of nonlinear filtering problems, the following two ideas have mainly been adopted in the design of recursive filters using the EKF structure. The first is to transform the truncation error into the form of a product of the scaling matrix and a time-varying positive definite matrix related to the problem, then set the upper bound of the covariance matrix, obtain the upper bound of the minimum covariance matrix by numerically solving the Ricatti equation, and finally use recursive filtering, which cannot accurately analyze the scaling matrix, before obtaining the covariance matrix through the Ricatti equation. The inverse matrix process is more complicated, and in occasions with high dimensions the computational burden is greatly increased to the point that real-time requirements cannot be met. The second is to use the idea of Taylor expansion to expand the high-order term while using the original equation to expand the dimension, design new state and measurement equations, use recursive filtering to obtain the required state estimation value, and avoid first above method. The "dimension disaster" caused by inversion in the dimensional situation, as well as the expansion order, can be reasonably selected according to the nonlinear strong and weak characteristics in the actual problem, which takes into account the requirements of computing power and estimation accuracy and has a good real-time filtering feature. A block diagram of the algorithm is shown in Figure 1.

Taylor Expansion of Nonlinear System and Introduction of Remainder Variables
Among the nonlinear filtering algorithms, the extended Kalman filter (EKF) algorithm is widely used. Its main idea is perform a first-order Taylor expansion of the nonlinear state equation around the filter estimated value. The main process is as follows.
Γx j (k); then, x i (k + 1) is subjected to a first-order Taylor expansion at the filter estimatex i (k|k) to obtain the following expression for x i (k + 1): where the state transition matrix is In the same way, the first-order Taylor expansion of the nonlinear observation function h i (x i (k + 1)) at the state prediction valuex i (k + 1|k) can be obtained: From this, we can see that in the original extended Kalman filtering algorithm, which uses only the first-order term information of the nonlinear function, the high-order term information and the truncation error are directly discarded, which makes the estimation accuracy low in occasions with a high degree of nonlinearity and leads to inaccurate estimation results. In certain extreme cases, the filtering results may even directly diverge. Here, we introduce remainder variables into the state equation and observation equation, respectively, to compensate for the loss of high-order information in the original filtering algorithm.

Dynamic Modeling of Remainder Variables and Establishment of Higher-Order Linear Systems
The remainder variable is obtained by integrating the constant term with the truncation error after Taylor expansion of the nonlinear state function and the observation function around the filter estimation value and the one-step prediction value, respectively. By using the remainder variable information to improve the estimation accuracy, the original state equation and observation equation need to be expanded and rewritten to establish a high-order Kalman filter. Lemma 1. The n-order Taylor formula of the binary function z = f (x, y) at (x 0 , y 0 ): Assuming that z = f (x, y) is continuous in a neighborhood of (x 0 , y 0 ) and has continuous partial derivatives up to order; if (x 0 + h, y 0 + h) is any point in this neighborhood, there are Lemma 2. Taylor expansion of a multivariate function at x k : 3.2.1. Pseudo-Linear Representation of State Function f i (x i,k ) According to the above lemma, the state transition function f i (x i,k ) in Formula (1) can be expressed as a polynomial form using Taylor expansion, as follows.
r) be a set of l-order hidden variables corresponding to the original system variable x(k). Theorem 2. The weight vector corresponding to the l-order latent variable vector is then We use the following two states as an example: . . .
tion is in a pseudo-linear form with higher-order latent variables x (l) (k)(l = 1, 2, · · · , r) and remainder variables β(x(k)). There is only a change in representation; as there is no essential difference, this is called pseudo-linearization.

Pseudolinear Representation of the Measurement Function
Similarly, the measurement function in formula (2) is expressed as a polynomial form by Taylor expansion, as follows: where γ(x(k + 1)) denotes the remainder variables of the measurement function: . . . . . .

Dynamic Modeling and Higher-Order Linearization of Remainder Variables Based on Taylor Expansion
In order to convert the pseudo-linearized model established above into a real highorder linearized form, it is necessary to expand the dimension and model using the highorder latent variables and remainder variables as new variables of the system.
In order to solve the problem in Remark 1, the high-order latent variable x (l) (k) and the remainder variables β(x(k)) and γ(x(k + 1)) are used as time-varying parameters and the linear coupling relationship between the l-order latent variable x (l) (k + 1) and the u-order latent variable x (u) (k) is established, along with the dynamic model relationship of the remainder variables β(x(k)) and γ(x(k + 1)): β(x(k + 1)) = A(β)β(x(k)) Then, the latent variables and remainder variables are used as the extension of the original state vector to realize the linear representation of the state model. Here, A (u) l (k) and β(x(k + 1)) can be identified according to the input information of the original state model; when there is no prior information, the settings are as follows: Then, The argmented state equation can be described as Based on the above-mentioned argmented state variables, the linearized description model of the observation equation is similarly described as follows: . . .
x (l) (k + 1) . . . Denote Then, the linearized representation of the nonlinear measurement model equation is as follows:

Projection Matrix Analysis
Rewrite the expanded state variables in Equation (46) as the original state variables and all higher-order variables and remainder variables in the following form: where K x (k + 1) and K β (k + 1) are the gain matrices corresponding to the original variables and the remainder augmented-dimensional variables, respectively, and h x (k + 1) and γ(h x (k + 1)) are the observation matrices corresponding to the original variables and the remaining augmented-dimensional variables, respectively. Let the projection operator be where I x and O β are matrices that match the dimensions of the original and remainder augmented variables, respectively. Then,x (k + 1|k + 1) = P e * X(k + 1|k Substituting Equation (48) into Equation (50), we obtain It is found that after projection through the projection matrix, only the original system state estimation value is retained, which reduces the complexity of the algorithm and the actual computational burden, and includes more information from the model, improving the estimation accuracy.

Covariance Matrix Analysis
Bring Equation (48) into Equation (47) and at the same time divide the covariance matrix into blocks according to the original variables and the remaining dimension expansion variables, then write as follows: P x (k + 1|k + 1) P xβ (k + 1|k + 1) P βx (k + 1|k + 1) P β (k + 1|k + 1) where P x (k + 1|k + 1) is the covariance matrix of the original variable, P xβ (k + 1|k + 1) is the covariance matrix of the original variable and the remainder variable, P βx (k + 1|k + 1) is the covariance matrix of the remainder variable and the original variable, and P β (k + 1|k + 1) is the covariance matrix of the remainder variable. Then, the original variable covariance matrix can be calculated by the projection matrix as follows: In the experiment, it was found that K x (k + 1) * γ(h x (k + 1)) * P xβ (k + 1|k) is a positivedefinite matrix. Compared with the original EKF algorithm, the posterior covariance matrix becomes smaller with the application of higher-order residual information, and the algorithm makes full use of the higher-order information of the model. Theoretically, the state estimation accuracy will be higher.

Simulation Experiments
For multi-robot formation systems , the state equation is as follows [24,35]: where (ξ i (k), ψ i (k)) and θ i (k) represent the position and direction of the i-th robot, respectively, and φ i (k) and δ i (k) represent the linear velocity and angular velocity, respectively.
In a multi-robot formation, Γ is a time-varying matrix related to the formation and Γ is the formation; if Γ is an upper triangular matrix, it means that the formation has a fixed movement formation, as shown in Figure 2. If Γ = Γ(k), it means that the formation is time-varying and has a random switching topology, as shown in Figure 3. We use the extended Kalman filter (EKF) and high-order Kalman filter methods to estimate the secondorder Kalman filter algorithm (SOKF), including the remainder extended Kalman filter algorithm (REKF), as well as the second-order remainder Kalman filter algorithm(SEORKF) for comparison.  As can be seen from the above figure, in a robot formation with a fixed formation, the robot maintains the formation at the first moment and continues to move with a unique nonlinear model. However, in a formation with a random switching topology, at the first moment it maintains a fixed formation structure. With the constraints of different communication or task division, the formation at each moment will be different. There are arbitrary formations or individual formations, and they maintain their own nonlinear motion.
The visual measurement is provided by the following equation: where (p i (k), q i (k)) defines the coordinates of the robot in the image plane, (d 1 , d 2 ) are its frame coordinates, z c f is the distance from the optical center of the camera to the robot, and γ u and γ v are variable pixel magnification factors; for a visual tracking system, the feature points are placed on the ceiling, (s x , s y ) are the coordinates of the feature points in the world frame, and (p 0 , q 0 ) is the camera image coordinate principal point; In the simulation, the following parameters of the visual tracking system are adopted: 1050,γ u = 902.13283,γ v = 902.50141,p 0 = 347.20436, q 0 = 284.34705. The process noise covariance matrix isQ i (k + 1) = diag{0.01, 0.01, 0.01} , and the measurement noise covariance matrix is R i (k) = diag{25 2 , 25 2 } (i = 1,2,3).
To illustrate the tracking performance of the proposed filter, we use the root mean square error (MSE) of the three robots. Over 50 Monte-Carlo runs were obtained, and the table below shows the mean squared error for each robot position: i (k)) represents the position estimate of the i-th robot at the m-th Monte Carlo run.

Trajectory Estimation of Multi-Robot in Fixed Formation
In a multi-robot formation, Γ is a time-varying matrix related to the formation, Γ is the formation, and when Γ is an upper triangular matrix, it means that the formation has a fixed movement formation to complete a certain task. Table 1 shows the multi-robot trajectory estimation MSE with fixed formation, where 'Improved' refers to the improved estimation accuracy compared with EKF method.
The original nonlinear filtering algorithms based on Taylor expansion, such as extended Kalman filtering and high-order Kalman filtering methods, directly discard the truncation error. The estimation accuracy has a great influence. Table 2 shows the truncation errors generated by the EKF and second-order Kalman filtering methods for estimating the fixed formation trajectories of multi-robot. It can be seen from the table below that both methods have considerable truncation errors.
In the high-order Kalman filtering method, the remainder variables are introduced into the algorithm to make full use of the nonlinear model information, which can effectively avoid the influence of the truncation error on the estimation results; however, there will always be rounding errors. Here, the rounding error results generated by the first-order remainder extended Kalman filtering algorithm and the second-order remainder Kalman filtering algorithm with a fixed formation of robots are compared and analyzed, as show in Table 3, with the result that the rounding error is much smaller than the truncation error.  As follows, in the fixed formation, Figure 4 shows the real trajectory of the fixed multi-robot formation; Figure 5a,b shows the R1 estimation error in X-Coordinate and Y-Coordinate respectively; Figure 6 shows the histogram of R1 estimation error; Figure 7a,b shows the R2 estimation error in X-Coordinate and Y-Coordinate respectively; Figure 8 shows the histogram of R2 estimation error; Figure 9a,b shows the R3 estimation error in X-Coordinate and Y-Coordinate respectively; Figure 10 shows the histogram of R3 estimation error.

Multi-Robot Formation Trajectory Estimation with Random Topology Structure
In the multi-robot formation, Γ is the time-varying matrix related to the formation, that is, Γ = Γ(k), indicating that the formation has a random switching topology. Table 4 shows MSE for multi-robot trajectory estimation and Table 5 shows the truncation error of trajectory estimation. The truncation error of the trajectory estimation extended Kalman algorithm and the second-order Kalman filtering algorithm in Table 6 shows the rounding error when the trajectory estimation is performed by introducing the remainder variables into the above two algorithms, and the results are similar to the trajectory estimation under the fixed formation.
As follows, in random topology formation, Figure 11 shows the real trajectory; Figure 12a,b shows the R1 estimation error in X-Coordinate and Y-Coordinate respectively; Figure 13 shows the histogram of R1 estimation error; Figure 14a,b shows the R2 estimation error in X-Coordinate and Y-Coordinate respectively; Figure 15 shows the histogram of R2 estimation error; Figure 16a,b shows the R3 estimation error in X-Coordinate and Y-Coordinate respectively; Figure 17 shows the histogram of R3 estimation error.  Figure 11. The real trajectory of the multi-robot swarm in random topology formation.          As shown in the figure above, we use a fixed formation formation and a random topology formation to respectively estimate the actual trajectories of the three robots, use the estimation error of each robot in the x-y direction of the two-dimensional plane as a measure, and use the EKF and the second-order Kalman filtering(SOKF) methods, respectively, to introduce the residual term variable in each of these two methods. From the results, it can be concluded that the estimation accuracy of EKF is the lowest, followed by the second-order Kalman filtering algorithm (SOKF), the remainder extended Kalman filtering Algorithm (REKF), and the remainder second-order Kalman filtering algorithm (SORKF). The other three filtering methods designed according to EKF structure, excepting EKF, significantly improve the estimation accuracy of EKF. Thus, we can see that for the robot in the motion model, the degree of non-linearity of the turning characteristics is relatively high. If first-order Taylor expansion is performed, the information contained in it cannot be completely extracted, which inevitably leads to a decrease in the estimation accuracy. The higher the degree of expansion performed, the higher the robot position estimation accuracy, and the more accurate the robot estimation estimation. In addition, we analyzed the truncation error of the EKF and second-order Kalman filtering(SOKF) algorithms. The truncation error has a great influence on the estimation accuracy. When the nonlinearity is strong, the model information cannot be fully utilized by directly discarding the truncation error, resulting in low estimation accuracy. After the variables are introduced, the rounding error of the estimation results is analyzed, and it is found that the rounding error is much smaller than the truncation error. The Taylor expansion of the nonlinear model and the use of remainder variables can make full use of the model information to improve the estimation accuracy and make the estimation more accurate.

Conclusions and Future Work
This paper studies the multi-robot trajectory estimation problem, and proposes a fusion estimation method based on the high-order Kalman filter algorithm. The extended Kalman filter (EKF) algorithm used in the existing robot position estimation only considers first-order expansion and ignores the high-order information. To solve the problem, a joint trajectory estimation method with a multi-robot formation based on the high-order Kalman filtering method was adopted, the Taylor expansion of the state equation and the observation equation was carried out, and the remainder variables were introduced on this basis, effectively improving the estimation accuracy. Through a simulation, we found that for robot fixed formations and random topology formations the introduction of remainder variables improves the accuracy of position estimation by more than 50%, and in certain cases even more, compared to the EKF algorithm. At the same time, the truncation error in the estimation process of the EKF and the high-order Kalman filter algorithms was analyzed and compared with the rounding error of the estimation algorithm after the introduction of the remainder variable. We found that the rounding error is much smaller than the truncation error. After the remainder variable is introduced during filtering, the algorithm makes full use of the model information, and the estimation accuracy is greatly improved.
In future research, we may encounter more complex models. If the state model has strong nonlinear characteristics, the observation model has super strong nonlinear characteristics, and the modeling error is non Gaussian white noise, our method will be difficult to achieve, and will be completed with the help of the characteristic function filtering method [36,37]; At the same time, for the neural network model, we can also introduce the high-order Kalman filtering method to update the real-time parameters [38]; Finally, for the state estimation of non-cooperative targets, there must be multiple targets and sensors. We can extend the method in this paper to distributed filtering. Similarly, for the distributed model of Federated learning, we will also applicable [39,40].

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Mathematical Symbols
The important mathematical symbols used in the text are described as follows: