# Global Path Planning for Unmanned Surface Vehicle Based on Improved Quantum Ant Colony Algorithm

Authors: Guoqing Xia; Zhiwei Han; Bo Zhao; Caiyun Liu; Xinwei Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

As a tool to monitor marine environments and to perform dangerous tasks instead of manned vessels, unmanned surface vehicles (USVs) have extensive applications. Because most path planning algorithms have difficulty meeting the mission requirements of USVs, the purpose of this study was to plan a global path with multiple objectives, such as path length, energy consumption, path smoothness, and path safety, for USV in marine environments. A global path planning algorithm based on an improved quantum ant colony algorithm (IQACA) is proposed. The improved quantum ant colony algorithm is an algorithm that benefits from the high efficiency of quantum computing and the optimization ability of the ant colony algorithm. The proposed algorithm can plan a path considering multiple objectives simultaneously. The simulation results show that the proposed algorithm’s obtained minimum was 2.1–6.5% lower than those of the quantum ant colony algorithm (QACA) and ant colony algorithm (ACA), and the number of iterations required to converge to the minimum was 11.2–24.5% lower than those of the QACA and ACA. In addition, the optimized path for the USV was obtained effectively and efficiently.


---

## Body

## 1. Introduction

An unmanned surface vehicle (USV) is a kind of autonomous marine vehicle. Determining the path of a USV is an important problem associated with its safety and efficiency ‎[1]. Depending on whether the environmental information is obtained from a digital map or sensors, path planning is divided into global and local stages ‎[2]. In this paper, a USV global path planning study is presented. Global path planning is the process of planning a path to connect the starting and destination points under a given planning space from a digital map and constraints according to the mission requirements. The indices for evaluating a path can be path length, energy consumption, path smoothness, and path safety.Obtaining a short path from the starting point to the destination point is one of the main objectives of global path planning. Planning the shortest path is an NP-hard problem ‎[3]. Existing methods take the path length as a single objective of the path planning, and neither energy consumption nor other indices are considered.The energy consumption during sailing determines the USV’s endurance and the duration of the mission. Since the environmental loads such as wind, waves, and ocean currents influence the performance of the USV, the calculation of the USV’s energy consumption is complex. Niu et al. ‎[4] considered the effect of the ocean current on the energy consumption of USVs. Lee et al. ‎[5] found a more economical path by considering the shallow water effect as well as tidal currents and wind for surface ship navigation. Most calculations of energy consumption have considered the effects of ocean currents on the USV without considering wind and waves.The smoothness of a path depends on the size and number of the turns that the USV makes while sailing along the planned path. The smoother path allows the USV to make fewer turns along the path, which reduces the mechanical wear on the steering actuators, such as rudders. Smooth paths can reduce unnecessary curvature discontinuities and possible stops. In a previous report ‎[6], the smoothness of a path was evaluated by summing the angles of each turn on the path that the vehicle follows. Ma et al. ‎[7] evaluated the turn angle set by adopting the maximum value of the turn angle set to assess the path smoothness for the USV.Obstacles such as islands and reefs affect the safety of the USV. Path safety means that the USV cannot collide with any obstacles while sailing. Ma et al. ‎[7] used circles that just covered the obstacles to identify the safe area.Since USV global path planning involves optimization algorithms, environmental models, and marine craft hydrodynamics, existing path planning algorithms have difficulty meeting the mission requirements. Intelligent optimization algorithms are widely used in global path planning, such as the genetic algorithm ‎[8], particle swarm algorithm ‎[9], NSGA-II ‎[10], and ant colony algorithm ‎[11]. With the development of quantum technology, the idea of combining quantum computing with intelligent optimization algorithms has been developed. Narayanan and Moore combined quantum mechanics principles and evolutionary computing methods for the first time ‎[12]. A quantum bit and superposition of states were proposed to solve the knapsack problem by a quantum-inspired evolutionary algorithm (QEA) ‎[13]. Based on the QEA with a quantum rotation gate strategy, an adaptive evolution-based quantum-inspired evolutionary algorithm (AEQEA) introduces an adaptive evolution mechanism ‎[14]. A new improved quantum evolution algorithm (IQEA) with a mixed local search procedure was proposed ‎[15]. Li et al. ‎[16] proposed a quantum ant colony algorithm (QACA) that combined quantum computing and the ant colony algorithm for continuous space optimization. You et al. ‎[17] proposed a novel parallel ant colony optimization algorithm based on a quantum dynamics mechanism (PQACO). An improved quantum ant colony algorithm was proposed for the optimization of evacuation paths from dangerous areas to safe areas ‎[18]. The quantum ant colony algorithm was used to determine campus path navigation ‎[19].In this paper, a global path planning algorithm for USV based on the improved quantum ant colony algorithm (IQACA) is proposed. The main contributions of the proposed approach are as follows:(1)
At present, most USV global path planning algorithms only search for a feasible path for one objective ‎[3–6]. In this paper, path planning was considered with multiple simultaneous objectives, which were path length, energy consumption, path smoothness, and path safety.(2)
The IQACA is a new optimization algorithm that combines quantum-inspired computing with the ant colony algorithm (ACA). The quantum bit (Q-bit) is used to encode the pheromone in the ACA to obtain the quantum pheromone, and the ant movement is determined based on the concentration of the quantum pheromone on the path. Compared to the existing QACA [16–19], the phase of the quantum ant colony is transformed by an adaptive quantum rotation gate, and the quantum pheromone is updated by local and global update rules in the IQACA.Simulation experiments in a complex environment with wind, waves, and ocean currents verified the effectiveness of the objective model, and we obtained a desired path based on the IQACA.The paper is organized as follows. In Section2, the USV path planning problem is established, and the USV kinetic model, environmental loads, and cost function of the path planning are described. In Section 3, the principles of the IQACA are provided, and we apply the IQACA to USV global path planning. In Section 4, the simulations for USV global path planning using the IQACA are presented. Conclusions are provided in Section 5.


## 2. Problem Statement

### 2.1. USV Kinetic Model

The kinetic model of a USV accounts for the forces, such as the control force and environmental loads, which cause USV motion. For the USV, the control force is mainly the thrust of each propeller. The environmental loads on the USV are generated by wind, waves, and ocean currents. The kinetic model of the USV, which was proposed previously ‎[20], is as follows:(1)Mν˙+Cνν+Dνν=τenv+τthr(2)τenv=τwind+τwave+τcurrentwhere M is the system inertia matrix, C(ν) is the Coriolis-centripetal matrix, D(ν)∈R3×3 is damping matrix. τwind, τwave, and τcurrent are wind, wave, and ocean current forces acting on the USV, respectively, and τthr is the thrust generated by the USV propulsion system. The generalized velocity ν=[u,v,r]T is obtained by (1), where the first two components (u,v) are the linear velocities of the surge and sway, and r is the angular velocity of the yaw.


### 2.2. Models of Environmental Loads

When planning a global path for USVs, it is necessary to consider the environmental effects on the vehicles. Thus, we need to analyze the impacts of wind, waves, and ocean currents on the USV. The planned area is a confined sea with some static obstacles, and the mission execution time is short. Therefore, it can be assumed that the environmental loads are basically stable in limited time and space.


