SafeZone*: A Graph-Based and Time-Optimal Cooperative 3D Printing Framework Free

Swarm manufacturing (SM) has garnered significant attention in recent years for its potential to transform traditional manufacturing processes by employing a swarm of robots to cooperate on manufacturing tasks, particularly in the field of additive manufacturing (AM) [1,2]. Cooperative 3D printing (C3DP), a specialized form of swarm manufacturing tailored to AM, uses multiple 3D printers (e.g., fused deposition modeling printers) to work on a part simultaneously [3]. The printers can be either stationary or mobile, depending on the system configuration and desired application, but they typically operate in a shared workspace. Although working in a shared workspace is necessary to ensure cooperation and leverage the benefits of SM, it introduces two new challenges, namely division of labor and collision avoidance. In the context of C3DP, division of labor is concerned with how to divide a part into sub-volumes and then assign them to each printer. Collision avoidance refers to the ways (e.g., through control, scheduling, communication, detection) that the printers cooperate without colliding either with each other or with the part being printed, which could occur due to printer’s nozzle geometry or kinematic chain, as in the case of C3DP with robotic manipulators. The approach to solving these problems also heavily depends on the dynamics of the printing robots, since C3DP has been achieved with drones [4], mobile robots [5], and robotic manipulators [6].

Minimizing the overall printing time (i.e., makespan) is a standard objective in C3DP [7]. Due to the inherent interdependencies between the different stages of C3DP, an optimal framework must co-consider division of labor and collision avoidance, ensuring that each informs the other. For instance, if work is evenly distributed among all robots, but the collision-avoidance methodology is simplistic, the process may incur unavoidable collisions. Conversely, if the collision-avoidance strategy is robust and accommodates any division of labor, it may lead to suboptimal printing times. Further complexity is added to the problem when physical constraints are enforced. Such constraints are particularly relevant in systems with robotic manipulators, which may have issues with reachability and positioning of the printers in the workspace.

Existing C3DP frameworks primarily focus on the geometric aspects of the division of labor (i.e., size of sub-volumes) when minimizing printing time, often by attempting to evenly distribute labor among robots. However, this approach overlooks the topological aspects of the division of labor (i.e., connectivity of sub-volumes), which can significantly impact scheduling and collision avoidance. In our previous work, we demonstrated the importance of integrating the topology of the division of labor into the design space when optimizing C3DP frameworks [8]. Building on this foundation, we introduce SafeZone* (SZ*), a Voronoi- and graph-based C3DP framework for collision-free, minimum-makespan printing. The key idea of SZ* is to abstract the complexity of all C3DP stages by formulating an optimization framework that depends solely on printer positioning. By establishing a structured mathematical relationship between printer placement, layer partitioning (geometry and topology), and collision-free scheduling, SZ* efficiently searches the design space to identify C3DP configurations that minimize makespan.

SZ* is novel in that it minimizes makespan by co-considering the geometric and topological implications of the division of labor (Fig. 1). Particularly, SZ* demonstrates how an integrated topological-geometric analysis, which relates Voronoi tessellations to partition graphs, enables not only the efficient generation of optimal printing schedules, but also fundamental insights into beneficial topologies for C3DP, based on the graph theory. Although SZ* is applicable to generic C3DP systems, it is tailored to robotic manipulator applications, where reachability and printer positioning are the primary concerns. Finally, SZ* is a hybrid approach: while the printing path and schedule are pre-computed, execution is coordinated in real-time by a centralized control hub that sends commands to the printers. As such, the main contributions of this work are as follows:

1. We establish a mathematical framework for representing and evaluating C3DP layer topology based on the graph theory, with a focus on graph coloring.

2. We introduce SZ*, a novel C3DP framework for makespan minimization that co-considers geometry and topology of a layer by relating Voronoi sites to partition graphs.

3. We numerically validate SZ* and establish the chromatic number as a key metric for the topology of a minimum-makespan C3DP partition.

4. We physically validate SZ* and demonstrate its robustness in constrained manufacturing environments by leveraging virtual robot positions (VRPs).

The remainder of the article is structured as follows. Section 2 reviews related research on C3DP and multi-robot systems that face similar problems. In Sec. 3, we establish the mathematical framework defining layer geometry, layer topology, and the key concept of buffer regions. The algorithmic framework of the SZ* for achieving collision-free and minimum-makespan printing is presented in Sec. 4. In Sec. 5, we show the numerical results using the framework in a variety of scenarios, including different numbers of printers and layer geometries. Additionally, we present a case study with four selective compliance articulated robotic arms (SCARA) printers as physical validation. Finally, we discuss limitations, future research directions, and topics of interest in Sec. 6.

While there has been extensive research on additive manufacturing [9], the concept of C3DP is a recent development [10,11]. There are several key challenges within C3DP that are not considered in additive manufacturing or robotics in general. With traditional single-robot or gantry-based 3D printing, collision avoidance between robots is not a concern. In many robotic tasks where multiple robots work together to complete a task, it is often preferable for robots to work as independently as possible, with little to no overlap in workspaces [12–14]. This can be seen in many coverage tasks where robots are required to cover the entire space, commonly through discovery and exploration [14,15]. In many cases, it is preferable that the robots interact as little as possible. In a work by Lin et al., two robots explore an unknown space, and in the ideal scenario, the paths of the robots never cross, and the robots, in general, have very little interaction with one another. A similar study demonstrates that efficient exploration of a known space can be achieved through the cellularization of the space, where each robot explores its own assigned region [16]. This is similar to what is done in C3DP, where robots are assigned known regions, and the method to segment the part into regions involves equal area distribution and, thus, nearly equal printing time. This approach is also seen in similar coverage works, such multi-robot car painting. In a work by Zbiss et al., several multi-axis robots work together to paint a car simultaneously. In this case, robots have shared workspaces, and collision avoidance is managed through careful ordering of the painting paths [17].

