Pathfinding algorithm based on rotated block AOR technique in structured environment

1.   Introduction

The surge of mobile robot navigation technology nowadays is rising by virtue of human desire. It is not only limited to transportation, but it is also very useful in other problems such as manufacturing, industrial, and so on. In general, the ability of a robot to navigate autonomously is a prerequisite and foundation for the development of intelligent robots. Pathfinding algorithms must be both efficient and practical to ensure the proper performance of mobile robots' various tasks.

Figure 1 illustrates the key phases related to the operation of a robot. The problem of robot pathfinding can be divided into four categories:

● Localization: Where do I go?

● Pathfinding: Where was I supposed to go? Which is the best path to be reached?

● Motion control: How can I move?

● Cognitive mapping: Where was I? What should I remember?

There is no doubt that pathfinding is a major issue in the navigation process. The pathfinding algorithm allows the robot to choose and identify the right path within the configuration space. The two major components of global path planning are the robot's representation of the world in configuration space and its ability to execute the algorithm. These two components are closely linked and have a significant impact on one another in choosing the optimum path and time for the robot to move around in the configuration space ^[2].

One of most difficult problem in robot applications is the construction of robust autonomous motion planning. To create a truly autonomous mobile robot, it must be able to compute its trajectory between locations quickly and accurately while also avoiding collisions with surrounding objects. The goal of this study is to use harmonic potential values to simulate point-robot pathfinding in a specified area using heat transfer theory. Laplace's equation ^[3] is used to model the heat transfer problem. Harmonic functions are the solutions to Laplace's equation, and they represent temperature values in the configuration space that will be exploited to simulate path generation.

Iterative methods are another way of generating harmonic functions by numerically solving the Laplace equation under Dirichlet or Neumann boundary conditions and discretizing the two-dimensional environment into a grid of discrete points ^[4,5,6,7,8]. Harmonic functions obey the min-max principle, thus there can be no spontaneous creation of local minima within the solution region. As a result, the only critical points that can occur are saddle points. A search in the neighbourhood of such a critical point would then reveal the exit from that point. Furthermore, any deviation or disturbance of route from such points results in a smooth path everywhere. The discretization of the Laplace equation yields a sparse linear system, in which the computing costs of iterative approaches are generally significantly greater than the analytical technique. The computational complexity is determined by the discretization resolution; hence, if a greater region of obstacles is occupied, fewer calculations and storage are required. Laplace's solutions can be numerically derived using the finite difference method. Laplace' equation is transformed into a linear system, and most of the matrix elements are zero when written in matrix notation. The resultant matrix is often extremely large and sparse, and typically requires very large computing resources. This study thereby solved the linear system by means of the iterative approach, since it avoids the need to store a very large matrix in memory.

For pathfinding method in this study, the Gradient Descent Search (GDS) technique utilizes the Laplacian potential values obtained using the iterative methods discussed in the next section as guidance during the path construction. From the current point, the GDS method examines the potential values of its eight neighbouring points and simply picks the node with the lowest Laplacian potential value. This procedure is repeated until it reaches the target point ^[4,9,10,11].

Recently, a real-time path planning method using harmonic functions based on the Proper Generalized Decomposition (PGD) has been suggested ^[12,13]. The proposed method was tested on two structured environments: an L-shaped corridor and a bug trap planning environment. The experiments showed encouraging results, in which a physical Lego Mindstorms robot was tested to validate the results. The experiments, however, were only conducted in a small and relatively simple environment. It was also stated that the offline phase of the PGD-based framework on a square environment of 150 x 150 nodes would take 14 hours using a high-performance computer to obtain the solution for every start and goal vehicle position. However, the online phase was faster and required minimal processing time.

