Unmatched uncertainty mitigation through neural network supported model
  predictive control
This paper presents a deep learning based model predictive control (MPC)
algorithm for systems with unmatched and bounded state-action dependent
uncertainties of unknown structure. We utilize a deep neural network (DNN) as
an oracle in the underlying optimization problem of learning based MPC (LBMPC)
to estimate unmatched uncertainties. Generally, non-parametric oracles such as
DNN are considered difficult to employ with LBMPC due to the technical
difficulties associated with estimation of their coefficients in real time. We
employ a dual-timescale adaptation mechanism, where the weights of the last
layer of the neural network are updated in real time while the inner layers are
trained on a slower timescale using the training data collected online and
selectively stored in a buffer. Our results are validated through a numerical
experiment on the compression system model of jet engine. These results
indicate that the proposed approach is implementable in real time and carries
the theoretical guarantees of LBMPC.
Mateus V. Gasparino, Prabhat K. Mishra, Girish Chowdhary
# Unmatched uncertainty mitigation through neural network supported model predictive control

###### Abstract

This paper presents a deep learning based model predictive control (MPC) algorithm for systems with unmatched and bounded state-action dependent uncertainties of unknown structure. We utilize a deep neural network (DNN) as an oracle in the underlying optimization problem of learning based MPC (LBMPC) to estimate unmatched uncertainties. Generally, non-parametric oracles such as DNN are considered difficult to employ with LBMPC due to the technical difficulties associated with estimation of their coefficients in real time. We employ a dual-timescale adaptation mechanism, where the weights of the last layer of the neural network are updated in real time while the inner layers are trained on a slower timescale using the training data collected online and selectively stored in a buffer. Our results are validated through a numerical experiment on the compression system model of jet engine. These results indicate that the proposed approach is implementable in real time and carries the theoretical guarantees of LBMPC.

 safety critical systems, deep learning, learning based MPC

## I Introduction

Machine learning and Model Predictive Control (MPC) complement each other by compensating drawbacks of each other and making their combination useful for safety critical applications. We refer readers to [1, 2] for excellent surveys on safe learning and robotics.

Learning based Model Predictive Control (LBMPC) [3] became popular due to its improvement over linear MPC in terms of transient response and overshoot with slight expense in processing time [4]. Several interesting applications such as autonomous driving [5, 6], heat, ventilation and air-conditioning systems [7], quad-copter [4], formation control [8], atmospheric pressure plasma jets [9], air-borne wind energy systems [10] etc., contributed in theoretical and practical advancements.

Due to the significant success of and progress in deep learning techniques, this is tempting to use a neural network as an oracle in LBMPC. However, boundedness and differentiability of the oracle is required to hold the results of [3]. Even though a bounded and differentiable neural network can be constructed, the estimation of their weights in real time is generally difficult [11]. Since training of DNN is a time demanding process, real time implementation of DNN supported LBMPC is challenging.

In this article, we demonstrate that the real time implementation of neural network based oracle is not only possible but also computationally efficient by two time scale training. In [12, 13], the training of the output layer and the hidden layers are separated in which the output layer is trained online through adaptation by considering the recently trained hidden layers as a feature basis function and the hidden layers are trained on a parallel machine by keeping the recently updated weights of the output layer fixed. We refer readers to excellent surveys on transfer learning [14, 15]. This method of two time scale training allows us to try different methods of training the output layer and hidden layers independently. In addition, different architectures such as Recurrent Neural Network and Convolutional Neural Network are also supported for the hidden layer. The theoretical guarantees mostly depend on last activation layer and the training mechanism of the linear output layer. The above methods [12, 13] are limited only to those uncertainties that enter into the dynamics through control channel. The present article extends the results of [12, 13] by generalizing the class of uncertainties. Our approach improves the performance of [3] by replacing the L2NW estimator by a DNN. Since PyTorch is popular for DNN implementation and CasAdi for MPC, we address some non-trivial issues related to their interface also.

This article is organized as follows. The problem statement is given in SSII. The architecture of neural network and its

Fig. 1: The neural network on the main loop with fixed hidden layers is connected with MPC and transmits the weights of output layer at time sequence \((T_{k})_{k\in\mathbb{Z}_{*}}\) to the neural network on the training thread. This neural network (on the training thread) transmits the weights of hidden layer to the neural network (on the main loop) at time instants \((t_{f})_{j\in\mathbb{Z}_{*}}\) after its \(j^{\text{th}}\) training.
training mechanism are explained in SSIII. LBMPC and its properties are presented in SSIV and SSV, respectively. We validate our results in SSVI and conclude in SSVII.

