# Minimum Energy Trajectory Optimization for Driving Systems of Palletizing Robot Joints

Authors: Ying He; Jiangping Mei; Zhiwei Fang; Fan Zhang; Yanqin Zhao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Palletizing robot is widely used in logistics operation. At present, people pay attention to not only the loading capacity and working efficiency of palletizing robots, but also the energy consumption in their working process. This paper takes MD1200-YJ palletizing robot as the research object and studies the problem of low energy consumption optimization of joint driving system from the perspective of trajectory optimization. Firstly, a multifactor dynamic model of palletizing robot is established based on the conventional inverse rigid body dynamic model of the robot, the Stribeck friction model and the spring balance torque model. And then based on joint torque, friction torque, motion parameter, and joule’s law, the useful work model, thermal loss model of joint motor, friction energy consumption model of joint system, and total energy consumption model of palletizing robot are established, and through simulation and experiment, the correctness of the multifactor dynamic model and energy consumption model is verified. Secondly, based on the Fourier series approximation method to construct the joint trajectory expression, the minimum total energy consumption as the optimization objective, with coefficients of Fourier series as optimization variables, the motion parameters of initial and final position, and running time constant as constraint conditions, the genetic algorithm is used to solve the optimization problem. Finally, through the simulation analysis the optimized Fourier series motion law and the 3-4-5 polynomial motion law are comprehensively evaluated to verify the effectiveness of the optimization method. Moreover, it provides the theoretical basis for the follow-up research and points out the direction of improvement.


---

## Body

## 1. Introduction

Palletizing robots have been widely applied in logistics operations such as storage, handling, and transportation of materials. At the present, energy saving and green manufacturing are emphasized. Energy consumption index of mechanical equipment has become a new standard to judge the advanced technology of industrial equipment [1, 2]. The energy consumption of industrial robots is approximately 8% of the total electrical energy consumed in production processes. Therefore, besides repeatability, accuracy, speed, efficiency, and load capacity of industrial robots, energy consumption indicator which is in the manufacturing and operation of industrial robots process also has been the focus of attention [3–5].There are many researches on reducing the energy consumption of robots. The main methods include optimization of robot structure [6–8], trajectory planning and so on. Optimization of robot structure method is especially suitable for robot structure design stage. Trajectory planning is more suitable to reduce the energy consumption of existing robots. This paper takes MD1200-YJ high speed and heavy load palletizing robot developed by Tianjin University in China as the research object and researches the minimum energy consumption of joint drive system which is the core of driving the manipulator to complete the motion control, and it is also the main link of energy consumption, from the perspective of trajectory planning.The driving torque of each joint drive system servo motor is an important factor affecting its energy consumption, and different motion laws will result in different joint driving torque under the same load conditions. Robot dynamics analysis is a theoretical method to determine the dynamic relationship between robot motion and joint drive torque. The commonly used methods of robot dynamics analysis include Newton-Euler method, Lagrange method, Kane equation method, and so on [9–12]. However, the conventional inverse rigid body dynamics model only contains the inertial term, gravity term, and centrifugal or Coriolis force term. And the energy consumption model derived from this dynamic model belongs to the useful work model, which is defined as the work done by the joints of robot in the process of carrying load in this paper. Nevertheless, in addition to the useful work, there are also various dissipative factors such as joint friction energy consumption and servo motor heat loss in the energy consumption of robot joint drive system, which cannot be ignored in energy consumption research. In addition, in order to reduce the load on the joint parts of heavy load robot, the balancing device, such as the balancing spring cylinder, is generally designed, which can effectively reduce the peak torque of the large arm drive joint [13], but this spring torque is not considered in the conventional rigid body dynamics model. Therefore, the conventional rigid body dynamic model cannot fully meet requirements of this research. A multifactor dynamic model of palletizing robot needs to be established; it is necessary to add friction torque and spring balance torque components on the conventional rigid body dynamic model basis. And combined with joule’s law, the heat loss model and energy consumption model of multijoint motor of palletizing robot necessary be established.Literatures [14, 15] point out that about 20 percent of the robot’s energy is used to overcome joint friction; therefore, friction is an important factor to consider when researching energy consumption problem of robot [16]. Friction model is a mathematical description of friction behavior and also is a theoretical basis for correct understanding of friction mechanism and effective prediction of friction behavior. The friction behavior is very complex, and it has been the focus of many scholars to build a friction model that can fully describe the frictional behavior. Up to now, there have been more than 30 friction models. Friction models are usually divided into two types: static friction models and dynamic friction models. The static friction model describes friction as a function of relative velocity, including Coulomb model, static friction-Coulomb-viscosity-Stribeck model, Karnopp friction model, and so on. The dynamic friction model describes the friction as a function of relative velocity and displacement and can describe comprehensively static and dynamic characteristics of friction. The classical dynamic friction model mainly includes Dahl model, mane model, LuGre model, generalized Maxwell model, and so on [17]. Compared with the static friction model, the dynamic friction model can describe the friction behavior more comprehensively, but the dynamic friction model is complex, and the difficulty of parameter identification is also significantly increased. Therefore, it should be reasonably selected according to the specific scope and precision of the research. Considering that this paper mainly researches the energy consumption of palletizing robot during high speed operation, the frictional influence in the dynamic process of low speed unsteady state can be ignored. Therefore, the friction model of each joint of palletizing robot is established based on the Stribeck model.Parameter identification of friction model is the process of collecting and processing relevant data by experimental method to obtain the unknown parameters of friction model. The accuracy of parameter identification determines the accuracy of friction model. Parameter estimation often uses least square method [18, 19].As for the optimization of minimum energy consumption of robot joint drive system, scholars at home and abroad have made some explorations and researches. The main method is to optimize the motion rule with the minimum energy consumption as the objective [20]. Wang zhe and Mei jiangping, respectively, researched the influence of 2 and 3 DOF parallel manipulator operating space and joint space trajectory planning on system energy consumption [21, 22]. Paes K et al. [4] to ABB’s IRB1600 industrial robot as the research object and proposes a scene recognition and time and energy optimal trajectory planning method, to work space, joint speed acceleration and dynamic parameters as constraint conditions, the minimum cost function based on the energy as the optimization goal. The sequential quadratic programming method was used to solve the problem and achieved good effect. Wigstrom O et al. [23] proposed a method to generate energy optimal trajectories using dynamic time scales. Without changing the time sequence and trajectory of the robot system, the optimal control time sequence of robot energy is obtained through dynamic programming. Compared with the traditional linear time scale, the energy consumption is significantly reduced. The above optimization strategies are to reduce the useful work by trajectory planning the robot. However, the dissipation energy, such as friction and heat loss, widely exists in the robot joint system should not be ignored. Some scholars have also studied the dissipation energy of robot joint system. Jin bo et al. [24] established the six-legged walking robot dynamics model. Based on the useful work and heat loss of the motor, combined with the quadratic programming method, the optimal distribution of each joint torque was carried out, and the energy consumption of the system was reduced by 22% to 39%. Pellegrinelli S et al. [25] researched the analysis and optimization of the energy consumption related to auxiliary robotic assembly processes, a frictional model considering the coulomb-viscous friction and inverse dynamic model were established, and an energy consumption model including constant power mechanical power and copper loss of motor was established, a method for automatic generation of collision-free trajectories was proposed to reduce the energy of the work cycle by 12%. Datouo R et al. [26] proposed an optimal motion planning method aiming to achieve minimum energy consumption for a three-wheeled omnidirectional mobile robot (TOMR) named Robotino. The A∗ path planner and the heuristic function integrating energy saving criterion were used to generate energy-saved paths. And the effectiveness of the proposed motion planning was demonstrated by performing a series of simulations. Izumi T [27] researched the energy consumption of the robot’s joint servo drive system under no-load conditions, clarified the effect of Coulomb friction, and proved the linear relationship between the optimal angular velocity cross time and the minimum energy consumption. Most of the above researches are focused on lightweight, light-loaded, and less joint mechanical arm. In this paper, the MD1200-YJ is a high speed and heavy load palletizing robot; it involved 4 joints and contains multiple constraint branches, which involves more factors, more variables, and more complicated problems.In this paper, the minimum energy trajectory optimization of the joint drive system of MD1200-YJ palletizing robot is researched. Firstly, based on the conventional inverse rigid body dynamic model of the robot, the Stribeck friction model, and the spring balance torque model, a multifactor dynamic model of the palletizing robot is established. And then by combining joint torque, friction torque, and motion parameter with joule’s law, the useful work model, thermal loss model of joint motor, friction energy consumption model of joint system, and total energy consumption model of palletizing robot are established. And through simulation and experiment, the correctness of the multifactor dynamic model and energy consumption model is verified. Secondly, the joint motion law is constructed by using the Fourier series approximation method, and taking the coefficients of Fourier series as optimization variables, the minimum energy consumption is taken as the optimization objective, the motion parameters of initial and final position and running time constant as constraint conditions, and the minimum energy consumption trajectory of the joint system is optimized by the genetic algorithm. Finally, the effectiveness of the optimization method is verified by simulation analysis.


## 2. Establishment of Energy Consumption Model of MD1200-YJ Palletizing Robot

Based on Kane equation method and the structural characteristics of MD1200-YJ palletizing robot, its rigid body dynamics model is established. Based on the Stribeck friction model and the parameter identification theory, the friction torque models of each joint of the robot are established. Based on the motion parameters of the robot shoulder joint combined with the structural parameters of the balance spring cylinder, the balance spring torque model is established. The multifactor dynamics model of MD1200-YJ palletizing robot is established based on the above model. The energy consumption model of the robot is established by the joule’s law by combining the multifactor dynamics model with the motion parameters. And the accuracy of the multifactor dynamics model and the energy consumption model is verified by experiments.


### 2.1. Introduction to the Mechanism of MD1200-YJ Palletizing Robot

The MD1200-YJ palletizing robot is a 4-DOF multijoint robot, load capacity: 120kg, maximum turning radius: 2400mm, repeat positioning accuracy: ±0.4mm, waist maximum rotation speed: 85°/s, belonging to high speed and heavy load palletizing robot, its 3D model shown in Figure1. It is mainly composed of machine base, waist joint, small arm joint, drive arm of small arm, balance spring cylinder, balancing weight, small arm drive link, attitude hold link1, attitude hold tripod, attitude hold link2, small arm link, big arm joint, big arm link, waist mounting bracket, end effector mounting bracket, and so on. The action of two parallelogram mechanical linkages (one consists of attitude hold link1, big arm link, attitude hold tripod, and waist mounting bracket, and the other consists of attitude hold tripod, attitude hold link2, small arm link, and end effector mounting bracket) makes undersurface of end effector keep parallel to the ground.Figure 1
3D model of MD1200-YJ palletizing robot.


### 2.2. Establishment of Multifactor Dynamic Model of MD1200-YJ Palletizing Robot

On the basis of the conventional rigid body inverse dynamics model of robot, the multifactor inverse dynamics model of palletizing robot with considering joint friction torque and spring balance torque is established.


#### 2.2.1. Rigid Body Dynamics Model

Rigid body dynamics model for MD1200-YJ palletizing robot is derived using Kane equation method in this paper. The basic idea of Kane equation method is to replace the generalized coordinate with the generalized velocity as the independent variable of the system. When using Kane method for dynamic modeling, it only involves the cross product operation of the vector’s dot product, which greatly improves the calculation efficiency. Meanwhile, the speed and acceleration of connecting link parts solution have the recursive nature, which is convenient for the computer to carry out numerical calculation.Literature [12, 13] described derivation ideas of Kane equation method in detail, and here we will not repeat them, but directly give the derivation results.Assuming thatτi represents the drive torque of jointi, and the rigid body inverse dynamics equation of the robot system is obtained,(1)τit=Hiθitθ¨it+Ciθit,θ˙it+Giθitwhere Hi(θi(t))θ¨i(t),Gi(θi(t)), and Ci(θi(t),θ˙i(t)), respectively, represent inertia term, gravity term, centrifugal term, and coriolis force term in the joint drive torque. According to the motion characteristics of the palletizing robot, the model is considered to be simplified properly without affecting the accuracy of the dynamic modeling of the robot. Ignore the rotatory joint (i.e., the wrist) of the end effector and fix the end effector directly to the end of the small arm. Therefore, the dynamic model only describes the dynamic characteristics of the waist joints, shoulder joints, and elbow joints, and corresponding subscripti =1, 2, 3, respectively.


#### 2.2.2. Establishment of Joint Friction Model for Palletizing Robot

