AutoTuning Environment for Static Obstacle Avoidance Methods Applied to USVs

Abstract

This work is focused on reactive Static Obstacle Avoidance (SOA) methods used to increase the autonomy of Unmanned Surface Vehicles (USVs). Currently, there are multiple approaches to avoid obstacles, which can be applied to different types of USV. In order to assist in the choice of the SOA method for a particular vessel and to accelerate the pretuning process necessary for its implementation, this paper proposes a new AutoTuning Environment for Static Obstacle Avoidance (ATESOA) methods applied to USVs. In this environment, a new simplified modelling of a LIDAR (Laser Imaging Detection and Ranging) sensor is proposed based on numerical simulations. This sensor model provides a realistic environment for the tuning of SOA methods that, due to its low load computation, is used by evolutionary algorithms for the autotuning. In order to analyze the proposed ATESOA, three SOA methods were adapted and implemented to consider the measurements given by the LIDAR model. Furthermore, a mathematical model is proposed and evaluated for using as USV in the simulation enviroment. The results obtained in numerical simulations show how the new ATESOA is able to adjust the SOA methods in scenarios with different obstacle distributions.

1. Introduction

Nowadays, the development of the USV is an active and growing field of research. The reason is that the applications of a USV are very broad, covering from environmental control, as well as many scientific and commercial applications, to national security and surveillance issues [1]. The competencies required a vessel in order to be considered autonomous can be grouped into three fields: Navigation, Guidance and Control, [1,2,3]. In more detail, these three fields can be decomposed into the subsystems shown in Figure 1. Where, firstly, the state vector of the vessel dynamics must be estimated from state observers and wave filters [4,5]. Once the state of the vehicle is known, the detection system processes the information received by the surrounding sensors in order to generate a model of the environment in which the USV is located [6,7,8,9,10]. Using this model of the environment, the guidance system (composed by algorithms of: obstacle avoidance [8,11,12,13,14,15,16], path following [4,5,17,18] and path planning [3,7,17,19,20]) demands the course and speed setpoints which guide the vessel safely to its goal. With the aim of ensuring that the vehicle reaches these setpoints, the control system commands the vessel’s actuators (propulsion system and steering machine) [4,5]. In this way, the degree of autonomy of the USV is determined by the level of development of each subsystem, as well as by the integration of the set [1,2,3]. From the subsystems shown in Figure 1, this work is focused on the obstacle avoidance systems applied to USVs [8,10,12,13,14,15,16,21,22,23,24,25,26,27,28,29,30,31]. Most of this work is focussed on the avoidance of dynamic obstacles in open sea situations [8,12,13,14,15,21,22,23,24,25,26,27,28]. And, although some of them also consider a specific approach to static obstacles [15,22,23,25], their study and evaluation is secondary. However, on a lesser scale, works such as [10,16,17,29,30,31] are centered on the capacity of a USV to avoid static obstacles.

All these obstacle avoidance methods have a common requirement, they need a simulation environment for their design, development and evaluation. In this aspect, there are very complete approaches [32,33,34,35], where the authors are focused on the development of simulation environments for the validation of USVs under realistic sea conditions (wind, current and waves), considering also the models of inertial and LIDARs sensors available in the GAZEBO framework [36]. On the other hand, most of the works focused on obstacle avoidance methods do not take into account the ambient disturbances that affect the USV and, moreover, they assume known: position, dimension and speed of the obstacles. As exceptions, in [12,15] the error present in the estimates of the environment model is taken into account. In addition, in [15,22] the external disturbances that affect the USV are also considered by modelling them as constant forces and moments referred to the earth axes. As an alternative to these extremely simplified simulation environments [10,13,21,23,24,25,26,27,28,29,30,31], or the highly complete environments proposed in [32,33,34,35,36], in this work a new simulation environment based on the simplified modelling of a LIDAR sensor is proposed. As a result, the time of numerical simulations is reduced and the methods of obstacle avoidance are exposed to a higher level of uncertainty present in the environment surrounding the vessel. This facilitates that the simulation environment proposed in this work can be used by iterative optimization algorithms, such as the evolutionary algorithms [37], with the advantage that it requires less computer load. In this way, the main contribution of this work is achieved, a new autotuning environment for SOA methods applied to USVs. Currently, there are different approaches for the autotuning of unmanned vehicles [38,39,40,41,42,43,44]. In [38] the authors use machine learning for autotuning the guidance system on a physical ground vehicle. As autotuning methods applied in the marine environment, in [39] a genetic algorithm is used to autotune the parameters of a sliding mode controller, which is applied to the mathematical model of an oil tanker vessel. On the other hand, the work carried out in [40] is focused on the autotuning of a dynamic positioning system in response to changes in the sea state. This is evaluated in numerical simulations with the dynamic model of an offshore supply vessel. With regard to underwater vehicles, in [41] the authors propose a autotuning method for heading control on an underwater vehicle and validate their approach with the micro-ROV vehicle. In relation to these works, the objective of the ATESOA approach proposed in this paper is to facilitate and accelerate the pretuning process necessary to implement SOA methods in USVs. In addition, in order to evaluate this autotuning environment, three methods of static obstacle avoidance have been adapted and implemented in this work [10,23,29,30], which have already been used in USVs. As a secondary contribution of this work, through the application of the generalized potential fields [45,46], the performance of potential fields method for USVs proposed in [23,30] has been improved. Finally, a USV model must be used to evaluate the new autotuning environment. In the current state of the art different mathematical models of USVs have been proposed [22,26,28,47,48,49], most of them for vessels of small length. This is mainly because one of the advantages of this kind of vehicles is that, as these do not need a crew, the size of the vehicle can be reduced, increasing its manoeuvrability and reducing costs [1,2]. Many of these models [22,26,28,47,48] simplify the dynamics of the vessel since they do not take into account the modelling of the actuators and the effect of the ocean current. In order to make realistic the simulation environment proposed in this work, a mathematical model is proposed for an USV of 9.2 m of length. This model has been synthesized from previous work in the field of marine surface vehicle control [4,5,28,49,50,51,52], whose parameters have been calculated through numerical simulations in order to obtain a characteristic behavior of marine surface vehicles [5,53,54,55,56,57]. Concretely, the USV model is composed of: model of the vessel and the actuators, current effect and course/speed controllers.

Major Contribution

This work is focused on the reactive obstacle avoidance methods, which form part of USV guidance systems. The main contribution of this work is a autotuning environment for SOA methods applied to USVs, the ATESOA environment. ATESOA can be applied to different SOA methods and vessels. In particular, it is based on the simplified modelling of a LIDAR sensor. This model takes into account the uncertainty present in the environment by considering the perspective effect of the sensor. In addition, a variation of the error present in the measurements delivered is also modelled as a function of the distance to the obstacle. Given its low computational load, this sensor model is used to perform a realistic autotuning of SOA methods using an evolutionary algorithm. In order to analyze the proposed ATESOA, the SOA methods: LROABRA (Local Reactive Obstacle Avoidance Based on Region Analysis), Vector Field Histogram+ and Potential Fields, have been implemented and adapted to the measures delivered by the sensor model. The results obtained have been analysed in detail in terms of performance and robustness.

The organization of this paper is as follows: Section 2 contains the mathematical models that form the new autotuning environment. The Section 3 presents the SOA methods adapted and implemented in this work to evaluate the autotuning environment. In Section 4, the setting of the genetic algorithm used for autotuning is presented, as well as the indicators that have been defined to evaluate the performance of each avoidance method. Section 5 contains the results obtained in numerical simulations. Lastly, in Section 6 the Conclusions are presented.

2. Autotuning Environment Based on Numerical Simulations

This section describes the mathematical models that form the simulation environment proposed in this work for the autotuning of SOA methods applied to USVs.

2.1. USV Model

The USV model is based on the development proposed in [4,5] for the modelling of three degrees of freedom applied to vessels. Specifically, the non-linear modelling proposed by Norrbin [5,49] is taken into account for the damping and the effect of the ocean currents is incorporated into the model in terms of relative velocities, considering that the current variation is significantly slower than the ship dynamics, [5]. Thus, the dynamic model of the vessel is defined as:

$$\begin{array}{cc}\hfill & \dot{u}(m-X_{\dot{u}})=mvr+m{x}_g{r}^2-Y_{\dot{v}}{v}_{r}r+X_u{u}_r+X_{\left|u\right|u}\left|u_r\right|u_r+{\tau }_x\hfill \\ \hfill & \dot{v}(m-Y_{\dot{v}})+m{x}_g\dot{r}=-mur-X_{\dot{u}}{u}_{r}r+Y_v{v}_r+Y_{r}r+Y_{\left|v\right|v}{v}_r|v_r|+Y_{\left|r\right|v}\left|r\right|v_r+Y_{\left|v\right|r}\left|v_r\right|r+{\tau }_y\hfill \\ \hfill & \dot{r}(I_z-N_{\dot{r}})+m{x}_g\dot{v}=-m{x}_{g}ur-N_{uv}{v}_r{u}_r+N_{r}r+N_v{v}_r+N_{\left|v\right|r}|v_r|r+N_{\left|r\right|v}\left|r\right|v_r+N_{\left|r\right|r}\left|r\right|r+{\tau }_n\hfill \\ \hfill & \underset{\dot{\mathit{\eta }}}{\underbrace{\left(\begin{array}{c}{\dot{x}}_E\\ {\dot{y}}_E\\ \dot{\psi }\end{array}\right)}}=\underset{\mathit{R}\left(\psi \right)}{\underbrace{\left(\begin{array}{ccc}cos\left(\psi \right)& -sin\left(\psi \right)& 0\\ sin\left(\psi \right)& cos\left(\psi \right)& 0\\ 0& 0& 1\end{array}\right)}}\underset{\mathit{\nu }}{\underbrace{\left(\begin{array}{c}u\\ v\\ r\end{array}\right)}}\hfill \end{array}$$

(1)

where $\mathit{\eta }$ represents the position vector of the vessel referenced to earth axes, $\mathit{\nu }$ and ${\mathit{\nu }}_{\mathit{r}}$ are the vectors of velocity and relative velocity referenced to body axes, $\mathit{R}\left(\psi \right)$ is the rotational matrix from Body-axes reference system to Earth-axes reference system (see Figure 2), ${\mathit{\tau }}_{\mathit{act}}={\left[{\tau }_x{\tau }_y{\tau }_z\right]}^T$ represents the forces and moment applied by the actuators and the viscous non-linear term $N_{uv}$ replaces $(X_{\dot{u}}-Y_{\dot{v}})$, which causes the destabilizing moment or Munk moment [5]. The relative velocity vector of the Model (1) is defined according to Equation (2), where $V_c$ is the module of the current velocity and ${\beta }_c$ its direction.

$${\mathit{\nu }}_{\mathit{r}}={\left[u_r{v}_{r}r\right]}^T={\left[\left(u-V_{c}cos({\beta }_c-\psi )\right)\left(v-V_{c}sin({\beta }_c-\psi )\right)r\right]}^T$$

