Hydrodynamic modeling with grey-box method of a foil-like underwater vehicle

In this study, a dynamic modeling method for foil-like underwater vehicles is introduced and experimentally verified in different sea tests of the Hadal ARV. The dumping force of a foil-like underwater vehicle is sensitive to swing motion. Some foil-like underwater vehicles swing periodically when performing a free-fall dive task in experiments. Models using conventional modeling methods yield solutions with asymptotic stability, which cannot simulate the self-sustained swing motion. By improving the ridge regression optimization algorithm, a grey-box modeling method based on 378 viscous drag coefficients using the Taylor series expansion is proposed in this study. The method is optimized for over-fitting and convergence problems caused by large parameter matrices. Instead of the PMM test data, the unsteady computational fluid dynamics calculation results are used in modeling. The obtained model can better simulate the swing motion of the underwater vehicle. Simulation and experimental results show a good consistency in free-fall tests during sea trials, as well as a prediction of the dive speed in the swing state.


Introduction
Recently, the shape designs of underwater robots have been diversified. Many well-known deep-sea underwater vehicles, such as Odyssey IV (Cooney, 2009) from MIT, Sentry (Jakuba, 2003) from WHOI, DEEPSEA CHAL-LENGER (Li et al., 2013), QIAN LONG II (Xu et al., 2015) from Shenyang Institute of Automation et al., have adopted foil-like bodies. In contrast to an underwater vehicle with a body of revolution, the design and modeling of an underwater vehicle with a foil-like body is complicated because no reference model is used during the design process.
The modeling of an underwater vehicle is the most challenging problem in the design of an underwater vehicle. The modeling is based on the existing reference model or performed using a grey-box modeling method, which considers part of the model as a white-box and the other part a black-box. The modeling of rigid-body dynamics is performed using the white-box method, which is not discussed in this paper. Therefore, it is difficult to model fluid forces, especially viscous forces. The most challenging problem for a foil-like underwater vehicle is the lack of a suitable refer-ence model for viscous forces. Gertler and Hagen (1967) proposed a standard submarine motion equation which is later widely used by underwater vehicles with a body of revolution as a reference model. However, this model required a similar shape of the underwater vehicle to the SUB-OFF model (Cao et al., 2016). Hence, the model could not be applied to vehicles of other different shapes. The hydrodynamic model adopted by Sentry AUV (Jakuba, 2003) and REMUS AUV (Prestero, 2001) used the theory of slender body to fit the empirical coefficient to determine the structure. Both of them used the white-box method to build the whole dynamic model. This modeling method can model the foil-like vehicles and the most specially shaped vehicles. However, the model operation is complex and difficult because it requires deep hydrodynamic foundations and extensive engineering experience. The hydrodynamic model of remotely operated vehicles (ROVs) is modeled by firstorder coefficients or first-order and diagonal second-order coefficients (Maalouf et al., 2015). The model is asymptotically stable in all directions of forces, and it cannot simulate the self-sustained swing motion. Martin and Whitcomb (2014) realized an online system identification method for an ROV, obtained a high precision in the experiment, and compared the identification results of many kinds of combined coefficients. However, this method required expensive sensors because it involved an on-line identification of the mass properties. In this study, the grey-box modeling method is used for dynamic modeling in which the rigidbody dynamics with a definite structure employs the whitebox approach and the dumping force lacking the reference model employs the black-box approach. The CFD method is utilized to obtain the virtual experimental data, and the improved ridge regression algorithm is used to optimize the large number of viscous coefficients, which can solve the above-mentioned problems.
Commonly used viscous force coefficients are obtained using the Taylor series expansion (Prestero, 2001); therefore, there are multicollinearity problems between the coefficients (Yoon and Rhee, 2003). The influence of multicollinearity on the identification of viscous force parameters is rarely discussed in the literature of common hydrodynamic studies. Martin and Whitcomb (2014) used an improved least-squares method to identify coefficients, which diminished the effect of multicollinearity on the identification results. However, this method is highly integrated and cannot judge the quality of data at the same time. Yoon and Rhee (2003) used the ridge regression method to solve the problem of multicollinearity. This method provided data quality judgment for users. However, the number of parameters used in their work was small. As the number of parameters increases, the simple ridge regression method demonstrates a serious decline in the goodness of fit. This study proposes an improved ridge regression algorithm to solve the problem of decline in the goodness of fit.
In this study, the Hadal Autonomous and Remotely operated Vehicle (Hadal ARV) is used as an example to verify the method. This ARV is a full-ocean-depth underwater vehicle developed by Shenyang Institute of Automation of the Chinese Academy of Sciences. The vehicle has a foillike shape as shown in Fig. 1. Owing to the work depth (~11000 m), the robot spends most of the time descending or ascending. The diving process shows a complex movement form. In the free-fall diving, the vehicle spontaneously exhibits the movement of a spiral swing. During the spiral swing motion, its dive speed greatly deviates from the resistance calculation results in the vertical direction. When the swing amplitude of the vehicle reaches its maximum, the attack angle increases. As a result, the resistance of water will increase rapidly. When the swing reaches the equilibrium point, the attack angle decreases and the vehicle's speed increases. In this case, using the CFD calculation results, the deviation between the simulation results and the actual situation is large. Hence, this kind of underwater vehicle should employ the dynamic model to predict its dive speed. In addition, because the vehicle exhibits the motion of a spiral swing, 2-D dynamic models are not appropriate for this kind of vehicle. A six-degree-of-freedom (6-DOF) dynamic model should be employed to accurately describe this motion. However, as described earlier, the conventional models such as the SUBOFF model, the white box model, and the ROV model are not appropriate for the modeling of the self-oscillation of the foil-like underwater vehicle. In addition, because the swing motion is difficult to explain, it is necessary to treat the damping force as a black-box model. This study uses a grey-box method to build the model in which the damping force is considered as a black-box part. The viscous force model is built using all 378 viscous drag coefficients which are obtained using the Taylor series expansion, and an improved ridge regression method is employed to solve the problem.
The rest of this paper is organized as follows. In Section 2, the modeling and optimization methods of the underwater vehicle are introduced. Section 3 focuses on the implementation of this method and results of the CFD calculation and the simulation results of the 6-DOF dynamic model are presented. The Hadal ARV is used as an example to build a hydrodynamic model by using this method. In Section 4, the results of the simulation are compared with the results of actual sea trials, and the velocity prediction capability is compared with that of the conventional free-fall dive model. The final section summarizes our approach.

