# Three-Dimensional Path Planning Method for Autonomous Underwater Vehicle Based on Modified Firefly Algorithm

Authors: Chang Liu; Yuxin Zhao; Feng Gao; Liqiang Liu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Path planning is a classic optimization problem which can be solved by many optimization algorithms. The complexity of three-dimensional (3D) path planning for autonomous underwater vehicles (AUVs) requires the optimization algorithm to have a quick convergence speed. This work provides a new 3D path planning method for AUV using a modified firefly algorithm. In order to solve the problem of slow convergence of the basic firefly algorithm, an improved method was proposed. In the modified firefly algorithm, the parameters of the algorithm and the random movement steps can be adjusted according to the operating process. At the same time, an autonomous flight strategy is introduced to avoid instances of invalid flight. An excluding operator was used to improve the effect of obstacle avoidance, and a contracting operator was used to enhance the convergence speed and the smoothness of the path. The performance of the modified firefly algorithm and the effectiveness of the 3D path planning method were proved through a varied set of experiments.


---

## Body

## 1. Introduction

AUVs are in strong demand within both military and civil fields. More and more research institutions are carrying out research into AUVs. Path planning is one of the key technologies in autonomous decision-making for AUVs, and it has become one of the most talked about issues in the AUV field. Many methods for AUV path planning have been put forward and most of the AUV operating environments were assumed to be two-dimensional. For example, some researchers proposed a path planning method using artificial potential field methods [1, 2]. Chen et al. brought forward a global path planning method for AUV based on the sparse A* search algorithm [3]. Some novel optimization algorithms, such as the genetic algorithm (GA) and ant colony optimization (ACO), were also used to solve the problems of AUV path optimization [4, 5]. Although these methods were proved effective in solving path planning problems, they inevitably face some problems in practical applications because the AUV working environment is a 3D marine space.3D path planning is a classical optimization problem. With the development of optimization algorithms, many researchers apply new optimization algorithms to solve the path planning problem in recent years. Swarm intelligence (SI) belongs to an artificial intelligence (AI) discipline that became increasingly popular over the last decade, and it refers to a research field that is concerned with a collective behavior within self-organized and decentralized systems. Nowadays, swarm intelligence has formed one of the hottest topics in the developments of new algorithms inspired by nature. And methods of swarm intelligence are used to solve traditional NP problems and their excellent performance and great potential have been proved in many practical applications [6]. In recent years, SI-based algorithms are also applied for solving 3D path planning problem. For example, a 3D path planning method based on the ACO was proposed by Liu et al. [7, 8]. Particle swarm optimization (PSO) was also used to solve the problem of 3D path optimization [9–11]. These works prove that it is feasible to solve 3D path planning problems using SI-based algorithm.Firefly algorithm (FA) is a more promising SI-based optimization algorithm based on firefly social features put forward by Dr. Yang in 2008. Although this algorithm is similar to other bioinspired algorithms, it is relatively simple both in theory and in implementation. One of its most prominent advantages is that global communications can run based on individual exercise at the same time as random motion [12]. The existing research results showed that the algorithm was very effective in dealing with a lot of optimization problems, and it had gradually been applied to some kinds of optimization fields with good results [13–21]. As a new optimization algorithm, the basic FA has some drawbacks in convergence speed and stability. Some improvements have been put forward recently. For example, Dr. Yang added Levy flight to the random part of the FA formulation and Levy-flight firefly algorithm (LFA) was constructed. Experimental result showed that the improved algorithm was more effective in searching for global best values [22, 23]. Farahani et al. made some improvements to the random motion part of FA to ensure the random step is larger in the initial period, in order to avoid falling into local optima. The random step was reduced after several iterations to make the algorithm converge to the global minimum rapidly. Meanwhile, the movement style of FA was improved and directed movement was proposed. Experimental results showed that the improved FA had better performance than the basic FA [24]. dos Santos Coelho and de Andrade Bernert used chaotic sequences to adjust the parameters of FA, and the improved algorithm had better results in optimizing benchmark test functions of reliability and redundancy allocation [25]. Liu et al. made some improvements to ensure parameters of FA adjusted according to the iteration numbers. Compared to the basic FA, the performance and convergence speed were improved [26]. Farahani et al. put forward a new FA that increases convergence speed using Gaussian distribution to move all fireflies to global best in each iteration. The modified algorithm was tested on five standard functions. And experimental results showed that the modified FA had better performance than the basic firefly algorithm [27]. Fister et al. proposed a comprehensive review of the FA. The status of applications within various application areas was summarized. At the same time, the review showed that FA was simple, flexible, and versatile, which was very efficient in solving a wide range of diverse real-world problems [20].In this paper, a modified FA and a 3D path planning method based on the new algorithm are proposed. The length of random steps and parameters of the algorithm are designed to be adjusted according to the distance between two fireflies and the iteration times. The autonomous flight strategy has been introduced to avoid invalid flight. An excluding operator and a contracting operator have been used to improve the effect of obstacle avoidance, the speed of convergence, and the smoothness of the path. The paper is organized as follows. Section2 analyzes the performance of basic FA. Section 3 develops a modified FA. A 3D path planning method is proposed in Section 4. Experiments in Section 5 prove the performance of the modified FA and the path planning method. Finally, Section 6 contains the main conclusions and describes some problems that warrant further research.


## 2. Basic Firefly Algorithm and Its Performance Analysis

