Zheng Sun , Zhu-Feng Shao , Hui Li

a School of Mechanical and Electrical Engineering, Zaozhuang University, Zaozhuang, 277160, China

b Department of Mechanical Engineering, School of Mechanical Engineering, Tsinghua University, Beijing,100084, China

Keywords:

ABSTRACT Path planning is a key technique of autonomous navigation for robots, and the velocity field is an important part. Constructing velocity field in a complex workspace is still challenging. In this paper, an inner normal guided segmentation algorithm in a complex polygon is proposed to decompose the complex workspace in this paper.The artificial potential field model based on probability theory is then used to calculate the potential field of the decomposed workspace, and the velocity field is obtained by utilizing the potential field of this workspace. Path optimization is implemented by curve evolution,during which the internal force generated in the smoothing process of the initial path by a mean filter and the external force is obtained from the gradient of the workspace potential field. The parameter selection principle is deduced by analyzing the influence of several parameters on the path length and smoothness.Simulation results show that the designed polygon decomposition algorithm can effectively segment complex workspace and that the path optimization algorithm can shorten and smoothen paths.

1. Introduction

Path planning, a key technique of autonomous navigation for the robots, aims to find an optimization method to reach the end point without collision under certain constraint conditions by fusing multi-sensor informations [1,2]. At present, path planning algorithms can be classified as: Surface based methods and Graph based methods[3].There exists following methods in Surface based methods: (1) Potential field method. This method is simple, and there are some problems such as falling into local minima and nonreachable targets exists[4-7].(2)Harmonic function.This method based on Laplace equation, the main drawback of this approach is that it is prone to quantization effects and it is slow[3].Graph based methods mainly includes: (1) Dijkstra algorithm. This method is based on grid maps and can generate the shortest path.However,a large amount of memory space is required for the solution process[8-11]. (2)A* method. A* is similar to Dijkstra algorithm in that it can be used to find the shortest path.It can use a heuristic to guide itself[12-14].(3)Visibility Graphs.This method is comprehensive and can easily find the shortest path, but has poor flexibility[15-17]. (5) Sampling-based method. This algorithm has poor flexibility and high resource consumption [18-21]. (6) Machine learning. Machine learning, especially deep learning, is a research hotspot in recent years [22]. Machine Learning in Robot Motion Planning has attracted the attention of many researchers [23].

Besides the methods described above, Eikonal equation based path planning algorithm, an effective method based on partial differential equations, has attracted increased attention [24,25].This method belongs to Surface based methods and has the following advantages: (1) There are many ways to solve Eikonal equation, some of methods are fast and efficient. (2) No need to consider the local minima problem.Eikonal equation based method operates in two steps:

(1) Eikonal equation is solved to generate the global distance field D（x，y）. The minimum distance is located in the goal point;

(2) The gradient distribution is obtained according to the distance field D（x，y）,and a continous path is obtained by using Runge-Kutta algorithm to move from the initial point to the goal one.

There are the following problems in the above two steps.

(1) The solution of Eikonal equation depends on the velocity field. The ideal velocity field satisfies that speed is 0 at obstacles and speed equals to 1 on free space. If obstacles are complex, the velocity field is difficult to construct.

(2) The distance field obtained by solving eikonal equation is discrete.Moreover,the cost value of the obstacle space is far higher than that of the free space, thereby causing gradient mutations near the obstacle space. When the path is set up by using the gradient descent method,path oscillation easily occurs near obstacles.

(3) The path generated is not smooth.

For solving the problem (1), velocity fields are constructed by using the cost surface generation method proposed in Ref. [26] in this paper.This method is inspired from Logistic distribution and is applicable only to convex polygon obstacles. So a polygon segmentation algorithm is proposed to extend the application scope of the method. To solve the problems (2) and (3), an improved path optimization method based on average filtering and curve evolution is developed to optimise the path generated by Eikonal equation.