#### 2.2.1. Wind Forces

The wind acts directly on the superstructure of the hull. As reported previously ‎[21], the wind forces are written as follows:(3)τwindX=12ρaAfVw2CwxαRτwindY=12ρaAsVw2CwyαRτwindN=12ρaAsVw2CwnαR·Lwhere ρa is the density of air, Af and As are the frontal and lateral projected areas, Cwx(αR), Cwy(αR), and Cwn(αR) are the empirical force coefficients, αR is the angle between the wind and the heading of the vessel, L is the length of the vessel, Vw is the relative wind speed, and τwindX, τwindY, and τwindN are the wind forces during the surge, sway, and yaw, respectively ‎[22].


#### 2.2.2. Wave Forces

When a vehicle is sailing on the sea, the interference of wave forces is complicated. The wave forces acting on the hull are first- and second-order wave forces. The second-order wave forces, which impact the heading and path of the USV, are proportional to the square of the wave height ‎[22]. The wave forces are simplified as follows:(4)τwaveX=Kw1ss2+2λ1ωe1s+ωe12w1+d1τwaveY=Kw2ss2+2λ2ωe2s+ωe22w2+d2τwaveN=Kw3ss2+2λ3ωe3s+ωe32w3+d3where wi(i=1,2,3) are Gaussian white noise processes, and τwaveX, τwaveY, and τwaveN are the wave forces during the surge, sway, and yaw, respectively. The amplitudes of τwaveX, τwaveY, and τwaveN are adjusted by choosing the constants Kwi(i=1,2,3), while the spectra are parameterized in terms of the pairs λi and ωei(i=1,2,3). The wave drift forces di(i=1,2,3) are usually modeled as slowly varying bias terms:(5)d˙1=w4d˙2=w5d˙3=w6where wi(i=4,5,6) are Gaussian white noise processes ‎[22].


#### 2.2.3. Ocean Current Forces

The ocean currents cause vessels sailing on the sea to change their positions and postures. The ocean current forces are given as follows:(6)τcurrentX=12ρAfVc2CXβτcurrentY=12ρAsVc2CYβτcurrentN=12ρAsVc2CNβ·Lwhere ρ is the density of the seawater, Af and As are the frontal and lateral projected areas below the waterline, respectively, CX, CY, and CN are the empirical force coefficients, Vc is the relative current speed, β is the angle between the ocean current and the heading of the vessel, L is the length of the vessel, and τcurrentX, τcurrentY, and τcurrentN are the ocean current forces during surge, sway, and yaw, respectively ‎[22].


### 2.3. Path Representation by Grids

The real task area is partitioned to reduce the modeling complexity. Visibility graphs ‎[23], Voronoi diagrams ‎[24], and grid maps ‎[25] are the most commonly used path planning algorithms. The grid map-based path planning algorithm is powerful in that it generates a path with the shortest computation time ‎[25]. To facilitate the calculation, the planned path is represented on grids. The area under consideration is discretized into grids. The information, such as the relative speed and direction of the wind, the amplitude and direction of the waves, the relative speed and direction of ocean current, and the position of the obstacles, is discretized in each grid. Stationary obstacles are encoded in a binary format on the grids. We assigned weights of 1 to all obstacle grids and weights of 0 to all free neighbor grids of them.


### 2.4. Objectives of USV Global Path Planning

Since USV global path planning is a multiobjective optimization problem, we should analyze the interrelated objectives and discuss the importance of each objective based on the requirements of the mission. A cost function can be constructed as a weighted sum of the objective functions. Finally, the cost function is used to evaluate the quality of the planned path.


#### 2.4.1. Path Length

Since the task area is modeled by grids, the planned path is represented on a rectangular grid. The path passes the centers of the grids. Thus, the distanceLi,i+1 between two adjacent waypoints pi=(xi,yi) and pi+1=(xi+1,yi+1) is equal to the Euclidean distance between the centers of the grids as follows:(7)Li,i+1=1,xi=xi+1oryi=yi+12,otherwiseThe positions ofpi and pi+1 are shown in Figure 1. If pi and pi+1 are adjacent in the horizontal or vertical direction, Li,i+1=1. If pi and pi+1 are adjacent in the diagonal direction, Li,i+1=2.Figure 1
The positions ofpi and pi+1.Therefore, the total length of the pathL is the sum of the distances between the adjacent waypoints:(8)L=∑i=1mLi,i+1where m is the number of path segments.


#### 2.4.2. Energy Consumption

In this paper, the energy consumption of the USV while sailing is derived from the propulsion system. Thus,E is the sum of the energy consumption of each segment along the entire path:(9)E=∑i=1mEi,i+1Supposing that the USV is sailing at a constant velocity betweenpi and pi+1, the energy consumption Ei,i+1 between pi and pi+1 equals the work done by the propulsion system to overcome the environmental loads, such that(10)Ei,i+1=τenv·v→usv·twhere t is the time for the USV to sail in Li,i+1.(11)t=Li,i+1v→outwhere v→usv is the magnitude of the velocity v→usv generated by the USV propulsion system, τenv is the resultant force of the environmental loads, and v→out is the magnitude of the velocity v→out of the USV moving in the horizontal plane. Since the headings of the USV in the grid are several fixed values, as shown in Figure 1, the angular velocity r caused by the yaw motion can be ignored when solving v→out. Hence, v→out is equal to(12)v→out=u2+v2where u and v are obtained by (1).It is known from (10) and (11) that the energy consumption is proportional to v→usv and 1/v→out, when the thrust τthr generated by the propulsion system is a fixed value. To reduce the energy consumption, it is necessary to adjust the USV’s heading to take advantage of the environmental loads to increase v→out.


#### 2.4.3. Path Smoothness

It is assumed that the current waypoint of the USV ispi=(xi,yi), the previous waypoint is pi-1=(xi-1,yi-1), and the next waypoint is pi+1=(xi+1,yi+1). Thus, the angle θi+1 of the vector pipi+1→ and the angle θi of the vector pi-1pi→ are(13)θi+1=arctan⁡yi+1-yixi+1-xi(14)θi=arctan⁡yi-yi-1xi-xi-1The differenceψi between θi+1 and θi is(15)ψi=absθi+1-θiθ i, θi+1, and ψi, are shown in Figure 2. Therefore, the cost function of the path smoothness Jsmooth is(16)Jsmooth=∑i=1Nψiwhere N is the number of differences ψi.Figure 2
The angles of the path segments.


#### 2.4.4. Path Safety