When multiple robots operate in a shared environment, preventing collisions is essential. While various methods exist, we explore a few relevant to our work including, scheduling, zone-based approaches, and real-time collision detection.

By carefully scheduling printing tasks, collision can be avoided before the printing process begins. This approach has been used in C3DP by scheduling printing paths [18,19] or sub-volumes [20–22]. In these cases, the relationship between part geometry and printing order must be carefully considered, as an overhanging region cannot before the region beneath it. One previously developed method uses dependency trees to determine printing order by leveraging the interfacial boundary of segments to construct the dependency tree [23]. The use of mobile robots and even aerial robots introduces additional dependencies, requiring careful coordination of printing sequences.

Another method for avoiding collisions is to classify regions as safe or unsafe, meaning that if printing occurs in one area, then all other robots either cannot print at all or are restricted to printing in specific regions [24–27]. While different articles use slightly varying terminology, the core idea remains the same: define regions where collisions are possible, typically along boundaries shared by multiple robots. These unsafe regions can then be managed by either printing them sequentially or restricting printing to a single robot at a time. In the case of Shen et al. [27], the division of labor was determined by optimizing the placement of cuts in an asymmetric part to evenly distribute work. Using this segmentation, an offset of the interfacial boundary was created to separate safe and unsafe regions [26,27]. Other approaches rely on predefined cellularization to define safe regions, where certain cells are labeled as safe or unsafe. In these methods, printing occurs sequentially in unsafe regions before transitioning to collaborative printing in safe regions (referred to as interfacing or non-interfacing in their terminology) [24,25].

Real-time collision avoidance can also be employed in C3DP, where robots continuously relay their current positions and/or trajectories to other robots. If a collision is imminent, one or more robots halt their work and resume once it is safe to proceed [26]. In such cases, optimizing the printing process becomes challenging since the path may not be predetermined, particularly when using a standard slicer for path generation. These methods can be beneficial in scenarios where the printing time or task completion time is unknown or variable, as other robots can continue working even if one is temporarily halted. However, scaling these methods to an arbitrary number of robots is difficult, as communication overhead and real-time distance calculations become increasingly costly with a larger number of robots.

While this is a brief overview of existing division of labor and collision-free scheduling methods, SZ* takes a fundamentally different approach by actively driving the partitioning process to ensure an optimal schedule. Unlike prior works that derive printing schedules from predefined partitions, SZ* simultaneously optimizes partitioning and scheduling through an integrated geometric-topological analysis of the C3DP layer. By leveraging Voronoi tessellations and partition graphs, it intelligently defines safe and unsafe regions while structuring them to enable collaborative, rather than strictly sequential, printing. This key distinction allows SZ* to minimize makespan while ensuring collision-free multi-robot coordination, setting it apart from existing methods.

In this section, we formulate the problem of interest. Given a layer geometry, in the form of boundary curves, and a set of printers, the question is: how to partition the layer geometry into sub-regions and assign them to each printer? Ideally, the choice of partition should facilitate collision-free movement of the printers and, consequently, minimize printing time. While the geometry of each sub-region (i.e., size and shape) is clearly an important aspect to be determined, SZ* also considers the topology of the partitioning (i.e., how the sub-regions are connected to each other). In fact, current literature shows that the partitioning topology of robots working in a shared workspace is understudied.

To better understand the implication of the topology of the layer partitioning, particularly in the context of C3DP, let us consider a layer with a circular boundary that is set to be printed by four printers. A symmetric partition would be a simple, yet naive way to create the sub-regions (Fig. 2(a)). In this case, the likelihood of collision is high as we approach the center of the circle since that neighborhood is shared among all four printers. However, we observe that there may be other partitions of the same circular layer such that only three robots share any existing boundary (Fig. 2(b)). In fact, it is even possible to create sub-regions, such that no more than two robots share a boundary (Fig. 2(c)). Knowing that we are interested in achieving time-optimal and collision-free cooperative printing, these fundamental differences between each of the aforementioned partitions reveal the criticality of assessing their topology and geometry in conjunction. Hence, the goal of SZ* is to tackle these two aspects simultaneously.

The partition topology of any layer geometry can be represented as a graph $GP$⁠, denoted as a partition graph, which is algebraically represented using an adjacency matrix $MP$⁠. The nodes of $GP$ represent the sub-regions, and the edges between the nodes represent the boundaries between sub-regions. Two sub-regions are said to be connected if they share a boundary, which can be a curve or a single point (Figs. 3(a)–3(c)). Generally, we can see that graphs that have nodes with lower degrees are better, in the sense that robots share fewer boundaries. However, we are interested in a more quantitative metric for judging whether a partition, topologically, is better than another. We approach this problem by relating the partition graph to the printing schedule, which, in turn, defines the printing time.

In this idealized scenario, topology becomes irrelevant, and geometry dominates the choice of partition. To consider a more realistic scenario, let us assume that, if two sub-regions share a boundary, their respective printers cannot work on them simultaneously. As shown in our literature review, this is a common approach for achieving collision-free printing [26,27]. Under this assumption, we notice that certain layer topologies can offer parallel printing opportunities. For example, if two sub-regions do not share a boundary, then they can be printed completely in parallel without concern of collision. This is where we can leverage the concept of graph coloring [28] to guide us in drawing fundamental conclusions about the goodness of the topology of a partition.

