Decentralized control for harmonic drive–based modular and reconfigurable robots with uncertain environment contact

Abstract

In this article, a decentralized control strategy is presented for harmonic drive–based modular and reconfigurable robots with uncertain environment contact. Unlike conventional methods that rely on robot–environment contact model or force/torque sensing, this article addresses the problem of controlling modular and reconfigurable robots in contact with uncertain environment using only encoder data of each joint module. By employing a control-oriented harmonic drive model, the dynamic model of modular and reconfigurable robot is formulated as a synthesis of interconnected subsystems, in which the interconnected joint couplings are with small magnitudes. Based on the integral sliding mode control technique and the adaptive super-twisting algorithm, the decentralized controller is designed to compensate model uncertainty in which the up-bound is unknown. The stability of the modular and reconfigurable robot system is proved using Lyapunov theory. Finally, simulations are conducted for 2-degree-of-freedom modular and reconfigurable robots with different configurations under the situations of dynamic contact and collision to investigate the advantage of the proposed approach.

Introduction

Modular and reconfigurable robots (MRRs) comprise the robot modules, which contain power supplies, processing systems, actuators, and sensors. These modules are assembled to desirable configurations with standard interfaces to satisfy the needs of various tasks with a complex working environment. Therefore, the MRRs are often utilized in unknown and dangerous environments, for example, space exploration, disaster assistance, and high-/low-temperature operations. Besides, in the face of a complex and uncertain environment, MRRs need suitable control systems to provide stability and accuracy after reconfigurations.

Besides the property of reconfigurability, an ideal MRR system should be designed to possess the ability of physically interacting with the external environment, for example, the dynamic contact and collision with time-varying environment constraint, and ensuring the accuracy of control performance in trajectory tracking. Some insightful investigations are presented for addressing the problems of robust stabilization and uncertainty compensation of mechatronic systems based on a novel robust integral of the sign of the error (RISE) feedback method^1,2 and extended state observer (ESO).^3,4 Besides, to tackle the concerns of raising interaction stability, numerous researches have focused on developing robot–environment contact models^5,6 which are implemented for investigating the cases of position control^7,8 and force control.^9,10 However, the limitation of the existing contact model–based control methods is that it relies on the exact knowledge of the robot and environment dynamics. Indeed, MRRs are always designed for satisfying various task requirements under uncertain environments; this means the dynamic model of MRRs will be changed under different configurations, and the contact model is also difficult to be obtained. To avoid using the enclosed forms of dynamic models and contact models, and seeking out high-performance control for MRRs, joint torque feedback techniques are attached with attentions of both robotic researchers and industrial manufacturers. Several successful applications of direct joint torque sensing techniques are published in the literatures.^11–13 However, using joint torque sensor, measurements are known directly with many drawbacks and may jeopardize the simplicity, reliability, and mechanical ruggedness of the modular robots. For one thing, the strain gauges, which are the most important components of the torque sensors, are known with intensive sensitivity to peripheral temperature variations; for another thing, the sensor data are contaminated with too much noise and unstable states, which are caused by the narrow bandwidth of the torque information. Therefore, an excellent control system for an MRR with uncertain environment contact should be able to eliminate the needs of contact models and force/torque sensing, and it relies only on encoder data of each joint module.

A significant property of MRRs is that the robot modules can be appended, removed, and replaced optionally without the requirement of readjusting the gains and parameters of the controllers in other modules, so that the MRRs can be reassembled to satisfy different requirements of various tasks. Decentralized control is a promising strategy in particular for implementing the control system of MRRs. The decentralized control design always arises from many complex conditions where there exist physical restrictions on information interchange among the subsystems for which there is an insufficient capacity of possessing a centralized controller. Furthermore, a decentralized control architecture based on local modular information not only provides the flexibility of structure, which is applicative for all the robot configurations, but also decrease the computation consumptions that are critical for the real-time operation of the MRR systems. A number of researches have been made toward the decentralized control method for MRRs. Based on generic locomotion algorithm, a decentralized control method is presented for a lattice-based reconfigurable robot. ¹⁴ A joint torque sensing–based decentralized control method is investigated to address the problems of fault detection and tolerant control of MRRs. ¹⁵ A decentralized reconfiguration strategy is proposed for MRRs assembling furniture-like structures from Roombots metamodules. ¹⁶ A decentralized robust control algorithm is presented for MRRs with harmonic drive transmission based on backstepping technique and Lyapunov theory. ¹⁷ In our previous researches, a stable decentralized control method is proposed for MRRs by employing adaptive fuzzy algorithm and sliding mode control technique for satisfying the conception of modular design; ¹⁸ building on these results, we also proposed a reinforcement learning-based decentralized optimal control method ¹⁹ and a decentralized integral nested sliding mode control scheme for MRRs under the situation of time-varying constraint. ²⁰ However, the main disadvantage of the above-mentioned methods is that the interconnected joint couplings, which include the Coriolis, centripetal torques, gravitational torques, and so on regarding all joint modules, are simply treated as external disturbances or model uncertainties with known up-bounds. Unfortunately, the magnitudes of the interconnected couplings between the joints in these strategies are very large, so that the decentralized compensating controllers, which are designed to handle the uncertainties, should be with higher control gains, and this may cause controller chattering and further spoiling the performance of the controllers. Besides, to the best of the authors’ knowledge, there are few types of research that focus on investigating the decentralized control strategy of the MRRs under uncertain environment contact. Therefore, it is meaningful to devise an appropriate decentralized control method for MRRs, particularly under the situations of uncertain environment contact, and the controllers should be with the capability of effectively compensating the interconnected joint couplings and reducing the chattering effect of the controllers.

In this article, we present a decentralized control strategy for harmonic drive–based MRRs with uncertain environment contact. First, by employing a control-oriented harmonic drive model that considers the compliance of both the wave generator and flexspline, the dynamic model of the MRR is formulated as a synthesis of interconnected subsystems, in which the magnitudes of the interconnected couplings between the joints are decreased significantly. Second, based on the integral sliding mode control (ISMC) technique^21,22 and the adaptive super-twisting algorithm (ASTA),^23,24 a decentralized controller is proposed which utilizes only encoder data of each isolate joint module to compensate model uncertainty and to reduce the chattering effect of the controller, as well as the stability of the closed-loop MRR system is proved using Lyapunov theory. Finally, simulations are conducted for 2-degree-of-freedom (DOF) MRRs with different configurations under the situations of dynamic contact and collision to investigate the advantage of the proposed approach.

The main contributions of this article can be summarized as follows:

• Unlike conventional methods that rely on robot–environment contact model or force/torque sensing, this article focuses on addressing the trajectory tracking control problems of MRRs with uncertain environment contact that uses only encoder data on the motor side and link side of each isolate joint module, and the joint modules can be used to construct various configurations without the requirement of readjusting the control gains and parameters.

• In this article, we consider the MRRs are with model uncertainties, which contain friction modeling error and uncompensated interconnected joint couplings, and then model uncertainties are compensated under ASTA-based ISMC, which can reduce the chattering effect of the controller and adjust the control gains and parameters adaptively while the up-bounds of the uncertainties are unknown.

The remainder of this article is organized as follows: section “Model and problem formulation” outlines the dynamic model of harmonic drive–based MRR. The proposed decentralized control strategy is presented in section “Decentralized ISMC design based on ASTA,” and the simulations are illustrated in section “Simulations.” Concluding remarks are given in section “Conclusion.”

Model and problem formulation

In this section, we formulated the dynamic model of harmonic drive–based MRRs with uncertain environment contact. Two parts are included in this section; in the first part, a control-oriented harmonic drive model is constructed which considers both wave generator compliance and flexspline compliance, and in the second part, the dynamic model of MRR with uncertain environment contact is formulated based on the harmonic drive model.

Harmonic drive model

Considering an MRR comprising n modules, and each module is integrated with a harmonic drive–based rotary joint. By referring to the literatures,^25,26 which are performed for analytically describing the structural property of harmonic drive, one concludes the following ideal kinematic relationship

$${\theta }_{\mathrm{wi}}={\gamma }_i{\theta }_{\mathrm{fi}}-({\gamma }_i+1){\theta }_{\mathrm{ci}}$$

(1)

where ${\gamma }_i$ indicates the gear ratio; ${\theta }_{\mathrm{wi}}$ , ${\theta }_{\mathrm{ci}}$ , and ${\theta }_{\mathrm{fi}}$ represent the angle positions of the wave generator, circular spline, and flexspline, respectively; note that subscript i denotes the ith module. The static force balance between the components is illustrated as follows

