# Multiobjective Trajectory Optimization and Adaptive Backstepping Control for Rubber Unstacking Robot Based on RFWNN Method

Authors: Le Liang; Yanjie Liu; Hao Xu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Multiobjective trajectory optimization and adaptive backstepping control method based on recursive fuzzy wavelet neural network (RFWNN) are proposed to solve the problem of dynamic modeling uncertainties and strong external disturbance of the rubber unstacking robot during recycling process. First, according to the rubber viscoelastic properties, the Hunt-Crossley nonlinear model is used to construct the robot dynamics model. Then, combined with the dynamic model and the recycling process characteristics, the multiobjective trajectory optimization of the rubber unstacking robot is carried out for the operational efficiency, the running trajectory smoothness, and the energy consumption. Based on the trajectory optimization results, the adaptive backstepping control method based on RFWNN is adopted. The RFWNN method is applied in the main controller to cope with time-varying uncertainties of the robot dynamic system. Simultaneously, an adaptive robust control law is developed to eliminate inevitable approximation errors and unknown disturbances and relax the requirement for prior knowledge of the controlled system. Finally, the validity of the proposed control strategy is verified by experiment.


---

## Body

## 1. Introduction

Countless rubber in the tire manufacturing and other industries has been applied, such as natural rubber and synthetic rubber. Due to the need for a large number of manual participation and taking into account the labor intensity and low efficiency, it is urgent to use robotic technology instead of artificial one to complete the rubber unstacking process. As the heat refining of the rubber block needs to be operated before the mixer processing, the rubber block will have a certain degree of melting and resulting in uneven surface of the rubber block which leads to the fact that adsorption type fixture cannot be used. Therefore, the way that the robot end effector directly is inserted into the rubber block is utilized to complete the unstacking process. During this process, the robot will be in sharp contact with the rubber, time-varying contact force will seriously interfere with the normal operation of the robot, and conventional trajectory planning and control methods are difficult to meet the system performance requirements.Efficiency, trajectory smoothness, and energy consumption are important performance indicators of rubber unstacking robot in trajectory planning. Particularly during recycling process, due to the viscoelastic characteristics of rubber, the faster the insertion velocity of the end effector, the greater the reaction force. According to the characteristics of robot operation, performance requirements, and corresponding constraints, multiobjective optimization model is established to obtain the optimal trajectory of the rubber unstacking robot during recycling process. Many researches have been done in this area [1–3]. Piazzi and Visioli developed an approach based on interval analysis to find the global minimum-jerk (MJ) trajectory of a robot manipulator within a joint space scheme using cubic splines [4]. Chettibi et al. discussed the problem of minimum cost trajectory planning for robotic manipulators, which put a few free via points uniformly distributed in time [5]. Korayem et al. converted dynamic load-carrying capacity problem into a trajectory optimization problem of cable-suspended parallel robots which is fundamentally a constrained nonlinear optimization problem [6]. Abu-Dakka et al. addressed an indirect method for trajectory planning for industrial robots operating in workspaces with obstacles using an evolutionary algorithm [7].In the aspect of control strategy, due to the influence of time-varying nonlinear contact force on the robot during the process of inserting the rubber, the modeling error caused by the stiffness of the robot bar, the friction of the joint, and the external disturbance are highlighted, which seriously affect the robot’s optimal trajectory tracking accuracy and control performance. The adaptive robust controllers based on hybrid neural networks can exhibit good properties as an attempt to cope with the problem of modeling uncertainty parameters and external disturbances. The neural networks can deal with the unknown dynamics problems of the rubber unstacking robot control system by exploiting their universal approximation ability [8–10]. At the same time, the approximate error of the neural network and the interference of the external system are estimated and compensated by adding the adaptive robust term in the controller [11–13]. Chairez applied differential neural network approach with activation functions described by wavelets in a state observation problem when the dynamic model of a plant contains uncertainties [14]. Yoo et al. proposed a method for the robust control of flexible-joint robots with model uncertainties which combined the adaptive dynamic surface control technique and the self-recurrent wavelet neural network [15]. Boukattaya et al. studied the trajectory tracking control problem of mobile manipulators subject to nonholonomic constraints, operating in task space, with the presence of external torque disturbances and dynamic uncertainties [16].The main objective of this work is related to multiobjective trajectory optimization and adaptive robust trajectory tracking control of the rubber unstacking robot during the recycling process. First, according to the nonlinearity, viscoelasticity, and anisotropy of the rubber, the Hunt-Crossley foundation model was used to establish the mechanical model while the end effector is inserted into the rubber block. Then, aimed at the performance indexes such as operation efficiency, running trajectory smoothness, and energy consumption, the problem of motion planning is transformed into multiobjective optimization problem. The position of the floating via points for the B-spline trajectory and the running time between adjacent points are optimized, so the optimal trajectory in joint space is obtained. Finally, a recursive fuzzy wavelet neural network is proposed to estimate the uncertainty model of the system by combining fuzzy and wavelet neural networks. At the same time, adaptive backstepping control method is adopted to eliminate the influence of uncertainties such as estimation error and external disturbance on the tracking accuracy of the system trajectory. The control performance is verified by simulation analysis.


## 2. Robot System and Dynamic Model

The robot unstacking system is shown in Figure1 which includes unstacking robot, visual system, ultrasonic system, six-axis force/torque sensor, and end effector with pointy tip. The dynamic model of robot unstacking process is shown as follows:(1)Mθθ¨+Cθ,θ˙θ˙+Gθ+JTθFs=τ,where θ, θ˙, and θ¨ are the robot joint angle, angular velocity, and angular acceleration, respectively, M(θ) is the inertia matrix, C(θ,θ˙) is Coriolis force item and centrifugal force item, G(θ) is gravity item, τ is the generalized input torque, J(θ) is Jacobian matrix for the conversion of joint space to Cartesian space, andFs represents the external force in Cartesian space when the robot is inserted into the rubber. Table 1 shows the definition of main structural physical parameters of the unstacking robot.Table 1
Definition of main structural physical parameters.
Physical meaning
Symbol
Value
Waist joint angle/°
θ
1
−165~165
Big arm joint angle/°
θ
2
5~130
Forearm joint angle/°
θ
3
60~200
Length of big arm/mm
l
1
945
Length of forearm/mm
l
2
1025
Length of parallel arm/mm
l
3
400Figure 1
Composition of the robot unstacking system.As can be seen from (1), M(θ), C(θ,θ˙), G(θ), and J(θ) can be calculated according to the robot kinematics and dynamics; however, it is difficult to describeFs with a simple linear mechanics model because the rubber has the characteristics of nonlinearity, viscoelasticity, and anisotropy. To obtain a more accurate viscoelastic model, Hunt and Crossley [17] made a nonlinear combination of spring and damping systems, and a nonlinear model was established as follows:(2)Fet=Kxβt+λxβx˙t,where K and λ are the elastic coefficient and the nonlinear damping coefficient, respectively, x(t) is the displacement inserted in the rubber, x˙(t) is insert speed, and β is the positive scalar.The mechanical model of the robot end effector insertion rubber block process was established on the basis of the selected HC nonlinear basic model and combined with the mechanical characteristics of rubber block. The insertion process is divided into two stages, namely, tip part insertion stage and rod part insertion stage. Figure2 showed the former mechanical model. Where R is the maximum bottom radius of the tip cone, H is the tip cone height, α is the tip cone angle, x˙f(t), x˙r(t) are the velocity component parallel and perpendicular to the direction of the conical bus, respectively,fT is the tip friction of robot end effector, fTh and fTr are, respectively, the axial and radial frictional forces,FT is tip positive pressure of the end effector (FT′ is the reaction force), fTh and fTr are, respectively, the axial and radial positive pressure components, andμ is the friction coefficient between the rubber and actuator.Figure 2
Mechanical model of the tip insertion stage.According to Figure2 and HC nonlinear viscoelastic model, the mechanical model equations of the tip part insertion stage xt≤H are as follows:(3)ri=xitan⁡αcos⁡α,x˙rt=x˙tsin⁡α,Fei=Kriβ+λriβx˙rt,FT=∫0xtFei2πxitan⁡αcos⁡αdxi=2πβ+2K+λx˙tsin⁡αtan⁡αcos⁡αβ+1xtβ+2,fT=μFT,Fst=fTcos⁡α+FTsin⁡α.Figure3 showed the mechanical model of the rod part insertion stage, in which fR is the rod friction of robot end effector, and FR is the rod positive pressure of the end effector (FR′ is the reaction force).Figure 3
Mechanical model of the rod insertion stage.According to Figure3 and HC nonlinear viscoelastic model, the mechanical model equations of the tip insertion stage (x(t)>H) are as follows:(4)FT=2πβ+2K+λx˙tsin⁡αtan⁡αcos⁡αβ+1Hβ+2,fT=μFT,FRi=KRβ,FR=∫0xt-HFRi2πRdxi=2πKRβ+1xt-H,fR=μFR,Fst=fTcos⁡α+FTsin⁡α+fR.The mechanical model of the robot end effector inserted into the rubber process can be obtained after solving (3) and (4):(5)Fst=2πβ+2tan⁡αcos⁡αβ+1μcos⁡α+sin⁡α·K+λsin⁡αx˙txtβ+2,xt≤H;2πβ+2tan⁡αcos⁡αβ+1μcos⁡α+sin⁡α·K+λsin⁡αx˙tHβ+2+2πμKRβ+1xt-H,xt>H,whereα,H,R, andμare known quantities and their values areα= 20°,H= 48 mm,R= 17.5 mm, andμ = 0.45,x(t) and x˙(t) are input variables,K,β, andλcan be obtained by experimental parameter identification, and their values areK = 16.32 KN/m,β= 1.15, andλ = 7.36 KN∗s/m2.


