ABSTRACT Title of dissertation: Optimality, Synthesis and a Continuum Model for Colle tive Motion Udit Halder, Do tor of Philosophy, 2019 Dissertation dire ted by: Professor P. S. Krishnaprasad Dept. of Ele tri al & Computer Engineering It is of importan e to study biologi al olle tives and apply the wisdom so a rued to modern day engineering problems. In this dissertation we attempt to gain insight into olle tive behavior where the main ontribution is twofold. First, a `bottom-up' approa h is employed to study individual level ontrol law synthesis and emergen e thereby of olle tive behavior. Three dierent problems, involv- ing single and multiple agents, are studied by both analyti al and experimental means. These problems arise from either a pra ti al viewpoint or from attempts at des ribing biologi ally plausible feedba k me hanisms. One result obtained in this ontext for a double agent s enario is that under a parti ular onstant bearing pursuit strategy, the problem exhibits ertain features ommon with the Kepler two body problem. Laboratory demonstrations of the solutions to these problems are presented. It is to be noted that these types of individual level ontrol prob- lems an help understand and onstru t building blo ks for group level behaviors. The se ond approa h is `top-down' in nature. It treats a olle tive as a whole and asks if its movement minimizes some kind of energy fun tional. A key goal of this work is to develop wave equations and their solutions for a natural lass of optimal ontrol problems with whi h one an analyze information transfer in o ks. Controllability arguments in innite dimensional spa es give strong sup- port to onstru t solutions for su h optimal ontrol problems. Sin e the optimal ontrol problems are innite dimensional in the state spa e and one annot simply expe t Pontryagin's Maximum Prin iple (PMP) to apply in su h a setting, the work has required are and attention to fun tional analyti onsiderations. In this work, it is shown that under a ertain assumption on nite o-dimensionality of a rea hable set, PMP remains valid. This assumption is then shown to hold true for the ase of a spe i ensemble of agents, ea h with state spa e as the Heisenberg group H(3). Moreover, analysis of optimal ontrols demonstrates the existen e of traveling wave solutions in that setting. Syn hronization results are obtained in a high oupling limit where deviation from neighbors is too ostly for every agent. The ombination of approa hes based on PMP and al ulus of variations have been fruitful in developing a solid new understanding of wave phenomena in olle tives. We provide partial results along these lines for the ase of a ontinuum of planar agents (SE(2) ase). Finally, a dierent top-down and data-driven approa h to analyze olle tive be- havior is also put forward in this thesis. It is known that the total kineti energy of a o k an be divided into several modes attributed to rigid-body translations, rotations, volume hanges, et . Flight re ordings of multiple events of European starling o ks yield time-signals of these dierent energy modes. This approa h then seeks an explanation of kineti energy mode distributions (viewed as o k- s ale de isions) by appealing to te hniques from evolutionary game theory and optimal ontrol theory. We propose the notion of ognitive ost that al ulates a suitably dened a tion fun tional and measures the ost to an event, resulting from temporal variations of energy mode distributions. Optimality, Synthesis and a Continuum Model for Colle tive Motion by Udit Halder Dissertation submitted to the Fa ulty of the Graduate S hool of the University of Maryland, College Park in partial fulllment of the requirements for the degree of Do tor of Philosophy 2019 Advisory Committee: Professor P. S. Krishnaprasad, Chair/Advisor Professor Steve Mar us Professor Nuno Martins Professor André Tits Professor Nikhil Chopra, Dean's Representative © Copyright by Udit Halder 2019 Dedi ated to my parents ii A knowledgments I would like to take this opportunity to thank my advisor Prof. P. S. Krish- naprasad for all the help and support he has generously provided over the past ve years. He has given me the freedom to think deeply about hallenging and fundamental problems. It is his onstant en ouragement and enthusiasm for s i- en e over ountless hours of enlightening dis ussions that positively propelled the work presented in this thesis. His persistent emphasis on attention to detail in every aspe t of resear h has helped me grow as a resear her. It has been an ab- solute pleasure having him as my thesis advisor. I would also like to thank Prof. Steve Mar us, Prof. Nuno Martins, Prof. André Tits and Prof. Nikhil Chopra for being in my dissertation ommittee and reviewing this do ument. The insightful omments from the members of the dissertation ommittee have helped improve the quality of this thesis. Collaboration with Dr. Uroš Kalabi¢ at the Mitsubishi Ele tri Resear h Laboratory has been very helpful. I also ollaborated with Dr. Eri Justh on a entral problem of my thesis and I thank him for all his feedba k and wise opinions. Over the past ve years, I have been fortunate to intera t with my olleagues in the Intelligent Servosystems Laboratory. I thank Ms. Vidya Raju, Dr. Biswadip Dey, Dr. Yunlong Huang, Prof. Kevin Galloway, Mr. Brent S hlotfeldt, and Mr. Kenneth Miltenberger for all the stimulating dis ussions we have had. I enjoyed a lot while ollaborating with Dr. Biswadip Dey, Mr. Brent S hlotfeldt, and Ms. Vidya Raju at several stages of my PhD. Spe ial thanks to Vidya who has been iii on a same boat as me from the start of my graduate study. I would also like to thank my friends outside the lab without whose ompany and support my life as a graduate student would have been in omplete. Biswadip, Vidya, Dipankar, Sid- dharth, Soham, Debdipta, Abhishek, Rajdeep, Ankit, Proloy, Soumyadip, Agniv, Jiaul, Kushal and others  thank you all for making this journey a memorable one. Finally, I would like to thank my parents who have been a sour e of onstant inspiration and en ouragement. I may not have a omplished this without their support. This resear h was supported by the Air For e O e of S ienti Resear h under AFOSR Grant FA9550-10-1-0250, an AFOSR FY2012 DURIP Grant No. FA2386-12-1-3002, the ARL/ARO MURI Program Grant No. W911NF-13-1- 0390, through the University of California Davis (as prime), the ARL/ARO Grant No. W911NF-17-1-0156, through the Virginia Polyte hni Institute and State University (as prime), and by Northrop Grumman Corporation. iv Table of Contents List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .viii List of Abbreviations and Notations.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Mathemati al Ba kground . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.1 Self Steering Parti le Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Experimental Setup. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 I Synthesis of Colle tive Motion: Bottom-up Approa hes 9 2 Feedba k Laws for Colle tive Motion .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Steering for Bea on Pursuit under Limited Sensing . . . . . . . . . . . . . . . . . . . 11 2.2.1 Tra king a Moving Bea on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2 Dynami s Restri ted to the Invariant Manifold . . . . . . . . . . . . . . . 15 π 2.2.3 Spe ial Case: α = . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2 2.2.4 The Limited Field of View (FOV) Problem . . . . . . . . . . . . . . . . . . . 23 2.2.5 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.2.6 Asso iated Lagrangian . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Biomimeti Algorithms for Coordinated Motion . . . . . . . . . . . . . . . . . . . . . . 33 2.3.1 Mutual Motion Camouage (MMC) . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3.2 Topologi al Velo ity Alignment (TVA) . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.3 Implementation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3 Optimal Steering of Agents on a Plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 Optimal Steering of a Uni y le . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.1 Optimal Control Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.2.2 Chara terizing the Types of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.2.3 On Time-optimality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.3 Optimal Control of a Colle tive of Agents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.4 Con luding Remarks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 v II Analysis of Colle tive Motion: Top-down Approa hes 65 4 Continuum Flo king and Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.1 Motivation.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 A Control System on a Loop Group.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3 Controllability.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4 Optimal Control Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.4.1 Cal ulus of Variations Approa h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.4.2 Maximum Prin iple Approa h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.5 Spe ial Case : G = H(3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.5.1 Controllability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.5.2 Equations of Optimal Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.5.3 Behavior of Optimal Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.5.4 Strong Coupling Limit, χ → ∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.5.5 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.6 A Continuum of Agents on the Plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.6.1 Equations of Optimal Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.6.2 Strong Coupling Limit, χ → ∞ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.6.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.7 Dis ussion and S ope of Future Resear h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5 Cognitive Cost of Flo king: A Geometri and Hamiltonian Per- spe tive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .120 5.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.2 Flo king Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.3 Data Smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.3.1 A linear generative model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.3.2 Data smoothing in the Eu lidean setting . . . . . . . . . . . . . . . . . . . . . . 125 5.4 Energy Modes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 5.5 Generative model on the 1-simplex and the data-smoothing problem131 5.6 Data Fitting Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.7 Dis ussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 6 Con lusions and Dire tions for Future Resear h . . . . . . . . . . . . . . . . . .148 A An Optimal Control Problem in an Innite Dimensional Setting152 A.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 A.2 Maximum Prin iple .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156 Bibliography.. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .167 vi List of Tables 5.1 Details of aptured o king events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.2 Hamiltonian Signature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 vii List of Figures 1.1 Mobile robot based experimental platform (Pioneer 3 DX) with two-wheel dierential and aster. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Illustration of s alar shape variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Phase portrait . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3 Null lines of (2.5) with u = 1, ν1 = 0.5, ν2 = 1. . . . . . . . . . . . . . . . . . . . . . . . . 20 2.4 Robot (Pioneer 3 DX) with Kine t mounted, and the orange one used as the bea on. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.5 Implementation results for the limited FOV problem .. . . . . . . . . . . . . . . . 29 2.6 Implementation results for MMC ontrol law. . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.7 Implementation results for TVA ontrol law .. . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1 Interse tion of Hamiltonian and Casimir surfa es . . . . . . . . . . . . . . . . . . . . . 51 3.2 Sample traje tories for dierent values of c, η and s1, s2 . . . . . . . . . . . . . . 53 4.1 Numeri al solution of (4.61) for experiment 1 . . . . . . . . . . . . . . . . . . . . . . . . . 99 4.2 Evolution of x3 for experiment 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 4.3 Numeri al solution of (4.61) for experiment 2 . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4 Evolution of x3 for experiment 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.5 Numeri al solutions of (4.61) for experiment 3 . . . . . . . . . . . . . . . . . . . . . . . . 103 4.6 Numeri al solution of (4.99) and state evolution for experiment 1. . . 112 4.7 State evolution for experiment 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.8 State evolution for experiment 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.9 Numeri al solution of (4.99) and state evolution for experiment 4. . . 116 4.10 Numeri al solution of (4.99) and state evolution for experiment 5. . . 117 5.1 Hamiltonian signatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.2 Optimal energy allo ation for event 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 5.3 Optimal energy allo ation for event 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141 5.4 Optimal energy allo ation for event 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.5 Optimal energy allo ation for event 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 5.6 Optimal energy allo ation for event 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 5.7 Optimal energy allo ation for event 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 5.8 Optimal energy allo ation for event 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 5.9 Optimal energy allo ation for event 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 viii List of Abbreviations and Notations ROS Robot Operating System CB Constant Bearing pursuit FOV Field-of-View MMC Mutual Motion Camouage TVA Topologi al Velo ity Alignment PMP Pontryagin's Maximum Prin iple ri Position of agent i {xi,yi} (planar) Natural Frenet Frame for agent i tr Tra e (of a matrix) 1n n× n identity matrix am z Amplitude z sn z Ellipti sine z, Ja obian ellipti fun tion n z Ellipti osine z, Ja obian ellipti fun tion dn z Delta amplitude z, Ja obian ellipti fun tion TmM Tangent spa e at m to manifold M G A matrix Lie group g Lie algebra of G S1 Cir le group G ∞ 1Loop group C (S ;G) L ∞ 1Loop algebra C (S ; g) H(3) Heisenberg group SE(2) Spe ial Eu lidean group h(3) Lie algebra of H(3) se(2) Lie algebra of SE(2) χ Coupling onstant ix Chapter 1 Introdu tion The last few de ades have witnessed an in rease in resear h eorts towards un- overing me hanisms behind olle tive motion [Nagy et al., 2010; Ballerini et al., 2008a; Cavagna et al., 2010; Inada and Kawa hi, 2002℄ and pursuit behavior [Ol- berg et al., 2000; Mizutani et al., 2003; Ghose et al., 2006; Chiu et al., 2010℄ in nature. Ranging from sh s hools to bird o ks, olle tive behavior is seen abundantly in nature. The on ept of safety in numbers is used in a omplishing variety of goals, from foraging food to evading predators. Re ent improvements in data olle tion and pro essing te hnology has enabled resear hers to study these natural o ks in more detail than ever before [Ballerini et al., 2008a,b℄. The driving question then be omes to answer how lo al intera tions between in- dividual agents in the olle tive give rise to group level ohesion and syn hrony. Although several attempts have been made to understand these behaviors [Cu ker and Smale, 2007; Mora and Bialek, 2011; Bialek et al., 2012; Young et al., 2013; Attanasi et al., 2014℄, the individual level me hanisms responsible for emergen e 1 of olle tive behavior remain mostly elusive to resear hers. It is therefore a sig- ni ant goal of this thesis to pursue su h questions. This thesis is distin tively divided into two parts, where we take two dier- ent approa hes to understand olle tive behavior. The rst approa h is alled a `bottom-up' approa h, i.e. instead of studying the o k as a whole, we on en- trate on dynami s of individual agents and analyze simple intera tion laws among small number of agents. Studying these intera tions are important sin e they an be used as a building blo k for group level motion. In a 1995 paper [Vi sek et al., 1995℄, a novel dis rete time self-driven parti le model was rst introdu ed to ad- dress self-ordered motion in a system of parti les. The on ept of self-steering parti les was developed in the following de ades [Justh and Krishnaprasad, 2004, 2006; Reddy et al., 2006; Mis hiati and Krishnaprasad, 2010, 2012; Galloway et al., 2013℄. We undertake the self-steering parti le model under gyros opi on- trol [Justh and Krishnaprasad, 2003, 2004℄ as the basi model for individuals in the o k. This model des ribes a traje tory of an individual as a urve, des ribed by the natural Frenet frame equations [Bishop, 1975℄ in the Eu lidean spa e; and the driving ontrols are given by speed and urvature of the urve. We show in Chapter 2, 3 that even in the single agent or double agent ase, interesting mo- tion patterns an be synthesized from arefully sele ting these ontrol inputs. The ontrol inputs an be generated from an underlying optimal ontrol problem or by applying some biologi ally plausible feedba k strategies. Parallel to the quest of 2 mathemeti al modeling, some groups in the roboti s ommunity have performed su essful implementation of various ontrol strategies [Thurrowgood et al., 2014; Vásárhelyi et al., 2014℄, and thereby demonstrated the power of a bio-inspired ap- proa h towards synthesizing olle tive motion. Our work is similar in spirit, and provides indoor demonstrations of problems raised in Chapter 2. Some of these problems were on eptualized from a pra ti al perspe tive and arry engineering value. The other approa h to study olle tive behavior is what an be alled as `top- down' view. Instead of spe ifying agent level ontrol laws, the idea is to infer those laws from solving a bigger problem that investigates the o k as a whole. Existing literature employs several methods su h as optimal ontrol [Justh and Krishnaprasad, 2015b,a℄, statisti al physi s [Mora and Bialek, 2011; Bialek et al., 2012℄ et . It is the framework of optimal ontrol [Justh and Krishnaprasad, 2015b℄ that we undertake and extend in this thesis. It has been observed from empir- i al data [Ballerini et al., 2008a℄ that intera tion among starlings in the o k is lo al, i.e. ea h bird intera t with six/seven neighbors during ight. Taking in- spiration from this idea, the entral on ept of [Justh and Krishnaprasad, 2015b℄ is to set up an optimal ontrol problem whi h penalizes ontrols of individual agents oupled with mismat h in ontrol with its `neighbors'. The neighbors are determined by a previously dened intera tion graph. We then let the number of agents in the o k to go to innity in order to propose a ontinuum model for 3 o king. Various ontinuum models have been studied for olle tives [Kudrolli et al., 2008; Topaz et al., 2006; Zhang et al., 2010℄. These models study a set of partial dierential equations that des ribe spatio-temporal evolution of the o k density. Our approa h is dierent in the sense that the system dynami s an be seen as an ordinary dierential equation in an appropriate innite dimensional Lie group setting. The oupling between birds are introdu ed through the mis- mat h term in the ost fun tional. A natural question of ontrollability of su h a system is addressed by using a generalized Chow-Rashevsky theorem for in- nite dimensional systems. This enables us to formulate the underlying optimal ontrol problem in an innite dimensional setting in whi h the usual Pontryagin's maximum prin iple fails in general without further assumptions. In Chapter 4, we invoke a maximum prin iple atered for this spe i setting. A spe i ex- ample of ontinuum of nonholonomi integrators is also studied in detail. This an be viewed as a ontinuum version of single agent Heisenberg ase [Justh and Krishnaprasad, 2016℄. It has been found that optimal ontrol solutions possess a traveling wave hara ter, whi h might enable information transfer in the o k. In addition to the Heisenberg ase, we provide optimal ontrol equations in the ase of a ontinuum o k of planar agents. Syn hronization results and numeri al simulations are presented for both the ases. In Chapter 5, we present another `top-down' approa h to the o king problem. This approa h is data-driven in nature. Kinemati energy modes of European 4 starling o ks are represented on a simplex whi h is then subje ted to des ription as traje tory of some evolutionary game dynami s. Solution of this data-tting problem on the simplex results in ontrol inputs that are interpreted as modulation of tness asso iated with the energy modes. We note that in ontrast to Chapter 4, where the ontrol inputs were individual agent-level (or `low-level') ontrols, the ontrols obtained by this data-driven approa h are o k-level (or `high-level') ontrols. The o k is on eptualized to apply these ontrols to optimally allo ate its kineti energy among dierent modes. 1.1 Mathemati al Ba kground 1.1.1 Self Steering Parti le Model We des ribe the parti le model that is the underlying generative model in all our subsequent analysis throughout this thesis. The traje tory of a single agent an be des ribed by a fun tion r : [0, T ] → R3, for some T > 0. We assume r(t) to be a regular urve, i.e. ṙ(t) =6 0, ∀t ∈ [0, T ]. Let s be the ar length parameter, ∫ t i.e. s(t) = ‖ṙ(σ)‖ dσ. Under the regularity assumption, s(t) is monotoni ally 0 in reasing and invertible fun tion of time. We an then reparametrize the urve r(t) by the ar length parameter s and the evolution equations an be expressed in terms of well known Fernet-Serret frames. However, this way of representation requires thri e dierentiability of the urve and need the urvature of the urve to be stri tly positive. To over ome these di ulties, we take an alternate approa h 5 for framing the urve, known as the Natural Frenet frame [Bishop, 1975℄. This approa h requires only twi e dierentiability and is well dened even when the se ond derivative vanishes. In 3D, to any point on the urve r(t), we atta h an orthonormal moving frame {x(t),y(t), z(t)}. The unit ve tor x(t) is tangent to the urve and points toward the heading of an individual. The unit ve tors {y(t), z(t)} are hosen in the plane normal to x(t). The evolution of these ve tors are given by the frame equations, ṙ(t) = ν(t)x(t) ẋ(t) = ν(t)(u(t)y(t) + v(t)z(t)) (1.1) ẏ(t) = −ν(t)u(t)x(t) ṙ(t) = −ν(t)v(t)x(t), where ν(t) is the speed (‖ṙ(t)‖) and (u(t), v(t)) are alled natural urvatures of the traje tory [Justh and Krishnaprasad, 2005℄. In a planar setting, we have the frame {x(t),y(t)} and the evolution equations are written as, ṙ(t) = ν(t)x(t) ẋ(t) = ν(t)u(t)y(t) (1.2) ẏ(t) = −ν(t)u(t)x(t). We an therefore treat ν and u variables as ontrol inputs to steer the individual on the plane, ν as the velo ity input and u as the urvature ontrol input. 6 Figure 1.1: Mobile robot based experimental platform (Pioneer 3 DX) with two- wheel dierential and aster. 1.2 Experimental Setup We provide a omprehensive des ription of the laboratory set up in the Intelli- gent Servosystems Lab, University of Maryland. All the laboratory experiments presented in this thesis are done under this setup. Our experimental test-bed is omprised of Pioneer 3 DX wheeled robots from Adept MobileRobots [Pioneer℄. These ompa t, dierential-drive mobile robots are equipped with reversible DC motors, high-resolution motion en oders and 19 m wheels, and the onboard om- putation is done via a 32-bit Renesas SH2-7144 RISC mi ropro essor, in luding the P3-SH mi ro ontroller with ARCOS. The sensors on the robot in lude eight forward-fa ing ultrasoni (sonar) sensors. ARIA [ROS-ARIA℄ provides an inter- fa e for ontrolling and re eiving data from the robot, and ommuni ation with the robot for sending ontrol ommands (forward velo ity and turning rate) is done via 802.11-b/g/n networking. The width of the robot is 380 mm and it has a swing radius of 260 mm. 7 Algorithm implementation (i.e, feedba k law omputation) has been done in C++ using ROS [ROS℄, along with ROS-ARIA [ROS-ARIA℄, as the interfa ing roboti s middleware. The experiments have been arried out in a laboratory en- vironment equipped with a sub-millimeter a urate Vi on motion apture system [Vi on℄. We use a Dell workstation to run ROS, and this omputer is onne ted to the Vi on server via a dedi ated Ethernet onne tion. The Vi on system aptures the motion of the robots and sends out the position and heading data to the omputer running ROS. The ontrol law program listens to this data, and transmits the individual velo ities and turning rates. Both of these operations are arried out at a frequen y of 25 Hz. The ontrol law program omputes the ontrols a ording to the strategy that is spe i to the problem onsidered. The omputed velo ity and urvature ontrol variables ν(t), u(t) an be translated to the turning rate ω(t) (in degrees/se ) as: ( ) 180 ω(t) = ν(t)u(t). (1.3) π 8 Part I Synthesis of Colle tive Motion: Bottom-up Approa hes 9 Chapter 2 Feedba k Laws for Colle tive Motion 2.1 Introdu tion Colle tive motion plays a ru ial role in modern day roboti s and engineering. It is be oming ommonpla e for a group of unmanned, remote ontrolled vehi les to be deployed to a omplish goals ranging from sear h and res ue to surveil- lan e. For the swarm of robots to fun tion in a harmonious manner, it is very important to ontrol them arefully. Natural olle tives are indeed an inspira- tion in this endeavor. On the other hand, a thorough study of those olle tives remain in omplete without understanding agent level intera tion laws. In this hapter therefore, we will build models for olle tive motion from `bottom-up', i.e. from individual level ontrol strategy to o k level syn hrony through in- tera tion among agents. This hapter presents a range of theoreti al results as well as laboratory demonstrations of ontrol laws that we propose or has been proposed before. This hapter has two main ontributory se tions, most of whi h 10 are taken verbatim from their respe tive publi ations. We start with a problem that studies a parti ular dyadi intera tion under the setting of a pursuit strategy alled onstant bearing pursuit [Halder et al., 2016℄. We gain interesting insight from this problem that onne ts to the Kepler two body problem. The obtained result of this problem is then used to solve a problem arising from a pra ti al robot maneuvering s enario. The last problem is purely experimental [Halder and Dey, 2015℄ whi h demonstrate another dyadi intera tion strategy potentially useful for surveillan e, and a o king strategy involving many agents. 2.2 Steering for Bea on Pursuit under Limited Sens- ing In this se tion, we will try to understand simple dyadi pursuit strategies (i.e. strategies based on pairwise intera tions), and exploit them as building blo ks for synthesis of omplex motion patterns for olle tives. In [Galloway et al., 2009, 2013℄, using symmetry prin iples and nonlinear dynami s, a spe i strat- egy, known as onstant bearing y li pursuit, is shown to produ e a ri h variety of behaviors for appropriate hoi es of parameters (bearing angles). In [Justh and Krishnaprasad, 2006℄ a biologi ally plausible feedba k ontrol law is inves- tigated that exe utes motion amouage, a type of stealthy pursuit asso iated with visually-guided ight in inse ts (e.g. hoveries and dragonies). Stealth arises from nulling opti ow in the visual eld of the target of pursuit, thereby 11 in reasing the han e of su ess in prey apture or territorial battle against a on- spe i . This type of dyadi intera tion is also exploitable in oordinated motion, for instan e see [Mis hiati and Krishnaprasad, 2010℄. The present work is similar to [Justh and Krishnaprasad, 2004℄ in motivation. We onsider a problem of two agents moving in a plane with onstant (not ne - essarily identi al) speeds and, one of them is free i.e. it assumes any open loop steering ( urvature) ontrol, while the other pursues it. The free agent may be onstrued as a bea on and the pursuer's task is to rea h a safe vi inity of the bea on and ir ulate around it. In the interesting ase when the bea on is sta- tionary, but the pursuer has a sensor with limited eld of view (FOV) to dete t the bea on, the ir ling law proposed in [Justh and Krishnaprasad, 2004℄ may be foiled. One goal of this work is to devise a prin ipled approa h ( ontrol algorithm) for this problem that opes with sensor limitation. We do this via a two-step pro- ess. We rst analyze a slightly dierent problem of tra king a (slowly) moving bea on assuming that: (a) the bea on tra k is of onstant urvature (i.e. on a straight line or on a ir le); and (b) the sensor on-board the pursuer has no FOV limitation. For the hoi e of a onstant bearing pursuit feedba k ontrol law, one obtains a ri h dynami s. The phase portrait in turn suggests the se ond step  a feedba k law modi ation that is appli able to the setting of stationary bea on, and limited FOV. In this ase, one needs an additional ingredient  an estimator using odometry to tra k the bea on when it has fallen out of the FOV. The idea here is to use the odometry-based estimate in the feedba k law as if it is exa t 12 PSfrag repla ements y1 x1 x2 κ1 y2 κ2 r1 r12 = r1 − r2 r2 Figure 2.1: Illustration of s alar shape variables (ρ, κ1, κ2) used to parametrize the shape spa e. (a type of ertainty equivalen e) and use dire t observation of the bea on when it is re- aptured in the FOV. Su h an estimator is used in an implementation of bea on tra king in a laboratory test-bed with a range amera (Kine t [Kine t℄) as the sensor mounted on a mobile ground robot. A high pre ision (marker-based) motion apture system (Vi on [Vi on℄) is used to determine ground truth and analyze the performan e. In addition, the un onstrained tra king of the moving bea on problem is revisited and it is shown that the resulting dynami s an be identied with motion of a harged parti le in an ele tromagneti eld. Moreover, at a parti ular value of the bea on urvature, the ombined dynami s is exa tly same as the Kepler problem of two bodies. 2.2.1 Tra king a Moving Bea on Let us onsider two agents moving on a plane, ea h abiding the self steering par- ti le equations of motion. We assume that both their speeds ν1, ν2 are onstants. It is possible to represent the dynami s of the system of two agents by the help 13 of s alar shape variables ρ, κ1, κ2 (Fig. 2.1) as ρ̇ = −ν1 cosκ1 − ν2 cosκ2 − 1κ̇1 = ν1u1 + (ν1 sin κ1 + ν2 sin κ2) (2.1) ρ − 1κ̇2 = ν2u2 + (ν1 sin κ1 + ν2 sin κ2). ρ This is rather a straightforward al ulation. We view agent 1 as a slowly moving bea on to whi h agent 2 pays attention. Let us make the following assumption. A-1: The speed of agent 1 is less than the speed of agent 2, i.e. ν1 < ν2. We pi k the feedba k ontrol law for agent 2 as follows: ( ) ( ) − r21 1 r21u2 = µ̃ R(α)y2 · | | − ⊥ r ν | · ṙr | |r | 21 , (2.2)21 2 21 21 ⊥ for some µ̃ > 0. Here we denote a = R(π/2)a, for any ve tor a in the plane of motion, R(·) is the planar rotation matrix. Note that this ontrol law is a standard onstant bearing (CB) pursuit law [Galloway et al., 2013℄ with parameter α. The feedba k ontrol law an be expressed in terms of the s alar shape variables as 1 u2 = µ̃ sin(κ2 − α) + (ν1 sin κ1 + ν2 sin κ2) ν2ρ The losed loop dynami s of (2.1) then takes the form ρ̇ = −ν1 cosκ1 − ν2 cos κ2 1 κ̇1 = −ν1u1 + (ν1 sin κ1 + ν2 sin κ2) (2.3) ρ κ̇2 = −µ̃ν2 sin(κ2 − α). A fundamental result [Galloway et al., 2013℄ for the CB strategy tells us that under the a tion of the ontrol law (2.2), the manifold Mα + 1 1CB = {(ρ, κ1, κ2) ∈ R × S × S : κ2 = α} 14 is an attra tive invariant manifold for all initial onditions ex ept κ2(0) = α + π. The invarian e follows dire tly from the losed loop dynami s (2.3). The attra tiveness an be proved by dening Λ(t) = − cos(κ2(t)− α). Thus Λ(0) =6 1 implies Λ(t) → −1 as t → ∞ or equivalently κ2 onverges to α. 2.2.2 Dynami s Restri ted to the Invariant Manifold At this stage, we are ready to make another assumption: A-2: We onsider the urvature of the bea on to be onstant, i.e. u1 = u, for some u ∈ R onstant. Now we fo us our analysis on the dynami s on the invariant manifold ( alled M + 1Shape = R × S ) whi h may be expressed as ρ̇ = −ν1 cos κ1 − ν2 cosα (2.4) 1 κ̇1 = −ν1u+ (ν1 sin κ1 + ν2 sinα). ρ It is of interest to hara terize the solutions of the restri ted dynami s (2.4) on the invariant manifold. Note that given ν1, ν2, u and α, (2.4) might have at most (ρ∗, κ∗ ν sinκ∗+sinα two equilibrium points 1), with cosκ ∗ = − cosα ρ∗ 11 and = , whereν νu ν = ν1we denote < 1. Existen e of su h equilibrium points is guaranteed if ν ≥ ν2 | cosα| and (ν sin κ∗1 + sinα)u > 0. Linearizing (2.4) around su h an equilibrium point gives the Ja obian matrix   0 ν sin κ∗  1 ν  2 ,   −νu − cosα ρ∗ ρ∗ 15 with asso iated Eigenvalues ( ) √ − cosα± cos2 α− 4ν2ρ∗u sinκ∗ λ = ν 12 . 2ρ∗ Depending on α, following ases will arise I. α ∈ (−π,−π/2) ∪ (π/2, π] In this region, cosα < 0, whi h makes the equilibrium points unstable. So all traje tories tend to blow up in (ρ, κ1) plane. II. α ∈ (−π/2, π/2) Here, cosα > 0 ∗, then (lo ally) stable equilibrium exists if u sin κ1 > 0, otherwise (ρ∗, κ∗1) is unstable whi h leads to eventual ollision. We note that [Davis, 1962℄ (pages 119-125) studies the same problem with α = 0. For the α = 0 ase, the existen e onditions of equilibrium read ν ≥ 1 and u sinκ∗1 > 0, whi h in turn means we will have a stable equilibrium only when ν ≥ 1. The urrent problem an be viewed as a generalization of that onsidered in [Davis, 1962℄. III. α ∈ {π/2,−π/2} In this ase, however, the dynami s (2.4) produ es a ri h behavior whi h we analyze next. We only provide the analysis for α = π/2 ase, the other ase being analogous. 16 π 2.2.3 Spe ial Case: α = 2 We rewrite (2.4) in this parti ular ase, ρ̇ = −ν1 cosκ1 (2.5) 1 κ̇1 = −ν1u+ (ν2 + ν1 sin κ1). ρ All the traje tories of (2.5) are losed and we will prove this in the same way as exploited in [Mis hiati and Krishnaprasad, 2012℄. We rst introdu e the following denitions and a theorem due to Birkho. Denition 2.1 (Involution). A dieomorphism F : M → M from a manifold M onto itself is said to be an involution if F 6= idM , the identity dieomorphism 2 and F = idM , i.e. F (F (m)) = m, ∀m ∈ M . Denition 2.2 (F-reversibility). A ve tor eld X dened over a manifold M is said to be F-reversible if there exists an involution F su h that F∗(X) = −X, i.e. F maps orbits of X to orbits of X , reversing the time parametrization. Here (F∗(X))(m) = (DF )F−1(m)X(F −1(m)), ∀m ∈ M is the push-forward of F . We all F the reverser of X . Theorem 2.2.1 (G.D. Birkho, [Birkho, 1915℄). Let X be a F-reversible ve tor eld on M and ΣF the xed-point set of the reverser F . If an orbit of X through a point of ΣF interse ts ΣF at another point, then it is periodi . See [Mis hiati and Krishnaprasad, 2012℄ for a detailed proof of this theorem. Based on these denitions and the theorem, we propose the following. 17 90 4 90 8 90 4 120 60 120 60 120 60 3 6 3 150 2 30 150 4 30 150 2 30 1 2 1 180 0 180 0 180 0 210 330 210 330 210 330 240 300 240 300 240 300 270 270 270 (a) u = −2 (b) u = 0 ( ) u = +2 Figure 2.2: Phase portrait (polar plot) of the system dynami s restri ted to the invariant manifold (2.5) with dierent values of u, keeping ν1 = 0.5, ν2 = 1 xed for all three ases. Theorem 2.2.2. (i) The quantity, 1 E(ρ, κ1) = ρ(ν2 + ν1 sin κ1)− ν1uρ2 = E(ρ(0), κ1(0)) (2.6) 2 is onserved along any traje tory of (2.5). (ii) Every solution of (2.5) is periodi . Proof. (i) Denote, χ = ρ(ν2 + ν1 sin κ1), then dχ = ρ̇(ν2 + ν1 sin κ1) + (ρν1 cosκ1)κ̇1 dt dρ = ν1uρ dt =⇒ 1χ(t) = ν1uρ2(t) + c, 2 where c = χ(0)− ν1uρ2(0)/2 = onstant, whi h, in turn implies 1 E(ρ, κ1) = ρ(ν2 + ν1 sin κ1)− ν1uρ2 = E(ρ(0), κ1(0)). 2 18 (ii) Step-1: Ve tor eld dened by (2.5) is F-reversible with reverser F (ρ, κ1) = (ρ, π − κ1). 2 Clearly, F is an involution sin e F (ρ, κ1) = (ρ, κ1). Next, (F∗(X))(ρ, κ1) = (DF )F−1(ρ,κ )X(F −1(ρ, κ1))1   1 0   =  X(ρ, π − κ1)   0 −1 = −X(ρ, κ1). Hen e, X is F-reversible. π Step-2: Fixed point set of F is given by ΣF = {(ρ, κ1) : ρ > 0, κ1 = ± }. So2 π every orbit of (2.5) rossing κ1 = ± line twi e is periodi . Now, depending on2 the value of u, dierent ases will arise. (a) u ≤ 0 : In this ase, we note that the assumption ν2 > ν1 is su ient to guarantee monotoni ity of κ̇1, in parti ular κ̇1 > 0 for all time. Hen e, any traje tory originating from any point on the κ1 = ±π/2 line (ex luding the origin) will travel ounter lo kwise until it hits the line again when κ1 gets in remented by an amount of π radian (see Fig. 2.2a, 2.2b). Note that the onserved energy, E prohibits any traje tory that starts with positive energy to go to the origin (with zero energy). (b) u > 0 : Be ause κ̇1 an assume any sign under this ase, we need a more serious argument for this ase. To determine the null line κ̇1 = 0, we ompute 19 PSfrag repla ements 90 4 120 60 A 3 150 2 30 1 III IV I II 180 0 O B NK 210 330 240 300 270 Figure 2.3: Null lines of (2.5) with u = 1, ν1 = 0.5, ν2 = 1. ρ = ν2+ν1 sinκ1 and ν2 > ν1 ensures existen e of a valid ρ for ea h value of κν u 1. The1 ρ̇ null line is simply the κ1 = ±π/2 line. This is illustrated in Fig. 2.3. It is ( ) ν2+ν1 π immediate that in this ase we have two equilibrium points, A = ,+ and ν1u 2 ( ) = ν2−ν1 ,−πB . The traje tories starting from either of those points are learly ν1u 2 periodi . Depending on the null lines, the whole spa e an be divided into four regions as shown in Fig. 2.3 and those regions are hara terized as Region I : ρ̇ < 0, κ̇1 < 0, Region III : ρ̇ > 0, κ̇1 > 0, Region II : ρ̇ > 0, κ̇1 < 0, Region IV : ρ̇ < 0, κ̇1 > 0. Now imagine traje tories starting on the line segment OA, ex luding both points. Sin e κ̇1 < 0, they will move into region III, whi h an produ e two out omes: (b1) It hits the OB line (ex luding both points). 20 (b2) It exits region III through NK and goes into region II. Now the onstan y of E gives the following observation, if a traje tory of (2.5) rosses NK at (1) (2) a ertain angle κ1 , then it must ross it again at a dierent angle κ1 and these angles are symmetri about κ1 = π/2. This means traje tories must not enter region III from region II. This leads all the traje tories in region II to hit the κ1 = +π/2 line beyond point A. Next, onsider traje tories starting at the boundary of region I and II with κ1 = π/2. κ̇1 < 0 gives rise to lo kwise motion into region I. Again, we need to analyze two s enarios: (b3) The traje tories rea h boundary between region I and II with κ1 = −π/2. (b4) They enter region IV through NK. Similar argument as in ase (b2) an be employed to prove they must rea h boundary between region III and IV with κ1 = π/2. Now, traje tories starting on OB line (ex luding both points) must move in region IV and hen e must hit OA line (ex luding both points). Finally, traje tories starting on the boundary between region I and II with κ1 = −π/2 must go into region II and must eventually hit the boundary between region I and II with κ1 = π/2 following a lo kwise path (and without entering region III). This ompletes the proof.  Remark 2.1. In a spe ial ase u = 0, (2.6) reads exa tly as the polar equation of 1 an ellipse = c̃(1 + ν cos θ), where the origin is pla ed at one fo us of the ellipse, ρ 21 the angle θ is measured from the origin with respe t to the major axis of the ellipse ν = ν1and is the e entri ity of the ellipse (re all A-1). It is of interest to note ν2 that the same ellipti equation omes from the analysis of the Kepler two-body problem [Goldstein et al., 2001℄. A result of the dynami s (2.5) regarding time period is immediate in the light of Proposition 2.2.2. Corollary 2.2.1. Every orbit of (2.4) has a period ∫ ρmax dρ T = 2 √ , (2.7) ρ 2 1min ν1 − ( ν1uρ+ E0 − ν 22 ρ 2) where ρmin and ρmax are solutions of the pair of equations 1 ρ(ν ± ν )− ν uρ22 1 1 = E(ρ(0), κ1(0)) =: E0 (2.8) 2 2πν2E0 In parti ular, for the ase u = 0, the time period be omes, T = 3 . (ν2 2 22−ν1 ) Remark 2.2. Note that for an admissible value of E0, the pair of equations (2.8) has only two solutions. For the spe ial ase, u = 0, we know the losed loop traje tories are des ribed by the ellipses ρ = E0/(ν2 + ν1 sin κ1) with semi-major ( ) a = 1axis, (ρmin + ρ E0 1 1 E0ν2 2 max ) = + 2 ν +ν ν − =ν ν2−ν2 . Then, from Corollary2 1 2 1 2 1 ( ) 2 4π2 3 2.2.1, we nd T = a . ν2E0 Remark 2.3. The ondition ν2 > ν1 is ne essary for the existen e of periodi orbits for the ase u = 0 while it is merely a su ient ondition for other values of u. 22 2.2.4 The Limited Field of View (FOV) Problem In this se tion, we will des ribe the problem of tra king a stationary bea on by a ontrolled agent equipped with a sensor (e.g. a depth amera like Kine t) with limited eld of view. As opposed to the problem dis ussed in the previous se tion, this problem is inspired by an implementation point of view. The ba kbone model (2.1) of the system stays the same. Agent 2 is supposed to sense the position of the bea on relative to its own position and use the sensed quantities to determine the ontrol a tion. Using the shape variables, it has a ess to the pair (ρ, κ2) (refer to Fig. 2.1), with the limitation that |κ2| ≤ κmax < π/2, whi h we all the eld of view onstraint. Although various feedba k ontrol laws have been proposed (for e.g. [Justh and Krishnaprasad, 2004℄) to en ir le a stationary bea on, permanent loss of the target (bea on) from the eld of view annot be avoided by those laws. More pre isely, the limited FOV problem boils down to en ir ling the bea on while being able to sense it (at least) periodi ally. Putting bea on speed, ν1 = 0 in (2.1) and ignoring κ1 dynami s, the equivalent shape spa e equations an be redu ed to ρ̇ = −ν2 cosκ2 (2.9) − 1κ̇2 = ν2u2 + ν2 sin κ2. ρ Remark 2.4. From (2.9), it is guaranteed that under the eld of view onstraint, the attempt of en ir ling the bea on would eventually make the bea on perma- nently invisible (as long as a ir ular orbit around the bea on is onsidered). Moreover, from ρ dynami s, meeting the onstraint |κ2| ≤ κmax < π/2 for all time 23 will lead to denite ollision. On e the bea on goes out of the FOV, the only hoi e for the vehi le would be to e iently estimate the position of the bea on and apply the ontrol based on those estimates. A epting the fa t mentioned in remark 2.4, one an only try to design the ontrol u2 in su h a way that the ontrol law provides some promise to bring the bea on ba k in the eld of view after losing it. The following proposition is meant to serve that purpose. Proposition 2.2.1. The feedba k ontrol law given by uFOV µ 2 = u0 − , u0 ≤ 0, µ > ν2, (2.10)ρν2 guarantees the periodi return of the bea on to the eld of view under ideal esti- mates. Proof. With the feedba k ontrol (2.10), the losed loop system be omes, ρ̇ = −ν2 cosκ2 (2.11) κ̇2 = − 1 ν2u0 + (µ+ ν2 sin κ2). ρ Noti ing that (2.11) is equivalent to (2.5), the laim follows dire tly from Theorem 2.2.2. From the polar phase portrait (Fig. 2.2), we see that the ondition u0 ≤ 0 ◦ is required for the angle variable κ2 to go through a full 360 rotation whi h is essential in order to bring the bea on ba k in the FOV.  As we will dis over next, the ondition on the parameter u0 an be relaxed to in lude positive values as well. From Theorem 2.2.2, the ondition µ > ν2 is 24 only su ient for any value of u0. It be omes ne essary for the parti ular ase of u0 = 0. Although the feedba k ontrol law (2.10) produ es all periodi orbits and it inherently takes are of ollision avoidan e problem (see Fig. 2.2), it la ks in the freedom of driving the vehi le to a parti ular desired orbit ( f. the ir ular orbits with desired radius as in [Justh and Krishnaprasad, 2004℄). Moreover, the periodi orbits are not orbitally asymptoti ally stable whi h makes them sus eptible to disturban es. To over ome these short omings, we propose the following. Proposition 2.2.2. Let Ed denote the admissible value of the desired energy, i.e. there exists a periodi orbit with E(ρ, κ2) = Ed. Here E(ρ, κ2) is as in (2.6) with the pair (ν1, ν2) interpreted as (ν2, µ) in present ontext. Then the ontrol law u = uFOV + uAD µ 2 2 2 = u0 − + kd(E(ρ, κ2)−Ed) cosκ2, (2.12)ρν2 with kd > 0 makes the orbit with energy Ed asymptoti ally stable with region of attra tion given by MShape \ {(ρ, κ2) : cosκ (µ+ν2 sinκ2)2 = 0, ρ = , u0 > 0}, whereν2u0 M = R+ 1Shape × S . Proof. Note that the traje tories of (2.9) with ontrol (2.12) will no longer be peri- AD odi be ause of the in lusion of the extra u2 term. Sin e E(ρ, κ2) is a ontinuous 2 fun tion of both ρ and κ2, it su es to prove that the quantity (E(ρ, κ2) − Ed) is monotoni ally de reasing. We obtain, d (E(ρ, κ2)−E 2d) = −2kdρν22(E(ρ, κ )− E )22 d cos2 κ2.dt 25 Clearly, with kd > 0, d(E(ρ, κ2)−Ed)2/dt ≤ 0, for all (ρ, κ2) su h that ρ > 0. The 2 largest invariant subset of {(ρ, κ2) : d(E(ρ, κ2)−Ed) /dt = 0} is indeed {(ρ, κ2) : E(ρ, κ2) = Ed}, whi h, in turn establishes the statement of the proposition.  In the light of Theorem 2.2.2, we see that the restri tion u0 ≤ 0 in Proposition 2.2.1 an be relaxed. In parti ular, for u0 > 0, one only has to hoose Ed su h that k2 ompletes full 360 ◦ rotation (for e.g. one might pi k Ed = E(ρ,−π/2), µ−ν2 with ρ > , see Fig. 2.2( )). ν2u0 2.2.5 Implementation In this laboratory implementation, we hose to use the newest Kine t model, whi h was reated for Mi rosoft's Xbox One. The Kine t primarily fun tions as a motion-sensing input devi e, enabling players to intera t with video games in ex iting ways. To a omplish this, the devi e is equipped with several sensors in luding an RGB sensor, 3D Depth Sensor, as well as Multi-array Mi rophones. The Kine t's RGB sensor has a 70.6 degree horizontal eld of view, and a 60 degree verti al eld of view (see Fig. 2.4). The Kine t operates at a rate of 30 Hz, and has an ee tive range between 0.5 meters, and 4.5 meters where a ura y is reliable. Despite it's original use ase as a video game ontroller, the Kine t has been studied re ently as a sensor for many roboti s appli ations, in luding autonomous vehi les [Oliver et al., 2012℄ and health are [Nghiem et al., 2012℄. In this experiment, the Kine t RGB amera a ts as a primary sensor for de- termining the distan e and relative heading of the bea on (i.e. ρ, κ2), whi h in 26 Figure 2.4: Robot (Pioneer 3 DX) with Kine t mounted, and the orange one used as the bea on. our experiment was an orange one. OpenCV [OpenCV℄ is used to perform a simple blob dete tion algorithm that al ulates the entroid of the one in pixel oordinates, and then uses the Kine t's oordinate mapping feature to transform the result into physi al, or amera spa e. To take advantage of these API features, we mount a laptop running the Windows operating system onto the robot, and utilize a ustom TCP/IP server to stream the oordinates ba k to the robot on- trol station. The ontrol station is a Dell omputer running ROS [ROS℄, and the algorithm implementation is done using the MATLAB ROS toolbox [MATLAB℄. Finally, the Vi on motion apture system [Vi on℄ is used to tra k the motion of the robot and bea on in the lab oordinate spa e to obtain ground truth results of the implementation. 2.2.5.1 Estimation In order to su essfully implement the proposed ontrol law, the robot has to be able to e iently determine the bea on position relative to its own position during 27 the periods of time when the bea on is not in FOV. Equivalently, the estimation problem then is to integrate the losed loop shape dynami s ((2.9) with ontrol (2.12)) given some initial ondition, whi h in this ase would be the last known (ρ, κ2) value before the robot loses sight of the bea on. Sin e there is a onserved quantity asso iated with (2.11) (Proposition 2.2.2), a mid-point based update rule performs better than the naive Euler rule [Austin et al., 1993℄. Denoting the estimate of (ρ, κ2) by (ρ̂, κ̂2) and the dis rete time step by ∆t, the update rule may be impli itly expressed as follows: ( ρ̂n+1 − ρ̂n κ̂n + κ̂n+1) = −ν 2 22 cos ∆t 2 (2.13) ( κ̂n+1 − κ̂n 2ν κ̂n + κ̂n+1)2 2 = − 2ν2un2 + sin 2 2 ,∆t ρ̂n+1 + ρ̂n 2 n n n 0 0 where u2 = u2(ρ̂ , κ̂2 ) as in (2.12) and (ρ̂ , κ̂2) is the last su essful measurement of the bea on position. We then solve the nonlinear equations (2.13) numeri ally (using MATLAB's fsolve) to produ e the ne essary estimates whenever the bea on is not in the eld of view of the sensor. This pro edure an be summarized in Algorithm 1. 2.2.5.2 Experimental Results To demonstrate our solution to the limited eld of view problem, we onstru ted an experiment for whi h the robot sees the orange one and attempts to en ir- le it using the des ribed ontrol me hanisms. The result is a traje tory that periodi ally brings the one ba k in its FOV so that the robot an fulll its net 28 Algorithm 1: Steering Law Computation for Limited FOV Problem Data: ρ : measured distan e to bea on, κ2 : measured angle to bea on; ρ̂ : estimate of ρ, κ̂2 : estimate of κ2 Parameters: u0, µ, kd, Ed begin while not stopped do if bea on visible then Compute u2 = u2(ρ, κ2) using equation (2.12) else 0 0 last last Initialize: (ρ̂ , κ̂2) ←− (ρ , κ2 ) Cal ulate (ρ̂, κ̂2) from (2.13) Determine u2 = u2 (ρ̂, κ̂2) using (2.12) Expe rimental 90 2000 Ideal Beacon 2000 120 60Robot 1500 ) 1500 150 1000 30 m m 1000 500 n i ( 500 180 0 PSfrag repla ements 0 −500 210 330 −1000 240 300 −2000 −1000 0 1000 2000 x 270(in mm) (a) Robot Traje tory (b) Phase Portrait Figure 2.5: Implementation results as re orded by Vi on (with ν2 = 200 mm/se , u −10 = 0 mm , µ = 1000 mm/se , kd = 5 × 10−9 −3 6mm se and Ed = 1.44 × 10 2 mm /se ). 29 y en ir ling of the bea on. The ground truth data was obtained using Vi on. We ran the experiment using a robot speed ν2 = 200 mm/se , and the parameters were u0 = 0 −1 mm , µ = 5ν −9 −32 = 1000 mm/se , kd = 5 × 10 mm se and Ed was taken to be energy orresponding to the orbit whi h maintains the minimum 6 2 distan e of 1200 mm from the bea on, in other words Ed = 1.44× 10 mm /se . Implementation results are shown in Fig. 2.5. The ground truth polar plot an be seen in omparison to the desired ellipse (sin e u0 was taken to be 0) in Fig. 2.5b. The mid-point rule estimation method results in a robust ontroller that a hieves the desired traje tory although it is slightly loser to the bea on than the theory predi ts. The error between these two orbits is observed (∼ 200 mm) to be within the size of the robot (∼ 400 mm). 2.2.6 Asso iated Lagrangian Here we will re-visit the problem of tra king a moving bea on as onsidered in Se tion 2.2.1. It is of spe i interest to ask whether the system dynami s admits some Lagrangian formulation. Without loss of generality, at this stage we take ν2 = 1 and denote ν1 = ν (note that A-1 translates to ν < 1). Writing r = r1−r2, we see that on the invariant manifold, the feedba k ontrol law (2.2) takes the form r · ṙ⊥ u2 = − | .r|2 30 Now it is a straightforward exer ise to see that the baseline ve tor satises the following se ond order ODE: r̈ = ν2uy1 − u2y2 ( )( ) 2 r · ṙ⊥ r= ν uy1 + | . (2.14)r|2 |r| Note that here we used the fa t that on the invariant manifold the bearing angle of the se ond agent to the bea on, κ2 = π/2 and hen e the y2 ve tor is aligned ( ) r with the baseline ve tor r, or in other words y2 = | | .r Proposition 2.2.3. On every level set of E, (E(ρ, κ1) = E0, as in (2.6)) the two dimensional system (2.14) is a tually the Euler-Lagrange equation of the La- grangian fun tion (of the type kineti energy−potential energy) L(r, ṙ) = K(r, ṙ)− V (r, ṙ) ( ) 1 | |2 − −E0 − 1= ṙ | | νu|r| − A(r) · ṙ , (2.15)2 r 2 1 ⊥ where A is dened as A(r) := − νur . 2 Proof. Note that the quantity, E0 = ρ(1 + ν sin κ) − 1νuρ2 is onserved and2 r · ṙ⊥ = −ρ(1 + ν sin κ). From here, we an rewrite (2.14) as ( )( ) r̈ = νu(νy1 − E0 1 r y2) + −| | + νur 2 2 |r| ( ) E0 1 = (ṙ×B)−∇ −| | − νu|r|r 2 = (ṙ×B) + E(r), (2.16) where we introdu e B := −νuẑ, a stati magneti eld in the dire tion per- pendi ular to the plane of motion of the agents, ẑ being the unit ve tor in 31 that dire tion. Also, we all E(r) := −∇Φ(r), the ele tri eld and Φ(r) := ( ) − E0 + 1| | νu|r| , the orresponding ele trostati potential.r 2 Now, the equation (2.16) resembles with that of the equation of motion of a harged parti le in an ele tromagneti eld and one an nd obvious similarity of the right hand side with that of Lorentz for e law. Then a standard result [Goldstein et al., 2001℄ in the theory of ele tromagnetism gives the Lagrangian formulation of (2.16). Sin e ∇ · B = 0, B an be written as url of a magneti ve tor potential A(t, r), i.e. B = ∇×A. Also, the ele tri eld E an be writ- ten as E = −∇Φ(r) − ∂A . Mat hing this with (2.16), we an see that A is a ∂t ve tor valued fun tion of r only. It is a straightforward exer ise to show that A = −1(r × B) = −1νur⊥ satises B = ∇ × A. Then the Lagrangian whi h 2 2 generates (2.16) is given by L(r, ṙ) = K(r, ṙ)− V (r, ṙ) 1 = |ṙ|2 − (Φ(r)−A(r) · ṙ) 2 1 = |ṙ|2 E0 1 1+ | | + νu|r| − νur ⊥ · ṙ 2 r 2 2  Remark 2.5. In essen e, Proposition 2.2.3 reveals that with open loop onstant urvature ontrol of one agent and with feedba k ontrol (2.2), namely onstant bearing pursuit law with parameter α = π/2, of the other, (on every level set of an invariant manifold) the oupled system behaves exa tly the same way as a harged parti le in a stati ele tromagneti eld. Moreover, the magneti eld, B is dependent on the onstant urvature (u) of the rst agent. Thus u = 0 implies 32 the absen e of the magneti eld and the agents are subje t to the ele trostati eld only. The eld E has familiar inverse square form and hen e the results for the spe ial ase u = 0 agree with Kepler's laws whi h are obtained from a two-body problem subje t to Newtonian gravitational for e. 2.3 Biomimeti Algorithms for Coordinated Mo- tion In this se tion, we will report implementation of two feedba k ontrol strategies on our laboratory test-bed. The rst of these two strategies is alled mutual mo- tion amouage (MMC) [Mis hiati and Krishnaprasad, 2012℄. Existing literature on dragonies [Corbet, 1999℄ provides qualitative analysis of territorial battles, wherein the traje tories display spiraling motion onsistent with the theoreti al predi tions [Mis hiati and Krishnaprasad, 2011℄. This parti ular bio-inspired on- trol algorithm inherits an appealing overage property through the me hanism of spa e lling urves, and our implementations are able to reprodu e overage pat- terns similar to the predi ted ones. Although there has been a long history of ontrol algorithms for o king, al- most every model of olle tive motion predi ts diusive transport of information. But, ontrary to the existing models, re ent ndings [Attanasi et al., 2014℄ from starling o ks suggest that dire tional information within a o k propagates with an almost onstant speed, and this linear growth of information an be explained 33 by models with wave-like aspe ts. In [Dey, 2015; Halder and Dey, 2015℄ a on- trol strategy alled topologi al velo ity alignment (TVA) was introdu ed, whi h onforms to this riterion and an explain how information about lo al neighbors an inuen e the agents in a o k to align their headings in a single ommon dire tion. Hen e it seems reasonable to use TVA strategy for olle tive motion synthesis. Furthermore, our implementation results in real robots have shown that redu tion in neighborhood size and external perturbation (similar to preda- tor atta k) an split a o k into smaller subgroups. 2.3.1 Mutual Motion Camouage (MMC) Here we onsider the mutual intera tion between two agents ea h applying the same pursuit law, while per eiving the other one as a target. As the dynami s of MMC in a planar setting has been studied earlier [Mis hiati and Krishnaprasad, 2012℄, we just reiterate some key results in order to have a omprehensive frame- work. Allowing dierent speeds for the agents, we begin with the following sym- metry: u1ν1 = u2ν2 = u. (2.17) Then the dynami s of the relative motion ve tors, namely r = r1 − r2, g = ṙ = ν x ⊥ ⊥1 1 − ν2x2 and h = g = ṙ , an be expressed as ṙ = g ġ = uh (2.18) ḣ = −ug. 34 Now we introdu e three s alar shape variables dened as ρ = |r|, γ = (r · g)/|r| and λ = (r ·h)/|r|. Then, a ording to [Justh and Krishnaprasad, 2006℄, we have ( ) ( ) − ru = µ ⊥ r|r| · ṙ = −µ | | · h = −µλ, (2.19)r where, µ > 0 denotes the feedba k gain. As shown earlier, the dynami s of relative motion (2.18) an be redu ed to yield a se ond order dynami s given by ρ̇ = γ (2.20) ( ) γ̇ = (1/ρ− µ) δ2 − γ2 , where, δ = |g| = |h| is onserved along any traje tory of (2.18). As detailed in the original work [Mis hiati and Krishnaprasad, 2012℄, individual traje tories an be re onstru ted from the solutions of (2.20). Moreover, the solutions of the redu ed dynami s (2.20) onstitute level sets for another onserved quantity, dened as E(ρ, γ) = ρ2(δ2 − γ2)e−2µρ = E(ρ0, γ0). (2.21) However, the absen e of damping in the redu ed dynami s (2.20) has poten- tial to deteriorate the performan e of the original MMC law (2.19). A modied feedba k law, with an added dissipative term to neutralize any deviation from the predi ted traje tories, an be expressed as ( ) utot = u+ udis = −µλ+ kdλγ E(ρ, γ)−Ed , (2.22) where Ed is set as the initial value of the onserved quantity E(ρ, γ), i.e. Ed = E(ρ0, γ0). Previous work [Mis hiati and Krishnaprasad, 2010℄ has shown that this modied ontrol law (2.22) with kd > 0 makes the periodi orbit (with energy Ed) 35 orbitally asymptoti ally stable, and the orresponding domain of attra tion is hara terized by {(ρ, γ) : ρ > 0,−δ < γ < δ, (ρ, γ) 6= (1/µ, 0)}. 2.3.2 Topologi al Velo ity Alignment (TVA) Here we formalize the strategy of topologi al velo ity alignment (TVA) [Halder and Dey, 2015; Dey, 2015℄, and assume that ea h member in a group of n-agents uses this strategy to move together while keeping its heading parallel to the neigh- borhood enter of mass velo ity. Letting Ni denote the neighborhood of the i-th agent, the enter of mass (COM) velo ity of this neighborhood is given by 1 ∑ v = ν x , (2.23) COM |N | j ji j∈Ni where |Ni| represents the number of neighbors inuen ing the i-th agent. Next, by assuming that v does not vanish to zero, we dene the dire tion of the COM enter of mass motion as v x = COMN | | . (2.24)i v COM It should be noted that xN is not well-dened over a thin set in the state spa e.i As the han e of getting into the thin set is very small, we an overlook this situation for all pra ti al purposes. Now we introdu e a ontrast fun tion 1 Θi = (xN − xi) · (xN − xi) = 1− xi · xN , (2.25) 2 i i i as a quantitative measure for the misalignment between the heading of an agent and the dire tion of motion of its neighborhood enter of mass. Clearly, this ontrast fun tion (Θi) assumes its minimum value (= 0) whenever the i-th agent's 36 velo ity is aligned with its neighborhood enter of mass velo ity, and in reases monotoni ally with in rease in the misalignment between them. Thus, Θi an be interpreted as a measure of departure from our goal of a hieving alignment. Next, by assuming a non-zero velo ity for the neighborhood enter of mass (v =6 0), the TVA law is given by [Dey, 2015℄ COM ( · )xN yi iui = µ , (2.26) νi where µ > 0 denotes a positive gain, and yi arries its usual meaning. Alterna- tively, lateral a eleration for this hoi e of ontrol laws (2.26) an be expressed as ( [ ] ) lat ai = µνi xN − xN · xi xi , (2.27)i i and this provides a physi al intuition behind (2.26) as the lateral a eleration is proportional to the proje tion of the normalized velo ity of its neighborhood enter of mass onto the transverse of its own dire tion of motion. Remark 2.6. Earlier works [Justh and Krishnaprasad, 2003, 2004℄ have onsid- ered a very similar form of this ontrol law with three omponents for attra - tion (while the agents are far away), repulsion (to avoid ollision) and velo ity alignment. However, the TVA ontrol law onsiders only velo ity alignment, but extends the s ope from a planar setting to a three dimensional environment. More- over it relaxes the assumption on uniform speed of the olle tive by allowing the agents non-uniform and time-varying speed proles. This relaxation plays an im- portant role in the ontext of applying this ontrol law to a group of heterogeneous agents. 37 It was shown in [Halder and Dey, 2015; Dey, 2015℄, for a two agent system it is possible to show that the TVA strategy (2.26) aligns the velo ity of the agents, only if their velo ity ve tors were not exa tly opposite to ea h other initially. [Dey, 2015℄ also provides further analyti al results for a general n agent system with a y li intera tion s enario. As analysis of an n-agent system with neighborhood dened as the set of K-nearest neighbors poses hard hallenges, we propose an algorithmi way (Algorithm 2) to implement TVA in a real system. We bring in an additional neighbor into onsideration whenever v be omes zero. Clearly, COM this provides a way to avoid ill-posedness asso iated with v being zero be ause COM non-zero speeds of individual agents ensure that onsidering an extra neighbor will make an otherwise zero v non-zero. COM 2.3.3 Implementation Results We present the implementation results of the two ontrol laws in our roboti test- bed. In this se tion, we are presenting results for whi h the speeds of all the individual agents are same, i.e. νi = νj , ∀i, j. Though it should be noted that both ontrol laws an be implemented with dierent speeds. 2.3.3.1 Implementation of MMC Here we will show some implementation results for MMC, and demonstrate the ee tiveness of using a dissipative ontrol term (2.22). Our analysis also in ludes a omparison between the observed traje tories and traje tories obtained from 38 Algorithm 2: Topologi al Velo ity Alignment Data: Initial Time - tinitial; Final Time - tfinal; Sampling Interval - ∆; n Number of Agents - n; Initial Position and Orientation - {gi}i=1; Neighborhood Size - K begin Initialize: tcurrent ←− tinitial ; for i = 1 to n do Initialize: State - Xi ←− gi ; while tcurrent ≤ tfinal do for i = 1 to n do Dene: Ni - the set of K-nearest neighbors ; Compute: Neighborhood enter of mass velo ity - v ; COM if v = 0 then COM Dene: Ni - the set of K + 1-nearest neighbors ; Compute: Neighborhood enter of mass velo ity - v ; COM Compute: Steering ontrol - ui; n Implement: Steering Control - {ui}i=1 ; n Update: State - {Xi}i=1 ; Update: Time - tcurrent ←− tcurrent +∆ ; theoreti al predi tions, obtained via integrating the redu ed system dynami s (2.20). Considering the prersen e of a onserved quantity (2.6) in the system, we used the method des ribed in [Austin et al., 1993℄ for integration instead of general ODE solver, whi h otherwise would not be able to keep the quantity E(ρ, γ) k k onstant to its initial value. Then, from the updated values of ρ and γ , k k k we re onstru t the traje tories along with their frame ve tors, i.e. ri , xi , yi , with i = 1, 2 kand k denoting the time indi es. At ea h time instan e t , the error 39 600 600 300 400 400 200 200 )200 e 0 0 s 100 / ) ) PSfrag repla ements -200 PSfrag repla ements PSfrag repla ements-200 m m m 0 m -400 -400 mm n n n i i -600 i -600 -100 ( ( ( -800 -800 -1000 -200-1000 -1200 -1200 -300 -500 0 500 1000 -500 0 500 1000 500 1000 1500 2000 x (in mm) x (in mm) ρ (in mm) (a) Traje tory (Experiment) (b) Traje tory (Simulation) ( ) Phase plot Figure 2.6: Performan e of the system signi antly improved upon addition of −1 the dissipative ontrol term (with µ = 0.001 mm , ν1 = ν2 = 200 mm/se and kd = 1× 10−15 −6 3mm se ) (ek k ki ) is omputed as: ei = |ri,expt − rki,ideal|. The plots of a sample run using the modied MMC feedba k law (2.22) are shown in Fig 2.6. This modied ontrol law has been applied with kd = 1 × 10−15 −6 3 −1mm se , and the parameters µ, ν1 and ν2 are sele ted as 0.001 mm , 200 mm/se and 200 mm/se , respe tively. The resulting performan e is quite satisfa tory as shown in Fig 2.6a (refer [YouTube℄ for implementation video). We have also observed that the error is bounded (∼ 250 mm) within the size of the robots (∼ 400 mm). 2.3.3.2 Implementation of TVA We implemented the TVA ontrol law (2.26) in a 2 dimensional setting (i.e. vi(t) is ignored). As the implementation is in dis rete time, we followed Algorithm 2 in our implementation in order to avoid the singular ase of |v | = 0. To demonstrate COM the performan e of TVA ontrol law, we designed three dierent experiments (refer 40 y y γ 1000 2000 1000 500 0 1000 0 −500 −1000 0 −1000 −2000 −1000 −1500 −3000 −2000 −2000 −4000 −2500 −3000 −2000 −1000 0 1000 2000 −2000 −1000 0 1000 2000 −5000 −4000 −3000 −2000 −1000 0 1000 2000 x (in mm) x (in mm) x (in mm) (a) Traje tories: 8 agents, o k- (b) Traje tories: 8 agents, split- ( ) Traje tories: 6 agents, pertur- ing ting bation 4.5 1.4 4 5 4 1.2 3.5 4 1 3 2.5 0.8 3 2 2 0.6 2 1.5 0.4 1 1 PSfrag repla ements PSfrag repla ements 0.2 PSfrag repla ements 0.5 0 0 0 0 0 10 20 30 40 50 0 5 10 15 20 0 20 40 60 80 Time (in sec) Time (in sec) Time (in sec) (d) Contrast Fun tion: 8 agents, (e) Contrast Fun tion: 8 agents, (f) Contrast Fun tion: 8 agents, o king splitting perturbation Figure 2.7: Robot traje tories and ontrast fun tions of TVA for (i) Experiment 1 [Fig (a),(d)℄ with 8 agents, demonstrates o king behavior (K = 3); (ii) Ex- periment 2 [Fig (b),(e)℄ with 8 agents, des ribes the splitting behavior due to low neighborhood size (K = 1), and (iii) Experiment 3 [Fig ( ),(f)℄ with 6 agents, shows that perturbation an ause a swarm to split, the traje tory of the per- turbing agent is not shown. (µ = 1 Hz and νi = 60 mm/se is kept xed for all experiments.) [YouTube℄ for implementation videos). In these experiments, the sonar sensors on the robots were a tivated to sense any obsta le in the dire tion of motion of the robots and if any robot an sense su h an obsta le, it will simply apply a sat maximum turning rate (ω ) to avoid ollision. The sonars are programmed to 41 y (in mm) Θ(t) y (in mm) Θ(t) y (in mm) Θ(t) Number of nearest neighbors dete t an obsta le only in lose proximity (∼ 300 mm) of the robots. In all our sat experiments ω is taken to be 50 rad/se , forward speeds of all of the robots are kept onstant at 60 mm/se and the value of the parameter µ is hosen to be 1 Hz. A system with eight agents is onsidered and we apply same TVA law to all of them. The neighborhood size is taken to be three (i.e. K = 3). The robots are initially pla ed in arbitrary positions and dire tions. The footprints of ∑ the robots and the orresponding ontrast fun tion, Θ(t) = i Θi(t) is plotted against time in Fig 2.7a, 2.7d. The initial and nal dire tions of the robots are shown using arrows and the nal positions of the robots are denoted using dots. It an be seen from Fig 2.7d that the ontrast fun tion de ays to zero very qui kly whi h indi ates perfe t velo ity alignment within the swarm. Next we de reased the neighborhood size and made it one (K = 1), so that every robot only ` ommuni ates' with its losest neighbor. We hose the initial positions in su h a way that they may form sub- lusters instead of moving as a single swarm. This behavior is alled `splitting ' in a swarm. From Fig 2.7b, we an learly see that the swarm of eight robots gradually split from ea h other and form three dierent lusters. It is to be noted that even if all the agents are not going in the same dire tion, the ontrast fun tion still onverges to zero (Fig 2.7e). This happens be ause ea h of the robots are aligned with their nearest neighbors and hen e ea h of the individual ontrast fun tions (Θi(t)) are zero. This experiment may explain the splitting phenomenon observable in nature. Lastly, we ombined the above two experiments, and ondu ted an experiment using six robots in a 42 swarm and another robot as a predator. A separate omputer was used for manual ontrol of the `predator' robot. At the beginning, neighborhood size is kept at K = 3, su h that the ` om- muni ation' graph among the robots stays onne ted and they move as an entire swarm in a ommon dire tion. When the swarm omes lose to the predator, the neighborhood size is de reased to one. As we are not using any onboard visual sensing and the sonar sensing is done only in very lose region (∼ 300 mm), the hange in neighborhood size is made manually. From Fig 2.7f, we an see that the hange in neighborhood size takes pla e at around 20 se onds and we an also see a tiny jump in the ontrast fun tion at that time. The predator then slowly approa hes to one of the agents in the swarm, whi h abiding to its ollision avoidan e rule, turns to avoid the predator. In Fig 2.7 , the traje tories of the agents are drawn in dashed lines before the o urren e of this event and in solid lines afterwards. The traje tory of the predator robot in not shown in the g- ure. After reating the initial perturbation, the predator is slowly moved through the swarm ausing some subsequent disturban es. These perturbations reate a noti eable impa t in the swarm. As the atta ked agent turns, its neighbor also tries to align itself with that agent and so does its neighbor. This goes on until the ommuni ation graph be omes dis onne ted and a split in the swarm is then observed [YouTube℄ just like in Experiment 2. As we an see in Fig 2.7 , the swarm is divided in two groups after the atta k of the predator. The jumps in plot of the ontrast fun tion in Fig 2.7f symbolize the perturbations aused by 43 the external agent. The ontrast fun tion eventually onverges to zero after the members are aligned with their neighbors within ea h subgroup. 44 Chapter 3 Optimal Steering of Agents on a Plane 3.1 Introdu tion The kinemati uni y le model is often used in path-planning for ground vehi les, sin e the onguration of a ground vehi le an often be represented by a point in a plane that is onstrained to move in the dire tion of the urrent heading [Bellaï he et al., 1998; LaValle, 2006℄. The state of this system an be represented as an element of the spe ial Eu lidean group SE(2), where the ontrol inputs are a urvature input whi h ontrols the rate of hange of the heading angle, and a velo ity input whi h ontrols the rate of hange of the uni y le position in the dire tion of the heading angle. Given the urrent onguration of the uni y le and a desired future ong- uration, an admissible path for moving the uni y le from an initial to a nal onguration an be determined via the minimization of some ost fun tional. There an be many variations to this problem depending on the hosen ost fun - 45 tional. A mu h elebrated problem is the problem of Euler's elasti a [Euler, 1744℄ where the minimum urvature path joining two given ongurations on the plane is onsidered. Parti ularly elegant and other well-known variations in lude the minimum-time solutions of [Dubins, 1957℄ and [Reeds and Shepp, 1990℄. Optimal paths of the minimum-time problem onsist only of straight-line and ir ular-ar segments whi h, when pat hed together, reate dis ontinuities in the path ur- vature and ause potential di ulty in implementation sin e abrupt hanges in urvature are hard to tra k. Proposed modi ations that alleviate this problem enfor e that the urvature stay ontinuous, e.g., [Frai hard and S heuer, 2004℄, yet it is also possible to penalize the total urvature along the path in the expe - tation that the optimal urvature will be ontinuous. [Halder and Kalabi , 2017℄ takes the latter approa h, onsidering the minimization of the urvature along a path onne ting initial and nal uni y le ongurations with free nal time. In this hapter, we will present a problem that penalizes both the urvature and speed ontrols in maneuvering a uni y le from initial to desired nal ong- uration. This helps both the urvature and speed ontrols to be smooth along an optimal traje tory. These optimal traje tories losely resemble to those obtained in [Halder and Kalabi , 2017℄. Our solution is obtained using geometri optimal ontrol, where the ne essary onditions for optimality are obtained via the Pon- tryagin Maximum Prin iple (PMP), and Lie-Poisson redu tion [Krishnaprasad, 1993; Ohsawa, 2013℄. Using geometri optimal ontrol on SE(2) to nd solutions to path-planning problems has also been onsidered by [Sussmann and Tang, 1991; 46 Krishnaprasad, 1993; Agra hev and Sa hkov, 2004; Dey and Krishnaprasad, 2014; Justh and Krishnaprasad, 2015b℄. Studying a uni y le on the plane is important sin e often a olle tive of N agents exhibits single agent behavior under syn hronization. In the later part of this hapter, we will present a framework for analyzing su h a olle tive. This framework is based on [Justh and Krishnaprasad, 2015b℄, where every agent was assumed to have onstant speed. Here we will, however, onsider a ost fun tion that penalizes both speed and urvature ontrols. In addition to the individual ontrol osts, there is one ost that is attributed to the `mismat h in ontrol' of an agent with its `neighbors'. The neighbors of any agent is di tated by a xed graph of intera tion. The strength of su h intera tion is aptured by a oupling parameter. It is shown that in extreme ases (no oupling and high oupling) the optimal olle tive ontrols are dire tly asso iated with optimal ontrols for a single agent problem. This framework is what we use in the later hapters of this thesis where we onsider optimal ontrol problems for a ontinuum of agents. 3.2 Optimal Steering of a Uni y le In this se tion, we onsider minimizing the urvature and speed ontrol osts of a path in SE(2) onne ting an initial uni y le onguration g0 with its desired nal onguration gT at time T . We formulate this optimization as a geometri optimal ontrol problem and derive the ne essary onditions using PMP and Lie- Poisson redu tion. From the ne essary onditions, we show that there are two 47 onstants of motion: the Hamiltonian and the Casimir. We show that there are three possible families of solutions depending on the values of Casimir and the Hamiltonian. In the rst ase, the motion onsists of segments of a U-turn; in the se ond ase, the motion onsists of segments of parallel parking traje tories; in the third ase, the motion onsists of straight lines or asymptoti approa hes thereto. Consider the uni y le kinemati equations, ẋ(t) = v(t) cos θ(t), (3.1a) ẏ(t) = v(t) sin θ(t), (3.1b) θ̇(t) = u(t), (3.1 ) where (x(t), y(t)) ∈ 2R is the position of the uni y le on the Cartesian plane, θ(t) ∈ S1 is the heading of the uni y le, v(t) is the uni y le speed ontrol, and u(t) is the steering ontrol, equal to the rate of hange of the heading θ(t). The onguration of the uni y le an be represented as an element of the matrix Lie group SE(2). Let g(t) ∈ SE(2) where,   cos θ(t) − sin θ(t) x(t)       g(t) = sin θ(t) cos θ(t) y(t) . (3.2)       0 0 1 Then the equations (3.1) an be written in left-invariant form, ġ(t) = g(t)ξ(u(t), v(t)), (3.3) where, ξ(u(t), v(t)) = u(t)X1 + v(t)X2, (3.4) 48 and,     0 −1 0 0 0 1             X 1 = 1 0 0 , X2 = 0 0 0 . (3.5)             0 0 0 0 0 0 The matri es X1 and X2 are elements of the Lie algebra se(2). Together with X3 = [X1, X2] = X1X2 −X2X1, X1 and X2 form a basis for se(2). Without loss of generality, we an assume that g0 = I3, sin e g(t) an always −1 be redened a ording to g(t) := g0 g(t). Given a nal time T > 0 and a nal state gT , ost fun tion to be minimized is, ∫ 1 T ( ) u(t)2 + v(t)2 dt. (3.6) 2 0 3.2.1 Optimal Control Solution In order to solve the problem we form the pre-Hamiltonian, H = 〈 1p, gξ(u, v)〉 − (u2 + v2), (3.7) 2 where p(t) ∈ SE(2)∗ is the adjoint variable. To simplify the Hamiltonian, we perform Lie-Poisson redu tion, introdu ing the variable µ(t) ∈ se(2)∗ satisfying the translation to identity, 〈µ, ξ(u, v)〉 = 〈p, gξ(u, v)〉. As an element of the dual spa e, µ(t) an be represented as µ(t) = µ X♭1 1+µ ♭ 2X2+ µ X♭ {X♭ ♭3 3, where 1, X2, X♭3} are the basis ve tors dual to {X1, X2, X3}. 49 The pre-Hamiltonian therefore be omes, 1 H(µ, u, v) = 〈µ, ξ(u, v)〉 − (u2 + v2), 2 = 〈µ1X♭1 + µ2X♭2 + µ3X♭3, uX1 + vX2〉 − 1(u2 + v2), 2 1 = µ1u+ µ2v − (u2 + v2). 2 ∗ ∗ A ording to the PMP, the optimal ontrol (u , v ) satises, H(µ∗, u∗, v∗) = max H(µ∗, u, v). (3.8) (u,v)∈R2 Therefore the optimal ontrols are given by, u∗1 = µ1, (3.9) u∗2 = µ2, (3.10) The redu ed Hamiltonian is therefore, 1 ( ) H = H(µ, u∗1, u ∗ 2) = µ 2 2 2 1 + µ2 , (3.11) and the dynami s of the µi variables are given by [Krishnaprasad, 1993℄,     µ̇1 0 −µ   3 µ2         ∂h µ̇  =  µ 0 0  , (3.12)  2  3  ∂(µ1, µ2, µ3)         µ̇3 −µ2 0 0 µ̇1 = −µ2µ3, (3.13a) µ̇2 = µ1µ3, (3.13b) µ̇3 = −µ1µ2. (3.13 ) 50 (a) c = 0.5 (b) c = 1.5 ( ) c = 1 Figure 3.1: M (bla k) plotted as an interse tion of C (red) and H (blue) for h = 1 and three values of the Casimir c Let, c = µ22 + µ 2 3. (3.14) This variable is alled the Casimir and it is a onstant of motion, implying that √ the variables µ2 and µ3 evolve on ir le of radius c. Along with the Casimir, the Hamiltonian (3.11) is also a onserved quantity of motion. For onvenien e of subsequent al ulations, we will work with the following s aled Hamiltonian, h = 2H = µ21 + µ 2 2. (3.15) 3.2.2 Chara terizing the Types of Motion A ording to (3.14) and (3.15), the dynami s (3.13) evolve on the manifold M where, M = C ∩ H, C = {(µ1, µ2, µ3) ∈ 3R : c = µ22 + µ23}, (3.16) H = {(µ1, µ2, µ3) ∈ 3R : h = µ2 21 + µ2}. 51 The manifold M is one-dimensional and equal to the interse tion of the ylinders C and H. The shape of M is determined by the Casimir c, while the shape of H is determined by h. The motion of µ evolves on M. Due to ontinuity, it must evolve on a onne ted omponent of M, so it is important to onsider the types of possible interse tions, of whi h there are three orresponding to three dierent ases: in Case 1, C is stri tly smaller than H; in Case 2, C is stri tly larger than H; in Case 3, C is equal in size to H. To perform a ase-by- ase ategorization of M √, we note that, a ording to (3.14), the variable µ2 is restri ted to ± c and, √ √ √ a ording to (3.15), µ2 is restri ted to ± h. In Case 1, c < h, so |µ2| ≤ c < h. Therefore the motion evolves on a onne ted omponent of M where µ1 does not 2 hange sign sin e µ1 = h − µ22 ≥ h − c > 0. Similarly, in Case 2, c > h, so √ |µ2| ≤ √ h < c. The motion evolves on a onne ted omponent of M where µ3 2 2 does not hange sign sin e µ3 = c − µ2 ≥ c − h > 0. In Case 2, c = h, so the two onstraints agree at the extremes. Instead of having M as a one onne ted omponent, it a tually has four dis onne ted omponents. These omponents √ meet ea h other asymptoti ally at the extremes µ1 = µ3 = 0, µ2 = ± h. See Fig. 3.1 for a visualization of the three ases. In the following, we study the three types of motion in further detail. Case 1: 0 < c < h The equations (3.13) admit the following expli it solutions by means of Ja obian 52 6 6 6 5 4 4 4 2 2 3 0 0 2 -2 -2 1 -4 -4 0 -6 -6 -0.5 0 0.5 1 1.5 2 -2 0 2 4 -1 0 1 2 (a) c = 0.5 (solid bla k) and 0.75 (b) c = 1.5 (solid bla k) and 2 ( ) c = 1 (dot-dashed red) (dot-dashed red) Figure 3.2: Extremal traje tories for dierent ases (solid bla k)  in reasing c produ es the dot-dashed red urve; in reasing η produ es the solid blue urve; hanging s1 or s2 produ es the dashed green urve. h = 1 for all three ases. ellipti fun tions [Davis, 1962; Byrd and Friedman, 1971℄, √ (√ ) µ1 = s1 h dn h(t+ η), k , (3.17a) √ (√ ) µ2 = s1 c sn h(t+ η), k , (3.17b) √ (√ ) µ3 = c n h(t+ η), k , (3.17 ) 2 c where the modulus k of the ellipti fun tions is given by k = . The parameters h s1 ∈ {1,−1} and η ∈ R do not depend on c and h, but on the endpoint onstraints. 4K√(k) Note that these solutions are periodi with the period given by T1 := , whereh K(k) denotes the omplete ellipti integral of the rst kind. Proposition 3.2.1. Let 0 < c < h. Assume η = 0, then any extremal traje tory 53 is given by, s ( (√ ))1 x(t) = 1− dn ht, k (3.18a) k 1 (√ ( (√ ) )) y(t) = ht− E am ht, k , k (3.18b) k (√ ) ( (√ )) ( (√ )) θ(t) = s1 · am ht, k = s1 cos−1 n ht, k = s1 sin−1 sn ht, k , (3.18 ) where E(·, ·) denotes the in omplete ellipti integral of the se ond kind. Proof. We have from (3.9) and (3.17a), √ (√ ) u1 = µ1 = s1 h dn ht, k . Integrating the equation θ̇ = u1 gives (√ ) ( (√ )) ( (√ )) θ(t) = s1 · am ht, k = s cos−1 ht, k = s sin−11 n 1 sn ht, k . Sin e the optimal speed ontrol is given by (3.10) and (3.17b), √ (√ ) u2 = µ2 = s1 c sn ht, k , we may now integrate x and y dynami s to obtain position variables as fun tions of time. We have, ∫ t x(t) = u2 cos(θ(t))dt 0 ∫ t √ (√ ) (√ ) = s1 c sn ht, k n ht, k dt 0 √ − c 1 [ (√ )]t = s1√ · dn ht, k h k2 t=0 s ( (√ ))1 = 1− dn ht, k (3.19) k 54 Similarly, ∫ t y(t) = u2 sin(θ(t))dt 0 ∫ t√ (√ ) = c 2sn ht, k dt √0 c 1 [√ ( (√ ) )]t = √ · ht− E am ht, k , k h k2 t=0 1 (√ ( (√ ) )) = ht−E am ht, k , k (3.20) k  Studying the extremal traje tories in η = 0 ase is important sin e any other traje tory an be expressed by means of these traje tories after a suitable trans- lation and rotation. This is demonstrated in the following proposition. Proposition 3.2.2. Assume 0 ≤ c < h. Let us denote an extremal traje tory belonging to η = 0 ase by (x0(t), y0(t), θ0(t)) a ording to (3.18). Then any extremal traje tory generated by (3.17) an be expressed as,       x(t) x (t+ η) x (η)    0   0    = R(−ψ) −   , (3.21)       y(t) y0(t+ η) y0(η) θ(t) = θ0(t + η)− θ0(η), (3.22) where R(·) is the planar rotation matrix and ψ is dened as, (√ ) ψ := s · am h η = θ0(η). (3.23) 55 Proof. Indeed, equation (3.22) is readily a hieved by integrating (3.17a). Further, noti e that √ (√ ) ( (√ ) ) ẋ(t) = s1 c sn h(t+ η), k cos s1 · am h(t+ η), k − ψ √ ( (√ ) (√ ) = c s1 · cos(ψ) · sn h(t+ η), k n h(t+ η), k (√ )) + sin(ψ) · 2sn h(t + η), k , and, √ (√ ) ( (√ ) ) ẏ(t) = s1 c sn h(t + η), k sin s1 · am h(t+ η), k − ψ √ ( (√ ) (√ ) = c −s1 · sin(ψ) · sn h(t+ η), k n h(t + η), k (√ )) +cos(ψ) · 2sn h(t+ η), k . Compa tly, we have       ẋ(t) cos(ψ) sin(ψ) ẋ0(t + η)        =   ·   . (3.24)       ẏ(t) − sin(ψ) cos(ψ) ẏ0(t + η) Integration of (3.24) yields (3.21).  Case 2: c > h This ase admits solutions analogous to those of Case 1 and we will pro eed in a similar way. To start, we write down solutions to (3.13) as, √ (√ ) µ1 = h n c(t+ η), k , (3.25a) √ (√ ) µ2 = s2 h sn c(t + η), k , (3.25b) √ (√ ) µ3 = s2 c dn c(t+ η), k , (3.25 ) 56 2 h with the modulus k = . s2 ∈ {1,−1} and η ∈ R are similar parameters as inc ase 1. These solutions are periodi as well, with period given by T := 4K√(k)2 .c Proposition 3.2.3. Let h > 0 and c > h. Assume η = 0, then any extremal traje tory is given by, ( (√ )) x(t) = s2 · k 1− n c t, k (3.26a) (√ ( (√ ) )) y(t) = s2 c t− E am c t, k , k (3.26b) θ(t) = cos− ( (√ )) ( (√ )) 1 dn c t, k = sin−1 k sn c t, k , (3.26 ) where E(·, ·) denotes the in omplete ellipti integral of the se ond kind. Proof. The proof is similar to that of Proposition 3.2.1 and uses elementary inte- grals of Ja obi ellipti fun tions.  Analogous to Proposition 3.2.2, we have the following result. Proposition 3.2.4. Assume h > 0 and c > h. Let us denote an extremal traje - tory belonging to η = 0 ase by (x0(t), y0(t), θ0(t)) a ording to (3.26). Then any extremal traje tory generated by (3.17) an be expressed as,       x(t) x0(t+ η) x0(η)        = R(−ψ) −   , (3.27)       y(t) y0(t+ η) y0(η) θ(t) = θ0(t + η)− θ0(η), (3.28) where R(·) is the planar rotation matrix and ψ is dened as, − ( (√ )) ( (√ ))ψ := cos 1 c η, k = sin−1dn k sn c η, k = θ0(η). (3.29) 57 Proof. The proof is essentially the same to that of Proposition 3.2.2. Note the dierent denition of ψ in (3.29)  Case 3: c = h 2 This ase is transitional between ase 1 and ase 2. Putting the modulus k = c = 1 in either of the equations (3.17) or (3.25), we get the following solutions h √ (√ ) µ1 = s1 c se h c(t + η) , (3.30a) √ (√ ) µ2 = s2 c tanh c(t + η) , (3.30b) √ (√ ) µ3 = s3 c se h c(t + η) , (3.30 ) with s1, s2, s3 ∈ {1,−1} and η ∈ R. We readily obtain the following result. Proposition 3.2.5. Let c = h > 0. Assume η = 0, then any extremal traje tory is given by, ( (√ )) x(t) = s2 1− se h c t , (3.31a) (√ (√ )) y(t) = s1s2 c t− tanh c t , (3.31b) ( (√ )) θ(t) = s1 tan −1 sinh c t . (3.31 ) Proof. Indeed, integrating θ̇ = u1 = µ1 yields (3.31 ). Sin e u1 = µ2, we get √ (√ ) (√ ) ẋ = s2 c se h c t tanh c t , √ (√ ) ẏ = s1s2 c tanh 2 c t , whi h in turn gives (3.31a)(3.31b) after integration.  58 In the same spirit as before, we write any general extremal traje tories in terms of these traje tories. Proposition 3.2.6. Assume c = h > 0. Let us denote an extremal traje tory belonging to η = 0 ase by (x0(t), y0(t), θ0(t)) a ording to (3.31). Then any extremal traje tory generated by (3.30) an be expressed as,       x(t) x (t+ η) x (η)    0   0    = R(−ψ) −   , (3.32)       y(t) y0(t+ η) y0(η) θ(t) = θ0(t + η)− θ0(η), (3.33) where R(·) is the planar rotation matrix and ψ is dened as, ψ := s tan− ( (√ )) 1 1 sinh c η = θ0(η). (3.34) This ase onsists of two types of solutions: a straight line solution, orrespond- ing to the sub ase where the nal ondition lies on the x-axis, and an asymptoti solution, whi h asymptoti ally approa hes a straight line with slope cotψ. See Fig. 3.2 for a graphi . 3.2.3 On Time-optimality From an engineering perspe tive, it seems appealing to onsider a similar problem where nal time T is free. Therefore we want to rea h from initial onguration g 2 20 to the nal onguration gT in a minimal time so that the ontrol ost (u +v ) is minimized along the optimal traje tory. It also makes sense to add a penalty 59 to the time it takes for the uni y le to omplete the maneuver. The ost fun tion an be expressed as ∫ 1 T min (a+ u(t)2 + v(t)2)dt, (3.35) T,u(·),v(·) 2 0 for some time-penalty parameter a > 0. Note that without the penalty on time, a solution ould orrespond to a prohibitively large nal time whi h is not desired from a pra ti al viewpoint. This time-optimal version is just a spe ial ase of what we have onsidered in se tion 3.2. To see this, we ompute the Hamiltonian 1 2 2 (3.11) as, H = (µ1 + µ2 − a). Sin e time-optimality requires the Hamiltonian to2 2 2 be identi ally zero, we have spe ial ase of h = µ1 + µ2 = a ( .f. (3.15)). This √ also gives the bounds of the optimal ontrols u1 = µ1, u2 = µ2 to be within ± a. We an, therefore, use the parameter a to set a desired bound on the ontrols that is permitted by physi al onstraints. A losely related problem was studied in [Halder and Kalabi , 2017℄, where the speed ontrol v was assumed to be of onstant magnitude, and the minimum time problem asso iated with minimum urvature path was onsidered. 3.3 Optimal Control of a Colle tive of Agents Now we onsider a olle tive of N agents moving on the plane. Motion of ea h agent an be modeled by the uni y le dynami s (3.1). As seen before, this dynam- i s an be equivalently expressed as a ontrolled dynami s in SE(2), ġk = gkξk(uk), 60 where ξk(uk) = uk1X1 + uk2X2, k = 1, 2, ..., N. (3.36) We suppose these agents intera t among themselves dire ted by a xed adja en y matrix A = [a ] ∈ RN×Nij . aij = 1 if agent i and agent j intera t and aij = 0 otherwise. The nature of the intera tion will be made pre ise shortly. LetD be the degree matrix, i.e. the diagonal matrix where the i-th diagonal entry represents the number of agents the i-th agent intera t with. Then the graph Lapla ian is dened by B := D − A. In this setup, we seek to minimize ∫ T L = L(ξ1(u1(t)), ..., ξN(uN(t)))dt 0 ( ) ∫ 1 T N N N ∑ ∑∑ = |ξ 2 2k| + χ akj |ξk − ξj| dt, (3.37) 2 0 k=1 k=1 j=1 for some onstant χ ≥ 0 and the xed endpoint onditions gk(0) = gk0, gk(T ) = √ gkT , k = 1, ..., N T. Note that we used the tra e norm |ξ| = tr(ξ ξ). The param- eter χ is alled a oupling onstant sin e it a ts as a weight to the se ond term in the ost fun tional (3.37). Without the oupling term, this problem simplies to solving N opies of the single agent problem as onsidered in Se . 3.2. The oupling term penalizes agent k through the `mismat h in ontrol' with the agents that it is intera ting with (i.e. nonzero entries of k-th row of the matrix A). This type of ost fun tional is aimed to apture the `allelomimeti behavior' or the tenden y to opy neighbors in a natural olle tive. [Justh and Krishnaprasad, 2015b℄ studies a very similar problem where the speeds of ea h agent is assumed to be onstant. Here the speed ontrols (uk2, k = 1, 2, ..., N) are to be determined 61 by solving the optimal ontrol problem (3.37). Sin e the underlying optimal ontrol problem is essentially the same, we will use the results from [Justh and Krishnaprasad, 2015b℄ to derive rst order opti- mality onditions using the Pontryagin's Maximum Prin iple and the Lie-Poisson ∗ redu tion te hnique. Denote µk ∈ se (2) by, 3 ∑ µ = µ X♭k ki i , (3.38) i=1 [ ]T ♭ ♭ ♭ where {X1, X2, X3} is dual basis to {X1, X2, X3}. If we dene µ̃k = µk1 µk2 µ ,k3 rst order optimality (PMP) yields     u µ̃  1  1         .   = Ψ . . . , (3.39) . .             uN µ̃N where Ψ = ((IN + 2χB)⊗ −1I2) = (IN + 2χB)−1 ⊗ I2. (3.40) ∗ N Here ⊗ denotes the Krone ker produ t. The redu ed Hamiltonian in (se (2)) takes the form   µ̃  1 [ ]   1   h = · · · Ψ .µ̃T µ̃T  . . (3.41)2 1 N .      µ̃N The Lie-Poisson redu ed dynami s is then expressed as follows. Dene µ = 62 [ ]T µT · · · µT . Then,1 N µ̇ = Λ(µ)∇h, (3.42) where Λ(µ) = −diag(Ω1, · · · ,ΩN), with   0 µk3 −µ k2     Ωk = −µ  , k = 1, ..., N. (3.43)  k3 0 0      µk2 0 0 [ ]T [ ]T Also we have ∇h = (∇h) · · · (∇h) , (∇h)k = ∂h ∂h1 N 0 , k =∂µk1 ∂µk2 1, ..., N . Note that along with the Hamiltonian (3.41), there are N Casimirs whi h are onstants of motion. The Casimirs are dened as, c 2 2k = µk2 + µk3. (3.44) In both the extreme ases (i) no oupling (χ = 0) and (ii) high oupling (χ → ∞), the olle tive optimal problem simplies to studying the single agent problem as onsidered in Se tion 3.2. The χ = 0 ase is immediate. The details of the high oupling limit χ → ∞ is worked out in [Justh and Krishnaprasad, 2015b℄ (see se tion 3( )). Dening the quantities N N N 1 ∑ 1 ∑ 1 ∑ α1 = µj1, α2 = µj2, α3 = µj3, (3.45) N N N j=1 j=1 j=1 we obtain the following dierential equations α̇1 = −α2α3, α̇ = α (3.46)2 1α3, α̇3 = −α1α2. 63 These are the same equations we already obtained and analyzed in details for a single agent ase (3.13). 3.4 Con luding Remarks In this hapter, we studied an optimal ontrol problem of a uni y le on the plane in detail. First order ne essary onditions are obtained by using the Pontryagin's Maximum Prin iple and the Lie-Poisson redu tion te hnique. All possible mo- tion types are properly ategorized by the relative values of the Casimir and the Hamiltonian. In the later part, we presented a framework for studying a lass of optimal ontrol problems involving many agents. These agents intera t with ea h other by a pre-determined intera tion graph. The intera tion enters into the optimal ontrols of the agents through the additive ` ontrol-mismat h' term in ost fun tional. This type of ost has been used in literature [Justh and Krish- naprasad, 2015a,b℄ to apture the `allelomimeti behavior' in natural o ks. The single agent ase, onsidered in this hapter, emerges naturally in a syn hroniza- tion limit of the olle tive model. This olle tive framework gives us the starting point to on eptualize a ontinuum o k where we study the limiting ase of N → ∞ under a spe i intera tion graph. These topi s are des ribed in detail in later hapters of this thesis. 64 Part II Analysis of Colle tive Motion: Top-down Approa hes 65 Chapter 4 Continuum Flo king and Control 4.1 Motivation It is a ommon pra ti e in lassi al me hani s to onsider a ontinuous des ription of a physi al system. Appli ations in lude vibrating rods, vibrating membranes, uid me hani s et . [Goldstein et al., 2001; Chorin et al., 1990℄ The transition from dis rete-parti le system to a ontinuum enables a ompa t des ription of the system, often leading to partial dierential equations that reveal deep insights into the system, e.g. wave-like phenomena, whi h may be too obs ure or inelegant in the dis rete ounterpart. In the same spirit, we attempt to on eptualize a ontinuum o k and address its optimal maneuvering properties. Biologi al o ks are known to show remarkable response to predator atta ks. In the ase of an atta k, the whole o k seems to divert away from the predator. They an perform these tasks my means of propagating information (in this ase, threat) through the o k at a mu h higher speed than the o king speed [Attanasi et al., 2014℄. 66 A key goal of our approa h is to apture this phenomenon, i.e. to un over the wave-like aspe ts of o king. An optimal ontrol problem for a o k of nite agents is presented in Se tion 3.3. We take the same framework and formulate a general ontinuum version of the problem. Consider an ensemble of N identi al, self-steering parti les, ea h obeying a drift free left invariant system on a matrix Lie group G, ġi = giξi, gi ∈ G, ξi ∈ g, i = 1, ..., N, (4.1) where g is the Lie algebra of G. These parti les intera t with ea h other by a pre-dened graph of intera tion. The olle tive behavior of su h a system was extensively studied in [Justh and Krishnaprasad, 2015b℄ by solving an appropriate optimal ontrol problem. Borrowing the notations of [Justh and Krishnaprasad, 2015b℄, we write the ost fun tional as self-energy term oupled with a mismat h in steering term ∫ T J = L(ξ1(t), ..., ξN(t))dt, (4.2) 0 where, ( ) N N N 1 ∑ ∑∑ L(ξ1, ..., ξ ) = ‖ξ ‖2N i + χ aij ‖ξi − ξ 2j‖ , (4.3) 2 i=1 i=1 j=1 where binary valued aij 's populate the adja en y matrix that denes the graph of intera tion and χ ≥ 0 is a oupling onstant. Note that the inner produ t √ 〈ξ, η〉 = (ξTtr η), and the orresponding tra e norm, ‖ξ‖ = 〈ξ, ξ〉 are in ee t, where ξ, η ∈ g. 67 In an attempt to extend this view, we onsider an innite number of parti les, i.e. the limiting ase of N → ∞. Here we onsider a one-dimensional ontinuum 1 of parti les, i.e. ea h parti le is labeled by a point on a ir le S . This way, the agents are thought as a virtual lament. Moreover, we onsider a y li intera tion graph, i.e. ea h parti le is thought to be intera ting with the `next' parti le on the ir le. We introdu e the maps, g : R × S1 → G and ξ : R × S1 → g. The mismat h in steering term an then be written as the gradient of ξ in the limiting ase. In other words, in ontinuum limit of N → ∞, the Lagrangian in (4.3) take the form, ( ) ∥ ∥ ∫ 2 1 2π ∥∂ξ(t, θ)∥ L(ξ) = ‖ξ(t, θ)‖2 + χ ∥ ∥ dθ. (4.4) 2 ∥ ∂θ ∥0 Note that the summations over the number of parti les in (4.3) have been repla ed by integral over the ir le in the ontinuum setting in (4.4). Let n be the dimension of the Lie algebra g and {A1, A2, ..., An} denote an or- 1 thonormal basis of g. We introdu e the ontrols ui : R × S → R, i = 1, ..., m, so that ξ an be written as, m ∑ ξ(t, θ) = ui(t, θ)Ai, (4.5) i=1 where m < n. With this substitution, the Lagrangian in (4.4) an be rewritten as, ( ) m ( )∫ 2 1 2π ∑ ∂ui L(u1, ..., um) = u 2 + χ dθ. (4.6) 2 i0 ∂θi=1 68 Finally, we attempt to minimize the ost fun tional, ∫ T J = L(u1, ..., um)dt, (4.7) 0 ( ) m ∂g(t,θ) ∑ subje t to the group dynami s, = g(t, θ)ξ(t, θ) = g(t, θ) ui(t, θ)Ai and∂t i=1 the xed end-point onstraints, g(0, θ) = g0(θ) and g(T, θ) = gT (θ). We will note that this optimal ontrol problem an be ast in a more onvenient setting of loop groups, the group of smooth fun tions from the ir le to the Lie group G. In Se tion 4.2, we will develop a general framework for su h optimal ontrol problems in loop group setting. Controllability results will be dis ussed in Se tion 4.3. This helps us to des ribe optimal ontrol solutions in Se tion 4.4. Ne essary onditions will be derived by both al ulus of variations and Pontrya- gin's Maximum Prin iple approa h. An example of ontinuum of nonholonomi integrators will be studied in detail in Se tion 4.5. Se tion 4.6 will present deriva- tion of optimal ontrol equations in the SE(2) ase along with their numeri al treatment. This is a joint work with Dr. E. Justh [Halder et al., 2019a℄. 4.2 A Control System on a Loop Group Let G be a nite dimensional matrix Lie group and g be its Lie algebra of di- 1 mension n. We will study spa es of smooth maps from the ir le S to G and g, ∞ 1 G = C (S ;G), L = C∞(S1; g). 69 We an onstru t Sobolev ompletions of G and L as done in [Krishnaprasad et al., 1983℄. We an always view the Lie algebra g as a subalgebra of the general linear algebra gl(r,R) ∞ 1for some r > n. Dening the spa e R = C (S ; gl(r,R)), we have G ⊂ R, L ⊂ R. For any f ∈ R, use the Sobolev k-norm (k ≥ 1) k ∣ ∣∫ l 2 ∑ ∣ ∣ ‖ ‖ df ∣ ∣k = f(θ) dθ, (4.8) ∣dθl ∣ 1 l=0S where |f |2 = tr(fTf). Let the ompletions of G , L , and R in this norm be denoted as Gk, Lk, and Rk, respe tively. By Proposition 3.1 of [Krishnaprasad et al., 1983℄, Gk is a tually a Lie group under pointwise multipli ation operation (g1g2)(θ) = g1(θ) · g2(θ), g1, g2 ∈ G 1k, θ ∈ S . Moreover, Lk is the Lie algebra of Gk under pointwise Lie bra ket dened as [f1, f2](θ) = [f1(θ), f2(θ)], f1, f2 ∈ Lk, θ ∈ S1. The spa es Gk and Lk are alled loop groups and loop algebras [Pressley and Segal, 1986℄. Similar to the nite dimensional Lie groups, we introdu e (pointwise) left a tion by Lg : Gk → Gk, h 7→ gh, the left translation by g ∈ Gk. The tangent map of Lg is then given as ThLg : ThGk → TghGk. We now dene a left invariant ve tor eld on Gk as follows. A ve tor eld X : Gk → TGk, h 7→ X(h) is alled left invariant if ThLg(X(h)) = X(gh), ∀h ∈ Gk. (4.9) 70 Re ognizing the Lie algebra Lk as the tangent spa e at identity e of Gk (e = {f ∈ Gk|f(θ) ≡ eG}, where eG is the identity element of G), i.e. Lk = TeGk, we an dene a left invariant ontrol system as, dg(t) = TeLg(t)(ξ(t)) = g(t) · ξu(t), (4.10) dt where a given ontrol input u(t) determines a ontrolled ve tor ξu(t) in the Lie algebra Lk. Note that the loop algebra Lk an be identied with the tensor ∞ 1 produ t spa e g ⊗ F , where F is the ring of real valued C fun tions on S . Choose a basis of g as {A1, A2, . . . An}. Then, any ξ ∈ Lk an be written as, ξ(θ) = ξ1(θ)A1 + · · ·+ ξn(θ)A 1n, θ ∈ S , where ea h of ξk's (k = 1, . . . , n) are smooth fun tions on the ir le. We will now limit ourselves to the study of ontrol ve tors ξu of the form, ξu(t) = u1(t)A1 + · · ·+ um(t)Am, (4.11) where m < n and the ontrol input u(·) = (u1(·), . . . , um(·)) belongs to the set U of pie ewise ontinuous U valued fun tions, where U is ve tor spa e of Rm valued smooth fun tions on the ir le, i.e. U := {u(·) : u is pie ewise ontinuous in t, u(t) ∈ U = C∞(S1;Rm)}. 4.3 Controllability Having onstru ted the ontrol system on the loop group Gk, it is natural to ask the question of ontrollability or a essibility, i.e. given any two points g1 and 71 g2 in Gk, if they an be onne ted by a pie ewise dierentiable urve, onsisting of possibly nitely many pie es, ea h pie e being an integral urve of a left in- variant ve tor eld dened by hoosing a ontrol u(·) ∈ U . In nite dimensional 1 analogue of this question, i.e. where we shrink the ir le S down to a point, the ontrollability question is answered by the well known Chow-Rashevsky theorem [Wei-Liang, 1939; Rashevsky, 1938℄. In innite dimensional ases, however, it is not immediate if the Chow-Rashevsky theorem remains valid. There is a body of literature [Heintze and Liu, 1999; Salehani and Markina, 2014℄ that attempts to atta k this problem. It is the result of [Heintze and Liu, 1999℄ that we use in this se tion. This result addresses the ontrollability question in a weaker sense whi h we will make expli it. Let M be a omplete onne ted Hilbert manifold and let X(M) denote the set of all smooth ve tor elds dened on M. Let F ⊂ X(M) be a given family of smooth ve tor elds onM. LetRF(x) be the set of points inM that an be joined from x ∈ M by means of a pie ewise dierentiable urve, ea h pie e of whi h is an integral urve of a ve tor eld in F . Let Lie F be the Lie subalgebra of X(M) generated by F , and Liex F = {X(x) : X ∈ Lie F} - the evaluation of Lie F at x ∈ M. If M is nite dimensional, the lassi al Chow-Rashevsky theorem holds: if Liex F = TxM for ea h x ∈ M, then RF (x) = M, for every x ∈ M. In a general Hilbert manifold M, the following generalized Chow-Rashevsky theorem holds: 72 Theorem 4.3.1 ([Heintze and Liu, 1999℄). Let M be a omplete onne ted Hilbert manifold and F a family of smooth ve tor elds dened on M. If Liex F is dense in TxM for all x ∈ M, then RF (x) is dense in M for all x ∈ M. Theorem 4.3.1 is a weaker statement than the one for nite dimensional ase. Here we make pre ise the strong and weak notions of ontrollability. Consider the ontrol system onstru ted in (4.10)(4.11). Note that the loop group Gk an be given a stru ture of a smooth Hilbert manifold [Eells Jr, 1966; Ebin and m Marsden, 1970℄. In this ase, the family F ∈ X(Gk) is given by {Xi}i=1 , where Xi(g(t)) = g(t) · (ui(t)Ai), for g(t) ∈ Gk. Denition 4.1. (Strong Controllability) The ontrol system (4.10)(4.11) is said to be strongly ontrollable if RF = Gk, i.e. given any two points g1, g2 ∈ Gk, we an nd a ontrol input that will transfer the system from g1 to g2. Denition 4.2. (Weak Controllability) The ontrol system (4.10)(4.11) is said to be weakly ontrollable if RF is dense in Gk, i.e. given any two points g1, g2 ∈ Gk, we an nd a ontrol input that will transfer the system from g1 to a state that is arbitrarily lose to g2. The set {A1, ..., Am} is said to be bra ket generating if the iterated bra kets of its elements span the Lie algebra g. In the nite dimensional analogue of the ontrol system dened in (4.10)(4.11), the (strong) ontrollability ondition m a ording to Chow-Rashevsky theorem is equivalent to having the set {Al}l=1 bra ket generating in g. We will now try to establish (weak) ontrollability of 73 the innite dimensional loop group ase by assuming that {A1, ..., Am} is bra ket generating in g. Theorem 4.3.2. Consider the ontrol system (4.10)(4.11) on the loop group Gk. m Assume that the set {Al}l=1 is bra ket generating in g. Then the system is weakly ontrollable. Proof. The denition of left invariant ve tor elds on Gk (4.9) an also be made expli it by means of smooth fun tions on Gk. Let D be the set of smooth real val- ued fun tions on Gk. Then given an element ξ ∈ Lk, we an dene a dierentiable ve tor eld Xξ : D → D as, (Xξf)(g) = (Df)g · gξ, f ∈ D, (4.12) whereD denotes the dierential operator. Given two ve tor eldsXξ, Xη ∈ X(Gk), we an al ulate their Ja obi-Lie bra ket dened as, [Xξ, Xη]f = Xξ(Xηf)−Xη(Xξf), f ∈ D. We ompute, Xξ(Xηf)(g) = (Xξ((Df)g · gη))(g) = D((Df)g · gη)g · gξ = (D2f)g · (gη, gξ) + (Df)g · (D(gη)g · gξ) = (D2f)g · (gη, gξ) + (Df)g · (gξη). Similarly, Xη(X 2 ξf)(g) = (D f)g · (gξ, gη) + (Df)g · (gηξ). 74 2 The symmetry of the se ond dierential operator D yields [Xξ, Xη]f(g) = (Df)g · (g(ξη − ηξ)) = X[ξ,η]f(g), where [ξ, η] is the usual (pointwise) Lie bra ket on Lk. This leads us to a detailed study of the Lie bra ket of the loop algebra Lk. It is immediate that Lk is mr generated by the generators {Pr }m ∈Z, r∈{1,...,n} dened as,r Pmr imrθr = e Ar. (4.13) r Let the stru ture onstants of g be denoted by Γpq. Then, n ∑ [Pmp, Pmq ] = Γr mp+mqp q pqPr . (4.14) r=1 With this notation, a ontrolled ve tor ξ(t) = ξu(t) ∈ Lk an be expressed as ( ) m ∑ ∑ ξ mr mru(t) = ur (t)Pr , (4.15) r=1 mr∈Z mr where for ea h r = 1, ..., m, ur (t) ∈ C's are the Fourier oe ients of the ontrol ur(t) and umrr (t) = u −mr r (t) (sin e the ontrols are real). We now dene a family mr of ve tor elds on Gk as F = {Xr }m ∈Z, r∈{1,...,m}, wherer Xmrr f(g) = (Df) · gPmrg r , f ∈ D. Taking bra ket of any two ve tor elds from the family F yields another ve tor led whi h is governed by the ommutator relationship in (4.14). Note that sin e the set {A1, ..., Am} is bra ket generating in g, for ea h l ∈ {m+ 1, ..., n}, we are guaranteed to get the item ∑ a m Pml := P r r rl l , a ∈ +r Z ∪ {0}, mr ∈ Z, r ∈ {1, ..., m}, 75 at some depth of iterated bra kets from the family F . By hoosing mr ∈ Z, r ∈ {1, . . . , m}, we an then a hieve any ml ∈ Z. We have thus proved that if the set {Ar}r∈{1,...,m} is bra ket generating in g, Lieg F is dense in the tangent spa e TgGk at ea h g ∈ Gk. The generalized Chow's theorem 4.3.1 then provides the required (weak) ontrollability result.  4.4 Optimal Control Problems We start with the left invariant ontrol system on the loop group Gk as in (4.10)- (4.11). Now for a given T > 0, onsider the following xed end-point optimal ontrol problem, ∫ T min J(u) = L(g(t), u(t))dt u∈U ( 0PG) (4.16) subje t to: ġ = g · ξ, g(0) = g0, g(T ) = gT ; g0, gT ∈ Gk. We are interested in deriving ne essary onditions for optimality for su h optimal ontrol problems. Spe ial are needs to be taken sin e the problem is posed in an innite dimensional setting. We provide two dierent approa hes for doing that. 4.4.1 Cal ulus of Variations Approa h ∞ 1 r Let x(t) = (x1(t), ..., xr(t)) ∈ C (S ,R ) =: X denote a ve tor that an be used to represent the omponents of g(t) ∈ G, for some r ≥ n. The group dynami s ġ = g · ∑ξ(u) = g · ( i uiAi) ee tively lets us write the ontrol u(t) as a fun tion of (x(t), ẋ(t)). The xed endpoint onstraints in g an be translated to some 76 nonholonomi onstraints of the form ∫ T Φ(x, ẋ) = φ(x(t), ẋ(t))dt = 0, (4.17) 0 where φ(x(t), ẋ(t)) ∈ C∞(S1,Rl) =: Z, for some l < r. Then the problem (PG) an be written as, ∫ T min J = L(x(t), ẋ(t))dt, (4.18) 0 subje t to the nonholonomi onstraints (4.17). Here we have to keep in mind that the variations in x are to be both in t and θ. This is a well known problem alled the `Lagrange problem' in al ulus of variations. The one-dimensional Lagrange problem is well studied [Gelfand and Fomin, 1963; Elsgol , 2012℄. However, the theory behind multidimensional problem is more ompli ated and less omplete [Giaquinta and Hildebrandt, 1996; Bliss, 1946℄. The di ulty arises sin e not all the ẋ are freely variable. A ording to [Giaquinta and Hildebrandt, 1996, p. 112℄, there exist a Lagrange multiplier λ ∈ C∞(R× S1;Rl), su h that we an nd the free extremals of the augmented Lagrangian in an usual way. Moreover, sin e the onstraints (4.17) are of isoperimetri type, λ does not depend on t [Rund, 1966, p. 349℄. We an write the augmented Lagrangian as L̃ = L+ 〈λ, φ〉Z , (4.19) ∞ 1 l where λ = (λ1, λ2, ..., λl) ∈ C (S ;R ). Furthermore, we are interested in a spe ial stru ture of the Lagrangian, namely ∫ 2π L(x, ẋ) = L(x, ẋ) dθ, (4.20) 0 77 where the fun tional L is alled the Lagrangian density. The augmented ost fun tion take the form, ∫ T ∫ 2π l ∫ 2π ∑ J̃(x, ẋ) = L(x, ẋ) dθdt+ λj(θ)φj(x, ẋ) dθ 0 0 j=1 0 ∫ T ∫ 2π = L̃(x, ẋ) dθdt, 0 0 where the redened Lagrangian density is given by, l ∑ L̃ = L+ λjφj. (4.21) j=1 We will relabel L̃ by L in subsequent analysis for onvenien e. By invoking ∂x ∂x ∂2x notations xt = , xθ = , xtθ = = x et ., we an write L = L(x, x , x∂t ∂θ ∂t∂θ θt θ t, xtθ). In order to optimize this ost fun tional, the variational prin iple requires that, ∫ T ∫ 2π r ( ) ∑ ∂L ∂L ∂L ∂L δJ̃ = δxi + δxi,θ + δxi,t + δxi,tθ dθdt = 0, 0 0 ∂xi=1 i ∂xi,θ ∂xi,t ∂xi,tθ where δy denotes variation of the quantity y that vanishes at the endpoints of t and θ. Using integration by parts, for ea h i, we may write, ∣2π ∫ 2π ∂L ( )∫∂L ∣ 2π− ∂ ∂L∣δxi,θdθ = δxi δxidθ∣ 0 ∂xi,θ ∂xi,θ ∣ 0 ∂θ ∂xi,θ 0 ∫ 2π ( ) − ∂ ∂L= δxidθ, 0 ∂θ ∂xi,θ Similarly, ∫ T ∂L ∫ T (∂ ∂L ) δxi,tdt = − δxidt, 0 ∂xi,t 0 ∂t ∂xi,t ∫ T ∫ 2π L ∫ T ∫ 2π ( L )∂ ∂2 ∂ δxi,tθdθdt = δxi dθdt. 0 0 ∂xi,tθ 0 0 ∂t∂θ ∂xi,tθ 78 Hen e, T 2π ( L ( L ) (∫ ∫ ∑ ∂ − ∂ ∂ − ∂ ∂L ) ( )) ∂2 ∂L δJ̃ = + δxi dθdt 0 0 ∂xi ∂θ ∂xi,θ ∂t ∂xi,t ∂t∂θ ∂xi i,tθ The Euler-Lagrange equations an then be expressed as, ( ) ( ) ( ) ∂L ∂ ∂L ∂ ∂L ∂2− − ∂L+ = 0, i = 1, 2, ..., r. (4.22) ∂xi ∂θ ∂xi,θ ∂t ∂xi,t ∂t∂θ ∂xi,tθ 4.4.2 Maximum Prin iple Approa h In this se tion we provide a brief exposure to Pontryagin's Maximum Prin iple (PMP) type argument in innite dimensional spa es. It is to be noted that PMP does not automati ally hold in general innite dimensional optimal ontrol prob- lems, one requires some more assumptions for it to work. A detailed study on this subje t is done in Appendix A. Here we only dene some notations and assumptions to state the ne essary theorem. We onsider an abstra t dierential equation in a Hilbert spa e X , dx(t) = f(t, x(t), u(t)) a.e. in [0, T ], (4.23) dt where x(t) ∈ X , u(·) ∈ U , and T > 0. Here X is alled the state spa e and U is the set of all measurable fun tions u(·) : [0, T ] → U , where U is a separable metri spa e alled the ontrol spa e. With this setup, we formulate the following optimal ontrol problem (P), ∫ T min J(u) = L(t, x(t), u(t))dt u∈U 0 (P) (4.24) subje t to: ẋ = f(t, x, u) a.e. in [0, T ], x(0) = x0, x(T ) = xT . 79 We assume that both the fun tions f(·, ·, ·) and L(·, ·, ·) are Bo hner integrable in t ∈ [0, T ] and Lips hitz ontinuous in x(t) ∈ X , with onstant K. Furthermore, ′ we require the existen e and ontinuity of the Fré het derivatives fx(t, x, u) and L′x(t, x, u) ′ ′ . We also assume the fun tions f, L and their derivatives fx, Lx to be bounded, i.e. there exists an M > 0, su h that ‖f(t, x, u)‖ ≤ M, ‖f ′x(t, x, u)‖ ≤ M, ‖L(t, x, u)‖ ≤ M, ‖L′x(t, x, u)‖ ≤ M, for all (t, x(t), u(t)) ∈ [0, T ]×X×U . Note that these hypotheses ensure a ontin- uous and unique solution of (4.23) to exist [Avez, 1986℄. The following te hni al details is one of the key ingredients in the proof of the PMP. Denition 4.3. (Finite Codimensionality) [Fattorini, 1987℄ A subset S of a Hilbert spa e Z is alled to be nite odimensional in Z, if there exists a losed subspa e Zc ⊆ Z of nite odimension su h that Sc = Πc( o(S)), has nonempty interior in Zc, where Πc denotes the orthogonal proje tion from Z onto Zc and o means losed onvex hull. We will now make a key assumption to derive a nontrivial maximum prin iple. ∗ Let a solution of problem (P) exist and the optimal ontrol be denoted by u ∈ U ∗ and let the orresponding optimal traje tory be denoted as x (·). Then dene the `rea hable set' as, { ∫ t ′ ∗ ∗ R := z(T ) ∈ X | z(t) = fx(s, x (s), u (s)) · z(s)ds 0 } ∫ t + (f(s, x∗(s), v(s))− f(s, x∗(s), u∗(s))) ds, for some v(·) ∈ U (4.25) 0 80 (A1) The set R is nite odimensional in X . Remark 4.1. In general, it is not lear whether there exists a relationship between ontrollability (strong or weak) of the system and the nite odimensionality assumption of the `rea hable set' R. We will, however, prove that in a spe ial ase of G = H(3), the Heisenberg group, the strong ontrollability implies nite odimensionality of R. It is of future onsideration to address this question in a general ase. Using usual formalism, we invoke the pre-Hamiltonian fun tion H : R×X × U ×R×X∗ → R as, H(t, x(t), u(t), p0, p(t)) = p0L(t, x(t), u(t)) + 〈p(t), f(t, x(t), u(t))〉 , (4.26) p(t) ∈ X∗where is alled the ostate variable. Then the PMP an be written as, ∗ Theorem 4.4.1. (Maximum Prin iple) Let u ∈ U be an optimal ontrol for ∗ problem (P) and x (t) be the orresponding optimal traje tory. Then, there exist (p∗a pair 0, p ∗(t)) ∈ R × X∗, t ∈ [0, T ] ∗ ∗ ∗ ∗, su h that (p0, p ) ≡6 (0, 0), p0 ≤ 0, p (·) satises the dierential equation, ṗ∗(t) = − (f ′ (t, x∗ ⋆x (t), u∗(t)) p∗(t)− p∗0L′x(t, x∗(t), u∗(t)), (4.27) ⋆ where by A we denote the adjoint operator of the operator A. The pointwise maximization of the pre-Hamiltonian holds, H(t, x∗(t), u∗(t), p∗0, p ∗(t)) = maxH(t, x∗(t), v, p∗, p∗(t)), (4.28) v∈ 0U 81 for a.e. t ∈ [0, T ] ∗ ∗. Moreover, x and p satisfy Hamilton's anoni al equations, i.e. dx∗ δH = ∗ (t, x ∗, u∗, p∗ ∗ dt δp 0 , p ) (4.29) dp∗ δH = − ∗ (t, x ∗, u∗, p∗0, p ∗). (4.30) dt δx A proof of this theorem is rather ompli ated and is given in Appendix A. We now use this result to state a maximum prin iple for the loop group ase. ∗ Theorem 4.4.2. (Maximum Prin iple in loop group setting) Let u ∈ U ∗ be an optimal ontrol for problem (PG) and g (t) be the orresponding optimal traje tory. Assume the nite odimensionality ondition (A1). Denote Rk, the ∞ 1 Hilbert spa e of k-Sobolev ompletion of the spa e R = C (S , gl(r,R)), for some r > n (p∗, p∗(t)) ∈ R×R , t ∈ [0, T ] (p∗, p∗. Then, there exist a pair 0 k , su h that 0 ) 6≡ (0, 0), p∗0 ≤ 0, p∗(·) satises the dierential equation ṗ∗(t) = −p∗(t) · ξ(u∗(t))T, (4.31) and the pointwise maximization of the pre-Hamiltonian holds, H(g∗(t), u∗(t), p∗, p∗0 (t)) = maxH(g ∗(t), v, p∗, p∗(t)), (4.32) v∈ 0U ∗ ∗ for a.e. t ∈ [0, T ]. Moreover, g and p satisfy Hamilton's equations, i.e. dg∗ δH = ∗ (g ∗, u∗, p∗0, p ∗) dt δp (4.33) dp∗ −δH= ∗ (g ∗, u∗, p∗, p∗). dt δg 0 Proof. It is almost immediate that under a nite odimensionality assumption like (A1), we an state a maximum prin iple like Theorem 4.4.1 for problem (PG). The 82 only aveat is that the state spa e Gk is not a Hilbert spa e and hen e Theorem 4.4.1 annot be applied dire tly. However, adopting an `enlargement' te hnique [Bro kett, 1973; Justh and Krishnaprasad, 2015a℄, we an state an analogous maximum prin iple. We re ognize that the loop group Gk is a subset of Rk. The spa e Rk an then be regarded as the `raised' state spa e. The dynami s (4.10), 1 along with the initial ondition g(0, θ) = g0(θ) ∈ G for all θ ∈ S , ensures that g(t, θ) remains in G for all (t, θ) ∈ [0, T ] × S1. Endow the spa e gl(r,R) with ( ) the tra e inner produ t and an indu ed norm, i.e. 〈A,B〉 (r, ) = Ttr A B andgl R √ ‖A‖ Tgl(r) = tr (A A), for A,B ∈ gl(r,R). We now dene the pre-Hamiltonian H : Rk × U ×R× Rk → R as, H(g(t), u(t), p0, p(t)) = 〈p(t), g(t)ξ(u(t))〉 + p0L(u(t)), (4.34)Rk where the duality pairing in the denition of H an be expli itly written as, 2π k 〈 〉∫ i ∑ 〈 d d i p(t), g(t)ξ(u(t))〉 = p(t, θ), (g(t, θ)ξ(u(t, θ))) dθ Rk dθi i0 dθi=0 gl(r,R) ∫ 2π k ( )i ∑ d di = tr p(t, θ)T · (g(t, θ)ξ(u(t, θ))) dθ. 0 dθ i dθi i=0 ∗ We are now all set to apply Theorem 4.4.1. If u ∈ U is an optimal ontrol, then we have, H(g∗(t), u∗(t), p∗, p∗0 (t)) = maxH(g ∗(t), v, p∗, p∗0 (t)). (4.35) v∈U δH (g∗, u∗, p∗, p∗) = g∗ξ(u∗) = ġ∗It is obvious that 0 . We an also derive for anyδp∗ g̃(t) ∈ Rk, (suppressing other arguments) δH · H(g ∗ + αg̃)−H(g∗) 〈 〉 ∗ g̃ = lim = 〈p ∗, g̃ξ(u∗)〉 = p∗ξ(u∗)T, g̃ , δg α→0 α Rk Rk 83 whi h implies the adjoint equation to (4.10) is, ṗ∗(t) = −p∗(t) · ξ(u∗(t))T. (4.36)  4.5 Spe ial Case : G = H(3) The previous se tion on optimal ontrol provides a on rete foundation in whi h we an state the maximum prin iple for the onsidered optimal ontrol problem. In this se tion, we will explore a spe ial ase where we take the Lie group, G as the Heisenberg group, H(3). Note that the nite number of parti les ase of this problem has been onsidered in [Justh and Krishnaprasad, 2016℄ and hen e this work an be thought as a ontinuum ounterpart of it. In H(3), g(t, θ) ∈ H(3) has the stru ture,   1 x1 x x1x2  3 + 2      g = 0 1 x  ,  2      0 0 1 that satisfy the group evolution equation, ∂ g(t, θ) = g(t, θ) · (u1(t, θ)A1 + u2(t, θ)A2) , (4.37) ∂t where,     0 1 0 0 0 0             A 1 = 0 0 0 , A2 = 0 0 1 ,             0 0 0 0 0 0 84 along with A3 = [A1, A2], form an orthonormal basis for the asso iated Lie algebra h(3). We attempt to address the optimal ontrol problem formulated before, under this Heisenberg group setting, i.e. we solve the following, ∫ T ∫ T ∫ 2π min L(u)dt = L(u)dθdt 0 0 0 ( ( )) T 2π ( )2 ( )∫ ∫ 21 ( )2 2 ∂u1 ∂u2= u1 + u2 + χ + dθdt,2 0 0 ∂θ ∂θ (4.38) where L is alled the Lagrange density fun tion. 4.5.1 Controllability It is a dire t exer ise of the generalized Chow-Rashevsky theorem 4.3.1 to show ∞ 1 weak ontrollability in this ase. The loop algebra C (S , h(3)) is spanned by {eim1θA1, eim2θA2, eim3θA3}. The family F of left invariant ve tor elds that is m1 m2 hosen by means of ontrol inputs is given by F = {X1 , X2 }, where Xmrr f(g) = (Df) · geimrθg Ar, f ∈ D, r = 1, 2. Sin e the only non-vanishing bra kets in h(3) are [A1, A2] = A3 = −[A2, A1], the ∞ 1 set Lieg F would span the tangent spa e at every point g ∈ C (S ,H(3)). However, we an provide an argument that establishes the strong ontrol- lability in this ase. Here we des ribe an approa h to onstru t a andidate 0 T smooth ontrol given any endpoint onditions xi (θ) and xi (θ), i = 1, 2, 3. With- 0 out loss of generality, we may assume xi (θ) = 0 for all i. Now sin e, xi(t, θ) = 85 ∫ t ui(t, θ)dt, i = 1, 2, we an hoose smooth ontrols vi(t, θ), t ∈ [0, t̄ ], i = 1, 2, for0 some t̄ < T , su h that x1 and x2 rea h their nal endpoints. At time t = t̄, let the T `error' in x3 variable be denoted as ∆x3(θ) = x3 (θ)− x3(t̄, θ). Note that without loss of generality, we may assume that ∆x3(θ) > 0 for all θ. We know that in a single nonholonomi integrator ase, if we omplete a loop in time for (x1, x2) variables, the hange in x3 variable will be given by the area of the loop. We may use the same idea in the ontinuum ase to onstru t smooth ontrols. We an generate the following ir ular loops (in time) in (x1, x2) variables,       x1(t, θ) x T 1 (θ)− r(θ) cos(ω(t− t̄))        =  + r(θ)  ,       x T2(t, θ) x2 (θ) sin(ω(t− t̄)) (4.39) 2π ω = − , t ∈ (t̄, T ],T t̄ where r(·) is a smooth fun tion in θ. The ontrols required to generate these loops are given by     ṽ1(t, θ) − sin(ω(t− t̄))      = r(θ)ω  , (4.40)     ṽ2(t, θ) cos(ω(t− t̄)) We then ompute, ∫ 1 T x3(T, θ) = x3(t̄, θ) + (x1(t, θ)ṽ2(t, θ)− x2(t, θ)ṽ1(t, θ))dt 2 t̄ = xT3 (θ)−∆x3(θ) + πr2(θ). Sin e ∆x3(·) is smooth, we an always hoose smooth fun tion r(·) su h that 86 ∆x3(θ) = πr 2(θ), and thus the smooth ontrols    vi(t, θ), t ∈ [0, t̄ ] ui(t, θ) = , i = 1, 2, (4.41)   ṽi(t, θ), t ∈ (t̄, T ] make the required state transitions possible. This shows that in the Heisenberg ase we have strong ontrollability. Remark 4.2. We an assume ∆x3 > 0 for all θ be ause if it was not the ase, we ould add another pie e of ontrols v̂i's before applying the ontrols ṽi's. The purpose of the ontrols v̂i's would be to make the states (x1, x2) undergo a ir ular 2 loop of radius r̂, for all θ. This will produ e a hange in the x3 variable by πr̂ 2 for all θ. Hen e the new error an be written as ∆x3(θ) = ∆x3(θ) + πr̂ . We an always pi k a r̂ so that ∆x3(θ) > 0 for all θ. 4.5.2 Equations of Optimal Control Ne essary onditions for optimality an be derived by various methods. We will present two su h approa hes to solve the optimal ontrol problem (4.38). 87 4.5.2.1 Cal ulus of Variations Approa h Note that the group dynami s in (4.37) an also be expressed by the dierential equations, ∂x1 = u1, (4.42) ∂t ∂x2 = u2, (4.43) ∂t ∂x3 1 = (x1u2 − x2u1). (4.44) ∂t 2 ∫ T Sin e we have an integral onstraint, namely, ∆x3(θ) = 1 (x u 2 0 1 2 − x2u1)dt, we invoke a Lagrange multiplier, λ ∈ C∞(S1,R) to write the augmented Lagrangian density fun tion as, 1 [( ) ( ) ]L = x2 21,t + x2,t + χ x2 21,θt + x2,θt + λ (x1x2,t − x2x1,t) . (4.45)2 2 where we adopted the notation onventions x = ∂xii,θ , x ∂ xi ∂θ i,tθ = = x ∂t∂θ i,θt et . Appli ation of the Euler-Lagrange equations (4.22) gives the following equations, ∂2x1 ∂x2 ∂ 4x1 = λ + χ ∂t2 ∂t ∂θ2∂t2 (4.46) ∂2x 42 − ∂x1 ∂ x2= λ + χ , ∂t2 ∂t ∂θ2∂t2 whi h yields the evolution equations for optimal ontrols, A∂u1 = λu2 ∂t (4.47) A∂u2 = −λu1, ∂t where we denote A := (1− χ∆), ∆ being the Lapla ian. ∞ 1 Remark 4.3. A is a positive denite self adjoint operator in C (S ,R), having eigenvalues αn = 1 + χn 2 inθ with asso iated eigenve tors en = e for n ∈ Z. 88 Lemma 4.5.1. The quantities, ( ) ∫ 2π ∂nu ∂nu ∂n n1A 1 u2A∂ u2hn = + dθ, n = 0, 1, 2, ... (4.48) ∂tn ∂tn ∂tn ∂tn0 are onserved along any optimal traje tories satisfying (4.47). dhn Proof. It is easy to establish that for ea h n, = 0, whi h follows dire tly from dt the way optimal ontrols behave in (4.47) and the fa t that the operator A is self adjoint.  4.5.2.2 PMP Approa h ∞ 1 3 We introdu e the ostate variable p(t) = (p1(t), p2(t), p3(t)) ∈ C (S ;R ), t ∈ [0, T ]. The pre-Hamiltonian an be written as ( onsidering only normal extremals, i.e. where p0 6= 0 and an be normalized to −1), H = 〈ẋ, p〉 − L 2π ( )∫ 1 1 ( ) χ ( ) = u1p1 + u2p2 + (x1u2 − x 2 2 2 22u1)p3 − u + u − u + u dθ. 0 2 2 1 2 2 1,θ 2,θ (4.49) Remark 4.4. Note that the nite odimensionality assumption is satised in this ase. To see this, note that the omponents of members of the `rea hability set' (4.25) an be expressed as, ∫ T z1(T ) = w1(t)dt 0 ∫ T z2(T ) = w2(t)dt 0 ∫ 1 T ∫ 1 T z3(T ) = (z1(t)u ∗ 2(t)− z2(t)u∗1(t)) dt+ (x∗1(t)w ∗2(t)− x2(t)w1(t)) dt,2 0 2 0 89 where w(t) = (w1(t), w2(t)) is any arbitrary ontrol input. If we denote traje tories u orresponding to any input u(·) as x (t) = x(t, u(·)) z (T ) = xw, then, 1 1 (T ), z2(T ) = xw2 (T ). Furthermore note that, [∫ ∫ 1 T T z3(T ) = (x ∗u∗1 2 − x∗ ∗2u1)dt+ (xw w1 w2 − x2 w1)dt2 0 0 ∫ T ] − ((x∗ − xw)(u∗ − w )− (x∗ w1 1 2 2 2 − x2 )(u∗1 − w1))dt 0 = x∗3(T ) + x w 3 (T )− x̃3(T ), ∗ w ∗ where x̃ = x − x . We may now hoose w = u . This makes x̃(T ) = 0 and (z1(T ), z2(T ), z3(T )) = (x ∗ 1(T ), x ∗ 2(T ), 2x ∗ 3(T )). Sin e the Heisenberg ase is strongly ontrollable, the `rea hability set' R spans the whole of the state spa e X , making it trivially nite odimensional in X . We an now dire tly apply Theorem 4.4.1 to derive ne essary optimality on- ditions. The maximum prin iple would require us to maximize (4.49) pointwise over the ontrols, i.e. we are attempting to nd the Hamiltonian as, H(x, p) = sup H(x, {v }2i i=1, p) v ∞ 11,v2∈C (S ;R) ( ∫ 2π 1 = sup v1p1 + v2p2 + (x1v2 − x2v1)p3 v ,v ∈C∞ 1 21 2 (S ;R) 0 ) 1 ( ) ( )− 2 2 − χv1 + v 2 22 v1,θ + v2,θ dθ (4.50)2 2 This maximization results in two Euler-Lagrange equations, ( ) ∂H − ∂ ∂H = 0, i = 1, 2, (4.51) ∂vi ∂θ ∂vi,θ 90 that yield the optimal ontrols, Au1 = p1 − 1 x2p3, 2 (4.52) A 1u2 = p2 + x1p3, 2 with the usual denition of A = (1− χ∆). The Hamiltonian an be read as, ∫ 1 2πH = (u1Au1 + u2Au2) dθ 2 0 [( ) ( ) ∫ 1 2π − 1 1= p1 x −12p3 A p1 − x2p3 2 0 2 2 ( ) ( )] 1 + p2 + x1p3 A−1 1 p2 + x1p3 dθ (4.53) 2 2 The dynami s of the ostate variable p an be al ulated from Hamilton's equa- ∂p = − δH δHtion, , where denotes the fun tional derivative of H with respe t ∂t δx δx to x. Note that the Hamiltonian fun tion H is smooth in x, so we an take this 〈 〉 δH derivative. This an be dened as, DH(xi) · σ = , σ , i = 1, 2, 3, whereδxi DH(xi) · σ is the Fré het derivative of H at xi in the dire tion of σ. This is dened as, ∣ d ∣ DH(x ∣i) · σ = H(x1 + ǫσ) . dǫ ∣ǫ=0 We may al ulate, for i = 1, ∫ 1 d 2π [( ) ( ) H · − 1D (x ) σ = p −1 11 1 x2p3 A p1 − x2p3 2 dǫ 0 2 2 ( ) ( )] ∣ 1 1 ∣ + p2 + (x1 + ǫσ)p A−13 p2 + (x1 + ǫσ)p ∣3 dθ 2 2 ∣ǫ=0 ( ( )) ∫ 1 2π 1 = σp A−13 p2 + (x1 + ǫσ)p3 dθ 2 0 2 2π ( ( ))∫1 = σp3A−1 1 p2 + x1p3 dθ 2 0 2 〈 ( ) 〉 p3 = A−1 1p2 + x1p3 , σ , (4.54) 2 2 91 so that we may write, H ( )δ p3A−1 1= p2 + x1p3 . (4.55) δx1 2 2 ( ) δH We an similarly al ulate = −p3A−1 p − 11 x2p δH3 and = 0. The evolutionδx2 2 2 δx3 of p an then be expressed as, ( ) ∂p1 −p3= A−1 1p2 + x1p3 ∂t 2 2 ( ) ∂p2 p3 = A−1 1p1 − x2p (4.56)3 ∂t 2 2 ∂p3 = 0. ∂t From (4.52) and (4.56), we noti e that, ∂ Au1 = −p A−13 u2 ∂t (4.57) ∂ Au −12 = p3A u1 ∂t Hen e, re ognizing p3 as the negative of Lagrange multiplier λ in the previous se tion, we redis over (4.47). 4.5.3 Behavior of Optimal Control Denote, z(t, θ) := u1(t, θ) + iu2(t, θ). Then, (4.47) an be expressed as, A∂z = −iλz. (4.58) ∂t Sin e u1, u2 are periodi fun tions in θ with period 2π, they have a Fourier series ∞ ∑ (ν) (ν) representation, uν(t, θ) = un (t)e inθ, ν = 1, 2, where un 's are the Fourier n=−∞ 92 oe ients of uν . Moreover, sin e uν 's are real valued fun tions of θ, we have, (ν) (ν) u−n = un . Transforming (4.58) in Fourier domain, we get, dž A = −iΛž, (4.59) dt ( ∞ ) ∞ (1) where we denote, A = diag {αn}n=−∞ , ž(t) = [zn(t)]−∞. zn(t) = un (t) + (2) iun (t), i.e. zn(t)'s are Fourier oe ients of z(t, θ). Λ is the innite Toeplitz matrix generated by the Fourier oe ients of λ, i.e.   λ λ λ  0 −1 −2 · · ·      λ   1 λ0 λ−1 · · · Λ =   (4.60)    λ λ λ · · ·  2 1 0      · · · · · · · · · · · · Sin e λn's are Fourier oe ients of real valued fun tion λ, we have λ−n = λn. ∗ This leads to the observation that Λ is (innite) Hermitian matrix, i.e. Λ = Λ . 4.5.3.1 Trun ation of Fourier Coe ients Here we will onsider rst N+1 Fourier oe ients of z and provide an analysis of T (4.59) in the trun ated nite dimensional ase. We write, žN = (z−N · · · z0 · · · zN ) ∈ 2N+1 C . Then, the trun ated version of (4.59) an be written as, džN AN = −iΛN žN , (4.61) dt where AN and ΛN 's are appropriately trun ated matri es from A and Λ, respe - 1/2 −1 tively. Note that AN ≻ 0. Let us denote, ẑ = AN žN and BN = AN , then we get, 93 dẑ = −iΛ̂N ẑ, (4.62) dt ( )∗ 1/2 1/2 where Λ̂N = BN ΛNBN = Λ̂N . Sin e the matrix −iΛ̂N is skew Hermitian, all its eigenvalues are on the imaginary axis. Denote them by −iσn, σn ∈ R, for n = −N, ..., N . There exists a unitary matrix V that diagonalizes −iΛ̂N , i.e. −iΛ̂N = V ∗DV, D = diag({−iσn}). We perform another oordinate hange by, 1/2 z̃ = V ẑ = V AN žN , (4.63) whi h yields the de oupled equations, dz̃ = Dz̃, (4.64) dt i.e. dz̃n = −iσnz̃n, dt =⇒ z̃ (t) = e−iσntn z̃n(0), n = −N, · · · , N. (4.65) Performing the inverse Fourier operation, we see that, N ∑ z̃(t, θ) = z̃n(t)e inθ n=−N N ∑ =⇒ z̃(t, θ) = ei(nθ−σnt)z̃n(0). (4.66) n=−N 94 This is an equation of superposition of 2N + 1 traveling waves with n being the σn wave number and vn = is the speed of propagation asso iated with n-th moden of the wave. 4.5.3.2 Velo ity of Propagation We know that in the wave equation of (4.66), the velo ity vn of propagation σn orresponding to n-th frequen y is determined as, vn = , where σn's are (real)n 1/2 1/2 1/2 1/2 eigenvalues of the Hermitian matrix Λ̂N = BN ΛNBN . Sin e BN ΛNBN ∼ ΛNBN ∼ BNΛN (similar matri es), the eigenvalues of Λ̂N an be hara terized by those of BN and ΛN . N 1 1 We have eig(BN ) = {βn}−N , where βn = = 2 . Now, ΛN is a Toeplitzαn 1+χn 2N Hermitian matrix formed by the Fourier oe ients {λn}−2N . Given those oe- ients, it is in general not possible to write down losed form representation of its eigenvalues. However, the bounds of eigenvalues of su h a matrix is well known. We will make a little detour to state these results. 95 4.5.3.3 Toeplitz Matri es and Eigenvalues Let f be a periodi fun tion over the interval [0, 2π) and {fn} are its Fourier oe ients. Let us denote Tn(f), the n×n Toeplitz Hermitian matrix dened as,   f0 f−1 f−2 · · · f −(n−1)      f   1 f0 f−1 · · · f−(n−2)     T n(f) = f2 f1 f0 · · ·  . (4.67)      . . .   . . . . . .       f(n−1) f0 We also dene, mf = ess inf(f), (4.68) Mf = ess sup(f). (4.69) Let the eigenvalues of Tn(f) be denoted by τn,k, k = 1, ..., n. Then, mf ≤ τn,k ≤ Mf . (4.70) Note that, max |τn,k| ≤ max(|mf | , |Mf |) ≤ M|f |. We re all another useful result k here. Lemma 4.5.2. Let P,Q be Hermitian positive denite matri es of same order. If τ(X) denote eigenvalues of X, then, τmax(PQ) ≤ τmax(P ) · τmax(Q) (4.71) τmin(PQ) ≥ τmin(P ) · τmin(Q) (4.72) 96 4.5.3.4 Bounds on the Velo ity To get bounds on σn's, we need an useful assumption that will be apparent shortly. A-1 mλ > 0, i.e. ΛN = T2N+1(λ) ≻ 0. Sin e µn = τ(BNT2N+1(λ)), Lemma 4.5.2 gives the following bound, for n = −N, ..., N , τmin(BN)τmin(T2N+1(λ) ≤ σn ≤ τmax(BN)τmax(T2N+1(λ)), (4.73) ⇒ mλ= ≤ σ ≤ M . 1 + χN2 n λ A ordingly, the velo ity is bounded by, mλ ≤ ≤ Mλvn (4.74) n(1 + χN2) n Remark 4.5. The assumption A-1 an be extended to in lude the ase Mλ < 0, i.e. ΛN ≺ 0 as well. 4.5.3.5 Spe ial Cases 1. Case - I: Constant λ : λ(θ) = λ0 6= 0, a onstant. We may assume that λ0 > 0. In this ase, (4.59) an be expli itly solved and the solution an be expressed as, ∞ ( )βnλ0 ∑ in θ− t z(t, θ) = e n z −jλ0tn(0) + e z0(0). (4.75) n=−∞ n=6 0 λ0 Here the velo ities vn an be written as, vn = , n ∈ Z \ {0}.n(1+χn2) 97 2. Case - II: Band-limited λ : Here, we onsider only one frequen y omponent of λ, i.e. λ(θ) = λ0 + λ eiθ1 + λ e −iθ −1 . Let A-1 hold, i.e. ΛN = T2N+1(λ) ≻ 0. In this ase, ΛN is a tri-diagonal Toeplitz matrix. The eigenvalues of su h a matrix are known to take the following form. Lemma 4.5.3. ( ) kπ τk(ΛN) = λ0 + 2 |λ1| cos , k = 1, ..., 2N + 1. (4.76) 2N + 2 This, ombined with Lemma 4.5.2, we get the following bound, ( ) λ0 + 2 |λ1| cos (2N+1)π2N+2 ( )≤ σ ≤ λ + 2 |λ | cos πn 0 1 (4.77) 1 + χN2 2N+2 4.5.4 Strong Coupling Limit, χ → ∞ It is interesting to note that in the limit χ → ∞, the equations (4.47) take simple ∞ 1 form. To see this, note that for some z ∈ C (S ;R), we may express ∞ ∑ A−1z = βn 〈φn, z〉 φn, n=−∞ 1 −1 where βn = 2 are the eigenvalues of A , and φn = einθ. This implies,1+χn ( ) ∑ A−1 1lim z = lim 〈z, φ 〉φ χ→∞ χ→∞ n n1 + χn2 n = 〈z, 1〉 ∫ 2π =⇒ lim A−1z = zdθ. χ→∞ 0 98 t = 0 s t = 6.67 s t = 13.33 s t = 20 s 1 1 1 1 PSfrag repla ements 0.5 0.5 0.5 0.5 0 0 0 0 -0.5 -0.5 -0.5 -0.5 -1 -1 -1 -1 0 5 0 5 0 5 0 5 θ θ θ θ 3 3 3 3 2 2 2 2 1 1 1 1 0 0 0 0 -1 -1 -1 -1 -2 -2 -2 -2 -3 -3 -3 -3 -2 0 2 -2 0 2 -2 0 2 -2 0 2 x1 x1 x1 x1 Figure 4.1: Numeri al solution of (4.61) for experiment 1. Evolution of u1 (blue) and u2 (red) is given in the rst row; the ontrols form a simple traveling wave along the θ domain. The se ond row shows evolution of x1, x2 variables. With this realization, we may rewrite (4.47) in the strong oupling limit as, ∫ ∂u 2π1 = λu2dθ ∂t 0 (4.78) ∫ ∂u 2π2 = − λu1dθ. ∂t 0 It is lear from (4.52) that both u1 and u2 are independent of θ. Then the equations (4.78) an be equivalently written as u̇1 = λ̃u2 (4.79) u̇2 = −λ̃u1, 99 u1, ux 22 x2 x2 x2 t = 0 s t = 4 s t = 8 s 1 0 0 0.5 -0.5 -1 0 -1 -2 -0.5 -1 -1.5 -3 1 1 1 1 1 1 0 0 0 PSfrag repla ements 0 0 0 -1 -1 -1 -1 -1 -1 t = 12 s t = 16 s t = 20 s 0 0 0 -1 -2 -2 -2 -4 -4 -3 -6 -4 -6 -8 1 1 1 1 1 1 0 0 0 0 0 0 -1 -1 -1 -1 -1 -1 Figure 4.2: Evolution of x3 for experiment 1. The blue loop represents the ir le S1 while the height of ea h point on the red loop is given by the value of x3 at the orresponding point of θ. ∫ 2π where λ̃ = λdθ. Equations (4.79) are optimal ontrol evolution equations for 0 a single agent [Justh and Krishnaprasad, 2016℄. Thus in the strong oupling limit χ → ∞, the optimal ontrol solutions for the ontinuum of agents ollapses to that of a single agent. This is alled the syn hronization of the o k, where every agent in the o k behaves the same way. 100 t = 0 s t = 6.67 s t = 13.33 s t = 20 s 1 1 1.5 1.5 PSfrag repla ements 0.5 1 1 0.5 0.5 0 0.5 0 0 -0.5 0 -0.5 -0.5 -1 -0.5 -1 -1 -1.5 -1 -1.5 0 5 0 5 0 5 0 5 θ θ θ θ 3 3 33 2 2 22 1 1 11 0 0 00 -1 -1 -1-1 -2 -2 -2-2 -3 -3 -3-3 -2 0 2 -2 0 2 -2 0 2 -2 0 2 x1 x1 x1 x1 Figure 4.3: Numeri al solution of (4.61) for experiment 2. Evolution of u1 (blue) and u2 (red) is given in the rst row. The se ond row shows evolution of x1, x2 variables. 4.5.5 Simulation Results We simulate the evolution of optimal ontrols u1 and u2 governed by the linear partial dierential equations (4.47) by means of Fourier analysis as presented in the se tion 4.5.3. In parti ular, here we present the solutions of the trun ated ordinary dierential equations (4.61), where we only keep tra k of rst N + 1 Fourier oe ients of ea h variable. Note that λ is assumed to have less than N + 1 oe ients. We will now present numeri al solutions in dierent ases by 101 u1, ux 22 x2 x2 x2 t = 0 s t = 4 s t = 8 s 1 0 0 0.5 -1 -0.5 0 -2 -1 -0.5 -3 -1 -1.5 -4 1 1 1 1 1 1 0 0 0 PSfrag repla ements 0 0 0 -1 -1 -1 -1 -1 -1 t = 12 s t = 16 s t = 20 s 0 0 0 -2 -2 -4 -5 -4 -6 -6 -8 -10 1 1 1 1 1 1 0 0 0 0 0 0 -1 -1 -1 -1 -1 -1 Figure 4.4: Evolution of x3 for experiment 2. The blue loop represents the ir le S1 while the height of ea h point on the red loop is given by the value of x3 at the orresponding point of θ. varying the initial onditions and parameter values. For all the experiments pre- sented here, N is taken to be 30 and the nal time T is set as 20 se onds and four snapshots of the optimal ontrols are shown. The Hamiltonian is veried to be staying a onstant (up to ma hine pre ision) for all of the experiments. The evolution of state variables is also re orded. For all the experiments presented here, (x1, x2) is set to start from a unit ir le and x3 is initially zero for all values of θ. 102 PSfrag repla ements t = 0 s t = 6.67 s t = 13.33 s t = 20 s 4 1.5 1 1.5 1 0.5 1 3 0.5 0 0.5 2 0 -0.5 0 -0.5 -1 -0.5 1 -1 -1.5 -1 0 -1.5 -2 -1.5 0 5 0 5 0 5 0 5 θ θ θ θ PSfrag repla ements t = 0 s t = 6.67 s t = 13.33 s t = 20 s 4 3 1.5 2 1 2 1.5 3 0.5 1 1 0 2 0.5 0 -0.5 0 1 -1 -1 -1.5 -0.5 0 -2 -2 -1 0 5 0 5 0 5 0 5 θ θ θ θ PSfrag repla ements t = 0 s t = 6.67 s t = 13.33 s t = 20 s 4 4 3 4 3 3 3 2 2 2 2 1 1 1 1 0 0 0 0 -1 -1 -1 0 5 0 5 0 5 0 5 θ θ θ θ Figure 4.5: Numeri al solutions of (4.61) for experiment 3 for (a) rst row, χ = 0.1, (b) se ond row, χ = 1, ( ) third row, χ = 10. In all plots, u1 is shown in blue and u2 is shown in red. The speed of information propagation de reases as the oupling onstant in reases. 103 u1, u2 u1, u2 u1, u2 The simplest set of initial onditions u1(0, θ) = cos(θ), u2(0, θ) = sin(θ), λ(θ) = 5, χ = 1 generate a traveling wave solution, a ording to equation (4.75). This an be seen from Fig. 4.1. The states (x1, x2, x3) are also integrated from an initial ondition for (x1, x2) forming a ir le on the plane and x3 being identi ally zero for all θ. The evolution of (x1, x2) an be seen from Fig. 4.1. The shape of the (x1, x2) ir le did not hange, although its size varied over time. The evolution of x3 is shown in Fig. 4.2, whi h appeared to de rease steadily for all θ. Next, for experiment 2, we onsider a band-limited λ, i.e. λ(θ) = 5 + cos(θ). Keeping all other onditions same as in experiment 1, we get Fig. 4.3-4.4. Here both size and shape of the (x1, x2) ir le hanged over time. The value of x3 de reased in this ase as well but more asymmetri ally than in experiment 1. In experiment 3, we show how a lo alized disturban e gets spread in the ontinuum. For this experiment, we let the ontrol u1 is initially zero everywhere, u1(0, θ) = 0. However, u2 has a lo alized peak at a ertain spatial point. We took the example of a Gaussian form, 1 − (θ−π) 2 u2(0, θ) = 2ρ2√ e , (4.80) 2πρ2 with ρ = 0.1. We then plot the solutions in three dierent settings of χ values, χ = 0.1, 1, 10 in Fig. 4.5. λ is taken to be a onstant in all these ases, λ(θ) = 5. It is dis overed in the previous se tion that the speed of traveling wave de reases as χ in reases. This an be seen learly from Fig. 4.5 as the disturban e is seen to be not well propagated for higher values of χ. 104 4.6 A Continuum of Agents on the Plane As an extension of the problem onsidered in Se tion 3.3, we will explore the ase where we take the underlying group as the spe ial Eu lidean group, SE(2). This ase then an be seen as a ontinuum ounterpart of [Justh and Krishnaprasad, 2015b℄. Every g(t, θ) ∈ SE(2) an be represented as,   cosx3 − sin x 3 x1     g =  sin x3 cosx x  ,  3 2     0 0 1 with group evolution dynami s, ∂ g(t, θ) = g(t, θ) · (u1(t, θ)A1 + u2(t, θ)A2) , (4.81) ∂t where,     0 −1 0 0 0 1             A1 = 1 0 0 , A =  2 0 0 0 ,             0 0 0 0 0 0 along with A3 = [A1, A2], form a basis for the asso iated Lie algebra se(2). Note that sin e {A1, A2} is bra ket generating in se(2), similar argument as in the Heisenberg ase would provide (weak) ontrollability result in this ase as well. We then seek to minimize the same ost fun tional as in (4.38), ( ( )) T 2π ( )2 ( )∫ ∫ 21 J = (u2 + u2 ∂u1 ∂u2 ) + χ + dθdt 2 1 20 0 ∂θ ∂θ 105 4.6.1 Equations of Optimal Control 4.6.1.1 Cal ulus of Variations Approa h The system dynami s an be equivalently expressed as, ∂x1 = u2 cosx3 ∂t ∂x2 = u (4.82)2 sin x3 ∂t ∂x3 = u1. ∂t Note that the ontrols u1 and u2 an be expressed by the xj variables and their derivatives as, u1 = x3,t, u2 = x1,t cos x3 + x2,t sin x3. Sin e we have two integral onstraints in this ase, namely, ∫ T ∫ T ∆x1(θ) = u2 cosx3dt = (x1,t cosx3 + x2,t sin x3) · cosx3dt 0 0 ∫ T ∫ T ∆x2(θ) = u2 sin x3dt = (x1,t cosx3 + x2,t sin x3) · sin x3dt, 0 0 ∞ 1 we invoke Lagrange multipliers λ, µ ∈ C (S ;R) and the augmented Lagrangian density an be read as, L 1 [( ) ( )] = x2 2 23,t + (x1,t cos x3 + x2,t sin x3) + χ x3,tθ + (x1,tθ cos x3 + x 2 2,tθ sin x3) 2 + λ (x1,t cosx3 + x2,t sin x3) · cosx3 + µ (x1,t cosx3 + x2,t sin x3) · sin x3. (4.83) We re all the Euler-Lagrange equations from (4.22), L ( L ) ( ) ( )∂ ∂ ∂ 2− − ∂ ∂L ∂ ∂L+ = 0, ∂xi ∂θ ∂xi,θ ∂t ∂xi,t ∂θ∂t ∂xi,θt 106 for i = 1, 2, 3. We note the following quantities, u2 = x1,t cosx3 + x2,t sin x3, u2,θ = x1,tθ cosx3 + x2,tθ sin x3, ( ∂ ∂L ) = u2,t cos x3 − u1 (u2 sin x3 + λ sin 2x3 − µ cos 2x3) , ∂t ∂x1,t ( ∂2 ∂L ) = χ (u2,θθt cosx3 − x3,θ (u2,θt sin x3 + u2,θu1 cos x3) ∂θ∂t ∂x1,θt − sin x3 (u2,θθu1 + u2,θu1,θ)) , ( ) ∂ ∂L = u2,t sin x3 + u1 (u2 cosx3 + λ cos 2x3 + µ sin 2x3) , ∂t ∂x2,t ( ) ∂2 ∂L = χ (u2,θθt sin x3 + x3,θ (u2,θt cosx3 − u2,θu1 sin x3) ∂θ∂t ∂x2,θt +cosx3 (u2,θθu1 + u2,θu1,θ)) , ( ) ∂ ∂L = u1,t, ∂t ∂x3,t ( L )∂2 ∂ = χu1,θθt. ∂θ∂t ∂x3,θt Subsequently, the Euler-Lagrange equations for the Lagrangian (4.83) take the form, u2,t cos x3 − u1u2 sin x3 − u1 (λ sin 2x3 − µ cos 2x3) = χ (u2,θθt cosx3 − sin x3 (u2,θθu1 + u2,θu1,θ)− x3,θ (u2,θt sin x3 + u2,θu1 cosx3)) , u2,t sin x3 + u1u2 cos x3 + u1 (λ cos 2x3 + µ sin 2x3) = χ (u2,θθt sin x3 + cosx3 (u2,θθu1 + u2,θu1,θ) + x3,θ (u2,θt cosx3 − u2,θu1 sin x3)) , u1,t − χu1,θθt = −u2 (λ sin x3 − µ cosx3 − χu2,θx3,θ) , (4.84) 107 whi h after some readjustments yield, ∂ (1− χ∆)u1 = −u2 (λ sin x3 − µ cosx3 − χx3,θu2,θ) ∂t ∂ (1− χ∆)u2 = u1 (λ sin x3 − µ cosx3 − χx3,θu2,θ) (4.85) ∂t ∂ (λ sin x3 − µ cosx3 − χu2,θx3,θ) = −u1(1− χ∆)u2. ∂t Denoting A := (1 − χ∆) and µ3 := λ sin x3 − µ cosx3 − χu2,θx3,θ, we then an express (4.85) as, ∂ Au1 = −µ3u2 ∂t ∂ Au (4.86)1 = µ3u1 ∂t ∂µ3 = −u1Au2. ∂t 4.6.1.2 PMP Approa h ∞ 1 3 We introdu e the ostate variable p(t) = (p1(t), p2(t), p3(t)) ∈ C (S ;R ), t ∈ [0, T ]. The pre-Hamiltonian ( onsidering only normal extremals) an be written as, H = 〈ẋ, p〉 − L 2π ( )∫ 1 ( ) χ ( ) = u2 (p 2 2 2 2 1 cosx3 + p2 sin x3) + u1p3 − u + u − u 2 1 2 2 1,θ + u2,θ dθ. 0 (4.87) Remark 4.6. Here we only onsider normal extremals, i.e. when p0 6= 0 and an be normalized to −1. Note that in this ase, the emptiness (i.e. the full ostate being identi ally zero) of the PMP would not o ur. It is of future eort to investigate whether the nite odimensionality ondition is satised in the SE(2) ase. 108 The maximum prin iple (Theorem 4.4.1) would require us to maximize (4.87) pointwise over the ontrols, i.e. we are attempting to nd the Hamiltonian as, H = sup H v1,v2∈C∞(S1;R) ∫ 2π = sup (v2 (p1 cosx3 + p2 sin x3) + v1p3 v1,v ∈C∞2 (S1;R) 0 ) 1 ( ) χ ( )− v21 + v22 − v21,θ + v22 2 2,θ dθ (4.88) This maximization results in two Euler-Lagrange equations, ( ) ∂H − ∂ ∂H = 0, i = 1, 2, (4.89) ∂vi ∂θ ∂vi,θ that yields the optimal ontrols, Au1 = p3, (4.90) Au2 = p1 cosx3 + p2 sin x3, where we denote A = (1− χ∆), as usual. The Hamiltonian an be read as, ∫ H 1 2π = (u1Au1 + u2Au2) dθ 2 0 ∫ 1 2π ( ) = p3A−1p3 + (p1 cosx3 + p2 sin x −13)A (p1 cosx3 + p2 sin x3) dθ 2 0 (4.91) The dynami s of the ostate variable p an be al ulated from Hamilton's equa- ∂p = − δH δHtion, , where denotes the fun tional derivative of H with respe t to ∂t δx δx x. Expli it al ulations yield, ∂p1 = 0 ∂t ∂p2 = 0 (4.92) ∂t ∂p3 = −u2 (−p1 sin x3 + p2 cos x3) . ∂t 109 If we denote µ3 := −p1 sin x3 + p2 cosx3, then from (4.90) and (4.92), we an redis over (4.86). Renaming µi = Aui, i = 1, 2, (4.86) an also be written as, ∂µ1 = −µ A−13 µ2 ∂t ∂µ2 = µ A−1µ (4.93)3 1 ∂t ∂µ3 = −µ A−12 µ1. ∂t 4.6.2 Strong Coupling Limit, χ → ∞ Similar to the Heisenberg ase, it an be shown that syn hronization is a hieved in ∞ 1 strong oupling limit in this ase as well. We know that for some z ∈ C (S ;R), − ∫1 2π we have lim A z = zdθ. Thus, we may rewrite (4.93) in the strong oupling χ→∞ 0 limit as, ( ) ∫ ∂µ 2π1 = −µ3 µ2dθ ∂t 0 ( ) ∫ ∂µ 2π2 = µ3 µ1dθ (4.94) ∂t 0 ( ) ∫ ∂µ 2π3 = −µ2 µ1dθ . ∂t 0 If we dene the variables, ∫ 2π αi = µidθ, i = 1, 2, 3, (4.95) 0 then these variables will evolve a ording to, α̇1 = −α2α3 α̇2 = α α (4.96)1 3 α̇3 = −α1α2. 110 Equations (4.96) are the equations for a single agent s enario and they are studied in detail in Se tion 3.2. This indi ates the syn hronization phenomenon in the planar ontinuum o k. Note that α1 and α2 are essentially the optimal ontrols u1 and u2, respe tively. The Hamiltonian (4.88) simply be omes, 1 ( ) h∞ = α 2 1 + α 2 2 . (4.97)2 4.6.3 Simulation Results While it has not been possible to hara terize general solutions of (4.93) analyti- ally, here we demonstrate numeri al solutions. To numeri ally solve the evolution equations of µi variables, we used a nite dieren e method. We partitioning the spa e domain [0, 2π] uniformly in M points, 0 = θ1, . . . , θM = 2π, so that the 2π dieren e between two onse utive spa e points be ome δθ = . In this dis- M rete setting any z(t, θ) an be approximated as an M ve tor, z(t, θ) ≈ z(t) = [z1(t), z2(t), . . . , zM(t)] T with the onstraint z1(t) = zM(t) for all t to respe t the periodi ity property. Note also that in a se ond-order entral dieren e s heme, the double partial spa e derivative is expressed as, ∂2z(t, θj) zj+1 − 2zj + zj−1 = , ∂θ2 δθ2 for j = 1, . . . ,M with appropriate adjustments for the boundary points j = 1,M . The linear operatorA = (1−χ∆) an then be expressed by theM×M nonsingular 111 PSfrag repla ements t = 0 s t = 6.67 s t = 13.34 s t = 20 s 1 10 10 10 0.5 5 5 5 0 0 0 0 µ1µ2 µ3 -0.5 -5 -5 -5 -1 -10 -10 -10 0 5 0 5 0 5 0 5 θ θ θ θ 1.5 1 0.4 0.2 0.2 0 1 0.5 u1 0 -0.2 u2 0.5 0 -0.2 -0.4 0 -0.5 -0.4 -0.6 0 5 0 5 0 5 0 5 θ θ θ θ 0.01 0.5 1 2 0 0 0 0 -0.01 -0.5 -1 -2 0 10 5.5 6 6.5 6 7 8 4 6 8 x1 10 -3 x1 x1 x× 1 Figure 4.6: Numeri al solution of (4.99) (rst two rows) along with the state evolution (third row) for experiment 1. The arrows represent the dire tion of movement of the parti le at that point. matrix AM ,   ( ) ( ) ( ) 1 + 2χ − χ 0 . . . − χ 0  δθ2 δθ2 δθ2      ( ) ( ) ( )  − χ 1 + 2χ − χ 0 . . . 0   δθ2 δθ2 δθ2      ( ) ( ) ( ) AM =  0 − χ 1 + 2χ χ  . (4.98)  δθ2 δθ2 − 2 . . . 0δθ      . . . .   . . . . . . . .       ( ) ( ) ( ) 0 − χ2 . . . 0 − χ2 1 + 2χδθ δθ δθ2 112 x2 u µ , µ , µ1, u2 1 2 3 With this notation, the partial dierential equations (4.93) an be expressed as a system of ordinary dierential equations (ODE), µ̇ = −µ A−11 3 M µ2 µ̇ = µ A−1µ (4.99)2 3 M 1 µ̇3 = −µ −12AM µ1. These ODEs (4.99) are then solved using a mid-point based ODE solver in MAT- LAB. The optimal ontrols ui's an be derived from the µ variables by the relation, ui = A −1 M µi, i = 1, 2 whi h are used in the quadrature of the state variables xj 's. Here we present results of some experiments with varying initial onditions. The nal time T and spa e dis retization fa tor M is kept xed at T = 20 se onds and M = 128 for all the experiments. A high value of M is hosen for a faithful al ulation of the spatial derivatives. In the subsequent experiments we try to investigate the behavior of a simple loop under the optimal ontrols generated by (4.99), i.e. we take, x1(0, θ) = 0.01 cos(θ) and x2(0, θ) = 0.01 sin(θ) so that initially the parti les start on a ir le. The remaining initial onditions of x3 and µi variables and the parameter χ is varied in the following experiments. It is to be noted that the Hamiltonian and the Casimir variables are validated to be onstant in ea h of the experiments. Experiment 1 We take a simple example where ea h agent start moving in the positive x axis and with unit speeds, i.e. x3(0, θ) = 0, µ2(0, θ) = 1. The initial urvature of ea h 113 th th th th 0.02 10 10 0.01 10 0 0 0 0 PSfrag repla ements -0.01 -10 -0.02 -10 -10 0 -10 0 10 -10 0 10 -10 0 10 x1 x1 x1 x1 Figure 4.7: State evolution is only shown for experiment 2. The arrows represent the dire tion of movement of the parti le at that point. agent is taken to be zero, µ1(0, θ) = 0 whi h means every agent starts harmo- niously with same velo ity and urvature. The value of χ is taken as 1. Four snapshots of the µ, u and x variables are shown in Fig. 4.6. The urvature eld is seen to be forming two peaks in the spatial domain whi h gives rise to the twisted form of the initial ir le. Experiment 2 We keep all the initial onditions same as in experiment 1 ex ept the initial dire - tion of movement of the parti les. It is simulated that almost half the parti les try to go in one dire tion while the other half in the opposite dire tion. To write this as a ontinuous periodi fun tion, we take   ( ( ( ))) π   1 + tanh 100 θ − π , if 0 ≤ θ < π 2 2 x3(0, θ) = . (4.100)   ( ( ( ))) π  1− tanh 100 θ − 3π , if π ≤ θ ≤ 2π 2 2 This denition of the initial dire tion means that the parti les on the `east' part of the ir le initially go to the right and the parti les on the `west' part go to the left. The simulation results are shown in Fig. 4.7. Comparing to Fig. 4.6, 114 x2 th th th th 0.02 1010 10 0.01 5 0 0 0 0 PSfrag repla ements -0.01 -5 -0.02 -10 -10 -10 0 -10 0 10 -10 0 10 -5 0 5 x1 x1 x1 x1 Figure 4.8: State evolution is only shown for experiment 3. The arrows represent the dire tion of movement of the parti le at that point. the rst two rows are identi al sin e the initial onditions of µ variables did not hange and hen e they are omitted. What is interesting is that the ir le splits into two loops onne ted by very small number of parti les. Experiment 3 Similar to experiment 2, we try to investigate the ee t of hange in initial dire - tion of movement of the parti les. Here, we set the parti les to go on a radially outward path, i.e. x3(0, θ) = θ, with keeping all other onditions same as in ex- periment 2. The simulation results are shown in Fig. 4.8. The rst two rows are not shown sin e they are identi al with experiment 1. Experiment 4 In this experiment, again we x all the initial onditions and parameters same as in experiment 1, ex ept the initial urvature is given a lo al intensity. In other words we hoose this Gaussian fun tion, 1 (θ 2−π) µ1(0, θ) = √ e− 2σ2 , (4.101) 2πσ2 115 x2 PSfrag repla ements t = 0 s t = 6.67 s t = 13.34 s t = 20 s 10 10 20 10 5 5 5 10 0 µ1µ2 0 0 0 µ3 -5 -5 -5 -10 -10 0 5 0 5 0 5 0 5 θ θ θ θ 1.1 1.5 1.5 1.5 1 1 1 1.05 u1 0.5 0.5 0.5 u2 1 0 0 0 0.95 -0.5 -0.5 -0.5 0 5 0 5 0 5 0 5 θ θ θ θ 2 1 0.01 1 1 0 0.5 0 0 -0.01 0 -1 -1 0 10 3 3.5 4 5 6 4 4.5 5 5.5 x1 -3 x x×10 1 1 x1 Figure 4.9: Numeri al solution of (4.99) (rst two rows) along with the state evolution (third row) for experiment 4. The arrows represent the dire tion of movement of the parti le at that point. with σ = 0.05. The purpose of hoosing this initial ondition is to see whether a lo alized information gets spread a ross the ontinuum or not. The simulation results are shown in Fig. 4.9. Experiment 5 Finally, we demonstrate the ee t of strong oupling limit by taking a large value 116 x2 u1, u2 µ1, µ2, µ3 PSfrag repla ements t = 0 s t = 6.67 s t = 13.34 s t = 20 s 2 10 20 20 5 10 10 1 0 0 0 µ1µ2 0 µ3 -5 -10 -10 -1 -10 -20 -20 0 5 0 5 0 5 0 5 θ θ θ θ 1.5 1 0.75 0.5 1 0 u1 0.5 0.5 0.7 u2 -0.5 0 -0.5 0 0.65 -1 0 5 0 5 0 5 0 5 θ θ θ θ 0.02 1.67 3.9 0.01 4.491.66 3.89 4.48 0 1.65 3.88 4.47 -0.01 1.64 3.87 4.46 0 10 6.46 6.47 6.48 12.5 12.515 17.4 17.42 x1 10 -3 x1 x x× 1 1 Figure 4.10: Numeri al solution of (4.99) (rst two rows) along with the state evolution (third row) for experiment 5. The arrows represent the dire tion of movement of the parti le at that point. The u1, u2 solutions are almost `at', indi ating single agent solution or syn hronization. of χ = 1000. We note that even in the ase, µ1(0, θ) = cos(θ), µ2(0, θ) = 1 + 0.2 cos(θ), µ3(0, θ) = sin(θ), x3(0, θ) = π/4, we essentially get the system derived by the optimal ontrols that are spatially non varying. 117 x2 u1, u2 µ1, µ2, µ3 4.7 Dis ussion and S ope of Future Resear h In this hapter, we have presented a general framework for a ontinuum des rip- tion of a o k. We are interested in solving optimal ontrol problems to explain olle tive movement of su h a o k. We re ognize this is a hallenging problem that naturally provides several open questions for further resear h. We itemize few su h possibilities. • It is shown that under a ertain nite odimensionality of a rea hable set, the PMP remains valid in a general Hilbert spa e setting. One might want to dis over its relationship with the ontrollability ondition. In parti ular, if a system on the loop group is strongly ontrollable, does the PMP ondition satisfy automati ally? We have been able to show this to be true in the Heisenberg ase. Does this remain valid if we only have weak ontrollability? • It is of interest to extra t meaningful features of the optimal solutions of the SE(2) ase. While we have not been able to solve (4.93) analyti ally, we want to answer few questions about it. For example, do these equations possess a traveling wave like solution (like the Heisenberg ase)? If so, what is the speed of those waves? The answer might give an insight toward the information transfer in biologi al swarms. It is the inherent nature of the numeri al study presented here that there exist many possibilities by varying initial ongurations whi h makes it parti ularly di ult. We explored only a tiny fra tion of possible variations in the initial ongurations in this 118 do ument. A future work ould perform more extensive numeri al study of these partial dierential equations (4.93). • We have presented the results in this hapter under the ase of a xed y li intera tion topology. A more general, possibly state dependent (hen e time dependent, too) intera tion s heme an be modeled and subje ted under similar questions. • It will be an interesting future work to establish ontinuum parallel of the Lie-Poisson redu tion in the loop group ase. 119 Chapter 5 Cognitive Cost of Flo king: A Geometri and Hamiltonian Perspe tive 5.1 Introdu tion It has been an appealing question to resear hers from several elds to address how natural olle tives fun tion at a fundamental level. Many theories have been proposed to des ribe this phenomenon over the past few de ades. The la k of a urate motion apturing te hnologies had limited the study of natural olle tives for many years. However, as motion apture be ame more sophisti ated, more movement data of these olle tives were re orded. This enabled resear hers to un over several underlying me hanisms behind o king [Ballerini et al., 2008a,b; Attanasi et al., 2014; Cavagna et al., 2018; Nagy et al., 2010℄. These studies shed light on how individual agents intera t with its neighboring agents or how information may be propagated through the o k. 120 Continuing the spirit of the `top-down' view of the o k, we will present a novel perspe tive for analyzing olle tive motion data. The o k movement results in a time-series of its kineti energy, whi h an be divided into several energy modes. Normalized modes dene a urve in some appropriate dimensional simplex whi h we attempt to des ribe by an evolutionary game dynami s. Individual energy modes are onsidered as pure strategies of su h a game. An optimal ontrol problem is proposed to best t the data on the simplex, where the ontrol inputs modulate the tness asso iated with the strategies. This is in ontrast to the optimal ontrol problem posed in Chapter 4, where the ontrols are `low-level' i.e. individual ontrol inputs are determined post-solution of the optimal ontrol problem. In the present ontext however, the ontrol inputs are `high-level'. The olle tive itself is thought to be de iding the optimal allo ation of its energy among several dierent modes during a ight event. A notion of ognitive ost is introdu ed to denote the optimal ost for the olle tive to perform this allo ation. This work brings together several key ingredients for this data-driven approa h. In se tion 5.2, the motion data of European o ks is detailed. This data is then subje ted to a linear data smoothing te hnique [Dey and Krishnaprasad, 2012℄ that re onstru ts smooth traje tory data of ea h bird in the o k. A nonlinear data smoothing te hnique [Dey and Krishnaprasad, 2014℄ is later used for the optimal energy allo ation problems. These smoothing te hniques are based on optimal ontrol theory and are des ribed in se tion 5.3. Se tion 5.4 ontains the geometri theory developed in [Mis hiati and Krishnaprasad, 2017℄ to ompute 121 dierent energy modes. Finally, in se tion 5.5, a generative model based on [Raju and Krishnaprasad, 2018℄ is des ribed to onstru t an optimal ontrol problem on a simplex. Numeri al solution of this optimal ontrol problem, as well as the idea of ognitive ost, are presented in se tion 5.6. This is a joint work with V. Raju [Halder et al., 2019b℄. 5.2 Flo king Data We are provided with ight data of European starlings that were taken by Dr. Andrea Cavagna and his ollaborators from the Colle tive Behaviour in Biologi al Systems (COBBS) group at the Institute for Complex Systems (ISC-CNR), Uni- versity of Rome La Sapienza". Starlings gather around urban areas during the winter months in order to get extra warmth from the ities. Flo ks of these kind of birds are well known to perform remarkable maneuvers, the purpose and me ha- nisms of whi h still elude resear hers. Equipped with modern imaging te hniques and sophisti ated algorithms for stereo re onstru tion, these group of resear hers managed to apture a series of ight events with dierent o k sizes in the winter months of 2011. See [Attanasi et al., 2014℄ for more details about the pro ess. We will study eight parti ular o king events, the details of whi h are given in table 5.1. 122 Flo king Flo k Size Duration Data Capture Rate Event (n) (se onds) (frames/se ond) 1 175 5.4875 80 2 123 1.8176 170 3 46 5.6118 170 4 485 2.3471 170 5 104 3.8824 170 6 122 4.1588 170 7 380 5.7353 170 8 194 1.7588 170 Table 5.1: Details of aptured o king events 5.3 Data Smoothing Given a time-indexed sequen e of sampled observations on a manifold, genera- tive models provide a meaningful way of apturing them through the use of an underlying dynami al system omplete with ontrol inputs having useful interpre- tations. The ontrol inputs are determined by solving an optimal ontrol problem, where the ost fun tion onsists of a tness term that penalizes mismat h between the generated traje tory and sampled data, and a smoothing term weighted by a parameter λ that ae ts the smoothness of the generated traje tory. We dis uss two generative models that have been proposed to solve this problem. 123 5.3.1 A linear generative model A rst approa h to solving the data smoothing problem, presented in [Dey and Krishnaprasad, 2012℄, is to formulate an optimal ontrol problem to minimize the jerk path integral, with intermediary state osts determining the t error. Suppose N 3 that {ri}i=0 denote the positions of the birds at ea h sampling time, with ri ∈ R . 3 In order to re over a traje tory t r(t) : [t0, tN ] → R , one an use the jerk-driven linear generative model, ṙ(t) = v(t) v̇(t) = a(t) ȧ(t) = u(t) (5.1) where v(t), a(t),u(t) denote the velo ity, a eleration and jerk (input) of the tra- je tory. The ost fun tional to be minimized is N tN∫ ∑ J = || 2 2l r(ti)− r(t)|| + λ ||u(t)|| dt (5.2) i=0 t0 where the minimization is over initial onditions r(t0),v(t0), a(t0) and the input u(t). Dening the state and output as   r(t)       x(t) = 9 3 (t)  ∈ R ,y(t) = x(t) ∈ Rv       a(t) 124 we obtain the linear state equations ẋ(t) = Ax(t) +Bu(t) y(t) = Cx(t), where     0 13 0 0            A =  0 0 1 , B =  0  , C = [13 0 0] (5.3)  3            0 0 0 13 Therefore, the problem of minimizing Jl subje t to (5.3) is a linear, quadrati optimal ontrol problem, whi h an be solved by a ompletion of squares of terms in the ost by invoking a path independen e lemma, or by applying the Pontryagin Maximum Prin iple as shown in [Dey and Krishnaprasad, 2012℄. This approa h has been used to smooth the starling o k data for all the events listed in table 5.1. 5.3.2 Data smoothing in the Eu lidean setting In this se tion, we present a general result on the Pontyagin Maximum Prin iple n based approa h for data smoothing on the Eu lidean spa e R . Suppose that { } xd N i denote the sampled data. For a generative model given by the dynami si=0 ẋ = f(x, u) n mon R , with the ontrol u ∈ R , the optimal ontrol problem an be formulated as: ∫ N λ tN ∑ min J(x(t0), u) = ‖u‖2 dt+ Fi(x(ti), xdi ), x(t0), u∈Rm 2 t0 i=0 (5.4) subje t to: ẋ = f(x, u), 125 where parameter λ > 0 is a regularization parameter, and Fi's are suitably dened t errors of the re onstru ted traje tories and sampled data at the sampling times. Using Pontryagin's Maximum Prin iple, the optimal ontrol values an be al ulated as a fun tion of the state and a o-state variable. The following result from [Dey and Krishnaprasad, 2014℄ states this pre isely. Theorem 5.3.1. (PMP for data smoothing [Dey and Krishnaprasad, 2014℄ ) Let u∗(·) ∗be an optimal ontrol input for (5.28), and let x (·) denote the orresponding n state traje tory. Then there exist a ostate traje tory p : [t0, tN ] → R , p 6= 0 , su h that ∗ ∂Hẋ = (t, x∗, p, u∗) ∂p (5.5) ∂H ṗ = − (t, x∗, p, u∗) ∂x during t ∈ (ti, ti+1), i = 0, 1, ..., N − 1, and the Hamiltonian is given as H(t, x∗, p, u∗) = max H(t, x∗, p, v), (5.6) v∈Rm for t ∈ [t0, tN ]\{t0, t1, ..., tN}, where the pre-Hamiltonian is dened asH(t, x, p, u) = 〈p, f(x, u)〉 − λ ‖u‖2. Moreover, jump dis ontinuities of the ostate variable an 2 be written as p(t−0 ) = 0, p(t+)− p(t− ∂Fi(x(ti))i i ) = , i = 0, 1, ..., N, (5.7)∂x(ti) p(t+N) = 0. The pie ewise ontinuous nature of the o-state traje tory due to jump on- ditions arising from mismat h between the sampled data points and the re on- 126 stru ted state must be noted here. The initial ondition x(t0) is identied by using the terminal ondition for the o-state, while the optimal value of λ is typi ally obtained through leave-one-out or ordinary ross validation. The re onstru ted traje tory is then obtained as the proje tion onto the state spa e of the solution of Hamilton's equations derived from the (maximized pre-) Hamiltonian. We refer the reader to [Dey, 2015℄ for a detailed treatment of these problems. This is the result that we will use in our data tting problem on a simplex. 5.4 Energy Modes Avian o ks display a variety of ight behaviors that may be hara terized as olle tive strategies su h as steady dire ted translation of enter of mass (whi h we denote by om), oherent rotation about enter of mass (rot), hange of form (ens), internal re-shuing of relative positions (dem), rapid expansion or on- tra tion of volume (vol) et . A o king event may display all of the mentioned strategies to varying degrees as governed by the time-dependent allo ation of ki- neti energy to ea h strategy. We take the viewpoint presented in [Mis hiati and Krishnaprasad, 2017℄ and study the fra tions of the total kineti energy of a o k allo ated to several `kinemati modes'  rigid translations, rigid rotations, inertia tensor transformations, expansion and ompression, in order to des ribe olle tive behavior. If the positions of the birds in a o k are denoted by {r1, r2, ..., rn}, the enter 127 of mass an be written as, n 1 ∑ r = ri, (5.8) om n i=1 where we treat every bird alike, i.e. their masses are taken to be equal. The ensemble inertia tensor is dened by n ∑ K = (ri − Tr ) (r − r ) . (5.9) om i om i=1 Let the velo ities of the birds be denoted as, {vr1, ...,vrn}, then the total kineti energy is, n 1 ∑ E = ‖ 2vri‖ . (5.10) 2 i=1 We an dene the position and velo ity ve tor with respe t to the enter of mass, i.e. c , [c1, ..., c ] ∈ R3×n 3×nn , where ci = ri − r ; vc , [vc1,vc2, ...,v ] ∈ R , om cn where vci = vri − v . Then, om n n 1 ∑ E = ‖ 2v ‖ , E ‖ 2, vci‖ . (5.11) om om rel 2 2 i=1 We thus have the splitting, E = E + E . As presented in [Mis hiati and om rel Krishnaprasad, 2017℄, instantaneous relative energy allo ations an be expressed 4 1 on a probability simplex (∆ ) by exploiting the ber bundle stru tures of the o k's total onguration spa e to split the total kineti energy using (i) ensemble bration or (ii) shape bration. (i) Ensemble Fibration: We note that the ensemble inertia tensor K (5.9) is a symmetri positive denite matrix. Hen e its eigende omposition an be 1 n Note that in this hapter we will use ∆ to denote the n-dimensional simplex. 128 written as, K = QΛQT, with Λ = diag(λ1, λ2, λ3), where λ1 ≥ λ2 ≥ λ3 > 0. Dene, F := cvT+v cTc and F̃ = [F̃ T ij ] = Q FQ. Then the following energyc modes an be al ulated, ( ) 1 F̃ 2 F̃ 2 F̃ 2 E , 12 + 13 + 23 ens.rot 2 λ1 + λ2 λ1 + λ3 λ2 + λ3 ( ) (5.12) 1 F̃ 211 F̃ 2 2 E , + 22 F̃ + 33 . ens.def 8 λ1 λ2 λ3 Furthermore, ( ) 1 2 Ttr cv E c, , (5.13) vol 2 tr(K) so that, E = E − E . We may also al ulate E = E − ens.res ens.def vol dem rel E − E . Hen e, in this bration we have the following splitting of ens.rot ens.def the kineti energy, { } E E E E E om dem ens.rot vol ens.res , , , , ∈ ∆4 (5.14) E E E E E (i) Shape Fibration: Dene n ∑ J = (ci × vci) , i=1 (5.15) n ∑ ( ) Ic = ‖ci‖2 T1− cici . i=1 Then the rotational energy E an then be al ulated as, rot 1 E , JTI−1J, (5.16) rot 2 c The shape residual energy is given by E = E −E −E , whi h shp.res rel rot end.def provides the splitting in this bration as below { } E E E E E om rot shp.res vol ens.res , , , , ∈ ∆4 (5.17) E E E E E 129 While we an split the kineti energy in 5 dierent modes (5.14),(5.17), many o king events show a predominant allo ation of nearly onstant energy of rigid translation (E ). We ex lude this omponent from the total E in our analysis, om and onsider the allo ation of the remaining energy E to obtain a time dependent rel tra e of ea h event on a lower dimensional simplex. In parti ular, we apture the tra e generated by the following de omposition of E using ensemble bration rel on the 1-simplex by two dierent methods, { } E E dem ens ( 1ENS-I) , ∈ ∆ , (5.18) E E rel rel { } E E − E ens.rot rel ens.rot (ENS-II) , ∈ ∆1, (5.19) E E rel rel where E = E−E , and E = E −E = E +E +E . Similarly, rel om ens rel dem ens.rot vol ens.res a one dimensional simplex des ription using shape bration may be given by two ways, { E E − }E shp.res rel shp.res ( 1SHP-I) , ∈ ∆ , (5.20) E E rel rel { } E E rot shp (SHP-II) , ∈ ∆1, (5.21) E E rel rel where E = E − E = E + E + E . shp rel rot shp.res vol ens.res In this way, moment-to-moment de isions made by individuals in a o k, tak- ing a ount of the de isions of their neighbors, ontribute to o k-s ale strategies as aptured by su h time dependent tra es on the probability simplex. Treating the strategy prevalen e as being given by the respe tive energy fra tions, we resort to a generative evolutionary game dynami s to model the ompetition between the o k-s ale strategies. 130 5.5 Generative model on the 1-simplex and the data-smoothing problem Sin e we are interested in des ribing the evolution of two o k strategies as in eqs. (5.18) and (5.19) for ensemble bration or eqs. (5.20) and (5.21) for shape bration, we apture the tra e of o king events via a generative model on the 1-simplex. We onsider an evolutionary game model, namely repli ator dynami s equipped with a multipli ative ontrol, in order to des ribe their evolution in the interior (0, 1) of the one-dimensional simplex. The hoi e of repli ator dynami s is inuen ed by its universality in des ribing simplex-preserving dynami s, and by virtue of being an extremal for a variational problem [Svirezhev, 1972; Raju and Krishnaprasad, 2018℄. Presently, with the in lusion of a ontrol variable, we onsider a dierent variational problem that aims to perform data smoothing using regularization as in [Dey and Krishnaprasad, 2014℄. To see this, let x = [x x ]T ∈ ∆11 2 where xi, i = 1, 2 denote the prevalen e of strategies i (to be spe ied) on the simplex with the natural onstraint x1 + x2 = 1. xi = 1, i = 1, 2 orrespond to allo ation of Erel entirely to one of the two pure strategies. Suppose that the frequen ies asso iated with the strategies are updated a ording to the rule f i(x) xi(t+ 1) = xi(t) (5.22) f̄(x) 131 i where the tness f (x) = Ax and f̄ = x1f 1(x) + x f 22 (x) 2 . Here, A = [aij ] ∈ R th th denes a payo matrix with aij denoting the payo of the i strategy against j strategy. In the ase that the payos do not depend on the strategy j of against whi h it is mat hed up, the olumns of A are identi al. In the ode limit of (5.22), after an inhomogeneous time-s ale hange, we get the mean eld equations: ẋi(t) = xi(t)(f i(x)− f̄(x)), i = 1, 2 (5.23) It an be readily veried that (5.23) is simplex-preserving, leaving the pure strate- gies invariant. Sin e addition of the same term to ea h omponent of the tness keeps the dynami s (5.23) un hanged, by subtra ting a21 and a12 from the rst and se ond olumn elements of A respe tively, we get the equivalent payo matrix   a  11 − a21 0  Ã =   (5.24)   0 a22 − a12 We introdu e a ontrol input ũ that s ales the tness, and hoose the parameters of the matrix su h that a11 − a21 = −(a22 − a12) = 1 so that the tness an be rewritten as:   1 0   f(x) = ũ  x (5.25)   0 −1 Due to the simplex onstraint, (5.23) is ompletely des ribed using x = x1: ẋ(t) = ũ(t)x(t)(1 − x(t))(f 1(x)− f 2(x)) (5.26) with x = 0, 1 orresponding to the pure strategies 2 and 1 respe tively. Due to 1 2 our hoi e of the payo matrix parameters, f − f is a onstant. This allows 132 1 2 us to adopt a time-s ale hange by the fa tor f − f to arrive at our generative model: ẋ(t) = u(t)x(t)(1− x(t)) (5.27) This dynami s results in asymptoti onvergen e to the pure strategy x = 1 in the absen e of ontrol, that is, when u(t) ≡ 1. However, the time-varying on- trol variable u serves to model hanging preferen es for the o k strategies by appropriate hanges in its sign and magnitude. Su h a temporal modulation of the tness ensures feasibility of apturing arbitrary tra es in the interior of the simplex. {xd d d dGiven a set of data points 0, x1, ..., xN} with ea h xk ∈ (0, 1), k = 0, 1, ..., N , at time instants {t0, t1, ..., tN}, we formulate the optimal ontrol problem, ∫ N λ tN ∑ min J(x(t0), u) = u 2dt+ Fi(x(ti)), x(t0), u∈R 2 t0 i=0 (5.28) subje t to: ẋ = ux(1− x), where the t errors Fi's are given by the Kullba k-Leibler divergen e measure of mismat h between the data and the state, ( ) ( ) xd 1− xd Fi(x) = x d i i log + (1− xd) log ii − , i = 0, 1, ..., N. (5.29)x 1 x We an dire tly appeal to Pontryagin's Maximum Prin iple (PMP) and theo- rem (5.3.1) to write ne essary onditions for optimality. We an write the pre- Hamiltonian as, λ H(x, p, u) = upx(1− x)− u2. (5.30) 2 133 The Hamiltonian maximization ondition (5.6) yields an optimal ontrol in ea h time interval t ∈ (ti, tt+1), i = 0, 1, ..., N − 1, 1 u = px(1− x), (5.31) λ with Hamiltonian given by, H 1(x, p) = p2x2(1− x)2. (5.32) 2λ Hamilton's equations (5.5) read, 1 ẋ = px2(1− x)2 λ (5.33) 1 ṗ = − p2x(1− x)(1− 2x). λ The jump onditions for p (5.7) an be written as, p(t−0 ) = 0, x(t di)− x p(t+i )− p(t−) = i , i = 0, 1, ..., N, (5.34)i x(ti)(1− x(ti)) p(t+N) = 0. du Remark 5.1. Note that the optimal ontrol is pie ewise onstant sin e = 0 dt for ea h of these time intervals t ∈ (ti, ti+1), i = 0, 1, ..., N − 1. Therefore, denoting xk = x(tk), k = 0, 1, ..., N , any optimal ontrol an be des ribed by a ve tor (u0, u1, ..., uN) with the onditions 1 u0 = (x0 − xd λ 0 ), 1 uk − uk−1 = (x dk − xk), k = 1, 2, ..., N (5.35)λ uN = 0. 134 Pie ewise onstan y of the ontrol input allows us to write the solution to the state equation (5.26) expli itly. Suppose the sampling time of the tra e is uniform, i.e. ∆t := tk+1 − tk, ∀k ∈ {0, ..., N − 1}, integrating the state equation (5.26) in (tk, tk+1), we an write x euk∆tk xk+1 = − , k = 0, 1, ..., N − 1. (5.36)1 + xk (euk∆t 1) By iteration, we an in turn write every xk as a fun tion of x0 and u0, u1, ..., uk−1, x0e (u0+u1+···+uk )∆t−1 xk = xk(x0) = , k = 1, 2, ..., N. (5.37) 1 + x (e(u0+u1+···+uk 1)∆t0 − − 1) The endpoint ondition (uN = 0) an then be written as, x + x + · · ·+ x = xd0 1 N 0 + xd1 + · · ·+ xdN , (5.38) where the left hand side of (5.38) is a fun tion of x0. Solving the optimal ontrol problem (5.28) thus boils down to solving (5.38) for x0 ∈ (0, 1). Remark 5.2. The value of the regularization parameter λ is usually hosen through ross validation te hnique. We do not employ any su h te hniques here. The value of λ is hosen su h that the root nding algorithm for solving (5.38) onverges for all events. For λ = 0.2, the roots were found with reasonably good −5 a ura y with value of the fun tion at the root being of the order of 1 × 10 or lower for all events. For lower λ however, the problem be omes stier and left hand side of (5.38) demonstrates `ee tive dis ontinuity' in x0. This poses seri- ous problem in solving (5.38). It is to be noted that the original aptured ight data was subje ted to data-smoothing to obtain smooth traje tories [Dey, 2015℄. 135 0.5 1 PSfrag repla ements PSfrag repla ements ENS-I ENS-I SHP-I 0.4 0.8 SHP-I ENS-II ENS-II SHP-II SHP-II 0.3 0.6 0.2 0.4 0.1 0.2 0 0 0 2 4 6 8 0 2 4 6 8 Event Event (a) Time Averaged Hamiltonian Integral (b) Time Averaged Total Cost Figure 5.1: Hamiltonian signatures The data-smoothing problem in [Dey, 2015℄ onsidered a linear generative model as in Se tion 5.3.1 and used ordinary ross validation for traje tory of ea h bird to determine the appropriate weight to the regularization term. This generated smooth traje tories with suppressed level of noise ompared to the original data. d d We then take the sampled data {x0, · · · , xN} from these smooth traje tories. This an justify taking same value of λ a ross all the events. As a future step, ross validation ould be employed to arrive at a good value of λ in the range where (5.38) an be solved. 5.6 Data Fitting Results For all 8 events, we solve the optimal ontrol problem (5.28) and report the results here. The value of the regularization weight λ is taken to be 0.2 and 100 data samples at regular time intervals are taken for all events. Given the data 136 ∫ 1 Hdt T J/T ve tor, we solve equation (5.38) for x0 ∈ (0, 1). In Table 5.2, we report time averaged Hamiltonian integrals and time averaged total osts for all the dierent games that we onsider in eqs. (5.18) to (5.21). These time averaged Hamiltonian integrals are thought of as ognitive osts of the events. As seen from Table 5.2, the trend of (ENS-I) losely follow the game (SHP-I), while the other two games seem to follow ea h other. This is graphi ally represented in Fig. 5.1. Optimal ontrol solutions for the games ENS-I (5.18) and SHP-II (5.21) for individual events are shown in Fig. 5.25.9. We note that more variation in the energy time signal results in higher ognitive ost (in both measures). This is interpreted as the olle ting having to `think' more to properly allo ate the modes, in urring higher osts. These ognitive osts for a parti ular game an thus indi ate overall physi al behavior of the o k. For example, in the games (ENS-II) or (SHP-II) where a rotational energy is onsidered as one of the pure strategies, relatively higher ognitive osts for event 2, 5 indi ate that the o ks went through more rotations than the other events during the ight periods. On the other hand, low ost for event 4 is justied by almost re tilinear overall motion. Similar on lusions an be drawn for the other set of games (ENS-I) and (SHP-II), where the respe tive ognitive ost will stipulate nature of variation of the demo rati (reshuing within the o k) energy. The higher the ost is, more aggressively the relative positions of the birds within the o ks are hanged, leading to a more omplex ight event. 137 ∫ ∫Hdt J(∫x0,u) Duration dt dt (se onds) (ENS-I) (SHP-I) (ENS-II) (SHP-II) (ENS-I) (SHP-I) (ENS-II) (SHP-II) 5.4875 0.1232 0.1263 0.0976 0.1077 0.1981 0.1975 0.1454 0.1499 1.8176 0.1432 0.1018 0.2210 0.1760 0.2227 0.1619 0.3769 0.3118 5.6118 0.2735 0.2392 0.0613 0.1073 0.4595 0.4092 0.1557 0.2495 2.3471 0.1021 0.1270 0.0107 0.0190 0.2440 0.2702 0.0594 0.0610 3.8824 0.0779 0.2699 0.1587 0.1383 0.0896 0.3692 0.3001 0.3041 4.1588 0.1809 0.1634 0.0846 0.1105 0.2799 0.2706 0.2063 0.2090 5.7353 0.0804 0.1293 0.0576 0.0619 0.1127 0.2079 0.1087 0.1221 1.7588 0.4569 0.4069 0.0731 0.1090 0.8037 0.8361 0.2074 0.3810 Table 5.2: Hamiltonian Signature 5.7 Dis ussion In this hapter, we have brought together several results from geometry, optimal data-tting and evolutionary game theory to asso iate a ognitive aspe t to o k- ing. The ight data of Starling o ks give rise to time-signals of energy mode distributions. Here, the whole o k is on eptualized to making de isions about how to optimally allo ate its energy in several modes. The dierent energy modes are thought as pure strategies of an evolutionary game and their tness is mod- ulated by some de ision or ontrol variables. These ontrols are then determined by optimally tting this model to the observed energy mode distributions in the data. The ost to this data-tting are referred to as ognitive ost for the o k. In this work, we have only onsidered splitting energy into two modes. In this 138 setting, the optimal ontrol solutions present interesting hara teristi s. It will be an important dire tion to onsider energy splitting in several energy modes, hen e solving the tting problem in a higher dimensional simplex. It will also be of interest to interpret the ognitive osts in su h s enarios. 139 (a) Flo k Traje tory 1 3 0.9 2.5 0.8 2 0.7 1.5 0.6 1 0.5 0.5 0.4 0 Trajectory fit 0.3 Data -0.5 0.2 -1 0.1 -1.5 0 -2 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 2 Trajectory fit 0.9 Data 1.5 0.8 1 0.7 0.5 0.6 0.5 0 0.4 -0.5 0.3 -1 0.2 -1.5 0.1 0 -2 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.2: Event 1, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel E rot (d)-(e) x = (SHP-II) E rel 140 x, x x, x d d u u (a) Flo k Traje tory 1 1.5 0.9 1 0.8 0.5 0.7 0 0.6 -0.5 0.5 -1 0.4 -1.5 0.3 0.2 -2 0.1 -2.5Trajectory fit Data 0 -3 0 0.5 1 1.5 2 0 0.5 1 1.5 2 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 2 Trajectory fit 0.9 Data 1.5 0.8 1 0.7 0.5 0.6 0 0.5 -0.5 0.4 -1 0.3 0.2 -1.5 0.1 -2 0 -2.5 0 0.5 1 1.5 2 0 0.5 1 1.5 2 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.3: Event 2, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel E rot (d)-(e) x = (SHP-II) E rel 141 x, x x, x d d u u (a) Flo k Traje tory 1 3 0.9 2 0.8 1 0.7 0 0.6 Trajectory fit Data 0.5 -1 0.4 -2 0.3 -3 0.2 -4 0.1 0 -5 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 2.5 Trajectory fit 0.9 Data 2 0.8 1.5 0.7 1 0.6 0.5 0.5 0 0.4 -0.5 0.3 -1 0.2 -1.5 0.1 -2 0 -2.5 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.4: Event 3, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel E rot (d)-(e) x = (SHP-II) E rel 142 x, x x, x d d u u (a) Flo k Traje tory 1 2.5 0.9 2 0.8 1.5 0.7 1 0.6 0.5 0.5 0.4 0 0.3 -0.5 0.2 -1 0.1 Trajectory fit Data 0 -1.5 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 0.8 Trajectory fit 0.9 Data 0.6 0.8 0.4 0.7 0.2 0.6 0.5 0 0.4 -0.2 0.3 -0.4 0.2 -0.6 0.1 0 -0.8 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.5: Event 4, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel (d)-(e) x = Erot (SHP-II) E rel 143 x, x x, x d d u u (a) Flo k Traje tory 1 1.5 Trajectory fit 0.9 Data 1 0.8 0.7 0.5 0.6 0 0.5 -0.5 0.4 0.3 -1 0.2 -1.5 0.1 0 -2 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 2 Trajectory fit 0.9 Data 1.5 0.8 1 0.7 0.5 0.6 0 0.5 -0.5 0.4 -1 0.3 0.2 -1.5 0.1 -2 0 -2.5 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.6: Event 5, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel E rot (d)-(e) x = (SHP-II) E rel 144 x, x x, x d d u u (a) Flo k Traje tory 1 3 0.9 2 0.8 0.7 1 0.6 0.5 0 0.4 -1 0.3 0.2 -2 0.1 Trajectory fit Data 0 -3 0 1 2 3 4 5 0 1 2 3 4 5 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 1 Trajectory fit 0.9 Data 0.5 0.8 0 0.7 -0.5 0.6 0.5 -1 0.4 -1.5 0.3 -2 0.2 -2.5 0.1 0 -3 0 1 2 3 4 5 0 1 2 3 4 5 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.7: Event 6, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel x = Erot(d)-(e) (SHP-II) E rel 145 x, x x, x d d u u (a) Flo k Traje tory 1 1.5 0.9 1 0.8 0.5 0.7 0 0.6 0.5 -0.5 0.4 -1 0.3 -1.5 0.2 -2 0.1 Trajectory fit Data 0 -2.5 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 2 Trajectory fit 0.9 Data 1.5 0.8 0.7 1 0.6 0.5 0.5 0 0.4 0.3 -0.5 0.2 -1 0.1 0 -1.5 0 1 2 3 4 5 6 0 1 2 3 4 5 6 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.8: Event 7, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel E rot (d)-(e) x = (SHP-II) E rel 146 x, x x, x d d u u (a) Flo k Traje tory 1 5 0.9 4 0.8 3 0.7 2 0.6 0.5 1 0.4 0 0.3 -1 0.2 -2 0.1 Trajectory fit Data 0 -3 0 0.5 1 1.5 2 0 0.5 1 1.5 2 Time (s) Time (s) (b) Traje tory Fit ( ) Optimal Control 1 2 Trajectory fit 0.9 Data 1.5 0.8 0.7 1 0.6 0.5 0.5 0.4 0 0.3 0.2 -0.5 0.1 0 -1 0 0.5 1 1.5 2 0 0.5 1 1.5 2 Time (s) Time (s) (d) Traje tory Fit (e) Optimal Control E dem Figure 5.9: Event 8, λ = 0.2, Number of samples = 100, (b)-( ) x = (ENS-I), E rel E rot (d)-(e) x = (SHP-II) E rel 147 x, x x, x d d u u Chapter 6 Con lusions and Dire tions for Future Re- sear h We have made an attempt to explain olle tive behavior in natural o ks in this thesis. Possible roboti appli ations in this ontext are also presented. The thesis is distin tively divided into two parts depending on the underlying approa h  either olle tive behavior is viewed as an emergen e of intera tions between small number of agents in a `bottom-up' fashion or those intera tions are inferred in a `top-down' way. We summarize below the ontributions of this dissertation along with dire tions in whi h this line of resear h an be ontinued. In Chapter 2, we explored inter-agent intera tion strategies from both theoret- i al and implementation perspe tives. First, we onsider a two-agent s enario in whi h one agent pursues the other using onstant-bearing (CB) pursuit law. The pursued agent behaves like a moving bea on whose movement is independent to the other. It is then shown that under parti ular parameter setting of the CB law 148 and onstant urvature paths of the bea on, the ombined equations of motion of the system resemble to that of a ouple of gravitating parti les. Periodi orbits are shown to exist, ea h orresponding to a xed energy. This result is used in a roboti appli ation subsequently. We have ee tively utilized the results of this problem in the problem of en ir ling a stati bea on that is sensed visually by a mobile robot by means of a amera with limited eld-of-view (FOV). Proper feed- ba k law for the robot is advised to make a desired losed loop in the phase spa e asymptoti ally stable. This guarantees intermittent appearan e of the bea on in the amera's FOV. Laboratory demonstration of this problem in orporates online estimation of the bea on's position when it falls out of the FOV. Se ondly, labo- ratory implementations of two biologi ally plausible feedba k laws are presented. These laws in lude another dual-agent law alled Mutual Motion Camouage and a multi-agent onsensus type law alled Topologi al Velo ity Alignment. In this hapter we have shown how omplex olle tive motion patterns an emerge from simple intera tions among the agents in a o k. We study the problem of optimal steering of a single agent in Chapter 3. The agent is driven from initial to nal onguration on the plane while minimizing the ontrol ost that penalizes both speed and urvature ontrol. Optimal ontrol solution is obtained by using Pontryagin's Maximum Prin iple (PMP) and Lie- Poisson redu tion te hnique. Optimal traje tories are ategorized by the values of the Hamiltonian and another onserved quantity alled Casimir. This problem is then extended to apture the s enario of a o k of agents moving on the plane. 149 These agents intera t with ea h other by a predened graph. The individual ontrol osts are oupled with mismat h in ontrol with the neighbors. This hapter forms a bridge between the two parts of the thesis. In Chapter 4, we take the problem of optimal steering of many agent ase and onsider its ontinuum limit. A goal of su h an approa h is to develop wave equa- tions that may explain observed phenomenon of information transfer in natural o ks. We only onsider the y li graph of intera tion that enables us to present the results in a ompa t way. A general optimal ontrol problem in the loop group ase is stated. General ontrollability result in innite dimensional setting is shown to be helpful to onstru t su h optimal ontrol problem. The ne essary onditions for optimality, namely the Pontryagin's Maximum Prin iple (PMP) in Hilbert spa e setting is only valid under a ondition of nite o-dimensionality of a rea hable set. Two spe ial ases of this problem are studied. The ase in whi h the underlying group is the Heisenberg group H(3), i.e. a ontinuum of nonholo- nomi integrators is studied in detail. We have shown that the optimal ontrol solutions possess traveling wave hara ter. Moreover, a syn hronization result is obtained in whi h the innite oupling strength prohibits every agent in the o k to behave dierently. The ase of planar ontinuum, i.e. agents moving in the spe ial Eu lidean group SE(2) is also onsidered. Optimal ontrol evolution equa- tions are obtained by both al ulus of variations and PMP approa h. Similar to the H(3) ase, syn hronization result is obtained. Numeri al simulations for both these ases are presented. However, we have not been able to perform a thorough 150 analyti al study of the partial dierential equations obtained in the SE(2) ase. It is one of the dire tions in whi h future resear h ould be ondu ted. A riti al question will answer whether traveling wave solutions exist in this ase. Further- more, dierent intera tion graphs an be onsidered to un over more interesting details about this problem. Chapter 5 presents a data-driven analysis of ight data of European Starling o ks, aptured in Rome. This data gives rise to temporal signals of the o k's energy distribution in several energy modes. We use an optimal ontrol based data-tting te hnique to explain this data as the out ome of an evolutionary game on a simplex. We all the data-tting ost fun tionals of the underlying optimal ontrol problem as ` ognitive ost' that measures the ognitive eort of the o k to allo ate its energy in dierent modes. In our work, we have only onsidered energy splitting into two modes so as to onsider a simple game on the one-dimensional simplex. This an be extended to higher dimensions where multiple energy modes are onsidered. 151 Appendix A An Optimal Control Problem in an Innite Dimensional Setting A.1 Introdu tion After Pontryagin provided his method for solving optimal ontrol problems in - nite dimensional setting [Pontryagin et al., 1962℄, there have been many attempts to try and prove similar prin iples in innite dimensions. However, the ounterex- ample of Egorov [Egorov, 1963℄ posed a serious hallenge to that pro ess. This ounterexample showed that the Pontryagin's maximum prin iple does not gen- erally hold in innite dimensions. In parti ular, the ostate variable an be ome identi ally zero, making the maximum prin iple empty. The advan ements in the following de ades [Ekeland, 1979; Fattorini, 1987; Li and Yong, 2012; Krastanov et al., 2011℄ showed that it is possible to state PMP in some ases where some additional assumptions are made. In this work, we adopt a similar path to prove 152 the maximum prin iple set in a mu h friendlier setting. We onsider an abstra t dierential equation in a Hilbert spa e, dx(t) = f(t, x(t), u(t)), a.e. in [0, T ], (A.1) dt where x(t) ∈ X , u(·) ∈ U , and T > 0. Let X be a Hilbert spa e alled the state spa e and U be the set of all measurable fun tions u(·) : [0, T ] → U , where U is a separable metri spa e alled the ontrol spa e. With this setup, we formulate the following optimal ontrol problem (P), ∫ T min J(u) = L(t, x(t), u(t))dt u∈U 0 (P) (A.2) subje t to: ẋ = f(t, x, u), a.e. in [0, T ], x(0) = x0, x(T ) = xT . We assume that both the fun tions f(·, ·, ·) and L(·, ·, ·) are Bo hner integrable in t ∈ [0, T ] and Lips hitz ontinuous in x(t) ∈ X , with onstant K. Further- ′ more, we require the existen e of the ontinuous Fré het derivatives fx(t, x, u) and L′x(t, x, u) ′ ′ . We also assume the fun tions f, L and their derivatives fx, Lx to be bounded, i.e. there exists an M > 0, su h that ‖f(t, x, u)‖ ≤ M, ‖f ′x(t, x, u)‖ ≤ M, (A.3) ‖L(t, x, u)‖ ≤ M, ‖L′x(t, x, u)‖ ≤ M, for all (t, x(t), u(t)) ∈ [0, T ]× X × U . Note that these hypotheses ensure a on- tinuous and unique solution of (A.1) to exist [Avez, 1986℄. Let the spa e U be endowed with the distan e fun tion, d(u, v) = meas{t : u(t) 6= v(t)}, (A.4) 153 where meas{·} denotes the usual Lebesgue measure in [0, T ]. Then, a ording to Theorem 5.3 of [Fattorini, 1987℄, the spa e U is omplete with respe t to the distan e d. A dire t onsequen e of the assumptions stated above leads to an important result [Krastanov et al., 2011℄. Lemma A.1.1. Let u1(·), u2(·) be any two arbitrary members of U . Denote the state traje tories asso iated with these ontrols by, xi(·) = x(·, ui(·)), i = 1, 2. Then there exist positive onstants C1, C2 su h that, sup ‖x1(t)− x2(t)‖ ≤ C1d(u1, u2), (A.5) t∈[0,T ] |J(u1)− J(u2)| ≤ C2d(u1, u2). (A.6) Proof. Let S ⊂ [0, T ] be the set where the ontrols u1 and u2 dier, i.e. d(u1, u2) = meas{S}. We know that (A.1) an also be written as, ∫ t x(t) = x0 + f(s, x(s), u(s))ds. 0 Then, ∫ t x1(t)− x2(t) = (f(s, x1(s), u1(s))− f(s, x1(s), u1(s))) ds 0 ∫ t = (f(s, x1(s), u1(s))− f(s, x2(s), u1(s))) ds 0 ∫ t + (f(s, x2(s), u1(s))− f(s, x2(s), u2(s))) ds 0 ∫ t = (f(s, x1(s), u1(s))− f(s, x2(s), u1(s))) ds 0 ∫ + (f(s, x2(s), u1(s))− f(s, x2(s), u2(s))) ds. [0,t]∩S 154 Taking norm of both sides we get, ∫ t ‖x1(t)− x2(t)‖ ≤ K ‖x1(s)− x2(s)‖ ds+ 2Md(u1, u2), 0 where the Lips hitz property of f in x is used in the rst term and the uni- form boundedness is used in the se ond term. By the use of Gronwall-Bellman inequality, we arrive at (A.5). The result (A.6) an be derived analogously.  We now pro eed to solve this problem by a maximum prin iple based approa h. Before we state the maximum prin iple, let us give some additional te hni al details that are going to be essential for the proof of the maximum prin iple. Denition A.1. (Finite Codimensionality) [Fattorini, 1987℄ A subset S of a Hilbert spa e Z is alled to be nite odimensional in Z, if there exists a losed subspa e Zc ⊆ Z of nite odimension su h that Sc = Π( o(S)), has nonempty interior in Zc, where Πc denotes the orthogonal proje tion from Z onto Zc and o means losed onvex hull. We will now make a key assumption to derive a nontrivial maximum prin iple. ∗ Let a solution of problem (P) exist and that optimal ontrol is denoted as u ∈ U ∗ and let the orresponding optimal traje tory be denoted as x (t). Then dene the `rea hable set' as, { ∫ t R := z(T ) ∈ X | z(t) = f ′ (s, x∗x (s), u∗(s)) · z(s)ds 0 t }∫ + (f(s, x∗(s), v(s))− f(s, x∗(s), u∗(s))) ds, for some v(·) ∈ U (A.7) 0 (A1) The set R is nite odimensional in X . 155 A.2 Maximum Prin iple Using usual formalism, we invoke the pre-Hamiltonian fun tion H : R×X ×U × R×X∗ → R as, H(t, x(t), u(t), p0, p(t)) = p0L(t, x(t), u(t)) + 〈p(t), f(t, x(t), u(t))〉 , (A.8) where p(t) ∈ X∗ is alled the ostate variable. Intuitively, we want to make the following statement of the maximum prin iple that needs to be validated. ∗ Theorem A.2.1. (Maximum Prin iple) Let u ∈ U be an optimal ontrol for ∗ problem (P) and x (t) be the orresponding optimal traje tory. Then, there exist (p∗, p∗(t)) ∈ R × X∗a pair 0 , t ∈ [0, T ] ∗ ∗, su h that (p0, p ) ≡6 (0, 0), p∗0 ≤ 0 p∗, (·) satises the dierential equation, ṗ∗(t) = − (f ′ ⋆(t, x∗(t), u∗(t)) p∗(t)− p∗L′ (t, x∗(t), u∗x 0 x (t)), (A.9) ⋆ where by A we denote the adjoint operator of the operator A. The pointwise maximization of the pre-Hamiltonian holds, H(t, x∗(t), u∗(t), p∗0, p ∗(t)) = maxH(t, x∗(t), v, p∗0, p ∗(t)), (A.10) v∈U ∗ ∗ Moreover, x and p satisfy Hamilton's anoni al equations, i.e. dx∗ δH = ∗ (t, x ∗, u∗, p∗ ∗0, p ) (A.11)dt δp dp∗ −δH= (t, x∗, u∗∗ , p ∗, p∗0 ). (A.12)dt δx 0 Proof. At the outset, we begin by introdu ing the variable, x (t) ∈ R that obeys   x0(t)   0 0   the dynami s, ẋ = L(t, x, u), x (0) = 0. Dene, y(t) = ∈ R × X, so   x(t) 156 that   dy(t) L(t, x, u)  =   =: g(t, y, u), a.e. in [0, T ]. (A.13) dt   f(t, x, u) The ore of the argument in proving the maximum prin iple will follow re- sults of Ekeland [Ekeland, 1979℄, and te hniques developed in [Fattorini, 1987; Krastanov et al., 2011; Li and Yong, 2012℄. We will now state a result known as Ekeland variational prin iple [Ekeland, 1979℄. Lemma A.2.1. (Ekeland Variational Prin iple) Let V be a omplete metri spa e with respe t to the distan e fun tion d(·, ·) and let F : V → R ∪ {+∞} be lower semi ontinuous and bounded below with F 6≡ +∞. Let ǫ > 0 and u ∈ V be su h that F (u) ≤ inf{F (w) : w ∈ V }+ ǫ. (A.14) Then there exists v ∈ V su h that √ d(u, v) ≤ ǫ (A.15) √ F (w)− F (v) ≥ − ǫ d(w, v), ∀w ∈ V. (A.16) Let us pro eed by assuming that an optimal ontrol to the problem (P) exists ∗ ∗ 0,∗ ∗ and is denoted by u and let y = (x , x ) be the orresponding optimal traje - ∗ tory. We write the minimum ost by η0, i.e. η0 = J(u ). Now, for given ǫ > 0, we onsider the fun tion Jǫ : U → R, √ Jǫ(u) = (J(u)− η0 + ǫ)2 + ‖x(T )− xT ‖2X . (A.17) 157 It is evident that Jǫ(u) > 0 for all u ∈ U and ǫ > 0. Moreover, J (u∗ǫ ) = ǫ ≤ inf Jǫ(u) + ǫ, u∈U ǫ whi h by the Ekeland variational prin iple yields the existen e of u ∈ U su h that √ d(uǫ, u∗) ≤ ǫ, (A.18) √ J (w)− J (uǫǫ ǫ ) ≥ − ǫ d(w, uǫ), ∀w ∈ U . (A.19) ǫ Next, we introdu e a variation in ontrol u what is known as needle variations". For any v(·) ∈ U , let h : [0, T ] → R×X ,   L(t, xǫ(t), v(t))− L(t, xǫ(t), uǫ(t))   h(t) = (g(t, xǫ(t), v(t))− g(t, xǫ(t), uǫ(t))) =   .   f(t, xǫ(t), v(t))− f(t, xǫ(t), uǫ(t)) (A.20) Then, a ording to Corollary 3.9 (p. 144) of [Li and Yong, 2012℄ , for any ρ ∈ (0, 1], there is a measurable set Fρ ⊂ [0, T ] su h that meas{Fρ} = ρT and ∫ t ∫ ρ h(s)ds = h(s)ds+ o(ρ), (A.21) 0 Fρ∩[0,t] o(ρ) where → 0 as ρ ↓ 0, uniformly in t ∈ [0, T ]. The perturbed ontrol is then ρ dened as,    uǫ(t), t ∈/ Fρ, uǫρ(t) = . (A.22)   v(t), t ∈ Fρ ǫ It is of interest to express the perturbation in traje tory when the ontrol uρ is ǫ ǫ applied, i.e. we want a Taylor like expansion of yρ(t) = y(t, uρ) with respe t to ρ 158 at ρ = 0 ǫ, i.e. at y (t) = y(t, uǫ). Let us write, yǫρ(t) = y ǫ(t) + ρψǫ(t) + o(ρ). Then, yǫρ(t)− yǫ(t) ψǫ(t) = lim ρ↓0 ρ ( ) ∫ t ∫1 t = lim g(s, yǫ(s), uǫρ ρ(s))ds− g(s, yǫ(s), uǫ(s))ds ρ↓0 ρ 0 0 ∫ t g(s, yǫ(s), uǫρ ρ(s))− g(s, yǫ, uǫρ(s)) = lim ds ρ↓0 0 ρ ∫ g(s, yǫ(s), v(s))− g(s, yǫ, uǫ(s)) + lim ds ρ↓0 ρ [0,t]∩Fρ ∫ t ∫ t = g′y(s, y ǫ(s), uǫ(s)) · ψǫ(s)ds+ (g(s, yǫ(s), v(s))− g(s, yǫ, uǫ(s))) ds, 0 0 (A.23) ′ where the se ond term follows from (A.21). gy is the Fré het derivative of g with respe t to y and an be de omposed as,   L′x(t, x ǫ(t), uǫ(t)) · q g′   (t, yǫy (t), u ǫ(t)) · q̄ =   , for any q̄ = (q0, q) ∈ R×X.   f ′x(t, x ǫ(t), uǫ(t)) · q ∣ ǫ ǫ ǫ ǫ d ǫ ∣ Let's write ψ (t) = (z0(t), z (t)). In parti ular, we have z0(T ) = J(u ) anddρ ρ ρ=0 ∣ zǫ(t) = d xǫ ∣ρ(t)) , ea h of whi h an be spelled out separately from (A.23),dρ ρ=0 ∫ T ∫ T zǫ(T ) = L′0 x(s, x ǫ(s), uǫ(s)) · zǫ(s)ds+ (L(s, xǫ(s), v(s))− L(s, xǫ(s), uǫ(s))) ds, 0 0 (A.24) ∫ t ∫ t zǫ(t) = f ′x(s, x ǫ(s), uǫ(s)) · zǫ(s)ds+ (f(s, xǫ(s), v(s))− f(s, xǫ(s), uǫ(s))) ds. 0 0 (A.25) 159 ǫ ǫ Now, as a next step, we derive ne essary onditions for the pair (y (t), u (t)) to be suboptimal. We do this by using the Ekeland variational prin iple and letting ρ tend to zero. In (A.19), we set w = uǫ ǫ ǫρ. Note that the ontrols u and uρ dier only in the set Fρ, whi h has a measure ρT . Then, √ Jǫ(u ǫ ρ)− Jǫ(uǫ) ≥ − ǫd(uǫρ, uǫ) ≥ −√ǫρT. i.e. Jǫ(u ǫ ǫ ρ)− Jǫ(u ) ≥ − √T ǫ. (A.26) ρ Now, note that, ∣ Jǫ(u ǫ ρ)− J ǫǫ(u ) dJ ǫ ∣ǫ(uρ) ∣ lim = ∣ ρ↓0 ρ dρ ∣ ρ=0 [ 1 dJ(uǫρ) = 2(J(uǫ )− η0 + ǫ) 2Jǫ(uǫρ) ρ dρ ∣ ( ǫ )] ∣ ∥ ∥ ∥ ∥ dxρ(T ) ∣ +2 ∥xǫ ǫ∥ ∥ ∥ρ(T )− xT xρ(T )− xT · ∣dρ ∣ ρ=0 〈 〉 (J(uǫ)− η0 + ǫ) xǫ(T )− xT = zǫ0(T ) + , z ǫ(T ) . (A.27) J (uǫǫ ) Jǫ(uǫ) Thus, taking the limit in (A.26), we an write, √ ξǫzǫ0 0(T ) + 〈ξǫ, zǫ(T )〉 ≥ −T ǫ, (A.28) ǫ ǫ ξǫ = (J(u )−η0+ǫ) ξǫ = x (T )−xT ∈ X∗where 0 ǫ and ǫ . Note additionally that,Jǫ(u ) Jǫ(u ) (ξǫ)2 + ‖ξǫ‖20 = 1. (A.29) 160 (yǫEquation (A.28) an be regarded as the ne essary onditions for (t), uǫ(t)). Finally, we will let ǫ ∗ ∗tend to zero to obtain ne essary onditions for (y (t), u (t)) to be optimal. Dene, ∫ T ∫ T z0 := L ′ x(s, x ∗(s), u∗(s)) · z(s)ds+ (L(s, x∗(s), v(s))− L(s, x∗(s), u∗(s))) ds, 0 0 (A.30) ∫ t ∫ t z(t) := f ′x(s, x ∗(s), u∗(s)) · z(s)ds + (f(s, x∗(s), v(s))− f(s, x∗(s), u∗(s))) ds. 0 0 (A.31) Sin e v(·) is any arbitrary element in U , z(T ) ∈ R, .f. (A.7). Lemma A.2.2. The following results hold true. lim |zǫ0(T )− z0| = 0, ǫ↓0 (A.32) lim sup ‖zǫ(t)− z(t)‖ = 0, ǫ↓0 t∈[0,T ] ǫ ∗ Proof. Let us denote Sǫ = {t ∈ [0, T ] : u (t) 6= u (t)}. Then, meas{Sǫ} = √ d(uǫ, u∗) ≤ ǫ, by (A.18). From the denition of z(t), we nd, ∫ t ‖z(t)‖ ≤ M ‖z(s)‖ ds+ 2MT, 0 ′ where boundedness of both f and fx have been used. Applying the Gronwall- MT Bellman inequality, we get ‖z(t)‖ ≤ 2MTe , for all t ∈ [0, T ]. We now write, ∫ t zǫ(t)− z(t) = f ′ (s, xǫ(s), uǫ(s)) · (zǫx (s)− z(s))ds 0 ∫ t + (f ′ (s, xǫx (s), u ǫ(s))− f ′x(s, x∗(s), u∗(s))) · z(s)ds 0 ∫ t + (f(s, xǫ(s), v(s))− f(s, x∗(s), v(s))) ds 0 ∫ t + (f(s, xǫ(s), uǫ(s))− f(s, x∗(s), u∗(s))) ds. 0 161 ′ Taking norm on both sides and utilizing boundedness of fx and z(t), we obtain, ∫ t ‖zǫ(t)− z(t)‖ ≤ M ‖zǫ(s)− z(s)‖ ds 0 ∫ t + 2MTeMT ‖f ′x(s, xǫ(s), uǫ(s))− f ′ (s, x∗x (s), u∗(s))‖ ds 0 ∫ t + ‖f(s, xǫ(s), v(s))− f(s, x∗(s), v(s))‖ ds 0 ∫ t + ‖f(s, xǫ(s), uǫ(s))− f(s, x∗(s), u∗(s))‖ ds. (A.33) 0 The last term in (A.33) an be written as, ∫ t ‖f(s, xǫ(s), uǫ(s))− f(s, x∗(s), u∗(s))‖ ds 0 ∫ = ‖f(s, xǫ(s), u∗(s))− f(s, x∗(s), u∗(s))‖ ds [0,t]\Sǫ ∫ + ‖f(s, xǫ(s), uǫ(s))− f(s, x∗(s), u∗(s))‖ ds [0,t]∩Sǫ ∫ ≤ K ‖xǫ(s)− x∗(s)‖ ds+ 2Md(uǫ, u∗) [0,t]\Sǫ ≤ √ ǫ↓0(KC ǫ ∗1T + 2M)d(u , u ) ≤ (KC1T + 2M) ǫ −→ 0. Note that we have used the Lips hitz property of f and result of Lemma A.1.1. The se ond and third term an be treated in a similar fashion to show they are of o(1) whi h goes to 0 as ǫ tends to 0. Note that, instead of Lips hitz ontinuity, we ′ would use ontinuity of fx in x in order to use appropriate upper bound. Hen e, (A.33) an be written as, ∫ t ‖zǫ(t)− z(t)‖ ≤ M ‖zǫ(s)− z(s)‖ ds+ o(1), 0 162 whi h by Gronwall-Bellman inequality yields, ‖zǫ(t)− z(t)‖ ≤ eMT −ǫ↓0o(1) → 0, uniformly in t ∈ [0, T ] ǫ. The onvergen e of z0(T ) an be shown analogously.  Using (A.28) and the Cau hy-S hwarz inequality, we may now write, ξǫ0z0 + 〈ξǫ, z(T )〉 = ξǫ0zǫ0 + 〈ξǫ, zǫ(T )〉 − ξǫ ǫ0 (z0 − z0)− 〈ξǫ, zǫ(T )− z(T )〉 ≥ − √T ǫ− |ξǫ0| |zǫ0 − z | − ‖ξǫ0 ‖ ‖zǫ(T )− z(T )‖ ≥ − √T ǫ− |zǫ0 − z0| − ‖zǫ(T )− z(T )‖ . (A.34) ǫ 2 ǫ 2 ǫ ǫ The last inequality follows, sin e (ξ0) + ‖ξ ‖ = 1, both |ξ0| ≤ 1 and ‖ξ ‖ ≤ 1. √ ǫ ǫ ǫ Denote, κ = −T ǫ − |z0 − z0| − ‖z (T )− z(T )‖ and by the onvergen e results ǫ (A.32), we see that κ → 0 as ǫ ↓ 0. Thus, (A.34) an be expressed as, ξǫz + 〈ξǫ0 0 , z〉 ≥ −κǫ, ∀z0 ∈ R, z ∈ R, (A.35) ǫ ǫ↓0 where, κ −→ 0. Now the assumption (A1) that the set R is nite odimensional ǫ ǫ in X is going to be useful in proving nontriviality of the limit of the pair (ξ0, ξ ) as ǫ goes to 0. Here we state the following lemma from [Fattorini, 1987℄, as a onsequen e of nite odimensionality. Lemma A.2.3. Let Q be a set of nite odimension in a Hilbert spa e Z and let {zn} be a sequen e of ve tors in Z su h that 0 < c ≤ ‖zn‖ ≤ C. 163 Assume additionally that, 〈zn, q〉 ≥ −θn, for q ∈ Q and θn → 0 as n → ∞. Then there exists a subsequen e of {zn} that onverges weakly to z ∈ Z, and z 6= 0. Now hoose a sequen e {ǫ(n)} su h that ǫ(n) → 0 as n → ∞. Sin e both ǫ(n) sequen es {ξ0 } {ξǫ(n) ǫ(nk)and } are bounded, there exist subsequen es {ξ } and {ξǫ(nk)} ∗that onverge weakly to some ξ̄0 ∈ R and ξ̄ ∈ X . For simpli ity, let the subsequen es be denoted by themselves. Showing nontriviality of the pair (ξ̄0, ξ̄) is a ru ial step in proving maximum prin iple in innite dimensional ase. Re all ξǫ = (J(u ǫ)−η0+ǫ) that, 0 ǫ , so that ξ ǫ > 0, ∀ǫ > 0. Hen e, we may only have ξ̄ ≥ 0. J 0ǫ(u ) 0 If ξ̄0 6= 0 ǫ, we are done proving that (ξ̄0, ξ̄) 6= (0, 0). Otherwise, let ξ0(n) → 0 as ( ) ∥ ∥ → ∞ ǫ(n) 2 2 ǫ(n) n ∥ ∥. Then from the relation (A.29), we get, 1 ≥ ξ = 1 − ξ0 ≥ 1 − δ > 0, for some δ > 0, for n large enough. Finally, by the lemma A.2.3, we get ξ̄ 6= 0 in the ase ξ̄0 = 0. Hen e we on lude that, ∗ (ξǫ ǫ0, ξ ) ⇀ (ξ̄0, ξ̄) 6= (0, 0), ξ̄0 ≥ 0. (A.36) Then, nally taking the limit ǫ ↓ 0 in (A.35), we get for any z ∈ R and z0 as spe ied in (A.30), there exists a pair R×X∗ ∋ (ξ̄0, ξ̄) 6= (0, 0), ξ̄0 ≥ 0, so that, 〈 〉 ξ̄0z0 + ξ̄, z ≥ 0. (A.37) ∗ ∗ Now, let us introdu e the ostate variable p (t) ∈ X , that obeys the following dierential equation, ṗ∗ ⋆ (t) = − (f ′x(t, x∗(t), u∗(t)) p∗(t)− p∗0L′x(t, x∗(t), u∗(t)), (A.38) 164 ∗ with p (T ) = −ξ̄ and p∗0 = −ξ̄0 ≤ 0 ⋆. Here by A we denote the adjoint operator of the operator A. Now sin e z(0) = 0, we have, 〈p∗(T ), z(T )〉 = 〈p∗(T ), z(T )〉 − 〈p∗(0), z(0)〉 ∫ T d = 〈p∗(t), z(t)〉 dt 0 dt ∫ T [〈 = − (f ′ ⋆ 〉 x(t, x ∗(t), u∗(t)) p∗(t)− p∗L′0 x(t, x∗(t), u∗(t)), z(t) 0 + 〈p∗(t), f ′ ∗ ∗x(t, x (t), u (t))z(t) + f(t, x∗(t), v(t))− f(t, x∗(t), u∗(t))〉] dt ∫ T = −p∗0 〈L′x(t, x∗(t), u∗(t)), z(t)〉 dt 0 ∫ T + (〈p∗(t), f(t, x∗(t), v(t))− f(t, x∗(t), u∗(t))〉) dt (A.39) 0 This, ombined with the denition of z0 (A.30) and equation (A.37) yields, 〈p∗(T ), z(T )〉+ p∗0z0 ∫ T = [H(t, x∗(t), v(t), p∗ ∗0, p (t))−H(t, x∗(t), u∗(t), p∗0, p∗(t))] dt ≤ 0. (A.40) 0 Sin e the ontrol set U is separable, the similar argument as in [Krastanov et al., 2011; Li and Yong, 2012℄ would give the pointwise maximization riterion of the pre-Hamiltonian, H(t, x∗(t), v(t), p∗0, p ∗(t)) ≤ H(t, x∗(t), u∗(t), p∗0, p∗(t)), a.e. in [0, T ], ∀v ∈ U . (A.41) From the denition of the Hamiltonian, we an nally ompute its derivatives. In what follows, the appropriate arguments will be suppressed for simpli ity and the notation |∗ will imply the fun tion has been evaluated at optimal parameters. We 165 then nd for any (x̃, p̃) ∈ X ×X∗, ∣ ∗ − ∗ ∣δH ∣ · H(p + σp̃) H(p ) ∣ ∣ p̃ = lim ∣ δp∗ ∣ ∣∗ σ→0 σ ∗ ∣ d ∗ ∣= H(p + σp̃)∣ dσ ∣σ=0,∗ = 〈p̃, f |∗〉 and, ∣ δH H(x∗∣ + σx̃)− ∣H(x∗) ∣ ∣ ∣ δx∗ · x̃ = lim ∣ σ→0 σ ∣∗ ∗ ∣ d ∣ = H(x∗ + σx̃)∣ dσ ∣σ=0,∗ ∣ d ∣ = (〈p∗, f(x∗ + σx̃)〉+ p∗0L(x∗ + σx̃)) ∣dσ ∣σ=0,∗ = 〈p∗, f ′x|∗ · x̃〉+ p∗0 〈L′x|∗, x̃〉 〈 ′ | ⋆ 〉 = (fx ∗) p ∗ + p∗L′0 x|∗, x̃ . Thus, we may write the anoni al Hamilton's equations of motion, dx∗ δH = ∗ (t, x ∗, u∗, p∗ ∗ dt δp 0 , p ), (A.42) dp∗ −δH= ∗ (t, x ∗, u∗, p∗, p∗). dt δx 0 This ompletes the proof of the maximum prin iple.  166 Bibliography Agra hev, A. and Sa hkov, Y. 2004. Control Theory from the Geometri Viewpoint. Springer, Berlin. Attanasi, A., Cavagna, A., Del Castello, L., Giardina, I., Grigera, T. S., Jeli¢, A., Melillo, S., Parisi, L., Pohl, O., Shen, E., and Viale, M. 2014. Information transfer and behavioural inertia in starling o ks. Nature physi s 10:691696. Austin, M. A., Krishnaprasad, P. S., and Wang, L.-S. 1993. Almost poisson integration of rigid body systems. Journal of Computational Physi s 107:105117. Avez, A. 1986. Dierential al ulus. John Wiley & Sons In . Ballerini, M., Cabibbo, N., Candelier, R., Cavagna, A., Cisbani, E., Giardina, I., Le omte, V., Orlandi, A., Parisi, G., Pro a ini, A., Viale, M., and Zdravkovi , V. 2008a. Intera tion ruling animal olle tive behavior depends on topologi al rather than metri distan e: Eviden e from a eld study. Pro eedings of the National A ademy of S ien es 105:12321237. Ballerini, M., Cabibbo, N., Candelier, R., Cavagna, A., Cisbani, E., Giardina, I., Orlandi, A., Parisi, G., Pro a ini, A., Viale, M., et al. 2008b. Empiri al investigation of starling o ks: a ben hmark study in olle tive animal behaviour. Animal behaviour 76:201215. Bellaï he, A., Jean, F., and Risler, J. J. 1998. Geometry of nonholonomi systems. In J.-P. Laumond (ed.), Robot Motion Planning and Control. Springer, London. Bialek, W., Cavagna, A., Giardina, I., Mora, T., Silvestri, E., Viale, M., and Wal zak, A. M. 2012. Statisti al me hani s for natural o ks of birds. Pro eedings of the National A ademy of S ien es 109:47864791. Birkhoff, G. D. 1915. The restri ted problem of three bodies. Rendi onti del Cir olo Matemati o di Palermo 39:265334. Bishop, R. L. 1975. There is more than one way to frame a urve. Ameri an Mathemati al Monthly pp. 246251. 167 Bliss, G. A. 1946. Le tures on the al ulus of variations, phoenix s ien e series. The University of Chi ago Press, Chi ago 37:16. Bro kett, R. 1973. Lie theory and ontrol systems dened on spheres. SIAM Journal on Applied Mathemati s 25:213225. Byrd, P. F. and Friedman, M. D. 1971. Handbook of Ellipti Integrals for Engineers and S ientists. Springer Berlin Heidelberg. Cavagna, A., Cimarelli, A., Giardina, I., Parisi, G., Santagati, R., Stefanini, F., and Viale, M. 2010. S ale-free orrelations in starling o ks. Pro eedings of the National A ademy of S ien es 107:1186511870. Cavagna, A., Conti, D., Giardina, I., and Grigera, T. S. 2018. Propagat- ing speed waves in o ks: a mathemati al model. Physi al Review E 98:052404. Chiu, C., Reddy, P. V., Xian, W., Krishnaprasad, P. S., and Moss, C. F. 2010. Ee ts of ompetitive prey apture on ight behavior and sonar beam pattern in paired big brown bats, Eptesi us Fus us. Journal of Experimental Biology 213:33483356. Chorin, A. J., Marsden, J. E., and Marsden, J. E. 1990. A mathemati al introdu tion to uid me hani s, volume 3. Springer. Corbet, P. S. 1999. Dragonies: Behavior and E ology of Odonata. Comsto k Publishing Asso iates. Cu ker, F. and Smale, S. 2007. Emergent behavior in o ks. IEEE Transa - tions on automati ontrol 52:852862. Davis, H. T. 1962. Introdu tion to nonlinear dierential and integral equations, volume 971.Dover Publi ation, In . Dey, B. 2015. Re onstru tion, Analysis and Synthesis of Colle tive Motion. PhD thesis, University of Maryland, College Park. Dey, B. and Krishnaprasad, P. S. 2012. Traje tory smoothing as a linear optimal ontrol problem. In Pro eedings of Annual Allerton Conferen e on Communi ation, Control, and Computing, pp. 14901497, Urbana-Champaign, IL. IEEE. Dey, B. and Krishnaprasad, P. S. 2014. Control-theoreti data smoothing. In Pro eedings of 53rd IEEE Conferen e on De ision and Control, pp. 5064 5070. IEEE. Dubins, L. E. 1957. On urves of minimal length with a onstraint on aver- age urvature, and with pres ribed initial and terminal positions and tangents. Ameri an Journal of Mathemati s 79:497516. 168 Ebin, D. G. and Marsden, J. 1970. Groups of dieomorphisms and the motion of an in ompressible uid. Annals of Mathemati s pp. 102163. Eells Jr, J. 1966. A setting for global analysis. Bulletin of the Ameri an Mathemati al So iety 72:751807. Egorov, Y. V. 1963. The optimal ontrol in bana h spa e. Dokl. Akad. Nauk SSSR 150:241244. Ekeland, I. 1979. Non onvex minimization problems. Bulletin of the Ameri an Mathemati al So iety 1:443474. Elsgol , L. D. 2012. Cal ulus of variations. Courier Corporation. Euler, L. 1744. Methodus inveniendi lineas urvas maximi minimive proprietate gaudentes sive solutio problematis isoperimetri i latissimo sensu a epti. Fattorini, H. 1987. A unied theory of ne essary onditions for nonlinear non- onvex ontrol systems. Applied Mathemati s and Optimization 15:141185. Frai hard, T. and S heuer, A. 2004. From Reeds and Shepp's to ontinuous- urvature paths. IEEE Transa tions on Roboti s 20:10251035. Galloway, K., Justh, E., and Krishnaprasad, P. 2009. Geometry of y li pursuit. In Pro eedings of 48th IEEE Conferen e on De ision and Control, IEEE, pp. 36293643. Galloway, K. S., Justh, E. W., and Krishnaprasad, P. S. 2013. Sym- metry and redu tion in olle tives: y li pursuit strategies. Pro eedings of the Royal So iety A 469:20130264. Gelfand, I. and Fomin, S. 1963. Cal ulus of variations. Revised English edition translated and edited by Ri hard A. Silverman. Prenti e-Hall In ., Englewood Clis, NJ. Ghose, K., Horiu hi, T. K., Krishnaprasad, P. S., and Moss, C. F. 2006. E holo ating bats use a nearly time-optimal strategy to inter ept prey. PLoS Biology 4:865873. Giaquinta, M. and Hildebrandt, S. 1996. Cal ulus of Variations: Cal ulus of Variation. I. Springer-Verlag. Goldstein, H., Poole, C., and Safko, J. 2001. Classi al me hani s.Addison- Wesley. Halder, U. and Dey, B. 2015. Biomimeti algorithms for oordinated motion: theory and implementation. In Pro eedings of International Conferen e on Roboti s and Automation, pp. 54265432. IEEE. 169 Halder, U., Justh, E. W., and Krishnaprasad, P. S. 2019a. Continuum o king and ontrol. In preparation for submission. Halder, U. and Kalabi , U. 2017. Time-optimal solution for a uni y le path on se2 with a penalty on urvature. IFAC-PapersOnLine 50:63206325. Halder, U., Raju, V., Krishnaprasad, P. S., and et al. 2019b. Cognitive ost of o king: A geometri and hamiltonian perspe tive. In preparation for submission. Halder, U., S hlotfeldt, B., and Krishnaprasad, P. S. 2016. Steer- ing for bea on pursuit under limited sensing. In Pro eedings of 55th IEEE Conferen e on De ision and Control, pp. 38483855. IEEE. Heintze, E. and Liu, X. 1999. Homogeneity of innite dimensional isopara- metri submanifolds. Annals of mathemati s 149:149181. Inada, Y. and Kawa hi, K. 2002. Order and exibility in the motion of sh s hools. Journal of Theoreti al Biology 214:371  387. Justh, E. W. and Krishnaprasad, P. S. 2003. Steering laws and ontinuum models for planar formations. In Pro eedings of 42nd IEEE Conferen e on De ision and Control, pp. 36093615. IEEE. Justh, E. W. and Krishnaprasad, P. S. 2004. Equilibria and steering laws for planar formations. Systems & Control Letters 52:25  38. Justh, E. W. and Krishnaprasad, P. S. 2005. Natural frames and intera ting parti les in three dimensions. In Pro eedings of 44th IEEE Conferen e on De ision and Control, pp. 28412846. IEEE. Justh, E. W. and Krishnaprasad, P. S. 2006. Steering laws for motion amouage. Pro eedings of Royal So iety A 462:36293643. Justh, E. W. and Krishnaprasad, P. S. 2015a. Enlargement, geodesi s, and olle tives. In International Conferen e on Networked Geometri S ien e of Information, pp. 558565. Springer. Justh, E. W. and Krishnaprasad, P. S. 2015b. Optimality, redu tion and olle tive motion. Pro eedings of Royal So iety A 471:20140606. Justh, E. W. and Krishnaprasad, P. S. 2016. Subriemannian geodesi s for oupled nonholonomi integrators. In Pro eedings of Ameri an Control Conferen e, pp. 72817288. IEEE. Kine t. Mi rosoft kine t. http://www.xbox. om/en-US/xbox-one/ a essories/kine t-for-xbox-one. [Online℄. 170 Krastanov, M. I., Ribarska, N., and Tsa hev, T. Y. 2011. A pontryagin maximum prin iple for innite-dimensional problems. SIAM Journal on Control and Optimization 49:21552182. Krishnaprasad, P., Mar us, S. I., and Hazewinkel, M. 1983. Current algebras and the identi ation problem. Sto hasti s: An International Journal of Probability and Sto hasti Pro esses 11:65101. Krishnaprasad, P. S. 1993. Optimal ontrol and Poisson redu tion. Te hni al Report 93-87, Institute for Systems Resear h, University of Maryland. Kudrolli, A., Lumay, G., Volfson, D., and Tsimring, L. S. 2008. Swarm- ing and swirling in self-propelled polar granular rods. Physi al review letters 100:058001. LaValle, S. M. 2006. Planning Algorithms. Cambridge university press. Li, X. and Yong, J. 2012. Optimal ontrol theory for innite dimensional systems. Springer S ien e & Business Media. MATLAB. MATLAB Roboti s Toolbox. http://www.mathworks. om/ produ ts/roboti s/. [Online℄. Mis hiati, M. and Krishnaprasad, P. S. 2010. Motion amouage for ov- erage. In Pro eedings of Ameri an Control Conferen e, pp. 64296435. IEEE. Mis hiati, M. and Krishnaprasad, P. S. 2011. Mutual motion amouage in 3D. In Pro eedings of the 18th IFAC World Congress, pp. 44834488. Mis hiati, M. and Krishnaprasad, P. S. 2012. The dynami s of mutual motion amouage. Systems & Control Letters 61:894903. Mis hiati, M. and Krishnaprasad, P. S. 2017. Geometri de ompositions of olle tive motion. Pro eedins of Royal So iety A 473:20160571. Mizutani, A., Chahl, J. S., and Srinivasan, M. V. 2003. Inse t behaviour: Motion amouage in dragonies. Nature 423:604604. Mora, T. and Bialek, W. 2011. Are biologi al systems poised at riti ality? Journal of Statisti al Physi s 144:268302. Nagy, M., Akos, Z., Biro, D., and Vi sek, T. 2010. Hierar hi al group dynami s in pigeon o ks. Nature 464:890893. Nghiem, A. T., Auvinet, E., and Meunier, J. 2012. Head dete tion using Kine t amera and its appli ation to fall dete tion. In Pro eedings of 11th International Conferen e on Information S ien e, Signal Pro essing and their Appli ations, pp. 164169. IEEE. 171 Ohsawa, T. 2013. Symmetry redu tion of optimal ontrol systems and prin ipal onne tions. SIAM Journal on Control and Optimization 51:96120. Olberg, R. M., Worthington, A. H., and Venator, K. R. 2000. Prey pursuit and inter eption in dragonies. Journal of Comparative Physiology A 186:155162. Oliver, A., Kang, S., Wüns he, B. C., and Ma Donald, B. 2012. Using the Kine t as a navigation sensor for mobile roboti s. In Pro eedings of the 27th Conferen e on Image and Vision Computing, pp. 509514. ACM. OpenCV. http://open v.org/. [Online℄. Pioneer. Pioneer P3-DX. http://www.mobilerobots. om/resear hrobots/ pioneerp3dx.aspx. [Online℄. Pontryagin, L. S., Boltyanskiy, V. G., Gamkrelidze, R. V., and Mish henko, E. F. 1962. Mathemati al theory of optimal pro esses. In- ters ien e. Pressley, A. and Segal, G. B. 1986. Loop Groups. Clarendon Press. Raju, V. and Krishnaprasad, P. S. 2018. A variational problem on the probability simplex. In Pro eedings of 57th IEEE Conferen e on De ision and Control, pp. 35223528. IEEE. Rashevsky, P. 1938. About onne ting two points of a ompletely nonholonomi spa e by admissible urve. U h. Zapiski Ped. Inst. Libkne hta 2:8394. Reddy, P. V., Justh, E., and Krishnaprasad, P. S. 2006. Motion amou- age in three dimensions. In Pro eedings of 45th IEEE Conferen e on De ision and Control, pp. 33273332. IEEE. Reeds, J. and Shepp, L. 1990. Optimal paths for a ar that goes both forwards and ba kwards. Pa i Journal of Mathemati s 145:367393. ROS. Robot Operating System. http://www.ros.org. [Online℄. ROS-ARIA. http://wiki.ros.org/ROSARIA. [Online℄. Rund, H. 1966. The Hamilton-Ja obi theory in the al ulus of variations: its role in mathemati s and physi s. Krieger Pub Co. Salehani, M. K. and Markina, I. 2014. Controllability on innite-dimensional manifolds: a howrashevsky theorem. A ta Appli andae Mathemati ae 134:229246. Sussmann, H. J. and Tang, G. 1991. Shortest paths for the Reeds-Shepp ar: A worked out example of the use of geometri te hniques in nonlinear optimal ontrol. Te hni al Report 91-10, Rutgers Center for Systems and Control. 172 Svirezhev, Y. M. 1972. Optimum prin iples in geneti s. Studies on theoreti al geneti s pp. 86102. Thurrowgood, S., Moore, R. J., So ol, D., Knight, M., and Srini- vasan, M. V. 2014. A biologi ally inspired, vision-based guidan e system for automati landing of a xed-wing air raft. Journal of Field Roboti s 31:699727. Topaz, C. M., Bertozzi, A. L., and Lewis, M. A. 2006. A nonlo al ontin- uum model for biologi al aggregation. Bulletin of mathemati al biology 68:1601. Vásárhelyi, G., Virágh, C., Somorjai, G., Tar ai, N., Szörényi, T., Nepusz, T., and Vi sek, T. 2014. Outdoor o king and formation ight with autonomous aerial robots. In Pro eedings of IEEE/RSJ International Conferen e on Intelligent Robots and Systems, pp. 38663873. IEEE. Vi on. Vi on motion apture system. http://www.vi on. om. [Online℄. Vi sek, T., Czirók, A., Ben-Ja ob, E., Cohen, I., and Sho het, O. 1995. Novel type of phase transition in a system of self-driven parti les. Physi al Review Letters 75:12261229. Wei-Liang, C. 1939. Uber systeme von linearen partiellen dierentialglei hungen erster ordnung. Math. Ann 117:98105. Young, G. F., S ardovi, L., Cavagna, A., Giardina, I., and Leonard, N. E. 2013. Starling o k networks manage un ertainty in onsensus at low ost. PLoS omputational biology 9:e1002894. YouTube. Implementation videos. http://ter.ps/ICRA2015. [Online℄. Zhang, H.-P., Be’er, A., Florin, E.-L., and Swinney, H. L. 2010. Colle tive motion and density u tuations in ba terial olonies. Pro eedings of the National A ademy of S ien es 107:1362613630. 173