We let \(\mathds{R}\) denote the set of real numbers, \(\mathds{N}\) 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 0\), \(\left\|v\right\|_{M}^{2}\) is used to denote \(v^{\top}Mv\). For a given matrix \(A\), the trace, the largest eigenvalue and pseudo-inverse are denoted by \(\operatorname{tr}(A)\), \(\lambda_{\max}(A)\) and \(A^{\dagger}\), respectively. By notation \(\left\|A\right\|\) and \(\left\|A\right\|_{\infty}\), we mean the \(2-\)norm and \(\infty-\)norm when \(A\) is a vector; induced \(2-\)norm and \(\infty-\)norm when \(A\) is a matrix, respectively. A vector or a matrix with all entries \(0\) is represented by \(\mathbf{0}\) and \(I\) is the identity matrix of appropriate dimensions. We let \(M^{(i)}\) denote the \(i^{\text{th}}\) column of a given matrix \(M\).

## II Problem setup

Let us consider a discrete time dynamical system

\[x_{t+1}=Ax_{t}+Bu_{t}+h(x_{t},u_{t});\quad t\in\mathds{N}, \tag{1}\]

1. \(x_{t}\in\mathcal{X}\subset\mathds{R}^{d}\), \(u_{t}\in\mathds{U}\subset\mathds{R}^{m}\); \(d,m\in\mathds{Z}_{+}\),
2. \(h:\mathcal{X}\times\mathds{U}\rightarrow\mathds{R}^{d}\) is a continuous and (possibly) non-linear function, which represents state-action dependent unmatched uncertainty (or modeling error),
3. \(h(x,u)\in\mathds{W}\) for every \(x\in\mathcal{X}\) and \(u\in\mathds{U}\),
4. \(\mathcal{X},\mathds{U}\) and \(\mathds{W}\) are polytopes,
5. the matrix pair \((A,B)\) is stabilizable.

The matrices \(A\) and \(B\) represent our domain knowledge or prior knowledge about the system dynamics, and the continuous function \(h\) represents the unknown component of the system dynamics1. LBMPC [3] has been developed to improve the closed-loop performance of tube-based robust MPC by modifying the cost function in the underlying optimization problem. The cost function is modified by learning (or estimating) the unknown function \(h\) with the help of data. In this article, we address issues associated with the use of DNN as an estimator of \(h\). Our main focus is to use DNN in such a way that it can be implemented in real time on a hardware with limited computational power. The general problem description of this article is as follows:

Footnote 1: For example, when a non-linear dynamics is linearized by the matrix pair \((A,B)\), \(h\) represents the linearization error.

**Problem Statement 1.** Present a stabilizing, robust and real-time implementable control framework for (1), which respects physical constraints (1-a), optimizes a given performance index, and reduces the effect of unmatched uncertainties by using a trainable DNN.

## III Deep neural network

Any continuous function \(h\) on a compact set \(\mathcal{X}\times\mathds{U}\) can be approximated by a multi-layer network with number of layers \(L\succeq 2\) such that

\[h(x,u)=W_{L}^{\top}\psi_{L}\left[W_{L-1}^{\top}\psi_{L-1}\left[\cdots[\psi_{ 1}(x,u)]\right]\right]+\varepsilon^{*}(x,u), \tag{2}\]

where \(x\in\mathcal{X},u\in\mathds{U},\psi_{i},W_{i}\) for \(i=1,\ldots,L\), are the activation functions in the \(i^{\text{th}}\) layer and the corresponding ideal weights, respectively [16, SS7.1]. We can represent \(h(x_{t},u_{t})\) with the help of a neural network with a desired accuracy.

Let us define \(\phi^{*}(x,u)\coloneqq\psi_{L}\left[W_{L-1}^{\top}\psi_{L-1}\left[\cdots[\psi_ {1}(x,u)]\right]\right]\) and \(W^{*}\coloneqq W_{L}\), then

\[h(x_{t},u_{t})=W^{*\top}\phi^{*}(x_{t},u_{t})+\varepsilon^{*}(x_{t},u_{t}), \tag{3}\]