The Stribeck models of the waist joint, shoulder joint, and elbow joint of the MD1200-YJ palletizing robot are established and their parameters are identified.(1) Establishment of Joint Friction Torque Model of Robot. A single joint servo system mainly composed of AC servo motor and speed reducer is taken as an example, and its abstract model is shown in Figure 2. The speed reducer is divided into two parts according to the high-speed and low-speed stages. The motor output shaft is connected with the high-speed stage, which is called the high-speed shaft. The low speed stage is connected to the load, called the low-speed shaft.Figure 2
Single joint system abstract model.In Figure2,n is the deceleration ratio of the joint speed reducer; τmotor is the output torque of the motor (N·m); τconnect is output torque of high-speed level (N·m); n·τconnect is input torque of low-speed level (N·m);Jh,Jl are, respectively, the rotational inertia of the high-speed part and the low-speed part (kg·m2); θ¨h, θ¨l are, respectively, the angular acceleration of the high-speed part and the low-speed part (rad/s2); τfh,τfl are, respectively, the friction torque of the high-speed part and the friction torque of the low-speed part (N·m); τload is load torque (N·m).Here, Lagrange equation is used for dynamic modeling of low-speed shaft and high-speed shaft, respectively [17].Dynamic Model of High-Speed Shaft. As shown in Figure 2, the kinetic energy and potential energy of the high-speed part composed of motor and high-speed shaft are, respectively,(2)Kh=12Jhθ˙h2(3)Ph=a(4)Lh=Kh-Ph=12Jhθ˙h2-awhere,Kh,Ph, andLh are, respectively, the kinetic energy, the potential energy, and the Lagrange function of the high-speed part, andais a constant.Calculate, respectively,∂Lh/∂θ˙h, (d/dt)∂Lh/∂θ˙h and ∂Lh/∂θh, results are as follows:(5)∂Lh∂θ˙h=Jhθ˙h(6)ddt∂Lh∂θ˙h=Jhθ¨h(7)∂Lh∂θh=0The external force acting on this part has motor driving torqueτmotor, friction torque τfh, and restoring torque of the low-speed part τconnect, according to the Lagrange method and has (8)τmotor-τfh-τconnect=ddt∂Lh∂θ˙h-∂Lh∂θhSubstitute equations (5), (6), and (7) into(8), and the simplification is as follows:(9)τmotor=Jhθ¨h+τfh+τconnectStribeck model of friction torqueτfh is as follows [16]:(10)τfh=Mfch+Mfsh-Mfche-θ˙h/θ˙sh2sgn⁡θ˙h+σ2hθ˙hwhere Mfch,Mfsh,θ˙sh and σ2h are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient for high-speed part.Dynamic Model of Low-Speed Shaft. As shown in Figure 2, the kinetic energy and potential energy of the low-speed part composed of load and low-speed shaft of the speed reducer are, respectively,(11)Kl=12Jlθ˙l2(12)Pl=b(13)Ll=Kl-Pl=12Jlθ˙l2-bwhere,Kl,Pl andLl are, respectively, the kinetic energy, the potential energy, and the Lagrange function of the low-speed part, andbis a constant.Calculate, respectively,∂Ll/∂θ˙l, (d/dt)∂Ll/∂θ˙l, and ∂Ll/∂θl, and the results are as follows:(14)∂Ll∂θ˙l=Jlθ¨l(15)ddt∂Ll∂θ˙l=Jlθ¨l(16)∂Ll∂θl=0The external force acting on this part has restoring torque of the low-speed partn·τconnect, friction torque τfl, and load torque τload, according to the Lagrange method and has(17)nτconnect-τfl-τload=ddt∂Ll∂θ˙l-∂Ll∂θlSubstitute equations (14), (15), and (16) into(17), and the simplification is as follows:(18)nτconnect=Jlθ¨l+τfl+τloadStribeck model of friction torqueτfl is as follows:(19)τfl=Mfcl+Mfsl-Mfcle-θ˙l/θ˙sl2sgn⁡θ˙l+σ2lθ˙lwhere Mfcl,Mfsl,θ˙sl, and σ2l are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient for low-speed part.The expression for the load moment is as follows:(20)τload=Jloadθ¨+Twhere Jload is the rotational inertia of the load;T is torque generated by gravity, friction, and so no.The combination of equation (9) and equation (18) is obtained,(21)τmotor=Jhθ¨h+1nJlθ¨l+τfh+1nτfl+1nτloadwhere friction term τfh+(1/n)τfl is described by total friction τf, mean(22)τf=τfh+1nτfl(2) Simplified Method for Parameter Identification of Friction Model. In literature [28], the high-speed shaft and low-speed shaft of the reducer are separately used for friction modeling, which brings great difficulties to the identification of friction parameters. On the one hand, it is necessary to detect the friction in the gear reducer; this detection has difficulty. On the other hand, the number of friction parameters needed to be identified is large and the difficulty of identification is also increased.Literature [29] points out that surface topography, sliding speed, and contact pressure are important factors influencing friction behavior of sliding contact interface. Considering that material and lubrication conditions of high-speed shaft and low-speed shaft of the speed reducer are almost the same and friction behaviors are similar, in order to reduce the difficulty of parameter identification, both of them take the same static and dynamic friction parameters, mean(23)σ2h=σ2l;Mfch=Mfcl;Mfsh=Mfsl;θ˙sh=θ˙slThe friction parameters of four equations in equation (23) are expressed as uniform symbols, which are expressed as follows σ2,Mfc,Mfs, and θ˙S.Friction torque is a function of the angular velocity of the joint. Therefore, the corresponding frictional torque of the system at different steady state velocities can be measured, and frictional model can be fitted with the least square method, so as to obtain its identification parameters.Assuming that the observed value of the total friction torque isτf, predictive value is τ^f, definition identification error is(24)e=τf-τ^fDefinition the objective function for parameter identification:(25)Je=12∑i=1εei2Using the least square method, the friction parameters are obtained by fitting the measured speed-friction torque curve, and the minimum value of objective optimization functionJe can be obtained.(3) Establishment of Waist Joint Friction Model. Abstract model of waist joint system is shown in Figure 3. The speed reducer is composed of the first speed reduction part and the second speed reduction part. The first speed reduction part is composed of the input gear and the center gear, and its reduction ratio is n11=3.91. The second speed reduction part is input from spur gear to RV reducer. The RV reducer adopts the mode of shaft fixed and enclosure output, and its reduction ratio isn12=35.61. In addition,τmotor1 is output torque of waist motor (N·m); τconnect1 is output torque of input gear (N·m). τf1 is the total friction torque to the motor shaft, τf11 is the friction torque on the high-speed shaft, τf1Z is the total friction torque on the intermediate shaft, τf12 is the friction torque on the intermediate shaft, and τf13 is the friction torque on the second speed reduction part.Figure 3
Abstract model of the waist joint system.(1) Stribeck Friction Model of the Waist Joint System. Combined equations (10), (22), and (23) establish friction models for the first and second speed reduction parts, respectively.Friction model of the second speed reduction part:(26)τf1Z=τf12+1n12τf13=-Mfc1+Mfs1-Mfc1e-n12θ˙1/θ˙s12sgn⁡θ˙1-σ21n12θ˙1+Mfc1+Mfs1-Mfc1e-θ˙1/θ˙s12sgn⁡θ˙1+σ21θ˙1n12Friction model of the first speed reduction part, namely, the total frictional moment model of the waist joint:(27)τf1=τf11+1n11τf1Z=Mfc1+Mfs1-Mfc1e-n11n12θ˙1/θ˙s12sgn⁡θ˙1+σ21n11n12θ˙1+1n11τf1Zwhere θ˙1 is angular velocity of the waist joint; Mfs1, Mfc1, θ˙s1, and σ21 are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient, and these four parameters need to be identified.(2) Parameter Identification of Waist Joint Friction Model. The waist joint rotates at a constant velocity, the inertial force is zero, and the shaft of the joint is parallel to the direction of gravity; the motor output torque is equal to the system friction torquel therefore, the friction torque can be obtained by measuring the motor output torque, namely,(28)τmotor1=τf1Therefore, it is necessary to measure the motor output torque at different steady speed by experiment.The experiment mainly collects the motion parameters and output torque of the servo motor of the waist joint of MD1200-YJ palletizing robot. The main equipment of the experiment is shown in Table1. Figure 4 shows the experimental environment and relations of data stream between the experimental units. And the following experiments are generally conducted based on this process.Table 1
Experiment equipment.
Experimental equipment Type Manufacturer Industrial Personal Computer(IPC) BECKHOFF C6920 Beckhoff Automation GmbH Driver Cooldrive A8 Tsino Dynatron Data acquisition software DriveStarter V1.0 Tsino Dynatron Upper computer programming software TwinCAT 2 Beckhoff Automation GmbH PC ThinkPad E420 LenovoFigure 4
Experiment scene and relations of data stream between the experimental units.
(a)
Experiment scene (b)
Relations of data stream between the experimental unitsThe basic processes of robot control, data communication, and data acquisition are as follows:(a) The IPC, as the upper computer, writes the motion control commands of robot joint servo motors by programming software.(b) The upper computer transfers control commands to the servo driver, and the servo driver drives the motion of each joint servo motor according to the commands.(c) During the operation, the rotary encoder feedback data signals such as angular position and angular velocity to the servo driver.(d) Run data acquisition software on PC, establish communication, and collect the required physical data.In the experiment, only the waist joint moves; other joints do not move. The range of angular velocity of the waist joint is [0°/s, 100°/s]. In [0°/s, 1°/s], take a point every 0.2°/s, take a point every 1°/s in [1°/s, 5°/s], take a point every 10°/s in (5°/s, 80°/s), and measure the corresponding friction torque. The relationship between the two is shown in Figure5. Then, the friction parameters are identified by using the least square method, and the results are shown in Table 2.Table 2
Result of friction parameters of waist joint identification.
M f c 1 (N∙m) M f s 1 (N∙m) θ ˙ s 1 (rad∙s -1 ) σ 21 (N∙m∙rad -1 ) 0.827 1.276 15.271 0.056Figure 5
Relation between friction torque and angular velocity of waist joint.(4) Establishment of Shoulder and Elbow Joint Friction Models. The system abstract model of shoulder and elbow joints is shown in Figure 2. The reduction ratio of their reducer isn2=n3=171.(1) Establishment of Shoulder Joint Friction Model. Combining (10), (22) with (23), the shoulder joint friction model is established:(29)τf2=Mfc2+Mfs2-Mfc2e-n2θ˙2/θ˙s22sgnθ˙2+σ22n2θ˙2+Mfc2+Mfs2-Mfc2e-θ˙2/θ˙s22sgnθ˙2+σ22θ˙2n2where θ˙2 is angular velocity of the shoulder joint; Mfs2, Mfc2, θ˙s2, and σ22 need to be identified; τf2 is friction torque of the shoulder joint.(2) Parameter Identification of Shoulder Joint Friction Model. Unlike the waist joint, the shaft of the shoulder joint is perpendicular to the direction of gravity. Therefore, when shoulder joint moves at a constant velocity, it needs to overcome both the friction torque and the gravitational torque of the load. In order to measure the data of friction torque, the reciprocating constant velocity motion test method is adopted here [30]; namely, the following two groups of experiments are conducted.Experiment I. The shoulder joint constant speed moves fromθ2 toθ2’, and the output torqueu1 of the motor is measured.Experiment II. The shoulder joint constant speed moves fromθ2’ toθ2, and the output torqueu2 of the motor is measured.Assumption. The output torque of the motor in Experiment I is(30)u1=τgθ2+Ffθ˙2where τg(θ2) is gravitational torque; Ff(θ˙2) is friction torque.The output torque of the motor in Experiment II is(31)u2=τgθ2+Ff-θ˙2Because every joint of MD1200-YJ palletizing robot adopts high precision RV reducer; the friction difference in the positive and negative directions is very small. Therefore,(32)Ff-θ˙2=-Ffθ˙2From equations (30) ~ (32), it can be obtained that the friction torque and the weight torque are, respectively:(33)Ffθ˙2=12u1-u2(34)τgθ2=12u1+u2According to the experimental data, the friction torque at different angular velocities can be obtained from equation (33).In the experiment, only the shoulder joint moves, while the other joints do not move. Within the angular velocity range of the shoulder joint [0°/s, 75°/s], a point is taken every 0.2°/s within [0°/s, 1°/s], and a point is taken every 1°/s within [1°/s, 5°/s]. In the (5°/s, 75°/s), a point is taken every 10°/s, and the corresponding friction torque is measured. The relationship between the two is shown in Figure6. The friction parameters are then identified, and the results are shown in Table 3.Table 3
Result of static parameters identification.
M f c 2 (N∙m) M f s 2 (N∙m) θ ˙ s 2 (rad∙s -1 ) σ 22 (N∙m∙rad -1 )  0.565 3.031 26.103 0.057Figure 6
Relation between friction torque and angular velocity of shoulder joint.(3) Establishment of Elbow Joint Friction Model. The elbow joint friction model is as follows:(35)τf3=Mfc3+Mfs3-Mfc3e-n3θ˙3/θ˙s32sgn⁡θ˙3⁡+σ23n3θ˙3+Mfc3+Mfs3-Mfc3e-θ˙3/θ˙s32sgn⁡θ˙3+σ23θ˙3n3where θ˙3 is angular velocity of the shoulder joint; Mfs3, Mfc3, θ˙s3, and σ23 need to be identified; τf3 is friction torque of the shoulder joint.The shaft of elbow joint is also perpendicular to the direction of gravity. The measurement of friction torque and the recognition of friction parameters are the same as that of shoulder joint. Here, the experimental data is directly given, as shown in Figure7 and friction parameters, as shown in Table 4.Table 4
Result of static parameters identification.
M f c 3 (N∙m) M f s 3 (N∙m) θ ˙ s 3 (rad∙s -1 ) σ 23 (N∙m∙rad -1 )  0.015 0.259 16.703 0.0066Figure 7
Relation between friction torque and angular velocity of elbow joint.(5) The Dynamic Model of MD1200-YJ Palletizing Robot considering Joint Friction Is Established. Because the model of the friction torque derived in this paper is to convert the friction torque of the joint system to the motor shaft; therefore, it needs to be converted to the output end of the joint; there is(36)FSiθ˙it=τfi·niF S i ( θ ˙ i ( t ) ) is joint friction torque; ni is reduction ratio of each joint.On the basis of equation (1), the friction torque term FSi(θ˙i(t)) is added to obtain the dynamic model considering friction of jointi: (37)τi∗t=Hiθitθ¨it+Ciθit,θ˙it+Giθit+FSiθ˙itwhere τi∗ is joint driving torque considering joint friction.


#### 2.2.3. Establishment of Torque Model of Balance Spring Cylinder

The balancing spring cylinder of MD1200-YJ palletizing robot has inner and outer two springs, and the spring parameters are shown in Table5.Table 5
Main parameters of spring.
Spring Wire diameter d t (mm) Mean diameter D (mm) Pitch p (mm) Effective number n Free length L 0 (mm) Single ring stiffness K (N/mm) Inner spring 10 70 24.2 22 543 288 Outer spring 14 100 34.8 16 580 375The force of the spring cylinder is shown in Figure8. The cylinder barrel of the balance spring cylinder is connected with the upper big arm by hinge (at point B), and its cylinder rod is connected with the waist mounting bracket by hinge (at point A). The distance between the spring cylinder mounting hole on the waist mounting bracket and the shoulder joint shaft isd (230mm), and the angle between the big arm and the vertical direction isθ2. When the big arm moves from the OB’ vertical position to the OB position, the balance spring cylinder extends with the rotation of the big arm (that is, the increase ofθ2), while the internal spring is compressed, generating the balance spring force FT.Figure 8
Force analysis on spring balance mechanism.Spring stiffnessK:(38)KI=288nI(39)KO=375nOTwo springs are installed in parallel, the total spring stiffnessK: (40)K=KI+KOSpring preload:(41)F0I=KI·X0I(42)F0O=KO·X0OTwo springs are installed in parallel, the total preloadF0:(43)F0=F0I+F0OWorking length of the spring LAB:(44)LAB=d2+L22-2dL2cos⁡θ2Spring deformation ΔX:(45)ΔX=LAB-L2-dThen, spring tension force generated in the working processF is (46)F=F0+K·ΔXForce arm of spring tensionlf is(47)lf=L2dsin⁡θ2LABBalance torque of the springMT is(48)MTθ2=F·lf=F0+Kd2+L22-2dL2cos⁡θ2-L2+dL2dsin⁡θ2d2+L22-2dL2cos⁡θ2Therefore, the dynamic model of shoulder joint considering joint friction and spring moment is established.(49)τ~2t=H2θ2tθ¨2t+C2θ2t,θ˙2t+G2θ2t+FS2θ˙2t+MTθ2t


### 2.3. Establishment of Energy Consumption Model of MD1200-YJ Palletizing Robot

In this paper, the energy consumption of the robot joint is divided into the useful work to drive the load and dissipation energy due to friction and heat loss of the motor.Since this study is carried out in the no-load condition, compared with the other three joints, the load of the wrist joint is small and its movement speed is not high, so its energy consumption is less [18]. In order to simplify problem, energy consumption of the wrist joint is ignored. Therefore, the total energy consumption is only the sum of the energy consumption of the waist joint (i =1), shoulder joint (i =2), and elbow joint (i =3).The following energy consumption model can be established for the motion of jointi system int0 totf time:(50)WSi=∫t0tfτmotori·θ˙i·ni+ii2Ridtwhere WSi, τmotori, θ˙i, ni,ii, and Ri are, respectively, energy consumption, motor torque, angular velocity, reduction ratio of reducer, effective value of motor current, torque constant of motor, and equivalent resistance of motor of jointi.The motor torque model of the jointi can be obtained from equations (37) and (49)(51)τmotori=τiniAndτmotori has the following relationship with the effective currentii of the motor(52)τmotori=ki·iiwhereki is the torque constant of the motor of jointi.The thermal power model of joint motori can be obtained by joule’s law combined with equation (52)(53)PRi=τmotoriki2RiThus, the total energy consumption model of palletizing robot int0 totf time is established.(54)WS=∑i=13WSi=∑i=13∫t0tfτmotori·θ˙i·ni+τmotoriki2RidtFriction energy model:(55)Pfi=τfi·θ˙i·niThe total dissipative energy consumption model of palletizing robot int0 totf time is established:(56)Wloss=∑i=13∫t0tfPfi+PRidt=∑i=13∫t0tfτfi·θ˙i·ni+τmotoriki2RidtThe mechanical power of jointi(57)PEi=τmotori·θ˙i·niThe total mechanical energy model of palletizing robot int0 totf time is established:(58)WE=∑i=13∫t0tfPEidt=∑i=13∫t0tfτmotori·θ˙i·nidt


### 2.4. Validation and Analysis of Dynamic Model and Energy Consumption Model

Based on the example of 3-4-5 polynomial motion law, the multifactor dynamic model and energy consumption model are verified and analyzed by experiment and simulation.


#### 2.4.1. Description of Palletizing Robot Parameters and Calculation Examples

The main parameters of MD1200-YJ palletizing robot are shown in Tables6–8.Table 6
Main function size and material of palletizing robot.
Part name Main size parameter Material Big arm 1200 mm Zl110 Forearm 1200 mm Zl110 Drive arm of forearm 360 mm QT-5007 Driving link of forearm 1200 mm Q235A Horizontal holding link 1 1200 mm Q235A Horizontal holding link 2 1200 mm Q235ATable 7
Parameters of motors and reducer.
Joint Servo motor Reducer Waist joint Rated power: 5.5 kWTorque constant: 2.37 N·m/A Type: RV320CSpeed ratio: 139.1 Shoulder joint Rated power: 5.5 kWTorque constant: 2.37 N•m/A Type: RV160ESpeed ratio: 171 Elbow joint Rated power: 5.5 kWTorque constant: 2.37 N•m/A Type: RV160ESpeed ratio: 171Table 8
Mass of the main part of palletizing robot.
Part name Mass (kg) Part name Mass (kg) Waist joint 486.085 Big arm 107.105 Forearm 67.202 Drive arm of forearm 13.513 End effector 31.10716 Driving link of forearm 31.107Problem Description of Calculation Example. In the basic coordinate system, the end position and posture of the palletizing robot are represented by coordinates (X,Y,Z,Φ). The torque and mechanical power of the joint motor at the end of the joint from the base coordinateS0 (1200mm, 1800mm, 800mm, 0rad) to the pose coordinateSf (-1200mm, 1200mm, 500mm, 0rad) are discussed. The interval of starting and ending time of motion is [0s, 3.5s]. Through the forward and inverse kinematics solution, the joint coordinates corresponding toS0 andSf, and the motion parameters at the start and stop time, as shown in Tables 9 and 10, the whole process of no-load operation is conducted.Table 9
Joint position ofS0.
Joint Angular position(rad) Angular velocity(rad/s) Angular acceleration(rad/s 2 ) θ 1(t0) 0.9828 0 0 θ 2(t0) 1.10251 0 0 θ 3(t0) 2.6202 0 0Table 10
Joint position ofSf.
Joint Angular position(rad) Angular velocity(rad/s) Angular acceleration(rad/s 2 ) θ 1(tf) 2.3562 0 0 θ 2(tf) 1.1519 0 0 θ 3(tf) 2.2028 0 0According to the conditions given in this example, the trajectory planning of 3-4-5 polynomial motion law commonly used in engineering is carried out, and the expression is as follows:(59)θ1=0.9828+0.3203t3-0.1373t4+0.0157t5θ2=1.0251+0.0296t3-0.0127t4+0.0014t5θ2′=2.6202-0.0974t3+0.0417t4-0.0048t5t∈03.5The curves of end trajectory and the angular position of waist joint, shoulder joint, and elbow joint are shown in Figures9 and 10.Figure 9
End trajectory of 3-4-5 polynomial.Figure 10
Angular displacement trajectory of 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow joint


