\section{Introduction}
Model Predictve Control (MPC) is a well-established and widely used method to control complex dynamical systems. MPC requires the real-time solution of a sequence of optimal control problems (OCPs) on a finite (but large) time horizon, see e.g.\ \cite{grune2017, rawlings2019}, which can be computationally demanding. This is for example the case when the model is the result of the (spatial) discretization of a Partial Differential Equation (PDE) or in the simulation of interaction particle systems. It is therefore natural to use numerically efficient approximations to speed up the solution of the finite-horizon OCPs that need to be solved in MPC. 
One recently-proposed numerically-efficient approximation method is the Random Batch Method (RBM) \cite{jin2020,rbm},  which is closely related to the stochastic algorithms like Stochastic Gradient Descent (SGD) used in machine learning. 
The effectiveness of the resulting combination of MPC with RBM (RBM-MPC) has been demonstrated for interacting particle systems in \cite{ko2021model}, but a rigorous stability and convergence analysis of the proposed method was still lacking. 


The idea for MPC can be traced back to Propoi \cite{propoi1963} and is already found in the classical book by Lee and Markus, \cite{lee1967}, but the interest in MPC only started to grow since the late 1970's when sufficient computational power for the implementation of MPC became available. Since then, various analysis methods for MPC have been proposed. Often terminal constraints or terminal costs are imposed to guarantee stability, see, e.g., \cite{mayne2000}, \cite{grune2017}, and  \cite{rawlings2019} for an overview. Because this paper considers RBM-MPC for unconstrained linear systems, the stability and convergence of MPC can be understood simply from the convergence of the solution of the Algebraic Riccati Equation (ARE) to the solution of the Riccati Differential Equation (RDE) here, see also \cite{mpc}. 

Random Batch Methods (RBMs) have been introduced recently in \cite{jin2020} and are inspired by the successes of stochastic algorithms such as SGD  in supervised learning. In the first formulation of RBM, the simulation of interacting particle systems is sped up by considering not all interactions between particles simultaneously but only a randomly chosen subset at each time instant \cite{jin2020}. Although initially proposed for the simulation of the forward dynamics, the same ideas have also been applied to speed up the computation of optimal controls for interacting particle systems \cite{ko2021model}. A first rigorous analysis of RBMs in linear-quadratic optimal control has been provided in \cite{rbm}. The same paper also shows that original idea of the RBM for interacting particle systems can be extended to a more general operator-splitting setting. 

In this paper, we provide the first rigorous analysis of the RBM-MPC algorithm in the unconstrained linear quadratic setting. The obtained error estimates clearly demonstrate the influence of the different tuning parameters on the expected performance of the proposed algorithm and the obtained convergence rates are validated in numerical experiments.

The remainder of this paper is structured as follows. 
The proposed RBM-MPC algorithm is presented in Section \ref{sec:rbm_mpc} and the main stability and convergence results are presented in Section \ref{sec:main_results}. Section \ref{sec:preliminaries} summarizes the existing results about the analysis of MPC and RBM from \cite{mpc} and \cite{rbm} that form the basis for the proofs of the main result in Section \ref{sec:proofs}. The obtained convergence rates and the effectiveness of the proposed method are demonstrated in numerical experiments in Section \ref{sec:numerical}. Finally, conclusions and perspectives for future work are presented in Section \ref{sec:conclusions}. 

\section{The RBM-MPC algorithm} \label{sec:rbm_mpc}

\subsection{The infinite-horizon problem}
The RBM-MPC algorithm analyzed in this paper is a way to approximate the control $\bm{u}^*_\infty(t)$ that minimizes
\begin{equation}
    J_\infty(\bm{u}) = \int_0^\infty \left((\bm{x}(t))^\top Q \bm{x}(t) + (\bm{u}(t))^\top R \bm{u}(t) \right) \ \mathrm{d}t, \label{eq:Jinf}
\end{equation}
subject to the dynamics
\begin{equation}
    \dot{\bm{x}}(t) = A \bm{x}(t) + B \bm{u}(t), \qquad \qquad \bm{x}(0) = \bm{x}_0,  \label{eq:dyn_x}
\end{equation}
where $\bm{x}(t)$ evolves in $\mathbb{R}^n$ starting from the initial condition $\bm{x}_0$, the control $\bm{u}(t)$ evolves in $\mathbb{R}^m$, $Q \in \mathbb{R}^{n \times n}$ is symmetric positive semi-definite, $R \in \mathbb{R}^{m \times m}$ is symmetric positive definite, $A \in \mathbb{R}^{n\times n}$, and $B \in \mathbb{R}^{n\times m}$. It is assumed that $(A,B)$ is stabilizable and $(A,Q)$ is detectable.  

\begin{remark}
It is well-known that the problem \eqref{eq:Jinf}--\eqref{eq:dyn_x} admits a unique minimizer of feedback form, see e.g.\ \cite{sage1968}, which can be computed efficiently by solving an Algebraic Riccati Equation (ARE) when $n$ is not very large. MPC and RBM-MPC thus only offer a computational advantage when the ARE cannot be solved efficiently (i.e., when $n$ is large) or when the optimal control is not given by a linear feedback law (i.e., in constrained and/or nonlinear optimal control). 
\end{remark}

\begin{remark}
Although we present the RBM-MPC algorithm here in the LQ setting, the same approach can be used for problems with constraints and/or nonlinear dynamics, see \cite{ko2021model}. The stability and convergence analysis for the LQ problem in this paper is the first step towards a rigorous analysis of the RBM-MPC algorithm in these more general settings. 
\end{remark}

\subsection{Model Predictive Control} \label{ssec:mpc}
One method to approximate the solution of the infinite-horizon OCP \eqref{eq:Jinf}--\eqref{eq:dyn_x} is Model Predictive Control (MPC). 
In MPC, two time scales arise: the prediction horizon $T$ and the shorter control horizon $\tau$. Set $\tau_i := i \tau$ (with $i \in \mathbb{N}$) and let $\bm{u}^*_{T}(t;\bm{x}_i,\tau_i)$ and $\bm{x}^*_{T}(t;\bm{x}_i,\tau_i)$ denote the control and state trajectory that minimize
\begin{align}
\begin{split}
     J_{T}&(\bm{u}_{T}; \bm{x}_i, \tau_i) = \\ &  \int_{\tau_{i}}^{\tau_{i} + T} \bm{x}_{T}(t)^TQ\bm{x}_{T}(t) + \bm{u}_{T}(t)^TR\bm{u}_{T}(t)\,dt, \label{eq:JT}
\end{split}
\end{align}
where $\bm{x}_{T}(t)$ (for $t \in [\tau_i, \tau_i + T]$) fulfills 
\begin{align}
    \dot{\bm{x}}_{T}(t) = A\bm{x}_{T}(t) + B\bm{u}_{T}(t),\;\;\;\;\;\;\;\;\bm{x}_{T}(\tau_i) = \bm{x}_i. \label{eq:dyn_xT}
\end{align}
The MPC control $\bm{u}_M(t)$ and corresponding state trajectory $\bm{x}_{M}(t)$ are now computed as follows, see, e.g.\ \cite{mpc, azmi2016, azmi2018}.
\begin{enumerate}
    \item Initialize $\bm{x}_{M}(0) = \bm{x}_0$ and $i = 0$.
    \item Compute $\bm{u}_T^{*}(t;\bm{x}_{M}(\tau_i), \tau_i)$ and $\bm{x}_T^{*}(t;\bm{x}_{M}(\tau_i), \tau_i)$ on $[\tau_i, \tau_i + T]$.
    \item Set $\bm{u}_M(t) = \bm{u}^*_T(t;\bm{x}_{M}(\tau_i), \tau_i)$ and $\bm{x}_M(t) = \bm{x}^*_{T}(t,\bm{x}_{M}(\tau_i), \tau_i)$ on $[\tau_i, \tau_{i+1}]$.
    \item Set $i = i + 1$ and got to Step 2. 
\end{enumerate}
The MPC algorithm is illustrated in Figure \ref{fig:mpc}. 


\begin{figure}
\includegraphics[width=0.98\columnwidth]{Fig1}
\caption{Illustration of the MPC algorithm. The controls $u_T^*(t;\bm{x}_i, \tau_i)$ are computed on $[\tau_i, \tau_i+T]$ (dashed black lines) but are only applied on $[\tau_i, \tau_{i+1}]$ (solid red lines) for $i \in \{1,2 \}$, yielding $x_M(t)$.}
\label{fig:mpc}
\end{figure}

\subsection{The Random Batch Method} \label{ssec:rbm}
The Random Batch Method (RBM) can speed up the solution of the OCP \eqref{eq:JT}--\eqref{eq:dyn_xT}. 
We start by splitting the system matrix $A$ into submatrices $A_m$
\begin{align}
    A = \sum_{m=1}^M A_m. \label{eq:Asplit}
\end{align}
The time interval $[\tau_i, \tau_i + T]$ is divided into $K$ subintervals of length $h$. The grid points are denoted by
\begin{equation}
    h_{i,k} = \tau_i + k h = i \tau + k h, \qquad i \in \mathbb{N}, k \in \{0,1,\ldots, K \}.
\end{equation}

Note that there are $2^M$ subsets of the indices $m \in \{1,2, \ldots, M \}$, which are denoted by $S_1, S_2, \dots, S_{2^M}$. In each time interval $[h_{i,k-1}, h_{i,k})$, one subset $S_{\omega_{i,k}}$ (with $\omega_{i,k} \in \{1,2,\ldots,2^M\}$) is chosen according to certain chosen probabilities $p_\omega$ (with $\omega \in \{1,2,\ldots,2^M\}$). In particular, $p_\omega \in [0,1]$ is the probability that $S_\omega$ is selected in a certain time interval and $\sum_{\omega} p_\omega = 1$. The indices for one interval $[\tau_i, \tau_i + T]$ are collected in a vector 
\begin{align}
    \bm{\omega}_i = (\omega_{i,1}, \omega_{i,2},\dots, \omega_{i,K}) \in \{1,2,\dots,2^M\}^K. \label{eq:def_omegai}
\end{align}
This vector $\bm{\omega}_i$ defines a time-dependent matrix on $[\tau_i, \tau_i + T]$
\begin{align}
    A_R(\bm{\omega}_i,t; \tau_i) = \sum_{m \in S_{\omega_{i,k}}} \frac{A_m}{\pi_m}, \quad \quad t\in [h_{i,k-1},h_{i,k}), \label{eq:def_AR}
\end{align}
where $\pi_m$ denotes the probability of having the index $m$ in the selected subset
\begin{equation}
    \pi_m := \sum_{\substack{\omega \in \{1,2,\dots,2^M\} \\ m \in S_\omega}} p_\omega. 
\end{equation}
Observe that the definition of $A_R(\bm{\omega}_i, t; \tau_i)$ requires that $\pi_m > 0$ for all $m \in \{1,2, \ldots, M \}$. 
The distance between $A_R$ and $A$ is measured by
\begin{align}
    \mathrm{Var}[A_R] &:= \sum_{\omega = 1}^{2^M} \left\| A - \sum_{m \in S_{\omega}} \frac{A_m}{\pi_m}  \right\|^2 p_\omega \label{eq:def_varA}
