A novel non-collision paths planning strategy for multi-manipulator cooperate manufacturing system

Analogy to the definition of Human-robot Interaction (HRI), the case of multiple manipulators with shared workspace, non-simultaneous manufacturing tasks and separate objects is named as multi-manipulator cooperation, which is becoming more widely employed in modern industrial manufacturing system and requiring non-collision path planning as a key issue in terms of safety and efficiency. In this paper, a novel method called Sampling based Position Space Map Search (SbPSMS) method which combines the map search method with the time-sampling based method will be proposed, including a minimum distance prediction method based on PSO-BP neural network for collision detection and two candidate position determination methods for search map establishing of all manipulators. After the specific search map simplification process, the local path fragments during each sampling time interval can be determined via cost function, which will be glued together to generate the final collision-free paths. The simulation results not only show that the PSO-BP hybrid algorithm has more accurate of nearly 2mm than the standard BP neural network in minimum distance prediction, but also demonstrated that our proposal can successfully achieve collision avoidance of dual manipulators system whilst meeting the real-time requirements for multi-manipulator cooperate assembling scenarios. The further satisfactory simulation results of triple manipulators suggesting our algorithm can be extended to applications of multiple manipulators cooperate manufacturing.


Introduction
At present, six degree-of-freedom (6-DOF) industrial manipulators play an increasingly important role in automated manufacturing, due to a number of advantages including the provision of tireless repetitive labour, faster moving speeds and higher accuracy performance [1]. Thus, tasks requiring multiple workers can undoubtedly be accomplished by multiple manipulators cooperation， which means that multiple manipulators respective manufacturing will not only work side by side, but as dyads and teams. Analogy to the definition of Human-robot Interaction (HRI) [2,3], the case of dual or multiple manipulators manufacturing system can be divided into cooperation and collaboration according to the tasks arrange-ment ( Fig.1). For the case of collaborate manufacturing system, all manipulators together with the executed mechanism constitute a complete multi-degree-of-freedom closed-loop or parallel manufacturing system. Thus, the non-collision manufacturing control strategy under this circumstance can be transformed into the currently matured internal obstacle avoidance strategy [4]. However, for the scenario of multiple manipulators cooperate manufacturing system (e.g. assembling as shows in the right figure of Fig.1), problems and limitations arise in small overlapping workspaces accommodating numerous cooperative manipulators, whereby collisions can occur if no related countermeasures are put in place. In this case, finding a reliable strategy of simultaneously collision-free paths planning for all the cooperative manipulators is a sensible choice to overcome the upper-mentioned limitations. Fig.1 The schematic diagram of multi-manipulator cooperate and collaborate manufacturing system. Multi-manipulator collaborate manufacturing

ORIGINAL ARTICLE
system (Left): the case of multiple manipulators with shared workspace, simultaneous manufacturing tasks and common object; Multi-manipulator cooperate manufacturing system (Right): the case of multiple manipulators with shared workspace, non-simultaneous tasks and separate object; 6-DOF industrial manipulator implies a serial chain robot with six revolving joints, which are structurally much more complex than mobile robot. For modelling of the 6-DOF in-dustrial manipulator, a description of the position space, posture space, execution space, and configuration space are given as descripted in Fig.2. is the base coordinate system of the manipulator; is the direction vector of the central axis of the sixth joint;

Relative work
is the position coordinate vector in . All position coordinates that the manipulator EEF can reach within the space set are gathered and designated as the position space. = ( , , ) is the direction vector of the central axis of the sixth joint, which can also be used to represent the pose of manipulator EEF. All possible direction vectors that manipulator EEF can reach are collected in a set, referred to as the posture space. The combination of the position space and posture space constitutes the execution space. The value of each joint can be denoted as ( = 1~6). Then = ( 1 , 2 , 3 , 4 , 5 , 6 ) can be recognized as a point in the 6-DOF configuration space (C-space). The relations between the execution space, C-space, position space, and posture space are shown in (b). Every point in the C-space can determine a only point in execution space via forward kinematic calculation, but every point in execution space can correspond to more than one point in the C-space for the inverse kinematic calculation of 6-DOF manipulator has multiple solutions.
Since single manipulator can be modelled as a particle in its own C-space, the particle-based path planning methods [5] can be used for single manipulator path-planning directly. However, considering the non-collision path planning of multi-manipulator cooperation scenarios, the main encumbrance that prevent us from directly applying particle-based path planning algorithms can be summarized in the following three points: (a) The dimensions of C-space are equivalent to the joint number of the manipulator. As a consequence, using particle-based path planning methods in high-dimensional space will greatly increase computation time and reduce the efficiency of the algorithm [6]. (b) The mapping relationship between C-space and execution space is nonlinear, which means the optimal path in C-space may not be optimal in execution space [7]. As the actual workspace of the manipulator is execution space, finding an optimal execution space path rather than C-space path is the final target of our exertion. (c) The cooperate manipulators are themselves regarded as obstacles to each other. Therefore, all manipulators must be projected into the same C-space. Due to the irregular structures of manipulators, this process is very complicated and timeconsuming [8].
Limited progress has been made for non-collision path planning of multiple cooperating manipulators in C-space. Obstacles were modelled as simple geometries by Jia et al. [9], which will then be projected into the C-space of a manipulator. This method is not effective for multiple manipulators due to their complex structures. Following on from this, Yu et al. [10] divided the C-space into multiple layers based on certain dimensions, and implemented a rapidly created Cspace grid map on each layer, however, due to the huge amount of elements in the maps, the whole algorithm tends to be time-consuming. Li et al. [11] planned collision-free paths for 2D horizontally articulated dual-arm robots by dividing the C-space of each robot into multiple blocks, marking the obstacle space and free space at different time intervals, and mapping the search method to seek out an optimal path in free space. The method was too computationally expensive for real-time applications, and moreover, it is difficult to project a 6-DOF manipulator body into the C-space of another 6-DOF manipulator. Harada et al. [12] proposed a general manipulation planner for a dual-arm industrial manipulator, which combined the C-space of the arms and obstacle space into one overall search map. The method mainly focused on the trajectory arrangement and movement order of each arm to allow a target object to pass from starting point to end point, however, the possibility of collision between the two working arms was not considered.
Hence, how to explore a simple, effective and reliable execution space path planning algorithm has become one of the more and more enthusiastic focuses in recent decades. An online collision-free trajectory generation algorithm for dualarm robots was established by Lee et al. [13], by setting up a Virtual Road Map (VRM) in the execution space, and refreshing the map with a new collision-free path. However, the execution time is long due to difficulties in setting up the VRM, and extending the method to multiple robots is challenging. Larsen et al. [14] divided the working area into different parts, and set up obstacle models within each part. A master-slave method was then used to choose a free path; however dynamic environments cannot be accommodated since the obstacle must be static. Cohen et al. [15] combined the C-space and execution space to build a "motion space", and used the Lazy Weighted A* method to plan a non-collision path in an environment containing N-manipulators, but only path planning for one robot is allowed at a time and the planning process is performed off-line. While the ARA* method was also tested to search for the optimal path in the workspace (execution space) of a manipulator, the study focused only on dual arms performing the same task and did not consider the coordination of different tasks.
Additional methods have been proposed for the non-collision path planning of multiple manipulators. Chiddarwar et al. [16] used the A* method to plan a collision-free path in Cspace without considering the motion of other robots, and a Path Modification Sequences (PMS) method was then proposed to arrange the moving sequence of each robot, resulting in a much longer execution time. Another method proposed by Afaghani et al. [17] used a collision map for detecting the collisions between two robots in order to avoid deadlocks, which can occur if one robot becomes an obstacle to another. Unfortunately, the method simply delays the movement of the robot to avoid the collision. Rodríguez et al. [18] suggested an approach based on a variation of the Probabilistic Road Map, called the Probabilistic Road Map with Obstacles (PRMwO). The method does not exclude collision samples with removable objects, but instead classifies them as collided obstacle(s), and allows the search for free paths b highlighting which objects must be removed from the workspace to make a valid path. While this approach removes any obstacles along the path of the working manipulator, it does not enable the manipulator to actually bypass obstacles. Habibnejad et al. [19] used the Artificial Potential Field (APF) method to plan the paths of multiple cooperative manipulators on mobile bases, however the study focused solely on the path of the end effector (EEF) and did not consider the entire arm, thus lacking important systemic considerations.
A preliminary conclusion can be drawn from the above review that when it comes to the path planning of dual or multiple industrial manipulators, the negative influence, due to the above characteristics such as model complexity, high dimensionality and complex planning space, will become even more critical. Instantly multi-manipulator path planning algorithms can be divided into map search based method, timesampling based method, mathematical model method and others [20]. Map search based method [21][22][23] is a global optimization algorithm which will cost huge time on planning map construction. While time-sampling based method [24][25][26] divide the whole planning period into several orderly but isolated intervals, and splice all the calculated outcomes of each interval into a complete result path or planning map in order.

