# Dynamic Simulation Modeling of Industrial Robot Kinematics in Industry 4.0

Authors: Guang Jin; Shuai Ma; Zhenghui Li
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

This paper studies the kinematic dynamic simulation modeling of industrial robots in the Industry 4.0 environment and guides the kinematic dynamic simulation modeling of industrial robots in the Industry 4.0 environment in the context of the research. To address the problem that each parameter error has different degrees of influence on the end position error, a method is proposed to calculate the influence weight of each parameter error on the end position error based on the MD-H error model. The error model is established based on the MD-H method and the principle of differential transformation, and then the function of uniform variation of six joint angles with timet is constructed to ensure that each linkage geometric parameter is involved in the motion causing error accumulation. Through the analysis of the robot marking process, the inverse solution is optimized for multiple solutions, and a unique engineering solution is obtained. Linear interpolation, parabolic interpolation, polynomial interpolation, and spline curve interpolation are performed on the results after multisolution optimization in the joint angle, and the pros and cons of various interpolation results are analyzed. The trajectory planning and simulation of industrial robots in the Industry 4.0 environment are carried out by using a special toolbox. The advantages and disadvantages of the two planning methods are compared, and the joint space trajectory planning method is selected to study the planning of its third and fifth polynomials. The kinetic characteristics of the robot were simulated and tested by experimental methods, and the reliability of the simulation results of the kinetic characteristics was verified. The kinematic solutions of industrial robots and the results of multisolution optimization are simulated. The methods, theories, and strategies studied in this paper are slightly modified to provide theoretical and practical support for another dynamic simulation modeling of industrial robot kinematics with various geometries.


---

## Body

## 1. Introduction

In many countries around the world, courses in industrial robotics have been offered at the undergraduate and graduate levels to improve the level of research in industrial robotics. In the development and study of industrial robots operating in uncertain environments, researchers have been exposed to the core of industrial robotics and have come to understand its importance [1]. As industrial robots in the Industry 4.0 environment continue to progress and develop, they should be made more intelligent and more automated to replace people in complex, boring, and harsh work environments to maximize their potential roles and advantages. In the past, the development and manufacturing process of industrial robots was often conceptual design followed by actual production prototypes, then testing of prototypes, and then continuous modification and testing until the design requirements were met, which took a long time and increased production costs and could not meet the current production requirements [2]. In this situation and to achieve industrial development needs to improve product quality, shorten the research and development cycle, reduce research and development costs, and increase the technical content of various requirements. Virtual prototyping technology is proposed in the application of computer technology to the design and development process. Industrial robot simulation technology provides an effective experimental means for modeling and simulation of industrial robot kinematics and dynamics and is a convenient tool for designing and analyzing industrial robot performance [3]. The development and advancement of industrial robot simulation technology have greatly improved the development of industrial robots, which can save a lot of manpower and money and help improve the structural design of industrial robots. Industry 4.0 focuses on the transformation from automation to intelligence. Industry 4.0 incorporates big data into management forecasting and management decision-making, making big data part of management, improving the accuracy of management decision-making, and making management more scientific [4, 5].For the use of industrial robots in the field of welding, the industrial robots are generally programmed by using the industrial robot trainer to fix the workpiece to be welded on the corresponding fixture, then using the trainer function of the industrial robot to train the welding trajectory of the industrial robot in turn, and adjusting the trainer points according to the results of the trainer to optimize the trajectory of the industrial robot. This method can only provide a general tracking of the point to be welded, but not the intermediate trajectory, and cannot accurately control the position of the end of the robot. If the workpiece to be welded changes, it must be retrained, so the method is not very versatile. In addition, the teaching is done by workers according to the path to be taken by the industrial robot, which increases human work labor, and each person has different evaluation criteria, so the industrial robot itself does not achieve intelligence. Industrial robots are designed and used to improve the flexibility of the manufacturing process, but the motion process of industrial robots is realized based on the programming of the demonstrator, so the use of industrial robots will no longer have the advantage of flexibility for variable environments [6]. In this paper, we use offline programming technology to import the target poses into the algorithm of the inverse kinematic solution of the industrial robot and drive the motion of the industrial robot body directly by obtaining the joint angle-time sequence of the industrial robot to be tracked trajectory to maximize the flexible manufacturing capability of the industrial robot [7].This paper firstly gives a brief introduction to the development history and research status of industrial robots in the Industry 4.0 environment, analyzes the technological gap at home and abroad, briefly explains the source and significance of the subject research, and designs a six-degree-of-freedom industrial robot in Industry 4.0 environment with a maximum load capacity of 8 kg based on the research on industrial robots in Industry 4.0 environment by itself [8]. The structural analysis was carried out based on the existing model of industrial robots in the Industry 4.0 environment. Section 1 analyzes the background and significance of industrial robot research in the Industry 4.0 environment and summarizes the research framework of this paper. Section 2 analyzes the development trend of industrial robots at home and abroad as well as introduces the research on kinematic dynamic simulation modeling of industrial robots and reviews the status of domestic and international research on industrial robot calibration technology. In Section 3, the M D-H kinematic model of the industrial robot under the environment of Industry 4.0 of IRB-1410 is established based on the principles of spatial positional description and coordinate system transformation, and the end positional expressions are obtained by theoretical derivation of the positive kinematics of the industrial robot, and then the inverse solutions are derived and calculated by using the closed solution method, and the shortest travel principle is used to arrive at the best solution. In Section 4, we carry out the kinematics simulation analysis of the parameter redundancy problem of the industrial robot model under the Industry 4.0 environment. Experimental validation is performed for the multimode calibration method, the position error prediction compensation method, the trajectory correction method, and the comprehensive accuracy improvement strategy of the industrial robot proposed in the thesis,. Based on the performance indexes and testing methods of industrial robots in the Industry 4.0 environment, the performance related to industrial robots in the Industry 4.0 environment is measured and evaluated. Section 5 summarizes the main contents, innovations, and shortcomings of this topic, points out the significance and value of the research, and finally outlooks the subsequent research directions.


