A Decomposition and a Scheduling Framework for Enabling Aerial 3D Printing

Aerial 3D printing is a pioneering technology yet in its conceptual stage that combines frontiers of 3D printing and Unmanned aerial vehicles (UAVs) aiming to construct large-scale structures in remote and hard-to-reach locations autonomously. The envisioned technology will enable a paradigm shift in the construction and manufacturing industries by utilizing UAVs as precision ﬂying construction workers. However, the limited payload-carrying capacity of the UAVs, along with the intricate dexterity required for manipulation and planning, imposes a formidable barrier to overcome. Aiming to surpass these issues, a novel aerial decomposition-based and scheduling 3D printing framework is presented in this article, which considers a near-optimal decomposition of the original 3D shape of the model into smaller, more manageable sub-parts called chunks. This is achieved by searching for planar cuts based on a heuristic function incorporating necessary constraints associated with the interconnectivity between subparts, while avoiding any possibility of collision between the UAV’s extruder and generated chunks. Additionally, an autonomous task allocation framework is presented, which determines a priority-based sequence to assign each printable chunk to a UAV for manufacturing. The efﬁcacy of the proposed framework is demonstrated using the physics-based Gazebo simulation engine, where various primitive CAD-based aerial 3D constructions are established, accounting for the nonlinear UAVs dynamics, associated motion planning and reactive navigation through Model predictive control.


Introduction
Advancements in 3D printing technology are anticipated to significantly transform the perception of the operational procedures for the manufacturing and construction industries [1,2], while it has already gained widespread use in various fields, such as healthcare [3], agriculture [4], construction [5], automotive and aerospace industries [6].Recent initiatives explore the feasibility of using 3D printing to construct full-scale infrastructures, such as emergency shelters, in disaster-stricken areas [7,8], remote and challenging environments exposed to extreme weather conditions, unstructured terrains, and distant locations.However, constructing such full-scale infrastructure using the conventional 3D printing approach demands large-scale manufacturing machines, making it difficult to realize in hard-to-reach and remote areas [9].In this context, the use of aerial 3D printing recently introduced in [10] is a viable alternative.In principle, the concept of collaborative aerial manufacturing is inspired by nature, where builders like wasps and bees construct nests in groups.Considering the significant recent advancements in autonomous UAV technologies that can enable precise motion planning, navigation, and task execution in a group for material depositing using UAVs.A feasibility study for UAV-based masonry construction of real-scale structures is presented in [11,12].Hence, the use of aerial robots to deposit materials for 3D printing is an innovative and still a relatively unexplored field, emerging as a promising direction for autonomous scalable manufacturing.To the best of the authors' knowledge, only a few research are reported, most of which are still at the primitive conceptual stage.This article presents a more general and flexible aerial 3D printing framework aimed at seamlessly integrating the simultaneous operation of multiple UAVs in a collaborative manner.

Related Work
The current approach to aerial 3D printing is limited by its primitive design strategy, which follows a 2.5 dimensional formalism commonly used in additive manufacturing.This restricts the motion of the printing nozzle to a horizontal (X − Y ) plane, with limited variability along the Z axis [10,12].An extensive review recently reported in [13] has identified that the limited manoeuvrability results in poor bonding strength between layers, staircase-like effects on surface quality and external support are often necessary.A more advanced printing strategy involving flexible multiple degrees of freedom (DoF) motion planning is recommended to address these issues.
Aerial robots possess agile maneuverability with six degrees of freedom (DoF) motion, which can expedite the construction of complex shapes with minimal external support.This is an essential aspect of 3D printing in general, as support structures are typically required to prevent the object from deforming or collapsing due to gravity [14].In order to minimize the requirement of support structures, various optimal slicing mechanisms [15,16] are being explored that distribute the object to be printed into smaller sub-parts enabling multi-axis layered manufacturing.However, in the context of aerial additive manufacturing, it is necessary to design the sub-parts in a manner that incorporates printability constraints to accommodate the extruder's shape, while also ensuring seamless connectivity between them.In a collaborative manufacturing scenario, a scheduling mechanism [17] establishes the order in which the sub-parts are to be printed in a distributed manner to build the object from the bottom up.To the best of our knowledge, no advanced general framework utilizing multiaxis-layered manufacturing in the context of aerial 3D printing has been described in the existing literature.