Paper organization
This paper focus on the simultaneously non-collision path planning for the scenario of multi-manipulator cooperate manufacturing. Under this circumstance, a path planning algorithm that can skirt round the constraints of mechanical structure, simplify the planning space, and determine the paths of all manipulators participating in cooperation system in real time will be the kernel of this paper.
Thus, a novel method called Sampling based Position Space Map Search (SbPSMS) method which combines the map search method with the time-sampling based method will be proposed. Equivalently speaking, a local planning map in each time interval will be built and a local optimal path will be determined according to the local planning map, then all path fragments that be planned in chronological order are glued together to generate collision-free paths, which is very efficient at the cost of abandoning the global optimum.
The first step in building a planning map is modelling and collision detection. In order to further shorten the execution time of these two steps, we try to use neural network algorithm to merge these two steps into one step. Then, two different candidate position determination method as well as specific search map simplification process will be descripted in detail to set up the local planning map of all manipulators. After path planning based on map search algorithms, the remaining work is to repeat the same work during each sampling time interval until all manipulators reaches their respective targets.
In this paper, a brief path planning method review of multiple manipulator cooperation is shown in Section 1. Section 2 provides a detailed framework to describe how to predict the shortest distance between two industrial manipulators in real time based on the PSO-BP neural network. Section 3~5 are the detail content of SbPSMS method. In section 3, two candidate position points selection methods are determined for a single manipulator, providing the nodes of the search map for all the cooperative manipulators. In section 4, we describe the search map simplification method for dual and multiple manipulators. In section 5, formulae for the cost function are introduced as an evaluation index of the optimal path. Section 6 presents the simulation results and relevant analyses. Finally, in section 7, some conclusions are presented.

Background
Usually, in order to determine the potential collision, a safety control strategy of collision detection by comparing the shortest distance between manipulators with threshold is ap-plied. In the early stage, the minimum distance was calculated through the 3D geometric models and many efficient algorithms were well developed [27,28]. But because of its computational cost, hierarchical bounding volume which employs geometrical primitives to approximate the objects is adopted to overcome the shortage in real application, including axis-aligned bounding box [29], object-oriented bounding box [30,31], sphere-swept volume [32,33], ellipsoid [34], sphere [35] and convex hull [36]. Basically, most bounding volumes are struggling to get a trade-off between the model accuracy and algorithm complexity, computation efficiency and the tightness of fitting the underlying object. Due to the rotation invariance, sphere is a popular choice of bounding volume. But the difficulties lie in the optimization of the radii and the number of the sphere in order to obtain a relative high accuracy. There are also other techniques that focus not only on the minimum distance between the objects, but also take the system state into consideration to evaluate risk index [34]. Through the methodology described above to approximate the objects, minimum distance is equivalent to calculate the shortest distance among geometrical primitives. But the algorithm calculating the shortest distance between two geometry objects is usually complex to implement and many different cases need to consider.
In the developed approach, to simplify the implementation of the algorithm, potential collision detection of the two manipulators is predicted by adopting PSO-BP neural network method. BP neural network is very powerful in nonlinear mapping and can approximates nonlinear function in any precision [37]. Besides it's popular in prediction, its drawback is easy to fall into local minimum. Nevertheless, PSO algorithm perform well in global optimization [38]. The improved PSO with two optimization stages integrating k-mean clustering to keep population variety and described in the subsequent section is applied before training BP neural network. With the optimized values as initial weights, the BP neural network can decrease the training time and generate a better training result of accuracy than the standard BP neural network with respect to given training and test data set, as investigated in [39,40].

