Tracking Control of Autonomous Underwater Vehicles with Acoustic Localization and Extended Kalman Filter

: The absence of global positioning system (GPS) signals and the inﬂuence of ocean currents are two of the main challenges facing the autonomy of autonomous underwater vehicles (AUVs). This paper proposes an acoustic localization-based tracking control method for AUVs. Particularly, three buoys that emit acoustic signals periodically are deployed over the surface. Times of arrivals of these acoustic signals at the AUV are then obtained and used to calculate an estimated position of the AUV. Moreover, the uncertainties involved in the localization and ocean currents are handled together in the framework of the extended Kalman ﬁlter. To deal with system physical constraints, model predictive control relying on online repetitive optimizations is applied in the tracking controller design. Furthermore, due to the different sampling times between localization and control, the dead-reckoning technique is utilized considering detailed AUV dynamics. To avoid using the highly nonlinear and complicated AUV dynamics in the online optimizations, successive linearizations are employed to achieve a trade-off between computational complexity and control performance. Simulation results show that the proposed algorithms are effective and can achieve the AUV tracking control goals.


Introduction
Autonomous underwater vehicles (AUVs) have been widely used for various purposes, e.g., military, commercial, and marine scientific survey applications [1]. For different purposes and scenarios, the requirement for the localization and control accuracy of AUVs ranges from meters to hundreds of meters. Compared to manned underwater vehicles, AUVs have no human pilots on board and thus avoid being exposed to dull and dangerous underwater environments. Compared to the tethered unmanned underwater vehicles, i.e., remotely operated vehicles (ROVs) [2], AUVs are tether free and thus can reach wider ocean spaces without the technical issues brought by the tether cable. The perception-planaction cycle [3] is achieved by AUVs without persistent human instructions in order to accomplish the programmed missions. However, in the underwater environment, the attenuation rate of the global positioning system (GPS) signals is high, which makes GPS not applicable. Moreover, the ubiquitous ocean currents pose great influences on the relatively slow-moving AUVs. These factors impede the realization of the autonomy of AUVs.
In general, there are two types of localization methods underwater, i.e., distance independent methods and distance dependent methods [4]. Distance independent methods estimate the positions of underwater targets based on the topology of several beacon nodes. The most commonly used distance dependent methods include received signal strength indicator, time difference of arrival (TDOA), and time of arrival (TOA). They differ in the distance calculation processes. Multipath and shadow fading effects are usually concerned with the received signal strength indicator method. TOA methods use fewer states, and the uncertainties due to the ocean current, an EKF approach is proposed to provide an estimate of the system's full state. With the estimated system states, MPC is applied to track a reference trajectory. System physical constraints and control performance are explicitly considered. Since the localization and control have different sampling times, the dead-reckoning technique is utilized with detailed AUV dynamics. In the MPC controller design, to relieve the possible computational burden, the AUV dynamics are simplified through successive linearizations. Briefly, this paper contributes to the state of the art in the following aspects: The more precise definition of the body-fixed AUV velocity is the velocity of the point O b with respect to {n} expressed in {b}. The ocean current velocity vector ν c is also defined in the body-fixed coordinate system {b}, as to be specified later. The controlled forces and moments are τ = ξ Y uuδ u 2 δ N uuδ u 2 δ T where ξ and δ, i.e., the system inputs, are the propeller thrust and rudder angle, respectively, and Y uuδ and N uuδ are the corresponding hydrodynamic parameters. The rotation matrix R(ψ) relates the motions in {n} and {b}: where ψ is the heading angle, and R satisfies d dt {R(ψ)} =ψR(ψ)S with In Equation (2), ν r = ν − ν c = u r v r r r T is the relative velocity in {b} between the AUV hull and the ocean current. For simplification, we consider constant irrotational ocean current which is denoted as b = [V c , β c ] T where V c and β c are the current speed and angle in the inertial coordinate system {n}, respectively. Therefore, β c = 0 means that the direction of the ocean current is aligned with the x n axis. Then Therefore,ν r =ν − rS T ν c .
Rigid-body and added mass matrices are indicated as M RB ∈ R 3×3 and M A ∈ R 3×3 , respectively M RB and M A are given as constant matrices here, and the inverse (M RB + M A ) −1 exists for the given matrices. Similarly, C RB and C A are the rigid-body and added Coriolis and centripetal matrices, respectively. The damping effects are modeled by a D L and a nonlinear part D NL . The added and damping matrices are due to the hydrodynamic effects of marine crafts moving in fluid. The detailed matrix forms and paramter values can be found in the Appendix A. Physical limitations are imposed on AUV speeds, and propeller and rudder working ranges as:

Problem Statement
The AUV is required to start from an origin and track a time parameterized path to arrive at the destination. In the specified origin-destination region where the AUV operates, the ocean current information, i.e., V c and β c , usually can be either predicted by regional ocean models or roughly measured by the sensors carried with the AUV [10]. However, these predictions or measurements are not precise and are mostly stochastically uncertain. Therefore, we assume that the ocean current could in principle be obtained by predictions, but with small prediction biased errors b and uncertainties following a normal distribution as, i.e., b ∼ , andV c ,β c are the predicted mean current speed and angle in the inertial coordinate system, respectively; is the corresponding covariance that reflect the prediction accuracy where diag means diagonal matrices; b = [ v , β ] with v and β are the corresponding prediction biased errors, respectively. Note that the biased errors are unknown in the design process, and are modeled only for simulation purposes. In the simulations, we will explore how the ocean current biased errors influence the AUV motions as well as the robustness to small magnitude biased errors due to the feedback from acoustic localization. Furthermore, V c and β c are assumed to be independent and identically distributed random variables.
For the controller design, full states need to be known precisely for feedback in general. However, the sensors such as GPS that are used for the state measurement of ground or airborne vehicles are not applicable in the underwater environment. Therefore, in the next section, we first propose a distance dependent localization method to "measure" the AUV position.

Acoustic Localization
Taking advantage of the finite propagation speed of sound in water and the availability of GPS or global navigation satellite system (GNSS) signals [26] over the ocean surface, we deploy three surface buoys with known positions [x n , y n , z n ] T , n = 1, 2, 3, as shown in Figure 2. However, the underwater velocity of sound is generally not constant but varies with the water pressure, temperature, and conductivity. Therefore, the propagation path of the sound is actually curved, which leads to non-random deviations of signal arrival times. To overcome this difficulty, the ray tracing technique [27] with isogradient sound speed profiles (SSP) is utilized in the AUV acoustic localization problem.
Consider the scenario that each of the three buoys carries a pinger and emits acoustic signals every T l s seconds. Denote the position of AUV as [x(t), y(t), z(t)] T at time t. Since the depth of underwater vehicles can be accurately measured with a depth sensor [28], we consider the AUV to move at a known constant depth of z meters. The clocks of the three buoys and the AUV can be synchronized with either GPS or lower power atomic clocks. For the ray tracing problem between the surface buoys and the AUV, the SSP is dependent on the water depth as: where c is the sound speed, a is a constant depending on the underwater environment, and b is the sound speed at the surface. Then, the TOAs of the signals from the three buoys at the AUV, denoted as ∆t n , n = 1, 2, 3, can be calculated in simulations as [27]: Y n (t) = L n (t)/X n (t) (11) α n (t) = arctan(Y n (t)) (12) β n (t) = arctan(X n (t)) (13) θ B n (t) = β n (t) + α n (t) (14) θ A n (t) = β n (t) − α n (t) (15) ∆t where d n stands for the actual horizontal distance between buoy n and the AUV; L n , X n and Y n are the intermediate variables; α n is the angle between the actual sound ray and the straight line path for buoy n; β n is the angle between the horizontal direction and the straight line path for buoy n. Lastly, θ B n and θ A n are the grazing ray angles at buoy n and the AUV, respectively. In simulations, TOAs ∆t n (t), n = 1, 2, 3 can be calculated as Equation (16), and in practice, they can be measured with synchronized clocks. Then, with the known TOAs, AUV depth z, and the fixed buoy locations, the horizontal distanced n could be estimated reversely using algorithms such as the binary search [27]. The main idea is to compare the arrival time of the median point with the actual propagation time repeatedly so that the search area could be narrowed down until a precise enough horizontal distance estimatê d n is found.
Next, we are able to obtain the position of the AUV withd n . If in one signal emitreceive cycle, the AUV is moving slowly and the position is assumed to be the same when the three signals from buoys 1, 2, and 3 are received, then we can get a set of nonlinear equations as: by solving which we get the acoustically estimated position (x(t),ŷ(t)) of the AUV. Note that in practice, only two equations are necessary to be solved for (x(t),ŷ(t)). Available functions such as fsolve using trust-region algorithm in Matlab [29] can be applied to solve Equation (17).
If the AUV horizontal heading is measured by on board sensors such as the inertial measurement unit asψ(t), the AUV horizontal pose can be obtained as: which can be seen as the "measured" system output. In the underwater environment, signal-to-noise ratio (SNR) is usually used to describe the quality of the received signals.
High SNR of the received signals means more accurate TOAs, which leads to high localization accuracy. Therefore, due to the uncertainties in measured TOAs, the surface buoy positions [x n , y n , z n ] T , and the yaw angle measurement system, the "measured" AUV outputŷ(t) is assumed to be subject to random noises ε ∼ N(0, Σ) where ε is an R 3×1 vector, and the covariance with σ x , σ y , and σ ψ being the corresponding standard deviations. (1) and (2) with state feedback information from acoustic localization, we propose the trajectory tracking controller with MPC and EKF. MPC uses a system prediction model to obtain the system trajectory over a future finite time interval. Optimization problems are formulated based on the predicted trajectories and are solved repetitively online in a receding horizon way, and thus also called receding horizon control. To relieve the possible computational burden using Equations (1) and (2) as the prediction models, successive linearizations are applied in the framework of MPC. Moreover, an EKF approach is proposed to handle the uncertainties in the system and acoustic measurements. Corrected system state information is obtained and used in the predictive tracking controller design. During the control time steps when the updated localization information is not available, the dead-reckoning technique is utilized.

