Admittance-Based Adaptive Cooperative Control for Multiple Manipulators With Output Constraints

This paper proposes a novel adaptive control methodology based on the admittance model for multiple manipulators transporting a rigid object cooperatively along a predefined desired trajectory. First, an admittance model is creatively applied to generate reference trajectory online for each manipulator according to the desired path of the rigid object, which is the reference input of the controller. Then, an innovative integral barrier Lyapunov function is utilized to tackle the constraints due to the physical and environmental limits. Adaptive neural networks (NNs) are also employed to approximate the uncertainties of the manipulator dynamics. Different from the conventional NN approximation method, which is usually semiglobally uniformly ultimately bounded, a switching function is presented to guarantee the global stability of the closed loop. Finally, the simulation studies are conducted on planar two-link robot manipulators to validate the efficacy of the proposed approach.


I. INTRODUCTION
In recent decades, robots have been widely employed in a variety of applications, such as entertainment, manufacturing industry and medical service [1], [2], etc. Among these applications, robot manipulators play a significant role and hence has attracted a lot of attention. Subsequently, large amount of efforts have been dedicated to the advance of robot manipulator technology [3]- [8]. An increasing complexity of tasks performed by robots demand higher manipulation skills and has rendered single manipulators ineffective in many situations. Hence, the research about coordinated multiple manipulators is becoming progressively significant [9]- [11]. For instance, the problem of transporting a long, heavy object has been studied in [12], where the task was accomplished easily with the cooperation of multiple robots.
However, the cooperative control of multi-manipulators is much more complicated than that of single one. The complex dynamics due to the presence of a large number of degrees of freedom (DOFs) of multiple manipulators in a system, result in a closed holonomic constrained chain mechanism [13]. In [14], a compound position/force control strategy was investigated under the assumption that an object was carried without relative motion between the object and manipulators. However, the relative motion between manipulators during the cooperative movement exists for a number of applications. A self-tuning control for cooperative manipulators was studied in [15], where the closed kinematic chain was formed by manipulators in the presence of uncertainty of kinematics model. Two cascaded estimators were adopted to change the kinematic parameters online to attain a good tracking performance. A flexible payload transforming problem for multiple collaborative agents was addressed in [16], which was regarded as a formation control problem by modelling the contact forces as the gradients of nonlinear potentials.
For multiple manipulators transporting a rigid object, some properties of grasp need to be taken into consideration. In grasp planning, there are two main classes of grasps known as "form closure" and "force closure" grasps [17]. These terminologies were first used by Reuleaux [18] in 1875 when he investigated the mechanism of some machines. Form closure is a pure geometry property, which describes the capability that the contacts can prevent all motions of a grasped object [17]. Reuleaux found that to obtain the form closure property, at least four contact points are needed in a planar case and Somov reported the number for a general spatial case is at least seven [19]. A further detailed classification of form closure can be seen in [20]. The significant difference between force-closure and form-closure is whether the friction forces are considered or not [21]. The manner in which the contact forces are being exerted on the object and the kinematics of the manipulators are contained by Force closure.
Investigations of interaction control mainly include force control and impedance control. Impedance control was introduced in [22], with the kernel idea of taking the mechanical impedance into consideration and mapping from state of the system to the interaction forces. The feasible, robust impedence control can be seen in [13], [23], [24]. On the contrary, another control method called admittance can be regarded as an inverse process of impedance control, mapping from forces to desired trajectory [25]. In conventional admittance control system, with predefined desired trajectory and the interaction forces exerted on the object, a virtual desired trajectory can be obtained. In this paper, an admittance model with interaction force between manipulators and object (measured by pressure sensors) is utilized to generate reference trajectories online for each manipulator.
In practical applications, there are still considerable difficulties to determine the precise actual kinematic and dynamic models of robots due to the uncertainties and environmental disturbances. Therefore, neural networks based methods have been widely employed in designing adaptive controller in the absence of availability of a precise model of the system [26]- [32]. The radial basis function Neural Networks (RBFNN) is a highly effectively tool for obtaining a model due to its universal approximation ability, which means it can approximate a smooth nonlinear function, and hence has been extensively utilized in designing controller for uncertain or unknown system [6]. In [33], the NN approximation strategy was used to compensate for the uncertain dynamics of the manipulated object and the robot manipulator. A class of multiinput-multi-output (MIMO) nonlinear manipulators adaptive control problem was investigated in [34]. In this work, two RBFNN namely, a critic NN and an actor NN are employed to achive the optimal control. Control method based on adaptive RBFNN to guarantee the global stability was developed in [35]- [38]. To deal with the unknown system paramters and complex couplings among several subsystems, an ANN was presented in [39].
It is noted that the conventional NN approximation method can only guarantee the semi-globally uniformly ultimately bounded (SGUUB) stability of the closed-loop system, which means the NN approximator is only valid in its active domain. Therefore it is essential to design NN controller with globally uniformly ultimately bounded (GUUB) stability [40], [41]. In [41], a novel switching scheme was introduced by combining a robust controller with an adaptive neural controller in the approximation domain, to guaranteed that the states in closedloop system are GUUB. The actual control system always faces some limitations in practical situations, such as state constraints, input and output constraints, etc. Barrier Lyapunov Function (BLF) was proposed due to its advantage of achieving a good tracking performance and satisfying the constraints simultaneously [42]. By utilizing some barrier functions whose value is infinite at corresponding limits, these methodologies keep the barriers not to be violated. As a consequence, the constraints are guaranteed to be valid all the time. BLFs are used for designing controller in several forms, for instance, Brunovsky form [43], feedback output-constrained system [44], stateconstrained system [45]. Compared with the generalized BLF whose constraints tend to be more conservative, an improved approach named integral barrier Lyapunov function (iBLF) [46] for nonlinear system allows the mixture of the initial state constraints and the errors. In this paper, iBLF is also presented in the analysis of system stability.
Inspired by the aforementioned discussion, a control scheme for a multiple homogeneous manipulators system grasping an object and tracking a predefined trajectory without knowing the precise system model is designed in this paper. For clarity, the main contributions of this paper are highlighted as follows: 1) The admittance model is introduced to work as a trajectory generator, which also depicts the interaction between the object and manipulators. Once the desired trajectory of the carried object is given, the reference trajectory for each manipulator could be generated online. 2) A novel iBLF function is presented in this paper. The iBLF function is of great significance in dealing state constraints. 3) An NN control scheme is proposed by combing a switch function with an adaptive neural network. Such a scheme could achieve a global approximation for the uncertainties in the system dynamics. Then a global RBFNN based adaptive control approach is specially designed for each manipulator to track the corresponding reference trajectories. By learning from system states and applying updating law for the weight matrix, optimal estimated weight matrix of the neural networks is obtained.
The structure of this paper is organised as follows. Section II formulates the problem and introduces some preliminaries. In section III, the multiple manipulators model is presented and the trajectory generator is designed. Section IV proposes the control strategy and the proof of system stability. In Section V, the simulation studies are conducted to verify the proposed methodology. The conclusion is drawn in Section VI.

