<figure><img src="image_1.png" /></figure>

HAL
open science

# Directional Critical Load Index: a Distance-to-Instability Metric for Continuum Robots

Federico Zaccaria, Edoardo Idà, Sébastien Briot

<figure><img src="image_2.png" /></figure>

## To cite this version:

Federico Zaccaria, Edoardo Idà, Sébastien Briot. Directional Critical Load Index: a Distance-to-Instability Metric for Continuum Robots. IEEE Transactions on Robotics, In press. hal-04639378

HAL Id: hal-04639378

https://hal.science/hal-04639378v1

Submitted on 9 Jul 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Directional Critical Load Index: a
Distance-to-Instability Metric for Continuum Robots

Federico Zaccaria<fnref n="1" />,<fnref n="2" />, Edoardo Idà<fnref n="1" />, and Sébastien Briot<fnref n="3" />

***Abstract***—Equilibrium stability assessment is a primary issue in continuum robots (CRs). The possible stable-to-unstable transitions that CRs may admit complicate the use of CRs in tasks where safety and human-robot interactions are mandatory. In this context, metrics measuring the distance from instability are essential but rarely developed. Existing metrics are frequently based on the evaluation of matrices involving mixed units, thus resulting in unit-dependent metrics. Moreover, the physical meaning of existing metric is hard to interpretate. This paper proposes to use the magnitude of a force that brings instability to the CR equilibrium as a measure of the distance to instability. The major advantages of this metric are the intrinsic physical meaning, the practical interpretation of the results, and the well-defined unit of the measurements. The proposed metric (named directional critical load index) is based on the linearization of the eigenvalues of the reduced Hessian matrix of the total potential energy, which can be achieved regardless of the employed discretization technique. Three different case studies illustrate and demonstrate the main results of this paper.

***Index Terms***—Continuum Robots, Performance Index, Equilibrium Stability.

## I. INTRODUCTION

CONTINUUM robots (CRs) are a new class of manipulators developed by researchers to improve the reduced interaction capacities of rigid-link robots [1]. CRs are usually made by the assembly of flexible components, and the motion is obtained by controlled deformations of the constitutive elements of the robot [2]. The intrinsic flexibility of CRs enables the possibility to use them in many fields where human-robot interaction is fundamental, such as minimally invasive surgery [3] and collaborative tasks [4]. Still, also inspection tasks benefit from the CRs ability to reach complex shapes and to work in confined environments [5].

Continuum robotics has been a fast-growing research community in the last decade, and the increasing use of CRs in different technological fields motivated researchers to propose various CRs designs. Serial-like designs [6] include pneumatically-actuated [7], tendon-driven [8], magnetic [9] or concentric tube [10] CRs (CTRs). Alternatively, parallel designs were proposed to increase the load capacity of CRs [11] and to obtain higher accuracy [12]. Parallel designs include continuum parallel robots (CPRs) made by the assembly of several passive slender links arranged in parallel [13], with the motion obtained by the use of traditional motors placed

at the CPRs base, or parallel continuum robots (PCRs) made by several serial CRs mounted in parallel and interactively displacing a shared end-effector (EE) [14].

A well-known limitation of serial and parallel CR designs is their equilibrium stability and consequently its assessment: the elastic structure of CRs enables possible stable-to-unstable transitions, which ultimately depends on how the CR is loaded. For instance, CTRs exhibit instabilities with significant impact on the usability and controllability of the robot [15]. As the tubes rotate and translate with respect to each other, elastic potential energy accumulates until an unstable configuration is met, and the energy is released with a dangerous snapping [16]. To avoid instability, the tube curvature can be optimized [17], but also design anisotropies [18] reduce the risk of snapping. Alternatively, the work of [19] proposes a criterion that, for given tube lengths, ensures the stable behaviour of the CTR no matter the tubes' base orientation is selected. In the same fashion, CPR designs displayed stable to unstable transition [20], which limits the CPR motion abilities [21]. It should be mentioned that stability assessment is a relevant problem not only in robotics: the Euler's buckling load [22] established the foundations for the elastic stability assessment of rigid beam structures. Successive works focused on different aspects of equilibrium stability, such as secondary bifurcations [23], post-buckling instabilities [24], and investigation of stability bifurcations [25].

### A. Related works on Stability Analysis

Stability analysis of mechanical systems made of thin beams is a challenging topic due to the severe nonlinearities resulting from large deformations, and it is an active field of research [26], [27]. Energetic considerations are necessary to characterise equilibrium stability because stable CR configurations are associated with a minimum of the total CR potential energy [28], [29]. When considering the stability of quasi-static systems, two prevalent approaches have been identified in the literature: optimal control approaches [30] and the analysis of the Hessian matrix of the CR energy [29]. Stable-to-unstable transitions are due to the presence of saddle points of the potential energy [29]. In dynamics, other phenomena may appear, like the presence of limit cycles [31], but their investigation is beyond the scope of our work.

When continuous (not discrete) CRs modelling approaches are used (e.g. [13], [32]), equilibrium stability is frequently studied using optimal control approaches (OC) [15]. These approaches derive stability conditions through non-discretized CR equations, and the resulting numerical test, based on the

<footnote n="1">¹DIN, University of Bologna, Bologna, Italy, {federico.zaccaria3,edoardo.ida2}@unibo.it</footnote>

<footnote n="2">²École Centrale de Nantes, Laboratoire des Sciences du Numérique de Nantes (LS2N), Nantes, France</footnote>

<footnote n="3">³CNRS, Laboratoire des Sciences du Numérique de Nantes (LS2N), Nantes, France, Sebastien.Briot@ls2n.fr</footnote>
integration of differential equations, determines the robot's stability. OC approaches provide a rigorous approach minimally affected by discretization issues at the cost of increased mathematical complexity. State-of-the-art methods on CTRs determine the equilibrium stability by the use of OC approaches [15], [30], and OC is also used in [20] to show that CPRs admits stable-to-unstable transitions. Equilibrium stability of planar CPRs also received significant attention, and OC theory is preliminary used in [33], and on a family of three-actuated-DoF planar CPRs in [34]. In [35], OC is used for stability analysis of tendon-driven continuum robots. Indeed, other robotic systems benefit from OC approaches for the equilibrium stability assessment. Sagging cables of cable-driven parallel robots are analogous to long and slender flexible beams frequently modelled with Irvine's model (a particular subcase of the Cosserat beam's model) and OC can be used to assess the equilibrium stability of sagging cable-driven parallel robot [36].

Alternatively, when the *CRs* configuration is described by a finite number of variables (e.g. when using discretized *CRs* equations), equilibrium stability can be characterized by evaluating the positive definiteness of the Hessian matrix of the potential energy [29]. In contrast to *OC* approaches, stability assessment based on the Hessian matrix provides intuitive mathematical derivation and simplicity, but the accuracy depends on the number of discretization coordinates. The approach proposed in [37] established a numerical method (based on the reduced Hessian matrix evaluation) to assess the *CTRs* stability, follow stable paths by continuation, and determine bifurcations where instabilities occur. Equilibrium stability assessment plays a crucial role in the workspace evaluation of *CRs* [38], and unstable configurations define the attainable workspace boundaries [21]. During the workspace computation, the equilibrium stability of each configuration is to be verified, and the positive definiteness of the reduced Hessian matrix is a straightforward and effective strategy that avoids differential equation integrations proper of *OC* approaches. The same approach is used for different problems, such as equilibrium stability assessment of underconstrained cable-driven robots [39] (assuming cables as straight), or stable equilibrium continuation of elastic beams [25].

While characterizing if a configuration is stable or not is rather straightforward by using any of the methods mentioned above, measuring a “distance to instability” is less obvious. However, such stability metrics are essential in many cases, such as when planning stable *CRs* trajectories [40], for defining controllers avoiding instability, or for designing robots with a large stable workspace. However, stability metrics were rarely investigated, and only a few works were proposed in this direction. Based on the *OC* framework, a metric for the equilibrium stability measurement of *CTRs* is given in [30]. Except for the two-tubes case, the proposed metric is based on the use of a matrix determinant which cannot meaningfully indicate the closeness to non-optimality in general [41].<fnref n="1" />

One may use the condition number as an alternative to the matrix determinant to solve this issue. However, when the matrix is made of non-homogeneous units, the condition number depends on the units selection, providing unit-dependent performance measurements. The work of [42] proposed to measure the equilibrium stability of *CTRs* from a different perspective. *S-curves*, that describes how tubes base rotations are transmitted to the *CTR* tip for fixed tubes lengths, are used to describe the input-output relation of a two-tube *CTR*. The observation that stable *CTRs* *S-curve* do not exhibit a negative slope suggests using the *S-curve* slope as stability metric. However, this approach is limited to the two-tubes case. In the case of *CPRs*, a heuristic metric is proposed in [43] within the *OC* framework. The proposed metric is based on an equivalent integration length which indicates the so-called conjugate points where instability should occur. The proposed metric is heuristic, and its use is limited by the fact that there is no straightforward proof that the mechanism reaches the limit of stability when the conjugate points appears at the abscissa zero. Additionally, the sensitivity of this metric to small changes in other model parameters could be high, and the metric should be used with caution [20]. Except for the specific case of the two-tube *CTRs*, and the case where the heuristic metric provided in [43] could be applied, the generic distance-to-instability metrics are based on the computation of determinant or condition number of matrices which may suffer from the same issues. On the one hand, for a given configuration of a robot under study, at a given "physical distance" from the limit of stability, changing the number of variables for parameterizing any discretized model will lead to different values of the matrix determinant or condition number. On the other hand, those matrices mix terms with different units with some terms that are dimensionless, whatever the technique used for modeling the system. Using millimeters instead of meters in such kind of matrices may lead to a drastic change in any metrics based on the computation of the determinant or condition number. Therefore, there is a necessity to provide metrics that are not sensitive to such kind of issues, and that would reflect more the physics rather than numerical issues.

Indeed, according to [44], a good performance metric should i) not involve mixed units of measurement, ii) admit an analytical expression (to be used for optimization), iii) be bounded in magnitude, and iv) attain a physical meaning to enable effective comparison and quantification. However, except for the two-tube *CTR* case in [42], none of the previously discussed metrics ([30], [20]) satisfy these requirements. In this paper, the focus is directed toward a methodology for the measurement of *CRs* distance to instability. A criterion for said measurement is introduced, which is based on the magnitude of the force that brings instability to the *CR* equilibrium; such criterion has a physical meaning and a well-defined unit (Newton) intuitively, and it can be analytically approximated.

## B. Contributions and Outline

The metric proposed in this paper, named the *directional critical load index (DCLI)*, estimates the magnitude of an

<footnote n="1">¹As an example, consider the identity matrix of dimension three, whose determinant is unitary. However, multiplying the matrix by an arbitrarily small constant (e.g. $10^{-5}$) drastically reduces the determinant, but still, the matrix is full rank.</footnote>
external load that will cause a stable-to-unstable transition of the CR when applied even quasi-statically in a given direction. The derivation of the DCLI is based on the energy-based modelling formalism, and its formulation can be applied to CRs of different architectures, serial and parallel alike, as it will be shown in the Case Studies in Sec. IV. The formulation of DCLI does not depend on the selected discretization strategy, and a finite-differences modelling technique [45] is used for its efficiency in computing successive model-equations derivatives. However, any discretized modelling approaches could also be used. A major advantage of the DCLI is its well-defined measurement unit (Newton) and its intuitive lower and critical bound zero: these characteristics lead to practical physical interpretations of the index results. Although it would be beneficial to have an upper bound of DCLI, there may exist no force in a given direction, namely a force of infinite magnitude, to cause the CR instability, and this is further confirmed in Sec. IV.

In this paper, DCLI is defined by considering only external forces as a cause of the CR instability, as they cover many practical cases (e.g. gravitational loads applied to the EE or at a specific CR location, contacts forces). Moreover, the derivation of DCLI assumes to know the application point and direction of the force that causes instability. Although this may seem a limiting factor, guidelines for selecting the application point and direction of the force may come from the practical scope of the CR. For instance, when the CR is used for manipulation tasks, a possible source of instability may be an additional EE load aligned with the gravity. Finally, the DCLI is computed by investigating the influence of an external force on the equilibrium stability when the motors are fixed, considering a practical case where the robot is moved by imposing motor values. In a straightforward manner, the DCLI should be found solving iteratively a nonlinear system of equations, whose solution is time-consuming to compute and which even may admit no solution. Therefore, in practice, we propose to compute an approximate solution derived via a linearization of the eigenvalues of the reduced Hessian matrix of the potential energy. This solution, even if less accurate that the one that may be found by the iterative procedure, is fast to compute. Accordingly, an analytically computable approximation of the external load for which instability occurs, i.e. for which an eigenvalue vanishes, can be established. As an additional benefit, the index derivation is based on algebraic computations only, and differential equations integration typical of OC approaches are avoided [30], [43]. On the other hand, as DCLI is found by a first-order Taylor approximation of the true critical load, its value is only approximated. The approximation, though is reduced when near an instability: thus, the value of the DCLI can be practically exploited in performance-driven design iterations, where a more rapid metric may be favorable over a more accurate one.

