11384 FACTA UNIVERSITATIS Series: Mechanical Engineering https://doi.org/10.22190/FUME230215017Z © 2020 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper MODELING AND SIMULATION OF INDUSTRIAL ROBOT ARMS USING SIMSCAPE MULTIBODY Cheng Zhang1, 2, Songxiao Li3, Zhuo Zhang4 1Faculty of Innovation Engineering, Macau University of Science and Technology, China 2School of Electronic and Information Engineering, Zhuhai College of Science and Technology, China 3Department of Mathematics, Shantou University, China 4College of Robotics, Gaungdong Polytechnic of Science and Technology, China Abstract. The dynamic simulation modeling problem of industrial robot arm is solved, and the trajectory planning dynamic simulation is performed in this paper. In response to the lack of trajectory planning and motion controller interfaces in the robotic modelling study, including the lack of dynamic simulation visualization, a Simscape Multibody-based method for building a dynamic model of industrial robot arm is proposed and the effectiveness of the model is verified through dynamic simulation. The simulation model integrates the robotic arm trajectory planning, motion controller and data acquisition module. It has a clear structure and the parameters are easy to modify. It can reasonably simulate the structure and parameters of the research object and facilitate the subsequent research of related algorithms. It provides an innovative and open-source research and development platform for the dynamic simulation study of the robot arm. Keywords: Robotic Dynamic Model, Industrial Robot Arms, Robot Simulation System, Simscape Multibody 1. INTRODUCTION Robots are known as the "pearl at the crown of manufacturing". At present, the world is experiencing the wave of "machine replacement", and artificial intelligence has become a strategic technology leading a new round of scientific and technological revolution and industrial transformation. As a representative of artificial intelligence, industrial robots integrate advanced technologies in the fields of machinery, electronics, sensors, wireless Received: February 15, 2023 / Accepted May 28, 2023 Corresponding author: Songxiao Li Department of Mathematics, Shantou University, 515063, Shantou, Guangdong, China E-mail: jyulsx@163.com 2 C. ZHANG, S. LI, Z. ZHANG communication, sound recognition, image processing and artificial intelligence, involving multiple disciplines. It is a country's level of science and technology development and an important symbol of national economy modernization and informatization. At present, various industries have integrated industrial robots into the production line, and made use of industrial robot technology to meet the actual needs of the development of circular economy and improve the development of time efficiency, making the unique advantages of higher production efficiency and lower cost. Compared with traditional industrial processing auxiliary equipment, industrial robots have lower maintenance costs, higher degree of intelligence and better flexibility to meet the requirements of modern production. It can better adapt to a variety of varieties, small batch, on-site processing requirements. On the whole, industrial robots have become a wide range of tools for various tasks in modern society, forming a new generation of "intelligent robot manufacturing"[1]. Under these circumstances, it is particularly important to further research on industrial robot related technology. Many scholars have carried out research and achieved many achievements, including but not limited to the problems of industrial robot grasping, assembly, process control and industrial man-machine cooperation in typical operations. A group of authors [2-4] studied the trajectory planning strategy of industrial robot arm and optimized it. Özyer [5], Yilmaz et al. [6] and Gharib [7] carried out studies based on algorithms such as adaptive neural network control, adaptive fuzzy control and robust control for industrial robot arm. Rojas-Garcia et al. [8] and Gutierrez-Giles et al. [9] also proposed optimization methods for industrial robot arm control strategy from the perspective of force control. The premise of the various robot algorithms is how to build a precise, flexible model. The effectiveness and accuracy of the dynamic model of the industrial robot arm are crucial for its practical application. Usually, robot manufacturers only provide kinematic data, but this can only guarantee a certain level of localizer accuracy, without than insight into its dynamic features. Therefore, how to establish the modeling of industrial robot arm is a very important research content in industrial robot technology. It is the primary premise of all follow-up research, and also one of the hot spots and key points in the field of robot research. In recent years, scholars have studied the modeling of industrial robot arm. In general, the common methods for dynamic modeling of industrial robot arm include Newton-Euler method, Lagrange method, robot modeling with MATLAB / SIMIULINK or robot toolbox, robot modeling with multi-body system, and parameter identification method based on acceleration, speed, position and voltage, etc. Urrea and Pascal [10] introduced the acquisition process of the dynamic equation of the 6R articulated industrial robot, and obtains the expression of the generalized dynamic force of the robot to represent the equation of the inverse dynamic model. Fernández et al. [11] used the Euler-Lagrange equation to model the robot. Bugnar et al. [12] studied the transfer moments of the industrial robot arm joints, developed a mathematical model of the 6-DOF robot, and determined the kinematic and kinetic equations. Urrea and Pascal [13] derived a direct kinematic mathematical model, an inverse kinematic equation, and a dynamic model of the effector torque in the joint. Šulekić et al. [14] studied the forward and reverse kinematics of the DOBOT kinematic robot and simulated it in the Matlab environment. Ho et al. [15] proposed a fast inverse kinematics algorithm, which has a closed solution for a specific 6-DOF arm. Bien [16] proposed the solution algorithm for a closed inverse kinematic model of a human two-arm robot, and verified the effectiveness of the algorithm through Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 3 simulation. Wang et al. [17, 18] proposed dynamic model for articulated industrial robots suitable for machining applications, and determined the link and rotor inertia, joint stiffness and damping parameters. Alam et al. [19] and Nguyen et al. [20] respectively designed the parameter identification model of industrial robots, which is formulated based on acceleration, speed, position and voltage. Harfoush et al. [21] proposed an algorithm based on the position of a set of measured end effectors to identify the sources of error included in the kinematic model. Huynh et al. [22] describes the data-driven modeling approach and its application in controlling a new set of series elastic actuators. Segota et al. [23] used the multilayer perceptron (MLP) algorithm to model the robot to obtain the best combination of output parameters. Truc et al. [24] defined an implementation method for solving the optimal geometric structure model of the plane z-type mechanism, and verified the method by using the geometric modeling software. Barrientos-Diez et al. [25] presented a fast and low-cost fencing robotic arm modeling system. Ogbemhe et al. [26] provides a workspace determination and analysis based on spatial and plane parallel manipulator graphics technology, and completes its implementation in the CATIA workbench using CAD software. Aboulissane et al. [27] introduced a kinematic model of robot in drive and configuration space. After the comparison and analysis of the above modeling methods, the following remarks can be made: (1) Some modeling methods require special measuring of the physical parameters such as the center of mass and mass of the robot. Especially when the model parameters change, the measurement is cumbersome and the accuracy is not always sufficiently good. (2) Some studied robot models are not intuitive and convenient for subsequent trajectory planning and control strategy data. (3) Some modeling methods only discuss the simple axial range of axis motion, not including the case of physical interference between joints, while the coupling between joints is a relatively common case. The main contribution of this paper is to establish a dynamic simulation model of an industrial robot arm, and to simulate the motion trajectory planning of the robot involving all the above mentioned challenging factors. The parameters of the center and mass of the model are automatically calculated by the multi-body system, which is accurate and convenient, and can satisfy the mathematical description of various parameters of the robot. At the same time, the parameters coupling is solved and realized for the structured robot. The model is set up with sufficiently rich data interface to carry out the research and verification of trajectory planning and control algorithm. The study protocol of this paper is shown in Fig. 1. It includes the following work: (1) analyzing the kinematic structure and kinematic model of industrial robots; (2) reconstructing the CAD model of the robot; (3) establishing the dynamic model of the robot with MBS and setting up external interfaces such as controller data; (4) conducting the motion trajectory planning simulation to prove the validity and rationality of the model. 4 C. ZHANG, S. LI, Z. ZHANG Fig. 1 Research scheme for the dynamic modeling of the industrial robot arm 2. KINEMATIC ANALYSIS 2.1. Kinematic Characteristics DOBOT MG400 (hereinafter referred to as MG400) is the subject of this paper, which is an industrial super small desktop manipulator developed and released by DOBOT in December 2020. It combines the characteristics of the series mechanism and also the parallel mechanism. The quadrilateral link of the robot arm restricts its own partial working space, it has good stiffness, precision and speed. MG400 can be widely used in loading and unloading, stacking, visual sorting, distribution, viscose, automatic testing and other industrial scenarios [28-30]. The paper deals with modeling and simulation of MG400. The method used is also applicable to other similarly structured robotic arms. The motion mechanism of MG400 can be modeled as a planar mechanism with four degrees of freedom, as shown in Fig. 2. The definitions of joints J1 to J4 refer to MG400 hardware user manual (hereinafter referred to as the manual) officially published by DOBOT. The definition of each component is shown in Fig. 3. The degree of freedom F can be calculated by formula, so that n is the number of active components, PL is the number of low pairs, and PH is the number of high pairs. Then the degree of freedom F of the mechanism is given as: 3 2 3 10 2 13 0 4 L H F n P P         (1) Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 5 Fig. 2 Joints of MG400 Fig. 3 Components of MG400 In Fig. 3, the red part represents the moving component driven by J1, the blue part represents the moving component driven by J2, the yellow part represents the moving component driven by J3, and the green part is the moving component driven by J4. J1 and J4 are parallel to each other in space, and J2 and J3 are parallel to each other in space. Since there are two parallelogram mechanisms composed of L2, L2a, L2b, L33, L33a and L33b, the assembly L4 remains perpendicular to the horizontal plane, independent of the joint angles of J2 and J3. J1, J2 and J3 are orthogonal in space. J1 only affects the rotation angle between the Axz1 vertical plane (xoz plane of the J1 reference coordinate system) and the Axz0 vertical plane (xoz plane of the base coordinate system). J2 and J3 only affect the position of the execution point at the end of the robot arm on the Axz1 vertical plane. The manual does not give the dimensional parameters of each component of the robot arm, and it is challenging to obtain the dimensional parameters through actual 6 C. ZHANG, S. LI, Z. ZHANG measurement. The official support center of DOBOT released the CAD model of MG400 in November 2021. There is a certain error between the default coordinate system of L0, L1 and L4 part models and the parameters given in the manual. This error can be eliminated by creating a reference coordinate system in Solidworks. For the convenience of analysis, the parameters used in this paper are obtained by measuring the default coordinate system of the CAD model. The definition of each component reference coordinate system is shown in Fig. 4, where Ai represents the reference coordinate system of component Li during modeling (base is L0). Fig. 4 The definition of reference coordinate system for each component of MG400 When the joint angle from J1 to J4 is 0°, the conversion parameters of each coordinate system under A0 are shown in Tab. 1. The Euler angles (Rx, Ry, Rz) under the A0 reference are obtained by specifying the order of the three basic rotations. The order is x-y-z. Table 1 The conversion parameters of each coordinate system based on A0 coordinate system Coordinate system Tx (mm) Ty (mm) Tz (mm) Rx (deg) Ry (deg) Rz (deg) A0 0 0 0 0 0 0 A1 -5 0 123 0 0 0 A2 38.5 0 228 -90 0 -90 A2a -0.47 0 250.5 -90 0 -90 A2b 38.5 0 403 -90 0 0 A31 38.5 0 228 -90 0 -90 A32 -4.5 0 228 -90 0 -90 A33 -4.5 0 403 -90 0 0 A33a 67.43 0 437.47 -90 0 0 A33b 213.5 0 403 -90 0 0 A4 279.5 0 367 0 0 0 Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 7 2.2. Forward Kinematics Forward kinematics is the process of obtaining the position and pose of a point at the end of a robot arm through the known angle values of each joint, which is given as: 1 2 33 4 [ , , , ] [ , , , ] x y z z P P P R f     (2) In order to simplify the analysis, the motion mechanism of the robot arm can be equivalent to the series mechanism as shown in Fig. 5, in which the key parameters such as D-H joint reference coordinate system have been marked. Fig. 5 The equivalent D-H model of MG400 According to Fig. 5 and Tab. 1, the robot arm standard D-H parameter tab (SDH) is shown in Tab. 2, where {i} represents the number of the movable connecting rod corresponding to the D-H joint, and Link{i} represents the connecting rod corresponding to the {i} joint. The range of motion of θ2 and θ3 are mutually influenced, that is 2 3 ( ) ( 25,105)    . The joint of L33b is a driven joint, so the joint speed has no upper limit. Angle θ4 is structurally connected to L33b. The initial angle represents the angle value of θi in the initial state of the robot arm. Table 2 MG400 D-H parameter table {i} Link {i} θi (deg) di (mm) ai (mm) αi (deg) 1 L1 θ1 105 43.5 -90 2 L2 θ2 0 175 0 33 L33 θ33 0 175 0 33b L33b -θ2-θ33 0 66 90 4 L33b θ4 -36 0 0 8 C. ZHANG, S. LI, Z. ZHANG In order to express homogeneous transformation matrix conveniently, let Eq. (Fehler! Verweisquelle konnte nicht gefunden werden.) hold. According to Tab. Fehler! Verweisquelle konnte nicht gefunden werden. and SDH homogeneous transformation matrix formula [31-33], the homogeneous transformation matrix between adjacent connecting rods is easily obtained as: cos( ) sin( ) cos( ) sin( ) i i i i i j i j i j i j c s c s                   (3) The coordinates of the base center of the robot arm are defined as the base coordinate origin of the robot arm. Based on the parameters A0 and A1 in Table 1, the conversion equation from A0 to A1 is obtained: 0 1 0 0 5 0 1 0 0 0 0 1 123 0 0 0 1 B T             (4) Therefore, the homogeneous transformation matrix of the robot arm is given as: 1 4 1 4 1 2 33 2 1 4 1 4 1 2 33 20 1 2 33 33 4 0 1 2 33 33 4 2 2 33 0 ( (350 350 219)) / 2 0 ( (350 350 219)) / 2 0 0 1 175 175 192 0 0 0 1                          B B b b c s c c s s c s c s T T T T T T T c s (5) Columns 1 to 3 in 4 B T are the direction vectors of the robot arm end execution point under the reference of the base coordinate system xb (x0), yb (y0) and zb (z0) axes respectively, and column 4 is the position vector of the robot arm end execution point under the reference of the base coordinate system. From the above derivation, it can be seen that the position of the end execution point of the robot arm is only related to the three joints J1, J2 and J3, but not J4. The attitude of the end execution point of the robot arm (the attitude after rotating θ1+θ4 around zb) is only related to J1 and J4, but not J2 and J3. Assuming that the coordinates of the end execution point of the robot arm under the reference of the base coordinate system are [Px, Py, Pz] and the attitude angle is Rz, the forward kinematics equation of the robot arm is given as follows: 1 2 33 2 1 2 33 2 2 2 33 1 4 cos (175 cos( ) 175sin 109.5) sin (175 cos( ) 175sin 109.5) 175 cos 175sin( ) 192 x y z z P P P R                               (6) Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 9 2.3. Inverse Kinematics Inverse kinematics is opposite to forward kinematics. That is to calculate the angle of each joint through the known pose, which is shown as Eq. (7). 1 2 33 4 [ , , , ] [ , , , ] x y z z f P P P R     (7) (1) Calculate θ1 Let the robot arm end execution point homogeneous transformation matrix be shown as Eq. (8). 4 1 4 1 4 1 2 33 2 1 4 1 4 1 2 33 2 2 2 33 0 0 0 1 0 ( (350 350 219)) / 2 0 ( (350 350 219)) / 2 0 0 1 175 175 192 0 0 0 1 x x x x y y y yB z z z z n o a P n o a P T n o a P c s c c s s c s c s c s                                     (8) Multiply Eq. (8) left and right by 0 1 1 1 0 B T T   . Then we get Eq. (9). 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 4 4 2 33 2 2 33 2 4 4 43.5 228 0 0 0 1 0 175 175 66 0 0 1 175 175 36 0 0 0 0 0 1 x y x y x y x y z z z z y x y x y x y x n c n s o c o s a c a s p c p s n o a p n c n s o c o s a c a s p c p s c s c s s c s c                                           (9) Since Eq. (9) holds, the elements on the left and right of the equation must be equal. According to the position elements of matrix (3, 4) in Eq. (9), we get 1 1 0 y x p c p s  (10) Therefore, 1 tan y x p p   (11) Eq. (11) can be further written as: 1 atan2( , ) y x p p  (12) 10 C. ZHANG, S. LI, Z. ZHANG (2) Calculate θ3 According to the matrix (1, 4) and (2, 4) position elements in Eq. (9), we get Eq. (13). 1 1 2 33 2 2 33 2 43.5 175 175 66 228 175 175 36 x y z p c p s c s p s c             (13) Therefore, we have Eq. (14). 2 2 2 1 1 33 2 2 175 ( 109.5) (192 ) sin 2 175 x y z p c p s p          (14) According to the initial angle of θ33 and the height relationship between the level of A33b and the level of A2, we have Eq. (15). 33 33 1 33 33 33 1 1 sin (3, 4) (3, 4) cos 1 sin (3, 4) (3, 4) B B B B T T T T            (15) In Eq. (15), 33 B T is the homogeneous transformation matrix from A0 to A33b, 1 B T is the homogeneous transformation matrix from A0 to A2, then we have Eq. (16). 33 33 33 atan2(sin , cos )   (16) (3) Calculate θ2 Transforming Eq. (13), we obtain Eqs. (17-19): 2 2 1 1 33 2 33 33 2 33 ( 109.5)( 1) 175 ( 1) (192 ) 175        x y z p c p s s s s p c s c (17) 2 2 2 33 2 33 33 1 1 33 175 ( 1) 175 (192 ) ( 109.5)( 1)        z x y s s s c p c p c p s s (18) 2 2 2 33 33 33 1 1 33 (175( 1) 175 ) (192 ) ( 109.5)( 1)        z x y s s c p c p c p s s (19) By further transformation, the following equations are obtained: 33 1 1 33 2 2 2 33 33 (192 ) ( 109.5)( 1) sin 175( 1) 175 z x y p c p c p s s s c          (20) 2 2 cos 1 sin   (21) Then we have Eq. (22). 2 2 2 atan2(sin , cos )   (22) (4) Calculate θ4 According to the matrix (1,1) and (2,1) position elements in Eq. (Fehler! Verweisquelle konnte nicht gefunden werden.), we get the following equation: 1 4 1 4 x y n c n s      (23) Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 11 Since atan2( , ) z y x R n n (24) According to Eq. (6) and Eq. (24), the following equation is obtained: 4 1 atan2( , ) y x n n   (25) The influence of the initial state of the robot arm joint and the motion trajectory of the corresponding initial state is not considered, because it solves the minimum angle value of θ4 when the robot arm is in a certain pose, and we have Eq. (26). 1 4 4 1 4 atan2( , ) 360 180 atan2( , ) 360 180 y x y x n n n n                 (26) According to the above calculation, the inverse kinematics equation of the robot arm is shown as Eq. (27). 1 33 2 2 2 1 1 2 33 2 33 1 1 33 2 2 33 33 2 4 1 atan2( , ) atan2 2 175 ( 109.5) (192 ) ( 2 175 , 1 sin ) atan2 (192 ) ( 109.5)( 1) ( 175( 1) 175 , 1 sin ) atan2( , ) 360 y x x y z z x y y x p p p c p s p p c p c p s s s c n n                                          (27) 2.4. Kinematic Verification In order to verify the correctness of Eq. (6) and Eq. (27), the kinematics solution function is written by MATLAB. Ten groups of joint angles that conform to the motion range of the joints of the robot arm are randomly generated, and Eq. (6) is used to calculate the forward kinematics solution. Then, according to the obtained position and posture of the terminal execution point, the inverse kinematics solution is obtained by using Eq. (27). The correctness of Eq. (6) and Eq. (27) is verified by comparing the results of the two calculations. The specific results of forward and inverse kinematics solutions are shown in Tab. 3 and Tab. 4. 12 C. ZHANG, S. LI, Z. ZHANG Table 3 Forward kinematics conditions and calculation results of MG400 No. θ1(deg) θ2(deg) θ3(deg) θ4(deg) Px(mm) Py(mm) Pz(mm) Rz(deg) 1 -1.86 60.70 7.26 290.68 327.62 -10.66 115.45 -71.19 2 125.10 11.76 48.32 -217.58 -133.65 190.19 211.66 -92.48 3 -150.23 56.85 -16.85 -14.46 -338.58 -193.69 175.21 -164.68 4 129.51 42.09 13.21 258.80 -207.68 251.84 178.01 28.31 5 97.76 38.44 -39.66 -187.25 -53.08 389.66 332.80 -89.49 6 123.68 -21.85 40.03 -239.09 -116.82 175.27 299.81 -115.41 7 153.18 53.40 -13.34 -20.82 -342.62 173.24 183.72 132.36 8 -140.92 50.02 -69.50 -308.56 -317.17 -257.55 362.82 -89.48 9 6.93 -14.36 67.61 228.63 169.56 20.60 221.32 -124.44 10 71.18 -8.52 50.46 13.39 68.96 202.33 248.11 84.57 Table 4 Inverse kinematics conditions and calculation results of MG400 No. Px(mm) Py(mm) Pz(mm) Rz(deg) θ1(deg) θ2(deg) θ3(deg) θ4(deg) 1 327.62 -10.66 115.45 -71.19 -1.86 60.70 7.26 -69.32 2 -133.65 190.19 211.65 -92.48 125.10 11.76 48.32 142.42 3 -338.58 -193.69 175.21 -164.68 -150.23 56.85 -16.85 -14.46 4 -207.68 251.84 178.01 28.31 129.51 42.09 13.21 -101.20 5 -53.08 389.66 332.80 -89.49 97.76 38.44 -39.66 172.75 6 -116.82 175.27 299.81 -115.41 123.68 -21.85 40.03 120.91 7 -342.62 173.24 183.72 132.36 153.18 53.40 -13.34 -20.82 8 -317.17 -257.55 362.82 -89.48 -140.92 50.02 -69.50 51.44 9 169.56 20.60 221.37 -124.44 6.93 -14.36 67.61 -131.37 10 68.96 202.33 248.11 84.57 71.18 -8.51 50.46 13.39 When solving the inverse solution of the robot arm, the minimum value of the rotation angle of each joint of the robot arm is obtained. When the rotation angle of θ4 exceeds the range of [-pi, pi], the algorithm will let θ4 map to the range of [-pi, pi], and select a smaller inverse solution under the condition of ensuring the attitude of the robot arm. According to Tab. 3 and Tab. 4, it can be concluded that the kinematic analysis represented by Eq. (6) and Eq. (27) is correct. 3. CONSTRUCTION OF MBS MODEL 3.1. Build Solid Module Mechanism of MBS Model To establish the solid module mechanism of the MBS model of the robot arm, it requires understanding the parameters of each part of the robot arm, but it is difficult to directly obtain the specific parameters of each part of the robotic arm. Therefore, the relevant parameters in this article are generally obtained through CAD model analysis. The Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 13 solid module mechanism of the established MBS model is shown in Fig. 6 (a). Solidworks is used to establish the reference coordinate system shown in Tab. Fehler! Verweisquelle konnte nicht gefunden werden. in the official CAD model, and all Solidworks related entities are output as step format files under the corresponding reference coordinate system. (a) (b) Fig. 6 a) Solid module mechanism of MBS model and b) its workspace The solid module model is used for the dynamic simulation of the system, and the working space of each joint needs to be accurately consistent with the real robot arm, so that the various parameters of the model can be directly applied to the real robot arm. The workspace of the robot arm model system is obtained in joint space as shown in Fig. 6 (b). Comparing with the MG400 manual, it is obtained that the model established can accurately simulate the working state and parameters of the actual robot, and is reliable and feasible. 3.2. Build Robot Arm MBS Model A robot arm MBS system has been established, mainly including trajectory planning, controller, robot MBS model, and simulation data module [34-36]. The MBS model is composed of joint module and environment module. The joint module is divided into rotation and translation according to the solid module motion mode, as shown in Fig. 7 which includes solid module coordinate system transformation, joints and part solid modules. The joint module is used to simulate the motion state of each joint of the robot arm and the components connected on the joint, and can also simulate the motion state of the outer shaft of the robot arm. The environment module is shown in Fig. 8. 14 C. ZHANG, S. LI, Z. ZHANG Fig. 7 Joint module of the MBS model Fig. 8 Environment module of the MBS model In order to make the model consistent with the equivalent model, the position of L33 rotating joint is coincident with A2b. This causes L2a, L2b, L31, L32, L33a and L33b in the robot arm model to be driven components. According to Tab. 1 and Tab. 2, the solid transform parameters of each joint module can be obtained, as shown in Tab. 5. Table 5 The solid transform parameters of joint module Coordinate System Conversion Module Tx (mm) Ty (mm) Tz (mm) Rx (deg) Ry (deg) Rz (deg) RT 0_1 -5 0 123 0 0 0 RT 1_2 43.5 0 105 -90 0 -90 RT 1_2a -0.47 0 250.5 -90 0 -90 RT 2_2b 175 0 0 0 0 90 RT 1_31 43.5 0 105 -90 0 -90 RT 31_32 0 -43 0 0 0 0 RT 2_33 175 0 0 0 0 90 RT 2b_33a 28.93 -34.47 0 0 0 0 RT 33_33b 175 0 0 0 0 0 RT 33b_4 66 36 0 90 0 0 In Tab. 5, RT {i}_{j} represents the conversion parameters from A{i} to A{j}, and {i} and {j} represent the number of the robot arm assembly and the corresponding coordinate system. In this model, the corresponding joint modules of L1, L2, L33 and L4 are used as the Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 15 active components of the robot arm, and the input can be the target angle (j{i}) or torque (t{i}). The visualization of the robot arm is obtained by importing the parts of models shown in Fig. 6 into the solid module. As the MG400 is made of aluminum alloy, the model parameter density is set to 2700kg/m3, and the physical properties of the components such as mass, center of mass coordinates and moment of inertia will be automatically calculated according to the geometric characteristics of the model, as shown in Tab. 6. Table 6 The physical properties of solid module modules Solid Module Mass (kg) Centroid coordinates (cm) Moment of inertia (g·cm2) L1 4.781 [0.462, -0.006, 5.699] [163602, 156446, 200109] L2 0.234 [5.963, -0.286, -0.049] [1049.56, 9144.35, 8310.61] L2a 0.179 [7.378, -1.967, 0.001] [1009.9, 5914.17, 5415.16] L2b 0.099 [-0.511486, -2.702, 0.002] [408.52, 663.483, 634.498] L31 0.012 [-0, -1.728, -0.003] [43.717, 3.788, 45.756] L32 0.025 [8.934, 0.204, -0.453] [3.752, 966.085, 968.624] L33 0.101 [8.21, 0.487, 0.011] [68.215, 4436.23, 4433.04] L33a 0.122 [8.78, -0.668, -0.001] [541.489, 3825.9, 3563.67] L33b 0.431 [5.526, 0.337, 0.003] [2676.58, 2873.23, 4239.39] L4 0.004 [0, 0, 0] [3.824, 3.824, 0.223] Connecting each module which are shown in Fig. 2 to get the MBS model of the robot arm, as shown in Fig. 9. The physical quantities in MBS model are simulated by its basic equations, which have original physical effects and characteristics. The motion state of the system is solved by the basic motion equation at each sampling time, and the visual images are synchronously generated in mechanics explorers. Fig. 9 The MBS model of robotic arm 16 C. ZHANG, S. LI, Z. ZHANG 3.3. MBS Dynamic Simulation Modeling The MBS dynamic simulation model after adding the trajectory planning module, control algorithm module and data acquisition module is shown in Fig. 10. These interfaces all greatly facilitate the follow-up research and development of the robot. In addition, the MBS model retains the interface between the robot arm and the robot outer axis and the tool, enabling the simulation system to be used repeatedly and effectively in subsequent studies of the robot arm, in this way to facilitate subsequent ongoing development. Fig. 10 The dynamic simulation model of MG400 4. SIMULATION RESULTS AND ANALYSIS The dynamic simulation of the MBS model system by the above. In this paper, the palletizing operation task of the industrial mechanical arm is selected as the trajectory planning target. According to the operation requirements, given the coordinates of Target points such as Initial point (A), Lift hand point (B1, B2) , Drop hand point (C1, C2) and Task point (D1, D2). The schematic diagram of the movement of the robotic arm operation is shown in Fig. 11. Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 17 Fig. 11 The schematic diagram of the movement of the robotic arm operation After setting the target points for J1, J2, J33 and J4 joints, trajectory planning simulations were performed using three methods: the Quintic Polynomial Trajectory, the B-Spline and the Minimum Jerk Polynomial Trajectory. Generate the target point parameters for the robotic arm, and the driving signal of each joint module is obtained by solving the inverse solution. Then the signal is directly transmitted to the MBS model to obtain the real-time torque or force. All simulation data is available via the rotating joints and transformation sensor modules. Similarly, a control algorithm module, added to the system, using torque as a drive signal, can be transmitted to the motion controller module to verify the effect of the control algorithm. The paths and their process dynamic maps planned by the above three trajectory planning algorithms are obtained, as shown in Figs. 12-14. After the path planning, 3D dynamic animation can also be generated in the mechanical explorer, which provides a nice and straightforward option to observe the trajectory. A few screenshots of the obtained trajectories are shown in Fig. 15. (a) (b) Fig. 12 Path and its dynamic map planned by Quintic Polynomial Trajectory 18 C. ZHANG, S. LI, Z. ZHANG (a) (b) Fig. 13 Path and its dynamic map planned by B-Spline (a) (b) Fig. 14 Path and its dynamic map planned by Minimum Jerk Polynomial Trajectory Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 19 Fig. 15 3D simulation visualization of robot arm The obtained joint angles, angular velocities and angular accelerations are shown in Fig. 16, Fig. 17 and Fig. 18, respectively, and the output torques are shown in Fig. 19, including Quintic Polynomial Trajectory, B-Spline, and Minimum Jerk Polynomial Trajectory. As can be observed from these figure, in the motion data obtained from the model, there is an obvious proportional relationship between the torque and acceleration of J1. J2 and J33, which are greatly affected by gravity, and J4 is not equipped with any fixtures, so its torque is almost zero. It can be obtained by comprehensive comparison that the Quintic Polynomial Trajectory is an algorithm involving both the starting point and the end point, yielding the smoothest path between the two points. The B-Spline can be inserted into multiple pathway points, those on the path between the start and the end. The path generated by the B-Spline algorithm does not pass through the path points. It has the advantages of a smooth curve and short path, but it will cause discontinuity in acceleration and torque. The Minimum Jerk Polynomial Trajectory can insert several pathway points between the starting point and the end point. Unlike the B-Spline, the path planned by this method is through the pathway point. The advantage is that the path is the smoothest, while ensuring the complete continuity of angle, angular velocity, angular acceleration and angular acceleration. The above data well simulate the actual operation of the robot arm. 20 C. ZHANG, S. LI, Z. ZHANG Fig. 16 Movement angle of each joint of MBS model Fig. 17 Angular velocity of each joint of MBS model Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 21 Fig. 18 Angular acceleration of each joint of MBS model Fig. 19 Torque of each joint of MBS model 22 C. ZHANG, S. LI, Z. ZHANG The error between the actual simulated trajectory and the ideal trajectory at the end execution point of the robot arm model is shown in Fig. 20 and the error statistics are shown in Tab. 7. The average error of simulation is as low as 10-11 mm, and the minimum error is 0, which is a very good accuracy. To sum up, the developed robot arm model can simulate well the actual kinematics and dynamics of the robot arm. Fig. 20 Error between simulation trajectory and ideal trajectory of robot arm Table 7 Error between simulation trajectory and ideal trajectory of the end execution point of robot arm MBS model Method Maximum Error (mm) Minimum Error (mm) Average Error (mm) Quintic Polynomial Trajectory 3.8730e-09 0 9.0446e-11 B-Spline 2.1104e-09 0 4.4756e-11 Minimum Jerk Polynomial Trajectory 4.4408e-09 0 1.2206e-10 5. CONCLUSION The MBS modeling and Simulation of DOBOT robot arm are carried out by using Simscape Multibody in this paper. First, the kinematic properties of the robot arm are analyzed and the forward and inverse kinematic equations of the robot arm are established. Secondly, Similink is used to establish the MBS model logic framework of the robot arm, and then the solid model is imported into the corresponding solid module in the model, calculate the corresponding physical properties of the corresponding components, and obtain the visualization of the robot arm. Third, an offline dynamic simulation model of the robot arm was established by adding a trajectory planning module, a control algorithm module, and a data acquisition module to the MBS model. Fourth, the motion trajectory of Modeling and Simulation of Industrial Robot Arms Using Simscape Multibody 23 the simulation model is planned and simulated in the joint space. The simulation results show that the simulation model can accurately simulate the motion state of the robot arm, the visual simulation animation is easily obtained to intuitively observe the motion process of the robot arm, and the output simulation data is convenient for the study and analysis of the subsequent control algorithm. The peripheral mechanisms can also be extended by adding joint modules or related follow-up studies by embedding the code into the trajectory planning module or control algorithm module. Based on the dynamic simulation function of the simulation model, it is more economical and easier to study the robot arm with the same structure as MG400. Acknowledgement: The paper is supported by the 'Three Levels' Talent Construction Project of Zhuhai College of Science and Technology, the Higher Education Teaching Reform Project of Zhuhai College of Science and Technology (ZLGC20220529) and the Science and Technology Projects in Guangzhou (202102021282). REFERENCES 1. Kouritem, S.A., Abouheaf, M.I., Nahas, N., Hassan, M., 2022, A Multi-Objective Optimization Design of Industrial Robot Arms, Alexandria Engineering Journal, 61(12), pp. 12847-12867. 2. Damaševičius, R., Maskeliūnas, R., Narvydas, G., Narbutaitė, R., Połap, D., Woźniak, M., 2020, Intelligent Automation of Dental Material Analysis Using Robotic Arm with Jerk Optimized Trajectory , Journal of Ambient Intelligence and Humanized Computing, 11(12), pp. 6223-6234. 3. Kim, J., Croft, E.A., 2019, Online near Time-Optimal Trajectory Planning for Industrial Robots, Robotics and Computer-Integrated Manufacturing, 58, pp. 158-171. 4. Rout, A., BBVL, D., Biswal, B.B., Mahanta, G.B., 2021, Optimal Trajectory Planning of Industrial Robot for Improving Positional Accuracy, Industrial Robot-the International Journal of Robotics Research and Application, 48(1), pp. 71-83. 5. Özyer, B., 2020, Adaptive Fast Sliding Neural Control for Robot Manipulator, Turkish Journal of Electrical Engineering and Computer Sciences, 28(6), pp. 3154-3167. 6. Yilmaz, B.M., Tatlicioglu, E., Savran, A., Alci, M., 2021, Adaptive Fuzzy Logic with Self-Tuned Membership Functions Based Repetitive Learning Control of Robotic Manipulators, Applied Soft Computing, 104, pp. 107-183. 7. Gharib, M.R., 2020, Comparison of Robust Optimal QFT Controller with TFC and MFC Controller in a Multi-Input Multi-Output System, Reports in Mechanical Engineering, 1(1), pp. 151-61. 8. Rojas-Garcia, L., Bonilla-Gutierrez, I., Mendoza-Gutierrez, M., Chavez-Olivares, C., 2022, Adaptive Force/Position Control of Robot Manipulators with Bounded Inputs, Journal of Mechanical Science and Technology, 36(3), pp. 1497-1509. 9. Gutierrez-Giles, A., Evangelista-Hernandez, L.U., Arteaga, M. A., Cruz-Villar, C.A., Rodriguez-Angeles, A., 2021, A Force/Motion Control Approach Based on Trajectsory Planning for Industrial Robots with Closed Control Architecture, IEEE Access, 9, pp. 80728-80740. 10. Urrea, C., Pascal, J., 2021, Design and Validation of a Dynamic Parameter Identification Model for Industrial Manipulator Robots, Archive of Applied Mechanics, 91(5), pp. 1981-2007. 11. Fernández, J.G., Yu, B.B., Bargsten, V., Zipper, M., Sprengel, H., 2020, Design, Modelling and Control of Novel Series-Elastic Actuators for Industrial Robots, Actuators, 9(1), pp. 414-420. 12. Bugnar, F., Nedezki, C.M., Moholea, I.F., 2019, Dynamic Model of 6R Articulated Industrial Robot Using Newton-Euler Formulation, Acta Technica Napocensis Series-Applied Mathematics Mechanics and Engineering, 62(1), pp. 169-174. 13. Urrea, C., Pascal, J., 2021, Dynamic Parameter Identification Based on Lagrangian Formulation and Servomotor-Type Actuators for Industrial Robots, International Journal of Control Automation and Systems, 19(8), pp. 2902-2909. 14. Šulekić, I., Osman, K., Lavrić, Z., 2022, Implementation of a Mathematical Model and Comparison of Control Algorithms for an Industrial Robot Arm, Polytechnic and design, 10(3), pp. 211-221. 24 C. ZHANG, S. LI, Z. ZHANG 15. Ho, T., Kang, C.G., Lee, S., 2012, Efficient Closed-Form Solution of Inverse Kinematics for a Specific Six-Dof Arm, International Journal of Control Automation and Systems, 10(3), pp. 567-73. 16. Bien, D.X., 2023, Investigation of Driving Torques at the Joints of Industrial Robot Arms Based on the Topology Optimization Technique, Journal of Applied and Computational Mechanics, 9(3), pp. 820-833. 17. Wang, J.W., Xu, J.X., 2021, Kinematic Modeling and Simulation of Dual-Arm Robot, Journal of Robotics, Networking and Artificial Life, 8(1), pp. 56-59. 18. Wang, J.W., Xu, J.X., 2021, Kinematic Modeling and Simulation of Humanoid Dual-Arm Robot, Proc. 26th International Conference on Artificial Life and Robotics, 26, pp. 381-384. 19. Alam, M.M., Ibaraki, S., Fukuda, K., 2021, Kinematic Modeling of Six-Axis Industrial Robot and Its Parameter Identification: A Tutorial, International Journal of Automation Technology, 15(5), pp. 599-610. 20. Nguyen, V., Marvel, J.A., 2022, Modeling of Industrial Robot Kinematics Using a Hybrid Analytical and Statistical Approach, Journal of Mechanisms and Robotics-Transactions of the Asme, 14(5), 051009. 21. Harfoush, A., Hossam, M., 2020, Modelling of a Robot-Arm for Training in Fencing Sport, International Journal of Intelligent Robotics and Applications, 4(1), pp. 109-121. 22. Huynh, H.N., Assadi, H., Riviere-Lorphevre, E., Verlinden, O., Ahmadi, K., 2020, Modelling the Dynamics of Industrial Robots for Milling Operations, Robotics and Computer-Integrated Manufacturing, 61, 101852. 23. Segota, S.B., Andelic, N., Car, Z., Sercer, M., 2022, Prediction of Robot Grasp Robustness Using Artificial Intelligence Algorithms, Tehnicki Vjesnik-Technical Gazette, 29(1), pp. 101-107. 24. Truc, L.N., Nguyen, T.L., 2020, Quasi-Physical Modeling of Robot IRB 120 Using Simscape Multibody for Dynamic and Control Simulation, Turkish Journal of Electrical Engineering and Computer Sciences, 28(4), pp. 1949-1964. 25. Barrientos-Diez, J., Dong, X., Axinte, D., Kell, J., 2021, Real-Time Kinematics of Continuum Robots: Modelling and Validation, Robotics and Computer-Integrated Manufacturing, 67, 102019. 26. Ogbemhe, J., Mpofu, K., Mokakabye, M., 2022, Robot Dynamic Model: Freudenstein-Based Optimal Trajectory and Parameter Identification, Cogent Engineering, 9(1), 2046682. 27. Aboulissane, B., Bakkali, L.E., Bahaoui, J.E., 2020, Workspace Analysis and Optimization of the Parallel Robots Based on Computer-Aided Design Approach, Facta Universitatis-Series Mechanical Engineering, 18(1), pp. 79-89. 28. Brito, T., Lima, J., Braun, J., Piardi, L., Costa, P., 2020, A Dobot Manipulator Simulation Environment for Teaching Aim with Forward and Inverse Kinematics, Proc. 14th APCA International Conference on Automatic Control and Soft Computing, Braganca, PORTUGAL, vol. 695, pp. 303-312. 29. Miki, K., Nagata, F., Furuta, K., Arima, K., Shimizu, T., Ikeda, T., Kato, H., Watanabe, K., Habib, M. K., 2022, Development of a Hyper Cls Data-Based Robotic Interface for Automation of Production-Line Tasks Using an Articulated Robot Arm, Artificial Life and Robotics, 27(3), pp. 547-553. 30. Liri, F., Luu, A., George, K., Angulo, A., Dittloff, J., 2022, Real-Time Dynamic Object Grasping with a Robotic Arm: A Design for Visually Impaired Persons, Proc. IEEE World AI IoT Congress, Seattle, WA, pp. 344-350. 31. Guo, F.Y., Cai, H., Ceccarelli, M., Li, T., Yao, B.T., 2020, Enhanced D-H: An Improved Convention for Establishing a Robot Link Coordinate System Fixed on the Joint, Industrial Robot-the International Journal of Robotics Research and Application, 47(2), pp. 197-205. 32. Singh, A., Singla, A., Soni, S., 2015, Extension of D-H Parameter Method to Hybrid Manipulators Used in Robot-Assisted Surgery, Proceedings of the Institution of Mechanical Engineers Part H-Journal of Engineering in Medicine, 229(10), pp. 703-712. 33. Gupta, K.C., 1986, Kinematic Analysis of Manipulators Using the Zero Reference Position Description, International Journal of Robotics Research, 5(2), pp. 5-13. 34. Li, Y., Wang, J., Ji, Y.M., 2020, Function Analysis of Industrial Robot under Cubic Polynomial Interpolation in Animation Simulation Environment, International Journal of Interactive Multimedia and Artificial Intelligence, 6(4), pp. 105-112. 35. Ozakyol, H., Karaman, C., Bingul, Z., 2019, Advanced Robotics Analysis Toolbox for Kinematic and Dynamic Design and Analysis of High-Dof Redundant Serial Manipulators, Computer Applications in Engineering Education, 27(6), pp. 1429-1452. 36. Baglioni, S., Cianetti, F., Braccesi, C., Micheli, D.M., 2016, Multibody Modelling of N Dof Robot Arm Assigned to Milling Manufacturing. Dynamic Analysis and Position Errors Evaluation, Journal of Mechanical Science and Technology, 30(1), pp. 405-420.