A two-nozzle cable-driven parallel robot for 3D printing building construction: path optimization and vibration analysis

This paper proposes a novel two-nozzle cable-driven parallel robot (CDPR) for 3D printing building construction. This system has two independently moving nozzles mounted on the existing printing head. Thus, the printing time can be reduced dramatically with this system as the travel path of the printing head can be reduced to almost half thanks to those two nozzles that print almost half of the printing contour. To fully take advantage of two nozzle structures effectively, the path of the printing head is optimized to secure the minimum travel length of both the printing head and two nozzles. The smoothness of the optimal path is secured by applying the non-uniform rational B-splines (NURBS). In addition, free vibration of the proposed CDPR printer’s structure is analyzed to improve the printing quality and help the control of the proposed CDPR by using a finite element formulation of cables of the proposed robot. From the analysis results, the critical design guide—stiffness of the proposed structure could be raised by increasing the minimum cable tension in the tension distribution optimization algorithm—was obtained as well.


Introduction
In recent years, 3D printing technology has been creating lots of applications in many industries, from mechanical manufacturing to medical applications. 3D printing has now expanded its application area even to the building construction. Due to the large-scale of the required work in construction, a conventional 3D printer based on the gantry or x-y-z Cartesian mechanism might not be suitable because they are heavy, has limitation in dimension expansion, and their construction (especially z-axis) cost are expensive. Thus, cable-driven parallel robots (CDPRs) have been proposed for 3D printing building construction [1][2][3][4][5] for their benefits in terms of large workspace, low cost, lightweight, flexible reconfiguration, and easy transportation [2]. Along with such benefits, CDPRs could be more competitive for 3D printing building constructions if construction time can be reduced more as the construction time is a key factor for the fulfillment of the construction deadline and the cost. Reducing the time for the 3D printing building construction can be achieved by several options such as using rapid-solidifying printing material and reducing the printing path length. In this research, we focus on reducing the printing path length approach. Also, as a viable option for path length reduction, we propose a CDPR with a multi-nozzle printing head. With the multi-nozzle printing head, path length of the printing head can be reduced dramatically as each nozzle covers the portion of the required printing path.
Considering the fact that the contour of the most buildings to be printed has a closed-shape, a number of nozzles are selected as two in this research. To the best knowledge of the authors, there has been no such two-nozzle CDPRs reported in the literature. Note that, even for the single nozzle CDPR, there have been only a few reported cases [1][2][3][4][5] for 3D printing building construction.
As can be imagined, applying extra nozzles to the existing single nozzle printing head to print certain portions of the building component would require a change in the whole process of 3D printing including the path planning. Especially for printing an asymmetric part using two nozzles, it requires independent motion control of each nozzle to print the object. Hence, there is a need for a specific path planning for the integrated system which has multiple nozzles working together [6].
Typical path planning aims to find a collision-free path from a starting point to a target point and to satisfy some criteria such as distance, time, or energy. Among them, the distance is one of the most commonly adopted criteria [7].
To find the optimal path for the two-nozzle CDPR, the following three constraints are considered: (1) The path of the printing head should be positioned in the middle of the contour of the parts to be printed to minimize the time required for the movement of both nozzles. This will minimize the difference between the path lengths of two nozzles; thus, both nozzles move the same length as much as possible which would make the motion control of two nozzles plain (2) The path length of the printing head should be minimum to minimize the printing time. (3) The path shall be as smooth as possible and at least satisfy the continuity so that it can be implemented physically. The first two constraints are related not only with time, but also with the energy of motion control system of the CDPR with twonozzle printing head. These three constraints form the basis for the objective function for the path optimization problem.
To find out the optimal path of the printing head, optimization algorithms need to be determined to solve such objective functions for the path planning. For path planning, numerous optimization algorithms have been proposed, and they can be divided into two groups: (1) classical methods and (2) evolutionary methods. Some of the most used classical methods are cell decomposition (CD), potential field (PF), roadmap, and subgoal network (SN) [8]. However, the main drawbacks of the classical methods are that they are difficult to find an optimum path for the complex system [9]. In addition, the time for the determination of feasible collision-free path increased, and these algorithms might be trapped in the local optimal solution [7]. Hence, evolutionary algorithms, such as genetic algorithm (GA) [10,11], particle swarm optimization (PSO) [12,13], and differential evolution (DE) [14,15], which randomly find solutions in a given search space can be reasonable alternatives for complex problems. Hence, GA, DE, and PSO algorithms are applied in this research to find the optimal paths of the printing head of the two-nozzle CDPR 3D printer.
For the smooth path generation, the non-uniform rational B-Spline (NURBS) is applied. The NURBS curves have been noticed by research communities for robotic path planning [16][17][18] because it has a high level of flexibility and a great capacity to produce natural smooth curves [18].
Regarding the two-nozzle CDPRs for 3D printing building construction, in addition to the path planning, there is another important issue that needs to be handled. Due to the inevitably flexible characteristics of cables, the two-nozzle CDPRs for 3D printing building construction are highly capable of vibrating in both axial and transversal directions [19]. It is obvious that the vibration gives an adverse effect on the performance of the two-nozzle CDPR. Hence, vibration analysis for the two-nozzle CDPR for 3D printing building construction should be conducted. In this research, the finite element method which was reported in [5,21] is applied to study the vibration of the two-nozzle CDPR 3D printer for the first time. The natural frequency which is a good parameter to evaluate the stiffness and dynamic performance of CDPRs is computed, and its validity is confirmed through the comparison with those of commercial software, SAP2000.