Successively Linearized Prediction Models
The highly nonlinear and uncertain AUV dynamics in Equations (1) and (2) render the online optimization based MPC control difficult. The conventional one step linearization approach was proposed for autonomous surface vessels in [30]. The computational time is shown to be much shorter than using the nonlinear prediction models directly. However, the one step linearized model could lead to large linearization errors over the whole prediction horizon and thus degrade the control performance. Therefore, the successive linearization approach proposed in [31] is applied to achieve a trade-off between computational complexity and control performance. The basic idea is to utilize the whole sequence of control inputs from a previous time step and pre-calculate a shifted system trajectory for linearizations over all prediction steps. For completeness while avoiding repetitions, we briefly describe the implemented successive linearization algorithm for the uncertain AUV dynamics.
Generalize the AUV dynamics Equations (1) and (2) as: where f : is the control input vector; b is defined as before. For numerical simulations, continuous time model Equation (19) is discretized with zero-order-hold assumption as: where time step k relates with the real time t as t = kT c s with control sampling time T c s . Note that the control sampling time T c s is not necessarily the same with the localization sampling time T l s . At each control time step k, successive linearizations are primarily implemented as: 1.
Shift the optimal control input sequence from the previous step, u * (i|k − 1), i = 0, 1, . . . , N p − 1 with N p being the discrete horizon steps, to obtain a seed control input trajectory u 0 (i|k):

2.
Apply u 0 (i|k) to Equation (20) and obtain a seed state trajectory x 0 (i|k), i = 0, 1, . . . , N p . The initial state is set as x 0 (0|k) = x(k) with x(k) being the current system state.
The seed disturbance values are set as whereb(i|k) is constant as [V c ,β c ] for simplicity without loss of generality. The proposed algorithm is still applicable and works effectively when time-varying values for [V c ,β c ] are available. The EKF that will be introduced in the next subsection is also based on linearizations at the mean values of random disturbances.

3.
Linearize the nonlinear dynamics at the seed trajectories where are small perturbations around x 0 (i|k), u 0 (i|k) and b 0 (i|k), respectively; A(i|k), B(i|k), and E(i|k) are the corresponding Jacobian matrices, respectively, of the nonlinear dynamics function Equation (20) evaluated at x 0 (i|k), u 0 (i|k), b 0 (i|k) .