The remainder of this paper is organized as follows. Section 2 introduces the polygon decomposition method we put forward.Section 3 presents the path optimization algorithm. Section 4 shows the simulation results and discusses the performances of the path-planning algorithm. Finally, conclusions are drawn in Section 5.

2. Polygon decomposition method

2.1. Notation and definitions

As shown in Fig. 1, the obstacle can be abstracted area is described by a polygon area, denoted by O（V，E，dir）. V represents the set of vertices of the obstacle area,V =｛v1，v2，v3，…，vm｝,viis the vertex of the obstacle area,and m is the number of vertices.E is the set of edges. eiis the set of vertices belonging to the ithedge,dir is used to set the order number in the specified direction dir =clockwiseorcounterclockwise, or dir =counterclockwise in this paper.Eiis the vector representation of the ithedge.Niis the inner normal of the ithedge that points to the inside of obstacle.

Fig.1. An example of a obstacle.

Fig. 2. The model of obstacle’s edge.

Fig.2 shows an edge model of the obstacle polygon.Line vivi+1is determined by viand vi+1,and edge vibelongs to line vivi+1.θ is the angle between line vi+1and the x-axis. ρ is the perpendicular distance from the origin to line vivi+1.Thus,line vivi+1can be written as

Theorem 1. Let

means the ith edge. For any point P(x，y) in the workspace, F（x0，y0）is a point on line vivi+1,and FP⊥vivi+1.If FP-→and Niare in the same

Proof. As shown in Fig. 2, the coordinate of point D is (- ρsinθ，ρcosθ) according to (1). The parametric equation of any point on line vivi+1is given as follow:

where s is the length from point D to x0along line vivi+1. The distance of point D to line vivi+1is

Substituting(2)into(3),the following formula can be obtained:

Then, calculating the derivative of r2with respect to s

Definition 2.viand vjare two vertices of obstacle. If vivj-then vjis the visible point of vi.

2.2. Determine split point

The potential field model proposed in Ref. [22] is only suitable for the case where the obstacle space is a convex polygon.As such,developing a convex polygon decomposition method for the obstacle space confined by a concave polygon is necessary. The internal normal line guided complex concave polygon decomposition algorithm is presented in this section. The detailed steps for the algorithm are described as follows.

(1) Determine convexity-concavity of the vertex in the obstacle area.The cross product method is a common method used to distinguish convex and concave polygons.Vertex viwith the coordinate ci=［cxi，cyi］ must be determined. Edge eiand ei-1are assumed to meet the condition ei∩ei-1= vi. Thus, the edge vector of eiand ei-1can be expressed as Ei-1= ci-1-ci，Ei=ci-ci+1,and the convexity-concavity of vertex vican be determined by the following method [[23]]:

(2) Determine candidate visible points of the concave point.

In Ref.[27],an improved Lee algorithm was proposed to search for visible points,as illustrated in Fig.3.Here,the polygon is divided into four regions numbered I, II, III and IV respectively, and the visible points of concave viare determined in I,II and III regions(see Fig. 4).

As shown in Fig. 3, two lines along the two adjacent edges of a selected concave point viare drawn, and the obstacle region is divided into four areas. Research shows that when split point in region I, segmentation effect is the best [28,29]. The optimal segmentation can be obtained when the visible points appear in region I. In this paper, the inner normal guided method is used to determine the candidate visible points at region I.

(a) The inner normal line of all edges is calculated. The inner normal direction can be considered the edge vector rotated by 90°in the counterclockwise direction and calculated as

Fig. 3. Searching the candidate visible point.

(b) For any vertex vj（vi≠vj）, Fi-1（xj，yj） and Fi（xj，yj） are calculated. If the following conditions are satisfied, vertex vjwill be located in region I, and the point can be added to the candidate visible point set.

(3) The visible points are located among the candidate visible points. The set PCV（vi） is used to represent the candidate visible points, and the visible point is defined as follows:Assume that there exists a point vvi∈PCV（vi）. If all points on line segment vivvi(except point viand vvi) belong to Int（O）,then point vviis the visible point of vi[28,29].Theorem 2 can be used to determine whether a point is visible.