Let printer $Pi$ be responsible for area $Ai$⁠. Let $S$ be the printing schedule, which we define as a sequence of tuples $τj$⁠. Each $τj$ may consist of one or more $Ai$ which can be printed in parallel and do not share edges with one another. If a $τj$ has multiple $Ai$⁠, it simply means that those regions are printed in parallel during that step of the schedule. For instance, in the idealized scenario we just discussed (Eq. (1)), S = (⁠$τ1$⁠), where $τ1=(A1,A2,A3,A4)$⁠, meaning that everything is parallelized.

Assuming that we are given a generic partition, we can represent it as a partition graph $GP$ and use any graph coloring technique to assign a color to each node such that no adjacent nodes have the same color. The coloring process is the same as determining a collision-free printing schedule. Thus, in our context, each color represents a tuple $τj$ in the schedule $S$⁠. The number of nodes with a given color represents the size of $τj$⁠, which is the number of sub-regions printed in parallel at the $jth$ step of $S$⁠.

We established, based on the printing schedule, and assuming that the order of the graph is equal to the number of printers, that the optimal topology generates a partition graph $GP$ with the lowest possible chromatic number $χ(GP)$⁠. From graph theory, we know that, when considering the family of all possible $n$ partitions for connected geometries, $χ(GP)$ is bounded between 2 and $n$⁠. If $χ(GP)$ is equal to $n$⁠, it means that the resulting partition graph is complete $(Kn)$⁠. This is equivalent to the symmetric case (Fig. 4(a)), which produces the worst possible schedule for $n$ printers (Eq. (2)). If $χ(GP)$ is equal to 2, $GP$ becomes a path graph $(Pn$⁠), which is equivalent to the third example provided (Fig. 4(c)).

Although topologies with $Pn$ partition graphs have the smallest chromatic number, which would make them theoretically optimal from a printing scheduling perspective, they are generally not viable in C3DP systems with robotic manipulators due to scalability and reachability issues. These are two key constraints that SZ* considers in the formulation of its partitioning framework (Sec. 4.1). For instance, such a topology could produce sub-regions that stretch very far from the original positions of the printers. Additionally, as more printers are added, and the part of interest is not slender, the workspace may become too cluttered as the distances between the robot positions decreases (Fig. 5). However, there are special cases in which we can leverage other families of partition graphs that have $χ(GP)=2$⁠, in feasible ways (Sec. 6.1) but they are difficult to generalize to any layer geometry.

As such, there are several ways in which a layer geometry can be partitioned to allow for a desired topology. However, there are three main concerns to be addressed specifically for C3DP with robotic manipulators. First, we need to minimize the likelihood of collision between robots that share a common boundary. Second, the method should provide a means to avoid reachability issues. Finally, the method to generate the partitions should offer a natural way to control the partition topology, especially the degree-3 condition. Therefore, a systematic way of determining and controlling the geometry of the sub-regions $Ai$ is necessary.

We note that constructing a Voronoi tessellation of the planar domain using the robot locations as the Voronoi sites naturally alleviates all of these concerns (Fig. 7). Given a set of points ${pi}∈Rn$⁠, a Voronoi tessellation is a mutually exclusive and exhaustive set of polygonal regions ${Ai}$ such that any point $q∈Ai$ is closest to $pi$⁠. If ${pi}$ are robot positions, then using such a partition to segment the layer geometry will automatically result in printable sub-regions that are also closest to the nearest robot. This addresses the first two concerns, namely minimizing the likelihood of collision between robots that share a common boundary in the partition and naturally avoiding reachability concerns.

Besides the closest-site property, Voronoi tessellations offer a natural way to control the topology of the partitioning at the vertices. To understand this, we consider the Delaunay diagram, which is the dual of the Voronoi tessellation (i.e., the vertices of Delaunay diagram are the sites of the Voronoi tessellation). The Delaunay diagram can consist of convex polygons that can all be triangulated to obtain the well-known Delaunay triangulation. Also note that all polygons in a Delaunay diagram are cyclic polygons (i.e., each triangle in a Delaunay triangulation has a unique circumcircle). In the general case, a Delaunay diagram may not have a unique triangulation. Consider a case where a quadrilateral cell can be split into two triangles in two different ways (by choosing to connect either of the diagonal vertices with an edge). However, when there are no more than three Voronoi sites that are cyclic, we always obtain a unique Delaunay triangulation [29]. This also means that the resulting Voronoi tessellation will always contain degree-3 vertices. This offers a principled approach to generating robot positions that ensures degree-3 partitioning. Therefore, we employ Voronoi tessellations to simultaneously control the geometry and the topology of the partition, as it naturally alleviates the key concerns in C3DP with robotic manipulators.

Creating a collision-free printing schedule based solely on the sub-regions $Ai$ is too conservative. Collisions can generally only occur in the vicinity of shared boundaries of the partition, meaning robots can safely print in areas away from these interfaces. Hence, preventing two adjacent sub-regions from being printed in parallel due to a small (relative to $Ai$⁠) shared boundary area is not ideal. This problem is intensified in C3DP applications due to the large size of the $Ai$ sub-regions, as the ratio of the shared boundary area to the total area of $Ai$ becomes smaller. To address this issue, we introduce the idea of buffer regions to improve the original printing schedule.

