# Local Path Planning for Unmanned Surface Vehicle Collision Avoidance Based on Modified Quantum Particle Swarm Optimization

Authors: Guoqing Xia; Zhiwei Han; Bo Zhao; Xinwei Wang
Journal: Complexity
 (2020)

---

## Abstract

An unmanned surface vehicle (USV) plans its global path before the mission starts. When dynamic obstacles appear during sailing, the planned global path must be adjusted locally to avoid collision. This study proposes a local path planning algorithm based on the velocity obstacle (VO) method and modified quantum particle swarm optimization (MQPSO) for USV collision avoidance. The collision avoidance model based on VO not only considers the velocity and course of the USV but also handles the variable velocity and course of an obstacle. According to the collision avoidance model, the USV needs to adjust its velocity and course simultaneously to avoid collision. Due to the kinematic constraints of the USV, the velocity window and course window of the USV are determined by the dynamic window approach (DWA). In summary, local path planning is transformed into a multiobjective optimization problem with multiple constraints in a continuous search space. The optimization problem is to obtain the USV’s optimal velocity variation and course variation to avoid collision and minimize its energy consumption under the rules of the International Regulations for Preventing Collisions at Sea (COLREGs) and the kinematic constraints of the USV. Since USV local path planning is completed in a short time, it is essential that the optimization algorithm can quickly obtain the optimal value. MQPSO is primarily proposed to meet that requirement. In MQPSO, the efficiency of quantum encoding in quantum computing and the optimization ability of representing the motion states of the particles with wave functions to cover the whole feasible solution space are combined. Simulation results show that the proposed algorithm can obtain the optimal values of the benchmark functions and effectively plan a collision-free path for a USV.


---

## Body

## 1. Introduction

An unmanned surface vehicle (USV) is an autonomous marine vehicle that has emerged as a viable tool to perform tasks that are dangerous or unsuitable for manned vessels. Planning the path of a USV is important as regards safety and efficiency [1]. USV path planning can be divided into two stages: global path planning based on prior environmental information and local path planning based on sensor information [2]. Since unexpected conditions can occur during USV sailing, such as other vessels and obstacles on its path, the global path must be adjusted locally to avoid collision. As part of local path planning, USV collision avoidance for moving obstacles is carried out online based on real-time information of shipborne sensors such as radars. This paper studies local path planning for USV collision avoidance.All USVs must have the ability to autonomously avoid obstacles, including land masses, watercraft, low-hanging obstacles, submerged shallow obstacles, interface obstacles, and submerged obstacles for USVs that tow systems [3].There are some results for local path planning in the literature. An obstacle-avoidance solution based on rules and criteria was proposed, in which the collision avoidance rules and criteria are determined according to the target’s collision priority [4]. A method to calculate the collision probabilities of vessels was established under the constraints of mission space and number of vessels [5]. An alternative method assessed the collision risk for surface ships in close-range encounters that is compliant with the International Regulations for Preventing Collisions at Sea (COLREGs) [6].The velocity obstacle (VO) method is a local path planning method that considers the velocity of obstacles [7]. The method has been used to study the case of a single vehicle avoiding some obstacles moving along known linear trajectories. Using the relative velocity between the vehicle and each obstacle, the dynamic problem is transformed into a number of static problems. These are converted to a single problem by vector transformation, and a set of velocity vectors are calculated to avoid obstacles. To avoid collisions, both the current position and the velocities of vehicles are used to compute their future trajectories [8]. An improved VO method was proposed to choose an optimal velocity by introducing a cost function, which consists of a pass time and clearance [9]. A conflict-free maneuver that aims to avoid obstacles while limiting the deviation from the original route was generated by the selective VO method [10]. A finite-time velocity obstacle (FVO) method was used to plan the motion of a mobile robot and complete a crossing through an unknown environment [11]. Based on the nonconvex nature of the FVO constraints, a parabolic function and the approximated function were used to calculate the optimal next-step velocity of a mobile robot to ensure collision-free motion.In [12, 13], considering the kinematic constraints of the robot, a dynamic window approach (DWA) was proposed, which calculates the velocity that the robot can reach within a given time interval. The reachable velocities constitute a set of velocity space, known as dynamic window. The path planner selects the optimal feasible velocity vector from the dynamic window by the objective functions. DWA was proposed to plan an energy efficient local path for navigation of omnidirectional battery-powered mobile robots in dynamic environments with incorporating a cost function based on energy consumption [14].The COLREGs are maritime traffic regulations established by the International Maritime Organization (IMO) to prevent collisions between vessels. The COLREGs specify several collision situations during a voyage: head-on, crossing, and overtaking [15]. It is reasonable for a USV to follow COLREGs during navigation. Some researchers have applied COLREGs to actual collision avoidance. When a ship was detected, the range of its position during a given time frame was estimated, and the path planning system, considering COLREGs, produced a new, safe, and smooth path in real time [16]. Based on the motion velocity model and COLREGs, a reverse eccentric expansion method was designed to deal with dynamic obstacles [17]. A multiobjective optimization approach for path planning of an autonomous surface vehicle combined COLREGs with good seaman’s practice and stratified the objectives [18].Some intelligent optimization algorithms have been applied to the path planning of vehicles. These include ant colony optimization [19], genetic algorithms [20], fuzzy algorithms [21], and particle swarm optimization (PSO) [22]. The algorithms have global optimization ability and overcome the shortcomings of many traditional path planning algorithms. As a method to solve optimization problems, PSO has been developed rapidly in recent years and widely used in many fields, such as combinatorial optimization [23], scheduling problems [24], and neural networks [25].With the development of quantum theory, to combine quantum computing with intelligent optimization algorithms has been proposed. Quantum-behaved particle swarm optimization (QBPSO) was first proposed to improve the global search ability of the original PSO [26]. The position of a particle in PSO was sampled by the quantum delta potential well model around the mean best position [27]. The basic principle of QBPSO is to regard the optimization process as a moving process of particles in the potential field of quantum mechanics to the lowest point of potential energy (the center of the potential well) [28]. Compared with PSO, the iterative equation of QBPSO needs no velocity vectors of particles, is easier to tune, and requires fewer parameters to implement [29].Quantum mechanics principles and evolutionary computing methods were first combined in [30]. A quantum bit and superposition of states were used to solve the knapsack problem by a quantum-inspired evolutionary algorithm (QEA) [31]. Based on the QEA with a quantum rotation gate strategy, an adaptive evolution-based quantum-inspired evolutionary algorithm (AEQEA) introduced an adaptive evolution mechanism [32]. A quantum ant colony algorithm was used to plan the global path of a USV [33, 34]. A quantum evolutionary algorithm was integrated into particle swarm optimization [35]. The positions of particles were encoded by probability amplitudes of quantum bits, and the movements of particles were performed by quantum rotation gates to achieve particle searching. Quantum-inspired particle swarm optimization (QIPSO) has a stronger search ability and quicker convergence speed, benefiting from quantum computing theory, self-adaptive probability selection, and chaotic sequence mutation [36].This study proposes a local path planning algorithm for USV collision avoidance based on modified quantum particle swarm optimization (MQPSO). The main contributions are as follows:(1)
Based on the model proposed in [7, 17], a USV collision avoidance model is established by the VO method considering the constraints of USV kinematic model, the rules of COLREGs, and the uncertainty of obstacle’s velocity. The collision avoidance model not only considers the velocity and course of the USV but also handles the variable velocity and course of an obstacle. According to the collision avoidance model, the USV needs to adjust its velocity and course simultaneously to avoid collision and minimize its energy consumption. Due to the kinematic constraints of the USV, the velocity window and course window of the USV are determined by DWA. Therefore, the local path planning for USV collision avoidance is transformed into a multiobjective optimization problem with multiple constraints in a continuous search space.(2)
In order to solve the optimal value in a continuous search space that satisfies a series of constraints, an optimization algorithm is needed. Because the particles in QBPSO are encoded by real numbers in [26–29], the search efficiency is low. The motion states of particles are represented by a vector of velocity and position [36]. To avoid premature convergence and local optima, a new optimization algorithm called modified quantum particle swarm optimization (MQPSO), combining QBPSO and quantum computing, is proposed in this study. The quantum bits (Q-bits) are used to encode the positions of particles, and the motion states of particles are determined by wave functions. A quantum nongate is introduced to increase the population’s diversity. Since USV local path planning is completed in a short time, it is essential that the optimization algorithm can quickly obtain the optimal value. MQPSO meets that requirement.The remainder of this study is organized as follows. In Section2, the model of USV local path planning is established, and the USV kinematic model, the rules of COLREGs, the VO method, the velocity window and course window of the USV, the collision avoidance model with considering uncertainties, and the cost function of path planning are described. In Section 3, the principles of MQPSO are provided, and it is applied to plan a USV local path. In Section 4, simulations for testing the performance of MQPSO and USV local path planning using MQPSO are presented. Conclusions are provided in Section 5.


## 2. Problem Statement

In this section, the USV kinematic model, rules of COLREGs, velocity obstacle method, and cost function of path planning are introduced.


### 2.1. USV Kinematic Model

A three-degree-of-freedom model is usually adopted for USV path planning. The kinematic equation can be expressed according as [16, 37](1)η˙=Rψν+Vc,where η=x,y,ψT∈ℝ3 is the vector including the positions and heading angle of the USV in the north-east coordinate system. ν=ur,vr,rT is the relative velocity between the USV and the ocean current in the body-fixed reference frame, where ur=u−uc and vr=v−vc. u and v are the velocities in surge and sway, respectively, of the USV. Let Vc=Vcx,Vcy,0T∈ℝ3 represent the velocity of the ocean current in the north-east coordinate system. The body-fixed ocean current velocities uc,vc and north-east current velocities Vcx,Vcy satisfy uc,vcT=RTψVcx,VcyT, and the rotation matrix Rψ is(2)Rψ=cosψ−sinψ0sinψcosψ0001.Therefore, the USV moving in surge, sway, and yaw can be described as(3)x˙=urcosψ−vrsinψ+Vcx,y˙=ursinψ−vrcosψ+Vcy,ψ˙=r.The relative velocity of the USV is(4)U=ur2+vr2.The USV at positionx,y moving in the sea has velocity(5)vusv=x˙2+y˙2.


### 2.2. Collision Avoidance Model

#### 2.2.1. VO-Based Modeling