Most path planning methods in the literature that are based on the harmonic potentials rely on the used of the standard Laplacian operator. This study presents several to evaluate the capability of a rotated block Laplacian operator, also known as Explicit Decoupled Group (EDG), in the computation of harmonic potentials. A linear system of rotating approximation equations is developed and computed via EDG Accelerated Over Relaxation (EDGAOR) by employing the rotated block operator. The results of the implemented EDGAOR algorithm are compared and analysed in terms of the number of iterations and CPU time against the existing approaches. The findings showed that the EDGAOR outperformed the existing methods.

2.   Pathfinding technique

2.1.   Global pathfinding

Global planning stores the map of its environmental and then utilises that map to create a feasible path. Accordingly, global pathfinding seeks to search for a collision-free route by which the robot may move from its original position toward its intended position without colliding with any obstacles. Figure 2 illustrates the global pathfinding framework. A model of the robot workspace that includes position of the robot, obstacles and vacant spaces, must first be made. The vacant area in the workspace is then discretized to form a graph representing the connectivity of the space. The path finding algorithm is then employed to determine a feasible route for the robot to trace using the graph of connectivity in order to reach the goal.

The process of pathfinding construction in this study comprises the following steps:

Begin

Step 1. Create the map of the robot's environment.

Step 2. Establish the formulation and modelling of the iterative schemes using finite difference method.

Step 3. Develop and implement the algorithms of the proposed iterative schemes.

Step 4. Run the numerical simulations to obtain the solutions.

Step 5. Conduct the performances evaluation and complexity analysis of the algorithms.

End

2.2.   Harmonic function

Let Laplace's equation below, be satisfied by a harmonic function in a region $\Omega \subset {\Re ^n}$

with ${x_i}$ as the i-th coordinates in the Cartesian and the dimension is $n$. For the robot path creation, the region boundary is made up of walls, obstacle boundaries, and the target location. Since the harmonic potential complies with the min-max criterion ^[4] in the specified domain, it does not have local minima except the target position, and generally produces smooth path, thus making harmonic potential approach a feasible and very attractive option for robot navigation. The Laplace equation is commonly solved using standard approaches ^[14,15,16]. In this study, Eq (1) was solved using fast rotated block iterative technique to speed up the computational execution.

A global approach is used to compute harmonic potentials of the robot workspace by employing the analogy of heat distribution in areas in which all obstacles and boundaries (with the highest potential value) represent heat sources and target position (with the highest potential value) as heat sink. Laplace's equation is used to model this heat transfer activity and is solved numerically to obtain the heat distribution that represents the harmonic potential for each node in the graph. Utilizing the property of heat distribution that flow from high temperature to lower temperature, a gradient search can be used to generate path from any start point with high potential value to target position with the lowest potential value. The pathfinding algorithm utilizes the gradient descent search to determine a feasible route for the robot to navigate the workspace safely from the start point to reach the goal position.

This research seeks to replicate the aforementioned paradigm for pathfinding, describing the solution of Laplace's Eq (1) through the analogy of temperature (for the potential) and heat flow (for the pathway). The experiment is conducted with a two-dimensional domain with various forms of obstacles, including the inner and external walls. To solve Eq (1) and obtain the potential values for each node, the EDGAOR iterative approach is utilised. For performance comparison purposes, the existing Explicit Group Successive Over Relaxation (EGSOR), Explicit Group AOR (EGAOR), and EDGSOR were also investigated.

3.   Materials and methods

3.1.   Rotated block technique

Abdullah ^[17] was the first to present the EDG iterative approach that employed block half-sweep technique on rotated grid to solve the two-dimensional Poisson equation. This approach is also used in ^{[18,19,20,21,22,23,24,25,26]} to solve partial differential equations. A modified version of this approach was investigated for the solution of the diffusion equation ^[24]. Meanwhile, ^[25] discussed early efforts on integrating SOR with other techniques. In the robotics literature, the classic GS ^[9] and SOR ^[3,26] algorithms have been utilised to solve Eq (1). To obtain the solutions of Laplace's Eq (1), this work developed a robust numerical solver using rotated block iteration approach that employ half-sweep technique.