## 3. Multiobjective Trajectory Optimization Algorithm

### 3.1. B-Spline Trajectory Planning in Joint Space

Robot joint space trajectory planning is based on a series of key constraint points to plan robot motion trajectory which areθmi=θmiti+1,i=0,1,…,n,j=1,2,…,M, where θmi(t) represents theith trajectory segment of joint m. θmi is the angle at the i+1 via point of jointm that is the angle at ti+1, andM is the number of robot joints. The number of constraint points for any joint trajectory curve is n+1, which are starting point, end point, and n-1 via points, respectively.The unstacking robot joint space trajectory is described using the nonuniform B-spline curve. Based on determining the orderk and the node vectorU= u0,u1,…,un+k+1, B-spline curve recursive equation using the de Boer algorithm is achieved as follows [18]:(6)pu=∑j=i-k+1idjlNj,k-1u=⋯=dik,ui≤u<ui+1,djl=dj,l=01-αjldj-1l-1+αjldjl-1,l=1,2,…,k;j=i-k+l,…,i,αjl=u-ujuj+k+1-l-uj,where di∈RM×1 is the control vertices of the joint spline trajectory, u∈ui,ui+1 is the normalized time node vector of spline curve in segmenti, and Ni,k(u) is k-order specification B-spline basis function. According to the local support character of the B-spline curve, ther-order derivative of p(u) is obtained(7)pru=∑j=i-k+ridjrNj,k-ru,ui≤u<ui+1,djl=dj,l=0k+1-ldjl-1-dj-1l-1uj+k+1-l-uj,l=1,2,…,r;j=i-k+l,…,i.According to (7), robot joint space velocity, acceleration, and jerk can be obtained as follows:(8)θ˙mit=p′u=∑j=i-k+1idj1Nj,k-1u,θ¨mit=p′′u=∑j=i-k+2idj2Nj,k-2u,θ⃛mit=p′′′u=∑j=i-k+3idj3Nj,k-3u.


### 3.2. Optimization Model

During the robot recycling process, the trajectories of each joint are designed as segment B-spline curves. By using the same time interval parameterhi to couple the joint tracks to each other, whereθmi(t) represents theith trajectory of jointm, θmi, θ˙mi, and θ¨mi are, respectively, the angle, angle velocity, and acceleration at thei via point of jointm. The main constraints of the optimization process are set as follows.(
1) Position Constraint. According to the recycling characteristics, the rubber unstacking robot usually adopts the “door-shaped” trajectory to operate as shown in Figure 4. In order to avoid interference, the AB, CE, and GE segments are set as linear trajectories. In some former researches [1–7], the trajectory is usually optimized according to a number of fixed teaching points and the way of inserting fixed via points. The optimization method overreliance on the selection of the initial trajectory is mainly to optimize the segmentation trajectory between the fixed via points. However, since all the via points have been set in advance, the optimization range of the Cartesian space trajectory is limited and it is difficult to obtain the optimal trajectory.Figure 4
Recycling process trajectory before optimization.The trajectory optimization method proposed in this paper contains two movement forms: linear trajectories AB and GH which are limited by obstacle avoidance and recycling form; the free trajectory passes through the barrier point B, the highest point D, and the insertion point G which are optimized by the rectangular space composed of B, D, and G points; meanwhile, to avoid interference, it is necessary to ensure that the robot end effector is always moving in the space. The positional coordinates of the end effector in the Cartesian space along they- and z-axes can be solved by positive kinematics. It is only necessary to verify whether the limit position of each trajectory satisfies the constraint condition as compared with the coordinate value of the whole motion time period. When the velocity of the rubber unstacking robot along y- and z-axes is zero, the corresponding position is the limit position of the B-spline segment, which can be calculated as follows:(9)y˙it=cos⁡θ1itl1cos⁡θ2it-l3cos⁡θ3it+S1θ˙1it-l1sin⁡θ1itsin⁡θ2itθ˙2it+l3sin⁡θ1itsin⁡θ3itθ˙3it,z˙it=l1cos⁡θ2itθ˙2it-l3cos⁡θ3itθ˙3it.The limit position ofyi∗ and zi∗ of the spline trajectory of theith segment is obtained by the above equations, and the spatial position constraints are as follows:(10)yMIN≤yi∗≤yMAX,i=b,…,g,zMIN≤zi∗≤zMAX,i=b,…,g,whereb and g are the spline segments of the transition point B and G, respectively. Since the linear trajectory segments AB and GE are in the z-axis direction, the position constraint is as follows:(11)yi∗-ySTRAIGHT≤δMAX,i=1,…,b&g,…,n.(
2) Velocity, Acceleration, and Jerk Constraints. Under the condition of the drive system which has been set, the joint velocity and acceleration of the rubber unstacking robot are constrained to a certain extreme range and, at the same time, due to excessive acceleration, can cause the robot mechanical system resonance and increase the wear of the mechanical system, so the jerk is also used as one of the constraints. Any section trajectory of any joint should satisfy the constraint conditions as follows:(12)θ˙mLB≤θ˙mit≤θ˙mUB,θ¨mLB≤θ¨mit≤θ¨mUB,θ⃛mit≤JMAX,i=0,…,n,j=1,…,m,where θ˙mLB, θ˙mUB, θ¨mLB, and θ¨mUB are, respectively, the velocity and acceleration limit value of jointm and JMAX is the maximum jerk. In the process of optimization, the constraint will be transformed into a constraint on the first three derivatives of the B-spline control vertex dm as shown in (7).(
3) Torque Constraints. The joint output torque is limited by the robot drive system and needs to meet the corresponding constraints. According to (1) one can get(13)τ=Mθθ¨+Cθ,θ.θ.+Gθ+JTθFs≤τUB,where τUB is the output torque limit, which is defined as the 80% of the maximum torque of the joint motor.(
4) Normalized Torque Constraints. In order to ensure the smoothness of the robot motion, the normalized torque is taken as the constraint condition as follows:(14)∑m=1M1/H∫0Hτmit2dt1/N∑k=1NτmikΔT-τmik-1ΔT2≤MUP,whereN is the number of servo control cycles, ΔT is the robot servo cycles, and MUP is the maximum value.In order to improve the working efficiency, take the work time of the rubber unstacking robot during recycling process as the optimization goal. At the same time, taking into account the trajectory smoothness and energy consumption, the motion planning problem is transformed into multiobjective optimization problem. Select the running time (hi) between the via points and the joint angle θmi at the floating via points is used as an optimization variable. The optimization model is established as shown in (15), whereX is the optimized variables,SJ is the joint average acceleration which is used to measure the smoothness of the robot motion trajectory, andSq is the joint total energy consumption. The relevant parameters of the optimization model are set as follows: yMIN = −1 m, yMAX = 1 m, zMIN = 0.6 m, zMAX = 1 m, δMAX = 1 mm, θ˙mLB= −30°/s, θ˙mUB = 30°/s, θ¨mLB= −150°/s2, θ¨mUB = 150°/s2, JMAX = 500°/s3, τUB = 8000 Nm, and MUP = 100.(15)minXH=∑i=1nhiSJ=∑m=1M1H∫0Hθ⃛mit2dtSq=∑i=0n∑m=1M∫hihi+1τmitdθmis.t.yMIN≤yi∗≤yMAX;zMIN≤zi∗≤zMAX;yi∗-ySTRAIGHT≤δMAX;θ˙mLB≤θ˙mit≤θ˙mUB;θ¨mLB≤θ¨mit≤θ¨mUB;θ⃛mit≤JMAX;τ≤τUB;∑m=1M1/H∫0Hτmit2dt1/N∑k=1NτmikΔT-τmik-1ΔT2≤MUP.


### 3.3. Analysis of Optimization Results

The unstacking robot is recycled in accordance with the standard trajectory before the optimization, the total time-consuming is 3.5 seconds, the robot joint trajectory is constructed with 7 order B-spline curve, and the velocity, acceleration, and jerk are set to zero at the start and stop times. For multiobjective optimization problem, genetic algorithm is one of the most mature methods, and this paper chose NSGA II multiobjective evolutionary algorithm [19]. According to the actual experience, the number of population is 200, the number of iterations is 1000, the probability of parent crossing is 0.5, and the probability of individual population variation is 0.01. The Pareto frontier of the optimization model is shown in Figure 5.Figure 5
The Pareto frontier of the optimization model.It can be seen in Figure5 that robot efficiency index and track smoothness indicators and energy consumption indicators conflict with each other; the higher the efficiency, the greater the energy consumption and the less smooth the trajectory. In order to shorten the no-load running time of the robot during recycling process, take the efficiency as the main objective and energy consumption and smoothness as reference index. To get the best results from the Pareto frontier surface of the optimization model, add constraints to the reference index as SJ≤ 800°/s2 and Sq≤3000 J. So the optimal results can be obtained as shown in Table 2.Table 2
Trajectory optimization results.
Joint 1
Joint 2
Joint 3
Optimal position sequence (°)
[−45  −45  −45  −45  −45  −45  −44.852  −42.154  −25.503   −16.098  −0.967   16.616  32.361  40.608  44.449  45   45  45  45  45  45]
[55.2078  55.253  55.803  59.071  61.206  63.373  66.314  76.713  91.353   93.209  94.417  92.086  87.013  79.314   68.553  62.888  59.210  55.221  51.230   50.733  50.644]
[116.812  116.827  117.050  118.214  119.520  121.034  123.033  129.159   133.419  133.504  133.464  133.490   131.504  125.722  119.897  118.047   116.636  114.876  112.666  112.302   112.314]
Optimal time series (s)
[0  0.028  0.084  0.214  0.273  0.322  0.378  0.543  0.844  0.942  1.106  1.320  1.516  1.691  1.887  1.982  2.048  2.143  2.332  2.423  2.479]
Optimal object
H
=
2.479 s, SJ=794.14°/s2, Sq=2785.64 JAccording to the optimization results, the Cartesian spatial trajectory comparison and the corresponding optimal trajectory in joint space are obtained in Figure6. The velocity, acceleration, and output torque of each joint are shown in Figure 7.Figure 6
Cartesian spatial trajectory comparison and the corresponding optimal trajectory in joint space.
(a)
Cartesian spatial trajectory comparison
(b)
Optimal trajectory in joint spaceFigure 7
Comparison of each joint before and after optimization.
(a)
Comparison of the velocity of each joint before optimization
(b)
Comparison of the velocity of each joint after optimization
(c)
Comparison of the acceleration of each joint before optimization
(d)
Comparison of the acceleration of each joint after optimization
(e)
Comparison of the output torque of each joint before optimization
(f)
Comparison of the output torque of each joint after optimizationIt can be seen in Figure6 that time required for the robot recycling process is reduced from 3.5 s to 2.479 s, and the production efficiency is improved by 29.17%. The total energy consumption is reduced from 3186.15 J to 2785.64 J, and the joint average acceleration is reduced from 943.71°/s2 to 794.14°/s2; that is, the robot energy consumption and track smoothness were reduced by 12.57% and 15.85%. It can be seen from Figures 7(a)–7(d) that although the time required for the recycling process is greatly shortened by optimization, the angular velocity and acceleration of each joint are not significantly increased, and the acceleration of the starting point and the ending point are zero, so as to avoid causing infinite jerk and the impact on the system. It can be seen in Figures 7(e)–7(f) that the joint torque values are within the allowable range after optimization, especially when the robot end effector is inserted into the rubber. Joints 2 and 3 are subjected to significantly impact torque and the resistance is large in the process before optimization. Through the optimization of the floating via points, the time of end effector inserted into the rubber process is shortened, and at the same time the impact force and resistance during the insertion process are significantly reduced.


## 3.1. B-Spline Trajectory Planning in Joint Space

Robot joint space trajectory planning is based on a series of key constraint points to plan robot motion trajectory which areθmi=θmiti+1,i=0,1,…,n,j=1,2,…,M, where θmi(t) represents theith trajectory segment of joint m. θmi is the angle at the i+1 via point of jointm that is the angle at ti+1, andM is the number of robot joints. The number of constraint points for any joint trajectory curve is n+1, which are starting point, end point, and n-1 via points, respectively.The unstacking robot joint space trajectory is described using the nonuniform B-spline curve. Based on determining the orderk and the node vectorU= u0,u1,…,un+k+1, B-spline curve recursive equation using the de Boer algorithm is achieved as follows [18]:(6)pu=∑j=i-k+1idjlNj,k-1u=⋯=dik,ui≤u<ui+1,djl=dj,l=01-αjldj-1l-1+αjldjl-1,l=1,2,…,k;j=i-k+l,…,i,αjl=u-ujuj+k+1-l-uj,where di∈RM×1 is the control vertices of the joint spline trajectory, u∈ui,ui+1 is the normalized time node vector of spline curve in segmenti, and Ni,k(u) is k-order specification B-spline basis function. According to the local support character of the B-spline curve, ther-order derivative of p(u) is obtained(7)pru=∑j=i-k+ridjrNj,k-ru,ui≤u<ui+1,djl=dj,l=0k+1-ldjl-1-dj-1l-1uj+k+1-l-uj,l=1,2,…,r;j=i-k+l,…,i.According to (7), robot joint space velocity, acceleration, and jerk can be obtained as follows:(8)θ˙mit=p′u=∑j=i-k+1idj1Nj,k-1u,θ¨mit=p′′u=∑j=i-k+2idj2Nj,k-2u,θ⃛mit=p′′′u=∑j=i-k+3idj3Nj,k-3u.


## 3.2. Optimization Model

During the robot recycling process, the trajectories of each joint are designed as segment B-spline curves. By using the same time interval parameterhi to couple the joint tracks to each other, whereθmi(t) represents theith trajectory of jointm, θmi, θ˙mi, and θ¨mi are, respectively, the angle, angle velocity, and acceleration at thei via point of jointm. The main constraints of the optimization process are set as follows.(
1) Position Constraint. According to the recycling characteristics, the rubber unstacking robot usually adopts the “door-shaped” trajectory to operate as shown in Figure 4. In order to avoid interference, the AB, CE, and GE segments are set as linear trajectories. In some former researches [1–7], the trajectory is usually optimized according to a number of fixed teaching points and the way of inserting fixed via points. The optimization method overreliance on the selection of the initial trajectory is mainly to optimize the segmentation trajectory between the fixed via points. However, since all the via points have been set in advance, the optimization range of the Cartesian space trajectory is limited and it is difficult to obtain the optimal trajectory.Figure 4
Recycling process trajectory before optimization.The trajectory optimization method proposed in this paper contains two movement forms: linear trajectories AB and GH which are limited by obstacle avoidance and recycling form; the free trajectory passes through the barrier point B, the highest point D, and the insertion point G which are optimized by the rectangular space composed of B, D, and G points; meanwhile, to avoid interference, it is necessary to ensure that the robot end effector is always moving in the space. The positional coordinates of the end effector in the Cartesian space along they- and z-axes can be solved by positive kinematics. It is only necessary to verify whether the limit position of each trajectory satisfies the constraint condition as compared with the coordinate value of the whole motion time period. When the velocity of the rubber unstacking robot along y- and z-axes is zero, the corresponding position is the limit position of the B-spline segment, which can be calculated as follows:(9)y˙it=cos⁡θ1itl1cos⁡θ2it-l3cos⁡θ3it+S1θ˙1it-l1sin⁡θ1itsin⁡θ2itθ˙2it+l3sin⁡θ1itsin⁡θ3itθ˙3it,z˙it=l1cos⁡θ2itθ˙2it-l3cos⁡θ3itθ˙3it.The limit position ofyi∗ and zi∗ of the spline trajectory of theith segment is obtained by the above equations, and the spatial position constraints are as follows:(10)yMIN≤yi∗≤yMAX,i=b,…,g,zMIN≤zi∗≤zMAX,i=b,…,g,whereb and g are the spline segments of the transition point B and G, respectively. Since the linear trajectory segments AB and GE are in the z-axis direction, the position constraint is as follows:(11)yi∗-ySTRAIGHT≤δMAX,i=1,…,b&g,…,n.(
2) Velocity, Acceleration, and Jerk Constraints. Under the condition of the drive system which has been set, the joint velocity and acceleration of the rubber unstacking robot are constrained to a certain extreme range and, at the same time, due to excessive acceleration, can cause the robot mechanical system resonance and increase the wear of the mechanical system, so the jerk is also used as one of the constraints. Any section trajectory of any joint should satisfy the constraint conditions as follows:(12)θ˙mLB≤θ˙mit≤θ˙mUB,θ¨mLB≤θ¨mit≤θ¨mUB,θ⃛mit≤JMAX,i=0,…,n,j=1,…,m,where θ˙mLB, θ˙mUB, θ¨mLB, and θ¨mUB are, respectively, the velocity and acceleration limit value of jointm and JMAX is the maximum jerk. In the process of optimization, the constraint will be transformed into a constraint on the first three derivatives of the B-spline control vertex dm as shown in (7).(
3) Torque Constraints. The joint output torque is limited by the robot drive system and needs to meet the corresponding constraints. According to (1) one can get(13)τ=Mθθ¨+Cθ,θ.θ.+Gθ+JTθFs≤τUB,where τUB is the output torque limit, which is defined as the 80% of the maximum torque of the joint motor.(
4) Normalized Torque Constraints. In order to ensure the smoothness of the robot motion, the normalized torque is taken as the constraint condition as follows:(14)∑m=1M1/H∫0Hτmit2dt1/N∑k=1NτmikΔT-τmik-1ΔT2≤MUP,whereN is the number of servo control cycles, ΔT is the robot servo cycles, and MUP is the maximum value.In order to improve the working efficiency, take the work time of the rubber unstacking robot during recycling process as the optimization goal. At the same time, taking into account the trajectory smoothness and energy consumption, the motion planning problem is transformed into multiobjective optimization problem. Select the running time (hi) between the via points and the joint angle θmi at the floating via points is used as an optimization variable. The optimization model is established as shown in (15), whereX is the optimized variables,SJ is the joint average acceleration which is used to measure the smoothness of the robot motion trajectory, andSq is the joint total energy consumption. The relevant parameters of the optimization model are set as follows: yMIN = −1 m, yMAX = 1 m, zMIN = 0.6 m, zMAX = 1 m, δMAX = 1 mm, θ˙mLB= −30°/s, θ˙mUB = 30°/s, θ¨mLB= −150°/s2, θ¨mUB = 150°/s2, JMAX = 500°/s3, τUB = 8000 Nm, and MUP = 100.(15)minXH=∑i=1nhiSJ=∑m=1M1H∫0Hθ⃛mit2dtSq=∑i=0n∑m=1M∫hihi+1τmitdθmis.t.yMIN≤yi∗≤yMAX;zMIN≤zi∗≤zMAX;yi∗-ySTRAIGHT≤δMAX;θ˙mLB≤θ˙mit≤θ˙mUB;θ¨mLB≤θ¨mit≤θ¨mUB;θ⃛mit≤JMAX;τ≤τUB;∑m=1M1/H∫0Hτmit2dt1/N∑k=1NτmikΔT-τmik-1ΔT2≤MUP.


## 3.3. Analysis of Optimization Results

The unstacking robot is recycled in accordance with the standard trajectory before the optimization, the total time-consuming is 3.5 seconds, the robot joint trajectory is constructed with 7 order B-spline curve, and the velocity, acceleration, and jerk are set to zero at the start and stop times. For multiobjective optimization problem, genetic algorithm is one of the most mature methods, and this paper chose NSGA II multiobjective evolutionary algorithm [19]. According to the actual experience, the number of population is 200, the number of iterations is 1000, the probability of parent crossing is 0.5, and the probability of individual population variation is 0.01. The Pareto frontier of the optimization model is shown in Figure 5.Figure 5
The Pareto frontier of the optimization model.It can be seen in Figure5 that robot efficiency index and track smoothness indicators and energy consumption indicators conflict with each other; the higher the efficiency, the greater the energy consumption and the less smooth the trajectory. In order to shorten the no-load running time of the robot during recycling process, take the efficiency as the main objective and energy consumption and smoothness as reference index. To get the best results from the Pareto frontier surface of the optimization model, add constraints to the reference index as SJ≤ 800°/s2 and Sq≤3000 J. So the optimal results can be obtained as shown in Table 2.Table 2
Trajectory optimization results.
Joint 1
Joint 2
Joint 3
Optimal position sequence (°)
[−45  −45  −45  −45  −45  −45  −44.852  −42.154  −25.503   −16.098  −0.967   16.616  32.361  40.608  44.449  45   45  45  45  45  45]
[55.2078  55.253  55.803  59.071  61.206  63.373  66.314  76.713  91.353   93.209  94.417  92.086  87.013  79.314   68.553  62.888  59.210  55.221  51.230   50.733  50.644]
[116.812  116.827  117.050  118.214  119.520  121.034  123.033  129.159   133.419  133.504  133.464  133.490   131.504  125.722  119.897  118.047   116.636  114.876  112.666  112.302   112.314]
Optimal time series (s)
[0  0.028  0.084  0.214  0.273  0.322  0.378  0.543  0.844  0.942  1.106  1.320  1.516  1.691  1.887  1.982  2.048  2.143  2.332  2.423  2.479]
Optimal object
H
=
2.479 s, SJ=794.14°/s2, Sq=2785.64 JAccording to the optimization results, the Cartesian spatial trajectory comparison and the corresponding optimal trajectory in joint space are obtained in Figure6. The velocity, acceleration, and output torque of each joint are shown in Figure 7.Figure 6
Cartesian spatial trajectory comparison and the corresponding optimal trajectory in joint space.
(a)
Cartesian spatial trajectory comparison
(b)
Optimal trajectory in joint spaceFigure 7
Comparison of each joint before and after optimization.
(a)
Comparison of the velocity of each joint before optimization
(b)
Comparison of the velocity of each joint after optimization
(c)
Comparison of the acceleration of each joint before optimization
(d)
Comparison of the acceleration of each joint after optimization
(e)
Comparison of the output torque of each joint before optimization
(f)
Comparison of the output torque of each joint after optimizationIt can be seen in Figure6 that time required for the robot recycling process is reduced from 3.5 s to 2.479 s, and the production efficiency is improved by 29.17%. The total energy consumption is reduced from 3186.15 J to 2785.64 J, and the joint average acceleration is reduced from 943.71°/s2 to 794.14°/s2; that is, the robot energy consumption and track smoothness were reduced by 12.57% and 15.85%. It can be seen from Figures 7(a)–7(d) that although the time required for the recycling process is greatly shortened by optimization, the angular velocity and acceleration of each joint are not significantly increased, and the acceleration of the starting point and the ending point are zero, so as to avoid causing infinite jerk and the impact on the system. It can be seen in Figures 7(e)–7(f) that the joint torque values are within the allowable range after optimization, especially when the robot end effector is inserted into the rubber. Joints 2 and 3 are subjected to significantly impact torque and the resistance is large in the process before optimization. Through the optimization of the floating via points, the time of end effector inserted into the rubber process is shortened, and at the same time the impact force and resistance during the insertion process are significantly reduced.


## 4. Adaptive Backstepping Control Strategy Based on RFWNN Method

The adaptive backstepping control is a systematic and recursive design methodology for the feedback control of nonlinear systems with parametric uncertainties. Unlike the feedback linearization method with the problems such as the precise model requirement and the cancellation of useful nonlinear terms, the backstepping approach offers a choice of design tools for accommodation of uncertainties and nonlinearities and can avoid wasteful cancellations.


### 4.1. Backstepping Control System Planning