(2)

On the other hand, the actuators are the elements that allow the steering of the vessel [4,5,49,50,51]. For the USV model proposed in this work, an outboard motor is considered as the actuator, which generates the propulsion force and the steering torque. This actuator is modelled as the set formed by propeller and rudder [49,50,51]. In this way, the force vector produced by the outboard located on the vessel’s centreline is defined according to Equation (3).

$${\mathit{\tau }}_{\mathit{act}}=\left(\begin{array}{c}{\tau }_x\\ {\tau }_y\\ {\tau }_n\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\\ 0& -l_x\end{array}\right)\left(\begin{array}{c}{T}_{act}-D_{act}\\ L_{act}\end{array}\right),\begin{array}{cc}\hfill & T_{act}=T_{\left|n\right|n}\left|n\right|n-T_{\left|n\right|u}\left|n\right|u_r\hfill \\ \hfill & D_{act}=D_{\left|\delta \right|}\left|\delta \right||u_r|u_r\hfill \\ \hfill & L_{act}=(L_{\delta }\delta -L_{\left|\delta \right|\delta }\left|\delta \right|\delta \left)\right|u_r|u_r\hfill \end{array}$$

(3)

where $l_x$ is the distance of the actuator to the vessel’s center of gravity, n represents the revolutions per minute (RPM) of the propeller, $\delta$ is the rudder angle relative to the vessel’s centerline, $T_{act}$ represents the thrust generated by the propeller and, finally, $D_{act}$ and $L_{act}$ are the drag and transversal forces generated by the rudder, respectively. In addition, the dynamics of the outboard motor is taken into account as:

$$\begin{array}{cc}\hfill & \dot{\delta }=(1/{\tau }_{\delta })({\delta }_c-\delta ),\delta \in [-{\delta }_{lim},{\delta }_{lim}],\dot{\delta }\in [-{\dot{\delta }}_{lim},{\dot{\delta }}_{lim}]\hfill \\ \hfill & \dot{n}=(1/{\tau }_n)(n_c-n),n\in [n_{min},n_{max}]\hfill \end{array}$$

(4)

where ${\delta }_c$ and $n_c$ represent the setpoints for the rudder angle and the propeller RPM, respectively. The time constants ${\tau }_{\delta }$ and ${\tau }_n$ define the dynamics of the outboard motor.

Once the mathematical model has been defined, in order to adjust its parameters, the turning radius $\left(K_r=\sqrt{Iz/m}\right)$ and the Nomoto model identified in [28] for the vessel Viknes 830, as well as its maximum speed, have been taken as references to synthesize the proposed mathematical model and its parameter values. The main features of Viknes 830 are listed in

Table 1. Viknes 830 parameters expressed in the international system of units [28].

Symbol	Description	Value	Symbol	Description	Value
$B_v$	Beam	2.97	$m_v$	Mass	3980
$L_v$	Length	8.52	$I_{z_v}$	z-axis moment of inertia	19,703
$D_v$	Draft	0.82	$u_{max}$	Maximum forward speed	10.5
$K_{Nomoto-Model}$	Gain $r/\delta$	0.5	${\tau }_{Nomoto-Model}$	Time constant	4
$r_{max}$	Maximum angular rate	0.35	${\delta }_{max}$	Maximum rudder angle	0.26

.

The parameters of the USV model proposed in this paper are collected in

Table 2. Parameters of the vessel model proposed in this work.

Symbol	Value	Symbol	Value	Symbol	Value	Symbol	Value
L	9.2 m	B	3 m	D	0.7 m	$x_g$	0.5 m
$l_x$	4.5 m	m	3705.78 kg	$I_z$	18,345.50 kg m${}^2$	$X_{\dot{u}}$	−185.23 kg
$Y_{\dot{v}}$	−7272.34 kg	$N_{\dot{r}}$	−51,433.21 kg m${}^2$	$\rho$	1027 kg/m${}^3$	$X_u$	−512.5
$Y_v$	−2200	$Y_r$	−500	$N_r$	−15,000	$N_v$	−311.1
$X_{\left|u\right|u}$	−59.5	$Y_{\left|v\right|v}$	−5000	$Y_{\left|r\right|v}$	−7000	$N_{\left|r\right|r}$	−4000
$N_{\left|v\right|r}$	−50	$Y_{\left|v\right|r}$	−40	$N_{\left|r\right|v}$	−1700	$N_{uv}$	−22
$T_{\left|n\right|n}$	1.2 $\times {10}^{-3}$	$T_{\left|n\right|u}$	3.27 $\times {10}^{-2}$	$L_{\delta }$	1.3886	$L_{\left|\delta \right|\delta }$	1.28 $\times {10}^{-2}$
$D_{\left|\delta \right|}$	0.71	${\tau }_n$	0.5 s	${\tau }_{\delta }$	1 s	$n_{max}$	3300 RPM
$n_{min}$	800 RPM	${\delta }_{lim}$	35$°$	${\dot{\delta }}_{lim}$	10$°$/s

. Assuming a homogeneous mass distribution and taking as a reference the Viknes 830, the mass of the new model and its moment of inertia are obtained according to the Equation (5). In addition, the added masses that affect the vessel are estimated according to [52], see Equation (6).

$$m=\left({\displaystyle \frac{LBD}{L_v{B}_v{D}_v}}\right)m_v,I_z=K_{r_v}^{2}m$$

(5)

where L, B y D represent the length, the beam and the draft of the vessel, respectively.

$$X_{\dot{u}}=-0.05m,Y_{\dot{v}}=-(1/2)\rho \pi D^{2}L,N_{\dot{r}}=-(1/24)(0.1m{B}^2+\rho \pi D^2{L}^3)$$

(6)

where $\rho$ defines the density of seawater.

The other parameters are adjusted to obtain stationary values for forward speed and angular rate close to those shown in Table 1. In order to verify this setting, a zero current velocity is assumed ($V_c=0$) and, once the operating range of the actuators ($n_c$ and ${\delta }_c$) has been discretised, the points of equilibrium of the system are obtained (stationary values of the velocity vector $\nu$). As can be seen in Figure 3, the forward speed increases with propoller revolution (where $u_{max}=10.5$ m/s) and falls when the rudder is engaged. Furthermore, the relation of the stationary value of r with the rudder angle is consistent with the behaviour of surface marine vehicles exposed in [5,53,54,55,56].

On the other hand, the parameters listed in Table 2 must ensure a stable behavior of the Model (1) at all its operating points. In order to check its dynamic behavior, Model (1) is linearized at each operating point shown in Figure 3. As a result, the effective time constants obtained for the velocity vector $\nu$ are shown in Figure 4. As can be seen, the system is stable over its entire operating range and its dynamics becomes faster in function of the vessel’s speed [5,53,54,55,56]. Moreover, with a zero rudder angle, the effective time constant ${\tau }_r$ varies, depending of the speed, around the 4 s obtained by [28] in the Nomoto model identified for the Viknes 830. To conclude with the analysis of the Model (1), it is evaluated through numerical simulations by carrying out different sea tests [57]. The Runge-Kutta numerical integration method of fourth order is used with an integration step of 0.01 s, holding $n_c=2430$ RPM and $V_c=0$ m/s. The Turning Circle Tests can be seen in Figure 5a. On the other hand, in Figure 5b,c the spiral (or Dioudonne) manoeuvre and zig-zag manoeuvres are shown.

2.2. Course and Speed Controllers

Most obstacle avoidance systems for USVs [10,12,13,15,21,22,23,24,25,26,27,29,30,31] require that the vessel be equipped with course and speed controllers. These controllers guarantee, within a range of defined sea states, that the vessel reaches the course and speed requested by the avoidance system in a finite time. To this end, there are currently a multitude of control methods [4,5,52,53,58]. As the objective of this work is not focused on the design of the controllers, and in order to limit the problem, two PID structures are used to control the course and speed, see Equation (7). Moreover, due to the strong non-linearity of the Model (1), the adaptive control technique Gain Scheduling is used [53,58]. For this purpose, the controller parameters are set experimentally for three operating points: low ($u_{low}=3$ m/s), medium ($u_{medium}=6$ m/s) and high ($u_{high}=9$ m/s) speed, see

Table 3. Parameters of the course and speed controllers.

	${\mathit{K}}_{{\mathit{P}}_{\mathit{\psi }}}$	${\mathit{K}}_{{\mathit{I}}_{\mathit{\psi }}}$	${\mathit{K}}_{{\mathit{D}}_{\mathit{\psi }}}$	${\mathit{K}}_{{\mathit{P}}_{\mathit{u}}}$	${\mathit{K}}_{{\mathit{I}}_{\mathit{u}}}$	${\mathit{K}}_{{\mathit{D}}_{\mathit{u}}}$
$u_{low}$	−40	−0.01	−30	40	2	6
$u_{medium}$	−16	−0.01	−32	35	2	12
$u_{high}$	−9.5	−0.01	−17	45	2	18

. For intermediate speeds, linear interpolations are carried out between the PIDs. Finally, these controllers have been discretised using the Euler backwards method with a sampling period $T_m=0.05$ s.

$$\begin{array}{cc}\hfill {\delta }_c\left(k\right)& =\underset{{\delta }_P\left(k\right)}{\underbrace{K_{P_{\psi }}{e}_{\psi }\left(k\right)}}+\underset{{\delta }_I\left(k\right)}{\underbrace{T_m{K}_{I_{\psi }}{e}_{\psi }\left(k\right)+{\delta }_I(k-1)}}+\underset{{\delta }_D\left(k\right)}{\underbrace{{\displaystyle \frac{K_{D_{\psi }}}{T_m}}\left(e_{f_{\psi }}\left(k\right)-e_{f_{\psi }}(k-1)\right)}}\hfill \\ \hfill n_c\left(k\right)& =\underset{n_P\left(k\right)}{\underbrace{K_{P_u}{e}_u\left(k\right)}}+\underset{n_I\left(k\right)}{\underbrace{T_m{K}_{I_u}{e}_u\left(k\right)+n_I(k-1)}}+\underset{n_D\left(k\right)}{\underbrace{{\displaystyle \frac{K_{D_u}}{T_m}}\left(e_{f_u}\left(k\right)-e_{f_u}(k-1)\right)}}\hfill \\ \hfill e_{f_i}\left(k\right)& ={\displaystyle \frac{1}{c_f+1}}\left(e_i\left(k\right)+c_f{e}_{f_i}(k-1)\right),e_i\left(k\right)=s{p}_i\left(k\right)-i\left(k\right)\hfill \end{array}$$

(7)

where $i\in \{u,\psi \}$, $c_f={\tau }_f/T_m$ and the time constants have been fixed as ${\tau }_f=0.1(K_D/K_P)$.

Due to the effect of external disturbances, a compensation is carried out on the error signals of the controllers (7) [5,52]. In this way, the course angle ($\chi =\psi +\beta$, see Figure 2) is controlled instead of the heading angle ($\psi$), as well as the module of the velocity vector (U) of the USV.