$$\begin{array}{c}\hfill {\tau }_{\mathrm{wi}}=-\frac{1}{{\gamma }_i+1}{\tau }_{\mathrm{ci}}=\frac{1}{{\gamma }_i}{\tau }_{\mathrm{fi}}\hfill \end{array}$$

(2)

where ${\tau }_{\mathrm{wi}}$ , ${\tau }_{\mathrm{ci}}$ , and ${\tau }_{\mathrm{fi}}$ refer to the torque at the wave generator, circular spline, and flexspline, respectively. Besides, since the circular spline is fixed leaving the wave generator and flexspline as the input and output, respectively, so that one can conclude the relation ${\theta }_{\mathrm{ci}}={\omega }_{\mathrm{ci}}=0$ , and the reaction torque ${\tau }_{\mathrm{ci}}$ is not of concern.

However, the experimental measurements of the kinematic relationship between the wave generator and the flexspline provided in the literatures^27,28 clearly show that the link-side position measurement is not linearly related to the motor-side position measurement. The causation of this nonlinearity is due to the compliance in the harmonic drive components, and the kinematic error may attribute to gear meshing and machining errors. Given this ideal kinematic relationship which describes the motion and force constraints present in harmonic drives, the remaining effects can be incorporated by introducing compliance and kinematic error properties.

By combining the link-side and motor-side position measurements, the torsional compliance of the flexspline and wave generator is described as follows

$$\Delta {\theta }_{\mathrm{fi}}={\theta }_{\mathrm{fOi}}-{\theta }_{\mathrm{fIi}},\Delta {\theta }_{\mathrm{wi}}={\theta }_{\mathrm{wOi}}-{\theta }_{\mathrm{wIi}}$$

(3)

where ${\theta }_{\mathrm{fOi}}$ and ${\theta }_{\mathrm{fIi}}$ indicate the angle positions of the flexspline at the load side and the gear side in the ith joint, respectively; ${\theta }_{\mathrm{wOi}}$ and ${\theta }_{\mathrm{wIi}}$ represent the angle positions of the wave generator outside part and the center part in the ith joint, respectively. Only the flexible output position (joint’s angle) ${\theta }_{\mathrm{fOi}}$ and the wave generator input position (motor’s angle) ${\theta }_{\mathrm{wIi}}$ can be measured using the link-side and motor-side encoders. These obtainable positions are used to compute the torsional angle of the harmonic drive in the ith joint utilizing the following relation

$$\Delta {\theta }_i={\theta }_{\mathrm{fOi}}-\frac{{\theta }_{\mathrm{wIi}}}{{\gamma }_i}$$

(4)

By adding and subtracting the terms ${\theta }_{\mathrm{wOi}}/{\gamma }_i$ and ${\theta }_{\mathrm{fIi}}$ into equation (5), we can get that

$$\begin{array}{c}\multicolumn{1}{c}{\Delta {\theta }_i={\theta }_{\mathrm{fOi}}-{\theta }_{\mathrm{fIi}}+(\frac{{\theta }_{\mathrm{wOi}}}{{\gamma }_i}-\frac{{\theta }_{\mathrm{wIi}}}{{\gamma }_i})+({\theta }_{\mathrm{fIi}}-\frac{{\theta }_{\mathrm{wOi}}}{{\gamma }_i})}\\ \multicolumn{1}{c}{=\Delta {\theta }_{\mathrm{fi}}+\frac{\Delta {\theta }_{\mathrm{wi}}}{{\gamma }_i}+\Delta {\theta }_{\mathrm{erri}}}\end{array}$$

(5)

where the kinematic error of the harmonic drive torsional angle $\Delta {\theta }_{\mathrm{erri}}$ is defined as

$$\Delta {\theta }_{\mathrm{erri}}={\theta }_{\mathrm{fIi}}-\frac{{\theta }_{\mathrm{wOi}}}{{\gamma }_i}$$

(6)

To compensate the kinematic error $\Delta {\theta }_{\mathrm{erri}}$ , a kinematic error model of the harmonic drive is developed as follows