Hydrodynamic modeling method
This study focuses on synthesizing the viscous force equation with 378 viscous hydrodynamic coefficients obtained using the Taylor series expansion. An improved ridge regression method is utilized to solve the problem of multicollinearity when the coefficient matrix is huge. The implementation of this modeling method includes the following steps.

Rigid-body dynamics
In this study, the north-east-down (NED) coordinate system is chosen as the geographic coordinate system. The body coordinate system is defined as shown in Fig. 2, where B is the buoyancy of the vehicle and G is the gravity of the vehicle. In Fig. 2, the x direction is the forward direction that indicates the sailing direction; the y direction is the traverse direction, and there is no active motion along this direction; and the z direction is the vertical direction that shows the direction of ascent or descent.
The object of this study is a foil-like underwater vehicle with a large metacentric height. A 6-DOF rigid-body dynamic model with its gravity center away from the buoyancy center is mentioned in another study . As the center of gravity of the vehicle is located at a farther position from the center of the buoyancy, the vehicle has a larger attitude restoration capability. The 6-DOF motion equation can be expressed as follows: where is the inertial matrix, is the Coriolis force matrix and represents the force exerted by the water on the vehicle. The external force is shown as . is the velocity and angular velocity of the vehicle in 6-DOF.
can be divided into three sections.
where is added mass matrix, is the Coriolis force matrix induced by the added mass force, is the viscous force matrix, is the effect induced by gravity and buoyancy, and the mass of the Hadal ARV is m.
, and are more easily to compute .