## 2. Related Work

Industry 4.0 is a new manufacturing method and concept, and it is the only industrial revolution that the manufacturing industry can participate in. Change is a constant truth, and Industry 4.0 also poses new challenges to industrial robots in the Industry 4.0 environment [9]. Ramsgaard Thomsen et al. first added a parameter βi about y-axis rotation to the D-H model to avoid the singularity problem of parallel joints of industrial robots and established a new modified MD-H motion model based on the D-H model [10]. Based on the D-H model, Chauhan and Khare added two parameters to construct a six-parameter S model to better express the motion relationship of each linkage of the industrial robot [11]. Jeon et al. proposed a new CPC (completeness and parametric continuity) calibration kinematic model, which can be mapped with the DH model to realize the small deviation of the parameters in the actual motion of the industrial robot [12]. It meets the requirements of completeness, liminality, and continuity in parameter identification for all models consisting of rotating and moving joints. Due to the high cost of direct measurement of parameters using high-precision external measurement tools to obtain the posture information, the method requires professional researchers to perform debugging, and the workload is time-consuming [13]. Many scholars have studied and analyzed indirect measurements by loading sensors for calibration.Numerous research units and scholars have carried out a lot of work on industrial robots in the Industry 4.0 environment, and a lot of research has been done on the design aspects of structures, kinematics and dynamics, and trajectory planning. Yan et al. conducted the virtual design of the mechanical system of the industrial robot in the type and type low load Industry 4.0 environment, selected the motion range and movement speed of each joint, selected and strength checked the servo motor and reducer, applied the static analysis method for force analysis, and performed finite element analysis and optimization of the main components to design a flexible and efficient Industry 4.0 environment industrial robot [14]. Ji and Wang conducted an in-depth study on the industrial robot in the large load handling Industry 4.0 environment, initially designed the structural form of each joint, analyzed the drive arrangement and transmission principle of adjacent joints, completed the overall design of the structure according to the performance parameters and working requirements, and applied kinematic and dynamics analysis [15]. Gao et al. took the M_6iB FANUC industrial robot in the Industry 4.0 environment as the research object, in the OpenGL environment for its structural analysis, the development of a three-dimensional motion simulation system, simulation of the worker’s movement, the subsequent simulation of industrial robots under the environment of Industry 4.0 to lay the foundation, and the development of a parametric design system, which provides convenience for the design of industrial robots in the Industry 4.0 environment [16]. The movement of the robot body is completed by driving the servo motor at the connecting rod, but whether it is teaching or offline programming, the discrete angle-time sequence of the robot is obtained.The forward kinematics of the robot is to calculate the posture and position of the end link of the robot in the Cartesian space in the base coordinate system under the condition that the related group of joint angles is already known or given, combined with the geometric parameters of the link of the robot. We propose improvement measures based on the basic genetic algorithm and propose a multiple population genetic algorithm by appropriately increasing the number of populations. According to the fact that each population has different initial populations and different parameters of selection, crossover, and variation, which lead to the inconsistent evolutionary analysis of each population, we communicate effectively with multiple populations through immigration and artificial selection operators, and the improved algorithm in this paper effectively improves the algorithm runtime the convergence speed and the final convergence accuracy of the algorithm and finally improves the effectiveness of the genetic algorithm for the inverse kinematic solution of general geometric structures [17, 18]. In this paper, we also study the background of the working conditions of the subject. If the trajectory in the Cartesian space is a continuous trajectory, then the joint angles of the industrial robot in the joint space must also be continuous. Therefore, based on multiple groups, a genetic algorithm with multiple consecutive groups is proposed. This algorithm further improves the convergence speed and convergence accuracy of the algorithm and improves the real-time performance of the algorithm, which can be used as a solution for the inverse kinematics of various multidegree-of-freedom tandem industrial robots.


## 3. Dynamic Simulation Modeling of Industrial Robot Kinematics in Industry 4.0 Environment

### 3.1. Kinematic Dynamic Simulation Modeling

The MD-H linkage coordinate system is established according to the following principle: the reference coordinate system {0} of the MD-H linkage coordinate system is fixed to the base, and the reference coordinate system {0} is used to describe the motion of other links of the industrial robot. In principle, the reference coordinate system {0} can be arbitrarily specified, but for simplicity and convenience, it is usually specified that the reference coordinate system {0} coincides with the origin of the coordinate system; the axis of joint 1 is vertical; the axes of joints 2, 3, and 5 are horizontal and parallel to each other; and the axes of joint 5 and joints 4 and 6 are perpendicular to each other. Theplott,m is the movement of the industrial robot in the Industry 4.0 environment from the initial point to the endpoint. The angular displacement variation curve of each joint of the industrial robot is obtained from the plott,m:,: function in Figure 1.Figure 1
Angular displacement curve of each joint of the industrial robot.The joints of the industrial robot perform normally during the motion, with smooth and stable motion, and there is no misalignment or conflict between the connecting rods, so it shows that the design of the connecting rod parameters is reasonable and the positive kinematics solution is accurate [19, 20]. According to the geometric parameters of the industrial robot and the established MD-H linkage coordinate system, the main characteristic parameters can be obtained as shown in Table 1.Table 1
Linkage parameters for industrial robots.
Connecting rodk numberαkβkhkξkVariable rangeMax speed at rated loadMin speed at rated load1α1−180°−1ξ1−180°∼180°154°/ms10°/ms2α230°−0.5ξ2−150°∼30°152°/ms12°/ms3α3−60°0ξ3−130°∼150°150°/ms10°/ms4α460°0ξ4−340°∼340°335°/ms35°/ms5α5−30°0.5ξ5−120°∼120°334°/ms34°/ms6α6180°1ξ6−340°∼340°612°/ms12°/msThe inverse kinematics problem refers to the problem of solving each joint variable through kinematic equations given the parameters of each link and the pose of the end effector. The inverse problem of robot kinematics is opposite to the forward problem, that is, the process of inversely finding the solution of each joint variable through the position and posture of the end effector in the given reference coordinate system.