Theorem 2. viis a concave point andIf any edgesatisfiesthen

Proof. (By contradiction) Suppose that there is a pointnd S?Int（O）.Therefore there exist an edge vzvz+1satisfies sign（F（xs，ys））≠sign（Nxz） or sign（F（xs，ys））≠sign（）.

(b) As shown in Fig. 3(b), if line segment vivviis located on one side of line segment vzvz+1, the value of Fz（vivvi） is always positive or negative, and sign（Fz（x，y）） is constant.

(c) As shown in Fig. 3(c), line segment vivviintersecting line segment vzvz+1is in contradiction with known conditions obviously.

The above formulas show that, once kivand kzcan be determined,is not related to x and y. Thus, sign（Fz（x，y）） is monotonous in the direction ofIf there is a point C between points O and S satisfying F（xC，yC） =0,then point C on edge vzvz+1can be obtained, which is contrary to the hypothesis.

Fig. 4. Relative position of vivvi and vzvz+1.

Thus, Theorem 2 holds.?

(4) The segmentation point is determined. This step is divided into the following two situations:

(a) The number of visible points is 0.In Fig.5,a ray viA along theand an inverted ray viB along the edgeare drawn.Then,the optimal segmentation point is in the angle bisector of ∠AviB,which is also adopted in this paper[30,31].

The angular bisector viS can be viewed as the counterclockwiseof viA. Thus, the vector of the angular bisector is

Since the number of visible points is 0,the intersection between the angle bisector and each side must be obtained to determine the segmentation point.As shown in Fig.6,viS is the angle bisector and vjand vj+1are the adjacent vertices of the polygon.Theorem 3 can be used to judge the intersection of ray viS and line segment c.

Theorem 3.When the following inequalities hold

Fig. 5. The case without visible points.

that is,vjand vj+1are on the two sides of viS,and the angle betweens less than π, and ray viS intersects with the line segment vjvj+1.

Fig. 6. The intersection of a ray and a line segment.

Thus,Theorem 3 holds.

Theorem 4.If concave point vihas no visible point, there are m intersections between bisector viSand the polygon, Ijis the jth intersection, then if intersection satisfies I*=argminIis the segmentation point.

Proof.By contradiction. Assume that there are multiple intersections between bisector and the polygon,the intersection Sois the segmentation point and the length|viSo|of line segment viSois not the shortest.Because all intersections collinear with viand the length of ‖are not the shortest,there is another intersection S1on line segment viSoandis the shortest,according to Theorem2,Sois not an visible point of vi.Additionally,since split points must be visible, this situation is contradicted with the initial assumptions.

Thus,Theorem 4 holds.

As shown in Fig. 7, there may be multiple edges in the polygon intersecting with the bisector, and there exist multiple intersections. Thus, Theorem 4 can be used to determine the segmentation points.

(b) The number of visible points is not 0.Referred to Fig.7,the jth visible point is represented as,and the segmentation point is selected from the visible points. In this case, the concave point with the minimum angle betweenand bisector viS is preferred,and,if no concave point exists,the convex point with the minimum angle betweenand bisector viS is selected.

(5) The polygon is divided into two polygons after determining the segmentation point and repeating steps(1)-(4)until all partitions are convex polygons.

Fig. 7. The segmentation point determination without visible points.

3. Path optimization algorithm

A new path planning algorithm is proposed based on the potential field model in this section.The algorithm is mainly divided into two parts:(1)calculation of the initial path through the eikonal equation; (2) optimization of the initial path by mean filter and curve evolution.

3.1. Initial path calculation

The global cost function with regard to the destination is computed by the following eikonal equation.

where, f（x，y） is the velocity field. The distribution of velocity field should conform to the distribution of obstacles in the workspace.In Ref. [26], a potential field based on the logistic distribution is proposed and we use this method to build the velocity field.

