Over the last decade, it has become clear that autonomous marine vehicles will revolutionize ocean sampling. Several researchers have investigated approaches for ASVs and AUVs for adaptive ocean sampling Yuh (2000)-Smith et al. (2010c) and fundamental marine sampling techniques for ASVs and AUVs are discussed in Singh et al. (1997). Besides control algorithms for Oceanic Sampling, an alternative approach is to use static sensor placements to maximize information gathering Zhang and Sukhatme (2008).

Our work connects also with research on control design for AUVs for adaptive ocean sampling, Yoerger and Slotine (1985); Frazzoli et al. (2002); Low et al. (2009); Rudnick and Perry (2003); Yuh (2000); Frank and Jónsson (2003); Graver (2005); Barnett et al. (1996); Carreras et al. (2000); Ridao et al. (2000); Rosenblatt et al. (2002); Turner and Stevenson (1991); Whitcomb et al. (1999, 1998); McGann et al. (2008b), McGann et al. (2008a), McGann et al. (2008c), Yoerger and Slotine (1985)-McGann et al. (2008c). Applications of ocean sampling techniques for autonomous vehicles are discussed in Singh et al. (1997)-Eriksen et al. (2001). This body of research differs from the proposed research in that we plan to utilize predictive models in the form of Partial Differential Equations (PDE) to enable effective sampling, navigation, and localization within dynamic features.

Reinforcement learning in marine robotics, especially model-free methods, is an attractive alternative to finding plans for several reasons. First, executing marine robotics experiments and deployments is expensive, time-consuming, and often risky; controllers learned through RL can represent significant time and cost savings and shorten the time to deployment. Second, system identification can sometimes be challenging in marine environments due to several factors such as unmodeled dynamics and environment’s unknowns; for that reason, model-free RL approaches can be an alternative in these scenarios. Examples of approaches that have used RL for ASVs or AUVs include path planning Yoo and Kim (2016), control Cui et al. (2017) and tracking Martinsen et al. (2020).

Our ideas are also connected to the use of Machine Learning models in the context of Partial Differential Equations. Due to their usefulness and impact in several domains, there have been efforts to use modern machine learning techniques to solve high dimensional PDEs Han et al. (2018), find appropriate discretizations Han et al. (2018), and control them Farahmand et al. (2017).

Partial differential equations (PDEs) have been used to model water features of interest such as pH, temperature, turbidity, salinity, and chlorophyll-A. Depending on the nature of their motion, they can be modeled through diffusion, advection or a combination of both. It is important to evaluate how they behave given certain initial conditions to understand their evolution in time. We model the ocean features of interest as a scalar field f : R 2 × [ 0 , ∞ ) → R .

Once we chose a PDE as a model, it is crucial to estimate the parameters of the PDE to get a reliable model. This problem belongs to the family of inverse problems since those parameters are sensitive to the observations and given initial conditions Richard et al. (2019) Antman et al. (2006). Because of this sensitivity, it is computationally expensive to find the PDE parameters. There are optimization-based techniques to solve this problem. These techniques problems balance the fitting parameter to the observations and the model sensitivity to those parameters. One of most used methods is the Tikhonov regularization technique Nair and Roy (2020) Bourgeois and Recoquillay (2018). It comprises solving a regularized optimization problem to get a regularized solution. It can be highly efficient depending on the regularization norm (especially if the L ² norm is used). However, it depends on the regularization constant to achieve good results.

Other approaches to solving the PDE estimation problem take advantage of Bayesian theory Xun et al. (2013). In this case, bayesian learning is connected to regularization since the regularization problem coincides with the maximization of the likelihood of the parameters given the observations Bishop (1995). Therefore, Machine Learning techniques have been proposed to take advantage of the capability of the models to discover hidden relationships between the input data and the final estimation Jamili and Dua (2021). Most of those models use the fact that the samples are given in advance. This work proposes a learning mechanism to select samples that can reasonably estimate the model without exploring the complete domain. This principle has been used in numerical integration problems resulting in several quadrature rules, such that Gauss–Kronrod, Gauss-Legendre, or Newton cotes Kincaid et al. (2009). Those methods have proven to be more efficient since they can give reliable estimations using few points. In this work, we employ an intelligent agent capable of sampling the environment, searching for reliable samples, and using them to compute the parameters of a PDE. Also, this allows estimating the ocean feature behavior in the domain according to Eq. 2.

We modeled the marine environment as a 2-D water layer (representing, for example, the surface) denoted as W ⊂ R 2 where W is an open and bounded set. The obstacle-free state space for our robot is represented by S = W \ O , where O represents the set of locations that are not accessible to the robot.

To estimate the flow field, we define a scheme of fitting problems based on the known initial conditions of Eq. 1 h(x) and the current samples acquired by the agent. First, we expect to collect samples y _i at the location x _i and time t _i for i = 1, … , n such that the field minimizes the mean square error of the collected samples. Taking advantage of the closed solution described in the homogeneous version of Eq. 1, the fitting error function e _f ( b ) is defined as e f b ; x 1 , … , x n = 1 n ∑ i = 1 n f x i , t i − y i 2 = 1 n ∑ i = 1 n h x i − t i b − y i 2 (3) where n is the number of collected samples. Next, we define the fitting error e _f (x ₁, … , x _n ) associated to the locations e f x 1 , … , x n = min b ∈ R 2 e f b ; x 1 , … , x n . (4)

