An iterative technique for solving path planning in identified environments by using a skewed block accelerated algorithm

1.   Introduction

Path planning is a computational problem to find a sequence of valid configurations that move the object of interest from an arbitrary source to a specific destination. The term path planning has been used widely, mostly in the robotics area. Good path planning for a mobile robot means it should be executed while avoiding walls and dodging any obstacles in between. Autonomous robot path planning is very impactful, as it is extremely beneficial not only for transportation, but also for households, pharmaceuticals, industries, manufacturing, aerospace and others. In general, the ability of a robot to navigate autonomously is a prerequisite and foundation for the development of intelligent robots. To ensure the successful completion of the many tasks that mobile robots are designed to accomplish, path-planning algorithms must be both effective and practical.

On the whole, navigation problems are divided into four main categories, i.e., localization, path planning, motion control and cognitive mapping ^[1]. In essence, the localization of a robot is telling where the robot is currently located, whereas path planning reveals in which direction the robot is meant to travel and what is the best route to arrive at the destination. Robot motion control gives instructions on the manner in which the robot is able to move. Cognitive mapping concerns the previous, current and next positions of the robot, as well as to what extent the robot should remember. Hence, from all of these category's functions, it is not wrong to say that path planning is a crucial issue to be solved in achieving intelligent robots. With the aid of a path-planning algorithm, the robot should be able to select and recognize the ideal route inside the given configuration region.

It is common knowledge that there are two types of navigation, i.e., global and local approaches. This work considers global navigation, where prior knowledge of the configuration region is given, which is also known as an off-line mode for path planning. The robot's ability to represents the real world and to execute the algorithm are the two key elements of global path planning. These two elements are interconnected and significantly influence one another while determining the optimal route and minimum duration for the robot to travel across the known environment ^[2]. The main objective of this study is to utilize the harmonic potential in simulating point-robot path-planning in the identified environment following the resemblance of heat transfer. The model of the aforementioned heat transfer problem was designed by using Laplace's equation ^[3]. Harmonic functions are literally the solutions to Laplace's equation, which are also called Laplacian potentials. These Laplacian potentials portray the analogy of temperature values from heat transfer in the configuration region that will be used to simulate the path generation. Heat transfer's capacity to surpass the local minima enigma is one of its most crucial characteristics, making it very encouraging for robot navigation control. The following sections provide further explanations. Numerous approaches can be used to achieve harmonic functions, but the most prevalent approaches rely on numerical techniques due to the availability of rapid processing machines, as well as on their elegance and competence in solving the problem ^[4,5,6].

Furthermore, the iterative approach could be implemented to produce the Laplacian potential by computationally solving Laplace's equation and discretizing the environment through the use of mesh grid points ^[7,8,9,10]. The cost of computation for iterative techniques is substantially higher than analytical techniques since the discretization process produces a sparse linear system. However, the discretization resolution determines the computational complexity. Using a finite difference approximation, the Laplace equation can be translated to a linear system. When this linear system is documented in matrix notation, the majority of the elements are zero, and the matrix size of the end product is usually large and sparse. Therefore, the iterative technique is employed to solve the linear system to prevent the higher requirement for memory storage of those large and sparse matrices. Several experiments were executed in this work to measure the competency of the skewed block accelerated iterative approach in terms of the computation of the Laplacian potential when generating robot-point paths. Fundamentally, a complete process of the path-planning construction phase in this study consists of the following steps:

Begin

Step 1: Mapping of the robot's configuration region.

Step 2: Formulation and modeling of the finite-difference approximation of the proposed iterative schemes.

Step 3: Algorithms of the proposed iterative schemes.

Step 4: Numerical simulations.

Step 5: Evaluation and analysis.

End.

2.   Materials and methods

2.1.   Harmonic function

Consider Laplace's partial differential equation in two dimensions:

where $U\left({x, y} \right)$ is some unknown function of two variables. A harmonic function in the region $\Omega \subset {\Re ^n}$ should be able to satisfy the generalization of Eq (1), as below.