Modelling & Minimum distance calculation
The manipulators are composed of irregular and complex parts, which results in difficulties in computing the minimum distance. In this paper, sphere-swept volume is applied to approximate the linkages and joints in order to calculate the minimum distance between each two manipulators for the neural network training process, as shown in figure 3.
According to the modelling, the minimum distance calculation process is converted to determine the shortest distance between the surfaces of two capsule-like geometrical objects. Assume 1 , 2 are the two points sets in Cartesian space, 1 ( 1 , 1 , 1 ) ∈ 1 , 2 ( 2 , 2 , 2 ) ∈ 2 , then the minimum distance between 1 , 2 can be expressed as: Right: Minimum distance of two capsule models In figure 3, define the surfaces of manipulator A and B by ∈ 3 , ∈ 3 respectively; the center line segments of the cylinders and the center points of the spheres are defined by ∈ 3 , ∈ 3 ; , are the radiis of the hemispheres ( ∈ (1,2,3,4)) .Thus, the minimum distance function can be summarized as [41]: This fact indicates that the complex calculation of two cylinder surfaces is replaced by the calculation of minimum distance of the two centre line segments of the cylinders and then subtracting the radii of hemispheres A and B [41]. Meanwhile, in order to determine the position and orientation of line segment, forward kinematic [42] is applied which is complex and time consuming. This paper try to propose an algorithm constructed by PSO-BP neural network [43] to predict the minimum distance with respect to the given joint values of all cooperated manipulators.

Detail of PSO-BP neural network
In this section, a hybrid method of PSO and BP neural network is developed to predict the shortest distance in a given position and posture. In the presented hybrid method, the PSO algorithm is utilized to optimize the initial connection weights of neurons and threshold values of network. With the introduction of an extra learning item evaluated by the fitness-distance-ratio [44] and k-mean clustering, the PSO can avoid premature. In this case, the hybrid algorithm combining the PSO and BP neural network has higher generalization and accuracy [39,40].
1) PSO algorithm PSO is a global optimization algorithm based on swarm intelligence and evolutionary computation, which searches the global optimization solution via the competition and learning among all particles. The original PSO formulas proposed by Kennedy and Eberhart [45] were modified [46] with inertial parameter which was empirically improve the overall performance.
The core idea of PSO is define every particle in the swarm with size of as a potential solution problem in a -dimensional space, with the i-th particle represented as = ( 1 , 2 , 3 , … , ). Each particle has two important properties of position and velocity, it preserves its previous best position ( ) and velocity ( ) along each dimension, represented as: = ( 1 , 2 , 3 , … , ) = ( 1 , 2 , 3 , … , ) In each generation, the of the particles with the best fitness named as gbest which can be denoted as = ( 1 , 2 , 3 , … , ), together with the pbest of the current particle, are used to adjust the velocities of each dimension. The updated velocities are then used to generate a new position for the particle. In mathematics, the position and velocity of each particle are computed by the following update formulae: where is the serial number of the particle; is the ℎ dimension of the particle; ω = − − × is the inertial weight and dynamically changes with respect to iteration; = 1,2,3, … , ( is the maximum iteration of PSO) is the current iteration of the algorithm; 1 , 2 are learning factors or acceleration coefficients determining the relative influence of the cognitive and social components, which are usually equal to 2; 1 , 2 ∈ [0,1] are random numbers.
2) Fitness-distance-ratio In order to avoid the premature convergence observed in many applications of PSO, Peram T, Veeramachaneni K and Mohan C K [44] proposed FDR-PSO in which the particle moves towards nearby particles with higher fitness. By using the ratio of the fitness difference to the one-dimension distance, the velocity of the ℎ particle in ℎ dimension is also updated via a nearby best particle with the highest ratio, denoted as : In this case, the particle's velocity is influenced by three particles: previous best particle, global best particle and the nearby best particle. Considering the item of nearby particle, the formulae (4) can be changed into: where 3 is also a constant quantity representing the relative influence of the nearby best particle; 3 ∈ [0,1] is a random number.
3) -mean clustering -mean clustering algorithm is an unsupervised learning algorithm that intends to divide the particles into subgroups with the shortest Euclidean distance to the clustering centroids. -mean clustering method is introduced in order to improve the globalization optimization and population diversity of PSO algorithm. This algorithm aims to minimize the cost function with respect to the cluster centroids using following formulas: where ‖… ‖ denotes the Euclidean distance and is the number of the particles; represents the cluster centroid of the ℎ subgroup; is the number of the particles in the ℎ group. 4) BP neural network BP neural network is an artificial neural network that simulates the cerebrum information processing. BP neural network with hidden layer can approximate linear or nonlinear function with any precision. As the consequence of its powerful ability of nonlinear mapping, it has been widely used in pattern recognition, image processing, automatic control and some other field.
The basic structure of BP neural network is input layer, hidden layer and output layer. Input layer accepts the training data set, namely 12 joint values in this paper; hidden layer maps the input data to the output layer through connection weights and activation functions; output layer linear-combine the outputs of hidden layer to computes the predicted result. The model of BP neural network with layers and one output is shown in figure 4.

Fig.4 Model of BP neural network
For a neural network with one output, the mean square error of all training sample is used to evaluates the performance of neural network and can be expressed in formula (9).
In this learning algorithm, the training process of the network terminates when the network error is less than the given tolerance or reach the maximum epoch. Back propagation (BP) algorithm and gradient descent are combined to update the connection weights of neuron throughout the whole training process. In the output layer, the error can be denoted as: where is the error of the output layer; is the output of the network and is the training data. When the error backs propagation into the hidden layers, it needs to take the connection weights into consideration as shown in equation (11).
In the formula, refers to the error of the ℎ hidden layer( = 2, … , − 1 ); +1 is the weight matrix of the ( + 1) ℎ layer; is the input of ℎ hidden layer; ( ) ′ is the derivative of activation function with respect to and ° means Hadamard product of matrix. With the gradient descent, the weights of neural network can be updated using the following formulas: is the learning rate which has significant influence on the convergence speed and the accuracy of the network. But usually it is difficult to assign a specific value from the beginning in real application. Larger numbers will result in oscillations in the minimum value, while a small number will cause a long training time. Using adaptive learning rate [47] with respect to the error difference is a feasible solution of making a balance between these two aspects:

5) PSO-BP neural network
Standard BP neural network uses gradient descent for training which indicates that the training result and accuracy are sensitive to the initial weight values and threshold. It is a feasible way to introduce PSO into BP neural network. The basic idea is to use the PSO algorithm to optimize the weight values, and then the optimized solution of the PSO is transferred to BP neural network as its initial weight. A PSO algorithm containing two stages is utilized to optimize the initial weights. In the first stage, all the particles that are divided into subgroups with -mean clustering search the global optimization together until the maximum iteration; in the second stage, PSO algorithm works on a new population that is composed of global best particle of subgroups. After these two steps, the searching space of BP neural network narrows down and the hybrid method has better performance in prediction than standard BP neural network. The details are shown in the Fig.5.

Candidate position determination for single manipulator
The candidate positions for a single manipulator consist of all possible positions of the EEF that the manipulator can reach within the position space during the next scanning pe-riod. Determining these candidate positions is the cornerstone of proposal method, and the key to building a simple and effective position space search map (PSSM).
During the movement process of an industrial manipulator, the absolute value of EEF speed is generally held uniform throughout the entire movement period to maintain stability of the system, which means the step length of the ma-nipulator EEF in every scanning interval must be constant. Therefore, the set of all candidate positions of a single manipulator constitutes a spherical surface, in which the current space position point of the EEF, , is the sphere center and is the radius, and we call this set as Search Feasible Region (SFR). Theoretically, there can be an infinite number of candidate position points within SFR, however infinite continuous points cannot be used to build a node map, thus selecting a specific number of points is necessary. Here, we propose finding three non-coplanar vectors as the reference direction vectors and determining all candidate positions of a single manipulator, referring to the method of establishing a "unit cell" in chemistry. In this section, two separate methods are proposed to find the three reference direction vectors: the first is based on the base coordinate system of each manipulator, and the second considers the current velocity of each manipulator EEF.
However, before talking about how to determine the three reference direction vectors, the concepts of main movement axis, secondary movement axis, displacement axis, and main movement plane are proposed.
Definition: Assuming that a point moves from start position to endpoint under a coordinate system as shown in Fig. 4, and the coordinate system can be arbitrary, the vector from to is = ( , , ), which indicates the size and direction of the movement. Displacements in the positive , and directions are , , and . The values ‖ ‖, ‖ ‖, and ‖ ‖ can be compared and sorted in a descending order, such that the axis corresponding to the minimum value is set as the displacement axis. Moreover, the main movement plane is defined as the plane composed of the other two axes, whereby the main axis of the plane is the axis corresponding to the maximum number, and the remaining axis is the secondary movement axis.
As an example scenario in Fig.6, if ‖ ‖ ≥ ‖ ‖ ≥ ‖ ‖ then the -axis is the main movement axis, theaxis is the secondary movement axis, the -axis is the displacement axis, and the -plane is the main movement plane.

Coordinate system based method
The principle of coordinate system based method is to project the base coordinate system of the single manipulator onto the global coordinate system to establish the three reasonable reference direction vectors by using the motion attribute of the manipulator while moving along a given path in Cartesian space. In a multi-manipulator system, the motion of every single manipulator can be calibrated conveniently by each manipulator's own base coordinate system, but the relative motion among multiple manipulators can only be calibrated under the global coordinate system. Besides, when applying the PSMS method on the non-collision path planning of multi-manipulator system, we need to combine the motion of every single manipulator and that of all the manipulators together to judge how to reduce the number of nodes in the PSSM in order to increase the efficiency of the algorithm. The common base coordinate system for a single manipulator is shown as Fig.7. When the manipulator is in its initial position (which means that all the values of six joints are zero), the − plane of the base coordinate system is the bottom surface of the manipulator base, the origin point of base coordinates system is the center point of the bottom surface of the base. From the origin point, the axis of the base coordinates system is the axis perpendicular to the bottom surface of the base and point to the direction of the manipulator mechanical arm, the axis is the axis in the − plane and point to the position of EEF, and then the axis can be determined by the right-hand rule. Then the base coordinate system of the manipulator is built up. is the global coordinate system. is the global coordinate system as well as is the base coordinate system. The posture of base coordinate system is arbitrary. Take the scenario of Fig.4 for example, the reference direction vector is the projection vector of axis on the displacement axis , is the projection vector of axis on the main movement plane , and is close to axis and perpendicular to both the reference direction vector and .
The execution steps of the coordinate system based method are shown as follows: Step 1: determine the main movement axis, secondary movement axis, displacement axis, and main movement plane of the manipulator under the global coordinate system; Step 2: move the origin point of the global coordinate system to the origin point of the base system coordinate system. Since the base coordinate system can be arbitrary, the base coordinate system and the global coordinate system may not coincide; Step 3: project the axis of the base coordinate system on the displacement axis to be the reference direction vectors , project axis of the base coordinate system onto the main movement plane to be the reference direction vectors , then the reference direction vectors will be determined via the following formula: = × ( × ) (15) which means the reference direction vector is close to the axis and perpendicular to both the and . However, in some cases like the axis of the base coordinate system will in the main movement plane which means the axis is perpendicular to the displacement axis, or the axis coincides with the displacement axis which means the projection of axis on the main movement plane is the origin point itself, the aforementioned process will lose efficacy. In these cases, we only need to interchange the axis and the axis of the base coordinate system, and then the three reference direction vectors can be got according to the processes aforementioned.
It is obvious that these three reference direction vectors ( , , ) are perpendicular to each other as the base coordinate system and global coordinate system do.