Contributions
The article introduces a novel collaborative aerial 3D printing framework that facilitates seamless integration of distributed aerial 3D printing with multiple UAVs in a multi-axis layered manufacturing environment.The primary contributions of the article are twofold.Firstly, it presents a generic chunk formation mechanism that effectively decomposes any arbitrary geometric shape to be printed into multiple sub-parts suitable for collaborative aerial 3D printing that extends the work done in [10].The formation of the chunks is executed based on multiple planar cuts, which incorporate printability constraints while adhering to the extruder's shape, thus enabling smooth connectivity between neighboring sub-parts as an addition to [15].Secondly, it establishes an innovative task scheduling mechanism, which generates a workable sequence and priority for parallel execution of distributed aerial 3D printing in a coordinated manner by deploying multiple UAVs extending the study presented in [17] towards aerial platforms.Finally, the proposed aerial 3D printing mechanism is evaluated in a Gazebo-based simulation framework incorporating multiple UAVs to print a complex, large-scale 3D structure.

Organization of the Article
The rest of the article is structured as follows.In Section 2, an overview of the problem formulation is presented, while in Section 3, the novel mechanism for the construction of the feasible chunks is proposed.The execution of aerial printing with UAVs and associated trajectory generation and a path-following controller is presented in Section 4. The overall efficiency of the proposed framework is validated with realistic simulation results in Section 5, while the article is concluded in Section 6.

Methodology
The first step towards the realization of the novel aerial 3D Printing is a mesh decomposition process in which the original mesh given by the user is decomposed into smaller sub-parts called chunks.This is done by near-optimally selecting a set of planes, which will act as planar cuts for the mesh.This selection is executed through a search-based approach in which possible cutting planes are sampled and evaluated.The best possible combinations of planes proceed to the next iteration, which will be extended with another cut.The evaluation is implemented by introducing a heuristic function aiming to reduce the dispersion between the volumes of the chunks and increase the parallelability of chunk printing tasks.A Binary Space Partitioning Tree (BSP) [18] is constructed to keep track of the planar cuts and the resulting structure of the produced chunks.After the search phase, the BSP tree serves as a printing task priority sequencer, as well due to the inherited parent-child relations of its structure, which imposes the dependencies between the chunks.The priority sequence during printing must be followed to successfully recreate the expected original mesh.Finally, every chunk is assigned to each UAV in a reactive fashion, while each chunk is sliced through a commercial open-source slicer and transformed into a trajectory which will be later executed by the assigned UAV.An overview of the proposed approach can be seen in Fig. 1.

Iterative Search-based Chunks Generation
The chunk generation process is carried out by intersecting the mesh with a set of planes, fragmenting it into smaller  printable sub-parts.The cutting planes are denoted by j i ∈ R 3 , defined with the surface-normal vector n i ∈ R 3 , perpendicular to the plane, originated at p j ∈ R 3 and lies on respective cutting surfaces.The proposed chunk generation framework utilizes an iterative beam search approach to generate a BSP tree that decomposes the mesh into multiple parts guided by a heuristic function that ensures uniformity of chunks in volumetric distribution, while accounting for constraints to assess the feasibility and printability of the generated chunks with respect to collaborative aerial 3D printing.A comprehensive description of each of the components culminating in their association with the proposed iterative search-based chunk generation framework is presented next.

Organization of Chunks with BSP Tree
The executing planes fragment the overall mesh into multiple chunks, which must be organized systematically, while preserving their associations and connectivity to mechanize the printing with an orderly sequence.A Binary Space Partitioning (BSP) Tree T is introduced in the process of the chunk generation to keep track of both the sequence of cuts and the resulting chunks.Each node of the BSP tree preserves information regarding the plane with which a cut is being executed and the associated mesh to which the cut is being applied.The original mesh is located at the root node, while all the generated chunks C i are located at the leaves of T .Whenever a new planar cut i is imposed, all the chunks C j affected by it are found recursively and are extended further into children nodes by executing the planar cut.The intermediate children nodes are denoted with chunks C − j,i and C + j,i , which are essentially the result of the cut between i and C j .Out of those two, the C + j,i lies on the half-space that is along with the direction of the normal n i and is called a positive chunk, while the other C − j,i is called a negative chunk.For this approach, the convention of placing the negative chunk C − j,i as the left child of the parent node is followed in this article.It is important to note that the topology of the BSP tree may change depending on the order in which the planar cuts are executed.However, this will not affect the resulting sets of chunks, which will remain the same.

