MAPS 2 : Multi-robot autonomous motion planning under signal temporal logic specifications

Abstract

This article presents MAPS²: a distributed algorithm that allows multi-robot systems to deliver coupled tasks expressed as Signal Temporal Logic (STL) constraints. Classical control theoretical tools addressing STL constraints either adopt a limited fragment of the STL formula or require approximations of min/max operators. Meanwhile, works maximising robustness through optimisation-based methods often suffer from local minima, thus relaxing any completeness arguments due to the NP-hard nature of the problem. Endowed with probabilistic guarantees, MAPS² provides an autonomous algorithm that iteratively improves the robots’ trajectories. The algorithm selectively imposes spatial constraints by taking advantage of the temporal properties of the STL. The algorithm is distributed in the sense that each robot calculates its trajectory by communicating only with its immediate neighbours as defined via a communication graph. We illustrate the efficiency of MAPS² by conducting extensive simulation and experimental studies, verifying the generation of STL satisfying trajectories.

Introduction

Autonomous robots can solve significant problems when provided with a set of guidelines. These guidelines can be derived from either the physical constraints of the robot, such as joint limits, or imposed as human-specified requirements, such as pick-and-place objects. An efficient method of imposing such guidelines is by using logic-based tools, which enable reasoning about the desired behaviour of robots. These tools help us describe the behaviour of a robot at various levels of abstraction, such as interactions between its internal components to the overall high-level behaviour of the robot (Lamport, 1983). This strong expressivity helps us efficiently encode complex mission specifications into a logical formula. Recent research has focused on utilising these logic-based tools to express requirements on the behaviour of robots. Once these requirements are established, algorithms are developed to generate satisfying trajectories. Such is the focus of our work.

Examples of logic-based tools include formal languages, such as Linear Temporal Logic (LTL), Metric Interval Temporal Logic (MITL), and Signal Temporal Logic (STL). The main distinguishing feature between these logics is their ability to encode time. While LTL operates in discrete-time and discrete-space domain, MITL operates in the continuous-time domain but only enforces qualitative space constraints. On the other hand, STL allows for the expression of both qualitative and quantitative semantics of the system in both continuous-time and continuous-space domains (Maler and Nickovic, 2004). STL thus provides a natural and compact way to reason about a robot’s motion since it operates in a continuously evolving space-time environment. Additionally, STL is accompanied by a robustness metric which allows us to determine the extent of satisfaction compared to only absolute satisfaction (Donzé, 2013).

Another important property of autonomous robots is their ability to coordinate and work in teams. The use of multiple robots is often necessary in situations where a single robot is either insufficient, the task is high-energy demanding, or unable to physically perform certain tasks. However, multi-robot systems present their own set of challenges, such as communication overload, the need for a central authority for commands, and high computational demands. The challenge, therefore, is to derive solutions for multi-robot problems utilising logic-based tools, ensuring the achievement of specified high-level behaviour.

In this article, we propose MAPS² – Multi-Robot Autonomous Motion Planning under Signal Temporal Logic Specifications – to address the multi-robot motion-planning problem subject to coupled STL constraints. The algorithm encodes these constraints into an optimisation function and selectively activates them based on the temporal requirements of the STL formula. While doing so, each robot only communicates with its neighbours and iteratively searches for STL satisfying trajectories. The algorithm ensures distributed trajectory generation to satisfy STL formulas that consist of coupled constraints for multiple robots. The article’s contributions are summarised in the following attributes:

• The algorithm’s effectiveness lies in its ability to distribute STL planning for multiple robots and in providing a mechanism to decouple the STL formula among robots, thereby facilitating the distribution of tasks.

• As opposed to previous work, it covers the entire STL formula and is not limited to a smaller fragment. It reduces conservatism by eliminating the need for approximations of max/min operators and samples in continuous time to avoid abstractions.

• It incorporates a wide range of coupled constraints (both linear and nonlinear) into the distributed optimisation framework, enabling the handling of diverse tasks such as pick-and-place operations and time-varying activities like trajectory tracking.

• We present extensive simulation and hardware experiments that demonstrate the execution of complex tasks using MAPS².

Additionally, the algorithm presented is sound, meaning that it produces a trajectory that meets the STL formula and is probabilistically complete, meaning that it will find such a trajectory if one exists.

In our prior study (Sewlia et al., 2023), we addressed the STL motion-planning problem for two coupled agents. There, we extended the conventional Rapidly exploring Random Trees (RRT) algorithm to sample in both the time and space domains. Our approach incrementally built spatio-temporal trees through which we enforced space and time constraints as specified by the STL formula. The algorithm employed a sequential planning method, wherein each agent communicated and waited for the other agent to build its tree. In contrast, the present work addresses the STL motion-planning problem for multiple robots. Here, our algorithm adopts a distributed optimisation-based approach, where spatial and temporal aspects are decoupled to satisfy the STL formula. Instead of constructing an incremental tree, as done in the previous work, we introduce a novel metric called the validity domain and initialise the process with an initial trajectory. In the current research, we only incorporate the STL parse tree and the Satisfaction variable tree from our previous work. Additionally, we present experimental validation results and introduce a novel STL verification architecture.

The rest of the paper is organised as follows. The related work is presented next. Then preliminaries and problem formulation are discussed, followed by the decomposition of the STL formula into temporal and spatial constraints. Next, the main algorithm MAPS² with its analyses is presented. Afterwards, simulations of various robotics tasks are shown, followed by experimental validation on a real multi-robot setup. Finally, the conclusion is presented.

Related work

In the domain of single-agent motion planning, different algorithms have been proposed to generate safe paths for robots. Sampling-based algorithms, such as CBF-RRT (Yang et al., 2019), have achieved success in providing a solution to the motion-planning problem in dynamic environments. However, they do not consider high-level complex mission specifications. Works that impose high-level specifications in the form of LTL, such as Ayala et al. (2013); Bhatia et al. (2010); Fainekos et al. (2009); Vasile and Belta (2013), resort to a hybrid hierarchical control regime resulting in abstraction and explosion of state-space. While a mixed-integer program can encode this problem for linear systems and linear predicates (Wolff et al., 2014), the resulting algorithm has exponential complexity, making it impractical for high-dimensional systems, complex specifications, and long duration tasks. To address this issue, Kurtz and Lin (2022) propose a more efficient encoding for STL to reduce the exponential complexity in binary variables. Additionally, Lindemann and Dimarogonas (2017) introduce a new metric, discrete average space robustness, and composes a Model Predictive Control (MPC) cost function for a subset of STL formulas.

In multi-agent temporal logic control, works such as Verginis and Dimarogonas (2018) and Kress-Gazit et al. (2009) employ workspace discretisation and abstraction techniques, which we avoid in this article due to it being computationally demanding. Some approaches to STL synthesis involve using mixed-integer linear programming (MILP) to encode constraints, as previously explored in Belta and Sadraddini (2019), Raman et al. (2014), and Sadraddini and Belta (2015). However, MILPs are computationally intractable when dealing with complex specifications or long-term plans because of the large number of binary variables required in the encoding process. The work in Sun et al. (2022) encodes a new specification called multi-agent STL (MA-STL) using mixed-integer linear programs (MILP). However, the predicates here depend only on the states of a single agent, can only represent polytope regions, and finally, temporal operations can only be applied to a single agent at a time. In contrast, this work explores coupled constraints between robots and predicates are allowed to be of nonlinear nature.

As a result, researchers have turned to transient control-based approaches such as gradient-based, neural network-based, and control barrier-based methods to provide algorithms to tackle the multi-robot STL satisfaction problem (Kurtz and Lin, 2022). Such approaches, at the cost of imposing dynamical constraints on the optimisation problem, often resort to using smooth approximations of temporal operators at the expense of completeness arguments or end-up considering only a smaller fragment of the syntax (Lindemann et al., 2017; Charitidou and Dimarogonas, 2021; Chen and Dimarogonas, 2022; Lindemann and Dimarogonas, 2018). STL’s robust semantics are used to construct cost functions to convert a synthesis problem to an optimisation problem that benefits from gradient-based solutions. However, such approaches result in non-smooth and non-convex problems and solutions are prone to local minima (Gilpin et al., 2020). In this work, we avoid approximations and consider the full expression of the STL syntax. The proposed solution adopts a purely geometrical approach to the multi-robot STL planning problem. Our current focus is directed towards the planning problem, specifically the generation of trajectories that fulfil STL constraints, rather than the dynamical constraints or the precise control techniques used to execute the trajectory.

Notations: The set of natural numbers is denoted by $\mathbb{N}$ and the set of real numbers by $\mathbb{R}$ . With $n\in \mathbb{N}$ , ${\mathbb{R}}^n$ is the set of n-coordinate real-valued vectors and ${\mathbb{R}}_{+}^n$ is the set of real n-vector with non-negative elements. The cardinality of a set A is denoted by |A|. If $a\in \mathbb{R}$ and $[b,c]\in {\mathbb{R}}^2$ , the Kronecker sum is defined as $a\oplus [b,c]=[a+b,a+c]\in {\mathbb{R}}^2$ . We further define the Boolean set as $\mathbb{B}=\{\top ,\perp \}$ (True, False). The acronym DOF stands for degrees of freedom.

Preliminaries and problem formulation

In this section, we start by introducing STL and STL parse tree, followed by the problem formulation.

Signal Temporal Logic (STL)

Let $\mathbf{x}:{\mathbb{R}}_{+}\to {\mathbb{R}}^n$ be a continuous-time signal. Signal temporal logic (Maler and Nickovic, 2004) is a predicate-based logic with the following syntax:

$\phi =\top |{\mu }^h|\neg \phi |{\phi }_1{\mathcal{U}}_{[a,b]}{\phi }_2|{\phi }_1\wedge {\phi }_2$

