Kinematic Calibration of a Cable-Driven Parallel Robot for 3D Printing

Three-dimensional (3D) printing technology has been greatly developed in the last decade and gradually applied in the construction, medical, and manufacturing industries. However, limited workspace and accuracy restrict the development of 3D printing technology. Due to the extension range and flexibility of cables, cable-driven parallel robots can be applied in challenging tasks that require motion with large reachable workspace and better flexibility. In this paper, a cable-driven parallel robot for 3D Printing is developed to obtain larger workspace rather than traditional 3D printing devices. A kinematic calibration method is proposed based on cable length residuals. On the basis of the kinematic model of the cable-driven parallel robot for 3D Printing, the mapping model is established among geometric structure errors, zero errors of the cable length, and end-effector position errors. In order to improve the efficiency of calibration measurement, an optimal scheme for measurement positions is proposed. The accuracy and efficiency of the kinematics calibration method are verified through numerical simulation. The calibration experiment based on the motion capture system indicates that the position error of end-effector is decreased to 0.6157 mm after calibration. In addition, the proposed calibration method is effective and verified for measurement positions outside optimal positions set through experiments.


Introduction
Cable-driven parallel robots (CDPRs) are known as a type of parallel robots. In CDPRs, the end-effector is suspended by several flexible cables, taking the place of rigid links in traditional rigid-link parallel robots [1][2][3]. Compared with traditional rigid-link parallel robots, CDPRs have much smaller inertia and higher payload to weight ratio, which provides high speed and acceleration of the end-effector [4]. In addition, due to the extension range and flexibility of cables, CDPRs can be applied in challenging tasks that require motion with large reachable workspace and better flexibility as well [5]. Three-dimensional (3D) printing technology has been greatly developed in the last decade and gradually applied in the construction, medical, and manufacturing industries. Limited workspace and accuracy restrict the development of 3D printing technology [6,7]. In this paper, a cable-driven parallel robot for 3D Printing is developed, in order to obtain larger workspace than traditional 3D printing devices.
The accuracy of CDPRs is mainly affected by the kinematic parameter errors. The kinematic parameter errors consist of the geometric parameter errors and the zero errors of the robot joints. It is necessary to compensate the geometric parameter errors and zero errors in order to realize the accurate position control of CDPRs [8]. The kinematic calibration is an effective method to improve the end-effector position accuracy by improving the precision of kinematic model. In general, the kinematic calibration can be divided into error modeling, error measurement, error identification, and error compensation [9][10][11]. The error measurement and error identification are two important steps in the calibration process. The error measurement is the process of measuring end-effector positions by means of the sensors, and comparing the measured values with the theoretical values. The error identification is a process of analyzing error measurement results based on a predetermined parameter identification equation, aiming to obtain kinematic parameter errors. Thus, the accuracy of the error identification is determined by error measurement, including the accuracy of measuring sensors and the selection of measuring positions. There are some areas in the workspace where the kinematic parameter errors have little influence on the end-effector position errors, while the influence of the measurement sensors noises on the accuracy of the error measurement is greatly increased. Therefore, the accuracy of error measurement can be improved obviously by selecting the end-effector positions with large position errors in the workspace as the measurement positions. According to the difference of measurement methods, the calibration methods can be divided into self-calibration method [12][13][14][15] and external calibration method [16][17][18][19].
In general, additional sensors in joints are necessary to obtain joint variables with self-calibration methods. In CDPRs, the end-effector is suspended by several flexible cables instead of rigid links in traditional rigid-link parallel robots, the installation of the sensors in robot joints will decrease the transmission accuracy of cables. Therefore, the external calibration method is more widely used for CDPRs. For the kinematic calibration of robots, many scholars have done a lot of research on improving the precision and efficiency of kinematic calibration, and proposed many efficient calibration methods. Joshi and Surianarayan [20] developed a six degrees of freedom cable-driven parallel device and proposed a kinematic calibration method based on the rotary angle errors of end platform. During the kinematic calibration, the assembly error of the inclinometers is considered and error angle on the perpendicularity of the two inclinometers is taken as the parameter to be calibrated. Ren et al. [21] has proposed a new orientation constraint that keeps the two orientation of the end-effector invariant, and deduces a calibration algorithm using different angle combinations. Lau et al. [22] has proposed a kinematic calibration method based on relative cable length measurement that is used to calibrate the initial cable length and end-effector poses of the CDPR. Duan et al. [23] employed the calibration of a cable-driven parallel manipulator by the six degrees of freedom laser tracker, meanwhile they verified the correctness and to what extent the calibration improves the system through the trajectory planning and motion control. Mustafa [24] developed a novel cable-driven manipulator with seven degrees of freedom, according to the characteristics of its redundant drives, a self-calibration scheme based on the cable joint variables is proposed.
The above calibration methods can achieve an effective calibration, however most of them need to measure a large number of redundant measurement poses in order to solve nonlinear problem in the calibration process. This is a very time-consuming operation process. Therefore, it is necessary to carry out the optimal selection of the measurement poses, which can reach the acceptable accuracy with fewer measurement configurations. The observability of the measurement configurations was used for this purpose. Borm and Menq [25] put forward the concept of observability index, which is used to evaluate the observability of the identification Jacobian matrix for kinematic errors, thus providing a theoretical basis for measurement positions optimization. Scholars have different opinions on the specific value of observability index. Borm and Menq [26] proposed the observability index based on the product of all non-zero singular values of identification Jacobian matrix. Driels and Pathre [27] proposed the condition number which equals the ratio of the largest singular value to the smallest non-zero singular value. Regarding the inverse of the condition number as the observability index, Nahvi et al. [28] regarded the minimum value of the singular value of the identification Jacobian matrix as the observability index. Jia et al. [29] proposed an optimization model with two observability indexes, which can take the partial features and the global ones of the singular values into account at the same time. With the observability index determined, the additional challenge is to find an effective algorithm to construct a configuration set whose observability index is maximal. Li et al. [30] and Zhou et al. [31] adopted DETMAX algorithm to construct the configuration set. Wang et al. [32]  proposed the modified annealing algorithm for the polishing robot. These algorithms improved convergence performance but often led to a low convergence rate. Many recent methods used a large pool of randomly selected configuration candidates, and focused on finding a certain number of optimal configurations within the pool [33,34]. These methods usually make improvements on DETMAX algorithm. Because the size of the candidate pool is usually large, a huge number of computations is required to check all candidates in the pool. Therefore, many scholars put forward different methods to choose the smaller candidate pool. Li et al. [30] identified the candidate pool of a six DOF serial robot and compared the characteristics of the optimal selection schemes based on different observability indices. Nategh et al. [35] demonstrated the distribution of the candidate pool of a rigid parallel mechanism, and calibrated the six degrees of freedom rigid parallel mechanism using a minimum number of measurement configurations. Wang et al. [36] proposed an efficient configuration search method that is based on the closed-form mapping from configuration perturbations to singular-value variations. The above researches on the optimal selection of measurement poses provides a systematic implementation method for the optimal poses selection for parallel or serial rigid robots. Due to the unilateral actuating property of cables, some widely used calibration methods cannot be applied in CDPRs directly, which must be modified to meet the special property of cables.
In this paper, a cable-driven parallel robot for 3D Printing is developed. A kinematic calibration method is proposed based on cable length residuals. On the basis of the kinematic model of the cable-driven parallel robot for 3D Printing, the mapping model is established among geometric structure errors, zero errors of cable, and end-effector position errors. In order to improve the efficiency of calibration measurement, an optimal scheme for measurement positions is proposed. The accuracy and efficiency of the kinematics calibration method are verified through simulation. Furthermore, the calibration experiments based on the motion capture system are carried out to demonstrate the high performance and effectiveness of the proposed calibration method.

