\section{Introduction}\label{sec:intro}
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~\cite{Chung2018,Chowdhury2015,Alasseur2020}.
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~\cite{Villani2003} between the empirical distribution based on the state of the agents and the target distribution over dynamical \blu{systems}.
The OT problem over dynamical systems consists of finding an assignment of agents to targets and control inputs that drive the agents to the assigned targets in order to minimize the total cost of interest.
The difficulty of this problem lies in the large scale of the collective dynamics.





{\it Literature review:}
The assignment problem has been extensively studied in the context of combinatorial optimization, and many methods to find the optimal assignment have been proposed such as the well-known Hungarian algorithm~\cite{Kuhn1955} and auction algorithm~\cite{Bertsekas1992}.
These algorithms have been applied to multi-agent assignment problems; see e.g., \cite{Yu2014,Mosteo2017} and references therein.
In the literature, the dynamics of agents are simplified as the single integrator dynamics, and easily computable assignment costs, e.g., distance-based cost, are considered in general.
On the other hand, when considering more general dynamics and cost functions for the stabilization to targets, it is difficult to obtain the associated assignment costs and optimal controls. This is because, in most cases, infinite horizon optimal control (OC) problems stabilizing agents to desired targets are computationally intractable.
A promising approach to overcome this problem is model predictive control~(MPC)~\cite{Mayne2014}, in which the current control input is determined by solving, at each sampling instant, a finite horizon OC problem using the current state as the initial state. For example \subm{in \cite{Morgan2016},} MPC is used to solve a finite horizon assignment problem over \subm{dynamical systems} in real time.
Now it is important to emphasize that when performing MPC for \subm{a dynamic} OT problem, \subm{it is desirable to update the target assignment for agents at each time as well as control inputs.}
However, when the number of the agents is large, solving the assignment problem at each sampling instant is computationally very expensive even with \subm{the Hungarian algorithm}.
\subm{Even worse, the changes of the assignment along the controlled state trajectories are not continuous, and this makes it difficult to ensure the stability of the dynamics under MPC.}