Hydrodynamics model
The modeling of fluid forces requires the support of experimental data. Towing tank testing and plane motion mechanism (PMM) testing are commonly used simulation test methods. However, during the design phase of the robot, the shape of the vehicle is often changed, and it is uneconomical to repeatedly test the new model using a towing tank or PMM test . In Eq. (2), and can be obtained by the Hess-Smith method (Antes and Panagiotopoulos, 1992).
is difficult to express directly in a linear matrix form. In most cases, is written as a Taylor series expansion, as shown in Eq. (3).
has 63 parameters in each direction, totaling 378 parameters on the six force and torque directions in F x , F y , F z , T x , T y , and T z . To reduce the impact of a manual selection on the number of parameters, all 378 parameters will participate in optimization and modeling.
In the dynamic modeling of underwater vehicles, the number of viscous force parameters is generally small, and each number of force parameters in one direction is often smaller than 20. The collinearity between the input data is unobvious and the fitting and data acquisition processes are also easy. For instance, among the 6-DOF standard submarine equations of motion (Gertler and Hagen, 1967), the pitch moment equation has the maximum numbers of damping force parameters, in a total of 14 items. Maalouf et al. (2015) used only two damping force parameters to characterize the damping force in each degree of freedom. When the number of parameters increased to 63, the choice of an optimization method will become significant. Owing to the large number of parameters to be fitted and the problem of multicollinearity between the input data, it is difficult to obtain valid results by using unbiased estimators such as the conventional least squares method. The direct use of least squares method brings a poor prediction effect and simulation divergence. The ridge regression method with crossvalidation is utilized to improve the prediction accuracy. In addition, the gradient method of adding a correlation analysis is used to solve the problem of weak convergence caused by multicollinearity.
Furthermore, the CFD method is utilized to obtain the modeling data. In this work, the ridge regression method is used, which can draw a ridge trace. The quality of modeling data can be evaluated visually based on the ridge trace . This can guide the designer to select appropriate CFD data to improve the modeling accuracy.

Regression method
For a least-squares linear regression problem, the linear system is represented by Eq. (4) The optimal objective function is arg , The least-squares estimate iŝ u|u| u|u| Because of the strong collinearity between the columns in X, for example, the output of the parameters u, u 2 , and are very close when u<1. u 2 and are approximately linear. In this case, the gradient approximation leads to the singularity of the result of , and furthermore, leads to the estimation result sensitive to the noise of y, resulting in inaccurate estimates . When the ill-posed problem occurs, the mean squared error of the least-squares result is large. The ridge regression method can be used to solve the problem. The regression coefficient, obtained using the ridge regression method, is added to the optimization objective function, which is shown as The estimation result is expressed as: where the ridge regression parameter is k; the adjustment of k can effectively reduce the singular result. The parameter k is determined through the cross-validation method; 80% of the data is considered training data and 20% the validation data. To improve the prediction accuracy of the model during the cross-validation process, validation data should not participate in the training. In addition, a ridge trace is drawn, and k is set to be the minimum value of the crossvalidation result.
u|u| Another problem encountered in the optimization process is the lack of data incentives and the widening of residuals caused by ridge regression. For example, if a vehicle cannot move backward, the parameters u 2 and of this vehicle will be the same and cannot be distinguished. The ridge regression result has higher residuals than the least squares result due to the introduction of regularization terms, which means the predicted value of the regression parameters obtained in the simulation is always smaller than the actual value. In the optimization process, introduction of prior knowledge can better improve the following two problems by referencing pre-obtained correlation coefficients into gradient method. This method can produce the effect of artificial intervention and does not affect the condition that the algorithm can converge to a minimum value. This study calculates the correlation coefficient (Zou et al., 2003) between each group of force and each velocity term before optimization. Let the square of the correlation coefficient of each group data be . This work constructs the following optimization objective function.
arg min When C is not all 0, Eq. (9) is positive definite, and the optimal solution exists. The derivative of Eq. (9) yields β Hydro j = c 2 j β j From Eq. (9), it can be seen that the actual hydrodynamic coefficient value can be substituted into Eq. (8) and Eq. (9), the following formula is obtained.
The comparison between Eq. (8) and Eq. (10) shows that the introduction of the correlation coefficient is essentially equivalent to the re-allocation of various parameters in the regularization terms of the ridge regression. However, avoid the arithmetic error caused by the division and improve the operation precision when c j tends to be infinitely small. In addition, from Eq. (10), it can be seen that the parameter with a larger correlation reduces the contribution of the regularization term to the gradient due to the large value of c j , so that the result is closer to the result of unbiased estimation. For parameters with less correlation, the contribution of the regularization term to the gradient can be greatly improved, and the coefficient can be prevented from being larger.

Modeling implementation
The modeling method described in this paper includes three steps. The first step is to obtain the mass parameter through the white-box method. The second step is to use the CFD method to obtain the modeling data. The third step is to use the CFD data and the improved ridge regression algorithm to build the grey-box model of the viscous force.
Most full-ocean-depth underwater vehicles apply a freefall strategy to touch the seabed (Murashima et al., 2009). In this study, the CFD data are used to build the model, and the simulation results are compared with the free-fall data of the Hadal ARV on the sea trails.