#### 2.4.2. Validation of Dynamic Model and Energy Consumption Model

According to the 3-4-5 polynomial motion law and the known motion parameters, the motion of each joint of the palletizing robot is controlled, and the data acquisition software is used to collect the torque and mechanical power data of each joint motor within a motion period.It should be noted that since the data acquisition software cannot directly measure the heat loss of the motor and friction energy consumption, therefore, the total energy consumption of the palletizing robot cannot be directly measured. However, according to equation (54), the accuracy of the energy consumption model can be indirectly verified by comparing the experimental value of mechanical power with the simulated value.The experimental data collected show a trend of high frequency oscillation, but the overall trend is regular. Therefore, the fitting curve of the experimental data is obtained by the polynomial curve fitting method, as shown in Figure11, and the fitting curve of experimental data is compared with the simulation curve of the theoretical model, as shown in Figure 12.Figure 11
Experimental fitting curves of joint motor torque.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 12
Comparison of joint motor torque experimental curve with simulation curve.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointIt can be seen from Figure12 that the variation trend of the motor torque simulation curve of the three joints is similar to that of the experimental curve, which indicates that the multifactor dynamics model has higher accuracy. However, it can be seen that there is a certain error between the two, which may be caused by (1) the scale parameters and inertia parameters of all parts of the simulation model are measured by the 3D model, which are different from the corresponding parts of the actual prototype; (2) there will be errors in the experimental curve fitting.The influence of friction torque and spring torque on output torque of shoulder joint motor is analyzed, as shown in Figure13. It shows that the friction torque and spring torque have obvious effects on the output torque of the motor. The conventional rigid body dynamics model cannot fully reflect components of the joint torque. The addition of friction torque and spring torque makes the output torque of the shoulder joint closer to the real situation.Figure 13
Comparison of shoulder joint motor torque experimental curve with simulation curve under different conditions.
(a)
No spring and no friction (b)
No spring and consider friction (c)
Consider spring and no friction (d)
Consider spring and frictionNext, the mechanical power data of each joint are analyzed. The mechanical power experimental data of waist joint, shoulder joint, and elbow joint are obtained through data acquisition software. The phenomenon of high frequency oscillation still exists, and the fitting curve of the experimental data is obtained by the polynomial curve fitting method, as shown in Figure14, and compared with the simulation curve of mechanical power, as shown in Figure 15.Figure 14
Experimental fitting curves of joint mechanical power.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 15
Comparison of joint mechanical power experimental curve with simulation curve.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointAs can be seen from Figure15, the variation trend of the mechanical power simulation curve of the three joints is similar to that of the experimental curve. As can be seen from Table 11, the theoretical value of mechanical power consumption differs from the experimental value by 80.2988 J, with an error of about 7.69%, which proves the accuracy of the energy consumption model.Table 11
Mechanical power of joints (J).
Comparison object Waist joint Shoulder joint Elbow joint Total mechanical power consumption Simulation 837.0672 103.0638 104.6278 1044.8 Experimental 776.6805 82.0542 105.7665 964.5012 Error (%) 7.21408 20.38504 1.088334 7.685567


## 2.1. Introduction to the Mechanism of MD1200-YJ Palletizing Robot

The MD1200-YJ palletizing robot is a 4-DOF multijoint robot, load capacity: 120kg, maximum turning radius: 2400mm, repeat positioning accuracy: ±0.4mm, waist maximum rotation speed: 85°/s, belonging to high speed and heavy load palletizing robot, its 3D model shown in Figure1. It is mainly composed of machine base, waist joint, small arm joint, drive arm of small arm, balance spring cylinder, balancing weight, small arm drive link, attitude hold link1, attitude hold tripod, attitude hold link2, small arm link, big arm joint, big arm link, waist mounting bracket, end effector mounting bracket, and so on. The action of two parallelogram mechanical linkages (one consists of attitude hold link1, big arm link, attitude hold tripod, and waist mounting bracket, and the other consists of attitude hold tripod, attitude hold link2, small arm link, and end effector mounting bracket) makes undersurface of end effector keep parallel to the ground.Figure 1
3D model of MD1200-YJ palletizing robot.


## 2.2. Establishment of Multifactor Dynamic Model of MD1200-YJ Palletizing Robot

On the basis of the conventional rigid body inverse dynamics model of robot, the multifactor inverse dynamics model of palletizing robot with considering joint friction torque and spring balance torque is established.


### 2.2.1. Rigid Body Dynamics Model

Rigid body dynamics model for MD1200-YJ palletizing robot is derived using Kane equation method in this paper. The basic idea of Kane equation method is to replace the generalized coordinate with the generalized velocity as the independent variable of the system. When using Kane method for dynamic modeling, it only involves the cross product operation of the vector’s dot product, which greatly improves the calculation efficiency. Meanwhile, the speed and acceleration of connecting link parts solution have the recursive nature, which is convenient for the computer to carry out numerical calculation.Literature [12, 13] described derivation ideas of Kane equation method in detail, and here we will not repeat them, but directly give the derivation results.Assuming thatτi represents the drive torque of jointi, and the rigid body inverse dynamics equation of the robot system is obtained,(1)τit=Hiθitθ¨it+Ciθit,θ˙it+Giθitwhere Hi(θi(t))θ¨i(t),Gi(θi(t)), and Ci(θi(t),θ˙i(t)), respectively, represent inertia term, gravity term, centrifugal term, and coriolis force term in the joint drive torque. According to the motion characteristics of the palletizing robot, the model is considered to be simplified properly without affecting the accuracy of the dynamic modeling of the robot. Ignore the rotatory joint (i.e., the wrist) of the end effector and fix the end effector directly to the end of the small arm. Therefore, the dynamic model only describes the dynamic characteristics of the waist joints, shoulder joints, and elbow joints, and corresponding subscripti =1, 2, 3, respectively.


### 2.2.2. Establishment of Joint Friction Model for Palletizing Robot

The Stribeck models of the waist joint, shoulder joint, and elbow joint of the MD1200-YJ palletizing robot are established and their parameters are identified.(1) Establishment of Joint Friction Torque Model of Robot. A single joint servo system mainly composed of AC servo motor and speed reducer is taken as an example, and its abstract model is shown in Figure 2. The speed reducer is divided into two parts according to the high-speed and low-speed stages. The motor output shaft is connected with the high-speed stage, which is called the high-speed shaft. The low speed stage is connected to the load, called the low-speed shaft.Figure 2
Single joint system abstract model.In Figure2,n is the deceleration ratio of the joint speed reducer; τmotor is the output torque of the motor (N·m); τconnect is output torque of high-speed level (N·m); n·τconnect is input torque of low-speed level (N·m);Jh,Jl are, respectively, the rotational inertia of the high-speed part and the low-speed part (kg·m2); θ¨h, θ¨l are, respectively, the angular acceleration of the high-speed part and the low-speed part (rad/s2); τfh,τfl are, respectively, the friction torque of the high-speed part and the friction torque of the low-speed part (N·m); τload is load torque (N·m).Here, Lagrange equation is used for dynamic modeling of low-speed shaft and high-speed shaft, respectively [17].Dynamic Model of High-Speed Shaft. As shown in Figure 2, the kinetic energy and potential energy of the high-speed part composed of motor and high-speed shaft are, respectively,(2)Kh=12Jhθ˙h2(3)Ph=a(4)Lh=Kh-Ph=12Jhθ˙h2-awhere,Kh,Ph, andLh are, respectively, the kinetic energy, the potential energy, and the Lagrange function of the high-speed part, andais a constant.Calculate, respectively,∂Lh/∂θ˙h, (d/dt)∂Lh/∂θ˙h and ∂Lh/∂θh, results are as follows:(5)∂Lh∂θ˙h=Jhθ˙h(6)ddt∂Lh∂θ˙h=Jhθ¨h(7)∂Lh∂θh=0The external force acting on this part has motor driving torqueτmotor, friction torque τfh, and restoring torque of the low-speed part τconnect, according to the Lagrange method and has (8)τmotor-τfh-τconnect=ddt∂Lh∂θ˙h-∂Lh∂θhSubstitute equations (5), (6), and (7) into(8), and the simplification is as follows:(9)τmotor=Jhθ¨h+τfh+τconnectStribeck model of friction torqueτfh is as follows [16]:(10)τfh=Mfch+Mfsh-Mfche-θ˙h/θ˙sh2sgn⁡θ˙h+σ2hθ˙hwhere Mfch,Mfsh,θ˙sh and σ2h are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient for high-speed part.Dynamic Model of Low-Speed Shaft. As shown in Figure 2, the kinetic energy and potential energy of the low-speed part composed of load and low-speed shaft of the speed reducer are, respectively,(11)Kl=12Jlθ˙l2(12)Pl=b(13)Ll=Kl-Pl=12Jlθ˙l2-bwhere,Kl,Pl andLl are, respectively, the kinetic energy, the potential energy, and the Lagrange function of the low-speed part, andbis a constant.Calculate, respectively,∂Ll/∂θ˙l, (d/dt)∂Ll/∂θ˙l, and ∂Ll/∂θl, and the results are as follows:(14)∂Ll∂θ˙l=Jlθ¨l(15)ddt∂Ll∂θ˙l=Jlθ¨l(16)∂Ll∂θl=0The external force acting on this part has restoring torque of the low-speed partn·τconnect, friction torque τfl, and load torque τload, according to the Lagrange method and has(17)nτconnect-τfl-τload=ddt∂Ll∂θ˙l-∂Ll∂θlSubstitute equations (14), (15), and (16) into(17), and the simplification is as follows:(18)nτconnect=Jlθ¨l+τfl+τloadStribeck model of friction torqueτfl is as follows:(19)τfl=Mfcl+Mfsl-Mfcle-θ˙l/θ˙sl2sgn⁡θ˙l+σ2lθ˙lwhere Mfcl,Mfsl,θ˙sl, and σ2l are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient for low-speed part.The expression for the load moment is as follows:(20)τload=Jloadθ¨+Twhere Jload is the rotational inertia of the load;T is torque generated by gravity, friction, and so no.The combination of equation (9) and equation (18) is obtained,(21)τmotor=Jhθ¨h+1nJlθ¨l+τfh+1nτfl+1nτloadwhere friction term τfh+(1/n)τfl is described by total friction τf, mean(22)τf=τfh+1nτfl(2) Simplified Method for Parameter Identification of Friction Model. In literature [28], the high-speed shaft and low-speed shaft of the reducer are separately used for friction modeling, which brings great difficulties to the identification of friction parameters. On the one hand, it is necessary to detect the friction in the gear reducer; this detection has difficulty. On the other hand, the number of friction parameters needed to be identified is large and the difficulty of identification is also increased.Literature [29] points out that surface topography, sliding speed, and contact pressure are important factors influencing friction behavior of sliding contact interface. Considering that material and lubrication conditions of high-speed shaft and low-speed shaft of the speed reducer are almost the same and friction behaviors are similar, in order to reduce the difficulty of parameter identification, both of them take the same static and dynamic friction parameters, mean(23)σ2h=σ2l;Mfch=Mfcl;Mfsh=Mfsl;θ˙sh=θ˙slThe friction parameters of four equations in equation (23) are expressed as uniform symbols, which are expressed as follows σ2,Mfc,Mfs, and θ˙S.Friction torque is a function of the angular velocity of the joint. Therefore, the corresponding frictional torque of the system at different steady state velocities can be measured, and frictional model can be fitted with the least square method, so as to obtain its identification parameters.Assuming that the observed value of the total friction torque isτf, predictive value is τ^f, definition identification error is(24)e=τf-τ^fDefinition the objective function for parameter identification:(25)Je=12∑i=1εei2Using the least square method, the friction parameters are obtained by fitting the measured speed-friction torque curve, and the minimum value of objective optimization functionJe can be obtained.(3) Establishment of Waist Joint Friction Model. Abstract model of waist joint system is shown in Figure 3. The speed reducer is composed of the first speed reduction part and the second speed reduction part. The first speed reduction part is composed of the input gear and the center gear, and its reduction ratio is n11=3.91. The second speed reduction part is input from spur gear to RV reducer. The RV reducer adopts the mode of shaft fixed and enclosure output, and its reduction ratio isn12=35.61. In addition,τmotor1 is output torque of waist motor (N·m); τconnect1 is output torque of input gear (N·m). τf1 is the total friction torque to the motor shaft, τf11 is the friction torque on the high-speed shaft, τf1Z is the total friction torque on the intermediate shaft, τf12 is the friction torque on the intermediate shaft, and τf13 is the friction torque on the second speed reduction part.Figure 3
Abstract model of the waist joint system.(1) Stribeck Friction Model of the Waist Joint System. Combined equations (10), (22), and (23) establish friction models for the first and second speed reduction parts, respectively.Friction model of the second speed reduction part:(26)τf1Z=τf12+1n12τf13=-Mfc1+Mfs1-Mfc1e-n12θ˙1/θ˙s12sgn⁡θ˙1-σ21n12θ˙1+Mfc1+Mfs1-Mfc1e-θ˙1/θ˙s12sgn⁡θ˙1+σ21θ˙1n12Friction model of the first speed reduction part, namely, the total frictional moment model of the waist joint:(27)τf1=τf11+1n11τf1Z=Mfc1+Mfs1-Mfc1e-n11n12θ˙1/θ˙s12sgn⁡θ˙1+σ21n11n12θ˙1+1n11τf1Zwhere θ˙1 is angular velocity of the waist joint; Mfs1, Mfc1, θ˙s1, and σ21 are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient, and these four parameters need to be identified.(2) Parameter Identification of Waist Joint Friction Model. The waist joint rotates at a constant velocity, the inertial force is zero, and the shaft of the joint is parallel to the direction of gravity; the motor output torque is equal to the system friction torquel therefore, the friction torque can be obtained by measuring the motor output torque, namely,(28)τmotor1=τf1Therefore, it is necessary to measure the motor output torque at different steady speed by experiment.The experiment mainly collects the motion parameters and output torque of the servo motor of the waist joint of MD1200-YJ palletizing robot. The main equipment of the experiment is shown in Table1. Figure 4 shows the experimental environment and relations of data stream between the experimental units. And the following experiments are generally conducted based on this process.Table 1
Experiment equipment.
Experimental equipment Type Manufacturer Industrial Personal Computer(IPC) BECKHOFF C6920 Beckhoff Automation GmbH Driver Cooldrive A8 Tsino Dynatron Data acquisition software DriveStarter V1.0 Tsino Dynatron Upper computer programming software TwinCAT 2 Beckhoff Automation GmbH PC ThinkPad E420 LenovoFigure 4
Experiment scene and relations of data stream between the experimental units.
(a)
Experiment scene (b)
Relations of data stream between the experimental unitsThe basic processes of robot control, data communication, and data acquisition are as follows:(a) The IPC, as the upper computer, writes the motion control commands of robot joint servo motors by programming software.(b) The upper computer transfers control commands to the servo driver, and the servo driver drives the motion of each joint servo motor according to the commands.(c) During the operation, the rotary encoder feedback data signals such as angular position and angular velocity to the servo driver.(d) Run data acquisition software on PC, establish communication, and collect the required physical data.In the experiment, only the waist joint moves; other joints do not move. The range of angular velocity of the waist joint is [0°/s, 100°/s]. In [0°/s, 1°/s], take a point every 0.2°/s, take a point every 1°/s in [1°/s, 5°/s], take a point every 10°/s in (5°/s, 80°/s), and measure the corresponding friction torque. The relationship between the two is shown in Figure5. Then, the friction parameters are identified by using the least square method, and the results are shown in Table 2.Table 2
Result of friction parameters of waist joint identification.
M f c 1 (N∙m) M f s 1 (N∙m) θ ˙ s 1 (rad∙s -1 ) σ 21 (N∙m∙rad -1 ) 0.827 1.276 15.271 0.056Figure 5
Relation between friction torque and angular velocity of waist joint.(4) Establishment of Shoulder and Elbow Joint Friction Models. The system abstract model of shoulder and elbow joints is shown in Figure 2. The reduction ratio of their reducer isn2=n3=171.(1) Establishment of Shoulder Joint Friction Model. Combining (10), (22) with (23), the shoulder joint friction model is established:(29)τf2=Mfc2+Mfs2-Mfc2e-n2θ˙2/θ˙s22sgnθ˙2+σ22n2θ˙2+Mfc2+Mfs2-Mfc2e-θ˙2/θ˙s22sgnθ˙2+σ22θ˙2n2where θ˙2 is angular velocity of the shoulder joint; Mfs2, Mfc2, θ˙s2, and σ22 need to be identified; τf2 is friction torque of the shoulder joint.(2) Parameter Identification of Shoulder Joint Friction Model. Unlike the waist joint, the shaft of the shoulder joint is perpendicular to the direction of gravity. Therefore, when shoulder joint moves at a constant velocity, it needs to overcome both the friction torque and the gravitational torque of the load. In order to measure the data of friction torque, the reciprocating constant velocity motion test method is adopted here [30]; namely, the following two groups of experiments are conducted.Experiment I. The shoulder joint constant speed moves fromθ2 toθ2’, and the output torqueu1 of the motor is measured.Experiment II. The shoulder joint constant speed moves fromθ2’ toθ2, and the output torqueu2 of the motor is measured.Assumption. The output torque of the motor in Experiment I is(30)u1=τgθ2+Ffθ˙2where τg(θ2) is gravitational torque; Ff(θ˙2) is friction torque.The output torque of the motor in Experiment II is(31)u2=τgθ2+Ff-θ˙2Because every joint of MD1200-YJ palletizing robot adopts high precision RV reducer; the friction difference in the positive and negative directions is very small. Therefore,(32)Ff-θ˙2=-Ffθ˙2From equations (30) ~ (32), it can be obtained that the friction torque and the weight torque are, respectively:(33)Ffθ˙2=12u1-u2(34)τgθ2=12u1+u2According to the experimental data, the friction torque at different angular velocities can be obtained from equation (33).In the experiment, only the shoulder joint moves, while the other joints do not move. Within the angular velocity range of the shoulder joint [0°/s, 75°/s], a point is taken every 0.2°/s within [0°/s, 1°/s], and a point is taken every 1°/s within [1°/s, 5°/s]. In the (5°/s, 75°/s), a point is taken every 10°/s, and the corresponding friction torque is measured. The relationship between the two is shown in Figure6. The friction parameters are then identified, and the results are shown in Table 3.Table 3
Result of static parameters identification.
M f c 2 (N∙m) M f s 2 (N∙m) θ ˙ s 2 (rad∙s -1 ) σ 22 (N∙m∙rad -1 )  0.565 3.031 26.103 0.057Figure 6
Relation between friction torque and angular velocity of shoulder joint.(3) Establishment of Elbow Joint Friction Model. The elbow joint friction model is as follows:(35)τf3=Mfc3+Mfs3-Mfc3e-n3θ˙3/θ˙s32sgn⁡θ˙3⁡+σ23n3θ˙3+Mfc3+Mfs3-Mfc3e-θ˙3/θ˙s32sgn⁡θ˙3+σ23θ˙3n3where θ˙3 is angular velocity of the shoulder joint; Mfs3, Mfc3, θ˙s3, and σ23 need to be identified; τf3 is friction torque of the shoulder joint.The shaft of elbow joint is also perpendicular to the direction of gravity. The measurement of friction torque and the recognition of friction parameters are the same as that of shoulder joint. Here, the experimental data is directly given, as shown in Figure7 and friction parameters, as shown in Table 4.Table 4
Result of static parameters identification.
M f c 3 (N∙m) M f s 3 (N∙m) θ ˙ s 3 (rad∙s -1 ) σ 23 (N∙m∙rad -1 )  0.015 0.259 16.703 0.0066Figure 7
Relation between friction torque and angular velocity of elbow joint.(5) The Dynamic Model of MD1200-YJ Palletizing Robot considering Joint Friction Is Established. Because the model of the friction torque derived in this paper is to convert the friction torque of the joint system to the motor shaft; therefore, it needs to be converted to the output end of the joint; there is(36)FSiθ˙it=τfi·niF S i ( θ ˙ i ( t ) ) is joint friction torque; ni is reduction ratio of each joint.On the basis of equation (1), the friction torque term FSi(θ˙i(t)) is added to obtain the dynamic model considering friction of jointi: (37)τi∗t=Hiθitθ¨it+Ciθit,θ˙it+Giθit+FSiθ˙itwhere τi∗ is joint driving torque considering joint friction.


