arXiv:1609.06283v1 [cs.SY] 20 Sep 2016  Robotic Swarm Control from Spatio-Temporal Specifications  Iman Haghighi, Sadra Sadraddini, and Calin Belta  Abstract — In this paper, we study the problem of controlling a two-dimensional robotic swarm with the purpose of achieving high level and complex spatio-temporal patterns. We use a rich spatio-temporal logic that is capable of describing a wide range of time varying and complex spatial configurations, and develop a method to encode such formal specifications as a set of mixed integer linear constraints, which are incorporated into a mixed integer linear programming problem. We plan trajectories for each individual robot such that the whole swarm satisfies the spatio-temporal requirements, while optimizing to- tal robot movement and/or a metric that shows how strongly the swarm trajectory resembles given spatio-temporal behaviors. An illustrative case study is included.  I. I NTRODUCTION  Robotic swarms have received a lot of attention from the robotics research community in recent years [1]. Large teams of robots are suitable for a broad range of applica- tions such as distributed task allocation [2], area patrolling and coverage [3], [4], search and rescue missions [5], and simultaneous localization and mapping (SLAM) [6]. With the recent technological developments, producing a large number of inexpensive robots that are equipped with sophisticated sensing, computation and communication tools has become a reality. Describing complex spatial specifications for swarms is a non-trivial task. The existing methods rely on spatial config- urations generated from simple geometrical shapes, potential fields or sets of target points [7], [8], [9], [10]. However, it is practically easier to specify collective spatial behaviors of a swarm as opposed to specifying trajectories for each individual robot. The authors in [11], [12] introduced a method for controlling the abstract behavior of swarms based on the first and second moments of their spatial distribution. This is a useful approach to specify some simple patterns such as ellipsoids and boxes. However, there is a necessity for a more powerful framework of pattern specification that is not only easily definable and interpretable by the user, but is also rich enough to capture a wide range of complex spatial patterns that are not expressible by merely statistical moments. For this reason, we propose to use a formal spatial logic [13] that is capable of describing high level global behaviors in multi agent systems.  This work was partially supported by the National Science Foundation under grants NRI-1426907 and CMMI-1400167. I. Haghighi is with the Division of Systems Engineering, Boston Univer- sity, Boston MA, USA   haghighi@bu.edu  S. Sadraddini is with the Department of Mechanical Engineering, Boston University, Boston MA, USA   sadra@bu.edu  C. Belta is with the Department of Mechanical Engineering, Boston University, Boston MA, USA   cbelta@bu.edu  • • • • • •• • • •  Fig. 1.   An example of spatio-temporal patterning requirements for a swarm: While avoiding the unsafe zone (red) at all times, attain the following formations in any order within 30 seconds: 1) form a checkerboard pattern (green) by populating every other cell on the north east quadrant of the workspace. 2) Populate one of the grey squares in the south west quadrant. After completing both tasks, gather in one of the L shaped upload regions (cyan).   All the tasks must be completed   within 40 seconds   and each formation must be maintained for at least 3 seconds.  In this paper, we consider square-shaped two dimensional workspaces that are gridded to equal-sized cells as illus- trated in Fig. 1. A user can express spatial requirements for the swarm by defining shapes that can be formed by unions of cells in the grid. The user can give the swarm choices between distinct patterns. Furthermore, they can also specify certain requirements for how these patterns should evolve over time. An example of such a spatio-temporal specification is given in Fig. 1. Such specifications involve logical reasoning and provide different choices for the swarm movement. A   wide   variety   of   complex patterns can   be defined in this framework that are not easily expressible by earlier work in the literature. However, specifications of this type can be naturally expressed as spatial temporal logic (SpaTeL) [14] formulas. We formulate the swarm motion planning problem corre- sponding to a SpaTeL specification as a mixed integer linear programming (MILP) problem. A feasible solution to this problem provides a high-level plan for the movements of the swarm. We are also able to optimize a cost function. For instance, we are able to minimize the total swarm movement (energy), or find a plan that maximizes the satisfaction of a specification (robustness), or a combination of both. Finally, we develop a low-level strategy to move each individual robot according to the high level plan. Two different solutions for the specification given in Fig. 1 corresponding to two different initial conditions are given in Fig. 5 and 6. It can be seen that the optimal solution in Fig. 5 involves populating
the grey region before forming the checkerboard pattern and populating the right cyan region at the end, while the optimal solution in Fig. 6 forms the checkerboard pattern first and populates the left cyan region. Although temporal logic specifications in the context of mobile robot control and motion planning have been recently explored in the literature, there is very limited prior work in which complex requirements are expressed in   both   space and time. The authors in [12] attempt to solve a similar problem, but spatial specifications are limited to statistical moments of the swarm in their work and thus complex spatial patterns are not easily expressible. The authors in [15] introduced a   procedure to specify emergent spatial behaviors in swarms by linear temporal logic and used model checking techniques to verify such behaviors in swarms, but the control problem is not discussed in that work. As opposed to linear temporal logic multi-robot motion planning [16], [17], [18], [19], our solution is optimization-based which is advantageous in the following ways. First, we are able to optimize a general cost function, which is difficult to formalize in automata-based approaches. We are also able to deal with infeasibility by minimizing the   distance   of the swarm trajectory from satisfaction. Furthermore, under some relaxations, the complexity of our approach is independent of the size of the swarm. Therefore, our approach is easily applicable to large swarms. This paper is organized as follows. First, we provide the necessary background on SpaTeL specifications in Sec. II. Next, the problem is formulated in Sec. III. The solution and the technical details are explained in Sec. IV. Finally, an illustrative case study is presented in Sec. V. II. P RELIMINARIES  A. Quad Transition Systems  A   quad transition system   (QTS) [13], [14] is a tree data structure defined as the tuple   Q ( t ) = ( V ,   E , v ι , V f   , μ,   L , l ) , where:   V   is the set of nodes (vertices).   E ⊂ V × V   is the set of directed edges (transitions). We say that   v 2   is a   child  of   v 1   if and only if   ( v 1 , v 2 )   ∈ E ;   v ι   is the root (the only node which is not a child of another node);   V f   is the set of leaves (nodes without children);   μ   :   V ×   R ≥ 0   →   R +   is the valuation function, designating each node a real value at any given time   t   ≥   0 ;   L   is a finite set of labels;   l   :   E → L   is the labeling function that maps each edge to a label. Let   A ( t )   ∈   R 2 D   × 2 D  represent a time-varying   2 D   ×   2 D  matrix, where   D   ∈   N   is the   depth   of the matrix and   t   ∈   R ≥ 0  is time. We construct a QTS from   A ( t )   as follows. We let the  root   node   v ι   represent   A ( t ) . Next, we partition the matrix into four   2 D − 1   ×   2 D − 1   sub-matrices, where each sub-matrix is represented by a child of   v ι . We label each edge with a  directional   label from the set   L   =   { N W, N E, SE, SW   } , where   N W   represents north west,   SE   represents south east, etc (see Fig. 2). Next, we execute the same procedure for each child until the leaf nodes are obtained, i.e. each leaf is a single element matrix. Note that   | V f   |   = 2 2 D   and   |V|   =  D ∑  i =0  2 2 D . For each leaf node   v f   ∈ V f   , we let   μ ( v f   , t )   to be the value of the corresponding element in   A ( t ) . The valuation function for other nodes is recursively defined as the sum of the valuations of its children, i.e.  μ ( v, t ) =   ∑  ( v,v c   ) ∈E  μ ( v c , t )   ∀ v   ∈ V \   V f   .   (1) An example of a QTS construction is given in Fig. 2. Given a subset of labels   B ⊆ L , a   labeled path   of a QTS is defined as a function that maps a vertex to a set of infinite sequences of nodes:  λ B   ( v 0 ) :=   { ( v 0 , v 1 , v 2   · · ·   ,   v f   ) ∣ ∣ ( v i , v i +1 )   ∈ E , v f   ∈   V f   , l ( v i , v i +1 )   ∈ B , i   ∈   N ≥ 0 } ,  (2) where   v f   denotes infinite repetitions of leaf node   v f   . The  i -th element of a labeled path   π   ∈   λ B ( v )   is denoted by   π i . For example, in Fig. 2,   ( v 0 , v 1 ,   v 5 )   and   ( v 0 , v 1 ,   v 6 )   are both members of   λ { N W,N E } ( v 0 ) . A   QTS signal   starting at time   t   is defined as   Q t   =  { Q ( τ   ) | τ   ≥   t } .  B. Spatial Temporal Logic (SpaTeL)  SpaTeL formulas are defined by nesting tree spatial su- perposition logic (TSSL) specifications [13] inside tempo- ral operators of signal temporal logic (STL) [20]. Formal definitions of SpaTeL syntax and semantics can be found in   [14].   Informally, SpaTeL formulas are STL   formulas in which linear predicates over signals are replaced with spatial formulas over quad transition systems and allow for describing how spatial patterns change over time. A TSSL formula is recursively formed by linear predicates over the valuation function (1), spatial operators, and boolean operators ( ∧ , ∨ , ¬ ). For example,   φ   :=   μ   ∼   c , where   ∼∈ {≥  , < } , is a very simple TSSL formula consisting of a single predicate that indicates that the function   μ   of (1) at the initial node   v ι   must have a value of larger (smaller) than threshold  c . If the predicate is true, we write   Q v ι   | =   φ   (read as QTS   Q  satisfies formula   φ   at node   v ι ). Spatial operators are used to define specifications at nodes located in lower tree levels. For instance,   ∃ B   ©   φ   ( B ⊆ { N W, N E, SW, SE } ) is the spatial “there exists next” operator which means that the spatial formula   φ   has to be satisfied for at least one of the children of the initial node with directional label   l ( v ι , v ′ )   ∈ B . Fur- thermore,   ∀ B   © φ , read as “for all next”, indicates that   φ   must be satisfied by all such children. Specifications at deeper tree levels can be expressed similarly by nesting several spatial next operators. In addition to spatial next, TSSL is equipped with spatial until operators ( ∃ B φ 1 U κ φ 2 ,   ∀ B   φ 1 U κ φ 2 ). Formal definitions for all these operators are presented in [13]. A SpaTeL formula is recursively formed by TSSL for- mulae, temporal operators, and boolean operators. Three common temporal operators are   eventually ( F I   ),   always ( G I   ), and until ( U I   ), where   I   = [ t 1 , t 2 )   is a time interval. For instance, A QTS signal   Q t   | =   F I   φ   ( read as   Q t   satisfies  F I   φ ) if there   ∃ τ   ∈   [ t   +   t 1 , t   +   t 2 )   such that   Q τ   | =   φ   and  Q t   | =   G I   φ   if   ∀ τ   ∈   [ t + t 1 , t + t 2 )   Q τ   | =   φ . A SpaTel formula is satisfied by a QTS signal   Q t   if and only if it is satisfied by the QTS at time t ( Q ( t ) ).
16  v ι  0  0 0 0 0  NW   NE   SW   SE  2  1 1 0 0  NW   NE   SW   SE  7  3 2 1 1  NW   NE   SW   SE  7  0 0 4 3  NW   NE   SW   SE  NW  NE   SW  SE  A   =       3   4   1   1  0   0   2   3  0   0   0   0  1   1   0   0       NW   NE SW   SE  v 1   v 2   v 3   v 4  v 5   v 6   v 7   v 8   v 9   v 10   v 11   v 12   v 13   v 14   v 15   v 16   v 17   v 18   v 19   v 20  Fig. 2.   The QTS corresponding to the matrix   A .  Example 1:   Consider a   4   ×   4   matrix with the requirement that every other entry is zero (thus forming a checkerboard pattern). There are two different realizations for this pattern, that can be specified by the following TSSL formulas:  φ c 1   =   ∀ L   ©   ( ∀ { NW,SE }   ©   ( μ   = 0))  φ c 2   =   ∀ L   ©   ( ∀ { NE,SW   }   ©   ( μ   = 0)) ,  where   μ   is the valuation function in the definition of QTS. Now consider a spatio-temporal requirement that the checkerboard pattern periodically switches tiles. The follow- ing SpaTeL formula specifies this requirement:  Φ c   =   G [0 ,t 1   ) ( F [0 ,t 2 ) φ c 1   ∧   F [0 ,t 2   ) φ c 2 ) .   (3) SpaTeL is equipped with quantitative semantics. Quanti- tative valuation (robustness)   ρ (Φ , Q t )   of a SpaTeL formula  Φ   with respect to QTS signal   Q t   is calculated recursively:  ρ ( φ, Q t ) =   ρ s ( φ, v ι ) , ρ ( ¬ Φ , Q t ) =   − ρ (Φ , Q t ) , ρ (Φ 1   ∧   Φ 2 , Q t ) = min( ρ (Φ 1 , Q t ) , ρ (Φ 2 , Q t )) , ρ (Φ 1   ∨   Φ 2 , Q t ) = max( ρ (Φ 1 , Q t ) , ρ (Φ 2 , Q t )) , ρ ( F I 1   ,I 2 ) Φ , Q t ) = sup t ′   ∈ [ t + I 1 ,t + I 2   )   ρ (Φ , Q t ′   ) , ρ ( G I 1   ,I 2 ) Φ , Q t ) = inf t ′   ∈ [ t + I 1 ,t + I 2   )   ρ (Φ , Q t ′   ) , ρ (Φ 1 U [ I 1 ,I 2   ) Φ 2 , Q t ) = sup t ′ ∈ [ t + I 1   ,t + I 2 ) ( min( ρ (Φ 2 , Q t ′   ) ,   inf t ′′   ∈ [ t + I 1 ,t ′ )   ρ (Φ 1 , Q t ′′   ))) ,  (4) where   ρ s ( φ, v )   is the robustness of TSSL formula   φ   with respect to node   v   ∈ V :  ρ s ( ⊤ , v ) = 1 , ρ s ( μ   ∼   c, v ) = ( μ   −   c )   if   ( ∼   is   ≥ ) ,   ( c   −   μ )   if   ( ∼   is   ≤ ) , ρ s ( ¬ φ, v ) =   − ρ s ( φ, v ) , ρ s ( φ 1   ∧   φ 2 , v ) = min( ρ s ( φ 1 , v ) , ρ s ( φ 2 , v )) , ρ s ( φ 1   ∨   φ 2 , v ) = max( ρ s ( φ 1 , v ) , ρ s ( φ 2 , v )) , ρ s ( ∃ B   ©   φ, v ) = max π ∈ λ B   ( v )   ρ s ( φ, π 1 ) , ρ s ( ∀ B   ©   φ, v ) = min π ∈ λ B   ( v )   ρ s ( φ, π 1 ) , ρ s ( ∃ B   φ 1 U k φ 2 , v ) = sup π ∈ λ B   ( v ) ,i ∈ (0 ,k ] ( min( ρ s ( φ 2 , π i ) ,   inf j ∈ [0 ,i )   ρ s ( φ 1 , π j   ))) , ρ s ( ∀ B   φ 1 U k φ 2 , v ) = inf π ∈ λ B   ( v ) ,i ∈ (0 ,k ] ( min( ρ s ( φ 2 , π i ) ,   inf j ∈ [0 ,i )   ρ s ( φ 1 , π j   ))) .  (5) Positive and negative robustness indicate satisfaction and violation, respectively.  ρ (Φ , Q t )   >   0   ⇒   Q t   | = Φ , ρ (Φ , Q t )   <   0   ⇒   Q t   6   | = Φ .   (6) The absolute robustness value can be viewed as a measure of ”distance to satisfaction”. In other words, a higher absolute value for robustness indicates stronger satisfaction (violation) of a specification. We use the definition of robustness ((4), (5)) in subsequent sections to translate SpaTeL specifications into mixed integer constraints.  Example 2 (Example 1 continued):   Consider a stationary QTS   Q ( t ) =   Q , where   Q   is the QTS depicted in Fig. 2. By computing quantitative semantics from (5), it is straightfor- ward to verify that specification   Φ c   in (3) is violated by   Q 0  with a robustness of   − 4 . The horizon of a SpaTeL formula is defined similar to STL [21]. Intuitively, the horizon   T   of a SpaTeL formula  Φ   is the maximum time for which some specification in   Φ  must be checked against   Q t . For instance, the time horizon of   G [0 , 20) F [0 , 5) ∀ L   ©   ( μ   ≥   1)   is   T   = 25 . III. P ROBLEM   F ORMULATION AND   A PPROACH  Consider   N   homogenous planar robots with negligible sizes in a two-dimensional space. The position of robot   r   at time   t   is denoted by   x r   ( t )   ∈ X   ,   r   = 1 ,   · · ·   , N   , where   X ⊂  R 2   is the workspace of the robots, which is assumed to be the following square:   X   := [ −   a  2   ,   a  2   ] × [ −   a  2   ,   a  2   ] ,   where   a   is the length of the square. Note that any rectangular workspace can be normalized to meet this assumption. We denote the state of the swarm by   x ( t ) =   ( x 1 ( t ) T   , x 2 ( t ) T   ,   · · ·   , x N   ( t ) T   ) T   . The   kinematics of   each   individual robot is   assumed   as follows:   ̇ x r   ( t ) =   u r ( t ) ,   r   = 1 ,   · · ·   , N,   (7) where   u r ( t )   ∈ U   is the control applied to robot   r   at time  t   and   U   =   { u r  ∣ ∣   ‖ u ‖ 2   ≤   u m  } , where   u m   is the maximum speed that a robot can attain.  X   is partitioned into   2 D   × 2 D   number of equal-sized cells, where   D   is the depth of the grid. A user expresses desirable patterns by defining shapes that are formed by unions of cells in the workspace, defining thresholds for the number
• • • • • • • • • • •   • • • • • • • • • •        1   1   0   4 2   0   1   0 1   3   2   0 0   1   2   3       Fig. 3.   (Left) A swarm of   21   robots in a square region. The square is gridded into 16 cells. (Right) The matrix representing the number of the robots in each cell.  of agents populating each shape, and expressing temporal requirements for those pre-identified patterns. The objective in this paper is to synthesize a control policy for (7) such that the spatio-temporal requirements expressed by the user are met. We will provide a formal formulation for this problem later in this section. We construct the matrix   N   ( t )   ∈   N 2 D   × 2 D  , where the value of each element is the number of robots in the corresponding cell, as illustrated in an example in Fig. 3. We construct the time varying QTS   Q ( t )   from   N   ( t )   using the procedure outlined   in   Sec.   II-A.   Note   that   the   shapes   defined   by unions of cells can be easily expressed using the spatial next operator in tree spatial superposition logic (see Example 3). Consequently, a SpaTeL specification   Φ   can be automatically generated from the input specification.  Remark 1:   A supervised learning algorithm is proposed in [13] for automatically learning TSSL formulae that are satisfied by a set of positive spatial configurations (images) and violated by a set of negative images. This method can be used to learn TSSL (and SpaTeL) formulas describing more complex high level patterns (circular clusters, ellipsoids, etc). Two sets of images, one illustrating the desirable pattern and   the   other   lacking the   pattern, should   be   artificially created. The learning algorithm generates a TSSL formula by using the generated images as a training set. Although this is a very effective method to find SpaTeL descriptors for arbitrary patterns, the resulting formulas are often too long and complex. Therefore, The mixed integer linear program- ming problems that result from the framework presented in subsequent sections become unsolvable by existing solvers. As a result, the input spatio-temporal requirements in this paper are limited to unions of squares. As explained in Sec. II-B, such requirements can be intuitively translated into TSSL/SpaTeL specifications and the resulting formulae are small and manageable.  Example 3:   The specification that was introduced in Sec. I (Fig. 1) is formalized by the following SpaTeL formula:  Φ d   =   G [0 , 40) ( ¬ φ 1 )   ∧   ( F [0 , 30) G [0 , 3) φ 2 )  ∧ ( F [0 , 30) G [0 , 3) φ 3 )   ∧   ( F [30 , 40) φ 4 ) ,   (8) where   φ i   are TSSL formulas describing different patterns illustrated in Fig. 1:   φ 1   represents the red danger zone in Fig. 1,   φ 2   specifies formation of a checkerboard pattern in the north west quadrant,   φ 3   specifies gathering inside one of the grey cells in the south west quadrant, and   φ 4   represents populating one of the L-shaped cyan regions. These formulas are automatically generated by representing each cell in the gridded workspace using appropriate spatial next operators of TSSL.  φ 1   =   ∀ SE   © ∀ N W   ©   ( μ   ≤   0) ∧ ∀ SW   © ∀ N E   © ∀ { N W,N E }   ©   ( μ   ≤   0) ∧ ∀ SW   © ∀ N W   © ∀ N E   ©   ( μ   ≤   0) , φ 2   =   ∀ N E   ©   ( ∀ L   © ∀ { N W,SE }   ©   ( μ   ≥   γ 1 )) , φ 3   =   ∀ SW   ©   ( ∀ SW   © ∃ L   ©   ( μ   ≥   γ 2 )) , φ 4   =   ∀ N W   ©   ( φ 5   ∨   φ 6 ) , φ 5   =   ( ∀ N E   © ∀ { N W,N E,SE }   ©   ( μ   ≥   γ 3 )) ∧  ( ∀ SE   © ∀ N E   ©   ( μ   ≥   γ 4 )) , φ 6   =   ( ∀ SW   © ∀ { N W,SW,SE }   ©   ( μ   ≥   γ 5 )) ∧  ( ∀ N W   © ∀ SW   ©   ( μ   ≥   γ 6 )) ,  (9) where   L   =   { N W, N E, SW, SE } , and   μ   is the the number of robots residing in a subregion of the workspace identified by spatial operators and   γ 1 − 6   are thresholds for the minimum number of robots that are required to populate each pattern. We wish to find a control strategy that steers the swarm such that   Φ   is satisfied. Such a policy is not usually unique. Therefore, we choose a policy that optimizes a cost function. For instance, we can minimize the total number of robot displacements (one displacement is defined as moving one robot from its current location to a neighboring cell). In addition, a natural candidate for optimization is maximizing the SpaTeL robustness. The problem that we consider in this paper is formulated as follows:  Problem 1:   Given a swarm of   N   agents with initial po- sitions at   x (0)   and a SpaTeL formula   Φ   that describes time varying spatial requirements of the user, find an   optimal and correct   control strategy such that:  u r   ( t )  r =1 , ···   ,N,t ∈ [0 ,T   ]  =   argmin   − α ρ (Φ , Q 0 ) +   J f   ( x ( T   )) +   ∫   T  0   J r   ( x ( t ) , u ( t )) dt, s.t.   Q 0   | = Φ ,  (10) where   ρ   is the SpaTeL robustness,   Q 0   is the QTS signal starting at time   0 ,   J f   :   R 2 N   →   R ,   J r   :   R 2 N   ×   U N   →   R , are the endpoint cost and the running cost (Lagrangian), respectively. The end time   T   is the time horizon of the SpaTeL formula and   α   is a positive constant designating a weight for SpaTeL robustness. Our approach to problem 1 can be summarized as follows. First, we find   N   ( t ) ,   0   ≤   t   ≤   T   such that   Q 0   | = Φ . It is known that this problem is undecidable in continuous time [22]. Therefore, we (approximately) solve the problem in discrete time assuming that at each time step, each robot can be displaced by one cell to its right, left, up or down. Therefore, we choose a sampling time such that:   ∆ t   ≥  a  2 D − 1   u m   .   We assume that the time intervals of the temporal operators of   Φ   are multiples of   ∆ t . This assumption can be matched by increasing   D   such that the time intervals can be reasonably approximated by multiples of   ∆ t . We also denote the last discrete time by   K   :=   T  ∆ t   . The matrix   N   at time  t   =   k ∆ t   is denoted by   N   [ k ] . We construct a discrete time model for the evolution of   N   [ k ] . Next, we find the required values at each time such that the SpaTeL specification is
4 ,   1  3 ,   1  2 ,   1  1 ,   1   1 ,   2  2 ,   2  3 ,   2  4 ,   2   4 ,   3  3 ,   3  2 ,   3  1 ,   3   1 ,   4  2 ,   4  3 ,   4  4 ,   4  Fig. 4.   A flow network with   16   cells. The robots are only able to move from one cell to a neighboring cell in one time step.  satisfied using a MILP-based approach that is explained in Sec. IV. Finally, we find continuous time controls for each individual robot such that the number of each cell at time  t   =   k ∆ t   matches its corresponding value in   N   [ k ] . IV. S OLUTION  A. Swarm Flow in Discrete Time  In this section, we develop a discrete time model that characterizes the evolution of   N   [ k ] . At each time step, each robot is only able to remain at its current cell or move to an adjacent cell (The cell to its right, left, top, or down). All the robots move synchronously during one time step. The flow of the robots between the cells can be thought as a network as depicted in Fig. 4. The index of each cell is represented by  [ i, j ] , where   i   is the row and   j   is the column of the element in the matrix   N   [ k ] . The set of cells that are adjacent to   [ i, j ]   is denoted by   Ω([ i, j ]) . We denote the number of robots in the cell   [ i, j ]   at time step   k   by   N [ i,j ] [ k ] . The number of robots that move from cell   [ i, j ]   to an adjacent cell   [ i ′ , j ′ ]   ∈   Ω([ i, j ])  during time   [ k, k   + 1]∆ t   is denoted by   f   [ i ′ ,j ′   ] [ i,j ]   [ k ] , which is a non-negative integer. The total number of robots that move out from cell   [ i, j ]   at time step   k   is:  f   out  [ i,j ] [ k ] :=   ∑  [ i ′   ,j ′   ] ∈ Ω([ i,j ])  f   [ i ′ ,j ′   ] [ i,j ]   [ k ] .   (11) We add the following constraint:  N [ i,j ] [ k ]   ≥   f   out  [ i,j ] [ k ] ,   (12) which indicates that the number of robots moving out from a cell can not be more than the number of robots in the cell. The number of robots that enter cell   [ i, j ]   at   k   is:  f   in  [ i,j ] [ k ] :=   ∑  [ i ′   ,j ′   ] ∈ Ω([ i,j ])  f   [ i,j ] [ i ′ ,j ′   ] [ k ] .   (13) The discrete time evolution of   N   [ k ]   is:  N [ i,j ] [ k   + 1] =   N [ i,j ] [ k ]   −   f   out i,j   [ k ] +   f   in i,j   [ k ] ,   (14) which is a function of decisions made on the values of  f   [ i ′ ,j ′   ] [ i,j ]   [ k ] . In a compact form, we define the decision variable  f   [ k ]   as the set:  f   [ k ] =  {  f   [ i ′ ,j ′   ] [ i,j ]  ∣ ∣ ∀ [ i ′ , j ′ ]   ∈   Ω([ i, j ]) ,   ∀ [ i, j ]  }  ,   (15) and the discrete time evolution of   N   [ k ]   is written as:  N   [ k   + 1] =   F   ( N   [ k ] , f   [ k ]) .   (16)  B. Mixed-Integer Formulation of SpaTeL Specifications  In this section, we explain how to recursively transform a SpaTeL formula into a set of mixed-integer constraints. Our method is inspired by the binary mixed-integer encoding of STL formulas presented in [22]. For a predicate of a SpaTeL formula   σ   = ( μ   ≥   c ) , a set of binary variables   z σ [ v, k ]   ∈ { 0 ,   1 } ,   v   ∈ V ,   0   ≤   k   ≤  K , is associated such that values 1 and 0 indicate   T rue  and   F alse , respectively. The corresponding mixed integer constraints are:  {   μ [ v, k ]   −   M z σ [ v, k ]   ≤   c, μ [ v, k ] +   M   (1   −   z σ [ v, k ])   ≥   c,   (17) where   M   is a sufficiently large positive number. Mixed integer constraints for all predicates in the form of   σ ′   = ( μ   ≤  c )   are defined similarly. For encoding a SpaTeL formula, The following rules are used to map boolean, temporal, and spatial operators into mixed integer constraints. These rules are derived from the definition of SpaTeL robustness (4),(5).  •   Negation :  Ψ =   ¬ Φ   →   z Ψ [ v, k ] = 1   −   z Φ [ v, k ];  •   Conjunction :  Ψ =  m ∧  i =1  Φ i   →      z Ψ [ v, k ]   ≤   z Φ i   [ v, k ] , i   = 1 ,   · · ·   , m, z Ψ [ v, k ]   ≥   1   −   m   +   m ∑  i =1  z Φ i   [ v, k ];  •   Disjunction :  Ψ =  m ∨  i =1  Φ i   →      z Ψ [ v, k ]   ≥   z Φ i   [ v, k ] , z Ψ [ v, k ]   ≤   m ∑  i =1  z Φ i   [ v, k ];  •   There exists spatial next :  Ψ =   ∃ B   ©   Φ   →   z Ψ [ v, k ] =   ∨  π ∈ λ B   ( v )  z Φ [ π 1 , k ];  •   For all spatial next :  Ψ =   ∀ B   ©   Φ   →   z Ψ [ v, k ] =   ∧  π ∈ λ B   ( v )  z Φ [ π 1 , k ];  •   There exists spatial until :  Ψ =   ∃ B Φ 1 U κ Φ 2   →  z Ψ [ v, k ] =   ∨  π ∈ λ B   ( v ) ,i ∈ (0 ,κ ]  ( z Φ 2   [ π i , k ]   ∧   ∧  j ∈ [0 ,i )  z Φ 1   [ π j   , k ]);  •   For all spatial until :  Ψ =   ∀ B Φ 1 U κ Φ 2   →  z Ψ [ v, k ] =   ∧  π ∈ λ B   ( v ) ,i ∈ (0 ,κ ]  ( z Φ 2   [ π i , k ]   ∧   ∧  j ∈ [0 ,i )  z Φ 1   [ π j   , k ]);  •   Temporal eventually :  Ψ =   F [ k 1   ∆ t,k 2 ∆ t ) Φ   →   z Ψ [ v, k ] =   ∨  k ′   = k 1   , ···   ,k 2  z Φ [ v, k ′ ];
•   Temporal always :  Ψ =   G [ k 1   ∆ t,k 2 ∆ t ) Φ   →   z Ψ [ v, k ] =   ∧  k ′   = k 1   , ···   ,k 2  z Φ [ v, k ′ ];  •   Temporal until :  Ψ = Φ 1 U [ k 1 ∆ t,k 2 ∆ t ) Φ 2   →  z Ψ [ v, k ] =   ∨  k ′   = k 1 , ···   ,k 2  ( z Φ 2   [ v, k ′ ] ∧   ∧  k ′′   = k 1 , ···   ,k ′  z Φ 1   [ v, k ′′ ]) .  Note that   z Ψ [ v, k ]   ∈   [0 ,   1]   is not required to be declared an integer since it is automatically enforced to take binary values. Finally, the problem of satisfying a general SpaTeL formula,   Q 0   | = Φ , reduces to the following constraint:  z Φ ( v ι ,   0) = 1 ,   (18) where   v ι   is the root node of quad transition system.  C. Robustness-Based Encoding  In this section, we briefly explain how to incorporate Spa- TeL robustness into the mixed-integer encoding. The method is much in spirit of the method in [23], where the authors characterize the changes in the satisfaction of the specifi- cation with respect to the changes in the predicates. For a predicate in the form of   σ   = ( μ   ∼   c ) , it is straightforward to see from (4) and (5) that   ∂ρ (Φ ,Q 0   )  ∂c   ∈ { 0 ,   1 }   (non-decreasing) or   ∂ρ (Φ ,Q 0   )  ∂c   ∈ { 0 ,   − 1 }   (non-increasing), depending on the operators preceding the predicate. Therefore, by increas- ing (decreasing) the value of   c   for a non-increasing (non- decreasing) predicate, a constraint is   tightened . Therefore, we alter the values of   c   in the predicates as follows:  {  c   ←   c   +   ̺   ∂ρ (Φ ,Q 0   )  ∂c   ∈ { 0 ,   − 1 } , c   ←   c   −   ̺   ∂ρ (Φ ,Q 0   )  ∂c   ∈ { 0 ,   1 } .   (19) Next, we add the constraint   ̺   ≥   0   to ensure satisfaction of   Φ . It is easy to show that the maximum   ̺   that renders   Q 0   | = Φ  is equal to   ρ (Φ , Q 0 ) .  D. High Level Planning  In the previous sections, we formulated the dynamics and SpaTeL objectives as mixed-integer constraints. We formu- late the discrete-time version of Problem 1 as the following MILP:  f   [ k ]  k =0 , 1 , ···   ,K  =   argmin   − α ̺   +   J f   ( N   [ K ]) +   ∑ K  0   J r   ( N   [ k ] , f   [ k ]) , s.t.   N   [ k   + 1] =   F   ( N   [ k ] , f   [ k ]) , z Φ ( v ι ,   0) = 1 , ̺  ≥   0 ,  (20) where   J f   and   J r   are the discrete time versions of the end- point and running cost, respectively. Note that we assume the costs are linear functions. In this paper, we are particularly interested in the following cost:  J r   ( N   [ k ] , f   [ k ]) =   ∑   f   [ k ] ,   (21) which corresponds to the total number of robot displacements (energy). Note that all the values of   f   are non-negative. In case the MILP above is infeasible, no control strategy is able to satisfy the SpaTeL formula. In this case, we relax the last constraint   ̺   ≥   0 , and choose a very large value for   α  (or remove the other costs). Therefore, the resulting solution solely maximizes the SpaTeL robustness, which is a negative value. In other words, the SpaTeL violation is minimized.  E. Low Level Control Policy  The decision variables   f   [ i ′ ,j ′   ] [ i,j ]   [ k ]   are obtained from the solution to (20). The only remaining problem is to choose the set of individual robots that must be moved from cell   [ i, j ]   to adjacent cells. For this purpose, we choose   f   [ i ′   ,j ′   ] [ i,j ]   [ k ]   number of robots that are closest to the edge between   [ i, j ]   and   [ i ′ , j ′ ]  and move them with a constant velocity on a straight line. In other words, if   R [ i,j ] [ k ]   is the the set of robot indices that are located inside cell   [ i, j ]   at time step   k   and   R [ i ′   ,j ′   ] [ i,j ]   [ k ]   is the set of robot indices that are supposed to be moved from   [ i, j ]  to   [ i ′ , j ′ ]   ∈   Ω([ i, j ])   at time step   k , the control law would be:  u r ( t ) =   a  2 D   ∆ t   .  (   j ′   −   j i   −   i ′  )  ,   r   ∈ R [ i ′   ,j ′   ] [ i,j ]   [ k ] , k ∆ t   ≤   t <   ( k   + 1)∆ t.  (22) Algorithm 1 presents the procedure to determine   R [ i ′   ,j ′   ] [ i,j ]   [ k ] .  Data :   [ i, j ]   (cell index),   R [ i,j ] [ k ]   (robots inside that cell),   f   [ k ]   (flow variables)  Result :   R [ i ′   ,j ′   ] [ i,j ]   [ k ]   ( set of robots that are moved from  [ i, j ]   to a neighbor)  while   ∑  [ i ′   ,j ′   ] ∈ Ω([ i,j ])   f   [ i ′ ,j ′   ] [ i,j ]   [ k ]   >   0   do  find robot   r   ∈ R [ i,j ] [ k ]   that has the minimum distance to the edges of   [ i, j ] ; add   r   to   R [ i ′   ,j ′   ] [ i,j ]   [ k ]   where   [ i ′ , j ′ ]   is the neighbor of  [ i, j ]   that   r   was closest to; remove   r   from   R [ i,j ] [ k ] ;  f   [ i ′ ,j ′   ] [ i,j ]   [ k ]   ←   f   [ i ′   ,j ′   ] [ i,j ]   [ k ]   −   1 ;  end Algorithm 1:   How to assign robots to move from a cell to its neighbors  Remark 2:   As   mentioned   earlier,   we   do   not   consider physical sizes for robots in this paper. In practice, careful strategies are required for collision avoidance among the robots. This issue   will   be   further investigated in   future work. As a preliminary solution, we propose the following approach. Let   n cap   be   the   maximum number of robots that can populate one cell without physically overlapping one another. We basically add the specification that for all cells the number of robots should not exceed   n cap , which guarantees there is always enough empty space for swarm movements. Localized policies such as the methods in [24] can be used for guaranteeing collision avoidance.
a)   t   = 0   b)   t   = 3   c)   t   = 7   d)   t   = 11   e)   t   = 16   f)   t   = 18  g)   t   = 21   h)   t   = 24   i)   t   = 27   j)   t   = 30   k)   t   = 32   l)   t   = 40  Fig. 5.   Case study: Snapshots of the optimal swarm movement satisfying SpaTel formula (8) starting from the initial condition shown in figure a). First the robots are gathered in one cell in the grey region to satisfy   φ 3 , then robots move toward forming the checkerboard pattern, satisfying   φ 2 . Finally robots move to the populate the upper L-shaped pattern, satisfying   φ 5 .  a)   t   = 0   b)   t   = 1   c)   t   = 3   d)   t   = 6   e)   t   = 10   f)   t   = 13  g)   t   = 15   h)   t   = 18   i)   t   = 20   j)   t   = 22   k)   t   = 31   l)   t   = 40  Fig. 6.   Case study: Snapshots of the optimal swarm movement satisfying SpaTel formula (8) starting from the initial condition shown in figure a). First the robots are forming the checkerboard pattern   φ 2 , then they gather in one grey cell to satisfy   φ 3 , and finally robots move to the populate the lower L-shaped pattern, satisfying   φ 6 .  F. Complexity  The worst case complexity of the framework outlined in previous sections depends on the complexity of the MILP formulation in Sec. IV-B. The complexity of a MILP problem grows exponentially, in the worst case, with respect to the number of integer variables and polynomially with respect to the continuous variables. It is worth to note that for large swarms, the flow variables can be approximated as continuous numbers and   be   rounded off   after   obtaining the MILP solution. This approximation significantly reduces the computational complexity without significantly altering the optimal solution and makes the complexity independent from the total number of robots. The number of integer variables in (20) is   KP   |V f   | , where   K   is the total number of time steps,   P   is the number of linear predicates in   Φ , and   |V f   |   is the number of cells in the gridded workspace. This suggests that the framework might not be scalable for grids with high resolutions or extremely complicated and long spatio-temporal requirements. However, as illustrated in the examples in Sec. V, quite complicated patterns are achievable in practice with relatively low computation time. V. C ASE   S TUDY  640   robots in a workspace partitioned into a   8   ×   8   grid. The SpaTeL formula of (8) corresponding to the specification of Fig. 1 is the   target, where we   set   γ 1   =   80 , γ 2   = 640 , γ 3 − 6   = 160 . The cost function that is minimized is the total   number of robot displacements given by (21). We demonstrate the results for two different initial con- ditions.   A   movie   illustrating   both   cases   is   available   on  https://youtu.be/x-uI8N9iN3I .  A. Case 1  We set the initial configuration of robots to be in the uni- formly distributed in the   SW   quadrant of the   SE   quadrant (see Fig. 5 a). We formulate (20) as a MILP, which we solve
using Gurobi   1 . The MILP is solved in 54 seconds on a 3GHz Dual core Macbook Pro. Next, we move the robots according to the plan obtained from the solution of the MILP. Fig. 5 shows snapshots of the swarm movement during its completion of the mission described by (8). It is seen that the swarm first satisfies   φ 3 , then   φ 2   and then   φ 5 .  B. Case 2  Now we set the initial condition to be uniformly dis- tributed in the   N W   quadrant of the   N E   quadrant (see Fig. 6 a).The MILP is solved in 43 seconds. The snapshots of the swarm movement are shown in Fig. 6. This time, the optimal plan is to satisfy   φ 2   first, then   φ 3   and finally   φ 6 . VI. C ONCLUSION AND   F UTURE   W ORK  We presented a fully automated framework to synthesize controls that steer a fully actuated robotic swarm while sat- isfying high-level spatio-temporal requirements. Such spec- ifications are naturally expressed as spatial temporal logic formulae. A computationally efficient framework was intro- duced to determine the high level plan from which a low level control strategy executes the swarm movements. Directions for future research include extending the frame- work to under actuated swarms and developing distributed control strategies for coordination of   movements among robots. We also plan to incorporate machine learning meth- ods from [13] in order to synthesize control policies for more complex spatial patterns which are automatically learned from training data. Furthermore, we plan to create a graphical user interface in which a potential user can define required patterns for execution. The user would draw the patterns that they want to emerge and specify time requirements. The interface will use machine learning techniques to generate SpaTeL formulas for those patterns. These formulas will then be used by algorithms presented in this paper to synthesize control policies for the swarm. R EFERENCES [1]   M. Mesbahi and M. Egerstedt,   Graph theoretic methods in multiagent networks .   Princeton University Press, 2010. [2]   N. Michael, M. M. Zavlanos, V. Kumar, and G. J. Pappas, “Distributed multi-robot task assignment and formation control,” in   Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on . IEEE, 2008, pp. 128–133. [3]   J. Cortes, S. Martinez, T. Karatas, and F. Bullo, “Coverage control for mobile sensing networks,” in   Robotics and Automation, 2002. Proceedings. ICRA’02. IEEE International Conference on , vol. 2. IEEE, 2002, pp. 1327–1332. [4]   F. Bullo, J. Cort ́ es, and S. Martinez,   Distributed control of robotic net- works: a mathematical approach to motion coordination algorithms . Princeton University Press, 2009. [5]   G. Kantor, S. Singh, R. Peterson, D. Rus, A. Das, V. Kumar, G. Pereira, and J. Spletzer, “Distributed search and rescue with robot and sensor teams,” in   Field and Service Robotics .   Springer, 2006, pp. 529–538. [6]   S.   Thrun   and   Y.   Liu,   “Multi-robot   SLAM   with   sparse   extended information filers,” in   Robotics Research. The Eleventh International Symposium .   Springer, 2005, pp. 254–266. [7]   Y. Chen and Z. Wang, “Formation control:   a review   and a new consideration,” in   Intelligent Robots and Systems, 2005.(IROS 2005). 2005 IEEE/RSJ International Conference on .   IEEE, 2005, pp. 3181– 3186.  1 www.gurobi.com [8]   L. C. Pimenta, N. Michael,   R. C. Mesquita,   G. A. Pereira,   and V. Kumar, “Control   of swarms based on hydrodynamic   models,” in   Robotics and Automation, 2008. ICRA 2008. IEEE International Conference on .   IEEE, 2008, pp. 1948–1953. [9]   M. Egerstedt and X. Hu, “Formation constrained multi-agent control,”  IEEE TRANSACTIONS ON ROBOTICS AND AUTOMATION , vol. 17, no. 6, p. 947, 2001. [10]   S. G. Lee and M. Egerstedt, “Multi-robot control using time-varying density functions,”   arXiv preprint arXiv:1404.0338 , 2014. [11]   P. Yang, R. Freeman, K. M. Lynch, and Others, “Multi-agent coor- dination by decentralized estimation and control,”   Automatic Control, IEEE Transactions on , vol. 53, no. 11, pp. 2480–2496, 2008. [12]   M. Kloetzer and C. Belta, “Temporal logic planning and control of robotic swarms by hierarchical abstractions,”   Robotics, IEEE Trans- actions on , vol. 23, no. 2, pp. 320–330, 2007. [13]   E. A. Gol, E. Bartocci, and C. Belta, “A formal methods approach to pattern synthesis in reaction diffusion systems,” in   Decision and Control (CDC), 2014 IEEE 53rd Annual Conference on .   IEEE, 2014, pp. 108–113. [14]   I. Haghighi, A. Jones, Z. Kong, E. Bartocci, R. Gros, and C. Belta, “SpaTeL: a novel spatial-temporal logic and its applications to net- worked systems,” in   Proceedings of the 18th International Conference on Hybrid Systems: Computation and Control .   ACM, 2015, pp. 189– 198. [15]   A. F. Winfield, J. Sa, M.-C. Fern ́ andez-Gago, C. Dixon, and M. Fisher, “On formal specification of emergent behaviours in swarm robotic systems,”   International journal of advanced robotic systems , vol. 2, no. 4, pp. 363–370, 2005. [16]   Y. Chen, X. C. Ding, A. Stefanescu, and C. Belta, “Formal approach to   the   deployment   of   distributed   robotic   teams,”   Robotics,   IEEE Transactions on , vol. 28, no. 1, pp. 158–171, 2012. [17]   J. Tumova and D. V. Dimarogonas, “A receding horizon approach to multi-agent planning from local LTL specifications,” in   American Control Conference (ACC), 2014 .   IEEE, 2014, pp. 1775–1780. [18]   A. Ulusoy, S. L. Smith, X. C. Ding, and C. Belta, “Robust multi-robot optimal path planning with temporal logic constraints,” in   Robotics and Automation (ICRA), 2012 IEEE International Conference on .   IEEE, 2012, pp. 4693–4698. [19]   Y. Diaz-Mercado, A. Jones, C. Belta, and M. Egerstedt, “Correct-by- construction control synthesis for multi-robot mixing,” in   Decision and Control (CDC), 2015 IEEE 54th Annual Conference on .   IEEE, 2015, pp. 221–226. [20]   O. Maler and D. Nickovic, “Monitoring temporal properties of contin- uous signals,” in   Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems .   Springer, 2004, pp. 152–166. [21]   A. Dokhanchi, B. Hoxha, and G. Fainekos, “On-line monitoring for temporal logic robustness,” in   Runtime Verification .   Springer, 2014, pp. 1–20. [22]   V. Raman, A. Donz ́ e, M. Maasoumy, R. M. Murray, A. Sangiovanni- Vincentelli, and S. A. Seshia, “Model predictive control with signal temporal logic specifications,” in   Decision and Control (CDC), 2014 IEEE 53rd Annual Conference on .   IEEE, 2014, pp. 81–87. [23]   S. Sadraddini and C. Belta, “Robust Temporal Logic Model Predictive Control,”   53rd Annual Conference on Communication, Control, and Computing (Allerton) , 2015. [24]   J. Van Den Berg, S. J. Guy, M. Lin, and D. Manocha, “Reciprocal n-body collision avoidance,” in   Robotics research .   Springer, 2011, pp. 3–19.