2.1. The two-link manipulator dynamic model

As shown in Fig. 1, it is a physical picture of the collaborative manipulator, which has three degrees of freedom (three rotating joints), each joint is mainly controlled by the internal frameless torque motor and servo drive. Attached to the base is the first joint with no limit of rotation angle, followed by the second and third joints. Considering the difficulty of algorithm verification and the importance of affecting motion accuracy in practical applications, we selectively fixed the first joint and only considered the second and third joints.

Figure 1.

The collaborative robot plotform.

Figure 2 is a simplified collaborative robotic arm modeling where $m_1$ and $m_1$ are the masses of link 2 and link 3, $m_3$ is the mass of the second joint module, and $m_4$ is the mass of the third joint joint module motor. $l_1$ and $l_2$ represent the lengths of links two and three, respectively, and $l_{c1}$ and $l_{c2}$ are the lengths from the center of mass to the center of rotation of each joint. alpha1 is the angle between joint 2 and the X-axis, and alpha2 is the angle between joint 3 and joint 2.

Figure 2.

The simplified model of two-DOF collaborative robot.

A simplified model of a robotic arm has been established, and we now proceed to introduce the Lagrange energy method for conducting dynamic modeling of the robotic arm. The expression for a general mechanical system is as follows:

(1) \begin{align} \frac{d}{dt}\left(\frac{\partial L (p, \dot{p}, t)}{\partial \dot{p}_{i}}\right)-\frac{\partial L (p, \dot{p}, t)}{\partial p_{i}}=Q \nonumber\\[-25pt] \end{align}

(2) \begin{align} Q=\tau _i(t) -\tau _{fi}(t), \quad i=1,2, \cdots, n \end{align}

In the above equation, $L(p, \dot{p}, t)$ is the Lagrangian operator, which is defined as the difference between the kinetic and potential energy of a mechanical system. $Q$ corresponds to the generalized force, where is the difference between the input torque and the external disturbance torque:

(3) \begin{align} L(p, \dot{p}, t)=T(p, \dot{p}, t)-V(p, t) \end{align}

and $T(p, \dot{p}, t)$ means for system kinetic energy and $V(p, t)$ means for system potential energy. $p_i \in R^n$ is the generalized coordinates of the system and $\dot{p}_i \in R^n$ is the generalized velocity of the system. $\tau _i \in R^n$ and $\tau _{fi} \in R^n$ mean the input torque and the external interference torque, respectively. Because the motor of joint 2 is built at the origin of the coordinate system, $m_3$ can be ignored in the derivation of the equation. By combining equation (1), equation (2), and equation (3), we can get the dynamics model of the collaborative manipulator:

(4) \begin{equation} \begin{aligned} \tau (t) = &H(p,t)\ddot p + C(p, \dot p, t)\dot p + F(\dot p, t)+ G(p, t), \;\;\;{\mkern 1mu} p \in{{\rm{R}}^n} \end{aligned} \end{equation}

$t$ denotes time and $p=[p_1,p_2]^T$ is the angular displacement of two joints and also an output parameter of the system. $\dot{p}=[\dot{p_1},\dot{p_2}]^T$ and $\ddot{p}=[\ddot{p_1},\ddot{p_2}]^T$ represent the angular acceleration and angular velocity of the joints, respectively. $\varepsilon$ is an uncertain term. $\tau (t)$ represents the joint torque of the two-axis arm, and $H({\cdot})\in{\mathrm{R}^n}$ represents the inertia matrix. Matrix $C({\cdot})$ represents the Coriolis and centrifugal force terms, while matrix $C({\cdot})$ represents frictional forces and external disturbances. Matrix $G({\cdot})$ represents the gravitational terms.

During continuous operation of collaborative robots, friction becomes a critical factor that cannot be overlooked. Inaccurately identifying the friction force can cause structural damage, joint jitter, and instability in the end effector’s motion and accuracy. Therefore, it is essential to precisely calculate the friction force as the robot moves. Joint friction performance, however, is especially complicated. Most extant literature assumes that friction is a function of velocity. However, after performing a parameter identification analysis, we discovered that Coulomb friction force and viscous friction force have the biggest effects on the movement of collaborative robot joints. As a result, we can derive a specific model of force of friction as follows:

(5) \begin{align} {F_i}(\dot p ) ={\delta _{vi}}{\dot p _i} +{\delta _{ci}}sign\left ({{{\dot p }_i}} \right ),\;\;\;{\mkern 1mu} i = 1, \ldots n \end{align}

Here, $\delta _{vi}$ and $\delta _{ci}$ represent the viscous friction matrices and the coefficients of Coulomb friction, respectively.