The FA is derived from the simplification and simulation of firefly group behavior. It has three idealized constraints which are derived from firefly features. (1) The firefly is not gender-specific. And they will fly to more attractive and larger brightness companion regardless of its gender. (2) Firefly attractive size is proportional to its brightness. And its brightness decreases with the distance between individuals. If there is no brighter or more attractive one, then it will flight randomly. (3) The brightness or attractiveness of a firefly is determined by the specified value of the objective function [12]. For the maximization problem, the luminous intensity is proportional to the value of the objective function. Based on the above three rules, the execution process of firefly algorithm is shown in Figure 1.Figure 1
The execution process of FA.The core of FA is that the absolute brightness of fireflies represents the value of objective function and the position of the fireflies represents the solution of the problem to be solved. The relative brightness is gained by comparing two fireflies and it is associated with attraction. One firefly is attracted by a brighter one and adjusts its position according to this attraction. Only the brightest firefly moves randomly. The formula of FA is as follows:(1)x→it+1=x→it+βijrijx→it-x→jt+αε→i,(2)βijrij=β0e-γrij2,
where β0 stands for the biggest attraction and γ is the absorption coefficient which controls the change in light intensity and determines the speed of convergence. α is the random coefficient which controls the range of movement. γij is the artesian distance between firefly i and firefly j·ε→i is the random vector decided by Gauss distribution [12].As described in formula1, the changing position of a firefly depends on its attraction to brighter fireflies and their random movement. At the initial period of algorithm running, fireflies are randomly distributed within the search area. The distance between each of them is greater, while βij(rij) is small and the range of movement is small, so the exploring ability of FA is insufficient with slower convergence. At the later stage of the algorithm running, fireflies gather around the optimal insect, the distance between each other is smaller, and βij(rij) is bigger. At the same time, random motion is still taking place which is not useful for the convergence of the algorithm [13].In addition, the position updating of one firefly is affected by other brighter fireflies, some of which may not be necessarily useful for solving the problem, and this may make the new location of the firefly far away from the optimal solution (we call this behavior invalid flight). Although invalid flight may be improved by the influence of other fireflies in the following movements, it will slow down the search process and reduce the speed of convergence. Thus, in order to improve the ability of basic FA to solve practical problems, some changes are necessary.


## 3. A Modified Firefly Algorithm

There are two important parameters in the formula for the position updating: the absorption coefficientγ and the random coefficient α. The smaller the γ, the greater the attraction between two fireflies and the faster the convergence speed. The range of random movement is larger and the convergence speed is slower, with a bigger α. In practical optimization problems, generally global search helps the algorithm converge to an area quickly and then a local search helps obtain a high precision solution. Reference [26] proposed an adaptive FA in which α and γ were adjusted with the iteration times, and its performance was better than that of the basic FA. Because the random movement step is 1 and the changing range of the random part is −0.5 to 0.5, the convergence speed is very slow in solving problems with a large range. In order to solve this problem, we proposed a new FA in which the random movement step was designed to be the distance between two fireflies. The position updating formula of the new FA is as follows:
(3)x→it+1=x→it+βijrij·x→it-x→jt+αt·rij·ε→i;βij(rij), α(t), and γ(t) can be obtained as follows:
(4)βijrij=β·e-γtrij2,γt=γb+tN·γe-γb,αt=αb+tN·αe-αb,
where x→i and x→j stand for the spatial positions of fireflies i and j; N is the whole iteration time; γe>γb, where γb and γe are the initial and ultimate value of γ, respectively; αe<αb, where αb and αe are the initial and ultimate value of α, respectively.As seen in formula3, the attraction between two fireflies is small when the distance between them is great. Attraction therefore does not play an important role in position updating. However, as the random movement step is big, fireflies move within a big range to explore a larger search space. When the distance among fireflies is small, the random movement step reduces and the random movement range is small. Attraction therefore begins to play a more important role in position updating. Fireflies move to the brightest one quickly and the algorithm comes to convergence at a great speed.In addition, in order to avoid invalid flight, we introduced the autonomous flight strategy (AFS). The new position of a firefly should be judged before it moves. The firefly will move and update the position if the new position is useful for solving the problem. Otherwise, the firefly will stay at the current position.


## 4. 3D Path Planning for AUV Using the Modified Firefly Algorithm

### 4.1. Coding of FA

In order to use FA to solve 3D path planning problem, the coding of FA should be done first. If a path hasn route points, it must be constituted by n-1 path segments. We assume that each firefly represents a candidate path in the firefly group, and the dimensions of position for each firefly correspond to the route points. So the number of fireflies corresponds to the candidate path number. So if a firefly group contains m members and each firefly in the group has n dimensions, the (i)th candidate path Path(i) represented by the (i)th firefly can be described as follows:
(5)Pathi=Pi,j∣j=1,…,n,
where Pi,j=(xi,j,yi,j,zi,j) stands for the (j)th route point and it corresponds to the (j)th dimension of the position for the (i)th firefly.In addition, as the brightness represents the quality of a firefly, it corresponds to the path quality in the path planning problem. The brightest firefly represents the optimal path.


### 4.2. Cost Function

There are many criteria for evaluating a candidate path. Generally, the smaller the path length is, the better it is. Additionally, in order to ensure the safety of AUV, the path should be away from obstacles, so we used the path lengthEL and risk value ED to evaluate the candidate path. The cost function of path planning can be calculated as follows:
(6)E=ω1EL+ω2ED,
where ω1+ω2=1.So the aim of path planning is to find a path which can satisfymin⁡(E). The calculating methods of every cost part are as follows.


#### 4.2.1. Path Length

If a path hasn route points, its length can be calculated by the following formula:
(7)EL=LSP1+∑i=1n-1ΔLi+LPnE,
where (xi,yi,zi),1≤i<n, stands for the coordinates of the (i)th route point,LSP1stands for the distance between the start point and the first route point, ΔLi=(xi-xi+1)2+(yi-yi+1)2+(zi-zi+1)2 represents the distance between the(i)th route point and the (i+1)th route point, and LPnE stands for the distance between the destination point and the last route point.


#### 4.2.2. Risk Value

A path holds risks if it passes through obstacles. So the risk value can be judged by how much of a path contains obstacles. In order to calculate the risk value easily, we selectedm evaluating points in every path segment. The whole risk value can therefore be judged by calculating how many route points and evaluating points which are in obstacles. For example, if there are k1 route points and k2 evaluating points in obstacles, the risk value can be calculated as follows:
(8)ED=Ck1+k2,
where C>0 stands for the cost coefficient.


### 4.3. Optimization Operator

#### 4.3.1. Excluding Operator