The collision avoidance model is established as follows. Figure1 shows the schematic of obstacle-avoidance modeling. Take the current position of the USV, that is, point V, whose coordinates are (xv, yv), as the origin of the north-east coordinate system, and point O, whose coordinates are (xo, yo), as the current position of the obstacle. The USV and the obstacle are moving. vusv and vobs are the velocities of the USV and obstacle, respectively. Thus, Δv=vusv−vobs is the velocity difference between the USV and the obstacle. α is the angle from the X-axis to vusv, which is represented as the course angle of the USV. β is the angle from the X-axis to vobs. θ is the angle from the X-axis to the line VO. ϕ is the angle from Δv to vusv. Although the shape of the obstacle is unknown, the obstacle’s collision area can be expanded to a circle. With the point O as the center and robs as the radius, the circle of the collision area is constructed and represented by ⊙O. The lines VM and VN are two tangents of the collision area ⊙O. γ is the angle from the line VO to Δv. μ is the angle from the line VO to VM, which is called the minimum safety angle.Figure 1
Schematic of collision avoidance modeling.Based on the current velocity of the obstacle, to avoid collision, the vectorΔv cannot be in the polygon VMON in the next time interval. Therefore, γ+Δγ>μ is the collision avoidance criterion.The angleγ is related to Δv. Δv is decomposed into velocity components Δvo, directed to the obstacle, and Δve, perpendicular to Δvo. Thus, Δvo and Δve are represented as(6)Δve=vusvsinα−θ−vobssinβ−θ,Δvo=vusvcosα−θ−vobscosβ−θ.As shown in Figure1, the angle γ is represented as(7)γ=arctanvusvsinα−θ−vobssinβ−θvusvcosα−θ−vobscosβ−θ.The difference ofγ in the time interval Δt can be calculated as(8)Δγ=−vobssinα−βvusv2+vobs2+2vusvvobscosα−βΔvusv+vusv2−vusvvobscosα−βvusv2+vobs2+2vusvvobscosα−βΔα+−vusvsinα−βvusv2+vobs2+2vusvvobscosα−βΔvobs+vobs2−vusvvobscosα−βvusv2+vobs2+2vusvvobscosα−βΔβ.Since there is a relationship betweenvusv, vobs, and Δv, as shown in Figure 2, the following equations hold:(9)vobssinϕ=Δvsinα−β,vobssinϕ=vusvsinα−β+ϕ,vusv=vobscosα−β+Δvcosϕ,vobs=vusvcosα−β−Δvcosα−β+ϕ,Δv2=vusv2+vobs2−2vusvvobscosα−β.Figure 2
Positional relationship betweenvusv, vobs, and Δv.Substituting (9) in (8) yields(10)Δγ=−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ,where α, β, and vobs are measured by the shipborne sensors. Δvobs and Δβ can also be measured. μ and ϕ are calculated as follows:(11)μ=arcsinrobsxv−xo2+yv−yo2,(12)ϕ=arcsinvobssinα−βvusv−vobs.As shown in Figure1, if γ<μ, then the USV will collide with the obstacle, so Δγ must be adjusted to satisfy γ>μ, that is,(13)γ+Δγ>μ.Therefore, the adjustment range ofΔγ is(14)Δγ>μ−γ,γ>0,Δγ<−μ−γ,γ≤0.It can be concluded from (8) and (10) that Δγ is related to Δvusv, Δα, Δvobs, and Δβ. But only Δvusv and Δα are control variables of the USV. Thus, the local path planning for USV collision avoidance is transformed to solve the velocity variation Δvusv and course variation Δα of the USV.As mentioned above, the USV adjusts its velocityvusv and the course α within a given time interval Δt to avoid collision. However, due to the constraints of USV kinematics, the velocity and course of the USV will not change much in a given time interval. DWA was introduced to handle the kinematic constraints of the USV. The DWA considers the kinematics and calculates the velocity and course that the USV can achieve within a given time window, which is called the dynamic window VDW, respectively:(15)VDW=vusv,αvusv∈vcur−ΔvusvΔt,vcur+ΔvusvΔt,α∈αcur−ΔαΔt,αcur+ΔαΔt,where vcur is the current velocity of the USV and αcur is the current course of the USV [12].The velocityvusv and course α reached by USV within a given Δt can be determined by (15). The elements in the set VDW are continuous. Therefore, the local path planner needs to use optimization algorithm to obtain the optimized Δvusv and Δα that meet the optimization objectives and constraints within the feasible continuous space of velocity and course, and then a new USV velocity vector vusv can be obtained to avoid collision.


#### 2.2.2. Uncertainty of the Obstacle

In actual collision avoidance, the path planner must consider various types of uncertainties. The uncertainty in USV local path planning is mainly divided into two categories: the uncertainty of sensor detection and the uncertainty of obstacle motion.It is assumed that moving obstacles are detected and tracked by shipborne sensors. Performance characteristics of the sensors in turn affect the noise and state estimation errors for the tracked obstacle. Another uncertainty is from the imprecise motion of the moving obstacle [38]. The velocities of moving obstacle are assumed to be constant in traditional VO, but in reality they do not move with constant velocities. Considering the uncertainty of the obstacle, the velocity of the obstacle is calculated as(16)vobs=v¯obs+δobs,where v¯obs is the nominal velocity (i.e., the expected velocity that is estimated by the obstacle tracker) and δobs represents the uncertainties of the obstacle’s velocity [38]. We assumed that the uncertain component of the velocity lies in a set δobs∈⊙wobs, where ⊙wobs is a bound set and treated as a constant. Since Δv=vusv−vobs, the vector Δv is also affected by the uncertainty and cannot be in the polygon VMON in the next time interval. Due to the uncertainty of vobs, the collision area ⊙O with the worst case uncertainty ⊙WO is defined as(17)⊙WO=⊙O⊕⊙wobs,where ⊕ is the Minkowski sum operation.Figure3 shows the same situation as shown in Figure 2 and takes into account the uncertainties. The circle drawn with a dotted line is represented as ⊙WO. The lines VM′ and VN′ are two tangents of the area ⊙WO. Therefore, considering the uncertainties, the vector Δv will not be in the polygon VM′ON′ to avoid collision in the next time interval.Figure 3
Schematic of collision avoidance modeling considering uncertainty.


### 2.3. Collision Avoidance Rules Based on COLREGs

Rules 13, 14, and 15 of the COLREGs [15] make provisions for three types of collisions that may occur during the sailing of a vessel, that is, head-on, crossing, and overtaking, respectively. The rules are as follows:(1)
Overtaking: any vessel overtaking any other shall keep out of the way of the vessel being overtaken.(2)
Head-on: when two power-driven vessels are meeting on reciprocal or nearly reciprocal courses so as to involve risk of collision, each shall alter her course to starboard so that each shall pass on the port side of the other.(3)
Crossing: when two power-driven vessels are crossing so as to involve risk of collision, the vessel which has the other on her own starboard side shall keep out of the way and shall, if the circumstances of the case admit, avoid crossing ahead of the other vessel.However, the rules are only behavioral constraints on operation and do not explicitly state the angle and scope in practical applications. In this study, we assume that there is no communication between the USV and the obstacle. Due to the special nature of the unmanned operation, the USV is regarded as the give-way vessel and the obstacle is regarded as the stand-on vessel. Therefore, the USV must avoid the obstacle [38]. According to the empirical data of collision avoidance, the angle ranges of the four encounters of overtaking, head-on, left crossing, and right crossing are determined [39], as shown in Figure 4.(1)
Overtaking rule: when the relative course angle between the USV and the obstacle is within112.5°,180°∪−180°,−112.5°, that is,(18)α−β∈112.5°,180°∪−180°,−112.5°,the USV will pass along the port side of the obstacle, as shown in Figure 5(a).(2)
Head-on rule: when the relative course angle between the USV and the obstacle satisfies−15°,15°, that is,(19)180°−α−β∈−15°,15°,the USV will sail along the starboard side of the obstacle, as shown in Figure 5(b).(3)
Crossing from left rule: when the relative course angle between the USV and the obstacle satisfies−112.5°,−15°, that is,(20)α−β∈−112.5°,−15°,the USV will adjust the course to its port side and sail along the tail of the obstacle, as shown in Figure 5(c).(4)
Crossing from right rule: when the relative course angle between the USV and the obstacle satisfies15°,112.5°, that is,(21)α−β∈15°,112.5°,the USV will adjust the course to its starboard side and sail along the tail of the obstacle, as shown in Figure 5(d).Figure 4
Possible encounter situations between USV and moving obstacle.Figure 5
Possible encounter situations and collision avoidance paths between the USV and the obstacle. (a) Overtaking. (b) Head-on. (c) Crossing from left. (d) Crossing from right.
(a)
(b)
(c)
(d)


### 2.4. Optimization Model of Local Path Planning

The mathematical model of the optimization objectives and constraints of local path planning for USV collision avoidance is constructed.


#### 2.4.1. Objective Functions

According to Section2.2, the USV must adjust its velocity and course to avoid collisions. To adjust these simultaneously, it is necessary to obtain the optimal values of the two variables under constraints. Hence, the local path planning for USV collision avoidance is transformed to a multiobjective optimization problem. The energy consumption while sailing determines the USV’s endurance and duration [34]. Therefore, the energy consumption of the USV is considered as one of the objectives of local path planning.A USV sails on a preplanned global path until it detects unsafe obstacles, such that local adjustment must be made to the global path to avoid collision. To ensure the efficiency of the USV’s task execution, it is expected that the motion of the USV should be changed as little as possible during local path planning, that is, the variations of the USV’s velocity and course should be small.In local path planning for USV collision avoidance, the velocity variation, course variation, and energy consumption of the USV during time intervalΔt are used as the optimization objectives.Suppose the velocity of the USV at timet is vusvt, and after time interval Δt, its velocity is vusvt+Δt. The velocity variation Jv is defined as(22)Jv=Δvusv=vusvt+Δt−vusvt.Suppose the course of the USV at timet is αt, and after time interval Δt, its course is αt+Δt. Consequently, the course variation Jα is(23)Jα=Δα=αt+Δt−αt.The energy consumption of the USV while sailing is derived from the propulsion system. The resistance determines the effective power of the vessel [40]. The hydrodynamic drag Fd is calculated as(24)Fd=0.5ρU2CDA,where ρ is the density of the water, CD is the drag coefficient, A is the reference area, and U is obtained by (4).Suppose the USV’s energy consumption after time intervalΔt is(25)Je=FdUΔt=0.5ρU3CDAΔt.In summary, the cost function of USV local path planning is established as(26)J=minw1Jv+w2Jα+w3Je,where w1, w2, and w3 are the weights of Jv, Jα, and Je, respectively.


#### 2.4.2. Constraints of Local Path Planning

To plan an optimal local path for a USV, some constraints, such as kinematic constraints, COLREGs rules, and obstacle-avoidance constraints, must be met. The kinematic constraints are the limitations of the USV’s velocity and course variation. The obstacle-avoidance constraints require the USV to be outside the collision area, and these are obtained from (14).Therefore, the constraints of USV local path planning are as follows:(27)0<vmin≤vusv≤vmax,(28)Δvusv≤a,(29)Δα≤r·Δt,(30)−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ>μ−γ,γ>0,(31)−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ<−μ−γ,γ≤0,where vmax and vmin are, respectively, the upper and lower bounds of vusv and a is the limit of the size of Δvusv. (27), (28), and (29) are constraints of USV kinematics. (30) and (31) are collision avoidance constraints. In addition, (18)–(21) are the constraints based on COLREGs.In summary, the local path planning for USV collision avoidance is transformed into a multiobjective optimization problem with multiple constraints in a continuous search space.


## 2.1. USV Kinematic Model

A three-degree-of-freedom model is usually adopted for USV path planning. The kinematic equation can be expressed according as [16, 37](1)η˙=Rψν+Vc,where η=x,y,ψT∈ℝ3 is the vector including the positions and heading angle of the USV in the north-east coordinate system. ν=ur,vr,rT is the relative velocity between the USV and the ocean current in the body-fixed reference frame, where ur=u−uc and vr=v−vc. u and v are the velocities in surge and sway, respectively, of the USV. Let Vc=Vcx,Vcy,0T∈ℝ3 represent the velocity of the ocean current in the north-east coordinate system. The body-fixed ocean current velocities uc,vc and north-east current velocities Vcx,Vcy satisfy uc,vcT=RTψVcx,VcyT, and the rotation matrix Rψ is(2)Rψ=cosψ−sinψ0sinψcosψ0001.Therefore, the USV moving in surge, sway, and yaw can be described as(3)x˙=urcosψ−vrsinψ+Vcx,y˙=ursinψ−vrcosψ+Vcy,ψ˙=r.The relative velocity of the USV is(4)U=ur2+vr2.The USV at positionx,y moving in the sea has velocity(5)vusv=x˙2+y˙2.


## 2.2. Collision Avoidance Model

### 2.2.1. VO-Based Modeling