The conventional friction models are frequently irregular and discontinuous in character, and Coulomb friction is frequently represented using a symbolic function. In order to solve this issue, the sign function is approximated using the function of hyperbolic tangent, as shown in equation (6) [Reference Kovincic, Müller, Gattringer, Weyrer, Schlotzhauer and Brandstötter32]:

(6) \begin{align} sign(\dot p _i)\approx tanh(\gamma \dot p _i), \quad i=1,\ldots n \end{align}

To ensure that the function of hyperbolic tangent closely approximates the sign function, the value of “ $ \gamma$ ” must be chosen to be sufficiently large.

The Lagrange equation can be further expanded:

(7) \begin{align} H(p,t) = \left[ \begin{array}{l@{\quad}l} H_{11} & H_{12}\\[3pt] H_{21} & H_{22}\end{array}\right] \end{align}

where

(8) \begin{equation} \left \lbrace \begin{aligned} H_{11}=\,&I_1+I_2+(m_2+m_4)l_2^2+m_3l_1^2+m_2l_{c2}^2+m_4l_2^2+2(m_2l_1l_{c2}+m_4l_1l_2)cosp_2\\[3pt] H_{12}=\,&H_{21}=I_2+m_2l_{c2}^2+m_4l_2^2 +(m_2l_1l_{c2}+m_4l_1l_2)cosp_2\\[3pt] H_{22}=\,&I_2+m_2l_{c2}^2+m_4l_2^2 \end{aligned} \right . \end{equation}

and

(9) \begin{align} C(p,\dot{p},t) = \left[ \begin{array}{l@{\quad}l} C_{11} & C_{12}\\[3pt]C_{21}&C_{22}\end{array}\right] \end{align}

where

(10) \begin{equation} \left \lbrace \begin{aligned} C_{11}=&-2(m_2l_1l_{c2}+m_4l_1l_2)sinp_2\ddot p_2\\[3pt] C_{12}=&-(m_2l_1l_{c2}+m_4l_1l_2)sinp_2\ddot p_2\\[3pt] C_{21}=&(m_2l_1l_{c2}+m_4l_1l_2)sinp_2\ddot p_1\\[3pt] C_{22}=&0. \end{aligned} \right . \end{equation}

and

(11) \begin{align} G(p,\dot{p},t)= \left[\begin{array} {l} G_{11}\\[3pt] G_{12}\end{array} \right] \nonumber\\[-25pt] \end{align}

(12) \begin{align} \left \lbrace \begin{aligned} G_{21}=&(m_2l_1l_{c2}+m_4l_1l_2)sinp_2\ddot p_1\\[3pt] G_{22}=&0. \end{aligned} \right . \end{align}

In the actual working environment, it is difficult to obtain accurate system parameters, and it is also necessary to consider that the system will be affected by some uncertain factors such as mechanical structure and unknown disturbance. Therefore, we introduce the parameter $\sigma \in R^2$ to describe these uncertainties, and the kinetic equation can be written as:

(13) \begin{align} H(p, \sigma, t)\ddot{p}(t)+C(p, \dot{p}, \sigma, t)\dot{p}(t)+F(p, \dot{p}, \sigma, t)= \tau (t) \end{align}

2.2. Constraint equations in mechanical systems

In order to simplify the analysis process, we first consider the constraint equation of the nondeterministic mechanical system, that is, the fully smooth, ideal state. According to the U–K equation:

(14) \begin{align} H(p,t)\ddot{p}=U(p,\dot{p},t) \end{align}

In the above equation, $H(p,t)\in R^n$ represents the symmetric positive definite inertia matrix, and $U(p,\dot{p},t)\in R^n$ stand for the generalized active moment matrix which acts on the system to release constraint forces. Then, we can get

(15) \begin{align} \ddot{p}=H^{-1}U(p,\dot{p},t)=\alpha (p,\dot{p},t) \end{align}

The above expression is the generalized acceleration expression for unconstrained systems. Based on the UK equations and without considering friction, the motion equations of a general unconstrained system without uncertainty can be derived. The representation of these control parameters can be expressed in matrix form as follows:

(16) \begin{align} E(p,t)=a(p,t) \end{align}

where $E=[E_{ji}]_{m\times n}$ and $a=[a_1,a_2,\ldots, a_m]^{m\times n}$ . We first take the derivative of the time t with respect to the constraint equation, which can give us its second-order form. Then we analyze the constraint equation so that we can obtain the performance of the second-order constraint form. The specific process is as follows:

Consider the servo restriction as follows:

(17) \begin{align} \sum \limits _{i=1}^{n} E_{ji}(p,t)=a_j(p,t),\quad j=1,\ldots, m \end{align}

We take the derivative of equation (17) with respect to time t and get:

(18) \begin{align} \sum \limits _{i=1}^{n} \frac{d}{dt} E_{ji}(p,t)=\frac{d}{dt} a_j(p,t) \end{align}

We rewrite Eq. (18) as:

(19) \begin{align} \sum \limits _{i=1}^{n} D_{ji}(p,t)\dot p=d_j(p,t) \end{align}

Eq. (19) provides the constraint’s first-order form. Differentiating Eq. (19) with respect to time, we obtain the following equation:

(20) \begin{align} \sum \limits _{i=1}^{n} \left(\frac{d}{dt}D_{ji}(p,t)\right)\dot p_i+\sum \limits _{i=1}^{n} D_{ji}(p,t)\ddot p_i=\frac{d}{dt} d_j(p,t) \end{align}

Rewrite Eq. (20) as:

(21) \begin{align} \sum \limits _{i=1}^{n} D_{ji}(p,t)\ddot p_i=c_j(p,\dot p,t) \end{align}

The second-order form of the constraint is represented by Eq. (21). The constraint equation is then stated as follows after this equation:

(22) \begin{align} D(p,t)\dot p=d(p,t) \nonumber\\[-25pt] \end{align}

(23) \begin{align} D(p,\dot p,t)\ddot p=c(p,\dot p,t) \end{align}

where $A=[A_{ji}]_{m\times n}$ , $d=[d_1,d_2,\ldots, d_m]^T$ , and $c=[c_1,c_2,\ldots, c_m]^T$ . When the uncertainty is known by [Reference Chen14], we can now obtain the constraint force.

3.1. In the absence of uncertainty, constraint force

When there is no uncertainty in the system and the initial conditions are compatible, we analyze the constraints within the system.

Theorem 1. Taking into account both the dynamics equation (13) and the constraint equation (23), the expression for the constraint force of the system can be derived under the Assumption 1 and 2 as follows:

(24) \begin{gather} U^c(p,\dot p,t)=\bar{H}^{1/2}(p,t)(A(p,t)\bar{H}^{-1/2}(p,t))^+[c(p,\dot p,t) \notag \\[3pt] +D(p,t)\bar{H}^{-1}(p,t)(\bar{C}(p,\dot p,t)\dot p+\bar{G}(p,t))+\bar{F}(p,t))] \end{gather}

In the above equation, the symbol $\bar{({\cdot})}$ represents a nominal term, and the symbol $()^{+}$ represents the Moore–Penrose generalized inverse. This constraint equation constitutes the first part of the controller.

3.2. Robust servo control design

In practical mechanical systems, we need to consider the uncertainty of the controlled object. We divide the terms in the dynamic equation into deterministic and stochastic parts, as shown below:

(25) \begin{gather} H(p,\sigma, t)=\bar{H}(p,t)+\varDelta H(p,\sigma, t) \notag \\[3pt] U(p,\dot p,\sigma, t)=\bar{U}(p,\dot p,t)+\varDelta U(p,\dot p,\sigma, t) \notag \\[3pt] U^c(p,\dot p,\sigma, t)=\bar{U^c}(p,\dot p,t)+\varDelta U^c(p,\dot p,\sigma, t) \end{gather}

where $\bar{H}, \bar{U}, \bar{U^c}$ are the nominal portions, $\varDelta \bar{H}, \varDelta \bar{U}, \varDelta \bar{U^c}$ stand for the uncertain parts in the constraint system. And $\bar{H}, \bar{U}, \bar{U^c}$ and $\varDelta{H}, \varDelta{U}, \varDelta{U^c}$ are continuous. We let

(26) \begin{gather} E(p,t)=\bar{H}^{-1}(p,t) \notag \\[3pt] \varDelta E(p,\sigma, t)=H^{-1}(p,\sigma, t)-\bar{H}^{-1}(p,t) \notag \\[3pt] K(p,\sigma, t)=\bar{H}(p,t)H^{-1}(p,\sigma, t)-I \end{gather}

Therefore,

(27) \begin{gather} \varDelta E(p,\sigma, t)=E(p,t)K(p,\sigma, t) \notag \\[3pt] H^{-1}(p,\sigma, t)=E(p,t)+\varDelta E(p,\sigma, t) \end{gather}

Assumption 4. Considering the condition stated in Assumption 3 , let $P\in R^{n\times n}$ be a positive definite matrix:

(28) \begin{gather} W(p,\sigma, t)=PD(p,t)E(p,t)K(p,\sigma, t)\bar{H}(p,t)D^T(p,t)(D(p,t)D^T(p,t))^{-1}P^{-1} \end{gather}

For every $(p,t)\in R^n\times R$ , $\hat{\rho }_E(\bullet ): R^n\times R\longrightarrow ({-}1,\infty )$ is real:

(29) \begin{align} \frac{1}{2}\mathop{min}\limits _{\sigma \in \sum }\lambda _m\!\left(W(p,\sigma, t)+W^T(p,\sigma, t)\right)\geq \hat{\rho }_E(p,t) \end{align}

where $min(\bullet )$ represents the smallest integer and $\lambda _m(\bullet )$ represents the smallest eigenvalue.

Remark 3. Due to the unknown nature of the uncertainty bound, the constant $\hat{\rho }_E(\bullet )$ remains unknown as well. In the absence of uncertainty in the system, where $\bar H$ is equal to $H$ , $K$ is equal to 0, and $W$ is equal to 0, it is possible to choose $\hat{\rho }_E(\bullet )$ = 0. Therefore, by continuity, this assumption implies that the presence of uncertainty affects the potential deviations of $H$ from $\bar H$ within a specific threshold. It is crucial to highlight that this threshold is unidirectional, indicating that there are no restrictions or limitations in one particular direction.

Based on the analysis and assumptions mentioned above, we propose the following RC:

(30) \begin{gather} \tau (t)=y_1(p(t),\dot p(t),t)+y_2(p(t),\dot p(t),t)+y_3(p(t),\dot p(t),t) \end{gather}

where

(31) \begin{gather} y_1\!\left(p(t),\dot p(t),t\right)=\bar{H}^{1/2}(p,t)\left(A(p,t)\bar{H}^{-1/2}(p,t)\right)^+\!\left[c(p,\dot p,t) \right. \notag \\[3pt] \left. +D(p,t)\bar{H}^{-1}(p,t)\left(\bar{C}(p,\dot p,t)\dot p+\bar{G}(p,t)\right)+\bar{F}(p,t))\right] \end{gather}

that is, $y_1(p,\dot p,t)=Q^c(p,\dot p,t)$ .

(32) \begin{gather} y_2\!\left(p(t),\dot p(t),t\right)=-k\bar{H}(p,t)D^T(p,t)\left(D(p,t)D^T(p,t)\right)^{-1} \notag \\[3pt] P^{-1}\eta (p,\dot p,t) \end{gather}

(33) \begin{gather} y_3\!\left(p(t),\dot p(t),t\right)= {-} \left[\bar{H}(p,t)D^T(p,t)\left(D(p,t)D^T(p,t)\right)^{-1}P^{-1}\right] \notag \\[3pt] \gamma (p,\dot p,t)\delta (p,\dot p,t)\rho (p,\dot p,t) \end{gather}

where $\epsilon, k\gt 0$ ,

(34) \begin{align} \gamma (p,\dot p,t)=\left \{ \begin{array}{rcl} \dfrac{(1+\hat{\rho }_E(p,t))^{-1}}{\lVert \delta (p,\dot p,t)\rVert } & &{\lVert \delta (p,\dot p,t)\rVert \gt \epsilon }\\[3pt] \dfrac{(1+\hat{\rho }_E(p,t))^{-1}}{\epsilon } & &{\lVert \delta (p,\dot p,t)\rVert \leq \epsilon }\\[3pt] \end{array} \right . \nonumber\\[-25pt]\end{align}

(35) \begin{align} \delta (p,\dot p,t)&=\eta (p,\dot p,t)\rho (p,\dot p,t), \notag \\[3pt] \eta (p,\dot p,t)&=D(p,t)\dot p-d(p,t) \end{align}

The function $\rho (\bullet ):R^n\times R^n\times R\longrightarrow R_+$ is seleced so that

(36) \begin{align} \rho (p,\dot p,t)&\geq \max \limits _{\sigma \in \varSigma }\lVert PD(p,t)\varDelta E(p,t)({-}U\dot p+y_1(p,\dot p,t) \notag \\[3pt] &+y_2(p,\dot p,t))+PD(p,t)E(p,t)\varDelta Q(p,\dot p,t))\rVert \end{align}

Eq. (28) gives the upper bound of the expression on the right.

Proof. A valid Lyapunov function candidate is constructed as demonstrated below:

(37) \begin{equation} V(\eta )=\eta ^TP\eta \end{equation}

Take the devirative of $V$ with respect to t, then we get

