Model-Based Learning of Turbulent Flows using Mobile Robots

12/10/2018 ∙ by Reza Khodayi-mehr, et al. ∙ Duke University 12

In this paper we consider the problem of model-based learning of turbulent flows using mobile robots. The key idea is to use empirical data to improve on numerical estimates of time-averaged flow properties that can be obtained using Reynolds-Averaged Navier Stokes (RANS) models. RANS models are computationally efficient and provide global knowledge of the flow but they also rely on simplifying assumptions and require experimental validation. In this paper, we instead construct statistical models of the flow properties using Gaussian Processes (GPs) and rely on the numerical solutions obtained from RANS models to inform their mean. We then utilize Bayesian inference to incorporate empirical measurements of the flow into these GPs, specifically, measurements of the time-averaged velocity and turbulent intensity fields. Our formulation accounts for model ambiguity and parameter uncertainty via hierarchical model selection. Moreover, it accounts for measurement noise by systematically incorporating it in the GP models. To obtain the velocity and turbulent intensity measurements, we design a cost-effective mobile robot sensor that collects and analyzes instantaneous velocity readings. We control this mobile robot through a sequence of waypoints that maximize the information content of the corresponding measurements. The end result is a posterior distribution of the flow field that better approximates the real flow and also quantifies the uncertainty in the flow properties. We present experimental results that demonstrate considerable improvement in the prediction of the flow properties compared to pure numerical simulations.



There are no comments yet.


page 1

page 7

page 11

page 12

page 14

page 16

page 18

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Knowledge of turbulent flow properties, e.g., velocity and turbulent intensity, is of paramount importance for many engineering applications. At larger scales, these properties are used for the study of ocean currents and their effects on aquatic life, [1, 2, 3], meteorology, [4], bathymetry, [5], and localization of atmospheric pollutants, [6], to name a few. At smaller scales, knowledge of flow fields is important in applications ranging from optimal HVAC of residential buildings for human comfort, [7], to design of drag-efficient bodies in aerospace and automotive industries, [8]. At even smaller scales, the characteristics of velocity fluctuations in vessels are important for vascular pathology and diagnosis, [9] or for the control of bacteria-inspired uniflagellar robots, [10]. Another important application that requires global knowledge of the velocity field is chemical source identification in advection-diffusion transport systems, [11, 12, 13].

The most cost-effective way to estimate a flow field is using numerical simulation and the Navier-Stokes (NS) equations. These equations are a set of Partial Differential Equations (PDEs) derived from first principles, i.e., conservation of mass, momentum, and energy, that describe the relationship between flow properties,

[14]. They can be simplified depending on the nature of the flow, e.g., laminar versus turbulent or incompressible versus compressible, to yield simpler, more manageable equations. Nevertheless, when the flow is turbulent, solving the NS equations becomes very challenging due to the various time and length scales that are present in the problem. Turbulent flow is characterized by random fluctuations in the flow properties and is identified by high Reynolds numbers. In many of the applications discussed before, especially when air is the fluid of interest, even for very small velocities the flow is turbulent. In this case, Direct Numerical Solution (DNS) of the NS equations is still limited to simple problems. A more effective approach is to solve the time-averaged models called Reynolds-Averaged Navier Stokes equations, or RANS for short. RANS models are useful since in many engineering applications we are interested in averaged properties of the flow, e.g., the mean velocity components, and not the instantaneous fluctuations. Note that despite its chaotic nature, turbulent flow has structure; it consists of eddies that are responsible for generation and dissipation of kinetic energy in a cascade that begins in large or integral scales and is dissipated to heat in small scales called Kolmogorov scales. Understanding these structures is essential for proper modeling of turbulence; see [15] for more details.

Averaging the flow properties in RANS models introduces additional unknowns, called Reynolds stresses, in the momentum equation that require proper modeling to complete the set of equations that are solved. This is known as the ‘closure problem’. Depending on how this problem is addressed, various RANS models exist that are categorized into two large classes, Eddy Viscosity Models (EVMs) and Reynolds Stress Models (RSMs); see [15] for more details. The former assume that the Reynolds stress is proportional to the strain defining, in this way, an effective turbulent viscosity, while in the latter a transport equation is derived for the Reynolds stress components and is solved together with the averaged flow equations. EVMs have an assumption of isotropy built-into them meaning that the Reynolds stresses are direction-independent. Although this assumption is not valid for flow fields that have strong asymmetries, there are many benchmark problems for which the results from EVMs are in better agreement with empirical data than RSMs; see e.g. [16]. In general, solutions provided by these RANS models are not necessarily compatible with each other and more importantly with the real world and require experimental validation, [17]. Furthermore, precise knowledge of the Boundary Conditions (BCs) and domain geometry is often unavailable which contributes to even more inaccuracy in the predicted flow properties. Finally, appropriate meshing is necessary to properly capture the boundary layer and get convergent solutions from RANS models, which is not a straight-forward process.

