Analytical derivation of singularity-free tubes in the constant-orientation workspace of 6-6 Stewart platform manipulators

Abstract

This paper presents the novel concept of a singularity-free tube (SFT) in the constant orientation workspace of a spatial parallel manipulator. The concept is developed and demonstrated in the context of a $6$-$6$ spatial parallel manipulator, namely, the semi-regular Stewart platform manipulator. Given two points in the said workspace, the SFT is a tubular volume which contains these points and is free of gain-type or forward-kinematic singularities. The purpose of identifying such regions in space is to allow abundant freedom to the path-planner to connect the said points by a path, which can be free of gain-type singularities simply by remaining inside the SFT at all times. To demonstrate the concept, two smooth paths obtained by formulating two different optimisation problems have been presented as examples. The SFT can be of great help in singularity-free path-planning in many similar manipulators.

Nomenclature

Abbreviations

CAS

computer algebra system

DoF

degree(s)-of-freedom

LSFS

largest singularity-free sphere

SFS

singularity-free sphere

SFT, $\mathcal{T}$

singularity-free tube

SPM

Stewart platform manipulator

SRSPM

semi-regular Stewart platform manipulator

List of Symbols

${\boldsymbol{p}}_{\textrm{c}}(t)$

The centre-curve of the SFT

${\boldsymbol{p}}_{\textrm{c}_i}$

centre of $s_{i}$

${\boldsymbol{p}}_{\textrm{p}}(t)$

path connecting ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ , where $t$ is the path parameter

${\boldsymbol{p}}_{\textrm{o}}(t)$

path connecting ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ and having minimal length

$r_i$

radius of $s_{i}$

$\mathcal{G}$

geodesic curve connecting ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$

${\boldsymbol{G}}$

set of all points ${\boldsymbol{q}}_{i}$ sampled from $\mathcal{G}$

${\boldsymbol{q}}_{i}$

ith point sampled from $\mathcal{G}$

$s_{i}$

ith LSFS tangential to $\mathcal{S}$ at ${\boldsymbol{q}}_{i}$

${\boldsymbol{q}}_{\textrm{s}}$ , ${\boldsymbol{q}}_{\textrm{d}}$

projection of ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ onto $\mathcal{S}$ , respectively

$\mathcal{S}$

singularity surface in the constant-orientation space of the SRSPM

${\boldsymbol{p}}_{\textrm{s}},\,{\boldsymbol{p}}_{\textrm{d}}$

given source and destination points in the path-planning problem, respectively

$s_{\textrm{s}}$ , $s_{\textrm{d}}$

SFS centred at ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ , respectively

$\mathcal{W}$

constant-orientation workspace of the SRSPM

1. Introduction

Planning paths for spatial six degrees of freedom (DoF) manipulators has inherent technical challenges since $\mathbb{SE}(3)$ consists of $\mathbb{SO}(3)$ and ${\mathbb{R}}^3$ , which do not admit the same metrics globally. In the case of parallel manipulators, the problem gets complicated further due to additional issues of singularities of the gain-type ¹, (see, e.g., refs. [1] and [2]) which occur inside their workspaces. Various approaches have been proposed in the literature to solve this problem, specifically for manipulators of practical interest, such as the Stewart platform manipulator (SPM). A majority of these methods may be broadly categorised into two groups: (a) analytical and (b) randomised.

The work of Dasgupta and Mruthyunjaya reported in ref. [3] is a typical example of the first type, wherein a constrained optimisation-based analytical approach is used to find safe via-points in the workspace of the SPM through which a continuous singularity-free path may be constructed. A similar approach is presented in ref. [4], where a singularity-free path is planned by considering the limitations of the actuator forces. However, these approaches rely upon the determinants of certain Jacobian/transformation matrices in their formulations, which are known to be devoid of any physical significance since such matrices are dimensionally inhomogeneous. Further, a popular remedy to this problem,that is, homogenising such a matrix using a “characteristic length scale”, is also futile, since no such unique scale exists in $\mathbb{SE}(3)$ . This problem has been alleviated in ref. [5] through the use of an object-oriented metric and demonstrated in the context of linear pentapods having simple singular surfaces.

The randomised approaches, on the other hand, follow a common procedure of sampling the workspace in a discrete manner and planning a path inside it by considering the singularities within it as obstacles. For instance, in ref. [6], the workspace of the SPM is discretised and analysed using flood-filling algorithms and safe via-points for the singularity-free path are identified employing the probabilistic road map (PRM) algorithm. Similarly, in ref. [7], the workspace is discretised using slicing techniques and a singularity-free path is found by a local routing method based on line geometry. These techniques typically suffer from two major drawbacks: (a) large computational expenses, as they involve handling huge datasets generated as a consequence discretisation of the workspace and (b) the need for post-processing steps, such as enhancing the continuity of the singularity-free initial path, which may not guarantee the singularity-free nature of the final path.

While the above-mentioned works address the problem of planning a path between any two given points in $\mathbb{SE}(3)$ , a considerable amount of research has been done on planning such paths in the constant-orientation workspace ² of the SPM. In ref. [8], the geometric and algebraic properties of the singularity surface presented in ref. [9] have been used to plan a singularity-free path between two points in the constant-orientation workspace of the semi-regular Stewart platform manipulator (SRSPM). However, this method is only applicable when the terminal points of the desired path share the same $Z$ -coordinate.

Path planning in the constant-orientation workspace of the SPM finds practical applications in contour machining operations where the tool attached to the moving platform (MP) is required to maintain a constant orientation throughout the process. The SPMs used for such machining purposes are commonly known as hexapods, and examples of path-planning algorithms for these hexapods, specifically within their constant-orientation workspaces, may be found in refs. [10, 11]. An instance of path planning in the constant-orientation workspace of other 6-DoF spatial parallel manipulators may be found in ref. [12], wherein a method based on Hermite polynomials is presented for determining a singularity-free path of the least length for a 6-RSS parallel manipulator.

In most of the methods discussed above, an unique path is determined for each (successful) run of the algorithm. In other words, the entire set of computations needs to be repeated, in most cases, to find a path afresh, should either the source, or destination, or both be changed. Intuitively, it appears that if one could delimit all such singularity-free paths between two given points, it would provide the analyst with much greater flexibility in choosing their desired paths based on any combination of additional objectives and constraints. This thought process leads to the concept of safe-tubes, which serve as the bounding geometric object mentioned above.

The idea of a safe-tube appears to be fairly novel, as the authors could locate only two records in the existing literature wherein a similar concept has been presented, albeit not in the context of parallel manipulators. In ref. [13] the authors introduced a series of cubes termed as “bubbles”, that are free of singularities and obstacles, and computed them inside the configuration space of a serial manipulator. The union of these bubbles constituted the safe-tube. A safe tube bounded by a two-parametric surface has been constructed in ref. [14], once again, within the task space of a serial manipulator. Therefore, as far as parallel manipulators are concerned, there is no known mention of a safe tube of any kind in the existing literature (to the best of the authors’ knowledge).

In the present paper, a unique approach is proposed to compute a safe-tube (hereafter referred to as the singularity-free tube and abbreviated as SFT) containing the source and destination points, such that any path contained within this SFT would be free of gain-type singularities. It is to be noted that the SFT is described by a one-parameter family of spheres analytically, in the closed form. This method utilises the knowledge of the gain-type singularity surface, $\mathcal{S}$ , described implicitly as an algebraic surface in the variables chosen to parametrise $\mathbb{SE}(3)$ , as described in ref. [9]. In this paper, however, the formulation and results pertain to a given point in $\mathbb{SO}(3)$ since the orientation of MP is considered to be fixed. Once the SFT is obtained, any number of singularity-free paths may be obtained using the SFT as the container of the feasible paths and imposing additional constraints in accordance with the intended application. To illustrate the process two smooth mono-parametric paths from the source to the destination are presented, (a) path lies inside the SFT and is as far as possible from the boundary of the SFT and (b) path lies inside the SFT and has minimal length.

The derivation of the SFT involves computing a series of largest singularity-free spheres (LSFSs), which are tangential to the singularity surface at points on a geodesic curve connecting a pair of points obtained by projecting the given source and the destination points onto $\mathcal{S}$ . Further, mono-parametric cubic polynomials are fitted to the centres of these LSFSs, which are described by ${{\boldsymbol{p}}_c}(t)=\{x_{\textrm{c}}(t),y_{\textrm{c}}(t),z_{\textrm{c}}(t)\}$ , where $t \in [0,1]$ is the path parameter. The coefficients of these polynomials are obtained by solving an ordinary least-squares (OLS) (for more details, refer to p. $355$ of ref. [15]) problem after incorporating two boundary constraints, that is., the path originates from and terminates at the centres of the first and the last LSFS, respectively. Another polynomial in the same parameter $t$ , denoted by $r_{\textrm{c}}(t)$ , is obtained for computing the safe radii $\forall t \in (0,1)$ . The polynomials $\{x_{\textrm{c}}(t),y_{\textrm{c}}(t),z_{\textrm{c}}(t),r_{\textrm{c}}(t)\}$ collectively present a mono-parametric description of a family of spheres centred at ${{\boldsymbol{p}}_c}(t)$ , having radius $r_{\textrm{c}}(t)$ for $t\in [0,1]$ . This mono-parametric family of spheres represents the SFT (denoted by $\mathcal{T}$ ) analytically. The symbolic computations involved in these steps are performed in the computer algebraic system (CAS) Mathematica v13.0 [16].

The remainder of this paper is organised as follows: Section 2 describes the geometry of the SRSPM and singularity surface in the constant-orientation space of the SRSPM. Section 3 presents the steps followed in computing the SFT in the constant-orientation space of the SRSPM. Section 4 presents the numerical examples demonstrating the construction of the SFT and its utility in obtaining two singularity-free paths using different optimisation objectives. The main results and the contributions are summarised in Section 5, and the conclusions are presented in Section 6.

2. Geometry of the SRSPM and formulation of the singularity surface

The geometry of the SRSPM and the algebraic equation describing its singularities, adopted from ref. [9], are presented in the following.

2.1. Geometry of the SRSPM

A schematic of the SRSPM is shown in Figure 1. The manipulator consists of an MP driven by six UPS (i.e., a prismatic actuator with a universal joint with the fixed platform (FP) and a spherical joint with the MP) legs attached to an FP. The FP and MP are hexagonal in shape, with alternate sides having equal lengths. The prismatic joint of each leg is actuated, and its state is defined by its length $l_i, i=1,2,..,6$ . In some cases, the manipulator may have SPS (i.e., a prismatic actuator with spherical joints at both ends) legs. However, this leg architecture is kinematically equivalent to that of the UPS legs, as the roll motion of the SS links corresponds to an idle DoF.

Figure 1. Kinematic description of the SRSPM.

