Abstract
1. Introduction
Symmetrically distributed limbs are a trait of higher primates. This fact can be observed as an optimum achieved by a long evolutionary process, or as a result of an optimization process. The limbs enable humans, for example, to complete tasks, which require the use of both hands and cannot be completed by using only one hand. In this line of thought, placing two robots next to each other enables the robots to be more flexible, to help each other and mutually increase their scope of work. Moreover, by placing robots so that their workspaces overlap, their DOFs are shared. Thus, the flexibility of the system increases. A standard industrial robot has a six DOF structure, whereas a human hand has 27 degrees of freedom in total, with 8 of them located in the wrist. The disparity is obvious and the consequence is that today robots, even the most sophisticated ones, have problems with opening a door, for example.
One could ask why a robot or a robotic hand with more than six degrees of freedom has not been built. It has, but the problems associated with controlling such a system are complex and limit its application. Problems with accuracy and repeatability limit implementation even more seriously [1, 2]. Building robots with a higher number of degrees of freedom is not the only way to increase their work scope. We adopt an anthropomatic approach, where a system consists of more than one robot, but the kinematics of the robots is relatively simple. This way, the control of each robot is simpler, but favourable traits of high flexibility, increased scope of work and autonomy result from the combination of several simple structures.
This paper deals with the problem of motion coordination for two standard six DOF industrial robots whose workspaces are shared (Figure 1).

Robot setup. Spheres indicate the corresponding workspaces.
The idea behind arranging robots this way is to increase the work scope of the robotic system with respect to the set of actions that can be performed by two separate robots. It is immediately possible to conclude that two robots in such a setup cannot work completely independently. The occurrence of collision is almost certain if one robot does not take into account the motion of the other one. The simplest way to tackle this problem is to ensure a safety envelope for each robot and to ensure that the other robot never enters the safety envelope of the first one. This is actually a very common way to solve this problem in industrial applications but limits the performance of the robots as a system working together and achieving synergy. Another common approach is to give a role of “master” to one of the robots and declare the other robot a “slave”. The slave robot constantly tracks the position of the master and adapts its movement so that no interference between the robots' links occurs. The problem is that such a scenario results in intermittent motion of the robots and time delays. There are also approaches that utilize simple breakaway screws or external shock sensors. In such applications, collision actually occurs, with all the physical consequences, replacement and restart procedures and the results in a production delay.
The issues mentioned above have resulted in an increased focus by the research community on finding suitable solutions for autonomous motion planning for a single robot, or in recent times, in the multi-robot domain. A comprehensive overview on robot motion planning can be found in [3]. If a path planning problem is presented in the optimization context, robust optimization techniques, such as evolutionary algorithms [4], swarm intelligence concepts [5] etc. have proven suitable, even though the problem is NP-complete and PSPACE-hard even in its simplest form. Its complexity increases exponentially with the dimensions of the robot's configuration space [6]. In [7] a special genetic algorithm (GA) for optimized robot trajectories is proposed. The main characteristics of this algorithm are the use of dynamic chromosomes structures and a modified crossover operator called an analogous crossover. The goal of the proposed GA is to minimize the accumulative deviation between the actual and the desired path. A multi-objective GA for the evolution of joint-space strings for manipulator configurations is proposed in [8]. Five indices, namely the travelling distance of the manipulator joints, the joint velocity, the Cartesian distance, the Cartesian velocity and energy are used to qualify the evolving trajectories. In [9] the preliminary results of the Constructive Solid Geometry-based approach to path planning of multiple robot arms are presented. The authors use a two phase GA to obtain a plan for robotic arms by using a strategy that combines the exploration of free space, while looking for the target position from each previously explored area.
The majority of papers dealing with the application of evolutionary algorithm-based methods for path planning consider significant simplifications, such as a single robot in an environment without obstacles or a robot surrounded by a set of stationary obstacles.
An approach to the path planning of a dual arm reconfigurable robot is presented in [10]. The paper deals with a SCARA configuration dual arm robot controlled by a single controller. The algorithm is based on the configuration space exploration for the given initial and final configuration of two hands, carried out in a horizontal plane.
A coupled dual arm system is a system composed of two robots with a point of contact between them. The contact is usually realized between the end-effectors of the two robots. In this scenario, impedance control-based methods may be used to solve the coordination of the two robots [11].
This paper presents a co-evolutionary approach to solving the path planning of two robots that share their workspace and are controlled by separate controllers. The input variables are the current positions of the two robots, given by the end-effector coordinates and the desired positions of the two robots, again given by the reference point coordinates in a reference coordinate system. The main characteristics of the proposed algorithm are: real value encoded chromosomes with joint angles on the loci, two coevolving populations (each population represents a set of potential solutions to path planning for one robot), the Pareto-based selection used for the evaluation of a multi objective inter-population fitness function and the hall of fame procedure for the preservation of the cross-population best collaborators (since the quality of each individual from one population depends on individuals from the other population).
The performance of the algorithm is tested through a number of simulation experiments and compared to the generic co-evolutionary algorithm without the Pareto-based selection and the hall of fame preservation. After the fine-tuning of the parameter space, the algorithm is implemented in a real robotic system, at this point limited to planning in a vertical or horizontal plane.
2. The proposed co-evolutionary algorithm
The purpose of this paper is not to outperform the existing heuristic approaches or to find better solutions to the problem of dual arm robot motion planning. Its purpose is to show how co-evolutionary algorithms may be efficiently applied to this problem and to describe how to do it.
Co-evolutionary algorithms offer great potential for concurrent multiagent domains [12–14]. Concepts of co-evolution were also successfully coupled with other heuristic optimization techniques [15–17] to solve complex optimization problems when search-spaces are connected. The main problems reported concerning co-evolutionary algorithms are their game-theoretic background and resulting pathologies, namely cyclic dynamics, loss of fitness gradient and evolutionary forgetting [18].
The motivation for using co-evolutionary algorithms to solve the defined problem is the distributed nature of the problem itself. Prior to this research, we experimented with standard evolutionary algorithms and realized [19] that the explicit fitness function for the evaluation of a system composed of a single robot with several optimization criteria becomes increasingly complex. The consequence is that the time required for finding feasible solutions is unacceptably long if the algorithm is to be implemented with real robots.
Co-evolution is by nature a distributed process, which can be applied to two separate controllers, making the computation parallel and thus reducing the time required for finding feasible solutions. Moreover, the definition of co-evolution is the adaptation of one species triggered by a change in the other species, which is directly applicable to the adaptation of the behaviour of one robot to the behaviour of another one. The abstract pseudo-code of a co-evolutionary algorithm is presented in Figure 2. The interaction between the two coevolving populations naturally mimics the solution of the problem described in this paper.

