Adaptive backstepping controller based on a novel framework for dynamic solution of an ankle rehabilitation spherical parallel robot

Abstract

This research offers an adaptive model-based methodology for autonomous control of 3-RRR spherical parallel manipulator (RSPM) based on a novel modeling framework. RSPM is an overconstrained parallel mechanism that has a variety of applications in medical procedures such as ankle rehabilitation because of its precision and accuracy. However, obtaining a complete explicit dynamic model of these mechanisms for tracking purposes has been a problematic challenge due to their inherent singularities, coupling effects of the limbs, and redundant constraints imposed by the intermediate joints. This paper presents a novel algorithm to obtain the analytical kinematic solutions of RSPMs based on the closed-loop vector method, which includes constraint analysis. By incorporating constrained kinematics into the dynamic model, a comprehensive explicit dynamic solution of the non-overconstrained version 3-RCC of RSPM is developed in task space, based on screw theory and the linear homogeneous property of algebraic equations on the manipulator twist. Based on the proposed computational framework, a robust self-tuning backstepping control (STBC) strategy is applied to the robot to overcome the effect of external disturbances and time-varying uncertainties. Furthermore, an observer-based compensation (OBC) method is presented for dealing with the nonlinear hysteresis loops of the ankle during trajectory tracking purposes. The closed-loop stability of the whole system including STBC and OBC is theoretically performed by Lyapunov methods. The proposed methodologies are validated by realistic co-simulations in different scenarios. For instant, in the presence of external disturbances, the maximum tracking error norm of STBC is 37.5% less than the sliding mode approach.

1. Introduction

1.1. Ankle rehabilitation devices and matching the ankle joint complex

Parallel robots (PRs) are extensively utilized in robotic rehabilitation, specifically ankle rehabilitation. Due to their excessive rigidity and high accuracy, PRs are widely used when performing active and passive ankle rehabilitation exercises. Numerous different parallel robots, including exoskeletons, soft robots, and platform robots, have been proposed and controlled for ankle rehabilitation therapy. A general summary of the latest ankle rehabilitation robots and their properties and structural types has been presented in ref. [1–3].

The ankle joint complex (AJC) consists of the three primary bones, the tibia, fibula, and talus. The mobility of AJC can be visualized by three separate rotations, which are illustrated in Fig. 1. The rotational axes do not precisely coincide from an anatomical perspective [5], which causes residual translational motions of the foot [6]. Six generalized spherical parallel manipulators (GSPMs) are built in ref. [7], taking advantage of the module combination configuration synthesis approach and ankle motion fitting models, in order to accommodate the complexity of the ankle joint. The ankle models can be fitted correctly by the GSPMs. Despite the fact that the suggested GSPMs offer a useful starting point for the development and study of ankle rehabilitation robots, the complexity of these GSPMs poses significant control issues. But as Fig. 2 illustrates, the mechanism in our study considers that the AJC’s center of rotation can be permanently located [8–10] at the talocrural joint, which is the joint between the tibia and the talus bone in the ankle [11]. This assumption reduces the mechanism’s complexity and this can result in a reduction in control difficulty. So, we are able to implement more advanced control techniques.

Figure 1. Morphological mobilities of ankle-foot complex.

Figure 2. Rehabilitation device and its adjustment with the ankle joint complex [4].

1.2. Mathematical modeling challenges and control strategies

The presented mechanism with three degrees of freedom is shown in Fig. 2, which offers the ankle three rotational degrees of freedom and has an appropriate adaptation to the kinematics of the human ankle [8]. Although 3-RRR spherical parallel manipulator (RSPM) is an overconstrained mechanism, non-overconstrained versions of the manipulator, such as 3-RCC, 3-RRS, 3-RUU, 2-RSC-RRR, and RSC-RRU-RRC [12], can be used for the evaluation of dynamic analysis and control design methods. This modification leads to consequences such as different position and orientation accuracies. In this study, the 3-RCC assembly version of the manipulator is used for applying the proposed methods because of its highest accuracy [13] in terms of orientation of MP, and we assume that the displacements caused by manufacturing errors are minor displacements. So far, most studies have focused on the kinematic issues of RSPMs, while quite a few have referred to the dynamic control of these robots. For RSPMs, several studies have been reported on a variety of relevant problems, which include modeling of workspace in joint and Cartesian space [14–17], inverse and forward kinematic analysis to obtain analytically unique real-time solutions [18–25], design and optimization [26–32], design and robustness [33], singularity analysis, and derivation of Jacobian matrices [34, 35]. However, the constrained kinematic analysis has gone unnoticed in the literature. As will be shown in this paper, the constrained kinematic analysis can be useful for the extraction of the robot’s closed-form model. Furthermore, the explicit dynamic model of RSPM is needed for model-based controller design in various control strategies for rehabilitation exercises [36]. The complexity of RSPM, with its closed-loop structure and kinematic constraints, compels researchers to employ control techniques that are independent of the exact dynamic model of the robot. However, the problem of dynamic analysis of RSPMs has been addressed in various works using Newton-Euler, the Lagrange method [37], the recursive matrix method [38], the principle of virtual work [39, 40], and the D’Alembert principle approach [41]. A Newton-Euler approach based on screw theory has been proposed to obtain actuator torques and constraint wrenches [42]. In ref. [43], a general solution was proposed for the inverse and forward dynamics of parallel robots in joint space using dynamic models of legs and a moving platform. The dynamic modeling techniques in the literature are not convenient when designing model-based controllers for RSPMs.

The foundation of these controllers is referred to as the computed torque control (CTC) and passivity-based control schemes [44]. In the presence of model uncertainties, the CTC approach is unable to provide high trajectory tracking performance [45]. To solve this issue, various adaptive and robust control methods have been proposed. A novel robust adaptive approach when sufficient extraction of the kinematic and dynamic model of the robot is unavailable has been proposed in ref. [46] and is implemented on an RSPM. This technique provides bounded control input and output signals. Despite supplying excessive stabilization in the presence of uncertainties and external disturbances, trajectory tracking functionality has no longer been provided through this method. Li et al. [47] have proposed a spherical parallel mechanism, which is a special case of RSPMs, with additional sliding actuators and four moving links. Trajectory tracking control of this manipulator has been performed using a robust adaptive switching learning algorithm. A promising structure of the adaptive control approach, which requires many processes, resorts to a linear regression model of the robot through the use of a linearly parameterized dynamic model of a manipulator [48]. Moreover, time delay controller [49] and time delay estimation (TDE) [50] techniques can be employed for estimating the unmodeled dynamics and uncertainties through the usage of all state variables and their first and second delayed derivatives. As these approaches call for information on the acceleration of state variables, they might cause non-negligible noise in the estimation process [41].

As sufficient system identification for the robot is available, it is more practical and feasible to use model-based controllers and fuse them with an adaptation law. We have used the backstepping approach for the known dynamic model, and in order to compensate for model uncertainties, an adaptation law has been combined with backstepping method using Lyapunov stability analysis.

Regarding the control strategy, it is worth mentioning that controllers can be implemented in joint [51] or Cartesian [48] space. Cartesian space control is more efficient than joint space control in parallel robots due to its parallel kinematic chains and complex constraints. Because the position of MP is determined by the state of the robot and the ankle orientation, Cartesian space control theoretically provides state feedback control, resulting in better accuracy than joint space control, which is no longer state feedback control [52]. From the literature, many proposed control approaches for PRs are implemented in joint space [53], and despite their adequate performance, they can be ameliorated by implementing them in Cartesian space.

1.3. Hysteresis loops of the ankle-resistance torques

Since RSPM is used for ankle rehabilitation applications, to achieve high position tracking accuracy and training safety during passive rehabilitation exercises for patients with increased ankle stiffness due to a stroke or muscle injury, the ankle-resistance torque (ART) [54] applied to the robot must be considered in the control design process. ART, which is caused by stiffness in the ankle joint, is depicted in Fig. 3. The design of passive controllers for ankle rehabilitation robots requires enough knowledge of the ART dynamic model as well as the robot. It should be considered that ankle stiffness varies across its range of motion during robotic rehabilitation training and depends on individual characteristics such as weight, height, gender, age, and type of injury. Therefore, it is important that the robot provides assistive torque with adjustable compliance. An estimation of passive ankle stiffness is useful to compensate for ankle resistance in real time. Passive ankle stiffness has been extensively measured and reported in both healthy and neurologically impaired subjects [54, 55]. In addition, some mechanical computational ankle models [56] have been developed to assess ART, including muscles with the feature of subject-specific adjustments [56, 57]. Despite their ability to simulate spasticity, computational models may not be suitable for control purposes because of their complexity [58].

Figure 3. Applied ankle resistance torque $\mathcal{F}_{r}$ to the robot in Cartesian space.

Although any unknown dynamics from environmental sources could be considered external disturbances, if these disturbances are very large, the system may become unstable or its performance may decline significantly. In this circumstance, one way is to take steps to identify the source of the disturbances and reduce their effects by adding additional components to the controller. However, identifying the exact model for ART is not available, and using a high-gain proportional-derivative (PD) controller may be proposed in order to compensate for the error between the precise and estimated values of disturbances. As shown in ref. [36], various passive motion control schemes can also be used for stiffness control of the robot, given the controller parameters. Thus, using only a high-gain controller increases the stiffness of the robot and the Cartesian wrench applied to the patient’s ankle joint and amplifies the oscillations and vibrations. These phenoniums can affect the comfort and safety of the patient during the rehabilitation process.

1.4. Article’s motivation and organization

The main contribution of this paper is the development of a novel computational framework for an explicit dynamic solution of the 3-RCC spherical parallel mechanism and the design of an adaptive robust controller with trajectory tracking capability based on the established dynamic model for ankle rehabilitation applications. Motivated by the abovementioned points, the following is the paper’s structure.

In Sect. 2, a complete analytical inverse kinematic model of the robot is derived using the closed-loop vector method. Then the Jacobian matrices of the robot are extracted. In Sect. 3, based on the screw theory, the extraction of the explicit dynamic model of the RSPM in Cartesian space is presented. In Section 4, a self-tuning backstepping controller (STBC) is designed in task space by providing the Lyapunov stability analysis. Time-varying uncertainties have been taken into consideration within the design process. In Sect. 5, a predictive model based on the established robot’s dynamic and a linear predictive model of ankle stiffness is proposed to compensate for the ART’s dynamic. In Sect. 6, the 3-RCC version of the manipulator is assembled in a multibody environment, and co-simulations are conducted in different scenarios on a virtual prototyping simulation to confirm the proposed explicit dynamic solution and controller. Prior to creating a prototype, the co-simulations can provide a computer-based benchmark for validating the control algorithm and improving the closed-loop system. As a result, it offers a more intuitive simulation than the majority of current ones, which rely solely on proposed algorithms.

Figure 4. The overconstrained 3-RRR spherical parallel manipulator.

2. Architecture and kinematics description