Beam Search Based Planar Cuts
The primary objective of the search module is to find the optimal combination of planes j i that will serve as planar cuts to the original mesh and decompose it into a set of chunks with desired characteristics to ensure collaborative printability with multiple UAVs.In order to achieve this goal, an iterative beam-search-based approach [15] is considered in this article.The candidate planes for the search-based chunk generation are uniformly sampled in the following approach.The possible normal vectors are initially sampled in a collection denoted as N .The vectors are sampled along the spherical space S ∈ R 3 = {r , φ sample , θ sample } with respect to the constraints r = 1, φ sample ∈ [−φ max sample , φ max sample ] and θ sample ∈ [0, 2π).The number of sampled normals in each dimension is determined by the constants M φ and M θ for φ and θ angle correspondingly.An instance of this sampling is shown in Fig. 2.
For each normal vector n i ∈ N a family of planes F = { j i }, j = 1, . . ., K is calculated where j i = ( n i , p j ).All the origin vectors p j are sampled along the direction defined by the normal vector n i by taking the minimum and maximum of the collection M = { j ∈ R : j = pr oj(v k , n i )} which contains all the projection of the vertices v i of the mesh to the normal vector n i .An interpolation between j min and j max is executed by fixed step δ ∈ R , resulting to the set P = { p j } where p j = ( j min + l • δ) • n i , l = 0, 1, . . ., (l max − l min )/δ.A family F of plans of planes for a given normal vector n i is shown in Fig. 3.
All the possible planar cuts j i are executed to the given tree T that needs to be extended and the newly generated trees after the extension are returned into a list.The whole process regarding the extension of the tree along with the evaluation of the new trees is shown in Algorithm 1.
This search is being executed iteratively, and in each iteration, all the trees that do not satisfy the termination condition are appended in a set T containing the most promising ones up to that point based on the heuristic function.For each tree in this set, an investigation is carried out for further extension with new planes.Every tree is extended by all the candidate planes sampled, but only the best W inner are chosen to be inserted into the list with all the extended trees of this search iteration.This task is executed iteratively for all the nonfeasible trees and after the list with the newly extended trees

Heuristic
As presented in Algorithm 1, each intermediate BSP tree containing the chunks is evaluated using a heuristic function, which initially aims to spread uniformity among the volumes of the generated chunks.This leads to shorter waiting times between each print since all the printing tasks take approx-imately the same duration to be fulfilled.Additionally, the uniformity in volume results in an evenly distributed arrangement of the chunks, with their geometric centres spaced apart sufficiently.This enables the feasibility of collaborative printing by allowing for collision-free operation of UAVs in a spread-out workspace.The volumetric uniformity is integrated into the heuristic function by penalizing the dispersion c v of the volumes of the generated chunks.Assuming that the whole volume of the original mesh is notated by V and the volume of the i-th chunk by V i then V = N i=1 V i , where N is the number of chunks.Towards calculating the terms c v , the mean μ = N i=1 V i /N along with the standard deviation 2 /N of the volumes are extracted and combined as follows By penalizing this value the search is led to a uniform distribution of chunk volumes.
During the search phase, another heuristic value that is considered is the total number of seed chunks present in a BSP tree.Following the notation described in [19], a seed chunk is defined as a chunk where all its faces resulting from the planar cuts are positive, as illustrated in Fig. 4.These types of chunks are desired to be as more as possible in the final tree since their only requirement is to be printed on the ground or on top of a flat chunk and then all of their faces act as supports for the next chunks in the sequence.This property leads to an increase in the parallelability of printing tasks since no further dependencies are imposed.
In order to identify if a chunk C k is a seed or not, a set C = { i , i = 1, 2, . . ., M} is gathered where i are all the planes that resulted in cuts affecting this chunk and are co-planar with the chunk.Having these planes leads to the next step of identifying whether the chunk is negative or positive to the corresponding plane i .This is extracted by taking a random (non-lying on the investigated face) vertice v k j ∈ R 3 of the chunk's C k mesh and calculating its projection pr oj n i ( v k j ) to the plane normal n i as follows: where p i is the origin of the plane i .The final classification is evaluated by the value s k ∈ R where So, a chunk C k is a seed whenever s k = 1, and this property is rewarded by adding the term s k in the heuristic function.In addition to the previous, a further extension is also introduced where the number of positive faces of a chunk C k is rewarded, provided that this chunk is a seed.A positive face of a seed chunk C k is defined as the face that is coplanar with a plane i and is characterized by the term f k+ i, j .
Consequently, for calculating the number of these faces of a given chunk, the term p k f aces is introduced as a reward in the heuristic.
The two aforementioned heuristics can be incorporated into a common heuristic function g(C k ) ∈ R which adds them scaled by some arbitrary weights G part and G f aces that are tuned based on the objectives of the application Finally, the total heuristic function h(T ) for the evaluation of the tree T is given by where M is the number of chunks of the specific BSP tree T .
It must be noted that the sum of the functions g(C k ) referring to each chunk are being negated due to the fact that these values need to be rewarded and not penalized.