in which ${\nabla ^2}$ is the n-dimensional Laplace operator (or Laplacian as it is often called), the ${x_i}$ is the i-th coordinate in a Cartesian system and $n$ is the dimension of the region. This work only focuses on two-dimensional problems with $x - y$ notation represented as $i - j$. The region boundary of $\Omega$ comprises internal and external walls, boundaries of obstaclesand the goal position. According to ^[7], the existence of the spontaneous formation of local minima in the solution region cannot arise because harmonic functions obey the min-max principle. Thus, the only critical points that can happen are saddle points. The escape from such a critical point would then be found by searching the neighborhood. Any path deviation or disturbance from the critical point will eventually cause a smooth path throughout the region. For this reason, harmonic potential provides a highly beneficial option in navigation, as it offers a complete path-planning algorithm. The equation of Laplace can be solved with a range of numerical techniques, i.e., Jacobi, Gauss-Seidel and successive over-relaxation (SOR) are the standard approaches ^[11,12]. The development of the SOR iterative method has sparked interest in examining and resolving several issues. Since its establishment, the generalization of over-relaxation iterative approaches has been the subject of extensive research, and it has had promising outcomes. Therefore, this work will solve Eq (2) by exploiting the combination of skewed block accelerated iterative approach with the accelerated over-relaxation (AOR) and two-parameter over-relaxation (TOR) methods, for rapid computation.

In the configuration region for this model, the robot is described as a node point and the designated area is defined in the mesh grid pattern. As stated before, the off-line mode approach is used to compute the Laplacian operator of the point-robot designated area by utilizing the analogy of heat transfer. All wall boundaries and obstacles (with the highest potential value) stand for heat sources, while the goal position (with the lowest potential value) represents the heat sink. This heat transfer activity is modeled from Laplace's equation and then numerically solved to gain the heat distribution, which represents the harmonic potential for each nodal point in the mesh grid. By making use of the heat distribution property which flows from higher to lower temperatures, a gradient search can be used to generate the path from any starting position with a high potential value to the goal which has the lowest potential value. The path-planning algorithm utilizes the gradient descent search (GDS) to determine a feasible route for the robot to travel around the environment safely from the start point to reach the goal position. From the current point, the GDS method examines the potential values of its eight neighboring points on the finite-difference grid and simply picks the node with the lowest Laplacian potential value. This procedure is repeated until it reaches the target point ^[7,13,14,15].

This work aims to imitate the aforesaid paradigm for path planning, describing the solution of Laplace's Eq (2) over the analogy of temperature (for the potential) and heat flow (for the pathway). The experiment was carried out in a two-dimensional domain with walls and various forms of obstacles. A new technique called the skewed block two-parameter over-relaxation (SBTOR) iterative approach is proposed to solve Eq (2) and obtain the potential values for each node. For the purpose of performance comparison, the existing block successive over-relaxation (BSOR), block accelerated over-relaxation (BAOR), block two-parameter over-relaxation (BTOR), skewed block successive over-relaxation (SBSOR), and skewed block accelerated over-relaxation (SBAOR) techniques are also investigated.

2.2.   Skewed block technique

The pioneering work of the skewed block technique (originally known as explicit decoupled group, EDG) arose from Abdullah's ^[16] focus on solving two-dimensional Poisson equations. The EDG iterative scheme essentially engaged the half-sweep technique in a block frame on the rotated mesh grid. By virtue of the proficiency of this technique, it has been widely used to solve numerous differential equations ^{[17,18,19,20,21,22,23,24,25]}. In the robotics literature, the standard Gauss-Seidel ^[18,19] and SOR ^[3,20] schemes have been employed to solve Eq (2). This work, however, generalizes a reliable numerical solver, when addressing the solutions of Laplace's equation, by utilizing a skewed block technique.