The paper is structured as follows. Section II recalls the energy-based modelling approach. Section III is devoted to the stability metric derivation, and case studies are proposed in Section IV to verify the capability of DCLI to measure the distance to the instability and to quantify the closeness of DCLI to the exact critical load. Finally, conclusions and

<figure><img src="image_4.png" /><figcaption>Fig. 1: Continuous parametrization of a flexible beam.</figcaption></figure>

limitations are highlighted in Section V.

## II. MODELLING

This section describes the energy-based modelling approach of this paper. Sec. II-A presents the CRs architecture and relevant variables. Then, Sec. II-B,II-C derive the deformation energy and external load energy, respectively. The total CR energy and constraints are discussed in Sec. II-D. Finally, the discretization process is introduced to derive the CR geometrico-static model in Sec. II-E.

### A. Description of the Continuum Robot

According to [3], a CR is an actuable architecture whose constitutive material forms curves with continuous tangent vectors. In this work, a CR made by $n$ flexible beams is considered. For instance, $n = 1$ may represent an isolated beam, while the case $n > 1$ includes CTRs [16] or CPRs [11]. The CR is actuated by $m$ motors, and the variables associated with the motor actions are grouped into the vector $q_a \in \mathbb{R}^m$. Then, a frame $\mathcal{F}_p$ is rigidly attached to a specific CR location, assumed to be the EE reference point. The EE pose is described by $q_p = [\mathbf{p}_p; \alpha_p] \in \mathbb{R}^{n_c}$, where $n_c = 3$ for the planar case, $n_c \ge 6$ for the spatial case, $\mathbf{p}_p$ describes the EE position, and $\alpha_p$ is a set of orientation parameters that define the rotation matrix $\mathbf{R}_p$.

Assuming the same number of controlled and actuated variables, the vector $q_c \in \mathbb{R}^m$ stacks the controlled variables. Typically, $q_c$ is a subset of $q_p$, and the vector $q_u$ collects the remaining variables of $q_p$ not included in $q_c$.

### B. Deformation Energy

Let us consider a flexible beam as represented in Fig. 1: the beam is assumed to be one of the $n$ flexible components of the CR, and the index $i = 1, \dots, n$ denotes the $i$-th beam. Being $L_i$ the length of the $i$-th beam, the coordinate $s_i \in [0, L_i]$ parameterizes the beam's centerline. A frame $\mathcal{F}_i(s_i)$ is attached at each beam's cross-section and the pose of the cross-section is defined by $\mathbf{T}_i(s_i) \in SE(3)$:

$$ \mathbf{T}_i(s_i) = \begin{bmatrix} \mathbf{R}_i(s_i) & \mathbf{p}_i(s_i) \\ \mathbf{0} & 1 \end{bmatrix} \quad (1) $$

where $\mathbf{R}_i(s_i) \in SO(3)$, $\mathbf{p}_i(s_i) \in \mathbb{R}^3$ represent the rotation matrix and the position of $\mathcal{F}_i(s_i)$ with respect to (w.r.t.) a
fixed frame $\mathcal{F}_0$. To get the deformation energy of the beam, let us first evaluate the strain vector $\xi_i \in \mathbb{R}^6$ at $s_i$:

$$ \hat{\xi}_i(s_i) = \mathbf{T}_i^{-1}(s_i)\mathbf{T}_i'(s_i) \quad (2) $$

with $(\cdot)' = \frac{d}{ds}, \xi_i = [\mathbf{u}_i, \mathbf{v}_i]$ and $\hat{\xi}_i \in se(3)$ is defined as:

$$ \hat{\xi}_i(s_i) = \begin{bmatrix} \hat{\mathbf{u}}_i(s_i) & \mathbf{v}_i(s_i) \\ \mathbf{0} & 0 \end{bmatrix} \quad (3) $$

where $\hat{\mathbf{u}}_i \in so(3)$ is the skew-symmetric matrix obtained by $\mathbf{u}_i$. Moreover, $\mathbf{u}_i \in \mathbb{R}^3$ represents bending and torsion of the beam, and $\mathbf{v}_i \in \mathbb{R}^3$ describes shear and extensibility. Assuming linear isotropic elasticity, the deformation energy of the beam is given by [46]:

$$ V_{ei} = \frac{1}{2} \int_{0}^{L_i} (\xi_i(s_i) - \xi_i^*(s_i))^T \mathbf{K}_i (\xi_i(s_i) - \xi_i^*(s_i)) ds \quad (4) $$

where $(\cdot)^*$ denotes the undeformed strain configuration, and $\mathbf{K}_i \in \mathbb{R}^{6 \times 6}$ is the material stiffness matrix. Frequently, and if the links cross-section has circular symmetry, $\mathbf{K}_i = \text{diag}(EI_x, EI_y, GI_z, GA, GA, EA)$ where $E$ is the Young's modulus, $G$ is the shear modulus, $I_x, I_y, I_z$ are the principal inertia moments of the cross-section, and $A$ is the cross-section area.

### C. External Loads Energy

Let us compute the potential energy due to external conservative loads independent of the robot system and variables$^2$. First, the potential energy due to a distributed force $\mathbf{f}_{di}$ applied at the $i$-th beam is given by:

$$ V_{di} = - \int_{0}^{L_i} \mathbf{f}_{di}(s_i)^T \mathbf{p}_i(s_i) ds \quad (5) $$

Then, concentrated loads$^3$ are considered. The contribution of a concentrated force $\mathbf{f}$ applied at $\mathbf{p}_{app}$ is:

$$ V_f = -\mathbf{f}^T \mathbf{p}_{app} \quad (6) $$

where $\mathbf{p}_{app}$ is the application point of $\mathbf{f}$. For instance, when the force is applied to the EE, $\mathbf{p}_{app} = \mathbf{p}_p$ while, if $\mathbf{f}$ is applied at the coordinate $s_{app}$ of the $i$-th beam, $\mathbf{p}_{app} = \mathbf{p}_i(s_{app})$.

For the following derivation, it is convenient to represent $\mathbf{f}$ as its magnitude $f$ and its direction $\mathbf{d}_f$ [47]. Thus, Eq. (6) becomes:

$$ V_f = -\mathbf{f}^T \mathbf{p}_{app} = -f \mathbf{d}_f^T \mathbf{p}_{app} = -fh \quad (7) $$

with $h = \mathbf{d}_f^T \mathbf{p}_{app}$.

<figure><img src="image_5.png" /><figcaption>Fig. 2: Graphical illustration of the geometric-static model variables.</figcaption></figure>

### D. Continuum Robot Energy and Constraints

The CR total potential energy is obtained as the sum of the deformation energy and external load contributions, that is:

$$ V_{tot} = \sum_{i=1}^{n} (V_{ei} + V_{di}) + V_f \quad (8) $$

Geometric constraints frequently have to be considered in CRs. For instance, the closure-loop geometric constraints of CPRs involve leg and rigid-platform variables [29]. In serial CRs, when a set of independent variables $\mathbf{q}_p$ is used to represent the EE position, geometric constraints are introduced as well. Additionally, when using quaternions to represent orientations, unitary constraints must be enforced [25]. Without loss of generality, the constraints are represented by:

$$ \Phi(\mathbf{q}_a, \mathbf{q}_p, \xi_1(s_1), \dots, \xi_n(s_n)) = \mathbf{0} \quad (9) $$

where a vector $\Phi \in \mathbb{R}^{n_\Phi}$ is introduced to stack all the constraints in homogeneous form.

### E. Discretization and Geometrico-Static Modelling

CRs equilibrium configurations are static configurations corresponding to critical points of $V_{tot}(\mathbf{q}_a, \mathbf{q}_p, \xi_1(s_1), \dots, \xi_n(s_n))$ subject to the constraints imposed by $\Phi(\mathbf{q}_a, \mathbf{q}_p, \xi_1(s_1), \dots, \xi_n(s_n)) = \mathbf{0}$ [48]. The use of shooting-like approaches frequently solves continuous models [13], providing computational performances [43] and widespread applicability [49]. On the other side, discretization of the potential energy equations through a finite number of variables brings simple but effective CRs analysis, in particular when dealing with equilibrium stability assessment. Discretization techniques include piecewise constant curvature [50], finite differences [45], piecewise constant strains [51], piecewise linear strains [52], assumed strain modes [53], and the interested reader is addressed to specialized review papers [48]. Moreover, discretized CR equations offer the possibility to verify the equilibrium stability by checking the positive definiteness of the reduced Hessian matrix of the total potential energy. For these reasons, this paper is based on the use on discretization techniques.

After the discretization process, $V_{tot}$ (defined in Eq. (8)) is discretized with a finite set of variables $\mathbf{q}_{ei} \in \mathbb{R}^{N_{fi}}$, with $N_{fi}$ being the number of discretization variables. Then, being

<footnote>$^2$The same reasoning could be applied if specific conservative loads explicitly depending on the robot variables are considered, such as linear elastic springs attached to the frame and to the robot links</footnote>

<footnote>$^3$Three-dimensional pure moments, which are non-conservative, are assumed not to appear.</footnote>
$N_e = \sum_{i=1}^{n} N_{fi}$, the vector $\mathbf{q}_e = [\mathbf{q}_{e1}, \dots, \mathbf{q}_{en}] \in \mathbb{R}^{N_e}$ stacks all the discretization variables. For later convenience, the vectors $\mathbf{q}_d = [\mathbf{q}_u, \mathbf{q}_e]$, and $\mathbf{x} = [\mathbf{q}_d, \mathbf{q}_c]$ are introduced. Figure 2 proposes a summary of the variables to ease the reader comprehension. The CRs configuration is thus defined by $\mathbf{q}_a, \mathbf{x}$: the total potential energy and the constraints become a discrete function of $\mathbf{q}_a, \mathbf{x}$, that is:

$$V_{tot} = V_{tot}(\mathbf{q}_a, \mathbf{x}); \quad \mathbf{\Phi} = \mathbf{\Phi}(\mathbf{q}_a, \mathbf{x}) \quad (10)$$

In the static case and for a fixed value of $\mathbf{q}_a$, feasible CRs configurations are determined by the minimum of the total potential energy $V_{tot}$ [28], [29]. However, constraints $\mathbf{\Phi}(\mathbf{q}_a, \mathbf{x})$ should be enforced, and first-order Lagrange conditions determine critical points of $V_{tot}$ [54]. Assuming $\nabla_x \mathbf{\Phi}$ full rank, $\mathbf{x}$ is associated with a critical point of $V_{tot}$ if there exist a set of Lagrange multipliers $\lambda \in \mathbb{R}^{n_\Phi}$ such as [54]:

