**Multiple-Target Homotopic Quasi-Complete Path Planning Method for Mobile Robot Using a Piecewise Linear Approach**

**Gerardo Diaz-Arango 1, Hector Vazquez-Leal 2,3,***<sup>∗</sup>* **, Luis Hernandez-Martinez 4, Victor Manuel Jimenez-Fernandez 2, Aurelio Heredia-Jimenez 5, Roberto C. Ambrosio 6, Jesus Huerta-Chua 7, Hector De Cos-Cholula <sup>4</sup> and Sergio Hernandez-Mendez <sup>1</sup>**


Received: 4 May 2020; Accepted: 27 May 2020; Published: 8 June 2020

**Abstract:** The ability to plan a multiple-target path that goes through places considered important is desirable for autonomous mobile robots that perform tasks in industrial environments. This characteristic is necessary for inspection robots that monitor the critical conditions of sectors in thermal, nuclear, and hydropower plants. This ability is also useful for applications such as service at home, victim rescue, museum guidance, land mine detection, and so forth. Multiple-target collision-free path planning is a topic that has not been very studied because of the complexity that it implies. Usually, this issue is left in second place because, commonly, it is solved by segmentation using the point-to-point strategy. Nevertheless, this approach exhibits a poor performance, in terms of path length, due to unnecessary turnings and redundant segments present in the found path. In this paper, a multiple-target method based on homotopy continuation capable to calculate a collision-free path in a single execution for complex environments is presented. This method exhibits a better performance, both in speed and efficiency, and robustness compared to the original Homotopic Path Planning Method (HPPM). Among the new schemes that improve their performance are the Double Spherical Tracking (DST), the dummy obstacle scheme, and a systematic criterion to a selection of repulsion parameter. The case studies show its effectiveness to find a solution path for office-like environments in just a few milliseconds, even if they have narrow corridors and hundreds of obstacles. Additionally, a comparison between the proposed method and sampling-based planning algorithms (SBP) with the best performance is presented. Furthermore, the results of case studies show that the proposed method exhibits a better performance than SBP algorithms for execution time, memory, and in some cases path length metrics. Finally, to validate the feasibility of the paths calculated by the proposed planner; two simulations using the pure-pursuit controlled and differential drive robot model contained in the Robotics System Toolbox of MATLAB are presented.

**Keywords:** robot motion; path planning; piecewise linear approximation; multiple-target path planning; autonomous mobile robot; homotopy based path planning

#### **1. Introduction**

Autonomous mobile robot is an entity capable of performing a wide variety of tasks that involve displacement in its workspace, such as home service, pickup and delivery assistance in offices, monitoring factories, and so forth. An autonomous robot develops its assignment without human intervention because it commonly implies a certain degree of risk for a human being or environmental conditions make the use of teleoperation impossible. To execute its task, the robot must be capable to plan a path from an initial position to a target-position. In principle, path planning is a geometric process in which an autonomous agent must find a collision-free path in its workspace, without considering its kinematics and dynamics restrictions [1–4]. Then, once a path is specified, another process is executed to calculate the motion plan using the kinodynamic properties of the robot. Usually, for a path planning process, the robot is considered as a point in the configuration space (*C*space), which is the space generated by all feasible positions that it can reach [3,5,6]. Then, *C*space is divided into free configuration space (*C*free) for valid positions and obstacles space (*C*obs) for all forbidden configurations.

#### *1.1. Planning Algorithms*

Finding a collision-free path in *C*free may seem like an easy task for a human agent, nevertheless, it is a very complex issue for artificial intelligence. This has been the main motivation for researchers and developers community which have worked in this area during the last five decades [2–4,6–11]. These efforts have generated diverse planning algorithms with particular characteristics that define its degree of completeness, approach, and configuration of problems that each planner is capable of solve. In general, the planners can be classified into three categories.


contains all probabilistic-complete algorithms such as: probabilistic road maps (PRM) [7,29], expansive spaces trees (EST) [30], rapidly-exploring random trees (RRT) [3,7,31], bidirectional RRT (Bi-RRT, also named RRT-Connect) [32], RRT combined with a shortest path approach A (RRT\*) [7], kynodinamic motion planning by interior-exterior cells of exploration (KPIECE) [1], and collision checking efficient algorithms (LazyPRM and LazyRRT) [33,34] which have the common modules; uniform distribution samples generator, collision checker, local planner, and smoothing post-processing algorithm [2,3,35]. (e) Homotopy continuation based planners are a new category that contains some variants of the original method introduced in Reference [19]. These variants have been proposed to improve performance, minimize computation time, and obtain the shortest path (reported in References [10,20,36,37]).

3. Planners are classified according to their ability to reuse pre-processed data for solving problems into: (a) single-query and (b) multiple-query. In this way, the algorithms and methods mentioned previously, as well as their variants, fall in one of these categories. On the one hand, multiple-query algorithms are commonly applied to solve static environments and the roadmap generated can be reused as many times as needed. Therefore, queries are very fast, nevertheless, the computational cost and time to generate the roadmap are impractical for dynamic environments; algorithms like PRM and graph-based (Visibility graphs and A\*) have this property. On the other hand, for single-query planners, roadmap generation and extend function are developed in parallel to reduce the high computational cost of analyzing the entire environment. This characteristic makes these algorithms faster, nevertheless, the resulting roadmap is only useful for the current query. Some algorithms contained in this category are RRT, EST, KPIECE, artificial potential fields method, HPPM, among others (see Table 1).

Table 1 shows the main characteristics of the most used collision-free path planners. The last row presents the advantages of HPPM (the method of interest in this work) which makes it one of the most versatile and reliable planners, based on its single-query execution, deterministic approach, and quasi-complete success. Furthermore, this table remarks the main issue, or bottleneck, found in each sampling-based algorithm, APF method, visibility graph approach, A\* algorithm, and HPPM. On the one hand, for SBP algorithms, the bottleneck is the computing time spent during the collision query procedure which consumes about sixty percent of the total time [9]. On the other hand, the main feature of LazzyRRT and LazzyPRM over RRT and PRM is the lower time for the collision query. Nevertheless, for these planners, the number of samples (density of nodes) has a high impact on its success to find a feasible path. For A\* algorithm, its success lies in the resolution of the discretized workspace (better resolution implies higher memory consumption and execution time). For AFP, the main issue is the local minimum, which confine the robot and prevents the method to continue. For the visibility graph approach, the bottleneck is the calculation time, which increases according to the number of obstacles. As for HPPM, the main issue lies in the selection of the repulsion parameter for obstacles. This parameter plays an important role in this method because it determines the length of the path and, in consequence, the execution time. Although, the repulsion parameter selection for HPPM is considered solved for particular applications in References [38,39] and References [20,21,36,37]. A generalized strategy, at this moment, it is currently considered an open problem. One of the main contributions in the present work is the implementation of a systematic criterion to select the repulsion parameter to generate an optimal path, in terms of length and execution time, for HPPM.


**Table 1.** Main characteristics of collision-free path planners.

#### *1.2. Multiple-Target Path Planning*

Commonly, the single-query and multiple-query collision-free planners are developed to find a valid path from a starting point to a single target point. Nevertheless, for applications like: robot vacuum cleaner, service at home, pickup and delivery services in office, industrial process and inspection, victim rescue, museum guidance, and land mine detectors, robots need to move through a sequence of target points [40–45]. The multiple-target collision-free path planning problem lacks deeper study, since it is considered a trivial task by using a point-to-point strategy. It operates through a set of sub-paths which creates a full multiple-target path, which connects each target point to another. While this strategy operates properly, its main disadvantage is the unnecessary and redundant turnings generated by connecting two or more segments; causing the sub-paths to cross each other and, in consequence, the robot will cross through the same place more than once [46] (see Figure 1b). In this sense, a complete path found through the point-to-point strategy is not optimal in terms of length.

**Figure 1.** Efficient route and unnecessary turnings for multiple-target problem. (**a**) Efficient path for multiple-target problem. (**b**) Unnecessary turnings in multiple-target path (point-to-point strategy).

Figure 1a shows an example of desirable path (the shortest path) for a multiple-target problem. It can be noticed that the path is smooth without redundant configurations and pass very close to every target. On the other hand, Figure 1 depicts an example of point-to-point strategy for multiple-target problem. The unnecessary turnings and crossings between sub-paths are noticeable, making the trajectory inefficient in terms of distance. Figure 2 shows an example of a common office-like environment, this floor plan corresponds to a museum hall and the robot performs guidance tasks. In this case, it is considered that the robot calculates a path in some milliseconds, then, people and other dynamic agents could be represented as quasi-static obstacles (circular shapes in the environment map). The tasks of autonomous robots are to guide a group of people, provide relevant information about the displayed art, and visit all rooms while avoiding any collision with objects.

**Figure 2.** Multiple-target path planning issue in a complex environment (museum floor plan).

The museum environment in Figure 2 contains three common challenges for path planning algorithms. First, the map contains many obstacles that represent an issue for planners like BUG algorithm, artificial potential field method, and visibility graph methods. For these, the complexity and computing time increases accordingly to obstacles density. Second, the map displays several narrow corridors between walls and obstacles that are difficult to solve by SBP algorithms; while for artificial potential fields method this configuration exhibits local minimum problems. Third, for the layout in this map, the path should visit fourteen places, that is, the robot system must be capable to plan a path for multiple-target points. Although the first two issues have been thoroughly studied, multiple-targets is still considered as an open problem [40,44]. In this paper, a capable planner based on homotopy continuation to deal with narrow corridors, environment maps with hundreds of obstacles, and multiple-target issues is proposed. This paper is organized as follows. In Section 2, the bases of the homotopic path planning method are explained. The spherical algorithm applied to trace homotopy curves is presented in Section 3. Chua's canonical piecewise linear model is shown in Section 4. Section 5 provides the main contributions and formulation of our proposed multiple-target HPPM. Furthermore, in Section 6, a HPPM with visibility approach is explained. Some case studies and performance comparison between the proposed methodology and SBP algorithms are presented in Section 7. Path tracking examples applying the pure-pursuit algorithm are provided in Section 8. Finally, the concluding remarks are given in Section 9.

#### **2. Homotopic Path Planning Method**

Homotopic path planning method (HPPM) is a single-query planner capable to obtain a collision-free path by using a mathematical model emanated from the configuration space, initial and final positions [10,19]. This planner is based on homotopy continuation methods (HCM) which, frequently, are employed to find multiple solutions of non-linear algebraic equation system (NAES) of the form

$$f(X) = 0: \quad \mathbb{R}^n \longrightarrow \mathbb{R}^n. \tag{1}$$

The homotopy operates through a deformation of *f*(*X*) by adding a function *G*(*X*) and a homotopy parameter *λ*, such that

$$H(X,\lambda) = \lambda f(X) + (1-\lambda)G(X) = 0,\tag{2}$$

where *<sup>H</sup>* (*f*(*X*)) <sup>=</sup> *<sup>H</sup>*(*X*, *<sup>λ</sup>*) : <sup>R</sup>*n*+<sup>1</sup> −→ <sup>R</sup>*n*, *<sup>X</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* , *<sup>λ</sup>* <sup>∈</sup> [0, 1]; *<sup>G</sup>*(*X*) is a function with a trivial or known solution. The homotopy system has the following properties


HPPM process takes the system of equations that models the configuration space and, through the application of Newton's homotopy, a free collision path is obtained. For Newton's homotopy, the auxiliary function is *G*(*X*) = *f*(*X*) − *f*(*X*0); where *X*<sup>0</sup> is the known starting point. Then, HPPM employs the spherical algorithm (SA) to properly track the homotopy curve. For this system of equations, the curve represents a sequence of points that describes a continuous motion from a starting point to a target-point avoiding collisions with obstacles [10,19]. The configuration for a 2-D environment is presented as NAES *f*1(*x*, *y*) = 0 and *f*2(*x*, *y*) = 0. For both equations the unique solution lies on the target-point (*xT*, *yT*) [10,19]. The configuration space is modeled by the following equations

$$f\_1(\mathbf{x}, \mathbf{y}) = L\_1(\mathbf{x}, \mathbf{y}) = \mathbf{0},\tag{3}$$

$$L f\_2(\mathbf{x}, y) = L\_2(\mathbf{x}, y) + \mathbf{g}(\mathbf{x}, y) = \mathbf{0},\tag{4}$$

where *g*(*x*, *y*) model the obstacles on the map, which creates singular regions around them. *L*1(*x*, *y*) and *L*2(*x*, *y*) are two auxiliary straight lines that intersect at the target-point, these are represented by

$$L\_i(\mathbf{x}, y) = -y + m\_i \mathbf{x} + b\_i = 0; \; i = 1, 2,\tag{5}$$

then, applying Newton's homotopy to Equation (3) and Equation (4), the NAES is transformed into

$$H = \begin{cases} H\_1(f\_1(\mathbf{x}, y), \lambda) = f\_1(\mathbf{x}, y) - (1 - \lambda)f\_1(\mathbf{x}\_0, y\_0) = 0, \\ H\_2(f\_2(\mathbf{x}, y), \lambda) = f\_2(\mathbf{x}, y) - (1 - \lambda)f\_2(\mathbf{x}\_0, y\_0) = 0, \end{cases} \tag{6}$$

where (*x*0, *y*0) is the initial point. The obstacles in the configuration space can be modeled by circumferences, ellipsoids, and other closed curves. For this work, only circular and quasi-rectangular shapes modeled using circle and ellipsoid equations [19], respectively are employed. Furthermore, the solid obstacles representation proposed in Reference [36] is used, such that each circular obstacle is defined by

