Clarifying clairvoyance: Analysis of forecasting models for near-sinusoidal periodic motion as applied to AUVs in shallow bathymetry

Ocean Engineering 190 (2019) 106385

Contents lists available at ScienceDirect

Ocean Engineering journal homepage: www.elsevier.com/locate/oceaneng

Clarifying clairvoyance: Analysis of forecasting models for near-sinusoidal periodic motion as applied to AUVs in shallow bathymetry Jon Woolfrey ∗, Wenjie Lu, Teresa Vidal-Calleja, Dikai Liu Centre for Autonomous Systems, University of Technology Sydney, 81 Broadway Ultimo, 2007, NSW, Australia

ARTICLE Keywords: Wave Motion Prediction AUV Robot

INFO

ABSTRACT This paper shows that Gaussian Process Regression (GPR) with a periodic kernel has better mean prediction accuracy and uncertainty bounds than time series or Fourier series when forecasting motion data of underwater vehicles subject to wave excitation. Many robotic systems, such as autonomous underwater vehicles (AUVs), are required to operate in environments with disturbances and relative motion that make task performance difficult. This motion often exhibits periodic, near-sinusoidal behaviour. By predicting this motion, control strategies can be developed to improve accuracy. Moreover, factoring in uncertainty can aid the robustness of these predictive control methods. Time series and Fourier series have been applied to several predictive control problems in a variety of fields. However, there are contradictory results in performance based on parameters, assessment criteria, and application. This paper seeks to clarify these discrepancies using AUV motion as a case study. GPR is also introduced as a third candidate for prediction based on previous applications to time series forecasting in other fields of science. In addition to assessing mean prediction accuracy, the ability of each model to adequately bound prediction error is also considered as a key performance indicator.

1. Introduction Robotic systems have the ability to perform intervention tasks in environments that humans are otherwise incapable of. This could be due to environmental hazards, or a lack of physical aptitude. Often, this involves operating within the context of natural environs which produce disturbances that affect performance. In certain situations, these disturbances manifest themselves as periodic, near-sinusoidal motion. The movement of marine vessels under wave excitation is exemplary of this kind of problem (e.g. Fig. 1). In such cases, conventional methods such as feedback control may not be able to adequately compensate for this undesired motion imposed on the system. Methods such as Model Predictive Control (MPC), on the other hand, formulate a sequence of control actions across a finite control horizon to minimize an objective function (typically task error). It follows that, if the motion disturbances can be forecast for this horizon, it would be possible to negate them allowing for more precise control. These periodicities and patterns in the disturbance lend themselves to mathematical modelling, which can then be used to forecast said disturbances. Indeed, these types of phenomena have already been modelled and applied to predictive control schemes in a variety of marine robotic systems (Riedel and Healy, 1998; From et al., 2011; Medagoda and Williams, 2012; Kormushev and Caldwell, 2013; Woolfrey et al., 2016; Fernandéz and Hollinger, 2016). Similar strategies

have also been applied to surgical applications. For instance, Becker et al. (2008) attempted to predict the tremor of a surgeon’s hand during teleoperated microsurgery. Heart motion prediction was also explored by Yuen et al. (2008), from which Bowthorpe and Tavakoli (2016) successfully applied a predictive control method for robot-assisted beatingheart surgery. Despite applying similar prediction models and control strategies, there still exists discrepancies in performance within and between fields. 1.1. Motion prediction of marine vessels Wave theory asserts that the undulation of fluids can be modelled via a superposition of sinusoids (Dean and Dalrymple, 1991). This assumption has formed the basis for several prediction modelling methods in literature. For instance, Chung et al. (1990) used a discrete Fourier series to model and predict the motion of a ship in waves for the purpose of improving the accuracy of on-board targeting systems. From et al. (2011) compared this method to time series known as autoregression (AR) for ship motion prediction in sea swell. They concluded the latter to be the better method. Moreover, they noted that AR had lower prediction error variance than Fourier series. In contrast, Yang et al. (2008) stated that AR had poor long-term prediction accuracy

∗ Corresponding author. E-mail address: [email protected] (J. Woolfrey).

https://doi.org/10.1016/j.oceaneng.2019.106385 Received 28 December 2018; Received in revised form 27 August 2019; Accepted 27 August 2019 Available online 17 September 2019 0029-8018/© 2019 Elsevier Ltd. All rights reserved.

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

1.3. Motion prediction in related fields Periodic motion compensation has also been a topic of research in other fields. For instance, Yuen et al. (2008) compared different models for heart motion prediction with the intended use for teleoperated, beating-heart surgery. They also tested Fourier series due to the periodic motion of the heart. However, unlike Chung et al. (1990) and From et al. (2011), they applied an Extended Kalman Filter (EKF) which can track the changing periodicity. This method was compared to an AR(30) model. Their conclusion was that the Fourier series with EKF outperformed time series. This research was then successfully applied to a predictive control method by Bowthorpe and Tavakoli (2016). By using the motion prediction, it was shown that tracking error between a surgical tool and the surface of the heart could be reduced. Time series forecasting has also been applied to modelling and predicting unwanted hand tremor during microsurgery (Becker et al., 2008). The authors implemented both AR, and autoregressive moving averages (ARMA). Their conclusion was that an ARMA model performed better at predicting this form of periodic motion than a simpler AR format.

Fig. 1. The Submerged Pile Inspection Robot (SPIR) may be required to perform under waves and turbulence in and around bridge structures.

for ship motion prediction. Instead, they developed a method that used past observations of ship motion plus wave excitation input to produce predictions for maritime flight operations. Love et al. (2004) showed that the inertial forces produced by sea swell on ship-mounted manipulators could be significantly large such that feedback control was insufficient for compensation. Instead, they used Repetitive Learning Control (RLC) to estimate and counteract these dynamic disturbances. However, they noted that this method failed when periodicity was not constant. This is an important factor to consider if operating in the natural environment. In the same line of research, From et al. (2009) showed that the dynamic properties of the ship have no effect on the manipulator dynamics. All that is required is knowledge of the ship motion. In that regard, they showed that knowledge of the future ship motion could be used in motion planning to reduce joint torque of ship-mounted manipulators. This work was extended by From et al. (2011) in which actual predictions were used. Furthermore, factoring in joint state uncertainty improved the robustness of the predictive motion planning.