A. Problem Formulation
Consider a system with multiple homogeneous robot manipulators, the coordinated task is for m manipulators to grasp and move an object along a predefined trajectory, as shown in Fig.1. The forces between the manipulators and the object can be measured by the pressure sensors. By employing admittance model [25] which mapping from forces to trajectories, the reference trajectories are obtained and tracked. In order to successfully fulfil the task, the following control objectives need to be accomplished: (i) the reference trajectory can be generated by using admittance method.
(ii) the end-effectors of the manipulators could track the reference end-effectors trajectories and the actual trajectories should be bounded within a predefined range due to the practical limitation.
(iii) the contact force between the object and manipulators could be restricted within an allowable range, in case of damaging to the object.
Before discussing the control of multiple manipulators, we need to take the contacts on the rigid object into consideration. Therefore, the following assumption is considered as a premise of the proposed method.
Assumption 1: The holding of the rigid object by multiple manipulators is form-closured, which means the contacts between manipulators and object prevent all kinds of motions of the rigid object.
Assumption 2: The end-effector is rigidly attached to the object, i.e., there is no relative movement between the endeffector and the object. In this paper, a class of rigid robot manipulators with ndegree of freedom (DOF) has been considered. The dynamics of the ith manipulator [47] is represented as follows: where q i ,q i ,q i ∈ R Ni are the joint position, joint velocity and joint acceleration, respectively, M e,i (q i ) ∈ R Ni×Ni denotes the symmetric positive definite inertia matrix, C e,i (q i ,q i ) ∈ R Ni×Ni is the Coriolis-centrifugal torque matrix, G e,i (q i ) ∈ R Ni denotes the gravity torque vector; τ e,i ∈ R Ni stands for the control input torque vector, F i is the force vector exerted on the end-effector of the ith manipulator, N i (i = 1, · · · , m) stands for the number of the ith manipulator's DOF, and m is the number of manipulators. Considering the position of the end-effector of the ith manipulator x i ∈ R ni in the Cartesian space, the kinematics of the ith manipulator is x i = φ i (q i ). Differentiate the kinematics with regards to time yieldṡ where J e,i (q i ) ∈ R ni×Ni is the Jacobian matrix. The dynamics of ith manipulator in Cartesian space can be represented as below The rigid robot manipulators depicted in (1)have the following properties [47] : Property 2: The symmetric and positive definite inertia matrix M i (q i ) is uniformly bounded, there exists lower limit constant m i > 0 and upper limit constant m i > 0, and M i (q i ) satisfies the following inequality k ci ||q i ||, and ||G i (q i )|| k gi , respectively, where k ci and k gi are positive constants.