In general, the EDG technique computes only half of all points in the region, resulting in a significant reduction in total processing time. The application of EDG iteration, which examines the interior grid points (black grid points), is compared to the conventional full-sweep (alternative name is Explicit Group, EG) technique, which iterates every grid point in the region (see Figure 3). The half-sweep iteration technique was first used in robotics in ^[18]. In this study, the EDG iteration is essentially subjected to the rotated 5-points finite difference approximation (5-FDA) equation to solve Eq (1). The key advantage of this iterative technique is that it reduces the computing complexity by evaluating only a portion of the total number of node grid points. Figure 4 depicts the computational molecules of the iterative techniques of full-sweep (FS) and half-sweep (HS).

All the algorithms in the Block iterative techniques family indicate that acquiring four Laplacian potentials per calculation speeds up the evaluation. In the event of computing groups of points that are close to the boundary, it will be managed as a set of a single point and two points, as illustrated in Figure 3. After all, for Block iterative techniques, the approximate values of the residual node points are computed directly by using direct methods ^[17,27]. It should be emphasised that all algorithms in this study are executed so long as the results meet a predefined convergence criteria denoted as $\varepsilon$. The stopping criterion in this study is based on $\left\| {{u^{\left({k + 1} \right)}} - {u^{\left(k \right)}}} \right\| \leqslant \varepsilon$.

In general, the executions of all iterative schemes are imposed on black node points in Figure 3 until the convergence test criterion is satisfied. As depicted in Figure 3 (left), all formulations employing the EG iterative techniques compute a group of four nodes at once during the iteration process (except for certain cases near the boundary). Meanwhile, the EDG is derived from a rotated 5-FDA ^[17]. As shown in Figure 3 (right), the solution domain for the EDG method is divided into two types of node points, $\cdot$ and $\circ$. The solutions of each node point group can be executed by pairing the points of an equal type and can be carried out independently for each of the pairings. Due to this independency, almost half of the iteration over the solution domain is performed in either type of the points, hence saves the performance time. Furthermore, the direct method ^[28] is computed via Eq (3), as demonstrated in the next section, to compute sets that are close to the border, i.e., with only one or two points.

3.2.   5-point finite difference approximation

Consider the 2-dimensional equation of Laplace in Eq (1), which is rewritten as

The Laplacian operator is denoted as ${\nabla ^2}$. To solve Eq (2), the system needs to be discretised via finite difference method before it can be computed efficiently using the numerical method. The simplest 5-FDA for the 2D Laplacian, ${\nabla ^2}$, is given as

where $U$ is a function to satisfy the equation of Laplace, $u$ is the potential node at point $\left({x, y} \right)$, and $h$ is the distance between node points for each direction. The above equation can be termed as a second-order equation since its error term is of order two, ${h^2}$. By using the central difference formula, it will then be deduced using Taylor series approximation into the second derivative, and be expressed as follows (which represents FS iteration case)

Meanwhile, the approximation that is based on the cross-orientation operator is formed upon rotating the $x - y$ plane by ${45^ \circ }$ in the clockwise direction ^[29,30]. This produces the rotated (skewed) 5-point approximation, also known as HS iteration, denoted as

This also can be shown in matrix form as

and

for FS and HS iteration cases respectively.

By letting ${U_{ij}}$ approaches the solution of u for Laplace's Eq (2) along the grid point $\left({{x_i}, {y_j}} \right)$, the discretization of these Laplace equation via conventional 5-point stencil is

for the FS case. Meanwhile, for HS iteration holds that

Then, by employing Eqs (3) and (4) to the 2-dimensional Laplace's problem based on Eq (2), it will eventually generate a linear system. This system can be formed in matrix notation and produce large and sparse matrices. The aforementioned linear system is displayed as

The coefficient $A$ and $b$ are both identified, where $A$ is in matrix form and $b$ is a vector. Whereas coefficient u is an unknown vector. For the FS iteration case, the matrix $A$ is said to be a block triangle matrix with order $\left({N - 1} \right) \times \left({N - 1} \right)$, and signified as