1.4. Contribution and outline Evidently, the ability to forecast near-sinusoidal, periodic motion has wide applications for predictive control strategies in robotics. Although they have been applied successfully, there are still several questions that remain which this paper will attempt to address: • Did the AR model perform better than Fourier series for ship motion prediction by From et al. (2011) because it is better suited to the application? Or would the implementation of the Fourier series with EKF used by Yuen et al. (2008) improve results? • And would an ARMA model, like that used by Becker et al. (2008), be better than AR for motion prediction of marine vessels as used by From et al. (2011) and Woolfrey et al. (2016)? • Furthermore, how well do each of the prediction models adequately capture uncertainty about the future state? Forecast uncertainty is often overlooked in literature on predictive control, despite its utility in planning and optimization. As Silver (2012) writes:

1.2. Disturbance prediction for AUVs Some of the earliest work for wave compensation in autonomous underwater vehicles (AUVs) was conducted by Riedel and Healy (1998). A Kalman Filter (KF) was used to estimate the external sea effects. This was then used with Sliding Mode Control (SMC) to improve the accuracy of station keeping under oscillating forces. In later work, Medagoda and Williams (2012) used an Acoustic Doppler Current Profiler (ADCP) to sense the current flow beneath an AUV descending in a water column. MPC was then used to proactively counter the anticipated disturbances and reduce trajectory tracking error. Riedel and Healy (1998) had actually argued against the use of remote sensing technology in order to reduce weight and monetary costs in the AUV design. This necessitates that some form of on-board estimation be developed. In more recent work, Fernandéz and Hollinger (2016) applied MPC for stationkeeping to a simulated, 2D AUV. They showed reduced position error over a PID feedback controller. However, simulation data was given to the MPC controller; no forecasting method was presented. The ability to predict disturbances can also be utilized to minimize energy consumption. For instance, Kormushev and Caldwell (2013) modelled disturbances to an AUV using a superposition of sinusoids. Using machine learning, the amplitudes and frequencies were adapted to best approximate the local conditions. This was then used in a control strategy that enabled the system to reduce energy expenditure. Forecasts of sea currents were utilized in research by Huynh et al. (2015) and Jones and Hollinger (2017) to plan efficient trajectories over long distances. Factoring prediction uncertainty was shown to improve robustness.

The words predict and forecast are largely used interchangeably today, but in Shakespeare’s time, they meant different things ... Making a forecast typically implied planning under conditions of uncertainty. It suggested having prudence, wisdom, and industriousness, more like the way we now use the word foresight. From et al. (2011) had opted for the AR model over Fourier series in part because the former had lower prediction error variance than the latter. This uncertainty was used as part of their optimization problem. However, it is not evident that low uncertainty in a forecast is desirable per se. On the contrary, overconfidence may lead to sub-optimal or unsafe actions. An interesting piece of research by Fang and Williams (2014) showed that statistical distributions beyond the conventional Gaussian distribution can provide better uncertainty bounds around trajectories. This paper provides a systematic comparison of the forecasting abilities between time series, Fourier series, and Gaussian Process Regression (GPR) using AUV motion data as a case study. This includes both the mean accuracy and uncertainty. GPR is considered as a third candidate due to its application in time series prediction for wind power generation (Yan et al., 2016), physiology (Dürichen et al., 2015), and many other natural phenomena (Roberts et al., 2013). Section 2 covers time series modelling, with some statistical analysis of motion data to guide future model development. Section 3 outlines the Fourier series and EKF method. Section 4 covers GPR and the model 2

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

selection for the kernel functions. In Section 5, different metrics for assessing forecast quality are examined. In particular, Theil’s Inequality Coefficient (Theil, 1966) is opted for as a dimensionless metric to better compare linear and angular motion. The standard error of the forecasts is also considered as a means of assessing forecast uncertainty. Section 6 covers the analysis of forecasting motion data for two different AUVs under wave excitation. The three aforementioned methods are assessed, with 2 variations on the parameters for each. 2. Time series forecasting 2.1. Structure of a time series model In time series modelling, the expected value of a state observation y(t) ̂ is expressed as a weighted sum of previous observations and forecast errors: y(t) ̂ =

p ∑ i=1

𝛼i y(t − i) −

q ∑

𝛽j 𝜀(t − j)

Fig. 2. Correlograms for a sample of AUV motion data, showing statistically significant autocorrelation back in time.

(1)

j=1

where: 2.3. Producing forecasts with time series

• p is the number of autoregression (AR) terms, • q is the number of moving average (MA) terms, and • 𝜀 = ŷ − y ∼  (0, 𝜎𝜀2 ) is the prediction error.

After optimizing the coefficients in Eq. (1), From et al. (2011) and Woolfrey et al. (2016) used an iterative process to calculate predictions using an AR model:

Eq. (1) is referred to as an ARMA(p,q) model, which assumes stationarity in the data, i.e. it exhibits a long-term mean. This equation can be modified to account for nonstationary data. A thorough reference on time series can be found in Box et al. (2016). The coefficients 𝛼i for i ∈ {1, … , p} and 𝛽j for j ∈ {1, … , q} can optimized via linear least squares. These can be fitted at each time step in a control loop using a moving window of 𝑙 historical state observations. This approach was taken by From et al. (2011) and Woolfrey et al. (2016) so that the values best reflected the current conditions.

y(t ̂ + i|t) =

p ∑

𝛼j y(t ̂ + i − j|t),

for i ∈ {1, … , L}.

(2)

j=1