\end{align}
Note that $\mathrm{Var}[A_R]$ depends on the chosen splitting \eqref{eq:Asplit} and the probability distribution $p_\omega$ and that $\mathbb{E}[\| A_R(t) - A \|^2] = \mathrm{Var}[A_R]$ for all $t \in [\tau_i, \tau_i + T]$. The analysis of the RBM also requires a uniform growthbound $\mu_R \geq 0$ for $A_R(\bm{\omega}_i,t; \tau_i)$, i.e.\
\begin{equation}
\bm{x}^\top A_R(\bm{\omega}_i, t; \tau_i) \bm{x} \leq \mu_R |\bm{x}|^2, \label{eq:def_muR}
\end{equation}
for all $\bm{x} \in \mathbb{R}^n$, $\bm{\omega}_i \in \{1,2, \ldots, 2^M \}^K$, and $t \in [\tau_i, \tau_i + T]$. When $\mu_R > 0$, the RBM-error grows exponentially in time.  When $A$ is dissipative, it is often possible to choose all submatrices $A_m$ in \eqref{eq:Asplit} dissipative and achieve $\mu_R = 0$. This is the case for the example in Section \ref{sec:numerical} and \cite[Section 4]{rbm}. 

The idea is now that replacing the matrix $A$ in \eqref{eq:dyn_xT} by $A_R(\bm{\omega}_i, t; \tau_i)$ will lead to a significant reduction in computational cost and that the dynamics generated by $A_R(\bm{\omega}_i, t; \tau_i)$ is (in expectation) close to the dynamics generated by the original matrix $A$ when $h \mathrm{Var}[A_R]$ is sufficiently small. 

It therefore makes sense to replace the original matrix $A$ in the OCP \eqref{eq:JT}--\eqref{eq:dyn_xT} by $A_R(\bm{\omega}_i, t; \tau_i)$ and to introduce $\bm{u}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ and $\bm{x}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ as the control and state trajectory that minimize
\begin{align}
\begin{split}
     & J_R(\bm{u}_R; \bm{\omega}_i, \bm{x}_i, \tau_i) = \\ & \int_{\tau_i}^{\tau_i+T} \bm{x}_{R}(\bm{\omega}_i,t)^TQ\bm{x}_{R}(\bm{\omega}_i,t) + \bm{u}_{R}(t)^TR\bm{u}_{R}(t)\,dt, \label{eq:JR}
\end{split}
\end{align}
where $\bm{x}_{R}(t)$ (for $t \in [\tau_i, \tau_i+T]$) fulfills 
\begin{multline}
    \dot{\bm{x}}_{R}(\bm{\omega}_i, t) = A_R(\bm{\omega}_i, t; \tau_i)\bm{x}_{R}(\bm{\omega}_i, t) + B\bm{u}_{R}(t),\\
    \bm{x}_{R}(\bm{\omega}_i,\tau_i) = \bm{x}_i. \label{eq:dyn_xR}
\end{multline}
Note that $\bm{x}_R(\bm{\omega}_i, t)$ depends on $\bm{\omega}_i$ through $A_R(\bm{\omega}_i, t; \tau_i)$. 

\begin{remark} \label{rem:RBMGD}
When $n$ is large, the OCP \eqref{eq:JR}--\eqref{eq:dyn_xR} is usually solved by a gradient-based algorithm in which the gradient is computed as
\begin{equation}
    \nabla J_R(\bm{\omega}_i,t) = 2B^\top\bm{\varphi}(\bm{\omega}_i,t) + 2 R \bm{u}_R(t),
\end{equation}
where the adjoint state $\bm{\varphi}(\bm{\omega}_i,t)$ satisfies
\begin{equation}
    \dot{\boldsymbol \varphi}(\bm{\omega}_i,t) = (A_R(\bm{\omega}_i,t;\tau_i))^\top \bm{\varphi}(\bm{\omega}_i,t) + Q \bm{x}_R(\bm{\omega}_i,t), \label{eq:dyn_phiR}
\end{equation}
and the final condition $\bm{\varphi}(\bm{\omega}_i,\tau_i + T) = 0$. Note that $\bm{\omega}_i$ is fixed during this iteration process and that the matrix $A_R(\bm{\omega}_i,t;\tau_i)$ in \eqref{eq:dyn_xR} is the same as the matrix $A_R(\bm{\omega}_i,t;\tau_i)$ in \eqref{eq:dyn_phiR}. An approach that does not require to store $\bm{\omega}_i$ in memory would be to use a different realization of the matrix $A_R(\bm{\omega}_i,t;\tau_i)$ in every forward pass \eqref{eq:dyn_xR} and backward pass \eqref{eq:dyn_phiR}. The analysis and the formulation of good termination conditions for this algorithm are interesting topics for future research. 
\end{remark}

\subsection{The RBM-MPC algorithm} \label{ssec:rbmmpc}
The controls $\bm{u}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ will be used to control the original dynamics \eqref{eq:dyn_x}. Introduce therefore $\bm{y}_R^*(\bm{\omega}_i,t; \bm{x}_i, \tau_i)$ as the state trajectory that results from applying the control $\bm{u}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ to the plant \eqref{eq:dyn_xT}, i.e.
\begin{multline}
    \dot{\bm{y}}_R^*(\bm{\omega}_i, t) = A \bm{y}_R^*(\bm{\omega}_i, t) + B \bm{u}^*_R(\bm{\omega}_i, t; \bm{x}_i, \tau_i), \\ \bm{y}_R^*(\bm{\omega}_i, \tau_i) = \bm{x}_i. \label{eq:dyn_yR}
\end{multline}
where $\bm{y}_R^*(\bm{\omega}_i, t)$ denotes $\bm{y}_R^*(\bm{\omega}_i,t; \bm{x}_i, \tau_i)$, again for brevity. 

The RBM-MPC algorithm computes the control $\bm{u}_{R-M}(t)$ and corresponding state trajectory $\bm{x}_{R-M}(t)$ as follows. 
\begin{enumerate}
    \item Initialize $\bm{x}_{R-M}(0) = \bm{x}_0$ and $i = 0$.
    \item Select indices $\omega_{i,k}$ for $k \in \{1,2, \ldots, K \}$ according to the probability distribution $p_\omega$ and build $\bm{\omega}_i$ as in \eqref{eq:def_omegai}. 
    \item Compute $\bm{u}_R^{*}(\bm{\omega}_i, t;\bm{x}_{R-M}(\tau_i), \tau_i)$ and \\ $\bm{y}_R^{*}(\bm{\omega}_i, t;\bm{x}_{R-M}(\tau_i), \tau_i)$ on $[\tau_i, \tau_i + T]$.
    \item Set $\bm{u}_{R-M}(t) = \bm{u}_R^{*}(\bm{\omega}_i, t;\bm{x}_{R-M}(\tau_i), \tau_i)$ and \\ 
    $\bm{x}_{R-M}(t) = \bm{y}^*_R(t,\bm{x}_{R-M}(\tau_i), \tau_i)$ on $[\tau_i, \tau_{i+1}]$.
    \item Set $i = i + 1$ and got to Step 2. 
\end{enumerate}
Note that $\bm{x}_{R-M}(\tau_{i+1})$ depends on all indices $\omega_{j,k}$ with $j \in \{1, \ldots, i\}$ and $k \in \{ 1,2, \ldots, K \}$ which are collected in 
\begin{equation}
    \Omega_{i} = (\bm{\omega}_0, \bm{\omega}_1, \bm{\omega}_2, \ldots, \bm{\omega}_{i}) \in \{1,2, \ldots, 2^M \}^{(i+1)K}. \label{eq:def_Omegai}
\end{equation}

\begin{remark}
The proposed algorithm can be viewed as an MPC strategy with an imperfect plant model constructed by the RBM. The feedback nature of MPC will naturally create some robustness against these imperfections. 
\end{remark}

\begin{remark}
The computational advantage of the RBM-MPC strategy depends on the splitting of $A$ into submatrices $A_m$ in \eqref{eq:Asplit}. For example, when the original matrix $A$ is not tridiagonal and the submatrices $A_m$ are (up to a permutation of the rows and columns), the computational cost for one time step in an implicit time discretization scheme reduces from $O(n^3)$ to $O(n)$, see also Section \ref{sec:numerical}. 
A further discussion on the reduction of the computational cost achievable in RBM-constrained OCPs can be found in \cite[Examples 4 and 5]{rbm}. 
\end{remark}

\begin{remark} \label{rem:MPCRBM} 
The proposed algorithm is the result of first applying MPC and then the RBM (first-MPC-then-RBM). Another variant of RBM-MPC is obtained by first applying the RBM and then MPC (first-RBM-then-MPC). This means that we first construct a randomized $A$-matrix $\tilde{A}_R : \{1,2,\ldots,2^M \}^\mathbb{N} \times [0,\infty) \rightarrow \mathbb{R}^{n \times n}$ as in Subsection \ref{ssec:rbm}. The plant \eqref{eq:dyn_x} is then controlled by an MPC controller that uses the imperfect plant model based on $\tilde{A}_R$. When $\tau = T$, first-MPC-then-RBM and first-RBM-then-MPC result in the same algorithm, but for $\tau < T$ this is not the case. This phenomenon is similar to the discussions on first-optimize-then-discretize (FOTD) versus first-discretize-then-optimize (FDTO) in the numerics of optimal control problems, see, e.g.,  \cite{hinze2008, ervedoza2013}.
\end{remark}

\section{Main Results}
\label{sec:main_results}
Main Result \ref{res:main_1} concerns the stability of the RBM-MPC trajectory $\bm{x}_{R-M}(t)$.
\begin{Main_Result}\label{res:main_1}
There exist positive constants $M_\infty$, $\mu_\infty$, $C_{\mathrm{MPC}}$, $C_{[A,B,Q,R,T]}$ such that for $t \geq 0$
\begin{align}\label{eq:mr_stab}
    \mathbb{E}[|\bm{x}_{R-M}(t)|] \leq M_\infty e^{-t\,\mu_{R-M}} |\bm{x}_0|,
\end{align}
where
\begin{multline}
    \mu_{R-M} = \mu_\infty - C_{\mathrm{MPC}} e^{-\mu_\infty(T-\tau)} \\ - C_{[A,B,Q,R,T]}e^{2\mu_RT} \sqrt{h \mathrm{Var}[A_R]} . \label{eq:muRM}
\end{multline}
\end{Main_Result}

\begin{remark}
The constants $M_\infty$ and $\mu_\infty$ depend only on $A$, $B$, $Q$, and $R$ and are related to the infinite-horizon OCP \eqref{eq:Jinf}--\eqref{eq:dyn_x}, see \eqref{eq:expon_stab_infty} in Subsection \ref{ssec:Riccati}. The constant $C_{\mathrm{MPC}}$ also depends only on $A$, $B$, $Q$, and $R$ and also appears in stability analysis of MPC in Lemma \ref{lem:MPC} in Subsection \ref{ssec:mpc_ana}. 
\end{remark}