$$\left\{ \begin{array}{l} \nabla_x \mathcal{L} = 0 \\ \nabla_\lambda \mathcal{L} = 0 \end{array} \right. \qquad (11)$$

where $\mathcal{L}$ is the Lagrangian function defined as:

$$\mathcal{L} = V_{tot} + \mathbf{\Phi}^T \mathbf{\lambda} \quad (12)$$

Equations (11) form the geometric-static model of a CR, and it is an undetermined set of $n_c + N_e + n_\Phi$ equations in $m + n_c + N_e + n_\Phi$ unknowns. By fixing $m$ variables to desired values, a square system is obtained. In particular, the forward geometric-static (FGSP) consists in finding $\mathbf{x}$ and $\mathbf{\lambda}$ for given $\mathbf{q}_a$ and external loads:

$$\mathbf{F} = \left\{ \begin{array}{l} \nabla_x V_{tot} + \nabla_x \mathbf{\Phi}^T \mathbf{\lambda} = 0 \\ \mathbf{\Phi} = 0 \\ \mathbf{q}_a - \mathbf{q}_a^d = 0 \end{array} \right. \quad (13)$$

with $\mathbf{q}_a^d$ desired motor values. Equation (13) forms a square system of $m + n_c + N_e + n_\Phi$ nonlinear equations in the unknowns $\mathbf{y} = [\mathbf{q}_a, \mathbf{x}, \mathbf{\lambda}]$ that can be solved by using root-finding techniques, such as the Newton method. The Jacobian matrix $\mathbf{J} = \partial \mathbf{F}/\partial \mathbf{y}$, that can be supplied to the numerical solver of Eq.(13) to speed up the computation, is structured as follows:

$$\mathbf{J} = \frac{\partial \mathbf{F}}{\partial \mathbf{y}} = \begin{bmatrix} \mathbf{A}_1 & \mathbf{U}_1 & \mathbf{P}_1 & \mathbf{G}^T \\ \mathbf{A}_2 & \mathbf{U}_2 & \mathbf{P}_2 & \mathbf{0} \\ \mathbf{I}_m & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix} \quad (14)$$

where $\mathbf{I}_m$ is the identity matrix of dimension $m$, and:

* $\mathbf{A}_1 = \nabla_{\mathbf{q}_a} (\nabla_x \mathcal{L})$, $\mathbf{U}_1 = \nabla_{\mathbf{q}_d} (\nabla_x \mathcal{L})$

* $\mathbf{P}_1 = \nabla_{\mathbf{q}_c} (\nabla_x \mathcal{L})$, $\mathbf{G} = (\nabla_\lambda (\nabla_x \mathcal{L}))^T = \nabla_x \mathbf{\Phi}$

* $\mathbf{A}_2 = \nabla_{\mathbf{q}_a} \mathbf{\Phi}$, $\mathbf{U}_2 = \nabla_{\mathbf{q}_d} \mathbf{\Phi}$, $\mathbf{P}_2 = \nabla_{\mathbf{q}_c} \mathbf{\Phi}$

All the terms included in $\mathbf{J}$ can be computed analytically from $\mathbf{F}$, and their expressions depend on the specific discretization technique employed (the expressions for a finite-difference modelling approach are reported in [55]). As shown in [29], matrix $\mathbf{J}$ contains the necessary information to investigate singularity conditions and to assess the CR equilibrium stabil-

ity. In particular, when a Type-2 singularity<fnref n="4" /> is encountered, matrix $\mathbf{J}$ becomes singular.

### III. DIRECTIONAL CRITICAL LOAD INDEX

This section derives the DCLI focus of this paper. First, equilibrium stability conditions are derived in Sec. III-A. Then, Sec. III-B proposes a numerical approach to compute the exact load value that causes instability, and it discusses its limitations. Sec. III-C presents the distance-to-instability index introduced in this paper. Then, the index computation is discussed in Sec. III-D.

#### A. Equilibrium Stability Assessment

Solutions of Eq. (13) are CRs equilibrium configurations that correspond to critical points of $V_{tot}$. Second-order conditions establish sufficient conditions for a critical point to be a local minimizer of $V_{tot}$ [54] and, consequently, a stable equilibrium configuration.

To assess equilibrium stability, the Hessian matrix of $\mathcal{L}$, $\mathbf{H}$, is computed as follows:

$$\mathbf{H} = [\mathbf{U}_1 \quad \mathbf{P}_1] = \frac{\partial^2 \mathcal{L}}{\partial \mathbf{x} \partial \mathbf{x}}; \quad \mathbf{H} \in \mathbb{R}^{(n_c+N_e) \times (n_c+N_e)} \quad (15)$$

Then, the matrix $\mathbf{Z}$ spanning the right nullspace of $\mathbf{G}$, is defined as:

$$\mathbf{G}\mathbf{Z} = \mathbf{0}; \quad \mathbf{Z} \in \mathbb{R}^{(N_e+n_c) \times n_z} \quad (16)$$

with $n_z = N_e+n_c-n_\Phi$. According to second-order Lagrange conditions, a robot equilibrium configuration is stable if the reduced Hessian matrix $\mathbf{H}^r$ is positive definite [54], where $\mathbf{H}^r$ is:

$$\mathbf{H}^r = \mathbf{Z}^T \mathbf{H} \mathbf{Z}; \quad \mathbf{H}^r \in \mathbb{R}^{n_z \times n_z} \quad (17)$$

In practice, equilibrium stability is assessed by verifying that all eigenvalues of $\mathbf{H}^r$ are strictly positive. Being $\sigma_k$ the k-th eigenvalue of $\mathbf{H}^r$, the configuration is stable if the minimum eigenvalue $\sigma_{min}$ is positive:

$$\sigma_{min} = \min_k (\sigma_k) > 0 \quad (18)$$

#### B. Exact Directional Critical Load Computation

As previously mentioned, this paper proposes to measure the distance from instability by investigating the influence of an external load $\mathbf{f}$ on the equilibrium stability of a configuration obtained by the solution of Eq. (13), and thus for fixed motor values<fnref n="5" />. The force direction $\mathbf{d}_f$ and the application point $\mathbf{p}_{app}$

<footnote n="4">⁴As detailed in [29]: (i) Type 1 singularities are singularities of the inverse geometric-static problem (FGSP) and are limits of the end-effector space: a non-null motion of the actuators leads to no motion of the end-effector; (ii) Type 2 singularities are singularities of the forward geometric-static problem (FGSP) and are limits of the active-joint space: a non-null motion of the end-effector leads to no motion of the actuators. Moreover, it has also been proven that they are the zones delimiting the stable configuration domain, i.e. that they are the limits of stability. In such a configuration, the stiffness of the robot is null along at least one direction.</footnote>

<footnote n="5">⁵It is possible also to measure the influence of an external load on the equilibrium stability of configurations obtained by the solution of the inverse geometric-static problem. However, there is little practical interest in the investigation of the case where the values of some controlled variables are assigned: maintaining fixed the value of controlled variables when varying $\mathbf{f}$ in practical cases is not trivial as it requires complex force-estimation techniques [56]. Instead, it is more frequent for practical applications the case of assigned motor values. Thus, only the FGSP solutions are considered in the following.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 3: Illustration of when solving $\sigma_{min} = 0$ with a root-finding technique may not be effective.</figcaption></figure>

are considered as *known and fixed* (see Eq. (7)), and the goal is to compute the force magnitude $f$ that causes instability in a given direction. The selection of $d_f, p_{app}$ is guided by the scope of the robot in order to consider possible scenarios where an external load is applied to the robot (e.g. an EE load during a pick and place operation, or a possible load on the CPR leg during a contact).

As a distance-to-instability metric, it is desirable to compute the exact critical load that causes instability: this could be done by varying $f$ and solving the FGSP until a value of $f$ that causes instability is met, that is, the value of $f$ for which $\sigma_{min} = 0$. For the scope, a numerical approach may be devised, aiming at an efficient computation of $f_{CRIT}$, based on a nonlinear root-finding technique (e.g. the Newton method): a solution to Eq. (18), where $f$ is the unknown and $\sigma_{min} = 0$, is sought. Unfortunately though, this approach requires solving the FGSP in Eq. (13) at each procedure and, although this approach seems to be more efficient, it may be ineffective. Let us consider for simplicity a Newton method for the solution of $\sigma_{min} = 0$. Starting from $f_0$, the force update may select a new value of $f_{new}$ that corresponds to unstable regions where no FGSP solution is numerically reachable (see Fig. 3). As illustrated in [29], passing from positive to negative $\sigma_{min}$, there exists a value of $f$ for which $\sigma_{min} = 0$, that is, $f_{CRIT}$. This value of $f$ is a stable-to-unstable transition that may correspond to a Type-2 singularity [29], also defining a limit of the FGSP solvability for the same assembly mode of the robot under consideration. In fact, if a type-2 singularity is crossed, the correct solution to the FGSP would be different than the one under consideration, leading to completely different stability performances, and ultimately wrongly assessing the value of $f_{CRIT}$ for the considered assembly mode and equilibrium configuration. In order to avoid these issues, approaches that avoid overpassing the singular configuration like [57], [58] could be used.

In the present work, we prefer to use a different numerical approach for the computation of $f_{CRIT}$, less time efficient but simpler to implement, that employs a constant and limited update of $f$. This approach is illustrated in Fig. 4. Starting from $f = 0$, the value of $f$ is gradually incremented with a fixed increment $\delta f$ until instability occurs. First, the FGSP is solved with a given initial guess $y_0$ and $f_0 = 0$ to get

the configuration $y$. Then, $H_r$ is built from the output of the FGSP, and $\sigma_{min}$ is computed. If $\sigma_{min} \neq 0$, $f_0$ is incremented of a user-defined quantity $\delta f$, that is:

$$f_{new} = f_0 + \delta f \qquad (19)$$

After the force update, the algorithm restarts by repeating the FGSP solution, and an initial guess for the robot configuration is required. The previous FGSP solution may be used as an initial guess, but it is convenient to update better $y$ accordingly to $f_{new}$. In this way, at the next iteration, the convergence of the solver is faster since the given initial guess is in accordance with the new value of $f$. For the scope, the tangent vector $t_f$ is defined as:

$$t_f = \frac{\partial y}{\partial f} / \left\| \frac{\partial y}{\partial f} \right\| \qquad (20)$$

where $\frac{\partial y}{\partial f}$ is obtained by using Eq. (30). Then, $t_f$ is used to update $y$ accordingly to $\delta f$:

$$y_{new} = y + t_f \delta f \qquad (21)$$

The values $y_{new}, f_{new}$ are used as $y_0, f_0$ to solve again the FGSP. The algorithm is repeated until a value of $f$ is found so that $\sigma_{min} = 0$, or the maximum number of allowed iterations is reached.

The proposed numerical approach is iterative and, depending on the choice of $\delta f$, several FGSP solutions are required with a consequent increase in the overall computational cost for the identification of $f_{CRIT}$. Moreover, using a constant and small update $\delta f$ poses complexities in the selection of $\delta f$. Since $f_{CRIT}$ is unknown, $\delta f$ should be selected as small to obtain accuracy in the computation of $f_{CRIT}$. However, small values of $\delta f$ considerably increases the computational cost of the algorithm, and a fast (and reliable) estimation of $f_{CRIT}$ becomes desirable. Thus, the next Section proposes a computationally-efficient approximation of $f_{CRIT}$ to be used as distance-to-instability index.

### C. Distance to Instability Index

The scope of this section is to provide an index approximating the exact value of $f_{CRIT}$ with reduced computational effort. To do this, let us consider Eq. (18). The value of $\sigma_k$ for a given $f = f_k^*$, namely $\sigma_k^*(f_k^*)$, can be approximated by first-order Taylor's expansion of $\sigma_k$ around a generic $f$:

$$\sigma_k^*(f_k^*) \approx \sigma_k(f) + S_k(f)(f_k^* - f) \qquad (22)$$

where $S_k$ can be obtained as [59]:

$$S_k(f) = \frac{d\sigma_k}{df} = \nu_k(f)^T \frac{d\mathbf{H}^r(f)}{df} \nu_k(f) \qquad (23)$$

the proof of Eq.(23) is given in Appendix B. The term $\nu_k$ is the normalized eigenvector associated with the $k$-th eigenvalue of $\mathbf{H}^r$, namely $\sigma_k$. Please note the dependence of $\sigma_k, S_k, \nu_k$ on $f$: these values are computed at the current robot configurations and $\sigma_k^*$ is estimated by using Eq.(22). To estimate the force magnitude for which instability occurs, let
<figure><img src="image_7.png" /><figcaption>Fig. 4: Schematics of the numerical approach for the computation of the exact critical load.</figcaption></figure>

us compute the value of $f_k^*$ for which $\sigma_k^* = 0$. This is obtained by fixing $\sigma_k^* = 0$ in Eq. (22), and rearranging its terms as:

$$f_k^* = f + \left( -\frac{\sigma_k(f)}{S_k(f)} \right) \quad (24)$$

$f_k^*$ represents a value of $f$ for which, approximately, at least one eigenvalue of $H^r$ is zero, and thus the matrix is not positive definite anymore. The quantity $\sigma_k/S_k$ estimates the additional load magnitude to be applied to cause $\sigma_k = 0$, while $S_k$ provides information on the direction of additional force to be applied: as $\sigma_k$ is always positive for stable equilibrium, $\sigma_k/S_k$ could be either positive or negative. Negative values of $S_k$ represent additional forces oriented as $d_f$, while positive values of $S_k$ are associated with an additional load oriented opposite to $d_f$. By using Eq. (24) for each $k, n_z$ different values of $f_k^*$ are obtained. To measure the distance to instability, the smallest additional magnitude should be used as metric. It seems legitimate to verify only the value of $\sigma_k/S_k$ associated with $\sigma_{min}$ but, depending on the value of $S_k$, any eigenvalue could be the one associated with the smaller value of $f_k^*$, and all the $n_z$ eigenvalues must be checked. Moreover, $S_k$ admits both positive and negative values, but only loads oriented as $d_f$ are of interest. Thus, the directional critical load index (DCLI) is defined as follows:

$$DCLI(f) = \min_{\forall S_k < 0} \left( -\frac{\sigma_k(f)}{S_k(f)} \right) \quad (25)$$

DCLI represents the smallest additional force magnitude of $f$ with respect to nominal external actions, i.e. the robot weight, that causes at least a zero eigenvalue and, consequently, a limit of the stable equilibrium<fnref n="6" />. DCLI has a well-defined unit (Newton), and it can be used to measure and physically understand the distance from the instability, where larger values indicate greater distance. Moreover, it should be stressed that DCLI is directional since the direction $d_f$ is known and fixed, and DCLI represents the additional load applied at $p_{app}$ in a given direction $d_f$ that causes instability.

## D. Index Computation

This Section discusses how to compute the DCLI in practice, as several steps are necessary, and differentiating $H^r$ is

<caption>**Algorithm 1: DCLI computation.**</caption>

1. [y, J] = Solve Geometrico-Static Problem;
2. Extract G, H from J;
3. Z = NullspaceComputation(G);
4. Compute H^r = Z^T H Z;
5. [σ, V] = Eigenvalue decomposition of H^r;
6. **if** Equilibrium is Stable then
7. &nbsp;&nbsp;&nbsp;&nbsp;Compute C and ∂y/∂f = -J⁻¹C;
8. &nbsp;&nbsp;&nbsp;&nbsp;Compute ∂H/∂f, ∂Z/∂f, ∂H^r/∂f;
9. &nbsp;&nbsp;&nbsp;&nbsp;**for** k = 1: n_z **do**
10. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;σ_k = σ(k), ν_k = V(:, k);
11. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute S_k = ν_k^T ∂H^r/∂f ν_k;
12. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;f^*(k) = σ_k/S_k
13. &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end**
14. **DCLI** = min_k(f^*)
15. **else**
16. &nbsp;&nbsp;&nbsp;&nbsp;**DCLI** = 0;
17. **end**

not straightforward. The required steps to calculate the DCLI are schematically summarized in Alg. 1, and the detailed methodology for computing its terms is shown in the following. For a given external load, the solution of Eq. (13) gives the CR configuration. However, the resulting configuration depends on the value of $f$ that is, $y = y(f)$. Additionally, after a solution to Eq. (13) is found, matrix J is obtained from Eq. (14), and G and H, are extracted as blocks of J without further computations. Then, Z can be computed from G, and $H^r$ is obtained from Eq. (17). It is noteworthy that even though there is an infinite possibility of computing Z, a specific one is required for properly computing DCLI, as it will detailed at the end of this Section. To check the positive definiteness of $H^r$, an eigenvalue decomposition is performed to get i) the vector $\sigma$ collecting $n_z$ eigenvalues, and ii) the matrix V whose columns are $n_z$ eigenvector. Suppose the equilibrium is stable (verified by Eq. (18)): in that case, the index computation continues. If the configuration is unstable, the DCLI is not defined.