Description of Experimental Prototype
In this paper, a novel CDPR for 3D printing is developed, which is designed for large-scale 3D printing. The structure of the CDPR for 3D Printing is shown in Figure 1. The extension range and flexibility of cables can significantly improve the workspace of 3D printing device at a lower manufacturing cost. The CDPR for 3D Printing adopts the Fused Deposition Molding (FDM) technology. The motion of the end-effector is mainly divided into scanning motion in each horizontal layer and lifting movement in vertical direction. Therefore, in order to realize 3D printing, the end-effector should have at least three degrees of freedom (3-DOF). At the same time, for the sake of ensuring the scanning accuracy of the device on each print layer, the rotational motion of the end-effector must be avoided and then the planeness of each printing layer can be ensured. These algorithms improved convergence performance but often led to a low convergence rate. Many recent methods used a large pool of randomly selected configuration candidates, and focused on finding a certain number of optimal configurations within the pool [33,34]. These methods usually make improvements on DETMAX algorithm. Because the size of the candidate pool is usually large, a huge number of computations is required to check all candidates in the pool. Therefore, many scholars put forward different methods to choose the smaller candidate pool. Li et al. [30] identified the candidate pool of a six DOF serial robot and compared the characteristics of the optimal selection schemes based on different observability indices. Nategh et al. [35] demonstrated the distribution of the candidate pool of a rigid parallel mechanism, and calibrated the six degrees of freedom rigid parallel mechanism using a minimum number of measurement configurations. Wang et al. [36] proposed an efficient configuration search method that is based on the closed-form mapping from configuration perturbations to singular-value variations. The above researches on the optimal selection of measurement poses provides a systematic implementation method for the optimal poses selection for parallel or serial rigid robots. Due to the unilateral actuating property of cables, some widely used calibration methods cannot be applied in CDPRs directly, which must be modified to meet the special property of cables.
In this paper, a cable-driven parallel robot for 3D Printing is developed. A kinematic calibration method is proposed based on cable length residuals. On the basis of the kinematic model of the cabledriven parallel robot for 3D Printing, the mapping model is established among geometric structure errors, zero errors of cable, and end-effector position errors. In order to improve the efficiency of calibration measurement, an optimal scheme for measurement positions is proposed. The accuracy and efficiency of the kinematics calibration method are verified through simulation. Furthermore, the calibration experiments based on the motion capture system are carried out to demonstrate the high performance and effectiveness of the proposed calibration method.

Description of Experimental Prototype
In this paper, a novel CDPR for 3D printing is developed, which is designed for large-scale 3D printing. The structure of the CDPR for 3D Printing is shown in Figure 1. The extension range and flexibility of cables can significantly improve the workspace of 3D printing device at a lower manufacturing cost. The CDPR for 3D Printing adopts the Fused Deposition Molding (FDM) technology. The motion of the end-effector is mainly divided into scanning motion in each horizontal layer and lifting movement in vertical direction. Therefore, in order to realize 3D printing, the endeffector should have at least three degrees of freedom (3-DOF). At the same time, for the sake of ensuring the scanning accuracy of the device on each print layer, the rotational motion of the endeffector must be avoided and then the planeness of each printing layer can be ensured.   Considering the above two requirements, the end-effector is driven by three cable groups. Each cable group consists of two parallel cables. A sketch of the CDPR for 3D Printing is shown in Figure 2. According to parallelogram principle, the end-effector cannot rotate around the coordinate axis under the constraint of three cable groups. Therefore, the end-effector only has three translational degrees of freedom, and the stable translation of the end-effector along each coordinate axis are realized under the constraint of the six cables and the follow-up spring, which is used to keep the tension on cables. The upper end of the follow-up spring is connected to the slide rail, which can guarantee that the follow-up spring is always in the vertical direction. Undoubtedly, the elasticity of the cables will affect the printing precision. It is necessary to choose the cables which has the high tensile strength. The modulus of the Kevlar49 cable is 861.9 cN/dtex and the modulus of the stainless steel cable is about 254.4 cN/dtex. Therefore, choosing Kevlar49 cables as the driven cables that will effectively reduce the effect of cable elasticity on printing accuracy.
As shown in Figure 2, one end (B 1 ∼ B 6 ) of six cables (L 1 ∼ L 6 ) is connected to the end-effector, and the cable winds around the fixed pulley set, through the cable outlets (A 1 ∼ A 6 ), the other end is connected with the slider of the linear module which is made up of slide rail and ball screw. Three linear modules are used to pull six cables. The coordinate system of the CDPR for 3D printing is shown in Figure 2. OXYZ is the global coordinate frame which fixed to the center of the base, OX 1 Y 1 Z 1 is the local coordinate frame fixed to the end-effector. Considering the above two requirements, the end-effector is driven by three cable groups. Each cable group consists of two parallel cables. A sketch of the CDPR for 3D Printing is shown in Figure 2. According to parallelogram principle, the end-effector cannot rotate around the coordinate axis under the constraint of three cable groups. Therefore, the end-effector only has three translational degrees of freedom, and the stable translation of the end-effector along each coordinate axis are realized under the constraint of the six cables and the follow-up spring, which is used to keep the tension on cables. The upper end of the follow-up spring is connected to the slide rail, which can guarantee that the follow-up spring is always in the vertical direction. Undoubtedly, the elasticity of the cables will affect the printing precision. It is necessary to choose the cables which has the high tensile strength. The modulus of the Kevlar49 cable is 861.9 cN/dtex and the modulus of the stainless steel cable is about 254.4 cN/dtex. Therefore, choosing Kevlar49 cables as the driven cables that will effectively reduce the effect of cable elasticity on printing accuracy.
As shown in Figure  connected with the slider of the linear module which is made up of slide rail and ball screw. Three linear modules are used to pull six cables. The coordinate system of the CDPR for 3D printing is shown in Figure 2. O X Y Z is the global coordinate frame which fixed to the center of the base,