where \(W^{*}\in\mathds{R}^{(n_{L}+1)\times d}\) denotes the weights of the output layer. There are \(n_{L}\) number of neurons in the last hidden layer. The first row of \(W^{*}\) represents the bias term in the output layer and the first element of \(\phi^{*}\in\mathds{R}^{n_{L}+1}\) is \(1\).

The major challenge is associated with real-time implementation of DNN because its training takes much more time than the sampling interval of fast hardware like quad-copter. Therefore, we address this issue with the help of two DNN's as shown in Fig. 1. Both DNNs have same architecture as in Fig. 2. The DNN in the main loop is located on the main machine and the other DNN is located on some secondary (or remote) machine. We update the weights of the output layer on the main machine in real time at each time instant with the help of a weight update law while keeping the weights of hidden layers fixed. The hidden layers are trained on a secondary machine by using the approach [17] in which the weights of the output layer are copied from the main machine at the start of the training and remain fixed during the training. Once the training of DNN on a secondary machine is complete, new weights of hidden layers are updated on the main machine and remain fixed until the new set of weights are again obtained from the secondary machine. Training details about the output layer and hidden layers are provided in SSIII-A and SSIII-B, respectively.

Let \((t_{j})_{j\in\mathds{Z}_{+}}\) represent an increasing time sequence of instants when weights of hidden layers are updated on main machine. Therefore, we get the following expression:

\[h(x_{t},u_{t})=W^{*\top}\phi_{j}(x_{t},u_{t})+\varepsilon_{j}(x_{t},u_{t}), \tag{4}\]

where \(\varepsilon_{j}(x_{t},u_{t})=W^{*\top}\left(\phi^{*}(x_{t},u_{t})-\phi_{j}(x_{t },u_{t})\right)+\varepsilon^{*}(x_{t},u_{t})\). We can assume that \(\left\|\phi_{j}(x,u)\right\|\) will be bounded for each \(j\), \(u\in\mathds{R}^{m}\)

Fig. 2: The neural network architecture can use an arbitrary number of layers and neurons. The input layer has the same number of neurons as the experience data. The output layer has as many neurons as system states. The hidden layers can have any architecture provided suitable adjustments are made to facilitate their training.
and \(x\in\mathds{R}^{d}\) due to the presence of bounded activation layer in Fig. 2 consisting bounded neurons, i. e. sigmoidal, tanh, etc. Since \(h\) is bounded due to (1-c), we can assume that ideal weights \(W^{*}\) in the output layer are also bounded. We make the following assumption:

**Assumption 1**.: There exist \(\bar{W}_{i}>0\) for \(i=1,\ldots,d\), and \(\sigma,\bar{\epsilon}>0\) such that \(\left\|W^{*(i)}\right\|\leqslant\bar{W}_{i}\), for \(i=1,\ldots,d\), and \(\left\|\phi_{j}(x,u)\right\|\leqslant\sigma\) for every \(x\in\mathcal{X},u\in\mathds{U}\) and \(j\in\mathds{N}\).

The above assumption is standard in literature [19, 12, 18]. 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 [16, SS7.1]. A priori knowledge about the bounds on the ideal weights \(W^{*}\) of the output layer is useful to avoid parameter drift phenomenon.

At \(t_{0}=0\), weights of DNN on both machines are randomly initialized with desired bound on the output layer weights as per the Assumption 1. Therefore, for \(j\in\mathds{N}\), we have

\[\hat{h}_{t}(x_{t},u_{t})\coloneqq K_{t}^{\top}\phi_{j}(x_{t},u_{t})\text{ for }t\in\{t_{j},t_{j}+1,\ldots,t_{j+1}-1\}. \tag{5}\]

We update \(K_{t}\) in an unsupervised manner on main machine while collecting the training data for the DNN on secondary machine. In the next subsections we provide the relevant details of the training of DNN.

### _Adaptive learning of \(W^{*}\) on the main machine_

For \(t\in\{t_{j},t_{j}+1,\ldots,t_{j+1}-1\}\), the output of DNN is given by (5) and the bounded features are given by \(\phi_{j}(x_{t},u_{t})\) at time \(t+1\)2. We get the estimated state

