Experimental Research on Motion Analysis Model and Trajectory Planning of GLT Palletizing Robot

Abstract

To improve wood structure processing efficiency, a palletizing robot suitable for loading and unloading glued laminated timber (GLT) has been developed. The robot comprises a six-axis connecting rod mechanism and a sponge sucker as a grasping actuator, which can enable the intelligent automatic loading and unloading and palletizing operations for small-sized GLT. Matlab robotics was used to construct the kinematic model of the GLT loading and unloading robot. Based on Matlab and Monte Carlo methods, the robot workspace was simulated and analyzed to determine the scope of the robot workspace. Using the high-order quintic and sixtic polynomial curve interpolation method, the trajectory of wood structure parts in the process of loading and unloading operations was planned, respectively, under the two conditions of staying and not staying. Tests verified that the simulation results of the pose of the end-effector were consistent with the actual pose of the robot. The robot’s working range could be analyzed intuitively and effectively. The robot’s operation trajectory planning provides data support and a parameter basis for the automatic control and program design of a loading, unloading and palletizing robot.

1. Introduction

Glued laminated timber (GLT) is one of the important products in the ecological timber industry chain, which can be widely used in wood structure construction, furniture, decoration and bridges. The GLT industry also plays a vital role in China’s employment-demand consumption industry and the national economy [1]. Several targeted industrial development plans have been issued with the acceleration of China’s ecological civilization process and the sustainable development of low-carbon industry. Applying GLT in the innovative development of prefabricated wood structure buildings represents a new opportunity to optimize the industrial development pattern of green buildings in China. However, with the diversified demands of wood structure customization, the quantity of glulam processing exponentially increases and the processing becomes more complex, as well as having low output and efficiency [2]. The development of intelligent automatic processing technology and equipment for glulam is comparatively lagging behind, which has become the main bottleneck to the promotion and application of modern wood structure buildings. Especially in the intelligent manufacturing and production of wood structures, the loading, unloading and palletizing links of GLT are mainly manual operations, which affects the production efficiency to a large extent.

Structural glued laminated timber has good performance, diverse types and a perfectly relevant standard system, which is usually used for load-bearing components in building structures. Domestic scholars have studied its physical structural properties [3,4], drilling processing properties [5], production process [6,7] and flexible drilling control system [8], which has accelerated the application of structural glued laminated timber in wood structure buildings in China. Scholars at home and abroad have explored and studied the integrated digital information manufacturing platform and software system of GLT [9], which provides a theoretical basis for intelligent manufacturing. With the wide application and promotion of industrial robots in industrial automation production, industrial robots have been applied in the intelligent manufacturing of GLT. By analyzing the process characteristics of plate processing, the layout of equipment and the running tracks of industrial robots are planned. Machine vision is used to quickly identify the position of plates, and an automatic loading and unloading system for a woodworking center is built [10]. The intelligent automatic loading and unloading palletizing technology and equipment of GLT is a key component of the intelligent manufacturing of wood structures, which can remove the need for a large portion of manual labor from heavy repetitive loading, unloading and palletizing operations, reduce labor intensity and improve production efficiency. Therefore, the study of the trajectory planning method and motion analysis model of a palletizing robot is of great significance to realize the intelligent manufacturing of wood structures.

In this paper, a six-axis robot was used, with a sponge sucker as the integrated material to grasp the actuator, for intelligent, automatic loading and unloading to improve the automation level and work efficiency. By constructing a robot kinematic model, the workspace of the robot was analyzed using the Monte Carlo method and Matlab based on a kinematic analysis. The robot’s automatic loading and unloading operation trajectory planning provides data support and a parameter basis for the automatic control and program design of a loading, unloading and palletizing robot.

2. Overview of the GLT Palletizing Robot

GLTs are mostly wood boards and wood strips of various kinds and sizes. To meet the requirements of the integrated material loading and unloading palletizing operation, a 6-DOF, 6-axis robot with a multidimensional and comprehensive working space was equipped with a fixture. There are knots and air permeability on the surface of soft GLT, and a mechanical clamping fixture can easily cause damage to the surface of the material. A traditional vacuum suction cup cannot effectively form a confined space on the surface of the material and easily discharges during operation, which cannot sufficiently adapt to the automatic loading and unloading of GLT. Therefore, through the comprehensive analysis of various vacuum suction cup performances, principles and vacuum system control methods [11,12], the internal check valve sponge suction cup vacuum generator was chosen as the integrated material loading and unloading capture actuator.

The wooden structure palletizing robot was composed of a base, a waist connecting rod, a large arm connecting rod, an elbow connecting rod, a small arm connecting rod, a wrist connecting rod, a rotary flange and an end effector. Each connecting rod was composed of a servo motor, an encoder and a transmission mechanism. The end-effector was composed of a connecting plate, a spring, a screw and a sponge sucker. The structure is shown in Figure 1.

3. Kinematic Analysis

3.1. Establishment of Coordinate System and Determination of D-H Parameters

In order to facilitate the complex geometry of each link of the robot, it was necessary to set a link coordinate system on each link to describe the relationship of the coordinate system between these links. The relative relationship between the joints, as well as the function between the end-effector position and the operation base frame, were studied.

The link was numbered from the fixed base (link 0) of the robot. The first movable link was named link 1, and using this method, the terminal link of the robot was used as link 6. The convention was established according to Craig’s reference coordinate system [13,14], and the coordinate Z-axis and origin of each robot arm were established on the axis in front of the arm. Except for the first and sixth links, the axis of each link was normal, which is the common normal of the front and rear adjacent links. The robot’s connecting rod coordinate system is shown in Figure 2:

According to the Denavit–Hartenberg method [15], D-H parameters of each connecting rod were determined, and the D-H parameters are shown in