The two-nozzle cable driven parallel robot
This section presents a novel conceptual design of the twonozzle CDPR which could reduce the printing-time dramatically. In addition, the dynamic and kinematic analysis of the printing head and nozzles are also derived.

Conceptual design
A two-nozzle CDPR for constructional 3D printing saving printing time and cost of construction simultaneously is proposed in Fig. 1. A number of nozzles are selected as two considering the fact that most of the contour of the buildings to be printed are closed shape. Two-nozzles are attached to the printing head as shown in Fig. 1. In the figure, an embodiment of a fully constrained two-nozzle cable-driven parallel robot includes a printing head 1, two-nozzles 2 and 3, the sliding bar that extends the nozzle length 12, four vertical support columns 11, drums 8 and 9 mounted on each vertical columns, lower cables 6, and upper cables 7. The pair of cables 6 and 7 are twisted around the drums 9 and 8, respectively, and are driven by actuators 4 and 5. These motors are placed outside the workspace and are placed on the ground. Each drum 9 is mounted on a bracket 10. The bracket can move up and down along the vertical support column by using a motor. Thus, it makes the present embodiment reconfigurable and helps avoid the interference between its system and building parts being printed. Overall, this system has eight cables to suspend the printing head. The poses of the printing head and the nozzles are adjusted by changing the lengths of the cables routed through drums. Figure 2 shows a front view of the two nozzles 2 and 3. These nozzles can move along the axis of the sliding bar 12 and can move independently to print the asymmetric parts at the same time. The sliding bars 12 can rotate around the printing head to flexibly print various printing shapes.
In the case of printing the same-shaped multiple parts or symmetric parts like the walls 15 and 16 as in Fig. 3a, the lengths of the nozzles are fixed. In contrast, for printing an un-symmetric part 17 as shown in Fig. 3b, each nozzle is controlled independently and differently. Hence, proper path of the printing head is determined first, and then path and movement speed of each nozzle is computed subsequently. Feed rate of constructional printing materials such as concrete or foam for each nozzle is controlled in proportion to the movement speed of each nozzle by two separate pumps 14, and they are supplied by two separate supply tubes 13 in Fig 1.

Dynamic equation of the printing head
This section presents the dynamic equation of the printing head. The dynamic behavior of the printing head can be considered the dynamics of the end effector of the general cable robot. As reported in [19,20], a general scheme of the The first kinematic equation is given as follows: Next, differentiation of the previous equation with respect to time and arranging the nth equations into a matrix form yields: where, In these above equations, ̇ and are the linear velocity vector of point C and the angular velocity vector of the end effector, respectively; J is the [n × 6] Jacobian matrix of the cable robot. Thus, the general dynamic equation of cable robot can be presented as follows [19,20]: where m is the mass of the end effector; I P is the inertia tensor of the end effector about point C in the frame O 0 ; I 3x3 is the 3x3 identity matrix; g is the gravity acceleration vector; f e and e are the external force and moment vector applied to the printing head, respectively; and T is the vector of cable forces. It should be noted that cables cannot be compressed, and the vector of cable force must be satisfied with the following force closure equation: where w is the resultant wrench applied at the end effector defined as follows: Fig. 4 A general scheme of the printing head [19] The determination of the cable tension is a key to the operation of the CDPR. The end-effector of the proposed CDPR operates with eight cables, but the motion of the end effector is governed by only six dynamic equations of motion defined in Eq. (4). This implies an under-determined problem; hence, eight cable tensions need to be solved by some form of optimization algorithms. Hence, a nonlinear programming form optimization problem is posed as follows for the determination of the cable tension distribution: where C and c are the weighting factors for the objective function in Eq. (7) and T max and T min are the minimum and the maximum tension allowed to cables, respectively.
Note that 2-norm is used as the objective function of the quadratic programming because discontinuities could result when 1-norm is used as reported in [1].
The weighting factors are set as follows: C = I for the quadratic part and c = 0 for nullifying the linear part.