In order to address the shortcomings of RANS models and complement the associated numerical methods, various experimental techniques have been developed that directly measure the flow properties. These methods can also be used as stand-alone approaches alleviating the need for accurate knowledge of the BCs and geometry. A first approach of this kind utilizes hot-wire constant temperature anemometers. A hot-wire anemometer is a very thin wire that is placed in the flow and can measure flow fluctuations up to very high frequencies. The fluctuations are directly related to the heat transfer rate from the wire, [18]. Hot-wire sensors are very cost effective but require placing equipment in the flow that might interfere with the pattern and is not always possible. A next generation of approaches are optical methods that seed particles in the flow and optically track them. Two of the most important variants of optical methods are Particle Imaging Velocimetry (PIV) and Laser-Doppler Anemometry (LDA). In PIV, two consecutive snapshots of the particles in the flow are obtained and their spatial correlation is used to determine the shift and the velocity field, [19]. On the other hand, in LDA a pair of intersecting beams create a fringe in the desired measurement location. The light from a particle passing through this fringe pulsates and the frequency of this pulse together with the fringe spacing determine the flow speed, [20]. Optical methods although non-intrusive, require seeding particles with proper density and are expensive. Another sensing approach utilizes Microelectromechanical Systems (MEMS) for flow measurements. Because of their extremely small size and low cost, these devices are a very attractive alternative to aforementioned methods, [21].

Empirical data from experiments are prone to noise and error and are not always available for various reasons, e.g., cost and accessibility, scale of the domain, and environmental conditions. On the other hand, as discussed before, numerical simulations are inaccurate and require validation but are cost-effective and can provide global information over the whole domain. In this paper we propose a model-based learning method that combines numerical solutions with empirical data to improve on the estimation accuracy of flow properties. Specifically, we construct a statistical model of the time-averaged flow properties using Gaussian Processes (GPs) and employ the physics of the problem, captured by RANS models, to inform their prior mean. We construct GPs for different RANS models and realizations of the random BCs and other uncertain parameters. We also derive expressions for measurement noise and systematically incorporate it into our GP models. To collect the empirical data, we use a custom-built mobile robot sensor that collects and analyzes instantaneous velocity readings and extracts time-averaged velocity and turbulent intensity measurements. Specifically, we propose an online path planning method that guides the robot through optimal measurement locations that maximize the amount of information collected about the flow field. Then, using Bayesian inference and hierarchical model selection, we obtain the posterior probabilities of the GPs and the flow properties for each model after conditioning on the empirical data. The proposed robotic approach to learning the properties of turbulent flows allows for fewer, more informative measurements to be collected. Therefore, it is more efficient, autonomous, and can be used for learning in inaccessible domains.

Ii Related Work

Ii-a Machine Learning in Fluid Dynamics

Machine Learning (ML) algorithms have been used to help characterize flow fields. One way to do this is by model reduction, which speeds up the numerical solutions, [22, 23, 24, 25]. Particularly, [22]

utilizes Neural Networks to reduce the order of the physical model and introduce a correction term in the Taylor expansion of the near wall equations that improves the approximations considerably.


uses Deep Learning to construct an efficient, stable, approximate projection operator that replaces the iterative solver needed in the numerical solution of the NS equations.

[24] proposes an alternative approach to traditional numerical methods utilizing a regression forest to predict the behavior of fluid particles over time. The features for the regression problem are selected by simulating the NS equations. A different class of methods utilize ML to improve accuracy instead of speed. For instance, [26] presents a modification term in the closure model of the RANS equations which is learned through data collected experimentally or using DNSs. The modified closure model replaces the original one in the RANS solvers. In a different approach, [17] utilizes classification algorithms to detect regions in which the assumptions of the RANS models are violated. They consider three of the main underlying assumptions: (i) the isotropy of the eddy viscosity, (ii) the linearity of the Boussinesq hypothesis, and (iii) the non-negativity of the eddy viscosity.