Velocity based method
According to particle kinematics and dynamics, when the manipulator EEF passes point at time instant with speed , it will move at a speed which is as close as possible to in the next time instant according to the principle of minimum energy. When the current velocity of the EEF is not parallel to any of the three axes of the manipulator base coordinate system, the direction of the velocity vector, as determined by the base coordinate method, is not necessarily consistent with the current velocity direction. Therefore, we propose a new method of three reference direction vectors determination based on the current velocity.
Firstly, a new coordinate system called velocity coordinate system will be set up to calibrate the motion of every single manipulator, in which is the current position point of manipulator EEF. As shown on the right half side of Fig.  6, the velocity of the manipulator EEF at point is at the current time instant. Thus, the direction of vector can be taken as the direction of the reference -axis. Then, the EEF position and the base of the manipulator can be connected to construct the space vector ⃗⃗⃗⃗⃗⃗ , which is projected onto the bottom surface of the manipulator base in the same coordinate system and serves as our reference -axis. Finally, the -axis of the candidate position points can be calculated from the -and -axis according to: = × (16) Then how to determine the three reference direction vectors are the same as those presented in Section 4.1.
After setting up the three reference direction vectors, the determination of candidate positions for a single manipulator will be an easy job. As shown in Fig. 8, we assume the coordinate of the center point (black point) of the EEF at time is = ( , , ), and ∆ = (∆ , ∆ , ∆ ) is the step gain required for the EEF to reach the target position along the three reference direction vectors, where ∆ = ∆ = ∆ . That is to say, if the axis coordinate of point at time instant is , a possible coordinate for +1 must be one of the following: − ∆ , , and + ∆ . The same can be applied for +1 and +1 .
According to the assumptions presented above, 27 different coordinate combinations can be obtained using the following formula: = ( +1 , +1 , +1 ) where is the set flag and the operator "×" denotes "combination".
These coordinate combinations can set up as 27 points (marked ~; yellow points in Fig. 8) in Cartesian space, and form a cube with side length 2∆ about point P. The 27 points are referred to as the candidate positions of point P. However, inconsistencies exist between the displacement values from point P to the 27 candidate positions, which results in difficulties in setting up the cost function.
In order to make every displacement fixed, an inscribed sphere of radius ∆ is placed into each cube, as shown in Fig. 8. From point P to ~, there will be 27 radials intersecting the inscribed sphere resulting in 27 intersection points marked as ′~′ (green points in Fig. 8). The intersection points are the candidate positions with constant displacement in Cartesian space for a single manipulator during the next scanning period. After the process above, 27 choices are available for one 6-DOF industrial manipulator, and can be uniformly marked in matrix . The 27 candidate positions ( ) also constitute the PSSM of a single manipulator at time .

Search map simplification process
For multi-manipulator system, after setting up the "unit cell" structure of a single manipulator, the next step is to set up a reasonable PSSM for the entire system. Since the number of nodes in the "unit cell" can be reduced according to the physical and kinematic characteristics of the multi-manipulator system, a system of dual manipulators and one with no less than three manipulators will have different physical space position features. Thus, the search map simplification (SMS) of these two systems must be described separately.

SMS for dual-manipulator system
The two manipulators are marked as A and B in Fig. 9. The global coordinate system of the dual-manipulator is = ( , , ) , while = ( , , ) and = ( , , ) are the base coordinate systems of manipulators A and B, respectively. Since the two manipulators are mounted on the same plane, we can define the -axis of each base coordinate system in the same direction as the -axis of the global coordinate system. As mentioned in Section 4, the candidate positions of manipulators A and B can be expressed as and , with each containing 27 elements in their own three reference direction vectors. Every element represents a position point under global coordinate system, can be set as the horizontal coordinate and as the vertical coordinate, in order to compose a PSSM of 27 × 27 nodes. Finding an optimal node in this map can provide an optimal solution for the dual-manipulator cooperation in the next step. If the time variable is taken into consideration, the time dimension can be added into the PSSM to produce a time-position space map (TPSM). Thus, planning of a collision-free path can be converted into finding an optimal path in the TPSM (Fig.10). However, a PSSM with 27 × 27 nodes is not efficient enough to meet real-time requirements. Hence, additional measures must be taken to reduce the number of nodes. One possible method is to reduce the number of candidate positions of a single manipulator.
Every element of the candidate positions for a single manipulator consists of 3 axis coordinates according to their own reference direction vectors, in which ( , , ) represents three non-coplanar reference direction vectors.
The above example can be used to demonstrate how to reduce the elements of the PSSM. The search choices can be reduced along the -, -, and -axis based on the following movement analyses: a) Along -axis Based on the previous assumption, the coordinate value of point along the -axis at time is and the physical meaning of three possible values of +1 ( − ∆ , and + ∆ ) represent three choices: "step backward towards the starting point", "remain stationary " and "step forward towards the target". Since the "step backward to the starting point" option does not promote efficient movement of a manipulator, this choice can be abandoned. Thus, movement along main movement axis can be reduced from 3 choices to 2.
b) Along -axis The movement in the main movement plane of manipulator A (Fig. 6) is shown in Fig.11.   Fig. 11 Movement in the main movement plane of manipulator A "Movement towards the manipulator base" (green arrow) and "remain stationary" options are the only two reasonable avoidance movements along -axis, since movement towards the base of the other manipulator would increase the possibility of collision. In the example scenario, manipulator A abandons moving along positive -axis, while manipulator B abandons along the negative -axis. In a similar way, the choice of movement along the secondary movement axis can also be reduced from 3 to 2 choices. c) Along -axis In the example scenario, the -axes of dual manipulators are the same and both the same with the -axis of global coordinate system, the coordinate values at time for the two manipulators under the global coordinate system are and , respectively. Two boundary values ε and -ε are set, where ε > 0. If − > , the reasonable movements of manipulators A and B along their respective -axes are in the positive (towards the target point) and negative (away from the target point) directions, respectively. If − < − , the reasonable movements of manipulators A and B along their respective -axes are in the negative and positive directions, respectively. If ‖ − ‖ ≤ , manipulators A and B only need to move in the main movement plane of each manipulator.
In summary, movement choices along the displacement axis can be reduced from 3 to either 2 or 1, depending on the actual position coordinates.
The results of the choice reduction for manipulators A and B of the example scenario are summarized in Table 1.  However, if the -axes of dual manipulators are both along the a same axis but the axis is not the -axis of global coordinate system, then we take the along-axis of the global coordinate system as a new " -axis" and do as step c; if the -axes of dual manipulators are not along the a same axis, then manipulator A and B only need to move in the main movement plane of each manipulator.
The above analysis shows that the number of elements in the PSSM can be reduced from 27 × 27 to 7 × 7 , or in some cases 3 × 3.