As the path planning space is 3D and contains many obstacles, it is inevitable that some points on the path will come into contact with obstacles in the process of searching. Sometimes these points would come out of the obstacles at last, but it will cost much time. In order to reduce searching time and increase the success rate of obstacle avoidance, it is necessary to help these points avoid obstacles quickly. An excluding operator is used here to make this easier. A path point inside an obstacle will be excluded to the surface of obstacle’s plane [11]. Figure 2 shows how the excluding operator works.Figure 2
Sketch map of excluding operator.In Figure2, ABC is a plane on the surface of one obstacle and n→ stands for the normal direction of ABC. Pi is a path point inside the obstacle. In order to ensure safety, Pi should move out of the obstacle using an excluding operator. Pi′ is on plane ABC which stands for the new position of Pi after it is excluded from obstacles. The position of Pi′ and the displacement of Pi can be calculated using the following formula:
(9)Pi′=Pi+PiPi′→,PiPi′→=PiPi′→•n→n→.


#### 4.3.2. Contracting Operator

The generation of a safe path should be done through the cooperation of every path point. Without taking path knowledge into account, the planning process will take a long time, especially in path planning with many dimensions. A path generated by FA only may be unsmooth, so a contracting operator is used to manage the cooperation of route points and to make the path smoother [11]. Figure 3 shows how the contracting operator works.Figure 3
Sketch map of contracting operator.In Figure3, Pi-1′ is the route point before Pi and Pi+1′ is the route point after Pi. On the assumption that Pi-1′ and Pi+1′ are fixed, Pi will be pulled by these two points and the direction of the resultant force is PiPi-1′→+PiPi+1′→. Theoretically, with the pulling force, Pi will move to Pi′′. In fact, because it is affected by air resistance, Pi can only move to Pi′. The displacement of Pi and its real position can be calculated by the following formula:
(10)Pi′=Pi+SPi→,SPi→=kPiPi-1′→+PiPi+1′→,
where SPi→ is the displacement vector of Pi with pulling forces from Pi-1′ and Pi+1′ and k∈(0,1) is the resistance coefficient.The path planning flow based on the modified FA is as follows: firstly, initializing path parameters and modeling the operating environment and then searching for the optimum path using modified FA, using an excluding operator to avoid obstacles and a contracting operator to optimize the path. These steps should be repeated until a satisfactory path is found. The detailed workflow is shown in Figure4.Figure 4
The flowchart for 3D path planning.


## 4.1. Coding of FA

In order to use FA to solve 3D path planning problem, the coding of FA should be done first. If a path hasn route points, it must be constituted by n-1 path segments. We assume that each firefly represents a candidate path in the firefly group, and the dimensions of position for each firefly correspond to the route points. So the number of fireflies corresponds to the candidate path number. So if a firefly group contains m members and each firefly in the group has n dimensions, the (i)th candidate path Path(i) represented by the (i)th firefly can be described as follows:
(5)Pathi=Pi,j∣j=1,…,n,
where Pi,j=(xi,j,yi,j,zi,j) stands for the (j)th route point and it corresponds to the (j)th dimension of the position for the (i)th firefly.In addition, as the brightness represents the quality of a firefly, it corresponds to the path quality in the path planning problem. The brightest firefly represents the optimal path.


## 4.2. Cost Function

There are many criteria for evaluating a candidate path. Generally, the smaller the path length is, the better it is. Additionally, in order to ensure the safety of AUV, the path should be away from obstacles, so we used the path lengthEL and risk value ED to evaluate the candidate path. The cost function of path planning can be calculated as follows:
(6)E=ω1EL+ω2ED,
where ω1+ω2=1.So the aim of path planning is to find a path which can satisfymin⁡(E). The calculating methods of every cost part are as follows.


### 4.2.1. Path Length

If a path hasn route points, its length can be calculated by the following formula:
(7)EL=LSP1+∑i=1n-1ΔLi+LPnE,
where (xi,yi,zi),1≤i<n, stands for the coordinates of the (i)th route point,LSP1stands for the distance between the start point and the first route point, ΔLi=(xi-xi+1)2+(yi-yi+1)2+(zi-zi+1)2 represents the distance between the(i)th route point and the (i+1)th route point, and LPnE stands for the distance between the destination point and the last route point.


### 4.2.2. Risk Value

A path holds risks if it passes through obstacles. So the risk value can be judged by how much of a path contains obstacles. In order to calculate the risk value easily, we selectedm evaluating points in every path segment. The whole risk value can therefore be judged by calculating how many route points and evaluating points which are in obstacles. For example, if there are k1 route points and k2 evaluating points in obstacles, the risk value can be calculated as follows:
(8)ED=Ck1+k2,
where C>0 stands for the cost coefficient.


## 4.2.1. Path Length

If a path hasn route points, its length can be calculated by the following formula:
(7)EL=LSP1+∑i=1n-1ΔLi+LPnE,
where (xi,yi,zi),1≤i<n, stands for the coordinates of the (i)th route point,LSP1stands for the distance between the start point and the first route point, ΔLi=(xi-xi+1)2+(yi-yi+1)2+(zi-zi+1)2 represents the distance between the(i)th route point and the (i+1)th route point, and LPnE stands for the distance between the destination point and the last route point.


## 4.2.2. Risk Value

A path holds risks if it passes through obstacles. So the risk value can be judged by how much of a path contains obstacles. In order to calculate the risk value easily, we selectedm evaluating points in every path segment. The whole risk value can therefore be judged by calculating how many route points and evaluating points which are in obstacles. For example, if there are k1 route points and k2 evaluating points in obstacles, the risk value can be calculated as follows:
(8)ED=Ck1+k2,
where C>0 stands for the cost coefficient.


## 4.3. Optimization Operator

### 4.3.1. Excluding Operator

As the path planning space is 3D and contains many obstacles, it is inevitable that some points on the path will come into contact with obstacles in the process of searching. Sometimes these points would come out of the obstacles at last, but it will cost much time. In order to reduce searching time and increase the success rate of obstacle avoidance, it is necessary to help these points avoid obstacles quickly. An excluding operator is used here to make this easier. A path point inside an obstacle will be excluded to the surface of obstacle’s plane [11]. Figure 2 shows how the excluding operator works.Figure 2
Sketch map of excluding operator.In Figure2, ABC is a plane on the surface of one obstacle and n→ stands for the normal direction of ABC. Pi is a path point inside the obstacle. In order to ensure safety, Pi should move out of the obstacle using an excluding operator. Pi′ is on plane ABC which stands for the new position of Pi after it is excluded from obstacles. The position of Pi′ and the displacement of Pi can be calculated using the following formula:
(9)Pi′=Pi+PiPi′→,PiPi′→=PiPi′→•n→n→.


