Multi-robot Path Planning and Scheduling via Model Predictive Optimal Transport (MPC-OT)

TL;DR

This paper introduces a novel multi-robot path planning method combining optimal transport and model predictive control to generate non-overlapping, cost-efficient trajectories in complex environments.

Contribution

It develops a new approach that integrates optimal transport with MPC for collision-free, optimal multi-robot navigation, addressing path overlaps and dynamic constraints.

Findings

Provides a computational complexity analysis for the proposed method.

Demonstrates the ability to generate non-overlapping paths efficiently.

Shows how to incorporate robot dynamics and re-planning into the transport framework.

Abstract

In this paper, we propose a novel methodology for path planning and scheduling for multi-robot navigation that is based on optimal transport theory and model predictive control. We consider a setup where $N$ robots are tasked to navigate to $M$ targets in a common space with obstacles. Mapping robots to targets first and then planning paths can result in overlapping paths that lead to deadlocks. We derive a strategy based on optimal transport that not only provides minimum cost paths from robots to targets but also guarantees non-overlapping trajectories. We achieve this by discretizing the space of interest into $K$ cells and by imposing a ${K\times K}$ cost structure that describes the cost of transitioning from one cell to another. Optimal transport then provides \textit{optimal and non-overlapping} cell transitions for the robots to reach the targets that can be readily deployed without any scheduling considerations. The proposed solution requires $\unicode{x1D4AA}(K^3\log K)$ computations in the worst-case and $\unicode{x1D4AA}(K^2\log K)$ for well-behaved problems. To further accommodate potentially overlapping trajectories (unavoidable in certain situations) as well as robot dynamics, we show that a temporal structure can be integrated into optimal transport with the help of \textit{replans} and \textit{model predictive control}.

Equations

$${\dot{\mathbf{x}}}_n(t)$$

$${\mathbf{y}}_n(t)$$

$$\underset{n,m}{\min }\sum _{l=1}^L{\int }_0^{T}c(q_{nm}(t))dt,$$

$$\mu (T^{-1}(B)):=\mu (\{a\in \mathcal{A}|T(a)\in B\})=\nu (B),\forall B\subseteq \mathcal{B},$$

$$\underset{T}{\inf }{\int }_{\mathcal{A}}c(a,T(a))d\mu (a)\text{mbox}s.t.T♯\mu =\nu .$$

$$\underset{\pi (a,b)\in \Pi (\mu ,\nu )}{\inf }{\int }_{\mathcal{A}\times \mathcal{B}}c(a,b)d\pi (a,b).$$

$$\mu (a)=\sum _{i=1}^n{\mu }_i{\delta }_{{\mathbf{p}}_i}(a),\nu (b)=\sum _{j=1}^m{\nu }_j{\delta }_{{\mathbf{r}}_j}(b),$$

$$\underset{\Pi =\{{\pi }_{ij}\}}{\inf }\sum _i\sum _j{\pi }_{ij}\cdot c({\mathbf{p}}_i,{\mathbf{r}}_j),$$

$$\sum _{j=1}^m{\pi }_{ij}={\mu }_i,\forall i,\sum _{i=1}^n{\pi }_{ij}={\nu }_j,\forall j,\sum _{i=1}^n{\mu }_i=\sum _{j=1}^m{\nu }_j.$$