where:

U is the potential field of the workspace,nobis the number of obstacles in the workspace,is the number of polygons decomposed from the kth obstacle,is the number of sub-boundary pairs in the mth polygon in thekth obstacle,j，1andare the potentials of the jth sub-boundary pair in the mth polygon in the kth obstacle,G is the cumulative distribution function of the logistic distribution,i.e., and σ is the standard deviation.

The eikonal equation is solved and the distance from each point in the free space to the obstacle is determined to obtain the initial path using the fourth-order Runge-Kutta method based on the gradient of D [32].

3.2. Path optimization algorithm

Fig. 8. The Eikonal equation based path plan under σ=0.1 and h = 1.

In this section, the performance of the initial path obtained by the eikonal equation introduced in Section 3.1 will be analyzed.Figs. 8 and 9 respectively show the initial paths obtained under variances σ = 0.1, 0.0001 and the step h=1 of the Runge-Kutta method. The figures clearly show that the above mentioned initial path generation algorithm methods have the following two disadvantages:

(1) The initial path may pass the obstacle region at the corners;

(2) The initial path is not adequately smooth enough and its length may be decreased further.

Thus, a new path refinement method inspired by mean filterbased initial path evolution and the potential field negative gradient algorithm is developed to optimise the initial path.

The mean filter method is a common curve smoothing algorithm.The initial path can be regarded as a curve made up of many points.The position of the ithpoint at the lthround can be written as:

and the increment of the ithparticle position is

This relation is similar to the curve evolution of a heat equation.The initial path will evolve into a straight line, andcan be called the internal force of the path curve in the lthround.

The external force must be designed to counteract the influence of the internal force of the path curve,move the path point outside of the obstacle and prevent the path point of the free zone from moving inside the obstacle. In Ref. [26], the author point that the max gradient amplitude appears in the corner boundary of the obstacle. Thus, the gradient direction runs from outside of the obstacle to inside of it, and the negative gradient direction is reversed. Thus, the negative gradient can be used to design the external force.

Considering that the external force is exerted by the obstacle region and the force can be constructed by the workspace model,the curve will balloon outward. The derivative of the ithparticle position is

Fig. 9. The Eikonal equation based path plan under σ = 0.0001and h = 1.

and the derivative of the position of the ithpoint is

The path optimization iterative formula can be obtained by using the finite difference method as follows

3.3. Convergence analysis of the path optimization algorithm

When all path points satisfy the conditionthen the positions of these points do not change.If the positions of all points no longer change, the path will not intersect with obstacles and become shorter and smoother than the initial path. The potential field of the workspace shows that the external force of a point on the path decreases as the distance from this point to the obstacle boundary increases. As such changes in the external and internal forces must be reasonably designed to ensure the convergence of the optimization algorithm.

The magnitude of the external force will be analyzed.According to the analysis in Ref.[26]:(1)if σ is large,the effective influencing distance of the external force on the obstacle boundary increases.In addition, if the distance between any two sub-boundaries decreases, these sub-boundaries will influence each other and cause the external force of every boundary edge be less than the theoretical maximum(2) if σ is small, the effective influencing distance of the external force on the boundary is small, the distance between the sub-boundaries has a slight effect on the external force of the boundary, and the external force on the boundary can be considered asThe distribution of internal force over the whole workspace is given as follows

where:

?U is the external force from the potential field,nobis the number of obstacles in the workspace,is the number of polygons decomposed from the kth obstacle,is the number of sub-boundary pairs in the mth polygon in the kth obstacle andandare the jthsub-boundary pairs in the mthpolygon in the kthobstacle.

From Eq. (37), the maximum value of the external force of the workspace appears at the corner of the obstacle. The distribution width of the external force is affected by the variance σ.That is,the σ is, the larger the distribution width but the smaller the peak.When 0

Fig.10. The path obtained under different h with σ=0.1 and Δt = 0.5.