C. RBFNN Constructure
According to the Weierstrass high order Approximation Theorem [48], given sufficient basis nodes, every continuous function F (Z) : Ω z → R over a compact set Ω Z ⊂ R Nin can be approximated as closely as desired by utilizing a basis set {s(z)} [49]. In this paper, Gaussian radial basis function is applied to approximate nonlinear function. The RBFNN structure is represented as follows: .., s Ns (Z)] T is the regressor vector with Gaussian radial basis function s i (·), ε(Z) ∈ Ω z is the approximation error and |ε(Z)| < ε * with constant ε * > 0 for all Z ∈ Ω Z , N in and N o are control input and output dimension, respectively, N s is the number of neural nodes.
It is noted that the ideal weight matrix W * is usually unknown in (7). In practice, the estimated weightŴ , which can be trained by a weight updating law, is often used to replace W * to approximate an nonlinear function, thus RBFNN in (7) can be represented as:

III. MODELLING PROCEDURE A. Dynamics of Multiple Manipulators
The dynamics of the object is: where F o ∈ R n is the resulting force exerted on the object, F d stands for the envirmonment force vector exerted on the object. From Assumption 2, let J i (x o ) ∈ R ni×n be the Jacobian matrix relating the the position of the end-effector x i ∈ R ni with the mass centre of the object x o ∈ R n , as shown in Fig.1, (3), the dynamics of m manipulators can be written in a compact form we can rewrite the dynamics of the ith manipulator in the following form: where By multiplying left side of the equation (13) by J T (x o ), integrating equation (10) and (11), we have In this paper, a damping-stiffness system model ( Fig.2) is considered: where C e and K e are damping and stiffness parameters of the interaction system predefined by the experimenter, and F e is the impedance force exerted on the end-effector. In general, a target admittance model in the Cartesian Space is depicted as below: where x ri and x di are reference trajectory and desired trajectory of the ith manipulator, respectively, as shown in Fig.2(a).
Applying Laplace transformation on (16), we can obtain: and we can further obtain that Therefore, once the impedance force F ei and the desired trajectory x di are obtained, the reference trajectory x ri could easily produced by the trajectory generator, the process can be seen in Fig.2(b). Remark 1: As shown in Fig.2(b), the inputs of trajectory generator are force F ei (S) and desired trajectory X di (S), the output is reference trajectory X ri (S). Therefore, the generator can be regarded as a filter, which means the reference trajectory is obtained by filtering the impedance force and desired trajectory.
According to the task that multi-manipulators system need to perform, the desired trajectory of the object x d in the Cartesian Space can be determined. Then the desired trajectory for each manipulator x di could be specified by the experimenter in advance, which locates on the object and closer to the centre than the actual contact point between manipulator and the object does. It can be seen clearly in Fig.2(a), as the reference trajectory x ri and the desired trajectory x di is different, the force between the manipulator and the object F ei = 0. Once the force obtained from pressure sensors mounted on the manipulators, the reference trajectory x ri can be derived by using admittance model. The relationship of desire object trajectory x d , the desired manipulator trajectory x di and the reference manipulator trajectory x ri can be seen in Fig.1.
Denoting e rdi = x ri − x di , we havė According to the knowledge of linear system, we can deduce that where x 0i is the initial state of the ith robot manipulator. Then we could get the reference trajectory: Assumption 3: The parameters of the ith admittance model is known, that means parameter matrix A i and B i also are known.
Assumption 4: The external forces exerted on the ith manipulator F ei can be measured by the pressure sensor mounted on the manipulator.
Remark 2: According to the assumption 3 and 4, considering the predetermined initial state x 0i and desired trajectory x di , we can conclude that the reference trajectory x ri could be generated by employing (21).