### 4.3.2. Contracting Operator

The generation of a safe path should be done through the cooperation of every path point. Without taking path knowledge into account, the planning process will take a long time, especially in path planning with many dimensions. A path generated by FA only may be unsmooth, so a contracting operator is used to manage the cooperation of route points and to make the path smoother [11]. Figure 3 shows how the contracting operator works.Figure 3
Sketch map of contracting operator.In Figure3, Pi-1′ is the route point before Pi and Pi+1′ is the route point after Pi. On the assumption that Pi-1′ and Pi+1′ are fixed, Pi will be pulled by these two points and the direction of the resultant force is PiPi-1′→+PiPi+1′→. Theoretically, with the pulling force, Pi will move to Pi′′. In fact, because it is affected by air resistance, Pi can only move to Pi′. The displacement of Pi and its real position can be calculated by the following formula:
(10)Pi′=Pi+SPi→,SPi→=kPiPi-1′→+PiPi+1′→,
where SPi→ is the displacement vector of Pi with pulling forces from Pi-1′ and Pi+1′ and k∈(0,1) is the resistance coefficient.The path planning flow based on the modified FA is as follows: firstly, initializing path parameters and modeling the operating environment and then searching for the optimum path using modified FA, using an excluding operator to avoid obstacles and a contracting operator to optimize the path. These steps should be repeated until a satisfactory path is found. The detailed workflow is shown in Figure4.Figure 4
The flowchart for 3D path planning.


## 4.3.1. Excluding Operator

As the path planning space is 3D and contains many obstacles, it is inevitable that some points on the path will come into contact with obstacles in the process of searching. Sometimes these points would come out of the obstacles at last, but it will cost much time. In order to reduce searching time and increase the success rate of obstacle avoidance, it is necessary to help these points avoid obstacles quickly. An excluding operator is used here to make this easier. A path point inside an obstacle will be excluded to the surface of obstacle’s plane [11]. Figure 2 shows how the excluding operator works.Figure 2
Sketch map of excluding operator.In Figure2, ABC is a plane on the surface of one obstacle and n→ stands for the normal direction of ABC. Pi is a path point inside the obstacle. In order to ensure safety, Pi should move out of the obstacle using an excluding operator. Pi′ is on plane ABC which stands for the new position of Pi after it is excluded from obstacles. The position of Pi′ and the displacement of Pi can be calculated using the following formula:
(9)Pi′=Pi+PiPi′→,PiPi′→=PiPi′→•n→n→.


## 4.3.2. Contracting Operator

The generation of a safe path should be done through the cooperation of every path point. Without taking path knowledge into account, the planning process will take a long time, especially in path planning with many dimensions. A path generated by FA only may be unsmooth, so a contracting operator is used to manage the cooperation of route points and to make the path smoother [11]. Figure 3 shows how the contracting operator works.Figure 3
Sketch map of contracting operator.In Figure3, Pi-1′ is the route point before Pi and Pi+1′ is the route point after Pi. On the assumption that Pi-1′ and Pi+1′ are fixed, Pi will be pulled by these two points and the direction of the resultant force is PiPi-1′→+PiPi+1′→. Theoretically, with the pulling force, Pi will move to Pi′′. In fact, because it is affected by air resistance, Pi can only move to Pi′. The displacement of Pi and its real position can be calculated by the following formula:
(10)Pi′=Pi+SPi→,SPi→=kPiPi-1′→+PiPi+1′→,
where SPi→ is the displacement vector of Pi with pulling forces from Pi-1′ and Pi+1′ and k∈(0,1) is the resistance coefficient.The path planning flow based on the modified FA is as follows: firstly, initializing path parameters and modeling the operating environment and then searching for the optimum path using modified FA, using an excluding operator to avoid obstacles and a contracting operator to optimize the path. These steps should be repeated until a satisfactory path is found. The detailed workflow is shown in Figure4.Figure 4
The flowchart for 3D path planning.


## 5. Experimental Setup

### 5.1. Experiments Using Modified FA

Genetic algorithm (GA) and particle swarm optimization (PSO) are two classic algorithms with extreme influence in optimization field, and research on them is relatively mature. FA has some qualities in common with these two algorithms. In order to verify the performance of the modified FA (AMFA) proposed in this paper, we used four algorithms (basic FA, AMFA, PSO, and GA) to solve six classic test functions. The four classic test functions are as follows.Sphere:f1(x→)=∑i=1nxi2; the function has its global optimum f*=0 at s*=(0,…,0).Rosenbrock:f2(x→)=∑i=1n-1[100(xi2-xi+1)2+xi-1)2;, the function has global optimum f*=0 at s*=(1,…,1).Rastrigin:f3(x→)=∑i=1nxi2-10cos⁡2πxi+10;, the function has global optimum f*=0 at s*=(0,…,0).Griewank:f4(x→)=∑i=1n-1[100(xi2-xi+1)2+(xi-1)2]; the function has global optimum f*=0 at s*=(0,…,0).Ackley:f5=-20exp⁡[-0.21/d∑i=1dxi2]-exp⁡[1/d∑i=1dcos⁡(2πxi)]+(20+e); the function has global optimum f*=0 at s*=(0,…,0).Zakharo:f6=∑i=1Dsi2+(1/2∑i=1Disi)2+(1/2∑i=1Disi)4; the function has global optimum f*=0 at s*=(0,…,0).Parameters of these four algorithms are set as follows.FA:β0=1.0, γ=1.0, and α=1.0. rand was the uniformly distributed random number in the range of 0-1.GA: binary encoding was used and the binary digit was 20. The roulette wheel selection method and single point crossover were used. Crossover probability was defined aspc=0.6, mutation probability was described as pm=1/n, and n was the dimension of the optimization function.PSO: using parameters recommended by Li et al. [10], ω=0.729 and c1=c2=1.494. The maximum speed was set to the top of the search range, Vmax⁡=Xmax⁡.AMFA:β0=1.0, γe=2.0, γb=1.0, αe=0.8, and αb=2.0. rand was the uniformly distributed random number in the range of 0-1.In order to compare the optimization performance of these algorithms fairly and sufficiently, they were used to solve the four test functions in both small and wide ranges. The population size of every algorithm was set to 40 and the dimension of the testing function was set to 30. Each algorithm was run 30 times independently, and the number of fitness evaluations was limited to 60,000.