$$\{\begin{array}{c}\hfill \Delta {\theta }_{\mathrm{erri}}={\tilde{\theta }}_{\mathrm{wi}}+{\tilde{\theta }}_{\mathrm{fi}}\hfill \\ \multicolumn{1}{c}{{\tilde{\theta }}_{\mathrm{wi}}=a_{w0}+a_{w1}\cos ({\omega }_w{\theta }_{\mathrm{wIi}})+b_{w1}\sin ({\omega }_w{\theta }_{\mathrm{wIi}})}\\ \multicolumn{1}{c}{+a_{w2}\cos (2{\omega }_w{\theta }_{\mathrm{wIi}})+b_{w2}\sin (2{\omega }_w{\theta }_{\mathrm{wIi}})}\\ \multicolumn{1}{c}{{\tilde{\theta }}_{\mathrm{fi}}=a_{f0}+a_{f1}\cos ({\omega }_f{\theta }_{\mathrm{fOi}})+b_{f1}\sin ({\omega }_f{\theta }_{\mathrm{fOi}})}\end{array}$$

(7)

where ${\tilde{\theta }}_{\mathrm{fi}}$ and ${\tilde{\theta }}_{\mathrm{wi}}$ are the kinematic errors of the flexspline and wave generator, respectively. The determined parameters of the kinematic error model in equation (7) are obtained utilizing the experimental data in the literature ²⁹ and are formulated in Table 1.

Table 1.

Parameters of the dynamic model.

Name	Value	Name	Value	Name	Value
$I_{\mathrm{mi}}$	120	$a_0$	0	$b_{w2}$	90.94e−5
$c_w$	83.5	$a_{w1}$	15.94e−5	$b_{f1}$	−43.09e−4
$c_f$	8.9e−2	$a_{w2}$	32.88e−5	${\omega }_w$	17.34e−3
$k_{w0}$	1.33	$a_{f1}$	−14.76e−5	${\omega }_f$	17.31e−3
$k_{f0}$	8.3e+3	$b_{w1}$	−18.91e−4	${\gamma }_i$	100

According to the typical stiffness property of a harmonic drive, which is proposed by Zhang et al., ³⁰ it is obvious that the local elastic coefficient increases as the flexspline torque increases. Therefore, one can define the local elastic coefficient of the flexspline $k_{\mathrm{fL}}$ as

$$k_{\mathrm{fL}}=\frac{d{\tau }_{\mathrm{fi}}}{d\Delta {\theta }_{\mathrm{fi}}}$$

(8)

Considering the symmetric performance of harmonic drive stiffness, the local elastic coefficient $k_{\mathrm{fL}}$ is approximated by

$$k_{\mathrm{fL}}=k_{f0}(1+{(c_f{\tau }_{\mathrm{fi}})}^2)$$

(9)

where $c_f$ and $k_{f0}$ are the constants to be determined. If $k_{f0}\ne 0$ , one can calculate the flexspline torsion $\Delta {\theta }_{\mathrm{fi}}$ as

$$\Delta {\theta }_{\mathrm{fi}}=\underset{0}{\overset{{\tau }_{\mathrm{fi}}}{\int }}\frac{d{\tau }_{\mathrm{fi}}}{k_{\mathrm{fL}}}=\frac{\arctan (c_f{\tau }_{\mathrm{fi}})}{c_f{k}_{f0}}$$

(10)

Besides, the deformation range of harmonic drive drops to zero sharply at rated torque, which implies that the wave generator stiffness increases sharply. For the purpose of replicating the hysteresis shape of this stiffness property, the local elastic coefficient of wave generator is modeled as

$$k_{\mathrm{wL}}=k_{w0}{e}^{c_w\left|{\tau }_{\mathrm{wi}}\right|}$$

(11)

where $c_w$ and $k_{w0}$ are the constants to be determined. If $k_{w0}\ne 0$ , one can calculate the torsional angle of wave generator using the following relation

$$\Delta {\theta }_{\mathrm{wi}}=\underset{0}{\overset{{\tau }_{\mathrm{wi}}}{\int }}\frac{d{\tau }_{\mathrm{wi}}}{k_{\mathrm{wL}}}=\frac{\mathrm{sgn}({\tau }_{\mathrm{wi}})}{c_w{k}_{w0}}(1-e^{-c_w\left|{\tau }_{\mathrm{wi}}\right|})$$

(12)

where $\mathrm{sgn}(·)$ is a sign function, defined as

$$\mathrm{sgn}({\tau }_{\mathrm{wi}})=\{\begin{array}{c}\hfill 1\mathrm{for}{\tau }_{\mathrm{wi}}>0\hfill \\ \multicolumn{1}{c}{0\mathrm{for}{\tau }_{\mathrm{wi}}=0}\\ \multicolumn{1}{c}{-1\mathrm{for}{\tau }_{\mathrm{wi}}<0}\end{array}$$

(13)

Substituting equations (10) and (12) into equation (5), define the total torsional angle of a harmonic drive as follows

$$\begin{array}{c}\multicolumn{1}{c}{\Delta {\theta }_i=\Delta {\theta }_{\mathrm{fi}}+\frac{\Delta {\theta }_{\mathrm{wi}}}{{\gamma }_i}+\Delta {\theta }_{\mathrm{erri}}}\\ \multicolumn{1}{c}{=\frac{\arctan (c_f{\tau }_{\mathrm{fi}})}{c_f{k}_{f0}}+\frac{\mathrm{sgn}({\tau }_{\mathrm{wi}})}{{\gamma }_i{c}_w{k}_{w0}}(1-e^{-c_w\left|{\tau }_{\mathrm{wi}}\right|})+\Delta {\theta }_{\mathrm{erri}}}\end{array}$$

(14)

Deforming the equation above, one obtains that

$${\tau }_{\mathrm{fi}}=\frac{1}{c_f}\tan \left(c_f{k}_{f0}\left(\begin{array}{c}\hfill \Delta {\theta }_i-\Delta {\theta }_{\mathrm{erri}}\hfill \\ \multicolumn{1}{c}{-\frac{\mathrm{sgn}({\tau }_{\mathrm{wi}})(1-e^{-c_w\left|{\tau }_{\mathrm{wi}}\right|})}{{\gamma }_i{c}_w{k}_{w0}}}\end{array}\right)\right)$$

(15)

In equation (15), ${\tau }_{\mathrm{fi}}$ denotes the flexspline torque of harmonic drive including the environment constrained torque and the load torque; ${\tau }_{\mathrm{wi}}$ denotes the wave generator torque which is obtained using the motor torque command.

Dynamic model of harmonic drive–based MRR

According to the torque sensing–based robot manipulator dynamic modeling method proposed by Imura et al., ³¹ and the encoder data–based flexspline torque estimation method given in equation (15), in this part, the dynamic model of harmonic drive–based MRR system is formulated as a synthesis of interconnected subsystems, in which the ith subsystem dynamic model is formulated as follows

$$\begin{array}{c}\hfill I_{\mathrm{mi}}{\gamma }_i{\stackrel{··}{\theta }}_i+f_i({\theta }_i,{\stackrel{·}{\theta }}_i)+I_{\mathrm{mi}}\sum _{j=1}^{i-1}{z}_{\mathrm{mi}}^T{z}_{\mathrm{lj}}{\stackrel{··}{\theta }}_j\hfill \\ \multicolumn{1}{c}{+I_{\mathrm{mi}}\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{z}_{\mathrm{mi}}^T(z_{\mathrm{lk}}\times z_{\mathrm{lj}}){\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j+\frac{{\tau }_{\mathrm{fi}}}{{\gamma }_i}={\tau }_i}\end{array}$$

(16)

where $I_{\mathrm{mi}}$ is the moment of inertia of the rotor about the axis of rotation, ${\tau }_i$ indicates the motor output torque, $I_{\mathrm{mi}}\sum _{j=1}^{i-1}{z}_{\mathrm{mi}}^T{z}_{\mathrm{lj}}{\stackrel{··}{\theta }}_j+I_{\mathrm{mi}}\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{z}_{\mathrm{mi}}^T(z_{\mathrm{lk}}\times z_{\mathrm{lj}}){\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j$ is the interconnected joint couplings, $f_i({\theta }_i,{\stackrel{·}{\theta }}_i)$ denotes the friction torque; $z_{\mathrm{mi}}$ and $z_{\mathrm{li}}$ denote the unity vectors along the rotation axis of the ith rotor and joint, respectively; and ${\tau }_{\mathrm{fi}}/{\gamma }_i$ reflects the dynamic information of the external environment contact torque and the load torque. Then, according to the dynamic model analysis in our previous research, ³² model uncertainties of the subsystem dynamic model equation (16) may exist in terms of $f_i({\theta }_i,{\stackrel{·}{\theta }}_i)$ , $I_{\mathrm{mi}}\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{z}_{\mathrm{mi}}^T(z_{\mathrm{lk}}\times z_{\mathrm{lj}}){\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j$ , and $I_{\mathrm{mi}}\sum _{j=1}^{i-1}{z}_{\mathrm{mi}}^T{z}_{\mathrm{lj}}{\stackrel{··}{\theta }}_j$ .

For the friction torque term $f_i({\theta }_i,{\stackrel{·}{\theta }}_i)$ , by referring to the investigation, ³³ one can define it as a function with regard to the position and velocity of the ith joint, which is given as follows

$$f_i({\theta }_i,{\stackrel{·}{\theta }}_i)=b_{\mathrm{fi}}{\stackrel{·}{\theta }}_i+(f_{\mathrm{ci}}+f_{\mathrm{si}}{e}^{(-f_{\tau i}{\stackrel{·}{\theta }}_i^2)})\mathrm{sgn}({\stackrel{·}{\theta }}_i)+f_{\mathrm{pi}}({\theta }_i,{\stackrel{·}{\theta }}_i)$$

(17)

where $b_{\mathrm{fi}}$ is the viscous friction coefficient, $f_{\mathrm{ci}}$ is the relevant parameter of Coulomb friction, $f_{\mathrm{si}}$ represents the relevant parameter of static friction, $f_{\tau i}$ is a position parameter for the Stribeck effect, and $f_{\mathrm{pi}}({\theta }_i,{\stackrel{·}{\theta }}_i)$ represents the position dependence friction and other friction modeling errors. According to the linearization scheme proposed by Liu, ³⁴ assuming the nominal values of $f_{\mathrm{si}}$ and $f_{\tau i}$ are close to their actual values, we can linearize the Stribeck effect $f_{\mathrm{si}}{e}^{(-f_{\tau i}{\stackrel{·}{\theta }}_i^2)}$ at nominal parameter values ${\hat{f}}_{\mathrm{si}}$ and ${\hat{f}}_{\tau i}$ , and friction torque $f_i({\theta }_i,{\stackrel{·}{\theta }}_i)$ can be approximated by

$$\begin{array}{c}\multicolumn{1}{c}{f_i({\theta }_i,{\stackrel{·}{\theta }}_i)\approx {\hat{b}}_{\mathrm{fi}}{\stackrel{·}{\theta }}_i+({\hat{f}}_{\mathrm{ci}}+{\hat{f}}_{\mathrm{si}}{e}^{(-{\hat{f}}_{\tau i}{\stackrel{·}{\theta }}_i^2)})\mathrm{sgn}({\stackrel{·}{\theta }}_i)}\\ \multicolumn{1}{c}{+f_{\mathrm{pi}}({\theta }_i,{\stackrel{·}{\theta }}_i)+Y({\stackrel{·}{\theta }}_i){\tilde{F}}_i}\end{array}$$

(18)

where ${\hat{b}}_{\mathrm{fi}}$ , ${\hat{f}}_{\mathrm{ci}}$ , ${\hat{f}}_{\mathrm{si}}$ , and ${\hat{f}}_{\tau i}$ denote the estimated values of the friction parameters; ${\tilde{F}}_i$ , which indicates the parametric uncertainty of the friction, is defined as

$${\tilde{F}}_i=[\begin{array}{cccc}\hfill b_{\mathrm{fi}}-{\hat{b}}_{\mathrm{fi}}\hfill & \hfill f_{\mathrm{ci}}-{\hat{f}}_{\mathrm{ci}}\hfill & \hfill f_{\mathrm{si}}-{\hat{f}}_{\mathrm{si}}\hfill & \hfill f_{\tau i}-{\hat{f}}_{\tau i}\hfill \end{array}{]}^T$$

(19)

and $Y({\stackrel{·}{\theta }}_i)$ is defined as

$$Y({\stackrel{·}{\theta }}_i)=[{\stackrel{·}{\theta }}_i,\mathrm{sgn}({\stackrel{·}{\theta }}_i),e^{(-{\hat{f}}_{\tau i}{\stackrel{·}{\theta }}_i^2)}\mathrm{sgn}({\stackrel{·}{\theta }}_i),-{\hat{f}}_{\mathrm{si}}{\stackrel{·}{\theta }}_i^2{e}^{(-{\hat{f}}_{\tau i}{\stackrel{·}{\theta }}_i^2)}\mathrm{sgn}({\stackrel{·}{\theta }}_i)]$$

(20)

Besides, for terms $I_{\mathrm{mi}}\sum _{j=1}^{i-1}{z}_{\mathrm{mi}}^T{z}_{\mathrm{lj}}{\stackrel{··}{\theta }}_j$ and $I_{\mathrm{mi}}\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{z}_{\mathrm{mi}}^T(z_{\mathrm{lk}}\times z_{\mathrm{lj}}){\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j$ in equation (16), one can rewrite them as

$$\begin{array}{c}\multicolumn{1}{c}{I_{\mathrm{mi}}\sum _{j=1}^{i-1}{z}_{\mathrm{mi}}^T{z}_{\mathrm{lj}}{\stackrel{··}{\theta }}_j=I_{\mathrm{mi}}\sum _{j=1}^{i-1}{D}_j^i{\stackrel{··}{\theta }}_j}\\ \multicolumn{1}{c}{=\sum _{j=1}^{i-1}[\begin{array}{cc}\hfill I_{\mathrm{mi}}{\hat{D}}_j^i\hfill & \hfill I_{\mathrm{mi}}\hfill \end{array}][\begin{array}{cc}\hfill {\stackrel{··}{\theta }}_j\hfill & \hfill {\tilde{D}}_j^i{\stackrel{··}{\theta }}_j\hfill \end{array}{]}^T}\\ \multicolumn{1}{c}{=\sum _{j=1}^{i-1}{I}_j^i{U}_j^i}\end{array}$$

(21)

$$\begin{array}{c}\multicolumn{1}{c}{I_{\mathrm{mi}}\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{z}_{\mathrm{mi}}^T(z_{\mathrm{lk}}\times z_{\mathrm{lj}}){\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j=I_{\mathrm{mi}}\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{\Theta }_{\mathrm{kj}}^i{\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j}\\ \multicolumn{1}{c}{=\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}[\begin{array}{cc}\hfill I_{\mathrm{mi}}{\hat{\Theta }}_{\mathrm{kj}}^i\hfill & \hfill I_{\mathrm{mi}}\hfill \end{array}][\begin{array}{cc}\hfill {\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j\hfill & \hfill {\tilde{\Theta }}_{\mathrm{kj}}^i{\stackrel{·}{\theta }}_k{\stackrel{·}{\theta }}_j\hfill \end{array}{]}^T}\\ \multicolumn{1}{c}{=\sum _{j=2}^{i-1}\sum _{k=1}^{j-1}{J}_{\mathrm{kj}}^i{V}_{\mathrm{kj}}^i}\end{array}$$

(22)

In equations (21) and (22), one obtains the relationships of $D_j^i=z_{\mathrm{mi}}^T{z}_{\mathrm{lj}}$ , ${\Theta }_{\mathrm{kj}}^i=z_{\mathrm{mi}}^T(z_{\mathrm{lk}}\times z_{\mathrm{lj}})$ , ${\hat{D}}_j^i=D_j^i-{\tilde{D}}_j^i$ , and ${\hat{\Theta }}_{\mathrm{kj}}^i={\Theta }_{\mathrm{kj}}^i-{\tilde{\Theta }}_{\mathrm{kj}}^i$ , where ${\hat{D}}_j^i$ is the dot product of $z_{\mathrm{mi}}$ and $z_{\mathrm{lj}}$ and ${\tilde{D}}_j^i$ denotes the alignment error. Similarly, ${\hat{\Theta }}_{\mathrm{kj}}^i$ denotes the dot product of $z_{\mathrm{mi}}$ and $z_{\mathrm{lk}}\times z_{\mathrm{lj}}$ ; ${\tilde{\Theta }}_{\mathrm{kj}}^i$ is the alignment error.

Then, on the basis of the decomposition scheme of model uncertainty, which is represented by Liu and Goldenberg, ³⁵ the model uncertain terms of ${\tilde{F}}_i$ , $U_j^i$ , and $V_{\mathrm{kj}}^i$ in equations (19), (21), and (22) are decomposed as the summations of constant parts and variable parts, that is

$$\begin{array}{c}\hfill {\tilde{F}}_i={\tilde{F}}_{\mathrm{ic}}+{\tilde{F}}_{\mathrm{iv}},U_j^i=U_{\mathrm{jc}}^i+U_{\mathrm{jv}}^i,V_{\mathrm{kj}}^i=V_{\mathrm{kjc}}^i+V_{\mathrm{kjv}}^i\hfill \end{array}$$

(23)

where superscripts c and v denote the constant and variable parts, respectively. By analyzing the parameters and variables in the dynamic model equation (16) and the uncertain terms in equation (23), one obtains that model uncertainties ${\tilde{F}}_i$ , $U_j^i$ , and $V_{\mathrm{kj}}^i$ possess unknown up-bounds. Then, one can classify these model uncertainties into the following three categories.

1. Constant parametric uncertainties, which exist in terms of ${\tilde{F}}_{ic}$ , $U_{\mathrm{jc}}^i$ , and $V_{\mathrm{kjc}}^i$ , do not vary along with the reconfiguration of MRRs.

2. Variable parametric uncertainties, which exist in terms of $U_{\mathrm{jv}}^i$ and $V_{\mathrm{kjv}}^i$ , result in the reconfiguration of MRRs. Besides, ${\tilde{F}}_{\mathrm{iv}}$ may vary along with the MRR moving to different positions.

3. Nonparametric model uncertainty contains the nonparametric friction term $f_{\mathrm{pi}}({\theta }_i,{\stackrel{·}{\theta }}_i)$ . In the following contents, we use $f_{\mathrm{pi}}$ instead of $f_{\mathrm{pi}}({\theta }_i,{\stackrel{·}{\theta }}_i)$ to simplify the expression.

Substituting equations (21) and (22) into equation (16), one can rewrite the dynamic model of ith subsystem as follows

$${\stackrel{··}{\theta }}_i=-B_i\left(\frac{{\tau }_{fi}}{{\gamma }_i}+f_i({\theta }_i,{\stackrel{.}{\theta }}_i)+{\displaystyle \sum _{j=1}^{i-1}{I}_j^i}{U}_j^i+{\displaystyle \sum _{j=2}^{i-1}{\displaystyle \sum _{k=1}^{j-1}{J}_{kj}^i}}{V}_{kj}^i-{\tau }_i\right)$$

(24)

where $B_i=(1/I_{\mathrm{mi}}{\gamma }_i)\in R^{+}$ .

Define the system state vector $x_i=[\begin{array}{cc}\hfill x_{i1}\hfill & \hfill x_{i2}\hfill \end{array}{]}^T=[\begin{array}{cc}\hfill {\theta }_i\hfill & \hfill {\stackrel{·}{\theta }}_i\hfill \end{array}{]}^T$ and control input $u_i={\tau }_i$ , $i=1,2,\dots ,n$ , then the state space of ith subsystem is formulated as follows

$$S_i:\{\begin{array}{c}{\stackrel{.}{x}}_{i1}=x_{i2}\hfill \\ {\stackrel{.}{x}}_{i2}=-B_i\left(f_i(x_i,{\stackrel{.}{x}}_i)+{\displaystyle \sum _{j=1}^{i-1}{I}_j^i}{U}_j^i+{\displaystyle \sum _{j=2}^{i-1}{\displaystyle \sum _{k=1}^{j-1}{J}_{kj}^i}}{V}_{kj}^i+\frac{{\tau }_{fi}}{{\gamma }_i}\right)+B_i{u}_i\hfill \\ y=x_{i1}\hfill \end{array}$$

(25)

The control object of this article is to design a decentralized control strategy for harmonic drive–based MRR with uncertain environment contact that uses only encoder data of each joint module and to make the joints and end-effector follow the desired trajectories.

Decentralized ISMC design based on ASTA

According to the dynamic model of harmonic drive–based MRR in equation (16) and the state space in equation (25), in this section, a decentralized ISMC method is proposed based on ASTA to compensate model uncertainty in which the up-bound is unknown and to reduce controller chattering.

Let the desired position, velocity, and acceleration of the ith joint be $x_{\mathrm{id}}$ , ${\stackrel{·}{x}}_{\mathrm{id}}$ , and ${\stackrel{··}{x}}_{\mathrm{id}}$ , respectively, which are known and bounded function. Then, one can define the integral sliding surface as

$$s_i={\stackrel{·}{e}}_i(t)+k_i{e}_i(t)+z_i-{\stackrel{·}{e}}_i(t_0)-k_i{e}_i(t_0)$$

(26)

where $e_i(t)=x_{i1}(t)-x_{\mathrm{id}}(t)$ denotes the tracking position error, $e_i(t_0)$ is the position error in the initial time, $k_i$ denotes a positive constant to be determined, $z_i$ is a potential function integral variable, and its time derivative ${\stackrel{·}{z}}_i$ is defined as

$${\stackrel{·}{z}}_i=\{\begin{array}{c}\hfill \eta \sin (\pi e_i/2\eta )\left|e_i\right|<\eta \hfill \\ \multicolumn{1}{c}{\eta ,e_i\ge \eta }\\ \multicolumn{1}{c}{-\eta ,e_i\le -\eta }\end{array}$$

(27)

where $\eta$ is a determined positive constant. Obviously, for the potential function ${\stackrel{·}{z}}_i$ , when $\left|e_i\right|<\eta$ , ${\stackrel{·}{z}}_i$ is monotone increasing strictly; when $\left|e_i\right|\ge \eta$ , ${\stackrel{·}{z}}_i$ is saturated.

Remark 1

The nonlinear potential function ${\stackrel{·}{z}}_i$ is applicative for reducing the chattering effect of the sliding mode system, since the desired position error behavior is obtained by choosing a suitable $\eta$ . Compared with the conventional sliding surface design, the proposed integral sliding surface possesses stronger robustness.

Combining equations (16), (18), (21), and (22) with equation (26), the time derivative of $s_i$ is given as

$$\begin{array}{l}{\stackrel{.}{s}}_i={\stackrel{··}{e}}_i(t)+k_i{\stackrel{.}{e}}_i(t)+{\stackrel{.}{z}}_i-{\stackrel{··}{e}}_i(t_0)-k_i{\stackrel{.}{e}}_i(t_0)\\ =-B_i\left(\begin{array}{c}{\widehat{b}}_{fi}{\stackrel{.}{x}}_{i1}+\left({\widehat{f}}_{ci}+{\widehat{f}}_{si}{e}^{(-{\widehat{f}}_{\tau i}{\stackrel{.}{x}}_{i1}^2)}\right)\mathrm{sgn}({\stackrel{.}{x}}_{i1})+f_{pi}\\ +\frac{{\tau }_{fi}}{{\gamma }_i}+Y({\stackrel{.}{x}}_{i1}){\tilde{F}}_i+{\displaystyle \sum _{j=1}^{i-1}{I}_j^i}{U}_j^i+{\displaystyle \sum _{j=2}^{i-1}{\displaystyle \sum _{k=1}^{j-1}{J}_{kj}^i}}{V}_{kj}^i\end{array}\right)\\ -{\stackrel{··}{x}}_{id}+k_i\left({\stackrel{.}{e}}_i(t)-{\stackrel{.}{e}}_i(t_0)\right)-{\stackrel{··}{e}}_i(t_0)+{\stackrel{.}{z}}_i+B_i{u}_i\end{array}$$

(28)

In equation (28), we can rewrite the decentralized control law $u_i$ as $u_i=u_{i1}+u_{i2}$ to compensate the certain and uncertain terms of ${\stackrel{·}{s}}_i$ .

First, design the control law $u_{i1}$ as

$$\begin{array}{c}\multicolumn{1}{c}{u_{i1}=\frac{{\tau }_{\mathrm{fi}}}{{\gamma }_i}+{\hat{b}}_{\mathrm{fi}}{\stackrel{·}{x}}_{i1}+({\hat{f}}_{\mathrm{ci}}+{\hat{f}}_{\mathrm{si}}{e}^{(-{\hat{f}}_{\tau i}{\stackrel{·}{x}}_{i1}^2)})\mathrm{sgn}({\stackrel{·}{x}}_{i1})}\\ \multicolumn{1}{c}{-B_i^{-1}({\stackrel{·}{z}}_i+k_i({\stackrel{·}{e}}_i(t)-{\stackrel{·}{e}}_i(t_0))-{\stackrel{··}{x}}_{\mathrm{id}}-{\stackrel{··}{e}}_i(t_0))}\end{array}$$

(29)

where $u_{i1}$ is utilized for compensating the flexspline torque, modeled joint friction, and other certain terms in equation (28).

Second, design the control law $u_{i2}$ to compensate the model uncertain terms in ${\stackrel{·}{s}}_i$ . In order to provide a description of the model uncertainties supported by the ASTA, according to equation (23) and the categories of model uncertainties, one can reformulate model uncertainties as

$$\begin{array}{l}{\psi }_i=Y({\stackrel{.}{x}}_{i1}){\tilde{F}}_i+{\displaystyle \sum _{j=1}^{i-1}{I}_j^i}{U}_j^i+{\displaystyle \sum _{j=2}^{i-1}{\displaystyle \sum _{k=1}^{j-1}{J}_{kj}^i}}{V}_{kj}^i+f_{pi}\\ ={\psi }_{ic}+{\psi }_{iv}\\ {\psi }_{ic}=Y({\stackrel{.}{x}}_{i1}){\tilde{F}}_{ic}+{\displaystyle \sum _{j=1}^{i-1}{I}_j^i}{U}_{jc}^i+{\displaystyle \sum _{j=2}^{i-1}{\displaystyle \sum _{k=1}^{j-1}{J}_{kj}^i}}{V}_{kjc}^i\\ \le {\varsigma }_{i1}(x_i,t)\left|{\sigma }_{i1}(s_i)\right|\\ {\psi }_{iv}=Y({\stackrel{.}{x}}_{i1}){\tilde{F}}_{iv}+{\displaystyle \sum _{j=1}^{i-1}{I}_j^i}{U}_{jv}^i+{\displaystyle \sum _{j=2}^{i-1}{\displaystyle \sum _{k=1}^{j-1}{J}_{kj}^i}}{V}_{kjv}^i+f_{pi}\\ \le {\displaystyle \underset{0}{\overset{t}{\int }}{\varsigma }_{i2}(x_i,t)\left|{\sigma }_{i2}(s_i)\right|dt}\end{array}$$

(30)

where the continuous functions ${\varsigma }_{i1}(x_i,t)$ and ${\varsigma }_{i2}(x_i,t)$ are unknown up-bounds of model uncertainties. ${\sigma }_{i1}(s_i)$ and ${\sigma }_{i2}(s_i)$ are designed based on ASTA and represented as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\sigma }_{i1}(s_i)={\left|s_i\right|}^{\frac{1}{2}}\mathrm{sgn}(s_i)+{\kappa }_{\mathrm{ki}}{s}_i}\\ \multicolumn{1}{c}{{\sigma }_{i2}(s_i)=\frac{1}{2}\mathrm{sgn}(s_i)+\frac{3}{2}{\kappa }_{\mathrm{ki}}{\left|s_i\right|}^{\frac{1}{2}}\mathrm{sgn}(s_i)+{\kappa }_{\mathrm{ki}}^2{s}_i}\end{array}$$

(31)

where ${\kappa }_{\mathrm{ki}}$ is a determined positive constant.

Substituting equations (29) and (30) into equation (28), the time derivative of sliding surface ${\stackrel{·}{s}}_i$ is rewritten as follows

$${\stackrel{·}{s}}_i=-B_i({\psi }_{\mathrm{ic}}+{\psi }_{\mathrm{iv}}-u_{i2})$$

(32)

Then one can design the ASTA-based control law $u_{i2}$ as follows

$$u_{i2}={\epsilon }_{i1}{\sigma }_{i1}(s_i)+\underset{0}{\overset{t}{\int }}{\epsilon }_{i2}{\sigma }_{i2}(s_i)\mathrm{dt}$$

(33)

where ${\epsilon }_{i1}$ and ${\epsilon }_{i2}$ are the adaptive update laws, which are designed as follows

$$\begin{array}{c}{\stackrel{·}{\epsilon }}_{i1}=\{\begin{array}{cc}\hfill {\kappa }_{\mathrm{ai}}s{g}_i(s_i,{\kappa }_{\mathrm{uj}})\hfill & \hfill {\epsilon }_{i1}>{\kappa }_{\mathrm{vi}}\hfill \\ \multicolumn{1}{c}{{\kappa }_{\mathrm{ri}}\left|s_i\right|}& \hfill {\epsilon }_{i1}\le {\kappa }_{\mathrm{vi}}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\epsilon }}_{i2}=2{\kappa }_{\mathrm{ti}}{\stackrel{·}{\epsilon }}_{i1}}\end{array}$$

(34)

where $s{g}_i(s_i,{\kappa }_{\mathrm{uj}})=\mathrm{sgn}(s_i-{\kappa }_{\mathrm{ei}}){\left|s_i\right|}^{{\kappa }_{\mathrm{uj}}}$ , $j=1,2$ . In addition, ${\kappa }_{\mathrm{ai}},{\kappa }_{\mathrm{uj}},{\kappa }_{\mathrm{ei}},{\kappa }_{\mathrm{vi}},{\kappa }_{\mathrm{ri}}$ , and ${\kappa }_{\mathrm{ti}}$ are the positive constants to be determined.

Combining equation (29) with equation (33), the proposed decentralized control law $u_i$ for the ith joint of MRR is designed as follows

$$\begin{array}{c}\multicolumn{1}{c}{u_i=u_{i1}+u_{i2}}\\ \multicolumn{1}{c}{=({\hat{f}}_{\mathrm{ci}}+{\hat{f}}_{\mathrm{si}}{e}^{(-{\hat{f}}_{\tau i}{\stackrel{·}{x}}_{i1}^2)})\mathrm{sgn}({\stackrel{·}{x}}_{i1})}\\ \multicolumn{1}{c}{-{\stackrel{··}{x}}_{\mathrm{id}}-B_i^{-1}({\stackrel{·}{z}}_i+k_i({\stackrel{·}{e}}_i(t)-{\stackrel{·}{e}}_i(t_0))-{\stackrel{··}{e}}_i(t_0))}\\ \multicolumn{1}{c}{+\frac{{\tau }_{\mathrm{fi}}}{{\gamma }_i}+{\hat{b}}_{\mathrm{fi}}{\stackrel{·}{x}}_{i1}+{\epsilon }_{i1}{\sigma }_{i1}(s_i)+\underset{0}{\overset{t}{\int }}{\epsilon }_{i2}{\sigma }_{i2}(s_i)\mathrm{dt}}\end{array}$$

(35)

Theorem

Consider an MRR with uncertain environment contact comprising n modules, with the dynamic model of the ith subsystem as given in equation (16) and model uncertainties as in equations (19), (21), and (22) with unknown up-bounds. The trajectory tracking errors of each robotic joint are bounded under the decentralized control law designed by equation (35).

Proof

To prove the theorem, at first, we can reformulate the time derivative of the sliding surface in equation (32) as follows

$${\stackrel{·}{s}}_i=B_i({\epsilon }_{i1}{\sigma }_{i1}(s_i)-{\psi }_{\mathrm{ic}}+K_i)$$

(36)

where $K_i$ is defined as

$$K_i=-{\psi }_{\mathrm{iv}}+\underset{0}{\overset{t}{\int }}{\epsilon }_{i2}{\sigma }_{i2}(s_i)\mathrm{dt}$$

(37)

Consider the Lyapunov function candidate as the following form

$$V=V_s+V_a$$

(38)

First, define one part of Lyapunov function candidate $V_s$ as follows

$$V_s=\sum _{i=1}^n{V}_{\mathrm{si}}(s_i,K_i)=\sum _{i=1}^n{\xi }_i^{T}P{\xi }_i$$

(39)

where ${\xi }_i^T$ and P are given as

$$\begin{array}{c}\multicolumn{1}{c}{{\xi }_i^T=\left[\begin{array}{cc}\hfill {\sigma }_{i1}(s_i)\hfill & \hfill K_i\hfill \end{array}\right]}\\ \multicolumn{1}{c}{P=\left[\begin{array}{cc}\hfill p_1\hfill & \hfill p_2\hfill \\ \multicolumn{1}{c}{p_2}& \hfill p_3\hfill \end{array}\right]}\end{array}$$

(40)

where $p_1=\mathrm{diag}(4{\lambda }_i+4k_{\mathrm{ti}}^2)$ , $p_2=\mathrm{diag}(-2k_{\mathrm{ti}})$ , and $p_3=2I$ , ${\lambda }_i>0$ . Rewriting equation (30) yields $h_{\mathrm{ic}}={\vartheta }_{i1}(x_i,t){\sigma }_{i1}(s_i)$ and $h_{\mathrm{iv}}={\int }_0^t{\vartheta }_{i2}(x_i,t){\sigma }_{i2}(s_i)\mathrm{dt}$ with the functions of $|{\vartheta }_{i1}(x_i,t)|\le {\varsigma }_{i1}(x_i,t)$ and $|{\vartheta }_{i2}(x_i,t)|\le {\varsigma }_{i2}(x_i,t)$ . Then one obtains the derivative of ${\xi }_i$ as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{\xi }}_i=\left[\begin{array}{c}\hfill {\sigma \prime }_{i1}(s_i)B_i({\epsilon }_{i1}{\sigma \prime }_{i1}(s_i)-{\psi }_{\mathrm{ic}}+K_i)\hfill \\ \multicolumn{1}{c}{{\epsilon }_{i2}{\sigma \prime }_{i2}(s_i)-(d/\mathrm{dt}){\psi }_{\mathrm{iv}}}\end{array}\right]}\\ \multicolumn{1}{c}{={\sigma \prime }_{i1}(s_i)B_i\left[\begin{array}{cc}\hfill {\epsilon }_{i1}-{\vartheta }_{i1}(x_i,t)\hfill & \hfill 1\hfill \\ \multicolumn{1}{c}{{\epsilon }_{i2}-{\vartheta }_{i2}(x_i,t)}& \hfill 0\hfill \end{array}\right]{\xi }_i}\\ \multicolumn{1}{c}{={\sigma \prime }_{i1}(s_i)B_{i}A(x_i,t){\xi }_i}\end{array}$$

(41)

So that the time derivative of $V_s$ is obtained as follows

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_s=\sum _{i=1}^n{\stackrel{·}{V}}_{\mathrm{si}}(s_i,K_i)}\\ \multicolumn{1}{c}{=\sum _{i=1}^n{\sigma \prime }_{i1}(s_i)B_i{\xi }_i^T(A^T(x_i,t)P+\mathrm{PA}(x_i,t)){\xi }_i}\\ \multicolumn{1}{c}{=\sum _{i=1}^n-{\sigma \prime }_{i1}(s_i)B_i{\xi }_i^T\Omega {\xi }_i}\end{array}$$

(42)

where

$$\Omega =\left[\begin{array}{cc}\hfill {\Omega }_1\hfill & \hfill {\Omega }_2\hfill \\ \multicolumn{1}{c}{{\Omega }_2}& \hfill {\Omega }_3\hfill \end{array}\right]$$

and ${\Omega }_1$ , ${\Omega }_2$ , and ${\Omega }_3$ are represented as

$$\{\begin{array}{c}{\Omega }_1=2p_1{\epsilon }_{i1}+2p_2{\epsilon }_{i2}-2{\kappa }_{\mathrm{ti}}I\hfill \\ \multicolumn{1}{c}{{\Omega }_2=p_2{\epsilon }_{i1}+p_3{\epsilon }_{i2}-p_{1}I}\\ \multicolumn{1}{c}{{\Omega }_3=-2p_{2}I-2{\kappa }_{\mathrm{ti}}I}\end{array}$$

(43)

In equation (43), note that the following three relations are satisfied: (a) ${\Omega }_1>0$ , (b) ${\Omega }_3>0$ , and (c) ${\Omega }_1-{\Omega }_2{\Omega }_3^{-1}{\Omega }_2^T>0$ to ensure that the symmetric matrix $\Omega$ is positive definite. Then, one obtains that

$${\stackrel{·}{V}}_s\le -\sum _{i=1}^{n}2{\kappa }_{\mathrm{ti}}\sigma {\prime }_{i1}(s_i)B_i{\xi }_i^T{\xi }_i\le -\sum _{i=1}^{n}2{\kappa }_{\mathrm{ti}}\left(\frac{1}{2{\left|s_i\right|}^{\frac{1}{2}}}\right)B_i{\xi }_i^T{\xi }_i$$

(44)

Since

$$\begin{array}{c}\multicolumn{1}{c}{{\lambda }_{min}\{P\}\Vert {\xi }_i{\Vert }_2^2\le {\xi }_i^{T}P{\xi }_i}\\ \multicolumn{1}{c}{\le {\lambda }_{max}\{P\}\Vert {\xi }_i{\Vert }_2^2}\end{array}$$

(45)

where ${\lambda }_{min}\{P\}$ and ${\lambda }_{max}\{P\}$ are the minimum and maximum eigenvalues of P, respectively, and

$$\begin{array}{c}\multicolumn{1}{c}{\Vert {\xi }_i{\Vert }_2^2={\xi }_{i1}^2+{\xi }_{i2}^2}\\ \multicolumn{1}{c}{=\left|s_i\right|+2{\kappa }_{\mathrm{ki}}{\left|s_i\right|}^{\frac{3}{2}}+{\kappa }_{\mathrm{ki}}^2{s}_i^2+K_i^2}\end{array}$$

(46)

where $\Vert {\xi }_i{\Vert }_2^2$ is the Euclidean norm of ${\xi }_i$ , and

$$\left|{\xi }_{i1}\right|\le \Vert {\xi }_i{\Vert }_2\le \frac{V_{\mathrm{si}}^{\frac{1}{2}}(s_i,K_i)}{{\lambda }_{min}^{\frac{1}{2}}\{P\}}$$

(47)

From equations (44) to (47), one obtains

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_s\le -\sum _{i=1}^n\frac{2{\kappa }_{\mathrm{ti}}{B}_i{P}^{\frac{1}{2}}}{P}{V}_{\mathrm{si}}^{\frac{1}{2}}(s_i,K_i)+\frac{4{\kappa }_{\mathrm{ti}}{B}_i{\kappa }_{\mathrm{ki}}}{{\lambda }_{max}\{P\}}{V}_{\mathrm{si}}(s_i,K_i)}\\ \multicolumn{1}{c}{\le -\sum _{i=1}^n{\Gamma }_{i1}{V}_{\mathrm{si}}^{\frac{1}{2}}(s_i,K_i)-\sum _{i=1}^n{\Gamma }_{i2}{V}_{\mathrm{si}}(s_i,K_i)}\end{array}$$

(48)

where ${\Gamma }_{i1}$ and ${\Gamma }_{i2}$ are positive definites and are represented as

$${\Gamma }_{i1}=\frac{2{\kappa }_{\mathrm{ti}}{B}_i{\lambda }_{min}^{\frac{1}{2}}\{P\}}{{\lambda }_{max}\{P\}},{\Gamma }_{i2}=\frac{4{\kappa }_{\mathrm{ti}}{B}_i{\kappa }_{\mathrm{ki}}}{{\lambda }_{max}\{P\}}$$

(49)

Second, define the other part of Lyapunov function candidate $V_a$ as follows

$$\begin{array}{c}\multicolumn{1}{c}{V_a=\sum _{i=1}^n{V}_{\mathrm{ai}}({\epsilon }_{i1},{\epsilon }_{i2})}\\ \multicolumn{1}{c}{=\frac{1}{2}\sum _{i=1}^n{\tilde{\epsilon }}_{i1}^2+{\tilde{\epsilon }}_{i2}^2}\end{array}$$

(50)

where ${\tilde{\epsilon }}_{i1}={\epsilon }_{i1}-{\overline{\epsilon }}_{i1}$ and ${\tilde{\epsilon }}_{i2}={\epsilon }_{i2}-{\overline{\epsilon }}_{i2}$ . ${\overline{\epsilon }}_{i1}$ and ${\overline{\epsilon }}_{i2}$ are positive constants with the relationship of $\left|{\tilde{\epsilon }}_{i1}\right|\le {\overline{\epsilon }}_{i1}$ and $\left|{\tilde{\epsilon }}_{i2}\right|\le {\overline{\epsilon }}_{i2}$ , respectively. So that one can obtain the time derivative of $V_a$ as

$$\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{V}}_a=-\sum _{i=1}^n\left|{\tilde{\epsilon }}_{i1}\right|{\stackrel{·}{\epsilon }}_{i1}+\left|{\tilde{\epsilon }}_{i2}\right|{\stackrel{·}{\epsilon }}_{i2}}\\ \multicolumn{1}{c}{\le \sum _{i=1}^n{\varpi }_{i1}(\left|{\tilde{\epsilon }}_{i1}\right|{\stackrel{·}{\epsilon }}_{i1}+\left|{\tilde{\epsilon }}_{i2}\right|{\stackrel{·}{\epsilon }}_{i2})-\sum _{i=1}^n\left|{\tilde{\epsilon }}_{i1}\right|{\stackrel{·}{\epsilon }}_{i1}}\\ \multicolumn{1}{c}{-\sum _{i=1}^n{\varpi }_{i1}({\tilde{\epsilon }}_{i1}^2+{\tilde{\epsilon }}_{i2}^2)+\sum _{i=1}^n{\varpi }_{i2}(\left|{\tilde{\epsilon }}_{i1}\right|+\left|{\tilde{\epsilon }}_{i2}\right|)}\\ \multicolumn{1}{c}{-\sum _{i=1}^n\left|{\tilde{\epsilon }}_{i2}\right|{\stackrel{·}{\epsilon }}_{i2}-\sum _{i=1}^n{\varpi }_{i2}(\left|{\tilde{\epsilon }}_{i1}\right|+\left|{\tilde{\epsilon }}_{i2}\right|)}\end{array}$$

(51)

where ${\varpi }_{i1}$ and ${\varpi }_{i2}$ are positive constants. Combining equations (48) and (51) with ${\epsilon }_{i1}$ in equation (34), there are two cases of the following analysis. For the case of ${\epsilon }_{i1}>{\kappa }_{\mathrm{vi}}$ , one obtains the time derivative of the Lyapunov function candidate $\stackrel{·}{V}$ as follows

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{V}={\stackrel{·}{V}}_s+{\stackrel{·}{V}}_a}\\ \multicolumn{1}{c}{\le -\sum _{i=1}^n{\Gamma }_{i1}{V}_{\mathrm{si}}^{\frac{1}{2}}(s_i,K_i)-\sum _{i=1}^n{\Gamma }_{i2}{V}_{\mathrm{si}}(s_i,K_i)}\\ \multicolumn{1}{c}{-\sum _{i=1}^n{\varpi }_{i1}({\tilde{\epsilon }}_{i1}^2+{\tilde{\epsilon }}_{i2}^2)-\sum _{i=1}^n{\varpi }_{i2}(\left|{\tilde{\epsilon }}_{i1}\right|+\left|{\tilde{\epsilon }}_{i2}\right|)}\\ \multicolumn{1}{c}{-\sum _{i=1}^n\left|{\tilde{\epsilon }}_{i1}\right|({\kappa }_{\mathrm{ai}}(s{g}_i(s_i,{\kappa }_{u1})+s{g}_i(s_i,{\kappa }_{u2}))-{\varpi }_{i2}-{\varpi }_{i1}{\overline{\epsilon }}_{i1})}\\ \multicolumn{1}{c}{-\sum _{i=1}^n\left|{\tilde{\epsilon }}_{i2}\right|(2{\kappa }_{\mathrm{ti}}{\kappa }_{\mathrm{ai}}(s{g}_i(s_i,{\kappa }_{u1})+s{g}_i(s_i,{\kappa }_{u2}))-{\varpi }_{i2}-{\varpi }_{i1}{\overline{\epsilon }}_{i2})}\end{array}$$