### 2.2.3. Establishment of Torque Model of Balance Spring Cylinder

The balancing spring cylinder of MD1200-YJ palletizing robot has inner and outer two springs, and the spring parameters are shown in Table5.Table 5
Main parameters of spring.
Spring Wire diameter d t (mm) Mean diameter D (mm) Pitch p (mm) Effective number n Free length L 0 (mm) Single ring stiffness K (N/mm) Inner spring 10 70 24.2 22 543 288 Outer spring 14 100 34.8 16 580 375The force of the spring cylinder is shown in Figure8. The cylinder barrel of the balance spring cylinder is connected with the upper big arm by hinge (at point B), and its cylinder rod is connected with the waist mounting bracket by hinge (at point A). The distance between the spring cylinder mounting hole on the waist mounting bracket and the shoulder joint shaft isd (230mm), and the angle between the big arm and the vertical direction isθ2. When the big arm moves from the OB’ vertical position to the OB position, the balance spring cylinder extends with the rotation of the big arm (that is, the increase ofθ2), while the internal spring is compressed, generating the balance spring force FT.Figure 8
Force analysis on spring balance mechanism.Spring stiffnessK:(38)KI=288nI(39)KO=375nOTwo springs are installed in parallel, the total spring stiffnessK: (40)K=KI+KOSpring preload:(41)F0I=KI·X0I(42)F0O=KO·X0OTwo springs are installed in parallel, the total preloadF0:(43)F0=F0I+F0OWorking length of the spring LAB:(44)LAB=d2+L22-2dL2cos⁡θ2Spring deformation ΔX:(45)ΔX=LAB-L2-dThen, spring tension force generated in the working processF is (46)F=F0+K·ΔXForce arm of spring tensionlf is(47)lf=L2dsin⁡θ2LABBalance torque of the springMT is(48)MTθ2=F·lf=F0+Kd2+L22-2dL2cos⁡θ2-L2+dL2dsin⁡θ2d2+L22-2dL2cos⁡θ2Therefore, the dynamic model of shoulder joint considering joint friction and spring moment is established.(49)τ~2t=H2θ2tθ¨2t+C2θ2t,θ˙2t+G2θ2t+FS2θ˙2t+MTθ2t


## 2.2.1. Rigid Body Dynamics Model

Rigid body dynamics model for MD1200-YJ palletizing robot is derived using Kane equation method in this paper. The basic idea of Kane equation method is to replace the generalized coordinate with the generalized velocity as the independent variable of the system. When using Kane method for dynamic modeling, it only involves the cross product operation of the vector’s dot product, which greatly improves the calculation efficiency. Meanwhile, the speed and acceleration of connecting link parts solution have the recursive nature, which is convenient for the computer to carry out numerical calculation.Literature [12, 13] described derivation ideas of Kane equation method in detail, and here we will not repeat them, but directly give the derivation results.Assuming thatτi represents the drive torque of jointi, and the rigid body inverse dynamics equation of the robot system is obtained,(1)τit=Hiθitθ¨it+Ciθit,θ˙it+Giθitwhere Hi(θi(t))θ¨i(t),Gi(θi(t)), and Ci(θi(t),θ˙i(t)), respectively, represent inertia term, gravity term, centrifugal term, and coriolis force term in the joint drive torque. According to the motion characteristics of the palletizing robot, the model is considered to be simplified properly without affecting the accuracy of the dynamic modeling of the robot. Ignore the rotatory joint (i.e., the wrist) of the end effector and fix the end effector directly to the end of the small arm. Therefore, the dynamic model only describes the dynamic characteristics of the waist joints, shoulder joints, and elbow joints, and corresponding subscripti =1, 2, 3, respectively.


## 2.2.2. Establishment of Joint Friction Model for Palletizing Robot

The Stribeck models of the waist joint, shoulder joint, and elbow joint of the MD1200-YJ palletizing robot are established and their parameters are identified.(1) Establishment of Joint Friction Torque Model of Robot. A single joint servo system mainly composed of AC servo motor and speed reducer is taken as an example, and its abstract model is shown in Figure 2. The speed reducer is divided into two parts according to the high-speed and low-speed stages. The motor output shaft is connected with the high-speed stage, which is called the high-speed shaft. The low speed stage is connected to the load, called the low-speed shaft.Figure 2
Single joint system abstract model.In Figure2,n is the deceleration ratio of the joint speed reducer; τmotor is the output torque of the motor (N·m); τconnect is output torque of high-speed level (N·m); n·τconnect is input torque of low-speed level (N·m);Jh,Jl are, respectively, the rotational inertia of the high-speed part and the low-speed part (kg·m2); θ¨h, θ¨l are, respectively, the angular acceleration of the high-speed part and the low-speed part (rad/s2); τfh,τfl are, respectively, the friction torque of the high-speed part and the friction torque of the low-speed part (N·m); τload is load torque (N·m).Here, Lagrange equation is used for dynamic modeling of low-speed shaft and high-speed shaft, respectively [17].Dynamic Model of High-Speed Shaft. As shown in Figure 2, the kinetic energy and potential energy of the high-speed part composed of motor and high-speed shaft are, respectively,(2)Kh=12Jhθ˙h2(3)Ph=a(4)Lh=Kh-Ph=12Jhθ˙h2-awhere,Kh,Ph, andLh are, respectively, the kinetic energy, the potential energy, and the Lagrange function of the high-speed part, andais a constant.Calculate, respectively,∂Lh/∂θ˙h, (d/dt)∂Lh/∂θ˙h and ∂Lh/∂θh, results are as follows:(5)∂Lh∂θ˙h=Jhθ˙h(6)ddt∂Lh∂θ˙h=Jhθ¨h(7)∂Lh∂θh=0The external force acting on this part has motor driving torqueτmotor, friction torque τfh, and restoring torque of the low-speed part τconnect, according to the Lagrange method and has (8)τmotor-τfh-τconnect=ddt∂Lh∂θ˙h-∂Lh∂θhSubstitute equations (5), (6), and (7) into(8), and the simplification is as follows:(9)τmotor=Jhθ¨h+τfh+τconnectStribeck model of friction torqueτfh is as follows [16]:(10)τfh=Mfch+Mfsh-Mfche-θ˙h/θ˙sh2sgn⁡θ˙h+σ2hθ˙hwhere Mfch,Mfsh,θ˙sh and σ2h are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient for high-speed part.Dynamic Model of Low-Speed Shaft. As shown in Figure 2, the kinetic energy and potential energy of the low-speed part composed of load and low-speed shaft of the speed reducer are, respectively,(11)Kl=12Jlθ˙l2(12)Pl=b(13)Ll=Kl-Pl=12Jlθ˙l2-bwhere,Kl,Pl andLl are, respectively, the kinetic energy, the potential energy, and the Lagrange function of the low-speed part, andbis a constant.Calculate, respectively,∂Ll/∂θ˙l, (d/dt)∂Ll/∂θ˙l, and ∂Ll/∂θl, and the results are as follows:(14)∂Ll∂θ˙l=Jlθ¨l(15)ddt∂Ll∂θ˙l=Jlθ¨l(16)∂Ll∂θl=0The external force acting on this part has restoring torque of the low-speed partn·τconnect, friction torque τfl, and load torque τload, according to the Lagrange method and has(17)nτconnect-τfl-τload=ddt∂Ll∂θ˙l-∂Ll∂θlSubstitute equations (14), (15), and (16) into(17), and the simplification is as follows:(18)nτconnect=Jlθ¨l+τfl+τloadStribeck model of friction torqueτfl is as follows:(19)τfl=Mfcl+Mfsl-Mfcle-θ˙l/θ˙sl2sgn⁡θ˙l+σ2lθ˙lwhere Mfcl,Mfsl,θ˙sl, and σ2l are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient for low-speed part.The expression for the load moment is as follows:(20)τload=Jloadθ¨+Twhere Jload is the rotational inertia of the load;T is torque generated by gravity, friction, and so no.The combination of equation (9) and equation (18) is obtained,(21)τmotor=Jhθ¨h+1nJlθ¨l+τfh+1nτfl+1nτloadwhere friction term τfh+(1/n)τfl is described by total friction τf, mean(22)τf=τfh+1nτfl(2) Simplified Method for Parameter Identification of Friction Model. In literature [28], the high-speed shaft and low-speed shaft of the reducer are separately used for friction modeling, which brings great difficulties to the identification of friction parameters. On the one hand, it is necessary to detect the friction in the gear reducer; this detection has difficulty. On the other hand, the number of friction parameters needed to be identified is large and the difficulty of identification is also increased.Literature [29] points out that surface topography, sliding speed, and contact pressure are important factors influencing friction behavior of sliding contact interface. Considering that material and lubrication conditions of high-speed shaft and low-speed shaft of the speed reducer are almost the same and friction behaviors are similar, in order to reduce the difficulty of parameter identification, both of them take the same static and dynamic friction parameters, mean(23)σ2h=σ2l;Mfch=Mfcl;Mfsh=Mfsl;θ˙sh=θ˙slThe friction parameters of four equations in equation (23) are expressed as uniform symbols, which are expressed as follows σ2,Mfc,Mfs, and θ˙S.Friction torque is a function of the angular velocity of the joint. Therefore, the corresponding frictional torque of the system at different steady state velocities can be measured, and frictional model can be fitted with the least square method, so as to obtain its identification parameters.Assuming that the observed value of the total friction torque isτf, predictive value is τ^f, definition identification error is(24)e=τf-τ^fDefinition the objective function for parameter identification:(25)Je=12∑i=1εei2Using the least square method, the friction parameters are obtained by fitting the measured speed-friction torque curve, and the minimum value of objective optimization functionJe can be obtained.(3) Establishment of Waist Joint Friction Model. Abstract model of waist joint system is shown in Figure 3. The speed reducer is composed of the first speed reduction part and the second speed reduction part. The first speed reduction part is composed of the input gear and the center gear, and its reduction ratio is n11=3.91. The second speed reduction part is input from spur gear to RV reducer. The RV reducer adopts the mode of shaft fixed and enclosure output, and its reduction ratio isn12=35.61. In addition,τmotor1 is output torque of waist motor (N·m); τconnect1 is output torque of input gear (N·m). τf1 is the total friction torque to the motor shaft, τf11 is the friction torque on the high-speed shaft, τf1Z is the total friction torque on the intermediate shaft, τf12 is the friction torque on the intermediate shaft, and τf13 is the friction torque on the second speed reduction part.Figure 3
Abstract model of the waist joint system.(1) Stribeck Friction Model of the Waist Joint System. Combined equations (10), (22), and (23) establish friction models for the first and second speed reduction parts, respectively.Friction model of the second speed reduction part:(26)τf1Z=τf12+1n12τf13=-Mfc1+Mfs1-Mfc1e-n12θ˙1/θ˙s12sgn⁡θ˙1-σ21n12θ˙1+Mfc1+Mfs1-Mfc1e-θ˙1/θ˙s12sgn⁡θ˙1+σ21θ˙1n12Friction model of the first speed reduction part, namely, the total frictional moment model of the waist joint:(27)τf1=τf11+1n11τf1Z=Mfc1+Mfs1-Mfc1e-n11n12θ˙1/θ˙s12sgn⁡θ˙1+σ21n11n12θ˙1+1n11τf1Zwhere θ˙1 is angular velocity of the waist joint; Mfs1, Mfc1, θ˙s1, and σ21 are, respectively, maximum static friction force, Coulomb friction force, Stribeck characteristic velocity, and viscosity friction coefficient, and these four parameters need to be identified.(2) Parameter Identification of Waist Joint Friction Model. The waist joint rotates at a constant velocity, the inertial force is zero, and the shaft of the joint is parallel to the direction of gravity; the motor output torque is equal to the system friction torquel therefore, the friction torque can be obtained by measuring the motor output torque, namely,(28)τmotor1=τf1Therefore, it is necessary to measure the motor output torque at different steady speed by experiment.The experiment mainly collects the motion parameters and output torque of the servo motor of the waist joint of MD1200-YJ palletizing robot. The main equipment of the experiment is shown in Table1. Figure 4 shows the experimental environment and relations of data stream between the experimental units. And the following experiments are generally conducted based on this process.Table 1
Experiment equipment.
Experimental equipment Type Manufacturer Industrial Personal Computer(IPC) BECKHOFF C6920 Beckhoff Automation GmbH Driver Cooldrive A8 Tsino Dynatron Data acquisition software DriveStarter V1.0 Tsino Dynatron Upper computer programming software TwinCAT 2 Beckhoff Automation GmbH PC ThinkPad E420 LenovoFigure 4
Experiment scene and relations of data stream between the experimental units.
(a)
Experiment scene (b)
Relations of data stream between the experimental unitsThe basic processes of robot control, data communication, and data acquisition are as follows:(a) The IPC, as the upper computer, writes the motion control commands of robot joint servo motors by programming software.(b) The upper computer transfers control commands to the servo driver, and the servo driver drives the motion of each joint servo motor according to the commands.(c) During the operation, the rotary encoder feedback data signals such as angular position and angular velocity to the servo driver.(d) Run data acquisition software on PC, establish communication, and collect the required physical data.In the experiment, only the waist joint moves; other joints do not move. The range of angular velocity of the waist joint is [0°/s, 100°/s]. In [0°/s, 1°/s], take a point every 0.2°/s, take a point every 1°/s in [1°/s, 5°/s], take a point every 10°/s in (5°/s, 80°/s), and measure the corresponding friction torque. The relationship between the two is shown in Figure5. Then, the friction parameters are identified by using the least square method, and the results are shown in Table 2.Table 2
Result of friction parameters of waist joint identification.
M f c 1 (N∙m) M f s 1 (N∙m) θ ˙ s 1 (rad∙s -1 ) σ 21 (N∙m∙rad -1 ) 0.827 1.276 15.271 0.056Figure 5
Relation between friction torque and angular velocity of waist joint.(4) Establishment of Shoulder and Elbow Joint Friction Models. The system abstract model of shoulder and elbow joints is shown in Figure 2. The reduction ratio of their reducer isn2=n3=171.(1) Establishment of Shoulder Joint Friction Model. Combining (10), (22) with (23), the shoulder joint friction model is established:(29)τf2=Mfc2+Mfs2-Mfc2e-n2θ˙2/θ˙s22sgnθ˙2+σ22n2θ˙2+Mfc2+Mfs2-Mfc2e-θ˙2/θ˙s22sgnθ˙2+σ22θ˙2n2where θ˙2 is angular velocity of the shoulder joint; Mfs2, Mfc2, θ˙s2, and σ22 need to be identified; τf2 is friction torque of the shoulder joint.(2) Parameter Identification of Shoulder Joint Friction Model. Unlike the waist joint, the shaft of the shoulder joint is perpendicular to the direction of gravity. Therefore, when shoulder joint moves at a constant velocity, it needs to overcome both the friction torque and the gravitational torque of the load. In order to measure the data of friction torque, the reciprocating constant velocity motion test method is adopted here [30]; namely, the following two groups of experiments are conducted.Experiment I. The shoulder joint constant speed moves fromθ2 toθ2’, and the output torqueu1 of the motor is measured.Experiment II. The shoulder joint constant speed moves fromθ2’ toθ2, and the output torqueu2 of the motor is measured.Assumption. The output torque of the motor in Experiment I is(30)u1=τgθ2+Ffθ˙2where τg(θ2) is gravitational torque; Ff(θ˙2) is friction torque.The output torque of the motor in Experiment II is(31)u2=τgθ2+Ff-θ˙2Because every joint of MD1200-YJ palletizing robot adopts high precision RV reducer; the friction difference in the positive and negative directions is very small. Therefore,(32)Ff-θ˙2=-Ffθ˙2From equations (30) ~ (32), it can be obtained that the friction torque and the weight torque are, respectively:(33)Ffθ˙2=12u1-u2(34)τgθ2=12u1+u2According to the experimental data, the friction torque at different angular velocities can be obtained from equation (33).In the experiment, only the shoulder joint moves, while the other joints do not move. Within the angular velocity range of the shoulder joint [0°/s, 75°/s], a point is taken every 0.2°/s within [0°/s, 1°/s], and a point is taken every 1°/s within [1°/s, 5°/s]. In the (5°/s, 75°/s), a point is taken every 10°/s, and the corresponding friction torque is measured. The relationship between the two is shown in Figure6. The friction parameters are then identified, and the results are shown in Table 3.Table 3
Result of static parameters identification.
M f c 2 (N∙m) M f s 2 (N∙m) θ ˙ s 2 (rad∙s -1 ) σ 22 (N∙m∙rad -1 )  0.565 3.031 26.103 0.057Figure 6
Relation between friction torque and angular velocity of shoulder joint.(3) Establishment of Elbow Joint Friction Model. The elbow joint friction model is as follows:(35)τf3=Mfc3+Mfs3-Mfc3e-n3θ˙3/θ˙s32sgn⁡θ˙3⁡+σ23n3θ˙3+Mfc3+Mfs3-Mfc3e-θ˙3/θ˙s32sgn⁡θ˙3+σ23θ˙3n3where θ˙3 is angular velocity of the shoulder joint; Mfs3, Mfc3, θ˙s3, and σ23 need to be identified; τf3 is friction torque of the shoulder joint.The shaft of elbow joint is also perpendicular to the direction of gravity. The measurement of friction torque and the recognition of friction parameters are the same as that of shoulder joint. Here, the experimental data is directly given, as shown in Figure7 and friction parameters, as shown in Table 4.Table 4
Result of static parameters identification.
M f c 3 (N∙m) M f s 3 (N∙m) θ ˙ s 3 (rad∙s -1 ) σ 23 (N∙m∙rad -1 )  0.015 0.259 16.703 0.0066Figure 7
Relation between friction torque and angular velocity of elbow joint.(5) The Dynamic Model of MD1200-YJ Palletizing Robot considering Joint Friction Is Established. Because the model of the friction torque derived in this paper is to convert the friction torque of the joint system to the motor shaft; therefore, it needs to be converted to the output end of the joint; there is(36)FSiθ˙it=τfi·niF S i ( θ ˙ i ( t ) ) is joint friction torque; ni is reduction ratio of each joint.On the basis of equation (1), the friction torque term FSi(θ˙i(t)) is added to obtain the dynamic model considering friction of jointi: (37)τi∗t=Hiθitθ¨it+Ciθit,θ˙it+Giθit+FSiθ˙itwhere τi∗ is joint driving torque considering joint friction.