The collision avoidance model is established as follows. Figure1 shows the schematic of obstacle-avoidance modeling. Take the current position of the USV, that is, point V, whose coordinates are (xv, yv), as the origin of the north-east coordinate system, and point O, whose coordinates are (xo, yo), as the current position of the obstacle. The USV and the obstacle are moving. vusv and vobs are the velocities of the USV and obstacle, respectively. Thus, Δv=vusv−vobs is the velocity difference between the USV and the obstacle. α is the angle from the X-axis to vusv, which is represented as the course angle of the USV. β is the angle from the X-axis to vobs. θ is the angle from the X-axis to the line VO. ϕ is the angle from Δv to vusv. Although the shape of the obstacle is unknown, the obstacle’s collision area can be expanded to a circle. With the point O as the center and robs as the radius, the circle of the collision area is constructed and represented by ⊙O. The lines VM and VN are two tangents of the collision area ⊙O. γ is the angle from the line VO to Δv. μ is the angle from the line VO to VM, which is called the minimum safety angle.Figure 1
Schematic of collision avoidance modeling.Based on the current velocity of the obstacle, to avoid collision, the vectorΔv cannot be in the polygon VMON in the next time interval. Therefore, γ+Δγ>μ is the collision avoidance criterion.The angleγ is related to Δv. Δv is decomposed into velocity components Δvo, directed to the obstacle, and Δve, perpendicular to Δvo. Thus, Δvo and Δve are represented as(6)Δve=vusvsinα−θ−vobssinβ−θ,Δvo=vusvcosα−θ−vobscosβ−θ.As shown in Figure1, the angle γ is represented as(7)γ=arctanvusvsinα−θ−vobssinβ−θvusvcosα−θ−vobscosβ−θ.The difference ofγ in the time interval Δt can be calculated as(8)Δγ=−vobssinα−βvusv2+vobs2+2vusvvobscosα−βΔvusv+vusv2−vusvvobscosα−βvusv2+vobs2+2vusvvobscosα−βΔα+−vusvsinα−βvusv2+vobs2+2vusvvobscosα−βΔvobs+vobs2−vusvvobscosα−βvusv2+vobs2+2vusvvobscosα−βΔβ.Since there is a relationship betweenvusv, vobs, and Δv, as shown in Figure 2, the following equations hold:(9)vobssinϕ=Δvsinα−β,vobssinϕ=vusvsinα−β+ϕ,vusv=vobscosα−β+Δvcosϕ,vobs=vusvcosα−β−Δvcosα−β+ϕ,Δv2=vusv2+vobs2−2vusvvobscosα−β.Figure 2
Positional relationship betweenvusv, vobs, and Δv.Substituting (9) in (8) yields(10)Δγ=−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ,where α, β, and vobs are measured by the shipborne sensors. Δvobs and Δβ can also be measured. μ and ϕ are calculated as follows:(11)μ=arcsinrobsxv−xo2+yv−yo2,(12)ϕ=arcsinvobssinα−βvusv−vobs.As shown in Figure1, if γ<μ, then the USV will collide with the obstacle, so Δγ must be adjusted to satisfy γ>μ, that is,(13)γ+Δγ>μ.Therefore, the adjustment range ofΔγ is(14)Δγ>μ−γ,γ>0,Δγ<−μ−γ,γ≤0.It can be concluded from (8) and (10) that Δγ is related to Δvusv, Δα, Δvobs, and Δβ. But only Δvusv and Δα are control variables of the USV. Thus, the local path planning for USV collision avoidance is transformed to solve the velocity variation Δvusv and course variation Δα of the USV.As mentioned above, the USV adjusts its velocityvusv and the course α within a given time interval Δt to avoid collision. However, due to the constraints of USV kinematics, the velocity and course of the USV will not change much in a given time interval. DWA was introduced to handle the kinematic constraints of the USV. The DWA considers the kinematics and calculates the velocity and course that the USV can achieve within a given time window, which is called the dynamic window VDW, respectively:(15)VDW=vusv,αvusv∈vcur−ΔvusvΔt,vcur+ΔvusvΔt,α∈αcur−ΔαΔt,αcur+ΔαΔt,where vcur is the current velocity of the USV and αcur is the current course of the USV [12].The velocityvusv and course α reached by USV within a given Δt can be determined by (15). The elements in the set VDW are continuous. Therefore, the local path planner needs to use optimization algorithm to obtain the optimized Δvusv and Δα that meet the optimization objectives and constraints within the feasible continuous space of velocity and course, and then a new USV velocity vector vusv can be obtained to avoid collision.


### 2.2.2. Uncertainty of the Obstacle

In actual collision avoidance, the path planner must consider various types of uncertainties. The uncertainty in USV local path planning is mainly divided into two categories: the uncertainty of sensor detection and the uncertainty of obstacle motion.It is assumed that moving obstacles are detected and tracked by shipborne sensors. Performance characteristics of the sensors in turn affect the noise and state estimation errors for the tracked obstacle. Another uncertainty is from the imprecise motion of the moving obstacle [38]. The velocities of moving obstacle are assumed to be constant in traditional VO, but in reality they do not move with constant velocities. Considering the uncertainty of the obstacle, the velocity of the obstacle is calculated as(16)vobs=v¯obs+δobs,where v¯obs is the nominal velocity (i.e., the expected velocity that is estimated by the obstacle tracker) and δobs represents the uncertainties of the obstacle’s velocity [38]. We assumed that the uncertain component of the velocity lies in a set δobs∈⊙wobs, where ⊙wobs is a bound set and treated as a constant. Since Δv=vusv−vobs, the vector Δv is also affected by the uncertainty and cannot be in the polygon VMON in the next time interval. Due to the uncertainty of vobs, the collision area ⊙O with the worst case uncertainty ⊙WO is defined as(17)⊙WO=⊙O⊕⊙wobs,where ⊕ is the Minkowski sum operation.Figure3 shows the same situation as shown in Figure 2 and takes into account the uncertainties. The circle drawn with a dotted line is represented as ⊙WO. The lines VM′ and VN′ are two tangents of the area ⊙WO. Therefore, considering the uncertainties, the vector Δv will not be in the polygon VM′ON′ to avoid collision in the next time interval.Figure 3
Schematic of collision avoidance modeling considering uncertainty.


## 2.2.1. VO-Based Modeling

The collision avoidance model is established as follows. Figure1 shows the schematic of obstacle-avoidance modeling. Take the current position of the USV, that is, point V, whose coordinates are (xv, yv), as the origin of the north-east coordinate system, and point O, whose coordinates are (xo, yo), as the current position of the obstacle. The USV and the obstacle are moving. vusv and vobs are the velocities of the USV and obstacle, respectively. Thus, Δv=vusv−vobs is the velocity difference between the USV and the obstacle. α is the angle from the X-axis to vusv, which is represented as the course angle of the USV. β is the angle from the X-axis to vobs. θ is the angle from the X-axis to the line VO. ϕ is the angle from Δv to vusv. Although the shape of the obstacle is unknown, the obstacle’s collision area can be expanded to a circle. With the point O as the center and robs as the radius, the circle of the collision area is constructed and represented by ⊙O. The lines VM and VN are two tangents of the collision area ⊙O. γ is the angle from the line VO to Δv. μ is the angle from the line VO to VM, which is called the minimum safety angle.Figure 1
Schematic of collision avoidance modeling.Based on the current velocity of the obstacle, to avoid collision, the vectorΔv cannot be in the polygon VMON in the next time interval. Therefore, γ+Δγ>μ is the collision avoidance criterion.The angleγ is related to Δv. Δv is decomposed into velocity components Δvo, directed to the obstacle, and Δve, perpendicular to Δvo. Thus, Δvo and Δve are represented as(6)Δve=vusvsinα−θ−vobssinβ−θ,Δvo=vusvcosα−θ−vobscosβ−θ.As shown in Figure1, the angle γ is represented as(7)γ=arctanvusvsinα−θ−vobssinβ−θvusvcosα−θ−vobscosβ−θ.The difference ofγ in the time interval Δt can be calculated as(8)Δγ=−vobssinα−βvusv2+vobs2+2vusvvobscosα−βΔvusv+vusv2−vusvvobscosα−βvusv2+vobs2+2vusvvobscosα−βΔα+−vusvsinα−βvusv2+vobs2+2vusvvobscosα−βΔvobs+vobs2−vusvvobscosα−βvusv2+vobs2+2vusvvobscosα−βΔβ.Since there is a relationship betweenvusv, vobs, and Δv, as shown in Figure 2, the following equations hold:(9)vobssinϕ=Δvsinα−β,vobssinϕ=vusvsinα−β+ϕ,vusv=vobscosα−β+Δvcosϕ,vobs=vusvcosα−β−Δvcosα−β+ϕ,Δv2=vusv2+vobs2−2vusvvobscosα−β.Figure 2
Positional relationship betweenvusv, vobs, and Δv.Substituting (9) in (8) yields(10)Δγ=−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ,where α, β, and vobs are measured by the shipborne sensors. Δvobs and Δβ can also be measured. μ and ϕ are calculated as follows:(11)μ=arcsinrobsxv−xo2+yv−yo2,(12)ϕ=arcsinvobssinα−βvusv−vobs.As shown in Figure1, if γ<μ, then the USV will collide with the obstacle, so Δγ must be adjusted to satisfy γ>μ, that is,(13)γ+Δγ>μ.Therefore, the adjustment range ofΔγ is(14)Δγ>μ−γ,γ>0,Δγ<−μ−γ,γ≤0.It can be concluded from (8) and (10) that Δγ is related to Δvusv, Δα, Δvobs, and Δβ. But only Δvusv and Δα are control variables of the USV. Thus, the local path planning for USV collision avoidance is transformed to solve the velocity variation Δvusv and course variation Δα of the USV.As mentioned above, the USV adjusts its velocityvusv and the course α within a given time interval Δt to avoid collision. However, due to the constraints of USV kinematics, the velocity and course of the USV will not change much in a given time interval. DWA was introduced to handle the kinematic constraints of the USV. The DWA considers the kinematics and calculates the velocity and course that the USV can achieve within a given time window, which is called the dynamic window VDW, respectively:(15)VDW=vusv,αvusv∈vcur−ΔvusvΔt,vcur+ΔvusvΔt,α∈αcur−ΔαΔt,αcur+ΔαΔt,where vcur is the current velocity of the USV and αcur is the current course of the USV [12].The velocityvusv and course α reached by USV within a given Δt can be determined by (15). The elements in the set VDW are continuous. Therefore, the local path planner needs to use optimization algorithm to obtain the optimized Δvusv and Δα that meet the optimization objectives and constraints within the feasible continuous space of velocity and course, and then a new USV velocity vector vusv can be obtained to avoid collision.


## 2.2.2. Uncertainty of the Obstacle

In actual collision avoidance, the path planner must consider various types of uncertainties. The uncertainty in USV local path planning is mainly divided into two categories: the uncertainty of sensor detection and the uncertainty of obstacle motion.It is assumed that moving obstacles are detected and tracked by shipborne sensors. Performance characteristics of the sensors in turn affect the noise and state estimation errors for the tracked obstacle. Another uncertainty is from the imprecise motion of the moving obstacle [38]. The velocities of moving obstacle are assumed to be constant in traditional VO, but in reality they do not move with constant velocities. Considering the uncertainty of the obstacle, the velocity of the obstacle is calculated as(16)vobs=v¯obs+δobs,where v¯obs is the nominal velocity (i.e., the expected velocity that is estimated by the obstacle tracker) and δobs represents the uncertainties of the obstacle’s velocity [38]. We assumed that the uncertain component of the velocity lies in a set δobs∈⊙wobs, where ⊙wobs is a bound set and treated as a constant. Since Δv=vusv−vobs, the vector Δv is also affected by the uncertainty and cannot be in the polygon VMON in the next time interval. Due to the uncertainty of vobs, the collision area ⊙O with the worst case uncertainty ⊙WO is defined as(17)⊙WO=⊙O⊕⊙wobs,where ⊕ is the Minkowski sum operation.Figure3 shows the same situation as shown in Figure 2 and takes into account the uncertainties. The circle drawn with a dotted line is represented as ⊙WO. The lines VM′ and VN′ are two tangents of the area ⊙WO. Therefore, considering the uncertainties, the vector Δv will not be in the polygon VM′ON′ to avoid collision in the next time interval.Figure 3
Schematic of collision avoidance modeling considering uncertainty.


## 2.3. Collision Avoidance Rules Based on COLREGs