(52)

According to equations (38) and (52), one can conclude the following inequality relationship

$$\stackrel{·}{V}+{\omega }_{i1}V+{\omega }_{i2}{V}^{1/2}\le 0$$

(53)

where ${\omega }_{i1}=min({\Gamma }_{i1},{\varpi }_{i1})$ and ${\omega }_{i2}=min({\Gamma }_{i2},{\varpi }_{i2})$ .

For the case of ${\epsilon }_{i1}\le {\kappa }_{\mathrm{vi}}$ , one can obtain that

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{V}\le -\sum _{i=1}^n{\Gamma }_{i1}{V}_{\mathrm{si}}^{\frac{1}{2}}(s_i,K_i)-\sum _{i=1}^n{\Gamma }_{i2}{V}_{\mathrm{si}}(s_i,K_i)}\\ \multicolumn{1}{c}{-\sum _{i=1}^n\left|{\tilde{\epsilon }}_{i1}\right|{\kappa }_{\mathrm{ri}}\left|s_i\right|+2\left|{\tilde{\epsilon }}_{i2}\right|{\kappa }_{\mathrm{ti}}{\kappa }_{\mathrm{ri}}\left|s_i\right|}\\ \multicolumn{1}{c}{\le 0}\end{array}$$

(54)

Therefore, according to equations (53) and (54), on the basis of Lyapunov theory, the MRR system is stable. This concludes the proof of the theorem.