Recall that the collision-free schedule we derived for the symmetric partition of the circular layer (Fig. 2(a)), has the same printing time as if only one printer was used. This is a clear example of how the printing schedule based solely on $Ai$ can perform poorly due to its conservativeness. Therefore, we can improve the schedule, and thus the layer printing time for any partition, regardless of topology or geometry, by introducing buffer regions $Bi$⁠. The idea is that the buffer regions will act as buffers between the original $Ai$ sub-regions, such that they are the only regions of the layer geometry where collisions between robots may occur. This is guaranteed if the printer that was originally assigned to print sub-region $Ai$ is also assigned to its corresponding buffer $Bi$⁠, and the buffer region geometry is properly defined based on the printer dynamics. In general, determining the geometry of buffer regions requires knowledge of the mechanical model of the printers, that is, whether they are drones, mobile robots, or robotic arms [26]. However, regardless of the mechanical model, at a minimum one must consider the geometry of the nozzle when testing for collisions.

This process can be repeated for every $Ai$ until all necessary buffer regions are created. Note that if $∂Ai$ is disjoint, so is $Bi$⁠, although the buffer region is still treated as one singular node in the partition graph. Additionally, there is the possibility that $Bi$ is empty, which can happen if the resulting intersection itself is empty. In this case, the node representing that particular $Bi$ still exists in the partition graph, but it has an area value of zero.

The algebraic form of the resulting adjacency matrix $MP*$ (Eq. (7)) is constant irrespective of the topology of the original partition graph. This result is important to SZ* as it allows the incorporation of an efficient heuristic-based scheduler to the current framework (Sec. 4.1.3) since the sub-graph formed by the $Bi$ regions will always be three-colorable, provided that the original partition graph is also three-colorable. Therefore, once the buffer regions are efficiently printed using the three-coloring schedule, all $Ai*$ can be safely printed in parallel during the final step of the schedule.

The overarching framework for SZ* consists of two core components

1. Co-solving of partitioning and scheduling

2. Printer communication

The framework is modular in that it also allows for user-defined solutions. For example, if the user provides their own partitioning scheme by pre-defining the $Ai$ sub-regions, SZ* is only responsible for creating the buffer regions (if desired), scheduling the printing, and overseeing the process through the established communication protocol using wait-and-go commands [26] (Fig. 10). However, in general, SZ* is responsible for co-solving partitioning and scheduling by co-considering layer geometry and topology in its design space exploration. The user can be interested in an user-defined partition for many reasons, but mainly because of the design of the part. For instance, if there is a multi-material process (PLA and ABS), then the partitions could be determined by the intended material deposition areas. The user might also want to avoid certain partitions that have boundaries in mechanically important regions, such as critical features that will be subjected to elevated stress and strain concentrations.

The design space is defined to be the set of Voronoi sites ${pi}$⁠, which correspond to the printer locations (Sec. 3.2). The locations of the Voronoi sites directly influence the structure of geometry and topology at each layer through their associated Voronoi tessellation and resulting partition graph (Sec. 3.3). As discussed, we are interested in minimizing the makespan (⁠$Tl$⁠) of a layer, which can be extracted from the schedule of the partition graph. Thus, we can frame the choice of Voronoi sites as a constrained optimization problem, in which the objective function is the makespan function $C(pi)$ of the partition graph. Before formally defining the optimization problem, we first describe the nature of the existing constraints.

Within the context of C3DP, it is not uncommon for some robots to have constraints on physical locations and cannot be arbitrarily placed, limiting the placement of the Voronoi sites. As the location of the printers is directly related to the resulting partitioning, we call these constraints partitioning constraints. SZ* considers three such constraints when exploring the design space.

• Reachability: Given a radius $rp,i$⁠, defining a two-dimensional circular workspace $Cp,i$ centered at the Voronoi site $pi∈R$⁠, and given a sub-region $Ai$⁠, constructed as the intersection of the Voronoi cell of $pi$ and the layer geometry, then the reachability constraint states that $Ai$ must be completely enclosed in $Cp,i$

  $∀i,Ai⊆Cp,i$

  (10)

  guaranteeing that the printers are able to reach their assigned sub-regions. Note that printers are also guaranteed to reach their buffer regions $Bi$⁠, if they are created, since each $Bi$ is enclosed within their $Ai$⁠.
• Separation: Given a radius $rs,i$⁠, defining a circle $Cs,i$ centered at the Voronoi site $pi∈R$⁠, then the separation constraint states that none of the circles $Cs,i$ should intersect with one another

  $∀i≠j,Cs,i∩Cs,j=∅$

  (11)

  thus, we can tune the radii $rs,i$ to create safe distances between each printer. This prevents printers to be placed too close to each other, which could result in collisions even before the printing process starts due to their spatial occupancy.
• Configuration Space: Given the configuration spaces $Ci∈R2$ defining feasible values for the Voronoi sites $pi$ coordinates, and given their corresponding obstacle space $Oi∈R2$ defined by the layer geometry, then the configuration space constraint enforces the following:

  $∀i,pi∈Ci∖Oi$

  (12)

  avoiding solutions in which the Voronoi sites $pi$ are inside of the layer geometry. If the layer geometry is not a genus-0 surface, meaning that it has inner holes, this constraint does not avoid the printer locations to be contained in those holes. However, if the user does not wish to consider those solutions, then the obstacle space $O$ can be taken to be the union of the layer geometry and any associated holes.

