# Direct Kinematic Singularities and Stability Analysis of Sagging Cable-driven Parallel Robots

Sébastien Briot, Jean-Pierre Merlet

► To cite this version:

Sébastien Briot, Jean-Pierre Merlet. Direct Kinematic Singularities and Stability Analysis of Sagging Cable-driven Parallel Robots. IEEE Transactions on Robotics, 2023, 39 (3), pp.2240-2254. hal-04007182

HAL Id: hal-04007182

https://hal.science/hal-04007182

Submitted on 4 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Direct Kinematic Singularities and Stability Analysis of Planar Sagging Cable-driven Parallel Robots

Sébastien Briot<fnref n="1" /> and Jean-Pierre Merlet<fnref n="2" />

## I. INTRODUCTION

This technical report is associated with the paper [1]. It aims at giving the equations for the computation of the geometrico-static model equations of planar CDPRs based on the Irvine's model, as well as the equations for their stability analysis. It ends with a case study.

## II. GEOMETRICO-STATIC MODEL OF PLANAR SAGGING CDPRs, AND STABILITY ANALYSIS.

In this section, we provide the geometrico-static model of planar sagging CDPRs, and we derive their stability conditions.

### A. Formalization of the problem

Let us consider a planar CDPR with sagging cables moving in a plane $(O, y_0, z_0)$ (Fig. 1). The robot is made of $n$ cables connected to a rigid platform. The cable $i$ ($i = 1, \dots, n$) is connected at point $B_i$ on the platform by a passive revolute joint and at point $A_i$ on the base, a fixed point which is the output point of the winch system<fnref n="1" />. The state vector of the cable $i$ can be defined as:

$$ \mathbf{x}_i(p_i) = [y_i(p_i) \quad z_i(p_i)]^T \in \mathbb{R}^2 \quad (1) $$

where $p_i \in [0 \ell_i]$, $p_i$ is the curvilinear abscissa along the cable, $\ell_i$ is the cable length of the undeformed cable, $y_i(p_i)$ is the location of a point of the cable at the abscissa $p_i$ along the $y_0$ axis, and $z_i(p_i)$ is the location of a point of the cable at the abscissa $p_i$ along the $z_0$ axis. Moreover, the vector $\mathbf{x}_i(p_i = 0) = \mathbf{x}_{i0}$ represents the coordinates of point $A_i$ in the world frame. The state vector derivative, with respect to the initial curvilinear abscissa $p_i$, is given by:

$$ \frac{\partial \mathbf{x}_i}{\partial p_i}(p_i) = \begin{bmatrix} \frac{\partial y_i(p_i)}{\partial p_i} \\ \frac{\partial z_i(p_i)}{\partial p_i} \end{bmatrix} = \mathbf{f}_i(\mathbf{x}_i, \mathbf{u}_i, p_i) \quad (2) $$

<figure><img src="image_1.png" /><figcaption>Fig. 1. Schematic representation of a sagging planar CDPR.</figcaption></figure>

In what follows, it will be necessary to replace again the variable $p_i \in [0 \ell_i]$ associated with the cable $i$ by a normalized variable $s \in [0 1]$ defined by $s = p_i/\ell_i$. Therefore we have

$$ \mathbf{x}'_i(s) = \ell_i \frac{\partial \mathbf{x}_i}{\partial p_i} = \ell_i \mathbf{f}_i(\mathbf{x}_i, \mathbf{u}_i, s) \quad (3) $$

The platform state vector is given by $\mathbf{x}_p = [y_p z_p \varphi_p]^T \in \mathbb{R}^3$, where $y_p$ is the platform center of mass (COM, denoted as point $P$) position along $y_0$, $z_p$ is the platform COM position along $z_0$ and $\varphi_p$ is the angle characterizing the rotation of the platform local frame around $\mathbf{x}_0$. The location of points $B_i$ in the platform frame $\mathcal{F}_p$ are given by the constant vector $^p\overrightarrow{PB_i}$. As a result, $^0\overrightarrow{PB_i} = \mathbf{R}_p \overrightarrow{PB_i}$ where

$$ \mathbf{R}_p = \begin{bmatrix} \cos \varphi_p & -\sin \varphi_p \\ \sin \varphi_p & \cos \varphi_p \end{bmatrix} \quad (4) $$

Finally, the geometric loop-closure constraints of the robot can be written under the form $\beta = [\beta_1^T \dots \beta_n^T]^T = \mathbf{0}$ where

$$ \beta_i(\mathbf{x}_p, \mathbf{x}_{if}) = \begin{bmatrix} y_p \\ z_p \end{bmatrix} + \mathbf{R}_p \overrightarrow{PB_i} - \begin{bmatrix} y_{if} \\ z_{if} \end{bmatrix} = \mathbf{0} \quad (5) $$

Let us assume that all cables are extensible and, for reasons of brevity, that they have all the same cross-section area $A$ and that they are made with the same material having a linear density $\rho$ and a Young's modulus $E$ (assuming a Hookean

<footnote>*This work was partially supported by the French ANR project COSSE-ROOTS (ANR-20-CE33-0001).</footnote>

<footnote n="1">¹ S. Briot is with the Laboratoire des Sciences du Numérique (LS2N) at the Centre National de la Recherche Scientifique (CNRS), 44321, Nantes, France. Sebastien.Briot@ls2n.fr</footnote>

<footnote n="2">² J.-P. Merlet is with the Centre de Recherche Inria Sophia Antipolis, 06902 Sophia Antipolis, France. Jean-Pierre.Merlet@inria.fr</footnote>

<footnote n="1">¹For reasons of brevity, we do not consider here and in the other sections any pulley model. We assume that the point $A_i$ is fixed on the base. The following results could be extended to include pulleys.</footnote>