On the other hand, recently, a different approach to solve a dynamical assignment problem using OT theory has attracted much attention~\cite{Chen2018,Bakshi2020,Krishnan2018}.
In this approach, a large population limit is considered, and infinitely many agents are represented as a probability density of the state of a single system. 
Then, the dynamical assignment problem boils down to a density control problem~\cite{Chen2017,de2021,Ito2022bridge} finding a feedback control law that steers an initial \subm{state} density to a target density with minimum cost.
Consequently, this approach can avoid the difficulty due to the large scale of the collective dynamics.
\blu{Nevertheless}, it has the drawback that even for linear systems, the density control requires to solve a nonlinear partial differential equation such as the Monge-Amp\`{e}re equation or the Hamilton-Jacobi-Bellman equation, which are generally difficult to solve.



{\it Contributions:}
With this in mind, we deal with the collective dynamics directly without taking the number of agents to infinity, but utilizing the results of computational OT.
Specifically, in \cite{Cuturi2013}, \red{several favorable computational properties of an entropy-regularized version of OT are highlighted. In particular,} entropy-regularized OT problems can be solved efficiently by an iterative algorithm called the Sinkhorn algorithm.
Inspired by this, we propose a dynamical transport algorithm integrating MPC and the Sinkhorn algorithm, which we call {\em Sinkhorn MPC}.
This method incorporates the Sinkhorn iterations into MPC as a dynamic controller and can be seen as simultaneously solving an assignment problem while executing control actions.
The contributions of this paper coming from the introduction of Sinkhorn MPC are as follows:
\begin{enumerate}
	\item[1)] By combining MPC and the Sinkhorn algorithm, the computational effort for determining destinations of agents at each time is reduced substantially;
 	\item[2)] Thanks to the smoothing effect of the entropy regularization, we reveal the global convergence property of Sinkhorn MPC with a \subm{sufficiently} large number of \las{Sinkhorn} iterations;
  \item[3)] For a quadratic control cost, we show the ultimate boundedness and the local asymptotic stability for Sinkhorn MPC without the assumption of the number of \blu{Sinkhorn} iterations.
\end{enumerate}



Compared to a preliminary version of this work~\cite{Ito2022sinkhorn}, the most notable improvement of the current paper is that we consider continuous-time systems and derive 2) the global convergence property for Sinkhorn MPC, which is one of the most crucial properties of dynamical transport algorithms. 
In addition, we provide several illustrative examples, which describe the usefulness of Sinkhorn MPC, and we give the proof of Lemma~\ref{lemma:exponential_convergence}, which is omitted in the preliminary version.




\textit{Organization:}
The remainder of this paper is organized as follows. In Section~\ref{sec:optimal_transport}, we introduce OT between discrete distributions. In Section~\ref{sec:formulation}, we provide the problem formulation.
In Section~\ref{sec:proposed}, we describe the idea of Sinkhorn MPC.
\modi{In Section~\ref{sec:example}, numerical examples illustrate the utility of the proposed method.
Section~\ref{sec:global} is devoted to the global convergence analysis of the proposed method.
In Section~\ref{sec:linear_system}, for a quadratic control cost, we investigate \blu{fundamental} properties of Sinkhorn MPC, \blu{such as local asymptotic stability}.}
Some concluding remarks are given in Section~\ref{sec:conclusion}.

\textit{Notation:}
Let $ \bbR $ denote the set of real numbers.
The set of all positive (resp. nonnegative) vectors in $ \bbR^n $ is denoted by $ \bbR_{> 0}^n $ (resp. $ \bbR_{\ge 0}^n $). We use similar notations for the set of all real matrices $ \bbR^{m\times n} $ and integers $ \bbZ $, respectively.
The set of integers $ \{1,\ldots,N\} $ is denoted by $ \bbra{N} $.
The Euclidean norm is denoted by $ \|\cdot \| $.
For a positive semidefinite matrix $ A $, denote $\|x\|_A := (x^\top A x)^{1/2}$.
The identity matrix of size $n$ is denoted by $I_n$ \subm{or $ I $ when its size is clear in the context}. The matrix norm induced by the Euclidean norm is denoted by $ \| \cdot \|_2 $.
\rev{For vectors $ x_1,\ldots,x_m \in \bbR^n $, a collective vector $ [x_1^\top \ \cdots \ x_m^\top]^\top \in \bbR^{nm} $ is denoted by $ [x_1 ; \ \cdots \ ; x_m] $.}
For $ A = [a_1  \ \cdots \ a_n] \in \bbR^{m\times n}$, we write $ {\rm vec} (A) := \rev{[a_1; \ \cdots \ ;a_n]} $.
For $ \alpha = [\alpha_1 \ \cdots \ \alpha_N ]^\top \in \bbR^N $, the diagonal matrix with diagonal entries $ \{\alpha_i\}_{i=1}^N $ is denoted by \rev{$ \alpha^\diagbox $}.
The block diagonal matrix with diagonal entries $ \{A_i\}_{i=1}^N, A_i \in \bbR^{m\times n} $ is denoted by \rev{$ \{A_i\}_i^\diagbox $}.
Especially when $ A_i = A, \forall i $, $ \{A_i\}_i^\diagbox $ is also denoted by \rev{$  A^{\diagbox,N} $}.
Let $ (\subm{\calM},d) $ be a metric space.
The open ball of radius $ r > 0 $ centered at $ x \in \subm{\calM} $ is denoted by $ B_r(x) := \{y\in \subm{\calM} : d(x,y) < r \} $.
The element-wise division of $ a, b \in \bbR_{>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 $ {\bf 1}_N $.
The gradient of a function $ f $ with respect to the variable $ x $ is denoted by $ \nabla_x f $.
\rev{For $ x, x' \in \bbR_{> 0}^n $, define an equivalence relation $ \sim $ on $ \bbR_{> 0}^n $ by $ x \sim x' $ if and only if $ \exists r> 0, x = r x' $.}







\section{Background on optimal transport}\label{sec:optimal_transport}
Here, we briefly review OT between discrete distributions $ \mu := \sum_{i=1}^{N} {\mathbf a}_i \delta_{x_i}, \nu := \sum_{j=1}^{M} {\mathbf b}_j \delta_{y_j} $ where ${\mathbf a} \in \Sigma_N := \{ p\in \bbR_{\ge 0}^N : \sum_{i=1}^{N} p_i = 1 \}, {\mathbf b}\in \Sigma_M $, $ x_i, y_j \in \bbR^n $, and $ \delta_x $ is the Dirac delta at $ x $.
Given a cost function $ c : \bbR^n \times \bbR^n (\ni (x,y)) \rightarrow \bbR $, 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 $ \subm{{\bf T}} : \{x_1,\ldots,x_N\} \rightarrow \{y_1,\ldots,y_M\}$ that solves
\begin{equation}\label{prob:Monge}
	\begin{aligned}
	&\underset{\subm{{\bf T}}}{\rm minimize} \ \sum_{i \in \bbra{N}} c(x_i, \subm{{\bf T}}(x_i)) \\
	&\text{subject to} \ {\mathbf b}_j = \sum_{i : \subm{{\bf T}}(x_i) = y_j} {\mathbf a}_i, \ \forall j \in \bbra{M} .
	\end{aligned}
\end{equation}
Especially when $ M = N $ and $ {\mathbf a} = {\mathbf b} = {\bf 1}_N /N $, the optimal map $ \subm{{\bf T}} $ gives the optimal assignment for transporting agents with the initial states $ \{x_i\}_i $ to the desired states $ \{y_j\}_j $, \red{and then for example, the Hungarian algorithm can be adopted to solve \eqref{prob:Monge}.}
\red{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 \subm{program}:
\begin{equation}\label{prob:Kantorovich}
	\underset{P\in \calT({\mathbf a},{\mathbf b})}{\rm minimize} \ \sum_{i\in \bbra{N},j\in \bbra{M}} C_{ij} P_{ij}
\end{equation}
where $ C_{ij} := c(x_i, y_j) $ and
\[
	\calT({\mathbf a},{\mathbf b}) := \left\{ P \in \bbR_{\ge0}^{N\times M} : P {{\bf 1}}_M = {\mathbf a}, \ P^\top {\bf 1}_N = {\mathbf b}   \right\}.
\]
A matrix $ P \in \calT({\mathbf a},{\mathbf 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 $ {\mathbf a} = {\mathbf b} = {\mathbf 1}_N/N $, there exists an optimal solution from which we can reconstruct an optimal map for \subm{Monge's problem}~\eqref{prob:Monge} \cite[Proposition~2.1]{Peyre2019}.
\red{However, similarly to \eqref{prob:Monge}, for a large number of agents and destinations, the problem~\eqref{prob:Kantorovich} with $ N M $ variables is challenging to solve.}


In view of this, \cite{Cuturi2013} employed entropy regularization to \eqref{prob:Kantorovich}:
\begin{equation}\label{prob:Sinkhorn}
		\underset{P\in \calT({\mathbf a},{\mathbf b})}{\rm minimize} \ \sum_{i\in \bbra{N},j\in \bbra{M}} C_{ij} P_{ij} - \varepsilon {\rm H}(P),
\end{equation}
where $ \varepsilon > 0 $ is a regularization parameter and the entropy of $ P $ is defined by ${\rm H}(P) := - \sum_{i,j} P_{ij} (\log (P_{ij}) - 1)$. 
Define the Gibbs kernel $ K $ associated \subm{with} the cost matrix $ C = (C_{ij}) $ as
\[
	K = (K_{ij}) \in \bbR_{> 0}^{N\times M}, \ K_{ij} := \exp\left(  - C_{ij}/\varepsilon \right) .
\]
Then, a unique solution of the entropic OT problem~\eqref{prob:Sinkhorn} has the form
\begin{equation}\label{eq:opt_form}
	P^* =  (\alpha^*)^\diagbox K  (\beta^*)^\diagbox,
\end{equation}
where the two scaling variables $ (\alpha^*,\beta^*) \in \bbR_{>0}^N \times \bbR_{>0}^M $ are determined by
\begin{equation}\label{eq:opt_scaling_condition}
	\alpha^* = {\mathbf a} \oslash [K\beta^*], \ \beta^* = {\mathbf b} \oslash [K^\top \alpha^*] .
\end{equation}
The variables $ (\alpha^*, \beta^*) $ can be efficiently computed by the Sinkhorn algorithm:
\begin{equation}\label{eq:Sinkhorn_algorithm}
	\alpha[k+1] = {\mathbf a} \oslash [K\beta[k]], \ \beta[k] = {\mathbf b} \oslash \left[K^\top \alpha[k]\right]
\end{equation}
where
\[
	\lim_{k\rightarrow \infty} \alpha[k+1]^\diagbox K \beta  [k]^\diagbox  =  P^*, \ \forall \alpha[0] = \alpha_0 \in \bbR_{>0}^N. 
\]

\rev{Now, let us introduce Hilbert's projective metric}
\begin{equation}\label{eq:hilbert_metric}
	{d_\calH} (\beta, \beta') := \log \max_{i,j\in \bbra{M}} \frac{\beta_i \beta'_j}{\beta_j \beta'_i}, \ \beta, \beta' \in \bbR_{>0}^M ,
\end{equation}
\rev{which is a distance on the projective cone $ \bbR_{>0}^M / {\sim} $ (see the Notation in Section~\ref{sec:intro} for $ \sim $)} and is useful for the convergence analysis of the Sinkhorn algorithm; see~\cite[Remark~4.12 \rev{and 4.14}]{Peyre2019}.
Indeed, for any $ (\beta, \beta') \in (\bbR_{> 0}^M )^2 $ and any $ \bar K \in \bbR_{>0}^{N\times M} $, it holds
\begin{equation}\label{ineq:hilbert_metric}
		{d_\calH} (\bar K\beta, \bar K\beta') \le \lambda (\bar K) {d_\calH} (\beta, \beta')
\end{equation}
where
\begin{equation*}
	\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}} .
\end{equation*}
Then it follows from \eqref{ineq:hilbert_metric} that
\begin{align*}
	&{d_\calH} (\beta[k+1], \beta^*) = {d_\calH} \left({\mathbf b}\oslash [K^\top \alpha[k+1]], {\mathbf b}\oslash [K^\top \alpha^*]\right) \\
	&={d_\calH} (K^\top \alpha[k+1], K^\top \alpha^*) \\
	&\le \lambda (K) {d_\calH} (\alpha [k+1], \alpha^*) \le \lambda^2 (K) {d_\calH} (\beta[k], \beta^*) \subm{,}
\end{align*}
which implies $ V_{\subm{\calH}} (\beta) := {d_\calH} (\beta, \beta^*) $ \subm{serves as} a Lyapunov function of \eqref{eq:Sinkhorn_algorithm}, and $ \lim_{k\rightarrow \infty} \beta [k] = \beta^* \in \bbR_{>0}^M / {\sim} $.





\section{Problem formulation}\label{sec:formulation}
In this paper, we consider the problem of \rev{stabilizing agents efficiently to a given discrete distribution over dynamical systems.
This can be formulated as Monge's OT problem.
\begin{problem}
	\label{prob:main_continuous}
	Given initial and desired states $\{x_i^0\}_{i=1}^N, $ $ \{x_j^{\sf d}\}_{j=1}^N  \in (\mathbb{R}^{n})^N$,
	find control inputs $\{u_i\}_{i=1}^N$ and a permutation $ \sigma : \bbra{N} \rightarrow \bbra{N} $ that solve
	\begin{align}
		&\underset{\sigma }{\rm minimize} ~~ \sum_{i\in \bbra{N}} c_\infty^i (x_i^0, x_{\sigma(i)}^{\sf d} ) .  \label{eq:infinite_horizon_cost}
	\end{align}
	Here, the cost function $ c_\infty^i : \bbR^n \times \bbR^n \rightarrow \bbR $ is defined by
	\begin{align}
		& &&\hspace{-2cm} c_\infty^i (x_i^0, x_j^{\sf d}) := \min_{u_i} \ \int_0^\infty  \ell_i (x_i (t), u_i (t) ; x_j^{\sf d}) {\rm d} t \label{eq:value_infty}\\
		&\hspace{0.5cm} \text{\rm subject to} &&\red{ \dot{x}_i(t)  = A_i x_i (t)  + B_i u_i (t) }, \label{eq:linear_dynamics_conti}\\
		& &&x_i (t) \in \bbX_i \subseteq \bbR^n , \ \forall t \ge 0 , \label{eq:state_constraint} \\ 
		& &&u_i (t) \in \bbU_i \subseteq \bbR^m , \ \forall t \ge 0 , \label{eq:input_constraint} \\ 
		& &&x_i(0) = x_{i}^0, \label{eq:constraint_initial} \\
		& &&\lim_{t\rightarrow \infty} x_i(t) =  x_j^{\sf d},\label{eq:constraint_infty}
	\end{align}
	where $ x_i(t) \in \bbR^n $ denotes the state of the agent $ i $, and $ A_i\in \bbR^{n\times n}, B_i \in \bbR^{n\times m} $.
	\hfill $ \diamondsuit $
\end{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} $. \rev{Throughout this paper, we assume the existence of an optimal solution of OC problems. \red{In addition, we assume that there exists a \blu{constant} input $ \bar{u}_{\blu{ij}} $ under which $ x_i = x_j^{\sf d} $ is an equilibrium of \eqref{eq:linear_dynamics_conti}.}
\red{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} $ \las{and at least one such input} $ u_{i} = \bar{u}_{\blu{ij}} $, there is not a cost incurred, i.e., $ \ell_i(x_j^{\sf d}, \bar{u}_{\blu{ij}} ; x_j^{\sf d}) = 0 $.}
\blu{For instance,} if $ B_i $ is square and invertible, $ \bar{u}_{\blu{ij}} =  - B_i^{-1} A_ix_j^{\sf d} $ makes $ x_i = x_j^{\sf d} $ an equilibrium.
}


In most cases, the infinite horizon OC problem $ c_\infty^{\magenta{i}} (x_i^0,x_j^{\sf d}) $ is computationally intractable. To avoid this difficulty, we use MPC, which solves a tractable finite horizon OC problem with a prediction horizon $ T_{\rm h} > 0 $ at each time:
\begin{align}
	c_{T_{\rm h}}^{\magenta{i}} (\check{x}_i, x_j^{\sf d}) := &\min_{u_i} \ \int_{0}^{T_{\rm h}} \ell_i (x_i(t), u_i (t); x_j^{\sf d}) {\rm d} t\label{eq:cost_finite_conti}\\
	&\text{subj. to  \eqref{eq:linear_dynamics_conti}\text{--}\magenta{\eqref{eq:input_constraint}},}~~ x_i(0) = \check{x}_i, \  x_i (T_{\rm h}) = x_{j}^{\sf d} . \nonumber
\end{align}
Denote the first control in the optimal sequence of the above problem by $ u_i^{\rm MPC} (\check{x}_i, x_j^{\sf d}) $.
Also for $ \check{x} = [\check{x}_1 ; \cdots ; \check{x}_N] \in \bbR^{nN} $, denote by $ \sigma(\cdot; \check{x}) $ the optimal permutation of the following problem:
\begin{align}
	&\underset{\sigma }{\rm minimize} ~~ \sum_{i\in \bbra{N}} c_{T_{\rm h}}^i (\check{x}_i, x_{\sigma(i)}^{\sf d} ) .  \label{eq:finite_horizon_OT}
\end{align}
Then the dynamics \eqref{eq:linear_dynamics_conti} under MPC for Problem~\ref{prob:main_continuous} is given by
\begin{equation}\label{eq:permutation_dynamics}
  \dot{x}_i (t) = A_i x_i (t) + B_i u_i^{\rm MPC} \left(x_i(t), x_{\sigma(i;x(t))}^{\sf d} \right) , \ \forall i\in \bbra{N} ,
\end{equation}
where $ x(t) := [x_1(t);\cdots ; x_N(t)] $.
Note that along the trajectory $ x(t) $, at several times, the permutation $ \sigma(\cdot;x(t)) $ changes in general.
The state-dependent permutation $ \sigma(\cdot; x(t)) $ is expected to \blu{reduce the cost accumulated during the transport} \las{more} than the permutation $ \sigma(\cdot; x^0) $ that is fixed at the initial time $ t=0 $.
Despite the merit, the state-dependency of the permutation poses the following computational and theoretical difficulties:
\begin{itemize}
	\item Solving the assignment problem~\eqref{eq:finite_horizon_OT} at each time leads to the high computational burden when $ N $ is large;
 	\item The optimal permutation $ \sigma(\cdot; x) $ is not continuous in $ x $. That is, the target states $ \{x_{\sigma(i;x(t))}^{\sf d} \}_i $ for the agents change discontinuously along the trajectory $ x(t) $, and this makes \blu{it} difficult to ensure the \blu{convergence} of the dynamics~\eqref{eq:permutation_dynamics}.
\end{itemize}
In the remainder of this paper, we reveal that entropy regularization mitigates the above problems.








\section{MPC with entropy-regularized optimal transport}\label{sec:proposed}
\subsection{Introduction of the entropy regularization to MPC}
Now, to avoid the issues observed in the previous section, we employ the entropy regularization.
To this end, we first consider the linear program:
\begin{equation}\label{eq:equi_LP}
	\underset{P\in \calT(\comm{\frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N},\frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N}})}{\rm minimize} \ \sum_{i,j\in \bbra{N}} C_{ij} (x) P_{ij},
\end{equation}
where $ C_{ij} (x) := c_{T_{\rm h}}^i (x_i,x_j^{\sf d}), \ x = [x_1;\cdots;x_N] $. Then as mentioned in Section~\ref{sec:optimal_transport}, the optimal permutation $ \sigma $ can be obtained by the optimal permutation matrix $ P^\sigma $ of \eqref{eq:equi_LP} satisfying $ P_{ij}^\sigma = 1/N $ if $ j= \sigma(i) $, and $ 0 $, otherwise.
Next, we introduce the entropy regularization to \eqref{eq:equi_LP} as in \eqref{prob:Sinkhorn}.
Then, based on the optimal coupling $ P^* $ of the entropic OT problem, we determine a target state for each agent.
Specifically, we introduce a map $ x_i^{\rm tmp} : \bbR_{\ge 0}^{N\times N} \rightarrow \calX (\subset \bbR^n) $ as a policy to determine a temporary target $ x_i^{\rm tmp} (P^*) $ for agent $ i $.
\comm{We call $ x_i^{\rm tmp} $ a navigator function.}
A typical \comm{navigator function} to approximate Monge's OT map from a coupling matrix $ P $ is the so-called barycentric projection~\cite[Remark~4.11]{Peyre2019}:
\begin{equation}\label{eq:barycentric_target}
x_i^{\rm tmp} (P) = N \sum_{j=1}^{N} P_{ij} x_j^{\sf d}, \ P \in \bbR_{\ge 0}^{N\times N} .
\end{equation}
Note that, for a permutation matrix $ P^\sigma $, it holds $ N\sum_{j=1}^N P_{ij}^\sigma x_j^{\sf d} = x_{\sigma(i)}^{\sf d} $.
\comm{Fig.~\ref{fig:image} illustrates the states of three agents $ \{x_i\} $, destinations $ \{x_{\subm{j}}^{\sf d}\} $, and \blu{temporary} targets $ \{x_i^{\rm tmp}(P)\} $ determined by the barycentric projection~\eqref{eq:barycentric_target} for a given coupling matrix $ P $.}


Now, we propose to use the control law \las{$ u_i^{\rm MPC} \bigl(x_i(t), $ $ x_i^{\rm tmp} (P^*(x(t)))\bigr) $} where
\begin{equation}
	P^* (x) := \mathop{\rm arg~min}\limits_{P\in \calT(\frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N}, \frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N})} \ \sum_{i,j\in \bbra{N}} C_{ij} (x) P_{ij} - \varepsilon {\rm H}(P). \label{eq:P_opt_sink}
\end{equation}
In summary, for any given \subm{navigator function}~$ x_i^{\rm tmp} $ and $ \varepsilon > 0 $, the dynamics of the agents are written as
\begin{align}
  &\dot{x}_i (t) = A_i x_i (t) + B_i u_i^{\rm MPC} \left(x_i(t), x^{\rm tmp}_{i} \Bigl(P^*(x(t))\Bigr) \right) , \nonumber \\
  &\hspace{6cm} \forall i\in \bbra{N} , \label{eq:mpc_dynamics_conti}\\
  &x_i (0) = x_i^0, \ \forall i \in \bbra{N} . \nonumber 
\end{align}
The entropy regularization enables to use the Sinkhorn algorithm~\eqref{eq:Sinkhorn_algorithm}, which contributes to \blu{reducing} the computational burden of determining target states at each time.
In addition, we will see that the entropy regularization also enables to analyze the global convergence property of \eqref{eq:mpc_dynamics_conti} in Section~\ref{sec:global}.




\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.4]{figures/image3.eps}
	\caption{\comm{\subm{Three} agents $ \{x_i\} $, desired states $ \{x_{\subm{j}}^{\sf d}\} $, and temporary targets $ \{x_i^{\rm tmp} (P)\} $ determined by the barycentric projection~\eqref{eq:barycentric_target} for a given coupling matrix $ P $.}}
	\label{fig:image}
\end{figure}



\subsection{Integrating MPC and the Sinkhorn algorithm}\label{subsec:proposed}
In the previous subsection, \subm{it was implicitly} assumed that at each time, the optimal coupling $ P^*(x(t)) $ is available for determining temporary targets for agents.
The Sinkhorn algorithm achieves a speed-up in the computation of an optimal coupling.
However, in challenging situations in which the number of agents is very large and the sampling time is small, only a few Sinkhorn iterations are allowed. In such a case, an available approximate coupling matrix may not be close enough to the optimal coupling, and therefore the above assumption is not valid.
To address this issue, in this subsection, we propose to integrate MPC and the Sinkhorn algorithm.
Since the Sinkhorn algorithm works in discrete time, we consider a time-discretized version of \eqref{eq:linear_dynamics_conti} obtained by e.g., a zero-order hold discretization\footnote{Throughout this paper, we use bold symbols for discrete-time systems while we use italic letters for continuous-time systems.}:
\begin{equation}\label{eq:linear_dynamics_disc}
	{\boldsymbol x}_i [k+1]  = {\boldsymbol A}_i {\boldsymbol x}_i [k]  + {\boldsymbol B}_i {\boldsymbol u}_i [k], \ k\in \bbZ_{\ge 0} .
\end{equation}
Then, the cost function $ {\boldsymbol c}_{\tau_{\rm h}}^i $ with a finite horizon $ \tau_{\rm h} \in \bbZ_{>0} $ is defined by
\begin{align}
	& &&\hspace{-1cm} {\boldsymbol c}_{\tau_{\rm h}}^i (\check{x}_i, x_j^{\sf d}) := \min_{{\boldsymbol u}_i} \ \sum_{k=0}^{\tau_{\rm h}-1}  {\boldsymbol \ell}_i \left({\boldsymbol x}_i [k], {\boldsymbol u}_i [k] ; x_j^{\sf d} \right) \label{eq:cost_disc}\\
	&\hspace{-0.1cm} \text{\subm{subj. to}} &&\eqref{eq:linear_dynamics_disc},\nonumber\\
	& &&{\boldsymbol x}_i[k] \in \bbX_i, \ {\boldsymbol u}_i [k] \in \bbU_i , \ \forall k \in \subm{\bbra{\tau_{\rm h} -1}\cup \{0\}} , \nonumber \\ 
	& &&{\boldsymbol x}_i[0] = \check{x}_{i}, \label{eq:constraint_initial_disc} \\
	& &&{\boldsymbol x}_i[\tau_{\rm h}] =  x_j^{\sf d}. \label{eq:constraint_end_disc}
\end{align}
Denote the first control in the optimal sequence of the above problem by $ {\boldsymbol u}_i^{\rm MPC} (\check{x}_i, x_j^{\sf d}) $. Let $ {\boldsymbol x}[k] := [{\boldsymbol x}_1[k];\cdots; {\boldsymbol x}_N[k]] $ and
\begin{equation}
	{\boldsymbol P}^* (x) := \mathop{\rm arg~min}\limits_{P\in \calT(\frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N}, \frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N})} \ \sum_{i,j\in \bbra{N}} {\boldsymbol C}_{ij} (x) P_{ij} - \varepsilon {\rm H}(P), \label{eq:P_opt_sink_disc}
\end{equation}
where $ {\boldsymbol C}_{ij} (x) := {\boldsymbol c}_{\tau_{\rm h}}^i (x_i, x_j^{\sf d}) $.
Note that if $ {\boldsymbol c}_{\tau_{\rm h}}^i (x_{\blu{i}},x_j^{\sf d}) $ is continuous in $ x_{\blu{i}} $ for all $ i\in \bbra{N} $, then from the relations~\eqref{eq:opt_form},~\eqref{eq:opt_scaling_condition}, $ {\boldsymbol P}^* $ is continuous.
Hence, it is expected that if we take a sampling time for \eqref{eq:linear_dynamics_disc} appropriately so that the difference between $ {\boldsymbol x}_i [k+1] $ and $ {\boldsymbol x}_i [k] $ is small, then the scaling variables $ (\alpha^*[k+1], \beta^*[k+1]) $ for $ {\boldsymbol P}^*({\boldsymbol x}[k+1]) $ are close to the variables $ (\alpha^*[k],\beta^*[k]) $ for $ {\boldsymbol P}^* ({\boldsymbol x}[k]) $. This implies that $ (\alpha^*[k],\beta^*[k]) $ yield good initial estimates of $ (\alpha^*[k+1], \beta^*[k+1]) $.


Based on this observation, we present a dynamical transport algorithm integrating MPC and the Sinkhorn algorithm.
Let $ S[k] \in \bbZ_{>0} $ be the number of \blu{Sinkhorn} iterations at time $ k $.
\rev{For any given \comm{navigator function} $ x_i^{\rm tmp} $ and $ \varepsilon > 0 $,} the proposed algorithm, which we call {\em Sinkhorn MPC} is given as the following dynamics where the Sinkhorn algorithm behaves as a dynamic controller.\\
{\bf Sinkhorn MPC:}
\begin{align}
&{\boldsymbol x}_i[k+1] = {\boldsymbol A}_i {\boldsymbol x}_i [k] + {\boldsymbol B}_i {\boldsymbol u}_i^{\rm MPC} \bigl({\boldsymbol x}_i [k], x_i^{\rm tmp} \left(P[k] \right)  \bigr),\nonumber\\
&\hspace{6cm} \forall i \in \bbra{N}, \label{eq:SMPC_x} \\
&P[k] = \alpha\left[k,S[k]+1\right]^\diagbox {\boldsymbol K}({\boldsymbol x}[k]) \beta\left[k,{S}[k]\right]^\diagbox, \label{eq:SMPC_P} \\
&\text{Sinkhorn iterations:} \nonumber\\
&\begin{cases}
\alpha\left[k, {l+1}\right] = {\bf 1}_N/N \oslash \left[{\boldsymbol K}({\boldsymbol x}[k])\beta[k,l] \right], \\%\label{eq:SMPC_alpha} \\
\beta[k, l] = {\bf 1}_N/N \oslash \left[{\boldsymbol K}({\boldsymbol x}[k])^\top \alpha[k,l] \right] ,
\end{cases} \ l \in \bbra{S[k]}, \label{eq:SMPC_sink}\\
&\alpha[k+1,1] = \alpha[k,S[k]+1],  \label{eq:SMPC_next} \\ 
&{\boldsymbol x}_i [0] = x_i^0, \ \alpha[0,1] = \alpha_0, \nonumber 
\end{align}
where
\begin{align*}
{\boldsymbol K}_{ij}(x) := \exp\left(- \frac{{\boldsymbol c}_{\tau_{\rm h}}^{\magenta{i}} (x_i, x_{j}^{\sf d} )}{\varepsilon} \right), \ x = [x_1;  \cdots  ;x_N] ,
\end{align*}
and the initial value \las{$ \alpha_0 \in \bbR_{>0}^N $} is arbitrary.
\hfill $ \diamondsuit $


The important point here is that for the initial value ${ \alpha[k+1,1] }$ of the Sinkhorn iterations at time $ k+1 $, we use the final value $ \alpha[k,{S[k] + 1}] $ of the iterations at the previous time $ k $.
When the difference between $ {\boldsymbol x}[k+1] $ and $ {\boldsymbol x}[k] $ is small, $ \alpha[k,{S[k] + 1}] $ will be a good initial estimate of the scaling variable $ \alpha^*[k+1] $ for $ {\boldsymbol P}^*({\boldsymbol x}[k+1]) $ even if $ S[k] $ is small.
A convenient way to determine $ S[\cdot] $ is to fix it to a suitable constant $ \hat{S}\in \bbZ_{>0} $ in terms of computation time.

The pseudocode of Sinkhorn MPC is described in Algorithm~\ref{algo:sinkhorn_mpc}.
Note that, of course, the proposed method can be applied to systems that are not discretizations of continuous-time systems and originally evolve in discrete time.



\begin{algorithm}
	\caption{Sinkhorn MPC}\label{algo:sinkhorn_mpc}
	\begin{algorithmic}[1]
		\renewcommand{\algorithmicrequire}{\textbf{Input:}}
		\renewcommand{\algorithmicensure}{\textbf{Output:}}
		\Require $ \{x_i^0\}_{i=1}^N, \{x_j^{\sf d}\}_{j=1}^N, \alpha_0, \ \varepsilon, \ \tau_{\rm h}, \ S[\cdot] $\\
		\textit{Initialization}: $ {\boldsymbol x}_i[0] := x_i^0,  \forall i, \ \alpha := \alpha_0$
		\For{$ k = 0, 1, 2,\ldots$}
		\For{$ i,j \in \bbra{N}$ ({\it run in parallel})}
		\State Compute $ {\boldsymbol c}_{\tau_{\rm h}}^i ({\boldsymbol x}_i[k], x_j^{\sf d}) $ defined in \eqref{eq:cost_disc}
		\State $ {\boldsymbol K}_{ij} := \exp\left(- \frac{{\boldsymbol c}_{\tau_{\rm h}}^{\magenta{i}} ({\boldsymbol x}_i [k], x_{j}^{\sf d} )}{\varepsilon} \right) $
			\EndFor
			\For{$ l = 1,2,\ldots,S[k] $}
			\State $ \beta := {\bf 1}_N/N \oslash [{\boldsymbol K}^\top \alpha] $
			\State $ \alpha := {\bf 1}_N/N \oslash [{\boldsymbol K} \beta] $
		\EndFor
		\State $ P: = \alpha^\diagbox {\boldsymbol K} \beta^\diagbox $
		\For{$ i \in \bbra{N} $ ({\it run in parallel})}
			\State Compute $ {\boldsymbol u}_i^{\rm MPC} ({\boldsymbol x}_i[k], x_i^{\rm tmp} (P)) $ and apply it to the agent $ i $
		\EndFor
		\EndFor
	\end{algorithmic}
\end{algorithm}


















\section{Illustrative examples}\label{sec:example}
\subsection{Quadratic cost}\label{subsec:quad_invertible}
This section gives examples for Sinkhorn MPC.
First we consider \blu{a} quadratic cost
\begin{equation}\label{eq:energy}
	{\boldsymbol \ell}_i (x_i, u_i; x_j^{\sf d}) = \|u_i - {\boldsymbol B}_i^{-1} (x_j^{\sf d} - {\boldsymbol A}_i x_j^{\sf d})\|^2 ,
\end{equation}
where we assumed the invertibility of $ {\boldsymbol B}_i $, and let $ {\bbX_i = \bbU_i = \bbR^n }$.
Note that for a constant input $ \bar{{\boldsymbol u}}_{\blu{ij}} := {\boldsymbol B}_i^{-1} (x_j^{\sf d} - {\boldsymbol A}_i x_j^{\sf d}) $, under which $ {\boldsymbol x}_i = {\boldsymbol x}_j^{\sf d} $ is an equilibrium of \eqref{eq:linear_dynamics_disc}, it holds $ {\boldsymbol \ell}_i (x_j^{\sf d}, \bar{{\boldsymbol u}}_{\blu{ij}} ; x_j^{\sf d}) = 0 $.
Then, the dynamics under Sinkhorn MPC can be written as follows \cite[\rev{Section~2.2, pp.~37-39}]{Lewis2012}:
\begin{align}
&{\boldsymbol x}_i[k+1] = \bar{{\boldsymbol A}}_i {\boldsymbol x}_i [k] +  (I - \bar{{\boldsymbol A}}_i) x_i^{\rm tmp} (P[k]), \label{eq:SMPC_linear_x} \\
&{\boldsymbol u}_i^{\rm MPC} (x_i, \red{\hat{x}}) = - {\boldsymbol B}_i^\top ({\boldsymbol A}_i^\top)^{\tau_{\rm h} - 1} {\boldsymbol G}_{i,\tau_{\rm h}}^{-1} {\boldsymbol A}_i^{\tau_{\rm h}} (x_i - \red{\hat{x}} ) \nonumber\\
&\hspace{2.3cm} + {\boldsymbol B}_i^{-1} ( \red{\hat{x}} - {\boldsymbol A}_i \red{\hat{x}}), \ \forall i \in \bbra{N}, \ \red{\forall \subm{x_i}, \hat{x}\in \bbR^n} \nonumber
\end{align}
with \eqref{eq:SMPC_P},~\eqref{eq:SMPC_sink} where
\begin{align*}
&{\boldsymbol K}_{ij}(x)  = \exp \left(- \frac{  \|x_i - x_j^{\sf d}\|_{\calG_{i}}^2 }{\varepsilon}   \right), \\
&\calG_{i} := ({\boldsymbol A}_i^{\tau_{\rm h}})^\top {\boldsymbol G}_{i,\tau_{\rm h}}^{-1}  {\boldsymbol A}_i^{\tau_{\rm h}}, \ {\boldsymbol G}_{i,\tau_{\rm h}} := \sum_{k=0}^{\tau_{\rm h} -1} {\boldsymbol A}_i^k {\boldsymbol B}_i {\boldsymbol B}_i^\top ({\boldsymbol A}_i^\top)^k, \\
&\bar{{\boldsymbol A}}_i := {\boldsymbol A}_i - {\boldsymbol B}_i {\boldsymbol B}_i^\top ({\boldsymbol A}_i^\top)^{\tau_{\rm h} -1} {\boldsymbol G}_{i,\tau_{\rm h}}^{-1} {\boldsymbol A}_i^{\tau_{\rm h}} .
\end{align*}
In the examples below, we use the barycentric target~\eqref{eq:barycentric_target} \subm{as a navigator function}.

First, consider \eqref{eq:linear_dynamics_conti} with
\begin{equation}\label{eq:ex2d_conti}
	A_i = 
	\begin{bmatrix}
		2 & 1.3\\
		-0.5 & 1
	\end{bmatrix}, \
	B_i = I_2, \ \forall i\in \bbra{N} .
\end{equation}
By using the Euler method with a step size $ 0.02 $, we obtain 
\begin{equation}\label{eq:ex2d}
	{\boldsymbol A}_i = 
	\begin{bmatrix}
		1.04 & 0.026\\
		-0.01 & 1.02
	\end{bmatrix}, \
	{\boldsymbol B}_i = 0.02 I_2, \ \forall i\in \bbra{N} .
\end{equation}
Set $ N = 120, \ \varepsilon = 2.0, \ \tau_{\rm h} = 100, \ \alpha_0 = {\bf 1}_N $.
Here, we consider the case where the optimal coupling $ {\boldsymbol P}^* $ is available at each time.
Specifically, rather than using $ S[\cdot] $ fixed beforehand for Sinkhorn MPC, we employ the stopping criterion for the Sinkhorn iterations~\cite[Remark~4.14]{Peyre2019} given by
\begin{equation}\label{ineq:stop}
	\| P[k] {\bf 1}_N - {\bf 1}_N/N \|_1 + \|P[k]^\top {\bf 1}_N - {\bf 1}_N/N \|_1 < 0.005,
\end{equation}
where $ \| \cdot \|_1 $ denotes the $ \ell^1 $-norm.
For given initial and desired states, the trajectories of the agents governed by \eqref{eq:SMPC_linear_x} with \eqref{eq:SMPC_P}--\eqref{eq:SMPC_next},~\subm{\eqref{eq:ex2d}}, \eqref{ineq:stop} are illustrated in Figs.~\ref{fig:2d_trajectory},~\ref{fig:3d_trajectory}.
It can be seen that the agents converge sufficiently close to the target states. We will study the convergence property in Sections~\ref{sec:global},~\ref{sec:linear_system}.
The number of \blu{Sinkhorn} iterations $ \bar{S}[k] $ satisfying \eqref{ineq:stop} at each time $ k $ is shown in Fig.~\ref{fig:iteration}. The number of \blu{iterations} is drastically reduced from $ \bar{S}[0] \simeq 520 $ to $ \bar{S}[1] \simeq 100 $ in one time step, and $ \bar{S} $ continues to decrease as $ k $ increases.
This clarifies that the optimal scaling variables $ (\alpha^*[k],\beta^*[k]) $ can be used for good initial estimates of $ (\alpha^*[k+1], \beta^*[k+1]) $ as expected in Subsection~\ref{subsec:proposed}.

The computation time for one Sinkhorn iteration and the number of \blu{Sinkhorn} iterations~$ \bar{S}[0] $ at the initial time with different $ N $ are shown in Table~\ref{table:computation}. The algorithm has been implemented in MATLAB on MacBook~Pro with Apple~M1~Pro.
Table~\ref{table:computation} also shows the computation time for solving an optimal assignment problem to obtain \subm{the permutation} $ \sigma(\cdot;{\boldsymbol x}[k]) $ by the Hungarian algorithm~\cite{Cao2022}.
As can be seen, the Hungarian algorithm is not scalable and thus not suitable for MPC.
Hence, introducing the entropy regularization to MPC contributes to reducing the computational burden.








\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.45]{figures/2d_trajectory.eps}
	\caption{Trajectories $ {\boldsymbol x}_i[k] = [{\boldsymbol x}_i^{(1)}[k] \ {\boldsymbol x}_i^{(2)} [k]]^\top $ of 120 agents for \eqref{eq:ex2d} (solid), initial states (filled triangles), steady states (filled circles), and desired states (black circles).}
	\label{fig:2d_trajectory}
\end{figure}


\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.35]{figures/3d_trajectory.eps}
	\caption{Trajectories of $ 14 $ agents out of $ 120 $ agents for \eqref{eq:ex2d} (solid), initial states (magenta crosses), and desired states (black circles).}
	\label{fig:3d_trajectory}
\end{figure}



\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.45]{figures/iteration_semilog.eps}
	\caption{Trajectories $ {\boldsymbol x}_i^{(2)} [k] $ of 14 agents (top) and the number of \blu{Sinkhorn} iterations at each time $ \bar{S}[k] $ (bottom, \res{semi-log plot}) for \eqref{eq:ex2d}.}
	\label{fig:iteration}
\end{figure}




\renewcommand{\arraystretch}{1.3}
\begin{table*}[hbtp]
  \caption{Computation time for one Sinkhorn iteration and the Hungarian algorithm, and the number of \blu{iterations} $ \bar{S}[0] $.}
  \label{table:computation}
  \centering
  \begin{tabular}{lrrrr}
    \hline
      & \begin{tabular}{c}Computation time for\\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}} Sinkhorn iteration \end{tabular}  & \begin{tabular}{c}$ \bar{S}[0] $ \\with $ \varepsilon = 2.0 $ \end{tabular}  &\begin{tabular}{c} $ \bar{S}[0] $\\with $ \varepsilon = 4.0 $ \end{tabular} & Hungarian algorithm~\cite{Cao2022}\\
    \hline 
    $ N  =120 $  & $ 0.0060 $~ms  & $ 527 $  & $ 267 $ & $ 0.1 $~s\\
    $ N = 500 $  & $ 0.21 $~ms   & $ 381 $ & $ 189 $ & $ 0.6 $~s\\
		$ N = 1000 $  & $ 0.83 $~ms   & $ 379 $ & $ 186 $ & $ 3.0 $~s\\
		$ N = 3000 $  & $ 3.0 $~ms   & $ 256 $ &  $ 129 $&  $ 57.0 $~s\\
    \hline
  \end{tabular}
\end{table*}
\renewcommand{\arraystretch}{1.0}





\subsection{Effect of the number of \blu{Sinkhorn} iterations and the regularization parameter on Sinkhorn MPC}
Here, we describe how the number of \blu{Sinkhorn} iterations $ S $ affects the bahavior of Sinkhorn MPC.
To this end, consider a simple case for~\eqref{eq:energy} with $ N = 14, \ \tau_{\rm h} = 20, \ \varepsilon = 0.1 $, and
\begin{equation}\label{eq:ex1d}
	{\boldsymbol A}_i = 1, \ {\boldsymbol B}_i = 0.1, \ \forall i \in \bbra{N}.
\end{equation}
Then the trajectories of the agents with $ S[k] = 1, \forall k $ and $ S[k] = 5, \forall k $ are illustrated in Fig.~\ref{fig:iteration_compare}. 
Also, the trajectories without the regularization $ (\varepsilon = 0) $ following the discretized version of \eqref{eq:permutation_dynamics} are shown with the black dotted lines.
In this example, for all the cases, $ {\boldsymbol x}[k] $ converges to almost the same point close enough to the desired distribution.
The total energy cost $ \sum_{i,k} \|{\boldsymbol u}_i[k] \|^2 $ for the case without the regularization is $ 21.6 $.
It can be seen from Fig.~\ref{fig:iteration_compare} that one iteration per time step is not enough to determine an appropriate destination for each agent while performing control and results in the total energy cost $ 34.7 $.
On the other hand, the trajectories for five iterations are similar to the trajectories without the regularization, and the total energy cost is $ 19.1 $. Note that since we use MPC, the total cost for Sinkhorn MPC can be smaller than for the case without the regularization as in this example.



\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.35]{figures/iteration_compare.eps}
	\caption{Trajectories $ {\boldsymbol x}_i [k] $ of 14 agents for \eqref{eq:ex1d} with $ S[k] \equiv 1 $ (chain) and $ S[k] \equiv 5 $ (solid), respectively, and desired states (black circles). The black dotted lines indicate trajectories without \blu{the} regularization ($ \varepsilon = 0 $) following \eqref{eq:permutation_dynamics}.}
	\label{fig:iteration_compare}
\end{figure}





Next, we investigate the effect of the regularization parameter $ \varepsilon $ on Sinkhorn MPC.
We continue to consider \eqref{eq:ex1d} and set $ N = 7, S[k] = 1,  \forall k $.
Then the trajectories of the agents with $ \varepsilon = 0.4, 0.9 $ are shown in Fig.~\ref{fig:ex1d}.
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 fact, it is known that the convergence of the Sinkhorn algorithm deteriorates as $ \varepsilon \rightarrow + 0 $~\cite[Remark~4.14]{Peyre2019}; see Table~\ref{table:computation}. This degrades the transient behaviors under Sinkhorn MPC.
Taking larger $ S[k] $ remedies this issue.
The steady states $ \lim_{k\rightarrow \infty} {\boldsymbol x}_i [k] $ for a fixed set of initial states as a function of $ \varepsilon $ are illustrated in Fig.~\ref{fig:steady_state}. Note that for different initial states, we obtained the same result or the one flipped upside down.
The obtained behavior is due to the fact that as $ \varepsilon $ becomes larger, the optimal coupling of the entropic OT problem is more blurred to the maximum entropy coupling $ {\bf 1}_N {\bf 1}_N^\top /N^2 $~\cite[Proposition~4.1]{Peyre2019}.
The behavior of the equilibrium points for Sinkhorn MPC as $ \varepsilon \rightarrow + 0 $ will be revealed in Lemma~\ref{lemma:exponential_convergence} in Section~\ref{sec:linear_system}.
Although we have considered the simple setting~\eqref{eq:ex1d}, the above observations apply to the general case.


\begin{figure}[t]
	\begin{minipage}[b]{0.53\linewidth}
	
		\hspace{-0.45cm}
		\includegraphics[keepaspectratio, scale=0.3]
		{figures/epsilon_compare.eps}
		\subcaption{}\label{fig:ex1d}
	\end{minipage}
	\begin{minipage}[b]{0.46\linewidth}
	
	
		\includegraphics[keepaspectratio, scale=0.3]
		{figures/steady_state.eps}
		\subcaption{}\label{fig:steady_state}
	\end{minipage}


	\caption{(a) Trajectories $ {\boldsymbol x}_i[k] $ of 7 agents for \eqref{eq:ex1d} with $ \varepsilon = 0.9 $ (chain) and $ \varepsilon = 0.4 $ (solid), respectively, \rev{trajectories without the regularization (dotted)}, and desired states (black circles).  \\
	\rev{\subm{(b)} Semi-\blu{log} plot of steady states $ \lim_{k\rightarrow\infty} {\boldsymbol x}_i [k] $ with different $ \varepsilon \in [10^{-1},10^2] $ (solid) and desired states (dashed).}
\end{figure}


\subsection{Non-quadratic cost}
Lastly, we investigate the behavior of Sinkhorn MPC for \las{the} non-quadratic cost and \las{the} bounded input spaces:
\begin{align}
	&{\boldsymbol \ell}_i (x_i, u_i; x_j^{\sf d}) = \|x_i - x_j^{\sf d} \|^2  +  \|u_i - {\boldsymbol B}_i^{-1} (x_j^{\sf d} - {\boldsymbol A}_i x_j^{\sf d})\|_1 \subm{,}  \label{eq:l1_cost} \\
	&\bbU_i =  \subm{\left\{ u_i = [u_i^{(1)} \ u_i^{(2)}]^\top \in \bbR^{2} : |u_i^{(j)}| \le 5,\  j = 1,2 \right\} } , \nonumber\\
	&\subm{\bbX_i = \bbR^2}, \hspace{5cm} \forall i\in \bbra{N} .\nonumber
\end{align}
Here $ \ell^1 $-norm encourages $ u_i $ to be identically equal to $ {\boldsymbol B}_i^{-1} (x_j^{\sf d} - {\boldsymbol A}_i x_j^{\sf d}) $ for long time and is used for sparse optimal control~\cite{Nagahara2016,Ito2021sparse}.
For the computation of $ {\boldsymbol c}_{\tau_{\rm h}} $ and $ {\boldsymbol u}_i^{\rm MPC} $, we used {\tt cvx} package with MATLAB~\cite{cvx}.
Set $ N  = 15, \varepsilon = 2.0, \tau_{\rm h} = 50, S [k] = 20, \forall k $ and consider \eqref{eq:ex2d}.
Then, Fig.~\ref{fig:lasso} shows the trajectories under Sinkhorn MPC with the barycentric target~\eqref{eq:barycentric_target}. As can be seen, the agents have achieved the desired distribution.
This result shows that Sinkhorn MPC achieves the transport also for the non-quadratic cost.




\begin{figure}[tb]
	\centering
	\includegraphics[scale=0.35]{figures/lasso.eps}
	\caption{Trajectories of $ 15 $ agents for \eqref{eq:l1_cost} with $ \varepsilon = 2.0 $ (solid) and desired states (black circles).}
	\label{fig:lasso}
\end{figure}









\section{Global convergence property of Sinkhorn MPC}\label{sec:global}
In the remainder of this paper, we investigate the fundamental properties of Sinkhorn MPC.
In this section, we consider the case where the sampling time for obtaining \eqref{eq:linear_dynamics_disc} is small and $ S $ is large enough so that the dynamics~\eqref{eq:mpc_dynamics_conti} is well approximated by the discretized system~\eqref{eq:SMPC_x}.
Then, we analyze the global convergence property of the dynamics~\eqref{eq:mpc_dynamics_conti}.
To this end, we suppose the following condition holds.
\begin{assumption}\label{ass:positive_cost}
  For all $ i\in \bbra{N} $, $ \bbX_i = \bbU_i = \bbR^n $, $ B_i $ is invertible, and
	the function $ \ell_i $ satisfies
  \begin{equation}\label{eq:positive_cost}
    \ell_i (x_i, u_i; x_j)
    \begin{cases}
      =0  & \text{if} \ x_i =x_j \ \text{and} \ u_i = -B_i^{-1} A_i x_i, \\
      >0  & \text{otherwise}.
    \end{cases}
  \end{equation}
  \hfill $ \diamondsuit $
\end{assumption}
Then, a point $ x^{\sf e} = [x_1^{\sf e};\cdots;x_N^{\sf e}] \in \bbR^{nN} $ satisfying
\begin{equation}\label{eq:equi_general}
  x_i^{\sf e} = x^{\rm tmp}_{i} (P^*(x^{\sf e})), \ \forall i \in \bbra{N}
\end{equation}
is an equilibrium of \eqref{eq:mpc_dynamics_conti}.
Indeed, by \eqref{eq:positive_cost} and \eqref{eq:equi_general}, the constant input $ u_i (t) \equiv -B_i^{-1} A_i x_i^{\sf e} $, under which the state process starting from $ x_i(0) = x_i^{\sf e} $ is the constant $ x_i(t) \equiv x_i^{\sf e} $, is the unique optimal solution of the OC problem \eqref{eq:cost_finite_conti} with $ \check{x}_i = x_i^{\sf e}, x_j^{\sf d}  = x^{\rm tmp}_{i} (P^*(x^{\sf e})) = x_i^{\sf e} $. 
Therefore,
\[
  u_{i}^{\rm MPC} \left(x_i^{\sf e}, x_i^{\rm tmp} (P^*(x^{\sf e})) \right) = -B_i^{-1} A_i x_i^{\sf e} ,
\]
and
\[
  A_i x_i^{\sf e} + B_i u_i^{\rm MPC} \left(x_i^{\sf e}, x_i^{\rm tmp} (P^*(x^{\sf e}))\right) = 0 .
\]


The following proposition ensures the existence of a point satisfying \eqref{eq:equi_general}.
\begin{proposition}\label{prop:fixed_point_general}
  Assume that the codomain $ \calX $ of $ x_i^{\rm tmp} $ is a compact convex set and for all $ i\in \bbra{N} $, $ x_i^{\rm tmp} : \bbR_{\ge 0}^{N\times N} \rightarrow \calX $ is continuous.
	Assume further that for all $ i,j\in \bbra{N} $, $ c_{T_{\rm h}}^i (x_i,x_j^{\sf d}) $ is continuous in $ x_i $.
  Then, the set ${ \calR := \{ x^{\sf e} \in \bbR^{nN} : x_i^{\sf e} =  x^{\rm tmp}_{i} (P^*(x^{\sf e})), \ \forall i \in \bbra{N}\} }$ is non-empty.
\end{proposition}
\begin{pf}
	Define a map $ h : \bbR^{nN} \rightarrow \bbR^{nN} $ as
	\begin{equation}
	h(x) := \left[x_1^{\rm tmp} (P^*(x)); \cdots ; x_N^{\rm tmp} (P^*(x))\right]
	, \ x \in \bbR^{nN} . \label{eq:func_fixed_point}
	\end{equation}
	It is obvious from \eqref{eq:opt_form},~\eqref{eq:opt_scaling_condition}, and the continuity of $ c_{T_{\rm h}}^i $ that $ P^* $ is continuous.
  Since $ x_i^{\rm tmp} $ and $ P^* $ are continuous, $ h $ is also continuous.
	The set of all fixed points of $ h $ coincides with $ \calR $.
	For brevity, we abuse notation and regard $ \calX^N $ as a subset of $ \bbR^{nN} $.
	Let $ h_{\calX} : {\calX}^N \rightarrow {\calX}^N $ be the restriction of $ h $ to $ {\calX}^N $.
	Now we use Brouwer's fixed point theorem (\comm{see e.g.,~\cite[Corollary~1.1.1]{Florenzano2003}}). That is, since $ h_{\calX} $ is a continuous map from a compact convex set $ \calX^N $ into itself, there exists a point $ x^{{\sf e}} \in {\calX}^N $ such that $ x^{{\sf e}} = h(x^{\sf e}) $.
  \hfill $ \Box $
\end{pf}



Next, as a tool for the convergence analysis of \eqref{eq:mpc_dynamics_conti}, we consider the entropic OT cost
\begin{equation}\label{eq:transport_cost_conti}
	\calE(x,x^{\sf d}) := \min_{P\in \calT(\frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N}, \frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N})} \ \sum_{i,j\in \bbra{N}} C_{ij} (x) P_{ij} - \varepsilon {\rm H}(P) \subm{,}
\end{equation}
\subm{where $ x^{\sf d} := [x_1^{\sf d}; \cdots ; x_N^{\sf d}] $.}
Assume that for any $ i,j\in \bbra{N} $, $ c_{T_{\rm h}}^i (x_i, x_j^{\sf d}) $ is continuously differentiable with respect to $ x_i $.
Then, thanks to the regularization, $ \calE(x,x^{\sf d}) $ is continuously differentiable~\cite[Eq. (9.6)]{Peyre2019} with respect to $ x $\subm{,} and
\begin{align}
	\nabla_{x_i} \calE(x, x^{\sf d}) = \sum_{j=1}^N P_{ij}^* (x) \nabla_{x_i} c_{T_{\rm h}}^i (x_i, x_j^{\sf d}) . \label{eq:grad_E}
\end{align}
This is in clear contrast to the case without the entropy regularization ($ \varepsilon = 0 $), in which the optimal coupling $ P^*(x) $ is not continuous similarly to the optimal permutation $ \sigma(\cdot;x) $, and thus $ \calE $ with $ \varepsilon = 0 $ is not differentiable. This difference is crucial for analyzing the \subm{global convergence property} of \eqref{eq:mpc_dynamics_conti} as shown in Theorem~\ref{thm:lasalle}.
If \comm{navigator functions} $ \{x_i^{\rm tmp}\} $ are designed appropriately, then it is expected that the state $ x(t) $ following~\eqref{eq:mpc_dynamics_conti} moves in a direction where the cost $ \calE(x(t),x^{\sf d}) $ decreases.
In fact, the following result shows that this is indeed the case and, as a result, ensures the convergence to the set of equilibria.
We say that $ x(t) $ converges to a set $ \calM \subset \bbR^{nN} $ as $ t\rightarrow\infty $ if for each $ \epsilon > 0 $, there exists $ \tau > 0 $ such that $\inf_{p\in \calM} \| x(t) - p \| < \epsilon $ for all $t \ge \tau $.
The proof of Theorem~\ref{thm:lasalle} is shown in Appendix~\ref{app:lasalle}.

\begin{theorem}\label{thm:lasalle}
  Suppose that Assumption~\ref{ass:positive_cost} holds.
  Assume that for any $ i\in \bbra{N} $ and $ \hat{x}\in \bbR^n $, $ c_{T_{\rm h}}^i (x_i, \hat{x}) $ is continuously differentiable with respect to $ x_i $ and $ T_{\rm h} $. Also assume that for any $ i,j\in \bbra{N} $,
  \begin{equation}\label{eq:unbounded}
    c_{T_{\rm h}}^i (x_i, x_j^{\sf d}) \rightarrow +\infty ~~ \text{as} \ \|x_i\| \rightarrow +\infty .
  \end{equation}
  Assume further that for any $ i\in \bbra{N} $, there exists a constant $ a_i > 0 $ such that for any $ x = [x_1;\cdots;x_N] \in \bbR^{nN} $,
	\begin{equation}\label{eq:tmp_condition}
		\sum_{j=1}^{N} P_{ij}^* (x) \nabla_{1} c_{T_{\rm h}}^i (x_i, x_j^{\sf d}) = a_i \nabla_{1} c_{T_{\rm h}}^i \left(x_i, x_i^{\rm tmp} (P^* (x))\right) ,
	\end{equation}
	where $ \nabla_1 c_{T_{\rm h}}^i $ denotes the gradient of $ c_{T_{\rm h}}^i $ with respect to the first variable.
	Then, for any initial state $ x(0) = x^0 \in \bbR^{nN} $, \blu{the solution} $ x(t) $ \blu{of \eqref{eq:mpc_dynamics_conti}} converges to the set $ \calR = \left\{x^{\sf e} \in \bbR^{nN} : x_i^{\sf e}  = x_i^{\rm tmp} (P^* (x^{\sf e})), \forall i \in \bbra{N} \right\} $ as $ t\rightarrow \infty $.
	\hfill $ \diamondsuit $
\end{theorem}






\begin{rmk}
  In the proof of Theorem~\ref{thm:lasalle}, the linearity of the system~\eqref{eq:linear_dynamics_conti} is not used at all.
  Hence, the same proof works for nonlinear systems of the form:
  \[
    \dot{x}_i = f_i (x_i) + g_i(x_i) u_i ,
  \]
  where $ g_i(x_i) $ is square and invertible for any $ x_i\in \bbR^n $.
  \hfill $ \diamondsuit $
\end{rmk}
\begin{rmk}
  Assume that $\calR $ consists of only isolated points. Then by Theorem~\ref{thm:lasalle}, for any initial state, $ x(t) $ converges to one of the equilibrium points in $ \calR $ as $ t\rightarrow \infty $.
  \hfill $ \diamondsuit $
\end{rmk}
\begin{rmk}\label{rmk:nonuniform}
	Let us consider the general case where the number of agents $ N $ is not equal to the number of targets $ M $, and the agents and targets have mass distributions that are not necessarily uniform, i.e., $ {\mathbf a} \neq {\bf 1}_N/N, \ {\mathbf b} \neq {\bf 1}_M/M $. Then, $ P^* $ in the dynamics~\eqref{eq:mpc_dynamics_conti} is replaced by
	\begin{equation}
	\bar{P}^* (x) := \mathop{\rm arg~min}\limits_{P\in \calT({\mathbf a}, {\mathbf b})} \ \sum_{i\in \bbra{N}, j\in \bbra{M}} C_{ij} (x) P_{ij} - \varepsilon {\rm H}(P) .
	\end{equation}
	Even in this case, the same proof as in Theorem~\ref{thm:lasalle} works. That is, under the assumptions in Theorem~\ref{thm:lasalle}, for any initial state $ x(0) = x^0 \in \bbR^{nN} $, $ x(t) $ converges to $ \calR $ as $ t\rightarrow \infty $.
	In other words, Sinkhorn MPC can be applied to general OT problems whereas MPC with e.g., the Hungarian algorithm works only for OT problems that are equivalent to assignment problems. This is one of the advantages of the proposed method.
	\hfill $ \diamondsuit $
\end{rmk}

The condition \eqref{eq:tmp_condition} gives a guideline for the design of $ x_i^{\rm tmp} $.
However, it is not necessary for $ x_i^{\rm tmp} $ to satisfy \eqref{eq:tmp_condition} to ensure the convergence. In fact, in Section~\ref{sec:example}, we observed that for the non-quadratic cost~\eqref{eq:l1_cost}, the barycentric projection~\eqref{eq:barycentric_target}, which is not desinged based on \eqref{eq:tmp_condition} achieves the transport to the target distribution.



Next, as a specific example, we consider a quadratic cost
\begin{equation}\label{eq:quadratic}
	\ell_i (x_i, u_i; x_j^{\sf d}) = \|u_i + B_i^{-1} A_i x_j^{\sf d}\|^2 .
\end{equation}
Then the corresponding control law $ u_i^{\rm MPC} $ and the cost function $ c_{T_{\rm h}}^i $ can be written as follows~\cite[Section~3.3,~pp.~138-140]{Lewis2012}:
\begin{align}
  &u_i^{\rm MPC} (x_i, \red{\hat{x}}) = - B_i^\top \scrG_i (x_i - \red{\hat{x}} ) - B_i^{-1} A_i \hat{x}, \label{eq:mpc_quad_conti}\\
    &c_{T_{\rm h}}^i (x_i, \hat{x}) = \| x_i - \hat{x} \|_{\scrG_i}^2, ~~ \forall x_i, \hat{x} \in \bbR^n,\ \forall i \in \bbra{N}, \label{eq:quad_value_func_conti}
  \end{align}
  where
  \begin{align*}
  & \scrG_i := \Ee^{A_i^\top T_{\rm h}} G_{i,T_{\rm h}}^{-1} \Ee^{A_i T_{\rm h}}, \ G_{i,T_{\rm h}} := \int_0^{T_{\rm h}} \Ee^{A_i t} B_i B_i^\top  \Ee^{A_i^\top t} {\rm d} t .
  \end{align*}
  Thus, the condition \eqref{eq:tmp_condition} is rewritten as
  \begin{align*}
    \frac{2}{N} \scrG_i \biggl(x_i - N\sum_j P_{ij}^* (x) x_j^{\sf d} \biggr) = 2a_i \scrG_i \left(x_i - x_i^{\rm tmp} (P^*(x))\right) .
  \end{align*}
  This is fulfilled by $ a_i = 1/N $ and the barycentric projection~\eqref{eq:barycentric_target}.

The cost function \eqref{eq:quadratic} does not satisfy Assumption~\ref{ass:positive_cost} because it does not depend on the variable $ x_i $. Nevertheless, similarly to Theorem~\ref{thm:lasalle}, the following convergence result holds. The proof is given in Appendix~\ref{app:lasalle}.
\begin{corollary}\label{cor:convergence}
  Consider the quadratic cost \eqref{eq:quadratic} with $ {\bbX_i = \bbU_i = \bbR^n} $ for all $ i\in \bbra{N} $ and the barycentric target~\eqref{eq:barycentric_target}.
  Then, for any initial state $ x(0) = x^0 \in \bbR^{nN} $, \blu{the solution} $ x(t) $ \blu{of \eqref{eq:mpc_dynamics_conti}} converges to the set \blu{$ \calR = \{x^{\sf e} \in \bbR^{nN} : x_i^{\sf e}  = N \sum_{j=1}^N P_{ij}^*(x^{\sf e}) x_j^{\sf d}  , \forall i \in \bbra{N}\} $} as $ t\rightarrow \infty $.
	\hfill $ \diamondsuit $
\end{corollary}



The above result justifies that the barycentric projection~\eqref{eq:barycentric_target}, which \subm{is} typically used to approximate Monge's OT map from a coupling matrix, gives an appropriate direction where the cost $ \calE(x(t), x^{\sf d}) $ becomes smaller under Sinkhorn MPC for the quadratic cost~\eqref{eq:quadratic}.




\res{Theorem~\ref{thm:lasalle} and Corollary~\ref{cor:convergence} ensure the global convergence of the proposed method that uses the optimal coupling $ P^*(x(t)) $ at each time.
Hence, in terms of the convergence, it is desirable to perform a sufficiently large number of Sinkhorn iterations at each time to obtain a coupling close enough to $ P^*(x(t)) $.
On the other hand, as observed in Fig.~\ref{fig:iteration}, once we perform a sufficiently large number of iterations at some point, we can obtain a coupling close enough to $ P^*(x(t)) $ at later times by a smaller number of iterations. 
This implies that we can reduce the computational burden at later times while maintaining the convergence property.}









\section{Boundedness and local asymptotic stability for Sinkhorn MPC with a quadratic cost}\label{sec:linear_system}
In this section, we consider the general case where the number of the Sinkhorn iterations $ S $ is not necessarily large, and thus $ P[k] $ may not be close enough to the optimal coupling $ P^* ({\boldsymbol x}[k]) $.
Then, we elucidate the fundamental properties of Sinkhorn MPC on discrete-time systems~\subm{\eqref{eq:linear_dynamics_disc}}.
Specifically, we reveal that even when $ S $ is small, the ultimate boundedness and the local asymptotic stability for Sinkhorn MPC hold for \modi{the quadratic cost \eqref{eq:energy}} and $ \bbX_i = \bbU_i = \bbR^n $. \rev{Hereafter, we assume the invertibility of $ {\boldsymbol B}_i $.}
In addition, for notational simplicity, we deal only with the case where just one Sinkhorn iteration is performed at each time, i.e., $ S[k] = 1, \forall k $. Nevertheless, by similar argument, all \rev{of} the results in this section are still valid when more iterations are performed.
For convenience, we recall the dynamics under Sinkhorn MPC with $ S \subm{[k]}  \equiv 1 $:
\begin{align}
	&{\boldsymbol x}_i[k+1] = \bar{{\boldsymbol A}}_i {\boldsymbol x}_i [k] +  (I - \bar{{\boldsymbol A}}_i) x_i^{\rm tmp} (P[k]),\ \forall i \in \bbra{N}, \label{eq:SMPC_linear_x_1} \\
	&P[k] = \alpha[k+1]^\diagbox {\boldsymbol K}({\boldsymbol x}[k]) \beta[k]^\diagbox, \label{eq:SMPC_P_1} \\
	&\alpha[k+1] = {\bf 1}_N/N \oslash \left[{\boldsymbol K}({\boldsymbol x}[k])\beta[k] \right], \label{eq:SMPC_alpha_1} \\
	&\beta[k] = {\bf 1}_N/N \oslash \left[{\boldsymbol K}({\boldsymbol x}[k])^\top \alpha[k] \right], \label{eq:SMPC_beta_1}   \\
	&{\boldsymbol x}_i [0] = x_i^0, \ \alpha[0] = \alpha_0 . \nonumber 
\end{align}



\subsection{Ultimate boundedness for Sinkhorn MPC}\label{subsec:bounded}
Here, we assume that for the codomain $ \calX $ of $ x_i^{\rm tmp} $, there exists a constant $ \bar{r} > 0 $ such that
\begin{equation}\label{ineq:x_upp}
	\|x \| \le  \bar{r} , \ \forall x \in \calX .
\end{equation}
For example, if $ \calX $ is the convex hull of $ \{x_j^{\sf d}\}_j $, we can take $ \bar{r} = \max_j \|x_j^{\sf d}\| $.
It is known that, under the assumption that $ {\boldsymbol B}_i $ is invertible, $ \bar{{\boldsymbol A}}_i $ is stable, i.e., \blue{the spectral radius $ \rho_i $ of $ \bar{{\boldsymbol A}}_i $ satisfies $ \rho_i < 1 $} \cite[Corollary~1]{Kwon1975}.
Using this fact, we derive the ultimate boundedness of \eqref{eq:SMPC_linear_x_1} with \eqref{eq:SMPC_P_1}--\eqref{eq:SMPC_beta_1}.

\subm{
\begin{proposition}\label{prop:bounded}
	Assume that there exists a constant $ \bar{r} > 0 $ satisfying \eqref{ineq:x_upp}.
	Then, for any $ \delta > 0, \{x_i^0\}_i $, and $ \{\nu_i\}_i $ satisfying $ \nu_i > 0, \rho_i + \nu_i < 1, \forall i\in \bbra{N} $, there exist $ \kappa_i(\nu_i) > 0, i\in \bbra{N} $ and $ \tau(\delta, \{x_i^0\},\{\nu_i \} ) \in \bbZ_{>0} $ such that the solution $ \{{\boldsymbol x}_i\}_i $ of \eqref{eq:SMPC_linear_x_1} with \eqref{eq:SMPC_P_1}--\eqref{eq:SMPC_beta_1} satisfies
	\begin{equation}\label{ineq:ultimate_bound}
		\|{\boldsymbol x}_i [k] \| < \delta  +  \frac{ \kappa_i \bar{r} \|I- \bar{{\boldsymbol A}}_i\|_2}{1 - (\rho_i + \nu_i ) }  , \ \forall k \ge \blue{\tau}, \ \forall i \in \bbra{N}.
	\end{equation}
\end{proposition}
\begin{pf}
	Let $ \tilde{{\boldsymbol u}}_i [k] := ( I - \bar{{\boldsymbol A}}_i) x_i^{\rm tmp} (P[k]) $. 
	Then, it follows from \eqref{ineq:x_upp} that
	\begin{align*}
	\| \tilde{{\boldsymbol u}}_i [k] \| \le \bar{r} \| I -  \bar{{\boldsymbol A}}_i \|_2, \ \forall k\in \bbZ_{\ge 0} .
	\end{align*}
	By \cite[Lemma~7.3.2]{Golub2013}, for any $ \nu_i > 0 $, there exists $ \kappa_i (\nu_i) > 0 $ such that
	\[
		\| \bar{{\boldsymbol A}}_i^k \|_{\subm{2}} \le \kappa_i (\rho_i + \nu_i )^k, \ \forall k \in \bbZ_{\ge 0} . 
	\]
	Hence, the desired result is straightforward from
	\begin{align*}
	\| {\boldsymbol x}_i [k] \| \le \| \bar{{\boldsymbol A}}_i^k \|_2 \| x_i^0 \| + \sum_{s = 1}^{k} \| \bar{{\boldsymbol A}}_i^{s-1} \|_2 \| \tilde{{\boldsymbol u}}_i [k-s]\| .
	\end{align*}
	\hfill $ \Box $
\end{pf}
}



\rev{We emphasize that Proposition~\ref{prop:bounded} holds for any \comm{navigator function} $ x_i^{\rm tmp} $ \fin{whose codomain $ \calX $ satisfies \eqref{ineq:x_upp}.}}



\subsection{Existence of the equilibrium points}\label{subsec:equilibrium}
In the \rev{remainder} of this section, we focus on the barycentric target~\eqref{eq:barycentric_target}.
For $ (x,\beta) \in \bbR^{nN} \times \bbR_{>0}^N $ and $ X^{\sf d} := [x_1^{\sf d} \ \cdots \ x_N^{\sf d}] \in \bbR^{n\times N} $, define
\begin{align*}
&f_1 (x,\beta) := \{\bar{{\boldsymbol A}}_i\}_i^\diagbox x \nonumber\\
&\hspace{1.6cm}+ N  \{I_n - \bar{{\boldsymbol A}}_i\}_i^\diagbox  (X^{\sf d})^{\diagbox,N} {\rm vec}(\tilde{P}(x,\beta)), \\
&f_2 (x,\beta) := {\bf 1}_N/N \oslash \left[ {\boldsymbol K}(f_1 (x,\beta))^\top ({\bf 1}_N/N \oslash [{\boldsymbol K}(x) \beta]) \right], \\
&\tilde{P}(x,\beta) :=  \left( {\bf 1}_N/N \oslash [{\boldsymbol K}(x)\beta]  \right)^\diagbox {\boldsymbol K}(x)  \beta^\diagbox. 
\end{align*}
Then, the collective dynamics \eqref{eq:SMPC_linear_x_1}--\eqref{eq:SMPC_beta_1} is \blu{rewritten as}
\begin{align}
&{\boldsymbol x}[k+1] = f_1 ({\boldsymbol x}[k], \beta[k])  , \label{eq:x_global}\\
&\beta[k+1] = f_2({\boldsymbol x}[k], \beta[k]) . \label{eq:beta_global}
\end{align}




A point $ x^{{\sf e}} = \rev{[x_1^{\sf e}; \ \cdots \ ;x_N^{\sf e} ]} \in \bbR^{nN} $ is an equilibrium of \eqref{eq:x_global},~\eqref{eq:beta_global} if and only if
\begin{align}
&\hspace{-0.2cm}(I_n - \bar{{\boldsymbol A}}_i) \biggl( x_{i}^{{\sf e}} - N \sum_{j=1}^{N} {\boldsymbol P}_{ij}^* (x^{\sf e}) x_j^{\sf d} \biggr) = 0 , \ \forall i \in \bbra{N}. \nonumber
\end{align}
Here note that $ {\boldsymbol P}^* $ satisfies
\begin{align}
&\hspace{-0.2cm} {\boldsymbol P}_{ij}^* (x) =  \alpha_{i}^* {\boldsymbol K}_{ij}(x) \beta_{j}^*,\ \alpha^*, \beta^* \in \bbR_{>0}^N,\\
&\hspace{-0.2cm}\alpha^* = {\bf 1}_N/N \oslash \left[ {\boldsymbol K}(x) \beta^*	\right],  \beta^* = {\bf 1}_N/N \oslash \left[ {\boldsymbol K}(x)^\top \alpha^*	\right]. \label{eq:optimal_beta}
\end{align}
The stability of $ \bar{{\boldsymbol A}}_i $ implies that it has no eigenvalue \rev{equal to} $ 1 $, and therefore $ I_n - \bar{{\boldsymbol A}}_i $ is invertible. Thus, the necessary and sufficient condition for the equilibria is given by
\begin{equation}\label{eq:equilibrium_iff}
x_{i}^{{\sf e}} - N\sum_{j=1}^{N} {\boldsymbol P}_{ij}^* (x^{\sf e}) x_j^{{\sf d}} = 0, \ \forall i\in \bbra{N} ,
\end{equation}
which coincides with \eqref{eq:equi_general}.
Similarly to Proposition~\ref{prop:fixed_point_general}, we show the existence of an equilibrium.
\begin{corollary}\label{prop:fixed_point}
	The dynamics \eqref{eq:x_global}, \eqref{eq:beta_global} has at least one equilibrium point $ (x^{\sf e}, \beta^{\sf e}) \in \bbR^{nN} \times (\bbR_{> 0}^N / {\sim}) $.
\end{corollary}
\begin{pf}
	Note that if a point $ x^{{\sf e}}\in \bbR^{nN} $ satisfies \eqref{eq:equilibrium_iff}, the corresponding $ \beta^{\sf e} \in \bbR_{> 0}^N / {\sim} $ is uniquely determined by $ {\beta^{\sf e} = \beta^* }$ in \eqref{eq:optimal_beta} with $ x= x^{\sf e} $~\cite[Theorem~4.2]{Peyre2019}.
	Note also that for any $ i \in \bbra{N}$ and any $ x \in \bbR^{nN} $, $ N\sum_{j=1}^N {\boldsymbol P}_{ij}^* (x) x_j^{\sf d} $ belongs to the convex hull $ \calX $ of $ \{x_j^{\sf d}\}_j $.
	Then, by the same argument as in the proof of Proposition~\ref{prop:fixed_point_general}, we obtain the desired result.
	\hfill $ \Box $
\end{pf}

Sometimes, in order to emphasize the dependence of $ (x^{\sf e}, \beta^{\sf e}) $ on $ \varepsilon $, we write $ (x^{\sf e} (\varepsilon), \beta^{\sf e} (\varepsilon)) $.







\subsection{Local asymptotic stability for Sinkhorn MPC}
Next, we analyze the stability of the equilibrium points.
For this purpose, the following lemma is crucial when $ \varepsilon $ is small.
The proof is shown in Appendix~\ref{app:convergence}.
\begin{lemma}\label{lemma:exponential_convergence}
	Assume that $ x_i^{\sf d} \neq x_j^{\sf d} $ for all $(i,j), \ i\neq j $, \blu{and $ {\boldsymbol A}_i $ is invertible for all $ i\in \bbra{N} $.}
	For a permutation $ \sigma : \bbra{N} \rightarrow \bbra{N} $, define $ x^{{\sf d}} (\sigma) := \rev{[x_{\sigma(1)}^{{\sf d}}; \ \cdots \ ;x_{\sigma(N)}^{{\sf d}} ]} $ and a permutation matrix $ P^\sigma = (P_{ij}^\sigma)$ as $ P_{ij}^\sigma := 1/N $ if $ j= \sigma(i) $, and $ 0 $, otherwise.
	Then \subm{for any permutation $ \sigma $,} there exists an equilibrium $ (x^{\sf e}(\varepsilon), \beta^{\sf e}(\varepsilon)) $ of \eqref{eq:x_global}, \eqref{eq:beta_global} such that $ x^{\sf e} (\varepsilon) $ and $ P^* (x^{\sf e} (\varepsilon)) $ converge exponentially to $ x^{{\sf d}} (\sigma) $ and $ P^\sigma $, respectively, as $ \varepsilon \rightarrow +0 $, i.e., there exists $ \zeta > 0 $ such that
	\[
		\lim_{\varepsilon \rightarrow +0} \frac{\|\eta(\varepsilon)  \|_2}{\exp (-\zeta/\varepsilon)} = 0
	\]
	for $ \eta (\varepsilon) = x^{\sf e} (\varepsilon) - x^{\sf d} (\sigma) $ and $ \eta (\varepsilon) = P^*(x^{\sf e} (\varepsilon)) - P^\sigma$.
	\hfill $ \diamondsuit $
\end{lemma}

Denote by $ {\rm Exp}(\sigma) $ the set of all equilibria \fin{$ (x^{\sf e} (\cdot), \beta^{\sf e} (\cdot)) $} of \eqref{eq:x_global}, \eqref{eq:beta_global} having the property in Lemma~\ref{lemma:exponential_convergence} for a permutation $ \sigma $.


For $ \bar{P} \in \bbR^{N\times N} $ and $ x  = \rev{[x_1; \ \cdots \ ;x_N]} \in \bbR^{nN} $, define
\begin{align*}
	V_{\subm{\bar{P}}}(x) := \sum_{i=1}^{N} \Bigl\| x_i - N \sum_{j=1}^{N} \bar{P}_{ ij} x_j^{\sf d} \Bigr\|_{\calG_i}^2 .  
\end{align*}
Then, $ V_{\subm{\bar{P}}} $ is a Lyapunov function of \eqref{eq:x_global} where $ \tilde{P}(x,\beta) $ is fixed by $ \bar{P} $~\cite{Mayne2000}.
Indeed, we have
\begin{align*}
	&V_{\subm{\bar{P}}} ({\boldsymbol x}[k+1]) - V_{\subm{\bar{P}}} ({\boldsymbol x}[k]) \le - \sum_{i=1}^N W_{1,i} ({\boldsymbol x}_i[k], \bar{P}) \\
	&W_{1,i} (x_i, \bar{P}) := \Bigl\| {\boldsymbol B}_i^\top ({\boldsymbol A}_i^\top)^{\tau_{\rm h} -1} {\boldsymbol G}_{i,\tau_{\rm h}}^{-1} {\boldsymbol A}_i^{\tau_{\rm h}} \\
	&\hspace{4cm} \times \bigl(x_i - N\sum_{j=1}^N \bar{P}_{ij} x_j^{\sf d} \bigr) \Bigr\|^2 . 
\end{align*}
Given an equilibrium $ (x^{\sf e}, \beta^{\sf e}) $, let us take the optimal coupling $ P^{\sf e} := {\boldsymbol P}^* (x^{{\sf e}}) $ as $ \bar{P} $, and for $ \gamma > 0 $, define
\begin{equation}\label{eq:V}
V(x, \beta) :=  V_{\subm{P^{\sf e}}} (x) +  \gamma {d_\calH} ( \beta, \beta^{\sf e} ), \ (x,\beta)\in \bbR^{nN} \times  (\bbR_{>0}^N / {\sim} ) . 
\end{equation}
The following theorem follows from the fact that, for sufficiently small or large $ \varepsilon > 0$ and large $ \gamma > 0 $, $ V $ behaves as a Lyapunov function of \eqref{eq:x_global},~\eqref{eq:beta_global} with respect to $ (x^{\sf e}, \beta^{\sf e}) $.
We give the proof in Appendix~\ref{app:convergence}.

\begin{theorem}\label{thm:stability}
	Assume that for all $ i\in \bbra{N} $, $ {\boldsymbol A}_i $ is invertible.
	Then the following hold:
	\begin{itemize}
		\item[(i)] Assume that $ (x^{\sf e}, \beta^{\sf e}) $ is an isolated equilibrium of \eqref{eq:x_global}, \eqref{eq:beta_global}. Then, for a sufficiently large $ \varepsilon > 0 $, $ (x^{\sf e}, \beta^{\sf e}) $ is locally asymptotically stable.
		\item[(ii)]  
		Assume that $ x_i^{\sf d} \neq x_j^{\sf d} $ for all $(i,j), \ i\neq j $. \fin{Assume further that for some $ \varepsilon' > 0 $, $ (x^{\sf e} (\varepsilon'), \beta^{\sf e}(\varepsilon')) $ is an isolated equilibrium of \eqref{eq:x_global}, \eqref{eq:beta_global} and for some permutation $ \sigma $, $ (x^{\sf e} (\cdot), \beta^{\sf e}(\cdot)) \in {\rm Exp}(\sigma) $.}
		Then, for sufficiently small $ \varepsilon >0$, $ (x^{\sf e} (\varepsilon), \beta^{\sf e}(\varepsilon)) $ is locally asymptotically stable. \hfill$ \diamondsuit $
	\end{itemize}
\end{theorem}




\subsection{Interpretation of Sinkhorn MPC as an alternating descent and ascent method}
Lastly, we give an interpretation of Sinkhorn MPC.
We continue to work on the quadratic cost~\eqref{eq:energy} with the barycentric target~\eqref{eq:barycentric_target}.
	First, similarly to \eqref{eq:transport_cost_conti} for the continuous-time systems, consider the entropic OT cost
	\begin{equation}\label{eq:OT_cost_disc}
		E(x,x^{\sf d}) := \min_{P\in \calT(\frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N}, \frac{\scalebox{0.6}{\bf 1}_{\scalebox{0.4}{$N$}}}{N})} \ \sum_{i,j\in\bbra{N}} {\boldsymbol C}_{ij} (x) P_{ij} - \varepsilon {\rm H}(P),
	\end{equation}
	where $ {\boldsymbol C}_{ij}(x) = \|x_i - x_j^{\sf d} \|_{\calG_i}^2 $.
	Since it holds
	\begin{align*}
		\nabla_{x_i} E(x, x^{\sf d}) = 2\calG_i \sum_{j=1}^{N} {\boldsymbol P}_{ij}^* (x) (x_i - x_j^{\sf d}) ,
	\end{align*}
	\blu{a point $ x $ satisfying the following is a stationary point of $ E $.}
	\begin{equation}\label{eq:stationary_point}
		x_i - N \sum_{j=1}^{N} {\boldsymbol P}_{ij}^* (x) x_j^{\sf d} = 0, \ \forall i\in \bbra{N} .
	\end{equation}
	This is exactly the condition \eqref{eq:equilibrium_iff} for the equilibrium points $ x^{{\sf e}} $ of \eqref{eq:x_global}, \eqref{eq:beta_global}.
	\subm{Hence, Sinkhorn MPC can be viewed as a cost-effective search method to find the stationary points of the associated entropic OT cost $ E $.}


	Next, we introduce the dual problem associated \subm{with} \eqref{eq:OT_cost_disc}:
	\begin{align*}
		& E(x,x^{\sf d}) = \max_{{\bf f},{\bf g}\in \bbR^N} \ Q({\bf f},{\bf g}; x), \\
		&Q({\bf f},{\bf g}; x) := {\bf f}^\top ({\bf 1}_N / N) + {\bf g}^\top ({\bf 1}_N / N) - \varepsilon (\Ee^{{\bf f}/\varepsilon})^\top {\boldsymbol K}(x) \Ee^{{\bf g}/\varepsilon} ,
	\end{align*}
	where $ \Ee^{{\bf f}/\varepsilon} \in \bbR^{N} $ denotes the element-wise exponential of $ {\bf f}/\varepsilon $.
	Let $ {\bf f}[k] = \varepsilon \log (\alpha[k]), {\bf g}[k] = \varepsilon \log (\beta[k]) $ for scaling variables of the Sinkhorn algorithm~\eqref{eq:Sinkhorn_algorithm} where $ (\log(\alpha))_i = \log (\alpha_i) $. Then, the Sinkhorn iterations \eqref{eq:Sinkhorn_algorithm} are equivalent to a block coordinate ascent~\cite[Remark~4.21]{Peyre2019}, which updates alternatively $ {\bf f} $ and $ {\bf g} $ to cancel the respective gradients
	\begin{align*}
		&\nabla_{\bf f} Q({\bf f},{\bf g}; x) = {\bf 1}_N/N - \Ee^{{\bf f}/\varepsilon} \odot ({\boldsymbol K}(x)\Ee^{{\bf g}/\varepsilon}), \\
		&\nabla_{\bf g} Q({\bf f},{\bf g}; x) = {\bf 1}_N/N - \Ee^{{\bf g}/\varepsilon} \odot ({\boldsymbol K}(x)^\top \Ee^{{\bf f}/\varepsilon}) ,
	\end{align*}
	where $ \odot $ denotes element-wise multiplication.
	On the other hand, the gradient with respect to $ x_i $ is
	\begin{align}
		&\nabla_{x_i} Q({\bf f}[k+1], {\bf g}[k] ; x) \nonumber\\
		&= \sum_{j\in \bbra{N}} \alpha_i[k+1] \beta_j[k] \exp\left( - \frac{{\boldsymbol c}_{\tau_{\rm h}}^i (x_i,x_j^{\sf d})}{\varepsilon}  \right) \nabla_{x_i} {\boldsymbol c}_{\tau_{\rm h}}^i (x_i,x_j^{\sf d}) \nonumber\\
		&= \sum_{j\in \bbra{N}} P_{ij} [k] \nabla_{x_i} {\boldsymbol c}_{\tau_{\rm h}}^i (x_i, x_j^{\sf d}),
	\end{align}
	which has the same form as \eqref{eq:grad_E}. Now, let us consider the case when the dynamics~\eqref{eq:linear_dynamics_conti} is well approximated by the discretized system~\eqref{eq:linear_dynamics_disc}. 
	Then using the same derivation as for \eqref{eq:lyap_decrease} in the proof of Theorem~\ref{thm:lasalle}, along the trajectory $ x(t) $ following
	\begin{equation}
		\dot{x}_i (t) = A_i x_i (t) + B_i u_i^{\rm MPC} \left(x_i(t), x^{\rm tmp}_{i} (P[k]) \right) ,
	\end{equation}
	we have
	\begin{align}
		&\frac{{\rm d}}{{\rm d} t} \calQ\left({\bf f}[k+1], {\bf g}[k] ; x(t)\right) \nonumber\\
		&\qquad\begin{cases}
			< 0, & x_i(t)  \neq x_i^{\rm tmp} (\subm{P[k]}), \ \exists i\in \bbra{N}, \\
			= 0, & x_i(t)  = x_i^{\rm tmp} (\subm{P[k]}), \ \forall i\in \bbra{N} ,
		\end{cases}
	\end{align}
	where $ \calQ $ is the continuous-time version of $ Q $:
	\begin{align*}
		\calQ({\bf f},{\bf g}; x) &:= {\bf f}^\top ({\bf 1}_N / N) + {\bf g}^\top ({\bf 1}_N / N) - \varepsilon (\Ee^{{\bf f}/\varepsilon})^\top K(x) \Ee^{{\bf g}/\varepsilon} \\
		&\simeq Q({\bf f},{\bf g}; x).
	\end{align*}
	Therefore, if the sampling time is small, the state trajectory $ {\boldsymbol x}[k] $ moves in a direction where $ Q({\bf f}[k+1], {\bf g}[k]; {\boldsymbol x}[k]) $ decreases.
	Note that the above argument applies to the general cost under the assumptions in Theorem~\ref{thm:lasalle}.


	In summary, Sinkhorn MPC can be interpreted as an alternating descent and ascent method to seek a solution of the minimax problem
	\[
		\min_{x\in \bbR^{nN}} \max_{{\bf f},{\bf g}\in \bbR^N} \ Q({\bf f},{\bf g};x) ,
	\]
	where the minimizers satisfy \eqref{eq:stationary_point}.





\section{Conclusion}\label{sec:conclusion}
\magenta{
In this paper, we presented the concept of Sinkhorn MPC, which integrates MPC and the Sinkhorn algorithm to achieve scalable, cost-effective transport over dynamical systems.
The numerical examples described the usefulness of the proposed method.
Moreover, thanks to the entropy regularization, under some assumptions, we ensured the global convergence for Sinkhorn MPC, which is one of the most important properties of transport algorithms.
Furthermore, for linear systems with a quadratic cost, we analyzed the ultimate boundedness and the local asymptotic stability for Sinkhorn MPC based on the stability of the constrained MPC and the conventional Sinkhorn algorithm.

On the other hand, in the numerical example, we observed that the regularization parameter plays a key role in the trade-off between the stationary and transient behaviors for Sinkhorn MPC.
Hence, an important direction for future work is to investigate the design of a time-varying regularization parameter to balance the trade-off.
In addition, although we focused on the case where an OT problem is equivalent to an assignment problem, in Remark~\ref{rmk:nonuniform}, we mentioned that the convergence result for Sinkhorn MPC still holds in more general settings. Hence, it is also interesting to explore applications of Sinkhorn MPC for general OT problems.
\res{In this paper, for simplicity, we assumed the invertibility of $ B_i $ for the convergence analysis. Possible relaxation of this assumption will be reported in a future publication.}

}




\begin{ack}
	This work was supported in part by JSPS KAKENHI Grant Numbers JP21J14577, JP21H04875, JST, ACT-X Grant Number JPMJAX2102, and the joint project of Kyoto University and Toyota Motor Corporation, titled ``Advanced Mathematical Science for Mobility Society.''
\end{ack}