Footnote 2: This is important for implementation to note that \(u_{t}\) is computed at time \(t\) but we need to use \(\phi_{j}(x_{t},u_{t})\) and \(\hat{h}(x_{t},u_{t})\) at time \(t+1\).

\[\hat{x}_{t+1}=Ax_{t}+Bu_{t}+\hat{h}_{t}(x_{t},u_{t}). \tag{6}\]

We can compute the error \(\bar{x}_{t+1}=\hat{x}_{t+1}-x_{t+1}\). For the online training of the output layer, we employ the projection based robust weight update law and refer readers to [20, Chapter 10] for other methods. For a given learning rate \(0<\gamma<1\), we first modify \(K_{t}\) by taking \(\bar{x}_{t+1}\) into account to get \(\bar{K}_{t+1}\) and then project \(\bar{K}_{t+1}\) to ensure its boundedness as follows:

\[\bar{K}_{t+1} =K_{t}-\gamma\frac{\phi_{j}(x_{t},u_{t})}{\left\|\phi_{j}(x_{t},u _{t})\right\|^{2}}\bar{x}_{t+1}^{\top}, \tag{7}\] \[K_{t+1}^{(t)} =\text{Proj}\bar{K}_{t+1}^{(t)}=\begin{cases}\bar{K}_{t+1}^{(t)}& \text{if }\left\|\bar{K}_{t+1}^{(t)}\right\|\leqslant\bar{W}_{i}\\ \frac{\bar{W}_{i}}{\left\|\bar{K}_{t+1}^{(t)}\right\|}\bar{K}_{t+1}^{(t)}& \text{otherwise}.\end{cases}.\]

The implications of (7) are discussed in SSV.

### _Supervised learning of \(\phi^{*}\) on a secondary machine_

In this section, the training data and the loss function required for the supervised training of the DNN on secondary machine are explained. For a given state-action pair \((x_{t},u_{t})\) as input at time \(t\), the label \(h(x_{t},u_{t})\) is computed at time \(t+1\) by the relation:

\[h(x_{t},u_{t})=x_{t+1}-Ax_{t}-Bu_{t}. \tag{8}\]

These pairs \((x,u)\) and corresponding labels \(h(x,u)\) are stored in a training buffer until the buffer is full. In particular, a full buffer consists of \(\left((x^{i},u^{i}),h(x^{i},u^{i})\right)\) for \(i=1,\ldots,p_{\max}\).

Once the buffer is full, new data is stored by replacing some old data. Several methods are proposed in the literature to increase the richness of the training buffer; see [21, 22, 23] and references therein. We refer readers to [24] for different methods of designing replay buffer.

Let \((T_{k})_{k\in\mathcal{Z}_{u}}\) be an increasing time sequence. At time \(T_{k}\), the weight of output layer of the primary neural network are copied in the secondary neural network. During the training, the output layer weights in the secondary network remain fixed. Therefore, we are interested in finding the weights \(W_{1:L-1}\coloneqq W_{1},\ldots,W_{L-1}\) which minimize the following cost for a given input \((x,u)\) and label \(h(x,u)\):

\[\ell\left((x,u),W_{1:L-1}\right)\coloneqq\\ \left\|h(x,u)-K_{T_{k}}^{\top}\psi_{L}\left[W_{L-1}^{\top}\psi_{L-1} \left[\cdots\left[\psi_{1}(x,u)\right]\right]\right]\right\|^{2}.\]

Let \(M\) represent the number of training samples and \(\mathcal{D}_{k}\coloneqq\left((x^{i},u^{i}),h(x^{i},u^{i})\right)_{i=1}^{M}\) is training data consisting \(M\) data points randomly sampled from the buffer for the \(k^{\text{th}}\) training. The following loss function is considered for the training of DNN:

\[\mathcal{L}(\mathcal{D}_{k},W_{1:L-1})=\frac{1}{M}\sum_{i=0}^{M}\ell\left((x^{i },u^{i}),W_{1:L-1}\right).\]

## IV Model predictive controller

We first fix an optimization horizon \(N\in\mathds{Z}_{+}\). Let \(x^{r},u^{r}\) be a reference (or equilibrium) state-action pair. For given positive definite matrices \(Q,R>0\), \(P>0\) is the solution of the following Lyapunov equation:

\[(A+BK)^{\top}P(A+BK)-P=-(Q+K^{\top}RK), \tag{9}\]