EKF with Acoustic Localization Information
Kalman filter is a widely used correct-predict formulation to deal with velocity estimation and stochastic environmental disturbances [15]. EKF further extends the system dynamics to nonlinear types by linearizations to propagate an approximation of the state estimation and covariance. With the models Equations (18), (20) and (22), we are ready for the EKF design.
Assuming that T l s = nT c s and n ≥ 2 is an integer. Then, at control time step k while the new acoustic localization information is available, i.e., the localization step k/n, the AUV predicted state is corrected bŷ where is the output matrix;x(k) is the predicted state by ignoring the uncertainties in the nonlinear dynamics Equation (20), i.e.,x where f d is the discrete function that generalizes Equation (20). Note the the parameters with a horizontal bar overhead such asb andx all indicate the predicted system information. The Kalman gain matrix K(k) is derived by The covariance matrixP(k) for the prediction error x(k) −x(k) is computed by the recursion where P(k) is the covariance matrix of the estimation error x(k) −x(k). Through this standard Kalman filter procedure, the predicted system statex(k) can be corrected by the acoustically "measured" outputs y(k) whenever this information is available. The corrected statesx(k) can then be used for feedback design in MPC. However, as mentioned earlier, the control cycle is generally shorter than the localization cycle, i.e., T c s ≤ T l s , which means y(k) is not guaranteed to be there for each control time step k. When the acoustic measurement is not available, the predicted statex(k) from Equation (28) is used for feedback, which is commonly known as the dead-reckoning technique [14].

Predictive Tracking Control
Briefly, at each time step k, MPC computes a future control sequence over the prediction horizon N p based on the current system state. The predicted system trajectory is then driven close to the given reference. Only the first element of the control sequence is applied to the system. At the next time step k + 1, this process is repeated over the shifted prediction horizon with newly measured system states. MPC is based on solving optimization problems online repetitively and mostly numerically, which naturally and explicitly handles system constraints. Therefore, MPC is especially suitable to systems of which analytical control laws are difficult to derive, and suitable to systems with constraints. For the AUV tracking problem, system constraints Equations (4)-(6) are imposed. We next propose the MPC approach to solve the AUV tracking problem with the corrected information from EKF or the predicted information from dead-reckoning.
Consider at time step k, given the current estimatedx(k) when the acoustic estimated statex(k) is not available, future predicted trajectories x(i|k) over N p are then made following Equation (22). Moreover, given a reference path determined by an origin (x o , y o ) and a destination (x d , y d ) in the inertial coordinate system, a constant reference speed u d is set to parameterize path. The reference trajectory over a future prediction horizon N p at k is calculated as where ψ d = arctan((y d − y 0 )/(x d − x 0 )). The tracking problem is then represented by solving a constrained finite horizon linear programming problem subject to Equation (22) and u min u(i|k) u max , i = 0, 1, · · · , N p , (34) whereȳ(i|k) = Cx(i|k). The cost function Equation (32) aims to minimize the tracking error between the predicted state trajectory and the reference trajectory with optimal control efforts. The three terms are the terminal tracking error cost, the stage tracking error cost, and the control input cost, respectively; Q f , Q and R are positive definite weighting matrices for the respective cost terms and can be adjusted to set priority between convergence speed and control efforts. Due to the constraints on system states and control inputs, the optimization problem does not have an analytical solution, and needs to be solved numerically repetitively online with newly "measured" system states and time-varying prediction models. However, due to the successive linearization, the above optimization problem becomes convex and can be solved efficiently.
Overall, the predictive tracking control algorithm with EKF and acoustic localization information or dead-reckoning can be summarized as Algorithm 1.