The computer-aided design (CAD) model of the robot is shown in Fig. 4. The moving platform and the base platform (BP) of the robot are connected to each other by three equally spaced parallel chains. The chains are denoted as being counterclockwise (e.g., i = 1, 2, 3). Each chain consists of two rigid, curved links. The robot has three revolute joints in each limb. The joint connecting the first and second links of each limb is named the passive joint. The first revolute joint of each limb is actuated by a fixed motor. Unit vectors defined along the axes of the active and passive joints are denoted by $\mathbf{u}_{i}$ , $\mathbf{w}_{i}$ , respectively. Moreover, the unit vectors fixed to the MP are denoted by $\mathbf{v}_{i}$ , $i=1,2,3$ . These vectors also intersect at a point called the center of rotation, and the MP has an arbitrary rotation with respect to that point. The arc angles of the proximal and distal links are indicated by $\alpha _{1}$ and $\alpha _{2}$ , respectively. Moreover, $\gamma$ is the internal angle between the axis perpendicular to BP and the axis of the first revolute joint of the robot. The right-hand reference coordinate system is chosen such that the origin is at the robot’s center of rotation. The mentioned kinematic parameters of the robot are shown in Fig. 5. Now, following the earlier definition of the reference coordinate system, the unit vectors $\mathbf{u}_{\mathbf{i}}$ , $i=1,2,3$ are defined as

(1) \begin{equation}\mathbf{u}_{\mathbf{i}}=\left[\begin{array}{c@{\quad}c@{\quad}c} -{sin } \left(\eta _{i}\right){sin }\ \gamma, & -{cos } \left(\eta _{i}\right){sin }\ \gamma, & {cos }\ \gamma, \end{array}\right]^{T} \eta _{i}=\frac{2(i-1)}{3}\pi\end{equation}

The unit vectors of the intermediate joints are obtained in terms of the input joint angles described in the fixed frame as follows:

(2) \begin{equation}\mathbf{w}_{i}=\left[\begin{array}{c} -\textrm{s}\alpha _{1}\left(c\eta _{i}c\Theta _{i}-c\gamma s\eta _{i}s\Theta _{i}\right)-c\alpha _{1}s\eta _{i}s\gamma \\[3pt] \textrm{s}\alpha _{1}\left(c\Theta _{i}s\eta _{i}+c\eta _{i}c\gamma s\Theta _{i}\right)-c\alpha _{1}c\eta _{i}s\gamma \\[3pt] c\alpha _{1}c\gamma +s\alpha _{1}s\gamma s\Theta _{i} \end{array}\right]\end{equation}

where c indicates cos and s indicates sin, and $\Theta _{i}, i=1,2,3$ define the angular states of the respective active joints. Moreover, $\mathbf{v}_{\mathbf{i}}$ only depends on the orientation of MP and is denoted as

(3) \begin{equation}\mathbf{v}_{\mathbf{i}}=\mathbf{Q}\mathbf{v}_{\mathbf{i}\mathbf{0}}\end{equation}

where the rotation matrix Q describes the orientation of the MP with respect to the reference coordinate system. Matrix Q can be chosen using the XYZ-Euler angle conventions [59], and the reference configuration of $\mathbf{v}_{\mathbf{i}}$ is chosen such that:

\begin{equation*} \mathbf{v}_{\mathbf{10}}=\mathbf{u}_{\mathbf{3}}\quad \mathbf{v}_{\mathbf{20}}=\mathbf{u}_{\mathbf{1}}\quad \mathbf{v}_{\mathbf{30}}=\mathbf{u}_{\mathbf{2}} \end{equation*}

3. Inverse kinematic

3.1. Analytical solution for active joints

For the closed-loop chain of the parallel manipulator, the following constraint equations hold as the loop closure:

(4) \begin{equation}\mathbf{w}_{\mathbf{i}}.\mathbf{v}_{\mathbf{i}}=\cos \alpha _{2}\end{equation}

Substituting Eqs. (2) and (3) into Eq. (4), inverse kinematic solutions can be solved analytically by the following uncoupled nonlinear algebraic equations for active joint angles $\Theta _{i}$ :

(5) \begin{equation}A_{i}{T_{i}}^{2}+B_{i}T_{i}+T_{i}=0\quad i=1,2,3\end{equation}

with

(6) \begin{equation}T_{i}={\tan } \left(\frac{\Theta _{i}}{2}\right)\end{equation}

In order to obtain unique inverse kinematic solutions, the positive roots of Eq. (6) are considered due to the accepted initial robot configuration where all three proximal links are rotated in the positive direction of the actuated joints.

Figure 5. Limbs, vector of active and passive joints, links, kinematic parameters, and reference coordinate system of RSPM.

Figure 6. Relative degree between proximal and distal links.

3.2. Analytical solution for passive joints

Let $\unicode{x1D6C3}=[\begin{array}{c@{\quad}c@{\quad}c} \beta _{1}, & \beta _{2}, & \beta _{3} \end{array}]^{T}$ denote the vector of the intermediate joints, and $\varphi$ , $\theta$ , and $\psi$ denote the MP motion variables. To calculate the passive joint angles $\beta _{i}, i=1,2,3$ , and according to Fig. 6, two coordinate systems are used, the z-axes of which are directed along the intermediate joint axes, but the y-axis corresponds to the first coordinate system $(\mathbf{y}_{1,i})$ is perpendicular to the plane in which the proximal link is located, and the y-axis of the second coordinate system $(\mathbf{y}_{2,i})$ is perpendicular to the plane in which the distal link is located, i.e.,

\begin{equation*} \mathbf{y}_{2,i}\left(\varphi,\theta,\psi,\Theta _{i}\right)=\mathbf{v}_{\mathbf{i}}\times \mathbf{w}_{\mathbf{i}} \end{equation*}

\begin{equation*} \mathbf{y}_{1,i}\left(\Theta _{i}\right)=\mathbf{w}_{\mathbf{i}}\times \mathbf{u}_{\mathbf{i}} \end{equation*}

The dot product of these vectors determines the cosine of the relative angle that proximal and distal links have to each other. Thus, the closed-loop vector equation for passive joints can be written as follows:

(7) \begin{equation}\cos \beta _{i}\left(\varphi,\theta,\psi,\Theta _{i}\right)=\frac{\mathbf{y}_{1,i}.\mathbf{y}_{2,i}}{\left\| \mathbf{y}_{1,i}\right\| \left\| \mathbf{y}_{2,i}\right\| }\end{equation}

Eq. (7) is proposed to represent the constrained kinematics that hold for passive joints. By using this equation, the inverse kinematic problem of the robot in terms of the passive joint angles can be analytically obtained.

4. Jacobian matrices of the manipulator

4.1. Jacobian matrix in terms of the active joint rates

Let $\unicode{x1D6DA}$ denote the angular velocity vector of the MP, $\unicode{x1D6C9}=[\begin{array}{c@{\quad}c@{\quad}c} \Theta _{1}, & \Theta _{2}, & \Theta _{3} \end{array}]^{T}$ denote the vector of the actuated joint angles, and $\mathbf{x}=[\begin{array}{c@{\quad}c@{\quad}c} \varphi, & \theta, & \psi \end{array}]^{T}$ denote the vector of MP variables. The relationship between the vector of active joint rate $\dot{\unicode{x1D6C9}}$ and the angular velocity of MP can be written as [33]:

(8) \begin{equation}\dot{\unicode{x1D6C9}}=\mathbf{J}_{a}\unicode{x1D6DA}\end{equation}

where the Jacobian matrix $\mathbf{J}_{a}$ can be expressed as:

(9) \begin{equation}\mathbf{J}_{a}=-{\mathbf{J}_{\theta }}^{-1}\mathbf{J}_{x}\end{equation}

in which,

\begin{equation*} \mathbf{J}_{\theta }=diag(\mathbf{w}_{1}\times \mathbf{u}_{1}.\mathbf{v}_{1}, \mathbf{w}_{2}\times \mathbf{u}_{2}.\mathbf{v}_{2},\mathbf{w}_{3}\times \mathbf{u}_{3}.\mathbf{v}_{3}), \end{equation*}

and

\begin{equation*} \mathbf{J}_{x}=\left[\begin{array}{c} \left({\mathbf{w}_{1}}\times {\mathbf{v}_{1}}\right)^{T}\\ \left({\mathbf{w}_{2}}\times {\mathbf{v}_{2}}\right)^{T}\\ \left({\mathbf{w}_{3}}\times {\mathbf{v}_{3}}\right)^{T} \end{array}\right] \end{equation*}

For fully parallel manipulators, $\mathbf{J}_{\theta }\ \textrm{and}\ \mathbf{J}_{x}$ are both $n\times n$ matrices [36] ( $n=3$ for RSPM). It should be noted that in order to compute the angular velocity of MP, a rotation matrix E, which is a function of Euler angles, should be used as follows [59]:

(10) \begin{equation}\unicode{x1D6DA}=\mathbf{E}\dot{\mathbf{x}}\end{equation}

Thus, Eq. (8) can be rewritten as:

(11) \begin{equation}\dot{\unicode{x1D6C9}}=\mathbf{J}_{\sigma }\dot{\mathbf{x}}\end{equation}

in which the $3\times 3$ matrix $\mathbf{J}_{\sigma }$ is defined as:

(12a) \begin{equation}\mathbf{J}_{\sigma }=\mathbf{J}_{a}\mathbf{E}=\left[\begin{array}{c@{\quad}c@{\quad}c} {J_{\sigma }}_{11} & {J_{\sigma }}_{12} & {J_{\sigma }}_{13}\\ {J_{\sigma }}_{21} & {J_{\sigma }}_{22} & {J_{\sigma }}_{23}\\ {J_{\sigma }}_{31} & {J_{\sigma }}_{32} & {J_{\sigma }}_{33} \end{array}\right]=\left[\begin{array}{c} \unicode{x1D6D4}_{1}\\ \unicode{x1D6D4}_{2}\\ \unicode{x1D6D4}_{3} \end{array}\right]_{3\times 3}\end{equation}

As a result, $\unicode{x1D6D4}_{i}$ is a $1\times 3$ vector that is defined as:

(12b) \begin{align}\unicode{x1D6D4}_{i}=\left(\begin{array}{c} {J_{\sigma }}_{i1}\\[3pt] {J_{\sigma }}_{i2}\\[3pt] {J_{\sigma }}_{i3} \end{array}\right)^{T}, i=1,2,3\end{align}

4.2. Constrained kinematics

In order to derive a relation that maps the MP motion variable rates to the vector of passive joint rates $\dot{\unicode{x1D6C3}}$ , each passive joint angle is considered to be a function of MP Euler angles and active joint variables, namely,

\begin{equation*} \beta _{i}=\beta _{i}\left(\Theta _{i}, \varphi, \psi, \theta \right), i=1,2,3 \end{equation*}

Therefore, the time derivative of $\beta _{i}$ can be written as follows:

(13) \begin{equation}\dot{\beta }_{i}=\nabla {\beta _{i}}\left(\Theta _{i}, \theta, \psi, \varphi \right)_{1\times 4}\left[\begin{array}{c} \dot{\Theta }_{i}\\[3pt] \dot{\theta }\\[3pt] \dot{\psi }\\[3pt] \dot{\varphi } \end{array}\right]_{4\times 1}=\frac{\partial \beta _{i}}{\partial \Theta _{i}}\dot{\Theta }_{i}+\frac{\partial \beta _{i}}{\partial \theta }\dot{\theta }+\frac{\partial \beta _{i}}{\partial \psi }\dot{\psi }+\frac{\partial \beta _{i}}{\partial \varphi }\dot{\varphi }, i=1,2,3\end{equation}

Writing Eq. 13 three times for each limb and substituting the value of $\dot{\Theta }_{i}$ from Eq. (11) results in the $3\times 3$ matrix that relates the passive joint rates and the rates of Euler angles.

(14) \begin{equation}\dot{\unicode{x1D6C3}}=\mathbf{J}_{\mathbf{p}}\dot{\mathbf{x}}\end{equation}

where

\begin{equation*} \mathbf{J}_{\mathbf{p}}=\left[\begin{array}{c@{\quad}c@{\quad}c} \dfrac{\partial \beta _{1}}{\partial \Theta _{1}}{J_{\sigma }}_{11}+\dfrac{\partial \beta _{1}}{\partial \varphi } & \dfrac{\partial \beta _{1}}{\partial \Theta _{1}}{J_{\sigma }}_{12}+\dfrac{\partial \beta _{1}}{\partial \theta } & \dfrac{\partial \beta _{1}}{\partial \Theta _{1}}{J_{\sigma }}_{13}+\dfrac{\partial \beta _{1}}{\partial \psi }\\[15pt] \dfrac{\partial \beta _{2}}{\partial \Theta _{2}}{J_{\sigma }}_{21}+\dfrac{\partial \beta _{2}}{\partial \varphi } & \dfrac{\partial \beta _{2}}{\partial \Theta _{2}}{J_{\sigma }}_{22}+\dfrac{\partial \beta _{2}}{\partial \theta } & \dfrac{\partial \beta _{2}}{\partial \Theta _{2}}{J_{\sigma }}_{23}+\dfrac{\partial \beta _{2}}{\partial \psi }\\[15pt] \dfrac{\partial \beta _{3}}{\partial \Theta _{3}}{J_{\sigma }}_{31}+\dfrac{\partial \beta _{3}}{\partial \varphi } & \dfrac{\partial \beta _{3}}{\partial \Theta _{3}}{J_{\sigma }}_{32}+\dfrac{\partial \beta _{3}}{\partial \theta } & \dfrac{\partial \beta _{3}}{\partial \Theta _{3}}{J_{\sigma }}_{33}+\dfrac{\partial \beta _{3}}{\partial \psi } \end{array}\right] \end{equation*}

(15a) \begin{equation}=\left[\begin{array}{c@{\quad}c@{\quad}c} {J_{\rho }}_{11} & {J_{\rho }}_{12} & {J_{\rho }}_{13}\\[4pt] {J_{\rho }}_{21} & {J_{\rho }}_{22} & {J_{\rho }}_{23}\\[4pt] {J_{\rho }}_{31} & {J_{\rho }}_{32} & {J_{\rho }}_{33} \end{array}\right]=\left[\begin{array}{c} \unicode{x1D6D2}_{\mathbf{1}}\\[4pt] \unicode{x1D6D2}_{\mathbf{2}}\\[4pt] \unicode{x1D6D2}_{\mathbf{3}} \end{array}\right]_{3\times 3}\end{equation}

So, $1\times 3$ vector $\unicode{x1D6D2}_{i}$ could be written as below:

(15b) \begin{equation}\unicode{x1D6D2}_{i}=\left(\begin{array}{c} \dfrac{\partial \beta _{i}}{\partial \Theta _{i}}{J_{\sigma }}_{i1}+\dfrac{\partial \beta _{i}}{\partial \varphi }\\[15pt] \dfrac{\partial \beta _{i}}{\partial \Theta _{i}}{J_{\sigma }}_{i2}+\dfrac{\partial \beta _{i}}{\partial \theta }\\[15pt] \dfrac{\partial \beta _{i}}{\partial \Theta _{i}}{J_{\sigma }}_{i3}+\dfrac{\partial \beta _{i}}{\partial \psi } \end{array}\right)^{T}, i=1,2,3\end{equation}

5. Explicit dynamics in task space

In this study, the dynamic modeling of the manipulator is derived based on the Newton-Euler formulation and screw notation [60]. According to this approach, the dynamic model of all moving links of the robot can be written as follows:

(16) \begin{equation}\mathbf{M}\dot{\mathbf{t}}+\mathbf{WMt}-{}^{G}{\mathbf{w}}{}={}^{A}{\mathbf{w}}{}+{}^{C}{\mathbf{w}}{}\end{equation}

where

(17) \begin{equation}\mathbf{t}=[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} \mathbf{t}_{\mathbf{1}} & \mathbf{t}_{\mathbf{2}} & \ldots & \mathbf{t}_{\mathbf{7}} \end{array}]^{T},\quad {}^{l}{\mathbf{w}}{}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {}^{l}{\mathbf{w}_{1}}{} & {}^{l}{\mathbf{w}_{2}}{} & \ldots & {}^{l}{\mathbf{w}_{7}}{} \end{array}\right]^{T}\end{equation}

and

(18) \begin{equation}\mathbf{M}=diag\left(\mathbf{M}_{1},\ldots,\mathbf{M}_{7}\right)\end{equation}

(19) \begin{equation}\mathbf{W}=diag\left(\mathbf{W}_{1},\ldots,\mathbf{W}_{7}\right)\end{equation}

in which $\mathbf{t}_{j}$ denotes the twist array of the jth moving link. Moreover, ${}^{l}{\mathbf{w}_{j}}{}, l=A,G,C$ represents the working, gravitational, and non-working constraint wrench exerted on body j, where it is assumed that torques are applied at the center of mass of the body [61], and $\mathbf{M}_{j}$ is the 6 × 6 inertia matrix of body j. Also, $\mathbf{W}_{j}$ is the 6 × 6 angular velocity matrix of the same body. The foregoing matrices are defined below for the jth body,

\begin{equation*}\mathbf{t}_{j}=\left[\begin{array}{c} \unicode{x1D6DA}_{j}\\[5pt] \mathbf{v}_{j} \end{array}\right] \quad \mathbf{M}_{j}=\left[\begin{array}{c@{\quad}c} \mathbf{I}_{j} & \mathbf{0}_{3\times 3}\\[5pt] \mathbf{0}_{3\times 3} & m_{j}\mathbf{1}_{3\times 3} \end{array}\right]\quad \mathbf{W}_{j}=\left[\begin{array}{c@{\quad}c} \boldsymbol{\Omega }_{j} & \mathbf{0}_{3\times 3}\\[5pt] \mathbf{0}_{3\times 3} & \mathbf{0}_{3\times 3} \end{array}\right],\quad j=1,2,\ldots,7\end{equation*}

where $\unicode{x1D6DA}_{j}$ indicates the angular velocity and $\mathbf{v}_{j}$ denotes the linear velocity of the center of mass of each link. Also, $\mathbf{I}_{j}$ and $\boldsymbol{\Omega }_{j}$ denote the inertia matrix and the cross-product matrix of $\unicode{x1D6DA}_{j}$ , respectively [61]. Moreover, $m_{j}$ denotes the mass of link j. The robot actuator-wrench array ${}^{A}{\mathbf{w}}{}$ can be defined as a function of actuator torques by a transformation, as shown below [61]:

(20) \begin{equation}{}^{A}{\mathbf{w}}{}={}^{\mathbf{A}}{\mathbf{T}}{}\unicode{x1D6D5}\end{equation}

where ${}^{\mathbf{A}}{\mathbf{T}}{}$ is the 6n × k actuator-wrench shaping matrix of the robot and τ denotes a k-dimensional array of actuator torques [62] ( $n=7, k=3$ for RSPM and $\unicode{x1D6D5}=[\begin{array}{c@{\quad}c@{\quad}c} \tau _{1} & \tau _{2} & \tau _{3} \end{array}]^{T}$ ). Since no actuator wrench is applied to bodies number 2, 3, 5, or 7, one has

(21) \begin{equation}{}^{A}{\mathbf{w}_{2}}{}={}^{A}{\mathbf{w}_{3}}{}={}^{A}{\mathbf{w}_{5}}{}={}^{A}{\mathbf{w}_{7}}{}=\mathbf{0}_{6\times 1}\end{equation}

As each actuator applies a torque on the proximal links about the axis of rotation of the first joint in each chain of the robot, the actuator-wrench shaping matrix of links 1, 4, and 6 can be defined as:

(22) \begin{equation}{}^{A}{\mathbf{w}_{1}}{}=\left[\begin{array}{c} \mathbf{u}_{1}\tau _{1}\\[4pt] \mathbf{0}_{3\times 1} \end{array}\right],\quad {}^{A}{\mathbf{w}_{4}}{}=\left[\begin{array}{c} \mathbf{u}_{2}\tau _{2}\\[4pt] \mathbf{0}_{3\times 1} \end{array}\right],\quad {}^{A}{\mathbf{w}_{6}}{}=\left[\begin{array}{c} \mathbf{u}_{3}\tau _{3}\\[4pt] \mathbf{0}_{3\times 1} \end{array}\right]\end{equation}

Hence, by resorting to Eq. (21) and (22), ${}^{A}{\mathbf{T}}{}$ is obtained as follows:

(23) \begin{equation}{}^{\mathbf{A}}{\mathbf{T}}{}=\left[\begin{array}{c@{\quad}c@{\quad}c} \mathbf{u}_{\mathbf{1}} & \mathbf{0}_{\mathbf{3}\times \mathbf{1}} & \mathbf{0}_{\mathbf{3}\times \mathbf{1}}\\[3pt] \mathbf{0}_{\mathbf{15}\times \mathbf{1}} & \mathbf{0}_{\mathbf{15}\times \mathbf{1}} & \mathbf{0}_{\mathbf{15}\times \mathbf{1}}\\[3pt] \mathbf{0}_{\mathbf{3}\times \mathbf{1}} & \mathbf{u}_{\mathbf{2}} & \mathbf{0}_{\mathbf{3}\times \mathbf{1}}\\[3pt] \mathbf{0}_{\mathbf{9}\times \mathbf{1}} & \mathbf{0}_{\mathbf{9}\times \mathbf{1}} & \mathbf{0}_{\mathbf{9}\times \mathbf{1}}\\[3pt] \mathbf{0}_{\mathbf{3}\times \mathbf{1}} & \mathbf{0}_{\mathbf{3}\times \mathbf{1}} & \mathbf{u}_{\mathbf{3}}\\[3pt] \mathbf{0}_{\mathbf{9}\times \mathbf{1}} & \mathbf{0}_{\mathbf{9}\times \mathbf{1}} & \mathbf{0}_{\mathbf{9}\times \mathbf{1}} \end{array}\right]\end{equation}

where ${}^{\textrm{A}}{\mathbf{T}}{}$ is a $42\times 3$ matrix that transfers the actuator torques to the center of mass of each link. Hence, by using Eqs. (16) and (20), the equations of motion of the parallel manipulator can be rewritten as follows:

(24) \begin{equation}\mathbf{M}\dot{\mathbf{t}}+\mathbf{WMt}-{}^{\mathbf{G}}{\mathbf{w}}{}={}^{\mathbf{C}}{\mathbf{w}}{}+{}^{\mathbf{A}}{\mathbf{T}}{}\unicode{x1D6D5}\end{equation}

In the following, the derivation of dynamic equations is considered in explicit form. For this purpose, the robot’s twist array can be written as a function of the rates of generalized coordinates through a Jacobian matrix [61], namely,

(25) \begin{equation}\mathbf{t}=\mathbf{J}\dot{\mathbf{x}}\end{equation}

(26) \begin{equation}\dot{\mathbf{t}}=\dot{\mathbf{J}}\dot{\mathbf{x}}+\mathbf{J}\ddot{\mathbf{x}}\end{equation}

$\mathbf{J}$ is a Jacobian matrix that maps the robot’s independent variables to the twist array of all the robot’s moving links. In order to derive this Jacobian, the twist array of each moving link should be written as a function of the main independent coordinate rates. Using Eqs. (10), (11), (12), (14), and (15) and the vectors shown in Fig. 7, the angular velocity vector of moving links can be written as follows:

(27a) \begin{equation}\unicode{x1D6DA}_{1}=\mathbf{u}_{1}\unicode{x1D6D4}_{1}\dot{\mathbf{x}}\end{equation}

(27b) \begin{equation}\unicode{x1D6DA}_{2}=\mathbf{u}_{1}\unicode{x1D6D4}_{1}\dot{\mathbf{x}}+\mathbf{w}_{1}\unicode{x1D6D2}_{1}\dot{\mathbf{x}}\end{equation}

(27c) \begin{equation}\unicode{x1D6DA}_{3}=\mathbf{E}\dot{\mathbf{x}}\end{equation}

(27d) \begin{equation}\unicode{x1D6DA}_{4}=\mathbf{u}_{2}\unicode{x1D6D4}_{2}\dot{\mathbf{x}}\end{equation}

(27e) \begin{equation}\unicode{x1D6DA}_{5}=\mathbf{u}_{2}\unicode{x1D6D4}_{2}\dot{\mathbf{x}}+\mathbf{w}_{2}\unicode{x1D6D2}_{2}\dot{\mathbf{x}}\end{equation}

(27f) \begin{equation}\unicode{x1D6DA}_{6}=\mathbf{u}_{3}\unicode{x1D6D4}_{3}\dot{\mathbf{x}}\end{equation}

(27g) \begin{equation}\unicode{x1D6DA}_{7}=\mathbf{u}_{3}\unicode{x1D6D4}_{3}\dot{\mathbf{x}}+\mathbf{w}_{3}\unicode{x1D6D2}_{3}\dot{\mathbf{x}}\end{equation}

Figure 7. The vectors connected to proximal and distal links and moving platform.

Then, using some algebraic manipulations, the linear velocity vector of each moving link can be written in a suitable form as follows:

(28a) \begin{equation}\mathbf{v}_{1}=-CPM\left(\mathbf{c}_{1}\right).(\mathbf{u}_{1}\unicode{x1D6D4}_{1}).\dot{\mathbf{x}}\end{equation}

(28b) \begin{equation}\mathbf{v}_{2}=-CPM\left(\mathbf{p}_{1}\right).(\mathbf{u}_{1}\unicode{x1D6D4}_{1}).\dot{\mathbf{x}}-CPM\left(\unicode{x1D6CC}_{1}\right).\left(\mathbf{u}_{1}\unicode{x1D6D4}_{1}+\mathbf{w}_{1}\unicode{x1D6D2}_{1}\right).\dot{\mathbf{x}}\end{equation}

(28c) \begin{equation}\mathbf{v}_{3}=-CPM\left(\unicode{x1D6DD}\right).\mathbf{E}\dot{\mathbf{x}}\end{equation}

(28d) \begin{equation}\mathbf{v}_{4}=-CPM\left(\mathbf{c}_{2}\right).(\mathbf{u}_{2}\unicode{x1D6D4}_{2}).\dot{\mathbf{x}}\end{equation}

(28e) \begin{equation}\mathbf{v}_{5}=-CPM\left(\mathbf{p}_{2}\right).\left(\mathbf{u}_{2}\unicode{x1D6D4}_{2}\right).\dot{\mathbf{x}}-CPM\left(\unicode{x1D6CC}_{2}\right).\left(\mathbf{u}_{2}\unicode{x1D6D4}_{2}+\mathbf{w}_{2}\unicode{x1D6D2}_{2}\right).\dot{\mathbf{x}}\end{equation}

(28f) \begin{equation}\mathbf{v}_{6}=-CPM(\mathbf{c}_{3}).(\mathbf{u}_{3}\unicode{x1D6D4}_{3}).\dot{\mathbf{x}}\end{equation}

(28g) \begin{equation}\mathbf{v}_{7}=-CPM\left(\mathbf{p}_{3}\right).\left(\mathbf{u}_{3}\unicode{x1D6D4}_{3}\right).\dot{\mathbf{x}}-CPM\left(\unicode{x1D6CC}_{3}\right).\left(\mathbf{u}_{3}\unicode{x1D6D4}_{3}+\mathbf{w}_{3}\unicode{x1D6D2}_{3}\right).\dot{\mathbf{x}}\end{equation}

in which $CPM()$ indicates the cross-product matrix [61] of a vector. Thus, the Jacobian matrix $\mathbf{J}$ can be driven as follows:

(29) \begin{equation}\boldsymbol{J}=\left[\begin{array}{c} \boldsymbol{u}_{\mathbf{1}}\unicode{x1D6D4}_{\mathbf{1}}\\[2pt] -CPM(\boldsymbol{c}_{1}).\boldsymbol{u}_{1}\unicode{x1D6D4}_{1}\\[2pt] \boldsymbol{u}_{1}\unicode{x1D6D4}_{1}+\boldsymbol{w}_{1}\unicode{x1D6D2}_{1}\\[2pt] -CPM\left(\boldsymbol{p}_{1}\right).(\boldsymbol{u}_{1}\unicode{x1D6D4}_{1})-CPM\left(\unicode{x1D6CC}_{1}\right).\left(\boldsymbol{u}_{1}\unicode{x1D6D4}_{1}+\boldsymbol{w}_{1}\unicode{x1D6D2}_{1}\right)\\[2pt] \boldsymbol{E}\\[2pt] -CPM(\unicode{x1D6DD})\boldsymbol{E}\\[2pt] \boldsymbol{u}_{2}\unicode{x1D6D4}_{2}\\[2pt] -CPM\left(\boldsymbol{c}_{2}\right).\boldsymbol{u}_{2}\unicode{x1D6D4}_{2}\\[2pt] \boldsymbol{u}_{2}\unicode{x1D6D4}_{2}+\boldsymbol{w}_{2}\unicode{x1D6D2}_{2}\\[2pt] -CPM\left(\boldsymbol{p}_{2}\right).\left(\boldsymbol{u}_{2}\unicode{x1D6D4}_{2}\right)-CPM\left(\unicode{x1D6CC}_{2}\right).\left(\boldsymbol{u}_{2}\unicode{x1D6D4}_{2}+\boldsymbol{w}_{2}\unicode{x1D6D2}_{2}\right)\\[2pt] \boldsymbol{u}_{3}\unicode{x1D6D4}_{3}\\[2pt] -CPM(\boldsymbol{c}_{3}).\boldsymbol{u}_{3}\unicode{x1D6D4}_{3}\\[2pt] \boldsymbol{u}_{3}\unicode{x1D6D4}_{3}+\boldsymbol{w}_{3}\unicode{x1D6D2}_{3}\\[2pt] -CPM\left(\boldsymbol{p}_{3}\right).\left(\boldsymbol{u}_{3}\unicode{x1D6D4}_{3}\right)-CPM\left(\unicode{x1D6CC}_{3}\right).\left(\boldsymbol{u}_{3}\unicode{x1D6D4}_{3}+\boldsymbol{w}_{3}\unicode{x1D6D2}_{3}\right) \end{array}\right]_{42\times 3}\end{equation}

$\mathbf{u}_{\textrm{i}}$ and $\mathbf{w}_{i}$ are unit vectors of active and passive joints, respectively. Moreover, $\mathbf{c}_{i}$ , $\mathbf{p}_{i}$ , $\unicode{x1D6CC}_{i}$ , and $\unicode{x1D6DD}$ are vectors shown in Fig. 7 and defined in the reference coordinate system. Moreover, $\unicode{x1D6D4}_{i}$ and $\unicode{x1D6D2}_{i}$ are defined in Eq. (12b) and Eq. (15b) and $i=1,2,3$ signifies the number of each limb. According to linear homogeneous property on the array of manipulator twist, and as the non-working constraint wrench ${}^{\mathbf{C}}{\mathbf{w}}{}$ produces no work on the robot and its only function is to hold the adjacent links together, it can be concluded that the power generated by this wrench on the twist of the system is zero for any feasible motion of the robot [61], i.e.,

(30) \begin{equation}\mathbf{t}^{\mathbf{T}}{}^{\,\mathbf{C}}{\mathbf{w}}{}=\mathbf{0}\end{equation}

By substituting Eq. (25) in Eq. (30), one has the following relation:

\begin{equation*} \dot{\mathbf{x}}^{\mathbf{T}}\mathbf{J}^{\mathbf{T}}{}^{\,\mathbf{C}}{\mathbf{w}}{}=\mathbf{0} \end{equation*}

Since $\dot{\mathbf{x}}$ cannot be zero, the following equation is obtained:

\begin{equation*} \mathbf{J}^{\mathbf{T}}{}^{\,\mathbf{C}}{\mathbf{w}}{}=\mathbf{0} \end{equation*}

By multiplying $\mathbf{J}^{T}$ in Eq. (24), the dynamic robot equation takes the following form:

(31) \begin{equation}\mathbf{J}^{T}\mathbf{M}\dot{\mathbf{t}}+\mathbf{J}^{T}\mathbf{WMt}-\mathbf{J}^{T}{}^{\,\mathbf{G}}{\mathbf{w}}{}=\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{w}}{}+\mathbf{J}^{T}{}^{\mathbf{C}}{\mathbf{w}}{}=\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\unicode{x1D6D5}\end{equation}