$$\displaystyle\mu_{k}=\left\{\begin{array}[]{cc}1,&\mbox{cell~$k$ is not a target},\\ 0,&\mbox{cell~$k$ is a target}.\end{array}\right.$$

$$\displaystyle\nu_{k}=\left\{\begin{array}[]{cc}1,&\mbox{cell~$k$ is not a robot},\\ 0,&\mbox{cell~$k$ is a robot}.\end{array}\right.$$

$$\displaystyle C=\left[\begin{array}[]{cccccc}0&1&10^{2}&1&1&10^{3}\\ 1&0&1&1&1&1\\ 10^{2}&1&0&10^{3}&1&1\\ 1&1&10^{3}&0&1&10^{2}\\ 1&1&1&1&0&1\\ 10^{3}&1&1&10^{2}&1&0\end{array}\right],$$

$$\displaystyle\begin{array}[]{cccccc}\hskip-21.33955pt\bm{\mu}\backslash\bm{\nu}~~\>\>0&~1&~~\>1&\!~~~1&~~~1&~~1\end{array}$$

$$\displaystyle\left.\begin{array}[]{c}1\\ 1\\ 1\\ 1\\ 1\\ 0\end{array}\right|$$

$$\text{mbox}P1:\underset{\Pi =\{{\pi }_{ij}\}}{\inf }\sum _i$$

$$\Pi \cdot 1_K=\mathbf{\mu },1_K^{\top }\Pi ={\mathbf{\nu }}^{\top },$$

$$\underset{\Pi =\{{\pi }_{ij}\}}{\inf }\sum _i\sum _j{\pi }_{ij}\cdot c_{ij},i,j=\{1,\dots ,K\},$$

$$\Pi \cdot 1_K\le \mathbf{\mu },1_K^{\top }\Pi \le {\mathbf{\nu }}^{\top },1_K^{\top }\Pi 1_K=m.$$

$$\displaystyle\begin{array}[]{cccccc}\hskip-21.33955pt\bm{\mu}\backslash\bm{\nu}~~\>\>0&0&~~1&~~1&~~1&~~~1\end{array}$$

$$\displaystyle\left.\begin{array}[]{c}1\\ 1\\ 1\\ 1\\ 1\\ 0\end{array}\right|$$

$$\begin{array}[]{l}\min\limits_{\mathbf{u}_{n}}J_{n}:=\int_{t}^{t+T}(\mathbf{y}_{n}(t)-\mathbf{y}^{d}_{n}(r,t))^{\top}Q_{1}(\mathbf{y}_{n}(t)-\mathbf{y}^{d}_{n}(r,t))dt\\ \hskip 42.67912pt+\int_{t}^{t+T}\mathbf{u}_{n}(t)^{\top}Q_{2}\mathbf{u}_{n}(t)dt,\end{array}$$

$$\begin{array}[]{l}\dot{\mathbf{x}}_{n}(t)=f_{n}(\mathbf{x}_{n}(t))+g_{n}(\mathbf{x}_{n}(t))\cdot\mathbf{u}_{n}(t)\\ \mathbf{y}_{n}(t)=h_{n}(\mathbf{x}_{n}(t))\\ \mathbf{u}_{n}^{-}\leq\mathbf{u}_{n}(t)\leq\mathbf{u}_{n}^{+}\\ \mathbf{x}_{n}^{-}\leq\mathbf{x}_{n}(t)\leq\mathbf{x}_{n}^{+}\\ ||\mathbf{e}_{n}(t_{k+1})||_{P}\leq\alpha_{n}||\mathbf{e}_{n}(t_{k})||_{P},~~\alpha_{n}\in(0,1),\end{array}$$

$${\dot{p}}_{x,n}=v_n\cos {\theta }_n,{\dot{p}}_{y,n}=v_n\sin {\theta }_n,{\dot{\theta }}_n={\omega }_n,$$

Full text

1 INTRODUCTION

Robot-to-target assignment and navigation is an important problem in multi-robot swarms, automated warehouses, and multi-armed robotic manipulations Agaskar et al. [2025]. In these problems, $N$ robots are tasked to navigate to $M$ targets in a common space with obstacles. This problem is also known as MAPF (multi-agent path finding) Stern [2019], where robot-target pairs are given, and the task is to navigate the entire fleet of robots to their pre-assigned targets. A variant, anonymous MAPF Ali and Yakovlev [2023], is where any robot-target pairing is valid and we seek the pairing with minimal cost (travel, makespan) for the entire fleet. Anonymous MAPF (AMAPF) is tractable under certain assumptions, akin to solving a linear program, whereas general MAPF is significantly more challenging, requiring an integer linear program due to the terminal constraints.

In this paper, we study the AMAPF problem, which is highly relevant in scenarios where it is sufficient for any robot to reach any target. A plausible approach for AMAPF is to solve an assignment problem over an ${N\times M}$ cost matrix, whose elements are the minimum costs to travel between all robot-target pairs. This approach, although simple, has several disadvantages:

(i) needs pre-computed minimum travel costs between all robot-target pairs;

(ii) each min-cost path is computed independently and the procedure is susceptible to collisions, bottlenecks, and deadlocks;

(iii) the search space is restricted to exactly one min-cost path for every pair and thus ignores other paths of potentially the same cost.

In contrast, we propose an optimal transport (OT) framework, formulated over a $K$-cell discretization of the underlying space $\Omega$, that does not rely on pre-computed paths. In particular, we construct a ${K\times K}$ cost matrix $C$ of single cell transitions, that captures every possible path (sequence of transitions) between all robot-target pairs. The resulting transport $\Pi$ is optimized over $C$ and therefore over all possible paths between all robot-target pairs. The transport $\Pi$ must adhere to the marginal distribution constraints, which lead to non-overlapping paths. In situations where overlapping paths are unavoidable, we propose replans where the OT problem is solved again on a finer discretization of $\Omega$, equivalent to increasing travel capacity of the cells. Combining replans with MPC (model-predictive control), we show that each robot is able to track the OT trajectories, within a bounded error, and eventually converges exponentially to the final OT trajectory leading them to the desired targets.

While OT has a long history, it has gained renewed attention due to its applications in machine learning and data science Villani [2003], Peyré and Cuturi [2019], Figalli [2022]. Its flexibility and mathematical rigor have also motivated exploration in multi-agent path planning (PP) problems. Existing work includes Bandyopadhyay et al. [2014], de Badyn et al. [2018] that consider density control of networked multi-agent systems. Refs. Krishnan and Martínez [2018], Frederick et al. [2022] prescribe transport PDEs and gradient flows for mult-robot swarms, while Kachar and Gorodetsky [2022] integrates robot dynamics and LQR control costs in the ${N\times N}$ OT cost matrix. Related work also includes: Emerick and Bamieh [2023], on modeling swarms as continuum densities; Le et al. [2023], on entropic-regularized formulations; and more recently Liang et al. [2024], which utilizes diffusion models in continuous spaces. See also Stern [2019] for a comprehensive survey on MAPF, and Ma and Koenig [2016], Andreychuk et al. [2022, 2023], Fine et al. [2023] for some recent results on anonymous MAPF. In contrast, the proposed OT-based approach is formulated such that the entire space $\Omega$ is reconfigured from a starting distribution to the desired distribution, while naturally resulting in non-overlapping trajectories due to the marginal distribution constraints on the transport plan.

We further enhance our approach by incorporating real-time replanning, refining the space to account for cell travel capacity, and integrating temporal trajectory tracking using model predictive control (MPC) Rossiter [2004], Benosman [2016]. The receding horizon nature of MPC interleaved with the OT planner constitutes an attempt to integrate time constraints in OT, which is a notoriously hard problem Shi [2025]. Moreover, the proposed OT formulation opens up new directions for MAPF by e.g., adding robot dynamics and collision avoidance via appropriate control actions, while also offering the flexibility to shape paths with varying geometry and density through the cost matrix or additional constraints. Furthermore, it bridges MAPF with recent advances in OT, such as low-rank approximations, entropic smoothing, and gradient flow formulations, which can be employed to develop scalable solutions and other modern variations.

We now describe the rest of the paper. Section 2 provides the problem formulation and recaps OT. Section 3 describes our main contribution on using OT for joint matching and path planning. We consider replans and MPC-based OT in Sections 4 and 5, respectively. Finally, Section 6 provides numerical simulations and Section 7 concludes the paper.

2 PROBLEM FORMULATION

Consider $N$ robots and $M$ targets located in a region of interest ${\Omega\subset\mathbb{R}^{2}}$ with ${\kappa\subset\Omega}$ being the set of obstacles. The starting (centroid) location of robot $n$ is denoted by $\mathbf{p}_{n}$ and the $m$-th target location is given by $\mathbf{r}_{m}$, with $n\in[1,\ldots,N],m\in[1,\ldots,M]$. We assume that each robot $n$ has the following dynamics (${t\geq 0}$):

[TABLE]

where $\mathbf{x}_{n}(t)$ is the state of robot $n$, $\mathbf{u}_{n}(t)$ is the control input, and $\mathbf{y}_{n}(t)$ is its location, with ${\mathbf{y}_{n}(0)=\mathbf{p}_{n}}$. Let $q_{nm}(t)$ be a path that starts at robot $n$ and ends at target $m$. Let ${L=\min(N,M)}$, our goal is to design $L$ control inputs that map distinct robot-target pairs in a cost-optimal way:

[TABLE]

where $c(\cdot)$ denotes the cost to traverse a path and may include both the path length and the control cost. Since we are dealing with AMAPF, the minimization is over all robot-target pairs. Ideally, we would like to design trajectories that are non-overlapping and therefore each robot can charter its own course without any scheduling considerations.

This is a challenging problem overall as it requires minimizing the total cost (path length plus control inputs) over all distinct robot-target pairs. In this paper, we break this problem into three phases that are discussed in the following three sections:

(i) **Section 3–OT-based matching and path planning:**Find non-overlapping minimum cost paths that match $L$ distinct robots to $L$ distinct targets such that no other robot-target pair could result in a lower cost. We achieve this with the help of optimal transport without including the control cost in the first phase;

(ii) **Section 4–replans:**When overlapping paths are acceptable and can be implemented, e.g., through scheduling, we propose finer discretizations, effectively increasing travel capacity in certain lanes. To implement this, we decompose the paths obtained from OT into shorter segments and solve a series of OT problems for each segment;

(iii) **Section 5–Integrating MPC:**Finally, we combine controller design based on the robot kinematics (1)-(2) with the paths obtained from OT replans in order to add both physical and temporal constraints.

Before we proceed, we briefly recap optimal transport.

2.1 Optimal Transport

In 1781, the French mathematician Gaspard Monge formulated the problem of transporting mass from a source destination to a target destination that can be cast in terms of probability measures. Consider a source measure ${\mu:\mathcal{A}\rightarrow\mathbb{R}_{+}}$ and a target measure ${\nu:\mathcal{B}\rightarrow\mathbb{R}_{+}}$, where $\mu(A)$ tells us how much mass we have in the set ${A\subseteq\mathcal{A}}$, with ${\mu(\mathcal{A})=\nu(\mathcal{B})=1}$. Optimal transport seeks a map ${T:\mathcal{A}\rightarrow\mathcal{B}}$ such that

[TABLE]

where $\mu(T^{-1}(B))$ is the push-forward measure $T\sharp\mu$ of $\mu$ through $T$ Villani [2003], Peyré and Cuturi [2019], Figalli [2022]. Monge’s OT formulation is given by

[TABLE]

An issue with Monge’s formulation is that a feasible map $T$ may not exist because it does not allow mass to be split. For example, consider the source measure $\mu$ to be a delta function and the target measure $\nu$ to be some continuous measure, then $T$ can only push-forward the source delta to another delta and thus a feasible solution does not exist.

In 1939, Leonid Kantorovich provided an alternate formulation that does not consider a transport map but instead seeks a transport plan that allows mass to be split and go to different places. A transport plan $\Pi(A,B)$ describes how much mass goes from ${A\subseteq\mathcal{A}}$ to ${B\subseteq\mathcal{B}}$, which in fact is a measure on the product space ${\mathcal{A\times B}}$. We would like mass conservation that is given by ${\Pi(A,\mathcal{B})=\mu(A)},{\Pi(\mathcal{A},B)=\nu(B)}$, $\forall{A\subseteq\mathcal{A}},{B\subseteq\mathcal{B}}$. These constraints focus only on those measures, denoted by $\Pi(\mu,\nu)$, that have the correct marginals instead of all possible joint measures $\Pi(\mathcal{A,B})$. We now consider a cost $c(a,b)$ of going from ${a\in\mathcal{A}}$ to ${b\in\mathcal{B}}$ weighted further by the amount of mass moving from $a$ to $b$. Formally, Kantorovich’s transport problem is given by

[TABLE]

When both the source and target measures are discrete resources (as in robots being matched to targets), OT takes a very convenient form. Consider two discrete densities:

[TABLE]

where ${\delta_{\mathbf{p}_{i}}(a)=+\infty}$, when ${a=\mathbf{p}_{i}}$, and [math] otherwise. A transport plan is now a matrix ${\Pi=\{\pi_{ij}\}}$ with elements ${\pi_{ij}\geq 0}$ that tells us how much mass from $\mathbf{p}_{i}$ is moved to the location $\mathbf{r}_{j}$, such that mass conservation is not violated. Discrete OT is formally stated as

[TABLE]

The first constraint says that all mass at each $\mathbf{p}_{i}$ is appropriately transported; the second constraint is that the demand at each $\mathbf{r}_{j}$ is met; while the last constraint is mass conservation.

3 OT-BASED MATCHING and PATH PLANNING

In this section, we describe the main idea of this paper: to use optimal transport (OT) for joint matching and path planning. The discrete OT formulation described earlier has a direct relation to the robot-target matching problem. In its simplest form, the transport plan $\Pi$ can be chosen as an ${N\times M}$ matrix with elements in $\{0,1\}$ such that each row and each column sums to $1$; giving us the minimum cost robot-target matching. However, this implementation does not consider all available paths from the robots to the targets, does not explicitly account for obstacles, and therefore the resulting robot-target pairings may result in collisions, bottlenecks, and deadlocks. Our goal is not only to match robots to targets but to also compute optimal, non-overlapping paths that each robot could navigate. To achieve this, we cast a novel and unique OT formulation on a discretization of the underlying space $\Omega$.

Discretization of $\Omega$: Let $\mathcal{D}(\Omega)$ be a discretization of $\Omega$ into $K$ cells; the cell size/placement is such that no robot or target can occupy more than one cell. An example configuration of $\Omega$ and the grid-discretization is shown in Fig. 1; we choose a uniform, equally-spaced discretization for the sake of simplicity, which can be relaxed as we discuss later in Section 4.1.

Cost Matrix $C$: We now design a ${K\times K}$ cost matrix ${C=\{c_{ij}\}}$ that describes the transition cost from one cell to another such that moving to the adjacent cells is cheaper and incurs a cost of, e.g., $1$, while moving to any other distant cell has a higher cost and staying in the same cell has [math] cost. Similarly, the cost to move into an obstacle is chosen arbitrarily high. An example cost is shown in Fig. 1 (middle). With the help of this space discretization and cost structure, a particular novelty of our approach is that the robots are incentivized to travel to the targets by successively sliding into nearby cells instead of making a long-distance one-step match (jump) to the targets. The following OT formulation formally describes our approach.

OT Formulation: Optimal transport provides the optimal strategy to move mass from a source distribution to a target distribution. Let ${\bm{\mu}=\{\mu_{k}\}\in\mathbb{R}^{K}}$ be a $\{0,1\}$ vector of the source distribution such that

[TABLE]

Similarly, let ${\bm{\nu}=\{\nu_{k}\}\in\mathbb{R}^{K}}$ be a $\{0,1\}$ vector of the target distribution such that

[TABLE]

This unconventional encoding of the source and target distributions is then used in OT to move the mass at $\bm{\mu}$ to $\bm{\nu}$. Note that, by design, both $\bm{\mu}$ and $\bm{\nu}$ include a virtual one-unit mass on the entire free space where we have neither robots nor targets. Since the cost of staying in place is zero, this virtual mass is displaced only if it facilitates a lower-cost sequence of transitions to transport the non-virtual mass in the source (robots) to the non-virtual mass in the target distribution. An alternate viewpoint of this mass transport is that, instead of matching the robots to targets (and solving an ${N\times M}$ assignment), the proposed setup assigns mass to the entire space $\Omega$ such that mass transport is achieved by reconfiguring $\Omega$ from a source configuration to a desired target configuration in a cost-optimal way. The following example explains this concept.

*Example: *Consider a simple one robot and one target case in a region of interest $\Omega$ (Fig. 2). We discretize $\Omega$ into a $3\times 2$ grid, with ${K=6}$ cells. With this discretization, the robot occupies cell [math], while the target occupies cell $5$. An example cost matrix for this problem is

[TABLE]

where $c_{ij}$ is the cost of moving mass from cell $i$ to cell $j$. The source and target distribution vectors are ${\bm{\mu}^{\top}=[1~1~1~1~1~0]}$ and ${\bm{\nu}^{\top}=[0~1~1~1~1~1]}$, respectively. The assignment matrix ${\Pi=\{\pi_{ij}\}}$, shown below, is such that $\pi_{ij}$ represents how much of the mass $\mu_{i}$, at location $i$, is transported to meet the demand $\nu_{j}$, at location $j$. In other words, we would like the row sum of $\Pi$ to match $\bm{\mu}$ and the column sum of $\Pi$ to match $\bm{\nu}$ ensuring that the source mass is successfully transported to the target. Since the cost of keeping mass at its location is zero, i.e., ${c_{i=j}=0}$, choosing $\Pi$ to be identity could result in the smallest (zero) cost. However, since no mass can be brought to cell [math] (column-sum constraint), we must have ${\pi_{00}=0}$, i.e., mass at cell [math] cannot stay in place. Consequently, mass at cell [math] must incur a non-zero transition cost to ensure that it is transported out (row-sum constraint). In other words, the robot must be displaced to an adjacent cell $1,3$, or $4$ (each with a unit cost). Suppose, ${\pi_{04}=1}$, i.e., the robot slides into cell $4$ leading to ${\pi_{44}\neq 1}$ (column-sum constraint); mass at cell $4$ therefore moves out to satisfy the row-sum constraint, which can be achieved by ${\pi_{45}=1}$. The rest of the assignment is ${\pi_{11}=\pi_{22}=\pi_{33}=1}$, each with zero cost. The entire assignment is shown by underlined variables below and can be verified that it matches the constraints. The non-zero costs describe the path $0\rightarrow 4\rightarrow 5$ taken by the robot to the target.

[TABLE]

3.1 Balanced Case $(N=M)$

The above exposition can be cast in a precise mathematical framework. In particular, the OT problem is to find a transport plan $\Pi$, among all possible plans, that results into the minimum cost such that the source vector $\bm{\mu}$ in (8) is mapped to the target vector $\bm{\nu}$ in (11), i.e.,

[TABLE]

where $\mathbf{1}_{K}$ is a column vector with $K$ ones; the constraints ensure that the transport has the correct marginals. Clearly, multiple solutions may exist and the minimum cost path is not necessarily unique (both ${0\rightarrow 4\rightarrow 5}$ and ${0\rightarrow 1\rightarrow 5}$ in Fig. 2). Optimal transport provides one such min-cost transport as noted by the $\inf$ in P1. To describe the main results and guarantees of P1, we next provide a key concept.

Definition 1**.**

A practically feasible discretization $\overline{\mathcal{D}}(\Omega)$ of $\Omega$, with $K$ cells, is such that there exists a transport plan $\Pi$ on $\overline{\mathcal{D}}(\Omega)$ with the following properties:

1. (i)

$\mathbf{1}_{K}^{\top}(\Pi\odot C)\mathbf{1}_{K}$* is finite;* 2. (ii)

the row- and column-sum constraints are satisfied by $\Pi$; 3. (iii)

if $\pi_{ij}=1$, for any $i$ and $j$, then $c_{ij}\leq 1$.

Definition 1 formalizes all possible discretizations of $\Omega$, where the OT problem has a well-behaved solution, achievable by at-most one-step transitions (we assume for simplicity that one-step transitions that are physically realizable have a cost of $\leq 1$). In other words, $\Pi$ does not require a cell transition where a robot jumps more than one cell. Scenarios that violate Def. 1 will be considered in the next section. Key implications and guarantees regarding P1 are presented next.

Lemma 1**.**

Let ${N=M}$ robots and targets be situated in $\Omega$ with a practically feasible discretization $\overline{\mathcal{D}}(\Omega)$. The transport $\Pi^{*}$ obtained from the OT Problem in P1, described on $\overline{\mathcal{D}}(\Omega)$, has the following properties:

1. (i)

$\Pi^{*}$* has $\{0,1\}$ elements, i.e., the robots do not split on their paths to the targets;* 2. (ii)

$\Pi^{*}$* provides non-overlapping paths;* 3. (iii)

The total mass is conserved by the transport, i.e., all robots reach a distinct target; 4. (iv)

$\Pi^{*}$* is optimal under $\overline{\mathcal{D}}(\Omega)$ and the OT constraints;* 5. (v)

The complexity of P1 ranges from $\mathcal{O}(K^{3}\log K)$ to $\mathcal{O}(K^{2})$ for well-behaved problems.

Proof.

• It can be verified that P1 is a linear program (LP) with linear constraints and is well-defined over $\overline{\mathcal{D}}(\Omega)$, leading to the following.

(i) The marginal constraints result in a constraint matrix that is totally unimodular, ensuring that the underlying LP has integral solutions;

(ii) We can prove this by contradiction. Assume that two paths intersect at cell $j$, then masses from two different cells $i_{1}$ and $i_{2}$ arrive at cell $j$, i.e., ${\pi_{i_{1}j}=\pi_{i_{2}j}=1}$, which violates the column-sum constraint since ${\bm{\mu},\bm{\nu},\Pi\in\{0,1\}}$. The same argument applies to other situations.

(iii) We have ${\bm{\mu}^{\top}\mathbf{1}_{K}=\bm{\nu}^{\top}\mathbf{1}_{K}}$, since ${N=M}$, and the mass conservation constraint is satisfied, which in addition to (i)-(ii) ensures that all robots reach a distinct target;

(iv) P1is a feasible LP and therefore finds one of the optimal solutions.

(v) LPs can be solved using interior point methods, at $\mathcal{O}(K^{3}\log K)$, or with auction algorithms at $\mathcal{O}(K^{2})$ under favorable conditions.

∎∎

Assertion (i) represents a key theoretical contribution of this work, i.e., it establishes that the proposed OT formulation inherently preserves integrality, ensuring that mass is never split during the entire transport. This result leverages classical unimodularity results for integral solutions to LPs Schrijver [1998], but its application in this setting is both novel and essential. Assertions (ii)–(iv) then build on this foundation and ensure non-overlapping, optimal paths that cover all robots and targets, while assertion (v) follows from standard arguments on algorithmic and computational aspects Veinott and Dantzig [1968].

3.2 Unbalance Case $(N\neq M$)

The classical OT setup described above in P1 returns a transport plan $\Pi$ that preserves the total mass between the source and target, ensuring that each robot reaches a distinct target. In scenarios with unequal robots and targets, the formulation in P1 can be modified with the help of unbalanced optimal transport (UOT) as described next Chizat et al. [2018], Chapel et al. [2020]. Let ${m=\min(\|\bm{\mu}\|_{1},\|\bm{\nu}\|_{1})}$ be the least amount of mass that must be transported. The UOT formulation is given by

[TABLE]

This formulation is best explained with a simple example. Consider the same setup as in Fig. $2$ but with an additional robot at cell $1$. The cost matrix remains the same as before. The transport plan $\Pi$ with the UOT constraints in P2 is

[TABLE]

where the additional robot appears as ${\nu_{i1}=0},\forall i$. Following the column-sum constraint, the first two columns are now zeros implying that mass cannot be brought to either cell [math] or cell $1$. The UOT inequality constraints allow for a zero solution, which is avoided by the final constraint that $\Pi$ must have exactly ${m=\min(5,4)=4}$ ones. Among several possible plans $\Pi$, two examples are $\{{0\rightarrow 4\rightarrow 5},{2\rightarrow 2},{3\rightarrow 3}\}$ and $\{{1\rightarrow 5},{2\rightarrow 2},{3\rightarrow 3},{4\rightarrow 4}\}$, where the latter has a cost of $1$ and is returned by UOT.

4 REPLANS with OT

The OT and UOT formulations discussed up to now assume a practically feasible discretization $\overline{\mathcal{D}}(\Omega)$. As a consequence, the linear program implemented by OT is guaranteed to find a min-cost solution, with non-overlapping paths that comprise of at-most one cell transitions. We now consider scenarios where a preferred discretization $\underline{\mathcal{D}}(\Omega)$, not necessarily feasible, is given and the goal is to design robot trajectories within the confines of $\underline{\mathcal{D}}(\Omega)$. For a given $\underline{\mathcal{D}}(\Omega)$, however, implementing trajectories obtained from OT may not be ideal, see Fig. 3:

(i) $\underline{\mathcal{D}}(\Omega)$is practically feasible (Def. 1) but results in a longer trajectory;

(ii) non-overlapping trajectories with at most one-cell transitions do not exist.

These are consequences of the fact that the OT plans do not account for timing and scheduling possibilities. We emphasize that OT plans are still meaningful when ${\underline{\mathcal{D}}(\Omega)=\overline{\mathcal{D}}(\Omega)}$, as they can be implemented without any temporal considerations, as in Fig. 3 (left).

**Simple replans: **One approach to address potentially longer (but non-overlapping) OT trajectories is to sort them in an ascending-cost order and execute the shorter paths first (composed by unit cost transitions). The OT problem is then re-solved with a reduced set of robots and targets. Clearly, this approach handles both scenarios in Fig. 3, without bringing in path conflicts. For example, the shorter path in Fig. 3 (left) is executed first in $2$ steps and then OT is solved again to find the shortest path for the remaining robot. Another approach is to find a set of overlapping trajectories and bring timing and scheduling considerations to avoid conflicts. We describe this next.

4.1 Discretization (revisited)

Consider the scenario in Fig. 3 again with a given $\underline{\mathcal{D}}(\Omega)$. Let the robots and targets be point masses and assume a finer discretization $\underline{\underline{\mathcal{D}}}(\Omega)$ of $\underline{\mathcal{D}}(\Omega)$, such that $\underline{\underline{\mathcal{D}}}(\Omega)$ is practically feasible, i.e., OT returns lower cost on a higher resolution grid; see Fig. 4. In other words, the high resolution in $\underline{\underline{\mathcal{D}}}(\Omega)$ is equivalent to increasing the travel capacity of the corresponding cells. Clearly, from the perspective of OT in P1 or P2, the solutions in Fig. 4 are feasible as no constraints are violated because the trajectories traverse distinct spaces. These trajectories can be projected back in $\underline{\mathcal{D}}(\Omega)$, where a transport plan is originally sought, however, they cannot be readily implemented and require scheduling considerations with the help of replans as discussed in the next section. The following lemma characterizes the finer discretizations $\underline{\underline{\mathcal{D}}}(\Omega)$.

Lemma 2**.**

Consider ${N}$ robots and $M$ targets in $\Omega$, ${N,M<\infty}$, with a preferred discretization $\underline{\mathcal{D}}(\Omega)$. Assume that each robot has a path (in $\Omega$) to each target that does not cross through an obstacle, another robot, or another target. Then there exists a finer discretization $\underline{\underline{\mathcal{D}}}(\Omega)$ of $\underline{\mathcal{D}}(\Omega)$, which is practically feasible.

A formal proof is beyond the scope of this paper, but we provide a brief sketch of the argument. If each robot has at least one unblocked path to a target, feasible paths can be accommodated even in narrow regions by using finer discretizations to resolve conflicts. The condition stated in Lemma 2 avoids situations where a robot or a target permanently blocks another robot from reaching its goal, for example, by occupying the only traversable corridor. This requirement is related to the concept of valid infrastructures introduced in Cáp et al. [2015], which ensure that agents can maneuver without creating deadlocks. We note that non-uniform or adaptive discretizations, using finer resolution only in tight or congested areas, are also possible; see Fig. 4. While finer discretizations increase the likelihood of feasible, schedulable paths, they also incur greater computation, both in computing OT and in generating and implementing schedules. However, the OT complexity remains $\mathcal{O}(K^{3}\log K)$, provided the number of subdivisions does not scale with $K$.

In the next section, we use the set of paths generated by OT and design MPC strategies for each robot to track these paths while respecting the robot dynamics and control costs. The proper integration of OT replans with a receding horizon strategy allows us to extend the OT planning to the time-varying case, e.g., when the obstacles are dynamic.

5 MPC-OT: TRAJECTORY TRACKING

Up to now, we have neglected the robots’ dynamics in the OT formulations, however, in most real-life applications the path planned by OT will be executed by robots with physical constraints, given by (1). To include these dynamics in the overall multi-robot system control, we propose to close the loop on these trajectories via decentralized model predictive control (MPC) Rossiter [2004], Benosman [2016]. To this aim, for each robot ${n\in[1,\ldots,N]}$, let $\mathbf{y}_{n}^{d}(r,t)$ be the desired trajectory generated by the $r$th OT replan, obtained by breaking the OT paths, resulting for example from the aforementioned simple replans, into smaller segments. Let $[t_{k},t_{k+1}]$ denotes the time interval horizon of the MPC applied to the system, such that ${t_{k}=t}$, ${t_{k+1}=t+T_{c}}$, ${T_{c}<T}$, with $T>0$ denoting the MPC integration horizon. For each robot $n$, we solve the following contractive MPC problem:

[TABLE]

under the constraints

[TABLE]

where ${\mathbf{e}_{n}(t)=h_{n}(\mathbf{x}_{n}(t))-\mathbf{y}^{d}_{n}(r,t)}$ is the error process in tracking the trajectory, the matrices $Q_{1},Q_{2},P$ are arbitrary, positive definite, and $||\mathbf{e}_{n}(t)||_{P}=(\mathbf{e}_{n}(t)^{\top}P\mathbf{e}_{n}(t))^{1/2}$. Note that the first two constraints are the physical dynamics (1), the third and fourth constraints enforce that the state and control remains inside a bounded interval, while the last constraint ensures that the error process $\mathbf{e}_{n}(t)$ is contractive.

We will now prove boundedness of the tracking errors and exponential convergence to the target locations for all robots, when controlled by the MPC-OT presented in Algorithm 1. Before we proceed, we provide the assumptions next.

Assumption 1**.**

For any ${n\in[1,\ldots,N]}$, the system in (1) is controllable along its OT target output trajectories.

Assumption 2**.**

For any ${n\in[1,\ldots,N]}$ and any $r\in\{1,2,\ldots\}$, there exists a $\rho_{nr}\in\>]0,\infty[$ such that for all ${\mathbf{x}_{t_{0}}\in\mathbf{B}_{\rho_{nr}}}\stackrel{{\scriptstyle\triangle}}{{=}}\{{\mathbf{x}_{t_{0}}\in\mathbb{R}^{n}}|\;||\mathbf{e}_{n}(\mathbf{x}_{t_{0}})||_{P}\leq\rho_{nr}\}$, the MPC problem (12), (13) has a solution.

Assumption 3**.**

For any ${n\in[1,\ldots,N]}$ and any $r\in\{1,2,\ldots\}$, there exist a constant ${\beta_{nr}\in\>]0,\infty[}$, such that $||\mathbf{e}_{n}(t)||_{P}\leq\beta_{nr}||\mathbf{e}_{n}(t_{k})||_{P},\;\forall t\in[t_{k},\;t_{k+1}],\;k=0,1,\dots$.

Assumption 4**.**

For any ${n\in[1,\ldots,N]}$, and any $r\in\{1,2,\ldots\}$, there exists a $(t^{*},r^{*})$, s.t., $\forall{t>t^{*}},\;{r>r^{*}}$ the OT generates the same optimal trajectories $\mathbf{y}^{d^{*}}_{n}(r^{*},.)$, i.e., the OT replanning leads to the same target trajectories.

Assumptions 1 and 2 are classical in MPC, allowing us to ensure the feasibility of the MPC problem. The third assumption on the solutions’ boundedness is straightforward from the box constraints on the inputs and the continuity of the system’s model. The fourth assumption is needed to allow the system to settle on a final optimal trajectory to track. This assumption is needed to prove exponential stability of the error dynamics w.r.t. to the final optimal trajectory.

Theorem 1**.**

Under Assumptions 1–4, the MPC-OT Algorithm 1 applied to the multi-robot system (1), leads to bounded tracking errors $\mathbf{e}_{n}(t)$ and exponential convergence to the target locations $\mathbf{y}^{target}_{n}(T),\forall n=[1,\ldots,N]$.

The proof of Theorem 1 will be provided elsewhere. In short, it shows that the robots track, within a bounded error, the trajectories generated by the iterations of the OT planner, and eventually converge exponentially to the final OT trajectory leading them to the desired targets.

6 Simulations

In this section, we evaluate the effectiveness of the proposed MPC-OT algorithm with the help of numerical simulations. We consider a region $\Omega$ discretized into a $75\times 50$ grid, yielding ${K=3,\!750}$ cells; see Fig. 5. A team of ${N=20}$ unicycle robots is deployed in this environment, each governed by the following dynamics:

[TABLE]

where $p_{x,n}$ and $p_{y,n}$ denote the position of robot $n$, $\theta_{n}$ is its orientation, and $v_{n}$, $\omega_{n}$ are its forward and angular velocities. The state, control input, and spatial location of each robot are denoted by $\mathbf{x}_{n}=[p_{x,n},p_{y,n},\theta_{n}]$, $\mathbf{u}_{n}=[v_{n},\omega_{n}]$, and $\mathbf{y}_{n}=[p_{x,n},p_{y,n}]$, respectively. We set the number of targets equal to the number of robots (${M=N}$) and randomly generate obstacles with varying shapes, sizes, and locations. Robots and targets are randomly placed in the obstacle-free space and are marked as red triangles and blue circles, respectively. The optimal transport problem P1 is then solved to compute minimum-cost, non-overlapping trajectories that assign each robot to a corresponding target.

To obtain reference trajectories for model predictive control (MPC), we apply cubic time interpolation to the discrete OT-generated paths, assuming a fixed transition time of $1$ second between neighboring cells. These continuous reference trajectories are tracked using the MPC strategy described in Sec. 5. Each robot’s orientation in its starting cell is aligned with the direction of its motion. As shown in Fig. 5, all robots successfully follow their assigned OT paths and reach their respective targets. In dynamic environments, the set of occupied cells may change over time due to moving obstacles or new sensor information, necessitating online replans. To simulate such scenarios, we relocate one obstacle during MPC execution, that is effectively captured by the OT replans. As illustrated in Fig. 6, the robot’s reference trajectory is updated in response to the displaced obstacle, and the MPC controller successfully adapts to guide the robot to its target without collision.

7 Conclusions

In this paper, we address the problem of multi-robot target assignment and path planning using optimal transport (OT). Rather than directly matching individual robots to targets, we formulate the task as a mass reconfiguration problem, transforming the entire space $\Omega$ from an initial to a final distribution. We develop conditions under which cost-optimal, non-overlapping paths exist, and introduce modifications to handle cases where such paths cannot be realized. We incorporate temporal considerations through periodic replanning and receding horizon MPC controllers. The proposed methodologies offer a novel framework for multi-agent path finding (MAPF) combined with optimal transport, laying the groundwork for incorporating control-based constraints, modulating the spatial structure and distribution of paths, and leveraging modern advances in optimal transport.