IV. CONTROL DESIGN
The objective of the NN control is to track the reference trajectory generated from section III-B. The framework of the multi-manipulator controller is shown in Fig.3.
Define the state variables x 1 = x o , x 2 =ẋ o , thus the kinetic equation in (14) can be rewritten as feedback form: A. Controller Design using iBLF The error variables are defined as follows where error vector z 1 = [z 1i , · · · , z 1n ] T , z 2 = [z 2i , · · · , z 2n ] T , α represents virtual control aiming to let the tracking error z 1 converge to a small neighbourhood of zero. Instead of using logarithmic Barrier Lyapunov Functions (BLF), a modified BLF method is employed in this paper, i.e. integral BLF, which allows the states constraints mixed with errors. For system (22), consider the iBLF candidate where k ci is positive constant satisfying |x i | < k ci . Then V is positive definite and continuously differentiable in the set |α i | < k ci for i = 1, · · · , n. We denote the set X := {x ∈ R n : |x i | < k ci , i = 1, · · · , n} ⊂ R n . Remark 3: The iBLF candidate V satisfies where |x di | A i < k ci , which is helpful for the proof of global stability. Additionally, a more useful conclusion can be given The proof of (26) can be seen in APPENDIX A. We use backstepping method to design the controller. It consists of several steps.