The vertices of the MP, denoted by ${\boldsymbol{t}}_1,{\boldsymbol{t}}_2,\dots, {\boldsymbol{t}}_6$ , lie on a circle of radius $r_{\textrm{t}}$ . The vertices of the FP, namely, ${\boldsymbol{b}}_1,{\boldsymbol{b}}_2,\dots, {\boldsymbol{b}}_6$ , also lie on a circle, whose radius is scaled to unity without any loss of generality. The angular spacings between the adjacent pairs of legs are given by $2{\gamma _{\textrm{b}}}$ and $2{\gamma _{\textrm{t}}}$ for the FP and MP, respectively. The global reference frame ${\boldsymbol{O}}$ - ${\boldsymbol{X}}{\boldsymbol{Y}}{\boldsymbol{Z}}$ , denoted by $\{1\}$ , is attached to the centroid of the FP, that is, ${\boldsymbol{O}}$ , such that the axis ${\boldsymbol{X}}$ passes through ${\boldsymbol{b}}_1$ . Similarly, the reference frame ${\boldsymbol{o}}$ - ${\boldsymbol{x}}{\boldsymbol{y}}{\boldsymbol{z}}$ , denoted by $\{2\}$ , is attached at the centroid of the MP, ${\boldsymbol{o}}$ , and ${\boldsymbol{x}}$ -axis passes through ${\boldsymbol{t}}_1$ . This particular choice of the coordinate frames results in vanishing of two of the coordinates of the first vertices of the FP and MP respectively (in the respective frames). This reduction of architecture parameters is found to help in reducing the complexity of subsequent symbolic computations. The axes ${\boldsymbol{Z}}$ and ${\boldsymbol{z}}$ are normal to the FP and MP, respectively. Under these specifications of the frames $\{1\}$ and $\{2\}$ , the angular displacement of the six legs, measured counter-clockwise (CCW) from the axes ${\boldsymbol{X}}$ and ${\boldsymbol{x}}$ , are given by $\boldsymbol{\gamma }_{\text{t}} = [0, 2{\gamma _{\textrm{t}}}, 2\pi /3, 2\pi /3+2{\gamma _{\textrm{t}}}, 4\pi /3, 4\pi /3+2{\gamma _{\textrm{t}}}]^{\top }$ and $\boldsymbol{\gamma }_{\text{b}} = [0, 2{\gamma _{\textrm{b}}}, 2\pi /3, 2\pi /3+2{\gamma _{\textrm{b}}}, 4\pi /3, 4\pi /3+2{\gamma _{\textrm{b}}}]^{\top }$ for the FP and MP, respectively. The states of the prismatic actuators, denoted by ${\boldsymbol{l}} = [l_1, l_2, l_3, l_4, l_5, l_6]^{\top }$ , define the active/actuated variables. The pose of the MP w.r.t. $\{1\}$ is defined by the coordinates of its centroid, namely, ${\boldsymbol{o}}=[x, y, z]^{\top }$ , and the orientation of $\{2\}$ w.r.t. $\{1\}$ , defined by the rotation matrix ${{}^1_2{\boldsymbol{R}}} \in{\mathbb{SO}(3)}$ . In the present work, ${}^1_2{\boldsymbol{R}}$ is expressed in terms of Rodrigues parameters (see, e.g., ref. [17]) of $\mathbb{SO}(3)$ , denoted by the variables $\{c_1, c_2, c_3\}$ . Thus, the task-space coordinates of the MP in frame $\{1\}$ are $\{x, y, z, c_1, c_2,c_3\}$ .

2.2. Algebraic description of the singularity surface

The singularity surface of the SRSPM has been derived analytically in the closed-from in ref. [9]. The equation defining the singularity surface is reproduced below:

(1) \begin{align} \notag f(x, y, z) \;:\!=\;\,&e_1 x^2 z + e_2 x^2 + e_3 x y z + e_4 x y + e_5 x z^2 + e_6 x z + e_7 x + e_8 y^2 z + e_9 y^2 + e_{10} y z^2+ e_{11} y z \\ &+ e_{12} y + e_{13} z^3 + e_{14} z^2 + e_{15} z + e_{16} = 0. \end{align}

The coefficients $e_i, i=1,2,\dots, 16,$ appearing in Eq. (1), are closed-form functions of the orientation variables $c_1, c_2, c_3$ and the architecture parameters ${r_{\textrm{t}}},{\gamma _{\textrm{b}}},{\gamma _{\textrm{t}}}$ . In this work, an SFT containing the source and destination points is constructed for an SRSPM whose architecture parameters and the orientation of its MP, that is, ${}^1_2{\boldsymbol{R}}$ , are known. Hence, in the context of this article, the coefficients $e_i$ may be treated as constants, and therefore, Eq. (1) describes the gain-type singularity surface denoted by $\mathcal{S}$ , which is the locus of all the singular points ${\boldsymbol{o}} = [x, y, z]^{\top } \in{\mathbb{R}}^3$ in the constant-orientation space of the SRSPM.

Figure 2. Sinularity-free path, ${\boldsymbol{p}}_{\textrm{p}}(t)$ connecting ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ , which is also verified to be contained in $\mathcal{W}$ post its computation. (Note: these pictures are for illustration purposes only).

3. Mathematical formulation of the singularity-Free tube

The objective of the present work is to enable singularity-free path-planning between the source, and the destination, denoted by ${{\boldsymbol{p}}_{\textrm{s}}} = [x_{\textrm{s}}, y_{\textrm{s}}, z_{\textrm{s}}]^{\top }$ and ${{\boldsymbol{p}}_{\textrm{d}}} = [x_{\textrm{d}}, y_{\textrm{d}}, z_{\textrm{d}}]^{\top }$ , respectively, such that, ${{\boldsymbol{p}}_{\textrm{s}}},{{\boldsymbol{p}}_{\textrm{d}}}\in{\mathcal{W}}$ , where, $\mathcal{W}$ denotes the constant-orientation workspace of the SRSPM. The central idea is to introduce the concept of the SFT, which is free of gain-type singularities and may be described analytically. This would enable the analyst to specify various objectives for the path-planning problem while treating $\mathcal{T}$ as the container of singularity-free paths. For the sake of completeness, this paper presents a smooth mono-parametric candidate path (denoted by ${\boldsymbol{p}}_{\textrm{p}}(t)\in{\mathcal{T}}$ ) connecting a given pair of source and destination points, which is obtained by maximising the distance from the boundary of $\mathcal{T}$ , along with the constraints ${\boldsymbol{p}}_{\textrm{p}}(0)={{\boldsymbol{p}}_{\textrm{s}}}$ and ${\boldsymbol{p}}_{\textrm{p}}(1)={{\boldsymbol{p}}_{\textrm{d}}}$ . Section 3.1 explains, in brief, the steps followed to find $\mathcal{T}$ in a parametric form. The mathematical details are presented subsequently.

Figure 3. Plot of $\mathcal{S}$ and $\mathcal{W}$ in the constant-orientation space of the SRSPM, for example, in Section 4.1. The linear path (shown in red colour) between ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ intersects $\mathcal{S}$ .

Figure 4. The projections of the starting and destination points and the geodesic curve joining these.

3.1. Steps followed in finding a singularity-free tube, $\mathcal{T}$

The proposed method builds upon some of the concepts used in ref. [18]. However, the main focus of the paper is to develop a parametric description of the SFT containing the source and destination points, which would be a key enabler for singularity avoidance during the path-planning process. In addition to finding the SFT, a candidate path joining the source and destination, which lies inside the SFT and is described in the form of mono-parametric cubic polynomials, is also presented. The steps are summarised below.

1. 1. Given ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ , first, it is ensured that these lie inside the workspace and on the same side of the singularity surface, $\mathcal{S}$ (this is ensured by checking that $f({{\boldsymbol{p}}_{\textrm{s}}})$ and $f({{\boldsymbol{p}}_{\textrm{d}}})$ are non-zero and have the same sign, that is, $f({{\boldsymbol{p}}_{\textrm{s}}})f({{\boldsymbol{p}}_{\textrm{d}}})\gt 0$ ). This is done very easily, by performing inverse kinematics and confirming that all the actuator lengths are within the given limits (depending on the specifications of the given SRSPM at hand). Next, the existence of the trivial solution, that is, a non-singular linear path connecting the ${\boldsymbol{p}}_{\textrm{d}}$ to the ${\boldsymbol{p}}_{\textrm{s}}$ , is examined. Such a path is given by:

   (2) \begin{align} {{\boldsymbol{p}}_{\textrm{p}}(\alpha )} ={{\boldsymbol{p}}_{\textrm{s}}} + \alpha ({{\boldsymbol{p}}_{\textrm{d}}} -{{\boldsymbol{p}}_{\textrm{s}}}), \quad \alpha \in [0, 1]. \end{align}
   Equation (2) presents the desired solution if ${{\boldsymbol{p}}_{\textrm{p}}(\alpha )}\notin{\mathcal{S}}\,\forall \alpha \in [0, 1]$ . If the linear path in Eq. (2) intersects $\mathcal{S}$ at one or more $\alpha \in [0, 1]$ (see, e.g., Figure 3(b)), then it fails to serve its purpose. In such cases, the general procedure is to be followed, as described below.

2. 2. To find a feasible path, firstly, ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ are projected onto $\mathcal{S}$ by computing the singularity-free spheres (SFSs), namely, $s_{\textrm{s}}$ and $s_{\textrm{d}}$ , centred at ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ , respectively. This is done by following the method presented in refs. [19, 20]. The point of tangency between $\mathcal{S}$ and $s_{\textrm{s}}$ , denoted by ${\boldsymbol{q}}_{\textrm{s}}$ , is the projection of ${\boldsymbol{p}}_{\textrm{s}}$ onto $\mathcal{S}$ (see Figure 4(a)). Similarly, ${\boldsymbol{q}}_{\textrm{d}}$ is the projection³ of ${\boldsymbol{p}}_{\textrm{d}}$ onto $\mathcal{S}$ .

   

   Figure 5. Largest SFS $s_{i}$ centred at ${\boldsymbol{p}}_{\textrm{c}_i}$ and with radius $r_i$ that touches $\mathcal{S}$ at ${\boldsymbol{q}}_{i}$ and ${\boldsymbol{p}}_{\textrm{b}_i}$ .

   

   Figure 6. Graphical depiction of the step $5$ followed in obtaining the gain-type singularity-free tube, $\mathcal{T}$ .

3. 3. The geodesic curve on $\mathcal{S}$ (denoted by $\mathcal{G}$ ) between ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$ is computed next. This step helps in (a) finding a continuous path connecting ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$ , and (b) guiding the subsequent steps towards computing the SFT containing the source and destination points. As $\mathcal{G}$ is computed numerically, only a set of $n$ discrete points on $\mathcal{G}$ (where, $n$ is an even positive integer), denoted by ${{\boldsymbol{G}}} = \{{{\boldsymbol{q}}_{i}} \in{\mathcal{G}}, i=1,2,\dots, n\}$ , are obtained, by sampling $\mathcal{G}$ uniformly (as per the length of $\mathcal{G}$ ), as shown in Figure 4(b).

4. 4. For each of the points ${\boldsymbol{q}}_{i}$ , an SFS is computed such that (a) it is tangential to $\mathcal{S}$ at ${\boldsymbol{q}}_{i}$ , and (b) it has the maximum possible radius. Such an SFS, denoted by $s_{i}$ , is called the largest SFS, and abbreviated as the LSFS. The LSFS $s_{i}$ tangential to $\mathcal{S}$ at ${\boldsymbol{q}}_{i}$ , computed using the method described in Section 3.4, is shown in Figure 5. The radius and centre of $s_{i}$ are given by $R_i$ and ${\boldsymbol{p}}_{\textrm{c}_i}$ , respectively. At this stage, a series of LSFSs $s_i$ that are tangential to $\mathcal{S}$ at ${\boldsymbol{q}}_{i}$ and ${\boldsymbol{p}}_{\textrm{b}_i}, i=1,2,\dots, n$ , are found.