$$\begin{array}{cc}\hfill e_{\psi }\left(k\right)& =s{p}_{\psi }\left(k\right)-\chi \left(k\right),\chi \left(k\right)=\psi \left(k\right)+atan2(v\left(k\right),u\left(k\right))\hfill \\ \hfill e_u\left(k\right)& =s{p}_u\left(k\right)-U\left(k\right),U\left(k\right)=\sqrt{u{\left(k\right)}^2+v{\left(k\right)}^2}\hfill \end{array}$$

(8)

2.3. Simplified Modelling of a LIDAR Sensor

As a novel contribution made in this work, and with the aim of increasing the realism of the numerical simulations carried out, a LIDAR sensor has been modelled. As a reference, the characteristics of the LIDAR Ultra Puck of the company Velodyne have been used [59]. In addition, as in [10,12,13,15,21,22,23,24,25,26,27,28,29,30,31], a 2D Cartesian space ($x_E$, $y_E$) is used to evaluate the obstacle avoidance systems applied to USVs. Therefore, only one of the vertical channels of the LIDAR sensor is taken into account. Thus, the computation of the sensor model is decreased, reducing simulation times and allowing its use in iterative optimization algorithms, as evolutionary algorithms. In particular, from the information provided by the manufacturer in [59], the following data are taken:

• Distance Measurement Range of 200 m ($d_{range}=200$ m).
• Horizontal Field of View of 360$°$ ($h_{range}=360°$).
• Horizontal Angular Resolution between 0.1$°$ and 0.4$°$. Looking for the worst case, $h_{res}=0.4°$.
• Rotation Rate between 5 Hz to 20 Hz. Looking for a balance between computation and refresh rate of the environment information $\left(T_m^L\right)$, a rotation frequency $f_m^L=10$ Hz is set.
• Range Accuracy up to $\pm 3$ cm. The minimum value for the standard deviation of the WGN (White Gaussian Noise) that affects the measurement is taken as ${\sigma }_{min}^L=0.03$ m.

Before proceeding with the modelling of the LIDAR sensor, it is necessary to define the obstacle scenarios in which the USV is located. In this work, each obstacle $Obs{t}_E$ is defined as a polygon formed by four concatenated segments, which are referenced to Earth-axes. Consequently, a scenario ($S_E$) is the result of the union of a finite number of obstacles ($n_O$), see Equation (9).

$$S_E=\bigcup _{i_O=1}^{n_O}Obs{t}_E^{i_O},Obs{t}_E^{i_O}=\left\{\begin{array}{c}se{g}_{1E}^{i_O}=\left(x_{1E}^{i_O},y_{1E}^{i_O},x_{2E}^{i_O},y_{2E}^{i_O}\right),se{g}_{2E}^{i_O}=\left(x_{2E}^{i_O},y_{2E}^{i_O},x_{3E}^{i_O},y_{3E}^{i_O}\right),\hfill \\ se{g}_{3E}^{i_O}=\left(x_{3E}^{i_O},y_{3E}^{i_O},x_{4E}^{i_O},y_{4E}^{i_O}\right),se{g}_{4E}^{i_O}=\left(x_{4E}^{i_O},y_{4E}^{i_O},x_{1E}^{i_O},y_{1E}^{i_O}\right)\hfill \end{array}\right\}$$

(9)

where $se{g}_{1E\cdots 4E}^{i_O}$ represents the concatenated segments that form each obstacle and $x_{1E\cdots 4E}^{i_O}$ are the coordinates of its four points.

It is important to highlight that we have chosen to model each obstacle as a rectangle due to the facility for generating scenarios from these polygons. However, this LIDAR sensor model can also be used in the presence of irregular polygons such as the ones that define the marine environments [16,17,19]. For this purpose, two approaches can be taken. Firstly, Delaunay Triangulation [60] can be applied to polygons that form the coastlines [61]. This triangulation decomposes the polygon into a finite number of triangles ($n_O$). Thus, each triangle is defined as a polygon formed by three concatenated segments. Secondly, the scenario can be defined as the set of segments collected in .shp files such as the ones collected in [61].

Defined the LIDAR parameters, as well as the obstacle scenario ($S_E$), the measurement of this sensor can be modelled as the set of segments $L_B$ (10), which is referenced to Body-axes.

$$\begin{array}{c}{L}_B=\left\{se{g}_B^{i_L}\left(x_{1B}^{i_L},y_{1B}^{i_L},x_{2B}^{i_L},y_{2B}^{i_L}\right)|x_{1B}^{i_L}=y_{1B}^{i_L}=0,x_{2B}^{i_L}=cos\left(i_L{h}_{res}\right)d_{range},y_{2B}^{i_L}=sin\left(i_L{h}_{res}\right)d_{range}\right\}\hfill \\ I_L=\left\{i_L|i_L\in \mathbb{N},i_L<n_L\right\},n_L\in \mathbb{N}\hfill \end{array}$$

(10)

where $x_{1B,2B}^{i_L}$ and $y_{1B,2B}^{i_L}$ represent the coordinates of the two points that form a segment and $n_L$ is the number of segments (LIDAR beams) that compose the set $L_B$ ($n_L=h_{range}/h_{res}=900$).

In order to model the distances measured by the LIDAR in its horizontal field of view ($h_{range}$), the points of intersection between the two sets of segments $S_E$ and $L_B$ are taken. For this purpose, both sets of segments must refer to the same reference system. Therefore, once the set $L_B$ has been expressed in matrix form, it is rotated and translated to Earth-axes at each simulation instant (k):

$$\begin{array}{cc}\hfill & \underset{{\mathit{L}}_{\mathit{E}}}{\underbrace{\left(\begin{array}{cccc}{x}_{1E}^0& y_{1E}^0& x_{2E}^0& y_{2E}^0\\ ⋮& ⋮& ⋮& ⋮\\ x_{1E}^{i_L}& y_{1E}^{i_L}& x_{2E}^{i_L}& y_{2E}^{i_L}\\ ⋮& ⋮& ⋮& ⋮\\ x_{1E}^{n_L-1}& y_{1E}^{n_L-1}& x_{2E}^{n_L-1}& y_{2E}^{n_L-1}\end{array}\right)}}=\underset{{\mathit{L}}_{\mathit{B}}}{\underbrace{\left(\begin{array}{cccc}{x}_{1B}^0& y_{1B}^0& x_{2B}^0& y_{2B}^0\\ ⋮& ⋮& ⋮& ⋮\\ x_{1B}^{i_L}& y_{1B}^{i_L}& x_{2B}^{i_L}& y_{2B}^{i_L}\\ ⋮& ⋮& ⋮& ⋮\\ x_{1B}^{n_L-1}& y_{1B}^{n_L-1}& x_{2B}^{n_L-1}& y_{2B}^{n_L-1}\end{array}\right)}}\left(\begin{array}{cc}{\mathit{R}}_{\mathit{p}}^T& \mathit{0}\\ \mathit{0}& {\mathit{R}}_{\mathit{p}}^T\end{array}\right)+\hfill \\ \hfill & +{\left(\begin{array}{ccccc}1& \cdots & 1& \cdots & 1\end{array}\right)}^T\left(\begin{array}{cccc}{x}_E\left(k\right)& y_E\left(k\right)& x_E\left(k\right)& y_E\left(k\right)\end{array}\right),{\mathit{R}}_{\mathit{p}}=\left(\begin{array}{cc}cos\left(\psi \left(k\right)\right)& -sin\left(\psi \left(k\right)\right)\\ sin\left(\psi \left(k\right)\right)& cos\left(\psi \left(k\right)\right)\end{array}\right)\hfill \end{array}$$

(11)

The set of segments that form the scenario is also expressed in matrix form $S_E$:

$${\mathit{S}}_{\mathit{E}}=\left(\begin{array}{c}{\mathit{Obst}}_{\mathit{E}}^{\mathit{1}}\\ ⋮\\ {\mathit{Obst}}_{\mathit{E}}^{{\mathit{i}}_{\mathit{O}}}\\ ⋮\\ {\mathit{Obst}}_{\mathit{E}}^{{\mathit{n}}_{\mathit{O}}}\end{array}\right),{\mathit{Obst}}_{\mathit{E}}^{{\mathit{i}}_{\mathit{O}}}=\left(\begin{array}{c}se{g}_{1E}^{i_O}\\ se{g}_{2E}^{i_O}\\ se{g}_{3E}^{i_O}\\ se{g}_{4E}^{i_O}\end{array}\right)=\left(\begin{array}{cccc}{x}_{1E}^{i_O}& y_{1E}^{i_O}& x_{2E}^{i_O}& y_{2E}^{i_O}\\ x_{2E}^{i_O}& y_{2E}^{i_O}& x_{3E}^{i_O}& y_{3E}^{i_O}\\ x_{3E}^{i_O}& y_{3E}^{i_O}& x_{4E}^{i_O}& y_{4E}^{i_O}\\ x_{4E}^{i_O}& y_{4E}^{i_O}& x_{1E}^{i_O}& y_{1E}^{i_O}\end{array}\right)$$

(12)

As the problem of the intersection of segments can be solved from the field of computational geometry [60], in this work the matrices (11) and (12) are used as inputs for the function proposed in [62]. This function calculates the intersections between two sets of segments, considering coincident and parallel segments. Between other outputs, it returns the matrices ${\mathit{F}}_{\mathit{Inter}}$ and ${\mathit{ND}}_{{\mathit{L}}_{\mathit{E}}}$, both of dimension $n_L\times n_O$. Where ${\mathit{F}}_{\mathit{Inter}}$ defines which segments of the sets (11) and (12) are intersected (${\mathit{F}}_{\mathit{Inter}}(i_L,i_O)=1$) and ${\mathit{ND}}_{{\mathit{L}}_{\mathit{E}}}$ returns the normalized distances between the intersection points and the starting points of the segments of the set $L_E$. From the matrix ${\mathit{F}}_{\mathit{Inter}}$, in the Equation (13) the vector ${\mathit{f}}_{\mathit{Inter}}^{\mathit{L}}$ is obtained, which indicates that beams of the LIDAR sensor have detected an obstacle.

$$\begin{array}{cc}\hfill & {\mathit{f}}_{\mathit{Inter}}^{\mathit{L}}={\left(\begin{array}{ccccc}{f}_{Inter}^0& \cdots & f_{Inter}^{i_L}& \cdots & f_{Inter}^{n_L-1}\end{array}\right)}^T,f_{Inter}^{i_L}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{\mathit{F}}_{\mathit{Inter}}^{{\mathit{i}}_{\mathit{L}}}·{\mathit{F}}_{\mathit{Inter}}^{{\mathit{i}}_{\mathit{L}}}T\ne 0\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \\ \hfill & {\mathit{F}}_{\mathit{Inter}}^{{\mathit{i}}_{\mathit{L}}}ϵM_{1\times n_O}\left(\mathbb{R}\right)\hfill \end{array}$$

(13)