Using the safety cost of the nodes on the grids cannot accurately represent the threat impact of each path segment. First, three sampling points are selected on a path segment and the average Euclidean distance between the three sampling points and the center of the obstacle is calculated. The schematic diagram of the calculation of the path safety is shown in Figure3. lk is the length of the kth path segment between the waypoint pi and pi+1. For the kth path segment, three sampling points are taken at lk/6, lk/2, and 5lk/6, respectively. The average Euclidean distance between the three sampling points and the center of the obstacle is(17)DT,k=13dT,kjlk6+dT,kjlk2+dT,kj5lk6where dT,kj() is the Euclidean distance from the sampling point on the kth path segment to the center of the obstacle Tj.Figure 3
Schematic diagram of the calculation of the path safety.The path safety cost between waypointpi and pi+1 denoted as Jsafei,i+1 is(18)Jsafei,i+1=0,d>dsafe_max1DT,ki,i+1,dsafe_min≤d≤dsafe_max1,d<dsafe_minwhere d is the distance between the USV and the obstacle’s center, dsafe_max is the radius of the obstacle’s affected area, and dsafe_min is the radius of the no-sail zone. JT,ki,i+1 is obtained using (17).Thus, the entire path safety cost functionJsafe is(19)Jsafe=∑i=1NJsafei,i+1where N is the number of the waypoints of the planned path.


#### 2.4.5. Cost Function

In summary, the cost function of the USV global path planning was established as(20)min⁡J=w1·L+w2·E+w3·Jsmooth+w4·Jsafewhere L, E, Jsmooth, and Jsafe are obtained by (8), (9), (16), and (19), respectively. w1, w2, w3, and w4 represent the weights of the path length, energy consumption, path smoothness, and path safety in the cost function, respectively, subject to(21)L≤Lmax0≤v→out≤vmaxwhere Lmax is the maximum voyage distance of the USV and vmax is the maximum speed of the USV.


## 2.1. USV Kinetic Model

The kinetic model of a USV accounts for the forces, such as the control force and environmental loads, which cause USV motion. For the USV, the control force is mainly the thrust of each propeller. The environmental loads on the USV are generated by wind, waves, and ocean currents. The kinetic model of the USV, which was proposed previously ‎[20], is as follows:(1)Mν˙+Cνν+Dνν=τenv+τthr(2)τenv=τwind+τwave+τcurrentwhere M is the system inertia matrix, C(ν) is the Coriolis-centripetal matrix, D(ν)∈R3×3 is damping matrix. τwind, τwave, and τcurrent are wind, wave, and ocean current forces acting on the USV, respectively, and τthr is the thrust generated by the USV propulsion system. The generalized velocity ν=[u,v,r]T is obtained by (1), where the first two components (u,v) are the linear velocities of the surge and sway, and r is the angular velocity of the yaw.


## 2.2. Models of Environmental Loads

When planning a global path for USVs, it is necessary to consider the environmental effects on the vehicles. Thus, we need to analyze the impacts of wind, waves, and ocean currents on the USV. The planned area is a confined sea with some static obstacles, and the mission execution time is short. Therefore, it can be assumed that the environmental loads are basically stable in limited time and space.


### 2.2.1. Wind Forces

The wind acts directly on the superstructure of the hull. As reported previously ‎[21], the wind forces are written as follows:(3)τwindX=12ρaAfVw2CwxαRτwindY=12ρaAsVw2CwyαRτwindN=12ρaAsVw2CwnαR·Lwhere ρa is the density of air, Af and As are the frontal and lateral projected areas, Cwx(αR), Cwy(αR), and Cwn(αR) are the empirical force coefficients, αR is the angle between the wind and the heading of the vessel, L is the length of the vessel, Vw is the relative wind speed, and τwindX, τwindY, and τwindN are the wind forces during the surge, sway, and yaw, respectively ‎[22].


### 2.2.2. Wave Forces

When a vehicle is sailing on the sea, the interference of wave forces is complicated. The wave forces acting on the hull are first- and second-order wave forces. The second-order wave forces, which impact the heading and path of the USV, are proportional to the square of the wave height ‎[22]. The wave forces are simplified as follows:(4)τwaveX=Kw1ss2+2λ1ωe1s+ωe12w1+d1τwaveY=Kw2ss2+2λ2ωe2s+ωe22w2+d2τwaveN=Kw3ss2+2λ3ωe3s+ωe32w3+d3where wi(i=1,2,3) are Gaussian white noise processes, and τwaveX, τwaveY, and τwaveN are the wave forces during the surge, sway, and yaw, respectively. The amplitudes of τwaveX, τwaveY, and τwaveN are adjusted by choosing the constants Kwi(i=1,2,3), while the spectra are parameterized in terms of the pairs λi and ωei(i=1,2,3). The wave drift forces di(i=1,2,3) are usually modeled as slowly varying bias terms:(5)d˙1=w4d˙2=w5d˙3=w6where wi(i=4,5,6) are Gaussian white noise processes ‎[22].


### 2.2.3. Ocean Current Forces

The ocean currents cause vessels sailing on the sea to change their positions and postures. The ocean current forces are given as follows:(6)τcurrentX=12ρAfVc2CXβτcurrentY=12ρAsVc2CYβτcurrentN=12ρAsVc2CNβ·Lwhere ρ is the density of the seawater, Af and As are the frontal and lateral projected areas below the waterline, respectively, CX, CY, and CN are the empirical force coefficients, Vc is the relative current speed, β is the angle between the ocean current and the heading of the vessel, L is the length of the vessel, and τcurrentX, τcurrentY, and τcurrentN are the ocean current forces during surge, sway, and yaw, respectively ‎[22].


## 2.2.1. Wind Forces

The wind acts directly on the superstructure of the hull. As reported previously ‎[21], the wind forces are written as follows:(3)τwindX=12ρaAfVw2CwxαRτwindY=12ρaAsVw2CwyαRτwindN=12ρaAsVw2CwnαR·Lwhere ρa is the density of air, Af and As are the frontal and lateral projected areas, Cwx(αR), Cwy(αR), and Cwn(αR) are the empirical force coefficients, αR is the angle between the wind and the heading of the vessel, L is the length of the vessel, Vw is the relative wind speed, and τwindX, τwindY, and τwindN are the wind forces during the surge, sway, and yaw, respectively ‎[22].


## 2.2.2. Wave Forces

When a vehicle is sailing on the sea, the interference of wave forces is complicated. The wave forces acting on the hull are first- and second-order wave forces. The second-order wave forces, which impact the heading and path of the USV, are proportional to the square of the wave height ‎[22]. The wave forces are simplified as follows:(4)τwaveX=Kw1ss2+2λ1ωe1s+ωe12w1+d1τwaveY=Kw2ss2+2λ2ωe2s+ωe22w2+d2τwaveN=Kw3ss2+2λ3ωe3s+ωe32w3+d3where wi(i=1,2,3) are Gaussian white noise processes, and τwaveX, τwaveY, and τwaveN are the wave forces during the surge, sway, and yaw, respectively. The amplitudes of τwaveX, τwaveY, and τwaveN are adjusted by choosing the constants Kwi(i=1,2,3), while the spectra are parameterized in terms of the pairs λi and ωei(i=1,2,3). The wave drift forces di(i=1,2,3) are usually modeled as slowly varying bias terms:(5)d˙1=w4d˙2=w5d˙3=w6where wi(i=4,5,6) are Gaussian white noise processes ‎[22].


## 2.2.3. Ocean Current Forces

