Optimal MultiRobot Path Planning on Graphs: Complete Algorithms and Effective Heuristics
Abstract
We study the problem of optimal multirobot path planning on graphs (MPP) over four distinct minimization objectives: the makespan (last arrival time), the maximum (singlerobot traveled) distance, the total arrival time, and the total distance. In a related paper Yu and LaValle (2015), we show that these objectives are distinct and NPhard to optimize. In this work, we focus on efficiently algorithmic solutions for solving these optimal MPP problems. Toward this goal, we first establish a onetoone solution mapping between MPP and networkflow. Based on this equivalence and integer linear programming (ILP), we design novel and complete algorithms for optimizing over each of the four objectives. In particular, our exact algorithm for computing optimal makespan solutions is a first such that is capable of solving extremely challenging problems with robotvertex ratio as high as . Then, we further improve the computational performance of these exact algorithms through the introduction of principled heuristics, at the expense of some optimality loss. The combination of ILP model based algorithms and the heuristics proves to be highly effective, allowing the computation of optimal solutions for problems containing hundreds of robots, densely populated in the environment, often in just seconds.
I Introduction
In this paper, we study the problem of optimal multirobot path planning on graphs (MPP), focusing on the design of complete algorithms and effective heuristics. In an MPP instance, the robots are uniquely labeled (i.e., distinguishable) and are confined to an arbitrary connected graph. The robots may move from a vertex to an adjacent one in one time step in the absence of collision, which may occur when two robots simultaneously move to the same vertex or along the same edge in different directions. A distinguishing feature of our MPP formulation is that we allow robots on fully occupied cycles to rotate synchronously. Such a formulation, more appropriate for multirobot applications, has not been widely studied (except, e.g., Standley (2010); Standley and Korf (2011)). Over the basic MPP formulation, we look at four commonly studied minimization objectives: the makespan (last arrival time), the maximum (singlerobot traveled) distance, the total arrival time, and the total distance. We note that these global objectives have direct relevance toward real world multirobot applications such as autonomous warehouse systems Wurman et al. (2008). For example, minimizing makespan is equivalent to minimizing the task completion time, whereas minimizing total distance is applicable to minimizing the fuel consumption of the entire fleet of robots.
In a related paper Yu and LaValle (2015), we show that these objectives are all distinct and NPhard to optimize, suggesting that efforts on solving optimal MPP should be directed at finding effective, nearoptimal algorithms. In this work, we make an attempt toward this goal and propose a novel yet general framework for solving MPP optimally. By examining space and time dimensions jointly, we observe a onetoone mapping between a solution for an MPP instance and that for a multicommodity flow problem derived from the MPP problem. Based on the equivalence relationship, we can then translate the MPP problem into an integer linear programming (ILP) model, which can be subsequently solved using an ILP solver. The generality of ILP allows the encoding of all four objectives to yield complete algorithms for optimizing these objectives. From here, we take a further step and introduce several heuristics to boost the algorithmic performance at a slight loss of solution optimality. Our method is especially effective in computing nearoptimal minimum makespan solutions, capable of computing optimal solutions for hundreds of robots, densely populated on the underlying graph, often in just seconds.
Related work. Multirobot path planning problems, in its many formulations, have been actively studied for decades Erdmann and LozanoPérez (1986); Zelinsky (1992); LaValle and Hutchinson (1998); Silver (2005); Kloder and Hutchinson (2006); Ryan (2008); Jansen and Sturtevant (2008); Surynek (2009); Luna and Bekris (2011); Standley and Korf (2011); van den Berg and Overmars (2005); van den Berg et al. (2009); Solovey and Halperin (2012); Yu and LaValle (2013a); Turpin et al. (2014); Solovey et al. (2015). As a universal subroutine, collisionfree path planning for multiple robots finds applications in tasks spanning assembly Halperin et al. (2000); Nnaji (1992), evacuation Rodriguez and Amato (2010), formation control Balch and Arkin (1998); Poduri and Sukhatme (2004); Shucker et al. (2007); Smith et al. (2009); Tanner et al. (2004), localization Fox et al. (2000), micro droplet manipulation Ding et al. (2001); Griffith and Akella (2005), object transportation Matarić et al. (1995); Rus et al. (1995), search and rescue Jennings et al. (1997), and so on. We refer the readers to Choset et al. (2005); Latombe (1991); LaValle (2006) and the references therein for a more comprehensive review on the general subject of multirobot path and motion planning.
The algorithmic study of graphbased multirobot path planning problems, the focus of this work, can be traced to as least 1879 Story (1879), in which Story makes the observation that the feasibility of the puzzle Loyd (1959) is decided by the parity of the game. The puzzle is a restricted MPP instance moving labeled game pieces on a grid, from some initial configuration to some goal configuration. The restriction is that only a single game piece near the only empty vertex may move to the empty vertex in a step–multiple mobile robots could move simultaneously. A generalization of the puzzle is introduced in Wilson (1974), extending the problem from game pieces on a grid to labeled pebbles on an vertex, connected graph. It is shown, together with an implied algorithm, that an instance is always feasible if the graph is nonbipartite. When the graph is bipartite (such as the puzzle), all pebble configurations are split into two groups of equal size such that any two configurations in the same group form a feasible instance. A further generalization is introduced in Kornhauser et al. (1984), allowing pebbles on a graph with vertices. For this problem, an algorithm is provided to solve an instance or decide that the instance is infeasible.
As computer games and multirobot systems gain popularity, concurrent movements are introduced and pebbles are replaced with robots (or agents). On the feasibility side, the MPP problem studied in this paper is shown to be solvable also in time Yu and Rus (2014). To distinguish the formulations, we denote the formulation that does not allow cyclic rotations of robots along fully occupied cycles as cyclefree MPP. Until recently, the majority of algorithmic study on MPP is on the cyclefree case. Since the problem is shown to be tractable Kornhauser et al. (1984), most algorithmic study of cyclefree MPP put some emphasis on optimality. Through the clever use of primitive operations, algorithms from Ryan (2008); Surynek (2009); Luna and Bekris (2011); Sajid et al. (2012); de Wilde et al. (2014) could quickly solve difficult problems with some form of completeness guarantees. These algorithms do not have optimality guarantees but the produced solutions are often of much better quality than the bound given by Kornhauser et al. (1984). For more discussion and references on suboptimal methods, see Sajid et al. (2012); de Wilde et al. (2014).
Pushing more toward the optimality side, most algorithmic results explore ways to limit the exponential search space growth induced by multiple robots. One of the first such algorithms, Local Repair A* (LRA*) Zelinsky (1992), plans robot paths simultaneously and performs local repairs when conflicts arise. Focusing on fixing the (locality) shortcomings of LRA*, Windowed Hierarchical Cooperative A* (WHCA*) proposes using a spacetime window to allow more choices for resolving local conflicts while simultaneously limiting the search space size Silver (2005). Iterative deepening technique is shown to be applicable to MPP problems in Sharon et al. (2012). A technique called subdimensional expansion is shown to perform well in complex environment Wagner and Choset (2011); the robot density is however relatively low ( cells per robot according to the paper). In Standley (2010); Standley and Korf (2011), instead of agnostic dissection of an instance, the natural idea of independence detection (ID) is explored to only consider multiple robots jointly (the source of exponential search space growth) as necessary. With operator decomposition (OD) that treats each legal move as an “operator”, the authors produced algorithms (ID,OD+ID, and related variants) that prove to be quite effective in computing total time or distanceoptimal solutions. We point out that ID and OD+ID have support for handling cycles (i.e., they apply to MPP instead of cyclefree MPP). Algorithms designed for minimizing makespan have also been attempted, e.g., Surynek (2012), but the solution quality degrades rapidly as the robotvertex ratio increases.
Many approaches have also been proposed for solving multirobot path planning problems in the continuous domain. A representative method called velocityobstacles Kant and Zucker (1986); van den Berg et al. (2008, 2011) explicitly examines velocitytime space for coordinating robot motions. In Griffith and Akella (2005), mixed integer programming (MIP) models are employed to encode the robot interactions. A method based on the spacetime perspective, similar to ours, is explored in I. Karamouzas (2012). In Peasgood et al. (2008), an A*based search is performed over a discrete roadmap abstracted from the continuous environment. In Solovey et al. (2014), discreteRRT (dRRT) is proposed for the efficient search of multirobot roadmaps. Algorithms for discrete MPP, cyclefree or not, have also helped solving continuous problems Solovey and Halperin (2012); Krontiris et al. (2012).
Contributions. We study the optimal MPP formulation allowing up to robots on a vertex connected graph, which we believe is better suited for multirobot applications. The formulation is not widely studied, perhaps due to the inherent difficulty in handling cyclic rotations of robots. Besides the novelty of the problem, this work brings several algorithmic contributions. First, based on the equivalence relationship between MPP and network flow, we establish a general and novel solution framework allowing the compact encoding of optimal MPP problems using ILP models. We show that the framework readily produces complete algorithms for minimizing the makespan (last arrival time), the maximum (singlerobot traveled) distance, the total arrival time, and the total distance, which are perhaps four most common global objectives for MPP. The resulting algorithms, in particular the one for computing minimum makespan, are highly effective in computing challenging problems instances with robotvertex ratio up to . Second, we introduce several principled heuristics, in particular a way split heuristic that divides an MPP instance over the time horizon, to give the exact algorithms a sizable performance boost at the expense of some solution optimality loss. With these heuristics, we are able to extend our algorithms to tackle problems with several hundred robots that are extremely densely populated, while at the same time maintain solution optimality. Last but not least, our successful exploit of ILP to attack optimal MPP shows that integer linear programming method is competitive with direct search methods in this problem domain, especially when the number of robots becomes large. This is surprising because ILP solvers are not designed specifically for MPP.
The rest of the paper is organized as follows. In Section II, we define the optimal MPP problems and provide background material on network flow. We then establish the equivalence relationship between MPP and multicommodity flow in Section III. We derive complete algorithms in Section IV based on the equivalence and continue to describe the performance boosting heuristics in Section V. We evaluate the algorithms in Section VI and conclude in Section VII. This paper is partly based on Yu and LaValle (2013a, b, c).^{1}^{1}1Yu and LaValle (2013c) is a preliminary poster presentation. In comparison to Yu and LaValle (2013a, b), besides demonstrating significantly improved computational performance due to the addition of the way split heuristic, we have substantially extended the generality of our ILPbased algorithmic framework, which now supports all common global time and distancebased objectives.
Ii Preliminaries
We now define MPP and the optimality objectives studied in this paper. Following the problem statements, we gave a brief review on network flow.
Iia MultiRobot Path Planning on Graphs
Let be a connected, undirected, simple graph, with being the vertex set and the edge set. Let be a set of robots. The robots move at discrete time steps (i.e., at ). At time step , each robot occupies a distinct vertex of . In general, at any time step , the robots assume a configuration that is an injective map from to . The start (initial) and goal configurations of the robots are denoted as and , respectively. Fig. 1(a) shows a possible configuration of robots on a grid graph. Fig. 1(b) shows a possible goal configuration, in which the robots are ordered based on rowmajor ordering ^{2}^{2}2In this paper, we generally use shaded discs to mark start locations of robots and discs without shading for goal locations..
(a)  (b) 
During a discrete time step, each robot may either remain stationary or move to an adjacent vertex. To formally describe a plan, let a scheduled path be a map , in which . A scheduled path is feasible if it satisfies the following properties: 1) . 2) For each , there exists a smallest such that . 3) For any , . 4) For any , or (if , robot stays at vertex between the time steps and ). We say that two paths are in collision if there exists such that (meet collision) or (headon collision). As an illustration, Fig. 2 shows the feasible and infeasible moves for two robots during a single time step ^{3}^{3}3We assume that the graph only allows “meet” or “headon” collisions. The assumption is mild. For example, a (arbitrary dimensional) grid with unit edge distance is such a graph for robots of with radii of no more than .. The multirobot path planning on graph (MPP) problem is defined as follows.
(a)  (b)  (c) 
Problem 1 (Multirobot Path Planning on Graphs)
Given a tuple , find a set of paths such that ’s are feasible paths for respective robots ’s and no two paths are in collision.
For example, Fig. 1(a) and Fig. 1(b) define an MPP problem on the grid. We call this particular problem the 9puzzle problem (a variant of the 15puzzle Ratner and Warmuth (1990)), which readily generalizes to puzzles.
Remark. With a few exceptions (e.g., Standley and Korf (2011)), most existing studies on discrete multirobot path planning problems require empty vertices as swap spaces. In these formulations, in a time step, a nonintersecting chain of robots may move simultaneously only if the head of the chain is moving into a previously unoccupied vertex. In contrast, our MPP formulation allows synchronized rotations of robots along fully occupied cycles (see, e.g., Fig. 1 and 3). This implies that even when the number of robots equals the number of vertices, robots can still move on disjoint cycles. We note that MPP can be solved in polynomial time with feasibility test taking only linear time Yu and Rus (2014).
IiB Optimal Formulations
Let be an arbitrary feasible solution to some fixed MPP instance. For a path , let denote the length of the path , which is increased by one each time when the robot passes an edge. A robot, following a path , may visit the same edge multiple times. Recall that denotes the arrival time of robot . In the study of optimal MPP formulations, we examine four common objectives with two focusing on time optimality and two focusing on distance optimality. Below, each objective is formally defined.
Objective 1 (Makespan)
Compute a path set that minimizes .
Objective 2 (Maximum Distance)
Compute a path set that minimizes .
Objective 3 (Total Arrival Time)
Compute a path set that minimizes .
Objective 4 (Total Distance)
Compute a path set that minimizes .
IiC Network Flow Review
A network consists of a directed graph with as the maps specifying capacities and costs over edges, respectively, and as the set of sources and sinks. Let , with denoting the set of source vertices, denoting the set of sink vertices, and . For a vertex , let (resp., ) denote the set of edges of going to (resp., leaving) . A feasible (static) flow on this network is a map that satisfies edge capacity constraints,
(1) 
the flow conservation constraints at non terminal vertices,
(2) 
and the flow conservation constraints at terminal vertices,
(3) 
The quantity is called the value of the flow . The classic (singlecommodity) maximum flow problem asks the following question: Given a network , what is the maximum that can be pushed through the network? The minimum cost maximum flow problem further requires the flow to have a minimum total cost among all maximum flows. That is, we want to find a flow among all maximum flows that also minimizes the quantity
(4) 
The network flow formulation described so far only considers a single commodity, corresponding to all robots being interchangeable. For general MPP formulations, the robots are distinct and must be treated as different commodities. Such problems can be captured with the Multicommodity flow or simply multiflow. Instead of having a single flow function , we have a flow function for each commodity . The constraints (1), (2), and 3 become
(5) 
(6) 
(7) 
Maximum flow and minimum cost flow problems may also be posed under a multiflow setup; we omit the details. Our review on network flow only touches aspects pertinent to the current work; for a thorough coverage on the subject of network flow, see Aronson (1989); Ahuja et al. (1993) and the references therein. Note that the multiflow model stated here is sometimes also referred to as integer multiflow because must have integer values.
Iii From MultiRobot Path Planning to Multiflow
A close algorithmic connection exists between optimal MPP and network flow problems. Maximum (single commodity) flow problems generally admit efficient (lowdegree polynomial time) algorithmic solutions Ahuja et al. (1993), whereas maximum multiflow is a wellknown NPhard problem, difficult to even approximate Andrews and Zhang (2005). Mirroring the disparity between single and multicommodity flows, in the domain of MPP problems, if there is a single group of interchangeable robots (here, for a group, it does not matter which robot goes to which goal as long as all goal locations assigned to the group are occupied by robots from the same group), then many optimal formulations admit polynomial time algorithms Yu and LaValle (2013a). However, as soon as a single group of robots splits into two or more groups, finding optimal paths for these robots become intractable Yu and LaValle (2015). The apparent similarity between optimal MPP and multiflow is perhaps best explained through a graphbased reduction from MPP problems to network flow problems. The reduction will also form the basis of our algorithmic solution.
To describe the reduction, we use as an example the undirected graph in Fig. 4(a), with start vertices and goal vertices . An instance of MPP is given by . We will reduce the problem to a network flow problem .
(a)  (b) 
The reduction proceeds by constructing a network that is a timeexpanded version of the graph , which then allows the explicit consideration of the interactions among the robots over space and time. To carry out this expansion, a time horizon must first be decided. For different optimality objectives, the expansion time horizon, some natural number , is generally different; for now we assume that is fixed.
To begin building the network, we create copies of ’s vertices, with indices , as shown in Fig. 5. For each vertex , denote these copies . For each edge and time steps , , the mergesplit gadget shown in Fig. 4(b) is added between and (arrows from the gadget are omitted from Fig. 5 since they are small). For the gadget, we assign unit capacity to all edges, unit cost to the horizontal middle edge, and zero cost to the other four edges. The mergesplit gadget ensures that two robots cannot travel in opposite directions on an edge in the same time step, which prevents headon collision between two robots. To finish the construction of Fig. 5, for each vertex , we add one edge between every two successive copies (i.e., we add the edges ). These correspond to the green and blue edges in Fig. 5. For all green edges, we assign them unit capacity and cost; for all blue edges, we assign them unit capacity and zero cost. The green edges allow robots to stay at a vertex during a time step whereas blue edges ensure that each vertex holds at most one robot, enforcing the meet collision constraint.
Fig. 5 (with the exception of edges and , which will become relevant shortly), the timeexpanded network, is the desired . For the set , we may simply let and . is now complete; we have reduced MPP to an integer multiflow problem on with each robot from as a single type of commodity.
Theorem 1
Let be an MPP instance. There is a bijection between its solution set (with a maximum number of time steps up to ) and the integer maximum multiflow solutions of flow value on the timeexpanded network constructed from with time steps.
Proof. Injectivity. Assume that ( is the number of robots) is a solution to an instance of MPP. For each and every time step , we mark the copy of and (recall that corresponds to a vertex of ) at time step in the timeexpanded graph . Connecting these vertices of sequentially (there is a unique way to do this) yields one unit of flow on (after connecting to appropriate source and sink vertices in , which is trivial). It is straightforward to see that if two paths are not in collision, then the corresponding flows on are vertex disjoint paths and therefore do not violate any flow constraint. Since any two paths in are not in collision, the corresponding set of flows is feasible and maximal on .
Surjectivity. Assume that is an integer maximum multiflow on the network with for all ’s. First we establish that any pair of flows are vertex disjoint. To see this, we note that (both are unit flows) cannot share the same source or sink vertices due to the unit capacity structure of enforced by the blue edges. If share some nonsink vertex at time step , both flows then must pass through the same blue edge (see Fig. 4(b)) with being either the head or tail vertex, which is not possible. Thus, are vertex disjoint on . We can readily convert each flow to a corresponding path (after deleting extra source vertex, sink vertices, vertices in the middle of the gadgets, and tail vertices of blue edges) with the guarantee that no will collide due to a meet collision. By construction of , the gadget we used ensures that a headon collision is also impossible. The set is then a solution to MPP. ∎
Remark. A multiflow problem can be reduced to an MPP problem (on a directed graph) as well. In particular, if all edges in the network have unit capacity and there is a single unit of each type of commodity, then a multiflow problem is a multirobot path planning problem, often known as the edge disjoint path problem Robertson and Seymour (1995).
Iv Complete, Integer Linear ProgrammingBased Algorithms for Optimal Mpp Problems
Because optimizing MPP solutions over Objectives 14 are computationally intractable, reducing MPP to multiflow problems does not make these optimal MPP problems any easier. However, with a network flow formulation (see Secion IIC), it becomes possible to establish integer linear programming (ILP) models for optimal MPP formulations. These ILP models can then be solved with powerful linear programming packages. In comparison to Abased algorithms augmented with heuristics which often target important but a limit set of problem structures, ILPbased algorithms proposed here are agnostic to specific problem structures. As such, ILPbased algorithms appear more capable at attacking a wider range of MPP problems and in particular difficult MPP instances in which the robotvertex ratio is high. In this section, we build ILP models for each of Objectives 14.
Iva Minimizing the Makespan
A minimum makespan solution to an MPP instance can be computed using a maximum multiflow formulation. Fixing a time span , let be the timeexpanded network for , a set of loopback edges are added to by connecting each pair of corresponding start and goal vertices in , from the goal to the start. We use ’s to denote edges of and let the loopback edges take the first indices, with being the edge connecting the goal vertex of to the start vertex of . For example, for the in Fig. 5, and are the loopback edges for and , respectively. The lookback edges have unit capacities and zero costs. Next. for each edge (including the lookback edges), create binary variables corresponding to the flow through that edge, one for each robot. The is, if and only if robot passes through in . The variables ’s must satisfy two edge capacity constraints and one flow conservation constraint,
(8) 
(9) 
The objective function is
(10) 
For each fixed , a solution to the above ILP model with an objective equaling means that a feasible solution to MPP is found. We are to find the minimal that yields such a feasible solution. To do this, start with being the maximum over all robots the shortest possible path length for each robot, ignoring all other robots. An ILP model for this is then built and tested for a feasible solution. If the model is not feasible, is increased and the procedure repeated. The first feasible is the optimal . Once a feasible model with the minimum expansion time horizon is found, the robots’ paths can be extracted based on the proof of Theorem 1. The algorithm is in fact complete: Since the problem is discrete, there is only a finite number of possible states. Therefore, for some sufficiently large , there must either be a feasible solution or we can pronounce that none can exist ( is big enough Yu and Rus (2014)). Denoting the resulting algorithm as MinMakespan we have shown the following.
Proposition 2
Algorithm MinMakespan is a complete algorithm for finding minimum makespan solutions for MPP.
IvB Minimizing the Maximum SingleRobot Traveled Distance
For minimizing the maximum distance traveled by any single robot, the network and variable creation remains the same as the minimum makespan setup; constraints (8) and (9) remain unchanged. Because we want to send all robots to their goal, a maximum flow may be forced through the additional constraint
(11) 
To encode the minmax objective function, we introduce an additional integer variable and add the constraint
(12) 
for all . For a fixed , the left side of (12) represents the distance traveled by robot . The objective function is then simply
(13) 
Denoting the algorithm as MinMaxDist, we have
Proposition 3
Algorithm MinMaxDist is a complete algorithm for finding solution to MPP that minimize the maximum distance traveled by a single robot.
Remark. Assuming that we have a minimum makespan solution, we can better bound the time horizon needed for MinMaxDist. Let be the minimum makespan computed by MinMakespan, setting is then sufficient for finding the . To see that this is true, we first note that because the minimum makespan solution requires robots to to synchronize their moves, which may force some robots to travel unnecessarily. Then, the robots cannot move a total distance of more than because no robot may travel more than edges and therefore cannot use more than time in total.
IvC Minimizing the Total Arrival Time
The ILP model for minimizing the total arrival time is more involved in the way the objective function is constructed. First, the network and all variables from the minimum makespan ILPmodel are inherited. We also inherit constraints (8), (9), and (11). To represent the objective function, for each time step and each , , we create a binary variable . Then, we give new indicies to some variables that are already created. Recall that for each edge (e.g., the four extra bold blue edges in Fig. 5), a variable is created. There are a total of such variables. Here, we give these variable a second index . That is, is the binary variable indicating whether edge , , is used by robot .
Given a network with a fixed , if constraints (8), (9), and (11) can be satisfied, then there is a feasible solution to the original MPP problem. In this case, for . We let . Then, each , is defined recursively over and as
(14) 
The constraint (14) effectively performs the logical and operation over and and stores the result in . In the end, the smallest for which is the time robot reaches its goal (and stops). Therefore, for each , is the number of time steps from the time arrives at its goal until time . is then the time spent by . To minimize the total arrival time, the objective function can be expressed as
(15) 
Denoting the resulting algorithm (with a sufficiently large ) as MinTotalTime, we have
Proposition 4
Algorithm MinTotalTime is a complete algorithm for finding minimum total arrival time solutions for MPP.
IvD Minimizing the Total Distance
From the ILP model for minimizing the maximum distance, we only need to change the objective function for computing a minimum total distance solution. We do not need the variable and simply update the objective function to
(16) 
Denoting the algorithm as MinTotalDist, we have
Proposition 5
Algorithm MinTotalDist is a complete algorithm for finding minimum total distance solutions for MPP.
Again, is sufficient for building a network that contains a minimum total distance solution, if one exists.
V Heuristics for Effective Computation of NearOptimal Solutions
In Section IV, in constructing the ILP models for Objectives 14, our goal is to show the universal applicability of the network flow model (e.g., Fig. 5) toward many different optimal MPP formulations. The resulting ILPbased algorithms are complete and always produce the optimal solution in principle. As will be shown in Section VI, such algorithms perform very well in handling relatively small but extremely challenging problems. Nevertheless, as the problem size grows (i.e., as the graph and the number of robots get larger), the computation time needed by ILP solvers grows rapidly. From a practical point of view, it may be far more desirable to quickly compute a good quality but suboptimal solution than waiting forever for the optimal solution. In this section, we introduce several heuristics to accomplish just that, with a particular focus on computing solutions with minimum makespan.
Va Building More Compact Models
To get the best performance out of a solver, it is beneficial to have a lean model (i.e., fewer columns and rows). So far, our focus has been to provide a general network flow based framework so that the ILP models can be easily built. When it comes to translating the models to an ILP solver, they can be further simplified. The heuristics discussed in this subsection aim at making the representation of the constraints (8) and (9) more compact in the resulting ILP models. As such, they apply to all the optimality objectives.
Better encoding of the collision constraints
Recall that in building the network flow model (e.g., Fig. 5), we used a mergesplit gadget (Fig. 4(b)) for enforcing the headon collision constraint and extra time steps (e.g., the blue edges in Fig. 5) for avoiding meet collisions. When we translate this into linear constraints, these structures can be simplified to yield the more compact structure illustrated in Fig. 6.
In the newer structure, each mergesplit gadget now has two edges instead of five. Also, the blue edges are removed. The updated gadget for an edge between time steps and is shown in Fig. 7 (note that due to the removal of the blue edges, vertices such as is no longer needed).
That is, instead of five, only two variables are needed for each robot. Denoting these binary variables as and for a robot , the headon collision constraint for a single gadget can be readily encoded as
(17) 
Then, to enforce the meet constraint, for example at a vertex , we simply require that at most one outgoing edge from may be used, i.e.,
(18) 
Overall, the newer ILP model is roughly half of the size of the original model.
Reachability analysis
In the timeexpanded graph, there are redundant binary (edge) variables that can never be true because some edges are never reachable. For example, in Fig. 6, at , the only outgoing edges that can possible be used are those originates from and . The rest can be safely removed. In general, for each robot , based on its reachability from its start vertex and to its goal vertex, a sizable number of binary variables ’s can be deleted.
VB DivideandConquer Over Time Domain
In evaluating the ILP modelbased algorithm for optimal makespan computation, we observe that the ILP solver running time appears to grow exponentially as the size of the model grows. This prevents the algorithm from performing well over instances with more than a few tens of robots. The observation, while hampering the effectiveness of the exact algorithm, turns out to offer a useful insight toward a highly effective heuristic. We find that, when the overall size of the ILP model is relatively small and the robotvertex ratio is not approaching , even when there are a large number of robots, the model usually does not present much challenge for ILP solvers. To apply the ILP modelbased method to more challenging problems (e.g., solving problems with hundreds of robots quickly), we simply limit the size of the individual ILP model fed to the solver. One way to achieve this is through divideandconquer over the time domain. We use a simple example (see Fig. 8) to illustrate the idea.
(a)  (b) 
In Fig. 8(a), we have a simple planning problem for two robots on a grid. To execute the heuristic, we first compute a shortest path for every pair of start and goal locations. In this case, we get the orange and green paths for robots and , respectively. Then, if we decide to split the problem into two smaller problems, for each of the paths, it is split into two (generally) equal length pieces and the middle node is set as the intermediate goal. In our example, we may do this for robot by setting the intermediate goal location at from the topleft corner (the brown disc labeled in Fig. 8(b)). For robot , because the middle location coincides with that of robot , we pick an alternative location that is not already occupied as the intermediate goal for robot , in this case from the topleft corner. The intermediate goals for the first instance will also serve as the start locations of the second instance. This yields two child instances with both requiring a time expansion with steps each, effectively making the individual ILP model roughly half the size of the original one, which requires a time expansion with steps. In general, we may divide a problem into arbitrarily many smaller instances in the time domain.
If a problem is divided in this manner to sub problems, we call the resulting heuristic a way split. Because the division is over time, there is no interaction between the individual, smaller instances. Once we obtain the solution for each child instance, the solutions can be glued together by concatenating the results. In practice, this simple heuristic dramatically improves algorithm performance without heavy negative impact on path optimality in terms of makespan; we observe a consistent speedup in computational experiments.
Remark. The way split heuristic, by design, is particularly suitable for the makespan objective. This is the case due to the additive nature of the makespan objective over the split subproblems. Besides makespan, the heuristic also applies to distance objectives (i.e., Objectives 2 and 4) quite well, as long as the time horizon required for finding distance optimal solution does not differ greatly from the time horizon required for minimum makespan solution. The heuristic does not directly apply to Objective 3 because total time is not additive over the split subproblems. As an example, suppose that a 2way split is carried out with each subproblem having a time horizon of . If a robot does not move in the solution to the first subproblem (i.e., ), it contributes to the total distance. However, if moves even a single step in the solution to the second subproblem (i.e., ), then will contribute at least to the total arrival time. Nevertheless, way split is still helpful in this case as we may use it to quickly compute an initial for performing the time expansion.
Vi Experimental Evaluation
In this section, we evaluate the performance of our optimal and nearoptimal MPP algorithms with an emphasis on MinMakespan. Our performance evaluation covers a broad spectrum of typical problem settings. For each setting, we push the limit on the robotvertex ratio–to as high as . To the best of our knowledge, the majority of the settings with high robotvertex ratio have never been attempted with much success prior to our study.
When applicable, we also compare our results with the stateoftheart found in the literature. In particular, we have examined OD+ID, ID, WHCA*, and Cobopt Surynek (2012), among others. OD+ID and ID support cycles whereas WHCA*appears to be designed for cyclefree MPP as it could not solve any puzzle. The problem definition for Cobopt suggests it solves MPP but it employs a cyclefree subroutine for finding feasible solutions. OD+ID, ID, and WHCA* are designed for optimizing total time and total distance optimal objectives, and do not naturally extend to makespan computation. However, the associated makespan produced by these algorithms are usually of good quality. Among these three, our experiments show that IDbased anytime algorithm is the most versatile due to its IDA*like incremental structure. On the other hand, OD+ID and WHCA* do not scale well when the robotvertex ratio goes beyond , depending on the particular problem setting.^{4}^{4}4Some of these algorithms were evaluated without requiring all robots reach their goals. For example, in the WHCA* work Silver (2005), if an instance with robots is solved for robots, the problem is counted as partially solved. We require each instance to be fully solved to be counted as a success. Cobopt is designed for makespan computation.
We implemented all algorithms (MinMakespan, MinMaxDist, MinTotalTime, and MinTotalDist) in the Java programming language. We take advantage of multicore CPUs when the way split heuristic is being used. Also, Gurobi Gurobi Optimization (2014), the ILP solver used in our implementation, can engage multiple cores automatically for hard problems. We ran all our tests on a MacBook Pro laptop computer (Intel Core i74850HQ, 16GB memory). We thank Trevor Standley for sharing the C code implementing OD+ID, ID, and WHCA*, among others. We modified (the original code supports only grid) and compiled the code as a 64bit windows executable under MSVC 2010 with all speed optimization flags turned on. The comparison to Cobopt uses the result provided in Surynek (2012), which covers only and grids.
Via Performance of MinMakespan and way Split Heuristic
We begin our experimental evaluation focusing on the MinMakespan algorithm and the way split heuristic. For this purpose, we use as the based graph a grid with varying number of vertices () removed to simulate obstacles. The connectivity of the graph is always maintained. We note that with vertices removed, the graph is already sparsely connected at some places (see e.g., Fig. 9),
making solving these problems optimally a challenging task. In presenting the computational results, each data point in the figures is an average over sequentially, randomly generated instances. For each obstacle percentage, we start from the lowest number of robots (usually or ) and allocate a maximum of seconds for each problem instance. If an instance takes more than seconds to produce a result, the instance is failed and we move to the next obstacle percentage. This also means that a data point is given only if each of the instances is completed within seconds. Over the same set of problem instances, the MinMakespan algorithm is executed in the exact manner (which produces optimal makespan solutions) and with the way split heuristic.
The exact makespan computation result is summarized in Fig. 10. For all obstacle settings, the MinMakespan algorithm computes optimal makespan solutions consistently for up to robots with an average computation time of no more than seconds. notably, for robots and obstacles up to , the MinMakespan algorithm is able to complete in about seconds in all cases. From the top plot of Fig. 10, we observe that for each fixed obstacle percentage, the computation time appears to grow exponentially with respect to the number of robots. The computational difficulty of a particular problem instance also heavily depends on the actual optimal makespan. For example, a problem instance in the case of and robots has a particularly long makespan (see the bottom plot of Fig. 10), resulting an unexpected jump of the computation time.
Whereas the exact MinMakespan algorithm is reasonably efficient, the way split heuristic brings a significant performance boost, allowing a much higher robotvertex density in general. In our tests, we are often able to double or triple the supported robotvertex density. For the grid, we evaluated the way split heuristic for up to . The way split performance is illustrated in Fig. 11.
In the figure, we measure optimality using a conservatively estimated optimality ratio. To obtain this, we divide the objective value returned by the optimizer over a conservative estimate (a lower bound). For makespan, this lower bound estimate is obtained by first computing the shortest path for each robot ignoring all other robots. The minimum makespan estimate is obtained by taking the maximum length over all these shortest paths. Clearly, the optimality ratio obtained this way is an underestimate.
We make three comments over Fig. 11. First, from the top plot, we observe that our method is highly effective in terms of computation time, capable of computing minimum makespan solutions for up to robots, which translates to a maximum robotvertex ratio of . The majority of the cases are solved in under seconds. Even when there are obstacles, we could solve the problem consistently for up to robots in about seconds. Second, we again observe an exponential relationship between computation time and the number of robots. Third, all computed solutions are very close to being optimal, with all but one case having an optimality ratio of below . The average minimum makespan for these instances is about .
The rest of the way split evaluation is presented in Fig. 12, in which the computation time and optimality ratio are shown side by side without the axis labels. We also omit the key of the plots, which is the same as those from Fig. 10, representing different obstacle percentages. In conjunction with Fig. 10 and Fig. 11, as increases, we observe a general trend of reduced computation time at the expense of some optimality loss. With split, MinMakespan can solve problems with robots, corresponding to a robotvertex ratio of .
For comparison, we ran IDbased anytime algorithm over the same set of instances with a second time limit (we also attempted OD+ID and WHCA*, which do not go past robots under the same setup). The result is plotted in Fig. 13. ID actually performs quite well for up to robots, which can be attributed to its A* root with minimum overhead as compared to our method. However, the performance of ID degrades faster–it does not scale well beyond robots in our tests. MinMakespan, with way split, readily outperforms ID when there are or more robots. Conceivably, it may be possible to combine ID and way split to make it run faster. However, adding way split to ID will inevitably make the overall makespan more suboptimal.
ViB Minimum Makespan Solution to puzzles
Next, we evaluate the efficiency of the algorithm MinMakespan for finding minimum makespan solutions to the puzzle for and . These problems have a robotvertex ratio of , making them highly constrained and extremely challenging. Note that an puzzle instance is always solvable for (see the Appendix); this means that all states are connected in the state (search) space. We ran Algorithm MinMakespan on 100 randomly generated puzzle instances for . For the 9puzzle, computation on all instances completed successfully with an average computation time of 0.46 seconds per instance. To compare the computational result, we implemented a (optimal) BFS algorithm. The BFS algorithm is heavily optimized: For example, cycles of the grid are precomputed and hard coded to save computation time. Since the state space of the 9puzzle is small, the BFS algorithm is capable of optimally solving the same set of 9puzzle instances with an average computation time of about 1.08 seconds per instance.
Once we move to the 16puzzle, the power of general ILP solvers becomes evident. MinMakespan solved all 100 randomly generated 16puzzle instances with an average computation time of 4.2 seconds. On the other hand, the BFS algorithm with a priority queue that worked for the 9puzzle ran out of memory after a few minutes. As our result shows that an optimal solution for the 16puzzle generally requires 6 time steps, it seems natural to also try bidirectional search, which cuts down the total number of states stored in memory. To complete such a search, one side of the bidirectional search generally must reach a depth of 3, which requires storing over states (the branching factor is over 1000), each taking 64 bits of memory. This translates into over 4GB of raw memory and over 8GB of working memory, which is more than the JavaVM can handle: A bidirectional search ran out of memory after about 10 minutes in general. We also experimented with C++ implementation using STL libraries, which yields similar result (i.e., ran out of memory before reaching a search depth of 3).
For the 25puzzle, without a good heuristic, bidirectional search cannot explore a tiny fraction of the fully connected state space with about states. On the other hand, MinMakespan again consistently solves the 25puzzle, with an average computational time of 391.6 seconds over 100 randomly created problems. Fig. 14 shows one of the solved instances with a 7step solution given in Fig. 15. Note that 7 steps is obviously the least
possible since it takes at least 7 steps to move robot 10 to its desired goal. We also briefly tested MinMakespan on the 36puzzle. While we had some success here, MinMakespan generally does not seem to solve a randomly generated instance of the 36puzzle within 24 hours, which has states and a branching factor of well over .
As comparison, WHCA* can not solve puzzle. OD+ID, and ID and can only solve puzzle consistently and cannot solve any puzzles in seconds.
ViC Minimum Makespan on , , and Grids
In this section, we evaluate MinMakespan with underlying graphs that are grids, grids, and grids with obstacles. In addition to further demonstrating the effectiveness of MinMakespan, this allows us to better compare our results. and grids are used as the environment for evaluation in Surynek (2012). grids with obstacles are used for evaluation in Silver (2005); Standley and Korf (2011).
For and grids, the instances are constructed using the same procedure stated in Section VIA. Again, each data point is an average over sequentially randomly created instances. Given the size of and grids, we limit to ; using way split can solve more instances but incur an average solution optimality between  optimal. Each instance is given a time limit of seconds. The outcome of these experiments is plotted in Fig. 16, along with the result from running ID. OD+ID and WHCA* cannot consistently solve the instances with robots in seconds. Makespan, instead of optimality ratio, is used in Fig. 16 for easy comparison with the results from Surynek (2012).
We observe that, over the grid, with way split heuristic, MinMakespan can solve problems with robots to almost true optimal solutions in just seconds. With way split, MinMakespan can further push to robots (robotvertex ratio of ) with solutions that are within optimal. ID can only handle up to robots. As reported in Surynek (2012), Cobopt generally takes more than half an hour to produce its final solution when there are or more robots.^{5}^{5}5We did not directly run Cobopt over our randomly created instances. However, the instances in Surynek (2012) are created in an identical manner. Therefore, given that similarly powered computers are used, the computation time and solution makespan are directly comparable between ours and those from Surynek (2012). The solution quality also degrades quickly as the number of robots increases. For example (Fig. 2 and Fig. 3 in Surynek (2012)), at robots, Cobopt takes over an hour to produce a solution with a makespan of over , whereas our way split heuristic yields a nearoptimal makespan of in just seconds.
Over the grid, MinMakespan is able to handle instances with (robotvertex ratio of ) robots with way split while at the same yielding solutions that are always less than optimal. When switched to way split, MinMakespan can consistently solve problems with up to robots to no more than optimal. In comparison, ID can solve instances with up to robots to relatively good quality. Taking on average an hour of computation, Cobopt can handle up to robots; the solution quality is quite poor. For example (Fig. 4 and Fig. 5 in Surynek (2012)), for about robots, the computed makespan by Cobopt is at whereas the optimal makespan is about . Across and grids, we generally observe a speedup of over times when MinMakespan (with the way split heuristic) is compared with Cobopt. At the same time, our method yields solutions with much smaller makespan.
A classical test scenario is the connected grid with vertices randomly removed.^{6}^{6}6Some work (e.g., Standley and Korf (2011)) also adopts an connected model. That is, each vertex is on the grid is connected to its neighborhood. This causes unit cost to be assigned to all edges, although a diagonal edge should have length times that of a nondiagonal edge. Since we are modeling robots in this work, we do not discuss the connected model here. However, we mention that our algorithms easily extend to connected model. Our tests show that we can in fact compute nearoptimal makespan for robots on grids assuming connectivity. For completeness, we also perform a brief evaluation of this setup. We randomly generate the instance as before, run the test, and plot the result in Fig. 17. From the figure, we observe a pattern consistent with experiments on the , and grids.
ViD Algorithm Performance over All Objectives
Last, we evaluate the performance of our algorithms at optimizing all Objectives 14. The result on MinMakespan is already presented in Section VIA, which we also use as the solution to optimizing Objective 2 (i.e., we simply use MinMakespan in place of MinMaxDist due to its superior performance). Note that no change to the plots are needed here because, on one hand, we can use a (near)optimal makespan solution as a solution to minimizing maximum singlerobot traveled distance. This is true because the makespan of a solution is always no less than the minimum maximum singlerobot distance. On the other hand, the lower bound estimate for minimum makespan is the same as that for minimum maximum singlerobot distance.
Before moving to MinTotalTime and MinTotalDist, we note that these algorithms possess some properties of an anytime algorithm, which is of practical importance. In solving these ILP models, the solver generally uses variations of the branchandbound algorithm Land and Doig (1960). For computing total time (and distance) optimal solutions, a branchandbound algorithm always try to find a feasible solution first and then iteratively improve over the feasible solution. This naturally leads to improving solution quality commonly observed in an anytime algorithm. The anytime property of MinTotalTime and MinTotalDist allows us to set a desired suboptimal threshold to reduce the computation time. Note that the same cannot be said for MinMakespan because a feasible solution here is the optimal solution.
Our next set of results focuses on the MinTotalTime algorithm (Fig. 18). The general setup is the same as that used in Section VIA. In particular, the same set of problem instances is used. In out experiment, we limit both time ( seconds) and required suboptimality threshold (automatically adjusted) to achieve a balanced performance. The lower bound estimate for computing optimality ratio is obtained by summing over the individual shortest path lengths. The actual optimal total time is about times the number of robots (i.e., for robots and for robots), regardless of the percentage of obstacles. We observe that the MinTotalTime algorithm is fairly effective, capable of computing optimal solutions for up to robots in the allocated time.
Similar outcome is also observed in the performance evaluation of the MinTotalDist algorithm with the way split heuristic (Fig. 19). The optimal total distance is again about times the number of robots. In comparison with the total time optimal case, due to the way split heuristic, the MinTotalDist algorithm is faster but produced solutions that are more suboptimal but still quite good.
For comparison, we run MinTotalTime, MinTotalDist, and ID (total time and total distance versions) on grids with  obstacles with a maximum time limit of seconds. The setting is slightly different from that used in obtaining Fig. 18 and 19; we do not set a suboptimal threshold here. The result is plotted in Fig. 20. In the case of total time (Objective 3), ID could solve more instances. For the instances that are solved, the achieved optimality is similar. For total distance (Objective 4), we observe similar outcome. Here, ID could produce one more data point; the achieved optimality by both methods are again comparable (and good). Overall, MinTotalTime, and MinTotalDist are competitive with ID.
Vii Conclusion
In this paper, we propose a general algorithmic framework for solving MPP problems optimally or nearoptimally. Through an equivalence relationship between MPP and network flow, we provide ILP modelbased algorithms for minimizing the makespan (last arrival time), the maximum (singlerobot traveled) distance, the total arrival time, and the total distance. In conjunction with additional heuristics, our algorithmic solutions are highly effective, capable of computing nearoptimal solutions for hundreds of robots in seconds in scenarios with very high robotvertex ratio. In pushing for high performance algorithms aim at solving MPP optimally or nearoptimally, our eventual goal is to apply these algorithms to multirobot path planning problems in continuous domains–our preliminary work toward this goal has begun to show promising results, producing algorithms that can compute solutions for around a hundred disc robots in 2D environments with holes.
References
 Yu and LaValle (2015) J. Yu and S. M. LaValle, “Optimal multirobot path planning on graphs: Structure and computational complexity,” 2015, reference to be updated.
 Standley (2010) T. Standley, “Finding optimal solutions to cooperative pathfinding problems,” in Proceedings AAAI National Conference on Artificial Intelligence, 2010, pp. 173–178.
 Standley and Korf (2011) T. Standley and R. Korf, “Complete algorithms for cooperative pathfinding problems,” in Proceedings International Joint Conference on Artificial Intelligence, 2011, pp. 668–673.
 Wurman et al. (2008) P. R. Wurman, R. D’Andrea, and M. Mountz, “Coordinating hundreds of cooperative, autonomous vehicles in warehouses,” AI Magazine, vol. 29, no. 1, pp. 9–19, 2008.
 Erdmann and LozanoPérez (1986) M. A. Erdmann and T. LozanoPérez, “On multiple moving objects,” in Proceedings IEEE International Conference on Robotics & Automation, 1986, pp. 1419–1424.
 Zelinsky (1992) A. Zelinsky, “A mobile robot exploration algorithm,” IEEE Transactions on Robotics & Automation, vol. 8, no. 6, pp. 707–717, 1992.
 LaValle and Hutchinson (1998) S. M. LaValle and S. A. Hutchinson, “Optimal motion planning for multiple robots having independent goals,” IEEE Transactions on Robotics & Automation, vol. 14, no. 6, pp. 912–925, Dec. 1998.
 Silver (2005) D. Silver, “Cooperative pathfinding,” in The 1st Conference on Artificial Intelligence and Interactive Digital Entertainment, 2005, pp. 23–28.
 Kloder and Hutchinson (2006) S. Kloder and S. Hutchinson, “Path planning for permutationinvariant multirobot formations,” IEEE Transactions on Robotics, vol. 22, no. 4, pp. 650–665, 2006.
 Ryan (2008) M. R. K. Ryan, “Exploiting subgraph structure in multirobot path planning,” Journal of Artificial Intelligence Research, vol. 31, pp. 497–542, 2008.
 Jansen and Sturtevant (2008) R. Jansen and N. Sturtevant, “A new approach to cooperative pathfinding,” in In International Conference on Autonomous Agents and Multiagent Systems, 2008, pp. 1401–1404.
 Surynek (2009) P. Surynek, “A novel approach to path planning for multiple robots in biconnected graphs,” in Proceedings IEEE International Conference on Robotics & Automation, 2009, pp. 3613–3619.
 Luna and Bekris (2011) R. Luna and K. E. Bekris, “Push and swap: Fast cooperative pathfinding with completeness guarantees,” in Proceedings International Joint Conference on Artificial Intelligence, 2011, pp. 294–300.
 van den Berg and Overmars (2005) J. van den Berg and M. Overmars, “Prioritized motion planning for multiple robots,” in Proceedings IEEE/RSJ International Conference on Intelligent Robots & Systems, 2005.
 van den Berg et al. (2009) J. van den Berg, J. Snoeyink, M. Lin, and D. Manocha, “Centralized path planning for multiple robots: Optimal decoupling into sequential plans,” in Robotics: Science and Systems, 2009.
 Solovey and Halperin (2012) K. Solovey and D. Halperin, “color multirobot motion planning,” in Proceedings Workshop on Algorithmic Foundations of Robotics, 2012.
 Yu and LaValle (2013a) J. Yu and S. M. LaValle, “Multiagent path planning and network flow,” in Algorithmic Foundations of Robotics X, Springer Tracts in Advanced Robotics (STAR). Springer Berlin/Heidelberg, 2013, vol. 86, pp. 157–173.
 Turpin et al. (2014) M. Turpin, K. Mohta, N. Michael, and V. Kumar, “CAPT: Concurrent assignment and planning of trajectories for multiple robots,” International Journal of Robotics Research, vol. 33, no. 1, pp. 98–112, 2014.
 Solovey et al. (2015) K. Solovey, J. Yu, O. Zamir, and D. Halperin, “Motion planning for unlabeled discs with optimality guarantees,” in Robotics: Science and Systems, 2015, to appear.
 Halperin et al. (2000) D. Halperin, J.C. Latombe, and R. Wilson, “A general framework for assembly planning: The motion space approach,” Algorithmica, vol. 26, no. 34, pp. 577–601, 2000.
 Nnaji (1992) B. Nnaji, Theory of Automatic Robot Assembly and Programming. Chapman & Hall, 1992.
 Rodriguez and Amato (2010) S. Rodriguez and N. M. Amato, “Behaviorbased evacuation planning,” in Proceedings IEEE International Conference on Robotics & Automation, 2010, pp. 350–355.
 Balch and Arkin (1998) T. Balch and R. C. Arkin, “Behaviorbased formation control for multirobot teams,” IEEE Transactions on Robotics & Automation, vol. 14, no. 6, pp. 926–939, 1998.
 Poduri and Sukhatme (2004) S. Poduri and G. S. Sukhatme, “Constrained coverage for mobile sensor networks,” in Proceedings IEEE International Conference on Robotics & Automation, 2004.
 Shucker et al. (2007) B. Shucker, T. Murphey, and J. K. Bennett, “Switching rules for decentralized control with simple control laws,” in American Control Conference, Jul 2007, pp. 1485–1492.
 Smith et al. (2009) B. Smith, M. Egerstedt, and A. Howard, “Automatic generation of persistent formations for multiagent networks under range constraints,” ACM/Springer Mobile Networks and Applications Journal, vol. 14, no. 3, pp. 322–335, Jun. 2009.
 Tanner et al. (2004) H. Tanner, G. Pappas, and V. Kumar, “Leadertoformation stability,” IEEE Transactions on Robotics & Automation, vol. 20, no. 3, pp. 443–455, Jun 2004.
 Fox et al. (2000) D. Fox, W. Burgard, H. Kruppa, and S. Thrun, “A probabilistic approach to collaborative multirobot localization,” Autonomous Robots, vol. 8, no. 3, pp. 325–344, Jun. 2000.
 Ding et al. (2001) J. Ding, K. Chakrabarty, and R. B. Fair, “Scheduling of microfluidic operations for reconfigurable twodimensional electrowetting arrays,” IEEE Transactions on Computeraided Design of Integrated Circuits and Systems, vol. 20, no. 12, pp. 1463–1468, 2001.
 Griffith and Akella (2005) E. J. Griffith and S. Akella, “Coordinating multiple droplets in planar array digital microfluidic systems,” International Journal of Robotics Research, vol. 24, no. 11, pp. 933–949, 2005.
 Matarić et al. (1995) M. J. Matarić, M. Nilsson, and K. T. Simsarian, “Cooperative multirobot box pushing,” in Proceedings IEEE/RSJ International Conference on Intelligent Robots & Systems, 1995, pp. 556–561.
 Rus et al. (1995) D. Rus, B. Donald, and J. Jennings, “Moving furniture with teams of autonomous robots,” in Proceedings IEEE/RSJ International Conference on Intelligent Robots & Systems, 1995, pp. 235–242.
 Jennings et al. (1997) J. S. Jennings, G. Whelan, and W. F. Evans, “Cooperative search and rescue with a team of mobile robots,” in Proceedings IEEE International Conference on Robotics & Automation, 1997.
 Choset et al. (2005) H. Choset, K. M. Lynch, S. Hutchinson, G. Kantor, W. Burgard, L. E. Kavraki, and S. Thrun, Principles of Robot Motion: Theory, Algorithms, and Implementations. Cambridge, MA: MIT Press, 2005.
 Latombe (1991) J.C. Latombe, Robot Motion Planning. Boston, MA: Kluwer, 1991.
 LaValle (2006) S. M. LaValle, Planning Algorithms. Cambridge, U.K.: Cambridge University Press, 2006, also available at http://planning.cs.uiuc.edu/.
 Story (1879) E. W. Story, “Note on the ‘15’ puzzle,” American Journal of Mathematics, vol. 2, pp. 399–404, 1879.
 Loyd (1959) S. Loyd, Mathematical Puzzles of Sam Loyd. New York: Dover, 1959.
 Wilson (1974) R. M. Wilson, “Graph puzzles, homotopy, and the alternating group,” Journal of Combinatorial Theory (B), vol. 16, pp. 86–96, 1974.
 Kornhauser et al. (1984) D. Kornhauser, G. Miller, and P. Spirakis, “Coordinating pebble motion on graphs, the diameter of permutation groups, and applications,” in Proceedings IEEE Symposium on Foundations of Computer Science, 1984, pp. 241–250.
 Yu and Rus (2014) J. Yu and D. Rus, “Pebble motion on graphs with rotations: Efficient feasibility tests and planning,” in Proceedings Workshop on Algorithmic Foundations of Robotics, 2014.
 Sajid et al. (2012) Q. Sajid, R. Luna, and K. E. Bekris, “Multiagent pathfinding with simultaneous execution of singleagent primitives.” in Fifth Annual Symposium on Combinatorial Search, 2012.
 de Wilde et al. (2014) B. de Wilde, A. W. ter Mors, and C. Witteveen, “Push and rotate: a complete multiagent pathfinding algorithm,” Journal of Artificial Intelligence Research, pp. 443–492, 2014.
 Sharon et al. (2012) G. Sharon, R. Stern, A. Felner, and N. Sturtevant, “ConflictBased Search for Optimal MultiAgent Path Finding,” in Proc of the TwentySixth AAAI Conference on Artificial Intelligence, 2012.
 Wagner and Choset (2011) G. Wagner and H. Choset, “M*: A complete multirobot path planning algorithm with performance bounds,” in Proceedings IEEE/RSJ International Conference on Intelligent Robots & Systems, 2011, pp. 3260–3267.
 Surynek (2012) P. Surynek, “Towards optimal cooperative path planning in hard setups through satisfiability solving,” in Proceedings 12th Pacific Rim International Conference on Artificial Intelligence, 2012.
 Kant and Zucker (1986) K. Kant and S. Zucker, “Towards efficient trajectory planning: The path velocity decomposition,” International Journal of Robotics Research, vol. 5, no. 3, pp. 72–89, 1986.
 van den Berg et al. (2008) J. van den Berg, M. C. Lin, and D. Manocha, “Reciprocal velocity obstacles for realtime multiagent navigation,” in Proceedings IEEE International Conference on Robotics & Automation, 2008, pp. 1928–1935.
 van den Berg et al. (2011) J. van den Berg, J. Snape, S. J. Guy, and D. Manocha, “Reciprocal collision avoidance with accelerationvelocity obstacles,” in Proceedings IEEE International Conference on Robotics & Automation, 2011, pp. 3475–3482.
 I. Karamouzas (2012) A. F. v. d. S. I. Karamouzas, R. Geraerts, “Spacetime group motion planning,” in Proceedings Workshop on Algorithmic Foundations of Robotics, 2012.
 Peasgood et al. (2008) M. Peasgood, C. Clark, and J. McPhee, “A complete and scalable strategy for coordinating multiple robots within roadmaps,” IEEE Transactions on Robotics, vol. 24, no. 2, pp. 283–292, 2008.
 Solovey et al. (2014) K. Solovey, O. Salzman, and D. Halperin, “Finding a needle in an exponential haystack: Discrete RRT for exploration of implicit roadmaps in multirobot motion planning,” in Proceedings Workshop on Algorithmic Foundations of Robotics, 2014.
 Krontiris et al. (2012) A. Krontiris, Q. Sajid, and K. Bekris, “Towards using discrete multiagent pathfinding to address continuous problems,” Workshop on Multiagent Pathfinding (in conjunction with AAAI), 2012.
 Yu and LaValle (2013b) J. Yu and S. M. LaValle, “Planning optimal paths for multiple robots on graphs,” in Proceedings IEEE International Conference on Robotics & Automation, 2013, pp. 3612–3617.
 Yu and LaValle (2013c) ——, “Fast, nearoptimal computation for multirobot path planning on graphs,” in Proceedings AAAI National Conference on Artificial Intelligence, 2013, late breaking papers.
 Ratner and Warmuth (1990) D. Ratner and M. Warmuth, “The puzzle and related relocation problems,” Journal of Symbolic Computation, vol. 10, pp. 111–137, 1990.
 Aronson (1989) J. E. Aronson, “A survey on dynamic network flows,” Annals of Operations Research, vol. 20, no. 1, pp. 1–66, 1989.
 Ahuja et al. (1993) R. K. Ahuja, T. L. Magnanti, and J. B. Orlin, Network Flows: Theory, Algorithms, and Applications. Upper Saddle River, NJ, USA: PrenticeHall, Inc., 1993.
 Andrews and Zhang (2005) M. Andrews and L. Zhang, “Hardness of the undirected edgedisjoint paths problem,” in Proceedings of the thirtyseventh annual ACM symposium on Theory of computing. ACM, 2005, pp. 276–283.
 Robertson and Seymour (1995) N. Robertson and P. D. Seymour, “Graph minors. XIII. the disjoint paths problem,” Journal of combinatorial theory, Series B, vol. 63, no. 1, pp. 65–110, 1995.
 Gurobi Optimization (2014) I. Gurobi Optimization, “Gurobi optimizer reference manual,” 2014. [Online]. Available: http://www.gurobi.com
 Land and Doig (1960) A. H. Land and A. G. Doig, “An automatic method of solving discrete programming problems,” Econometrica, vol. 28, no. 3, pp. 497–520, 1960.
