\section{Introduction}

Electrocardiogram (ECG) remains the golden standard to study the physiology of the heart, to date. Nevertheless, the limitations of ECG, e.g., placement of electrodes on patient's chest and limbs by a trained professional, requirement of no chest movement during the experiment, etc., have inspired the researchers to explore and develop alternate methods as well. To this end, photoplethysmography (PPG) is a low-cost, non-invasive, optical technique that has recently gained widespread acceptance for (do-it-yourself) monitoring of cardiac activity. The off-the-shelf PPG devices (pulse oximeters, wearable smart bands, smart watches) are typically attached to the peripheral sites (e.g., wrist, finger, ear, etc.) to monitor various important physiological parameters, e.g., heart rate, heart rate variability, blood Oxygen saturation level (SpO2), and respiratory rate~\cite{A2}. More recently, the PPG signal has been used to aid in the diagnosis of several cardiac vascular, respiratory, and neurological diseases~\cite{A3}.

A PPG device utilizes light at two different wavelengths (typically, red and infrared) to measure the volumetric changes in blood flow inside blood vessels at the distal sites. That is, a PPG device is composed of a light emitting diode (LED) that emits (red and infrared) light, and a photodetector (PD) that collects the fraction of light reflected from the skin tissue. Two kinds of PPG devices are prevalent: transmissive and reflective. In transmissive type (oximeters), LED and PD are on the opposite side of the skin tissue, the LED impinges the light onto the skin tissue while the PD collects the light that penetrates through the tissue. In reflective type (smart watches), the PD is placed besides the LED on the same side of the skin tissue (see Fig. \ref{fig:1} a). More recently, video-PPG has started to attract the attention of the research community whereby one places his/her fingertip on the smartphone camera in order to record a small video which is then processed to extract the PPG signal (see Fig. \ref{fig:1} b).

\begin{figure}
    \centering
     \subfigure[]
    {\includegraphics[width=0.24\textwidth]{fig1a.png}}
   \subfigure[]
   {\includegraphics[width=0.24\textwidth]{fig1b.png}}

    \caption{(a) Oximeter-based PPG, (b) Smartphone-based video-PPG. }
    \label{fig:1}
\end{figure}

Having motivated the PPG techniques, there are certain challenges associated with the PPG waveform too, e.g., person-specific penetration depth of the light incident onto the skin tissue~\cite{A5}, variation in accuracy from one body-site to another, interference from nearby power-lines, baseline wander noise ~\cite{A6}, motion artefacts due to ambient light and due to voluntary or involuntary movement of the subject (e.g., during exercise, sleep).
Therefore, in order to make sure that the accuracy of the body vitals returned by a PPG device is not compromised, there is a pressing need to design efficient denoising algorithms to remove unwanted noise and motion-induced artefacts (e.g., baseline drift and spikes) from a noisy PPG signal. This is precisely the agenda of this work.

{\bf Contributions.}
The proposed method consists of two distinct phases. During the decomposition phase, the fast wavelet transform (FWT) is utilized to extract multiple approximation and detailed coefficients from a noisy PPG signal. During the reconstruction phase, a custom feedforward neural network (FFNN) systematically selects the (noise-free) wavelet sub-signals outputted by the inverse-FWT in order to help synthesize a clean PPG signal.
To train the FFNN in a supervised manner, gradients of mean squared error (MSE) loss function are computed for the backpropagation purposes.

{\bf Outline.}
The rest of this paper is organized as follows. Section II summarizes the selected related work. Section \ref{NN} provides key details of our proposed DL-assisted hybrid denoising method. Section IV provides selected numerical results. Finally, Section \ref{conc} concludes the paper.

\section{Related Work}

Due to space-constraint, this section summarizes only selected related works, as presented below.\footnote{For instance, there are works that design adaptive filters, perform singular value decomposition, or use the multi-scale principle component analysis for PPG denoising. These are not referenced due to the lack of space. }