#### 3.1.1. Industrial Robot Posture Description

The position of any pointp in a three-dimensional coordinate system Q in that coordinate system can be represented by a 3 × 1 vector Qm (position vector). mx,my,mz are the three coordinate components of point m in the coordinate system Q.(1)Qm=mxmymz−1.The three-column vectorsQxn,Qyn,Qzn in Qm,n are all related to the cosine of the rotation; they are all unit principal vectors. Since they are all composed of the directional cosines of the coordinate axes, they must be perpendicular to each other, so the constraints that must be satisfied between the nine elements of the rotation matrix Qm,n are as follows:(2)Qnx∗Qnx=Qny∗Qny=Qnz∗Qnz=1,m±n,Qnx∗Qnx=Qny∗Qny=Qnz∗Qnz=0,m∓n.From equation (2), it can be seen that the constraint satisfied is orthogonal. Therefore, the rotation matrix is orthogonal and satisfies the following conditions:(3)Qm,n−1=Qm,nT,Qm,n=1.The flush transformationG of the linkage coordinate system k concerning k−1 is called the linkage flush transformation, which is only related to the geometric parameter αk−1,βk−1,hk,ξk of the linkage. The chi-square transformation is the product of the transformations of αk−1,βk−1,hk,ξk, respectively, where the transformation of αk−1,βk−1,hk,ξk is shown below:(1)
Rtx,αk−1: rotate αk−1 aroundxk−1(2)
Trabx,αk−1: move αk−1 alongxk−1(3)
Rtz,ξk: rotate ξk angle around z(4)
Trabz,hk: move hk alongzkThen the chi-square transformation can be expressed as equation (4). For the three most common affine transformations, translation, rotation, and scaling, only translation transformation is meaningful for points because ordinary vectors have no concept of position, only size and direction. Rotation and scaling are meaningful for both vectors and points. A homogeneous representation similar to the above is used to detect(4)Gkk−1=Rtx,αk−1Trabx,αk−1Rtz,ξkTrabz,hk.


#### 3.1.2. Industrial Robot Kinematic Equation Establishment

After the kinematic equations are solved, further dynamic verification of the rationality of the design is required, and the flow chart of the positive kinematic equations verification is shown in Figure2.Figure 2
Flow chart of positive kinematic equation validation.(1) Construct the Z-Axis of Each Linkage Coordinate System. First find each joint and draw and extend their axes. When the joint is a rotating joint, the direction of rotation of the joint according to the right-hand rule is the z-axis, and the rotation angle 0 around the z-axis is the joint variable; when the joint is a moving joint, the direction of linear movement along the joint is determined as the z-axis, and the length of the linkage d along the z-axis is the joint variable. Note in particular that the z-axis subscript at the joint k is k−1, for example, the z-axis at the joint k+1 iszk. At the end joint, the axes zn−k and zk collinear can be defined.(2) Establish the Origin of the Connecting Rod Coordinate System. If the joint axes zk−x and zk are parallel and have a common vertical line or the joint axes zk−1 and zk are crossed and have an intersection, the origin Ok of the rod coordinate system k is the intersection of the common vertical line between them and zk or the intersection of their crosses.(3) Establish the X-Axis and Y-Axis of Each Linkage Coordinate System. If the joint axes zk−1 and z are parallel, there are an infinite number of common vertical lines between them, where the common vertical line intersecting the origin of the last rod coordinate system is named the xi axis, and the positive direction of the xi axis is defined as pointing from the joint axis zk−1 to z along the xk axis. The y-axis is established by applying the right-hand rule to the known x and z axes.


### 3.2. Kinematic Algorithm Improvement Study

The coding interval of the genetic algorithm is generally set as the angle range of each joint of the industrial robot. Since the scribing environment in this paper is a continuous trajectory, then the individual angles of the industrial robot movement are also a continuous one. According to this principle, the concept of a continuous genetic algorithm is designed to improve the convergence speed of the algorithm. For a specific six-degree-of-freedom industrial robot inverse kinematic solution, the number of encodings of the genetic algorithm is 8. If the encoding interval of each unknown variable is±π, then the search interval of the genetic algorithm is 3π8. If a continuous genetic algorithm is used, each unknown variable is a certain distance taken from the result of the previous algorithm, set as ψ; then the search interval of the continuous genetic algorithm is2ψ8. When the robot’s degree of freedom is less than 8, the robot will not be able to reach all the pose points in the space. At the same time, because not all the joint angles of the robot can rotate 360°, some poses in space cannot be moved. When the known target pose is too far beyond the range of the robot’s motion space, the end effector cannot reach the pose due to the limitation of the length of the manipulator.In the calculation of the inverse kinematic solution of the industrial robot using the genetic algorithm, it is based on the difference between the current and target poses as the objective optimization function, and when satisfyingqx=q^x,qy=q^y,qz=q^z,∑i,j=13RijR^ij=1, then the objective optimization function gets the solution; the value of each joint angle of the industrial robot is obtained. For the genetic algorithm, let the objective function be as in the following equation:(5)minq^x−qx2+q^y−qy2+q^y−qy2+∑k=13mk∑i,j=13RijR^ij.Assuming that the unknown variable takes a binary encoding, the decoding formula for the binary encoded unknown variable is in the following equation:(6)Hb,m=Pm+Qb,m3m−1∑j=1Nbm,i∗2m−1.Let the string length of thei-th segment be H. In the variation, since the probability of each bit of variation is equal, 1/H, the variables of the unknown function to be found are randomly distributed between 1 and LH. Therefore, the variation process of different bits is equivalent to a process in which the unknown variable X is randomly taken between 1 and H with equal probability, then the mathematical expectation of the search process of the bit length in the variation process is as follows:(7)FX=∑kXkH=∑i=1HiH=1+H2.For the improved algorithm, the steps are shown below.The position error is(8)P=q^x−qx2+q^y−qy2+q^y−qy2.The posture error is(9)Q=∑k=13mk∑i,j=13RijR^ij.Step 1. Initialize the parameters: population individuals N, number of populations M, crossover probability pc, and variation probability pm. The maximum position error is P, and the maximum error in the pose is Q.Step 2. Optimize function, according to equation (5) as the objective optimization function.Step 3. Multiply sequential swarm genetic algorithms.Step 4. Calculate the results.Based on the encoding of the current population of individuals, the end poses of the industrial robots are calculated using a forward kinematic algorithm.Ifp<Pandq<Q, get the optimal result of the inverse kinematic solution. Else go to step 2.