Equivalent of Analytical Structure
Before analyzing the kinematics, equivalent model of the CDPR for 3D printing can be established for simplicity. The schematic of the equivalent model is shown in Figure 3. For the presented 3-DOF CDPR, the center B of the end-effector is regarded as the equivalent end point. A . Therefore, it can be obtained from the principle of parallelogram that:

Equivalent of Analytical Structure
Before analyzing the kinematics, equivalent model of the CDPR for 3D printing can be established for simplicity. The schematic of the equivalent model is shown in Figure 3. For the presented 3-DOF CDPR, the center B of the end-effector is regarded as the equivalent end point. The six cables L i  According to the above parallel conditions, we can simplify the kinematic inverse operation of six cables to the kinematic inverse operation of three cables without losing the generality, and then the computational complexity of inverse kinematic analysis is greatly reduced.
According to the above parallel conditions, we can simplify the kinematic inverse operation of six cables to the kinematic inverse operation of three cables without losing the generality, and then the computational complexity of inverse kinematic analysis is greatly reduced.

Kinematic Error Modeling
There are several factors that reduce the printing accuracy of the CDPR for 3D Printing, including kinematic error, transmission error, nonlinear error of the control system, deformation error, the measurement error of the sensors, and thermal error. The kinematic errors caused by the manufacturing and assembly are the main error sources. The error model is established in order to describe the relationship between the position errors of the end-effector, the geometric parameter errors of the mechanism, and the zero errors of the cable length.
Kinematic analysis is carried out based on the equivalent model shown in x y z ], therefore, the length of cables can be obtained as follows: where ka L are the cable lengths from virtual cable outlets to equivalent end point, B , ko L are the cable lengths from virtual cable outlets to the equivalent end point B at the initial position, kr L are cable lengths variable, which are defined as the differences between ko L and ka L .
The implicit function form of Equation (1) can be written as: One can obtain from Equation (2)

Kinematic Error Modeling
There are several factors that reduce the printing accuracy of the CDPR for 3D Printing, including kinematic error, transmission error, nonlinear error of the control system, deformation error, the measurement error of the sensors, and thermal error. The kinematic errors caused by the manufacturing and assembly are the main error sources. The error model is established in order to describe the relationship between the position errors of the end-effector, the geometric parameter errors of the mechanism, and the zero errors of the cable length.
Kinematic analysis is carried out based on the equivalent model shown in Figure 4. B(x, y, z) can be represented as the reference point of the local coordinate frame. The Euler angle of the moving coordinate frame is [ 0 0 0 ], the global coordinates of the virtual cable outlets A ka are represented as [x ka , y ka , z ka ], therefore, the length of cables can be obtained as follows: where L ka are the cable lengths from virtual cable outlets to equivalent end point, B, L ko are the cable lengths from virtual cable outlets to the equivalent end point B at the initial position, L kr are cable lengths variable, which are defined as the differences between L ko and L ka . The implicit function form of Equation (1) can be written as: One can obtain from Equation (2) that the kinematic errors of the device consist of the position error of virtual cable outlets that is represented as ( dx ka dy ka dz ka ), and the error of initial cable length (i.e., the zero errors of the cable length) that can be represented as dL k , k = (1, 2, 3) . Thus, a total of 12 geometric error sources can be obtained.
Equation (2) can be differentiated to obtain the following equation: Equation (3) can be expanded and written into matrix form. The mapping relationship of geometric parameter errors, zero errors of the cable length and the end-effector position errors can be obtained as follows: where: where δ is the position error of the end-effector, and ∆q is kinematic error of the CDPR for 3D printing. The matrix C are reversible in nonsingular positions, thus, Equation (4) can be rewritten as follows: where where J 0 is the identification Jacobian matrix.

Optimal Selection
During the process of measuring position, the noise of the measurement sensors is inevitable. The observability is considered as a criterion to find out optimal measurement positions where the position errors of the end-effector are most distinguishable. The noise of the measurement sensors slightly affect the error measurement results in optimal measurement positions.
Optimal selection for measurement positions in this section proceeds in two stages. Through