Simulations

Simulation setup

To verify the advantage of the proposed decentralized control method, in this section, two 2-DOF MRRs with uncertain environment contact are used to conduct the simulations. The simulation setup of the MRRs is illustrated in Figure 1, where l₁, l₂, and l₃ denote the length of the links, respectively, and l₄ denotes the length between the base module and the time-varying constraint. By referring to our previous work, ¹⁹ the friction model parameters, which are used to conduct the simulations, are given as follows

$$\begin{array}{c}\multicolumn{1}{c}{f_{\mathrm{ci}}=3.5+0.7\sin (10{\theta }_i)\mathrm{N}\mathrm{m}}\\ \multicolumn{1}{c}{f_{\mathrm{si}}=5+\sin (10{\theta }_i)\mathrm{N}\mathrm{m}}\\ \multicolumn{1}{c}{f_{\tau i}=100+20\sin (10{\theta }_i){\mathrm{s}}^2/\mathrm{ra}{\mathrm{d}}^2}\\ \multicolumn{1}{c}{b_{\mathrm{fi}}=1.5+0.3\sin (10{\theta }_i)\mathrm{N}\mathrm{ms}/\mathrm{rad}}\end{array}$$

and the other parameters of the dynamic model and controller are given in Tables 1 and 2, respectively. The initial positions of the two MRRs are chosen as $x_1(0)=0$ and $x_2(0)=-2$ , and the desired trajectories of the two configurations are defined as follows:

• Configuration A

  $$\begin{array}{l}{x}_{1d}=0.5\cos (t)+0.2\sin (3t)\\ x_{2d}=\arcsin \left(\frac{l_4}{l_2}\sin \left(\alpha (t)\right)-\frac{l_1}{l_2}\sin \left(\alpha (t)\right)-\frac{l_3}{l_2}\sin \left(\alpha (t)-\frac{\pi }{2}\right)\right)-\alpha (t)\end{array}$$

• Configuration B

  $$\begin{array}{c}\multicolumn{1}{c}{x_{1d}=0.1\cos (4t)+0.2\sin (3t)}\\ \multicolumn{1}{c}{x_{2d}=0.2\cos (t)+0.3\sin (2t)}\end{array}$$

Figure 1.

Simulation setup of the MRRs with external environment constraint and collision: (a) configuration A and (b) configuration B.

Table 2.

Parameters of the controller.

Name	Value	Name	Value	Name	Value
${\hat{b}}_{\mathrm{fi}}$	1.2	$k_i$	100	${\kappa }_{\mathrm{ti}}$	0.2
${\hat{f}}_{\mathrm{si}}$	4	${\kappa }_{\mathrm{ki}}$	1	${\kappa }_{u1}$	1.6
${\hat{f}}_{\tau i}$	80	${\kappa }_{\mathrm{ri}}$	0.6	${\kappa }_{u2}$	0.6
${\hat{f}}_{\mathrm{ci}}$	3.0	${\kappa }_{\mathrm{vi}}$	0.01	${\kappa }_{\mathrm{ai}}$	0.02
${\lambda }_i$	1	$\eta$	0.05	${\kappa }_{\mathrm{ei}}$	0.1