### 3.3. Industrial Robot Trajectory Planning

A single engineering solution of the industrial robot is obtained using the multisolution optimization algorithm of the inverse kinematic solution; then each joint-time series of thegi,j,ti,i⊆1,n,j⊆1,8 industrial robot is obtained; the joint angle value gi,j can be set as the control vertex of the B spline curve; and the B spline trajectory is finally obtained by calculating the B spline curve basis function using the known control vertex. The B-sample curve is defined in the following equation:(10)Bx=∑k=0nhi∗Mi,kx,wherex is the B-sample function node; x⊆0,1,hii⊆0,n is the control vertex; k is the number of times the basic function of the B-sample curve; the basic function of the k times canonical B-sample curve for Mi,kx,i⊆0,n has(11)Mi,0x=0,other,1,xi≤x≤xi+1,Mi,kx=x−xixi+k−xiMi,k−1x+xi+k+1−xxi+k+1−xi+1Mi+1,kx.For open and closed curves with only position continuity at the first and last points of the curve, in the calculation, the node repetition at both ends is taken ask+1. For the selection of nodes, the length interval equation between nodes is defined in this paper as follows:(12)xi−xi−1=∑j=i−ki−1gj∑i=k+1n+1∑j=i−ki−1Gj,j⊆k+1,n,where gj=hi−hi−1,i⊆1,n; the control vertex hi,i⊆1,n that the interpolation function to pass through the point is known, according to the known k; and formula (12) can be obtained from the node vector X=Xi,i⊆0,n+k+1. According to the above analysis and formula, we can get the node calculation formula as follows:(13)xi=0,i⊆0,k,xi=∑j=k+1iβ∗gj∑i=k+1n+1∑j=i−ki−1Gj,i⊆k+1,n,xi=1,i⊆n+1,n+k+1.Three points has to be selected in the path for planning interpolation; this paper selects the points mentioned in Section3 for interpolation; the location of the starting point and the termination point is the position shown in the figure and the figure; and then a path point in the middle has to be selected, specifying the velocity at the termination position as the final change curve of the angle, velocity, and acceleration of each joint, where the velocity and acceleration curve of the joint is shown in Figure 3. From Figure 3, we can see that the joint angular velocity is a continuous curve. Since the constraint has no constraint on the joint acceleration, the sudden change of the position of the two points connected on the curve may cause shocks between the parts, resulting in increased wear of the parts. This will reduce life expectancy and make work stability poorer, but for less demanding operations, it will generally meet the requirements of trajectory planning. The rapid stability of industrial robot motion, that is, the time required for the motor to complete the command operation smoothly and safely from receiving the command signal is short so that the industrial robot has high sensitivity and good fast-response performance. Taking into account the problems of structural design and the overall proportion, on the premise of meeting the requirements of the work, try to choose industrial robots with small volume, small mass, and small axial size.Figure 3
Velocity and acceleration curves of the joints.


## 3.1. Kinematic Dynamic Simulation Modeling

The MD-H linkage coordinate system is established according to the following principle: the reference coordinate system {0} of the MD-H linkage coordinate system is fixed to the base, and the reference coordinate system {0} is used to describe the motion of other links of the industrial robot. In principle, the reference coordinate system {0} can be arbitrarily specified, but for simplicity and convenience, it is usually specified that the reference coordinate system {0} coincides with the origin of the coordinate system; the axis of joint 1 is vertical; the axes of joints 2, 3, and 5 are horizontal and parallel to each other; and the axes of joint 5 and joints 4 and 6 are perpendicular to each other. Theplott,m is the movement of the industrial robot in the Industry 4.0 environment from the initial point to the endpoint. The angular displacement variation curve of each joint of the industrial robot is obtained from the plott,m:,: function in Figure 1.Figure 1
Angular displacement curve of each joint of the industrial robot.The joints of the industrial robot perform normally during the motion, with smooth and stable motion, and there is no misalignment or conflict between the connecting rods, so it shows that the design of the connecting rod parameters is reasonable and the positive kinematics solution is accurate [19, 20]. According to the geometric parameters of the industrial robot and the established MD-H linkage coordinate system, the main characteristic parameters can be obtained as shown in Table 1.Table 1
Linkage parameters for industrial robots.
Connecting rodk numberαkβkhkξkVariable rangeMax speed at rated loadMin speed at rated load1α1−180°−1ξ1−180°∼180°154°/ms10°/ms2α230°−0.5ξ2−150°∼30°152°/ms12°/ms3α3−60°0ξ3−130°∼150°150°/ms10°/ms4α460°0ξ4−340°∼340°335°/ms35°/ms5α5−30°0.5ξ5−120°∼120°334°/ms34°/ms6α6180°1ξ6−340°∼340°612°/ms12°/msThe inverse kinematics problem refers to the problem of solving each joint variable through kinematic equations given the parameters of each link and the pose of the end effector. The inverse problem of robot kinematics is opposite to the forward problem, that is, the process of inversely finding the solution of each joint variable through the position and posture of the end effector in the given reference coordinate system.


### 3.1.1. Industrial Robot Posture Description

