Sinkhorn MPC: Model predictive optimal transport over
dynamical systems*

The problem of controlling a large number of agents has become a more and more important area in control theory with a view to applications in sensor networks, smart grids, intelligent transportation systems, and systems biology, to name a few. One of the most fundamental tasks in this problem is to stabilize a collection of agents to a desired distribution shape with minimum cost. This can be formulated as an optimal transport (OT) problem [1] between the empirical distribution based on the state of the agents and the target distribution over a dynamical system.

In [2, 3, 4, 5], infinitely many agents are represented as a probability density of the state of a single system, and the problem of steering the state from an initial density to a target one with minimum energy is considered. Although this approach can avoid the difficulty due to the large scale of the collective dynamics, in this framework, the agents must have the identical dynamics, and they are considered to be indistinguishable. In addition, even for linear systems, the density control requires us to solve a nonlinear partial differential equation such as the Monge-Ampère equation or the Hamilton-Jacobi-Bellman equation. Furthermore, it is difficult to incorporate state constraints due to practical requirements, such as safety, into the density control.

With this in mind, we rather deal with the collective dynamics directly without taking the number of agents to infinity. This straightforward approach is investigated in multi-agent assignment problems; see e.g., [6, 7] and references therein. In the literature, homogeneous agents following single integrator dynamics and easily computable assignment costs, e.g., distance-based cost, are considered in general. On the other hand, for more general dynamics, the main challenge of the assignment problem is the large computation time for obtaining the minimum cost of stabilizing each agent to each target state (i.e., point-to-point optimal control (OC)) and optimizing the destination of each agent. This is especially problematic for OT in dynamical environments where control inputs need to be determined immediately for given initial and target distributions.

For the point-to-point OC, the concept of model predictive control (MPC) solving a finite horizon OC problem instead of an infinite horizon OC problem is effective in reducing the computational cost while incorporating constraints [8]. On the other hand, in [9], several favorable computational properties of an entropy-regularized version of OT are highlighted. In particular, entropy-regularized OT can be solved efficiently via the Sinkhorn algorithm. Based on these ideas, we propose a dynamical transport algorithm combining MPC and the Sinkhorn algorithm, which we call Sinkhorn MPC. Consequently, the computational effort can be reduced substantially. Moreover, for linear systems with an energy cost, we reveal the fundamental properties of Sinkhorn MPC such as ultimate boundedness and asymptotic stability.

Organization: The remainder of this paper is organized as follows. In Section II, we introduce OT between discrete distributions. In Section III, we provide the problem formulation. In Section IV, we describe the idea of Sinkhorn MPC. In Section V, numerical examples illustrate the utility of the proposed method. In Section VI, for linear systems with an energy cost, we investigate the fundamental properties of Sinkhorn MPC. Some concluding remarks are given in Section VII.

Notation: Let ${\mathbb{R}}$ denote the set of real numbers. The set of all positive (resp. nonnegative) vectors in ${\mathbb{R}}^{n}$ is denoted by ${\mathbb{R}}_{>0}^{n}$ (resp. ${\mathbb{R}}_{\geq 0}^{n}$). We use similar notations for the set of all real matrices ${\mathbb{R}}^{m\times n}$ and integers ${\mathbb{Z}}$, respectively. The set of integers $\{1,\ldots,N\}$ is denoted by $[\![N]\!]$. The Euclidean norm is denoted by $\|\cdot\|$. For a positive semidefinite matrix $A$, denote $\|x\|_{A}:=(x^{\top}Ax)^{1/2}$. The identity matrix of size $n$ is denoted by $I_{n}$. The matrix norm induced by the Euclidean norm is denoted by $\|\cdot\|_{2}$. For vectors $x_{1},\ldots,x_{m}\in{\mathbb{R}}^{n}$, a collective vector $[x_{1}^{\top}\ \cdots\ x_{m}^{\top}]^{\top}\in{\mathbb{R}}^{nm}$ is denoted by $[x_{1};\ \cdots\ ;x_{m}]$. For $A=[a_{1}\ \cdots\ a_{n}]\in{\mathbb{R}}^{m\times n}$, we write ${\rm vec}(A):={\color[rgb]{0,0,0}{[a_{1};\ \cdots\ ;a_{n}]}}$. For $\alpha=[\alpha_{1}\ \cdots\ \alpha_{N}]^{\top}\in{\mathbb{R}}^{N}$, the diagonal matrix with diagonal entries $\{\alpha_{i}\}_{i=1}^{N}$ is denoted by $\alpha^{\boxbslash}$. The block diagonal matrix with diagonal entries $\{A_{i}\}_{i=1}^{N},A_{i}\in{\mathbb{R}}^{m\times n}$ is denoted by $\{A_{i}\}_{i}^{\boxbslash}$. Especially when $A_{i}=A,\forall i$, $\{A_{i}\}_{i}^{\boxbslash}$ is also denoted by $A^{\boxbslash,N}$. Let $(M,d)$ be a metric space. The open ball of radius $r>0$ centered at $x\in M$ is denoted by $B_{r}(x):=\{y\in M:d(x,y)<r\}$. The element-wise division of $a,b\in{\mathbb{R}}_{>0}^{n}$ is denoted by $a\oslash b:=[a_{1}/b_{1}\ \cdots\ a_{n}/b_{n}]^{\top}$. The $N$-dimensional vector of ones is denoted by $1_{N}$. The gradient of a function $f$ with respect to the variable $x$ is denoted by $\nabla_{x}f$. For $x,x^{\prime}\in{\mathbb{R}}_{>0}^{n}$, define an equivalence relation $\sim$ on ${\mathbb{R}}_{>0}^{n}$ by $x\sim x^{\prime}$ if and only if $\exists r>0,x=rx^{\prime}$.