Fixed two-nozzles mounted to the main printing head
When the two additional nozzles are fixed at the printing head to print the same-shaped objects, these nozzles and the printing head become one rigid body. Then, two nozzles will result in the external forces (f 1 ) and external moments (τ 1 ) to the printing head by the following two equations: where m 1 and m 2 are the mass of the two nozzles, respectively, and d 1 and d 2 are the relative positions of the two with respect to the printing head, respectively.

Moving-nozzles mounted to the main printing head
Generally, it is assumed that the end effector of a CDPR for the 3D printing construction moves with small velocity and constant speed in most cases. Hence, we can neglect the dynamic effect of the printing head on the motion of the moving nozzles. The moving nozzles exert the external forces (f 1 ) and external moments ( 1 ) to the printing head as in Eq. (10), and they can be controlled by two linear motors attached to the printing head. As the sliding bar on which two nozzles slide is assumed to have very low friction, the dynamic equation of two moving nozzles are given as follows: where d 1 and d 2 are the acceleration of the two nozzles and F 1 and F 2 are the applied forces by the nozzle-driving actuators to the two nozzles.

Free vibration analysis
This section presents a finite element formulation to analyze the vibration of the proposed system. Then, some numerical analyses are performed to investigate its vibration behavior.  (N) [60,80,100,120] The maximum tension T max (N) 1000

Cable formulation
For vibration analysis, cables attached between drums and the printing head are divided into small cable elements. Each cable element contains two nodes, i and j, with the length of l ij and the applied tension of T ij . Then, the general equation of motion of cable element is as given follows [5,21]: where M is the mass matrix; and u are the acceleration and displacement vector, respectively; K L is the conventional stiffness matrix; and K G is the geometric stiffness matrix. They are defined below: = u xi u yi u zi u xj u yj u zj T where ij is the weight per unit length; k s is the elastic stiffness with k s = EA ; and G is the transformation matrix which transforms the element stiffness matrix of each cable in the local frame (oxy) to the global frame (OXY) and is given as follows:

Numerical examples
This section presents a free vibration analysis of the proposed CDPR with two nozzles. Some assumptions that should be noted for modeling of the proposed CDPR with two-nozzles are given as follows: • In 3D printing construction, the printing head moves slowly and the dynamic effect of the printing head on the motion of moving nozzles could be neglected. • The nozzles are assumed to be fixed to the printing head hence the nozzles, and the printing head are considered to be one rigid body. • All cables are connected at the centroid of the printing head. • The total mass of the printing head and the additional nozzles is the point mass at the connecting point • Applied cable tensions are calculated by tension distribution algorithms explained in Sect. 2.2.
For the free vibration analysis, the general equation of motion for the cable robot's structure can be given as follows: where M and K are the global mass matrix and the global stiffness matrix, respectively.
The equation for the eigenvalue problem is as follows: where ω and are the natural frequency and the associated displacement vector. To get the nontrivial solutions of the eigenvalue problem, the following determinant must be equal to zero: The vibration analysis presented here is implemented by using a FEM model built by the CDRP model presented in Fig. 1 in the Matlab R2018b program on a desktop computer with Intel Core(TM) i7-7700 CPU @ 3.60 GHz (8 CPUs), 8.00 GB RAM of memory, and Windows 10 Professional with 64-bit operating system. The dimension and the parameter values used for the analysis are summarized in Tables 1  and 2. Fig. 5 shows the SAP2000 model of the proposed systems and parameters for one of the cables. Table 3 shows the first  Table 3 and Fig. 6, it is revealed that there is a tendency that the stiffness of the proposed CDPR is increasing with the minimum cable tension T min . For example, when T min increases from 60 to 120 N, the first natural frequency pulls up its values from 13.7 Hz to 16.4 Hz. Since the stiffness and the vibration resistance of the structure is related to the natural frequency, the stiffness and the vibration resistance of the proposed CDPR 3D printer can be raised by increasing the value of T min. In Table 3, FEM analysis results were compared with the result of commercial software SAP2000 for the verification. It is observed that the values obtained by FEM analysis are almost identical to the results by SAP2000 and the minimum and maximum values of mean errors of 10 modes are only 0.3988 and 0.2183, respectively, thereby, the validity of the FEM analysis is now confirmed. Here, it should be pointed out that modeling a CDPR and making a variation  for model change for other types of CDPRs using the commercial software requires a lot of time. While the FEM model specifically designed for a CDPR can be easily modified for other comprehensive CDPR models. Also, FEM model provides general formulations of cables for modeling of CDPRs which SAP2000 cannot [5]. Moreover, FEM is more economic in terms of computation time as shown in [5]. Figure 7a shows the circular trajectory for the proposed CDPR 3D printer with a radius of 0.75 m on the plane Az = 1.25 m. Figure 7b illustrates the first natural frequency of the proposed CDPR 3D printer along the circular trajectory with two different values of T min . In the case of T min = 100 N, at every pose of the trajectory, the proposed CDPR 3D printer always has the higher first natural frequency than the one in the case of T min = 80 N. It implies that the stiffness of the proposed CDPR 3D printer can be improved by increasing the values of T min . Figure 8 illustrates the first natural frequency of the proposed CDPR 3D printer for different positions in 3D space. From the figure, it is shown that as the printing head moves up, its natural frequency gets increased, and at the center of the workspace, the natural frequency reaches its maximum. Based on this observation, it is inferred that the structural stiffness of the CDPR at a low position in the workspace is low and vulnerable to vibration. In contrast, if the end effector moves towards the top center of the workspace, the structural stiffness increases; thus, vibration resistance is achieved.

Path optimization
For a two-nozzle CDPR, the path optimization of the printing head is conducted based on three constraints: (1) minimize the path lengths difference between two nozzles in the printing process to balance the required motion for two nozzles; (2) achieve the shortest path length to save the printing time; and (3) the path shall be smooth based on the nonuniform rational basis spline (NURBS). Figure 9 shows the 2D contour of an arbitrary shape of the asymmetric part. As mentioned, to print this part, an optimal path of the printing head of Fig. 1 must be determined ahead with three following criteria being satisfied: (1) the path for the printing head should be positioned in the middle of the contour to save and balance the time required to Fig. 11 The flow chart of path optimization move both nozzles. This means that the difference between the path lengths of two nozzles shall be minimized and both nozzles will travel as similar length as possible. (2) The path length of the printing head shall be minimum to save the printing time; (3) the path shall be as smooth as possible by using a NURBS curve.

The objective function for path optimization
In Fig. 9, the contour is constructed with a set of discretized points-the left contour 20 consisting of discretized points 24 of L 1 , L 2 ,…, L n+1 and the right contour 21 consisting of discretized points 25 of R 1 , R 2 , ..., R n+1 . The optimal path 22 of the printing head should go through the start point 18, target point 19, and inside points 23 of P 1 , P 2 ,…, P n+1 , and it should satisfy the first two conditions of the optimal path.
With the given conditions, the optimal path of the printing head can be determined by solving the optimization problem minimizing the following objective function: where n is the numbers of discretized points (which are defined by users); P i is the ith inside points; L i is the ith node on the left boundary; R i is the ith node on the right boundary; and α is the penalty parameter and set to be 10 2 to ensure that two parts in Eq. (22) are equally weighted.