From matrix $A$, each element of block $T$ is represented as

while $I$ is an identity matrix with the same order. Meanwhile, coefficients u and $b$ from Eq (5) may be defined as

where $k = 1, 2, \ldots, N - 2, N - 1$.

Meanwhile, for the HS iteration case, matrix $A$ from the linear system Eq (5) can be described as

where

Vectors u and $b$ from Eq (5) for the HS iteration case are defined respectively as

with $i = 1, 3, 5, \ldots, N - 3, N - 1$ and $j = 2, 4, 6, \ldots, N - 4, N - 2$.

As stated previously, the system of linear Eq (5) generates a sparse and large matrix. Consequently, the iterative methods can be solved using either Point or Block iterative techniques because the proposed methods are ideal for solving this problem ^[17]. Hence, considering the finite difference of Eqs (3) and (4), the Gauss-Seidel iterative schemes for FS and HS cases can be rewritten and denoted respectively as follows,

As previously stated, this study concentrates on the family of Block iterative schemes for the proposed solvers, i.e., SOR and AOR, using EG and EDG iterative methods. With the intention of implementing the Block scheme, the Point iteration scheme for corresponding methods is required. The FSSOR method is, in fact, the standard SOR method that applies traditional FS iteration via a 5-point discretization scheme ^[31]. The iterative scheme for the Point SOR method is given below. It is referring to Eq (6), by inserting a weighted parameter, $\omega$ via SOR ^[31].

Meanwhile, the HSSOR iterative method is formulated by employing the HS concept ^[17]. From Eq (7), the Point rotated SOR iterative scheme is given as

The difference between Eqs (6) and (7) to Eqs (8) and (9) are clearly on the addition of weighted parameters $\omega$. The basic idea of the SOR iterative method comes from Gauss-Seidel iterative schemes with weighted parameters ^[32,33]. Note that if the value of $\omega$ is equal to 1, then the SOR scheme is effectively simplified to the standard Gauss-Seidel scheme.

The AOR method belongs to the family of over-relaxation methods. It is a generalisation of the SOR method with two parameters. For this study, these two relaxation parameters are denoted as $r$ and $\omega$. Both parameters can be utilised to generate iterative methods that can speed up the convergence rates, and AOR are more flexible and suitable than any other methods in this family. To perform the AOR iterative scheme proposed in ^[34], replacing $u_{i - 1, j - 1}^{\left({k + 1} \right)}$ and $u_{i + 1, j - 1}^{\left({k + 1} \right)}$ with $u_{i - 1, j - 1}^{\left(k \right)}$ and $u_{i + 1, j - 1}^{\left(k \right)}$, respectively, as well as appending the nodes of $\frac{{r\left({u_{i - 1, j - 1}^{\left({k + 1} \right)} - u_{i - 1, j - 1}^{\left(k \right)}} \right)}}{4}$ and $\frac{{r\left({u_{i + 1, j - 1}^{\left({k + 1} \right)} - u_{i + 1, j - 1}^{\left(k \right)}} \right)}}{4}$ into Eq (8) are required. Thus, the Point AOR iterative scheme is

Whereas, based on the explanation above and implementing every replacement of terms into Eq (9), the rotated Point AOR schemes are stated as

Based on Laplace Eq (2), all of the proposed iterative techniques simply replace each node's value from the average of four neighbours of the corresponding node. The node values that indicate the boundaries, obstacles, and target position in this study remains constant.

3.3.   Explicit group accelerated over relaxation scheme

To facilitate the formulation of the Block EGSOR method, observe a set of four node points $\left({4 \times 4} \right)$ in Figure 3 (left). By examining Eqs (3) and (8), the scheme generally can be expressed as ^[17,28,30].

where