The key idea of the backstepping design is to select recursively some appropriate state variables as virtual inputs for lower dimension subsystems of the overall system and the Lyapunov functions are designed for each stable virtual controller. Therefore, the finally designed actual control law can guarantee the stability of total control system. Combined with (1), the dynamic model of the rubber unstacking robot during recycling process which considered the robot joint friction and the system unknown disturbance is shown as follows:(16)Mθθ¨+Cθ,θ˙θ˙+Gθ+fθ˙+JTθFs+τd=τ,where f(θ˙)∈RM×1 is the joint friction term and τd∈RM×1 is the system unknown disturbance term.θ
d
(
t
), θ˙d(t), and θ¨d(t) are the optimal trajectory parameters calculated in Section 3. Define M-θ, C-(θ,θ˙), G-(θ), and f-(θ˙) as the actual values with uncertainties in the nominal values Mθ, C(θ,θ˙), G(θ), and f(θ˙), respectively. According to the above definition, the robot dynamics model is reconstructed as follows:(17)Mθθ¨+yθ,θ˙+JTθFs+τd=τ,where yθ,θ˙ is the uncertainty term of the robot dynamics system which can be expressed as follows:(18)yθ,θ˙=MθM--1θC-θ,θ˙θ˙+G-θ+fθ˙.Define state variablesX1=θ and X2=θ˙; dynamics model (17) is reconstructed as follows:(19)X˙1=X2,X˙2=M-1X1τ-JTθFs-yθ,θ˙-τd.According to (19), the goal of the control system is to design an adaptive control law so that the state vectorX1 can precisely track the desired trajectory θd. Define tracking error E1(t) as(20)E1t=X1t-θdtand its derivative is(21)E˙1t=X˙1t-θ˙dt=vt-θ˙dt,where vt=X˙1t is called the virtual control function and the stabilizing function st is defined as(22)st=-ς1E1t+θ˙dt,whereς1 is a positive definite diagonal matrix. The first Lyapunov function is chosen as(23)V1t=12E1TtE1tand its derivative is(24)V˙1t=E1TtE˙1t=E1TtX˙1t-θ˙dt=E1Ttvt-st-ςE1t.According to the Lyapunov stability criterion, the position tracking error E1t converges asymptotically when the virtual control function v(t) is equal to the stability functions(t). Define E2t=v(t)-s(t)-E˙1t, and then the derivative of E2t is expressed as(25)E˙2t=v˙t-s˙t=X˙2t+ς1E˙1t-θ¨dt=M-1X1τ-JTX1Fs-yX1,X2-τd+ς1E˙1t-θ¨t. To design the backstepping control system, the Lyapunov function is defined as(26)V2t=V1t+12E2TtE2t.And its derivative can be derived as follows:(27)V˙2t=V˙1t+E2TtE˙2t=E1TtE2t-ς1E1t+E2Tt×M-1X1τ-JTX1Fs-yX1,X2-τd+ς1E˙1t-θ¨t.From (27), the backstepping control law that satisfies the Lyapunov stability condition is designed as follows:(28)τ=MX1θ¨t-Λt+JTX1Fs+yX1,X2+τd,where Λt=ς1E˙1t+ς2E2t+E1t is the defined system synthesis error andς2 is a positive definite diagonal matrix. The backstepping control law shown in (28) is difficult to achieve due to the effect of the uncertainty term of the robot dynamic model y(X1,X2) and the time-varying disturbance τd. Therefore, adaptive robust control technology is needed, which can not only identify the uncertain dynamic of the robot but also compensate for the external time-varying disturbance to ensure the tracking accuracy of the robot trajectory. An adaptive robust online learning control strategy based on recursive fuzzy wavelet neural network (RFWNN) is proposed. RFWNN is used to estimate the uncertainty term of the robot dynamics system. The robust compensator compensates for the external disturbance and the estimation error of the former.


### 4.2. RFWNN Structure

The RFWNN structure consists of a recursive structure and an FWNN structure as shown in Figure8, FWNN structure also includes TSK fuzzy system based on fuzzy rules and wavelet neural network (WNN) structure. Each fuzzy rule is considered as a sub-WNN space, and each sub-WNN contains a wavelet with a specified dilation value (resolution). Therefore, sub-WNNs of different resolution levels can capture different behaviors of the approximate functions, and the role of fuzzy sets is to determine the contribution of the sub-FWNN space to input/output of the FWNNs. Therefore, with the main feature stated above, the RFWNN can achieve the same approximation accuracy with a smaller size in comparison with the FNNs/NNs. In addition, feedback loop with the time delay at the FWNN fuzzy layer will make the system become a dynamic structure that have sufficient capability to handle the high dynamics of the robot control system.Figure 8
RFWNN structure.The RFWNN structure consists of four layers, the first layer is the input layer, input signalsx1,x2,…,xn are transported to the fuzzy layer by neurons, each of which is represented by its corresponding fuzzy membership function μAij(xi), and the equation is as follows:(29)μAijxi=e-dji2xi-cji2,where dji is scale parameters, cji is translation parameter,p is rule number, and j=1,…,p,i=1,…,n,(p,n)∈N. A local feedback unit with a real-time delay section is added to the layer, so the recursive input of the layer is represented as follows:(30)xrit=xit+πriμAijxit-T,where the membership function μAijxit-T represents the delay value of μAij(xi) within the time intervalT and πri is the recursive weight vector of the feedback unit.The third layer is a fuzzy rule layer, each neuron in the layer is represented as a rule, the prerequisite for the execution of the neuron is matched with the relevant rule, and the output of the layer is multiplied by(31)ωj=∏iwAijμAijxri,where wAij is the weight between the fuzzy layer and the rule layer and the multiplication operator ∏ is used to determine the excitation intensity in the fuzzy reasoning mechanism. The excitation intensity of rulej is simplified by a combination of ωj and ∏iζji(xri). The fuzzy wavelet basis function is as follows:(32)Φjxr=ωj∏iζjixri,where ζji(xri)=1-dji2xri-cji2, xri=xr1,…,xrnT∈Rpn×1, and r=1,…,p. The wavelet basis function Φj(xr) is expressed as a multidimensional wavelet function as follows:(33)Φjxr=∏iϕjxri,where ϕj(xri)=1-dji2xri-cji2e-dji2xri-cji2 is the Mexican Cap Wavelet function.The fourth layer is a fuzzy output layer, where each node represents the output variable and is calculated by summing all the input signals(34)yl=∑jwljΦjxr,where wlj is the weight between the rule layer and the output layer, l=1,…,m,m∈N,m is RFWNN number of output variables, and the above equation is expressed as a vector form as follows:(35)yx,d,c,π,W=WTΦxr,d,c,π;Φxr,d,c,π=Φ1xr,Φ2xr,…,ΦpxrT∈Rp×1;WT=w11w12⋯w1pw21w22⋯w2p⋮⋮⋱⋮wm1wm2⋯wmp∈Rm×p,where d∈Rnp×1, x∈Rn×1, and y∈Rm×1 are scaling parameter vector, c∈Rnp×1 is translation parameter vector, and π∈Rnp×1 is recursive weight parameter vector.The RFWNN structure is mainly used to online estimate the uncertain dynamics model of the robot. Based on the approximate error analysis, there is an ideal RFWNN structure with the optimal parameters [20](36)yxt=W∗TΦ∗xrt,d∗,c∗,π∗+Δxrt,where W∗, d∗, c∗, and π∗ are the optimal parameters corresponding to W, d, c, and π and Δ(xr(t)) is the estimation error vector.


### 4.3. Adaptive Backstepping Controller Design

