Abstract:
A collaborative control method for tracking Lagrangian coherent structures (LCSs) and manifolds on flows employs at least three autonomous underwater vehicles (AUVs) each equipped with a local flow sensor. A first flow sensor is a tracking sensor and the other sensors are herding sensors for controlling and determining the actions of the tracking sensor. The AUVs are deployed in a body of water whereby the tracking sensor is positioned with respect to the herding sensors such that the herding sensors maintain a straddle formation across a boundary. A local flow velocity measurement is obtained from each AUV; and based on the local flow velocity measurements a global flow structure that is useful for plotting an optimal course for a vessel between two or more locations is predicted.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
       [0001]    This Application claims the benefit of U.S. Provisional Application 61/638,694 filed on Apr. 26, 2012, and incorporated herein by reference. 
     
    
     FIELD OF THE INVENTION 
       [0002]    The invention is directed to tracking coherent structures and manifolds on flows, and in particular to deploying autonomous underwater vehicles (AUVs) to track stable/unstable manifolds of general 2D conservative flows through local sensing alone. 
       BACKGROUND OF THE INVENTION 
       [0003]    In realistic ocean flows, time-dependent coherent structures, or Lagrangian coherent structures (LCS), are similar to separatrices that divide the flow into dynamically distinct regions. LCS are extensions of stable and unstable manifolds to general time-dependent flows (see, e.g., G. Haller and G. Yuan, “Lagrangian coherent structures and mixing in two-dimensional turbulence,” Phys. D, vol. 147, pp. 352-370 (December 2000)) and they carry a great deal of global information about the dynamics of the flows. 
         [0004]    For two-dimensional (2D) flows, LCS are analogous to ridges defined by local maximum instability, and can be quantified by local measures of Finite-Time Lyapunov Exponents (FTLE) (S. C. Shadden, F. Lekien, and J. Marsden, “Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows,” Physica D: Nonlinear Phenomena, vol. 212, no. 3-4, pp. 271-304 (2005)). Recently, LCS have been shown to coincide with optimal trajectories in the ocean which minimize the energy and the time needed to traverse from one point to another (see, e.g., T. Inane, S. Shadden, and J. Marsden, “Optimal trajectory generation in ocean flows,” in Proceedings of the 2005 American Control Conference, pp. 674-679, 2005, and C. Senatore and S. Ross, “Fuel-efficient navigation in complex flows,” in Proceedings of the 2008 American Control Conference, pp. 1244-1248, 2008). Furthermore, to improve weather and climate forecasting, and to better understand various physical, chemical, and geophysical processes in the ocean, there has been significant interest in the deployment of autonomous sensors to measure a variety of quantities of interest. One drawback to operating sensors in time-dependent and stochastic environments like the ocean is that the sensors will tend to escape from their monitoring region of interest. Since the LCS are inherently unstable and denote regions of the flow where more escape events may occur (see, e.g., E. Forgoston, L. Billings, P. Yecko, and I. B. Schwartz, “Set-based corral control in stochastic dynamical systems: Making almost invariant sets more invariant,” Chaos, vol. 21, 013116, 2011), knowledge of the LCS are of paramount importance in maintaining a sensor in a particular monitoring region. 
         [0005]    Existing work in cooperative boundary tracking for robotic teams that relies on one-dimensional (1D) parameterizations include C. Hsieh, Z. Jin, D. Marthaler, B. Nguyen, D. Tung, A. Bertozzi, and R. Murray, “Experimental validation of an algorithm for cooperative boundary tracking,” in Proceedings of the 2005 American Control Conference, pp. 1078-1083,2005, S. Susca, S. Martinez, and F. Bullo, “Monitoring environmental boundaries with a robotic sensor network,” IEEE Trans. on Control Systems Technology, vol. 16, no. 2, pp. 288-296, 2008, I. Triandaf and I. B. Schwartz, “A collective motion algorithm for tracking time-dependent boundaries,” Mathematics and Computers in Simulation, vol. 70, pp. 187-202, (2005) and V. M. Goncalves, L. C. A. Pimenta, C. A. Maia, B. Dutra, and G. A. S. Pereira, “Vector fields for robot navigation along time-varying curves in n-dimensions,” IEEE Trans. on Robotics, vol. 26, no. 4, pp. 647-659 (2010), for static and time-dependent cases respectively. Formation control strategies for distributed estimation of level surfaces and scalar fields in the ocean are presented in F. Zhang, D. M. Fratantoni, D. Paley, J. Lund, and N. E. Leonard, “Control of coordinated patterns for ocean sampling,” Int. Journal of Control, vol. 80, no. 7, pp. 1186-1199 (2007), K. M. Lynch, P. Schwartz, I. B. Yang, and R. A. Freeman, “Decentralized environmental modeling by mobile sensor networks,” IEEE Trans. on Robotics, vol. 24, no. 3, pp. 710-724 (2008), and W. Wu and F. Zhang, “Cooperative exploration of level surfaces of three dimensional scalar fields,” Automatica, the IFAC Journal, vol. 47, no. 9, pp. 2044-2051 (2011), and pattern formation for surveillance and monitoring by robot teams is discussed in J. Spletzer and R. Fierro, “Optimal positioning strategies for shape changes in robot teams,” in Proceedings of the IEEE Int. Conf. on Robotics &amp; Automation, Barcelona, Spain pp. 754-759, 2005, S. Kalantar and U. R. Zimmer, “Distributed shape control of homogeneous swarms of autonomous underwater vehicles,” Autonomous Robots (intl. journal), 2006, and M. A. Hsieh, S. Loizou, and V. Kumar, “Stabilization of multiple robots on stable orbits via local sensing,” in Proceedings of the Int. Conf. on Robotics &amp; Automation (ICRA), 2007). 
       BRIEF SUMMARY OF THE INVENTION 
       [0006]    According to the invention, a collaborative control method for tracking Lagrangian coherent structures (LCSs) and manifolds on flows employs at least three autonomous underwater vehicles (AUVs) each equipped with a local flow sensor. A first flow sensor is a tracking sensor and the other sensors are herding sensors for controlling and determining the actions of the tracking sensor. The AUVs are deployed in a body of water whereby the tracking sensor is positioned with respect to the herding sensors such that the herding sensors maintain a straddle formation across a boundary. A local flow velocity measurement is obtained from each AUV; and based on the local flow velocity measurements a global flow structure that is useful for plotting an optimal course for a vessel between two or more locations is predicted. 
         [0007]    The invention advantageously uses cooperative robots to find coherent structures without requiring a global picture of the ocean dynamics, and enables a team of robots to track the stable/unstable manifolds of general 2D conservative flows through local sensing alone. The invention provides tracking strategies for mapping LCS in the ocean using AUVs, using nonlinear dynamical and chaotic system analysis techniques to create a tracking strategy for a team of robots. The cooperative control strategy leverages the spatio-temporal sensing capabilities of a team of networked robots to track the boundaries separating the regions in phase space that support distinct dynamical behavior. Additionally, boundary tracking relies solely on local measurements of the velocity field. The method of the invention may be generally applied to any conservative flow. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0008]      FIG. 1  shows three robots tracking the stable structure B s  in a conservative vector field according to the invention; 
           [0009]      FIGS. 2A-B  show trajectories of 3 robots tracking a sinusoidal boundary ( FIG. 2A ) and a star-shaped boundary ( FIG. 2B ) according to the invention; 
           [0010]      FIGS. 3A-B  show trajectories of a 3 robot team tracking a star shape ( FIG. 3A ) and a snapshot of the multi-robot experiment ( FIG. 3B ) according to the invention; 
           [0011]      FIG. 4  shows a phase portrait of a time-independent double-gyre model according to the invention; and 
           [0012]      FIGS. 5A-5H  show snapshots of the trajectories of the team of 3 robots tracking Lagrangian coherent structures according to the invention. 
       
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
       [0013]    We consider the problem of controlling a team of N planar AUVs to collaboratively track the material lines that separate regions of flow with distinct fluid dynamics. This is similar to the problem of tracking the stable (and unstable) manifolds of a general nonlinear dynamical system where the manifolds separate regions in phase space with distinct dynamical behaviors. We assume the following 2D kinematic model for each of the AUVs: 
         [0000]        dx   i   /dt=V   i  cos θ i   +u   i    (1a)
 
         [0000]        dy   i   /dt=V   i  sin θ i   +v   i ;   (1b)
 
         [0000]    where x i =[x i ,y i ] T  is the vehicle&#39;s planar position, V i  and θ i  are the vehicle&#39;s linear speed and heading, and u i =[u i , v i ] T  is the velocity of the fluid current experienced/measured by the i th  vehicle. Additionally, we assume each agent can be circumscribed by a circle of radius r, i.e., each vehicle can be equivalently described as a disk of radius r. 
         [0014]    In this work, u i  is provided by a 2D planar conservative vector field described by a differential equation of the form 
         [0000]        dx/dt=F ( x ).   (2)
 
         [0000]    In essence, u i =F x (x i ) and v i =F y (x i ). Let B S  and B U  denote the stable and unstable manifolds of Eq. (2). In general, B S  and B U  are the separating boundaries between regions in phase space with distinct dynamics. For 2D flows, B* are simply one-dimensional curves where * denotes either stable (S) or unstable (U) boundaries. For a small region centered about a point on B*, the system is unstable in one dimension. Finally, let ρ(B*) denote the radius of curvature of B* and assume that the minimum of the radius of curvature ρ min (B*)&gt;r. This last assumption is needed to ensure the robots do not lose track of the B* due to sharp turns. 
         [0015]    The objective is to develop a collaborative strategy to enable a team of robots to track B* in general 2D planar conservative flow fields through local sampling of the velocity field. While the focus is on the development of a tracking strategy for B S , the method can be easily extended to track B U  since B U  are simply stable manifolds of Eq. (2) for t&lt;0. 
         [0016]    The PIM Triple Procedure 
         [0017]    The method of the invention originates from the Proper Interior Maximum (PIM) Triple Procedure, H. E. Nusse and J. A. Yorke, “A procedure for finding numerical trajectories on chaotic saddles,” Physica D Nonlinear Phenomena, vol. 36, pp. 137-156, 1989 (hereinafter “Nusse et al.”)—a numerical technique designed to find stationary trajectories in chaotic regions with no attractors. While the original procedure was developed for chaotic dynamical systems, the approach can be employed to reveal the stable set of a saddle point of a general nonlinear dynamical system. The procedure consists of iteratively finding an appropriate PIM Triple on a saddle straddling line segment and propagating the triple forward in time. 
         [0018]    Given the dynamical system described by Eq. (2), let D ∈ R 2  be a closed and bounded set such that D does not contain any attractors of Eq. (2). Given a point x ∈ D, the escape time of x, denoted by T E (x), is the time x takes to leave the region D under the differential map given by Eq. (2). 
         [0019]    Let J be a line segment that crosses the stable set B S  in D, i.e., the endpoints of the J are on opposite sides of B S . Let {x L ,x C ,x R } denote a set of three points in J such that x C  denotes the interior point. Then {x L ,x C ,x R } is an Interior Maximum triple if T E (x C )&gt;max{T E (x L ),T E (x R )}. Furthermore, {x L ,x C x R } is a Proper Interior Maximum (PIM) triple if it is an interior maximum triple and the interval [X L , x R ] in J is a proper subset of J. Then the numerical computation of any PIM triple can be obtained iteratively starting with an initial saddle straddle line segment J 0 , Let x t0  and X R0  denote the endpoints of J 0  and apply an ε 0 &gt;0 discretization of J 0  such that x L0 =q 0 &lt;q 1 &lt; . . . &lt;q M =x R0 . For every point qi, determine T E (q i ) by propagating q i  forward in time using Eq. (2). Then the PIM triple in J 0  is given by the points {q k−1 ,q k ,q k+1 } where q k =argmax i=1, . . . , M T E (q i ). This PIM triple can then be further refined by choosing J 1  to be the line segment containing {q k−1 ,q k ,q k+1 } and reapplying the procedure with another ε 1 &gt;0 discretization where ε 1 &lt;ε 0 . 
         [0020]    Given an initial saddle straddling line segment J 0 , it has been shown that the line segment given by any subsequent PIM triple on J 0  is also a saddle straddling line segment [H. E. Nusse and J. A. Yorke, “A procedure for finding numerical trajectories on chaotic saddles,” Physica D Nonlinear Phenomena, vol. 36, pp. 137-156, 1989.]. Furthermore, if we use a PIM triple x(t)={x L ,x C ,x R } as the initial conditions for the dynamical system given by Eq. (2) and propagate the system forward in time by Δt, then the line segment containing the set x(t+Δt), J t+Δt , remains a saddle straddle line segment. As such, the same numerical procedure can be employed to determine an appropriate PIM triple on J t+≢t . This procedure can be repeated to eventually reveal the entire stable set B S  and unstable set B U  within D if time was propagated forwards and backwards respectively. Furthermore, since the procedure always begins with a valid saddle straddling line segment, by construction, the procedure always results in a non-empty set. 
         [0021]    Building upon the PIM Triple Procedure, as described below the invention utilizes a cooperative saddle straddle control strategy for a team of N≧3 robots to track the stable (and unstable) manifolds of a general conservative time-independent flow field F(x). The invention differs from the PIM procedure where it relies solely on information gathered via local sensing and shared through the network. In contrast, a straight implementation of the PIM Triple Procedure necessitates global knowledge of the structure of the system dynamics throughout a given region given its reliance on computing escape times. 
         [0022]    Controller Synthesis 
         [0023]    Consider a team of three robots and identify them as robots {L,C,R}. While the robots may be equipped with similar sensing and actuation capabilities, we propose a heterogeneous cooperative control strategy. 
         [0024]    Let x(0)=[x L   T (0),x C   T (0),x R   T (0)] T  be the initial conditions for the three robots. Assume that x(0) lies on the line segment J 0  where J 0  is a saddle straddle line segment and {x L (0),x C (0),x R (0)} constitutes a PIM triple. Similar to the PIM Triple Procedure, the objective is to enable the robots to maintain a formation such that a valid saddle straddle line segment can be maintained between robots L and R. Instead of computing the escape times for points on J 0  as proposed by the PIM Triple Procedure, robot C must remain close to B S  using only local measurements of the velocity field provided by the rest of the team. As such, we refer to robot C as the tracker of the team while robots L and R maintains a straddle formation across the boundary at all times. Robots L and R may be thought of herding robots, since they control and determine the actions of the tracking robot. 
       Straddling Formation Control 
       [0025]    The controller for the straddling robots consists of two discrete states: a passive control state, U P , and an active control state, U A . The robots initialize in the passive state U P  where the objective is to follow the flow of the ambient vector field. Therefore, V i =0 for i=L,R. Robots execute U P  until they reach the maximum allowable separation distance d Max  from robot C. When ∥x i −x C ∥&gt;d max  robot i switches to the active control state, U A , where the objective is to navigate to a point p i  on the current projected saddle straddle line segment Ĵ i  such that ,∥p i −p C |=d Min  and p C  denotes the midpoint of Ĵ i . When robots execute U A , V i =∥p i −x i −u i ∥ and θ i (t)=α i (t) where α i  is the angle between the desired, (p i −x i ), and current heading, u i , of robot i as shown in  FIG. 1 . In summary, the straddling control strategy for robots L and R is given by 
         [0000]    
       
         
           
             
               
                 
                   
                     V 
                     i 
                   
                   = 
                   
                     { 
                     
                       
                         
                           0 
                         
                         
                           
                             
                               if 
                                
                               
                                   
                               
                                
                               
                                 d 
                                 Min 
                               
                             
                             &lt; 
                             
                                
                               
                                 
                                   x 
                                   i 
                                 
                                 - 
                                 
                                   x 
                                   C 
                                 
                               
                                
                             
                             &lt; 
                             
                               d 
                               Max 
                             
                           
                         
                       
                       
                         
                           
                              
                             
                               
                                 ( 
                                 
                                   
                                     p 
                                     i 
                                   
                                   - 
                                   
                                     x 
                                     i 
                                   
                                 
                                 ) 
                               
                               - 
                               
                                 u 
                                 i 
                               
                             
                              
                           
                         
                         
                           
                             otherwise 
                             , 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     3 
                      
                     a 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     θ 
                     i 
                   
                   = 
                   
                     { 
                     
                       
                         
                           0 
                         
                         
                           
                             
                               if 
                                
                               
                                   
                               
                                
                               
                                 d 
                                 Min 
                               
                             
                             &lt; 
                             
                                
                               
                                 
                                   x 
                                   i 
                                 
                                 - 
                                 
                                   x 
                                   C 
                                 
                               
                                
                             
                             &lt; 
                             
                               d 
                               Max 
                             
                           
                         
                       
                       
                         
                           
                             α 
                             i 
                           
                         
                         
                           
                             otherwise 
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     3 
                      
                     b 
                   
                   ) 
                 
               
             
           
         
       
     
         [0000]    We note that while the primary control objective for robots L and R is to maintain a straddle formation across B S , robots L and R are also constantly sampling the velocity of the local vector field and communicating these measurements and their relative positions to robot C. Robot C is then tasked to use these measurements to track the position of B S . 
         [0026]    Manifold Tracking Control 
         [0027]    Let û L (t), û C , and û C (t) denote the current velocity measurements obtained by robots L, C, and R at their respective positions. Let d(•,•) denote the Euclidean distance function and assume that d(x C ,B S )&lt;ε such that ε&gt;0 is small. Given the straddle line segment J t  such that x L (k) and x R (k) are the endpoints J t , we consider an ε t &lt;ε discretization of J t  such that x L =q 1 &lt;q 2 &lt; . . . &lt;q M =x R . The objective is to use the velocity measurements provided by the team to interpolate the vector field at the points q 1 , . . . , q M . Since Eq. (2) has C l  continuity and if x C  is ε-close to B S , then the point q B =argmax k=1, . . . , M u(q k ) T û C (t) should be δ-close to B S  where ε&lt;δ&lt;A and A is a small enough positive constant. 
         [0028]    While there are numerous vector field interpolation techniques available (J. C. Agui and J. Jimenez, “On the performance of particle tracking,” Journal of Fluid Mechanics, vol. 185, pp. 447-468, 1987, and E. J. Fuselier and G. B. Wright, “Stability and error estimates for vector field interpolation and decomposition on the sphere with rbfs,” SIAM J. Numer. Anal., vol. 47, pp. 3213-3239, 2009), we employ the inverse distance weighting method. For a given set of velocity measurements û i (t) and corresponding position estimates {circumflex over (x)} i (t), the velocity vector at some point q k  is given by 
         [0000]    
       
         
           
             
               u 
                
               
                 ( 
                 
                   q 
                   k 
                 
                 ) 
               
             
             = 
             
               
                 ∑ 
                 j 
                 
                     
                 
               
                
               
                   
               
                
               
                 
                   ∑ 
                   
                     i 
                     = 
                     1 
                   
                   N 
                 
                  
                 
                     
                 
                  
                 
                   
                     
                       w 
                       ij 
                     
                      
                     
                       
                         
                           u 
                           ^ 
                         
                         i 
                       
                        
                       
                         ( 
                         j 
                         ) 
                       
                     
                   
                   
                     
                       ∑ 
                       j 
                       
                           
                       
                     
                      
                     
                         
                     
                      
                     
                       
                         ∑ 
                         
                           i 
                           = 
                           1 
                         
                         N 
                       
                        
                       
                         w 
                         ij 
                       
                     
                   
                 
               
             
           
         
       
     
         [0000]    where w ij =∥{circumflex over (x)} i (j)−q i ∥ −2 . Rather than rely solely on the current measurements provided by the three robots, it is possible to include the recent history of û i (t) to improve the estimate of u(q k ), i.e., û i (t−ΔT), û i (t−2ΔT), and so on, where ΔT is the sampling period and i={L,C,R}. Thus, the control strategy for the tracking robot C is given by 
         [0000]        V   C =∥[( q   B   +bû   B )− x   C   ]−u   C ∥  (4a)
 
         [0000]      θ C =β C    (4b)
 
         [0000]    where β C  denotes the difference in the heading of robot C and the vector (q B −û B ) and b&gt;r is a small number. The term bû B  is included to ensure that the control strategy aims for a point in front of robot C rather than behind it. As such, the projected saddle straddle line segment Ĵ t  at each time step is given by p c =q C +bu C  with Ĵ t  orthogonal to B S  at q C  and ∥Ĵ t ∥ chosen to be in the interval [2d Min ,2d Max ]. 
         [0029]    Analysis 
         [0030]    Regarding the implementation of the saddle straddle control strategy, we begin with the following key assumption on the robots&#39; initial positions. 
         [0031]    Assumption 1 Given a team of three robots {L, C, R}, assume that d(x C ,(0),B S )&lt;ε for a small value of ε&gt;0, ∥x L −x C ∥=∥|x R −x∥=d Min  with d Min &gt;2r, and robots L and R are on opposite sides of B S . 
         [0032]    In other words, assume that the robots initialize in a valid PIM triple formation and their positions form a saddle straddle line segment orthogonal to B S . Our main result concerns the validity of the saddle straddle control strategy. 
         [0033]    Theorem 1 Given a team of 3 robots with kinematics given by Eq. (1) and u i  given by Eq. (2), the feedback control strategy Eq. (3) and Eq. (4) maintains a valid saddle straddle line segment in the time interval [t,t+Δt] if the initial positions of the robots, x(t), is a valid PIM triple. 
         [0034]    The above theorem guarantees that for any given time interval [t,t+Δt] the team maintains a valid PIM triple formation. As such, the iterative application of the proposed control strategy leads to the following proposition. 
         [0035]    Proposition 1 Given a team of 3 robots with kinematics given by Eq. (1) and u i  given by Eq. (2), the feedback control strategy results in an estimate of B S , denoted as {circumflex over (B)} S , such that &lt;B S , {circumflex over (B)} S &gt; L2 &lt;W for some W&gt;0 where &lt;•,•&gt; L2  denotes the inner product (which provides an L 2  measure between the B S  and {circumflex over (B)} S  curves). 
         [0036]    From Theorem 1, since the team is able to maintain a valid PIM triple formation across B S  for any given time interval [t,t+Δt], this ensures that an estimate of B S  in the given time interval also exists. Applying this reasoning in a recursive fashion, one can show that an estimate of B S  can be obtained for any arbitrary time interval. Preferably, one also determines the bound on W such that {circumflex over (B)} S  results in a good enough approximation since W depends on the sensor and actuation noise, the vector interpolation routine, the sampling frequency, and the time scales of the flow dynamics. 
         [0037]    Results 
         [0038]    Simulations: 
         [0039]    We illustrate the proposed control strategy given by Eq. (3) and Eq. (4) with the following simulation results.  FIG. 2A  shows the trajectories of three robots tracking a sinusoidal boundary while  FIG. 2B  shows the team tracking a 1D star-shaped boundary. We note that throughout the entire length of the simulation, the team maintains a saddle straddle formation across the boundary. 
         [0040]    In both examples, u=−a∇φ−b∇×ψ where a,b&gt;0 and •(x) is an artificial potential function such that Φ(x)=0 for all x ∈ B• and Φ(x)&lt;0 for any x ∈ R 2 /B•. The vector ψ is a 3×1 vector whose entries are given by [0,0,γ(x, y)] T  where γ(x,y) is the curve describing the desired boundary. Lastly, the estimated position of the boundary is given by the position of the tracking robot, i.e., robot C. In these examples, we filtered the boundary position using a simple first-order low pass filter. 
         [0041]    Experiments 
         [0042]    We also implemented the control strategy on our multi-robot testbed. The testbed consisted of three mSRV-1 robots in a 4.8×5.4 meter workspace. The mSRV-1 are differential-drive robots equipped with an embedded processor, color camera, and 802.11 wireless capability. Localization for each robot was provided via a network of overhead cameras.  FIG. 3A  shows the trajectories of the robots tracking a star shaped boundary.  FIG. 3B  is a snapshot of the experimental run. 
         [0043]    Extension to Periodic Boundaries 
         [0044]    Next, we consider the system of 3 robots with kinematics given by Eq.(1) where u i  is determined by the wind-driven double-gyre flow model with noise 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       x 
                       . 
                     
                     = 
                     
                       
                         
                           - 
                           π 
                         
                          
                         
                             
                         
                          
                         A 
                          
                         
                             
                         
                          
                         
                           sin 
                            
                           
                             ( 
                             
                               π 
                                
                               
                                 
                                   f 
                                    
                                   
                                     ( 
                                     
                                       x 
                                       , 
                                       t 
                                     
                                     ) 
                                   
                                 
                                 s 
                               
                             
                             ) 
                           
                         
                          
                         
                           cos 
                            
                           
                             ( 
                             
                               π 
                                
                               
                                 y 
                                 s 
                               
                             
                             ) 
                           
                         
                       
                       - 
                       
                         μ 
                          
                         
                             
                         
                          
                         x 
                       
                       + 
                       
                         
                           η 
                           1 
                         
                          
                         
                           ( 
                           t 
                           ) 
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   
                     5 
                      
                     a 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     
                       y 
                       . 
                     
                     = 
                     
                       
                         π 
                          
                         
                             
                         
                          
                         A 
                          
                         
                             
                         
                          
                         
                           cos 
                            
                           
                             ( 
                             
                               π 
                                
                               
                                 
                                   f 
                                    
                                   
                                     ( 
                                     
                                       x 
                                       , 
                                       t 
                                     
                                     ) 
                                   
                                 
                                 s 
                               
                             
                             ) 
                           
                         
                          
                         
                           sin 
                            
                           
                             ( 
                             
                               π 
                                
                               
                                 y 
                                 s 
                               
                             
                             ) 
                           
                         
                          
                         
                           
                              
                             f 
                           
                           
                              
                             x 
                           
                         
                       
                       - 
                       
                         μ 
                          
                         
                             
                         
                          
                         y 
                       
                       + 
                       
                         
                           η 
                           2 
                         
                          
                         
                           ( 
                           t 
                           ) 
                         
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   
                     5 
                      
                     b 
                   
                   ) 
                 
               
             
             
               
                 
                   
                     f 
                      
                     
                       ( 
                       
                         x 
                         , 
                         t 
                       
                       ) 
                     
                   
                   = 
                   
                     
                       ɛ 
                        
                       
                           
                       
                        
                       
                         sin 
                          
                         
                           ( 
                           
                             
                               ω 
                                
                               
                                   
                               
                                
                               t 
                             
                             + 
                             ψ 
                           
                           ) 
                         
                       
                        
                       
                         x 
                         2 
                       
                     
                     + 
                     
                       
                         ( 
                         
                           1 
                           - 
                           
                             2 
                              
                             ɛ 
                              
                             
                                 
                             
                              
                             
                               sin 
                                
                               
                                 ( 
                                 
                                   
                                     ω 
                                      
                                     
                                         
                                     
                                      
                                     t 
                                   
                                   + 
                                   ψ 
                                 
                                 ) 
                               
                             
                           
                         
                         ) 
                       
                        
                       
                         x 
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   
                     5 
                      
                     c 
                   
                   ) 
                 
               
             
           
         
       
     
         [0045]    When ε=0, the double-gyre flow is time-independent, while for ε≠0, the gyres undergo a periodic expansion and contraction in the x direction. In Eq. (5a-c), A approximately determines the amplitude of the velocity vectors, ω/2π gives the oscillation frequency, ε determines the amplitude of the left-right motion of the separatrix between the gyres, ψ is the phase, μ determines the dissipation, s scales the dimensions of the workspace, and η i (t) describes a stochastic white noise with mean zero and standard deviation σ=√{square root over (2I)}, for for noise intensity I. In this work, η i (t) can be viewed as either measurement or environmental noise.  FIG. 4  shows the phase portrait of the time-independent double-gyre model. 
         [0046]      FIGS. 5A-5H  show trajectories of the team of 3 robots tracking Lagrangian coherent structures of the system described by Eq. (5a-c) with A=10, μ=0.005, ε=0.1, ψ=0, I=0.01, and s=50. The trajectories of the straddling robots are shown in black and the estimated LCS is shown in white. 
         [0047]    While the present invention has been described with respect to exemplary embodiments thereof, it will be understood by those of ordinary skill in the art that variations and modifications can be effected within the scope and spirit of the invention.