$$|\mathcal{C}\_i(\mathbf{x}, \mathbf{y})| + \mathcal{C}\_i(\mathbf{x}, \mathbf{y}) = \mathbf{0},\tag{7}$$

where *C*(*x*, *y*) is the general equation of the circumference, modeled by

$$\mathcal{C}\_{i}(\mathbf{x}, \mathbf{y}) = (\mathbf{x} - \mathbf{x}\_{\mathbb{C}, i})^2 + (\mathbf{y} - \mathbf{y}\_{\mathbb{C}, i})^2 - r\_{\mathbb{C}, i}^2 = \mathbf{0},\tag{8}$$

*i* = 1, 2, 3, ... , *k*, *k* is the number of circular obstacles in the map; (*xC*,*i*, *yC*,*i*) represents the center of the *i*-th circular obstacle, and *rC*,*<sup>i</sup>* its radius. The solid circle obstacle representation has solution for all points inside and at the contour of the circumference, that is, <sup>∀</sup> (*x*, *<sup>y</sup>*) <sup>∈</sup> <sup>R</sup><sup>2</sup> such that (*x* − *xC*,*i*)<sup>2</sup> + (*y* − *yC*,*i*)<sup>2</sup> ≤ *rC*,*i*. Then, the expression representing all circular obstacles in the map is

$$\mathcal{W}\_{\mathbb{C}}(\mathbf{x}, \mathbf{y}) = \sum\_{i=1}^{i=k} \frac{p\_{\mathbb{C}, i}}{|\mathbb{C}\_{i}(\mathbf{x}, \mathbf{y}) + |\mathbb{C}\_{i}(\mathbf{x}, \mathbf{y})|} \, \text{} \, \tag{9}$$

here, *pC*,*<sup>i</sup>* is the repulsion parameter of the *i*-th circular obstacle [10,19]. Likewise, this principle is applied to ellipsoidal obstacles using

$$R\_{\hat{\jmath}}(\mathbf{x}, \mathbf{y}) = \left(\frac{\mathbf{x} - \mathbf{x}\_{R, \hat{\jmath}}}{n}\right)^{2\eta} + \left(\frac{\mathbf{y} - \mathbf{y}\_{R, \hat{\jmath}}}{\beta}\right)^{2\nu} - 1 = 0,\tag{10}$$

where, the shape of each rectangular obstacle is defined by *R*(*x*, *y*) in Equation (10); *j* = 1, 2, 3, ... , *l*, *l* is the number of rectangular obstacles in the map; (*xR*,*j*, *yR*,*j*) is the center for the *j*-th rectangular obstacle; *α* and *β* are width and length, respectively; *η* and *ν* are two integers that define the sharp of the vertex [19,36]. The solid representation for ellipsoidal obstacles is modeled as

$$\left| R\_{\dot{\jmath}}(\mathbf{x}, \mathbf{y}) \right| + R\_{\dot{\jmath}}(\mathbf{x}, \mathbf{y}) = \mathbf{0},\tag{11}$$

in a similar way to circular obstacles, this representation has a solution for all points inside and at the contour of the ellipsoid. Then, the expression that represents all ellipsoidal obstacles in the map is

$$\mathcal{W}\_{\mathbb{R}}(\mathbf{x}, \mathbf{y}) = \sum\_{j=1}^{j=1} \frac{p\_{\mathbb{R}, j}}{\mathcal{R}\_{j}(\mathbf{x}, \mathbf{y}) + \left| \mathcal{R}\_{j}(\mathbf{x}, \mathbf{y}) \right|} \, \text{} \tag{12}$$

*pR*,*<sup>j</sup>* is the repulsion parameter for the *j*-th ellipsoidal obstacle. It is important to notice that the points inside and at the contour of the obstacles result in divisions by zero for Equations (9) and (12). These divisions by zero creates singularities in the map, thus, HPPM uses them to automatically avoid collisions with the obstacles. The effect of all obstacles, whether they are ellipsoid or circular, are contained in the following expression

$$\mathcal{W}(\mathbf{x}, \mathbf{y}) = \mathcal{W}\_{\mathbb{C}}(\mathbf{x}, \mathbf{y}) + \mathcal{W}\_{\mathbb{R}}(\mathbf{x}, \mathbf{y}), \tag{13}$$

where *WC*(*x*, *y*) represents the solid circular obstacles and *WR*(*x*, *y*) represents the solid ellipsoidal obstacles. The term *g*(*x*, *y*) in Equation (4) is redefined as

$$\mathcal{g}(\mathbf{x}, \mathbf{y}) = \mathcal{W}(\mathbf{x}, \mathbf{y}) - \mathcal{W}(\mathbf{x}\_T, \mathbf{y}\_T), \tag{14}$$

here, *W*(*xT*, *yT*) is added to the NAES to cancel the effect of the obstacles at the target-point.

#### **3. Homotopy Path Tracking Scheme**

Spherical algorithm (SA) is a tool applied to track continuous curves resulting from a non-homogeneous NAES. Previous works have demonstrated the effectiveness of SA to trace homotopy curves *<sup>γ</sup>* [10,36,47–50]. The homotopy is an operator such that *<sup>H</sup>* (*f*(*X*)) : <sup>R</sup>*n*+<sup>1</sup> −→ <sup>R</sup>*n*, that is, homotopy system (6) is non-homogeneous. SA operates by adding an *n*-dimensional sphere equation *Si*(*X*, *λ*) (the dimension of sphere or hypersphere depends on the number of variables in the system of equations) that guarantees, at least, two cross points with the homotopy curve (6), as long as its center is at *H*−1(0) [47]. The homotopy process is described as follows: first step, the center of first sphere (*S*1) is placed at *O*<sup>1</sup> = (*x*0, *y*0, *λ* = 0) (see Figure 3). Then, the first root on the right is calculated by a predictor-corrector scheme; the root represents the intersection between the homotopy curve and the first hypersphere for values of *λ* ≥ 0. The next step consists of the assignation of the calculated root on the right as the center (*O*2) of a new hypersphere (*S*2), as shown in Figure 3. Numerical path tracking scheme is based on