The Fourier series analysis and bandpass filtering \cite{A7}, iterative singular spectrum-based method \cite{A9}, least mean square algorithm \cite{A12}, empirical mode decomposition combined with spectrum subtraction \cite{A13}, spectral matrix decomposition \cite{A14}, independent component analysis combined with wavelet \cite{A22}, discrete wavelet transform (DWT) \cite{A15}, multi-level DWT \cite{A16}, DWT along with adaptive filtering \cite{A17}, and periodic averaging filtering \cite{A8} are some of the techniques used for noise and artefact removal from PPG signals.

Recently, machine and deep learning based works on PPG denoising have emerged. For instance, \cite{A20} utilizes the principal component analysis along with an artificial neural network, while \cite{A1} considers a cycle generative adversarial network for noise removal from PPG signals.

\section{System Model and Problem Formulation}
\label{SM}

\begin{figure}
\centering
\includegraphics[width=3.4in]{Figures/SM.PNG}
\caption{System Model}
\centering
\label{fig:SM}
\end{figure}

We have designed multi-level fast wavelet decomposition with neural network to remove the noise from PPG signal. An input noisy  signal of length N with the window size of 8sec is given to wave decomposition block. In this the signal is decomposed into detail $cD_n$ and approximate coefficients $cA_n$. The approximate coefficient at each level helps obtain the next level of coefficients.

The same noisy input is imported into neural network which generates the neurons with the same output size as of decomposition. Dot product of the coefficients obtained from wave decomposition and values of output neuron is computed to obtain the reconstructed signal as shown in  Figure \ref{fig:SM}. Network is trained to update the neurons which helps cancels the effect of those coefficient which introduces artifact in the signal. Sharp spikes and noise will be eliminated from output signal.

\subsection{Case 1: Single Antenna without Interference}

\section{The Proposed Method}\label{NN}

This section first introduces the PPG denoising problem as an optimization program. Then, the key details of the proposed deep learning (DL) and wavelets-based hybrid method are outlined, as well as the architecture of the custom DL model.

\subsection{Problem Formulation}

The raw PPG signal acquired by a PPG device is typically contaminated by: the respiratory-induced baseline, hardware-induced high-frequency noise, artefacts due to motion and ambient light. Thus, in order to denoise the raw noisy PPG signal, we begin by taking the FWT of the noisy PPG signal in order to decompose the signal into ($L$) detail coefficients and one approximate coefficient (see Fig. \ref{fig:SM}).

\begin{figure}
\centering
\includegraphics[width=3.5in]{Figures/SM.PNG}
\caption{The proposed deep learning \& the wavelets-based hybrid method. }
\centering
\label{fig:SM}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=2.7in]{Figures/Recnstrctd_sub_sig.png}
\caption{Reconstruction of sub-signals by the wavelet reconstruction block.}
\centering
\label{fig:sub_sig}
\end{figure}

The output of the wavelet reconstruction block consists of the sub-signals $\mathbf{s}_i$ that together satisfy the following equation:
\begin{equation}
\label{eq:2}
    \textbf{y} =  \textbf{s}_{1}  + \textbf{s}_{2} +... + \textbf{s}_{L+1},
\end{equation}
where  $\textbf{s}_i \in \mathcal{R}^{N \times 1}$  is the sub-signal reconstructed from the $i$-th detailed coefficient (while setting all other coefficients to 0). $\textbf{s}_{L+1}$ is the sub-signal obtained using the approximate coefficient only (with all detailed coefficients being 0). For a quick illustration, Fig. \ref{fig:sub_sig} shows a few such sub-signals ($\mathbf{s}_1$ to $\mathbf{s}_7$) outputted by the wavelet reconstruction block.

Equation (\ref{eq:2}) suggests that the noisy PPG signal $\textbf{y}$ can be denoised by selectively suppressing some sub-signals. This implies that the denoised PPG signal can be written as follows:
\begin{equation}
\label{eq:3}
    \hat{\textbf{x}} = a_{1} \textbf{s}_{1}  +a_{2} \textbf{s}_{2} +... +a_{L+1} \textbf{s}_{L+1},