Constraints
During the sampling phase of the normal vectors, the sampling surface S is bounded due to constraints derived from the connectivity between the chunks and the mitigation of collisions between the extruder and the previously printed chunks.

Connectivity Constraint
Following the mesh decomposition phase, the generated chunks need to be printed by depositing material on top of previously printed chunks.Hence, it is necessary to ensure sufficient printable overlapping areas between the faces of the adjacent surfaces to establish a robust connection to stick together effectively.This requirement is set as the connectivity constraint for the search algorithm and is incorporated during the sampling of the normals n i .Specifically, the angle φ defined by the horizontal axis and the normal vector n i is directly tied to the printable overlapping area, as depicted in Fig. 5. Hence, an upper-bound constraint φ conn max is imposed to restrict steepness.This bound is set to 45 • , which is usually a rule of thumb for 3D printing scenarios [20] in order to ensure 50% overlap between layers.It is worth mentioning that this limit may vary depending on the printing material and the extruder's freedom of motion so that more steep angles could be reached.

Extruder Collision Constraint
In addition to the connectivity constraint, the chunk generation process also considers the possibility of collision between the extruder nozzle/head and previously printed chunks.To prevent such collisions, the maximum allowable cutting slope is determined based on the geometry of the extruder nozzle and head.A representative visualization is depicted in Fig. 6.Considering that a generic shape of the extruder head can be approximated with a bounding rectangle, the maximum allowable cutting slope is determined as where, h represents the height of the nozzle, and l denotes the length of the extruder nozzle to the outer surface of the extruder's head.
The aforementioned connectivity and collision constraints are integrated into a common maximum angle φ max sample , supplementary to the minimum of the φ conn max and φ extr max since it is referring to spherical coordinates (with respect to the vertical z-axis).
The constraint on φ max sample is imposed while determining the sample space for the planar cuts in the beam search algorithm.Similarly, the sorting approach is used for the chunk volumes, which are described by the set C = [c 1 , c 2 , . . ., c k ].Before initiating the search algorithm, a primary feasibility check is required to ensure that the quantity of material carried by the UAVs is larger than the volume of the entire original mesh, represented as

Feasibility Condition
where, V denote the volume of the original mesh.The feasibility condition is necessary to initiate the chunk generation.

Termination Condition
The formation of a BSP tree is considered to be terminated when all the chunks generated at its leaf node can be printed by at least one of the available UAVs.A UAV can print multiple chunks if its carrying capacity allows the required material volume.However, the reverse scenario, i.e., multiple UAVs printing one chunk, is considered to be avoided to eliminate the associated complexity of task planning and coordination of multiple UAVs in a shared workspace.Hence, a tree is considered feasible to be printed when a non-empty set exists and is defined as: The sets C, D and S are updated recursively after the construction of each chunk.