Optimal Selection
During the process of measuring position, the noise of the measurement sensors is inevitable. The observability is considered as a criterion to find out optimal measurement positions where the position errors of the end-effector are most distinguishable. The noise of the measurement sensors slightly affect the error measurement results in optimal measurement positions.
Optimal selection for measurement positions in this section proceeds in two stages. Through analyzing the position errors distribution in the workspace of the CDPR for 3D printing, the preliminary area of the measurement positions are determined. The coordinates of each optimal position are accurately determined through the calculating the observable index of each point in the preliminary selection region.

Analysis of Workspace and Preliminary Selection
The workspace of the CDPR for 3D printing is the position set that the end-effector can reach within the structural frame. An important characteristic of CDPRs is well known, as cables can only be driven by positive tension in order to keep the straight line shape, rather than negative compression. Therefore, the constraint should be satisfied that cables are all in tension (t k > 0) in the process of calculating the workspace.
Due to the low inertia of the CDPRs rather than rigid link robots, the end-effector can obtain higher acceleration. Therefore, the dynamic constraints caused by acceleration should be taken into account in the analysis of the workspace. The specific mathematical description of acceleration constraints is as follows: where a max is the scalar value of the maximum acceleration of the end-effector, A = ( u 1 u 2 u 3 ) is the structure matrix of the robot, u k k = (1, 2, 3) is unit directional vector of each cable, and T = (t 1 , t 2 , t 3 ) T is the tension value of each cable, W is the external force on the end-effector, the value of which is about 15 N, m is the weight of end-effector, a is the acceleration scalar for end-effector, and c is the unit directional vector of end-effector acceleration. As described in Section 3, the motion of the end-effector is equivalent to the movement under the constraint of three cables and the follow-up spring. The preload of the follow-up spring is far greater than the self-weight of the end-effector. The follow-up spring can be regarded as a special cable whose tension changes with the z value of the end-effector. The direction of the spring is along the z axis of the global coordinate. Therefore, the force balance equation of the CDPR for 3D printing can be expressed as follows: where b = ( 0 0 1 ) T is the unit directional vector of the follow-up spring tension, k is the elastic coefficient of the follow-up spring, z k is the z value of the end-effector in the initial state, Z is the z value of the end-effector, and g is the acceleration of gravity. Equation (7) can be rewritten as follows: According to Equation (8), the mathematical description of acceleration constraints can be rewritten as follows: ∀a < a max ∃T > 0 : The solution of equation can be expressed as: where A + is generalized inverse matrix of A , X + ( ma max ma max ma max ) T is a vector in null space of A , according to matrix theory, r · [X + ( ma max ma max ma max ) T ] is still a vector in null space of A , while r is an arbitrary constant, therefore, Equation (10) can be transformed into: The inertial force generated by the acceleration in specific direction will produce a large resolution force on the cables. Hence, it is necessary to apply the solution method of the force-closure workspace [37] to the solution of the workspace that is constrained by the acceleration. It can be obtained from Equation (11) that if ∃X > 0, T > 0 can be satisfied while r trend to be infinite. According to matrix theory, if ∃V ∈ ker(A ), ∃X > 0 can be satisfied while V > 0. Therefore, the solution conditions of the workspace can be expressed as: According to the need of actual printing, in addition to guaranteeing the tension of cables are all under acceleration constraints, the following constraints should be satisfied during the simulation: (1) The rigid frame constraint on the range of the workspace. Specific constraints can be expressed as follows: where l 0 is the side length of a regular triangle with three virtual cable outlets as apexes, and l z is the initial height of the end-effector.
(2) The cable length constraints on the range of the workspace. The device is driven by the fixed length cables. When the end-effector is in the initial position, the cable length from the cable outlets to the end-effector is the maximum l max . The length of each cable at the remaining positions is less than or equal to l max . Hence, the specific constraints can be expressed as follows: (3) The constraint of cable inclination relative to the horizontal plane. When the inclination of the cables α k k = (1, 2, 3) is too small, it can be known from the decomposition of the force that the force on each cable will increase dramatically. Too large cable tension will affect the performance of the motors. The specific constraints are expressed as follows: (15) As shown in Figure 4, the whole frame can be sketched as a tri-prism structure. The bottom of the tri-prism is a regular triangle with the side length 600 mm, and the height of the tri-prism is 650 mm.
According to the constraints and structural parameters listed above, the workspace of the end-effector is obtained as Figure 5: As shown in Figure 4, the whole frame can be sketched as a tri-prism structure. The bottom of the tri-prism is a regular triangle with the side length 600 mm, and the height of the tri-prism is 650 mm. According to the constraints and structural parameters listed above, the workspace of the end-effector is obtained as Figure 5: Figure 5. The workspace of the CDPR for 3D Printing.
The shape of the workspace can be sketched as a tri-prism. The bottom of the tri-prism is a regular triangle with the side length 180 mm and the height of the tri-prism is 240 mm. During the experimental stage, in order to ensure the stability of the CDPR for 3D Printing, the CDPR for 3D Printing is driven by the fixed length cables. However, it results in a smaller workspace. As mentioned above, the cable length and the rigid frame are the main constraints of the workspace. Therefore, the workspace can be expanded through adopting the winches to drive the CDPR for 3D printing or enlarging the scale of the rigid frame. Due to the extension range and flexibility of cables, the workspace can be expended at a low engineering cost. The shape of the workspace can be sketched as a tri-prism. The bottom of the tri-prism is a regular triangle with the side length 180 mm and the height of the tri-prism is 240 mm. During the experimental stage, in order to ensure the stability of the CDPR for 3D Printing, the CDPR for 3D Printing is driven by the fixed length cables. However, it results in a smaller workspace. As mentioned above, the cable length and the rigid frame are the main constraints of the workspace. Therefore, the workspace can be expanded through adopting the winches to drive the CDPR for 3D printing or enlarging the scale of the rigid frame. Due to the extension range and flexibility of cables, the workspace can be expended at a low engineering cost.
In order to reduce the search range of optimal measurement positions and obtain the preliminary positions with greater observability, it is necessary to analyze the distribution of position errors in the workspace. Firstly, random geometric errors (from −1 mm to 1 mm) is added to all desired kinematic parameters. The construction volume errors of the end-effector in the workspace are obtained via forward kinematics. The construction volume error of parallel 3D printing device at different height in workspace can be obtained, the results are depicted in Figure 6. The construction volume errors in the workspace of the CDPR for 3D printing can be expressed as follows: where ∆x is the position error of end-effector in x direction, ∆y is the position error of end-effector in y direction, and ∆z is the position error of end-effector in z direction. As shown in Figure 6, the z values of the end-effector are 110 mm, 150 mm, 190 mm, and 230 mm, respectively. It can be seen that the maximum value of the construction volume errors is distributed at the boundary of the each horizontal plane and the structural volume error of the end-effector decreases with the printing height. Based on the principle of selecting the positions with the larger structural volume errors, the boundary region of the z = 110 mm plane in the workspace is regarded as the preliminary selection area of the measurement positions. Sensors 2018, 18, x FOR PEER REVIEW 10 of 23 Figure 6. The workspace of the CDPR for 3D Printing.
In order to reduce the search range of optimal measurement positions and obtain the preliminary positions with greater observability, it is necessary to analyze the distribution of position errors in the workspace. Firstly, random geometric errors (from −1 mm to 1 mm) is added to all desired kinematic parameters. The construction volume errors of the end-effector in the workspace are obtained via forward kinematics. The construction volume error of parallel 3D printing device at different height in workspace can be obtained, the results are depicted in Figure 6. The construction volume errors in the workspace of the CDPR for 3D printing can be expressed as follows: x  is the position error of end-effector in x direction, y  is the position error of endeffector in y direction, and z  is the position error of end-effector in z direction.
As shown in Figure 6, the z values of the end-effector are 110 mm, 150 mm, 190 mm, and 230 mm, respectively. It can be seen that the maximum value of the construction volume errors is distributed at the boundary of the each horizontal plane and the structural volume error of the end-effector decreases with the printing height. Based on the principle of selecting the positions with the larger structural volume errors, the boundary region of the z = 110 mm plane in the workspace is regarded as the preliminary selection area of the measurement positions.