Here, we briefly review OT between discrete distributions $\mu:=\sum_{i=1}^{N}a_{i}\delta_{x_{i}},\nu:=\sum_{j=1}^{M}b_{j}\delta_{y_{j}}$ where $a\in\Sigma_{N}:=\{p\in{\mathbb{R}}_{\geq 0}^{N}:\sum_{i=1}^{N}p_{i}=1\},b\in\Sigma_{M}$, $x_{i},y_{j}\in{\mathbb{R}}^{n}$, and $\delta_{x}$ is the Dirac delta at $x$. Given a cost function $c:{\mathbb{R}}^{n}\times{\mathbb{R}}^{n}(\ni(x,y))\rightarrow{\mathbb{R}}$, which represents the cost of transporting a unit of mass from $x$ to $y$, the original formulation of OT due to Monge seeks a map $T:\{x_{1},\ldots,x_{N}\}\rightarrow\{y_{1},\ldots,y_{M}\}$ that solves

		$\displaystyle\underset{T}{\rm minimize}\ \sum_{i\in[\![N]\!]}c(x_{i},T(x_{i}))$		(1)
		$\displaystyle\text{subject to}\ b_{j}=\sum_{i:T(x_{i})=y_{j}}a_{i},\ \forall j\in[\![M]\!].$

Especially when $M=N$ and $a=b=1_{N}/N$, the optimal map $T$ gives the optimal assignment for transporting agents with the initial states $\{x_{i}\}_{i}$ to the desired states $\{y_{j}\}_{j}$, and the Hungarian algorithm [10] can be used to solve (1). However, this method can be applied only to small problems because it has $O(N^{3})$ complexity.

On the other hand, the Kantorovich formulation of OT is a linear program (LP):

$$\underset{P\in{\mathcal{T}}(a,b)}{\rm minimize}\ \sum_{i\in[\![N]\!],j\in[\![M]\!]}C_{ij}P_{ij}$$

(2)

where $C_{ij}:=c(x_{i},y_{j})$ and

$${\mathcal{T}}(a,b):=\{P\in{\mathbb{R}}_{\geq 0}^{N\times M}:P1_{M}=a,\ P^{\top}1_{N}=b\}.$$

A matrix $P\in{\mathcal{T}}(a,b)$, which is called a coupling matrix, represents a transport plan where $P_{ij}$ describes the amount of mass flowing from $x_{i}$ towards $y_{j}$. In particular, when $M=N$ and $a=b=1_{N}/N$, there exists an optimal solution from which we can reconstruct an optimal map for (1) [11, Proposition 2.1]. However, similarly to (1), for a large number of agents and destinations, the problem (2) with $NM$ variables is challenging to solve.

In view of this, [9] employed an entropic regularization to (2):

$$\underset{P\in{\mathcal{T}}(a,b)}{\rm minimize}\ \sum_{i\in[\![N]\!],j\in[\![M]\!]}C_{ij}P_{ij}-\varepsilon H(P),$$

(3)

where $\varepsilon>0$ and the entropy of $P$ is defined by $H(P):=-\sum_{i,j}P_{ij}(\log(P_{ij})-1)$. Define the Gibbs kernel $K$ associated to the cost matrix $C=(C_{ij})$ as

$$K=(K_{ij})\in{\mathbb{R}}_{>0}^{N\times M},\ K_{ij}:=\exp\left(-C_{ij}/\varepsilon\right).$$