\begin{remark}
For $T - \tau$ sufficiently large and $h\mathrm{Var}[A_R]$ and sufficiently small, $\mu_{R-M} > 0$. When $\mu_{R-M} > 0$, the RBM-MPC strategy is stabilizing with probability 1. To see this, note that Markov's inequality \cite{rohatgi2015} shows that for any $\varepsilon > 0$
\begin{equation}
\mathbb{P}[|\bm{x}_{R-M}(\Omega_i, t)| \geq \varepsilon] \leq \frac{\mathbb{E}[|\bm{x}_{R-M}(t)|]}{\varepsilon}. 
\end{equation}
Because $\mathbb{E}[|\bm{x}_{R-M}(t)|] \rightarrow 0$ for $t \rightarrow \infty$ if $\mu_{R-M} > 0$, the probability that $|\bm{x}_{R-M}(\Omega_i, t)| < \varepsilon$ approaches 1 for $t \rightarrow \infty$. 
\end{remark}

\begin{remark} \label{rem:asymptotics}
For clarity, the estimates only include the leading terms in the limit $T-\tau \rightarrow \infty$ and $h \mathrm{Var}[A_R] \rightarrow 0$. In particular, by rescaling time, we may assume without loss of generality that $T \geq 1$ and $\tau \leq 1$. Also, terms proportional to $h \mathrm{Var}[A_R]$ are omitted because they are negligible to terms proportional to $\sqrt{h \mathrm{Var}[A_R]}$ in the limit $h \mathrm{Var}[A_R] \rightarrow 0$. 
\end{remark}

The second main result shows that for $h \mathrm{Var}[A_R] \rightarrow 0$, the RBM-MPC state trajectory $\bm{x}_{R-M}(\Omega_i, t)$ and control $\bm{u}_{R-M}(\Omega_i, t)$ converge (in expectation) to the MPC state trajectory $\bm{x}_M(t)$ and control $\bm{u}_M(t)$, respectively. 

\begin{Main_Result} \label{res:main_2}
Suppose that $\mu_{R-M} > 0$, then there exists a constant $C_{[A,B,Q,R,T]}$ such that
\begin{multline}
	\mathbb{E}[|\bm{x}_{R-M}(t) - \bm{x}_M(t)|] + \mathbb{E}[|\bm{u}_{R-M}(t) - \bm{u}_M(t)|] \leq \\
	C_{[A,B,Q,R,T]}e^{2\mu_R T}\sqrt{h \mathrm{Var}[A_R]}\frac{e^{-\mu_{R-M}(t-\tau)}(t+1)}{\mu_{R-M}}|\bm{x}_0|.
\end{multline}
\end{Main_Result}

\begin{remark}
The estimates in this paper lead to a constant $C_{[A,B,Q,R,T]}$ proportional to $T^{5/2}$. 
The numerical example in Section \ref{sec:numerical} suggests that $C_{[A,B,Q,R,T]}$ increases at a lower rate. 
\end{remark}

\begin{remark}
The proof also clearly indicates that the constant $C_{[A,B,Q,R,T]}$ decreases when $\tau$ is smaller. However, it does not hold that $C_{[A,B,Q,R,T]} \rightarrow 0$ for $\tau \rightarrow 0$. The dependence of $C_{[A,B,Q,R,T]}$ on $\tau$ was therefore omitted. 
\end{remark}

Main Result \ref{res:main_2} has the following corollary that relates the RBM-MPC strategy to the infinite horizon problem \eqref{eq:Jinf}--\eqref{eq:dyn_x}. 

\begin{corollary} \label{corr:RBMMPCconv}
Suppose that $\mu_{R-M} > 0$, then there exist constants $C_{[A,B,Q,R]}$ and $C_{[A,B,Q,R,T]}$ such that
\begin{multline} \label{eq:mr_conv}
    \mathbb{E}[|\bm{x}_{R-M}(t) - \bm{x}^*_\infty(t)|] + \mathbb{E}[|\bm{u}_{R-M}(t) - \bm{u}_\infty^*(t)|] \leq \\C_{[A,B,Q,R,T]}\sqrt{h \mathrm{Var}[A_R]}\frac{e^{-\mu_{R-M}(t-\tau)}(t+1)}{\mu_{R-M}}|\bm{x}_0| \\
    +C_{[A,B,Q,R]}e^{-2\mu_\infty(T-\tau)} \frac{e^{-\mu_{R-M}t}(t+1)}{\mu_{R-M}}|\bm{x}_0|. 
\end{multline}
\end{corollary}

\begin{proof}
The estimate on the state follows because
\begin{multline}
\mathbb{E}[|\bm{x}_{R-M}(t) - \bm{x}^*_\infty(t)|] \leq \\ \mathbb{E}[|\bm{x}_{R-M}(t) - \bm{x}_M(t)|] + |\bm{x}_M(t) - \bm{x}^*_\infty(t)|. 
\end{multline}
A bound for the first term is given in Main Result \ref{res:main_2}. The bound for the second term follows from the analysis of MPC in \cite[Theorem 2.7]{mpc}, see also Lemma \ref{lem:MPCconv} in Section \ref{sec:preliminaries}. The bound on the controls follows analogously. 
\end{proof}

The form of the estimates in Main Results \ref{res:main_1} and \ref{res:main_2} clearly indicate a tuning strategy for the main parameters $T$, $\tau$, and $h$ in the RBM-MPC strategy. First, $T-\tau$ should be chosen such that the MPC strategy is stabilizing, i.e.\ such that $\mu_\infty - C_{\mathrm{MPC}}e^{-\mu_\infty(T-\tau)} > 0$. After that $h$ can be chosen such that $\mu_{R-M} > 0$ and such that the RBM-MPC strategy leads to a sufficiently good approximation of the MPC strategy. 


\section{Preliminaries} \label{sec:preliminaries}

\subsection{Notation}
For a symmetric matrix $M$, $M \succcurlyeq 0$ or $M \succ 0$ indicates that $M$ is positive semi-definite or positive definite, respectively. For two symmetric matrices $M$ and $N$, $M \succcurlyeq N$ means that $M - N \succcurlyeq 0$. The norm of a vector $\bm{x} \in \mathbb{R}^n$ is $|\bm{x}| = \sqrt{\bm{x}^\top\bm{x}}$. The operator norm of a matrix $A$ is defined as
\begin{equation}
\| A \| = \sup_{|\bm{x}| = 1} |A\bm{x}|. 
\end{equation}
Furthermore, $C_{[a,b,\ldots,c]}$ denotes a constant depending only on the parameters $a,b, \ldots, c$ that may vary from expression to expression. We may for example write $a^b C_{[a,b]} \leq C_{[a,b]}$. 

Let $V$ be a vector space. The expected value of a random variable $X : \{1,2,\ldots,2^M \}^K \rightarrow V$ depending on $\bm{\omega}_i$ is
\begin{equation}
     \mathbb{E}_i[X] = \sum_{\bm{\omega}_i \in \{1,2, \ldots, 2^M \}^K} X(\bm{\omega}_i) p(\bm{\omega}_i),
\end{equation}
where
\begin{equation}
     p(\bm{\omega}_i) = p_{\omega_{i,1}} p_{\omega_{i,2}} \cdots p_{\omega_{i,K}}. 
\end{equation}
The expected value of a random variable $X : (\{1,2,\ldots,2^M \}^K)^{i+1} \rightarrow V$ depending on $\Omega_i$ in \eqref{eq:def_Omegai} is denoted by
 \begin{multline}
     \mathbb{E}[X] = \sum_{\bm{\omega}_0, \bm{\omega}_1, \ldots , \bm{\omega}_i \in \{1,2, \ldots, 2^M \}^K} X(\bm{\omega}_0, \bm{\omega}_1, \ldots , \bm{\omega}_i) \times \\ p(\bm{\omega}_0) p(\bm{\omega}_1) \cdots p(\bm{\omega}_i).
 \end{multline}
For the expected value of a random variable $X(\Omega_i)$ w.r.t.\ only the last $\bm{\omega}_i$, we will write $\mathbb{E}_i[X(\Omega_{i-1})]$ to indicate that the result only depends on $\Omega_{i-1}$. 
 
For two random variables $X(\Omega_i)$ and $Y(\Omega_i)$, the Cauchy-Schwartz inequality in the probability space shows that 
\begin{equation}
\mathbb{E}[XY] \leq \sqrt{\mathbb{E}[X^2]\mathbb{E}[Y^2]}, \qquad \mathbb{E}[X] \leq \sqrt{\mathbb{E}[X^2]}.
\end{equation}
where the second equality follows after setting $Y = 1$. 
Similar expressions hold for the expectation $\mathbb{E}_i$ introduced above. 

\subsection{Riccati Theory} \label{ssec:Riccati}
This subsection recalls several well-known results from Riccati theory that can for example be found in \cite{sage1968}. 

First of all, the state trajectory $\bm{x}^*_\infty(t)$ minimizing the infinite horizon cost \eqref{eq:Jinf}--\eqref{eq:dyn_x} satisfies
\begin{equation}
    \dot{\bm{x}}_\infty^*(t) = A_\infty \bm{x}_\infty^*(t), \qquad \bm{x}^*_\infty(0) = \bm{x}_0, \label{eq:xinfty}
\end{equation}
with 
\begin{equation}
    A_\infty = A-BR^{-1}B^\top P_\infty, \label{eq:Ainfty}
\end{equation}
where $P_\infty$ is the solution of the (unique) symmetric positive-definite solution of the Algebraic Riccati Equation (ARE)
\begin{align}
    A^\top P_{\infty} + P_{\infty}A - P_{\infty}BR^{-1}B^\top P_{\infty} + Q = 0. \label{eq:ARE}
\end{align}
Comparing \eqref{eq:xinfty} and \eqref{eq:dyn_x}, it follows that
\begin{equation}
    \bm{u}^*_\infty(t) = -R^{-1}B^\top P_\infty \bm{x}^*_\infty(t). \label{eq:uinf}
\end{equation}
The matrix $P_\infty$ also has the property that 
\begin{equation}
    \min_{\bm{u}_\infty} J_\infty(\bm{u}_\infty) = J_\infty(\bm{u}_\infty^*) = \bm{x}_0^\top P_\infty \bm{x}_0. \label{eq:minJinfty}
\end{equation}
The observability of the pair $(A,Q)$ implies that $A_\infty$ is stable, i.e. there exist $M_\infty \geq 1$ and $\mu_\infty > 0$ such that
\begin{equation}
\| e^{A_\infty t} \| \leq M_{\infty} e^{-\mu_{\infty}t}. \label{eq:expon_stab_infty}
\end{equation}

Similarly, the state trajectory $\bm{x}_T^*(t; \bm{x_i}, \tau_i)$ minimizing the finite horizon cost \eqref{eq:JT}--\eqref{eq:dyn_xT} satisfies
\begin{equation}
    \dot{\bm{x}}^*_T(t) = (A - B R^{-1}B^\top P_T(t - \tau_i)) \bm{x}^*_T(t), \quad \bm{x}^*_T(\tau_i) = \bm{x}_i, \label{eq:xT}