<footnote n="2">²In what follows, an index “p” will denote a vector in the platform frame while an index “0” will be for the world frame.</footnote>
deformation law). Moreover, we consider that the CDPR is subject to a gravity field $\mathbf{g} = -g\mathbf{z}_0$ ($g > 0$), and that a planar conservative (constant) wrench $\mathbf{w}_p = [F_y F_z M_x]^T \in \mathbb{R}^3$ is applied at point P ($\mathbf{w}_p$ includes the gravity force on the platform). Note that, as we are in statics, $\mathbf{w}_p$ can however take any value. The total potential energy of the CDPR is given by the functional:

$$U = \sum_{i=1}^{n} \ell_i \int_{0}^{1} \left( \frac{EA}{2} (\epsilon_{ix} - 1)^2 + g\rho z_i(s) \right) ds - \mathbf{w}_p^T \mathbf{w}_p \quad (6)$$

where

$$\epsilon_{ix} = \|\mathbf{u}_i(s)\| = (\ell_i \|\mathbf{x}'_i(s)\|) \quad (7)$$

is the strain in the cable i. The term $-\mathbf{w}_p^T \mathbf{w}_p$ in the energy U corresponds to the potential energy of the planar wrench $\mathbf{w}_p$. It should be mentioned that, in this integrand,

$$U_{ei} = \ell_i \frac{EA}{2} \int_{0}^{1} (\epsilon_{ix} - 1)^2 ds \quad (8)$$

is the (elastic) deformation energy of cable i,

$$U_{gi} = \ell_i \int_{0}^{1} g \rho z_i(s) ds \quad (9)$$

is the potential energy due to gravity on cable i, and

$$U_w = -\mathbf{w}_p^T \mathbf{w}_p \quad (10)$$

is the potential energy due to the constant wrench applied on the platform.

The static (stable) configuration of the CDPR must be a local minimum of the potential energy U. As a result, the following fixed-time Bolza optimization problem must be solved:

$$\begin{array}{ll} \min_{\mathbf{u}(s)} & U \\ \text{subject to} & \mathbf{x}'_i = \ell_i \mathbf{f}_i(\mathbf{x}_i, \mathbf{u}_i, s) \quad \text{for } i = 1, \dots, n \\ & \mathbf{x}_i(s_0) = \mathbf{x}_{i0} \\ & \beta_i(\mathbf{x}_p, \mathbf{x}_{if}) = \mathbf{0} \end{array} \quad (11)$$

where $\mathbf{u} = [\mathbf{u}_1^T \dots \mathbf{u}_n^T]^T$.

## B. First-order conditions of optimality

For solving the problem (11), let us define the augmented cost function:

$$U' = G(\mathbf{x}_p, \mathbf{x}_{1f}, \dots, \mathbf{x}_{nf}) + \sum_{i=1}^{n} \int_{0}^{1} \left( H_i(\mathbf{x}_i, \mathbf{u}_i, \boldsymbol{\lambda}_i, s) + \boldsymbol{\lambda}_i^T \mathbf{x}'_i \right) ds \quad (12)$$

where

$$G = -\mathbf{w}_p^T \mathbf{w}_p + \beta^T \boldsymbol{\nu} = -\mathbf{w}_p^T \mathbf{w}_p + \sum_{i=1}^{n} \beta_i^T \boldsymbol{\nu}_i \quad (13)$$

with $\boldsymbol{\beta} = [\beta_1^T \dots \beta_n^T]^T \in \mathbb{R}^{2n}$, $\boldsymbol{\nu} = [\nu_1^T \dots \nu_n^T]^T \in \mathbb{R}^{2n}$, and

$$H_i = \ell_i \left( \frac{EA}{2} (\|\mathbf{u}_i\| - 1)^2 + g\rho z_i - \boldsymbol{\lambda}_i^T \mathbf{f}_i \right) \quad (14)$$

Taking the first variation of $U'$ leads to:

$$\begin{aligned} \delta U' = G_{\mathbf{x}_p} \delta \mathbf{x}_p &+ \sum_{i=1}^{n} G_{\mathbf{x}_{if}} \delta \mathbf{x}_{if} + \beta^T \delta \boldsymbol{\nu} \\ &+ \sum_{i=1}^{n} \int_{0}^{1} \left( H_{i\mathbf{x}_i} \delta \mathbf{x}_i + H_{i\mathbf{u}_i} \delta \mathbf{u}_i \right. \\ & \qquad \left. + (H_{i\boldsymbol{\lambda}_i} + \mathbf{x}_i^{\prime T}) \delta \boldsymbol{\lambda}_i + \boldsymbol{\lambda}_i^T \delta \mathbf{x}'_i \right) ds \end{aligned} \quad (15)$$

The last term in the integrand can be integrated by parts to obtain:

$$\begin{aligned} \delta U' = G_{\mathbf{x}_p} \delta \mathbf{x}_p &+ \sum_{i=1}^{n} (G_{\mathbf{x}_{if}} + \boldsymbol{\lambda}_{if}^T) \delta \mathbf{x}_{if} + \beta^T \delta \boldsymbol{\nu} \\ &+ \sum_{i=1}^{n} \int_{0}^{1} \left( (H_{i\mathbf{x}_i} - \boldsymbol{\lambda}_{i}^{\prime T}) \delta \mathbf{x}_i + H_{i\mathbf{u}_i} \delta \mathbf{u}_i \right. \\ & \qquad \left. + (H_{i\boldsymbol{\lambda}_i} + \mathbf{x}_i^{\prime T}) \delta \boldsymbol{\lambda}_i \right) ds \end{aligned} \quad (16)$$

The necessary condition for local optimality which involves that $\delta U' = 0$ are now expressed by [2], for $i = 1, \dots, n$:

$$\mathbf{x}'_i = \ell_i \mathbf{f}_i(\mathbf{x}_i, \mathbf{u}_i, s) \quad (17)$$

$$\boldsymbol{\lambda}'_i = H_{i\mathbf{x}_i}^T \quad (18)$$

$$H_{i\mathbf{u}_i} = \mathbf{0} \quad (19)$$

$$\boldsymbol{\lambda}_{if} = -G_{\mathbf{x}_{if}}^T \quad (20)$$

$$G_{\mathbf{x}_p} = \mathbf{0} \quad (21)$$

$$\boldsymbol{\beta} = \mathbf{0} \quad (22)$$

Equations (17) to (19) are the ODEs of the Irvine's model [3], leading to the fact that $\boldsymbol{\lambda}_i = \boldsymbol{\tau}_i$ the tension in the cable i (see also [1]). Equation (22) is the vector of the geometry loop-closure constraint equations for the CDPR. Moreover, from (20) and (21), we get:

$$G_{\mathbf{x}_{if}} = -\boldsymbol{\nu}_i^T \Rightarrow \boldsymbol{\lambda}_{if} = \boldsymbol{\tau}_{B_i} = \boldsymbol{\nu}_i \quad (23)$$

$\boldsymbol{\tau}_{B_i}$ being the tension in the cable i at point $B_i$, and

$$\begin{aligned} G_{\mathbf{x}_p} = -\mathbf{w}_p &+ \sum_{i=1}^{n} \left[ \frac{\partial \boldsymbol{\beta}_i}{\partial \mathbf{x}_p} \right]^T \boldsymbol{\nu}_i \\ &- \mathbf{w}_p + \sum_{i=1}^{n} \left[ \frac{\partial \boldsymbol{\beta}_i}{\partial \mathbf{x}_p} \right]^T \boldsymbol{\tau}_{B_i} = \mathbf{0} \end{aligned} \quad (24)$$

where

$$\frac{\partial \boldsymbol{\beta}_i}{\partial \mathbf{x}_p} = \left[ \mathbf{I}_2 \quad (\mathbf{E}^0 \overrightarrow{\mathbf{P}\mathbf{B}_i}) \right] \in \mathbb{R}^{2 \times 3} \quad (25)$$

in which $\mathbf{I}_2$ is the identity matrix of dimension 2 and

$$\mathbf{E} = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix} \quad (26)$$

Note that the equation (24) is nothing else than the equation of the static equilibrium on the platform.
Summarizing the results, the geometrico-static model of the planar CDPR is governed by the equations:

$$
\begin{bmatrix} \mathbf{x}'_i \\ \boldsymbol{\tau}'_i \end{bmatrix} = \mathbf{h}_i = \ell_i \begin{bmatrix} \epsilon_{ix}(s) \cos \theta_i(s) \\ \epsilon_{ix}(s) \sin \theta_i(s) \\ 0 \\ g\rho \end{bmatrix} \quad \text{for } s \in [0 \ 1] \quad (27)
$$

$$
\text{with } \mathbf{x}_i(0) = \mathbf{x}_{i0}, \quad \boldsymbol{\tau}_i(0) = \boldsymbol{\tau}_{i0}, \quad i = 1, \dots, n
$$

$$
\theta_i(s) = \tan^{-1}(\tau_{i2}(s) / \tau_{i1}(s))
$$

$$
\epsilon_{ix}(s) = \frac{\|\boldsymbol{\tau}_i(s)\|}{EA} + 1 \quad (28)
$$

$$
\text{and } \mathbf{b} = \begin{bmatrix} \beta \\ -\mathbf{w}_p + \sum_{i=1}^{n} \left[ \frac{\partial \beta_i}{\partial \mathbf{x}_p} \right]^T \boldsymbol{\tau}_{B_i} \end{bmatrix} = \mathbf{0} \quad (29)
$$

The vector **b** = **b**(*x*<sub>p</sub>, *x*<sub>f</sub>, *τ*<sub>B</sub>) ∈ ℝ<sup>2n+3</sup>, where **x**<sub>f</sub> = [**x**<sub>1f</sub>,..., **x**<sub>nf</sub>]<sup>T</sup> and **τ**<sub>B</sub> = [**τ**<sub>B1</sub>,..., **τ**<sub>Bn</sub>]<sup>T</sup>, is the vector
of the boundary conditions for the system of Differential Algebraic Equations (DAEs). The part **β** = **0** of **b** represents
the geometric constraints on the locations of the extremities
**B**<sub>i</sub> of the cables. The second part of **b** is the mechanical
equilibrium condition. Obviously, when solving (27), we get
the explicit solution of the Irvine's model for a single cable
provided by:

$$
y_i(s) = y_{i0} + \frac{F_{iy} \ell_i s}{EA} + \frac{|F_{iy}|}{\rho g} \left\{ \sinh^{-1} \left[ \frac{F_{iz} - \rho g \ell_i (1-s)}{F_{iy}} \right] - \sinh^{-1} \left[ \frac{F_{iz} - \rho g \ell_i}{F_{iy}} \right] \right\} \quad (30)
$$

$$
\begin{equation}
\begin{split}
z_i(s) = z_{i0} &+ \frac{F_{iz} \ell_i s}{EA} + \frac{\rho g \ell_i^2 s}{EA} \left(\frac{s}{2} - 1\right) \\
& + \frac{1}{\rho g} \left\{ \sqrt{F_{iy}^2 + [F_{iz} - \rho g \ell_i (1-s)]^2} \right. \\
& \qquad \left. - \sqrt{F_{iy}^2 + [F_{iz} - \rho g \ell_i]^2} \right\}
\end{split}
\tag{31}
\end{equation}
$$