Fundamentally, merely half of the node points in the mesh region were computed while using a skewed block approach. In contrast, the explicit group (EG) scheme is basically a full-sweep (FS) technique in a block frame, evaluated every node point in the mesh region. Figure 1 provides an example of each EG and EDG scheme in a $10 \times 10$ mesh region. To obtain the Laplacian potentials in the configuration region, only the black node points (see Figure 1) will be computed throughout the iteration process, at least until the convergence condition is satisfied. The convergence criteria in this work are denoted as $\varepsilon$, and the stopping condition is subject to $\left\| {{u^{(k + 1)}} - {u^{(k)}}} \right\| \leqslant \varepsilon$. The rest of the residual nodes, or the white node points (see Figure 1(b)), will be calculated independently using a direct technique ^[16,26,27]. Moreover, the block iterative techniques involve four Laplacian potentials simultaneously per calculation, thus indirectly speeding up the computation. As shown in Figure 1, a set of one-node point and two-node points can be earmarked to measure groups of nodal points that are adjacent to the boundary.

From Figure 1(a), each and every formulation that practices the EG iterative approaches shall compute a set of four nodes in a group all at once throughout the iteration process (excluding the unique set adjacent to the boundary). Meanwhile, the EDG is essentially derived from a skewed 5-point finite-difference approximation (5-FDA) ^[16]. As illustrated in Figure 1(b), the configuration domain for the EDG technique is distributed with two different types of nodal points, i.e., the black node ($\bullet$) and the white node ($\circ$). By pairing the identical node, the solutions of each group could be executed independently for each pair. On account of this independency, roughly half of the iterations across the solution domain are to be executed using either type of node, reducing performance time as well as computation complexity. Additionally, to compute the unique sets adjacent to the boundary, the direct method ^[27] is computed by Eq (1). One of EDG's main advantages is that these approaches diminish the computational complexity by measuring only half of all node points. The FS and half-sweep (HS) iterative approaches' computational stencils are shown in Figure 2, where $h$ is the distance between node points for each direction. These molecules can be displayed in matrix form for FS and HS cases, respectively, as

and

2.3.   5-point finite-difference approximation

Laplace's equation is a special case of Poisson's equation ${\nabla ^2}U = f$, in which the function $f$ is equal to zero. Hence, a two-dimensional Laplace equation (2D-Le) could be written using Eqs (1) and (2) altogether.

The system must be discretized using the finite-difference approach in order to solve the 2D-Le before it can be computed effectively using the numerical technique. The simplest second derivative of 5-FDA (also known as an FS iteration) is

in which $U$ is a function to satisfy Laplace's equation, while $u$ is the potential node at the point $\left({x, y} \right)$. In addition, the $x - y$ plane of the configuration region is rotated by ${45^ \circ }$ clockwise to provide the approximation that is based on the cross-orientation operator ^[28,29]. As a result, a skewed 5-point approximation, commonly referred to as the HS iteration, is produced and denoted by

Using the notation ${U_{ij}}$ to represent the solution of the potential node for 2D-Le through the mesh grid point $\left({{x_i}, {y_j}} \right)$, the Laplace equation for the FS case is discretized using a standard 5-point stencil, as follows:

whereas the discretization of a rotated case via the 5-point stencil is

A linear system is then generated by applying Eqs (3) and (4) to the Laplace's problem subject to the 2D-Le. Matrix notation can be used to build this linear system, which results in large and sparse matrices, and is represented as

in which $A$ represents known coefficients that are expressed in matrix form and $b$ is a known vector, while $u$ is the unknown vector. The expansion of each coefficient for both the FS and HS cases is discussed in depth in ^[25]. Since the linear system of Eq (5) creates large and sparse matrices, iterative approaches can be used to solve this problem using conventional point or block techniques ^[16]. In order to elucidate the block technique, the fundamental concept behind the point technique must first be described. Thus, from the second-order central finite difference of Eqs (3) and (4), the point Gauss-Seidel iterative schemes for FS and HS can be reformulated and designated respectively as below.

By integrating a weighted parameter ^[30], the point SOR iterative technique for FS and HS cases is given respectively as

Normally, the optimal weighted parameter $\omega$ is determined in the range of $1 \leqslant \omega < 2$ ^[31]. If the parameter value is equal to 1, the approach is essentially simplified to the conventional Gauss-Seidel method. The SOR and Gauss-Seidel methods are actually rather comparable, with the exception that SOR uses a scaling factor to reduce the approximation error.

