University of Southern California Dissertations and Theses

Physics-based and data-driven models for bio-inspired flow sensing and motion planning

(USC Thesis Other)

Physics-based and data-driven models for bio-inspired flow sensing and motion planning

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Physics-based and data-driven models for bio-inspired ow sensing and motion planning A Dissertation Presented to the Faculty of the USC Graduate School in Partial Fulllment of the Requirements of the Degree Doctor of Philosophy Aerospace Engineering Brendan Colvert Advisor: Dr. Eva Kanso Department of Aerospace and Mechanical Engineering University of Southern California August 2018 Abstract From naval applications to space exploration, there is a growing need to develop engineering solu- tions for autonomous submersibles. If subsurface vehicles are going to be capable of high levels of autonomy, they will need the ability to sense and decode information embedded in the uid ows that surround them and make decisions based on the cues they extract. Inspired by the capabilities of many aquatic creatures that leverage various ow sensing modalities to achieve critical tasks such as navigation, mating, and predation, we use both physics-based and data-driven techniques to model how biological and bioinspired sensory systems can extract and react to relevant cues in the hydrodynamic environment. With a physics-based framework, we address the problem of optimizing sensory layouts, develop a local ow characterization scheme, and test hypotheses about behavioral responses to ow stimuli. We then use data-driven techniques to identify and decipher which cues in the hydrodynamic environment are most relevant for achieving classication of un- known objects in the ow eld. These results serve as a framework for addressing more complex and pressing challenges such as detection and localization of unknown objects in the far eld and autonomous navigation. i To Gavin, Danielle, Gabrielle, and Jamie. ii Acknowledgements The successful completion of this work would not have been possible without the aid of many people over the course of my life. Although it is impossible to thank every person who has helped me along the way, I would like to take this opportunity to name a few individuals whose guidance and support has been particularly impactful. My academic success has been in large part due to the great number of excellent teachers who have shaped my mind. All of them delivered a high quality of instruction and challenged me to push myself academically. However, in particular, I regard the eorts Richard Quaglieri as having a uniquely transformative impact on my academic work ethic. `Mr. Q.', as we called him during my seventh and eighth grade math classes, taught me that success does not come from simply accomplishing what you've been asked to do, but instead pushing yourself to work harder and gain a deep understanding of the material that you are tasked with learning. Without Mr. Q., I am not sure whether I would have developed the work ethic that has allowed me to complete this dissertation. I would also like to thank my undergraduate research advisor, Dr. Veronica Eliasson. Her passion for research as well as her innovative and energetic presentation of its merits in the under- graduate classroom were responsible for inspiring my own interest in academic research. Veronica oered me the opportunity during my time as an undergraduate to work in her lab and to gain experience. Her passion for her work ignited my desire to pursue a Ph.D. and I am certain that her inspiration played a central role in my success today. The completion of my Ph.D. would not have been possible had I not had the guidance and support of my graduate advisor, Dr. Eva Kanso. From the rst time we met four years ago, I knew that Eva was the person who would nurture and guide my academic mind to this moment. Eva showed me the beauty of the interplay between mathematics, engineering, and biology. Her adherence to a singular standard of excellence when it comes to communication and presentation of scientic knowledge is both an inspiration to me and a habit that will serve me well for the rest of my professional life. Eva repeatedly gave me opportunities to grow and mature as an academic, even quite early on in my tenure as a graduate student, and these chances to show my knowledge and skills allowed me to become the researcher that I am today. I will be forever grateful for the mentorship and friendship of Eva that shaped my time as a graduate student at USC. A word of credit is also due to those with whom I have collaborated with over the years including Nicholas Foster, Tom Peaslee, Monica Nguyen, Xingtian Tao, Mohamad Alsalman, Dr. Kevin Chen, Dr. Geo Spedding, Dr. Geng Liu, and Dr. Haibo Dong. I would also like to thank the members of the Kanso Biodynamical Systems Laboratory, whose insight and companionship have proved invaluable over the past four years. iii iv The greatest credit is due to my family for raising and supporting me as I grew into the person that I am today. I am forever indebted to my parents for the hardship they endured to provide me with the highest quality of education. My family is my bedrock and without them, I would not have accomplished anything. I am also thankful for the support and companionship of my anc ee, Jamie, during the past ve years. Many hours of studying, working, and traveling have been made bearable with the knowledge that I am working to build a life for the two of us and for our future family. Funding Sources The work of Brendan Colvert and Eva Kanso is partially supported by the Of- ce of Naval Research (ONR) grants N00014-14-1-0421 and N00014-17-1-2287 (to E.K.). Brendan Colvert also acknowledges support from the Department of Defense (DoD) through the National Defense Science & Engineering Graduate Fellowship (NDSEG) Program. The physics-based ow characterization work was, in part, carried out in collaboration with Kevin K. Chen, who was sup- ported by the Viterbi Postdoctoral Fellowship provided by the Viterbi School of Engineering at the University of Southern California. The data-driven work was, in part, carried out in collaboration with Mohamad Alsalman, whose work is partially supported by the Army Research (ARO) grant W911NF-16-1-0074 (to E.K.). Contents 1 Introduction and Motivation 1 2 Literature Review 9 3 Flow Sensing: Physics-Based Approach 13 3.1 Bio-inspired ow sensing system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2 Measuring shear ows using velocity sensors . . . . . . . . . . . . . . . . . . . . . . . 16 3.3 Characterizing ow types using velocity sensors . . . . . . . . . . . . . . . . . . . . . 22 3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4 Flow Sensing: Data-Driven Approach 32 4.1 Vortex wakes of pitching airfoils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Neural networks for wake classication . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3 Network performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5 Flow Sensing in Fish-Inspired Behavior 47 5.1 Potential ow model for pressure sensing . . . . . . . . . . . . . . . . . . . . . . . . . 48 5.2 Optimal sensor layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.3 Response to perceived pressure dierence . . . . . . . . . . . . . . . . . . . . . . . . 52 5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 6 Physics-Based Source Tracking Using Flow Sensing 56 6.1 Sensory system and signal eld model . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.2 Analytical performance of sensor in nding source . . . . . . . . . . . . . . . . . . . 61 6.3 Locating a pitching airfoil by following its wake . . . . . . . . . . . . . . . . . . . . . 69 6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 7 Data-Driven Source Tracking Using Visual Sensing 74 7.1 Mobile agent with visual sensing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 7.2 Data-driven controller design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 7.3 Performance of trained agent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 7.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 8 Conclusions and Outlook 86 Bibliography 88 v List of Figures 1.1 Von K arm an vortex street in the wake behind Heard Island in the southern Indian Ocean. This image is a modied version of the original captured by the Moderate Resolution Imaging Spectroradiometer (MODIS) and is in the public domain courtesy of NASA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Illustration of the dierence between local (surface pressure or velocity) and global (entire ow eld) sensing. An airfoil in uniform ow at (left) 0 angle of attack, (center) 15 angle of attack, and (right) 15 angle of attack with a vortex on the leading edge. We plot the (top) pressure distribution and (middle) velocity distri- bution on the surface of the airfoil. (bottom) Then, we plot the entire pressure eld and the associated streamlines. Can global ow characteristics be determined from measuring only local ow properties? . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Schematic block diagram for autonomous submersible. Sensing { How does the agent extract information from the hydrodynamic environment? Signal Processing { How does the agent extract relevant cues from measured data? Actuation { What actions does the agent take in response to these cues? Swimming Dynamics { How do these actions translate into motion? Hydrodynamic Environment { What information is encoded in the ow eld by the movement of the platform as well as other agents? How is this information measured by the sensors? . . . . . . . . . . . . . . . . . . . 6 1.4 (a) Systems of interest have similar forms { an input x is transformed into an output y by the mapping y = f(x). (b) The job of the modeler is to approximate this map- ping in one of two ways: Physics-Based Modeling uses physical intuition to construct a model f phys (x) which estimates the output ^ y from input x. Data-Driven Model- ing constructs a model f phys (x;), then uses known input-output pairs to optimize parameters , which allows estimation of the output ^ y from input x. . . . . . . . . 7 3.1 Problem setup and summary of approach. (a) The goal is to develop a framework to detect select ow features by a local sensory array. (b) The uid velocity eld is linearized around the sensory array. (c) A simple bio-inspired sensory system consists of two sensors and one sensory measurement, the dierence between the two sensors. (d) The sensory measurements are decoded to detect the desired ow properties. . . 14 3.2 The Taylor decomposition of (a) a full ow eld into (b) a constant component, (c) a linearly varying component and (d) higher order terms, where the velocity magnitude is shown in color, with blue to red representing low to high speed. (e) Constant and linear components are combined into linearized ow for comparison with full nonlinear ow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 vi LIST OF FIGURES vii 3.3 Two-sensor system with sensory array extent l. Velocity is measured in body frame (b 1 ; b 2 ). Orientation with respect to inertial frame is parameterized by the angle . The relative rotation between angles n and m is given by m;n = n m . . . . 16 3.4 (a) Sensors that measure s 1 = u 1 (x l ) u 1 (x r ) are called orthogonal sensors and (b) sensors that measure s 2 = u 2 (x l ) u 2 (x r ) are called parallel sensors. (c) Full velocity sensors have access to both u 1 and u 2 . By denition, u 1 = ub 1 and u 2 = ub 2 . 17 3.5 (a) Symmetries in the orthogonal measurement. Measurements at = o , o , o , and + o all produce the same sensory measurement. (b) This is due to the re ective symmetries in the measurement, with axes of re ection at = 0 and (horizontal axis) as well as ==2 and 3=2 (vertical axis). . . . . . . . . . . . . . 18 3.6 (a) Symmetries in the parallel measurement. Measurements at = o , + o , =2 o , and=2 o all produce the same sensory measurement. (b) This is due to the re ective symmetries in the measurement, with axes of re ection at ==4 and3=4 as well as = 3=4 and =4 (both diagonal axes). . . . . . . . . . . . . . 19 3.7 Measured value over true value of for orthogonal sensors. (a) `First' branch on the left for > 0 and (b) `second' branch on the right, also for > 0. If < 0, then the two branches switch. Black lines delineate contours where = 0, which is where the two branches of the solution exchange forms. . . . . . . . . . . . . . . . . . . . . . . 20 3.8 Measured value over true value of for parallel sensors. Solutions exchange sign at black circles, when n1;n ==2 mod =2. . . . . . . . . . . . . . . . . . . . . . . . 21 3.9 Variation of ow type with parameter for (a) counter-clockwise rotation, (b) counter- clockwise shear, (c) extensional ow, (d) clockwise shear, and (e) clockwise rotation. 23 3.10 (a{d) Lamb{Oseen vortex, (e{h) Gaussian shear and (i{l) cylinder in uniform ow at Re D = 60. (a,e,i) Velocity vector eld u and vorticity eld !. (b,f,j ) Flow type, with the black contours depicting =1 (c,g,k) Flow intensity , (d,h,l) Flow orientation r 1 vector eld, with the colormap showing the orientation angle =. The cylinder data are taken from [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.11 (a,c) Flow type and (b,d) ow intensity measurements for varying sensor pa- rameters in nondimensionalized (a,b) Gaussian shear prole and (c,d) Lamb-Oseen vortex. Analytical solutions shown in dashed black lines. First plots show xed sensor length l = 1 and varying =f0:1; 0:3; 0:5g. Second plots show xed = 0:1 and varying l =f0:5; 1; 5g. . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.12 Velocity gradients and ow parameters as function of time for sensory array placed at x o = 9e 1 + 0e 2 in the cylinder ow, as shown in gure 3.10. True values for velocity gradients and parameters are computed via nite dierence and shown as dashed black lines. Computed values using M + full are shown for (a) , (b) , (c) !, (d) , and (e) . Sensor system parameters are = 0:05 and l =f0:1; 0:5; 1g. The data are adapted from [12]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.1 Classication of ow patterns behind an oscillating airfoil: (left) proposed approach using only local measurements, and (right) classic approach relying on global inspec- tion of the spatiotemporal patterns in the ow eld. . . . . . . . . . . . . . . . . . . 33 LIST OF FIGURES viii 4.2 Wake behind an oscillating airfoil for A = 1:656 and St ranging from 0:038 to 0:079 as shown in (a). Three wake types are observed, 2P+4S, 2P+2S, and 2S, shown respectively by plotting contours of vorticity eld !(x;t) for (b) St = 0:044, (c) St = 0:064 and (d) St = 0:075. Time traces of vorticity at downstream locations I, II, and III are plotted in gure 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.3 Time traces of vorticity in the wake of an oscillating airfoil for wake types 2P+4S, 2P+2S, and 2S. Point-measured time traces are taken at several locations down- stream along the lines indicated by I, II and III in Figure 4.2. . . . . . . . . . . . . 36 4.4 Architecture of the neural network consisting of four successive layers: a linear layer, a nonlinear activation layer, another linear layer, and a nonlinear normalization layer. 37 4.5 A neural network is trained by optimizing a loss function over the network parameters subject to a training set of snapshots s m and labels y m . The optimization problem is solved using stochastic gradient descent. . . . . . . . . . . . . . . . . . . . . . . . 38 4.6 Error incurred by ve dierent trained networks in predicting the wake type or class of the test simulations. The error is plotted as a function of the number of features N f for (a) all classes, as well as for the individual classes (b) 2P+4S, (c) 2P+2S, and (d) 2S. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.7 Accuracy of the trained network in predicting the ow type of the test simulations. Diagonal entries show percentages of correct classication and o-diagonal entries show percentages of wrong classication and to which class they were erroneously assigned. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 4.8 Samples of the time traces of vorticity obtained after altering the data by (a) adding white noise, (b) truncating the total time interval of measurement below the oscil- lation period of the airfoil, and (c) decreasing the sampling rate of the sensor below the computational rate 1=t, where t is the simulation time step. . . . . . . . . . . 41 4.9 Accuracy of the trained network in predicting the ow type of the test simulations with (a) noisy measurements, (b) reduced measurement intervals, and (c) reduced sampling rate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.10 (a) Salience matrix of the N f = 34 features of the model for N c = 3 classes. (b) The ve most positively and negatively salient features for each class are plotted as functions of time and ranked by salience. Features that occur as a pair of positive and negative identiers for dierent classes are highlighted using the same color. These features are also marked by colored boxes above their corresponding columns in the salience matrix. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.11 Comparison of the vorticity eld and the performance of the neural network in cor- rectly classifying the wake as a function of the positionx m where the vorticity eld is probed. (left) Contours of the vorticity eld !(x;t) repeated from Figure 4.2 to facilitate comparison and (right) the spatial distribution of the likelihood ^ y of correct classication for wake types (a) 2P+4S atSt = 0:042, (b) 2P+2S atSt = 0:066, and (c) 2S at St = 0:073. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.1 (a) Spatial distribution of the lateral-line canal neuromasts in Astyanax fasciatus (Blind Mexican cave sh), inspired from [74]. (b) Canal density in various species of sh, adapted from [88]. (c) Lateral-line canal seen from the top and corresponding model system proposed here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 LIST OF FIGURES ix 5.2 Schematic of the model system (a) in the physical plane, (b) in a body-xed frame, with dashed lines corresponding to level sets of constant elliptic coordinates; and (c) conformally mapped to the circle plane, with dashed lines corresponding to level sets of constant polar coordinates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3 (a) Pressure distributions along the body for b=a = 0:266, = 5 , V = 0. (b) Pressure dierence p in (5.5) as a function of for = 5 ; 10 ; and 20 andV = 0. (c) Stimulation (p=) for = 5 ; 10 ; and 20 and its theoretical limit obtained by using small in (5.5) compared with the average canal density of Figure 5.1(b). (d) Pressure dierence as as a function of for V = 5; 10; 20; 40 and = 5 . (e) p=V for same parameter values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 5.4 (a) Pressure dierence p as a function of body orientation , swimming speed V and sensor location . Black contour lines show p = 0. (b) p = 0 are xed points of the governing equations for which the body does not change orientation. The stability of these xed directions depend on the swimming speed V . The ellipse aspect ratio is set to b=a = 0:3, but the results hold qualitatively for all 0<b=a< 1. 51 5.5 Trajectories of swimmer in uniform ow for (a) V = 1:44 and (b) V = 11:5. Initial orientations are (0) = =6;=3;=2; 2=3; 5=6. Parameter values are b=a = 0:3, V cr = 4:33, and G = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 5.6 Basins of attraction: (a) For V = 0, the set of all initial orientations are equally split into a white subset that tends to align in the direction opposite to the incoming ow, and a grey subset that tends to align in the same direction as the incoming ow. (b) For 0 < V < V cr , the grey region shrinks until at V = V cr , a pitchfork bifurcation occurs making = unstable. (c) For V >V cr , all trajectories, except those parallel to , tend to align in the direction opposite to the ow. . . . . . . . . 53 6.1 (a) A mobile sensor, located at x s and equipped with body-xed frame (b 1 ; b 2 ), moves at constant speed V and changes its orientation with angular velocity in response to local sensory measurements of the signal eld f(x;t). Near the sensor, f(x;t) is expanded into a series of traveling waves with modal amplitudes a n and b n , frequencies w n , wavenumberskk n k and directions k n =kk n k indicating the local direction of signal propagation. (b) Relation between the inertial frame (e 1 ; e 2 ), body-xed frame (b 1 ; b 2 ), where b 1 points in the heading direction of the sensor, and body-xed frame (c 1 ; c 2 ), where c 1 is in the direction pointing from the sensor to the source. The pointing or alignment error n is the angle between c 1 and r n =kr n k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6.2 Block diagram depicting the sensory and control laws of the mobile sensor. The time history of the signal eldf(x s ;t) is measured locally at x s over one oscillation period of the signal and transformed to the frequency domain ^ f(x s ;w n ). The magnitude m n (x) and phase elds n (x) are computed from the spectral information to obtain the sensory output s and determine the gain function G of the feedback controller =Gs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.3 Radially-symmetric signal eld in (6.19) for ` = 6:5. (a) Snapshot of the signal eld f(x;t) at time t = 0, (b) the magnitude eld m(x), (c) the phase eld (x), (d) the direction eldr. Note that the associated alignment error (x) is identically zero. 62 LIST OF FIGURES x 6.4 Behavior of the mobile sensor in the radially-symmetric signal eld in (6.19). The left panel depicts trajectories in the (x;y) plane starting at positions located successively further away from the source. In the middle panel, trajectories are shown in the phase space (r cos ;r sin ), with xed points depicted as black dots. Relative equilibria along sin = 0 are depicted as solid black lines. The right panel depicts the solution curves in the (r; _ r) space for various values of the integral of motionQ. For all cases V = 1 and = 2:0. (a) Static gain. Two center-type xed points occur. Unconditional convergence is achieved for all initial conditions. (b) Proportional gain with ` = 6:5. Two center-type and two saddle-type xed points occur. The homoclinic orbits associated with the saddle xed points are highlighted in black. Conditional convergence is achieved for initial conditions inside the closed-loop of these orbits. (c) Proportional gain with ` = 5:4. No xed points occur and no convergence. (d) Inverse gain with ` = 5:8. Two center-type xed points occur. Unconditional convergence is achieved for all initial conditions. . . . . . . . . . . . . 64 6.5 Behavior of the mobile sensor in the wake of an oscillating airfoil in a uniform incom- ing ow of strengthU. The airfoil oscillates at dimensionless frequency St (Strouhal number) and amplitude A. A wake forms downstream of the oscillating airfoil, with regions of spatiotemporally-correlated coherent vortical structures. (a) Snapshot of vorticity eld eld!(x;t) for Re = 500, St = 0:075 and A = 1:66. (b) the magnitude eld m(x) of the spectral vorticity, (c) the phase eld (x), (d) the direction eld r=krk and the associated alignment error (x). . . . . . . . . . . . . . . . . . . 70 6.6 Trajectories of the mobile sensor in the wake of the oscillating airfoil for the param- eter values shown in gure 6.5: (a) static gain law, (b) proportional gain law, and (c) inversely proportional gain law. In all cases, the parameters of the mobile sensor are set to G o = 1 and V = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 7.1 The agent located at x = xe 1 +ye 2 is tasked with nding a target located at the origin. A body frame (b 1 ; b 2 ) is axed to the agent and the agent's heading is denoted . The distance r and bearing are computed relative to the agent. The agent's motion is controlled by changing its forward speed V and its rotation rate . 75 7.2 The perception space of the agent. The eld of view and success boundaries are de- picted as red and green, respectively. The perception space is discretized into states S =fs i g corresponding to the cells depicted with black lines. A given states i (high- lighted in black) is acted upon by action a k yielding s k i (highlighted in blue). The compute the intersection of the mapped state with all reference states to determine a discrete probability transition matrix. . . . . . . . . . . . . . . . . . . . . . . . . . 76 7.3 Top row: Trajectories sampled from policies at various stages of training (N train = 100; 1000; 4000). Bottom row: Value function after full training. (a) Model-based policy optimization using Value Iteration. Agent takes path directly to target. (b) Model-free policy optimization using SARSA. In early stages of training, agent takes circuitous path and sometimes does not reach target. As training progresses, the target is reached every time. (c) Model-free policy optimization using Q-Learning. Policy is learned early on and does not deviate much from optimal policy. . . . . . 83 LIST OF FIGURES xi 7.4 Variation in performance parameters as a function of training episodesN train . (a) Prob- ability of failure p fail for SARSA (blue) and Q-Learning (red), with smoothing window size 1000. Optimal value of p fail = 0:0494 shown as a dotted line. (b) Relative value function errorE value . (c) Fraction of suboptimal state-action pairsF suboptimal . . . . 84 List of Tables 3.1 RMS errors for strain rates (, ), vorticity (!), ow type (), and intensity () for time series in Figure 3.12. The sensor separation distancel is varied between 0.1 and 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.1 Fifteen ow simulations were performed at dierent values of the Strouhal number St resulting in three distinct classes of wakes: 2P+4S, 2P+2S, and 2S. Each simulation was assigned an ID value that identies its wake type and Strouhal number. . . . . 35 4.2 The neural network is trained on ve distinct data sets. In each data set, four simulations from each wake type are selected for training the neural network and one simulation is reserved for testing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 xii Chapter 1 Introduction and Motivation There is, of late, tremendous interest in the engineering community { and indeed in the population at large { in developing systems with high levels of autonomy. The Encyclopdia Britannica denes automation as \the application of machines to tasks once performed by humans or, increasingly, to tasks that would otherwise be impossible" [41]. The eld of automation began with mechanization, the use of machines to supplement or replace human labor, allowing humans to increase strength, speed, and productivity in a variety of industrial and practical applications. With the invention of the computer, humans automated the task of calculation. When sensing technology was developed, we invented feedback control systems. Within the last few decades, advances in computing power have accelerated the drive to develop every increasingly anthropomorphic engineering systems, which promise to relieve humans from performing tasks that are mundane or dangerous. Some of the earliest developments of `robots' were computer numerically controlled (CNC) ma- chine tools developed in the 1950s. These machines took advances in feedback control and computer computation to convert sets of programmable instructions to machine controls to automate various industrial processes [83]. In the 1950s, the U.S. Navy developed the Remote Underwater Manip- ulator (RUM) as a means of replacing human divers with submersible `telechirics' { robots that could manipulate objects underwater with similar capabilities as the human hand. In the 1970s, the Navy began to study building `teleoperators' (submersible robots) with various forms of articial intelligence, to reduce risk to life and cost of using divers. A technology report from that era stated that although developments at the time in the elds of articial intelligence had \implications" for underseas operations, these were \less well understood" [98]. Over the years, varying attempts at higher levels of autonomy have been attempted. In the aviation industry, autopilot technologies have been in use for nearly a century. Indeed, automated systems for steering sea-faring vehicles with varying levels of capability have been available since the late 19th century. Early attempts simply reduced the workload of the helmsman by keeping the ship on course. More recently, there have been a number of developments leveraging various forms of `intelligent control' to more eectively cope with the diculties of ship-steering in a manner that is at least as eective as the eorts of an experienced helmsman [89]. The technology to engineer systems with higher forms of autonomy with ever-more anthropo- morphic qualities is in a state of rapid expansion. Whether or not programmed machines can ever achieve the type of `intelligence' that human beings enjoy is a subject of great debate (see, e.g. [115, 55]). These questions aside, there is no doubt that many advances in technology have al- 1 CHAPTER 1. INTRODUCTION AND MOTIVATION 2 low us to emulate with greater and greater delity some of the higher-level tasks that we observe in the natural world. These include, among others, deep learning, computer vision, adaptive feedback control, and the advances in data science that have enabled `machine learning'. In addition, there is a tremendous in ux of nancial and human capital both within and outside the scientic com- munity dedicated to further developing machines capable of assisting or replacing the human. This eort, in recent years, has been directed primarily at the development of autonomous automobiles. Whether or not these cars become technologically and politically viable, these will undoubtedly advance our capability to develop submersible vehicles with high levels of autonomy. Motivation There are a number of applications for autonomous seafaring vehicles. Some of the most famous early examples of unmanned research vehicles are those built by the Woods Hole Oceanographic Institute during the latter parts of the 20th century. The rst vehicle, the Acoustically Navigated Geological Underwater Survey, ANGUS for short, and later Argo, were involved in the famous discovery of the wreckage of the Titanic on the sea oor of the Atlantic Ocean. These vehicles were fairly primitive when it came to automation and had to be tethered to `mother ships', but they allowed the exploration of regions of the ocean inaccessible to humans and are responsible for many discoveries that have fundamentally changed our scientic understanding of the deep ocean [48]. More recently, The Boeing Company has developed a unmanned underwater vehicle (UUV) that has a range of 7500 miles and uses a hybrid rechargeable battery. The Echo Voyager can be used to inspect underwater infrastructure, map the sea oor, or help with ocean exploration [24]. Addi- tionally, Boston Engineering is developing a bioinspired UUV called GhostSwimmer that resembles a large sh both in shape and in swimming style. This new type of UUV is well suited for surveil- lance and reconnaissance missions as well as hull inspection tasks in harbors, in particular due to its novel combination of autonomy and bioinspired propulsion [42]. Finally, the National Aeronautics and Space Administration (NASA) is developing an autonomous submersible to explore the liquid hydrocarbon oceans of Saturn's moon, Titan. NASA wants to send a submarine to Titan's largest sea, Kraken Mare, to carry out scientic investigations under the surface. The Titan Submarine will be fully autonomous, because real-time interplanetary command and control is neither feasible nor desirable [81]. These are just some among many instances where various engineering entities, both public and private, are recognizing the need to develop autonomous underwater vehicles. The Department of Defense itself is also highly interested in developing autonomous ships to operate at sea. One such concept that is reaching the later stages of development is an autonomous surface ship christened the \Sea Hunter" by the Navy. The Anti-Submarine Warfare Continuous Trail Unmanned Vessel (ACTUV) is an unmanned platform whose purpose is to track enemy submarines in the open ocean [82]. The advantages in terms of range, endurance, and reduced risk to Navy personnel have motivated the development of a system capable of operating completely without human input. The vehicle is said to be capable of \safe navigation, autonomous system management for operational reliability, and autonomous interactions with an intelligent adversary". Although, for obvious reasons, the full extent of the ship's capabilities are not public knowledge, it is known that the system will use `non-conventional sensor technologies' for tracking quiet submarine targets in the open ocean [117]. Flow sensing The modern revolution in autonomous systems is, at least in part, due to the number of dierent and novel sensing modalities available. An autonomous submersible might be equipped with several dierent types, such as visual, acoustic, chemical, or hydrodynamic sensors. CHAPTER 1. INTRODUCTION AND MOTIVATION 3 Figure 1.1: Von K arm an vortex street in the wake behind Heard Island in the southern Indian Ocean. This image is a modied version of the original captured by the Moderate Resolution Imaging Spectroradiometer (MODIS) and is in the public domain courtesy of NASA. In terms of visual sensing, a camera that operates in the visual spectrum, or perhaps a forward- looking infrared (FLIR) sensitive camera might be used. The sub would likely be equipped with external lighting to actively illuminate the environment. Submarines have a long history of being tted with acoustic sensing capabilities, both active and passive. Sonar is widely used by submarines to survey the underwater environment. More recently, there has been a great deal of interest in outtting autonomous submersibles with chemical sensors. This is largely driven by the mission concept for these vehicles: tracking plumes of chemicals released into the environment (see e.g. [32]). These three sensing modalities have existing sensor technology and their use is fairly widespread. As mentioned in the context of the Navy's autonomus sub hunting ship, there is a renewed interest in using `non-conventional sensor technologies' to probe the underwater environment. Dur- ing the Cold War, the Soviet Union lacked the capability to detect submarines in the open ocean using passive acoustic measurements, unlike the United States which had developed the U.S. Sound Surveillance System (SOSUS) [79]. To counter this, the Soviets developed the \System Obnaruje- nia Kilvaternovo Sleda" or \wake object detection system," which in part relied upon detecting the wake generated by a submarine as it moved through the water [46]. Although the U.S. felt that \[i]t is unlikely any of these methods will enable detection of submarines at long ranges" in a 1974 assessment of the technology [87], much of the research on these technologies remains redacted in the intelligence documentation that has recently (2017) been released to the public. The fact that much of this information remains classied suggests that the U.S. Department of Defense maintains a strong interest in developing these technologies and has a reasonable hope that hydrodynamic sensing has great promise as a means of probing the ocean for vehicle signatures. The primary reason why there is great hope for using hydrodynamics as a means of probing the ocean is because in recent decades there has been a number of discoveries within the uid dynamics community that support the idea that hydrodynamic signatures are both strongly coherent (well- organized) and very persistent (long-lasting). Any uid dynamicist knows that ow over a body generates a discernible and coherent ow pattern. This phenomenon occurs over many orders of magnitude in length scale. At the scale of atmospheric ows, there are often von K arm an vortex streets in the wake behind islands as the wind ows over them. An example of this phenomenon is the wake of Heard Island in the southern Indian Ocean, as shown in Figure 1.1. At a more CHAPTER 1. INTRODUCTION AND MOTIVATION 4 intermediate scale, it is well known that ships generate coherent and persistent wakes as they pass through the ocean. Finally, at a smaller scale, when sh swim through the water, the oscillatory motions they make to propel their bodies through the water generate coherent vortical structures, see e.g. [122, 28]. The classical approach to identication, visualization, and quantication of coherent ow struc- tures relies on global reconstruction of the ow eld by an external observer, using either experi- mental or computational methods, and identication of its topological structure (see e.g. [121, 56, 43, 95]). Using a global approach, there are a number of methods available to identify ow patterns or structures that are important for various engineering applications such as material transport, propulsive performance improvement, and active ow control. One example of this is the discovery of Lagrangian Coherent Structures (LCS) from ow data (see e.g. [45, 96, 97, 57, 44]). This method identies ow patterns that are important for transport processes from a Lagrangian perspective. Alternatively, there are the techniques from modal analysis which identify Eulerian coherent struc- tures. These techniques can be either based on governing physical equations, such as Resolvent (see e.g. [69, 66]) and Global Stability Analyses (see e.g. [77]) or driven by data from empirical or numerical experiments such as Proper Orthogonal Decomposition (see e.g [6]) and Dynamic Mode Decomposition (see e.g. [93, 90]). These methods are widely used for identifying coherent structures that can be used for low order models for ow control. Local sensing Global ow characterization yields valuable insight for engineering design, but it does not provide a suitable framework for an organism or an engineered vehicle with local sensory capabilities. In such applications, one has access to pointwise and discrete sensory measurements in a small subdomain of the ow eld. The objective is to use these local sensory measurements to compute desired ow properties. Two general questions are important to answer here: what are the desired ow properties that need to be computed and how can they be computed? The how part consists of deciding which physical quantities the sensors should measure to obtain enough information, as well as the processing algorithms that decode the desired ow properties from the sensory measurements. To illustrate what distinguishes local sensing from the global techniques, consider the airfoils shown in Figure 1.2. Here, we subject a Joukowski airfoil (see e.g. [129] and references contained therein) to three dierent ow conditions: uniform ow at 0 angle of attack, 15 angle of attack, and 15 angle of attack with a point vortex near the leading edge. Using potential ow theory, we compute the velocity eld subject to the uniform ow condition in the far-eld and the no- penetration boundary condition on the surface of the airfoil. Then, using Bernoulli's equation, we compute the pressure eld from velocity. In the rst row of Figure 1.2, we plot the pressure distribution and in the second row, we plot the velocity distribution, each on the surface of the airfoil. Each surface measurement corresponds to a given ow condition. To obtain the surface measurements, we solved the forward problem: computation of the sensory measurements given the ow properties (boundary conditions, vortex locations, etc.); this process was deterministic and unique. The inverse problem of computing ow properties from sensory measurements is markedly more dicult to solve because it can involve incomplete information and could have zero or multiple solutions. For a well-trained eye, it may be possible to match surface measurements in the rst two rows to the ow elds plotted in the third row. However, any attempt of this type most likely ad hoc and based on heuristics. We are thus inclined to ask: are there more robust and systematic tools for solving this inverse problem? More precisely, we want to ascertain if there is a mathematical mapping from the surface measurements to the ow eld properties. If CHAPTER 1. INTRODUCTION AND MOTIVATION 5 Figure 1.2: Illustration of the dierence between local (surface pressure or velocity) and global (entire ow eld) sensing. An airfoil in uniform ow at (left) 0 angle of attack, (center) 15 angle of attack, and (right) 15 angle of attack with a vortex on the leading edge. We plot the (top) pressure distribution and (middle) velocity distribution on the surface of the airfoil. (bottom) Then, we plot the entire pressure eld and the associated streamlines. Can global ow characteristics be determined from measuring only local ow properties? it does exist, then how do we go about obtaining it and how well-conditioned is it? These are the types of questions we must ask when solving inverse problems. Despite this apparent diculty with solving these ow sensing problems, there is a substan- tial and growing body of evidence shows that identiable wake signatures from islands, sh, and plankton are detected by organisms that live in the water. Detection of hydrodynamic signals is well documented in aquatic organisms at various length and time scales. These organisms `feel' the surrounding uid through specialized sensory modalities such as pressure or velocity sensors and respond accordingly; see, e.g. [100, 16, 17] and references therein. Harbor seals are known to track a moving object by following its wake [25], and sh respond to a wide variety of ow stimuli in behaviors ranging from upstream migration to predator evasion [18, 30, 73, 88]. For copepods, a group of microscopic crustaceans that are the most abundant multicellular organisms on the planet, the ability to respond to hydrodynamic stimuli is essential for their foraging, mating and escape maneuvers [105, 125, 126]. The details of how these organisms measure the relevant hydrodynamic signals, process them, and respond to cues embedded therein depend on complex and intertwined sensory, neural, and actuation networks. In this disseration, we do not attempt to directly model the animals themselves; CHAPTER 1. INTRODUCTION AND MOTIVATION 6 Hydrodynamic Environment Sensing Swimming Dynamics Actuation Signal Processing Other Agents Locomotion Flow Sensing Figure 1.3: Schematic block diagram for autonomous submersible. Sensing { How does the agent extract information from the hydrodynamic environment? Signal Processing { How does the agent extract relevant cues from measured data? Actuation { What actions does the agent take in response to these cues? Swimming Dynamics { How do these actions translate into motion? Hydrodynamic Environment { What information is encoded in the ow eld by the movement of the platform as well as other agents? How is this information measured by the sensors? instead, we are interested in developing problem solving strategies suited for use in engineered devices. However, these aquatic creatures leverage various ow sensing modalities to achieve critical tasks such as navigation, mating, and predation. Because ow sensing is essential for their survival, we believe that it is a viable paradigm for extracting valuable information from uid ows that remains largely unexplored in the engineering community. For these reasons, we are inspired by the abilities of these remarkable organisms to develop the basic scientic principles to allow for localized ow sensing in autonomous submersibles. Physics-based and data-driven modeling strategies In order to design autonomous sub- mersibles with ow sensing capabilities, we will need to model all of the components that make up the larger systems, see Figure 1.3. First, we will need a model of the sensing systems, including their physical layout, what physical quantities they can measure, and the quality of those measure- ments. Then, we will need to build signal processing algorithms to take the sensory measurements and decipher the relevant cues encoded therein. We will spend a major component of this thesis focusing on the development of signal processing algorithms. Next, we will need to determine how we will actuate our system in response to the perceived cues to achieve the desired behavior. We will then need to couple this actuation to the swimming dynamics, producing the locomotive behavior of the system. Finally, we will close the loop through the hydrodynamic environment. We place the hydrodynamic environment at the center of Figure 1.3 because all of the elements of the system couple to it. CHAPTER 1. INTRODUCTION AND MOTIVATION 7 system of interest (a) (b) model model input output input estimate input estimate parameters physics-based data-driven Figure 1.4: (a) Systems of interest have similar forms { an input x is transformed into an output y by the mapping y =f(x). (b) The job of the modeler is to approximate this mapping in one of two ways: Physics-Based Modeling uses physical intuition to construct a model f phys (x) which estimates the output ^ y from input x. Data-Driven Modeling constructs a modelf phys (x;), then uses known input-output pairs to optimize parameters, which allows estimation of the output ^ y from input x. The model of each application and subsystem will have a similar structure. First, each subsystem will have an input or measurement. The model will then take that input and produce an output. This paradigm is illustrated by Figure 1.4(a). The modeling process consists of creating a mapping from input to output that produces the desired eect or behavior. Each subsystem and application will also admit a dierent modeling strategy. We can split the modeling approaches into two categories: physics-based and data-driven, as shown in Figure 1.4(b). In physics-based approach, the input-output model is constructed using our intuition about which physical phenomena are relevant to the situation. For example, this might include using governing equations of the dynamics or leveraging our knowledge about the system to choose a signal pro- cessing scheme. This strategy is appropriate when the physical insight is available and the problem is tractable. On the other hand, there may be cases where the physical insight is unavailable or, perhaps, our knowledge of the physics of the problem is insucient and the problem becomes in- tractable. In these cases, if we can build a sizable library of known input-output pairs, we can use a data-driven modeling approach. For this strategy, we employ a parametric mapping then compute a set of parameters which map our input data to the outputs in some optimal sense of our choosing. Both of these approaches may be useful in dierent situations and it is an active area of research to determine when and why certain approaches are useful for which applications. We employ both physics-based and data-driven modeling in this disseration. Chapter overview In this introductory chapter, we have demonstrated that their is a growing desire and need to develop autonomous engineering systems and, in particular, for autonomous submersibles. These systems will, courtesy of the inevitable synergies with other industries, likely be capable of high levels of autonomy. Due to some of the unique challenges associated with oper- ating in the underwater environment, autonomous submarines will need to leverage unconventional sensing modalities such as hydrodynamic tracking. Although there is a great deal of research as- sociated with global ow pattern detection, new paradigms will need to be developed to support CHAPTER 1. INTRODUCTION AND MOTIVATION 8 localized ow sensing. These new approaches will undoubtedly be inspired by the many organisms that leverage ow sensing daily for their survival. In Chapter 2, we will brie y survey the literature to see what attempts have been made to create engineering solutions for local ow sensing problems. Then, in Chapters 3 and 4, we will cover two case studies of ow sensing, one physics-based and one data driven. Next, in Chapters 5 and 6, we will build closed-loop models with ow sensing and motion planning. Then, in Chapter 7, we will explore new paradigms for control of nonlinear and distributed systems by examining a case study in reinforcement learning. Finally, in Chapter 8, we will conclude the thesis and oer outlook on the future directions for this body of research. Chapter 2 Literature Review In recent years, there have been a number of attempts to incorporate techniques from ow modeling, control theory, signal processing, and machine learning to solve ow sensing problems. Here, we provide a brief survey of a selection of these works which are relevant to the problems addressed in this thesis. Local ow characterization Recent attempts have been made to use ow sensing as a means to characterize ows. Local ow characterization using physics-based methods have included various attempts to dene vortical ow by measuring gradients in the velocity eld. A great deal of eort has been dedicated to identifying vortices, in part because of their signicance as coherent structures in many ows of engineering importance. In the 1950s, Truesdell introduced the kinematic vorticity number [107] as a dimensionless measure of the `quality of rotation' which has been used in uid dynamics (see e.g. [70]), geology (see e.g. [104]), and atmospheric ow science (see e.g. [92]). Later, Hunt et al. [54] introduced the widely-used Q-criterion as a means of delineating vortices. The 2 -criterion was introduced by Jeong and Hussain in 1995 [56]. These denitions satised many requirements of vortex identication but failed in noninertial frames of reference, driving Haller to introduce his objective dention of a vortiex [43]. All of these serve as local measures of ow character based on measurements of the velocity eld. In Chapter 3 of this thesis, we develop similar measures of local ow character and also examine how they can be measured by a sensory array living in the ow and sampling the local velocity eld. Global ow characterization There also have been a number of attempts to use local mea- surements to classify far-eld ow conditions. These include the eorts of Fernandez, where both physics-based and data-driven techniques informed studies of an articial lateral line consisting of an array of pressure sensors. There, the sensory setup is used to perform both shape identication and vortex tracking [33]. Because these model-based results were built on inviscid theory, ow separation made the method limited. However, Maertens and Triantafyllou showed that although they diminished the eectiveness of the ideal ow model, the form of the disturbances due to vis- cous eects are specic to the shape that generated them. Therefore, they proposed that unsteady viscous eects could, in fact, be leveraged to improve object identication [67]. In another series of experiments, Klein and Bleckmann demonstrated that an articial lateral line can be used to esti- mate the presence, size, and position of an upstream cylinder by measuring the the pressure trace 9 CHAPTER 2. LITERATURE REVIEW 10 of the vortex wake that it generates [58]. Finally, in a series of experiments with another articial lateral line, Venturelli et al. demonstrated that von Ka rma n vortex streets could be distingushed from uniform ows. In addition, the orientation and position of the sensory array with respect to the ow and the the vortex street could be estimated [112]. These results served as early motivation for the work contained in this thesis because they demonstrated that characteristics of the far eld could be determined from measurements of local information. Machine learning More recently, there has been an eort to identify wake structures using data-driven methods. Wang and Hemati used an ideal ow model to measure time series of velocity on the surface of a shlike body immersed in a variety of exotic von-K arm an-type vortex street wakes. Then, taking the spectral information from these time series and using a Gaussian t to form `features,' a library of data corresponding to dierent wakes was formed. Then, using a k- nearest-neighbors (KNN) machine learning algorithm, they demonstrated that their wake detection protocol can have accuracy rate of over 95% in dierentiating the wake types they considered [118]. In a study with similar motivations, Tu performed numerical experiments on pitching and plunging plates at low Reynolds number to build a library of wake data. Then, using a variety of statistics of the time series of local velocity measurements in the wake, and using the machine learning technique of linear discriminant analysis (LDA), it was shown that using only point measurements, the wake of a pitching plate can be distinguished from that of a plunging plate [109]. In Chapter 4 of this thesis, we use the machine learning tool of neural networks to classify vortex wakes using local time series of vorticity. Lateral line-inspired ow sensing The sh ability to sense minute water motions in their vicinity is attributed to the lateral-line mechanosensory system [30]. The lateral-line system consists primarily of two types of sensors: supercial and canal neuromasts. Supercial neuromasts are located in the skin and are functionally suited for measuring velocity [61]. Canal neuromasts are recessed in bone-like canals and are thought to measure pressure [21]. The layout of the canal sensory system is highly conserved across diverse species and evolutionary time [20]; namely, it extends over much of the body but is concentrated near the head [99]. This basic arrangement also holds in the blind Mexican cavesh, albeit with more sensory receptors in comparison to their surface-dwelling relatives who retained a functional visual system [128]. The universality of the lateral-line architecture suggests a functional advantage in having a higher density of sensors in the anterior segment of the sh. The fact that sensors are concentrated on the surface of a sh in areas where dierences in pressure are greatest supports the notion that these locations are optimized for ow sensing [119]. In fact, recent experiments by Ristroph et al. correlate anatomical measurements of canal density with pressure measurements from ow experiments [88]. They suggest that the relevant sensory output is the dierence in pressure between the left and right canals, rather than the absolute pressure itself. Because behaviors such as predator/prey detection [116] and rheotaxis [73] are mediated by the lateral lines, there is a great deal of interest in understanding how these sensory organs can be used to provide sensory feedback during maneuvering of the sh. Once the relevant sensory information can be obtained, the desired behavior must be eected by actuating the system appropriately. Salumae et al. used two pressure sensors on the head of a robotic sh to measure lateral pressure dierences to estimate orientation to the incoming ow. This sensory information was fed into a Braitenberg controller, and they showed that the robot is able control its orientation with respect to the ow [91]. Gao and Triantafyllou equipped a towed CHAPTER 2. LITERATURE REVIEW 11 underwater vehicle with an articial lateral line of pressure sensors. They used an unscented Kalman ltering technique informed by a panel-based ow model to estimate the ow direction. This was fed into a feedback control law that allowed the vehicle to control its orientation [36]. DeVries et al. developed a robotic shlike vehicle with an articial lateral line and used estimation and control strategies to allow the vehicle to compute ow properties necessary for feedback control. Using a potential ow model, they developed an algorithm to estimate the owspeed, angle of attack, and the relative position of an upstream obstacle. Then, this information is fed back into the control law to perfom rheotaxis and station-holding [26]. These results demonstrated that shlike behaviors could be achieved by using sensory information inspired by the lateral lines of shes. In Chapter 5 of this thesis, we propose a mathematically tractable model system that allows us to investigate the optimal sensor placement that maximizes a lateral-line inspired sensory output. We then use a kinematic model to investigate how rheotactic behavior emerges in light of this sensory information. Trajectory planning and source seeking Sensing disturbances in a signal eld to locate the source of those disturbances is an often-studied type of problem in the control community called source-seeking. Algorithms developed for source-seeking can be applied to trajectory planning prob- lems (see e.g.[13, 35]) or can be used in plume tracing (see e.g. [32]). Although these foundational works are either agnostic to the signal eld type or assume a chemical signal eld, we are interested in examining situations where the signal eld is embedded in the hydrodynamics itself. One such application of this is formation ight and schooling to determine a follower's position relative to the wake of a leader. Using a lifting-line model and an array of pressure sensors, Hemati et al. demonstrated that various estimation techniques can be used to compute relative positions within the wake [51]. DeVries and Paley used a similar modeling approach and a Bayesian algorithm to estimate the wake parameters of the leader aircraft. Then, using optimal control theory, they devel- oped a controller that drives the follower to track the leader while following a path that maintains the observability (in the control theory sense) of the leader's relative position [27]. These results demonstrated that relatively simple ow eld models can be coupled to techniques from optimal estimation and control to allow a follower to localize a wake-generator by sensing the hydrodynamic signature contained in its wake. Cochran and Krstic developed a source-seeking algorithm that, although not directly tuned to hydrodynamic signals, did establish a framework for understanding how periodic probing of a signal eld could be used to measure the local signal gradient. They used a nonholonomic unicycle model equipped with an extremum seeking controller to steer a vehicle to the source of a scalar signal eld. In Chapter 6 of this thesis, inspired by this work, we will develop a physics-based strategy to allow an agent to track the wake of an airfoil and follow the signal to its source. Reinforcement learning for motion planning of swimmers For problems in distributed actuation and sensing, the relationship between what is sensed and how the system should be actuated is often not known or nding it is not tractable a priori. Recent advances in computer science have introduced a promising new paradigm for control of these complex systems called reinforcement learning [103]. Reinforcement learning methods have been used very successfully in the gaming industry, but only recently began to make their way to control of swimming robotic systems. Gazzola et al. used a system of coupled nite dipoles to model the hydrodynamic interactions of a large school of swimmers. Using a reinforcement learning algorithm, they demonstrated that swimmers can use adaptive decision-making to adjust their gaits in order maintain their relative CHAPTER 2. LITERATURE REVIEW 12 positions within a school-like formation [38]. Building on this work, Novati et al. used reinforcement learning to eciently synchronize the motion of two swimming bodies. They demonstrated that when the follower uses reinforcement learning to adapt its gait, up to 30% reduction in energy ex- penditure for the follower and up to 20% increase in its swimming-eciency can be observed. [78]. Further advancing these ideas, Verma et al. showed that sh can improve their propulsive e- ciency by following vortices in the wake of a leading swimmer. Using deep reinforcement learning, the \smart-swimmer" can adapt its position and gait to intercept the vortices from the wake of the leader [114]. These studies demonstrate that reinforcement learning can be used to develop actuation strategies for swimmers in ow elds where the appropriate strategy is not known a priori. In Chapter 7 of this thesis, we consider the simpler problem of pursuit and target-locating with visual sensing. Applications of these types of problems include rescue missions with vehicles operating under water, in caves or in urban environments. In many cases, Global Positioning System (GPS) is not available and detailed information about the position of the target and the pursuing agent are absent. Using a reinforcement learning approach, we train a nonholonomic agent to steer towards the target. These results demonstrate that reinforcement learning is a paradigm that allows the appropriate behavior to be learned from experience by an agent. Chapter 3 Flow Sensing: Physics-Based Approach The classical approach to ow characterization relies on global reconstruction of the ow eld by an external observer, using either experimental or computational methods, and identication of its global topological structure (e.g. [121, 95]). Another means of ow characterization is to identify the ow `type'. Locally, all incompressible ows lie on a spectrum from extensional to rotational [68, 63]. These types of characterizations are used to delineate coherent ow structures [56, 43]. Global ow characterization yields valuable insight for engineering design, but it does not provide a suitable framework for an organism or an engineered vehicle with local sensory capabilities. In such applications, one has access to pointwise and discrete sensory measurements in a small subdomain of the velocity eld. The objective is to use these local sensory measurements to computed desired ow properties. In this chapter, we develop a framework suitable for a sensor to measure the local velocity eld and provide a local ow `characterization.' Given lack of denite knowledge of which velocity component is sensed in various aquatic organisms, we consider here three cases where the velocity sensors have access to (i) only the orthogonal component of the ambient velocity relative to a frame attached to the sensory array; (ii) only the parallel component of the ambient velocity; and (iii) both orthogonal and parallel components of the velocity, in other words, the full velocity vector. To begin, we consider an array of two sensors separated by a nite distancel. The sensory system has access to one measurement: the dierence in velocity between the two sensors. We couple these sensors to a general incompressible velocity eld, see Figure 3.1. Under the assumption that the ow eld evolves at a larger spatial scale than the dimensions of the sensory array, we linearize the velocity eld around the sensory array. Then, we one-way couple the sensory array to the ow (the sensors do not aect the ow) and calculate the local strain rates and vorticity eld. Next, we consider the case where we know a priori that the local ow eld takes the form of a simple shear ow and we ask whether the sensory system is capable of measuring the shear rate. We consider this simple ow s a rst step that allows us to probe the fundamentally nonlinear nature of the inverse problem. We develop, within the context of this simple sensor-environment model, decoding algorithms that would allow the sensory system to decipher the shear rate and orientation from primitive sensory measurements. We nd that only in the case of full velocity sensors, the sensory system is capable of uniquely determining the characteristics of the ambient shear. 13 CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 14 (a) Global flow features (b) Linearized flow (c) Sensory layout s l s r (d) Sensory system Sensory measurements Decoding algorithms Flow properties Figure 3.1: Problem setup and summary of approach. (a) The goal is to develop a framework to detect select ow features by a local sensory array. (b) The uid velocity eld is linearized around the sensory array. (c) A simple bio-inspired sensory system consists of two sensors and one sensory measurement, the dierence between the two sensors. (d) The sensory measurements are decoded to detect the desired ow properties. Then, we consider a fully general two-dimensional, incompressible ow eld of unknown charac- ter and we ask whether the sensory system is capable of characterizing the local ow and measuring its intensity. Inspired by [63], we introduce the notion of ow `type' and ow `intensity' and show how these ow properties are related to the local strain rates and vorticity. We nd that, under certain conditions on the sensory measurements, we can derive decoding algorithms that allow for reconstruction of the local strain rates and vorticity eld. We then go beyond reconstructing the strain rates and vorticity to identifying the ow type on the spectrum from extensional to rotational ows, thus providing a local characterization of the ow eld. To illustrate these ideas, we virtually place this sensory system in several canonical ow elds and test its performance. Numerical results show that the sensory system and processing scheme we develop accurately characterize the local ow. 3.1 Bio-inspired ow sensing system Consider a uid environment with velocity eld u (x;t), where x is the Euclidean position vector measured relative to a xed inertial frame and t is the time. At a given instant in time, we expand the velocity eld u about a particular location x o , which we will later take to be the location of the sensory system, using the Taylor series u = u o + (x x o )ru + O kx x o k 2 : (3.1) Here, u o is the velocity evaluated at x o , and the velocity gradient tensorru is evaluated at x 0 . In Figure 3.2, we show an example of the spatial decomposition of a nonlinear velocity eld into a uniform component, a spatially-linear component, and a residual eld representing the higher order terms. We discard the higher order terms and employ the linear approximation u u o + (x x o )ru: (3.2) The velocity gradient tensorru is canonically decomposed into a symmetric strain rate tensor S and a skew-symmetric vorticity tensor (see e.g. [56, 43]) ru = S + ; (3.3) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 15 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 Figure 3.2: The Taylor decomposition of (a) a full ow eld into (b) a constant component, (c) a linearly varying component and (d) higher order terms, where the velocity magnitude is shown in color, with blue to red representing low to high speed. (e) Constant and linear components are combined into linearized ow for comparison with full nonlinear ow. where S = 1 2 ru +ru T and = 1 2 ruru T : (3.4) The superscript () T denotes the transpose operation. Given an incompressible, two-dimensional ow eldr u = 0, we can write S = 1 2 and = 1 2 0 ! ! 0 ; (3.5) where and are the strain rates and ! is the vorticity. We then construct a sensory system consisting of two velocity sensors, located at x r and x l a distance l apart, as depicted in Figure 3.3. Let (b 1 ; b 2 ) be an orthonormal frame attached to the system such that x l x r = lb 2 ; let denote the orientation of the b 1 -axis measured from the inertial x-direction. Each sensor samples the velocity eld u in the (b 1 ; b 2 ) frame such that the data from both sensors are combined to produce one sensory measurement s = u(x l ) u(x r ). By virtue of (3.2), s is given by s = u(x l ) u(x r ) =lb 2 ru: (3.6) The goal of the sensory system is to invert this computation to determine the components ofru from sensory output s. The inverse problem, however, is not algebraically closed with one set of measurements. The vector s = s 1 b 1 + s 2 b 2 contains at most two independent quantities s 1 and s 2 , whereasru is a tensor containing three independent quantities (, , and !). For this reason, it is necessary to take more than one set of independent measurements by changing the orientation of the sensory system. Let n and m denote the orientation of the sensory system when the nth andmth samples are taken, respectively. We dene the angle m;n = n m that parameterizes the relative rotation of the sensory array between the two angles (see Figure 3.3). It is worth noting here that the sensory system has access to its relative rotation m;n but is not aware of its orientation n and m in inertial space. Let (b 1 ; b 2 ) n and (b 1 ; b 2 ) m denote the body frames associated with n and m , respectively. Let (ru) n , (s) n and (ru) m , (s) m denote the velocity gradient tensor and sensory measurements in the n and m frames, respectively. The sensory output in the n frame is given by CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 16 Figure 3.3: Two-sensor system with sensory array extentl. Velocity is measured in body frame (b 1 ;b 2 ). Orientation with respect to inertial frame is parameterized by the angle . The relative rotation between angles n and m is given by m;n =nm. (3.6) which, when expressed in component form, yields the following system of linear equations (s 1 ) n (s 2 ) n = l 2 0 1 1 1 0 0 2 4 n n ! n 3 5 : (3.7) Clearly, given the strain rates and vorticity, one can directly evaluate the sensory measurement. However, we emphasize that the inverse problem { computing ( n ; n ;! n ) from ((s 1 ) n ; (s 2 ) n ) { is not algebraically closed. Thus, we need additional sensory measurements ((s 1 ) m ; (s 2 ) m ) at m and we need to relate ( n ; n ;! n ) to ((s 1 ) m ; (s 2 ) m ). The velocity gradients (ru) m and (ru) n are related by a change of coordinates (ru) m = m;n (ru) n T m;n ; (3.8) where m;n is the rotation matrix between n and m frames m;n = cos m;n sin m;n sin m;n cos m;n : (3.9) The sensory measurement ((s 1 ) m ; (s 2 ) m ) in them frame is obtained by substituting (3.8) into (3.6), yielding s m =lb 2 ( m;n (ru) n T m;n ): (3.10) which, upon further simplication, yields the sensory measurement in component form (s 1 ) m (s 2 ) m = l 2 sin 2 m;n cos 2 m;n 1 cos 2 m;n sin 2 m;n 0 2 4 n n ! n 3 5 : (3.11) 3.2 Measuring shear ows using velocity sensors The local approximation of the velocity eld u(x) takes the form of (3.2). To allow us to probe the fundamental mathematical nature of this problem, we begin by assuming that the velocity eld CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 17 (a) orthogonal (b) parallel (c) full velocity Figure 3.4: (a) Sensors that measure s 1 = u 1 (x l )u 1 (xr ) are called orthogonal sensors and (b) sensors that measure s 2 = u 2 (x l ) u 2 (xr ) are called parallel sensors. (c) Full velocity sensors have access to both u 1 and u 2 . By denition, u 1 =u b 1 and u 2 =u b 2 . can be locally approximated by a simple shear such that ru = 0 0 0 ; (3.12) where is the shear rate. We remind the reader that, in general, the velocity gradient tensor of an incompressible velocity eld has three independent components (, , and !). Consider a sensory system located at x o and let (b 1 ; b 2 ) be an orthonormal frame attached to the sensory system. Each sensor samples the velocity eld u. The data from both sensors is combined to produce one sensory measurement s, namely, the dierence in velocity between the left and right sensors s = u(x l ) u(x r ); (3.13) where the velocity at the left and right sensors are denoted u(x l ) and u(x r ), respectively. In component form, one has s = s 1 b 1 + s 2 b 2 and from (3.11), s 1 s 2 = 1 2 l(cos 2 + 1) l sin 2 ; (3.14) where and are measured with respect to the inertial frame. We distinguish three types of sensors: orthogonal sensors that measure s 1 only, parallel sensors that measure s 2 only, and full velocity sensors that measure both s 1 and s 2 , as shown in Figure 3.4. Next, we discuss the sensory measurement for each of these sensor types. Sensory measurements We consider rst the forward problem, namely, the computation of the sensory measurements, s 1 , s 2 or both, given , l and . From (3.14), both s 1 and s 2 scale linearly with bothl and but depend nonlinearly on orientation. Therefore, we examine more closely their dependence on . Orthogonal sensors { Orthogonal sensors measure s 1 = 1 2 l(cos 2 + 1). Clearly, s 1 depends sinusoidally on with period , see Figure 3.5(a). This induces re ection symmetries about the directions 0k 1 =2, where k 1 is an integer, as shown in Figure 3.5(b). In other words, the same sensory measurement is obtained at = o , = + o , and = o , for any given orientation o of the sensory system. Further, one may think of the periodicity in the sensory measurement as a re ection symmetry about an axis oriented at o +=2, CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 18 axis of symm. axis of symm. axis of symm. ax. of symm. Figure 3.5: (a) Symmetries in the orthogonal measurement. Measurements at =o,o, o, and +o all produce the same sensory measurement. (b) This is due to the re ective symmetries in the measurement, with axes of re ection at = 0 and (horizontal axis) as well as ==2 and 3=2 (vertical axis). which rotates with the sensory system. This non-uniqueness in sensory measurements can be viewed as a limitation of the sensory system and will have important implications on the inverse problem of reconstructing ow properties from sensory measurements as we will discuss in the Decoding section. Parallel sensors { Parallel sensors measure s 2 = 1 2 l sin 2, which also depends sinusoidally on with period; see Figure 3.6(a). However, this dependence induces re ection symmetries about the directionsk 2 =4, wherek 2 6= 0 is an integer, as shown in Figure 3.6(b). Similar to the orthogonal case, there are re ection symmetries that are characterized by two perpendic- ular axes of symmetry, however, in this case they are are diagonal, see Figure 3.6(b). That is to say, for a given orientation o , the same sensory measurement is produced at = + o , ==2 o , and = =2 o . Here too, the periodicity in the sensory measurement induces a re ection symmetry about an axis oriented at o +=2, which rotates with the sensory system. Full velocity sensors { The axes of symmetry in the orthogonal and parallel sensors are at a =4 phase relative to each other. This phase dierence between the two sensory components breaks the re ection symmetries about the 0k 1 =2 andk 2 =4 directions in the case of full velocity sensors. However, because both sensory measurements have the same periodicity , the re ection symmetry associated with o +=2 remains in the full velocity measurement. Decoding the sensory measurements The goal of the sensory system is to invert (3.14) to identify the properties of the shear ow, namely, to compute from sensory measurement. The decoding of sensory measurement is what we refer to as the inverse problem. The solution of the inverse problem requires the inversion of the sinusoidal functions in (3.14). Due to this inherent nonlinearity in the forward measurements, there may be some degree of ambiguity in the solutions we obtain for the inverse. We solve the fully nonlinear inverse problem to explore how the nonlinear nature of the measurements impacts the ability to nd solutions to the inverse problem. We consider that the sensory array does not have knowledge of its orientation relative to the inertial frame but is aware of any rotation it makes relative to its body frame. To invert (3.14) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 19 axis of symm. axis of symm. axis of symm. axis of symm. ax. of symm. Figure 3.6: (a) Symmetries in the parallel measurement. Measurements at =o, +o,=2o, and=2o all produce the same sensory measurement. (b) This is due to the re ective symmetries in the measurement, with axes of re ection at ==4 and3=4 as well as = 3=4 and =4 (both diagonal axes). and compute from s 1 , s 2 or both, the sensory system needs to circumvent its lack of knowledge of . One approach is to take multiple sensory measurements at dierent orientations. Algebraic closure of the inverse problem requires a minimum of two independent measurements. Let n and m denote the orientation of the sensory array when the nth and mth samples are taken, repsectively, and let m;n = n m . Substituting into (3.14), we get (s 1 ) m (s 2 ) m = 1 2 l(cos(2 n 2 m;n ) + 1) l sin(2 n 2 m;n ) : (3.15) Expanding the trigonometric functions (s 1 ) m (s 2 ) m = 1 2 l(sin 2 n sin 2 m;n + cos 2 n cos 2 m;n + 1) l(cos 2 n sin 2 m;n sin 2 n cos 2 m;n ) : (3.16) The inversion of (3.16) to compute the shear intensity and sensor orientation n from local sensory measurements will depend on the type of sensors, namely, whether the sensors are orthogonal, parallel or full velocity sensors. In each case, algebraic closure of the inverse problem necessitates at least two independent equations to eliminate n and solve for given the sensory measurements and m;n . Orthogonal sensors { Consider the sensory measurements (s 1 ) n and (s 1 ) n1 taken at n and at n1 = n . Note that we dropped the subscript in n;n1 for simplicity. Using the rst row of (3.16), we have that (s 1 ) n (s 1 ) n1 = 1 2 l(cos 2 n + 1) l(sin 2 n sin 2 + cos 2 n cos 2 + 1) : (3.17) We combine the two rows to eliminate n and solve for , which yields s = (s 1 ) n1 + (s 1 ) n 2 p (s 1 ) n1 (s 1 ) n cos 2 l sin 2 : (3.18) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 20 Figure 3.7: Measured value over true value of for orthogonal sensors. (a) `First' branch on the left for > 0 and (b) `second' branch on the right, also for > 0. If < 0, then the two branches switch. Black lines delineate contours where = 0, which is where the two branches of the solution exchange forms. where s is the value of obtained from the sensory system. Clearly, the inverse mapping from (s 1 ) n , (s 1 ) n1 , and to is not unique. The inverse solution has two branches: a branch associated with the positive sign and another associated with the negative sign. To examine how s deviates from its actual value , we substitute (3.17) into (3.18) to get s = 1 2 sin 2 (cos 2( n ) + cos 2 n + 2)j(cos 2( n ) + cos 2 n + cos 2 + 1)j ; (3.19) where the branches are again chosen by selecting the positive or negative version of. For notational convenience, let =(cos 2( n ) + cos 2 n + cos 2 + 1). Depending on the sign of , we can simplify (3.19) further. When > 0, the two branches of the solution take the form s = 2 cos 2( n ) + 2 cos 2 n + cos + 3 2 sin 2 d and s = 1; (3.20) respectively, for 6= 0 mod . When < 0, the two branches switch places while keeping the same functional forms. Thus, for all, the value of s depends both upon the measurement angle n as well as the rotation angle , as is shown in Figure 3.7. Parallel sensors { Consider the sensory measurements (s 2 ) n and (s 2 ) n1 taken at n and at n1 = n . Using the second row of (3.16), we have that (s 2 ) n (s 2 ) n1 = 1 2 l sin 2 n l(cos 2 n sin 2 sin 2 n cos 2) : (3.21) We combine the two rows to eliminate n and solve for , which yields s = 2 q (s 2 ) 2 n1 + (s 2 ) 2 n 2(s 2 ) n1 (s 2 ) n cos 2 l sin 2 : (3.22) Again, the inverse mapping from (s 2 ) n , (s 2 ) n1 , and to is not unique; it has two branches, each associated with the choice of the positive or negative sign in (3.22). CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 21 Figure 3.8: Measured value over true value of for parallel sensors. Solutions exchange sign at black circles, when n1;n ==2 mod =2. As in the case of orthogonal sensors, we examine how s diers from by substituting (3.21) into (3.22), we nd that s =jj j sin 2j sin 2 : (3.23) For > 0 and 6==2 mod =2, the rst and second branches of the solution take the form s = 1 and s =1: (3.24) For < 0, the rst and second solutions switch place. Thus, the sensed value s captures the correct magnitude of the ambient shear but cannot determine the sign of the shear unambigu- ously. Figure 3.8 illustrates the branches of the solution and how each branch switches sign at = =2 mod =2, implying that each branch may have the wrong sign of the ambient shear. Full velocity sensors { If the sensory system has access to full velocity information, algebraic closure of the system can be achieved from one sensory sample taken at n . To this end, the sensory measurement is given (s 1 ) n (s 2 ) n T . Using both rows of (3.16), we have that (s 1 ) n (s 2 ) n = 1 2 l(cos 2 n + 1) l sin 2 n : (3.25) We combine both rows to eliminate n and solve for , which yields s = (s 1 ) 2 n + (s 2 ) 2 n l(s 1 ) n : (3.26) Unlike the orthogonal and parallel sensors, the inverse in (3.26) is unique. To verify it has the correct functional form, we substitute (3.25) into (3.26), yielding s = 2 cos 2 n cos 2 n + 1 ; (3.27) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 22 which reduces to s = 1; (3.28) for n 6= =2 mod . Thus, the value of s is correct in both magnitude and sign and the sensory system is capable of reconstructing the correct value of the ambient shear. 3.3 Characterizing ow types using velocity sensors In the previous section, we considered the measurement of shear ow using a local sensory system. This measurement of shear is a proxy to more complex local ows and served as a starting point to examine the fully nonlinear nature of the inverse problem. In this section, we now examine a fully general, two-dimensional incompressible ow eld. At a given instant in time, we consider the locally linear approximation of (3.2). The velocity gradient tensorru is canonically decomposed into S and as by (3.3). Given an incompressible, two-dimensional ow eldr u = 0, we can write the components of S and in the (;;!) parameterization given in (3.5). LetkSk = p tr (SS T ) andk k = p tr ( T ) denote the Frobenius norms of S and , respec- tively; then, kSk = 1 p 2 p 2 + 2 and k k = 1 p 2 j!j: (3.29) Note that the Frobenius norm of a tensor is frame invariant, that is to say,kRAk =kAk for any matrix A and rotation matrix R. It follows that any quantities derived from the norm of the tensor are also frame invariant. Since S is symmetric, it can be diagonalized by a rotation matrix R = r 1 r 2 composed of its normalized eigenvectors r 1 = 2 6 6 6 6 4 + p 2 + 2 q 2 ( 2 + 2 ) + 2 p 2 + 2 q 2 ( 2 + 2 ) + 2 p 2 + 2 3 7 7 7 7 5 ; r 2 = 2 6 6 6 6 4 p 2 + 2 q 2 ( 2 + 2 ) 2 p 2 + 2 q 2 ( 2 + 2 ) 2 p 2 + 2 3 7 7 7 7 5 : (3.30) The diagonalized strain rate tensor takes the form ~ S = 1 p 2 kSk 0 0 kSk ; (3.31) where ~ () denotes a quantity written in the eigenbasis of S. Meanwhile, in the case of two- dimensional ows, is invariant under rotation and can be generally written in the form ~ = sgn (!) 1 p 2 0 k k k k 0 : (3.32) To this end, when expressed in the eigenbasis of S, the velocity gradient tensor takes the form r~ u = 1 p 2 kSk sgn (!)k k sgn (!)k k kSk : (3.33) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 23 (a) (b) (c) (d) (e) Rotational Shear Extensional Shear Rotational Figure 3.9: Variation of ow type with parameter for (a) counter-clockwise rotation, (b) counter-clockwise shear, (c) extensional ow, (d) clockwise shear, and (e) clockwise rotation. The (;;!) parameterization of the velocity gradient tensor describes the magnitude of various velocity gradients in the ow: stretch, shear, and vorticity. However, each one fails independently to provide an easy means by which one can identify key characteristics of the ow such as its type. We then pose the question: is there a dierent parameterization where the parameters represent more meaningful quantities with regard to the local ow characterization? We introduce two parameters, the ow type and ow intensity = sgn (!) k k kSk and = p 2 p k k 2 +kSk 2 ; (3.34) where 2 (1;1) describes the ratio of vorticity to strain and 2 [0;1) describes the intensity of velocity gradients. Figure 3.9 illustrates the parameterization of ow type by on a spectrum from extensional to rotational. The velocity gradient tensor can be rewritten in terms of and as r~ u = 2 2 6 6 4 1 p 2 + 1 p 2 + 1 p 2 + 1 1 p 2 + 1 3 7 7 5 : (3.35) The (;;!) parameterization of the velocity gradient tensor for a two-dimensional incompressible ow is a three-parameter description of the gradients in the ow. We translate these parameters into type and intensity by = ! p 2 + 2 and = p 2 + 2 +! 2 : (3.36) In addition, we also consider the rotation angle that diagonalized the strain rate tensor S, where tan = + p 2 + 2 : (3.37) The (;; ) parameterization has a distinct advantage over that of (;;!) in that each parameter has a intuitive signicance that can be ascertained independently of the others. Perhaps more importantly, the values of and are frame independent, since they are derived fromkSk and k k, whereas (;;!) are fundamentally linked to a given frame of reference. In Figure 3.10, we demonstrate the advantage of this parameterization by examining the spatial variation of (;; ) for three canonical uid ows: a Lamb{Oseen vortex, a Gaussian shear prole, CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 24 -2 0 2 -2 0 2 -2 0 2 -2 0 2 -2 0 2 -2 0 2 -2 0 2 -2 -1 0 1 2 0 1 2 3 0 2 4 6 8 10 -0.5 0 0.5 / Figure 3.10: (a{d) Lamb{Oseen vortex, (e{h) Gaussian shear and (i{l) cylinder in uniform ow at Re D = 60. (a,e,i) Velocity vector eld u and vorticity eld !. (b,f,j) Flow type, with the black contours depicting =1 (c,g,k) Flow intensity , (d,h,l) Flow orientation r 1 vector eld, with the colormap showing the orientation angle =. The cylinder data are taken from [12]. and a circular cylinder in uniform free stream at Re y = 60. The Lamb{Oseen vortex is given by the velocity eld u = (1 exp(r 2 =r 2 c ))A=r, where u is the azimuthal velocity, r = p x 2 +y 2 is the radial dimension, r c = 2 is the vortex core size, and A =3 is the vortex strength, with all quantities dened in polar coordinates. The Gaussian shear ow is dened be the velocity eld u = B exp(y 2 =y 2 c ), where u is the velocity in the x direction, y c = 5 is the wake width, and B = 5 is the wake intensity. The cylinder data is a Navier-Stokes simulation of uniform ow past a circular cylinder at Re y = 60, originally published in [12]. The rst row in Figure 3.10 depicts the velocity vectors and vorticity elds. The second and third rows show that provides an intensity- independent measure of ow type and, accordingly, provides a type-independent measure of ow intensity. The fourth row depicts the values of and the associated eigenvector r 1 , which provides the local orientation in the direction of maximum stretch. As an aside, for a two dimensional ow, the type and intensity can be related to the canonical CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 25 \Q-criterion" [54]. The Q-criterion is often cited as a simple means of identication of vortical structures in a ow and is computed by Q = 1 2 k k 2 kSk 2 : (3.38) In the type-intensity parameterization, Q is given by Q = 2 4 2 1 2 + 1 : (3.39) There are two important items to note regarding the (;; ) parameterization and theQ-criterion. First, the numerical value ofQ is inherently tied to both the type and intensity of the ow, whereas and independently reveal the type and intensity. Secondly, since Q depends on 2 rather than directly, the Q criterion reveals type of ow up to a sign. It cannot distinguish the sense of the vorticity, however. We have shown that (;; ) provides a convenient parameterization of the velocity gradient tensorru and that (;;!) can be used to compute this parameterization. Next, we devise a scheme by which a sensory system that measures local velocity dierences can characterize the ow type and intensity. Sensory measurements and decoding Consider a sensory system consisting of two velocity sensors, located at x r and x l a distance l apart, as depicted in Figure 3.3. Let (b 1 ; b 2 ) be an orthonormal frame attached to the system such that x l x r = lb 2 ; let denote the orientation of the b 1 -axis measured from the inertial x-direction. Each sensor samples the velocity eld u in the (b 1 ; b 2 ) frame such that the data from both sensors are combined to produce one sensory measurement s = u(x l ) u(x r ). Let n and m denote the orientation of the sensory system when the nth and mth samples are taken, respectively. We dene the angle m;n = n m that parameterizes the relative rotation of the sensory array between the two angles (see Figure 3.3). It is worth noting here that the sensory system has access to its relative rotation m;n but is not aware of its orientation n and m in inertial space. Let (b 1 ; b 2 ) n and (b 1 ; b 2 ) m denote the body frames associated with n and m , respectively. Let (ru) n , (s) n and (ru) m , (s) m denote the velocity gradient tensor and sensory measurements in the n and m frames, respectively. We nd that, s is given in component form by (3.11). We next discuss how to combine (3.7) and (3.11) to compute the desired ow characteristics from sensory measurements. This problem will depend on the capabilities of the sensory system, namely whether the velocity sensors can access s 1 , s 2 , or both. Orthogonal sensors { Consider sensors that measure s 1 only. We refer to these sensors as orthogonal sensors because they measure the component of velocity orthogonal to the sen- sory array (see Figure 3.4(a)). Barring rank deciency, the algebraic closure of the system necessitates, at a minimum, three sensory samples are taken at n, n 1 and n 2. That is to say, our sensory information consists of the vector (s 1 ) n (s 1 ) n1 (s 1 ) n2 T . Using the rst row of (3.7) and (3.11), we have that 2 4 (s 1 ) n (s 1 ) n1 (s 1 ) n2 3 5 = M ortho 2 4 n n ! n 3 5 ; (3.40) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 26 where the measurement matrix M ortho is M ortho = l 2 2 4 0 1 1 sin 2 n1;n cos 2 n1;n 1 sin 2 n2;n cos 2 n2;n 1 3 5 : (3.41) If n1;n ; n2;n 6= 0 mod and n1;n 6= n2;n , M ortho has full column rank. Thus, since we can compute the inverse (omitted here), we nd that the local strain rates n and n and the local vorticity ! n can be determined by the system. However, we note that the determination here relies on measurements at three distinct orientations, which may be disadvantageous in a time-dependent velocity eld. Parallel sensors { We now consider sensors that measure s 2 only. We denote these sensors as parallel sensors because they measure the component of velocity parallel to the sensory array (see Figure 3.4(b)). Once again, barring rank deciency, the algebraic closure requires that our sensory information consists of the vector (s 2 ) n (s 2 ) n1 (s 2 ) n2 T . Using the second row of (3.7) and (3.11), we have that 2 4 (s 2 ) n (s 2 ) n1 (s 2 ) n2 3 5 = M parallel 2 4 n n ! n 3 5 ; (3.42) where the the measurement matrix M parallel is M parallel = l 2 2 4 1 0 0 cos 2 n1;n sin 2 n1;n 0 cos 2 n2;n sin 2 n2;n 0 3 5 : (3.43) M parallel has maximum column rank 2, and therefore, the velocity gradient cannot be uniquely determined. Here, we are unable to determine ! n in the case of parallel sensors because it lies in the null space of M parallel . Note that this result holds true regardless of how many measurements are taken. Full velocity sensors { If the sensory system has access to full velocity data, we will see that we only need to use two sample measurements atn andn1. We concatenate the orthogonal and parallel sensory data at n and n 1 into the vector (s 1 ) n (s 2 ) n (s 1 ) n1 (s 2 ) n1 T . For simplicity, we dene n1;n =. Using both rows of (3.7) and (3.11), we have that 2 6 6 4 (s 1 ) n (s 2 ) n (s 1 ) n1 (s 2 ) n1 3 7 7 5 = M full 2 4 n n ! n 3 5 ; (3.44) where the measurement matrix M full is M full = l 2 2 6 6 4 0 1 1 1 0 0 sin 2 cos 2 1 cos 2 sin 2 0 3 7 7 5 : (3.45) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 27 For 6= 0, M full has full column rank. If the ow eld u were truly linear, then n , n , and ! n could be uniquely determined from these sensory data. In the general case, however, the sensory data will also include O(kx x o k 2 ) terms from the velocity eld (3.1), and the linear system (3.44) is therefore overdetermined. Using the Moore{Penrose pseudoinverse [40], we nd that the least-squares solution is given by 2 4 n n ! n 3 5 = M + full 2 6 6 4 (s 1 ) n (s 2 ) n (s 1 ) n1 (s 2 ) n1 3 7 7 5 ; (3.46) where M + full = (M T full M full ) 1 M T full , yielding M + full = 1 2l 2 4 sin 2 cos 2 3 sin 2 cos 2 + 1 cos 2 + 1 (cos 2 + 1) cot cos 2 1 (cos 2 3) cot 2 2 cot 2 2 cot 3 5 : (3.47) Thus, for 6= 0, the local strain rates n and n and the local vorticity ! n can be obtained from the sensory system. We now demonstrate how the sensory information can be used to characterize the ow locally. To this end, we will test the ability of the algorithm described in (3.47) to compute ow properties, while varying sensor parameters, in several canonical ow regimes. Gaussian shear { The nondimensional ow prole for a Gaussian Shear is described by the velocity eld u(x) = exp(y 2 )e 1 : (3.48) Accordingly, the ow type and intensity are given by =sgn(y) and = 2 p 2jyj exp(y 2 ): (3.49) We plot ow type in Figure 3.11(a) and intensity in Figure 3.11(b) as a function of the vertical coordinatey. We then apply the algorithm described in (3.47), computing these ow properties rst for xed l and varying , and then vice versa. We see that the algorithm performs very well in capturing both ow type and intensity, and in fact, changing has very little impact on the performance. We also see that yet again, forl O(1), the algorithm performs very well, and for larger l, the performance deteriorates. This result is quite understandable, given the fact that our primary assumption in the development of (3.47) was that the ow was locally linear, and for larger l, the assumption of locality is violated. However, since we did not place any restrictions on the size of , we should expect that even for nite values of , we can accurately recover the ow properties. Lamb-Oseen vortex { The nondimensional ow prole for a Lamb-Oseen vortex is described by the velocity eld u(x) = 1 exp(r 2 ) r 2 (ye 1 +xe 2 ); (3.50) CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 28 -5 0 5 0 1 2 Gaussian Shear -10 0 -10 0 0 1 2 -10 0 10 -10 0 10 Lamb-Oseen V ortex 10 10 (a) (b) (c) (d) -5 0 5 Figure 3.11: (a,c) Flow type and (b,d) ow intensity measurements for varying sensor parameters in nondimen- sionalized (a,b) Gaussian shear prole and (c,d) Lamb-Oseen vortex. Analytical solutions shown in dashed black lines. First plots show xed sensor length l = 1 and varying =f0:1; 0:3; 0:5g. Second plots show xed = 0:1 and varying l =f0:5; 1; 5g. where r = p x 2 +y 2 . Accordingly, the ow type and intensity are given by = r 2 r 2 exp(r 2 ) + 1 and = 2 exp(r 2 ) r 2 p 2r 4 2r 2 exp(r 2 ) + 2r 2 + exp(2r 2 ) + 1: (3.51) We plot ow type in Figure 3.11(c) and intensity in Figure 3.11(d) as a function of the vertical coordinatey. We then apply the algorithm described in (3.47), computing these ow properties rst for xed l and varying , and then vice versa. As was the case with the Gaussian shear, we see that the algorithm performs very well. We conrm that in this canonical ow, the impact of changing is minimal, as well. Much the same, we also see that the algorithm performs well for nontrivial l, but for l > 1, the performance deterioriates. Cylinder in uniform ow { To test the performance in a time-varying ow, we use the com- putational model system of a cylinder in uniform ow shown in Figure 3.10. We place the sensory system consisting of full velocity sensors in its wake. More specically, we place two full velocity sensors at a separation distance l at a distance 9 units downstream from the cylinder's center, along the wake's axis. We allow the ow to evolve in time and measure the velocity at each sensor and calculate the velocity dierence. At each time step, we obtain the sensory measurements at orientation and. We then feed these sensory measurements to the pseudoinverse M + full and plot its output values for strain rates and vorticity. Results are shown in Figure 3.12(a), (b) and (c) for = 0:05. Making use of (3.36) and (3.37), we see that from n , n , and ! n , we can compute n to characterize the local ow type and CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 29 Figure 3.12: Velocity gradients and ow parameters as function of time for sensory array placed atxo = 9e 1 +0e 2 in the cylinder ow, as shown in gure 3.10. True values for velocity gradients and parameters are computed via nite dierence and shown as dashed black lines. Computed values using M + full are shown for (a) , (b) , (c) !, (d) , and (e) . Sensor system parameters are = 0:05 and l =f0:1; 0:5; 1g. The data are adapted from [12]. l ! 0.1 0.0492 0.0669 0.0667 0.379 0.0549 0.5 0.0471 0.0689 0.0638 0.389 0.0522 1 0.0407 0.0755 0.0559 0.429 0.0467 Table 3.1: RMS errors for strain rates (, ), vorticity (!), ow type (), and intensity () for time series in Figure 3.12. The sensor separation distance l is varied between 0.1 and 1. n to determine its local intensity (Figure 3.12(d) and (e)). We also plot the true values of , , !, and evaluated directly at the location of the sensory array. By comparing these results, we see that the pseudoinverse provides reasonable estimates for all quantities. For quantiative comparison, we compute the RMS error (e.g. = 1=T R T 0 ( sensor ) 2 dt) for each time series for each value ofl, with errors tabulated in Table 3.1. We emphasize that the ow type and intensity are frame-independent quantities characteristic of the ow itself and not the frame of reference of the sensory array. CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 30 3.4 Discussion We considered a bioinspired sensory system that measures dierences in local velocities. We placed the sensory system in simple shear and derived conditions under which the system is able to de- termine the shear ow properties from local sensory measurements. In particular, we examined orthogonal, parallel and full velocity sensors that respectively measure the orthogonal, parallel or full components of the velocity vector in a frame of reference attached to the sensory system. We found that the orthogonal sensors cannot determine the shear properties and the parallel sensors can determine the magnitude of the ambient shear but not its direction. Full velocity sensors are able to unambiguously identify the properties of the ambient shear. Then, we developed a framework for locally characterizing general two-dimensional, incompress- ible ows by examining the velocity gradient tensorru. By re-parameterizingru in terms of ow type and intensity , we showed that all ows can be locally characterized on a spectrum from extensional to rotational. We illustrated the use of these parameters in three canonical ows: a Lamb{Oseen vortex, a Gaussian shear prole, and a circular cylinder in uniform free stream at Re y = 60. We applied our sensory setup to the problem of characterizing the local ow eld and examined its performance when measuring either the orthogonal, parallel, or full components of the velocity. We found that parallel sensors are not capable of uniquely determining all ow properties, irrespec- tive of how many times they sample the ambient ow. More specically, we found that using parallel velocity sensors, the sensory system is unable to calculate the local vorticity. We then showed that, when the sensors have access to the orthogonal component or full velocity vector, the sensory system is able to unambiguously compute all three ow properties, however in the case of the full velocity vector, this takes only two independent measurements. We derived a `lter' that solves the inverse problem to obtain the local strain rates and vorticity from the sensory measurements. Finally, we demonstrated how this lter can be used in conjunction with the (;; ) parameterization of the velocity gradient tensor. To illustrate the capabilities of this sensory system, we used it to probe several canonical ow elds. We showed that the sensory scheme we developed provided good agreement with direct measurements of the velocity gradients, types, and intensities of the ow. This framework provides an eective method of measuring strain rates and vorticity from local sensory information and for computing the local ow type. However, we must take careful note of several limitations. First, we assumed that the sensor does not disturb the ambient ow. We also linearized the ow eld around the sensor. This assumption is only valid in certain ow regimes when the dimensions l of the submerged sensory system are small relative to that of the ow structures. It is worth noting that, to generate the numerical results in Figure 3.12, we relaxed this requirement by taking the sensor dimension to be equal to 1 (the same as the cylinder's radius) and obtained good agreement between the actual and sensed ow properties. However, the results for the Gaussian shear and Lamb-Oseen vortex demonstrate the eective performance of the sensory system whenl is up to the same order of magnitude as the ow structures, demonstrating that the use of l.O(1) is reasonable. Further, there is an inherent assumption regarding the timescale of the sensory system and the timescale at which the ow eld is evolving. Because we sample the ow twice at n and n to reconstruct the ow properties, the frequency at which we sample must be much greater than that at which the velocity eld uctuates. Finally, the relative angle must be properly tuned to the ow eld that is being sensed. When choosing , one must be aware that if is too small, the computation of strain rates and vorticity becomes ill-conditioned. In sum, careful selection of the sensory length scale l and relative sampling angle will be important to CHAPTER 3. FLOW SENSING: PHYSICS-BASED APPROACH 31 ensure the sensed ow properties are close to the actual ones. We posit that, for a given ow eld, these parameters can be nely tuned to provide optimal ow information and that the parameter studies we performed with the canonical ows can provide guidance on these choices. Although the analyses performed here were illustrated with two-dimensional ows, the use of and to uniquely characterize ows holds for the three-dimensional case. For any velocity eld, there are three tensor invariants that can be computed forru. For incompressible ow, one of these is always 0, since S is traceless. The other invariants can be uniquely mapped to and . For both the two- and three-dimensional cases, the remaining independent parameters ofru are due to the orientations of the principle axes of strain and rotation. For the 2D case, we denote this , the angle of the principle axes of strain { the principle axis or rotation is always orthogonal to the plane. For the 3D case, there are six independent parameters describing the orientations of the principle axes of strain and rotation. In both these cases, we see that ow type and intensity uniquely characterize all ows, up to their principle axes of strain and rotation. In future studies, we will develop quantiable metrics for assessing the performance of a given sensory array in relation to the specic ow eld it senses. Dening quantitative metrics of sensory performance will be an indispensable step for assessing and comparing the performance of distinct sensory arrays and sensing strategies (static versus dynamic), as well as for optimizing the layout of the sensory platform. In dynamically deployed sensors, as is the case in biological swimmers, we will investigate how the sensory information can be incorporated in a feedback loop to determine the motion of the sensory array. In particular, it will become increasingly important to develop strategies to more intelligently guide the dynamic deployment of the sensory array, using techniques from estimation theory (e.g. [113]), machine learning (e.g. [38]), and compressive sensing (e.g. [11]). Chapter 4 Flow Sensing: Data-Driven Approach The spatiotemporal organization of vorticity in the unsteady wake of a moving body is intimately linked to the parameters of the body and its motion. Classical examples include the von K arm an vortex street behind a stationary cylinder in a uniform free stream and the more \exotic" vortex wakes behind oscillating cylinders [121] and airfoils [9, 10, 39, 59, 60, 64, 95], as well as the wake structures of swimming sh [76, 110, 129] and apping wings [85]. Vortical wakes are often described by the number of vortices shed per oscillation or apping period [121]; 2S refers to wakes in which two vortices of opposite sign are shed per period whereas in 2P wakes, two vortex pairs are shed per period. Classifying the wake type (2S, 2P, etc.) as a function of the body parameters is relevant to many applications including the design of oshore structures and the analysis of propulsive forces in aquatic and aerial organisms and bio-inspired vehicles. The classical approach to wake classication relies on global inspection of the ow eld by an external observer in order to identify global topological patterns [95, 102, 121] or Lagrangian coherent structures [43, 56]. This global approach provides valuable insight for engineering design and analysis but it is not suited for applications where one has access to local information only. The objective of this study is to classify global ow patterns using local ow measurements. To address this problem, we use supervised machine learning and neural networks. Neural networks are versatile and powerful tools that have been applied to many branches of engineering and science ranging from autonomous vehicles [23] to skin cancer classication [31]. This versatility is due to the fact that neural networks, with as few as one hidden layer, can approximate any continuous function with arbitrary accuracy [52]. In addition, recent advances in computational capabilities and fast algorithms [50, 120] made it relatively straighforward to design, construct and train neural networks. In this section, we design a neural network that classies the wake type from local time mea- surements of the vorticity eld in the wake of an oscillating airfoil; see Figure 4.1. Specically, we train the neural network to distinguish between three wake types: 2S, 2P+2S, and 2P+4S. We analyze the performance of the trained network in terms of its accuracy in predicting the wakes of test simulations not used for training and we discuss the salient features relevant to classify each wake type. 32 CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 33 visual inspection measure time series identify coherent structures neural network Global Local wake classification Figure 4.1: Classication of ow patterns behind an oscillating airfoil: (left) proposed approach using only local measurements, and (right) classic approach relying on global inspection of the spatiotemporal patterns in the ow eld. 4.1 Vortex wakes of pitching airfoils We consider a symmetric airfoil of chord C and diameter D undergoing a simple pitching motion of the form (t) = max sin 2ft; (4.1) where is the angle of attack, max is the maximum angle of the airfoil, f is the frequency of oscillation, and t is time; see Figure 4.1. This system can be characterized by three dimensionless parameters: the Strouhal number St, which is a dimensionless measure of the oscillation frequency, the dimensionless amplitude A of the airfoil trailing edge, and the Reynold number Re, St = fD U ; A = 2C sin max D ; Re = UD : (4.2) Here, and are the uid density and dynamic viscosity, and U is the velocity of the uniform free stream. For Reynolds numbers Re between 10 2 and 10 4 , the qualitative structure of the wake is independent of Re and depends on the Strouhal number St and dimensionless oscillation amplitude A [95]. The Strouhal number is inversely proportional to the oscillation period, which is typically an integer multiple of the vortex shedding period. At higher St, the number of vortices shed per period is generally two, a set of oppositely-signed vortices, and the wake type is 2S. As St decreases, the number of vortices shed per period increases and more exotic wakes are observed, including 2P+2S and 2P+4S wakes [95]. To obtain numerical data of the wake, we solve the incompressible Navier-Stokes equations, @u @t + uru =rp + 1 Re u; r u = 0; (4.3) in the uid domain surrounding the oscillating airfoil. Here, the uid velocity eld u and the pressure eld p depend on space, parameterized by the position vector x, and time t. Further, u is subject to the no-slip condition at the airfoil boundary. Since the airfoil motion is prescribed, CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 34 Figure 4.2: Wake behind an oscillating airfoil for A = 1:656 and St ranging from 0:038 to 0:079 as shown in (a). Three wake types are observed, 2P+4S, 2P+2S, and 2S, shown respectively by plotting contours of vorticity eld !(x;t) for (b) St = 0:044, (c) St = 0:064 and (d) St = 0:075. Time traces of vorticity at downstream locations I, II, and III are plotted in gure 4.3. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 35 Wake Type 2P+4S 2P+2S 2S Case ID St Case ID St Case ID St A1 0.038 B1 0.063 C1 0.071 A2 0.040 B2 0.064 C2 0.073 A3 0.042 B3 0.065 C3 0.075 A4 0.044 B4 0.066 C4 0.077 A5 0.046 B5 0.067 C5 0.079 Table 4.1: Fifteen ow simulations were performed at dierent values of the Strouhal number St resulting in three distinct classes of wakes: 2P+4S, 2P+2S, and 2S. Each simulation was assigned an ID value that identies its wake type and Strouhal number. the uid-structure interactions are one-way coupled: the airfoil aects the ow but not vice-versa. To solve for u(x;t), we use a numerical algorithm based on the Immersed Boundary Method, initially devised for fully-coupled uid-structure interactions [84]. The algorithm is described in detail in [72], and has been implemented, optimized and tested extensively by the group of Haibo Dong; see, e.g., [111, 8]. We then calculate the vorticity eld ! = (r u) e 3 , where e 3 = e 1 e 2 is the unit normal vector to the plane. The vorticity can be expressed as a scalar eld !(x;t) given the two-dimensional nature of the problem. We conducted a total of fteen ow simulations by varying St from 0:038 to 0:079 while xing A = 1:656 and Re = 500; the values of St are shown in Table 4.2. Three dierent types of periodic wakes are observed: 2P+4S, 2P+2S, and 2S, each occurring at ve distinct values of St as shown in Figure 4.2(a). These wake types are identied by global inspection of spatiotemporal patterns of the vorticity eld!. Representative examples of these wakes are shown by plotting the contours of the vorticity eld in Figures 4.2(b)-(d). Local ow measurements We measure the vorticity eld !(x;t) at a position x m in the wake using a sampling interval t. Examples of the time evolution of the local vorticity are shown in Figure 4.3. We represent the local vorticity measurements as a column vector s m of lengthN, where N is the total number of samples. In other words, the n th entry of the measurement vector s m is equal to !(x m ;nt). Our goal is to train a neural network that takes as input the measurement vector s m , also called a snapshot, and classies the wake type into 2P+4S, 2P+2S or 2S. We construct a training dataset S2R NM of M snapshots S = 2 4 j j j s 1 s 2 s M j j j 3 5 : (4.4) We then assign to each snapshot s m a label y m that identies the wake type or class. The three classes 2P+ 4S, 2P+ 2S, and 2S are represented, respectively, by the one-hot encoding y m = 2 4 1 0 0 3 5 ; 2 4 0 1 0 3 5 ; and 2 4 0 0 1 3 5 : (4.5) Here, y m is a column vector of length N c = 3, where N c is the number of classes. The class label CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 36 -20 0 20 -20 0 20 0 1 2 -20 0 20 0 1 2 0 1 2 Figure 4.3: Time traces of vorticity in the wake of an oscillating airfoil for wake types 2P+4S, 2P+2S, and 2S. Point-measured time traces are taken at several locations downstream along the lines indicated by I, II and III in Figure 4.2. matrix Y2R NcM for the training set can be written as Y = 2 4 j j j y 1 y 2 y M j j j 3 5 : (4.6) For each ow simulation, we probe the wake at select locations x m chosen according to an equally-spaced grid of 51121 points. Ten snapshots are recorded at each grid point as follows. First, the time series of the local vorticity eld is sampled at t = 0:012 forN = 210, resulting in a total time of approximately 2:5 units. We then shift the time at which the rst sample is recorded by about 0:25 time units to create a total of ten distinct snapshots. That is to say, from each ow simulation, we extract a total of 5112110 = 61,710 snapshots, resulting in 1561,710 = 925,650 snapshots from all simulations. 4.2 Neural networks for wake classication A neural network is formally dened as a mapping ^ y =f(s;) from an input s2R N to an estimate ^ y2R Nc of the wake type. Here, denotes the mapping parameters. Note that the hat is used to distinguish the estimate ^ y from the actual label y. Also note that we temporarily drop the subscript () m for clarity. The mapping f is constructed by the successive application of various functions, each one is called a layer. The particular functional form for each layer and the number of layers depend on the specic application [94]. In this work, we use the feedforward neural network shown in Figure 4.4. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 37 Neural Nework Linear Linear Activation Normalization Figure 4.4: Architecture of the neural network consisting of four successive layers: a linear layer, a nonlinear activation layer, another linear layer, and a nonlinear normalization layer. Network architecture We design a feedforward neural network of four layers. The rst layer is an ane transformation q = W s + b; (4.7) where q2R N f is the output vector, W2R N f N is the weight matrix, b2R N f is the bias vector, and N f is the number of features. We interpret the rows of W as features of the vorticity time series that are relevant for dierentiating between classes. Consequently, Ws is a projection of the input s onto these features. The second layer is called an activation layer and is dened by the nonlinear function ~ q = max(0; q) that sets negative entries of q to zero and `activates' only positive entries. In component form, one has ~ q j = ( q j q j 0; 0 q j < 0; j = 1;:::;N f : (4.8) The third layer is another ane transformation ~ ~ q = ~ W ~ q + ~ b; (4.9) where ~ ~ q2 R Nc , ~ W2 R NcN f , and ~ b2 R Nc . Finally, the last layer is a normalization layer ^ y = softmax( ~ ~ q), which can be written in component form as ^ y c = exp( ~ ~ q c ) P c exp( ~ ~ q c ) ; c = 1;:::;N c : (4.10) This layer normalizes the output vector ^ y such that its entries take values between 0 and 1 and its L 1 -normk^ yk = 1. The output ^ y is a probability distribution vector, where the c th entry ^ y c represents the likelihood that s belongs to class c. The mapping ^ y = f(s;) is constructed by successive application of these four layers. The parameter vector consists of all the entries of the weight matrices W and ~ W and bias vectors b and ~ b. That is, 2R Np is a column vector of dimension N p =N f (N c +N + 1) +N c : (4.11) CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 38 training inputs Neural Network training outputs Loss Function output estimates Compute Gradient Stochastic Gradient Descent Figure 4.5: A neural network is trained by optimizing a loss function over the network parameters subject to a training set of snapshots sm and labels y m . The optimization problem is solved using stochastic gradient descent. Network training The training algorithm, depicted schematically in Figure 4.5, optimizes the values of subject to a loss function of the form l(y; ^ y) = Nc X c=1 y c ln ^ y c (1 y c ) ln(1 ^ y c ): (4.12) This loss function is zero only when ^ y = y and nonnegative otherwise. The logarithmic form of the function means that the loss increases dramatically as the estimate ^ y deviates further from y. Given the training dataset of inputs s m 2 S and outputs y m 2 Y. The optimization problem can be written as = argmin 1 M M X m=1 l m (y m ;f(s m ;)); (4.13) where l m is the loss associated with s m as in (4.12) and the total loss is an average of the losses over all snapshots in the training set. To solve (4.13), we use a stochastic gradient descent method, new = old @ @ 1 M X m l(y m ;f(s m ; old ) ! ; (4.14) where @()=@ is the gradient of the loss function and is a scalar called the learning rate. Equa- tion (4.14) is applied iteratively until the change in loss reaches an acceptable threshold. In practice, in order to ensure that the solution to (4.13) is statistically signicant, we need to have N p M, so that the problem is suciently overdetermined. Numerical implementation We implement the neural network using Wolfram Mathematica version 11.1.1. The network is created using the NetChain function, where the number of layers as well as N, N f , N c are dened. The linear layers are created using the LinearLayer function, the activation layer using the ElementwiseLayer function with the Ramp option, and the normalization layer using SoftmaxLayer. This network is then trained using the NetTrain function and the training set S together with its label set Y. After training, we evaluate the accuracy of the trained network on a distinct test dataset using the function ClassifierMeasurements. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 39 0 10 20 30 40 10 -3 10 -2 10 -1 10 0 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 Data Set 1 Data Set 2 Data Set 3 Data Set 4 Data Set 5 Figure 4.6: Error incurred by ve dierent trained networks in predicting the wake type or class of the test simulations. The error is plotted as a function of the number of featuresN f for (a) all classes, as well as for the individual classes (b) 2P+4S, (c) 2P+2S, and (d) 2S. 4.3 Network performance We train the neural network on ve distinct training sets obtained by arranging the fteen ow simulations into ve combinations; for each wake type, four simulations are chosen for training and one simulation is reserved for testing. The simulations reserved for testing are highlighted in Table 4.2. In other words, in each combination, 80% of the data was used for training and the remaining 20% was used to analyze the performance of the trained network. Network accuracy We evaluate the performance of the trained network for each combination by measuring its accuracy in correctly classifying the data reserved for testing. Figure 4.6 shows the error (= 1accuracy=100) in log-scale as a function of the number of features N f for all combinations. The overall error for all classes is depicted in Figure 4.6(a) and the error per class in Training Data Test Data Data Set 1 A2, A3, A4, A5 A1 B1, B2, B4, B5 B3 C1, C2, C3, C4 C5 Data Set 2 A1, A2, A4, A5 A3 B1, B2, B3, B4 B5 C1, C3, C4, C5 C2 Data Set 3 A1, A2, A3, A4 A5 B1, B2, B3, B5 B4 C2, C3, C4, C5 C1 Data Set 4 A1, A3, A4, A5 A2 B2, B3, B4, B5 B1 C1, C2, C3, C5 C4 Data Set 5 A1, A2, A3, A5 A4 B1, B3, B4, B5 B2 C1, C2, C3, C5 C4 Table 4.2: The neural network is trained on ve distinct data sets. In each data set, four simulations from each wake type are selected for training the neural network and one simulation is reserved for testing. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 40 0.12 0.89 91.59 1.35 97.21 0.74 98.53 1.90 7.67 Data Set 1 2P+4S 2P+2S 2S Actual 2S 2P+2S 2P+4S 1.23 2.82 98.55 1.95 90.99 0.02 96.81 6.19 1.43 Data Set 2 2S 2P+2S 2P+4S 0.16 0.02 96.60 4.01 98.01 0.34 95.83 1.97 3.06 Data Set 3 2S 2P+2S 2P+4S Predicted 0.53 1.20 99.86 0.95 91.51 0.02 98.52 7.29 0.12 Data Set 4 2S 2P+2S 2P+4S 0.41 0.70 99.68 1.09 98.19 0.08 98.50 1.11 0.24 Data Set 5 2S 2P+2S 2P+4S Figure 4.7: Accuracy of the trained network in predicting the ow type of the test simulations. Diagonal entries show percentages of correct classication and o-diagonal entries show percentages of wrong classication and to which class they were erroneously assigned. Figures 4.6(b)-(d). Combination 5 exhibits the best overall performance for N f 10, with errors as small as 0.01 and smallest error at N f = 34. For comparison purposes, a purely random guess among three classes would incur an error of about 2=3. For the cases where the network fails, we track the fraction of cases assigned to each of the erroneous wake types or classes. The results are depicted in confusion matrices in Figure 4.7 for all combinations. The diagonal entries of these matrices show the accuracy of the network in correctly classifying to wake type and correspond to the results in Figures 4.6(b)-(d). The o-diagonal entries quantify the degree of confusion between classes in cases of erroneous classication. For example, the rst row, third column of the confusion matrix shows the percentage of cases that belong to 2S but were misclassied as 2P+4S. Obviously, in all combinations and for all classes the accuracy of predicting the correct class is above 90% and can be as high as 98 or 99%. Note that higher degrees of confusion are observed between 2S, 2P+4S wakes and 2P+2S, 2P+4S wakes than between 2S, 2P+2S wakes. Also, in nearly every combination, there is a higher percentage of 2P+2S cases incorrectly classied as 2P+4S than 2S. From a physical standpoint, the vortex laments in 2S and 2P+2S wakes (see Figure 4.2(c) and (d)) can be potentially confused with those of the 2P+4S wakes (see Figure 4.2(b)). Further, since 2P+2S and 2P+4S wakes are produced at decreasing frequencies of airfoil oscillations, the overall intensity of the vorticity eld is lower, leading to potentially higher confusion between these wake types than with 2S wakes. Network robustness We evaluate the robustness of the network to three corruptions of the sensory measurements: (i) we add measurement noise, (ii) we truncate the vorticity time series so that the measurement duration is smaller than the period of airfoil oscillation, and (iii) we reduce the sampling rate of the sensory measurements; Illustrative examples of the modied measurements are shown in Figure 4.8. To add noise, we compute the magnitude of the signal at a given sensor location as the root mean square value of local vorticity over time, s rms = v u u t 1 N N X n=1 s 2 n : (4.15) We then average this value over all the time series from a given simulation to compute s mean . We CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 41 0 0.5 1 1.5 2 -9 -6 -3 0 3 6 9 Original Noisy 0 0.5 1 1.5 2 Original Truncated 0 0.5 1 1.5 2 Original Downsampled Figure 4.8: Samples of the time traces of vorticity obtained after altering the data by (a) adding white noise, (b) truncating the total time interval of measurement below the oscillation period of the airfoil, and (c) decreasing the sampling rate of the sensor below the computational rate 1=t, where t is the simulation time step. add measurement noise to the local signal s noisy = s + w; (4.16) where w2R N is a random vector whose elements are drawn from a Gaussian distribution of zero mean and variance given by Var(w) = ((NSR)s mean ) 2 ; (4.17) where NSR is the noise-to-signal ratio. As the noise-to-signal ratio increases, we observe a monotonic loss in the network prediction accuracy as shown in Figure 4.9(a). As the noise-to-signal ratio approaches unity, the network performance degrades to about 80% accuracy, indicating that the neural network is still able to decode the underlying patterns in the data that make the classes discernible. Next, we truncate the time series by decreasing the time length N, providing the network with less information. We quantify this by dening the signal truncation ratio (NN new )=N and let it vary from 0 to 0.5. The motivation is to consider how much information from the full oscillation period the network needs to correctly identify the wake class. Figure 4.9(b) shows that the network maintains an accuracy above 90% until approximately one third of the time period of oscillation is removed. Finally, we vary the sampling rate by dening the downsampling factor t new =t, and allow it to vary from 1 to 21. As we downsample more agressively, the vorticity is sampled at progressively lower resolution providing the network with less information. Figure 4.9(c) shows that the accuracy decreases slowly indicating that the essential information is maintained even at lower sampling rates. Even at 10% of the original sampling rate, we observe a 90% accuracy which demonstrates that relevant patterns in the data are preserved and the network is capable of identifying them. Salient features To analyze the features relevant for dierentiating between classes, we consider the trained network based on Data Set 5 and N f = 34. The network features are the rows of W. We refer to the k th row W kn , n = 1;:::;N, as feature k. Examples of the features arrived at by the trained network are plotted versus time in Figure 4.10. While some features bear resemblance to the snapshots shown in Figure 4.3, most features are not recognizable by direct inspection of the vorticity time series. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 42 noise-to-signal ratio signal truncation ratio downsampling factor accuracy/100 1.0 0.9 0.8 0.7 0.6 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.1 0.2 0.3 0.4 0.5 1 5 9 13 17 21 Figure 4.9: Accuracy of the trained network in predicting the ow type of the test simulations with (a) noisy measurements, (b) reduced measurement intervals, and (c) reduced sampling rate. The rst layer of the trained network projects a snapshot onto the features of the network, thus quantifying the prevalence of each feature in the snapshot. A feature acts as either a positive or negative identier. If a feature is a positive identier for a given class, its prevalence in a snapshot indicates a higher likelihood that the snapshot belongs to that class. Conversely, if a feature is a negative identier for this class, its prevalence in a snapshot indicates a lower likelihood that the snapshot belongs to that class. The likelihood ^ y c of class c is given by (4.10). Before normalization, its value is determined by Numerator(^ y c ) = exp X ~ W ck ~ q k + ~ b c = Y k exp( ~ W ck ~ q k ) exp( ~ b c ) (4.18) From (4.7) and (4.8), we know that ~ q k = max(0; W kn s n + b k ) is zero or positive, re ecting which features are active in ~ q k . By substituting into (4.18), it is straightforward to see that the weights ~ W ck dictate how important feature k is for classc. With this understanding we dene the salience matrix ck , ck = ~ W ck q P n W 2 kn max k ~ W ck q P n W 2 kn ; (4.19) where ck 2 [1; 1] and q P n W 2 kn is the L 2 -norm of feature k. For similar measures of salience, see [5]. The saliency matrix of combination 5 is shown in Figure 4.10(a). In Figure 4.10(b), we rank the features by their salience and explore which features are most salient for each class. In particular, we plot ten features per class: ve corresponding to highest positive salience and ve to lowest negative salience. Some features such as F25 or F32 can be physically interpreted as counting the number of vortices passing over the vorticity sensor, but others do not lend themselves to such intuitive interpretation. Further, features that are most positively salient for one class are also most negatively salient for another class. This duality reinforces the notion that features act as positive and negative identiers; if a feature is very CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 43 F33 2P+4S 0.98 F13 0.97 F26 0.87 F4 0.85 F14 0.79 F22 -0.85 F34 -0.93 F20 -0.96 F18 -0.99 F8 -1.00 F19 2P+2S 0.55 F24 0.53 F31 0.53 F25 0.51 F8 0.44 F17 -0.64 F14 -0.67 F5 -0.73 F32 -0.77 F26 -1.00 F11 2S 0.42 F32 0.38 F27 0.37 F29 0.35 F17 0.35 F28 -0.59 F25 -0.62 F12 -0.62 F4 -0.89 F9 -1.00 2P+4S 2P+2S 2S F5 F10 F15 F20 F25 F30 1 0 -1 1 0 -1 1 0 -1 1 0 -1 1 0 -1 1 0 -1 1 0 -1 1 0 -1 1 0 -1 1 0 -1 0 1 2 0 1 2 0 1 2 (b) 1 0 -1 (a) Figure 4.10: (a) Salience matrix of the N f = 34 features of the model for Nc = 3 classes. (b) The ve most positively and negatively salient features for each class are plotted as functions of time and ranked by salience. Features that occur as a pair of positive and negative identiers for dierent classes are highlighted using the same color. These features are also marked by colored boxes above their corresponding columns in the salience matrix. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 44 important for identifying a snapshot as belonging to a certain class, it is also likely that it is important for identifying the snapshot as not belonging to the other classes. Figure 4.10(b) shows that for the class 2P+4S, most of the highly salient features, with salience values close to 1 and -1, do not appear with opposite level of salience in other classes, explaining why this class resulted in the most confusion. Mapping onto the physical space We examine the performance of the trained network on the wakes reserved for testing by mapping the network accuracy onto the physical space. In Figure 4.11, we plot the probability of correct classication as a function of space for the entire domain of the wake. To be precise, we test on a larger region of the wake than the region we trained on, showing that the data that we used to train is somehow characteristic of the entire wake itself, thus supporting our hypothesis that local measurements contain enough information to classify the wake. We also show that, save for some small areas downstream of the airfoil, the network has a very high probability of correct prediction, showing that the relevant information about vortical coherent structures is embedded in the time history of local vorticity throughout the wake. The performance is degraded in the large regions (the dark conic-shaped zones around the airfoil) where there is essentially no information; the vorticity time series are always zero in these regions. 4.4 Discussion How should one process local measurements to infer global information about the ow? This question leads to a multi-valued inverse problem that is dicult to solve in general. Here, we used neural networks to develop a processing algorithm that decodes the wake type from local vorticity measurements behind an oscillating airfoil. One of the advantages of this data-driven approach is that it requires no a priori knowledge of the specic features of the decoding algorithm. The patterns in the time series of the local vorticity (Figure 4.3) that are most useful for dierentiating wake types are not readily recognizable by direct inspection. Yet, the neural network uncovers the salient features of these wakes (Figure 4.10) and maps a time series of the local vorticity to the overall wake type; the resulting mapping is shown to be both remarkably accurate and robust to measurement degradation. A few comments on the limitations of the feedforward neural network employed here are in order. This simple network is not shift-independent. Namely, if a certain pattern occurs in a time series, the network will not recognize that it is the same pattern if it gets shifted in time. To overcome this issue, convolutional neural networks will be used in future extensions of this work. Another extension that would allow us to capture more complex relationships between input and output is to increase the number of hidden linear layers beyond one, creating a so-called deep neural network. We deliberately relied on simple local measurements, in the form of vorticity time series, in order to focus on the classication algorithm. In the biological systems that motivated this work, the physiological modalities involved in ow sensing are more complex, with potentially distributed sensors tuned to the hydrodynamic cues relevant for the particular organism. These cues could be in the form of changes in the uid velocity, pressure or density. Flow sensors could independently or jointly probe one or more of these elds. For example, the lateral line system of sh consists of two types of sensors that measure velocity [61] and pressure [21] along the sh body. Models of uid- structure interactions in individual and schools of sh as well as behavior-based sh models often CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 45 Figure 4.11: Comparison of the vorticity eld and the performance of the neural network in correctly classifying the wake as a function of the positionxm where the vorticity eld is probed. (left) Contours of the vorticity eld!(x;t) repeated from Figure 4.2 to facilitate comparison and (right) the spatial distribution of the likelihood ^ y of correct classication for wake types (a) 2P+4S at St = 0:042, (b) 2P+2S at St = 0:066, and (c) 2S at St = 0:073. CHAPTER 4. FLOW SENSING: DATA-DRIVEN APPROACH 46 ignore this ow sensing ability; see, e.g., [4, 34, 108] and references therein. A notable exception is our recent work that uses pressure sensors and postulates a behavior-based model coupled to a physics-based ow model to study sh rheotaxis, or alignment to an oncoming ow eld, [18]. It will be of major interest to train neural networks to detect ows using pressure and velocity sensors and to implement the trained networks in behavior- and physics-based models of sh swimming in isolation and in groups. Chapter 5 Flow Sensing in Fish-Inspired Behavior The sh ability to sense minute water motions in their vicinity is attributed to the lateral-line mechanosensory system [30]. The lateral line system consists primarily of two types of sensors: supercial and canal neuromasts. Supercial neuromasts are located in the skin and are functionally suited for measuring velocity [61]. Canal neuromasts are recessed in bone-like canals and are thought to measure pressure [21]. The layout of the canal sensory system is highly conserved across diverse species and evolutionary time [20]; namely, it extends over much of the body but is concentrated near the head [99]; see Figure 5.1. This basic arrangement also holds in the blind Mexican cavesh, albeit with more sensory receptors in comparison to their surface-dwelling relatives who retained a functional visual system [128]. The universality of the lateral line architecture suggests a functional advantage in having higher density of sensors in the anterior segment of the sh. Experiments by Ristroph, Liao and Zhang correlate anatomical measurements of canal density with pressure measurements from ow exper- iments [88]. They suggest that the relevant sensory output is the dierence in pressure between the left and right canals, rather than the absolute pressure itself. They further show that the canal system is concentrated at locations on the body that experience strong variations in this pressure dierence. While these ndings provide valuable insights into the relevant quantity that is sensed, a deeper understanding of how behavior emerges from sensory information requires, in addition to identifying the sensory output, a knowledge about how such sensory information is processed and translated into decision making and behavior. Discerning how sensory signals translate into behav- ior in live sh is a complicated task. Behavior is the result of intertwined complex mechanisms, where sensory feedback and neural control, muscular actuation and pre exes, and hydrodynamics are combined to produce the remarkable ability of sh to respond to hydrodynamic cues. In this chapter, we propose a mathematically-tractable model system that allows us to shed light on some of these questions. Namely, we investigate the optimal sensor placement that maximizes the sensory output, that is, the dierence in pressure between the left and right sensors, and how behavior emerges in light of this sensory information. Our model consists of a planar elliptic body swimming at a constant speed along its major axis, and equipped with pressure sensors located laterally on the two sides of the ellipse in a symmetric conguration about its major axis, see Figure 5.1(c). The sensors provide one sensory output, namely, the dierence in pressure between 47 CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 48 Side view Canal Neuromasts Top view Model left canal right canal Fish Position along body Catfish Cave Fish Gar Herring Sea Bass Trout Average (c) Canal density (b) 0.0 0.1 0.2 0.3 0.4 0.5 (a) 0.6 0.4 0.2 0.0 Figure 5.1: (a) Spatial distribution of the lateral-line canal neuromasts in Astyanax fasciatus (Blind Mexican cave sh), inspired from [74]. (b) Canal density in various species of sh, adapted from [88]. (c) Lateral-line canal seen from the top and corresponding model system proposed here. the body's two sides, which is consistent with experimental ndings of [88]. We postulate that the body responds to the sensory output by adjusting its orientation. Our goal is to test whether this hypothesis or `control law' is sucient to capture the rheotactic behavior exhibited by many sh species. To this end, we place the sh model in oncoming uniform ows in an otherwise unbounded domain of inviscid and incompressible uid. We investigate optimal sensor placement in these ows and examine the response of the body to pressure cues. Interestingly, the body responds by aligning with or against the ow as observed in sh rheotaxis. This rheotactic behavior depends crucially on the orientation of the shlike body and its swimming speed relative to the oncoming ow. We conclude by commenting on the relevance of these ndings to understanding both the sensory response in sh and the conceptual design of engineered underwater sensory systems. 5.1 Potential ow model for pressure sensing Consider a planar elliptic body with semi-major axis a and semi-minor axis b, submerged in an otherwise unbounded domain of inviscid and incompressible uid. Let (x c ;y c ) denote the position of the center of the body with respect to a xed inertial frame with Cartesian coordinates (x;y). Let denote the orientation angle measured from the positive x-direction to the ellipse's major axis as shown in Figure 5.2(a). The linear and angular velocities of the body are given by ( _ x c ; _ y c ) and _ , respectively. Here, the dot _ ( ) corresponds to derivative with respect to time t. To emulate the sh lateral line, we equip the body with two pressure sensors located at (x l ;y l ) and (x r ;y r ) in a symmetric conguration about the major axis of the ellipse. These sensors capture the local value of the pressure eld p(x;y) to provide one sensory output, namely, the dierence in pressure p =p(x l ;y l )p(x r ;y r ). We investigate the dependence of the sensory output p on body orientation and swimming velocityV in uniform ows of speedU in the negativex-direction. Analytic expressions for the uid velocity and pressure elds are obtained by conformally mapping the innite uid domain outside the elliptic body into the domain outside a circle of radius R = (a +b)=2. For convenience, we introduce the complex coordinate z =x + iy, where i = p 1 is the imaginary unit, in the physical CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 49 physical plane body-fixed coordinates map conformally to circle pane circle plane transform to body-fixed frame left sensor right sensor (a) (b) (c) Figure 5.2: Schematic of the model system (a) in the physical plane, (b) in a body-xed frame, with dashed lines corresponding to level sets of constant elliptic coordinates; and (c) conformally mapped to the circle plane, with dashed lines corresponding to level sets of constant polar coordinates. plane, such that the location of the center of the body is given by z c =x c + iy c and the locations of the sensors are z r = x r + iy r and z l = x l + iy l . These coordinates are then transformed to a body-xed frame, as in Figure 5.2(b), prior to mapping the ellipse plane conformally to the circle plane, as in Figure 5.2(c). The conformal mapping from the circle to the physical plane takes the form (see, e.g., [71, 129]), z =z c + + 2 e i : (5.1) Here, = +i is the complex variable in the circle plane and = ( p a 2 b 2 )=2 is a measure of the ellipse's eccentricity, whereas e () denotes the exponential function. Let u(z) = u x (x;y) + iu y (x;y) denote the complex uid velocity in the physical plane. The corresponding complex potential function can be written asF (z) = (x;y)+i (x;y), with being the real potential function and the stream function. The existence of and is guaranteed by the irrotationality and incompressibility conditions; see, e.g., [65]. By linearity of the problem, the complex potentialF (z) can be written as a superposition of two contributions, F =F b +F u , where F b (z) =V (R 2 2 )= is due to the motion of the elliptic body andF u (z) =Ue i Ue i R 2 = is due to the external uniform ow U. Note that we neglected the eect of the turning rate _ on F b , considering only quasi-static changes in body orientation to best emulate the experimental procedure of [88]. We introduce non-dimensional counterparts to the above quantities using the characteristic length scale and time scale =U. The non-dimensional radius R of the circle is given byR = (a+b)=(2) = p (a +b)=(ab), and all parameters and variables are now considered in non-dimensional form. To this end, the dimensionless stream function takes the form F =F b +F u =V 1R 2 e i R 2 e i : (5.2) The complex conjugatesu() =u iu andu(z) =u x iu y of the uid velocity eld in the circle and physical planes are given by u() = dF () d ; u(z) = dF d d dz : (5.3) The pressure eld p is obtained from the steady Bernoulli equation pp 1 = 1 2 juj 2 ; where p 1 is the pressure at innity and is undetermined. The sensory output is given by the dierence in CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 50 left right Average canal density Small angle limit Large limit 40 20 10 5 0 -1 3 2 1 0.0 0.2 0.4 0.15 0.10 0.05 0.0 0.1 0.2 0.3 0.4 0.5 12 8 4 0 0.4 0.2 0.0 0.2 0.4 0.1 0.2 0.3 0.4 0.5 0.2 0.4 0.5 0.2 0.4 0.5 Stimulation Figure 5.3: (a) Pressure distributions along the body for b=a = 0:266, = 5 , V = 0. (b) Pressure dierence p in (5.5) as a function of for = 5 ; 10 ; and 20 and V = 0. (c) Stimulation (p=) for = 5 ; 10 ; and 20 and its theoretical limit obtained by using small in (5.5) compared with the average canal density of Figure 5.1(b). (d) Pressure dierence as as a function of for V = 5; 10; 20; 40 and = 5 . (e) p=V for same parameter values. pressure between the left and right sensors l and r , p =p( l )p( r ) = 1 2 juj 2 r 1 2 juj 2 l : (5.4) The location of the right/left sensors in the circle plane is expressed as r =Re i and l =Re i , with being the sensors angular position, see Figure 5.2(c). Alternatively, one could use the arc- lengths normalized by the lengthL, as done in [88], see Figure 5.1(a). In the context of the elliptic model, one can put the normalized arc-length s=L in one-to-one correspondence with , as noted in Figure 5.3. Hereafter, we use to refer to the location of the sensors. 5.2 Optimal sensor layout We examine the dependence of the sensory output p on the location of the sensors as the body changes orientation relative to the oncoming ow. Ristroph et al. [88] showed by experiment that the sensor locations on the body of a sh correspond to areas where pressure changes are greatest, suggesting they are optimized to sensing lateral dierences in pressure. Here, we explore these ideas within the context our our analytical model and we compare our ndings to those experimental results. One of the advantages of the modeling approach employed here is that it is amenable to analytic expressions for the pressure (expression not shown) and the pressure dierence p. One gets, upon substituting (5.2) into (5.3) and using the resulting expression in (5.4), that p = 2R 2 sin(2) R 4 2R 2 cos(2) + 1 h V (R 2 1) sin() +R 2 sin(2) i : (5.5) Note that the dependence of the function p on is separable from its dependence on V and . This property has important consequences on optimal sensor placement. Figure 5.3(a) shows the pressure distribution for the right and left sensors as a function of when the body is oriented at = 10 . The ow-facing (right) side experiences higher pressure than the leeward (left) side, consistent with [88, Figure 3]. Figure 5.3(b) depicts p as a function of CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 51 6 4 2 0 0.5 0.3 0.1 2.0 1.5 1.0 0.5 0.0 8 6 4 2 0 -2 -4 -6 -8 0 0.5π π 1.5π 2π 6 0 1 2 3 4 5 pitchfork bifurcation stable stable unstable stable unstable unstable stable stable (a) (b) Figure 5.4: (a) Pressure dierence p as a function of body orientation, swimming speedV and sensor location. Black contour lines show p = 0. (b) p = 0 are xed points of the governing equations for which the body does not change orientation. The stability of these xed directions depend on the swimming speed V . The ellipse aspect ratio is set to b=a = 0:3, but the results hold qualitatively for all 0<b=a< 1. for three dierent body orientations = 5; 10 and 20 . The pressure dierence increases with increasing but the dependence on follows a similar trend for all, which is to be expected from the separability property of p in (5.5). The increase in the sensory signal p as increases can be interpreted by the swimmer as an indication of the degree of misalignment with the oncoming ow, as noted in [112, 88]. The qualitative agreement of the results in Figure 5.3(a,b) with the experimental ndings in gure 3(a,b) of [88] is remarkable (the gures of Ristroph et al. are not reproduced here). To further explore the accord between our model and the experimental results, we plot in Figure 5.3(c) the pressure stimulation parameter, dened in Ristroph et al. as the ratio of the pressure dierence to orientation angle p=. The curves in Figure 5.3(c) collapse into one curve for small as in [88, Figure 3]. Indeed, for small , p= is independent of as can be readily veried from (5.5) and, thus, the sensory output p scales linearly with . The collapsed curve shows a remarkable resemblance to the average canal density itself, shown in Figure 5.1(b) and overlaid on the stimulation curves in Figure 5.3(d) for direct comparison. This correspondence conrms the ndings of [88] that the lateral line canals are concentrated at locations of strong hydrodynamic signals p. We now go beyond the results of [88] to explore the eect of the swimming speed V on p, as shown in Figure 5.3(d) for four dierent speeds V = 5; 10; 20 and 40. Clearly, p follows a similar trend for all V such that the curves p=V (for non-zero V ) collapse onto one curve for large V . The fact that for large V , p=V is independent of V follows immediately from (5.5). It implies that for swimming speeds V much larger than the speed of the background ow (here normalized to 1), the sensory output p scales linearly with V . Further, the value of V does not alter the sensors location for which p is maximum. To compute analytically the optimal sensor location opt for which p takes its maximum value, we set @p=@ = 0 to nd that opt = 1 2 cos 1 2R 2 R 4 + 1 : (5.6) CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 52 -2 -1 0 1 2 x 1 y 0 -1 (a) V< V cr V> V cr -2 -1 0 1 2 x 1 y 0 -1 (b) U U = 1 Figure 5.5: Trajectories of swimmer in uniform ow for (a) V = 1:44 and (b) V = 11:5. Initial orientations are (0) ==6;=3;=2; 2=3; 5=6. Parameter values are b=a = 0:3, Vcr = 4:33, and G = 1. It follows immediately from the separability property in (5.4) that opt is independent of V and and depends only on R = p (a +b)=(ab). Figure 5.4(a) presents a depiction of the pressure dierence p as a function of body orientation , swimming speedV and sensor location. This pictorial depiction conrms that the optimal sen- sor location opt for which p is maximum/minimum is independent ofV and. Meanwhile, there exist orientations for which the sensory output is identically zero irrespective of sensor placement . These blind spots, orientations for which the sensory output is zero, are highlighted in black in Figure 5.4(a) and calculated analytically by setting p = 0 in (5.5) and solving for , which yields = 0; = cos 1 V R 2 1 2R 2 ; =; (mod 2): (5.7) That is, p is zero for = 0 and when the body is aligned opposite to or in the same direction as the oncoming uniform ow. A third blind spot exists only when1V (R 2 1)=2R 2 1, that is, when the swimming speed V is below a critical value V cr = 2R 2 =(R 2 1). 5.3 Response to perceived pressure dierence We assume that the shlike body swims at a constant speed V in the direction of its major axis of symmetry, as shown in Figure 5.2(a), and responds to sensory information by adjusting its orientation depending on the value of the pressure dierence p. Namely, we describe the body response to sensory output using the following kinematic model _ x c =V cos; _ y c =V sin; _ =Gp: (5.8) Here,G is a constant `gain' parameter. The pressure output p depends on as well as the location of the pressure sensors and the swimming speedV . Equations (5.8) can be thus viewed as a model of the neuro-mechanical control or response of the sh to sensory information p. We examine the behavior when the body changes its orientation in response to the local pressure dierence p according to (7.1). Note that p and, consequently, the right-hand side of (7.1), depends explicitly on , but not on the position (x c ;y c ) of the body. This fact is due to the translational symmetry of motion in uniform ows. One can thus `reduce' equations (7.1) from three equations for (x c ;y c ;) to one equation in terms of only in the sense that one can solve CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 53 stable stable stable stable un- stable unstable unstable stable unstable unstable (a) (b) (c) Figure 5.6: Basins of attraction: (a) ForV = 0, the set of all initial orientations are equally split into a white subset that tends to align in the direction opposite to the incoming ow, and a grey subset that tends to align in the same direction as the incoming ow. (b) For 0<V <Vcr, the grey region shrinks until atV =Vcr, a pitchfork bifurcation occurs making = unstable. (c) For V > Vcr, all trajectories, except those parallel to , tend to align in the direction opposite to the ow. the orientation equation _ =Gp separately. Once the orientation is known, the trajectory in the (x c ;y c ) plane can be readily reconstructed. Orientations for which p = 0 are xed points or, more precisely, relative equilibria of (7.1). Along these directions, the body moves with constant speedV without changing orientation. These directions are linearly stable when @p=@ < 0 and unstable otherwise. A summary of the three families of xed points and their stability type is given in Figure 5.4(b). For swimmers moving at speeds V < V cr , both = 0 and = are stable, whereas = cos 1 (V=V cr ) is unstable. At V =V cr , the system undergoes a subcritical pitchfork bifurcation where the two unstable branches of = cos 1 (V=V cr ) collide with the stable branch =, making it unstable. For V >V cr , one has only one stable solution = 0 where the body aligns in the opposite direction to the oncoming ow. In other words, for low swimming speeds V relative to the normalized background ow, the swimmer orients either with or against the ow. As the swimming speed increases, the swimmer orients only in the opposite direction to the ow. Figure 5.5 depicts representative trajectories for both V < V cr and V > V cr . These results are consistent with observations of sh rheotaxis where sh are seen to respond to water motion either by eeing high- ow regions or by aligning their bodies in the direction opposite to the oncoming ow for upsteam migration and position holding [2, 73]. As time increases, the linearly stable direction = 0 is an attractor for all V and = is an attractor for V <V cr . For each stable direction , the basin of attraction is dened as the set of all initial values of for which the body eventually aligns with . Figure 5.6 depicts the basins of attraction for three representative cases corresponding to (a) zero swimming speed V = 0, (b) swimming at subcritical speeds V <V cr , and (c) swimming at supercritical speeds V >V cr . In (a) and (b), the set of all orientations is divided into two regions: the white region denes the basin of attraction for the stable orientation = 0, where the body aligns in the opposite direction to the ow, and the grey region denes the basin of attraction for = for which the body aligns with the ow. For V > V cr as in (c), the latter orientation becomes unstable and all initial conditions tend to align in the direction opposite to the oncoming ow. CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 54 5.4 Discussion The lateral line sensory system, which contains canal receptors that measure hydrodynamic pres- sure, is thought to be responsible for the ability of sh to sense minute water motions in their vicinity, and, thus, to play an important role in the sh response to ambient ows such as sh rheotaxis. Recent ndings based on sh anatomical measurements and ow experiments suggest that the arrangement of the lateral line system is correlated with locations along the sh body that experience strong variations in the hydrodynamic signal [88]. In this chapter, we presented a mathematical model consisting of an elliptic body equipped with laterally-arranged pressure sensors that emulate the sh lateral line sensory system. We studied the sensory output (the dierence in hydrodynamic pressure) in oncoming uniform ows as a function of the sensors location, swimming speed and body orientation. We found that the optimal placement of sensors that maximizes the sensory output is independent of the swimming speed and body orientation and depends only on the body geometry. Our results are consistent with the experimental ndings of [88] and support a physical and hydrodynamic explanation of the highly conserved architecture of the lateral line sensory system across many sh species { the sensory system extends over much of the body but is concentrated near the head. Most sh species have a body plan characterized by high curvature at the head. This body plan leads to a pressure prole that is greatest at the head, making it particularly sensitive to variations in the hydrodynamic signal, such as those due to changes in body orientation or swimming speed. For small deviations from the oncoming ow directions, the sensory output varies linearly with the body orientation. This linear dependence holds well up to = 20 . The pressure output also changes linearly with the swimming speed V when V is much larger than the background ow U = 1, namely, when V=U > 40. Such linear encoding of changes in sensory signal as the sh swimming parameters change could greatly simplify the neural processing involved in signal decoding. However, this linear dependence does not hold for nite body orientation and relatively moderate swimming speeds, implying that more complex computations are needed to decipher the sensory information in these regimes. We proposed that the shlike body responds to these sensory cues by adjusting its orientation proportionally to the perceived dierence in pressure. This modeling approach can be viewed as an abstraction, or `meta-model', of the coupling between the sensory output and the neural and muscular drive underlying sh behavior. Our working hypothesis is that these intertwined sensory, neural and actuation mechanisms produce an emergent rule or law that underlies sh rheotaxis. The behavior obtained from our proposed model is remarkably similar to the sh rheotactic behavior. At relatively low swimming speeds, sh go with the ow if they are initially oriented in a way that they cannot hold out against and turn into the oncoming ows. Otherwise, they align into the oncoming ows. This modeling framework couples the sensory and response mechanisms underlying sh rheo- taxis, albeit in the context of a potential ow model.It establishes the response law to sensory signals, regardless of the physiological mechanisms that bring this about. Here, a few comments on the model limitations are in order: (i) We considered a planar model. But an ellipsoidal body in uniform potential ows leads to similar results on optimal sensor location, equilibrium orientations, and stability. In uniform ows, the ow topology and pressure distribution around the centerline of an ellipsoidal body are similar to those around an ellipse { this is the reason why our results agree with the experimental ndings of [88]. CHAPTER 5. FLOW SENSING IN FISH-INSPIRED BEHAVIOR 55 (ii) We ignored the eects of viscosity and separation. In the experiments of [88] where the sh model is placed in uniform ows, separation occurs at the trailing edge of the sh { thus it does not aect the optimal location of sensors. Viscosity and separation will play prominent roles in non-uniform ows and complex sh maneuvers. Under these conditions, sh need multiple sensory outputs from their distributed lateral line to decipher the location of ow separation. (iii) The control law in (7.1) is linear. This is in part due to the lack of experimental evidence to support a particular choice of linear versus nonlinear control law. More importantly, a linear control law simplies the neural processing involved in responding to a given sensory output. Therefore, it is reasonable to conjecture that a linear law is preferable. We showed that the linear model leads the proper rheoactic behavior. More complex sh behavior might require more complex control laws. (iv) The model in (7.1) is kinematic. That is to say, it does not explicitly account for forces and moments, although they can be calculated posteriori for each trajectory using balance of momenta. Kinematic models are common in studying emergent behavior such as in sh schooling [22] where details of the coupling among the sensory feedback, neural control, and muscular actuation is not fully understood. Our goal here was to test whether the rule `turn in response to sensory output' reproduces shlike rheotaxis. We will link this rule-based model to the mechanics of motion in ambient ows in future studies. The two main results obtained from the simplied model { the optimal sensor placement that maximizes the sensory output and how rheotactic behavior emerges in light of this sensory infor- mation { serve to direct future research. We focused here on single sensory output and rheotactic behavior in uniform ows. But future work will address the interplay of multiple sensory outputs and how such sensory information could be used to avoid collision with nearby objects, detect ambi- ent vorticity, and inform changes in the sh body undulations [33, 86, 123]. Future models will also account for three-dimensional and viscous eects [49, 29, 37], which we expect to play a prominent role in non-planar sh motions and maneuvers as well as in non-uniform and high-speed ows, where ow separation is paramount. These models could also serve as a design tool for developing optimal sensory layouts and deployment strategies of underwater vehicles [33, 106, 112, 124]. The development of autonomous vehicles that detect hydrodynamic cues is important for their deploy- ment in unknown underwater environments and is relevant for many environmental, military, and industrial applications including ocean exploration, navigation, surveillance and monitoring. Chapter 6 Physics-Based Source Tracking Using Flow Sensing Several algorithms have been proposed for driving a mobile \agent" towards the source of a chemical signal. Algorithms that follow the gradient of the signal are particularly convenient when the signal itself is suciently high and smooth. In many situations, signal detection algorithms need to be combined with gradient-based algorithms; see, for example, [53, 32]. The group of Naomi Leonard used gradient-based methods to cooperatively control a network of mobile sensors to climb or descend an environmental gradient; see, for example, [3, 80, 7, 75]. Another variant of the gradient- based approach was developed in the group of Miroslav Krstic, by combining a nonholonomic mobile sensor with an extremum-seeking controller that drives the sensor towards the source of a static signal eld; see, for example, [14, 13]. Time-dependent signals, advected by a background ow, but in one dimension, were considered in [47]. In [113], an ingenious `infotaxis' approach was developed that is particularly suited to applications with time-dependent, spatially-sparse signals. The signal generated by the source is thought of as a message sent from the source and received by the mobile sensor through a noisy communication channel. The sensor reconstructs the message using a Bayesian approach conditioned on prior information gathered by the sensor. The sensor then maneuvers in a manner that seeks to balance its need to gain more information (exploration) versus its desire to drive towards the source (exploitation). In this study, we extend the notion of source seeking to problems where the relevant signal is contained in the ow eld itself, as opposed to a chemical scalar eld that is carried by the ow. The information content of various ows (e.g., shear rate, vorticity eld, vortex wake patterns, etc.) and how to extract such information from local sensory measurements of the pressure and velocity elds was explored in a series of papers. Specically, in [17, 19], we used a physics-based approach to decipher the local ow character (extensional, rotational, or shear ow) from velocity sensors, while in [15], we trained a neural network to identify the global ow type (wake pattern of an oscillating airfoil) from local sensory measurements. Here, we develop a local sensory and control approach that is designed to operate in time-periodic signals of \traveling wave"-like character. This includes the wakes generated by stationary and self-propelled bodies at intermediate Reynolds numbers. We posit that a mobile sensor moving in the opposite direction to the direction of propagation of the signal will locate its generating source. We measure the signal strength locally over one period of signal oscillation and transform its time history to the frequency domain to determine the direction 56 CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 57 signal approximated as traveling wave agent source (a) (b) Figure 6.1: (a) A mobile sensor, located at xs and equipped with body-xed frame (b 1 ;b 2 ), moves at constant speedV and changes its orientation with angular velocity in response to local sensory measurements of the signal eld f(x;t). Near the sensor, f(x;t) is expanded into a series of traveling waves with modal amplitudes an and bn, frequencies wn, wavenumberskknk and directions kn=kknk indicating the local direction of signal propagation. (b) Relation between the inertial frame (e 1 ;e 2 ), body-xed frame (b 1 ;b 2 ), whereb 1 points in the heading direction of the sensor, and body-xed frame (c 1 ;c 2 ), wherec 1 is in the direction pointing from the sensor to the source. The pointing or alignment error n is the angle between c 1 andrn=krnk. of propagation of the signal and its intensity. The mobile sensor then changes its orientation in response to this sensory information so as to orient opposite to the signal propagation. The coupling between the response of the mobile sensor and the signal eld makes the system amenable to analysis only for relatively simple signals. Employing radially-symmetric signal elds and using techniques from local stability analysis and nonlinear dynamical systems theory, we rigorously prove unconditional convergence of the mobile sensor to an \orbit-like" attractor around the source under certain assumptions on the control gain. To assess the behavior of the sensor in more realistic ows, we conduct careful numerical simulations of uid ow past an oscillating airfoil and we probe the response of the mobile sensor in these more complex ow elds. 6.1 Sensory system and signal eld model We consider a mobile underwater sensor operating in an unknown hydrodynamic environment. The mobile sensor is tasked with locating the source of a periodic ow disturbance in an unbounded two-dimensional domain. The sensor measures the local signal contained in the uid ow, then adjusts its heading accordingly. Our goal is to design a controller that will drive the mobile sensor to the source of the ow. Mobile sensor We model the mobile sensor as a nonholonomic \agent" with a sensor located at its center. The sensor moves at a constant translational speed V but is capable of changing its heading angle in response to its local sensory measurement (see Figure 6.1). The equations of motion of the sensor are given by _ x s =V cos e 1 +V sin e 2 ; _ = (s); (6.1) where x s is the position vector of the sensor with respect to an inertial frame (e 1 ; e 2 ) and is measured from the positive e 1 direction. The angular velocity changes according to the local sensory measurement s. CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 58 Kinematics of Mobile Sensor Time History of Signal Field Spectral Transform Signal Processing Figure 6.2: Block diagram depicting the sensory and control laws of the mobile sensor. The time history of the signal eld f(xs;t) is measured locally at xs over one oscillation period of the signal and transformed to the frequency domain ^ f(xs;wn). The magnitude mn(x) and phase elds n(x) are computed from the spectral information to obtain the sensory output s and determine the gain function G of the feedback controller =Gs. Signal eld We consider an oscillating source located at the origin of the inertial frame. The source generates a hydrodynamic signal in the form of a scalar eld f(x;t) that propagates away from the source. As the source oscillates with period T , we assume that the signal eld reaches a steady-state limit cycle such that f(x;t) = f(x;tT ) is also time periodic with period T . This admits a broad class of signal elds including vortex wakes generated by biological and bioinspired engineered devices. The exact shape of f(x;t) is not known to the sensor. Our goal is to steer the mobile sensor towards the origin using only time history of the signal at the sensor location. In light of the assumed structure of the signal eld, we can locally (in the vicinity of the sensor at x s ) decompose f(x;t) as a series of traveling waves f(x;t) = 1 X n=0 a n (x) cos (k n (x) (x x s )w n t) +b n (x) sin (k n (x) (x x s )w n t); (6.2) where n is a nonnegative integer indicating the mode index, a n (x) and b n (x) are the amplitudes, w n are the angular frequencies, and k n (x) are the wavenumber vectors of the n th modes. Next, we establish a framework that allows the mobile sensor to adjust its heading direction relative to the direction of propagation of the signal in order to locate the source position. This framework relies on the premises that the signal eld propagates in a traveling wave fashion and that the source- seeking agent is capable of deciphering salient information contained in these traveling waves, as depicted in Figure 6.2 and explained next. Local ow sensing We dene the frequency spectrum ^ f(x;w) of the signal eldf(x;t) as follows ^ f(x;w) = lim !1 !w 1 2 Z f(x;t) exp(it)dt; (6.3) where i = p 1 is the imaginary unit and w is the frequency at which the spectrum is be- ing evaluated. We substitute (6.2) into (6.3), recalling Euler's trigonometric identities cosu = CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 59 [exp(iu) + exp(iu)]=2 and sinu = [exp(iu) exp(iu)]=2i. We get that at w =w n , ^ f(x;w n ) = 1 2 1 X n=0 (a n + ib n ) exp(ik n (x x s )): (6.4) Ecient and practical estimation of the relevant frequencies w n contained within a signal is the subject of a great deal of established literature (see, e.g., [101]). Here, we remain agnostic to the method of spectral analysis and we assume that the period T of the signalf(x;t) can be estimated oine and that w n = 2n=T . For each mode n, we dene the spectral magnitude eld m n (x) and the spectral phase eld n (x), m n (x) =k ^ f(x;w n )k and n (x) =\ ^ f(x;w n ); (6.5) wherekk and \() denote the magnitude and argument of a complex number, respectively. In light of (6.4), these elds are given by m n (x) = 1 2 p a 2 n (x) +b 2 n (x); (6.6) and n (x) = tan 1 b n (x) a n (x) k n (x) (x x s ) mod 2: (6.7) We expand the phase eld (x) about x s using Taylor series expansion, n (x) = n (x s ) +r n (x) xs (x x s ) +O(kx x s k 2 ); (6.8) wherer() xs denotes the spatial gradient evaluated at x s . Comparing (6.7) and (6.8) and neglect- ing higher order terms, we obtain that n (x s ) tan 1 b n a n and r n (x s )k n : (6.9) The wavenumber vector k n describes the direction of propagation of the signal in (6.2). Therefore, if the phase eld n and its gradientr n can be measured, the latter can be used as a tool to determine the direction of signal propagation. With this in mind, we dene the spectral direction eldr n =kr n k, which is a unit vector eld describing the local direction of signal propagation for the n th mode. A few remarks regarding the practical computation of the phase gradient are in order. The assumption of direct access to the gradient is feasible because one can use multiple sensors on the sensory platform to produce a nite-dierence approximation of the gradient. A more subtle issue lies in the fact that n (x) is a scalar eld in the frequency domain whereas our mobile sensor lives in the time domain and needs to respond to the measured signal in real time. Measuring r n (x s ) requires access to the time history of the f(x;t) over one period of oscillations. Here, we assume that the mobile sensor moves at a much slower time scale compared to the time scale of the signal oscillations. This translates to a constraint on the speed V of the mobile sensor such that V w n =kk n k where w n =kk n k is the wavespeed of mode n. In this regime, we can use the quasi- steady assumption that the sensor measures a full time-period of the signal eld \instantaneously" relative to its swimming speed. CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 60 Sensory information and control law Since the rst mode n = 1 in the ambient signal is likely to be among the most energetic modes and therefore the most easily measured with accuracy in practical applications, we assume that it is sucient for the sensor to measure the rst mode only. Hereafter, we drop the index n with the understanding that n = 1. We ax a body frame (b 1 ; b 2 ) to the agent such that b 1 is pointing in the heading direction (b 1 e 1 = cos as shown in Figure 6.1). We consider the measurable sensory signal s to be the lateral component of the phase gradient at the position x s of the sensor, namely, s = r(x s ) kr(x s )k b 2 (): (6.10) We introduce the feedback controller = G(m(x s ))s(x s ;), where G(m(x s )) is a proportional gain function that may dependent on the spectral magnitude eld m(x s ) evaluated at the sensor location. The objective of this controller is to steer the mobile sensor to follow the local gradient of the spectral phase eld. The details of both the sensory system and control model are summarized schematically in Figure 6.2. Equations of motion for sensor platform We substitute the control law =G(m(x s ))s(x s ;) into equations (6.1) to obtain a closed system of equations that govern the motion of the mobile sensor in response to the local sensory signal s it perceives, _ x s =V b 1 (); _ =G r(x s ) kr(x s )k b 2 (): (6.11) In component form, one gets three scalar, rst-order dierential equations, _ x =V cos; _ y =V sin; _ =G 2 6 6 6 6 4 @ @x sin + @ @y cos s @ @x 2 + @ @y 2 3 7 7 7 7 5 : (6.12) Here, the spectral magnitude m and phase of the signal (as well as @=@x and @=@y) are all evaluated at the location x s (x;y) of the sensor. For a given choice of G(m), these coupled equations, together with the local sensory measurements of the signal, constitute a closed-form system that can be solved to obtain the trajectory of the mobile sensor (x(t);y(t)) and its heading angle (t). To analyze the ability of the mobile sensor to track the signal eld and locate its source, it is convenient to introduce a right-handed orthonormal frame (c 1 ; c 2 ) that is attached to the sensor and pointing towards the source such that c 1 =x=kxk. It is also convenient to dene the alignment error (x) as the angle measured from c 1 (the unit vector pointing from the sensor directly to the source) to the direction of propagation of the signalr(x)=kr(x)k such that r krk = cos c 1 + sin c 2 : (6.13) The alignment error eld(x) is a local measure of the degree of alignment of the spectral direction eldr=krk with the direction towards the source. When = 0,r points directly toward the CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 61 source. For < 0 or > 0,r points to the left or to the right of the source, respectively. Since we are driving the mobile sensor to follow the local spectral direction eld under the premise that it will lead to the source, the alignment error provides a measure of how well the mobile sensor will head towards the source from a given location x in the signal eld. We now rewrite the equations of motion (7.2) in the (c 1 ; c 2 ) frame. To this end, note that (c 1 ; c 2 ) is related to (b 1 ; b 2 ) via the orthogonal transformation c 1 c 2 = cos sin sin cos b 1 b 2 ; (6.14) and to the inertial frame (e 1 ; e 2 ) via the transformation c 1 c 2 = cos sin sin cos e 1 e 2 : (6.15) Here, the angle is given by the polar transformation from (x;y) to (r;) x =r cos; y =r sin; (6.16) with r 2 =x 2 +y 2 . The angle is given by = (). In Figure 6.1(b) is a depiction of the all three reference frames (e 1 ; e 2 ), (b 1 ; b 2 ) and (c 1 ; c 2 ), together with the polar coordinate angle and the angles and . It is important to emphasize here that does not contain any information about the heading direction of the sensor but does. By denition of (r;) and (c 1 ; c 2 ), we get that x s =rc 1 and _ x s = _ rc 1 r _ c 2 . Substituting these expressions into (6.11) and using (6.13) and (6.14), we get that _ r =V cos ; r _ =V sin ; r _ =V sin rG (cos sin + sin cos ): (6.17) Here, the alignment error (r;) is expressed as a function of the polar coordinates (r;). This is rooted in the fact that the signal eld f, and consequently its magnitude m and phase , are functions of space only, which can be expressed as functions of (r;) only. By choice, we will consider three types of gain functions: (i) static gain of constant valueG o , (ii) gain proportional to the magnitude of the signalG =G o m(r;), and (iii) inversely-proportional gain G =G o =m(r;). It will be convenient for later to introduce the constant dimensionless parameter =V=G o , which can be thought of as a gain length scale. In all three cases, the gainG is function of the sensor position only. The coupled equations in (6.17) together with the local sensory measurement of the signal eld and the choice of the gain function constitute a closed-form system that describe the motion of the sensor in terms of (r;; ). They are equivalent to equations (6.12), except at the origin. 6.2 Analytical performance of sensor in nding source Canonical signal eld We consider a simplied signal eld that locally approximates the signal generated by an oscillating source. In particular, we consider a scalar signal eld f(r;;t)f(r;t) propagating in a radially symmetric fashion (independent of ) from the origin, f(r;t) = 2 exp r ` cos (krwt): (6.18) CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 62 Figure 6.3: Radially-symmetric signal eld in (6.19) for` = 6:5. (a) Snapshot of the signal eldf(x;t) at timet = 0, (b) the magnitude eldm(x), (c) the phase eld(x), (d) the direction eldr. Note that the associated alignment error (x) is identically zero. Here, , k and w are the signal magnitude, wavenumber and frequency, respectively, and ` is a length scale that determines the spatial decay of the signal strength. To simplify our analysis we use the time scale 1=w, the length scale 1=k, and the signal intensity scale to dene the dimensionless quantities x y =xk; y y =yk; r y =rk; t y =tw; ` y =`k; f y = f : For clarity, we drop the () y notation and rewrite (6.18) in nondimensional form, f(r;t) = 2 exp r ` cos(rt): (6.19) A snapshot of this signal is shown in Figure 6.3(a). Using (6.3), the spectrum of (6.19) evaluated at w = 1 is given by ^ f(r) = exp r ` ir : (6.20) Substituting into (6.5), one gets the associated magnitude and phase elds, m(r) = exp r ` and (r) =r: (6.21) These elds are shown in Figures 6.3(b) and (c), respectively. The normalized gradient of (r) is given byr=krk = c 1 and the alignment error is identically zero, = 0, everywhere in the (r;)-plane, meaning that the spectral direction eldr always points towards the source; the normalized gradient eld and contours of the alignment error are depicted in Figure 6.3(d). CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 63 We substitute = 0 into (6.17) to get the equations of motion of the sensor in the simplied signal eld f(r;t), _ r =V cos ; r _ =V sin ; r _ = (VrG(r)) sin ; (6.22) with the understanding that V and G are dimensionless. Due to the rotational symmetry of the signal eld, the equations of motion are invariant under superimposed rotations on the eld-sensor system. As a result, the equations for r and in (6.22) decouple from the equation for and can be solved independently. The results can then be substituted into the equation r _ = V sin to reconstruct (t) = R t 0 dV sin( ())=r(). The two coupled equations in (r; ) can be reduced further by noting that the phase gradientr associated with the signal eld in (6.19) is invariant under translations from r to rt, which gives rise to an integral of motionQ. To computeQ, it is convenient to rewrite the equations for r and in dierential form, r VrG d sin = dr V cos ; (6.23) which after straightforward manipulation yields the exact dierential, dQ = d(sin ) sin + 1 r G V dr = 0: (6.24) A direct integration of the above dierential givesQ = r sin exp R G(r)dr=V = constant. Without loss of generality, we can scaleQ by the constant 1= =G o =V to get Q(r; ) = r sin exp Z G(r) V dr = constant: (6.25) We solve for in terms ofQ and r and substitute back into the equation for r in (6.22) to get _ r =V v u u u t 1 Q 2 r 2 2 exp 2 V R G(r)dr : (6.26) This equation can be readily integrated to get r(t). Next, we analyze the closed-loop performance of the mobile sensor in the simplied signal eld in the context of three types of gain functions: (i) static gain of constant valueG o , (ii) gain proportional to the magnitude of the signal G = G o m(r), and (iii) inversely-proportional gain G = G o =m(r). In each case, we show the behavior of the sensor in the (x;y) plane by numerically integrating the nonlinear equations of motion expressed in cartesian coordinates in (6.12). We then analyze the behavior in the phase space (r cos ;r sin ) associated with the system of equations expressed in (6.22) after eliminating the equation. Lastly, we determine the success of the sensory-control law in following the signal eld and converging towards its source. Static gain Consider the case of a static gain and setG to a constant valueG o . We solve the sys- tem of nonlinear equations in (6.12) numerically forV = 1 andG o = 0:5. Figure 6.4(a) (left) shows representative trajectories in the (x;y) plane for distinct initial conditions, at progressively larger distances from the source, marked as `' on the plot. In each case, the mobile sensor approaches and loops around the target, providing multiple opportunities to intercept the source. CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 64 Figure 6.4: Behavior of the mobile sensor in the radially-symmetric signal eld in (6.19). The left panel depicts trajectories in the (x;y) plane starting at positions located successively further away from the source. In the middle panel, trajectories are shown in the phase space (r cos ;r sin ), with xed points depicted as black dots. Relative equilibria along sin = 0 are depicted as solid black lines. The right panel depicts the solution curves in the (r; _ r) space for various values of the integral of motionQ. For all cases V = 1 and = 2:0. (a) Static gain. Two center- type xed points occur. Unconditional convergence is achieved for all initial conditions. (b) Proportional gain with ` = 6:5. Two center-type and two saddle-type xed points occur. The homoclinic orbits associated with the saddle xed points are highlighted in black. Conditional convergence is achieved for initial conditions inside the closed-loop of these orbits. (c) Proportional gain with` = 5:4. No xed points occur and no convergence. (d) Inverse gain with ` = 5:8. Two center-type xed points occur. Unconditional convergence is achieved for all initial conditions. CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 65 We analyze this behavior in the phase space (r; ) where, from (6.22), one has _ r =V cos ; _ = 1 r 1 V sin : (6.27) This system of equation admits two xed points (r ; ) for which _ r = _ = 0 for all time, r =; = 2 : (6.28) These xed points are marked as `' in Figure 6.4(a) (middle) in the phase space (r cos;r sin ). We determine the linear stability of these equilibria by computing the Jacobian matrix J of (6.27) and evaluating it at (r ; ) J ; 2 =V 0 1 2 0 : (6.29) The corresponding eigenvalues =iG o are purely imaginary; thus, the xed points are centers by nature and are marginally stable. We also observe a line of degenerate relative equilibria along sin = 0, for which = 0 mod() and _ = 0 and _ r =V for all time. Along this line, the sensor heads directly towards the source. We mark this line in Figure 6.4(a) (middle) and note that it separates the sensor trajectories that orbit clockwise around the source from those that orbit in a counter-clockwise fashion. To complete the phase portrait of the mobile agent, we superimpose onto Figure 6.4(a) (middle) the contours of the integral of motionQ. For static gain, one has Q = r sin exp r = constant: (6.30) Generic contour lines are shown in gray. The specic contours corresponding to the trajectories shown in Figure 6.4(a) (left) are highlighted in color. Lastly, in Figure 6.4)(a) (right), we plot the radial phase velocity _ r as a function of r for the ve dierent values ofQ corresponding to the trajectories in 6.4a (left). Here, we use the radial dynamics in (6.26), which for static gain becomes _ r =V v u u u t 1 Q 2 r 2 2 exp 2r : (6.31) In each case we mark the initial condition and phase velocity by `'. We also label each curve by its corresponding value ofQ. Clearly, for each value ofQ, the radial excursion is bounded from below and above. To prove this for all initial conditions, we note that physically-speaking, since r is a positive scalar, equation (6.31) is well-dened only when the right-hand side is real-valued, that is to say, Q r exp r 1: (6.32) This condition can be used to determine the upper and lower bounds on the distance r, W 0 (jQj) r W 1 (jQj); (6.33) CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 66 whereW p () is thep th Lambert product logarithm function. Thus, we have shown thatr is bounded from both below and above for all values ofQ; in other words, r remains nite for all time and all initial conditions. The implication of this result is that regardless of the initial conditions, the trajectory of the mobile sensor will circle around the source within a bounded distance r from the source. These periodic loops around the source have the practical advantage that the sensor will have multiple opportunities to approach the source. We say that the mobile sensor exhibit unconditional convergence to the source. Proportional gain We now consider a gain strategy that is proportional to the magnitude of the signal eld G =G o exp(r=`). Figure 6.4(b) (left) shows the behavior of the mobile sensor in the (x;y) plane for V = 1, G o = 0:5, and ` = 6:5 and same initial conditions as in Figure 6.4(a). Unlike in Figure 6.4(a) (left), here, three trajectories (shown in green, purple, and yellow) get closer to the source once and then travel o to innity, whereas two trajectories (shown in red and blue) make multiple loops around the source. This implies a change in the stability characteristics of the system as the gain changes, with some initial conditions resulting in unstable trajectories, drifting o tor!1 as time progresses. To understand this change in behavior, we examine the dynamics in the phase space (r cos ;r sin ), given by the equations of motion _ r =V cos ; _ = 1 r 1 exp r ` V sin : (6.34) Here, _ is inversely proportional to an exponential function of r. The physical interpretation is that, for this proportional gain, the mobile sensor will turn more sharply as it gets closer to the source. Equations (6.34) admit four xed points for which _ r = _ = 0 for all time, r =`W 1 ` ; = 2 ; (6.35) and r =`W 0 ` ; = 2 : (6.36) These xed points are marked as `' in Figure 6.4(b) (middle). To examine the linear stability of these xed points, we evaluate the Jacobians J =V 0 1 X 0 ; (6.37) where Xj (r ; ) = 1 ` 2 W 1 ` 0 @ 1 W 1 ` + 1 1 A ; (6.38) and Xj (r ; ) = 1 ` 2 W 0 ` 0 @ 1 W 0 ` + 1 1 A : (6.39) CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 67 This leads to the eigenvalue pairs ` p V = v u u u t 1 W 1 ( ` ) 0 @ 1 W 1 ( ` ) + 1 1 A ; (6.40) and ` p V =i v u u u t 1 W 0 ( ` ) 0 @ 1 W 0 ( ` ) + 1 1 A ; (6.41) for (r ; ) and (r ; ), respectively. The xed points (r ; ) are unstable (saddle type) while (r ; ) are linearly stable (centers). As before, = 0 mod () corresponds to a degenerate family of relative equilibria. We complete the phase portrait in Figure 6.4(b) (middle) by plotting the contour lines ofQ given by Q = r sin exp ` exp r : (6.42) Notice the homoclinic-type separatrices, highlighted with a thickened black curve, associated with the saddle xed points. Trajectories that start inside the homoclinic orbit are stable and those outside are unstable. This distinction explains the dierence in behavior we observed in Fig- ure 6.4(b) (left). In Figure 6.4(b) (right), we plot the radial phase velocity _ r as a function of r for the ve dierent values ofQ corresponding to the trajectories in 6.4(b) (left). Here, we use the radial dynamics in (6.26), which for proportional gain becomes _ r =V v u u u t 1 Q 2 r 2 2 exp 2 ` exp r ` : (6.43) We dene the critical conserved quantityQ cr as the value ofQ corresponding to the separatrices associated with the saddle points (r ;=2). That is to say, we evaluateQ at (r ;=2) to get Q cr = ` W 1 ` exp 2 4 1 W 1 ` 3 5 : (6.44) IfjQj <jQ cr j, the sensor will approach the source only once, as seen in the green, purple, and yellow curves. On the other hand, ifjQj>jQ cr j, the trajectory will either loop around the source indenitely or will diverge to innity depending on whether the initial radial position is greater than or less than r , respectively. Given that the radial excursion is bounded from above for only a subset of the initial conditions, we will call this conditional convergence. Let us consider the case where ` = 5:4. Keeping V and G o the same, we solve the system of nonlinear equations in (6.1) numerically, generating the trajectories in the (x;y) plane which we plot in Figure 6.4(c) (left). Clearly, all of the the trajectories make one close-in pass of the source and then move away to innity. CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 68 To understand this behavior, we reexamine (6.34) with the lower value of` and nd that there are no xed points. This is because a saddle-node bifuraction occurs when ` =e and the saddle- and center-type xed points collide and annihilate. Upon further examination of Figure 6.4(c) (middle), we conrm that all trajectories regardless of initial condition are unbounded. Finally, the radial dynamics in Figure 6.4(c) (right) conrm thatr is unbounded for all initial conditions. The sensor is said to diverge. Inversely proportional gain Lastly, we consider the gain strategy G = G o exp(r=`) that is inversely proportional to the magnitude eld. That is, the gain decreases with signal strength. In Figure 6.4(d) (left), we plot the motion of the sensor in the (x;y) plane for V = 1, G o = 0:5, and ` = 5:8. In this case, the trajectories resemble those in Figure 6.4(a) (left), but exhibit longer loops with tighter turns around the source. As before, we examine the dynamics expressed in the phase space (r cos ;r sin ), _ r =V cos ; _ = 1 r 1 exp r ` V sin : (6.45) Here, the rate of change of the heading angle _ is proportional to exp(r=`), implying that the mobile sensor turns more sharply when it is further away from the source. This is in contrast to the proportional gain strategy _ exp(r=`), where the gain increases with signal strength and the sensor turns more sharply closer to the source. Equations (6.45) admit two xed points for which _ r = _ = 0 for all time, r =`W 0 ` ; = 2 ; (6.46) which are depicted in Figure 6.4(d) (middle) as `'. To examine the linear stability of these xed points, we evaluate the Jacobian J r ; 2 =V 0 1 X 0 ; (6.47) where X = 1 ` 2 W 0 ` 0 @ 1 W 0 ` + 1 1 A : (6.48) This leads to the eigenvalue pairs ` p V =i v u u u t 1 W 0 ` 0 @ 1 W 0 ` + 1 1 A : (6.49) Since these eigenvalue pairs are purely imaginary and equal in magnitude, the xed point is a center and is linearly stable. We complete the phase portrait in Figure 6.4(d) (middle) by plotting the contour lines of the conserved quantity, Q = r sin exp ` exp r ` ; (6.50) CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 69 The contour lines are closed orbits around the xed points, in agreement with the fact that these points are of center type. Importantly, all orbits are periodic and bounded, meaning that the mobile sensor does not travel o towards r!1. The boundedness of the radial distance from the source r can be best seen by examining the radial equation of motion _ r =V v u u u t 1 Q 2 r 2 2 exp 2r ` ; (6.51) We plot the radial velocity _ r versus r in Figure 6.4(d) (right) for all ve dierent trajectories highlighted in Figure 6.4(d) (left). The radial excursion r is bounded from both below and above for all initial conditions and hence the sensor dynamics is said to be unconditionally convergent. 6.3 Locating a pitching airfoil by following its wake Now that we analyzed the convergence properties of the mobile sensor in a simplied signal eld, we seek to demonstrate its performance in a more complex signal eld arising from the ow past an oscillating airfoil. To this end, we consider an airfoil-shaped body with diameter D and chord C in an otherwise uniform ow of speed U. We control the pitch angle (t) sinusoidally with frequency f such that (t) = max cos(2ft), as depicted in Figure 6.5. The parameters can be reduced to three dimensionless parameters: the Strouhal number St =fD=U, the amplitude ratio A = 2C sin max =D, and the Reynolds number Re = UD=, where and are the uid density and viscosity. These dimensionless parameters characterize the structure of the resultant wake; see, e.g., [95, 15]. The ow around the airfoil is described by the uid velocity vector eld u(x;t) and the scalar pressure eld p(x;t). The time evolution of these elds is governed by the incompressible Navier- Stokes equations @u @t =ururp + 1 Re r 2 u; r u = 0: (6.52) together with the no-slip condition at the airfoil boundary. We solve these equations numerically for Re = 500 and St = 0:075 using an algorithm based on the Immersed Boundary Method [84]. The algorithm is described in detail in [72], and has been implemented, optimized and tested extensively by the group of Haibo Dong; see, e.g., [111, 8]. We calculate the vorticity eld!(x;t) = (ru)e 3 where e 3 = e 1 e 2 is the unit vector normal to the plane. Figure 6.5(a) shows contours of vorticity for a snapshot in time. We take !(x;t) as the signal eld. We compute its Fourier transform ^ !(x;w n ), evaluated at the rst Fourier mode n = 1, and calculate its magnitude m(x) and phase (x), as depicted in Figures 6.5(b) and (c), respectively. The direction eldr=krk and the alignment error(x) are plotted in Figure 6.5(d). Although the structure of the vorticity signal eld !(x;t) and its spectral elds in Figure 6.5 is far more complex than the radially-symmetric signal given in (6.19) and plotted in Figure 6.3, the two signals have some similarities. First, the magnitude m(x) of the signal increases as we approach the airfoil. Second, the direction eldr=krk points in the upstream direction over a large area downstream of the airfoil. These two components lead us to believe that the control law proposed in this chapter will drive a mobile sensor initially positioned downstream to follow the hydrodynamic signal and converge to the airfoil. In fact, downstream of the airfoil, the alignment error (x) is nonzero; to the left of the centerline, the alignment error is largely positive, and the CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 70 wake airfoil Figure 6.5: Behavior of the mobile sensor in the wake of an oscillating airfoil in a uniform incoming ow of strengthU. The airfoil oscillates at dimensionless frequency St (Strouhal number) and amplitude A. A wake forms downstream of the oscillating airfoil, with regions of spatiotemporally-correlated coherent vortical structures. (a) Snapshot of vorticity eld eld !(x;t) for Re = 500, St = 0:075 and A = 1:66. (b) the magnitude eld m(x) of the spectral vorticity, (c) the phase eld (x), (d) the direction eldr=krk and the associated alignment error (x). CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 71 vector eldr=krk points towards the centerline. To the left of the centerline, the alignment error is largely negative, andr=krk points towards the centerline. Therefore, a mobile sensor initially positioned in the wake of the airfoil and re-orienting according tor=krk will be directed towards the source. Following the analysis done for the simplied signal, we test this prediction for three types of gain functions: static, proportional and inversely-proportional. Static gain We solve (6.12) numerically using the signal from the vorticity eld for various initial conditions located downstream of the airfoil at successively larger distances from the source. Results are shown in Figure 6.6(a) in the (x;y) plane. The mobile sensor largely oscillates about the centerline until it reaches the airfoil. For the trajectories shown in green, red, and blue, the mobile sensor ends up very close to the airfoil. However, for the purple and yellow trajectories, the mobile sensor's nal position is further o the centerline. In particular, the yellow trajectory is an example of an overshoot, where the mobile sensor did not turn sharply enough to stay in the wake. Overall, the trajectories end up in the vicinity of the airfoil but do not concentrate along the centerline { instead, they are distributed transversally to the wake. This is due to the fact that the gain strategy does not take into consideration any information about the signal strength, which is largely concentrated along the centerline. Proportional gain We incorporate the signal strength by choosing a gain strategyG =G o m(x) proportional to the signal magnitude. We again solve (6.12) for the same initial conditions. The resulting trajectories in the (x;y) plane are plotted in Figure 6.6(b). For this choice of initial conditions, the mobile sensor approaches then tracks the centerline of the wake. The proportional gain strategy causes the mobile sensor to turn more sharply in response to higher signal magnitude. Since the magnitude of the signal is concentrated along the centerline of the wake, the mobile sensor stays more closely aligned to it once it reaches the centerline. Therefore, the trajectories track the airfoil and approach it very closely. Inversely-proportional gain Lastly, we consider the gain strategy whereG is inversely propor- tional to the signal magnitude such that G = G o =m(x). Figure 6.6(a) shows the trajectories of the sensor for the same initial conditions considered in Figures 6.6(a,b). Clearly, the mobile sensor oscillates about the centerline until it reaches the airfoil. The sensor tends to undergo more lateral excursions as it approaches the airfoil. In the regions of the wake where the magnitude eld is vanishingly small, the gain function becomes singular and the control law is ill-posed. Therefore, this algorithm is not practical for online implementation given that in real applications, there is no guarantee that the signal eld will be always non-zero. 6.4 Discussion We presented a novel framework for analyzing the information embedded in time-periodic ow elds in the context of a source-seeking problem, where a mobile sensor, with only local information, is required to track the ow and locate its generating source. Our approach takes advantage of the \traveling wave"-like character of the signal eld in order to reorient the mobile sensor in the direction opposite to the direction of propagation of the signal. Specically, we measure the direction of signal propagation in the frequency domain and map back to the time domain to CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 72 Figure 6.6: Trajectories of the mobile sensor in the wake of the oscillating airfoil for the parameter values shown in gure 6.5: (a) static gain law, (b) proportional gain law, and (c) inversely proportional gain law. In all cases, the parameters of the mobile sensor are set to Go = 1 and V = 1. CHAPTER 6. PHYSICS-BASED SOURCE TRACKING USING FLOW SENSING 73 reorient the sensor. This method is best suited for systems where the time scale associated with the signal is much smaller than the time scale associated with the sensor. We rigorously analyzed the convergence of the sensor to the source in the context of a simplied signal eld. Then, through carefully laid out simulations of uid ow past an oscillating airfoil, we demonstrated the ecacy of this sensory-control system in tracking the ow signal and locating the airfoil. A few comments on the limitations of the algorithm are in order. The transformation of time domain information at one position in space to the frequency domain requires in principle that the mobile sensor remains stationary for one period of the signal. However, waiting to collect the signal may not be feasible or practical to enforce in certain applications. In this case, we would need to extend the algorithm to perform true real-time spectral estimation, for which there are a variety of useful algorithms in the spectral analysis literature (see e.g. [1, 62]). Another limitation of the model is that the mobile sensor is one-way coupled to the signal eld | that is, it can probe the signal but it does not aect it. This condition was used to simplify the analysis and demonstrate the eectiveness of the sensing scheme. Further research is required to understand the coupling of ow sensing to self-generated ow disturbances as done in [36, 35, 127]. Our results, that the time history of local ow quantities contain information that can be used to track the ow and locate its source, can be seen as a rst step towards laying a foundation for source seeking in uid ows. Indeed, in [15], we demonstrated, using neural networks, that local time histories of the ow contain information relevant to categorizing the ow eld and hence its source, whereas here we showed, using a deterministic control law, that local time histories contain information relevant to nding the source. Future extensions of these works will unify these two strategies, by combining tools from articial intelligence and machine learning with rigorous dynamical systems and control theoretic methods, in a robust manner to allow mobile \agents" equipped with hydrodynamic sensing capabilities to identify and localize other objects or agents in the ow eld. Chapter 7 Data-Driven Source Tracking Using Visual Sensing In this chapter, we consider the problem of locating a target using a mobile sensor with limited `eld of view.' We employ a reinforcement learning approach to steer the vehicle to the target. Our control strategy keeps the forward velocity at a constant value and nds an optimal controller for the angular velocity, a setting applicable to most autonomous vehicles. Reinforcement learning methods (see, e.g., [103]) have been used very successfully in the gaming industry, but only recently began to make their way to physics-based and robotic applications. See, for example, [38, 78, 114]. Here, we develop a simple method to map a time-continuous dynamical system into a discrete Markov process. We then discuss three algorithms for nding optimal control strategies: a model- based algorithm that requires a priori knowledge of the dynamics and reward over the whole space of states and actions, and two model-free algorithms that learn from experience. We nd that in all algorithms, the vehicle learns a policy that allows it to reach the target with high success rates. This framework can be readily extended to many applications in autonomous vehicle control. 7.1 Mobile agent with visual sensing Consider an agent whose position is described by x = xe 1 +ye 2 where (e 1 ; e 2 ) is an inertial frame. The agent is tasked with nding a target located at the origin of this frame. Let denote the heading angle of the agent, measured counter-clockwise from e 1 . We introduce a body frame (b 1 ; b 2 ) attached to agent such that represents the angle from e 1 to b 1 as depicted in Fig 7.1. Equations of motion The motion of the agent is given by _ x =V cos; _ y =V sin; _ = ; (7.1) where V is a translational speed and is a rotation rate. The position of the target relative to the agent is dened via the relative position vector r =xe 1 ye 2 . The coordinates (x;y) can be expressed in body frame using the nonlinear transformationx =r cos( +) andy =r sin( +), where r is a positive distance (r 2 = x 2 +y 2 ) and is the bearing angle depicted in Fig 7.1. In body frame, one has r =r cos b 1 +r sin b 2 . The goal of the agent to drive r! 0. 74 CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 75 agent target field of view inset Figure 7.1: The agent located at x =xe 1 +ye 2 is tasked with nding a target located at the origin. A body frame (b 1 ;b 2 ) is axed to the agent and the agent's heading is denoted . The distance r and bearing are computed relative to the agent. The agent's motion is controlled by changing its forward speed V and its rotation rate . Equations (7.1) can be rewritten in body frame in the form _ r =V cos ; _ = V r sin : (7.2) Further, introducing r 1 =r cos and r 2 =r sin , we can rewrite (7.2) as a linear dynamical system, _ r = 0 0 r V 0 ; (7.3) where r = r 1 r 2 T and the superscript [] T denotes the transpose operation. We approximate the time derivative in (7.3) using the implicit nite dierence _ r = (r n+1 r n )=t, wheren is a discrete time index. To this end, we get r n+1 = A(u) r n + b(u); (7.4) where u = V T and A = 1 1 + ( t) 2 1 t t 1 ; b = V t 0 : (7.5) Visual sensing model We view the (r; ) plane as the perception space because r and dene the position of the target as perceived by the agent. We prescribe a conic-like boundary in the perception space which denes the `eld of view' of the agent (see inset of Fig. 7.1). Outside this boundary, the agent has no knowledge of the target's position and the agent has failed to nd it. We also prescribe another circular boundary around the origin which denes the `success region'. When the target enters the success region, the agent has found the target. We discretize the perception space such that the boundaries of eld of view and success region are dened as polygons in the discrete perception space, as depicted in Fig. 7.2. Basically, we lay down a set of pointsf i g and compute the Voronoi cells associated with these points. By denition, CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 76 -5 0 5 10 15 -10 -5 0 5 10 Figure 7.2: The perception space of the agent. The eld of view and success boundaries are depicted as red and green, respectively. The perception space is discretized into statesS =fs i g corresponding to the cells depicted with black lines. A given states i (highlighted in black) is acted upon by action a k yieldings k i (highlighted in blue). The compute the intersection of the mapped state with all reference states to determine a discrete probability transition matrix. thei th Voronoi cell of point i includes all points of the perception space that are closer to it than any other point j . LetS cells be the set of all voronoi cells;S cells is a discrete represention of the perception space. Discrete state and action spaces As the agent moves about, the target enters various cells in the discrete perception space. We dene the state of the agent at a given time as the cell in the perception space where the target is located. It will be necessary to dene two more states, s success ands fail , as the states into which the agent transitions when it enters the success region and when it exits the eld of view of the agent, respectively. We group these states intoS terminal =fs success ;s fail g because these are the states where we will terminate an episode of the agent's experience. We can now dene the state spaceS by computing the union of the voronoi cells of the perception space and the terminal statesS =S cells [S terminal . We dene the action spaceA as the set of control inputs u that the agent can choose from. For simplicity, we discretize the action space byA =fa k g, where a k = (V k ; k ) and V k and k correspond to a particular choice of translational speed and rotation rate. The particular choice of actions can be arbitrary as long asA is nite. Discrete dynamics The deterministic dynamics of (7.5) map onto the state space in the form of a stochastic process that arises entirely from the discretization of the perception space. At each time step n, we denote the cell that the agent is in as s n and the action it chooses as a n . To this CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 77 end, we dene the stochastic transition matrix via its componentsP k ij , P k ij =P[s n+1 =s j js n =s i ; a n =a k ]: (7.6) Here,P k ij or P[s n+1 = s j js n = s i ; a n = a k ] is the conditional probability of transitioning from states i tos j given a choice of actiona k . To calculateP k ij , we let the ane transformation in (7.5) act on cells i given actiona k to get the image polygons k i and compute the intersection of the image of s k i with cells s j to yield P k ij = area [s k i \s j ] areas i : (7.7) Meanwhile, the probability of transitioning from any state inS terminal to any other state is zero. Reward We compute the discrete outcome function o n by o n = 8 > < > : 1 s n =s success ; 1 s n =s fail ; 0 otherwise; (7.8) which is a variable used to determine whether the target is within the bounds of the agent's eld of view or if the agent has reached the target. We provide the agent with feedback on its performance in the form of a reward signal R n R n =o n R o ; (7.9) where R o > 0 is a positive constant. Namely, the agent receives a positive reward when it reaches the target, a negative reward it the target leaves its eld of view, and no intermediate rewards. We dene the reward functionR k i , R k i =E[R n+1 js n =s i ; a n =a k ]; (7.10) which is the expected reward given that the agent is in state s i and takes action a k . This can be computed given R k i = X j P k ij o j R o ; (7.11) whereR o o j is the reward given when transitioning into state s j given by (7.8) and (7.9). Together, P k ij andR k i form a model of how the agent transitions and what rewards it should expect given a particular choice of action. 7.2 Data-driven controller design In order to control the agent's behavior, we need to dene a policy (ajs) given by (a k js i ) =P[a n =a k js n =s i ]; (7.12) which is the probability of choosing action a k given that the agent is in states i . The policy can be thought of as a controller that maps states to control actions. The control problem is to design a CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 78 policy that optimizes the behavior of the agent. In order to optimize, we need to dene a measure of `goodness' for a policy as follows. First, we dene the returnG n asG n = P 1 m=0 m R n+m+1 , where2 [0; 1] is the discount factor. This is the sum of future discounted rewards. We discount to encode a sense of time where 1 is a `far-sighted' evaluation whereby we look far into the future to compute the return, whereas 0 is a `myopic' evaluation where we only consider short-term rewards. Note that, by denition, G n =R n+1 +G n+1 : (7.13) This will be a convenient relation later. Next, we dene the state value functionV (s) as V (s i ) =E [G n js n =s i ]; (7.14) whereE indicates an expectation taken assuming the agent follows policy . Similarly, we dene the state-action value functionQ (s;a) as Q (s i ;a k ) =E [G n js n =s i ; a n =a k ]: (7.15) The value functions for a given policy are related by V (s i ) = X k (a k js i )Q (s i ;a k ): (7.16) We can use the state value function to dene a partial ordering of policies such that 0 i V (s i )V 0(s i );8s i 2S: (7.17) In this way, we can dene the optimal policy and its state value functionV as = argmax V (s) and V (s) = max V (s): (7.18) If we can solve for the optimal value function, then we have found an optimal control scheme for our agent. We have now described all the elements of a Markov Decision Process (MDP), consisting of the tuplehS;A;P;R;i. A policy is considered a solution to the problem of an MDP, and an optimal solution is the optimal policy. Since the state and action spaces of our MDP are nite, there exists a unique optimal policy [103]. First, let us dene the greedy action a greedy (s) = argmax a Q(s;a); (7.19) which is the action that maximizes the value function. The greedy policy is then (a k js) = ( 1 a k =a greedy (s); 0 otherwise: (7.20) The optimal policy is the greedy one with respect to the optimal value functionQ . The task for the agent is to learn the optimal value function then to act greedily with respect to it. We use two approaches to solve the MDP system described above: a model-based approach, that has knowledge ofhP;Ri over the whole space and computes the optimal value function using dynamic programming, and a model-free approach that learns the optimal value function directly from experience. CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 79 Algorithm 1 Iterative Policy Evaluation Input policy , stopping criterion m = 0 Initialize ^ V 0 arbitrarily repeat for all s i 2S do for all a k 2A do ^ Q(s i ;a k ) =R k i + P j P k ij ^ V m (s j ) end for ^ V m+1 (s i ) = P k (a k js i )Q(s i ;a k ) end for =k ^ V m+1 ^ V m k m =m + 1 until < Model-based approach The solution entails computing the value functionV for a given policy, and optimizing the policy itself. It may seem dicult to do both simultaneously, but for a greedy policy, we can take advantage of some simplications. First, we describe how to compute the value function for a given policyV . A fundamental property of the value function is that it satises the Bellman equation [103] V (s i ) = X k (a k js i ) 2 4 R k i + X j P k ij V (s j ) 3 5 : (7.21) Basically, if we have a modelhP;Ri and a policy, we can compute the value functionV by solving the set of linear equations in (7.21). In practice, we solve these equations iteratively following the update rule ^ V m+1 (s i ) = X k (a k js i ) 2 4 R k i + X j P k ij ^ V m (s j ) 3 5 ; (7.22) where ^ V m is an incremental estimate of the value function. Notice thatV is a xed point of (7.22), such that lim m!1 ^ V m =V : Therefore, we can iteratively compute the value function using Algo- rithm 1. The optimal policy corresponds to greedily choosing actions with respect to Q , therefore from (7.21), we get the Bellman optimality equation V (s i ) = max k 0 @ R k i + X j P k ij V (s j ) 1 A : (7.23) This nonlinear set of equations cannot be solved directly. However, based on this equation, we can use a dynamic programming algorithm for `value iteration' ^ V m+1 (s i ) = max k 0 @ R k i + X j P k ij ^ V m (s j ) 1 A : (7.24) CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 80 Algorithm 2 Value Iteration Input stopping criterion m = 0 Initialize ^ V 0 arbitrarily repeat for all s i 2S do for all a k 2A do ^ Q(s i ;a k ) = (R k i + P j P k ij ^ V m (s j )) end for ^ V m+1 (s i ) = max k Q(s i ;a k ) end for =k ^ V m+1 ^ V m k m =m + 1 until < This algorithm converges exponentially fast in the parameter to the optimal value functionV , namely, lim m!1 ^ V m =V . This algorithm is described in Algorithm 2. Notice that by construction, we obtainQ (s;a), with respect to which we act greedily. Model-free approach Suppose that we do not have available or that it is impractical to compute hP;Ri and we want to learn the optimal policy in an online way. We are again faced with the issue of both learning the value function for a given policy as well as optimizing that policy itself. Here, our task is dicult because we do not explicitly know the model and therefore cannot estimate the value function by computing it over all states and actions. Instead, we have to experience various state and action pairs and build our map of the value function incrementally. For this reason, we learn the optimal state-action value functionQ (s;a) online, as it allows us to learn the optimal policy | recalling that (ajs) = greedyQ (s;a) | without ever using a model. One approach to learnQ is called temporal dierence (TD) learning, which we will explain as follows. We rst consider an algorithm to learn the value function for a given policy (ajs). Suppose we are in a states n and take an actiona n according to(a n js n ), resulting in rewardR n+1 and landing us in state s n+1 , at which point we select action a n+1 according to the policy. We start with an estimate of the value function ^ Q(s;a) from which we now compute the TD error, dened as n =R n+1 + ^ Q(s n+1 ;a n+1 ) ^ Q(s n ;a n ): (7.25) We can now introduce an update rule ^ Q n+1 (s n ;a n ) = ^ Q n (s n ;a n ) + n ; (7.26) where is a step size parameter, as described in Algorithm 3. If all state-action pairs are visited suciently often and satises the Robbins-Monroe stochastic estimation conditions [103], then the value function converges to the optimal value function, lim m!1 ^ Q m (s;a) =Q (s;a): If is kept constant, then the policy converges on average [103]. Now that we have an online way of learning the value function, we consider how to also optimize the policy. Since at every time step of TD we have available to us ^ Q, we can read our policy to act greedily with respect to it; we call this `exploitation' because we exploit our estimate of the value CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 81 Algorithm 3 TD Policy Evaluation Input policy , step size m = 0 Initialize ^ Q 0 arbitrarily repeat for each episode n = 0 Initialize s 0 repeat Select a n given by for s n Take a n , observe R n+1 and s n+1 Compute n ^ Q m+1 (s n ;a n ) = ^ Q m (s n ;a n ) + n n =n + 1 m =m + 1 until s n is terminal until satisfactory performance function to act more closely to optimal. But, because the value function estimate is non-optimal, we also want to `explore' by sometimes choosing the non-greedy action. Therefore, we dene an -greedy policy with respect toQ(s;a) as (a k js i ) = ( 1 1 1 jAj ; a k = argmax a l Q(s i ;a l ); jAj ; otherwise; (7.27) where < 1 is a small parameter determining the probability of acting non-greedily. Now, we simply apply TD learning and adjust the policy in an -greedy way online and we have the SARSA algorithm, described in Algorithm 4. The SARSA algorithm is called an on-policy algorithm because we compute n usinga n+1 where a is chosen from our current-greedy policy. What if we do not want to update based on this policy, but what we know to be the optimal policy | that is, acting truly greedily with respect to the value function? In this case, we can modify our denition of n as follows n =R n+1 + max k ^ Q(s n ;a k ) ^ Q(s n ;a n ): (7.28) When we use this expression to calculate the TD error, we get the Q-Learning algorithm, described in Algorithm 5. 7.3 Performance of trained agent We discretized the perception space as depicted in Fig. 7.2 resulting in a total number of statesjSj = 155 including the terminal states. We then discretized the action space by setting the linear velocity of the agent to a constant,V =V o , and its angular velocity =f1:5;1;0:5; 0; 0:5; 1; 1:5g o , leading to a total number of actionsjAj = 7. We choseV o = 1 and o = 0:5, leading to a minimum turning radius of V o =1:5 o = 0:667. The time step was set to t = 0:1. CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 82 Algorithm 4 SARSA Policy Optimization Input exploration parameter , step size m = 0 Initialize ^ Q 0 arbitrarily repeat for each episode n = 0 Initialize s 0 Choose a 0 from -greedy ^ Q(s 0 ;a) repeat Take a n , observe R n+1 and s n+1 Choose a n+1 from -greedy ^ Q(s n+1 ;a) Compute n ^ Q m+1 (s n ;a n ) = ^ Q m (s n ;a n ) + n n =n + 1 m =m + 1 until s n is terminal until satisfactory performance Algorithm 5 Q-Learning Policy Optimization Input exploration parameter , step size m = 0 Initialize ^ Q 0 arbitrarily repeat for each episode n = 0 Initialize s 0 repeat Choose a n from -greedy ^ Q(s n ;a) Take a n , observe R n+1 and s n+1 Compute n ^ Q m+1 (s n ;a n ) = ^ Q m (s n ;a n ) + n n =n + 1 m =m + 1 until s n is terminal until satisfactory performance Value Iteration We computed the optimal policy using the model-based algorithm (Algorithm 2). Sample trajectories are shown in Fig. 7.3(a) when setting the stopping criterion = 10 3 . The agent follows Dubins-like paths independent of initial conditions. The corresponding optimal value functionV is shown in Fig. 7.3(b). In the eld of view, the value function is high because, by following the optimal policy, the agent should nearly always expect to reach the target. In the area to the left of the target region, the value function is lower, because there are some initial positions within these states from which the agent will be unable to reach the target due to the minimum turning radius. We can evaluate the performance of this value function by assuming a greedy policy and sampling CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 83 -10 -5 0 5 10 -10 -5 0 5 10 -10 -5 0 5 10 -10 -5 0 5 10 -10 -5 0 5 10 -10 -5 0 5 10 -5 0 5 10 15 -10 -5 0 5 10 -5 0 5 10 15 -10 -5 0 5 10 -5 0 5 10 15 -10 -5 0 5 10 -100 -50 0 50 100 (a) (c) (e) (b) (d) (f) Figure 7.3: Top row: Trajectories sampled from policies at various stages of training (N train = 100; 1000; 4000). Bottom row: Value function after full training. (a) Model-based policy optimization using Value Iteration. Agent takes path directly to target. (b) Model-free policy optimization using SARSA. In early stages of training, agent takes circuitous path and sometimes does not reach target. As training progresses, the target is reached every time. (c) Model-free policy optimization using Q-Learning. Policy is learned early on and does not deviate much from optimal policy. episodes from it. One performance metric is the probability of failing to reach the target p fail . Because there exist some locations in the state space from which the agent will not be able to reach the target, even the optimal policy will have a non-zero p fail . We compute an average of p fail over 10,000 samples with initial conditions distributed uniformly over the state space, with the asymptotic value of p fail = 0:0494;. This value is depicted in Fig. 7.4(a) as a dashed line. SARSA To learn from experience, we implement the SARSA algorithm. We rst discretize the state space as before, and then we step forward the dynamics in (7.1). At each time step, we check which state the agent is in by computing the target's location (r; ) in the perception space. We then choose an action according to our current estimate of the value functionQ, acting -greedy ( = 0:01). We update the value function according to the SARSA update rule, and continue stepping the dynamics until the episode ends. At the end of an episode, we reinitialize the agent at a random location in the state space, and start the process again. This process is repeated over 20,000 episodes. Sample trajectories are shown in Fig. 7.3(c) after various rounds of training. Unlike in the model-based learning, the agent's path are more sinuous. The learned optimal value functionV is CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 84 0 25000 50000 0.01 0.1 1 0 25000 50000 0.1 0.5 1 0 25000 50000 0.5 0.7 1 Figure 7.4: Variation in performance parameters as a function of training episodes N train . (a) Probability of failure p fail for SARSA (blue) and Q-Learning (red), with smoothing window size 1000. Optimal value of p fail = 0:0494 shown as a dotted line. (b) Relative value function errorE value . (c) Fraction of suboptimal state-action pairs F suboptimal . shown in Fig. 7.3(d). During training, we can compute several measures of performance. First, we can compute a moving average of p fail , which we show in Fig. 7.4(a) with a smoothing window of 1000. The probability of failure goes down signicantly as training progresses. Compared to model-based optimal p fail , we see that on average the SARSA agent achieves the optimal failure rate. Another performance metric we use is to compare the computed value functionV to the optimal oneV . We do this by deningE value E value = 1 jSj X i V(s i )V (s i ) V (s i ) ; (7.29) which is the average normalized error over the states. In Fig. 7.4(b), we plotE value as a function of N train . As N train increases,E value converges to the average value of 0:25. That is, the SARSA agent has learned about 75% of the optimal value function. This is sucient to allow the agent to reach the target at least as often as in the optimal case. Finally, we compute how often the agent acts optimally. Using (7.19), we evaluate the learned greedy action ^ a greedy (s) for ^ Q and the optimal greedy action a greedy (s) forQ . We then compute F suboptimal = 1 jSj X i H ^ a greedy (s i )a greedy (s i ) ; (7.30) whereH() is the Heaviside function. That is,F suboptimal represents the fraction of states in which the greedy agent chooses a suboptimal action. In Fig. 7.4(c), we plotF suboptimal versus N train . Surprisingly, even thoughF suboptimal converges to 0.75, that is, the agent chooses a suboptimal action in 75% of the states, it still achieves a low failure rate of p fail 0:05, as shown in Fig. 7.4(a). Q-Learning We repeat the same learning process using the Q-Learning algorithm to update the value function. Sample trajectories and the learned value function are shown in Fig. 7.3(e,f). The performance measures are shown in Fig. 7.4 in solid red lines. While the details of the individual CHAPTER 7. DATA-DRIVEN SOURCE TRACKING USING VISUAL SENSING 85 trajectories are distinct, the statistical performance of the Q-learning are comparable to the SARSA algorithm, with slight improvement in convergence rate. 7.4 Discussion We demonstrated the ecacy of using reinforcement learning as a paradigm for learning the optimal control strategy in the goal-based problem of target nding. We used a relatively simple example, which has a known solution in terms of Dubins path, as a means of introducing and applying reinforcement learning techniques to solve adaptive optimal control problems. In particular, we found that an agent is capable of learning the optimal policy from experience. These results can be viewed as a rst step towards understanding the applicability of reinforcement learning to pursuit problems. Chapter 8 Conclusions and Outlook There is great need to develop autonomous submersibles. For example, the Navy is building un- manned underwater vehicles that can carry out reconnaissance tasks and deliver payloads to targets while reducing the risk to personnel. NASA is developing an autonomous submersible for explo- ration of the methane oceans on Saturn's moon, Titan. Along with the amount of nancial and human capital that been directed towards the development of autonomous aircraft and autonomous motor vehicles, both engineering corporations and government agencies are realizing there is an un- met need in the submersible market. These autonomous submersibles will likely leverage many of the technologies that have been developed for these other applications, however there are a unique set of technological challenges associated with operating underwater. The need for new basic science research to address these challenges is readily apparent. These submersibles will be equipped with a suite of sensors leveraging a variety of sensing modalities. One such modality is hydrodynamic sensing. Because uid ows exhibit a high degree of coherence and structure, and these ow structures persist in the uid for remarkably long times, it is likely that relevant cues from the environment are embedded in the uid ow itself. When we look to nature, we nd evidence that many organisms are capable of extracting information from the hydrodynamic environment that surrounds them. An autonomous submersible would be at a great advantage if it were capable of decoding these cues and making decisions based on them. With the aim of building a framework suitable for implementation in engineering systems, we explored both physics-based and data-driven modeling approaches for some cases where ow sensing can be used. Physics-based modeling consisted of taking concepts and equations that described our understanding of the relevant physics and using them to couple together the dierent parts of the ow sensing problem. Data-driven modeling, involved the acquisition of an extensive dataset of measurements and behaviors that we wanted to model. We then used a number of parameterized models and their associated tools to construct a model that ts the data optimally. When had knowledge of the relevant physics, we used it to inform our modeling of the coupling between sensory data and ow properties. In these instances, we we able to validate our intuition by comparing our results with examples from nature and from experiment. We also validated our ideas by performing numerical experiments and analytical investigations. In all cases, we demonstrated that knowledge of the relevant physical phenomena can provide a great deal of insight and inform powerful and robust models of ow sensing. In other instances, the governing equations did not provide a tractable means of building the 86 CHAPTER 8. CONCLUSIONS AND OUTLOOK 87 model, or they were simply not available to us. In these cases, data-driven modeling approaches proved useful. The use of data-driven modeling was necessitated by the fact that local ow measure- ments depend in a nonobvious way on global ow characteristics, meaning that physical intuition and heuristics failed to posit these relationships a priori. A most striking example of this phe- nomenon was the work where we classied wakes from local measurements of vorticity time series. The relationship between local information and global wake type was not apparently a priori, but a well-trained neural network was able to uncover the salient features for performing this task and did so with a remarkable degree of accuracy. However, we know that these models are only as good as the underlying data on which they are trained. In the case of wake classication, we provided a large, rich, and clean dataset of vorticity time series from numerical simulations. How these results extend to empirical data is not clear yet. However, we have great condence in our approach based on our investigations of the robustness of the neural network and the many instances of their power and exibility in other applications in the engineering world. We must take great care, however, to be cautious about where and when we use these data- driven models. As is the case with many modern machine-learning tools, how the neural network was capable of uncovering the relevant patterns in data was barely interpretable by the human user. This is because the structure of the neural network is so high-dimensional and abstract. This is what allows the networks to capture strong nonlinearities in the data, but it is also what makes them dicult to interpret. We aimed to demonstrate by our work that these results serve as a framework for addressing more complex and pressing challenges such as detection and localization of unknown objects in the far eld and autonomous navigation. Going forward, it is evident that we need to leverage both our knowledge of ow physics and modern tools in data science to construct informative and parsimonious models for ow sensing. Physics-based methods are promising for developing low- order models which capture relevant relationships in some cases. These models then enjoy the benet of our extensive analytical tools because we have the models in equation form. However, there are a number of instances where data is abundant and of reasonable quality, which may mean that using data-driven methods may be an appropriate modeling strategy. As the eld of data science progresses across all engineering applications, the use of data-driven modeling will become more prevalent in ow sensing. It will serve us well to keep abreast of the advances in data science while being sure to ground our studies rmly in the framework where the physics is well known, taking care to leverage the advantages of both modeling approaches to develop robust and powerful tools for ow sensing. Bibliography [1] J. B. Allen and L. R. Rabiner. A unied approach to short-time fourier analysis and synthesis. Proceedings of the IEEE, 65(11):1558{1564, 1977. [2] G. Arnold. Rheotropism in shes. Biological Reviews, 49(4):515{576, 1974. [3] R. Bachmayer and N. E. Leonard. Vehicle networks for gradient descent in a sampled envi- ronment. In Decision and Control, 2002, Proceedings of the 41st IEEE Conference on, pages 112{117, 2002. [4] A. D. Becker, H. Masoud, J. W. Newbolt, M. Shelley, and L. Ristroph. Hydrodynamic schooling of apping swimmers. Nature Communications, 6, 2015. [5] L. M. Belue and K. W. Bauer. Determining input features for multilayer perceptrons. Neu- rocomputing, 7(2):111{121, 1995. [6] G. Berkooz, P. Holmes, and J. L. Lumley. The proper orthogonal decomposition in the analysis of turbulent ows. Annual Review of Fluid Mechanics, 25(1):539{575, 1993. [7] E. Biyik and M. Arcak. Gradient climbing in formation via extremum seeking and passivity- based coordination rules. Asian Journal of Control, 10(2):201{211, 2008. [8] M. Bozkurttas, R. Mittal, H. Dong, G. Lauder, and P. Madden. Low-dimensional models and performance scaling of a highly deformable sh pectoral n. Journal of Fluid Mechanics, 631:311{342, 2009. [9] J. B. Bratt. Flow patterns in the wake of an oscillating aerofoil. Technical report, Ministry of Supply { Aeronautical Research Council, 1953. [10] J. H. Buchholz and A. J. Smits. The wake structure and thrust performance of a rigid low-aspect-ratio pitching panel. Journal of Fluid Mechanics, 603:331{365, 2008. [11] E. J. Candes. The restricted isometry property and its implications for compressed sensing. Comptes Rendus Mathematique, 346(9):589{592, 2008. [12] K. K. Chen, J. H. Tu, and C. W. Rowley. Variants of dynamic mode decomposition: boundary condition, Koopman, and Fourier analyses. J. Nonlinear Sci., 22:887{915, 2012. [13] J. Cochran, E. Kanso, S. D. Kelly, H. Xiong, and M. Krstic. Source seeking for two nonholo- nomic models of sh locomotion. IEEE Transactions on Robotics, 25(5):1166{1176, 2009. 88 BIBLIOGRAPHY 89 [14] J. Cochran and M. Krstic. Nonholonomic source seeking with tuning of angular velocity. IEEE Transactions on Automatic Control, 54(4):717{731, 2009. [15] B. Colvert, M. Alsalman, and E. Kanso. Classifying vortex wakes using neural networks. Bioinspiration & Biomimetics, 2018. [16] B. Colvert, K. K. Chen, and E. Kanso. Bioinspired sensory systems for shear ow detection. J. of Nonlinear Sci., Accepted Jan. 2017., 2017. [17] B. Colvert, K. K. Chen, and E. Kanso. Local ow characterization using bioinspired sensory information. Journal of Fluid Mechanics, 818:366{381, 2017. [18] B. Colvert and E. Kanso. Fishlike rheotaxis. J. Fluid Mech., 793:656{666, 2016. [19] B. Colvert, G. Liu, H. Dong, and E. Kanso. How can a source be located by responding to local information in its hydrodynamic trail? In Decision and Control (CDC), 2017 IEEE 56th Annual Conference on, pages 2756{2761, 2017. [20] S. Coombs, J. Janssen, and J. F. Webb. Diversity of lateral line systems: evolutionary and functional considerations. In Sensory biology of aquatic animals, pages 553{593. Springer, 1988. [21] S. Coombs and S. Van Netten. The hydrodynamics and structural mechanics of the lateral line system. In Fish Biomechanics, volume 23 of Fish Physiology, pages 103{139. Academic Press, 2005. [22] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and N. R. Franks. Collective memory and spatial sorting in animal groups. Journal of theoretical biology, 218(1):1{11, 2002. [23] H. Dahlkamp, A. Kaehler, D. Stavens, S. Thrun, and G. Bardski. Self-supervised monocular road detection in desert terrain. In Robotics: Science and Systems, volume 38. Philadelphia, 2006. [24] A. Davies. Boeing's monstrous underwater robot can wander the ocean for 6 months. Wired, March 2016. [25] G. Dehnhardt, B. Mauck, W. Hanke, and H. Bleckmann. Hydrodynamic trail-following in harbor seals (Phoca vitulina). Science, 293:102{104, 2001. [26] L. DeVries, F. D. Lagor, H. Lei, X. Tan, and D. A. Paley. Distributed ow estimation and closed-loop control of an underwater vehicle with a multi-modal articial lateral line. Bioinspiration & Biomimetics, 10(2):025002, 2015. [27] L. D. DeVries and D. A. Paley. Wake estimation and optimal control for autonomous aircraft in formation ight. In Proceedings of the AIAA Guidance, Navigation, and Control Conference, page 4705, 2013. [28] E. G. Drucker and G. V. Lauder. Locomotor forces on a swimming sh: three-dimensional vortex wake dynamics quantied using digital particle image velocimetry. Journal of Experi- mental Biology, 202(18):2393{2412, 1999. BIBLIOGRAPHY 90 [29] J. D. Eldredge. Dynamically coupled uid?body interactions in vorticity-based numerical simulations. J. Comput. Phys., 227:9170{9194, 2008. [30] J. Engelmann, W. Hanke, J. Mogdans, and H. Bleckmann. Neurobiology: Hydrodynamic stimuli and the sh lateral line. Nature, 408(6808):51{52, 11 2000. [31] A. Esteva, B. Kuprel, R. A. Novoa, J. Ko, S. M. Swetter, H. M. Blau, and S. Thrun. Dermatologist-level classication of skin cancer with deep neural networks. Nature, 542(7639):115{118, 2017. [32] J. A. Farrell, S. Pang, W. Li, and R. Arrieta. Chemical plume tracing experimental results with a remus auv. In OCEANS, volume 2, pages 962{968, 2003. [33] V. I. Fernandez. Performance analysis for lateral-line-inspired sensor arrays. PhD thesis, Massachusetts Institute of Technology, 2011. [34] A. Filella, F. Nadal, C. Sire, E. Kanso, and C. Eloy. Hydrodynamic interactions in uence sh collective behavior. arXiv:1705.07821, 2107. [35] B. A. Free and D. Paley. Model-based observer and feedback control design for a rigid joukowski foil in a k arm an vortex. Bioinspiration & Biomimetics, 2018. [36] A. Gao and M. Triantafyllou. Bio-inspired Pressure Sensing for Active Yaw Control of Un- derwater Vehicles. PhD thesis, Massachusetts Institute of Technology, 2013. [37] M. Gazzola, P. Chatelain, W. van Rees, and P. Koumoutsakos. Simulations of single and mul- tiple swimmers with non-divergence free deforming geometries. J. Comput. Phys., 230:7093{ 7114, 2011. [38] M. Gazzola, A. A. Tchieu, D. Alexeev, A. de Brauer, and P. Koumoutsakos. Learning to school in the presence of hydrodynamic interactions. Journal of Fluid Mechanics, 789:726{749, 2016. [39] R. Godoy-Diana, J.-L. Aider, and J. E. Wesfreid. Transitions in the wake of a apping foil. Physical Review E, 77(1):16308, 2008. [40] G. H. Golub and C. F. Van Loan. Matrix computations. Johns Hopkins U. Press, 3 edition, 1996. [41] M. P. Groover. Automation. Encyclopdia Britannica, March 2018. [42] E. Guttierrez. Navy tests new unmanned underwater vehicle at jeblc-fs. Navy News, December 2014. [43] G. Haller. An objective denition of a vortex. Journal of Fluid Mechanics, 525:1{26, 2005. [44] G. Haller. Lagrangian coherent structures. Annual Review of Fluid Mechanics, 47:137{162, 2015. [45] G. Haller and G. Yuan. Lagrangian coherent structures and mixing in two-dimensional tur- bulence. Physica D: Nonlinear Phenomena, 147(3-4):352{270, 2000. [46] D. Hambling. How the soviet union snooped waters for enemy subs|without sonar. Popular Mechanics, 2017. BIBLIOGRAPHY 91 [47] A. Hamdi. The recovery of a time-dependent point source in a linear transport equation: application to surface water pollution. Inverse Problems, 25(7):075006, 2009. [48] A. L. Hammond. Project famous: Exploring the mid-atlantic ridge. Science, 187(4179):823{ 825, March 1975. [49] E. S. Hassan. Mathematical description of the stimuli to the lateral line system of sh derived from a three-dimensional ow eld analysis. Biological Cybernetics, 66(5):443{452, 1992. [50] R. Hecht-Nielsen. Theory of the backpropagation neural network. In International Joint Conference on Neural Networks, pages 593{605. IEEE, 1989. [51] M. S. Hemati, J. D. Eldredge, and J. L. Speyer. Wake sensing for aircraft formation ight. Journal of Guidance, Control, and Dynamics, 37(2):513{524, 2014. [52] K. Hornik, M. Stinchcombe, and H. White. Multilayer feedforward networks are universal approximators. Neural Networks, 2:35{366, 1989. [53] Y. Huang, J. Yen, and E. Kanso. Detection and tracking of chemical trails in bio-inspired sensory systems. European Journal of Computational Mechanics, 26(1-2):98{114, 2017. [54] J. C. R. Hunt, A. A. Wray, and P. Moin. Eddies, streams, and convergence zones in turbulent ows. In Proceeings of the Summer Program 1988, pages 193{208, Stanford, CA, USA, 1988. Center for Turbulence Research. Article no. N89-24555. [55] N. S. Jankel. Ai vs. human intelligence: Why computers will never create disruptive innova- tions. Hungton Post, April 2015. [56] J. Jeong and F. Hussain. On the identication of a vortex. Journal of Fluid Mechanics, 285:69{94, 1995. [57] J. Kasten, C. Petz, I. Hotz, H.-C. Hege, B. R. Noack, and G. Tadmor. Lagrangian feature extraction of the cylinder wake. Physics of Fluids, 22(9):091108, 2010. [58] A. Klein and H. Bleckmann. Determination of object position, vortex shedding frequency and ow velocity using articial lateral line canals. Beilstein Journal of Nanotechnology, 2:276{283, 2011. [59] M. Koochesfahani. Wake of an oscillating airfoil. Physics of Fluids, 29(9):2776{2776, 1986. [60] M. M. Koochesfahani. Vortical patterns in the wake of an oscillating airfoil. AIAA Journal, 27(9):1200{1205, 1989. [61] A. B. Kroese and N. A. Schellart. Velocity- and acceleration-sensitive units in the trunk lateral line of the trout. Journal of Neurophysiology, 68(6):2212{2221, 1992. [62] T. Kuo and S. Chan. Continuous, on-line, real-time spectral analysis of systemic arte- rial pressure signals. American Journal of Physiology-Heart and Circulatory Physiology, 264(6):H2208{H2213, 1993. [63] R. R. Lagnado. The stability of two-dimensional linear ows. PhD thesis, California Institute of Technology, 1985. BIBLIOGRAPHY 92 [64] J. Lai and M. Platzer. Jet characteristics of a plunging airfoil. AIAA Journal, 37(12):1529{ 1537, 1999. [65] H. Lamb. Hydrodynamics. Cambridge University Press, 1932. [66] M. Luhar, A. S. Sharma, and B. J. McKeon. Opposition control within the resolvent analysis framework. Journal of Fluid Mechanics, 749:597{626, 2014. [67] A. P. Maertens and M. S. Triantafyllou. The boundary layer instability of a gliding sh helps rather than prevents object identication. Journal of Fluid Mechanics, 757:179{207, 2014. [68] G. Marrucci and G. Astarita. Linear, steady, two-dimensional ows of viscoelastic liquids. AIChE Journal, 13:931{935, 1967. [69] B. McKeon and A. Sharma. A critical-layer framework for turbulent pipe ow. Journal of Fluid Mechanics, 658:336{382, 2010. [70] M. V. Melander and F. Hussain. Polarized vorticity dynamics on a vortex column. Physics of Fluids A: Fluid Dynamics, 5(8):1992{2003, 1993. [71] S. Michelin and S. G. L. Smith. An unsteady point vortex method for coupled uid{solid problems. Theo. and Comp. Fluid Dynamics, 23(2):127{153, 2009. [72] R. Mittal, H. Dong, M. Bozkurttas, F. Najjar, A. Vargas, and A. von Loebbecke. A versatile sharp interface immersed boundary method for incompressible ows with complex boundaries. Journal of Computational Physics, 227(10):4825{4852, 2008. [73] J. C. Montgomery, C. F. Baker, and A. G. Carton. The lateral line can mediate rheotaxis in sh. Nature, 389(6654):960{963, 1997. [74] J. C. Montgomery, S. Coombs, and C. F. Baker. The mechanosensory lateral line system of the hypogean form of astyanax fasciatus. Environmental Biol. of Fishes, 62(1-3):87{96, 2001. [75] B. J. Moore and C. Canudas-de Wit. Source seeking via collaborative measurements by a circular formation of agents. In American Control Conference 2010, Proceedings of the, pages 6417{6422. IEEE, 2010. [76] U. K. M uller, J. G. van den Boogaart, and J. L. van Leeuwen. Flow patterns of larval sh: undulatory swimming in the intermediate ow regime. Journal of Experimental Biology, 211(2):196{205, 2008. [77] B. R. Noack and H. Eckelmann. A global stability analysis of the steady and periodic cylinder wake. Journal of Fluid Mechanics, 270:297{330, 1994. [78] G. Novati, S. Verma, D. Alexeev, D. Rossinelli, W. M. van Rees, and P. Koumoutsakos. Syn- chronisation through learning for two self-propelled swimmers. Bioinspiration & Biomimetics, 12(3):036001, 2017. [79] D. of Science and Technology. Intelligence report: Soviet antisubmarine warfare: Current capabilities and priorities. Technical report, Central Intelligence Agency, 1972. BIBLIOGRAPHY 93 [80] P. Ogren, E. Fiorelli, and N. E. Leonard. Cooperative control of mobile sensor networks: Adaptive gradient climbing in a distributed environment. IEEE Transactions on Automatic Control, 49(8):1292{1302, 2004. [81] S. Oleson. Titan submarine: Exploring the depths of kraken. June 2014. [82] D. Outreach. Actuv \sea hunter" prototype transitions to oce of naval research for further development. Defense Advanced Research Projects Agency, 2018. [83] W. Pease. An automatic machine tool. Scientic American, 187(3):101{115, September 1952. [84] C. S. Peskin. Numerical analysis of blood ow in the heart. Journal of Computational Physics, 25(3):220{252, 1977. [85] M. F. Platzer, K. D. Jones, J. Young, and J. C. Lai. Flapping-wing aerodynamics: progress and challenges. AIAA Journal, 46(9):2136, 2008. [86] Z. Ren and K. Mohseni. A model of the lateral line of sh for vortex sensing. Bioinspiration and Biomim., 7(3):036016, 2012. [87] I. I. Report. Prospects for soviet success in improving detection of submarines. Technical report, 1974. [88] L. Ristroph, J. C. Liao, and J. Zhang. Lateral line layout correlates with the dierential hydrodynamic pressure on swimming sh. Phys. Rev. Lett., 114:018102, Jan 2015. [89] G. Roberts, R. Sutton, A. Zirilli, and A. Tiano. Intelligent ship autopilots{a historical per- spective. Mechatronics, 13:1091{1103, 2003. [90] C. W. Rowley, I. Mezi c, S. Bagheri, P. Schlatter, and D. S. Henningson. Spectral analysis of nonlinear ows. Journal of Fluid Mechanics, 641:115{127, 2009. [91] T. Salum ae, I. Ran o, O. Akanyeti, and M. Kruusmaa. Against the ow: A braitenberg controller for a sh robot. In Robotics and Automation (ICRA), 2012 IEEE International Conference on, pages 4210{4215. IEEE, 2012. [92] L. Schielicke, P. N evir, and U. Ulbrich. Kinematic vorticity number{a tool for estimating vortex sizes and circulations. Tellus A: Dynamic Meteorology and Oceanography, 68(1):29464, 2016. [93] P. J. Schmid. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics, 656:5{28, 2010. [94] J. Schmidhuber. Deep learning in neural networks: An overview. Neural Networks, 61:85{117, 2015. [95] T. Schnipper, A. Andersen, and T. Bohr. Vortex wakes of a apping foil. J. Fluid Mech., 633:411{423, 2009. [96] S. C. Shadden, F. Lekien, and J. E. Marsden. Denition and properties of lagrangian coherent structures from nite-time lyapunov exponents in two-dimensional aperiodic ows. Physica D: Nonlinear Phenomena, 212(3-4):271{304, 2005. BIBLIOGRAPHY 94 [97] S. C. Shadden, F. Lekien, J. D. Paduan, F. P. Chavez, and J. E. Marsden. The correla- tion between surface drifters and coherent structures based on high-frequency radar data in monterey bay. Deep Sea Research Part II: Topical Studies in Oceanography, 56(3-5):161{172, 2009. [98] T. B. Sheridan and W. L. Verplank. Human and computer control of undersea teleoperators. Technical report, Massachusetts Institute of Technology { Man-Machine Systems Laboratory, 1978. [99] Y. I. Siregar. Morphology and growth of lateral line organs of the rainbow trout (oncorhynchus mykiss). Acta Zoologica, 75(3):213{218, 1994. [100] G. R. Spedding. Wake signature detection. Annual Review of Fluid Mechanics, 46:273{302, 2014. [101] P. Stoica and R. L. Moses. Spectral Analysis of Signals, volume 452. Prentice Hall, 2005. [102] M. A. Stremler and S. Basu. On point vortex models of exotic blu body wakes. Fluid Dynamics Research, 46(6):061410, 2014. [103] R. S. Sutton and A. G. Barto. Reinforcement learning: An introduction. MIT Press Cam- bridge, 2 edition, 2017. [104] B. Tiko and H. Fossen. The limitations of three-dimensional kinematic vorticity analysis. Journal of Structural Geology, 17(12):1771{1784, 1995. [105] J. Titelman. Swimming and escape behavior of copepod nauplii: implications for predator- prey interactions among copepods. Marine Ecology-Progress Series, 213:203{213, 2001. [106] M. S. Triantafyllou, G. D. Weymouth, and J. Miao. Biomimetic survival hydrodynamics and ow sensing. Annual Review of Fluid Mechanics, 48:1{24, 2016. [107] C. Truesdell. Two measures of vorticity. Journal of Rational Mechanics and Analysis, 2:173{ 217, 1953. [108] A. C. H. Tsang and E. Kanso. Dipole interactions in doubly periodic domains. Journal of Nonlinear Science, 23(6):971{991, 2013. [109] J. H. Tu. Machine learning for classication of vortex patterns generated by pitching and plunging plates. In Bulletin of the American Physical Society. American Physical Society { Division of Fluid Dynamics, 2017. [110] E. D. Tytell and G. V. Lauder. The hydrodynamics of eel swimming. Journal of Experimental Biology, 207(11):1825{1841, 2004. [111] A. Vargas, R. Mittal, and H. Dong. A computational study of the aerodynamic performance of a dragon y wing section in gliding ight. Bioinspiration & Biomimetics, 3(2):026004, 2008. [112] R. Venturelli, O. Akanyeti, F. Visentin, J. Je zov, L. D. Chambers, G. Toming, J. Brown, M. Kruusmaa, M. W. Megill, and P. Fiorini. Hydrodynamic pressure sensing with an articial lateral line in steady and unsteady ows. Bioinspiration and Biomim., 7(3):036004, 2012. BIBLIOGRAPHY 95 [113] M. Vergassola, E. Villermaux, and B. I. Shraiman. `infotaxis' as a strategy for searching without gradients. Nature, 445(7126):406{409, 2007. [114] S. Verma, G. Novati, and P. Koumoutsakos. Ecient collective swimming by harnessing vortices through deep reinforcement learning. arXiv preprint arXiv:1802.02674, 2018. [115] V. Vinge. Signs of the singularity. IEEE Spectrum, June 2008. [116] C. Von Campenhausen, I. Riess, and R. Weissert. Detection of stationary objects by the blind cave shanoptichthys jordani (characidae). Journal of Comparative Physiology, 143(3):369{ 374, 1981. [117] A. M. G. Walan. Anti-submarine warfare (asw) continuous trail unmanned vessel (actuv). Defense Advanced Research Projects Agency. [118] M. Wang and M. S. Hemati. Detecting exotic wakes with hydrodynamic sensors. arXiv preprint arXiv:1711.10576, 2017. [119] J. F. Webb. Lateral line morphology and development and implications for the ontogeny of ow sensing in shes. In Flow Sensing in Air and Water. Springer, 2014. [120] P. Werbos. Backpropagation through time: What it does and how to do it. Proceedings of the IEEE, 78(19):1550{1560, 1990. [121] C. H. K. Williamson and A. Roshko. Vortex formation in the wake of an oscillating cylinder. J. Fluids Struct., 2:355{381, 1988. [122] M. Wolfgang, J. Anderson, M. Grosenbaugh, D. Yue, and M. Triantafyllou. Near-body ow dynamics in swimming sh. Journal of Experimental Biology, 2-2(17):2303{2327, 1999. [123] K. Yanase, N. A. Herbert, and J. C. Montgomery. Disrupted ow sensing impairs hydrody- namic performance and increases the metabolic cost of swimming in the yellowtail kingsh, seriola lalandi. The J. of Experimental Biol., 215(22):3944{3954, 2012. [124] Y. Yang, J. Chen, J. Engel, S. Pandya, N. Chen, C. Tucker, S. Coombs, D. L. Jones, and C. Liu. Distant touch hydrodynamic imaging with an articial lateral line. Proc. of the Nat. Acad. of Sci., 103:18891{18895, Dec. 2006. [125] J. Yen, P. H. Lenz, D. V. Gassie, and D. K. Hartline. Mechanoreception in marine copepods: electrophysiological studies on the rst antennae. Journal of Plankton Research, 14(4):495{ 512, 1992. [126] J. Yen, M. J. Weissburg, and M. H. Doall. The uid physics of signal perception by mate- tracking copepods. Philos. Trans. R. Soc. London B: Biol. Sci., 353:787{804, 1998. [127] W.-K. Yen, D. M. Sierra, and J. Guo. Controlling a robotic sh to swim along a wall using hydrodynamic pressure feedback. IEEE Journal of Oceanic Engineering, 2018. [128] M. Yoshizawa, W. Jeery, S. Van Netten, and M. McHenry. The sensitivity of lateral line receptors and their role in the behavior of mexican blind cavesh (astyanax mexicanus). The J. of Experimental Biol., pages 886{895, 2013. BIBLIOGRAPHY 96 [129] A. Ysasi, E. Kanso, and P. K. Newton. Wake structure of a deformable joukowski airfoil. Physica D: Nonlinear Phenomena, 240(20):1574{1582, 2011.

Abstract (if available)

Abstract From naval applications to space exploration, there is a growing need to develop engineering solutions for autonomous submersibles. If subsurface vehicles are going to be capable of high levels of autonomy, they will need the ability to sense and decode information embedded in the fluid flows that surround them and make decisions based on the cues they extract. Inspired by the capabilities of many aquatic creatures that leverage various flow sensing modalities to achieve critical tasks such as navigation, mating, and predation, we use both physics-based and data-driven techniques to model how biological and bioinspired sensory systems can extract and react to relevant cues in the hydrodynamic environment. With a physics-based framework, we address the problem of optimizing sensory layouts, develop a local flow characterization scheme, and test hypotheses about behavioral responses to flow stimuli. We then use data-driven techniques to identify and decipher which cues in the hydrodynamic environment are most relevant for achieving classification of unknown objects in the flow field. These results serve as a framework for addressing more complex and pressing challenges such as detection and localization of unknown objects in the far field and autonomous navigation.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

PDF

Evaluating sensing and control in underwater animal behaviors

PDF

Synchronization in carpets of model cilia: numerical simulations and continuum theory

PDF

Mechanical and flow interactions facilitate cooperative transport and collective locomotion in animal groups

PDF

Transitions in the collective behavior of microswimmers confined to microfluidic chips

PDF

Aerodynamics at low Re: separation, reattachment, and control

PDF

Modeling and analysis of parallel and spatially-evolving wall-bounded shear flows

PDF

RANS simulations for flow-control study synthetic jet cavity on NACA0012 and NACA65(1)412 airfoils.

PDF

Multiscale modeling of cilia mechanics and function

PDF

Modeling and simulation of circulating tumor cells in flow. Part I: Low-dimensional deformation models for circulating tumor cells in flow; Part II: Procoagulant circulating tumor cells in flow

PDF

Boundary layer and separation control on wings at low Reynolds numbers

PDF

Model based design of porous and patterned surfaces for passive turbulence control

PDF

Physics-based data-driven inference

PDF

Feature and model based biomedical system characterization of cancer

PDF

Novel techniques for analysis and control of traffic flow in urban traffic networks

PDF

Speeding up trajectory planning for autonomous robots operating in complex environments

PDF

Tracking and evolution of compressible turbulent flow structures

PDF

Feature learning for imaging and prior model selection

PDF

Optimization-based whole-body control and reactive planning for a torque controlled humanoid robot

PDF

Exploiting mechanical properties of bipedal robots for proprioception and learning of walking

Asset Metadata

Creator Colvert, Brendan Thomas (author)

Core Title Physics-based and data-driven models for bio-inspired flow sensing and motion planning

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Aerospace Engineering

Publication Date 07/31/2018

Defense Date 05/07/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag bioinspired engineering,control theory,fish behavior,fluid dynamics,machine learning,OAI-PMH Harvest

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Kanso, Eva (committee chair), Luhar, Mitul (committee member), Newton, Paul (committee member), Savla, Ketan (committee member), Spedding, Geoffrey (committee member)

Creator Email brendancolvert@gmail.com,colvert@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-44559

Unique identifier UC11672614

Identifier etd-ColvertBre-6579.pdf (filename),usctheses-c89-44559 (legacy record id)

Legacy Identifier etd-ColvertBre-6579.pdf

Dmrecord 44559

Document Type Dissertation

Format application/pdf (imt)

Rights Colvert, Brendan Thomas

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA