*Article*

## **Cable Tension Analysis Oriented the Enhanced Stiffness of a 3-DOF Joint Module of a Modular Cable-Driven Human-Like Robotic Arm**

**Kaisheng Yang 1,2,3,\*, Guilin Yang 2,\*, Chi Zhang 2, Chinyin Chen 2, Tianjiang Zheng 2, Yuguo Cui 1 and Tehuan Chen 1,4**


Received: 18 November 2020; Accepted: 7 December 2020; Published: 11 December 2020

**Abstract:** Inspired by the structure of human arms, a modular *cable-driven human-like robotic arm* (CHRA) is developed for safe human–robot interaction. Due to the unilateral driving properties of the cables, the CHRA is redundantly actuated and its stiffness can be adjusted by regulating the cable tensions. Since the trajectory of the *3-DOF joint module* (3DJM) of the CHRA is a curve on Lie group SO(3), an enhanced stiffness model of the 3DJM is established by the covariant derivative of the load to the displacement on SO(3). In this paper, we focus on analyzing the how cable tension distribution problem oriented the enhanced stiffness of the 3DJM of the CHRA for stiffness adjustment. Due to the complexity of the enhanced stiffness model, it is difficult to solve the cable tensions from the desired stiffness analytically. The problem of *stiffness-oriented cable tension distribution* (SCTD) is formulated as a nonlinear optimization model. The optimization model is simplified using the symmetry of the enhanced stiffness model, the rank of the Jacobian matrix and the equilibrium equation of the 3DJM. Since the objective function is too complicated to compute the gradient, a method based on the genetic algorithm is proposed for solving this optimization problem, which only utilizes the objective function values. A comprehensive simulation is carried out to validate the effectiveness of the proposed method.

**Keywords:** cable-driven robots; human-like robotic arms; human–robot interactions; stiffness adjustment; cable tension analysis

## **1. Introduction**

Unlike conventional robots that work in structured environments, safe human–robot interactions have been a key element for the robots that work in unstructured and unpredictable environments. Inspired by the structure of human arms, a modular *cable-driven human-like robotic arm* (CHRA) is developed for safe human–robot interaction, which employs cables to mimic the functionality of the human muscles. As shown in Figure 1, the CHRA consists of a shoulder joint, an elbow joint and a wrist joint in series, where the shoulder joint and the wrist joint have *three degrees of freedom* (3-DOF) and the elbow joint has *one degree of freedom* (1-DOF). The CHRA employs lightweight cables to drive the rigid links and the cables can be wound into winches mounted onto the base of the CHRA. A *variable-stiffness device* (VSD) is designed and placed along with each driving cable to increase the

flexibility of the cables. With these arrangements, the CHRA has flexibility advantages [1], low moving mass [2], large workspace [3] and high payload-to-weight ratio [4]. Due to these advantages, the CHRA is intrinsically safe for human–robot interactions. The proposed CHRA and its joint modules are one kind of *cable-driven mechanisms* (CDMs). In the last decades, various CDMs have been designed for various applications, such as medical robots [5–7], rehabilitation robots [8–12], inspection and repair [13–15], and moving payloads [16–19].

**Figure 1.** Concept design of the modular cable-driven human-like robotic arm: (**a**) structure of a human arm, (**b**) structure of the *cable-driven human-like robotic arm* (CHRA).

Unlike rigid links, the cables have the unilateral driving property (i.e., they can only pull, but can not push). Hence, the CHRA and its joint modules are redundantly actuated. For a given pose, a lot of cable tension solutions are feasible for the CHRA. Furthermore, the stiffness of the CHRA can be adjusted by regulating the cable tensions. The characteristic of variable stiffness increases the flexibility and safety of the CHRA.

The problem of *stiffness-oriented cable tension distribution* (SCTD) of a CDM aims at finding the optimal cable tensions to achieve a desired stiffness. However, in the last decades, most studies on cable tension distribution have been carried out to minimize the cable tensions to reduce the energy consumption [20–27]. In [28], the SCTD problem of a CDM was studied by formulating it as an optimization model, and a gradient-projection-based algorithm was presented to solve the optimization problem. However, it utilized the determinant of the stiffness matrix as the objective function. With this method, the desired stiffness cannot be achieved accurately. In [29], the SCTD problem of a 3-DOF cable-driven spherical mechanism was studied and it established the optimization model with all entries of the stiffness matrix, not only the determinant of the stiffness matrix. However, it employed the conventional stiffness model of the CDM, which was derived by the conventional differential formula of the load to the displacement. Since the trajectory of the 3-DOF cable-driven spherical mechanism is a curve on Lie group SO(3) and SO(3) is nonlinear, the stiffness model based on the conventional differential formula is not exactly accurate on SO(3). Its stiffness should be evaluated by the variation of its load against its displacement on SO(3).

In order to study the SCTD problem of the proposed CHRA, the SCTD problems of the 1-DOF and 3-DOF joint modules should be studied first. Due to the simple structure of the 1-DOF joint module, its SCTD problem can be solved easily. For the *3-DOF joint module* (3DJM) of the CHRA, its trajectory is a curve on SO(3). We derive the stiffness model of the 3DJM by the covariant derivative of the load to the displacement on SO(3), and call it the enhanced stiffness model of the 3DJM. In this paper, we focus on investigating how the cable tension distribution problem oriented the enhanced stiffness of the 3DJM. Due the complexity of the enhanced stiffness model, it is difficult to solve the cable tensions from the desired stiffness of the 3DJM analytically. We formulate the SCTD problem of the 3DJM as a nonlinear constrained optimization problem. Due to the symmetry of the enhanced stiffness model

of the 3DJM, the desired stiffness matrix can be transformed to a diagonal matrix by an orthogonal transformation. Based on the analysis of the rank of the Jacobian matrix of the 3DJM, the six cables can be divided into two groups: one group with three cables for position adjustment by regulating the cable lengths, and another group with the remaining three cables for stiffness adjustment by regulating the cable tensions. In this manner, the position and stiffness of the 3DJM can be adjusted simultaneously. Furthermore, three cable tensions for position adjustment can be expressed by another three cable tensions for stiffness adjustment using the equilibrium equation of the 3DJM. That means the six decision variables of the optimization model can be reduced to three. Since the objective function of the optimization model is too complicated to compute the gradient, a direct optimization method based on the genetic algorithm is proposed for solving the optimization problem, which only utilizes the objective function values. A comprehensive simulation is carried out to validate the effectiveness of the proposed method. The results show that the proposed method provides an accurate and efficient way to adjust the stiffness of the 3DJM by regulating the cable tensions. In summary, the main contribution of this paper is solving the cable tension distribution problem with the enhanced stiffness model of the 3DJM using a method based on the genetic algorithm.

#### **2. Enhanced Stiffness Model of the 3DJM**

As shown in Figure 2, the 3-DOF joint module (3DJM) of the CHRA was made up with a base, a moving-platform and a passive spherical joint connecting them. Six cables were employed to drive the moving-platform and a *variable-stiffness device* (VSD) was placed along with each driving cable. For routing the cables, six small holes were drilled in the moving-platform and the base, denoted by *Ai* and *Bi* (*i* = 1, 2, ... , <sup>6</sup>), respectively. In this design, *A*2*A*3 = *A*4*A*5 = *A*6*A*1 = *lA*, *B*1*B*2 = *B*3*B*4 = *B*5*B*6 = *lB*, *A*1*A*2 = *A*3*A*4 = *A*5*A*6 = *eA*, *B*2*B*3 = *B*4*B*5 = *B*6*B*1 = *eB*, where *eA* ≈ 0 and *eB* ≈ 0. The distance between the moving-platform and the spherical joint is denoted as *hA* ∈ R, and the distance between the base and the spherical joint is denoted as *hB* ∈ R. Each cable was actuated by a cable driving unit.

To describe the pose of the 3DJM, two frames were attached to the moving-platform (named moving frame {A}) and the base (named base frame {B}), which were both located at the center of the joint *O*. With the two frames, the pose of the 3DJM can be represented by a rotation matrix *R* ∈ *SO*(3). Thus, the motion trajectory of the 3DJM was a parameterized curve *R*(*t*) on *SO*(3). The enhanced stiffness model of the 3DJM was established by using the covariant derivative of the load to the displacement on SO(3) [30].

 **Figure 2.** Concept design of the *3-DOF joint module* (3DJM).

#### *2.1. Displacement and Load of the 3DJM on SO(3)*

According to the exponential formula [31], the trajectory curve of the 3DJM *R*(*t*) yields

$$\mathcal{R}(t) = \epsilon^{\sharp(t)},\tag{1}$$

where ˆ *ζ* ∈ *so*(3) is an element of Lie algebra *so*(3) and satisfies

$$\mathcal{J} = \begin{pmatrix} \mathbb{\mathcal{J}}\_1 \\ \mathbb{\mathcal{J}}\_2 \\ \mathbb{\mathcal{J}}\_3 \end{pmatrix} \to \dot{\mathfrak{F}} = \begin{pmatrix} 0 & -\mathbb{\mathcal{J}}\_3 & \mathbb{\mathcal{J}}\_2 \\ \mathbb{\mathcal{J}}\_3 & 0 & -\mathbb{\mathcal{J}}\_1 \\ -\mathbb{\mathcal{J}}\_2 & \mathbb{\mathcal{J}}\_1 & 0 \end{pmatrix}. \tag{2}$$

Here, *ζ* ∈ R<sup>3</sup> is called the coordinates of ˆ*ζ* ∈ *so*(3). The derivative of *R*(*t*) to time *t*, i.e., *R*˙ (*t*), is an element of the tangent space of *SO*(3) at the point *<sup>R</sup>*(*t*), denoted as *<sup>T</sup>R*(*t*)*SO*(3). According to (1), *R* ˙ (*t*) satisfies the following equation

ˆ

$$
\dot{\mathcal{R}}(t) = \mathcal{R}(t)\dot{\mathcal{J}}(t) = \mathcal{R}(t)\dot{\mathcal{w}}(t), \tag{3}
$$

where *ω*(*t*) = ˙ *ζ*(*t*) ∈ *so*(3) is the angular velocity of the 3DJM. Given *e*ˆ*j* (*j* = 1, 2, 3) as the basis of *so*(3), where *e*1 = (1, 0, <sup>0</sup>)*<sup>T</sup>*, *e*2 = (0, 1, 0)*<sup>T</sup>* and *e*3 = (0, 0, <sup>1</sup>)*<sup>T</sup>*, then *R*˙ (*t*) can be expressed as

$$\dot{\mathbf{R}}(t) = \mathbf{R}(t)\hat{\boldsymbol{\omega}}(t) = \mathbf{R}(t)\sum\_{j=1}^{3} \boldsymbol{\omega}^{j}\dot{\mathbf{e}}\_{j} = \sum\_{j=1}^{3} \boldsymbol{\omega}^{j}\mathbf{R}(t)\mathbf{\dot{e}}\_{j} = \sum\_{j=1}^{3} \boldsymbol{\omega}^{j}L\_{j}|\_{\mathbf{R}(t)}.\tag{4}$$

where *<sup>L</sup>j*|*R*(*t*) = *<sup>R</sup>*(*t*)*e*<sup>ˆ</sup>*j* (*j* = 1, 2, 3) are taken as the basis of *<sup>T</sup>R*(*t*)*SO*(3) and *<sup>L</sup>j*|*R*(*t*) are always written as *Lj* for short. From (4), it can be concluded that *<sup>T</sup>R*(*t*)*SO*(3) is isomorphic to *so*(3).

Since *R* ˙ (*t*) describes the velocity of the 3DJM, the instantaneous displacement of the 3DJM can be represented by *δR*, yielding

$$
\delta \mathcal{R}(t) = \dot{\mathcal{R}}(t) \delta t. \tag{5}
$$

That means the instantaneous displacement of the 3DJM can be described on the tangent space of *SO*(3). Similarly, since the moment load *τ*ˆ is an element of *so*<sup>∗</sup>(3), the dual space of *so*(3), by analogy with *ω*<sup>ˆ</sup> ∈ *so*(3) and *R*˙ (*t*) ∈ *<sup>T</sup>R*(*t*)*SO*(3), we define *F*(*t*) = *<sup>R</sup>*(*t*)*τ*<sup>ˆ</sup>(*t*) as an element of the dual space of the tangent space *<sup>T</sup>R*(*t*)*SO*(3) at *<sup>R</sup>*(*t*), called *cotangent space* and denoted as *<sup>T</sup>*<sup>∗</sup>*R*(*t*)*SO*(3). The element of the cotangent space is named as the cotangent vector. Given *λ* ˆ *k* (*k* = 1, 2, 3) as the basis of *so*<sup>∗</sup>(3) and **<sup>Λ</sup>***<sup>k</sup>*|*R*(*t*) = *R*(*t*)*λ*<sup>ˆ</sup> *k* (*k* = 1, 2, 3) as the basis of *<sup>T</sup>*<sup>∗</sup>*R*(*t*)*SO*(3), then *F*(*t*) ∈ *<sup>T</sup>*<sup>∗</sup>*R*(*t*)*SO*(3) satisfies

$$F(t) = \mathcal{R}(t)\hat{\mathbf{r}}(t) = \sum\_{k=1}^{3} \tau\_k \mathcal{R}(t)\hat{\boldsymbol{\lambda}}^k = \sum\_{k=1}^{3} \tau\_k \boldsymbol{\Lambda}^k|\_{\mathcal{R}(t)}.\tag{6}$$

It shows that *<sup>T</sup>*<sup>∗</sup>*R*(*t*)*SO*(3) is isomorphic to *so*<sup>∗</sup>(3). Thus, *F*(*t*) can be employed to represent the load applied on the 3DJM. That means the load applied on the 3DJM can be described on the cotangent space of *SO*(3).

#### *2.2. Enhanced Stiffness Model of the 3DJM on SO(3)*

According to the definition, the stiffness of the 3DJM is evaluated by the variation of its load against its displacement. However, the loads and displacements for different poses of the 3DJM are in different cotangent spaces and tangent spaces of SO(3). Considering the property of the motion of the 3DJM, an affine connection called *Levi-Civita* connection is introduced into *SO*(3) to connect different tangent spaces on SO(3) [32]. Given a curve *R*(*t*) and a vector field *V*(*t*) on *SO*(3), the affine connection specifies how a vector *<sup>V</sup>*(*<sup>t</sup>*0) in the tangent space at point *<sup>R</sup>*(*<sup>t</sup>*0) can be mapped to another vector *Vt*0 (*<sup>t</sup>*1) of the tangent space at some other point *<sup>R</sup>*(*<sup>t</sup>*1). The vector *Vt*0 (*<sup>t</sup>*1) is called the parallel

transport of *<sup>V</sup>*(*<sup>t</sup>*0) along the curve *R*(*t*) [33]. Then, a differentiation operation can be defined on SO(3) as below

$$\left. \nabla\_{\vec{\mathcal{R}}} V \right|\_{\mathcal{R}(t\_0)} = \lim\_{t\_1 \to t\_0} \frac{V^{t\_0}(t\_1) - V(t\_0)}{t\_1 - t\_0},\tag{7}$$

where ∇*R***˙** *<sup>V</sup>R*(*<sup>t</sup>*0) is called the *covariant derivative* of *V*(*t*) along the curve *<sup>R</sup>*(*t*). For two vectors, *V*1 and *V*2, in the tangent spaces of SO(3), ∇*<sup>V</sup>*2*V*<sup>1</sup> is employed to represent the covariant derivative of *V*1 in the direction *V*2. For a real-valued function *f* on *SO*(3), the covariant derivative ∇*V f* is usually written as *V* ◦ *f* . According to the property of covariant derivative, the covariant derivative of the load *F*(*t*) to the instantaneous displacement *δR*(*t*) can be written as

$$\nabla\_{\delta \mathcal{R}(t)} \mathcal{F} = \nabla\_{\mathcal{R}(t) \delta t} \mathcal{F} = \delta t \nabla\_{\sum\_{j=1}^{3} \omega^{j} L\_{j}} \mathcal{F} = \delta t \sum\_{j=1}^{3} \omega^{j} \nabla\_{L\_{j}} \mathcal{F} = \delta t \mathbf{K} \omega = \mathbf{K} \delta \mathsf{J}\_{\prime} \tag{8}$$

and *K* = (∇*<sup>L</sup>*1*F*, ∇*<sup>L</sup>*2*F*, <sup>∇</sup>*<sup>L</sup>*3*<sup>F</sup>*) ∈ R3×<sup>3</sup> is employed as the stiffness of the 3DJM. The components of the stiffness matrix *K*, i.e., *Kjk* (*j*, *k* = 1, 2, <sup>3</sup>), yield

$$\mathbb{K}\_{jk} = \langle \nabla\_{\mathbf{L}\_k} \mathbf{F}, \mathbf{L}\_j \rangle. \tag{9}$$

According to the property of the Levi–Civita connection [30], we have

$$\mathbb{K}\_{\vec{j}k} = \mathbb{L}\_k \diamond \langle \mathbf{F}, \mathbb{L}\_{\vec{j}} \rangle - \langle \mathbf{F}, \nabla\_{\mathbb{L}\_k} \mathbb{L}\_{\vec{j}} \rangle = \mathbb{L}\_k \diamond \tau\_{\vec{j}} - \frac{1}{2} \sum\_{r=1}^3 \tau\_r \gamma\_{kj^\*}^r \tag{10}$$

where the coefficients *γrkj*(*j*, *k*,*<sup>r</sup>* = 1, 2, 3) ∈ R are zero except

$$
\gamma\_{12}^3 = \gamma\_{31}^2 = \gamma\_{23}^1 = 2, \quad \gamma\_{21}^3 = \gamma\_{13}^2 = \gamma\_{32}^1 = -2.\tag{11}
$$

Furthermore, since the 3DJM is a conservative mechanical system, its stiffness matrix *K* can be proved to be symmetric at every pose, even if there is an external load applied on it [30].

#### *2.3. Parametric Stiffness Formulation of the 3DJM*