Therefore, this set of partitioning constraints allows the SZ* framework to generate optimal solutions while maintaining feasibility in practical C3DP applications. Once the Voronoi sites are selected, the buffer regions $Bi$ are created and the new partition graph $GP*$ is constructed (Sec. 3.3). Finally, the objective function must be evaluated, which involves extracting the layer printing time (makespan) from the partition graph schedule.

We can treat the nodes of $GP*$ as tasks with time duration equal to the time taken for a printer to finish printing a sub-region (i.e., a task) that a node represents. As we are interested in deriving the optimal schedule $S*$ from the partition graph, which is in itself another optimization problem, we must first define how to evaluate the makespan function $C(pi)$ and which physical constraints are involved.

Overall, SZ* considers two constraints in scheduling, which are derived from our previous layer topology analysis (Sec. 3.1), namely (1) that the printers cannot be working on two adjacent nodes at the same time, and (2) that each printer can only work on its associated safe region $Ai*$ and buffer region $Bi$⁠. This is a standard graph scheduling problem that can be solved with a mixed-integer linear program (MILP). The problem is mixed-integer since there are continuous decision variables defining the starting time of the tasks, which can be any real positive number, and also binary decision variables defining the order in which the tasks are sequenced. Within the SZ* framework, we solve this inner optimization problem using the computational infrastructure for operations research branch and cut (COIN-OR CBC) solver [30].

Alternatively, we can leverage the analysis for the layer printing time $Tl$ when adding buffer regions $Bi$ (Sec. 3.3) and simplify the scheduling problem using a heuristic-based approach. The idea is to divide the scheduling problem into two smaller sub-problems. The first sub-problem refers to scheduling only the buffer regions $Bi$⁠, while the second sub-problem involves scheduling the $Ai*$ sub-regions (Fig. 9). The second sub-problem is trivial as, by construction, all sub-regions $Ai*$ can be printed in parallel, and so no computational effort is needed. The first sub-problem can now be solved using the MILP optimizer, but considering only half the number of nodes when compared to the non-heuristic-based approach. Finally the schedules of the first and second sub-problems are appended, providing a new definition for $C(pi)$⁠. Although this approach does not guarantee an optimal solution, it is computationally efficient, and approximates the optimal makespan well when the size of the buffer region is small (Eq. (9)).

As discussed in the previous section, the makespan function can be the optimal makespan of the partition graph, which could be found using any MILP solver, or it could be heuristic-driven, such that the buffer regions $Bi$ are first printed using graph coloring, and then the safe regions $Ai*$ are printed in the final step. Irrespective, of how the makespan function $C(pi)$ is defined, the formulation of (Eq. (13)) remains the same.

To solve the optimization problem above, we employed a genetic algorithm (GA) formulation, where the gene is represented by the set of Voronoi site coordinates (which create the partition). The rationale for using a non-gradient-based optimizer is that the design space is not guaranteed to be continuous. Small perturbations in the location of the Voronoi sites may lead to partition graphs with distinct topologies creating discontinuities and a non-smooth makespan function. Therefore, a gradient-free, meta-heuristic optimizer such as the GA is an appropriate choice to solve (Eq. (13)). To facilitate convergence, the partitioning constraints $C1$ can be construed as penalty functions and added to C(⁠$pi$⁠). For instance, the reachability constraint can be redefined as a penalty function $w(d)$⁠, where $d$ is the amount of area assigned to printer $i$ that is outside of $Cp,i$⁠, guaranteeing a smoother makespan function for better convergence. A similar rationale follows for the remaining partitioning constraints.

While other gradient-free and meta-heuristic optimizers are possible, such as simulated annealing (SA) or particle swarm optimization (PSO), we chose GA as it is more well-suited to high-dimensional, non-smooth design spaces such as ours. SA is a single-point optimizer that can have difficulties in exploring high-dimensional design spaces. Additionally, the determination of the cooling function is not trivial, complicating the criteria for convergence and increasing the chances of local-minima solutions. PSO was originally conceived for smoother problems, as the velocity-driven updates may struggle in rugged design spaces, leading to premature convergence. Hence, while a separate study can be conducted comparing these algorithms for our specific design problem, we defer it to future work.

The last step in the SZ* framework is to establish the communication protocol during the printing process. The purpose of the protocol is to certify that the optimal schedule $S*$ is followed. We approach this problem from a centralized perspective, where a single control hub communicates independently with each printer and sends only two types of commands: wait or go [26]. The control hub follows the optimal schedule $S*$ generated through the partition graph $GP$⁠, which is a product of the Voronoi tessellation, to coordinate the sequence of wait and go commands. The steps are as follows:

1. Iterate over each tuple of sub-regions $τ$ in $S*$⁠. Stop if there are no more $τ$⁠.

2. Send the go command to each printer that has a region in $τ$⁠.

3. Once a printer finishes its region, send the wait command.

4. Once all printers that have regions in $τ$ receive the wait command, go back to step 1.

Since the regions in each $τ$ can be printed in parallel, this simple communication protocol guarantees a collision-free and optimal printing process. Note that SZ* is a framework that uses layer synchronization, meaning that the printers are not allowed to move on to the next layer until it is finished. This property is automatically enforced in our framework since each $τ$ can only contain regions in the same layer. Layer synchronization is applied because errors between predicted and actual time to print can cause uncertainty (e.g., accumulated error) in planning that is not considered in the scheduling. This preserves the collision-free and layer-optimal scheduling properties outlined in Sec. 4.