Two types of external environment contacts are considered in the simulations including continuous time-varying environment constraint (configuration A) and collision at random time point (configuration B). The environment constraint of configuration A is considered as a kind of column, with length of 0.5 m and mass of 4 kg, and rotating as a time-varying angle variable

$$\alpha (t)=0.75\pi +0.2\sin \left(\frac{t}{2}\right)$$

For configuration B, the constant collision force of 50 N is upright exerted at the end-effector at the time points of 5 and 6 s. The implemented time-varying constrained force and collision force curves are given in Figure 2. Note that the contact force and the collision force information are used only for setting up the the simulations but the control designs, so that for the controllers of the MRR systems, the environment contact forces are unknown all the time.

Figure 2.

Environment contact force curves: (a) time-varying constraint force of configuration A and (b) collision force of configuration B.

Simulation results

The simulation results are obtained by employing the constructed simulation setup of the MRRs with different configurations to investigate the variation in position tracking, position error, flexspline torque, and control torque of each joint module, as well as the performance of position tracking of the end-effectors. Two kinds of control method including the conventional super-twisting algorithm (STA)-based sliding mode control (SMC)^36,37 and the proposed ASTA-based ISMC are implemented to compare the control performance and to verify the advance of the proposed method in this article.

Joint position tracking and position error

Figures 3–6 illustrate the joint trajectory tracking curves and joint position error curves under different MRR configurations using STA-based SMC and ASTA-based ISMC. From Figures 3 and 4, we obtain that the desired trajectories are tracked effectively with the time-varying constraint under both STA-based SMC and ASTA-based ISMC. However, when the collision force is exerted at the end-effector of the MRR, a 0.2-s time interval is required for STA-based SMC to restabilize the robot joints. Besides, the robot joints can be restabilized within only 0.1 s under ASTA-based ISMC. From Figures 5(a) and 6(a), we observe that the convergence of the position error is achieved under conventional STA-based SMC, and the maximum position error values of each joint are close to 4.0 e−3 rad. Besides, since the control gains and parameters cannot be adjusted adaptively when the up-bound of the uncertainty is unknown, the periodic position error is relatively obvious under this control. From the simulation results in Figures 5(b) and 6(b), one can conclude that the tracking performances of joint positions are improved obviously since model uncertainties have been compensated accurately under ASTA-based ISMC, and the sliding manifold is reached in less than 0.2 s; after that, the MRR system is in sliding manifold.

Figure 3.

Joint trajectory tracking curves of configuration A: (a) STA-based SMC and (b) ASTA-based ISMC.

Figure 4.

Joint trajectory tracking curves of configuration B: (a) STA-based SMC and (b) ASTA-based ISMC.

Figure 5.

Joint position error curves of configuration A: (a) STA-based SMC and (b) ASTA-based ISMC.

Figure 6.

Joint position error curves of configuration B: (a) STA-based SMC and (b) ASTA-based ISMC.

Flexspline torque

Figures 7 and 8 show the flexspline torque curves of each joint module, which are estimated using the motor-side and link-side position measurements. From the simulation results, one obtains that the flexspline torques of each joint are estimated appropriately under the situations of external environment contact and collision, and this may eliminate the needs of joint torque sensors, which can hardly withstand the instantaneous collisions. However, in Figures 7(a) and 8(a), one observes that the estimated flexspline torques are with the chattering effects obviously under STA-based SMC, because the chattered position measurements are employed in the harmonic drive model to estimate the flexspline torques. Moreover, improved simulation results of the joint flexspline torques are illustrated in Figures 7(b) and 8(b), since the chattering effects of the flexspline torques are reduced obviously under ASTA-based ISMC. Besides, the joint flexspline torque curves are smoother than the former ones with unknown environment contact, especially in the intervals of the initial time and collision time.

Figure 7.

Flexspline torque curves of configuration A: (a) STA-based SMC and (b) ASTA-based ISMC.

Figure 8.

Flexspline torque curves of configuration B: (a) STA-based SMC and (b) ASTA-based ISMC.

Control torque

Figures 9 and 10 show the control torque curves of each joint under STA-based SMC and ASTA-based ISMC, respectively. In Figures 9(a) and 10(a), one concludes that the desired sliding surface is reached within 0.2 s under STA-based SMC, and the chattering effects of the control torques are obvious relatively because the switching nonideality causes a high-frequency oscillation to appear in the neighborhood of the sliding manifold; moreover, another source of controller chattering is that the flexspline torques, which are used for controller design, are also with chattering effects. Besides, in Figures 9(b) and 10(b), one observes that the desired sliding surface is reached in the initial time and the control torque curves are smoother obviously than the former ones since the chattering effects are reduced effectively under ASTA-based ISMC.

Figure 9.

Control torque curves of configuration A: (a) STA-based SMC and (b) ASTA-based ISMC.

Figure 10.

Control torque curves of configuration B: (a) STA-based SMC and (b) ASTA-based ISMC.

Three-dimensional end-effector trajectory tracking

Figure 11 shows the three-dimensional (3D) end-effector trajectory tracking curves under different configurations of MRR with dynamic contact and collision under ASTA-based ISMC. From the simulation results, one obtains that when the external environment contact and collision force are exerted on the end-effector of an MRR with different configurations, the end-effector desired trajectories can be tracked accurately under the proposed ASTA-based ISMC.

Figure 11.

3D end-effector trajectory tracking curves under ASTA-based ISMC: (a) configuration A and (b) configuration B.

From the aforementioned simulation results, we can conclude that the proposed decentralized controller can provide accuracy and stability for MRRs to satisfy the requirements of various tasks with complex working environments.

Conclusion

This article focuses on investigating harmonic drive–based MRRs with uncertain environment contact and addresses the problems of trajectory tracking control using only encoder data. The dynamic model of an MRR is formulated based on a control-oriented harmonic drive model, and the magnitudes of the interconnected joint couplings are decreased significantly in this model. A decentralized control strategy is proposed based on ISMC and ASTA to compensate model uncertainty in which the up-bound is unknown and to reduce the chattering effect. The stability of the MRR system is proved by Lyapunov theory. Finally, simulations are performed to verify the advantage of the proposed method.