Simulation Results and Discussions
Simulations are run to demonstrate the effectiveness of the proposed tracking control algorithms for AUVs. The proposed integrated localization and control algorithm is applied in following a reference path. Such applications can be seen in sea mapping scenarios where only waypoints in the concerned ocean region are given. The AUV is positioned at the initial point [100, 100, 300] m and is controlled to reach the destination [600, 600, 300] m along the following curved reference path x d (t) = 600 − 500 cos(πt/1200), y d (t) = 100 + 500 sin(πt/1200), for t ∈ [0, 600] s. The three surface buoys are positioned at [0, 0, 0] m, [500, 800, 0] m, and [1000, 0, 0] m, respectively. The isogradient SSP parameters in (7) are set as: a = −0.05 and b = 1540 m/s. In this region, ocean currents exist. The average speed is predicted to be 0.5 m/s and the angle is predicted to be π/6. However, there might be prediction biased errors. In the simulations, we will first validate the performance of the proposed integrated localization and tracking control algorithm with zero prediction biased errors, i.e., v = 0 and β = 0. Then, to illustrate the inherent robustness (to a certain degree) of the feedback due to acoustic localization, we also consider non-zero prediction biased errors with v = 0.15 m/s and β = π/20 rad. Note that the biased errors are not considered explicitly in the algorithm design, and are only implemented for simulation purposes to show the inherent robustness due to the feedback from the acoustic localization. Moreover, there exist uncertainties with the predicted values following the normal distribution with covariance Σ b = diag([0.1 m 2 /s 2 0.01 rad 2 ]), which means low prediction accuracy in the ocean current. It should be noted that AUVs with any design may fail at those worst predictions of the sea. The uncertain acoustic measurement is with covariance Σ = diag([0.1 m 2 0.1 m 2 0.01 rad 2 ]). The AUV depth is 300 m, and the buoy depths are all 0 m. Controller parameters are set as follows: prediction horizon N p = 10; weight parameter Q f = Q = diag([100 100]) and R = diag([1 180/π]); sampling time T c s = 0.5 s and T l s = 1 s, which means the system feedback states are corrected by the acoustic localization results every two control time steps, i.e., n = 2 in Algorithm 1. Physical system constraints are: u max = −u min = 2 m/s, ξ max = −ξ min = 86 N, and δ max = −δ min = 13.6π/180 rad. The REMUS AUV system parameter values are referred to [25]. Note that the practical parameter values might be different from the above set values. For adverse environmental situations, the algorithm might not be applicable. All the algorithms are implemented in MATLAB 2016b [29] with solver Cplex [32] on a platform with Intel(R) Core(TM) i3-7100 CPU @ 3.70 GHz.

Localization Results
For the considered scenario, the three buoys are localized at the surface as shown in Figure 3a. The TOAs of the signals from the three buoys calculated by Equations (8)- (16) considering the isogradient SSP Equation (7) are plotted in Figure 3b. Since T l s = 1 s, there are a total 600 TOAs being calculated for each buoy. As the AUV moves from the initial position to the destination, the distance between the AUV and buoy 1 is increasing while the distances between the AUV and buoys 2 and 3 are decreasing. Therefore, it is observed in Figure 3b that the TOA trajectory of buoy 1 ascends while the other two trajectories descend. Note that since the AUV approaches towards buoy 2 faster than buoy 3, the TOA trajectory of buoy 2 declines also more rapidly than that of buoy 3.    Figure 4 further shows the distance errors between the real horizontal distances d i (t) and the estimated horizontal distancesd i (t) as calculated with a binary search algorithm [27]. Small distance errors are observed. The micrometer magnitude in Figure 4 is because the threshold specified in the binary search algorithm for the estimated horizontal distance is set to 1 × 10 −5 . The binary search algorithm makes trade-offs in terms of speed and accuracy.

Tracking Control Results without Prediction Biased Errors
For the MPC controller design, system state feedback is necessary at each time step as in Equation (33). However, due to the uncertainties in system dynamics and measurements, only estimated system states from EKF or predicted system states from dead-reckoning with nominal system dynamics can be utilized. Figure 5 shows that the proposed algorithm can achieve the trajectory tracking goal and control the AUV moving from the initial position to the destination along the path. In Figure 5, five trajectories are plotted in total:

1.
The reference path is determined by connecting the initial position and the destination; 2.
The real AUV trajectory is x(t) recorded from the simulations by evolving Equation (20) based on Equations (1) and (2) with known uncertain ocean currents. 3.
The estimated AUV trajectoryx(t) is the trajectory from the EKF, i.e., from Equation (26). This trajectory is obtained by correcting the dead-reckoning trajectory using the acoustic localization feedback information, i.e., the "measured" output,ŷ(t).