Appendix
a Properties of the puzzle
The example problem from Fig. 1 easily extends to an grid; we call this class of problems the puzzle. Such problems are highly coupled: No robot can move without at least three other robots moving at the same time. At each step, all robots that move must move synchronously in the same direction (per cycle) on one or more disjoint cycles (see e.g., Fig. 3). To put into perspective the computational results on puzzles that follow, we make a characterization of the state structure of the puzzle for (the case of is trivial).
Proposition 6
All states of the 9puzzle are connected via legal moves.
Proof. We show that any state of a 9puzzle can be moved into the state shown in Fig. 1(b). From any state, robot 5 can be easily moved into the center of the grid. We are left to show that we can exchange two robots on the border without affecting other robots. This is possible due to the procedure illustrated in Fig. 21. ∎
Larger puzzles can be solved recursively: We may first solve the top and right side of the puzzle and then the left over smaller square puzzle. For a 16puzzle, Fig. 22 outlines the procedure, consisting of six main steps:

Move robots 1 and 2 to their respective goal locations, one robot at a time (first 1, then 2).

Move robots 3 and 4 (first 3, then 4) to the lower left corner (topmiddle figure in Fig. 22).

Move robots 3 and 4 to their goal location together via counterclockwise rotation along the cycle indicated in the topmiddle figure in Fig. 22.