In (28), the y part is approximated by the RFWNN network, the process will inevitably produce approximate error, and the τd part contains the unstructured uncertainties interference of the mechanical model established by the robot insertion process. Thus, this section estimates and compensates for the approximate error and τd by adding robust terms to the controller. System torque control law is(37)τ=Mθ¨-Λ+JTFs+y^+τ^d,where τ^d is a robust term used to compensate for uncertainties, such as approximate error, unknown interference, and unspecified parts of the robot control system, and y^ is an approximate term for the equation y, expressed as(38)y^xt=W^TΦxrt,d^,c^,π^,where x=θT,θ˙T,θdT,θ˙dT,θ¨dTT and W^, d^, c^, and π^ are, respectively, the approximate value of W∗, d∗, c∗, and π∗.The trajectory tracking control rate shown in (37) is introduced into the robot dynamics equation shown in (17) to obtain the closed-loop control system of the rubber unstacking robot(39)-MΛ=τd+y~-τ^d,where y~ is as follows:(40)y~=W∗TΦ~+W~TΦ^+Δ,where W~=W∗-W^, Φ^≡Φ(xr,d^,c^,π^), Φ∗≡Φxr,d∗,c∗,π∗, and Φ~=Φ∗-Φ.In order to achieve good tracking performance of the controller, the parameters such as weight, scaling, and translation in the RFWNN neural network need to be adjusted online. The nonlinear output of the RFWNN is transformed into a partial linear form by using the linearization technique, which can be extended by applying the Lyapunov theorem. Thus,Φ~ is expanded in the Taylor series to get the following form:(41)Φ~=PTd∗-d^+KTc∗-c^+HTπ∗-π^+Θd∗-d^,c∗-c^,π∗-π^,where vector Θ is the higher order term in the Taylor series expansion, P=∂Φ1/∂d,…,∂Φp/∂dd=d^, K=∂Φ1/∂c,…,∂Φp/∂cc=c^, H=∂Φ1/∂π,…,∂Φp/∂ππ=π^ are bounded normal vectors, and(42)∂Φj∂d=0,…,0︸j-1n,∂Φj∂dj1,…,∂Φj∂djn,0,…,0︸p-jnT,∂Φj∂c=0,…,0︸j-1n,∂Φj∂cj1,…,∂Φj∂cjn,0,…,0︸p-jnT,∂Φj∂π=0,…,0︸j-1n,∂Φj∂πj1,…,∂Φj∂πjn,0,…,0︸p-jnT.Defined~=d∗-d^, c~=c∗-c^, and π~=π∗-π^, and (41) is as follows:(43)Φ~=PTd~+KTc~+HTπ~+Θd~,c~,π~.According to (40) and (43) one can obtain the following:(44)y~+τd=W~TΦ^-PTd^-KTc^-HTπ^+W^TPTd~+KTc~+HTπ~+Γu,where Γu=W∗TΦ~+PTd^+KTc^+HTπ^-W^TPTd∗+KTc∗+HTπ∗+Δ+τd.Substituting (44) into (39),(45)-MΛ=W~TΦ^-PTd^-KTc^-HTπ^-τ^d+Γu+W^TPTd~+KTc~+HTπ~.Set the positive real numberWB, dB, cB, πB, ΔB, and τdB to satisfy the condition W∗≤WB, d∗≤dB, c∗≤cB, π∗≤πB, Δ≤ΔB, and τd≤τdB, and it can be obtained as follows:(46)Γu=W∗TΦ~+ITd^+KTc^+HTπ^-W^TITd∗+KTc∗+HTπ∗+Δ+τd≤W∗TΦ~+Δ+τd+W∗TITd^+W∗TKTc^+W∗THTπ^+ITd∗+KTc∗+HTπ∗W^.AddαWB2/4+αdB2/4+αcB2/4+απB2/4 to both ends of inequality (46),(47)Γu+αWB24+αdB24+αcB24+απB24≤β∗Tυ,where υT=1,d^,c^,π^,W^,α is positive constant, β∗T=β1,β2,β3,β4,β5,β1,β2,β3,β4, andβ5 are positive constant and the uncertainty of the following term: W∗TΦ~+Δ+τd+αWB2/4+αdB2/4+αcB2/4+απB2/4,  W∗TPT,  W∗TKT,  W∗THT, and PTd∗+KTc∗+HTπ∗.In order to ensure that the dynamic closed-loop system shown in (35) is stable, the robust term τ^d must cancel the indeterminate term Гu. Therefore, using τ^d to estimate the uncertainty limit β∗Tυ, the expression of the robust term τ^d is given as follows:(48)τ^d=E1E1β^Tυ+bsE1E12+δ,where bs and δ are positive constants and β^ is the online estimate of β∗. According to the above analysis, an adaptive online learning algorithm based on RFWNN control system and a robust term expression are proposed as follows:(49)W^˙=KWΦ^-PTd^-KTc^-HTπ^E1T-αKWE1W^,d^˙=KdPW^E1-αKdE1d^;c^˙=KcKW^E1-αKcE1c^,π^˙=KπHW^E1-αKπE1π^;β^˙=E1Kβυ,where KW, Kd, Kc, Kπ, and Kβ are diagonal normal matrix. Based on the above analysis, the controller structure is shown in Figure 9.Figure 9
Control block diagram of the adaptive backstepping controller.Combined with the control block diagram, the design procedure of the proposed control system is detailed in the following:(1)
Specify the structure of the RFWNN(n,p,m), and set constant parameters (ς1,ς2,KW,Kd,Kc,Kπ,Kβ and α).(2)
Initial the parameters of the RFWNNW^,c^,d^,π^,μAijxi0.(3)
Update the RFWNN inputs (x=θT,θ˙T,θdT,θ˙dT,θ¨dTT) and memory data W^,c^,d^,π^,μAijxit-T.(4)
Compute tracking errorE1(t) and system synthesis error Λ(t), compute the recurrent inputs xri(t) via (30), and compute the membership functions μAij(xri(t)) via (29).(5)
Compute the outputsωj of the fuzzy rule layer via (31), and then compute fuzzy wavelet basic functions Φj(xr) via (32) and (33).(6)
Compute the outputs of the RFWNNy^(x(t)) via (34), compute the robust term τ^d via (48), and then construct the control inputτ via (37).(7)
Adjust the weights, dilations, and translations and the recurrent values of the RFWNNs via (49).(8)
Save these data in the memory and return to step (3).


## 4.1. Backstepping Control System Planning

The key idea of the backstepping design is to select recursively some appropriate state variables as virtual inputs for lower dimension subsystems of the overall system and the Lyapunov functions are designed for each stable virtual controller. Therefore, the finally designed actual control law can guarantee the stability of total control system. Combined with (1), the dynamic model of the rubber unstacking robot during recycling process which considered the robot joint friction and the system unknown disturbance is shown as follows:(16)Mθθ¨+Cθ,θ˙θ˙+Gθ+fθ˙+JTθFs+τd=τ,where f(θ˙)∈RM×1 is the joint friction term and τd∈RM×1 is the system unknown disturbance term.θ
d
(
t
), θ˙d(t), and θ¨d(t) are the optimal trajectory parameters calculated in Section 3. Define M-θ, C-(θ,θ˙), G-(θ), and f-(θ˙) as the actual values with uncertainties in the nominal values Mθ, C(θ,θ˙), G(θ), and f(θ˙), respectively. According to the above definition, the robot dynamics model is reconstructed as follows:(17)Mθθ¨+yθ,θ˙+JTθFs+τd=τ,where yθ,θ˙ is the uncertainty term of the robot dynamics system which can be expressed as follows:(18)yθ,θ˙=MθM--1θC-θ,θ˙θ˙+G-θ+fθ˙.Define state variablesX1=θ and X2=θ˙; dynamics model (17) is reconstructed as follows:(19)X˙1=X2,X˙2=M-1X1τ-JTθFs-yθ,θ˙-τd.According to (19), the goal of the control system is to design an adaptive control law so that the state vectorX1 can precisely track the desired trajectory θd. Define tracking error E1(t) as(20)E1t=X1t-θdtand its derivative is(21)E˙1t=X˙1t-θ˙dt=vt-θ˙dt,where vt=X˙1t is called the virtual control function and the stabilizing function st is defined as(22)st=-ς1E1t+θ˙dt,whereς1 is a positive definite diagonal matrix. The first Lyapunov function is chosen as(23)V1t=12E1TtE1tand its derivative is(24)V˙1t=E1TtE˙1t=E1TtX˙1t-θ˙dt=E1Ttvt-st-ςE1t.According to the Lyapunov stability criterion, the position tracking error E1t converges asymptotically when the virtual control function v(t) is equal to the stability functions(t). Define E2t=v(t)-s(t)-E˙1t, and then the derivative of E2t is expressed as(25)E˙2t=v˙t-s˙t=X˙2t+ς1E˙1t-θ¨dt=M-1X1τ-JTX1Fs-yX1,X2-τd+ς1E˙1t-θ¨t. To design the backstepping control system, the Lyapunov function is defined as(26)V2t=V1t+12E2TtE2t.And its derivative can be derived as follows:(27)V˙2t=V˙1t+E2TtE˙2t=E1TtE2t-ς1E1t+E2Tt×M-1X1τ-JTX1Fs-yX1,X2-τd+ς1E˙1t-θ¨t.From (27), the backstepping control law that satisfies the Lyapunov stability condition is designed as follows:(28)τ=MX1θ¨t-Λt+JTX1Fs+yX1,X2+τd,where Λt=ς1E˙1t+ς2E2t+E1t is the defined system synthesis error andς2 is a positive definite diagonal matrix. The backstepping control law shown in (28) is difficult to achieve due to the effect of the uncertainty term of the robot dynamic model y(X1,X2) and the time-varying disturbance τd. Therefore, adaptive robust control technology is needed, which can not only identify the uncertain dynamic of the robot but also compensate for the external time-varying disturbance to ensure the tracking accuracy of the robot trajectory. An adaptive robust online learning control strategy based on recursive fuzzy wavelet neural network (RFWNN) is proposed. RFWNN is used to estimate the uncertainty term of the robot dynamics system. The robust compensator compensates for the external disturbance and the estimation error of the former.


## 4.2. RFWNN Structure