Mass properties
Because of the complexity of the mass property calculation of the vehicle, the experimental data are combined to obtain its mass properties through the equivalent model. The metacentric height of the ARV is measured by the pool experiment, and the simplified model is designed with a similar metacentric height. This vehicle has many open areas that contain a large amount of seawater. The internal water will influence the moment of inertia of the vehicle during the rotation. These parts of the internal water are incorporated into the calculation of mass properties. When the vehicle is working undersea, it is in the state of neutral buoyancy. According to the shape of the outer wrap surface, the volume of the drainage volume can be acquired, and the Solidworks model of the whole vehicle is divided into two parts. Then, the density of the two parts is adjusted with density modifiable materials. Fig. 3 shows the equivalent model with a mass and metacentric height close to those of the actual ARV when the density of the upper part is 160 kg/m 3 , and the density of the lower part is 2000 kg/m 3 . The following mess parameters can be obtained by Solidworks with the equivalent model above-mentioned.
where J is the moment of inertia, m is the equivalent mass of the robot in water, V is the volume of the robot drainage, and is the position of the mass center in the body coordinate system. The mass parameters in Eq. (11) are not exactly the same as the actual parameters of the robot. They are also more stable in the simulation and similar to those obtained from the theoretical calculation results.

CFD calculation
CFD simulations are utilized instead of PMM tests to collect modeling data. Owing to the swing motion, the results obtained by using the steady solver are inaccurate. Therefore, an unsteady solver is utilized for the calculations. As shown in Fig. 4, the size of the model calculation is the same as the actual size of the robot. Meanwhile, all the appendages participated in the CFD calculations. The outer wrap surface for the appendages is created, and then the surface mesh and the volume mesh are drawn.
A cubic region of 4 m×4 m×6 m with a total of 550000 polyhedral cells is chosen. The distance from the top of the robot to the boundary is set at 3.3 m, and the distance from the boundary to the other side of the robot is about 1.5 m. The K-Epsilon model is utilized as the turbulence model, the DFBI function of StarCCM+ is utilized to perform unsteady-state simulations and all the motion states in the actual motion of the ARV are considered by the simulation movement as far as possible (Yang et al., 2015). The input parameter is configured as a random number in the vicinity of the actual operating conditions to allow the simulation to go through as many states as possible. Therefore, the mass properties are set to an identity matrix in the DFBI and an external force as shown in Eq. (12).

Hydrodynamic model identification
A total of 1600 sets of CFD data were obtained in 16 s in the last subsection, using the method in Section 2 for identification and modeling. The results of the model's matching precision and cross-validation accuracy are shown in Table 1.
The goodness of fit of 80% forces is greater than that of the classic ridge regression method, and the value of k is larger than that of the classic ridge regression. This also shows that the parameters obtained, which are not essential, can be suppressed to a lower order of magnitude to achieve the purpose of a manual intervention.
Although the method still needs all 378 dumping coefficients in the simulation, this method has a better parameter reduction effect than the least squares method. Out of the 378 dumping coefficients obtained by this method, 91 parameters have an absolute value smaller than 1, and 22 have an absolute value smaller than 0.01. The results obtained using the standard least-squares method are 35 and 9. It can be concluded that the method described in this paper can re-  LIU Xin-yu et al. China Ocean Eng., 2017, Vol. 31, No. 6, P. 773-780 777 duce the number of the effective parameters.