Optimal Positions Selection
After the preliminary selection area of the measurement positions is determined, the optimal positions selection can be divided into two steps. Firstly, the candidates of measurement positions

Optimal Positions Selection
After the preliminary selection area of the measurement positions is determined, the optimal positions selection can be divided into two steps. Firstly, the candidates of measurement positions should be selected from z = 110 mm. Secondly, an effective algorithm should be applied to construct an optimal position set, the observability index of which is maximal.
As shown in Figure 7a, the z = 110 mm plane in the workspace is an equilateral triangle region with sides equal to 170 mm. In order to facilitate the selection of the candidates of measurement positions, the boundary of the equilateral triangle region is discretized. Through moving the sides of the equilateral triangle inward by 1mm, an embedded equilateral triangle can be obtained. The candidates of measurement positions are distributed along the edges of the two equilateral triangles. The distribution of the position candidates is shown in Figure 7b, 60 candidates of measurement positions are selected.
with sides equal to 170 mm. In order to facilitate the selection of the candidates of measurement positions, the boundary of the equilateral triangle region is discretized. Through moving the sides of the equilateral triangle inward by 1mm, an embedded equilateral triangle can be obtained. The candidates of measurement positions are distributed along the edges of the two equilateral triangles. The distribution of the position candidates is shown in Figure 7b, 60 candidates of measurement positions are selected. As mentioned in Section 3.2, the CDPR for 3D printing has 12 independent kinematic parameter errors, and each measurement position has 3 kinematic constraint equations. In order to accomplish the error identification, the number of measurement positions K is at least 4. Therefore, the optimal selection of the measurement positions in this paper aimed at selecting 4 optimal positions with the maximum observability index from the 60 candidate positions. Four kinds of observability index are mentioned in the introduction. In this paper, the reciprocal of conditional numbers 2 O , is chosen as the observability index. Among four kinds of observability index, the calibration algorithms will be led to a high convergence rate while regarding 2 O as observability index [30]. Regarding the nonzero singular value of the identification Jacobian 0 J , as (1) c Q is the candidate position set, there are C = 60 candidate points in the initial state.  As mentioned in Section 3.2, the CDPR for 3D printing has 12 independent kinematic parameter errors, and each measurement position has 3 kinematic constraint equations. In order to accomplish the error identification, the number of measurement positions K is at least 4. Therefore, the optimal selection of the measurement positions in this paper aimed at selecting 4 optimal positions with the maximum observability index from the 60 candidate positions. Four kinds of observability index are mentioned in the introduction. In this paper, the reciprocal of conditional numbers O 2 , is chosen as the observability index. Among four kinds of observability index, the calibration algorithms will be led to a high convergence rate while regarding O 2 as observability index [30]. Regarding the non-zero singular value of the identification Jacobian J 0 , as σ m ≤ . . . ≤ σ 1 , the observability index O 2 , can be expressed as: where σ m is the maximal non-zero singular value of the identification Jacobian, and σ 1 is the minimal non-zero singular value of the identification Jacobian. By calculating the observability index O 2 , of the identification Jacobian J 0 , in the above 60 candidate positions, the K(K = 4) positions with the maximum observability index O 2 , among the 60 candidate positions are selected. The method of optimal positions selection is presented as follows: (1) Q c is the candidate position set, there are C = 60 candidate points in the initial state.
(2) N l is the optimal measurement positions set, it concludes a number of l = 0 positions in the initial state.
(3) Search the O 2 maximum point q l , in Q c , using the extremum seeking function.
(4) Add q l to N l , l = l + 1, meanwhile, remove q l from Q c , C = C − 1 (5) Repeat steps 3, 4, until l = K, stop the cycle. After the above steps, the result is that the maximum observation index of the optimal positions set is O 2max = 13.99 × 10 2 , the coordinates of each optimal position in N l are shown in Table 1.