(38) \begin{equation} \dot{V}=2\eta ^TP(D\ddot{p}-c) \end{equation}

According to the previous analysis, we can obtain $\ddot{p}=H^{-1}(U+y_1+y_2+y_3)$ . Substituting it into the equation and combining with formulas (25) and (27), we can get

(39) \begin{align} \dot{V}&=2\eta ^TP(D\ddot{p}-c) \notag \\[3pt] &=2\eta ^TP\!\left\{D[H^{-1}(U+y_1+y_2+y_3)-b\right\} \notag \\[3pt] &=2\eta ^TP\!\left\{D\!\left[(E+\varDelta E)(U+y_1+y_2+y_3)\right]-b\right\} \notag \\[3pt] &=2\eta ^TP\!\left[DE(\bar U + y_1)-a\right] \notag \\[3pt] &=2\eta ^TPD\!\left[E\varDelta{\bar U}+\varDelta E(U+y_1+y_2)\right] \notag \\[3pt] &=2\eta ^TPDEy_2 \notag \\[3pt] &=2\eta ^TPD\!\left(E+\varDelta E\right)y_3 \notag \\[3pt] &=\sum _{i = 1}^{4} \varGamma _i(\eta, t) \end{align}

where

(40) \begin{equation} \left \lbrace \begin{aligned} \varGamma _1(\eta, t) & = 2\eta ^TP\!\left[DE(\bar{U}+y _1)-a\right] \\[3pt] \varGamma _2(\eta, t) & = 2\beta ^TPD\!\left[E\varDelta \bar{U}+\varDelta E(U+y _1+y _2)\right] \\[3pt] \varGamma _3(\eta, t) & = 2\beta ^TPDEy _2 \\[3pt] \varGamma _4(\eta, t) & = 2\beta ^TPD(E+\varDelta E)y _3 \end{aligned} \right . \end{equation}

Decompose $\dot{V}(\eta, t)$ into four parts as above and prove and analyze them separately.

When not considering uncertainty, the equation of motion of the system is

(41) \begin{equation} \bar{H}\ddot p=\bar{U}+U^c \end{equation}

where $U^c=y_1$ , that is,

(42) \begin{equation} \ddot p=\bar{H}^{-1}(\bar{U}+y_1) \end{equation}

Combined with the constraint equation Eq. (23), it can be obtained as:

(43) \begin{equation} D\!\left[\bar{HJ}^{-1}(\bar{U}+y_1)\right]-c=0 \end{equation}

that is,

(44) \begin{equation} \varGamma _1(\eta, t)=2\eta ^TP\!\left[DE(\bar{U}+y_1)-a\right]=0 \end{equation}

For $\varGamma _2(\eta, t)$ , and by Eq. (36), there is

(45) \begin{equation} \begin{aligned} \varGamma _2(\eta, t) &= 2\eta ^TPD\!\left[E\varDelta U+\varDelta E(U+y_1+y_2)\right] \\[3pt] &\leq 2 \lVert \eta \rVert \left\lVert PD\!\left[E\varDelta U + \varDelta E(U+y_1+y_2)\right]\right\rVert \\[3pt] &\leq 2\lVert \eta \rVert \rho \end{aligned} \end{equation}

Next, for $\varGamma _3(\eta, t)$ , based on Eq. (26) and Eq. (32)

(46) \begin{equation} \begin{aligned} \varGamma _3(\eta, t) &=2\eta ^TPDEy_2 \\[3pt] &=2\eta ^TPDE\!\left({-}k\bar{H}D^T(DD^T)^{-1}P^{-1}\eta \right) \\[3pt] &=-2k\eta ^T\eta \\[3pt] &=-2k\lVert \eta \rVert ^2 \end{aligned} \end{equation}

Finally, for $\varGamma _4(\beta, t)$ , by Eq. (25)–(29), Eq. (33), and Eq. (35), we can get

(47) \begin{equation} \begin{aligned} \varGamma _4(\eta, t) &=2\eta ^TPD(E+\varDelta E)y_3 \\[3pt] &=-2\eta ^TPD(\bar{H}^{-1}+EK) \left[\bar HD^T(DD^T)^{-1}P^{-1}\right]\gamma \delta \rho \\[3pt] &=-2\gamma \delta ^T\delta -2\delta ^T\!\left[PDEK\bar HD^T(DD^T)^{-1}P^{-1}\right]\gamma \delta \\[3pt] &=-2\gamma \delta ^T\delta -\gamma \delta ^T(W +W ^T)\delta \\[3pt] &\leq -2\gamma \lVert \delta \rVert ^2-\gamma \lambda _m(W +W ^T)\lVert \delta \rVert ^2 \\[3pt] &\leq -2\gamma (1+\hat{\rho }_E)\lVert \delta \rVert ^2 \end{aligned} \end{equation}

Combining Eq. (44)–(47), there is

(48) \begin{equation} \dot V(\eta, t)\leq 2\lVert \eta \rVert \rho -2k\lVert \eta \rVert ^2-2\gamma (1+\hat{\rho }_E)\lVert \delta \rVert ^2 \end{equation}

when $\lVert \delta \rVert \gt \epsilon$ , that is, $\gamma =(1+\hat{\rho }_E)^{-1}/{\lVert \delta \rVert }$

(49) \begin{equation} \begin{aligned} \dot{V}(\eta, t) &\leq -2k\lVert \eta \rVert ^2 -2\frac{(1+\hat{\rho }_E)^{-1}}{\lVert \delta \rVert }(1+\hat{\rho }_E)\lVert \delta \rVert ^2+ 2\lVert \eta \rVert \rho \\[3pt] &=-2k\lVert \eta \rVert ^2 \end{aligned} \end{equation}

when $\lVert \delta \rVert \leq \epsilon$ , that is, $\gamma =(1+\hat{\rho }_E)^{-1}/{\epsilon }$

(50) \begin{equation} \begin{aligned} \dot{V}(\eta, t) &\leq 2\lVert \eta \rVert \rho -2k\lVert \eta \rVert ^2 -2\frac{(1+\hat{\rho }_E)^{-1}}{\epsilon }(1+\hat{\rho }_E)\lVert \delta \rVert ^2 \\[3pt] &= 2\lVert \delta \rVert -2k\lVert \eta \rVert ^2 -2\frac{\lVert \delta \rVert ^2}{\epsilon }\\[3pt] &\leq \frac{\epsilon }{2} -2k\lVert \eta \rVert ^2 \end{aligned} \end{equation}

Finally, based on Eq. (49) and Eq. (50), we conclude that

(51) \begin{equation} \dot V\leq -2k\lVert \eta \rVert ^2+\epsilon /2 \end{equation}

It is clear from the preceding analysis that the controlled system demonstrates uniform boundedness:

(52) \begin{equation} d(r)=\left \lbrace \begin{aligned} r\left [\frac{\lambda _{M}(P)}{\lambda _{m}(P)}\right ]^{1/2} &\quad{r}=R \\[3pt] R\left [\frac{\lambda _{M}(P)}{\lambda _{m}(P)}\right ]^{1/2} & \quad r \leq R \\[3pt] \end{aligned} \right . \end{equation}

where

(53) \begin{equation} R= \left [ \frac{\epsilon }{4k }\right ] ^{1/2} \end{equation}

Additionally, the verification of uniform ultimate boundedness is conducted. That is,

(54) \begin{equation} \bar{d}\gt \underline{d}=R\left [\frac{\lambda _{M}(P)}{\lambda _{m}(P)}\right ]^{1/2} \end{equation}

(55) \begin{equation} T(\bar{d},r)=\left \lbrace \begin{aligned} &0 &\quad &{r}\leq \bar{d} \left [\frac{\lambda _{M}(P)}{\lambda _{m}(P)}\right ]^{1/2} \\[3pt] &\frac{\lambda _{M}(P) r^{2}-\left (\lambda _{m}^{2}(P) / \lambda _{M}(P)\right ) \bar{d}^{2}}{2 k \bar{d}^{2}\left(\lambda _{m}(P) / \lambda _{M}(P)\right)-(\epsilon / 2)} &\quad &otherwise \end{aligned} \right . \end{equation}

The appropriate design parameter $\epsilon$ can be used to predetermine the size of the ultimate boundedness ball $\bar{d}$ .According to the above analysis and derivation, the design process of RC can be summarized as the following steps:

1. 1. First, we disregard the constraint and obtain its equation of motion through Newton–Lagrange mechanics.

2. 2. Then consider the existing constraints, such as holonomic constraints, nonholonomic constraints, hardening constraints, etc., exist alone or together. Obtain the constraint equation(22).

3. 3. Write the constraint equation as second order to get equation (23) and then get the binding force based on equation (24). Then the RC is designed according to equations (31)–(33).

4.1. Parameters selection

The selection of different parameters can have different effects on the performance of the controller, so choosing appropriate parameters is an important part of the experiment. In practical working environments, the collaborative robotic arm will be affected by various factors. Therefore, it is necessary to thoroughly consider the selection of optimal parameters to compare the performance of various algorithms. Table I shows the dynamic parameters, such as mass and arm length, obtained from the three-dimensional model of the collaborative robotic arm. As discussed earlier, the F term in the dynamic equation represents frictional force. Therefore, the selection of suitable friction type and parameters can also affect the performance of the collaborative robotic arm. Based on literature review and analysis, we have listed the reference values of Coulomb coefficient and viscosity coefficient for joint friction of the collaborative robotic arm in Table II.

Table I.

The parameters and variables of a robotic manipulator system with two links.

Table II.

The initial parameters of a collaborative manipulator.

4.2. Computer simulations

In this section, we use computer modeling to simulate the robotic arm system and demonstrate the feasibility of the proposed controller through simulation results.

The structure of the two-link manipulator can be similar to the human arm. Compared with the complex multi-degree-of-freedom collaborative manipulator, the manipulator has simple structure and convenient operation. Fig. 2 shows the simplified model and overall structure of the two-degree-of-freedom manipulator. The dynamic equation of the manipulator is established using Newton–Lagrange method. The following table shows the parameters and variables of each joint of the manipulator.

We know that a typical PID control system consists of three variables: Kp, Ki, and Kd, which correspond to the proportional, integral, and derivative control parameters, respectively. They respectively affect the response speed, steady-state error, and disturbance rejection capability of the mechanical system. Therefore, for the RC we proposed, we can also control it from these three aspects. Based on the previous analysis and assumptions, p is the given parameter matrix of the RC controller, which is directly proportional to the system’s response speed. However, if it is too large, it may lead to instability and increased cost. K is the gain coefficient used to solve the problem of incompatible initial conditions. An excessively large K will reduce the tracking accuracy of the mechanical system. The parameter a represents the size of the ultimate boundedness radius, and it can be designed to improve the overall performance of the system (stability, convergence speed, and robustness). A value of a that is too small may result in system instability. In this case, the p and K parameters of the RC12 controller are the same as those of the RC controller.

Table III presents the parameter values for different controllers that exhibit relatively good control performance. In the process of parameter selection, we first adjust the parameters that have a significant impact to make the response signal follow the desired output. Then, we fine-tune the other parameters to improve tracking performance. Throughout this process, we conducted multiple experiments and adjustments, simulating different algorithms with various parameters, and eventually obtained a set of parameters that yielded a favorable overall performance.

Table III.

Selecting control parameters for a controller.

4.3. Numerical simulations and result

To validate the effectiveness of the controller, we applied a sine signal to link 1, causing it to oscillate back and forth in a regular manner. Similarly, link 2 also oscillates back and forth, with the sum of their angles satisfying the condition of being zero. Thus, it is assumed that the two-axis arm meets the following servo constraints:

(56) \begin{gather} p_1+p_2=0 \notag \\[3pt] p_1=\frac{3}{2}sin\!\left(\frac{\pi }{6}t\right) \end{gather}

By taking the second derivative of Eq. (51) with respect to time, we obtain

(57) \begin{equation} \left \lbrace \begin{aligned} \dot p_1&+\dot p_2= 0 \\[3pt] \dot p_1&= \frac{\pi }{4}cos\!\left(\frac{\pi }{6}t\right) \end{aligned} \right . \end{equation}

Taking the derivative of the above equation with respect to time t, we obtain

(58) \begin{gather} \ddot p_1+\ddot p_2=0 \notag \\[3pt] \ddot p_1=-\frac{\pi ^2}{24}sin(\frac{\pi }{6}t) \end{gather}

Once the constraint equations are determined, we can use the previous formulas to obtain the solution:

(59) \begin{equation} \begin{aligned} D= \left [\begin{array}{l@{\quad}l} 1 & 1\\[3pt] 1 & 0 \end{array}\right ], \quad d=\left [\begin{matrix} 0 \\[3pt] \frac{\pi }{4}cos(\frac{\pi }{6}t) \end{matrix}\right ], \quad c=\left [\begin{matrix} 0 \\[3pt] -\frac{\pi ^2}{24}sin(\frac{\pi }{6}t) \end{matrix}\right ] \end{aligned} \end{equation}

The computer simulation results are shown below. In this simulation, we introduce variations in the mass of the connecting rod and the motor to mimic uncertain parameters within the system. By adjusting these parameters, we can assess the robustness and adaptability of the control algorithm in the face of such uncertainties. The change rule is $ m_{actual} = m_{nominal} + 0.01cos(t)$ , which indicates that the mass of the mechanical arm changes with the change of time t. The simulation results demonstrate a significant impact of this change on the system’s stability, thereby enabling a better verification of the proposed controller’s robustness. Taking into account the situation where the initial conditions of the system are incompatible, we set the initial angles of the system to be nonzero. Joint 1 is given angle $\frac{\pi }{9}$ , and joint 2 is given angle $-\frac{\pi }{9}$ .

The figures display the simulation results. Figs. 3–6 represent the simulation result graphs. Figure 3(a) and (b) show the angular displacement variation curve and error of link 1. Fig. 4(a) and (b) depict the angular displacement variation curve and error of link 2. Figure 5 represents the sum of angular displacement variation curves for link 1 and link 2. Figure 6(a) and (b), respectively, illustrate the control torque output curves for the two joints.

Figure 3.

The experimental chart of collaborative robot joint 1.

Figure 4.

The experimental chart of collaborative robot joint 1.

Figure 5.

Sum of link 1 and link 2’s angular displacements.

Figure 6.

The experimental chart of collaborative robot joint 1.

From Figs 3–6, it can be seen that when only y1 is applied, when ignoring the issue of incompatible initial conditions and uncertain influences, the collaborative robotic arm completely fails to produce the desired tracking trajectory, and the error increases with time. However, when both y1 and y2 are applied simultaneously, the robotic arm is able to track the trajectory within a certain range of error, showing some superiority compared to the PID controller. If we also consider y1, y2, and y3 simultaneously, the collaborative robotic arm, under the control of servo constraints, achieves the desired trajectory with a tracking error less than 0.05.

According to the above analysis, the proposed controller can track the desired trajectory quickly and accurately, has good robustness, and still has small steady-state error under the interference of uncertainty. Therefore, the controller we designed is reasonable and effective.

4.4. Experimental platform introduction

As shown in Fig. 7, the experimental apparatus for this experiment consists of a load, MATLAB/Simulink software, a collaborative three-axis robotic arm, and a CSPACE control system. The collaborative robotic arm has three rotational degrees of freedom, and in this experiment, the last two axes are utilized. The joints are directly driven by motors and provide real-time joint position feedback. The graphical user interface (GUI) and the CSPACE control platform make up the CSPACE control system simultaneously. The CSPACE controller is equipped with a powerful signal processor that can be employed to validate the control performance of algorithms.

Figure 7.

Experimental platform of collaborative robot.

The design and control process of the algorithm for this control platform are shown in Fig. 8, and the approximate steps for development and validation are as follows:

• First, use MATLAB/Simulink to model the collaborative robotic arm.

• Based on the DSP development board, convert the dynamic model into C code.

• Download the compiled files to the CSPACE controller through a host computer.

• The GUI allows for real-time monitoring of the mechanical object’s dynamic performance. Control parameters can be continuously adjusted to achieve the desired results.

• Store and handle data from experiments.

Figure 8.

The flow chart for developing algorithms using CSPACE control system.

4.5. Experimental verification

In order to verify the control effectiveness of this algorithm, we choose a three-axis robot arm to fix the robot arm connected to the base and use two arms at the end to verify. The control parameters used in this part are consistent with the simulation process. According to the simulation, it can be seen that the control algorithm using y1 alone cannot meet the tracking performance and may have unpredictable risks. Therefore, in the experiment, we only compared the dynamic performance of y1y2 and y1y2y3 control algorithms and compared them with PID.

We give the system a fixed sine wave and then use different control algorithms to control the robot arm to track this trajectory, so as to verify the control performance of different algorithms. Figure 9 show the angular displacement tracking performance and error of joint 1, and Fig. 10 show the tracking performance and error of joint 2. According to Figs. 9 and10, it can be seen that under the y1y2y3 control algorithm, the tracking performance of joint 1 is the best, with an error of about $\pm 0.05$ , followed by the y1y2 algorithm, which is slightly worse than the y1y2y3 algorithm, but its error can also be maintained at about $\pm 0.1$ . Compared with the two, the error of PID algorithm is relatively large. At the same time, in order to verify the robustness of the controller, we apply different loads to the manipulator system, respectively, and make the manipulator track the desired trajectory under different loads. As can be seen from Figs. 11(a) –11(f), loads of 0 kg, 0.5 Kg, and 1.5 Kg are applied to the system, respectively. According to the experimental data, it can be seen that with the increase of load, the experimental error also increases gradually, and the fluctuation of joint 2 is larger. The results show that the controller proposed in this paper can track the trajectory well under different loads, indicating that the controller has strong robustness.

Figure 9.

The experimental chart of collaborative robot joint 1.

Figure 10.

The experimental chart of collaborative robot joint 2.

Figure 11.

Position error curves of collaborative manipulator under different loads.