\section{Introduction}

Multi-agent systems with non-interacting agents have a wide range of potential applications \cite{marine_search, agbots_Chowdhary, agriculture}. In such systems, individual agents can rely only on the commands issued by the central node because they do not have access to the information about their neighbors. This framework is advantageous when individual agents do not have enough computational power and the central node has exact model of each agent. Since modeling errors and environmental uncertainties are unavoidable in practice, purely centralized and model based controllers exhibit some unwanted behaviors.
\par Learning based methods became extremely successful in approximating the complex features, which is otherwise very difficult. Since these methods can exhibit some unsafe behavior during the learning phase, their implementation on \emph{safety critical systems} is prohibited \cite{survey_LMPC, robust_action_governor, safe_rl_Krause, robust_regression, safe_exploration, safe_exploration_Krause, safe_RL_Tomlin}. Some newly developed approaches like \emph{system level synthesis} \cite{SLS_safe}, stochastic learning \emph{model predictive control} (MPC) \cite{learning_MPC_Mesbah}, reachable tube learning \cite{deepreach} and iterative learning MPC for autonomous car racing  \cite{learning_autonomous_car, LMPC_Borrelli}, employ MPC to ensure safety.
\par The constraint satisfaction capability of MPC is combined with function approximation capability of \emph{deep neural networks} (DNN) in several recent works; please see \cite{Deep_MPC_Gopaluni, Deep_MPC_Lucia_18, rapid_MPC_Borreli, minimax_Deep_MPC, Allgower_MPC_approximate} and references therein. A linear tube based MPC is employed along with a \emph{neural network} (NN) in \cite{NNMPC_Findeisen} in which MPC is augmented by NN to deal with nonlinearities arising due to the linearization around the operating point. In order to satisfy constraints, first MPC is designed, then NN is constrained to the remaining action space. Alternatively, a robust action governor is utilized in \cite{robust_action_governor}, which projects the learning based control to the safe action space.
\par When (local) uncertainties are structured, they can be simply represented in terms of (possibly) high dimensional feature basis functions. Local agents, equipped with some learning mechanism, act on their respective disturbances by themselves and transmit only low dimensional estimates of their actions to the central node instead of high dimensional features. The computationally powerful central node ensures the satisfaction of underlying constraints by using only the received estimates. Recently, this framework is utilized in \cite{MWGC} with a central MPC equipped with local adaptation for structured uncertainties. In this current article, we extend the results of \cite{MWGC} by considering the unstructured uncertainties.
\par Generally, when uncertainties are unstructured, neuro-adaptive controllers are utilized in the framework of \emph{model reference adaptive control} (MRAC) \cite{MRAC_Nguyen}. One such method is \cite{HHN_08} in which the inner layer weights of DNN are randomly assigned and only the outer layer is trained by using the MRAC type weight update law. Although adaptation of only the outer layer weights gives some acceptable performance, it cannot learn complex features \cite{Chowdhary_GP_14, GPMRAC_18, DMRAC}. In \cite{DMRAC}, authors propose DNN based adaptive architecture for continuous time linear systems in which the outer layer weights are updated in real time while the inner layer weights are trained on a slower time scale using batch updates.

\par In this article, we follow the method proposed in \cite{DMRAC} to train the inner layer weights for a control affine non-linear discrete time system. The outer layer weights are updated according to discrete adaptive weight update law. We employ tube based MPC to ensure that the uncertain state trajectory stay within a tube around a nominal reference trajectory.  In contrast to \cite{NNMPC_Findeisen, robust_action_governor}, we employ bounded activation functions in the outermost layer and project their weights in a bounded set to ensure boundedness of the learning based controller. Then we assign the remaining control authority to MPC for the sake of optimality under constraints.
\par This manuscript is organized as follows: We present a problem setup in \secref{s:setup}. The formulation of adaptive controller and the MPC controller are given in \secref{s:adaptive} and \secref{s:MPC}, respectively. Our overall algorithm and its stability is discussed in \secref{s:stability}. We validate our theoretical results with the help of a numerical experiments in \secref{s:experiment} and conclude in \secref{s:epilogue}.

	\section{Problem setup}\label{s:setup}
	We let $\mathds{R}$ denote the set of real numbers, $\mathds{N}_0$ the set of non-negative integers and $\mathds{Z}_{+}$ the set of positive integers. For a given vector $v$ and positive (semi)-definite matrix $M \succeq \mathbf{0}$, $\norm{v}_M^2$ is used to denote $v \transp M v$. For a given matrix $A$, the trace, the largest eigenvalue, pseudo-inverse and Frobenius norm are denoted by $\trace(A)$, $\lambda_{\max}(A)$, $A^{\dagger}$ and $\norm{A}_F$, respectively. By notation $\norm{A}$ and $\norm{A}_{\infty}$, we mean the standard $2-$norm and $\infty-$norm, respectively, when $A$ is a vector, and induced $2-$norm and $\infty-$norm, respectively, when $A$ is a matrix. A vector or a matrix with all entries $0$ is represented by $\mathbf{0}$ and $I$ is an identity matrix of appropriate dimensions.
	\par Let us consider a discrete time dynamical system
	\begin{equation}\label{e:system}
	x_{t+1} = f(x_t) + g(x_t)\left( u_t + h(x_t) \right),
	\end{equation}
	where
	\begin{enumerate}[leftmargin = *, nosep, label=(1-\alph*), widest = b]
		\item \label{e:constraints} $x_t \in \mathcal{X} \subset \mathds{R}^d$, $ u_t \in \mathds{U} \coloneqq \{v \in \mathds{R}^m \mid \norm{v}_{\infty} \leq u_{\max} \}$,
		\item $\mathcal{X}\subset \mathds{R}^d$ is a compact set,
		\item system function $f:\mathds{R}^d \rightarrow \mathds{R}^d $, control influence function $g:\mathds{R}^d \rightarrow \mathds{R}^m$ are given Lipschitz continuous functions,
		\item $h(x_t)$ is the state dependent matched uncertainty at time $t$ such that $g(x_t)h(x_t) \in \mathds{W} \coloneqq \{v \in \mathds{R}^d \mid \norm{v}_{\infty} \leq w_{\max} \}$, $h$ is continuous and $\norm{g(x_t)} \leq \delta_g$ for some $\delta_g> 0$ for every $x_t \in \mathcal{X}$.
	\end{enumerate}