\end{equation}
where $\bm{x}_T^*(t)$ denotes $\bm{x}_T^*(t; \bm{x_i}, \tau_i)$ for brevity and $P_T(t)$ is the solution of the Riccati Differential Equation (RDE)
\begin{align}
-\dot{P}_T(t) & = A^\top P_T(t) + P_T(t)A \nonumber \\ & - P_T(t)BR^{-1}B^\top P_T(t) + Q, \quad P_T(T) = 0, \label{eq:RDE}
\end{align}
on $[0,T]$. 
Comparing \eqref{eq:dyn_xT} and \eqref{eq:xT}, it follows that
\begin{equation}
    \bm{u}^*_T(t; \bm{x_i}, \tau_i) = -R^{-1}B^\top P_T(t - \tau_i) \bm{x}^*_T(t; \bm{x}_i, \tau_i), \label{eq:uTstar}
\end{equation}
and it holds that
\begin{equation}
    \min_{\bm{u}_T} J_T(\bm{u}_T; \bm{x}_i, \tau_i) = J_T(\bm{u}_T^*; \bm{x}_i, \tau_i) = \bm{x}_i^\top P_T(0) \bm{x}_i. \label{eq:minJT}
\end{equation}
Comparing \eqref{eq:minJT} with $\bm{x}_i = \bm{x}_0$ to \eqref{eq:minJinfty}, it follows that
\begin{equation}
    \bm{x}_0^\top P_T(0) \bm{x}_0 \leq J_T(\left. \bm{u}_\infty^* \right|_{[0,T]}; \bm{x}_0, 0) \leq \bm{x}_0^\top P_\infty \bm{x}_0, 
\end{equation}
where $\left. \bm{u}_\infty \right|_{[0,T]}$ denotes the restriction of $\bm{u}_\infty^*(t)$ to $[0,T]$. 
Since this equation holds for all $\bm{x}_0$, it follows that $0 \preccurlyeq P_T(0) \preccurlyeq P_\infty$, and thus that
\begin{equation}
    \| P_T(0) \| \leq \| P_\infty \|. \label{eq:est_PT0}
\end{equation}
Because the definition of $P_T(t)$ in \eqref{eq:RDE} shows that $P_T(t) = P_{T-t}(0)$, \eqref{eq:est_PT0} implies that (for all $t \in [0,T]$)
\begin{equation}
    \| P_T(t) \| \leq \| P_\infty \|. \label{eq:boundPT}
\end{equation}

Analogously, the state trajectory $\bm{x}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ that minimizes the finite horizon randomized cost \eqref{eq:JR}--\eqref{eq:dyn_xR} satisfies
\begin{multline}
    \dot{\bm{x}}^*_R(\bm{\omega}_i, t) = \\ (A_R(\bm{\omega}_i,t;\tau_i) - B R^{-1}B^\top P_R(\bm{\omega}_i, t - \tau_i)) \bm{x}^*_R(\bm{\omega}_i, t), \\ \bm{x}^*_R(\bm{\omega}_i, \tau_i) = \bm{x}_i, \label{eq:xR}
\end{multline}
where $\bm{x}_R^*(\bm{\omega}_i, t)$ denotes $\bm{x}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ for brevity and $P_R(\bm{\omega}_i, t)$ is the solution of the Randomized Riccati Differential Equation (RRDE)
\begin{align} \label{eq:RRDE}
&\dot{P}_R(\bm{\omega}_i, t) = 
A_R(\bm{\omega}_i, t; \tau_i)^\top P_R(\bm{\omega}_i, t) + \\ 
&\ P_R(\bm{\omega}_i, t)A_R(\bm{\omega}_i, t; \tau_i)
-P_R(\bm{\omega}_i,t)BR^{-1}B^\top P_R(\bm{\omega}_i, t) + Q, \nonumber
\end{align}
and the final condition $P_R(\bm{\omega}_i, T) = 0$ on $[0,T]$. Comparing \eqref{eq:dyn_xR} and \eqref{eq:xR}, it follows that
\begin{equation}
    \bm{u}^*_R(\bm{\omega}_i, t) = -R^{-1}B^\top P_R (\bm{\omega}_i, t - \tau_i) \bm{x}^*_R(\bm{\omega}_i, t). \label{eq:uRstar}
\end{equation}
where $\bm{u}_R^*(\bm{\omega}_i, t)$ denotes $\bm{u}_R^*(\bm{\omega}_i, t; \bm{x_i}, \tau_i)$ for brevity. It also holds that
\begin{multline}
    \min_{\bm{u}_T} J_R(\bm{\omega}_i, \bm{u}_T; \bm{x}_i, \tau_i) = J_R(\bm{\omega}_i, \bm{u}_R^*(\bm{\omega}_i); \bm{x}_i, \tau_i)  \\
    = \bm{x}_i^\top P_R(\bm{\omega}_i,0) \bm{x}_i. \label{eq:minJR}
\end{multline}


\subsection{Analysis of Model Predictive Control} \label{ssec:mpc_ana}
In this subsection, the main ideas from the analysis of MPC in \cite{mpc} are summarized. We note that several other analysis methods for MPC have been proposed, see, e.g., \cite{mayne2000, grune2017, rawlings2019}. 

A fundamental result is the convergence of the solution $P_T(t)$ of the RDE \eqref{eq:RDE} to the solution $P_\infty$ of ARE \eqref{eq:ARE}. 
\begin{lemma} \label{lem:Pconv}
If $(A,B)$ is controllable, $(A,Q)$ is observable, and $\mu_\infty$ is as in \eqref{eq:expon_stab_infty}, there exists a constant $C_{\mathrm{MPC}}'$ only depending on $A,B,Q$ and $R$ such that
    \begin{equation}\label{P_infty - P_T 0}
        ||P_T(t) - P_{\infty}|| \leq C_{\mathrm{MPC}}' e^{-2\mu_{\infty}(T-t)}.
    \end{equation}
\end{lemma}

\begin{proof}
See \cite[Lemma 2.2]{mpc}, \cite{porretta2013}, or \cite{callier1994}. 
\end{proof}

\begin{remark}
Lemma \ref{lem:Pconv} also holds under the weaker assumptions that $(A,B)$ is stabilizable and $(A,Q)$ is detectable. 
\end{remark}

The form of the optimal state trajectory $\bm{x}^*_T(t;\bm{x}_i, \tau_i)$ in \eqref{eq:xT} and the MPC algorithm from Subsection \ref{ssec:mpc} imply that
\begin{equation}
    \dot{\bm{x}}_{M}(t) = A_{T,\tau}(t) \bm{x}_M(t), \qquad \bm{x}_M(0) = \bm{x}_0. \label{eq:dyn_xMPC}
\end{equation}
where $A_{T,\tau}(t)$ is the $\tau$-periodic matrix
\begin{align}
    A_{T,\tau}(t) &= A - BR^{-1}B^\top P_{T,\tau}(t), \label{eq:def_ATtau} \\
    P_{T,\tau}(t) &= P_{T}(t \text{ mod } \tau). \label{eq:def_PTtau}
\end{align}

Lemma \ref{lem:Pconv} and \eqref{eq:dyn_xMPC} lead to the following two lemmas on the stability and convergence of MPC. The first lemma shows that MPC is stabilizing for $T-\tau$ sufficiently large. 

\begin{lemma} \label{lem:MPC}
Let $M_\infty \geq 1$ and $\mu_\infty > 0$ be as in \eqref{eq:expon_stab_infty} and let $C_{\mathrm{MPC}}'$ be as in Lemma \ref{lem:Pconv}, then
\begin{equation}
    |\bm{x}_M(t)| \leq M_\infty e^{-\mu_{T-\tau} t} |\bm{x}_0|, 
\end{equation}
where
\begin{align}
    \mu_{T-\tau} &= \mu_\infty - C_{\mathrm{MPC}}e^{-2\mu_\infty (T-\tau)}, \\
    C_{\mathrm{MPC}} &= M_\infty \| BR^{-1}B^\top \| C_{\mathrm{MPC}}'. \label{eq:CMPC}
\end{align}
\end{lemma}

\begin{proof}
See \cite[Theorem 2.6]{mpc}. 
\end{proof}

The second lemma shows that $\bm{x}_M(t)$ and $\bm{u}_M(t)$ converge to $\bm{x}^*_\infty(t)$ and $\bm{u}^*_\infty(t)$ for $T-\tau \rightarrow \infty$. 

\begin{lemma} \label{lem:MPCconv}
There exists a constant $C_{[A,B,Q,R]}$ such that
\begin{multline}
    |\bm{x}_M(t) - \bm{x}_\infty^*(t)| + |\bm{u}_M(t) - \bm{u}^*_\infty(t)| \leq \\
    C_{[A,B,Q,R]} e^{-2\mu_\infty (T-\tau)} e^{-\mu_{T-\tau}} (t+1) |\bm{x}_0|. 
\end{multline}
\end{lemma}

\begin{proof}
See \cite[Theorem 2.7]{mpc}. 
\end{proof}

\subsection{Convergence of the Random Batch Method}\label{RBM}

This subsection contains the important results from \cite{rbm} that will be used in the analysis of RBM-MPC in the next section. 

Let $S_R(\bm{\omega}_i, t,s; \tau_i)$ (for $\tau_i \leq s \leq t \leq \tau_i + T$) denote the semigroup generated by $A_R(\bm{\omega}_i,t; \tau_i)$, i.e. $S_R(\bm{\omega}_i, t,s; \tau_i)\bm{x}= \bm{x}(\bm{\omega}_i, t)$ where $\bm{x}(\bm{\omega}_i, t)$ satisfies
\begin{equation}
\dot{\bm{x}}(\bm{\omega}_i,t) = A_R(\bm{\omega}_i, t; \tau_i) \bm{x}(\bm{\omega}_i,t), \qquad  \bm{x}(\bm{\omega}_i,s) = \bm{x}. 
\end{equation}

The following lemma from \cite{rbm} then shows that $S_R(\bm{\omega}_i,t,s)$ is (in expectation) close to $e^{A(t-s)}$ when $h \mathrm{Var}[A_R]$ is small. 
\begin{lemma} \label{lem:semigroup}
Let $\mathrm{Var}[A_R]$ and $\mu_R$ be as in \eqref{eq:def_varA} and \eqref{eq:def_muR}, then
\begin{equation}
\mathbb{E}_i[\| S_R(t,s) - e^{A(t-s)} \|^2 ] \leq h \mathrm{Var}[A_R]f(t-s) .
\end{equation}
where 
\begin{equation}
f(t) = (2t + \|A \|t^2)e^{2\mu_Rt}. \label{eq:def_f}
\end{equation}
\end{lemma}
\begin{proof}
See \cite[Theorem 1 and Corollary 1]{rbm}. 
\end{proof}