Rules 13, 14, and 15 of the COLREGs [15] make provisions for three types of collisions that may occur during the sailing of a vessel, that is, head-on, crossing, and overtaking, respectively. The rules are as follows:(1)
Overtaking: any vessel overtaking any other shall keep out of the way of the vessel being overtaken.(2)
Head-on: when two power-driven vessels are meeting on reciprocal or nearly reciprocal courses so as to involve risk of collision, each shall alter her course to starboard so that each shall pass on the port side of the other.(3)
Crossing: when two power-driven vessels are crossing so as to involve risk of collision, the vessel which has the other on her own starboard side shall keep out of the way and shall, if the circumstances of the case admit, avoid crossing ahead of the other vessel.However, the rules are only behavioral constraints on operation and do not explicitly state the angle and scope in practical applications. In this study, we assume that there is no communication between the USV and the obstacle. Due to the special nature of the unmanned operation, the USV is regarded as the give-way vessel and the obstacle is regarded as the stand-on vessel. Therefore, the USV must avoid the obstacle [38]. According to the empirical data of collision avoidance, the angle ranges of the four encounters of overtaking, head-on, left crossing, and right crossing are determined [39], as shown in Figure 4.(1)
Overtaking rule: when the relative course angle between the USV and the obstacle is within112.5°,180°∪−180°,−112.5°, that is,(18)α−β∈112.5°,180°∪−180°,−112.5°,the USV will pass along the port side of the obstacle, as shown in Figure 5(a).(2)
Head-on rule: when the relative course angle between the USV and the obstacle satisfies−15°,15°, that is,(19)180°−α−β∈−15°,15°,the USV will sail along the starboard side of the obstacle, as shown in Figure 5(b).(3)
Crossing from left rule: when the relative course angle between the USV and the obstacle satisfies−112.5°,−15°, that is,(20)α−β∈−112.5°,−15°,the USV will adjust the course to its port side and sail along the tail of the obstacle, as shown in Figure 5(c).(4)
Crossing from right rule: when the relative course angle between the USV and the obstacle satisfies15°,112.5°, that is,(21)α−β∈15°,112.5°,the USV will adjust the course to its starboard side and sail along the tail of the obstacle, as shown in Figure 5(d).Figure 4
Possible encounter situations between USV and moving obstacle.Figure 5
Possible encounter situations and collision avoidance paths between the USV and the obstacle. (a) Overtaking. (b) Head-on. (c) Crossing from left. (d) Crossing from right.
(a)
(b)
(c)
(d)


## 2.4. Optimization Model of Local Path Planning

The mathematical model of the optimization objectives and constraints of local path planning for USV collision avoidance is constructed.


### 2.4.1. Objective Functions

According to Section2.2, the USV must adjust its velocity and course to avoid collisions. To adjust these simultaneously, it is necessary to obtain the optimal values of the two variables under constraints. Hence, the local path planning for USV collision avoidance is transformed to a multiobjective optimization problem. The energy consumption while sailing determines the USV’s endurance and duration [34]. Therefore, the energy consumption of the USV is considered as one of the objectives of local path planning.A USV sails on a preplanned global path until it detects unsafe obstacles, such that local adjustment must be made to the global path to avoid collision. To ensure the efficiency of the USV’s task execution, it is expected that the motion of the USV should be changed as little as possible during local path planning, that is, the variations of the USV’s velocity and course should be small.In local path planning for USV collision avoidance, the velocity variation, course variation, and energy consumption of the USV during time intervalΔt are used as the optimization objectives.Suppose the velocity of the USV at timet is vusvt, and after time interval Δt, its velocity is vusvt+Δt. The velocity variation Jv is defined as(22)Jv=Δvusv=vusvt+Δt−vusvt.Suppose the course of the USV at timet is αt, and after time interval Δt, its course is αt+Δt. Consequently, the course variation Jα is(23)Jα=Δα=αt+Δt−αt.The energy consumption of the USV while sailing is derived from the propulsion system. The resistance determines the effective power of the vessel [40]. The hydrodynamic drag Fd is calculated as(24)Fd=0.5ρU2CDA,where ρ is the density of the water, CD is the drag coefficient, A is the reference area, and U is obtained by (4).Suppose the USV’s energy consumption after time intervalΔt is(25)Je=FdUΔt=0.5ρU3CDAΔt.In summary, the cost function of USV local path planning is established as(26)J=minw1Jv+w2Jα+w3Je,where w1, w2, and w3 are the weights of Jv, Jα, and Je, respectively.


### 2.4.2. Constraints of Local Path Planning

To plan an optimal local path for a USV, some constraints, such as kinematic constraints, COLREGs rules, and obstacle-avoidance constraints, must be met. The kinematic constraints are the limitations of the USV’s velocity and course variation. The obstacle-avoidance constraints require the USV to be outside the collision area, and these are obtained from (14).Therefore, the constraints of USV local path planning are as follows:(27)0<vmin≤vusv≤vmax,(28)Δvusv≤a,(29)Δα≤r·Δt,(30)−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ>μ−γ,γ>0,(31)−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ<−μ−γ,γ≤0,where vmax and vmin are, respectively, the upper and lower bounds of vusv and a is the limit of the size of Δvusv. (27), (28), and (29) are constraints of USV kinematics. (30) and (31) are collision avoidance constraints. In addition, (18)–(21) are the constraints based on COLREGs.In summary, the local path planning for USV collision avoidance is transformed into a multiobjective optimization problem with multiple constraints in a continuous search space.


## 2.4.1. Objective Functions

According to Section2.2, the USV must adjust its velocity and course to avoid collisions. To adjust these simultaneously, it is necessary to obtain the optimal values of the two variables under constraints. Hence, the local path planning for USV collision avoidance is transformed to a multiobjective optimization problem. The energy consumption while sailing determines the USV’s endurance and duration [34]. Therefore, the energy consumption of the USV is considered as one of the objectives of local path planning.A USV sails on a preplanned global path until it detects unsafe obstacles, such that local adjustment must be made to the global path to avoid collision. To ensure the efficiency of the USV’s task execution, it is expected that the motion of the USV should be changed as little as possible during local path planning, that is, the variations of the USV’s velocity and course should be small.In local path planning for USV collision avoidance, the velocity variation, course variation, and energy consumption of the USV during time intervalΔt are used as the optimization objectives.Suppose the velocity of the USV at timet is vusvt, and after time interval Δt, its velocity is vusvt+Δt. The velocity variation Jv is defined as(22)Jv=Δvusv=vusvt+Δt−vusvt.Suppose the course of the USV at timet is αt, and after time interval Δt, its course is αt+Δt. Consequently, the course variation Jα is(23)Jα=Δα=αt+Δt−αt.The energy consumption of the USV while sailing is derived from the propulsion system. The resistance determines the effective power of the vessel [40]. The hydrodynamic drag Fd is calculated as(24)Fd=0.5ρU2CDA,where ρ is the density of the water, CD is the drag coefficient, A is the reference area, and U is obtained by (4).Suppose the USV’s energy consumption after time intervalΔt is(25)Je=FdUΔt=0.5ρU3CDAΔt.In summary, the cost function of USV local path planning is established as(26)J=minw1Jv+w2Jα+w3Je,where w1, w2, and w3 are the weights of Jv, Jα, and Je, respectively.


## 2.4.2. Constraints of Local Path Planning

To plan an optimal local path for a USV, some constraints, such as kinematic constraints, COLREGs rules, and obstacle-avoidance constraints, must be met. The kinematic constraints are the limitations of the USV’s velocity and course variation. The obstacle-avoidance constraints require the USV to be outside the collision area, and these are obtained from (14).Therefore, the constraints of USV local path planning are as follows:(27)0<vmin≤vusv≤vmax,(28)Δvusv≤a,(29)Δα≤r·Δt,(30)−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ>μ−γ,γ>0,(31)−sinϕΔvΔvusv+vusvcosϕΔvΔα+sinα−β+ϕΔvΔvobs−vobscosα−β+ϕΔvΔβ<−μ−γ,γ≤0,where vmax and vmin are, respectively, the upper and lower bounds of vusv and a is the limit of the size of Δvusv. (27), (28), and (29) are constraints of USV kinematics. (30) and (31) are collision avoidance constraints. In addition, (18)–(21) are the constraints based on COLREGs.In summary, the local path planning for USV collision avoidance is transformed into a multiobjective optimization problem with multiple constraints in a continuous search space.


## 3. Optimization Algorithm

As mentioned above, in order to obtain the velocity variationΔvusv and course variation Δα in a continuous search space that satisfies a series of constraints, an optimization algorithm is needed. In this section, we introduce the MQPSO algorithm, which is used to plan the local path for USV collision avoidance.


### 3.1. Particle Swarm Optimization

PSO is an important method in swarm intelligence. The basic algorithm is as follows [41]. In the n-dimensional space, a population of m particles is represented as X=x1,x2,…,xm. The position and velocity vectors of particle i can be represented as xi=xi1,xi2,…,xinT and vi=vi1,vi2,…,vinT, respectively. The n-dimensional best previous position vector is pi=pi1,pi2,…,pinT, and the position vector of the best among all particles in the population is pg=pg1,pg2,…,pgnT. We substitute xi in the objective function to calculate its fitness. Particle i changes its velocity and position according to the following equations:(32)vidt+1=vidt+c1r1pidt−xidt+c2r2pgdt−xidt,(33)xidt+1=xidt+vidt+1,where d=1,2,…,n; i=1,2,…,m, and t is the current iteration number. r1 and r2 are random numbers in 0,1, and c1 and c2 are positive constants. By the iteration of (32) and (33), each particle in the population gradually approaches the global optimal solution.


### 3.2. Quantum-Behaved Particle Swarm Optimization

In PSO, the velocity obtained by (32) is limited. In (33), the particles’ movement range is also limited, such that particles cannot cover the whole feasible solution space. Therefore, PSO cannot guarantee a global optimal solution [42].QBPSO uses the principle of quantum mechanics to determine the states of particles by wave functions, and the particles are bound by the center of the attractive potential generated by the delta potential well. In this way, particles in a quantum bound state can appear anywhere in the search space with a certain probability density. Therefore, the particles in QBPSO can be searched in the whole feasible solution space and approach global optimality.In QBPSO, each particle is treated as spinless and moving in quantum space, and an individual particle is assumed to move in a delta potential well. In quantum mechanics space, the motion state of a particle is represented by a wave functionψ, which replaces the method described by the velocity vector and position vector of the particle. The probability density that a particle appears at a certain point in the search space is represented by ψ2. After determining the probability distribution function of the particles, a particle’s position is determined by the Monte Carlo method. In the n-dimensional space, a population is composed of m particles, and QBPSO updates related parameters through the following equations:(34)pijt=φjtpbestijt+1−φjtgbestjt,φjt∼U0,1,where 1≤i≤m and 1≤j≤n; pij is the center position of a delta potential well; pbestij is the best previous position of particle i and is called the personal best (pbest) position; and gbestj is the position of the best particle among all particles and is called the global best (gbest) position [26, 43].The mean best positionCt defined by the average of the best previous positions of all the particles is(35)Ct=C1t,C2t,…,Cnt=1m∑i=1mPit=1m∑i=1mPi1t,1m∑i=1mPi2t,…,1m∑i=1mPint.Thus, the position of particlei updates as follows:(36)Xijt+1=pijt±αCjt−Xijtln1uijt,uijt∼U0,1,where α is a contraction-expansion coefficient to adjust the algorithm.


### 3.3. Modified Quantum Particle Swarm Optimization

In QBPSO, particles with quantum behavior search in the whole feasible solution space. However, since the particles in QBPSO are encoded by real numbers, the search is only in one-dimensional space, and the search efficiency is low. Quantum bits are introduced to encode the positions of particles to improve the performance of QBPSO. Since the USV local path planning must be completed in real time, it is essential that the optimization algorithm can quickly obtain the optimal value. Therefore, MQPSO is proposed.