5.
The measured outputs y(t) are obtained from Equation (18) by adding stochastic uncertainties at the computed results from the TOAs according to Equation (17) at each localization step.
Overall, the differences among these trajectories are small. Therefore, a zoom-in around position (300, 300) m is added. The root-mean-square-error for tracking the reference path is 5.59 m which is qualified in applications such as large-area ocean mapping. Note that the estimated trajectory and the measured trajectory are plotted with fewer position points than the dead-reckoning trajectory since the control cycle is faster than the acoustic localization cycle. However, for the speed trajectory, the dead-reckoning speed trajectory satisfies the constraints well while small violations of the real and estimated trajectories with respect to the speed limits are observed. This is because both the optimization problem Equation (32) subject to Equations (22), (33)-(36) and the dead-reckoning trajectory use the nominal system dynamics without uncertainties. If all the constraints Equations (22), (33)-(36) are satisfied in the optimization problems, the dead-reckoning trajectory will also satisfy the constraints. However, compared to the dead-reckoning trajectory, the real trajectory further incorporates the uncertainties in ocean currents, and the estimated trajectory further incorporates the uncertainties in measurements. These uncertainties push the dead-reckoning trajectories that are already on the constraint boundary out of the feasible zone. Approaches such as buffer zone design in [31] can be applied to deal with this issue. Note that the differences between the estimated and dead-reckoning speed trajectories are small, and the root-mean-square-root value between these two trajectories is 0.03 m/s. This difference demonstrates the correction capability of the EKF as in Equation (26). To relieve the potential computational burden of the MPC online optimization problems, successively linearized prediction models Equation (22) are used. Therefore, only linear programming online problems need to be solved. Figure 8 shows the solver times of the proposed tracking controller. Mostly throughout the simulation, the solver times are between 10 ms to 20 ms, and are smaller than the system control sampling time T c s = 0.5 s, which indicates the potential of the proposed algorithm for real-time applications.

Tracking Control Results with Prediction Biased Errors
When there are biased errors in the predicted ocean current information, due to the inherent robustness of the MPC feedback control by using the updated states from acoustic localization, the AUV can still follow the reference path, as shown in Figure 9. However, compared to Figure 5 where there are no prediction biased errors, the AUV trajectory sees larger fluctuations. The root-mean-square-error for the reference tracking increases to 6.17 m compared to 5.59 m in the case of no biased errors. The larger tracking errors are partly due to the modeling errors, i.e., the error between the model used for the controller design and the model used for simulations. Furthermore, we also run a simulation where no ocean current information is used in the controller design, i.e., the predicted mean current speed, angle, and their covariances are all set to zero in the controller. In this case, infeasibilities occur and the simulation could not continue. Therefore, a proper consideration of the influence of ocean currents is necessary.

Tracking Control without Acoustic Localization
In order to demonstrate the role of acoustic localization for the AUV tracking control, we also run a set of simulations without acoustic localization. The system states used for MPC feedback are purely provided by the dead-reckoning module. This is the case when the AUV is not equipped with any positioning devices in an underwater environment. In this case, the AUV can never even reach the destination. Therefore, we set the simulation time the same as in the case with acoustic localization, and the tracking results are shown in Figure 10. It can be observed that the trajectory deviates far away from the reference path. The tracking task fails which means the tracking control algorithm with purely deadreckoning is not acceptable in practical applications. Compared with Figure 5, Figure 10 further demonstrates the effectiveness of the proposed integrated acoustic localization and tracking control approach. Note that since no acoustic localization feedback is used, there is no measured trajectory and the estimated trajectory coincides with the dead-reckoning trajectory. Moreover, since now the controller is purely based on the dead-reckoning states which are subject to no disturbances, the dead-reckoning trajectory can follow the reference path well.

Conclusions and Future Work
For the two challenges that face the autonomy of AUVs, i.e., localization and ocean currents, this paper proposes a systematic approach. Firstly, the localization problem is solved with the aid of three surface buoys. By leveraging the different latencies of the three surface signals arriving at the AUV, the position of AUV can be calculated every localization cycle. The acoustically calculated position can then be utilized to correct the predicted and possible inaccurate position in dead-reckoning with EKF. Moreover, the uncertain ocean currents and full state reconstruction can also be incorporated in the framework of EKF. An MPC controller is designed to handle system constraints and control metrics based on the estimated full states from EKF. Simulation results demonstrate that the proposed integrated acoustic localization and tracking control approach outperforms the tracking control without acoustic localization. The potential of applying the proposed algorithm to AUVs is illustrated especially in addressing the challenges of underwater localization and ocean currents. Future research will consider extending the AUV dynamics to three dimensional space and the possibility of cooperative localization and control with multiple AUVs.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
The involved REMUS AUV model matrices and parameter values are presented as follows [14,25].
The rigid-body and added mass matrices are The rigid-body and added Coriolis and centripetal matrices are Hydrodynamic derivatives follow the SNAME nomenclature. For instance, the hydrodynamic added mass force X along the x axis due to an accelerationu in the x direction is written as which implies {M A } 11 = −Xu. The REMUS AUV parameter values are listed in Table A1.