{'fraction_non_alphanumeric': 0.02574396755101196, 'fraction_numerical': 0.008940027316750134, 'mean_word_length': 4.769340974212034, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic System Identification of Underwater Vehicles Using Multi-Output Gaussian Processes

Non-parametric system identification with Gaussian Processes for underwater vehicles is explored in this research with the purpose of modelling autonomous underwater vehicle (AUV) dynamics with low amount of data. Multi-output Gaussian processes and its aptitude to model the dynamic system of an underactuated AUV without losing the relationships between tied outputs is used. The simulation of a first-principles model of a Remus 100 AUV is employed to capture data for the training and validation of the multi-output Gaussian processes. The metric and required procedure to carry out multi-output Gaussian processes for AUV with 6 degrees of freedom (DoF) is also shown in this paper. Multi-output Gaussian processes are compared with the popular technique of recurrent neural network show that Multi-output Gaussian processes manage to surpass RNN for non-parametric dynamic system identification in underwater vehicles with highly coupled DoF with the added benefit of providing a measurement of confidence.


Introduction
Dynamic modelling of unmanned underwater vehicles (UUVs) has been a subject of interest among researchers since the early days of underwater exploration. Nowadays, UUVs are extensively employed in research, industry and military applications. Modelling of autonomous underwater vehicles (AUVs) is an important step for design of mission, control and navigation systems. Thus, accurate modelling and adaptability of such systems is an important issue due to such extensive applications. The most common methodologies use a mathematical model which is derived from Newtonian-Lagrange mechanics. This mathematical model is composed of a series of coefficients that need to be calculated to obtain an accurate model. Difference between the obtained model and the reality is usually treated in the literature as noise and in most cases, is modelled as a Gaussian distribution. A Gaussian distribution can be extended arXiv:2006.02194v1 [eess.SY] 3 Jun 2020 to the calculation of an approximation to a real model of a vehicle with higher exactitude and adaptability than a mathematical model [1].
Over the years, multiple methods for the calculation of coefficients of underwater vehicles have been proposed. One way to obtain the hydrodynamic coefficients is to perform a series of captive model tests such as rotating arm and planar motion mechanics [2,3,4]. Another common technique for the hydrodynamic coefficient calculation is the use of computational fluid dynamics (CFD) [5]. However, for the successful application, CFD still requires verification of results with experiments [6].
Nevertheless, research has probed the variability of mathematical models for AUVs as the vehicle operates in proximity to objects [7], near surface [8] and most commercial available underwater vehicles with a modular architecture involving variable geometric and mass. Furthermore, in certain applications, the precision of some coefficients is require to be within 5% of accuracy [9]. Therefore, the variability of coefficients and the high precision required make it cumbersome or even impossible to acquire an exact analytical system model based on physical rules.
Another procedure to obtain coefficients from a model of an underwater vehicle is the use of observers. Common observers applied to obtain the hydrodynamic coefficients of AUVs from measured data are least-squares [10,11], nonlinear Kalman filters such as extended (EKF) and unscented Kalman filter (UKF) [12]. The EKF requires the linearization at each time step for the approximation of non-linearities which can be difficult to regulate and implement. A method to overcome this is the use of UKF which applies the unscented transform over a set of methodically chosen samples to model the system nonlinearity [13]. Other common methodologies are frequency domain identification [14]), neural networks (NN) [15] and support vector machines (SVM) [16]. The latter two methodologies are machine-learning algorithms and are more commonly used in online learning of the coefficients and also provide system adaptability.
The adaptation of mathematical model has inherent defects such as the dependency of initial values, small quantity of coefficients to be updated, ill-conditioned matrix and drift.
Machine learning algorithms are not limited to the calculation of hydrodynamic coefficients as they can learn to behave as part of the system or the complete system. Multiple applications have taken advantage of this ability and used NN [17] and SVM [18] to learn the damping model for the system which is placed in parallel to a well-known partial mathematical model. Other applications have used pure machine learning algorithms to identify a complete underwater vehicle as a black-box model with the use of nonlinear autoregressive model with exogenous (NARX) architecture. [19] have used multiple architectures of NN for the regression of an AUV model in a NARX architecture and used the learned model for model predictive control. Their study shows that recurrent NN (RNN) provides higher faithfulness to the plant.
Recently, [20] compared different machine learning algorithms for the system regression of underwater vehicles, i.e. NN, SVM, Gaussian Process Regression (GPR) and Kernel Ridge regression (KRR). Their results show that the machine learning algorithms could model an AUV from onboard sensor data in comparison to a least squares approach. Nevertheless, in their study, a structure for dynamic system identification has not been employed and each degree of freedom was treated as separate element. This can be problematic in AUVs as the outputs are strongly coupled. In the specific case of modelling with Gaussian Processes (GPs) [1], research shows that the dynamic regression of a system with GPs can produce better results than other methodologies. The most common methodology for Multi-Input-Multi-Output systems is to model each DoF as a separate system [21]. More advance methodologies for Dynamic system identification have been proposed in [22] and [23] and specific methodologies are introduced for the identification of multi-output GPs based on the use of variation of dependent GPs.
GPs is a well-established methodology in fields such as geostatistics, where the method is called 'kriging' [24]. In GPs-based system identification, the model is built over input-output data and a covariance function is used to characterize the vehicle behaviour. The advantage of GPs is their ability to work with small quantities of data and noisy data. The predicted results consist of a mean and variance value which can be used for other purposes as well such as control, navigation and model based fault detection as it contains a measure of confidence.
Multi-output GPs are a special case of GPs with the capability to model the nonlinear behaviour and coupling among outputs of a multi-output system [25], both which are important for AUV dynamics. In this study, a non-parametric dynamic system identification with Multi-Output GPs architecture employed by the authors for ships [26] is extended to AUVs. The output from the algorithm will be a predictive value and a measure of confidence of the predictive value. The present implementation was made over data obtained from a non-conventional test with variable frequency of a nonlinear simulation model of a REMUS 100 AUV. Multiple sample times and data length were tested to find the best metric that can describe an AUV. A RNN was employed as a comparison to measure the effectiveness of the proposed method.