## 2.2.3. Establishment of Torque Model of Balance Spring Cylinder

The balancing spring cylinder of MD1200-YJ palletizing robot has inner and outer two springs, and the spring parameters are shown in Table5.Table 5
Main parameters of spring.
Spring Wire diameter d t (mm) Mean diameter D (mm) Pitch p (mm) Effective number n Free length L 0 (mm) Single ring stiffness K (N/mm) Inner spring 10 70 24.2 22 543 288 Outer spring 14 100 34.8 16 580 375The force of the spring cylinder is shown in Figure8. The cylinder barrel of the balance spring cylinder is connected with the upper big arm by hinge (at point B), and its cylinder rod is connected with the waist mounting bracket by hinge (at point A). The distance between the spring cylinder mounting hole on the waist mounting bracket and the shoulder joint shaft isd (230mm), and the angle between the big arm and the vertical direction isθ2. When the big arm moves from the OB’ vertical position to the OB position, the balance spring cylinder extends with the rotation of the big arm (that is, the increase ofθ2), while the internal spring is compressed, generating the balance spring force FT.Figure 8
Force analysis on spring balance mechanism.Spring stiffnessK:(38)KI=288nI(39)KO=375nOTwo springs are installed in parallel, the total spring stiffnessK: (40)K=KI+KOSpring preload:(41)F0I=KI·X0I(42)F0O=KO·X0OTwo springs are installed in parallel, the total preloadF0:(43)F0=F0I+F0OWorking length of the spring LAB:(44)LAB=d2+L22-2dL2cos⁡θ2Spring deformation ΔX:(45)ΔX=LAB-L2-dThen, spring tension force generated in the working processF is (46)F=F0+K·ΔXForce arm of spring tensionlf is(47)lf=L2dsin⁡θ2LABBalance torque of the springMT is(48)MTθ2=F·lf=F0+Kd2+L22-2dL2cos⁡θ2-L2+dL2dsin⁡θ2d2+L22-2dL2cos⁡θ2Therefore, the dynamic model of shoulder joint considering joint friction and spring moment is established.(49)τ~2t=H2θ2tθ¨2t+C2θ2t,θ˙2t+G2θ2t+FS2θ˙2t+MTθ2t


## 2.3. Establishment of Energy Consumption Model of MD1200-YJ Palletizing Robot

In this paper, the energy consumption of the robot joint is divided into the useful work to drive the load and dissipation energy due to friction and heat loss of the motor.Since this study is carried out in the no-load condition, compared with the other three joints, the load of the wrist joint is small and its movement speed is not high, so its energy consumption is less [18]. In order to simplify problem, energy consumption of the wrist joint is ignored. Therefore, the total energy consumption is only the sum of the energy consumption of the waist joint (i =1), shoulder joint (i =2), and elbow joint (i =3).The following energy consumption model can be established for the motion of jointi system int0 totf time:(50)WSi=∫t0tfτmotori·θ˙i·ni+ii2Ridtwhere WSi, τmotori, θ˙i, ni,ii, and Ri are, respectively, energy consumption, motor torque, angular velocity, reduction ratio of reducer, effective value of motor current, torque constant of motor, and equivalent resistance of motor of jointi.The motor torque model of the jointi can be obtained from equations (37) and (49)(51)τmotori=τiniAndτmotori has the following relationship with the effective currentii of the motor(52)τmotori=ki·iiwhereki is the torque constant of the motor of jointi.The thermal power model of joint motori can be obtained by joule’s law combined with equation (52)(53)PRi=τmotoriki2RiThus, the total energy consumption model of palletizing robot int0 totf time is established.(54)WS=∑i=13WSi=∑i=13∫t0tfτmotori·θ˙i·ni+τmotoriki2RidtFriction energy model:(55)Pfi=τfi·θ˙i·niThe total dissipative energy consumption model of palletizing robot int0 totf time is established:(56)Wloss=∑i=13∫t0tfPfi+PRidt=∑i=13∫t0tfτfi·θ˙i·ni+τmotoriki2RidtThe mechanical power of jointi(57)PEi=τmotori·θ˙i·niThe total mechanical energy model of palletizing robot int0 totf time is established:(58)WE=∑i=13∫t0tfPEidt=∑i=13∫t0tfτmotori·θ˙i·nidt


## 2.4. Validation and Analysis of Dynamic Model and Energy Consumption Model

Based on the example of 3-4-5 polynomial motion law, the multifactor dynamic model and energy consumption model are verified and analyzed by experiment and simulation.


### 2.4.1. Description of Palletizing Robot Parameters and Calculation Examples

The main parameters of MD1200-YJ palletizing robot are shown in Tables6–8.Table 6
Main function size and material of palletizing robot.
Part name Main size parameter Material Big arm 1200 mm Zl110 Forearm 1200 mm Zl110 Drive arm of forearm 360 mm QT-5007 Driving link of forearm 1200 mm Q235A Horizontal holding link 1 1200 mm Q235A Horizontal holding link 2 1200 mm Q235ATable 7
Parameters of motors and reducer.
Joint Servo motor Reducer Waist joint Rated power: 5.5 kWTorque constant: 2.37 N·m/A Type: RV320CSpeed ratio: 139.1 Shoulder joint Rated power: 5.5 kWTorque constant: 2.37 N•m/A Type: RV160ESpeed ratio: 171 Elbow joint Rated power: 5.5 kWTorque constant: 2.37 N•m/A Type: RV160ESpeed ratio: 171Table 8
Mass of the main part of palletizing robot.
Part name Mass (kg) Part name Mass (kg) Waist joint 486.085 Big arm 107.105 Forearm 67.202 Drive arm of forearm 13.513 End effector 31.10716 Driving link of forearm 31.107Problem Description of Calculation Example. In the basic coordinate system, the end position and posture of the palletizing robot are represented by coordinates (X,Y,Z,Φ). The torque and mechanical power of the joint motor at the end of the joint from the base coordinateS0 (1200mm, 1800mm, 800mm, 0rad) to the pose coordinateSf (-1200mm, 1200mm, 500mm, 0rad) are discussed. The interval of starting and ending time of motion is [0s, 3.5s]. Through the forward and inverse kinematics solution, the joint coordinates corresponding toS0 andSf, and the motion parameters at the start and stop time, as shown in Tables 9 and 10, the whole process of no-load operation is conducted.Table 9
Joint position ofS0.
Joint Angular position(rad) Angular velocity(rad/s) Angular acceleration(rad/s 2 ) θ 1(t0) 0.9828 0 0 θ 2(t0) 1.10251 0 0 θ 3(t0) 2.6202 0 0Table 10
Joint position ofSf.
Joint Angular position(rad) Angular velocity(rad/s) Angular acceleration(rad/s 2 ) θ 1(tf) 2.3562 0 0 θ 2(tf) 1.1519 0 0 θ 3(tf) 2.2028 0 0According to the conditions given in this example, the trajectory planning of 3-4-5 polynomial motion law commonly used in engineering is carried out, and the expression is as follows:(59)θ1=0.9828+0.3203t3-0.1373t4+0.0157t5θ2=1.0251+0.0296t3-0.0127t4+0.0014t5θ2′=2.6202-0.0974t3+0.0417t4-0.0048t5t∈03.5The curves of end trajectory and the angular position of waist joint, shoulder joint, and elbow joint are shown in Figures9 and 10.Figure 9
End trajectory of 3-4-5 polynomial.Figure 10
Angular displacement trajectory of 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow joint


### 2.4.2. Validation of Dynamic Model and Energy Consumption Model

According to the 3-4-5 polynomial motion law and the known motion parameters, the motion of each joint of the palletizing robot is controlled, and the data acquisition software is used to collect the torque and mechanical power data of each joint motor within a motion period.It should be noted that since the data acquisition software cannot directly measure the heat loss of the motor and friction energy consumption, therefore, the total energy consumption of the palletizing robot cannot be directly measured. However, according to equation (54), the accuracy of the energy consumption model can be indirectly verified by comparing the experimental value of mechanical power with the simulated value.The experimental data collected show a trend of high frequency oscillation, but the overall trend is regular. Therefore, the fitting curve of the experimental data is obtained by the polynomial curve fitting method, as shown in Figure11, and the fitting curve of experimental data is compared with the simulation curve of the theoretical model, as shown in Figure 12.Figure 11
Experimental fitting curves of joint motor torque.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 12
Comparison of joint motor torque experimental curve with simulation curve.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointIt can be seen from Figure12 that the variation trend of the motor torque simulation curve of the three joints is similar to that of the experimental curve, which indicates that the multifactor dynamics model has higher accuracy. However, it can be seen that there is a certain error between the two, which may be caused by (1) the scale parameters and inertia parameters of all parts of the simulation model are measured by the 3D model, which are different from the corresponding parts of the actual prototype; (2) there will be errors in the experimental curve fitting.The influence of friction torque and spring torque on output torque of shoulder joint motor is analyzed, as shown in Figure13. It shows that the friction torque and spring torque have obvious effects on the output torque of the motor. The conventional rigid body dynamics model cannot fully reflect components of the joint torque. The addition of friction torque and spring torque makes the output torque of the shoulder joint closer to the real situation.Figure 13
Comparison of shoulder joint motor torque experimental curve with simulation curve under different conditions.
(a)
No spring and no friction (b)
No spring and consider friction (c)
Consider spring and no friction (d)
Consider spring and frictionNext, the mechanical power data of each joint are analyzed. The mechanical power experimental data of waist joint, shoulder joint, and elbow joint are obtained through data acquisition software. The phenomenon of high frequency oscillation still exists, and the fitting curve of the experimental data is obtained by the polynomial curve fitting method, as shown in Figure14, and compared with the simulation curve of mechanical power, as shown in Figure 15.Figure 14
Experimental fitting curves of joint mechanical power.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 15
Comparison of joint mechanical power experimental curve with simulation curve.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointAs can be seen from Figure15, the variation trend of the mechanical power simulation curve of the three joints is similar to that of the experimental curve. As can be seen from Table 11, the theoretical value of mechanical power consumption differs from the experimental value by 80.2988 J, with an error of about 7.69%, which proves the accuracy of the energy consumption model.Table 11
Mechanical power of joints (J).
Comparison object Waist joint Shoulder joint Elbow joint Total mechanical power consumption Simulation 837.0672 103.0638 104.6278 1044.8 Experimental 776.6805 82.0542 105.7665 964.5012 Error (%) 7.21408 20.38504 1.088334 7.685567


## 2.4.1. Description of Palletizing Robot Parameters and Calculation Examples