Apparently, Eq (12) may also be converted into the form of a system of linear Eq (5). It will then be translated upon establishing the inverse of a matrix to coefficient $A$, as shown below

where

Now, the Block EGSOR scheme for Eq (13) is indicated as

Meanwhile, the formulation of the Block EGAOR method also considers the node points in Figure 3 (left), but assessing the approximation equation from Eq (3) and (10). Then, the AOR block scheme of AOR is stated as ^[35]

with

Again, Eq (15) may also be transformed into a system of linear Eq (5) and translated as below, by initiating the inverse of matrix $A$

with

Finally, the Block EGAOR iterative scheme for Eq (16) is presented as

3.4.   Explicit decoupled group accelerated over relaxation scheme

The EDG approximation equation is on the basis of the cross-orientation operator, which allows it to be formed upon rotating or turning the $i - j$ plane by ${45^ \circ }$ in clockwise direction. Now, to simplify the formulation of the Block EDGSOR scheme, a set of four node points (see Figure 3 (right)) is examined in order to construct a $\left({4 \times 4} \right)$ algebraic linear system ^[17,27,36], expressed by Eq (15), where

and

The linear system equation in matrix form as in Eq (15), can be reduced to Eqs (18) and (19) and solved individually using two $\left({2 \times 2} \right)$ linear algebraic equations as follows

The Block EDGSOR scheme for both Eqs (20) and (21) can be defined independently as shown below, respectively,

The algorithm for Block EDGSOR iterative scheme can be implemented by using either Eq (22) or (23). Both equations can be used to obtain solution.

Applying similar approach to EGAOR, by using the rotated AOR formula (7) to the Block scheme gives rise to a new formula called EDGAOR ^[37]. As a result, the AOR scheme for the EDG iteration is stated as

where

Contrasting from SOR, no generic formula exists for determining the optimum $r$ and $\omega$ values that result in the minimal iteration counts. According to ^[34], the $r$ value is often set to be close to the SOR $\omega$ value, and the numerical analysis is then repeated for a certain range of $\omega$, with $1 \leqslant \omega \leqslant 2$. As some values do not converge in certain cases, $r$ and $\omega$ values for each full- sweep and half-sweep cases are different in order to obtain the optimal value. The effect of complexity on discovering parameter values to overall computation does not change, since the value of each parameter is set before the execution/computation. The optimal values used throughout the experiments are listed in Table 1.

The implementation of Block EDGAOR scheme based on Eq (24) into solving the two-dimensional Laplace's problem, as expressed in Eq (2), is stated in Algorithm 1.

4.   Results and discussion

The implementation of the EDG iterative method computes approximately half of the group of node points in the block scheme (see Figure 3) during the iteration process. Therefore, it will drastically reduce the computational complexity, i.e., approximately 50%. With the purpose of measuring the performance of these algorithms, the simulations are conducted on four configuration regions in three different sizes, i.e., $300 \times 300$, $600 \times 600$, and $900 \times 900$. Their performances are examined using three criteria:

(i). Number of iteration count,

(ii). CPU time, and

(iii). Success of pathfinding.

The number of iterations and CPU time required by each algorithm is captured when the tolerance rate of the iteration process is met. The computed Laplacian potentials are then used by the GDS algorithm to ensure that a path from the initial to the target position can be successfully created.