Then by multiplying $({\mathbf{J}^{\mathbf{T}}}{}^{\mathbf{A}}{\mathbf{T}}{})^{-\mathbf{1}}$ in Eq. (31), the following equation is derived:

(32) \begin{equation}\unicode{x1D6D5}=({\mathbf{J}^{T}}{}^{\mathbf{A}}{\mathbf{T}}{})^{-1}\mathbf{J}^{T}\mathbf{M}\dot{\mathbf{t}}+({\mathbf{J}^{T}}{}^{\mathbf{A}}{\mathbf{T}}{})^{-1}\mathbf{J}^{T}\mathbf{WMt}-\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}{}^{\,\mathbf{G}}{\mathbf{w}}{}\end{equation}

Now, by replacing Eqs. (25) and (26) in Eq. (32), the closed-form dynamic formulation for the manipulator in the Cartesian space is obtained as follows:

(33) \begin{equation}\unicode{x1D6D5}=\left[({\mathbf{J}^{T}}{}^{\mathbf{A}}{\mathbf{T}}{})^{-1}\mathbf{J}^{T}\mathbf{MJ}\right]\ddot{\mathbf{x}}+[\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}\mathbf{WMJ}+\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}\mathbf{M}\dot{\mathbf{J}}]\dot{\mathbf{x}}-\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}{}^{\,\mathbf{G}}{\mathbf{w}}{}\end{equation}

Here $\mathbf{J}$ is the Jacobian matrix, which was given in Eq. (29). Eq. (33) is very convenient for applying multiple advanced control schemes. For sake of simplicity, Eq. (33) can be written in the following standard form in task space in the absence of external disturbances:

(34) \begin{equation}\mathbf{H}(\mathbf{x})\ddot{\mathbf{x}}+\mathbf{C}(\mathbf{x},\dot{\mathbf{x}})\dot{\mathbf{x}}+\mathbf{G}(\mathbf{x})=\unicode{x1D6D5}\end{equation}

in which

(35a) \begin{equation}\mathbf{H}\left(\mathbf{x}\right)=({\mathbf{J}^{T}}{}^{\mathbf{A}}{\mathbf{T}}{})^{-1}\mathbf{J}^{T}\mathbf{MJ}\end{equation}

(35b) \begin{equation}\mathbf{C}\left(\mathbf{x},\dot{\mathbf{x}}\right)=\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}\mathbf{WMJ}+\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}\mathbf{M}\dot{\mathbf{J}}\end{equation}

(35c) \begin{equation}\mathbf{G}\left(\mathbf{x}\right)=-\left(\mathbf{J}^{T}{}^{\mathbf{A}}{\mathbf{T}}{}\right)^{-1}\mathbf{J}^{T}{}^{\,\mathbf{G}}{\mathbf{w}}{}\end{equation}

In this formulation, $\mathbf{H}(\mathbf{x})$ , $\mathbf{C}(\mathbf{x},\dot{\mathbf{x}})$ , and $\mathbf{G}(\mathbf{x})$ signify the mass matrix, the Coriolis and centrifugal matrices, and the gravity vector, respectively. The algorithm flowchart of the entire method of deriving the manipulator’s explicit dynamic is depicted in Fig. 8.

Figure 8. Proposed computational framework for deriving the explicit dynamic solution of RSPM in task space by using the analytical solutions of inverse kinematic problem.

6. Self-tuning backstepping controller

To achieve a suitable control performance in the presence of dynamic uncertainties, the control input can be computed as follows:

(36) \begin{equation}\unicode{x1D6D5}=\hat{\mathbf{H}}(\mathbf{x})\mathbf{a}_{r}+\hat{\mathbf{C}}(\mathbf{x},\dot{\mathbf{x}})\dot{\mathbf{x}}+\hat{\mathbf{G}}(\mathbf{x})\end{equation}

In which,

(37) \begin{equation}\mathbf{a}_{r}=\ddot{\mathbf{x}}_{d}+\mathbf{K}_{d}\dot{\mathbf{e}}_{x}+\mathbf{K}_{p}\mathbf{e}_{x}+\unicode{x1D6C5}_{a}\end{equation}

where $\hat{\mathbf{H}}(\mathbf{x})$ , $\hat{\mathbf{C}}(\mathbf{x},\dot{\mathbf{x}})$ , and $\hat{\mathbf{G}}(\mathbf{x})$ denote the estimated values of the inertia matrix, the Coriolis and centrifugal matrix, and the gravity vector of the robot, respectively. $\mathbf{K}_{p}$ and $\mathbf{K}_{d}$ denote the PD controller gains. $\mathbf{x}_{d}$ is vector of the desired orientation of the MP. It is assumed that $\mathbf{x}_{d}$ is a twice continuously differentiable signal. Moreover, $\mathbf{e}_{\mathbf{x}}=\mathbf{x}_{d}-\mathbf{x}$ is the tracking error vector, and $\dot{\mathbf{e}}_{x}=\dot{\mathbf{x}}_{d}-\dot{\mathbf{x}}$ indicates the first derivative of the tracking error vector in task space. $\unicode{x1D6C5}_{\mathbf{a}}=[\begin{array}{c@{\quad}c@{\quad}c} {\delta _{a}}_{1} & {\delta _{a}}_{2} & {\delta _{a}}_{3} \end{array}]^{T}$ signifies an additional term, which is designed to compensate for time-varying uncertainties. By defining the error as the difference between the estimated and the accurate value of the corresponding term as $(\,\tilde{}\,).$ Hence,

\begin{equation*} \tilde{\mathbf{H}}=\hat{\mathbf{H}}-\mathbf{H}\quad \tilde{\mathbf{C}}=\hat{\mathbf{C}}-\mathbf{C}\quad \tilde{\mathbf{G}}=\hat{\mathbf{G}}-\mathbf{G} \end{equation*}

Similarly, $(\,\tilde{}\,)$ notation may be applied to the motion variables as:

\begin{equation*} \tilde{\mathbf{x}}=\mathbf{x}-\mathbf{x}_{d}=-\mathbf{e}_{x} \quad\dot{\tilde{\mathbf{x}}}=\dot{\mathbf{x}}-\dot{\mathbf{x}}_{d}=-\dot{\mathbf{e}}_{x} \end{equation*}

Now substituting Eq. (36) in Eq. (34) and adding and subtracting $\mathbf{H}\mathbf{a}_{r}$ , it yields,

(38) \begin{equation}\ddot{\mathbf{x}}=\mathbf{a}_{\mathbf{r}}+\unicode{x1D6C8}\end{equation}

in which

(39) \begin{equation}\unicode{x1D6C8}=\mathbf{H}^{-1}\left(\tilde{\mathbf{H}}\mathbf{a}_{r}+\tilde{\mathbf{C}}\dot{\mathbf{x}}+\tilde{\mathbf{G}}\right)\end{equation}

where $\unicode{x1D6C8}$ signifies model uncertainties. Now rewrite Eq. (37) using the $(\,\tilde{}\,)$ notion, it yields,

(40) \begin{equation}\mathbf{a}_{r}=\ddot{\mathbf{x}}_{d}-\mathbf{K}_{d} \dot{\tilde{\mathbf{x}}}-\mathbf{K}_{p}\tilde{\mathbf{x}}+\unicode{x1D6C5}_{a}\end{equation}

Substituting Eq. (40) into Eq. (38) yields,

(41) \begin{equation}\ddot{\tilde{\mathbf{x}}}=-\mathbf{K}_{d}\dot{\tilde{\mathbf{x}}}-\mathbf{K}_{p}\tilde{\mathbf{x}}+\unicode{x1D6C5}_{a}+\unicode{x1D6C8}\end{equation}

The closed-loop dynamics outlined in Eq. (41) can be written as a state space form through pure mathematical manipulation by changing variables,

\begin{equation*} \mathbf{z}_{1}=\tilde{\mathbf{x}} \end{equation*}

\begin{equation*} \mathbf{z}_{2}=\dot{\tilde{\mathbf{x}}} \end{equation*}

which results in the following first-order differential equations (ODEs),

(42) \begin{equation}\dot{\mathbf{z}}_{1}=\mathbf{z}_{2}\end{equation}

(43) \begin{equation}\dot{\mathbf{z}}_{2}=-\mathbf{K}_{d}\mathbf{z}_{2}-\mathbf{K}_{p}\mathbf{z}_{1}+\unicode{x1D6C5}_{a}+\unicode{x1D6C8}\end{equation}

Consider the system shown in Eqs. (42) and (43). For designing a backstepping controller, the following steps have been followed. In the first step, $\mathbf{z}_{2}$ can be used as control input for the system presented in Eq. (42). So, $\mathbf{z}_{2}$ is defined as a function of $\mathbf{z}_{1}$ as follows:

(44) \begin{equation}\mathbf{z}_{2}=\unicode{x1D6DF}\left(\mathbf{z}_{1}\right)=-\boldsymbol{\mathcal{B}}\mathbf{z}_{1}\end{equation}

As a result, Eq. (42) is rewritten as follows:

(45) \begin{equation}\dot{\mathbf{z}}_{1}=-\boldsymbol{\mathcal{B}}\mathbf{z}_{1}\end{equation}

where $\boldsymbol{\mathcal{B}}$ is a diagonal positive-definite matrix defined by the user,

(46) \begin{equation}\boldsymbol{\mathcal{B}}=\left[\begin{array}{c@{\quad}c@{\quad}c} \mathcal{B}_{1} & 0 & 0\\ 0 & \mathcal{B}_{2} & 0\\ 0 & 0 & \mathcal{B}_{3} \end{array}\right]\end{equation}

By defining a Lyapunov function as $\mathbf{V}_{{z_{1}}}$ , it is obvious that the origin of Eq. (45) is asymptotically stable if $\dot{\mathbf{V}}_{{z_{1}}}\leq \mathbf{0}$ . Now, let’s define a new variable as,

(47) \begin{equation}\overline{\mathbf{z}}=\mathbf{z}_{2}-\unicode{x1D6DF}\left(\mathbf{z}_{1}\right)\end{equation}

Eq. (45) holds if $\overline{\mathbf{z}}$ converges to zero. For ensuring this convergence, Eq. (42) is written as below by using Eq. (47) and some manipulations,

(48) \begin{equation}\dot{\mathbf{z}}_{1}=\overline{\mathbf{z}}-\boldsymbol{\mathcal{B}}\mathbf{z}_{1}\end{equation}

Then by differentiating Eq. (47), the following equation is obtained:

(49) \begin{equation}\dot{\overline{\mathbf{z}}}=\boldsymbol{\mathcal{B}}\overline{\mathbf{z}}-\mathbf{K}_{d}\mathbf{z}_{2}-(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2})\mathbf{z}_{1}+\unicode{x1D6C5}_{a}+\unicode{x1D6C8}\end{equation}