#### 5.1.1. Experiments in a Small Range

The range of every test function is shown in Table1. Table 2 shows the results of simulation for every algorithm. In order to observe the iterative process of the algorithm, the curves of the six test functions are shown in Figures 5 to 10.Table 1
Test functions.
Function
Sphere
Rosenbrock
Rastrigin
Griewank
Ackley
Zakharo
Range
[−5.12, 5.12]
[−2.048, 2.048]
[−5.12, 5.12]
[−8, 8]
[−2.768, 2.768]
[−5, 10]Table 2
Experimental results in a small range.
Function
Term
AMFA
FA
PSO
GA
Sphere
Best
0
5.20E − 09
5.05E − 41
9.64E − 07
Mean
5
.58E
−
13
2.76E
−
06
1.75E
−
36
4.91E − 06
STD
1.73E − 12
3.58E − 06
9.19E − 36
4.56E − 06
“0” times
22
0
0
0
Rosenbrock
Best
0
3.66E − 07
1.27E − 02
7.8507
Mean
2.38E − 06
6.70E − 05
2.17E − 01
51.3280
STD
4.34E − 07
1.07E − 04
0.75
29.7729
“0” times
21
0
0
0
Rastrigin
Best
00E + 00
1.98E − 07
1.00
43.78
Mean
00E + 00
9.81E − 06
4.88
58.04
STD
00E + 00
1.70E − 04
2.24
8.28
“0” times
30
0
0
0
Griewank
Best
00E + 00
1.35E − 10
00E + 00
1.14E − 07
Mean
3.20E − 11
5.36E − 07
5.6E − 03
2.30E − 02
STD
1.12E − 10
8.51E − 07
8.3E − 03
4.40E − 02
“0” times
26
0
7
0
Ackley
Best
0
3.16E − 05
1.42E − 14
1.34E + 00
Mean
1.75E − 07
3.12E − 04
3.02E − 01
2.34E + 00
STD
7.23E − 07
2.69E − 04
5.68E − 01
3.71E − 01
“0” times
21
0
0
0
Zakharo
Best
00E + 00
1.7E − 07
4.6E − 06
9.8E + 01
Mean
1.6E − 06
2.4E − 03
6.7E − 05
1.8E + 02
STD
9.0E − 06
3.8E − 03
6.8E − 05
4.6E + 01
“0” times
4
0
0
0Figure 5
Test results of Sphere.Experimental results in Table2 show that the precision of solution obtained by AMFA is much better than that of FA. Meanwhile, Figures 5–10 show that the convergence speed of AMFA is quicker than that of FA. So we can conclude that the performance of AMFA is better than that of FA and the method proposed here is effective in improving the performance of basic FA. The performance and convergence speed of AMFA are better than those of PSO and GA. Although the solution precision of Sphere obtained by AMFA is not the best, it is valid. Meanwhile, the solution of Sphere obtained by AMFA is found to be zero 22 times. We can therefore conclude that AMFA is an effective algorithm for solving problems with a small range.Figure 6
Test results of Rosenbrock.Figure 7
Test results of Rastrigin.Figure 8
Test results of Griewank.Figure 9
Test results of Ackley.Figure 10
Test results of Zakharo.


#### 5.1.2. Experiments in a Wide Range

The range of every test function is shown in Table3. Table 4 shows the results of simulation for every algorithm. In order to observe the iterative process of the algorithm, the curves of the six test functions are shown in Figures 11–16.Table 3
Test functions.
Function
Sphere
Rosenbrock
Rastrigin
Griewank
Ackley
Zakharo
Range
[−100, 100]
[−100, 100]
[−100, 100]
[−600, 600]
[−23.768, 23.768]
[−50, 100]Table 4
Experimental results in wide range.
Function
Term
AMFA
FA
PSO
GA
Sphere
Best
00E + 00
5.87E − 09
5.23E − 39
3.42E − 04
Mean
0
0E
+
00
4.00E
−
03
8.73E
−
35
1.71E
−
03
STD
00E + 00
2.15E − 02
2.16E − 34
1.52E − 03
“0” times
30
0
0
0
Rosenbrock
Best
00E + 00
5.95E − 07
2.58E − 02
2.88E + 00
Mean
00E + 00
6.00E + 04
8.61E − 01
7.67E + 01
STD
00E + 00
3.28E + 05
1.46E + 00
4.98E + 01
“0” times
30
0
0
0
Rastrigin
Best
00E + 00
6.11E − 05
2.98E + 00
1.99E + 01
Mean
00E + 00
8.13E + 01
6.40E + 00
6.08E + 01
STD
00E + 00
2.54E + 02
2.38E + 00
3.15E + 01
“0” times
30
0
0
0
Griewank
Best
00E + 00
4.05E + 02
00E + 00
4.85E − 04
Mean
00E + 00
5.72E + 02
2.30E − 02
5.45E − 02
STD
00E + 00
5.52E + 01
2.85E − 02
5.30E − 02
“0” times
30
0
2
0
Ackley
Best
0
3.57E − 05
2.13E − 14
4.10E − 03
Mean
00E + 00
2.30E − 03
1.63E + 00
7.7E − 03
STD
0
4.40E − 03
1.12E + 00
3.1E − 03
“0” times
30
0
0
0
Zakharo
Best
0
2.0E − 06
4.8E + 01
1.5E + 04
Mean
1.5E − 322
1.5E − 01
3.3E + 02
3.0E + 04
STD
0
7.6E − 01
2.8E + 02
6.2E + 03
“0” times
5
0
0
0Figure 11
Test results of Sphere.Figure 12
Test results of Rosenbrock.Figure 13
Test results of Rastrigin.Figure 14
Test results of Griewank.Figure 15
Test results of Ackley.Figure 16
Test results of Zakharo.Experimental results in Table4 show that the precision of solution obtained by AMFA is the most effective for solving all six test functions. Additionally, it can find the global optimum every time. Figures 11 to 16 show that the convergence speed of AMFA is the quickest of the four algorithms. So we can conclude that the method proposed here is effective in improving the performance of basic FA and AMFA is powerful in solving problems with a wide range.