The AOR contributes a supplementary relaxation parameter that, in this work is denoted as $r$, along with the weighted parameter $\omega$ from the SOR technique. Both parameters are used to generate iterative schemes that are able to accelerate the convergence rates. AOR is a simple yet powerful scheme, due to the fact that two parameters are involved, for the larger linear system first introduced by Hadjidimos ^[32]. The point AOR schemes for standard and skewed approaches are shown respectively as

According to Hadjidimos, the $r$ value is typically selected to be close to the corresponding $\omega$ value. Additionally, he claimed that the uncertain optimal values of these two parameters did not impose any limitations on the ability in to obtain the minimum iterations number.

The TOR iterative method is basically a generalization of the AOR method with another additional parameter, $r'$. The TOR scheme reduces to the AOR iterative scheme provided the value of $r'$ is equal to $r$ ^[33]. As stated in ^[33], the additional parameter $r'$ could provide faster convergence if an appropriate value is chosen. The point TOR iterative schemes for FS and HS are stated respectively as

Similar to the AOR approach, all of the parameters are in the range of $[1, 2)$ and they are chosen to be near to the value of the corresponding SOR parameters ^[32,34]. In accordance with the 2D-Le, all of the proposed iterative schemes simply substitute each node's value with the average of its four neighbors' values. In this work, the node values that represent the walls, obstacles and target position remain constant.

2.4.   Explicit group two-parameter over-relaxation scheme

As discussed earlier, the block iterative technique obtains four Laplacian potentials per computation. As illustrated in Figure 1(a), all formulations using EG iterative methods compute a group of four nodes at once during the iteration process. Basically, the EG iteration is based on a group of a small number of points, and it is derived using the standard 5- FDA. By analyzing Eq (3) and the point SOR iterative scheme in Eq (8), the block SOR scheme ^[16,28,31] can be written as

for

For this scheme, it is also possible to generate a linear system of the same form as Eq (5). The scheme can later be translated to the coefficient $A$ after determining the inverse of its matrix, as represented below.

with

Based on the point SOR concept using 5-FDA, the block SOR (BSOR) iterative scheme for Eq (15) is now denoted as

Besides, by examining Eq (3) and the point AOR approximation in Eq (10), the formulation of the block AOR method is stated as ^[35]

with

Again, Eq (17) may also be transformed into a linear system in the form of Eq (5) and translated as in Eq (15). Hence, the block AOR (BAOR) iterative scheme can also be expressed as Eq (16).

The general expression of the block TOR (BTOR) iterative scheme, considering Eq (3) and the point TOR iterative method in Eq (12), can be written as

where

2.5.   Explicit decoupled group two-parameter over-relaxation scheme

As noted in the previous section, the EDG is derived from skewed 5-FDA ^[16]. The solution domain for the EDG technique is shown in Figure 1(b). From the solution domain, the block EDGSOR formulation may be expressed in matrix form, as Eq (14), with

and

Both Eqs (19) and (20) are deduced from the system of linear equations in Eq (14). These equations can be solved independently, as shown below.

hence, it can be defined as follows (respectively):

The algorithm for SBSOR scheme can be enforced by using either Eq (23) or (24). Both equations lead to an equivalent solution.

Similarly, adopting the skewed AOR formula from Eqs (7) and (11) in the block scheme yields a new formula known as SBAOR ^[36]. As a result, the SBAOR iteration scheme is described as

where

The same goes for the rotated block TOR, which applies the skewed TOR formula from Eq (7), as well as Eq (13), to the block scheme, and it provides the SBTOR iterative scheme as follows:

in which