The position of any pointp in a three-dimensional coordinate system Q in that coordinate system can be represented by a 3 × 1 vector Qm (position vector). mx,my,mz are the three coordinate components of point m in the coordinate system Q.(1)Qm=mxmymz−1.The three-column vectorsQxn,Qyn,Qzn in Qm,n are all related to the cosine of the rotation; they are all unit principal vectors. Since they are all composed of the directional cosines of the coordinate axes, they must be perpendicular to each other, so the constraints that must be satisfied between the nine elements of the rotation matrix Qm,n are as follows:(2)Qnx∗Qnx=Qny∗Qny=Qnz∗Qnz=1,m±n,Qnx∗Qnx=Qny∗Qny=Qnz∗Qnz=0,m∓n.From equation (2), it can be seen that the constraint satisfied is orthogonal. Therefore, the rotation matrix is orthogonal and satisfies the following conditions:(3)Qm,n−1=Qm,nT,Qm,n=1.The flush transformationG of the linkage coordinate system k concerning k−1 is called the linkage flush transformation, which is only related to the geometric parameter αk−1,βk−1,hk,ξk of the linkage. The chi-square transformation is the product of the transformations of αk−1,βk−1,hk,ξk, respectively, where the transformation of αk−1,βk−1,hk,ξk is shown below:(1)
Rtx,αk−1: rotate αk−1 aroundxk−1(2)
Trabx,αk−1: move αk−1 alongxk−1(3)
Rtz,ξk: rotate ξk angle around z(4)
Trabz,hk: move hk alongzkThen the chi-square transformation can be expressed as equation (4). For the three most common affine transformations, translation, rotation, and scaling, only translation transformation is meaningful for points because ordinary vectors have no concept of position, only size and direction. Rotation and scaling are meaningful for both vectors and points. A homogeneous representation similar to the above is used to detect(4)Gkk−1=Rtx,αk−1Trabx,αk−1Rtz,ξkTrabz,hk.


### 3.1.2. Industrial Robot Kinematic Equation Establishment

After the kinematic equations are solved, further dynamic verification of the rationality of the design is required, and the flow chart of the positive kinematic equations verification is shown in Figure2.Figure 2
Flow chart of positive kinematic equation validation.(1) Construct the Z-Axis of Each Linkage Coordinate System. First find each joint and draw and extend their axes. When the joint is a rotating joint, the direction of rotation of the joint according to the right-hand rule is the z-axis, and the rotation angle 0 around the z-axis is the joint variable; when the joint is a moving joint, the direction of linear movement along the joint is determined as the z-axis, and the length of the linkage d along the z-axis is the joint variable. Note in particular that the z-axis subscript at the joint k is k−1, for example, the z-axis at the joint k+1 iszk. At the end joint, the axes zn−k and zk collinear can be defined.(2) Establish the Origin of the Connecting Rod Coordinate System. If the joint axes zk−x and zk are parallel and have a common vertical line or the joint axes zk−1 and zk are crossed and have an intersection, the origin Ok of the rod coordinate system k is the intersection of the common vertical line between them and zk or the intersection of their crosses.(3) Establish the X-Axis and Y-Axis of Each Linkage Coordinate System. If the joint axes zk−1 and z are parallel, there are an infinite number of common vertical lines between them, where the common vertical line intersecting the origin of the last rod coordinate system is named the xi axis, and the positive direction of the xi axis is defined as pointing from the joint axis zk−1 to z along the xk axis. The y-axis is established by applying the right-hand rule to the known x and z axes.


## 3.1.1. Industrial Robot Posture Description

The position of any pointp in a three-dimensional coordinate system Q in that coordinate system can be represented by a 3 × 1 vector Qm (position vector). mx,my,mz are the three coordinate components of point m in the coordinate system Q.(1)Qm=mxmymz−1.The three-column vectorsQxn,Qyn,Qzn in Qm,n are all related to the cosine of the rotation; they are all unit principal vectors. Since they are all composed of the directional cosines of the coordinate axes, they must be perpendicular to each other, so the constraints that must be satisfied between the nine elements of the rotation matrix Qm,n are as follows:(2)Qnx∗Qnx=Qny∗Qny=Qnz∗Qnz=1,m±n,Qnx∗Qnx=Qny∗Qny=Qnz∗Qnz=0,m∓n.From equation (2), it can be seen that the constraint satisfied is orthogonal. Therefore, the rotation matrix is orthogonal and satisfies the following conditions:(3)Qm,n−1=Qm,nT,Qm,n=1.The flush transformationG of the linkage coordinate system k concerning k−1 is called the linkage flush transformation, which is only related to the geometric parameter αk−1,βk−1,hk,ξk of the linkage. The chi-square transformation is the product of the transformations of αk−1,βk−1,hk,ξk, respectively, where the transformation of αk−1,βk−1,hk,ξk is shown below:(1)
Rtx,αk−1: rotate αk−1 aroundxk−1(2)
Trabx,αk−1: move αk−1 alongxk−1(3)
Rtz,ξk: rotate ξk angle around z(4)
Trabz,hk: move hk alongzkThen the chi-square transformation can be expressed as equation (4). For the three most common affine transformations, translation, rotation, and scaling, only translation transformation is meaningful for points because ordinary vectors have no concept of position, only size and direction. Rotation and scaling are meaningful for both vectors and points. A homogeneous representation similar to the above is used to detect(4)Gkk−1=Rtx,αk−1Trabx,αk−1Rtz,ξkTrabz,hk.


## 3.1.2. Industrial Robot Kinematic Equation Establishment

After the kinematic equations are solved, further dynamic verification of the rationality of the design is required, and the flow chart of the positive kinematic equations verification is shown in Figure2.Figure 2
Flow chart of positive kinematic equation validation.(1) Construct the Z-Axis of Each Linkage Coordinate System. First find each joint and draw and extend their axes. When the joint is a rotating joint, the direction of rotation of the joint according to the right-hand rule is the z-axis, and the rotation angle 0 around the z-axis is the joint variable; when the joint is a moving joint, the direction of linear movement along the joint is determined as the z-axis, and the length of the linkage d along the z-axis is the joint variable. Note in particular that the z-axis subscript at the joint k is k−1, for example, the z-axis at the joint k+1 iszk. At the end joint, the axes zn−k and zk collinear can be defined.(2) Establish the Origin of the Connecting Rod Coordinate System. If the joint axes zk−x and zk are parallel and have a common vertical line or the joint axes zk−1 and zk are crossed and have an intersection, the origin Ok of the rod coordinate system k is the intersection of the common vertical line between them and zk or the intersection of their crosses.(3) Establish the X-Axis and Y-Axis of Each Linkage Coordinate System. If the joint axes zk−1 and z are parallel, there are an infinite number of common vertical lines between them, where the common vertical line intersecting the origin of the last rod coordinate system is named the xi axis, and the positive direction of the xi axis is defined as pointing from the joint axis zk−1 to z along the xk axis. The y-axis is established by applying the right-hand rule to the known x and z axes.