NURBS interpolation for paths
For the optimal path to satisfy the third condition-the smoothness-a non-uniform rational B-spline (NURBS) is used because they have been noticed by research communities for robotic path planning [16][17][18] for its high level of flexibility and a great capacity to produce natural smooth curves [18]. In addition to such features, it is often used for the optimal path generation as it has ability to represent the geometrical shapes in a compact form as reported in [18]. A NURBS curve is defined by its order, a set of weighted control points, and a knot vector. Its schematic view is shown  Fig. 12 The shape of the printing object in example 1 in Fig. 10. It is defined as Eq. (23) and is used for the optimal path: where P i are the n control points (n > 2) and R i,p (u) are the rational B-spline basis functions and are defined as follows: N i,p (u) denotes the i th B-spline basis function of parameter u and degree p, i , stands the weights of control points.
The optimization problem to find the inside points which compose an optimal printing head path is solved by evolutionary algorithms such as GA, DE, and PSO which are the most-used evolutionary algorithms in engineering optimization and path optimization. The sequence of path optimization is shown in the flow chart of Fig. 11.
A numerical simulation is performed in Sect. 4.3 to show the feasibility of GA, DE, and PSO in finding the optimal path of the two-nozzle CDPR.

Computer simulation of path planning
In the following simulation, we assume that the initial and target points are already determined. The main objective is to find the optimal paths that have a minimum length. To find out the points of the optimal paths, genetic algorithm (GA), particle swarm optimization (PSO), and differential evolution (DE) are used. Parameters for these algorithms are provided in Table 4, where GE is the number of generations, NP is the number of populations, CR is the crossover probability, MU is the fraction of mutation, SE is the percentages of selection, θ is the inertia factor and is equal to 0.6, and ε 1 and ε 2 are learning parameters or acceleration constants, respectively, and are set to be 2 [22].
Also, NURBS curves in which all weights of control points set to one are used for the optimal paths of the printing head of a two-nozzle CDPR 3D printer [18].
Let us consider an arbitrary shaped building shown in Fig. 12. The contour of the building to be printed is divided into two sides, the left-and right-half contours. The coordinates of 7-discretized points along these contours are given in Table 5.
The path optimization by GA, PSO, and DE is implemented using Matlab R2018b program on a desktop computer with Intel Core(TM) i7-7700 CPU @ 3.60 GHz (8 CPUs), 8.00 GB RAM of memory, and Windows 10 Professional with 64-bit operating system. The optimal points computed by GA, PSO, and DE are shown in Table 6, and the optimal paths obtained by NURBS based on those points are shown in Fig. 13. From the figure, it is confirmed that three algorithms generated almost similar optimal paths and the optimal paths are not only located in the middle of the printing parts but also forms the shortest paths. Note that the smoothness requirement of paths is satisfied by using NURBS curves. Table 7 shows the travel distances of both nozzles obtained from those three algorithms, and it confirms that the optimal path is located in the middle. Since they travel almost same distances, the speed of two nozzles would be quite similar. The similar speed patterns of both nozzles help suppress the possible vibration that could be incurred when the nozzle speed on both sides unbalances and large acceleration/deceleration is required to compensate the speed unbalance.
For the comparison of the effectiveness of the algorithms, convergence rate of the fitness functions of those algorithms are illustrated Fig. 14.
From the figure, PSO and DE have faster convergence speeds; thus, they require a shorter computation time than GA; hence, they fit better for the real time application. Note that the cost in the figure represents the length of the optimal path computed from Eq. (22).  From the above results, evolutionary algorithms such as GA, DE, and PSO are validated to be effective for the path planning of the two-nozzle CDPR which has never been handled before.

Conclusion
A novel two-nozzle CDPR 3D printer is proposed for the reduction of printing time and the construction cost. The kinematics and dynamics of the printing head and nozzles were analyzed. Due to slow and steady motion of the printing head in 3D printing construction, dynamic effect of the nozzles on the printing head was neglected. For the path generation of the two-nozzle printing head, three constraints-minimizing the path lengths difference between the two nozzles, keeping the path length shortest to save the printing time, and maintaining the smoothness of the path using the NURBS-were imposed. Three evolutionary algorithms (GA, DE, PSO) were used as optimizers for the path optimization problem having above mentioned constraints, and they successfully generated shortest paths which is suitable to control both nozzles mounted on the printing head.
To investigate the vibration of the two-nozzle CDPR, its natural frequencies were computed by using by using the FEM model implemented in the Matlab R2018b program, and its results were verified with a comparison with the results of commercial software (SAP2000). Through natural frequency analysis, it was found that the higher stiffness or the higher natural frequencies could be achieved by increasing the minimum cable tension constraint used in the tension distribution optimization algorithm. It is also revealed that the natural frequency goes up as the printing head moves up towards the top center of the workspace. These features could serve as a design guide for two-nozzle CDPR for 3D printing building construction.