Pseudo-code for a sequential co-evolutionary algorithm
The setup used for the purpose of experimental evaluation in this paper is given in Figure 3. The two robots are arranged so as to ensure the overlapping of corresponding workspaces.

System architecture
Motion planning is based on two coevolving populations. Each population consists of a set of potential solutions to the problem of path planning from the given initial to the given final configuration for the left and right robot, respectively.
An individual in the population encodes a set of transitional configurations of the robot. At this time the planning is executed in the vertical or horizontal plane, limiting the search space only to two DOF for each robot. An individual in a population is represented as a chromosome containing a real-valued vector in the joint space of the robot (Equation 1):
where
Adopting the form of an individual in the population according to (1), which is a discrete form of an individual resulting in discrete information given to the robot controller, there is a trade off between the granularity of the discretization and the real time computational requirements. There is an important question - how can we decide what is the appropriate length of the chromosome representing the individual and defining the granularity of discretization?
It is assumed for this calculation that the robot moves as shown in Figure 4: if the second joint of the robot is static, the joint displacement of the first joint is Δφ1; if the first joint of the robot is static, the joint displacement of the second joint is Δφ2. In both cases the distance travelled by the reference point equals Δ

Calculation of the chromosome length
The following formula can be derived: Δφ1·(
where
where
Since the evaluation of populations composed of chromosomes whose length is 32 is computationally complex and therefore is impossible to perform in real time, we keep the length at
3. Fitness evaluation
It is obvious that fitness evaluation has to take into account the number of collisions between the links of the robots and the trajectory length. Both of these criteria have to be minimized. After the initial simulation it was determined that these two parameters are not enough and that additional criteria are required to obtain satisfactory behaviour of the robots, so velocity profile and total joint rotation angle are added. By adding these two parameters smooth motions of the robots, free of unnecessary joint rotations were obtained. The criteria are described in details as follows:
Collision penalty
The collision penalty, depending on a collision between Robot 1 and Robot 2 in corresponding configurations, is given by:
where:
Since the initial and final configurations are known in advance and by definition, a collision cannot exist, to save calculation time, these two configurations are omitted from the search, thus
The total distance of the end-effector movement is as follows:
where
Total angle of rotation
Since the robots are redundant systems even in this simple form of 2
where αi, is the angle between the link of the robot and the positive horizontal axis in the discrete time step is
End-effector velocity distribution
To enable an even distribution when going through points
where
Regarding other evolutionary operators, a simple single point crossover showed good results. Standard values of crossover in the range
The mutation operator replaces one allele with a given probability using Eq. 9.
The size of the populations is also an important part of successful convergence. It is said that one part of the problem of any evolutionary algorithm is creating the algorithm itself [20]. The other part of the problem is fine-tuning the algorithms parameters, whose interconnections are nonlinear, so often a significant effort is necessary to make the algorithm work well. This surely holds also for the co-evolutionary algorithms that have even more complex natures compared to standard evolutionary algorithms [21]. The problem with large populations is the time required to find satisfactory solutions. Note that in this problem, the best possible solution is not always searched for, since time is a critical factor. Rather any feasible solution that solves the problem in terms of collision free trajectory is accepted if there is not enough time to continue searching. Therefore, the sizes of the populations are 100 individuals. To increase the speed of the convergence, we modify the co-evolution to approach the real-time usability of the system, developing methods of selecting the best collaborators based on a combination of the Pareto front exploration and the hall of fame creation to save the best collaborators from the co-evolving populations.
The Pareto front exploration aims to identify the best individuals from a given population. The individuals forming the Pareto front are also known as non-dominated solutions. This front is a 4 -dimensional one, partial fronts are illustrated in Figure 5. It is a very significant and positive trait of the proposed algorithm that it is able to identify the Pareto fronts.
We want to employ this trait to keep the best solutions that form the Pareto front and allow these solutions unconditioned transition to the following generation.
These solutions might serve as a good seed to direct the evolution in the desired direction in a shorter period of time.
This is analogous to the elitism in a simple evolutionary algorithm, but different in the way that the elitism is not directly employable in co-evolutionary algorithms since the performance of each individual from one population depends on each individual from the other, co-evolving population, as already described.
It is obvious that the proposed algorithm can find solutions on the Pareto front in an early phase of evolution, after only 30 generations (Figure 5). This is an advantage of the algorithm, with the consequence of an increased likelihood of finding the global-optimal area of the multi objective search space. The red dots indicate individuals on the Pareto front, or non-dominated solutions, whereas blue dots indicate inferior, dominated solutions. Red, non-dominated solutions are the highest quality ones, i.e., the ones that are kept in the hall of fame until replaced with fitter ones.

Pareto fronts in an early phase of evolution (
The hall of fame is an operator that determines how many of the individuals are to be kept and directly copied to the following generation. The depth of the hall of fame is a percentage of the individuals in respect to the size of the original population that are directly copied to the next generation. The members of the hall of fame are only replaced when better (fitter) solutions on the Pareto front are identified.
The problem analysed in this paper represents a multi objective optimization problem with conflicting criteria. All the criteria are to be minimized, but they are conflicting in nature.
The time required to find a feasible solution becomes a critical parameter once the algorithm is implemented in the real physical system. In this paper, the convergence time is accepted if it is in the range of <10s. However, since the algorithm is heuristic, it is possible that no feasible solution is found in the predefined time range (<10s). The good thing is that it is easy to check whether the algorithm has converged, based on the critical criterion of a collision. If no convergence has occurred, and collision is present, the algorithm generates new populations and the search starts again. Since planning is done only once, at the beginning of the motion, a pause of 10
Regarding the number of generations required for convergence, we have compared several variations of the standard approach based on the roulette wheel selection, methods I-IV, Fig. 6. To increase the speed of convergence, we introduced the Pareto-based selection. The idea is to keep not only the best collaborators, as pairs of values from both populations in the hall of fame, but also to include the whole Pareto front of the current generation for evaluation. The methods based on the Pareto domination are the labelled methods V-VIII in Fig. 6. It is visible that a significant impact on the speed of convergence is achieved by introducing the new selection method.

Impact of the modified co-evolution on the convergence speed
The differences between methods are as follows: I – IV are individuals of the standard roulette wheel selection process, with different values of evolutionary parameters, such as size of populations, mutation and crossover probabilities.
Methods V–VIII are based on the Pareto selection, with different parameters of population sizes, crossover and mutation operator values. The methods are evaluated with 20 runs for each method and the corresponding results are illustrated, with standard deviation, in the upper part of Fig. 6 and the box-plot in the lower part of the same figure. We experimented with the depth of the hall of fame procedure. Namely, if the hall of fame is not changed during the process of evolution, an opposite effect in terms of convergence speed can occur, since the number of good individuals tends to grow. A consequence is that more calculations per generation are required. In an early phase of evolution, this is not important since the hall of fame is empty, or contains a small number of solutions. The parameter depth of the hall of fame indicates the memory effect, which determines after how many generations the new members will start to replace the old ones. For our purposes, we obtained the best results when a hard limit of 15% of the population size was set to limit the size of the hall of fame. It is important to note that the size of the front changes during evolution, since the global optimum for some criteria might be found, which means the front is actually reduced to a single point.
Regarding the mutation, as mentioned already, we experimented with a simple, constant mutation probability. It was determined that problems with a constant mutation rate occur in the later phase of evolution, when higher mutation rates have a devastating impact on the individuals. To overcome this issue the following procedure was designed and implemented. Eq. 10: higher mutation rates are beneficiary in an early phase of a search in order to quickly sample the search space. Once “good” parts of the search space are identified, the mutation should decrease in order to not destroy good individuals.
where
The mutation rate is variable, according to Eq. 10, when roulette wheel selection is implemented and individuals are encoded as real-valued vectors with joint angles at the loci. Since the number of evaluations per generation is large, it is important to identify the operation that is most-time consuming. It was identified that the collision check procedure is the most time-consuming operation. At first, we used a simple method of calculating intersections between robot links defined as lines. For a population of size 100 and two robots with two links each, the number of calculations for only these criteria is 100×100×2×2 and the same number of calculations has to be carried out for each discrete time step. This is a large number with a high impact on calculation time. We have analysed the problem and were able to significantly reduce this time by introducing dimensionless parameters and by transposing the calculation to only the angle of the rotation space. The procedure is as follows: Dimensionless parameters μ and
The dimensionless parameters are defined as

Analysis of intersection
Solving the following equation:
we can obtain the values for parameters μ and
The consequence of this procedure is a significantly les complex calculation than the calculation of th intersections of a set of lines, with the beneficia consequence of increasing the speed of the calculaton Analogously, the collisions between two arbitrary links o the robots can be checked.
4. Implementation and experimental results
The proposed algorithm was examined through multipl simulations in Matlab (Fig. 8 and 9) and finally on the

Motion of two robots in the Matlab simulation environment
This is the reason why a blank space is visible between the set of points describing consecutive positions of the end-effector positions of each robot. The picture itself describes the trail of the robot, starting from the initial configuration and ending with the final one.
Fig. 9 describes the interpolated values for the discrete joint rotation values obtained as a result of the co-evolutionary algorithm. The interpolation is performed using a polynomial of the third order to ensure second order derivability and thus smooth transitions. Also, it is necessary to check the gradient of the interpolating polynomial to physically enable the robot to perform the motion.

Interpolated joint rotation values for discrete values obtained by the co-evolutionary algorithm
The PC that runs Matlab performs all the calculations upon receiving the actual and desired configurations of the robot. To enable communication protocols in the system, we assigned the following IP addresses: Robot 1 (left robot): 192.168.123.26; Robot 2: 192.168.123.25 and PC: 192.168.123.40. Two communication channels were opened. The first channel is for sending data from a robot to the PC. In that case, the robot has the role of the server, while the PC is a client receiving the desired data. The second channel is for sending the response from the PC to the robot controller.
Now the PC has the role of the server and the robot controller is the client receiving the desired data. Messages that are sent are coded and can have different forms. They indicate whether the motion is possible or not (i.e., whether the intersection of the robot links in the initial or final configurations). If the motion is possible, a list of coordinates and going-through positions of the end-effector is sent.
It is important to ensure that the time delay between the start of the motion of both robots is kept to minimum. This is done through directly connected I/O signals between the two robot controllers.
The role of these signals is to trigger the motion at a desired moment in both robots. The algorithm was tested for various scenarios: different initial and final conditions, one robot moving while the other one is waiting, stationary obstacles present in the workspace, etc.
Fig. 10 and 11 present the final checks prior to implementation in the real robot setup. The pictures are from the Roboguide physical simulator developed by Fanuc Robotics. Fig. 10 is for two Lr Mate robots, for a simple scenario when one robot is stationary. In this case, the problem boils down to avoiding a stationary obstacle.

Motion execution in the Roboguide physical simulator

Safety factor for the end-effector
This is a simple scenario, but a problem was noticed after implementation in the physical simulator. Namely, since the robots are defined as a set of lines, the collision can actually occur due to the physical dimensions of the real robot, even if the solution from the simulator is collision-free. To tackle this problem, we described each link of the robot as a set of two boundary lines.
A circle with a variable radius size parameter is defined around the TCP to enable tuning according to the current tool the robot is carrying (as illustrated in Fig. 11). Since the robots are to have various tools attached to the flange, we immediately noticed the problem of the fixed geometry of the robot as implemented in Matlab.
The three trajectories in Fig. 11, implemented in the M10 robot, in this case illustrate how the safety factor works, forcing the robot to move away from a potentially dangerous zone by various amounts depending on the robot's current setup. In order to implement the algorithm in our M10 robot, we had to adapt the lengths of the links and global coordinate systems in the Matlab simulation environment. Note also that this motion was executed in the horizontal plane, providing robots with only the desired position of the end-effector and using robots internal kinematics model to calculate everything else.
The distance between the lines describing the robot and the sides of the rectangle were experimentally determined to fit the configurations of the robots used in our laboratory. This way the problems of inconsistencies between the model and the real robot were eliminated.
5. Conclusions
In industrial environments where several autonomous robots work in a shared environment, it is essential to find feasible, collision-free motions.
In this paper, we have presented a method based on co-evolutionary algorithms that finds optimal or near–optimal trajectories for two industrial robots that share a workspace. Exact values for the granularity of calculation, through chromosome length, are calculated, based on the physical limitations of the robotic system used for experimental verification. This approach is scalable in terms of adapting it to different dimensions of robots. Also, it is general and implies the possibility of adding new degrees of freedom or new robots to the system. Standard approaches to co-evolutionary algorithms are implemented and evaluated and modifications, including the Pareto-based selection and the transition of collaborations to the new generations, are presented. The methods are compared and experimentally verified. The result of the modifications is the increase in speed of the convergence in the case of the modified algorithm. A formal model analysis of the dual arm robot is performed to simplify the time-consuming operation of collision detection. The computational complexity is reduced by means of a collision check conducted only by corresponding angles of rotation, instead of using the exhaustive line intersections procedure.
The algorithm developed in this research is expected to serve as a low–level procedure ensuring the safe motion of robots requested to perform high–level tasks. To solve the problem, we have developed a simulation environment that analyses movements of simplified robot models. The next step is implementation in the physical simulation environment, where the interaction between physical constraints and the simplified model occurs and the tuning of the model can be performed to satisfy real-world constraints. Finally, the full implementation is performed on two Fanuc Lr Mate 200
We anticipate further work in adding more robots and more degrees of freedom to the model. The model is developed to be scalable and allows the inclusion of additional features, although important questions regarding the computational power are expected. The main bottleneck of the approach remains its computational complexity, despite the effort to include algorithmic simplification.
A possible way to overcome this problem is to apply the co-evolutionary process to the distributed computer architecture and to make independent evolutionary operators perform in parallel.
At this point in the research, planning is performed only once, at the beginning of the motion of the two robots. This is fine if nothing unexpected occurs during the motion of the two robots. If one robot unexpectedly stops during the motion, there would be a possibility of a collision with the other robot, since the robots are not aware of each other after the motion starts. This is a limitation, but future work will consider permanent motion sampling for the two robots. If a possibility of collision is detected, the planning routine can be invoked more than once to avoid collisions after initial planning. Features for web-based, real time monitoring and control of such a complex robotic setup, see for example [22, 23] are also considered for future research. This would further increase the safety of a system tending to be of a high level of autonomy.