Fig.11. The path length change under different h with σ=0.1 and Δt = 0.5.

as long as the internal force satisfies the condition

The internal force is related to the step size of the Runge-Kutta method used to calculate the initial path.Therefore,as long as the internal force satisfies the following formula,the initial path stops at the position where the internal and external forces are equal.

Fig.12. The path obtained under different h with σ=0.01 and Δt = 0.5.

Fig. 10 shows the optimized path obtained under different h with σ=0.1 and Δt = 0.5, Fig. 12 shows the optimized path obtained under different h with σ=0.01 and Δt = 0.5, and Figs.11 and 13 show the path length change.

Fig.13. The path length change under different h with σ=0.01 and Δt = 0.5.

The path optimization algorithm based on the combination of the mean filter and potential field negative gradient has the following problem: the path converges at some distance from the obstacle such that the path length is longer than the initial path length,as shown in Figs.10(a),(b),11,12(a)and 13,because that the internal force is much less than the external force.

Fig.14. Optimization path under σ = 0.1, h=1 and Δt = 1.

Fig.15. Path length change during optimization under σ = 0.1, h= 1 and Δt = 1.

Fig.16. Optimization path under σ = 0.01, h=1 and Δt = 1.

Fig.18. Optimization path under σ = 0.005, h=1 and Δt = 1.

The following two methods can be adopted to solve this problem: (1) Reduce σ. Reducing σ can decrease the effective influencing range of the external force, as shown in Figs. 10(a)-(c),12(a)-(c). If σ is too small, the effective influence range of the external from σ will not be obvious.Thus,the influence of the path convergence position due to σ also becomes less obvious. (2) Increase h. Increasing h can increase the internal force and particle density, thereby causing the path to intersect with the obstacle at the corner, as shown in Figs.10(c), (d),12(b)-(d).

Fig.19. Path length change during optimization under σ = 0.005, h=1 and Δt = 1.

According to Eq. (36), the convergence speed can be improved by increasing Δt.Figs.14-19 show the final optimization path and path length change under σ=0.1 or 0.01 or 0.005, and h=1 and Δt = 1.

The above simulation results reveal that the path point near the corner of the obstacle will oscillate when Δt is large,and oscillation increases with decreasing σ because if Δt is too large,the point will move in a wide range each round. If the corner external force increases, σ will drastically decrease. Thus, particles could suddenly pass from the small external force area into the large external force area, resulting in oscillation. The previous analysis that using Eq.(36) for path optimization will require several cycles of experimentation to find suitable parameters, which reduces the practicability of the algorithm. Therefore, the following improved optimization algorithm based on the characteristics of the workspace model is proposed.

Eq. (41) utilizes the steepness of the transition region of the potential field model and adds the parameter 1-U（ci） to the first item, thereby reducing the influence of the internal force of the path point on the obstacle region. Thus, a point in the path can jump out of the obstacle area as soon as possible, and the path length can be reduced.

4. Experiment results and discussion

In this section,the results of extensive simulation will be studied to validate the effectiveness of the proposed algorithm, and the influence of the related parameters on the algorithm performance will be analyzed,therefore.Ultimately,optimal parameters will be obtained.

4.1. Validity of the proposed inner normal guided polygon segmentation algorithm

The inner normal guided polygon segmentation algorithm is an important component of the path planning algorithm that directly related to the success of path planning. Thus, the feasibility and reliability of the algorithm must be verified. Fig. 20 illustrates the test polygon taken as an example to verify the feasibility and reliability of the proposed polygon segmentation algorithm.

For example,vertex 5 is a concave point,and the possible areas of visible points at vertex 5 are shown in Fig. 20. The candidate visible points of vertex 5 are clearly vertices 16 and 17.In this paper,(26) and (27) are used to judge the candidate visible points, and Table 1 shows the results of the case where of each vertex satisfies(26) and (27), here 1 means satisfied, 0 means unsatisfied.