As stated before, there is no generic formula for defining the optimal value for all weighted parameters, which results in a minimum number of iterations. Consistent with ^[32], the values of $r$ and $r'$ are commonly appointed as close to the corresponding $\omega$ value. Note that not every parameter value leads to convergence. The optimal values for all three relaxation parameters for each FS and HS case are different. Due to each relaxation value being predetermined before execution, the complexity related to determining the relaxation values for the overall computation is unaffected. Table 1 provides a list of the optimal relaxation values used throughout the experiments. Sensitivity analysis has been used to specify each relaxation parameter's value. The values are the same for every N-size grid, although they differ from other grid sizes and are not significant enough to be revealed; they vary by roughly 0.0001.

Algorithm 1 describes the implementation of the SBTOR scheme in accordance with Eq (26) to solve the two-dimensional Laplace problem.

3.   Results and discussion

Following the aim of this study, an experiment has been conducted by using a robot 2D simulator ^[37] built by the authors to analyze the competency of the proposed algorithm, i.e., SBTOR scheme, as a tool to solve the path-planning problem. In the block mesh grid domain (see Figure 1) with the execution of the skewed block iterative method, approximately half of the group of node points are computed throughout the iteration process. This will logically result in a substantial reduction of the computational complexity by nearly 50% of the entire operation. The simulations are carried out on four different configuration spaces with five different grid sizes in order to assess the performance of the proposed algorithms. The aforementioned performances were analyzed with three criteria: (i) number of iterations, (ii) CPU time in seconds and (iii) success of path planning. The amount of CPU time and iteration count required by every algorithm is recorded once the tolerance rate for the iteration process is reached. The GDS approach then makes use of the computed Laplacian potentials to guarantee that a path from the starting point to the goal position is successfully constructed.

In the designated areas, several obstacles with different sizes and shapes are positioned. The Dirichlet boundary condition was utilized in the initial configuration, in which the walls and obstacles are set at the high potential value (the highest temperature). The start position had no initial value specified, whereas the goal location was set at the low potential value (viewed as the lowest temperature). The remainder of mesh grid points were initialized to zero. A machine running a 2.50-GHz processor with 8 GB of RAM has been used to carry out the simulation experiments. Numerical computation was performed iteratively until the stopping condition was met. Provided the Laplacian potential values give no further changes and the difference of current iterations (k) with the next iterations (k+1) is extremely small, e.g., ${1.0^{ - 15}}$, the loop should be stopped. This accuracy level was necessary to prevent any point on the surface which has critical point slopes that result in an abrupt absence of descending gradient patterns.

Each iterative scheme compared in the experiments is shown in Table 2, along with the iteration counts (k) and the execution time (t). Both the block (EG) and skewed block (EDG) schemes evidently show that the TOR iterative approach surpassed its generalization techniques. The computational complexity of the skewed block iterative approach will decrease by approximately half of that relative to the HS process. Since only half of the total nodes in the mesh grid are examined, the skewed schemes surpass the standard block schemes that utilize the FS process. The size of the mesh grid and pattern of the region plays a significant role in producing the outputs of k and t. Obviously, the larger grid size will require more iterations and a longer time to be executed. The improvement ratio between each over-relaxation technique for the smallest grid size (300) to the largest (1500) varies. In terms of the iteration count, on average, the BTOR efficiently decreased by 3-6% compared to BAOR, and 12-20% next to BSOR, whereas SBTOR reduced it by more than 3-8% compared to SBAOR, and roughly 19-32% of SBSOR. However, in terms of execution time, BTOR is between 3-7% lower relative to BAOR, and 14-18% in comparison with BSOR, while SBTOR effectively depreciated approximately 1-4% from SBAOR, and 15-23% in comparison with SBSOR. This study still contributes good results since previous studies ^[3,13,38] only did the experiments on extremely small grid meshes, i.e., 50 by 50 and 70 by 70 grid. Figures 3 and 4 show the results of Table 2 graphically for the number of iterations and CPU time, respectively. From the graphs, again, the skewed block TOR undoubtedly offers the optimum performance in contrast to its predecessor.