Move robot 8 to its goal location.

Move robots 12 and then 16 to the lower left corner.

Rotate robots 12 and 16 to their goal locations.
It is straightforward to see that larger puzzles can be solved similarly. We have thus outlined the essential steps for proving Proposition 7 below; a more generic proof can be written using generators of permutation groups, which we omit here due to its length. Proposition 7 implies that, for , all instances of puzzles are solvable. The constructive proofs of Proposition 6 and 7 lead to recursive algorithms for solving any puzzle (clearly, the solution is not time/distance optimal in general).
Theorem 7
All states of an puzzle, are connected via legal moves.
Corollary 8
All instances of the puzzle, , are solvable.
By Proposition 7, since all states of a puzzle for are connected via legal moves, the state space of searching an puzzle equals factorial. For 16puzzle and 25puzzle, . Large state space is one of the three reasons that make finding a time optimal solution to the puzzle a difficult problem. The second difficulty comes from the large branching factor at each step. For a 9puzzle, there are 13 unique cycles, yielding a branching factor of 26 (clockwise and counterclockwise rotations). For the 16puzzle, the branching factor is around 500. This number balloons to over for the 25puzzle. This suggests that on typical commodity personal computer hardware (assuming a 1GHz processor), a baisc breadth first search algorithm will not be able to go beyond depth of 3 for the 16puzzle and depth 2 for the 25puzzle in reasonable amount of time. Moreover, enumerating these cycles is a nontrivial task. The third difficulty is the lack of obvious heuristics: Manhattan distances of robots to their respective goals prove to be a bad one. For example, given the initial configuration as that in Fig. 1(a), the first step in the optimal plan from Fig. 3 gets robots 1, 3, 4, 6, 8, 9 closer to their respective goals while moving robots 2, 7 farther. On the other hand, rotating counterclockwise along the outer cycle takes robots 1, 3, 4, 5, 6, 8, 9 closer and only moves robot 7 farther. However, if we instead take this latter first step, the optimal plan afterward will take 5 more steps.