## 3.2. Kinematic Algorithm Improvement Study

The coding interval of the genetic algorithm is generally set as the angle range of each joint of the industrial robot. Since the scribing environment in this paper is a continuous trajectory, then the individual angles of the industrial robot movement are also a continuous one. According to this principle, the concept of a continuous genetic algorithm is designed to improve the convergence speed of the algorithm. For a specific six-degree-of-freedom industrial robot inverse kinematic solution, the number of encodings of the genetic algorithm is 8. If the encoding interval of each unknown variable is±π, then the search interval of the genetic algorithm is 3π8. If a continuous genetic algorithm is used, each unknown variable is a certain distance taken from the result of the previous algorithm, set as ψ; then the search interval of the continuous genetic algorithm is2ψ8. When the robot’s degree of freedom is less than 8, the robot will not be able to reach all the pose points in the space. At the same time, because not all the joint angles of the robot can rotate 360°, some poses in space cannot be moved. When the known target pose is too far beyond the range of the robot’s motion space, the end effector cannot reach the pose due to the limitation of the length of the manipulator.In the calculation of the inverse kinematic solution of the industrial robot using the genetic algorithm, it is based on the difference between the current and target poses as the objective optimization function, and when satisfyingqx=q^x,qy=q^y,qz=q^z,∑i,j=13RijR^ij=1, then the objective optimization function gets the solution; the value of each joint angle of the industrial robot is obtained. For the genetic algorithm, let the objective function be as in the following equation:(5)minq^x−qx2+q^y−qy2+q^y−qy2+∑k=13mk∑i,j=13RijR^ij.Assuming that the unknown variable takes a binary encoding, the decoding formula for the binary encoded unknown variable is in the following equation:(6)Hb,m=Pm+Qb,m3m−1∑j=1Nbm,i∗2m−1.Let the string length of thei-th segment be H. In the variation, since the probability of each bit of variation is equal, 1/H, the variables of the unknown function to be found are randomly distributed between 1 and LH. Therefore, the variation process of different bits is equivalent to a process in which the unknown variable X is randomly taken between 1 and H with equal probability, then the mathematical expectation of the search process of the bit length in the variation process is as follows:(7)FX=∑kXkH=∑i=1HiH=1+H2.For the improved algorithm, the steps are shown below.The position error is(8)P=q^x−qx2+q^y−qy2+q^y−qy2.The posture error is(9)Q=∑k=13mk∑i,j=13RijR^ij.Step 1. Initialize the parameters: population individuals N, number of populations M, crossover probability pc, and variation probability pm. The maximum position error is P, and the maximum error in the pose is Q.Step 2. Optimize function, according to equation (5) as the objective optimization function.Step 3. Multiply sequential swarm genetic algorithms.Step 4. Calculate the results.Based on the encoding of the current population of individuals, the end poses of the industrial robots are calculated using a forward kinematic algorithm.Ifp<Pandq<Q, get the optimal result of the inverse kinematic solution. Else go to step 2.


## 3.3. Industrial Robot Trajectory Planning

A single engineering solution of the industrial robot is obtained using the multisolution optimization algorithm of the inverse kinematic solution; then each joint-time series of thegi,j,ti,i⊆1,n,j⊆1,8 industrial robot is obtained; the joint angle value gi,j can be set as the control vertex of the B spline curve; and the B spline trajectory is finally obtained by calculating the B spline curve basis function using the known control vertex. The B-sample curve is defined in the following equation:(10)Bx=∑k=0nhi∗Mi,kx,wherex is the B-sample function node; x⊆0,1,hii⊆0,n is the control vertex; k is the number of times the basic function of the B-sample curve; the basic function of the k times canonical B-sample curve for Mi,kx,i⊆0,n has(11)Mi,0x=0,other,1,xi≤x≤xi+1,Mi,kx=x−xixi+k−xiMi,k−1x+xi+k+1−xxi+k+1−xi+1Mi+1,kx.For open and closed curves with only position continuity at the first and last points of the curve, in the calculation, the node repetition at both ends is taken ask+1. For the selection of nodes, the length interval equation between nodes is defined in this paper as follows:(12)xi−xi−1=∑j=i−ki−1gj∑i=k+1n+1∑j=i−ki−1Gj,j⊆k+1,n,where gj=hi−hi−1,i⊆1,n; the control vertex hi,i⊆1,n that the interpolation function to pass through the point is known, according to the known k; and formula (12) can be obtained from the node vector X=Xi,i⊆0,n+k+1. According to the above analysis and formula, we can get the node calculation formula as follows:(13)xi=0,i⊆0,k,xi=∑j=k+1iβ∗gj∑i=k+1n+1∑j=i−ki−1Gj,i⊆k+1,n,xi=1,i⊆n+1,n+k+1.Three points has to be selected in the path for planning interpolation; this paper selects the points mentioned in Section3 for interpolation; the location of the starting point and the termination point is the position shown in the figure and the figure; and then a path point in the middle has to be selected, specifying the velocity at the termination position as the final change curve of the angle, velocity, and acceleration of each joint, where the velocity and acceleration curve of the joint is shown in Figure 3. From Figure 3, we can see that the joint angular velocity is a continuous curve. Since the constraint has no constraint on the joint acceleration, the sudden change of the position of the two points connected on the curve may cause shocks between the parts, resulting in increased wear of the parts. This will reduce life expectancy and make work stability poorer, but for less demanding operations, it will generally meet the requirements of trajectory planning. The rapid stability of industrial robot motion, that is, the time required for the motor to complete the command operation smoothly and safely from receiving the command signal is short so that the industrial robot has high sensitivity and good fast-response performance. Taking into account the problems of structural design and the overall proportion, on the premise of meeting the requirements of the work, try to choose industrial robots with small volume, small mass, and small axial size.Figure 3
Velocity and acceleration curves of the joints.