We assume that the multi-agent system consists of a family of finite number of agents with each agent having the dynamics in the form \eqref{e:system} but $f,g,h, \mathcal{X}$, and $\mathds{U}$ may be different for each agent. One crucial assumption is that agents are uncoupled and non-interacting. The extension of the present approach for interacting agents with coupled dynamics needs separate treatment. So we do not consider that in this article. Some recent results on distributed MPC \cite{Lygeros_distributedMPC_19} may be helpful in this context. Since the agents are non-interacting, we consider the dynamics of only one agent for notational simplicity and without loss of generality \cite[Remark 1]{MWGC}. The problem description of this article is as follows:
	\begin{pstatement}\rm{
			Present a stabilizing control framework for \eqref{e:system}, which can be generalized to non-interacting, multi-agent systems, respects physical constraints \ref{e:constraints}, optimizes a given performance index, and reduces the effect of unstructured uncertainties by using neural networks for each agent.}
	\end{pstatement}

	\par Our proposed solution is based on constraint satisfaction and cost minimization capabilities of MPC, and universal approximation property of neural networks.
	We break the appied control $u_t$ such that
	\begin{equation}\label{e:total_conrol}
	u_t = u_t^a + u_t^m,
	\end{equation}
	where $u_t^a$ and $u_t^m$ are the neuro-adaptive and MPC components, respectively. The centralized MPC controller employs only the nominal dynamics of \eqref{e:system}, which is given below for easy reference
	\begin{equation}\label{e:nominal}
	x_{t+1} = f(x_t) + g(x_t) u_t^m  \coloneqq \bar{f}(x_t, u_t^m).
	\end{equation}
	Therefore, the dynamics \eqref{e:system} can be written as
	\begin{equation}\label{e:actual_agent}
	x_{t+1} = \bar{f}(x_t, u_t^m) + g(x_t) \left( u_t^a + h(x_t)\right).
	\end{equation}

	\par Please notice that MPC experiences the term $g(x_t) \left( u_t^a + h(x_t)\right)$ as a disturbance. In a broader sense, the MPC component $u_t^m$ is responsible for input-to-state stability (ISS) of closed-loop states in the presence of bounded disturbances, and the neuro-adaptive control component $u_t^a$ acts on $h(x_t)$. In particular, the job of $u_t^a$ is to approximate $-h(x_t)$ and stay within a known bound during the learning phase so that MPC can always experience a bounded disturbance with a known bound.
	\par In the following sections \secref{s:adaptive} and \secref{s:MPC}, we design the adaptive controller $u_t^a$ and MPC controller $u_t^m$, respectively.
	\section{Learning based agent}\label{s:adaptive}
	 Any continuous function $h$ on a compact set $\mathcal{X}$ can be approximated by a multi-layer network with number of layers $L \geq 2$ such that
	\begin{equation}
	h(x) = W_L\transp \psi_L \left[ W_{L-1}\transp \psi_{L-1} \left[ \cdots \left [ \psi_1(x)\right] \right] \right] + \varepsilon_1(x),
	\end{equation}
	where $x \in \mathcal{X}$, $\psi_i, W_i$ for $i=1, \ldots , L$, are activation functions and ideal weights, respectively, in the $i$th layer. The reconstruction error function $\varepsilon_1$ is bounded by a known constant $\bar{\varepsilon}_1 > 0$ for each $x \in \mathcal{X}$, i.\ e.\ $\norm{\varepsilon_1(x)} \leq \bar{\varepsilon}_1$. Therefore, we can represent $h(x_t)$ with the help of a neural network with a desired accuracy. If the neural network is not minimal then the ideal weights may not be unique. However, for the neural-adaptive controller design only the existence of ideal weights is assumed, which is always guaranteed when $h$ is a continuous function on a compact set \cite[\S 7.1]{Lewis_NN_99}. Let us define $\phi^{\ast}(x) \coloneqq \psi_L \left[ W_{L-1}\transp \psi_{L-1} \left[ \cdots \left[ \psi_1(x)\right] \right] \right]$ and $W^{\ast} \coloneqq W_{L-1} $, then
	\begin{equation}\label{e:NN_structure}
	h(x_t) = W^{\ast \top}\phi^{\ast}(x_t) + \varepsilon_1(x_t),
	\end{equation}
	where $W^{\ast} \in \mathds{R}^{(n_L + 1) \times m}$ and there are $n_L$ number of neurons in the last layer. The ideal inner layer weights defining $\phi^{\ast}(\cdot)$ are neither known nor unique. But for a given family $\mathcal{F}$ such that $\phi^{\ast} \in \mathcal{F}$, we can assume that there exist $\varepsilon_{\phi}: \mathds{R}^d \rightarrow \mathds{R}^m$ for each $\phi \in \mathcal{F}$ such that $\phi^{\ast}(x) = \phi(x) + \varepsilon_{\phi}(x)$ for each $x \in \mathcal{X}$. Therefore, \eqref{e:NN_structure} becomes
	\begin{equation}\label{e:NN_structure_general}
	h(x_t) = W^{\ast \top}\phi(x_t) + \varepsilon(x_t),
	\end{equation}
	where $\varepsilon(x_t) = \varepsilon_1(x_t) + W^{\ast \top}\varepsilon_{\phi}(x_t)$ is the overall reconstruction error.  Please notice that even when the hidden layer weights are randomly assigned as in ELM \cite{ELM_06}, the universal approximation property of the neural network allows us to make the overall reconstruction error $\varepsilon(\cdot)$ as small as desired by increasing the width of the network. A network with trained hidden layers can capture several useful features, which in turn results in performance improvement \cite{DMRAC}. We train the outermost layer in real time with the help of weight update law and the inner layers by using the approach \cite{DMRAC}. In the next subsections we provide the relevant details of the training of DNN.
	\subsection{Training of the outermost layer}\label{s:outer layer training}
	\par If the activation function in the last layer is bounded, i.\ e.\ sigmoidal, tanh, etc., then we can assume that $\norm{\phi(x)}$ is bounded \cite{HHN_08, NN_15}. Further, if we have a priori knowledge about the bounds on the ideal weights of the last layer then we can use this knowledge to avoid parameter drift phenomenon. We make the following assumption:
	\begin{assumption}\label{as:bounds_uncertainty}
		\rm{
			There exist $\bar{W}_{i} > 0$ for $i = 1, \ldots , m$, and $\sigma, \bar{\varepsilon} > 0$ such that $\norm{W^{\ast (i)}} \leq \bar{W}_{i}$, for $i = 1, \ldots , m$, and $\norm{\phi(x)} \leq \sigma, \sup_{\phi \in \mathcal{F}} \norm{\varepsilon(x)} \leq \bar{\varepsilon}$ for every $x \in \mathcal{X}$. }
	\end{assumption}
	The above assumption is standard in literature \cite{HHN_08, NN_15, Joshi_DMRAC}. When the feature function $\phi$ varies with time due to the training of inner layer weights, a more rigorous analysis may be performed along the lines of \cite{time-varying_features}.
	We employ \[ u_t^a = -K_t\transp \phi(x_t) \] as an adaptive control at time $t$, where $K_t$ is the weight of the outer layer which is trained according to the adaptive weight update law and $\phi \in \mathcal{F}$ is a feature basis function obtained from the inner layers of deep neural network. Let us initialize $K_0$ such that $\norm{K_0^{(i)}} \leq \bar{W}_i$; $i = 1, \ldots , m$. For a given learning rate $0< \theta < \frac{1}{\sigma^2}$, we employ the following weight update law:
	\begin{equation}\label{e:update_law}
	\bar{K}_{t+1} = K_t - \theta \phi(x_t)(x_{t+1} - \bar{f}(x_t, u_t^m))\transp g(x_t)^{\dagger},
	\end{equation}
	where $g(x_t)^{\dagger}$ represents pseudo-inverse of $g(x_t)$. When $g(x_t)^{\dagger}$ does not exist, we use $\bar{K}_{t+1} = K_t$. We employ the discrete projection method to ensure boundedness of $K_t^{(i)}$ for $i= 1, \ldots , m$, as follows:
	\begin{equation}\label{e:projection}
	K_t^{(i)} = \proj \bar{K}_t^{(i)} = \begin{cases} \bar{K}_t^{(i)} & \text{ if } \norm{\bar{K}_t^{(i)}} \leq \bar{W}_{i} \\
	\frac{\bar{W}_{i}}{\norm{\bar{K}_t^{(i)}}} \bar{K}_t^{(i)} & \text{ otherwise. } \end{cases}
	\end{equation}
	One important property of the projection is the following \cite[\S 4.11.3]{Ioannou_Fidan_tutorial}, \cite[page 341]{Landau_Karimi_11}:
	\begin{equation}\label{e:effect_projection}
	(K_t^{(i)} - W^{\ast (i)})\transp (\bar{K}_{t}^{(i)} - K_{t}^{(i)})  \geq 0 \text{ for each } i = 1, \ldots , m.
	\end{equation}
	Let us recall the following definition:
	\begin{definition}[\cite{Ioannou_Fidan_tutorial}, page 117]\label{def:small_disturbance}
		\rm{
			The sequence vector $(s_t)_{t \in \mathds{N}_0}$ is called $\mu$ small in mean square sense if it satisfies $\sum_{t=k}^{k+N-1}\norm{s_t}^2 \leq Nc_0 \mu + c_0^{\prime}$ for all $k \in \mathds{Z}_{+}$, a given constant $\mu \geq 0$ and some $N \in \mathds{Z}_{+}$, where $c_0, c_0^{\prime} \geq 0$. }
	\end{definition}