By choosing the controller gains $\mathcal{B}_{1}, \mathcal{B}_{2}, \mathcal{B}_{3}$ such that the matrix $\boldsymbol{\mathcal{B}}$ becomes Hurwitz, a symmetric positive-definite matrix P exists to satisfy the Lyapunov equation for any arbitrary symmetric positive-definite matrix $\boldsymbol{\Upsilon }$ as follows:

(50) \begin{equation}\boldsymbol{\mathcal{B}}^{T}\mathbf{P}+\mathbf{P}\boldsymbol{\mathcal{B}}^{T}=-\boldsymbol{\Upsilon }\end{equation}

On the other hand, consider the error dynamic of the uncertainty vector as follows:

\begin{equation*} \mathbf{E }=\hat{\unicode{x1D6C8}}-\unicode{x1D6C8} \end{equation*}

\begin{equation*} \dot{\mathbf{E }}=\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}} \end{equation*}

in which $\hat{\unicode{x1D6C8}}$ and $\unicode{x1D6C8}$ are the estimated and the unknown actual uncertainty vectors, respectively. Now it is desirable to design a robust backstepping controller with the ability to estimate time-varying uncertainties $\unicode{x1D6C8}$ . For this purpose, by choosing $\mathbf{V}_{{z_{1}}}=\frac{1}{2}{\mathbf{z}_{1}}^{T}\mathbf{z}_{1}$ , the following Lyapunov function is proposed:

(51) \begin{equation}\mathbf{V}_{C}=\frac{1}{2}{\mathbf{z}_{1}}^{T}\mathbf{z}_{1}+\overline{\mathbf{z}}^{T}\mathbf{P}\overline{\mathbf{z}}+\frac{1}{2}\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\mathbf{E }\end{equation}

where $\mathbf{G}_{a}$ is a diagonal positive-definite constant matrix, which is adaptive gain control. The time derivative of the Lyapunov function $\mathbf{V}_{C}$ is obtained as follows:

\begin{align*} \dot{\mathbf{V}}_{C} & =-{\mathbf{z}_{1}}^{T}\boldsymbol{\mathcal{B}}\mathbf{z}_{1}+\dot{\overline{\mathbf{z}}}^{T}\mathbf{P}\overline{\mathbf{z}}+\overline{\mathbf{z}}^{T}\mathbf{P}\dot{\overline{\mathbf{z}}}+\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\big(\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}}\big)\\& =-{\mathbf{z}_{1}}^{T}\boldsymbol{\mathcal{B}}\mathbf{z}_{1}+\left(\boldsymbol{\mathcal{B}}\overline{\mathbf{z}}\right)^{T}\mathbf{P}\overline{\mathbf{z}}+\left(-{\mathbf{K}_{d}}{\mathbf{z}_{2}}-\left({\mathbf{K}_{p}}+{\boldsymbol{\mathcal{B}}^{2}}\right){\mathbf{z}_{1}}+\unicode{x1D6C8}+{\unicode{x1D6C5}_{a}}\right)^{T}\mathbf{P}\overline{\mathbf{z}}\\& \quad +\overline{\mathbf{z}}^{T}\mathbf{P}\boldsymbol{\mathcal{B}}\overline{\mathbf{z}}+\overline{\mathbf{z}}^{T}\mathbf{P}\left(-\mathbf{K}_{d}\mathbf{z}_{2}-\left(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2}\right)\mathbf{z}_{1}+\unicode{x1D6C8}+\unicode{x1D6C5}_{a}\right)+\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\big(\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}}\big)\\& =\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}\left(-\mathbf{K}_{d}\mathbf{z}_{2}-\left(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2}\right)\mathbf{z}_{1}+\unicode{x1D6C8}+\unicode{x1D6C5}_{a}\right)+\overline{\mathbf{z}}^{T}\left(\boldsymbol{\mathcal{B}}^{T}\mathbf{P}+\mathbf{P}\boldsymbol{\mathcal{B}}^{T}\right)\overline{\mathbf{z}}+\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\big(\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}}\big) \end{align*}

For the sake of simplicity, denote $\overline{\mathbf{z}}^{T}\mathbf{P}$ by a vector $\mathbf{h}$ ,

\begin{equation*} \mathbf{h}=\mathbf{P}^{T}\overline{\mathbf{z}} \end{equation*}

By using Eq. (50), $\dot{\mathbf{V}}_{C}$ can be written as follows:

(52) \begin{equation}\dot{\mathbf{V}}_{C}=\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}(-\mathbf{k}_{d}\mathbf{z}_{2}-(\mathbf{k}_{p}+\boldsymbol{\mathcal{B}}^{2})\mathbf{z}_{1})+\mathbf{2}\mathbf{h}^{T}(\unicode{x1D6C8}+\unicode{x1D6C5}_{\textrm{a}})-\overline{\mathbf{z}}^{T}\boldsymbol{\Upsilon }\overline{\mathbf{z}}+\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\big(\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}}\big)\end{equation}

Then, the following control effort $\unicode{x1D6C5}_{a}$ and the adaptation law are chosen:

(53) \begin{equation}\unicode{x1D6C5}_{a}=-\mathbf{h}-\hat{\unicode{x1D6C8}}\left\{\begin{array}{l@{\quad}l} \dfrac{\mathbf{h}}{\left\| \mathbf{h}\right\| } & if \left\| \mathbf{h}\right\| \gt \epsilon \\[15pt] \dfrac{\mathbf{h}}{\epsilon } & if \left\| \mathbf{h}\right\| \leq \epsilon \end{array}\right.\end{equation}

and

(54) \begin{equation}\dot{\hat{\unicode{x1D6C8}}}=\mathbf{2}\mathbf{G}_{a}\mathbf{h}\end{equation}

in which $\epsilon$ is a threshold width on the $\mathbf{h}$ variable. If $\epsilon$ is suitably chosen to be larger than the measurement noise amplitude, this approximation significantly reduces the output chattering. Therefore, the time derivative of the Lyapunov function becomes,

(55) \begin{equation}\dot{\mathbf{V}}_{C}\leq -\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}\left(\mathbf{K}_{d}\mathbf{z}_{2}+\left(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2}\right)\mathbf{z}_{1}\right)-2\mathbf{h}^{T}\mathbf{h}-\overline{\mathbf{z}}^{T}\boldsymbol{\Upsilon }\overline{\mathbf{z}}-\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\dot{\unicode{x1D6C8}}\end{equation}

By choosing $\mathbf{K}_{d}$ and $\mathbf{K}_{p}$ as positive-definite matrices, $-\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}(\mathbf{K}_{d}\mathbf{z}_{2}+(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2})\mathbf{z}_{1})$ becomes negative definite. In the following, the stability analysis is performed under two different assumptions.

Assumption 1: Assuming that the uncertainties are randomly large and change slowly with time, $\dot{\unicode{x1D6C8}}$ can be omitted.

\begin{equation*} \dot{\mathbf{V}}_{C}\lt \mathbf{0} \end{equation*}

Assumption 2: Assuming that the uncertainties are arbitrarily large and quickly change over time. In this case, the following scenarios can occur:

\begin{equation*} if\ \mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\dot{\unicode{x1D6C8}}\gt 0\Rightarrow \dot{\mathbf{V}}_{C}\lt 0 \end{equation*}

In the worst case, i.e., $\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\dot{\unicode{x1D6C8}}\lt 0$ , inequality (55) can be written as follows:

\begin{equation*} \dot{\mathbf{V}}_{C}\leq -\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}\left(\mathbf{K}_{d}\mathbf{z}_{2}+\left(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2}\right)\mathbf{z}_{1}\right)-2\mathbf{h}^{T}\mathbf{h}-\overline{\mathbf{z}}^{T}\boldsymbol{\Upsilon }\overline{\mathbf{z}}+\boldsymbol{\Pi } \end{equation*}

in which $\boldsymbol{\Pi }=-\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\dot{\unicode{x1D6C8}}\gt 0$ . In this case, $\boldsymbol{\Pi }$ may make $\dot{\mathbf{V}}$ positive definite. A circumstance in which $\mathbf{z}\neq \mathbf{0}$ , by choosing controller parameters $\mathbf{K}_{p}, \mathbf{K}_{d}$ , and $\boldsymbol{\mathcal{B}}$ as large as possible, positive term $\boldsymbol{\Pi }$ can become negligible and asymptotical stability of the closed-loop system is guaranteed. On the other hand, when $\mathbf{z}=\mathbf{0}$ , by choosing the adaptation gain $\mathbf{G}_{a}$ large enough, the positive term $\boldsymbol{\Pi }$ becomes very small. It should be noted that larger $\mathbf{G}_{a}$ results in a faster adaptation, but also increases the control effort.

The estimated uncertainty vector can be computed by integrating Eq. (54) as follows:

\begin{equation*} \hat{\unicode{x1D6C8}}=2\mathbf{G}_{a}\left(\int \mathbf{h}dt\right) \end{equation*}

Therefore, the proposed controller (PC) is written as follows:

(56) \begin{equation}\unicode{x1D6D5}=\hat{\mathbf{H}}\left(\ddot{\mathbf{x}}_{d}+\mathbf{K}_{d}\dot{\mathbf{e}}_{x}+\mathbf{K}_{p}\mathbf{e}_{x}\right)+\hat{\mathbf{H}}\left(-\mathbf{P}^{T}\overline{\mathbf{z}}-\mathbf{2}\mathbf{G}_{a}\left(\int \mathbf{h}dt\right)sign\left(\mathbf{h}\right)\right)+\hat{\mathbf{C}}\dot{\mathbf{x}}+\hat{\mathbf{G}}\end{equation}

where