Priority Based Task Allocator
The BSP tree is introduced to organize the planar cuts during the search phase.Nonetheless, allocating the generated chunks into the leaves of the BSP tree retains the advantage of establishing an orderly sequence in which the chunks can be seamlessly manufactured to create the desired overall mesh.Limiting the slope of the executing plane by imposed constraints φ max ensures its normal vectors n i retain a positive component in the vertical direction (i.e., along z axis).As a result, the two resulting chunks produced after the execution of the cut will have a dependency between them regarding their manufacturability and structure.Specifically, the chunk that lies on the half-space in the direction of the normal vector (referred as the positive chunk denoted with C + j,i ) will have all its vertices and faces over the plane i of the cut in contrast with the negative chunk C − j,i that will have them below it, resulting to the fact that C + j,i being on top of C − j,i .Considering the characteristic of additive manufacturing, which is typically performed by constructing from the bottom to top layered fashion, leads to an imposed priority on printing the negative chunks over the positive ones to enable seamless integration of the overall topology cohesively.In view of that, the left child nodes of the BSP tree T and all of their subordinates are prioritized over the right one for each node.Expanding this property recursively to the whole tree leads to the priority given by a Depth-First-Search (DFS) pre-order traversal and keeping only the leaf nodes since the chunks are located there.Following this pattern, although the correct sequence that the chunks need to be manufactured is given, this sequence depends on the order of the planes selected during the search phase.Towards eliminating this unwanted result and following the overall approach of manufacturing chunks in a uniform way from the bottom to top, the set P of the planar cuts contained in the resultant BSP tree T is extracted and is reordered based on the z-coordinate of each origin point p, starting from lower to higher.This leads to a new set P or d and all the cuts are executed again by following Fig. 7 Chunks manufacturing sequence before and after reordering based on the vertical component of origin vector p j of each planar cut j i the newly ordered sequence.An illustration of this ordering can be seen in Fig. 7.
The assigned priority enforces a manageable and seamless deposition process, devoid of any infringement upon the structural interdependencies between the chunks.Additionally, the ordered printing of chunks is arranged such that each subsequent piece is supported by the preceding one, while maintaining the potential for parallelization.

Chunk Printing Procedure
Once the planar cuts and chunk generation have been completed, the chunks are allocated to individual UAVs according to the sequence P or d dictated by the topology of the BSP tree and the procedure mentioned in the previous section.Whenever a chunk is assigned to a UAV, a series of tasks are carried out in sequence to produce a trajectory for the UAV to follow.This process is depicted in Fig. 8.

Formation of Extruder Path using Slicer
In this study, each chunk is preserved as an individual mesh, consisting of triangles connecting the vertices of the mesh.To determine the path necessary for successful manufacturing of each chunk using the 3D extruder, the open-source software 'Cura' slicer [21] is used, which is typically interfaced with commercial 3D printers.In this context, the configuration is tailored with precise parameters for aerial 3D printing, including the UAV and extruder dimensions and properties, while integrating variables related to material deposition, such as the height of each layer, the diameter of the extruder, and the infill percentage between the walls of the chunk, which are adjusted to suit the requirements.Once the appropriate parameters are set, the slicer generates a g code path in linear segments that define the movements of the extruder.Typically, a specific material extrusion rate and temperature are also determined for each movement.However, the present work primarily focuses on the motion planning of the UAVs to enable accurate tracking of the directed path generated by the slicer, while the material extraction rate and temperature are not given significant importance for the simulation-based demonstration of the proposed framework.The utilization of the 'Cura' slicer provides a comprehensive and customizable Fig. 9 UAV-extruder setup.A gimbal of length l g is attached below the UAV, and at the other end of it, an extruder of length l ex is mounted approach, enabling successful manufacturing with UAVsassisted aerial 3D printing.

G-Code To UAV Trajectory Transformation
During the process of aerial 3D printing, the UAV controller cannot handle the "G-Code" format, which is generated by the 'Cura' slicer.Thus, a post-processing task is required to discard unnecessary information and store the path in a simpler format that can be fed to the UAV.This is achieved by assuming that the UAV moves at a constant velocity and that all the segments extracted by the G-code are converted into equivalent waypoints that UAV needs to pass through.The resulting trajectory P ex is in a format that UAV controller can be interfaced with.However, it should be noted that this trajectory refers to the path of the extruder, not the path of the UAV.Therefore, a transformation is needed to convert the frame of the extruder tip to the body frame of the UAV.In this context, a linear gimbal with a length of l g and an extruder with a length of l ex are used for the printing action, with the extruder tip pointing downwards vertically towards the ground.An illustration of the setup is presented in Fig. 9.
Finally, the UAV trajectory P U AV is obtained by transforming the extruder path to the UAV body frame, with the where, cθ g sφ g sθ g cφ g sθ g 0 0 cφ g −sφ g 0 −sθ g sφ g cθ g cφ g cθ g −l ex 0 0 0 1 Here, φ g and θ g , denote the gimbal angles around x-axis and y-axis of the gimbal frame G respectively, and sine and cosine of angles are abbreviated with s and c.