Once the beams that have detected obstacles are known, the measured distance vector ${\mathit{D}}_{\mathit{m}}^{\mathit{L}}$ is defined by the Equation (14). This LIDAR modelling considers the perspective effects that affect to this kind of sensors [1,21,26,32], since from each beam only the closest distance ($d_{min}^{i_L}$) is taken, while the later points of intersection with other segments (obstacles) are ignored.

$$\begin{array}{cc}\hfill & {\mathit{D}}_{\mathit{m}}^{\mathit{L}}={\left(\begin{array}{ccccc}{d}_m^0& \cdots & d_m^{i_L}& \cdots & d_m^{n_L-1}\end{array}\right)}^T,d_m^{i_L}=\left\{\begin{array}{cc}{d}_{min}^{i_L}+{\eta }^{i_L}\hfill & \mathrm{if}{f}_{Inter}^{i_L}=1\hfill \\ d_{range}\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \\ \hfill & d_{min}^{i_L}=min\left({\mathit{ND}}_{{\mathit{L}}_{\mathit{E}}}^{{\mathit{i}}_{\mathit{L}}}\right),{\mathit{ND}}_{{\mathit{L}}_{\mathit{E}}}^{{\mathit{i}}_{\mathit{L}}}ϵM_{1\times n_O}\left(\mathbb{R}\right)\hfill \end{array}$$

(14)

where ${\eta }^{i_L}$ is the WGN that affects to each distance measurement, whose standard deviation (${\sigma }^{i_L}$) depends of $d_{min}^{i_L}$ accoding to the Equation (15). Thus, the accuracy of each measure of the sensor is modelled with continuous limits and depends on the distance to the obstacle.

$${\sigma }^{i_L}=\left(tanh\left({\displaystyle \frac{10d_{min}^{i_L}}{d_{range}}}-5\right)+1\right)\left({\displaystyle \frac{{\sigma }_{max}^L-{\sigma }_{min}^L}{2}}\right)+{\sigma }_{min}^L$$

(15)

where the parameter ${\sigma }_{max}^L$ must be set according to the sea states. Thereby, the performance of the sensor depends on the distance of the obstacles as well as of the environmental conditions. As maximum noise level, in this work ${\sigma }_{max}^L=5$ m has been set.

Together with the measured distance matrix (${\mathit{D}}_{\mathit{m}}^{\mathit{L}}$) and intersection flags (${\mathit{f}}_{\mathit{Inter}}^{\mathit{L}}$), the sensor model also provides the matrix of detected points (16). Figure 6 shows the measurements delivered by the modelling proposed in this work for a LIDAR sensor.

$${{\mathit{P}}_{\mathit{obst}}^{\mathit{L}}}_{\mathit{B}}={\left(\begin{array}{ccccc}{x_{obst}^0}_B& \cdots & {x_{obst}^{i_L}}_B& \cdots & {x_{obst}^{n_L-1}}_B\\ {y_{obst}^0}_B& \cdots & {y_{obst}^{i_L}}_B& \cdots & {y_{obst}^{n_L-1}}_B\end{array}\right)}^T,\begin{array}{c}{x_{obst}}_B^{i_L}=cos\left(i_L{h}_{res}\right)d_m^{i_L}\hfill \\ {y_{obst}}_B^{i_L}=sin\left(i_L{h}_{res}\right)d_m^{i_L}\hfill \end{array}$$

(16)

where the reference system of the LIDAR sensor coincides with the Body-axes reference system.

This modelling, with respect to other works [10,12,13,15,21,22,23,24,25,26,27,28,29,30,31] that also evaluate avoidance systems for USVs through numerical simulations, is a novel contribution. Because it exposes the avoidance systems to a higher level of the uncertainty present in the environment surrounding the vessel, without reaching the complexity required when a complete modelling of the environment is realized [32,33,34,35,36]. This facilitates its application in autotuning environments based on evolutionary algorithms, due to the large number of simulations that these algorithms must carry out. As a visual demonstration, the video [63] shows the LROABRA algorithm adapted to the approach for the LIDAR sensor model proposed in this work. As can be seen in this video, this modelling adds measurement error to the detected points in function of their distance to the vehicle. In addition, this model also takes into account the uncertainty present in the environment, since it does not detect hidden obstacles due to the effects of perspective.

2.4. Scenarios for Static Obstacle Avoidance Methods Testing

Before the description of the obstacle scenarios, it is necessary to emphasize that the SOA algorithms studied in this work are reactive methods [11]. These methods are characterized by the use of environmental information processed in real time from the measurements delivered by different sensors. Therefore, these methods do not ensure that the vehicle reaches the target if there are local minima. For this reason, in order to provide a complete solution to the problem of autonomous navigation [1,2,3], it is necessary to combine obstacle avoidance methods or reactive methods [8,10,12,13,14,15,21,22,23,24,25,26,27,28,29,30,31] with global or path planning algorithms [3,7,17,19,20]. In this sense, the obstacle scenarios studied in this work are focused only on the reactive part, unknown scenarios that the autonomous vehicle discovers while navigating. Without considering scenarios with local minimums such as those that can be found in marine environments (in which it would be necessary to have a guidance system that also incorporates a global algorithm). Specifically, the scenarios used to carry out the autotuning and to evaluate the SOA methods studied in this work are shown in Figure 7. A scenario is defined by the obstacles (9), the goal point $P_{goal}$ and the initial position and speed vectors ($\mathit{\eta }$, ${\mathit{\nu }}_{\mathit{r}}$) of the USV. In this work, five scenarios are proposed. With scenario 1 the implementation of the SOA methods is verified. The second scenario, mainly, is used to evaluate the stability of the avoidance algorithm in situations where there are close obstacles without risk of collision. As well as the subsequent response of these methods when a collision risk is presented. Scenario 3 places the USV in a zig-zag scenario where all directions, at greater or lesser distances, lead to collision. Thus, the behaviour of SOA methods in closed environments is evaluated. Finally, scenarios 4 and 5 are densely occupied and allow a wide range of alternative routes. These scenarios are used to study the robustness of SOA methods that, while these have been adjusted for a different scenario, must solve new geometric situations not contemplated in their initial adjustment. In all scenarios the vehicle starts moving at the goal speed $u_{goal}=u=7$ m/s and, looking for the most unfavourable case in which the current leads the USV towards the obstacles, $V_c=0.5$ m/s and ${\beta }_C=0$${}^{\circ }$ has been set.

3. Static Obstacle Avoidance Methods

In order to evaluate the new autotuning environment proposed in this work, three obstacle avoidance methods have been adapted and implemented to the LIDAR sensor Model (10).

3.1. LROABRA Method

The LROABRA method is proposed in [29] as a static obstacle avoidance system that allows USVs to navigate at high speeds. In this approach, the obstacles are represented in a system of polar coordinates referenced to Body-axes, where each obstacle is modelled as a circle of radius $r_i$ and position $O_i({\rho }_i,{\theta }_i)$. From these circles, the non-permitted directions are defined:

$$B_{obs}=\bigcup _{i=1}^{n_O}{R}_{osb}^i,R_{osb}^i=\left\{\begin{array}{cc}\left[\psi \left(k\right)+{\theta }_i-{\varphi }_i-\sigma ,\psi \left(k\right)+{\theta }_i+{\varphi }_i+\sigma \right]\hfill & \mathrm{if}{\rho }_i\le d_{near}\hfill \\ \varnothing \hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(17)

where $\sigma$ is a safety angle used to expand the angular space occupied by obstacles, $d_{near}$ is the threshold distance from which the obstacles are considered and ${\varphi }_i=asin(r_i/{\rho }_i)$.

Based on [64], the guidance directions ($\theta$) and speed commands ($s{p}_u$) that can be reached by the USV are limited to the following dynamic window:

$$\begin{array}{cc}\hfill V_T& =\left\{s{p}_u|s{p}_uϵ\left[u\left(k\right)-{\dot{u}}_{max}{T}_m^a,u\left(k\right)+{\dot{u}}_{max}{T}_m^a\right]\right\}\hfill \\ \hfill V_{Head}& =\left\{\theta |\theta ϵ\left[\psi \left(k\right)+r\left(k\right)T_m^a-{\displaystyle \frac{1}{2}}{\dot{r}}_{max}{T_m^a}^2,\psi \left(k\right)+r\left(k\right)T_m^a+{\displaystyle \frac{1}{2}}{\dot{r}}_{max}{T_m^a}^2\right]\right\}\hfill \end{array}$$

(18)

where ${\dot{u}}_{max}$ and ${\dot{r}}_{max}$ represent the maximum linear and angular acceleration of the vessel, respectively, and $T_m^a$ is the sample time of the algorithm.

Once the sets of forbidden ($B_{obs}$) and reachable ($V_{Head}$) orientations have been defined, the Heuristics (19) is maximized in order to obtain the angle ${\theta }_{optimal}$.

$$\begin{array}{cc}\hfill & F_{ϵ}\left(\theta \right)=ϵ\underset{f_{head}\left(\theta \right)}{\underbrace{\left(1-{\displaystyle \frac{|{\theta }_{goal}-\theta |}{\pi }}\right)}}+(1-ϵ)\underset{f_{safe}\left(\theta \right)}{\underbrace{\left({\displaystyle \frac{1}{n_O}}\sum _{i=1}^{n_O}{\displaystyle \frac{|\theta -({\theta }_i+\psi \left(k\right)\left)\right|}{\pi }}\right)}}\hfill \\ \hfill & \theta \in V_{Head},\theta \notin B_{obs},\theta \in [-\pi ,\pi ]\hfill \end{array}$$

(19)

where $ϵ$ is the heuristic adjustment parameter and ${\theta }_{goal}$ represents the goal direction.

Due to the fact that the angle ${\theta }_{optimal}$ defines a safe and reachable direction for the USV, it is used to insert a temporary guidance waypoint $P_{insert}(x_{insert},y_{insert})$ at a distance $D_{insert}$.

$$\begin{array}{cc}\hfill & x_{insert}=x_E\left(k\right)+D_{insert}cos\left({\theta }_{optimal}\right)\hfill \\ \hfill & y_{insert}=y_E\left(k\right)+D_{insert}sin\left({\theta }_{optimal}\right)\hfill \end{array},s{p}_{\psi }\left(k\right)=atan2\left(y_{insert}-y_E\left(k\right),x_{insert}-x_E\left(k\right)\right)$$

(20)

where $P_{insert}$ is updated when the USV reaches it or if the guidance angle is far from the optimal angle:

$$\begin{array}{cc}\hfill & |x_{insert}-x_E\left(k\right),y_{insert}-y_E\left(k\right)|\le d_{gate}\hfill \\ \hfill & |s{p}_{\psi }\left(k\right)-{\theta }_{optimal}\left(k\right)|\ge {\theta }_{gate}\hfill \end{array}$$

(21)

where $d_{gate}$ and ${\theta }_{gate}$ represent the linear and angular distance thresholds that set the $P_{insert}$ update.

On the other hand, the speed setpoint $s{p}_u$ is calculated according to the Equation (22), which is limited to the dynamic window $V_T$ (18).

$$s{p}_u\left(k\right)=u_{goal}{\eta }_L\left(k\right)\left(1-{\displaystyle \frac{|{\theta }_{\delta }\left(k\right)|}{\pi }}\right),{\theta }_{\delta }\left(k\right)=s{p}_{\psi }\left(k\right)-\psi \left(k\right)$$

(22)

where the variable ${\eta }_L$ depends of the minimum distance to obstacles which lie within the impact zone of the USV, see Equation (23).

$${\eta }_L=\left\{\begin{array}{cc}{\eta }_1={\displaystyle \frac{r\left(k\right)}{max({\theta }_{\delta }\left(k\right),r\left(k\right))}}\hfill & \mathrm{if}{d}_{min}\le d_{near}\hfill \\ {\eta }_2={\displaystyle \frac{d_{min}^2+2r_{max}{d}_{min}}{d_{near}^2+r_{max}{d}_{near}}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.,d_{min}=min\left\{d_{obs}^i\le d_{near}\wedge \left|{\theta }_i\right|<{\theta }_v\right\}$$

(23)

where $r_{max}$ is the maximum angular rate of the vessel, $d_{obs}^i$ is the distance to the center of each obstacle and ${\theta }_v$ defines the angular impact zone of the USV.

In order to expose the LROABRA algorithm to a higher level of uncertainty in the environment, in this work the set of non-permitted directions is modified. Thus, $B_{osb}$ is generated from the measurements delivered by the LIDAR sensor model, considering each beam from the Model (10) as an obstacle, see Equation (24). Furthermore, with the purpose that the algorithm considers the direction of the velocity vector, in the Equations (18), (19) and (22), $\psi \left(k\right)$ is replaced by $\chi \left(k\right)$. One of the executions of the LROABRA algorithm implemented in this work is shown in Figure 8.

$$B_{obs}=\bigcup _{i=1}^{n_L}{R}_{osb}^i,R_{osb}^i=\left\{\begin{array}{cc}\left[\chi \left(k\right)+i_L{h}_{res}-\sigma ,\chi \left(k\right)+i_L{h}_{res}+\sigma \right]\hfill & \mathrm{if}{d}_m^{i_L}\le d_{near}\hfill \\ \varnothing \hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(24)

where the term ${\theta }_i\pm {\varphi }_i$ of the Equation (17) has been replaced by $i_L{h}_{res}$.

3.2. Potential Fields Methods

The potential field (PF) method, due to its simplicity, has been widely used in different problems of mobile robotics [7,11,16,19,23,26,30,45,65]. A disadvantage of the classical potential fields used as a reactive obstacle avoidance method is that it can produce oscillatory trajectories [46]. In order to improve their performance in USVs, the classical potential fields [65] are modified by different authors in [23,30]. In their work, the forces of attraction and repulsion are modelled according to Equation (25).

$$\begin{array}{cc}\hfill & f_{att}=m_{att}{\alpha }_{att}{|P_{goal}-P_{USV}\left(k\right)|}^{m_{att}-1}\hfill \\ \hfill & f_{rep}^i=\underset{f_{rep1}^i}{\underbrace{{\eta }_{rep}\left({\displaystyle \frac{1}{p_s^i}}-{\displaystyle \frac{1}{p_o}}\right){\displaystyle \frac{1}{{p_s^i}^2}}{|P_{USV}\left(k\right)-P_{goal}|}^{n_{rep}}}}+\underset{f_{rep2}^i}{\underbrace{{\displaystyle \frac{n_{rep}}{2}}{\eta }_{rep}{\left({\displaystyle \frac{1}{p_s^i}}-{\displaystyle \frac{1}{p_o}}\right)}^2{|P_{USV}\left(k\right)-P_{goal}|}^{n_{rep}-1}}}\hfill \end{array}$$

(25)

where $P_{USV}\left(k\right)=\left[x_E\left(k\right)y_E\left(k\right)\right]$, $P_{goal}$ is the position of the goal waypoint, ${\alpha }_{att}$ and $m_{att}$ are positive parameters used to tune the magnitude of the attraction potential, $p_o$ is the distance that limits the range of influence of the obstacles, $p_s^i$ is the distance to each obstacle, and, finally, the magnitude of the repulsion potential is tuned with the positive parameters ${\eta }_{rep}$ and $n_{rep}$.

The repulsive forces of Equation (25), as shown in [23], are generated from the vertices of the obstacles which form the coastline. In addition, the distances considered in [23] are in the range of kilometres. In order to expose the algorithm to a higher level of uncertainty in the environment and evaluate it in proximity situations ($d_{range}=200$ m), which increases the difficulty of avoidance scenarios, in this work a repulsive force is generated by each of the beams of the LIDAR sensor (10). In this way, the total repulsive force, the attraction force and the resulting force, all of them referenced to Body-axes, are defined according to Equation (26).

$$\begin{array}{cc}\hfill & F_{rep}\left(k\right)=\sum _{i_L=1}^{n_L}{f}_{rep}^{i_L}\left(\begin{array}{cc}cos\left(i_L{h}_{res}-\zeta \left(k\right)\right)& sin\left(i_L{h}_{res}-\zeta \left(k\right)\right)\end{array}\right)\hfill \\ \hfill & F_{att}\left(k\right)=f_{att}\left(\begin{array}{cc}cos\left({\theta }_{goal}\left(k\right)-\chi \left(k\right)\right)& sin\left({\theta }_{goal}\left(k\right)-\chi \left(k\right)\right)\end{array}\right)\hfill \\ \hfill & F_R\left(k\right)=F_{rep}\left(k\right)+F_{att}\left(k\right)=\left(F_{Rx}{F}_{Ry}\right)\hfill \end{array}$$

(26)

where $\zeta \left(k\right)=\psi \left(k\right)-\chi \left(k\right)$, ${\theta }_{goal}\left(k\right)=atan2\left(y_{goal}-y_E\left(k\right),x_{goal}-x_E\left(k\right)\right)$ and the module of each repulsive force ($f_{rep}^{i_L}$) corresponds to term $f_{rep}^i$ of Equation (25), in which $p_s^i$ has been replaced by the measured distance $d_m^{i_L}$, which is defined in Equation (14).

The course setpoint $s{p}_{\psi }$ is obtained from the resulting force ($F_R$) according to the Equation (27), [46]. As can be seen, a low-pass filter is applied to reduce the oscillations [66].

$$\begin{array}{cc}\hfill {\theta }_{FR}\left(k\right)& =atan2\left(F_{Ry},F_{Rx}\right)+\chi \left(k\right)\hfill \\ \hfill s{p}_{\psi }\left(k\right)& ={\displaystyle \frac{1}{c_F+1}}\left({\theta }_{FR}\left(k\right)+c_{F}s{p}_{\psi }(k-1)\right)\hfill \end{array}$$

(27)

where $c_F={\tau }_F/T_m^a$, ${\tau }_F$ is the effective time constant of the filter applied to the resulting force and $T_m^a$ is the sample time of the obstacle avoidance method.

On the other hand, in [23,30] the authors do not consider changes in the speed of the vessel. In order to ensure that the potential fields have the same manoeuvrability over the vessel as the other SOA methods studied in this paper, the approach proposed in [66] is adopted to calculate $s{p}_u$. Thus, the speed setpoint depends on the resulting force and the speed of the vessel, see Equation (28).

$$\begin{array}{cc}\hfill & s{p}_u\left(k\right)=u_{goal}(1-\left|cos\left(\theta \right)\right|)\hfill \\ \hfill & cos\left(\theta \right)=\left\{\begin{array}{cc}0\hfill & \mathrm{if}|F_R|=0\hfill \\ {\displaystyle \frac{u{F}_{Rx}+v{F}_{Ry}}{|U_B\left|\right|F_R|}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(28)

where $U_B(u,v)$ represents the vessel’s velocity vector referenced to Body-axes.

As will be shown in the results section, despite the favourable results obtained in [23,30] when applying the modified potential fields (PFs) (for distances of several kilometres), if the detection range is limited ($d_{range}=200$ m [59]) the PFs [23,30] do not provide acceptable results. As a secondary contribution of this work, and based on the concept of generalized potential fields (GPFs) proposed in [45] and studied in [46], we propose to use potential fields that consider the relative direction between the obstacle and the vessel’s velocity vector, see Equation (29). For this purpose, the module of repulsive forces $\left(f_{rep}^{i_L}\right)$ calculated in (26) is replaced by:

$$f_{repG}^{i_L}=f_{rep}^{i_L}\left|cos{\left({\varphi }^{i_L}\right)}^{{\beta }_{gpf}}\right|,{\varphi }^{i_L}=\left\{\begin{array}{cc}\pi /2\hfill & \mathrm{if}|i_L{h}_{res}-\zeta \left(k\right)|>\pi /2\hfill \\ i_L{h}_{res}-\zeta \left(k\right)\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(29)

where ${\beta }_{gpf}$ is a positive tuning parameter that decreases the effect of obstacles towards which the USV does not go.

The attraction force in the GPFs method is calculated as in the PFs proposed in [23].Therefore, any change in the guidance direction of one method with respect to another depends only on the repulsive forces. A comparative between the repulsive forces generated when applying the PFs proposed in [23,30] and those obtained when applying the GPFs proposed in this paper is shown in Figure 9. As can be seen in the bottom left graph, in the case of GPFs, only the obstacles towards which the USV is moving produce a repulsive force (blue circles in Figure 9). While, in the case of PFs, all obstacles under the range of influence ($p_o$) produce repulsive forces (brown circles in Figure 9). As a consequence, at the simulation instant shown, the resultant force obtained when applying GPFs ($F{t}_{GPF}$) guides the USV towards the goal, while if the resultant force obtained by the PFs ($F{t}_{PF}$) were applied, this would drive the USV away from the obstacle zone (although these do not represent a real danger) and, consequently, away from the goal (see bottom right graph). In the simulation shown, the GPFs method is the one that guides the USV, which has been manually adjusted in an iterative procedure.

3.3. VFH+ Method

The VFH+ method is proposed in [67] as an improved version of the VFH (Vector Field Histogram) [68]. Later, it has been used as an obstacle avoidance system for an USV in [10]. In this work, the block of this algorithm available in Simulink is used [69]. As model of the environment, the VFH+ method uses The Primary Polar Histogram (${\mathit{H}}_{\mathit{P}}$), [67]. To generate it, the version of the method available in [69] is fully compatible with the model of the LIDAR sensor proposed in this work. As inputs, this block uses the distance measurements (14) and their angular positions (30). Furthermore, as in the other SOA methods, and in order to consider the effect of environmental disturbances [5,52], in this work the angular positions of the measured distances are referenced to the course angle $\chi \left(k\right)$.

$${\mathit{\theta }}_{\mathit{m}}^{\mathit{L}}={\left(\begin{array}{ccccc}{\theta }_m^0& \cdots & {\theta }_m^{i_L}& \cdots & {\theta }_m^{n_L-1}\end{array}\right)}^T,{\theta }_m^{i_L}=i_L{h}_{res}-atan2(v\left(k\right),u\left(k\right))$$

(30)

In the same way that the LROABRA method increases the angular occupation of obstacles with the parameter $\sigma$, and in contrast to the potential fields which do not consider the size of the USV, the VFH+ method approximates the shape of the vehicle to a circle of radius $r_{width}$. Thus, the angular occupation of each detected point is extended according to ${\gamma }_r^{i_L}$, see Equation (31). In this way, each bar of ${\mathit{H}}^{\mathit{P}}$ is defined according to the Equation (32).

$${\mathit{\gamma }}_{\mathit{r}}^{\mathit{L}}={\left(\begin{array}{ccccc}{\gamma }_r^0& \cdots & {\gamma }_r^{i_L}& \cdots & {\gamma }_r^{n_L-1}\end{array}\right)}^T,{\gamma }_r^{i_L}=\left\{\begin{array}{cc}asin\left(r_{width}/d_m^{i_L}\right)\hfill & \mathrm{if}{f}_{Inter}^{i_L}=1\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(31)

$$\begin{array}{cc}\hfill & {\mathit{H}}_{\mathit{P}}={\left(\begin{array}{ccccc}{h}_P^0& \cdots & h_P^{i_L}& \cdots & h_P^{n_L-1}\end{array}\right)}^T,h_P^{i_L}=\sum _{i_L=1}^{n_O}f\left(d_m^{i_L}\right)b^{i_L}\hfill \\ \hfill & b^{i_L}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{i}_L{h}_{res}\in [{\theta }_m^{i_L}-{\gamma }_r^{i_L},{\theta }_m^{i_L}+{\gamma }_r^{i_L}]\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(32)

where $f\left(d_m^{i_L}\right)$ defines the occupation density as a function of distance, see [67,69].

Obtained ${\mathit{H}}_{\mathit{P}}$, this is converted into a Binary Polar Histogram (${\mathit{H}}_{\mathit{B}}$) that classifies the directions between occupied or safe. To ensure that the state of the directions does not change several times, which generates an oscillatory behavior that can drive the system to dangerous situations [67], the hysteresis thresholds ${\tau }_{high}$ and ${\tau }_{low}$ are applied:

$$\begin{array}{cc}\hfill & {\mathit{H}}_{\mathit{B}}={\left(\begin{array}{ccccc}{h}_B^0& \cdots & h_B^{i_L}& \cdots & h_B^{n_L-1}\end{array}\right)}^T,h_B^{i_L}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}{h}_P^{i_L}\left(k\right)>{\tau }_{high}\hfill \\ 0\hfill & \mathrm{if}{h}_P^{i_L}\left(k\right)<{\tau }_{low}\hfill \\ h_B^{i_L}(k-1)\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(33)

In addition, in order to take into account the dynamic and kinematic of the vehicle, the VFH+ method simplifies its trajectory to circular arcs of radius $r_{ship}$,67]. These circular trajectories limit the safe directions, applying a mask to the ${\mathit{H}}_{\mathit{B}}$ and generating The Masked Polar Histogram:

$$\begin{array}{c}{\mathit{H}}_{\mathit{M}}={\left(\begin{array}{ccccc}{h}_M^0& \cdots & h_M^{i_L}& \cdots & h_M^{n_L-1}\end{array}\right)}^T,h_M^{i_L}=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{h}_B^{i_L}\left(k\right)=0and{\theta }_m^{i_L}\in \left\{[{\varphi }_r,\chi ],[\chi ,{\varphi }_l]\right\}\hfill \\ 1\hfill & \mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(34)

where ${\varphi }_l$ and ${\varphi }_r$ are the safe angles to port and starboard, respectively. The definition of these angles is recorded in [67].

Defined ${\mathit{H}}_{\mathit{M}}$, the course setpoint ($s{p}_{\psi }$) is chosen between the safe or collision-free directions obtained ($h_M^{i_L}=0$). For this purpose, the following cost function is minimized:

$$g\left({\theta }_m^{i_L}\right)={\mu }_1\mathsf{\Delta }\left({\theta }_m^{i_L},{\theta }_{goal}\right)+{\mu }_2\mathsf{\Delta }\left({\theta }_m^{i_L},\chi \left(k\right)\right)+{\mu }_3\mathsf{\Delta }\left({\theta }_m^{i_L},s{p}_{\psi }(k-1)\right)$$

(35)

where the parameters ${\mu }_1$, ${\mu }_2$ and ${\mu }_3$ are weighing factors that take into account: alignment with the goal direction, cost associated with the change of direction and alignment with the previous course setpoint (memory effect), respectively. The operator $\delta (c1,c2)$ calculates the angular distance between the directions $c1$ and $c2$.

Finally, the speed setpoint $s{p}_u$ is calculated according to the Equation (36), [68]. Figure 10 shows one execution of the VFH+ method implemented in this work.

$$s{p}_u\left(k\right)=u_{goal}\left(1-{\displaystyle \frac{\left|r\right(k\left)\right|}{r_{max}}}\right)$$

(36)

where $r_{max}$ is a tuning parameter which represents the maximum angular rate of the USV.

4. Simple Autotuning Approach for Obstacle Avoidance Methods

Currently there are a great variety of obstacle avoidance methods for USVs [8,10,12,13,14,15,21,22,23,24,25,26,27,28,29,30,31], as well as a wide range of surface marine vehicles [4,5,8,10,12,21,22,26,28,29,47,48,49,52] in which these can be applied. For this reason, and as main novelty contribution of this work, an AutoTuning Environment for Static Obstacle Avoidance (ATESOA) methods applied to USVs is proposed. This environment is flexible to any type of marine surface vehicle (that is equipped with course and speed controllers) and can be applied to different SOA methods. It should be noted that, starting from the limitations of each avoidance method, its performance to avoid obstacles and to reach the goal are determined by the tuning of its parameters (${\mathsf{\Theta }}^{tuning}$). In addition, this tuning must be carried out on a specific obstacle scenario. However, in the problems of obstacle avoidance there can be an infinite number of geometric combinations that define the scenario. Therefore, even if the optimum value of ${\mathsf{\Theta }}^{tuning}$ is obtained for a defined function, this parameter vector would only be optimal for that particular scenario. For this reason, the study carried out in this work does not pretend to find the optimal parameter vector for each avoidance method applied to a specific vessel. Instead, its purpose is to facilitate the pretuning process necessary to implement SOA methods in USVs, as well as to assist in their choice for a particular vessel, through the use of the new ATESOA proposed in this paper.

In this work, a Genetic Algorithm is used as optimization technique to perform autotuning. This election is due to its established use in vehicle autotuning [37,39,42,44] and to the fact that the methods of obstacle avoidance [8,10,12,13,14,15,21,22,23,24,25,26,27,28,29,30,31] usually include decision factors based on thresholds (that convert them into non-smooth functions, complicating the application of conventional least-squares optimization methods [43,53]). Thus, each tuning parameter of an avoidance method is considered as a gene and, consequently, its tuning parameter vector ${\mathsf{\Theta }}^{tuning}$ represents an individual of the population. In addition, to quantify the performance achieved by an avoidance method, with a particular tuning vector ${\mathsf{\Theta }}^{tuning}$, the following indicators are proposed:

$$\begin{array}{cc}\hfill \mathsf{\Delta }s{p}_{\psi }\left(k\right)& ={\displaystyle \frac{|s{p}_{\psi }\left(k\right)-s{p}_{\psi }(k-1)|}{\pi }}+\mathsf{\Delta }s{p}_{\psi }(k-1)\hfill \\ \hfill \mathsf{\Delta }s{p}_u\left(k\right)& ={\displaystyle \frac{|s{p}_u\left(k\right)-s{p}_u(k-1)|}{u_{max}}}+\mathsf{\Delta }s{p}_u(k-1)\hfill \\ \hfill dis{t}_D\left(k\right)& =|(x_E\left(k\right)-x_E(k-1)),(y_E\left(k\right)-y_E(k-1))|+dis{t}_D(k-1)\hfill \\ \hfill time\left(k\right)& =k{T}_m^A\hfill \\ \hfill dis{t}_G\left(k\right)& =|(x_E\left(k\right)-x_{goal}),(y_E\left(k\right)-y_{goal})|\hfill \end{array}$$

(37)

where $\mathsf{\Delta }s{p}_{\psi }$ and $\mathsf{\Delta }s{p}_u$ are used to quantify the control effort solicited by the avoidance system, $dis{t}_D$ is the distance traveled by the USV, $time$ is the travel time and $dis{t}_G$ is the distance up to $P_{goal}$.

Once the indicators (37) have been defined, the fitness function that quantifies the performance of an avoidance method for each individual ${\mathsf{\Theta }}^{tuning}$ is defined according to:

$$J=\left\{\begin{array}{cc}{J}_s={\displaystyle \frac{time\left(k\right)}{tim{e}^{lim}}}+{\displaystyle \frac{dis{t}_D\left(k\right)}{dis{t}_D^{lim}}}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathsf{\Delta }s{p}_{\psi }\left(k\right)}{\mathsf{\Delta }s{p}_{\psi }^{lim}}}+{\displaystyle \frac{\mathsf{\Delta }s{p}_u\left(k\right)}{\mathsf{\Delta }s{p}_u^{lim}}}\right)\hfill & \mathrm{if}{f}_{stop}=1\hfill \\ J_f=3+{\displaystyle \frac{dis{t}_G\left(k\right)}{dis{t}_G\left(0\right)}}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(38)

where $tim{e}^{lim}$, $dis{t}_D^{lim}$, $\mathsf{\Delta }s{p}_{\psi }^{lim}$ and $\mathsf{\Delta }s{p}_u^{lim}$ set the maximum allowed values for the indicators (37) in a particular scenario and $dis{t}_G\left(0\right)$ represents the initial distance between the USV and $P_{goal}$. On the other hand, $f_{stop}$ defines two situations: the USV has reached its target or not,

$$f_{stop}=\left\{\begin{array}{cc}1\hfill & \mathrm{if}dis{t}_G\left(k\right)\le d_{wp}\hfill \\ 0\hfill & \mathrm{if}\left(n\left(k\right)=n^{lim}\vee min\left({\mathit{D}}_{\mathit{m}}^{\mathit{L}}\right)\le d_{collision}\right)\hfill \end{array}\right.,n\in \left\{time,dis{t}_D,\mathsf{\Delta }s{p}_{\psi },\mathsf{\Delta }s{p}_u\right\}$$

(39)

where $d_{wp}$ and $d_{collision}$ represent the distance thresholds, which define whether the USV has reached $P_{goal}$ or has collided, respectively.

In this way, if the parameter vector ${\mathsf{\Theta }}^{tuning}$ drives the USV to the goal, within the established limits (39), the index $J_s$ is used to minimize: time, distance and control effort in the travel, see Equation (38). On the other hand, if the safety distance ($d_{collision}$) or any of the established limits is exceeded, the numerical simulation is stopped and the index used is $J_f$. The value of $J_f$ will always be higher than that of $J_s$ and, furthermore, it depends on the final distance between the USV and $P_{goal}$. This last facilitates that, if in the initial populations a parameter vector that successfully guides the USV to $P_{goal}$ is not obtained, the GA prioritizes those vectors ${\mathsf{\Theta }}^{tuning}$ that drive the USV closest to the goal. Once the fitness function (38) has been described, the tuning parameter vectors of the four SOA methods studied in this work are defined in the Equation (40).

$$\begin{array}{cc}\hfill & {\mathsf{\Theta }}_{LROABRA}^{tuning}={\left(\begin{array}{ccccccccc}{\dot{u}}_{max}& {\dot{r}}_{max}& ϵ& D_{insert}& (d_{near}-D_{insert})& \sigma & {\theta }_{gate}& r_{max}& {\theta }_v\end{array}\right)}^T\hfill \\ \hfill & {\mathsf{\Theta }}_{PF}^{tuning}={\left(\begin{array}{cccccc}{m}_{att}& {\alpha }_{att}& {\eta }_{rep}& n_{rep}& p_o& {\tau }_{F_R}\end{array}\right)}^T,{\mathsf{\Theta }}_{GPF}^{tuning}={\left(\begin{array}{cc}{{\mathsf{\Theta }}_{PF}^{tuning}}^T& {\beta }_{GPF}\end{array}\right)}^T\hfill \\ \hfill & {\mathsf{\Theta }}_{VFH+}^{tuning}={\left(\begin{array}{ccccccc}{\tau }_{low}& ({\tau }_{high}-{\tau }_{low})& r_{ship}& {\mu }_1& {\mu }_2& {\mu }_3& r_{max}\end{array}\right)}^T\hfill \end{array}$$

(40)

For the implementation of the AG, the Toolbox available in Matlab [70] has been used. In order to take into account the physical or adjustment limits established for the vectors ${\mathsf{\Theta }}^{tuning}$, these have been limited considering the guidelines exposed in [23,29,30,66,67,68,69], see Equation (41). Consequently, the option mutationadaptfeasible has been established as a mutation function [70] and the initial population has been limited to $x^{low}$ and $x^{hihg}$, see Equation (41). On the other hand, in order to limit the simulation times, we have established: a population of 200 individuals, a maximum of 20 generations, and the parameter $MaxStallGenerations$ as $25\%$ of the maximum generations. Finally, with the purpose of avoiding early convergence [37,70], elite individuals have been limited to $1.5\%$ of the population and the crossover factor has been decreased (thus increasing the mutation factor) to 0.7. The other parameters are the ones established by default in [70].

$$\begin{array}{c}{\mathsf{\Theta }}_{PF}^{low}={\left(\begin{array}{cccccc}0.1& 50& 8000& 0.1& 50& 1\end{array}\right)}^T,{\mathsf{\Theta }}_{PF}^{high}={\left(\begin{array}{cccccc}5& 1000& 2.0\times {10}^5& 5& 200& 6\end{array}\right)}^T\hfill \\ {\mathsf{\Theta }}_{LROABRA}^{low}={\left(\begin{array}{ccccccccc}0.2& 0.3& 0.5& 60& 0& \pi /18& \pi /4& 0.1& \pi /36\end{array}\right)}^T,{\mathsf{\Theta }}_{GPF}^{low}={\left(\begin{array}{cc}{{\mathsf{\Theta }}_{PF}^{low}}^T& 2\end{array}\right)}^T\hfill \\ {\mathsf{\Theta }}_{LROABRA}^{up}={\left(\begin{array}{ccccccccc}2& 2.5& 0.9& 100& 100& \pi /6& 3\pi /4& 1& \pi /6\end{array}\right)}^T,{\mathsf{\Theta }}_{GPF}^{high}={\left(\begin{array}{cc}{{\mathsf{\Theta }}_{PF}^{high}}^T& 15\end{array}\right)}^T\hfill \\ {\mathsf{\Theta }}_{VFH+}^{low}={\left(\begin{array}{ccccccc}50& 10& 20& 10& 0& 0& 0.1\end{array}\right)}^T,{\mathsf{\Theta }}_{VFH+}^{up}={\left(\begin{array}{ccccccc}500& 1000& 100& 20& 10& 10& 0.5\end{array}\right)}^T\hfill \end{array}$$

(41)

where the respective magnitudes are expressed in the international system of units.

5. Results

This section shows the results obtained by the new ATESOA applied to the USV model proposed in this work. For numerical simulations, Runge-Kutta numerical integration method of fourth order is used with an integration step of 0.01 s and the SOA methods have a sampling period of $T_m^a=1$ s. Specifically, four SOA methods (LROABRA, PF, GPF and VFH+) have been automatically adjusted for the USV formed by the vessel Model (1), which is governed by the Controllers (7). In addition, in order to expose the algorithms to a higher level of uncertainty in the environment, these have been adapted to the measurements provided by the LIDAR sensor Model (10) proposed in this work, see Section 3. On the other hand, the scenarios 1, 2 and 3, which are defined in Section 2.4, have been used to carry out the autotuning of the SOA methods. As a result, the parameter vectors of each obstacle avoidance method (automatically adjusted for a scenario following the procedure given in Section 4) are shown in the

Table 4. Parameter vectors obtained by ATESOA in the scenarios 1, 2 and 3.

MethodScenary	Parameter Vectors
LROABRAS1	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccccc}1.83& 1.87& 0.68& 62.28& 79.44& 0.19& 0.86& 0.77& 0.51\end{array}\right)}^T$
LROABRAS2	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccccc}1.43& 1.84& 0.63& 67.13& 41.91& 0.19& 1.17& 0.66& 0.11\end{array}\right)}^T$
LROABRAS3	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccccc}1.07& 1.57& 0.56& 60.98& 60.35& 0.18& 2.18& 0.68& 0.09\end{array}\right)}^T$
PFS1	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{cccccc}1.36& 833.16& 1.32\times {10}^5& 1.68& 107.95& 1.54\end{array}\right)}^T$
PFS2	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{cccccc}2.06& 550.53& 2.02\times {10}^4& 1.75& 100.59& 1.50\end{array}\right)}^T$
PFS3	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{cccccc}2.85& 454.25& 1.03\times {10}^5& 4.20& 172.67& 1.91\end{array}\right)}^T$
GPFS1	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}2.62& 909.49& 1.89\times {10}^5& 4.48& 161.7& 2.27& 14.49\end{array}\right)}^T$
GPFS2	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}1.07& 245.22& 1.85\times {10}^5& 0.34& 141.87& 3.28& 10.64\end{array}\right)}^T$
GPFS3	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}0.75& 447.15& 1.55\times {10}^5& 0.50& 120.94& 1.06& 13.75\end{array}\right)}^T$
VFH+S1	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}93.71& 714.1& 50.94& 15.81& 7.58& 0.30& 0.45\end{array}\right)}^T$
VFH+S2	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}305.97& 176.25& 72.96& 13.01& 2.12& 9.71& 0.43\end{array}\right)}^T$
VFH+S3	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}407.43& 40.24& 70.08& 19.84& 3.13& 3.18& 0.47\end{array}\right)}^T$