The ocean currents cause vessels sailing on the sea to change their positions and postures. The ocean current forces are given as follows:(6)τcurrentX=12ρAfVc2CXβτcurrentY=12ρAsVc2CYβτcurrentN=12ρAsVc2CNβ·Lwhere ρ is the density of the seawater, Af and As are the frontal and lateral projected areas below the waterline, respectively, CX, CY, and CN are the empirical force coefficients, Vc is the relative current speed, β is the angle between the ocean current and the heading of the vessel, L is the length of the vessel, and τcurrentX, τcurrentY, and τcurrentN are the ocean current forces during surge, sway, and yaw, respectively ‎[22].


## 2.3. Path Representation by Grids

The real task area is partitioned to reduce the modeling complexity. Visibility graphs ‎[23], Voronoi diagrams ‎[24], and grid maps ‎[25] are the most commonly used path planning algorithms. The grid map-based path planning algorithm is powerful in that it generates a path with the shortest computation time ‎[25]. To facilitate the calculation, the planned path is represented on grids. The area under consideration is discretized into grids. The information, such as the relative speed and direction of the wind, the amplitude and direction of the waves, the relative speed and direction of ocean current, and the position of the obstacles, is discretized in each grid. Stationary obstacles are encoded in a binary format on the grids. We assigned weights of 1 to all obstacle grids and weights of 0 to all free neighbor grids of them.


## 2.4. Objectives of USV Global Path Planning

Since USV global path planning is a multiobjective optimization problem, we should analyze the interrelated objectives and discuss the importance of each objective based on the requirements of the mission. A cost function can be constructed as a weighted sum of the objective functions. Finally, the cost function is used to evaluate the quality of the planned path.


### 2.4.1. Path Length

Since the task area is modeled by grids, the planned path is represented on a rectangular grid. The path passes the centers of the grids. Thus, the distanceLi,i+1 between two adjacent waypoints pi=(xi,yi) and pi+1=(xi+1,yi+1) is equal to the Euclidean distance between the centers of the grids as follows:(7)Li,i+1=1,xi=xi+1oryi=yi+12,otherwiseThe positions ofpi and pi+1 are shown in Figure 1. If pi and pi+1 are adjacent in the horizontal or vertical direction, Li,i+1=1. If pi and pi+1 are adjacent in the diagonal direction, Li,i+1=2.Figure 1
The positions ofpi and pi+1.Therefore, the total length of the pathL is the sum of the distances between the adjacent waypoints:(8)L=∑i=1mLi,i+1where m is the number of path segments.


### 2.4.2. Energy Consumption

In this paper, the energy consumption of the USV while sailing is derived from the propulsion system. Thus,E is the sum of the energy consumption of each segment along the entire path:(9)E=∑i=1mEi,i+1Supposing that the USV is sailing at a constant velocity betweenpi and pi+1, the energy consumption Ei,i+1 between pi and pi+1 equals the work done by the propulsion system to overcome the environmental loads, such that(10)Ei,i+1=τenv·v→usv·twhere t is the time for the USV to sail in Li,i+1.(11)t=Li,i+1v→outwhere v→usv is the magnitude of the velocity v→usv generated by the USV propulsion system, τenv is the resultant force of the environmental loads, and v→out is the magnitude of the velocity v→out of the USV moving in the horizontal plane. Since the headings of the USV in the grid are several fixed values, as shown in Figure 1, the angular velocity r caused by the yaw motion can be ignored when solving v→out. Hence, v→out is equal to(12)v→out=u2+v2where u and v are obtained by (1).It is known from (10) and (11) that the energy consumption is proportional to v→usv and 1/v→out, when the thrust τthr generated by the propulsion system is a fixed value. To reduce the energy consumption, it is necessary to adjust the USV’s heading to take advantage of the environmental loads to increase v→out.


### 2.4.3. Path Smoothness

It is assumed that the current waypoint of the USV ispi=(xi,yi), the previous waypoint is pi-1=(xi-1,yi-1), and the next waypoint is pi+1=(xi+1,yi+1). Thus, the angle θi+1 of the vector pipi+1→ and the angle θi of the vector pi-1pi→ are(13)θi+1=arctan⁡yi+1-yixi+1-xi(14)θi=arctan⁡yi-yi-1xi-xi-1The differenceψi between θi+1 and θi is(15)ψi=absθi+1-θiθ i, θi+1, and ψi, are shown in Figure 2. Therefore, the cost function of the path smoothness Jsmooth is(16)Jsmooth=∑i=1Nψiwhere N is the number of differences ψi.Figure 2
The angles of the path segments.


### 2.4.4. Path Safety

Using the safety cost of the nodes on the grids cannot accurately represent the threat impact of each path segment. First, three sampling points are selected on a path segment and the average Euclidean distance between the three sampling points and the center of the obstacle is calculated. The schematic diagram of the calculation of the path safety is shown in Figure3. lk is the length of the kth path segment between the waypoint pi and pi+1. For the kth path segment, three sampling points are taken at lk/6, lk/2, and 5lk/6, respectively. The average Euclidean distance between the three sampling points and the center of the obstacle is(17)DT,k=13dT,kjlk6+dT,kjlk2+dT,kj5lk6where dT,kj() is the Euclidean distance from the sampling point on the kth path segment to the center of the obstacle Tj.Figure 3
Schematic diagram of the calculation of the path safety.The path safety cost between waypointpi and pi+1 denoted as Jsafei,i+1 is(18)Jsafei,i+1=0,d>dsafe_max1DT,ki,i+1,dsafe_min≤d≤dsafe_max1,d<dsafe_minwhere d is the distance between the USV and the obstacle’s center, dsafe_max is the radius of the obstacle’s affected area, and dsafe_min is the radius of the no-sail zone. JT,ki,i+1 is obtained using (17).Thus, the entire path safety cost functionJsafe is(19)Jsafe=∑i=1NJsafei,i+1where N is the number of the waypoints of the planned path.


### 2.4.5. Cost Function

In summary, the cost function of the USV global path planning was established as(20)min⁡J=w1·L+w2·E+w3·Jsmooth+w4·Jsafewhere L, E, Jsmooth, and Jsafe are obtained by (8), (9), (16), and (19), respectively. w1, w2, w3, and w4 represent the weights of the path length, energy consumption, path smoothness, and path safety in the cost function, respectively, subject to(21)L≤Lmax0≤v→out≤vmaxwhere Lmax is the maximum voyage distance of the USV and vmax is the maximum speed of the USV.


## 2.4.1. Path Length

Since the task area is modeled by grids, the planned path is represented on a rectangular grid. The path passes the centers of the grids. Thus, the distanceLi,i+1 between two adjacent waypoints pi=(xi,yi) and pi+1=(xi+1,yi+1) is equal to the Euclidean distance between the centers of the grids as follows:(7)Li,i+1=1,xi=xi+1oryi=yi+12,otherwiseThe positions ofpi and pi+1 are shown in Figure 1. If pi and pi+1 are adjacent in the horizontal or vertical direction, Li,i+1=1. If pi and pi+1 are adjacent in the diagonal direction, Li,i+1=2.Figure 1
The positions ofpi and pi+1.Therefore, the total length of the pathL is the sum of the distances between the adjacent waypoints:(8)L=∑i=1mLi,i+1where m is the number of path segments.