\begin{align*} sign\left(\mathbf{h}\right)=\left\{\begin{array}{l@{\quad}l} \dfrac{\mathbf{h}}{\left\| \mathbf{h}\right\| } & if \left\| \mathbf{h}\right\| \gt \epsilon \\[15pt] \dfrac{\mathbf{h}}{\epsilon} & if \left\| \mathbf{h}\right\| \leq \epsilon \end{array}\right. \end{align*}

7. Observer-based compensator of the ankle-resistance torques

7.1. The hysteresis loop of the ankle-resistance torque

In a case where the robot is applied to passive ankle rehabilitation, the dynamic model can be written as follows [60]:

(57) \begin{equation}\mathbf{H}(\mathbf{x})\ddot{\mathbf{x}}+\mathbf{C}(\mathbf{x},\dot{\mathbf{x}})\dot{\mathbf{x}}+\mathbf{G}(\mathbf{x})-\unicode{x1D6D5}_{r}=\mathbf{H}(\mathbf{x})\ddot{\mathbf{x}}+\mathbf{C}(\mathbf{x},\dot{\mathbf{x}})\dot{\mathbf{x}}+\mathbf{G}(\mathbf{x})-{\mathbf{J}_{a}}^{-T}\boldsymbol{\mathcal{F}}_{r}=\unicode{x1D6D5}\end{equation}

where $\unicode{x1D6D5}$ denotes the control input, and $\boldsymbol{\mathcal{F}}_{r}$ and $\unicode{x1D6D5}_{r}$ are the ARTs in the task space and joint space, respectively. As shown in ref. [51], the torque-angle curves of the ankle joint behave like a hysteresis loop (Fig. 9). Thus, by using a linear predictive model (LPM), an estimate of the ART exerted on the robot’s MP during passive exercises can be written as:

(58) \begin{equation}\mathcal{F}_{rl}=\begin{cases} K _{{rl_{df}}}\theta _{a}+\mathcal{F}_{{rl_{0}}}\theta _{a}\gt 0\\ K _{{rl_{pf}}}\theta _{a}+\mathcal{F}_{{rl_{0}}}\theta _{a}\lt 0 \end{cases}\end{equation}

where $\mathcal{F}_{rl}$ is the linear model of $\mathcal{F}_{r}$ . Moreover, $K _{{rl_{df}}}$ and $K _{{rl_{pf}}}$ represent the ankle stiffness during dorsiflexion and plantarflexion movements, respectively. $\theta _{a}$ represents the sagittal plane angular displacement of the ankle, and $\mathcal{F}_{{rl_{0}}}$ is a bias term giving the resistance torque when the angle of the ankle joint is zero. By adjusting the $K _{{rl_{df}}}$ , $K _{{rl_{pf}}}$ , and $\mathcal{F}_{{rl_{0}}}$ coefficients [61], it is possible to get a primary estimate of each patient’s resistance torque. Therefore, when the robot is used for rehabilitation, an estimate of the required actuator torque to compensate for the ART can be obtained as follows:

(59) \begin{equation}\unicode{x1D6D5}_{rl}={\mathbf{J}_{a}}^{-T}\left(\begin{array}{c} \mathcal{F}_{rl}\\ \delta \\ \delta \end{array}\right)\end{equation}

The main challenge of this modeling is its simplicity, and determining the appropriate parameters ( $K _{{r_{df}}}$ , $K _{{r_{pf}}}$ , and $\mathcal{F}_{{r_{0}}}$ ) for each patient can become an erroneous task and make the modeling error non-negligible. Moreover, even by choosing suitable LPM parameters, there will always be a mismatch due to the kinematic uncertainties in ${\mathbf{J}_{a}}^{-T}$ . To compensate for this problem, it is proposed to exploit an observer together with one in Eq. 59 as a robustifying term. This observer computes the assistive torque required to overcome the error between the estimated resistant torque $\unicode{x1D6D5}_{rl}$ and its actual value $\unicode{x1D6D5}_{r}$ indicated by $\unicode{x1D6D5}_{\varrho }$ , i.e., $\unicode{x1D6D5}_{\varrho }=\unicode{x1D6D5}_{rl}-\unicode{x1D6D5}_{r}$ . So, in the next step, estimating unknown bounded torque $\unicode{x1D6D5}_{\varrho }$ becomes of interest. The proposed observer in this research has the ability to estimate $\unicode{x1D6D5}_{\varrho }$ and convergence time of the estimation process can be adjusted by tuning the gain of the observer. This method is a model-based approach originally proposed in refs. [62, 63] for serial manipulators and based on the robot dynamic model. The observer is proposed as follows:

(60) \begin{equation}\unicode{x1D6CF}_{\varrho }=\mathbf{K}\left(\hat{\mathbf{H}}\dot{\mathbf{x}}-\int _{0}^{t}\left(\unicode{x1D6D5}-\hat{\mathbf{C}}\dot{\mathbf{x}}-\hat{\mathbf{G}}+\frac{d\hat{\mathbf{H}}}{dt}\dot{\mathbf{x}}+\unicode{x1D6CF}_{\varrho }\right)ds\right)\end{equation}

where $\unicode{x1D6CF}_{\varrho }$ is an estimation of the actual resistant torque $\unicode{x1D6D5}_{\varrho }$ , and $\mathbf{K}$ is the gain of the observer that is designed to be a positive diagonal matrix, i.e., $\mathbf{K}=diag(k_{1},k_{2},k_{3})\gt \mathbf{0}$ . Vector $\unicode{x1D6CF}_{\varrho }$ can be computed by using the measured $\mathbf{X}$ and $\dot{\mathbf{X}}$ and the commanded torque control $\unicode{x1D6D5}$ .

Figure 9. Representative torque-angle curves (hysteresis loops).

7.2. Stability analysis of the observer and the proposed controller

Consider the following dynamic equations of the robot in the presence of $\unicode{x1D6D5}_{\varrho }$ :

(61) \begin{equation}\mathbf{H}\left(\mathbf{x}\right)\ddot{\mathbf{x}}+\mathbf{C}\left(\mathbf{x},\dot{\mathbf{x}}\right)\dot{\mathbf{x}}+\mathbf{G}\left(\mathbf{x}\right)=\unicode{x1D6D5}+\unicode{x1D6D5}_{\varrho }\end{equation}

The ODEs of the closed-loop system can be written as:

(62) \begin{equation}\dot{\mathbf{z}}_{1}=\mathbf{z}_{2}\end{equation}

(63) \begin{equation}\dot{\mathbf{z}}_{2}=-\mathbf{K}_{d}\mathbf{z}_{2}-\mathbf{K}_{p}\mathbf{z}_{1}+\unicode{x1D6C5}_{a}+\unicode{x1D6C8}+\unicode{x1D6D5}_{\varrho }\end{equation}

For the convergence analysis, let’s consider the following variable, which denotes the error dynamics of $\unicode{x1D6D5}_{\varrho }$ :

(64) \begin{equation}\dot{\tilde{\unicode{x1D6CF}}}=\dot{\unicode{x1D6D5}}_{\varrho }-\dot{\unicode{x1D6CF}}_{\varrho }\end{equation}

By substituting $\unicode{x1D6D5}$ from Eq. (61) into Eq. (60), the dynamic relation between the resistance torque $\unicode{x1D6D5}_{\varrho }$ and its estimated value $\unicode{x1D6CF}_{\varrho }$ is obtained as follows:

(65) \begin{equation}\dot{\unicode{x1D6CF}}_{\varrho }=\mathbf{K}\left(\unicode{x1D6D5}_{\varrho }-\unicode{x1D6CF}_{\varrho }\right)=\mathbf{K}\tilde{\unicode{x1D6CF}}\end{equation}

Thus, by taking a suitable value for $\mathbf{K}$ , the estimated value $\unicode{x1D6CF}_{\varrho }$ converges to the accurate value $\unicode{x1D6D5}_{\varrho }$ . For the stability analysis of both the observer and controller, the following positive-definite Lyapunov function is considered:

(66) \begin{equation}\mathbf{V}_{T}=\frac{1}{2}\tilde{\unicode{x1D6CF}}^{T}\tilde{\unicode{x1D6CF}}+\frac{1}{2}{\mathbf{z}_{1}}^{T}\mathbf{z}_{1}+\overline{\mathbf{z}}^{T}\mathbf{P}\overline{\mathbf{z}}+\frac{1}{2}\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\mathbf{E }\end{equation}

Differentiating Eq. (66) with respect to time yields,

(67) \begin{align}\dot{\mathbf{V}}_{T} & =\tilde{\unicode{x1D6CF}}^{T}\dot{\tilde{\unicode{x1D6CF}}}-{\mathbf{z}_{1}}^{T}\boldsymbol{\mathcal{B}}\mathbf{z}_{1}+\dot{\overline{\mathbf{z}}}^{T}\mathbf{P}\overline{\mathbf{z}}+\overline{\mathbf{z}}^{T}\mathbf{P}\dot{\overline{\mathbf{z}}}+\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\big(\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}}\big)\nonumber\\&=\tilde{\unicode{x1D6CF}}^{T}\left(\dot{\unicode{x1D6D5}}_{\varrho }-\dot{\unicode{x1D6CF}}_{\varrho }\right)+\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}\left(-\mathbf{K}_{d}\mathbf{z}_{2}-\left(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2}\right)\mathbf{z}_{1}+\unicode{x1D6C8}+\unicode{x1D6C5}_{a}\right)\nonumber\\& +\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P} \unicode{x1D6D5}_{\varrho }+\overline{\mathbf{z}}^{T}(\boldsymbol{\mathcal{B}}^{T}\mathbf{P}+\mathbf{P}\boldsymbol{\mathcal{B}}^{T})\overline{\mathbf{z}}+\mathbf{E }^{T}{\mathbf{G}_{a}}^{-1}\big(\dot{\hat{\unicode{x1D6C8}}}-\dot{\unicode{x1D6C8}}\big)\end{align}

Then by choosing the control law and adaptation rule in Eq. (53) and Eq. (54), respectively, and substituting Eq. (65) into Eq. (67), $\dot{\mathbf{V}}_{T}$ is obtained as

(68) \begin{equation}\dot{\mathbf{V}}_{T}\leq \tilde{\unicode{x1D6CF}}^{T}\dot{\unicode{x1D6D5}}_{\varrho }+\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P} \unicode{x1D6D5}_{\varrho }-\tilde{\unicode{x1D6CF}}^{T}\mathbf{K}\tilde{\unicode{x1D6CF}}-\mathbf{2}\overline{\mathbf{z}}^{T}\mathbf{P}\left(\mathbf{K}_{d}\mathbf{z}_{2}+\left(\mathbf{K}_{p}+\boldsymbol{\mathcal{B}}^{2}\right)\mathbf{z}_{1}\right)-2\mathbf{h}^{T}\mathbf{h}-\overline{\mathbf{z}}^{T}\boldsymbol{\Upsilon }\overline{\mathbf{z}}+\boldsymbol{\Pi }\end{equation}

Assuming that the actual resistance torque $\unicode{x1D6D5}_{\varrho }$ varies slowly with time, i.e., $\dot{\unicode{x1D6D5}}_{\varrho }\approx 0$ , and the observer provides an appropriate estimation of $\unicode{x1D6D5}_{\varrho }$ , $\dot{\mathbf{V}}_{T}$ can become negative definite by choosing suitable control parameters. If the $\dot{\unicode{x1D6D5}}_{\varrho }\approx 0$ condition is not met, the Lyapunov function can become negative definite by choosing $\mathbf{K}$ as large as possible. However, the condition $\dot{\unicode{x1D6D5}}_{\varrho }\approx 0$ holds as the hysteresis loops of the ARTs do not change rapidly. As a result, the following robust linear predictive model (RLPM) is proposed to overcome the ARTs:

(69) \begin{equation}\unicode{x1D6D5}_{rlpm}=\unicode{x1D6CF}_{\varrho }+{\mathbf{J}_{\mathbf{a}}}^{-T}\boldsymbol{\mathcal{F}}_{rl}\end{equation}

8. Results and discussion

Extensive computer co-simulations have been accomplished to study the overall efficiency of the proposed algorithms for passive ankle rehabilitation purposes. For this purpose, simulations have been done without considering the ARTs in the first step. After ensuring the stability and high performance of the PC, the robot-human interactions arising from passive rehabilitation exercises have been taken into account, and the stability and performance of the proposed schemes have been discussed.