as well as the expressions of the tension $\tau_i(s) =$ $[\tau_{i1}(s) \ \tau_{i2}(s)]^T$ all along the cable:

$$
\tau_{i1}(s) = F_{iy} \tag{32}
$$

$$
\tau_{i2}(s) = F_{iz} - \rho g \ell_i (1 - s) \quad (33)
$$

Based on this explicit relationships and the boundary con-
ditions **b** = **0**, the geometric model can be solved by
using a numerical solver. In such a case, the unknown to be
found would be either the cable lengths are rest ℓᵢ and
the cable tensions τ<sub>Bᵢ</sub> (i = 1, . . . , n) for the inverse geometric
model [4], or the platform pose **x**<sub>p</sub> and the tensions τ<sub>Bᵢ</sub> for the
forward geometric model [5]. Note that these problems may
have both multiple solutions, and finding all their solutions is
a rather complicated task.

C. Second-order conditions of optimality

Here, we study the stability conditions of the planar sagging
CDPR.

1) Legendre-Clebsch conditions: We need to check the positive-definiteness of the block-diagonal matrix Huu whose expression is:

$$
\mathbf{H}_{\text{uu}} = \begin{bmatrix} \mathbf{H}_{\text{u}_1 \text{u}_1} & \cdots & \mathbf{0} \\ \vdots & \ddots & \vdots \\ \mathbf{0} & \cdots & \mathbf{H}_{\text{u}_n \text{u}_n} \end{bmatrix} \quad (34)
$$

where

$$
\mathbf{H}_{\mathbf{u}_i \mathbf{u}_i} = 2 \begin{bmatrix} \frac{\epsilon_{ix}^3 - u_{i2}^2}{\epsilon_{ix}^3} & \frac{u_{i1} u_{i2}}{\epsilon_{ix}^3} \\ \frac{u_{i1} u_{i2}}{\epsilon_{ix}^3} & \frac{\epsilon_{ix}^3 - u_{i1}^2}{\epsilon_{ix}^3} \end{bmatrix} \quad (35)
$$

with $u_{ij}$ the jth component of $\mathbf{u}_i$ ($j=1,2$, $i=1,...,n$).
$n$ eigenvalues of $\mathbf{H}_{uu}$ are equal to 2 while the $n$ others are
equal to:

$$
\sigma_i = 2 \frac{\epsilon_{ix} - 1}{\epsilon_{ix}} \quad \text{for } i = 1, \dots, n \qquad (36)
$$

From Equation (28), $\epsilon_{ix} > 1$, except if the cable tension $\tau(t)$ is null, which is unfeasible in practice on Earth, as long as we are not in the water. Therefore, we always have $\mathbf{H}_{uu} \succ 0$.

2) Jacobi conditions: Let us now compute the matrix necessary for assessing the existence of conjugate points. For verifying it, let us express the variation δ**b**:

$$
\delta \mathbf{b} = \mathbf{0} = \mathbf{b}_{x_p} \delta \mathbf{x}_p + \mathbf{b}_{x(s)} \delta \mathbf{x}(s) + \mathbf{b}_{\tau(s)} \delta \mathbf{\tau}(s) \quad (37)
$$

where $\delta \mathbf{x}(s) = [\delta \mathbf{x}_1(s)^T \dots \delta \mathbf{x}_n(s)^T]^T \in \mathbb{R}^{2n}$,
$\delta \mathbf{\tau}(s) = [\delta \mathbf{\tau}_1(s)^T \dots \delta \mathbf{\tau}_n(s)^T]^T \in \mathbb{R}^{2n}$, $\mathbf{b}_{\mathbf{x}(s)} \in \mathbb{R}^{(2n+3)\times 3}$,
$\mathbf{b}_{\mathbf{x}(s)} \in \mathbb{R}^{(2n+3)\times 2n}$ and $\mathbf{b}_{\mathbf{\tau}(s)} \in \mathbb{R}^{(2n+3)\times 2n}$. The expression
of $\mathbf{b}_{\mathbf{\tau}(s)}$ (also equal to $\mathbf{b}_{\lambda(s)}$) is detailed hereafter, as well
as the expression of $\mathbf{b}_{\mathbf{x}_p}$. The expression of $\mathbf{b}_{\mathbf{x}(s)}$ has no
importance in what follows; therefore its calculation is
skipped.

**Computation of the matrix b<sub>τ(s)</sub>:** In order to compute the matrix **b**<sub>τ(s)</sub>, three different types of functions must be calculated, as shown in [1]: ∂**b**/∂**x**<sub>f</sub>, ∂**b**/∂**τ**<sub>B</sub> and the transition matrix **Φ**(s, s<sub>f</sub> = 1) of the DAE (27), which relates ∂**x**<sub>f</sub> and ∂**τ**<sub>B</sub> to ∂**x**(s) and ∂**τ**(s) such that [2]:

$$
\begin{bmatrix} \delta \mathbf{x}(s) \\ \delta \mathbf{\tau}(s) \end{bmatrix} = \mathbf{\Phi}(s, s_f = 1) \begin{bmatrix} \delta \mathbf{x}_f \\ \delta \mathbf{\tau}_B \end{bmatrix} \quad (38)
$$

where

$$
\Phi(s, s_f = 1) = \begin{bmatrix} \frac{\partial \mathbf{x}(s)}{\partial \mathbf{x}_f} & \frac{\partial \mathbf{x}(s)}{\partial \tau_B} \\ \frac{\partial \tau(s)}{\partial \mathbf{x}_f} & \frac{\partial \tau(s)}{\partial \tau_B} \end{bmatrix} \quad (39)
$$

and **x**<sub>f</sub> = **x**(s = 1) and **τ**<sub>B</sub> = **τ**(s = 1). The expression of the vector **b** being an explicit function of **x**<sub>f</sub> and **τ**<sub>B</sub>, it directly comes that:

$$
\mathbf{b}_{\mathbf{x}_f} = \frac{\partial \mathbf{b}}{\partial \mathbf{x}_f} = \begin{bmatrix} -\mathbf{I}_{2n} \\ \mathbf{0}_{3 \times 2n} \end{bmatrix} \quad (40)
$$

$$
\mathbf{b}_{\tau_B} = \frac{\partial \mathbf{b}}{\partial \tau_B} = - \begin{bmatrix} \mathbf{0}_{2n \times 2n} \\ \left[ \left[ \frac{\partial \beta_1}{\partial x_p} \right]^T \quad \dots \quad \left[ \frac{\partial \beta_n}{\partial x_p} \right]^T \right] \end{bmatrix} \quad (41)
$$

where **I**<sub>2n</sub> is the identity matrix of dimension 2n and the
matrix ∂β<sub>i</sub>/∂**x**<sub>p</sub> is defined in (25).
From [1], [2], the matrix $\Phi(s, s_f = 1)$ is nothing else than:

$$ \Phi(s, s_f) = \begin{bmatrix} \frac{\partial \mathbf{x}}{\partial \mathbf{x}_f} & \frac{\partial \mathbf{x}}{\partial \tau_B} \\ \frac{\partial \tau}{\partial \mathbf{x}_f} & \frac{\partial \tau}{\partial \tau_B} \end{bmatrix} \quad (42) $$

Noticing that, by expressing (30) and (31) at $s = 1$, it comes that $\partial \mathbf{x}_f / \partial \mathbf{x}_0 = \mathbf{I}_{2n}$, then the blocks of the matrix $\Phi(s, s_f = 1)$ can be obtained by differentiating the equations (30)–(33) by $\mathbf{x}_f$ and $\tau_B$:

$$ \frac{\partial \mathbf{x}}{\partial \mathbf{x}_f} = \frac{\partial \mathbf{x}}{\partial \mathbf{x}_0} = \mathbf{I}_{2n}, \quad \frac{\partial \tau}{\partial \mathbf{x}_f} = \mathbf{0}, \quad \frac{\partial \tau}{\partial \tau_B} = \mathbf{I}_{2n} \quad (43) $$

Moreover, $\partial \mathbf{x} / \partial \tau_B$ is a block-diagonal matrix whose expression is detailed in Appendix A.

Finally, from [1], we know that:

$$ \mathbf{b}_{\tau(s)} = \mathbf{b}_{\mathbf{x}_f} (\mathbf{\Phi}^{-1})_{12} + \mathbf{b}_{\tau_B} (\mathbf{\Phi}^{-1})_{22} \quad (44) $$

where the matrices $(\mathbf{\Phi}^{-1})_{12} \in \mathbb{R}^{2n \times 2n}$ and $(\mathbf{\Phi}^{-1})_{22} \in \mathbb{R}^{2n \times 2n}$ are the upper right and lower right blocks of the matrix $\mathbf{\Phi}^{-1}(s, s_f)$, respectively.

**Computation of the matrix $\mathbf{b}_{\mathbf{x}_p}$:** When differentiating the matrix $\mathbf{b}_{\mathbf{x}_p}$ with respect to $\mathbf{x}_p$, its expression takes the form:

$$ \mathbf{b}_{\mathbf{x}_p} = \begin{bmatrix} \frac{\partial \beta_1}{\partial \mathbf{x}_p} \\ \vdots \\ \frac{\partial \beta_n}{\partial \mathbf{x}_p} \\ \left[ \mathbf{0}_{3 \times 2} - \sum_{i=1}^{n} (\overrightarrow{\mathbf{P} \mathbf{B}_i})^T \tau_{B_i} \right] \end{bmatrix} \quad (45) $$

**Computation of the Jacobi conditions:** Finally, let us find the expression of the matrix $\bar{\mathbf{S}}(s)$ which allows checking the presence of conjugate points. For this, let us first define the matrix $\mathbf{Z}$ which spans the nullspace of the matrix $\mathbf{b}_{\mathbf{x}_p}^T$, i.e.

$$ \mathbf{b}_{\mathbf{x}_p}^T \mathbf{Z} = \mathbf{0} \quad (46) $$

Assuming that $\mathbf{b}_{\mathbf{x}_p}$ is of full rank, which is true as long as all points $B_i$ are not all superposed<fnref n="3" />, then $\mathbf{Z} \in \mathbb{R}^{(2n+3) \times 2n}$. Left-multiplying (37) by $\mathbf{Z}^T$, we get

$$ \mathbf{0} = \mathbf{Z}^T \mathbf{b}_{\mathbf{x}(s)} \delta \mathbf{x}(s) + \mathbf{Z}^T \mathbf{b}_{\tau(s)} \delta \tau(s) \quad (47) $$

where $\mathbf{Z}^T \mathbf{b}_{\mathbf{x}(s)} \in \mathbb{R}^{2n \times 2n}$ and $\mathbf{Z}^T \mathbf{b}_{\tau(s)} \in \mathbb{R}^{2n \times 2n}$. We then deduce that

$$ \delta \tau(s) = -(\mathbf{Z}^T \mathbf{b}_{\tau(s)})^{-1} (\mathbf{Z}^T \mathbf{b}_{\mathbf{x}(s)}) \delta \mathbf{x}(s) \quad (48) $$

and that, from [1], [2]

$$ \bar{\mathbf{S}}(s) = -(\mathbf{Z}^T \mathbf{b}_{\tau(s)})^{-1} (\mathbf{Z}^T \mathbf{b}_{\mathbf{x}(s)}) \quad (49) $$

