Verifiable swarm engineering with limited communication

University of Strathclyde Department of Mechanical and Aerospace Engineering Verifiable Swarm Engineering with Limited Communication by Giuliano Punzo A thesis presented in fulfilment of the requirements for the degree of Doctor of Philosophy 2013 Declaration of authors rights This thesis is the result of the authors original research. It has been composed by the author and has not been previously submitted for examination which has led to the award of a degree. The copyright of this thesis belongs to the author under the terms of the United Kingdom Copyright Acts as qualified by University of Strathclyde Regulation 3.50. Due acknowledgement must always be made of the use of any material contained in, or derived from, this thesis. Acknowledgements This thesis contains the work developed over 3 years of intense efforts. These efforts inevitably influenced my life beyond the research job and the one of people around me. My first and more sincere thank you goes to those people, in particular to those not involved in this thesis directly, for tolerating my behaviour. Elena, my Fiancée, standing always next to me in the crucial moments of my PhD, proved how our life together will be bright and reliable in a more convincing way I did in this work about the future of swarming systems. My family in Italy accepted my permanence in the United Kingdom not without sorrow, but I hope I managed to convince them that was the best thing for myself. Either way I owe them gratitude. I also want to thank those who keep on having my name on their close friend list without seeing me more than once per year. To my supervisor Malcolm Macdonald I acknowledge an insane bravery for having allowed me the freedom of developing my research the way I wanted, even deviating from what I suspect the original plans were. As he already knows this, the only thing it is left to say is thank you. Derek Bennet knows how important he was for the continuous support provided throughout three years of “wandering in the science”; if he didn’t, then this acknowledgement will make him know. I cannot forget the long coffee trail that linked my PhD to the one of Phil Karagiannakis, and the great work we managed to produce with the help of Stephan Weiss, another person who doesn’t deserve to be neglected here. We are all looking forward to see where this trail will take. A special thanks goes to Gordon Dobie, Charles Macleod and Rahul Summan sentinels and keepers of the secrets of the dungeon of the Electrical and Electronic ii Engineering department at Strathclyde University. I have to acknowledge them and the whole CUE research group for experimental work developed in this thesis. I would then like to thank Jules Simo for his inventive and original contribution to our work on swarm fragmentation. Naomi Leonard in Princeton will probably remember me every now and then when thinking to the fast manoeuvring or when thinking to take a city break in Glasgow. I aim to bother her together with George Young and Darren Pais again with the latest old news about eigenvectors and avian flocks. For the time being they fully deserve my gratitude and my best wishes. Finally I want to acknowledge the Advanced Space Concept Laboratory and the Department of Mechanical and Aerospace Engineering at the University of Strathclyde that provided me with the support to develop my research. Abstract The new paradigms of swarm engineering, distributed architectures and autonomous multi-agent systems, are foreseen to redefine the way many engineering problems are approached. The affirmation of these new concepts requires the complete understanding of complex dynamics by the designers. That is, any system whose concept departs from the monolithic architecture must deliver its tasks in a predictable way and be controlled in a safe manner, while keeping the maximum possible autonomy. This work aims to span the gap between a complete foreseeable behaviour and system autonomy using precise mathematical descriptions of the dynamics and control of multi-agent systems. Dynamical system theory, Lyapunov stability, linear algebra and graph theory are used to rigorously frame the problem and delineate the characteristics of such systems in relation to a number of applications and performance parameters. The work first considers multi-agent systems as multi-particle systems in a physics fashion to draw fundamental results about the robustness to fragmentation when the individuals do not benefit from all-to-all communications. The exploitation of limited communications together with artificial potential functions is shown to be an effective way to shape formations of agents in a range of applications for future engineering, and in particular this scheme is proved to be effective for space-based communications through the autonomous deployment of antenna arrays. In this context, application to robotics is explored through laboratory tests exploiting wheeled robots with possible applications to structural inspection or planetary exploration. A stable fractal formation is proved to emerge out of a number of agents whose interaction network presents a recursive layout whereby relative motion is driven by artificial potential functions. Finally, the fast manoeuvring problem is covered together with one of allocating resources in an efficient way to track an external signal for the benefit of the group as a whole. Through an algebraic approach, the tracking capabilities are distributed amongst the agents producing advantages at group level for the tracking of an external signal. This also translates into fast reaction to threats. Contents Contents v List of Figures viii List of Tables xii 1 Introduction 1 1.1 Towards Decentralised Architectures . . . . . . . . . . . . . . . . 1.1.1 Biological Aspects . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Technology Aspects . . . . . . . . . . . . . . . . . . . . . . 3 6 7 1.2 Self Organization, Group or Swarm Intelligence . . . . . . . . . . 1.2.1 Concept of Emergence . . . . . . . . . . . . . . . . . . . . 7 8 1.3 1.4 Challenges in modelling, simulating, implementing group Behaviours 9 From Particles Systems to Intelligent Cooperating Agents . . . . . 10 1.4.1 Stochastic Approach . . . . . . . . . . . . . . . . . . . . . 11 1.4.2 1.4.3 Deterministic Approach . . . . . . . . . . . . . . . . . . . Graph Theoretic Methods . . . . . . . . . . . . . . . . . . 12 13 1.5 1.4.4 Behavioral Rules and Algorithms . . . . . . . . . . . . . . 1.4.5 Other Approaches . . . . . . . . . . . . . . . . . . . . . . . Aims and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . 15 16 16 1.6 1.7 Thesis Layout and Methodology . . . . . . . . . . . . . . . . . . . Papers Authored . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 19 2 Control Problem, Theory and Tools 2.1 21 2.2 Dynamical system theory . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 System Energy and Hamiltonian . . . . . . . . . . . . . . . Control Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 22 23 2.3 Stability Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 v 2.3.1 2.4 2.5 2.6 2.7 Linear Stability . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3.2 Nonlinear Stability - Lyapunov stability . . . . . . . . . . Graph Theoretic Methods . . . . . . . . . . . . . . . . . . . . . . Nonsmooth Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 27 29 33 Swarm Modelling and Control . . . . . . . . . . . . . . . . . . . . 2.6.1 Artificial Potential Functions . . . . . . . . . . . . . . . . 34 35 2.6.2 Connection Networks . . . . . . . . . . . . . . . . . . . . . 2.6.3 The Combined Effect of APFs and Connection Network . . Schooling, Flocking and Swarming . . . . . . . . . . . . . . . . . . 39 43 44 2.7.1 2.7.2 Biological Studies . . . . . . . . . . . . . . . . . . . . . . . Crystalline Patterns . . . . . . . . . . . . . . . . . . . . . 44 46 2.7.3 2.7.4 Continuous Swarming Systems . . . . . . . . . . . . . . . . Consensus, a Common Problem . . . . . . . . . . . . . . . 49 50 3 The Cohesion of Emerging Patterns and the Emergence of Bot52 tlenecks 3.1 3.2 3.3 The Need for a Cohesive Swarm . . . . . . . . . . . . . . . . . . . A Minimum Number of Links for Cohesion . . . . . . . . . . . . . The Limit Case and the Cheeger Constant . . . . . . . . . . . . . 52 54 58 3.4 3.5 Effectiveness of the Model . . . . . . . . . . . . . . . . . . . . . . Final Considerations . . . . . . . . . . . . . . . . . . . . . . . . . 65 69 4 Fractal Fractionated Architectures and Arrays - Application in Space 4.1 Fractionated Spacecraft and Formation Flight . . . . . . . . . . . 4.1.1 Artificial Potential Functions . . . . . . . . . . . . . . . . 4.1.2 4.1.3 4.2 4.3 4.4 71 72 74 Adjacency Matrix . . . . . . . . . . . . . . . . . . . . . . . Distributed Fractal Antennas . . . . . . . . . . . . . . . . 75 78 The Control Law and the Fractal Antenna . . . . . . . . . . . . . 4.2.1 Central Symmetry Emergence . . . . . . . . . . . . . . . . 4.2.2 APF Coefficient Definition . . . . . . . . . . . . . . . . . . 81 81 87 4.2.3 4.2.4 Stability of Control Law . . . . . . . . . . . . . . . . . . . Fractal Antenna Array Design and Analysis . . . . . . . . 87 89 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . Final Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 97 5 Pattern Creation - Experimental Results 102 5.1 Real World Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.1.1 Navigation . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.1.2 Actuators . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.2 Swarm Robotic Experiments . . . . . . . . . . . . . . . . . . . . . 106 5.2.1 Testbed and Methodology . . . . . . . . . . . . . . . . . . 106 5.3 5.2.2 Guidance and Control . . . . . . . . . . . . . . . . . . . . 109 5.2.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . 115 Real and Simulated Real World . . . . . . . . . . . . . . . . . . . 119 6 Fast Consensus and Manoeuvring 123 6.1 6.2 The Need for Fast Manoeuvring . . . . . . . . . . . . . . . . . . . 123 Fast manoeuvring out of Network Structure . . . . . . . . . . . . 126 6.3 6.4 6.5 Eigenvectors and Eigenvalues in a Network System . . . . . . . . 127 Distributing Tracking Capabilities . . . . . . . . . . . . . . . . . . 130 A Good Choice with Limited Resources . . . . . . . . . . . . . . . 132 6.6 The Significance of the First Left Eigenvector and the Effects of Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.6.1 Limited Resources and Small Perturbations . . . . . . . . 139 6.7 6.8 6.9 6.6.2 Fixed Trace, a Misleading Effect . . . . . . . . . . . . . . . 142 Solution of the First Order Linear Dynamics . . . . . . . . . . . . 143 6.7.1 6.7.2 6.7.3 The Regular Lattice Case . . . . . . . . . . . . . . . . . . 144 The Ring Case . . . . . . . . . . . . . . . . . . . . . . . . 145 The Small World Network Case . . . . . . . . . . . . . . . 146 6.7.4 6.7.5 The Random Network Case . . . . . . . . . . . . . . . . . 146 Performance of First Order System . . . . . . . . . . . . . 151 The Second Order and the Nonlinear Dynamics . . . . . . . . . . 152 Fast Manoeuvring - Final Considerations . . . . . . . . . . . . . . 160 7 Conclusions and Future Work 7.1 A Summary of the Topics Covered 7.2 7.3 163 . . . . . . . . . . . . . . . . . 163 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 Bibliography 170 List of Figures 1.1 1.2 1.3 (a) The game of life as example of cellular automata as reported in Conway’s original work. (b) Beni’s original concept of cellular robotic systems as illustrated in its work of 1988 . . . . . . . . . . Examples of emergent patterns in nature . . . . . . . . . . . . . . 4 9 Illustration for the problem of the seven bridges of Königsberg and the associated graph . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1 (a) The control scheme of a monolithic system and (b) the control 25 2.2 scheme of a distributed scheme. . . . . . . . . . . . . . . . . . . . Time history for a single degree of freedom, second order linear 2.3 2.4 system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Original and modified Morse potentials . . . . . . . . . . . . . . . Examples of potential functions . . . . . . . . . . . . . . . . . . . 28 38 40 (a) Graph representation of a small world network on 20 nodes with in (c) the correspondent average path lengths. (b) Scale free network on 60 nodes with in (d) its node degree distribution. . . . 42 2.6 2.7 Vortex Formation . . . . . . . . . . . . . . . . . . . . . . . . . . . Distorted Formations . . . . . . . . . . . . . . . . . . . . . . . . . 48 49 3.1 The illustration shows how for k = hN/2i the swarm must be connected . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Swarm systems in a two-dimensional space relaxing to different 56 2.5 3.2 shapes as the number of agents N increases while the number of connections allowed per agent is held constant at k = 60. . . . . . viii 57 3.3 Dumbbell emergence in a three-dimensional case. (a) A 130 agent swarm in a dumbbell shape due to the number of connections per agent being limited to 65 and (b) the corresponding adjacency matrix obtained associating each node to one agent, after having them sorted from one end to the other of the dumbbell. Dots represent nonzero entries. . . . . . . . . . . . . . . . . . . . . . . . 3.4 3.5 3.6 3.7 58 Circulation definition . . . . . . . . . . . . . . . . . . . . . . . . . Adjacency matrices for spontaneously relaxed (a) and idealised connected (b) 60 agent swarm. . . . . . . . . . . . . . . . . . . . . 61 Adjacency matrices for idealised connected 61 agent swarm. . . . Comparison of Cheeger constant obtained by the analytic expres- 64 62 sion and the one obtained by making a swarm of agents relax through pairwise interactions for even number of agents. Average values for each potential and extreme values obtained over all the 3.8 numerical tests are reported as well. . . . . . . . . . . . . . . . . . Comparison of Cheeger constant obtained by the analytic expres- 66 sion and the one obtained by making a swarm of agents relax through pairwise interactions for odd number of agents. Average values for each potential and extreme values obtained over all the numerical tests are reported as well. Some of the mamimum values have been cut off the plot to make the curves more clear. . . . . . 4.1 4.2 4.3 4.4 4.5 4.6 67 Connections between groups of fully connected agents. . . . . . . Adjacency matrix for an ensemble of 25 agents. The nonzero entries are represented by dots. . . . . . . . . . . . . . . . . . . . . . 76 77 Adjacency matrix for an ensemble of 125 agents. The self-similarity of the matrix can be noticed. . . . . . . . . . . . . . . . . . . . . . 77 Node degrees as number of links belonging to each node. A selfsimilar scheme can be observed with nodes in central position being the most connected ones. In this scheme the maximum number of connections per node is 28. . . . . . . . . . . . . . . . . . . . . . . Homogeneous square lattice array configuration . . . . . . . . . . 5 agents arranged in a homogeneous formation due to all-to-all 78 79 potential with the same coefficients for all the agents . . . . . . . 82 4.7 4.8 Contours for the potential sensed by agent in (0,0) in the case the other agents have all the same value of the repulsive potential scale length Lr (i) and in the case one of the other agents has a repulsive scale distance Lr ′ < Lr . . . . . . . . . . . . . . . . . . . . . . . . 84 Cross pattern emerging by shrinking the repulsive potential scale length broadcasted for the agent in the centre. . . . . . . . . . . . 85 4.9 First three stages of growth of the Purina fractal array . . . . . . 4.10 Directivity plots for the first three stages of Purina fractal arrays shown in Figure 4.9, with an assumed direction of the main beam 90 towards broadside (θ = 0) . . . . . . . . . . . . . . . . . . . . . . 91 4.11 Groups of agents controlled at each step throughout the duty cycle. 95 4.12 Formation deployment in GEO. Snapshots taken at (i) t=0 s, (ii) t=60 s, (iii) t=1 hr, (iv) t=24 hr . . . . . . . . . . . . . . . . . . 97 4.13 Errors in design positioning after 1 day from release of the formation 98 5.1 5.2 Wheeled robots used for swarming test. Not to scale . . . . . . . 107 One of the robots used in the tests. . . . . . . . . . . . . . . . . . 108 5.3 5.4 Testbed architecture . . . . . . . . . . . . . . . . . . . . . . . . . 109 The testbed arena where tests are performed, with Vicon T160 cameras in red. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.5 A formation of 5 agents arranged in (a) pentagon shape due to all-to-all potential with the same coefficients for all the robots and (b) in a cross shape due to one agent being less “repulsive” with respect to the others. Attitude is not considered. Formation is defined just by the middle point of the wheel axle . . . . . . . . . 112 5.6 Trajectories followed by the 5 agents while arranging in a pentagon formation first (stars) and then in a cross formation (circles). . . . 118 5.7 A formation of 5 agents arranged in (a) pentagon formation and then (b) in cross formation. The pictures were extracted from the video taken by one of the Vicon cameras. . . . . . . . . . . . . . . 119 5.8 Trajectories followed by the 5 agents while arranging in a cross formation and rotating about its centre. . . . . . . . . . . . . . . 120 A formation of 5 agents arranged in (a) four point star and (b) 5.9 in a three point star that excludes an agent with a failure. The pictures were extracted from the video taken by one of the Vicon cameras. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.1 Eigenvalues of a negated Laplacian matrix according to the Ger- 6.2 shgorin circle theorem. . . . . . . . . . . . . . . . . . . . . . . . . 128 Sketch of the lattice network. Each node in the inner section has 3 outgoing and 3 incoming edges. Nodes at the left hand side have more outgoing than incoming edges, while nodes at the right hand side end have more incoming than outgoing edges. Note that an 6.3 outgoing edge implies that the node where the edge is originated “observes” the node where the edge ends. . . . . . . . . . . . . . . 145 Consensus in a lattice, performance comparison. . . . . . . . . . . 147 6.4 6.5 Consensus in a ring, performance comparison. . . . . . . . . . . . 148 Consensus in a Small World network, performance comparison. . . 149 6.6 6.7 6.8 Consensus in a random network, performance comparison. . . . . 150 Manoeuvring in three-dimensional space. . . . . . . . . . . . . . . 158 Error time history with respect to the desired heading. . . . . . . 159 6.9 Detrimental effect of inter-agent forces in manoeuvring. . . . . . . 160 6.10 Detrimental effect of inter-agent forces on the error time history with respect to the desired heading. . . . . . . . . . . . . . . . . . 161 List of Tables 2.1 2.2 Conditions for Lyapunov asymptotic stability. . . . . . . . . . . . Examples of connected and disconnected graphs. . . . . . . . . . . 29 32 4.1 Numerical values of coefficients used in numerical simulations. . . 96 5.1 Numerical values of coefficients used in numerical simulations re- 5.2 ferred to a cross formation with a arm of 700 mm. . . . . . . . . . 115 Results of the static relative positioning tests for a cross formation with a 700 mm arm. Tolerance and error are referred to the design inter-agent distance . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.1 Graph 2 example . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 6.2 6.3 Graph 1 example . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 First eigenvalues of the system matrix (L + D) where, L is the network graph Laplacian and D is a diagonal matrix whose nonzero entries are either the first left eigenvector (F. L. E.), a uniform vector or a vector numerically optimised to maximise the corre- 6.4 6.5 spondent eigenvalue. All of these have unitary L1 norm. . . . . . 151 Coefficients used in the APF for the three dimensional manoeuvring simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 Desired velocity vectors for the swarm to produce the changes in heading. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 xii Chapter 1 Introduction This chapter deals with the different aspects of aggregation of individuals. It separates the case of groups of individuals from structures composed of cooperating elements, and the case of spontaneous aggregation in biological groups from the emerging aggregation in engineering applications. It presents the concepts of swarm engineering and swarm intelligence, defines the concepts of emergence, illustrates some emerging behaviours and stresses their advantages and the way they can be possibly exploited in engineering or other scientific disciplines. Definitions The engineering systems which are the subject of this dissertation are, in turn, composed of engineering systems themselves. This characteristic makes the multiagent systems susceptible to confusion. Terms such as units are used to indicate parts of a more complex device, yet they may be autonomous on their own and, in most cases, they are engineered devices themselves. In order to avoid misunderstandings due to terminology as much as possible, a glossary of possibly misleading terms is provided as an aid. - The term part refers generically to elementary entities, or groups of them, which belong to a set of interest. - A system is a set to which several parts belong. - A component is a part of a system. Components and can be either functional and dependent on one another, or autonomous and independent. 1 CHAPTER 1. INTRODUCTION 2 - A unit is an individual entity that does not constitute a system on its own, at least not one of our interest, but is a part of a system composed of interacting components. - A multi-agent system is an engineering system composed of a number of autonomous, yet interacting, units. - An agent is a unit that is capable of autonomous behaviour. - The term vehicle is used in place of agent when the attention is put on the mobility characteristics of it. - The term swarm is used to indicate a multi-agent system characterised by relative motion of its agents. Despite its biological connotation, it is used to address either biological and non biological groups. - The term ensemble is used as synonym of group of units. Only if this units are autonomous, hence are agents, ensemble is a synonym of swarm. Throughout this dissertation the term ensemble does not have any statistical connotation. - The term architecture is used instead of system to stress the characteristic of being designed with particular emphasis on the connections of its parts and their interplay. - The term distributed indicates the characteristic of multi-agent systems of performing some tasks in absence of central coordination. - The term state indicates a set of time dependent parameters that characterise the units and, hence, the system. In the most common use through the dissertation it indicates position and velocity, but other parameters, such as internal energy, may be included in general. Obviously this glossary does not complete the set of terms that are used in the following, yet it provides some reference to which this work will remain coherent in its development. CHAPTER 1. INTRODUCTION 1.1 3 Towards Decentralised Architectures The idea of achieving challenging tasks through the cooperation of elementary units, rather than through a dedicated, centralised system, is relatively new in engineering. Natural science, however, has been looking at the dynamics arising from the interaction of a number of elementary parts since the mid-20th century [1]. This particular class of problems has been analysed through in the fields of chaos, complexity and collective behaviors within disciplines such as particle physics, mathematics, biology, sociology and many more. The behavior of such systems cannot be explained by looking at only the single components, instead the interplay amongst these single components must be considered [2]. While science is in general interested in modelling those systems belonging to the natural world, to get some characteristics and predict their behaviour, engineering tackles the problem from the opposite angle: it models and designs the system such as to get a desired behaviour. In this sense the design of complex systems starting from the autonomous units that form them is awkward. Biology and science in general have been fascinated by some collective behaviours and have studied them for over a century [3]. Conversely, in engineering, the idea of multiple, autonomous units designed to cooperate and to achieve a result through their interaction is more recent. The seminal work of John von Neumann about Cellular Automata in the pioneering era of computer science can be considered the first step towards the shaping of a new branch of engineering. Cellular automata (CA) are cells packed in arrays capable of processing information [4]. Each element, (automaton) is characterised by a finite number of states, specifically “ON” and “OFF”, and it switches between them depending on the information it receives from the neighbouring units. The switch happens synchronously at discrete time intervals; the result is the creation of evolving patterns in the array. The popular Conway’s Game of Life is simply an application of CA [5]. Cells in an array are turned on when there is a “birth” event, kept on if they “survive” and turned off when they “die” as Figure 1.1(a) illustrates. In the figure it can be seen how each cell is adjacent to 8 others in the array. At each time interval, also called generation, the “birth” of a new cell happens wherever an empty cell is adjacent to 3 existing ones, while any living cell dies whenever it is adjacent to less than 2, or more than 3 existing others. The patterns propagate from the initial state corresponding to the first column to successive states. As the rules for “birth”, “survive” and “die” events are given, then the outcome is determined univocally by the initial 4 CHAPTER 1. INTRODUCTION configuration. CA allowed for a rigorous mathematical description of their dynamics but present limits in the mainly theoretical framework in which they develop and their self contained dynamics. The information on the neighbour states is processed by the Automaton and returns its state which is in turn another piece of information for the other neighbouring CA. This limitation was pointed out and overcome, at least in a theoretical framework, by Gerardo Beni in 1988 [6]. Beni introduced the concept of “Cellular Robotic System” (CRS). CRS are different from CA as the units, although always arranged on a grid, are able to move and operate physically on the surrounding environment. In Figure 1.1(b) the grid arrangement is visible together with the capability of the CRS to operate on the environment by moving objects in the grid. The leap forward provided by CRS with respect to CA can be summarised in the possibility of processing information and matter as well. This is the characterising feature of robotic systems and, with its multiple agent architecture, the CRS can be considered the first example of “swarm robotics”. The term “swarm” in this context was first proposed at “Il Ciocco” conference while attempting to find a name that could refer to the CRS just presented at (a) (b) Figure 1.1: (a) The game of life as example of cellular automata as reported in [5]. (b) Beni’s original concept of cellular robotic systems as illustrated in its work of 1988[6]. The picture refers to the case of moving objects in the field. CHAPTER 1. INTRODUCTION 5 the same event in an appealing way. An immediate link to something popular such as the biological swarms was considered a key for the success of that new idea as reported by Beni himself [7]. Indeed, the concept of agents designed to achieve a global task through local interactions only has several similarities to biological swarms. The number of agents involved is too large for an approach that considers each of them singularly, and too small for applying statistical methods. Moreover, all the agents are quasi-identical and autonomous so that a decentralised control scheme emerges naturally. Finally, a consequence of the shift from Automata to Robotic Agents is the possibility of exploiting “stigmergy”. Stigmergy is the capability of an agent of passing and processing information by altering the surrounding environment. It is well know how ants exploit this by releasing pheromone along the trails from the nest to the food sources. This way the colony is able to trace the trail back to the food source [8]. CRS are agents with robotic characteristics, that is to say, they are able to process matter and hence to operate on the environment. Due to this, CRS are potentially able to exploit stigmergy. Although the name “Cellular Robots” was used by Fukuda [9] a few months before the seminal paper of Beni, yet Beni was the first to spot the potential of exploiting the collective behaviour and hence to start the swarm robotics. One year before the CRS was presented to the scientific community, Craig Reynolds, in a context of computer graphics, presented an algorithm to reproduce the flocking behaviour of birds [10]. The algorithm was based on three simple rules applied locally for each virtual bird, he called “boid ”. These are: Separation: each boid shall keep a minimum distance from any of its neighbours; Cohesion: each boid shall stay within a given distance from any of its neighbours; Alignment: each boid shall align its velocity to that of its neighbours. The neighbours of each boid were its flock mates within a given sensing distance. As simple as it was, Reynolds’ flock model, not only was able to reproduce qualitatively the behaviour of birds flocking, but also was a source of inspiration for the scientific community that in a few months would have been introduced to swarm robotics. Initially, the local interactions in groups of swarming robots were defined exploiting heuristic rules, in the same way as Reynolds’ rules operate. Nonetheless, the engineering of such a new concept could not remain abstract from the certainty of behaviour for long. By its own nature, engineering needs to rely on certain rules and its product must deliver what it is designed for, reliably and in safe CHAPTER 1. INTRODUCTION 6 conditions, regardless the level of autonomy it is provided with. In this respect a system where the components are designed not to directly deliver the task, while the totality does so by exhibiting a behaviour different from the one of the components, presents several hazards for which verification is needed. From the contrast between the level of autonomy of a swarm robotic system and the need for predictable behaviour of an engineered product, Alan Winfiled et al. developed the most modern concept of Swarm Engineering [11]. In their work Winfield et al. point out how verification of a swarm system, from the engineering point of view, should include verifying both the single component level and the group level behaviours for the liveliness and safety. These characteristics correspond to the property of exhibiting a desirable behaviour and avoiding any undesirable behaviour. Verification of the behaviour should then rely on mathematical proofs. Hence the “Direct Lyapunov Method” is identified as a way to prove stability of the control system at agent level and statistical methods to model the evolution of the swarm and hence prove that it is stable for each configuration it is likely to take. In a swarm where internal interactions are based on neighbouring, and where, by definition, internal local interactions produce effects on the state at global level, the evolution of the neighbouring must be considered ahead of any stability analysis. On the other hand, predicting all possible neighbouring states is not feasible from an engineering point of view for a system that processes information, decides and performs in real time. The system is then to be designed such that, whatever the relative states of the single agents are, it results in a stable and predictable behaviour that exhibits the characteristics of liveliness and safety. This is the starting point of this dissertation. 1.1.1 Biological Aspects Social aggregation is a characteristic found in many living organisms in nature. Wolves developed strategies for hunting in packs that have been found to be based upon stigmergy [12]. Even elementary organisms get benefit from aggregation: bacteria, that are normally found in colonies, benefit from this property in a number of situations such as foraging and increasing their resistance to antibiotics. In doing this bacteria are able to process and transfer information amongst them not only in a diffusive way, but establishing an efficient network of connections CHAPTER 1. INTRODUCTION 7 that also in some cases presents features of hierarchical organization [13]. Finally, popular behaviours that are worth mentioning are those exhibited by fireflies and ants. Fireflies synchronise their flashing in order to make the signal stronger for potential mating partners [14], while ants construct social interactions exploiting stigmergy to locate and optimise the foraging path [8]. 1.1.2 Technology Aspects Multi-agent systems are becoming more and more popular for intelligence and military applications, see for example [15], mainly because of their flexibility in adapting to fast changing scenarios and robustness against failures or attacks [16]. Integrated systems that include space, aerial and ground based segments have been suggested as well for search and rescue or military applications although these are far from being realised [17]. More advanced projects are ongoing for tasks such as assisting human operators, firefighters in the case studied in [18], or structural or machinery inspection as analysed in [19]. One of the few operative examples of multi-agent system is the Onyx parachute system [20] that is able to operate a number of parachutes to land multiple payloads launched from a cargo aircraft using swarming algorithms for collision avoidance and precision landing. Very popular projects, now in their test phases (see for example [21] and the SARTRE project [22]), concern the self-driving vehicles which can be considered as an application of multi-agent systems. These vehicles do not need to be designed to collaborate actively, but operate simply considering each other as non collaborative obstacles. In this sense even just the introduction of an autonomous vehicle in the normal traffic makes it comparable to a multi-agent system. 1.2 Self Organization, Group or Swarm Intelligence The definition of the concept of intelligence is beyond the scope of this thesis, nor is of interest of that part of literature that this work refers to. Beni [7] pointed out that one of the characterising aspects of an intelligent behaviour is producing an improbable outcome, something ordered, hardly achievable simply by chance. The other is the impossibility of predicting it, i.e. it must be something original, CHAPTER 1. INTRODUCTION 8 different from the output of any non-intelligent machine such as a printer, that will always output something more tidy than ink stains but yet it will always print what it is commanded to print. Nevertheless, according to Beni, an intelligent swarm can be formed out of non-intelligent machines. Strogatz [23] pointed out how ordered behaviours and in particular synchronization are emergent phenomena in nature for which just small amount or even no intelligence is needed. Mechanical coupled metronomes were observed to sync as well as people clapping in a theatre do, although not previously trained for that. Of course intelligence, in the most common meaning of the word, is not a characteristic of the metronomes. It can be argued that not only a certain kind of intelligence, often identified as swarm intelligence, generates emerging behaviour, but that the swarm intelligence can be considered an emergent characteristic itself. 1.2.1 Concept of Emergence The word “emergence” was used already in this introduction. It refers to the process that leads to self organization of a system in its macroscopic properties through micro-scale interactions. According to Goldstein, it can be defined as: “the arising of novel and coherent structures, patterns and properties during the process of self-organization in complex systems” [24]. The patterns arising in the Conway’s Game of Life and the corresponding behaviour in the Beni’s Cellular Robotic Systems are examples of emerging behaviours. Indeed, as for swarm intelligence, emergent behaviours do not need processing to take place. Snow crystals form out of interactions of water molecules and a number of configurations are possible depending on the thermodynamic characteristics of the surrounding environment [25], yet no water molecule does anything but acquiring a minimum energy state. Fish schools are often found forming swirling vortices as result of keeping the same relative positions and orientation with their local neighbours while swimming [26]. Moreover, Parrish et al. noted that particular configurations, such as whorls, funnel or tori, can quite often be found in nature as results of abiotic interactions as well; popular examples include galaxies, hurricanes and tornadoes. Some cases are illustrated in Figure 1.2.1. In many respects, the natural world relies on emergent behaviours and its dynamics, at different scales, resemble very closely a natural cellular automata. One of the most striking examples is the case of the self-regulation of gas exchange in plants, consequence of CHAPTER 1. INTRODUCTION 9 Figure 1.2: Examples of emergent patterns in nature: a satellite image of an hurricane in the Atlantic Ocean (NASA), a fish school forming a vortex (National Geographic) and a snowflake from [25]. a distributed calculation perfectly reproducible by a CA [27]. 1.3 Challenges in modelling, simulating, implementing group Behaviours The wide variety of group behaviours observable in nature and achievable in engineered swarming systems increases the difficulty of studying these systems with a rigorous scientific approach. On one hand the multitude pushes the analyst towards finding a common ground on which a critical analysis can be carried out, on the other the features of the swarming agents so deeply characterise the system that it is impossible to disregard them. The motion constraints, the sensing capabilities and the environmental conditions make, for instance, a flock of birds exhibit dynamics different from those of a school of fish; similarly a swarm of wheeled robots will have completely different dynamics from one of Unmanned Aerial Vehicles (UAVs). A common approach to modelling the basic features of most swarming system is considering them as particle systems (see for example [28–32]). In this case the agents are considered as physical particles moving in a force field that can be generated by the reciprocal interaction, the surrounding environment, or both. The evolution of the system is then described through the dynamical system theory. Considering agents as particles, statistical methods are applicable and particularly suitable for large swarms [33, 34]. In a very similar physical-like fashion, some researchers exploited entropy measurements to capture the emerging ordered state [35, 36] that has occasionally created some confusion within the literature as to the physical meaning of the entropy [35, 37]. These approaches often do not consider any constraints on the agents, yet they can be claimed to be valid modelling strategies for particular applications such CHAPTER 1. INTRODUCTION 10 as underwater or space vehicles [38], or to analyse emerging behaviors on their own (see for example [30, 32]). When the number of agents is high, modelling them as particles, all driven by the same force field, is unrealistic. While particles exchange forces, agents exchange information and produce forces as control actions. As different from forces, information need dedicated communication channels, which are limited in number. To account for that, the presence of a communication network is considered in this dissertation consistently with the relevant literature in the field of network modelling. In this framework graph theoretic methods are used to study networked systems through rigorous mathematical tools (see for example [39–42]). When peculiar aspects of the agents are considered, such as particular interaction networks or motion constraints, the modelling looses its universality. It is then important to find a balance between looking at a wide horizon and focussing on few specific cases. In this dissertation the principal methodologies found within the literature will be reviewed, however the approach taken goes in the direction of preserving the generality of the system rather than deeply characterise it. Laboratory tests are then used to confirm part of the theoretic outcomes. 1.4 From Particles Systems to Intelligent Cooperating Agents The term swarm is found in several contexts. It can refer to the most popular biological meaning or to particle systems animated by internal motion, or indicating a multi-agent engineering system. This variety of fields has produced a number of approaches to the swarm modelling and control problem, each capturing different aspects. While statistical methods apply better when swarms are composed of a very large number of units, they are often considered not effective from a design point of view. Here the various ways of approaching swarm behaviours are reviewed considering their significance in different fields. The number of assumptions and the level of detail that make a model more peculiar to a given phenomenon must be balanced against the extent of its applicability. CHAPTER 1. INTRODUCTION 11 Hence, modelling self propelled particles, able to move freely in the space, has more comprehensive features than a model that considers unit actuation capabilities and environmental constraints, but the latest returns more valuable results when a specific application is targeted. 1.4.1 Stochastic Approach Statistical physics has been extensively used to model the performance of mobile agents in large groups. Considering an agent as a self-propelled particle, it is possible to study the evolution of the particle system state on a statistic basis. Despite the appealing characteristics and the wide range of possibilities opened by the stochastic methods applied to swarm systems, in the present thesis these methodologies are not considered as a completely deterministic approach is instead preferred. The statistical physics approach is more valid as more particles are considered. Schweitzer [33] considers some emergent behaviours in particle swarms and relates them to system thermodynamic properties. This way, rather complex behaviour can be traced back and linked to some global properties of the system, such as energy, and their rate of change. The presence of noise is considered as well and the system state variables are described in terms of their statistical distribution. Ebeling et al [34] [43] consider the state of a particle swarm subject to thermal noise in terms of its variable distribution in statistical terms. The velocity distribution across the swarm and its mean square displacement is derived. In this, as well as in similar works, results from numerical simulations match the predictions made by means of statistical physics tools when the number of particles in the swarm is reasonably high to provide the data with statistical significance. Martinoli [44] considers a different approach by modelling the evolution of a number of coordinated robots with the task of gathering some targets. The robots’ neighbourhood is defined on a statistical base. A series of stochastic events are triggered, each depending on the success of the previous one. The first event consists in assigning each robot a random position which also corresponds to assigning a neighbourhood. This includes the target position and the positions of the other robots in the field, with respect to the first one. If in its position the CHAPTER 1. INTRODUCTION 12 robot finds the object to catch, then the second event is triggered, that is the gathering and so on. By showing good matching between the performance of real gathering robots and the outcome produced by the series of stochastic events, Martinoli argued that problems such as gathering were easily tractable on a statistical base. This is to say that considering the odds, from time to time, that a certain event useful towards the accomplishment of the gathering task happens, is sufficient to predict the success or the failure of the task and in which timescale. This approach rules out, in theory, some simulation campaigns based on physical modelling of robots and environment that usually require several repetitions to account for all possible variables. An original approach to control problems with a stochastic perspective has been considered by linking a certain behaviour of a robot in a group of agents to a probability distribution influenced by the external environment. That is, each agent has a probability distribution associated to all possible behaviour it can pick, with the expected outcome corresponding to the most appropriate behaviour in response to the external stimulus as in [45]. This way the possibility of an agent not picking the correct behaviour is accounted on statistical base. For a more extensive review of the stochastic methods specifically applied to robotics the reader can refer to [46] where extended information is provided. 1.4.2 Deterministic Approach An engineering approach to the swarm problems requires reliable and verifiable behaviour. This means that proving a given system exhibits some global behaviours, either described through the probability distribution of its states or through the probability certain events will occur, is not always considered satisfactory from the verification point of view. As it was already pointed out in Section 1.1, the system should respond in an intelligent way, that is, showing both liveliness and safety characteristics in every circumstance. For this reason many authors looked at a more deterministic approach, where the group behaviour can be studied both at global and at agent level. In particular, the group level is addressed by constructing an attractor for a system of differential equations, each of them describing the dynamics of a single agent. CHAPTER 1. INTRODUCTION 13 McInnes proposed the Artificial Potential Function (APF) method, which is extensively described in Section 2.6.1, in the context of space formation flying as an elegant way to coordinate the deployment of several spacecraft in a given formation, where the final arrangement is guaranteed through the minimization of an energy function [47]. Bennet [38, 48–50] exploited the APF method and the bifurcation of the equilibrium to produce a pattern transition between two formations in swarms of autonomous mobile agents. The approach chosen in this case is not confined to a particular vehicle architecture, although the applications to spacecraft and unmanned aerial vehicles (UAVs) are considered with real world issues coming into play such as actuator dynamics and saturation. Mabrouk [51–53] exploited provable vortex behaviour to drive a swarm of robotic agents capable to overcome some local minima problems, which are peculiar of the APF approach, by relying on their internal state. Using similar methodologies and a focus on space applications, Badawy [54, 55] implemented hyperbolic functions to control vehicles flying in formation and for on-orbit assembly. In a comparable way, but exploiting more geometrical rather than energetic considerations, Spears and Gordon proposed the Artificial Physics method for the control of a formation of autonomous robots [56, 57]. Leonard [58, 59] proposed an approach for the coordination of groups of vehicles based on Lie algebra and Lyapunov stability. The latter will be addressed in detail in Section 2.3.2 as it is used to draw some of the results of this thesis. These studies, together with other ones about motion synchronisation (see for example [60–62]) were fundamental for the deployment of a fleet of underwater vehicles in Monterey bay, California, USA, for ocean sampling purposes [63–65]. 1.4.3 Graph Theoretic Methods Graph theory is extensively used when dealing with multi-agent systems and swarms as a powerful tool for their modelling and analysis. Graph Theory studies the mathematical objects called graphs. A graph is a pair G{V, E} where V is a set of vertices and E is a set of edges. More rigorous definitions and property descriptions are provided in Chapter 2, while in the remaining of this section a general overview of the graph theory practical utility is presented. 14 CHAPTER 1. INTRODUCTION The first example of graph study in the modern sense can be found in the 1736 paper by Leonhard Euler on the seven bridges of Königsberg problem [66, 67]. Euler’s work focussed on the problem of crossing all seven bridges of the city of (a) (b) Figure 1.3: (a) Illustration for the problem of the seven bridges of Königsberg (source [66]). (b) Graph associated with the problem Königsberg in a single trip crossing each of them only once. This is represented in Figure 1.3 where the city areas are modelled as nodes of the graph and the bridges as edges. This makes the problem tractable in a rigorous way excluding from it the information about the path to follow within the city areas. A seminal contribution to the modern Graph Theory was given by Paul Erdos with his model of random graphs. These particular graphs are constructed from a group of nodes that are linked to each other in a random way. Depending on the number of nodes and edges, it is then possible to study the characteristics of the network from a statistical point of view [68–70]. The works by Erdos were of fundamental importance for the study of complex networks and, more refined, later models largely benefited from them. In particular the model by Barabasi CHAPTER 1. INTRODUCTION 15 and Albert [71] and the model proposed by Watts and Strogatz [72] departed from the complete randomness of edge structuring. Both the models are widely used in works concerning multi-agent systems as they provide robustness and short average paths between nodes compared to lattice networks respectively. Moreover both have been found to be representative of several biological and social systems and able to capture some complex behaviours emerging from these. In multi-agent systems graph theory provides a means to track the network of connections amongst the agents and to extract its properties, such as its connectivity, which determines the stability of the system and its response as it will be explained in Section 2.4. In particular it is possible to model a multi-agent system as a graph whose nodes are the single agents and whose edges are the connections between any pair of them. The distance between any two agents can then be accounted as the number of edges through which they are connected (provided that such a path exists) and hence it is possible to predict the capability of the system to express coherent behaviour and its dynamical properties. This kind of tracking of system characteristics has also been performed in experimental tests as in [73]. A line graph coupled with inter-agent potentials is at the base of the control scheme that produces the ”Tetunagi-oni ” game described by Yamaguchi and Arai [74]. In this a number of agents, linked by pairwise communication in a line graph, encircle a single agent not belonging to the line graph. Connected agents produce an always narrower ark closing around the isolated one. This behaviour emerges out of the connection scheme and the control laws that are enforced in it. Graph theory can be used, beside to track the interactions in a group of agents, as control means itself. Controllability in this sense was investigated in [75], while the problem of converging to consensus, that is to a common behaviour, in absence of external drives, is extensively covered in [76]. 1.4.4 Behavioral Rules and Algorithms A popular approach to control problems in robotics in general is the use of behavioural rules; this approach is commonly found also in swarm robotics. When behavioural rules are in place, agents are driven by instructions of the kind ”if... then...”. Although the results obtainable are remarkable in terms of coordination and effectiveness in task achievement, the performances of the robots cannot be described in rigorous mathematical terms, hence a proof of the system responding CHAPTER 1. INTRODUCTION 16 positively to the requirements of safety and liveliness is not obtainable. Reynolds’ work, described in Section 1.1, can be considered a precursor of behavioural rules applied to swarm robotics. Lewis and Bekey [77] proposed this behavioural control method to be applied to a colony of microrobots for surgical purposes and supported their claims through numerical simulation on a grid in a discrete time, which mimic the CA approach described in Section 1.1. Alternatively this kind of approach can be validated through extensive simulations and testing campaigns as in [78]. A mathematical validation of a behavioural algorithm is given in [79] although no proof is provided. Argumentations based on an approximated mathematical model and numerical simulations support the claim of the validity and effectiveness of the methodology proposed. 1.4.5 Other Approaches The examples mentioned are not exhaustive of the extensive literature dealing with multi-agent systems, yet they provide an idea of the general approach to this kind of problem from an engineering research perspective. More extended analyses of some of the methodologies and the approaches taken are provided later, when relevant to the topic covered. For a more general overview about swarm robotics literature, which is beyond the scopes of this work, the reader can refer to the reviews in [80–82]. 1.5 Aims and Objectives This work focuses on controlling a swarm of agents in analytically provable ways. This is achieved in this dissertation through the use of graph theoretic methods and artificial potential functions. The work aims to develop a new approach to distributed control of swarms considering the reduction of the communication links as a way to enhance the control rather than a limitation. This looks towards the practical implementation of autonomous swarms where the ideal all-to-all communication network is often infeasible. In particular this thesis investigates: ♦ The influence of limited communication on the onset of fragmentation in swarms of autonomous agents; ♦ The leveraging of limited communication to shape a formation of autonomous agents; CHAPTER 1. INTRODUCTION 17 ♦ The issues connected with the implementation of artificial potential function based methods to real robotic agents; ♦ The connection between the sensing network and the fast maneuvering performance in biological swarms together with the possibility of applying it to engineered agents; ♦ The exploitation of swarming systems to real world applications such as space-based telecommunication arrays, structural inspection and defence. Considering the real world implementation, the application of the theoretical control algorithms produced to real robots will also be investigated. A clear contribution to knowledge and towards the affirmation of swarm systems is presented here through meeting practical control problems with analytically provable methodologies. 1.6 Thesis Layout and Methodology This thesis presents a quite extended first part over the first two chapters aimed to introduce the current state of research in swarm engineering and control but also to introduce the different fields which are touched within it. In this first part the general theoretical background for the methodologies used later is presented, clarifying aspects such as emergence, provable behaviours and consensus. The results achieved are then presented from Chapter 3 to Chapter 6. Conclusions and future developments in the last chapter then contain the author’s observations on the contents presented. This chapter has outlined the complexity of the problems related to the interaction of many autonomous units. The challenges presented by an always expanding technology paradigm were presented noting how the balance between an autonomous and a foreseeable behaviour is to be searched by relying on analytically provable methodologies. This would allow multi-agent systems to be truly autonomous with the resolution of the technology gap that separates us from this new engineering. CHAPTER 1. INTRODUCTION 18 The approach based on considering swarming agents as particle systems is discussed. This will be used often in this dissertation but, as different from some particle physics works, the methodologies adopted here will not make use of stochastic techniques. Methodologies more widely applied in engineering are introduced as well as graph theory, which will be exploited extensively in this thesis. The approach used is then detailed in Chapter 2 which presents the mathematical tools that are exploited in this work. Emphasis is put on the application of control theory as well as on way these tools apply to swarm control. The control methodologies used are presented in detail and reference to appropriate literature is made marking the point from which this thesis departs to produce the original achievements stated in Section 1.5. In Chapter 3 the problem of cohesion in a swarm with an incomplete communication graph is considered. The first result of this thesis is drawn by reversing the problem from a graph theoretical perspective. While a swarm is required to be connected in order to be controllable, here conditions to have it connected at all time are drawn avoiding the problem of fragmentation and presenting a measure for the weakness of the swarms with respect to fragmentation. While Chapter 3 presents valid results with a wide theoretical applicability, but does not target any particular application, in Chapter 4 the design of the communication graph in a swarm is used to produce a space based architecture aimed to telecommunications. Exploiting the emergence of a recursive shape and the resonances that this produces when composed of radiating elements, the concept of a distributed fractal antenna is presented. This is done within a control theory framework aiming to produce a formation flying architecture. Chapter 5 considers the application of the control laws produced onto hardware. A coherent behaviour emerges out of the actions of autonomous wheeled robots controlled by the same means used in Chapters 3 and 4, proving their flexibility in being transferred onto real hardware and adapting to some specific challenges. Chapter 6 presents the results obtained in terms of consensus achievement and its consequences on the fast manoeuvring of swarms of mobile vehicles. A theoretical study is presented and supported by numerical simulations that point to CHAPTER 1. INTRODUCTION 19 the algebraic characteristics of the communication network to improve the swarm performance in terms of reactivity. Conclusions and possible future developments for the results presented are finally detailed in Chapter 7. 1.7 Papers Authored Giuliano Punzo, Derek J. Bennet, Malcolm Macdonald. Swarm shape manipulation through connection control. In: Towards Autonomous Robotic Systems 2010. Lecture Notes in Artificial Intelligence, 6856 (1). Springer. Giuliano Punzo, Derek J. Bennet, Malcolm Macdonald. Enhancing selfsimilar patterns by asymmetric artificial potential functions in partially connected swarms. In: 11th Conference Towards Autonomous Robotic Systems 2011. 2011-08-31 - 2011-09-02, Sheffield, UK. Lecture Notes in Artificial Intelligence, 6856 (1). Institution of Engineering and Technology. ISBN 978-3-642-23231-2. Giuliano Punzo Fractal patterns in fractionated spacecraft. In: 62nd International Astronautical Congress 2011, 2011-10-03 - 2011-10-07, Cape Town, South Africa. Giuliano Punzo, Derek J. Bennet, Malcolm Macdonald. A fractally fractionated spacecraft. In: 62nd International Astronautical Congress 2011, 2011-10-03 - 2011-10-07, Cape Town, South Africa. Gareth Pierce, Giuliano Punzo, Gordon Dobie, Rahul Summan, Charles N. MacLeod, Colin McInnes, James Biggs, Malcolm Macdonald, Derek J. Bennet. Reconfigurable robotic platforms for structural health monitoring. In: 6th European Workshop on Structural Health Monitoring & 1st European Conference of the Prognostics and Health Management (PHM) Society, 2012-07-03 - 2012-07-06, Dresden. Punzo, Giuliano, Gordon Dobie, Derek J. Bennet, Jonathan Jamieson, Malcolm Macdonald. Low-cost, multi-agent systems for planetary surface exploration. In: 63rd International Astronautical Congress, 2012-10-01 2012-10-05, Naples, Italy. CHAPTER 1. INTRODUCTION 20 Giuliano Punzo, Philippos Karagiannakis, Derek J. Bennet, Malcolm Macdonald, Stephan Weiss. Enabling and Exploiting Self similar formations. IEEE Transaction on Aerospace and Electronic Systems. 2013 in press. Giuliano Punzo, Jules Simo, Derek Bennet, Malcolm Macdonald. Characteristics of Swarms on the Edge of Fragmentation. In preparation for resubmission to Physical Review E. Philippos Karagiannakis, Stephan Weiss, Giuliano Punzo, Malcolm Macdonald, Jamie Bowman, Robert Stewart. Impact of Purina Fractal Array Geometry on Beamforming Performance and Complexity. 21st European Signal Processing Conference. 2013-09-09 - 2013-09-13, Marrakech, Morocco. Chapter 2 Control Problem, Theory and Tools In this chapter the tools for the design, modelling and control of multi-agent systems, that will be used in this dissertation, are illustrated. The two main theoretical pillars upon which the mathematical foundations of this thesis are developed are dynamical system theory and graph theory. These will be detailed in the remainder of this chapter. Moreover, these tools will be linked to modelling and control of swarms, linking them to different behaviours observable in nature and reproducible through engineering systems. 2.1 Dynamical system theory Dynamical system theory is at the base of the analytical modelling of multi-agent systems regardless of whether the agents considered are particles, biological or engineered units. It provides the tools to deal with N-agent systems regardless their size and the physical space considered. The system dynamics are modelled through differential equations, with partial derivatives in the most general case. Given a system, its state is identified by the value associated to a number of variables, namely the state variables which are associated with its degrees of freedom. To make this clear, let’s consider an N-particle system moving on a plane, hence a two dimensional space; consider furthermore that the particles undergo acceleration depending on their positions and speed. This dynamical system has 2 × N degrees of freedom as each particle can move in the plane and a state composed of 2 × 2 × N variables as each particle state is characterised by its position and 21 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 22 velocity, which are in turn composed of two components due to the planar motion. In this dissertation first and second order dynamical systems are considered. A first order dynamical system can be expressed as dx = ẋ = f (x, t) dt (2.1) where, the vector x is composed by the state variables of the system, hence is called the state vector of the system. The time derivative of a given variable or function is often indicated as the same variable dotted. This kind of notation is used as well through this dissertation. A second order dynamical system is expressed through a vectorial equation or a system of equations in the form d2 x = ẍ = f (x, ẋ, t) dt2 (2.2) where, the state vector is in general composed of the vector x and its derivative ẋ. Given a system of m differential equations, each of n-th order, it is possible to turn it into a system of n × m differential equations of first order with n × m variables by defining (n − 1) × m auxiliary variables and associating them to the time derivatives of the original state variables. For the case of a second order system this turns into ẋ = v (2.3) v̇ = f (x, t) (2.4) The time dependance is included in the above equations for sake of completeness although the systems studied in the following do not have any explicit dependance on time. Such systems are said “autonomous” as opposed to the ones where the time is an explicit variable and are called “non-autonomous”. 2.1.1 System Energy and Hamiltonian When multi-particle systems are used as tool to model multi-agent systems, thermodynamic analogies are considered. In such a case canonical coordinates are de- CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 23 fined for the phase space. These are the generalised coordinates q and generalised momenta p and are linked to the Hamiltonian function H by the Hamiltonian Equations ∂H ∂q ∂H q̇ = ∂p ṗ = − (2.5) (2.6) Hamiltonian function H is defined as H = p · q̇ − L(q, q̇) (2.7) where, (·) indicates the scalar product and L is the Lagrangian, which is directly linked to the energy of the system and defined as L= T −U (2.8) with T and U being the kinetic and the potential energy respectively. The energy approach above is particularly important for the development and stability verification of control methodologies such as the artificial potential functions. 2.2 Control Theory Control theory provides the framework and the methods to design control algorithms, where such algorithms are meant as the formulas and procedures to determine the control actions. Control actions are executed by devices able to process information and translate actions on the system to control [83]. In classical monolithic systems, a controller is in charge of determining control actions of all the degrees of freedom of the system. In distributed systems, the degrees of freedom of the systems are the sum of the degrees of freedom of all the agents. In this case the control is distributed as well, that is, each agent has a controller that operates on the agent’s state only, although it exploits the information of other agents and the surrounding environment as well. CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 24 In a classic closed loop control scheme, for a monolithic system, the readings of the sensors that define the state of the system are compared to a reference state or guidance law; the control action decided is a consequence of the difference between the actual state of the system and the reference one. In the distributed control scheme, the agent’s own state is not only checked against the reference but also the states of other agents contribute to define the action to be taken. Agents in the system can be seen as a part of the changing environment surrounding the agent to control [11]. The controller of that agent has to consider the other units when driving it, determining the evolution of the distributed system towards its final state. A comparison between the control of a monolithic system and that of a distributed system is illustrated in Figure 2.1 where, it is shown how the control actions driving a classic monolithic system are driven by the system state, the external environment sensed and the reference or guidance law, whereas in a distributed system the states of the other agents are used as well to define the control for any unit and hence for the whole system. In many cases the Guidance and Control functions are not distinguished from one another and referred to generically as “Control”. Indeed if the system does not aim to any specific state in an absolute sense, but instead it aims to achieve a target defined on the base of the agent relative states, then the control provides the guidance information at each time step and tracks them as well. A complementary case occurs when the actuators are not considered or are idealised. What then comes from the comparison of the reference with the actual state is the feedback the system receives. This is often used as a simplifying assumption as it does not take into account how the control action is produced. In refined modelling of engineered systems, the control input is taken by a low level controller that maps it into a driving signal for the actuators, which are peculiar to the system considered. 2.3 Stability Theory Checking the stability of dynamical systems means verifying the long time behaviour converges to either periodic orbits in the phase space or to a fixed point. The stability verification of the control system is carried out considering the high level control. What matters for the control stability is that the control action commanded is in the direction of decreasing the error without considering whether 25 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS (a) (b) Figure 2.1: (a) The control scheme of a monolithic system and (b) the control scheme of a distributed system. 26 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS the low level controller and the actuators are able to achieve it. In this dissertation the control of swarming systems is verified for both linear and nonlinear stability. Here some insights about stability verification are given. 2.3.1 Linear Stability Linear stability for dynamical systems is verified through considering the eigenvalues of the Jacobian matrix. This is defined as the matrix of the partial derivatives of the system with respect to its variables. For instance considering a first order system in m variables            dx1 dt dx2 dt = f1 (x1 , x2 ...xm ) = f1 (x1 , x2 ...xm ) ... ... dxm = fm (x1 , x2 ...xm ) dt (2.9) with x = (x1 , x2 ...xm ), the Jacobian matrix is defined as  ∂f 1 (x1 ,x2 ,...) ∂x  ∂f2 (x1 ,x1 2 ,...)  ∂(f1 , f2 , ...fm ) ∂x1 J= = ..  ∂(x1 , x2 , ...xm )  . ∂fn (x1 ,x2 ,...) ∂x1 ∂f1 (x1 ,x2 ,...) ∂x2 ∂f2 (x1 ,x2 ,...) ∂x2 .. . ... ... ... .. . ... ∂f1 (x1 ,x2 ,...) ∂xm ∂f2 (x1 ,x2 ,...)   ∂xm  .. . ∂fn (x1 ,x2 ,...) ∂xm  (2.10)   For linear systems the analysis of the Jacobian is sufficient to prove stability. Stability of nonlinear systems can be analysed too through linearisation but just to some extent and in local sense only, that is the results hold just in a neighbourhood of the equilibrium point. The linearisation makes the Jacobian a constant matrix. If the linearisation is performed about a fixed point x0 , the Jacobian matrix can be calculated in that point as J0 = ∂(f1 , f2 , ...fm ) ∂(x1 , x2 , ...xm ) x=x0 (2.11) The Hartman-Grobman theorem states that if all the eigenvalues have nonzero real part (the fixed point is said hyperbolic), then the trajectories of the nonlinear system can be mapped to trajectories of the local linearised system in a neighbourhood of the equilibrium. The nonlinear stability can then be determined locally through the spectrum of the Jacobian. In particular the following behaviours can be related to the eigenvalues as exemplified in Figure 2.2 for a CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 27 linear damped oscillator and to a linearised one: - if all the eigenvalues are purely imaginary, that is the real part is zero, a linear system exhibits undamped oscillations (Figure 2.2.a); In this case, nothing can be said about the stability of a nonlinear system. - if all the eigenvalues have negative real parts, and some have nonzero imaginary parts, the system (both linear and linearised) stabilises according to an exponential decay that envelopes the oscillating behaviour damping it out eventually (Figure 2.2.b); - if at least one eigenvalue has a positive real parts, and some have imaginary parts the system (both linear and linearised) will exhibit oscillations of exponentially increasing magnitude (Figure 2.2.c); - if all the eigenvalues are real and negative the system (both linear and linearised) is exponential stable, that is it stabilises to the equilibrium according to an exponential decay (Figure 2.2.d); - if all the eigenvalues are real and at least one is positive the system (both linear and linearised) is exponential unstable, that is it exhibits an unbounded motion according to an exponential increase (Figure 2.2.e). It can be concluded that, for linear systems, an asymptotically stable behaviour requires eigenvalues with negative real parts, whereas a stable behaviour simply excludes the presence of eigenvalues with positive real parts. Moreover if the system dynamics is governed by a matrix whose eigenvalues have negative real part, the system will exhibit a stable behaviour. This kind of matrix is also known as Hurwitz matrix. For linearised system, instead, nothing can be said if at least one eigenvalue has zero real part. 2.3.2 Nonlinear Stability - Lyapunov stability In this dissertation the stability of nonlinear systems is analysed using the Lyapunov method, although other criteria to assess nonlinear stability exist for which the reader is addressed to the literature dealing with nonlinear dynamical systems and their stability (see for example [84–86]). Lyapunov method can assess 28 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 4 1 1.5 0.5 1 3 x 10 0 −0.5 2 Amplitude Amplitude Amplitude 2.5 0.5 1.5 1 0.5 0 0 −1 0 10 time −0.5 0 20 (a) 10 time 20 −0.5 0 5 10 time (b) 15 20 (c) 8 2.5 0.8 Amplitude Amplitude 1 0.6 0.4 0.2 0 0 5 10 time (d) 15 20 x 10 2 1.5 1 0.5 0 0 5 10 time 15 20 (e) Figure 2.2: Time history for a single degree of freedom, second order linear system in case of (a) purely imaginary eigenvalues, (b) complex eigenvalues with negative real part, (c) complex eigenvalues with positive real part, (d) a single negative, real eigenvalue, (e) a single positive, real eigenvalue stability for nonlinear systems characterised by smooth dynamics, although it has also been extended to non-smooth dynamics as will be illustrated in Section 2.5. In his book of 1892 “The General Problem of Stability” [87] Alexandr Lyapunov defined a method to verify stability of nonlinear dynamical systems that does not require solving the system. Before describing how this works, some definitions must be provided; these are mainly given accordingly to [86]. First the equilibrium point has to be defined: Given a dynamical system as the one in Equation 2.1, with initial conditions x(t0 ) = x0 a point x∗ is said to be an equilibrium point of the system if f (x∗, t) = 0 for all t ≥ 0. Without loss of generalities and to comply with the most spread formulation in the literature, from this point on consider x∗ = 0, that is, the origin. Stability in the sense of Lyapunov prescribes that x = x∗ = 0 is a stable equilibrium point for the system in equation 2.1 with initial conditions x(t0 ) = x if for all t0 ≥ 0 and ǫ > 0, there exists δ(t0 , ǫ) such that |x0 | < δ(t0 , ǫ) → |x(t)| < ǫ ∀t > t0 (2.12) where x(t) is the solution of Equation 2.1 starting from x0 at t0 and ǫ and δ are scalars. The stability is called uniform if δ can be chosen independent of t0 . The origin CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 29 x = 0 is an asymptotically stable equilibrium point for the system in equation 2.1 with initial conditions x(t0 ) = x if it is a stable equilibrium point for the system and it is attractive, that is for all t0 ≥ 0, there exists δ(t0 ) such that |x0 | < δ → lim |x(t)| = 0 t→∞ (2.13) The origin x = 0 is a Globally Asymptotically Stable equilibrium point if it is stable and limt→∞ |x(t)| = 0 for all x0 ∈ Rn with n the dimension of the state space. This means that an equilibrium point is Globally Asymptotically Stable if the system will always evolve towards it for any set of initial conditions considered. The method then foresees the verification of the stability by finding a scalar function V (x), defined over the variables’ space of the system, that takes zero value at the equilibrium and is positive definite everywhere else; moreover its time derivative is zero at the equilibrium and negative definite elsewhere. The conditions are summarised in Table 2.1 Table 2.1: Conditions for Lyapunov asymptotic stability. V (x) > 0 V̇ (x) < 0 V (x) = 0 for x 6= x0 for x = x0 V̇ (x) = 0 The function V is called Lyapunov function and finding it is the challenge associated with the method as no a-priori indication is given about its definition. Usually the energy of the system can be taken as a candidate for the definition of a Lyapunov function. The stability can be global or just local depending on whether the Lyapunov function satisfies the conditions in table 2.1 respectively in a neighbourhood of x0 or everywhere. 2.4 Graph Theoretic Methods Graph theory was introduced in Section 1.4. In this Section the terminology, used later in this dissertation, is defined. Moreover, some algebraic graph prop- CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 30 erties are illustrated for cases that will be discussed later in this dissertation. The reader can refer to a number of textbooks for a deeper understanding of graph theory. The textbooks [88–91] have been used as reference. A graph G{V, E} is a mathematical object defined as the pair composed of a set of nodes or vertices V and a set of ordered pairs of vertices E that are called edges. • A graph is said to be directed if all the edges are ordered pairs of vertices. A directed graph is also named digraph. A graph is said to be undirected if all the edges are unordered pairs of vertices. A graph with ordered and unordered pairs of vertices is said to be mixed, although it can just be considered as a special case of directed graph. • Any two vertices i and j are said to be adjacent if there is an oriented edge from vertex i to vertex j. • A path is a subset of edges arranged in a sequence such that any two of them, which are consecutive in the sequence, connect adjacent vertices. • A node i of a graph is said to be reachable by another node j if a directed path from node j to node i exists. Node i is said to be globally reachable if it is reachable by any other node of the graph. • A graph is said to be connected if and only if it has at least one globally reachable node. • A graph is said to be strongly connected if there exists a directed path between any two nodes of it. • A directed graph is said to be weakly connected if it can be made connected by replacing its directed edges with undirected ones. • A directed graph is said to be complete if and only if any two distinct vertices of the graph are the end-points of an edge of the graph. • The degree of a vertex u in an undirected graph is the number of edges which include u. The in-degree (out-degree) of a vertex u in a directed graph is the number of edges entering (exiting from) u. A graph is said to be regular if all the vertices have the same degree. CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 31 In the following only “sensing” graphs are considered, that is the oriented edges in the graph point to the nodes which are “observed” and have origin in the nodes that “observe” or “sense” the first ones. The term degree can then be used instead of out-degree in the context of directed graphs. The connectivity properties of a graph are related to the algebraic properties of two fundamental matrices: the Adjacency Matrix and the Laplacian Matrix. The adjacency matrix of a graph G on N vertices denoted by A(G) is an N ×N matrix (square matrix of size N), having rows and columns labeled by the vertices of G, and ij th entry, aij , defined as ( aij = 1 ⇔ ui , uj ∈ E aij = 0 otherwise (2.14) where, ui , uj ∈ V. It is important to notice that the way the adjacency matrix is defined implies that if an agent i observes an agent j, then there is a nonzero entry at the j-th column of the i-th row of the adjacency matrix. The Laplacian matrix of a graph on N vertices can be defined on the base of its adjacency matrix as L= D−A where, D is a diagonal degree matrix, that is its i-th diagonal element is the out-degree of the node i, that is the sum along the rows of the adjacency matrix. It is immediate to verify that L is positive semidefinite and, because of the zero row sum, it always has at least one zero eigenvalue. It can be proved that a graph is connected if and only if the zero eigenvalue of the Laplacian matrix is simple, that is, it has multiplicity equal to 1. When a graph is disconnected it is common to refer to its connected subsets as components of the graph. When the graph is connected, hence one only component is present, it can be referred to as the giant component. Two examples of graphs are reported in Table 2.2 with their adjacency matrices and Laplacian matrices. It can be seen that an edge connecting two generic nodes i, j corresponds to a nonzero entry in the i − th row at the j − th column of the adjacency matrix. The Laplacian matrix is obtained from considering a negated adjacency matrix where the i − th entry along the main diagonal corresponds to the outdegree of node i. 32 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS Table 2.2: Examples of connected and disconnected graphs. ✲ ❞1 2❞ (a) 3 ❅ ■ ❅ ❅ ❅ ❅ ❄ ❅ ✲ ❞4 ❞ 2❞ ❅ ■ ❅ (b) 3 ❞✠ ✒ ❅ ❞5  0 1 A= 0 0 0 0 0 1  0 0  1 0 0 1 0 0  0 0 0 0 −1 2 −1 0   L= 0 0 1 −1 0 −1 0 1  ❞1 ❅ ❅ ❘❞4 ❅  0 0  A= 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1  0 0  0  0 0  0 0 0 0 0 0 0 0   0 0 0 L= 0 0 0 0 0 −1 −1 −1 −1  0 0  0  0 4 Example 1. In the first example, (a), the graph is connected as node 1 is globally reachable. If node 3 is removed the graph would be on just 3 nodes, but it would still be connected. If instead the edge 3,4 is removed, the 4 node graph would be disconnected as there would be no path from node 3 to node 1, that hence would not be globally reachable anymore. Example 2. In the second example, (b), the graph is disconnected, but it is also weakly connected. This graph is also known as “exploding star”. If all the existing edges are made non-directed, the graph becomes connected as node 5 becomes globally reachable. The same result can be obtained by simply reversing all the edges, such that these all point to node 5. The graph so produced is also known as “imploding star”, that, as different from the exploding star, is connected. Adjacency matrix, and Laplacian in turn, can be weighted, that is the entries can attain values other than just 0 and 1. The entries in this case would represent the strength of an interaction beside its existence between two nodes. The linearisation of the interactions amongst agents in a networked system leads to a matrix representation where the system reduces to its Laplacian, or to a matrix proportional to it. Consider for instance a generic networked system with nonlinear interactions. Consider that each node or agent is characterised by a one-dimensional state variable xi . The dynamics of a generic agent i then evolves as consequence of the interactions with the other agents it is connected into the network, as a function of the difference between their states, as CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS ẋi = X j f (xj − xi ) . 33 (2.15) Considering the system in equilibrium for its state corresponding to the vector ξ, said x̃i the state linearisation about ξ, this is x̃˙i = K  X j ∂f f (ξj − ξi ) + ∂xj xj =ξj (xj − ξj − (xi − ξi ))  (2.16) where, the sum extends to the K nodes node i is connected with. As ξ is an equilibrium point, Equation 2.16 reduces to x̃˙i = K  X ∂f j ∂xj xj =ξj (xj − ξj − (xi − ξi ))  (2.17) that, in vector form becomes x̃˙ = −LJ0T (x − ξ) (2.18) where, J0T is the Jacobian matrix referred to the equilibrium point as defined in Section 2.3.1 and the superscript indicates the transpose. The linear stability can then be studied considering the spectrum of the LJ0T matrix. 2.5 Nonsmooth Analysis Nonsmooth analysis is a subfield of nonlinear analysis. It studies non differentiable phenomena, framing them within the rigorous boundary of mathematics or, conversely, expanding the classic nonlinear analysis to non differentiable functions [92]. When considering a time changing graph, or intermittent interactions in a swarm of agents, discontinuities are found in the time history of the signal that the agents receive from the other units in the swarm. In the particular case of the APF method, the derivative of the potential corresponding to acceleration or velocity commanded to the agents cannot be defined within the classical analysis whenever a discontinuity is found in the right hand side of the equations. Without the appropriate tools is then impossible to assert that the system is not going to output an unbounded response at the discontinuity, that is, stability cannot be proved. In such cases nonsmooth analysis is used to bound the value CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 34 of the derivative to a set, which can then be used to prove system stability. The first important contribution to the solution of this kind of problems was given by Filippov [93] who later provided the analytic tools to study equilibria in this framework through using smooth Lyapunov functions [94]. The extension to nonsmooth Lyapunov functions was provided by Shevitz and Paden [95]. In this dissertation nonsmooth analysis is considered only as a tool for future developments, however, it is important to address its role in the study of swarms and multi-agent problems. In recent papers Tanner [96, 97] studied the performance of groups of mobile agents linked in a network based on proximity. As the agents are able to interact with neighbours only within a given radius and the relative motion of the agents changes the local neighbour in a continuous fashion, the network of interaction changes dynamically. In this case nonlinear attractive-repulsive interactions govern the motion of the agents that can be proved to align their velocity along a common direction, hence producing a stable behaviour, by means of the nonsmooth analysis. This is pursued by proving that the generalised derivative at the discontinuity is such that the system does not destabilise because of the continuous jumps the agents experience in the control input. The requirement such systems have to be able to prove their stability is that the graph representing the interaction network remains connected over all the configurations the network takes throughout the switches. This hypothesis is fundamental to guarantee that the system evolves to the consensus. It is the case to note that being continuously connected is a strong assumption as the connectivity of the graph of interaction cannot be guaranteed by simply relying on the interactions within the local neighbours. 2.6 Swarm Modelling and Control Mathematical handling of swarming problems requires the definition of interactions between the agents from an analytical point of view. This is done for the interactions in terms of signals and actions exchanged between pairs of communicating agents, and the network that spans the swarm, which is mainly responsible for its cohesion. In this section, the method of the artificial potential functions CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 35 is explained with reference to the Morse potential. Also there is the need to consider the role of the network of interactions in the swarm and the control issues depending on the choice of both the network and the Artificial Potential Functions. 2.6.1 Artificial Potential Functions Artificial Potential Functions (APFs) are a popular control means for multi-agent systems, although it is also quite often used to produce a time variant reference guide rather than a control output. In APF methods, given a system modelled as a dynamical system ẋ = f (x, t), a scalar potential is defined as a function of the system state, namely U = U(x). The control output ẋc is then defined by the negative gradient of the potential as ẋc (t) = −∇U(x) (2.19) where ∇ indicates the gradient. For a generic function ϕ(χ1 , χ2 ...χn ) in the generic n variables χ1 , χ2 ...χn the gradient is defined as the vector field defined by the partial derivatives of ϕ, that is ∇(ϕ) = ∂ϕ ∂ϕ υ1 + υ2 + ... ∂χ1 ∂χ2 + ∂ϕ υn ∂χn (2.20) with υ1 , υ2 , ...υn being the unit vectors forming the canonical base of the ndimensional space. The control output is zero in case of stationary points in the potential. The system is therefore driven towards the minimization of the potential that corresponds to the target state. In order to ensure stability, the potential function has to be convex about the equilibrium point. The APFs can either be dependent upon the state of the system in absolute sense, that is, the position and velocity of each agent with respect to a given pair of parameters, or in relative sense, that is the potential for each agent is defined on the base of its state compared to that of one or some of the others. In the second case APFs are sometimes referred to as based on pairwise interactions, or more briefly pairwise APFs, which is the case of the systems analysed in this dissertation. APFs can be used to define the guidance law for an engineering system. In the case of multi-agent mobile systems, they can define the velocity field according to CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 36 Equation 2.19. The control system has to track this velocity field by producing the accelerations through the actuators. When the APFs are used to model dynamical systems this can be done as ẍ(t) = ẍc (t) = −∇U(x) . (2.21) Equation 2.21 can then be integrated twice to get the time history for a given set of initial conditions on position and velocity. In this second case the APFs can be used to model particle systems. When used to control engineered agents some issues arise. By their own nature, physical systems which undergo the action of the potential only are conservative, that is the energy provided by the potential field to the system is transferred to the system as kinetic energy and returned to the potential field when the system slows down without leakages. Hence the APFs alone would not be suitable in general to drive a system from an initial state to a final state characterised by a different level of total energy. For this reason the APFs are often complimented by an orientation function that, in the general case, pumps in or out virtual energy to the system. The Hamiltonian representation of the system, with the notation introduced in Section 2.1.1, is then ṗ = − ∂H ∂H − g(H) ∂q ∂p (2.22) where, g(H) is some function in general able to reduce or increase the energy level of the system. If g(H) is an always positive function, then the system will always reduce its energy level. Conversely, if it is an always negative function, the system will increase its energy level. This can be easily seen by computing the time derivative of the Hamiltonian function as X  ∂H 2 dH = −g(H) dt ∂pi i (2.23) In the simplest case of viscous-like effects g(H) is often taken constant and the system will eventually damp out all its kinetic energy relaxing to a minimum of its potential energy. Morse Potential Morse potential is named after Philip M. Morse who defined it in 1929 [98] for the purpose of modelling the vibrational structure of diatomic molecules. It is a CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 37 scalar potential composed of the sum of two exponential functions and its first formulation included explicitly the position and the value of the minimum for the potential associated to the bond between the atoms of a diatomic molecule. It is E(r) = De−2a(r−r0 ) − 2De−a(r−r0 ) (2.24) where, D and a are constants that define the magnitude of the potential as well as its attraction and repulsion range. The potential has a finite value at zero distance and vanishes asymptotically for r → ∞. In this formulation the potential attains a unique minimum of −D at a distance r = r0 . This is illustrated in Figure 2.3.a where, the attractive and repulsive contributions to the potential are shown as well as their sum. When modelling particle and swarming systems, a more convenient formulation is commonly used for pairwise interactions between two generic particles or agents i and j, see for example [28, 32]. This is UijM orse = −C a e− |xi −xj | La + C r e− |xi −xj | Lr (2.25) where, C a and C r are constants that tune the magnitude of the potential while La and Lr regulate the range of attraction and repulsion. This formulation uses vector notation for the positions of two generic agents i and j with | · | indicating here and thereafter the euclidean norm when applied to vectors, the determinant when applied to matrices and the absolute value when applied to scalars. The attractive and repulsive parts are clearly identified by a the superscripts a and r on the coefficients. This version of the potential presents a unique minimum at La Lr ln xi − xj = r L − La  C a Lr C r La  (2.26) where the potential attains the value of M orse Umin = −C a  C a Lr C r La r  LaL−L r + Cr  C a Lr C r La a  LaL−L r (2.27) Figure 2.3.b shows the attractive and repulsive functions that together form the Morse potential. These can be easily compared to the ones in Figure 2.3.a which refer to the original formulation of the Morse potential. 38 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 10 1 Morse potential Attractive potential Repulsive potential 0.5 U E(r) 5 0 0 −0.5 −5 −10 0 Morse potential Attractive potential Repulsive potential 1 2 r 3 4 (a) −1 0 1 2 xi−xj 3 4 (b) Figure 2.3: Original and modified Morse potentials. (a) A plot of the Morse potential as it was published in [98], according to Equation 2.24 with D = a = r0 = 1. The two exponential functions that shape the potential are shown as well and (b) a plot of the Morse potential according to Equation 2.25 with C a = C r = La = 1 and Lr = 0.5. Lennard-Jones and other Potentials The use of APFs is not limited to the Morse potential. In particular, when studying molecular dynamics a popular potential is that of Lennard-Jones. It was originally introduced to model the interaction between neutral atoms in gas molecules in correlation with the study of the relationship between the gas viscosity and the temperature [99]. Distinct from the Morse potential, the Lennard-Jones potential presents a singularity at the origin, that is, it produces an unbound force over particles at very close distance. The Lennard Jones potential sensed by a particle i in the vicinity of a particle j is λ2 λ1 − (2.28) UijLJ = n |xi − xj | |xi − xj |m where, λ1 and λ2 are free parameters that can be used to set a desired separation distance between particles, n and m are exponents that, in the original formulation and quite often in literature, take value 6 and 12 respectively; however, these can be tuned to better design the interactions between pairs of agents. In literature it is easy to find examples of potential functions other than Morse and Lennard-Jones. Badawy, McInnes and Bennet (see for example [49, 55]) exploited the property of hyperbolic potentials that have almost constant gradient far from CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 39 equilibrium. The hyperbolic potential U hyp is described by the equation Uijhyp = Ch [(|xi − xj | − d)2 + 1]0.5 (2.29) where, d is the desired distance between the agents and Ch is a scaling parameter. This potential produces a control law where saturation problems can be more easily tackled by scaling the actuator outputs between zero (at equilibrium) and the maximum gradient achieved by the potential. Moreover the quasi-linearity of the potential far from equilibrium allows for superposition of the effects easily tractable in analytic modelling. Some authors, see for example [59, 97, 100], proposed artificial potentials that present a discontinuity in the gradient at a designed cut-off distance, that is a distance within which two agents interact, while they do not anymore if further apart. In this case the control input experiences a jump when two agents enter each-others range. This should be accounted for in any stability analysis and tackled through non-smooth analysis. Three examples of potential functions quoted in this sections are reported in Figure 2.4 where some characteristics can be noted. The singularities of the Lennard-Jones potential visible in the first two panels are responsible for unbounded control actions at close distance. The linear trend of the hyperbolic potential far from equilibrium corresponds to a bounded derivative, which translates into a constant acceleration commanded to the agents when far from their design arrangement. Finally, in the rightmost panels, the non-smooth blending of the potential at a cut-off distance produces a clear discontinuity in the derivative, that translates to a sudden change from a continuous to a null control action. Lastly, it is worth noting how, even in cases where the potential does not present a cut-off distance, some of the functions used in literature present a vanishing gradient at large distances making the interaction faint and, practically, ineffective beyond a given radius. This is the case of Lennard-Jones potentials and Morse potential for instance which are often claimed to affect the agents relative neighbouring only. 2.6.2 Connection Networks Beside the interactions between pairs of agents, the other fundamental constituent of distributed control is the network of connections which link all the agents. All- 40 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 3 60 2.5 40 −0.4 U(r) −0.2 U(r) U(r) 0 2 20 1.5 0 −0.6 0 1 r 2 1 0 3 −3 3 1 2 r 3 4 −20 0 5 1 r 2 3 −3 x 10 1 x 10 5 4 0.5 1 3 0 f(r) 2 f(r) f(r) 2 0 1 0 −0.5 −1 −1 −2 −2 0 1 (a) r 2 3 −1 0 1 2 r 3 4 5 −3 0 0.5 (b) 1 1.5 2 2.5 3 r (c) Figure 2.4: Examples of potential functions in the first line and their gradient, below, in one dimension. (a) Lennard-Jones potential; (b) Hyperbolic potential; (c) Potential presenting a gradient discontinuity as found in [100] to-all connections are commonly used in works motivated by the investigation of particle systems. Attraction and repulsion are experienced by each particle because of all the other particles in the system: the corresponding graph is hence complete and equilibrium conditions for a single particle depends on the states of all the others. In works considering the engineering aspects of multi-agent systems, the infeasibility of all-to-all communications is often considered. Proper networks are hence introduced which can be classified in two main categories based on their graph, that is Symmetric and Asymmetric. In the symmetric case an edge between two generic nodes i and j implies the existence of another edge between node j and node i that is oriented in the opposite direction. In this case the graph is non-directed and its adjacency matrix is symmetric. The symmetry of the network does not imply symmetry of interactions. In a symmetric network the existence of a generic non-oriented edge i − j does not imply that the action of node i on node j is the same of the action produced by node j over node i. On the other hand, non-symmetric networks exclude the possibility of symmetric interactions as in this case the corresponding graph is directed. In the case of agents interacting only in case of proximity, that is within a certain distance, these are often modelled considering pairwise forces which are null beyond a given range; the network arises then as consequence of this. If the range CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 41 is the same for all the agents, the network that is formed is symmetric but there is no guarantee that it is or stays connected as the swarm evolves. Particular networks have been proposed to link autonomous agents in coherent swarms while enhancing the overall performances by exploiting some of the network characteristics. In particular the models of Scale Free networks by Barabasi and Albert [71] and the Small World network by Watts and Strogatz [72], already introduced in Section 1.4.3, deserve further discussion for the attention they received in the literature. Scale free networks are formed by connecting nodes to an initial cluster, with each new node having a probability of being connected to an existing node proportional to the number of connections the latter has. The network emerging out of this process has a node degree distribution (number of nodes featuring a given number of connections) that presents a small portion of nodes, also named hubs, with high number of connections while most nodes feature just few connections. This arrangement produces advantages in terms of controllability as controlling just the hubs it is possible to control the whole network [101]. The benefits of small world networks in multi-agent systems have been stressed with an emphasis comparable to scale free networks. Small world networks are developed by rewiring randomly a small amount of links in a regular lattice. This is a network where each node is linked to a number of neighbouring nodes and neighbouring property is determined by the length of the path between the nodes in the graph. By rewiring a small number of them, shortcuts are created so that two nodes separated by a long path in the lattice find themselves much closer after the rewiring [72]. The creation of a small world network, out of a regular lattice, enhances faster information routing in the network. This has been exploited in a number of applications which feature the presence of a network, including the controlled deployment of sensor networks and consensus problems, see for example [102–104]. In Figure 2.5 the appearance and main features of small world and scale free networks are illustrated. While the network can significantly change the performance of a multi-agent system, its connectivity is a critical issue for the stability of the distributed control. This mostly affects networks based on proximity, including the case in which interactions present a cut-off distance, and 42 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS (a) (b) 3 25 20 2 Nr of Nodes Average Path Length 2.5 1.5 1 10 5 0.5 0 0 15 5 10 Node (c) 15 20 0 1 2 3 4 Nr of Edges 5 6 (d) Figure 2.5: (a) Graph representation of a small world network on 20 nodes. (b) Graph representation of a scale free network on 60 nodes. (c) The comparison between the average path lengths for the small world network, in red, and the ones of regular lattice, in yellow, show the characteristics of small world networks of lowering the average path length of each node with respect to any other in comparison with a regular lattice with the same number of edges. (d) The node degree distribution of the scale free network in (b) is typical of this kind of network as characterised by few nodes with many edges while most of nodes have just few edges. CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 43 dynamically changing networks in general whenever a control on the connectivity is not considered. The chance of having isolated components in the graph is often ruled out by hypothesis in literature dealing with stability of distributed control, but it is a practical issue to account for when actual designs are considered. 2.6.3 The Combined Effect of APFs and Connection Network The dynamics of swarms are influenced by both the connection network and the pairwise relations that govern inter-agent interactions in a different way. They are both responsible for system failures either for what concerns fragmentation of the swarm, or instabilities in the control respectively. On the other hand, the use of the connection network and the control law operating along its links can boost the performance of a multi-agent system. Badawy proposed the exploitation of the network to route the action of the artificial potential functions towards the agents to joint in a self-assembling structure by modifying the adjacency matrix of the graph [54]. This way a number of structural configurations are possible by changing the adjacency matrix and hence defining which are the agents to lock together. An alternative was taken by the Dynamical Control System Laboratory in Princeton. This considers the concept of tensegrity structures [100, 105] for shaping arbitrary formations of vehicles. Tensegrity structures are composed of a network of linking elements that meet in the nodes of the structure. The position of the nodes in the structure is uniquely determined by the force in the structural elements connecting them. In turn this force is determined for each element depending on a rigidity constant and on the base of its length. Depending on the difference between the rest length and the actual length of the elements, these ones can either be “strings”, if they exert a pulling force, or “struts” if they exert a pushing force. This way a multi-agent system can be viewed as a set of masses interconnected by springs and its equilibrium in a given formation can be found by assigning “spring constants” accordingly. Also in this case there is no need to connect all the agents in a complete graph and the stability can be analysed considering the energy of the system. Indeed the symmetry of the network and the interactions in it allow for stability proof following a Lyapunov approach and CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 44 relying on the energy of the system as shown in Sections 2.3.2 and 2.6.1. Tensegrity formations, however, exploit the precise tuning of the virtual forces amongst the agents in the formation. This requires either off-line computations before the deployment or distributed, heavy, on-line computations. Moreover, because of the unique choice of the spring constants for a given equilibrium, the failure of one agent would lead to sensible distortions in the formation. The limitations highlighted by tensegrity structures are overcome in this dissertation by exploiting some asymmetries in the interactions amongst agents. Moreover the hierarchical arrangement of the agents discussed in Chapter 4 returns more flexibility in the agent arrangements, as it will be shown later on. 2.7 Schooling, Flocking and Swarming Schooling, flocking and swarming are three emergent group behaviours found in nature. Although it is recognised that some individuals are more influential than others, there is no superior intelligence that coordinates the motion in animal groups. Since Beni [6], the engineering systems designed taking inspiration out of such phenomena have been named Swarms, although the biological connotation of swarming originally only referred to the behaviour observed in bee colonies when a new queen leaves the hive to form a new one. In this section these terms will be explained, considering their importance for modelling animal behaviours that are inspirational for engineering design of multi-agent systems. The conditions that distinguish a schooling behaviour from a flocking one in terms of engineering design parameters are extrapolated to then be used to model non-biological systems. Importance is given to the problem of cohesion that, although natural in biological system, might result of difficult implementation when the agents rely only on a limited number of communication channels. 2.7.1 Biological Studies Schooling refers to be behaviour often found in large groups of fish which swim keeping fixed relative positions and orientations. Flocking instead refers to a group motion in which individuals, although remaining in close proximity, do not keep the same positions with respect to one another; this is often associated to avian groups. Both these behaviours act to better protect the group from CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 45 predators by increasing the capabilities to quickly become aware of an attack and confuse the predators, making it difficult to isolate a single target. This way animals, which are smaller and weaker than the predator on their own, increase their possibility of surviving. In schools and in flocks, individuals keep track of a fixed number of nearest neighbours [106, 107]. This is understandable as senses are limited as well as the amount of resources individuals can invest in social behaviours. It was claimed, for instance, that the number of neighbours considered by each bird in a flock is on average between six and seven and this is not dependent upon the distance between the flock members [106]. This result was found through the analysis of footage of large starling (Sturnus vulgaris) flocks (see [108, 109]). The same set of videos provided evidence for a number of behavioural models and a way to validate their numerical implementation. It was found that the velocities of the birds in a flock are correlated in a scale-free fashion, that is the correlation does not depend on the size of the flock [110]. Again it was possible to obtain important data about the reaction time and the manoeuvre speed of starlings under predator attack [111], providing evidence in support of the Trafalgar effect phenomenon [112] and the Chorus line hypothesis [113]. These phenomena will be better detailed in Chapter 6 together with the importance for their engineering spin off. Schooling behaviour has also received noticeable attention, see for example [114, 115]. Biologists attentively studied the mechanism producing the emergence of tidy patterns. They found maintenance of alignment and lateral distance are mainly due to the sensing organs schooling fish are provided with. In particular, tidiness in schooling is possible due to the presence of the so-called “lateral line”, a string of pressure sensors that, operating in differential fashion, allow the sensing of movement. This is used independently of visual sensing, helping towards schooling and fast manoeuvring in response to lateral threats [116–119]. As mentioned at the beginning of this section, swarming has a meaning associated to the animal world as well as flocking and schooling. When swarming bees move in search of a new hive location, the group dynamics is more similar to flocking than to schooling being the cohesion maintained without necessarily CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 46 keeping alignment and relative separation. Several studies have been conducted to investigate the emergent behaviours that arise in colonies and the hierarchical structure that appears to be a result of genetic development of individuals. Individuals swarming in search of a nest site are bigger than their swarm mates and communicate with visual signals produced by an oscillatory movement. The final decision for the nest site becomes a “democratic” process in which all the swarm gives credit to the dancing scouts in the site proximity [120, 121]. The swarming behaviour is far more complex than a chaotic flight in close proximity and bees are inspirational for engineers in the design of complex systems [122]. 2.7.2 Crystalline Patterns The control of relative positions within a set of agents is at the base of advances in the technology for autonomous formation flights and self-driving cars. In schooling behaviour as well as in formation flying of some species, such as geese or ducks, keeping fixed positions in a given way is useful towards the reduction of the energy the individuals have to supply to the system. Engineering swarming systems get advantages from arrangement in fixed relative positions for a number of cooperative tasks. For example, surveying missions or, more generally, remote sensing tasks can be carried out by several autonomous agents in place of just one if these keep fixed relative positions. At the same time agents that keep fixed relative positions can rely on the navigation capabilities of just some of them in the position of leaders, leaving the followers more resources to utilise in mission related tasks. The task allocation is something observed in the natural world as well, with migratory events generally led by few informed individuals in the group. Distributed control of swarms in crystalline formation can be achieved by defining a virtual energy function depending on relative positions that attains a minimum in correspondence of the desired distance between any two agents. This can be done for both first order controllers and second order ones, that is, either producing the desired velocity or acceleration components as output of the controller. From an energy point of view, a swarm whose agents stabilise in fixed relative positions achieves a minimum, that is a static, stable equilibrium, that corresponds to a given level of energy. This is in general different from the one the swarm is provided with before the control operates. Controllers based on APFs leverage the level of total energy by adding or subtracting some virtual potential CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 47 energy when an unwanted quantity of real kinetic energy is respectively missing or over-abundant into the swarm. In order to achieve a crystalline pattern the controller has to damp out the velocity components that tend to make agents drift away from a given pattern or collide against each other. Damping out the kinetic energy does not imply achieving a static swarm, which is cancelling the motion of the group while staying in a pattern. In [32] self propelled particles arrange in a rotating vortex driven by a pairwise potential and a steering function that produce the equations of motion dxi = vi dt dvi m = −∇Ui − ui dt (2.30) , (2.31) for the i−th particle. In Equation 2.31 U is a Morse potential of the kind reported in Equation 2.25 while ui is a steering function in the form ui = Co X j (vij · x̂ij )exp( |xij | )x̂ij lo (2.32) where, {·} represent the scalar product, x̂ represent the unitary vector in a given direction, Co and lo are scalar constants and the double subscript indicates the difference between the variables associated with the subscripts i and j. Under the effect of the above potential and steering functions, particles then arrange in a crystalline formation rotating about its centre. This is shown in Figure 2.6 where the time history of the swarm’s total energy is reported as well and it can be seen that it is an always decreasing function as expected for this particular model. Excluding external drives or agents inputting energy into the system, in case of fixed relative positions the total energy stabilises to a minimum as all kinetic components are damped completely or in part. Results quoted so far in this section refer to systems where agents are either able to track the positions of all the swarm members or, in case of particles, they are subject to the potential generated by all the other particles in the swarm. In this framework it is possible to deal with the system in terms of total energy and global characteristics such as linear and angular momentum with a physics approach. When the number of connections each agent can keep, that is the number of agents it can keep track of, is limited, several issues arise. Figure 2.7 depicts the arrangement obtained 48 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS −330 0.4 −340 0.3 Total Energy 0.2 0.1 0 −0.1 −0.2 −350 −360 −370 −380 −0.3 −0.4 −0.5 −390 −0.4 −0.3 −0.2 −0.1 0 (a) 0.1 0.2 0.3 0.4 0.5 5 10 Time 15 20 (b) Figure 2.6: Vortex formation. (a) A crystalline arrangement emerging out of the interactions of 40 agents interacting with APFs and (b) the time history of system total energy, which is an always decreasing function as per [32]. with the same control law that produced the results in [32] previously shown, but, in this case, the number of connections per agent is limited and the connections are done producing a small world network and a scale free one (see Section 2.6.2). Clear distortions in the pattern are visible: these are due to the missing communication links between agents which are physically close but do not have awareness of each other. In Chapter 4 reciprocal and symmetric interactions will be shown to be sufficient to keep on considering the system in a physics oriented framework, with the possibility of achieving useful pattern for some particular space applications. The physics approach fails, or has to be heavily corrected, when the connections are not static, that is the case, for instance, of distance dependent interactions. It is still possible to prove that, as long as the system is cohesive, the energy will eventually decrease and stabilise to a minimum level, as per [96, 97], but this requires a nonsmooth analysis approach. To avoid this, the reciprocal influence can be considered a continuous and vanishing function of the distance, rather than a discontinuous one presenting a cut-off distance. This was done, for example, in [123]. Finally, the achievement of fixed relative positions can be put in relation with the organizational entropy [36]. It was observed that organizational entropy undergoes a phase transition when the free parameters, governing the system’s Hamiltonian, are tuned to pass from a chaotic behaviour to a stabilised one characterised 49 CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 0.8 0.8 0.6 0.6 0.4 0.4 0.2 Y Y 0.2 0 0 −0.2 −0.2 −0.4 −0.4 −0.6 −0.8 −0.6 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 −1 −0.5 0 X (a) 0.5 1 X (b) Figure 2.7: Distorted formations. The vortex arrangement produced in [32] with 40 agents either (a) fails to form or (b) presents dramatic distortions. In (a) a scale free network is implemented in which agents connect in turn to the ones which are already connected. In (b) the connections are restricted for each agent to its 8 nearest neighbours and some of this are rewired randomly. The probability of rewiring is set to 0.1. by the agents in fixed relative positions. Similar results were observed by Vicsek et al. [30] who also accounted for the role of the noise in the transition between more equilibrium configurations. 2.7.3 Continuous Swarming Systems In Section 2.7.2 the parallels between the natural behaviour of schooling and the motion of vehicles, keeping fixed relative positions, was presented. Another recurrent parallel with the natural world is the one between the flocking and swarming behaviours in birds and bees and the possible application for flight in close proximity but without keeping relative positions. This is most often meant by the word “swarming”. Such a system is rarely found in literature as the flocking behaviour is inefficient in terms of energy and agents, when provided with a stable dynamics, will evolve towards a common velocity with no fluctuation. From an engineering and physical point of view, this minimises the energy dispersions. It emerges as result of a process that, through mutual influence, drives the agents towards the so-called “consensus achievement”, which will be detailed in the next section. CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 50 It is in theory possible to produce continuous fluctuations about the common velocity by just inserting noise in each agent dynamics or getting rid of the velocity alignment term, that is ui in Equation 2.31. In practical terms, this would not find justification in any engineering application where shuffling agent positions requires resources and, as such, has to be justified by a necessity. In birds flocking this would correspond to the need of confusing predators or the attempt of some birds to pass alerting messages by flying against the main direction of the flock. The latter behaviour can be seen as an information diffusion. The difference in the behaviour of few birds and the rest of the flock can be considered to be caused by the information gradient across the flock. A bird that deviates from the main stream does so alerting the others and producing a manoeuvre wave with consequent position shuffle. This was observed in the context of the Chorus Line hypothesis [113] and is the basis on which the results on Chapter 6 will be constructed. 2.7.4 Consensus, a Common Problem Regardless of whether the system is similar to a flock of birds or to a school of fish, consensus problems remain in both cases as convergence to a given response must be guaranteed. Couzin studied extensively the achievement of a common will in large groups of animals, where individuals manage to share information despite uncertainty, see for example [124–127]. Biologists are not interested into consensus the same way engineers are, their interest is often focussed in linking corporal characteristics of the species to their behaviours in a group. On the other hand, they offer a number of hints to engineers who look for the defining characteristics the single agents have to possess to achieve consensus. A trivial example is the number of group members each individual tracks in a biological group to ensure cohesion and continuous information passage: knowing that this number is usually much smaller than the group size made engineers concentrate on the problem of achieving an always connected communication network rather than having a large number of connections per agents. A stable control algorithm drives the agents of a connected swarm to converge towards a common, stable behaviour. Consensus is hence a form of stability: in particular it guarantees that a stable behaviour for a group is achieved whereby stability of the control law for the single elements is provided. CHAPTER 2. CONTROL PROBLEM, THEORY AND TOOLS 51 The consensus problem is strictly related to the connectivity network of interactions. If in a swarm there are groups of agents which are isolated, in the sense that they might have visibility of some neighbours but none of these belongs to the rest of the swarm, these cannot participate to the group dynamics. Hence an isolation in terms of network will translate into physical isolation when the group moves. Indeed a single, globally reachable node in the graph is needed to ensure the achievement of a consensus as proved for example in [39] and [128]. This stability and consensus achievement characteristic was extended to asynchronous swarms under the hypothesis of connectedness by Gazi [129], taking into account also possible delays in the communication channels. The challenge of achieving consensus quickly in a swarm, with obvious benefits for the reactivity of the system, has been attentively investigated. As was discussed in Section 2.6.2, particular networks can enhance the performance in terms of consensus achievement speed. Olfati-Saber and Murray [130], while pointing out the centrality of connectivity to reach consensus, related a particular class of directed graphs to the property of achieving consensus quickly considering also the effect of delays in communications. The peculiar case of small world networks has been considered to this end as well, for example in [103, 104]. A review of the most popular consensus seeking strategies in networked systems was recently produced by Olfati-Saber, Fax and Murray in [76]. In this dissertation the problem of quick consensus achievement will be analysed in Chapter 6 in terms of the design of a distributed controller, where the resources needed for its implementation are provided by the agents on the basis of their position in the network. Chapter 3 The Cohesion of Emerging Patterns and the Emergence of Bottlenecks The possibility of shaping a swarm by means of the network of interactions is explored here for the case of a simple reduction in the number of connections each agent can establish. A particle system approach is adopted where agents move in a three-dimensional space by means of pairwise Morse potential and leak energy because of a viscous damping term. The main results of this chapter have been reported in [131, 132]. 3.1 The Need for a Cohesive Swarm Within the context of swarming systems it is clear how coherent behaviour of the swarm depends upon reciprocal interactions amongst the units. These interactions are numerically restricted in biological systems and must be similarly in any practical implementation of a large scale swarm. As illustrated in the preceding chapters in these cases, the determinability of the behaviour of the swarm becomes highly complex. It is thus important to determine conditions in which the behaviour of the swarm becomes more or less coherent due to restricted interactions. The effect of limited interactions has been addressed in the literature with stability analyses that rely on the swarm staying connected and numerical simulations in support of the behaviours [61, 96]. At the same time some works have ad52 CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 53 dressed the problem of keeping the swarm connected by producing control inputs in the direction of increasing the connectivity of the graph of interactions, see for example [133, 134]. However, it is fundamental to assess what characteristics the swarm needs at agent level in order to establish a communication network that stays connected as the swarm evolves, even without an active control on the connectivity. Associating agents to nodes of a graph, it is possible to track mutual interactions as edges of the graph and thus conclude characteristics of the system. Consider agents that interact according to pairwise potentials. In this chapter, three potentials will be used to extend the validity of the results obtained beyond the particular form of interaction of the agents. The Morse potential was already introduced in Section 2.6.1 and it is repeated for sake of completeness as, UijM orse     |xij | |xij | r = −C exp − a + C exp − r L L a (2.25) where, xij is the relative position vector of an agent i with respect to an agent j. C a , C r represent the strength of the potential while La and Lr govern the range of the potentials. The condition La > Lr for C a = C r provide the potential with a convex shape, which in turn makes its stationary point a stable equilibrium point, as shown in [28, 32]. The hyperbolic potential was introduced in Section 2.6.1 as, Uijhyp = Ch [(|xi − xj | − d)2 − 1]0.5 . (2.29) where, Ch is a scaling constant and d is a desired inter-agent distance. Finally, the third potential considered is a simple quadratic potential in the form Uijq = (|xij | − d)2 (3.1) where, d is again the desired distance between any two agents. The quadratic potential provides an inter-agent attraction force increasing linearly as the distance between two connected agents increases beyond d and decreasing in the same fashion when the distance shrinks below the desired distance. The equations of motions for the i−th agent of the swarm, regardless the potential CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 54 used, are expressed as, dxi = vi dt dvi m = −∇Ui − σvi dt (3.2) (3.3) where, σvi introduces a velocity dependent damping term in the dynamics, with σ being a positive constant. Throughout the analysis presented in this chapter, σ takes the value of 0.7 as this guarantees a convenient relaxation time in the P numerical simulations. Ui = j (aij Uij ), with aij being the entry of the adjacency matrix that takes values “1” or “0” depending respectively on whether the agents are interacting or not. The presence of a viscous damping guarantees the achievement of a crystalline formation as per the physical system illustrated in Section 2.6.1 for an always positive g(H) function. When this happens, as it will be clear in the following, the connection network doesn’t undergo any further change, as it is established on the base of inter-agent distances. 3.2 A Minimum Number of Links for Cohesion The limit to the number of connections per agent necessary to ensure cohesion can be deduced for some cases depending on the rule agents follow to establish connections. Consider a swarm of N agents whose number of interactions is strictly limited. Assume a simple connection rule: suppose that an agent can sense the potential of at most k other agents, and so it does by sensing the k closest ones in terms of physical positions. As seen in Section 2.7.1, this is consistent with both biological studies and the performances expected from an engineered multi-agent system where a limit to the number of communication channels is to be considered. In the remainder of this chapter the connection scheme produced by the existence of this limit is referred to as the k Nearest Neighbours Rule (k − NNR). When representing the connection network into a sensing graph, the result is a directed graph as it can be understood easily through the following considerations. Two generic agents, call them A and B, are each other neighbours in the graph, if both of them have less than k other agents closer in distance to each of them. If this is the case, then an undirected link is established between the two CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 55 agents. Conversely, if agent A has k or more neighbours closer in distance, then agent B will not be included amongst the neighbours of A; hence, although there might still be a directed edge from B to A, there will be no edge from A to B, so the link is directed. Each node of the graph then has k outgoing edges, that is it senses the action of the closest k neighbours, while it can have a number of incoming edges, that is, it can be sensed by up to the total number of agents in the graph, excluding itself, i.e. N − 1. A general consideration is that, for a graph to be connected, the total number of interactions, hence of connections within the swarm, must be at least N − 1; this is, amongst others, the case of a line graph or a star graph. In the particular case presented, when the connections depend upon the relative positions, and the positions of the agents depend in turn upon the forces exchanged through the interactions, the pairwise potential presented tends to cluster the interacting agents. This prevents the spontaneous formation of such edge-saving connected graphs. When each agent can keep a number of connections greater or equal to the total number of agents in the swarm, that is k ≥ N − 1 excluding self connections, the graph is connected, and in particular it is complete. As more agents join the system, N increases and eventually becomes greater than k. At this point the graph of the interactions is no longer complete, yet it can still be connected. It is straightforward to understand that k = hN/2i (where h·i rounds down to the nearest integer) for each agent is a sufficient, though not necessary, condition for the graph to be connected. Hence, the swarm will remain cohesive. Further, in a dispersed swarm there will be at least one subgroup composed of n agents, with n ≤ hN/2i = k. As the number of connections per agent is greater than the number of agents in the subgroup, there must be k − n + 1 edges for each node in this group connecting to some of the other N − n nodes. As connected agents gather together, driven to relaxation by the inter-agent forces, this will produce a cohesive group. This circumstance is illustrated in Figure 3.1 where 8 agents establish 4 connections each: if the agents found themselves split into isolated subgroups, at least in one subgroup there would not be enough agents to complete all the connections. This does not hold anymore when k < hN/2i. For the case k = hN/2i two clustered, but yet connected, groups arise. This is shown in Figure CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 56 Figure 3.1: The illustration shows how for k = hN/2i the swarm must be connected. Each of the agents in the subgroup on the left must complete its 4 connections by joining the group on the right. 3.2 for a planar case where the emergence of a bottlenecked swarm is evident when the number of agents increases, while the number of interactions per agent is held constant. As the connections that an agent does not establish in its own cluster are established always on the base of closeness, these will be with some agents on the closest region of the other cluster. This gives rise to the dumbbell shape where the agents in the central, narrower part bond the two clusters together and, by symmetry, have the same number of connections to both sides of the dumbbell. Meanwhile, they are sensed by all N agents. Despite the fact that each agent still keeps a constant number of connection, the bottleneck becomes problematic when information have to travel across the two halves of the dumbbell. While the maximum path length is al small as 2 (two hops are sufficient to move from any node to any other), all the paths that put in communication any two agents not directly connected pass through the central agents which may experience high traffic volumes leading to possible chocking problems and lags in the information transfer. In Figure 3.3 the final arrangement after relaxation of a 130 agent swarm with 65 connections per agent is shown together with a graphical representation of the corresponding adjacency matrix for a three-dimensional case. The dumbbell shape is reflected in the adjacency matrix of the graph as well if each node is associated to an agent and the nodes are labelled from one-end of the dumbbell to the other. CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS (N = 60) (N = 80) (N = 90) (N = 100) (N = 110) (N = 120) (N = 121) (N = 122) 57 Figure 3.2: Swarm systems in a two-dimensional space relaxing to different shapes as the number of agents N increases while the number of connections allowed per agent is held constant at k = 60. 58 CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 0 20 row index i 40 60 80 100 120 0 20 (a) 40 60 80 column index j 100 120 (b) Figure 3.3: Dumbbell emergence in a three-dimensional case. (a) A 130 agent swarm in a dumbbell shape due to the number of connections per agent being limited to 65 and (b) the corresponding adjacency matrix obtained associating each node to one agent, after having them sorted from one end to the other of the dumbbell. Dots represent nonzero entries. 3.3 The Limit Case and the Cheeger Constant Numerical simulations and logical arguments presented in the previous Section show dramatically how hN/2i is a critical value for the number of connections per node. In the following, the bottleneck characteristics of the system for the critical case of k = hN/2i are studied considering the Cheeger constant of the graph. In graph theory the Cheeger constant or Isoperimetric number is a numerical value associated to a connected graph that identifies the groups of vertex which are more poorly connected to the rest of the graph. A small value of the Cheeger constant translates to having a group of vertices in the graph that are reachable through a smaller number of paths compared to the total number of paths that link the other vertices in the graph. In order to define the Cheeger constant, the definitions provided by Fan Chung in [135] are used. Consider a graph G{V, E}, according to the symbols defined is Section 1.4.3. Let S be any subset of V in the graph and S is its complement, while ∂S is the subset of all edges connecting a node in S to one in S, that is ∂S = {{i, j} ∈ E(G) : i ∈ S and j ∈ S} . (3.4) CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 59 For any subset S the volume of S is defined as, vol(S) = X dx (3.5) x∈S where, dx is the degree of the generic node x, equal to the sum of the weights of all the edges exiting from x in case of G being a weighted graph. Relying on the above definitions, the Cheeger constant referred to a subset S ∈ G is defined in [135] as, hG (S) = |∂S| min{vol(S), vol(S)} . (3.6) For an undirected graph, the Cheeger constant is the minimum value achieved by Equation 3.6 over all possible partitions S of G, that is h(G) = min hG (S) S⊂V . (3.7) When the Isoperimetric number is small, there is at least a partition S of the graph that is hardly reachable, which means there is a bottleneck in the network of paths that connect the nodes of the graph. When the graph G is direct, the sum of edge weights cannot be used to calculate the Cheeger constant that has to be obtained, for an oriented graph G, as [136], h(G) = inf S⊂V F (∂S) min{F (S), F (S)} (3.8) where, F = [fij ] is a function known as “circulation” that takes a nonzero value on each directed edge of the graph. F can be defined on the basis of the probability distribution matrix P = [pij ] and its dominant left eigenvector, φ. This is identified, here and thereafter, as the eigenvector corresponding to the rightmost eigenvalue on the complex plane. Circulation can then be defined as, fijφ = φi pij (3.9) where, i and j are indexes corresponding to generic nodes in the graph while φi is the i − th entry of the dominant left eigenvector. P = pij is a matrix whose generic entry i, j gives the probability of moving from vertex i to vertex j based on the number of links departing from i, accounted through the entries of the CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS adjacency matrix, that is aij pij = P j aij . 60 (3.10) Equation 3.9 does not directly provide information on the behaviour of the Cheeger constant as a function of the number of nodes. However, other circulations can be used to define the Isoperimetric number as long as they satisfy the condition, X i i→j fij = X fjw (3.11) w j→w for the generic nodes i, j and w, that is the flow in one node of the graph is null. Equation 3.11 is to be interpreted as balance between the flow entering in node j from all the neighbouring nodes accounted by the indices i and flow leaving the node towards all the neighbours accounted by the indices j. For the matrix F this translates in having the sum over the rows equal to the sum over the columns, that is the sum of the elements of the first row is equal to the sum of the elements of the first column, the sum of the elements of the second row is equal to the sum of the elements of the second column and so on. Therefore, instead of using a circulation based on the dominant left eigenvector as in [136], a new circulation is derived. This leads to an analytic expression for the Cheeger constant, which is linked directly to the number of agents in the swarm being dependent on the number of nodes in the graph. The new circulation can be defined by constant values associated to the edges of the graph grouped as follows. Two edges will be associated to the same value of the circulation if the nodes to which they point have the same number of incoming connections, that is are observed by the same number of agents. This implies that if a node has s incoming edges, then all s will be associated with the same value of the circulation, that is also the same for all edges entering nodes with a total of s incoming edges. This is illustrated by the sketch in Figure 3.4. This new definition for the circulation requires that the number of nonzero elements in each column of the adjacency matrix is known. This corresponds to considering known the number of incoming edges for each node. For this reason the adjacency matrix is idealised as composed of two blocks plus some linking rows. To be consistent with the connection rule, each row of the idealised adjacency matrix shall still present hN/2i nonzero entries. The cases of even and odd CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 61 Figure 3.4: Circulation definition. The edges of the same colour will be associated with the same value of the circulation as they all enter in node with the same number of entering edges. This circulation is composed of just 2 values, one represented by the colour blue for edges entering nodes with just one incoming edge, and the other by the colour red, for edges entering nodes with a total of 3 incoming edges. numbers of nodes can be detailed as follows. Even Number of Nodes Two agents, by symmetry, will find their position in the centre of the dumbbell, as shown in Figure 3.2 for N = 120. The agents in the centre are sensed by both the clusters and their columns in the adjacency matrix do not have any zero entries, except along the main diagonal. This idealised approximation is shown in Figure 3.5 for a graph composed of 60 nodes compared to the one resulting from a spontaneous relaxation of 60 agents following the dynamics earlier described, with initial conditions randomly distributed in a unit cube. The idealisation of the adjacency matrix in Figure 3.5.b represents the connection structure that produces the smallest possible bottleneck compliant with the physical restriction of having at least 2 central communicating agents for N even. This leads the corresponding Cheeger constant to be the minimum achievable, although, given the spontaneous relaxation of the agents, a spatial configuration that gives rise to this minimum value is not guaranteed to occur as shown in Figure 3.5.a. Being defined only for existing edges, the circulation used has a matrix representation with the same nonzero entries as the adjacency matrix. The actual value for each entry is found by requiring that entries belonging to columns with the same number of nonzero elements take the same value. This, together with the condition expressed by Equation 3.11, that is a zero net flow at the node, provides a number of relations that are sufficient to fully define the circulation for the matrix in Figure 3.5.b. This circulation takes four possible 62 CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS values, as there are four different column lengths in the matrix; call these a, b, c and d. Four equations are hence derived by requiring the sums along the rows and the columns of the circulation matrix to be equal, as per Equation 3.11, one for each column length as            N 2 N 2

− 2 a + b + N4 − 1 c + 2d

− 1 a + N4 − 1 c + 2d

− 1 a + b + N4 − 2 c + 2d
(N − 1) d = b + 2 N4 − 1 c + d
−2 a =
−1 b = N c = 2 N 4 N 4 N 4 (3.12) where, again N is the number of agents/nodes in the graph. To allow a solution other than the zero solution, a is considered known and a ∈ ℜ+ . The solution of the system in Equation 3.12 can be expressed as function of a as,   b =    c =     d = N −2 a N N 4 a−b−2d N −1 4 N−2 − N +N 2 N +2 = 2−N N +N −2 4 N 4 2−N + N 2 N 2+N −1 (3.13) a a 0 0 10 10 20 20 row index i row index i The Cheeger constant can thus be derived from its definition of Equation 3.7 as 30 30 40 40 50 50 60 0 10 20 30 40 column index j (a) 50 60 60 0 10 20 30 40 column index j 50 60 (b) Figure 3.5: Adjacency matrices for spontaneously relaxed (a) and idealised connected (b) 60 agent swarm. CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 63 follows, h(G) = h(N) = N 4 −1
N 2
−2 a+ N d 2 N 2
+ N4 − 1 c

− 1 b + N4 − 1 N2 c + (N − 1) d . (3.14) This, after long but straightforward algebraic manipulation, can be reduced to h(G) = h(N) = 2N (N 2 − 4N + 6) N 4 − 2N 3 − 6N 2 + 20N + 8 . (3.15) Equation 3.15 is obtained in the hypothesis of N being a multiple of 4, otherwise more complicated expressions are to be defined that take into account the shifting of the rows of one position depending on the value of N. In this case Equation 3.14 should be used considering the minimum achievable for each configuration rounding up or down the values of N/4 that are used is the equation. Equation 3.14 can be proved to be a decreasing function of N. It is composed of a numerator which is a first degree polynomial in N and a denominator which is a second degree polynomial in N. It follows that for N positive and integer this is an always decreasing function. In particular, for the same reason, h(N) tends to zero as N approaches infinity, that is, lim h(N) = 0 N →∞ (3.16) The Cheeger constant so defined does not depend on any of the coefficients as they are all proportional to the variable a, which eventually cancels out by substituting Equation 3.13 into Equation 3.14, that makes all the terms proportional to a. Odd Number of Nodes The same procedure can be followed for an odd number of nodes. The minimum number of bridging units, in this case, is one, with the two halves of the dumbbell comprising hN/2i = N 2−1 agents as shown in Figure 3.2 for N = 121. In this case, a circulation can be composed of as few as three values, as three is the number of different column length found in the adjacency matrix describing the minimal link arrangement. When hN/2i is odd (e.g. N = 31, 51 etc...) the positions of CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 64 the nonzero entries of the central row of can be shifted one left or right of one cell. This is understandable considering the fact that when hN/2i is odd the central agent cannot split its connections equally between the two halves. The following derivations are performed considering hN/2i being even as and the central agents splitting evenly its connections between the two halves as in the case shown in Figure 3.6 for a swarm of 61 agents. Under these condition, hN/4i = N 4−1 holds. Following the same derivation described in Section 3.3, in the hypothesis just 0 10 row index i 20 30 40 50 60 0 10 20 30 40 column index j 50 60 Figure 3.6: Adjacency matrices for idealised connected 61 agent swarm. introduced, the values of the circulation associated to the adjacency matrix in Figure 3.6 are obtained as, 
N N N   h 2 i − 1 a = h 4 ia + h 4 ib + c = h N4 ia + h N4 ib + c h N2 ib   (N − 1) c = h N2 ib. (3.17) Considering a known and a ∈ ℜ+ as for the case of N even, the solution of the system in Equation 3.17 is then,   b =          c = ia hN 4 hN i 2 hN i−h N i− N−1 2 4 hN i 2 N −1 hN ia 4 (3.18) hN i 2 i−h N i− N−1 hN 2 4 . CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS Considering hN/2i = N −1 2 and hN/4i = ( b = c = N −1 4 N −1 N −3 N −1 2(N −3) 65 these further simplify as . (3.19) The Cheeger constant can thus be derived from its definition of Equation 3.7 as follows, h N2 ic
h(G) = h(N) = N . (3.20) h 4 i h N2 i − 1 a + h N4 ih N4 ib In the hypothesis hN/2i being even, Equation 3.20 can be algebraically manipulated to obtain, N −1 (3.21) h(G) = h(N) = 2 N − 4N + 5 As before, it can be easily verified that lim h(N) = 0 N →∞ (3.22) being the denominator a polynomial of higher degree of the numerator. Also in this case the Cheeger constant tends to zero as the number of agents increases towards infinity. 3.4 Effectiveness of the Model The results obtained in the previous section are checked against numerical simulations of particle swarms relaxing under the action of the potential functions described in Section 3.1. A comparison is made between the above analytically determined Cheeger constant as a function of only the number of nodes, and the Cheeger constant obtained through the spontaneous relaxation of agents driven by pairwise potential, that is to say creating a spontaneous network of links based on the k − NNR. With reference to Equations 2.25, 2.29 and 3.1, the parameters used to shape the potentials are C a = C r = La = Ch = 1, Lr = 0.2. The desired distance in the hyperbolic and quadratic potentials d is set to 0.1 in order to obtain a separation of the same order of magnitude of the one obtained with the Morse potential. Finally, reference to Equations 3.2 and 3.3, the velocity dependent damping constant sigma is set equal to 0.7 to guarantee a convenient relaxation time in the simulations, while the mass of each agent is set unitary. 66 CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS The simulations are performed providing the agents with random initial positions and velocities in the unit sphere; the simulations are run for 35 simulated seconds and the state and the integration is performed using an explicit Euler scheme with a time step of 10−3 seconds. The values of the Cheeger constants obtained by averaging 100 tests are presented in Figure 3.7 for an even number of agents. A comparison for an odd number of agents is also presented here and shown in Figure 3.8. The tests are compared to analytic detemined Cheeger constant also for odd values of hN/2i by developing the adaquate expression needed through the same procedure described in Section 3.3. 0.45 Average Morse potential Average quadratic potential Average hyperbolic potential Upper bound Morse potential Upper bound quadratic potential Upper bound hyperbolic potential Lower bound Morse potential Lower quadratic Morse potential Lower bound hyperbolic potential Analytic lower bound prediction 0.4 0.35 h(G) 0.3 0.25 0.2 0.15 0.1 0.05 0 20 40 60 80 100 120 140 Nr of Agents Figure 3.7: Comparison of Cheeger constant obtained by the analytic expression and the one obtained by making a swarm of agents relax through pairwise interactions for even number of agents. Average values for each potential and extreme values obtained over all the numerical tests are reported as well. It is demonstrated that the Cheeger constant obtained from the adjacency matrix made up of two blocks and one or two linking rows tracks closely the numerical data, although being always below, or at most equal, to the least value. This result holds independently from the kind of potential used but differences are appreciable in the clustering produced by the difference in the interactions. In 67 CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 0.1 Analytic lower bound prediction Average Morse potential Average quadratic potential Average hyperbolic potential Upper bound Morse potential Upper bound quadratic potential Upper bound hyperbolic potential Lower bound Morse potential Lower bound quadratic potential Lower bound hyperbolic potential 0.09 0.08 0.07 h(G) 0.06 0.05 0.04 0.03 0.02 0.01 0 20 40 60 80 100 120 140 Nr of Agents Figure 3.8: Comparison of Cheeger constant obtained by the analytic expression and the one obtained by making a swarm of agents relax through pairwise interactions for odd number of agents. Average values for each potential and extreme values obtained over all the numerical tests are reported as well. Some of the mamimum values have been cut off the plot to make the curves more clear. CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 68 particular, it can be noticed how the Morse potential guarantees an higher degree of robustness compared to the other two as it provides an arrangement with more links across the two halves of the dumbbell. The comparison between quadratic and hyperbolic potentials reveals how for these two interactions, a more or less robust behaviour is dependent on the number of agents. Overall, the analytic curve, which has been proven to converge to zero as the swarm size increases towards infinite, provides a lower bound for the prediction of the Cheeger constant and can be easily determined for very large swarms, where calculation of isoperimetric number using eigenvalues becomes problematic. Results show that, on average, the graph structure the swarm creates keeps a number of crossing links between the two clusters, in excess of those strictly needed to ensure cohesion. However, while the actual value of the Cheeger constant depends on the final configuration of the swarm and the form of inter-agent potential, its lower bound only depends on the number of agents in the swarm, and, although representing a limit case, it is sometimes achieved by the Cheeger constant. Moreover, as the lower bound of the Cheeger constant is a decreasing function of the number of agents, it can be argued that the links between the two halves of the swarm are an always smaller minority of the total. That is, when swarms grow larger, the number of links that help towards cohesion reduces as a function of the total swarm size if compared to the total number of interactions. Thus as the Cheeger constant is a measure of the bottleneck characteristics of a graph, the results show how a swarm of agents that interact on the base of the k − NNR, with k = hN/2i, tends to become more and more bottlenecked as the number of agents increases. It can be concluded that an increase in the number of agents is not entirely compensated by an increase in the number of cross-links between the two clusters of the dumbbell that the system eventually relaxes into. The Cheeger constant can then be bound from below just by knowing the number of agents and by considering the approximation of the network configuration produced by the k − NNR. Thinking towards the engineering of multi-agent systems, in the case of a very large number of agents, possible consequences arising from the behaviour described are even more incisive. Emergence of a bottleneck restricts sensing and information flow, hence updates of the system’s state, which agents need for cohesion, is delayed, which in-turn directly results in fragmentation into sub-groups even in the case of k close to, but still greater than, hN/2i. CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 69 A systematic study would be required to better understand the effects of the Cheeger constant reduction on the engineering of multi-agent systems. This is not addressed here and represents a possible direction for future developments. 3.5 Final Considerations The butterfly or dumbbell shape emerging here, as a result of decreasing the connections per agent to hN/2i, is not common in biological swarms, but it is of great relevance in engineering systems. The simple reduction in the number of available connections leads to a shaping of the formation in a predictable way. This is different from the distortions in the crystalline patterns observable when a particular network is imposed as shown in Figure 2.7. The shaping towards a more elongated formation starts appearing already when the number of connections is larger than hN/2i, however, it is just at k = hN/2i that the limit case presents itself. This produces the narrowest bottleneck achievable for spontaneously relaxing agents driven by inter-agent forces. With a further reduction of the number of connections per agents, the swarm cannot guarantee cohesion without actively controlling it. The dumbbell shape is a predictable outcome of the connection reduction. Beyond the physical shape, the risk of bottlenecks in communications has to be considered. In an APF approach, as the one just discussed, the shape is the most visible consequence of the connection reduction, but even more alarming is the reduction in the Cheeger constant which indicates how a small number of agents become responsible for the cohesion of the whole group. These can be expected to experience a large volume of information flow to allow communications between the two halves of the group. This in turn sets the design parameter for the communication capabilities of the agents, or alternatively, elects as central agents those ones which can guarantee the best performance in terms of data relaying. When N increases, k = hN/2i becomes a large number of connections to keep for single agents, however this is a condition to guarantee cohesion without actively controlling it. From the design point of view a trade-off should be made to understand whether a multi-agent system would benefit more from a active cohesion control as in [133] or, as in this case, from a cohesiveness guaranteed regardless the control. CHAPTER 3. EMERGING PATTERNS AND BOTTLENECKS 70 While the cohesion is not dependent upon the control law chosen, the shape, under some respect, is so. The different behaviour of the Cheeger constant with different potential witness how a change of agent interactions may reflect into slight changes in the final arrangement. In particular, More potential seems to cluster agent in the proximity of the bottleneck, more than the quadratic and hyperbolic potentials. This is confirmend by the higher values achieved on average by the Cheeger constant, providing a beneficial effect towards keeping system safe from fragmentation. Similar difference can be expected when the tuning parameters of the potential ar tuned within the stability region outlined in Section 3.1, although this was not shown in this chapter. However, despite the fact different potentials return slightly different results, the physical limit of the Cheeger constant is confirmed beyond the particular potential used, which is the main result reported in this chapter. Chapter 4 Fractal Fractionated Architectures and Arrays Application in Space In Chapter 3 it was shown how topographic characteristics of the network reflect the final arrangement of the agent crystalline formation. The possibility of shaping a crystalline formation through changing the network of interactions can be exploited in the design of formation flying architectures. In particular, the design of the final arrangement of a number of agents can be worked out starting from the arrangement of agents in smaller groups, and then, instead of increasing the group size, the number of groups can be increased. The same control laws that drive the arrangement of the agents in a group can be used to drive the arrangement of the groups in larger formations, with the object of the relative positioning control shifted from agents to groups. This technique can be used iteratively to produce always larger formations with a well-structured communication network. The use of the APF method allows for stable behaviour, which is analytically provable. In this chapter the control of a formation of satellites in Earth orbit is designed to deploy a fractionated antenna. The self-similar characteristics of the control law and of the connection network provide a self-similar arrangement of the spacecraft and hence of the antenna elements carried by them. This enhances the performance of the antenna array in terms of directivity. The main results of this chapter have been published or presented in [137–140]. 71 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 4.1 72 Fractionated Spacecraft and Formation Flight A fractionated spacecraft is composed of several independent modules that, flying in close proximity, interact with each other as the subsystems of classic, monolithic spacecraft do. Having the functional units hosted on board several platforms, rather than a single one, provides an increased degree of flexibility and task sharing. A single subsystem can be used in this sense by more than one spacecraft. Although not strictly necessary, fractionated spacecraft concepts often make use of formation flying to rely on precise positioning of the different spacecraft parts beyond the property of being in a close range to allow interactions. Formation flight techniques are gaining popularity for space science, remote sensing and telecommunication applications [141–145]. Formation flying concepts proposed so far have been based on a relatively low number of cooperating spacecraft, as in the case of Lisa, Proba-3 or StarLight missions [146–148]. As the number of spacecraft increases, system complexity increases as well, with significant consequences for the operations and control aspects. Nonetheless, the exploitation of a fractionated architecture and formation flight with an increased number of elements, which maintains an acceptable level of system complexity, can be pursued through the control of autonomous and independent agents as a single group entity [142, 149]. Formation flying techniques in space have been proposed in the area of communications or for applications where electromagnetic radiations have a primary role, see for example the Olfar project [150]. The possibility of producing complex patterns using spacecraft in a reliable way will enable the potential of grouping a number of antenna elements into a cooperative structure. This has long been known and applied in antenna array theory [151, 152] and the move to spacebased architectures appears extremely promising. The key point in the exploitation of formation flying techniques for the deployment of an antenna array is the performance of a homogeneous pattern of array elements can be matched or surpassed by fractal geometries as per [153] and [154]. Fractal geometries as defined by [153] can be considered self-similar structures propagated from a core initiator through a number of stages of growth by CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 73 an identical generator. Application of fractal geometries in antenna array design has mainly focussed on single structures, that is to say one device housing the antenna array. In this chapter a new scenario is considered where, each satellite houses an antenna that contributes to form the fractal pattern. Hence, the problem turns into producing a fractal pattern from a formation of spacecraft, which provide a platform for a number of array elements able to exploit the fractal pattern characteristics. From a control point of view this can be realised through Artificial Potential Functions. The way to obtain complex formations through APF, while maintaining a high-degree of reliability and analytically provable characteristics, can be revealed through the design of a limited connection network. As seen in Section 2.7, network characteristics reflect on the final pattern deployed through APF acting along its edges. In particular when the connection network presents selfsimilarity characteristics, i.e. the same network structure repeats for nodes and groups of nodes, this impacts not only on the final formation but also on the stability and robustness properties, which can be analysed in the same fashion for single spacecraft or groups of those. The overall control architecture result is scalable and with a good degree of fault tolerance. The most inner nodes belong to more than one subgroup of spacecraft and, as such, are more strongly embedded into the network structure and more resistant to connection loss. The most peripheral nodes just contribute marginally to the formation shaping and a fault at this level is not catastrophical for the formation. This will be clearly explained later on in this chapter. From the array point of view self-similarity and sparseness, that is the low number of elements per unit area, lead to a number of benefits: similar performance in operation across a number of frequencies becomes possible due to the repetitive nature of the array pattern; additionally array performance degrades gracefully with element failure and finally equivalent performance can be achieved for a fraction of the number of elements used in square lattice arranged arrays [155]. An ensemble of N spacecraft is considered divided in subgroups of n agents such that N = nk for some integer k > 0. It is assumed that each spacecraft carries an element of the array where the pair spacecraft-array element will be CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 74 named as agent. Spacecraft and array element will be instead used when referring to these parts in any of the agents. The agents are connected according to a non-directed graph described by the adjacency matrix A of dimension N. The spacecraft are controlled through APFs that act only along the edges of the graph with a pairwise scheme. There is no global position or orientation of the agent formation, but within the formation, relative positions are considered for agents and groups of agents while relative orientation is considered for groups of agents only. This implies that the single array elements are assumed to be correctly pointed or, which is what is done here, they are assumed as isotropic sources. In this section it is shown how it is possible to obtain a self-similar formation starting from mutually interacting agents and how the array performances can be analysed for such a system. Artificial potential function characteristics and communication graph topology are described as well. The fundamental concept of applying fractal geometries to the design of antenna arrays using a self-scaling method is described. Only planar configurations are considered but the same argumentations can also apply to linear as well as 3D formations. 4.1.1 Artificial Potential Functions The spacecraft are controlled through artificial potential functions operating along the edges of a communication network. The APFs operate on a pairwise basis, that is they do not depend on position or velocity of the agents but only on their positions relative to the other spacecraft which they are connected with; in particular the Morse potential is used. This was already introduced in Section 2.6.1 and used in Chapter 3, but here, differently from what done previously, the tuning parameters are not the same for all the agents: they depend on the pair of agents they refer to. The parameters C a , C r , La and Lr are redefined and associated to a pair of agents. For the generic pair of agents ij the attractive and repulsive components of the potential are respectively defined as, Uija Uijr   |xij | = exp − a Lij   |xij | r = Cij exp − r Lij −Cija (4.1) (4.2) where, Cija , Cijr , Laij and Lrij are constants whose values shape potentials sensed by agent i because of interaction with agent j; xij is the relative position vector of CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 75 agent i with respect to agent j. The control law is completed by a virtual viscouslike damping in the form σvi . This control law together with the hypothesis of no external disturbances and idealised sensing and actuation capabilities results in the following equations of motion, ẋi = vi (4.3) mv̇i = −∇Uia − ∇Uir − σvi (4.4) where, Uia = X (aij Uija ) Uir = j X (aij Uijr ) (4.5) j and aij is the entry of the adjacency matrix as defined in Section 4.1.2. This approach for controlling spacecraft formations was already undertaken in literature, see for example [38, 47, 50, 55] 4.1.2 Adjacency Matrix As reported in Section 4.1.1 agents communicate through a network of links. In general in a network system studied through graph theory, an adjacency matrix is a matrix which presents a nonzero entry in the ij location whenever there is a directed edge from node i to node j, which corresponds to a communication link between the two agents represented by nodes i and j. Moreover, the matrix is not weighted, that means its elements only take 0 or 1 values; the “weight”, i.e. the strength of the interactions, is provided by the APF used. The adjacency matrix here proposed is symmetric, hence the graph is not directed but this does not imply that the virtual interactions amongst the agents are symmetric. Along the diagonal the edges belonging to fully connected n-agent groups form a block-like matrix. In each group of n2 agents the groups of n agents are directly linked through 2 linking agents (4 directed edges). These account for both position and orientation as shown in Figure 4.1. Extending this scheme to more agents, it can be defined that for any group of nk agents for any integer k > 0 the connections with groups of equal number of agents are ensured through 2nk−1 agents. The adjacency matrices for 25 and 125 agents are given in Figures 4.2 and 4.3. In Figure 4.4 the node degree is reported for the adjacency matrix of dimension 125, CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 76 Figure 4.1: Connections between groups of fully connected agents. Four groups of four agents each keep relative distances and orientations by means of links, in red, connecting agents of different groups. As positions of connected agents can be adjusted and kept, this results into the positions and relative orientation of the groups being uniquely determined. that is the number of links each node is connected to. Nodes are sorted from the central to the peripheral ones. The network is designed such that the peripheral nodes are weaker than the central ones. This means that loss of control of one node due to loss of link is more likely for nodes that belong to peripheral region of the formation, hence they do not act as a bridge between large portions of the ensemble. This implies that the loss of some links is more likely to produce the disconnection of a smaller and peripheral portion of the network than of a large portion. Each node is in any case at least connected to n − 1 other nodes. When the number of generators increases, those groups which were end-points for the previous generator become embedded and more firmly bonded into the larger pattern. This ensures that in the most critical scenario the loss of at least n−1 links is needed for fragmentation to occur. CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 77 0 row index j 5 10 15 20 25 0 5 10 15 column index i 20 25 Figure 4.2: Adjacency matrix for an ensemble of 25 agents. The nonzero entries are represented by dots. 0 20 row index j 40 60 80 100 120 0 20 40 60 80 column index i 100 120 Figure 4.3: Adjacency matrix for an ensemble of 125 agents. The self-similarity of the matrix can be noticed. 78 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 30 25 Degree 20 15 10 5 0 0 20 40 60 80 100 120 Node Index Figure 4.4: Node degrees as number of links belonging to each node. A self-similar scheme can be observed with nodes in central position being the most connected ones. In this scheme the maximum number of connections per node is 28. 4.1.3 Distributed Fractal Antennas The advantage of grouping individual antenna elements into arrays consists in the possibility to produce a highly directional radiation pattern, that is the resulting antenna has a very narrow beam and concentrates the radiated energy in a given direction minimizing its dispersion elsewhere. This can be quantified by the directivity, which is a measure of how efficient an antenna is at radiating its energy in a desired direction. Using a group of smaller distributed antenna elements to steer the main beam, especially in the context of space communications, offers a number of benefits including cost reduction and risk mitigation. Fractal electrodynamics studies the arrangement of these radiating elements into fractal patterns. It is defined as the combined study of fractal geometries with electromagnetic theory and provides methods for the theoretical analysis and synthesis of fractal antenna arrays. When an array presents radiating elements arranged into a fractal geometry, then it is termed as a fractal array. It provides the CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 79 Figure 4.5: Homogeneous square lattice array configuration opportunity to merge the benefits offered by periodic and random arrays, while mitigating some of the side effects. In deterministic fractals the pattern can be achieved by switching elements of a fully populated (square lattice) array on/off until the desired fractal pattern emerges. Following this procedure, the thinned generating sub-array can be copied, scaled and translated to produce the final array. Deterministic fractal arrays created in this manner can be viewed as arrays of arrays due to the recursive nature of the development procedure. The effects of several antenna elements combined in an array are mainly visible in the overall radiation pattern, which changes from the one of the single elements. These changes can be modelled using the so called array factor : it quantifies the effect of combining radiating elements in an array without the element specific radiation pattern being taken into account. Figure 4.5 shows a symmetric planar array with uniformly spaced elements, dx and dy apart in the x− and y−directions, where each element is considered an isotropic source; spacing of elements is commonly measured in fractions of wavelengths. In terms of mathematical modelling the array factor Γ in a given direction, identified by the azimuthal angle φ and the elevation angle θ, enters the directivity in the sense that it can be used to quantify the amount of energy radiated in a certain direction, which is then to be compared with the total energy radiated. Directivity is CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES hence: Γ2P (θ, φ) DP (θ, φ) = 1 4π R2π Rπ 0 0 80 (4.6) Γ2P (ϑ, ϕ) sin(ϑ) dϑ dϕ where, the array factor in the (θ, φ) direction is compared against its integral. The integral sweeps an hemispherical surface because of the planar configuration of the array. The subscript P refers to the level of growth of the fractal as ΓP attains different values depending on the stage of growth of the fractal. In particular, defining δ as the growth or expansion factor, that is the number of elements on the side of a square lattice needed to produce the initiator, the basic structure of the fractal, then the array factor is ΓP (u) = P Y Γfrac (δ p−1u) , (4.7) p=1 where, u = (θ, φ). δ controls how much the array grows with each application of the generating sub-array and corresponds to the number of elements on the side of the symmetric planar array before the thinning. The frequency at which the array operates scales with the fractal configuration starting from the frequency f0 at which the individual elements are designed to operate fP = f0 δ −P . (4.8) Conventional symmetric arrays are typically designed to operate at a chosen frequency f0 in the sense that they are expected to focus most of their energy in the main lobe rather than radiate their energy evenly in all directions. Fractal arrays are known to exhibit multiband characteristics as demonstrated by [153] and [154], which concluded that infinite fractal arrays have the same array factor at any infinite number of bands. For bandlimited realisation similar properties only exist for as many bands as growth stages, i.e. P . Fractal distributed arrays present a number of appealing features for space architecture to produce their deployment. The use of some characteristics of APF allows a Purina fractal arrangement to be produced in a relatively simple way. While fractal arrangement of array elements have been investigated in the context of fractal electrodynamics, the following sections are aimed to illustrate the innovative integration of this fractal shaped antenna in a formation flying architecture controlled through the APF method. A more detailed analysis of the array characteristics emerging out CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 81 of the fractal deployment is contained in the work the author developed with Philippos Karagiannakis et al. [137]. 4.2 The Control Law and the Fractal Antenna In this section the characteristics of the control technique used to drive an ensemble of agents towards the formation of a fractal pattern and the issues related to the design of a fractal shaped antenna array are considered. It is first shown how an asymmetry in attraction-repulsion potential leads necessarily to a central symmetry configuration. It is then shown how the APF coefficients are calculated in order to get the desired distance between agents. Analysis of the control law is completed by considering the nonlinear stability characteristics. The characteristics of the fractal antenna deployed through the satellite formation are finally illustrated for the case of a Purina fractal antenna array [153]. In the remainder of this chapter just the case of an initiator of 5 elements is considered, that is, with reference to Section 4.1, n = 5. 4.2.1 Central Symmetry Emergence Central symmetry emerges at initiator level by means of an asymmetry between the interactions of an agent with the group. This is obtained through a different value of the Lrij parameter along the directed edges connecting the agent to the other 4 in the initiator structure. This is herein explained by finding the conditions that make the artificial potential derivatives null along two orthogonal axes centred on the agent considered. Considering a spatial arrangement of n agents the first derivative of the artificial potential sensed by the generic i − th agent can be calculated as in Equations 4.9 and 4.10. n   Cija |xi − xj | exp − − Laij Laij   n  |xi − xj | ∂Ui X Cija exp − − = ∂yi Laij Laij j=1 ∂Ui X = ∂xi j=1    Cijr |xi − xj | xi − xj exp − (4.9) r r Lij Lij |xi − xj |   Cijr |xi − xj | yi − yj exp − (4.10) r r Lij Lij |xi − xj | Excluding the trivial case for Lrij = Laij and Cijr = Cija , Equations 4.9 and 4.10 can be made null while satisfying the stability conditions in [32] Lrij < Laij . From here on, only changes in the Lrij parameter are considered, where i, j refers to CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 82 their indexing within the 5 agent group, while Laij , Cija and Cijr are considered independent from the pair of agents, that is, these take the same value for every index i, j, hence the indices will be dropped. Consider the planar formation in Figure 4.6, which will be used extensively throughout this chapter. Agent 1, by symmetry, is in equilibrium for the y component of the potential. Figure 4.6: 5 agents arranged in a homogeneous formation due to all-to-all potential with the same coefficients for all the agents The equilibrium for the y-axis holds for all possible distances d either in case Lrij = Lr for all ij, that is, it takes the same values along all the edges, or in the case agent 1 has a different repulsive scale distance. In the reference frame chosen any point along the x-axis has null y derivative of the potential because of the symmetry of the formation about the x-axis. Equilibrium along the x-axis for agent 1 does not lead to an explicit expression for the equilibrium distance, nonetheless the x derivative of the potential referring to any agent can be calculated. Due to the homogeneity of the configuration, any agent can be taken to analyse the artificial potential field. In particular for agent 1 ∂U1 ∂x1 where,       Ca d d2 = 2 a exp − a (cos α) + exp − a (cos β) L L L pentagon       r d d2 C −2 r ′ exp − r ′ (cos α) + exp − r ′ (cos β) L L L (4.11) CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES d d2 = 2 s 1 tan α + cos α 2 + 1 = kd 83 (4.12) can be determined, where k > 1. This is considered as an initial equilibrium scenario for some equilibrium distance d and for Lr = Lr ′ that is the same repulsive scale distance sensed by all the agents. In this scenario, Equation 4.11 must return zero but if Lr 6= Lr ′ and in particular Lr ′ < Lr the separation distance must shrink, that is the equilibrium distance reduces as the scale separation distance shrinks. This can be verified by differentiating Equation 4.11 with respect to Lr ′ . This returns dU1 dx1       Cr d kd = 2 r ′2 exp − r ′ (cos α) + exp − r ′ (cos β) L L L    r    r r C d kd C C −2 r ′ exp − r ′ (cos α) + ′2 exp − r ′ (cos β) . L Lr ′2 L Lr L pentagon dLr ′ (4.13) Equation 4.13 can be shown to be negative, as given in Equation 4.15, that is verifying that a reduction of Lr ′ produces an acceleration on agent 1 in the direction of positive x-axis, hence a reduction of its equilibrium distance:       Cr d kd 2 r ′2 exp − r ′ (cos α) + exp − r ′ (cos β) L L L  r      r r C C d kd C exp − r ′ (cos α) + r ′2 exp − r ′ (cos β) < 0 −2 r ′ L L L Lr ′2 L (4.14)           d d kd kd ∴ 1 − r′ exp − r ′ (cos α) + 1 − r ′ exp − r ′ (cos β) < 0 . L L L L (4.15) ′ This is always verified for d > Lr . The sufficient condition d > Lr ′ can be obtained by a wide choice of system parameters and understood by inspecting the equilibrium distance for the simple case of two agents. This is obtained by summing up and setting equal to zero the derivatives of Equations 4.1 and 4.2 for |xij | = d, and then solving for d C a Lr La Lr ln r a > Lr d= r a L −L C L . (4.16) 84 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES y y 5 5 5 -5 10 x x 5 -5 -5 10 -5 (i) (ii) Figure 4.7: Contours for the potential sensed by agent in (0,0) in the case the other agents have all the same value of the repulsive potential scale length Lr (i) and in the case one of the other agents has a repulsive scale distance Lr′ < Lr In particular, for C a = C r this is verified as long as La 6= Lr , but as stability imposes La > Lr , to make the potential function convex about equilibrium, it can be concluded that this is always verified in this condition and possible to achieve for other choices of C a and C r parameters. The other agents in the group considered tend to keep the same relative distance with respect to agent 1. This produces the new equilibrium configuration that sees the agent with reduced separation distance finding its equilibrium position in the centre of the 5-agent group while fulfilling also equilibrium conditions for the other agents. A contour plot of the potential which agent 1 senses is shown in Figure 4.7 for both equilibrium and non-equilibrium parameter choices. The same could be stated for parameter C r as Equation 4.11 is linear in C r . Here parameter Lr ′ is used to impose the central symmetry configuration over the pentagon one, while parameter C r is used to produce the desired inter-agent distance only. The cross configuration generated by the asymmetry in the potential repulsive scale length is sketched in Figure 4.8 where the reference frame is rotated by 45 degrees compared to Figure 4.6. Considering that interactions amongst agents are only along the edges of the adjacency matrix, a representation of the repulsive and attractive scale parameter, as well as of the other coefficients influencing Equations 4.1 and 4.2, can be given CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 85 Figure 4.8: Cross pattern emerging by shrinking the repulsive potential scale length broadcasted for the agent in the centre. in terms of matrix that has the same structure as the adjacency matrix described in Section 4.1.2. An extract from repulsive distance matrix is reported in the Equation 4.17  0 Lr  r′  L 0   r′  L Lr  r′  L Lr   Lr ′ Lr   r′  L2 0   0 Lr  3   0 0  Lr Lr Lr Lr Lr Lr Lr2 0 0 Lr3 0 0 0 Lr Lr Lr 0 Lr Lr Lr 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Lr Lr Lr ′ 0 Lr 0 0 .. . 0 Lr ′ Lr 0        ...             .. . , (4.17) where, zeros are in the same positions in the adjacency matrix in Figure 4.2 and 4.3, and where the coefficients regulating the interactions among nodes that are centres of two different 5-agent groups are denoted by Lr2 . Finally Lr3 is used to indicate the value along the edges connecting peripheral agents across different 5-agent groups. Hence, coefficients Lr , La , C r and C a can be arranged in square CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 86 matrices of dimension N as they are referred to the edge and take a different value depending on which agents the edge connects. It is important to notice how the pentagon arrangement is not guaranteed by the condition Lr ′ = Lr , that is, the cross configuration is an equilibrium configuration also in the case Lr ′ = Lr . Conversely, having Lr ′ 6= Lr excludes an equilibrium configuration in the shape of a pentagon. It can be concluded that the pentagon arrangement, representing a local minimum configuration, is excluded by the choice of a reduced Lr ′ . As the cross configuration in Figure 4.8 can be obtained for both the choices of Lr ′ considered, the reduction of the Lr parameter produces the exclusion of one of the two configurations, hence, it can be seen as a method for escaping a local minimum configuration. A consequence of this is once the cross configuration is achieved through the reduced Lr ′ , returning this parameter to its original value Lr does not produce any change on the formation. This is because, due to the symmetry, the central agent is in equilibrium whatever choice of Lr parameter is done as long as it guarantees a stable minimum of the potential. The symmetry of the arrangement translates to two pairs of equal and opposite terms for the sums in Equations 4.9 and 4.10 making both equations trivially null. Equilibrium conditions for the surrounding agents according to the scheme of Figure 4.8 is only determined by Equation 4.9 as the y-component is null by symmetry. Equilibrium distance d as reported in Figure 4.8 is found by solving for d the Equation Cr Lr exp  √ !!  √ 2d −2d + exp + 2 exp = r L Lr     √ Ca −d −2d 2 exp exp + exp + La La La −d Lr   √ !! − 2d La . (4.18) Equation 4.18 is obtained by expanding Equation 4.9 for the present case; however, this is not analytically solvable. On the other hand, equilibrium can be found for a given d, by tuning C a and C r parameters. This is better explained in Section 4.2.2. CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 4.2.2 87 APF Coefficient Definition The coefficients of the APF acting along the edges of the graph are calculated such as to set the desired distance amongst the spacecraft. Just the C r coefficient is calculated as function of the others which are set arbitrarily as long as suitable to produce an equilibrium configuration according to the criteria discussed in Section 4.2.1. The change of C r parameter only or, more precisely, the change in the ratio C r /C a is sufficient to modify the position of the minimum of the potential, hence the design distance, for the APF used. In particular, an interaction between two spacecraft belonging to two different n-agent groups is considered, with a design distance dd ; C r coefficient can hence be calculated by manipulating Equation 4.16 as   Cr Lr La − Lr = a exp dd a r Ca L L L . (4.19) Equation 4.19 can be reversed to calculate the equilibrium distance once the coefficients are set. When more than 2 agents are involved, an analytic expression for the equilibrium distance cannot be defined, but given a desired distance, it is always possible to get an expression for the value of the ratio C r /C a that produces that separation. In particular for a fully connected group of 5 agents C r /C a ratio can be calculated equating to zero the gradient of the potential for the formation according to the scheme in Figure 4.8. As the y-component is trivially null, C r /C a can be calculated considering just the x-component of the gradient in Equation 4.9. It is then possible to solve for the C r /C a ratio and obtain  √ 

√ 2dd dd 2 exp − L2da d + exp − + exp − La La L C  √  = a Ca L exp − dd
+ exp − 2dd
+ √2 exp − 2dd Lr Lr Lr r r . (4.20) This tuning method can be extended to the other links of the adjacency matrix; by defining the coefficients in this way the desired self-similar pattern is produced. 4.2.3 Stability of Control Law Asymptotic stability can be proved by adopting the procedure introduced in [32]. Consider the time derivative of the energy as sum of artificial potential and real CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES kinetic energy: dEt dKt dUt = + dt dt dt where Ut = , 88 (4.21) 1 XX aij Uij 2 i j (4.22) is the total potential energy per unit mass with Uij = Uija + Uijr and Kt = , (4.23) 1X 1X Ki = (vi · vi ) 2 i 2 i (4.24) is the total kinetic energy per unit mass. Expanding 4.21,   ∂Kt dEt X = ∇Ut · vi + · v̇i dt ∂v i i (4.25) where, the operator ∇(·) is the one defined in Equation 2.20. Substituting Equation 4.4 and Equation 4.23 into Equation 4.25   dEt X dKt = ∇Ut · vi + · (−∇Ui − σvi ) dt dv i i ∴  dEt X  = (∇Ut · vi − ∇Ui · vi ) − σ|vi |2 dt i (4.26) . (4.27) As the potential depends upon pairwise interactions, xi derivative is not null for both the Uij and Uji potentials that form the total potential Ut . If the agents interacted in a symmetric way, this would cancel out with the gradient ∇Ui , but as the sum of the potential derivatives upon any agent includes asymmetric terms, this does not occur. Nevertheless the difference between the gradients can always be damped by the artificial viscous damping. Hence, it can be concluded that ∃σ>0: X i  (∇Ut · vi − ∇Ui · vi ) − σ|vi |2 ≤ 0 . (4.28) This is enabled by the property that the Morse potential and its gradients are bounded functions. CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 89 As total energy time derivative can be made a negative semi-defined function, this can be compared to a Lyapunov-like function whose derivative is always proved to be negative and zero at equilibrium, corresponding to null speed. Thus, the system will leak energy and stabilise eventually into a static formation which corresponds to a local minimum of the total energy. 4.2.4 Fractal Antenna Array Design and Analysis The fractal geometry of the formation provides the possibility of reducing the number of elements of a lattice array without significantly diminishing the performance of the array. In a formation flying architecture this is a noteworthy factor as, especially in close proximity, the challenges in position-keeping increase with the number of agents, as the presence of each agent in the group introduces position constraints for the others. Basing antenna array formations on fractal geometries provides not only the potential to reduce the number of elements but also offers the possibility to operate across a range of frequencies and the selfreplicating nature of fractal patterns extends to their performance characteristics. Following the analysis methodology described in Section 4.1.3, a thinned planar array based on the Vicsek or Purina fractal can be considered; this fractal has the following simple sub-array at growth scale one (P = 1):  1 0 1  .   S1 =  0 1 0  1 0 1 (4.29) The array fractal pattern SP at an arbitrary growth scale P ∈ N, P ≥ 2 is given by SP = S1 ⊗ SP −1 , (4.30) with ⊗ denoting the Kronecker product. Hence, the growth factor δ is 3 and the array will increase three times its linear size for each stage of growth. This first stage of growth starts off with a 9-element generating sub-array of uniformly spaced elements; the center and corner four elements are switched on (“1”), while the remaining elements are switched off (“0”) or can be considered completely removed. The Purina fractal array is generated recursively by repeat- 90 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES edly applying this sub-array over P -scales of growth. This is done by replacing each “1” or switched on element by a generating sub-array and replacing each “0” by an all zero array the size of the generating sub-array. The result of the process is illustrated in Figure 4.9. (a) (b) (c) Figure 4.9: First three stages of growth of the Purina fractal array for (a) P = 1, (b) P = 2, and (c) P = 3. Credits P. Karagiannakis (from [137]). The theory developed by Werner [153] about fractal arrays describes in details the above steps illustrating how the combination of the antenna elements has the potential to alter the radiation characteristics of an ensemble of antennas and results in a steerable and highly directive beam. Providing an in-depth analysis of the radiation pattern and the electrodynamics properties of this fractal antenna array is beyond the targets of this thesis. However it is important to understand what the advantages of such a configuration are. This can be done by plotting Equation 4.6 for the arrangements shown in Figure 4.9. This is done in Figure 4.10 where it is evident how as the stages of growth increase, the beam is more concentrated, that is, a greater percentage of the radiated energy is concentrated into the prescribed direction, resulting in a better signal strength for arrays with the same total area. 4.3 Simulation Results The control method illustrated in Section 4.2 is used to simulate a possible operative scenario in which a spacecraft formation is used to form a distributed array in Earth orbit. A geostationary orbit is chosen to simulate the dynamics although the application is not specifically aimed at terrestrial telecommunications. Deployment of a fractal antenna array is simulated where the system is composed 91 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 5 0 −10 −15 P |D (θ,φ) | / [dB] −5 10log 10 −20 −25 −30 P=1 P=2 P=3 −35 −40 −80 −60 −40 −20 0 20 40 60 80 o θ/[ ] Figure 4.10: Directivity plots for the first three stages of Purina fractal arrays shown in Figure 4.9, with an assumed direction of the main beam towards broadside (θ = 0) of 125 radiating elements. The requirements of the system drive selection of actuators and dictate to a certain degree choices regarding agent selection and separation. The method of control and the possibilities offered by reducing the size of individual radiating elements, while maintaining an overall large aperture, drive towards the selection of a satellite in the size range of pico- or nano-satellite suitable for a separation in the order of 1m. This is the separation chosen as the inter-spacecraft distance is still small enough to control motion through mutually exchanged electromagnetic forces and far apart enough to allow for relatively coarse accuracy, in particular at the release from a carrier spacecraft or launcher. The 125 unitary mass agents reproduce the shape of a Purina fractal at a growth stage of P = 3; they are deployed in 25 groups of 5-agent subgroups which is the elementary unit of the formation (N = 125, n = 5). The dynamics of the spacecraft formation is based on Clohessy-Wilthshire (CW) [156] linearised equations in an orbiting reference frame. This is so oriented: 92 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES • the x-axis is tangent to the orbit and parallel to orbital velocity vector; • the y-axis is parallel to angular momentum vector; • the z-axis is orthogonal to the first two and pointing towards the Earth centre. CW equations in the reference frame chosen are ẍ = −2ν ż ÿ = −ν 2 y z̈ = −2ν ẋ − 3ν 2 z (4.31) where, ν is the orbital frequency. Initial conditions were set such as each spacecraft had an initial position randomly picked within a sphere centred on its final position and radius equal to 1.5 times the distance to its nearest neighbour to account for possible initially swapped positions between near agents; initial relative velocities are null. This corresponds to assuming that a carrier spacecraft or launcher releases the agents with coarse accuracy i.e. not completely randomly. Attitude for the single spacecraft is not considered while overall formation attitude control for rotation around x and y axis is guaranteed by positioning control through a parabolic potential that flatten the formation on the x − y plane. Sensors are idealised, that is, the exact position of each agent is known without delay by all the agents it is linked with. In terms of the actuators, although these are not modeled, some characteristics related to the possible use of electromagnetic forces are considered. In particular, actuators of the kind proposed in [157] and [158] are considered. In, particular, in [157] the effectiveness of Coulomb actuators is assessed. The magnitude of the Coulomb forces attainable by two spacecraft in a geostationary orbit, lies between 10−3 and 10−5 newtons for separation distances between 10 and 50 meters. By inspection of Equations 4.31, considering that the orbital frequency ν has an order of magnitude of 10−5 for the geostationary orbit chosen, a unitary ż (corresponding to the drift of one agent out of the plane of the formation at the speed of 1 ms−1 , that is extremely high velocity) would fall within the actuation capabilities of the system. The actuators would then be able to correct the error in this case that represents the most critical scenario to face in a realistic setting. In fact, the z acceleration is the z direction is of the same magnitude 93 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES of the x component while the y component is 5 orders of magnitude smaller while . Electromagnetic actuators, in particular the ones based on Coulomb forces, cannot be used all at the same time to avoid interferences, a duty cycle is set up and the ensemble is split into a number of groups so that any two groups which are active at the same time are relatively far apart. This allows interferences to be neglected. Each group is controlled across a time period of the duty cycle. Over the whole duty cycle each group of agents is controlled for the same amount of time. As a consequence, agents belonging to more than one group (e.g. linking agents between groups) are controlled for longer. The frequency of the duty cycle needs to be high enough to prevent spacecraft drifting away between control periods. This can be bounded from below by considering a linearised version of the control law and computing the frequency of the associated harmonic oscillator. Considering the APF only, the x-component of the acceleration provided to the generic agent i, that is ẍi by the control can be linearised about equilibrium as mẍ˜i = X  Cija j    Cijr −dij −dij exp − r exp − Laij Laij Lij Lrij ) "  #   Cijr Cija −dij −dij − r 2 exp (xi − dij ) exp Laij Lrij Laij 2 Lij  , (4.32) where, it was assumed that the equilibrium position is at a distance d from the neighbouring agents and that these agents are fixed in their positions. The ycomponent of the acceleration is not considered for symmetry reasons. The sum is extended to all the neighbouring agents acting along one axis. As an example, considering the central agent of Figure 4.8, this means that only 2 agents contribute to its oscillatory motion along the orthogonal axes. Equation 4.32 is in the form of a linearised harmonic oscillator perturbed by a constant acceleration. The frequency associated with this system is: v     uX a Cij Cijr u −dij −dij ωi = t − r 2 exp exp Laij Lrij Laij 2 Lij j . The frequency chosen for the control to operate shall not be smaller than (4.33) 94 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES sup ωi (4.34) i obtained by considering all the sets of values defining the control of the groups. For the case under consideration the whole duty cycle lasts 2 seconds, which corresponds to a control frequency of 0.5 Hertz. The 125 spacecraft are considered as belonging to 9 groups, these are reported in the following list while Figure 4.11 identifies them by highlighting their position in the formation. • the 5 5-agent groups at the centre of 25-agent groups (Figure 4.11.a); • the 5 5-agent groups at the top of 25-agent groups (Figure 4.11.b); • the 5 5-agent groups at the bottom of 25-agent groups (Figure 4.11.c); • the 5 5-agent groups at the left of 25-agent groups (Figure 4.11.d); • the 5 5-agent groups at the right of 25-agent groups (Figure 4.11.e); • the agents linking the centres of the 5-agent groups in 25 agent groups (Figure 4.11.f); • the agents bonding the 5-agent side by side in the 25-agent groups (Figure 4.11.g); • the agents bonding the centres of the 25-agent groups (Figure 4.11.h); • the agents bonding the sides of the 25-agent groups (Figure 4.11.i). The connections between each group (consisting of 25 agents) are ensured by pairs of agents instead of groups of agents. This allows a reduction of the computational efforts for each agent and a reduction of the computational resources needed for the simulation. On the other hand, this reduces the control power and slows down the deployment of the formation. Table 4.1 reports the values of the coefficients used. The agent at the centre of the formation (say agent 1) is the only one linked to the centre of reference frame by a quadratic potential in the form Uc = ζ|x1 |2 with ζ as a weighting parameter set to 0.1. The control force is provided through the negative gradient of Uc as before plus a dissipative term in the form −σ ẋi CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 4.11: Groups of agents controlled at each step throughout the duty cycle. 95 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 96 Table 4.1: Numerical values of coefficients used in numerical simulations. Fully connected groups (f.c.g.) Centres of f.c.g. Peripheral between adjacent f.c.g. Centres of 25-agent groups Peripheral of 25-agent groups σ = 0.1 for all the agents Ca 4 Cr La 3.94722 2 Lr 1 Lr’=0.5 1 0.99596 4.5 4 Lr’=2 0.8925 1 2 0.5 2500 2505.3 10 9.9 Lr’=4.5 69.96 70.0 3 2.9 with σ also set to 0.1. This is to provide a kind of orbit tracking capability or, in practical terms, the possibility to stay anchored to centre of the reference frame. The value of 0.1 is chosen in order to make the orbit tracking effective without superimposing disturbing dynamics over the inter-agent APFs. Also, this suggests that the task of tracking the orbit can potentially be carried out by only a single agent, while the others just track their relative position with respect to the central agent. This is not necessarily the one in the centre as it is done here for simplicity. The control law is applied for just x and y axis of the orbital reference frame with control on z-axis performed through the gradient of a simple parabolic potential plus a dissipative term in the form z̈i = −∇Uzi −σ żi , for i = 1...N, with Uzi = ζ|zi|2 where, ζ and σ have the same roles and values as used previously, that flatten the formation on the plane z = 0. Four snapshots from the deployment are reported in Figure 4.12. It can be noticed that after one day the deployment presents some distortions in particular from peripheral groups. Finally in Figure 4.13 errors on the designed relative positions after one day are plotted. More precisely the error is considered as the difference between the actual distance of each spacecraft from the centre of the formation and the ideal design distance; this is then plotted as a percentage of the desired spacing. It can be seen that the greatest error is lower than 5%. Both snapshots and error evaluations are considered after a maximum of 24 hours; this is sufficient to show the self arranging capabilities of the control technique. After a further 24 hours the magnitude of the maximum error is halved. Theoretically a complete 97 10 10 5 5 Y [m] Y [m] CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 0 0 −5 −5 −10 −10 −15 −10 −5 0 5 10 15 −15 −10 −5 X [m] 10 10 5 5 0 10 15 5 10 15 0 −5 −5 −10 −10 −10 5 (ii) Y [m] Y [m] (i) −15 0 X [m] −5 0 5 10 15 −15 −10 X [m] (iii) −5 0 X [m] (iv) Figure 4.12: Formation deployment in GEO. Snapshots taken at (i) t=0 s, (ii) t=60 s, (iii) t=1 hr, (iv) t=24 hr relaxation with no positioning errors is possible but only after an infinite period of time due to the viscous-like damping. On the other hand, the error after finite time can be improved by either optimising the damping constant or replacing the viscous damping with more sophisticated means of virtual energy dissipation. However, both these improvements are beyond the scope of the study presented here althugh represent possible development directions. 4.4 Final Remarks The idea of meeting needs for highly directional wideband antenna arrays through a space based fractionated architecture is constructed around the possibility of locating a number of spacecraft, each carrying an antenna element, according to 98 CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 0.6 1 0.5 0.8 ∆ Y [m] ∆ X [m] 0.4 0.6 0.3 0.4 0.2 0.2 0 0 0.1 20 40 60 AGENT 80 100 120 0 0 20 40 (a) 60 AGENT 80 100 120 (b) 5 DISTANCE IN % 4 3 2 1 0 0 20 40 60 AGENT 80 100 120 (c) Figure 4.13: Errors in design positioning after 1 day from release of the formation. Distances are computed with respect to the agent at the centre of the formation (a) along the x axis, (b) along the y axis and (c) in percentage of the design distance from it. CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 99 a precise fractal scheme. This improves overall antenna performance and capabilities while using a reduced number of elements. In turn the possibility of using small spacecraft enables the formation of a wideband fractionated antenna, but requires accurate spacing between the elements. Orientation is not considered here for single agents as they are assumed to be isotropic sources. Thus, in the case of an antenna array as described above, the relative agent positions within the whole array is the key requirement as this influences the performance of the array. Hence considering just coarse attitude control for single agents, a description of the system characteristics in a global sense is possible as long as relative positions are precisely known. Utilising this knowledge, directivity through array phasing is achievable: at group-level for compensation of global attitude errors and at agent-level to accommodate misalignment of the single elements. From a control point of view, the need for precise close formation flying can be tackled through using reliable techniques and implementing these on relatively small agents. In this respect APFs are particularly suited to the task as their stability characteristics are analytically provable, hence they do not need extensive montecarlo test campaigns to validate their behaviour. Moreover, APFs allow for highly nonlinear control through quite straightforward computation due to their smoothness. As the amount of information needed is just the relative position of a number of neighbours, the connection network presented here has the double advantage of shaping the formation on one side and reducing the number of connections on the other. These combined characteristics make small spacecraft, even with reduced computation capabilities, able to carry out the task of arranging into a formation through exclusively inter-agent interaction in a decentralised way. APFs provide collision avoidance capabilities to the agents directly connected into the network. However, some agents might be found to fly, at least in some phases, close to others which they are not connected with. This means that these agents would not sense each other as they are unaware of their relative positions. The problem can be tackled considering avoidance manoeuvres to operate just in case of close proximity, revealed through any sensor scanning of the local neighbours. Any avoidance action should then be properly designed so as not to introduce persistent instability in the control of the agents already linked through the net- CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 100 work. Moreover, as the communication architecture presented requires each agent to keep a relatively small number of communication channels open, temporary links can be established between spacecraft in close proximity to negotiate on a case-by-case basis a suitable avoidance manoeuvre. The concept of self-similar arrangement is quite general and can be exported from formation flying for the exploitation of electromagnetic phenomena, to other applications. However, in this chapter the space flight case was illustrated, hence the dynamics of the formation in the orbit environment imposes to consider specific orbit parameters and suitable actuators, at least in the definition of the simulation scenario. A geostationary orbit was considered although agents are not specifically targeted at telecommunication purposes. The choice not to detail the actuator modelling in depth was made to mainly concentrate on the control aspects and the benefits the distributed architecture can get out of the hierarchical communication scheme. Actuator characteristics were considered only in part. Although the response of the actuators was not included, their choice took into account the close proximity scenario. The use of inter-agent electromagnetic forces was proposed rather than thrusters as it avoids possible plume impingement problems considered the close proximity at which the agents operate. Moreover the APF methods drive the system through an oscillatory stage before the achievement of the equilibrium configuration during which residual energy (both virtual potential and real kinetic) is dissipated. This translates into fuel wasting when considering the use of thrusters. The introduction of a duty cycle in the control operation is a consequence of the choice of actuators. Also the advantage of having actuators that mimic the virtual inter-agent action of the artificial potential makes the performances of control techniques more evident as less dependent on the actuators. The duty cycle just applies to inter-agent forces, which are supposed of Coulomb type [157], while no hypothesis is done on forces that act along z-axis or that bond the centre of the formation to the origin of the reference frame. For what concerns z-axis, the use of Lorentz forces as in [158] might be considered, although their effectiveness is to be investigated further in relation to the magnetic environment. CHAPTER 4. FRACTAL FRACTIONATED ARCHITECTURES 101 The communication network was intended in the first place for control purposes only, but the need for task assignment in the fractionated architecture as well as array phasing can be carried out through the same architecture. In particular, the system inherits a structured hierarchical network, where the ranking of the agents depends on the number of links they are connected to. This does not imply that the resulting architecture is centralised, but allows the task assignment to be carried out on the basis of the hierarchy of the agents. For instance the guidance for the whole formation can be carried out by a number of spacecraft within the group which communicate in an all-to-all scheme in order to share the computational efforts (e.g. the centres of the 25 agent groups shown in Figure 4.11.h), and then passed to another module able to compare this to the navigation to eventually generate a control input for the whole formation. This is different from the guidance, navigation and control (GNC) functions that each spacecraft carries out: while each spacecraft should find its position in a distributed architecture, the whole system follows a guidance law that enables the mission task achievement. It is worthwhile stressing how the position of each agent is not pre-determined in a strict sense. The links of each agent are pre-assigned, but this does not prevent agents, or groups of agents belonging to the same level, to swap their positions. Finally, despite the planarity of the formation, nothing prevents to produce the emergence of a central symmetry and the building up of several hierarchical levels in a self-similar fashion also in 3D formations. This can be done considering an initiator composed of a different number of agents as well. An analysis on the ground of fractal electrodynamics can be carried out in three dimensions using the same methodology described in this chapter. Chapter 5 Pattern Creation - Experimental Results In Chapter 3 and Chapter 4 the formation shaping and control for swarms of agents was analysed from a theoretical point of view, demonstrating the soundness and the benefits of techniques based on the APF methods, with advantages from modification and design of the communication network. Practical aspects such as latencies in the process, accuracy of the navigation sensors, actuator saturation and dynamics, or possible constraints due to their nature, were not taken into account so far. In Section 4.4 it was pointed out how the choice of actuators may heavily influence the dynamics of the system, but for the purpose of the control technique analysis carried out in Chapter 4, a deeper analysis of the “real world ” issues was avoided. In this chapter, the implementation on real hardware of the control technique proposed in Chapter 4 provides insight into some real world issues that are either difficult to model through simulations, or requiring a number of assumptions that diminish the validity of the modelling. The implementation of innovative technologies on hardware is a fundamental step when climbing the Technology Readiness Level (TRL) ladder, from mathematical modelling up to the point the technology is made available for the final users. The main results of this chapter have been reported in [159, 160]. 102 CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 5.1 103 Real World Issues When departing from theoretical models and looking into specific applications, APF methods are widely used to define the guidance law, that is the velocity field, rather than the forces to produce through the actuators to drive the system to its desired state; this is done, for example, in [161]. APFs are indeed suitable for defining the desired, time varying state of the system, and less appropriate for driving it there. The oscillatory motion produced about the equilibrium is characteristic of this control method and, although the virtual damping presented in Chapter 4 can be improved by a certain amount, using the APFs to define the accelerations of the agents returns improvable relaxation times and energy exploitation. Another fundamental reason that prevents the application of the APF to the definition of the required action from the actuators is the lack of full actuation in a real system. The APF method was applied in the previous chapter to produce a force dependent on the positions of a number of agents. This cannot be generally done in real systems: primarily because actuators cannot produce an action in any direction at any time. Constraints apply to actuator operations, which are due to their nature and to the physical characteristics of the system and the environment. Consider specifically the ground mobility case: wheels are the most widely used actuators and they are extremely limiting if compared to the assumption of free translation along any direction which was made in Chapter 4. Excluding particular and complex designs, wheeled agents are prevented from producing pure lateral shifting. Curvature radius also has to be taken into account for any manoeuvre. As different from the idealised free space case, APF cannot be simply used to produce the motion of the vehicles as a direct response to the potential gradient, which may be in a direction perpendicular to the heading of the vehicle, precluded in the motion of the wheeled system. APFs are very effective in designing the guidance law, that is where the vehicle should move and with which velocity. However, a controller different from the APF that produced the guidance law has to be implemented. This maps the error generated by the difference between the actual state and the desired one produced by the APFs into a command for the actuators at the agent-level. Although this second controller could exploit APFs as well, the advantages of the APF method would be greatly diminished when used to control a single agent rather than a swarm in distributed fashion. For this reason, in the experiments presented in CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 104 Section 5.2, the low level control, that is the control of the actuators, is performed through a quasi-linear function limiting the use of the artificial potential to the definition of the target position to achieve. 5.1.1 Navigation Navigation is the determination of position and direction of motion of an object. In a swarm, navigation must also provide information about the relative position of an agent with respect to the others so to provide a meaningful input to the guidance and control systems which drive the agent motion relying on relative position and velocity. In numerical simulations of swarms, agents are virtual objects for whose motion all information can be known. In the real world swarming agents must be able to estimate their own position and velocity and the ones of their neighbours. This is a demanding task for the agents, and, as it happens also for the control, is the object of specific studies oriented either to the definition of the state of a single agent or to its state relatively to the other agents of the swarm. A detailed analysis of the navigation problems is beyond the scope of this dissertation, hence, the problem of navigation is not faced directly. The results reported in this chapter are obtained using wheeled robots and the navigation problems, which cannot be avoided as in numerical simulations, are overcome by tracking the agents through external sensors and then controlling them in distributed fashion by exploiting this information. 5.1.2 Actuators The problems connected with actuators were partially anticipated at the beginning of this section. They can be summarised into 3 main concepts: - Actuators have their own dynamics that reflects onto the swarm behaviour mainly through delaying and limiting the motion within the swarm; - Actuators are subject to problems such as saturation, deadbands and locking: control system should be robust enough to compensate for these; CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 105 - Actuators interface with the external environment and often exploit its feature to produce forces and torques on the system. Actuators dynamics can be seen in terms of rising time from the moment the input is fed to the achievement of the level of actuation requested, provided that this is feasible. When the input is decided on the base of a time varying guidance law, reaction time of the actuators should be small enough to keep up with it. This is indeed the case of swarming systems where at each instant, desired states for each agent depend upon the continuously changing states of the others. Failing to keep up may result in the loss of convergence to consensus, that is agents do not manage to coordinate their action and produce a coherent behaviour. Saturation, deadband and locking are undesirable states which are present in all actuators. Avoiding operative conditions close to, or with the possibility of evolving into, those states relies on the controller characteristics. Trivially, the control system should not command the actuators to output within their deadband (which would result in null output) or close to the maximum output achievable (which would leave the system with no control margin for contingency manoeuvres, beside stressing consistently the actuators). Finally, most actuators use the external environment to input forces and torques onto the system. For example, wheels use a solid surface and the grip with it to roll while propellers use a fluid flow to accelerate vehicles. The interaction with the external environment is difficult to model as very little is usually known about the environment where the actuators are operating. The design of the control relies on a partial knowledge of its characteristics: the torque to the wheels in a road vehicle should be, for instance, commanded considering the static friction provided by the ground and the characteristics of the tire. There are then characteristics that are peculiar of particular actuators. For a wheeled vehicle, control system shall not command, for instance, out-of-plane motions or lateral shifting. By the same argument, the autonomous controller of a fixed wing aircraft should not command a climbing angle beyond the stalling point of the vehicle. It has instead to work out the best way to minimize the error between the actual state and the desired one by relying on the physical characteristics of the system and its actuators. All these aspects need to be con- CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 106 sidered carefully when passing from a numerical simulation to the implementation of swarm control onto real robotic agents. 5.2 Swarm Robotic Experiments In this section tests performed on a number of robots are described. First the adaptation of the APF control method to the task is presented with a description of the experimental apparatus. Tests of formation achievement, the switching of the formation from a pentagon to a cross shape, and the case of a cross formation rotating about its centre are then described. 5.2.1 Testbed and Methodology Hardware setup The robots used were designed at the Centre for Ultrasonic Engineering and manufactured within the facilities of the University of Strathclyde. These are shown in Figure 5.1. They are composed of two differentially driven motors mounted on an aluminium chassis, a ball bearing is used as third support. They are approximately 175 × 124 × 80 mm with a wheel diameter of 54 mm. The robots are controlled by an embedded Linux computer (720MHz ARM CortexA8 with 512MB of RAM) and powered by a Lithium Polymer battery (11.1V, 2Ahr, providing 4 hour run time). Guidance and control functions were performed centrally on a host computer using a C# Graphical User Interface (GUI) that was interfaced with the robot’s Software development kit (SDK). The SDK interfaces with the robots using WiFi. The robots were tracked by a 6 camera VICON T160 [162] positioning system, which provides 1 millimetre, 6 degree of freedom accuracy at 100 Hz over a 3.8×4 m2 area. Individual robots are identified through a unique pattern of reflective targets (14 mm spheres) affixed to the chassis, these targets are visible in Figure 5.2. In these experiments the robots were driven on the floor of the cell, so only X,Y positions and yaw angle were required to fully characterise the robot position. A more detailed description of the robotic agents used can be found in [163, 164]. Control actions were defined on the basis of the VICON positional estimate, CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS Figure 5.1: Wheeled robots used for swarming test. Not to scale 107 CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 108 Figure 5.2: One of the robots used in the tests. which was passed directly to the robots SDK. The computer emulated distributed control: inputs were produced for each agent in turn on the basis of its position relative to the others. This ensured that each robot was controlled individually hence any group behaviours emerge out of singles’ actions. Virtual robots can be added to the test. These are simulated agents that, through the GUI, can be made part of the experiment increasing the number of agents considered. This feature was not used for the tests presented in the following. A scheme of how the testbed works is presented in Figure 5.3 while the VICON [162] positioning cell is sketched in Figure 5.4. Software The guidance and control laws were coded using C# and, as previously explained, run on a dual core 2.5 GHz, 2 GB RAM, Windows XP desktop computer rather than on the single robots. The results are fed to the guidance law determined by the Artificial Potential Functions that in turn was fed to the control of the wheel speed as it will be explained in the next section. This is done in turn for each object robot considered by the GUI, that automatically labels the robots consecutively from 1 to 5. Robot state, as measured by VICON, was exported to file for post-processing and visualisation in MATLAB R . CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 109 Figure 5.3: Testbed architecture. The Vicon cameras track the robots in the test area and send their positions and attitude to the centralised computer that can also include virtual agents. The control is computed and input to the robots’ wheels are passed through the wireless network. 5.2.2 Guidance and Control The argumentations illustrated in Section 5.1 drive to the conclusion that the guidance must be separated from the control when including experiments with real hardware rather than just simulations. The Guidance Functions Morse artificial potential functions, already used in this thesis, produce the desired velocity for the system ẋdi = (ẋdi , ẏid ). As before, these are composed of exponential functions that provide an attractive and repulsive component defined as, |xij | = exp − a Lij   |xij | r r Uij = Cij exp − r Lij Uija −Cija   (4.1) (4.2) with constants Cija , Cijr , Laij and Lrij shaping the formation and its size. Differently CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 110 Figure 5.4: The testbed arena where tests are performed, with Vicon T160 cameras in red. CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 111 from what done in Chapters 3 and 4, this time APF gradients are used to obtain the desired velocity for any robot i as, ẋdi = −∇Uia − ∇Uir (5.1) where, ∇(·) = d(·) dxi Uia = X j (Uija ) Uir = X (Uijr ) . (5.2) j The order of the dynamic equation drops from second to first and, compared to Equation 4.4. In Equation 5.1 the virtual viscous damping is missing as the first order dynamics does not produce overshoots to damp out. Up to this point, relative position was used to produce the guidance law. The velocity field can be modified considering global translational and rotational terms. Here just a rotational one, Hi , is considered. It is defined as Hi = H(xi ) = kr K(xi − X) (5.3) where, X is the position of the centre of mass of the formation; K is a square matrix of dimension 2 that produces a rotation θ of the velocity vector. The value θ = π2 K is considered here to produce a velocity component in a direction orthogonal to the position vector of each robot with respect to the centre of mass of the formation, that is producing a pure rotation of the swarm. The matrix K is hence defined as " # " # cos(θ) − sin(θ) 0 −1 K= = . (5.4) sin(θ) cos(θ) 1 0 Matrix K produces a component of the desired velocity orthogonal to the position vector xi − X of each agent with respect to the centre of mass of the formation. A general form can include a parameter kr that tunes the magnitude of the tangential velocity. For the experiments presented in the proceeding, this was taken as unitary. Additionally the rotational term returns a general form of the guidance law as, ẋdi = −∇Uia − ∇Uir + Hi . (5.5) The process that leads a group of 5 agents to achieve a pentagon rather than a cross shape has already been described in Section 4.2.1. The group of 5 robots, CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS (a) 112 (b) Figure 5.5: A formation of 5 agents arranged in (a) pentagon shape due to all-to-all potential with the same coefficients for all the robots and (b) in a cross shape due to one agent being less “repulsive” with respect to the others. Attitude is not considered. Formation is defined just by the middle point of the wheel axle controlled through pairwise artificial potential, can arrange in two different cluster configurations that correspond to stable equilibria. Also here, by tuning one single parameter Lrij along the directed edges connecting one single agent to the other 4, the pentagon formation is no longer an equilibrium one. The formation is pictured in Figure 5.5 and it can be analysed considering the gradient of the artificial potential sensed by any of the robots as already done in Equations 4.9 and 4.10, for the initiator of the fractal formation presented in Chapter 4. For what concerns the stability of the control law, satisfying the condition Lrij < Laij guarantees the presence of a minimum (i.e. convexity property) at the equilibrium [32]. Null derivative can be obtained as a function of Cijr /Cija ratio. For agent 1, as done already in Chapter 4, exploiting the convenient reference frame of Figure 5.5(a), dropping the indices, the x-derivative becomes ∂U1 = ∂x1 pentagon       Ca d d2 = 2 a exp − a (cos α) + exp − a (cos β) L L L       r d d2 C −2 r exp − r (cos α) + exp − r (cos β) L L L where, (5.6) CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS d d2 = 2 s 1 tan α + cos α 2 + 1 = kd . 113 (5.7) If the repulsive scale distance considered just for agent 1 shrinks to a value Lr ′ such that Lr ′ < Lr the pentagon configuration is not an equilibrium one anymore. Indeed, the derivative of Equation 5.6 with respect to Lr ′ can be proved to be always negative for d > Lr ′ , as per Equations 4.13-4.15, or otherwise can be made always negative by setting correctly the other parameters of the equation. Again a reduction of the repulsive scale length produces an increase in velocity along the positive x − axis. When this is done for just one robot, it will collapse in the centre of the formation, which will rearrange in the cross shape. Moreover, as Equation 5.6 is linear with respect to C r , as pointed out previously, the ratio C r /C a can be used to scale the physical size of the formation. The Control Functions The low level controller uses the desired velocity output of APF to provide the rotational speed values for the wheels. This was developed considering the design of the wheeled robots used. The desired velocity vector components are mapped to the wheel mean speed (WMS) and the wheel speed difference (WSD). These are calculated, for a generic robot i, as function of the magnitude of the artificial potential derivative and of the error in heading, W MSi = |ẋdi | Sr Υ i |∇Ui∗ | W SDi = ∆θi Sr π (5.8) where, Sr is the maximum speed value allowed by the motor in the non-dimensional range [-100; 100], ∆θ is the error in heading measured as difference between the desired heading and the actual one, ∆θi = θid − θi θid = tan−1 ẏid ẋdi . (5.9) while θi is provided by the tracking system. U ∗ is the value of the artificial potential measured at a desired stand-off distance. The function Υi damps the magnitude of the speed commanded to the wheel in case of large heading error. CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS It is defined as Υi = Υ(∆θi ) = cos  ∆θi 2  . 114 (5.10) This function, positive semi-definite in the range [−π/2, π/2], weights the error in heading more than the ones in position. It ensures that the correct heading is acquired before producing a high WMS which might be demanded by the APF in case for large errors in position that correspond to large absolute value of artificial potential. The actual command to the speed is 2 element vector where the first element is the left wheel speed and the second one is the right wheel speed. It is defined as wi = " W MSi − W SDi /2 W MSi + W SDi /2 # . As the motors used do not output any torque in the interval ]-13;13[ a deadband was set imposing a rotational speed corresponding to level 13 for any input in the interval with exclusion of zero. Saturation of the motors can only be partially tackled by the scaling parameter U ∗ in Equation 5.2.2 that can be considered valid just for the case of 2 robots. When more robots are considered, it potential |∇U ∗ | is heuristically defined as |∇Ui∗ | = (N − 1)  Ca exp La  −d∗ La  Cr − r exp Li  −d∗ Lri  (5.11) where, just the index of Lr parameter is considered being the other parameters the same for all the agents. The constant d∗ = 500mm is a separation distance between two generic robots, considered isolated, and N is the number of robots. d∗ is chosen considering twice the maximum distance between axle centre (around which the robot rotates) and the farthest point of the robot chassis from this one. Finally the sensitivity to positioning errors is defined as the error on the magnitude of the artificial potential gradient for each agent. Once again in the case of 2 or 3 robots it is relatively easy to scale up the threshold on the potential gradient below which the group of robots is considered in a formation, that is, when the potential gradient magnitude is directly related to the error in relative positioning. When there are more than 3 robots, because any possible equilibrium does not produce the same separation between any two robots the computation of the threshold becomes more difficult. For these reasons the tolerance to relative 115 CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS positioning errors, that maps into tolerance on the APF gradient magnitude was heuristically scaled to account for a variable number of robots. The value of the threshold was taken in terms of potential gradient as ∇U(dd + dtoll ) where, dd is the design distance and δd is the required precision, or the tolerance. The associated potentia is ∆Uitoll N −1 = N +1  Ca exp La  −(dd + dtoll ) La  Cr − r exp Li  −(dd + dtoll ) Lri  . (5.12) The achievement of the desired spacing is obtained by solving Equation 4.18 for Cijr /Cija as a function of the other parameters, which are assigned, with the distance amongst the robots considered as a design parameter. The same values were used for parameter Laij , for all (i, j) to calculate Cijr /Cija . Lrij parameter was instead assumed to be the same for all the interactions amongst the robots, except for the interactions sensed by the robot taking the central position. As the testbed labels the robots consecutively (i.e. from 1 to 5), agent number 1 was given a reduced value for Lrij , that is Lr ′ . Values used, calculated using Equation 4.20 for a cross arm design size of 700 mm, are reported in Table 5.1. Table 5.1: Numerical values of coefficients used in numerical simulations referred to a cross formation with a arm of 700 mm. Ca 99.9996 Cr 100 La 700 Lr 698 Lr ′ =69.8 A series of tests were performed to validate the switching and the rotational behaviours that the system exhibits when the value of Lr parameter sensed by one agent drops to Lr ′ and when the coupling to the centre of mass of the formation is enabled through matrix H respectively. Here observations made are reported with two tests described in detail as examples. 5.2.3 Experimental Results In this section tests performed are described. The arrangement in a cross formation is tested and outputs for 3 requested positioning distances are illustrated, CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 116 beside the capability of rotating the formation or switching from a pentagon to a cross formation. Formation Acquisition and Static Positioning Accuracy The success in the achievement of a static formation is strongly dependent upon the level of accuracy in positioning that is requested. For the 700 mm cross case an accuracy of 20 mm was difficult to obtain. Relaxing the requirement to 30 mm acquiring a static formations took on average 324 seconds but this figure shrinks dramatically as soon as the precision requirement is relaxed further. To test the ability of achieving the cross formation within a prescribed tolerance, 30 tests were performed, 10 tests for each required tolerance on the inter-agent distance (30, 50 and 100 mm). The robots were given random initial positions within the testing area and zero initial velocities. As the final configuration is achieved by a minimization of the inter-agent potential that just depends on relative positions and the final positions with respect to the external reference frame are not assigned, the accuracy in getting a precise positioning must be mapped into a threshold value of the potential gradient sensed by each rover. Results obtained are summarised in Table 5.2. Table 5.2: Results of the static relative positioning tests for a cross formation with a 700 mm arm. Tolerance and error are referred to the design inter-agent distance Tolerance Average Error Error Variance Average time taken 30mm 18.3mm 364mm2 324s 50mm 32.7mm 1238mm2 107s 100mm 61.9mm 3965mm2 94s A critical issue was spotted in the frequency at which the system is updated, which for 5 agents and the testbed architecture described previously, was 10 Hz. This is ten times lower than the frequency at which position data are streamed by the Vicon; hence just one every ten measurements from Vicon was used to generate the commands to the actuators. As each command keeps on executing CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 117 till the new command is fed in, robots may end up passing from an error in positioning on one side to an error on the other side while trying to correct the first one. Moreover, the control frequency scales down linearly with the increase in the number of robots, while the number of local maxima and minima in the artificial potential field increases. Difficulties in the achievement of very precise positioning are also partially due to the noise in measurements and the misalignment between the centre of the robots perceived by the tracking system and the actual centre about which the robot rotates. This range of problems is closely linked to the testbed architecture which includes external tracking and centralised computation to mimic on-board intelligence. As such, these limitations can be overcome through improving the testing technology available but do not diminish the validity of the results obtained in terms of pattern formation using the APF method. Some issues can be spotted in the design of the robots and their suitability for the use in conjunction with the APF method. The limited mobility due to the differential drive is a key issue to account for when considering the capabilities of the system to quickly achieve a precise configuration. It was noted how, despite the cosine term in the low level control function, the final equilibrium position was missed several times when the rovers were committed to sharp bends in the final approach to the equilibrium positions. Actuators other than differentially driven wheels can be used to improve the design of the robots and make them more suitable for the APF control method. The deadband of the actuators, that is, the lower bound on the available spinning rate of the wheels, was found to contribute significantly, in negative sense, to the coarse precision in the final approach. The performance of the actuators should be compliant with the action requested by the controller, which, in turn, may be improved and optimised to bridge the APF velocity field with the achievable performance of the actuators. Formation Switching In the formation switching tests a group of 5 robots randomly deployed in the test area at the beginning of the test acquires the pentagon formation described, then, the Lr parameter is manually switched to Lr ′ through the GUI for the agent that the system labelled as nr 1 according to the automatic labelling introduced CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 118 400 200 Y [mm] 0 −200 −400 −600 −800 −2800 −2600 −2400 −2200 −2000 −1800 −1600 −1400 −1200 −1000 X [mm] Figure 5.6: Trajectories followed by the 5 agents while arranging in a pentagon formation first (stars) and then in a cross formation (circles). in Section 5.2.1. This collapses into the centre forcing the others to adjust their distance accordingly. The two formations are pictured in Figures 5.6 and 5.7. Rotating Formation and Failure Tolerance The rotating formation was tested starting from a random arrangement and both including the rotation since the beginning of the test, and switching the rotational term on when the formation was already stabilized in a cross shape. Just the cross shape was tested although, in theory, nothing would prevent the pentagon formation to undergo rotation by the same means. Results obtained are illustrated in Figures 5.8 and 5.9. Although the tracks of the robots in the rotating formation do not overlap exactly, the maximum distance between 2 trajectories is always below 40 mm. The loss of one agent was simulated as well in this context. When this happens, agents are sorted again and the one with Lrij = Lr ′ is identified always as the CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS (a) 119 (b) Figure 5.7: A formation of 5 agents arranged in (a) pentagon formation and then (b) in cross formation. The pictures were extracted from the video taken by one of the Vicon cameras. corresponding to the first sorted agent. The formation will then rotate about the agent newly labelled 1 that now presents the reduced value of repulsive scale distance. The rotation about the central agent is a consequence of the rotation about the centre of mass of the formation, which coincides with the position of the agent number 1. The pattern reconfiguration confirms that the emergence of the central symmetry is not dependent upon the number of agents involved. It is the equilibrium configuration that robotic agents achieve to compensate for the asymmetry of the potential with respect to one agent. 5.3 Real and Simulated Real World The effective implementation of APF techniques towards autonomous robotics is to be regarded as a step in the direction of increasing the technology readiness level (TRL) of this control method. This has been done in this chapter in a all-toall connection scheme, that is, each robot has access to the information about the state of every other robot and uses this information to compute the desired state and evolve towards it. However, this is not always possible, in particular when CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 120 200 0 Y [mm] −200 −400 −600 −800 −1000 −1200 −2800 −2600 −2400 −2200 −2000 −1800 −1600 −1400 −1200 −1000 X [mm] Figure 5.8: Trajectories followed by the 5 agents while arranging in a cross formation and rotating about its centre. Each robot is tracked with a different colour. The initial positions are marked by the diamonds while the final positions by circles. The picture shows the robustness to the loss of an agent. The robot tracked with the red line stops operating while the formation is rotating about the robot tracked in blue. The system then sorts again the agents forming a 3 point star rotating about its centre where the agent in black finds its place. CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS (a) 121 (b) Figure 5.9: A formation of 5 agents arranged in (a) four point star and (b) in a three point star that excludes an agent with a failure. The pictures were extracted from the video taken by one of the Vicon cameras. the number of robots increases significantly (see also [105, 165]). Findings related to incomplete or switching communication graphs should also be verified in an experimental environment. The main contribution given by the tests presented is the experimental validation of one particular emergent behaviour, that is the collapse of a given formation into one of central symmetry by non symmetric changes in the guidance potential. Testing similar behaviour in a not fully connected system is the next step for this research. This kind of multi-agent system can find application in structural inspection because of the advantages of having devices able to perform the task in an autonomous and flexible fashion. These may include hard to access areas, or areas where access may have consequences for the plant. Concerns may arise because of the complex geometries and tough conditions that are usually found in real world industrial environments, while in the experiments robots were moving on a plane, hazard-free surface. This is the main difference, and the step that separates the real world from a simulated real world. As such, to achieve a further improvement on the TRL ladder, this issue has to be considered. However, the tests aimed to reproduce emergent behaviours experimentally that had previously been seen in simulations, or in mathematical models. Thus, they CHAPTER 5. PATTERN CREATION - EXPERIMENTAL RESULTS 122 consolidate the control techniques developed while keeping them unattached from any particular operational theatre. Beside this, the test campaign included the evaluation of the effect that the failure of one robot can have on the whole system, which is often the case in real operational scenarios. Chapter 6 Fast Consensus and Manoeuvring In the two previous chapters the focus was on formations, here the point is remaining cohesive and exploiting the network of connections to boost manoeuvring performance. The content of this chapter, even more than the others, is inspired by the biological world. The starting point is represented by concepts well known amongst naturalists such as the “Chorus line hypothesis” and the “Trafalgar Effect”, previously mentioned in Section 2.7. This chapter turns the problem of consensus and fast manoeuvring in avian flocks into one of designing an efficient distributed control scheme for swarms. This, taking advantage of the swarm characteristics, drives the agents to consensus and to manoeuvre quickly in response to external stimuli. 6.1 The Need for Fast Manoeuvring The break through of autonomous systems, especially in defence and intelligence applications, imposes that complex manoeuvres, up to now left to the human ability, be performed autonomously. This applies to swarms as well and, even more, to swarms characterised by quick dynamics. As an example the DARPA F6 project aimed to produce a highly responsive space fractionated architecture, with the capability of quickly dispersing and regathering [15]. Works confirming the growing interest in the direction of the F6 project can also be found in the literature, see for example [166]. In the design of a control architecture able to produce fast responses, nature can be extremely inspirational. In avian flocks, birds manage to share informa123 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 124 tion about migratory routes by influencing each other direction of flight and at the same time they defend themselves efficiently from predators. In a flock the chances of spotting predators increase together with the capabilities of confusing them when attacking. The quicker the decision - manoeuvre - regathering chain is triggered and performed, the better the chances for the flock. As inspired by natural world, swarm engineering aims to turn the spontaneous behaviour of biological swarm, optimised by several years of Darwinian evolution, into design rules able to return the same or comparable behaviour for engineered swarming systems. Performance of such systems are strongly dependent on the capabilities of communicating and sharing tasks. Naturalists observed a number of dynamics contributing to the fast response in animal groups. Some of them are inspiring and worth mentioning. In 1981 Theherne and Foster isolated a dynamic in the group motion of marine insects able to disorient the approaching predators. They named this dynamic as the “Trafalgar Effect” [112]. Predators approaching a side of the group produce changes in the locomotory behaviour of the individuals at the periphery. These ones increase their activity with rapid turnings that disorient the attacker. This produces more, apparently random, encounters which are responsible for the information to be transmitted throughout the group. This kind of behaviour triggers an avoidance manoeuvre faster than the approaching speed of the predator preserving the insects from the attack. Information is then passed on a neighbouring basis, in the sense that, individuals aware of the danger alert neighbours through their movements. A few years after the study by Theherne and Foster, Wayne Potts postulated the even more interesting “Chorus line hypothesis” [113]. Through footage analysis of avian flocks, Potts noticed that the manoeuvres, including those to avoid predators, are initiated by individuals at the edges and propagated inward with increased speed. Flock’s members realise of the approaching manoeuvre wave and start propagating it before the wave reaches them, producing an increase in the propagation speed. The manoeuvre wave then travels throughout the flock up to four time faster than the speed achievable just by passing the information on neighbouring basis. This reduces the chances of the predators to catch any of the flock’s members. The propagation of manoeuvring waves was later found to be CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 125 unrelated to the flock size. More precisely, it exhibits scale-free characteristics, whereby the speed of propagation in flocks of larger sizes increases compensating for the increased sizes and hence requiring approximately the same time for a manoeuvring wave to be transmitted throughout the flock [110]. Transferring the capabilities of animal groups onto engineered systems raises the question: how is it possible to ensure that information flows quickly and sufficient attention is paid to the threats of the external world as well as the targets the system is designed for? The answer has to consider two aspects. The first one is translating the concept of “paying attention” typical of mindful creatures into something that can be used by devices with a very low-level of, or no intelligence at all. The second one is ensuring this task is carried out by the swarm efficiently as a whole, exploiting its own features. Flocking birds manoeuvre in response to predator attacks or follow the group motion initiated elsewhere in the flock; in both cases they respond to an external stimulus but, as different from monolithic systems, the group response emerges out of the responses of the single flock members to the stimulus and to the group motion itself. The idea of “paying attention” to something becomes then weighting the external drives against the group dynamics. In terms of engineered swarms, an external drive can easily be seen as a signal to follow, regardless this signal contains the route to the accomplishment of a task or is generated by an approaching threat. Providing all agents with the capability of tracking the external drives would mean not taking advantage of the group structure, that is having a number of independent agents, each one looking after itself only. On the other hand, exploiting the connections in the group it is possible to allocate the tracking duties more efficiently amongst all the swarm members. The amount of resources each individual has to provide for tracking is then strongly connected to the position of the agent within the system. This is only partially related to the physical position of the agents as in the case of marine insects through which the Trafalgar Effect was postulated; the position of an agent in the swarm has to be considered as its position in the network of agents. This will be discussed in the following sections. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.2 126 Fast manoeuvring out of Network Structure In Section 2.6.2 it was shown how particular networks make quick information spreading possible, and drive multi-agent systems towards consensus and coherent responses. When the spatial arrangement of the agents determines the network, as in the case of nearest neighbours networks, and the network determines the dynamics of the agents, the swarm is not guaranteed to evolve towards a configuration that boosts the consensus speed or the information spreading such as the small world network. The lack of control on network arrangement gives cause to find alternative methods to achieve quick consensus. The problem is then reversed: given any network of interaction, how can this be exploited to produce fast responses? In a swarm where all agents have similar characteristics, the network can be used to rank them according to their influence on the rest of the swarm. Consider an agent observed by a large number of swarm members. This agent will be highly influential as its dynamics influences a wide portion of the ensemble. Ideally this agent should be visible by its local neighbours (other agents at close distance) and, possibly, by some distant ones. This would push the swarm towards a small world arrangement and would allow the agent to influence both its local neighbourhood and more remote regions of the swarm. It would be clever then to provide this agent with more power to follow an external signal or the ability to spot threats for the whole group. Indeed, even without having such an advantageous network of connections, it is always possible to provide the most observed agents with more tracking capabilities. Several measurements are available in the literature that aim to define the degree of centrality of a node in a network. These include, though being not limited to, Betweenes Centrality, Closeness Centrality, Node Certainty [167, 168]. Although none of them is detailed in the following, the idea of leveraging a “reciprocal influence” measurement can be used to improve the swarm reactivity. In the next section it will be shown how algebraic network characteristics can be fed back into the dynamics to influence it directly. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.3 127 Eigenvectors and Eigenvalues in a Network System In Section 2.4 it was stated how the number of zero eigenvalues of the Laplacian indicates the number of components in the graph. If zero is a simple eigenvalue, then the graph is connected. It can be proved that the spectrum of the Laplacian matrix is positive semi-definite. Consider a complex, square matrix A of order P n, with entries aij , for i, j ∈ {1, . . . , n}. Define Ri = j6=i |aij | as the sum of the absolute values of the non-diagonal entries in the i-th row. For such a matrix a Gershgorin disc D(aii , Ri ) is defined as the closed disc on the Re − Im plane centred at aii with radius Ri . Consider then Gershgorin circle theorem: Theorem 1. Every eigenvalue of A lies within at least one of the Gershgorin discs D(aii , Ri ). The proof of the Gershgorin Circle Theorem can be found in many linear algebra textbooks, see for example [169], hence is here omitted. As the diagonal elements of the Laplancian are either positive or null, and each row, excluding the diagonal entry, sums up to the negated diagonal element, all the eigenvalues of the Laplacian matrix are either positive or null. As discussed in Section 2.4, the negated Laplacian can be used to model a multi-agent system in a linear fashion (see for example [76, 103, 104]). In this case the modelling takes the form ẋ = −Lx (6.1) As long as the zero eigenvalue is a simple eigenvalue, as state in Section 2.4, the graph is connected and this model guarantees the achievement of a consensus. In particular, the system converges to a consensus value as the negated Laplacian, considering the Gershgorin theorem, will have all its eigenvalues in the left part of the Re − Im plane, which includes the origin. The Re − Im plane for the eigenvalues of a negated Laplacian is illustrated in Figure 6.1 where, Gershgorin disks are centred on three diagonal entries of the Laplacian. As the zero eigenvalue always belong to the spectrum, the convergence speed of a linear system modelled through the negated Laplacian matrix is expressed by the second smallest eigenvalue in magnitude. The larger the magnitude of the second eigenvalue, the faster the system will converge. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 128 Figure 6.1: Eigenvalues of a negated Laplacian matrix according to the Gershgorin circle theorem. Every eigenvalue lies within at least one of the Gershgorin discs; these are centred on the matrix diagonal entry and have radius equal to the row sum. All the eigenvalues will have non-positive real part In a matrix the first left eigenvector, here identified by c, is the one associated with the largest eigenvalue in magnitude. For the Perron-Frobenius theorem this coincides with the spectral radius for positive semidefinite matrices [170]. A negated Laplacian matrix −L can be turned into a positive semi-definite matrix M by adding its spectral radius ρ(L) to the diagonal entries, i.e. M = −L + ρ(L)I (6.2) where, I is the identity matrix of appropriate dimensions. This way the whole spectrum is shifted in the positive half of the Re − Im plane and the zero eigenvalue coincides with the spectral radius. Consequently, the left eigenvector corresponding to the zero eigenvalue of the Laplacian becomes the first left eigenvector of the augmented Laplacian M, without undergoing any change as a constant diagonal perturbation does not change the eigenvectors. As a consequence of this, all the properties that apply to the dominant eigenvalue and its eigenvector of positive semidefinite matrices translate for the negated Laplacian for its zero eigenvalue and the corresponding eigenvectors (left and right) equally. This includes the properties obtainable as a consequence of the Perron-Frobenius theorem. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 129 For the Laplacian and the adjacency matrices, the eigenvector associated to the zero eigenvalue reveals some fundamental characteristics of the network. In particular, it is a measure of centrality for the nodes of the network as it quantifies the number of paths that pass through a given node compared to all possible paths. In order to understand this, consider the power method which can be used to approximate the first left eigenvector. Given a generic matrix A, the power method provides an approximation of its first eigenvector c, closer to the actual eigenvector as higher powers are applied to the matrix A. This is multiplied by the uniform vector w(0) = (1, 1, . . . , 1); at the l − th iteration the method returns w(l) = (1, 1, , 1)A(l) . (6.3) When A is the adjacency matrix, the first iteration w(1) returns the row sum of A, that is the outdegree of each node, i.e. w(1) = (1, 1, , 1)A = (d1 , d2, ...., dn ) . (6.4) At the second iteration the i − th entry of vector w(2) corresponds to the sum of the outdegrees of all the neighbours of node i, that is the number of possible walks of length 2 departing from node i. Propagating this, it is easy to understand how, for high l, w(l) provides the number of paths of any length departing from each node. As the power method provides an estimate of the Perron vector, it can then be concluded that the i − th entry of vector c is the property of node i to be the origin of a number of arbitrary long paths in the graph. An important consequence of this is the first left eigenvector ranks the nodes based on their ability to start a successful information spreading across the network. This characteristic goes beyond the use done here of the Perron vector, as applications can be easily targeted in the spread of epidemic diseases, fades or emerging behaviour in fields other than engineering [171, 172]. When the graph’s Laplacian matrix can be taken as a representation of the associated dynamical system, the characteristic of ranking the most influential nodes allows use of the first left eigenvector to define the final, steady state of the unperturbed system subject to some initial conditions. This is easily proved considering that the inner product of the Perron left eigenvector and the state vector is an CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 130 invariant of the system as Equation 6.5 shows: d(c · x) = c · ẋ = −cLx = 0 dt . (6.5) Hence, the expression c · x(0) = c · x(t) holds and considering a generic vector norm k · k, it follows lim x(t) = t→∞ c · x(0) (1, 1, ...1)T kck . (6.6) Furthermore, the Perron-Frobenius theorem reveals that all the entries of the first eigenvector are nonnegative, which will be useful in the following sections. Finally, it is worth noting how the zero row sum of the Laplacian causes the right eigenvector, call it e, corresponding to the zero eigenvalue, to be a uniform vector, that is e ∈ span(1) 6.4 . Distributing Tracking Capabilities In a swarm designed to accomplish some mission tasks, the amount of computational resources used for signal pursuing is to be considered relative to the total amount of resources on which the swarm relies to carry out its mission. Allocating these resources then becomes central for the accomplishment of the swarm tasks. The problem of signal pursuing carried out by a group of individuals connected through a network of sensing links is then analysed here. When the resources are limited, their allocation amongst the group members for signal pursuing has to be considered. In particular the performance of fast convergence in the direction of the stimulus is strongly influenced by the way the swarm pursues the stimulus. Consider a generic, linear, first order dynamics as ẋ = Ax + Bu (6.7) where, x is the system state, a stacked vector of order n containing single agent states, A ∈ Rn×n is the system matrix defining the linear system dynamics, that is how the states affect each other; u is the input vector of order m that contains CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 131 the external stimuli and B ∈ Rn×m is the input matrix that determines how the inputs affect the system states. In this case the relative importance of the external drives, compared to the group dynamics, is given by the relative weight on the system of the Ax and the Bu terms. Looking at the same system’s representation, the problem of tracking an external signal efficiently turns then into designing the input matrix B, defining which state is influenced by which input, or, thinking about a swarm system, which agent pays attention to which input and how much it does so compared to the other agents. For the purpose of this analysis u is constant as it is used to evaluate the response of the system to a simple perturbation. It is possible to associate u to a step input for the system and use it to characterise its performance. However, later on in this chapter, u will change with time becoming piecewise constant. The case of a continuous time changing reference signal u(t) is not considered in this chapter and considered towards possible future developments In multi-agent systems the linear relations are often expressed using the Laplacian matrix, as explained is Section 2.4. The Laplacian L hence will replace matrix A in Equation 6.7. Matrix B is then to be determined. Let’s assume that each agent in the swarm is able to track the same signal, although the strength sensed depends, for each agent, on the amount of resources this allocates to the tracking task. This is modelled using an input vector u ∈ span {1}, that is, u = u(1, 1, ....1)T which has dimension n and B ∈ Rn×n . In particular B can be taken diagonal, that is each state is influenced by one single entry of u, in this case the entries of B can be organised into an n-dimensional vector b which contains the signals for each of the n states. In a multi-agent system where the Laplacian is used to model inter-agent interactions and a driving signal from outside is fed to each agent, the dynamics of the i-th agent can then be modelled as ẋi = X j∈Ni (xi − xj ) + bi (u − xi ) (6.8) where, xi is the state of the i − th agent, u is the reference signal that the agents pursue with different strength, according to the gains bi . Finally, Ni is the set of neighbours of agent i in the graph, that is the set comprising of those 132 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING agents which agent i observes. Equation 6.8 considers an unweighted graph. The external signal u is compared to the actual state as each agent can only see the external signal as deviation from its actual state. For instance if the external signal is a heading direction to follow, then the agent can only measure the new heading starting from the actual one, hence as a deviation from that. Each agent i is able to track the signal as long as the corresponding bi is nonzero. In the next section the model in Equation 6.8 will be used to show how bi coefficients can be chosen to drive the swarm to quick consensus about the external signal, i.e. all the states will converge to u. 6.5 A Good Choice with Limited Resources The choice of the bi coefficients in Equation 6.8 determines the amount of effort each agent makes to track the signal u, as previously explained. Considering the role of the first left eigenvector in ranking the nodes, based on their reciprocal influence, a possible choice is leveraging the Perron vector of the Laplacian, i.e. the one associated with the zero eigenvalue. This is a positive vector as per the Perror-Frobenius theorem. The magnitude of c is determined by the way the eigenvectors are scaled, but, when fed into a dynamical system, it has to be compared with the magnitude of the underlying dynamics, that is, the Laplacian. A large magnitude of c implies a considerable amount of resources allocated globally to follow an external signal. Although the following arguments apply with any scaling of the eigenvectors, it is useful at this point to introduce the definition of L1 vector norm, also known as Frobenius norm. This will be important in Section 6.6. Definition 1. Given a generic vector v ∈ RN its L1 norm is defined as the sum of the absolute value of its components, that is, kvkF = kvk1 = N X i=1 |vi | . (6.9) At this point the model can be written in a more elegant way using vector notation as ẋ = −Lx + C(u − x) (6.10) CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 133 where, u is the signal which is the same for all the agents. C is a diagonal matrix whose entries are the entries of the first left eigenvector c of the Laplacian matrix L of the system, Equation 6.10 can be further transformed to a more suitable version as ẋ = −(L + C)x + Cu . (6.11) A change of coordinates can be performed to eliminate the term Cu from Equation 6.11. Consider y = x − (L + C)−1 Cu . (6.12) With u constant in time, the derivative of y is equal to the derivative of x. Isolating x from Equation 6.12 and substituting in Equation 6.11, this becomes ẏ = ẋ = −(L + C)(y + (L + C)−1 Cu) + Cu (6.13) where, most of terms cancel out turning the equation into ẏ = −(L + C)y . (6.14) This change of coordinates shifts the origin of the reference frame to the value of the signal u, that hence is associated with the value zero. The first left eigenvector allows the most influential nodes to be identified and allocated with a better signal pursuing capability. Moreover, it allocates no pursuing capabilities to those node that are not observed by any other nodes, and hence cannot contribute to make the system converge towards the signal. This is confirmed by the following theorem and lemmas. Theorem 2. Let L ∈ Rn×n be the Laplacian matrix of a connected graph and let the matrix C being defined as C = diag{c} with c being the first left eigenvector of L satisfying cT L = 0. Then the matrix −L̃ = −L − C is Hurwitz. In order to prove this theorem the following Lemmas are needed. Lemma 1. Let L be the Laplacian matrix of a connected digraph and let c be its first left eigenvector corresponding to the zero eigenvalue. Then c has a zero component ci if and only if node i has only outgoing edges or for each edge j − i entering node i with nonzero out-degree, cj = 0 holds. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 134 The lemma implies that all nodes i that correspond to ci = 0 are not globally reachable. Proof. Consider the i − th component of the vector cT L [cT L]i = ci Lii + X cj Lji = 0 . (6.15) j6=i By hypothesis the j = i term, that was taken out of the sum, is zero. For the Perron Froebenius theorem, elements of c are nonnegative, hence for the sum in Equation 6.15 to be zero, it must be cj = 0 for each j corresponding to an existing edge j − i. Considering the same equation, the reversed implication becomes straightforward. The lemma also implies that if a node is globally reachable, it must have ci 6= 0. Indeed, looking at Equation 6.15, a globally reachable node must have at least one incoming edge Lji corresponding to a cj 6= 0, hence the equation must be satisfied for some nonzero, hence positive, value of ci Lii . Lemma 2. Let L be the Laplacian matrix of a connected digraph and let c be its first left eigenvector corresponding to the zero eigenvalue. Then c has a zero component ci = 0 if and only if node i is not globally reachable. Proof. The first implication is a consequence of Lemma 1. More precisely, if ci = 0 Equation 6.15 can be satisfied either for Lji = 0 for all j 6= i, or for cj = 0 for Lji 6= 0. In the first case the node would not be reachable, i.e. observed, by any other node; in the second case all the nodes connecting to node j would have either no incoming connections or all connections belonging to nodes that correspond to zero components of the eigenvector. As there must be at least one nonzero component of c, that means that the nodes with cj = 0 are either not reachable or reachable amongst them, but not globally. For the second implication, in order to have non-globally reachable nodes, the graph must not be strongly connected, hence the Laplacian matrix L is reducible and, through permutations, it is possible to arrange it in the form L= " L0 B [0] LGR # (6.16) CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 135 where, L0 ∈ Rl×l and LGR ∈ Rk×k with n = l + k square matrices whose rows refer respectively to non-globally reachable and globally reachable nodes. [0] is a matrix of consistent dimensions whose entries are all zeros and B is a matrix of consistent dimensions with no relevant characteristics. Note that L0 is not a Laplacian, while LGR still is a Laplacian, being its row sum equal to zero. The eigenvalues of L are the union of the spectra of L0 and LGR , and because the graph is connected L has one single eigenvalue equal to zero that must belong to the spectrum of LGR as this is a Laplacian itself. Consider to partition the first left eigenvector accordingly, that is the first l and the second m entries as c= h cN GR cGR iT . (6.17) As cTN GR L0 = 0 must hold, and zero is not an eigenvalue for L0 , this means that cN GR = 0. Lemma 3. Let L be the Laplacian matrix of a connected digraph and let c be its first left eigenvector corresponding to the zero eigenvalue. If c has a zero component ci = 0 then the diagonal element of the Laplacian matrix Lii 6= 0. Proof. The lemma can be proved by contradiction by noting that if both ci = 0 and Lii = 0 then the node i would have no outgoing edges and all its incoming edges would start from non-globally reachable nodes, hence the node i would belong to an isolated component and the graph would be disconnected. The following examples show the importance of the lemmas stated above. Example 3. Consider the digraph already discussed in table 2.2 and here reported again in table 6.1 with its Laplacian matrix and its first left eigenvector. Table 6.1: Graph 2 example 2❡ ❅ ■ ❅ ❅ ❄ 3❡ ✲ ❡1 ❅ ❅ ❅ ❅ ❅ ✲ ❡4  0 0 0 0 −1 2 −1 0   L= 0 0 1 −1 0 −1 0 1    1      0 c= 0      0 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 136 Nodes 2, 3 and 4 are not globally reachable (they cannot be reached from node 1), but still they have both incoming and outgoing edges. The second, third and fourth entries of the first left eigenvector are then zeros. Example 4. Consider the digraph in table 6.2 with its Laplacian matrix and its first left eigenvector. Table 6.2: Graph 1 example 2 ❡✛ ❅ ■ ❅ ❅ ❄ 3❡ ❡1 ❅ ❅ ❅ ❅ ❅ ✲ ❡4  1 −1 0 0 0 1 −1 0   L= 0 0 1 −1 0 −1 0 1    0       0.33 c= 0.33      0.33 Reorienting the edge 2-1 makes vertices 2, 3 and 4 globally reachable, while node 1 is no longer globally reachable, nor simply reachable from any other node. This is reflected into the first left eigenvector which is scaled to 1 considering the L1 norm. It is now possible to prove theorem 2: Proof. First suppose L to be reducible and consider its partition as per Equation 6.16. The first l columns of L correspond to those nodes with index i for which ci = 0. All incoming edges to those nodes must have origin in nodes indexed by some j for which cj = 0 too. So all the nonzero entries of the first l columns must be in the top left partition L0 . As the zero eigenvalue is found in LGR , −L0 is Hurwitz. It follows that the matrix −L̃ = −L − C is Hurwitz too because the zero entries of c sum up the diagonal of L0 and the positive ones make the matrix −LGR − Ck Hurwitz, where Ck is the diagonal matrix composed of the last k components of c. The case of L being irreducible, can be considered as a particular one where l = 0 A consequence of Theorem 2 is shown through the following Theorem 3. Let L ∈ Rn×n be the Laplacian matrix of a connected graph and let c be its first left eigenvector, corresponding to the zero eigenvalue. Let the matrix D ∈ Rn×n be a diagonal matrix whose entries are the elements of the vector d CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 137 and are all nonnegative. Then the matrix −L∗ = −L − D is Hurwitz if and only if the scalar product hd, ci 6= 0. Proof. As both d, c are nonnegative, for their scalar product to be nonzero it must be that di , ci 6= 0 for at least one index i. If this is the case, considering the matrix L partitioned as in Equation 6.16, the nonzero elements of D would sum up the lower partition LGR making at least one row sum of it positive. The proof is completed considering the Theorem III in [173] and here reported as Theorem 4 with the notation used so far for sake of completeness. Theorem 4. ([173]) Let L ∈ Rn×n such that Lii > 0 and Lij 6 0 for j 6= i. Assume moreover that P Lii > j6=i Lij and that L is irreducible. The determinant of L then vanishes if P and only if nj=1 Lij = 0 for i = 1, 2, ..., n. Because the determinant does not vanish, the zero eigenvalue of the Laplacian disappears and Gershgorin circle theorem ensures that all the eigenvalues of −L∗ have then negative real part. The result of Theorem 3 can be explained considering that, whenever a driving signal is fed to one of the globally observable nodes, this will affect the whole group making it converge towards the signal. It can be concluded that a diagonal matrix, whose nonzero entries are the elements of the first left eigenvector, has the property of stabilising the system once summed up to the Laplacian of the graph. Up to now this chapter showed that the first left eigenvector of the Laplacian matrix is able to stabilise and drive to consensus a swarm of agents whose dynamics is modelled by the Laplacian itself. However, the effectiveness of it compared to other possible choices was not assessed. This assessment is produced in the next section. 6.6 The Significance of the First Left Eigenvector and the Effects of Scaling In a number of cases, the first left eigenvector represents the best possible choice of resource allocation over the nodes of the system for improving fast pursuing CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 138 of an external signal. The gain used to track the signal can be introduced in the form of a diagonal perturbation of the Laplacian. This is valid when the perturbation is of small magnitude with respect to the Laplacian. The amount of resources, which corresponds to the scaling of the perturbation, is a fundamental parameter. Consider again the system in Equation 6.18 ẏ = −(L + D)y (6.18) where, as different from Equation 6.14, matrix D represents a generic diagonal perturbation to the Laplacian. The spectrum of −(L+D) is shifted by a constant quantity k if a uniform diagonal perturbation kI is added, where I is the identity matrix of appropriate dimensions. This was already done in Section 6.3. Assume now that the matrix −(L + D) is Hurwitz and that the scalar k is positive and larger than the largest eigenvalue in magnitude of the −(L + D) matrix. This makes [−L − D + kI] a strictly positive matrix with the smallest eigenvalue in magnitude of −(L + D), call it λ1 , contributing to the spectral radius of [−L − D + kI], namely ρ([−L − D + kI]) = λ1 (−L − D) + k. The spectral radius received more attention in the literature than the smallest eigenvalue considered so far. In particular, given a generic matrix A ∈ Rn×n , the spectral radius is an always increasing function of each entry of A, and in particular the expression ∂ρ(A) ci ej = ∇ρ(A) = T ∂ij e c (6.19) holds where, e and c are, respectively, the first right and left eigenvectors as reported in [174]. In [175] an expression is given for the spectrum of a matrix subject to a small perturbation, that is λ(A + D) = λ(A) + cT De eT c . (6.20) This have important consequences. The first is that the spectral radius is an always increasing function of the matrix entries. Consider the matrix M = −L + kI: its spectrum is the spectrum of Laplacian matrix L shifted by the positive quantity k, and its eigenvectors are exactly the ones of L. This implies that, as long as a perturbation of small magnitude is introduced, the spectral radius ρ(M − D) = ρ(−L − D + kI) can be approximated by means of the first left and right eigenvectors of the Laplacian. This is detailed in the following subsection. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.6.1 139 Limited Resources and Small Perturbations The magnitude of the diagonal perturbation to the system dynamics can be measured in terms of vector norm. Consider a perturbation D = −C = −diag{c}: in particular let c and e be scaled such that their L1 norm is equal to some generic positive scalars γ and ǫ, that is, kck1 = γ kek1 = ǫ X X i i |ci | = γ |ei | = ǫ X X ci = γ i ei = ǫ . i Furthermore, assume that the system has high coupling gains, that is kLk ≫ kDk for some norm k·k. Then, for what was previously stated, P cT Ce γ 2 ǫ i c2i P ρ(M − C) = ρ(M) − T = ρ(M) − = ρ(M) − γkck22 e c γǫ i ci . (6.21) This implies that the smallest eigenvalue in magnitude of the perturbed negated Laplacian is λ1 (−L − C) = −γkck22 . (6.22) Note that here the Euclidean and the Frobenius norm are indicated with the classical p-norm notations to avoid confusion. An important consequence is that the first left eigenvector represents a gradient descending direction for the spectral radius at −L that is for the unperturbed system. This local property also holds in a neighbourhood of C = [0]. As the perturbation increases in magnitude, the minimum of the spectral radius is not said to be obtainable through a diagonal perturbation composed of the first left eigenvector. Considering a fixed Frobenius norm for the perturbation, the same argumentations used to find the diagonal matrix that maximizes the spectral radius in [174] can be adopted to find the minimum of it. Consider the gradient of the spectral radius in Equation 6.19 and consider the eigenprojection E of the generic matrix A as defined in [176]. Definition 2. The eigenprojection of a matrix A corresponding to the eigenvalue 0 or, for brevity, the eigenprojection of A, is the projection, i.e. the idempotent 140 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING matrix, Z such that R(Z) = N (Aν ) and N (Z) = R(Aν ) where, R indicates the Range and N the Null space. ν is the smallest k such that rank(Ak+1 ) = rank(Ak ) k ∈ {0, 1, 2....} . As the idempotent matrix is uniquely defined, so is the eigenprojection. It can be stated that the gradient in Equation 6.19 corresponds to the transpose of the eigenprojection E, that is E= ecT cT e (6.23) ceT . (6.24) cT e The following lemma shows how a diagonal perturbation ED composed of the ∇ρ(A) = E T = Perron vector attains the minimum of the spectral radius for the matrix (−L − ED + KI) amongst the diagonal perturbations with fixed Frobenius norm, in the hypothesis of small magnitude of it. Lemma 4. Let L be the Laplacian matrix of a directed graph on N nodes. Let e and c be respectively the right and left Perron vectors associated with the zero eigenvalue. Let K be a positive scalar greater than the largest eigenvalue in magnitude of L and ED be the diagonal matrix defined as ED = diag e ◦ c cT e (6.25) where “◦” indicates the product element by element. Then there exists a t∆ > 0 such that for all 0 < t ≤ t∆ the minimum of the spectral radius of (−L−t∆+KI) with ∆ belonging to the space of diagonal matrices of unitary Frobenius norm is achieved for ∆ = ED /kED kF . Proof. Consider the derivative of the spectral radius of the matrix M = −L + KI in the direction of ∆ and ED /kED kF , respectively 141 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING ρ (M − t∆) − ρ (M) t→0 t T = −E · ∆ = −trace (E∆)   ρ M − t kEEDDkF − ρ (M) (M) = lim t→0 t E D = −kED kF = −E T · kED kF ρ′∆ (M) = lim ρ′ED /kED kF (6.26) (6.27) where the inner product between matrices is defined as P · Q = trace (P Q) which provides a definition for the Frobenius norm of matrices as kP kF = Now, for the Cauchy-Schwarz inequality [174] trace(E∆) = kE T ∆k2F ≤ kEE T kF k∆T ∆kF = kEkF p trace (P T P ) (6.28) that implies ρ′ED /kED kF (M) ≤ ρ′∆ (M) . (6.29) Due to the definition of ED in Equation 6.25, and due to the fact that M has the same eigenvectors of L, e is always uniform and the gradient of the spectral radius depends upon the left eigenvector only. In the special case of a balanced graph, the Laplacian matrix is line sum symmetric, that is rows and columns sum up to the same values, though not necessary the same for every row or column. In this special case both the first left and the first right eigenvectors are uniform. Given a norm of the perturbation vector, the property of achieving the minimum of the spectral radius locally extends to all possible perturbation magnitudes, as the matrix M − C keeps the same eigenspace of M regardless of the magnitude of C. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.6.2 142 Fixed Trace, a Misleading Effect The following theorem and lemma, given and proved in [177], provide a useful means to find the best control matrix D for the system previously described, in the space of diagonal matrices with fixed trace. In particular, the theorem provides a way to calculate the matrix D with zero trace which returns the lowest spectral radius possible and the lemma extends the results to all possible traces. Note that when elements of D are all nonnegative, the trace corresponds to the L1 norm of the vector composed of the same nonzero entries of D. In order to make the following theorem and lemma understandable, the following notation is needed: - µ(A; t) is the minimum of the spectral radius for a generic matrix A subject to a diagonal perturbation with trace t. In particular for zero trace perturbation, it holds µ(A) , µ(A; 0). - Dq indicates a diagonal matrix whose nonzero entries are the elements of vector q. - The diagonal similarity of a matrix A is defined as A1 = Dq−1 ADq where q is a positive vector. The diagonal similarity of a nonnegative matrix A is nonnegative too and the value of µ(A) is not effected by the transformation. Theorem 5. Let A be an n-by-n essentially nonnegative matrix, and suppose that P is a permutation matrix such that P T AP is in Frobenius normal form. Let A0 be the direct sum of irreducible matrices that is obtained from P T AP by replacing all entries in off-diagonal blocks with 0’s, and let B be the line sum symmetric diagonal similarity of A0 . Then µ(A) = µ(A0 ) = (1/n)1T B1, and µ(A) = min{λmax (A + Dq ) : 1T q = 0} is achieved only by q = P [µ(A)1 − B1]. Lemma 5. Let A be an n-by-n essentially nonnegative matrix and let t be a scalar. Then µ(A; t) = µ(A) + t/n. Furthermore, if the minimum µ(A) is achieved only by the vector q with 1T q = 0, then the minimum µ(A; t) is achieved only by a vector q + t/n1 with 1 = (1, 1, 1....1)T . The calculation of the matrix B from Theorem 5 is not immediate. It consists in finding the positive diagonal matrix Dq that makes matrix A0 line sum symmetric. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 143 This corresponds to finding the solution of the system of equations Dq−1 Aq − Dq AT δ = 0 (6.30) where δ is a vector whose elements are the inverse of the elements of q. This can be solved through iterative procedures such as the Newton-Raphson method. The above theorem and lemma provide useful way to calculate the best diagonal perturbation for fixed trace. Unfortunately the fixed trace condition produces, in many cases, the effect of “fouling” the less influential nodes to allocate more influence to the most influential ones. As the arithmetic sum of the diagonal element must be the same, some elements, those corresponding to less observed nodes, can be taken negative so to allow other elements to achieve larger values, even beyond the amount which all the elements have to sum up to. This “mathematical trickery” does not match the initial intent of distributing resources in an efficient way as a negative element, that corresponds to passing a node wrong information about the driving signal, requires a resource allocation by the node anyway. It is then the absolute value of the elements, that corresponds to the actual amount of resources used, to be considered, and this is the case of L1 norm which the present chapter is mostly concerned on. 6.7 Solution of the First Order Linear Dynamics Numerical simulation for a network of 20 nodes connected in different ways were carried out and are reported here. The networks used are a non-symmetric lattice, a periodic lattice or ring, a random network and a Small World network. The system is of first order, of the kind described by Equation 6.10 and consensus speed is measured by the smallest eigenvalue in magnitude of the dynamical system matrix L + C or L + D, depending on the kind of diagonal matrix used to stabilise the system. Performances obtained through a uniform distribution of tracking capabilities and one corresponding to the first left eigenvector were compared with the speed of consensus obtained for a capability distribution obtained through numerical optimisation. Due to the continuous nature of the problem, the optimisation is performed using a gradient-based method implemented through the MATLAB R built-in optimiser fmincon. This performs the CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 144 gradient descending search using one of four algorithms automatically selected time to time depending on the problem’s peculiarities. Its objective is finding the diagonal perturbation D that maximises the smallest eigenvalue in magnitude of the matrix (L + D). The optimisation is constrained to the vectors of unitary L1 norm to produce a meaningful comparison with the first left eigenvector, which is scaled in the same way. Each node is seen as a mono dimensional variable “x” whose initial condition is zero and is pushed towards convergence to 1, being vector u in Equation 6.10 defined as u = (1, 1, ....1)T . The nodes are hence in agreement at the beginning of the simulation and converge towards the unitary signal in different fashions depending on how the tracking resources are allocated. In the following subsections the cases referred to the different networks considered are detailed. In subsection 6.7.5 a summary of the results is given, which includes a comparison of the performance achieved for the different networks with the numerical optimiser and the first left eigenvector. The comparison is made using the value achieved for the smallest eigenvalue of the perturbed Laplacian. 6.7.1 The Regular Lattice Case The lattice structure considered is non-symmetric in the sense that a node in general does not observe and is observed by the same subset of nodes, hence the corresponding graph adjacency and Laplacian matrices are non-symmetric. The generic node i “observes” the two in front of it, i + 1 and i + 2 and just one behind it, i − 1. Node 1 observes just nodes 2 and 3; the second-last node observe just the last node and the third-last while the last node just observes the second-last node. The asymmetry in the graph is wanted as a symmetric graph would produce a uniform first left eigenvector. This case is considered separately. A graphical representation of the lattice used is reported in Figure 6.2. This structure of the lattice graph finds justification in the preferential direction of view found in animal groups. The attention of individuals in large group focalises on a restricted group of individuals found in particular areas within their sensing range. Here it was chosen to give preference to the forward sensing more than the backward one within a maximum sensing range of 2 subsequent individuals, hence the asymmetry of the network. Moreover, a lattice is a plausible, CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 145 Figure 6.2: Sketch of the lattice network. Each node in the inner section has 3 outgoing and 3 incoming edges. Nodes at the left hand side have more outgoing than incoming edges, while nodes at the right hand side end have more incoming than outgoing edges. Note that an outgoing edge implies that the node where the edge is originated “observes” the node where the edge ends. though quite rigid, representation of the sensing links in groups where sensing range is limited. The values for the first eigenvalues achieved in the three cases considered for this graph are λ1 = −0.2994 when the first left eigenvector is used, λ1 = −0.3253 for the numerically optimised vector and λ1 = −0.05 for the case of uniform resource allocation. The difference between the performance using the optimised vector and using the first left eigenvector is very small compared to the consensus speed attained using a uniform distribution. In the latter case all nodes move rigidly towards the consensus value, but the motion is much slower than the one obtained with a non uniform allocation of resources. This is well represented in Figure 6.3. 6.7.2 The Ring Case The case of a periodic lattice builds upon the aperiodic lattice illustrated in the previous section. The difference consists in having agents at the extremities connected so to have the same number of inbound and outbound edges for each node, regardless its index. The resulting adjacency and Laplacian matrices are said to be “circulant”. Circulant matrices are balanced matrices and, as well as symmetric matrices, are a subset of normal matrices. In normal matrices the left and right eigenspaces present no difference and in particular the first left and right eigenvectors are the same. Thus, the first eigenvectors of the Laplacian for this network are both uniform. As shown in Section 6.6.1 the choice that maximise the speed of consensus, for a given magnitude, is obtained for a uniform distribution of the weights. The plots in Figure 6.4 illustrate the outcomes of the numerical simulation for the 20 node, first order, linear system for a ring network. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.7.3 146 The Small World Network Case The Small World network is often formed with a ring network with a number of rewired connections. Here the rewiring is kept but the original network is a nonperiodic lattice, as in Section 6.7.1. The rewiring is performed with probability psw = 0.1. The plots in Figure 6.5 show the dynamics of the 20 node, first order, linear system for a Small World network. 6.7.4 The Random Network Case Distinct from the three previous cases of lattice, ring and Small World networks, a random network does not have an underlying or base structure that is then modified. Each node connects to the others with a certain probability, here arbitrarily set to pr = 0.3. The resulting graph might be not connected but the results presented here do refer to a connected graph. The plots in Figure 6.6 show the dynamics of the 20 node, first order, linear system for a random network. The resource allocation presented in the figure is an example of how a numerical optimiser can encounter problems due to local minima. Because of this the first left eigenvector can surpass the performance of the optimised dynamics. However it should be noted how this is strongly dependent on the initial guess provided to the optimiser. Providing the first left eigenvector or a uniform vector as initial guess would have probably produced a different outcome. In particular, when the first left eigenvector is taken, the optimiser is likely to return the first left eigenvector as optimal choice. This is a consequence of the local optimal conditions guaranteed by this choice. 147 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 0 2 0.5 Optimised vector Uniform vector First left eigenvector 4 Resource allocation Node index 6 8 10 12 14 0.4 0.3 0.2 16 0.1 18 20 0 5 10 Node index 15 0 0 20 (a) 5 10 Node index 15 20 (b) 1 0.9 0.8 0.7 X 0.6 0.5 0.4 0.3 First left eigenvector Uniform vector Optimised vector First left eigenvector average Optimised vector average 0.2 0.1 0 0 5 10 t [s] 15 20 (c) Figure 6.3: Consensus in a lattice, performance comparison. (a) The adjacency matrix of the lattice where the dots represent nonzero entries. (b) Comparison between the values of the diagonal entries of matrix D for the matrix composed of the first left eigenvector, a vector obtained through numerical optimisation and one obtained by distributing evenly the tracking characteristics. All these vectors are scaled to have a unitary L1 norm. (c) Time evolution of the first order system driven to consensus about x = 1 by the diagonal perturbation of the Laplacian. 148 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 2 0.05 4 0.045 6 0.04 Resource allocation Node index 0 8 10 12 14 0.035 0.03 0.025 Optimised vector Uniform vector First left eigenvector 0.02 0.015 16 0.01 18 0.005 20 0 5 10 Node index 15 0 0 20 (a) 5 10 Node index 15 20 (b) 0.7 0.6 0.5 X 0.4 0.3 0.2 First left eigenvector Uniform vector Optimised vector 0.1 0 0 5 10 t [s] 15 20 (c) Figure 6.4: Consensus in a ring, performance comparison. (a) The adjacency matrix of the ring where the dots represent nonzero entries. (b) Comparison between the values of the diagonal entries of matrix D for the matrix composed of the first left eigenvector, a vector obtained through numerical optimisation and one obtained by distributing evenly the tracking characteristics. All these vectors are scaled to have a unitary L1 norm. (c) Time evolution of the first order system driven to consensus about x = 1 by the diagonal perturbation of the Laplacian. As the three distributions coincide the plots overlap. 149 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 0 2 0.5 Optimised vector Uniform vector First left eigenvector 4 Resource allocation Node index 6 8 10 12 14 0.4 0.3 0.2 16 0.1 18 20 0 5 10 Node index 15 0 0 20 (a) 5 10 Node index 15 20 (b) 1 0.9 0.8 0.7 X 0.6 0.5 0.4 0.3 First left eigenvector Uniform vector Optimised vector First left eigenvector average Optimised vector average 0.2 0.1 0 0 5 10 t [s] 15 20 (c) Figure 6.5: Consensus in a Small World network, performance comparison. (a) The adjacency matrix of the network where the dots represent nonzero entries. (b) Comparison between the values of the diagonal entries of matrix D for the matrix composed of the first left eigenvector, a vector obtained through numerical optimisation and one obtained by distributing evenly the tracking characteristics. All these vectors are scaled to have a unitary L1 norm. (c) Time evolution of the first order system driven to consensus about x = 1 by the diagonal perturbation of the Laplacian. 150 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 0.4 0 2 0.35 Optimised vector Uniform vector First left eigenvector 4 0.3 Resource allocation Node index 6 8 10 12 14 0.25 0.2 0.15 16 0.1 18 0.05 20 0 5 10 Node index 15 0 0 20 (a) 5 10 Node index 15 20 (b) 1 0.9 0.8 0.7 X 0.6 0.5 0.4 0.3 First left eigenvector Uniform vector Optimised vector First left eigenvector average Optimised vector average 0.2 0.1 0 0 5 10 t [s] 15 20 (c) Figure 6.6: Consensus in a random network, performance comparison. (a) The adjacency matrix of the network. The dots represent nonzero entries. (b) Comparison between the values of the diagonal entries of matrix D for the matrix composed of the first left eigenvector, a vector obtained through numerical optimisation and one obtained by distributing evenly the tracking characteristics. All these vectors are scaled to have a unitary L1 norm. (c) Time evolution of the first order system driven to consensus about x = 1 by the diagonal perturbation of the Laplacian. 151 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.7.5 Performance of First Order System The dominant eigenvalue can be taken as a measure of the speed of the system dynamics. For a system modelled through a Hurwitz matrix this eigenvalue is the smallest in magnitude. It is then possible to summarise the results obtained in terms of speed of response for the linear dynamics depending on the diagonal perturbation added to the Laplacian. This is done in table 6.3 for the cases analysed within this Section. Table 6.3: First eigenvalues of the system matrix (L+D) where, L is the network graph Laplacian and D is a diagonal matrix whose nonzero entries are either the first left eigenvector (F. L. E.), a uniform vector or a vector numerically optimised to maximise the correspondent eigenvalue. All of these have unitary L1 norm. D diagonal values Lattice Ring Small World Random F. L. E. 0.2994 0.05 0.2416 0.2430 Uniform 0.05 0.05 0.05 0.05 Optimised 0.3256 0.05 0.2466 0.0823 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 6.8 152 The Second Order and the Nonlinear Dynamics The dynamics observed for the first order system, where linear interactions represented by the Laplacian of the connected graph are augmented by the first left eigenvector of it, can be translated for a second order system where nonlinear relations can be put in place. This has the advantage of attaining a less abstract representation of the dynamical system composed of a swarm of mobile agents that manoeuvre in response to external stimuli while staying cohesive. Consider the second order dynamics: ẍi = Cid X j∈Ni (xj − xi − d) + Civ X j (ẋj − ẋi ) + Ciz (zi − xi ) + Ciu (ui − ẋi ) (6.31) where, the symbols are the same used for first order system; moreover Cid , Civ , Ciz , Ciu are constants that weight the relative importance of the terms of the equation, zi is a reference for the actual position xi and d is a desired gap between the positions of any two agents i and j connected in the network. Equation 6.31 can be used to model the dynamics of a number of agents aligned along a line. The difference between the positions of any two of them can be positive or negative depending on one preceding the other or viceversa. The difference xj − xi in the first term of Equation 6.31 is positive if the position of the j − th agent is further down the line than agent i. In this case it makes sense to compare it to the positive scalar d. Equation 6.31 can then be used as it is to model the dynamics of a number of forward-looking agents along the line. This corresponds to a line graph, either periodic (a ring) or not. This simple one-dimensional model can be used to show that some properties of the first order dynamics hold for the second order dynamics as well, where nonlinearities can be introduced to make the model better suited to reproduce the motion of a swarm of autonomous agents. First consider that each agent does not have a preference, or a driving input for its position in absolute terms, that is Ciz = 0∀i ∈ {1, ...N}. Equation 6.31 can be written in vector form as ẍ = −Cd Lx + d − Cv Lẋ + Cu u − ẋ (6.32) 153 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING where, C d , C v and C u are diagonal matrices and d is a uniform vector stating the desired distance each agent aims to keep with respect to its neighbours in the graph. This second order model can be translated into the state space becoming ( ẇ = y ẏ = (−C v L − C u )y − Cd Lw + Cu u − Cd Ld . (6.33) A coordinate change can be now operated to get rid of the terms related to the external input. Define two new variables ξ and θ such that ( ξ = w+d θ = y − (Cv L + Cu )−1 Cu u . (6.34) Using this transformation, Equation 6.33 becomes ( ξ̇ = θ θ̇ = −(C v L + C u )θ − Cd Lξ . (6.35) Finally Equation 6.35 can be turned into matrix form as ( ) " #( ) ξ˙ [0] I ξ = d v u θ̇ −C L −(C L + C ) θ (6.36) where, I is the identity matrix of appropriate dimensions. The entries of C d , C v and C u matrices are to be chosen consistently with the constant (unitary) total amount of resources the swarm uses to track an external signal and to produce a coherent swarm behaviour, that is Cid + Civ + Ciu = 1 ∀i ∈ {1, ...N} (6.37) For positive C d , C v and C u the system matrix is not Hurwitz as it keeps the zero eigenvalue but, for a connected swarm, this is a simple eigenvalue and all the other eigenvalues have negative real parts. This can be proved in the following Lemma 6. Let L be the Laplacian matrix of a connected directed graph on N nodes. Let C d , C v and C u be diagonal matrices with positive entries and in particular let C u have nonzero elements along the diagonal where the first left eigenvector, corresponding to the zero eigenvalue of the Laplacian matrix, presents 154 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING nonzero entries. Then all the eigenvalues of matrix S= " [0] I −C d L −(C v L + C u ) # (6.38) have non-positive real part. Moreover the spectrum of submatrix −(C v L + C u ) belongs to the spectrum of matrix S. Proof. Consider the characteristic equation for matrix S, that is det(S − λI) = −λI I d v −C L −(C L + C u + λI) . (6.39) As −λI and −C d L are permutable, from [178], the determinant can be calculated as det(S − λI) = |λI(C v L + C u + λI) + C d L| = 0 . (6.40) Consider to reduce all the matrices, which are not already diagonal, to their diagonal entries, that is their eigenvalues, and indicate these reduced matrices with subscript D, then det(S − λI) = |λI(C v LD + C u + λI) + C d LD | = 0 . (6.41) The characteristic equation reduces then to the determinant of a diagonal matrix. As this is given by the product of the diagonal terms, the eigenvalue problem can be treated as a scalar one and turned into finding the λs that make a generic diagonal element null. As there are only diagonal matrices, it is possible to easily multiply the terms in Equation 6.41 to obtain, for the generic i − th diagonal element, λ2i + (Civ Li + Ciu )λi + Cid Li = 0 , (6.42) which is a second degree equation. The λ that satisfies it can be found as λi = −(Civ Li + Ciu ) ± p (Civ Li + Ciu )2 − 4Cid Li 2 . (6.43) As Cid Li ≥ 0 the real part of the Equation 6.43 will always be non-positive. Moreover, the zero eigenvalue belongs to the spectrum of S as it belongs to the spectrum of C d L. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 155 The second part of the Lemma can be proved by considering the Schur formula [178] according to which the determinant det(S − λI) can also be calculated as det(S − λI) = | − λI − (C v L + C u + λI)C d L|| − C v L − C u − λI| = 0 (6.44) from which it is clear how the spectrum S includes the eigenvalues of −(C v L + C u ). The linear system is then stable and the presence of the zero eigenvalue reflects the fact that agents do not achieve consensus about a unique position as a consequence of C z = 0. However, they maintain distances with respect to each other. Extending the model to more dimensions, the design distance d has to be compared with a norm to map the relative positions between pairs of connected agents into a scalar. As an Euclidean norm is a nonlinear function of the relative positions, the extension beyond the one-dimensional model is made together with the introduction of nonlinear interactions in the network in the form of artificial potential functions. As the artificial potential functions represent a viable way to control relative positions in a highly nonlinear fashion they are suitable for the needs herein, and can be introduced considering the background outlined in the previous chapters. The introduction of nonlinearities is presented here although most results would require deeper analyses to be fully understood. In the following, results based on numerical tests are reported that open a wide horizon for prospective developments. Once more, the Morse potential is used and the model takes then form dx dt dv dt = v = −∇U − Lv + D(vdes − v) (6.45) where, ∇U is a time dependent vector representing the gradient of the Morse potential used in Chapters 3, 4 and 5. The numerical values for the coefficients of the potential are summarised in Table 6.4. The small magnitude of C a and C r coefficient is chosen not to let the inter-agent distance keeping dominate over the consensus to a common velocity. The signal to track is introduced into the equation as vdes while L is the graph Laplacian and D is a diagonal matrix containing the weights through which each agent tracks the signal. In Figure 6.7 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 156 it is shown how manoeuvring in a three-dimensional space is different depending on the values of D despite the fact that in both cases D has unitary Frobenius norm. The network used is generated at the beginning of any new manoeuvre as a directed graph in which each node observes the 5 nodes which are physically closer to it. This is similar to the connection rule described in Chapter 3, although the number of neighbours is significantly smaller, hence the connectivity of the graph cannot, in principle, be guaranteed. However, this choice is coherent with the structure of sensing networks in animal groups, as previously explained. As already verified for the first order case, when the diagonal of D is composed of the elements of the first left eigenvector, the swarm accelerates towards a reference direction quicker, covering longer distances in the same time. In Figure 6.7 this is compared to the case in which the weights are the same for all the agents, that is, the diagonal entries of D are all the same. The difference in performance for the change in heading can be better understood looking at the time history of the errors between the actual and the desired headings, shown in Figure 6.8. The error is calculated as the instantaneous angle ∆Θ between each agent’s velocity vector and the desired velocity vector. It can be seen how the error reduces quicker for the swarm where the tracking resources are allocated according to the first left eigenvector. Moreover, comparing panels (b) and (d) of Figure 6.8 it can be noticed how, for approximately the same change in heading, the swarm with uniform distribution reacts differently, increasing the gap in performance. This is most likely due to the increased inertia effects in panel (d) respect to panel (b) for the swarm with uniform distribution of the weights. This, in turn, is caused by an increased modulus of the velocity at which the change of heading is performed. In fact, the manoeuvre in panel (b) takes place at lower speed than in panel (d) due to the initial zero velocity condition, that, as understandable from the shorter distance covered (Figure 6.9), influenced more the swarm with uniform distribution of the weights. In the last manoeuvre the influence of the initial condition becomes negligible and the difference in performance appears more clear. The heading requested to the agents, after the initial static relaxation, are reported in Table 6.5. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 157 Table 6.4: Coefficients used in the APF for the three dimensional manoeuvring simulations. Ca 10−5 Cr 10−5 La 10 Lr 5 Table 6.5: Desired velocity vectors for the swarm to produce the changes in heading. The table shows the three components of the reference velocity vector which each agent is supposed to follow in the simulations shown in this section. 1 2 3 4 Vx Vy Vz √ √ √ 1 0 3 3 0 √ 2 2 3 3 0 √ 2 2 √ − 22 − 3 3 − √ 2 2 0 Avoiding the Undesirable Effects of Distance Keeping The use of nonlinear inter-agent forces for relative position keeping, that is the APFs in this case, was limited to preserve the group from dispersing. As the graph is not symmetric, so the inter-agent forces are, and this produces unbalanced drifts into the swarm of the same kind of the ones used in Chapter 4 to make one agent overcome the potential barrier and change the configuration of the swarm. This doesn’t always return a stable configuration as in Chapter 4, and the unbalanced interactions may produce an effect counteracting the desired convergence towards a common velocity. Increasing by a factor 102 the magnitude of the APFs, that is having C r = C a = 10−3 , with all the other parameters unchanged, introduces fluctuations in the velocity that reflect on the final trajectory and on the speed of manoeuvre. In Figure 6.9 the final trajectories followed by the swarm with the uniform tracking capability distribution (in blue) and the one with a tracking capability distribution reflecting the first left eigenvector are shown. It can be seen how dramatically the difference builds up. In Figure 6.10 instead the comparisons of the heading errors 158 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING Uniform vector First left eigenvector 10 Z [m] 5 0 −5 10 40 0 20 0 Y [m] −10 −20 X [m] Figure 6.7: Manoeuvring in three-dimensional space. Agents track a step changing external signal that provides the vectorial reference velocity. The agents sensing the signal weighted by the first left eigenvector (in red) achieve the velocity target sooner than the ones sensing the signal uniformly (in blue), hence cover more distance. are shown. The set of reference velocity used here is the same previously used and reported in Table 6.5. It can be noted how a uniform distribution of the weights produces bounces in the error time history. These are due to the effect of the asymmetries in the graph that generate asymmetric distance-keeping forces, and when the effect of these become predominant over the consensus to a common velocity, the consensus achievement slows down and the error increases again. When the consensus is sought through the first left eigenvector, the detrimental effect of cohesion forces looks to be minimised and a more coherent convergence to consensus emerge. This can be in part due to the fact that the first left eigenvector awards the nodes that are mostly observed with more tracking capability than the nodes that are mostly observers. As the observers tend to keep a given relative distance from the observed, the distribution provided by the first left eigenvector makes the observed node to drive the swarm relying on the distance keeping of the observers to follow up on the motion. However, these are just preliminary studies whose outcome look promising for the future expansion of the research branch here outlined. 159 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 60 150 Uniform vector First left eigenvector 100 40 ∆ θ [deg] ∆ θ [deg] Uniform vector First left eigenvector 50 50 30 20 10 0 20 25 30 t [s] 35 0 40 40 45 (a) 50 t [s] 55 60 (b) 60 100 Uniform vector First left eigenvector Uniform vector First left eigenvector 50 80 ∆ θ [deg] ∆ θ [deg] 40 60 40 30 20 20 0 60 10 65 70 t [s] (c) 75 80 0 80 85 90 t [s] 95 100 (d) Figure 6.8: Error time history with respect to the desired heading. Before starting the manoeuvres the agents relax for the first 20 seconds aiming to a null velocity. (a) The first manoeuvre is a straight acceleration that sets the agents in motion aiming to get unitary speed; (b) the first change in heading is performed with the swarm weighting the signal using the first left eigenvector (in red) surpassing the swarm were sensing is uniform distributed (in blue) in achieving the new heading; (c) a similar difference can be observed for the second change in direction; (d) the gap in performance increases in the last change in direction. Although the angular gap is similar to the one in case (a), the increased effects of the inertia penalise the swarm characterised by a uniform distribution more than the other. 160 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING Uniform vector First left eigenvector 5 Z [m] 0 −5 −10 10 0 5 −10 Y [m] −20 −5 10 15 0 X [m] Figure 6.9: Detrimental effect of inter-agent forces in manoeuvring. Agents track a step changing external signal that provides the vectorial reference velocity. The agents sensing the signal weighted by the first left eigenvector (in red) are less affected by the contrasting behaviour of keeping relative distances than the agents using a uniform weight distribution (in blue). 6.9 Fast Manoeuvring - Final Considerations The first left eigenvector can be used to create a diagonal control matrix that stabilises a multi-agent system with dynamics governed by the Laplacian matrix. In the case of a small magnitude of the control matrix with respect to the coupling strength, the first left eigenvector represents the optimal choice for maximising the magnitude of the system smallest eigenvalue, that is ensuring fast convergence towards an external input signal. This finds evidence into the local maximisation of the magnitude of the first eigenvalue of the system, that is the smallest in magnitude for a system with dynamics described by a Hurwitz matrix. In the case where the connection graph of the system is balanced, the first left eigenvector represents the best possible choice of allocating resources amongst the system’s nodes for tracking the external signal, regardless its magnitude. In this case the first left eigenvector belongs to span{1}. 161 CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 200 100 Uniform vector First left eigenvector Uniform vector First left eigenvector 50 ∆ θ [deg] ∆ θ [deg] 150 100 50 0 20 0 −50 −100 25 30 t [s] 35 −150 40 40 45 (a) 55 60 (b) 100 200 Uniform vector First left eigenvector Uniform vector First left eigenvector 150 80 100 50 ∆ θ [deg] ∆ θ [deg] 50 t [s] 0 −50 −100 60 40 20 −150 −200 60 65 70 t [s] (c) 75 80 0 80 85 90 t [s] 95 100 (d) Figure 6.10: Detrimental effect of inter-agent forces on the error time history with respect to the desired heading. Before starting the manoeuvres the agents relax for the first 20 seconds aiming to a null velocity. (a) The effect of the inter-agent forces can be noticed since the first acceleration: as the swarm with uniform distributed tracking gets close to the target velocity, some agents experience a braking effect that can be visualised as reversing the direction of motion of 180 ◦ C; (b) the first change in heading highlights the difference in behaviour between the two swarms; (c) a similar difference can be observed for the second change in direction; (d) the swarm driven to consensus through uniform weight distribution seems to produce a quick minimisation of the error but, it bounces back towards larger errors before achieving the target heading. CHAPTER 6. FAST CONSENSUS AND MANOEUVRING 162 When the system is not balanced nor the control matrix is of small magnitude, the first left eigenvector cannot be said to drive the best choice possible for tracking resource allocation. It still leads to a minimum of the first eigenvalue of the system though just in local sense. The global minimum can be analytically found for the set of diagonal control matrices with fixed trace although this has less applicability to a problem of resource allocation including negative entries in the control gains. Through linking the scaling of the eigenvector to the amount of resources available for pursuing a given external signal, it is possible to build up a meaningful model of resources allocation in distributed architectures. The fundamental problem of calculating the first left eigenvector in a distributed fashion, that is referred to a network for which each node does not have full visibility, can be overcome by using a distributed estimate as proposed by Qu et al. [179]. However, it is beyond the objectives of this dissertation to assess the performance of the swarm when distributed calculations are performed. Relying on the graph algebraic properties to improve the reactivity of a swarm is advantageous with respect to a numerical optimisation for the definition of the amount of control each node has to provide. Local minima can be avoided and results are close to the optimal ones, if not optimal. The extension towards second order nonlinear dynamics introduces a number of challenges to account for. This is just partially done here, outlining the possibilities to expand the research in this direction. It is important to stress how, relying on the communication network rather than on the agents own dynamics, this methodology is widely applicable across a number of distributed systems, which are not restricted to swarms of mobile agents. Information networks may get benefits through a quicker diffusion of new data, power grids can potentially be monitored against cascade failures in a more efficient way or computer virus infections spreading across the internet can be detected and targeted quicker by providing just some nodes with the required capabilities. Chapter 7 Conclusions and Future Work 7.1 A Summary of the Topics Covered In this Thesis a general introduction to a rather complex field such as swarm engineering was first provided through Chapter 1. The approach based on considering agents as particles was discussed in as well as this method is widely used throughout the Thesis and allows for obtaining results of wide applicability. The first chapter also stated also the targets this thesis aimed and introduced its development. The mathematical tools used and their application to the specific problem of distributed control were presented in Chapter 2 with a review of the previous works in the area of swarming featuring the elements common to this thesis presented in Section 2.7. This was done by connecting the present work to the inspiring biological fields. The way group behaviours found in nature can be achieved with the techniques presented in the previous chapters is introduced. Chapter 3 leads to the first of the results of this Thesis. The problem of fragmentation in particle swarms, which are more in general representative of agent swarms, was considered by finding the conditions that guarantee the swarm to be connected. Reversing the popular approach which requires a swarm to be connected to guarantee its stability, in Chapter 3 conditions that lead to swarm fragmentation are highlighted and the Cheeger constant is introduced as metric for the swarm fragility. A technique requiring minimal computational efforts to have a measure of the swarm fragility was finally presented. 163 CHAPTER 7. CONCLUSIONS AND FUTURE WORK 164 Chapter 3 also introduces the possibility of formation shaping offered by the reduction of links. These are more extensively exploited in Chapter 4 where a formation flying concept is designed around this emerging characteristic of the swarming system. It was shown how, by exploiting asymmetric potential functions and a particular communication network, the deployment of a fractal arrangement of agents can emerge out of the swarm. This arrangement finds application in the electromagnetic field as antenna array, which is provided with enhanced radiating properties by its peculiar shape. The asymmetry in the artificial potential functions was then exploited again in a test campaign aimed to validate the emerging behaviour on real hardware. This is described in Chapter 5 where, the performance of a group of wheeled robots are evaluated when controlled through the artificial potential functions aforementioned. This clearly goes in the direction of increasing the level of confidence in swarming systems towards real world applications. Lastly, a noteworthy point covered in this Thesis in the area of consensus achievement was discussed in Chapter 6. Manoeuvring performance of a swarm of mobile agents are analysed considering the amount of resources the swarm can allocate to tracking an external signal or staying in watch for threats. It was proved how a uniform distribution of resources is able to provide the swarm with good performance only in few cases. The results of Chapter 6 indicate a criterium for selection of the agents mainly deputed to signal following in a swarm depending on their position in the network. These results, inspired by flocking behaviours of birds, are here applied to the manoeuvres in swarms of mobile agents. However the extension of these results into more general area of networked systems was suggested considering the capital importance of the reaction speed to external inputs for interconnected systems. 7.2 Conclusions This dissertation has investigated the problems of modelling and controlling a swarm of agents in an effective and reliable way departing from the all-to-all CHAPTER 7. CONCLUSIONS AND FUTURE WORK 165 communication schemes and leveraging the communication network in a swarm to produce emergent behaviours. This has been done in a rigorous analytic way with the scope of addressing real problems of distributed control with effective and reliable tools. In particular this dissertation has provided the following novel contributions to knowledge. ♦ The problem of fragmentation in large groups of agents with reduced communication capabilities was linked to the formation of a bottleneck in the communication network. This reflects onto the final spatial arrangement for swarms controlled through APFs. A fundamental limit to the number of links each agent has to keep was identified and connected to the flow of a generic quantity across the network. The decreasing trend of the flow with the increase in the number of agents highlights how in large engineering multi agent systems the loss of cohesion can happen also when agents are supposed to keep an high number of connections. ♦ The possibilities disclosed in the area of formation shaping by the clever design of the communication network were explored and the application to space-based telecommunication arrays discussed. An original contribution to knowledge was provided in the exploitation of a fractal shaped satellite formation for the deployment of a distributed array. Calculations show how such a telecommunication architecture surpassed in performance both large monolithic antennas and distributed, regularly spaced, array lattices. ♦ The artificial potential function method was used to validate an emergent segregation behaviour and explore the implications of applying APF methods to real robotic agents. Although experimental tests match validations, these highlight a strong dependance of the performance on the particular actuators used. ♦ The fast manoeuvring performance in swarms were linked with the capability to cleverly exploit the connection network characteristics and allocate limited tracking resources across the swarm members. Quick consensus and consequent fast manoeuvring are achieved without the need for numerical optimisation. These four topics cover open problems in swarm engineering and provide original results for the development of this discipline. CHAPTER 7. CONCLUSIONS AND FUTURE WORK 7.3 166 Future Work Multi-agent systems, and swarms in particular, have to overcome a certain resistance to their affirmation which is mainly related to their autonomy. Having one single autonomous vehicle interacting with the real world means achieving a very high level of confidence in the autonomous controls. This is penetrating into the technological paradigm nowadays. Having a swarm of vehicles mutually interacting and then, as a whole, interacting with the real world, in a fully autonomous way, is probably something that will be hardly seen in the near future. Nevertheless, both the always increasing interests in this technology raised within the research environments, and the evolution of our society towards an always more interconnected world suggest this is the direction our technological paradigm will develop. Hence the need, for the scientific community, to increase the knowledge about the multi-agent systems getting ready for the time in which our technological level and social structure will be ready to welcome them. One of the paths that leads to this increase of technological awareness starts from this Thesis and develops throughout a number of steps which can be summarised as follows. A fundamental step to take consists in the increasing of the TRL of swarming systems by constructing and testing fully autonomous swarms of robots, UAV, (micro) spacecraft and so on. Initially this has to be done in an isolated environment. Swarm systems are being tested nowadays in laboratories but few of them are ready lo leave this ”simulated real environment” to then operate in the real world. Before this, swarm systems need to be tested in the absence of external tracking and far from human presence and intervention. This can be achieved by employing these systems in tasks that do not require the presence of humans in the neighbourhood. Autonomous mining, structural inspections, the exploration of remote regions, or space applications, as the ones suggested in this work, can be a valuable validation for technologies that, at some point will join our everyday life. The research about autonomy can progress through considering the possibility of evolving the control laws here presented with time varying potentials. This should include more than one emergent behaviour and the switching from one to the other in an autonomous way as well. Here the introduction CHAPTER 7. CONCLUSIONS AND FUTURE WORK 167 of an asymmetry in the artificial potential was proved to produce a switch from one configuration to the other. In the future having systems that are able to switch their configuration depending on the external stimuli, but in the same predictable and reliable way as it was done here, would significantly increase the level of autonomy. This of course is to be seen in a context even wider than formation shaping. A refinement to the theoretical results of this Thesis would be the introduction of noise in the driving signal and the capability of the swarm of filtering out the disturbance through sharing their understanding of the signal. To the author’s knowledge this development direction is being pursued in several research institutions and represents another steps towards meeting the challenges of the real world. A further step into the refining of the theoretical models presented may include the problems related with the discontinuities when the network changes. This would shift the problem into the field of nonsmooth analysis. The results presented here mainly exploited the heterogeneity of the interactions, that was used to rank and group the agents. Future swarming system are likely to include agents heterogeneous by their physical characteristics. Wheeled vehicles may interact with aerial ones and also include computing nodes with no capability of processing matter but just information as it was in the early Cellular Automata. The concept of distributed fractal antenna presented here requires almost continuous actuation to contrast the orbital dynamics and force the agents in frozen relative positions. Nevertheless if the emergent behaviour can be blended with the natural motion given by the orbital dynamics, a more refined system can be produced. Moreover the array produced can be improved significatively by considering three-dimensional formations. This would imply a deeper integration of the swarm dynamics with the operational scenario towards a more defined, although less open-ended, applica- CHAPTER 7. CONCLUSIONS AND FUTURE WORK 168 tion. Although the number makes swarming systems intrinsically safe from mission critical point of view, their interactions with humans has to be faced at some point. This will produce a shift from mission-critical to safety-critical in the design of swarming systems. Hence research has to be addressed in the direction of human-machine interaction, with the machine being something different from what it is nowadays meant. However, what the machine should be is not fully defined at present. The technology that produced self-driving cars looks to be ready for the public in the next future. Although this is not exactly an example of swarm engineering, yet it witness its feasibility or, at least, the will to produce the shift to the next technological paradigm. From fast flying cars visionary foresaw to self driving vehicles on the doorstep the paradigm has somehow already shifted from “speed” to “control”. Drawing physical devices out of the paradigm and making them available to anyone is likely to produce the next social shift. Bibliography [1] W. Weaver. Science and complexity. American scientist, (36):536–44, 1948. [2] A. Ilachinski. Land Warfare and Complexity, Part I: Mathematical Background and Technical Sourcebook. AD-a362 620. Center for naval analyses, Alexandria, VA, 1996. [3] I. Couzin. Collective minds. Nature, 445(7129):715, 2007. [4] J. Von Neumann. The general and logical theory of automata. Cerebral mechanisms in behavior, pages 1–41, 1951. [5] M. Gardner. The fantastic combinations of John Conway’s new solitaire game “life”. Scientific American, 223:120–123, October 1970. [6] G. Beni. The concept of cellular robotic system. In Intelligent Control, 1988. Proceedings., IEEE International Symposium on, pages 57 –62, aug 1988. [7] G. Beni. From swarm intelligence to swarm robotics. In Proceedings of the 2004 international conference on Swarm Robotics, SAB’04, pages 1–9, Berlin, Heidelberg, 2005. Springer-Verlag. [8] R. Beckers, S. Goss, J. L. Deneubourg, and J.M. Pasteels. Colony size, communication and ant foraging strategy. Psyche, (96):239–256, 1989. [9] T. Fukuda and S. Nakagawa. Approach to the dynamically reconfigurable robotic system. Journal of Intelligent and Robotic Systems, 1:55–72, 1988. [10] C. W. Reynolds. Flocks, herds and schools: A distributed behavioral model. SIGGRAPH Comput. Graph., 21(4):25–34, August 1987. 170 171 BIBLIOGRAPHY [11] A. F. T. Winfield, C. J. Harper, and J. Nembrini. Towards dependable swarms and a new discipline of swarm engineering. volume 3342, pages 126 – 142, 2005. [12] K. A. J. White, J. D. Murray, and M. A. Lewis. Wolf-deer interactions: A mathematical model. Proceedings of the Royal Society of London. Series B: Biological Sciences, 263(1368):299–305, 1996. [13] E. Ben-Jacob. Bacterial self-organization: co-enhancement of complexification and adaptability in a dynamic environment. Philosophical Transactions: Mathematical, Physical and Engineering Sciences, 361(1807):1283– 1312, 2003. [14] A. Moiseff and J. Copeland. Firefly Synchrony: A Behavioral Strategy to Minimize Visual Clutter. Science, 329(5988):181+, July 2010. [15] DARPA. Broad agency announcement system f6, 2010. [16] O. Brown, P. Eremenko, and B. A. Hamilton. Fractionated spacearchitectures: A vision for responsive space. In Proceedings of the 4th Responsive Space Conference. AiAA, 2006. [17] H. Kitano. Robocup rescue: a grand challenge for multi-agent systems. In MultiAgent Systems, 2000. Proceedings. Fourth International Conference on, pages 5 –12, 2000. [18] J. Saez-Pons, L. Alboul, J. Penders, and L. Nomdedeu. Multi-robot team formation control in the GUARDIANS project. Industrial Robot: An International Journal, 37(4):372–383, 2010. [19] N. Correll and A. Martinoli. A challenging application in swarm robotics: The autonomous inspection of complex engineered structures. Bulletin of the Swiss Society for Automatic Control, pages 15–19, 2007. [20] Atair Aerospace. Atair - onyx parachute system. http://www.atair.com/onyx/. Accessed 25 March, 2013. [21] C. Urmson and W. Whittaker. Self-driving cars and the urban challenge. Intelligent Systems, IEEE, 23(2):66 –68, march-april 2008. BIBLIOGRAPHY 172 [22] C. Bergenhem, E. Hedin, and D. Skarin. Vehicle-to-vehicle communication for a platooning system. Procedia - Social and Behavioral Sciences, 48(0):1222 – 1233, 2012. [23] S. H. Strogatz. Sync: The Emerging Science of Spontaneous Order. Hyperion, 1st edition, 2003. [24] J Goldstein. Emergence as a Construct: History and Issues. Emergence: Complexity and Organization, (1):1–9, 1999. [25] J. Gravner and D. Griffeath. Modeling snow-crystal growth: a threedimensional approach. Physical Review E, 75:1–18, 2009. [26] J.K. Parrish, S.V. Viscido, and D. Grnbaum. Self-organized fish schools: an examination of emergent properties. Biology Bulletin, 202(3):206–305, june 2002. [27] D. Peak, J. D. West, S. M. Messinger, and K. A. Mott. Evidence for complex, collective dynamics and emergent, distributed computation in plants. Proceedings of the National Academy of Sciences of the United States of America, 101(4):918–922, January 2004. [28] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes. SelfPropelled Particles with Soft-Core Interactions: Patterns, Stability, and Collapse. Physical Review Letters, 96(10), 2006. [29] H. Levine, W. Rappel, and I. Cohen. Self-organization in systems of selfpropelled particles. Phys. Rev. E, 63:017101, Dec 2000. [30] T. Vicsek, A. Czirók, E. BenJacob, I. Cohen, and O. Shochet. Novel type of phase transition in a system of self-driven particles. Phys. Rev. Lett., 75:1226–1229, Aug 1995. [31] V. Gazi and K.M. Passino. Stability analysis of swarms. Automatic Control, IEEE Transactions on, 48(4):692 – 697, april 2003. [32] C. R. McInnes. Vortex formation in swarms of interacting particles. Phys. Rev. E, 75:032904, Mar 2007. BIBLIOGRAPHY 173 [33] F. Schweitzer, W. Ebeling, and B. Tilch. Statistical mechanics of canonicaldissipative systems and applications to swarm dynamics. Phys. Rev. E, 64:021110, Jul 2001. [34] Werner E. Nonequilibrium statistical mechanics of swarms of driven particles. Physica A: Statistical Mechanics and its Applications, 314(14):92 – 96, 2002. [35] H. Van Dyke Parunak and S. Brueckner. Entropy and self-organization in multi-agent systems. In Proceedings of the fifth international conference on Autonomous agents, AGENTS ’01, pages 124–130, New York, NY, USA, 2001. ACM. [36] W. A. Wright, R. E. Smith, M. Danek, and P. Greenway. A generalisable measure of self-organisation and emergence. In Proceedings of the International Conference on Artificial Neural Networks, ICANN ’01, pages 857–864, London, UK, 2001. Springer-Verlag. [37] F. L. Lambert. Shuffled cards, messy desks, and disorderly dorm rooms examples of entropy increase? nonsense! Journal of Chemical Education, 76(10):1385, 1999. [38] D. J. Bennet. Pattern Formation in Swarming Systems using Bifurcating Potential Fields. PhD thesis, 2010. [39] W. Ren and R.W. Beard. Consensus seeking in multiagent systems under dynamically changing interaction topologies. Automatic Control, IEEE Transactions on, 50(5):655 – 661, may 2005. [40] R. P. Wiegand, M. A. Potter, D. A. Sofge, and W. M. Spears. A generalized graph-based method for engineering swarm solutions to multiagent problems. In Proceedings of the 9th international conference on Parallel Problem Solving from Nature, PPSN’06, pages 741–750, Berlin, Heidelberg, 2006. Springer-Verlag. [41] A. Jadbabaie, J. Lin, and A.S. Morse. Coordination of groups of mobile autonomous agents using nearest neighbor rules. Automatic Control, IEEE Transactions on, 48(6):988 – 1001, june 2003. 174 BIBLIOGRAPHY [42] L. Moreau. Stability of multiagent systems with time-dependent communication links. Automatic Control, IEEE Transactions on, 50(2):169 – 182, feb. 2005. [43] W. Ebeling and U. Erdmann. Nonequilibrium statistical mechanics of swarms of driven particles. 8(4):23–30, 2003. [44] A. Martinoli, A.J. Ijspeert, and F. Mondada. Understanding collective aggregation mechanisms: From probabilistic modelling to experiments with real robots. Robotics and Autonomous Systems, 29(1):51 – 63, 1999. [45] E. Sahin, T. H. Labella, V. Trianni, J. Deneubourg, P. Rasse, D. Floreano, L. Gambardella, F. Mondada, S. Nolfi, and M. Dorigo. Swarm-bot: Pattern formation in a swarm of self-assembling mobile robots. volume 4, pages 145 – 150, Yasmine Hammamet, Tunisia, 2002. [46] K. Lerman, A. Martinoli, and A. Galstyan. A review of probabilistic macroscopic models for swarm robotic systems. volume 3342, pages 143 – 152, Santa Monica, CA, United states, 2005. [47] C. R. McInnes. Autonomous ring formation for a planar constellation of satellites. Journal of Guidance, Control, and Dynamics, 18(5):1215 – 1217, 1995. [48] D. J. Bennet and C. R. McInnes. Verifiable control of a swarm of unmanned aerial vehicles. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 223(7):939–953, 2009. [49] D. J. Bennet and C.R. McInnes. Spacecraft formation flying using bifurcating potential fields. In Proceeding of the 59th International Astronautical Congress, Glasgow, Scotland, UK, 2008. [50] D. J. Bennet and C. R. McInnes. Pattern transition in spacecraft formation flying via the artificial potential field method and bifurcation theory. In 3rd International Symposium on Formation Flying, Missions and Technologies, European Space Agency, (Special Publication) ESA SP, ESA-ESTEC, The Netherlands, April 2008. European Space Agency. BIBLIOGRAPHY 175 [51] M. H. Mabrouk and C. R. McInnes. Nonlinear stability of vortex formation in swarms of interacting particles. Physical Review E (Statistical, Nonlinear, and Soft Matter Physics), 78(1), 2008. [52] M.H. Mabrouk and C.R. McInnes. Wall following to escape local minima for swarms of agents using internal states and emergent behaviour. London, England, UK, 2008. [53] M. M. Mabrouk, C. W. Murray, K. Johnstone, and C.R. McInnes. Internal agent states: experiments using the swarm leader concept. [54] A. Badawy and C. McInnes. On-orbit assembly using superquadric potential fields. Journal of Guidance, Control and Dynamics, 31(1):30–43, 2008. [55] A. Badawy and C.R. McInnes. Small spacecraft formation-flying using potential functions. Acta Astronautica, 65(11-12):1783–1788, 2009. [56] D.F. Gordon, W.M. Spears, O. Sokolsky, and L. Insup. Distributed spatial control, global monitoring and steering of mobile agents. In Information Intelligence and Systems, 1999. Proceedings. 1999 International Conference on, pages 681 –688, 1999. [57] W.M. Spears, R. Heil, D.F. Spears, and D. Zarzhitsky. Physicomimetics for mobile robot formations. In Autonomous Agents and Multiagent Systems, 2004. AAMAS 2004. Proceedings of the Third International Joint Conference on, pages 1528 –1529, july 2004. [58] T. R. Smith, H. Hanmann, and N. E. Leonard. Orientation control of multiple underwater vehicles with symmetry-breaking potentials. In 40th IEEE Conference on Decision and Control (CDC), December 4, 2001 - December 7, 2001, volume 5 of Proceedings of the IEEE Conference on Decision and Control, pages 4598–4603. Institute of Electrical and Electronics Engineers Inc., 2001. [59] N. E. Leonard and E. Fiorelli. Virtual leaders, artificial potentials and coordinated control of groups. In Decision and Control, 2001. Proceedings of the 40th IEEE Conference on, volume 3, pages 2968–2973 vol.3. BIBLIOGRAPHY 176 [60] R. Sepulchre, D. A. Paley, and N. E. Leonard. Stabilization of planar collective motion: All-to-all communication. IEEE Transactions on Automatic Control, 52(Compendex):811–824, 2007. [61] R. Sepulchre, D. A. Paley, and N. E. Leonard. Stabilization of planar collective motion with limited communication. IEEE Transactions on Automatic Control, 53(Compendex):706–719, 2008. [62] S. Nair and N. E. Leonard. Stable synchronization of rigid body networks. Networks and heterogeneous Media, 2(4):597–626, 2007. [63] E. Fiorelli, N. E. Leonard, P. Bhatta, D. Paley, R. Bachmayer, and D. M. Fratantoni. Multi-auv control and adaptive sampling in monterey bay. 2004 IEEE/OES Autonomous Underwater Vehicles, pages 134 – 147, 2004. [64] N. E. Leonard, D. A. Paley, F. Lekien, R. Sepulchre, D. M. Fratantoni, and R. E. Davis. Collective motion, sensor networks, and ocean sampling. Proceedings of the IEEE, 95(Compendex):48–74, 2007. [65] N. E. Leonard, D. A. Paley, R. E. Davis, D. M. Fratantoni, F. Lekien, and F. Zhang. Coordinated control of an underwater glider fleet in an adaptive ocean sampling field experiment in monterey bay. Journal of Field Robotics, 27(6):718–740, 2010. [66] L. Euler. Solutio problematis ad geometriam situs pertinentis. Commentarii academiae scientiarum Petropolitanae, 8:128–140, 1741. [67] N. Biggs, E. Lloyd, and R. Wilson. Graph Theory, 1736-1936. Oxford University Press, 1986. [68] P. Erds and A. Rnyi. On random graphs, i. Publicationes Mathematicae (Debrecen), 6:290–297, 1959. [69] P. Erdos and A. Renyi. On the evolution of random graphs. Publ. Math. Inst. Hung. Acad. Sci, 5:17–61, 1960. [70] P. Erdos and A. Rnyi. On the strength of connectedness of a random graph. Acta Mathematica Hungarica, 12(1):261–267, 1961. [71] A. Barabasi and R. Albert. Emergence of scaling in random networks. Science, 286(5439):509–512, 1999. BIBLIOGRAPHY 177 [72] D. J. Watts and S. H. Strogatz. Collective dynamics of ‘small-world’ networks. Nature, 393(6684):440–442, 1998. [73] N. E. Leonard, G. Young, K. Hochgraf, D. Swain, A. Trippe, W. Chen, and S. Marshall. In the dance studio: Analysis of human flocking. In American Control Conference (ACC), 2012, pages 4333 –4338, june 2012. [74] H. Yamaguchi and T. Arai. Distributed and autonomous control method for generating shape of multiple mobile robot group. In Intelligent Robots and Systems ’94. ’Advanced Robotic Systems and the Real World’, IROS ’94. Proceedings of the IEEE/RSJ/GI International Conference on, volume 2, pages 800 –807 vol.2, sep 1994. [75] A. Rahmani, M. Ji, M. Mesbahi, and M. Egerstedt. Controllability of multiagent systems from a graph-theoretic perspective. SIAM J. Control Optim., 48(1):162–186, February 2009. [76] R. Olfati-Saber, J.A. Fax, and R.M. Murray. Consensus and cooperation in networked multi-agent systems. Proceedings of the IEEE, 95(1):215 –233, jan. 2007. [77] M.A. Lewis and G.A. Bekey. The behavioral self-organization of nanorobots using local rules. In Intelligent Robots and Systems, 1992., Proceedings of the 1992 lEEE/RSJ International Conference on, volume 2, pages 1333 –1338, jul 1992. [78] T. Balch and R.C. Arkin. Behavior-based formation control for multirobot teams. Robotics and Automation, IEEE Transactions on, 14(6):926 –939, dec 1998. [79] L. Blazovics, K. Csorba, B. Forstner, and C. Hassan. Target tracking and surrounding with swarm robots. In Engineering of Computer Based Systems (ECBS), 2012 IEEE 19th International Conference and Workshops on, pages 135 –141, april 2012. [80] L. Bayindir and E. Sahin. A review of studies in swarm robotics. Turkish Journal of Electrical Engineering and Computer Sciences, 15(2):115 – 147, 2007. BIBLIOGRAPHY 178 [81] Y. Mohan and S.G. Ponnambalam. An extensive review of research in swarm robotics. pages 140 – 145, Coimbatore, India, 2009. [82] E. Bahceci, Soysal O., and E. Sahin. A review: Pattern formation and adaptation in multi-robot systems. Technical report, Robotics Institute, Carnegie Mellon University, Pittsburgh, PA, Oct 2003. [83] Z. Bubnicki. Modern Control Theory. Springer, Berlin Heidelberg, 2006. [84] D. P. Atherton. Stability of Nonlinear Systems. Research Study Press - John Wiley & Sons Ltd, 8 Willian Way, Letchworth, Herts SG6 2HG, England, 1981. [85] V Lakshmikantham, S. Leela, and A. A. Martynyuk. Stability Analysis of Nonlinear Systems. Marcel Dekker, Inc., 270 Madison Avenue, New York, New York 10016, 1989. [86] S. Sastry. Nonlinear Systems: Analysis, Stability and Control. SpringerVerlag, New York Berlin Heldelberg, 1999. [87] A. M. Lyapunov. The general problem of the stability of motion. International Journal of Control, 55:531–534, 1992. [88] B. Carre. Graphs and networks. Clarendon Press, 1979. [89] R. Gould. Graph theory. Benjamin/Cummings Pub. Co., Menlo Park, Calif., 1988. [90] J. A. Bondy and U. S. Murty. Graph Theory. Graduate texts in Mathematics, 244. Springer, New York, 2008. [91] N. Biggs. Algebraic Graph Theory. Cambridge University Press, 2nd edition, 1993. [92] F. H. Clarke. Nonsmooth analysis and control theory. Graduate texts in mathematics 178. Springer, New York, 1998. [93] A. Filippov. Differential Equations with Discontinuous Right-Hand Side. American Mathematica Society Translations, (2):191–231, 1964. [94] A. Filippov. Differential Equations with Discontinuous Right-Hand Side. Kluver, Boston, MA, 1988. BIBLIOGRAPHY 179 [95] D. Shevitz and B. Paden. Lyapunov stability theory of nonsmooth systems. In Proceedings of the 1993 ASME Winter Annual Meeting, November 28, 1993 - December 3, 1993, volume 53 of American Society of Mechanical Engineers, Dynamic Systems and Control Division (Publication) DSC, pages 53–60. Publ by ASME. [96] H. G. Tanner, A. Jadbabaie, and G. J. Pappas. Flocking in fixed and switching networks. Automatic Control, IEEE Transactions on, 52(5):863– 868, 2007. [97] H. G. Tanner, A. Jadbabaie, and G. J. Pappas. Stable flocking of mobile agents - part ii: Dynamic topology. 42nd Ieee Conference on Decision and Control, Vols 1-6, Proceedings, pages 2016–2021 6583, 2003. [98] P. M. Morse. Diatomic molecules according to the wave mechanics. ii. vibrational levels. Phys. Rev., 34(1):57–64, 07 1929. [99] J. E. Jones. On the Determination of Molecular Fields. II. From the Equation of State of a Gas. Royal Society of London Proceedings Series A, 106:463–477, October 1924. [100] B. Nabet and N. E. Leonard. Shape control of a multi-agent system using tensegrity structures. In 3rd IFAC Workshop on Lagrangian and Hamiltonian Methods in Nonlinear Control, LHMNLC 2006, July 19, 2006 - July 21, 2006, volume 366 LNCIS of Lecture Notes in Control and Information Sciences, pages 329–339. Springer Verlag. [101] X. F. Wang and G. Chen. Pinning control of scale-free dynamical networks. Physica A: Statistical Mechanics and its Applications, 310(3-4):521–531, 2002. [102] K.A. Hawick and H.A. James. Small-world effects in wireless agent sensor networks. International Journal of Wireless and Mobile Computing, 4(3):155 – 164, 2010. [103] R. Olfati-Saber. Ultrafast consensus in small-world networks. volume 4, pages 2371 – 2378, Portland, OR, United States, 2005. [104] F. Yu, D. Lv, and X. Liu. Optimal distributed consensus in small world network. In 25th Chinese Control Conference, CCC 2006, August 7, 2006 - BIBLIOGRAPHY 180 August 11, 2006, 2006 Chinese Control Conference Proceedings, CCC 2006, pages 333–336. Inst. of Elec. and Elec. Eng. Computer Society. [105] D. Pais, M. Cao, and N. E. Leonard. Formation shape and orientation control using projected collinear tensegrity structures. In 2009 American Control Conference, ACC 2009, June 10, 2009 - June 12, 2009, Proceedings of the American Control Conference, pages 610–615. Institute of Electrical and Electronics Engineers Inc. [106] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic. Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study. Proceedings of the National Academy of Sciences, 105(4):1232–1237, 2008. [107] J. K. Parrish, S. V. Viscido, and D. Grnbaum. Self-organized fish schools: An examination of emergent properties. Biol. Bull, 202:296–305, 2002. [108] A. Cavagna, I. Giardina, A. Orlandi, G. Parisi, and A. Procaccini. The starflag handbook on collective animal behaviour: 2. three-dimensional analysis. Animal Behaviour, 76(1):237–248, 2008. [109] A. Cavagna, I. Giardina, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic. The starflag handbook on collective animal behaviour: 1. empirical methods. Animal Behaviour, 76(1):217–236, 2008. [110] A. Cavagna, A. Cimarelli, I. Giardina, G. Parisi, R. Santagati, F. Stefanini, and M. Viale. Scale-free correlations in starling flocks. Proceedings of the National Academy of Sciences, 107(26):11865–11870, 2010. [111] A. Procaccini, A. Orlandi, A. Cavagna, I. Giardina, F. Zoratto, D. Santucci, F. Chiarotti, C. K. Hemelrijk, E. Alleva, G. Parisi, and C. Carere. Propagating waves in starling, sturnus vulgaris, flocks under predation. Animal Behaviour, 82(4):759–765, 2011. [112] J. E. Treherne and W. A. Foster. Group transmission of predator avoidance behaviour in a marine insect: The trafalgar effect. Animal Behaviour, 29(3):911–917, 1981. BIBLIOGRAPHY 181 [113] W. K. Potts. The chorus-line hypothesis of manoeuvre coordination in avian flocks. Nature, 309(5966):344–345, 1984. [114] B. L. Partridge. The structure and function of fish schools. Sci. Amer., 246(6):114–123, 1982. [115] T. J. Pitcher. Some ecological consequences of fish school volumes. Freshwater Biology, 10(6):539–544, 1980. [116] TJ Pitcher, BL Partridge, and CS Wardle. A blind fish can school. Science, 194(4268):963–965, 1976. [117] J. R. Hunter. Communication of velocity changes in jack mackerel (trachurus symmetricus) schools. Animal Behaviour, 17(Part 3):507–514, 1969. [118] J. K. Parrish. Re-examining the selfish herd: are central fish safer? Animal Behaviour, 38(6):1048–1053, 1989. [119] A. E. Magurran and A. Higham. Information transfer across fish shoals under predator threat. Ethology, 78(2):153–158, 1988. [120] P. K. Visscher and S. Camazine. Collective decisions and cognition in bees. Nature, 397(6718):400–400, 1999. 10.1038/17047. [121] M. R. Myerscough. Dancing for a decision: a matrix model for nest-site choice by honeybees. Proc Biol Sci, 270(1515):577–82, 2003. Myerscough, Mary R England Proceedings. Biological sciences / The Royal Society Proc Biol Sci. 2003 Mar 22;270(1515):577-82. [122] P. Lucic and D. Teodorovic. Transportation modeling: an artificial life approach. In Tools with Artificial Intelligence, 2002. (ICTAI 2002). Proceedings. 14th IEEE International Conference on, pages 216–223, 2002. [123] F. Cucker and S. Smale. Emergent behavior in flocks. Automatic Control, IEEE Transactions on, 52(5):852–862, 2007. [124] I. D. Couzin, J. Krause, N. R. Franks, and S. A. Levin. Effective leadership and decision-making in animal groups on the move. Nature, 433(7025):513– 516, 2005. 10.1038/nature03236. BIBLIOGRAPHY 182 [125] I. D. Couzin. Collective cognition in animal groups. Trends in Cognitive Sciences, 13(1):36–43, 2009. [126] N. O. Handegard, K. M. Boswell, C. C. Ioannou, S. P. Leblanc, D. B. Tjstheim, and I. D. Couzin. The dynamics of coordinated group hunting and collective information transfer among schooling prey. Current biology : CB, 22(13):1213–1217, 2012. [127] B. Nabet, N. E. Leonard, I. D. Couzin, and S. A. Levin. Dynamics of decision making in animal group motion. Journal of Nonlinear Science, 19(Compendex):399–435, 2009. [128] L. Zhiyun, B. Francis, and M. Maggiore. Necessary and sufficient graphical conditions for formation control of unicycles. Automatic Control, IEEE Transactions on, 50(1):121–127, 2005. [129] V. Gazi. Stability of an asynchronous swarm with time-dependent communication links. IEEE Transactions on Systems, Man, and Cybernetics, Part B: Cybernetics, 38(Compendex):267–274, 2008. [130] R. Olfati-Saber and R. M. Murray. Consensus problems in networks of agents with switching topology and time-delays. Automatic Control, IEEE Transactions on, 49(9):1520–1533, 2004. [131] G. Punzo, J. Simo, D. J. Bennet, and M. Macdonald. Characteristics of swarms on the edge of fragmentation. under review in Physical Review E. [132] G. Punzo, D. J. Bennet, and M. Macdonald. Swarm shape manipulation through connection control. In 12th Conference Towards Autonomous Robotic Systems 2010, August 2010. [133] E. Stump, A. Jadbabaie, and V. Kumar. Connectivity management in mobile robot teams. In Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on, pages 1525–1530, May 2008. [134] M.M. Zavlanos and G.J. Pappas. Controlling connectivity of dynamic graphs. In Decision and Control, 2005 and 2005 European Control Conference. CDC-ECC ’05. 44th IEEE Conference on, pages 6388–6393, December 2005. BIBLIOGRAPHY 183 [135] F. Chung. Spectral Graph Theory (CBMS Regional Conference Series in Mathematics, No. 92). American Mathematical Society, Providence, RI, 02940-6248, USA, 1997. [136] F. Chung. Laplacians and the cheeger inequality for directed graphs. Annals of Combinatorics, 9(1):1–19, 2005. [137] G. Punzo, P. Karagiannakis, D. J. Bennet, M. Macdonald, and S. Weiss. Enabling and exploiting self-similar central symmetry formations. IEEE Transaction on Aerospace and Electronic Systems, 2013 in press. [138] G. Punzo, D. J. Bennet, and M. Macdonald. A fractally fractionated spacecraft. In 62nd International Astronautical Congress 2011, pages Article IAC–11–D1.4.11, October 2011. [139] G. Punzo. Fractal patterns in fractionated spacecraft. In 62nd International Astronautical Congress 2011, pages Article IAC–11–E2.2.6, October 2011. [140] G. Punzo, D. J. Bennet, and M. Macdonald. Enhancing self-similar patterns by asymmetric artificial potential functions in partially connected swarms. In Towards Autonomous Robotic Systems, volume 6856 of Lecture Notes in Artificial Intelligence. Institution of Engineering and Technology, September 2011. [141] J. Leitner. Formation flying: The future of remote sensing from space. Technical report, NASA, Goddard Space Flight Center, MD, United States, 2004. [142] E. A. Barnes and B. E. Bishop. Utilizing navigational capability to control cooperating robotic swarms in reconnaissance-based operations. Proceedings of the 40th Southeastern Symposium on System Theory, pages 338–342, March 2008. [143] S Schwartz. Cross-scale - multi-scale coupling in space plasmas. Technical report, ESA, December 2009. [144] M. Bester, M. Lewis, B. Roberts, J. McDonald, D. Pease, J. Thorsness, S. Frey, D. Cosgrove, and D. Rummel. Themis operations. Space Science Reviews, 141(1-4):91–115, December 2008. BIBLIOGRAPHY 184 [145] P. R. Lawson, O. P. Lay, S. R. Martin, C. A. Beichman, K. J. Johnston, W. C. Danchi, R. O. Gappinger, S. L. Hunyadi, A. Ksendzov, B. Mennesson, R. D. Peters, D. P. Scharf, E. Serabyn, and S. C. Unwin. Terrestrial planet finder interferometer 2006-2007 progress and plans - art. no. 669308. Techniques and Instumentation for Detection of Exoplanets III, 6693:69308–69308, September 2007. [146] K. Danzmann. LISA - an ESA cornerstone mission for the detection and observation of gravitational waves. Advances in Space Research, 32(7):1233 – 1242, 2003. [147] F. Sjoberg, A. Karlsson, and B. Jakobsson. PROBA-3: A formation flying demonstration mission on the verge to be realised. Number 654 SP, April 2008. [148] G. Blackwood, O. Lay, A. Deininger, W.and Ahmed, R. Duren, C. Noecker, and B. Barden. The starlight mission: A formation-flying stellar interferometer. Proceedings of SPIE - The International Society for Optical Engineering, 4852(2):463 – 480, February 2002. [149] C. R. McInnes. Velocity field path-planning for single and multiple unmanned aerial vehicles. Aeronautical Journal, 107(1073):419–426, July 2003. [150] M. J. Bentum, C. J. M. Verhoeven, and A. J. Boonstra. Olfar - orbiting low frequency antennas for radio astronomy. In Proceedings of the ProRISC 2009, Annual Workshop on Circuits, Systems and Signal Processing, Veldhoven, pages 1–6, November 2009. [151] C. Balanis. Antenna theory: analysis and design. New York: Wiley, 1997. [152] D. H. Johnson. Array signal processing:concepts and techniques. Englewood Cliffs, NJ : P T R Prentice Hall, 1993. [153] D. H. Werner and R. Mittra. The theory and design of fractal antenna arrays. In Frontiers in Electromagnetics, chapter 3, pages 94–203. WileyIEEE Press, 1 edition, 2000. BIBLIOGRAPHY 185 [154] C. Puente-Baliarda and R. Pous. Fractal design of multiband and low sidelobe arrays. Antennas and Propagation, IEEE Transactions on, 44(5):730– 739, May 1996. [155] A. Maffett. Array factors with nonuniform spacing parameter. Antennas and Propagation, IRE Transactions on, 10(2):131 –136, March 1962. [156] H. Schaub and J.L. Junkins. Analytical mechanics of space systems. American Institute of Aeronautics and Astronautics, Reston, Va., 2003. [157] L.B. King, G.G. Parker, S. Deshmukh, and J. Chong. Study of interspacecraft coulomb forces and implications for formation flying. Journal of Propulsion and Power, 19(3):497 – 505, May/June 2003. [158] M. A. Peck, B. Streetman, Saaj C. M., and V. Lappas. Spacececraft formation flying using lorentz forces. Journal of the British Interplanetary Society, 60:263–267, 2007. [159] G. Punzo, G. Dobie, D. J. Bennet, J. Jamieson, and M. Macdonald. Lowcost, multi-agent systems for planetary surface exploration. In 63rd International Astronautical Congress, October 2012. [160] S. Pierce, G. Punzo, G. Dobie, R. Summan, C. N. MacLeod, C. McInnes, J. Biggs, M. Macdonald, and D. J. Bennet. Reconfigurable robotic platforms for structural health monitoring. In 6th European Workshop on Structural Health Monitoring & 1st European Conference of the Prognostics and Health Management (PHM) Society, July 2012. [161] D. Izzo and L. Pettazzi. Autonomous and distributed motion planning for satellite swarm. Journal of Guidance, Control, and Dynamics, 30(2):449– 459, March/April 2007. [162] Vicon t160 website, http://www.vicon.com/products/t160.html. [163] G.I. Dobie. Ultrasonic Sensor Platforms for Non-Destructive Testing. PhD thesis, University of Strathclyde, 2010. [164] S.G. Pierce, k. Worden, R. Summan, G. Dobie, and J.J. Hensman. Towards implementation of reconfigurable robotic strategies for structural health monitoring. In 5th European Workshop on Structural Health Monitoring EWSHM, pages 610 –615, Sorrento, Italy, July 2010. BIBLIOGRAPHY 186 [165] Y. Liu, J. Slotine, and A. Barabasi. Controllability of complex networks. Nature, 473(7346):167–173, 2011. 10.1038/nature10011. [166] S. Nag and L. Summerer. Behaviour based, autonomous and distributed scatter manoeuvres for satellite swarms. Acta Astronautica, 82(1):95–109, 2013. [167] L. C. Freeman. A Set of Measures of Centrality Based on Betweenness. Sociometry, 40(1):35–41, March 1977. [168] I. Poulakakis, L. Scardovi, and N.E. Leonard. Node certainty in collective decision making. In Decision and Control (CDC), 2012 IEEE 51st Annual Conference on, pages 4648–4653, 2012. [169] R. S. Varga. Matrix Iterative Analysis. Springer-Verlag, Berlin, 2009. [170] D. Serre. Matrices: theory and applications. Graduate texts in mathematics; v. 216. Springer, 2010. [171] Z. Qu, C. Li, and F. Lewis. Cooperative control based on distributed estimation of network connectivity. In American Control Conference (ACC), 2011, pages 3441–3446. [172] K. Klemm, M. Serrano, V. M. Eguı́luz, and M. S. Miguel. A measure of individual role in collective dynamics. Scientific Reports, 2, February 2012. [173] O. Taussky. A recurring theorem on determinants. The American Mathematical Monthly, 56(10):pp. 672–676, 1949. [174] L. Han, M. Neumann, and M. Tsatsomeros. Spectral radii of fixed frobenius norm perturbations of nonnegative matrices. SIAM J. Matrix Anal. Appl., 21(1):79–92, October 1999. [175] G. W. Stewart. Introduction to matrix computations [by] G. W. Stewart. Academic Press New York, 1973. [176] R. P. Agaev and P. Yu. Chebotarev. On determining the eigenprojection and components of a matrix. Autom. Remote Control, 63(10):1537–1545, October 2002. BIBLIOGRAPHY 187 [177] C. R. Johnson, D. P. Stanford, D. D. Olesky, and P. van den Driessche. Dominant eigenvalues under trace-preserving diagonal perturbations. Linear Algebra and Its Applications, 212-213:415–435, 1994. [178] F.R. Gantmacher. Matrix theory, Volume 1. Chelsea Publishing Company, New York, NY, 1960. [179] Z. Qu, J. Wang, and R. A. Hull. Cooperative Control of Dynamical Systems With Application to Autonomous Vehicles. Automatic Control, IEEE Transactions on, 53(4):894–911, 2008.