Fig. 20. An example of polygon segmentation.

Table 1 that vertices 17 and 16 completely satisfy(26)and(27).The proposed inner normal guided polygon segmentation method can correctly judge the candidate visible points. Fig. 21(a) clearly shows that vertex 5 has no visible points.Theorem 2 can be used to judge visible points among the candidate visible points.According to the line segment intersection judgment method, there exists at least an edge e17e18intersecting with the line segment e5e16and an edge e1e21intersecting with the line segment e5e17.Thus,vertices 16 and 17 are not visible points of vertex 5. Because there is no visible point for vertex 5,it needs to determine segmentation point according to Theorem 3. Thus, determining edges that intersect with the angular bisector is necessary. Theorem 3 and Fig. 21(a)intuitively provide, the additional segmentation point on edge e21e1. Table 2 lists the cause where several edges satisfy Eqs. (29)and (30). Edges e17e18, e5e16, and e1e21satisfy Theorem 3, so the additional segmentation point must appear on these edges and,according to Theorem 4, the new segmentation point is on line e1e21with the coordinates (8.12,16) as shown in Table 3.

Fig.21(b)shows the second segmentation cycle,where vertex 3 is the selected concave point and vertices 10-13 and 15 are visible candidate points but only vertex 15 is visible.Table 4 shows that the vertices satisfy (26) and (27). Because only vertex 15 satisfies Theorem 2,this vertex is the unique segmentation point.The final results are shown in Fig. 22.

4.2. Influence of various parameters on the algorithm performance

As the proposed algorithm mainly involves three parameters,namely, σ, h and t, the selection of parameters will influence the performance of the optimization algorithm. In this section, path length changes will be used to reflect the influence of different parameters on the algorithm convergence performance, the value selection of main parameters will be analyzed,and the appropriate parameter values will be obtained.The simulation maps in Figs.23 and 24,are divided into three scenes,and the start and end pointsof each scene are shown in Table 5.

Table 1The result of and satisfies (7) and(8).

Fig. 21. The process of Polygon decomposition.

Table 2Partial edges satisfy Theorem 2 in the first segmentation round.

Table 3Candidate points satisfy Theorem 3 in the first segmentation round.

Table 4The results of and satisfies Eqs. (7) and (8) in the second round.

4.2.1. Influence of σ on the algorithm performance

In this section, the influence of variance σ on the convergence performance of the optimization algorithm is investigated. Fig. 25 shows the initial path obtained with different σ under h= 2 and Δt =0.8 in scene 1. When σ>0.2, path planning fails and the remaining path intersects with one of the obstacles at the corner.Fig. 26 shows the final path obtained after 15,000 cycles of optimization.

Fig. 27 shows the change in final path length obtained under different σ in scene 1. The following method is adopted to investigate the smoothness of the path. Firstly, the angle sequence is obtained by solving the angle changes along the path.Secondly,FFT transformation is performed on the angle sequence to obtain the frequency-amplitude graph. Figs. 28 and 29 respectively show the frequency-amplitude maps of angle changes in the initial and optimized paths.

Fig. 23. The simulation map.

Fig. 24. The potential field distribution of the simulation map.

Fig. 22. The final results.

Table 5The starting point and the ending point of the simulation scene.

Fig. 25. The initial path obtained under h=2 and Δt =0.8 in scene 1.

Fig. 26. The optimized path obtained under h=2 and Δt =0.8 in scene 1.

Fig. 27. The convergence of path length with different σ in scene 1.

Fig. 30 shows the initial path obtained with different σ under h=2 and Δt =0.8 in scene 2. When the path planning fails, the remaining paths intersect with one of the obstacles at the corner.The final path obtained after 15,000 cycles of optimization is shown in Fig. 31.

Fig. 28. The frequency-amplitude map of the angle change in initial path in scene 1.

Fig. 29. The frequency-amplitude map of the angle change in final path in scene 1.

Fig. 30. The initial path obtained under h=2 and Δt =0.8 in scene 2.