(1)where φ₁, φ₂ are STL formulas and ${\mathcal{U}}_{[a,b]}$ encodes the operator until, with 0 ≤ a < b < ∞; μ ^h is a predicate of the form ${\mu }^h:{\mathbb{R}}^n\to \mathbb{B}$ defined by means of a vector-valued predicate function $h:{\mathbb{R}}^n\to \mathbb{R}$ as

${\mu }^h=\left\{\begin{array}{cc}\top \hfill & h(\mathbf{x}(t))\le 0\hfill \\ \perp \hfill & h(\mathbf{x}(t))>0\hfill \end{array}\right..$

(2)The satisfaction relation (x, t)⊧φ indicates that signal x satisfies φ at time t and is defined recursively as follows:

$\begin{array}{lll}\hfill & (\mathbf{x},t)\vDash {\mu }^h\hfill & \iff h(\mathbf{x}(t))\le 0\hfill \\ \hfill & (\mathbf{x},t)\vDash \neg \phi \hfill & \iff \neg ((\mathbf{x},t)\vDash \phi )\hfill \\ \hfill & (\mathbf{x},t)\vDash {\phi }_1\wedge {\phi }_2\hfill & \iff (\mathbf{x},t)\vDash {\phi }_1\wedge (\mathbf{x},t)\vDash {\phi }_2\hfill \\ \hfill & (\mathbf{x},t)\vDash {\phi }_1{\mathcal{U}}_{[a,b]}{\phi }_2\hfill & \iff \exists t_1\in [t+a,t+b]\text{\hspace{0.17em}s.t.\hspace{0.17em}}(\mathbf{x},t_1)\vDash {\phi }_2\hfill \\ \hfill & \hfill \wedge \forall t_2\in [t,t_1],(\mathbf{x},t_2)\vDash {\phi }_1.& \hfill \end{array}$

We also define the operators disjunction, eventually, and always as φ₁ ∨ φ₂ ≡ ¬(¬φ₁ ∧¬φ₂), ${\mathcal{F}}_{[a,b]}\phi \equiv \top {\mathcal{U}}_{[a,b]}\phi$ , and ${\mathcal{G}}_{[a,b]}\phi \equiv \neg {\mathcal{F}}_{[a,b]}\neg \phi$ , respectively. Each STL formula is valid over a time horizon defined as follows.

Definition 1

(Madsen et al., 2018). The time horizon th(φ) of an STL formula φ is recursively defined as

$\mathrm{t}\mathrm{h}(\phi )=\left\{\begin{array}{cc}0,\hfill & \text{if}\phi =\mu \hfill \\ \mathrm{t}\mathrm{h}({\phi }_1),\hfill & \text{if}\phi =\neg {\phi }_1\hfill \\ \max \{\mathrm{t}\mathrm{h}({\phi }_1),\mathrm{t}\mathrm{h}({\phi }_2)\},\hfill & \text{if}\phi ={\phi }_1\wedge {\phi }_2\hfill \\ b+\max \{\mathrm{t}\mathrm{h}({\phi }_1),\mathrm{t}\mathrm{h}({\phi }_2)\},\hfill & \text{if}\phi ={\phi }_1{\mathcal{U}}_{[a,b]}{\phi }_2.\hfill \end{array}\right.$

(3)

In this work, we consider only time bounded temporal operators, that is, when th(φ) < ∞. In the case of unbounded STL formulas, it is only possible to either falsify an always operator or satisfy an eventually operator in finite time, thus we consider only bounded time operators. Next, we state a common assumption regarding the STL formula.

Assumption 1

The STL formula is in positive normal form, that is, it does not contain the negation operator.

The above assumption does not cause any loss of expression of the STL syntax (1). As shown in Sadraddini and Belta (2015), any STL formula can be written in positive normal form by moving the negation operator to the predicate.

STL parse Tree

An STL parse tree is a tree representation of an STL formula (Sewlia et al., 2023). It can be constructed as follows:

• Each node is either a temporal operator node $\{{\mathcal{G}}_I,{\mathcal{F}}_I\}$ , a logical operator node {∨, ∧, ¬}, or a predicate node {μ ^h }, where $I\subset \mathbb{R}$ is a closed interval;

• temporal and logical operator nodes are called set nodes;

• a root node has no parent node and a leaf node has no child node. The leaf nodes constitute the predicate nodes of the tree.

A path in a tree is a sequence of nodes that starts at a root node and ends at a leaf node. The set of all such paths constitutes the entire tree. A subpath is a path that starts at a set node and ends at a leaf node; a subpath could also be a path. The resulting formula from a subpath is called a subformula of the original formula. In the following, we denote any subformula of an STL formula φ by $\overline{\phi }$ . Each set node is accompanied by a satisfaction variable $\tau :\overline{\phi }\to \{+1,-1\}$ and each leaf node is accompanied by a predicate variable π = μ ^h where h is the corresponding predicate function. A signal x satisfies a subformula $\overline{\phi }$ if τ = +1 corresponding to the set node where the subpath of $\overline{\phi }$ begins. Similarly, τ(root(φ)) = +1 ⇔ (x, t) ⊧ φ where root is the root node of φ. An analogous tree of satisfaction and predicate variables can be drawn, called satisfaction variable tree. The satisfaction variable tree borrows the same tree structure as the STL parse tree. Each set node from the STL parse tree maps uniquely to a satisfaction variable τ _i and each leaf node maps uniquely to a predicate variable π _i , where i is an enumeration of the nodes in the satisfaction variable tree. An example of construction of such trees is shown below.

Example 1

The STL parse tree and the satisfaction variable tree for the STL formula

$\phi ={\mathcal{F}}_{I_1}\left({\mu }^{h_1}\vee {\mathcal{G}}_{I_2}({\mu }^{h_2})\right)\wedge {\mathcal{G}}_{I_3}{\mathcal{F}}_{I_4}({\mu }^{h_3})\wedge {\mathcal{G}}_{I_5}({\mu }^{h_4}).$

(4)are shown in Figure 1. From the trees, one obtains the implications ${\tau }_2=+1\Rightarrow (\mathbf{x},t)\vDash {\mathcal{F}}_{I_1}\left({\mu }^{h_1}\vee {\mathcal{G}}_{I_2}({\mu }^{h_2})\right)$ , and ${\tau }_7=+1\Rightarrow (\mathbf{x},t)\vDash {\mathcal{G}}_{I_5}({\mu }^{h_4})$ .

Problem formulation

We consider a team of N robots, where each robot has state ${\mathbf{x}}_i\in {\mathcal{W}}_i\subset {\mathbb{R}}^{n_i}$ , i ∈ {1, …, N} and n _i is the number of degrees of freedom of robot i. The overall state vector is then $\mathbf{x}≔{[{\mathbf{x}}_1^{\top }\cdots {\mathbf{x}}_{\mathbf{N}}^{\top }]}^{\top }$ evolving in a workspace $\mathcal{W}={\mathcal{W}}_1\times \cdots \times {\mathcal{W}}_{\mathbf{N}}$ and we denote by n = n₁ + ⋯ + n _N the number of degrees of freedom of the multi-robot system. We consider the STL formula of the form (1) with a total of K predicates,

${\mu }^{h^{(k)}}=\left\{\begin{array}{cc}\top \hfill & h^{(k)}(\mathbf{x}(t))\le 0\hfill \\ \perp \hfill & h^{(k)}(\mathbf{x}(t))>0\hfill \end{array}\right.,k=1,\dots ,K.$

(5)Before we present the problem statement, we will introduce the multi-robot STL notation and the communication structure of the multi-robot system. In this direction, define a support set for each predicate function h^(k)(x) that captures all the robots upon which the predicate function imposes constraints. Define a projection matrix $E_i\in {\mathbb{R}}^{n\times n_i}$ such that $E_i^{\top }\mathbf{x}={\mathbf{x}}_i$ . The matrix E _i takes the form,

$E_i=\left[\begin{array}{c}\hfill {\mathbf{0}}_{n_1}\hfill \\ \hfill ⋮\hfill \\ \hfill I_{n_i}\hfill \\ \hfill ⋮\hfill \\ \hfill {\mathbf{0}}_{n_{\mathbf{N}}}\hfill \end{array}\right]\in {\mathbb{R}}^{n\times n_i},$

that is, $E_{i}v={\left[\begin{array}{c}\hfill 0,\dots ,0,v,0\dots ,0\hfill \end{array}\right]}^{\top }$ inserts any $v\in {\mathbb{R}}^{n_i}$ at the i-th position of the vector. The support set is then defined for each predicate function h^(k) as

$\begin{array}{l}\hfill {\mathbf{S}}_k≔\{i\in \mathcal{I}|\exists \mathbf{x}\in {\mathbb{R}}^n,v\in {\mathbb{R}}^{n_i},ϵ>0\\ \hfill :h^{(k)}(\mathbf{x}+ϵE_{i}v)\ne h^{(k)}(\mathbf{x})\}.\end{array}$

The support set thus captures all robots i for which some perturbation confined to their own state x _i can change the value of the predicate function h^(k). If a predicate function h^(k) imposes constraints on the states of multiple robots, that is, |S _k | ≥ 2, then we say that the predicate function is coupled.

Example 2

For an STL formula φ = (‖x₁ − x₂‖ > 5) ∧ (‖x₂ − x₃‖ < 2), h⁽¹⁾ = 5 − ‖x₁ − x₂‖ and h⁽²⁾ = ‖x₂ − x₃‖ − 2. Then S₁ = {1, 2} and S₂ = {2, 3}.

Let K _c ≤ K denote the number of coupled predicate functions, and let these be indexed as $h_j^c$ , where j ∈ {1, …, K _c }. For each j, there exists an index k _j ∈ {1, …, K} such that $h_j^c=h^{(k_j)}$ ; let ${\mathbf{S}}_{k_j}\subseteq \{1,\dots ,\mathbf{N}\}$ denote the support set of $h_j^c$ , that is, the set of robots whose states appear in $h_j^c$ . Each robot $i\in {\mathbf{S}}_{k_j}$ then uses a local copy $h_{i,j}^c≔h_j^c$ . This local copy is provided to the robots manually offline prior the start of the algorithm. The robots can also obtain the functions $h_{i,l}^d$ and $h_{i,j}^c$ that are coupled to their own states directly from φ. This could, for example, be done by defining a regular expression(regex) pattern and extracting predicate functions that involve the states x _i (Goyvaerts, 2016).

Next, let L _i be the number of independent predicate functions that involve only the states of robot i, with ${\sum }_{i=1}^{\mathbf{N}}{L}_i=K-K_c$ . Such predicate functions are indexed as $h_{i,l}^d$ , i ∈ {1, …, N} and l ∈ {1, …, L _i }.

Let ${\mathrm{\Pi }}_{k_j}$ be the projection that keeps the indicated state components of the robots in ${\mathbf{S}}_{k_j}$ only, that is, ${\mathrm{\Pi }}_{k_j}\mathbf{x}≔\{E_m^{\top }\mathbf{x}|m\in {\mathbf{S}}_{k_j}\}$ . The predicate function constraints for each robot i are then defined as follows:

$h_{i,l}^d({\mathbf{x}}_i)\le 0\text{and}{h}_{i,j}^c({\mathrm{\Pi }}_{k_j}\mathbf{x})\le 0$

(6)for all i ∈ {1, …, N}, l ∈ {1, …, L _i } and $j\in \{1,\dots ,K_c|i\in {\mathbf{S}}_{k_j}\}$ . The coupled predicate functions $h_{i,j}^c$ can reflect physical interactions between the robots if the constraint is such, for example, if $h_{i,j}^c$ specifies an obstacle avoidance constraint or an object handover task.

Example 3

Consider the STL formula,

$\begin{array}{ll}\hfill \phi =(\Vert x_1-x_2\Vert <1)& \wedge (\Vert x_2-x_3\Vert <1)\hfill \\ \hfill (\Vert x_3-x_4\Vert <1)& \wedge (\Vert x_4-x_1\Vert <1)\hfill \\ \hfill (\Vert x_1\Vert <1)& \wedge (\Vert x_2-x_3-x_4\Vert <1)\hfill \\ \hfill (\Vert x_2\Vert <1)& \wedge (\Vert x_5\Vert <1)\hfill \end{array}$

For the above STL formula, K = 8 and K _c = 5. The number of independent predicates for each robot are L₁ = 1, L₂ = 1, L₃ = L₄ = 0 and L₅ = 1. Table 1 depicts the labelled predicates for the above STL formula. The subscript i in the labels of the predicate functions $h_{i,l}^d$ and $h_{i,j}^c$ specifies the robot responsible for the predicate function.

Now we are ready to define the communication structure of the multi-robot system, which is dictated by a graph. Let the communication graph be given by G = (V, E) where V is the set of vertices corresponding to the indices of the robots and E is the set of edges. In particular, the edge (i, j) ∈ E indicates that robot i can communicate with robot j as the subsequent assumption states.

Table 1.

Labelled predicates for the example STL formula.

‖x1 − x2‖ − 1	$h_{\mathrm{1,1}}^c,h_{\mathrm{2,1}}^c$
‖x2 − x3‖ − 1	$h_{\mathrm{2,2}}^c,h_{\mathrm{3,1}}^c$
‖x3 − x4‖ − 1	$h_{\mathrm{3,2}}^c,h_{\mathrm{4,1}}^c$
‖x4 − x1‖ − 1	$h_{\mathrm{4,2}}^c,h_{\mathrm{1,2}}^c$
‖x1‖ − 1	$h_{\mathrm{1,1}}^d$
‖x2 − x3 − x4‖ − 1	$h_{\mathrm{2,3}}^c,h_{\mathrm{3,3}}^c,h_{\mathrm{4,3}}^c$
‖x2‖ − 1	$h_{\mathrm{2,1}}^d$
‖x5‖ − 1	$h_{\mathrm{5,1}}^d$

Assumption 2

If (i, j) ∈ E, then robots i and j can continuously communicate with each other.

We further assume that, every coupled predicate function induces an edge, ensuring that all state variables needed to compute the predicate function are locally available for each robot.

Assumption 3

If there exists k ∈ {1, …, K} such that i, j ∈ S _k , for i ≠ j, then (i, j) ∈ E.

Note that the aforementioned assumption implies that the graph is undirected, that is, (i, j) ∈ E implies (j, i) ∈ E. Additionally, based on G, the neighbourhood set ${\mathcal{N}}_i$ of a robot i is defined as ${\mathcal{N}}_i=\{j\in \mathbf{V}|(i,j)\in \mathbf{E}\}$ . We also assume that G is static, meaning that no new vertices are added and no edges are created or deleted. With the above assumptions, we are ready to define the distributed information flow:

Definition 2

An algorithm is called distributed, if it can be executed individually by each robot i (a local version of the algorithm) by using only information from its neighbours ${\mathcal{N}}_i$ .

This definition of distributed algorithm does not allow for any global information sharing among robots that are not neighbours with each other, thus a central computer cannot evaluate an STL formula. For example, consider the STL formula:

$\phi =\left(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert \le 1\right)\wedge \left(\Vert {\mathbf{x}}_2-{\mathbf{x}}_3\Vert \le 1\right)$

where x₁, x₂, and x₃ are the states of robot 1, 2, and 3 respectively. Then, Assumption 3 allows for a communication link between robot 1 and robot 2, and between robot 2 and robot 3. A distributed algorithm, in the sense of our work, does not allow for communication between robot 1 and robot 3.

We are now ready to formally state the problem addressed in this paper.

Problem 1. Given an STL formula φ that specifies tasks on a multi-robot system with N robots, design a distributed algorithm to find continuous time-varying trajectories ${\mathbf{y}}_i:[0,\mathrm{t}\mathrm{h}(\phi )]\to {\mathcal{W}}_i$ , starting at an initial configuration y _i (0) = x _i (0), such that (y, t) ⊧ φ, ∀t ∈ [0, th(φ)] with $\mathbf{y}≔{[{\mathbf{y}}_1^{\top },{\mathbf{y}}_2^{\top },\dots ,{\mathbf{y}}_{\mathbf{N}}^{\top }]}^{\top }$ .

It should be noted that we do not currently address the closed-loop stability of the underlying multi-robot system. Instead, we focus on the trajectory generation aspect and rely on existing low-level control approaches to track the generated trajectories. For more information, see Remark 3. The above problem is addressed assuming that at least one such solution exists. This will help us provide probabilistic completeness guarantees later on. Formally, we state the following assumption:

Assumption 4

There exists at least one y such that (y, t) ⊧ φ.

STL formula decomposition

In this section, we present how to retrieve spatial and temporal constraints from a given STL formula φ.

Spatial constraints

In Section ‘Problem Formulation', we provisioned each predicate function h^(k)(x) appearing in the STL formula φ over the complete multi-robot system to the corresponding robot i, denoting them as $h_{i,l}^d$ and $h_{ij}^c$ , depending om whether the robot is responsible for an independent or a coupled task.

For robot i, cast the constraints (6) into the cost function F ⁱ as

$F^i≔\sum _{l=1}^{L_i}\frac{1}{2}\max {\left(0,h_{i,l}^d\right)}^2+\sum _{\begin{array}{c}j\in \{1,\dots ,K_c\}\\ :i\in {\mathbf{S}}_j\end{array}}\frac{1}{2}\max {\left(0,h_{i,j}^c\right)}^2$

(7)Observe that $F^i:\mathcal{W}\to {\mathbb{R}}_{+}$ and F ⁱ = 0 if and only if all the constraints in (6) are satisfied. Then, enforcing conditions (6) is equivalent to finding x _i for a given ${\mathbf{x}}_j(j\in {\mathcal{N}}_i)$ such that F ⁱ = 0. This problem can be posed as

$\underset{{\mathbf{x}}_i\in {\mathcal{W}}_i}{\min }{F}^i$

(8)whose solution ${\mathbf{x}}_i^{\star }$ satisfies F ⁱ (x^⋆) = 0. In the cost function (7), to reduce computational costs, we only minimise h^(k)(x(t)) when h^(k)(x(t)) > 0 while leaving h^(k)(x(t)) ≤ 0 unchanged. This leads us to minimise: F ⁱ = max(0, h^(k)(x(t))), which results in F ⁱ = h^(k)(x(t)) when h^(k)(x(t)) > 0 and F ⁱ = 0 when h^(k)(x(t)) ≤ 0. Additionally, squaring the function penalises larger errors more than smaller ones. Other cost functions that enforce h^(k)(x(t)) ≤ 0 within the validity domain can also be considered. For example, the cost function $F^i≔{\sum }_{l=1}^{L_i}{h_{i,l}^d}^2+{\sum }_j{h_{i,j}^c}^2$ , j ∈ {1, …, K _c }: i ∈ S _j could also be used. However, it was not our first choice, as we aimed to minimise h^(k)(x(t)) only when h^(k)(x(t)) > 0, whereas this cost function would attempt to minimise h^(k)(x(t)) regardless of its sign. Additionally, our formulation is general and works for any type of STL formula as any type of objectives can be encoded in the STL formula, from which, we can extract the predicate functions and minimise the function F ⁱ .

The solution for finding the global minimum of a non-convex function is a subject of extensive research. We argue that employing gradient descent with random initialisations is adequate for addressing this problem, particularly since the initialisations are sampled from a compact set, ${\mathcal{W}}_i$ . Furthermore, using the knowledge that the minimum of the function, F ⁱ (x) = 0, acts as a stopping criterion and facilitates the attainment of the desired solution. We direct readers to the seminal work in Nedic and Ozdaglar (2009), which presents a distributed gradient descent algorithm for multi-agent systems. Additionally, under certain assumptions, Daneshmand et al. (2020) demonstrate that gradient descent with a constant step size avoids entrapment at saddle points. Gradient descent is also shown to efficiently manage most reach-avoid constraints without the need for re-initialisation, given that such constraints are expressible using norms. For our application of gradient descent, we utilise the function presented in Function 1. Function 1 implements the gradient descent

algorithm as described in Algorithm 9.3 of Boyd and Vandenberghe (2004), utilising initial conditions x _i , step size δ, maximum number of iterations L′, and activation variables λ _ij as inputs. The activation variables are presented later in Function 3. It returns the optimised states ${\mathbf{x}}_i^{\star }$ as output. In line 4, the function GradientComputation() computes the gradient, either analytically or numerically. The stopping criterion is met either when a feasible state is determined, indicated by F ⁱ = 0, or when the iteration count exceeds L′ (line 4), which may occur due to multiple conflicting predicates active within F ⁱ . This situation arises because the algorithm accounts for the possibility that the eventually operator may not be satisfied at every sampled point within its validity domain. This occurs, for example, if $\phi ={\mathcal{F}}_{[\mathrm{0,5}]}{\mathcal{G}}_{[\mathrm{0,5}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)\wedge {\mathcal{G}}_{[\mathrm{5,10}]}(g^{(2)}(\mathbf{x})\le {ϵ}_2)$ , and there is a conflict between h⁽¹⁾(x) ≔ g⁽¹⁾(x) − ϵ₁ and h⁽²⁾(x) ≔ g⁽²⁾(x) − ϵ₂ (i.e. $\nexists {\mathbf{x}}_i^{\star }\in {\mathcal{W}}_i$ such that $h^{(1)}(\mathbf{x})({\mathbf{x}}_i^{\star })\le 0\wedge h^{(2)}(\mathbf{x})({\mathbf{x}}_i^{\star })\le 0$ ). In such cases, it becomes necessary for h⁽¹⁾(x) ≤ 0 to be true exclusively within the interval [0,5][s] and for h⁽²⁾(x) ≤ 0 to be true exclusively within the interval [5,10][s].

The robots solve their respective optimisation problem cooperatively in a distributed manner via inter-neighbour communication. This makes the problem distributed, as every interaction between robots is part of the communication graph. Given the nature of the optimisation problem, there is a trade-off between robustness and optimisation performance since x^⋆ converges to the boundaries imposed by the STL formula constraints, making it vulnerable to potential perturbations. However, introducing a slack variable into the equation can enhance robustness, albeit at the cost of sacrificing completeness arguments. The example below shows how to construct the optimisation functions F ⁱ .

Example 4

Consider a system with 3 agents and the corresponding states {x₁, x₂, x₃}, and let the STL formula be: φ = (‖x₁ − x₂‖ > 5) ∧ (‖x₂ − x₃‖ < 2); then, the functions F ⁱ , for i ∈ {1, 2, 3}, are

$\begin{array}{ll}\hfill F^1& =\frac{1}{2}\max {(\mathrm{0,5}-\Vert x_1-x_2\Vert )}^2\hfill \\ \hfill F^2& =\frac{1}{2}\max {(\mathrm{0,5}-\Vert x_1-x_2\Vert )}^2+\frac{1}{2}\max {(0,\Vert x_2-x_3\Vert -2)}^2\hfill \\ \hfill F^3& =\frac{1}{2}\max {(0,\Vert x_2-x_3\Vert -2)}^2.\hfill \end{array}$

Now that spatial constraints are encoded into the optimisation problem, we are ready to encode temporal constraints in the following section, thus completing our STL decomposition into spatial and temporal constraints.

Temporal constraints

We now introduce the concept of validity domain, a time interval associated with every predicate and defined for every path in the STL formula. This interval represents the time domain over which each predicate applies and is defined as follows:

Definition 3

The validity domain $\mathrm{v}\mathrm{d}(\overline{\phi })$ of each path $\overline{\phi }$ of an STL formula φ, is recursively defined as

$\mathrm{v}\mathrm{d}(\overline{\phi })=\left\{\begin{array}{cc}0,\hfill & \text{if}\overline{\phi }={\mu }^h\hfill \\ \mathrm{v}\mathrm{d}({\overline{\phi }}_1),\hfill & \text{if}\overline{\phi }=\neg \overline{{\phi }_1}\hfill \\ [a,b],\hfill & \text{if}\overline{\phi }={\mathcal{G}}_{[a,b]}{\mu }^h\hfill \\ a\oplus \mathrm{v}\mathrm{d}({\overline{\phi }}_1),\hfill & \text{if}\overline{\phi }={\mathcal{G}}_{[a,b]}{\overline{\phi }}_1,{\overline{\phi }}_1\ne {\mu }^h\hfill \\ t^{\star }+T^{\star }\oplus \mathrm{v}\mathrm{d}({\overline{\phi }}_1),\hfill & \text{if}\overline{\phi }={\mathcal{F}}_{[a,b]}{\overline{\phi }}_1\hfill \end{array}\right.$

(9)where $T^{\star }≔\{t\in [a,b]|(\mathbf{x},t)\vDash {\mathcal{F}}_{[a,b]}\overline{\phi }\}$ is a time instant in [a, b] when the state x evaluated at t of a signal x(t) satisfies the eventually operator. The variable t^⋆ is initialised to 0, but takes the value t^⋆ = T^⋆ every time T^⋆ is updated and thus captures the last instance of satisfaction for the eventually operator.

The above definition of t^⋆ is necessary due to the redundancy of the eventually operator; we must ascertain the specific instances where the eventually condition is met to ensure finding a feasible trajectory. Additionally, we need to maintain the history of T^⋆ for nested temporal operators which require recursive satisfaction. The validity domain is determined for each path of an STL formula in a hierarchical manner, beginning at the root of the tree, and each path has a distinct validity domain. The number of leaf nodes in an STL formula is equal to the total number of validity domains. In Definition 3, we do not include the operators ∧ and ∨ because they do not impose temporal constraints on the predicates and thus inherit the validity domains of their parent node. If there is no parent node, operators ∧ and ∨ inherit the validity domains of their child node.

Remark 1

The validity domain is specially defined in the following cases. If a path contains only predicates, the validity domain of μ ^h is equal to the time horizon of φ (i.e. vd(μ ^h ) = th(φ)). Furthermore, if a path contains nested formulas with the same operators, such as $\overline{\phi }={\mathcal{G}}_{[\mathrm{1,10}]}{\mathcal{G}}_{[\mathrm{0,2}]}(\cdot )$ , then the validity domain of $\overline{\phi }$ is equal to the time horizon of the path i.e. $\mathrm{v}\mathrm{d}(\overline{\phi })=\mathrm{t}\mathrm{h}(\overline{\phi })$ . For example, $\mathrm{v}\mathrm{d}({\mathcal{G}}_{[\mathrm{1,10}]}{\mathcal{G}}_{[\mathrm{0,2}]}(\cdot ))=\mathrm{t}\mathrm{h}(\overline{\phi })=[\mathrm{1,12}]$ .

Example 5

Consider the following examples of the validity domain:

• ${\phi }_1={\mathcal{G}}_{[\mathrm{5,10}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ , then vd(φ₁) = [5, 10], which is the interval over which ${\mu }^{h^{(1)}}$ must hold. Here ${\mu }^{h^{(1)}}$ is the predicate corresponding to the predicate function h⁽¹⁾(x) = g⁽¹⁾(x) − ϵ₁.

• ${\phi }_2={\mathcal{F}}_{[\mathrm{5,10}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ , then t^⋆ is initialised to 0, T^⋆ ∈ [5, 10] and $\mathrm{v}\mathrm{d}({\mu }^{h^{(1)}})=0$ . Therefore, vd(φ₂) = T^⋆ ∈ [5, 10] is the instance when ${\mu }^{h^{(1)}}$ must hold.

• ${\phi }_3={\mathcal{F}}_{[\mathrm{5,10}]}{\mathcal{G}}_{[\mathrm{0,2}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ , then t^⋆ is initialised to 0, T^⋆ ∈ [5, 10], $\mathrm{v}\mathrm{d}({\mathcal{G}}_{[\mathrm{0,2}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1))=[\mathrm{0,2}]$ . Therefore, vd(φ₃) = 0 + T^⋆ ⊕ [0, 2] = [T^⋆, T^⋆ + 2] is the interval over which ${\mu }^{h^{(1)}}$ must hold such that φ₃ is satisfied.

• ${\phi }_4={\mathcal{G}}_{[\mathrm{2,10}]}{\mathcal{F}}_{[\mathrm{0,5}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ , then a = 2 and $\mathrm{v}\mathrm{d}({\phi }_4)=2\oplus \mathrm{v}\mathrm{d}({\mathcal{F}}_{[\mathrm{0,5}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1))=2+0+T^{\star }$ where T^⋆ ∈ [0, 5]. Suppose T^⋆ = 1, then vd(φ₄) = 3 is the time instance when ${\mu }^{h^{(1)}}$ must hold. Once ${\mu }^{h^{(1)}}=\top$ , then t^⋆ = T^⋆ and the new vd(φ₄) = 2 + 1 + T^⋆ where T^⋆ ∈ [0, 5].

• ${\phi }_5={\mathcal{F}}_{[\mathrm{0,100}]}{\mathcal{G}}_{[\mathrm{5,10}]}{\mathcal{F}}_{[\mathrm{0,1}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ , then t^⋆ = 0, T^⋆ ∈ [0, 100] and $\mathrm{v}\mathrm{d}({\phi }_5)=T^{\star }+a\oplus \mathrm{v}\mathrm{d}({\mathcal{F}}_{[\mathrm{0,1}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1))$ . Suppose T^⋆ = 50, then $\mathrm{v}\mathrm{d}({\phi }_5)=55\oplus \mathrm{v}\mathrm{d}({\mathcal{F}}_{[\mathrm{0,1}]}(g^{(1)}(\mathbf{x})\le {ϵ}_1))$ and so on.

Regarding the STL formula in equation (4), the validity domains are defined for the following paths: ${\mathcal{F}}_{I_1}\to {\mu }^{h^1},{\mathcal{F}}_{I_1}\to {\mathcal{G}}_{I_2}\to {\mu }^{h^2},{\mathcal{G}}_{I_3}\to {\mathcal{F}}_{I_4}\to {\mu }^{h^3}$ , and ${\mathcal{G}}_{I_5}\to {\mu }^{h^4}.$

We use the following notational convenience in this work: if a parent node of a leaf node of a path $\overline{\phi }$ is an eventually operator we denote the corresponding validity domain by vd ^F (), and, if the parent node of a leaf node of a path $\overline{\phi }$ is an always operator we denote the corresponding validity domain by vd ^G (). The notation vd ^F () indicates that the predicate of the respective leaf node needs to hold at some instance in the said interval, and vd ^G () indicates that the predicate of the respective leaf node needs to hold throughout the interval. The following lemma formalises the relation between the STL formula and its corresponding encoding as described above.

Lemma 1

Suppose $\mathbf{x}(t)=[{\mathbf{x}}_1^{\top },{\mathbf{x}}_2^{\top },\dots ]$ represents the states of all robots, and ${\overline{\phi }}_k$ encompasses all subformulas associated with the STL formula φ. Let ∑ _i F ⁱ (x(t)) = 0 for all $t\in {\bigcup }_k\mathrm{v}\mathrm{d}({\overline{\phi }}_k)$ . Then, it holds that y(t) ⊧ φ.

Proof

The proof follows from the construction of the optimisation function (7) and the validity domain. Notice that if the optimisation problem (8) converges to the desired minima at F ⁱ (x) = 0, then ${\mu }^{h_{i,l}^d}=\top$ and ${\mu }^{h_{i,j}^c}=\top$ for all l ∈ {1, …, L _i } and j ∈ {1, …, K _c }: i ∈ S _j . Next, by definition, the validity domain is defined for the STL formula and if F ⁱ is minimised during the validity domain, then y(t) ⊧ φ.

In the next Section, we present how to integrate the validity domain with the optimisation problem in (8), completing thus the spatial and temporal integration.

Main results

In this section, we present the algorithm for generating continuous trajectories that meet the requirements of a given STL formula φ. The algorithm is executed by the robots offline in a distributed manner, in the sense that they only communicate with their neighbouring robots. The algorithm builds a tree ${\mathcal{T}}_i=\{{\mathcal{V}}_i,{\mathcal{E}}_i\}$ for robot i where ${\mathcal{V}}_i$ is the vertex set and ${\mathcal{E}}_i$ is the edge set. Each vertex $z\in {\mathbb{R}}_{+}\times {\mathcal{W}}_i$ is sampled from a space-time plane. Until now, we denoted the states of robot i as x _i , but from here onward, we denote them as x ⁱ .

In what follows, we give a high-level description of the algorithm. The general idea is to start with an initial trajectory that spans the time horizon of the formula th(φ), then repeatedly sample random points along the trajectory and use gradient-based techniques to find solutions that satisfy the specification at these points. More specifically, the algorithm begins by connecting the initial and final points $z_0^i=\{t_0^i,{\mathbf{x}}_0^i\}$ and $z_f^i=\{t_f^i,{\mathbf{x}}_f^i\}$ with a single edge ${\mathcal{E}}_i=\{(z_0^i,z_f^i)\}$ . The initial conditions $z_0^i=\{t_0^i,{\mathbf{x}}_0^i\}$ depend on the robot’s initial position and time. The final conditions are chosen to be $z_f^i=\{\mathrm{t}\mathrm{h}(\phi )+ϵ,{\mathbf{x}}_f^i\}$ where ϵ > 0 and ${\mathbf{x}}_f^i\in {\mathcal{W}}^i$ . Let $t_0^i=0$ and $t_f^i=\mathrm{t}\mathrm{h}(\phi )+ϵ$ . The final states ${\mathbf{x}}_f^i$ can be randomly chosen since the states in the interval [0, th(φ)] will be determined by the algorithm based on the constraints imposed by φ. The algorithm then randomly selects a time instant t⁰ ∈ [0, th(φ)] and uses linear interpolation to determine the states of each robot at that time, denoted by x⁰. The robots then solve the distributed optimisation problem (8) to find new positions x^⋆ that meet the specification at time t⁰. The algorithm then repeats this process at a user-specified time density, updating the trajectories as necessary. The result is a trajectory that asymptotically improves the task satisfaction of the STL formula.

Example 6

Before we get into the technical details, let us consider an example of 4 agents, represented by the colours blue, green, yellow and magenta, to illustrate the procedure. Suppose, at a specific instance in time, say t⁰, the STL formula requires agent 1 () and agent 2 () to be more than 6 units apart and agent 3 () and agent 4 () to be closer than 6 units, that is, for ϵ > 0,

We begin the process by connecting the initial and final points $z_0^i$ and $z_f^i$ with an initial trajectory for all agents, as shown in Figure 2(a). Each agent’s vertex set is ${\mathcal{V}}_i$ and consists of the start and end points denoted by $z_0^i$ and $z_f^i$ respectively, while its edge set ${\mathcal{E}}_i$ contains only one edge connecting the start and end points. From the initial trajectory, the algorithm randomly selects a point at time instance t⁰ from the entire time domain and uses linear interpolation to determine the state of each agent at that time. The agents solve (8) using the initial position x⁰ to find new position x^⋆, as seen in Figure 2(b). As shown in Figure 2(c), the distributed optimisation problem (8) is solved, resulting in a solution x^⋆, in which agent 1 and agent 2 are positioned so that they are more than 6 units apart and agent 3 and agent 4 remain undisturbed. The latter is the result of using functions of the form $1/2\max {(0,h_{i,j}^c)}^2$ , and since agent 3 and agent 4 already satisfy the requirements, that is, $h_{i,j}^c<0$ , the function is valued 0. The newly determined positions of agents 1 and 2 are added to the tree, allowing the trajectory to be shaped to meet the requirements. The updated trajectory can be seen in Figure 2(d). This process of randomly selecting a point in time, determining the state of the agents and updating their positions is repeated for a user-defined number of times L, to ensure that the trajectory satisfies the STL formula φ throughout the time horizon.

Figure 1.

STL parse tree and satisfaction variable tree for the formula in (4).

MAPS²

The architecture of the algorithm MAPS² (short for ‘multi-robot anytime motion planning under signal temporal logic specifications’), is depicted in Figure 3. The algorithm, outlined in Algorithm 2, begins with an initial trajectory connecting $z_0^i$ and $z_0^f$ , along with a random seed and design constants as input (see lines 1-1). The random seed ensures that all robots select the same time instance. The algorithm proceeds by repeatedly sampling a time instance within the interval, interpolating states at the said time instance, applying gradient descent to minimise the function (7), and either adding or discarding the resulting optimal solution. This process is repeated until the total number of vertices, L, is reached, see lines 3-4. This is also illustrated in Figure 2.

Figure 2.

Illustration of the proposed algorithm.

These steps are implemented as follows: In line 4, the SearchSort() function separates the vertices ${\mathcal{V}}_i$ into two sets based on their time values: one set with time values lower than t⁰ (the vertex with the highest time in this set is indexed with ‘index’), and another with values greater than t⁰ (the vertex with the lowest time in this set is indexed with ‘index + 1’). The corresponding vertices are $z_{\text{index}}^i=\{t_{\text{index}}^i,{\mathbf{x}}_{\text{index}}^i\}$ and $z_{\text{index}+1}^i=\{t_{\text{index}+1}^i,{\mathbf{x}}_{\text{index}+1}^i\}$ . Then, the algorithm uses linear interpolation in line 4 via the function Interpolate() to obtain the vertex $z_{\text{inter}}^i=\{t^0,{\mathbf{x}}_{\text{inter}}^i\}$ . This is obtained by solving for ${\mathbf{x}}_{\text{inter}}^i$ element-wise as the solution of

${\mathbf{x}}_{\text{inter}}^i=\left(\frac{{\mathbf{x}}_{\text{index}+1}^i-{\mathbf{x}}_{\text{index}}^i}{t_{\text{index}+1}^i-t_{\text{index}}^i}\right)(t^0-t_{\text{index}}^i)+{\mathbf{x}}_{\text{index}}^i.$

The vertex $z_{\text{inter}}^i$ is the initial condition to solve the optimisation problem (8); and once a solution $z_{\text{opt}}^i$ is obtained, it is added to the vertex set ${\mathcal{V}}_i$ in line 4. The edge set ${\mathcal{E}}_i$ is reorganised to include $z_{\text{opt}}^i$ in lines 4-4.

Moreover, as a safeguard, if a solution remains undiscovered following L iterations, line 4 initiates a reset procedure. This involves setting the satisfaction variable for all eventually operators back to −1 and restarting the search. Since we assume that at least one viable solution always exists (refer to Assumption 4), the absence of a solution occurs solely when an eventually operator is satisfied at an impractical instance of time. Such an impractical instance of time affects the solution of the algorithm since there are redundancies in picking the satisfaction instance ( $\mathbf{x}(t)\vDash {\mathcal{F}}_{[a,b]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ if h⁽¹⁾(x) ≔ g⁽¹⁾(x) − ϵ₁ ≤ 0 at any single instance in [a, b]). By resetting these operators, the algorithm aims to locate a solution under feasible instances.

GradientDescent

The function is presented in Function 3 and computes the optimal value, $z_{\text{opt}}^i$ , by solving the problem presented in equation (8). This allows the robots to compute vertices that locally do not violate the STL formula. Once $z_{\text{opt}}^i$ is determined through Function 1, the satisfaction variables are updated in Function 4.

Based on the validity domain, the Function 3 determines which predicate functions are active in (7) at every sampled time instance t⁰. The Function ValidityDomain() in line 4 calculates the validity domains of each path $\overline{\phi }$ based on Definition 3. Let K _i be the total number of independent and coupled predicate functions associated with robot i, a binary variable λ _ij ∈ {0, 1}, j ∈ {1, …, K _i } is assigned to determine whether a predicate function is active or not. It is set to 1 if the predicate is active and 0 otherwise. For example,

• If ${\phi }_1={\mathcal{G}}_{[\mathrm{5,10}]}(\Vert x_1-x_2\Vert \le 2)$ , then λ₁₁ = λ₂₁ = 1 whenever t⁰ ∈ [5, 10] and 0 otherwise.

• If ${\phi }_2={\mathcal{F}}_{[\mathrm{10,15}]}(\Vert x_3\Vert \le 5)$ , then λ₃₁ = 1 whenever t⁰ ∈ [10, 15] and 0 otherwise. Once x₃(t) ⊧ φ₂, λ₃₁ = 0 ∀t.

The indices i and j in λ _ij and vd _ij refer to robot i and the jth predicate function associated with robot i, respectively. Here j ∈ {1, …, K _i }. We distinguish three cases: if the sampled point belongs to the validity domain of a single eventually operator and/or a single always operator, λ _ij = 1. If the sampled point belongs to the validity domain of multiple eventually operators, we activate only one of them at random, that is, λ _ij = 1 only for one of them. This avoids enforcing conflicting predicates as it can happen that multiple eventually operators may not be satisfied at the same time instance (e.g. $\phi ={\mathcal{F}}_{[\mathrm{0,1}]}(x>0)\wedge {\mathcal{F}}_{[\mathrm{0,1}]}(x<0)$ ); see lines 5-20.

In lines 21-34, the algorithm updates the satisfaction variable of all paths in the STL formula that impose restrictions on robot i’s states. The algorithm goes bottom-up, starting from the leaf node to the root node. First, it determines if $z_{\mathit{opt}}^i$ is the desired minimum (i.e. $F^i({\mathbf{x}}_{\mathit{opt}}^i)\le 0$ ) in line 27, and in lines 30-33, the algorithm updates the satisfaction variable of all nodes in the path $\overline{\phi }$ through the function SatisfactionVariable(). If $z_{\text{opt}}^i$ is not the desired minimum, then all the satisfaction variables of the path $\overline{\phi }$ are reset to −1 in line 35. This could result from conflicting predicates at the same time instance.

SatisfactionVariable

This function, presented in Function 4, updates the satisfaction variable tree, τ. The aforementioned procedure decides if the satisfaction variable corresponding to each node listed is +1 (satisfied) or −1 (not yet satisfied). The discussion of handling disjunction operators is deferred to Section ‘MAPS²:Branch-and-Pick for Disjunctions', as they are handled differently. Considering the premise that the predicate is true, as indicated in line 27 of Function 3, we evaluate the satisfaction variable as follows:

• ${\mathcal{F}}_I$ : The satisfaction variable of the eventually operator is updated along with t^⋆ = t⁰. This updated t^⋆ is used to determine the new validity domains in line 4 of Function 3; see Example 3 for an illustration of this procedure.

• ${\mathcal{G}}_I$ : Unlike the eventually operator, determining $\tau ({\mathcal{G}}_I)$ necessitates the computation of robustness over the entire validity domain of the operator. The function robust() uses the robust semantics of the STL presented in Maler and Nickovic (2004). Particularly, it samples a user-defined number of points in the interval ${\mathrm{v}\mathrm{d}}_{ij}^G()$ and computes ${\inf }_{t\in {\mathrm{v}\mathrm{d}}_{ij}^G}{h}_{i,l}^d({\mathbf{x}}^i(t))$ or ${\inf }_{t\in {\mathrm{v}\mathrm{d}}_{ij}^G}{h}_{i,j}^c({\mathbf{x}}^i(t))$ . If the robustness is non-negative, indicating satisfaction of the task, the value of $\tau ({\mathcal{G}}_I)$ is updated to +1.

• ∧: This set node returns the satisfaction variable as +1 since it does not impose spatial or temporal restrictions.

Branch-and-pick for disjunctions

In our approach, we address disjunctions as follows: Given an STL formula of the form φ =∨_i∈1,…,Kϕ _i , which can also be represented as φ = ∨(ϕ₁, ϕ₂, …, ϕ _K ), we divide it into K individual STL formulas. The agents then run Algorithm 2 separately for each φ = ϕ _i , where i ∈ 1, …, K. For instance, consider the STL formula represented as (4)

$\phi ={\mathcal{F}}_{I_1}\left({\mu }^{h_1}\vee {\mathcal{G}}_{I_2}({\mu }^{h_2})\right)\wedge {\mathcal{G}}_{I_3}{\mathcal{F}}_{I_4}({\mu }^{h_3})\wedge {\mathcal{G}}_{I_5}({\mu }^{h_4}).$

We branch it into two STL formulas: ${\varphi }_1={\mathcal{F}}_{I_1}{\mu }_1\wedge {\mathcal{G}}_{I_3}{\mathcal{F}}_{I_4}({\mu }_3)\wedge {\mathcal{G}}_{I_5}({\mu }_4)$ and ${\varphi }_2={\mathcal{F}}_{I_1}{\mathcal{G}}_{I_2}({\mu }_2)\wedge {\mathcal{G}}_{I_3}{\mathcal{F}}_{I_4}({\mu }_3)\wedge {\mathcal{G}}_{I_5}({\mu }_4)$ , as illustrated in Figure 4. The search terminates when any branch of the disjunction satisfies the condition τ(root) ≠ + 1, as specified on line 3 of Algorithm 1. We acknowledge that this naive method of handling disjunctions can result in exponential growth with the addition of more operators. An alternative approach, akin to the branch-and-bound method from optimisation (Morrison et al., 2016), involves evaluating the robustness of each ϕ _i for i ∈ 1, …, K and executing MAPS² only for the formulas that show a faster increase in satisfaction. However, this strategy might necessitate a higher level of communication among robots which goes beyond their existing communication network and possibly require a central authority to coordinate task fulfilment. For example, the STL formula

$\phi ={\mathcal{G}}_{[\mathrm{0,5}]}(x_1<5)\vee ({\mathcal{F}}_{[\mathrm{0,5}]}(|x_2-x_3|>2)).$

comprises disjunction between ${\varphi }_1={\mathcal{G}}_{[\mathrm{0,5}]}(x_1<5)$ and ${\varphi }_2={\mathcal{F}}_{[\mathrm{0,5}]}(|x_2-x_3|>2)$ . Observe that ϕ₁ requires no inter-robot communication, while ϕ₂ necessitates communication between robots 2 and 3. In the implementation of a method akin to branch-and-bound, we would branch into two formulas, ϕ₁ and ϕ₂, and repeatedly switch between them if we observe the robustness of one formula decaying faster compared to the other. This switching must be performed by a central authority that observes the decay in robustness. If the switching is decided among the robots, then robot 1 of ϕ₁ needs to communicate the robustness decay with the network of robots 2 and 3. This requires robot 1 to establish communication with the network of robots 2 and 3 in order to decide which branch to grow, thereby necessitating communication links where none existed before. Without such a communication link, both ϕ₁ and ϕ₂ would need to be satisfied using the naive approach presented in our work. This motivates our choice to use the naive approach.

Figure 3.

Architecture of the proposed algorithm.

Analysis

In this section, we analyse the proposed algorithm and arrive at proving the probabilistic completeness.

Let the set $\mathcal{S}\subseteq \mathcal{W}$ be a compact set where a trajectory $\mathbf{y}:[0,\mathrm{t}\mathrm{h}(\phi )]\to \mathcal{S}$ satisfies the STL formula. Along the lines of Kleinbort et al. (2019), let a trajectory y be located on the boundary of the set $\mathcal{S}$ , the satisfiable set, dividing $\mathcal{W}$ into a feasible set $\mathcal{S}$ and an infeasible set $\mathcal{W}\backslash \mathcal{S}$ .

Starting with an initial linear trajectory in the augmented time-space domain, each uniformly sampled time point t⁰ corresponds to a position x_inter either in $\mathcal{S}$ or $\mathcal{W}\backslash \mathcal{S}$ . If ${\mathbf{x}}_{\text{inter}}\in \mathcal{S}$ , we leave it unchanged as it meets the requirements. But if ${\mathbf{x}}_{\text{inter}}\notin \mathcal{S}$ , we use gradient descent to reach a point on y, since it lies on the boundary of the constraints’ set.

Next, divide the trajectory $\mathbf{y}:[0,\mathrm{t}\mathrm{h}(\phi )]\to \mathcal{S}$ into L + 1 points x _k , where 0 ≤ k ≤ L and y(th(φ)) = x _f = x _L by dividing the time duration into equal intervals of δ _t . Without loss of generality, assume that the points x _k and x_k+1 are separated by δ _t in time. With Lδ _t = th(φ), the probability of sampling a point in an interval of length δ _t can be calculated as $p=\frac{{\delta }_t}{\mathrm{t}\mathrm{h}(\phi )}$ . If δ _t ≪th(φ), then p < 1/2. Denote the sequential covering class ¹ of trajectory y as $Y_{{\delta }_t}({\mathbf{x}}_k)$ . The length of $Y_{{\delta }_t}({\mathbf{x}}_k)$ is δ _t in the time domain and is centred at x _k . See Figure 5 for reference.A trial is counted as successful if we sample a point t⁰ within the interval δ _t /2 on either side of x _k , that is, within $Y_{{\delta }_t}({\mathbf{x}}_k)$ . If there are L successful trials, the entire trajectory y is covered, and the motion-planning problem is solved. Consider m total samples, where m ≫ L, and treat this as m Bernoulli trials with success probability p since each sample is independent with only two outcomes. We are now ready to state the following lemma.

Figure 4.

Disjunction representation for disjunctive components using STL parse tree.

Lemma 2

Let a constant L and probability p such that $p<\frac{1}{2}$ . Further, let m represent the number of samples taken by the MAPS ² algorithm. Then, the probability that MAPS ² fails to sample a segment after m samples is at most $\frac{(m-L)p}{{(mp-L)}^2}$ .

Proof

The probability of not having L successful trials after m samples can be expressed as

$\mathbf{P}[X_m\le L]=\sum _{k=0}^{L-1}\left(\genfrac{}{}{0.0pt}{}{m}{k}\right)p^k{(1-p)}^{m-k}$

and according to Feller (1968), if $p<\frac{1}{2}$ , we can upper bound this probability as:

$\mathbf{P}[X_m\le L]\le \frac{(m-L)p}{{(mp-L)}^2}.$

As p and L are fixed and independent of m, the expression $\frac{(m-L)p}{{(mp-L)}^2}$ approaches 0 with as m increases, thus completing the proof.

Next, we present a final lemma which helps us prove the probabilistic completeness of the algorithm.

Lemma 3

No sampled point x _k is falsely labelled as satisfying the STL formula φ unless it actually does.

Proof

The algorithm initiates by setting all satisfaction variables, τ, to −1, as inputs to Algorithm 2. These variables are updated in Function 4 designed for evaluating whether τ meets the satisfaction criteria. The function adjusts τ in accordance with the definition of STL operators presented in the preliminaries section, ensuring that updates accurately reflect the satisfaction status. Furthermore, the update to τ(leaf) within Function 3 (referenced at line 20) occurs only when the condition F ⁱ ≤ 0 is met. This condition indicates that all active predicates are satisfied by definition. Thus, no satisfaction variable is incorrectly updated.

Next, the paper’s final result is presented, which states that the probability of the algorithm providing an STL formula satisfying trajectory (if one exists) approaches one as the number of samples tends to infinity. This is a desirable property for sampling-based planners and such algorithms are termed probabilistically complete.

Theorem 1

Algorithm 2 is probabilistically complete.

Proof

The proof follows from Lemmas 1, 2, and 3. From Lemma 1 and Lemma 3, we know that every sample added to the trajectory satisfies the STL formula. Thus, what needs to be shown is that the algorithm samples infinitely many times and covers the entire time horizon. From Lemma 2, we know that the probability of covering the entire time horizon is 1 − P[X _m ≤ L]. Suppose the Algorithm 2 reaches J = L′ samples without finding a feasible solution, then it discards J samples as seen in line 4 of Algorithm 2. Given Assumption 4, we have J < ∞, and since J is the number of discarded samples, we also have J ≤ m where m is the total number of samples sampled so far (including the discarded ones). Thus, the probability of the trajectory satisfying the STL formula is $1-\frac{((m-J)-L)p}{{((m-J)p-L)}^2}$ , which approaches one as m → ∞. Thus, the algorithm is probabilistically complete.

Remark 2

Our algorithm can be endowed in a post-processing stage with a module that smoothens the trajectory to avoid large accelerations. However, care needs to be taken since the smoothened paths may no longer satisfy the STL formula. One could also use more sophisticated approaches like B-splines to impose velocity and acceleration limits as shown in Lapandić et al. (2024).

Remark 3

At present, our approach does not incorporate kinematic or dynamic constraints. Incorporation of such constraints could be attempted by either deploying the kinodynamic version of the RRT algorithm (Webb and van Den Berg, 2012), or by using an existing low-level controller to track the generated open-loop trajectories. Some examples of such controllers include the Model Predictive Controller (Poignet and Gautier, 2000) and the input constrained Prescribed Performance Controller (Fotiadis and Rovithakis, 2024; Trakas and Bechlioulis, 2023). This incorporation is by no means straightforward but requires fusion with another type of methodological machinery that goes beyond the scope of the current work. Moreover, such controllers have been developed for a large variety of dynamical systems and hence the proposed algorithm is practical and applicable to a large class of robots.

Simulations

In this section, we present simulations of various scenarios encountered in a multi-robot system. Restrictions are imposed using an STL formula and MAPS² is utilised to create trajectories that comply with the STL formula. In the following we consider 4 agents, with δ = 0.1, η = 0.01 and L = L′ = 100. The simulations were run on an 8 core Intel^® Core^TM i7 1.9 GHz CPU with 16 GB RAM. ²

Collision avoidance

We begin with a fundamental requirement in multi-robot systems: avoiding collisions. In this scenario, it is assumed that all agents can communicate or sense each other’s positions. The following STL formula is used to ensure collision avoidance in the interval 20[s] to 80[s]:

$\phi ={\mathcal{G}}_{[\mathrm{20,80}]}(\Vert x_i-x_j\Vert \ge 1)$

where {i, j} ∈ {{1, 2}, {1, 3}, {1, 4}, {2, 3}, {2, 4}, {3, 4}}. As depicted in Figure 6a, all four agents maintain a distance of at least 1 unit from each other during the interval [20,80][s]. The maximum computation time by any agent is 0.1143[s].

Figure 5.

Illustration of $Y_{{\delta }_t}({\mathbf{x}}_k)$ .

Rendezvous

The next scenario is rendezvous. We use the eventually operator to express this requirement. The STL formula specifies that agents 1 and 3 must approach each other within 1 distance unit during the interval [40,60][s] and similarly, agents 2 and 4 must meet at a minimum distance of 1 unit during the same interval. The STL formula is

$\phi ={\mathcal{F}}_{[\mathrm{40,60}]}(\Vert x_1-x_3\Vert \le 1\wedge \Vert x_2-x_4\Vert \le 1).$

As seen in Figure 6b, agents 1 and 3 and agents 2 and 4 approach each other within a distance of 1 unit during the specified interval. It’s worth noting that the algorithm randomly selects the specific time t^⋆ within the continuous interval [40,60][s] at which the satisfaction occurs. The maximum computation time by any agent is 0.0637[s].

Stability

The last task is that of stability, which is represented by the STL formula ${\mathcal{F}}_{[a_1,b_1]}{\mathcal{G}}_{[a_2,b_2]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ . This formula requires that (g⁽¹⁾(x) ≤ ϵ₁) must always hold within the interval [t^⋆ + a₂, t^⋆ + b₂], where t^⋆ ∈ [a₁, b₁]. This represents stability, as it requires (g⁽¹⁾(x) ≤ ϵ₁) to always hold within the interval [t^⋆ + a₂, t^⋆ + b₂], despite any transients that may occur in the interval [a₁, t^⋆). Figure 6c presents a simulation of the following STL formula:

$\begin{array}{ll}\hfill \phi ={\mathcal{F}}_{[\mathrm{0,100}]}{\mathcal{G}}_{[\mathrm{0,20}]}& \left((1.9\le x_1\le 2.1)\wedge (3.9\le x_2\le 4.1)\right.\hfill \\ \hfill & \left.\wedge (5.9\le x_3\le 6.1)\wedge (7.9\le x_4\le 8.1)\right)\hfill \end{array}$

where t^⋆ = 63.97[s]. The maximum computation time by any agent is 0.0211[s].

Recurring tasks

The next scenario is that of recurring tasks. This can be useful when an autonomous vehicle needs to repeatedly survey an area at regular intervals, a bipedal robot needs to plan periodic foot movements, or a ground robot needs to visit a charging station at specified intervals. The STL formula to express such requirements is given by ${\mathcal{G}}_{[a_1,b_1]}{\mathcal{F}}_{[a_2,b_2]}(g^{(1)}(\mathbf{x})\le {ϵ}_1)$ , which reads as ‘beginning at a₁[s], g⁽¹⁾(x) ≤ ϵ₁ must be satisfied at some point in the interval [a₁ + a₂, a₁ + b₂][s] and this should be repeated every [b₂ − a₂][s]’. A simulation of the following task is shown in Figure 6d:

$\phi ={\mathcal{G}}_{[\mathrm{0,100}]}{\mathcal{F}}_{[\mathrm{0,20}]}(\Vert x_1-x_3\Vert \le 1).$

Every 20[s], the condition |x₁ − x₃| ≤ 1 is met. It’s worth noting that the specific time t^⋆ at which satisfaction occurs is randomly chosen by the algorithm. The maximum computation time by any agent is 0.2017[s].

In reference to Remark 2, an example of post-processing the trajectories is shown in Figure 7 for the STL formula,

$\phi ={\mathcal{G}}_{[\mathrm{0,100}]}{\mathcal{F}}_{[\mathrm{0,20}]}(\Vert x_1-x_2\Vert \le 1).$

(10)A 3rd order polynomial was applied using the Savitzky-Golay filter to smoothen the trajectory. Smoothening helps to avoid any large accelerations and sudden velocity changes, though it may come at the cost of potential STL violations.

Figure 6.

Simulation results of MAPS² with four agents

Multi-agent case study

In this case study, we design trajectories for a team of 100 agents that exist in a 100 × 100[m] space and [0,100][s] time span. The team needs to adhere to the following STL formula,

$\phi ={\mathcal{G}}_{[\mathrm{10,90}]}\left[\Vert x_i-x_j\Vert \ge 0.01\wedge \Vert x_i-(\mathrm{50,50})\Vert \le 5\right]$

(11)∀i, j ∈ {1, 2, …, 100} and i ≠ j. Note that the above STL formula has 5150 predicates. In the interval [10,90][s], the STL formula dictates every agent to be at least 0.01[m] apart from every other agent and to be at least 5[m] close to the centre point (50,50)[m]. The simulation results are shown in Figure 8 where the Figures 8(a)-8(c) are the trajectories before the start of the algorithm while Figures 8(d)-8(f) shows the trajectories at the end of j = 1000 iterations, as mentioned in Algorithm 2. The simulation took 17.84[s] to complete without parallelisation. The faster computation can be attributed to the nature of the design of the cost function in (8), which allows for points that do not violate the formula not to be changed. The robustness of the STL formula is shown in Figure 9, a negative robustness signifies task satisfaction. Here, the robustness converges to 0, because the robustness for an always operator reflects the worst-case scenario. It is important to note that computing the result for Figure 9 required 12 hours and 10 minutes of computation time since it had to be performed centrally.

Figure 7.

Non-smooth and smooth paths for the formula (10).

Figure 8.

Simulation of trajectory generation for 100 agents for the STL formula (11).

Overall case study

In this case study, we demonstrate the application of the aforementioned scenarios by setting up the following tasks:

• Agent 1 always stays above 8 units.

• Agents 2 and 4 are required to satisfy the predicate $x_2^2+x_4^2\le 2$ within the time interval [10,30][s].

• Agent 3 is required to track an exponential path within the time interval [20,60][s].

• Agent 2 is required to repeatedly visit Agent 1 and Agent 3 every 10 s within the interval [30,50][s].

• Agent 1 is required to maintain at least 1 unit distance from the other three agents within the interval [80,100][s].

The STL formula for the above tasks is as follows:

$\begin{array}{ll}\hfill \phi & =(x_1\ge 8)\wedge {\mathcal{G}}_{[\mathrm{10,30}]}(x_2^2+x_4^2\le 2)\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{20,60}]}(\Vert x_3-50\exp (-0.1t)\Vert \le 0.05)\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{30,50}]}{\mathcal{F}}_{[\mathrm{0,10}]}\left((\Vert x_2-x_1\Vert \le 0.5)\wedge (\Vert x_2-x_3\Vert \le 0.5)\right)\wedge \hfill \\ \hfill & {\mathcal{F}}_{[79.9,80.1]}{\mathcal{G}}_{[\mathrm{0,20}]}\left((\Vert x_1-x_2\Vert \ge 1)\wedge (\Vert x_1-x_3\Vert \ge 1)\right.\hfill \\ \hfill & \left.\wedge (\Vert x_1-x_4\Vert \ge 1)\right)\hfill \end{array}$

The parameter L was increased to 1000, and η was decreased to 0.001. In Figure 10, we show the resulting trajectories of each agent generated by MAPS² satisfying the above STL formula. The maximum computation time by any agent is 4.611[s].

Figure 9.

Robustness of the STL formula in (11).

Experiments

We now present an experimental demonstration of the proposed algorithm. The multi-robot setup involves three robots, as shown in Figure 11, and consists of 3 mobile bases and two 6-DOF manipulator arms. The locations of the three bases are denoted as ${\mathbf{x}}_1\in {\mathbb{R}}^2$ , ${\mathbf{x}}_2\in {\mathbb{R}}^2$ , and ${\mathbf{x}}_3\in {\mathbb{R}}^2$ , respectively. Base 2 and base 3 are equipped with manipulator arms, whose end-effector positions are represented as ${\mathbf{e}}_1\in {\mathbb{R}}^3$ and ${\mathbf{e}}_2\in {\mathbb{R}}^3$ , respectively.

Figure 10.

Overall case study.

Figure 11.

Experimental setup with three mobile bases and two 6-dof manipulators

The STL formula defining the tasks is the following,

$\begin{array}{ll}\hfill & \phi =\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert \ge 0.6\wedge \Vert {\mathbf{x}}_2-{\mathbf{x}}_3\Vert \ge 0.6\wedge \Vert {\mathbf{x}}_3-{\mathbf{x}}_1\Vert \ge 0.6\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{10,125}]}\Vert {\mathbf{x}}_1-1.8{[-\cos 0.0698t,\sin (0.0698t)]}^{\top }\Vert \le 0.05\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{30,70}]}\Vert {\mathbf{e}}_1-{[{\mathbf{x}}_1^{\top },0.35]}^{\top }\Vert \le 0.01\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{30,70}]}\Vert {\mathbf{x}}_2-1.1{[-\cos 0.0698t,\sin (0.0698t)]}^{\top }\Vert \le 0.05\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{80,120}]}\Vert {\mathbf{e}}_2-{[{\mathbf{x}}_1^{\top },0.35]}^{\top }\Vert \le 0.01\wedge \hfill \\ \hfill & {\mathcal{G}}_{[\mathrm{80,120}]}\Vert {\mathbf{x}}_3-1.1{[-\cos 0.0698t,\sin (0.0698t)]}^{\top }\Vert \le 0.05\wedge \hfill \\ \hfill & {\mathcal{F}}_{[\mathrm{180,200}]}\Vert {\mathbf{x}}_1-{[\mathrm{0,0}]}^{\top }\Vert \le 0.05\wedge \hfill \\ \hfill & {\mathcal{F}}_{[\mathrm{180,200}]}\left(\Vert {\mathbf{x}}_2-[1,-1]\Vert \le 0.05\wedge \Vert {\mathbf{e}}_1-[{\mathbf{x}}_2,0.6]\Vert \le 0.05\right)\wedge \hfill \\ \hfill & {\mathcal{F}}_{[\mathrm{180,200}]}\left(\Vert {\mathbf{x}}_3-[-\mathrm{1,1}]\Vert \le 0.05\wedge \Vert {\mathbf{e}}_2-[{\mathbf{x}}_3,0.6]\Vert \le 0.05\right).\hfill \end{array}$

The above task involves collision avoidance constraints that are always active given by the subformula ${\overline{\phi }}_1=(\Vert {\mathbf{x}}_1-{\mathbf{x}}_2\Vert \ge 0.6)\wedge (\Vert {\mathbf{x}}_2-{\mathbf{x}}_3\Vert \ge 0.6)\wedge (\Vert {\mathbf{x}}_3-{\mathbf{x}}_1\Vert \ge 0.6)$ . Next, in the duration [10,125][s], base 1 surveils the arena and follows a circular time-varying trajectory given by the subformula ${\overline{\phi }}_2=({\mathcal{G}}_{[\mathrm{10,125}]}\Vert {\mathbf{x}}_1-c_1(t)\Vert \le 0.05)$ where c₁(t) is the circular trajectory. In the duration [30,70][s], end-effector 1 tracks a virtual point 0.35[m] over base 1 to simulate a pick-and-place task, given by the subformula ${\overline{\phi }}_3={\mathcal{G}}_{[\mathrm{30,70}]}\Vert {\mathbf{e}}_1-{[{\mathbf{x}}_1^{\top },0.35]}^{\top }\Vert \le 0.01\wedge {\mathcal{G}}_{[\mathrm{30,70}]}\Vert {\mathbf{x}}_2-c_2(t)\Vert \le 0.05$ where c₂(t) is the circular trajectory. Similarly, in the duration [80,120][s], end-effector 2 takes over the task to track a virtual point 0.35[m] over base 1, given by the subformula ${\overline{\phi }}_4={\mathcal{G}}_{[\mathrm{80,120}]}\Vert {\mathbf{e}}_2-{[{\mathbf{x}}_1^{\top },0.35]}^{\top }\Vert \le 0.01\wedge {\mathcal{G}}_{[\mathrm{80,120}]}\Vert {\mathbf{x}}_3-c_2(t)\Vert \le 0.05$ . Finally, eventually in the duration [180,200][s], the robots assume a final position given by the subformula $\overline{{\phi }_5}={\mathcal{F}}_{[\mathrm{180,200}]}\Vert {\mathbf{x}}_1-{[\mathrm{0,0}]}^{\top }\Vert \le 0.05\wedge {\mathcal{F}}_{[\mathrm{180,200}]}\left(\Vert {\mathbf{x}}_2-[1,-1]\Vert \le 0.05\wedge \Vert {\mathbf{e}}_1-[{\mathbf{x}}_2,0.6]\Vert \le 0.05\right)\wedge {\mathcal{F}}_{[\mathrm{180,200}]}\left(\Vert {\mathbf{x}}_3-[-\mathrm{1,1}]\Vert \le 0.05\wedge \Vert {\mathbf{e}}_2-[{\mathbf{x}}_3,0.6]\Vert \le 0.05\right)$ .

The results are shown in Figure 12, where the x-axis represents time in seconds, and the y-axis represents the predicate functions defined by (5). The dashed line in the plots represents the predicate functions of the trajectories obtained by solving the optimisation problem (8), while the solid line represents the predicate functions of the actual trajectories by the robots. In the context of (5), negative values indicate task satisfaction. However, due to the lack of an accurate model of the robots and the fact that the optimisation solution converges to the boundary of the constraints, the tracking is imperfect, and we observe slight violations of the formula by the robots in certain cases. Nonetheless, the trajectories generated by the algorithm do not violate the STL formula. The coloured lines represent the functions that lie within the validity domain of the formula. Figure 12(a) shows that the collision constraint imposed on all 3 bases is not violated, and they maintain a separation of at least 60 cm. In Figure 12(b), base 1 tracks a circular trajectory in the interval [10, 125] seconds. In Figures 12(c) and 12(d), the end effectors mounted on top of bases 2 and 3 track a virtual point over the moving base 1 sequentially. In the last 20 seconds, the bases and end effectors move to their desired final positions, as seen in Figures 12(e) and 12(f). The maximum computation time by any robot is 3.611[s]. Figure 13 shows front-view and side-view at different time instances during the experimental run. ³

Figure 12.

Experimental verification of MAPS² with the setup in Figure 11.

Figure 13.

Front-view and side-view during experimental run with the setup in Figure 11.

Conclusion

This work proposed MAPS², a distributed planner that solves the multi-robot motion-planning problem subject to tasks encoded as STL constraints. By using the notion of validity domain and formulating the optimisation problem as shown in (8), MAPS² transforms the spatio-temporal problem into a spatial planning task, for which efficient optimisation algorithms already exist. Task satisfaction is probabilistically guaranteed in a distributed manner by presenting an optimisation problem that necessitates communication only between robots that share coupled constraints. Extensive simulations involving benchmark formulas and experiments involving varied tasks highlight the algorithms functionality. Future work involves incorporating dynamical constraints such as velocity and acceleration limits into the optimisation problem.

Supplemental Material