The approaches discussed above mainly focus on application of ML to improve on the speed or accuracy of the RANS models at the solver level without considering various sources of uncertainty that are not captured in these models, e.g., imperfect knowledge of the BCs and geometry. To the contrary, here we focus on validation, selection, and modification of the solutions obtained from RANS models using a data-driven approach. Particularly, we employ Gaussain Processes (GPs) to model the flow field and propose a supervised learning strategy to incorporate empirical data into numerical RANS solutions in real-time. GPs provide a natural way to incorporate measurements into numerical solutions and to construct posterior distributions in closed-form. Note that GPs have been previously used to model turbulent flows in other contexts. For instance,

[26] uses GPs to modify closure models and improve accuracy while [27] models turbulent fields via GPs and then use these GPs to obtain probabilistic models for the dispersion of particles in the flow field. Since turbulence is a non-Gaussian phenomenon, assuming Gaussianity is a major approximation, [28]

. Nevertheless, here we focus on mean flow properties which are Normally distributed regardless of the underlying distribution of the instantaneous properties,


Ii-B Active Learning for Gaussian Processes

Active learning using GPs has been investigated extensively in the robotics literature. This is because of the closed form of the posterior distribution that makes GPs ideal for inference and planning. Collecting informative data in active learning applications depends on defining appropriate optimality indices for planning with entropy and Mutual Information (MI) being the most popular ones. Particularly, MI is submodular and monotone and can be optimized using greedy approaches while ensuring sub-optimality lower-bounds; see [30]. Necessary conditions to obtain such sub-optimality lower bounds for MI are discussed in [31]. When used for planning, such optimality indices can guide mobile robots through sequences of measurements with maximum information content. For example, [32] proposes a suboptimal, non-greedy trajectory planning method for a GP using MI, while [2]

relies on Markov Decision Processes and MI to plan paths for autonomous under-water vehicles. On the other hand,

[33] collects measurements at the next most uncertain location which is equivalent to optimizing the entropy. In this paper, we employ a similar approach where the mobile sensor plans its path to minimize the entropy of unobserved areas in the flow.

Active learning using GPs has applications ranging from estimation of nonlinear dynamics to spatiotemporal fields. For instance, [33] presents an active learning approach to determine the region of attraction of a nonlinear dynamical system. The authors decompose the unknown dynamics into a known mean and an unknown perturbation which is modeled using a GP. Also, [34, 35] use similar approaches to estimate and control nonlinear dynamics using robust and model predictive control, respectively. [36] presents an active sensing approach using cameras to associate targets with their corresponding nonlinear dynamics where the target behaviors are modeled by Dirichlet process-GPs. On the other hand, [1] uses GPs to estimate the uncertainty in predictive ocean current models and rely on the resulting confidence measures to plan the path of an autonomous under-water vehicle. GPs have also been used for estimation of spatiotemporal fields, e.g., velocity and temperature. For instance, [37]

proposes a learning strategy based on subspace identification to learn the eigenvalues of a linear time-invariant dynamic system which defines the temporal evolution of a field.


addresses the problem of trajectory planning for a mobile sensor to estimate the state of a GP using a Kalman filter. Closely related are also methods for robotic state estimation and planning with Gaussian noise; see

[39, 40, 41].

The literature discussed above typically employs simple, explicit models of small dimensions and does not consider model ambiguity or parameter uncertainty. Instead, here we consider much more complex models of continuous flow fields that are implicit solutions of the NS equations, a system of PDEs. Solutions to the NS equations depend on uncertain parameters like BCs and geometry. We propose a model-based approach to learn the posterior distributions of these solutions along with the flow properties.

Ii-C Contributions

The contributions of this work can be summarized as follows: (i) To the best of our knowledge, this is the first model-based framework that can be used to learn flow fields using empirical data. To this date, ML methods have only been used to speed up numerical simulations using RANS models or to increase their accuracy. Instead, the proposed framework learns solutions of the NS equations by combining theoretical models and empirical data in a way that allows for validation, selection, and more accurate prediction of the flow properties. Compared to empirical methods to estimate flow fields, this is the first time that robots are used for data collection, while compared to relevant active learning methods in robotics, the flow models we consider here are much more complex, uncertain, and contain ambiguities that can only be resolved using empirical validation. (ii) We design a cost-effective flow sensor mounted on a custom-built mobile robot that collects and analyzes instantaneous velocity readings. We also propose a planning algorithm for the mobile robot based on the entropy optimality index to collect measurements with maximum information content and show that it out-performs similar planning methods that possess sub-optimality bounds. (iii) The proposed framework is able to select the most accurate models, obtain the posterior distribution of the flow properties and, most importantly, predict these properties at new locations with reasonable uncertainty bounds. We experimentally demonstrate the predictive ability of the posterior field and show that it outperforms prior numerical solutions.