NMPC-based Reactive Motion Planning
In order to execute the aerial 3D printing, each UAV needs to precisely traverse through the desired path P U AV , while depositing material uniformly.Towards this direction, the Nonlinear model predictive control (NMPC) based optimal path following controller [22] for accurate tracking of G code generated trajectory is adopted in this article.
The UAV dynamics is described using two coordinate frames, namely the world frame (W = {X W , Y W , Z W }) ∈ R 3 and the body-fixed frame The body-fixed frame is attached to the UAV's centre of mass, while the inertial frame is assumed to be attached at a point on the ground, as shown in Fig. 10.
The UAV's motion is governed by a set of non-linear dynamic equations presented as  the world frame, A x , A y , A z ∈ R represent linear damping terms, and g ∈ R the Earth gravitational acceleration.The attitude dynamics is approximated with a first-order system representation, characterized by the time constants τ φ and τ θ .The overall motion control system incorporates a lowlevel attitude controller with constant gains K φ and K θ that enforces to follow the reference orientation command φ re f ∈ R and θ re f ∈ R. Essentially, φ re f , and θ re f combined with the thrust magnitude T > 0 ∈ R constitute the control input, which is designed using NMPC-based high-level position controller to enable reactive path following accurately.
The UAV dynamics presented in Eq. 13 is discretized using the forward Euler method with sampling time δ t ∈ R + and expressed in a compact mathematical notation given as where, The discretized model is utilized as the receding horizon prediction model for the NMPC framework with a prediction horizon denoted as N .At a given time instant k, the predicted states and associated control input over the entire prediction horizon are presented as The objective of the control design is to find an optimal sequence of u k in each time instant that minimizes the cost function expressed as input smoothness cost (15) where The state cost penalizes the deviation of the predicted future position at each time step of the horizon with the commanded reference state.The input cost term enforces the controller to minimize deviation from reference control input u, and maintain as close as possible to hovering condition u re f = [g, 0, 0].The term input smoothness cost penalizes the deviations between the successive inputs to maintain smooth control action by minimizing the control rate.The cost function is minimized subjected to the system dynamics and input constraints u min ≤ u ≤ u max where u min , u max ∈ R 3 represents input bounds.In addition, the maximum bounds are also imposed on the rate of change of φ re f , θ re f , which are bounded to φ max , θ max and expressed as a set of equality constraint presented as Minimizing the control rate while enforcing its magnitude bound restricts aggressive maneuver and ensures sufficient smoothness on the path traversed by UAVs.This guarantees a uniform deposition of printing material throughout the trajectory, ensuring effective aerial 3D printing.The overall optimization problem in compact notation is presented as The solution of the optimization problem defined in Eq. 17, is obtained numerically using PANOC by employing a single-shooting approach [23] using the Optimization Engine (OpEn) platform [24].'OpEn' is an embedded nonconvex optimization framework that employs automatic differentiation to obtain the cost function gradient in CasADI [25].Compared to conventional approaches such as SQP and Interior-point methods, OpEn has a faster optimization time with higher accuracy [26].A representative numerical evaluation of the required computational time for obtaining the solution of the nonlinear optimization problem posed by the proposed NMPC-based navigation framework is depicted in Fig. 11.The computational times are recorded for obtaining the control command, consisting of an entire control window, at each time instant.In the present context, the computational time is evaluated in desktop computer with the processing speed of 2.1 Ghz, which can be fairly equivalent with the state of the art onboard computer for UAV, such as Intel ® NUC computer used for UAV as on board processing unit, with processing speed of 2.0 GHz.
The computational time remains consistently bounded, approximately within 30ms.Considering the practical implementation scenario for UAV control, typically the rate of control update operates around 20Hz, i.e., each control command computation should ideally be completed within the bound of 50ms.A similar NMPC-based waypoint tracking control including obstacle avoidance has been experimentally evaluated in [27].In view of that, the proposed MPCbased navigation framework allows a sufficient margin, as the computation time fall within the required timeframe.

Simulation Results and Discussion
The efficacy of the proposed collaborative aerial 3D printing framework is evaluated in multiple scenarios by constructing various primitive shapes with different UAV configurations using the chunk-based approach.The execution of the overall framework and the printing processes are assessed in the physics-based Gazebo simulation engine.A detailed description of simulation environments, execution steps, results and discussion are presented in the following sections.