Simulation
As mentioned above, the coordinates of the cable outlets and the initial cable length are the main kinematic parameters that determine the printing accuracy of the CDPR for 3D printing. In order to improve the printing accuracy, a kinematic calibration method based on cable length residuals is proposed. Through measuring the end-effector position errors and cable length variables synchronously by the motion capture system, the calibration method can calibrate the coordinates of the cable outlets and the initial cable length. This section mainly focuses on the simulation for the calibration to verify the effectiveness of this method.

Error Identification Model
The main aim of kinematic calibration is to reduce the position errors of end-effector by the way of kinematic parameters compensation. Therefore, it is the key problem of kinematic calibration to establish the error identification model. In this paper, the functional formula for the error identification is established, which realizes the error identification among the position errors of end-effector, geometric parameter errors, and zero errors of cable length.
The difference between the measured and actual values of the end-effector positions can be expressed as follows: where e k is the residuals between the measured and the true position coordinates of the end-effector, that is, the measurement noise of the measurement sensor, e m is the measured position coordinates of the end-effector, and e is the true position coordinates of the end-effector. According to Equation (1), the kinematic constraint equation of the CDPR for 3D printing is as follows: where L ko = L ko + ∆l k , L kr = L kr + ψ r , x = x + ∆x + e kx y = y + ∆y + e ky , z = z + ∆z + e kz , x ka = x ka + ∆x ka y ka = y ka + ∆y ka , z ka = z ka + ∆z ka , e = e + δ, q = q + ∆q where f (x) is the forward kinematics equation of the CDPR for 3D printing. Equation (20) can be rewritten as a functional form.
For the CDPR, the function equation based on the cable length residual is more convenient for the implementation of the nonlinear least square method. Hence, the functional equation for error identification is rewritten as follows: where ζ i is a set of cable length residuals caused by the measurement noise of the measurement sensors.
The error identification problems are usually solved through the nonlinear least square method. Hence, the equation of error identification can be written as follows:

Simulation Verification
The specific simulation calibration process is as follows: (1) The optimal measurement positions set N l , is regarded as the measurement positions set. The solving of inverse kinematic equation is performed based on the nominal values of the measurement positions. The cable length L kr,l , corresponding to the each measurement position can be obtained, which is used to simulate the cable length measured by motion capture system.
(2) By substituting the cable length L kr,l , into the forward kinematic model based on the true values of geometric parameters, the actual coordinates e l , of each measurement position will be calculated. The true values of the geometric parameters contain the true coordinate values of the cable outlets and the true values of the initial cable length, which are given ideal parameter during the simulation. The main purpose of this step is to simulate the measurement value of the motion capture system in the actual calibration process.
(3) The position errors of end-effector δ i , are calculated and substituted into the error identification equation Equation (23). The error identification values ∆q i , are calculated and compensated to each kinematic parameter q i = q i + ∆q i , i = i + 1 and the new nominal value, q i , of each parameter is obtained.
(4) Compare the absolute value of ∆q i with the given threshold ε, if ∆q i ≤ ε, the simulation calibration is terminated, otherwise i = i + 1, return to step (1).
Because the motion control system is based on the equivalent model of the CDPR for 3D printing mentioned in Section 3, the process of simulation calibration is based on the equivalent kinematics model. The true coordinates of the three virtual cable outlets are a 1 = (−260, −150.111, 78), a 2 = (260, −150.111, 78), and a 3 = (0, 300.222, 78), respectively, which are given ideal parameter during the simulation. The initial position of end-effector is specified as (0, 0, 330), the true values of the initial cable length can be determined using inverse kinematics as L ko = 392 mm, (k = 1, 2, 3). Before the kinematic calibration of the CDPR for 3D printing, preliminary measurement should be performed to obtain the nominal values of the kinematic parameters. According to the existing measuring tools, the errors of the preliminary measurement are generally within 3 mm. In the simulation, random errors in the range of ±3 mm are added to the true values of geometric parameters to obtain the nominal values of kinematic parameters. Therefore, it is assumed that the nominal values of the initial cable length are L 1o = 390 mm, L 2o = 391.5 mm, and L 3o = 389 mm, respectively. The nominal coordinates of each cable outlet are a 1 = (−258, −149, 79), a 2 = (263, −148, 77), and a 3 = (0, 301, 78.5), respectively.
Based on the data listed above, the threshold is set as ε = 1 × 10 −4 mm, the simulation calibration of a CDPR for 3D printing is performed. After the parameters iteration, the termination condition is satisfied while the parameter identification time is i = 4. The specific simulation results are shown in Table 2.
As shown in Figure 8a, before the calibration, the maximum construction volume error in the optimal measurement positions set is ∆V 0max = 14.4223 mm. After four times error identification, the maximum structural volume error in the optimal measurement positions set is ∆V 4max = 5.4321 × 10 −9 mm.    Figure 8b shows the variation of the ∆e x , ∆e y , and ∆e z with the error identification times. It can be seen from Figure 8b that before calibration, the average absolute values of the position errors in x, y, and z directions are ∆e x0 = 1.1772 mm, ∆e y0 = 1.9062 mm, and ∆e z0 = 14.1970 mm, respectively. After four times error identification, the mean absolute value of the position errors in the direction of each coordinate axis are reduced to ∆e x4 = 8.4842 × 10 −10 mm, ∆e y4 = 1.2539 × 10 −9 mm, and ∆e z4 = 1.5540 × 10 −9 mm. The experimental results show that this calibration method is of fast convergence speed and favorable calibration accuracy.

Calibration Experiment
The experimental platform of the CDPR for 3D printing is shown in Figure 9, the stepping motor drives the linear module to drive the fixed-length cables to realize the movement of the end-effector. Therefore, the displacements of the sliders on the linear module is equal to the change of cable length. The coordinates of the markers attached to the sliders can be measured with the motion capture system.
Sensors 2018, 18, x FOR PEER REVIEW 16 of 23 Figure 9. Experimental device of kinematic calibration.