. Where, for example, $P{F}^{S1}$ represents the PF method autotuning for the first scenario.

In order to validate the proposed autotuning environment and to compare the performance of each obstacle avoidance method studied in this work, Figure 11 and

Table 5. Indicators (37) achieved by each obstacle avoidance method in the scenarios 1, 2 and 3.

Method-Scenary/Indicators	${\mathit{f}}_{\mathit{stop}}$	time [s]	${\mathit{dist}}_{\mathit{D}}$ [m]	$\mathsf{\Delta }{\mathit{sp}}_{\mathit{\psi }}$ [rad]	$\mathsf{\Delta }{\mathit{sp}}_{\mathit{u}}$ [m/s]	${\mathit{dist}}_{\mathit{G}}$ [m]
LROABRAS1 - S1	✓	50	263.64	248.02	1.08	—
PFS1 - S1	✓	109	573.20	446.68	2.33	—
GPFS1 - S1	✓	59	294.23	427.20	1.11	—
VFH+S1 - S1	✓	47	265.50	162.14	1.29	—
LROABRAS2 - S2	✓	82	751.45	154.19	1.74	—
PFS2 - S2	✗	—	—	—	—	93.19
GPFS2 - S2	✗	—	—	—	—	91.62
VFH+S2 - S2	✓	78	739.36	131.86	1.23	—
LROABRAS3 - S3	✓	145	929.25	340.32	2.85	—
PFS3 - S3	✓	384	1993.10	1785.40	6.92	—
GPFS3 - S3	✓	226	1182.10	1837.20	1.05	—
VFH+S3 - S3	✓	148	949.69	428.77	2.28	—