Nonlinear Dynamic AUV Model
In [27] it was shown that the nonlinear dynamic equations of motion of an underwater vehicle can be expressed in vector notation defined by a state vector composed by the vector v of velocities on the body frame of the form [u, v, w, p, q, r] T and the vector η of position in the Earth fixed frame (Figure 1) of the form [ξ, η, ζ, φ, θ, ψ] T such that with the kinematic equationη where η position and orientation of the vehicle in Earth-fixed frame,  Equation (1) can be expanded into a more general equation of motion as has been shown in [28,29].The result of the expansion will be a system of six equation with 73 hydrodynamic coefficients. However, for a complete model the control surfaces must be modelled. In a general case, the resulting forces and moments of a control surface (thrusters and fins) can be expressed as [27] A more accurate thruster model can be found in [30] with the inclusion of the motor model and fluid dynamics. However, in this study, a more conservative model is used. Details of the Remus 100 AUV model used in this study are given in Section III.

Dynamic Identification with Multi-output GPs
GPs can be defined as a generalization of a multivariate Gaussian distribution. A multivariate Gaussian distribution is defined by it's mean and a covariance matrix. In the case of GPs is a distribution over functions rather a distribution over vectors. GPs is one of the methods based on kernel functions where the kernel function calculates the relationship between an input and an output point, and generates the covariance between them. The covariance determines how strongly linked (correlated) these two points are. In the case of multi-output GPs this is extended by the convolution of kernels to add not only the relationship between an input and an output but also the relationship between the outputs. The kernel is the key ingredient for the calculation of the covariance matrix that correlates inputs and outputs of training data [31].
The design of the algorithm for multi-output system identification with GPs is based on the previous work of Alvarez and Lawrence [32] and [1]. The dynamic identification problem can be defined as the search for relationship between a vector formed by delayed samples from the inputs u(k − 1) and outputs y(k − 1) and the future output values. The relationship can be expressed by the equation: where f (x(k), Θ) is a function that maps the sample data vector x(k) to the output space based on the hyperparameters Θ; v(k) accounts for the noise and error in the prediction of output y(k). In the case of dynamic system identification the discrete time variable (k) is presented as an embedded element in the regression process as it is accounted in the delayed samples.
A requirement for dynamic system identification of nonlinear systems is the selection of a nonlinear model structure such as nonlinear autoregressive model with exogenous input (NARX), nonlinear autoregressive (NAR), nonlinear output-error (NOE), nonlinear finite-impulse response (NFIR) and other structures. From all the possible structures, the simpler and most popular structure to implement is NARX as it only requires measurements of system output/s and input/s. In the case of an AUV, NARX is the most practical configuration since the measuring points are restricted to the available sensors [1].

Multi-output GPs
The non-linear dynamic system of an AUV (Eq. (1)) shows the level of coupling between the Newton-Lagrange equations of an AUV. The nonlinearity and coupling between outputs can be better represented by a multi-output GPs.
Multi-output GPs presented here are based on the work of [32]. Multi-output GPs are founded in the regression of data using convolving white noise process with a smoothing kernel function [33]. This was later introduced by [34] to the machine learning community by assuming multiple latent process defined over a space q . The dependency between two outputs is the model with a common latent process and their independency with a latent function, which does not interact with other outputs. If a set of functions {f q (x)} Q q=1 is considered, where Q is the Output Dimension for a N number of data points, where each function is expressed as the convolution between a smoothing kernel {k q (x)} Q q=1 and a latent function u(z) , This equation can be generalized for more than one latent function {u r (x)} R r=1 and includes a corruption function (noise) independent to each of the outputs The covariance between two different functions y q (x)and y s (x ) is: where If it is assumed that u r (z) is an independent white noise cov [u r (z) , u p (z )] = σ 2 ur δ rp δ z,z will become: The mean y with variance σ y of a predictive distribution at the point x given the hyperparameters Θ can be defined as A comprehensive description and implementation of the convolution process can be found in [32] and [35] respectively. In this study, the convolution of two square exponential kernels are used since squared exponential kernel is a universal kernel [36], provided that data is stationary and the function to be modelled is a smooth one. Furthermore, squared exponential kernel has small quantity of hyperparameters to be established.