B. Global NN Control
In this paper, the switching functions Q(Z) ∈ R n×n are defined as follows: where with positive constants 0 < d 1,k < d 2,k and positive constant ω ≥ 1.
Noted that there are uncertainties in G(x 1 ), C(x 1 ,ẋ 1 ), M(x 1 ), and F d therefore F * o cannot be obtained in a real system. RBFNN are used to approximate the uncertainties in terms G(x 1 ), C(x 1 ,ẋ 1 ), M(x 1 ), F d . Define where Assumption 5: The functions f i (Z), i = 1, · · · , n are unknown and bounded, and there exits known nonnegative smooth functions f U i (Z) satisfying is the approximation of f i (Z) given bŷ whereŴ i is the actual weight vector; S i (·) is the basis function introduced by (8),Ŵ * i is optimal weight vector, Then, the proposed adaptive neural network controller is designed as where K r = diag(k r1 , · · · , k rn ) > 0, Φ a and Φ b are designed as with being a positive parameter, F U (Z) = diag[f U 1 (Z), · · · , f U n (Z)], Γ i (·) = tanh(·). To improve the performance of the control system, the updating law is designed aṡ where Θ i is positive symmetric matrix; σ i is positive constant.
where i is approximating error satisfying i ¯ i ,¯ i is positive constant.
Remark 4: In the adaptive NN controller proposed in (42), the term K 2 z 2 providing the error feedback, NN approximation function Φ a i , robust item Φ b i and switching function Q i (Z) work together to ensure the global tracking performance. The scheme is shown in Fig.4. Ω 0 is the neural networks approximation within the admissible region. It is noted that the scale of smooth function m(·) is |m(·)| = 1 in the compact set Ω 1 and |m(·)| = 0 outside the set Ω 2 . Therefore, in the compact set Ω 1 (|x 1 | d 1 ), term Φ a i works to approximatê f i (Z). Outside the set Ω 2 (|x 1 | d 2 ), robust term Φ b i works to pull the state x 1 back to Ω 2 . When i will pull the state back to the compact set Ω 1 .

C. Stability Analysis
Choose a positive iBLF as Considering control input (42) and the derivative of (46), we haveV Let us define i as E i , i = 1, · · · , n for the interval t ∈ (0, +∞), where (·) is ith element of a vector. Therefore, we obtain E = [E 1 , · · · , E n ] T .
Substituting the weight updating laws into (47), we havė (49) and the gain matrix K r are designed to satisfy |E i | k ri , i = 1, · · · , n, we have z T 2 (K r sgn(z 2 ) − E) 0, thereforė where Therefore,V 3 is a negative definite function, and as time increases, V 3 will decay into a region close to zero.
Theorem 1: For the multiple manipulators system dynamics (14), under the NN controller (42) with weight updating law (44) and admittance trajectory generator (21), the system (14) would obtain satisfying control performance through environment-robot interaction force is imposed on (14). The errors will converge into a small neighbourhood near zero by designing appropriate controller parameters and updating laws. The system states still remain in the predefined region and the tracking error z 1 would converge into the compact Ω z1 := z 1 ∈ R n | |z 1i | √ 2B , i = 1, · · · , n. The tracking error z 2 converge into the compact set Ω z2 := z 2 ∈ R n | |z 2 |

V. SIMULATION STUDIES
In this section, to illustrate the efficacy of the proposed adaptive control method, two sets of simulation studies are conducted and the results are compared, and the model of planar 2-link manipulators are employed in this paper as shown in Fig.5. The manipulators with force sensor on the endeffector cooperate to move an object along a predefined desired trajectory.