are presented. On the one hand, Figure 11 shows the different trajectories described by the USV, on the three scenarios of autotuning, when each avoidance method is applied. In addition, Table 5 presents the Indicators (37) that define the performance achieved by each avoidance method when guiding the USV. As can be seen, with the exception of the potential fields applied to scenario 2, the new ATESOA achieves parameter vectors that guide the USV safely to $P_{goal}$. A situation like scenario 2 is not evaluated by the authors in [23,30]. Since, in the most similar scenario, the authors in [23] guide the model of the Mariner class vessel (171.8-meter vessel) proposed in [52] by a scenario where the space between obstacles is about 10 km, compared with the Model (1) of 9.2 m length and the 50 m of distance between obstacles defined in this paper. Furthermore, in [23], the goal point is contained between the obstacles, which does not cause the $F_{att}$ force to drive the USV into these obstacles (situation posed in Scenario 2). Additionally, with manual tuning, the PFs and GPFs methods have also failed in this scenario. This situation corresponds to the problem described in [19,71], in which the restriction of potential fields when passing through small passages is exposed. To solve it, it would be necessary to apply more advanced potential fields as the one proposed in [71]. On the other hand, the contribution made in this work through the application of the concept of generalized potential fields [45,46] to the improved PFs for USV proposed in [23,30], notably improves the guidance of the USV in the studied scenarios. This can be seen in scenarios 1 and 3, where travel times have been reduced to $54.13\%$ and $58.85\%$ and distances driven to $51.33\%$ and $59.31\%$, respectively. Nevertheless, the performance obtained, for the measurement range $d_{range}=200$ m and the USV Model (1, 7), by the methods based on potential fields is significantly lower than that obtained by the LROABRA and VFH+ methods, see Table 5. With respect to these, both present very similar indicators in the three scenarios. The performance obtained by the VFH+ is better in scenarios 1 and 2, while in scenario 3 the LROABRA method is superior.