8.1. Performance analysis and co-simulation setup

After determining the robot material (aluminum alloy), the physical parameters given in Tables I and II have been obtained in 3D CAD software. In addition, $\mathbf{c}_{i}$ , $\mathbf{p}_{i}$ , $\unicode{x1D6CC}_{i}$ , and $\unicode{x1D6DD}$ vectors have been determined in the initial configuration of the robot in the 3D CAD software. An XML file of non-over constrained version 3-RCC of RSPM is generated and then imported into a multibody environment. The co-simulation procedure for the PC and RLPM is shown in Fig. 10. The motion variables of the MP are sent to the control system by the sensor block in the multibody environment and sent to the controller blocks for running computer-based simulations. To discover the performances of the proposed schemes, a sequence of simulations is performed, and the PC is benchmarked by comparing its tracking performance in the presence of model uncertainties, noise measurements with an amplitude of %0.02 peak values of the original signals, and external disturbances with three different control topologies, namely, a CTC, sliding mode control (SMC), and adaptive sliding mode control using the time delay estimation technique [47] (SMC + TDE).

Table I. Link inertia ( $\textrm{gmm}^{2}$ ) taken at the center of mass and mass parameters (g) of the RSPM.

Table II. Dimensions of the RSPM.

Figure 10. Co-simulation processes for proposed schemes.

8.2. Simulation results and discussion

The performance and robustness of the proposed schemes are explored through the presentation and discussion of simulation results in this subsection. To do this, a desired trajectory in Cartesian space, as shown in Fig. 11, should be tracked by the MP in different scenarios. In the first step, the simulation has been performed under ideal conditions, where dynamic uncertainties, measurement noise, and disturbances are ignored. The controllers are tuned in such a way that they all have the same level of efficiency. The control parameters are given in Table III where $\unicode{x1D6C4}(\boldsymbol{z})$ denotes the estimated upper bound of the uncertainty vector and should be chosen to satisfy the inequality $\unicode{x1D6C4}(\boldsymbol{z})\geq \hat{\unicode{x1D6C8}}(\mathbf{z})$ . Moreover, $\mathbf{A }$ is a diagonal matrix, which is chosen as a positive, definite constant matrix and refers to the sliding surface $\unicode{x1D6D4}=\mathbf{A }\mathbf{z}_{1}+\mathbf{z}_{2}$ of the SMC method. To study the behavior of the closed-loop system in perturbed conditions where dynamic uncertainties, measurement noise, and external disturbances exist, a co-simulation is performed to command the RSPM to follow the same trajectory in Fig. 11 with the inertia parameters offset±30% and the kinematics parameters $\alpha _{1},\alpha _{2}, \gamma,\eta _{1},\eta _{2}$ , and $\eta _{3}$ offset ± 15% from their real values. In addition, noise measurement with an amplitude of % 0.01 peaks of the original signals is added to the control system, taking into account that all measured variables of vector x are contaminated. To guarantee the robustness of the controllers, the simulations were carried out in the presence of both lower and upper bounds of uncertainty. To compare the performance of the controllers, the norms of the error vectors are shown in Fig. 12, and when compared to other schemes, the PC provides a lower error norm.

For more detail, the errors of Cartesian space states are shown in Figs. 13 and 14 for upper and lower bound uncertainties, respectively. Moreover, as shown in Figs. 15 and 16, although the tracking error is not much increased compared to that of the system without noise, the required actuator torques of the TDE controller to carry out such a maneuver are very oscillatory. This is due to the fact that the TDE approach requires acceleration of Cartesian variables to accommodate the required tracking performance, and although the amplitude of noise is small, the signature of noise is apparently seen in the required actuator torques. In fact, generating such control input is infeasible in practice, and hence, reaching the tracking performance of the TDE method is a challenging task.

Table III. Controller parameter settings.

Figure 11. Desired trajectory for plantar/dorsi flexion configurations of MP.

Figure 12. Time histories of the norm of the tracking errors in the presence of (a) lower and (b) upper bound uncertainties.

Figure 13. The closed-loop tracking performance of the RSPM for different controllers in the presence of upper bound uncertainties.

Figure 14. The closed-loop tracking performance of the RSPM for different controllers in the presence of lower bound uncertainties.

Figure 15. Control input time histories (actuator torques) for different controllers during co-simulations in the presence of lower bound uncertainties and noise measurements.

Figure 16. Control input time histories (actuator torques) for different controllers during co-simulations in the presence of upper bound uncertainties and noise measurements.

Figure 17. (a) Time histories of disturbances applied to the robot; (b) Time histories of the norm of the tracking errors in the presence of external disturbances.

The disturbance suppression ability of controllers was evaluated by applying the disturbances shown in Fig. 17(a) to the robot in joint space. As the robot interacts with humans, disturbances with quick dynamic changes are more common and relevant because of human-induced impulses and human movement patterns. As shown in Fig. 17(b), the PC provides less tracking error and is more reliable. In order to demonstrate the superiority of the proposed algorithm, the performance measurements of all schemes are tabulated in Table IV. It demonstrates that the proposed method significantly decreases tracking errors in perturbed conditions, whereas the other controllers guarantee stability and have suitable performance. Overall, the closed-loop simulations confirm the power and effectiveness of the PC, and this result suggests that the proposed scheme could be applied to a wide range of manipulation tasks. In the next step, the desired trajectory, shown in Fig. 11, has been used for inversion and eversion movements. Fig. 18 shows a visual view of the animation for the RSPM movements enforced by the PC for dorsiflexion, plantarflexion, and inversion and eversion configurations.

Table IV. Comparative controller performances during trajectory tracking in different scenarios.

Figure 18. A virtual prototype model (3-RCC spherical parallel manipulator) of ankle exercise under passive exercise modes: (a) Inversion/eversion; (b) Dorsiflexion/plantarflexion.

In the following, co-simulations have been performed in four different scenarios. It is assumed that the robot should generate passive rehabilitation movements in the presence of ART, which is shown in Fig. 19. For this, it is necessary to overcome the weight and inertia of the ankle as well as the ART. According to [64], considering that the weight and inertia of the foot are about 1.5% of the human body weight, it can be considered the parametric uncertainty that has already been taken into account. The parameters of the linear model in Eq. (58) are chosen in four different cases, and the observer is expected to estimate the error between the mentioned linear model and the actual resistance torque applied to the robot in the joint space. Then, to evaluate the performance and necessity of the observer, the tracking error is compared in the presence and absence of the observer. The selection of parameters in four cases is shown in Table V, and to better understand the differences between each case, the selected parameters of the LPM in each case and the actual resistance torque of the ankle are shown in Fig. 19. The gain of the observer $\mathbf{K}$ is selected to be equal to 50. As shown in Fig. 20(a) in the first case, where the parameters of the LPM are chosen with good accuracy, the tracking error can be neglected, and there are not many differences between the LPM and the RLPM in terms of efficiency and accuracy. Since the exact selection of the ankle parameters in Eq. (58) is not possible for every patient, in the second case, the value of $\mathcal{F}_{{r_{0}}}$ is changed to six. This change is also made to measure the safety and robustness of the proposed method against instability. As shown in Fig. 20(b), the LPM stabilizes the closed-loop system, but the tracking error is not negligible. Obviously, this issue is due to the error of selecting $\boldsymbol{\mathcal{F}}_{{r_{0}}}$ from its actual value. However, as depicted in Fig. 20(b), the tracking error is reduced significantly by using RLPM, and according to Fig. 21, despite the measurement noise and parametric uncertainties, the predictive observer is able to estimate the torque required to overcome the resistance torque in the joint space.

Table V. Parameters of the linear predictive model of ankle resistance torque in different scenarios.

Figure 19. Selection of the range of ankle stiffness parameters and actual resistance torque of the ankle.

Figure 20. Comparison of the norm of the tracking error vector in (a) case 1 and (b) case 2.

Figure 21. Time histories of actual resistance torque $\unicode{x1D6D5}_{\varrho }$ and estimated resistance torque $\unicode{x1D6CF}_{\varrho }$ in (a) case 1 and (b) case 2.

In the third and fourth cases, the robot became unstable in the absence of the predictive observer. As it is clear from Fig. 22, the tracking error is less than 0.6 degrees by using the RLPM, and according to Fig. 23, the observer is able to estimate the resistance torque in these scenarios. The control efforts in all four cases are also shown in Fig. 24. By using the proposed method, the stability and high tracking efficiency of the closed-loop system can be guaranteed by avoiding the use of high-gain controllers, which lead to higher control effort, power consumption, and robot stiffness. Moreover, according to Fig. 19, the proposed method is largely robust to the error between the actual and selected parameters of the LPM.

Figure 22. Comparison of the norm of the tracking error vector in cases 3 and 4.

Figure 23. Time histories of actual resistance torque $\unicode{x1D6D5}_{\varrho }$ and estimated resistance torque $\unicode{x1D6CF}_{\varrho }$ in (a) case 3 and (b) case 4.

Figure 24. Control efforts in the presence of ankle resistance torque and different selections of joint ankle stiffness parameters in four cases.

9. Conclusion

In this paper, a spherical parallel robot is presented in aggregate with powerful STBC equipped with a robust linear predictive model for offering trajectory tracking capabilities in the presence of dynamic uncertainties, external disturbances, measurement noises, and ART. Considering the ART will become of interest when passive and assistive ankle rehabilitation exercises ought to be carried out. To attain the cited goals, the dynamic model of the robot has been extracted and validated using the Newton-Euler technique and screw theory. The advantage of this method is that, at the same time as being simple, it could be used for different spherical parallel manipulators with complicated kinematic structures. By using this dynamic model, advanced controllers such as CTC, SMC, and SMC + TDE are designed in Cartesian space, and so as to compare with the PC, which is a STBC, co-simulations are performed in different scenarios. Co-simulation outcomes show that the maximum tracking error norms have been decreased by 12% and 34%, respectively, with respect to the SMC and PC in comparison to the CTC in the presence of upper bound dynamic uncertainties. Also, in the presence of external disturbances, the maximum tracking error norms of the SMC and PC are 12% and 42%, much less than the CTC, respectively. The suitable robustness of the proposed approach toward uncertainties, measurement noise, and external disturbances indicates that the parallel robot can be effectively used in situations where trajectory tracking is a necessity.

The effectiveness of the proposed RLPM has been evaluated through passive rehabilitation exercises. The results show that, although different choices of ankle joint stiffness result in fast time-varying resistance torque in the joint space, the proposed technique is capable of providing a truthful prediction of the ART. With the assistance of this technique, not only is system stability guaranteed, but additionally, the tracking error may be significantly reduced in the presence of the ART. The validations show that RSPM can be successfully commercialized for utilization in ankle rehabilitation applications. In future work, the actuator’s dynamic behavior and passive joint friction will be extracted once a physical prototype is developed.