Then, Eq. (23) requires computing $dH^r/df$: a finite-difference approximation may be used as a straightforward solution (see [54], Chapter 8, Section 1), but an analytical formulation for $dH^r/df$ can also be derived when using discretized robot model equations, as shown in the following.

<footnote n="6">⁶In the case no $S_k < 0$ exists, DCLI is set as equal to $\infty$.</footnote>
Employing a finite-difference approximation for $d\mathbf{H}^r/df$ is simple but time-consuming, and, depending on the selected finite-difference approximation strategy, multiple FGSP solutions are required at the cost of higher computational time. Instead, an analytical formulation of $d\mathbf{H}^r/df$ is preferred when DCLI has to be computed several times, such as for workspace characterization. To obtain an analytical formulation of $d\mathbf{H}^r/df$, the use of the product derivative rule on Eq. (17) results in:

$$ \frac{d\mathbf{H}^r}{df} = \frac{d\mathbf{Z}^T}{df}\mathbf{H}\mathbf{Z} + \mathbf{Z}^T\frac{d\mathbf{H}}{df}\mathbf{Z} + \mathbf{Z}^T\mathbf{H}\frac{d\mathbf{Z}}{df} \quad (26) $$

The differentiation of $\mathbf{H}$ is addressed here first, and the computation of $d\mathbf{Z}/df$ is addressed later. Matrix $\mathbf{H}$ is computed after the solution of Eq. (13) and, in general, $\mathbf{H}$ depends on $y(f)$ and $f$, that is:

$$ \mathbf{H} = \mathbf{H}(y(f), f) \quad (27) $$

Consequently, the total derivative of $\mathbf{H}$ w.r.t. $f$ is obtained as the sum of two terms:

$$ \frac{d\mathbf{H}}{df} = \frac{\partial \mathbf{H}}{\partial f} + \sum_{i=1}^{n+n_c+m+n_{\phi}} \frac{\partial \mathbf{H}}{\partial y_i} \frac{\partial y_i}{\partial f} \quad (28) $$

Since $\mathcal{L} = V_{tot} + \mathbf{\Phi}^T\mathbf{\lambda}$ and $V_f$ only explicitly depends on $f$ (see Eq. (7)), the first term of Eq. (28) simplifies as follows:

$$ \frac{\partial \mathbf{H}}{\partial f} = \frac{\partial}{\partial f} \left( \frac{\partial^2 \mathcal{L}}{\partial \mathbf{x} \partial \mathbf{x}} \right) = - \left( \frac{\partial^2 h}{\partial \mathbf{x} \partial \mathbf{x}} \right) \quad (29) $$

where $h = \mathbf{d}_f^T \mathbf{p}_{app}$ is previously defined in Eq. (7). Instead, the term $\partial \mathbf{H}/\partial y_i$ can be computed analytically, and its expression depends on the specific discretization technique employed. Due to its lengthy expression, its formulation in the case of the finite-difference modelling approach is not reported in this paper for brevity, but the interested reader in addressed to a detailed technical report [60].

To compute the second term of Eq. (28), it is necessary to evaluate $\partial y/\partial f$, and the implicit functions theorem is used for the scope. Let us consider Eq. (13): $\mathbf{F}$ is a set of equations in the unknowns $y$ and dependent on the parameter $f$. $\mathbf{F}$ is assumed to be a set of continuous and differentiable functions w.r.t. $y, f$. Given a pair $(y, f)$ that satisfies $\mathbf{F}(y, f) = \mathbf{0}$, and assuming $\partial \mathbf{F}/\partial y$ full rank, there exists a unique function $\mathbf{r} : \mathbb{R} \to \mathbb{R}^{n+n_c+m+n_{\phi}}$ such as $y = \mathbf{r}(f)$. Moreover, the partial derivative of $\mathbf{r}$ w.r.t. $f$ (and thus $\partial y/\partial f$) is given by:

$$ \frac{\partial \mathbf{r}}{\partial f} = \frac{\partial y}{\partial f} = - \left( \frac{\partial \mathbf{F}}{\partial y} \right)^{-1} \frac{\partial \mathbf{F}}{\partial f} = - \mathbf{J}^{-1} \mathbf{C} \quad (30) $$

where $\mathbf{J}$ is defined in Eq. (14), and $\mathbf{C} = \partial \mathbf{F}/\partial f = [-\nabla_{\mathbf{x}}h; \mathbf{0}]$.

Let us now consider $\mathbf{Z}$: a specific computation methodology of $\mathbf{Z}$ is proposed, which allows for a streamlined derivation of $\partial \mathbf{Z}/\partial f$, that is ultimately needed in Eq. (26). Nullspace bases are frequently computed using numerical techniques, such as singular value decomposition. The resulting nullspace basis is orthonormal, that is, $\mathbf{Z}^T\mathbf{Z} = \mathbf{I}$ and $\mathbf{Z}$ is dense. However, since $\mathbf{Z}$ is obtained numerically, this approach leads to a cumbersome derivation of $\partial \mathbf{Z}/\partial f$, which depends on the specific numerical algorithm employed, and which may not

ensure continuity of $\mathbf{Z}$ with respect to variations of $f$ [61]. On the other side, fundamental nullspace basis [62] offers a non-orthonormal alternative that can be computed analytically. Thus, the analytical formulation of $\mathbf{Z}$ enables the possibility to calculate $\partial \mathbf{Z}/\partial f$ more easily.

To get a fundamental basis of $\mathbf{Z}$, it is necessary to permute the columns of $\mathbf{G}$ to obtain $\mathbf{G}_p = \mathbf{G}\mathbf{P}$, where $\mathbf{P} \in \mathbb{R}^{(m+n_c) \times (m+n_c)}$ is a permutation matrix. The scope of $\mathbf{P}$ is to get a full rank and well-conditioned matrix $\mathbf{G}_d \in \mathbb{R}^{n_{\phi} \times n_{\phi}}$ from:

$$ \mathbf{G}_p = [\mathbf{G}_d \quad \mathbf{G}_u] \quad (31) $$

with $\mathbf{G}_u \in \mathbb{R}^{n_{\phi} \times n_z}$. Matrix $\mathbf{P}$ can be found by inspecting $\mathbf{G}$ and identifying a set of $n_{\phi}$ linearly independent columns that form a well-conditioned $\mathbf{G}_d$. However, when the dimension of $\mathbf{G}$ increases, no trivial full-rank partition is available in general. As stated in [63], any choice of $\mathbf{P}$ ensuring $\mathbf{G}_d$ full-rank is adequate and, by randomly scanning the columns of $\mathbf{G}$ until a well-conditioned partition is found, matrix $\mathbf{P}$ can be computed. However, the computational cost of this strategy drastically increases with $n_{\phi}$. For instance, when using finite-difference techniques for the geometrico-static modelling, the quaternion-unitarity constraints must be enforced for each beam's cross-section, leading to a large value of $n_{\phi}$. To overcome this issue, a heuristic approach is proposed in Appendix A for the computation of $\mathbf{P}$. However, other approaches may be equivalently proposed for the scope.

As long as $\mathbf{G}_d$ is full rank, a fundamental nullspace basis of $\mathbf{G}_p$ is obtained as:

$$ \mathbf{Z}_p = \begin{bmatrix} \mathbf{Z}_d \\ \mathbf{I}_{n_z} \end{bmatrix} \quad (32) $$

with $\mathbf{Z}_d = -\mathbf{G}_d^{-1}\mathbf{G}_u$. Finally, $\mathbf{Z}$ is obtained by permutating the columns of $\mathbf{Z}_p$ as done for the rows of $\mathbf{G}_p$, that is $\mathbf{Z} = \mathbf{P}\mathbf{Z}_p$.

To analytically calculate $d\mathbf{Z}/df$, the first step requires to compute $d\mathbf{G}/df$. Since $\mathbf{G}$ depends on $y(f)$ and not on $f$ explicitly, $d\mathbf{G}/df$ is expressed as follows:

$$ \frac{d\mathbf{G}}{df} = \sum_{i=1}^{m+n_c+N_e+n_{\phi}} \frac{\partial \mathbf{G}}{\partial y_i} \frac{\partial y_i}{\partial f} \quad (33) $$

The expression of $\partial \mathbf{G}/\partial y_i$ depends on the specific modelling strategy, and its are reported in [60] for the interested reader. Then the derivative of $\mathbf{G}_d, \mathbf{G}_u$ w.r.t. $f$ are obtained by using $\mathbf{P}$ as follows:

$$ \frac{d\mathbf{G}_p}{df} = \frac{d\mathbf{G}}{df}\mathbf{P} = \begin{bmatrix} \frac{d\mathbf{G}_d}{df} & \frac{d\mathbf{G}_u}{df} \end{bmatrix} \quad (34) $$

The term $d\mathbf{Z}_p/df$ is obtained by deriving Eq. (32) w.r.t. $f$:

$$ \frac{d\mathbf{Z}_p}{df} = - \begin{bmatrix} \frac{d\mathbf{G}_d^{-1}}{df} \mathbf{G}_u + \mathbf{G}_d^{-1} \frac{d\mathbf{G}_u}{df} \\ \mathbf{0} \end{bmatrix} \quad (35) $$

where $d\mathbf{G}_d^{-1}/df = \mathbf{G}_d^{-1}(d\mathbf{G}_d/df)\mathbf{G}_d^{-1}$. Finally, $d\mathbf{Z}/df$ is recovered as $d\mathbf{Z}/df = \mathbf{P}d\mathbf{Z}_p/df$.
<figure><img src="image_8.png" /><figcaption>Fig. 5: Buckling of straight beams: (a) clamped beam with free distal section, and (b) clamped beam with pinned distal section. Then, the values of the *DCLI* (dotted line) and *σ*<sub>min</sub> (continuous line) are displayed by varying *f* for the scenario (a),(b) in (c),(d), respectively. The decreasing trend of *DCLI* (f) when approaching the critical load confirms the correctness of the distance-to-instability measurement.</figcaption></figure>