On the other hand, in the problems of obstacle avoidance there are infinite geometric combinations that define the possible scenarios that the autonomous vehicle may encounter. Therefore, once the tuning of an avoidance method for a specific scenario has been made, its robustness for other scenarios with different obstacle distributions must be evaluated. In this aspect, a robust behaviour in the guidance of the USV, which is realized by an avoidance method, that keeps its performance (based on the Indicators (37) or similar) is not the objective. Instead, the capability of the avoidance method to safely guide the USV to $P_{goal}$ in significantly different scenarios than those used for its tuning is analyzed. In this context, each avoidance method has been evaluated, with the parameter vectors collected in Table 4, in scenarios with different obstacle distributions than those used for their autotuning. As a result, Figure 12 shows the paths described by the USV when applying each avoidance method over the scenarios defined in this work for which it has not been adjusted. In addition, as support to the previous figure,

Table 6. Indicators (37) achieved in the guidance of the USV by the SOA methods (LROABRA, PF, GPF and VFH+) in scenarios different to those used for their autotuning.

Method-Scenary/Indicators	${\mathit{f}}_{\mathit{stop}}$	time [s]	${\mathit{dist}}_{\mathit{D}}$ [m]	$\mathsf{\Delta }{\mathit{sp}}_{\mathit{\psi }}$ [rad]	$\mathsf{\Delta }{\mathit{sp}}_{\mathit{u}}$ [m/s]	${\mathit{dist}}_{\mathit{G}}$ [m]
LROABRAS2 - S1	✓	71	350.42	573.68	1.67	—
LROABRAS3 - S1	✓	69	346.06	401.13	1.54	—
LROABRAS1 - S2	✗	—	—	—	—	408.43
LROABRAS3 - S2	✓	82	752.03	133.03	1.65	—
LROABRAS1 - S3	✓	167	951.04	509.81	3.83	—
LROABRAS2 - S3	✓	165	1002.80	638.77	3.83	—
LROABRAS1 - S4	✗	—	—	—	—	362.89
LROABRAS2 - S4	✗	—	—	—	—	357.68
LROABRAS3 - S4	✗	—	—	—	—	160.49
LROABRAS1 - S5	✓	109	780.98	344.93	1.60	—
LROABRAS2 - S5	✓	89	767.74	236.10	2.05	—
LROABRAS3 - S5	✓	88	758.31	232.56	2.01	—
PFS3 - S1	✗	—	—	—	—	103.32
PFS1 - S2	✗	—	—	—	—	419.99
PFS3 - S2	✗	—	—	—	—	406.31
PFS1 - S3	✗	—	—	—	—	646.32
PFS1 - S4	✗	—	—	—	—	422.80
PFS3 - S4	✗	—	—	—	—	418.08
PFS1 - S5	✗	—	—	—	—	303.24
PFS3 - S5	✗	—	—	—	—	139.99
GPFS3 - S1	✓	64	324.85	356.46	1.14	—
GPFS1 - S2	✗	—	—	—	—	323.60
GPFS3 - S2	✗	—	—	—	—	360.60
GPFS1 - S3	✗	—	—	—	—	222.95
GPFS1 - S4	✗	—	—	—	—	226.71
GPFS3 - S4	✗	—	—	—	—	228.07
GPFS1 - S5	✗	—	—	—	—	181.70
GPFS3 - S5	✗	—	—	—	—	177.98
VFH+S2 - S1	✓	48	270.62	164.54	1.33	—
VFH+S3 - S1	✓	47	267.96	220.53	1.27	—
VFH+S1 - S2	✓	80	748.90	130.65	1.39	—
VFH+S3 - S2	✓	78	743.54	129.58	1.21	—
VFH+S1 - S3	✓	177	1081.70	644.07	2.96	—
VFH+S2 - S3	✓	166	1018.60	600.59	3.39	—
VFH+S1 - S4	✓	78	758.63	194.63	1.44	—
VFH+S2 - S4	✓	76	754.82	157.85	1.30	—
VFH+S3 - S4	✓	75	750.52	151.35	1.23	—
VFH+S1 - S5	✓	93	791.22	400.88	1.96	—
VFH+S2 - S5	✓	92	790.23	275.57	1.83	—
VFH+S3 - S5	✓	92	793.59	287.66	1.84	—