## 2.4.2. Energy Consumption

In this paper, the energy consumption of the USV while sailing is derived from the propulsion system. Thus,E is the sum of the energy consumption of each segment along the entire path:(9)E=∑i=1mEi,i+1Supposing that the USV is sailing at a constant velocity betweenpi and pi+1, the energy consumption Ei,i+1 between pi and pi+1 equals the work done by the propulsion system to overcome the environmental loads, such that(10)Ei,i+1=τenv·v→usv·twhere t is the time for the USV to sail in Li,i+1.(11)t=Li,i+1v→outwhere v→usv is the magnitude of the velocity v→usv generated by the USV propulsion system, τenv is the resultant force of the environmental loads, and v→out is the magnitude of the velocity v→out of the USV moving in the horizontal plane. Since the headings of the USV in the grid are several fixed values, as shown in Figure 1, the angular velocity r caused by the yaw motion can be ignored when solving v→out. Hence, v→out is equal to(12)v→out=u2+v2where u and v are obtained by (1).It is known from (10) and (11) that the energy consumption is proportional to v→usv and 1/v→out, when the thrust τthr generated by the propulsion system is a fixed value. To reduce the energy consumption, it is necessary to adjust the USV’s heading to take advantage of the environmental loads to increase v→out.


## 2.4.3. Path Smoothness

It is assumed that the current waypoint of the USV ispi=(xi,yi), the previous waypoint is pi-1=(xi-1,yi-1), and the next waypoint is pi+1=(xi+1,yi+1). Thus, the angle θi+1 of the vector pipi+1→ and the angle θi of the vector pi-1pi→ are(13)θi+1=arctan⁡yi+1-yixi+1-xi(14)θi=arctan⁡yi-yi-1xi-xi-1The differenceψi between θi+1 and θi is(15)ψi=absθi+1-θiθ i, θi+1, and ψi, are shown in Figure 2. Therefore, the cost function of the path smoothness Jsmooth is(16)Jsmooth=∑i=1Nψiwhere N is the number of differences ψi.Figure 2
The angles of the path segments.


## 2.4.4. Path Safety

Using the safety cost of the nodes on the grids cannot accurately represent the threat impact of each path segment. First, three sampling points are selected on a path segment and the average Euclidean distance between the three sampling points and the center of the obstacle is calculated. The schematic diagram of the calculation of the path safety is shown in Figure3. lk is the length of the kth path segment between the waypoint pi and pi+1. For the kth path segment, three sampling points are taken at lk/6, lk/2, and 5lk/6, respectively. The average Euclidean distance between the three sampling points and the center of the obstacle is(17)DT,k=13dT,kjlk6+dT,kjlk2+dT,kj5lk6where dT,kj() is the Euclidean distance from the sampling point on the kth path segment to the center of the obstacle Tj.Figure 3
Schematic diagram of the calculation of the path safety.The path safety cost between waypointpi and pi+1 denoted as Jsafei,i+1 is(18)Jsafei,i+1=0,d>dsafe_max1DT,ki,i+1,dsafe_min≤d≤dsafe_max1,d<dsafe_minwhere d is the distance between the USV and the obstacle’s center, dsafe_max is the radius of the obstacle’s affected area, and dsafe_min is the radius of the no-sail zone. JT,ki,i+1 is obtained using (17).Thus, the entire path safety cost functionJsafe is(19)Jsafe=∑i=1NJsafei,i+1where N is the number of the waypoints of the planned path.


## 2.4.5. Cost Function

In summary, the cost function of the USV global path planning was established as(20)min⁡J=w1·L+w2·E+w3·Jsmooth+w4·Jsafewhere L, E, Jsmooth, and Jsafe are obtained by (8), (9), (16), and (19), respectively. w1, w2, w3, and w4 represent the weights of the path length, energy consumption, path smoothness, and path safety in the cost function, respectively, subject to(21)L≤Lmax0≤v→out≤vmaxwhere Lmax is the maximum voyage distance of the USV and vmax is the maximum speed of the USV.


## 3. Optimization Algorithm

In this section, we will introduce the optimization algorithm for the USV global path planning—the IQACA. The IQACA is a new optimization algorithm that combines quantum-inspired computing with ant colony optimization algorithm. We will introduce quantum code and a quantum rotation gate from quantum-inspired computing. Some rules based on the ant colony optimization algorithm are presented.


### 3.1. Quantum Code

The quantum bit (Q-bit) is the basic unit in quantum computing. A Q-bit is a system that has two possible states0 and 1. The state of a Q-bit φ is expressed as(22)φ=α0+β1where α and β are the probability amplitudes, which satisfy α2+β2=1. α2 and β2 are the probabilities in states 0 and 1, respectively. Thus, the state of the Q-bit φ is an uncertain superposition state between 0 and 1. When the number of Q-bits of an individual Xi is n, Xi is expressed as(23)Xi=αi1αi2…αinβi1βi2…βinwhere Xic=(αi1,⋯,αin) and Xis=(βi1,⋯,βin) are the two sets of solutions for individual Xi. Therefore, after quantum coding, every individual has two sets of solutions and the search space is doubled.In the IQACA, the quantum pheromone is obtained by encoding the pheromone left by the ants on the path in the ACA by the Q-bits. The transfer direction of the ants is selected by the quantum pheromone concentration on the path. Thus, the quantum pheromone concentration valueτijt of the ith ant on the jth point in the tth iteration is expressed as(24)τijt=αijtβijt


### 3.2. Adaptive Quantum Rotation Gate

In the quantum optimization algorithm, a quantum rotation gate is used to update the Q-bits. The update rule of a Q-bit is as follows:(25)αijt+1βijt+1=Uθtαijtβijtwhere [αijt,βijt]T represents the probability amplitude of the Q-bits in the tth iteration. U(θt) is the quantum rotation gate in the tth iteration(26)Uθt=cos⁡θt-sin⁡θtsin⁡θtcos⁡θtwhere θt is the rotation angle in the tth iteration. In a previous paper ‎[13], the rotation angle was obtained by looking it up in a table. In another paper ‎[26], the local and global updates of the pheromone concentration increments in the ACA were added to the rotation angle step function. In the IQACA, an adaptive adjustment strategy for the rotation angle is obtained by comparing the current solution and the global optimal solution currently being searched. Thus, the rotation angle θt in the tth iteration is(27)θt=-sgn⁡Ai·Δθiwhere -sgn⁡Ai is the direction of the rotation angle and Δθi is the size of the rotation angle. Ai is(28)Ai=α0α1β0β1where α0 and β0 are the probability amplitudes of the quantum pheromone corresponding to the global optimal solution currently searched and α1 and β1 are the probability amplitudes of the quantum pheromone corresponding to the current solution. Δθi is(29)Δθi=Je-Je-JkNmax·twhere Jk is the cost value of ant k in the current solution, Je is the cost value of the global optimal solution currently searched, and Nmax is the maximum number of iterations.