where \(K\) is such that \(A+BK\) is Schur stable. We define the following cost

\[\psi\left(z_{0:N+1},u_{0:N}\right)\coloneqq\left\|z_{N}-x^{r}\right\|_{P}^{2}+ \sum_{i=0}^{N-1}\left\|z_{i}-x^{r}\right\|_{Q}^{2}+\left\|v_{i}-u^{r}\right\|_{R }^{2}.\]

Let us define \(R_{i+1}=(A+BK)R_{i}\oplus\mathds{W}\) with \(R_{0}=\{0\}\). For \(i=0,\ldots,N-1\), we impose the following constraints:

\[\bar{z}_{i+1} =A\bar{z}_{i}+Bv_{i} \tag{10}\] \[\bar{z}_{i} \in\mathcal{X}\odot R_{i},\quad v_{i}\in\mathds{U}\odot KR_{i}\] \[\bar{z}_{N} \in\Omega\odot R_{N},\]

where \(\Omega\) is a disturbance invariant set. At each time \(t\), we measure the state \(x_{t}\) of the system (1) and solve the following optimization problem:

\[\min_{(c_{i})_{i=0}^{N-1}} \psi\left(z_{0:N+1},v_{0:N}\right)\] s. t. \[z_{0}=\bar{z}_{0}=x_{t}\] \[v_{i} =K\bar{z}_{i}+c_{i}\text{ for }i\in\mathds{Z}_{[0,N-1]}\] (11) \[z_{i+1} =Az_{i}+Bv_{i}+\hat{h}_{t}(z_{i},v_{i})\text{ for }i\in\mathds{Z}_{[0,N-1]}\] (10).

By solving the above problem, we get \(v_{0}=Kx_{t}+c_{0}\). We set \(u_{t}=v_{0}\) and apply \(u_{t}\) to the system (1).
## V Properties of LBMPC

For the purpose of analysis, we define \(\tilde{K}_{t}:=K_{t}-W^{*},\tilde{x}_{t+1}:=\tilde{x}_{t+1}-x_{t+1}=\tilde{h}_{t}( x_{t},u_{t})-h(x_{t},u_{t})=\tilde{K}_{t}^{\top}\phi_{J}(x_{t},u_{t})-\varepsilon_{J}(x _{t},u_{t})\) and \(\tilde{W}:=\sum_{t=1}^{m}\tilde{W}_{t}^{2}\). We recall the following definition:

**Definition 1** ([25], page 117).: The vector sequence \((s_{t})_{t\in\mathbb{N}}\) is called \(\mu\) small in mean square sense if it satisfies \(\sum_{t=k}^{k+N-1}\|s_{t}\|^{2}\leqslant Nc_{0}\mu+c_{0}^{\prime}\) for all \(k\in\mathbb{Z}_{+}\), a given constant \(\mu\geqslant 0\) and some \(N\in\mathbb{Z}_{+}\), where \(c_{0},c_{0}^{\prime}\geqslant 0\).

We make the following assumption:

**Assumption 2**.: There exists \(\bar{\varepsilon}>0\) such that

\[\left\|\varepsilon_{j}(x,u)\right\|\leqslant\bar{\varepsilon}\,\,\,\text{for each}\,\,\,(x,u)\in\mathcal{X}\times\mathbb{U}\,\,\text{and}\,\,\,j\in\mathbb{N}.\]

Since \(h\) is bounded due to (1-c), \(W^{*},\phi_{j}\) due to Assumption 1 and \(K_{t}\) due to (7), the above assumption is trivially satisfied. We have the following result that says that the estimation error \(\tilde{x}_{t+1}=\hat{h}_{t}(x_{t},u_{t})-h(x_{t},u_{t})\) is small in mean square sense.

**Lemma 1**.: _Consider the dynamical system (1), weight update law (7). Let the assumptions 1 and 2 hold. Let us define \(V_{a}(K_{t}):=\frac{1}{\gamma}\operatorname{tr}(\tilde{K}_{t}^{\top}\tilde{K} _{t})\). Then for all \(t\),_