Measurement
Before calibrating the coordinate values of each cable outlet and the initial cable length accurately, the initial measurement is carried out in order to obtain the reasonable structural parameters required by the control system. The upper surface of the base is the plane (z = 0 mm), the centroid of the regular triangular base is regarded as the origin of the global coordinate frame. The coordinate values of the cable outlets and the initial cable length measured in the global coordinate frame are shown in Table 3. Table 3. The initial measurement of kinematic parameters. As mentioned in Section 5.2, the coordinates of the cable outlets that should be calibrated are the coordinates of the three virtual cable outlets. According to the measured coordinates of each cable outlet, the coordinate values of virtual cable outlets are calculated based on the equivalent principle described in Section 3.1. The length of each cable group is the same as that of the corresponding virtual cables. The specific coordinate values of virtual cable outlets are shown in Table 4.

Measurement
Before calibrating the coordinate values of each cable outlet and the initial cable length accurately, the initial measurement is carried out in order to obtain the reasonable structural parameters required by the control system. The upper surface of the base is the plane (z = 0 mm), the centroid of the regular triangular base is regarded as the origin of the global coordinate frame. The coordinate values of the cable outlets and the initial cable length measured in the global coordinate frame are shown in Table 3. Table 3. The initial measurement of kinematic parameters.

Kinematic Parameters Parameters without Calibration (mm)
A As mentioned in Section 5.2, the coordinates of the cable outlets that should be calibrated are the coordinates of the three virtual cable outlets. According to the measured coordinates of each cable outlet, the coordinate values of virtual cable outlets are calculated based on the equivalent principle described in Section 3.1. The length of each cable group is the same as that of the corresponding virtual cables. The specific coordinate values of virtual cable outlets are shown in Table 4. Before the measurement, the markers should be adhere to the prototype to obtain the positions of each mobile component. The markers are distributed in three places, as shown in Figure 10. The first maker place is on the end-effector where the markers are triangle-distributed. The second maker place is triangle-distributed on the base. The centroid of the equilateral triangle is the coordinate origin of the global coordinate frame. The third maker place is located on the sliders of linear modules, the differences between the z values of the markers measured by the motion capture system and the z values of markers in the initial state are the cable length variables.  During the data measurement, the end-effector is moved to the above optimal positions (z = 110 mm), respectively. Meanwhile, the actual coordinate value of each position and the cable length are measured by the motion capture system synchronously.

Data Processing
In this calibration experiment, the measurement rate of the motion capture system is 60 frames per second. The effective test time for each location is 3 s. The coordinate of the markers is based on the world coordinate system determined by self-calibration of the motion capture system. However, the world coordinate system and the global coordinate system of the experimental device are difficult to be coincident though artificial setting. Therefore, in order to apply the measurement data to the error identification, the coordinate values measured by the motion capture system should be converted first.
The coordinate system A is set as the world coordinate system, the coordinate system B is set as the global system of experimental device, and the coordinate origin of coordinate system B in the coordinate system A is p is the centroid coordinate values of the regular triangle determined by the markers on the base, the y axis direction vector of the coordinate system B is also determined by the three attached markers on the base, which is the unit direction vector from a marker to the origin of coordinate system B. After obtaining the measurement data as described above, the rotation matrix A B R between the coordinate system A and the coordinate system B can be calculated out. The coordinate values of the point p that is measured by the motion capture system in coordinate system A is A p , the position of point p in coordinate system B is as follows: During the data measurement, the end-effector is moved to the above optimal positions (z = 110 mm), respectively. Meanwhile, the actual coordinate value of each position and the cable length are measured by the motion capture system synchronously.

Data Processing
In this calibration experiment, the measurement rate of the motion capture system is 60 frames per second. The effective test time for each location is 3 s. The coordinate of the markers is based on the world coordinate system determined by self-calibration of the motion capture system. However, the world coordinate system and the global coordinate system of the experimental device are difficult to be coincident though artificial setting. Therefore, in order to apply the measurement data to the error identification, the coordinate values measured by the motion capture system should be converted first.
The coordinate system A is set as the world coordinate system, the coordinate system B is set as the global system of experimental device, and the coordinate origin of coordinate system B in the coordinate system A is A p B0 . Through determining A p B0 and y axis direction vector of the coordinate system B in the coordinate system A, the conversion relationship between the two coordinate systems can be determined. A p B0 is the centroid coordinate values of the regular triangle determined by the markers on the base, the y axis direction vector of the coordinate system B is also determined by the three attached markers on the base, which is the unit direction vector from a marker to the origin of coordinate system B. After obtaining the measurement data as described above, the rotation matrix A B R between the coordinate system A and the coordinate system B can be calculated out. The coordinate values of the point p that is measured by the motion capture system in coordinate system A is A p, the position of point p in coordinate system B is as follows: where A B R −1 is the inverse matrix of rotation matrix A B R. After the data measurement, the coordinate values of three markers (F1, F2, F3) on the base in the coordinate system B are shown in Table 5. According to the data, it is calculated that A p B0 is (542.4536, 48.4085, 23.1651) and the specific value of rotation matrix A B R is as follows: In the motion control that is based on the initial measured kinematic parameters, the coordinate values of each optimal position are measured and the measured data are transformed by the mentioned method. The structural volume errors and the position errors of each optimal position are obtained as shown in Figure 11. It can be seen from Figure 11 that the error along the z axis is the main factor leading to the excessive position errors of the CDPR for 3D printing. Before calibration, the specific range of position errors and structural volume error of each optimal position are as shown in Table 6.  In the motion control that is based on the initial measured kinematic parameters, the coordinate values of each optimal position are measured and the measured data are transformed by the mentioned method. The structural volume errors and the position errors of each optimal position are obtained as shown in Figure 11. It can be seen from Figure 11 that the error along the z axis is the main factor leading to the excessive position errors of the CDPR for 3D printing. Before calibration, the specific range of position errors and structural volume error of each optimal position are as shown in Table 6.  The position errors of the end-effector reflects the accuracy of the CDPR for 3D printing. According to the analysis of the data in Figure 11 and Table 6, it can be seen that the position errors along z axis are more than 20 mm. The accuracy requirement for 3D printing can not be satisfied due to the position errors. In order to solve this problem, the error identification should be implemented based on the above measurement data. Through compensating the initial measured kinematic Δe (mm) Figure 11. The position errors of the each optimal position before calibration. The position errors of the end-effector reflects the accuracy of the CDPR for 3D printing. According to the analysis of the data in Figure 11 and Table 6, it can be seen that the position errors along z axis are more than 20 mm. The accuracy requirement for 3D printing can not be satisfied due to the position errors. In order to solve this problem, the error identification should be implemented based on the above measurement data. Through compensating the initial measured kinematic parameter with identification results, the printing accuracy will be improved.
In order to obtain valid error identification results, it is necessary to set a reasonable convergence threshold based on the measurement noise of the motion capture system. Through measuring the position coordinates of one fixed point repeatedly, the measurement noise of the motion capture system can be obtained. As shown in Figure 12, the fluctuation of the measurement noise is −0.12939 mm ≤e k ≤ 0.22407 mm. To reduce the influence of measurement noise on the error identification, a threshold as ξ = 0.80 mm is set. When the absolute value of the end-effector position errors meet the constraint |δ| ≤ ξ, the error identification cycle is terminated and the calibration is completed. After the parameter iteration, the position errors of the end-effector changing with the identification times are shown in Figure 13, after calibration, the errors of the end positions are as shown in Table 7. The identification errors of kinematic parameters are shown in Table 8.  Figure 13, after calibration, the errors of the end positions are as shown in Table 8. The identification errors of kinematic parameters are shown in Table 7.     Figure 13, after calibration, the errors of the end positions are as shown in Table 8. The identification errors of kinematic parameters are shown in Table 7.   Δe (mm) Figure 13. The relations between position errors and identification times. From the above experimental results, one can obtain that the construction volume error of the end-effector position is reduced from 23.4662 mm to 0.4740 mm through the calibration. Among the position errors in each direction, the kinematic calibration is obviously effective to reduce the position errors along the z axis. The average errors in the z direction are reduced from −23.0636 mm to 0.3594 mm. Different from the position errors in the x and y directions, the position errors along the z axis will accumulate on the printed components as the printing layers increases for the fused deposition modeling. The proposed calibration method can bring obvious improvement of the position accuracy in the z direction and therefore, improve the printing accuracy of the CDPR for 3D printing effectively.