where L is the number of lead terms or forecast steps. However, through recursive substitution of y(t − i) in Eq. (1) with its own ARMA model, it is possible to arrive at an equivalent state space form Box et al. (2016): ̂ + 1|t) = 𝝓(t)𝐲(t|t) ̂ 𝐲(t + 𝝍(t)𝜀(t),

(3)

where:

2.2. Determination of model parameters

y(t|t) ̂ ⎡ ⎤ ⎢ y(t ̂ + 1|t) ⎥ ⎢ ⎥ ̂ • 𝐲(t|t) =⎢ ⋮ ⎥, ⎢y(t ̂ + L − 1|t)⎥ ⎢ ⎥ y(t ̂ + L|t) ⎦ [⎣ ] 𝟎 𝐈 • 𝝓(t) = L×1 , 0 𝛼L−2 ⋯ 𝛼2 𝛼1 [ ]T • 𝝍(t) = 𝜓0 𝜓1 ⋯ 𝜓L−1 𝜓L ,

Several methods have been applied in literature to determine the number of parameters for p and q in Eq. (1). A popular method is to evaluate the root mean squared error (RMSE) of the predictions for various models (Becker et al., 2008; From et al., 2011; Woolfrey et al., 2016). However it is well known that linear regression methods can overfit data. To account for this, the Akaike Information Criterion (AIC) can be used which includes a penalty term for the number of model parameters, thereby promoting model parsimony. This was the approach taken by Yuen et al. (2008). Conversely, the AIC was criticized by Schwarz (1978) as been unreliable for large sample sizes, and instead proposed the Bayesian Information Criterion (BIC). This was used by Yang et al. (2008) for fitting parameters in their time series model for ship motion prediction. However, literature on time series modelling asserts that the model parameters are determined by the statistical properties of the data (Box et al., 2016). Namely, the number of p terms is determined by the number of statistically significant partial-autocorrelation values in the data, and the number of q terms is given by the number of statistically significant autocorrelation values. With this regard, the correlograms for a sample of IMU data for an AUV under wave excitation is shown in Fig. 2. There a numerous autocorrelation values in the data. But, an infinite MA process is equivalent to a finite AR process (Box et al., 2016). Moreover, the partial-autocorrelation values show a sharp decline, albeit with several significant values. These two observations suggest a high-order AR process. An ARMA(40,0) model is applied in this paper as it was found to give satisfactory prediction results with some balance for model parsimony. And, although the data do not indicate it, an ARMA(40,1) model may also help improve forecasts by correcting for the past prediction errors.

• 𝜓(𝛼, 𝛽) is a function of the coefficients in Eq. (1). The advantage to this method is that the state propagation of Eq. (3) gives all the predictions across the forecast horizon. Moreover, the application of a KF will correct both the current state estimate and the predictions when new sensor information is acquired. The uncertainty propagation of this model is then readily calculated using: 𝜮 y (t + 1|t) = 𝝓(t)𝜮 y (t|t)𝝓T + 𝜎𝜀2 𝝍(t)𝝍(t)T ,

(4)

𝜎𝜀2

where is the variance of the 1-step-ahead prediction errors. Eq. (4) gives both the current state uncertainty, as well as the forecast uncertainty. The observation matrix for the KF is given as: [ ] 𝐇 = 1 𝟎1×L , (5) since future states cannot be observed. 3. Fourier series 3.1. Superposition of sinusoids Wave theory asserts that the undulation of a body of water can be expressed as a superposition of sinusoids (Dean and Dalrymple, 3

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

1991). Indeed this assumption has been applied in many modelling and predictive control problems (Chung et al., 1990; From et al., 2011; Kormushev and Caldwell, 2013; Fernandéz and Hollinger, 2016). Using this method, the state equation for a system under periodic motion can be expressed as a discrete sum of m sinusoids: y(t) ̂ = c0 (t) + = c0 (t) +

m ∑ i=1 m ∑

ri (t) sin(𝜃i (t))

(6)

ai sin(i𝜔t) + bi cos(i𝜔t)

(7)

for i = {1, … , m}. For the initial state uncertainty, Yuen et al. (2008) and Bowthorpe and Tavakoli (2016) had opted for a diagonal matrix. In this paper, the uncertainty is determined by mapping Eq. (10) via the Jacobian: ( ) ( )T ̂ 𝜷̂ 𝜮 𝛽 𝜕 𝐱∕𝜕 ̂ 𝜷̂ + 𝜮 w , 𝜮 x (0) = 𝜕 𝐱∕𝜕 (12) where 𝜮 w is additional uncertainty to account for loss of information through linearization. The state propagation is given by:

i=1

𝜃i (t) = i

∫0

(8)

𝜔(𝜏)d𝜏 + 𝜙i (t).

where 𝐪 ∼  (𝟎, 𝜮 q ), and the state propagation matrix is defined as: ⎡𝐈m+2 ⎢ ⎢ 𝐀≜⎢ 𝟎 ⎢ ⎢ ⎣

Eq. (6) is the polar form, and Eq. (7) is the rectangular form arrived at through trigonometric identities. The coefficients in Eq. (7) can be initialized using a Best Linear Unbiased (BLU) estimate: ( )−1 T −1 𝜷̂ = 𝐗T 𝜮 −1 𝐗 𝜮y 𝐲 (9) y 𝐗 ( T −1 )−1 𝜮𝛽 = 𝐗 𝜮y 𝐗 , (10) [ • 𝜷̂ = ĉ 0 â 1 coefficients, [ • 𝐲 = y(t) ⋯

⋯

â m

̂b1

⋯

̂bm

]T

is the estimate of the

]T y(t − 𝑙) are prior sensor or state observations,

• 𝜮 y is the uncertainty of said observations, and ⎡ ⎤ 1 ⋯ 1 ⎢ sin(𝜔𝛥t) ⋯ sin(𝜔𝑙𝛥t) ⎥ ⎢ ⎥ ⋮ ⋯ ⋮ ⎢ ⎥ • 𝐗T = ⎢ sin(m𝜔𝛥t) ⋯ sin(m𝜔𝑙𝛥t) ⎥. ⎢ ⎥ ⋯ cos(𝜔𝑙𝛥t) ⎥ ⎢ cos(𝜔𝛥t) ⎢ ⎥ ⋮ ⋯ ⋮ ⎢ ⎥ ⎣cos(m𝜔𝛥t) ⋯ cos(m𝜔𝑙𝛥t)⎦

â j sin(𝜔i𝛥t) + ̂bj cos(𝜔i𝛥t),

1 0 ⋮ 0

𝟎 0 1 ⋮ 0

⋯ ⋯ ⋱ ⋯

0 0 ⋮ 1

⎤ ⎥ ⎥ ⎥. ⎥ ⎥ ⎦

(14)

The forecasts using the EKF method differ from Eq. (11) in that the intrinsic parameters and uncertainty must first be propagated: ̂ + i|t) = 𝐀i 𝐱(t|t) ̂ 𝐱(t 𝜮 x (t + i|t) = 𝐀i 𝜮 x (t|t)(𝐀i )T +

(17) i−1 ∑

𝐀j 𝜮 q (𝐀j )T ,

(18)

j=0

for i ∈ {1, … , L}. Then the predictions of the actual state outcome and uncertainty are calculated as:

Much like the time series application, authors like Chung et al. (1990) and From et al. (2011) optimized the coefficients for Eq. (7) using an unweighted form of Eq. (9). Then, forecasts were made recursively using: m ∑

𝛥t 2𝛥t ⋮ mΔt

For the EKF, the observation matrix is given by the partial-derivative of Eq. (6) with respect to the state vector: ( ) ̂ + i|t) 𝜕y 𝐱(t (15) 𝐇(t + i) ≜ ̂ + i|t) 𝜕 𝐱(t [1 sin(𝜃̂1 ) ⋯ sin(𝜃̂m ) = (16) 0 r̂1 cos(𝜃̂1 ) ⋯ r̂m cos(𝜃̂m )].

where:

y(t ̂ + i|t) = ĉ +

(13)

̂ + 1|t) = 𝐀𝐱(t|t) ̂ 𝐱(t + 𝐪,

t

(19)

̂ + i|t)) y(t ̂ + i|t) = f (𝐱(t 𝜎y2 (t

T

+ i|t) = 𝐇(t + i)𝜮 x (t + i|t)𝐇(t + i) .

(20)

4. Gaussian process regression

(11)

j=1

Both the time series and Fourier series methods use deterministic mathematical models to explain covariance between sets of data. This assumes that observations are generated according to a given mathematical model. Deviations (errors) are regarded as corruptions from noise. Gaussian Processes (GPs), on the other hand, model a distribution over functions. Inferences about an expected output are generated from the covariance between nearby observations. As can be seen from the time-series analysis in Section 2.2, the motion of a vessel under wave excitation exhibits statistically significant autocorrelation in time. This property has been exploited in using GPs to model time-series data in a variety of natural phenomena. Respiratory data was predicted by Brahim-Belhouari and Bermak (2004) for respiratory data, wind-power forecasting by Yan et al. (2016), and multivariate physiological phenomena by Dürichen et al. (2015). Roberts et al. (2013) also offer a variety examples for periodic motion such as tide heights to stellar light curves. The set of current state estimates y(t|t) ̂ and future (predicted) values y(t ̂ + 1|t) can be modelled via a joint Gaussian distribution: ]) [ ] ([ ] [ y(t|t) ̂ 𝜇(t) 𝜮 aa 𝜮 ab ∼ , , (21) y(t ̂ + 1|t) 𝜇(t + 1) 𝜮 Tab 𝜮 bb

for i ∈ {1, … , L}. One flaw with this method is that the fundamental frequency, 𝜔, is assumed to be known a priori. It is possible to use historical information to initialize this term, or even guess the parameter. Moreover, Eq. (9) only optimizes for the amplitudes and phase shifts and not the frequency 𝜔. Chung et al. (1990) had recognized this, and therefore devised a method to update the frequency based on an error gradient method. 3.2. Extended Kalman Filter (EKF) Parker and Anderson (1990) developed an EKF method to track hidden frequencies in noise based on the polar form of the Fourier series, Eq. (6). Their results showed that it was able to sufficiently track the underlying frequency, 𝜔, in a given data set. Yuen et al. (2008) extended this method to produce forecasts of heart motion. Then, Bowthorpe and Tavakoli (2016) successfully integrated these forecasts in to a predictive control strategy. This same method is applied in this paper to predict AUV motion under wave excitation. The parameters for the state vector can be initialized from Eq. (9): [ ]T ̂ • 𝐱(0) = ĉ (0) r̂ (0) 𝜔(0) ̂ 𝜃̂i (0) , ( 2 0 2 )1∕2 i • r̂i = â + ̂b , i

where 𝜇(⋅) is a mean function and the covariance matrices between state observations are: • 𝜮 aa = 𝑐𝑜𝑣(𝐚, 𝐚), • 𝜮 ab = 𝑐𝑜𝑣(𝐚, 𝐛),

i

• 𝜃̂i = atan2(̂bi , â i ), 4

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

• 𝜮 bb = 𝑐𝑜𝑣(𝐛, 𝐛),

forecasts of the GP and often optimization of the hyperparameters did not converge. As such, no further considerations are made.

and input spaces to the time-series GP are: [ ]T • 𝐚 = y(t ̂ − 1|t − 1) ⋯ y(t ̂ − 𝑙|t − 𝑙) [ ]T • 𝐛 = y(t|t) ̂ ⋯ y(t ̂ − 𝑙 + 1|t − 𝑙 + 1)

5. Forecast accuracy 5.1. Root Mean Squared Error (RMSE)

for given inputs 𝐚, 𝐛. A prediction about future expected states given a window of historical state observations can then be made using: ( ) y(t ̂ + 1|t) = 𝜇(t + 1) + 𝜮 Tab 𝜮 −1 ̂ − 𝜇(t) (22) aa y(t|t) 𝜎y2 (t + 1|t) = 𝜮 bb − 𝜮 Tab 𝜮 −1 aa 𝜮 ab .

A conventional approach to evaluate prediction accuracy is to measure the RMSE between the predicted state and the realized state observation: √ )2 ∑( y(t|t ̂ − 1) − y(t|t) ̂ . (28) RMSE = n−1

(23)

A KF can readily be applied to the state estimates from Eqs. (22) and (23).

However, there are two disadvantages to this metric: 4.1. Kernel functions 1. It does not give a measure of prediction performance against the minimum standard of a last-value predictor. Indeed, it is possible to perform worse than this. 2. It does not allow for comparison of a model across different dimensions, in this case linear acceleration, and angular velocity.

The crux of the GP is in the kernel functions of the covariance matrix. A variety of these can be found in work published by Rasmussen and Williams (2006) and Roberts et al. (2013). One useful property of the kernel functions is that additions and products can be constructed to model complex data.

5.2. Theil’s inequality coefficient 4.1.1. Matérn The Matérn function equates to a continuous-time AR process (Rasmussen and Williams, 2006; Roberts et al., 2013) and hence is of particular interest here. As can be seen from the correlogram in Fig. 2, AUV motion under wave excitation has many statistically significant partial-autocorrelation values. The Matérn kernel function is: ( √ ) ( √ ) 1 2 𝜈 ⋅ 𝛿 B𝜈 2 𝜈 ⋅ 𝛿 , (24) 𝑐𝑜𝑣M (a, b) = h2 𝜈−1 𝛤 (𝜈)2

In this paper, Theil’s Inequality Coefficient (Theil, 1966) is opted for over RMSE. This metric is given by the following equation: )2 ∑( y(t|t) ̂ − y(t|t ̂ − i) U2 (i) = ( (29) )2 , ∑ y(t|t) ̂ − y(t ̂ − i|t − i) which can be regarded as a ratio of the predicted error over the actual change. A value of U = 0 would indicate a perfect prediction, whereas U > 1 means the model performs worse than a na´’ive last-value predictor. And, unlike RMSE, this metric is dimensionless. Thus it is possible to compare a model’s ability to predict the linear motion of a vehicle against angular motion.

where 𝛿 = 𝜆−1 |a − b|, 𝛤 (⋅) is the Gamma function, and B(⋅) is the 2nd-order, modified Bessel function. 4.1.2. Rational quadratic This kernel function was combined with a periodic kernel by Roberts et al. (2013) to model changes in periodicity over time for stellar light curves. A similar formulation is thus employed in this paper. The kernel function is given as: )−𝛼 ( (a − b)2 𝑐𝑜𝑣RQ (a, b) = h2 1 + , (25) 𝛼𝜆2

5.3. Standard error of predictions Forecast uncertainty can provide useful information in motion planning and control. For instance, probabilistic bounds have been considered in trajectory generation for robotic systems (Fang and Williams, 2014). From et al. (2011) also attempted to minimize joint state uncertainty for predictive motion planning of manipulators on ship platforms. Evidently, sound information about forecast uncertainty can provide a better basis for planning and control. For example, if the forecast uncertainty underestimates the prediction error, this could lead to a false sense of security in applications that require collision avoidance. Conversely, an overly conservative estimate may hinder planning, exploration, and control. A useful metric is to take the average standardized prediction error, sometimes referred to as the Normalized Estimation Error Squared (NEES) (Li et al., 2001). ( )2 ∑ y(t ̂ + i|t + i) − y(t ̂ + i|t) 2 𝜎SE (i) = n−1 , (30) 𝜎 2 (t + i|t)

where h accounts for amplitude, 𝜆 is a length or time scale, and 𝛼 is known as the index parameter. 4.1.3. Periodic rational quadratic Periodicity can be introduced in to certain kernel functions by ( ) replacing (a − b)2 with sin2 𝜋(a − b)∕T with periodic parameter T. Substituting this in to the Rational Quadratic kernel (36) becomes: ( ))−𝛼 ( 1 |a − b| 𝑐𝑜𝑣PRQ (a, b) = h2 1 + (26) sin2 𝜋 | | 2 | T | 𝛼𝜆 4.2. Mean functions Mean functions influence the behaviour of the expected output in a GP model. The mean function may be useful in accounting for steady disturbances in addition to the periodicity, such as steady current flow. The zero mean, or even a constant mean function is the simplest: 𝜇(t) = c,

∀t.

where 𝜎 2 (t + i|t) is the forecast uncertainty of the 𝑖th step-ahead forecast. The square root of Eq. (30) will give the number of standard deviations that the realized state outcome y(t ̂ + i|t + i) is from the forecast mean y(t ̂ + i|t). Thus, when 𝜎SE > 1.96, then the state observation exceeds the 95% confidence interval of the forecast. This implies the prediction model is unable to sufficiently encompass the prediction errors.

(27)

Assuming zero mean or constant is reasonable for local predictions as the kernel functions can account for complex behaviour. However, distant predictions are dominated by the mean. It was found in this paper that other mean functions, such as linear, did not improve 5

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

6. Comparison of forecast models 6.1. Data collection process Data was collected from two different AUVs designed for intervention tasks in shallow bathymetry. The first is an AUV prototype being developed for exploration under Antarctic sea ice, dubbed the Antarctic-AUV (A-AUV), Fig. 3. The second is an AUV-manipulator system for high-pressure water cleaning of bridge piles and submerged infrastructure, dubbed the Submerged Pile Inspection Robot (SPIR), Fig. 6. Raw sensor data was collected from an inertial measurement unit (IMU) onboard each AUV whilst they were subject to wave motion. The linear motion is in m/s2 and the angular motion is in rad/s. Each of the AUVs was placed in a test tank at separate time, and left to float freely on the surface of the water. Waves were manually generated by using a large object to displace the water. The tank measures 6 m in length and 4 m in width, with the wavelength of 4 m matching the shorter side of the tank due to harmonic resonance. The data were recorded at 50 Hz.

Fig. 3. The Antarctic AUV (A-AUV) prototype for inspection under sea ice.

6.2. Selection of model parameters 6.2.1. Time series Model parameters for the ARMA(p,q) method were chosen based on analysing correlograms for samples of motion data, Fig. 2. Using values of p < 30 produced capricious results in the forecast when coupled with a KF. It was found that an ARMA(40,0) model produced stable forecasts up to 5 s ahead. A lag length of 𝑙 = 250 (5 s) of historical data was used to solve the coefficients in Eq. (1) each time step. This was then used to formulate the state propagation matrices in Eq. (3). Various combinations of an ARMA(40,q) model were also tested, though prediction results did not improve. An ARMA(40,1) is presented in this paper.

Fig. 4. A 30 s sample of the IMU data collected from the A-AUV.

By inspection, it can be seen that the GPR with either a Matérn or periodic kernel has better short-term mean predictions than either time series or the Fourier series methods. The long-term performance results of the ARMA(40,0) time series model has comparable performance to the GPR methods. Both of these also better captured the periodicity in the motion (1.2 s) than the Fourier series. Becker et al. (2008) had claimed that an ARMA model performed better than a simple AR for predicting the unwanted tremor of a surgeon’s hand during microsurgery. The results of this analysis showed that extending the ARMA(40,0) to an ARMA(40,1) did not improve results. In fact, it can be seen in Table 1 that the inclusion of the MA term over-corrects the 1-step-ahead predictions, leading to poorer results. This is to be expected from the analysis of the correlograms in Section 2, Fig. 2, which indicated a high-order AR(p,0) process. The Fourier series with m = 10 harmonics actually proved to have slightly better 1-step-ahead predictions than time series. This was also noted by From et al. (2011) in their analysis of ship motion prediction. However, the long term forecasts are not as accurate, and the method did not adequately track the periodicity in the data.

6.2.2. Fourier series The Fourier series with EKF was implemented with m = 5 and m = 10 harmonics. In general, it was found that a higher number of harmonics lead to the method overfitting the model to noise. To optimize the initial coefficients, Eq. (9) was solved using a lag length of 𝑙 = 400 data points (8 s) of historical data. 6.2.3. Gaussian process The hyperparameters for the kernel functions, Eqs. (24) and (26), can be fitted via gradient descent optimization via the log-likelihood function, Rasmussen and Williams (2006) and Roberts et al. (2013). These were fitted with 600 samples within the data collected. Initial values were randomized and several attempts at optimization were conducted to ensure convergence in the result. A moving window of 𝑙 = 200 (2 s) of historical observations were used for the input space for calculating the covariance matrices used in Eqs. (22) and (23). This was found to produce more accurate predictions than shorter windows.

6.3.2. Forecast uncertainty The standard error of the forecasts, Eq. (30), for the A-AUV motion data are plotted in Table 2 for 10-step intervals. The GPR with the periodic kernel function underestimates the prediction error for forecast horizons <1 s in length, with observations outside the 3𝜎 or 99.97% confidence interval. These generally improve for long-term predictions with forecasts under this 3𝜎 range. The Fourier series, overall, appears to have good short-term uncertainty bounds for m = 10 harmonics. Although the long-term uncertainty intervals encompass the realized state observation, this usually falls within the 1𝜎 range (68%) and hence could be regarded as too conservative

6.3. Antarctic AUV The A-AUV weighs approximately 17 kg. The three prediction methods were tested on 60 s worth of IMU data. Thirty (30) seconds worth of this data is shown in Fig. 4. Forecasts were made in 0.02 s increments, for a maximum of L = 250 steps equating to a 5 s forecast lead. 6.3.1. Mean forecast accuracy Results for Theil’s Inequality Coefficient, Eq. (29), were averaged for 10-step intervals up to 5 s in total length, which can be seen in Table 1. 6

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al. Table 1 Theil’s coefficient averaged for different forecast horizons, by model and degree-of-freedom for the A-AUV motion data.

Table 2 Standard error of the forecasts for different horizons, by model and degree-of-freedom for the A-AUV motion data.

The time series method overall had very poor prediction uncertainty in excess of 3𝜎, and thus state observations routinely fall outside the

99.97% confidence interval. From et al. (2011) had noted that the time series had low forecast error variance, and opted for this over Fourier 7

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

Fig. 5. Samples of the forecasts of the A-AUV motion data, with mean in red and 95% confidence interval in grey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

series for their MPC problem. However, the findings here would suggest the forecast uncertainty of the time series is unreliable for the data used in this paper.

As before, forecasts were made at 0.2 s intervals for L = 250 steps, equating to a maximum 5 s lead time. Theil’s Inequality Coefficient, Eq. (29), and the standard errors, Eq. (30), were calculated and averaged for every 10th step.

6.4. Submerged Pile Inspection Robot (SPIR)

6.4.1. Mean forecast accuracy Theil’s Inequality Coefficient for the SPIR motion prediction is shown in Table 4. It can be seen that this data proved much more difficult to predict for all the forecasting methods presented in this paper. In particular, Theil’s Inequality Coefficient U > 1 for the surge

The SPIR weighs approximately 45 kg, Fig. 6. As with the A-AUV, the different prediction methods were tested on 60 s worth of data. A 30 s sample is shown in Fig. 7. 8

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

Table 3 Summary of performance aptitudes for the forecasting methods. Time series Short term

Mean Uncertainty

Long term

Mean Uncertainty

Fourier series

GP

Periodicity

7. Discussion 7.1. Time series vs Fourier series From et al. (2011) stated that an AR time series model was better for predicting ship motion than Fourier series. Conversely, Yuen et al. (2008) concluded that Fourier series using an EKF was better at predicting heart motion than AR. The question was posed in the Introduction as to whether the application of the EKF could improve forecasts for marine vessels under wave excitation over time series. The results of this paper showed that an AR(40) model with state space form produced more accurate long-term predictions than the EKF method with the data presented. These results are to be expected, as the EKF method was originally proposed by Parker and Anderson (1990) to track hidden frequencies in periodic data, and not for producing forecasts per se.

Fig. 6. The Submerged Pile Inspection Robot (SPIR) designed for high-pressure cleaning of bridge support columns.

7.2. Model parameters Becker et al. (2008) had found that using an ARMA model was better at predicting the periodic tremors of a surgeon’s hand during teleoperated surgery, compared to a simple AR. Analysis of the statistical properties for some of the motion data indicated a high-order AR model. Moreover, the results of the mean prediction accuracy showed that an ARMA(40,1) model did not significantly improve motion prediction for AUVs over an AR(40,0) model. In fact, the additional MA term appeared to over-correct the 1-step-ahead predictions, Table 1. For the Fourier series, it was found that using m = 10 harmonics produced the best overall results for all DOF of vehicle motion with stationary data. These results are commensurate with the m = 8 harmonics that Yuen et al. (2008) used for heart motion prediction.

Fig. 7. A 30 s sample of the IMU data collected from the SPIR.

motion with all the prediction methods tested. This implies that, in practice, it would be better to use a last-value predictor. The GP methods had very accurate 1-step ahead predictions (U < 0.01), and the periodic kernel had generally good predictions for the dimensions other than surge for forecasts <2 s in length. The Matérn kernel performed poorly for long terms predictions. On further analysis of the SPIR motion data, it was found to be nonstationary. This may explain the poor results of the forecasting methods compared to the A-AUV data. Ironically, the Fourier series with m = 5 was better at forecasting surge motion than all the other models.

7.3. Forecast uncertainty The forecast uncertainty was considered a critical factor in this paper for evaluating a prediction model. From et al. (2011) had noted that the long-term prediction errors had lower variance than that of Fourier series. The results in this paper showed that the forecast uncertainty for AUV motion prediction is very small, leading to large standard errors for the predictions (Table 2, 5). The reason for this is clear. The uncertainty propagation for the time series model, Eq. (4), is a function of the 1-step-ahead prediction error variance 𝜎𝜀2 . Thus, a small variance in prediction error leads to overconfidence in the forecast uncertainty. Likewise, the GP method also greatly underestimated prediction error when applied to the nonstationary motion data of the SPIR. If these methods are to be applied in predictive control, the forecast uncertainty would need to be manually increased to account for this overconfidence.

6.4.2. Forecast uncertainty The results of the forecast uncertainty, Table 5, also prove to be much worse overall for the SPIR motion data than the A-AUV. Again, the time series is overconfident in its predictions with standard errors in excess of 3𝜎. And, both the Matérn and periodic kernel functions for the GP have significantly large standard errors for the predictions.

6.5. Performance summary Based on the analysis of the different forecasting methods, a summary of their different performance strengths is given in Table 3. Fourier series alone had the best short-term forecast uncertainties, although failed to capture periodicity. Time series appears to have good long-term mean predictions and estimates of the periodicity in the data. Overall, the GP with periodic kernel performed best with the stationary data.

7.4. Vehicle geometry and data stationarity Both the A-AUV and the SPIR motion data were collected in the same test tank. However, it was found that the A-AUV data exhibited 9

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al. Table 4 Theil’s coefficient averaged for different forecast horizons, by model and degree-of-freedom for the SPIR motion data.

Table 5 Standard error of the forecasts for different horizons, by model and degree-of-freedom for the SPIR motion data.

10

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

Fig. 8. Samples of the forecasts of the SPIR motion data, with mean in red and 95% confidence interval in grey. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

stationarity whereas the SPIR did not. As a consequence, all the prediction models presented performed poorly when forecasting the SPIR motion data. According to Dean and Dalrymple (1991), an unpropelled object floating in water is dominated by the particle accelerations and the pressure field. As such, the different geometry between the two vehicles may help explain the discrepancies in the observed motion. Both vehicles are neutrally buoyant and sit fully submerged in the water column. The length of the waves during the data collection was approximately 4 m due to resonance with the 4 m width of the tank. The wave height reached approximately 320 mm. Fig. 9 illustrates the

pressure gradients from the surrounding water on the A-AUV and SPIR from a passing wave. Because the body of the SPIR is longer than that of the A-AUV, a minute net force would be generated by a passing wave. This would cause it to drift over time with the propagation of the waves, which was also observed during data collection. The force in the 𝑥-direction from this pressure field on a rectangular object is given as Dean and Dalrymple (1991): ( ) −2𝜌gH sinh(kh) − sinh k(h − d) Fx = ⋅ ⋅𝛿 (31) k sin(𝜃) k cosh(kh) 11

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

Fig. 9. The body of the SPIR (top) is longer than that of the A-AUV (bottom), and hence is affected differently by passing waves.

Fig. 11. Underlying curvature in the SPIR motion data leading to nonstationarity.

7.5. Potential avenues to improve forecasting There are a couple of methods that can be used to deal with the nonstationary nature of the SPIR motion data. The first is to use an autoregressive integrated moving average (ARIMA) model (Box et al., 2016). Using this method, the data is differenced until stationarity is achieved. An ARIMA(p,1,q) model, for example, would predict 𝛥ŷ instead of y. ̂ Another option is to use GPR with a nonstationary kernel function. Brahim-Belhouari and Bermak (2004) showed good results for predicting respiratory data with a nonstationary kernel. This may be a point of future research. That said, the time series and GP with periodic kernel performed better than a last-value prediction for prediction horizons <2.5 s for all but the surge motion. Fig. 10. The forces acting on the SPIR in a wave field will differ from that of the A-AUV due to their geometry.

8. Conclusion There are many scenarios in which robotic systems are subject to disturbances which leads to periodic, near-sinusoidal behaviour. Both time series and Fourier series have been successfully applied to predictive control strategies in previous literature. However, performance results in the prediction accuracy of each method differs depending on derivation, parameters, and application. In this paper, time series, Fourier series, and Gaussian Process Regression were applied to motion prediction of two small AUVs designed for intervention tasks in shallow bathymetry. Mean forecast accuracy was assessed using Theil’s Inequality coefficient, which gives an agnostic measure of performance versus a naïve last-value predictor. Furthermore, the confidence interval of forecasts was studied for each model to determine if they could adequately capture uncertainty about the future state. The results of the analysis showed that a GP combined with a Periodic Rational Quadratic kernel function performed better than both the time series and Fourier series for the short-term, mean prediction error. And for the long-term, the GP had comparable performance to time series. Analysis of the forecast uncertainty revealed new and interesting insights. The time series method had low forecast uncertainty, as noted in prior research, and thus was overconfident in its prediction estimate. Conversely, the Fourier series was the most conservative of the three methods. The GP method was overconfident in the short term, but provided the best confidence bounds for long-term forecasts.

where: ( ) ( ) 𝛿 = sin 1∕2klx cos(𝜃) sin 1∕2kly sin(𝜃) sin(𝜔t)

(32)

and: • • • • • •

d is the draught (depth of submersion) for the vessel, g is the gravitational acceleration, h is the wavelength, H is the wave height, k is a constant, lx , ly is the length and width of the vessel respectively,

• 𝜌 is the fluid density, • 𝜃 is the angle of wave propagation w.r.t. the vehicle, and • 𝜔 is the wave frequency. The force being generated is a function of the vehicle heading relative to the wave propagation, 𝜃. Fig. 10 shows a diagram comparing the SPIR and the A-AUV in the x–y plane. The SPIR’s centre of mass does not align with the centre of volume. As such, any net wave forces on the body will cause a yaw moment and a subsequent change in 𝜃 over time. Fig. 11 shows how there is an underlying curvature to the yaw motion data in Fig. 8. Furthermore, if the angle of the wave propagation with respect to the robot, 𝜃, changes over time because of said yaw moment, then the surge and sway forces will also drift with time. This may account for the nonstationarity of the data.

Acknowledgements This work is supported in part by the Australian Research Council (ARC) Linkage Project (LP150100935), Roads and Maritime Services of NSW, Australia and the Centre for Autonomous Systems (CAS) at the University of Technology Sydney, Australia. The authors would like to thank Dr. Andrew To and Dr. Khoa Le for their generous assistance with data collection.

Conversely, the body of the A-AUV is cylindrical in shape and will thus not incur any large yaw moments from passing waves. The amplitude of the yaw motion data is mostly smaller than that of the SPIR. Moreover, it has a constant, zero mean (Fig. 5) and hence exhibits stationarity. 12

Ocean Engineering 190 (2019) 106385

J. Woolfrey et al.

References

Li, X., Zhao, Z., Jilkov, V., 2001. Practical measures and test for credibility of an estimator. In: Proc. Workshop on Estimation, Tracking, and Fusion—A Tribute to Yaakov Bar-Shalom. pp. 481–495. Love, L.J., Jansen, J.F., Pin, F.G., 2004. On the modeling of robots operating on ships. In: IEEE International Conference on Robotics and Automation, 2004 Proceedings. Vol. 3. ICRA’04. 2004, IEEE, pp. 2436–2443. Medagoda, L., Williams, S., 2012. Model predictive control of an autonomous underwater vehicle in an in situ estimated water current profile. In: OCEANS. pp. 1–8. Parker, P., Anderson, B., 1990. Frequency tracking of nonsinusoidal periodic signals in noise. Signal Process. 20, 127–152. Rasmussen, C.E., Williams, C.K.I., 2006. Gaussian Processes for Machine Learning. MIT Press. Riedel, J., Healy, A., 1998. Shallow water station keeping of AUVS using multi-sensor fusion for wave disturbance prediction and compensation. In: OCEANS 98. pp. 1064–1068. Roberts, S., Osborne, M., Ebden, M., Reece, S., Gibson, N., Aigrain, S., 2013. Gaussian processes for timeseries modelling. Philos. Trans. R. Soc. A 371 (1984), 20110550. Schwarz, G., 1978. Estimating the dimension of a model. Ann. Statist. 6 (2), 461–464. Silver, N., 2012. The Signal and the Noise: Why So Many Predictions Fail - But Some Don’T. The Penguin Press, Penguin Group (USA) Inc., 375 Hudson Street, New York 10014, U.S.A, p. 5. Theil, H., 1966. Applied Economic Forecasting. North-Holland Publishing Company, Amsterdam, Rand McNally & Company, Chicago. Woolfrey, J., Liu, D.K., Carmichael, M., 2016. Kinematic control of an autonomous underwater vehicle-manipulator system (AUVMS) using autoregressive prediction of vehicle motion and model predictive control. In: IEEE International Conference on Robotics and Automation, ICRA. Yan, J., Li, K., Bai, E., Yang, Z., Foley, A., 2016. Time series wind power forecasting based on variant Gaussian process and TLBO. Neurocomputing 189, 135–144. Yang, X., Pota, H., Garratt, M., Ugrinovskii, V., 2008. Ship motion prediction for maritime flight operations. IFAC Proc. Vol. 41 (2), 12407–12412. Yuen, S., Novotny, P., Howe, R., 2008. Quasiperiodic predictive filtering for robotassisted beating heart surgery. In: IEEE International Conference on Robotics & Automation, ICRA. pp. 3875–3880.

Becker, B., Tummala, H., Riviere, C., 2008. Autoregressive modeling of physiological tremor under microsurgical conditions. In: 30th Annual International IEEE EMBS Conference. pp. 1948–1951. Bowthorpe, M., Tavakoli, M., 2016. Generalized predictive control of a surgical robot for beating-heart surgery under delayed and slowly-sampled ultrasound image data. In: IEEE Robotics and Automation Letters (RA-L) Papers, presented at ICRA 2016. Box, G., Jenkins, G., Riensel, G., Ljung, G., 2016. Time Series Analysis: Forecasting and Control, fifth ed. John Wiley & Sons, Inc., Hoboken, New Jersey. Brahim-Belhouari, S., Bermak, A., 2004. Gaussian process for nonstationary time series prediction. Comput. Statist. Data Anal. 47, 705–712. Chung, J.C., Bien, Z., Kim, Y.S., 1990. A note on ship-motion prediction based on wave-excitation input estimation. IEEE J. Ocean. Eng. 15 (3), 244–250. Dean, R., Dalrymple, R., 1991. Water Wave Mechanics for Engineers and Scientists. In: Advanced Series in Ocean Engineering, World Scientific Publishing Co Pty Ltd., 5 Toh Tuck Link, Singapore, 596224. Dürichen, R., Pimentel, M., Clifton, L., Schweikard, A., Clifton, D., 2015. Multitask Gaussian process for multivariate phsyiological time-series analysis. IEEE Trans. Biomed. Eng. 62 (1), 314–322. Fang, C., Williams, C., 2014. General probabilistic bounds for trajectories using only mean and variance. In: IEEE International Conference on Robotics & Automation, ICRA. pp. 2501–2506. Fernandéz, D., Hollinger, G., 2016. Model predictive control for underwater robots in ocean waves. In: IEEE Interncational Conference on Robotics and Automation, ICRA. From, P.J., Duindam, V., Gravdahl, J.T., Sastry, S., 2009. Modeling and motion planning for mechanisms on a non-inertial base. In: 2009 IEEE International Conference on Robotics and Automation. IEEE, pp. 3320–3326. From, P., Gravahl, J., Lillehagen, T., Abbeel, P., 2011. Motion planning and control of robotic manipulators on seasborne platforms. Control Eng. Pract. 19, 809–819. Huynh, V.T., Dunbabin, M., Smith, R.N., 2015. Predictive motion planning for AUVs subject to strong time-varying currents and forecasting uncertainties. In: 2015 IEEE International Conference on Robotics and Automation, ICRA. IEEE, pp. 1144–1151. Jones, D., Hollinger, G.A., 2017. Planning energy-efficient trajectories in strong disturbances. IEEE Robot. Autom. Lett. 2 (4), 2080–2087. Kormushev, P., Caldwell, D., 2013. Improving the energy efficiency of autonomous underwater vehicles by learning to model disturbances. In: IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS. pp. 3885–3892.

13