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 [Reference Taghvaeipour, Angeles and Lessard59], 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 [Reference Taghirad33]:

(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 [Reference Staicu36] ( $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 [Reference Taghvaeipour, Angeles and Lessard59]:

(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 [Reference Saglia, Tsagarakis, Dai and Caldwell60]. 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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61], 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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61]. 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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61]:

(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 [Reference De Luca, Albu-Schaffer, Haddadin and Hirzinger62] ( $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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61], 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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61] 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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61], 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 [Reference Saglia, Tsagarakis, Dai and Caldwell60]:

(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. [Reference Chung, van Rey, Bai, Roth and Zhang51], 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 [Reference Roy, Krebs, Bever, Forrester, Macko and Hogan61], 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. [Reference De Luca, Albu-Schaffer, Haddadin and Hirzinger62, Reference Haddadin, De Luca and Albu-Schäffer63] 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}