As a result, assuming that $\mathbf{Z}^T \mathbf{b}_{\mathbf{x}(s)}$ is never rank deficient (which is true under mild conditions detailed in [1], [2], $\mathbf{S}^{-1}(s)$ becomes rank deficient if the matrix $\mathbf{Z}^T \mathbf{b}_{\tau(s)}$ is singular on the interval $s \in [0 \ 1]$. The rank deficiency of the matrix $\mathbf{Z}^T \mathbf{b}_{\tau(s)}$ on the interval $s \in [0 \ 1]$ is the condition for the appearance of the conjugate points.

Let us now deal with the general case: the spatial sagging CDPRs.

<figure><img src="image_2.png" /><figcaption>Fig. 2. Schematic representation of the planar three-cable CDPR under study (to scale).</figcaption></figure>

### III. CASE STUDY

In this section, we analyze the stability of a planar three-cable CDPR (Fig. 2). The cable properties are: Young's modulus $E = 210$ GPa, cable cross-section radius $r = 1$ mm, material density $\rho/A = 7800$ kg/m<fnref n="3" />. The platform mass is equal to $m_p = 0.25$ kg. There is no external force except gravity applied on the platform.

The anchor point positions on the base are: $\mathbf{x}_{10} = [-1 \ 0]^T \text{ m}$, $\mathbf{x}_{20} = [1 \ 0]^T \text{ m}$, $\mathbf{x}_{30} = [0 \ 0]^T \text{ m}$. The platform anchor points in the platform frame are parameterized by: $^p\overrightarrow{\mathbf{P}\mathbf{B}_1} = [-0.1 \ 0]^T \text{ m}$, $^p\overrightarrow{\mathbf{P}\mathbf{B}_2} = [0.1 \ 0]^T \text{ m}$, $^p\overrightarrow{\mathbf{P}\mathbf{B}_3} = [0 \ 0.1]^T \text{ m}$.

In order to analyze the stability of this robot, we first compute the end-effector configuration space at a given orientation $\varphi_p = -2.5$ rad. For this, we consider the poses $\mathbf{x}_p^g$ located at the nodes of a regular grid in the $yz$-plane and an initial configuration $S_{\text{init}} = (\ell_{\text{init}}, \mathbf{x}_p^g, \tau_B^g)$. Using a continuation process [6], we are able to calculate the configurations $(\ell^g, \mathbf{x}_p^g, \tau_B^g)$ which are connected to $S_{\text{init}}$ for each pose $\mathbf{x}_p^g$ of the grid.

For every configurations of the computed configuration space (for a platform orientation $\varphi_p = -2.5$ rad), the inverse condition numbers $\kappa_{\ell}^{-1}$ and $\kappa_p^{-1}$ of matrices $[\mathbf{b}_{\ell} \ \mathbf{b}_{\tau_0}]$ and $[\mathbf{b}_{\mathbf{x}_p} \ \mathbf{b}_{\tau_0}]$, respectively, are calculated. For a better display by color in Fig. 3, instead of showing directly $\kappa_{\ell}^{-1}$ and $\kappa_p^{-1}$, we provide their value at a given power: $\kappa_{\ell}^{-1.25}$ and $\kappa_p^{-0.6}$. Dark blue zones are zones with bad inverse condition numbers. It should be mentioned that these condition numbers are computed using matrices with components which have non-homogeneous units. This is valid because we want to analyze the degeneracy of the studied matrices, and not to characterize the physical performance of the robot [7].

The results show that:

* The matrix $[\mathbf{b}_{\ell} \ \mathbf{b}_{\tau_0}]$ loses its rank, i.e. the inverse condition number $\kappa_{\ell}^{-1}$ drops down to zero, near the end-effector configuration space boundaries (see the examples of configurations corresponding to these different zones in points $Q_4$, $Q_5$, $Q_6$ and $Q_7$ in Figs. 4 and 5). In these configurations, a variation of the cable length cannot bring any motion of the platform. Three main scenarios appear: either (a) some cables are too slack, (b) all cables are totally tensed, and (c) a cable has almost a null length.

<footnote n="3">³In this case, the platform is no more a solid, but is restricted to a point (i.e. the end-effector has only two degrees of freedom). As a result, there is no singularity of $\mathbf{b}_{\mathbf{x}_p}$ at all.</footnote>
<figure><img src="image_3.png" /><figcaption>(a) Inverse condition number of $[\mathbf{b}_{\ell} \quad \mathbf{b}_{\tau_0}]$ at the power $1.25$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Inverse condition number of $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ at the power 0.6</figcaption></figure>

<caption>Fig. 3. Inverse condition numbers for the matrices $[\mathbf{b}_{\ell} \quad \mathbf{b}_{\tau_0}]$ and $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ for the planar three-cable CDPR in its end-effector configuration space, for a platform orientation $\varphi_p = -2.5$ rad.</caption>

* The matrix $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ loses its rank, i.e. the inverse condition number $\kappa_p^{-1}$ drops down to zero, inside the end-effector configuration space (Fig. 3(b)). Singularities of $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ are suspected here.

In Figure 4, we show the robot end-effector configuration space, for the platform orientation $\varphi_p = -2.5$ rad. In this picture, both previously computed condition numbers are used to provide an overview of the kinemato-static model singularities in the end-effector configuration space: In particular, we highlight in red the areas where the inverse condition number of $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ is lower than $5 \cdot 10^{-4}$, i.e. near which singularities of $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ are foreseen to be present. These red zones separate the end-effector configuration space into two connected components. For assessing the presence of singularities in the red zones, we define a path between points $Q_1 = (0.78, -0.14)$ m and $Q_3 = (0.48, -0.22)$ m (Fig. 4) along which the stability criterion defined previously

<figure><img src="image_5.png" /><figcaption>Fig. 4. Stable and unstable end-effector configuration spaces of the planar three-cable CDPR, for a platform orientation $\varphi_p = -2.5$ rad. In green, the areas where the inverse condition number of $[\mathbf{b}_{\ell} \quad \mathbf{b}_{\tau_0}]$ is lower than $5 \cdot 10^{-4}$; in red, the areas where the inverse condition number of $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$ is lower than $5 \cdot 10^{-4}$.</figcaption></figure>

will be computed. Along this path, 100 points are defined. From Figs. 6(a) to 6(f), we observe the following things:

* From point $Q_1$ to point $Q_2$ (excluded) on the singularity loci of matrix $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$, there is no conjugate point

* At point $Q_2$, a conjugate point appears on the integration interval at $s = 0$, meaning that there is a singularity of the matrix $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$.

* From point $Q_2$ to point $Q_3$, a conjugate point exists on the interval $s \in [0 \ 1[$, meaning that the robot is unstable along this path.

Note that, in Fig. 6, the robot platform configurations have been voluntarily zoomed in order to show that it is impossible to understand by a direct visualization of the robot if its configuration is stable or not.

The red areas computed using the continuation algorithm provide only an indication of the presence of singularity curves for $[\mathbf{b}_{x_p} \quad \mathbf{b}_{\tau_0}]$. However, the change in the number of conjugate points actually proves the presence of these curves. Moreover, the existence of these singularity curves and the presence of no or 1 conjugate point allow also to conclude that:

* Zone A is the zone of the stable configurations,

* Zone B is the zone of the unstable configurations.

Because there is no possibility to check the stability of the solutions of the Irvine's model itself except applying the theory presented in this paper, in order to definitely validate these results on the stability prediction, we developed a *static* lumped model of our robot [8]. Each cable is now discretized with *N* elements. The robot potential energy becomes then an analytic function of the finite element variables plus the platform configuration variables, instead of a functional. As a result, checking the stability of the lumped model configuration remains to checking the positive-definiteness of the
<figure><img src="image_6.png" /><figcaption>(a) Configuration at point $Q_4$: a cable is slack. $\kappa_{\ell}^{-1} = 0.0003$ and $\kappa_p^{-1} = 0.0018$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Configuration at point $Q_5$: all cables are tensed. $\kappa_{\ell}^{-1} = 9.28 \cdot 10^{-5}$ and $\kappa_p^{-1} = 3.78 \cdot 10^{-7}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Configuration at point $Q_6$: a cable has almost a zero length. $\kappa_{\ell}^{-1} = 0.0008$ and $\kappa_p^{-1} = 0.001$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) Configuration at point $Q_7$: a cable is slack. $\kappa_{\ell}^{-1} = 0.004$ and $\kappa_p^{-1} = 0.014$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(a) Configuration at point $Q_1$ (stable): $\kappa_{\ell}^{-1} = 0.023$ and $\kappa_p^{-1} = 0.004$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Configuration at point $Q_2$ (singular): $\kappa_{\ell}^{-1} = 0.031$ and $\kappa_p^{-1} = 8.29 \cdot 10^{-6}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Configuration at point $Q_3$ (unstable): $\kappa_{\ell}^{-1} = 0.028$ and $\kappa_p^{-1} = 0.005$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) $\det(Z^T b_{\tau(s)})$ at point $Q_1$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(e) $\det(Z^T b_{\tau(s)})$ at point $Q_2$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(f) $\det(Z^T b_{\tau(s)})$ at point $Q_3$</figcaption></figure>

<caption>Fig. 6. Study of the determinant of the matrix $Z^T b_{\tau(s)}$ all along the path between points $Q_1$ and $Q_3$.</caption>

reduced Hessian of the energy, as it has been done in [9] for assessing the stability of continuum parallel robots, or in [10] for checking the stability of the straight-cable underactuated CDPRs.

First, in order to assert the lumped model validity, we check its configuration prediction accuracy. For doing this, we first compute the length of the cables based on the Irvine's model for each of the 100 configurations along the path between points $Q_1$ and $Q_3$. Then, we reperform the same computation, but with the lumped model. Results for the maximal cable length estimation error are summarized in Tab. I for several values of numbers of elements used per cables in the lumped model. We see that, with more than 30 elements per cable, this error is lower than $1.5 \cdot 10^{-5}$ m (all cables measuring more than dozens of centimeters). This error of prediction between the two models is negligible, and 30 elements per cables are

<table><caption>TABLE I
MAXIMAL ERROR OF MODELLING FOR THE LUMPED CABLE MODEL WITH RESPECT TO THE IRVINE'S MODEL, AS A FUNCTION OF THE NUMBER N OF ELEMENTS USED PER CABLE.</caption><tr><td>Err.</td><td>N = 2</td><td>N = 5</td><td>N = 10</td><td>N = 15</td><td>N = 30</td><td>N = 60</td></tr><tr><td>[mm]</td><td>3.789</td><td>0.601</td><td>0.149</td><td>0.065</td><td>0.015</td><td>0.003</td></tr></table>

used in the remaining of this section.

Then, we compute the smallest eigenvalue $\sigma_1$ of the reduced Hessian of the lumped model's potential energy and we check its positivity all along the path between configurations $Q_1$ and $Q_3$ (Fig. 7). Recall that a negative value for $\sigma_1$ corresponds to an unstable configuration of the model, while a positive value indicates a stable configuration [9]. We compare this value at each robot configuration with the maximal value of
<figure><img src="image_16.png" /><figcaption>Fig. 7. Comparison at each configuration all along the path of the stability prediction based on the Jacobi condition computed on the Irvine's model (computation of max($\det(Z^T b_{\tau(s)})$) on the integration interval $s \in [0, 1[$) with a stability prediction based on a lumped model (computation of the smallest eigenvalue $\sigma_1$ of the Hessian of the potential energy).</figcaption></figure>

$\det(Z^T b_{\tau(s)})$ on the interval $s \in [0, 1[$. We see that, up to the configuration 76 (point $Q_2$), the smallest eigenvalue $\sigma_1$ is positive, indicating the configurations are stable, which is confirmed by the fact there is no conjugate points for the Irvine's model (Fig. 6(d)). In point $Q_2$, $\sigma_1$ is null, which corresponds to the appearance of the conjugate point at $s = 0$ (Fig. 6(e)). Between $Q_2$ and $Q_3$, $\sigma_1$ is negative, showing that the robot modelled with the lumped cables is unstable, as it was also detected thanks to the Jacobi criterion on the Irvine's model (Fig. 6(f)).

Note that the Matlab code for computing the stability of the robot configuration is provided as an external material.

## APPENDIX

### A. Computation of the matrix $\partial\mathbf{x}/\partial\tau_B$ for the planar sagging CDPRs

The matrix $\partial\mathbf{x}/\partial\tau_B$ is a block-diagonal matrix whose $i$th block is given by the expression:

$$ \frac{\partial \mathbf{x}_i}{\partial \tau_{B_i}} = \begin{bmatrix} \phi_{i12}^{(1)} & \phi_{i12}^{(2)} \\ \phi_{i12}^{(3)} & \phi_{i12}^{(4)} \end{bmatrix} \qquad (50) $$

where

$$ \phi_{i12}^{(1)} = \frac{\ell_i s}{EA} + \frac{\mathrm{sign}(F_{iy})}{\rho g} \left( \frac{(F_{iz} - \rho g \ell_i)}{\sqrt{(F_{iz} - \rho g \ell_i)^2 + F_{iy}^2}} - \frac{F_{iz} - \rho g \ell_i (1 - s)}{\sqrt{(F_{iz} - \rho g \ell_i (1 - s))^2 + F_{iy}^2}} \right) + \frac{\mathrm{sign}(F_{iy})}{\rho g} \left( \sinh^{-1} \left[ \frac{F_{iz} - \rho g \ell_i (1 - s)}{F_{iy}} \right] - \sinh^{-1} \left[ \frac{F_{iz} - \rho g \ell_i}{F_{iy}} \right] \right) \quad (51) $$

$$ \phi_{i12}^{(2)} = \frac{\mathrm{sign}(F_{iy})}{\rho g} \left( \frac{1}{\sqrt{(F_{iz} - \rho g \ell_i (1 - s))^2 + F_{iy}^2}} - \frac{1}{\sqrt{(F_{iz} - \rho g \ell_i)^2 + F_{iy}^2}} \right) \quad (52) $$

$$ \phi_{i12}^{(3)} = \frac{F_{iy}}{\rho g} \left( \frac{1}{\sqrt{(F_{iz} - \rho g \ell_i (1 - s))^2 + F_{iy}^2}} - \frac{1}{\sqrt{(F_{iz} - \rho g \ell_i)^2 + F_{iy}^2}} \right) \quad (53) $$

$$ \phi_{i12}^{(4)} = \frac{\ell_i s}{EA} + \frac{F_{iy}}{\rho g} \left( \frac{F_{iz} - \rho g \ell_i (1 - s)}{\sqrt{(F_{iz} - \rho g \ell_i (1 - s))^2 + F_{iy}^2}} - \frac{F_{iz} - \rho g \ell_i}{\sqrt{(F_{iz} - \rho g \ell_i)^2 + F_{iy}^2}} \right) \quad (54) $$

where $\tau_{B_i} = [F_{iy} F_{iz}]^T$ and the function $\mathrm{sign}(x)$ takes the value 1 if $x > 0$, -1 if $x < 0$ and 0 if $x = 0$.

## REFERENCES

[1] S. Briot and J. Merlet, "Direct kinematic singularities and stability analysis of sagging cable-driven parallel robots," *IEEE Transactions on Robotics*, 2022, under review.

[2] D. Hull, *Optimal Control Theory for Applications*, F. Ling, Ed. New York, NY, USA: Oxford University Press: Springer, 2003.

[3] H. Irvine, *Cable Structures*. MIT Press, 1981.

[4] J.-P. Merlet and J. A. di Sandretto, "The forward kinematics of cable-driven parallel robots with sagging cables," in *Mechanisms and Machine Science*. Springer International Publishing, aug 2014, pp. 3–15.

[5] J.-P. Merlet, "On the inverse kinematics of cable-driven parallel robots with up to 6 sagging cables," in *2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2015)*. IEEE, sep 2015.

[6] K. G. Eugene L. Allgower, *Numerical Continuation Methods*. Springer Berlin Heidelberg, Oct. 2011. [Online]. Available: https://www.ebook.de/de/product/19302145/eugene_l_allgower_kurt_georg_numeric Continuation_methods.html

[7] J. Merlet, "Jacobian, manipulability, condition number, and accuracy of parallel robots," *ASME Journal of Mechanical Design*, vol. 128, no. 1, pp. 199–206, 2006.

[8] T. K. Mamidi and S. Bandyopadhyay, "Forward dynamic analyses of cable-driven parallel robots with constant input with applications to their kinetostatic problems," *Mechanism and Machine Theory*, vol. 163, p. 104381, 2021.

[9] S. Briot and A. Goldsztejn, "Singularity conditions for continuum parallel robots," *IEEE Transactions on Robotics*, pp. 1–19, 2021, in Early Access.

[10] M. Carricato and J.-P. Merlet, "Stability analysis of underconstrained cable-driven parallel robots," *IEEE Trans. on Robotics*, vol. 29, no. 1, pp. 288–296, 2013. [Online]. Available: http://www-sop.inria.fr/coprin/PDF/tro_carricato_2013.pdf