The main parameters of MD1200-YJ palletizing robot are shown in Tables6–8.Table 6
Main function size and material of palletizing robot.
Part name Main size parameter Material Big arm 1200 mm Zl110 Forearm 1200 mm Zl110 Drive arm of forearm 360 mm QT-5007 Driving link of forearm 1200 mm Q235A Horizontal holding link 1 1200 mm Q235A Horizontal holding link 2 1200 mm Q235ATable 7
Parameters of motors and reducer.
Joint Servo motor Reducer Waist joint Rated power: 5.5 kWTorque constant: 2.37 N·m/A Type: RV320CSpeed ratio: 139.1 Shoulder joint Rated power: 5.5 kWTorque constant: 2.37 N•m/A Type: RV160ESpeed ratio: 171 Elbow joint Rated power: 5.5 kWTorque constant: 2.37 N•m/A Type: RV160ESpeed ratio: 171Table 8
Mass of the main part of palletizing robot.
Part name Mass (kg) Part name Mass (kg) Waist joint 486.085 Big arm 107.105 Forearm 67.202 Drive arm of forearm 13.513 End effector 31.10716 Driving link of forearm 31.107Problem Description of Calculation Example. In the basic coordinate system, the end position and posture of the palletizing robot are represented by coordinates (X,Y,Z,Φ). The torque and mechanical power of the joint motor at the end of the joint from the base coordinateS0 (1200mm, 1800mm, 800mm, 0rad) to the pose coordinateSf (-1200mm, 1200mm, 500mm, 0rad) are discussed. The interval of starting and ending time of motion is [0s, 3.5s]. Through the forward and inverse kinematics solution, the joint coordinates corresponding toS0 andSf, and the motion parameters at the start and stop time, as shown in Tables 9 and 10, the whole process of no-load operation is conducted.Table 9
Joint position ofS0.
Joint Angular position(rad) Angular velocity(rad/s) Angular acceleration(rad/s 2 ) θ 1(t0) 0.9828 0 0 θ 2(t0) 1.10251 0 0 θ 3(t0) 2.6202 0 0Table 10
Joint position ofSf.
Joint Angular position(rad) Angular velocity(rad/s) Angular acceleration(rad/s 2 ) θ 1(tf) 2.3562 0 0 θ 2(tf) 1.1519 0 0 θ 3(tf) 2.2028 0 0According to the conditions given in this example, the trajectory planning of 3-4-5 polynomial motion law commonly used in engineering is carried out, and the expression is as follows:(59)θ1=0.9828+0.3203t3-0.1373t4+0.0157t5θ2=1.0251+0.0296t3-0.0127t4+0.0014t5θ2′=2.6202-0.0974t3+0.0417t4-0.0048t5t∈03.5The curves of end trajectory and the angular position of waist joint, shoulder joint, and elbow joint are shown in Figures9 and 10.Figure 9
End trajectory of 3-4-5 polynomial.Figure 10
Angular displacement trajectory of 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow joint


## 2.4.2. Validation of Dynamic Model and Energy Consumption Model

According to the 3-4-5 polynomial motion law and the known motion parameters, the motion of each joint of the palletizing robot is controlled, and the data acquisition software is used to collect the torque and mechanical power data of each joint motor within a motion period.It should be noted that since the data acquisition software cannot directly measure the heat loss of the motor and friction energy consumption, therefore, the total energy consumption of the palletizing robot cannot be directly measured. However, according to equation (54), the accuracy of the energy consumption model can be indirectly verified by comparing the experimental value of mechanical power with the simulated value.The experimental data collected show a trend of high frequency oscillation, but the overall trend is regular. Therefore, the fitting curve of the experimental data is obtained by the polynomial curve fitting method, as shown in Figure11, and the fitting curve of experimental data is compared with the simulation curve of the theoretical model, as shown in Figure 12.Figure 11
Experimental fitting curves of joint motor torque.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 12
Comparison of joint motor torque experimental curve with simulation curve.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointIt can be seen from Figure12 that the variation trend of the motor torque simulation curve of the three joints is similar to that of the experimental curve, which indicates that the multifactor dynamics model has higher accuracy. However, it can be seen that there is a certain error between the two, which may be caused by (1) the scale parameters and inertia parameters of all parts of the simulation model are measured by the 3D model, which are different from the corresponding parts of the actual prototype; (2) there will be errors in the experimental curve fitting.The influence of friction torque and spring torque on output torque of shoulder joint motor is analyzed, as shown in Figure13. It shows that the friction torque and spring torque have obvious effects on the output torque of the motor. The conventional rigid body dynamics model cannot fully reflect components of the joint torque. The addition of friction torque and spring torque makes the output torque of the shoulder joint closer to the real situation.Figure 13
Comparison of shoulder joint motor torque experimental curve with simulation curve under different conditions.
(a)
No spring and no friction (b)
No spring and consider friction (c)
Consider spring and no friction (d)
Consider spring and frictionNext, the mechanical power data of each joint are analyzed. The mechanical power experimental data of waist joint, shoulder joint, and elbow joint are obtained through data acquisition software. The phenomenon of high frequency oscillation still exists, and the fitting curve of the experimental data is obtained by the polynomial curve fitting method, as shown in Figure14, and compared with the simulation curve of mechanical power, as shown in Figure 15.Figure 14
Experimental fitting curves of joint mechanical power.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 15
Comparison of joint mechanical power experimental curve with simulation curve.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointAs can be seen from Figure15, the variation trend of the mechanical power simulation curve of the three joints is similar to that of the experimental curve. As can be seen from Table 11, the theoretical value of mechanical power consumption differs from the experimental value by 80.2988 J, with an error of about 7.69%, which proves the accuracy of the energy consumption model.Table 11
Mechanical power of joints (J).
Comparison object Waist joint Shoulder joint Elbow joint Total mechanical power consumption Simulation 837.0672 103.0638 104.6278 1044.8 Experimental 776.6805 82.0542 105.7665 964.5012 Error (%) 7.21408 20.38504 1.088334 7.685567


## 3. Minimum Energy Optimization of the Joint Drive System of Palletizing Robot

The MD1200-YJ palletizing robot is a kind of multijoint tandem robot, and there is complex coupling relationship between each joint. Single-joint system is the basic component of multijoint system, and its structure is relatively simple. For convenience of description, a single-joint system is taken as an example to illustrate the research process, and then this method is extended to the multijoint system of palletizing robot; in other words, based on the energy consumption optimization method of single-joint system, energy consumption optimization of each joint of palletizing robot is conducted, respectively.


### 3.1. Research Method of Minimum Energy Optimization for Palletizing Robot Joint System

Taking a single-joint system as an example, the theory and method of its minimum energy consumption optimization are illustrated, and provides the research basis for the multi-joint minimum energy consumption optimization of palletizing robot.


#### 3.1.1. Motion Law of Joints Are Constructed Based on Fourier Series Approximation Method

By (1), (37), (49), and (54), the palletizing robot joints system energy consumption is a function of motion parameters. Therefore, the minimum energy optimization of joint system can be attributed to the problem of finding a motion rule that minimizes the total energy consumption of joint system. BecauseWS is a nonlinear function of motion parameters, which is hard to obtain the analytic solution, therefore, a Fourier series approximation method is used to solve the problem.The method is as follows: assume the conclusionθ(t) obtained by optimization, but its expression has not yet been determined, can use Fourier series structure in advance a function expression of  θ^(t); in this way, the complex variational solving problem is transformed into the problem of finding the optimal parameter with known expression structure, namely, coefficients of Fourier series θ^(t) under the lowest energy consumption are optimized.There are numerous motion paths of the joint system from the angular positionθ0 at the moment oft0 to the angular positionθf at the moment oftf, and these paths are represented here by θ(t).The definitional domain ofθ(t) is [t0,tf], and θ(t) are extended to the continuous periodic function θ~(t) that satisfies the Dirichlet condition [30]. θ~(t) can be represented by Fourier series f(t). Within interval [t0,tf]:(60)θt=θ~t=ftThis moment, the energy consumption optimization problem of the joint driving system is transformed into the search for Fourier seriesf(t) to obtain the minimum energy consumption WSmin problem of the joint system in the time interval [t0,tf], and the coefficients of the Fourier series are the optimization variables.The expression ofθ(t) is(61)θt=a0+∑n=1∞ancos⁡nω0t+bnsin⁡nω0twhere an and bn are the coefficients of the Fourier series, a0 is the constant term, and ω0 is the angular frequency, ω0=π/(tf-t0).The constraint conditions that this formula needs to satisfy:(62)θt0=θ0,θtf=θfθ˙t0=0,θ˙tf=0θ¨t0=0,(63)θ¨tf=0The Fourier series is infinite series. In order to facilitate the practical operation, the finite terms are intercepted to approximateθ(t). When the number of intercepted terms is large enough, a better approximation can be obtained.Therefore, formula (62) is truncated into formula (64) to obtain the expression of the function representing any path:(64)θt=a0+∑n=1Mancos⁡nω0t+bnsin⁡nω0twhereM is the number of truncated terms (for convenience, this formula is also referred to as the Fourier series).The first and second derivatives ofθ(t) are solved to obtain the corresponding expressions of angular velocity and angular acceleration:(65)θ˙t=∑n=1Mnω0ancos⁡nω0t+π2+nω0bnsin⁡nω0t+π2(66)θ¨t=∑n=1Mn2ω02ancos⁡nω0t+π+n2ω02bnsin⁡nω0t+πTherefore, the specific expression of the function representing any path is obtained.


#### 3.1.2. Genetic Algorithm and Energy Optimization of Joint Drive System

Genetic Algorithm (GA) [31, 32], originated from Darwin’s theory of evolution, is an adaptive probabilistic optimization method for complex and nonlinear system optimization.The MD1200-YJ palletizing robot is a kind of nonlinear and strongly coupled complex system, and the minimum energy consumption optimization problem of its joint systems can be solved by genetic algorithm.Here, a single-joint drive system is used as an example to introduce the minimum energy consumption optimization method.MATLAB genetic algorithm toolbox is used to optimize the GA function call format as follows [33]:(67)x,fval,exitflag,output,population,scores=GA@fitnessfun,nvars,Aineq,Bineq,Aeq,Beq,LB,UB,nonlcon,options;The input parameters are shown in Table12.Table 12
Input parameters.
Parameter Meaning Parameter Meaning fitnessfun Fitness function nvars Number of optimized variables Aineq Linear inequality constraint coefficient matrix Bineq Linear constraint inequality conditional vector Aeq Equality constraint coefficient matrix Beq Equality constraint vector LB Lower limit of variable UB Upper limit of variable nonlcon Nonlinear constraint functions Options Algorithm properties set optionsThe optimization goal is energy consumptionWSmin, which is used as the fitness function of genetic algorithm. The optimized variable is the coefficient (a0,a1,a2,…,am,b1,b2,…,bm) of the Fourier series in formula (64), and the input parameter nvars value is 2×M+1.Formula (64), (65), and (66) are expressed as matrix forms:(68)A·x=Bwhere(69)A=cos⁡ω0t⋯cos⁡Mω0tsin⁡ω0t⋯sinMω0t1ω0cos⁡ω0t+π2⋯Mω0cos⁡Mω0t+π2ω0sin⁡ω0t+π2⋯Mω0sin⁡Mω0t+π20ω02cos⁡ω0t+π⋯Mω02cos2⁡Mω0t+πω02sin⁡ω0t+π⋯Mω02sin⁡Mω0t+π0(70)x=a1a2⋯aMb1b2⋯bMa0T(71)B=θtθ˙tθ¨tTCombining the equality constraint conditions of equation (62), the equality constraint coefficient matrix Aeq and the equality constraint vector Beq can be obtained.The upper and lower limits of the variables in Table12 are set to the default values of UB and LB.The “Options” in Table13 are set through the gaoptimset function; call format is:(72)options=gaoptimset(‘property1’,‘value1’,‘property2’,‘value2’,…,‘propertyn’,‘valuen’);Table 13
Property of GA.
Property Value P o p u l a t i o n S z i e 50 G e n e r a t i o n s 200 S t a l l G e n L i m i t 150 F i t n e s s L i m i n t 1e-6Some of the algorithm properties are set as shown in Table13.Run GA function and the output parameters are shown in Table14.Table 14
Output parameters.
Parameter Meaning Parameter Meaning X Final point Output Performance of each generation fval Final fitness function value Population Final population exitflag Stop flag Scores Final score


### 3.2. Minimum Energy Optimization of the Joint Drive System of Palletizing Robot

Based on the above theories and methods of the joint minimum energy optimization, and combining the multifactor dynamics model, energy consumption model, and corresponding motion law, the minimum energy consumption optimization of MD1200-YJ palletizing robot multiple joints system is studied.


#### 3.2.1. Establishment of Joint Motion Equation Based on Fourier Series Method

It is assumed that end effector of the palletizing robot moves from the positionx0,y0,z0,β0Tat the time oft0 to the position xf,yf,zf,βfTat the time oftf in the operating space. When starting and stopping, the corresponding relation between the palletizing robot’s operating space and the motion parameters under joint space is shown in Table 15.Table 15
Relationship between the operating space and the joint space.
Items Operating space Joint space Start position x 0 , y 0 , z 0 , β 0 T θ 10 , θ 20 , θ 30 , θ 40 T Stop position x f , y f , z f , β f T θ 1 f , θ 2 f , θ 3 f , θ 4 f T Start velocity x ˙ 0 , y ˙ 0 , z ˙ 0 , β ˙ 0 T θ ˙ 10 , θ ˙ 20 , θ ˙ 30 , θ ˙ 40 T Stop velocity x ˙ f , y ˙ f , z ˙ f , β ˙ f T θ ˙ 1 f , θ ˙ 2 f , θ ˙ 3 f , θ ˙ 4 f T Start acceleration x ¨ 0 , y ¨ 0 , z ¨ 0 , β ¨ 0 T θ ¨ 10 , θ ¨ 20 , θ ¨ 30 , θ ¨ 40 T Stop acceleration x ¨ f , y ¨ f , z ¨ f , β ¨ f T θ ¨ 1 f , θ ¨ 2 f , θ ¨ 3 f , θ ¨ 4 f TIn order to obtain the minimum energy consumption ofWSmin, the Fourier series of the four joints of the palletizing robot is constructed, andM finite terms are intercepted, and the motion equation of the joints is obtained as follows:(73)θ1t=α0+∑n=1Mαncos⁡nω0t+αn′sin⁡nω0tθ2t=β0+∑n=1Mβncos⁡nω0t+βn′sin⁡nω0tθ3t=η0+∑n=1Mηncos⁡nω0t+ηn′sin⁡nω0tθ4t=γ0+∑n=1Mγncos⁡nω0t+γn′sin⁡nω0twhere αn, αn′, βn, βn′, ηn, ηn′, γn and γn′, respectively, represent the coefficients of the Fourier series mentioned above, ω0=π/(tf-t0).The four joint motion equations must satisfy the constraint conditions:(74)θit0=θi0,θitf=θifθ˙it0=0,θ˙itf=0θ¨it0=0,θ¨itf=0i=1,2,3,4


#### 3.2.2. Optimization of Fourier Series Genetic Algorithm

The fitness function of genetic algorithm is energy consumptionWS, and the variables of this fitness function are αn, αn′, βn, βn′, ηn, ηn′, γn, and γn′, (n=1,2,3…,M), and four constant terms α0, β0, η0, γ0 and the number of variables is 8×M+4. Therefore, the complexity of optimization will depend on the number of Fourier series termsM in formula (73). The process of optimization is shown in Figures 16 and 17.Figure 16
Process of direct optimization.Figure 17
Process of the Calculation GA individual fitness value.For the convenience of genetic algorithm optimization, formula (73) is written into matrix form:(75)A′·x′=B′where(76)A′t=At0000At0000At0000AtHere,A(t) is shown in the equation (69).(77)x′=αi,αi′,α0,βi,βi′,β0,ηi,η′,η0,γi,γi′,γ0T(78)B′=θ1θ˙1θ¨1θ2θ˙2θ¨2θ3θ˙3θ¨3θ4θ˙4θ¨4TAccording to equation (74), the constraint conditions in the genetic algorithm are determined. Assuming the equality constraint matrix I and the equality constraint vector J, there are(79)I=A′t0A′tf,J=B′t0B′tfT


#### 3.2.3. Optimization Results and Evaluation

