Robust approximate constraint following control design for collaborative robots system and experimental validation

Abstract

The paper presents a novel control method aimed at enhancing the trajectory tracking accuracy of two-link mechanical systems, particularly nonlinear systems that incorporate uncertainties such as time-varying parameters and external disturbances. Leveraging the Udwadia–Kalaba equation, the algorithm employs the desired system trajectory as a servo constraint. First, the system’s constraints to construct its dynamic equation and apply generalized constraints from the constraint equation to an unconstrained system. Second, we design a robust approximate constraint tracking controller for manipulator control and establish its stability using Lyapunov’s law. Finally, we numerically simulate and experimentally validate the controller on a collaborative platform using model-based design methods.

1. Introduction

With the development of society, collaborative robots are playing an increasingly important role in the industrialization process. They have many advantages: lightweight and agile, strong perception capabilities, easy to program, user-friendly, etc., effectively achieving collaboration between humans and machines. Collaborative robots can help people engage in some repetitive, dangerous work, such as machining, spraying, and welding. Collaborative robots can be found not only in the industrial field but also in the fields of healthcare, education, and household. [1, 2].

Considering the use of collaborative robots, it is important to increase speed and accuracy. Collaborative robot systems have multiple inputs and outputs, and these variables are coupled with each other, making the collaborative robot system a complex nonlinear system. Moreover, deviations within the collaborative robot and disturbances in the environment can lead to system uncertainty, which poses a significant challenge to the design of the controller. Traditional control methods (such as proportional–integral–derivative (PID) control) are difficult to meet the control requirements [3, 4]; therefore, we must consider designing an effective nonlinear control algorithm.

Under such conditions, we propose a robust tracking servo constraint control algorithm [5]. This control method is to control the robotic manipulator system using a dynamic equation-based approach. The concept of constrained motion and its solution using Lagrange multipliers was initially introduced by Lagrange. As a result, numerous scholars have put forward various approaches to handle constrained systems, including the principles of Maggi, Boltzmann and Hamel, Gibbs and Appell, among others [6, 7]. Nevertheless, in practice, energy-based Lagrangian methods can rely on either the virtual displacement principle or D’Alembert’s principle, and to some extent, these two principles are equivalent in their application. Hence, all these approaches encounter a common limitation: while D’Alembert’s principle is applicable in most cases, it fails to address nonideal constraints.

Not until 1992 did Udwadia et al. introduce the Udwadia–Kalaba (U–K) equations as a solution, which expanded upon the principle of D’Alembert and the Lagrange equation [8–10]. They proposed a clear and concise set of discrete equations for dynamic systems. Since the motion equations follow the Gauss principle, this gave us a new perspective – that the Gauss principle can be used to handle constraints. UK combined the principle of D’Alembert with their previous research and concluded that the work done by constraint forces under virtual displacements is zero [11, 12]. In addition, they also extended the application range of the equations, making them suitable for constrained mechanical systems that do not follow the principle of D’Alembert, such as systems with nonideal constraints.

In addition to the aforementioned research efforts, there have been other notable contributions in the field of servo constraint control. For instance, researchers have explored the application of the U–K equation in various practical scenarios, providing valuable insights and theoretical foundations. Furthermore, Chen conducted further studies on adaptive robust control for systems with uncertainties [13–15]. Chen conducted an in-depth study of the servo constraint problem on the basis of the Maggi equation, and then he argued that a reasonable servo controller design can meet the required constraints [16]. Chen also conducted research in the area of adaptive robust control of uncertain systems, as stated in references [17, 18]. These investigations aimed to develop control strategies that can effectively handle system uncertainties, providing stability and performance guarantees. Another significant development came from Schutte, who delved into the theoretical analysis and simulation research of the U–K equation [19]. As a result, a state feedback controller was proposed, specifically tailored for systems with both holonomic and nonholonomic constraints. This controller proves beneficial in scenarios where these types of constraints are present in nonlinear systems. Overall, the research conducted by Udwadia, Chen, Bajodah, and Schutte, among others, has significantly advanced the understanding and application of servo constraint control. Their work has laid the groundwork for the development of effective control strategies for addressing dynamic inverse problems and trajectory tracking in nonlinear systems with constraints.

Mechanical arm systems are a common type of mechanical system, and many researchers are now looking for and exploring different control theories in various environments, hoping to achieve better dynamic control of mechanical arms. These control theories include classical PID control, derivative-enhanced PID control, neural network theory [20, 21], reinforcement learning methods [22], H-infinity control [23, 24], adaptive control [25, 26], and sliding mode theory [27, 28].