IV. CASE STUDIES

This Section proposes three different case studies: initially straight beams and their elastic buckling (Sec. IV-A), a two-tubes *CTR* with four controlled *DoF* (Sec. IV-B), and a spatial *CPR* with two controlled *DoF* (Sec. IV-C). These case studies are selected to illustrate how the proposed formulation of the *DCLI* can be applied to different scenarios (passive elements, serial *CRs*, and parallel *CRs*). Finite differences [29] are used as discretization techniques to obtain the geometric-static model of Eq. (13) for each case study. However, any other discretization technique can be used. Even though finite differences do not offer the best performances in terms of computational time [64], the analytical formulation of Eq. (13) considerably simplifies the computation of Equations (28) and (33). The details on the model's implementation and the explicit derivation of the equations are not reported in this section for brevity sake. However, detailed equations are available on the technical report [60] associated with this paper.

A. Buckling of beams

This Section proposes the analysis of initially straight beam instability and the comparison of the *DCLI* with Euler's critical load for the beam buckling to verify the correctness of the equilibrium stability prediction and the *DCLI*. For each beam, *DCLI* has been computed by performing a finite-differences discretization with $N_e = 200$. As shown in [45], a planar beam $N_e \ge 50$ ensures sufficient accuracy in the geometric-static problem solution. The early work of Euler [22] defined analytical conditions for the buckling of ideal elastic beams subjected to axial loads. For each considered beam, $L = 1$ m, the cross-section is circular with diameter 2 mm, and $E = 210$ GPa.

First, let us consider a clamped-free beam as illustrated in Fig. 5a: the beam is clamped at the proximal section, initially straight, and parallel to the fixed-frame *x* axis. Euler's buckling formula provides the value of the compressive axial force magnitude *f* to be applied at the tip of the beam to cause elastic instability. Under the assumptions that shear and extensibility are negligible, the Euler's critical load $f_{EUL}$ for
a clamped-free beam can be computed as:

$$f_{EUL} = \frac{EI\pi^2}{(2L)^2} \qquad (36)$$

for the selected beams parameters, $f_{EUL}$ results in 0.407 N. When no load is applied at the beam, $DCLI$ (0) = 0.407 N. $DCLI$ (0) is in accordance with $f_{EUL}$, and the difference between $DCLI$ (0) and $f_{EUL}$ is negligible up to four digits.

Then, the exact critical load $f_{CRIT}$ is computed by using the numerical algorithm proposed in Sec. III-B. The force applied at the beam tip is gradually increased with a $\delta f = 0.001$ N, and, for each step, $\sigma_{min}$ and $DCLI$ (f) are computed. The results of this computation are reported in Fig. 5c. The exact critical load value results in $f_{CRIT} = 0.407$ N, obtained in 407 steps (and thus 407 geometric-static problem solutions). In particular, $\sigma_{min}$ becomes negative and the equilibrium unstable when the Euler's load is reached, confirming the correctness of $f_{CRIT}$. Moreover, for each value of f, the sum $DCLI(f) + f \approx f_{EUL}$ is constant$^7$: as the beam remains straight and undeformed at each step, the value of $DCLI + f$ truly represents the critical load. Similar results were obtained in [43], where the index based on an equivalent integration length gives the exact beam length for which instability should occur for a given load.

Then, let us consider the beam of Fig. 5b: the beam is clamped at the proximal section, initially straight, parallel to the fixed-frame x axis, and pinned at the distal section. In this case, the Euler's critical load is obtained as:

$$f_{EUL} = \frac{EI\pi^2}{(0.699L)^2} \qquad (37)$$

and, with the selected beams parameters, $f_{EUL} = 3.322$ N and $DCLI = 3.322$ N. The difference between $DCLI$ and $f_{EUL}$ is negligible up to four digits. The exact critical load $f_{CRIT}$ is computed with the numerical algorithm proposed of Sec. III-B: the value of $f_{CRIT} = 3.322$ N is obtained with 3322 steps, with $\delta f = 0.001$ N. To confirm the correctness of $f_{CRIT}$, it is possible to note that $\sigma_{min}$ of Fig. 5d becomes negative when the Euler's load is reached. As before, $DCLI$ tends toward zero when $\sigma_{min}$ decreases, and $DCLI(f) + f$ is constant and equal to $f_{EUL}$.

## B. A two tubes CTR with four controlled DoF

This Section introduces the application of $DCLI$ for a two-tube $CTR$, a well-known class of CRs where instability occurs [19]. The two-tube case is a simple but effective benchmark to test the $DCLI$ since an analytical condition exists for the stability assessment, and is therefore used to showcase that our index performs at least as state-of-the-art ones. A $CTR$ made by $n = 2$ concentric tubes is considered (Fig. 6a). A fixed frame $\mathcal{F}_0$ is attached to the robot base, the $CTR$ centerline is parametrized with the coordinate $s$, and the index $i$ represents the $i$-th tube, where $i = 1$ is the inner tube, and $i = 2$ is the outer tube. Tubes are of length $L_i$ (measured from

$s = 0$) and actuated at $s = -\beta_i$. The tubes are actuated in translation and rotation: $\theta_{i0}$ is the rotation of the tube's base, and $\beta_i$ is called transmission length. The $CTR$ energy is obtained by considering shear-less and inextensible tubes, and the discretization process is performed by using finite differences, with 50 points for $CTR$ sections from 0 to $L_2$ and from $L_2$ to $L_1$.

As previously mentioned, the two-tube $CTR$ is a well-known situation where analytical conditions exist for the global stability of the $CTR$. Assuming planar precurvature only $(u_i^* = [u_{ix}^*, 0, 0])$, and no external load applied to the robot, the $CTR$ equilibrium is globally stable if the following inequality is verified [19]:

$$\zeta_{\gamma} = \frac{\cot(\gamma)}{\sqrt{\gamma}} < \zeta_{lim} \qquad (38)$$

where $\zeta_{\gamma}$ is computed by the knowledge of $\gamma$, defined as follows:

$$\gamma = L_2^2 u_{1x}^* u_{2x}^* \frac{k_{1b} k_{2b} (k_{1t} + k_{2t})}{k_{1t} k_{2t} (k_{1b} + k_{2b})} \qquad (39)$$

The term $k_{bi}$ is the flexural stiffness, and $k_{ti}$ is the torsional stiffness of the $i$-th tube while the term $\zeta_{lim}$ of Eq. (38) is equal to zero in the case $\beta_1$ and $\beta_2$ are assumed to be zero. Equation (38) determines conditions for the global stability of two-tube $CTR$s, and $S$-curves were introduced to practically visualize the $CTR$ motion abilities. $S$-curves describe the relationship between the base orientations and the resulting $CTR$ tip orientation, providing an intuitive representation of the robot motion abilities [19]. First, let us define $\phi_0 = (\theta_{10} - \theta_{20})$ as the rotation offset of the tubes at the base, while $\phi_L = (\theta_1(L_2) - \theta_2(L_2))$ is the rotation offset of the tubes at $L_2$. Then, $S$-curves are built by fixing $\theta_{20} = 0$, repeatedly solving the $FGSP$ with $\phi_0 \in [0, 2\pi]$ and measuring $\phi_L$ for each $\phi_0$.

First, the correctness of equilibrium stability prediction is verified by computing $S$-curves with no load applied to the robot ($f = 0$). The base orientation of the outer tube is fixed as zero, and the inner tube is rotated at the base counterwise by a small step. At each rotation, the $FGSP$ is solved with given base angles (thus fixed $\phi_0$) and the previous solution is used as initial guess. Then, the orientation of both tubes at the $CTR$ distal section is measured to compute $\phi_L$: as a consequence of the torsion over the tubes, these values differ for the orientations at the base.

Tubes parameters are selected as $E = 80$ GPa, $u_{1x}^* = 1/0.50$ m$^{-1}$, $u_{2x}^* = 1/0.70$ m$^{-1}$. The first tube has inner diameter $d_{1inn} = 1.0$ mm, and outer diameter $d_{1out} = 1.5$ mm, while the second tube has $d_{2inn} = 1.5$ mm, and $d_{2out} = 2.0$ mm. Computing the $S$-curves for $L_2 = 0.4$ m ensures globally stable equilibrium since $\zeta_{\gamma} = -0.164 < 0$. This is confirmed by the always positive value of $\sigma_{min}$, as shown in Fig. 7a. Instead, the $S$-curves for $L_2 = 0.6$ m violates the inequality of Eq. (38) ($\zeta_{\gamma} = +0.242 > 0$), and $\sigma_{min}$ goes toward negative values (Fig. 7b). Moreover, the observation that $S$-curves of stable $CTR$s are monotonic (the slope of the curve is always positive) motivated [42] to consider the slope

<footnote>$^7$As in any discretization process, the number of discretization coordinates influences the accuracy of the results. With the selected number of discretization coordinates ($N_c = 200$), the obtained numerical errors are lower that $10^{-4}$, thus negligible for the scope of $DCLI$</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 6: a two-tube CTR. (a) relevant dimensions and variables, (b) a cross-section of the tubes to highlight the torsion angles.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7: S-curves of a two-tube CTR. $\sigma_{min}$ is depicted over a globally stable S-curve, and on an unstable case (b). Instead, for the same unstable case, (c) displays the value of the stability index of [42].</figcaption></figure>

of S-curves as a stability metric. In particular, the stability metric based on the slope SL is defined as:

$$SL = \frac{\pi}{2} - \operatorname{atan}\left(\frac{\partial\phi_L}{\partial\phi_0}\right) \quad (40)$$

The index SL has a well-defined unit (radians), SL > 0, as long as the equilibrium is stable and larger values of SL indicate greater distance from instability. Fig. 7c illustrates the values of SL for the case L₂ = 0.6 m, and these values will be later used to verify the correctness of DCLI w.r.t. state-of-the-art indices, such as SL<fnref n="8" />

To compute DCLI, the influence of a tip load on the equilibrium stability is considered. The agreement of DCLI with the state-of-the-art is checked by selecting the unstable case L₂ = 0.6 m, and computing the S-curves highlighting DCLI over different directions d_f and no external load applied on the CTR. Figure 8 illustrates six S-curves at the value of DCLI (0): for each curve, d_f is aligned to one of the fixed-frame axes F₀, with positive and negative directions.

By looking at Fig. 8, it is possible to note that DCLI approaches zero when instability is reached<fnref n="9" />. This trend is in accordance with the values of SL (Fig. 7c). This comparison confirms the coherence of DCLI w.r.t. state-of-the-art results.

Then, the estimation provided by DCLI is compared with the exact value of f_CRIT by computing the absolute error e_abs = |f_CRIT - DCLI|. For instance, let us consider d_f = [1, 0, 0]e, and the corresponding values of DCLI over the S-curve of Fig. 8a. For each point of the S-curve, f_CRIT is computed with δf = +0.001 N, and a maximum number of steps equal to 2000, thus exploring forces up to 2 N. The resulting value of e_abs is displayed in Fig. 9a, and two scenarios can be identified. In a first scenario, there exist no f_CRIT between [0,2] N: this happens between φ₀ ∈ [0, 140]°, and in the upper branch of the S-curve. As an example, Fig. 9b reports the trend of σ_min when using the numerical algorithm of Sec. III-B for φ₀ = 60°: the value of σ_min varies when changing f but it never goes to negative values.

<footnote n="8">⁸Unstable solutions are not plotted in Fig. 7c, as well as in the following ones, as specialized numerical methods are needed to correctly track them, and this is not the aim of the current paper.</footnote>

<footnote n="9">⁹It should be considered that DCLI approaches zero in the proximity of a singularity and close to the instability. As DCLI is obtained by a linearization of the reduced Hessian matrix eigenvalues, small values of DCLI are obtained and this means that a small perturbation is required to cause instability.</footnote>
<figure><img src="image_11.png" /><figcaption>Fig. 8: For the same S-curve, different values of DCLI (0) are obtained by changing d_f.</figcaption></figure>