Chunk Based Mesh Decomposition
In the present simulation study, CAD models of three different shapes are being considered for multiple UAV-based collaborative aerial printing: a rectangular parallelepiped, a hollow rectangle, and a hemispherical dome.In order to execute the iterative search-based chunk generation, the bounding angles for the planar cuts sampling space S are set to φ max sample = 45 • and θ ∈ [0, 2π) while the weights of the heuristic function are designed as follows: G disp = 200, G part = 10 and G f aces = 20.
Initially, a solid parallelepiped is considered for chunk generation, which is decomposed into nine chunks, as shown that facilitates the printing process and reduces its associated complexity.The partitioning is designed to ensure adherence to the prescribed sequence, avoiding overlaps and infeasible prints.Furthermore, the resulting decomposition includes chunks that can be printed simultaneously (e.g.0 − 1, 3 − 4 and 6 − 8), enhancing overall efficiency.The BSP tree containing the planar cuts, the colour-coded chunks, and the sequence extracted by its topology is shown in Fig. 12b.
Another mesh considered in this study is a hollow rectangle, as shown in Fig. 13a.The volume dispersion of the chunks is obtained to be c v = 0.674, there are five seed chunks and the total negative faces are nine.In this case, again, a planar slice parallel to the ground is chosen and splits the mesh into two layers.The rest of the planar cuts generated from the iterative search-based algorithm, result in five seed chunks identified with the index by (0, 1, 3, 8 and 11).The generated chunks comply with the prescribed constraints and facilitate a pattern whereby each component supports the succeeding ones during printing, thereby ensuring the full assembly of the mesh upon adherence to the designated sequence.Note that, a small portion of chunks (indexed as 6 and 7) are obtained to be relatively smaller in volume, which contradicts the uniformity in the volumetric dispersion-based heuristic strategy.Such an outcome is expected since the tuning parameters associated with the heuristic function in Eq. 7 are selected to prioritize the number of seed chunks over their volumetric uniformity.However, these relatively smaller chunks are considered to be acceptable considering all practical purposes.In order to mitigate this phenomenon, the unification of these smaller chunks by combining them back to their parent node could be considered, which can be deemed as a post-processing step for chunk generation.However, such a scenario is not accounted for in the present study.
In addition, the proposed framework is evaluated for printing a hemispherical dome approximated with lowerresolution triangles as shown in Fig. 14a.The volumetric dispersion is obtained as c v = 0.592.A total number of three seed parts is obtained with a total of seven faces.Despite obtaining comparatively fewer seeds than in the aforementioned cases, the selected set of planar cuts results in seven positive faces.Similar to the previous cases, a parallel toground cut is executed, which renders all the chunks above it independent from each other in terms of their lower face.This Fig. 14 Color-coded chunking Results with the corresponding planar cuts for a dome, stacked together forming the decomposed mesh(a).BSP tree T visualizing the cutting planes 1 , . . ., 4 in the interme-diate nodes and the chunks C 1 , . . ., C 11 on the leaf nodes.The printing sequence graph P or d is extracted at the bottom (b) decreases the possibility of having unwanted dependencies, which would slow down the printing process in a parallel implementation framework of collaborative printing.Moreover, upon moving closer to the upper part of the dome, it becomes apparent that chunk C 11 does not require a support structure, as it is already supported by the surrounding previously printed chunks.The output tree containing all the cuts and chunks, along with the sequence graph for the dome mesh, can be seen in Fig. 14b.
Another series of experiments are conducted aiming to investigate the scaling and the behaviour of the algorithm for different initial UAV configurations.Different numbers of UAVs with various volumes of material, which are shown in Table 1, are used and tested with three different meshes (solid parallelepiped, hollow rectangular and hemispherical dome).As depicted in Table 2, the number of chunks generated for each mesh increases with the number of available UAVs which is expected since as the number of UAVs increases, their carrying material volume decreases.Thus, the printing task must be distributed in more chunks so each one of them is sufficiently small to be printed completely by one UAV.
Overall, the scalability of the approach is demonstrated by the consistent increase in the number of chunks generated as the number of UAVs and mesh complexity increase, showing that the method can effectively handle larger and more complex meshes.