### 5.2. Experiments of Path Planning

In order to simulate the real sailing environmental space of AUV, a maritime space with a horizontal area about2°×2° was selected as the experiment environment for AUV path planning. The terrain was generated by interpolation with elevation data extracted from electric chart. In order to observe the obstacle avoidance capability of the method, the actual lengths of the horizontal and vertical coordinates were reduced (1852×60)/1000 times. The resistance coefficient was 0.9553 and parameters of AMFA were the same as in the experiments described above. Because the population size of AMFA was 40 and the dimension of every firefly was set to 30, the number of candidate paths was 40 and every path should have 30 path points.In order to verify the performance of the new algorithm, two experiments were designed here. In the first experiment, depths of the starting point and destination are both 300 meters. The horizontal position of starting point is (111.10°E, 17.07°N), and the horizontal position of destination is (112.50°E, 16.69°N). The results of the path planning are shown in Figure17. In the second experiment, the start position and destination have different depths. The horizontal position of the starting point is (111.30°E, 15.87°N) and the depth is 230 meters. The horizontal position of the destination is (112.80°E, 17.11°N), and the depth is 450 meters. The results of the path planning are shown in Figure 18.Figure 17
Results of experiment 1.Figure 18
Results of experiment 2.Figures17 and 18 show that there are many obstacles in the environmental space. The line represents the path found by AMFA, the black dot is the starting point, and the red dot stands for the destination in the figures. From the results of the experiments we can conclude that the path planning method using AMFA can find the optimal path in the complex environmental space, and the paths in both experiments are smooth and do not path through obstacles. The experimental results prove that the path planning method proposed here is feasible and AMFA can be used in solving 3D path planning problem.


## 5.1. Experiments Using Modified FA

Genetic algorithm (GA) and particle swarm optimization (PSO) are two classic algorithms with extreme influence in optimization field, and research on them is relatively mature. FA has some qualities in common with these two algorithms. In order to verify the performance of the modified FA (AMFA) proposed in this paper, we used four algorithms (basic FA, AMFA, PSO, and GA) to solve six classic test functions. The four classic test functions are as follows.Sphere:f1(x→)=∑i=1nxi2; the function has its global optimum f*=0 at s*=(0,…,0).Rosenbrock:f2(x→)=∑i=1n-1[100(xi2-xi+1)2+xi-1)2;, the function has global optimum f*=0 at s*=(1,…,1).Rastrigin:f3(x→)=∑i=1nxi2-10cos⁡2πxi+10;, the function has global optimum f*=0 at s*=(0,…,0).Griewank:f4(x→)=∑i=1n-1[100(xi2-xi+1)2+(xi-1)2]; the function has global optimum f*=0 at s*=(0,…,0).Ackley:f5=-20exp⁡[-0.21/d∑i=1dxi2]-exp⁡[1/d∑i=1dcos⁡(2πxi)]+(20+e); the function has global optimum f*=0 at s*=(0,…,0).Zakharo:f6=∑i=1Dsi2+(1/2∑i=1Disi)2+(1/2∑i=1Disi)4; the function has global optimum f*=0 at s*=(0,…,0).Parameters of these four algorithms are set as follows.FA:β0=1.0, γ=1.0, and α=1.0. rand was the uniformly distributed random number in the range of 0-1.GA: binary encoding was used and the binary digit was 20. The roulette wheel selection method and single point crossover were used. Crossover probability was defined aspc=0.6, mutation probability was described as pm=1/n, and n was the dimension of the optimization function.PSO: using parameters recommended by Li et al. [10], ω=0.729 and c1=c2=1.494. The maximum speed was set to the top of the search range, Vmax⁡=Xmax⁡.AMFA:β0=1.0, γe=2.0, γb=1.0, αe=0.8, and αb=2.0. rand was the uniformly distributed random number in the range of 0-1.In order to compare the optimization performance of these algorithms fairly and sufficiently, they were used to solve the four test functions in both small and wide ranges. The population size of every algorithm was set to 40 and the dimension of the testing function was set to 30. Each algorithm was run 30 times independently, and the number of fitness evaluations was limited to 60,000.


### 5.1.1. Experiments in a Small Range

The range of every test function is shown in Table1. Table 2 shows the results of simulation for every algorithm. In order to observe the iterative process of the algorithm, the curves of the six test functions are shown in Figures 5 to 10.Table 1
Test functions.
Function
Sphere
Rosenbrock
Rastrigin
Griewank
Ackley
Zakharo
Range
[−5.12, 5.12]
[−2.048, 2.048]
[−5.12, 5.12]
[−8, 8]
[−2.768, 2.768]
[−5, 10]Table 2
Experimental results in a small range.
Function
Term
AMFA
FA
PSO
GA
Sphere
Best
0
5.20E − 09
5.05E − 41
9.64E − 07
Mean
5
.58E
−
13
2.76E
−
06
1.75E
−
36
4.91E − 06
STD
1.73E − 12
3.58E − 06
9.19E − 36
4.56E − 06
“0” times
22
0
0
0
Rosenbrock
Best
0
3.66E − 07
1.27E − 02
7.8507
Mean
2.38E − 06
6.70E − 05
2.17E − 01
51.3280
STD
4.34E − 07
1.07E − 04
0.75
29.7729
“0” times
21
0
0
0
Rastrigin
Best
00E + 00
1.98E − 07
1.00
43.78
Mean
00E + 00
9.81E − 06
4.88
58.04
STD
00E + 00
1.70E − 04
2.24
8.28
“0” times
30
0
0
0
Griewank
Best
00E + 00
1.35E − 10
00E + 00
1.14E − 07
Mean
3.20E − 11
5.36E − 07
5.6E − 03
2.30E − 02
STD
1.12E − 10
8.51E − 07
8.3E − 03
4.40E − 02
“0” times
26
0
7
0
Ackley
Best
0
3.16E − 05
1.42E − 14
1.34E + 00
Mean
1.75E − 07
3.12E − 04
3.02E − 01
2.34E + 00
STD
7.23E − 07
2.69E − 04
5.68E − 01
3.71E − 01
“0” times
21
0
0
0
Zakharo
Best
00E + 00
1.7E − 07
4.6E − 06
9.8E + 01
Mean
1.6E − 06
2.4E − 03
6.7E − 05
1.8E + 02
STD
9.0E − 06
3.8E − 03
6.8E − 05
4.6E + 01
“0” times
4
0
0
0Figure 5
Test results of Sphere.Experimental results in Table2 show that the precision of solution obtained by AMFA is much better than that of FA. Meanwhile, Figures 5–10 show that the convergence speed of AMFA is quicker than that of FA. So we can conclude that the performance of AMFA is better than that of FA and the method proposed here is effective in improving the performance of basic FA. The performance and convergence speed of AMFA are better than those of PSO and GA. Although the solution precision of Sphere obtained by AMFA is not the best, it is valid. Meanwhile, the solution of Sphere obtained by AMFA is found to be zero 22 times. We can therefore conclude that AMFA is an effective algorithm for solving problems with a small range.Figure 6
Test results of Rosenbrock.Figure 7
Test results of Rastrigin.Figure 8
Test results of Griewank.Figure 9
Test results of Ackley.Figure 10
Test results of Zakharo.