A variety of obstacles of various shapes are placed in the given region. In the preliminary configuration, the Dirichlet boundary condition was used, where the boundaries and obstacles are fixed at high-temperature values. The target position was placed at the lowest temperature reading, whereas no initial value was assigned at the start position. Other grid nodes remained at zero temperature. The execution had been performed using a computer, operating with a 2.50GHz processor and 8GB of RAM. The iterative procedure of numerical computation continued until the stopping condition was satisfied. The loop is supposed to be halted if potential values have not changed further and the changes within Laplacian potentials at the computation of current (k) and next (k+1) iterations are very small, i.e., ${1.0^{ - 10}}$. This level of precision was required to avoid flat areas in the solutions, often known as saddle points, which means the sudden non-existence of descending gradient patterns. As a result, the algorithm simply stops searching and gets stuck in flat areas at the lowest point. Tables 2 and 3 indicate the number of iterations along with CPU time (in seconds) for each numerical method compared in the tests. Clearly, the EDGAOR approach showed to be faster than other approaches tested. On average, the EDGAOR successfully reduced the number of iterations compared to EDGSOR by 19%, approximately 24% for EGAOR, and 33% for EGSOR. In terms of CPU time reduction, the EDGAOR averagely reduces the execution by 16%, 36%, and 40% against EDGSOR, EGAOR, and EGSOR respectively.

After determining the configuration region's potential values, a smooth path is formed by following the potentials distribution using the gradient descent method, in which the algorithm follows the descending gradient from the initial position to the next sequential points at lower potentials until it reaches the lowest potential (the target) point. Figure 5 depicts the pathways that were efficiently generated in a given stationary environment using the numerically computed potentials distribution. Every starting location (green/square point) reached the designated target location (red/circle point), avoiding every obstacle placed in the environments. The solution can be calculated and repeated as often as desired in a stationary situation where the target point and obstacles are set. It must be recalculated only if the obstacles or the target position are changed. The trajectory of the paths can be very quick, as they involve only the gradient evaluation of precomputed Laplacian potential ^[9]. The flow diagram of the pathfinding technique of this study is shown in Figure 6.

The performance graphs in terms of iteration numbers (see Figure 7) and CPU time (see Figure 8) are also illustrated. By referring to both graphs, it can be seen that the EDGAOR had outperformed their corresponding methods, either in terms of iteration counts or CPU time. The EDGAOR method gave the best performance. This idea can be clearly seen in Tables 2 and 3. As we can see from the table of results, the graphs for the number of iterations as well as CPU time, shows the same pattern.

5.   Conclusions

This study demonstrated that the paths can be obtained from the GDS algorithm by utilising the Laplacian potentials computed using the presented EDG iterative schemes. After succeeding in obtaining the potential values present in the environment, the path of the navigating robot can be constructed. It follows the principle of heat distribution by using the steepest descent method. The algorithms follow the temperature descent flow to the sink, which in this case is the goal point, thus providing a path line for the navigation. The path line can be generated using different starting and goal points by computing the Laplacian potentials numerically. The generated path line had successfully avoided any obstacles while navigating from the start green point to the target red point.

Harmonic functions offer a rapid approach to generating paths in a structured environment. The approach avoids the local minimum problem at some cost in computational speed for very large environments. However, given that the proposed solver is well suited for parallel processing, it can be applied to quickly solve Laplace's equation. The authors feel that this is not a serious drawback. Moreover, most inner and outer boundaries of an indoor structured environment remain static most of the time, so only the regions occupied by moving obstacles and people need to be re-computed. The computation of harmonic functions for these dynamic regions can be completed very quickly. In the future, the suggested Block iteration procedure such as EDGAOR can be extended using the log-space approach as demonstrated in ^[38], where the parallel implementation of the proposed approach showed great performance improvement.

The experiments conducted in this research revealed that numerical computation of Laplace's equation to solve the problem of robot pathfinding is notably appealing and practical due to recent developments and newly discovered numerical techniques, as well as the current availability of fast machines. The EDGAOR iterative approach proved to be much faster than the previous SOR and AOR approaches. An increase in the number of obstacles has no negative effect on performance; in fact, the calculation becomes faster considering the regions occupied by obstacles are ignored throughout the computation. In the future work, studies on half-sweep and quarter-sweep ^{[36,39,40,41]} iterations, would be explored in order to further improve the rate of convergence of the proposed approach.

Acknowledgments

The authors acknowledge the National Defence University of Malaysia for the funding of this article.

Conflict of interest

The researchers declare that there is no conflict of interest regarding the publication of this study.