In Figure 2, *ai* = −−→*OAi*, *bi* = −→*OBi* ∈ R<sup>3</sup> (*i* = 1, 2, ... , 6) are the position vectors of points *Ai* and *Bi*, respectively. The vector from *Ai* to *Bi* along the *ith* cable, denoted as *ci* = −−→*AiBi*, yields

$$\mathbf{c}\_{i} = \mathbf{b}\_{i} - \mathbf{a}\_{i} = \mathbf{c}\_{i}\mathbf{u}\_{i} \in \mathbb{R}^{3},\tag{12}$$

where *ci* = *<sup>c</sup>i* ∈ R is the length of the cable from *Ai* to *Bi* and *ui* = *<sup>c</sup>i*/*ci* ∈ R<sup>3</sup> is the unitary vector of *ci*.

Let *ti* = *ti<sup>u</sup>i* ∈ R<sup>3</sup> (*i* = 1, 2, ... , 6) be the cable tension vector of the *ith* cable, where *ti* = *<sup>t</sup>i* ∈ R<sup>3</sup> is the value of the cable tension. According to the equilibrium equation, the moment load *τ* applied on the moving platform with respect to the point *O* satisfies

$$
\mathbf{r} + \sum\_{i=1}^{6} \mathbf{a}\_{i} \times \mathbf{t}\_{i} = \mathbf{0}.\tag{13}
$$

Denoting *T* = (*<sup>t</sup>*1, *t*2, ... , *<sup>t</sup>*6)*<sup>T</sup>* ∈ R<sup>6</sup> as the vector of six cable tension values, the moment load *τ* can be written as

$$\mathbf{r} = -\sum\_{i=1}^{6} \mathbf{a}\_{i} \times \mathbf{t}\_{i} = -\sum\_{i=1}^{6} (\mathbf{a}\_{i} \times \mathbf{u}\_{i}) \mathbf{t}\_{i} = \mathbf{J}^{T} \mathbf{T} \tag{14}$$

where *J* ∈ R6×<sup>3</sup> is called the Jacobian matrix of the 3DJM and yields

$$J = -\left(\mathfrak{a}\_1 \times \mathfrak{u}\_1, \mathfrak{a}\_2 \times \mathfrak{u}\_2, \dots, \mathfrak{a}\_6 \times \mathfrak{u}\_6\right)^T. \tag{15}$$

Denoting *S* = *JT* ∈ R3×<sup>6</sup> and *Sji* (*j* = 1, 2, 3; *i* = 1, 2, ... , 6) as the components of *S*, the components of *τ* yield

$$
\pi\_j = \sum\_{i=1}^{6} S\_{ji} t\_i. \tag{16}
$$

Then, the expression *Lk* ◦ *τj* in (10) can be written as

$$L\_k \circ \tau\_{\dot{\gamma}} = \lim\_{\Delta t \to 0} \frac{\tau\_{\dot{\gamma}}|\_{\mathcal{R}(t+\Delta t)} - \tau\_{\dot{\gamma}}|\_{\mathcal{R}(t)}}{\omega^k \Delta t} = \frac{\partial \tau\_{\dot{\gamma}}}{\partial \zeta^k} \Big|\_{\mathcal{R}(t)} = \sum\_{i=1}^6 \left( \frac{\partial \mathcal{S}\_{\dot{\gamma}i}}{\partial \zeta^k} t\_i + S\_{\dot{\gamma}i} \frac{\partial t\_i}{\partial \zeta^k} \right). \tag{17}$$

where *<sup>τ</sup>jR*(*t*)(*j* = 1, 2, 3) are the components of *τ* when the 3DJM stays at the pose *<sup>R</sup>*(*t*).

Writing the six cable lengths as a vector *c* = (*<sup>c</sup>*1, *c*2, ... , *<sup>c</sup>*6)*<sup>T</sup>* ∈ R6, then we have the following equation according to the principle of virtual work, i.e.,

$$T^T \delta \mathfrak{c} = \mathfrak{x}^T \delta \mathfrak{y}.\tag{18}$$

Substituting (14) into (18), we have

$$
\delta \mathfrak{c} = f \delta \mathfrak{y}.\tag{19}
$$

It shows that *∂ci ∂ζ<sup>k</sup>* = *Jik* (*i* = 1, 2, ... , 6; *k* = 1, 2, <sup>3</sup>), where *Jik* ∈ R is the component of *J* ∈ R6×3. Writing *ki* = *∂ti ∂ci* as the stiffness of the *ith* cable, the expression *∂ti ∂ζ<sup>k</sup>* in (17) satisfies

$$\frac{\partial t\_i}{\partial \zeta^k} = \frac{\partial t\_i}{\partial c\_i} \frac{\partial c\_i}{\partial \zeta^k} = k\_i I\_{ik} \,. \tag{20}$$

Substituting (16), (17) and (20) into the stiffness model (10), *Kjk* can be written as

$$K\_{jk} = \sum\_{i=1}^{6} \left( \frac{\partial S\_{ji}}{\partial \zeta^k} t\_i + k\_i S\_{ji} J\_{ik} \right) - \frac{1}{2} \sum\_{r=1}^{3} \sum\_{i=1}^{6} \gamma\_{kj}^r S\_{ri} t\_i. \tag{21}$$

Since *S* = *J<sup>T</sup>*, the parametric formulation of the enhanced stiffness model (10) is given below

$$\mathbf{K} = \mathbf{D} + \mathbf{S} \mathbf{K}\_{\text{diag}} \mathbf{S}^T - \widehat{\mathbf{S}T}.\tag{22}$$

where *D* = *∂S ∂ζ*<sup>1</sup> *T*, *∂S ∂ζ*<sup>2</sup> *T*, *∂S ∂ζ*<sup>3</sup> *T* and *<sup>K</sup>*diag = **diag**{*k*1, *k*2, ··· , *k*6}.

#### **3. Modeling of the Cable Tension Distribution Problem Oriented the Enhanced Stiffness of the 3DJM**

For a given pose, the Jacobian matrix of the 3DJM is constant. The stiffness of the 3DJM only depends on the cable tensions according to (22). That means the stiffness of the 3DJM can be adjusted by regulating the cable tensions. Given a desired feasible stiffness matrix *K*des ∈ R3×<sup>3</sup> for the 3DJM at a given pose *R* with a load *τ*, the corresponding cable tensions should be solved from the following equation

$$D + \mathcal{SK}\_{\text{diag}} \mathcal{S}^T - \widehat{\mathcal{S}T} = \mathcal{K}\_{\text{des}}.\tag{23}$$

Due to the complexity of the stiffness formulation of the 3DJM, it is difficult to solve the cable tensions *T* ∈ R<sup>6</sup> from (23) analytically.

As the stiffness *K*des is a 3 × 3 real symmetric matrix in the frame {A}, there exists a real orthogonal matrix *Q* ∈ R3×<sup>3</sup> to transfer the stiffness matrix from the frame {A} to another frame {E}, in which the symmetric stiffness matrix can be represented by a diagonal matrix [34], i.e.,

$$\mathbf{Q}^T \mathbf{K}\_{\rm des} \mathbf{Q} = {}^E \mathbf{K}\_{\rm des} = \mathbf{diag}\{\mathbf{K}\_{d1}, \mathbf{K}\_{d2}, \mathbf{K}\_{d3}\},\tag{24}$$

where the orthogonal matrix *Q* ∈ *SO*(3) satisfies

$$\mathbf{Q}^T \mathbf{Q} = \mathbf{Q} \mathbf{Q}^T = \mathbf{I}\_{3 \times 3}. \tag{25}$$

and *<sup>E</sup>K*des represents the desired stiffness in the frame {E}. Here, *Kdi* (*i* = 1, 2, 3) is the diagonal component of the diagonal matrix *<sup>E</sup>K*des.

Since the desired stiffness is a diagonal matrix in the frame {E}, we can discuss the SCTD problem in the frame {E} to make it simple. On the other hand, as the actual stiffness matrix *K*act may not be exactly symmetric, a symmetric stiffness matrix *K*act can be obtained as below

$$\mathbf{K}\_{\rm act} = \frac{1}{2} (\mathbf{K}\_{\rm act} + \mathbf{K}\_{\rm act}^T). \tag{26}$$

In the frame {E}, the symmetric actual stiffness *K*act is represented by *<sup>E</sup>K*act, which yields

$$\mathbf{^E K}\_{\text{act}} = \mathbf{Q^T} \overline{\mathbf{K}}\_{\text{act}} \mathbf{Q} = \mathbf{Q^T} (\mathbf{D} + \mathbf{S} \mathbf{K}\_{\text{diag}} \mathbf{S}^T - \widehat{\mathbf{S}} \overline{\mathbf{T}}) \mathbf{Q}. \tag{27}$$

Denote *Kai* (*i* = 1, 2, 3) as the diagonal elements of *<sup>E</sup>K*act, the SCTD problem can be formulated as an optimization model

$$\begin{aligned} \text{Minimize } \qquad \mathcal{g}(T) &= \sqrt{\sum\_{i=1}^{3} (K\_{di}(T) - K\_{di})^2}, \\ \text{Subject to } \qquad \qquad ST &= \tau, \\ \qquad T\_{\text{min}} &\preceq T \preceq T\_{\text{max}}. \end{aligned} \tag{28}$$

where *Tmin* ∈ R<sup>6</sup> represents the minimum of the tension vector *T* to avoid the cable be slack and *Tmax* ∈ R<sup>6</sup> represents the maximum of the tension vector *T* to avoid the cable tensions exceeding the capability of the cable driving units. Here, *X Y* (*<sup>X</sup>*, *Y* ∈ R*n*) represents that each element of *X* is no more than the corresponding element of *Y*.

According to analysis of the rank of the Jacobian matrix of the 3DJM in the Appendix A, we have **rank**(*S*) = 3. Then, the column vectors of *S* ∈ R3×<sup>6</sup> can be divided into two parts, *<sup>S</sup>p*, *Ss* ∈ R3×3. Here *Sp* yields **rank**(*<sup>S</sup>p*) = 3 and it is called a basis of matrix *S*. Correspondingly, the cable tension vector *T* ∈ R<sup>6</sup> can be divided into two parts, *<sup>T</sup>p*, *Ts* ∈ R3. That means the six driving cables can be divided into two groups: one group with three cables for position adjustment by regulating the cable lengths, and another group with the remaining three cables for stiffness adjustment by regulating the cable tensions. In this manner, the position and stiffness of 3DJM can be adjusted simultaneously. According to equilibrium Equation (14), the tensions for the position adjustment cables, i.e., *<sup>T</sup>p*, can be represented by the tensions for the stiffness adjustment cables, i.e., *T<sup>s</sup>*, as follows

$$T\_p = \mathbf{S}\_p^{-1} \boldsymbol{\pi} - \mathbf{S}\_p^{-1} \mathbf{S}\_s T\_s. \tag{29}$$

Substituting *Tp* into the optimization model (28), the objective function *g*(*T*) is simplified to *h*(*<sup>T</sup>s*), i.e.,

$$\log(T) = \lg(T\_p, T\_s) = h(T\_s),\tag{30}$$

and the linear inequality constraints are written as

$$
\Phi \mathbf{T}\_8 \preceq w,\tag{31}
$$

where

$$\Phi = \begin{pmatrix} I\_{3 \times 3} \\ -I\_{3 \times 3} \\ -\mathcal{S}\_p^{-1} \mathcal{S}\_s \\ \mathcal{S}\_p^{-1} \mathcal{S}\_s \end{pmatrix} \tag{32}$$

and

$$\boldsymbol{\omega} = \begin{pmatrix} T\_{s-\max} \\ -T\_{s-\min} \\ -S\_p^{-1}\boldsymbol{\tau} + T\_{p-\max} \\ S\_p^{-1}\boldsymbol{\tau} - T\_{p-\min} \end{pmatrix} \tag{33}$$

Here, *<sup>T</sup>p*−*min*, *T<sup>s</sup>*−*min* ∈ R<sup>3</sup> represent the minimum of the tension vector *<sup>T</sup>p*, *T<sup>s</sup>*, respectively, and *<sup>T</sup>p*−*max*, *T<sup>s</sup>*−*max* ∈ R<sup>3</sup> represent the maximum of the tension vector *<sup>T</sup>p*, *T<sup>s</sup>*, respectively. Then, the optimization model (28) is simplified to an equivalent model, where six decision variables are reduced to three and the linear equality constraint are eliminated, i.e.,

$$\begin{array}{ll}\text{Minimize} & h(T\_s) = \sqrt{\sum\_{i=1}^{3} (K\_{ii}(T\_s) - K\_{di})^2}, \\\\ \text{Subject to} & \Phi T\_s - w \preceq 0. \end{array} \tag{34}$$

#### **4. Cable Tension Solution Based on the Genetic Algorithm for the 3DJM**

Since the objective function of the optimization model (34) is complicated and the stiffness of the VSD relative to the cable tension does not need to be different at each point, the widely used gradient-based algorithms for nonlinear optimization problems are not suitable for this optimization model. The direct optimization methods, which do not need the gradient of the objective function, can be employed for this optimization model, such as Complex method, Nelder–Mead algorithm [35] and genetic algorithm. In this paper, a generic method based on the genetic algorithm is proposed to solve the nonlinear constrained optimization model.

A genetic algorithm is inspired by biological evolutionary theory. It is an iterative procedure which usually operates on a population of constant size [36]. In order to apply the genetic algorithm, we revise the optimization model (34) as follows

$$\begin{array}{ll}\text{Maximize} & f(T\_s) = \frac{1}{h(T\_s)} '\\\text{Subject to} & \Phi T\_s - \mathfrak{w} \preceq \mathbf{0}.\end{array} \tag{35}$$

where *Ts* is taken as an individual (also called a *chromosomes*), and *f*(*<sup>T</sup>s*) is taken as the fitness function of the individual *T<sup>s</sup>*.

The genetic algorithm is a stochastic iterative algorithm, where each iteration step is also called a *generation*. Since the genetic algorithm cannot guarantee convergence [36], the termination condition of the proposed method is commonly triggered by finding an acceptable solution for the problem or by reaching a maximum number of generations. Here, we define a parameter *η* ∈ R to evaluate the closeness of two matrices *K*1, *K*2 ∈ R3×3,

$$\eta(\mathbb{K}\_1, \mathbb{K}\_2) = \frac{\mathbb{1}\_1 \cdot \mathbb{1}\_2}{||\mathbb{1}\_1|| \, ||\mathbb{1}\_2||} \times 100\% \tag{36}$$

where *νi* (*i* = 1, 2) is the vector form of *Ki* and *<sup>ν</sup>i* is the norm of *νi*. The iterative procedure of the proposed method is shown below, and it terminates when the following condition achieves

$$1 - \eta(\overline{\mathbf{K}}\_{\text{act}}, \mathbf{K}\_{\text{des}}) \le 0.0005. \tag{37}$$

**Step 1:** *Generate an initial population*

> Generate an initial population of individuals randomly or heuristically. Each individual can be represented by a binary string.

**Step 2:** *Evaluate all individuals*

> Compute *Tp* via (29), and evaluate the fitness function *f*(*<sup>T</sup>s*) for the individuals of the current population. If *Tp* does not satisfy *<sup>T</sup>p*−*min Tp <sup>T</sup>p*−*max*, set the fitness value as zero.

**Step3:** *Check termination condition*

> Check if the current population satisfies the termination condition. Stop the iterative procedure if it satisfies, and generate a new population if not.

**Step 4:** *Selective reproduction*

> Select the individuals of the current population (usually with a probability proportional to their relative fitness values) and produce offspring candidates.

**Step 5:** *Crossover and mutation*

> Perform two operators, named crossover and mutation, on the above offspring candidates for producing a new population. Execute **Step 3** for the new population.


The diagram of the proposed method is shown in Figure 3.

**Figure 3.** Diagram of the proposed method for the stiffness-oriented cable tension distribution problem of the 3DJM.

## **5. Simulation Examples**

In order to validate the proposed method, a simulation was carried out on a certain 3DJM. As shown in Figure 4, the dimensions of the 3DJM are given by *lA* = 0.10 [m], *lB* = 0.13 [m], *eA* = *eB* = 0.002 [m] and *hA* = *hB* = 0.08 [m].

**Figure 4.** CAD model of the 3DJM.

The VSDs were designed to fix onto the moving platform of the 3DJM to extend the range of stiffness variation of each cable. The CAD model of the VSD is shown in Figure 5a. According to the diagram of the VSD, as shown in Figure 5b, the length of the cable in the VSD yields

$$l = \sqrt{h^2 + d^2 - 2hd\cos\phi},\tag{38}$$

and the cable tension *tc* applied on the VSD yields

$$t\_c = \frac{lk\_s(\phi\_0 - \phi)}{hd\sin\phi},\tag{39}$$

where *h* ∈ R is the height of the revolute joint, *d* ∈ R is the length of the rigid-link, *ks* ∈ R is the stiffness of the spring, *φ* ∈ R is the angle of the rigid-link, and *φ*0 ∈ R is the initial value of *φ*. Then the stiffness of the VSD, denoted as *k*VSD ∈ R, can be described by

$$k\_{\rm VSD} = \frac{dt\_{\rm c}}{dl} = \frac{dt\_{\rm c}}{d\phi} \bigg/ \frac{dl}{d\phi}.\tag{40}$$

Given the parameters of the designed VSD, i.e., *ks* = 1.20 Nm/rad, *d* = 0.018 m, *h* = 0.03 m, and *φ*0 = 0.53 rad, the stiffness of the VSD can be approximated by a polynomial expression

$$k\_{\rm VSD} \approx 8.005t\_c^2 - 239.4t\_c + 5415. \tag{41}$$

The corresponding stiffness–tension curve is shown in Figure 5c. For this 3DJM, the lower limits of the cable tensions are given by

$$T\_{\min} = (10, 10, 10, 10, 10, 10)^T \, [N]\_\prime$$

and the upper limits of the cable tensions are given by

$$T\_{\text{max}} = \left(400,400,400,400,400,400\right)^T \text{ [N]}.$$

**Figure 5.** *variable-stiffness device* (VSD) for simulation: (**a**) CAD model of the VSD, (**b**) diagram of the VSD, (**c**) stiffness–tension curve of the VSD.

The simulation examples are implemented for two different desired poses *R*des1 and *R*des2 with different loads, where

$$\mathbf{R}\_{\rm des1} = \begin{pmatrix} 0.992 & -0.060 & -0.111 \\ 0.090 & 0.952 & 0.292 \\ 0.088 & -0.300 & 0.950 \end{pmatrix}, \quad \mathbf{R}\_{\rm des2} = \begin{pmatrix} 0.950 & 0.088 & -0.300 \\ -0.111 & 0.992 & -0.060 \\ 0.292 & 0.090 & 0.952 \end{pmatrix}. \tag{42}$$

For each desired pose, three desired stiffness matrices are given. Thus, the simulation examples can be divided into six sub-cases, as listed in Table 1.


**Table 1.** Six sub-cases of the simulation with different desired poses, applied loads and desired stiffness matrices.

Here, we take Case 1a as an example to perform the proposed method. The other five sub-cases are similar to Case 1a. In Case 1a, the rotation matrix *Q* is computed according to (24),

$$Q = \begin{pmatrix} 0.159 & -0.432 & 0.888 \\ -0.727 & 0.557 & 0.401 \\ 0.668 & 0.710 & 0.225 \end{pmatrix} \tag{43}$$

which transforms *K*des1a to a diagonal matrix *<sup>E</sup>K*des1a = **diag**{−975, −2937, −<sup>2756</sup>} [Nm/rad]. According to (15), the transpose of the Jacobian matrix *J*1 of the 3DJM at the given pose *R*des1 is computed

$$\mathbf{S}\_{1} = J\_{1}^{T} = \begin{pmatrix} 0.05 & -0.02 & -0.04 & -0.01 & 0.03 & 0.06 \\ -0.04 & -0.04 & 0 & 0.05 & 0.05 & 0.01 \\ -0.01 & 0.02 & -0.02 & 0.01 & -0.03 & 0.03 \end{pmatrix} . \tag{44}$$

The matrix *S*1 is divided into two parts as below

$$\mathbf{S}\_{p1} = \begin{pmatrix} 0.05 & -0.04 & 0.03 \\ -0.04 & 0 & 0.05 \\ -0.01 & -0.02 & -0.03 \end{pmatrix}, \quad \mathbf{S}\_{s1} = \begin{pmatrix} -0.02 & -0.01 & 0.06 \\ -0.04 & 0.05 & 0.01 \\ 0.02 & 0.01 & 0.03 \end{pmatrix}. \tag{45}$$

The proposed method is implemented to solve the optimization model (34) to find out the cable tension distribution for the desired stiffness *K*des1a. The iteration curve of the proposed method for Case 1a is shown in Figure 6a, which shows that the optimization process achieves the termination condition within 7 generations and the optimal cable tension distribution for the desired stiffness *K*des1a is

$$T\_{\rm opt1a} = \left(20.84, 72.86, 139.73, 83.88, 123.42, 248.62\right)^{\mathrm{T}} \left[N\right].$$

The corresponding actual stiffness *K*act1a and *K*act1a are computed as below

$$\mathbf{K}\_{\rm actla} = \begin{pmatrix} 1842.37 & 339.16 & 787.59 \\ 324.04 & 388.59 & -29.52 \\ 810.48 & -61.30 & 640.20 \end{pmatrix},\tag{46}$$

$$\mathbf{K}\_{\rm actla} = \begin{pmatrix} 1842.37 & 331.60 & 799.03 \\ 331.60 & 388.59 & -45.41 \\ 799.03 & -45.41 & 640.20 \end{pmatrix}.\tag{47}$$

According to (36), in Case 1a, *η*(*<sup>K</sup>*act1a, *<sup>K</sup>*des1a) = 99.9554%. It shows that, when the cable tensions achieve *<sup>T</sup>*opt1a, the actual stiffness approaches the desired stiffness *K*des1a.

The results of the simulation for all the six sub-cases are summarized in Table 2 and the corresponding iteration curves of the proposed method are given in Figure 6. The results show that the proposed method can achieve the desired stiffness accurately and efficiently. It is effective for the SCTD problem of the 3DJM.

**Figure 6.** Iteration curves of the proposed method for the six sub-cases: (**a**) Case 1a, (**b**) Case 1b, (**c**) Case 1c, (**d**) Case 2a, (**e**) Case 2b, (**f**) Case 2c.



## **6. Conclusions**

Inspired by the structure of human arms, a modular *cable-driven human-like robotic arm* (CHRA) was developed for safe human–robot interaction, since it has advantage of flexibility, low moving mass and intrinsic safety. Due to the unilateral driving properties of the cables, the CHRA is redundantly actuated and its stiffness can be adjusted by regulating the cable tensions. The cable tension distribution problem becomes a key element for the stiffness adjustment of the CHRA. Since the trajectory of the *3-DOF joint module* (3DJM) of the CHRA is a curve on Lie group SO(3), the stiffness of the 3DJM was evaluated by the covariant derivative of the load to the displacement on SO(3), called an enhanced stiffness model of the 3DJM. In this paper, we focus on analyzing how the cable tension distribution problem oriented the enhanced stiffness of the 3DJM. Since the enhanced stiffness model of the 3DJM is too complicated for solving the cable tensions from the desired stiffness analytically, the SCTD problem was formulated as a nonlinear optimization problem. By analyzing the rank of the Jacobian matrix and the equilibrium equation of the 3DJM, a variable elimination technique was employed to simplify the optimization model. A method based on the genetic algorithm was proposed to solve the optimization model, which only utilized the objective function values. The results of a comprehensive simulation show that the proposed method can solve the cable tension distribution from the desired stiffness accurately and efficiently.

**Author Contributions:** Conceptualization, K.Y. and G.Y.; supervision, G.Y.; project administration, C.Z. and Y.C.; writing—review and editing, C.C., T.Z. and T.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by the Public Welfare Technology Research Program of Zhejiang Province, China (Grant number: LGF19E050001), National Natural Science Foundation of China (Grant number:51705510), National-Zhejiang Joint Natural Science Foundation of China (Grant number: U1909215), Institute of robotics and intelligent manufacturing innovation, Chinese Academy of Science (Grant number: C2018005), Ningbo Natural Science Foundation (Grant number: 2019A610114), Open Research Project of the State Key Laboratory of Industrial Control Technology, Zhejiang University, China (Grant number: ICT20048) and Public Welfare Technology Research Program of Zhoushan, China (Grant number: 2019C31067).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Analysis of the Rank of the Jacobian Matrix of the 3DJM**

The transpose of the Jacobian matrix *J* ∈ R6×<sup>3</sup> is denoted as *S* ∈ R3×6, which is represented as below

$$S = f^T = -\left(\mathfrak{a}\_1 \times \mathfrak{u}\_1, \mathfrak{a}\_2 \times \mathfrak{u}\_2, \dots, \mathfrak{a}\_6 \times \mathfrak{u}\_6\right)^T. \tag{A1}$$

Assuming the column rank of *S* satisfies

$$\mathbf{rank}(S) < \mathfrak{Z} \tag{A2}$$

the mixed product of any three column vectors of the matrix *S* yields

$$\mathbf{S}\_{l} \cdot (\mathbf{S}\_{\circ} \times \mathbf{S}\_{k}) = \begin{vmatrix} s\_{l1} & s\_{\circlearrowright} & s\_{k1} \\ s\_{l2} & s\_{\circlearrowright} & s\_{k2} \\ s\_{l3} & s\_{\circlearrowright} & s\_{k3} \end{vmatrix} \equiv 0,\tag{A3}$$

where *i*, *j*, *k* = 1, 2, . . . , 6.

For the 3DJM of the CHRA, the mixed product of any three column vectors of the matrix *S* can be written as

$$\begin{aligned} &\mathbf{S}\_{i}\cdot(\mathbf{S}\_{j}\times\mathbf{S}\_{k})\\ &= (\mathbf{a}\_{i}\times\mathbf{u}\_{i})\cdot[(\mathbf{a}\_{j}\times\mathbf{u}\_{j})\times(\mathbf{a}\_{k}\times\mathbf{u}\_{k})] \\ &= \frac{[\mathbf{a}\_{i}\times(\mathbf{b}\_{i}-\mathbf{a}\_{i})]\cdot\{[\mathbf{a}\_{j}\times(\mathbf{b}\_{j}-\mathbf{a}\_{j})]\times[\mathbf{a}\_{k}\times(\mathbf{b}\_{k}-\mathbf{a}\_{k})]\}}{c\_{i}c\_{j}c\_{k}}\\ &= \frac{(\mathbf{a}\_{i}\times\mathbf{b}\_{i})\cdot[(\mathbf{a}\_{j}\times\mathbf{b}\_{j})\times(\mathbf{a}\_{k}\times\mathbf{b}\_{k})]}{c\_{i}c\_{j}c\_{k}}\\ &\equiv 0 \quad (i,j,k=1,2,\cdots,6). \end{aligned} \tag{A4}$$

In the design of the 3DJM, *eA* ≈ 0 and *eB* ≈ 0, we have *a*1 ≈ *a*2, *a*3 ≈ *a*4, *a*5 ≈ *a*6, and *b*1 ≈ *b*6, *b*2 ≈ *b*3, *b*4 ≈ *b*5. Then, the following mixed product can be represented as

$$\begin{aligned} &(\mathbf{a}\_1 \times \mathbf{b}\_1) \cdot [(\mathbf{a}\_2 \times \mathbf{b}\_2) \times (\mathbf{a}\_3 \times \mathbf{b}\_3)] \\ &\approx (\mathbf{a}\_1 \times \mathbf{b}\_1) \cdot [(\mathbf{a}\_1 \times \mathbf{b}\_3) \times (\mathbf{a}\_3 \times \mathbf{b}\_3)] \\ &= (\mathbf{a}\_1 \times \mathbf{b}\_1) \cdot (\mu \mathbf{b}\_3) \\ &= 0, \\ &(\mathbf{a}\_4 \times \mathbf{b}\_4) \cdot [(\mathbf{a}\_2 \times \mathbf{b}\_2) \times (\mathbf{a}\_3 \times \mathbf{b}\_3)] \\ &\approx (\mathbf{a}\_3 \times \mathbf{b}\_5) \cdot [(\mathbf{a}\_1 \times \mathbf{b}\_3) \times (\mathbf{a}\_3 \times \mathbf{b}\_3)] \\ &= (\mathbf{a}\_3 \times \mathbf{b}\_5) \cdot (\mu \mathbf{b}\_3) \\ &\equiv 0, \\ &(\mathbf{a}\_5 \times \mathbf{b}\_5) \cdot [(\mathbf{a}\_2 \times \mathbf{b}\_2) \times (\mathbf{a}\_3 \times \mathbf{b}\_3)] \\ &\approx (\mathbf{a}\_5 \times \mathbf{b}\_5) \cdot [(\mathbf{a}\_1 \times \mathbf{b}\_3) \times (\mathbf{a}\_3 \times \mathbf{b}\_3)] \\ &= (\mathbf{a}\_5 \times \mathbf{b}\_5) \cdot (\mu \mathbf{b}\_3) \\ &= 0, \end{aligned} \tag{A7}$$

where *μ* ∈ R and *μ* = 0.

Equation (A5) shows that vectors *b*3, *a*1 and *b*1 are in the same plane, (A6) shows that vectors *b*3, *a*3 and *b*5 are in the same plane, and (A7) shows that vectors *b*3, *a*5 and *b*5 are in the same plane.

In summary, we have that vectors *a*1, *a*3, *a*5, *b*1, *b*3 and *b*5 are in the same plane, which is impossible for the 3DJM. The assumption (A2) does not hold, hence

$$\mathbf{rank}(S) = \mathfrak{Z},\tag{A8}$$

i.e.,

$$\mathbf{rank}(f) = \mathbf{3}.\tag{A9}$$