#### 3.3.1. Quantum Code

The quantum bit (Q-bit) is the basic unit in quantum computing. A Q-bit is a system with two possible states,0 and 1. The state of a Q-bit φ is expressed as φ=cosθ·0+sinθ·1, where cosθ and sinθ are the probability amplitudes and θ is the phase angle of φ. cosθ2 and sinθ2 are the probabilities in state 0 and 1, respectively. The state of φ collapses to state 0 with probability cosθ2 and to state 1 with probability sinθ2. Thus, the state of φ is an uncertain superposition state between 0 and 1. When the dimension of an individual Xi is n, Xi is expressed as(37)Xi=cosθi1sinθi1cosθi2sinθi2⋯⋯cosθinsinθin,where Xic=cosθi1,cosθi2,…,cosθin and Xis=sinθi1,sinθi2,…,sinθin are two sets of solutions for Xi. Therefore, after quantum coding, in each iteration, every individual has two sets of solutions and the search space is doubled [35]. In this way, the ergodicity of the search space can be enhanced and the optimization process can be accelerated when the population size is same. In MQPSO, quantum coding is combined with quantum-behaved particle swarm optimization, that is, the positions of the particles are encoded by Q-bits according to (37).


#### 3.3.2. Solution Space Transformation

In MQPSO, each particle in the population contains 2n Q-bits’ probability amplitudes. Using linear transformation, these probability amplitudes can be mapped from n-dimensional unit space to the solution space of the optimization problem. Suppose the Q-bit of particle i on the jth dimension is cosθij,sinθijT. The corresponding solution space is(38)Xicj=bi1+cosθij+ai1−cosθij2,Xisj=bi1+sinθij+ai1−sinθij2,where the probability amplitude cosθij of quantum state 0 corresponds to Xicj and the probability amplitude sinθij of quantum state 1 corresponds to Xisj [35].


#### 3.3.3. Particle Evaluation

The two sets of solutionsXic,Xis corresponding to particle i are substituted in the cost function to calculate its cost value. For the minimum problem, the smaller the cost value of the particle, the better the corresponding fitness value. The best previous position of particle i on the jth dimension is determined by(39)pbestijt=Xijt,iffXijt<fpbestijt−1,pbestijt−1,iffXijt≥fpbestijt−1,where Xij is the current position of particle i on the jth dimension.The best particle’s positiongbestj among all particles in the population on the jth dimension is(40)g=argmin1≤i≤mfpbestijt,(41)gbestjt=Pgt,where g is the label of the particle in the global best position.


#### 3.3.4. Particle State Update

In MQPSO, the update of the particle’s position is converted to the update of the Q-bit probability amplitude’s phase angle. In this paper, the particle-updating mechanism in QBPSO is used to adjust the phase angleθ of the Q-bit to update the particle’s position. In the n-dimensional space, a population is composed of m particles. At the tth iteration, the best previous position of particle i and the position of the best among all particles in the population are expressed by probability amplitudes as(42)pbestijt=cosθijpt,sinθijptT,(43)gbestjt=cosθgjt,sinθgjtT,where 1≤i≤m and 1≤j≤n, respectively.The phase angle of the Q-bit in the delta potential well center positionCij is(44)θijct=φjtθijpt+1−φjtθgjtφjt∼U0,1,where θijpt and θgjt are, respectively, the phase angles of the best previous position of particle i and the position of the best among all particles.Therefore, at thet+1th iteration, the phase angle of the Q-bit of the position of the ith particle is(45)θijt+1=θijct±α1m∑i=1mθijpt−θijtln1uijtuijt∼U0,1,where m is the number of the particles.In summary, the position of particlei on the jth dimension at the t+1th iteration is expressed as(46)Xijt+1=cosθijt+1sinθijt+1.


#### 3.3.5. Mutation Operation

Because the population loses its diversity during optimization [27], the PSO algorithm tends to fall into local extreme values. So, the population mutation factor, that is, quantum nongate, is introduced to increase the population’s diversity [35]. First, the mutation probability Pi must be determined, and we choose a random number Rndi in 0,1 for each particle. If Rndi<Pi, then we randomly select n/2 Q-bits and switch the probability amplitudes of Q-bits by the quantum non-gate, such as(47)0110cosθijsinθij=sinθijcosθij=cosπ2−θijsinπ2−θij.Since the mutation operation is independent of the best previous position of the particle and the position of the best among all particles, it can increase the diversity of the population.


#### 3.3.6. Pseudocodes of MQPSO

Based on the above, the main steps of MQPSO to obtain the minimum value of the optimization problem are shown as Algorithm1.Algorithm 1: MQPSO.
Input: the number of the particles m,
the dimension of the problemn,
the maximum of iterationsTmax;
Output: optimal solution gbest;
(1)
Initialize the position of each particle,pbest, and gbest according to (37), (42), and (43), respectively;
(2)
t
←
0;
(3)
if
t
≤
T
max
,
then
(4)
t←t+1;
(5)
for each particlei, do
(6)
for each dimensionj, do
(7)
Calculate the particle’s fitness value according to (26);
(8)
Updatepbestijt according to (39);
(9)
Updategbestt according to (40) and (41);
(10)
Calculateθijct according to (44);
(11)
Updateθijt+1 according to (45);
(12)   Determine the new position of particlei according to (46);
(13)
Output the optimal solutiongbest.


### 3.4. Local Path Planning Based on MQPSO

As described in Section2, the local path planning for USV collision avoidance is transformed to a multiobjective optimization problem with multiple constraints. First, we formulate the optimization model of USV local path planning according to (18)–(21) and (22)–(31). Without loss of generality, it is assumed that the USV detects the obstacle that will affect its safe navigation at the waypoint P1 when t=0. According to the optimization model, MQPSO can obtain the velocity variation Δvusvt and course variation Δαt. Then, the USV sails at velocity vusvt+Δt=vusvt+Δvusvt and course αt+Δt=αt+Δαt to reach the waypoint P2. We repeat this process until the USV has completed collision avoidance. The waypoints are connected to form a local path for USV collision avoidance.


## 3.1. Particle Swarm Optimization

PSO is an important method in swarm intelligence. The basic algorithm is as follows [41]. In the n-dimensional space, a population of m particles is represented as X=x1,x2,…,xm. The position and velocity vectors of particle i can be represented as xi=xi1,xi2,…,xinT and vi=vi1,vi2,…,vinT, respectively. The n-dimensional best previous position vector is pi=pi1,pi2,…,pinT, and the position vector of the best among all particles in the population is pg=pg1,pg2,…,pgnT. We substitute xi in the objective function to calculate its fitness. Particle i changes its velocity and position according to the following equations:(32)vidt+1=vidt+c1r1pidt−xidt+c2r2pgdt−xidt,(33)xidt+1=xidt+vidt+1,where d=1,2,…,n; i=1,2,…,m, and t is the current iteration number. r1 and r2 are random numbers in 0,1, and c1 and c2 are positive constants. By the iteration of (32) and (33), each particle in the population gradually approaches the global optimal solution.


## 3.2. Quantum-Behaved Particle Swarm Optimization

In PSO, the velocity obtained by (32) is limited. In (33), the particles’ movement range is also limited, such that particles cannot cover the whole feasible solution space. Therefore, PSO cannot guarantee a global optimal solution [42].QBPSO uses the principle of quantum mechanics to determine the states of particles by wave functions, and the particles are bound by the center of the attractive potential generated by the delta potential well. In this way, particles in a quantum bound state can appear anywhere in the search space with a certain probability density. Therefore, the particles in QBPSO can be searched in the whole feasible solution space and approach global optimality.In QBPSO, each particle is treated as spinless and moving in quantum space, and an individual particle is assumed to move in a delta potential well. In quantum mechanics space, the motion state of a particle is represented by a wave functionψ, which replaces the method described by the velocity vector and position vector of the particle. The probability density that a particle appears at a certain point in the search space is represented by ψ2. After determining the probability distribution function of the particles, a particle’s position is determined by the Monte Carlo method. In the n-dimensional space, a population is composed of m particles, and QBPSO updates related parameters through the following equations:(34)pijt=φjtpbestijt+1−φjtgbestjt,φjt∼U0,1,where 1≤i≤m and 1≤j≤n; pij is the center position of a delta potential well; pbestij is the best previous position of particle i and is called the personal best (pbest) position; and gbestj is the position of the best particle among all particles and is called the global best (gbest) position [26, 43].The mean best positionCt defined by the average of the best previous positions of all the particles is(35)Ct=C1t,C2t,…,Cnt=1m∑i=1mPit=1m∑i=1mPi1t,1m∑i=1mPi2t,…,1m∑i=1mPint.Thus, the position of particlei updates as follows:(36)Xijt+1=pijt±αCjt−Xijtln1uijt,uijt∼U0,1,where α is a contraction-expansion coefficient to adjust the algorithm.


## 3.3. Modified Quantum Particle Swarm Optimization

In QBPSO, particles with quantum behavior search in the whole feasible solution space. However, since the particles in QBPSO are encoded by real numbers, the search is only in one-dimensional space, and the search efficiency is low. Quantum bits are introduced to encode the positions of particles to improve the performance of QBPSO. Since the USV local path planning must be completed in real time, it is essential that the optimization algorithm can quickly obtain the optimal value. Therefore, MQPSO is proposed.


### 3.3.1. Quantum Code

The quantum bit (Q-bit) is the basic unit in quantum computing. A Q-bit is a system with two possible states,0 and 1. The state of a Q-bit φ is expressed as φ=cosθ·0+sinθ·1, where cosθ and sinθ are the probability amplitudes and θ is the phase angle of φ. cosθ2 and sinθ2 are the probabilities in state 0 and 1, respectively. The state of φ collapses to state 0 with probability cosθ2 and to state 1 with probability sinθ2. Thus, the state of φ is an uncertain superposition state between 0 and 1. When the dimension of an individual Xi is n, Xi is expressed as(37)Xi=cosθi1sinθi1cosθi2sinθi2⋯⋯cosθinsinθin,where Xic=cosθi1,cosθi2,…,cosθin and Xis=sinθi1,sinθi2,…,sinθin are two sets of solutions for Xi. Therefore, after quantum coding, in each iteration, every individual has two sets of solutions and the search space is doubled [35]. In this way, the ergodicity of the search space can be enhanced and the optimization process can be accelerated when the population size is same. In MQPSO, quantum coding is combined with quantum-behaved particle swarm optimization, that is, the positions of the particles are encoded by Q-bits according to (37).


### 3.3.2. Solution Space Transformation

In MQPSO, each particle in the population contains 2n Q-bits’ probability amplitudes. Using linear transformation, these probability amplitudes can be mapped from n-dimensional unit space to the solution space of the optimization problem. Suppose the Q-bit of particle i on the jth dimension is cosθij,sinθijT. The corresponding solution space is(38)Xicj=bi1+cosθij+ai1−cosθij2,Xisj=bi1+sinθij+ai1−sinθij2,where the probability amplitude cosθij of quantum state 0 corresponds to Xicj and the probability amplitude sinθij of quantum state 1 corresponds to Xisj [35].


### 3.3.3. Particle Evaluation

The two sets of solutionsXic,Xis corresponding to particle i are substituted in the cost function to calculate its cost value. For the minimum problem, the smaller the cost value of the particle, the better the corresponding fitness value. The best previous position of particle i on the jth dimension is determined by(39)pbestijt=Xijt,iffXijt<fpbestijt−1,pbestijt−1,iffXijt≥fpbestijt−1,where Xij is the current position of particle i on the jth dimension.The best particle’s positiongbestj among all particles in the population on the jth dimension is(40)g=argmin1≤i≤mfpbestijt,(41)gbestjt=Pgt,where g is the label of the particle in the global best position.


### 3.3.4. Particle State Update