With Lemma \ref{lem:semigroup}, it is also possible to bound the difference between controlled state trajectories. In particular, let $\bm{u}_R: \{1,2, \ldots, 2^M \}^K \times [\tau_i, \tau_i + T] \rightarrow \mathbb{R}^q$ be a random control and let $\bm{x}_R(\bm{\omega}_i,t)$ and $\bm{y}_R(\bm{\omega}_i,t)$ be the state trajectories satisfying
\begin{align}
\dot{\bm{x}}_R(\bm{\omega}_i, t) &= A\bm{x}_R(\bm{\omega}_i, t) + B \bm{u}_R(\bm{\omega}_i,t), \label{eq:lem_dynxR} \\
\dot{\bm{y}}_R(\bm{\omega}_i, t) &= A_R(\bm{\omega}_i,t; \tau_i) \bm{y}_R(\bm{\omega}_i, t) + B \bm{u}_R(\bm{\omega}_i,t), \label{eq:lem_dynyR} \\
\bm{x}_R(\bm{\omega}_i,\tau_i) &= \bm{y}_R(\bm{\omega}_i,\tau_i) = \bm{x}_i(\bm{\omega}_i), \label{eq:lem_xRyRinit}
\end{align}
where $\bm{x}_i : \{1,2, \ldots, 2^M \}^K \rightarrow \mathbb{R}^n$ is a random initial condition. We then have the following lemma. 
\begin{lemma} \label{lem:xRyR}
Let $f(t)$ be as in Lemma \ref{lem:semigroup} and let $\bar{X}$ and $\bar{U}$ be such that
\begin{equation}
|\bm{x}_i(\bm{\omega}_i)| \leq \bar{X}, \qquad |\bm{u}_R(\bm{\omega}_i)|_{L^2(\tau_i, \tau_i+T)} \leq \bar{U}, \label{eq:def_barU}
\end{equation}
for all $\bm{\omega}_i \in \{1,2, \ldots, 2^M \}^K$, then
\begin{multline}
\mathbb{E}[|\bm{y}_R(t;\bm{x}_i, \tau_i) - \bm{x}_R(t;\bm{x}_i, \tau_i)|^2] \leq \\
f(t-\tau_i)h \mathrm{Var}[A_R] \left(\bar{X} + \| B \| \sqrt{t-\tau_i} \bar{U} \right)^2. 
\end{multline}
\end{lemma}
\begin{proof}
By a slight modification of \cite[Theorem 2]{rbm} in which the random initial condition was not considered. 
\end{proof}

Finally, we present a sharpened version of the no-gap condition from \cite[Theorem 3]{rbm} which bounds the difference between the minimal values of $J_R(\bm{\omega}_i, \cdot; \bm{x}_i, \tau_i)$ in \eqref{eq:JR} and $J_T(\cdot; \bm{x}_i, \tau_i)$ in \eqref{eq:JT}. Because these minimal values are related to the Riccati operators $P_R(\bm{\omega}_i,t)$ and $P_T(t)$ by \eqref{eq:minJR} and \eqref{eq:minJT}, the following result can be obtained.  
\begin{lemma} \label{lem:PRPT}
Let $P_R(\bm{\omega}_i,t)$ and $P_T(t)$ be the solutions of the RDEs \eqref{eq:RRDE} and \eqref{eq:RDE}, respectively, then 
\begin{align}
&\mathbb{E}_i[\|P_R(t) - P_T(t)\|] \leq \label{eq:lem_PRPT} \\
&C_{[A,B,Q,R]} \left( \sqrt{T^2f(T)h\mathrm{Var}[A_R]} + T^2 f(T) h \mathrm{Var}[A_R] \right). \nonumber
\end{align}
\end{lemma}
The proof of Lemma \ref{lem:PRPT} is in Appendix \ref{app:lem_PRPT}. 

\section{Proofs} \label{sec:proofs}

\subsection{Proof of Main Result \ref{res:main_1}}
Let $i$ be such that $t \in [\tau_{i}, \tau_{i+1})$. The RBM-MPC algorithm in Subsection \ref{ssec:rbmmpc} shows that $\bm{x}_{R-M}(\Omega_{i}, t)$ is $\bm{y}_R(\bm{\omega}_{i},t; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_{i})$. Because $\bm{y}_R(\bm{\omega}_i, t)$ satisfies \eqref{eq:dyn_yR}, it follows that
\begin{equation}
    \dot{\bm{x}}_{R-M}(\Omega_i,t) = A\bm{x}_{R-M}(\Omega_i,t; \tau_i) + B\bm{u}_R^*(\bm{\omega}_i,t), \label{eq:dyn_xRM}
\end{equation}
where $\bm{u}_R^*(\bm{\omega}_i,t; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_{i})$ is replaced by $\bm{u}_R^{*}(\bm{\omega}_{i},t)$ for brevity. The RHS of \eqref{eq:dyn_xRM} can be rewritten as
\begin{equation}
\dot{\bm{x}}_{R-M}(\Omega_{i},t) = A_\infty\bm{x}_{R-M}(\Omega_{i},t) + \bm{r}(\Omega_i, t), \label{eq:dyn_xRM1}
\end{equation}
with $A_\infty$ as in \eqref{eq:Ainfty} and the residual is given by
\begin{equation}
\bm{r}(\Omega_i,t) = B R^{-1}B^\top P_\infty \bm{x}_{R-M}(\Omega_i,t) + B \bm{u}^*_R(\bm{\omega}_i,t). \label{eq:def_r}
\end{equation}
The residual $\bm{r}(\Omega_i,t)$ can again be split into two parts
\begin{multline}
\bm{r}(\Omega_i,t) = B R^{-1}B^\top (P_\infty - P_T(t-\tau_i)) \bm{x}_{R-M}(\Omega_i,t) +  \\
B \bm{g}(\Omega_i,t), \label{eq:r_split}
\end{multline}
where 
\begin{equation}
\bm{g}(\Omega_i,t) = R^{-1}B^\top P_T(t-\tau_i) \bm{x}_{R-M}(\Omega_i,t) + \bm{u}^*_R(\bm{\omega}_i,t). \label{eq:def_g}
\end{equation}
The first term in \eqref{eq:r_split} represents the error introduced by MPC and can be bounded using Lemma \ref{lem:Pconv}. The following lemma bounds the expected value of $|\bm{g}(\Omega_i,t)|$ w.r.t.\ $\bm{\omega}_i$. 
\begin{lemma} \label{lem:bound_g}
Let $\bm{g}(\Omega_i,t)$ be as in \eqref{eq:def_g}, then there exists a constant $C_{[A,B,Q,R,T]}$ such that
\begin{multline}
\mathbb{E}_i[ |\bm{g}(\Omega_{i-1}, t)| ] \leq \\
C_{[A,B,Q,R,T]} \sqrt{h \mathrm{Var}[A_R]} |\bm{x}_{R-M}(\Omega_{i-1}, \tau_{i})|. 
\end{multline}
\end{lemma}
The proof of Lemma \ref{lem:bound_g} is postponed until Subsection \ref{ssec:proof_lem_bound_g}. 

Combining this result with Lemma \ref{lem:Pconv}, it follows that
\begin{multline}
\mathbb{E}_i[ |\bm{r}(\Omega_{i-1}, t)| ] \leq \frac{C_{\mathrm{MPC}}}{M_\infty} e^{-2\mu_\infty (T-\tau)} \mathbb{E}_i[|\bm{x}_{R-M}(\Omega_{i-1},t)|] \\
+ C_{[A,B,Q,R,T]} \sqrt{h \mathrm{Var}[A_R]} |\bm{x}_{R-M}(\Omega_{i-1}, \tau_{i})|.  \label{eq:r_bound}
\end{multline}

Main Result \ref{res:main_1} can now be proven as follows. Applying the variation of constants formula to \eqref{eq:dyn_xRM1} shows that
\begin{equation}
\bm{x}_{R-M}(\Omega_i,t) = e^{A_\infty t} \bm{x}_0 + \int_0^t e^{A_\infty (t-s)} \bm{r}(\Omega_i,s) \ \mathrm{d}s. 
\end{equation}
Taking the norm using \eqref{eq:expon_stab_infty} and then the expectation yields
\begin{multline}
\mathbb{E}[|\bm{x}_{R-M}(t)|] = M_\infty e^{-\mu_\infty t} \bm{x}_0 \\
+ M_\infty \int_0^t e^{-\mu_\infty (t-s)} \mathbb{E}[|\bm{r}(s)|] \ \mathrm{d}s. \label{eq:bound_xRM_step1}
\end{multline}
Write $i(s) = \lfloor s / \tau \rfloor$, so that $s \in [\tau_{i(s)}, \tau_{i(s)+1})$. By definition of the expectation
\begin{equation}
\mathbb{E}[|\bm{r}(s)|] = \mathbb{E}_1[\mathbb{E}_2[ \cdots \mathbb{E}_{i(s)-1}[\mathbb{E}_{i(s)}[|\bm{r}(s)|]]\cdots  ] ]. 
\end{equation}
Using \eqref{eq:r_bound} to bound the inner expectation and then applying the other expectation operators to the resulting estimate, it follows that
\begin{multline}
\mathbb{E}[|\bm{r}(s)|] \leq \frac{C_{\mathrm{MPC}}}{M_\infty}e^{-2\mu_\infty (T-\tau)} \mathbb{E}[|\bm{x}_{R-M}(t)|] \\
+ C_{[A,B,Q,R,T]} \sqrt{h \mathrm{Var}[A_R]} \mathbb{E}[|\bm{x}_{R-M}(\tau_{i(s)})|]. \label{eq:bound_expectation_g}
\end{multline}
After inserting this estimate back into \eqref{eq:bound_xRM_step1}, Main Result \ref{res:main_1} follows from the following variation of Gronwall's lemma. 
Note that $\tau \leq 1$ (see Remark \ref{rem:asymptotics}) so that the factor $e^{\mu_\infty \tau}$ that appears can be estimated as $e^{\mu_\infty \tau} \leq e^{\mu_\infty} \leq C_{[A,B,Q,R]}$. 
\begin{lemma} \label{lemma}
Suppose that $f(t) \geq 0$ and that there exist constants $\mu, C_1, C_2, C_3 \geq 0$ such that for all $t$
\begin{align}\label{assumption}
\begin{split}
    f(t) & \leq C_1e^{-\mu t} \\ & \;\;\; + \int_0^t e^{-\mu (t-s)}\bigg(C_2 f(\tau_{i(s)}) + C_3 f(s) \bigg)\,ds,
\end{split}
\end{align}
then
\begin{align}
    f(t) \leq C_1e^{(-\mu + C_2e^{\mu\tau} + C_3)\,t}.
\end{align}
\end{lemma}
The proof of Lemma \ref{lemma} can be found in Appendix \ref{app:lemma}. 

\subsection{Proof of Main Result \ref{res:main_2}}
Let $i$ again be such that $t \in [\tau_{i}, \tau_{i+1})$ and write
\begin{equation}
    \bm{e}_{R-M}(\Omega_i, t) := \bm{x}_{R-M}(\Omega_i, t) - \bm{x}_M(t). \label{eq:def_eRM}