Let $\tilde{K}_t \coloneqq K_t - W^{\ast}$ and $\tilde{u}_t \coloneqq u_t^a + h(x_t) = - \tilde{K}_t \transp \phi(x_t) + \varepsilon(x_t) $. Some straightforward arguments as in \cite[\S 4.11.3]{Ioannou_Fidan_tutorial} give us the following result:
\begin{lemma}\label{lem:adaptation}
	Consider the dynamical system \eqref{e:system}, weight update law \eqref{e:update_law} and the projection method \eqref{e:projection}. Let the Assumption \ref{as:bounds_uncertainty} hold. Then
$\tilde{u}_t$ is $\bar{\varepsilon}^2$ small in mean square sense with $c_0 = \frac{1}{1-\theta \sigma^2}$ and $c_0^{\prime} = \frac{4c_0}{\theta}\bar{W}$.
\end{lemma}
The Proof of Lemma \ref{lem:adaptation} is standard in literature. We provide it in the appendix for completeness.

Let $\bar{W} \coloneqq \sum_{i=1}^m \bar{W}_i^2$. It is evident that $\trace(K_t\transp K_t) \leq \bar{W}$ for all $t$ due to the projection.
Therefore, the neuro-adaptive control component $u_t^a$ is bounded, i.\ e.\
\begin{equation}
\norm{u_t^a}_{\infty} \leq \norm{u_t^a} \leq \sqrt{\bar{W}}\sigma \eqqcolon u_{\max}^a.
\end{equation}
The apparent disturbance term $g(x_t) \left( u_t^a + h(x_t)\right)$ in \eqref{e:actual_agent} is also bounded, i.\ e.\
\begin{equation}\label{e:overall_disturbance_bound}
\begin{aligned}
& \norm{g(x_t) \left( u_t^a + h(x_t)\right)}_{\infty} \leq \norm{g(x_t){u}_t^a} + \norm{g(x_t)h(x_t)} \\
& \quad \leq \delta_g\norm{u_t^a} + w_{\max} \leq \delta_g u_{\max}^a + w_{\max} \eqqcolon w_{\max}^{\prime}.
\end{aligned}
\end{equation}
\subsection{Training of the inner layers}\label{s:inner layer training}
\par Since the training of DNN is computationally demanding, we train inner layers separately with a slower rate than that of the outer layer. The inner layer weights are used to get feature basis function $\phi$ for the real time adaptation of the outer layer weights in \secref{s:outer layer training} and collection of data. The qualified data is used to train the inner layers, while $\phi$ in \secref{s:outer layer training} remains fixed until the training is complete.
\par The inner layers are trained by using the stochastic gradient descent method. In order to train them, independent and identically distributed samples of the ordered pair $\left( x, h(x) \right)$ are needed. Since the direct access of $h(x)$ for a given $x$ is not possible, we use a generative network to generate their estimates $\hat{h}(x)$ for each given $x$.
\par The generative network maintains a replay buffer in which $p_{max}$ number of pairs $(x, \hat{h}(x))$ are stored. In practice, the adaptive control component $u_t^a$ is used as a guess to estimate $-h(x_t)$ at time $t$. For $t > p_{\max}$, new pairs are added in the buffer only if they satisfy certain singular value maximization criterion \cite{svd_cl} after removing some less significant data from the buffer. Depending upon the computational power of the generative network, $p_0 < p_{\max}$ members of the buffer are randomly sampled to train the inner layers of the network. Once the training is complete, the feature basis function $\phi$ in \secref{s:outer layer training} is replaced by this new feature basis function.
\par In the next section, we provide a summary and important results of tube-based MPC relevant to the context of the present article.

	\section{Centralized controller}\label{s:MPC}
	 The discussion in this section closely follows \cite{Mayne_tube_NLMPC} with differences related to the adaptive controller discussed in \secref{s:adaptive}. Tube based MPC is based on ensuring that the closed-loop states stay within a tube around a reference trajectory. The reference trajectory is obtained by solving a reference governor problem \emph{offline} under the tightened constraints for regulation problems. Once a trackable reference trajectory is obtained by spending only a part of the available control authority, a reference tracking problem without state constraints is solved \emph{online} that utilizes full control authority.
	 \par Constraint tightening in the reference governor allows satisfaction of actual constraints by the actual states.  Knowledge of the exact bound on disturbance is needed to tighten the constraints. Although the disturbance in dynamical system \eqref{e:system} at time $t$ is $g(x_t)h(x_t)$, the disturbance experienced by MPC is $g(x_t)(u_t^a + h(x_t))$, which is uniformly bounded due to \eqref{e:overall_disturbance_bound}. Therefore, we re-define the disturbance set for MPC as follows:
	\[ \mathds{W}^{\prime} \coloneqq \{v \in \mathds{R}^d \mid \norm{v} \leq w_{\max}^{\prime} \} .\]
	We also notice that a part of control is utilized by the adaptive controller. Therefore, we need to re-define the control set as well.
	\[ \mathds{U}^{\prime} \coloneqq \{v \in \mathds{R}^m \mid \norm{v}_{\infty} \leq u_{\max} - u_{\max}^a \}. \]
	These modifications in tube-based MPC are already pointed out in \cite{MWGC, NNMPC_Findeisen} for similar context and in \cite{Bhasin_AMPC_19} for the linear systems with parametric uncertainties. We present a quick summary in this section for completeness. For a detailed discussion please refer to \cite{Mayne_tube_NLMPC, MWGC}.

	\par For some optimization horizon $N \in \mathds{Z}_{+}$, an offline reference governor is utilized to generate the reference trajectory
	 \begin{equation}\label{e:reference_signal}
	 \begin{aligned}
	 (x_t^r)_{t \in \mathds{N}_0} &\coloneqq \{ (x_t^r)_{t = 0}^{N-1}, 0, \ldots \}, \\
	 (u_t^r)_{t \in \mathds{N}_0} &\coloneqq \{ (u_t^r)_{t = 0}^{N-1}, 0, \ldots \} .
	 \end{aligned}
	 \end{equation}
	 such that $x_0^r = x_0$ and $x_i^r \in \mathcal{X}_r \subset \mathcal{X}, u_i^r \in \mathds{U}_r \subset \mathds{U}^{\prime}$, where $\mathcal{X}_r$ and $\mathds{U}_r$ are tightened state set and control set, respectively. The tightened constraint sets may be obtained by following the approach of \cite[Section 7]{Mayne_tube_NLMPC}.
	 \par In order to design the online reference tracking MPC, we first choose an optimization horizon $N \in \mathds{Z}_{+}$, which can be different from that chosen for the reference governor. Let $Q, R, Q_f \succ \mathbf{0}$ be given positive definite matrices of appropriate dimensions. Let
	 \[ c_{\mathrm{s}}(x_{t+i \mid t}, u_{t+i \mid t}) \coloneqq \norm{x_{t+i \mid t} - x_{t+i}^r}_Q^2 + \norm{u_{t+i \mid t} - u_{t+i}^r}_R^2\]
	 be the cost per stage at time $t+i$ predicted at time $t$ and let $c_{\mathrm{f}}(x) = x\transp Q_f x$ be the terminal cost. Let us define the terminal set $\mathcal{X}_f$ as in \cite{Mayne_tube_NLMPC} and make the following assumption:
	 \begin{assumption}\label{as:stability}
	 	\rm{
	 		There exists a control $u^{\prime} \in \mathds{U}^{\prime} $ such that the following holds
	 		\begin{equation}
	 		c_{\mathrm{f}} \left( \bar{f}(x, u^{\prime}) \right) + c_{\mathrm{s}}(x, u^{\prime}) -c_{\mathrm{f}}(x)  \leq 0
	 		\end{equation}
	 		for every $x \in \mathcal{X}_f$.
	 	}
	 \end{assumption}

 The online reference tracking MPC solves the following optimal control problem at each time instant:
 \begin{equation}\label{e:MPC}
 \begin{aligned}