SMS for multiple manipulators
According to the analysis in Section 5.1, additional processing procedures should be performed to simplify the PSSM in order to meet the requirements for real-time calcu-lation. In this section, the PSMS method is applied to the cooperation of multiple manipulators in the same way it can be applied to dual manipulator system. manipulators are marked as A, B, …, N. The transformation relationships of the position coordinates for each coordinate system are shown in Fig. 12. is the transformation matrix of coordinate system in relation to the global coordinate system , and similarly, to . According to Section 4, one manipulator has 27 candidate position points. Hence, for manipulators, the whole system produces a PSSM with 27 nodes, which undoubtedly hinders practical application.
If the coordinates for the starting points , , … , and target points , , … , of each manipulator under the global coordinate system are known, the main movement axis, secondary movement axis, and displacement axis of each manipulator can be identified. Thus, the number of elements in the PSSM can be reduced following the same logic used in Section 4.1.
Upon obtaining the coordinates for the starting points and target points of all manipulators, the displacement axes of all robots can be determined, and all manipulators with the same displacement axis can be collected as a set. For example, all manipulators with the same displacement axis are collected as the set of " displacement axis manipulators" (whose total number is ), and similarly, the set of " displacement axis manipulators" (whose total number is ) and " displacement axis manipulators" (whose total number is ( − − )) can be collected.
For each manipulator set, the candidate positions of each manipulator are defined following a certain set of rules: (a) Along the -axis there are two candidate positions ("remain stationary" or "move one step towards its target").
(b) Along the -axis there are two candidate positions ("remain stationary" or "move one step towards its base").
(c) Along the -axis is slightly more complex and must be classified based on the real-time situation. All the -axis coordinate values under the global coordinate system are sorted in descending order and be divided into three subsets, the "positive set", "negative set", and "zero set". Two boundary values are set, ε and −ε where ε > 0, such that values greater than ε are grouped into the "positive set", values lower than −ε are gathered to form the "negative set", and the remaining values constitute the "zero set". Then, every manipulator in the "positive set" has two candidate positions, remain stationary or move one step towards the positive direction; every manipulator in the "zero set" has one candidate position, remain stationary; and every manipulator in the "negative set" has two candidate positions, remain stationary or move one step in the negative direction.
Since remaining stationary along all the three reference direction vectors is not desirable in terms of maximizing manipulator movement efficiency, the nodes represented by remain stationary along all three reference direction vectors are removed. The above process can thus reduce the number of elements from 27 to 7 , or in some cases even 3 .
For manipulators in different sets, reductions in the number of candidate positions are independent. However, the above rules for element reduction may not be applicable under special conditions, for instance when one manipulator moves to its workspace boundary. In this situation, it is necessary to expand its number of candidate positions from 3 or 7 to 27 in order to find an optimal result for every manipulator during its next movement.
Next, we display the reduction results of the number of PSSM nodes for the cooperation of three manipulators. Assuming all -axes are -axes, all -axes are -axes, and all -axes are -axes under the global coordinate system, the results of the element reduction are given in Table 2.

Cost function and optimal solution selection
Formula representing the cost function for multiple manipulators scenarios are presented in this section. Assuming there exist manipulators, for the -th manipulator, the number of nodes in the PSSM at the current time instant is and the configuration of the -th node at current time instant are = ( 1 , 2 , 3 , 4 , 5 , 6 ) (18) The configuration of the -th node at time instant ( + 1) is + = ( 1 +1 , 2 +1 , 3 +1 , 4 +1 , 5 +1 , 6 +1 ) (19) The EEF position of the -th node at time instant ( + 1) under global coordinate system is + = ( +1 , +1 , +1 ) (20) The target position of the -th node under global coordinate system is = ( , , ) (21) Finally, the cost function using the classical A* algorithm is ( ) = ( ) + ( ) (22) where ( ) is the total cost value of the current node in the search map at the current time instant, ( ) is the cost value from the current node to candidate node, and ( ) represents the cost value from the current node to target node. The motion of a manipulator can be expressed in both C-space and execution space, therefore the cost value can also be calculated in these two spaces. However, the distances between the candidate nodes and the current node are all the same in the execution space, such that ( ) in the configuration space and ( ) in execution space can be calculated individually. The cost value for the -th node in the PSSM of the -th manipulator can be calculated as where ‖ ‖ is the secondary norm of vector , and and are cost coefficients, which can be predefined. The total cost value of the -th node for manipulators at the current time instant is (24) where = ( 1 , 2 , … , ) are the Boolean values (0 or 1) that determine whether the corresponding manipulator will execute the proposed path planning.
Further, the selection process of the optimal solutions can be defined as the selection of optimal nodes in the PSSM for all manipulators, and summarized as the following mathematical model: The dimension of each particle depending on the number of the connection weights and threshold can be calculated by the formula: = 0 1 + 1 + 1 2 + 2 + 2 3 + 3 .
And the population size of the particles is set to 60 which are divided into 3 subgroups to keep the swarm variety. The constant values in formula (6) are set to 1 = 2.8, 2 = 1.3, 3 = 2. As for the iteration of the PSO, in order to get a tradeoff between accuracy and computational cost, the values of 1 and 1 are set to 3,000 and 1,000 respectively. In order to eliminate variable type difference and saturation caused by large range of the input, the input data should be normalized into interval [-1,1] before training. The following formula is employed to pre-process the data.
where , are the minimum and maximum value of the inputs with respect to the exactly same feature, namely the values of the same joint in simulation and represents the training data.

Simulation result
A training data set with size of 300,000 samples in which each sample is composed of 12 values of joints and the minimum distance is used for training PSO-BP neural network, and conduct 10 trials to eliminate the random factors for a general result. The optimization process with respect to the two stages of PSO is showed in figure 13 and the training of BP neural network with optimized weights is presented in figure 14. In figure 13, it can be seen that the fitness of the particle decreased with the iteration increasing which shows a convergence process. The subfigure (a)(b)(c) corresponding to the three subgroups of the particles in the first optimization stage show that there's vibration of the fitness that does not appear in the subfigure (d). Because the particle with best fitness in a subgroup may be grouped into another one, and thus the best fitness of one of the subgroups might change and be larger than the original one which indicates the k-mean clustering algorithm works well in keeping the population diver-sity of the particles. In order to get a better optimization solution, the particles with best fitness in the subgroups form a new swarm used in the second stage. Subfigure (d) demonstrates that a relative optimal result is obtained in the second stage.