In MQPSO, the update of the particle’s position is converted to the update of the Q-bit probability amplitude’s phase angle. In this paper, the particle-updating mechanism in QBPSO is used to adjust the phase angleθ of the Q-bit to update the particle’s position. In the n-dimensional space, a population is composed of m particles. At the tth iteration, the best previous position of particle i and the position of the best among all particles in the population are expressed by probability amplitudes as(42)pbestijt=cosθijpt,sinθijptT,(43)gbestjt=cosθgjt,sinθgjtT,where 1≤i≤m and 1≤j≤n, respectively.The phase angle of the Q-bit in the delta potential well center positionCij is(44)θijct=φjtθijpt+1−φjtθgjtφjt∼U0,1,where θijpt and θgjt are, respectively, the phase angles of the best previous position of particle i and the position of the best among all particles.Therefore, at thet+1th iteration, the phase angle of the Q-bit of the position of the ith particle is(45)θijt+1=θijct±α1m∑i=1mθijpt−θijtln1uijtuijt∼U0,1,where m is the number of the particles.In summary, the position of particlei on the jth dimension at the t+1th iteration is expressed as(46)Xijt+1=cosθijt+1sinθijt+1.


### 3.3.5. Mutation Operation

Because the population loses its diversity during optimization [27], the PSO algorithm tends to fall into local extreme values. So, the population mutation factor, that is, quantum nongate, is introduced to increase the population’s diversity [35]. First, the mutation probability Pi must be determined, and we choose a random number Rndi in 0,1 for each particle. If Rndi<Pi, then we randomly select n/2 Q-bits and switch the probability amplitudes of Q-bits by the quantum non-gate, such as(47)0110cosθijsinθij=sinθijcosθij=cosπ2−θijsinπ2−θij.Since the mutation operation is independent of the best previous position of the particle and the position of the best among all particles, it can increase the diversity of the population.


### 3.3.6. Pseudocodes of MQPSO

Based on the above, the main steps of MQPSO to obtain the minimum value of the optimization problem are shown as Algorithm1.Algorithm 1: MQPSO.
Input: the number of the particles m,
the dimension of the problemn,
the maximum of iterationsTmax;
Output: optimal solution gbest;
(1)
Initialize the position of each particle,pbest, and gbest according to (37), (42), and (43), respectively;
(2)
t
←
0;
(3)
if
t
≤
T
max
,
then
(4)
t←t+1;
(5)
for each particlei, do
(6)
for each dimensionj, do
(7)
Calculate the particle’s fitness value according to (26);
(8)
Updatepbestijt according to (39);
(9)
Updategbestt according to (40) and (41);
(10)
Calculateθijct according to (44);
(11)
Updateθijt+1 according to (45);
(12)   Determine the new position of particlei according to (46);
(13)
Output the optimal solutiongbest.


## 3.3.1. Quantum Code

The quantum bit (Q-bit) is the basic unit in quantum computing. A Q-bit is a system with two possible states,0 and 1. The state of a Q-bit φ is expressed as φ=cosθ·0+sinθ·1, where cosθ and sinθ are the probability amplitudes and θ is the phase angle of φ. cosθ2 and sinθ2 are the probabilities in state 0 and 1, respectively. The state of φ collapses to state 0 with probability cosθ2 and to state 1 with probability sinθ2. Thus, the state of φ is an uncertain superposition state between 0 and 1. When the dimension of an individual Xi is n, Xi is expressed as(37)Xi=cosθi1sinθi1cosθi2sinθi2⋯⋯cosθinsinθin,where Xic=cosθi1,cosθi2,…,cosθin and Xis=sinθi1,sinθi2,…,sinθin are two sets of solutions for Xi. Therefore, after quantum coding, in each iteration, every individual has two sets of solutions and the search space is doubled [35]. In this way, the ergodicity of the search space can be enhanced and the optimization process can be accelerated when the population size is same. In MQPSO, quantum coding is combined with quantum-behaved particle swarm optimization, that is, the positions of the particles are encoded by Q-bits according to (37).


## 3.3.2. Solution Space Transformation

In MQPSO, each particle in the population contains 2n Q-bits’ probability amplitudes. Using linear transformation, these probability amplitudes can be mapped from n-dimensional unit space to the solution space of the optimization problem. Suppose the Q-bit of particle i on the jth dimension is cosθij,sinθijT. The corresponding solution space is(38)Xicj=bi1+cosθij+ai1−cosθij2,Xisj=bi1+sinθij+ai1−sinθij2,where the probability amplitude cosθij of quantum state 0 corresponds to Xicj and the probability amplitude sinθij of quantum state 1 corresponds to Xisj [35].


## 3.3.3. Particle Evaluation

The two sets of solutionsXic,Xis corresponding to particle i are substituted in the cost function to calculate its cost value. For the minimum problem, the smaller the cost value of the particle, the better the corresponding fitness value. The best previous position of particle i on the jth dimension is determined by(39)pbestijt=Xijt,iffXijt<fpbestijt−1,pbestijt−1,iffXijt≥fpbestijt−1,where Xij is the current position of particle i on the jth dimension.The best particle’s positiongbestj among all particles in the population on the jth dimension is(40)g=argmin1≤i≤mfpbestijt,(41)gbestjt=Pgt,where g is the label of the particle in the global best position.


## 3.3.4. Particle State Update

In MQPSO, the update of the particle’s position is converted to the update of the Q-bit probability amplitude’s phase angle. In this paper, the particle-updating mechanism in QBPSO is used to adjust the phase angleθ of the Q-bit to update the particle’s position. In the n-dimensional space, a population is composed of m particles. At the tth iteration, the best previous position of particle i and the position of the best among all particles in the population are expressed by probability amplitudes as(42)pbestijt=cosθijpt,sinθijptT,(43)gbestjt=cosθgjt,sinθgjtT,where 1≤i≤m and 1≤j≤n, respectively.The phase angle of the Q-bit in the delta potential well center positionCij is(44)θijct=φjtθijpt+1−φjtθgjtφjt∼U0,1,where θijpt and θgjt are, respectively, the phase angles of the best previous position of particle i and the position of the best among all particles.Therefore, at thet+1th iteration, the phase angle of the Q-bit of the position of the ith particle is(45)θijt+1=θijct±α1m∑i=1mθijpt−θijtln1uijtuijt∼U0,1,where m is the number of the particles.In summary, the position of particlei on the jth dimension at the t+1th iteration is expressed as(46)Xijt+1=cosθijt+1sinθijt+1.


## 3.3.5. Mutation Operation

Because the population loses its diversity during optimization [27], the PSO algorithm tends to fall into local extreme values. So, the population mutation factor, that is, quantum nongate, is introduced to increase the population’s diversity [35]. First, the mutation probability Pi must be determined, and we choose a random number Rndi in 0,1 for each particle. If Rndi<Pi, then we randomly select n/2 Q-bits and switch the probability amplitudes of Q-bits by the quantum non-gate, such as(47)0110cosθijsinθij=sinθijcosθij=cosπ2−θijsinπ2−θij.Since the mutation operation is independent of the best previous position of the particle and the position of the best among all particles, it can increase the diversity of the population.


## 3.3.6. Pseudocodes of MQPSO

Based on the above, the main steps of MQPSO to obtain the minimum value of the optimization problem are shown as Algorithm1.Algorithm 1: MQPSO.
Input: the number of the particles m,
the dimension of the problemn,
the maximum of iterationsTmax;
Output: optimal solution gbest;
(1)
Initialize the position of each particle,pbest, and gbest according to (37), (42), and (43), respectively;
(2)
t
←
0;
(3)
if
t
≤
T
max
,
then
(4)
t←t+1;
(5)
for each particlei, do
(6)
for each dimensionj, do
(7)
Calculate the particle’s fitness value according to (26);
(8)
Updatepbestijt according to (39);
(9)
Updategbestt according to (40) and (41);
(10)
Calculateθijct according to (44);
(11)
Updateθijt+1 according to (45);
(12)   Determine the new position of particlei according to (46);
(13)
Output the optimal solutiongbest.


## 3.4. Local Path Planning Based on MQPSO

As described in Section2, the local path planning for USV collision avoidance is transformed to a multiobjective optimization problem with multiple constraints. First, we formulate the optimization model of USV local path planning according to (18)–(21) and (22)–(31). Without loss of generality, it is assumed that the USV detects the obstacle that will affect its safe navigation at the waypoint P1 when t=0. According to the optimization model, MQPSO can obtain the velocity variation Δvusvt and course variation Δαt. Then, the USV sails at velocity vusvt+Δt=vusvt+Δvusvt and course αt+Δt=αt+Δαt to reach the waypoint P2. We repeat this process until the USV has completed collision avoidance. The waypoints are connected to form a local path for USV collision avoidance.


## 4. Evaluation and Simulation Studies

### 4.1. Performance Evaluation of MQPSO

To examine the effectiveness and efficiency of MQPSO, its performance is compared with that of PSO and QBPSO on five different kinds of test benchmark functions. The five benchmark functions are listed in Table1. They are the Ackley, Griewank, Rastrigin, Rosenbrock, and Schaffer functions as defined as follows.Table 1
Benchmark functions.
Benchmark function
Search limit
f
1
x
=
−
20
exp
−
0.2
1
/
N
∑
i
=
1
N
x
i
2
−
exp
1
/
N
∑
i
=
1
N
cos
2
π
x
i
+
20
+
exp
1
x
i
∈
−
8,8
f
2
x
=
∑
i
=
1
N
x
i
2
/
4000
−
∏
i
=
1
N
cos
x
i
/
i
+
1
x
i
∈
−
600,600
f
3
x
=
∑
i
=
1
N
x
i
2
−
10
cos
2
π
x
i
+
10
x
i
∈
−
5.12
,
5.12
f
4
x
=
∑
i
=
1
N
−
1
100
x
i
2
−
x
i
+
1
+
x
i
−
1
2
x
i
∈
−
8,8
f
5
x
1
,
x
2
=
0.5
+
sin
x
1
2
+
x
2
2
2
−
0.5
/
1
+
0.001
x
1
2
+
x
2
2
2
x
1
,
x
2
∈
−
10,10f
1
x is Ackley function. Ackley function is an n-dimensional function with one narrow global optimum basin and many minor local optima [44]. The function gets the minimum value of 0 at x1,x2,…,xn=0,0,…,0.f
2
x is Griewank function. Griewank function has a ∏i=1Ncosxi/i component causing linkages among variables, thereby making it difficult to reach the global optimum [44]. The function gets the minimum value of 0 at x1,x2,…,xn=0,0,…,0.f
3
x is Rastrigin function. Rastrigin function is a complex multimodel problem with a large number of local optima. There are about 10n local optimum in the range xi∈−5.12,5.12. When attempting to solve Rastrigin function, algorithms may easily fall into a local optimum [44]. The function gets the minimum value of 0 at x1,x2,…,xn=0,0,…,0.f
4
x is Rosenbrock function. Because Rosenbrock function provides little information for the optimization algorithm, it is difficult for the algorithm to identify the search direction and find the optimal solution [28]. The function gets the minimum value of 0 at x1,x2,…,xn=1,1,…,1.f
5
x
1
,
x
2 is Schaffer function. Schaffer function is a two-dimensional complex function. Because the function has strong oscillation characteristics, it is difficult to find the global optimal value [28]. The function gets the minimum value of 0 at x1,x2=1,1.In the test, the maximum number of iterations is 300, the number of particles is 50, and the dimensions of the benchmark functions are 30. The parameters of PSO, QBPSO, and MQPSO are consistent.Since all the above test benchmark functions are minimization problems, the smaller the final value of these functions, the better the performance of the corresponding algorithm. The test results, which contain optimal value, the number of iterations for optimization algorithm to converge to the minimum, the calculation time, and the average calculation time for per iteration, are listed in Table2. The 1st, 5th, and 9th rows of Table 2, respectively, show the optimal values obtained by PSO, QBPSO, and MQPSO for benchmark functions f1 to f5. The 2nd, 6th, and 10th rows show the number of iterations required for PSO, QBPSO, and MQPSO to converge to the minimum. The 3rd, 7th, and 11th rows, respectively, show the calculation time required for the benchmark functions f1 to f5 to perform the entire optimization calculation based on PSO, QBPSO, and MQPSO. The 4th, 8th, and 12th rows, respectively, show the average calculation time taken by each iteration of benchmark functions f1 to f5 based on PSO, QBPSO, and MQPSO.Table 2
Results of three algorithms on benchmark functions.
Algorithm
Performance
f
1
x
f
2
x
f
3
x
f
4
x
f
5
x
1
,
x
2
PSO
Optimal value
1.55
×
10
−
3
5.56
×
10
−
3
5.61
×
10
−
2
2.03
×
10
−
3
6.32
×
10
−
3
Iteration
228
245
258
206
261
Time (s)
1.8619
1.8428
1.9653
1.5378
1.6795
Average time per iteration (s)
0.0062
0.0061
0.0066
0.0051
0.0056
QBPSO
Optimal value
9.37
×
10
−
4
2.05
×
10
−
3
2.56
×
10
−
2
1.25
×
10
−
3
3.12
×
10
−
3
Iteration
215
231
254
197
263
Time (s)
1.4844
1.4688
1.4531
1.1563
1.2031
Average time per iteration (s)
0.0049
0.0049
0.0048
0.0039
0.0040
MQPSO
Optimal value
9.15
×
10
−
4
1.86
×
10
−
3
1.96
×
10
−
2
1.12
×
10
−
3
2.89
×
10
−
3
Iteration
149
158
136
132
142
Time (s)
1.5271
1.5163
1.5534
1.3123
1.4808
Average time per iteration (s)
0.0051
0.0051
0.0052
0.0044
0.0049The data in Table2 show that MQPSO is superior to PSO and QBPSO in both the optimal value and the number of iterations that converge to the minimum.The optimal values obtained by MQPSO and QBPSO are lower than obtained by PSO. Because the motion states of particles in MQPSO and QBPSO are determined by wave functions, the particles can cover the whole feasible solution space. MQPSO uses the above methods to avoid falling into local extreme values. And MQPSO further avoids falling into local extremes by using a quantum nongate to increase the population’s diversity. Therefore, the optimal values obtained by MQPSO are slightly lower than those obtained by QBPSO and PSO.The number of iterations required for MQPSO to converge to the minimum is obviously lower than for QBPSO and PSO. Because the positions of the particles are encoded by Q-bits, the search space is doubled, and the convergence speed is faster.In PSO, the strategy of first calculating the particle’s velocity at the next moment and then updating the positions of the particles is adopted, as shown in (32) and (33). So, the calculation time is long. QBPSO uses a method of directly updating the particles’ positions, which reduces the calculation time. The particle’s position update strategy of MQPSO is based on QBPSO. However, if the number of encoded Q-bits is n, each particle is updated n times in MQPSO. Therefore, the increasing calculation will lead to the longer calculation time. As the calculation time increases, the number of particle’s position updates is also increased, improving the optimization ability of the algorithm.Simulation results show that the optimization capability of MQPSO is indeed better than QBPSO and PSO.