The remainder of the paper is organized as follows. In Section III, we discuss properties of turbulent flows and develop a statistical model to capture the velocity field and corresponding turbulent intensity. Section IV describes the design of the mobile sensor, the proposed algorithms to collect and process the instantaneous velocity readings, the formulation of the measurement errors, and finally the proposed planning method for the mobile sensor to collect flow measurements. In Section V, we present experiments that demonstrate the ability of the posterior field to better predict the flow properties at unobserved locations. Finally, Section VI concludes the paper.

Iii Statistical Model of Turbulent Flows

In this section, we discuss the theoretical aspects of turbulence that are needed to model the flow field. We also present a statistical model of turbulent flows and discuss its components.

Iii-a Reynolds-Averaged Navier Stokes Models

Since turbulence is a 3D phenomenon, consider a 3D domain of interest

and the velocity vector field

defined over this domain where and . Following the analysis used in RANS models, we decompose the velocity field into its time-averaged value and a perturbation as



In order for this limit to exist, we assume that the flow is ‘stationary’ or ‘statistically steady’. Under this condition, the integral will be independent of .

Since turbulent flows are inherently random, we can view or equivalently as a Random Vector (RV). Let denote -th realization of this RV at location and time and consider the ensemble average defined as

where is the total number of realizations. Assuming ‘ergodicity’, is independent of time so that . This is a common assumption that allows us to obtain time-series measurements of the velocity field at a given location and use the time-averaged velocity as a surrogate for the ensemble average .

Next, consider the Root Mean Square (RMS) value of the fluctuations of the -th velocity component defined as

Noting that and using the ergodicity assumption, we have


denotes the variance of a Random Variable (RV). The average normalized RMS value is called turbulent intensity and is defined as


where is the normalizing velocity. This value is an isotropic measure of the average fluctuations caused by turbulence at point .

Finally, we consider the integral time scale of the turbulent field along the -th velocity direction, defined as


where is the autocorrelation of the velocity given by

and is the autocorrelation lag. We use the integral time scale to determine the sampling frequency in Section III-C3. For more details on turbulence theory, see [15].

Iii-B Statistical Flow Model

The RANS models use the decomposition given in (1) to obtain turbulence models for the mean flow properties. As discussed in Section I, depending on how they address the closure problem, various RANS models exist. The numerical results returned by these models depend on the geometry of the domain, the values of the Boundary Conditions (BCs), and the model that is used, among other things. Since there is high uncertainty in modeling the geometry and BCs and since the RANS models may return inconsistent solutions, numerical methods often require experimental validation. On the other hand, empirical data is prone to noise and this uncertainty needs to be considered before deciding on a model that best matches the real flow field. Bayesian inference provides a natural approach to combine theoretical and empirical knowledge in a systematic manner. In this paper, we utilize Gaussian Processes (GPs) to model uncertain flow fileds; see [42]. Using GPs and Bayesian inference we can derive the posterior distributions of the flow parameters conditioned on available data in closed-form. Note that while turbulence is a non-Gaussian phenomenon, [28], the time-averaged velocity and turbulent intensity fields are Normally distributed regardless of the underlying distribution of the instantaneous velocity field. Therefore, we define GP models for these time-averaged quantities, whose knowledge happens to also be the most important for many engineering applications.

We construct the proposed GP models based on numerical solutions as opposed to purely statistical models, i.e., empirical data are used only to modify model-based solutions and not to replace them. Particularly, given a mean function and a kernel function , we denote a GP by . The mean function is given by the numerical solution from a RANS model and we define the kernel function as


where encapsulates the uncertainty of the field and is the correlation function. Here we define this function as a compactly supported polynomial


where is the correlation characteristic length and we define the operator . The correlation function (5) implies that, two points with distance larger than

are uncorrelated, which results in a sparse covariance matrix. In the following we use a constant standard deviation, i.e., we set


Consider a measurement model with additive Gaussian noise and let denote a measurement at location given by