$$H\_S = \begin{cases} H\_1(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = 0, \\ H\_2(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = 0, \\ S\_i(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = 0, \end{cases} \tag{15}$$

for three dimensions (Bi-dimensional space and homotopy parameter *λ*), the sphere is represented as

$$S\_i(\mathbf{x}, y, \lambda) = (\mathbf{x} - \mathbf{c}\_x)^2 + (y - \mathbf{c}\_y)^2 + (\lambda - \mathbf{c}\_\lambda)^2 - r\_s^2 = 0,\tag{16}$$

where *rs* is the radius and *Oi* = (*cx*, *cy*, *cλ*) is the center of the sphere for each step in the spherical tracking.

**Figure 3.** 2-D representation of spherical tracking algorithm process.

Homotopy curves tracking is a complex task, because the sharp turning points makes the convergence of the corrector step difficult [10,19,48,49]. In this regard, a proper predictor scheme is employed to improve the performance of the SA; in this work the Euler's predictor is applied. Moreover, this scheme has shown adequate performance for tracking homotopy curves making the corrector scheme to converge faster [10,49]. This predictor scheme operates using the tangent vector to calculate an approximation of the next point in the tracking curve, as explained in the following section.

#### *3.1. Euler's Predictor Scheme*

Euler's scheme is a complement of SA that provides a predictor point close to the intersection between homotopy curve and the respective sphere at each step of the tracking. Work in References [10,49] explain its operation and show its effectiveness as a predictor scheme in SA. For system equations with three variables, the Euler predictor point is calculated using

$$\left(\mathbf{x}\_{p\prime}y\_{p\prime}\boldsymbol{\lambda}\_{p}\right) = \left(\mathbf{x}\_{i\prime}y\_{i\prime}\boldsymbol{\lambda}\_{i}\right) + r\_{\rm s}\left(\frac{\vec{v}\_{p}}{||\vec{v}\_{p}||}\right),\tag{17}$$

where - *xp*, *yp*, *λ<sup>p</sup>* is the predictor point; (*xi*, *yi*, *λi*) is the center of the sphere; *rs* is the radius of the sphere, and *v<sup>p</sup>* is the tangent vector. Figure 4 shows a 2-D view of the predictor scheme operation in the SA. The predictor point - *xp*, *yp*, *λ<sup>p</sup>* is located at the intersection between the vector point and the sphere *Si*.

**Figure 4.** Euler's predictor scheme.

*Sensors* **2020**, *20*, 3265

The tangent vector is obtained from the partial derivatives of the homotopy system (6), as explained in References [10,49]. This method uses partial derivatives with respect to an arbitrary parameter *ρ*, such that (*x*(*ρ*), *y*(*ρ*), *λ*(*ρ*)); then, the chain rule is employed to calculate the tangent vector for the homotopy system at the point (*xi*, *yi*, *λi*). Placing the tangent vector in Equation (17) the predictor point ( - *xp*, *yp*, *λ<sup>p</sup>* ) is obtained

$$f(\mathbf{x}\_{p}, \mathbf{y}\_{p}, \boldsymbol{\lambda}\_{p}) = (\mathbf{x}\_{i}, \mathbf{y}\_{i}, \boldsymbol{\lambda}\_{i}) + r\_{s} \frac{(\mathbf{x}\_{i}'(\boldsymbol{\rho}), \mathbf{y}\_{i}'(\boldsymbol{\rho}), \boldsymbol{\lambda}\_{i}'(\boldsymbol{\rho}))}{||(\mathbf{x}\_{i}'(\boldsymbol{\rho}), \mathbf{y}\_{i}'(\boldsymbol{\rho}), \boldsymbol{\lambda}\_{i}'(\boldsymbol{\rho}))||}\tag{18}$$

where,

$$(\mathbf{x}'(\rho), \mathbf{y}'(\rho), \lambda'(\rho)) = \left(\frac{\partial \mathbf{x}\_i(\rho)}{\partial \rho}, \frac{\partial \mathbf{y}\_i(\rho)}{\partial \rho}, \frac{\partial \lambda\_i(\rho)}{\partial \rho}\right),\tag{19}$$

here, the partial derivative vector of each variable in the system with respect to *ρ*, evaluated at (*xi*, *yi*, *λi*), represents the tangent vector [10,49].

#### *3.2. Broyden's Method as Corrector Scheme*

The corrector scheme is the core of the HCM process, capable to obtain one solution for homogeneous NAES from an initial condition. Commonly, Newton-like correctors are described by the expression

$$X\_{i+1} = X\_i - \left[ f(X\_i) \right]^{-1} f(X\_i),\tag{20}$$

here, *<sup>X</sup>* <sup>∈</sup> <sup>R</sup>*n*, *<sup>i</sup>* <sup>=</sup> 0, 1, 2, ... , *<sup>n</sup>*; *Xi*+<sup>1</sup> is the next point in the iterative process; *Xi* is the current point and represents the initial condition for *i* = 0, and [*J*(*Xi*)] <sup>−</sup><sup>1</sup> is the Jacobian inverse matrix of the NAES *f*(*X*). The iterative process ends when the stop criterion is achieved or the maximum number of iterations is reached; both parameters set by the user. Broyden's method is a quasi-Newton method employed to calculate the roots of non-linear algebraic equation systems. For this, the Jacobian matrix *J*(*Xi*) is replaced by an approximation *Ai* to reduce its calculation complexity. *Ai* is calculated using

$$A\_j = A\_{j-1} + \frac{f(X\_j) - f(X\_{j-1}) - A\_{j-1}[X\_j - X\_{j-1}]}{||X\_j - X\_{j-1}||^2} [X\_j - X\_{j-1}]^T,\tag{21}$$

for *j* = 1, 2, 3, . . . , *n*; Broyden's method is represented by

$$X\_{\circ+1} = X\_{\circ} - \left[A\_{\circ}(X\_{\circ})\right]^{-1} f(X\_{\circ}),\tag{22}$$

where, in the first iteration (*j* = 1), *X*<sup>1</sup> is calculated using Equation (20) and the matrix *A*<sup>0</sup> = *J*(*X*0), that is, the Jacobian matrix only needs to be calculated once. Furthermore, the Jacobian matrix calculation is also employed to obtain the predictor point, hence, this matrix is strongly used in two complementary SA processes. In this work, Broyden's Method stop criterion is set at *i* = 20 or when the next criterion is fulfilled

$$\|f(X\_i)\| < 1 \times 10^{-9}.\tag{23}$$

For the spherical algorithm, the initial point *X*<sup>0</sup> of Broyden's method procedure is the predictor point (*xp*, *yp*, *λp*) obtained using the Euler's scheme. Then, the Broyden's process is executed to solve the NAES system (15) for the current sphere (*Si*) until the stop criterion is fulfilled. Which implies that the next solution point (*xi*+1, *yi*+1, *λi*+1) has been found. Figure 5 depicts, in 2-D, the entire operation of the predictor-corrector scheme for SA.

**Figure 5.** Corrector scheme for spherical tracking.

#### **4. Canonical Piecewise Linear Representation**

A piecewise linear model can be defined as a mathematical representation that collects linear, related, descriptions to approximate a nonlinear function. The main reason that motivates the use of this type of model is the simplicity of their structure which allows an efficient implementation in algorithms. Although there are many piecewise linear models reported in the literature [51–55], due to its compact formulation, reduced number of parameters, and low computational requirements, the most popular is the so-called Chua's model; it is described by a compact global representation named canonical piecewise linear function, given by

$$y(\mathbf{x}) = a + b\mathbf{x} + \sum\_{i=1}^{\sigma} c\_i |\mathbf{x} - \mathbf{x}\_{B,i}| = \mathbf{0},\tag{24}$$

where *σ* is the number of breakpoints. Model parameters *a*, *b*, and *ci*, for *i* = 1, 2, ..., *σ*, can be determined as follows

$$a = y(0) - \sum\_{i=1}^{\sigma} c\_i |x\_{B,i}|\_{\sigma}$$

$$b = \frac{\mathcal{J}\_1 + \mathcal{J}\_{\sigma+1}}{2},$$

$$c\_i = \frac{\mathcal{J}\_{i+1} - \mathcal{J}\_i}{2}.$$

These parameters are strongly related to the graph of the piecewise linear function *y*(*x*). For instance, *b* and *ci* are both described in terms of J*<sup>i</sup>* which represents the slope of the *i*-th constitutive linear segment in the graph of *y*(*x*). Moreover, the parameter *a* is computed considering the *σ*-breakpoints (*xB*,*i*, *yB*,*i*), for *i* = 1, 2, ..., *σ* included in the graph of *y*(*x*). To graphically illustrate this relation, Figure 6 shows a single-valued piecewise linear function constituted by *σ* breakpoints and *σ* + 1 segments.

**Figure 6.** Single-valued piecewise linear function.

In this work, Chua's model serves as the base for the proposed homotopy scheme and the two-dimensional PWL model is described by

$$\text{PVIL}(\mathbf{x}, \mathbf{y}) = \mathbf{y} - \left( a + b\mathbf{x} + \sum\_{i=1}^{\sigma} c\_i |\mathbf{x} - \mathbf{x}\_{B,i}| \right), \tag{25}$$

where PWL(*x*, *y*) is an implicit representation of Equation (24). It can notice, that depends on *x* values. The breakpoints (*xB*,*i*, *yB*,*i*) can take any value in the *y*-axis but, for *x*-axis, the values must be incremental such that *xB*,*<sup>i</sup>* < *xB*,*i*+<sup>1</sup> < ... < *xB*,*i*+*n*.

#### **5. Multiple-Target Homotopic Path Planning Method**

This section is devoted to explain our proposed multiple-target homotopic path planning method (MTHPPM) and highlight the most relevant contributions of this work. MTHPPM is a planner based on the HPPM principle; for both, a set of auxiliary functions is necessary to generate a NAES in which the target-points are solutions. The multiple-target points in MTHPPM are provided by a set of two auxiliary functions that intersect in more of one point. In this work, target-points are generated by two single-valued PWL of the form Equation (24) with a similar formulation to the auxiliary functions in HPPM. These are represented by

$$f\_{\rm PWL\_1}(\mathbf{x}, y) = \rm PWL\_1(\mathbf{x}, y), \tag{26}$$

$$\begin{aligned} \left(\text{PWL}\_2(\mathbf{x}, \mathbf{y})\right) &=& \text{PWL}\_2(\mathbf{x}, \mathbf{y}) \\ &- \left(\text{g}(\mathbf{x}, \mathbf{y})\left(|\text{PWL}\_2(\mathbf{x}, \mathbf{y})| + |\text{PWL}\_1(\mathbf{x}, \mathbf{y})|\right)\right), \end{aligned} \tag{27}$$

where *g*(*x*, *y*) is the expression that models all obstacles in the proposed MTHPPM, represented by

$$\log(\mathbf{x}, \mathbf{y}) = \mathcal{W}\_{\mathbb{C}}(\mathbf{x}, \mathbf{y}) + \mathcal{W}\_{\mathbb{R}}(\mathbf{x}, \mathbf{y}), \tag{28}$$

*WC*(*x*, *y*) and *WR*(*x*, *y*) are the mathematical representation of circular and ellipsoidal obstacles, modeled by Equation (12) and Equation (9), respectively. Then, the homotopy system (6) is redefined by Equation (32), Equation (33), and Newton's homotopy formulation, obtaining

$$H = \begin{cases} H\_1(f\_{\text{PWM}\_1}(\mathbf{x}, y), \lambda) &= f\_{\text{PWM}\_1}(\mathbf{x}, y) \\ & -(1 - \lambda) f\_{\text{PWM}\_1}(\mathbf{x}\_0, y\_0) = 0, \\ H\_2(f\_{\text{PWM}\_2}(\mathbf{x}, y), \lambda) &= f\_{\text{PWM}\_2}(\mathbf{x}, y) \\ & -(1 - \lambda) f\_{\text{PWM}\_2}(\mathbf{x}\_0, y\_0) = 0, \end{cases} \tag{29}$$

where (*x*0, *y*0) is the initial point of the trajectory.

In addition to the PWL auxiliary lines, the second expression (27) is reformulated to keep the properties of obstacles for the multiple-target points approach. This reformulation guarantees that, when PWL1(*x*, *y*) = 0 and PWL2(*x*, *y*) = 0, the effect of the obstacles function *g*(*x*, *y*) over the homotopy curve is cancelled. In other words, the effect of the obstacles over the homotopy system progressively vanishes as the robot approaches to any target point, that is, this effect decreases in proportion to the difference between PWL1(*x*, *y*) and PWL2(*x*, *y*). The absolute values are employed in this formulation because any change in the sign of *g*(*x*, *y*) modifies the behavior of the homotopy curve, as it is presented in References [10,20]. It can be noticed that absolute value terms produce a singular Jacobian matrix during the iterative process. Then, in order to solve these issues and improve the performance of MGHPPM, a differentiable approximation of absolute value is used.

#### *5.1. Approximation of Absolute Value Function to Improve the Piecewise Linear Approach*

From definition of *f*(*x*) = |*x*|, *f*(*x*) is continuous for all values of *x*. Nevertheless, it is not differentiable when *x* = 0. Then, the absolute value function is not differentiable at the breakpoints, that is, for all points (*x*, *y*) such that *x* = *xB*,*<sup>i</sup>* (see expression (24)). In order to remove this mathematical issue in the procedure, an approximation of the absolute value function is employed. For this work, the approximation of the absolute value function presented in Reference [56] is used to warrant the continuity of the PWL derivative function. This formulation is defined by

$$|\mathbf{x}| \approx \xi(\mathbf{x}, \boldsymbol{\mathfrak{a}}) = \frac{1}{\mathfrak{a}} \left( \ln \left( 1 + \mathbf{e}^{-\mathbf{a} \cdot \mathbf{x}} \right) + \ln \left( 1 + \mathbf{e}^{\mathbf{a} \cdot \mathbf{x}} \right) \right), \tag{30}$$

where *α* is a parameter which reduces the error between |*x*| and *ξ*(*x*). Figure 7a shows the approximation of *ξ*(*x*, *α*) to |*x*| when *α* → ∞, this can also can be observed in Figure 7c. The first derivation graphs for *ξ*(*x*, *α*) and |*x*| are presented in Figure 7b. It shows that, like in the previous figure, the approximation *ξ*(*x*, *α*) is more similar to |*x*| when *α* → ∞.

**Figure 7.** Approximation of absolute value function. (**a**) Approximation of absolute value function. (**b**) First derivation of the approximation of the absolute value function. (**c**) Error (E) between the approximation of absolute value function and |*x*|.

By using this approximation of absolute value function, the canonical PWL(*x*, *y*) representation is modified to generate smooth and differentiable auxiliary PWL functions (PWL - (*x*, *y*)) for MGHPPM. The PWL - (*x*, *y*) representation is expressed as

$$\text{P}\overline{\text{W}}\overline{\text{L}}(\mathbf{x},\mathbf{y}) = \mathbf{y} - \left(\mathbf{a} + b\mathbf{x} + \sum\_{i=1}^{i=\sigma} \mathbf{c}\_{i} \left(\xi(\mathbf{x} - \mathbf{x}\_{B,i}, \mathbf{a})\right)\right), \tag{31}$$

Figure 8a presents a PWL - function generated with the canonical representation for *<sup>α</sup>* = 1, 5, 10, 50, 1 × 103, 1 × 109 in *<sup>ξ</sup>*(*x*, *<sup>α</sup>*). The first derivation of PWL - (*x*, *<sup>y</sup>*) is depicted in Figure 8b which shows the discontinuities present in PWL (*x*, *y*). It is important to note that for values of *α* ≤ 50, the behaviour of PWL - (*x*, *y*) is unsuitable for the MTHPPM formulation due to the inaccuracy between PWL(*x*, *y*) and PWL - (*x*, *y*) for these values of *α*, as it is depicted in Figure 8a,b. Nevertheless, a good approximation of PWL(*x*, *<sup>y</sup>*) is reached with *<sup>α</sup>* = <sup>1</sup> × 103 and *<sup>α</sup>* = <sup>1</sup> × 109 in PWL - (*x*, *<sup>y</sup>*) with the advantage that these have continuous differentiation, as it is shown in Figure 8b.

**Figure 8.** Smooth and differentiable piecewise linear function (PWL - ). (**a**) PWL - function generated with approximation of absolute value function. (**b**) First derivation of smooth PWL representation. -

For the MGHPPM, breakpoints are only used to generate intersections between auxiliary functions, in this sense, a good approximation that achieves this statement is when *<sup>α</sup>* = <sup>1</sup> × 103 (see Figure 8a). Then, this value is used for all case studies presented in the following section where, in order to simplify the notation, *<sup>ξ</sup>*(*x*, 1 × 103) is denoted by *<sup>ξ</sup>*(*x*) and expresions (26), (27) and (29) are redefined as

$$f\_{\bar{\text{PWL}}\_1}(\mathbf{x}, \mathbf{y}) = \mathbf{P}\bar{\text{W}}\mathbf{L}\_1(\mathbf{x}, \mathbf{y}), \tag{32}$$

$$\begin{aligned} f\_{\widehat{\rm PWL}\_2}(\mathbf{x}, y) &= \quad \widehat{\rm PWL}\_2(\mathbf{x}, y) \\ &- \quad \mathcal{g}(\mathbf{x}, y) \left( \xi \left( \widehat{\rm PWL}\_2(\mathbf{x}, y) \right) + \xi \left( \widehat{\rm PWL}\_1(\mathbf{x}, y) \right) \right), \end{aligned} \tag{33}$$

$$H = \begin{cases} H\_1(f\_{\overline{\rm PWL}\_1}(\mathbf{x}, y), \lambda) &= f\_{\overline{\rm PWL}\_1}(\mathbf{x}, y) \\ & -(1 - \lambda) f\_{\overline{\rm PWL}\_1}(\mathbf{x}\_0, y\_0) = 0, \\ H\_2(f\_{\overline{\rm PWL}\_2}(\mathbf{x}, y), \lambda) &= f\_{\overline{\rm PWL}\_2}(\mathbf{x}, y) \\ & -(1 - \lambda) f\_{\overline{\rm PWL}\_2}(\mathbf{x}\_0, y\_0) = 0, \end{cases} \tag{34}$$

#### *5.2. Breakpoints Selection for the Piecewise-Linear Functions*

A set of proper auxiliary PWL - functions generate the least number of unnecessary turns over the homotopy path. For this work, the way to obtain these are from two different initial points *y*(0) and the same breakpoints on the *x*-axis for both. This strategy guarantees that functions are different and intersect only at the desired points. Breakpoints are calculated from initial-points and target-points using the expressions

$$\mathbf{x}\_{B,i} = \mathbf{x}\_{T,i} + \frac{\mathbf{x}\_{T,i+1} - \mathbf{x}\_{T,i}}{\rho}, \ \rho = 2, 3, 4, \dots, <\infty \tag{35}$$

$$\mathcal{Y}\_{B,i} = \mathcal{J}\_i \ge\_{B,i} + \mathcal{Y}\_{B,i-1} - \mathcal{J}\_i (\mathbf{x}\_{B,i-1} - \mathbf{x}\_{B,i})\_\prime \tag{36}$$

$$\mathcal{J}\_i = \frac{y\_{T,i} - y\_{B,i-1}}{x\_{T,i} - x\_{B,i-1}},\tag{37}$$

where *<sup>ρ</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> is a proportionality parameter and determines the proximity between the *<sup>i</sup>*-th solution point and the *i*-th breakpoint (*i* = 1, 2, 3, ... , *σ*); for *σ* number of breakpoints there are *σ* + 1 target-points; (*xB*,*i*, *yB*,*i*) is the position of the *i*-th breakpoint; (*xT*,*i*, *yT*,*i*) is the position of the *i*-th target-point; J*<sup>i</sup>* is the slope of the *i*-th segment of the PWL - function, and (*xB*,0 = 0, *yB*,0 = *y*(0)) is the initial-point of the PWL - function. A simple environment map with two circular obstacles and two variants of PWL - <sup>1</sup> and PWL - <sup>2</sup> configurations is shown in Figure 9. For the PWL functions in Figure - 9a, the parameter is *ρ* = 2; for Figure 9b, the parameter *ρ* = 20 means that in this figure breakpoints are closer to the target-points.

**Figure 9.** Selection of breakpoints for the PWL - functions. (**a**) Breakpoints placed far from target-points. (**b**) Breakpoints placed close to target-points.

Figure 9b shows the similarity between PWL - <sup>1</sup> and PWL - 2, nevertheless, these only intersect at target-points. Figure 9a,b show the effect of PWL - <sup>1</sup> and PWL - <sup>1</sup> in the homotopy formulation through Equation (33), where, the similarity of these in the homotopy formulation produces path with a similar tendency than PWL - 's. Figure 9b shows two PWL - 's which are overlapped to obstacles, this configuration produces that homotopy path pass very close to obstacles. This characteristic could

be desirable in some cases, but when the map contains several grouped obstacles, the solution path could be inefficient in terms of length. In order to obtain paths very close to direct path between two points, for this work, the breakpoints are placed close to the target-points (similar form than Figure 9b).

#### *5.3. Technique for Successful Convergence and Avoid Reversal Effect*

Reversal effect is a phenomenon inherent to spherical path tracking and it is one of the most complex problems for the tracking techniques. Some works have proposed strategies to deal with this issue [10,47,48], nevertheless, these are inefficient when a PWL - approach is attempted. This section introduces a new strategy to avoid reversal effect during the tracking and improve the convergence of the corrector scheme. The proposed strategy is named double spherical tracking (DST); it is a spherical tracking embedded into another spherical tracking. DST is executed when a reversal effect or a non-convergence is detected. On the one hand, the reversal effect is recognized using the directional cosine strategy as explained in References [10,48]. On the other hand, a non-convergence is detected when the maximum number of iterations is reached without meeting the stop criterion. DST procedure is explained using the non-convergence case shown in Figure 10. Figure 10a shows a simple map with one obstacle and two target-points. Figure 10b shows the position of the predictor point (*xp*, *yp*, *λp*) which is far from the intersection between *H*1(*x*, *y*, *λ*), *H*2(*x*, *y*, *λ*), and *Si*(*x*, *y*, *λ*). DST formulation is based on the principle that <sup>∀</sup>(*x*, *<sup>y</sup>*, *<sup>λ</sup>*)*sol* <sup>∈</sup> *<sup>H</sup>*−<sup>1</sup> *<sup>S</sup>* (0) <sup>⇒</sup> (*x*, *<sup>y</sup>*, *<sup>λ</sup>*)*sol* <sup>∈</sup> *<sup>H</sup>*−<sup>1</sup> 1,*<sup>S</sup>* (0). Where, *<sup>H</sup>*−<sup>1</sup> *<sup>s</sup>* (0) is the solution set for the system of Equations (15) and *H*−<sup>1</sup> 1,*<sup>S</sup>* (0) is the solution set for the system of equations

$$H\_{1, \mathcal{S}} = \begin{cases} H\_1(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = \mathbf{0}, \\ S\_i(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = \mathbf{0}, \end{cases} \tag{38}$$

The procedure steps of DST are described as follows:

1. The system of equations *H*DST is established from *H*1(*x*, *y*, *λ*), *Si*(*x*, *y*, *λ*), and the sphere DST(*x*, *y*, *λ*). DST has been formulated to track the curve of intersection between *H*1(*x*, *y*, *λ*) and *Si*(*x*, *y*, *λ*) using the SA algorithm.

$$H\_{\rm DST} = \begin{cases} H\_1(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = \mathbf{0}, \\ S\_i(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = \mathbf{0}, \\ \rm{DST}\_k(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) = \mathbf{0}, \end{cases} \tag{39}$$

$$\begin{split} \text{DST}\_k(\mathbf{x}, \mathbf{y}, \boldsymbol{\lambda}) &= (\mathbf{x} - \mathbf{d}\mathbf{c}\_x)^2 + (\boldsymbol{y} - \mathbf{d}\mathbf{c}\_y)^2 \\ &+ (\boldsymbol{\lambda} - \mathbf{d}\mathbf{c}\_{\boldsymbol{\lambda}})^2 - r\_{\mathrm{dst}}^2 = \mathbf{0}, \end{split} \tag{40}$$

where *k* = 1, 2, 3, ... , *n*, *n* is the number of DST steps; *r*dst is the radius, and (dc*x*, dc*y*, dc*λ*) is the center of the sphere for every *k*-th step of the DST.


For his paper, the DST stop criterion is based on the distance between initial point of DST procedure (*xi*−1, *yi*−1, *λi*−1) and the current DST solution (dc*x*, dc*y*, dc*λ*)*<sup>k</sup>* as

$$\left\| \left( (\mathbf{dc}\_{x}, \mathbf{dc}\_{y}, \mathbf{dc}\_{\lambda})\_{k} - (x\_{i-1}, y\_{i-1}, \lambda\_{i-1}) \right) \right\| \prec^{r} r\_{\text{dst}} \tag{41}$$

here, the radius of DST sphere is proposed as *r*dst = (2*πrs*)/*n*dst, where *rs* is the radius of the sphere *Si* of the spherical tracking. *n*dst is the minimum number of steps of DST which is set to 24.

4. Finally, all points (dc*x*, dc*y*, dc*λ*)*k*, except (dc*x*, dc*y*, dc*λ*)1, are evaluated in Equation (15). Then, the point for which the evaluation is closer to zero, *HS*((dc*x*, dc*y*, dc*λ*)*k*) ≈ 0, is taken as a new predictor in the SA tracking of *HS*. Figure 10d shows the new predictor point (*xp*, *yp*, *λp*)new which is a point from the solution set *H*−<sup>1</sup> DST(0). Figure 10e provides a closer view where (*xp*, *yp*, *λp*)new is close to the intersection between *H*1(*x*, *y*, *λ*), *H*2(*x*, *y*, *λ*) and *Si*(*x*, *y*, *λ*). Then, this point is taken as the predictor point in the Broyden's corrector scheme. The position of the new point (*xi*+1, *yi*+1, *λi*+1) can be seen in Figure 10f.

**Figure 10.** Double spherical tracking strategy. (**a**) Environment map with two target-points. (**b**) Non-convergence issue. (**c**) Double spherical tracking *H*<sup>1</sup> and *Si*. (**d**) Double spherical tracking. (**e**) Double spherical tracking (zoom view). (**f**) New predictor-point.

It is important to remark that, DST is a backup technique for SA which is only applied when a non-convergence or reversal phenomenon is detected. Furthermore, the execution time of this technique is comparable with the execution time spent by one step of the SA. This because the *H*DST system of equations does not contain the obstacles representation in its formulation. For the map shown in Figure 10a, the DST strategy is employed only once for the non-convergence point found during the tracking. Figure 11 shows the 2-D and 3-D representation of the NAES and homotopy surfaces, the region where is located the non-convergence point is indicated.

**Figure 11.** 2-D and 3-D representation of the solution path. (**a**) Solution path 2-D representation. (**b**) Solution path 3-D representation (View 1). (**c**) Solution path 3-D representation (View 2).

Figure 11a shows the successful solution path obtained once that non-convergence is solved. Figure 11b,c show two different 3-D views of the surfaces which correspond at each equation of *HS*. In these figures, it can be noticed that the solution path represents the intersection between *H*1(*x*, *y*, *λ*) and *H*2(*x*, *y*, *λ*), crossing the two target-points. In this section, the DST technique for solving a non-convergence is explained and its ability to solve the reversion phenomenon is demonstrated.

#### *5.4. A Dummy Obstacle to Improve the Spherical Algorithm Performance*

Spherical tracking is an iterative process in which execution time and success depends on three main factors: predictor-corrector scheme, the radius of the sphere, and length of the path. First, the predictor-corrector scheme defines the total execution time, while the radius of the sphere and length of the path determines the number of predictor-corrector executions. For each SA step, the predictor is only employed one time, meanwhile, the number of iterations in Broyden's corrector depends on the predictor point. In other words, the success and fast convergence of the Broyden's corrector scheme depends on the quality of the predictor point. For quality, it means that proximity of the predictor point to the intersection between homotopy curve and the *i*-th sphere intersection (SA solution point). Then, the time for tracking an entire homotopy curve, from initial-point to target-point, is mostly determined by the total number of corrector scheme executions. In addition, the execution time for each corrector scheme iteration is the same for all SA steps; this implies a linear relationship between the total number of corrector iterations and the total computing time, as explained in Reference [10]. Second, if the radius of the sphere is reduced, the SA needs more steps to cover the same distance that could be traced by bigger spheres. Nonetheless, the sphere radius is restricted by the size of the obstacles, that is, if the sphere is bigger than the obstacles the SA could lose the path [10]. Third, the number of corrector executions increases as the homotopy curve length increases. The increase in length is due to the homotopy system formulation, as explained in Reference [10]. For homotopy with PWL - formulation, the length of the curve drastically increases caused by turnings when slopes changes. For the implementation in this work, the homotopy curves generated by PWL - formulation (34) are mapped in the space *x* ∈ (0, 1), *y* ∈ (0, 1) and *λ* ∈ (−∞, +∞). The plane *x* − *y* is delimited by the normalization of the workspace but *λ* is unlimited, nevertheless, for homotopy methods, the region of interest is *λ* ∈ (0, 1). Figure 11b,c show that displacement of the homotopy path ever the *λ*-axis is greater than over the *x* − *y* plane. Then, in order to minimize the *λ*-axis variations, a small perturbation in *H*1(*x*, *y*, *λ*) is proposed to generate a significant modification in the surface *H*1. This perturbation is performed by *Cd* which is a dummy obstacle placed very close to the initial point (*x*0, *y*0) with a radius smaller than the spheres employed in the SA. This obstacle generates a change of orientation of surface *H*1(*x*, *y*, *λ*), and their effect on the homotopy system is negligible for points far away from (*x*0, *y*0). By using the dummy obstacle, the surface *H*1(*x*, *y*, *λ*)

is forced to change its orientation. This change relates to the nature of Newton's homotopy which implies two properties in the new formulation: (a) its representation of the obstacles in the 3-D space are cylinders and rectangular prisms projections which are parallel to *λ*-axis. This can be concluded from the projections of Equations (7) and (10) expressions from 2-D to 3-D. Then, the new function that includes the dummy obstacle is represented by

$$\begin{split} f\_{\text{P}\overline{\text{W}}\mathbb{L}\_{1}}(\mathbf{x},y) &= \mathbf{P}\overline{\text{W}}\mathbb{L}\_{1}(\mathbf{x},y) \\ &\quad - \left(\frac{\mathbf{P}\_{\mathbb{C},d}}{\mathbb{C}\_{d}}\right) \left(\mathbb{J}^{\mathsf{x}}\left(\mathbf{P}\overline{\text{W}}\mathbb{L}\_{2}(\mathbf{x},y)\right) + \mathbb{J}\left(\mathbf{P}\overline{\text{W}}\mathbb{L}\_{1}(\mathbf{x},y)\right)\right), \end{split} \tag{42}$$

where, *pC*,*<sup>d</sup>* is the repulsion parameter of *Cd*. For this work, *pC*,*<sup>d</sup>* = 0.1, this value guarantees a minimal effect of the dummy obstacle in regions close to the starting point. b) The dummy obstacle integration in the homotopy system implies that PWL - <sup>2</sup>(*x*, *y*) must be modified to guarantee a continuous homotopy curve from initial-point to target-points. Figure 12 shows the effect of the *Cd* on the surface *H*1(*x*, *y*, *λ*). This surface is flattened and their variations over *λ*-axis is small for points away from the dummy obstacle. The use of a dummy obstacle allowed to reduce, substantially, the sweep on *λ*-axis for the homotopy system in the region away from this obstacle.

**Figure 12.** Map containing a dummy obstacle with isolated intersection zone. (**a**) Homotopy surfaces representation (View 1). (**b**)Homotopy surfaces representation (View 2).

The dummy obstacle represents a singularity in the flat surface, then for points closer to it, the variations over *λ*-axis are noticeable as shown in Figure 12a. It can observe that a non-intersection zone created an inclination on surface *H*2(*x*, *y*, *λ*), caused by its PWL - formulation. Figure 12b shows the isolated zone in which the spherical tracking should be confined. From the previous analysis, a solution for the isolated region is proposed and performed by the formulation of *Cd*

$$\mathbf{C}\_{d} = \left(\mathbf{x} - \mathbf{x}\_{\mathbb{C},d}\right)^{2} + \left(y - y\_{\mathbb{C},d}\right)^{2} - \left(r\_{\mathbb{C},d}\right)^{2},\tag{43}$$

$$\begin{aligned} \mathbf{x}\_{\mathbb{C},d} &= \mathbf{x}\_{0} - r\_{\mathbf{s},} \\ y\_{\mathbb{C},d} &= y\_{0} - \mathcal{J}\_{1}\left(\mathbf{x}\_{0} - \mathbf{x}\_{\mathbb{C},d}\right), \\ \mathcal{J}\_{1} &= \frac{y\_{\mathbb{C},1} - y\_{0}}{\mathbf{x}\_{\mathbb{C},1} - \mathbf{x}\_{0}}, \end{aligned}$$

here, the radius of dummy obstacle *rC*,*<sup>d</sup>* = *rs* <sup>100</sup> is small enough to be disregarded by the spherical tracking; (*x*0, *y*0) is the start point, and (*xG*,1, *yG*,1) is the first target-point. From the previous analysis, the initial point (*xB*,0, *yB*,0) of PWL - <sup>2</sup>(*x*, *y*) which guarantees a continuous solution path from start to every target-points is calculated using the next expression

$$(\mathbf{x}\_{B,0}, \mathbf{y}\_{B,0}) = \begin{pmatrix} 0, \mathbf{y}\_0 - \mathcal{J}\_1 \ \mathbf{x}\_0 \end{pmatrix},\tag{44}$$

Figure 13a depicts the dummy obstacle with center in (*xC*,*d*, *yC*,*d*) and using Equation (44) the point (*xB*,0, *yB*,0) for the PWL - <sup>2</sup>(*x*, *y*) is calculated. The successful solution path for multiple-target problem is shown in Figure 13b. Figure 13c,d show two views of the homotopy surfaces, *H*1(*x*, *y*, *λ*) and *H*2(*x*, *y*, *λ*). In these, it can notice the successful path begins at the start point and visits the two projected target-points. Figure 13d shows the effect of the dummy obstacle and the point (*xB*,0, *yB*,0) of PWL - <sup>2</sup>(*x*, *y*) to improve the homotopy path (in terms of length and number of SA stems). Furthermore, in this figure, the continuous curve denoted by the intersection between the homotopy surfaces *H*<sup>1</sup> and *H*<sup>2</sup> from the start point through all target-points is showed. The improvement of path (in terms of length) by using the dummy obstacle strategy can be validated from the visual comparison between the paths in the Figures 11 and 13.

**Figure 13.** Map with added dummy obstacle. (**a**) Environment map with dummy obstacle and two target-points. (**b**) 2-D representation solution path. (**c**) [3-D representation (View 1) solution path. (**d**) 3-D representation (View 2) solution path.

#### *5.5. Strategy to Simplified the Jacobian Matrix Based on Symbolic Manipulation*

To reduce the computational complexity and provide a simple implementation in *C*++ programming language, without the use of scientific specialized libraries, a symbolic manipulation of the Jacobian matrix (45) is developed. Jacobian matrix is the highest cost processes embedded in the SA; it is employed by predictor and corrector schemes for each step of the tracking. Furthermore, some terms also are employed on every step of DST.

$$J(\mathbf{x}, y, \lambda) = \begin{pmatrix} \frac{\partial H\_1(\mathbf{x}, y, \lambda)}{\partial \mathbf{x}} & \frac{\partial H\_1(\mathbf{x}, y, \lambda)}{\partial y} & \frac{\partial H\_1(\mathbf{x}, y, \lambda)}{\partial \lambda} \\\\ \frac{\partial H\_2(\mathbf{x}, y, \lambda)}{\partial \mathbf{x}} & \frac{\partial H\_2(\mathbf{x}, y, \lambda)}{\partial y} & \frac{\partial H\_2(\mathbf{x}, y, \lambda)}{\partial \lambda} \\\\ \frac{\partial S\_i(\mathbf{x}, y, \lambda)}{\partial \mathbf{x}} & \frac{\partial S\_i(\mathbf{x}, y, \lambda)}{\partial y} & \frac{\partial S\_i(\mathbf{x}, y, \lambda)}{\partial \lambda} \end{pmatrix} \tag{45}$$

The terms corresponding to the sphere *Si* and the partial derivatives to *λ* for *H*<sup>1</sup> and *H*<sup>2</sup> in the Jacobian matrix can be calculated by

1. From Equations (34) and (42)

$$\frac{\partial H\_1(\mathbf{x}, y, \lambda)}{\partial \lambda} = \overline{\text{PWL}}\_1(\mathbf{x}0, y0) - \left(\frac{\mathcal{P}\mathbb{C}\_d}{\mathbb{C}\_d}\right) \left(\overline{\xi} \left(\overline{\text{PWL}}\_2(\mathbf{x}0, y0)\right) + \overline{\xi} \left(\overline{\text{PWL}}\_1(\mathbf{x}0, y0)\right)\right),\tag{46}$$

2. From Equations (34) and (33)

$$\frac{\partial H\_{2}(\mathbf{x},\boldsymbol{y},\boldsymbol{\lambda})}{\partial \boldsymbol{\lambda}} = \widehat{\mathrm{PW}}\_{2}^{-}(\mathbf{x}\_{0},\boldsymbol{y}\_{0}) - \boldsymbol{g}(\mathbf{x}\_{0},\boldsymbol{y}\_{0}) \left( \boldsymbol{\xi}^{\mathrm{r}} \Big( \widehat{\mathrm{PW}}\_{2}^{-}(\mathbf{x}\_{0},\boldsymbol{y}\_{0}) \Big) + \boldsymbol{\xi}^{\mathrm{r}} \Big( \widehat{\mathrm{PW}}\_{1}^{-}(\mathbf{x}\_{0},\boldsymbol{y}\_{0}) \Big) \right), \tag{47}$$

3. From Equation (16)

$$\frac{\partial \mathcal{S}\_i(\mathbf{x}, y, \lambda)}{\partial \mathbf{x}} = 2(\mathbf{x} - \mathbf{c}\_x),\tag{48}$$

$$\frac{\partial S\_i(x, y, \lambda)}{\partial y} = 2(y - c\_y),\tag{49}$$

$$\frac{\partial S\_i(x, y, \lambda)}{\partial \lambda} = 2(\lambda - c\_{\lambda})\_{\prime} \tag{50}$$

Using symbolic manipulation and the chain rule to obtain the derivative of composite functions, and knowing that evaluation of *f* PWL - <sup>1</sup> (*x*0, *y*0) and *f* PWL - <sup>2</sup> (*x*0, *y*0) are constants, the rest of terms can be reduced as

1. From Equations (34) and (42)

$$
\begin{split}
\frac{\partial \left( H\_{1}(\mathbf{x}, \mathbf{y}, \lambda) \right)}{\partial \mathbf{x}} &=& \frac{\partial \left( \overline{\mathbf{P}} \overline{\mathbf{W}} \mathbf{L}\_{1}(\mathbf{x}, \mathbf{y}) \right)}{\partial \mathbf{x}} - \frac{\partial \left( \frac{p\_{\mathbb{C},d}}{\mathbf{C}\_{d}(\mathbf{x}, \mathbf{y})} \right)}{\partial \mathbf{x}} \left( \xi \left( \overline{\mathbf{P}} \overline{\mathbf{W}} \mathbf{L}\_{2}(\mathbf{x}, \mathbf{y}) \right) + \xi \left( \overline{\mathbf{P}} \overline{\mathbf{W}} \mathbf{L}\_{1}(\mathbf{x}, \mathbf{y}) \right) \right) \\ &- \frac{\partial \left( \frac{\eta}{\xi} \left( \overline{\mathbf{P}} \overline{\mathbf{W}} \mathbf{L}\_{2}(\mathbf{x}, \mathbf{y}) \right) + \frac{\eta}{\xi} \left( \overline{\mathbf{P}} \overline{\mathbf{W}} \overline{\mathbf{L}}\_{1}(\mathbf{x}, \mathbf{y}) \right) \right)}{\partial \mathbf{x}} \left( \frac{p\_{\mathbb{C},d}}{\mathbf{C}\_{d}(\mathbf{x}, \mathbf{y})} \right),
\end{split}
\tag{51}$$

$$\begin{split} \frac{\partial \left( H\_{1}(\mathbf{x}, y, \boldsymbol{\lambda}) \right)}{\partial y} &= \quad \frac{\partial \left( \overline{\mathrm{P}} \overline{\mathrm{W}} \mathrm{L}\_{1}(\mathbf{x}, y) \right)}{\partial y} \\ &- \frac{\partial \left( \frac{\mathrm{P}\_{\complement\mathcal{I}}}{\mathrm{C}\_{4}(\mathbf{x}, y)} \right)}{\partial y} \left( \overline{\mathrm{f}} \left( \overline{\mathrm{P}} \overline{\mathrm{W}} \mathrm{L}\_{2}(\mathbf{x}, y) \right) + \overline{\mathrm{f}} \left( \overline{\mathrm{P}} \overline{\mathrm{W}} \mathrm{L}\_{1}(\mathbf{x}, y) \right) \right) \\ &- \frac{\partial \left( \overline{\mathrm{f}} \left( \overline{\mathrm{P}} \overline{\mathrm{W}} \mathrm{L}\_{2}(\mathbf{x}, y) \right) + \overline{\mathrm{f}} \left( \overline{\mathrm{P}} \overline{\mathrm{W}} \mathrm{L}\_{1}(\mathbf{x}, y) \right) \right)}{\partial y} \left( \frac{p\_{\complement\mathcal{I}}}{\mathrm{C}\_{4}(\mathbf{x}, y)} \right), \end{split} \tag{52}$$

2. From Equations (34) and (33)

$$\begin{array}{rcl}\frac{\partial\left(H\_{2}(\mathbf{x},\boldsymbol{y},\boldsymbol{\lambda})\right)}{\partial\mathbf{x}} &=& \frac{\partial\left(\overline{\operatorname{PWD}}\_{2}(\mathbf{x},\boldsymbol{y})\right)}{\partial\mathbf{x}} - \frac{\partial\left(\operatorname{g}(\mathbf{x},\boldsymbol{y})\right)}{\partial\mathbf{x}} \left(\overline{\operatorname{g}}\left(\overline{\operatorname{PWD}}\_{2}(\mathbf{x},\boldsymbol{y})\right) + \overline{\operatorname{g}}\left(\overline{\operatorname{PWD}}\_{1}(\mathbf{x},\boldsymbol{y})\right)\right) \\ & \frac{\partial\left(\overline{\operatorname{g}}\left(\overline{\operatorname{PWD}}\_{2}(\mathbf{x},\boldsymbol{y})\right) + \overline{\operatorname{g}}\left(\overline{\operatorname{PWD}}\_{1}(\mathbf{x},\boldsymbol{y})\right)\right)}{\partial\mathbf{x}} \left(\operatorname{g}(\mathbf{x},\boldsymbol{y})\right), \end{array} \tag{53}$$

$$\begin{split} \frac{\partial \left( H\_{2}(\mathbf{x}, y, \lambda) \right)}{\partial y} &= \quad \frac{\partial \left( \mathbb{P} \overline{\text{W}} \mathbb{L}\_{2}(\mathbf{x}, y) \right)}{\partial y} - \frac{\partial \left( g(\mathbf{x}, y) \right)}{\partial y} \left( \mathbb{S} \left( \mathbb{P} \overline{\text{W}} \mathbb{L}\_{2}(\mathbf{x}, y) \right) + \mathbb{S} \left( \mathbb{P} \overline{\text{W}} \mathbb{L}\_{1}(\mathbf{x}, y) \right) \right)} \\ &- \frac{\partial \left( \mathbb{S} \left( \mathbb{P} \overline{\text{W}} \mathbb{L}\_{2}(\mathbf{x}, y) \right) + \mathbb{S} \left( \mathbb{P} \overline{\text{W}} \mathbb{L}\_{1}(\mathbf{x}, y) \right) \right)}{\partial y} \left( g(\mathbf{x}, y) \right), \end{split} \tag{54}$$

where *<sup>∂</sup>*(*g*(*x*,*y*)) *<sup>∂</sup><sup>x</sup>* and *<sup>∂</sup>*(*g*(*x*,*y*)) *<sup>∂</sup><sup>y</sup>* are

$$\frac{\partial \left( g(\mathbf{x}, \mathbf{y}) \right)}{\partial \mathbf{x}} = \sum\_{i=1}^{i=k} \frac{-p\_{O,i} \left( \frac{\partial (O\_i(\mathbf{x}, \mathbf{y}))}{\partial \mathbf{x}} \right)}{\left( |O\_i(\mathbf{x}, \mathbf{y})| + O\_i(\mathbf{x}, \mathbf{y}) \right)^2} \left( 1 + \frac{|O\_i(\mathbf{x}, \mathbf{y})|}{O\_i(\mathbf{x}, \mathbf{y})} \right), \tag{55}$$

$$\frac{\partial \left( g(\mathbf{x}, y) \right)}{\partial y} = \sum\_{i=1}^{i=k} \frac{-p\_{O,i} \left( \frac{\partial (O\_i(\mathbf{x}, y))}{\partial y} \right)}{\left( |O\_i(\mathbf{x}, y)| + O\_i(\mathbf{x}, y) \right)^2} \left( 1 + \frac{|O\_i(\mathbf{x}, y)|}{O\_i(\mathbf{x}, y)} \right), \tag{56}$$

where *Oi*(*x*, *y*) is the expression that describes the shape of the *i*-th obstacle; *Ci* for a circular obstacle, and *Ri* for an ellipsoidal obstacle. The parameter *pO*,*<sup>i</sup>* is the repulsion parameter of each obstacle; *pC*,*<sup>i</sup>* for a circular obstacle, and *pR*,*<sup>i</sup>* for an ellipsoidal obstacle.

$$\begin{split} \frac{\partial \left( \xi \left( \overrightarrow{\boldsymbol{m}\_{\mathrm{2}}} (\mathbf{x}, \boldsymbol{y}) \right) + \xi \left( \overrightarrow{\boldsymbol{m}\_{\mathrm{1}}} (\mathbf{x}, \boldsymbol{y}) \right) \right)}{\partial \mathbf{x}} &= \frac{1}{a} \frac{\partial \left( \overrightarrow{\boldsymbol{m}\_{\mathrm{2}}} (\mathbf{x}, \boldsymbol{y}) \right)}{\partial \mathbf{x}} \left( -\frac{a e^{-a \overrightarrow{\boldsymbol{P} \, \mathrm{W}} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}{1 + \,\_{\ast} \, -a \overrightarrow{\boldsymbol{P} \, \mathrm{W} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}}{1 + \,\_{\ast} \, a \overrightarrow{\boldsymbol{P} \, \mathrm{W} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}} \right) \\ &+ \frac{1}{a} \frac{\partial \left( \overrightarrow{\boldsymbol{m}\_{\mathrm{1}}} (\mathbf{x}, \boldsymbol{y}) \right)}{\partial \mathbf{x}} \left( -\frac{a e^{-a \overrightarrow{\boldsymbol{P} \, \mathrm{W}} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}{1 + \,\_{\ast} a \overrightarrow{\boldsymbol{P} \, \mathrm{W} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}} + \frac{a e^{a \overrightarrow{\boldsymbol{P} \, \mathrm{W} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}}}{1 + \,\_{\ast} a \overrightarrow{\boldsymbol{P} \, \mathrm{W} \mathrm{L}\_{1} (\mathbf{x}, \boldsymbol{y})}} \right), \end{split} \tag{57}$$

$$\begin{split} \frac{\partial \left( \xi \left( \widetilde{\boldsymbol{m}} \boldsymbol{\Sigma} \boldsymbol{\varepsilon} \mathbf{x} , \boldsymbol{y} \right) \right) + \xi \left( \widetilde{\boldsymbol{m}} \boldsymbol{\varepsilon} \mathbf{x} , \boldsymbol{y} \right) \rangle}{\partial \boldsymbol{y}} &= \frac{1}{a} \frac{\partial \left( \widetilde{\boldsymbol{m}} \boldsymbol{\varepsilon} \mathbf{x} , \boldsymbol{y} \right) \rangle}{\partial \boldsymbol{y}} - \left( \frac{a \boldsymbol{\varepsilon} - a \widetilde{\boldsymbol{\mathcal{W}} \mathbf{\tilde{W}} } \boldsymbol{\varepsilon} \left( \boldsymbol{x} , \boldsymbol{y} \right)}{1 + \boldsymbol{\varepsilon} \cdot \boldsymbol{a}^{\boldsymbol{\varepsilon} \mathbf{\tilde{W}} \mathbf{\tilde{L}} \boldsymbol{\varepsilon} \left( \mathbf{x} , \boldsymbol{y} \right)} \right) \\ &+ \frac{1}{a} \frac{\partial \left( \widetilde{\boldsymbol{m}} \boldsymbol{\varepsilon} \mathbf{\tilde{I}} \left( \mathbf{x} , \boldsymbol{y} \right) \right)}{\partial \boldsymbol{y}} \left( - \frac{a \boldsymbol{\varepsilon} - a \widetilde{\mathbf{P} \mathbf{\tilde{W}} } \boldsymbol{\varepsilon} \left( \mathbf{x} , \boldsymbol{y} \right)}{1 + \boldsymbol{\varepsilon} \cdot \boldsymbol{a}^{\boldsymbol{\varepsilon} \mathbf{\tilde{W}} \mathbf{\tilde{L}} \left( \mathbf{x} , \boldsymbol{y} \right)} + \frac{a \boldsymbol{\varepsilon}^{a} \widetilde{\mathbf{P} \mathbf{\tilde{W}} \mathbf{\tilde{L}}\_{1} \left( \mathbf{x} , \boldsymbol{y} \right)}{1 + \boldsymbol{\varepsilon} \cdot \boldsymbol{a}^{\boldsymbol{\varepsilon} \mathbf{\tilde{W}} \mathbf{\tilde{L}} \left( \mathbf{x} , \boldsymbol{y} \right)} \right)} . \end{split} \tag{58}$$

where, *<sup>α</sup>* = <sup>1</sup> × 103. Finally,

$$\frac{\partial \left( \text{P\"W}\_1(\mathbf{x}, y) \right)}{\partial \mathbf{x}} = - \left( b\_1 \mathbf{x} + \sum\_{i\_1=1}^{i\_1 = \sigma} c\_{i\_1} \left( \frac{\boldsymbol{\xi} \left( \mathbf{x} - \mathbf{x}\_{B, i\_1} \right)}{\mathbf{x} - \mathbf{x}\_{B, i\_1}} \right) \right), \tag{59}$$

$$\frac{\partial \left( \text{P} \overline{\text{W}} \text{L}\_2(\mathbf{x}, y) \right)}{\partial \mathbf{x}} = - \left( b\_2 \mathbf{x} + \sum\_{i\_2=1}^{i\_2=\sigma} c\_{i\_2} \left( \frac{\xi' \left( \mathbf{x} - \mathbf{x}\_{B, i\_2} \right)}{\mathbf{x} - \mathbf{x}\_{B, i\_2}} \right) \right), \tag{60}$$

$$\frac{\partial \left( \overline{\text{P}} \overline{\text{W}} \overline{\text{L}}\_1(\mathbf{x}, y) \right)}{\partial y} = 1,\tag{61}$$

$$\frac{\partial \left( \mathbf{P} \overline{\mathbf{W}} \overline{\mathbf{L}}\_2(\mathbf{x}, y) \right)}{\partial y} = 1,\tag{62}$$

where (*xB*,*ij* , *yB*,*ij* ) is the *i*-th breakpoint of the *j*-th PWL - functions PWL - <sup>1</sup> and PWL - 2, respectively. Every expression presented in this section has been implemented in our *C*++ program. The symbolic manipulation reduces the computing time because it only requires substitutions and evaluations.

#### *5.6. A Systematic Criterion to Select the Repulsion Parameter*

The selection of repulsion parameters for obstacles has been an open problem since the HPPM was proposed in Reference [19]. This is a complex task because the effect of each parameter has a great impact on the homotopy surfaces and, hence, in the solution curve. Advances in the characterization and creation of a criterion have been reported in previous works, where: A study of the effect produced by the sign and magnitude of each parameter over the homotopy curve was introduced in Reference [10], nevertheless, the criterion employed to set them was only based on empirical knowledge. Then, Reference [20] presents a criterion based on random selection and segmentation of the map to enhance the performance of the planner and find the shortest path. Also, a grouping obstacle strategy to reduce the number of repulsion parameters in the formulation was presented in Reference [20]. This strategy is applied to reduce the complexity of parameter selection tasks for maps with a large number of obstacles. Next, a first approach of systematic criterion based on the distance between obstacles and the diagonal path to properly select the repulsion parameter was presented in Reference [37]. Showed the effectiveness of a systematic criterion using examples, nonetheless, its performance takes a higher portion of the total execution time. This section is devoted to explaining a new systematic criterion based on results reported in Reference [10,20,37] and the properties of auxiliary PWL - functions proposed in this work. The main property of the formulation that has an effect over the repulsion parameter is denoted in the auxiliary function (33) which contains all mathematical representations of obstacles in the map. This function produces a scaling effect in the value of every repulsion parameter due to its absolute value terms. These reduce the effect of repulsion parameters in the vicinity of target-points and proportionally changing its value with the proximity between obstacles and any PWL function. From the manipulation of Equation ( - 33)

$$\begin{split} f\_{\text{P\overline{W}L}\_{2}}(\mathbf{x},\mathbf{y}) &= \mathbf{P}\bar{\text{W}L}\_{2}(\mathbf{x},\mathbf{y}) \\ &- \left( \sum\_{i=1}^{i=n} \frac{p\_{O\_{i}}\left(\xi\left(\mathbf{P}\overline{\text{W}L}\_{2}(\mathbf{x},\mathbf{y})\right) + \xi\left(\mathbf{P}\overline{\text{W}L}\_{4}(\mathbf{x},\mathbf{y})\right)\right)}{\mathcal{O}\_{i}(\mathbf{x},\mathbf{y}) + |\mathcal{O}\_{i}(\mathbf{x},\mathbf{y})|} \right), \end{split} \tag{63}$$

here, *i* = 1, 2, 3, ... *n*, *n* is the number of obstacles in the map; *Oi*(*x*, *y*) is the mathematical representation of any *i*-th circular obstacle, and *pO*,*<sup>i</sup>* represents its respective repulsion parameter. Then, the effective value of the parameter is denoted as

$$p\_{\rm eff,i} = p\_{O,i} \left( \zeta \left( \overline{\mathcal{P}} \overline{\text{W}} \overline{\text{L}}\_2(\mathbf{x}, y) \right) + \zeta \left( \overline{\mathcal{P}} \overline{\text{W}} \overline{\text{L}}\_1(\mathbf{x}, y) \right) \right), \tag{64}$$

where *p*eff,*<sup>i</sup>* is the effective value of repulsion parameter in Equation (33) for each *i*-th obstacle. It can be noticed that the value of *p*eff,*<sup>i</sup>* dynamically varies according to SA tracking the solution path, then, only one value of *p*base is selected for all circular obstacles. From the characterization presented in Reference [10], the next rank of magnitude for this parameter is employed in the case studies; 1 ≥ |*p*base| > 0 for circular obstacles and 100 ≥ |*p*base| > 1 for ellipsoidal obstacles. The expression (64) denotes that parameter magnitude of any obstacle decreases according to its proximity to any PWL function.

In this subsection, an environment map with three circular obstacles is used to explain the proposed systematic criterion for selecting the sign of repulsion parameters. As already explained above, solution path tends to follow the function PWL - <sup>2</sup>(*x*, *y*) as a consequence of using a dummy obstacle. In this sense, the criterion to set the sign of the repulsion parameter for any circular obstacle can be established from its center position with respect to PWL - <sup>2</sup>(*x*, *y*) (see Figure 14b). This criterion is based on that, for all points (*a*, *b*) located above PWL - <sup>2</sup>(*x*, *y*), the value of PWL - <sup>2</sup>(*a*, *b*) is greater than zero ; for points located below PWL - <sup>2</sup>(*x*, *y*), the value is less than zero (see Figure 14c). By using the dynamic value of repulsion parameter property provided by the formulation and relative position of the center ((*xC*,*i*, *yC*,*i*)) for every *i*-th circular obstacle with respect to PWL - <sup>2</sup>(*x*, *y*), a new definition of the repulsion parameter is described by

$$p\_{O,i} = \begin{cases} -p\_{\text{base}} \text{sgn} \left( \overline{\text{P}} \overline{\text{W}}\_2 (\mathbf{x}\_{\mathbb{C},i}, y\_{\mathbb{C},i}) \right) & \text{if} \quad \overline{\text{P}} \overline{\text{W}}\_2 (\mathbf{x}\_{\mathbb{C},i}, y\_{\mathbb{C},i}) \neq 0, \\ -p\_{\text{base}} & \text{if} \quad \overline{\text{P}} \overline{\text{W}}\_2 (\mathbf{x}\_{\mathbb{C},i}, y\_{\mathbb{C},i}) = 0, \end{cases} \tag{65}$$

where *pO*,*<sup>i</sup>* is the base parameter value with a sign that determines a lower or upper side path according to the *i*-th circular obstacle. Figure 14d shows the resulting solution path once the systematic criterion to set the repulsion parameter is executed. In this case study, the strategy is capable to properly set the signs of repulsion parameters which generate an efficient path with minimum length. This example validates the criterion effectiveness, which is addressed in the following sections.

The criterion introduced in this subsection has been formulated and characterized only for circular obstacles because of its symmetry. Nonetheless, it is assumed that ellipsoidal obstacles are only employed to represent walls or any other obstacle of big dimensions on the map. Then, the criterion to automatically assign the sign of repulsion parameter needs an additional step which will be explained in the next section.

**Figure 14.** Criterion to select the repulsion parameter. (**a**) Environment map. (**b**) Position of obstacles with respect to PWL - <sup>2</sup>(*x*, *y*, *λ*). (**c**) Positive and negative regions with respect to PWL - <sup>2</sup>(*x*, *y*, *λ*). (**d**) Solution path using systematic criterion to select repulsion parameter.

#### **6. Multiple-Target Homotopic Path Planning Method with Visibility Graph Approach**

Visibility graph (VG) is one of the most widely used roadmap methods to find the shortest path. This method employs a geometrical map representation to generate a roadmap that contains all links between vertices of the polygonal obstacles, start-point, and target-point [12–14]. Although this algorithm can find the shortest path for maps of structured environments, it is only employed to solve configurations with few polygonal obstacles, because its execution time and memory consumption depends on the number of obstacles [12]. In this work, the visibility graph algorithm is applied to pre-process the map and obtain a first approximation of the solution path for closed and office-like environments. In this way, only big polygonal objects such as walls and office furniture (static objects) are considered. Then, MTHPPM uses the path provided by VG (commonly the shortest path) to compute a collision-free path for the full environment containing non-polygonal obstacles. Figure 15 shows the operation process of MTHPPM using the visibility approach applying (MTHPPM\_VG) for an office-like environment with five hundred obstacles.

**Figure 15.** Visibility approach applied to Homotopic Path Planning Method. (**a**) Environment map with walls. (**b**) Visibility points and pre-processed path. (**c**) Visibility graph path without circular obstacles. (**d**) Solution path obtained by MTHPPM using the visibility approach.

For the example in Figure 15a, the main task of the mobile robot is to reach the target-point located on the opposite corner of the room, concerning to the initial position. Solving this example is done through the following steps; first, an approximation of the solution path is obtained by the visibility graph algorithm using a map, where only walls are considered. The configuration map without obstacles is reduced to a simple problem (see Figure 15b). The graph complexity is further reduced using the premise that PWL - formulation is a single-valued function. In this regard, every path with links that implies returns in *x*-axis is suppressed. Second, the first approximation of the path will be integrated by only links (corners of static obstacles and walls) in *x*-axis forward direction, as it is shown in Figure 15c. It can be noticed that, for environment maps with bug traps and maze layout, a path with visible vertices in forward direction might not exist. Two case studies with this problem are treated and solved in the next section. Third, the solution path provided by VG is employed to generate PWL - formulation which contains the information about the position of the ellipsoidal obstacles. Then, automatic sign assignation for this type of obstacle is executed, operating in the same way as circular obstacles. The function to obtain the sign of an ellipsoidal obstacle is

$$p\_{O,i} = \begin{cases} -p\_{\text{base}} \text{sgn} \left( \overline{\text{P}} \overline{\text{W}} \mathbb{L}\_2(\mathbf{x}\_{R,i}, y\_{R,i}) \right) & \text{if} \quad \overline{\text{P}} \overline{\text{W}} \mathbb{L}\_2(\mathbf{x}\_{R,i}, y\_{R,i}) \neq 0, \\ -p\_{\text{base}} & \text{if} \quad \overline{\text{P}} \overline{\text{W}} \mathbb{L}\_2(\mathbf{x}\_{R,i}, y\_{R,i}) = 0, \end{cases} \tag{66}$$

where *pO*,*<sup>i</sup>* is the base parameter value with a sign that determine if the solution path pass below or above the *i*-th ellipsoidal obstacle with center at point (*xR*,*i*, *yR*,*i*) and *p*base is the base value selected, arbitrarily, for all ellipsoidal obstacles; from empirical data presented in Reference [10], the value is in range 100 ≥ |*p*base| > 1. Finally, the MTHPPM is executed for the configuration map containing obstacles, walls, and vertices of the visibility path as target-points; solution is the shortest collision-free

path (see Figure 15d). The approach explained in this section can be used for problems with only one target-point or multiple-targets. For multiple-targets cases, the visibility points should be integrated in the set of targets as long as these do not represent a backward advance in the PWL - formulation, in which case these points must be removed.

The procedure of the proposed multiple-target planner can be summarised using the flow chart presented in Figure 16. If the planning problem has only a single target the method is named HPPM and for multiple targets MGHPPM. Furthermore, when the planning problem is performed on a structured workspace with a single target point, the procedure is the same and the visibility graph approach is used to calculate the first approximation of the solution, and the homotopy based method is named MGHPPM\_VG.

**Figure 16.** Multiple-target homotopic path planning method flow chart.

The flow chart (see Figure 16) contains all strategies explained previously and the addition of points generated by the visibility graph approach into the set of target-points. The proposed methodology considers that the environment map is known in advance or it is provided by sensors onboard (for real-time implementations) in the load environment map stage. In this sense, inputs in the proposed planner are a semi-algebraic representation of the environment map (*C*space), the initial position of the robot, and one or more target-points designed by the user or generated by a visibility-graph algorithm (because visibility graph planners have been treated strongly in others works such as References [12–14], in this work it is not integrated in the homotopic planner process, is only used as sub-process). Here, it is considered that the robot has a navigation module and control rules to guarantee the correct path following. Then, the iterative procedure of MTHPPM (for multiple-target points) or HPPM (for single target-point) starts. First, the mathematical model of the environment and data provided by the user is generated in the PWL - formulation stage. Second, the techniques to enhance convergence, length of path, and automatic assignation of repulsion parameters are performed (grouping obstacles, systematic criterion to select repulsion parameter, and dummy obstacle addition). Third, homotopy system of non-linear equations to represent the configuration space is formulated in NAES formulation stage. Then, the first predictor point for the HPPM (or MGHPPM) is calculated from the initial position of the robot (first point of the path). Fourth, the iterative procedure of spherical algorithm is executed; within it, the predictor-corrector scheme and reversion phenomenon (or non-convergence) check are performed for each step. This procedure is executed as explained in Section 3, where the Broyden's method is employed to calculate a new point of path which represents a solution of the system of homotopic equations and Euler's scheme is executed to obtain a new predictor point. During the Spherical algorithm procedure, if a reversion or non-convergence is detected, the DST technique is enabled and employed to solve this issue as already explained in the Section 5.3. Finally, the procedure finishes when the SA reaches a target-point or the maximum number of steps is achieved; then, returns a set of points that represent the solution path. In the next section, some case studies are presented to validate the utility of this methodology.

#### **7. Case Studies**

This section provides five case studies to provide certainty about the usefulness of the proposed methodology and how it solves path planning problems like bug traps and narrow corridors. The first three simulations are focused on how to solve planning problems in sceneries with narrow corridors. Simulation four demonstrates its effectiveness to solve maps with bug traps. On the one hand, for these four case studies, a comparison between the proposed methodology in its two variants (HPPM formulation and MTHPPM with visibility graph approach) and six sampling-based planning algorithms is provided. This comparison is based on important aspects like: memory consumption, execution time, percentage of feasible paths found, and path length. On the other hand, an example of multiple-target path planning for a closed environment with narrow corridors is shown in simulation five. All simulations were performed using the following set-up: all planners have the same step size (spherical tracking radius for MTHPPM and HPPM; collision checking resolution parameter for BSP algorithms), MTHPPM code in C++, and executed on PC (Intel i7 2.6 GHz processor, RAM 16 GB, and 64-bit Ubuntu 16.04 operating system). SBP algorithms ran on the same PC using the well known Open Motion Planning Library v.1.3.1 [57] and OMPL Planer Arena [58] to obtain performance data. To highlight advantages and weaknesses of the proposed homotopy based planner, it is compared to eight SBP algorithms: Expansive Space Trees (EST), Kinematic Planning by Interior-Exterior Cell Exploration (KPIECE1), Probabilistic Roadmap Method (PRM), Rapidly-exploring Random Trees (RRT), Bidirectional Rapidly-exploring Random Trees (Bi-RRT, also named RRT-Connect), Rapidly-exploring Random Trees with A\* approach (RRT\*), RRT and PRM with Lazy collision method (LazyRRT), and (Lazy PRM). To obtain a significant performance data, each SBP algorithm was run one hundred times.

#### *7.1. Case 1*

This case study is devoted to show the ability of MTHPPM to deal with the narrow corridor problem. Figure 17a depicts the visual comparison between HPPM and MTHPPM\_VG (MTHPPM with visibility graph approach) paths on a normalized office-like map in 2-D with five hundred circular obstacles. Here, the initial position of the robot is at point (0, 0) and target-point at (1, 1). The visual comparison between paths obtained by SBP algorithms is presented in Figure 17b.

**Figure 17.** Visual comparative between MTHPPM and HPPM against best run of each SBP algorithm for case study 1. (**a**) Paths obtained with HPPM and MTHPPM\_VG. (**b**) Paths obtained with SBP algorithms.

Figure 17a,b, the calculated paths by HPPM against MTHPPM\_VG and SBP algorithms (shortest path after 100 runs) are contrasted. It is important to note that the path obtained using HPPM is one of the longest because it tends is to round all obstacles in direction to the diagonal line between initial and target-point. On the other hand, the path obtained by MTHPPM with visibility approach is one of the shortest. Figure 17a shows the effectiveness of HPPM and MTHPPM (MTHPPM\_VG) to solve maps with narrow corridors and several obstacles. Furthermore, the ability of MTHPPM to calculate the shortest collision-free path using the visibility graph approach explained in Section 5 is validated in this case study. The box plots (box-and-whisker diagrams) in Figure 18 depicts a quantitative comparison between HPPM, MTHPPM\_VG, and SBP algorithms for percentage of feasible paths found, execution time, path length, and memory consumption.

Figure 18 shows the summarized results which denote the following characteristics: First, Figure 18a shows that percentage of fails for PRM is higher than twenty percent, while for LazyRRT is closer to that percentage. For HPPM and MTHPPM\_VG, the success rate is one hundred percent due to its deterministic formulation. Second, the execution time spent to solve the configuration map is presented in Figure 18c. In this box-and-whisker diagram, it can notice, that SBP algorithms are in the order of seconds, while MTHPPM\_VG and HPPM their time is in the order of milliseconds. Third, the results in Figure 18b show that the path obtained by MTHPPM\_VG is one of the shortest, RRT\* provided the best. Finally, memory consumption comparison exhibits a big gap between HPPM and MTHPPM against SBP algorithms. The difference is about three orders of magnitude, from KB to MB. This is because SBP algorithms store a roadmap of collision-free points, while homotopy based methods only store the path and obstacle positions.

**Figure 18.** Comparative results for case study 1. (**a**) Successful paths. (**b**) Path length. (**c**) Execution time. (**d**) Memory.

#### *7.2. Case 2*

This case study presents a path planning problem of an office-like environment with narrow corridors and five hundred obstacles. This example shows one of the future applications of HPPM and MTHPPM, to integrate it in the navigation system of a parcel service robot. The main task of the robot, in this example, is to collect and deliver packages or documents from one cubicle to another. Figure 19a provides the floor plan of an office represented on a normalized 2-D space, where, the initial and target points are located in opposites corners.

Figure 19a,b show a visual comparison of paths generated by HPPM against MTHPPM\_VG and between SBP algorithms, respectively. In these patterns, the difference between path lengths calculated by SBP algorithms, HPPM, and MTHPPM\_VG are depicted. It is important to note that, for SBP algorithms, the most optimistic simulation is taken, that is, the shortest path after one hundred runs. Figure 20 shows the box-and-whisker diagrams of the performance results for SBP algorithms and homotopy based planners (HPPM and MTHPPM).

**Figure 19.** Visual comparative for case study 2 with five hundred obstacles; best run of each SBP algorithm. (**a**) Paths obtained with HPPM and MTHPPM\_VG. (**b**) Paths obtained using SBP algorithms.

**Figure 20.** Comparative results for case study 2. (**a**) Successful paths. (**b**) Path length. (**c**) Execution time. (**d**) Memory.

From simulations, it can be concluded the following: First, RRT, PRM, RRT\*, and LazyRRT have a poor (less than eighty percent) performance of success percentage metric against the homotopy based planners. Second, like the results of the previous case study, HPPMs spent less time, significantly, than sampling-based planners (see Figure 20c). Third, meanwhile HPPM found one of the longest paths, MTHPPM\_VG obtained one of the three shortest; and the best compared to LazyPRM and

RRT\*. Finally, memory used by homotopy based planners is about a thousand times smaller than the sampling algorithms, this is due to the formulation of each approach.

#### *7.3. Case 3*

This case study presents an example in which some visibility points are selected to create a forward direction sequence of target-points. This is adequate to be integrated in a single-value piecewise linear formulation employed in MTHPPM, as already explained in Section 5. Figure 21a depicts the environment map for this case, contains walls, and two hundred circular obstacles. The operation of the proposed MTHPPM\_VG methodology for this configuration is similar to previous examples, furthermore, the visibility-points selection process is added. First, the approximation of the solution path is computed using the VG algorithm, then, an automatic criterion to discard nodes of the visibility path which implies backward advance in *x*-axis. Figure 21b shows nodes and paths after the visibility-point selection process is executed; the remaining nodes are employed as target-points for MTHPPM (dashed path). Figure 21c,d displays the resulting collision-free paths for HPPM against MTHPPM and between SBP algorithms (shortest path after one hundred runs for each SBP), respectively.

**Figure 21.** Comparative results for case study 3; best run of each SBP algorithm. (**a**) Path planning problem case study 3. (**b**) Visibility graph path without circular obstacles. (**c**) Paths obtained with HPPM and MTHPPM\_VG. (**d**) Paths obtained with SBP algorithms.

Figure 21 shows that the path obtained by MTHPPM\_VG is drastically shorter than the HPPM path. Figure 22 show the summarized results of HPPM, MTHPPM, and SBP algorithms for the most significant metrics.

**Figure 22.** Comparative results for case study 3. (**a**) Successful paths. (**b**) Path length. (**c**) Execution time. (**d**) Memory.

By analyzing the simulations, it is important to note: First, lazyRRT has the worst performance of all SBP algorithms concerning to percentage of feasible paths found; EST, KPIECE1, PRM, Bi-RRT, LazyPRM have one hundred percent of success. Second, using results from previous case studies, homotopy based planners spend, notoriously, the shortest time compared to sampling-based planners (see Figure 22c). Third, MTHPPM\_VG obtains a path very close to the shortest (RRT\* calculated the best). Finally, consumed memory by homotopy based planners are around one thousand times lower than SBP algorithms.

#### *7.4. Case 4*

This case study focuses on the bug trap problem, additionally, it deals with visibility nodes issue introduced in the previous case study. Figure 23a shows the configuration map with three hundred circular obstacles, walls, and narrow corridors. The initial position of the robot is inside of the bug trap, then, the planners need to plan a path that allows them to surround it, avoid the trap, and pass through two narrow corridors to reach target-point placed at the corner on the map. Similar to the previous example, for MGGHPPM, the VG algorithm found a first approximation of the path considering only the walls. Afterward, the discrimination process is executed to delete nodes of the VG path, which implies tracking in the backward direction as drawn in Figure 23b (dashed path). Finally, the full map is solved by MTHPPM using the remaining nodes of VG path as target-points. The resulting

collision-free paths for HPPM against MTHPPM and between SBP algorithms (shortest path after one hundred runs for each SBP) are presented in Figure 23c,d, respectively.

**Figure 23.** Comparative results for case study 4; best run of each SBP algorithm. (**a**) Path planning problem case study 4. (**b**) Visibility graph path without circular obstacles. (**c**) Paths obtained with HPPM and MTHPPM\_VG. (**d**) Paths obtained with SBP algorithms.

Figure 23c shows the ability of HPPM and MGHPPM\_VG to deal with a bug trap and narrow corridors. The tendency of HPPM to follow a direct path can be observed in Figure 23c. On the other hand, it also depicts an enhanced homotopy path since MTHPPM uses the points provided by the VG algorithm. Figure 24 shows a condensed performance information of the HPPM, MTHPPM, and SBP algorithms for time consumption, memory, percentage of feasible paths found, and length of the path.

**Figure 24.** Comparative results for case study 4. (**a**) Successful paths. (**b**) Path length. (**c**) Execution time. (**d**) Memory.

Analyzing the obtained results, can be summarized as follows: First, RRT and LazyRRT exhibit poor performance to deal with bug traps; unlike EST, KPIECE1, PRM, HPPM, and MTHPPM\_VG which have a one hundred percentage of success (see Figure 24a). Second, Figure 24c shows that homotopy based planners spent, noticeably, less execution time than sampling-based planners. Third, the HPPM creates a path with competitive length, while MTHPPM\_VG obtained a path close to the shortest (RRT\* calculated the best). For this metric, the results of RRT because the length of its best path corresponded to an unsuccessful run are discarded. Finally, homotopy planners exhibited execution times around three orders of magnitude lower than the best SBP algorithms.

#### *7.5. Case 5*

The utility of MTHPPM to solve the problem of multiple-target path planning in applications like pick-and-delivery activities, museum guidance robot, rescue task robot, and so forth, is presented. For this case study, target-points are previously selected by the user according to a sequence of rooms that the robot must visit. Figure 25a shows the environment with a sequence of fifteen milestones, including the start-position (point 1) and target-position (point 15). This map contains five hundred circular obstacles representing the configuration of people (considered as static obstacles for a given instant of time). The main objective of this case study is to show the ability of MTHPPM to calculate a collision-free path for multiple-targets in a map filled with obstacles. For this case, each SBP algorithm was executed one hundred times and applied the point-to-point strategy to generate a path to visit all fifteen target-points. Figure 25 depicts the solution paths of six SBP algorithms generated using the point-to-point technique. These figures show the best run (length of path) for each SBP algorithm between two points after one hundred runs.

**Figure 25.** Multiple-target path problem solution using SBP algorithms with point-to-point strategy. (**a**) Multiple-target point problem (case study 5). (**b**) Multiple-target path using EST and KPIECE1. (**c**) Multiple-target path using RRTConnect and PRM. (**d**) Multiple-target path using RRT and RRT\*.

Figure 25b shows that EST and KPIECE1 planners are not capable to solve the problem of the narrow corridor between target-points seven and eight. Here, it presents a visual comparison of full path length because benchmark results provided a similar outcome regarding memory consumption, time, and success percentage compared to previous experiments. It is important to remark that, the full path of each SBP algorithm is integrated by the shortest segments between two points obtained after one hundred executions. Figure 26a presents the single value PWL - function employed to calculate two multiple-target paths through MTHPPM. The fist path is calculated for the motion in the forward direction (path from point 1 to 8) and the second in the backward direction (path from point 8 to 15), as can be seen in Figure 26b.

**Figure 26.** Performance of MTHPPM for office-like environments. (**a**) Piecewise-linear representation of direct path for multiple-target approach. (**b**) Successful solution path.

Figure 26b shows the ability of MTHPPM to obtain a free collision path for a sequence of target-points. This planner calculates the full path in only two runs, one path in the forward and the other one in the backward direction. Although the path obtained by MTHPPM is not optimal, this case study shows the application of the proposed MGHPPM to solve multiple-target problems with hundreds of obstacles and narrow corridors. This characteristic is desirable for monitoring robots in oil platforms and forests or natural reserves.

#### **8. Pure-Pursuit Controller Using Matlab Robotics System Toolbox**

This section provides two simulations to validate the feasibility of paths obtained by HPPM, MTHPPM, and the pure-pursuit algorithm. The pure-pursuit algorithm procedure computes the angular velocity of the robot from the current position to reach a look-ahead point (see in Figure 27). This is not a common controller, although it is suitable to be employed for path tracking purposes [59,60]. Here, two tracking simulations of a waypoints set obtained by homotopy planners (HPPM and MTHPPM) are used. The tracking is simulated by the pure-pursuit class contained in MATLAB Robotics System Toolbox and the model of a differential drive robot. It is assumed that linear velocity is constant and angular velocity changes according to the instantaneous center of curvature (ICC). The ICC is calculated from the look-ahead parameter as explained in References [59,60]. Pure-pursuit operates as follows: First, the desired path is obtained by any path planning algorithm, a homotopy planner for this case. The path should be represented as a finite array of *n*-waypoints in the form (*xi*, *yi*), *i* = 1, 2, 3, ... , *n*. Second, using the look-ahead distance (value previously set), the algorithm calculates ICC and then the angular velocity to move the robot from one waypoint to another until the last point of the path is reached.

**Figure 27.** Pure-pursuit scheme.

Figure 27 presents the basic operation of the pure-pursuit controller. It can be noticed that the look-ahead distance shown in this figure is shorter than the distance between two waypoints. For these cases, the algorithm completes the desired path using segments of straight lines. For the next simulations, parameters of pure-pursuit controller are set as follows


In this example a configuration similar to case studies 3 and 4 in Section 7 is applied. Before the homotopy based planner is executed, the dimensions of the robot must be added to all obstacles, the robot dimensions are considered, as shown in Figures 28a and 29a. In these, the differential drive robot is modeled by a disc, then, its radius is added to dimensions of all obstacles as guard distance (dash line in Figures 28a and 29a), as explained in Reference [10]. The complete environment map, considering dimensions of the robot, is normalized to be introduced in HPPM. Figures 28b and 29b show the path followed by differential robot employing the pure-pursuit algorithm. These figures show that the difference between waypoints (homotopy path) and pure-pursuit path is minimal.

**Figure 28.** Pure-pursuit path tracking for the solution path of case study 3. (**a**) Floor plan considering dimensions of robot. (**b**) Path traced by differential drive robot.

**Figure 29.** Pure-pursuit path tracking for the solution path of case study 4. (**a**) Floor plan considering dimensions of the robot. (**b**) Path traced by differential drive robot.

For the example of Figure 28, the initial point is placed at (0, 0) and final point at (14.5 m, 14.5 m). Meanwhile, for the example of Figure 29, the initial point is placed at (0.5 m, 0.5 m) and final point at (14.5 m, 14.5 m).The apparent radius (dash line) of all circular obstacles is 0.26 m, which represents the sum of the obstacle radius (0.01 m) and the robot radius. The capability of HPPM to find solution of maps with narrow corridors can be observed in the example of Figure 29b.

Figures 28 and 29 validate that paths obtained by the proposed homotopy based planner (MGHPPM) can be applied to differential drive robots through the pure-pursuit controller. These simulations also denote that homotopy paths do not require post-processing or smoothed unlike SBP algorithms.

#### **9. Discussion**

This work presents a novel planner to obtain a collision-free path from a sequence of target-points which contains several improvements in the numerical implementation of the original homotopic path planner [19]. The main contributions of this work are outlined in four aspects.



**Table 2.** Number of reversal effects and non-convergences solved in SA for HPPM and MTHPPM, case studies 1–4, and using DST technique.


Table 3 shows the summarized results of case studies 1–4 for the maximum, minimum, and median values of performance metrics presented in its respective box diagrams. It is important to note that minimum length (dimensionless because map is normalized) for some SBP planners denotes the result of an unsuccessful run. Furthermore, this table shows the global execution time for all SBP algorithms, which considers the execution time of the algorithm and time spent to smooth the path (path simplification process). If the median value of each metric for every SBP planner is taken, then, it is possible to remark the following conclusions for the case studies: (I) The HPPM and MTHPPM have the best performance for memory consumption, about one-thousandth, compared to the amount spent by SBP algorithms; for instance, lazzyRRT and LazzyPRM for the second case study. (II) The execution of homotopy based planners is between ten and one hundred times faster than SBP algorithms. (III) The homotopy based planners presented in this paper have an adequate balance between the measured performance metrics; while for SBP algorithms, execution time and length of the path are inversely proportional parameters. In other words, faster planners (EST, KPIECE1,

and RRT) find longer paths and shortest path calculations consume the maximum amount of time (RRT\*). Results also show that the proposed methodology does not have any dependence on libraries, packages, or external functions because complex procedures are treated using symbolical analysis as it is presented in Section 5.5. Therefore, its implementation is simple and cheap in terms of computational resources. This allows an easy implementation in multiple platforms like embedded systems with microcontrollers, microprocessors up to implementations in PC's and FPGA's. In addition to these advantages, simulations presented in Section 8 show that the path obtained through homotopy based planners can be easily followed by a differential drive robot using the pure-pursuit algorithm.


**Table 3.** Summarized performance results for case studies 1–4. Metrics: memory consumption (M), execution time (T), and length of path (L).

Figures 30–33 show the minimum, median and maximum results for execution time, length of path, and memory of cases 1–4 (in number of times), according to the results of MGHPPM (MTHPPM\_VG). These figures represent a visual interpretation of Table 3 which helps to denote the advantages of the MGHPPM where, a) the memory used in all case studies by the SBP algorithms is between hundreds and thousands of times greater than MGHPPM; b) the path obtained through MGHPPM is very close

to the shortest path found by SBP algorithms, and c) the execution time of MGHPPM for all case studies is between five and one hundred of times smaller than every SBP algorithms.

**Figure 30.** Bar graphs for time, length of path, and memory in number of times respect to the MGHPPM results of case 1. (**a**) Execution time. (**b**) Memory. (**c**) Length of path.

**Figure 31.** Bar graphs for time, length of path, and memory in number of times respect to the MGHPPM results of case 2. (**a**) Execution time. (**b**) Memory. (**c**) Length of path.

**Figure 32.** Bar graphs for time, length of path, and memory in number of times respect to the MGHPPM results of case 3. (**a**) Execution time. (**b**) Memory. (**c**) Length of path.

**Figure 33.** Bar graphs for time, length of path, and memory in number of times respect to the MGHPPM results of case 4. (**a**) Execution time. (**b**) Memory. (**c**) Length of path.

#### **10. Conclusions and Future Work**

In this work, a path planning method with multiple-target applications has been presented, it is capable of solving complex maps of hundreds of obstacles. This method contemplates in its computational core a series of novel and effective tools such as Double Spherical Tracking, the integration of a dummy obstacle to improve performance and reduce computing time and the number of iterations. As well as, a scheme of multiple solutions of a NAES formulated by approximations of PWL functions. In the same way, two new schemes are proposed to avoid discontinuities in the spherical tracking using patches and bridges in the plane. Furthermore, a simple and efficient solution is proposed as a criterion for the automatic selection of the repulsion parameter of the obstacles representation for the homotopy based planner methods.

From the case studies it can be denoted first, the proposed planner (MTHPPM) is faster; between five and one hundred times than the average of one hundred runs for each SBP algorithms. Second, the calculated path is very close to the shortest path; the difference is between ten and twenty percent. Third, the success rate is one hundred percent for all case studies, while some SBP algorithms achieved around thirty-five percent of the failure rate. Fourth, the proposed planner uses just one-thousandth of the memory than the best SBP algorithms employs for every case study. Fifth, the proposed methodology does not have any dependence on libraries, packages, or external functions for which its implementation is simple as well as cheap in terms of memory and computational resources. Therefore, these characteristics allow its implementation for real-time applications in multiple platforms from embedded systems with microcontrollers of low resources, such as it is presented in Reference [10], and microprocessors to PCs and FPGAs, as it is presented in Reference [21].

As future work, it is left for further research on use of multivalued piecewise linear function formulation, presented in Reference [51], which could generate two research lines. One for multiple-target planning of robots with displacement in 3-D or more dimensions like: drones, underwater robots, computer animation, robot manipulators, and molecular simulations [3,8,46]. The other one, employing multivalued or parametrized piece-wise linear functions like the ones presented in Reference [61,62], could enhance the performance of MGHHPP\_VG in terms of path length, as it can be deducted from the cases studies, which would allow obtaining paths using visibility points in forward and backward movement on both axis. Besides, further work is needed to extend HPPM and MTHPPM to handle mechanical restrictions of non-holonomic robots, multi-agent systems, and path planning for dynamic environments with uncertainties. Finally, considering the advantages of the proposed homotopy based planner over SBP algorithms in terms of path length (close to the shortest path found by the best run of RRT\*), execution time (between five and one hundred times faster than the median of SBP algorithms), and memory consumption (about a thousandth of that used by the SBP algorithms) makes of this a good choice to be implemented in practical applications. Besides, the linear relation between the number of obstacles and the complexity (in terms of memory and execution time) for homotopic planners, as it was presented in References [10,21] against the exponential complexity increase of SBP algorithms [2,3,8,9,63–65] allows to conclude that MTHPPM is the best option to planning collision-free robot motion paths. Especially, when it is needed to be implemented in a system with limited resources (like embedded systems) or required to solve complex environments where time constraints are tight.

**Author Contributions:** Investigation: G.D.-A. and H.V.-L.; Supervision: L.H.-M. and J.H.-C.; Methodology: V.M.J.-F.; Validation; H.D.C.-C.; Resources: A.H.-J., R.C.A. and S.H.-M.; writing—original draft preparation: G.U.D.-A.; writing—review and editing: L.H.-M. and H.V.-L.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