As DCLI is obtained by linearization, it displays finite but high values. This means that the configurations is significantly far from the instability. Instead, between $\phi_0 \in [140, 190]^{\circ}$, there exist values of $f_{CRIT} \in [0, 2]$ N: in this region, when increasing $f$, the value of $\sigma_{min}$ reaches negative values, as illustrated in Fig. 9c for $\phi_0 = 160^{\circ}$. Thus, when there exists a load $f_{CRIT}$ that causes instability, DCLI displays reduced values and, when approaching the stability limits on the S-curve, the difference between DCLI and $f_{CRIT}$ reduces: this is confirmed by the value of $e_{abs}$ which decreases when approaching the instability.

Instead, when DCLI is low, it does not necessarily mean that $f_{CRIT}$ is low. This is evident in the upper part of Fig. 9a where no $f_{CRIT} \in [0, 2]$ N exists, but DCLI is low in the proximity of the instability. However, it should be considered that the S-curve of Fig. 8 (and thus the values of DCLI) are computed with $f = 0$, and the instability phenomenon is happening as a cause of the torsional energy accumulated in the CTR, as described in [15]. When varying $f$ during the computation of $f_{CRIT}$, the shape of the S-curves varies consequently, and it may happen that previously unreachable values of $\phi_0$ becomes accessible: the application of the external load may increase the stability of some configurations.

Then, the influence of the tip force direction $d_f$ on the values of DCLI is investigated. For the scope, $d_f$ is parameterized

by using spherical coordinates as follows:

$$d_f = \begin{bmatrix} \sin \alpha \cos \beta \\ \sin \alpha \sin \beta \\ \cos \alpha \end{bmatrix} \qquad (41)$$

where $\alpha$ is the inclination angle, and $\beta$ is the azimuth. For a given $\phi_0$ and given CTR configuration, a two-dimensional grid that discretizes uniformly $\alpha \in [0, 2\pi], \beta \in [0, 2\pi]$ is generated, and DCLI computed for each pair of $\alpha, \beta$.

First, let us consider the case of $\phi_0 = 60^{\circ}$ illustrated in Fig. 10a. By changing $\alpha, \beta$, the value of DCLI (0) modifies accordingly, and a minimum of DCLI = 0.149 N is found at $\alpha = 1.84, \beta = 5.17$ rad, corresponding to $d_f = [0.47, -0.86, -0.27]$. At the minimum, DCLI (0) and $f_{CRIT}$ are comparable, with 0.214 N and 0.149 N, respectively, and the absolute error is reduced. Instead, the previously investigated direction $d_f = [+1, 0, 0] \to \alpha = \pi/2, \beta = 2\pi$ is far from the minimum, high values of DCLI are displayed, and no $f_{CRIT}$ exists. When considering a configuration closer to the instability ($\phi_0 = 160^{\circ}$, Fig. 10b), the direction where DCLI (0) is minimum modifies to $\alpha = 1.663$ rad, $\beta = 6.098$ rad and DCLI (0) at the minimum is 0.136 N. $f_{CRIT}$ is 0.078 N for the same direction, displaying a reduced absolute error between $f_{CRIT}$ and DCLI (0). It should also be noted that, since polar coordinates are employed, the trend of DCLI (0) is periodic, and the same minimum is found twice in Fig. 10.
<figure><img src="image_12.png" /><figcaption>Fig. 9: Comparison between $f_{CRIT}$ and DCLI (0). Figure (a) displays the value of $e_{abs}$ for the S-curve computed with $L_2 = 0.6$ and $\mathbf{d}_f = +[1, 0, 0]$. In grey zones, no $f_{CRIT}$ exists between $[0, 2]$ N. Figure (b),(c) display the value of $\sigma_{min}$ when using the numerical algorithm of Sec. III-B, for $\phi_0 = 60^\circ, 160^\circ$, respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10: Influence of the inclination angle $\alpha$ and the azimuth angle $\beta$ of the tip force direction $\mathbf{d}_f$ on DCLI. Figure (a) is relative to $\phi_0 = 60^\circ$, while figure (b) corresponds to $\phi_0 = 160^\circ$. Minimums of DCLI are highlighted by a red circle.</figcaption></figure>

As a summary, DCLI effectively measures the distance to instability since i) DCLI goes to zero when the instability occurs, and ii) DCLI is in accordance with state-of-the-art indices. Moreover, near the instability when $f_{CRIT}$ is low, DCLI provides a good estimation of the instability force. Nevertheless, it should be considered that the DCLI can be used with systems with more than two tubes as well. As long as the discretized model equation for the CTR are derived as in Eq.(13), DCLI can be obtained as explained in Sec. III-C by computing successive derivatives of the model equation to obtain $d\mathbf{H}^r/df$ as in Eq.(26).

### C. A spatial CPR with two controlled DoF

This case study aims to illustrate the application of DCLI on a more complex continuum structure, that is, a CPR. In this context, the DCLI can be used to analyze the CPR limit payload, namely an additional gravitational force due to the additional weight on the EE, that would bring the robot instability even in the case of quasi-static motion. In particular, The two-controlled DoF CPR that was proposed in [65] and illustrated in Fig. 11 is used as a benchmark. This CPR is

considered because its workspace and its equilibrium stability limits have been experimentally validated. The CPR focus of this section is a *RFRF* robot, and it has two rotative motors (*R*) whose axes are coaxial and attached to the proximal section of two flexible chains. A transmission system transfers the motion from the two motors to the actuated coaxial joint *R*. Flexible chains, made by several slender beams arranged in parallel, are connected through a passive revolute joint (*R*) and the robot *EE* is coincident with the passive revolute joint *R*, with all the *R* joints being nominally parallel. With this beam arrangement, flexible link 1 and 3 (see Fig. 11a) are synchronously moved by the same motor, whereas the other actuator rotates link 2. Thus, the CPR has two actuated DoF: its end-effector moves in a vertical plane when no loads other than the ones belonging to the plane are applied on the robot. If out-of-plane loads are applied, this robot may deform in the space, as shown in [65]. In addition, some of its planar unstable configurations can only be predicted with spatial models.

An illustration of the prototype modelling framework is proposed in Fig. 11b: a frame $\mathcal{F}_p$ is attached to the robots *EE*, and variables $\mathbf{q}_p$ described the pose of $\mathcal{F}_p$ w.r.t. the fixed
<figure><img src="image_14.png" /><figcaption>Fig. 11: The *CPR* prototype of [65]. Figure (a) provides details on the prototype components, while Figure (b) illustrates the *CPR* modelling schematics.</figcaption></figure>

frame $\mathcal{F}_0$, As one motor actuated link 1 and 3, while the other motor moves link 2, $\mathbf{q}_a \in \mathbb{R}^2$, and the motors are moved to obtain a controlled $xy$ position of the *EE* (according to $\mathcal{F}_0$ of Fig. 11b) As noticeable in Fig. 11a, intermediate constraints are introduced in within each flexible chain to increase the prototype's stiffness in the direction orthogonal to the motion plane. The geometric-static modelling of these link, as well as the overall robot's modelling, have been discussed in [65] and, in the following, the *EE* mass $m = 0.216$ g is considered as a concentrated *EE* load aligned with the gravitational acceleration $\mathbf{g} = [0, -9.81, 0]$, and the links weight is included as distributed loads. The resulting *EE* force is $f_{EE} = 2.11$ N.

The evaluation of the equilibrium stability of the considered *CPR* is crucial since unstable configurations define the limits of the mobility of the robot. Let us first consider the prototype workspace (*WS*), experimentally validated in [65], and illustrated in Fig. 12a. Region (S) is a stable equilibrium region, where the assigned motor values correspond to attainable static robot configurations, while region (U) corresponds to an unstable region. The outer border of the *WS* is associated with a Type-1 singularity [29], which is a limit of the inverse geometric-static problem solvability. Instead, the inner border of (S) is associated with a leg singularity [29], and after crossing it, the robot equilibrium becomes unstable. In particular, when the equilibrium becomes unstable, the *EE* pose is not more controllable since an uncontrolled out-of-the-plane motion occurs.

This phenomenon is highly undesirable for practical applications of this prototype, and the goal is to measure the distance from the instability by using *DCLI*. To characterize the robots *WS* in terms of *DCLI*, the influence of a tip external load with direction aligned to the gravity is investigated ($\mathbf{d}_f = \pm[0, 1, 0]$ accordingly to the $\mathcal{F}_0$ of Fig. 11b). Results are illustrated in Fig. 12b, 12c, where the values of *DCLI* ($f_{EE}$) are displayed. In particular, the attention is directed to values of *DCLI* $\le 5$ N: as *DCLI* is used to measure the

distance to instability, and since $f_{EE} = 2.11$ N, variations of more than 200% of the *EE* load are not of practical interest, since the resulting *CPR* workspace may considerably differ from the one of Fig. 12a. Values of *DCLI* greater than 5 N are undisplayed (part in yellow of Fig. 12).

First, let us consider $\mathbf{d}_f = -[0, 1, 0]$, corresponding to $f$ aligned with gravity. As displayed in Fig. 12b, approaching the instability in the upper *WS* region causes *DCLI* values which tend to zero. Instead, the other *WS* regions display values greatly higher than $f_{EE}$, indicating a larger distance from the instability. On the other hand, when consider $\mathbf{d}_f = +[0, 1, 0]$ the resulting *DCLI* ($f_{EE}$) is illustrated in Fig. 12c. As approaching the instability at the lower *WS* region, *DCLI* ($f_{EE}$) goes to zero.

Then, the closeness of *DCLI* ($f_{EE}$) to $f_{CRIT}$ is quantified. To do this, at each step of the workspace computation *DCLI* and $f_{CRIT}$ and the absolute error $e_{abs} = |f_{CRIT} - \text{CLI}|$ are computed. In particular, to compare *DCLI* ($f_{EE}$) with $f_{CRIT}$, the computation of $f_{CRIT}$ is performed as follows. First, let us consider the workspace of Fig. 12a obtained with $f_{EE} = 2.11$ N aligned with $-y$. At each stable workspace point, the motor angles $\mathbf{q}_a$ are extracted from the *CPR* configuration and consider these values as fixed. Then, the numerical procedure of Sec. III-B is used: the *FGSP* is repeatedly solved with desired motor angles $\mathbf{q}_a$ by increasing the tip load until instability is met or the maximum allowed iterations number is reached. Since $f_{CRIT}$ is to be computed at each workspace point, $|\delta f|$ is selected as 0.01 N as a trade-off between accuracy and computational cost, and the maximum number of iterations is set as 500 to explore forces up to 5 N.

The results of this computation are displayed in Fig. 13a,13b $\mathbf{d}_f = -[0, 1, 0]$, $\mathbf{d}_f = +[0, 1, 0]$, respectively. As for the two-tubes *CTR* case, some regions exist where no $f_{CRIT} \in [0, 5]$ N exists. However, it is important to note that when *DCLI* is low, $f_{CRIT}$ is in accordance with *DCLI*, and $e_{abs}$ is reduced. In particular, regions close to the instability where *DCLI* ($f_{EE}$)
<figure><img src="image_15.png" /><figcaption>Fig. 12: The workspace of the CPR prototype. Figure (a) highlights the stable region S and the unstable region U. On the same workspace, the trend of DCLI ($f_{EE}$) with $d_f = [0, -1, 0]$ is reported in b) for a negative force sense, and c) for a positive force sense. Type-1 singularities are depicted in red, and leg singularities in black.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13: Comparison between $f_{CRIT}$ and DCLI ($f_{EE}$) for $d_f = [0, 1, 0]$. Figure (a) displays the value of $e_{abs}$ for a negative force sense, while figure (b) for a positive force sense. Grey zones are region where no $f_{CRIT}$ exists between [0, 5] N.</figcaption></figure>

$\le 1$ N, display reduced absolute errors $e_{abs}$ ($e_{abs} \le 0.1$ N). This confirms the capability of DCLI to measure the distance to the instability and, close to the instability, estimate the true load to be applied at the CPR to cause unstable transitions.

Finally, we measured the computational performances of DCLI w.r.t. state-of-the-art metrics. Focusing on Fig. 13a, we measured the computational cost for the determination of the determinant of $\mathbf{H}^r$ ($\det(\mathbf{H}^r)$), the conditioning number of $\mathbf{H}^r$, ($\mathrm{cond}(\mathbf{H}^r)$), DCLI and $f_{CRIT}$ at each workspace location of Fig. 13a where $f_{CRIT}$ exits. DCLI requires 0.0690 s, which is higher than the cost of $\det(\mathbf{H}^r)$ (0.0020 s) and $\mathrm{cond}(\mathbf{H}^r)$ (0.1340 s) but, still, its cost is reduced (lower than the average cost of the FGSP solution that is 0.118 s). However, DCLI provides a physical interpretability that $\det(\mathbf{H}^r)$ and $\mathrm{cond}(\mathbf{H}^r)$ cannot guarantee. Instead, $f_{CRIT}$ predicts accurate results but its cost is approximately 100 times higher than DCLI (6.11 s).