Field experiment result and analysis
The Hadal ARV performed a sea test in December 2015. In the test, the ARV was found to exhibit a similar elliptical swinging motion mode during the free-fall dive. The previously obtained dynamic model was configured according to the position and quality of load blocks during the sea test, and an open-loop simulation trace with a free fall was obtained as shown in Fig. 5.
A simulation program simulates the motion of the Hadal ARV from the surface after the release of a free-fall dive in the open-loop until a cyclical swing appears. In Fig. 5, the dashed line indicates the sea test data, and the solid line indicates the open-loop simulation data. In the simulation, this study used the quaternion kinematics model (Chen et al., 2013) to simulate the kinematic part of the vehicle. In this experiment, the TCM sensor utilized in this robot can only obtain the attitude information and not the angular velocity information (Ji et al., 2016). The simulations of the output of the roll, pitch and depth are compared with the sea test data of December 2, 2015 as shown in Fig. 5. The model identified by the CFD data shows the same vibration pattern as the sea test in the simulation, and the dive speed is basically the same. Although the kinetic model obtained the attitude angle and depth after two integrals, the simulation can still reproduce the more accurate data of the sea test. The amplitude error of the roll angle is smaller than 10° and the pitch angle is smaller than 2°. The accuracy of the vertical velocity is around 8% when the steady state is achieved. Therefore, the expected effect is achieved.
The frequency domain analysis of the simulation results and the sea trial data of December 2, 2015 are shown in Fig. 6. To fully express the response of the system in each modality, the data of the first 100 s are intercepted from a diving process. The data includes not only a stable state but also the process information about the system state that gradually changed from the initial state to the stable state. It can be seen from Fig. 6 that the roll data of the sea test and simulation data have similar trends in frequency. The main frequency components are 0.23 Hz and 0.22 Hz, and the frequency components near 0.06 Hz are decreased similarly. Z w|w| = −41.83 Fig. 7 shows the comparison between the modeling method described in this paper and the conventional drag coefficient modeling method (Liu et al., 2000) in the simulation. Common free-fall dive models are obtained by calculating the steady-state drag coefficient. The underwater drag coefficient of the no-swinging state of the Hadal ARV is . Substituting this coefficient into simulation, and the velocity trace is represented by the red dashed line in Fig. 7. The red solid line is the simulation result of the dynamic model outlined in this paper. The blue dotted line is the dive speed obtained by the difference of the CTD depth data in the sea test. As can be seen from the figure, the drag coefficient model does not consider the vehicle in the dive process of the swing motion on the vertical resistance changes, and thus, the forecast speed is high. The method described in this paper can well predict the change of the dive speed of the vehicle. The average speed of a steady state is about 0.7 m/s, which is in good agreement with the experimental results. It can simulate the free-fall dive process of the Hadal ARV well. The experimental data in Fig. 7 has several distinct noise data between 30 and 40 seconds due to the measurement error of the depth sensor and is magnified after the differentiation. This problem also   exists in other experiments and does not affect the comparison between simulation and experimental result. Furthermore, the model was utilized to predict the dive speed of the Hadal ARV on July 28, 2016 which reached a depth of 10767 m. In this dive, the mess of body increased considerably compared to the dive on December 2, 2015, but the increase in dive speed was not obvious. Fig. 8 shows the comparison of the data of a free-fall dive in the sea trial and the simulation forecast data. In this experiment, the simulation data of the dynamic model was utilized to predict the dive speed and the dive time before the test. The steadystate average dive velocity was about 0.8 m/s, which is in accordance with the simulation results.
In this test, the negative buoyancy increased by nearly 70%, but the steady-state average speed was smaller than 15%, while the predicted result of the drag coefficient model is that the dive speed should increase by 30%. It is also demonstrated that the dynamic model obtained by this method is reliable in simulating the motion characteristics of foillike underwater vehicles. Furthermore, as the mass of the dive load increases, the amplitude of the wobble in the dive expands, leading to the increase in the lateral flow area. This greatly increases the dive resistance. A model that directly uses a drag coefficient cannot represent this characteristic.

Conclusions
In this study, a modeling method was proposed for a foil-like underwater vehicle. The advantages of the method are as follows. Firstly, for an underwater robot without a reference 6-DOF dynamic model, such as a foil-like underwater robot, this method can build a 6-DOF dynamic model for them. Secondly, the modeling process considered the effect of overfitting on the prediction accuracy, and it had a better processing capability for multicollinearity data. In addition, the obtained dynamic model could simulate complex forms of motion. The simulation results of the open-loop with a free-fall dive were in good agreement with the experimental data.
The proposed method realized the fine modeling of the foil-like underwater vehicle and achieved a good prediction effect in the experiment. It simulated the swing characteristics of the vehicle and predicted the diving resistances in swinging. However, this method lacked the support of PMM experiments or other test data during the modeling phase. In the case of modeling using only CFD data, we found that if researchers lacked relevant experience in system identification, they could set erroneous velocities of movement and angular velocities during the simulation, resulting in inadequate system response and wrong modeling result. In future, the focus will be on obtaining more test data with high accuracy to improve the prediction accuracy.
In this method, the model parameters were large in quantity with a complex model structure. Hence, this model could not directly determine the handling characteristics of the vehicle from the parameters. Although this method has the effect of reducing the parameters involved. However, there is no study on how to select an appropriate reduction threshold. In future, the focus will be on optimizing the threshold selection method to achieve reduction estimates, and to reduce the number of parameters.