Calibration Result Verification
The calibrated kinematic parameters are applied to the control system, and the end-effector is controlled to move in four planes: Z = 175 mm, Z = 195 mm, Z = 215 mm, and Z=235 mm. On each plane, the ideal motion curve of the end-effector is the straight line from (20,20) (20,20). The static position errors measurement is carried out to eliminate the influence of delay of the controller, actuator, spring, etc. In the static case, the error of the end-effector mainly comes from the error of kinematic parameters.
One can obtain from Figure 14 that the measured positions covers four planes with different heights and the position errors of the calibrated CDPR for 3D printing in each direction of the end-effector remains within −0.80 mm ≤ δ ≤ 0.80 mm. The measurement results indicate that the proposed calibration method is effective and verified for measurement positions outside optimal positions set. Sensors 2018, 18, x FOR PEER REVIEW 21 of 23 Figure 14. The error distributions of the end-effector after calibration.

Conclusions
This paper presents a CDPR for 3D printing, which can improve the workspace of 3D printing device at a lower manufacturing cost, due to the extension range and flexibility of cables.
In order to solve the kinematic calibration problem of CDPR for 3D printing, a kinematic calibration method based on cable length residuals is proposed. The functional formula is established between the measured information of the motion capture system and the actual kinematic parameters. According to the structure characteristics of the CDPR for 3D printing, a coordinate system conversion method for data measurement of the motion capture system is proposed, which provides the guarantee for obtaining accurate measurement data. The accuracy and effectiveness of this calibration method are verified though simulation. Furthermore, the kinematic calibration experiment is carried out on the basis of synchronous measurement of end-effector position errors and cable length variables with motion capture system.
In order to further improve the efficiency of calibration and measurement, an optimal selection scheme for measurement positions is proposed. The simulation results of error modeling and analysis indicate that the accuracy of CDPR for 3D printing increases with printing height, and the maximum error on each horizontal plane is distributed close to the boundary of the workspace. The calibration experiments are carried out based on the optimized positions set, and the construction volume error of the end-effector is reduced from 23.4805 mm to 0.6157 mm. In addition, the proposed calibration method is effective and verified for measurement positions outside optimal positions set through experiments.

Position error (mm)
Position error (mm) Position error (mm) Position error (mm) Figure 14. The error distributions of the end-effector after calibration.

Conclusions
This paper presents a CDPR for 3D printing, which can improve the workspace of 3D printing device at a lower manufacturing cost, due to the extension range and flexibility of cables.
In order to solve the kinematic calibration problem of CDPR for 3D printing, a kinematic calibration method based on cable length residuals is proposed. The functional formula is established between the measured information of the motion capture system and the actual kinematic parameters. According to the structure characteristics of the CDPR for 3D printing, a coordinate system conversion method for data measurement of the motion capture system is proposed, which provides the guarantee for obtaining accurate measurement data. The accuracy and effectiveness of this calibration method are verified though simulation. Furthermore, the kinematic calibration experiment is carried out on the basis of synchronous measurement of end-effector position errors and cable length variables with motion capture system.
In order to further improve the efficiency of calibration and measurement, an optimal selection scheme for measurement positions is proposed. The simulation results of error modeling and analysis indicate that the accuracy of CDPR for 3D printing increases with printing height, and the maximum error on each horizontal plane is distributed close to the boundary of the workspace. The calibration experiments are carried out based on the optimized positions set, and the construction volume error of the end-effector is reduced from 23.4805 mm to 0.6157 mm. In addition, the proposed calibration method is effective and verified for measurement positions outside optimal positions set through experiments.

Conflicts of Interest:
The authors declare no conflict of interest.