\end{equation}
where $a_i$ equals $0$ for the sub-signals $\textbf{s}_i$ that correspond to the noise and artefacts, and $1$ for the sub-signals which (when summed together) yield the clean PPG signal. Thus, the denoising problem boils down to the optimal selection of $a_i$'s. Let us denote $\textbf{a} = [a_1, a_2, ..., a_{L+1}]^T$
and $\mathbf{S} = [\textbf{s}_1, \textbf{s}_2, ... , \textbf{s}_{L+1}]$,
then we can rewrite (\ref{eq:3}) as follows:
\begin{equation}
\label{eq:4}
    \hat{\textbf{x}} = \mathbf{S}\textbf{a}.
\end{equation}
To obtain the denoised signal $\hat{\textbf{x}}$, the optimal $\textbf{a}$ can be obtained by solving the following optimization problem:

\begin{equation}
\begin{aligned}
\min_{a_{1},a_{2} ... a_{L+1} } \quad & || \mathbf{S} \mathbf{a} -\mathbf{x} ||_2^2 \\
\textrm{s.t.} \quad & {a_{i}} \in \{0,1\}, i = 1,2,...,L+1 \\
\end{aligned}
\end{equation}
where $\mathbf{x}$ is the clean PPG signal. The binary integer constraint makes this optimization problem an integer linear program that is difficult to solve. Thus, we relax the integer constraint by allowing $a_i$ to take any value between $0$ and $1$. It helps since the constraint can now be satisfied by our DL model by using the sigmoid activation function at the output layer. Thus, we reach the following optimization problem (which also serves as the loss function for our DL model):

\begin{equation}
\begin{aligned}
\min_{a_{1},a_{2} ... a_{L+1} } \quad & || \mathbf{S} \mathbf{a} -\mathbf{x} ||_2^2 \\
\textrm{s.t.} \quad & 0 \leq {a_{i}} \leq 1, i = 1,2,...,L+1 \\
\end{aligned}
\end{equation}

\subsection{Proposed DL and Wavelet-based Hybrid Approach}
Under the proposed method, the noisy PPG signal goes as input to the wavelet decomposition block, while the wavelet reconstruction block returns the sub-signals $\mathbf{s}_i$ that satisfy (\ref{eq:2}) (see Fig. \ref{fig:SM}). The noisy PPG signal is also fed to the DL model in parallel, which generates the optimal weight vector $\textbf{a}$. The product of the matrix $\mathbf{S}$ (generated by the wavelet reconstruction block) with the vector $\textbf{a}$ (generated by the DL model) is the denoised signal. The MSE between the denoised signal and the clean PPG signal is utilized as a loss function to compute the gradients, which are then backpropagated in order to train the FFNN in a supervised manner.

\textit{The FFNN Architecture:}
The DL model used in the proposed method is an FFNN (see Fig. \ref{fig:NN}). It consists of an input layer of size $N$, three hidden layers, and an output layer. The number of neurons successively halve from one hidden layer to the next, with $L+1$ neurons in the output layer. Each hidden layer consists of a dense layer followed by batch normalization and activation (which in our case is ReLu). Finally, the output layer utilizes the sigmoid activation function.

\begin{figure}
\centering
\includegraphics[width=2.9in]{Figures/NN.PNG}
\caption{Architecture of the custom FFNN.}
\centering
\label{fig:NN}
\end{figure}

\subsection{Performance Metrics}
The performance of the proposed method is evaluated using the following two performance metrics:
MSE: $= \frac{1}{M} \sum_{j=1}^{\ M} (\textbf{x}_j-\hat{\textbf{x}}_j)^2$, and peak signal-to-noise ratio (PSNR): $= 10  {\log{_{10}}} (1/\textrm{MSE})$ in decibels, where $\textbf{x}_j$ is the clean PPG signal, $\hat{\textbf{x}}_j$ is denoised PPG signal, and $M$ is the batch size.

\section{Numerical Simulation}
\section{Performance Evaluation}
\label{performance}

We first provide key details about the training of the custom FFNN presented in Fig. 4, followed by a brief description of the two simulation experiments, and finally the results.

\textit{Training of the custom FFNN.}
All the trainable parameters are initialized by Xavier initializer. The model is optimized using Adam optimizer with an initial learning rate (LR) of 0.001. To avoid the over-fitting, an early stopping is applied. The details of other hyper-parameters are given in Table \ref{tab:table3}.