V. CONCLUSIONS

This paper proposed a criterion to measure the distance-to-instability of CRs. In contrast to state-of-the-art approaches, the DCLI does not involve the use of mixed units, and it provides the physical meaning of the results. As DCLI represents the additional load to be applied at a defined location of the robot to cause instability, DCLI possesses a well-defined unit (Newton). The applicability of DCLI was demonstrated over different case studies, namely the buckling of straight beams, instability of CTRs, and the stable-to-unstable transitions of an existing CPR prototype. No matter whether the case study is considered, as the instability is approaching, DCLI tends toward zero: this confirms the correctness of the distance-to-instability measurement. To further verify this, DCLI was compared with a state-of-the-art index for CTRs and, as DCLI goes to zero near instability, also the index of [40] vanishes. Finally, the accuracy of DCLI on the critical force estimation, and its computation times, are assessed by comparing DCLI with the exact force that causes instability, the latter computed
<table><caption>**Algorithm 2: Permutation Matrix Computation**</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Function [P] = GetPermutation(G):</td></tr><tr><td>2</td><td>[a, b] = size G;</td></tr><tr><td>3</td><td>P = 0, G<sub>d</sub> = ;</td></tr><tr><td>4</td><td>for i = 1: a do</td></tr><tr><td>5</td><td>R<sub>v</sub> = 0;</td></tr><tr><td>6</td><td>for j = 1: b do</td></tr><tr><td>7</td><td>G<sub>m</sub> = [G<sub>d</sub>(1: i, :), G(1: i, j)];</td></tr><tr><td>8</td><td>R<sub>v</sub>(j) = inverse conditioning of G<sub>m</sub>;</td></tr><tr><td>9</td><td>end</td></tr><tr><td>10</td><td>idx = column of G where max of R<sub>v</sub> occurs;</td></tr><tr><td>11</td><td>G<sub>d</sub> = [G<sub>d</sub>, G(:, idx)];</td></tr><tr><td>12</td><td>Set P(i, idx) = 1;</td></tr><tr><td>13</td><td>Set b = b - 1;</td></tr><tr><td>14</td><td>Remove G(:, idx) from G;</td></tr><tr><td>15</td><td>end</td></tr><tr><td>16</td><td>return</td></tr></tbody></table>

by an *ad hoc* numerical approach proposed in this paper. It is shown that (i) as the instability is approaching, the distance between *DCLI* and the exact load reduces and *DCLI* can be used to estimate the critical instability load, and (ii) *DCLI* is 100 times faster to compute. These characteristics can be positively exploited when iterative computations are required, such as when dealing with optimal design and non-linear optimization: a more rapid metric may be favorable over a more accurate one.

This work establishes the foundations for the future development of analysis tools for *CRs*. Authors believe that even though *DCLI* may only be used to measure the distance-to-instability in quasi-static scenarios, it may still be useful in contexts such as stable trajectory planning of *CRs* or pick-and-place task planning. Also, the applicability of *DCLI* in other underactuated robotics systems, and the extension of the index formulations to include inequality constraints (e.g. for cable-driven parallel robots) will be objective of future development. Nonetheless, the definition of *DCLI* also represents a starting point for more robust and precise indexes, for example accounting for how dynamic forces affect instability, or independent of the force application direction.

# APPENDIX

## A. Heuristic Algorithm for the computation of P

In the following, a heuristic algorithm (Alg. 2) is proposed for the computation of $\mathbf{P}$, aiming to create $\mathbf{P}$ that maximises the inverse conditioning of $\mathbf{G}_d$. The proposed approach requires scanning $n_\phi$ times the columns of $\mathbf{\Lambda}$, aiming to determine a matrix $\mathbf{P}$ that maximizes the inverse conditioning of $\mathbf{\Lambda}_d$. A pseudocode of the algorithm is reported in (Alg. 2). First, the algorithm starts by initializing $\mathbf{P} = \mathbf{0}$ and $\mathbf{\Lambda}_d = \emptyset$. The goal is to select $n_\phi$ columns of $\mathbf{\Lambda}$ to create a full-rank and well-conditioned $\mathbf{\Lambda}_d$. The algorithm starts by scanning the first row of $\mathbf{\Lambda}$ to select the term with the higher inverse conditioning, and the corresponding column of $\mathbf{\Lambda}$ is selected. This column (labelled with *idx*) is inserted in $\mathbf{\Lambda}_d$ and removed from $\mathbf{\Lambda}$ to avoid repetitions. Matrix $\mathbf{P}$ is updated accordingly to put the column *idx* as the first column of $\mathbf{\Lambda}_d$. Then, the second row is considered. For each column of $\mathbf{\Lambda}$, a $(2 \times 2)$

matrix $\mathbf{\Lambda}_m$ is obtained by collecting $\mathbf{\Lambda}_d$ and the considered column of $\mathbf{\Lambda}$ (see line 12 of Alg. 2). The second column to be put in $\mathbf{\Lambda}_d$ is the one that maximises the inverse conditioning of $\mathbf{\Lambda}_m$. Matrix $\mathbf{P}$ is updated to put the selected column as the second of $\mathbf{\Lambda}_d$. The algorithm proceeds in the same fashion for the next rows by building $\mathbf{\Lambda}_m$, selecting columns that maximise the inverse conditioning of $\mathbf{\Lambda}_m$, and creating $\mathbf{P}$ consequently. The algorithm stops when all the $n_\phi$ rows of $\mathbf{\Lambda}$ have been considered.

## B. Proof of Eq. (23)

This Appendix proposes the proof of Eq. (23). Let us consider $\mathbf{H}^r$ being a symmetric matrix with only real entries, and consider the following eigenproblem:

$$ \mathbf{H}^r \boldsymbol{\nu}_k = \sigma_k \boldsymbol{\nu}_k \qquad (42) $$

where $\sigma_k$ is the eigenvalue, and $\boldsymbol{\nu}_k$ the normalized eigenvector. Under these assumptions:

$$ \boldsymbol{\nu}_k^T \boldsymbol{\nu}_k = 1 \qquad \boldsymbol{\nu}_k^T \frac{d\boldsymbol{\nu}_k}{df} = 0 \qquad (43) $$

By differentiating Eq. (42) and left-multiplying by $\boldsymbol{\nu}_k^T$ we obtain:

$$ \boldsymbol{\nu}_k^T \frac{d\mathbf{H}^r}{df} \boldsymbol{\nu}_k + \boldsymbol{\nu}_k^T \mathbf{H}^r \frac{d\boldsymbol{\nu}_k}{df} = \frac{\sigma_k}{df} \boldsymbol{\nu}_k^T \boldsymbol{\nu}_k + \sigma_k \boldsymbol{\nu}_k^T \frac{d\boldsymbol{\nu}_k}{df} \qquad (44) $$

By considering Eq. (43), Eq. (44) simplifies as follows:

$$ \boldsymbol{\nu}_k^T \frac{d\mathbf{H}^r}{df} \boldsymbol{\nu}_k + \boldsymbol{\nu}_k^T \mathbf{H}^r \frac{d\boldsymbol{\nu}_k}{df} = \frac{\sigma_k}{df} \qquad (45) $$

Then, as long as $\mathbf{H}^r$ is symmetric:

$$ (\mathbf{H}^r \boldsymbol{\nu}_k)^T = (\sigma_k \boldsymbol{\nu}_k)^T \rightarrow \boldsymbol{\nu}_k^T \mathbf{H}^r = \sigma_k \boldsymbol{\nu}_k^T \qquad (46) $$

Thus, by inserting Eq. (46) into Eq. (45) and considering Eq. (43), it is possible to note that the term $\boldsymbol{\nu}_k^T \mathbf{H}^r \frac{d\boldsymbol{\nu}_k}{df}$ vanishes and we obtain the final expression of Eq. (23):

$$ \frac{d\sigma_k}{df} = \boldsymbol{\nu}_k^T \frac{d\mathbf{H}^r}{df} \boldsymbol{\nu}_k \qquad (47) $$

# REFERENCES

[1] G. Robinson and J. B. C. Davies, “Continuum robots-a state of the art,” in *Proceedings 1999 IEEE international conference on robotics and automation (Cat. No. 99CH36288C)*, vol. 4. IEEE, 1999, pp. 2849–2854.

[2] P. K. Singh and C. M. Krishna, “Continuum arm robotic manipulator: A review,” *Universal Journal of Mechanical Engineering*, vol. 2, no. 6, pp. 193–198, 2014.

[3] J. Burgner-Kahrs, D. C. Rucker, and H. Choset, “Continuum robots for medical applications: A survey,” *IEEE Transactions on Robotics*, vol. 31, no. 6, pp. 1261–1280, 2015.

[4] F. Campa, M. Diez, D. Diaz-Caneja, and O. Altzarra, “A 2 dof continuum parallel robot for pick & place collaborative tasks,” in *Advances in Mechanism and Machine Science: Proceedings of the 15th IFToMM World Congress on Mechanism and Machine Science 15*. Springer, 2019, pp. 1979–1988.

[5] D. A. Troncoso, J. A. Robles-Linares, M. Russo, M. A. Elbanna, S. Wild, X. Dong, A. Mohammad, J. Kell, A. D. Norton, and D. Axinte, “A continuum robot for remote applications: From industrial to medical surgery with slender continuum robots,” *IEEE Robotics & Automation Magazine*, 2022.

[6] S. Kolachalama and S. Lakshmanan, “Continuum robots for manipulation applications: A survey,” *Journal of Robotics*, vol. 2020, 2020.
[7] B. A. Jones and I. D. Walker, "Kinematics for multisection continuum robots," *IEEE Transactions on Robotics*, vol. 22, no. 1, pp. 43–55, 2006.

[8] M. Li, R. Kang, S. Geng, and E. Guglielmino, "Design and control of a tendon-driven continuum robot," *Transactions of the Institute of Measurement and Control*, vol. 40, no. 11, pp. 3263–3272, 2018.

[9] Q. Peyron, Q. Boehler, P. Rougeot, P. Roux, B. J. Nelson, N. Andreff, K. Rabenorosoa, and P. Renaud, "Magnetic concentric tube robots: introduction and analysis," *The International Journal of Robotics Research*, vol. 41, no. 4, pp. 418–440, 2022.

[10] P. E. Dupont, J. Lock, B. Itkowitz, and E. Butler, "Design and control of concentric-tube robots," *IEEE Transactions on Robotics*, vol. 26, no. 2, pp. 209–225, 2010.

[11] C. E. Bryson and D. C. Rucker, "Toward parallel continuum manipulators," in *2014 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2014, pp. 778–785.

[12] B. Mauzé, G. J. Laurent, R. Dahmouche, and C. Clévy, "Micrometer positioning accuracy with a planar parallel continuum robot," *Frontiers in Robotics and AI*, vol. 8, p. 706070, 2021.

[13] C. B. Black, J. Till, and D. C. Rucker, "Parallel continuum robots: Modeling, analysis, and actuation-based force sensing," *IEEE Transactions on Robotics*, vol. 34, no. 1, pp. 29–47, 2017.

[14] S. Lilge and J. Burgner-Kahrs, "Kinetostatic modeling of tendon-driven parallel continuum robots," *IEEE Transactions on Robotics*, 2022.

[15] J. Ha, F. C. Park, and P. E. Dupont, "Elastic stability of concentric tube robots subject to external loads," *IEEE Transactions on Biomedical Engineering*, vol. 63, no. 6, pp. 1116–1128, 2015.

[16] R. Xu, S. F. Atashzar, and R. V. Patel, "Kinematic instability in concentric-tube robots: Modeling and analysis," in *5th IEEE RAS/EMBS International Conference on Biomedical Robotics and Biomechatronics*. IEEE, 2014, pp. 163–168.

[17] J. Ha, F. C. Park, and P. E. Dupont, "Optimizing tube precurvature to enhance the elastic stability of concentric tube robots," *IEEE Transactions on Robotics*, vol. 33, no. 1, pp. 22–37, 2016.

[18] C. Rucker, J. Childs, P. Molaei, and H. B. Gilbert, "Transverse anisotropy stabilizes concentric tube robots," *IEEE Robotics and Automation Letters*, vol. 7, no. 2, pp. 2407–2414, 2022.

[19] R. J. Hendrick, H. B. Gilbert, and R. J. Webster, "Designing snap-free concentric tube robots: A local bifurcation approach," in *2015 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2015, pp. 2256–2263.