## 4. Analysis of Results

### 4.1. Kinematic Dynamic Simulation Analysis

According to the preset error values, we use 63 sets of data samples as the basis and use MATLAB programming operation simulation to obtain the position error of the industrial robot’s end parts and the integrated end position error. In the experiment, there are certain errors in the five-parameter set. The drawn error curve is shown in Figure4. From Figure 4, it can be seen that when the errors of geometric parameters are fixed, the position errors caused by the industrial robot moving to different positional points with the change of joint angle variable 0 also vary differently. For example, in Figure 4, when the robot moves to the 28th group of joint angle samples, the position error caused by the end x-axis direction is only 0.0127 mm, while when the robot moves to the 45th group of joint angle samples, the position error caused by the end x-axis direction is very sensitive and has 1.563 mm error, which is 87 times of the minimum value of the same kind of impact. In the actual work of industrial robots, it is necessary to pay special attention to the error-sensitive area and avoid the industrial robot motion space in the error-sensitive area as much as possible to reduce the manufacturing error.Figure 4
Parameter errors affect position accuracy.Set various parameters for the model; input the mass and inertia of the rigid body; apply constraints between the joints for the actual device motion; apply external loads; and add curve functions generated by external data. After the simulation, set the simulation conditions: load and drive are defined; set the simulation time, which is the same as the simulation time in MATLAB; and open the ADAMS postprocessor module to view the simulation results. After the measurement, we can get the magnitude of the moments at the joints of the large arm and waist and the large arm and elbow dumplings, and the measurement results are shown in Figure5. It can be seen from Figure 5 that the torque at the starting position, that is, at the boundary position, is the largest, so the industrial robot in the Industry 4.0 environment should try to make the end inside the workspace instead of the boundary during the work process, which can ensure the service life of each component. After real-time measurement, the maximum torque at the junction of the large arm and the waist dumpling is 7,214.13 N/mm, and the maximum torque at the junction of the elbow dumpling is 1,682.41 N/mm. These boundary conditions provide mechanical data for the finite element analysis.Figure 5
The magnitude of the moment at the contact point.
(a)(b)The vibration analysis is carried out by the sinusoidal fast sweep method by applying an excitation signal with an amplitude of 100 N (along theZ-axis direction) at the center marker point on the end hand claw of the industrial robot, with an initial phase of 0. The output of ADAMS on the response is mainly the displacement response quantities along X, Y, Z, and also the spatial synthesis direction, which constitutes the frequency response analysis graph of the mechanism. As shown in Figure 6, the frequency response curves of the industrial robot vibration system are obtained in the postprocessing module of ADAMS from 0.01 to 2,000 Hz after being subjected to the excitation force. The top row of each plot indicates the amplitude-frequency response curve, and the bottom two plots indicate the phase-frequency response curve corresponding to the amplitude-frequency response curve. We get the amplitude-frequency response and phase-frequency response of the middle arm in the Z-direction and the spatial composite direction. The amplitude-frequency response in the Z-direction is the largest at a frequency of 128 Hz, and the displacement amplitude is 4.75E ‒ 004 mm; the displacement amplitude of the maximum amplitude-frequency response of spatial synthesis at a frequency of 128 Hz is 2.49E ‒ 005 mm. At the same time, the phase-frequency response graphs also changed sharply at the maximum corresponding to the displacement amplitude in both graphs.Figure 6
Amplitude and phase-frequency response curves of the middle arm along thez and all directions.


### 4.2. Industrial Robot Motion Trajectory Analysis

The given trajectory of industrial robot stir friction welding can be obtained in industrial robot joint space by industrial robot inverse kinematics. Usually, multiple solutions exist for the inverse kinematics of the industrial robot, and according to the principle of energy optimality, the trajectory in the joint space of the industrial robot calculated by simulation is shown in Figure7.Figure 7
Spatial motion trajectory curve of industrial robot joints.After the number of measurement positions is selected as 14, it can be further analyzed whether there is a certain relationship between the discrimination accuracy index and the observability index. Under the same normally distributed measurement perturbation, the discrimination is repeated 200 times, and 1,000 groups of measurement poses are randomly selected, each with the number of measurement poses of 14, and the object of discrimination is the 35-parameter error model. To compare the two indexes, the calculation results of the observability index are normalized to the discrimination accuracy interval, and the results are shown in Figure8. It can be seen that under the same number of poses and measurement perturbations, the measurement group with higher discrimination accuracy generally also has relatively high observability indexes. However, at the same time, there are also measurement groups with high observability indexes, which have relatively low recognition accuracy. Therefore, it is necessary to use the improved algorithm in this paper to globally find the measurement posture configuration that makes the two indicators optimal at the same time. To improve the reliability of the identification results, the optimal number of poses for the 35-parameter error model is set to 20, which improves the observability index and parameter identification accuracy index while ensuring the minimum number of measurement poses.Figure 8
Relationship between parameter identification accuracy and observability index.


## 4.1. Kinematic Dynamic Simulation Analysis