Printing Process Simulation
The proposed framework is tested in the simulation environment of Gazebo by using the "RotorS" [28] package.The "Carma" [29] UAV is chosen as the application platform and is integrated into the system without having any actuation on the extruder joints (i.e., φ = θ = 0 • ), and the extruder length is equal to 0.45 m.The UAV's mass is m = 3.59The deposition action of the material is simulated by attaching an odometry sensor at the end-effector of the arm and keeping track of the pose (position and orientation) throughout the whole execution.Whenever a printing action is commanded from the UAV, the extrusion of the material occurs and is visualized by plotting small spherical markers placed right below the tip of the extruder.The behaviour and the physics of the deposited material are not taken into account in this study and the assumption that the material will stick to the previously overlayed surface is made.The initial search and chunking for the printing mission are conducted by a centralized entity.A centralized mission planner is also implemented and is responsible for assigning chunks to each UAV to be printed according to the sequence and priority is given by the priority graph P or d .
The printing procedure for a layer of one chunk can be seen in Figs. 15 and 16, where the reference path of a chunk's layer for a UAV along with the measured position of the extruder are depicted in X-Y axes and their time histories.The tracking of the given path is smooth, and precise and is considered sufficient for the purpose of the aerial 3D printing application.Even though there is a small deviation of the order of a few centimeters from the commanded trajectory persists, the diameter of the deposited material is sufficient to compensate for it.The position error is observed to be larger Fig. 17 Final mesh composed of all the chunks generated and printed sequentially from the available UAVs during the start of the movement for each linear segment, which is an expected phenomenon accounted for due to the momentary swinging motion of the UAV, which results in the deviation of the extruder from the reference position because of its distance from the center of rotation.
For the full-scale evaluation of the framework, the dome mesh is chosen to be printed, which is a sphere with a diameter of 1.5 meters approximated with triangles.The chunk generation for this shape is carried out as presented in Fig. 14.Each chunk is assigned to the UAVs sequentially based on the generated priority sequence, a snapshot of the simulation with the final mesh can be seen in Fig. 17 The sequence followed along with the printed chunks is shown in Fig. 18 where each step of the process is shown until the fully deposited mesh is completed.
A video where the whole printing process is simulated and visualized can be found in https://youtu.be/AM4kUXydpjc.

Conclusions
The article introduces a flexible and efficient large-scale aerial 3D printing framework, incorporating multi-DoF capabilities.The proposed framework introduces a generic automated optimal mesh-decomposition mechanism that can decompose any given complex geometric shape into smaller sub-parts called chunks.The mesh-decomposition mechanism is based on the beam-search algorithm, which evaluates planar cuts with respect to constraints that ensure that no collision between the extruder and previously deposited chunks occurs while also ensuring that the overlapping area between chunks is sufficient for them to bond.The cuts, along with the chunks, are managed by a Binary Space Partitioning (BSP) tree that facilitates the arrangement of the chunks.The BSP tree also inherently includes the priority sequence that the chunks must be printed to fully assemble the original mesh.The priority sequence is extracted by a task allocation mechanism, which assigns each chunk to a UAV in a reactive fashion.This mechanism sets the basis for enabling distributed, parallel aerial 3D printing with multiple agents.A Model Predictive Controller (MPC) for tracking the linear segments generated from the slicer is implemented ensuring smooth motion planning and increased accuracy during the execution.The proposed framework is demonstrated through a simulation performed in Gazebo Simulator.

Fig. 1
Fig. 1 Block diagram of the planes search and chunking phase

Fig. 2 Algorithm 2
Fig. 2 Instance of sampled normal vectors along the surface S with respect to the maximum angle φ max

Fig. 4
Fig. 4 Seed Chunk Demonstration.Both cutting faces of chunk C 1 generated by planes 1 and 2 are positive so the chunk is characterized as a seed and gets rewarded during the heuristic search

Fig. 5
Fig. 5 Connectivity constraint between chunks C 1 and C 2 aiming to have sufficient area overlapping between them

Fig. 8
Fig. 8 Block diagram of the generation of a path P U AV for a given chunk C i 3 represents the velocity of UAV in world frame.The UAV roll and pitch angles are defined by φ, θ ∈ [−π, π] respectively.R(φ, θ ) ∈ SO(3) denotes the rotation matrix that transforms the acceleration due to thrust acting on the UAV body frame to its equivalent components into

Fig. 10
Fig. 10 Illustration of the world W and body B frames of reference

Fig. 11 Fig. 13
Fig.11 Coputational time (ms.) for each optimization execution of the non-linear model predictive controller

Fig. 15
Fig. 15 Desired (blue) and Measured (orange) Extruder Path for one layer of a chunk given to the UAV .

Fig. 16
Fig.16 Snapshot with the reference (blue) and measured (orange) path of the extruder tip during the printing process

Fig. 18
Fig. 18 Sequential simulation snapshots of color-coded chunks being printed

Table 2
Search executions for three different meshes (Dome, Hollow Box and Solid Box) given three different initial UAV configurations (4, 6 and 8 UAVs)