In order to explore the law and characteristics of the method, the genetic algorithm is applied to optimize the item numberM0=3~600, respectively, and the energy consumption under the aforementioned 3-4-5 polynomial motion law is compared in the same conditions. The total energy consumption values corresponding to the optimization results of different item number are shown in Figure 18. Here, in order to simplify problem, the total energy consumption is still the sum of the energy consumption of the waist joint, shoulder joint, and elbow joint.Figure 18
Total Energy ofM0=3~600 with 3-4-5 polynomial.(1) Analysis of Optimization Results. The following can be seen from Figure 18.( 1 ) ExceptM0=3, the energy consumption values of the motion law obtained by the optimization of other terms are less than that of 3-4-5 polynomial motion law.( 2 ) As the value ofM0 increases, its total energy consumption decreases first and then increases and finally keeps fluctuating within a certain range; moreover, when the value ofM0=8, the total energy consumption reaches the minimum value of 917.5 J.(2) Comparison between the Optimal Motion Law of Fourier Series and the 3-4-5 Polynomial Motion Law. Here, Fourier series trajectories ofM0=8 (the smallest total energy consumption) andM0=600 (the largest number of intercepting terms) are compared with 3-4-5 polynomial motion law. The comparison of the terminal trajectories is shown in Figure 19.Figure 19
Trajectory ofM0=8,M0=600, and 3-4-5 polynomial.It can be seen from Figure19 that the terminal trajectory under theM0=600 motion law is similar to that under the 3-4-5 polynomial motion law, and the smoothness is higher than that under theM0=8 motion law. In the following part, the motion characteristics, energy consumption characteristics, and dynamic characteristics ofM0=8,M0=600, and 3-4-5 polynomial motion laws are evaluated, respectively. Here, the dynamic characteristics can be represented by the joint jerk, which is the first derivative of joint angular acceleration, and is an important indicator to evaluate the stability of robot motion. The smaller the joint jerk value, the smaller the change rate of angular acceleration, the better the stability of motion.Figures20–26 are the comparison diagrams of joint trajectory, mechanical power, dissipation power, total power, and jerk ofM0=8 andM0=600, and 3-4-5 polynomial motion law, respectively.Figure 20
Angular displacement trajectory ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 21
Angular velocity trajectory ofM0=8,M0=600, and 3-4-5 and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 22
Angular acceleration trajectory ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 23
Mechanical power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 24
Dissipation power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 25
Total power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 26
Jerk ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointThe energy consumption ofM0=8,M0=600, and 3-4-5 polynomial are compared, as shown in Table 16.Table 16
Comparison with energy consumption ofM0=8, M0=600, and 3-4-5 polynomial (J).
Motion law Total energy Dissipation energy Total energy of joints Dissipation energy of joints Waist Shoulder Elbow Waist Shoulder Elbow M 0=8 917.54 757.44 698.76 108.45 110.33 661.32 75.466 20.661 M 0=600 955.05 790.54 733.37 111.16 110.52 696.36 77.527 16.646 3-4-5 1055.9 894.33 841.11 109.50 105.34 799.20 77.210 17.925The following conclusions can be drawn from Figures20–26 and Table 16.(1) As can be seen from Figure20, the smoothness of the Fourier series trajectories of lower number of terms is not as good as that of the 3-4-5 polynomial trajectories, while the smoothness of the Fourier series trajectories of higher number of terms is significantly improved.(2) As can be seen from Figures22 and 26, as the number of terms increases, the number of high harmonics increases, the angular acceleration, and jerk of the Fourier series motion law show obvious oscillations, especially violent oscillation of jerk ofM0=600. Therefore, the practicability of the high term number Fourier series trajectory needs further consideration. The jerk ofM0=8 motion law is greater than 3-4-5 polynomial motion law, which indicates that its motion stability is lower. In the next step, further optimization research can be conducted by combining other optimization indexes or constraints.(3) It can be seen from Figures21, 23, 24, and 25 that they are significantly affected by angular velocity and have similar change trend. It can be seen from Table 16, compared with the 3-4-5 polynomial motion law, the energy consumption of the motion law obtained by the Fourier series optimization method is significantly reduced, among which the total energy consumption ofM0=8 reduced by 13.1% and its dissipation energy reduced by 15.3%. The effectiveness of the optimization method is verified. Since this research is carried out under no-load conditions, the dissipation energy consumption accounted for larger proportion of the total energy consumption, which also reflected the principle that “a big horse hauls a small gharry is wasting” commonly used in engineering. It can be speculated that when the palletizing robot drives a large load similar to the load capacity, the useful work will take a large proportion.


## 3.1. Research Method of Minimum Energy Optimization for Palletizing Robot Joint System

Taking a single-joint system as an example, the theory and method of its minimum energy consumption optimization are illustrated, and provides the research basis for the multi-joint minimum energy consumption optimization of palletizing robot.


### 3.1.1. Motion Law of Joints Are Constructed Based on Fourier Series Approximation Method

By (1), (37), (49), and (54), the palletizing robot joints system energy consumption is a function of motion parameters. Therefore, the minimum energy optimization of joint system can be attributed to the problem of finding a motion rule that minimizes the total energy consumption of joint system. BecauseWS is a nonlinear function of motion parameters, which is hard to obtain the analytic solution, therefore, a Fourier series approximation method is used to solve the problem.The method is as follows: assume the conclusionθ(t) obtained by optimization, but its expression has not yet been determined, can use Fourier series structure in advance a function expression of  θ^(t); in this way, the complex variational solving problem is transformed into the problem of finding the optimal parameter with known expression structure, namely, coefficients of Fourier series θ^(t) under the lowest energy consumption are optimized.There are numerous motion paths of the joint system from the angular positionθ0 at the moment oft0 to the angular positionθf at the moment oftf, and these paths are represented here by θ(t).The definitional domain ofθ(t) is [t0,tf], and θ(t) are extended to the continuous periodic function θ~(t) that satisfies the Dirichlet condition [30]. θ~(t) can be represented by Fourier series f(t). Within interval [t0,tf]:(60)θt=θ~t=ftThis moment, the energy consumption optimization problem of the joint driving system is transformed into the search for Fourier seriesf(t) to obtain the minimum energy consumption WSmin problem of the joint system in the time interval [t0,tf], and the coefficients of the Fourier series are the optimization variables.The expression ofθ(t) is(61)θt=a0+∑n=1∞ancos⁡nω0t+bnsin⁡nω0twhere an and bn are the coefficients of the Fourier series, a0 is the constant term, and ω0 is the angular frequency, ω0=π/(tf-t0).The constraint conditions that this formula needs to satisfy:(62)θt0=θ0,θtf=θfθ˙t0=0,θ˙tf=0θ¨t0=0,(63)θ¨tf=0The Fourier series is infinite series. In order to facilitate the practical operation, the finite terms are intercepted to approximateθ(t). When the number of intercepted terms is large enough, a better approximation can be obtained.Therefore, formula (62) is truncated into formula (64) to obtain the expression of the function representing any path:(64)θt=a0+∑n=1Mancos⁡nω0t+bnsin⁡nω0twhereM is the number of truncated terms (for convenience, this formula is also referred to as the Fourier series).The first and second derivatives ofθ(t) are solved to obtain the corresponding expressions of angular velocity and angular acceleration:(65)θ˙t=∑n=1Mnω0ancos⁡nω0t+π2+nω0bnsin⁡nω0t+π2(66)θ¨t=∑n=1Mn2ω02ancos⁡nω0t+π+n2ω02bnsin⁡nω0t+πTherefore, the specific expression of the function representing any path is obtained.


### 3.1.2. Genetic Algorithm and Energy Optimization of Joint Drive System

Genetic Algorithm (GA) [31, 32], originated from Darwin’s theory of evolution, is an adaptive probabilistic optimization method for complex and nonlinear system optimization.The MD1200-YJ palletizing robot is a kind of nonlinear and strongly coupled complex system, and the minimum energy consumption optimization problem of its joint systems can be solved by genetic algorithm.Here, a single-joint drive system is used as an example to introduce the minimum energy consumption optimization method.MATLAB genetic algorithm toolbox is used to optimize the GA function call format as follows [33]:(67)x,fval,exitflag,output,population,scores=GA@fitnessfun,nvars,Aineq,Bineq,Aeq,Beq,LB,UB,nonlcon,options;The input parameters are shown in Table12.Table 12
Input parameters.
Parameter Meaning Parameter Meaning fitnessfun Fitness function nvars Number of optimized variables Aineq Linear inequality constraint coefficient matrix Bineq Linear constraint inequality conditional vector Aeq Equality constraint coefficient matrix Beq Equality constraint vector LB Lower limit of variable UB Upper limit of variable nonlcon Nonlinear constraint functions Options Algorithm properties set optionsThe optimization goal is energy consumptionWSmin, which is used as the fitness function of genetic algorithm. The optimized variable is the coefficient (a0,a1,a2,…,am,b1,b2,…,bm) of the Fourier series in formula (64), and the input parameter nvars value is 2×M+1.Formula (64), (65), and (66) are expressed as matrix forms:(68)A·x=Bwhere(69)A=cos⁡ω0t⋯cos⁡Mω0tsin⁡ω0t⋯sinMω0t1ω0cos⁡ω0t+π2⋯Mω0cos⁡Mω0t+π2ω0sin⁡ω0t+π2⋯Mω0sin⁡Mω0t+π20ω02cos⁡ω0t+π⋯Mω02cos2⁡Mω0t+πω02sin⁡ω0t+π⋯Mω02sin⁡Mω0t+π0(70)x=a1a2⋯aMb1b2⋯bMa0T(71)B=θtθ˙tθ¨tTCombining the equality constraint conditions of equation (62), the equality constraint coefficient matrix Aeq and the equality constraint vector Beq can be obtained.The upper and lower limits of the variables in Table12 are set to the default values of UB and LB.The “Options” in Table13 are set through the gaoptimset function; call format is:(72)options=gaoptimset(‘property1’,‘value1’,‘property2’,‘value2’,…,‘propertyn’,‘valuen’);Table 13
Property of GA.
Property Value P o p u l a t i o n S z i e 50 G e n e r a t i o n s 200 S t a l l G e n L i m i t 150 F i t n e s s L i m i n t 1e-6Some of the algorithm properties are set as shown in Table13.Run GA function and the output parameters are shown in Table14.Table 14
Output parameters.
Parameter Meaning Parameter Meaning X Final point Output Performance of each generation fval Final fitness function value Population Final population exitflag Stop flag Scores Final score


## 3.1.1. Motion Law of Joints Are Constructed Based on Fourier Series Approximation Method

By (1), (37), (49), and (54), the palletizing robot joints system energy consumption is a function of motion parameters. Therefore, the minimum energy optimization of joint system can be attributed to the problem of finding a motion rule that minimizes the total energy consumption of joint system. BecauseWS is a nonlinear function of motion parameters, which is hard to obtain the analytic solution, therefore, a Fourier series approximation method is used to solve the problem.The method is as follows: assume the conclusionθ(t) obtained by optimization, but its expression has not yet been determined, can use Fourier series structure in advance a function expression of  θ^(t); in this way, the complex variational solving problem is transformed into the problem of finding the optimal parameter with known expression structure, namely, coefficients of Fourier series θ^(t) under the lowest energy consumption are optimized.There are numerous motion paths of the joint system from the angular positionθ0 at the moment oft0 to the angular positionθf at the moment oftf, and these paths are represented here by θ(t).The definitional domain ofθ(t) is [t0,tf], and θ(t) are extended to the continuous periodic function θ~(t) that satisfies the Dirichlet condition [30]. θ~(t) can be represented by Fourier series f(t). Within interval [t0,tf]:(60)θt=θ~t=ftThis moment, the energy consumption optimization problem of the joint driving system is transformed into the search for Fourier seriesf(t) to obtain the minimum energy consumption WSmin problem of the joint system in the time interval [t0,tf], and the coefficients of the Fourier series are the optimization variables.The expression ofθ(t) is(61)θt=a0+∑n=1∞ancos⁡nω0t+bnsin⁡nω0twhere an and bn are the coefficients of the Fourier series, a0 is the constant term, and ω0 is the angular frequency, ω0=π/(tf-t0).The constraint conditions that this formula needs to satisfy:(62)θt0=θ0,θtf=θfθ˙t0=0,θ˙tf=0θ¨t0=0,(63)θ¨tf=0The Fourier series is infinite series. In order to facilitate the practical operation, the finite terms are intercepted to approximateθ(t). When the number of intercepted terms is large enough, a better approximation can be obtained.Therefore, formula (62) is truncated into formula (64) to obtain the expression of the function representing any path:(64)θt=a0+∑n=1Mancos⁡nω0t+bnsin⁡nω0twhereM is the number of truncated terms (for convenience, this formula is also referred to as the Fourier series).The first and second derivatives ofθ(t) are solved to obtain the corresponding expressions of angular velocity and angular acceleration:(65)θ˙t=∑n=1Mnω0ancos⁡nω0t+π2+nω0bnsin⁡nω0t+π2(66)θ¨t=∑n=1Mn2ω02ancos⁡nω0t+π+n2ω02bnsin⁡nω0t+πTherefore, the specific expression of the function representing any path is obtained.


## 3.1.2. Genetic Algorithm and Energy Optimization of Joint Drive System

Genetic Algorithm (GA) [31, 32], originated from Darwin’s theory of evolution, is an adaptive probabilistic optimization method for complex and nonlinear system optimization.The MD1200-YJ palletizing robot is a kind of nonlinear and strongly coupled complex system, and the minimum energy consumption optimization problem of its joint systems can be solved by genetic algorithm.Here, a single-joint drive system is used as an example to introduce the minimum energy consumption optimization method.MATLAB genetic algorithm toolbox is used to optimize the GA function call format as follows [33]:(67)x,fval,exitflag,output,population,scores=GA@fitnessfun,nvars,Aineq,Bineq,Aeq,Beq,LB,UB,nonlcon,options;The input parameters are shown in Table12.Table 12
Input parameters.
Parameter Meaning Parameter Meaning fitnessfun Fitness function nvars Number of optimized variables Aineq Linear inequality constraint coefficient matrix Bineq Linear constraint inequality conditional vector Aeq Equality constraint coefficient matrix Beq Equality constraint vector LB Lower limit of variable UB Upper limit of variable nonlcon Nonlinear constraint functions Options Algorithm properties set optionsThe optimization goal is energy consumptionWSmin, which is used as the fitness function of genetic algorithm. The optimized variable is the coefficient (a0,a1,a2,…,am,b1,b2,…,bm) of the Fourier series in formula (64), and the input parameter nvars value is 2×M+1.Formula (64), (65), and (66) are expressed as matrix forms:(68)A·x=Bwhere(69)A=cos⁡ω0t⋯cos⁡Mω0tsin⁡ω0t⋯sinMω0t1ω0cos⁡ω0t+π2⋯Mω0cos⁡Mω0t+π2ω0sin⁡ω0t+π2⋯Mω0sin⁡Mω0t+π20ω02cos⁡ω0t+π⋯Mω02cos2⁡Mω0t+πω02sin⁡ω0t+π⋯Mω02sin⁡Mω0t+π0(70)x=a1a2⋯aMb1b2⋯bMa0T(71)B=θtθ˙tθ¨tTCombining the equality constraint conditions of equation (62), the equality constraint coefficient matrix Aeq and the equality constraint vector Beq can be obtained.The upper and lower limits of the variables in Table12 are set to the default values of UB and LB.The “Options” in Table13 are set through the gaoptimset function; call format is:(72)options=gaoptimset(‘property1’,‘value1’,‘property2’,‘value2’,…,‘propertyn’,‘valuen’);Table 13
Property of GA.
Property Value P o p u l a t i o n S z i e 50 G e n e r a t i o n s 200 S t a l l G e n L i m i t 150 F i t n e s s L i m i n t 1e-6Some of the algorithm properties are set as shown in Table13.Run GA function and the output parameters are shown in Table14.Table 14
Output parameters.
Parameter Meaning Parameter Meaning X Final point Output Performance of each generation fval Final fitness function value Population Final population exitflag Stop flag Scores Final score


## 3.2. Minimum Energy Optimization of the Joint Drive System of Palletizing Robot

Based on the above theories and methods of the joint minimum energy optimization, and combining the multifactor dynamics model, energy consumption model, and corresponding motion law, the minimum energy consumption optimization of MD1200-YJ palletizing robot multiple joints system is studied.


### 3.2.1. Establishment of Joint Motion Equation Based on Fourier Series Method

It is assumed that end effector of the palletizing robot moves from the positionx0,y0,z0,β0Tat the time oft0 to the position xf,yf,zf,βfTat the time oftf in the operating space. When starting and stopping, the corresponding relation between the palletizing robot’s operating space and the motion parameters under joint space is shown in Table 15.Table 15
Relationship between the operating space and the joint space.
Items Operating space Joint space Start position x 0 , y 0 , z 0 , β 0 T θ 10 , θ 20 , θ 30 , θ 40 T Stop position x f , y f , z f , β f T θ 1 f , θ 2 f , θ 3 f , θ 4 f T Start velocity x ˙ 0 , y ˙ 0 , z ˙ 0 , β ˙ 0 T θ ˙ 10 , θ ˙ 20 , θ ˙ 30 , θ ˙ 40 T Stop velocity x ˙ f , y ˙ f , z ˙ f , β ˙ f T θ ˙ 1 f , θ ˙ 2 f , θ ˙ 3 f , θ ˙ 4 f T Start acceleration x ¨ 0 , y ¨ 0 , z ¨ 0 , β ¨ 0 T θ ¨ 10 , θ ¨ 20 , θ ¨ 30 , θ ¨ 40 T Stop acceleration x ¨ f , y ¨ f , z ¨ f , β ¨ f T θ ¨ 1 f , θ ¨ 2 f , θ ¨ 3 f , θ ¨ 4 f TIn order to obtain the minimum energy consumption ofWSmin, the Fourier series of the four joints of the palletizing robot is constructed, andM finite terms are intercepted, and the motion equation of the joints is obtained as follows:(73)θ1t=α0+∑n=1Mαncos⁡nω0t+αn′sin⁡nω0tθ2t=β0+∑n=1Mβncos⁡nω0t+βn′sin⁡nω0tθ3t=η0+∑n=1Mηncos⁡nω0t+ηn′sin⁡nω0tθ4t=γ0+∑n=1Mγncos⁡nω0t+γn′sin⁡nω0twhere αn, αn′, βn, βn′, ηn, ηn′, γn and γn′, respectively, represent the coefficients of the Fourier series mentioned above, ω0=π/(tf-t0).The four joint motion equations must satisfy the constraint conditions:(74)θit0=θi0,θitf=θifθ˙it0=0,θ˙itf=0θ¨it0=0,θ¨itf=0i=1,2,3,4