Table 1. Improved D-H parameter table.

Joint	θi/rad	di/m	αi−1/rad	ai−1/m	Joint Rotation Range/rad
1	0	0	pi	0	−178/pi~178/pi
2	−pi/2	0	pi/2	149.078	−98/pi~158/pi
3	0	0	0	790.385	−178/pi~78/pi
4	pi	−860.711	pi/2	150.305	−400/pi~400/pi
5	pi	0	pi/2	0	−120/pi~120/pi
6	pi	0	pi/2	0	500/pi~500/pi

.

3.2. Forward Kinematics

Four homogeneous transformations were used from the position of connecting rod i−1 to rod i. Firstly, the {i−1} coordinate system was rotated around the x_i−1 axis α_i−1 angle so that the z_i−1 axis rotated in the same direction as the z_i axis. Then, the a_i−1 distance was translated along the x_i−1 axis and the coordinate system was moved to the i axis, and then, it rotated around the z_i axis θ_i angle. Finally, the d distance of translation along the zi axis coincided with the coordinate system {i} of i bar [16,17]. The transformation from coordinate {i−1} to coordinate {i} could be realized by using Equation (1).

$${}_i{}^{i-1}T=Rot\left(x,{\alpha }_{i-1}\right)Trans\left({\alpha }_{i-1},0,0\right)Rot\left(z,{\theta }_i\right)Trans\left(0,0,d_i\right)$$

(1)

Formula (2) was used, which could be calculated from four transformation matrices.

$${}_i{}^{i-1}T=\left[\begin{array}{cccc}c{\theta }_i& -s{\theta }_i& 0& a_{i-1}\\ s{\theta }_{i}c{\alpha }_{i-1}& c{\theta }_{i}c{\alpha }_{i-1}& -s{\alpha }_{i-1}& -d_{i}s{\alpha }_{i-1}\\ s{\theta }_{i}s{\alpha }_{i-1}& c{\theta }_{i}s{\alpha }_{i-1}& c{\alpha }_{i-1}& d_{i}c{\alpha }_{i-1}\\ 0& 0& 0& 1\end{array}\right]$$

(2)

where c = cos and s = sin.

Substituting the robot D-H parameters into Equation (2), we could obtain:

$${}_1{}^{0}T=\left[\begin{array}{cccc}c{\theta }_1& -s{\theta }_1& 0& 0\\ -s{\theta }_1& -c{\theta }_1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(3)

$${}_2{}^{1}T=\left[\begin{array}{cccc}c{\theta }_2& -s{\theta }_2& 0& 149.078\\ 0& 0& -1& 0\\ s{\theta }_2& c{\theta }_2& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(4)

$${}_3{}^{2}T=\left[\begin{array}{cccc}c{\theta }_3& -s{\theta }_3& 0& 790.385\\ s{\theta }_3& c{\theta }_3& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$$

(5)

$${}_4{}^{3}T=\left[\begin{array}{cccc}c{\theta }_4& -s{\theta }_4& 0& 150.305\\ 0& 0& -1& 860.711\\ s{\theta }_4& c{\theta }_4& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(6)

$${}_5{}^{4}T=\left[\begin{array}{cccc}c{\theta }_5& -s{\theta }_5& 0& 0\\ 0& 0& -1& 0\\ s{\theta }_5& c{\theta }_5& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(7)

$${}_6{}^{5}T=\left[\begin{array}{cccc}c{\theta }_6& -s{\theta }_6& 0& 0\\ 0& 0& -1& 0\\ s{\theta }_6& c{\theta }_6& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

(8)

$${}_{ 6}{}^{0}T={}_1{}^{0}T{}_2{}^{1}T{}_3{}^{2}T{}_4{}^{3}T{}_5{}^{4}T{}_6{}^{5}T$$

(9)

The robot base had a fixed transformation to Z relative to the workpiece reference system along the Z₁ axis d_B. The end of the axial robot tool had a fixed transformation E relative to the wrist’s end coordinate {6}. Firstly, dt was translated along the Z₆ axis and then rotated around the Y_e axis θ_e [18,19]. The transformation TT of the robot tool end relative to the workpiece reference system is:

$$\mathit{TT}=Z{}_6{}^{0}TE$$

(10)

$$Z=Trans\left(0,0,d_B\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& d_B\\ 0& 0& 0& 1\end{array}\right]$$

(11)

$$E=Rot\left(Y_e,pi\right)Trans\left(0,0,d_t\right)=\left[\begin{array}{cccc}C{\theta }_e& 0& S{\theta }_e& 0\\ 0& 1& 0& 0\\ -S{\theta }_e& 0& C{\theta }_e& 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& d_t\\ 0& 0& 0& 1\end{array}\right]$$

(12)

The robot kinematics equation X is:

$$TT=Z{}_6{}^{0}TE=\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& p_x\\ r_{21}& r_{22}& r_{23}& p_y\\ r_{31}& r_{32}& r_{33}& p_z\\ 0& 0& 0& 1\end{array}\right]$$

(13)

The homogeneous transformation matrix of Equations (3)–(6) was substituted into Equation (7), and then, Equation (8) could be obtained:

$$\begin{array}{l}{r}_{11}=s_6(c_4{s}_1+s_4(c_1{s}_2{s}_3-c_1{c}_2{c}_3\left)\right)-c_6\left(c_5\right(s_1{s}_4-c_4(c_1{s}_2{s}_3-c_1{c}_2{c}_3))+s_5(c_1{c}_2{s}_3+c_1{c}_3{s}_2\left)\right)\\ r_{12}=-c_6(c_4{s}_1+s_4(c_1{s}_2{s}_3-c_1{c}_2{c}_3\left)\right)-s_6\left(c_5\right(s_1{s}_4-c_4(c_1{s}_2{s}_3-c_1{c}_2{c}_3))+s_5(c_1{c}_2{s}_3+c_1{c}_3{s}_2\left)\right)\\ r_{13}= c_5(c_1{c}_2{s}_3 +c_1{c}_3{s}_2)-s_5(s_1{s}_4-c_4(c_1{s}_2{s}_3-c_1{c}_2{c}_3\left)\right)\\ r_{21}=s_6(c_1{c}_4+s_4(c_2{c}_3{s}_1-s_1{s}_2{s}_3\left)\right)-c_6\left(c_5\right(c_1{s}_4-c_4(c_2{c}_3{s}_1-s_1{s}_2{s}_3))-s5(c_2{s}_1{s}_3+c_3{s}_1{s}_2\left)\right),\\ r_{22}=-c_6(c_1{c}_4+s_4(c_2{c}_3{s}_1-s_1{s}_2{s}_3\left)\right)-s_6\left(c_5\right(c_1{s}_4-c_4(c_2{c}_3{s}_1-s_1{s}_2{s}_3))-s_5(c_2{s}_1{s}_3+c_3{s}_1{s}_2\left)\right)\\ r_{23}=-s_5(c_1{s}_4-c_4(c_2{c}_3{s}_1-s_1{s}_2{s}_3\left)\right)-c_5(c_2{s}_1{s}_3+c_3{s}_1{s}_2)\\ r_{31}= s_{23}{s}_4{s}_6-c_6(c_{23}{s}_5-s_{23}{c}_4{c}_5)\\ r_{32}=-s_6(c_{23}{s}_5-s_{23}{c}_4{c}_5) –s_{23}{c}_6{s}_4\\ r_{33}= c_{23}{c}_5+s_{23}{c}_4{s}_5\\ p_x= a_2{c}_1+a_4{c}_{23}{c}_1+d_4{s}_{23}{c}_1+a_3{c}_1{c}_2-d_{\mathrm{t}}{s}_{23}{c}_1{c}_5+d_{\mathrm{t}}{s}_1{s}_4{s}_5+d_{\mathrm{t}}{c}_1{c}_2{c}_3{c}_4{s}_5-d_{\mathrm{t}}{c}_1{c}_4{s}_2{s}_3{s}_5\\ p_y= d_t{s}_{23}{c}_5{s}_1-a_4{c}_{23}{s}_1-d_4{s}_{23}{s}_1-a_3{c}_2{s}_1-a_2{s}_1+d_t{c}_1{s}_4{s}_5-d_t{c}_2{c}_3{c}_4{s}_1{s}_5+d_t{c}_4{s}_1{s}_2{s}_3{s}_5\\ p_z= d_b +d_4{c}_{23} - a_4{s}_{23} - a_3{s}_2 - d_t{s}_{23}{s}_{45}/2-d_t{c}_{23}{c}_5+d_t s_{4-5}{s}_{23}/2\end{array}$$

where $c_i=\cos {\theta }_i$, $s_i=\sin {\theta }_i$, $c_{ij}=\cos ({\theta }_i+{\theta }_j)$, $s_{ij}=\sin ({\theta }_i+{\theta }_j), s_{i-j}=\sin ({\theta }_i-{\theta }_j)$

4. Matlab Simulation of Robot

4.1. Model Establishment

The Matlab Robotics Toolbox is a set of Matlab-based robotics toolboxes developed and maintained by the Australian scientist, Peter Corke; it contains many important functions such as kinematics, dynamics and trajectory planning in robotics research [20]. It can conduct a graphic simulation of a robot and analyze the test data results of real robot control. The Link function in Matlab robotics was used to build the robot connecting rod, and the six connecting rods were connected by the SerialLink function to build the simulation model of the robot [21]. Industrial robot modeling procedures are as follows:

clear;clc;clf;

th(1) = 0; d(1) = 0; a(1) = 0; alp(1) = pi;

th(2) = −pi/2; d(2) = 0; a(2) = 149.078; alp(2) = pi/2;

th(3) = 0; d(3) = 0; a(3) = 790.385; alp(3) = 0;

th(4) = pi; d(4) = −860.711; a(4) = 150.305; alp(4) = pi/2;

th(5) = pi; d(5) = 0; a(5) = 0; alp(5) = pi/2;

th(6) = pi; d(6) = 0; a(6) = 0; alp(6) = pi/2;

L1 = Link([th(1), d(1), a(1), alp(1), 0], ‘modified’);

L2 = Link([th(2), d(2), a(2), alp(2), 0], ‘modified’); L2.offset = −pi/2;

L3 = Link([th(3), d(3), a(3), alp(3), 0], ‘modified’);

L4 = Link([th(4), d(4), a(4), alp(4), 0], ‘modified’);L4.offset = pi;

L5 = Link([th(5), d(5), a(5), alp(5), 0], ‘modified’);L5.offset = pi;

L6 = Link([th(6), d(6), a(6), alp(6), 0], ‘modified’);L6.offset = pi;

robot = SerialLink([L1, L2, L3, L4, L5, L6]);

robot.name = ‘GREE Robot’;

robot.display();

robot.base = transl(0,0,449.822);

robot.tool = transl(0 ,0 ,-141.068)*troty(pi);

robot.teach();

The 3D simulation model of the robot is shown in Figure 3:

4.2. Forward Kinematics Simulation Verification

Within the allowable range of joint variables, three groups of joint variables were randomly set, and the robot model established in Matlab was used for the simulation to obtain the position (X, Y, Z) and attitude (R, P, Y) of the end effector. Three sets of variable values were calculated into the industrial robot equation, and the end effector’s position (X, Y, Z) and attitude (R, P, Y) were calculated. The forward kinematics equation was verified by comparing the three groups of variables. The comparison results of the kinematics equation and simulation calculation of the end effector’s attitude are shown in

Table 2. The posture teaching results and simulation comparison results of the end effector.

Group Number	Method	θ1/°	θ2/°	θ3/°	θ4/°	θ5/°	θ6/°	Result
								X/mm	Y/mm	Z/mm	R/°	P/°	Y/°
1	equation	15	30	45	60	90	160	783.71	−336.47	323.58	100.30	−43.68	90.22
	simulation							783.71	−336.47	323.58	100.30	−43.68	90.22
2	equation	90	60	40	25	15	300	−15.43	−775.88	−157.17	−156.35	−6.28	125.75
	simulation							−15.43	−775.88	−157.17	−156.35	−6.28	125.75
3	equation	0	107	−110	320	55	10	1842.03	74.28	329.72	−138.57	37.27	−12.42
	simulation							1842.03	74.28	329.72	−138.57	37.27	−12.42

.

4.3. Inverse Kinematics Simulation Verification

The inverse kinematics of a robots is a method to solve the motion equation. The desired posture was achieved by giving the geometric parameters of the robot link and the desired position and attitude (posture) of the robot’s end effector relative to the reference coordinate system and the joint variables. When the work task was described by the Cartesian coordinate system, the above provisions had to be transformed into a series of joint variables that the arm could drive. Therefore, the solution of joint variables to determine the position and posture of the arm was the solution of the motion equation. The problem of solving the robot motion equation, inverse kinematics, was problem synthesis. The transformation matrix T of the pose of the robot’s end effector corresponding to the reference coordinate system in Table 2 was verified via inverse kinematic simulation using the ikine function qi=ikine (robot, TT) in Matlab robotics. The triple transformation matrix is provided as Equations (14)–(16).

$$T{T}_1=\left[\begin{array}{cccc}-0.4670& -0.5522& -0.6907& 783.71\\ 0.3022& 0.6344& -0.7115& -336.47\\ 0.8310& -0.5410& -0.1294& 323.58\\ 0& 0& 0& 1\end{array}\right]$$

(14)

$$T{T}_2=\left[\begin{array}{cccc}-0.5808& -0.8067& -0.1094& -15.43\\ -0.7690& 0.4996& 0.3987& -775.88\\ -0.2670& 0.3157& -0.9105& -157.17\\ 0& 0& 0& 1\end{array}\right]$$

(15)

$$T{T}_3=\left[\begin{array}{cccc}0.7771& 0.1712& 0.6056& 1842.03\\ -0.2301& -0.8184& 0.5265& 74.28\\ 0.5858& -0.5485& -0.5966& 329.72\\ 0& 0& 0& 1\end{array}\right]$$

(16)

The simulation results of the inverse kinematic equations are shown in

Table 3. Results of inverse kinematics simulation.

Group Number	The End Effector Sits in the Base Coordinate Frame						Result
	X/mm	Y/mm	Z/mm	R/°	P/°	Y/°	θ1/°	θ2/°	θ3/°	θ4/°	θ5/°	θ6/°
1	783.71	−336.47	323.58	100.30	−43.68	90.22	15	30	45	60	90	200
2	−14.43	−775.83	−157.17	−156.35	−6.28	125.75	90	60	40	25	15	300
3	1842.03	74.28	329.72	−138.57	37.27	−12.42	0	107	−110	320	55	10

:

The rotation angles of each joint in the three groups of positions and postures of the robot’s end effector obtained from the simulation were consistent with the results of the joint variable parameters assumed in the forward kinematics simulation in Table 2, which validated the robot kinematics model.

5. Simulation Verification and Solution of Working Space

A robot’s working space is the geometric locus of all possible coordinates of the center point of the end effector, which is one of the main indicators of a robot’s design. It reflects the working ability when performing a specific task, the working space range of the attitude of the end effector, and also the basis for the robot’s position arrangement in the production line and the operator’s safe working space to avoid being injured by the robot [22,23,24].

According to Formula (4), the position point coordinates of the robot’s end effector [px,py,pz]^T can be obtained, which are only related to the relevant kinematic parameter, θ_i, and the point set of the robot’s working space is P.

$$P=\left\{\left(x,y,z\right)\left|\begin{array}{c}x=f_x\left({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6\right)\\ y=f_y\left({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6\right)\\ z=f_z\left({\theta }_1,{\theta }_2,{\theta }_3,{\theta }_4,{\theta }_5,{\theta }_6\right)\end{array}\right.\begin{array}{c}{\theta }_{min}\le {\theta }_i\le {\theta }_{max}\\ i=1,2,3,4,5,6\end{array}\right\}$$

(17)

After randomly selecting within the rotation range of each joint of the robot θ_i (i = 1,2,3,4,5,6) and substituting it into Formula (6), the working space point set of the center point of the robot’s end effector could be obtained.

The calculation of the robot kinematics equation was large, and the design and analysis were complicated. In order to solve this problem, Matlab robotics 2016b software was used for robot simulation modeling. According to the Monte Carlo method, the random variables of each joint were generated as well as the group of joint space vectors of the manipulator. Through calculating the forward kinematics solution, the working space distribution cloud map was drawn from the joint space to the working space at the end (Cartesian coordinate system) [25,26].

According to the Monte Carlo numerical analysis method, 30,000 groups within each joint range θi (i = 1,2,3,4,5,6) value were randomly selected by the Matlab rand function and then substituted into Equation (6) in turn to carry out 30,000 iterative operations. After that, the corresponding coordinate set of the central point of the corresponding end effector was obtained. The point set formed by drawing all points in the three-dimensional coordinates was the robot working space cluster.

The Matlab program code is as follows:

clear;clc;clf;

th(1) = 0; d(1) = 0; a(1) = 0; alp(1) = pi ;

th(2) = −pi/2; d(2) = 0; a(2) = 149.078; alp(2) = pi/2;

th(3) = 0; d(3) = 0; a(3) = 790.385; alp(3) = 0;

th(4) = pi; d(4) = −860.711; a(4) = 150.305; alp(4) = pi/2;

th(5) = pi; d(5) = 0; a(5) = 0; alp(5) = pi/2;

th(6) = pi; d(6) = 0; a(6) = 0; alp(6) = pi/2;

L1 = Link([th(1), d(1), a(1), alp(1), 0], ‘modified’); L1.qlim = [−178,178]/180*pi;

L2 = Link([th(2), d(2), a(2), alp(2), 0], ‘modified’); L2.offset = −pi/2; L2.qlim = [−98,158]/180*pi;

L3 = Link([th(3), d(3), a(3), alp(3), 0], ‘modified’); L3.qlim = [−178,78]/180*pi;

L4 = Link([th(4), d(4), a(4), alp(4), 0], ‘modified’);L4.offset = pi; L4.qlim = [−400,400]/180*pi;

L5 = Link([th(5), d(5), a(5), alp(5), 0], ‘modified’);L5.offset = pi; L5.qlim = [−120,120]/180*pi;

L6 = Link([th(6), d(6), a(6), alp(6), 0], ‘modified’);L6.offset = pi; L6.qlim = [−500,500]/180*pi;

robot = SerialLink([L1, L2, L3, L4, L5, L6]);

robot.name = ‘GREE Robot’;

robot.display();

robot.base = transl(0,0,449.822);

robot.tool = transl(0,0,−141.068)*troty(pi);

robot.teach();

num = 30,000;

P = zeros(num,3);

for i = 1:num

q1 = L1.qlim(1) + rand*(L1.qlim(2)−L1.qlim(1));

q2 = L2.qlim(1) + rand*(L2.qlim(2)−L2.qlim(1));

q3 = L3.qlim(1) + rand*(L3.qlim(2)−L3.qlim(1));

q4 = L4.qlim(1) + rand*(L4.qlim(2)−L4.qlim(1));

q5 = L5.qlim(1) + rand*(L5.qlim(2)−L5.qlim(1));

q6 = L6.qlim(1) + rand*(L6.qlim(2)−L6.qlim(1));

q = [q1 q2 q3 q4 q5 q6];

T = robot.fkine(q);

P(i,:) = transl(T);

end

plot3(P(:,1),P(:,2),P(:,3),’b.’,’markersize’ ,2);

xlabel(‘X/m’);

ylabel(‘Y/m’);

zlabel(‘Z/m’);

hold on;

The cloud map of the robot’s working space is shown in Figure 4:

6. Robot Planning for Loading and Unloading Path of GLT

Robot trajectory planning was used to plan the trajectory of the robot’s motion and generate the operation trajectory in joint space and Cartesian space based on the kinematics and dynamics of the manipulator. The so-called trajectory refers to the displacement, velocity and acceleration of the manipulator’s joints in the movement. Trajectory planning calculated the expected trajectory according to the task’s requirements. During the robot trajectory planning, the task, motion path and trajectory of the robot were firstly described. Then, the target pose of the gripper was given to determine the path point, duration, motion speed and other trajectory parameters to reach the target. Finally, for the trajectory described internally, the robot movement’s displacement, velocity and acceleration were calculated in real time to generate the trajectory. The manipulator is usually regarded as the movement of the tool coordinate system relative to the working coordinate system. The point mentioned in trajectory planning refers to the position and attitude of the manipulator’s tool coordinate system, which was used to describe the state of the starting position, the path position and the target position of the robot during the operation process. There are two common trajectory planning methods of manipulators: in the first method, the position, velocity and acceleration on the selected transition node (interpolation point) is given a set of explicit constraints to achieve the continuity and smoothness of the operation, the use of polynomial curve function, the node interpolation and meet the constraints. In the second method, the analytical expression of the motion path is given, such as the straight path in Cartesian coordinate space, and then, a trajectory is determined in joint space or Cartesian coordinate space to approximate the predetermined path. In the first method, constraint setting and trajectory planning are carried out in joint space; so, a collision with obstacles may occur. The path constraints of the second method are given in Cartesian coordinate space, while the joint driver is controlled in joint space. Therefore, to obtain a trajectory close to the given path, the Cartesian coordinate path constraints must be first converted into joint coordinate path constraints via the function approximation method, and then, the parameterized path should be determined and the joint path constraints should be satisfied.

6.1. The Robot Makes Route Planning for the Loading and Unloading Operation of GLT

The robot’s loading and unloading trajectory of GLT refers to the displacement, velocity and acceleration of each joint of the robot when the end effector (sponge sucker) moved the GLT from the grasping position to the unloading position. The loading and unloading trajectory planning calculated the expected motion trajectory based on the analysis of dynamics when the robot was loading and unloading so that the robot could complete the loading and unloading operation continuously, smoothly and accurately. Because the loading and unloading operation position of GLT is a point-to-point operation, there was no more requirement for the trajectory in the process; so, the trajectory planning was carried out in the joint space [27]. When the robot was loading and unloading, it grabbed GLT from the tray and then placed it on the conveyor line through the critical path. In the process of movement, due to the height difference between the tray and the conveying line and the large size of GLT, GLT could have easily collided with the conveying line. It was necessary to plan the robot’s loading and unloading operation path. The robot grabbed the GLT from the tray at point A and then moved along curve AB to point B. The schematic diagram of the path planning is shown in Figure 5. At point B, the position of GLT was higher than the conveying line, and there was no collision between the GLT edge and the conveying line. After passing point B, the GLT was placed at point C of the conveyor line along the BC curve to complete the loading operation of grabbing it from the tray and placing it on the conveyor line. In the unloading operation of GLT, the robot held the integrated material at point C and moved along the CBA curve to point A to complete the feeding operation. The loading operation was the same as the unloading operation path, but the direction of movement was opposite. The trajectory planning method was the same. Therefore, this paper only carried out the trajectory planning for the feeding operation.

The coordinates of the path’s key points A, B and C are shown in

Table 4. Coordinates of the path’s key points.

Key Points	End-Effector Coordinates in the Base Coordinate System						Joint Coordinate
	X/mm	Y/mm	Z/mm	R/°	P/°	Y/°	θ1/°	θ2/°	θ3/°	θ4/°	θ5/°	θ6/°
Target point A	0	−1280	55	90	0	180	90	53.75	11.77	0	24.48	0
Way point B	640	−1108.5	706	60	0	180	60	23.88	7.89	0	58.23	0
Target point C	1280	0	616	0	0	180	0	26.59	10.52	0	52.89	0

.

6.2. Loading Operation Trajectory Planning

Common robot trajectory planning currently includes trapezoidal programming, parabolic programming, polynomial programming and S-type programming [28]. Quintic polynomial programming is most commonly used in polynomial programming and is widely used in robot trajectory planning. The acceleration discontinuity of trapezoidal programming and S-type programming will have a mechanical impact on the robot. In contrast, quintic polynomial programming is relatively mature in the transition planning method, which can ensure the continuity of displacement, velocity, acceleration and acceleration during the robot movement [29]. Therefore, the high-order quintic polynomial interpolation trajectory method was used to complete the trajectory planning. Not only was the tracking accuracy of the planned trajectory ensured, but it also could ensure continuous smooth displacement, velocity and acceleration in the process of moving. Finally, the machine ran smoothly, reducing the impact vibration and effectively extending the machine’s service life [30]. However, when the quintic polynomial was used to interpolate the points through multiple paths, some acceleration points still had the problem of abrupt change, which could have easily caused vibration. Therefore, the high-order sixth-degree polynomial interpolation can be used to realize continuous flat smoothing with extra speed in multipoint piecewise path planning.

According to the actual work requirements of the loading and unloading operations, there are two main situations. One is that when the sponge sucker absorbs the aggregate material during the loading and unloading operation, it is necessary to stay at the passing point to avoid obstacles and align the loading and unloading position. That is, the speed and acceleration are zero at the critical nodes of each path. Second, when the sponge sucker absorbs the aggregate material for the loading and unloading operation, avoiding obstacles and aligning the taking and discharging position is unnecessary. Each passing point does not need to stay. That is, the speed and acceleration are not zero at the critical node of each path. In this paper, the quintic polynomial interpolation method and high-order sixtic polynomial interpolation method were used to study the GLT loading and unloading operation trajectory planning, respectively.

6.2.1. Loading Trajectory Planning Using High-Order Quintic Polynomial Interpolation

Because the sponge sucker will stay at the trajectory path point when gripping GLT for loading operation, the AB and BC sections can be independently planned using the quintic polynomial interpolation method for the loading and unloading operation path [31]. The general formula of trajectory planning for quintic polynomial interpolation is:

$${\theta }_i\left(t\right)=a_0+a_{1}t+a_2{t}^2+a_3{t}^3+a_4{t}^4+a_5{t}^5$$

(18)

The starting point of the feeding operation path is A, the passing point is B, and the target point is C. The pose state of the three points is shown in Table 4. Suppose the starting point position of each joint i is ${\theta }_A$, the path point position is ${\theta }_B$, the target point position is ${\theta }_C$, the starting point velocity is $v_A$, the path point velocity is $v_B$, the target point velocity is $v_c$, the starting point acceleration is $a_A$, the path point acceleration is $a_B$, the target point acceleration is $a_C$, the starting point time is $t_A$, the path point time is $t_B$ and the target point time is $t_C$. The kinematic parameters of each joint i in AB and BC are shown as Equations (19) and (20).

$$\left\{\begin{array}{c}{\theta }_A=a_{10}+a_{11}{t}_A+a_{12}{t}_A^2+a_{13}{t}_A^3+a_{14}{t}_A^4+a_{15}{t}_A^5\\ v_A=a_{11}+2a_{12}{t}_A+3a_{13}{t}_A^2+4a_{14}{t}_A^3+5a_{15}{t}_A^4\\ a_A=2a_{12}+6a_{13}{t}_A+12a_{14}{t}_A^2+20a_{15}{t}_A^3\\ {\theta }_B=a_{10}+a_{11}{t}_B+a_{12}{t}_B^2+a_{13}{t}_B^3+a_{14}{t}_B^4+a_{15}{t}_B^5\\ v_B=a_{11}+2a_{12}{t}_B+3a_{13}{t}_B^2+4a_{14}{t}_B^3+5a_{15}{t}_B^4\\ a_B=2a_{12}+6a_{13}{t}_A+12a_{14}{t}_A^2+20a_{15}{t}_A^3\end{array}\right.$$

(19)

$$\left\{\begin{array}{c}{\theta }_B=a_{20}+a_{21}{t}_B+a_{22}{t}_B^2+a_{23}{t}_B^3+a_{24}{t}_B^4+a_{25}{t}_B^5\\ v_B=a_{21}+2a_{22}{t}_B+3a_{23}{t}_B^2+4a_{24}{t}_B^3+5a_{25}{t}_B^4\\ a_B=2a_{22}+6a_{23}{t}_A+12a_{24}{t}_A^2+20a_{25}{t}_A^3\\ {\theta }_C=a_{20}+a_{21}{t}_C+a_{22}{t}_C^2+a_{23}{t}_C^3+a_{24}{t}_C^4+a_{25}{t}_C^5\\ v_c=a_{21}+2a_{22}{t}_C+3a_{23}{t}_C^2+4a_{24}{t}_C^3+5a_{25}{t}_C^4\\ a_C=2a_{22}+6a_{23}{t}_C+12a_{24}{t}_C^2+20a_{25}{t}_C^3\end{array}\right.$$

(20)

It is known that ${\theta }_A$,${\theta }_B$,${\theta }_C$,$v_A=v_B=v_c=0$,$a_A=a_B=a_C=0$. The time from point A to point B and from point B to point C is 2 s. By substituting each constraint into Equations (8) and (9), we can find that the quintic multiple of AB and BC is ${\theta }_1\left(t\right)$, ${\theta }_2\left(t\right)$. The specific expression is shown in Equation (21).

$$\left\{\begin{array}{c}{\theta }_1\left(t\right)={\theta }_A+\frac{5\left({\theta }_B-{\theta }_A\right)}{4}{t}^3-\frac{15\left({\theta }_B-{\theta }_A\right)}{16}{t}^4+\frac{3\left({\theta }_B-{\theta }_A\right)}{16}{t}^5 \\ {\theta }_2\left(t\right)={\theta }_B+\frac{5\left({\theta }_C-{\theta }_B\right)}{4}{t}^3-\frac{15\left({\theta }_C-{\theta }_B\right)}{16}{t}^4+\frac{3\left({\theta }_C-{\theta }_B\right)}{16}{t}^5\end{array}\right.$$

(21)

We substituted the coordinates and positions of the three key points of A, B and C of the loading operation into the quintic polynomial (21) to obtain the feeding trajectory and the displacement, velocity and acceleration curves of the six joints, as shown in Figure 6 and Figure 7. During the loading operation, the curve from points A to B and C was continuous and smooth without abrupt changes, but the entire ABC curve was continuous at point B, but not smooth. This is because during the loading process, it stayed at point B and then moved to point C, and the trajectory was planned by two different quintic polynomials in two segments.

As shown in Figure 7, starting from point A, the displacement of the six joints of the robot increased steadily with time and gradually tended to be stable, and the range of joint displacement was all within the action range of the robot. The velocity of the robot’s six joints accelerated from t = 0 s, reached the peak velocity at t = 1 s and t = 3 s, and then decreased to 0 at t = 2 s and t = 4 s. Each joint’s acceleration and deceleration time were equal, and the range of the velocity change was within the range of the speed of each joint. The acceleration of each robot joint first increased from 0 to the maximum value, then decreased to 0, then continuously decreased from 0 to the minimum acceleration, and finally increased to 0. The running time of the increasing and decreasing positive and negative acceleration stages was equal.

Although the ABC operation track curve was not smooth at point B in the three-dimensional space during the loading process, the displacement, velocity and acceleration curves of the robot’s six joints were continuous, smooth and flexible without sudden changes, jumps, discontinuities and other conditions. This indicated that the robot operated smoothly during loading, without adverse conditions such as impact and large vibration. Therefore, the higher-order quintic polynomial interpolation method applied to the loading operation trajectory planning under the condition of staying at path point B.

6.2.2. Loading and Unloading Trajectory Planning Using High-Order Six-Order Polynomial Interpolation

During the loading operation of the sponge sucker absorbing GLT, it did not stay at the passing point B. The AB section and BC section could be unified, and the six-order polynomial interpolation method could be used for the loading operation path planning. The general formula of six-degree polynomial interpolation trajectory planning is formula (22):

$${\theta }_i\left(t\right)=a_0+a_{1}t+a_2{t}^2+a_3{t}^3+a_4{t}^4+a_5{t}^5+a_6{t}^6$$

(22)

We substituted the constraints, such as ${\theta }_A$,${\theta }_B$,${\theta }_C$,$v_A=v_c=0$,$v_B$,$a_A=a_C=0$,$a_B$, and the time from point A to point B and from point B to point C into Equation (22), as shown in (23).

$$\left\{\begin{array}{c}{\theta }_A=a_0\\ v_A=a_1\\ a_A=2a_2 \\ {\theta }_B=a_0+a_1{t}_B+a_2{t}_B^2+a_3{t}_B^3+a_4{t}_B^4+a_5{t}_B^5+a_6{t}_B^6\\ {\theta }_C=a_0+a_1{t}_C+a_2{t}_C^2+a_3{t}_C^3+a_4{t}_C^4+a_5{t}_C^5+a_6{t}_C^6\\ v_c=a_1+2a_2+3a_3{t}_C^2+4a_4{t}_C^3+5a_5{t}_C^4+6a_6{t}_C^5\\ a_C=2a_2+6a_3{t}_C+12a_4{t}_C^2+20a_5{t}_C^3+30a_6{t}_C^4\end{array}\right.$$

(23)

Equation (23) can be written in the following matrix form, Equation (24):

$$\left[\begin{array}{c}{\theta }_A\\ v_A\\ a_A\\ {\theta }_B\\ {\theta }_C\\ v_c\\ a_C\end{array}\right]=\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 2& 0& 0& 0& 0\\ 1& t_B& t_B^2& t_B^3& t_B^4& t_B^5& t_B^6\\ 1& t_C& t_C^2& t_C^3& t_C^4& t_C^5& t_C^6\\ 0& 1& 2t_C& 3t_C^2& 4t_C^3& 5t_C^4& 6t_C^5\\ 0& 1& 2& 6t_C& 12t_C^2& 20t_C^3& 30t_C^4\end{array}\right]\times \left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6\end{array}\right]$$

(24)

Equation (24) can be simplified as follows:

$$\left[\theta \right]=\left[M\right]\times \left[a\right]$$

(25)

Equation (25) can be obtained from Equation (26):

$$\left[a\right]={\left[M\right]}^{-1}\times \left[\theta \right]$$

(26)

All sixth-degree polynomial coefficients were calculated by calculating [M]⁻¹.The coordinates and positions of the three key points of A, B and C of the feeding operation were substituted into the sixth-degree polynomial (10) to obtain the feeding trajectory and the displacement, velocity and acceleration curves of the six joints, as shown in Figure 8 and Figure 9.

Figure 8 shows that the whole process of the three-dimensional space trajectory during the loading operation from point A to point B to point C was continuous and smooth without sudden changes.

It can be seen from Figure 9 that during the whole loading operation from point A to point C through point B, the displacement, velocity and acceleration of the robot’s six joints were within the range of the robot’s motion constraints. The displacement, velocity and acceleration curves were continuous, smooth and flexible without sudden changes, jumps, discontinuities or other conditions. The results show that the robot ran smoothly without adverse impact and large vibrations.

From the analysis of Figure 8 and Figure 9, it can be concluded that the ABC operation track curve in the three-dimensional space was continuous and smooth during the whole loading process. The robot’s six joints’ displacement, velocity and acceleration curves were continuous, smooth and flexible, without sudden changes, jumps, discontinuities or other conditions, indicating that the robot operated smoothly during the loading process. The high-order six-order polynomial interpolation method can be applied to loading operation trajectory planning without stopping at path point B in stable conditions.

7. Test Verification

To verify the rationality of the robot model establishment and working space simulation analysis method, as well as the correctness of the simulation analysis results, a comparison test was carried out between the robot working space and the measured results. The test site is shown in Figure 10.

Ten groups were randomly selected within the rotation range of each joint axis of the robot Angle θi. The attitude value of the corresponding end effector was obtained through the robot measurement and robot simulation model. Finally, the accuracy of the comparison results was analyzed. The comparison test results are shown in

Table 5. The results of posture teaching and the end effector’s simulation comparison.

Group Number	Method	θ1/°	θ2/°	θ3/°	θ4/°	θ5/°	θ6/°	Result
								X/mm	Y/mm	Z/mm	R/°	P/°	Y/°
1	Test verification	0	−13.96	0	0	82.84	0	808.25	0	1438.84	0	21.13	180
	Simulation							808.26	0	1438.78	0	21.12	180
2	Test verification	145.19	−39.48	−39.40	9.30	90.46	−19.69	150.00	132.04	1905.55	130.60	−62.38	115.30
	Simulation							150.02	132.07	1905.46	130.60	−62.38	115.31
3	Test verification	−171.33	36.61	−39.40	9.30	90.46	−19.69	−1458.19	245.48	1137.23	−151.30	−0.88	−170.46
	Simulation							−1458.12	245.40	1137.198	−151.30	−0.87	−170.46

. The attitude values of the corresponding end effectors obtained by the robot measurement and the robot simulation model were consistent. This means that the established robot model was correct and could be used for the kinematics simulation analysis of the robot during the loading and unloading operation of GLT. The robot working space obtained via simulation analysis was consistent with the actual working space of the robot. The scope of the working space could meet the requirements for the palletization of GLT loading and unloading, providing a theoretical basis for the robot’s automatic loading and unloading operation trajectory planning and automatic control research.

8. Conclusions

(1)

Aiming to solve the problems of low automation levels of loading and unloading, palletizing and high labor cost in the process of wood structure production, this paper adopted a six-axis robot with a sponge sucker as the grasping actuator to carry out research on the intelligent, automatic loading and unloading of GLT, which can improve the intelligent manufacturing level of wood structures. The research conclusions are as follows:

(2)

According to the Craig reference coordinate system establishment agreement, the robot link coordinate system was established, the D-H parameters of each link were determined according to the Denavit–Hartenberg method, and the mathematical model was established. The forward kinematics of the robot was deduced and solved. The kinematics model of GLT loading and unloading robot was built using Matlab robots, and the simulation results were verified using simulation technology, proving the motion analysis’s correctness.

(3)

In the Matlab environment, the working space of robot was simulated and analyzed via the Monte Carlo method, and the scope of robot working space was determined. The analysis results show that the robot has a flexible structure and a wide range of working space, which can meet the requirements for palletizing and loading GLT and provide a theoretical basis for the robot’s automatic loading and unloading operation trajectory planning and automatic control research.

(4)

The simulation results show that the method of high-order quintic and sixth-order polynomial curve interpolation can be used to plan the trajectory of wood structure parts in the process of a loading and unloading operation under the two conditions of staying and not staying. It can ensure the tracking accuracy of the trajectory. At the same time, it can also provide continuous and stable movement without impacts, large vibrations or other adverse conditions.

A verification test was carried out on the test platform of the palletizing robot for laminated material loading and unloading. The results show that the kinematics model of the robot for structural part loading and unloading was correct, and the analysis method was feasible. The robot working space drawn using the Monte Carlo method can intuitively and effectively analyze the robot’s working range.

(5)

The robot’s automatic loading and unloading operation trajectory planning provides data support and parameter basis for the automatic control and program design of the loading, unloading and palletizing robot.

(6)

This palletizing robot is only suitable for loading and unloading small-sized GLT. Large-sized laminated materials require a truss-type manipulator for loading, unloading and palletizing operations.

(7)

Further research will be carried out in loading and unloading experiments according to the results regarding the robot’s automatic loading and unloading operation trajectory planning to prove it is correct and effective.