Then, a unique solution of the entropic OT (3) has the form $P^{*}=(\alpha^{*})^{\boxbslash}K(\beta^{*})^{\boxbslash}$ for two (unknown) scaling variables $(\alpha^{*},\beta^{*})\in{\mathbb{R}}_{>0}^{N}\times{\mathbb{R}}_{>0}^{M}$. The variables $(\alpha^{*},\beta^{*})$ can be efficiently computed by the Sinkhorn algorithm:

$$\alpha(k+1)=a\oslash[K\beta(k)],\ \beta(k)=b\oslash[K^{\top}\alpha(k)]$$

(4)

where

$$\lim_{k\rightarrow\infty}\alpha(k+1)^{\boxbslash}K\beta(k)^{\boxbslash}=P^{*},\ \forall\alpha(0)=\alpha_{0}\in{\mathbb{R}}_{>0}^{N}.$$

Now, let us introduce Hilbert’s projective metric

$${d_{\mathcal{H}}}(\beta,\beta^{\prime}):=\log\max_{i,j\in[\![M]\!]}\frac{\beta_{i}\beta^{\prime}_{j}}{\beta_{j}\beta^{\prime}_{i}},\ \beta,\beta^{\prime}\in{\mathbb{R}}_{>0}^{M},$$

(5)

which is a distance on the projective cone ${\mathbb{R}}_{>0}^{M}/{\sim}$ (see the Notation in Section I for $\sim$) and is useful for the convergence analysis of the Sinkhorn algorithm; see [11, Remark 4.12 and 4.14]. Indeed, for any $(\beta,\beta^{\prime})\in({\mathbb{R}}_{>0}^{M})^{2}$ and any $\bar{K}\in{\mathbb{R}}_{>0}^{N\times M}$, it holds

$${d_{\mathcal{H}}}(\bar{K}\beta,\bar{K}\beta^{\prime})\leq\lambda(\bar{K}){d_{\mathcal{H}}}(\beta,\beta^{\prime})$$

(6)

where

$$\lambda(\bar{K}):=\frac{\sqrt{\eta(\bar{K})}-1}{\sqrt{\eta(\bar{K})}+1}<1,\ \eta(\bar{K}):=\max_{i,j,k,l}\frac{\bar{K}_{ik}\bar{K}_{jl}}{\bar{K}_{jk}\bar{K}_{il}}.$$

Then it follows from (6) that

	$\displaystyle{d_{\mathcal{H}}}$	$\displaystyle(\beta(k+1),\beta^{*})={d_{\mathcal{H}}}(b\oslash[K^{\top}\alpha(k+1)],b\oslash[K^{\top}\alpha^{*}])$
		$\displaystyle={d_{\mathcal{H}}}(K^{\top}\alpha(k+1),K^{\top}\alpha^{*})$
		$\displaystyle\leq\lambda(K){d_{\mathcal{H}}}(\alpha(k+1),\alpha^{*})\leq\lambda^{2}(K){d_{\mathcal{H}}}(\beta(k),\beta^{*})$

which implies $V_{P}(\beta):={d_{\mathcal{H}}}(\beta,\beta^{*})$ is a Lyapunov function of (4), and $\lim_{k\rightarrow\infty}\beta(k)=\beta^{*}\in{\mathbb{R}}_{>0}^{M}/{\sim}$.

In this paper, we consider the problem of stabilizing agents efficiently to a given discrete distribution over dynamical systems. This can be formulated as Monge’s OT problem.

Note that the running cost $\ell_{i}$ depends not only on the state $x_{i}$ and the control input $u_{i}$, but also on the destination $x_{j}^{\sf d}$. Throughout this paper, we assume the existence of an optimal solution of OC problems. In addition, we assume that there exists a constant input $\bar{u}_{i}$ under which $x_{i}=x_{j}^{\sf d}$ is an equilibrium of (9). A necessary condition for the infinite horizon cost $c_{\infty}^{i}(x_{i}^{0},x_{j}^{\sf d})$ to be finite is that at $x_{i}=x_{j}^{\sf d}$ with $u_{i}=\bar{u}_{i}$, there is not a cost incurred, i.e., $\ell_{i}(x_{j}^{\sf d},\bar{u}_{i};x_{j}^{\sf d})=0$. If $B_{i}$ is square and invertible, $\bar{u}_{i}=B_{i}^{-1}(x_{j}^{\sf d}-A_{i}x_{j}^{\sf d})$ makes $x_{i}=x_{j}^{\sf d}$ an equilibrium.

The main difficulties of Problem 1 are as follows:

1. 1.

   In most cases, the infinite horizon OC problem $c_{\infty}^{{\color[rgb]{0,0,0}{i}}}(x_{i}^{0},x_{j}^{\sf d})$ is computationally intractable.

2. 2.

   In addition, given $c_{\infty}^{{\color[rgb]{0,0,0}{i}}}(x_{i}^{0},x_{j}^{\sf d}),\forall i,j\in[\![N]\!]$, the assignment problem needs to be solved, which leads to the high computational burden when $N$ is large.

To overcome these issues, we utilize the concept of MPC, which solves tractable finite horizon OC while satisfying the state constraint (10). Specifically, at each time, we address the OT problem whose cost function with the current state $\check{x}_{i}$ and the finite horizon $T_{h}\in{\mathbb{Z}}_{>0}$ is given by

	$\displaystyle c_{T_{h}}^{{\color[rgb]{0,0,0}{i}}}(\check{x}_{i},x_{j}^{\sf d}):=$	$\displaystyle\min_{u_{i}}\ \sum_{k=0}^{T_{h}-1}\ell_{i}(x_{i}(k),u_{i}(k);x_{j}^{\sf d})$
		$\displaystyle\text{subj. to \eqref{eq:nonlinear_dynamics}, {\color[rgb]{0,0,0}{\eqref{eq:state_constraint}}},}~{}~{}x_{i}(0)=\check{x}_{i},\ x_{i}(T_{h})=x_{j}^{\sf d}.$

Denote the first control in the optimal sequence by $u_{i}^{\rm MPC}(\check{x}_{i},x_{j}^{\sf d})$. From the viewpoint of the computation time for solving the OT problem, we relax the Problem 1 by the entropic regularization. It should be noted that, in challenging situations in which the number of agents is large and the sampling time is small, only a few Sinkhorn iterations are allowed. In what follows, we deal only with the case where just one Sinkhorn iteration is performed at each time. Nevertheless, by similar arguments, all of the results of this paper are still valid when more iterations are performed.

Based on the approximate solution $P(k)$ obtained by the Sinkhorn algorithm at time $k$, we need to determine a target state for each agent. Then, we introduce a set ${\mathbb{X}}\subset{\mathbb{R}}^{n}$ and a map $x_{\rm tmp}^{{\sf d},i}:{\mathbb{R}}_{\geq 0}^{N\times N}\rightarrow{\mathbb{X}}$ as a policy to determine a temporary target $x_{\rm tmp}^{{\sf d},i}(P(k))$ at time $k$ for agent $i$. Hereafter, assume that there exists a constant $r_{\rm upp}>0$ such that

$$\|x\|\leq r_{\rm upp},\ \forall x\in{\mathbb{X}}.$$

(13)

For example, if ${\mathbb{X}}$ is the convex hull of $\{x_{j}^{\sf d}\}_{j}$, we can take $r_{\rm upp}=\max_{j}\|x_{j}^{\sf d}\|$. A typical policy to approximate a Monge’s OT map from a coupling matrix $P$ is the so-called barycentric projection $x_{\rm tmp}^{{\sf d},i}(P)=N\sum_{j=1}^{N}P_{ij}x_{j}^{\sf d}$ [11, Remark 4.11].

For any given policy $x_{\rm tmp}^{{\sf d},i}$, we summarize the above strategy as the following dynamics where the Sinkhorn algorithm behaves as a dynamic controller.
Sinkhorn MPC:

	$\displaystyle x_{i}(k+1)=A_{i}x_{i}(k)+B_{i}u_{i}^{\rm MPC}\bigl{(}x_{i}(k),x_{\rm tmp}^{{\sf d},i}\left(P(k)\right)\bigr{)},$
	$\displaystyle\hskip 170.71652pt\forall i\in[\![N]\!],$		(14)
	$\displaystyle P(k)=\alpha(k+1)^{\boxbslash}K(x(k))\beta(k)^{\boxbslash},$		(15)
	$\displaystyle\alpha(k+1)=1_{N}/N\oslash\left[K(x(k))\beta(k)\right],$		(16)
	$\displaystyle\beta(k)=1_{N}/N\oslash\left[K(x(k))^{\top}\alpha(k)\right],$		(17)
	$\displaystyle x_{i}(0)=x_{i}^{0},\ \alpha(0)=\alpha_{0},$

where

$\displaystyle K_{ij}(x):=\exp\left(-\frac{c_{T_{h}}^{{\color[rgb]{0,0,0}{i}}}(x_{i},x_{j}^{\sf d})}{\varepsilon}\right),\ x=[x_{1};\cdots;x_{N}]\in{\mathbb{R}}^{nN},$