The fitting error expressed in Eq. 4 measures how well can the best fitted model prediction of the given samples (i.e., predict y _i given x _i and a parameter vector b . It is the “best” in the sense that is the minimum achievable error produced by the model given the samples x ₁, … , x _n ). Nevertheless, we can notice that if the locations x ₁, … , x _n are wrongly chosen, the fitting error can be low, but its capability of estimating the entire field may lead to over-fitting problems. To handle this issue, we add a new error term based on how well one sample can be predicted using the remaining ones. This is known as cross-validation. In this case, we propose the following cross-validation scheme. For each 1 ≤ i ≤ n let b i ∗ defined as b i ∗ = a r g min b ∈ R 2 e f b ; x 1 , … , x i − 1 , x i + 1 , … , x n . (5)

We define the cross validation error e _cv (x ₁, … , x _n ) as e c v x 1 , … , x n = 1 n ∑ i = 1 n h x i − t i b i ∗ − y i 2 , (6) and it measures on average how well the samples can fit a model, which is estimating the remaining sample. This avoids the over-fitting problems and allows to measure how reliable are the taken samples. Lastly, we define the total error e _total (x ₁, … , x _n ) or just e _total as e total x 1 , … , x n = e f x 1 , … , x n + e c v x 1 , … , x n . (7)

This error compound aims to have an equal trade off between the sample estimation measured by e _f (x ₁, … , x _n ) and the reliability of the samples measured by e _cv (x ₁, … , x _n ).

The agent is modeled as a rigid body that moves in R 2 and can be described by a non-linear system as x ̇ = f x , u z = o x , r (8)

Such that f (x, u) is the motion model of the vehicle, o (x, r) is the observation model of the vehicle, and r are additive, zero-mean noise to account for modeling errors and sensor imperfections.

Let S be the state space, i e., the set of all possible states x ∈ S and U be the action space, which represents the set of all possible actions. Therefore, a configuration of the vehicle can be described by x = x , y , ϕ u = u v , u ω (9)

In which (x, y) is the position of the vehicle and ϕ ∈ (−π/4, π/4) is the vehicle’s heading; the forward speed v and the angular velocity of the agent orientation ω can be set directly by the action variables u _v and u _ω , respectively. The kinematic model of the agent x ̇ = f ( x , u ) is described by Eq. 10. x ̇ = u v ⁡ cos ⁡ ϕ + v x y ̇ = u v ⁡ sin ⁡ ϕ + v y ϕ ̇ = u ω (10) where v _x and v _y account for the velocity components of the environment (flow field) in x and y directions.

Let x S ∈ S be the initial location of the agent. It is assumed that the agent takes advantage of ocean current dynamics as it drifts and moves forward with or against the currents and rotates clockwise or counterclockwise. Therefore, the action space is defined as U = 0 , v max / 2 , v max × ϕ min , ϕ max (11)

We discretize the action space to obtain a finite subset of U defined as A = v max , 0 , v max / 2,0 , v max , − ϕ , v max , + ϕ , v max / 2 , − ϕ , v max / 2 , + ϕ , 0,0 (12)

The description of the actions of the agent are summarized in Table 1.

For the observation model, we assume that the vehicle uses an IMU to measure its heading angle ϕ and has access to GPS at surface level. Also, the vehicle can observe its state with uncertainties due to sensor imperfections and the dynamic nature of the underwater environment. The observation space Z , the set of all possible sensor observations z ∈ Z , is given by Z = x min , x max × y min , y max × ϕ min , ϕ max . (13)

The observation model o(x) is represented by z = o x , r = I x + r (14)

Where r ∈ R 3 is noise distributed as r ∼ N ( 0 , Σ ) with Σ a diagonal covariance matrix to account for modeling errors and sensor imperfections and I is the identity matrix. It was also considered that the measurement noises of each sensor are uncorrelated and have constant covariance.

These elements allow us to formulate the following problem.

Problem: Given an aquatic environment W , the action set of the agent A , the state space S = W × ( ϕ min , ϕ max ) , the vehicle’s motion model, observations of a given ocean feature in several locations, estimate the flow field (and therefore the ocean feature distribution) by minimizing cross-validation error and error fitting within a given fixed number of steps.

In reinforcement learning problems, designing a reward function is not a trivial task, Ng et al. (1999). To avoid spurious exploration, we defined a terminal condition with a fixed number of observations taken to determine when to reset the environment for a new episode. To encourage the agent to minimize the fitting and cross-validation errors within a given number of steps, we provide a reward that is inversely proportional to the sum of the errors at the terminal state. For each episode, the agent collects 20 observations, and the reward function is defined as r s , a = 100 e total , if the number of observations = 20 c 2 , if the number of observations < 5 c 1 c 2 1.5 , otherwise (19) where c ₁ is the ratio between total error at previous and current step and c 2 = 1 − e total 20 0.4 .