5. 5. On obtaining the LSFSs, $s_{i}$ , one option for describing the SFT would be to obtain the envelope of $s_{i}$ . However, since $s_{i}$ are computed at discrete points, their envelope cannot be computed in a continuous manner. Therefore, portions of the SFT computed thus which fall between two consecutive LSFSs cannot be guaranteed to be free of singularities. It is desired to have a continuous algebraic description of the SFT. To obtain such a description, using the centres of the LSFSs $s_{i}$ , that is. ${\boldsymbol{p}}_{\textrm{c}_i}=\{x_{\textrm{c}_i},y_{\textrm{c}_i},z_{\textrm{c}_i}\}$ , an OLS problem is set up to estimate the coefficients of the set of cubic polynomials ${\boldsymbol{p}}_{\textrm{c}}(t)=\{x_{\textrm{c}}(t),y_{\textrm{c}}(t),z_{\textrm{c}}(t)\}$ , subject to two constraints ${\boldsymbol{p}}_{\textrm{c}}(0)={\boldsymbol{p}}_{\textrm{c}_0}$ and ${\boldsymbol{p}}_{\textrm{c}}(1)={\boldsymbol{p}}_{\textrm{c}_{\text{n}}}$ . The centre curve, that is, ${\boldsymbol{p}}_{\textrm{c}}(t)$ , may be expected to remain close to ${\boldsymbol{p}}_{\textrm{c}_i}$ , but need not pass through (all or any of) these. Thus, the radius of a sphere centred at each point on the centre curve would need to be recalculated, considering the deviation of the centre from the centres of the LSFSs. Thus, for each radius, $r_i$ of the $s_{i}$ , a corresponding $r_{\textrm{c}_i}$ is calculated as $r_{\textrm{c}_i}=r_i-\|{{\boldsymbol{p}}_{\textrm{c}_i}-{\boldsymbol{p}}_{\textrm{c}}(\frac{i}{n})}\|$ , such that the obtained sphere lies inside the corresponding $s_{i}$ (as shown in Figure 6(a)) and is, therefore, also free of singularities. Similar to ${\boldsymbol{p}}_{\textrm{c}}(t)$ , a constrained optimisation problem is set up using $r_{\textrm{c}_i}$ to obtain the coefficients of a cubic polynomial $r_{\textrm{c}}(t)$ . The four polynomials $\{x_{\textrm{c}}(t),y_{\textrm{c}}(t),z_{\textrm{c}}(t),r_{\textrm{c}}(t)\}$ describe the mono-parametric family of spheres which form the SFT, $\mathcal{T}$ , as visualised in Figure 6(b).

6. 6. The analytical description of the SFT is the main focus of the paper. However, for completeness, in Section 3.6, a smooth path, ${\boldsymbol{p}}_{\textrm{p}}(t)$ with ${\boldsymbol{p}}_{\textrm{p}}(0)={{\boldsymbol{p}}_{\textrm{s}}}$ and ${\boldsymbol{p}}_{\textrm{p}}(1)={{\boldsymbol{p}}_{\textrm{d}}}$ is obtained. The path is obtained by solving a constrained optimisation problem, that is, minimising the distance from the centre curve, ${\boldsymbol{p}}_{\textrm{c}}(t)$ while ensuring that the path lies inside the SFT. A Numerical example showing the computation of a gain-type singularity-free path obtained using the knowledge of the derived SFT is shown in Section 4.

3.2. Computation of the SFS

The details of the second step, namely, the computation of the SFS, are omitted as these may be found in ref. [19].

3.3. Computation of the geodesic curve, $\mathcal{G}$ on the singularity surface, $\mathcal{S}$

The computation of a geodesic curve is a fairly standard procedure, particularly when a parametric description of the surface is available (see, e.g., refs. [21], pp. 140–165, [22], pp. 127–136). As shown by Coste and Moussa in ref. [23], the singularity surface $\mathcal{S}$ does admit a rational parametrisation in terms of two parameters, $\{u, v\}$ . However, the denominator polynomial describes a rectangular hyperbola in the parameter space (see Appendix C). Therefore, the parameter space is split into disjoint subsets (see Figure 13). Consequently, such a parametrisation is not very suitable for computing the geodesic curves connecting an arbitrary pair of points on $\mathcal{S}$ . To circumvent this issue, a different formulation utilising the implicit form of $\mathcal{S}$ in Eq. (1) is considered, wherein the geodesic equations are derived by using the geometric property that the geodesic curvature vanishes at all points on the geodesic curve (see, e.g., ref. [22], p. 131).

