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. [Reference Shende, Patra, Prasad, Bandyopadhyay, Lenarčič and Siciliano17]) 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. [Reference Bandyopadhyay and Ghosal9]. 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.1. Steps followed in finding a singularity-free tube, $\mathcal{T}$

The proposed method builds upon some of the concepts used in ref. [Reference Prasad, Bandyopadhyay, Sen, Mohan and Ananthasuresh18]. 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. [Reference Nag, Reddy, Agarwal, Bandyopadhyay, Lenarčič and Merlet19, Reference Nag and Bandyopadhyay20]. 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 projectionFootnote ³ 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. [Reference Nag, Reddy, Agarwal, Bandyopadhyay, Lenarčič and Merlet19].

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. [Reference Kühnel21], pp. 140–165, [Reference Struik22], pp. 127–136). As shown by Coste and Moussa in ref. [Reference Coste and Moussa23], 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. [Reference Struik22], 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. [Reference Hyde, Ninham, Andersson, Larsson, Landh, Blum, Lidin, Hyde, Ninham, Andersson, Larsson, Landh, Blum and Lidin24]). 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. [Reference Struik22], 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. [Reference Struik22], 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. [Reference Lebedev, Cloud and Eremeyev25], 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. [Reference Nag, Reddy, Agarwal, Bandyopadhyay, Lenarčič and Merlet19, Reference Nag and Bandyopadhyay20]. 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:Footnote ⁴

(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 Footnote ⁵ 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 Footnote ⁶ 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 polynomialFootnote ⁷ 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. [Reference Boyd and Vandenberghe26], 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 Footnote ⁸ as described in ref. [Reference Boyd and Vandenberghe26], 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.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 [Reference Gough28].

• • 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. [Reference Nag and Bandyopadhyay20]. 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. [Reference Nag and Bandyopadhyay20]. 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. [Reference Mathews and Fink27], 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. [Reference Gough28]. 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 [Reference Guennebaud and Jacob29] 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.