The RFWNN structure consists of a recursive structure and an FWNN structure as shown in Figure8, FWNN structure also includes TSK fuzzy system based on fuzzy rules and wavelet neural network (WNN) structure. Each fuzzy rule is considered as a sub-WNN space, and each sub-WNN contains a wavelet with a specified dilation value (resolution). Therefore, sub-WNNs of different resolution levels can capture different behaviors of the approximate functions, and the role of fuzzy sets is to determine the contribution of the sub-FWNN space to input/output of the FWNNs. Therefore, with the main feature stated above, the RFWNN can achieve the same approximation accuracy with a smaller size in comparison with the FNNs/NNs. In addition, feedback loop with the time delay at the FWNN fuzzy layer will make the system become a dynamic structure that have sufficient capability to handle the high dynamics of the robot control system.Figure 8
RFWNN structure.The RFWNN structure consists of four layers, the first layer is the input layer, input signalsx1,x2,…,xn are transported to the fuzzy layer by neurons, each of which is represented by its corresponding fuzzy membership function μAij(xi), and the equation is as follows:(29)μAijxi=e-dji2xi-cji2,where dji is scale parameters, cji is translation parameter,p is rule number, and j=1,…,p,i=1,…,n,(p,n)∈N. A local feedback unit with a real-time delay section is added to the layer, so the recursive input of the layer is represented as follows:(30)xrit=xit+πriμAijxit-T,where the membership function μAijxit-T represents the delay value of μAij(xi) within the time intervalT and πri is the recursive weight vector of the feedback unit.The third layer is a fuzzy rule layer, each neuron in the layer is represented as a rule, the prerequisite for the execution of the neuron is matched with the relevant rule, and the output of the layer is multiplied by(31)ωj=∏iwAijμAijxri,where wAij is the weight between the fuzzy layer and the rule layer and the multiplication operator ∏ is used to determine the excitation intensity in the fuzzy reasoning mechanism. The excitation intensity of rulej is simplified by a combination of ωj and ∏iζji(xri). The fuzzy wavelet basis function is as follows:(32)Φjxr=ωj∏iζjixri,where ζji(xri)=1-dji2xri-cji2, xri=xr1,…,xrnT∈Rpn×1, and r=1,…,p. The wavelet basis function Φj(xr) is expressed as a multidimensional wavelet function as follows:(33)Φjxr=∏iϕjxri,where ϕj(xri)=1-dji2xri-cji2e-dji2xri-cji2 is the Mexican Cap Wavelet function.The fourth layer is a fuzzy output layer, where each node represents the output variable and is calculated by summing all the input signals(34)yl=∑jwljΦjxr,where wlj is the weight between the rule layer and the output layer, l=1,…,m,m∈N,m is RFWNN number of output variables, and the above equation is expressed as a vector form as follows:(35)yx,d,c,π,W=WTΦxr,d,c,π;Φxr,d,c,π=Φ1xr,Φ2xr,…,ΦpxrT∈Rp×1;WT=w11w12⋯w1pw21w22⋯w2p⋮⋮⋱⋮wm1wm2⋯wmp∈Rm×p,where d∈Rnp×1, x∈Rn×1, and y∈Rm×1 are scaling parameter vector, c∈Rnp×1 is translation parameter vector, and π∈Rnp×1 is recursive weight parameter vector.The RFWNN structure is mainly used to online estimate the uncertain dynamics model of the robot. Based on the approximate error analysis, there is an ideal RFWNN structure with the optimal parameters [20](36)yxt=W∗TΦ∗xrt,d∗,c∗,π∗+Δxrt,where W∗, d∗, c∗, and π∗ are the optimal parameters corresponding to W, d, c, and π and Δ(xr(t)) is the estimation error vector.


## 4.3. Adaptive Backstepping Controller Design

In (28), the y part is approximated by the RFWNN network, the process will inevitably produce approximate error, and the τd part contains the unstructured uncertainties interference of the mechanical model established by the robot insertion process. Thus, this section estimates and compensates for the approximate error and τd by adding robust terms to the controller. System torque control law is(37)τ=Mθ¨-Λ+JTFs+y^+τ^d,where τ^d is a robust term used to compensate for uncertainties, such as approximate error, unknown interference, and unspecified parts of the robot control system, and y^ is an approximate term for the equation y, expressed as(38)y^xt=W^TΦxrt,d^,c^,π^,where x=θT,θ˙T,θdT,θ˙dT,θ¨dTT and W^, d^, c^, and π^ are, respectively, the approximate value of W∗, d∗, c∗, and π∗.The trajectory tracking control rate shown in (37) is introduced into the robot dynamics equation shown in (17) to obtain the closed-loop control system of the rubber unstacking robot(39)-MΛ=τd+y~-τ^d,where y~ is as follows:(40)y~=W∗TΦ~+W~TΦ^+Δ,where W~=W∗-W^, Φ^≡Φ(xr,d^,c^,π^), Φ∗≡Φxr,d∗,c∗,π∗, and Φ~=Φ∗-Φ.In order to achieve good tracking performance of the controller, the parameters such as weight, scaling, and translation in the RFWNN neural network need to be adjusted online. The nonlinear output of the RFWNN is transformed into a partial linear form by using the linearization technique, which can be extended by applying the Lyapunov theorem. Thus,Φ~ is expanded in the Taylor series to get the following form:(41)Φ~=PTd∗-d^+KTc∗-c^+HTπ∗-π^+Θd∗-d^,c∗-c^,π∗-π^,where vector Θ is the higher order term in the Taylor series expansion, P=∂Φ1/∂d,…,∂Φp/∂dd=d^, K=∂Φ1/∂c,…,∂Φp/∂cc=c^, H=∂Φ1/∂π,…,∂Φp/∂ππ=π^ are bounded normal vectors, and(42)∂Φj∂d=0,…,0︸j-1n,∂Φj∂dj1,…,∂Φj∂djn,0,…,0︸p-jnT,∂Φj∂c=0,…,0︸j-1n,∂Φj∂cj1,…,∂Φj∂cjn,0,…,0︸p-jnT,∂Φj∂π=0,…,0︸j-1n,∂Φj∂πj1,…,∂Φj∂πjn,0,…,0︸p-jnT.Defined~=d∗-d^, c~=c∗-c^, and π~=π∗-π^, and (41) is as follows:(43)Φ~=PTd~+KTc~+HTπ~+Θd~,c~,π~.According to (40) and (43) one can obtain the following:(44)y~+τd=W~TΦ^-PTd^-KTc^-HTπ^+W^TPTd~+KTc~+HTπ~+Γu,where Γu=W∗TΦ~+PTd^+KTc^+HTπ^-W^TPTd∗+KTc∗+HTπ∗+Δ+τd.Substituting (44) into (39),(45)-MΛ=W~TΦ^-PTd^-KTc^-HTπ^-τ^d+Γu+W^TPTd~+KTc~+HTπ~.Set the positive real numberWB, dB, cB, πB, ΔB, and τdB to satisfy the condition W∗≤WB, d∗≤dB, c∗≤cB, π∗≤πB, Δ≤ΔB, and τd≤τdB, and it can be obtained as follows:(46)Γu=W∗TΦ~+ITd^+KTc^+HTπ^-W^TITd∗+KTc∗+HTπ∗+Δ+τd≤W∗TΦ~+Δ+τd+W∗TITd^+W∗TKTc^+W∗THTπ^+ITd∗+KTc∗+HTπ∗W^.AddαWB2/4+αdB2/4+αcB2/4+απB2/4 to both ends of inequality (46),(47)Γu+αWB24+αdB24+αcB24+απB24≤β∗Tυ,where υT=1,d^,c^,π^,W^,α is positive constant, β∗T=β1,β2,β3,β4,β5,β1,β2,β3,β4, andβ5 are positive constant and the uncertainty of the following term: W∗TΦ~+Δ+τd+αWB2/4+αdB2/4+αcB2/4+απB2/4,  W∗TPT,  W∗TKT,  W∗THT, and PTd∗+KTc∗+HTπ∗.In order to ensure that the dynamic closed-loop system shown in (35) is stable, the robust term τ^d must cancel the indeterminate term Гu. Therefore, using τ^d to estimate the uncertainty limit β∗Tυ, the expression of the robust term τ^d is given as follows:(48)τ^d=E1E1β^Tυ+bsE1E12+δ,where bs and δ are positive constants and β^ is the online estimate of β∗. According to the above analysis, an adaptive online learning algorithm based on RFWNN control system and a robust term expression are proposed as follows:(49)W^˙=KWΦ^-PTd^-KTc^-HTπ^E1T-αKWE1W^,d^˙=KdPW^E1-αKdE1d^;c^˙=KcKW^E1-αKcE1c^,π^˙=KπHW^E1-αKπE1π^;β^˙=E1Kβυ,where KW, Kd, Kc, Kπ, and Kβ are diagonal normal matrix. Based on the above analysis, the controller structure is shown in Figure 9.Figure 9
Control block diagram of the adaptive backstepping controller.Combined with the control block diagram, the design procedure of the proposed control system is detailed in the following:(1)
Specify the structure of the RFWNN(n,p,m), and set constant parameters (ς1,ς2,KW,Kd,Kc,Kπ,Kβ and α).(2)
Initial the parameters of the RFWNNW^,c^,d^,π^,μAijxi0.(3)
Update the RFWNN inputs (x=θT,θ˙T,θdT,θ˙dT,θ¨dTT) and memory data W^,c^,d^,π^,μAijxit-T.(4)
Compute tracking errorE1(t) and system synthesis error Λ(t), compute the recurrent inputs xri(t) via (30), and compute the membership functions μAij(xri(t)) via (29).(5)
Compute the outputsωj of the fuzzy rule layer via (31), and then compute fuzzy wavelet basic functions Φj(xr) via (32) and (33).(6)
Compute the outputs of the RFWNNy^(x(t)) via (34), compute the robust term τ^d via (48), and then construct the control inputτ via (37).(7)
Adjust the weights, dilations, and translations and the recurrent values of the RFWNNs via (49).(8)
Save these data in the memory and return to step (3).