and the initial value $\alpha_{0}$ is arbitrary. $\triangle$
Note that, under the assumption that for all $i\in[\![N]\!]$, $B_{i}$ is square and invertible, Sinkhorn MPC is obviously recursively feasible, and the dynamics (14) satisfies the state constraint (10) for all $i\in[\![N]\!]$.

This section gives examples for Sinkhorn MPC with an energy cost

$$\ell_{i}(x_{i},u_{i};x_{j}^{\sf d})=\|u_{i}-B_{i}^{-1}(x_{j}^{\sf d}-A_{i}x_{j}^{\sf d})\|^{2},$$

(18)

and ${\mathcal{X}}={\mathbb{R}}^{n}$. Then, the dynamics under Sinkhorn MPC can be written as follows [12, Section 2.2, pp. 37-39]:

	$\displaystyle x_{i}(k+1)=\bar{A}_{i}x_{i}(k)+(I-\bar{A}_{i})x_{\rm tmp}^{{\sf d},i}(P(k)),$		(19)
	$\displaystyle u_{i}^{\rm MPC}(x_{i},{\color[rgb]{0,0,0}{\hat{x}}})=-B_{i}^{\top}(A_{i}^{\top})^{T_{h}-1}G_{i,T_{h}}^{-1}A_{i}^{T_{h}}(x_{i}-{\color[rgb]{0,0,0}{\hat{x}}})$
	$\displaystyle\hskip 71.13188pt+B_{i}^{-1}({\color[rgb]{0,0,0}{\hat{x}}}-A_{i}{\color[rgb]{0,0,0}{\hat{x}}}),\ \forall i\in[\![N]\!],\ {\color[rgb]{0,0,0}{\forall\hat{x}\in{\mathbb{R}}^{n}}}$

with (15)–(17) where

	$\displaystyle K_{ij}(x)=\exp\left(-\frac{\|x_{i}-x_{j}^{\sf d}\|_{{\mathcal{G}}_{i}}^{2}}{\varepsilon}\right),$
	$\displaystyle{\mathcal{G}}_{i}:=(A_{i}^{T_{h}})^{\top}G_{i,T_{h}}^{-1}A_{i}^{T_{h}},\ G_{i,T_{h}}:=\sum_{k=0}^{T_{h}-1}A_{i}^{k}B_{i}B_{i}^{\top}(A_{i}^{\top})^{k},$
	$\displaystyle\bar{A}_{i}:=A_{i}-B_{i}B_{i}^{\top}(A_{i}^{\top})^{T_{h}-1}G_{i,T_{h}}^{-1}A_{i}^{T_{h}}.$

In the examples below, we use the barycentric target $x_{\rm tmp}^{{\sf d},i}(P)=N\sum_{j=1}^{N}P_{ij}x_{j}^{\sf d}$.

First, consider the case where

$$A_{i}=\begin{bmatrix}1.2&0.13\\ -0.05&1.1\end{bmatrix},\ B_{i}=0.1I_{2},\ \forall i\in[\![N]\!],$$

(20)

and set $N=150,\ \varepsilon=1.0,\ T_{h}=10,\ \alpha_{0}=1_{N}$. For given initial and desired states, the trajectories of the agents governed by (19) with (15)–(17) are illustrated in Fig. 1. It can be seen that the agents converge sufficiently close to the target states. The computation time for one Sinkhorn iteration is about $0.0063$ ms, $0.030$ ms, and $0.08$ ms for $N=150,500,800$, respectively, with MacBook Pro with Intel Core i5. On the other hand, solving the linear program (2) with MATLAB linprog takes about $0.12$ s, $6.4$ s, and $66$ s for $N=150,500,800$, respectively, and is thus not scalable. Hence, Sinkhorn MPC contributes to reducing the computational burden.

Next, we investigate the effect of the regularization parameter $\varepsilon$ on the behavior of Sinkhorn MPC. To this end, consider a simple case where $N=10,\ T_{h}=10$, and

$$A_{i}=1,\ B_{i}=0.1,\ \forall i\in[\![N]\!].$$

(21)

Then the trajectories of the agents with $\varepsilon=0.5,1.0$ are shown in Fig. 2. As can be seen, the overshoot/undershoot is reduced for larger $\varepsilon$ while the limiting values of the states deviate from the desired states. In other words, the parameter $\varepsilon$ reflects the trade-off between the stationary and transient behaviors of the dynamics under Sinkhorn MPC.

In this section, we elucidate ultimate boundedness and asymptotic stability of Sinkhorn MPC with the energy cost (18) and ${\mathcal{X}}={\mathbb{R}}^{n}$. Hereafter, we assume the invertibility of $B_{i}$.