In this section, we present a set of numerical studies using the proposed framework, as well as a physical validation demonstration using 3DOF robotic manipulators arranged in a discrete grid. We demonstrate how SZ* is not only a valuable tool for practical C3DP applications, but also a powerful aid in exploring fundamental research questions regarding optimal partitioning and scheduling of 3D-printed parts.

In order to demonstrate the viability of the SZ* framework, we conduct computational experiments using the proposed methodology and we postulate two hypotheses.

1. Optimal Voronoi sites generate three-colorable partition graphs

2. The areas of the $Ai$ sub-regions will be approximately equal

To show generality, we ran SZ* using a wide range of layer geometries, both convex and concave. Specifically, we tested four layer geometries: square, triangle, z-shape, and a cross section of the famous Stanford bunny. Additionally, we varied the number of printers $n$⁠, and thus the number of Voronoi sites used, starting at $n=4$ and going up to $n=7$⁠. The GA hyperparameters were kept constant at 50 generations, 100 population size, 20% mate rate, and 10% mutation rate. The partition constraint parameters were set at $rs=0.5$ and $rp=2.0$ units.

When generating results with this set of hyperparameters, we observed no issues in the convergence of the GA optimizer. However, due to the nature of the physical constraints, the user must be careful to ensure that the chosen physical scenario is realistic (e.g., adequate number of printers for the available space), or the GA may struggle to find a feasible solution. We have also provided the performance data for all numerical studies in  Appendix. Finally, all studies were performed on an HP Elite Tower 880 G9 computer with a 12th Gen Intel® Core™ i9-12900 2.40 GHz processor, with 64 GB of RAM.

For each combination of layer geometry and number of printers, and using the optimal MILP scheduler, we extracted the chromatic number of the resulting partition graph and analyzed the evolution of the areas of the sub-regions using a number of metrics.

For almost all of the scenarios, with one exception, the optimal solution had a partition graph with chromatic number $χ(G)=3$ (Fig. 11). In the scenario using the z-shaped printing part, with four printers, SZ* managed to find a solution with $χ(G)=2$⁠, which is theoretically better, while satisfying all the constraints. This solution was possible due to the slender and concave nature of the printing part, allowing a path graph $P8$ to be constructed. These results confirm our first hypothesis and the assertion made in the conceptual framework that, in general, topologies that are three-colorable, or less if possible, are the best for C3DP.

To understand if there are underlying patterns in the layer geometry of these three-colorable partitions, which would allow us to better guide the optimization process, we studied how the normalized areas of the sub-regions evolved with the generations. In particular, we looked at how the areas of the individual $Ai*$ and $Bi$ sub-regions evolve, their sum, and their difference. For difference metrics, we looked at the difference between the maximum and the minimum areas, the standard deviation, the mean squared error, and the sum of the absolute differences. However, all metrics behaved similarly and no notable distinction was observed. For both the difference and the sum of the areas, there was no monotonic relationship to the makespan. One interesting pattern observed is that the individual areas themselves tend to form symmetric clusters. For example, in the case of $n=4$⁠, either all areas are equal, or there exists a 3-to-1 ratio, or a 2-to-2 ratio. Similar ratios were observed for higher values of $n$⁠. Although the ratios change depending on the layer geometry, the areas were never observed to be randomly distributed (Fig. 12). The clustering pattern observed and the analysis of the areas are also seen for the other layer geometries considered. We posit that, in general, there is no simple partitioning geometry metric that relates to the makespan, implying that our original hypothesis regarding equality of areas is false. Therefore, minimizing the layer printing time does not necessarily yield a partition with sub-regions of equal area. Note that these conclusions are based on the partitioning strategy for SZ*, which uses $2n$ partitions and buffer regions, and may not necessarily hold for other families of partitions (Sec. 6.1).

To physically validate SZ* we consider a scenario with four printers, where the part of interest is a solid disk of radius 100 mm. Since all the layers of the part are equal, and our framework uses layer synchronization, we only print one layer of the disk, which has a height of 0.45 mm. For printing, we used four SCARA printers developed by AMBOTS in a symmetric arrangement, sharing a square printing area with 600 mm $×$ 600 mm dimensions (Fig. 13). The choice of a disk, as opposed to a more complex geometry like the Stanford bunny cross section, was made to simplify the visualization of layer partitioning and the printing sequence for demonstration purposes. However, our framework is fully capable of printing parts with complex geometries, as demonstrated by the numerical results.

The printers are equipped with a Duet controller board running the RepRap firmware. They are connected to a computer (control hub) through a local WiFi network. At the higher level, we have a python script running on the control hub that oversees the printing process by following the wait-and-go communication protocol (Sec. 4.2). However, any programming language can be used as long as it supports sending HTTP requests. No external sensors besides the ones necessary for printer calibration, such as the BLTouch, were used.

In this manufacturing scenario, we ran SZ* with four printers to determine the optimal printer positions, and their corresponding partitioning and schedule, considering the disk layer geometry. However, the optimal printer positions found were not feasible, due to the gridlike nature of the experimental setup, forcing the printers to be placed in specific locations (Fig. 13). However, we observed that it was possible to use the tessellation given by SZ*, while keeping the printers in their feasible locations. This is because the feasible locations also satisfied the partitioning constraints (Sec. 4.1.1) for the optimal tessellation, allowing us to map the optimized printer positions to feasible ones without altering the partitioning or schedule.

Since the printer locations used in the optimization process differ from those used in actual printing, we call them VRPs. These positions act as Voronoi sites in the optimizer but do not correspond to physical printer locations. Using VRPs is particularly beneficial in highly-constrained manufacturing scenarios, such as this one, as otherwise the design space would shrink significantly, restricting the types of partitions that SZ* can find. While the mapping between VRPs and feasible printer locations is trivial in this scenario, we recognize that this may not always be the case, warranting further exploration of this concept (Sec. 6.2).