### 3.3. Transfer Rule and Transition Probability

The ant colony optimization algorithm is a bionic intelligent algorithm inspired by the foraging behavior of ant colonies ‎[27]. During the foraging, ants produce a substance called a pheromone. The concentration of the pheromone, which is related to the path length, will determine the movement of other ants. If the path is shorter, the concentration of the pheromone left on the path is larger.To achieve multiobjective path planning, multiple pieces of heuristic information are used to determine the ant’s transfer rules and transition probabilities. The transfer rule of antk from point i to point j is(30)s=arg⁡maxs∈Sτijktαηijktβεijtγq≤q0s~q>q0where q is a random number in the range [0,1]. q0 is a constant within [0,1]. S is the set of points that ant k may reach by point i. s~ is the target waypoint selected by the following equation:(31)pijkt=τijktαηijktβμijktγ∑s∈allowediτijktαηijktβμijktγwhere τijk(t) is the pheromone on the path from point i to point j in the tth iteration and α(α>0) is the pheromone index. ηijk(t) is the multiple inspiration information on the path from point i to point j in the tth iteration, β(β>0) is the index of multiple inspiration information, μijk(t) is the quantum information strength on the path from point i to point j in the tth iteration, which is expressed as μijk(t)=1/αijk(t)2, and γ(γ>0) is the index of the quantum information strength.The multiple pieces of heuristic information include the path length heuristic informationϕijk(t), energy consumption heuristic information εijk(t), path smoothness heuristic information φijk(t), and path safety heuristic information ζijk(t).(32)ηijktβ=ϕijktaεijktbφijktcζijktd(33)ϕijt=1Lijεijt=1Eijφijt=1Jsmoothijζijt=1Jsafeijwhere Lij, Eij, Jsmoothij, and Jsafeij are obtained by (8), (9), (16), and (19), respectively. a, b, c, and d are the indices of the path length heuristic information, energy consumption heuristic information, path smoothness heuristic information, and path safety heuristic information, respectively.


### 3.4. Update Rules of Pheromone

After every ant completes a one-transfer, the pheromone on the path it passes is locally updated to avoid falling into a local optimum. When the current point of the ant ispi and the next point is pj, the pheromone local updating rule is(34)τpj=1-ρ1·τpi+ρ1·Δτijwhere τ(pi) is the pheromone of the current point, τ(pj) is the pheromone of the next point, ρ1(0<ρ1<1) is the pheromone local updating coefficient, and Δτij is the pheromone that every ant leaves on the path from pi to pj in this iteration, expressed as follows:(35)Δτij=∑k=1nΔτijk(36)Δτijk=QJk,thepathsegmentofthekthant0,elsewhere Q is a constant and Jk is the cost value of the kth ant’s path.After all the ants complete an iteration, the pheromone is globally updated to increase the pheromone concentration on the optimized path. The rules are as follows:(37)τpj=1-ρ2·τpj+ρ2·Δτbest,pj=s~1-ρ2·τpj,pj≠s~(38)Δτbest=QJe,ifi,jbelongstotheoptimalpathinthiscycle0,elsewhere Q is a constant, Je is the cost value of the optimal path in this iteration, ρ2(0<ρ2<1) is the pheromone global updating coefficient, and s~ is the global optimal solution currently being searched.


### 3.5. Global Path Planning Algorithm Based on IQACA

The flowchart of global path planning algorithm based on the IQACA is shown in Figure4.Figure 4
The flowchart of the global path planning algorithm based on the IQACA.The main steps are as follows:Step 1 (initialize the parameters).
The number of the ants in the colony isn. The maximum number of iterations is Nmax. The initial quantum pheromone concentration value of the ith ant on the jth waypoint is expressed as τijt=[αijt,βijt]T=[1/2,1/2]T, where t=0;Step 2.
The ants are placed at the starting point. The transfer rule of the ants is determined by (30), and the target waypoint is selected by (31).Step 3.
The phases of the Q-bits are updated by (25).Step 4.
The pheromone is locally updated by (34).Step 5.
After all the ants have passed by all the points in an iteration, the pheromone is globally updated by (37).Step 6.
The candidate solution selected by the ants is output and the path cost is calculated.Step 7.
If the iterationt>Nmax, the algorithm moves to Step 8; otherwise, it returns to Step 2.Step 8.
The waypoints of the optimized solution and the cost value of the path are output. The global optimized path is obtained by the waypoints of the optimized solution.Step 9.
The algorithm ends.


## 3.1. Quantum Code

The quantum bit (Q-bit) is the basic unit in quantum computing. A Q-bit is a system that has two possible states0 and 1. The state of a Q-bit φ is expressed as(22)φ=α0+β1where α and β are the probability amplitudes, which satisfy α2+β2=1. α2 and β2 are the probabilities in states 0 and 1, respectively. Thus, the state of the Q-bit φ is an uncertain superposition state between 0 and 1. When the number of Q-bits of an individual Xi is n, Xi is expressed as(23)Xi=αi1αi2…αinβi1βi2…βinwhere Xic=(αi1,⋯,αin) and Xis=(βi1,⋯,βin) are the two sets of solutions for individual Xi. Therefore, after quantum coding, every individual has two sets of solutions and the search space is doubled.In the IQACA, the quantum pheromone is obtained by encoding the pheromone left by the ants on the path in the ACA by the Q-bits. The transfer direction of the ants is selected by the quantum pheromone concentration on the path. Thus, the quantum pheromone concentration valueτijt of the ith ant on the jth point in the tth iteration is expressed as(24)τijt=αijtβijt


## 3.2. Adaptive Quantum Rotation Gate

In the quantum optimization algorithm, a quantum rotation gate is used to update the Q-bits. The update rule of a Q-bit is as follows:(25)αijt+1βijt+1=Uθtαijtβijtwhere [αijt,βijt]T represents the probability amplitude of the Q-bits in the tth iteration. U(θt) is the quantum rotation gate in the tth iteration(26)Uθt=cos⁡θt-sin⁡θtsin⁡θtcos⁡θtwhere θt is the rotation angle in the tth iteration. In a previous paper ‎[13], the rotation angle was obtained by looking it up in a table. In another paper ‎[26], the local and global updates of the pheromone concentration increments in the ACA were added to the rotation angle step function. In the IQACA, an adaptive adjustment strategy for the rotation angle is obtained by comparing the current solution and the global optimal solution currently being searched. Thus, the rotation angle θt in the tth iteration is(27)θt=-sgn⁡Ai·Δθiwhere -sgn⁡Ai is the direction of the rotation angle and Δθi is the size of the rotation angle. Ai is(28)Ai=α0α1β0β1where α0 and β0 are the probability amplitudes of the quantum pheromone corresponding to the global optimal solution currently searched and α1 and β1 are the probability amplitudes of the quantum pheromone corresponding to the current solution. Δθi is(29)Δθi=Je-Je-JkNmax·twhere Jk is the cost value of ant k in the current solution, Je is the cost value of the global optimal solution currently searched, and Nmax is the maximum number of iterations.


## 3.3. Transfer Rule and Transition Probability