\begin{table}
\caption{\label{tab:table3}Hyper-parameters for the custom FFNN.}
\begin{center}
\begin{tabular}{c*{2}{c}}
\\\hline
Hyper-parameter & Value \\\hline\hline
Training dataset size & $3873$\\\hline
Test dataset size & $240$ \\\hline
Training batch size & $100$ \\\hline
Validation batch size & $100$ \\\hline
Number of epochs & $500$ \\\hline
\end{tabular}
\end{center}
 \end{table}

We now describe each of the two experiments performed to validate the performance of the proposed method.

\textit{Experiment I.}
The first experiment is a controlled one. We fetch clean PPG signals from the BIDMC dataset \cite{BIDMC}, and corrupt them by adding four types of noises which are known to degrade the PPG signals quite often. Specifically, we consider the following noise types: Gaussian (zero-mean, variance$=0.05$), Poisson (with $\lambda =0.02$), uniform ($[ -0.1,0.1]$), and salt-and-pepper noise (with a density$=0.05$). Fig. \ref{fig:Noisy PPG} shows a snapshot of the resulting noisy PPG signals.

\begin{figure}
\centering
\includegraphics[width=3.2in]{Figures/Picture_1.PNG}
\caption{Some noisy PPG signals.}
\centering
\label{fig:Noisy PPG}
\end{figure}

Next, keeping in mind the multi-level decomposition capability of the FWT, we investigate the impact of the depth of decomposition on the denoising performance of the proposed method. Fig. 6 shows that the denoising performance of the proposed method improves with higher wavelet decomposition levels. This effect is most prominent on the PPG signal corrupted by the salt-and-pepper noise.
This suggests using the maximum level of decomposition (which is 8 in this work). Thus, in the subsequent simulations, we consider an 8-level decomposition. For Fig. 6, mother wavelet used was 'db4'.

Further, the unique morphology of PPG signals prompted us to study the denoising performance of the proposed method for various mother wavelets. Fig. \ref{fig:mse} shows that the MSE of the denoised signal is reduced significantly, compared to the noisy signal for all mother wavelets. However, 'db10' wavelet outperforms all other wavelets for all noise types.

 \begin{figure}
\centering
\includegraphics[width=2.8in]{Figures/mse.PNG}
\caption{Impact of the depth of decomposition on MSE. }
\centering
\label{fig:msedb4}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=2.8in]{Figures/MSE_impact_wavelets.PNG}
\caption{MSE of the proposed method under 10 different mother wavelets.}
\centering
\label{fig:mse}
\end{figure}

\textit{Experiment II.}
The second experiment is more ambitious. Here, we utilize clean PPG signals from the BIDMC dataset as before. Additionally, this time we add the artefact signals extracted from our custom video-PPG dataset (by doing FWT decomposition and retaining those detailed coefficients that capture the artefacts) to the training data. This mixed training dataset is then used to train the FFNN of Fig. 4. Note, that for the MSE computation, the clean PPG signal (for the video-PPG data) is obtained by doing the FWT decomposition of the PPG signal, and by discarding those detailed coefficients that captured the artefacts. Once the FFNN is trained, we evaluate its performance by denoise the video-PPG data during the testing phase. Fig. \ref{fig:NS} shows that our proposed method efficiently removes the baseline and high-frequency noise from the noisy PPG signal extracted from the video-PPG data.

\begin{figure}
\centering
\includegraphics[width=2.8in]{fig8.png}
\caption{Denoised video-PPG signal. }
\centering
\label{fig:NS}
\end{figure}

\section{Conclusion}
\label{conc}
Under the proposed method, the noisy PPG signal was simultaneously fed to a multi-level FWT decomposition block and a custom FFNN. The FFNN guided the FWT reconstruction phase by systematically setting the wavelet sub-signals corresponding to artefacts and noise to zero.
The FFNN was trained in a supervised manner (using the MSE loss function). Simulation results revealed that our proposed method reduces the MSE of the PPG signal significantly (compared to the MSE of the original noisy PPG signal): by 56.40\% for Gaussian noise, by 64.01\% for Poisson noise, 46.02\% for uniform noise, and by 72.36\% for salt-and-pepper noise (with 'db10' mother wavelet).


\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}