Fig.14 Training of BP neural network
The performance of training results is evaluated on a test set with a size of 5,000 through MSE and MAE from statistic perspective. Table 3 lists the MSE and MAE values of each trial respectively.  Table 3 demonstrates that the PSO-BP hybrid method gets a more accurate result than the standard BP neural network with 1.71 less average MAE and 64.94 2 less average MSE in ten trials. Considering the computational cost, if a larger size of population of particles and more iterations are applied in PSO, a better solution is supposed to attain in a high dimensional problem which is 393 in this paper.
The prediction results of the best trials which is the 3rd trial in hybrid method and 6th trial in BP neural network are selected to showed in figure 15. And it can be seen obviously that errors in (b) are larger and oscillate more intensely than those in (a). No matter what the algorithm is, it is impossible to predict the distance exactly without any error. Since when the line segments are in different relative position and orientation the computation formula is different in mathematics, the learning algorithm needs to compromise different situa-tion. Besides, the different radii of sphere-swept volume approximated the manipulator can lead to error while the learning algorithm intends to make general prediction, in which it needs to get a tradeoff among different radii.

Simulation of path planning
In this paper, two different methods were proposed to set up the PSSM for systems with multiple cooperating manipulators. In this section, the two methods are first tested for the dual manipulator scenario, and the simulation results are compared. Subsequently, the scenario of a triple manipulator system is used as an example to show how the method works in the situation of multiple manipulators. Finally, the simulation results are analyzed and some conclusions are drawn from the results. Moreover, the execution time of each method is assessed to determine whether it is possible meet real-time requirements.
Assuming manipulators in a system of multiple manipulator cooperating, all 6-DOF industrial manipulators have common D-H parameters, which are presented in Table 4.
The simulation results for dual-manipulator and triple-manipulator system will be analysed individually. For more than three manipulators, the results were relatively similar to those obtained for the triple-manipulator system, except for increases in number of the PSSM elements and execution time. All the algorithms were programmed in C++ and run on a 64bit Windows operation system with an i7-4790 CPU GHz with 8 GB RAM.

Simulation of dual-manipulator cooperate manufacturing system
For manipulators A and B, we can define 2 = 25 mm, 3 = 315 mm, 4 = 35 mm, and 4 = 365 mm (i=1,2), where A and B were given the same parameters and positioned face-to-face at a distance of = 800 mm. The minimum tolerance for the distance between the two manipulators was set as = 20 mm . The global coordinate system originated at the base of B, which means the coordinates of the base of B were (0,0,0) and for base of A were (0, −800,0). The step length of each manipulator and boundary values were predefined as = 10 mm and ε = 40 mm, respectively.
Three simple simulation scenarios were established as follows: ( 1) Results of simulation using coordinate system based method According to the principles of the coordinate system based method, the main movement axis, secondary movement axis, displacement axis, and main movement plane in the above three scenarios are the -axis, -axis, -axis and plane.
Results of the simulation using the coordinate system based method for each of the above scenarios are presented in Fig. 16. Movements along the --, and -axis for both manipulators are shown in Fig.17. Moreover, the paths of both manipulators are clearly divided into three distinct parts, marked as 1, 2, and 3 (Fig. 17) and named the "original path period", "avoidance period", and "forward period". The "original path period" indicates the manipulator moved along its original path without deviation, the "avoidance period" is the case where a manipulator path was re-planned in order to avoid a collision, and the "forward period" is the last situation for which the manipulator moved toward its target along a straight line.
For any scenario, it can be observed that when a collision is about to occur, it will be avoided by manipulator A in the direction and avoided by B in the + direction, indicating that both manipulators move toward their respective bases. In this way, the two manipulators bypass the Possible Collision Area (PCA) and achieve the goal of collision avoidance. However, the coordinate values of the -axis stay constant during the "avoidance period" for all the three scenarios. The reason for this is that for each of the three scenarios, when the two manipulators moved into the PCA, the coordinate values of the -axis of the two manipulators and met the condition of ‖ − ‖ ≤ , as described in Section 4.1. Here, simulation d was set up to observe how the coordinate values of the -axis would change when the condition ‖ − ‖ ≤ is untenable.
(d) Manipulator A moves from = (400, −350, 400) to = (−400, −350, 600) , while manipulator B simultaneously moves from = (−400, −450, 500) to = (400, −450, 500) . If no collision occurs, A will move from the start point to target point in a straight line of displacement along the main movement axis and displacement axis ( -and -axis), while B will move from the start point to target point in a straight line of displacement along the main movement axis ( -axis). According to the results presented in Fig. 18, when the two manipulators move into the PCA, ‖ − ‖ > 100 > , therefore B has the choice of "+ " and "remain stationary" and A has the choice of "− " and "remain stationary", according to the principles of the PSMS method. The results show that remained stationary while changed along + direction in the "avoid period".
Real-time execution is also a key assessment index of our algorithms. In order to get the per step execution time of the PSMS method, manipulator A was driven move from = (400, −450,600) to = (−400, −450,600) and B was moved from = (−400,450,600) to = (400,450,600) in the base coordinate system of each manipulator, simultaneously. Then, the base distance was changed from 1000 mm to 900 mm to 800 mm , and each change was tested 5 times. The simulation results are summarized in Fig. 19. The test results of execution time of whole path planning period are summarized in Table 5. As previously stated, the entire movement period can be divided into the "original period", "avoidance period" and "forward period", and a different path planning algorithm is required for each period. The execution times of each period are defined by the variables , , and , and by default, the three variables remain constant for each period. Fig.20 shows the movements of manipulators A and B (red and blue, respectively) along the -axis for the previously defined scenarios, and each of the three periods is demarcated.   Table 5, we can see that the errors are within the allowable range. The results suggest that path planning using the PSMS method for dual manipulators can be executed in less than 3 ms, and can therefore fully meet real-time requirements.
2) Results of simulation using velocity based method The simulation results for path planning, using the current velocity based method under scenarios a, b, and c, are presented in Fig. 21. The movements along the -, -, and -axis for both manipulators are shown in Fig. 22. From the three images shown in the lower half of Fig. 21, it can be observed that regardless of scenario, when a collision is about to occur, both A and B can plan a new path with different geometries to those presented in Fig. 16. However, A will usually avoid a collision by moving in thedirection, while B will move in the + direction, indicating both manipulators will move towards their own base. In this way, the two manipulators bypass the PCA and achieve the goal of collision avoidance. The coordinate values of the -axis remained constant during the "avoidance period" for all three scenarios, which were the same as those in the coordinate system based method. The simulation results (Fig. 23) show movements along -, -, and -axis under scenario d using the velocity-based method, clearly shown to be different along the -axis to those using the coordinate system based method (Fig. 18), mainly due to differences in the principles of each method.
To examine the per step execution time using the velocity based method, the same simulation scenario was established as for the coordinate system based method (Fig. 19). The simulation results are summarized in Fig. 24.
The results for execution time over the whole path-planning period are summarized in Table 6.  Table 6, errors were within the allowable range. The results suggest that path planning using the PSMS method for dual manipulators can be executed in less than 4 ms, which can fully meet real-time requirements.