For the purposes of physical validation, we partitioned the disk layer in two different ways. The first partition uses Voronoi sites that are co-located with the only feasible printer locations available, thus lying in a square, which is a cyclic polygon, producing a partition with a degree-4 vertex. When creating the second partition, we leverage the concept of VRPs and allow the GA optimizer to look for other solutions by decoupling the Voronoi sites from the printer locations. The resulting partition has a single degree-3 vertex, making it three-colorable. (Fig. 14). The final printing schedule was computed using the heuristic-based approach in which the buffer regions $Bi$ are printed first, followed by the safe regions $Ai*$⁠, and sent to the control hub for printing with the wait-and-go communication protocol. To facilitate visualization, we used different PLA filament colors for each printer: white, matte red, navy blue, and black. Regions $A1*$ and $B1$ correspond to the bottom left corner of the part and the numbering follows a counterclokwise order, thus, regions $A2*$ and $B2$ correspond to the bottom right corner, and so on (Figs. 15(a)–15(d)). The time taken to print each sub-region and their respective areas is shown in Table 1. Time-lapses of both cases are also available.²

Due to the presence of degree-4 nodes in the symmetric partition, we see that it lacks potential for parallelization when printing buffer regions, which ultimately impacts the layer printing time to be 2135 s or roughly 35 min. On the other hand, the degree-3 partition, which was created using our methodology, allows buffer regions $B2$ (matte red) and $B4$ (black) to be printed simultaneously, yielding a printing time of 1661 s or approximately 28 min. Considering that the single-printer layer printing time $Ts$ is about 3000 s (50 min), we see that the symmetric partition results in a printing time reduction of 28.83%, compared to 44.63% when using our degree-3 partition. Consequently, our degree-3 partition has a 22.20% reduction in printing time compared to the printing time of the symmetric case. It is important to note that the degree-3 partition has no reachability issues, which is enforced in our methodology by changing the form of the partitioning constraints with VRPs, thus making it collision-free, time-optimal, and applicable in real C3DP systems with robotic arms.

SZ* is a promising framework for C3DP and demonstrates good performance in physical systems, however, there are several points that require further investigation. The proposed methodology currently limits the number of sub-regions created in the partition to be at most twice the number of printers. This limitation impedes SZ* to look for other families of partitions which may be of interest to C3DP. Additionally, a more general approach to using VRPs could be construed, as SZ* currently requires the printer positions to be fixed. We discuss some of these issues in the following subsections. One important topic that we do not address here is that of the mechanical properties of cooperatively-printed parts. However, some of our previous research has partially investigated this area for sloped-chunking strategies [11] and topologically-interlocked parts [24].

By design, SZ* partitions space into at most $2n$ sub-regions, where $n$ is the number of printers used in the manufacturing process. This is always the case since the original $n$ partitions, created by the Voronoi tessellation, can generate at most an additional $n$ sub-regions during the buffer region creation process (Sec. 3.3). While it is possible to either modify SZ*, or devise new algorithms to consider other families of partitions, a number of challenges in the implementation must be considered. First, increasing the number of partitions also increases the complexity of the manufacturing process, as the printers have additional tasks that must be carefully coordinated. Second, more partitions increases the order of the partition graph $GP$⁠, which in turn increase the complexity of the graph scheduling problem. It is well-known that, due to the combinatorial nature of such graph scheduling problems, they do not scale well in the number of nodes. Finally, while our approach intelligently subdivides the original Voronoi cell into two, it is not trivial to devise a general partitioning and assigning rule that has the convenient properties of $SZ*$⁠, such as minimizing chances of collision and reducing issues with reachability. Challenges aside, we devote the remainder of this section to the analysis of interesting families of partitions, and their properties, which are of interest and may be explored in future research.

While using partitions that satisfy (Eq. (14)) could be an interesting approach for C3DP systems, one must be careful to create graphs that minimize the chances of collision. As illustrated in the given ideal graph (Fig. 16(a)), the likelihood of collision is high and significant control effort would be required to guarantee safety, regardless of printer placement. Furthermore, without a kinematically redundant system, which is not guaranteed, it would be challenging to avoid collisions between the printers working on $A1$ and $A3$ and the printer working on $A7$⁠. Hence, while this example may be well-suited for drones or mobile robots, it may lead to complications during implementation in systems with robotic manipulators. However, there may be other types of ideal graphs that have more favorable characteristics for robotic manipulators.

Some layer geometries may have holes in them, i.e., they may not be genus-0 surfaces. In such cases, if we consider disjoint sub-regions to be separate sub-regions, relaxing the assumption made in Sec. 4.2, then the underlying partition graph may change, potentially modifying its chromatic number, and allowing for a better printing schedule. Consider the case of a disk layer geometry, which is a genus-1 surface. Following the Voronoi tessellation, if we treat the buffer regions as separate we end up with 12 total sub-regions. In this case, buffer region $B1$ becomes $B1$ and $B2$⁠, buffer region $B2$ becomes $B3$ and $B4$⁠, and so on (Fig. 16(b)). This creates an even cycle graph $C12$ which is technically two-colorable. However, following the rule for ideal graphs (Eq. (14)), and given that the original partition considered four printers, we employ a three-coloring scheme, which under the same area assumptions (Sec. 6.1.1), yields the best possible schedule. The key insight here is that such a partition would not have been possible if the layer geometry was genus-0, and we did not treat the original buffer regions as disjoint. Finally, while we presented two different partitioning strategies that may be interesting for C3DP applications, we do not discard that there may be others, as this enumeration was not exhaustive, but we leave this to future work.

