A Workspace Visualization Method for a Multijoint Industrial Robot Based on the 3D-Printing Layering Concept

: The workspace of a robot provides the necessary constraint information for path planning and reliable control of the robot. In this paper, a workspace visualization method for a multijoint industrial robot is proposed to obtain a detailed workspace by introducing the 3D-printing layering concept. Firstly, all possible joint-angle groups of one pose in the joints’ ranges are calculated in detail according to the POE (product of exponential) theory-based forward-kinematics expressions of the multijoint industrial robot. Secondly, a multisolution selection method based on the key degree of the joint is proposed to select the appropriate joint-angle groups. The key degrees of all joints and their key order are obtained according to the sensitivity expressions of all joint angles, calculated from the Jacobian matrix of the robot. One principle based on the smallest di ﬀ erences of the nominal angle is established to select the possible solutions for one joint from the possible solutions for the joint with the smaller key order. The possible solutions for the joint with the highest key order are the appropriate joint-angle group. Thirdly, a workspace visualization method based on the layering concept of 3D printing is presented to obtain a detailed workspace for a multijoint industrial robot. The boundary formula of each layer is derived by forward kinematics, which is expressed as a circle or a ring. The maximum and minimum values of the radius are obtained according to the travel range of the joint angles. The height limitations of all layers are obtained with forward kinematics. A workspace boundary-extraction method is presented to obtain the array of path points of the boundary at each layer. The proposed postprocessing method is used to generate the joint-angle code of each layer for direct 3D printing. Finally, the e ﬀ ectiveness of the multisolution selection method and the workspace visualization method were veriﬁed by simulation and experiment.


Introduction
Industrial robots play an important role in manufacturing as manufacturing automation and intelligence continuously improve. The industrial robot industry is an important indicator of a country s manufacturing level and technology level. Industrial robots are widely used in different fields, such as automobile manufacturing, electrical and electronics, and aerospace. In addition, industrial robots are core pieces of equipment for the industrial intelligent manufacturing. The workspace of a robot is an important index of the flexibility of the robot. It provides the necessary constraint information for the path planning and control of the robot. It has important significance in the structural design and path planning.
Kinematics represent the relationship between the position and the pose of the end-effector and the joint variables. This lays the basis for the workspace of the multijoint industrial robots. The Denavit-Hartenberg (D-H) parametric method is a standard method of modeling [1][2][3]. Cumulative errors may be increased with the increased degrees of freedom when using this method. POE (product of exponential) theory is another option for kinematics modeling [4][5][6][7]. An et al. gave a generalized solution of a kinematics problem based on POE in order to analyze the kinematics problems of serial robots [8]. Ayiz et al. used POE theory to solve kinematics problems of industrial robots [9]. POE theory allows global description of rigid body motion and greatly simplifies the analysis of mechanisms. For inverse kinematics, the common methods are the algebraic method, geometric method, and numerical method. Yahya et al. proposed a geometric method to find the optimal solution of the inverse kinematics of redundant or hyperredundant robots from infinite solutions [10]. Ananthanarayanan et al. obtained the inverse kinematics solution by using iteration methods with the initial value of the analytical solution [11]. Modern intelligent algorithms are also widely used to solve the inverse kinematics of robots with complex structures. Baheshti et al. used the adaptive fuzzy logic method to determine the inverse kinematics solution of a 3-DOF (degree of freedom) planar robot [12]. Ayyıldız et al. obtained the solution of an inverse kinematics equation by using four different optimization algorithms [13]. These methods can produce accurate inverse kinematics solutions, but they need a lot of training. The inverse kinematics solution of robots is a multisolution problem. The criterion of shortest travel is generally adopted to select the appropriate solution. Wang et al. proposed a criterion of weighted shortest travel, and selected the appropriate solution by finding the minimum value of weighted Euclidean distance [14]. Liu et al. divided the result of inverse solution with the geometric attitude of the robot to facilitate selection [15].
The determination of the workspace boundary is generally an intermediate but critical step in analyzing multijoint industrial robots [16]. Traditionally, the workspace boundaries of multijoint industrial robots are usually expressed as numerical curves on the boundaries or analytical formulas of a certain kind of robot. Abdel-malek et al. used differential geometry and topology methods to obtain general formulas for the workspace recognition and visualization using the manifold layering method [17]. Wang et al. used the surface envelope overlay (SEO) method to identify and visualize the robot workspace [16]. Graphic methods are another way to find workspace boundaries [18]. Gan et al. took advantage of the arm length of a robot and the angle ranges of its joints to obtain section screenshots of the robot's workspace on the xz plane via the geometric drawing method [19]. Liu et al. discussed the geometric description of a robot's workspace and presented the design results of a delta robot for a given workspace [20]. This method is more intuitive, but it is not suitable for robots with more than 3 degrees of freedom. Most of the analysis methods try to extract boundary information from singular points of robots by analyzing the characteristics of the singular points. Based on this theorem, several analytical criteria for determining the workspaces boundaries were derived [21]. Kohli et al. analyzed a workspace and subspace by polynomial criterion and Jacobian matrix, respectively [22,23]. Yang et al. described singular behavior with the row rank defect of the Jacobian matrix and presented the boundary subsurface of the envelope of a workspace using the perturbation method by dividing the singular surface into several subsurfaces [24].
Most numerical methods of singularities or bounds are mapped by grids. All grid-based methods produce an approximate boundary, and the accuracy depends on the size of the grid. The Monte Carlo method is commonly used to generate sample points in a grid in order to identify the grid points near the workspace boundary via the statistical method, which bypasses the necessity of tightness. Cao et al. generated random values of the joint vector data by increasing the probability of the joint vector at the workspace boundary with the upper limit and lower limit of the joint angles [25]. Peidro et al. used the classical Monte Carlo method to generate an inaccurate workspace and used a Gaussian distribution to encrypt and extend the workspace until the boundary of the workspace was reached [26]. Liu et al. added random points locally at the boundary points of each layer in order to generate a clear workspace boundary [27]. However, this method cannot exclude the identification of false holes in the workspace. The layering concept in 3D printing can be introduced to obtain expressions of workspace boundaries, which makes the visualization of the workspace feasible.
In addition, interval-analysis-based methods have been used to analyze workspaces. Interval analysis can be used to evaluate the constraints, and branch-and-prune techniques used to characterize the constraint workspace. This is the basis of prescribed bounded velocity and force transmission factors. Merlet et al. proposed an interval-analysis-based approach to determine almost all the geometries of a simplified Gough platform for which the workspace included an arbitrary set of poses [28]. Chablat et al. used interval analysis to compute a dextrous workspace as well as the largest cube enclosed in this workspace [29]. Kaloorazi et al. proposed a systematic algorithm based on the interval-analysis concept in order to find the maximal singularity-free circle or sphere within the workspace of parallel mechanisms [30]. Viegas et al. used interval analysis to evaluate several performance indexes of parallel kinematic machine design [31]. However, interval-analysis-based methods are mostly used to analyze the workspace of parallel mechanisms.
This paper proposes a new workspace visualization method based on a 3D-printing layering concept in order to obtain and visualize the workspace accurately. At first, all possible solutions of joint angles are obtained in detail based on the forward kinematics and inverse kinematics of a multijoint industrial robot. Next, one multisolution selection method is proposed according to the sensitivity-based key degree of the joints to select the appropriate solutions. Next, a 3D-printing layering concept is introduced to establish one workspace visualization method. The forward kinematics can help to establish the boundary formula of each layer and obtain the limitation of the height of the workspace. The multisolution selection method can help to obtain the joint angles of each layer, which also can be used to print the workspace.
The remains of the paper are organized as follows: In Section 2, the kinematics model of a multijoint industrial robot based on POE theory is established, and all possible joint-angle groups of one pose in the joint range are presented in detail. In Section 3, a multisolution selection method based on the key degree of joints is proposed. Section 4 presents a workspace visualization method based on a 3D-printing layering concept. Section 5 describes the validation of the simulation and experiments that were carried out. Figure 1 shows the structure and structural parameters of a multijoint industrial robot, and Figure 2 shows the kinematic chain. BCS is the base coordinate system, TCS is the tool coordinate system, and WCS is the workpiece coordinate system. The structural parameters V A , V B , V C , V D , and V W represent the positions of the A-joint, B-joint, C-joint, D-joint, and the base relative to the end-effector, respectively. V U represents the position of the workpiece relative to the base. They were represented as