Collaborative robotic arm, as a typical effective nonlinear control system, has been a hotspot of research for many scholars [29, 30]. The U–K equation provides a new method for solving the servo constraint control problem, by which the trajectory tracking problem in nonlinear systems can be solved effectively. The notion of servo binding control was introduced by Chen [13, 31] utilizing the U–K equation. Establishing a dynamic model is the first step in designing control strategies for mechanical systems. This will help ensure that the control strategies can effectively meet constraints and handle uncertainties. Therefore, the model should accurately capture the dynamic performance of the system, especially the effects of constraints and uncertainties. Once the dynamic model is established, constraint equations can be transformed into a second-order form, making the development of control strategies more clear. At the same time, this method can also improve the accuracy and tracking capabilities of the system, enabling mechanical systems to achieve higher performance and reliability in practical applications.

This study makes three main contributions: i) first, the U–K equation is introduced, and the analytical expression of constraint forces for ideal systems is derived based on this equation. ii) Then, a novel robust controller (RC) for mechanical systems with uncertainties is proposed and proven to be stable. iii) Finally, we used this controller to resolve the trajectory tracking accuracy issue of the collaborative robots.

2. Dynamic model and constraint analysis

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) [32]:

(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 [14], we can now obtain the constraint force.

Remark 1. When creating a dynamic model, the state variables of the system must adhere to specific constraints as outlined in the constraint equations. These typically include holonomic constraints, rheonomic constraints, nonholonomic constraints, and other constraints. The constraint forces and active forces jointly influence the motion of the system.

3. Robust approximate constraint tracking control with uncertain parameters

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.

Assumption 1. For each $(p,t)\in R^n \times R$ , $\sigma \in \varSigma$ , $H(p,\sigma, t)\gt 0$ .

Definition 1. If equation ( 19 ) has at least one solution, the constraint equation ( 19 ) is said to be consistent.

Assumption 2. Constraint equation ( 23 ) is consistent.

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.

Remark 2. By utilizing the Lagrangian form of D’Alembert’s principle, the constraint force Equation (22) is identified as one among several potential alternative forces that possess the minimum norm. Additionally, these alternative forces are capable of satisfying Equation (23) as well.

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 3. $E(p,t)$ is full rank for every $(p,t)\in R^n\times R$ . Therefore, $D(p,t)D^T(p,t)$ is invertible.

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.

Remark 4. Equation (30) represents the RC that we ultimately propose, which consists of three components. The nominal control term y1, based on the U–K equation, is used to achieve the desired control of the system. By introducing y2, we can compensate for the actual initial conditions. The term y3 is used as a control measure to compensate for uncertainties in the controlled system.

Theorem 2. Assuming that Assumptions 1, 2, 3, and 4 is met, the control error arising from mechanical system uncertainty is represented by $\eta$ .

1. i UB: For any given $r\gt 0$ , there is a finite $d(r)$ so that if $\lVert \eta \!\left(p(t_0),\dot p(t_0),t_0\right)\rVert \leq r$ , then for all $t\geq t_0$ , it holds that $\lVert \eta\!\left (p(t),\dot p(t),t\right)\rVert \leq d(r)$ .

2. ii UUB: There exists a $\underline{d}\gt 0$ such that $\lVert \eta \!\left(p(t),\dot p(t),t\right)\rVert \leq d(r)$ for every $\bar{d}\gt \underline{d}$ as $t\geq t_0 T(\bar{d},r)$ , where $T(\bar{d},r)\lt \infty$ occurs for any $r\gt 0$ with $\lVert \eta\!\left(p(t_0),\dot p(t_0),t_0\right)\rVert \leq r$ . In addition, $\bar{d}\rightarrow 0$ as $\epsilon \rightarrow 0$ .

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. Simulation and experimental validation

During the numerical simulation and experimental validation stages, we compared our algorithm with other controllers to verify its effectiveness. We evaluated the dynamic tracking performance and steady-state error of the robotic arm to highlight the superiority of our controller. In the comparison process, we selected the classical PID controller as well as our proposed RC and its components (y1 (RC1) and y1+y2 (RC12)). Additionally, to validate the robustness of this control algorithm, we applied different loads to the robotic arm to observe the performance of the controller.

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.

5. Conclusions

The study proposes a novel robust control algorithm for mechanical systems with uncertainty and validates its stability and practicality through experiments on a two-axis arm platform.

The controller consists of three components. The initial component, y1, represents the fundamental aspect of the RC control algorithm. It is the analytical form of the nominal constraint force of the system, obtained using the UK equation with ideal constraints. The second component, y2, handles the incompatibility issue under initial conditions. The final component, y3, is presented to compensate for the affect of vulnerability on the framework when considering bounded instability. After designing the controller, the Lyapunov method is used to prove that the system is uniformly bounded and ultimately bounded.