Multi-robot cooperative 3D printing (C3DP) is a growing area of research. The motivation behind C3DP is to achieve efficiency, flexibility, and scalability in additive manufacturing by leveraging simultaneous deposition of material at multiple locations for a given part [1]. Despite the promise of C3DP, the so-called cooperation between different robots is inherently challenging due to several interlinked and often conflicting objectives such as collision avoidance, task scheduling, part segmentation, and so on. Accordingly, there have been several investigations of each of these problems individually or in combination with each other.
A critical gap in current research on C3DP is the lack of a principled approach that encapsulates the myriad challenges within a common formal framework. To address this gap, we present a topologicallyguided framework for spatial reasoning in C3DP. Our framework focuses on layer-wise cooperation across an arbitrary number of robots tasked with printing arbitrarily shaped parts. The core idea behind this framework is the notion of spatial parameterization induced by the partitions of a given layer assigned to each robotic printer. When applied to layers with internal cellular structures (as is typically the Fig. 1. The traveling Cellsman algorithm enables cooperation between robots through parametrically clustering the design space and allowing for optimization of a robots schedule. could vary ranging from maximizing cooperation between printers to the minimization of makespan under essential constraints such as collision avoidance. In many cases, the segmentation itself imposes restrictions or constraints on the printing schedule such as task dependency of sub-volumes [8]. Other layer-wise collaborative efforts create ''safezone-style'' regions where only one robot may print at a time [7,9,10].
We observe that the common element that connects these issues is that of determining the right spatial parameterization that simultaneously captures division of labor and scheduling and enables a good formalization of desirable objectives (e.g., makespan) under given constraints (e.g., collision avoidance). In our work, we specifically consider the scenario where multiple stationary robots cooperatively print in a layer-wise manner, and we assume a global-local control and coordination that enable robots to pause upon the detection of proximity with other robots. With these assumptions, our technical problem is to determine such a parameterization and also the methodology to subsequently use the parameterization to minimize the makespan (objective) while avoiding collision (constraints).
The main focus of this work is to develop a methodology to determine the sequence in which the cells in a given layer are printed. The design space of sequences, as such, is combinatorially large. We observe that our problem shares a remarkable resemblance with a variation of the TSP called the multiple Traveling Salesman Problem (MTSP) [15,16] (Fig. 2(a)). Unlike standard MTSP, where tasks can be freely reassigned to salesmen (robots in our case), our problem is more closely related to a colored MTSP approach where cells are preassigned to robots due to kinematic and positioning constraints [11,12,17]. For our problem, the coloring of cells is dependent on the partition of the printing space, only robots can print cells within their partition (Fig. 2(b)). To manage the tradeoff between combinatorial complexity and objective performance, additional works have further superimposed the notion of clustering that, in our case, would correspond to grouping cells in order to reduce the size of the design space when sequencing the tasks (Fig. 2(c)) [13,14]. Our approach in this work stems from the question: ''How should we cluster cells so as to determine the optimal sequence in which cells should be printed without compromising on the objective performance?' ' We begin by observing that the interfacing boundaries between subvolumes assigned to their respective printers naturally lend themselves to a parameterization that uniquely characterizes distances from those Fig. 2. Different variations of a TSP problem where tasks are not preassigned (a), preassigned to specific robots (colored) (b) [11,12], clustered into groups and preassigned to specific robots (colored-clustered) (c) [13,14], or clustered using a spatial parameterization (colored-parametrically clustered) (d).
boundaries. Formally, the level-sets of the distance field from the boundaries offer a direct way to determine the gross movement of each printing nozzle. Note that for layers subdivided into cells, this parameterization essentially translates to task scheduling by virtue of discretization. With this in view, our broader approach develops an interface-induced parameterization which can then be used as a representation to compute optimal sequence of cells. The crux of our approach is to use this parameterization as the means to generate clusters for colored clustered MTSP. We call this approach parametric clustering (Fig. 2(d)).
The first and fundamental contribution of this work is the development of interface-induced parameterization and clustering methodology for scheduling of tasks in layer-wise cooperative 3D printing (C3DP). This theoretically grounded methodology offers a natural representation for task scheduling for minimization of total printing time with C3DP. To this effect, our second contribution is the use of parametric clustering for the development of a novel computational framework for minimizing the total time to print. This framework, inspired by clustered multi-color traveling salesman problem is implemented using genetic algorithms (GA). Note that parametric clustering shares its strategy with min-max clustered MTSP [14], which is also formulated within the GA framework. However, the problem considered by Bao et al. [14] is conceived for topological graphs and does not offer any particular method for obtaining a clustering. In our case, the problem is situated in a metric space (specifically the Euclidean plane) and the method for clustering is therefore necessary. Not only that, our ''graph'' is induced from a tessellation of the underlying metric space. The key gap that our framework fills, therefore, is that it elegantly leverages our idea of parametric clustering to encode the printing sequence of cells in terms of two types of clusters, namely, radial and medial clusters of cells. Furthermore, we make a secondary contribution by developing a simplified area-based function for evaluating the print time to offer significant computational efficiency during GA-based optimization.
In terms of insights, our comprehensive numerical comparisons show that our parametric clustering arrives at a better optimal schedule much faster than a non-clustered sequencing approach in GA-based optimization. We also show that by changing the clustering type from medial to radial (or vise-versa), the time to print can be reduced. Through physical evaluation with three robots we found that an optimized strategy allows for more parallel printing than a non-optimized and a slicer-based approach. For example, for a long part (approximately 850 mm × 270 mm) cooperatively printed using 3 robots, our method helps save 8.5 min on average per layer leading to hours of time saved across the entire volume.
There has been extensive prior research on systems which have multiple robots which must collaborate in some way. The amount of collaboration between robots varies based on the system While it is also common to Some areas where collaborative systems are used are coverage tasks [18,19], assembly/disassembly in industrial robotics [20][21][22], and painting/spraying tasks [23,24].
While this paper is generally focused on multi-robot systems, we use collaborative 3D printing(C3DP) to validate where multiple robots must work together to share tasks and avoid collision. There are various different C3DP application systems with different robot geometries and setups with some using gantry systems, moving robots, robotic arms, welding wire arc robots, or flying robots [5,[25][26][27][28][29][30][31]. In these cases robots print in either planer [1] or non-planar [32] manner with multiple robots with various degrees of freedom. Robots printing tasks are within close proximity since they are working on a common part which can make collision avoidance a difficult challenge.
In many cases when robots work in close environments it is preferable for the robots to work as independently as possible [18,33,34]. These have been compared before to CPU systems where most/all tasks are parallelized without much difficulty [35]. Other methods exist for systems with temporal constraints [36], specific manufacturing constraints [37], or with cluttered environments [38].
Within C3DP systems tasks have been scheduled using various methods such as scheduling paths by creating subpaths [25,39,40], or by scheduling sub-volumes [1,41,42]. Other methods involve labeling specific regions as safe/unsafe and then sequentially printing the unsafe parts followed by collaborative printing of safe regions [6,7], others try to schedule the printing of unsafe and safe regions without pausing [9,10]. These are either limited by having large sub-regions which can be difficult to schedule well, or a low amount of collaborative printing.
One of the main limitations of prior work is that their predefined goal is to keep the number of regions low, in order to decrease the solution space, and therefore can solve most, if not all, potential schedules and determine the optimal solution [10]. In these cases, an increase in the number of regions (which is the case with cellular decompositions, as is common with infills) to print significantly increases the difficulty of determining the optimal schedule. For cellularized layers, prior works have adopted restrictive strategies to avoid collisions through efficient pausing/waiting [6,7]. Strategies where ordering of work is done independently to avoid collision by constantly checking robot position lead to high communication overhead and non-optimal printing [43]. In our work, we present an approach that is principled and generalizable, works for highly cellular layer geometry, and does so while avoiding the overhead caused by constant proximity checks between robots.
The traveling salesman problem (TSP) has been a popular problem in computer science [15] with applications in robot task scheduling, where tasks are treated as points and distance is minimized. With variations such as multiple TSP (MTSP) [16], where some have applied it to multi-robot task allocation and scheduling problems for optimization (Fig. 2(a)) [44][45][46].
However, in many MTSP problems with multiple robots, collision is not a consideration, and most problems can be reduced to a series of TSPs solved separately [13,35]. While this may be useful, it is often not the case that robots can work separately since they may share workspaces. One variation of the MTSP called the colored MTSP aims to reduce the possible solution space by assigning some, if not all, cities to certain salesmen (Fig. 2(b)) [11][12][13]47]. This is particularly applicable to many robotic task allocation jobs where only some robots can accomplish a task, due to manipulator requirements, tool availability, or reachability of the robots [12,48].
Since TSP and MTSP are combinatorial problems, various methods exists to obtain optimal or near-optimal solutions. Most approaches uses a metaheuristic optimization strategy such as simulated annealing, particle swarm optimization, or genetic algorithms [13,44,49,50]. In many cases, the design of the chromosome for genetic algorithms allows for the swapping of tasks between robots without increasing the design space [47,[50][51][52].
Another method for solving TSP is to cluster tasks into groups and then order the clusters, since the number of clusters is less than the number of tasks, the design space is reduced (Fig. 2(c)) [14,[53][54][55]. A common method for clustering is to use k-means clustering, which creates groups of points that are close to one another [56]. Another use of k-means clustering is to assign tasks in MTSP where the number of clusters matches the number of robots. After the assignment, the order for each robot is optimized independently [13]. While this is effective in grouping tasks, it groups tasks that are close to each other, which may not allow for collision avoidance.
There has also been work on solving TSP problems for 3D printing. One work uses TSP to order printing paths to minimize the number of retractions during printing which can take a long time, especially in DED or concrete printing [57,58]. These works place points within the object being printed and then use these points to traverse and optimize printing while minimizing the number of retractions needed. Other works use given paths and optimize the traversal between paths [59].
While the TSP and other scheduling methods have been applied before to collaborative robot systems, there currently exists a gap in literature for when there are a high number of tasks that must be scheduled without collision. When the number of tasks is reasonable, the number of combinations is reasonable to solve. In the case of 20 tasks, there are 2.43 × 10 18 combinations for a single robot or 7.25 × 10 6 combinations for two robots with equal distribution of tasks. However, if there are 100 tasks, then the problem is likely intractable. With single robot TSP, there are 9.33 × 10 157 combinations, and with two robot preassigned MTSP, there are 6.08 × 10 64 combinations. This is the gap that our work fills by parametrically clustering space using the coloring of the tasks so that an optimal solution can be found quickly (Fig. 2(d)). The closest terminology within the TSP space for this problem would be a Colored Parametrically-Clustered MTSP.
General TSP problems have a solution space of 𝑁! where 𝑁 is the number of tasks [13,47]. This can be reduced if the space is symmetric, meaning that traversal from point 1 to point 2 is identical to traversal from point 2 to point 1. This is not the case when a collision must be considered. Therefore, the solution space cannot be reduced [13]. For MTSP, where tasks are not pre-assigned, then the number of combinations is still at least 𝑁! depending on chromosome representation [11,52]. If tasks are preassigned to 𝑚 robots then the solution space is
where |𝑁 𝑘 | is the number of tasks assigned to the 𝑘th robot [47]. As mentioned before one method to reduce the solution space is to cluster tasks which reduces 𝑁 𝑘 if there are more than 1 task per cluster. This can be very useful when there are a high number of tasks. As an example, if there are 300 tasks, the solution space for TSP is 3.06 × 10 614 and with 3 robots and equal distribution of tasks the solution space becomes 8.128 × 10 473 . If each robot has 10 clusters, then the solution space becomes a more reasonable 4.778 × 10 19 possible combinations.
Our conceptual framework starts with a few assumptions about the robot-task space. The first assumption is that all tasks are located on surfaces (2-manifolds). This assumption is consistent with a layer-wise printing approach, usually seen in additive manufacturing applications. While we only show results for flat horizontal surfaces, it is possible for the same algorithm to work on curved orientable surfaces. The second assumption is that the tasks are located within some finite region. The robots are assumed to be connected with a home position fixed in space with kinematic constraints, in our case a two-linkage mechanism connects the home, and nozzle position.
We first define the mathematical framework for parametric clustering, which uses level sets of the partition to form medial or radial clustering. Then we discretize since our partitioning consists of cells to be printed and use the same approach in order to generate radial or medial clusters of space.
Consider the geometry 𝛺 of a layer and a number of robots with fixed locations. We assume that 𝛺 is a set of multiply connected components with potential holes (i.e. annular regions). A natural way to assign sub-regions of 𝛺 to the robots is to partition it into sub-regions where each sub-region is closest to its corresponding robot. Formally, this can be achieved by considering a Voronoi tessellation ( ) of the underlying space (R 2 in our case) where = { i } is the set of mutually exclusive and exhaustive faces/cells of . Further, = {𝝂 𝒊 }, and = {𝜀 𝑖 } denote the vertices and edges of , respectively. In this setup, the task of printing a sub-region to robot 𝑖 is equivalent to determining the intersection 𝛺 ∩ i . Fig. 3. The space can be parameterized using the partition (a) either radial (b) or medial (c) based on the distance from the partition vertices or edges respectively. This is the same in the discrete space with cells (bottom).
What is important to note is that the edges and vertices of this Voronoi tessellation define the boundaries of i where any edge 𝜀 ∈ is equidistant to two robots and any vertex 𝝂 ∈ is equidistant to three or more robots. As a result, and define the interface across the robots. This interface offers two natural ways for defining parameterization of the entire plane, namely medial and radial, based on the Euclidean distance function with respect to elements of and respectively.
Medial parameterization, in the continuous case, is characterized by the level sets of the function 𝑝(𝐱) = 𝑚𝑖𝑛 𝜀∈ 𝑑(𝐱, 𝜀) where 𝐱 ∈ R 2 is a point in the plane and 𝑑(𝐱, 𝜀) is the point-to-line segment distance between the point 𝐱 and an edge 𝜀 ∈ (Fig. 3(a)). Now, given a tessellation of 𝛺 into cells, the level sets of 𝑝(𝐱) can be used to cluster the cells for each robot with respect to the edges 𝐸 𝑃 of the partitioning (Fig. 3(a)).
Radial parameterization, in the continuous case, is characterized by the level sets of the function 𝑞(𝐱) = 𝑚𝑖𝑛 𝝂∈ 𝑑(𝐱, 𝝂) where 𝐱 ∈ R 2 is a point in the plane and 𝑑(𝐱, 𝝂) is the Euclidean distance between the point 𝐱 and a vertex 𝝂 ∈ . Now, given a tessellation of 𝛺 into cells, the level sets of 𝑞(𝐱) can be used to cluster the cells for each robot with respect to the branching points of the partitioning (Fig. 3(b)).
As such, interface-induced parameterization can directly be used to define robot paths. All that needs to be done is take the appropriate level sets of either field. However, in most practical additive manufacturing scenarios [2,6,7], it is reasonable (and many a time desirable) Fig. 4. The total number of combinations for solving the TSP using Medial, Radial, or no clustering.
to decompose the domain 𝛺 into mutually exclusive and exhaustive cellular regions; in other words a tessellation of 𝛺. Formally, we begin by redefining 𝛺 as a tessellation 𝛺(𝑉 𝛺 , 𝐸 𝛺 , 𝐹 𝛺 ) with 𝑉 𝛺 , , 𝐸 𝛺 , and 𝐹 𝛺 as the vertices, edges, and faces/cells of 𝛺. Consequently, the assignment of sub-regions of 𝛺 to a given robot 𝑖 ∈ [1, 𝑛] is represented as the partition
Here, 𝐜 𝐟 is the position of the centroid of a cell 𝑓 ∈ 𝐹 𝛺 and 𝐫 𝐣 is the position of robot 𝑗. This effectively produces a discretization of the original Voronoi tessellation leading to a discrete version of medial and radial parameterization.
The first step to obtain the medial and radial parameterization in the discrete case is to define the discrete counterparts 𝑉 ⊂ 𝑉 𝛺 and 𝐸 ⊂ 𝐸 𝛺 for the interfaces and of . Let 𝐹 𝛺 (𝐯) be the faces incident on a vertex 𝐯 ∈ 𝑉 𝛺 . Further, let 𝑃 (𝑓 ) denote the partition to which a cell 𝑓 belongs. The interfacing vertices 𝑉 are defined as 𝑉 = {𝐯 ∈ 𝑉 𝛺 | 𝑃 (𝑓 𝑖 ) ≠ 𝑃 (𝑓 𝑗 ) ∀ 𝑓 𝑖 , 𝑓 𝑗 ∈ 𝐹 𝛺 (𝐯), 𝑖 ≠ 𝑗}. In simpler terms, 𝑉 is the set of vertices of 𝑉 𝛺 shared by faces, none of which belong to the same partition. Similarly, to obtain 𝐸, let us first denote 𝐹 𝛺 (𝑒) to be the faces incident on the edge 𝑒 ∈ 𝐸 𝛺 . Unlike vertices, note that all edges (other than the boundary edges) are shared by exactly two faces. Therefore, the interfacing edges are defined as 𝐸 = {𝑒 ∈ 𝐸 𝛺 | 𝑃 (𝑓 1 ) ≠ 𝑃 (𝑓 2 ) , 𝑓 1 , 𝑓 2 ∈ 𝐹 𝛺 (𝑒)}. In other words, the interfacing edges are the ones shared by faces in two different partitions (i.e., to be printed by two different robots).
The sets 𝑉 and 𝐸 combined form the interfacing complex [6] and allow for a completely topological encoding for clustering the cells of 𝐹 𝛺 based on medial and radial parameterization. The key idea is to apply the moving fronts strategy as shown previously by Ebert et al. [6] and Krishnamurthy et al. [7]. In other words, This is due to a combination of generating a tessellation of 𝛺 into cells (𝐹 𝛺 ) and subsequently labeling the cells of 𝐹 𝛺 with respect to distance from the closest robot.
The first cluster using medial parameterization is formed by the cells which neighbor the partition edges which can be defined as
Effectively, the cells in cluster 𝐶 𝑖,𝑘 are the immediate neighbors of cells in cluster 𝐶 𝑖,𝑘-1 and do not belong to any previous cluster (Fig. 3(c)).
The first radial cluster is a set of cells that shares neighborhood with the interfacing vertices 𝑉 . Formally, the first cluster is defined by 𝐶 𝑖,1 = {𝑓 ∈ 𝑃 𝑖 | 𝑉 𝛺 (𝑓 ) ∩ 𝑉 ≠ ∅}. Note that once the first radial cluster is determines, the method for obtaining the subsequent clusters is identical to that for obtaining the medial clusters. Formally, the radial clusters following the first radial cluster are defined by
While the formation of the clusters through medial or radial parameterization allows for a smaller design space, that is, fewer possible combinations when ordering the tasks, the clustering does not inherently avoid collisions. Although collision can be avoided in real-time, a poor ordering of the clusters can result in long print times, as the number of pauses increases. Thus, we optimize the overall print time by using a genetic algorithm (GA) to select the ordering to visit the clusters. This is similar to many other TSP problems where the sequencing of cities/tasks can be modified depending on their impact on the improvement of distance/time.
To represent the ordering of clusters within the GA framework, we start by defining the clusters for all robots, either in a medial or radial fashion. The number of clusters for each robot is labeled as 𝐺 𝑖 for robot 𝑖 which is a positive integer. The clusters for each robot can be labeled according to their cluster number 𝑐, their respective robot 𝑖, and the order in which they will be printed 𝑝, by adopting the notation 𝐶 𝑝 𝑖,𝑐 (Fig. 5). As an example, 𝐶 5 1,3 indicates that cluster 3 of robot 1 is the 5th cluster to be printed. Using this notation, the schedule for a robot can be represented by 𝑆 𝑖 = [𝐶 1 𝑖,1 , 𝐶 2 𝑖,2 , …] which denotes the initial print order, i.e. cluster one is printed first, cluster two second, and so on (Fig. 5(b) and (c)). Furthermore, the entire schedule 𝑆 can be represented by 𝑆 = [𝑆 1 , 𝑆 2 , … , 𝑆 𝑚 ] (Fig. 5(c)) if there are 𝑚 robots with a solution size of ∏ 𝑚 𝑛=1 𝑆 𝑛 !. As an example, if there are three robots with 54, 53, and 56 tasks each, the total number of combinations for a nonclustered TSP is 7.02 × 10 215 (Fig. 4(c)). Using radial or medial clustering reduces the design space to 5.89 × 10 14 and 3.46 × 10 5 respectively (Fig. 4(a) and (b)). This chromosome representation is very similar to the setup described by Carter et al. which is described as a ''two-part chromosome'' [50] or the ''Path-city number'' chromosome described by Jiang et al. [60].
While the order of printing clusters may change, it is important to follow the rule that no two clusters for a robot can have the same schedule number, that is, 𝐶 𝑗 𝑖,𝑐 ≠ 𝐶 𝑗 𝑖,𝑐 ′ ∀𝑐 ≠ 𝑐 ′ . It follows that the scheduling order for all clusters must be integers between 1 and 𝐺 𝑖 (Fig. 6). Fig. 5. In order to represent the order in which clusters of tasks are scheduled, we first number all clusters (a) and then create a schedule of when clusters will be visited (b). Using this ordering per robot we can create a chromosome that describes the other in which clusters will be scheduled. Fig. 6. The order which robots print their clusters. Some clusters may take longer to print than others depending on the amount of printing required, this is an illustrative example.
Typical crossovers in the GA framework involve selecting two parent chromosomes and swapping parts of their genes. This would not work in our case, since it would not ensure that each cluster is selected only once. Instead, in our case, crossovers affect only the schedule 𝑆 𝑖 of a single robot. We first select a random integer between 1 and ⌊(𝐺 𝑖 ∕2)⌋, then select two non-overlapping portions of the 𝑆 𝑖 with the size of the random integer. These two non-overlapping portions of 𝑆 𝑖 can then be swapped. Finally, each switched portion can also be randomly flipped to form the chromosome (Fig. 7). This crossover methodology maintains that all numbers between 1 and 𝐺 𝑖 exist within 𝑆 𝑛 , meaning that all clusters are visited exactly once, while also allowing for different schedules. This is similar to only allowing for the ''in-route mutation/crossover'' described in prior works [52]. While another common method for modifying a chromosome within the GA is a mutation, in our problem formulation this is equivalent to only swapping two indices of schedule 𝑆 𝑖 since a single index cannot be modified alone without violating the property of all clusters being visited exactly once (Fig. 7(a) 𝑆 3 ). From our initial experiments, including a separate mutation function did not show a significant difference in results.
Through changing the print schedule, it is possible to avoid most collisions, but it is also likely that there are no solutions that completely avoid collision due to part geometry, robot locations, or robot kinematics. For example, if a small part is being printed by many robots, the robots will all be working around the same region, which can lead Fig. 7. How crossover happens within different schedules by either direct swapping (𝑆 1 ), swapping as well as reversing ordering (𝑆 2 ), or swapping the schedule order of just two clusters (𝑆 3 ).
to unsafe printing conditions. Because of these cases, the robots are always in communication with a central hub, which can determine if two robots are about to collide and send a pause command to one of the robots. This is similar to prior work in the area [43]. While there is an upper limit on the number of times robot locations may be checked (frequency) based on hardware, we assume for this work that Fig. 8. A representation of how changes in the clustering method can occur within a chromosome and how the change affects sequencing. The original chromosome is taken from Fig. 7.
the current location of robots may be checked 10 times a second and once told to pause to avoid collisions, robots stop their movement and pause within 10 ms. Once a robot receives a pause command, it must wait for 30 s before trying to continue. The safe distance for collision is 3 times the nozzle radius and the priority regarding which robot should pause in a two-robot collision situation is determined a priori.
In the description of the chromosome, the clustering methodology was chosen a priori. However, it is possible in some scenarios that it is not possible to reasonably intuit which may be preferable (medial or radial), especially with a high number of robots. Thus, we modify the chromosome to allow for changing the type of clustering within the chromosome.
We start by pre-computing the radial and medial clusters and storing them. For each robot, we then determine 𝐺 𝑖 for each robot by checking how many medial and radial clusters there are and taking the maximum value between them. The chromosome is then represented by an additional 𝑚 values which represent medial or radial selection (''M''= medial, ''R'' = radial). As an example, a chromosome for 3 robots can be represented by 𝑆 = [𝑆 1 , 𝑆 2 , 𝑆 3 , 𝑅, 𝑅, 𝑀] which means that the clusters of the first two robots are radially parameterized, and the clusters of the last robot are medially arranged (Fig. 8). Some clusters may contain no cells if there are unequal numbers of radial and medial clusters. As an example, if there are 5 medial clusters, but 6 radial clusters for robot 1 then 𝐶 1,6 contains no cells when medial, but contains tasks for radial clustering (Fig. 8 b).
In this chromosome representation, crossovers can include flipping the last 𝑚 values of 𝑆 to change the clustering method. From our initial testing, it was determined that making this change was rarely preferable, i.e. allow for within-cluster optimization and only changing clustering method less often. This is likely due to needing several generations to optimize a given clustering method and changing often does not allow this to occur.
In order to assess the chromosome choice, that is, the clustering and its associated ordering, we have to compute the time it would take to print. While most TSP problems only require distance to compute the optimum, robot tasks must also take into consideration collision avoidance strategies and division of work. For this reason, we simulate the time to print a given sequence using two different methods that replicate the physical setup to different degrees of accuracy and time computation, namely the path-based and the area-based evaluations.
The path-based evaluation relies on using the sliced path to compute the total time to print, which can be cumbersome. The path for each robot is computed based on the cluster scheduling to determine the cell printing order. The path for each cell can be computed by simply computing offsets of the cell to be printed, either completely filled in or with a set number of layers. Between the printing of each cell, the robot performs movements. Once the path for all robots is computed, each robot has its base position and the home position where it can return without collision. The printing process then begins with the robot moving from its safe position to the starting point of its path. Once that point is reached, we assume that it moves along the printing path with a set constant speed. At set intervals (similar to the robot communication setup), the robots check collision with other robots, which is limited to a frequency of 10 times per second to match the physical setup. The more often collision checks are made, the longer to compute since computing square roots frequently is computationally expensive. The most straightforward collision avoidance strategy is to check the distance between the end-effectors of any two robots, and assign a threshold value as a minimum acceptable distance. The choice of which robot is made to pause when a collision is detected is determined by a hierarchy (robot 1 pauses before robot 2 pauses etc.). More complicated collision-avoidance strategies can be used, such as arm-to-arm distance of robots to further verify collision avoidance. If a distance is below the given threshold, then a robot moves to the home position for a minimum of 30 s and returns to its path. Using this strategy, we compute the time a robot spends extruding material, moving between cells, pausing to avoid collision, and waiting to move to the next layer.
This method for print time evaluation can be costly to evaluate since distance checks for every robot, including arm-to-arm distance checks can be computationally expensive when evaluated multiple times a second.
Contrary to path-based evaluation, an area-based time evaluation can also be computed since the printing time of a cell is known and the position of a robot while printing a cell is constrained within specific bounds (i.e., the face of the cell). This process starts by computing the polygon-to-polygon distance for all cells and taking the smallest value (i.e., the most conservative approach). The arm-to-arm distance can also be computed for every cell based on the robot during the printing of the given cells using the boundary of the cell. The threshold for a pair of cells is then set to either the minimum cell distance or the armarm distance, whichever is smaller. The time to print a cell is computed based on the path that it is printed and the speed of the robot. Finally, the distance from every cell to every other cell can be computed which gives the time for a robot to move from one cell to any other cell given the speed of the robot.
The printing simulation can then begin by robots starting their assigned cells and checking for collision using precomputed cell values. As robots finish printing their cells, they move to the next cells using the precomputed values with no need to check for collision multiple times a second. The transition between clusters also consists of just transitions between cells. If the distance between two printed cells is lower than a threshold, then one robot moves to its safe position and waits for 30 s to continue. Each robot follows this process until all work is done while also keeping track of the time it spent printing, moving, pausing, and idling.
Note that most computationally expensive steps can be done a priori and then used during the optimization process to speed up convergence. Although this evaluation method is not as precise as the path-based method since the exact robot positions are not known, it allows for fast evaluation of printing time. A detailed comparison of the accuracy between the two methods is presented in later sections.
As a first implementation of our parametric clustering, we implemented our GA-based makespan minimization algorithm and conducted a series of systematic experiments using the gene encoding described earlier. The main goal of the experiments is to investigate (a) different fitness functions, (b) area vs. path based time evaluation, (c) the type of parameterization (radial/medial) on the makespan. For systematic evaluation, restricted the shape of the layer to a circular region with a radius of 600 mm and a two-tile tessellation generated using a cross section of the 3D truncated octahedral tessellation (Fig. 4).
Since a GA is a metaheuristic method, we use the following hyperparameters (unless stated otherwise) for our set of numerical results. To select the maximum number of generations, population per generation, and stall generations (number of generations without change before convergence is assumed), we carried out a series of studies using the area-based time evaluation. We found that while increasing the population size generally lowers the optimal makespan, it increases the computation time significantly for a marginal improvement (see Supp. Tables 5, 6 and, 7). For this reason, we choose a population size of 50, and a stopping generation of 50 since it terminates at a good stopping point in a reasonable amount of time. We also set the maximum number of generations to 5000.
The population for the genetic algorithm is initialized to a random sequence of clusters and each robot is randomly assigned to either a medial or radial clustering with equal probability. Both the sequence and the clustering method can change throughout the optimization process.
There are two other parameters that are necessary to establish for pausing a given robot, namely, the duration for which a robot pauses before resuming activity and a distance threshold for registering a collision. Based on preliminary physical experiments and previous works [10] that use the same physical setup as ours, we use a pausing duration of 30 s along with a safe distance of 3 times the nozzle radius. In order to confirm the sensitivity of these parameters on the time to print, we conducted numerical experiments with different ranges of pause durations and safe distances. Our analysis showed that increasing the enforced pause time has little effect on the overall time to print. So long as the pause duration is kept below 50 s (see Supp. Table 8 and 9), there is no effect. A large safe distance can make it difficult to find a good solution, and lead to unnecessary pausing. A safe distance of 3 times the nozzle radius, and a 30-s pause time is very close to the global minimum found (see Supp. Table 9)
The proposed chromosome representation for the printing schedule can be used directly in the GA framework and the subsequent evaluation of the printing time, either through the path-or area-based evaluators. However, the time to print is not a viable fitness function, as a large portion of the time spent in the process is due to the extrusion of the material, which cannot be minimized. We first start by tracking the how the robots spend their time during printing. There are 4 total activities that the robots can spend their time.
1. Extruding -The time robots spend completing their extruding/depositing material 2. Moving -The time a robot spends moving from one cell to another 3. Pausing -The time robots spend pausing to avoid collision with other robots, this includes time moving from its home position back to printing 4. Idling -The time a robot spends after all extruding tasks have been completed waiting for the other robots to finish so the next layer can be started Since the extruding time is only determined by the geometry of the printing part, it cannot be minimized, and therefore is not part of the fitness function. Furthermore, the time a robot spends idling is a function of the distribution of work between robots (equal distribution of extrusion time) and the pause/move time. As an example, if one robot has high pause/move times, other robots may have high idle time since the robots have finished their extrusion tasks. Since this work assumes near equal distribution of work, the idle time is only dependent on the pausing and moving of robots and, therefore, not included in the fitness function.
This leaves only the moving and pausing time as part of the optimization function. This makes sense as the moving time is typically the only objective for a TSP problem, while the pausing time is typical in a multi-robot motion planning optimization problem. However, calculating the moving and pause time is not trivial since it is different per robot. There are a few options for formulating the fitness function with movement and pause time by taking the maximum or the average value of all robots' move/pause times. We only evaluate maximum and average since we want to optimize the values towards zero with all robots having equal weight of punishment for poor sequences. To test this, we formulate 4 variations of the objective function to determine the best function (see Table 1).
In general, it was found that for printing a circle geometry with various numbers of robots, it is preferable to use the average values of move and pause time. Since the optimization method is heuristicdriven, it may not be possible to know exactly why one formulation is better than another. Using the maximum or average values tends to allow for fast convergence on a solution, and they both converge to similar values, although they may be different sequences. The main reason why we believe the average is the best is because it will take into account all robot's pause/move time as opposed to using the maximum value may cause dramatic shifts in magnitude, making it difficult to converge.
For the remainder of this work, anytime a fitness function is used to calculate an optimum, the average move and average pause time will be used for evaluation.
In Section 4.3 two different methods were discussed to evaluate the printing time, namely, area-based and path-based. The path-based method more closely imitates the real robot setup which involves traversal of the printing path to calculate collision but is therefore much more time intensive to calculate the time to print. On the other hand, the area-based method makes assumptions about robot printing time, collision avoidance, and relies on pre-computing for avoiding collisions (cell-to-cell distance rather than point-to-point). While the area-based method is inherently faster, it is important to determine whether they both reach the same maximum and do not result in a different time to print, meaning the optimal of one should be the optimal of the other.
The first test consists of running 30 different optimizations (different initial sequences) using the area-based time method and stopping the optimization early (after only 50 generations) which is before convergence. We then compare the time to print, pause, and moving time using path and area-based methods. By comparing the difference we can observe if the area-based is a reasonable estimation of the pathbased method. The hypothesis that will be testing is that the total time to print, pause time, and move time are statically similar and then when a linear regression is run will have a slope equal to or near 1.
Different fitness functions printing a circle with 3 robots using Medial clustering (robot setup shown in Fig. 9 right).
Fitness 3 Robots -Medial clustering 3 Robots -Radial clustering function Number of Time to Number of Time to generations print generations print Max(Move) + Max(Pause) 103 17.323 159 17.738 Max(Move) + Avg(Pause) 60 18.046 171 16.987 Avg(Move) + Max(Pause) 72 18.159 143 17.602 Avg(Move) + Avg(Pause) 57 17.387 132 17.176
Fig. 9. By computing starting the optimization process from 30 different starting populations we compare the time to print, move, and pause time using the area based (x-axis) and the path based (y-axis) with 2 and 3 robots which shows that they produce similar results. The dashed line on each graph is the ideal case if the two methods matched exactly. Note here that a high 𝑅 2 -score value is not desired, and the slope value is less than 1 (being conservative) is the goal in evaluating the performance.
We compare two cases where we are printing a circle with 2 and 3 robots with various cell sizes. The results show that, at least qualitatively, the time to print matches for 2 and 3 robots and similarly with the pause and print time, this can be seen since all points lie near the dotted line which is the ideal case if the two mirrored each other exactly (Fig. 9). One example is with two robots, all sequences resulted in zero pause time with both area based and path based which was exact matching of the 150 data points. Additionally it is seen that the move time generally increases for smaller cells because as the number of cells increases the traversal between cells must also increase, this is generally true for TSP problems if more cities are added.
While a qualitative comparison is good, linear regression can also be run on the data to determine the slope and the variance. With two robots and the total time to print the slope is 0.95 and an 𝑅 2 = 0.84, the move time has a slope of 1.0528 and an 𝑅 2 = 0.98, since the pause time data all lies on the same point linear regression cannot be run but the data matches exactly. For three robots, the total print time has a slope of 1.001 and an 𝑅 2 = 0.99, the pause time has a slope of 0.1577 and an 𝑅 2 = 0.09, and the move time has a slope of 1.0174 and an 𝑅 2 = 0.93. Note that lower 𝑅 2 values for the pause time indicates a conservative estimate in that it include false positives for collision avoidance (i.e. some cases where collision will not occur may also be included). This is primarily because of calculating distances between centers of cells rather than nozzle-to-nozzle distance in the path based time calculation. Overall, this works in our favor since optimal printing will result in robots staying further away from one another. Both the move and total time to print for three robots shows very good agreement between the area and the path based. From the two and three robot case, we fail to reject the null hypothesis, meaning the data does not provide strong enough evidence to conclude that the times differ significantly between the two cases. This test however does not provide evidence that if given the same starting sequence, both methods would converge on the same optimal sequence. While both methods do not converge on the same sequence, we show that they converge to similar sequences. In order to compare two sequences, we will use Kendall 𝜏 𝑎 measurement which compares two sequences on a scale -1 to 1 where 1 is exactly the same and -1 is nothing similar. Running the same setup as before we compare difference in cell printing (inserting the actual cells for printing into the cluster ordering) with 2 robots printing, Kendall's 𝜏 𝑎 = 0.6398 for Robot 1 and 𝜏 𝑎 = 0.1573 for robot 2. While these are not close to 1 and therefore we cannot say they are similar, they are also not completely dissimilar. From further testing it appears as though the sequences deviate in the first generation, meaning in the initial generation of 5 chromosomes, different sequences are preferred due to slight differences and mutated to continue. While it cannot be shown that the two different methods arrive at the same optimal sequence, the optimal they arrive at appear to also be near the optimal of another.
For the remainder of this paper, any time the time to print is used to optimize, the area based method is used unless otherwise specified.
One of the main contributions of this paper is the parameterization of tasks for multi-robot printing. To demonstrate the importance of parameterization, we compare the two different clustering strategies presented in this work, namely medial and radial, to a non-clustered strategy which would be a traditional MTSP problem without clustering. This will also be compared to a combination approach where one robot may use radial parameterization, and the other may use medial parameterization. The hypothesis that will be tested is that the combination parameterization approach will achieve the best optimal solution and that the non-clustered will preform the worst. This will Fig. 10. The convergence graph for different parameterization, radial only, medial only, and a combination of radial and medial all converge faster and to a better optimum than a non-clustered sequence.
be verified with various number of robots printing a circular layer geometry.
In the case of two robots, the non-clustered method took the predetermined maximum of 5000 generations and did not converge, the three parameterized methods took only a few hundred generations. The medial/radial combination and the radial strategies ended up at nearly identical optimal values, but the combination took 200 fewer generations (Fig. 10(a)). For three robots, the non-clustered does not converge (50 generations without change in optimism) at all, performing poorly. In this case, the combination reached a better optimum than the radial or medial parameterizations but took a few more generations (Fig. 10(b)). This trend continues for four robots, where the combination performs much better than the two parameterized models and the non-clustered takes the maximum of 5000 generations and did not converge as well as not reaching a value near the parameterized approaches (Fig. 10(c)).
From this testing, we can determine that the two parameterizations presented in this work perform much better than the non-clustered approach. Furthermore, combining these two parameterizations and allowing for the GA to change between them as needed generally results in better or equal results. While we only show with a single setup of geometry and a simple approach to robot positioning, future results reinforce this result.
The final numerical study conducted is to determine the effect of different cell sizes, which can vary depending on the tessellation of the printing part. Prior studies on C3DP show that changing the cell size from larger to small leads to an increase in printing time as the movement time of the robots also increases [6,7]. However, these works have different collision avoidance protocols which our work does not have. The main question that we are exploring is whether smaller cells allow for better sequencing such that pausing occurs less often, despite a potential increase in move time. This test will be conducted on a large part of size 1.5 m × 300 mm with 5 robots in various orientations around the part (Fig. 11). The robot setup is based on the constraints of our physical robot setup. The cell size will be sampled near the value of the radius of the nozzle (40 mm), which is the safe distance that the robots must maintain from each other. Although the parameterization test demonstrated that a combination of medial and radial is best, we will use the same four testing methods as before for completeness. We utilize several different orientations of robots to validate that the results are not due to a specific robot placement/positioning. To see further results, see the appendix (Supp. Fig. 1).
The first observation for the first orientation of robots is that, as the cell size decreases, the time to print increases for the nonclustered (Fig. 11 a). However, the medial, radial, and combination parameterization it generally decreases or stays the same as the cell size decreases. When the cell size decreases, the combination parameterization becomes the best by a wide margin, while non-clustered becomes the worst. After closer consideration this makes sense because with smaller cells there are more cells to sequence and non-clustered methods will have trouble converging while parameterized methods may be able to find sequences which do not have pause time.
The second orientation has a similar result. As the cell size decreases, the non-clustered becomes much worse in overall time to print, while the parameterized methods remain similar in overall print time. While the combination parameterization does not perform significantly better than the other methods, it is still the best in most cases.
Although this set of tests was done on a limited setup (e.g., geometry and number of robots), it demonstrates a result contrary to other cellularized work where smaller cells do not require a trade-off to higher print time. These results further demonstrate that overall parameterized clustering aids in finding a better optimal than a non-clustered methods.
The goal of our physical testing is twofold. First, we aim to determine the impact of parameterization for printing. Secondly, we physically validate that the algorithm and sequencing strategy can be utilized sequentially, on separate layers, to slice and print a multi-layer object with multiple robots.
The first test consists of physically printing a banana shape with dimensions of approximately 850 mm × 270 mm using 3 robots. We only print a single layer and compare the layer time differences. Three different methods for printing will be tested, namely:
1. Slicer-based: Using a standard slicer, contours are made from the outside of the part and working inside (Fig. 12a). 2. Cellular Non-Optimized: The cellularization approach presented in this paper and all cells are parameterized using the medial approach, but the printing sequence is not optimized. The printing order is the identity where the first cluster is printed first, second cluster is printed second, etc. (Fig. 12b). 3. Cellular Optimized: Similar to (b), however, in this case the sequence is determined by our GA optimizer, allowing some robots to use radial sequencing and others to use medial sequencing (Fig. 12c).
One layer of the banana-shaped part is printed using the robot setup described in prior work [1,8,61] (Fig. 14(a)) (see Table 2). Fig. 11. The time to print a large object (1500 mm × 300 mm) with 5 robots in two positions is shown. As the cell size decreases, the time to print for the parameterized methods decreases or remains similar. Fig. 12. The printing path for the robots using different strategies: the slicer based approach (a), the non-optimized identity cellular approach (b), and the optimized cellular approach which contains medial and radial contours.
The time taken for the different printing methods chosen in physical evaluation. The optimized cellular approach takes less time than the non-optimized cellular and the slicer-based contour method. Examples of the path are in Fig. 12.
Parameterization Time to method Print (min) Slicer-based contour method 95.29 Cellularization without optimization 91.66 Cellularization with optimization 86.77
The key result from the testing is that the circularized optimized result took the lowest overall time to print and was 8.5 min faster than the slicer-based method, and 4.89 min faster than the non-optimized approach. While this time saving may seem insignificant, as a percentage of the overall time to print, each robot has more than 43 min of extruding time which cannot be reduced. For a typical part with several hundreds of layers, say 200 layers, an average saving of 8.5 min per layer results in more than 24 hours of reduction in time to print. Of course, an exact extrapolation of time saved across layers must consider other factors such as inter-layer transitional issues, accumulation and correction of printing errors, and specific issues pertaining to the geometric and topological complexity of the part itself. Therefore, a potential time reduction on the order of magnitude of hours could be reasonably expected.
The time each robot spends extruding can also be compared. In the case of the slicer-based the printing occurred almost sequentially where the robots printing in red and in blue printed their work in parallel and then the robot printing in black finished its work while no other robots were working (Fig. 13a). This is much different from the other two cases where all three robots worked together for a significant portion of the print. In the non-optimized case, the robot printing in black has to wait until the red and blue robot are far enough away until it can start printing. In the case of the optimized all robots begin by printing together and then black robot pauses for a short time before resuming its work.
We further validate the algorithm by printing multi-layer objects. In order to complete the optimization on the objects, we first equally distribute the workload between robots by moving the part within the workspace, since workload changes between layers we sum the total work for each robot across all layers and move the part until work is equally distributed between robots. While this approach does not guarantee that all layers have equal distribution, it ensures equal work across all layers. This is also similar to the approach by NoodlePrint with their concurrence measure but without the penalty for interfacing cells [6]. The optimization can then be performed sequentially on every layer to determine the printing order for the clusters. Additionally, each robot may change between radial and medial, as shown in the chromosome representation.
A model of Einstein was chosen to print with approximate dimensions of 475 mm × 290 mm × 135 mm. During printing, it was observed that the segmentation method using truncated octahedron the interfaces were always planar but changed direction every few layers (Fig. 14(b)).
One of the main results of this work is demonstrating that parameterization of the space allows for better sequencing of tasks to avoid Fig. 14. The times when extrusion occurred during the physical printing of the banana per robot (red, blue, black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
collision and work in a time-optimal manner. This was first demonstrated with the numerical results, where the non-clustered performs poorly in overall time to print and struggles to converge on an optimal schedule (Section 5.5). The parameterized methods allow for much faster convergence. For printing a circular region with 3 robots, the combination approach found an optimal that takes 402.96 min less than the random, non-clustered approach. The combination takes 22.79 min and 25.79 min lesser than the medial-only and radial-only approaches respectively. This reduction in print time also increases with the number of robots. This was further demonstrated in the physical results where the slicer-based contour method took 9.52 min more to print a single layer than the cellularized optimized approach. If applied over an entire part this can potentially lead to reduction of print time on the order of magnitude of hours.
In order to reduce the computational effort to calculate total print time, the area-based evaluation method was introduced, which uses a conservative cell-to-cell distance calculation. The cell-to-cell distance can be precomputed and used throughout the optimization process, significantly reducing overall optimization time. It was also shown that computing the optimal using the path-based method was similar in overall time to print, idle time, and movement time compared to the path-based method. This result shows that future studies can reduce the computational effort when evaluating the fitness function by using the area-based method, which allows for more efficient optimization without compromising the accuracy.
Other works in C3DP that use a safe-area collision avoidance algorithmic approach make all other robots pause while one robot is printing an unsafe region and have scalability difficulties since there is little time for all robots to work together [6,7]. The Traveling Cellsman framework allows for all robots to start printing and working on their regions, while in some configurations printing cannot happen continuously, so the amount of pause time is minimized. This was demonstrated in the numerical tests, where it was shown that the cell size is relatively independent of the time to print since smaller cells have higher movement time but lower pause time due to better optimums. This result was unexpected since prior works have shown that smaller cells have higher print time due to higher movement time but no method to reduce the pause time [6,7]. We also demonstrated with 3 robots that can print a large multi-layer object by conducting optimization on every layer in order to allow for all robots to collaborate during printing.
One of the observations from physical printing was that there were gaps between the pieces printed by different robots which can lead to weak interface [6]. These gaps could have been caused by initial calibration or compounding small errors in the robot motor positioning throughout the printing process. While using some methods of segmentation this can lead to disastrous results, using the truncated octahedron allowed for maintaining cohesion between the parts printed by different robots. This is due the angled interface between truncated octahedron constantly change direction (Fig. 14(b)). Additionally, the truncated octahedron creates cells on every layer which are either squares, or hexagons which can be much easier to print without harsh angles.
While this work has focused on collaborative 3D printing, other multi-robot task scheduling applications could also be used with parameterization. As an example, colored MTSP exist in milling applications where multiple robots must work in shared spaces [12]. In this way, the parameterization shown here can be used to sequence contours, where each contour is the bit size. Other potential applications can also include welding, WAAM 3D printing, painting/spraying, manipulation/assembly, or surface prep applications where multiple robots share a working space. In these cases, optimizing the scheduling between tasks is critical since inefficient movements can result in wasted time or inefficient use of resources.
There are four important limitations of our work. First, our work does not address the important problem of determining the optimal placement of robots. This is important given that the location of the robots directly affects the partition of which leads to a completely different assignment, clustering, and workload of robots. For instance, the orientation of 5 robots, which results in the best time to print changes based on the cell diameter (Fig. 11). It is also possible to employ methods such as the ones by Stone et al. [10] to simplify our optimization. For example, for two adjacent robots printing medial clusters, they cannot print their innermost cluster at the same time. We believe this is a challenging problem that should be addressed in future work. As second related issue is that the positioning of robots is important to ensure equal (or nearly equal) workloads for further minimization of print time (having severe workload discrepancies leads to high idle time). We see interesting possibilities for future work where the restrictions to remain within one layer could be removed or three-dimensional parameterization could be introduced. Having mobile robots is yet another potential venue in this regard. The third limitation of our work is that we did not demonstrate our method for heterogeneity of robot kinematics. While a system of robots with different kinematic configurations does not effect our broader conceptual foundation defined by parametric clustering, it is bound to affect the optimization of the makespan as well as communication protocol owing to differences in path planning. Finally, our method currently assumes layers to be planar. We believe that generalizing the framework for the partitioning, parameterization, and clustering to curved layers is a particularly intriguing prospect for future research.
This work has that using parameterized clustering can be useful in creating the schedule for printing cells in cooperative 3D printing tasks. In order to create clusters efficiently, we introduced radial and medial clustering, which are formed based on the partitioning of space and its subsequent cellularization. The numerical results showed that clustering allowed for faster convergence to an optimal solution as opposed to a non-clustered approach. Additionally, allowing for different clustering methods (radial or medial) per robot allowed for better solutions. Finally, the physical results validated that the optimization process reduces the overall time to print as compared to a non-optimized ''moving-front'' style algorithm as well a slicer-based algorithm, traditionally used in C3DP applications.
The key advantage of our methodology is that it achieves strikes a balance between complexity and performance. As an example, by reducing the design space from 7.02 * 10 215 to 5.89 * 10 14 , the solution space can be more effectively navigated and optimized. This allows for saving nearly 10 min per layer. When extrapolated to 3D parts with a few hundreds of layers, this could potentially lead to time savings on the order of hours taking into account the time taken for inter-layer transitions. Additionally, while we only demonstrate our approach for C3DP applications, it can be applicable to many different types of multirobot applications where traversal between tasks must avoid collision while also minimizing total time moving between tasks. This is particularly relevant for multi-robot additive manufacturing applications, as our framework is agnostic to the process type (e.g., FDM or cold-spray).
Additive Manufacturing 112 (2025) 104987