A. Robot Manipulator Model
Consider two homogeneous manipulators with 2 revolute joints [50]. These manipulators with force sensor mounted on the end-effectors share same parameters as listed in TABLE I, where m i , l i and I i , i = 1, 2 are mass, length and inertia of the ith link, respectively.  The desired trajectory of the mass centre of the object x do is which is a circle with the centre at [0, 0] T and radius being 0.1m and there is no rotation during the translation of the object. x d1 and x d2 are the x-coordinate and y-coordinate of the object, respectively, θ is the orientation of the object. The parameters of dynamics for the object are given as Denote q 1 , q 2 as the first and the second joint angles of the first manipulator, respectively. The position of the endeffector of the first manipulator in the Cartesian space are x 1 and x 2 , respectively. Correspondingly, q 3 , q 4 , x 3 and x 4 are joint angles and position of the end-effector for the second manipulator. Therefore we have the following kinematic relationship , 0] T are the position of the bases of two manipulators. Therefore, the Jacobian matrix between the joint space of the manipulator and the corresponding Cartesian space is J e,1 (q) = −l 1 sinq 1 − l 2 sin(q 1 + q 2 ) −l 2 sin(q 1 + q 2 ) l 1 cosq 1 + l 2 cos(q 1 + q 2 ) l 2 cos(q 1 + q 2 ) .
The Jacobian matrix between the end-effectors and the object J o is The dynamic parameters of the ith(i = 1, 2) manipulator are where l c1 and l c2 are the centre of the first and the second link of the manipulator, respectively. Here q 11 , q 12 , q 21 and q 22 stand for q 1 , q 2 , q 3 and q 4 , respectively.  Fig. 6 -Fig. 10. Fig. 6 and 7 give the tracking performance of the PD controller and we can find that the tracking errors change with respect to time with a relative high vibration, as shown in Fig. 8 and 9. In Fig.  10, we can also find the control force is not stable and the amplitude is large.
2) Adaptive NN control: For the neural network, the number of NN nodes are chosen as 2 12 , the centre of NN nodes are evenly distributed in [−1, 1], and the variance is 50. The NN weights are all initialized as zeros, the gain of NN adaptive law are chosen as Θ 1 = 120, Θ 2 = 150, Θ 3 = 120, and σ 1 = σ 2 = σ 3 = 0.001. The designed parameters of the controller are selected as K 1 = diag (15,20,15), K 2 = diag (11,15,10) and    The simulation results are shown in Fig. 11 -Fig. 16. Fig.  11 and 12 show the tracking trajectories of the NN control. In the proposed iBLF function based control method, positive constant k ci stands for the output contraints. In our simulation, we set k c1 = k c2 = k c3 = 0.35. As shown in Fig. 12, the tracking trajectories remain in the predefined bound (green lines) during the whole simulation process, which validate the effectiveness of output constraints functionality of the proposed control scheme. From Fig. 13 and 14, the tracking errors converge to a small neighbourhood of zero quickly within 2 seconds, that is due to the learning process of the neural networks. In Fig.15, we can find that the control torque is stable at a relative low level when the NN weight matrix reach the optimal value. As depicted in Fig. 16, it takes less than 2 seconds for the estimated NN weight matrices to rapidly converge to stable values. In general, the NN control can make the controller more stable and efficient compared with PD control. These results verify the efficacy of the proposed method for the position tracking of the object carried by multiple robot manipulators.

VI. CONCLUSION
This paper presents an adaptive control methodology based on admittance model for multiple manipulators transporting an rigid object cooperatively along a predefined desired trajectory. The admittance model is utilized to generate reference trajectory online for each manipulator according to the desired path of the rigid object. A type of iBLF is introduced to tackle the constraints due to the physical and environmental limits. As for the uncertainties of the manipulator dynamics, adaptive neural networks are employed to approximate these uncertain terms. The switching function provides the required the global stability of the closed-loop. The simulation results carried on planar 2-link manipulators validate the efficacy of the presented method. Although the performance of the adaptive NN control is better than that of the conventional PD control when the weight vector W is stable, it is noted that it takes a longer time for NN control to catch up with the desired signal. To further validate the proposed method, conducting some experiment is one of our future work.

APPENDIX
A. Proof of (26) in Remark 3 Step 1 (the left part of the inequality) Denote g(z i ) = . It is noted that in the compact set X we have k 2 ci − x 2 i > 0. For case z i > 0, we have ∂g ∂zi < 0; For case z i < 0, we have ∂g ∂zi > 0. Since g(z i ) = 0 at z i = 0, therefore we can draw the conclusion that g i > 0 in the set X . That means Step 2 (the right part of the inequality) Define p i (σ, α i ) = σk 2 ci k 2 ci −(σ+αi) 2 , we have which is positive in the set |σ + α i | < k ci . Since p i (0, α i−1 ) = 0 for all |α i | < k ci and p i (σ, α i−1 ) is monotonically increasing with σ in the set |σ + α i | < k ci , we further have zi 0 p i (σ, α i−1 )dσ z i p i (σ, α i−1 ) for |σ + α i | < k ci , which leads to the right part of (26) after substituting for p i .