As shown in the physical validation of SZ* (Sec. 5.2), it can be useful to decouple the Voronoi sites from the physical printer locations by leveraging the concept of VRPs in highly-constrained manufacturing scenarios. However, the mapping between the optimal printing positions and any set of feasible printing positions is not trivial. This is due to two main reasons. First, the mapping must ensure that all the partitioning constraints are satisfied, but those constraints were originally defined with respect to the VRPs, so additional work is needed. Second, given that this mapping is not a well-defined function, meaning that a single set of VRPs can map to multiple sets of feasible printer locations, we must explicitly define how the algorithm should behave. In this section, we provide ideas for addressing these issues either assuming that there are user-defined fixed positions or not.

If we assume that the user determines the printer positions before the optimization process, then the mapping becomes easy, only leaving us with the constraint satisfaction problem. Since VRPs are not co-located with the printer positions, they have different partitioning constraints from the ones in the original framework (Sec. 4.1.1). Thus, if we have a user-defined set of fixed printer locations ${qi}$⁠, we can use the VRPs in the optimization framework by simply changing the reachability constraint to center the circular workspace $Cq,i$ around the fixed points ${qi}$⁠. In this case, it is expected that chosen locations already satisfy the separation and the configuration space constraints. With this new set of constraints, we can leverage the VRPs to search for a wider range of solutions (Fig. 17).

On the other hand, if we assume that the user does not fix the location of the printers, then the mapping must be first defined. This is because the map transforming a set of VRPs to a feasible set of physical robot positions is not a well-defined function, i.e., there could be multiple feasible printer locations associated with the optimal set of VRPs found. The simplest way to do this mapping is to choose the closest available printer location in the configuration space to that VRP. While this may work for the majority of cases, there can be scenarios in which doing so would result in a violation of the separation constraint, meaning that two printers are placed too closed to each other. Thus, we see that there is no simple rule that can define this mapping while guaranteeing a feasible set of physical locations. In fact, this mapping can be formally framed as a constraint satisfaction problem (CSP), in which we are interested in finding a feasible set of physical printer locations. In the framework of the CSP, the set of variables is the physical printer locations, the domain is given by the workspace of the manufacturing system, and the set of constraints is the partitioning constraints. While the CSP formulation can be efficiently solved using state-of-the-art algorithms, it would induce a higher computational burden in the SZ* framework. Therefore, a more detailed study regarding when VRPs should be used, given the level of freedom in printer positioning in a specific manufacturing system, is required.

SZ* is inherently a layer-by-layer framework, as layer geometry may vary between layers. Consequently, for a 3D-printed part with changing layer geometry, each printer may have a different optimal location depending on the layer number. While this poses no algorithmic challenge to SZ*, as the current framework assumes that the printers are mobile and can be repositioned before starting each new layer, we recognize that frequent movement may affect printer calibration procedures, which some users may wish to avoid. Therefore, this section explores how SZ* can be adapted to minimize inter-layer movement while preserving its optimal properties.

A key observation is that, for most 3D-printed parts, the layer-by-layer transitions are relatively smooth, meaning discontinuous changes in layer geometry are uncommon. Consequently, the optimal tessellation is likely to change only slightly between layers, resulting in minimal adjustments to printer positions. By leveraging the concept of VRPs (Sec. 6.2) and decoupling the Voronoi sites from printer positions, the same printer positions can be maintained across multiple optimal tessellations, reducing the need for repositioning between layers. Between layers where discontinuities occur, movement may be necessary. However, SZ* can be adapted to reuse printer positions as much as possible, ensuring that they satisfy reachability constraints for each optimal tessellation.

If printer movement between layers is not feasible, SZ* would need to be adapted to perform tessellation at the 3D level. This can be achieved by retaining the existing optimization framework while modifying the evaluation function to assess the schedule across all layers for a single tessellation. In this approach, an estimate of the total print makespan would guide tessellation decisions. However, this method would likely yield more suboptimal solutions, compared to the layer-by-layer approach, due to the lack of flexibility in printer repositioning. Additionally, it would be computationally expensive, as each solution would need to be evaluated across all layers of the printed part. Therefore, it is not trivial to design an efficient algorithm to find optimal tessellation at the 3D level for C3DP systems.

In this article, we developed and physically validated SafeZone*, a collision-free and time-optimal C3DP framework. The results show that the framework is viable, scalable, and has a lower printing time than traditional C3DP methods available in the literature. The key conceptual contribution of our algorithm is the co-solving of the partitioning and scheduling problems through an integrated geometric-topological analysis. However, based on the discussion topics presented, there are many interesting research directions that could further enhance our understanding of C3DP and the theory on which SafeZone* is grounded.

This work was partially supported by the National Science Foundation (NSF) Award #2048182 (Engineering Design and Systems Engineering Program) and the Army Research Lab (ARL) via Grant No. W911NF-24-2-0007. Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF or the ARL.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Additionally, we have also shown one example for the convergence graph of the fitness of our GA optimizer (Fig. 18). A similar pattern was observed for all numerical results, regardless of geometry and number of printers. Note that the theoretical best show in the figure is coincidentally equal to one, as the area of the z-beam is four units and there are four printers. The signs are flipped due to the internal implementation attempting to maximize the fitness function.