### 5.1.2. Experiments in a Wide Range

The range of every test function is shown in Table3. Table 4 shows the results of simulation for every algorithm. In order to observe the iterative process of the algorithm, the curves of the six test functions are shown in Figures 11–16.Table 3
Test functions.
Function
Sphere
Rosenbrock
Rastrigin
Griewank
Ackley
Zakharo
Range
[−100, 100]
[−100, 100]
[−100, 100]
[−600, 600]
[−23.768, 23.768]
[−50, 100]Table 4
Experimental results in wide range.
Function
Term
AMFA
FA
PSO
GA
Sphere
Best
00E + 00
5.87E − 09
5.23E − 39
3.42E − 04
Mean
0
0E
+
00
4.00E
−
03
8.73E
−
35
1.71E
−
03
STD
00E + 00
2.15E − 02
2.16E − 34
1.52E − 03
“0” times
30
0
0
0
Rosenbrock
Best
00E + 00
5.95E − 07
2.58E − 02
2.88E + 00
Mean
00E + 00
6.00E + 04
8.61E − 01
7.67E + 01
STD
00E + 00
3.28E + 05
1.46E + 00
4.98E + 01
“0” times
30
0
0
0
Rastrigin
Best
00E + 00
6.11E − 05
2.98E + 00
1.99E + 01
Mean
00E + 00
8.13E + 01
6.40E + 00
6.08E + 01
STD
00E + 00
2.54E + 02
2.38E + 00
3.15E + 01
“0” times
30
0
0
0
Griewank
Best
00E + 00
4.05E + 02
00E + 00
4.85E − 04
Mean
00E + 00
5.72E + 02
2.30E − 02
5.45E − 02
STD
00E + 00
5.52E + 01
2.85E − 02
5.30E − 02
“0” times
30
0
2
0
Ackley
Best
0
3.57E − 05
2.13E − 14
4.10E − 03
Mean
00E + 00
2.30E − 03
1.63E + 00
7.7E − 03
STD
0
4.40E − 03
1.12E + 00
3.1E − 03
“0” times
30
0
0
0
Zakharo
Best
0
2.0E − 06
4.8E + 01
1.5E + 04
Mean
1.5E − 322
1.5E − 01
3.3E + 02
3.0E + 04
STD
0
7.6E − 01
2.8E + 02
6.2E + 03
“0” times
5
0
0
0Figure 11
Test results of Sphere.Figure 12
Test results of Rosenbrock.Figure 13
Test results of Rastrigin.Figure 14
Test results of Griewank.Figure 15
Test results of Ackley.Figure 16
Test results of Zakharo.Experimental results in Table4 show that the precision of solution obtained by AMFA is the most effective for solving all six test functions. Additionally, it can find the global optimum every time. Figures 11 to 16 show that the convergence speed of AMFA is the quickest of the four algorithms. So we can conclude that the method proposed here is effective in improving the performance of basic FA and AMFA is powerful in solving problems with a wide range.


## 5.1.1. Experiments in a Small Range

The range of every test function is shown in Table1. Table 2 shows the results of simulation for every algorithm. In order to observe the iterative process of the algorithm, the curves of the six test functions are shown in Figures 5 to 10.Table 1
Test functions.
Function
Sphere
Rosenbrock
Rastrigin
Griewank
Ackley
Zakharo
Range
[−5.12, 5.12]
[−2.048, 2.048]
[−5.12, 5.12]
[−8, 8]
[−2.768, 2.768]
[−5, 10]Table 2
Experimental results in a small range.
Function
Term
AMFA
FA
PSO
GA
Sphere
Best
0
5.20E − 09
5.05E − 41
9.64E − 07
Mean
5
.58E
−
13
2.76E
−
06
1.75E
−
36
4.91E − 06
STD
1.73E − 12
3.58E − 06
9.19E − 36
4.56E − 06
“0” times
22
0
0
0
Rosenbrock
Best
0
3.66E − 07
1.27E − 02
7.8507
Mean
2.38E − 06
6.70E − 05
2.17E − 01
51.3280
STD
4.34E − 07
1.07E − 04
0.75
29.7729
“0” times
21
0
0
0
Rastrigin
Best
00E + 00
1.98E − 07
1.00
43.78
Mean
00E + 00
9.81E − 06
4.88
58.04
STD
00E + 00
1.70E − 04
2.24
8.28
“0” times
30
0
0
0
Griewank
Best
00E + 00
1.35E − 10
00E + 00
1.14E − 07
Mean
3.20E − 11
5.36E − 07
5.6E − 03
2.30E − 02
STD
1.12E − 10
8.51E − 07
8.3E − 03
4.40E − 02
“0” times
26
0
7
0
Ackley
Best
0
3.16E − 05
1.42E − 14
1.34E + 00
Mean
1.75E − 07
3.12E − 04
3.02E − 01
2.34E + 00
STD
7.23E − 07
2.69E − 04
5.68E − 01
3.71E − 01
“0” times
21
0
0
0
Zakharo
Best
00E + 00
1.7E − 07
4.6E − 06
9.8E + 01
Mean
1.6E − 06
2.4E − 03
6.7E − 05
1.8E + 02
STD
9.0E − 06
3.8E − 03
6.8E − 05
4.6E + 01
“0” times
4
0
0
0Figure 5
Test results of Sphere.Experimental results in Table2 show that the precision of solution obtained by AMFA is much better than that of FA. Meanwhile, Figures 5–10 show that the convergence speed of AMFA is quicker than that of FA. So we can conclude that the performance of AMFA is better than that of FA and the method proposed here is effective in improving the performance of basic FA. The performance and convergence speed of AMFA are better than those of PSO and GA. Although the solution precision of Sphere obtained by AMFA is not the best, it is valid. Meanwhile, the solution of Sphere obtained by AMFA is found to be zero 22 times. We can therefore conclude that AMFA is an effective algorithm for solving problems with a small range.Figure 6
Test results of Rosenbrock.Figure 7
Test results of Rastrigin.Figure 8
Test results of Griewank.Figure 9
Test results of Ackley.Figure 10
Test results of Zakharo.