The ant colony optimization algorithm is a bionic intelligent algorithm inspired by the foraging behavior of ant colonies ‎[27]. During the foraging, ants produce a substance called a pheromone. The concentration of the pheromone, which is related to the path length, will determine the movement of other ants. If the path is shorter, the concentration of the pheromone left on the path is larger.To achieve multiobjective path planning, multiple pieces of heuristic information are used to determine the ant’s transfer rules and transition probabilities. The transfer rule of antk from point i to point j is(30)s=arg⁡maxs∈Sτijktαηijktβεijtγq≤q0s~q>q0where q is a random number in the range [0,1]. q0 is a constant within [0,1]. S is the set of points that ant k may reach by point i. s~ is the target waypoint selected by the following equation:(31)pijkt=τijktαηijktβμijktγ∑s∈allowediτijktαηijktβμijktγwhere τijk(t) is the pheromone on the path from point i to point j in the tth iteration and α(α>0) is the pheromone index. ηijk(t) is the multiple inspiration information on the path from point i to point j in the tth iteration, β(β>0) is the index of multiple inspiration information, μijk(t) is the quantum information strength on the path from point i to point j in the tth iteration, which is expressed as μijk(t)=1/αijk(t)2, and γ(γ>0) is the index of the quantum information strength.The multiple pieces of heuristic information include the path length heuristic informationϕijk(t), energy consumption heuristic information εijk(t), path smoothness heuristic information φijk(t), and path safety heuristic information ζijk(t).(32)ηijktβ=ϕijktaεijktbφijktcζijktd(33)ϕijt=1Lijεijt=1Eijφijt=1Jsmoothijζijt=1Jsafeijwhere Lij, Eij, Jsmoothij, and Jsafeij are obtained by (8), (9), (16), and (19), respectively. a, b, c, and d are the indices of the path length heuristic information, energy consumption heuristic information, path smoothness heuristic information, and path safety heuristic information, respectively.


## 3.4. Update Rules of Pheromone

After every ant completes a one-transfer, the pheromone on the path it passes is locally updated to avoid falling into a local optimum. When the current point of the ant ispi and the next point is pj, the pheromone local updating rule is(34)τpj=1-ρ1·τpi+ρ1·Δτijwhere τ(pi) is the pheromone of the current point, τ(pj) is the pheromone of the next point, ρ1(0<ρ1<1) is the pheromone local updating coefficient, and Δτij is the pheromone that every ant leaves on the path from pi to pj in this iteration, expressed as follows:(35)Δτij=∑k=1nΔτijk(36)Δτijk=QJk,thepathsegmentofthekthant0,elsewhere Q is a constant and Jk is the cost value of the kth ant’s path.After all the ants complete an iteration, the pheromone is globally updated to increase the pheromone concentration on the optimized path. The rules are as follows:(37)τpj=1-ρ2·τpj+ρ2·Δτbest,pj=s~1-ρ2·τpj,pj≠s~(38)Δτbest=QJe,ifi,jbelongstotheoptimalpathinthiscycle0,elsewhere Q is a constant, Je is the cost value of the optimal path in this iteration, ρ2(0<ρ2<1) is the pheromone global updating coefficient, and s~ is the global optimal solution currently being searched.


## 3.5. Global Path Planning Algorithm Based on IQACA

The flowchart of global path planning algorithm based on the IQACA is shown in Figure4.Figure 4
The flowchart of the global path planning algorithm based on the IQACA.The main steps are as follows:Step 1 (initialize the parameters).
The number of the ants in the colony isn. The maximum number of iterations is Nmax. The initial quantum pheromone concentration value of the ith ant on the jth waypoint is expressed as τijt=[αijt,βijt]T=[1/2,1/2]T, where t=0;Step 2.
The ants are placed at the starting point. The transfer rule of the ants is determined by (30), and the target waypoint is selected by (31).Step 3.
The phases of the Q-bits are updated by (25).Step 4.
The pheromone is locally updated by (34).Step 5.
After all the ants have passed by all the points in an iteration, the pheromone is globally updated by (37).Step 6.
The candidate solution selected by the ants is output and the path cost is calculated.Step 7.
If the iterationt>Nmax, the algorithm moves to Step 8; otherwise, it returns to Step 2.Step 8.
The waypoints of the optimized solution and the cost value of the path are output. The global optimized path is obtained by the waypoints of the optimized solution.Step 9.
The algorithm ends.


## 4. Simulation Studies

In this section, the effectiveness and efficiency of the IQACA are validated. The section consists of two parts. The first subsection compares the performance of the ACA, QACA, and IQACA with the Traveling Salesman Problem (TSP). The second subsection deals with the USV global path planning based on the IQACA. To validate the proposed algorithm, simulations were conducted.


### 4.1. Performance Evaluation of IQACA

To validate the effectiveness of the IQACA presented in this paper, we compared the algorithm performance between ACA, QACA, and IQACA with the TSP. In this paper, RAND100 was selected from the TSPLIB standard library for the simulations. The maximum number of iterationsNmax=200, the number of the ants n=100, α=3, β=1, γ=2, and ρ1=ρ2=0.8. The obtained values are shown in Table 1, and the iterations are shown in Figure 5.Table 1
The obtained values from the simulations.
Algorithm Length Iteration ACA 8447 167 QACA 8062 142 IQACA 7891 126Figure 5
Convergence of the IQACA, QACA, and ACA.The known optimal value of RAND100 is 7891. From Table1, it was concluded that the path length of RAND100 obtained by the IQACA was 2.12% lower than that obtained by the QACA and 6.58% lower than that obtained by the ACA. The number of iterations required for the IQACA to converge to the minimum was 11.27% lower than the QACA and 24.55% lower than the ACA. The results show that the IQACA was superior to the QACA and the ACA in both the path length and iteration number. Since the algorithm uses the pheromone local and global updates, and the phases of the Q-bits are updated by the adaptive quantum rotation gate, the IQACA can avoid the local optimal solution. Because the pheromone is encoded by Q-bits, the search space is doubled, and the convergence speed is faster. Thus, the IQACA is an effective and efficient algorithm.


### 4.2. USV Global Path Planning Based on IQACA