## 5. Experimental Verification

A servo control system based on EtherCAT bus and 4-DOF unstacking robot body experimental platform were built in this work. Beckhoff IPC was used as industrial computer which included TwinCAT3 programming software. The RS2 series servo drives were applied in this work which supported EtherCAT bus communications. The actual output torque and actual position of the robot joints are fed back in real time by the NCtoPLC function of TwinCAT3 software. So real-time communication and robot control could be realized via high-speed EtherCAT bus. Figure10 shows the entire experimental system.Figure 10
Control system of unstacking robot.
(a)
Control system
(b)
Unstacking robotThe model of the unstacking robot is described by (16) with its main dynamics given as(50)Mθ=m11000m22m230m32m33,Cθ,θ˙=0C12θ˙2C13θ˙3C21θ˙10C23θ˙3C31θ˙1C32θ˙20,m11=4.4cos2θ3+0.26-45.1cos⁡θ3+0.26cos⁡θ2+234.4cos2θ3+345.9cos2θ2-394.6cos⁡θ2cos⁡θ3+192.5cos⁡θ2-94.6cos⁡θ3-10.8cos⁡θ3+0.26+52.6,m22=365.37,m33=242.26,m23=m32=22.6cos⁡θ2-θ3-0.26-194.2cos⁡θ2-θ3,C12=-2C21=45.1sin⁡θ2cos⁡θ3+0.26-345.9sin⁡θ2+394.6sin⁡θ2cos⁡θ3-192.5sin⁡θ3,C13=-2C31=10.8sin⁡θ3+0.26-4.4sin⁡2θ3+0.53+394.6cos⁡θ2sin⁡θ3+45.1sin⁡θ3+0.26cos⁡θ2+94.6sin⁡θ3-235.4sin⁡2θ3,C23=-C32=22.6sin⁡θ2-θ3-0.26-194.2sin⁡θ2-θ3.The RFWNN structure can be characterized byn=15, p=5, and m=3 nodes and the relevant parameters of the designed controller are set as follows:(51)σ=diag⁡120,60,25,bs=0.01,δ=0.001,Kc=Kd=Kπ=diag⁡35,KW=diag⁡60,Ks=diag⁡180,180,180,Kβ=diag⁡0.002,0.002,0.002,0.002,0.002.The control strategy is compared with PID control and PD synchronic control. The PID parameters are designed by a compromise between the superiority of control performance and the magnitude of control effort, and they are given asKp1=diag⁡(110,70,50), Ki1=diag⁡(0.6,0.4,0.5), and Kd1=diag⁡(1200,900,600). And the parameters of the PD synchronic controller are set as Kp2=diag⁡(120,50,80) and Kd2=diag⁡(1500,1200,800).In order to verify the effect of the viscoelastic contact between the robot and the rubber on the control performance of the system and the ability of the control strategy proposed in this paper to effectively overcome this effect, the experiment in this section is validated in two cases. In Case  1, the unstacking robot tracks only the optimal trajectory without physical contact with the rubber. In Case  2, the unstacking robot inserts the rubber with the optimal trajectory according to the actual operation requirements.In the first case, Figure11 represents the joint torques and tracking errors of the robot under PID control and PD synchronic control and the proposed control strategy. In the second case, the torque and tracking error of each robot joint under different control strategies are shown in Figure 12. Since the viscoelastic contact between the robot and the rubber occurs within the time 1.982 s ≤t≤ 2.479 s, the normalized mean square error (NMSE) of each joint trajectory is compared with the above two experimental cases during this time, as shown in Table 3.Table 3
NMSE comparison of each joint in different control strategies.
NMSE (×10−3)
Case  1 (rad)
Case  2 (rad)
Joint 1
Joint 2
Joint 3
Joint 1
Joint 2
Joint 3
PID
0.842
4.387
5.224
2.517
7.848
6.042
PD synchronic
0.955
2.282
2.897
1.338
3.635
4.055
Adaptive backstepping
0.225
1.512
2.199
0.301
2.167
2.552
N
M
S
E
=
∑
k
=
0
K
θ
k
+
1
-
θ
d
k
+
1
2
/
K and K is the total sampling instant.Figure 11
Experimental results of different control strategies in Case  1.
(a)
Joint torques under PID control
(b)
Joint tracking errors under PID control
(c)
Joint torques under PD sync control
(d)
Joint tracking errors under PD sync control
(e)
Joint torques under adaptive backstepping control
(f)
Joint tracking errors under adaptive backstepping controlFigure 12
Experimental results of different control strategies in Case  2.
(a)
Joint torques under PID control
(b)
Joint tracking errors under PID control
(c)
Joint torques under PD sync control
(d)
Joint tracking errors under PD sync control
(e)
Joint torques under adaptive backstepping control
(f)
Joint tracking errors under adaptive backstepping controlIt can be seen from Figure11 that because the robot does not directly contact with the rubber in Case  1, the robot runs smoothly without large torque ripple and the tracking error is within an acceptable range. Compared with PID control, both PD synchronic control and the proposed control strategy in this paper have better performance and can control the trajectory tracking error to a smaller extent.According to Figure12 and Table 3, it can be seen that, under the control of PID, the joint torque of the robot fluctuates greatly, and the impact is more intense especially during the robot contact with rubber in Case  2. In addition the trajectory tracking error increases significantly, which will seriously affect the normal operation of the robot. Although the PD synchronization control can relatively improve the tracking accuracy of the robot in Case  1, the controller parameters cannot be adjusted online adaptively, and the robustness to interference is not strong enough so that it still cannot meet the system requirements in Case  2. Compared with the above controllers, the RFWNN-based adaptive backstepping controller proposed in this paper can reduce the trajectory tracking error significantly when the robot is inserted into the rubber. In addition the torque peak and fluctuations of each joint are effectively suppressed, so the robot can precisely and accurately perform the rubber recycling process according to the optimal trajectory. According to Table 3, it can be concluded that compared with Case  1 the NMSE value increases significantly in Case  2. And compared with the PID and PD synchronization control, the control method proposed in this paper can decrease the NMSE value of the robot joints by more than 72% and 40%, respectively, in the impact process. In Case  2, the robust term of the control strategy and the online update adjustment of the estimated parameters are shown in Figure 13. It can be seen that the robust term τ^d and the estimated parameters W^, c^, π^, and d^ can be adaptively adjusted online according to the modeling error of the robot dynamics and the external disturbances to significantly enhance the system robustness.Figure 13
Robust term and online estimated parameters.
(a)
Robust termτ^d
(b)
Online estimated parameters


## 6. Conclusion

This work deals with multiobjective trajectory optimization and adaptive backstepping control strategy of the rubber unstacking robot during the recycling process; the main contents of the research include the following aspects:(1) Considering the mechanical analysis during the inserting process, the Hunt-Crossley nonlinear model based dynamic model of the rubber unstacking robot is deduced, which can be used for trajectory planning and controller design.(2) Aiming at the performance indicators of operational efficiency, trajectory smoothness, and energy consumption, the motion planning problem is transformed into multiobjective optimization problem. The position of the floating via points of the B-spline trajectory and the running time of each point are optimized, and the NSGA genetic algorithm is used to solve the problem. The optimization results show that compared to the standard trajectory the above indicators were raised by 29.17%, 15.85%, and 12.57%, respectively.(3) The proposed adaptive backstepping control method based on RFWNN is robust to external torque disturbances and can overcome the effects of the unknown dynamic parameters and approximation error. The backstepping strategy is designed based on Lyapunov synthesis, which can guarantee the stability of the control system. Through experimental analysis, compared with other control methods, the validity of the proposed control strategy is verified.(4) As the rubber in the heat refining process, adhesion phenomenon will be produced between the rubber blocks due to local melting, which would result in rubbers difficult to be separated. So the future work will be extended to deal with the control problem of the unknown viscoelastic force between rubber blocks during the rubber separation process.



---

*Source: 7974325-2018-01-02.xml*