## 5.1.2. Experiments in a Wide Range

The range of every test function is shown in Table3. Table 4 shows the results of simulation for every algorithm. In order to observe the iterative process of the algorithm, the curves of the six test functions are shown in Figures 11–16.Table 3
Test functions.
Function
Sphere
Rosenbrock
Rastrigin
Griewank
Ackley
Zakharo
Range
[−100, 100]
[−100, 100]
[−100, 100]
[−600, 600]
[−23.768, 23.768]
[−50, 100]Table 4
Experimental results in wide range.
Function
Term
AMFA
FA
PSO
GA
Sphere
Best
00E + 00
5.87E − 09
5.23E − 39
3.42E − 04
Mean
0
0E
+
00
4.00E
−
03
8.73E
−
35
1.71E
−
03
STD
00E + 00
2.15E − 02
2.16E − 34
1.52E − 03
“0” times
30
0
0
0
Rosenbrock
Best
00E + 00
5.95E − 07
2.58E − 02
2.88E + 00
Mean
00E + 00
6.00E + 04
8.61E − 01
7.67E + 01
STD
00E + 00
3.28E + 05
1.46E + 00
4.98E + 01
“0” times
30
0
0
0
Rastrigin
Best
00E + 00
6.11E − 05
2.98E + 00
1.99E + 01
Mean
00E + 00
8.13E + 01
6.40E + 00
6.08E + 01
STD
00E + 00
2.54E + 02
2.38E + 00
3.15E + 01
“0” times
30
0
0
0
Griewank
Best
00E + 00
4.05E + 02
00E + 00
4.85E − 04
Mean
00E + 00
5.72E + 02
2.30E − 02
5.45E − 02
STD
00E + 00
5.52E + 01
2.85E − 02
5.30E − 02
“0” times
30
0
2
0
Ackley
Best
0
3.57E − 05
2.13E − 14
4.10E − 03
Mean
00E + 00
2.30E − 03
1.63E + 00
7.7E − 03
STD
0
4.40E − 03
1.12E + 00
3.1E − 03
“0” times
30
0
0
0
Zakharo
Best
0
2.0E − 06
4.8E + 01
1.5E + 04
Mean
1.5E − 322
1.5E − 01
3.3E + 02
3.0E + 04
STD
0
7.6E − 01
2.8E + 02
6.2E + 03
“0” times
5
0
0
0Figure 11
Test results of Sphere.Figure 12
Test results of Rosenbrock.Figure 13
Test results of Rastrigin.Figure 14
Test results of Griewank.Figure 15
Test results of Ackley.Figure 16
Test results of Zakharo.Experimental results in Table4 show that the precision of solution obtained by AMFA is the most effective for solving all six test functions. Additionally, it can find the global optimum every time. Figures 11 to 16 show that the convergence speed of AMFA is the quickest of the four algorithms. So we can conclude that the method proposed here is effective in improving the performance of basic FA and AMFA is powerful in solving problems with a wide range.


## 5.2. Experiments of Path Planning

In order to simulate the real sailing environmental space of AUV, a maritime space with a horizontal area about2°×2° was selected as the experiment environment for AUV path planning. The terrain was generated by interpolation with elevation data extracted from electric chart. In order to observe the obstacle avoidance capability of the method, the actual lengths of the horizontal and vertical coordinates were reduced (1852×60)/1000 times. The resistance coefficient was 0.9553 and parameters of AMFA were the same as in the experiments described above. Because the population size of AMFA was 40 and the dimension of every firefly was set to 30, the number of candidate paths was 40 and every path should have 30 path points.In order to verify the performance of the new algorithm, two experiments were designed here. In the first experiment, depths of the starting point and destination are both 300 meters. The horizontal position of starting point is (111.10°E, 17.07°N), and the horizontal position of destination is (112.50°E, 16.69°N). The results of the path planning are shown in Figure17. In the second experiment, the start position and destination have different depths. The horizontal position of the starting point is (111.30°E, 15.87°N) and the depth is 230 meters. The horizontal position of the destination is (112.80°E, 17.11°N), and the depth is 450 meters. The results of the path planning are shown in Figure 18.Figure 17
Results of experiment 1.Figure 18
Results of experiment 2.Figures17 and 18 show that there are many obstacles in the environmental space. The line represents the path found by AMFA, the black dot is the starting point, and the red dot stands for the destination in the figures. From the results of the experiments we can conclude that the path planning method using AMFA can find the optimal path in the complex environmental space, and the paths in both experiments are smooth and do not path through obstacles. The experimental results prove that the path planning method proposed here is feasible and AMFA can be used in solving 3D path planning problem.


## 6. Conclusions

This paper presents a 3D path planning method based on a modified FA. The performance of the modified FA and the path planning method was verified through simulation. The simulation results show that the modified FA has a quick convergence speed and the path planning method based on it can find an effective path in a 3D environment. However, as a new SI-based algorithm, the firefly algorithm still has imperfections from a theoretical research point of view and its mathematical theory also needs to be studied and proved. Additionally, it has been assumed that the operating environment is static and obstacle information is known before the path planning method is begun. In reality, the marine environment is dynamic and AUVs will face many sudden threats while in action. Finally, many constraints for AUVs such as the maximal turning angle and the maximal sailing depth must be taken into account. Therefore further research is required to conclude how to plan a practical path for AUVs in the real dynamic marine environment.



---

*Source: 561394-2015-10-12.xml*