By definition of the geodesic, ${\mathcal{G}}: \,{{\boldsymbol{q}}(\beta )} ={[x(\beta ), y(\beta ), z(\beta )]}^{\top }$ represents a smooth curve on the implicit surface $f(x, y, z) = 0$ , parametrised by the curve-length $\beta \in [\beta _1, \beta _2]$ . At a generic point ${\boldsymbol{q}}(\beta )$ , the tangent, the normal, and the binormal vectors are given by ${{\boldsymbol{q}}'},\,{{\boldsymbol{n}}_{\textrm{g}}}$ and ${\boldsymbol{a}}_{\textrm{g}}$ , respectively, where:

(3) \begin{align} {{\boldsymbol{q}}'} \;:\!=\;{\frac{\textrm{d}{{\boldsymbol{q}}(\beta )}}{\textrm{d} \beta }}, \end{align}

(4) \begin{align} {{\boldsymbol{n}}_{\textrm{g}}} \;:\!=\; \nabla _{{\boldsymbol{q}}}f({{\boldsymbol{q}}(\beta )}), \end{align}

(5) \begin{align} {{\boldsymbol{a}}_{\textrm{g}}} \;:\!=\;{{\boldsymbol{q}}'}\times{{\boldsymbol{n}}_{\textrm{g}}}. \end{align}

It may be noted that Eq. (4) derives from the fact that at any point on a geodesic curve the normal to the curve coincides with the normal (i.e., the gradient) to the surface at the same point (see, e.g. ref. [24]). It is also important to note that for the purpose of this derivation $\mathcal{S}$ is considered to be devoid of singular points, that is, ${{\boldsymbol{n}}_{\textrm{g}}} \ne{\boldsymbol{0}}$ at all points. Indeed, this formulation breaks down when $\nabla _{{\boldsymbol{q}}} f({{\boldsymbol{q}}(\beta )}) ={\boldsymbol{0}}$ , which happens only at a finite number of isolated points on $\mathcal{S}$ . These points have been identified in Appendix B.1.

The vectors ${{\boldsymbol{q}}'},\,{{\boldsymbol{n}}_{\textrm{g}}}$ and ${\boldsymbol{a}}_{\textrm{g}}$ in Eqs. (3)–(5) form the axes of the Frenet-Serret local frame of reference (see, e.g., ref. [22], pp. 5–18) attached to ${\boldsymbol{q}}(\beta )$ (denoted by $\{F_1\}$ ) as shown in Figure 7. The osculating plane is spanned by ${\boldsymbol{q}}'$ and ${\boldsymbol{n}}_{\textrm{g}}$ (denoted by $\mathcal{A}_{\textrm{g}}$ , see Figure 7). The curvature vector at ${\boldsymbol{q}}(\beta )$ is given by ${\boldsymbol{q}}''$ :

(6) \begin{align} &{{\boldsymbol{q}}''} \;:\!=\;{\frac{\textrm{d}^2{{\boldsymbol{q}}(\beta )}}{\textrm{d} \beta ^2}}. \end{align}

As described in ref. [22], p. 131, since $\mathcal{G}$ is a geodesic curve on the surface $f(x, y, z) = 0$ the geodesic curvature at ${\boldsymbol{q}}(\beta )$ , denoted by ${\kappa _{\textrm{g}}}(\beta )$ , vanishes at each point on the curve:

(7) \begin{align} &{\kappa _{\textrm{g}}}(\beta ) \;:\!=\;{{\boldsymbol{a}}_{\textrm{g}}}^{\top }{{\boldsymbol{q}}''} = 0, \quad \forall \beta \in [\beta _1, \beta _2]. \end{align}

Equivalently, from Eq. (7), it can be seen that as ${\boldsymbol{q}}(\beta )$ is a point on the geodesic curve, the curvature vector at that point, ${\boldsymbol{q}}''$ , is contained in the osculating plane, $\mathcal{P}_{\textrm{g}}$ . Furthermore, ${\boldsymbol{q}}''$ is along ${\boldsymbol{n}}_{\textrm{g}}$ (see Figure 7), since, by definition, ${\boldsymbol{q}}''$ is orthogonal to ${\boldsymbol{q}}'$ :

(8) \begin{align} &\left ({{\boldsymbol{q}}'}\right )^{\top }{{\boldsymbol{q}}''} = 0. \end{align}

The vectors ${\boldsymbol{q}}'$ and ${\boldsymbol{q}}''$ also satisfy the second-order constraint:

(9) \begin{align} &{\frac{\textrm{d}^2 }{\textrm{d} \beta ^2 }} f({{\boldsymbol{q}}(\beta )}) = \left ({{\boldsymbol{q}}'}\right )^{\top }{{\boldsymbol{H}}_{\textrm{g}}}{{\boldsymbol{q}}'} +{{\boldsymbol{n}}_{\textrm{g}}}^{\top }{{\boldsymbol{q}}''} = 0, \quad \textrm{where,}\nonumber \\ &{{\boldsymbol{H}}_{\textrm{g}}} = \begin{bmatrix} f_{xx} & f_{xy} & f_{xz} \\ f_{yx} & f_{yy} & f_{yz} \\ f_{zx} & f_{zy} & f_{zz} \\ \end{bmatrix},\,\textrm{is the } \textit{Hessian} \textrm{of}\;f\; \textrm{w.r.t.}\,\{x, y, z\}\textrm{ and}\ f_x = \frac{\partial{f}}{\partial{x}}, \ f_{xx} ={\frac{\partial ^2 f}{\partial x^2}}, \ \textrm{etc.} \end{align}

Equations (7)–(9) form a system of three linear equations in ${\boldsymbol{q}}''$ that may be organised as:

\begin{align*} &{{\boldsymbol{A}}_{\textrm{g}}}{{\boldsymbol{q}}''} ={{\boldsymbol{b}}_{\textrm{g}}}, \quad \textrm{where} \ {{\boldsymbol{A}}_{\textrm{g}}} ={\left [{{\boldsymbol{a}}_{\textrm{g}}}, \, \,{{\boldsymbol{q}}'}, \, \,{{\boldsymbol{n}}_{\textrm{g}}} \right ]}^{\top },{{\boldsymbol{b}}_{\textrm{g}}} ={\left [0, \, \, 0, \, \, -\left ({{\boldsymbol{q}}'}\right )^{\top }{{\boldsymbol{H}}_{\textrm{g}}}{{\boldsymbol{q}}'} \right ]}^{\top }. \end{align*}

In Eq. (10), ${\boldsymbol{A}}_{\textrm{g}}$ defines a transformation from the global frame $\{1\}$ to the Frenet–Serret local frame $\{F_1\}$ . Solving Eq. (10), ${\boldsymbol{q}}''$ is found using the reciprocal basis (see, e.g., ref. [25], pp. 11–18) of ${\boldsymbol{A}}_{\textrm{g}}$ :

(10) \begin{align} {{\boldsymbol{q}}''} ={{\boldsymbol{A}}_{\textrm{g}}}^{-1}{{\boldsymbol{b}}_{\textrm{g}}}, \quad \textrm{where,} \end{align}

(11) \begin{align} {{\boldsymbol{A}}_{\textrm{g}}}^{-1} &= \frac{1}{D}\left [{{\boldsymbol{a}}_{\textrm{g}}}, \, \,{{\boldsymbol{n}}_{\textrm{g}}}\times{{\boldsymbol{a}}_{\textrm{g}}}, \, \,{{\boldsymbol{a}}_{\textrm{g}}}\times{{\boldsymbol{q}}'} \right ],\,\text{and} \nonumber \\[5pt] D &= \left ({{\boldsymbol{a}}_{\textrm{g}}}\times{{\boldsymbol{q}}'}\right )\cdot{{\boldsymbol{n}}_{\textrm{g}}} ={{\boldsymbol{a}}_{\textrm{g}}}^{\top }{{\boldsymbol{a}}_{\textrm{g}}} \ \text{(assumed to be nonzero)}; \end{align}

(12) \begin{align} \implies {{\boldsymbol{q}}''} = \left (\frac{{\left ({{\boldsymbol{q}}'}\right )}^{\top }{{\boldsymbol{H}}_{\textrm{g}}}{{\boldsymbol{q}}'}}{D}\right ){{\boldsymbol{q}}'}\times{{\boldsymbol{a}}_{\textrm{g}}}, \ \text{assuming}. \end{align}

Equation (12) defines a system of three ordinary differential equations that are satisfied by the geodesic curve $\mathcal{G}$ on the implicit surface $f(x, y, z) = 0$ . Therefore, given two points ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$ on $\mathcal{S}$ , the geodesic curve $\mathcal{G}$ between them may be obtained by solving Eq. (12) along with the boundary conditions ${\boldsymbol{q}}(0) ={{\boldsymbol{q}}_{\textrm{s}}}$ and ${\boldsymbol{q}}(1) ={{\boldsymbol{q}}_{\textrm{d}}}$ , assuming that $D \ne 0$ throughout the curve.

Figure 7. Depiction of the Frenet–Serret local frame $\{F_1\}$ attached to ${\boldsymbol{q}}(\beta )$ on the geodesic curve $\mathcal{G}$ lying on the implicit surface $f(x, y, z) = 0$ . The axes of $\{F_1\}$ are defined by the tangent, ${\boldsymbol{q}}'$ , normal, ${\boldsymbol{n}}_{\textrm{g}}$ , and binormal, ${\boldsymbol{a}}_{\textrm{g}}$ , vectors at ${\boldsymbol{q}}(\beta )$ . The curvature vector at ${\boldsymbol{q}}(\beta )$ , i.e., ${\boldsymbol{q}}''$ , is along ${\boldsymbol{n}}_{\textrm{g}}$ and the osculating plane, $\mathcal{A}_{\textrm{g}}$ , is spanned by ${\boldsymbol{q}}'$ and ${\boldsymbol{n}}_{\textrm{g}}$ .

If, on the other hand, $D = 0$ , then from Eqs. (5) and (11) ${\boldsymbol{a}}_{\text{g}} = {\boldsymbol{q}}'\times {\boldsymbol{n}}_{\text{g}} ={\boldsymbol{0}}$ . Further, since ${\boldsymbol{q}}'$ belongs to the tangent space of $\mathcal{S}$ at ${\boldsymbol{q}}(\beta )$ , it is orthogonal to ${\boldsymbol{n}}_{\textrm{g}}$ ; that is, ${\boldsymbol{q}}'^{\top } {\boldsymbol{n}}_{\text{g}} = 0$ . These two equations can only be satisfied simultaneously if one or both of ${\boldsymbol{q}}'$ and ${\boldsymbol{n}}_{\text{g}}$ vanish. Vanishing of ${\boldsymbol{q}}'$ at a point would imply that the geodesic curve passes through the point and suffers a singularity in itself at that point. Thus, vanishing of ${\boldsymbol{n}}_{\text{g}}$ is of greater significance, since in that case, it is the surface $\mathcal{S}$ that contains a singularity at such a point, and as a consequence thereof, no geodesic curve passing through that point can exist. Therefore, one must identify the points at which ${\boldsymbol{n}}_{\text{g}} ={\boldsymbol{0}}$ ; or, equivalently, the gradient vector $\nabla _{{\boldsymbol{q}}} f({\boldsymbol{q}}(\beta ))$ vanishes (see Eq. (4)), as the formulation breaks down at such singular points. These points are rare and finite in number, and they can occur only at specific combinations of architecture parameters and orientations of the MP, as shown in Appendix B. Examples of such points can be seen in Appendix B.1, Figures 12(a) and (b). As mentioned before, these points are to be avoided in the rest of this work, which is possible by identifying them a priori.

3.4. Determination of the largest singularity-free sphere (LSFS) given one point of Tangency with the singularity surface, $\mathcal{S}$

The geodesic curve computed above serves as a guide, as the SFT is essentially derived using the knowledge of the largest singularity-free spheres, or LSFSs, which are tangential to $\mathcal{S}$ at each point on the said curve. In the following, the computation of the LSFS is explained, which is a harder variant of the SFS problem presented in refs. [19, 20]. In these works, the centre of the SFS was specified; however, in the present problem, the centre of the LSFS is unknown, but the point of tangency between it and the $\mathcal{S}$ (denoted ${\boldsymbol{p}}_{\textrm{a}}$ ) is given.

Since the SFS should not intersect $\mathcal{S}$ , it may be grown hypothetically until it touches $\mathcal{S}$ at a point ${\boldsymbol{p}}_{\textrm{b}}$ , all the while remaining tangent to $\mathcal{S}$ at ${\boldsymbol{p}}_{\textrm{a}}$ . Given the nonlinear nature of $\mathcal{S}$ , there can exist multiple points of tangency between $\mathcal{S}$ and the LSFS other than ${\boldsymbol{p}}_{\textrm{a}}$ . Thus, the problem reduces to finding a sphere that touches the singularity surface at least at two distinct points, one of them being ${\boldsymbol{p}}_{\textrm{a}}$ . Among the solutions for the spheres that touch $\mathcal{S}$ at least at two distinct points, it is obvious that the LSFS is the sphere with the least radius since the other spheres would intersect $\mathcal{S}$ at other points.

3.4.1. Mathematical formulation of the largest singularity-free sphere (LSFS)

Let, a sphere, namely $s_{\textrm{g}}$ , of radius $R$ , centred at ${\boldsymbol{p}}_{\textrm{c}}$ , be tangential to $\mathcal{S}$ at a given point ${\boldsymbol{p}}_{\textrm{a}}$ , as well as another unknown point, ${{\boldsymbol{p}}_{\textrm{b}}} = [x, y, z]^{\top } \in{\mathbb{R}}^3$ . Such a sphere may be described by

(13) \begin{align} &g \;:\!=\; ({{\boldsymbol{p}}_{\textrm{b}}} -{{\boldsymbol{p}}_{\textrm{c}}})\cdot ({{\boldsymbol{p}}_{\textrm{b}}} -{{\boldsymbol{p}}_{\textrm{c}}}) - R^2 = 0. \end{align}

The unit normal vectors to $\mathcal{S}$ (given by Eq. (1)) and $s_{\textrm{g}}$ at ${\boldsymbol{p}}_{\textrm{a}}$ may be computed as, respectively:

(14) \begin{align} {\boldsymbol{n}}_1 \;:\!=\; \left . \frac{\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} f}{\|{\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} f}\|}\right \rvert _{{{\boldsymbol{p}}_{\textrm{a}}}}, \end{align}

(15) \begin{align} {\boldsymbol{n}}_2 \;:\!=\; \left . \frac{\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} g}{\|{\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} g}\|}\right \rvert _{{{\boldsymbol{p}}_{\textrm{a}}}} = \frac{({{\boldsymbol{p}}_{\textrm{a}}} -{{\boldsymbol{p}}_{\textrm{c}}})}{R}. \end{align}

Since ${\boldsymbol{p}}_{\textrm{a}}$ is a point of tangency between $\mathcal{S}$ and $s_{\textrm{g}}$ , the unit normal vectors in Eqs. (14) and (15) align; that is:

(16) \begin{align} &{\boldsymbol{n}}_1 = \pm {\boldsymbol{n}}_2. \end{align}

From Eqs. (14)–(16), the relation between ${\boldsymbol{p}}_{\textrm{c}}$ , ${\boldsymbol{p}}_{\textrm{a}}$ and $R$ may be derived as:

(17) \begin{align} &{{\boldsymbol{p}}_{\textrm{c}}} ={{\boldsymbol{p}}_{\textrm{a}}} + r{\boldsymbol{n}}_1, \quad \textrm{where,} \quad r = \pm R. \end{align}

Figure 8. Two spheres with same radius $R$ that touch $\mathcal{S}$ at ${\boldsymbol{p}}_{\textrm{a}}$ and lie on either side of $\mathcal{S}$ w.r.t. ${\boldsymbol{p}}_{\textrm{a}}$ .

Equation (17) shows that the centre of $s_{\textrm{g}}$ is constrained to lie on the line passing through ${\boldsymbol{p}}_{\textrm{a}}$ and directed along ${\boldsymbol{n}}_1$ , as the radius of $s_{\textrm{g}}$ , that is, $R$ , varies. Furthermore, for a given $R$ , there exist two spheres whose centres lie on either side of $\mathcal{S}$ , both tangential to it at ${\boldsymbol{p}}_{\textrm{a}}$ (see Figure 8). In further computations towards the derivation of the SFT, the sphere whose centre lies on the same side of $\mathcal{S}$ as ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ is used. Using Eq. (17), the equation of $s_{\textrm{g}}$ in Eq. (13) may be rewritten as:⁴

(18) \begin{align} &g \;:\!=\; ({{\boldsymbol{p}}_{\textrm{b}}} - ({{\boldsymbol{p}}_{\textrm{a}}} + r{\boldsymbol{n}}_1))\cdot ({{\boldsymbol{p}}_{\textrm{b}}} - ({{\boldsymbol{p}}_{\textrm{a}}} + r{\boldsymbol{n}}_1)) - r^2 = 0. \end{align}

Similar to Eq. (16), the condition for tangency between $\mathcal{S}$ and $s_{\textrm{g}}$ at the unknown point ${\boldsymbol{p}}_{\textrm{b}}$ may be written as

(19) \begin{align} &\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} f ({{{\boldsymbol{p}}_{\textrm{b}}}}) = \lambda _{{{\boldsymbol{p}}_{\textrm{b}}}} \nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} g({{{\boldsymbol{p}}_{\textrm{b}}}}), \quad \lambda _{{{\boldsymbol{p}}_{\textrm{b}}}} \in{\mathbb{R}}. \end{align}

The unknown multiplier $\lambda _{{{\boldsymbol{p}}_{\textrm{b}}}}$ may be eliminated by computing the cross-product of $\nabla f ({{{\boldsymbol{p}}_{\textrm{b}}}})$ with both sides of Eq. (19):

(20) \begin{align} & \nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} f({{{\boldsymbol{p}}_{\textrm{b}}}}) \times \nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} g ({{{\boldsymbol{p}}_{\textrm{b}}}}) = \mathbf{0}. \end{align}

Using Eq. (13), Eq. (20) may be written in terms of ${\boldsymbol{p}}_{\textrm{b}}$ and ${\boldsymbol{p}}_{\textrm{c}}$ as

(21) \begin{align} &\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} f({{{\boldsymbol{p}}_{\textrm{b}}}}) \times ({{\boldsymbol{p}}_{\textrm{b}}} -{{\boldsymbol{p}}_{\textrm{c}}}) = \mathbf{0}. \end{align}

Finally, substituting ${\boldsymbol{p}}_{\textrm{c}}$ from Eq. (17) in Eq. (21) yields:

(22) \begin{align} &\nabla _{{{\boldsymbol{p}}_{\textrm{b}}}} f({{{\boldsymbol{p}}_{\textrm{b}}}}) \times ({{\boldsymbol{p}}_{\textrm{b}}} -{{\boldsymbol{p}}_{\textrm{a}}} - r{\boldsymbol{n}}_1) = \mathbf{0}, \end{align}

Equation (22) comprises three polynomial equations in the variables ${{\boldsymbol{p}}_{\textrm{b}}} = [x, y, z]^{\top }$ and $r$ :

(23) \begin{align} &h_i(x, y, z, r) = 0, \quad i = 1,2,3. \end{align}

At the most, only two of the equations among $h_i = 0$ in Eq. (23) may be linearly independent and hence, in the general case, any two of the three may be considered. To simplify further analysis, the equations corresponding to the expressions $h_i$ that have the least degree in the variables $x, y, z, r$ and size ⁵ are considered. The degrees and sizes of the expressions $h_i$ in the variables $x, y, z, r$ are summarised in Table IV, where it can be seen that $h_3$ has the least degree and size among $h_i, i=1,2,3$ . Furthermore, since $h_1$ and $h_2$ are observed to contain identical power-products ⁶ in the variables $x, y, z, r$ , the expression with the smaller size among the two, that is, $h_1$ , is chosen. Thus, the two equations among $h_i = 0$ considered in the following are $h_1 = 0$ and $h_3 = 0$ .

Additionally, ${\boldsymbol{p}}_{\textrm{b}}$ satisfies the equation defining $\mathcal{S}$ (i.e., Eq. (1)) and ${{\boldsymbol{p}}_{\textrm{b}}}, r$ both satisfy the equation defining $s_{\textrm{g}}$ (in Eq. (18)). Considering all of these together, one obtains a system of four independent equations in terms of the variables $x, y, z$ , and $r$ , which may be expressed as

(24) \begin{align}f(x, y, z) = 0, \end{align}

(25) \begin{align} g(x, y, z, r) = 0, \end{align}

(26) \begin{align} h_1(x, y, z, r) = 0, \end{align}

(27) \begin{align} h_3(x, y, z, r) = 0. \end{align}

The above system of equations may be denoted by ${\boldsymbol{f}}={\boldsymbol{0}}$ , where ${\boldsymbol{f}}=\{f,g,h_1,h_3\}$ . The corresponding Jacobian matrix, denoted by ${\boldsymbol{J}}$ , is given by

(28) \begin{align} &{\boldsymbol{J}}= \frac{\partial{{\boldsymbol{f}}}}{\partial{\{x,y,z,r\}}}. \end{align}

Solving the system of polynomial equations ${\boldsymbol{f}}={\boldsymbol{0}}$ by replacing ${\boldsymbol{p}}_{\text{a}}$ by ${\boldsymbol{q}}_{i}$ yields the solution for the $s_{i}$ corresponding to ${{\boldsymbol{q}}_{i}}\in{{\boldsymbol{G}}}$ , with a radius $r_i$ . The corresponding centre ${\boldsymbol{p}}_{\textrm{c}_i}=\{x_{\textrm{c}_i}, y_{\textrm{c}_i}, z_{\textrm{c}_i}\}$ may be calculated using Eq. (17). The solution procedure is detailed in Appendix A.

3.5. Construction of the singularity-free tube (SFT)

The solution procedure presented in Appendix A is used to compute the first and last LSFSs, that is, $s_{0}$ and $s_{n}$ , that are tangential to $\mathcal{S}$ at the points ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$ , respectively. The first step in constructing the SFT is to obtain a series of $n$ LSFSs, denoted by $s_{i}$ , which are tangent to $\mathcal{S}$ at each point, ${\boldsymbol{q}}_i \in{{\boldsymbol{G}}},\,i=1,2,\dots, n$ . These are obtained using a Newton–Raphson method-based iteration scheme using the system of equations ${\boldsymbol{f}}={\boldsymbol{0}}$ and the Jacobian matrix, ${\boldsymbol{J}}$ , given in Eq. (28). The detailed procedure to obtain the solutions for ${\boldsymbol{p}}_{\text{b}_i}$ corresponding to each ${\boldsymbol{q}}_{i}$ is shown in Method 1. The procedure is repeated twice, once starting with $s_{0}$ and iterating for $i=1,2,\dots, n/2$ and a second time starting with $s_{n}$ and iterating for $i=n-1,n-2,\dots, n/2$ . To verify the accuracy, the radii and the centre of the $s_{n/2}$ obtained from both computations are compared and are found to be identical up to five decimal places. After ${\boldsymbol{p}}_{\text{b}_i}$ are obtained the corresponding ${\boldsymbol{p}}_{\textrm{c}_i}$ may be computed using Eq. (17).

Algorithm 1. Obtaining the series of LSFSs using Newton-Raphson iterations

3.5.1. A mono-parametric polynomial as the centre curve of the singularity-free tube (SFT)

Having computed a series of LSFSs, that is, $s_{i}$ centred at ${\boldsymbol{p}}_{\textrm{c}_i}$ with the corresponding radii $r_i,$ $i=1,2,\dots, n$ , the next goal is to describe a volume which is free of gain-type singularities. As mentioned in Section 3.1, such a volume is named the singularity-free tube (SFT) and is denoted by $\mathcal{T}$ .

As explained in Section 3.1 a polynomial description for the centre curve is chosen in this paper. Towards this, the coefficients of a mono-parametric cubic curve given by ${\boldsymbol{p}}_{\textrm{c}}(t)=\{x_{\textrm{c}}(t),y_{\textrm{c}}(t),z_{\textrm{c}}(t)\}$ are obtained by formulating an OLS problem with the objective of minimising the deviation from the centres ${\boldsymbol{p}}_{\textrm{c}_i}$ . This curve would form the centre curve of the SFT. Let the cubic polynomials be given by

(29) \begin{align}x_{\textrm{c}}(t)=a_{10} t^3 + a_{11} t^2 + a_{13} t + a_{14}, \end{align}

(30) \begin{align}y_{\textrm{c}}(t)=a_{20} t^3 + a_{21} t^2 + a_{23} t + a_{24}, \end{align}

(31) \begin{align}z_{\textrm{c}}(t)=a_{30} t^3 + a_{31} t^2 + a_{33} t + a_{34}, \quad a_{ij}\in{\mathbb{R}}\,\text{and}\,t\in [0,1]. \end{align}

The constraints ${\boldsymbol{p}}_{\textrm{c}}(0)={\boldsymbol{p}}_{\textrm{c}_0}$ and ${\boldsymbol{p}}_{\textrm{c}}(1)={\boldsymbol{p}}_{\textrm{c}_n}$ , may be expressed as

(32) \begin{align}x_{\textrm{c}}(0) = a_{14} = x_{\textrm{c}_0},\,\quad \,x_{\textrm{c}}(1) = a_{10} + a_{11} + a_{13} + x_{\textrm{c}_0} = x_{\textrm{c}_n}, \end{align}

(33) \begin{align}y_{\textrm{c}}(0) = a_{24} = y_{\textrm{c}_0},\,\quad \,y_{\textrm{c}}(1) = a_{20} + a_{21} + a_{22} + y_{\textrm{c}_0} = y_{\textrm{c}_n}, \end{align}

(34) \begin{align}z_{\textrm{c}}(0) = a_{34} = z_{\textrm{c}_0},\,\quad \,z_{\textrm{c}}(1) = a_{30} + a_{31} + a_{32} + z_{\textrm{c}_0} = z_{\textrm{c}_n}. \end{align}

By absorbing the constraints shown in Eqs. (32)–(34) into Eqs. (29)–(31), the number of unknowns may be reduced to two for each of the cubic polynomials and ${\boldsymbol{p}}_{\textrm{c}}(t)$ may be rewritten as shown in Eqs. (35)–(37):

(35) \begin{align}x_{\textrm{c}}(t) = a_{10} t^3 + a_{11} t^2 + (x_{\textrm{c}_n} -a_{10}-a_{11}-x_{\textrm{c}_0})t + x_{\textrm{c}_0}, \end{align}

(36) \begin{align}y_{\textrm{c}}(t) = a_{20} t^3 + a_{21} t^2 + (z_{\textrm{c}_n} -a_{20}-a_{21}-y_{\textrm{c}_0})t + y_{\textrm{c}_0}, \end{align}

(37) \begin{align}z_{\textrm{c}}(t) = a_{30} t^3 + a_{31} t^2 + (z_{\textrm{c}_n} -a_{30}-a_{31}-z_{\textrm{c}_0})t + z_{\textrm{c}_0}. \end{align}

Equations (35)–(37) may be rewritten as below by isolating the terms containing the unknowns:

(38) \begin{align} &{\boldsymbol{p}}_{\textrm{c}}(t)=\begin{bmatrix} x_{\textrm{c}}(t)\\ y_{\textrm{c}}(t)\\ z_{\textrm{c}}(t) \end{bmatrix}^{\top }=\left [ (t^3-t),(t^2-t) \right ]\boldsymbol{\theta } +\begin{bmatrix} (x_{\textrm{c}_n} -x_{\textrm{c}_0})t + x_{\textrm{c}_0}\\ (y_{\textrm{c}_n} -y_{\textrm{c}_0})t + y_{\textrm{c}_0}\\ (z_{\textrm{c}_n} -z_{\textrm{c}_0})t + z_{\textrm{c}_0} \end{bmatrix}^{\top }\text{where,}\,\boldsymbol{\theta }=\begin{bmatrix} a_{10} & a_{20} & a_{30}\\ a_{11} & a_{21} & a_{31} \end{bmatrix}^{\top }. \end{align}

The deviations of the centre curve from the centres of $s_{i}$ may be quantified as the sum of squared distances between the centres of the LSFSs, $s_{i}$ , and the corresponding uniformly sampled points on the centre curve, ${\boldsymbol{p}}_{\textrm{c}}(t)$ :

(39) \begin{align} &\|{{\boldsymbol{p}}_{\textrm{c}_i}-{\boldsymbol{p}}_{\textrm{c}}(i/n)}\|^2=(x_{\textrm{c}_i}-x_{\textrm{c}}(i/n))^2+(y_{\textrm{c}_i}-y_{\textrm{c}}(i/n))^2+(z_{\textrm{c}_i}-z_{\textrm{c}}(i/n))^2,\,i\in (1,2,\dots, n). \end{align}

The objective is to fit the curve ${\boldsymbol{p}}_{\textrm{c}}(t)$ , such that its deviation from the centres of $s_{i}$ is minimised:

(40) \begin{align} &\underset{\boldsymbol{\theta }}{\text{Minimise}}\,\sum _{i=0}^{n}\|{{\boldsymbol{p}}_{\textrm{c}_i}-{\boldsymbol{p}}_{\textrm{c}}(i/n)}\|^2. \end{align}

Using Eq. (39) the above minimisation problem may be solved by formulating three separate minimisation problems as shown in Eqs. (41)–(43), since the three polynomials $x_{\textrm{c}}(t),y_{\textrm{c}}(t)$ and $z_{\textrm{c}}(t)$ are independent, that is, do not necessarily share a common coefficient.

(41) \begin{align}\underset{[a_{10},a_{11}]}{\text{Minimise}}\,\sum _{i=0}^{n}(x_{\textrm{c}_i}-x_{\textrm{c}}(i/n))^2, \end{align}

(42) \begin{align}\underset{[a_{20},a_{21}]}{\text{Minimise}}\,\sum _{i=0}^{n}(y_{\textrm{c}_i}-y_{\textrm{c}}(i/n))^2, \end{align}

(43) \begin{align}\underset{[a_{30},a_{31}]}{\text{Minimise}}\,\sum _{i=0}^{n}(z_{\textrm{c}_i}-z_{\textrm{c}}(i/n))^2. \end{align}

Once the centres of the LSFSs are obtained, the above problems may be seen as standard OLS problems, the solutions to which are known to be unique. Hence, the details of the solution are not presented here.

The cubic curve obtained by substituting $\boldsymbol{\theta }$ in Eq. (38) with the solutions of Eqs. (41)–(43) is shown in both Figure 6(a) and (b) as a dot-dashed curve. The centre curve, ${{\boldsymbol{p}}_{\textrm{c}}(t)},\,t\in [0,1]$ , presents itself as a safe-path from ${\boldsymbol{p}}_{\textrm{c}_0}$ to ${\boldsymbol{p}}_{\textrm{c}_n}$ , as it is free of gain-type singularity.

3.5.2. The description of the singularity-free tube (SFT)

Since the centre curve ${\boldsymbol{p}}_{\textrm{c}}(t)$ may not pass through all the centres ${\boldsymbol{p}}_{\textrm{c}_i}$ , the radius of the safe-sphere at each point of ${\boldsymbol{p}}_{\textrm{c}}(t)$ is not the same as $r_i$ . Taking this into consideration, the radius of a sphere centred at ${\boldsymbol{p}}_{\textrm{c}}(i/n)$ may be computed using Eq. (44), such that the resulting sphere lies inside the corresponding $s_{i}$ , thus being singularity-free. The radius of the sphere centred at ${\boldsymbol{p}}_{\textrm{c}}(i/n)$ , which lies inside $s_{i}$ , is denoted by $r_{\textrm{c}_i}$ and hence this sphere singularity-free. As shown in Figure 6(a), the radius $r_{\textrm{c}_i}$ is computed as

(44) \begin{equation} r_{\textrm{c}_i}=r_i-\|{{\boldsymbol{p}}_{\textrm{c}}(i/n)-{\boldsymbol{p}}_{\textrm{c}_i}}\|,\,i\in [1,2,\dots, n-1]. \end{equation}

The computation of $r_{\textrm{c}_i}$ in the above manner ensures that the sphere centred at ${\boldsymbol{p}}_{\textrm{c}}(i/n)$ with the radius $r_{\textrm{c}_i}$ is the largest sphere tangential to $s_{i}$ while being inside $s_{i}$ . Once the set of radii $r_{\textrm{c}_i}$ are obtained for $i=0,1,\dots, n$ , mono-parametric polynomial description for the radii is to be obtained. For this purpose, a quintic polynomial⁷ is used, along with the constraints $r_{\textrm{c}}(0)=r_{\textrm{c}_0}$ and $r_{\textrm{c}}(1)=r_{\textrm{c}_{\text{n}}}$ , leading to:

(45) \begin{align} r_{\textrm{c}}(t) &= a_{45} t^5 + a_{44} t^4 + a_{43} t^3 + a_{42} t^2 + a_{41} t + a_{40}, \end{align}

(46) \begin{align} &=a_{45} t^5 + a_{44} t^4 + a_{43} t^3 + a_{42} t^2 +(r_{\textrm{c}_{\text{n}}} - a_{45} - a_{44} - a_{43} - a_{42}-r_{\textrm{c}_0})t+r_{\textrm{c}_0}, \end{align}

A constrained optimisation problem (defined in Eqs. (47) and (48)) is formulated to obtain the estimates of the coefficients of $r_{\textrm{c}}(t)$ , denoted by $\boldsymbol{\theta }_{\text{r}}= \left [a_{45},a_{44},a_{42},a_{41}\right ]^{\top }$ :

(47) \begin{align} \underset{\boldsymbol{\theta }_{\text{r}}}{\text{Minimise}} \, \, &\sum _{i=0}^{n}\|{r_{\textrm{c}_i}-r_{\textrm{c}}(i/n)}\|^2, \end{align}

(48) \begin{align} \text{subject to} \, \, &r_{\textrm{c}}(i/n)-r_{\textrm{c}_i}\lt 0,\,\text{and}\,-r_{\textrm{c}}(i/n)\lt 0,\,\forall \,i\in [1,2,\dots, n-1]. \end{align}

The constraint in the above optimisation problem ensures that the radii corresponding to all the points on the centre curve are positive and less than $r_{\textrm{c}_i}$ , which would, in turn, ensure that the SFT does not intersect $\mathcal{S}$ . Due to the constraint in Eq. (48) acting at discrete points one might question if the value of radii between two discrete points $(i/n)$ and $((i+1)/n)$ , is greater than $\texttt{max}(r_{\textrm{c}}(i/n),r_{\textrm{c}}((i+1)/n))$ , that is, the possibility of an inflection occurring in between two discrete points $(i/n)$ and $((i+1)/n)$ . A systematic explanation, as to how this problem is avoided is provided in Appendix D.

Since the objective function is a sum of squares, it is convex. Using this knowledge, the logarithmic barrier function, as described in ref. [26], pp. $562\text{-}564$ , is used to convert the constrained optimisation problem given in Eqs. (47) and (48) to the following unconstrained optimisation problem:

(49) \begin{equation} \underset{\boldsymbol{\theta }}{\text{Minimise}} \, \, \frac{n}{\epsilon }\sum _{i=0}^{n}\|{r_{\textrm{c}_i}-r_{\textrm{c}}(i/n)}\|^2 + \sum _{i=0}^{n} -\text{log}(r_{\textrm{c}_i}-r_{\textrm{c}}(i/n))+\sum _{i=0}^{n} -\text{log}(r_{\textrm{c}}(i/n)),\,\quad \,\text{where}\,\epsilon =20. \end{equation}

Equation (49) now represents an unconstrained optimisation problem with a convex objective function. Subsequently, the barrier method ⁸ as described in ref. [26], pp. $568\text{-}569$ , is used to obtain the optimal estimates of $\boldsymbol{\theta }_{\text{r}}$ .

A sphere of radius $r_i$ may be defined at each point on the centre curve, ${\boldsymbol{p}}_{\textrm{c}}(t)$ , which is guaranteed to be free of gain-type singularities. Consequently, the SFT may be represented by a mono-parametric family of spheres as shown in Eq. (50):

(50) \begin{equation} {\mathcal{T}}\;:\!=\;{\boldsymbol{T}}(x,y,z,t)=({\boldsymbol{x}}-{\boldsymbol{p}}_{\textrm{c}}(t))\cdot ({\boldsymbol{x}}-{\boldsymbol{p}}_{\textrm{c}}(t))-r^2_{\textrm{c}}(t)=0,\,\ \,\text{where,}\,{\boldsymbol{x}}=[x,y,z]^{\top }\,\text{and}\,t\in (0,1). \end{equation}

The SFT, namely $\mathcal{T}$ , is described as a volume embedded in $\mathbb{R}^3$ . The interior of $\mathcal{T}$ , which may be denoted by $\texttt{int}({\mathcal{T}})$ , consists of all the points lying inside the mono-parametric family of spheres described by Eq. (50), which also contains the source and destination points. The same has been depicted visually in Figure 6(b).

Developing a closed-form parametric description of the SFT is the main contribution of this work. The SFT would act as an enabler for various path-planning algorithms; for example, it would reduce the domain to be scanned in the case of exploration-based path-planning algorithms. Thus, how the knowledge of the SFT is used may vary based on the application and the desired objective. However, only for the sake of completeness, an example of an optimisation-based gain-type singularity-free path-planning approach using the knowledge of the SFT is presented in the following section.

3.6. Planning a singularity-free path using the SFT in the constant-orientation workspace of the SRSPM

3.6.1. Example 1

Similar to the centre curve, a mono-parametric quintic polynomial form, as shown in Eqs. (51)–(53), is chosen for the path connecting ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ :

(51) \begin{align} &x_{\textrm{p}}(t)=b_{10} t^5 + b_{11} t^4 + b_{12} t^3 + b_{13} t^4 + b_{14} t + b_{15}, \end{align}

(52) \begin{align} &y_{\textrm{p}}(t)=b_{20} t^5 + b_{21} t^4 + b_{22} t^3 + b_{23} t^4 + b_{24} t + b_{25}, \end{align}

(53) \begin{align} &z_{\textrm{p}}(t)=b_{30} t^5 + b_{31} t^4 + b_{32} t^3 + b_{33} t^4 + b_{34} t + b_{35}, \ \text{where,}\,b_{ij} \in{\mathbb{R}}, \ i = 1,2,3 \ \text{and} \ j=1,2,..,5. \end{align}

An optimisation problem is formulated to estimate the values of the unknown coefficients, such that the deviation of the obtained path from the directrix is minimised. The formulation and solution of this optimisation problem are explained in the following.

The polynomials in Eqs. (51)–(53) are subjected to two constraints, namely, ${{\boldsymbol{p}}_{\textrm{p}}(0)}={{\boldsymbol{p}}_{\textrm{s}}}$ and ${{\boldsymbol{p}}_{\textrm{p}}(1)}={{\boldsymbol{p}}_{\textrm{d}}}$ . On absorbing the constraints, the number of unknown coefficients is reduced by two (similarly as Eqs. (35)–(37)), and the polynomials may thus be written as

(54) \begin{align} &x_{\textrm{p}}(t) = b_{10} t^5 + b_{11} t^4 + b_{12} t^3 + b_{13} t^2 + (x_{\textrm{d}} - b_{10}-b_{11}-b_{12}-b_{13}-x_{\textrm{c}_0})t + x_{\textrm{s}}, \end{align}

(55) \begin{align} &y_{\textrm{p}}(t) = b_{20} t^5 + b_{21} t^4 + b_{22} t^3 + b_{23} t^2 + (y_{\textrm{d}} - b_{20}-b_{21}-b_{22}-b_{23}-y_{\textrm{c}_0})t + y_{\textrm{s}}, \end{align}

(56) \begin{align} &z_{\textrm{p}}(t) = b_{30} t^5 + b_{31} t^4 + b_{32} t^3 + b_{33} t^2 + (z_{\textrm{d}} - b_{30}-b_{31}-b_{32}-b_{33}-z_{\textrm{c}_0})t + z_{\textrm{s}}. \end{align}

The coefficients to be estimated may be organised in a vector given by $\boldsymbol{\theta }_{\textrm{p}}$ :

(57) \begin{equation} \boldsymbol{\theta }_{\textrm{p}}=\left [b_{10},b_{11},b_{12},b_{13},b_{20},b_{21},b_{22},b_{23},b_{30},b_{31},b_{32},b_{33}\right ]^{\top }. \end{equation}

For estimating the coefficients $\boldsymbol{\theta }_{\textrm{p}}$ , the centre curve ${\boldsymbol{p}}_{\textrm{c}}(t)$ is sampled uniformly at $m$ points, and a quadratically constrained quadratic objective problem is formulated as

(58) \begin{align} \underset{\boldsymbol{\theta }_{\textrm{p}}}{\text{Minimise}} \, \, &\sum _{i=0}^{m}\|{{{\boldsymbol{p}}_{\textrm{p}}(i/m)}-{\boldsymbol{p}}_{\textrm{c}}(i/m)}\|^2. \end{align}

(59) \begin{align} \text{subject to} \, \, &\|{{{\boldsymbol{p}}_{\textrm{p}}(i/m)}-{\boldsymbol{p}}_{\textrm{c}}(i/m)}\|^2-r_{\textrm{c}}^2(i/m)\lt 0,\,\forall \,i\in [1,2,\dots, m-1]. \end{align}

The optimisation problem in Eqs. (58) and (59) is solved using the QuadProg++ library in C++. The optimisation problem requires that the deviation from the centre curve, ${\boldsymbol{p}}_{\textrm{c}}(t)$ , is minimised, while ensuring that the obtained path, ${\boldsymbol{p}}_{\textrm{p}}(t)$ , lies inside $\mathcal{T}$ , and is thus free of any gain-type singularities.

3.6.2. Example 2.

In Section 3.6.1, the obtained path is ensured to lie inside the SFT, $\mathcal{T}$ , while being closer to the centre curve, ${\boldsymbol{p}}_{\textrm{c}}(t)$ . In this section, the objective is to obtain an optimal path connecting ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ which also lies inside $\mathcal{T}$ . Similar to Section 3.6.1, coefficients of a mono-parametric quintic polynomial path, ${\boldsymbol{p}}_{\text{o}}(t)\;:\!=\;[x_{\text{o}}(t),y_{\text{o}}(t),z_{\text{o}}(t)]$ are estimated by formulating an constrained optimisation problem. The polynomials, $x_{\text{o}}(t),y_{\text{o}}(t)$ and $z_{\text{o}}(t)$ , subject to two constraints, ${\boldsymbol{p}}_{\text{o}}(0)={{\boldsymbol{p}}_{\textrm{s}}}$ and ${\boldsymbol{p}}_{\text{o}}(1)={{\boldsymbol{p}}_{\textrm{s}}}$ , may be written as:

(60) \begin{align} &x_{\textrm{o}}(t) = u_{10} t^5 + u_{11} t^4 + u_{12} t^3 + u_{13} t^2 + (x_{\textrm{d}} - u_{10}-u_{11}-u_{12}-u_{13}-x_{\textrm{c}_0})t + x_{\textrm{s}}, \end{align}

(61) \begin{align} &y_{\textrm{o}}(t) = u_{20} t^5 + u_{21} t^4 + u_{22} t^3 + u_{23} t^2 + (y_{\textrm{d}} - u_{20}-u_{21}-u_{22}-u_{23}-y_{\textrm{c}_0})t + y_{\textrm{s}}, \end{align}

(62) \begin{align} &z_{\textrm{o}}(t) = u_{30} t^5 + u_{31} t^4 + u_{32} t^3 + u_{33} t^2 + (z_{\textrm{d}} - u_{30}-u_{31}-u_{32}-u_{33}-z_{\textrm{c}_0})t + z_{\textrm{s}}. \end{align}

In Eqs. (60)–(62), $u_{ij}\in{\mathbb{R}},\,i=1,2,3$ and $j=1,2,\dots, 5$ . These are organised in a vector as

(63) \begin{equation} \boldsymbol{\theta }_{\textrm{o}}=\left [u_{10},u_{11},u_{12},u_{13},u_{20},u_{21},u_{22},u_{23},u_{30},u_{31},u_{32},u_{33}\right ]^{\top }. \end{equation}

The following optimisation problem is formulated in order to estimate the optimal values of $\boldsymbol{\theta }_{\textrm{o}}$ , s.t. the length of the path is minimised while ensuring that the path lies inside $\mathcal{T}$ :

(64) \begin{align} \underset{\boldsymbol{\theta }_{\textrm{o}}}{\text{Minimise}} \, \, &\int _{0}^{1} \left (\left ({\frac{\textrm{d}{x_{\text{o}}}}{\textrm{d} t}}\right )^2+\left ({\frac{\textrm{d}{y_{\text{o}}}}{\textrm{d} t}}\right )^2+\left ({\frac{\textrm{d}{z_{\text{o}}}}{\textrm{d} t}}\right )^2\right )^{1/2} \text{d}t. \end{align}

(65) \begin{align} \text{subject to} \, \, &\|{{\boldsymbol{p}}_{\text{o}}(i/m)-{\boldsymbol{p}}_{\textrm{c}}(i/m)}\|^2-r_{\textrm{c}}^2(i/m)\lt 0,\,\forall \,i\in [1,2,\dots, m-1]. \end{align}

The NMinimise function in Mathematica may be used to solve the above optimisation problem and thus an optimal path from ${\boldsymbol{p}}_{\textrm{s}}$ to ${\boldsymbol{p}}_{\textrm{d}}$ , which lies inside $\mathcal{T}$ and is of minimal length.

It may be noted at this point that the procedure presented in this section merely serves the purpose of demonstrating that singularity-free paths may be found using the concept of SFT. The method for planning these paths is not novel, and no claim is made as regards either the superiority of the methods or the results thereof over other possibilities.

4. Numerical examples

The theoretical developments presented in Section 3 are illustrated via numerical examples in this section. The specific SRSPM chosen for this purpose is the MISTRAL 800 Hexapod, the architecture parameters of which are obtained from the website⁹ of its manufacturer, Symétrie, France. The radius of the circumcircle of the base platform is $r_{\textrm{b}} = 991.4$ mm. All the linear dimensions are scaled by $r_{\textrm{b}}$ , thus rendering them dimensionless, while all the angles are measured in radians. These are listed in Table I. An example of path planning corresponding to a fixed orientation of the MP is described in this section, demonstrating the methods and algorithms presented in the previous section. All the computations have been performed on a desktop computer with $128$ GB of RAM and a Ryzen $9$ $7950$ X CPU running at $4.5$ GHz, using a single core of the CPU.

Table I. Architecture parameters of MISTRAL 800 Hexapod after scaling the lengths with respect to $r_{\textrm{b}}$ (= $991.4$ mm).

4.1. Numerical example demonstrating the utility of SFT

In this example, the set of Rodrigues parameters that define the orientation of the MP are chosen to be: $c_1 = -0.9448$ , $c_2 = 0.0703$ , $c_3 = -0.4608$ . The corresponding $\mathcal{S}$ is a regular surface (as it is verified that the singularity condition, in the form of Eq. (86), is not satisfied at any point on $\mathcal{S}$ ). The singularity surface $\mathcal{S}$ and the workspace $\mathcal{W}$ are shown in Figure 3(a). The starting point, ${{\boldsymbol{p}}_{\textrm{s}}} = [-0.2000, -0.5000, 0.5000]^{\top }$ and the destination, ${{\boldsymbol{p}}_{\textrm{d}}} = [-0.0400, -0.2200, -0.3900]^{\top }$ , are both chosen arbitrarily, but inside $\mathcal{W}$ .

Table II. Values of the tolerances and step size used during numerical computations in C/C++, the other parameters being set to the default values specified in GSL [28].

• • Following the first step in Section 3.1, the existence of a singularity-free linear path connecting ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ is queried. As can be seen in Figure 3(b), the query fails as the line joining ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ intersects $\mathcal{S}$ .

• • Next, the SFSs centred at ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ are computed (as shown in Figure 4(a)) using the generalised eigenproblem approach presented in ref. [20]. The total number of finite solutions for the spheres centred at ${\boldsymbol{p}}_{\textrm{s}}$ or ${\boldsymbol{p}}_{\textrm{d}}$ is $22$ , which is in agreement with ref. [20]. Among these, only four solutions are found to be real for both ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ . Finally, the projections of ${\boldsymbol{p}}_{\textrm{s}}$ and ${\boldsymbol{p}}_{\textrm{d}}$ on to $\mathcal{S}$ are obtained using the least radius criterion, which are computed as ${{\boldsymbol{q}}_{\textrm{s}}} = [-0.3177, -0.6832, 0.2522]^{\top }$ , and ${{\boldsymbol{q}}_{\textrm{d}}} = [-0.0489, -0.2223, -0.3555]^{\top }$ , respectively.

• • Subsequently, the geodesic curve between ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$ , that is, $\mathcal{G}$ , is computed by solving the system of ordinary differential equations presented as Eq. (12), using the boundary conditions ${\boldsymbol{q}}(0) ={{\boldsymbol{q}}_{\textrm{s}}}$ and ${\boldsymbol{q}}(1) ={{\boldsymbol{q}}_{\textrm{d}}}$ . The shooting method is used to solve this boundary value problem (BVP), which is based on a scheme of iteration in which Eq. (12) are solved initially by guessing the tangent direction at ${\boldsymbol{q}}_{\textrm{s}}$ , say, ${\boldsymbol{v}}'_{\textrm{s}}$ , and subsequently improving ${\boldsymbol{v}}'_{\textrm{s}}$ via successive iterations until Eq. (12) yields the solution for ${\boldsymbol{q}}(\beta )$ that terminates sufficiently close to ${\boldsymbol{q}}_{\textrm{d}}$ . Newton’s iterative algorithm for root-finding is used to compute ${\boldsymbol{v}}'_{\textrm{s}}$ , such that the absolute error in reaching the destination, quantified as $\|{{\boldsymbol{q}}(1) -{{\boldsymbol{q}}_{\textrm{d}}}}\|$ , falls below a chosen tolerance of $\epsilon _{\textrm{z}}$ (see Table II, for the numerical values).

• • In each iteration, Eq. (12) is solved with ${\boldsymbol{q}}(0) ={{\boldsymbol{q}}_{\textrm{s}}}$ , $\left ({\frac{\textrm{d} {\boldsymbol{q}}}{\textrm{d} \beta }}\right )_{\beta =0}\,=\,{{\boldsymbol{v}}'_{\textrm{s}}}$ as the initial conditions using the explicit Runge–Kutta–Fehlberg $(4, 5)$ algorithm (see, e.g., ref. [27], pp. 458–474) with the initial step size and tolerances noted in Table II. This numerical scheme is implemented in +C/C+++ using the routines gsl_multiroot_fsolver_iterate and gsl_odeiv2_driver_apply for the said Newton’s iteration and Runge–Kutta algorithm, respectively, both available in the GNU Scientific Library (+GSL+) ref. [28]. Newton’s method is seen to converge in $30$ iterations in this case. The corresponding initial condition is given by ${{\boldsymbol{v}}'_{\textrm{s}}} = [0.3840, 0.5470,-0.5870]^{\top }$ . The obtained geodesic curve, $\mathcal{G}$ , is plotted in Figure 4(b). By design, for “good” initial guesses, this method is assured of attaining the solution for ${\boldsymbol{q}}(\beta )$ that satisfies Eq. (12) within the tolerances specified by $\epsilon _{\textrm{rel}}$ and $\epsilon _{\textrm{abs}}$ . It now remains to verify if indeed the points on $\mathcal{G}$ lie on $\mathcal{S}$ , that is, satisfy Eq. (1). For this purpose, $n = 250$ points are sampled uniformly from $\mathcal{G}$ , namely, ${{\boldsymbol{q}}_{i}},\,i=1,2,\dots, n$ , as shown in Figure 4(b) and are substituted into Eq. (1) to compute the residues and their maximum value, respectively, as defined below:

  (66) \begin{align} &e_{\textrm{r}_{\textrm{max}}} \;:\!=\; \max _{i}\{|f({{\boldsymbol{q}}_{i}})|\}, \quad i = 1,2,\dots, n. \end{align}
  The obtained geodesic curve $\mathcal{G}$ is deemed accurate if $e_{\textrm{r}_{\textrm{max}}}\lt{\epsilon _{\textrm{z}}}$ . For the computed ${\mathcal{G}},\,e_{\textrm{r}_{\textrm{max}}} = 2.5796\times 10^{-12} \lt{\epsilon _{\textrm{z}}}$ (see Table II) and hence, it is verified to lie on $\mathcal{S}$ . An animation, “geodesic_Animation.gif” showing the normal to $\mathcal{S}$ and the Frenet-Serret local frame at each point on the geodesic, $\mathcal{G}$ is attached herewith. In this animation, the alignment of the normal to $\mathcal{G}$ and the normal to $\mathcal{S}$ at each point on $\mathcal{G}$ may be seen.

• • The LSFSs $s_{0}$ and $s_{n}$ , which are tangential to $\mathcal{S}$ at ${\boldsymbol{q}}_{\textrm{s}}$ and ${\boldsymbol{q}}_{\textrm{d}}$ respectively, are computed next using the procedure presented in Section 3.4. The generalised eigenproblem in Eq. (76) is solved using the Eigen library [29] implemented in C/C++ with the same tolerance $\epsilon _{\textrm{z}}$ as above. Further, the iterative method shown in Method 1 is used to obtain the LSFSs $s_{i}$ , which are tangent to ${{\boldsymbol{q}}_{i}},\,i=1,2,\dots, n-1$ . For the sake of clarity, LSFS touching $\mathcal{S}$ at ${{\boldsymbol{q}}_{0}}={{\boldsymbol{q}}_{\textrm{s}}}$ and ${{\boldsymbol{q}}_{n}}={{\boldsymbol{q}}_{\textrm{d}}}$ , that is, $s_{0}$ and $s_{1}$ , respectively, are shown in Figure 10(b).

• • Subsequently, the SFT, $\mathcal{T}$ , is shown in Figure 10(c), by plotting a few of the spheres among the mono-parametric family of spheres. Figure 9 shows the continuous polynomial $r_{\textrm{c}}(t)$ satisfying the constraints given in Eq. (48). The polynomials $x_{\textrm{c}}(t),y_{\textrm{c}}(t),z_{\textrm{c}}(t)$ are cubic and $r_{\textrm{c}}(t)$ is quintic in this case, and they provide a continuous, analytical description of the SFT, (see Eq. (50)).

• • Using the knowledge of $\mathcal{T}$ , a mono-parametric quintic polynomial path (shown in Figure 10(d)), which lies inside the SFT, is obtained using the procedure explained in Section 3.6.1. The time taken to obtain ${\boldsymbol{p}}_{\textrm{p}}(t)$ with $m=300$ is about $0.8$ - $1.5$ ms using the QuadProg++ library in C++.

• • Alternatively, an optimal path, ${\boldsymbol{p}}_{\text{o}}(t)$ (also shown in Figure 10(d)), lying inside the SFT, $\mathcal{T}$ and which has minimal length, is obtained using the procedure described in Section 3.6.2. Thus, demonstrating the utility of the obtained SFT in obtaining constant orientation singularity-free optimal path for SRSPM. The external solver, NMinimize in Mathematica takes about $25$ s to obtain ${\boldsymbol{p}}_{\text{o}}(t)$ with $m=300$ . In order to validate if ${\boldsymbol{p}}_{\text{o}}(t)$ indeed lies inside the SFT, $\mathcal{T}$ , the shortest distance between the path and the boundary of the SFT is computed, which is found to be $0.06\, {\textrm {mm}}$ .

The step-by-step procedure followed in the construction of $\mathcal{T}$ and consequently obtaining the singularity-free path ${\boldsymbol{p}}_{\textrm{p}}(t)$ is shown in Figure 10.

Figure 9. The discrete values of radii, $r_{\textrm{c}_i},\,i=1,2,\dots, n$ (computed using Eq. (44)) and the corresponding continuous variation of $r_{\textrm{c}}(t)$ . As can be seen in the figure, at no point $r_{\textrm{c}_i}$ exceeds $r_{\textrm{c}}(t)$ .

Figure 10. Visualisation of the directrix, the SFT $\mathcal{T}$ , and the obtained path ${\boldsymbol{p}}_{\textrm{p}}(t)$ , computed for Example 1 in Section 4.1.

5. Discussions

The mathematical developments presented in this paper and the results thereof are revisited in this section with the intention of summarising the main contributions of this work and its potential applications and extensions.

The primary contribution of the paper is the introduction of the concept of SFT, specifically, in the context of parallel manipulators. Determination of the SFT helps significantly in the process of path-planning, as the SFT encloses not one but an infinite collection of paths between a given pair of points, all guaranteed to be free of singularities. Furthermore, since the SFT has two spheres as its end-caps, even if the source and destination points are varied within these spheres, any path joining them and remaining within the original SFT would still be free of singularities without necessitating the re-computation of the SFT. Thus, the SFT provides the means of efficiently studying parametric variations of the path by changing the end-points, if necessary from an operational perspective.

Several steps in computation of the SFT are entirely novel, to the extent of the knowledge of the authors. The first among these is the concept of computing the geodesic curve on the singularity surface. In general, in the existing works, the knowledge of the singularity surface, even when it is available, is used only to the extent of avoiding it, and therefore, such knowledge is not utilised to realise its full potential in application. In this paper, however, it is the geodesic curve on the singularity surface that finds a smooth path between the projections of the source and the destination points on the said surface, thus providing valuable guidance to the process of computing the actual path. It is a non-trivial task to find such a smooth connector via other methods, whether analytical or randomised. The geodesic curve has been computed using the implicit equation of the singularity surface, which obviates the need for a parametrisation of the same. This is likely to broaden the scope of application of such methods to other manipulators, wherein such singularity manifolds may have been computed but not explicitly parametrised. Furthermore, this formulation circumvents the issues of potential singularities in such parametrisations, even when it is present.

Since the SFT is derived as a mono-parametric family of spheres analytically, it is trivial and yet deterministic to ascertain if any given in space point falls inside it or not. Thus, any number of algorithms, irrespective of their nature, that is, analytical or randomised, may be used to plan singularity-free paths between these points just by remaining inside the SFT and using additional constraints and objectives, as desired. For instance:

• • Analytical algorithms for path planning would benefit from the parametric and closed-form description of the SFT, as it can be used directly for a constraint function (as exemplified in Section 3.6).

• • Algorithms involving randomised scans only need to scan inside the SFT as opposed to the entire workspace, which reduces the computational cost of employing such algorithms.

• • Since the centre curve of the SFT is always available as a feasible solution, its can used to bias the probabilistic algorithms such as RRT and RRT $^*$ to achieve their solutions faster.

• • The concept of SFT along with additional constraints may also be used in defining the risk/reward functions for reinforcement learning (RL)-based path planning/control of SPMs and other parallel manipulators.

Therefore, it may be stated that while the computation of the SFT is a mathematical exercise involving several concepts and tools, once derived, it can be used and reused by all forms of planning algorithms, without incurring the burden of having to consider singularities in the process.

It may be argued that gain-type singularities form only one of several considerations while planning the path of a parallel manipulator. Indeed, as described in ref. [30], there are four distinct issues to be considered, in general: workspace constraints of loss-type singularities,¹⁰ gain-type singularities, interference among the links, and physical limits on the permissible motions at the joints. It can also be argued that the problem of gain-type singularities is the hardest among these, since it involves demanding analytical computations involving the singularity surface. Once this problem is addressed satisfactorily, the others can be tackled incrementally by imposing additional constraints while planning a particular path inside the SFT. It also affords significant flexibility in choosing the methods to address these issues, whether analytical or numerical. For instance, if one uses a sampling based method to determine if a point on the path is free of link interference, the knowledge of the SFT helps in the process by delimiting the search space to its volume, as opposed to that of the workspace. Furthermore, the gain-type singularities and hence the SFT depend solely on the architecture of the manipulator, that is, its kinematic aspects, such as link lengths and the geometry of the platforms, as opposed to other physical details, such as the girths or lateral dimensions of the links which impact the last two considerations. Thus, the SFT introduces some modularity in the process, wherein one can consider the effects of kinematic and non-kinematic details of the manipulator in a decoupled manner.

Another important feature of this work is the computational efficiency. The computation of the SFT takes about $103$ ms (see Table III) using a single core of a contemporary CPU. At the level of applications, such an ease of computation will allow the analyst to evaluate multiple options, such as evaluating multiple fixed orientations for feasibility for computing a path between a given pair of points, or to choose one of these, based on some desired criteria. Such computational capabilities are not available for spatial parallel manipulators, as far as authors’ knowledge goes.

Table III. Time taken for computation in C/C++ in each step of obtaining the SFT. The average timing over $100$ trials is reported below.

Table IV. Degrees and sizes of the polynomials in Eqs. (25)–(27), Eqs. (68), (70)–(72).

In the future, this work will be extended to path planning in $\mathbb{SE}(3)$ , wherein the MP will be allowed to modify its orientation along the path.

6. Conclusions

This paper introduces the concept of a singularity-free tube, which is an important enabler for singularity-free path planning of spatial parallel manipulators. It incorporates several novel features, such as the use of the concept of the geodesic curves on the singularity surface itself for the initial task of exploration. The SFT also allows for a choice among an infinite number of possible smooth singularity-free paths from the source to the destination. These novel analytical developments are complemented by numerical optimisation schemes, which complete the solution process by obtaining a smooth mono-parametric analytical description of paths from the source to the destination. The methods presented are robust as well as computationally inexpensive, and as such, they may be suitable for a large number of applications in different scenarios as well as in other manipulators.