\end{equation}
Note that inserting \eqref{eq:r_split} into \eqref{eq:dyn_xRM1} and using the definition of $A_{T,\tau}(t)$ in \eqref{eq:def_ATtau} shows that
\begin{equation}
\dot{\bm{x}}_{R-M}(\Omega_i,t) = A_{T,\tau}(t) \bm{x}_{R-M}(\Omega_i,t) + B\bm{g}(\Omega_i,t). 
\end{equation}
Subtracting \eqref{eq:dyn_xMPC} from this equation shows that
\begin{equation}
\dot{\bm{e}}_{R-M}(\Omega_i,t) = A_{T,\tau}(t) \bm{e}_{R-M}(\Omega_i,t) + B\bm{g}(\Omega_i,t),
\end{equation}
and $\bm{e}_{R-M}(\Omega_i,0) = 0$. Applying the variation of constants formula yields
\begin{equation}
\bm{e}_{R-M}(\Omega_i,t) = \int_0^t e^{\int_s^t A_{A,\tau}(\sigma) \ \mathrm{d}\sigma} B\bm{g}(\Omega_i,s) \ \mathrm{d}s. \label{eq:bound_eRM_step0}
\end{equation}
Note that Lemma \ref{lem:MPC} implies that
\begin{equation}
\| e^{\int_s^t A_{A,\tau}(\sigma) \ \mathrm{d}\sigma} \| \leq M_\infty e^{-\mu_{T-\tau} (t-s)} \leq M_\infty,
\end{equation}
where the second inequality uses that $\mu_{T-\tau} > \mu_{R-M} > 0$. 
Taking norms and the expectation in \eqref{eq:bound_eRM_step0} thus yields
\begin{equation}
\mathbb{E}[|\bm{e}_{R-M}(t)|] \leq M_\infty \| B \| \int_0^t \mathbb{E}[|\bm{g}(s)|] \ \mathrm{d}s.  \label{eq:bound_eRM_step1}
\end{equation}
Using the bound from Main Result \ref{res:main_1} in the estimate from Lemma \ref{lem:bound_g} furthermore shows that
\begin{align}
\mathbb{E}[|\bm{g}(s)|] &\leq C_{[A,B,Q,R,T]} \sqrt{h \mathrm{Var}[A_R]} e^{\mu_{R-M} \tau} e^{-\mu_{R-M}s} |\bm{x}_0| \nonumber \\
&\leq C_{[A,B,Q,R,T]} \sqrt{h \mathrm{Var}[A_R]} e^{-\mu_{R-M}s} |\bm{x}_0|, \label{eq:bound_expectation_g2}
\end{align}
where the second estimate follows because $\tau \leq 1$ by assumption (see Remark \ref{rem:asymptotics}) so that $e^{\mu_{R-M} \tau} \leq e^{\mu_\infty}$. 
The estimate for the states now follows by inserting \eqref{eq:bound_expectation_g2} into \eqref{eq:bound_eRM_step1} and evaluating the remaining integral. 

For the bound on the controls, note that \eqref{eq:dyn_xMPC} shows that
\begin{align}
    \bm{u}_{R-M}(& \Omega_i, t) - \bm{u}_M(t) = \bm{u}_R^*(\bm{\omega}_i,t) + R^{-1}B^\top P_{T,\tau}(t) \bm{x}_M(t) \nonumber \\
    &=\bm{g}(\Omega_i,t) - R^{-1}B^\top P_{T,\tau}(t) \bm{e}_{R-M}(\Omega_i,t). 
\end{align}
where $\bm{u}^*_R(\bm{\omega}_i, t)$ denotes $\bm{u}^*_R(\bm{\omega}_i, t;\bm{x}_{R-M}(\Omega_{i-1},\tau_i), \tau_i)$ as before. 
The bound on the controls follows after taking the norm and the expected value using Lemma \ref{lem:bound_g} and the bound for $\mathbb{E}[|\bm{e}_{R-M}(t)|]$ obtained in the first part of the proof. 

\subsection{Proof of Lemma \ref{lem:bound_g}} \label{ssec:proof_lem_bound_g}
Let $i$ again by such that $t \in [\tau_i, \tau_{i+1})$. 
Note that $\bm{g}(\Omega_i, t)$ in \eqref{eq:def_g} can be rewritten as
\begin{multline}
    \bm{g}(\Omega_i,t) = R^{-1}B^\top \underbrace{P_T(t-\tau_{i})}_{\eqref{eq:boundPT}}\underbrace{(\bm{x}_{R-M}(\Omega_{i},t)-\bm{x}^*_R(\bm{\omega}_{i},t))}_{\mathrm{Lemma}\ \ref{lem:xRyR}} \\
    + R^{-1}B^\top\underbrace{(P_T(t-\tau_{i}) - P_R(\bm{\omega}_{i}, t-\tau_{i}))}_{\mathrm{Lemma}\ \ref{lem:PRPT}}\underbrace{\bm{x}^*_R(\bm{\omega}_{i},t)}_{\eqref{eq:def_muR}}, \label{eq:lem_g_step1}
\end{multline}
where \eqref{eq:uRstar} has been used to express $\bm{u}^*_R(\bm{\omega}_{i}, t)$ and $\bm{x}^*_R(\bm{\omega}_i, t)$ denotes $\bm{x}_R^*(\bm{\omega}_i, t; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_{i})$ for brevity. The main lemmas and equations are highlighted in \eqref{eq:lem_g_step1}. After taking norms and the expectation w.r.t.\ $\bm{\omega}_i$ in \eqref{eq:lem_g_step1}, it follows that
\begin{multline}
    \mathbb{E}_i[|\bm{g}(\Omega_{i-1},t)|] \leq \\
    C_{[A,B,Q,R]} \mathbb{E}_i[|\bm{x}_{R-M}(\Omega_{i-1},t)-\bm{x}^*_R(t)|] \\
     +\| R^{-1}B^\top \| \mathbb{E}_i[\| P_T(t-\tau_{i}) - P_R(t-\tau_{i})\|] \\ 
     \times \max_{\bm{\omega}_i} |\bm{x}^*_R(\bm{\omega}_{i},t)|. \label{eq:lem_g_step3}
\end{multline}
The bound for $\mathbb{E}_i[\| P_T(t-\tau_{i}) - P_R(t-\tau_{i})\|]$ is given in Lemma \ref{lem:PRPT}. 
It remains to find bounds for $\mathbb{E}_i[|\bm{x}_{R-M}(\Omega_{i-1},t)-\bm{x}^*_R(t)|]$ and $\max_{\bm{\omega}_i} |\bm{x}^*_R(\bm{\omega}_{i},t)|$. 

For $\mathbb{E}_i[|\bm{x}_{R-M}(\Omega_{i-1},t)-\bm{x}^*_R(t)|]$, apply Lemma \ref{lem:xRyR} with $\bm{u}_R(\bm{\omega}_i,t) = \bm{u}_R^*(\bm{\omega}_i,t; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_{i})$ and the initial condition $\bm{x}_i(\bm{\omega}_i) = \bm{x}_{R-M}(\Omega_{i-1},\tau_i)$. This makes
\begin{align*}
\bm{y}_R(\bm{\omega}_i,t) &= \bm{x}^*_R(\bm{\omega}_i,t; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_{i}), \\
\bm{x}_R(\bm{\omega}_i,t) &= \bm{y}_R^*(\bm{\omega}_i,t; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_{i}) = \bm{x}_{R-M}(\Omega_i,t).
\end{align*}
Clearly, $\bar{X} = |\bm{x}_{R-M}(\Omega_{i-1},\tau_i)|$ does not depend on $\bm{\omega}_i$. 
It remains to find $\bar{U}$. Because $R$ is positive definite and $\bm{u}^*_R(\bm{\omega}_i)$ minimizes $J_R(\bm{\omega}_i, \cdot)$, see \eqref{eq:minJR}, it follows that
\begin{align}
|\bm{u}^*_R(\bm{\omega}_i)|_{L^2}^2 &\leq C_{[R]} J_R(\bm{\omega}_i, \bm{u}_R^*(\bm{\omega}_i); \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_i) \nonumber \\ 
&\leq C_{[R]} J_R(\bm{\omega}_i, 0; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_i) \nonumber \\
&\leq C_{[Q,R]} |\bm{x}_R(\bm{\omega}_i; \bm{x}_{R-M}(\Omega_{i-1}, \tau_{i}), \tau_i)|_{L^2}^2. \label{eq:boundJR}
\end{align}
where $\bm{x}_R(\bm{\omega}_i,t; \bm{x}_i, \tau_i)$ satisfies \eqref{eq:dyn_xR} with $\bm{u}_R(t) = 0$. 
Now \eqref{eq:def_muR} shows that
\begin{multline}
\frac{d}{dt} |\bm{x}_R(\bm{\omega}_i,t)|^2 = 2\bm{x}_R^\top(\bm{\omega}_i,t) A_R(\bm{\omega}_i,t; \tau_i) \bm{x}_R(\bm{\omega}_i,t) \\ \leq 2\mu_R |\bm{x}_R(\bm{\omega}_i,t)|^2, 
\end{multline}
and it follows that
\begin{equation}
|\bm{x}_R(\bm{\omega}_i,t; \bm{x}_i, \tau_i)|^2 \leq e^{2\mu_R(t-\tau_i)}|\bm{x}_i|^2 \leq e^{2\mu_R T}|\bm{x}_i|^2.
\end{equation}
Inserting this bound into \eqref{eq:boundJR} yields
\begin{equation}
|\bm{u}^*_R(\bm{\omega}_i)|_{L^2(\tau_i, \tau_i+T)}^2 \leq C_{[Q,R]} T e^{2 \mu_R T} |\bm{x}_i|^2, \label{eq:bound_uRstar}
\end{equation}
which shows that $\bar{U}=\sqrt{T} e^{\mu_R T} |\bm{x}_{R-M}(\Omega_{i-1},\tau_i)|$. Lemma \ref{lem:xRyR} and the fact that $\mathbb{E}_i[|X|] \leq \sqrt{\mathbb{E}_i[|X|^2]}$ thus show that
\begin{multline}
\mathbb{E}_i[|\bm{x}_{R-M}(\Omega_{i-1},t)-\bm{x}^*_R(t)|] \leq \\ \sqrt{f(\tau) h \mathrm{Var}[A_R]} C_{[A,B,Q,R,T]} e^{\mu_RT} |\bm{x}_{R-M}(\Omega_{i-1},\tau_i)|, \label{eq:bound_xRMxR}
\end{multline}
for $t \in [\tau_i, \tau_{i+1}]$. 
This estimate is proportional to $\sqrt{\tau}$ because for $\tau \leq 1$, $f(\tau) \leq \tau (2 + \| A \|)e^{2\mu_R}$, see Remark \ref{rem:asymptotics} and \eqref{eq:def_f}. 