The ideology of path-planning flow within this experiment starts with establishing the initial start position and the target location. Next, we identify the optimal parameter according to the corresponding iterative schemes. Once the harmonic potential is generated from the selected algorithm, a complete smooth path is developed through the use of the GDS technique. This impression could describe the iterative scheme that tracks the descending gradient from its starting point to the next consecutive points with lower potentials from previous points, up to the target location (with the lowest potential value). These successful trails can be observed as demonstrated in Figure 5. In this simulation, the starting point is indicated by a green square point, while a red circle point is allocated for the target. As shown in Figure 5, each start location from every region has efficiently accomplished the route by arriving at the target position along while escaping any walls and obstacles in between (if any). The path trajectories can be really fast because they only involve the gradient evaluation of the precomputed Laplacian potential ^[13]. All four configuration regions in this experiment are relatively simple relative to those in ^[38].

There are various techniques to solve navigation problems. According to Patle et al. ^[1], these techniques can be divided into two main groups, i.e., classical approaches and reactive approaches. Due to the lack of artificial intelligence tools at the time, the classical approach was initially quite common as a tool to address robot navigational issues. When a task is performed utilizing a classical approach, it can be noticed that either a result will be acquired, or it will be confirmed that there is no result. This method is less suitable for real-time implementation due to its significant disadvantages, including its high processing cost and inability to respond to the uncertainty existing in the environment. Reactive techniques have surpassed traditional ways in recent years as the preferred method for mobile robot navigation. They are extremely adept at managing the environment's uncertainty. Meanwhile, Sanchez-Ibanez et al. ^[39] classified path planning into four categories with two major subcategories (see Figure 6). There are characteristics shared by two adjacent subcategories from distinct categories. The schematic also shows how some subcategories lean more toward global planning or local planning than others. From the reviewed papers ^[1,39], it is safe to say the proposed algorithm in this study has its own edge and drawback. Algorithm 1 has clearly implemented the improved version of potential field approaches. In essence, the goal and obstacles operate as charged surfaces, and the overall potential creates the imaginary forces on the robot. The robot is drawn toward the goal and kept clear of the obstacles by this imaginary force ^[8]. Later, the robot will travel along the negative gradient to avoid obstacles and arrive at its destination point. To prevent a local minimum problem, this work utilizes the harmonic function ^[7]. Moreover, Algorithm 1 exhibits superior computing execution when the skewed accelerated relaxation technique is implemented, as it performs much faster to obtain the solution of Laplace's equation to solve the path-planning problem.

4.   Conclusions

The iterative approaches investigated in this work are from over-relaxation families including the implementation of the block and skewed block schemes to enhance the execution performance as well as reduce the execution time. The introduction of accelerated weighted parameters to respective skewed nodes further improves the overall performance of the newly proposed SBTOR iterative schemes, yielding promising results. The experiments evidently demonstrate that the solution to the robot path-planning problem is feasibly solved by using numerical approaches owing to advanced algorithms and the accessibility of fast machines these days. From the table of results, the SBTOR scheme was found to significantly outperform its predecessor's techniques in terms of both the number of iterations and the time taken. The computational performance is unaffected by increasing the number of obstacles; in fact, it will complete more quickly because the calculation disregards the regions engaged by the obstacles. It can be said that the larger the area of obstacles taken, the less the required calculations and storage, or, in other words, the obstacles reduce computational complexity. Once again, the authors would like to highlight that the SBTOR outperformed the SBAOR (approximately by 3-8%) and SBSOR (by about 19-32%) in terms of the iteration count, while the SBTOR saved roughly 1-4% over SBAOR and 15-23% over SBSOR in terms of the processing time. The application of the skewed block TOR scheme families in robot pathfinding and in Algorithm 1 is another aspect of this study's originality or novelty. The authors intend to investigate the red and black strategy ^{[23,24,40,41,42]} by applying the proposed approaches in future work. It is believed that the aforementioned strategy will enhance and refine the overall computation.

Acknowledgments

This research was financially supported by Universiti Pertahanan Nasional Malaysia. The authors also acknowledge support from Science Foundation Ireland (SFI) under Grant Number SFI/16/RC/3918 (Confirm), and the Marie Skłodowska-Curie grant agreement No. 847577, co-funded by the European Regional Development Fund.

Conflict of interest

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