Learning Hyperparameters
There are two principal methods for learning the hyperparameters Θ namely: Bayesian model interference and marginal likelihood. Bayesian inference assumes that prior data of the unknown function to be mapped are known and a posterior distribution over the function is refined by incorporation of observations. The marginal likelihood method is based on the aspect that some hyperparameters are going to be more noticeable. Over this base, the posterior distribution of hyperparameters can be described with a unimodal narrow Gaussian distribution.
The learning of GPs hyperparameters Θ is normally carried out with the maximization of marginal likelihood. The marginal likelihood can be expressed as: where K is the covariance matrix, N is the number of input learning data points and y is a vector of learning output data of the form [y 1 ; y 2 ; · · · y N ] . To reduce the calculation complexity, it is preferred to use the logarithmical marginal likelihood that is obtained by the application of logarithmic properties to Eq. (14).
In order to find out a solution for the maximization of log-likelihood, there are multiple optimization methods that can be used such as particle swarm optimization, genetic algorithms, or gradient descent. For deterministic optimization methods, the computation of likelihood partial derivatives with respect to each hyperparameter is required. From [37] log-likelihood derivatives for each hyperparameter can be calculated by: Equation (14) gives us the learning process computational complexity, for each cycle the inverse of the covariance matrix of K has to be calculated. This calculation carries a complexity O(NM) 3 . After learning, the complexity of predicting the value y(k + 1) is O (NM) and to predict the mean value σ(k + 1) is O(NM) 2 . The higher order term O(NM) 3 is the major disadvantage of using multi-output GPs. If the number of data increases the complexity of learning the hyperparameters increases in a cubic form. Methods such as genetic algorithms, differential equations, and particle swarm optimization can be applied to avoid the calculation of the marginal likelihood partial derivatives and thereby reduce the computational time.

Experiment setup
The implementation of a mathematical model of an underactuated REMUS 100 AUV was used to generate the required identification data. The coefficients of [28] were used and adapted for simulation on Simulink with the addition of the thruster model from [38] as the original mathematical model produced by [28] has a constant thrust force. The resultant model was tested to mimic the original results obtained by [28] at a speed of 1.5 m/s. The AUV details can be found in Table 1. As shown in Fig. 2 a simulation setup was developed in MATLAB/Simulink to emulate the AUV behaviour. Figure 3 shows an example of input signals given for the rudder angle, thruster RPMs and elevator angle respectively. A total of 8 sets of data were produced by combining and initial chirp signal and after the first 1000 seconds the command signal change to a step function or ramp function. Simulation was carried out for 2000 seconds. A list of simulation done can be seen in Table 2. The objective of not using a standard test such as zigzag test or turning circle test is to test the ability of GPs under more drastic conditions. A sample data point was captured for each 1.5 seconds over the input and outputs. A total of 8000 points were captured over the six motion outputs and 4000 point over the three input signals (Propeller RPM, rudder angle and elevator angle). The data set was divided into two sets of points, the first set of points is used for the model learning, this data of training is equivalent to the section of chirp input signal, and the second set of points is used for learning validation. The validation data is purposely chosen to be beyond the range of training data and very different to the training data to test the ability of the method to predict beyond the training range.  Chirp+Ramp in All surfaces 5 Chirp+Step in Propeller 6 Chirp+Step in Rudder 7 Chirp+Step in Stern 8 Chirp+Step in All surfaces