### 4.2. Simulations of USV Local Path Planning

In this section, some simulations of USV local path planning based on MQPSO, QBPSO, and PSO considering the uncertainty of obstacle’s velocity are discussed.It is assumed that the USV can detect the motion states of the obstacle by the shipborne sensors and estimate its motion states, that is,Δvobs and Δβ in (10) can be measured. The initial motion information for the USV and the obstacles in Cases 1–3 is listed in Table 3. If the USV does not avoid collision when it detects an obstacle that affects its safe navigation at t=0s, then the USV will collide with the obstacle at t=20s. The time interval ΔT is 2s. In this way, the USV will plan a local path every 2 seconds, that is, 10 calculations will be performed during the entire collision avoidance.Table 3
Initial motion information for USV and obstacle.
Case 1
Case 2
Case 3
Initial location of USV
0,0
0
,
0
0
,
0
Location of goal
500,500
500
,
500
500
,
500
v
usv
0
18
m
/
s
18
m
/
s
18
m
/
s
α
0
45
°
45
°
45
°
Initial location of obstacle
150,150
500,400
450,40
v
obs
0
5
m
/
s
10
m
/
s
10
m
/
s
β
0
45
°
−
145
°
125
°
Δ
v
obs
0.6
m
/
s
0.8
m
/
s
0.8
m
/
s
Δ
β
0.8594
°
0.8594
°
−
0.8594
°
Encounter situation
Overtaking
Head-on
Crossing from right
Whether to follow COLREGs
Yes
Yes
Yes
Expected time of collision
20s
20s
20sThe radius of the collision area⊙O is 30m, and the radius of ⊙wobs is 9m. The parameters in (25) are set as CD=1 and A=300. In Cases 1–3, the weights w1,w2,w3 are 1000,100,1 in (26). The constraints in (27), (28), and (29), respectively, are 0m/s<vusv≤22m/s, Δvusv≤0.25m/s, and Δα/Δt≤2°.In the simulations, the maximum number of iterations is 200 and the number of particles is 50. A particle represents a candidate solution, so the particles are two-dimensional, where one dimension representsΔvusv and the other represents Δα. The parameters of MQPSO, QBPSO, and PSO are consistent.Figures6–8 show the planned local paths for USV collision avoidance based on MQPSO considering uncertainty in Cases 1–3, respectively. The location of the USV is represented by a blue pentagon every two seconds. The location of the obstacle is represented by a red triangle every two seconds. The location of the USV is represented by a black pentagon at t=20s. The collision area is represented as a dotted red line at t=20s, and the USV global path is represented by a dotted green line.Figure 6
Local path for USV collision avoidance based on MQPSO in Case 1.Figure 7
Local path for USV collision avoidance based on MQPSO in Case 2.Figure 8
Local path for USV collision avoidance based on MQPSO in Case 3.Table4 shows the cost values in each calculation, respectively, obtained by MQPSO, QBPSO, and PSO in Cases 1–3, where J1 represents the first calculation performed by local path planner when t=0 and so on. From the data in Table 4, we can conclude that the cost value obtained by the USV local path planner based on MQPSO is lower than those obtained by QBPSO and PSO. According to the evaluation criteria of (25), the paths obtained by MQPSO are better than the paths obtained by QBPSO and PSO. Therefore, it is reasonable to choose MQPSO as the optimization algorithm of USV local path planning in this study.Table 4
Cost values in each calculation.
Algorithm
J
1
J
2
J
3
J
4
J
5
J
6
J
7
J
8
J
9
J
10
Case 1
MQPSO
1319.75
1305.67
1328.42
1323.27
1322.91
1325.21
1525.87
1595.71
1698.57
1784.21
QBPSO
1406.68
1433.61
1407.88
1383.07
1357.21
1358.81
1600.38
1692.81
1784.61
1844.95
PSO
1420.98
1430.80
1443.26
1424.96
1413.08
1398.57
1634.82
1774.58
1877.89
1918.35
Case 2
MQPSO
1048.58
1037.63
1013.48
1007.50
967.91
937.69
917.59
890.68
845.15
836.10
QBPSO
1138.62
1130.76
1097.33
1071.18
1037.98
1006.69
979.31
952.58
908.22
888.58
PSO
1186.66
1167.49
1162.48
1118.84
1086.32
1058.32
1040.08
1005.78
967.77
949.83
Case 3
MQPSO
992.27
980.16
988.98
967.42
965.07
961.52
948.09
953.24
936.98
932.10
QBPSO
1017.63
1020.35
1022.21
1006.77
1018.66
996.09
984.16
986.49
986.19
972.12
PSO
1097.36
1083.25
1069.46
1063.22
1050.69
1055.83
1053.99
1041.82
1020.96
1021.23Figures6–8 show the paths obtained by MQPSO can avoid the moving obstacle in Cases 1–3. After t=20s, the USV returns to the global path as soon as possible. In Cases 1–3, the strategy of simultaneously adjusting the velocity and course of the USV is adopted.As shown in Table3, the encounter situation of Case 1 between the USV and the obstacle is overtaking. According to the rules in COLREGs and Section 2.3, the USV should pass along the port side of the obstacle. As shown in Figure 6, the local path planner based on MQPSO can obtain a collision-free path that meets the rules of COLREGs on overtaking. The adjusted velocity and course of the USV within 0–20 s in Case 1 are shown in Figure 9(a) and 9(b), respectively. As shown in Table 3, the encounter situation of Case 2 between the USV and the obstacle is head-on. According to the rules in COLREGs and Section 2.3, the USV should sail along the starboard side of the obstacle. As shown in Figure 7, the local path planner based on MQPSO can obtain a collision-free path that meets the rules of COLREGs on head-on. The adjusted velocity and course of the USV within 0–20 s in Case 2 are shown in Figure 9(c) and 9(d), respectively. As shown in Table 3, the encounter situation of Case 3 between the USV and the obstacle is crossing from right. According to the rules in COLREGs and Section 2.3, the USV should adjust the course to its starboard side and sail along the tail of the obstacle. As shown in Figure 8, the local path planner based on MQPSO can obtain a collision-free path that meets the rules of COLREGs on crossing from right. The adjusted velocity and course of the USV within 0–20 s in Case 3 are shown in Figure 9(e) and 9(f), respectively.Figure 9
Adjusted velocity and course based on MQPSO in Cases 1–3.
(a)
(b)
(c)
(d)
(e)
(f)The results verify that the proposed algorithm can effectively plan the local path for USV collision avoidance.


## 4.1. Performance Evaluation of MQPSO