1. \(V_{a}(K_{t})\leqslant\frac{4}{\gamma}\tilde{W}\)_,_
2. \(V_{a}(K_{t+1})-V_{a}(K_{t})\leqslant-\frac{1-\gamma}{\sigma^{2}}\left\|\tilde{ x}_{t+1}\right\|^{2}+\left\|\varepsilon_{j}(x_{t},u_{t})\right\|^{2}\)_,_
3. \(\tilde{x}_{t}\)_is_ \(\bar{\varepsilon}^{2}\) _small in mean square sense with_ \(c_{0}=\frac{\sigma^{2}}{1-\gamma}\) _and_ \(c_{0}^{\prime}=\frac{4c_{0}}{\gamma}\tilde{W}\) _as per the Definition_ 1_._

Proof.: (Proof of Lemma 1)

1. Since \(V_{a}(K_{t})=\frac{1}{\gamma}\operatorname{tr}(\tilde{K}_{t}^{\top}\tilde{K} _{t})=\frac{1}{\gamma}\sum_{i=1}^{m}\left\|K_{t}^{(i)}-W^{*(i)}\right\|^{2} \leqslant\frac{2}{\gamma}\sum_{i=1}^{m}\left\|K_{t}^{(i)}\right\|^{2}+\left\| W^{*(i)}\right\|^{2}\leqslant\frac{4}{\gamma}\sum_{i=1}^{m}\tilde{W}_{t}^{2}= \frac{4}{\gamma}\tilde{W}\).
2. By substituting \(\tilde{K}_{t+1}=\tilde{K}_{t+1}-W^{*}+K_{t+1}-\tilde{K}_{t+1}\) in \(V_{a}(K_{t+1})\) and defining \[\alpha_{t} \coloneqq(K_{t+1}-\tilde{K}_{t+1})^{\top}(K_{t+1}-\tilde{K}_{t+1})\] \[\quad+2(K_{t+1}-\tilde{K}_{t+1})^{\top}(\tilde{K}_{t+1}-W^{*})\] \[=-(K_{t+1}-\tilde{K}_{t+1})^{\top}(\tilde{K}_{t+1}-\tilde{K}_{t+1})\] \[\quad+2(K_{t+1}-K_{t+1})^{\top}(K_{t+1}-W^{*}),\] we get \[V_{a}(K_{t+1}) =\frac{1}{\gamma}\operatorname{tr}(\tilde{K}_{t+1}^{\top}\tilde{K }_{t+1})\] \[=\frac{1}{\gamma}\operatorname{tr}\left((\tilde{K}_{t+1}-W^{*})^{ \top}(\tilde{K}_{t+1}-W^{*})\right)+\frac{1}{\gamma}\operatorname{tr}\left( \alpha_{t}\right).\] One important property of the projection is the following [25, (4.61)]: \[(W^{*(i)}-K_{t}^{(i)})^{\top}(\tilde{K}_{t}^{(i)}-K_{t}^{(i)})\leqslant 0 \,\,\text{for each}\,\,i=1,\ldots,m.\] (12) Since \((K_{t+1}^{(i)}-\tilde{K}_{t+1}^{(i)})^{\top}(K_{t+1}^{(i)}-W^{*})\leqslant 0\) due to (12), we can ensure \(\operatorname{tr}(\alpha_{t})\leqslant 0\). Therefore, \[V_{a}(K_{t+1})\leqslant\frac{1}{\gamma}\operatorname{tr}\left(( \tilde{K}_{t+1}-W^{*})^{\top}(\tilde{K}_{t+1}-W^{*})\right)\] \[=V_{a}(K_{t})+\frac{1}{\left\|\phi_{J}(x_{t},u_{t})\right\|^{2}} \operatorname{tr}\left(\tilde{x}_{t+1}\tilde{x}_{t+1}^{\top}\right)\] \[\quad-\frac{1}{\left\|\phi_{J}(x_{t},u_{t})\right\|^{2}} \operatorname{tr}\left(\tilde{K}_{t}^{\top}\phi_{J}(x_{t},u_{t})\tilde{x}_{t+1} ^{\top}+\tilde{x}_{t+1}\phi_{J}(x_{t},u_{t})^{\top}\tilde{K}_{t}\right).\] By substituting \(\tilde{K}_{t}^{\top}\phi_{J}(x_{t},u_{t})=\tilde{x}_{t+1}+\varepsilon_{j}(x_{t}, u_{t})\) in the above inequality, we get \[V_{a}(K_{t+1})\leqslant V_{a}(K_{t})+\frac{1}{\left\|\phi_{J}(x_{t},u_{t}) \right\|^{2}}\left(\gamma\left\|\tilde{x}_{t+1}\right\|^{2}-2\operatorname{tr}( \tilde{x}_{t+1}\] \[\quad+\varepsilon_{j}(x_{t}))\tilde{x}_{t+1}^{\top}\right)\] \[=V_{a}(K_{t})+\frac{1}{\left\|\phi_{J}(x_{t},u_{t})\right\|^{2}} \left(\gamma\left\|\tilde{x}_{t+1}\right\|^{2}-2\tilde{x}_{t+1}^{\top} \varepsilon_{j}(x_{t},u_{t})\right)\] \[\leqslant V_{a}(K_{t})+\frac{1}{\left\|\phi_{J}(x_{t},u_{t}) \right\|^{2}}\left((\gamma-1)\left\|\tilde{x}_{t+1}\right\|^{2}+\left\| \varepsilon_{j}(x_{t},u_{t})\right\|^{2}\right)\] \[=V_{a}(K_{t})-\frac{1-\gamma}{\left\|\phi_{J}(x_{t},u_{t})\right\|^{2 }}\left\|\tilde{x}_{t+1}\right\|^{2}+\left\|\varepsilon_{j}(x_{t},u_{t}) \right\|^{2}\] \[\leqslant V_{a}(K_{t})-\frac{1-\gamma}{\sigma^{2}}\left\|\tilde{x}_{t+ 1}\right\|^{2}+\left\|\varepsilon_{j}(x_{t},u_{t})\right\|^{2},\] where the last inequality is due to \(1\leqslant\left\|\phi_{J}(x_{t},u_{t})\right\|^{2}\leqslant\sigma^{2}\). Therefore, \[V_{a}(K_{t+1})-V_{a}(K_{t})\leqslant-\frac{1-\gamma}{\sigma^{2}}\left\|\tilde{x}_{t+ 1}\right\|^{2}+\left\|\varepsilon_{j}(x_{t},u_{t})\right\|^{2}.\]
3. Consider Lemma 1-(ii) to get \[\frac{1-\gamma}{\sigma^{2}}\left\|\tilde{x}_{t+1}\right\|^{2} \leqslant-V_{a}(K_{t+1})+V_{a}(K_{t})+\left\|\varepsilon_{j}(x_{t},u_{t}) \right\|^{2}\] \[\leqslant-V_{a}(K_{t+1})+V_{a}(K_{t})+\bar{\varepsilon}^{2}.\] By summing from \(t=k\) to \(k+N-1\) in both sides, we get \[\frac{1-\gamma}{\sigma^{2}}\sum_{t=k}^{k+N-1}\|\tilde{x}_{t+1}\|^{2} \leqslant V_{a}(K_{k})+N\bar{\varepsilon}^{2},\] \[\leqslant\frac{4}{\gamma}\tilde{W}+N\bar{\varepsilon}^{2}.\] Therefore, \(\tilde{x}_{t}\) is \(\bar{\varepsilon}^{2}\) small in mean square sense with \(c_{0}=\frac{\sigma^{2}
Theorem 2]. The second result is a direct implication of [27, Lemma 3.5]. The recursive feasibility of (11) is due to [3, Theorem 1]. In particular, if \((c_{0},c_{1},\ldots,c_{N-1})\) is an optimizer of (11) at some time \(t\) then \((c_{1},\ldots,c_{N-1},0)\) will be a feasible solution at time \(t+1\). 

## VI Numerical experiment

We consider the compression system of a jet engine, which exhibits instabilities due to rotating stall and surge. The Moore-Greitzer compressor model is given by the following nonlinear dynamics:

\[\dot{z} =-y+z_{c}+1+\frac{3}{2}z-\frac{1}{2}z^{3} \tag{13}\] \[\dot{y} =\frac{1}{\beta^{2}}\left(z+1-r\sqrt{z}\right),\]

where \(z\) is mass flow, \(y\) is pressure rise, \(\beta>0\) is a constant and \(r\) is the throttle opening. Similar to [3], we assume that \(r\) is controlled by a second order actuation with transfer function \(r(s)=\frac{\omega_{n}^{2}}{s^{2}+2\zeta\omega_{n}s+\omega_{n}^{2}}u(s)\), where \(u\) is the input. Our simulation parameters are \(\beta=1,z_{c}=0,\zeta=\frac{1}{\sqrt{2}},\omega_{n}=10\sqrt{10}\). We have the following constraints:

\[z \in[0,1] y\in[1.1875,2.1875] \tag{14}\] \[r \in[0.1547,2.1547] \dot{r}\in[-20,20]\] \[u \in[0.1547,2.1547]\]

The state \(x\) of the system is represented by \(x=\begin{bmatrix}z&y&r\end{bmatrix}^{\top}\in\mathbb{R}^{4}\). The nonlinear model (13) is linearized around the equilibrium point \(x_{e}=[0.5&1.6857&1.1547&0]^{\top}\) and discretized with sampling time \(T=0.05\) seconds. We chose \(T=0.05\) to make it slightly larger than the solver time in one optimization problem corresponding to linear MPC for its online implementation. We use a learning rate of \(0.001\) to train the deep neural network in parallel.

We employed CasADI for the symbolic representation of the underlying optimization problem of MPC (11). This symbolic representation is done offline and therefore, helpful in online implementation of MPC by updating only the measured state. Since PyTorch has inbuilt tools for the training of neural networks, we want to use it along with CasADI. We design a neural network on main machine with the help of CasADI and that on the secondary machine by PyTorch.

Another difficulty is associated with the computation of terminal set \(\Omega\) and reachable sets \(R_{i}\)'s. The computation of \(R_{i}\) is very hard for large \(i\). MATLAB based tools like MPT are not available in python for polytopic manipulation. The python library pytone has very limited capability and is under development. Therefore, we followed the approach of [28]. Let the set \(\mathcal{X}\ominus\mathrm{W}\), \(\mathcal{X}\) and \(\mathrm{U}\) be given by

\[\mathcal{X}\ominus\mathrm{W} \coloneqq\{x\in\mathbb{R}^{d}\mid F_{p}x\leqslant h_{p}\}\] \[\mathcal{X} \coloneqq\{x\in\mathbb{R}^{d}\mid F_{xx}\leqslant h_{x}\}\] \[\mathrm{U} \coloneqq\{x\in\mathbb{R}^{m}\mid F_{ux}\leqslant h_{u}\},\]

where \(F_{p},F_{x},F_{u}\) are suitable matrices and \(h_{p},h_{x},h_{u}\) are suitable vectors required to provide the half-space representations of above sets. The set \(\Omega\) and \(R_{i}\) are approximated in [28].

We are interested in comparing our proposed approach with the state-of-the-art [3], which employs L2NW as an estimator. We implemented L2NW with the help of CasADI in terms of symbolic variables. One of the inputs for L2NW is a data buffer of fixed size, used as function parameters to learn the uncertainties. In an online implementation, this data buffer is empty at the beginning and its size increases with time. Therefore, we initialize a fixed size buffer in such a way that default values do not affect the outcome of the L2NW estimator. Our simulation parameters are same as in [3] except the sampling time \(T\). Figures 3 and 4 demonstrate that the proposed approach has faster transient response than that of linear MPC and smaller overshoot than that of [3]. It is demonstrated in Fig. 5 that the proposed approach and linear MPC both have comparable solver time but that in [3] is larger than the sampling time at the beginning.

Only drawback in [3] is the use of L2NW estimator, which makes the underlying optimization problem of LBMPC computationally demanding and therefore hard to implement on a small machines. Since the approach [3] gives choice to use any estimator, authors used linear estimators in [4, 28] instead of L2NW estimator. Our approach of using DNN as an estimator is not only computationally less demanding due to the parallel processing but surprisingly the underlying optimization problems have similar solver time as linear MPC (Fig. 5). Therefore, the proposed approach can be implemented on small machines.

## VII Conclusion

This article demonstrates that DNN can be used in the framework of LBMPC by dual time scale training and using bounded neurons in the last activation layer. The proposed approach is able to provide a faster solution because a DNN can be trained on a separate machine. Since LBMPC allows any estimator as an oracle, different DNN architectures can be investigated for different class of problems by following the proposed approach. Some interesting extensions may be possible along the lines of vision based navigation [29], stochastic MPC [30, 31], Bayesian Neural Network [32] and crystallization processes [33, 34].

Fig. 3: Proposed approach and [3] both have faster transient response than that of linear MPC but the proposed approach has smaller overshoot than that of [3].