For the bound on $|\bm{x}^*_R(\bm{\omega}_{i},t)|$, note that \eqref{eq:def_muR} implies that
\begin{multline}
\frac{d}{dt}|\bm{x}_R^*(\bm{\omega}_i, t)| 
= \frac{(\bm{x}_R^*(\bm{\omega}_i,t))^\top \dot{\bm{x}}_R^*(\bm{\omega}_i,t)}{|\bm{x}_R^*(\bm{\omega}_i,t)|} \\
\leq \mu_R |\bm{x}^*_R(\bm{\omega}_i,t)| + \| B \| |\bm{u}^*_R(\bm{\omega}_i,t)|.
\end{multline}
Integrating the above inequality from $\tau_i$ to $t$ and applying Gronwall's lemma shows that
\begin{multline}
|\bm{x}^*_R(\bm{\omega}_i,t)| \leq e^{\mu_R (t-\tau_i)}  |\bm{x}_{R-M}(\Omega_{i-1},\tau_i)| + \\
e^{\mu_R (t-\tau_i)} \| B \|\int_{\tau_i}^t |\bm{u}^*_R(\bm{\omega}_i, s)| \ \mathrm{d}s. \label{eq:lem_xR_bound1}
\end{multline}
Since $t \in [\tau_i, \tau_{i+1}]$ it follows that
\begin{multline}
\int_{\tau_i}^t |\bm{u}^*_R(\bm{\omega}_i,s)| \ \mathrm{d}s \leq \sqrt{\tau} |\bm{u}^*_R(\bm{\omega}_i)|_{L^2(\tau_i,\tau_i+T)} \\ 
\leq C_{[Q,R]} \sqrt{\tau T} e^{\mu_RT} |\bm{x}_{R-M}(\Omega_{i-1}, \tau_i)|, \label{eq:lem_xR_bound2}
\end{multline}
where the second inequality follows from \eqref{eq:bound_uRstar}. Inserting \eqref{eq:lem_xR_bound2} into \eqref{eq:lem_xR_bound1} now yields
\begin{multline}
|\bm{x}_R^*(\bm{\omega}_i,t)| \leq \\ 
\left(1 + C_{[B,Q,R]} \sqrt{\tau T} e^{\mu_RT} \right) e^{\mu_R\tau} |\bm{x}_{R-M}(\Omega_{i-1},\tau_i)| \\
\leq C_{[A,B,Q,R,T]}e^{\mu_RT}|\bm{x}_{R-M}(\Omega_{i-1},\tau_i)|. \label{eq:bound_xR}
\end{multline}
The result follows by inserting \eqref{eq:bound_xRMxR}, \eqref{eq:lem_PRPT}, and \eqref{eq:bound_xR} into \eqref{eq:lem_g_step3}. 

\section{Numerical Example} \label{sec:numerical}
The proposed RBM-MPC algorithm is applied to a heat equation on the unit circle $S^1$
\begin{equation}
    \frac{\partial z}{\partial t}(\theta,t) = 0.001\frac{\partial^2 z}{\partial \theta^2}(\theta,t) + b(\theta) u(t),
\end{equation}
where $\theta \in [0,1]$ and $\theta = 0$ is identified with $\theta  = 1$, $b(\theta)$ is $1/W$ for $\theta \in [0,W/2] \cup [1-W/2,1]$ and zero otherwise, and $u(t)$ is the control. We fix $W = 0.1$. The initial condition is
\begin{equation}
z(\theta,0) = 10 \exp\left(-2(1 + \sin(2 \pi \theta)) \right).
\end{equation}
Our objective is to minimize the functional
\begin{equation}
\mathcal{J}(u) = \frac{100}{2}\int_0^\infty \int_{S^1} (z(\theta,t))^2 \ \mathrm{d}\theta \ \mathrm{d}t + \frac{1}{2}\int_0^\infty (u(t))^2 \ \mathrm{d}t. 
\end{equation}

The unit circle $S^1$ is discretized with a uniform grid with $n$ grid points by finite differences. This leads to a problem of the form \eqref{eq:Jinf}--\eqref{eq:dyn_x} with
\begin{equation}
A = \frac{1}{\Delta\theta^2} \begin{bmatrix}
-2 & 1 & 0 & \cdots  & 0 & 1 \\
1 & -2 & 1 &  & 0 & 0 \\
0 & 1 & -2 &  & 0 & 0 \\
\vdots & & & \ddots & \\
0 & 0 & 0 &  & -2 & 1 \\
1 & 0 & 0 &  & 1 & -2 
\end{bmatrix}, 
\end{equation}
$B \in \mathbb{R}^{n \times 1}$ contains the nodal values of $b(\theta)$, $Q = \Delta \theta I_n$, $R = 1$, $\Delta \theta = 1/n$, and $I_n$ denotes the $n \times n$ identity matrix. 

Time is discretized by the Crank-Nicholson scheme on the truncated time interval $[0,200]$ with a fixed time step $\Delta t$. At every time step, a linear system involving the matrix $A - \tfrac{\Delta t}{2}I_n$ needs to solved. Because this matrix is not tridiagonal, the computational cost per time step is of order $O(n^3)$. The OCPs that appear in MPC and RBM-MPC are solved by a steepest descent algorithm. The gradients are computed based on the discretization of the adjoint state from \cite{Apel}. The stepsize minimizes the functional in the direction of the gradient and optimization procedure is stopped when the relative change in the control is below $10^{-5}$ or after 1000 iterations. 

\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{Fig2}
\caption{Spatial discretization of the unit circle. Removing one connection between two grid points transforms the unit circle into the line segment at the bottom. }
\label{fig:torus}
\end{figure}

To construct the randomized matrix $A_R(\bm{\omega}_i,t;\tau_i)$, note that $A$ can be written as the sum of $M = n$ interconnection matrices as in \eqref{eq:Asplit}, where the first $n-1$ interconnection matrices $A_m$ are zero except for a diagonal block of the form
\begin{equation}
\frac{1}{\Delta \theta^2} \begin{bmatrix}
-1 & 1 \\ 1 & -1
\end{bmatrix}, 
\end{equation}
and the last interconnection matrix is
\begin{equation}
A_n = \frac{1}{\Delta \theta^2}\begin{bmatrix}
-1 & 0 & \cdots & 0 & 1 \\
0 &   &        &   & 0 \\
\vdots & &     & & \vdots \\
0 &   &        &   & 0 \\ 
1 & 0 & \cdots & 0 & -1
\end{bmatrix}.
\end{equation}
Observe that adding up the first $n-1$ interconnection matrices leads to a tridiagonal matrix
\begin{equation}
\sum_{m=1}^{n-1}A_m = \frac{1}{\Delta\theta^2} \begin{bmatrix}
-1 & 1 & 0 & \cdots  & 0 & 0 \\
1 & -2 & 1 &  & 0 & 0 \\
0 & 1 & -2 &  & 0 & 0 \\
\vdots & & & \ddots & \\
0 & 0 & 0 &  & -2 & 1 \\
0 & 0 & 0 &  & 1 & -1 
\end{bmatrix}. \label{eq:Atridiag}
\end{equation}
Therefore, the computational cost for inverting the matrix $\sum_{m=1}^{n-1}A_m - \tfrac{\Delta t}{2}I_n$ is of order $O(n)$. In fact, the rotational symmetry implies that removing any of the interconnection matrices $A_m$ leads to a tridiagonal matrix of the form \eqref{eq:Atridiag} up to a permutation of the rows and columns. Randomly removing one of the interconnection matrices thus reduces the computational cost per time step from $O(n^3)$ to $O(n)$. The randomized matrix $A_R(\bm{\omega}_i, t; \tau_i)$ is constructed by assigning a probability $1/n$ to all subsets of $\{1,2, \ldots, n \}$ of size $n-1$. The probabilities $\pi_m$ that $m$ is an element of the selected subset is thus
\begin{equation}
\pi_m = \frac{n-1}{n}. 
\end{equation}
The grid spacing $h$ is chosen as small as possible, so $h = \Delta t$. 
For a given prediction horizon $T$ and control horizon $\tau$, the randomized matrix $A_R(\bm{\omega}_i,t;\tau_i)$ can now be defined as in \eqref{eq:def_AR}.
\begin{remark}
Note that $h$ can in principle be chosen as $N_h\Delta t$ with $N_h \in \mathbb{N}_{\geq 2}$, but that this is not advisable because it will not reduce the computational cost and increase the RBM-error.  
\end{remark}

Figure \ref{fig:results_time_control} compares 20 realizations of the RBM-MPC control $u_{R-M}(\Omega_i,t)$ to the MPC control $u_M(t)$ and the infinite horizon control $u^*_\infty(t)$ for $n = 100$ spatial grid points. As can be seen, $u^*_\infty(t)$ is smooth, the MPC control $u_M(t)$ jumps when $t$ is a multiple of $\tau = 10$, and the RBM-MPC controls contain high-frequent oscillations related to the grid spacing $\Delta t = h = 1$. Figure \ref{fig:results_time_state} shows that that despite the relatively large deviations of $u_{R-M}(\Omega_i,t) $ from $u_M(t)$, the norm of the resulting state trajectory $|\bm{x}_{R-M}(\Omega_i,t)|$ is very close to the norm of the MPC state trajectory $|\bm{x}_M(t)|$ for all 20 considered realizations $\Omega_i$. The RBM-MPC control thus leads to practically the same convergence rate as the MPC control in this example. Also note that $T = 15$ is not much larger than $\tau$ in this example, but that the MPC strategy is stabilizing. 

\begin{figure*}
\centering
\subfloat[The controls $u_{R-M}(\Omega_i,t)$, $u_M(t)$, and $u^*_\infty(t)$. \label{fig:results_time_control}]{%
\includegraphics[width=0.48\textwidth]{Fig3a}}
\hfill
\subfloat[The norm of the state trajectories $\bm{x}_{R-M}(\Omega_i,t)$, $\bm{x}_M(t)$, and $\bm{x}^*_\infty(t)$. \label{fig:results_time_state}]{%
\includegraphics[width=0.48\textwidth]{Fig3b}}
\caption{The RBM-MPC control and state trajectory $u_{R-M}(\Omega_i,t)$ and $\bm{x}_{R-M}(\Omega_i, t)$ for 20 realizations of $\Omega_i$ compared to the MPC control and state trajectory $u_M(t)$ and $\bm{x}_M(t)$ and the infinite horizon control and state trajectory $u^*_\infty(t)$ and $\bm{x}^*_\infty(t)$ for $n=100$, $h=1$, $\tau = 10$, and $T = 15$. The lines for $|\bm{x}_{R-M}(\Omega_i,t)|$ and $|\bm{x}_M(t)|$ in Figure \ref{fig:results_time_state} almost overlap. 
}
\label{fig:results_time}
\end{figure*}

Table \ref{tab:times_varyn} shows that RBM-MPC significantly reduces the computational cost compared to MPC, which is in turn more efficient than solving the OCP on $[0,200]$ directly. The numbers between round brackets in Table \ref{tab:times_varyn} indicate the estimated standard deviation of the running times based on 20 runs. For $n = 100$, MPC is almost 3 times faster than a classical optimal control approach, and RBM-MPC is again almost 3 times faster than MPC. For $n = 1000$, MPC is still approximately 3 times faster than solving the OCP directly, but RBM-MPC is 5 times faster than MPC. Note that the relative speed up of RBM-MPC compared to MPC is not of $O(n^2)$ as the theoretical estimates predict because it also includes overheads and because solving the RBM-constrained OCP sometimes requires a few more iterations than the original OCP. 

These observations are particularly interesting because Table \ref{tab:error_varyn} shows that the errors do not increase significantly when $n$ is increased. The numbers between round brackets in Table \ref{tab:error_varyn} indicate the estimated standard deviation based on 20 realizations of $\Omega_i$. The norm $\| \cdot\|_{L^\infty}$ is defined as
\begin{equation}
\| \bm{x} \|_{L^\infty} := \max_t \sqrt{\Delta \theta (\bm{x}(t))^\top \bm{x}(t)}.
\end{equation}

