We examine multi-robot path planning (MRPP, also known as multi-agent path finding or MAPF [45]) on two-dimensional grids, with potentially regularly distributed obstacles (see Fig. 1). The main objective of MRPP is to find a set of collision-free paths for routing many robots from a start configuration to a goal configuration. In practice, solution optimality is also of key importance; yet optimally solving MRPP is generally NP-hard [46,55], even in planar [52] and grid settings [10]. MRPP algorithms find many important largescale applications, including, e.g., in warehouse automation for general order fulfillment [51], grocery order fulfillment [35], and parcel sorting [50]. Other application scenarios include formation reconfiguration [38], agriculture [6], object transportation [40], swarm robotics [39,24], to list a few.
Motivated by applications including grocery fulfillment and parcel sorting, we focus on MRPP in which the underlying The authors are with the Department of Computer Science, Rutgers, the State University of New Jersey, Piscataway, NJ, USA. E-Mails: {teng.guo, jingjin.yu}@rutgers.edu. Fig. 1: (a) Real-world parcel sorting system (by JD.com) using many robots on a large grid-like environment with holes for dropping parcels; (b) A snapshot of a similar MRPP instance we can solve in polynomial-time with provable optimality guarantees. In practice, our algorithms scale to maps of size 450 × 300, supporting over 50K robots, and achieves 1.x-optimality (see, e.g., Fig. 11).
graph is an m 1 × m 2 grid, m 1 ≥ m 2 , with extremely high robot density. Whereas recent studies [53,10] have shown that such problems can be solved in polynomial time with O(1) optimality guarantees, the constant factor associated with the guarantee is generally prohibitively high (≫ 1) for these methods to be practical. In this research, we break this barrier by showing that, we can achieve (1 + δ)-makespan optimality for MRPP on large grids in polynomial-time in which δ ∈ (0, 0.5 + ε], ε → 0 as m 1 → ∞. Through the judicious application of a novel global object rearrangement method called Rubik Tables [48] together with many algorithmic techniques, and combined with careful analysis, we establish that in polynomial time:
• For m 1 m 2 robots, i.e., at maximum robot density, RTM (Rubik Table for MRPP) computes a solution for an arbitrary MRPP instance under a makespan of 7m 1 + 14m 2 ;
• For m1m2 3 robots and uniformly randomly distributed start/goal configurations, RTH (Rubik Tables with Highways) computes a solution with a makespan of m 1 + 2m 2 + o(m 1 ), with high probability. In contrast, such an instance has a minimum makespan of m 1 + m 2o(m 1 ) with high probability. This implies that, as m 1 → ∞, an optimality guarantee of m1+2m2 m1+m2 ∈ (1, 1.5] is achieved, with high probability;
• For m1m2 3 robots, for an arbitrary (i.e., not necessarily random) instance, a solution can be computed with a makespan of 3m 1 + 4m 2 + o(m 1 ) using RTH;
• For m1m2 2 robots, the same m1+2m2 m1+m2 optimality guarantee can be achieved with a slightly larger overhead using RTLM (Rubik Tables with Line Merge);
• The same m1+2m2 m1+m2 optimality guarantee may be achieved on grids with up to m1m2 9 regularly distributed obstacles together with 2m1m2 robots using RTH (e.g., Fig. 1(b)). Moreover, we have developed effective and principled heuristics to work together with RTH that further reduce the computed makespan by a large margin, i.e., for m1m2 3 robots, a makespan smaller than m 1 + 2m 2 can often be achieved. Demonstrated through extensive numerical evaluations, our methods are highly scalable, capable of solving instances with tens of thousands of robots in dense settings under two minutes. Simultaneously, the solution optimality approaches the 1-1.5 range as predicted theoretically. This level of scalability far exceeds what was possible. With the sub-1.5 optimality guarantee, our approach unveils a promising direction toward the development of practical, provably optimal multi-robot routing algorithms that runs in low polynomial time.
Related work. Literature on multi-robot path and motion planning [25,12] is expansive; here, we mainly focus on graph-theoretic (i.e., the state space is discrete) studies [56,45]. As such, in this paper, MRPP refers explicitly to graph-based multi-robot path planning. Whereas the feasibility question has long been positively answered for MRPP [26], the same cannot be said when it comes to securing optimal solutions, as computing time-or distance-optimal solutions are shown to be NP-hard in many settings, including for general graphs [16,46,55], planar graphs [52,2], and even regular grids [10], similar to the setting addressed in this study.
Nevertheless, given its high utility, especially in e-commerce applications [51,35,50] that are expected to grow significantly [9,1], many algorithmic solutions have been proposed for optimally solving MRPP. Among these, combinatorialsearch based solvers [27] have been demonstrated to be fairly effective. MRPP solvers may be classified as being optimal or suboptimal. Reduction-based optimal solvers solve the problem through reducing the MRPP problem to other problem, e.g., SAT [47], answer set programming [11], integer linear programming (ILP) [56]. Search-based optimal MRPP solvers include EPEA* [15], ICTS [42], CBS [43], M* [49], and many others. Due to the inherent intractability of optimal MRPP, optimal solvers usually exhibit limited scalability, leading to considerable interests in suboptimal solvers. Unbounded solvers like push-and-swap [31], pushand-rotate [8], windowed hierarchical cooperative A * [44], all return feasible solutions very quickly, but at the cost of solution quality. Balancing the running-time and optimality is one of the most attractive topics in the study of MRPP/MAPF. Some algorithms emphasize the scalability without sacrificing as much optimality, e.g., ECBS [3], DDM [22], EECBS [30], PIBT [36], PBS [34]. There are also learning-based solvers [7,41] that scales well in sparse environments. Effective orthogonal heuristics have also been proposed [18]. Recently, O(1)-approximate or constant factor time-optimal algorithms have been proposed, e.g. [53,10], that tackle highly dense instances. However, these algorithms only achieve lowpolynomial time guarantees at the expense of very large constant factors, rendering them theoretically interesting but impractical.
In contrast, with high probability, our methods run in low polynomial time with provable 1-1.5 asymptotic optimality. To our knowledge, this paper presents the first MRPP algorithms to simultaneously guarantee polynomial running time and 1.x solution optimality.
Organization. The rest of the paper is organized as follows. In Sec. II, we provide a formal definition of graph-based MRPP, and introduce the Rubik Table problem and the associated algorithm. RTM, a basic adaptation of the Rubik Table results for MRPP at maximum robot density which ensures a makespan upper bound of 7m 1 + 14m 2 , is described in Sec. III. An accompanying lower bound of m 1 +m 2 -o(m 1 ) for random MRPP instances is also established. In Sec. IV we introduce RTH for one third robot density achieving a makespan of m 1 + 2m 2 + o(m 1 ). Obstacle support is also discussed. In Sec. V, we show how one-half robot density may be supported with optimality guarantees similar to that of RTH. We thoroughly evaluate the performance of our methods in Sec. VI and conclude with Sec. VII. Given the amount of material included in the work, to provide a concentrated discussion, we refer the readers to [17] for proofs to theorems.
Graph-based multi-robot path planning (MRPP) seeks collision-free paths that efficiently route robots. Consider an undirected graph G(V, E) and n robots with start configuration S = {s 1 , . . . , s n } ⊆ V and goal configuration G = {g 1 , . . . , g n } ⊆ V . Each robot has start and goal vertices s i , g i . We define a path for robot i as a map P i : N → V where N is the set of non-negative integers. A feasible P i must be a sequence of vertices that connects s i and g i : 1)
With warehouse automation-like applications in mind, we work with G being 4-connected grids, aiming to minimize the makespan, i.e., max i {|P i |}. Unless stated otherwise, G is assumed to be an m 1 × m 2 grid with m 1 ≥ m 2 . Also, "randomness" in this paper always refers to uniform randomness. The version of MRPP we study is sometimes referred to as the one-shot MAPF problem [45]. We mention that our results also translate to guarantees on the life-long setting [45], which is briefly discussed in Sec. VII.
The Rubik Table problem (RTP) [48] formalizes the task of carrying out globally coordinated token swapping operations on lattices, with many interesting applications. The problem has many variations; in our study, we use the basic 2D form and the associated algorithms, which are summarized below. Problem 1 (Rubik Table Problem (RTP) [48]). Let M be an m 1 (row) × m 2 (column) table, m 1 ≥ m 2 , containing m 1 m 2 items, one in each table cell. The m 1 m 2 items are of m 2 colors with each color having a multiplicity of m 1 . In a shuffle operation, the items in a single column or a single row of M may be permuted in an arbitrary manner. Given an arbitrary configuration X I of the items, find a sequence of shuffles that take M from X I to the configuration where row i, 1 ≤ i ≤ m 1 , contains only items of color i. The problem may also be labeled, i.e., each item has a unique label in 1, . . . , m 1 m 2 .
A key result from [48], which we denote as the Rubik Table Algorithm (RTA), establishes that a colored RTP can be solved using m 2 column shuffles followed by m 1 row shuffles. Additional m 1 row shuffles then solve the labeled RTP.
Theorem 1 (Rubik Table Theorem [48]). An arbitrary Rubik Table problem on an m 1 ×m 2 table can be solved using m 1 + m 2 shuffles. The labeled Rubik Table problem can be solved using 2m 1 + m 2 shuffles.
We briefly illustrate how RTA works on an m 1 × m 2 table with m 1 = 4 and m 2 = 3 (in Fig. 2); we refer readers to [48] for more details. RTA operates in two phases. In the first phase, a bipartite graph B(T, R) is constructed based on the initial table configuration where the partite set T are the colors/types of items, and the set R are the rows of the table (Fig. 2(b)). An edge is added to B between t ∈ T and r ∈ R for every item of color t in row r. From B(T, R), a set of m 2 perfect matchings can be computed, as guaranteed by [21]. Each matching, containing m 1 edges, connects all of T to all of R, and dictates how a column should look like after the first phase. For example, the first set of matching in solid lines in Fig. 2(b) says that the first column should be ordered as yellow, cyan, red, and green as shown in Fig. 2(c). After all matchings are processed, we get an intermediate table, Fig. 2(c). Notice that each row of Fig. 2(a) can be shuffled to yield the corresponding row of Fig. 2(c); this is the key novelty of the RTA. After the first phase of m 1 row shuffles, the intermediate table (Fig. 2(c)) can then be rearranged with m 2 column shuffles to solve the colored RTP (Fig. 2(d)). Another m 1 row shuffles can then solve the labeled RTP (Fig. 2(e)). We note that it is also possible to perform labeled rearrangement using m 2 column shuffles followed by m 1 row shuffles and then followed by another m 2 column shuffles. The built-in global coordination capability of RTA naturally applies to solving makespan-optimal MRPP. Since RTA only requires three rounds of shuffles and each round involves either parallel row shuffles or parallel column shuffles, if each round of shuffles can be realized with makespan proportional to the size of the row/column, then a makespan upper bound of O(m 1 +m 2 ) can be guaranteed. This is in fact achievable even when all of G's vertices are occupied by robots, by recursively applying a labeled line shuffle algorithm [53], which can arbitrarily rearrange a line of m robots embedded in a grid using O(m) makespan.
Lemma 2 (Basic Simulated Labeled Line Shuffle [53]). For m labeled robots on a straight path of length m, embedded in a 2D grid, they may be arbitrarily ordered in O(m) steps. Moreover, multiple such reconfigurations can be performed on parallel paths within the grid.
The key operation is based on a localized, 3-step pair swapping routine, shown in Fig. 3. For more details on the line shuffle routine, see [53]. The basic simulated labeled line-shuffle algorithm, however, has a large constant factor. Borrowing ideas from parallel oddeven sort [4], we can greatly reduce the constant factor in Lemma 2. First, we need the following lemma. Proof of Lemma 3. Using integer programming [56], we exhaustively compute makespan-optimal solutions for arbitrary horizontal reconfiguration on 3×2 (8 possible cases) and 4×2 grids (16 possible cases), which confirms the claim.
As an example, it takes seven steps to horizontally "swap" all three pairs of robots on a 3 × 2 grid, as shown in Fig. 4. Lemma 4 (Faster Line Shuffle). For m robots on a straight path of length m, embedded in a 2D grid, they may be arbitrarily ordered in 7m steps. Moreover, multiple such reconfigurations can be performed simultaneously on parallel straight paths within the grid.
Proof. "Sorting" of m robots on a straight path of length m may be realized using parallel odd-even sort [4] in m -1 rounds, which only requires the ability to simulate potential pairwise "swaps" interleaving odd phases (swapping robots located at positions 2k + 1 and 2k + 2 on the path for some k) and even phases (swapping robots located at positions 2k + 2 and 2k + 3 on the path for some k). Here, it does not matter whether m is odd or even. To simulate these swaps, we can partition the grid embedding the path into 3 × 2 grids in two ways for the two phases, as illustrated in Fig. 5. A perfect partition requires that the second dimension of the grid, perpendicular to the straight path, be a multiple of 3. If this is not the case, some partitions at the bottom can use 4 × 2 grids. By Lemma 3, each odd-even sorting phase can be simulated using at most 7m steps. Clearly, shuffling on parallel paths is directly supported.
Combining RTA and fast line shuffle (Lemma 4) yields a polynomial time MRPP algorithm for fully occupied grids with a makepsan of 7m 1 + 14m 2 .
Theorem 5 (MRPP on Grids under Maximum Robot Density, Upper Bound). MRPP on an m 1 × m 2 grid, m 1 ≥ m 2 ≥ 3, with each grid vertex occupied by a robot, can be solved in polynomial time in a makespan of 7m 1 + 14m 2 .
We note that the case of m 2 = 2 can also be solved similarly except when m 1 = 2, with a slightly altered procedure since we can only use partitions of 2 × 3 grids. We omit the details for this minor case which readers can readily fill in.
The straightforward pseudo-code for RTM, the RTA based algorithm for MRPP on grids supporting the maximum possible robot density, is given in Alg. 1. The comments in the main RTM routine indicate the corresponding RTA phases. For MRPP on an m 1 × m 2 grid with row-column coordinates (x, y), we say robot i belongs to color 1 ≤ j ≤ m 1 if g i .y = j. Function Prepare() in the first phase finds intermediate states {τ i } for each robot through perfect matchings and routes them towards the intermediate states by (simulated) column shuffles. If the robot density is smaller than required, we may fill the table with "virtual" robots [23,53]. For each robot i we have τ i .y = s i .y. Function ColumnFitting() in the second phase routes the robots to their second intermediate states {µ i } through row shuffles where µ i .x = τ i .x and µ i .y = g i .y. In the last phase, function RowFitting() routes the robots to their final goal positions using additional column shuffles.
We now establish the optimality guarantee of RTM, assuming MRPP instances are randomly generated. For rearranging Proof. Without loss of generality, let the constant in Θ(m 1 m 2 ) be some c > 0, i.e., there are cm 1 m 2 robots. We examine the top left and bottom right corners of the m 1 × m 2 grid G. In particular, let G tl (resp., G br ) be the top left (resp., bottom right) αm 1 × αm 2 sub-grid of G, for some positive constant α ≪ 1. For u ∈ V (G tl ) and v ∈ V (G br ), assuming each grid edge has unit distance, then the Manhattan distance between u and v is at least (1 -2α)(m 1 + m 2 ). Now, the probability that some u ∈ V (G tl ) and v ∈ V (G br ) are the start and goal, respectively, for a single robot, is α 4 . For cm 1 m 2 robots, the probability that at least one robot's start and goal fall into G tl and G br , respectively, is
Because (1x) y < e -xy for 0 < x < 1 and y > 0 1 , p > 1-e -αfoot_1 cm1m2 . Therefore, for arbitrarily small α, we may choose m 1 such that p is arbitrarily close to 1. For example, we may let α = m -1 8 1 , which decays to zero as m 1 → ∞, then it holds that the makespan is
Comparing the upper bound established in Theorem 5 and the lower bound from Proposition 6 immediately yields Theorem 7 (Optimality Guarantee of RTM). For random MRPP instances on an m 1 × m 2 grid with Ω(m 1 m 2 ) robots, m 1 ≥ m 2 ≥ 3, as m 1 → ∞, RTM computes in polynomial time solutions that are 7(1 + m2 m1+m2 )-makespan optimal, with high probability.
We emphasize that RTM always runs in polynomial time and is not limited by any probabilistic guarantee; the high probability guarantee is only for solution optimality. The same is true for other algorithms' high probability guarantees proposed in this paper. We also note that high probability guarantees are stronger than and imply guarantees in expectation.
Though RTM runs in polynomial time and provides constant factor makespan optimality in expectation, the constant factor is still relatively large due to the extreme density. In practice, a robot density of around 1 3 (i.e., n = m1m2 3 ) is already very high. As it turns out, with n = cm 1 m 2 for some constant c > 0 and n ≤ m1m2 3 , which is assumed throughout this section, the constant factor can be dropped significantly by employing a "highway" heuristic to simulate the row/column shuffle operations.
Random MRPP instances. For the highway heuristics, we first work with random MRPP instances. Let us assume for the moment that m 1 and m 2 are multiples of three; we partition G into 3 × 3 cells (see, e.g., Fig. 1(b) and Fig. 6). We use Fig. 6, where Fig. 6(a) is a random start configuration and Fig. 6(f) is a random goal configuration, as an example to illustrate RTH-Rubik Table (for MRPP) with Highways, targeting robot density up to 1 3 . RTH involves two phases: anonymous reconfiguration and MRPP resolution with Rubik Table and highway heuristics.
In the anonymous reconfiguration phase, in which robots are treated as being indistinguishable or unlabeled, arbitrary start and goal configurations (under 1 3 robot density) are converted to intermediate configurations where each 3 × 3 cell contains no more than 3 robots. We call such configurations balanced configurations. With high probability, random MRPP instances are not far from being balanced. To establish this result (Proposition 9), we need the following.
Theorem 8 (Minimax Grid Matching [28]). Consider an m × m square containing m 2 points following the uniform distribution. Let ℓ be the minimum length such that there exists a perfect matching of the m 2 points to the grid points in the Theorem 8 applies to rectangles with the longer side being m as well (Theorem 3 in [28]).
Proposition 9. On an m 1 × m 2 grid, with high probability, a random configuration of n = m1m2 3 robots is of distance o(m 1 ) to a balanced configuration.
Proof. We prove for the case of m 1 = m 2 = 3m using the minimax grid matching theorem (Theorem 8); generalization to m 1 ≥ m 2 can be then seen to hold using the generalized version of Theorem 8 that applies to rectangles (Theorem 3 of [28], which in fact applies to arbitrarily simply connected region within a square region). Now let m 1 = m 2 = 3m. We may view a random configuration of m 2 robots on a 3m × 3m grid as randomly placing m 2 continuous points in an m×m square with scaling (by three in each dimension) and rounding. By Theorem 8, a random configuration of m 2 continuous points in an m × m square can be moved to the m 2 grid points at the center of the m 2 disjoint unit squares within the m × m square, where each point is moved by a distance no more than O(log 3 4 m), with high probability. Translating this back to a 3m × 3m gird, we have that m 2 randomly distributed robots on the grid can be moved so that each 3 × 3 cell contains exactly one robot and the maximum distance moved for any robot is no more than O(log 3 robots on an m 1 × m 1 gird can be moved so that each 3 × 3 cell contains exactly three robots and no robot needs to move more than a O(log three sets of reconfiguration paths will not cause an increase in the distance traveled by any robot (and will reduce it).
In the example, unlabeled reconfiguration corresponds to Fig. 6(a)→Fig. 6(b) and Fig. 6(f)→Fig. 6(e) (note that MRPP solutions are time-reversible). We simulated the process of anonymous reconfiguration for m 1 = m 2 = 300, i.e., on a 300 × 300 grids. For 1 3 robot density, the actual number of steps, averaged over 100 random instances, is less than 5. We call configurations like Fig. 6(b)-(e), which have all robots concentrated vertically or horizontally in the middle of the 3 × 3 cells, centered balanced configurations or simply centered configurations. Completing the first phase requires solving two unlabeled MRPP problems [54,32], easily doable in polynomial time.
In the second phase, RTA is applied with a highway heuristic to get us from Fig. 6(b) to Fig. 6(e), transforming between vertical centered configurations and horizontal centered configurations. To do so, RTA is applied (e.g., to Fig. 6(b) and (e)) to obtain two intermediate configurations (e.g., Fig. 6(c) and(d)). To go between these configurations, e.g., Fig. 6(b)→Fig. 6(c), we apply a heuristic by moving robots that need to be moved out of a 3 × 3 cell to the two sides of the middle columns of Fig. 6(b), depending on their target direction. If we do this consistently, after moving robots out of the middle columns, we can move all robots to their desired goal 3 × 3 cell without stopping nor collision. Once all robots are in the correct 3 × 3 cells, we can convert the balanced configuration to a centered configuration in at most 3 steps, which is necessary for carrying out the next simulated row/column shuffle. Adding things up, we can simulate a shuffle operation using no more than m + 5 steps where m = m 1 or m 2 . The efficient simulated shuffle leads to low makespan MRPP routing algorithms. It is clear that all operations take polynomial time; a precise running time is given at the end of this subsection.
Theorem 10 (Makespan Upper Bound for Random MRPP, ≤ 1 3 Density). For random MRPP instances on an m 1 × m 2 grid, where m 1 ≥ m 2 are multiples of three, for n ≤ m1m2 3
robots, an m 1 + 2m 2 + o(m 1 ) makespan solution can be computed in polynomial time, with high probability.
Proof. By Proposition 9, anonymous reconfiguration requires distance o(m 1 ) with high probability. By Theorem 1 from [53], this implies that a plan can be obtained for anonymous reconfiguration that requires o(m 1 ) makespan. For the second phase of MRPP resolution with Rubik Table and highway heuristics, by Theorem 1, we need to perform m 1 parallel row shuffles with row width of m 2 , followed by m 2 parallel column shuffles with column width of m 1 , followed by another m 1 parallel row shuffles with row width of m 2 . Simulating these shuffles require m 1 + 2m 2 + O(1) steps. All together, a makespan of m 1 + 2m 2 + o(m 1 ) is required, with very high probability.
Theorem 11 (Makespan Optimality for Random MRPP, ≤ 1 3 Density). For random MRPP instances on an m 1 × m 2 grid, where m 1 ≥ m 2 are multiples of three, for n = cm 1 m 2 robots with c ≤ 1 3 , as m 1 → ∞, a (1 + m2 m1+m2 ) makespan optimal solution can be computed in polynomial time, with high probability.
Since m 1 ≥ m 2 , 1 + m2 m1+m2 ∈ (1, 1.5]. In other words, in polynomial running time, RTH achieves (1 + δ) asymptotic makespan optimality for δ ∈ (0, 0.5], with high probability.
From the analysis so far, if m 1 and/or m 2 are not multiples of 3, it is clear that all results in this subsection continue to hold for robot density 1 3 -(m1 mod 3)(m2 mod 3)
, which is arbitrarily close to 1 3 for large m 1 . It is also clear that the same can be said for grids with certain patterns of regularly distributed obstacles (Fig. 1(b)), i.e., Corollary 12 (Random MRPP, 1 9 Obstacle and 2 9 Robot Density). For random MRPP instances on an m 1 × m 2 grid, where m 1 ≥ m 2 are multiples of three and there is an obstacle at coordinates (3k 1 + 2, 3k 2 + 2) for all applicable k 1 and k 2 , for n = cm 1 m 2 robots with c ≤ 2 9 , a solution can be computed in polynomial time that has makespan m 1 + 2m 2 + o(m 1 ) with high probability. As m 1 → ∞, the solution approaches 1 + m2 m1+m2 optimal, with high probability. Arbitrary MRPP instances. We now examine applying RTH to arbitrary MRPP instances under 1 3 robot density. If an MRPP instance is arbitrary, all that changes to RTH is the makespan it takes to complete the anonymous reconfiguration phase. On an m 1 × m 2 grid, by computing a matching, it is straightforward to show that it takes no more than m 1 + m 2 steps to complete the anonymous reconfiguration phase, starting from an arbitrary start configuration. Since two executions of anonymous reconfiguration are needed, this adds 2(m 1 + m 2 ) additional makespan. Therefore, we have Theorem 13 (Arbitrary MRPP, ≤ 1 3 Density). For arbitrary MRPP instances on an m 1 × m 2 grid, m 1 ≥ m 2 , for n ≤ m1m2 3 robots, a 3m 1 + 4m 2 + o(m 1 ) makespan solution can be computed in polynomial time. This implies that, for n = cm 1 m 2 ≤ m1m2 3 robots, a polynomial time can compute an asymptotic 3 + m2 m1+m2 makespan optimal solution, with high probability.
We now give the running time of RTM and RTH. 1 m 2 ) in deterministic time or O(m 1 m 2 log m 1 ) in expected time [14]. Anonymous MRPP may be tackled using the max-flow algorithm [13] in O(nm 1 m 2 T ) = O(nm 2 1 m 2 ) time, where T = O(m 1 + m 2 ) is the expansion time horizon of a time-expanded graph that allows a routing plan to complete.
Based on RTA, RTH has three simulated shuffle phases. Eventually, the makespan is dominated by the robot that needs the longest time to move, as a sum of moves for the robot in all three phases. As a result, the optimality of Rubik Table methods is determined by the first preparation phase. The matchings determine the intermediate states in all three phases. Finding arbitrary perfect matchings is fast but the process can be improved to reduce the overall makespan.
For improving matching, we propose two heuristics; the first is based on integer programming (IP). We create binary variables {x ri } where r represents the row number and i the robot. robot i is assigned to row r if x ri = 1. Define single robot cost as
We optimize the makespan lower bound of the first phase by letting λ = 0 or the third phase by letting λ = 1. The objective function and constraints are given by max r,i Eq. ( 1) is the summation of makespan lower bound of the first phase and the third phase. Note that the second phase cannot be improved through optimizing the matching. Eq. ( 2) requires that robot i be only present in one row. Eq. ( 3) specifies that each row should contain robots that have different goal columns. Eq. ( 4) specifies that each vertex (r, c) can only be assigned to one robot. The IP model represents a general assignment problem which is NP-hard in general. It has limited scalability but provides a way to evaluate how optimal the matching could be in the limit.
A second matching heuristic we developed is based on linear bottleneck assignment (LBA) [5], which takes polynomial time. LBA differs from the IP heuristic in that the bipartite graph is weighted. For the matching assigned to row r, the edge weight of the bipartite graph is computed greedily. If column c contains robots of color t, we add an edge (c, t) and its edge cost is
We choose λ = 0 to optimize the first phase. Optimizing the third phase (λ = 1) would give similar results. After constructing the weighted bipartite graph, an O( m 2.5 1 log m1 ) LBA algorithm [5] is applied to get a minimum bottleneck cost matching for row r. Then we remove the assigned robots and compute the next minimum bottleneck cost matching for next row. After getting all the matchings M r , we can further use LBA to assign M r to a different row r ′ to get a smaller makespan lower bound. The cost for assigning matching M r to row r ′ is defined as
The total time-complexity of using LBA heuristic for matching is O(
). We denote RTH with IP and LBA heuristics as RTH-IP and RTH-LBA, respectively. We mention that RTM, which uses the line swap motion primitive, can also benefit from these heuristics to re-assign the goals within each group. This can lower the bottleneck path length and improve the optimality.
The key design philosophy behind RTM and RTH is to effectively simulate row/column shuffles. With this in mind, we further explored the case of 1 2 robot density. Using a more sophisticated shuffle routine, 1 2 robot density can be supported while retaining most of the guarantees for the 1 3 density setting; obstacles are no longer supported.
To best handle 1 2 robot density, we employ a new shuffle routine called linear merge, based on merge sort, and denote the resulting algorithm as Rubik Table with Linear Merge heuristics or RTLM. The basic idea behind linear merge (shuffle) is straightforward: for m robots on a 2 × m grid, we iteratively sort the robots first on 2 × 2 grids, then 2 × 4 grids, and so on, much like how merge sort works. An illustration of the process on a 2 × 8 grid is shown in Fig. 7. We now show that linear merge is always feasible and has the desired properties.
Lemma 15 (Properties of Linear Merge). On a 2 × m grid, m robots, starting on the first row, can be arbitrarily ordered using m + o(m) steps. The motion plan can be computed in polynomial time.
Proof. We first show feasibility. The procedure takes ⌈log m⌉ phases; in a phase, let us denote a section of the 2 × m grid where robots are treated together as a block. For example, the left 2 × 4 grid in Fig. 7(b) is a block. It is clear that the first phase, involving up to two robots per block, is feasible (i.e., no collision). Assuming that phase k is feasible, we look at phase k + 1. We only need to show that the procedure is feasible on one block of length up to 2 k+1 . For such a block, the left half block of length up to 2 k is already fully sorted as desired, e.g., in increasing order from left to right. For the k + 1 phase, all robots in the left half block may only stay in place or move to the right. For these robots that stay, they must be all at the leftmost positions of the half block and will not block motions of any other robot. For the robots that do need to move to the right, their relative orders do not need to change, and therefore will not cause collisions among themselves. Because these robots that move in the left half block will move down on the grid by one edge, they will not interfere with any robot from left from the right half block. Because the same arguments hold for the right half block (except the direction change), the overall process of merging a block occurs without collision.
Next, we examine the makespan. For any single robot r, at phase k, suppose it belongs to block b and block b is to be merged with block b ′ . It is clear that the robot cannot move more than len(b ′ ) + 2 steps, where len(b ′ ) is the number of columns of b ′ and the 2 extra steps may be incurred because the robot needs to move down and then up the grid by one edge. This is because any move that r needs to do is to allow robots from b ′ to move toward b. Because there are no collisions in any phase, adding up all the phases, no robot moves more than m + 2(log m + 1) = m + o(m) steps.
Finally, it is clear that the merge sort-like linear merge shuffle primitive runs in O(m log m) time since it is a standard divide-and-conquer routine with log m phases.
With linear merge, the asymptotic properties of RTH for 1 3 robot density mostly carries over to RTLM.
Theorem 16 (Random MRPP, 1 2 Robot Density). For random MRPP instances on an m 1 × m 2 grid, where m 1 ≥ m 2 are multiples of two, for m1m2 3 ≤ n ≤ m1m2 2 robots, a solution can be computed in polynomial time that has makespan m 1 + 2m 2 + o(m 1 ) with high probability. As m 1 → ∞, the solution approaches an optimality of 1 + m2 m1+m2 ∈ (1, 1.5], with high probability.
In this section, we evaluate Rubik Table based algorithms and compare them with fast and near-optimal solvers, ECBS (w=1.5) [3] and DDM [22]. These two methods are, to our knowledge, two of the fastest near-optimal solvers for MRPP. We considered a state-of-the-art polynomial algorithm, pushand-swap [31], which gave fairly suboptimal results; the makespan optimality ratio is often above 100 for densities we examine. We also tested prioritized methods, e.g., [34,36,44], which faced significant difficulties in resolving deadlocks. Given the limited relevance and considering the amount of results we are presenting, we omit these methods in our comparison.
All experiments are performed on an Intel ® Core TM i7-9700 CPU at 3.0GHz. Each data point is an average over 20 runs on randomly generated instances, unless otherwise stated. A running time limit of 300 seconds is imposed over all instances. The optimality ratio is estimated as compared to conservatively estimated makespan lower bounds. The Rubik Table based algorithms are implemented in Python. The compared solvers are C++ based. As such, one can expect additional significant running time reductions from our algorithms implemented in C++. We choose Gurobi [20] as the mixed integer programming solver and ORtools [37] as the max-flow solver. The video of the simulations can be found at https://youtu.be/aphCjWFwfss.
We first evaluate the optimality achieved by RTM, RTLM, and RTH over randomly generated instances on their maximum designed robot density. That is, for RTM, the grid is fully occupied; for RTLM and RTH, the robot density is 1 2 and 1 3 , respectively. We test over three m 1 : m 2 ratios: 1 : 1, 3 : 2, and 5 : 1. The result is plotted in Fig. 8. Computation time is not listed (because we list the computation time later for RTH; the running times of RTM, RTLM, and RTH are similar). The optimality ratio is computed as the ratio between the solution makespan and the longest Manhattan distance between any pair of start and goal. Therefore, the estimate is an overestimate (i.e., the actual ratio may be lower/better). Fig. 8: Makespan optimality ratio for RTM, RTLM, and RTH for their maximum designed robot density, for different grid sizes and m1 : m2 ratios. We note that the largest problem has 90, 000 robots on a 300 × 300 grid.
We observe that RTM achieves 7-10.5+ makespan optimality ratio, which justifies the correctness of Theorem 7. Both and RTH achieve sub-2 optimality guarantee for most of the test cases, with result for RTH dropping below 1.5 on large grids. For all settings, as the grid size increases, there is a general trend of improvement of optimality across all methods/grid aspect ratios. This is due to two reasons: (1) the overhead in the shuffle operations becomes relatively smaller as grid size increases, and (2) with more robots, the makespan lower bound becomes closer to m 1 + m 2 . Lastly, as m 1 : m 2 ratio increases, the optimality ratio improves as predicted. For many test cases, the optimality ratio for the RTH for m 1 : m 2 = 5 setting is around 1.3.
We note that the performance of RTLM and RTH on optimality can be further improved using the heuristics described in Sec. IV-B. For the rest of the evaluations, we focus on RTH and its variants with additional heuristics.
1) Impact of grid size: For our first detailed comparative study of the performance of RTH, we set m 1 : m 2 = 3 : 2 and fix robot density at 1 3 . In Fig. 9, we compare the performance of ECBS [3], DDM [22], RTH, RTH-IP, and RTH-LBA, in terms of computation time and optimality ratio.
Despite the less efficient Python-based implementation, RTM, RTH, and RTH-LBA are faster than ECBS and DDM, due to their low polynomial running time guarantees. RTH and RTH-LBA can solve very large instances, e.g., on 450 × 300 grids with 45, 000 robots in about 100 seconds while neither ECBS nor DDM can. ECBS stopped working after m 2 = 30, though it shows better optimality on problems it can solve. DDM could handle up to m 2 = 90 but demonstrated poor optimality under the 1 3 density setting. The optimality ratio of RTH and RTH-LBA improves as the graph size increases, as predicted by our theoretical analysis. The optimality ratios of RTH and RTH-LBA reach as low as 1.49 and 1.26, respectively, agreeing with the ratio predicted by Theorem 11; here, the high probability asymptotic ratio is 1+ m2 m1+m2 = 1.4. RTH-LBA is able to do better than 1.4 because of the LBA heuristic. RTH-IP does slightly better on optimality in comparison to RTH and RTH-LBA, but its scalability is limited. On the other hand, RTH-LBA does nearly as well on optimality and remains competitive in terms of running time in comparison to RTH. As a consequence, we do not include further evaluation of RTH-IP.
For each method, the one standard deviation range is also shown in the figure. Because RTH and RTH-LBA are mostly deterministic and there are many robots, the change in optimality across different instances is small. We omit the inclusion of standard deviations from other plots as they are mostly similar in other tested settings.
We mention that we also evaluated ECBS with temporal splitting heuristics [19], which did not show significant difference in comparison to ECBS at the density we tried; so we did not include it here. We also evaluated push-and-swap [31], which runs fast but yields very poor optimality ratios (> 100 for many instances). We further evaluated prioritized methods, e.g., [36], which faced significant difficulties in resolving deadlocks. Given these, we did not include results from these methods in our comparatively study.
2) Impact of robot density: Next, we experiment the impact of different robot density on a 180 × 120 grid (Fig. 10). At density 1 3 , there are up to 7, 200 robots. For densities below 1 3 , we add "virtual robots" when the matching step is performed. ECBS does not appear because it cannot solve problems at this scale within 300 seconds.
RTH and RTH-LBA both return solutions around 10s on this graph for all instances. Both computation time and optimality ratio of DDM grow as the robot density increases, while the robot density has little impact on RTH and RTH-LBA. At lower density, DDM demonstrates better optimality. In environments with high densities, however, robots are highly coupled which causes more conflicts for structurallyagnostic approaches like DDM (and ECBS). RTH and RTH- LBA show improved optimality as the density increases, reaching 1.4, mostly due to the makespan of the instances getting larger.
3) Handling obstacles: RTH can also handle scattered obstacles and are especially suitable for cases where obstacles are regularly distributed. For instance, problems with underlying graphs like that in Fig. 1(b), where each 3×3 cell has a hole in the middle, can be natively solved without performance degradation. Such settings find real-world applications in parcel sorting facilities in large warehouses [50,29]. For this parcel sorting setup, we fix the robot density at 2 9 and test ECBS, DDM, RTH and RTH-LBA on graphs with varying sizes. The results are shown in Fig 11. Note that DDM can only apply when there is no narrow passage. So we added additional "borders" to the map to make it solvable for DDM. The results are similar as earlier ones; RTH and RTH-LBA run very fast and produce high-quality solutions, with conservatively estimated optimality ratio approaching 1.27. 4) Impact of grid aspect ratios: In this section, we fix m 1 m 2 = 90000 and vary the m 2 : m 1 ratio between 0 (nearly one dimensional) and 1 (square grids). We evaluated four algorithms, two of which are RTH and RTH-LBA. Now recall that RTP on an m 1 × m 2 table can also be solved using 2m 2 column shuffles and m 1 row shuffles. Adapting RTH and RTH-LBA with m 1 + 2m 2 shuffles gives the other two variants which we denote as RTH-LL and RTH-LBA-LL respectively, with "LL" suggesting two sets of longer shuffles are performed (each set of column shuffle work with columns of length m 1 ). The result is summarized in Fig. 12.
Interestingly but not surprisingly, the result clearly demonstrates the trade-offs between computation effort and solution optimality. RTH and RTH-LBA achieve better optimality ratio in comparison to RTH-LL and RTH-LBA-LL but require more computation time. Notably, the optimality ratios
Besides random start and goal settings, we also test RTH-LBA on many "special" instances; two are presented here (Fig. 13). For both settings, m 1 = m 2 . In the first, the "squares" setting, robots form concentric square rings and each robot and its goal are centrosymmetric. In the second, the "blocks" setting, the grid is divided into smaller square blocks (not necessarily 3 × 3) containing the same number of robots. robots from one block need to move to another random chosen block. RTH-LBA achieves optimality that is fairly close to 1.0 in the square setting and 1.7 in the block setting. The computation time is similar to that of Fig. 11; ECBS does well on optimality but scales poorly (only works on 30 × 30 grids). For some reason, DDM does very poorly on optimality and is not included.
In this study, we propose to apply Rubik Tables [48] to solving MRPP. A basic adaptation RTA, with a more efficient line shuffle routine, enables solving MRPP on grids at maximum robot density, in polynomial time, with previously unachievable optimality guarantee. Then, combining RTA, a highway heuristic, and additional matching heuristics, we obtain novel polynomial time algorithms that are provably asymptotically 1 + m2 m1+m2 makespan-optimal on m 1 × m 2 grids with up to 1 3 robot density, with high probability. Similar guarantees are also achieved with the presence of obstacles and at robot density up to one half. In practice, our methods can solve problems on graphs with over 10 5 number of vertices and 4.5 × 10 4 robots to 1.26 makespan-optimal (which can be better with larger m 1 : m 2 ratio). To our knowledge, no previous MRPP solvers provide dual guarantees on lowpolynomial running time and practical optimality.
Our study opens the door for many follow up research directions; we discuss a few here.
New line shuffle routines. Currently, RTLM and RTH only use two/three rows to perform a simulated row shuffle. Among other restrictions, this requires that the sub-grids used for performing simulated shuffle be well-connected (i.e. obstaclefree or the obstacles are regularly spaced so that there are at least two rows that are not blocked by static obstacles in each motion primitive to simulate the shuffle). Using more rows or even irregular rows in a simulated row shuffle, it is potentially possible to accommodate larger obstacles and/or support density higher than one half.
Better optimality at lower robot density. It is interesting to examine whether further optimality gains can be realized at lower robot density settings, e.g., 1 9 density or even lower, which are still highly practical. We hypothesize that this can be realized by somehow merging the different phases of RTA so that some unnecessary robot travel can be eliminated, after computing an initial plan.
Consideration of more realistic robot models. The current study assumes a unit-cost model in which a robot takes a unit amount of time to travel a unit distance and allow turning at every integer time step. In practice, robots will need to accelerate/decelerate and also need to make turns. Turning can be especially problematic and can cause significant increase in plan execution time, if the original plan is computed using the unit-cost model mentioned above. We note that RTH returns solutions where robots move in straight lines most of the time, which is advantageous in comparison to all existing MRPP algorithms, such as ECBS and DDM, which have large number of directional changes in their computed plans. it would be interesting to see whether the performance of RTA based MRPP algorithms will further improve as more realistic robot models are adapted.
Life-long MRPP settings. Currently, out RTA based MRPP solves are limited to a static setting whereas ecommerce applications of multi-robot motion planning often require solving life-long setting [33]. The metric for evaluating life-long MRPP is often the throughput, namely the number of goals reached per time step. We note that RTH also provides optimality guarantees for such settings, e.g., for the setting where m 1 = m 2 = m, we have Proposition 17 (Bound for Random Life-Long MRPP). The direct application of RTH to large-scale life-long MRPP on square grids yields an optimality ratio of 2 9 on throughput. Proof. We may solve life-long MRPP using RTH in batches. For each batch with n robots, RTH takes about 3m steps; the throughput is then T RT H = n 3m . As for the lower bound estimation of the throughput, the expected Manhattan distance in an m × m square, ignoring inter-robot collisions, is 2m 3 .
Therefore, the lower bound throughput for each batch is T lb = 3n 2m . The asymptotic optimality ratio is T RT H T lb = 2 9 . The 2 9 estimate is fairly conservative because RTH supports much higher robot densities not supported by known lifelong MRPP solvers. Therefore, it appears very promising to develop optimized Rubik Table inspired algorithms for solving life-long MRPP problems.
This is because log(1 -x) < -x for 0 < x < 1; multiplying both sides by a positive y and exponentiate with base e then yield the inequality.
m 1 ) steps, with high probability. We note that, because the robots are indistinguishable, overlaying