V_m(x_t) \coloneqq \min_{(u_{t+i \mid t})_{i=0}^{N-1}} & c_{\mathrm{f}}(x_{t+N \mid t}) + \sum_{i=0}^{N-1}c_{\mathrm{s}}(x_{t+i \mid t}, u_{t+i \mid t})  \\
 \sbjto \quad & x_{t \mid t} = x_t \\
 & \text{for } i = 0, \cdots, N-1  \\
 & x_{t+ i +1 \mid t} = \bar{f}(x_{t+i \mid t}, u_{t+i\mid t})   \\
 & u_{t+i \mid t} + u_t^a \in \mathds{U}.
 \end{aligned}
 \end{equation}
 The first control $u_{t\mid t}$ is called the MPC component $u_t^m$ and is applied along with $u_t^a$ to the system at time $t$, which in turn results in the satisfaction of control constraints. Please notice that the \emph{online reference tracking} problem does not have state constraints and it does not take into account the modified disturbance and control sets, $\mathds{W}^{\prime}, \mathds{U}^\prime$, respectively.
 \par Let $X_c(x_t^r)$ be the level set around $x_t^r$ of radius $c$ generated by $V_m(x_t)$ and $X_c$ be their union. The choice of the terminal cost functional $c_{\mathrm{f}}(\cdot)$ and the construction of terminal set ensure that $x_t \in X_c(x_t^r)+ \mathds{W}^{\prime}$ for every $t$, which in turn ensures the satisfaction of state constraints of the uncertain system. We recall the following results:
\begin{lemma}{\cite[Proposition 3]{Mayne_tube_NLMPC}}\label{lem:lipschitz_value}
		There exists $c_3>0$ such that $x_t \mapsto V_m(x_t)$ is Lipschitz continuous with Lipschitz constant $c_3$, when $x_t \in X_c(x_t^r)+ \mathds{W}^{\prime}$ for every $t$.
\end{lemma}
	\begin{lemma}{\cite[Proposition 2]{Mayne_tube_NLMPC}}\label{lem:bound_on_predicted_next}
		If Assumption \ref{as:stability} is satisfied, then for every $x_t \in X_c+\mathds{W}^{\prime}$ the following holds
		\begin{equation}\label{e:bound_on_predicted_next}
		V_m(x_{t+1 \mid t}) - V_m(x_t) \leq -c_{\mathrm{s}}(x_t, u_{t \mid t}) \text{ for all } t.
		\end{equation}
	\end{lemma}

\section{Algorithm and stability}\label{s:stability}
	In this section, we present our algorithm and provide our main result on stability. For a regulation problem, the reference governor is utilized at the beginning to generate a trackable reference trajectory. In case of reference tracking problem, the same reference governor can be invoked whenever the target state changes. At the beginning, the inner layer weights are randomly assigned and the outer layer weights are set to zero. A finite size replay buffer is used to store the state and adaptive action pair $(x_t, u_t^a)$, where $u_t^a$ is computed according to \secref{s:outer layer training}. The feature basis function $\phi$ needed to computed $u_t^a$ is time varying and represented by the inner layers of DNN. Once the sufficient data is stored in the replay buffer, a random subset of it is used to train the DNN with the help of generative network as described in \secref{s:inner layer training}. The replay buffer is updated by SVD maximization approach \cite{svd_cl} and DNN is trained after each fixed interval of time. Once the training of DNN is complete, $\phi$ is also updated. The adaptive control component $u_t^a$ is communicated with MPC, which solves \eqref{e:MPC} to compute $u_t^m$. The control $u_t = u_t^a + u_t^m$ is applied to the system. The same procedure is repeated after measuring the next state.
	Our algorithm is summarized in Algorithm \ref{a:algo}.
	\begin{algorithm}
		\caption{Deep MPC}
		\label{a:algo}
		\begin{algorithmic}[1]
			\Require  $x_0$
			\State choose $\theta < \frac{1}{\sigma^2}$
			\State Generate the reference trajectory \eqref{e:reference_signal} under the tightened constraints
			\State initialize $K_0 = \mathbf{0}$, $t=0$ and randomly assign weights of inner layers to get $\phi$
			\State initialize a generative network on a parallel machine
			\For{each $t$}
			\State update $\phi$ if training in the generative network is complete, otherwise move to the next step
			\State \label{al:first_step}compute $u_t^a = - K_t \transp \phi(x_t)$
			\State solve \eqref{e:MPC}, set $u_t^m = u_{t \mid t}$
			\State apply $u_t = u_t^m + u_t^a$ to the system and measure $x_{t+1}$
			\State compute $K_{t+1}$ by the weight update law \eqref{e:update_law} and \eqref{e:projection}
			\EndFor
		\end{algorithmic}
	\end{algorithm}

	\par We have the following main result:
	\begin{theorem}\label{th:stability}
		Consider the dynamical system \eqref{e:system} and let assumptions \ref{as:bounds_uncertainty} and \ref{as:stability} hold. Then under the control computed by the Algorithm \ref{a:algo} the closed-loop system \eqref{e:system} is ISS.
	\end{theorem}
	\begin{proof}
		Since $Q \succ 0$ and $f,g$ are Lipschitz continuous, by \cite[Lemma 3]{MWGC} there exist $c_1, c_2 > 0$ such that
		\begin{equation}
		c_1 \norm{x_t -x_t^r}^2 \leq V_m(x_t) \leq c_2 \norm{x_t -x_t^r}^2.
		\end{equation}
		We compute a bound on $\Delta_m(x_t) = V_m(x_{t+1})- V_m(x_t)$ as follows:
		\begin{align}\label{e:ISS}

		\Delta_m(x_t) & = V_m(x_{t+1 \mid t})- V_m(x_t) + V_m(x_{t+1}) - V_m(x_{t+1 \mid t}) \notag \\
		& \leq - c_{\mathrm{s}}(x_t, u_t^m) + V_m(x_{t+1}) - V_m(x_{t+1 \mid t}) \text{ by Lemma \ref{lem:bound_on_predicted_next}} \notag \\
		& \leq - c_{\mathrm{s}}(x_t, u_t^m) + c_3\norm{x_{t+1}-x_{t+1\mid t}} \quad \text{ by Lemma \ref{lem:lipschitz_value}} \notag \\
		& = - c_{\mathrm{s}}(x_t, u_t^m) + c_3\norm{g(x_t) \left( u_t^a + h(x_t)\right)}.
		\end{align}
		Since $\norm{g(x_t) \left( u_t^a + h(x_t)\right)}\leq w_{\max}^{\prime}$ is bounded for all $t$, the dynamical system \eqref{e:system} is ISS under the control computed by the Algorithm \ref{a:algo}.
	\end{proof}