### 3.2.2. Optimization of Fourier Series Genetic Algorithm

The fitness function of genetic algorithm is energy consumptionWS, and the variables of this fitness function are αn, αn′, βn, βn′, ηn, ηn′, γn, and γn′, (n=1,2,3…,M), and four constant terms α0, β0, η0, γ0 and the number of variables is 8×M+4. Therefore, the complexity of optimization will depend on the number of Fourier series termsM in formula (73). The process of optimization is shown in Figures 16 and 17.Figure 16
Process of direct optimization.Figure 17
Process of the Calculation GA individual fitness value.For the convenience of genetic algorithm optimization, formula (73) is written into matrix form:(75)A′·x′=B′where(76)A′t=At0000At0000At0000AtHere,A(t) is shown in the equation (69).(77)x′=αi,αi′,α0,βi,βi′,β0,ηi,η′,η0,γi,γi′,γ0T(78)B′=θ1θ˙1θ¨1θ2θ˙2θ¨2θ3θ˙3θ¨3θ4θ˙4θ¨4TAccording to equation (74), the constraint conditions in the genetic algorithm are determined. Assuming the equality constraint matrix I and the equality constraint vector J, there are(79)I=A′t0A′tf,J=B′t0B′tfT


### 3.2.3. Optimization Results and Evaluation

In order to explore the law and characteristics of the method, the genetic algorithm is applied to optimize the item numberM0=3~600, respectively, and the energy consumption under the aforementioned 3-4-5 polynomial motion law is compared in the same conditions. The total energy consumption values corresponding to the optimization results of different item number are shown in Figure 18. Here, in order to simplify problem, the total energy consumption is still the sum of the energy consumption of the waist joint, shoulder joint, and elbow joint.Figure 18
Total Energy ofM0=3~600 with 3-4-5 polynomial.(1) Analysis of Optimization Results. The following can be seen from Figure 18.( 1 ) ExceptM0=3, the energy consumption values of the motion law obtained by the optimization of other terms are less than that of 3-4-5 polynomial motion law.( 2 ) As the value ofM0 increases, its total energy consumption decreases first and then increases and finally keeps fluctuating within a certain range; moreover, when the value ofM0=8, the total energy consumption reaches the minimum value of 917.5 J.(2) Comparison between the Optimal Motion Law of Fourier Series and the 3-4-5 Polynomial Motion Law. Here, Fourier series trajectories ofM0=8 (the smallest total energy consumption) andM0=600 (the largest number of intercepting terms) are compared with 3-4-5 polynomial motion law. The comparison of the terminal trajectories is shown in Figure 19.Figure 19
Trajectory ofM0=8,M0=600, and 3-4-5 polynomial.It can be seen from Figure19 that the terminal trajectory under theM0=600 motion law is similar to that under the 3-4-5 polynomial motion law, and the smoothness is higher than that under theM0=8 motion law. In the following part, the motion characteristics, energy consumption characteristics, and dynamic characteristics ofM0=8,M0=600, and 3-4-5 polynomial motion laws are evaluated, respectively. Here, the dynamic characteristics can be represented by the joint jerk, which is the first derivative of joint angular acceleration, and is an important indicator to evaluate the stability of robot motion. The smaller the joint jerk value, the smaller the change rate of angular acceleration, the better the stability of motion.Figures20–26 are the comparison diagrams of joint trajectory, mechanical power, dissipation power, total power, and jerk ofM0=8 andM0=600, and 3-4-5 polynomial motion law, respectively.Figure 20
Angular displacement trajectory ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 21
Angular velocity trajectory ofM0=8,M0=600, and 3-4-5 and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 22
Angular acceleration trajectory ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 23
Mechanical power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 24
Dissipation power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 25
Total power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 26
Jerk ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointThe energy consumption ofM0=8,M0=600, and 3-4-5 polynomial are compared, as shown in Table 16.Table 16
Comparison with energy consumption ofM0=8, M0=600, and 3-4-5 polynomial (J).
Motion law Total energy Dissipation energy Total energy of joints Dissipation energy of joints Waist Shoulder Elbow Waist Shoulder Elbow M 0=8 917.54 757.44 698.76 108.45 110.33 661.32 75.466 20.661 M 0=600 955.05 790.54 733.37 111.16 110.52 696.36 77.527 16.646 3-4-5 1055.9 894.33 841.11 109.50 105.34 799.20 77.210 17.925The following conclusions can be drawn from Figures20–26 and Table 16.(1) As can be seen from Figure20, the smoothness of the Fourier series trajectories of lower number of terms is not as good as that of the 3-4-5 polynomial trajectories, while the smoothness of the Fourier series trajectories of higher number of terms is significantly improved.(2) As can be seen from Figures22 and 26, as the number of terms increases, the number of high harmonics increases, the angular acceleration, and jerk of the Fourier series motion law show obvious oscillations, especially violent oscillation of jerk ofM0=600. Therefore, the practicability of the high term number Fourier series trajectory needs further consideration. The jerk ofM0=8 motion law is greater than 3-4-5 polynomial motion law, which indicates that its motion stability is lower. In the next step, further optimization research can be conducted by combining other optimization indexes or constraints.(3) It can be seen from Figures21, 23, 24, and 25 that they are significantly affected by angular velocity and have similar change trend. It can be seen from Table 16, compared with the 3-4-5 polynomial motion law, the energy consumption of the motion law obtained by the Fourier series optimization method is significantly reduced, among which the total energy consumption ofM0=8 reduced by 13.1% and its dissipation energy reduced by 15.3%. The effectiveness of the optimization method is verified. Since this research is carried out under no-load conditions, the dissipation energy consumption accounted for larger proportion of the total energy consumption, which also reflected the principle that “a big horse hauls a small gharry is wasting” commonly used in engineering. It can be speculated that when the palletizing robot drives a large load similar to the load capacity, the useful work will take a large proportion.


## 3.2.1. Establishment of Joint Motion Equation Based on Fourier Series Method

It is assumed that end effector of the palletizing robot moves from the positionx0,y0,z0,β0Tat the time oft0 to the position xf,yf,zf,βfTat the time oftf in the operating space. When starting and stopping, the corresponding relation between the palletizing robot’s operating space and the motion parameters under joint space is shown in Table 15.Table 15
Relationship between the operating space and the joint space.
Items Operating space Joint space Start position x 0 , y 0 , z 0 , β 0 T θ 10 , θ 20 , θ 30 , θ 40 T Stop position x f , y f , z f , β f T θ 1 f , θ 2 f , θ 3 f , θ 4 f T Start velocity x ˙ 0 , y ˙ 0 , z ˙ 0 , β ˙ 0 T θ ˙ 10 , θ ˙ 20 , θ ˙ 30 , θ ˙ 40 T Stop velocity x ˙ f , y ˙ f , z ˙ f , β ˙ f T θ ˙ 1 f , θ ˙ 2 f , θ ˙ 3 f , θ ˙ 4 f T Start acceleration x ¨ 0 , y ¨ 0 , z ¨ 0 , β ¨ 0 T θ ¨ 10 , θ ¨ 20 , θ ¨ 30 , θ ¨ 40 T Stop acceleration x ¨ f , y ¨ f , z ¨ f , β ¨ f T θ ¨ 1 f , θ ¨ 2 f , θ ¨ 3 f , θ ¨ 4 f TIn order to obtain the minimum energy consumption ofWSmin, the Fourier series of the four joints of the palletizing robot is constructed, andM finite terms are intercepted, and the motion equation of the joints is obtained as follows:(73)θ1t=α0+∑n=1Mαncos⁡nω0t+αn′sin⁡nω0tθ2t=β0+∑n=1Mβncos⁡nω0t+βn′sin⁡nω0tθ3t=η0+∑n=1Mηncos⁡nω0t+ηn′sin⁡nω0tθ4t=γ0+∑n=1Mγncos⁡nω0t+γn′sin⁡nω0twhere αn, αn′, βn, βn′, ηn, ηn′, γn and γn′, respectively, represent the coefficients of the Fourier series mentioned above, ω0=π/(tf-t0).The four joint motion equations must satisfy the constraint conditions:(74)θit0=θi0,θitf=θifθ˙it0=0,θ˙itf=0θ¨it0=0,θ¨itf=0i=1,2,3,4


## 3.2.2. Optimization of Fourier Series Genetic Algorithm

The fitness function of genetic algorithm is energy consumptionWS, and the variables of this fitness function are αn, αn′, βn, βn′, ηn, ηn′, γn, and γn′, (n=1,2,3…,M), and four constant terms α0, β0, η0, γ0 and the number of variables is 8×M+4. Therefore, the complexity of optimization will depend on the number of Fourier series termsM in formula (73). The process of optimization is shown in Figures 16 and 17.Figure 16
Process of direct optimization.Figure 17
Process of the Calculation GA individual fitness value.For the convenience of genetic algorithm optimization, formula (73) is written into matrix form:(75)A′·x′=B′where(76)A′t=At0000At0000At0000AtHere,A(t) is shown in the equation (69).(77)x′=αi,αi′,α0,βi,βi′,β0,ηi,η′,η0,γi,γi′,γ0T(78)B′=θ1θ˙1θ¨1θ2θ˙2θ¨2θ3θ˙3θ¨3θ4θ˙4θ¨4TAccording to equation (74), the constraint conditions in the genetic algorithm are determined. Assuming the equality constraint matrix I and the equality constraint vector J, there are(79)I=A′t0A′tf,J=B′t0B′tfT


## 3.2.3. Optimization Results and Evaluation

In order to explore the law and characteristics of the method, the genetic algorithm is applied to optimize the item numberM0=3~600, respectively, and the energy consumption under the aforementioned 3-4-5 polynomial motion law is compared in the same conditions. The total energy consumption values corresponding to the optimization results of different item number are shown in Figure 18. Here, in order to simplify problem, the total energy consumption is still the sum of the energy consumption of the waist joint, shoulder joint, and elbow joint.Figure 18
Total Energy ofM0=3~600 with 3-4-5 polynomial.(1) Analysis of Optimization Results. The following can be seen from Figure 18.( 1 ) ExceptM0=3, the energy consumption values of the motion law obtained by the optimization of other terms are less than that of 3-4-5 polynomial motion law.( 2 ) As the value ofM0 increases, its total energy consumption decreases first and then increases and finally keeps fluctuating within a certain range; moreover, when the value ofM0=8, the total energy consumption reaches the minimum value of 917.5 J.(2) Comparison between the Optimal Motion Law of Fourier Series and the 3-4-5 Polynomial Motion Law. Here, Fourier series trajectories ofM0=8 (the smallest total energy consumption) andM0=600 (the largest number of intercepting terms) are compared with 3-4-5 polynomial motion law. The comparison of the terminal trajectories is shown in Figure 19.Figure 19
Trajectory ofM0=8,M0=600, and 3-4-5 polynomial.It can be seen from Figure19 that the terminal trajectory under theM0=600 motion law is similar to that under the 3-4-5 polynomial motion law, and the smoothness is higher than that under theM0=8 motion law. In the following part, the motion characteristics, energy consumption characteristics, and dynamic characteristics ofM0=8,M0=600, and 3-4-5 polynomial motion laws are evaluated, respectively. Here, the dynamic characteristics can be represented by the joint jerk, which is the first derivative of joint angular acceleration, and is an important indicator to evaluate the stability of robot motion. The smaller the joint jerk value, the smaller the change rate of angular acceleration, the better the stability of motion.Figures20–26 are the comparison diagrams of joint trajectory, mechanical power, dissipation power, total power, and jerk ofM0=8 andM0=600, and 3-4-5 polynomial motion law, respectively.Figure 20
Angular displacement trajectory ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 21
Angular velocity trajectory ofM0=8,M0=600, and 3-4-5 and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 22
Angular acceleration trajectory ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 23
Mechanical power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 24
Dissipation power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 25
Total power ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointFigure 26
Jerk ofM0=8,M0=600, and 3-4-5 polynomial.
(a)
Waist joint (b)
Shoulder joint (c)
Elbow jointThe energy consumption ofM0=8,M0=600, and 3-4-5 polynomial are compared, as shown in Table 16.Table 16
Comparison with energy consumption ofM0=8, M0=600, and 3-4-5 polynomial (J).
Motion law Total energy Dissipation energy Total energy of joints Dissipation energy of joints Waist Shoulder Elbow Waist Shoulder Elbow M 0=8 917.54 757.44 698.76 108.45 110.33 661.32 75.466 20.661 M 0=600 955.05 790.54 733.37 111.16 110.52 696.36 77.527 16.646 3-4-5 1055.9 894.33 841.11 109.50 105.34 799.20 77.210 17.925The following conclusions can be drawn from Figures20–26 and Table 16.(1) As can be seen from Figure20, the smoothness of the Fourier series trajectories of lower number of terms is not as good as that of the 3-4-5 polynomial trajectories, while the smoothness of the Fourier series trajectories of higher number of terms is significantly improved.(2) As can be seen from Figures22 and 26, as the number of terms increases, the number of high harmonics increases, the angular acceleration, and jerk of the Fourier series motion law show obvious oscillations, especially violent oscillation of jerk ofM0=600. Therefore, the practicability of the high term number Fourier series trajectory needs further consideration. The jerk ofM0=8 motion law is greater than 3-4-5 polynomial motion law, which indicates that its motion stability is lower. In the next step, further optimization research can be conducted by combining other optimization indexes or constraints.(3) It can be seen from Figures21, 23, 24, and 25 that they are significantly affected by angular velocity and have similar change trend. It can be seen from Table 16, compared with the 3-4-5 polynomial motion law, the energy consumption of the motion law obtained by the Fourier series optimization method is significantly reduced, among which the total energy consumption ofM0=8 reduced by 13.1% and its dissipation energy reduced by 15.3%. The effectiveness of the optimization method is verified. Since this research is carried out under no-load conditions, the dissipation energy consumption accounted for larger proportion of the total energy consumption, which also reflected the principle that “a big horse hauls a small gharry is wasting” commonly used in engineering. It can be speculated that when the palletizing robot drives a large load similar to the load capacity, the useful work will take a large proportion.


## 4. Conclusion

The minimum energy trajectory optimization problem of the joint drive system of MD1200-YJ palletizing robot is researched in this paper.(1) Rigid body dynamics model for MD1200-YJ palletizing robot is derived using Kane equation method; based on the Stribeck friction model and parameter identification theory, the friction torque model of each joint of the palletizing robot is established, the spring balance torque model is established by the motion parameters of the shoulder joint and the structure parameters of the balance spring cylinder. And then a multifactor dynamic model of the palletizing robot is established based on the above models. And by combining joint torque, friction torque, motion parameter with joule’s law, the useful work model, dissipation energy model, and total energy consumption model of the palletizing robot are established. Based on a calculating example under the 3-4-5 polynomial motion law the correctness of multifactor dynamic model and energy consumption model is verified through experiments.(2) The joint motion law is constructed by using the Fourier series approximation method, and taking the coefficients of Fourier series as optimization variables, the minimum energy consumption is taken as the optimization objective, the motion parameters of initial and final position, and running time constant as constraint conditions, and then the minimum energy consumption trajectory of the joint system is optimized by the genetic algorithm.(3) Through the calculating example simulation, the optimizedM0=8 andM0=600 Fourier series and the 3-4-5 polynomial motion law are comprehensively evaluated. The energy consumption of the optimized Fourier series is generally lower than that of the 3-4-5 polynomial motion law, especially energy consumption ofM0=8, which is total energy consumption reduced by 13.1% and dissipation energy reduced by 15.3% under the same conditions, the effectiveness of the optimization method is verified. However,M0=8 has a larger jerk value, and its motion stability is less than 3-4-5 polynomial motion law. Moreover, because the coefficients of the Fourier series (that is, the number of optimization variables) increase multiply with the number of terms, the computational complexity of the optimization solution is large, the convergence speed is slow, and the configuration requirements of the computing equipment are higher, to some extent, this method may be restricted in practical engineering application. Therefore, in the next step, further optimization research can be conducted: on the one hand, combine more optimization indexes or constraint conditions to reduce the jerk value of each joint under the Fourier series motion law; on the other hand, based on the Fourier series optimization method, some new optimization strategies need be explored to improve the calculation efficiency and reduce the difficulty of solving, to improve the practicability of this optimization method.



---

*Source: 7247093-2018-12-03.xml*