In this subsection, we will show a simulation of USV global path planning based on the proposed algorithm. Obstacles are black-colored. The coordinates of the starting point are(0.5,24.5), and the coordinates of the destination point are (29.5,0.5). The length of the side of a grid is 1km. We assumed that the wind, waves, and ocean current act on the vessel from the same direction, since in most cases the ocean current is the most significant environmental disturbance on the vessel. In this simulation, the direction of the disturbances is assumed as 240° in the Northeast coordinate system. The relative wind speed was 7.5m/s. The wave height was 2.5m. The relative ocean current was 2.0m/s. The thrust of the USV propulsion system was 500kN. The coefficients of the USV in this simulation are shown in Table 2. The maximum number of iterations was 500, the number of the ants was 100, α=3, γ=2, and ρ1=ρ2=0.8. The safety boundary of the obstacle was represented by a dotted red circle whose radius was 1.3 times the distance between the vertex and the center of the obstacle. The optimal path is represented by a solid blue line. The objectives of the global path planning were determined based on their respective weights. We considered four scenarios in this simulation. In the first scenario, we highly weighted the path length, so that the weights [w1,w2,w3,w4] in the cost function (20) were [0.5,0.1,0.3,0.1]. In the other three scenario, we weighted the energy consumption, path smoothness, and path safety, respectively, such that the weights [w1,w2,w3,w4] were [0.1,0.5,0.1,0.3], [0.3,0.1,0.5,0.1], and [0.1,0.3,0.1,0.5], respectively.Table 2
Coefficients of USV.
Property Value Length overall (m) 80.8 Breadth (m) 18.2 Draught (m) 5.0 Area of frontal projection above the waterline (m2) 330.9 Area of lateral projection above the waterline (m2) 874.8 Area of frontal projection below the waterline (m2) 91.0 Area of lateral projection below the waterline (m2) 323.4Figures6–9 show the planned paths in the four scenarios. The path length, energy consumption, path smoothness, and path safety index for each planned path are listed in Table 3. The results indicate that the proposed algorithm can plan feasible paths for the USV considering different objectives simultaneously. Moreover, by adjusting the weights of different objectives, the proposed algorithm can generate paths for different purposes.Table 3
Path data for different objectives.
Objective Path length(km) Energy consumption(1×107kJ) Path smoothness(rad) Path safety index Path length 41.87 12.38 18.85 26.89 Energy consumption 48.31 10.81 28.27 27.12 Path smoothness 43.63 11.95 14.14 15.51 Path safety 46.80 11.96 28.27 7.61Figure 6
Planned path focusing on the path length.Figure 7
Planned path focusing on the energy consumption.Figure 8
Planned path focusing on the path smoothness.Figure 9
Planned path focusing on the path safety.


## 4.1. Performance Evaluation of IQACA

To validate the effectiveness of the IQACA presented in this paper, we compared the algorithm performance between ACA, QACA, and IQACA with the TSP. In this paper, RAND100 was selected from the TSPLIB standard library for the simulations. The maximum number of iterationsNmax=200, the number of the ants n=100, α=3, β=1, γ=2, and ρ1=ρ2=0.8. The obtained values are shown in Table 1, and the iterations are shown in Figure 5.Table 1
The obtained values from the simulations.
Algorithm Length Iteration ACA 8447 167 QACA 8062 142 IQACA 7891 126Figure 5
Convergence of the IQACA, QACA, and ACA.The known optimal value of RAND100 is 7891. From Table1, it was concluded that the path length of RAND100 obtained by the IQACA was 2.12% lower than that obtained by the QACA and 6.58% lower than that obtained by the ACA. The number of iterations required for the IQACA to converge to the minimum was 11.27% lower than the QACA and 24.55% lower than the ACA. The results show that the IQACA was superior to the QACA and the ACA in both the path length and iteration number. Since the algorithm uses the pheromone local and global updates, and the phases of the Q-bits are updated by the adaptive quantum rotation gate, the IQACA can avoid the local optimal solution. Because the pheromone is encoded by Q-bits, the search space is doubled, and the convergence speed is faster. Thus, the IQACA is an effective and efficient algorithm.


## 4.2. USV Global Path Planning Based on IQACA

In this subsection, we will show a simulation of USV global path planning based on the proposed algorithm. Obstacles are black-colored. The coordinates of the starting point are(0.5,24.5), and the coordinates of the destination point are (29.5,0.5). The length of the side of a grid is 1km. We assumed that the wind, waves, and ocean current act on the vessel from the same direction, since in most cases the ocean current is the most significant environmental disturbance on the vessel. In this simulation, the direction of the disturbances is assumed as 240° in the Northeast coordinate system. The relative wind speed was 7.5m/s. The wave height was 2.5m. The relative ocean current was 2.0m/s. The thrust of the USV propulsion system was 500kN. The coefficients of the USV in this simulation are shown in Table 2. The maximum number of iterations was 500, the number of the ants was 100, α=3, γ=2, and ρ1=ρ2=0.8. The safety boundary of the obstacle was represented by a dotted red circle whose radius was 1.3 times the distance between the vertex and the center of the obstacle. The optimal path is represented by a solid blue line. The objectives of the global path planning were determined based on their respective weights. We considered four scenarios in this simulation. In the first scenario, we highly weighted the path length, so that the weights [w1,w2,w3,w4] in the cost function (20) were [0.5,0.1,0.3,0.1]. In the other three scenario, we weighted the energy consumption, path smoothness, and path safety, respectively, such that the weights [w1,w2,w3,w4] were [0.1,0.5,0.1,0.3], [0.3,0.1,0.5,0.1], and [0.1,0.3,0.1,0.5], respectively.Table 2
Coefficients of USV.
Property Value Length overall (m) 80.8 Breadth (m) 18.2 Draught (m) 5.0 Area of frontal projection above the waterline (m2) 330.9 Area of lateral projection above the waterline (m2) 874.8 Area of frontal projection below the waterline (m2) 91.0 Area of lateral projection below the waterline (m2) 323.4Figures6–9 show the planned paths in the four scenarios. The path length, energy consumption, path smoothness, and path safety index for each planned path are listed in Table 3. The results indicate that the proposed algorithm can plan feasible paths for the USV considering different objectives simultaneously. Moreover, by adjusting the weights of different objectives, the proposed algorithm can generate paths for different purposes.Table 3
Path data for different objectives.
Objective Path length(km) Energy consumption(1×107kJ) Path smoothness(rad) Path safety index Path length 41.87 12.38 18.85 26.89 Energy consumption 48.31 10.81 28.27 27.12 Path smoothness 43.63 11.95 14.14 15.51 Path safety 46.80 11.96 28.27 7.61Figure 6
Planned path focusing on the path length.Figure 7
Planned path focusing on the energy consumption.Figure 8
Planned path focusing on the path smoothness.Figure 9
Planned path focusing on the path safety.


## 5. Conclusion

This paper proposed a global path planning algorithm for the USV based on an improved quantum ant colony algorithm (IQACA).The IQACA is an optimization algorithm that combines quantum computing with the ACA. In IQACA, using Q-bits to encode the pheromone of the ants, the search space is doubled when the number of the ants is the same. The simulation results show that the proposed algorithm’s obtained minimum was 2.1–6.5% lower than those of the quantum ant colony algorithm (QACA) and ant colony algorithm (ACA), and the number of iterations required to converge to the minimum was 11.2–24.5% lower than those of the QACA and ACA. Based on the model of the kinetics of the USV and the marine environment, we defined the objectives of the path planning: the path length, energy consumption, path smoothness, and path safety. The simulation results showed that the proposed algorithm can consider several optimization objectives and generate paths satisfying these requirements.In the future, the following studies should be conducted in depth. First, the correlation between the multiple objectives should be calculated to determine the weight of each objective in the cost function to meet the actual mission requirements. Moreover, the kinetic and kinematic constraints of the USV should be added to the cost function. Finally, more practical environmental loads should be applied to calculate their effects on the path energy consumption of the USV.



---

*Source: 2902170-2019-04-24.xml*