shows the Indicators (37) achieved by each avoidance method in the guidance of the USV. In this robustness study, the indicator that presents the greatest interest is $f_{stop}$, which defines success or failure in the vehicle’s guidance. As can be seen, in the scenarios defined in this work, the PFs and GPFs methods do not present good results in the guidance of the USV over scenarios for which these have not been adjusted. Since only a favourable result is obtained by the GPFs method adjusted in Scenario 3 when it guides the USV over Scenario 1 (GPF^S3-S1). In the other scenarios, the methods based on potential fields lead the USV to a collision or generate trajectories that are proper of a limit cycle (PF^S1 in Scenarios 3 and 5). On the other hand, the LROABRA method successfully solves the majority of scenarios for which it has not been adjusted. With two exceptions: Scenario 2 with the autotuning performed for scenario 1 and scenario 4 with any parameter vector. The first case is due to the value $d_{near}=141.72$ m obtained in the autotuning, which constrains the passage between the two parallel obstacles. As a consequence, the avoidance method, without taking into account the USV’s dynamics, requires that the vessel turn completely to starboard, which causes a collision. In the second case, the numerical simulations carried out show that, in scenario 4, the collisions occur because the LROABRA method makes multiple switches in the guidance waypoint $P_{insert}$. These commutations, as established in [67], cause undecided behaviour that leads the USV to a collision. To correct this, in a similar way to the VFH+ method, a possible solution would consist in adding to the heuristics (19) a weighting factor that takes into account the alignment of the candidate directions $V_{Head}$ with the previous course setpoint $s{p}_{\psi }(k-1)$.

With respect to the VFH+ method, it successfully resolves all of the scenarios defined in this paper for which it has not been adjusted. Moreover, it presents a robust behaviour to the variations in its parameter vector, since the Indicators (37) obtained in the same scenario for different values of ${\mathsf{\Theta }}^{tuning}$ are very close, see Table 6. Therefore, for the defined scenarios, of the four methods of static obstacle avoidance applied to USVs that have been adapted and implemented in this work, the VFH+ method achieves the best results in terms of performance (see Table 5) and robustness (see Figure 12).

As shown in the results obtained in Figure 12 and in Table 6, the new ATESOA proposed in this work is flexible, such that it can be applied to different SOA methods used in USVs. Where, the robustness of the autotuning obtained depends on the scenario used, but mainly of the features specific of each obstacle avoidance method. Since, as shown below, the new ATESOA achieves a successful adjustment of the SOA methods that did not overcome the obstacle scenarios 4 and 5. In particular, the autotuning of the LROABRA method for scenario 4 and of the methods based on potentials fields for scenarios 4 and 5 is carried out. The parameter vectors obtained are collected in

Table 7. Parameter vectors obtained for the LROABRA, PF and GPF methods, scenarios 4 and 5.

MethodScenary	Parameter Vectors
LROABRAS4	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccccc}1.92& 0.83& 0.56& 74.94& 74.91& 0.21& 1.41& 0.85& 0.11\end{array}\right)}^T$
PFS4	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{cccccc}1.97& 209.81& 1.10\times {10}^5& 3.24& 34.02& 5.62\end{array}\right)}^T$
PFS5	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{cccccc}0.99& 525.02& 5.46\times {10}^4& 0.73& 42.79& 4.09\end{array}\right)}^T$
GPFS4	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}1.96& 867.47& 5.34\times {10}^4& 3.31& 76.05& 3.92& 11.45\end{array}\right)}^T$
GPFS5	${\mathsf{\Theta }}^{tuning}={\left(\begin{array}{ccccccc}0.70& 726.56& 1.70\times {10}^5& 1.73& 29.30& 5.99& 10.66\end{array}\right)}^T$

. As a result, Figure 13 shows the trajectories described by the USV when applying the LROABRA, PFs and GPFs methods autotuning in the new ATESOA. In addition,

Table 8. Indicators (37) achieved by the LROABRA, PF and GPF methods in the scenarios 4 and 5.

Method-Scenary/Indicators	${\mathit{f}}_{\mathit{stop}}$	time [s]	${\mathit{dist}}_{\mathit{D}}$ [m]	$\mathsf{\Delta }{\mathit{sp}}_{\mathit{\psi }}$ [rad]	$\mathsf{\Delta }{\mathit{sp}}_{\mathit{u}}$ [m/s]	${\mathit{dist}}_{\mathit{G}}$ [m]
LROABRAS4 - S4	✓	97	867.66	290.04	1.12	—
PFS4 - S4	✓	144	1143.20	511.33	4.36	—
GPFS4 - S4	✓	153	1063.40	915.83	2.04	—
PFS5 - S5	✓	137	1047.60	642.90	3.70	—
GPFS5 - S5	✓	95	815.77	443.22	2.67	—

presents the Indicators (37) obtained by each method in each scenario. As can be seen, all methods successfully solve the two scenarios. In scenario 4, the LROABRA method again improves performance with respect to potential fields. Regarding the potential fields [23,30], of the methods studied in this work, the results obtained indicate that these would be the worst option as an avoidance system for USVs in situations where the measurement range of the sensors is limited ($d_{range}\le 200$ m). Although in both scenarios, the GPFs method proposed in this work again produces a less oscillatory trajectory in the guidance than the PFs proposed in [23,30], which means shorter travel distances.

6. Conclusions

Due to the scientific and industrial interest, in the last years many works have been developed in order to provide a higher level of autonomy to the vessels. This has led to a diversity of methods for static obstacle avoidance applied to USVs. In order to design, development and validation of these methods, it is necesary an environment for computer simulation that allows the evaluation of the proposed algorithms. With this aim, the main contribution of this work is the new ATESOA, an autotuning environment for the design parameters associated with each static obstacle avoidance algorithm. This environment is adaptable to diverse SOA methods, which can be applied in different types of USVs. Due to the simplified model of the LIDAR sensor proposed in this paper, the simulation environment is more realistic than those used by most authors for the design of obstacle avoidance methods. With the advantage, over works focused on more sophisticated and complex simulation environments to validate USVs, that this model of the LIDAR sensor facilitates the application of evolutionary algorithms for autotuning (due to its lower computation load). To evaluate the new ATESOA, four SOA methods (LROABRA, PFs, GPFs and VFH+) have been adapted and implemented for using this sensor model. The results show how this new autotuning environment achieves a successful and, depending on the obstacle avoidance method, robust tuning of the parameters of the SOA methods for a range of measurement limited by the LIDAR [59] to 200 m. In addition, the parameters for the mathematical model of a USV of 9.2 m of length, whose modelling is based on previous work in the field of surface marine vehicle control, are proposed and evaluated. This model includes: non-linear modelling of the vessel and actuators dynamics, modelling of the effect of ocean currents and course/velocity controllers. In addition, as method of avoiding reactive obstacles, the performance of potential fields [23,30], applied to this USV model (length of 9.2 m) has been improved, through the adaptation of the generalized potential fields [45,46]. As future work, the new autotuning environment will be adapted for dynamic obstacle avoidance methods.