\begin{table}
\caption{Running times for a varying number of spatial grid points $n$ \\ ($h = 1$, $T = 15$, $\tau = 10$)}
\centering
\begin{tabular}{|l||r|r|r|} \hline
Running times [s]  & \multicolumn{1}{|c|}{$n=10$} & \multicolumn{1}{|c|}{$n=100$} & \multicolumn{1}{|c|}{$n=1000$} \\ \hline \hline
Optimal Control    & 12.9 ($\pm$1.35) & 32.7 ($\pm$1.53) & 218.5 ($\pm$4.36) \\ \hline
MPC          & 4.6 ($\pm$0.44) & 11.0 ($\pm$0.66) & 70.5 ($\pm$4.51) \\ \hline
RBM-MPC      & 2.0 ($\pm$0.23)  &  3.6 ($\pm$0.71) & 14.1 ($\pm$1.33)  \\ \hline
\end{tabular}
\label{tab:times_varyn}
\end{table}

\begin{table}
\caption{Errors for a varying number of spatial grid points $n$ \\ ($h = 1$, $T = 15$, $\tau = 10$)}
\centering
\begin{tabular}{|l||c|c|c|} \hline
Relative errors [-]  & $n=10$ & $n=100$ & $n=1000$ \\ \hline \hline
$| u_{R-M} - u^*_\infty |_{L^2}$ & 0.76 ($\pm$0.28) & 0.59 ($\pm$0.23) & 0.53 ($\pm$0.11) \\ \hline
$| u_{R-M} - u_M |_{L^2}$ & 0.63 ($\pm$0.30) & 0.41 ($\pm$0.27) & 0.33 ($\pm$0.16) \\ \hline
$| u_M   - u^*_\infty |_{L^2}$ & 0.41 ($\pm$0.00) & 0.37 ($\pm$0.00) & 0.39 ($\pm$0.00) \\ \hline \hline
$\| \bm{x}_{R-M} - \bm{x}^*_\infty \|_{L^\infty}$ & 0.35 ($\pm$0.09) & 0.28 ($\pm$0.08) & 0.25 ($\pm$0.04) \\ \hline
$\| \bm{x}_{R-M} - \bm{x}_M \|_{L^\infty}$ & 0.17 ($\pm$0.08) & 0.11 ($\pm$0.07) & 0.08 ($\pm$0.04)  \\ \hline
$\| \bm{x}_M   - \bm{x}^*_\infty \|_{L^\infty}$ & 0.24 ($\pm$0.00) & 0.22 ($\pm$0.00) & 0.22 ($\pm$0.00) \\ \hline
\end{tabular}
\label{tab:error_varyn}
\end{table}

The convergence rates from Main Result \ref{res:main_2} and Corollary \ref{corr:RBMMPCconv} are validated in Figure \ref{fig:state_conv} and \ref{fig:control_conv}. Figures \ref{fig:state_conv_h} and \ref{fig:control_conv_h} show that $\|\bm{x}_{R-M}(\Omega_i) - \bm{x}_M \|_{L^\infty}$ and $|u_{R-M}(\Omega_i) - u_M|_{L^2}$ decay as $\sqrt{h}$ for $h \rightarrow 0$, which is in agreement with Main Result \ref{res:main_2}. Also note that $\bm{x}_{R-M}(\Omega_i)$ and $u_{R-M}(\Omega_i)$ do not converge to $\bm{x}^*_\infty$ and  $u^*_\infty$ for $h \rightarrow 0$, as Main Result \ref{res:main_2} and Corollary \ref{corr:RBMMPCconv} indicate. 
Figures \ref{fig:state_conv_T} and \ref{fig:control_conv_T} show that $\| \bm{x}_M - \bm{x}^*_\infty \|_{L^\infty}$ and $|u_M - u^*_\infty|_{L^2}$ are proportional to $e^{-2\mu_\infty T}$, which is agreement with Corollary \ref{corr:RBMMPCconv}. Increasing $T$ increases $\| \bm{x}_{R-M}(\Omega_i) - \bm{x}_M \|_{L^\infty}$ and $|u_{R-M}(\Omega_i) - u_M|_{L^2}$, which confirms that the constant $C_{[A,B,Q,R,T]}$ in Main Result \ref{res:main_2} increases with $T$. 
Figures \ref{fig:state_conv_tau} and \ref{fig:control_conv_tau} show that varying $\tau$ does not affect $\|\bm{x}_{R-M}(\Omega_i) - \bm{x}_M \|_{L^\infty}$ and $|u_{R-M}(\Omega_i) - u_M|_{L^2}$ strongly and  $\|\bm{x}_{M} - \bm{x}^*_\infty \|_{L^\infty}$ and $|u_M - u^*_\infty|_{L^2}$ increase with $\tau$. This is in agreement with Main Result \ref{res:main_2} and Corollary \ref{corr:RBMMPCconv}. 

\begin{figure*}
\centering
\subfloat[Varying $h$, $T = 15$, $\tau = 10$. \label{fig:state_conv_h}]{%
\includegraphics[width=0.3\textwidth]{Fig4a}}
\hfill
\subfloat[$h = 1$, varying $T$, $\tau = 10$. \label{fig:state_conv_T}]{%
\includegraphics[width=0.3\textwidth]{Fig4b}}
\hfill
\subfloat[$h = 1$, $T = 40$, varying $\tau$. \label{fig:state_conv_tau}]{%
\includegraphics[width=0.3\textwidth]{Fig4c}}
\caption{Differences between the RBM-MPC state trajectory $\bm{x}_{R-M}(\Omega_i, t)$, the MPC state trajectory $\bm{x}_M(t)$, and the infinite horizon state trajectory $\bm{x}^*_\infty(t)$ for $n =100$. The error bars indicate the $2\sigma$ confidence intervals estimated based on 20 realizations of $\Omega_i$. }
\label{fig:state_conv}
\end{figure*}

\begin{figure*}
\centering
\subfloat[Varying $h$, $T = 15$, $\tau = 10$. \label{fig:control_conv_h}]{%
\includegraphics[width=0.3\textwidth]{Fig5a}}
\hfill
\subfloat[$h = 1$, varying $T$, $\tau = 10$. \label{fig:control_conv_T}]{%
\includegraphics[width=0.3\textwidth]{Fig5b}}
\hfill
\subfloat[$h = 1$, $T = 40$, varying $\tau$. \label{fig:control_conv_tau}]{%
\includegraphics[width=0.3\textwidth]{Fig5c}}
\caption{Differences between the RBM-MPC control $u_{R-M}(\Omega_i,t)$, the MPC control $u_M(t)$, and the infinite horizon control $u^*_\infty(t)$ for $n =100$. The error bars indicate the $2\sigma$ confidence intervals estimated based on 20 realizations of $\Omega_i$. }
\label{fig:control_conv}
\end{figure*}



The code used to generated the results in this section can be found on \url{https://github.com/DCN-FAU-AvH}. 


\section{Conclusion and Perspectives} \label{sec:conclusions}
This paper considers a randomized MPC strategy called RBM-MPC to efficiently approximate the solution of an infinite horizon linear-quadratic OCP. In RBM-MPC, the finite-horizon OCPs that need to be solved in each iteration of the MPC-algorithm are simplified by replacing the system matrix $A$ by a randomized system matrix. RBM-MPC can reduce the computational cost of MPC significantly. 

There are four critical parameters that determine the stability and convergence of the proposed strategy: the prediction horizon $T$, the control horizon $\tau$, the time scale for the randomization $h$, and variance in the randomization of the system matrix $\mathrm{Var}[A_R]$. The estimates in this paper demonstrate that the RBM-MPC strategy is stabilizing when $T-\tau$ is sufficiently large and $h \mathrm{Var}[A_R]$ is sufficiently small and that the RBM-MPC strategy converges in expectation to an infinite-horizon OCP when $T-\tau \rightarrow \infty$ and $h \mathrm{Var}[A_R] \rightarrow 0$. 

The computational advantage of the RBM-MPC strategy was demonstrated for a heat equation with periodic boundary conditions discretized with $n = 100$ spatial grid points. In this example, RBM-MPC is 9 times faster than solving the OCP direcly and 3 times faster than classical MPC. The relative speed-up of the RBM-MPC strategy increases with the number of spatial grid points $n$. 

The analysis in this paper opens up several extensions and interesting possibilities for future work. 
\begin{itemize}
\item \textbf{Higher-order moments.} Main Result \ref{res:main_1} bounds the first-order moment $\mathbb{E}[|\bm{x}_{R-M}(t)|]$ based on an estimate for $\mathbb{E}[|J_R(\bm{u}_R^*) - J_T(\bm{u}^*_T)|]$ from \cite{rbm}. For $\mathbb{E}[|\bm{x}_{R-M}(t)|^2]$ an estimate for $\mathbb{E}[|J_R(\bm{u}_R^*) - J_T(\bm{u}^*_T)|^2]$ should be proved. 

\item \textbf{Terminal constraint.} RBM-MPC with a terminal constraint $\bm{x}_R(\bm{\omega}_i,\tau_i +T;\bm{x}_i,\tau_i) = 0$ in the OCP \eqref{eq:JT}--\eqref{eq:dyn_xT} can be analyzed in a future work based on the literature on MPC with a terminal constraint, see e.g., \cite{mayne2000}. 

\item \textbf{Constraints.} Control constraints can be directly included in the RBM-MPC algorithm. Including state constraints is more challenging because the plant trajectory $\bm{y}_R(t)$ does not necessarily satisfy the constraints when the predicted trajectory $\bm{x}_R^*(t)$ does. 

\item \textbf{Nonlinear settings.} It follows as in \cite{mpc} that RBM-MPC based on a linear plant model is locally stabilizing for a nonlinear plant when the plant model is a sufficiently good approximation of the nonlinear plant near the origin. The RBM-MPC strategy has also been applied with a nonlinear plant model in \cite{ko2021model} but the analysis for this problem setting is still open. 

\item \textbf{Machine Learning.} Because the training of residual Deep Neural Networks (DNNs) can be seen viewed as a nonlinear OCP \cite{benning2019, esteve2020}, RBM-MPC may also be applied to speed up the training of DNNs. 

\item \textbf{Gradient descent for RBM-constrained OCPs.} A natural way to solve RBM-contrained OCPs is to use a different realization of the randomized $A$-matrix in each forward and backward pass, see also Remark \ref{rem:RBMGD}. The rigorous analysis of this algorithm is still missing. 

\item \textbf{First-RBM-then-MPC.} It is natural to also analyze the first-RBM-then-MPC algorithm and compare it to the first-MPC-then-RBM algorithm considered in this paper, see Remark \ref{rem:MPCRBM}. The stochastic dependencies in first-RBM-then-MPC will be somewhat more involved. 


\item \textbf{PDEs on Networks.} The numerical example in Section \ref{sec:numerical} demonstrates that removing cycles by the RBM can speed up computations significantly. These ideas can also be applied to the complex networks of (1-D) PDEs that e.g.\ appear in the modeling of gas transport \cite{herty2007, leugering2017}.

\item \textbf{Adaptive RBM-MPC} The computational cost for RBM-MPC may be reduced further by adapting the parameters $h$, $\tau$, and $T$ over the iterations. Developing such adaptation strategies for MPC and RBM-MPC and proving their stability is an interesting topic for future research. 

\end{itemize}



\bibliographystyle{IEEEtran}