To examine the effectiveness and efficiency of MQPSO, its performance is compared with that of PSO and QBPSO on five different kinds of test benchmark functions. The five benchmark functions are listed in Table1. They are the Ackley, Griewank, Rastrigin, Rosenbrock, and Schaffer functions as defined as follows.Table 1
Benchmark functions.
Benchmark function
Search limit
f
1
x
=
−
20
exp
−
0.2
1
/
N
∑
i
=
1
N
x
i
2
−
exp
1
/
N
∑
i
=
1
N
cos
2
π
x
i
+
20
+
exp
1
x
i
∈
−
8,8
f
2
x
=
∑
i
=
1
N
x
i
2
/
4000
−
∏
i
=
1
N
cos
x
i
/
i
+
1
x
i
∈
−
600,600
f
3
x
=
∑
i
=
1
N
x
i
2
−
10
cos
2
π
x
i
+
10
x
i
∈
−
5.12
,
5.12
f
4
x
=
∑
i
=
1
N
−
1
100
x
i
2
−
x
i
+
1
+
x
i
−
1
2
x
i
∈
−
8,8
f
5
x
1
,
x
2
=
0.5
+
sin
x
1
2
+
x
2
2
2
−
0.5
/
1
+
0.001
x
1
2
+
x
2
2
2
x
1
,
x
2
∈
−
10,10f
1
x is Ackley function. Ackley function is an n-dimensional function with one narrow global optimum basin and many minor local optima [44]. The function gets the minimum value of 0 at x1,x2,…,xn=0,0,…,0.f
2
x is Griewank function. Griewank function has a ∏i=1Ncosxi/i component causing linkages among variables, thereby making it difficult to reach the global optimum [44]. The function gets the minimum value of 0 at x1,x2,…,xn=0,0,…,0.f
3
x is Rastrigin function. Rastrigin function is a complex multimodel problem with a large number of local optima. There are about 10n local optimum in the range xi∈−5.12,5.12. When attempting to solve Rastrigin function, algorithms may easily fall into a local optimum [44]. The function gets the minimum value of 0 at x1,x2,…,xn=0,0,…,0.f
4
x is Rosenbrock function. Because Rosenbrock function provides little information for the optimization algorithm, it is difficult for the algorithm to identify the search direction and find the optimal solution [28]. The function gets the minimum value of 0 at x1,x2,…,xn=1,1,…,1.f
5
x
1
,
x
2 is Schaffer function. Schaffer function is a two-dimensional complex function. Because the function has strong oscillation characteristics, it is difficult to find the global optimal value [28]. The function gets the minimum value of 0 at x1,x2=1,1.In the test, the maximum number of iterations is 300, the number of particles is 50, and the dimensions of the benchmark functions are 30. The parameters of PSO, QBPSO, and MQPSO are consistent.Since all the above test benchmark functions are minimization problems, the smaller the final value of these functions, the better the performance of the corresponding algorithm. The test results, which contain optimal value, the number of iterations for optimization algorithm to converge to the minimum, the calculation time, and the average calculation time for per iteration, are listed in Table2. The 1st, 5th, and 9th rows of Table 2, respectively, show the optimal values obtained by PSO, QBPSO, and MQPSO for benchmark functions f1 to f5. The 2nd, 6th, and 10th rows show the number of iterations required for PSO, QBPSO, and MQPSO to converge to the minimum. The 3rd, 7th, and 11th rows, respectively, show the calculation time required for the benchmark functions f1 to f5 to perform the entire optimization calculation based on PSO, QBPSO, and MQPSO. The 4th, 8th, and 12th rows, respectively, show the average calculation time taken by each iteration of benchmark functions f1 to f5 based on PSO, QBPSO, and MQPSO.Table 2
Results of three algorithms on benchmark functions.
Algorithm
Performance
f
1
x
f
2
x
f
3
x
f
4
x
f
5
x
1
,
x
2
PSO
Optimal value
1.55
×
10
−
3
5.56
×
10
−
3
5.61
×
10
−
2
2.03
×
10
−
3
6.32
×
10
−
3
Iteration
228
245
258
206
261
Time (s)
1.8619
1.8428
1.9653
1.5378
1.6795
Average time per iteration (s)
0.0062
0.0061
0.0066
0.0051
0.0056
QBPSO
Optimal value
9.37
×
10
−
4
2.05
×
10
−
3
2.56
×
10
−
2
1.25
×
10
−
3
3.12
×
10
−
3
Iteration
215
231
254
197
263
Time (s)
1.4844
1.4688
1.4531
1.1563
1.2031
Average time per iteration (s)
0.0049
0.0049
0.0048
0.0039
0.0040
MQPSO
Optimal value
9.15
×
10
−
4
1.86
×
10
−
3
1.96
×
10
−
2
1.12
×
10
−
3
2.89
×
10
−
3
Iteration
149
158
136
132
142
Time (s)
1.5271
1.5163
1.5534
1.3123
1.4808
Average time per iteration (s)
0.0051
0.0051
0.0052
0.0044
0.0049The data in Table2 show that MQPSO is superior to PSO and QBPSO in both the optimal value and the number of iterations that converge to the minimum.The optimal values obtained by MQPSO and QBPSO are lower than obtained by PSO. Because the motion states of particles in MQPSO and QBPSO are determined by wave functions, the particles can cover the whole feasible solution space. MQPSO uses the above methods to avoid falling into local extreme values. And MQPSO further avoids falling into local extremes by using a quantum nongate to increase the population’s diversity. Therefore, the optimal values obtained by MQPSO are slightly lower than those obtained by QBPSO and PSO.The number of iterations required for MQPSO to converge to the minimum is obviously lower than for QBPSO and PSO. Because the positions of the particles are encoded by Q-bits, the search space is doubled, and the convergence speed is faster.In PSO, the strategy of first calculating the particle’s velocity at the next moment and then updating the positions of the particles is adopted, as shown in (32) and (33). So, the calculation time is long. QBPSO uses a method of directly updating the particles’ positions, which reduces the calculation time. The particle’s position update strategy of MQPSO is based on QBPSO. However, if the number of encoded Q-bits is n, each particle is updated n times in MQPSO. Therefore, the increasing calculation will lead to the longer calculation time. As the calculation time increases, the number of particle’s position updates is also increased, improving the optimization ability of the algorithm.Simulation results show that the optimization capability of MQPSO is indeed better than QBPSO and PSO.


## 4.2. Simulations of USV Local Path Planning

In this section, some simulations of USV local path planning based on MQPSO, QBPSO, and PSO considering the uncertainty of obstacle’s velocity are discussed.It is assumed that the USV can detect the motion states of the obstacle by the shipborne sensors and estimate its motion states, that is,Δvobs and Δβ in (10) can be measured. The initial motion information for the USV and the obstacles in Cases 1–3 is listed in Table 3. If the USV does not avoid collision when it detects an obstacle that affects its safe navigation at t=0s, then the USV will collide with the obstacle at t=20s. The time interval ΔT is 2s. In this way, the USV will plan a local path every 2 seconds, that is, 10 calculations will be performed during the entire collision avoidance.Table 3
Initial motion information for USV and obstacle.
Case 1
Case 2
Case 3
Initial location of USV
0,0
0
,
0
0
,
0
Location of goal
500,500
500
,
500
500
,
500
v
usv
0
18
m
/
s
18
m
/
s
18
m
/
s
α
0
45
°
45
°
45
°
Initial location of obstacle
150,150
500,400
450,40
v
obs
0
5
m
/
s
10
m
/
s
10
m
/
s
β
0
45
°
−
145
°
125
°
Δ
v
obs
0.6
m
/
s
0.8
m
/
s
0.8
m
/
s
Δ
β
0.8594
°
0.8594
°
−
0.8594
°
Encounter situation
Overtaking
Head-on
Crossing from right
Whether to follow COLREGs
Yes
Yes
Yes
Expected time of collision
20s
20s
20sThe radius of the collision area⊙O is 30m, and the radius of ⊙wobs is 9m. The parameters in (25) are set as CD=1 and A=300. In Cases 1–3, the weights w1,w2,w3 are 1000,100,1 in (26). The constraints in (27), (28), and (29), respectively, are 0m/s<vusv≤22m/s, Δvusv≤0.25m/s, and Δα/Δt≤2°.In the simulations, the maximum number of iterations is 200 and the number of particles is 50. A particle represents a candidate solution, so the particles are two-dimensional, where one dimension representsΔvusv and the other represents Δα. The parameters of MQPSO, QBPSO, and PSO are consistent.Figures6–8 show the planned local paths for USV collision avoidance based on MQPSO considering uncertainty in Cases 1–3, respectively. The location of the USV is represented by a blue pentagon every two seconds. The location of the obstacle is represented by a red triangle every two seconds. The location of the USV is represented by a black pentagon at t=20s. The collision area is represented as a dotted red line at t=20s, and the USV global path is represented by a dotted green line.Figure 6
Local path for USV collision avoidance based on MQPSO in Case 1.Figure 7
Local path for USV collision avoidance based on MQPSO in Case 2.Figure 8
Local path for USV collision avoidance based on MQPSO in Case 3.Table4 shows the cost values in each calculation, respectively, obtained by MQPSO, QBPSO, and PSO in Cases 1–3, where J1 represents the first calculation performed by local path planner when t=0 and so on. From the data in Table 4, we can conclude that the cost value obtained by the USV local path planner based on MQPSO is lower than those obtained by QBPSO and PSO. According to the evaluation criteria of (25), the paths obtained by MQPSO are better than the paths obtained by QBPSO and PSO. Therefore, it is reasonable to choose MQPSO as the optimization algorithm of USV local path planning in this study.Table 4
Cost values in each calculation.
Algorithm
J
1
J
2
J
3
J
4
J
5
J
6
J
7
J
8
J
9
J
10
Case 1
MQPSO
1319.75
1305.67
1328.42
1323.27
1322.91
1325.21
1525.87
1595.71
1698.57
1784.21
QBPSO
1406.68
1433.61
1407.88
1383.07
1357.21
1358.81
1600.38
1692.81
1784.61
1844.95
PSO
1420.98
1430.80
1443.26
1424.96
1413.08
1398.57
1634.82
1774.58
1877.89
1918.35
Case 2
MQPSO
1048.58
1037.63
1013.48
1007.50
967.91
937.69
917.59
890.68
845.15
836.10
QBPSO
1138.62
1130.76
1097.33
1071.18
1037.98
1006.69
979.31
952.58
908.22
888.58
PSO
1186.66
1167.49
1162.48
1118.84
1086.32
1058.32
1040.08
1005.78
967.77
949.83
Case 3
MQPSO
992.27
980.16
988.98
967.42
965.07
961.52
948.09
953.24
936.98
932.10
QBPSO
1017.63
1020.35
1022.21
1006.77
1018.66
996.09
984.16
986.49
986.19
972.12
PSO
1097.36
1083.25
1069.46
1063.22
1050.69
1055.83
1053.99
1041.82
1020.96
1021.23Figures6–8 show the paths obtained by MQPSO can avoid the moving obstacle in Cases 1–3. After t=20s, the USV returns to the global path as soon as possible. In Cases 1–3, the strategy of simultaneously adjusting the velocity and course of the USV is adopted.As shown in Table3, the encounter situation of Case 1 between the USV and the obstacle is overtaking. According to the rules in COLREGs and Section 2.3, the USV should pass along the port side of the obstacle. As shown in Figure 6, the local path planner based on MQPSO can obtain a collision-free path that meets the rules of COLREGs on overtaking. The adjusted velocity and course of the USV within 0–20 s in Case 1 are shown in Figure 9(a) and 9(b), respectively. As shown in Table 3, the encounter situation of Case 2 between the USV and the obstacle is head-on. According to the rules in COLREGs and Section 2.3, the USV should sail along the starboard side of the obstacle. As shown in Figure 7, the local path planner based on MQPSO can obtain a collision-free path that meets the rules of COLREGs on head-on. The adjusted velocity and course of the USV within 0–20 s in Case 2 are shown in Figure 9(c) and 9(d), respectively. As shown in Table 3, the encounter situation of Case 3 between the USV and the obstacle is crossing from right. According to the rules in COLREGs and Section 2.3, the USV should adjust the course to its starboard side and sail along the tail of the obstacle. As shown in Figure 8, the local path planner based on MQPSO can obtain a collision-free path that meets the rules of COLREGs on crossing from right. The adjusted velocity and course of the USV within 0–20 s in Case 3 are shown in Figure 9(e) and 9(f), respectively.Figure 9
Adjusted velocity and course based on MQPSO in Cases 1–3.
(a)
(b)
(c)
(d)
(e)
(f)The results verify that the proposed algorithm can effectively plan the local path for USV collision avoidance.


## 5. Conclusions

This study proposes a local path planning algorithm for USV collision avoidance based on MQPSO.First, a USV collision avoidance model is established, which is composed of USV kinematics, COLREGs rules, the uncertainty of the obstacle’s velocity, and a VO method, which not only considers the velocity and course of the USV but also handles variable velocities and courses of an obstacle. Due to the kinematic constraints of the USV, the velocity window and course window of the USV are determined by DWA. In this way, local path planning for USV collision avoidance is transformed into a multiobjective optimization problem with multiple constraints in a continuous search space. Thus, a USV can avoid collision by simultaneously adjusting its optimized velocity and course.Second, to solve the optimization problem, MQPSO is primarily proposed. MQPSO is an optimization algorithm combining quantum computing with QBPSO. It benefits from the high efficiency of quantum computing and the optimization ability of QBPSO. Simulation results show that MQPSO converges more quickly than PSO and QBPSO and avoids falling into local extreme values. The effectiveness of the proposed collision avoidance model is verified. Performance evaluation tests on benchmark functions show that MQPSO is superior to PSO and QBPSO in both the optimal value and the number of iterations required for convergence, and the proposed algorithm based on MQPSO can effectively and efficiently plan a collision-free USV path.In future work, the collision avoidance model proposed in this study should consider multiple obstacles that a USV may encounter during navigation. Moreover, the correlation between multiple objectives should be calculated, and MQPSO should accommodate multiobjective processing. The uncertainties of the obstacle’s velocity and course are also considered in future works.



---

*Source: 3095426-2020-04-13.xml*