Then is also a GP


with the following kernel function


where is the Dirac delta function. Given a vector of measurements at a set of locations , we can obtain the conditional distribution of (6) at a point

which is a Gaussian distribution

with mean and variance given by


where denotes the mean function evaluated at measurement locations and the entries of the covariance matrices and are computed using (7).

The GPs discussed above are defined individually for each one of the flow properties of interest, i.e., the 3D velocity field and the turbulent intensity. Specifically, denoting by the first mean velocity component at point , we can define the GP


for , with mean obtained from the numerical solution of a RANS model and the standard deviation in (4) defined as . This value is a measure of confidence in the numerical solution and can be adjusted to reflect this confidence depending on the model that is utilized and the convergence metrics provided by the RANS solver. GPs for the other velocity components and can be defined in a similar way. On the other hand, denoting by the turbulent intensity at point , we can define a GP


for , where, as before, we select the prior mean from the numerical solution of a RANS model and define the standard deviation in the kernel function (4) by a constant, i.e., .

Iii-C Measurement Model

Next, we discuss equation (6) used to obtain measurements of the velocity components and turbulent intensity. Specifically, at a given point , we collect a set of instantaneous velocity readings over a period of time and calculate the sample mean and variance. Note that by the ergodicity assumption, these time samples are equivalent to multiple realizations of the RV ; see Section III-A.

Iii-C1 Mean Velocity Components:

Consider the instantaneous first component of the velocity and a sensor model that has additive noise with standard deviation . An observation of this velocity at a given point and time is given by


where is the instantaneous measurement noise. Assume that we collect uncorrelated samples of at time instances for . Then, the sample mean of the first velocity component is given by


and is a random realization of the mean first velocity component .

Since the samples are uncorrelated, the variance of is given by


An unbiased estimator of

is the sample variance


Substituting this estimator in (13), we get an estimator for the variance of as


Note that as we increase the number of samples , becomes a more accurate estimator of .

Given the observation , we construct the measurement model for the mean first component of the velocity as


where is the perturbation with variance ; cf. equation (7). Note that the distribution of is Gaussian as long as the measurement noise in (11) is Gaussian; see [29]. In Section IV-B, we show that for typical off-the-shelf sensors, a Gaussian model is a good approximation for the measurement noise. Note that in addition to sensor noise , captured in the sample mean variance (14), other sources of uncertainty may contribute to . In Section IV-B, we derive an estimator for this variance taking into account these additional sources of uncertainty.

Next, we take a closer look at and the terms that contribute to it. Since the variance of the sum of uncorrelated RVs is the sum of their variances, from (11) we get . Furthermore, since the samples are uncorrelated, we have

and thus,




is the mean variance of the sensor noise. We derive an expression for the instantaneous variance in Section IV-B. Notice from equation (17) that the variance of is due to turbulent fluctuations and sensor noise . The former is a variation caused by turbulence and inherent to the RV whereas the latter is due to inaccuracy in sensors and further contributes to the uncertainty of the mean velocity component .

Iii-C2 Turbulent Intensity:

Referring to (2), turbulent intensity at a point is given as

Note that by the ergodicity assumption and similarly for the other velocity components; see Section III-A. Then, using equations (17) and (14), we have


where variances of the instantaneous velocity measurements and and the corresponding mean sensor noise variances and are computed using equations similar to (14) and (18), respectively.

Let denote the error in the turbulent intensity estimate. Then, we define the corresponding measurement as


Note that, unlike the mean velocity measurements, it is not straight-forward to theoretically obtain the variance of the turbulent intensity measurement. This is due to the nonlinearity in definition (2

). In general, estimating this variance requires the knowledge of higher order moments of the random velocity components; cf.

[29]. Specifically under a Gaussianity assumption, they depend on the mean values and which are not necessarily negligible. Consequently, we need to incorporate this uncertainty into our statistical model.

Here, we utilize the Bootstrap resampling method to directly estimate the variance . Assume that the samples are independent and consider the measurement set ; define and similarly. Furthermore, consider batches of size obtained by randomly drawing the same samples from , , and with replacement. Using (III-C2), we obtain estimates corresponding to the batches . Then, the desired variance can be estimated as


where is the mean of batch estimates ; see [43] for more details.

Iii-C3 Sampling Frequency:

In the preceding analysis, a key assumption was that the samples are independent. Since fluid particles cannot have unbounded acceleration, their velocities vary continuously and thus consecutive samples are dependent. As a general rule of thumb, in order for the samples to be independent, the interval between consecutive sample times and must be larger than where

is the maximum of the integral time scales (3) along different velocity directions at point .

In practice, we only have a finite set of samples of the instantaneous velocity vector over time and can only approximate (3). Let denote the discrete lag. Then, the sample autocorrelation of first velocity component is

This approximation becomes less accurate as increases since the number of samples used to calculate the summation decreases. Furthermore, the integral of the sample autocorrelation over the range is constant and equal to ; see [44]. This means that we cannot directly use (3) which requires integration over the whole time range. The most common approach is to integrate up to the first zero-crossing; see [45].

Iii-D Hierarchical Model Selection

Thus far, we have assumed that a single numerical solution is used to obtain the proposed statistical model. In this section, we show how to incorporate multiple numerical solutions in this statistical model. This is necessary since different RANS models may yield different and even contradictory approximations of the flow field. This effect is exacerbated by the uncertainty in the parameters contained in RANS models; see Section I. On the other hand, incorporating solutions from multiple RANS models in the proposed statistical model of the flow allows for a better informed prior mean and ultimately for model selection and validation using empirical data.

Specifically, consider a RV that collects all secondary variables that affect the numerical solution including uncertain BCs. Given an arbitrary distribution for each variable, we can construct a discrete distribution that approximates the continuous distributions using, e.g., a Stochastic Reduced Order Model (SROM); cf. [46] for details. Let denote the -th numerical solution corresponding to one combination of all these RVs and let the collection denote the discrete distribution, where is the number of discrete models. Given a vector of measurements at a set of locations , the posterior distribution over these discrete models can be easily obtained using Bayes’ rule as


where is the normalizing constant and is the likelihood of the data given model . Note that we have the marginal distributions in closed-form from the definition of the GP as

where and are short-hand notation for the mean and covariance of the GP constructed using model , see the discussion after equation (8). Noting that , we compute as


Given , we can easily compute the posterior distribution over the discrete models. Then, the desired mean velocity component fields and the turbulent intensity field are the marginal distributions , and after integrating over the discrete models. These marginal distributions are GP mixtures with their mean and variance given by


where denotes the posterior model probabilities obtained from (22). Equation (24b) follows from the principal that the variance of a RV is the mean of conditional variances plus the variance of the conditional means. The expressions for and are identical.

Iv Learning Flow Fields using Mobile Robot Sensors

In the previous section, we proposed a statistical model of turbulent flows based on GPs. Next, we develop a mobile robot sensor to measure the instantaneous velocity vector and extract the necessary measurements discussed in Section III-C. We also characterize the uncertainty associated with collecting these measurements by a mobile robot and formulate a path planning problem for it to maximize the information content of these measurements.

Iv-a Mobile Sensor Design

For simplicity, in what follows we focus on the in-plane components of the velocity field and take measurements on a plane parallel to the ground. Extension to the 3D case is straightforward. Particularly, we use the D6F-W01A1 flow sensor from OMRON; see Section V for more details. Figure 1 shows the directivity pattern for this sensor measured empirically.

Fig. 1: Directivity plot of the flow sensor used in the mobile robot where we measure the angle from the axis of the sensor. Above the measurements deviate from the theoretical prediction and above there is no detection. The bars on the experimental data show the standard deviation of the sensor samples caused by turbulent fluctuations and sensor noise.

In this figure, the normal velocity component is given as a function of the velocity magnitude and the angle of the flow from the axis of the sensor. It can be seen that there exists an angle above which the measurements deviate from the mathematical formula and for an even larger angle, no flow is detected at all. The numerical values for these angles are roughly and for the sensor used here. This directivity pattern can be used to determine the maximum angle spacing between sensors mounted on a plane perpendicular to the robot’s -axis, so that any velocity vector on the plane can be measured accurately. This spacing is , which means that eight sensors are sufficient to sweep the plane.

Figure 2 shows the mobile sensor comprised of these eight sensors separated by giving rise to a sensor rig that can completely cover the flow in the plane of measurements.

Fig. 2: Structure of the mobile sensor and arrangement of the flow velocity sensors. Eight sensors are placed in two rows in the center of the robot with an angle spacing of in each row and combined spacing of .