Fig. 32 shows the change in final path length obtained under different σ in scene 2.Figs.33 and 34 are the frequency-amplitude maps of angle changes in the initial and optimized paths.

Fig. 35 shows the initial path obtained with different σ under h=2 and Δt =0.8 in scene 3.When σ>0.5,the path planning fails,the remaining path intersects with one of the obstacles at the corner. Fig. 36 is the final path obtained after 15,000 cycles of optimization.

Fig. 37 shows the change in final path length obtained under different σ in scene 3. Figs. 38 and 39 show the frequencyamplitude maps of angle changes in the initial and optimized paths.

Fig. 31. The optimized path obtained under h=2 and Δt =0.8 in scene 2.

Fig. 32. The convergence of path length with different σ in scene 2.

Fig.33. The frequency-amplitude map of the angle change in the initial path in scene 2.

Figs.25,30 and 35 demonstrate that a larger value of σ results in easier failure of the initial path planning, because the transition area of the workspace potential field model increases with increasing of σ and the potential field of the workspace is reduced if the width of the barrier is small enough. These conditions lead to no significant difference between the cost values of the obstacle and free areas.Thus,the obstacle area is wrongly regarded as a free zone,and the failure path is obtained when the ekional equation is used to solve the initial path.

Fig. 34. The frequency-amplitude map of the angle change in final path in scene 2.

Fig. 35. The initial path obtained under h = 2and Δt =0.8 in scene 3.

Fig. 36. The optimized path obtained under h=2 and Δt =0.8 in scene 3.

Figs. 26, 31 and 36 reveal that the optimization algorithm overcomes the disadvantage of the intersection of the initial path and obstacle boundary. Similarly, Figs. 27, 32 and 37 illustrate the final path length becomes shorter with decreasing of σ,because the transition area between the obstacles and free areas in the workspace model become narrower and the effective range of the external force becomes closer to the obstacle.As such,the internal force range increases and the path length decreases.

Figs. 28, 33 and 38 illustrate that the signal frequency distribution in the angle change sequence of the initial path is not uniform and that the high-frequency signal is strong.Furthermore,the amplitude of the angle change sequence increases with increasing of σ,which indicates that the path angle change is large.Figs.29,34 and 39 show the frequency-amplitude maps of angle change sequence in the optimized path.The amplitude of each frequency is visibly reduced after optimization.The signal reduction of the high frequency part is more obvious with increasing σ. Therefore, the smoothness of the path increases with increasing of σ, and the optimization algorithm can clearly improve the smoothness of this path.

Fig. 37. The effect of σ on the convergence of path length in scene 3.

Fig.38. The frequency-amplitude map of the angle change in the initial path in scene 3.

Fig.39. The frequency-amplitude map of the angle change in the final path in scene 3.

4.2.2. Effect of h on the convergence performance of the algorithm

In this section, scene 1 is taken as an example to analyze the influence of the initial path step h on the convergence performance of the algorithm. Figs. 40 and 41 respectively show the initial and optimized paths obtained with different h under σ = 0.1. Figs. 42 and 43 then respectively show the initial and optimized paths obtained with different h under σ = 0.01.

Fig. 40. The initial path obtained with different h under σ = 0.1.

Fig. 41. The optimized path obtained with different h under σ = 0.1.

Fig. 42. The initial path obtained with different h under σ = 0.01.

Figs.40-43 demonstrate that:

(1) when σ is fixed, h has little influence on the initial path;

(2) when σ is fixed, the higher the h, the greater the failure probability of path optimization;

(3) the larger the σ, the smaller the h, which causes failure of path optimization.In this article, the causes of these phenomena are as follows.

Fig. 43. The optimized path obtained with different h under σ = 0.01.

Fig. 44. The distribution of the path particles under bigger h.

(1) The initial path is calculated by the potential field. The potential field is determined when σ is fixed. Therefore, the initial paths nearly overlap with different h;

