2.1. Environmental model

Grid maps offer several advantages in map modeling. They discretize the planning area into a finite number of grids, which can effectively store spatial information about the environment and achieve the desired map accuracy by utilizing grids of varying resolutions. Additionally, grid maps significantly simplify the search space, thereby enhancing the efficiency of path planning [Reference Zhang, Guo, Andriukaitis, Li, Królczyk and Li34, Reference Yang, Wang and Wang35].

In the grid map, the grid cell occupied by the USV serves as the central node, while the adjacent cells are designated as alternative nodes for subsequent movement. The USV must select a node that is free of obstacles and exhibits either the lowest energy consumption in windy conditions or the shortest distance to the target location in calm conditions. Additionally, if an obstacle does not completely occupy a grid cell, that node is included in the contraindication table to ensure sufficient safety margins for the USV, thereby allowing it to avoid potentially hazardous situations such as shallow waters or reefs. The relationship between the index value S and its corresponding coordinates $(x,y)$ within the grid environment is represented as follows:

(1) \begin{equation} \begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}\mathrm{mod}\left (S-1,n\right )\\-\mathrm{ceil}\left (S/n\right )\end{bmatrix}+\begin{bmatrix}0\\n\end{bmatrix}+\begin{bmatrix}0.5\\0.5\end{bmatrix} \end{equation}

where the modulus operation, denoted as ‘mod’, yields the remainder following the division of two numbers. Conversely, the ceiling function, denoted as ‘ceil’, rounds a given number upwards to the nearest integer that is greater than or equal to the original value.

2.2. Wind load model

The water surface resistance experienced by the USV during travel comprises friction resistance, eddy resistance, and wave resistance. Friction resistance is proportional to the square of the vehicle’s speed, while both eddy resistance and wave resistance increase with rising speed. Additionally, the water surface resistance is influenced by the structural parameters of the hull [Reference Lv, Yu, Chen, Zhao and Chi36]. Given that the water surface environment in the berthing area is relatively stable, the water surface resistance is treated as a variable that is speed-dependent in the energy consumption calculations [Reference Monex and Kumar37]. Taking into account the effects of wind load and the power of the USV, three degrees of freedom are analyzed: surge, sway, and yaw [Reference Wang and Li38]. The forces and moments associated with these three degrees of freedom from the incompressible Bernoulli equation [Reference Gu, Zhou, Wen, Xiao and Knoll39] are derived as follows:

(2) \begin{equation} \begin{cases}X_{wind}=\frac {1}{2}\rho _{a}S_{fp}V_{vw}^{2}A_{wx}\left (\theta _{w}\right )\\[3pt] Y_{wind}=\frac {1}{2}\rho _{a}S_{sp}V_{vw}^{2}A_{wy}\left (\theta _{w}\right )\\[3pt] N_{wind}=\frac {1}{2}\rho _{a}S_{sp}V_{vw}^{2}A_{wn}\left (\theta _{w}\right )L_{oa} \end{cases} \end{equation}

where $X_{wind}$ , $Y_{wind}$ , $N_{wind}$ represent the transverse force, longitudinal force, and overall force acting on the hull, respectively. $\rho _a$ denotes the air density, while $S_{fp}$ and $S_{sp}$ signify the projected area of the hull. $V_{vw}$ indicates the relative wind speed, and $L_{oa}$ represents the total length of the hull. Additionally, $A_{wx}(\theta _{w})$ , $A_{wy}(\theta _{w})$ , and $A_{wn}(\theta _{w})$ are the wind resistance coefficients in the transverse, longitudinal, and vertical directions, respectively. Lastly, $\theta _{w}$ is the angle between the apparent wind (relative wind) and the USV.

Figure 1.

Coordinate system.

Figure 2.

Data storage model.

Non-constant wind was employed as a disturbance factor, with varying wind speeds and directions at different locations. As illustrated in Figure 1, the true wind direction is defined, with $0^\circ$ representing north winds, $90^\circ$ representing east winds, and continuing in a clockwise manner, ranging from $0^\circ$ to $360^\circ$ . As illustrated in Figure 2, a three-dimensional data storage model is established using three channels to store environmental information, specifically obstacle data, wind speed, and wind direction. This approach enhances the fitting of distribution patterns across various wind speed intervals. The vector relationship between the true wind and the apparent wind is described as follows:

(3) \begin{equation} V_{vw} + V_{usv} = V_{rw} \end{equation}

where the apparent wind and true wind are represented by $V_{vw}$ and $V_{rw}$ , respectively, while the speed of the USV is denoted as $V_{usv}$ .

In the coordinate system illustrated in Figure 1, the heading angle of the USV on the map is represented by $\phi$ . The previous node serves as the starting point of the vector, while the current node acts as the endpoint. Consequently, the vector angle corresponds to the heading angle of the USV along that segment of the path. The motion state of the USV within the grid map is denoted as $(x_i,y_i,\phi _i)$ . The true wind is transformed from the inertial coordinate system to the attached coordinate system, resulting in the two components $u_{vw}$ and $v_{vw}$ of the apparent wind in the attached coordinate system, and the expression is as follows:

(4) \begin{equation} \left \{\begin{array}{llllllllll}{{u_{{\mathrm{rw}}}=-u+V_{{\scriptscriptstyle {rw}}}\cos (\theta _{{\mathrm{rw}}}-\varphi )}}\\ {{u_{vw}=-v+V_{{\scriptscriptstyle {rw}}}\cos (\theta _{{\mathrm{rw}}}-\varphi )}}\end{array}\right . \end{equation}

where $u$ and $v$ represent the USV’s speed in the attached coordinate system. $\theta _{rw}$ denotes the true wind direction.

It is important to note that $\theta _w$ is considered positive when the wind chord is oriented to the left. The magnitude of the apparent wind and the wind chord angle, $\theta _w$ , can be described as follows:

(5) \begin{align}&\qquad\qquad\qquad V_{vw}^2=u_{vw}^2+v_{vw}^2 \end{align}

(6) \begin{align}& \theta _w=\left \{\begin{array}{llllllllll}{{{\displaystyle \arctan \frac {V_{rw}}{u_{rw}}}}} & {{{\displaystyle u_{rw}\lt 0}}} \\ \\ {{{\displaystyle \arctan \frac {V_{rw}}{u_{rw}}-\pi \operatorname {sgn}(\nu _{rw})}}} & {{{\displaystyle u_{rw}\geq 0}}}\end{array}\right . \end{align}

2.3. Energy consumption modeling

In this paper, the energy consumption of the USV is categorized into two components: the energy required to overcome environmental resistance, specifically wind load, and the reaction force exerted by the water surface on the hull. The total energy consumption of the USV during travel can be determined using Eq. (7). To maximize the efficiency of energy use while ensuring task completion, it is essential to leverage natural forces that facilitate navigation, thereby effectively reducing overall energy consumption. Eqs. (8) and (9) indicate that the energy required for the USV to overcome water resistance is denoted as E1, while E2 represents the energy required to counteract wind resistance. The energy consumption data E1 for the subsequent drivable node at each position is stored in the adjacency matrix A1, whereas E2 is recorded in the adjacency matrix A2. The expressions are as follows:

(7) \begin{align}&\qquad\qquad\quad E=\alpha E_{1}+\beta E_{2} \end{align}

(8) \begin{align}& E_1=\int _0^tP_c(t)\;dt\qquad E_2=\int _0^tP_w(t)\;dt \end{align}

(9) \begin{align}&\qquad\qquad\,\,\, P(t)=F(t)V_{usv}(t) \end{align}

where $\alpha$ and $\beta$ represent the weighting parameter. $P_c$ and $P_w$ denote the power consumed by the USV to counteract water surface resistance and wind resistance, respectively. $F$ signifies the longitudinal force acting on the USV.

Figure 3.

Combined parameter analysis of energy consumption models.

Due to the gentle water current in the berthing area and the stable power output of the USV, it is assumed that the USV travels at a constant speed in this region, except when influenced by wind. Consequently, the frictional resistance between the water surface and the hull is treated as a constant value. As the USV maintains a constant speed, a shorter travel path results in a smaller value for E1. To reduce E2 in the presence of wind interference, it is essential to take into account the actual wind speed and direction, indicating that the shortest path may not be the optimal objective. However, for the computations to converge and enable the USV to reach its destination, it is necessary to employ the heuristic distance as a constraint. In the actual driving scenario, E1 is significantly larger than E2. To mitigate the dominance of E1 in the calculations, a weighted average approach is employed, incorporating weighted values and preprocessing of the data. This method enhances the differentiation of the fitness for each node, thereby benefiting the probabilistic search methods utilized by the algorithm. Both downwind and upwind paths are flexibly selected to optimize energy consumption, resulting in more reasonable routing. Furthermore, the value of $\alpha$ must be less than 0.5 to enhance the sensitivity of the USV to wind. Under these constraints, a one-factor experimental design method is implemented to investigate the impact of the parameter combinations on the experimental outcomes. As illustrated in Figure 3, Sailing Posture Changes and the Maximum Steering Angle of the USV are analyzed separately. A total of 20 samples are collected for averaging during periods when the USV can effectively harness the wind [Reference Gao, Li, Liu, Ge and Song40]. The experimental results indicate that with parameters $\alpha =0.35$ , $\beta =0.65$ the USV demonstrates increased sensitivity to wind and improved path smoothness.

3.1. Bacterial foraging optimization algorithm and function optimization

The BFO is a novel bionic class of algorithms proposed by K.M. Passino, inspired by the foraging behavior of E. coli. Bacteria tend to swim towards nutrient-rich areas and frequently flip when encountering unfavorable conditions. These behaviors are facilitated by flagella, which are integral to chemotaxis operation. The chemotaxis operation is central to the BFO, as it makes decisions based on comparisons, determining whether the bacteria should flip or swim. Therefore, the configuration of the chemotaxis operator is crucial for the effective convergence of the BFO. Upon the completion of the chemotaxis life cycle, when bacteria have sufficient nutrients, grow in size, and the environmental conditions are favorable, they undergo a process of binary fission, resulting in the formation of two new bacterial cells. This reproductive process not only regenerates a more energetic subset of the population but also tends to maintain the population size at an even number. Moreover, this propagation mechanism is linked to enhanced population quality and accelerated convergence performance among the bacteria. The elimination-dispersion operation is implemented in response to a sudden change in the environment or an attack on the colony that results in the death of some bacteria or their relocation to a new environment. This operation assesses the decision to relocate an individual by comparing a random number with the migration probability. Consequently, this operation enhances the algorithm’s capacity for global optimization [Reference Passino41].

The introduction of weight ratios in the energy consumption model alters the original fitness, which was solely based on distance. This modification encourages individuals to prioritize low-energy nodes, potentially leading them to deviate from the target node. When the E2 percentage is excessively high, USVs may continuously pursue the wind, resulting in slower convergence of computations or, in some cases, failure to converge altogether. Reward value serves as a feedback mechanism in the Markov decision-making process, where its magnitude is utilized to evaluate the quality of an intelligent agent’s action choices. To enhance result orientation, the reward function in reinforcement learning is introduced. In this context, the USV acts as the intelligent agent, and an energy consumption reward is devised based on the node where the USV is located after executing an action. If the distance from the USV’s current node to the target node is less than the distance from the previous node to the target node, a reward is granted; otherwise, a penalty is imposed. Consequently, the optimal path is identified as the one yielding the highest reward value. The incorporation of the reward function not only accelerates convergence but also ensures sample diversity. The reward function is formulated as follows:

(10) \begin{equation} r=\left \{\begin{array}{llllllllll}{{1/dist_i,\quad if dist_i\lt dist_{i-1}}}\\ {{-1/dist_i,\quad \quad otherwise}}\end{array}\right . \end{equation}

where $dist_i$ represents the distance from the USV’s current optional node to the target node, and $dist_{i-1}$ denotes the distance from the previous node to the target node.

Using a gradient strategy, among all the optional nodes in the subsequent step, the node that, after taking the action, has the shortest distance to the target node will receive the highest reward for this action. If $dist_i\lt dist_{i-1}$ , although this may not be the optimal choice, it will yield a relatively small reward value. The optimal path does not encourage the USV to travel in the opposite direction of the goal node; thus, if $dist_i\gt dist_{i-1}$ , this action will incur a penalty. Furthermore, in the context of global search, which also involves non-stationary rewards, receiving a smaller reward value when farther from the target node is advantageous for maintaining population diversity in the early stages of the algorithm, allowing for exploration of more unoccupied areas. As the USV approaches the target node and receives a larger reward value, this can enhance endpoint orientation and ultimately prompt the algorithm to converge.

Each step in this decision-making algorithm is linked to the fitness. Consequently, a new function for calculating the fitness is employed in the new application scenario. The fitness of the selectable nodes is derived from Eq. (11), and the initial paths are determined using a probabilistic selection method during the initialization phase. The expression is as follows:

(11) \begin{equation} fitness{{}}\left (i,j\right )=E(i,j)+\lambda r(i,j) \end{equation}

where $fitness$ is an indicator of individual merit. $\lambda$ denotes the equilibrium parameter. To ensure the proper functioning of the roulette wheel selection method during the initialization of the population, $\lambda$ is set to 20.

3.2. Trajectory constraints within the berth

While the natural environment near the berthing area is relatively favorable, uncertainties persist due to the movement of ships entering and leaving the harbor, as well as irregularities in the management of the berthing area. To mitigate the risk of accidents, the USV must reduce its traveling speed. In scenarios of insufficient energy, the power distribution between the propellers may not achieve optimal results, complicating challenging berthing maneuvers, such as continuous turning or executing large angle turns within a confined space. Consequently, it is essential to plan a simple and efficient berthing path.

When the USV approaches the berthing area, it is essential to impose constraints not only on energy consumption and path length but also on the change in the heading. Specifically, the steering angle of the USV must be kept below $45^\circ$ , which enhances the controller’s effectiveness and results in a smoother path. Upon reaching the designated docking point, the speed should ideally be reduced to zero. However, different routes will be planned based on varying practical situations within the berthing area. Consequently, the concept of equal proportional reduction, which is based on the number of nodes in the path under the same constraints, will ensure that the speed is reduced by a consistent amount across each segment of the path. The speed constraint of the USV can be expressed as follows:

(12) \begin{equation} V_{i}=\frac {N_{s}-N_{p i}}{N_{s}-N_{p k}}V_{k} \end{equation}

where $V_i$ denotes the speed of the USV upon reaching the $i-th$ node, $N_s$ represents the total number of nodes in the path, $N_{pi}$ indicates the number of nodes traveled so far, $N_{pk}$ signifies the number of nodes traversed upon entering the berthing area, and $V_k$ denotes the speed of the USV when it enters the berthing area.

Due to the constraints of heading and speed, a portion of the trajectory of the USV as it enters the berthing area from the upper left corner is illustrated in Figure 4. Both trajectories consist of four nodes and three path segments, with the velocity required to be reduced by one-third in each segment of the path.

Figure 4.

Bow restraint at berth.

The influence of natural factors on path planning is significant, particularly regarding the berthing capabilities of the USV in environments with obstacles. The effectiveness of the USV’s berthing ability is primarily determined by its obstacle perception and avoidance capabilities. By optimizing relevant parameters and imposing constraints on speed and heading, the proposed method enables the USV to smoothly berth at a designated position along a continuously curving path while successfully avoiding obstacles. This approach mitigates the risks associated with sharp turns or erratic pathways, thereby enhancing the stability and comfort of the USV during the berthing process. Furthermore, the constraints on speed and heading ensure that the USV adheres to a predetermined safe trajectory while berthing, maintaining a stable sailing posture and preventing swaying or course deviations caused by improper operation or external disturbances.

4.1. Generating paths

In the grid map, USV icon denotes the starting point, and USV dashed icon represents the endpoint. The red irregular objects symbolize obstacles, while the size and direction of the blue arrows indicate wind speed and direction, respectively. Figure 5 illustrates the lowest energy path, where the USV aims for the endpoint while searching for a route that leverages favorable winds, thereby conserving its energy output. Although the actual wind conditions are complex, a new fitness function can be assigned to each segment of the path to facilitate probabilistic selection. In Figure 5, the front part of the right-side route experiences a favorable wind direction, yet the wind speed is noticeably lower than that of the left-side route. Additionally, there exists a headwind section in the middle of the right-side route. Individuals often encounter challenges in selecting intersections, leading to increased confusion, and the initial population may include routes on the right side that are high in energy. However, after applying three operators and several constraints, optimal or relatively optimal individuals will emerge in the final population. The optimal path, referred to as MC-BFO, depicted in Figure 6, is determined by using the frequency of changes in the USV’s sailing posture as a filtering criterion within the final population. The pseudo-code for MC-BFO is provided below. The energy consumption value for the path illustrated in Figure 6 is represented by a blue pentagram in Figure 8. Although this path does not exhibit the lowest energy consumption overall, it is relatively more efficient. Additionally, this path is smoother than the one with the lowest energy consumption, making it more suitable for the safe travel of the USV, particularly in scenarios involving limited energy and reduced control effectiveness. The path’s design is therefore more reasonable.

Figure 5.

Minimum energy consumption berthing path.

Figure 6.

Considering the smoothness of the berthing path.

A comparison of the data between the paths shown in Figure 5 and Figure 6 is presented in Table I. Figure 7 depicts the shortest path, and it is calculated that the path in Figure 5 consumes 40.49% less energy than that in Figure 7, while the path in Figure 6 consumes 36.50% less energy than the path in Figure 7.

Table I.

Comparison of data between the optimal path and the lowest energy path.

Figure 7.

Shortest berthing path.

Figure 8.

Comparison of individual energy consumption of optimal and shortest paths.

Algorithm 1

A bacterial foraging optimization algorithm with multiple constraints (MC-BFO)

4.2. Experimental comparison

To verify the effectiveness of the proposed algorithm, simulation experiments were conducted using MATLAB 2016a (Windows 10, Intel(1R) Core (TM) i5–8300H CPU @ 2.30 GHz 2.30 GHz with 16G of RAM). Three experimental environments were established with grid map sizes of 20 $\times$ 20, 25 $\times$ 25, and 30 $\times$ 30. The ratio of obstacle area to map area remains consistent across these environments; however, the distribution of obstacles varies among them, as does the distribution of wind speed and direction. The paths planned by the MC-BFO, BFO, and GA are compared and analyzed under both downwind and upwind conditions.

Figure 9 illustrates the optimal paths generated by the three algorithms across three environments under both downwind and upwind conditions (with start and end positions exchanged). In this figure, the red solid line represents the path generated by the MC-BFO, the gold dashed line corresponds to the BFO, and the black dotted line indicates the GA. Table II presents data on the average number of sailing posture (ASP), average energy consumption (AEC), and average path length (APL) for the three algorithms in the three environments. These metrics were obtained from 20 repetitions of the experiment conducted under both downwind and upwind conditions. The results indicate that the MC-BFO path exhibits fewer sailing posture changes across all three environments, resulting in a simpler and more efficient trajectory. Specifically, under downwind conditions in Environment 1, the MC-BFO path demonstrates a reduction of 28.77% in AEC compared to the BFO path and 45.75% compared to the GA path. In upwind conditions, the MC-BFO path achieves a 30.93% reduction in AEC relative to the BFO path and 38.29% compared to the GA path. In Environment 2, the MC-BFO path saves 43.81% in AEC when compared to the BFO path and 44.56% relative to the GA path under downwind conditions, while it saves 8.07% over the BFO path and 5.69% over the GA path in upwind conditions. In Environment 3, the MC-BFO path results in a 58.73% reduction in AEC compared to the BFO path and 64.09% compared to the GA path under downwind conditions, and a 45.29% reduction over the BFO path and 45.95% over the GA path under upwind conditions. These findings suggest that the USV conserves more energy during downwind travel compared to upwind travel. While the analysis of energy consumption is intricately linked to the energy model, the parameters of the USV itself, and the wind conditions, the MC-BFO path consistently consumes the least energy under identical map models and parameter settings. Regarding APL, the differences among the paths range from 0.3 to 8.01%, with the APL of the MC-BFO path being shorter, smoother, and more reasonable in certain areas with complex obstacles. This characteristic is advantageous for the USV’s safe navigation and smooth berthing, particularly during periods of energy scarcity.

Figure 9.

Optimal berthing paths for different algorithms in each environment.

Table II.

Comparison of performance parameters of various algorithms.