In this configuration, one sensor in the rig will have the highest reading at a given time with one of its neighboring sensors having the next highest reading. Let and denote the magnitude and angle of the instantaneous flow vector and and denote the -th sensor’s reading and heading angle in the global coordinate system where . Assume that sensors and have the two highest readings. Then, we have


These equations are used to get the four-quadrant angle in global coordinates. Given , the magnitude is given by . Assume that sensor has the highest reading at a given time, then and or and . Any deviation from these conditions indicates inaccurate readings due to sensor noise or delays caused by the serial communication between the sensors and microcontroller onboard; see Section V.

Algorithm 1 is used by the mobile robot sensor to collect the desired measurements.

0:  Measurement location , sample number , and the orientation of the mobile sensor;
1:  Given , compute sensor headings for ;
2:  for  do
3:     Receive the readings from all sensors;
4:     Let and denote sensors with the highest readings;
5:     if  then
6:        warning: inaccurate sample!
7:     end if
8:     Let and compute:
9:     if  or  then
10:        warning: inaccurate sample!
11:     end if
12:     Velocity magnitude:
13:     First velocity component: ;
14:     Second velocity component: ;
15:  end for
16:  First mean velocity component: ;
17:  Second mean velocity component: ;
18:  Compute sample variances and via (34) and (37);
19:  Compute turbulent intensity measurement via (III-C2);
20:  Compute variance of turbulent intensity using (21);
21:  Return , , , , , and ;
Algorithm 1 Flow Measurements using Mobile Sensor

Given the measurement location and the heading of the mobile sensor , the algorithm collects instantaneous samples of the velocity vector field and computes the mean velocity components and as well as the turbulent intensity . Particularly, in line 4 it finds the sensor indices and with the two highest readings, respectively. As discussed above, these two sensors should be next to each other. If not, the algorithm issues a warning indicating an inaccurate measurement in line 5. In line 8, it computes the flow angle using the four-quadrant inverse Tangent function; cf. equation (25). This flow angle is validated in line 9. In lines 12-14, the algorithm computes the magnitude of the velocity and its components. The vectors and store samples at different times, i.e., is the -th sample of the instantaneous velocity component ; see also equation (11). In lines 16 and 17, mean denotes the mean function. Finally, in lines 19 and 20, the algorithm computes the turbulent intensity measurement and its corresponding variance.

Iv-B Measurement Noise

In Section III-C, we defined the sensor noise for the instantaneous first velocity component by and the corresponding measurement noise by ; see equations (11) and (16). In the following, we derive an explicit expression for and discuss different sources of uncertainty that contribute to . Particularly, we consider the effect of heading error and location error . Since the relation between these noise terms and the measurements is nonlinear, we employ linearization so that the resulting distributions can be approximated by Gaussians. Generally speaking, let denote an observation which depends on a vector of uncertain parameters . Then, after linearization where , , and denotes the nominal value of uncertain parameters. For independent parameters , we have where is the constant diagonal covariance matrix. Then

by the properties of the Gaussian distribution for linear mappings. Specifically, the linearized variance of the measurement given the uncertainty in its parameters, can be calculated as


where denotes the variance of the uncertain parameter . Consequently, since the sensor noise and heading error are independent, we can sum their contributions to get .

First, we consider the sensor noise for the flow sensor discussed in Section IV-A. As before, we assume this noise is additive and Gaussian. For sensors with a full-scale accuracy rating, this assumption is reasonable as the standard deviation is fixed and independent of the signal value. Let FS denote the full-scale error and denote the standard deviation of the sensor noise. The cumulative probability of a Gaussian distribution within is almost equal to one. Thus, we set


In order to quantify the effect of sensor noise on the readings of the velocity components, referring to equation (25), we have


where depends on the angle spacing between sensors and and denote the sensor headings. Then, and


from (26), where and denote the nominal headings of the sensors in the global coordinate system. Similarly,


and . Then,


Note that equations (28) and (30) are linear in the instantaneous sensor readings. Thus, the contribution of sensor noise to uncertainty in the instantaneous velocity readings is independent of those readings and only depends on the robot heading .

Next, we consider uncertainty caused by localization error, i.e., error in the heading and position of the mobile robot. This error is due to structural and actuation imprecisions and can be partially compensated for by relying on motion capture systems to correct the nominal values; see Section V. Specifically, let denote the distribution of the robot heading around the nominal value obtained from a motion capture system. Note that where is the relative angle of sensor in the local coordinate system of the robot; the same is true for . Then, from (28) we have