Training and validation
A script was written to implement the NARX structure with multi-output GPs and the implementation of multi-out GPs by [35] was employed. The multi-output GPs regressors were defined as: Where y is the vector of regressors [u k−1:3 , v k−1:3 , w k−1:3 , p k−1:3 , q k−1:3 , r k−1:3 ] T , the function f is a relation between the vector of regressors from the correspondent vehicle speeds (u, v, w, p, q, r) or the full vehicle state y , and the vector cthat content the regressors of the commanded signals (u rpm , u elevator , u rudder ) to the respective output of the system. The input signals were normalized between -1 and 1 to give all the inputs and outputs the same weight in the learning process.
For the training, a minimum search with the gradient descent method, in particular the interior-point algorithm, was used for the minimization of the negative logarithmical likelihood.  three terms of delays for the output to be feedback to the network and three terms delay of the inputs. The RNN 1 that was selected as relatively optimal for the task at hand used two hidden layers with logarithmic sigmoid functions for the hidden layer neurons and was trained with Levenberg-Marquardt backpropagation, this configuration and regressors was selected as provide the best results for our system. RNN 2( Figure 5) is a fully connected RNN with a single hidden layer. RNN 2 use the same logarithmic sigmoid functions for the hidden layer neurons and was trained with Levenberg-Marquardt backpropagation. A third step of simulation was carry out with the combination the full length of the data and feeding back after each step the output from the models. The neural network system was trained, validated, and tested with the same data used for the multi-output GPs. The complete implementation code can be found at the GitHub Repository 1 . Figure 6 presents the results of GPs training compared to the AUV simulator signals, and the error plots between the predicted and real systems. In all the figures, a 2σ variance is plotted. The variance values of the training data are in the expected value and encompass the error results.
The validation data consisted of the second part of the captured data in the vector form with the commanded inputs and the real output from the training data with the respective system delays. The segments of results from the validation with the second set of data are depicted in Fig. 7. The low validation errors show a good system prediction for the sway speed and yaw speed.

Simulation
With the objective to test the ability of the learning system, a simulation stage was implemented for the learned RNN and the multi-output Gaussian processes over the total length of the simulated data. Navigation applications as EKF, UKF, and control as model predictive control requires predicting the behaviour of the plant in a number steps ahead of the actual state to predict the correct position or control signals. In the case of the multi-output GPs the simulation is done by feeding back to the simulation the past inputs y i (k − n) , the initial position and control signals of rudder, elevator and forward speed where used, the naive simulation [1] covers training and validation data acquired from the original simulation in a close loop setup. Naive simulation provides an approximation where variances are not exactly the same, but provide a general guidance on uncertainty and this is enough for our study. If the variance is to be employed, such as in a control systems or a navigation problem, the uncertainty propagation can be included with the use a simulation based on Monte Carlo numerical approximation [1]. Figure 8 shows the results from the simulation of RNN 1, RNN 2 and the multi-output GPs compare to the original system. RNN 1 compare to RNN 2 show better performance with ramp input signals, and RNN 2 show better perfomance to simulate step functions in our simulations. However multi-output GPs can identify the system correctly and predict the behaviour of the system with chirp+ramp and chirp+step functions . The better capability of GPs to predict outside the training horizon from a number of difference variation from the initial training data is confirmed by the results of Tables 3 to 5. The mean value of the output root-mean-square error (RMSE) ,the predicted residual error sum of squares (PRESS) and mean absolute error (MAE) for GPs are smaller than the values of RNN 1 and RNN 2. seconds were carried out with commands signals composed of a chirp signal for half of the time and a ramp signal for the other half of the simulation. The first half of each data set was employed for training and the test simulation of the learned model was done over the complete extension of data. The RMSE ,PRESS and MAE were also measured over all simulation results. The result of the sensitivity analysis can be seen in Figs. 9 and 10. All the measurement of the sensitivity analysis show the same trend for each measured variable, the simulations with 1500 seconds of capture data for RNN1 and GPs show similar average results and RNN 1 over 4000 seconds of simulation can overpass the ability of GPs to simulate the system with a chirp+ramp signal for all inputs. RNN 2 and RNN 1 require higher quantity of data of rich data to be effective in the simulation of AUV outside of the learning horizon in comparison to multi-output GPs.

Conclusion
In this work, the use of multi-output GPs for the system identification of AUV dynamics was tested on a REMUS 100 AUV. It was demonstrated that the nonparametric multi-output GPs can model an AUV as well as RNN with the added value of a confidence measurement. In the simulations, GPs show a better ability than RNN to predict and simulate the behaviour of an AUV. In some cases, GPs performed better than RNN outside of the training horizons with the error between the GPs and the real system being relatively low as the convolution process is equivalent to represent the system through a differential equation. The GPs model obtained also has a smaller number of hyperparameters compared to the large number of coefficients of a mathematical model. The results of the sensitivity analysis show that multi-output GPs perform better that RNN with low quantities of data. RNN required a higher spectrum of data to be able to approximate the behaviour of the vehicle outside of the training horizon.
To improve further the capability of prediction of the system, more recent suggested techniques for GPs such as Recurrent GPs can be used. The simulation of GPs can also be improved if techniques such as Montercarlo and Taylor series can take advantage of the variance to increase the horizon of cover manoeuvres and the prediction accuracy. The next work will be devoted to the development of application for navigation and control using the obtained model. As the real world is a noisy environment that can be better described with Gaussian distributions, the use of GPs can show better performance in specific tasks such as navigation and control.