In reinforcement learning systems, feature construction is critical since it values each state of the agent. The main techniques for feature construction of linear methods are polynomial-based, Fourier basis and tile coding Sherstov and Stone (2005). As such, tile coding is a computationally effective feature design technique that divides the state space into divisions called tiles. Each element in the tiling is referred to as a tile. Different tilings are separated by a fixed-size fraction of the tile width Sutton and Barto (2018). If there are n tilings and each tiling has m × m tiles, the feature vector is x ( s ) ∈ R n × m × m . One of the main advantages of using tile coding with binary feature vectors is that the weighted sum in the approximate value function (Eq. 17) is easy to compute. Figure 2 shows an example of the representation of tile coding for two-dimensional continuous state space. In this case, x(s) is a feature vector with twelve components, one for each tile in each tiling. Each component of x(s) is inactive (zero-valued) except active components x ₀(s), x ₄(s) and x ₈(s) that corresponds to the current location states of the agent. As a consequence, there are n active features in x(s) because every position in state space falls into precisely one tile in each of the n tilings. Let the weight vector w = [ w 0 , … , w 11 ] ⊤ and the action space be A = {a ₀, a ₁, a ₂}. The feature vector regarding actions a ₀, a ₁ and a ₂ is x (s, a ₀) = x (s, a ₁) = x (s, a ₂) = [1,0,0,0,1,0,0,0,1,0,0,0] ^⊤ . Thus, the action-value function approximation q ̂ ( s , a , w ) described in Eq. 17 is computed as q ̂ s , a , w = ∑ i = 1 d w i (20) for each action in action space.

With tile coding, design issues for discrimination and generalization should be taken into account. The number and size of tiles, for example, affect the granularity of state discrimination, or how far the agent must move in state space to change at least one component of the feature vector. Aside from that, the shape of the tilings and the offset distance between them have an impact on generalization. As an example, if tiles are stretched along one dimension in state space, generalization will extend to states along that dimension as well Sutton and Barto (2018).

In problems with large state spaces, the eligibility trace is a technique to promote computational efficiency of reinforcement learning methods. The eligibility trace is a vector z t ∈ R d whose components maintain track of which components of the weight vector w _t have contributed to recent state values and temporarily records the occurrence of estimated events. Therefore, components of w _t that most frequently contribute to valuations of previous states are considered eligible for an update Singh et al. (1995). Eligibility trace components are updated based on the trace-decay parameter λ ∈ [0, 1], which specifies the pace at which the trace fades away exponentially. In contrast with n-step methods that perform action-value updates after a given number of steps, eligibility traces provide updates continually over the learning process. For this reason, agent behavior can be modified right after a new state is found rather than being delayed n steps.

The action-value return function G _t is a function approximation of the n-step return defined as G t : t + n = r t + 1 + ⋯ + γ n − 1 q ̂ s t + n , a t + n , w t + n − 1 , t + n < T (21) where γ is the discount rate that regulates the relative importance of near-sighted and far-sighted rewards. Thus, the λ-return G t λ is written as G t λ = 1 − λ ∑ n = 1 T − t − 1 λ n − 1 G t : t + n + λ T − t − 1 G t (22)

In this way, the update rule for the weight vector in Eq. 16 is modified as follows w t + 1 = w t + α G t λ − q ̂ s t , a t , w t ∇ q ̂ s t , a t , w t = w t + α δ t z t (23) where the action-value estimation error δ _t is defined as δ t = r t + 1 + γ q ̂ s t + 1 , a t + 1 , w t − q ̂ s t , a t , w t (24)

The action-value representation of the eligibility trace is defined as z − 1 = 0 z t = γ λ z t − 1 + ▽ q ̂ s t , a t , w t , 0 ≤ t ≤ T (25)

The complete algorithm for SARSA(λ) is presented in table 1, Sutton and Barto (2018).

Algorithm 1

SARSA(λ) with linear function approximation.

For future work, we consider the expansion of the proposed method for 3D environments. This can be accomplished by augmenting the vehicle model (state space, action space, observation space) and validating the proposed framework with deployments in aquatic environments such as in the Biscayne Bay area, Florida, United States. Besides that, it is possible to refine our estimation strategies with cooperative agents. A primary direction for future work is to incorporate a combination of heterogeneous agents in order to provide better estimates of the locations of the ocean feature. In this work, we assume known initial conditions for a given linear, constant flow field. A second direction for future work is to investigate how effective the proposed estimation framework is for time-varying flow fields and actual oceanic data from the Regional Ocean Modeling System (ROMS) Shchepetkin and McWilliams (2005). ROMS data set that provides current velocity prediction data consisting of three spatial dimensions (longitude, latitude, and depth) associated with time. Finally, tracking oceanic features, such as the Lagrangian coherent structures (LCS) contributes to a wide range of applications in ocean exploration Hsieh et al. (2012). Therefore, an additional direction of our work is to expand the current work towards an efficient method for LCS tracking using machine learning techniques.