where the last equality holds from (30). Evaluating (29) at the nominal values and and using (26), the contribution of the heading error to uncertainty in the measurement of the instantaneous first velocity component becomes


Following an argument similar to what used to derive (17), to obtain an expression for , and noting that the contributions of independent sources of uncertainty, i.e., sensor and heading errors, to the variance of the first mean velocity component are additive, cf. (26), we obtain an estimator for as


In (34), is given by (15) and can be obtained using an equation similar to (18). Note that unlike whose contribution to the uncertainty in can be reduced by increasing the sample size , the contribution of the heading error is independent of and can only be reduced by collecting multiple independent measurements. Similarly, for the second component of the velocity


Then using (26) and (35), we have




Finally, let denote a D distribution for the measurement location where is the nominal location and denotes its SROM discretization; see [46] for details. Without loss of generality, we assume , i.e., belongs to the set of SROM samples. From (6), the measurement at a point is Normally distributed as . Then, we can marginalize the location to obtain the expected measurement distribution as

This distribution is a Gaussian Mixture (GM) and we cannot properly model it as a normal distribution. Nevertheless, the probability of the measurement location being far from the mean (nominal) value drops exponentially. This means that corresponding to is larger than the rest of the weights and the distribution is close to a unimodal distribution. Noting that we can obtain the mean and variance of in closed-form, a Gaussian distribution can match up to two moments of the underlying GM. Particularly,


Considering equations (38), the following points are relevant: (i) For simplicity, we do not consider covariance between different measurements in this computation. This is reasonable as long as , where is the characteristic length of the correlation function (5). (ii) Equations (38) are only relevant at measurement locations. Thus, in equations (8), we replace the entry of at location with the corresponding mean value from (38a) and similarly, the diagonal entry of corresponding to with from (38b). (iii) The kernel function, defined in (7), depends on the heading error variances and of the velocity components through (34) and (37). These values depend on empirical data that are only available at the true measurement location and not at SROM samples . Thus, we use the same value for and at all SROM samples.

Iv-C Optimal Path Planning

Next, we formulate a path planning problem for the mobile sensor to collect measurements with maximum information content. Specifically, consider the domain and its D projection on the plane of mobile sensor discussed in Section IV-A. Let denote a discretization of and denote a discrete subset of points from that the mobile sensor can collect a measurement from. Our goal is to select measurement locations from the set , which we collect in the set , so that the entropy of the velocity components at unobserved locations , given the measurements , is minimized. With a slight abuse of notation, let denote this entropy. Then, we are interested in the following optimization problem

where denotes the cardinality of the measurement set . Noting that , we can equivalently write


The optimization problem (39) is combinatorial and NP-complete; see [31]. This makes finding the optimal set computationally expensive as the size of the set grows. Thus, in the following we propose a greedy approach that selects the next best measurement location given the current ones. Particularly, let denote the set of currently selected measurement locations, where

. By the chain rule of the entropy

Then, given , we can find the next best measurement location by solving the following optimization problem


To obtain the objective in (40), we note that for a continuous RV , the differential entropy is defined as

For a Gaussian RV , the value of the entropy is independent of the mean and is given in closed-form as


where and is defined in (8b). Particularly, here we consider the uncertainty in the velocity components and . Since these components are independent, we have

In order to predict the mean measurement variances and at a candidate location , we utilize the posterior turbulent intensity field given the current measurements, i.e., . To do so, we assume that the turbulent flow is isotropic, meaning that the variations are direction-independent. Then, the expression for the entropy defined in (41) can be approximated by

where in (8b) we approximate

see Section III-C1 for details.

Recall now the hierarchical model selection process, discussed in Section III-D, that allows to account for more than one numerical model. In this case, the posterior distribution is a GM, given by (24), and a closed-form expression for the entropy is not available; see [47]. Instead, we can optimize the expected entropy. Particularly, let denote the posterior probability of model , given the current measurements , obtained from (22). Then, the optimization problem (40) becomes


where is the added information between consecutive measurements corresponding to model and is given by


for entropy. Using this closed-form expression for the added information, the planning problem (42) can be solved very efficiently. Note that this problem depends on the current measurements through the posterior probabilities of the models as well as the posterior turbulent intensity field and must be solved online.

Algorithm 2 summarizes the proposed solution to the next best measurement problem (42) at step .

0:  Covariance matrix , the sets and , and current model probabilities ;
1:  for