Please notice that  tube based MPC is ISS by its design. The adaptive component in the presence of structured uncertainty can ensure convergence of states to origin \cite{MWGC}, which is difficult when uncertainties are unstructured. However, since the adaptive component learns a representable part of the uncertainty, we can expect an improvement in performance, without loosing the guarantee ensured by MPC.

	\section{Numerical experiment}\label{s:experiment}
	We consider two fixed-wing unmanned air vehicles (UAV) with limited processing power. The UAVs are subject to adverse lateral-longitudinal wing-rock dynamics, which are experienced at high angle of attack. These UAVs cannot interact with each other but they can communicate with a central base station.
	\par Wing-rock is a benchmark dynamical system, which has been widely used to evaluate adaptive controllers. Unlike existing adaptive control results \cite{yucelen2012command,chowdhary2010concurrent,chowdhary2013concurrent}, we require that the control and system states remain bounded in a pre-specified constraint sets. In this example, we assume that both agents suffer from bounded and unstructured uncertainties.
	\par Letting $\delta_t$ denote the roll angle in radian, and $p_t$ denote the roll rate in radian per second, the state of the wing-rock dynamics model is $x_t \coloneqq \bmat{\delta_t & p_t}\transp$ at time $t$. Let $x_t^{(i)}$ denote the state of the $i$th agent for $i= 1,2$, we consider the following discrete time dynamics:
	\begin{equation}
	x_{t+1}^{(i)} = Ax_t^{(i)} + B\Big(u_t^{(i)} + h(x_t^{(i)})\Big),
	\label{eq:disc-sys}
	\end{equation}
	where $A=\begin{bmatrix}
	1 & 0.05 \\ 0 & 1
	\end{bmatrix}$, $B=\begin{bmatrix}
	0 \\ 0.05
	\end{bmatrix}$,
	and $h(\cdot)$ is bounded uncertainty. In order to generate $h$ for the purpose of simulation, we use $ h(x_t^{(i)}) = V_t\transp \varsigma(x_t^{(i)}) + \omega$, with
	 $V_t = v_t  V_0$,
	\[ V_0\transp = \bmat{0.8 & 0.2314 & 0.6918 & -0.6245 & 0.0095 & 0.0214}, \] and  $\omega \in [-\omega_0, \omega_0]$ is a truncated normal random variable with $\omega_0 = 0.457$. The time varying multiplier $v_t$ is given by
	\begin{equation*}
	   v_t =  \begin{cases}
	        4 \text{ for } t \in [0,49] \cup [100, 149] \\
	        0  \text{ otherwise}.
	    \end{cases}
	\end{equation*}
	The function $\varsigma(\cdot)$ is saturated by a standard saturation function as $\varsigma(x) = \sat(\varsigma^{\prime}(x))$, where $\varsigma^{\prime}(x) = \begin{bmatrix}
	1 & \delta & p & \abs{\delta} p & \abs{p} p &\delta^3
	\end{bmatrix}\transp $ and $\sat(\cdot)$ is a standard saturation function with the threshold $\frac{\omega_0}{10}$.

	The controller is not aware of $\varsigma(\cdot)$ and $\omega$.
	\par The goal is to stabilize the wing-rock dynamics by driving states to the origin from arbitrarily chosen initial states $x_0^{(1)}=\bmat{\frac{\pi}{30} & \frac{5\pi}{60}}\transp$ and $x_0^{(2)}=\bmat{\frac{\pi}{30} & -\frac{\pi}{18}}\transp$, respectively. For both agents, constraint sets are given below:
	\[ \mathcal{X}  = \left[-\frac{\pi}{6},\frac{\pi}{6} \right] \times \left[ -\frac{\pi}{2},\frac{\pi}{2} \right], \text{ and }
	\mathds{U}  = \left[-\frac{\pi}{2},\frac{\pi}{2} \right]. \]
	\par We compare our proposed approach (Deep MPC computed according to the Algorithm \ref{a:algo}) with two controllers, namely tube MPC and shallow MPC. Tube MPC is designed according to \cite{Mayne_tube_NLMPC}. In order to design shallow MPC, we follow our algorithm \ref{a:algo} but we consider only a single layer neural network with $3$ neurons. Therefore, shallow MPC utilizes our adaptive weight update law \eqref{e:update_law} to update their weights but the training of NN is not required. To design the Deep MPC, we use a four layers network with sizes $[2, 5, 5, 3]$ respectively, where the first hidden layer has 5 neurons and the outermost layer has $3$ neurons. The outermost layer is updated with our adaptive weight update law \eqref{e:update_law}, while the remaining three inner layers are trained using SGD with momentum constant 0.9 and learning rate 0.1. We use nonlinear activation functions after each of the inner layers, and these functions are respectively $[ReLU, ReLU, tanh]$. We use the MATLAB based software package MPT 3.0 \cite{MPT3} to compute the polytopes needed for the constraint tightening. The simulation parameters are chosen to be same for both agents with optimization horizon $N=20$, $\theta = 0.9, Q = I$, and $R=0.1$. For the reference governor, we used a larger optimization horizon such that the terminal constraint on state can be satisfied.
	\par Fig.\  \ref{fig:roll_angle} and Fig.\ \ref{fig:roll_rate} depict the trajectories of the states of the system. Fig.\ ~\ref{fig:roll_angle} demonstrates that the roll angle trajectory obtained by tube MPC has some offset and oscillatory behavior for both agents. The roll rate trajectories as shown in Fig.\ \ref{fig:roll_rate} are noisy in which tube MPC has the highest fluctuations. It is evident that shallow MPC has reasonably good performance at the beginning in both states for both agents and keeps doing well when multiplier $v_t$ changes to $0$. At $t=100$, when multiplier changes again to $4$, shallow MPC is not able to adapt this change. However, deep MPC outperforms both of them and keeps improving its performance even after the sudden change in the multiplier $v_t$. A comparison in \emph{mean square error} (MSE) with respect to the reference trajectory is given in the following table:
\begin{center}
\begin{tabular}{ |c|c|c|c|c| }
\hline
Agent & state & tube MPC & shallow MPC & deep MPC \\
\hline
\multirow{2}{*}{First} & $\delta$ & 15.8961   & 3.7257 &   0.0611 \\
& $p$ & 4.5280    & 1.3179 &    0.1025 \\
\hline
\multirow{2}{*}{Second} & $\delta$ & 20.6749   & 3.6843 &   0.0049 \\
& $p$ & 5.2877  &  1.0997  &  0.0291 \\
\hline
\end{tabular}
\end{center}
	\begin{figure}
		\centering
		\begin{adjustbox}{width = \columnwidth}
			\includegraphics{images/roll_angle.pdf}
		\end{adjustbox}
		\caption{Roll angle (deg.) versus time steps}
		\label{fig:roll_angle}
	\end{figure}

	\begin{figure}
		\centering
		\begin{adjustbox}{width = \columnwidth}
			\includegraphics{images/roll_rate.pdf}
		\end{adjustbox}
		\caption{Roll rate (deg./sec.) versus time steps}
		\label{fig:roll_rate}
	\end{figure}
\par The above table shows that deep MPC has significant improvement over tube MPC and shallow MPC in terms of MSE.
	\section{Epilogue}\label{s:epilogue}
	A deep learning based control algorithm is presented for safety critical systems by combining the approaches of adaptive control, generative neural network and tube MPC. The proposed approach is advantageous over tube MPC when uncertainties are matched. It is demonstrated with the help of a numerical experiment that our approach with a single layer neural network has mean square error smaller than that of tube MPC but is not able to adapt in the presence of abrupt change in uncertainties, whereas the deep MPC with a very small network is able to adapt well. Our approach is conservative in the sense that constraint tightening is required to satisfy constraints. Therefore, it is worth investigating along the lines of \cite{deep_tube_learning, dynamic_tube_How}. Some interesting extensions of the present approach may incorporate output feedback, unreliable channels \cite{PDQ_intermittent, ref:PSC18} and/or interacting agents. Since MPC experiences a bounded disturbance, we are able to ensure input-to-state stability of the closed-loop states. Some modifications in MPC may ensure asymptotic stability for small gain type reconstruction error by choosing some suitable common Lyapunov function.

	\bibliographystyle{IEEEtran}