Forward Kinematics Based on POE Theory
Appl. Sci. 2020, 10, 5241 (1)  The clear geometric meaning of POE theory makes it an effective tool for kinematic modeling. One twist can represent the motion of one rigid body. The twist coordinates can be written as where ω is the unit vector of the rotational axis; ν is the linear velocity of the rotational motion; vector ν can be obtained as ν = p × ω; p is the position vector in the base for a point fixed on the axis; and the exponential matrix of the twist is represented as ( )  The clear geometric meaning of POE theory makes it an effective tool for kinematic modeling. One twist can represent the motion of one rigid body. The twist coordinates can be written as where ω is the unit vector of the rotational axis; ν is the linear velocity of the rotational motion; vector ν can be obtained as ν = p × ω; p is the position vector in the base for a point fixed on the axis; and the exponential matrix of the twist is represented as The clear geometric meaning of POE theory makes it an effective tool for kinematic modeling. One twist can represent the motion of one rigid body. The twist coordinates can be written as where ω is the unit vector of the rotational axis; ν is the linear velocity of the rotational motion; vector ν can be obtained as ν = p × ω; p is the position vector in the base for a point fixed on the axis; and the exponential matrix of the twist is represented as where e ωθ can be calculated by following the formula Using the POE formula, the forward kinematics of a multijoint industrial robot with n joints can be expressed as g st (θ) = e ξ 1 θ 1 e ξ 2 θ 2 · · · e ξ n θ n In order to obtain the position of the print nozzle (the end-effector) relative to the workpiece, and solve the inverse kinematics conveniently, the base coordinate system was set as the global coordinate system. The obtained kinematic POE formula was translated to the workpiece coordinate system in order to obtain the position of the print nozzle relative to the workpiece. A multijoint industrial robot can be seen as one open chain from the base to the printing nozzle. As shown in Figure 2, the order of the open chain is base→A-joint→B-joint→C-joint→D-joint→printing nozzle. Point R A is the intersection point of the rotation axis of the A-joint. The coordinates of R A in BCS are R A = −V W + V A . Point R B is the intersection point of the rotation axis of the B-joint. The coordinates relative to BCS are R B = −V W + V B . Point R C is the intersection point of the rotation axis of the C-joint. The coordinates relative to BCS are R C = −V W + V C . Point R D is the intersection point of the rotation axis of the D-joint. The coordinates relative to BCS are R D = −V W + V D . First, motion twists and exponential matrices of all axes were established according to the kinematic chain and POE theory. Unit twists of the A-joint, B-joint, C-joint, and D-joint were expressed as According to the order of the open chain, the transformation matrix of a multijoint industrial robot can be obtained with the exponential matrix of each joint. The movements of all joints relative to their local zero positions are positive relative to the global coordinate system (BCS). The POE formula of the tool relative to the base of the multijoint industrial robot was then expressed as b t T = e ξ b · e ξ 1 α · e ξ 2 β · e ξ 3 γ · e ξ 4 δ · e ξ t = a x P x n y o y a y P y n z o z a z P z 0 0 0 1 where ξ b = ξ t = [0,0,0,0,0,0] represents the twist of the base and the printing nozzle, respectively. The pose matrix is R = [n,o,a] and the position matrix is P = [P x ',P y ',P z '] T . However, it can be seen from Figure 1 that the fourth joint of the multijoint industrial robot modeled in this paper is a false joint. Due to this limitation of the mechanical structure, the end-effector is always parallel to the printing platform, meaning that δ = −β − γ can be obtained from the geometry, namely According to Equation (3), the position and pose of the multijoint industrial robot can be obtained. The last three joints of the common six-joint industrial robots are concentrated on the wrist to indicate pose and the first three joints indicate the position. Due to the limitation of the mechanical structure of the robot modeled in this paper and the existing 3D-printing methods, only three joints were considered. Thus, the forward kinematics equations of the printing nozzle relative to the workpiece were expressed as          P x = cos α[a 3 cos(β + γ) + a 2 cos β + a 1 + a 4 ] − U x P y = sin α[a 3 cos(β + γ) + a 2 cos β + a 1 + a 4 ] − U y P z = a 3 sin(β + γ) + a 2 sin β + d 1 − a 5 − U z (9) Appl. Sci. 2020, 10, 5241 6 of 20

All Possible Solutions of Joint Angle in Inverse Kinematics
For the kinematic modeling of robots established by POE theory, the Paden-Kahan subproblem is needed to solve the inverse kinematics. It requires the tool coordinate system relative to the base coordinate system. As shown in Figure 2, [P x ,P y ,P z ] T was converted to the base coordinate system and expressed as [P x ',P y ',P z '] T in the base coordinate system, where P x ' = P x + U x , P y ' = P y + U y , P z ' = P z + U z .
Firstly, ignoring the periodicity of trigonometric functions and the stroke range of the rotation axis, one rotation angle a, b, and c of the A-joint, B-joint, and C-joint was calculated as where r = P x 2 + P y 2 .
Secondly, considering the periodicity of trigonometric function and the travel range of rotation axis, according to the structure of the multijoint industrial robot, the stroke range of α was [−π, π], the stroke range of β was [−π/2, 3π/2], and the stroke range of γ was [−π, π]. The range of the inverse trigonometric function and the stroke range of the A-joint, B-joint, and C-joint are shown in Table 1. For the A-joint, when the range of α was [−π, π], the tangent function was as shown in Figure 3. Thus, the stroke range of equations were divided into four regions, and the dual solution of equations is shown in Table 2. Table 1. The range of inverse trigonometric functions and the stroke of rotation axes.

Solution of Rotation Angles Range of Inverse Trigonometric Functions Stroke of A-Joint, B-Joint, or C-Joint
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 22

Range of Inverse Trigonometric Functions
Stroke of A-Joint, B-Joint, or C-Joint Combining the solutions of joint angles of the A-joint, B-joint, and C-joint, the four solutions of (α, β, γ) were represented as follows. When P x = 0, r = (P x + U x ) 2 + P y + U y 2 , the solutions were represented as when Px 0, a = arctan P y + U y /(P When α, β, and γ are all within the feasible domain, there were four groups of solutions. As an example, when a > 0, the four groups of solutions were as follows.
For the multijoint robot shown in Figure 1, the A-joint is the first joint, B-joint is the second joint, and C-joint is the third joint. All possible groups of joint angles of one pose for this robot can be represented as follows. where q Q i represents the qth possible group of solutions of the point P i .; q θ h,i represents the joint angle of the hth joint in the qth possible group of solutions.

Multisolution Selection of Joint Angle Based on the Key Order of the Joint
Four groups of possible solutions can be obtained for one path point. The solutions should be compared and analyzed to obtain the appropriate solution for the actual control and the motion of the robot. This section proposes one multisolution selection method to obtain the appropriate group of joint angles from the four groups of possible solutions. The key order of the joints is determined by comparing the sensitivity of all joint angles. The possible solutions are selected by comparing the increment of joint angles one by one, according to the key degree of all joints.

Key Order of the Joint
When the multijoint industrial robot is running according to the specified path, the rotation or movement of different joints may have different effects on the end position, which may affect the positioning accuracy of the robot. The appropriate solutions can be selected based on the influences of the joints on the robot, which can evaluate the key degree of the joints.
Sensitivity is introduced to evaluate the influence of each joint on the end position. Generally, sensitivity refers to the ratio of independent variable change to dependent variable change in a system. For sensitivity analysis of the joint in this paper, the angles of each joint of the multijoint industrial robot were regarded as independent variables of the system, and the dependent variable was the displacement of the end of the robot in a single direction. The sensitivity component of the hth joint in one direction can be obtained by calculating the partial derivative of the angle of the hth joint with respect to the positon of the robot end-effector in this direction, that is, where P k is the position of the robot end-effector in the k direction; θ h is the joint angle of the hth joint, and m is the number of joints of the robot. Further, in order to understand the effect of rotation or movement of a single joint on the end position of the robot, the sensitivity of a joint to the end position can be expressed by S θ h , that is, where P = (P x , P y , P z ) is the path point of the robot. The sensitivity of each joint can be solved by Jacobian matrix. The Jacobian matrix J can be obtained from Equation (15).
The sensitivity can be used to evaluate the degrees of influence of the joints on the precision of the robot, so the sensitivity can be used to present the key degree of a joint. According to the sensitivity of all joints, joints are sorted from high key degree to low key degree. The sequence of each joint after sorting is the key order of the joint, which is from one to the number of the joints in the robot. For example, if the sequence of the hth joint after sorting is j, then the key order of the hth joint d h is j, and the hth joint also is jth key joint.

Multisolution Selection
Even though each solution is available for the kinematics in theory, it is necessary to select one appropriate solution of joint angles from the view of a whole tool path and precise machining. The magnitude of joint motion is directly related to the energy consumption and running time of the entire robot system, and it is also important to ensure that continuous changes in joint angles do not mutate. The incremental changes of each joint angle should be small enough that the relative motion amplitude of the joint between the two path points can be considered in path planning. The key joints and joint increments are used to select a suitable solution from all possible solutions in the inverse kinematics of a multijoint industrial robot.
For one point P i = [x i , y i , z i ] T in the trajectory of the end-effector of the robot, all possible groups of joint angles can be obtained based on Section 2. The previous point P i−1 in the trajectory and its corresponding joint angles of the robot T are used to select the appropriate solutions Q i of the point P i . The nominal increments of all joint angles from P i-1 to P i can be calculated with the Jacobian matrix by substituting in the joint angles Q i−1 of P i−1 as The nominal joint angles of point P i can be obtained with the joint angles Q i-1 as where ∆Q i represents the nominal increments of all joint angles of P i ; ∆θ k,i represents the nominal increments of the kth joint; J i−1 represents the Jacobian matrix by substituting into joint angles Q i−1 ; Q n i represents the nominal joint angles of point P i ; and θ n k,i represents the nominal joint angles of kth joint.
All joints are sorted based on their key degree, and the key order of all joints can be obtained. For example, if the key order of the hth joint d h = j, it means that the hth joint is the jth key joint. The nominal joint angles of the jth key joint for point P i can then be obtained as where θ d j,i represents the nominal joint angles of the jth key joint. One principle based on the smallest differences of the nominal angle for one joint is proposed to select the possible solutions for this joint. The possible joint angles of one joint q θ h,i in the possible groups of solutions obtained by the inverse kinematics can be compared with the nominal joint angles θ n h,i of this joint obtained by Equations (16) and (17). The differences can be expressed as q θ h,i − θ n h,i . The selection principle for one joint is to obtain the smallest differences of joint angles. In this way, the possible solutions for one joint are selected from all possible groups of solutions. The key order of one joint reflects the key degree of the influences of this joint. The possible solutions for the joints with low key order should be selected before selecting the solutions for the joints with high key order. Namely, the possible solutions for the jth key joint can be selected from the possible solutions for the (j−1)th key joint, and all possible groups of solutions obtained in Section 2 can be seen as the possible solutions for the 0th key joint.
The flow chart of the multisolution selection method based on the key order of joints for the appropriate solution of joint angles at one point in the trajectory of the robot is presented in Figure 4. First, all possible groups of joint angles of point P i are obtained based on Section 2 to serve as the possible solutions of the 0th key joint. Secondly, the Jacobian matrix is obtained by substituting in joint angles of point P i-1 . The sensitivities of all joints are calculated with the Jacobian matrix to sort the joints and obtain the key order of all joints. The nominal increments and the nominal joint angles of all joints are calculated based on Equations (16) and (17), with the Jacobian matrix. Thirdly, the possible solutions of all joints are selected according to the key order of all joints, applying the principle based on smallest differences of the nominal angles. When key order is j, the joint with the key order j is found by comparing the key orders of all joints with j in order to obtain its nominal joint angles. The differences between this nominal joint angle and the corresponding joint angles of all possible solutions for the (j−1)th key joint are calculated. The possible solutions for the jth key joint are the solutions with the smallest differences of the joint angles of this joint. When j is larger than the number of joints of the robot, the possible solutions for all joints are obtained. The possible solutions for the joint with largest key order are the appropriate solutions for point P i of the robot.
(j−1)th key joint, and all possible groups of solutions obtained in Section 2 can be seen as the possible solutions for the 0th key joint.
The flow chart of the multisolution selection method based on the key order of joints for the appropriate solution of joint angles at one point in the trajectory of the robot is presented in Figure 4. First, all possible groups of joint angles of point Pi are obtained based on Section 2 to serve as the possible solutions of the 0th key joint. Secondly, the Jacobian matrix is obtained by substituting in joint angles of point Pi-1. The sensitivities of all joints are calculated with the Jacobian matrix to sort the joints and obtain the key order of all joints. The nominal increments and the nominal joint angles of all joints are calculated based on Equations (16) and (17), with the Jacobian matrix. Thirdly, the possible solutions of all joints are selected according to the key order of all joints, applying the principle based on smallest differences of the nominal angles. When key order is j, the joint with the key order j is found by comparing the key orders of all joints with j in order to obtain its nominal joint angles. The differences between this nominal joint angle and the corresponding joint angles of all possible solutions for the (j−1)th key joint are calculated. The possible solutions for the jth key joint are the solutions with the smallest differences of the joint angles of this joint. When j is larger than the number of joints of the robot, the possible solutions for all joints are obtained. The possible solutions for the joint with largest key order are the appropriate solutions for point Pi of the robot.

Existing Boundary Extraction Method
A numerical Monte-Carlo-based method is adopted to obtain the shape and size of the workspace [32], as shown in the point cloud in Figure 5. This is a simple method to describe a workspace boundary from an engineering perspective. The principle of workspace generation is the mapping relationship between joint variables and workspace. The three joint variables are assigned the same number of random values in their respective motion ranges, and then mapped to the workspace through kinematics equations; thus, a three-dimensional workspace point cloud is formed.
A numerical Monte-Carlo-based method is adopted to obtain the shape and size of the workspace [32], as shown in the point cloud in Figure 5. This is a simple method to describe a workspace boundary from an engineering perspective. The principle of workspace generation is the mapping relationship between joint variables and workspace. The three joint variables are assigned the same number of random values in their respective motion ranges, and then mapped to the workspace through kinematics equations; thus, a three-dimensional workspace point cloud is formed. Assuming that the point cloud in the XY plane is extracted when Z = 200 mm, as shown in Figure  6, it may be found that the workspace shape described by the point cloud is not accurate enough. Figure 7 shows the boundary curve formed by extracting the boundary contour of the workspace in Figure 6. The boundary curve presents an irregular curve, which is quite different from the actual regular curve. This deviation has a negative impact not only on the shape, but also on the area and volume estimation of the workspace. This section presents a new method to obtain more accurate boundary. Assuming that the point cloud in the XY plane is extracted when Z = 200 mm, as shown in Figure 6, it may be found that the workspace shape described by the point cloud is not accurate enough. Figure 7 shows the boundary curve formed by extracting the boundary contour of the workspace in Figure 6. The boundary curve presents an irregular curve, which is quite different from the actual regular curve. This deviation has a negative impact not only on the shape, but also on the area and volume estimation of the workspace. This section presents a new method to obtain more accurate boundary.

Extraction of the Workspace Boundary Based on 3D-Layering Concept
In order to simplify the problem, the concept of 3D-printing layering was used to divide the workspace into many slices along the Z direction, so that the three-dimensional space problem can be converted to a two-dimensional plane. The set of point P is defined as the workspace of the robot, and the generalized joint variables are written as parametric equations, namely

Extraction of the Workspace Boundary Based on 3D-Layering Concept
In order to simplify the problem, the concept of 3D-printing layering was used to divide the workspace into many slices along the Z direction, so that the three-dimensional space problem can be converted to a two-dimensional plane. The set of point P is defined as the workspace of the robot, and the generalized joint variables are written as parametric equations, namely

Extraction of the Workspace Boundary Based on 3D-Layering Concept
In order to simplify the problem, the concept of 3D-printing layering was used to divide the workspace into many slices along the Z direction, so that the three-dimensional space problem can be converted to a two-dimensional plane. The set of point P is defined as the workspace of the robot, and the generalized joint variables are written as parametric equations, namely where θ i = (α i , β i , γ i ). From the kinematics Equation (5), the contour boundary equation can be derived as From the boundary equation, it can be seen that the boundary of the contour should be a circle or annulus in the XY plane, of which the center is and the radius is r = a 3 cos(β + γ) + a 2 cos β + a 1 + a 4 (22) In order to solve the minimum and maximum problems of r, γ is expressed as β by the kinematics Equation (8), namely The above expression is substituted into Equation (22), giving r = ± a 2 3 − (P z − a 2 sin β − d 1 + a 5 + U z ) 2 + a 2 cos β + a 1 + a 4 (24) As can be seen from the expression of radius r, the radius is only related to the variable β at a given height. According to the kinematics Formula (5) of the multijoint industrial robot, the relationship between β and γ can be obtained as − arctan a 3 sin γ a 2 + a 3 cos γ (25) Figure 8 is a contour of the workspace in the XZ plane. The contour is divided into four parts: FG, GM, ME, and EF. The maximum radius of the boundary of each layer of the workspace in the XY plane is on the contour segment EF, and the minimum radius is on the contour segments FG, GM, and ME. According to the structural parameters and the joint-angle stroke range of the multijoint industrial robot, the condition for forming the contour segment EF (i.e., the outer contour of the workspace) is that when the joint angle γ of the robot takes the maximum, that is, γ = γ max . β = arcsin a 3 (P z − d 1 + a 5 + U z ) a 2 2 + 2a 2 a 3 cos γ max + a 2 3 − arctan a 3 sin γ max a 2 + a 3 cos γ max (26) FG, GM, ME, and EF. The maximum radius of the boundary of each layer of the workspace in the XY plane is on the contour segment EF, and the minimum radius is on the contour segments FG, GM, and ME. According to the structural parameters and the joint-angle stroke range of the multijoint industrial robot, the condition for forming the contour segment EF (i.e., the outer contour of the workspace) is that when the joint angle γ of the robot takes the maximum, that is, γ = γmax.   (26) Similarly, the condition for forming the contour segment GM is that when the joint angle γ of the robot takes the minimum, that is, γ = γmin. Similarly, the condition for forming the contour segment GM is that when the joint angle γ of the robot takes the minimum, that is, γ = γ min .
− arctan a 3 sin γ min a 2 + a 3 cos γ min (27) The condition for forming the contour segment FG is that when the joint angle β of the robot takes the maximum, that is, β = β max . Similarly, the condition for forming the contour segment ME is that when the joint angle γ of the robot takes the minimum, that is, β = β min . Thus, the minimum and maximum of r in the XY plane can be obtained at a given height Z, as shown in Table 3. Table 3. The value of β when r takes the minimum and maximum at a given height Z.

The Scope of the P z
The Value of β when r min The Value of β when r max P z < P z (0, γ min ) (L 1 ∼ L 2 ) β min arcsin a3(Pz−d1+a5+Uz) a 2 2 +2a2a3 cos γ max +a 2 3 − arctan a3 sin γ max a2+a3 cos γ max P z (0, γ min ) ≤ P z < P z (β min , γ min ) (L 2 ∼ L 3 ) f(β, γ) = a 3 sin(β + γ) + a 2 sin β + d 1 − U z (28) The partial derivative is calculated: Through calculation, it is found that the critical point cannot be found, that is, (β, γ) cannot be found to make the above two equations simultaneously true. So the maximum and minimum of f(β, γ) has to be on the boundary of D, as follows: In conclusion, if the height Z of each layer is constant, then the motion of the robot is constrained to a plane, and accordingly, the kinematic expression of the robot degenerates to In fact, the method proposed in this paper is not only applicable to the robot used in this paper, but also applicable to other joint robots. Most of these joint robots have one common feature, that is, the first joint (some robots have other joints) is the joint with the rotation axis vertically upward. This causes the boundary of each layer of workspace to be a circle or arc. When the height Z is fixed, according to the kinematic formula, the maximum and minimum radius of the workspace at this time can be solved, so that the boundary of each layer of the workspace can be obtained.  (31) was imported into the modeling software to generate a 3D printing layering model of the entire workspace, as shown in Figure 9. The layer height was then set to 0.2mm and the 3D-printing codes for the entire workspace were exported for physical printing of the workspace. From the 3D-printing layering model, the boundary curves of each layer were obtained in the entire workspace, as shown in Figures 10 and 11.

Verification of the Extraction of Workspace Boundary Based on 3D-Layering Concept
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 according to the kinematic formula, the maximum and minimum radius of the workspace at this time can be solved, so that the boundary of each layer of the workspace can be obtained.  (31) was imported into the modeling software to generate a 3D printing layering model of the entire workspace, as shown in Figure 9. The layer height was then set to 0.2mm and the 3D-printing codes for the entire workspace were exported for physical printing of the workspace. From the 3D-printing layering model, the boundary curves of each layer were obtained in the entire workspace, as shown in Figures 10 and 11. Figure 9. The 3D-printing model of the entire workspace. Figure 9. The 3D-printing model of the entire workspace. Figure 9. The 3D-printing model of the entire workspace.  As in Section 4.1, in order to compare with other methods, especially the Monte Carlo method, the boundary curve of the workspace when Pz = 200 mm was extracted using the method in this paper, as shown in Figure 12. By comparing Figure 12 with Figure 7, it was observed that the boundary curve of the workspace extracted by Monte Carlo method was irregular, which affected the accuracy and efficiency of the solution of the workspace. The boundary curve extracted by the method proposed in this paper was closer to the real boundary and had a better fitting degree and smoothness compared to the actual workspace.

Verification of Multisolution Selection Method Based on the Key Degrees of the Joints
In order to verify the correctness of the multisolution selection method based on the key degrees of the joints, the workspace stratification model generated in Section 5.1 was used to select a layer As in Section 4.1, in order to compare with other methods, especially the Monte Carlo method, the boundary curve of the workspace when P z = 200 mm was extracted using the method in this paper, as shown in Figure 12. By comparing Figure 12 with Figure 7, it was observed that the boundary curve of the workspace extracted by Monte Carlo method was irregular, which affected the accuracy and efficiency of the solution of the workspace. The boundary curve extracted by the method proposed in this paper was closer to the real boundary and had a better fitting degree and smoothness compared to the actual workspace. As in Section 4.1, in order to compare with other methods, especially the Monte Carlo method, the boundary curve of the workspace when Pz = 200 mm was extracted using the method in this paper, as shown in Figure 12. By comparing Figure 12 with Figure 7, it was observed that the boundary curve of the workspace extracted by Monte Carlo method was irregular, which affected the accuracy and efficiency of the solution of the workspace. The boundary curve extracted by the method proposed in this paper was closer to the real boundary and had a better fitting degree and smoothness compared to the actual workspace.

Verification of Multisolution Selection Method Based on the Key Degrees of the Joints
In order to verify the correctness of the multisolution selection method based on the key degrees of the joints, the workspace stratification model generated in Section 5.1 was used to select a layer boundary for analysis. As shown in Figure 12

Verification of Multisolution Selection Method Based on the Key Degrees of the Joints
In order to verify the correctness of the multisolution selection method based on the key degrees of the joints, the workspace stratification model generated in Section 5.1 was used to select a layer boundary for analysis. As shown in Figure 12, two path points P 0 = (258.238, 12.483,200) and P 1 = (255.829, 37.333, 200) were taken from the printing path, and Q 0 = (2.7675, 99.0989, −117.2417) was known. Verification of the multisolution selection method was performed according to the following steps.
(1) According to the inverse kinematic formula in Section 2.2, four groups of solutions of joint angles at path point P 1 were solved, and then all possible groups of joint angles were obtained according to Equation (12), as shown in Table 4. They were seen as the possible solutions for the 0th key joint. (2) According to Section 3.1, the Jacobian matrix J 0 was obtainedat the path point P 0 . The sensitivity S q of each joint was calculated by the Jacobian matrix J 0 .
The key degrees of all joints were obtained and all joints were sorted according to their sensitivity. For S α > S γ > S β , the A-joint was the first key joint and its key order was one; the C-joint was the second key joint and its key order was two; the B-joint was the third key joint and its key order was three.
The nominal increments and the nominal joint angles of all joints were calculated based on Equations (16) and (17) with the Jacobian matrix. The displacement components of the end position from the path point P 0 to P 1 were calculated as ∆P x = −2.409 mm, ∆P y = 24.850 mm, and ∆P z = 0 mm. According to the Jacobian matrix J 0 and Equation (16), the nominal increments of all joint angles from P 0 to P 1 were calculated as

1.
For j = 1, the joint with the key order of 1 was found to be the A-joint. Its nominal joint angles θ n 1,1 were obtained as 8.2939 • based on Q n 1 . The differences q θ 1,1 − θ n 1,1 between θ n 1,1 and q θ 1,1 in Table 4 were calculated. The smallest difference of α-joint angles was 8.3025 • −8.2939 • =0.0086 • . The possible solutions for the first key order were 1 Q 1 and 2 Q 1 , and then, j = j + 1 = 2. The second key order was the C-joint, and its nominal joint angles θ n 3,1 were obtained as −117.65 • . The differences of θ n 3,1 were calculated with the joint angles of the C-joint in 1 Q 1 and 2 Q 1 . The possible solutions for the second key order were 1 Q 1 , which had the smallest difference of θ n 3,1 . For the third key order, the nominal joint angles were θ n 2,1 = 99.5755 and the possible solution was also 1 Q 1 . As a result, the appropriate solution of point P 1 was Q 1 = 1 Q 1 = [8.3025, 99.0992, −117.2420] T .

Experiments
The printing code of the workspace was generated with the proposed kinematics model and the multisolution method by setting V U = [0,0,0]T. The multijoint-industrial-robot-based 3D-printing system then printed its own workspace. The simulation are proposed with the printing code by setting the layer height as 10mm and 5mm, respectively as shown in Figure 13. It was found here that it was difficult to print the actual size of the workspace, because the end-effector approached singularities near the boundaries of its workspace, where it lost rigidity and precision. In addition, the layer height setting may make printing time be too long, which will waste time and printing materials.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 19 of 22 The printing code of the workspace was generated with the proposed kinematics model and the multisolution method by setting VU = [0,0,0] T . The multijoint-industrial-robot-based 3D-printing system then printed its own workspace. The simulation are proposed with the printing code by setting the layer height as 10mm and 5mm, respectively as shown in Figure 13. It was found here that it was difficult to print the actual size of the workspace, because the end-effector approached singularities near the boundaries of its workspace, where it lost rigidity and precision. In addition, the layer height setting may make printing time be too long, which will waste time and printing materials. The entire workspace was scaled by 10 times in equal proportion, so the robot could print its workspace. The origin of the base of the robot was chosen as the origin of the scaling, as shown in Figure 14. After scaling, the height dimensions became Zmax = 38.2924 mm and Zmin = −17.5947 mm. The total height was 55.8871 mm. For printing, the scaled workspace model was translated to the print platform, as shown in Figure 14. First, the model was translated from the origin of the base to the origin of the WCS with vector Va = [225, −70,12] T mm. At this time, the height of a part of the model was negative. The model was then translated with vector Vp = [70, 70, 5.5947] T mm to make the model on the print platform. The printing code with joint angles of the robot was generated based on the model and the proposed multisolution method. Vector VU = Va + Vp was used for the calculation of the all possible joint-angle groups. With the printing code, the scaled workspace was printed with this multijoint robot. The printing material was PLA. The layer height was set to 0.2 mm, and feed velocity was 50 mm/min. The printed model of the scaled workspace is shown in Figure 15. The results showed that the workspace visualization method produced the real workspace, and the proposed kinematics and multisolution method can be used for 3D printing with such a robot. The entire workspace was scaled by 10 times in equal proportion, so the robot could print its workspace. The origin of the base of the robot was chosen as the origin of the scaling, as shown in Figure 14. After scaling, the height dimensions became Z max = 38.2924 mm and Z min = −17.5947 mm. The total height was 55.8871 mm. For printing, the scaled workspace model was translated to the print platform, as shown in Figure 14. First, the model was translated from the origin of the base to the origin of the WCS with vector V a = [225, −70,12] T mm. At this time, the height of a part of the model was negative. The model was then translated with vector V p = [70, 70, 5.5947] T mm to make the model on the print platform. The printing code with joint angles of the robot was generated based on the model and the proposed multisolution method. Vector V U = V a + V p was used for the calculation of the all possible joint-angle groups. With the printing code, the scaled workspace was printed with this multijoint robot. The printing material was PLA. The layer height was set to 0.2 mm, and feed velocity was 50 mm/min. The printed model of the scaled workspace is shown in Figure 15. The results showed that the workspace visualization method produced the real workspace, and the proposed kinematics and multisolution method can be used for 3D printing with such a robot.
The proposed method can be used to obtain and visualize a detailed workspace of a multijoint industrial robot by using a 3D-printing layering concept. There are some possible industrial applications. For example, multijoint industrial robots play a great role in automatic production lines, such as in automation welding, automatic material delivery, and automatic assembly. Detailed workspace models of such robots can be a great help for the arrangement of the robots. The trajectory of the welding or the material delivery should be in the workspace of the robot. The workspace also provides a restricted area for path planning and reliable control of the robot in an automatic production line.
In addition, multijoint industrial robots can be used for 3D printing to manufacture complex freeform surfaces. The visual workspace can help to determine and optimize the locations of the printed parts. The printing trajectories are also limited by the workspace of the robot. The proposed method can be used to obtain and visualize a detailed workspace of a multijoint industrial robot by using a 3D-printing layering concept. There are some possible industrial applications. For example, multijoint industrial robots play a great role in automatic production lines, such as in automation welding, automatic material delivery, and automatic assembly. Detailed workspace models of such robots can be a great help for the arrangement of the robots. The trajectory of the welding or the material delivery should be in the workspace of the robot. The workspace also provides a restricted area for path planning and reliable control of the robot in an automatic production line. In addition, multijoint industrial robots can be used for 3D printing to manufacture complex freeform surfaces. The visual workspace can help to determine and optimize the locations of the printed parts. The printing trajectories are also limited by the workspace of the robot.

Conclusions
In this paper, a workspace visualization method of multijoint industrial robot based on 3D layering concept is proposed. First, a kinematic model of a multijoint industrial robot is established based on POE theory. According to the periodicity of trigonometric function and the stroke range of joint angles, all possible joint-angle groups of one pose are presented in detail.  The proposed method can be used to obtain and visualize a detailed workspace of a multijoint industrial robot by using a 3D-printing layering concept. There are some possible industrial applications. For example, multijoint industrial robots play a great role in automatic production lines, such as in automation welding, automatic material delivery, and automatic assembly. Detailed workspace models of such robots can be a great help for the arrangement of the robots. The trajectory of the welding or the material delivery should be in the workspace of the robot. The workspace also provides a restricted area for path planning and reliable control of the robot in an automatic production line. In addition, multijoint industrial robots can be used for 3D printing to manufacture complex freeform surfaces. The visual workspace can help to determine and optimize the locations of the printed parts. The printing trajectories are also limited by the workspace of the robot.

Conclusions
In this paper, a workspace visualization method of multijoint industrial robot based on 3D layering concept is proposed. First, a kinematic model of a multijoint industrial robot is established based on POE theory. According to the periodicity of trigonometric function and the stroke range of joint angles, all possible joint-angle groups of one pose are presented in detail.

Conclusions
In this paper, a workspace visualization method of multijoint industrial robot based on 3D layering concept is proposed. First, a kinematic model of a multijoint industrial robot is established based on POE theory. According to the periodicity of trigonometric function and the stroke range of joint angles, all possible joint-angle groups of one pose are presented in detail.
Secondly, a multisolution selection method based on the key order of the joint is proposed. According to the Jacobian matrix of the multijoint industrial robot, the sensitivity of each joint is calculated to determine the key degree of each joint. The key order of a joint with a higher key degree is lower. A principle based on the smallest nominal joint-angle difference is then established. The nominal increments and the nominal joint angles of all joints are calculated with the Jacobian matrix. The possible solutions for the joints with high key order are selected from the joints with low key order. The possible solutions for the joint with the highest key order are the appropriate joint-angle groups of the point.
Thirdly, a workspace visualization method based on a 3D-printing layering concept is proposed. The formula of each layer's contours are derived from the forward kinematics of the robot. It was proven that each layer of the workspace is a circle or ring. According to the structural parameters and the stroke range of joint angle, the minimum and maximum of the radius are obtained. The height limitations of the workspace are obtained based on the forward kinematics. According to the method of boundary extraction, the array of the boundary path points is obtained. Based on the postprocessing method proposed in this paper, the array is coded to obtain the joint-angle code of the workspace, which can be used for direct 3D printing. Finally, in order to verify the effectiveness of the proposed visualization method, an experimental multijoint-industrial-robot 3D-printing platform was used to print the whole workspace.