1) Comparison of simulation results
Based on the analysis and results of the previous sections, both methods can be used to plan non-collision paths for dual-manipulator cooperation systems in real time. However, a number of similarities and differences can be observed upon carefully comparing of the two methods: (a) The simulated paths of both methods are not smooth throughout, and there were "twists and turns" in the geometries. One reason may be that all nodes in the search map based on the PSMS method are special positions within the SFR. Another may be the influence of step length which will be simulated and analysed in Fig.26.     can be seen that the paths generated using the velocity based method had larger displacements along -and -axis than the coordinate system based method, however the planned stepped within the "avoidance period" were smaller by using velocity based method than that of the coordinate system based method.
(c) The simulated paths of the velocity based method (Figs. 21~25) had more apparent "twists and turns" than those of the coordinate system based method (Figs. 16~20). This is due to changes in the reference direction vectors throughout the entire planning period in the velocity-based method, whereas the reference direction vectors remain constant in the coordinate-based method.
In order to clearly identify the influence of step length on the flexibility of the simulated paths under different methods, simulations were performed using both methods according to scenario c ( Fig.19 and Fig.24), with step lengths of 10 mm, 5 mm, and 1 mm. The results are shown in Fig. 26. While a decrease in step length resulted in a smoother path using the coordinate system based method, there was little effect shown on paths generated using the current velocity based method (Fig. 26). The reason is not totally clear, however, one possibility is that the reference direction vectors of the velocity based method changes continuously throughout the planning period whereas the reference direction vectors of the coordinate based method remain constant.
Furthermore, the operating time required for the velocitybased method is longer than the coordinate system based method. As such, the coordinate system based method is more effective and demonstrates greater stability in comparison to the velocity-based method. Therefore, the coordinate system based method will be used to demonstrate use of the PSMS method for triple manipulator cooperation.  . The minimum tolerance for the distance between two manipulators was set as = 50

Simulation of triple-manipulator manufacturing system
. The global coordinate system originated at the base of B, which means the coordinates of the base of A, B, and C were (900,0,0), (0,0,0), and (450,450√3, 0).
During The simulation results based on the above scenario are presented in Fig. 28.
When a collision is about to occur, all manipulators move towards their respective bases, successfully bypassing the PCA and achieving the goal of collision avoidance, as shown in Fig. 23.
Again, to examine the per step execution time using the PSMS method for triple-manipulator cooperation, the same three simulation scenarios were used for which only the base distance D was varied from 950 mm to 900 mm to 850 mm, and each simulation was tested 5 times. Execution times for each scenario over the entire path-planning period are summarized in Table 7. Since the planning process of triple-manipulator cooperation is complex, the complete period of path planning can be divided into four parts: period of simultaneous triple-manipulator avoidance; period of simultaneous dual-manipulator avoidance; period of single manipulator avoidance; and noncollision period. For each of the defined periods, a distinct solution to the algorithm is required, which will lead to different execution times. As such, it is difficult to determine the exact execution time for each period. However, using the timer program embedded within the algorithm, it was possible to determine average values of execution time for triple manipulators, which were are all less than 10 ms.

Conclusion & Future research
In this paper, a novel path planning method called SbPSMS was proposed, which can successfully facilitate the planning of non-collision paths for all manipulators working under multiple-manipulator cooperate manufacturing scenarios according to the real-time circumstance at any time while all manipulators are working along their own tasks without pause.
Firstly, a method based on PSO-BP neural network to predict the potential collision of the dual manipulators system is presented. The hybrid algorithm utilizes PSO algorithm to optimize BP neural network initial weights to obtain a better result. Finally, the BP neural network outputs the predicted minimum distance and potential collision is detected if the minimum distance is less than the safety distance threshold. Through the given training data set, the MAE of the predicted minimum distance on the test set is about 1.71mm less than standard BP neural network.
Then, two different methods were presented to build up the three reference direction vectors in order to determine the candidate positions of all manipulators: coordinate system based method and velocity based method. Further, the process for decreasing the number of nodes of the PSSM, in order to meet real-time execution requirements, was described in detail. Finally, a series of simulation tests were performed, and demonstrated that the SbPSMS method can successfully achieve collision avoidance of manipulators whilst meeting the real-time requirements. However, different methods of setting up the PSSM may lead to different path geometries. Moreover, the method itself dictates the influence of step length on the smoothness of the path. Whilst decreasing the step length resulted in a smoother path for the coordinate system based method, minimal effects were seen when using the current velocity based method. Furthermore, the operation time of the current velocity based method is longer than that of the coordinate system based method. Therefore, it can be concluded that the coordinate system based method is more effective and demonstrated greater stability. The simulation of the coordinate system based method with triple manipulator cooperation produced satisfactory results, suggesting our algorithm can be extended to applications using multiple manipulators.
In future work, the maximum number of cooperating manipulators that the SbPSMS can deal with simultaneously should be clarified. In addition, planning of paths within entire SFR will be a major focus.

Author contribution
Chang Su conceived of the study, designed the study, and wrote the manuscript. All authors were involved in collecting and analyzing the data, also revisions