[20] J. Till and D. C. Rucker, "Elastic stability of cosserat rods and parallel continuum robots," *IEEE Transactions on Robotics*, vol. 33, no. 3, pp. 718–733, 2017.

[21] F. Zaccaria, E. Idá, and S. Briot, "A Boundary Computation Algorithm for the Workspace Evaluation of Continuum Parallel Robots," *Journal of Mechanisms and Robotics*, vol. 16, no. 4, pp. 1–19, 2023.

[22] L. Euler, *Methodus inveniendi lineas curvas maximi minimive proprietae gaudentes sive solutio problematis isoperimetrii latissimo sensu accepti*. Springer Science & Business Media, 1952, vol. 1.

[23] J. Maddocks, "Stability of nonlinearly elastic rods," *Archive for rational mechanics and analysis*, vol. 85, no. 4, pp. 311–354, 1984.

[24] A. H. Nayfeh and S. A. Emam, "Exact solution and stability of postbuckling configurations of beams," *Nonlinear Dynamics*, vol. 54, pp. 395–408, 2008.

[25] A. Lazarus, J. Miller, and P. M. Reis, "Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method," *Journal of the Mechanics and Physics of Solids*, vol. 61, no. 8, pp. 1712–1736, 2013.

[26] S. Jiajia, L. Luke, W. M. Ahmer, G. Leroy, P. Alberto, and G. R. M. J., "Probing in situ capacities of prestressed stayed columns: towards a novel structural health monitoring technique," *Phil. Trans. R. Soc. A.*, vol. 381, p. 20220033, 2023.

[27] A. Abbasi, T. Sano, D. Yan, and P. Reis, "Snap buckling of bistable beams under combined mechanical and magnetic loading," *Phil. Trans. R. Soc. A.*, vol. 381, p. 20220029, 2023.

[28] D. C. Rucker, R. J. Webster III, G. S. Chirikjian, and N. J. Cowan, "Equilibrium conformations of concentric-tube continuum robots," *The International journal of robotics research*, vol. 29, no. 10, pp. 1263–1280, 2010.

[29] S. Briot and A. Goldsztejn, "Singularity conditions for continuum parallel robots," *IEEE Transactions on Robotics*, vol. 38, no. 1, pp. 507–525, 2021.

[30] H. B. Gilbert, R. J. Hendrick, and R. J. Webster III, "Elastic stability of concentric tube robots: A stability measure and design test," *IEEE Transactions on Robotics*, vol. 32, no. 1, pp. 20–35, 2015.

[31] Y. A. Kuznetsov, *Elements of Applied Bifurcation Theory*, 3rd ed. Springer New York, NY, 2004.

[32] O. Altuzarra, D. Caballero, F. J. Campa, and C. Pinto, "Position analysis in planar parallel continuum mechanisms," *Mechanism and Machine Theory*, vol. 132, pp. 13–29, Feb. 2019.

[33] O. Altuzarra and F. J. Campa, "On singularity and instability in a planar parallel continuum mechanism," in *Advances in Robot Kinematics 2020*. Springer, 2021, pp. 327–334.

[34] O. Altuzarra, M. Urizar, M. Cichella, and V. Petuya, "Kinematic analysis of three degrees of freedom planar parallel continuum mechanisms," *Mechanism and Machine Theory*, vol. 185, p. 105311, 2023.

[35] Q. Peyron and J. Burgner-Kahrs, "Stability analysis of tendon driven continuum robots and application to active softening," *IEEE Transactions on Robotics*, vol. 40, pp. 85–100, 2024.

[36] S. Briot and J.-P. Merlet, "Direct kinematic singularities and stability analysis of sagging cable-driven parallel robots," *IEEE Transactions on Robotics*, 2023.

[37] Q. Peyron, K. Rabenorosoa, N. Andreff, and P. Renaud, "A numerical framework for the stability and cardinality analysis of concentric tube robots: Introduction and application to the follow-the-leader deployment," *Mechanism and Machine Theory*, vol. 132, pp. 176–192, 2019.

[38] F. Zaccaria, E. Idá, S. Briot, and M. Carricato, "Workspace computation of planar continuum parallel robots," *IEEE Robotics and Automation Letters*, vol. 7, no. 2, pp. 2700–2707, 2022.

[39] M. Carricato and J.-P. Merlet, "Stability analysis of underconstrained cable-driven parallel robots," *IEEE Transactions on Robotics*, vol. 29, no. 1, pp. 288–296, 2012.

[40] K. Leibrandt, C. Bergeles, and G.-Z. Yang, "Concentric tube robots: Rapid, stable path-planning and guidance for surgical use," *IEEE Robotics & Automation Magazine*, vol. 24, no. 2, pp. 42–53, 2017.

[41] J.-P. Merlet, "Jacobian, Manipulability, Condition Number, and Accuracy of Parallel Robots," *Journal of Mechanical Design*, vol. 128, no. 1, pp. 199–206, 2005.

[42] K. Leibrandt, C. Bergeles, and G.-Z. Yang, "Implicit active constraints for safe and effective guidance of unstable concentric tube robots," in *2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, 2016, pp. 1157–1163.

[43] J. Till, C. E. Bryson, S. Chung, A. Orekhov, and D. C. Rucker, "Efficient computation of multiple coupled cosserat rod models for real-time simulation and control of parallel continuum manipulators," in *2015 IEEE international conference on robotics and automation (ICRA)*. IEEE, 2015, pp. 5067–5074.

[44] S. Patel and T. Sobh, "Manipulator performance measures-a comprehensive literature survey," *Journal of Intelligent & Robotic Systems*, vol. 77, pp. 547–570, 2015.

[45] F. Zaccaria, S. Briot, M. T. Chikhaoui, E. Idá, and M. Carricato, "An Analytical Formulation for the Geometrico-static Problem of Continuum Planar Parallel Robots," in *Symposium on Robot Design, Dynamics and Control, ROMANSY2020, G.Ventuare et al. (Eds.)*, no. 512–520, 2020.

[46] S. Antman, *Nonlinear Problems of Elasticity*. Springer Verlag New York, 1995, vol. 107.

[47] S. A. Davis H.F, *Introduction to Vector Analysis*, 4th ed. Allyn and Bacon, Inc, Toronto, 1979.

[48] C. Armanini, F. Boyer, A. T. Mathew, C. Duriez, and F. Renda, "Soft robots modeling: A structured overview," *IEEE Transactions on Robotics*, 2023.

[49] A. Gotelli, F. Zaccaria, O. Kermorgant, and S. Briot, "A gazebo simulator for continuum parallel robots," in *Advances in Robot Kinematics 2022*. Springer, 2022, pp. 248–256.

[50] R. J. Webster III and B. A. Jones, "Design and kinematic modeling of constant curvature continuum robots: A review," *The International Journal of Robotics Research*, vol. 29, no. 13, pp. 1661–1683, 2010.

[51] F. Renda, F. Boyer, J. Dias, and L. Seneviratne, "Discrete cosserat approach for multisection soft manipulator dynamics," *IEEE Transactions on Robotics*, vol. 34, no. 6, pp. 1518–1533, 2018.

[52] H. Li, L. Xun, and G. Zheng, "Piecewise linear strain cosserat model for soft slender manipulator," *IEEE Transactions on Robotics*, 2023.

[53] F. Boyer, V. Lebastard, F. Candelier, and F. Renda, "Dynamics of continuum and soft robots: A strain parameterization based approach," *IEEE Transactions on Robotics*, vol. 37, no. 3, pp. 847–863, 2020.

[54] J. Nocedal and S. Wright, *Numerical Optimization*, 2nd ed. Springer, 2006.

[55] S. Briot and A. Goldsztejn, "Technical report associated with the paper 'singularity conditions for continuum parallel robots'," CNRS, Laboratoire des Sciences du Numérique de Nantes, Tech. Rep. 03210572, 2020. [Online]. Available: https://hal.science/hal-03210572v1
[56] V. Aloi, C. Black, and C. Rucker, "Stiffness control of parallel continuum robots," in *Dynamic Systems and Control Conference*, vol. 51890. American Society of Mechanical Engineers, 2018, p. V001T04A012.

[57] E. Riks, "An incremental approach to the solution of snapping and buckling problems," *International Journal of Solids and Structures*, vol. 15, no. 7, pp. 529-551, 1979. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0020768379900817

[58] M. Crisfield, "A fast incremental/iterative solution procedure that handles "snap-through"," *Computers and Structures*, vol. 13, no. 1, pp. 55-62, 1981. [Online]. Available: https://www.sciencedirect.com/science/article/pii/0045794981901085

[59] J. R. Magnus and H. Neudecker, *Matrix differential calculus with applications in statistics and econometrics*. John Wiley & Sons, 2019.

[60] F. Zaccaria, E. Ida, and S. Briot, "Technical report associated with the paper: Directional critical load index: a distance-to-instability metric for continuum number," *Laboratoire des Sciences du Numérique de Nantes (LS2N)*, Tech. Rep., 2023. [Online]. Available: https://www.dropbox.com/scl/fi/y8xx41uy1p50quajs8k7p/Tech_Report_TRO_Stability_Index.pdf?rlkey=h2ps55hmb20xguq8c5bya0vsi&dl=0

[61] R. H. Byrd and R. B. Schnabel, "Continuity of the null space basis and constrained optimization," *Mathematical Programming*, vol. 35, no. 1, pp. 32-41, 1986.

[62] T. F. Coleman and A. Pothen, "The null space problem ii. algorithms," *SIAM Journal on Algebraic Discrete Methods*, vol. 8, no. 4, pp. 544-563, 1987.

[63] E. Idá, S. Briot, and M. Carricato, "Natural oscillations of underactuated cable-driven parallel robots," *IEEE Access*, vol. 9, pp. 71660-71672, 2021.

[64] S. Briot and F. Boyer, "A geometrically exact assumed strain modes approach for the geometric-and kinemato-static modelings of continuum parallel robots," *IEEE Transactions on Robotics*, 2022.

[65] F. Zaccaria, E. Idá, and S. Briot, "Design and experimental equilibrium stability assessment of a RFRFR continuum parallel robot," *Mechtronics*, vol. 95, p. 103064, 2023.

<figure><img src="image_17.png" /></figure>

**Sébastien Briot** received the B.S. and M.S. degrees in mechanical engineering in 2004 and the Ph.D. degree in robotics in 2007, from the Institut National des Sciences Appliquées de Rennes, Rennes, France. He was a Postdoctoral Fellow with the Ecole de Technologie Supérieure, Montreal, QC, Canada, in 2008. In 2009, he has been recruited at CNRS as a researcher in the Laboratoire des Sciences du Numérique de Nantes, Nantes, France, where he has been the Head of the ARMEN Research Team since 2017. Since 2022, he is CNRS Director of Research

in the same lab. He has authored more than 50 referred journal papers, two books, and four inventions. His research interests include the design optimization of robots and the analysis of their performance, especially their singularities. Dr. Briot received the Best Ph.D. Thesis Award in Robotics from the French CNRS in 2007. In 2011, he received two other awards: the Award for the Best Young Researcher from the French Region Bretagne

<figure><img src="image_18.png" /></figure>

**Federico Zaccaria** received both his B.Sc. and M.Sc. degrees in Mechanical Engineering at the University of Bologna in 2017 and 2020, respectively, with a specialization curriculum in Mechanics of Automation and Robotics. He was a research fellow at CIRI-MAM, Advanced Application in Mechanical Engineering and Materials Technology, where he collaborated on the development of a robotic manipulation system for the COORSA project. Federico collaborates with the Laboratoire des Sciences du Numérique de Nantes, where he is enrolled at the

SIS Ph.D program of the Ecole Centrale de Nantes. Its research topics, as a Ph.D. student of Mechanics and Advanced Engineering Science degree, are connected to industrial robotics and continuum parallel robot robots, with a focus on design, modelling, workspace analysis, and performance evaluation of flexible parallel manipulators.

<figure><img src="image_19.png" /></figure>

**Edoardo Idá** received both B.Sc. and M.Sc. degrees (with honors) in Mechanical Engineering at the University of Bologna. He was the winner of the graduation award in memory of Professor Ettore Funaioli in 2017, of the best research paper award at CableCon 2021, and the Young Author Best Paper Award of the Italian IEEE Robotic and Automation Society in 2021. Head of the lab, after completing his Ph.D. in Mechanics and Advanced Engineering Science, he is now a Junior Assistant Professor, teaches courses in automation and robotics at the

bachelor and master levels, and focuses his research on cable-driven robotics systems, continuum robots, automation, and mechanism design.