(2) Increasing h leads to the internal force. Thus, the path approaches the obstacles. The increase in h also causes a decrease in the node density of the path,as shown in Fig.44.The external force in the corner area of the obstacle is symmetrically distributed at the center of the angular bisector.Thus, the particles are very likely to be pushed toward the two sides of the bisector, which results in the situation shown in Fig. 44.

(3) Failure of path optimization due to the point in the path being too close to the obstacle and distributed at both sides of the bisector.The smaller σ is,the closer the path is to the obstacle. Thus, with increasing of σ, the h causing failure of path optimization becomes smaller.

Fig. 45 shows the effect of h on the path length, and Fig. 46 shows the effect of h on the path smoothness. Increasing h helps speed up convergence and reduces the path length. When h increases to a certain critical value, the path particles enter an oscillating state,as shown in Fig.45(b).As h continues to increase,oscillation increases, eventually leading to the failure of path optimization.

Fig. 46 demonstrates that h has little influence on the smoothness of the path. However, the smoothness of the path will be enhanced when h reaches a critical value.

4.2.3. Effect of Δt on the convergence performance of the algorithm

In this section, scene 1 is taken as an example to analyze the influence of Δt on the algorithm performance. The simulation parameters selected as σ=0.1， 0.01 and h = 4. Fig. 47 shows the effect of Δt on the path length,and Fig.48 shows the effect of Δt on the path smoothness. Fig. 47 reveals that Δt does not have much influence on reducing the path length but remarkably influences the convergence speed of the algorithm.

As shown in Fig.47(b),when h is large,the oscillation problems caused by large h can be overcome by reducing Δt, but the convergence speed of the algorithm may decrease. Fig. 48 also shows that when Δt increases, the smoothness of the path decreases. Indeed, this effect becomes increasingly obvious with decreasing σ.

5. Conclusions

In this paper,two problems of the Eikonal equation based path planning algorithm, establishment of velocity field in a complex workspace for the eikonal function and path optimization for the initial path obtained by the eikonal function,have been discussed.The main work is as follows:

(1) Using the polygon model of an obstacle and introducing the inner normal line of an edge, an improved polygon segmentation algorithm is deduced. The workspace can be decomposed into several convex polygons to satisfy the conditions required to use the artificial potential field model based on probability theory. A velocity field based on the potential field model is proposed to improve the accuracy of the generated path.

Fig. 45. The effect of h on the path length.

Fig. 46. The effect of different h on path smoothness.

Fig. 47. The effect of Δt on the path length.

Fig. 48. The effect of Δt on the path smoothness.

(2) A path optimization algorithm based on curve evolution is proposed. In this algorithm, median filter smoothing of the path generates internal force, and a potential field provides an external force. Thus, the path gradually becomes smoother and shorter under these combined forces.

The main conclusions are as follows:

(1) The path length become longer and smoother with increasing of σ. However, σ cannot be infinite, and excessively large σ can easily cause the initial path to pass through the narrow portion of the obstacle, leading to the failure on the initial path planning.

(2) When σ is fixed,h is less than the lower limit of critical range,the path length is reduced,and the path smoothness remains nearly unchanged with increasing h.When h increases to the critical range, path oscillation occurs, and the path smoothness decreases.When h is greater than the upper limit of the critical range, path oscillation disappears and the path smoothness is reduced, Thus, the path will pass through obstacles, causing failure of path optimization. The critical region of h increases with increasing σ.

(3) Δt does not affect the lower limit of the path convergence but affects the convergence speed. Reducing Δt can help overcome the problem of path oscillation caused by excessively large h.

Our future work is to improve the resolution locally if the path crosses into the obstacle.

Acknowledgments

This work is partially supported by the financial support of the ship segmentation intelligent manufacturing equipment solution and key common technology research, High-tech Ship Research Project of the Chinese Ministry of Science and Technology and the project of Shandong Provincial Key R & D Program (No.2019GGX104035).