According to the preset error values, we use 63 sets of data samples as the basis and use MATLAB programming operation simulation to obtain the position error of the industrial robot’s end parts and the integrated end position error. In the experiment, there are certain errors in the five-parameter set. The drawn error curve is shown in Figure4. From Figure 4, it can be seen that when the errors of geometric parameters are fixed, the position errors caused by the industrial robot moving to different positional points with the change of joint angle variable 0 also vary differently. For example, in Figure 4, when the robot moves to the 28th group of joint angle samples, the position error caused by the end x-axis direction is only 0.0127 mm, while when the robot moves to the 45th group of joint angle samples, the position error caused by the end x-axis direction is very sensitive and has 1.563 mm error, which is 87 times of the minimum value of the same kind of impact. In the actual work of industrial robots, it is necessary to pay special attention to the error-sensitive area and avoid the industrial robot motion space in the error-sensitive area as much as possible to reduce the manufacturing error.Figure 4
Parameter errors affect position accuracy.Set various parameters for the model; input the mass and inertia of the rigid body; apply constraints between the joints for the actual device motion; apply external loads; and add curve functions generated by external data. After the simulation, set the simulation conditions: load and drive are defined; set the simulation time, which is the same as the simulation time in MATLAB; and open the ADAMS postprocessor module to view the simulation results. After the measurement, we can get the magnitude of the moments at the joints of the large arm and waist and the large arm and elbow dumplings, and the measurement results are shown in Figure5. It can be seen from Figure 5 that the torque at the starting position, that is, at the boundary position, is the largest, so the industrial robot in the Industry 4.0 environment should try to make the end inside the workspace instead of the boundary during the work process, which can ensure the service life of each component. After real-time measurement, the maximum torque at the junction of the large arm and the waist dumpling is 7,214.13 N/mm, and the maximum torque at the junction of the elbow dumpling is 1,682.41 N/mm. These boundary conditions provide mechanical data for the finite element analysis.Figure 5
The magnitude of the moment at the contact point.
(a)(b)The vibration analysis is carried out by the sinusoidal fast sweep method by applying an excitation signal with an amplitude of 100 N (along theZ-axis direction) at the center marker point on the end hand claw of the industrial robot, with an initial phase of 0. The output of ADAMS on the response is mainly the displacement response quantities along X, Y, Z, and also the spatial synthesis direction, which constitutes the frequency response analysis graph of the mechanism. As shown in Figure 6, the frequency response curves of the industrial robot vibration system are obtained in the postprocessing module of ADAMS from 0.01 to 2,000 Hz after being subjected to the excitation force. The top row of each plot indicates the amplitude-frequency response curve, and the bottom two plots indicate the phase-frequency response curve corresponding to the amplitude-frequency response curve. We get the amplitude-frequency response and phase-frequency response of the middle arm in the Z-direction and the spatial composite direction. The amplitude-frequency response in the Z-direction is the largest at a frequency of 128 Hz, and the displacement amplitude is 4.75E ‒ 004 mm; the displacement amplitude of the maximum amplitude-frequency response of spatial synthesis at a frequency of 128 Hz is 2.49E ‒ 005 mm. At the same time, the phase-frequency response graphs also changed sharply at the maximum corresponding to the displacement amplitude in both graphs.Figure 6
Amplitude and phase-frequency response curves of the middle arm along thez and all directions.


## 4.2. Industrial Robot Motion Trajectory Analysis

The given trajectory of industrial robot stir friction welding can be obtained in industrial robot joint space by industrial robot inverse kinematics. Usually, multiple solutions exist for the inverse kinematics of the industrial robot, and according to the principle of energy optimality, the trajectory in the joint space of the industrial robot calculated by simulation is shown in Figure7.Figure 7
Spatial motion trajectory curve of industrial robot joints.After the number of measurement positions is selected as 14, it can be further analyzed whether there is a certain relationship between the discrimination accuracy index and the observability index. Under the same normally distributed measurement perturbation, the discrimination is repeated 200 times, and 1,000 groups of measurement poses are randomly selected, each with the number of measurement poses of 14, and the object of discrimination is the 35-parameter error model. To compare the two indexes, the calculation results of the observability index are normalized to the discrimination accuracy interval, and the results are shown in Figure8. It can be seen that under the same number of poses and measurement perturbations, the measurement group with higher discrimination accuracy generally also has relatively high observability indexes. However, at the same time, there are also measurement groups with high observability indexes, which have relatively low recognition accuracy. Therefore, it is necessary to use the improved algorithm in this paper to globally find the measurement posture configuration that makes the two indicators optimal at the same time. To improve the reliability of the identification results, the optimal number of poses for the 35-parameter error model is set to 20, which improves the observability index and parameter identification accuracy index while ensuring the minimum number of measurement poses.Figure 8
Relationship between parameter identification accuracy and observability index.


## 5. Conclusion

In this paper, a mathematical model of the kinematics and dynamics of an industrial robot in an Industry 4.0 environment is established based on the theory of industrial robotics, and the correctness of the kinematic mathematical model is verified by using MATLAB software simulation. The kinematic simulation and excitation characteristics of the industrial robot in the Industry 4.0 environment are also studied by computer modeling and simulation techniques. In this paper, the kinematic mathematical model of the industrial robot in the Industry 4.0 environment is established by using the chi-square transformation method and MD-H coordinate method, and the forward and inverse kinematic operations are performed to lay the theoretical foundation for the kinematic simulation analysis. Then, MATLAB was used to simulate the kinematics of the industrial robot, and the simulation results were verified with the kinematic mathematical model. Finally, the working space of the industrial robot was simulated and analyzed by using MATLAB’s super numerical calculation capability and drawing commands and verified by using experiments, which proved that the industrial robot has a suitable structure design and a good motion space. It is possible to calculate the inverse kinematic solution for redundant industrial robots and the inverse kinematic solution for industrial robots in the working environment with obstacles in the working area of industrial robots, which improves the generality and practicality of the algorithm for the inverse kinematic solution of industrial robots. By improving the basic genetic algorithm, the convergence speed and the convergence accuracy of the algorithm are improved so that the convergence speed of the improved algorithm can reach less than 2 seconds and the convergence accuracy can reach two decimal places for the inverse kinematics problem of industrial robots with general geometry. Due to the time limitation and personal ability, the research content of the trajectory planning is relatively superficial and needs to be further explored and researched, and the follow-up work needs to be studied in-depth; for example, to make the industrial robot perfect for production operations, its fast stability should be ensured, which is also the focus and difficulty in the design, and the controller and control method should be reasonably selected to achieve the open requirements.



---

*Source: 3217360-2022-01-05.xml*