Error bounds for maxout neural network approximations of model
  predictive control
Neural network (NN) approximations of model predictive control (MPC) are a
versatile approach if the online solution of the underlying optimal control
problem (OCP) is too demanding and if an exact computation of the explicit MPC
law is intractable. The drawback of such approximations is that they typically
do not preserve stability and performance guarantees of the original MPC.
However, such guarantees can be recovered if the maximum error with respect to
the optimal control law and the Lipschitz constant of that error are known. We
show in this work how to compute both values exactly when the control law is
approximated by a maxout NN. We build upon related results for ReLU NN
approximations and derive mixed-integer (MI) linear constraints that allow a
computation of the output and the local gain of a maxout NN by solving an MI
feasibility problem. Furthermore, we show theoretically and experimentally that
maxout NN exist for which the maximum error is zero.
Dieter Teichrib, Moritz Schulze Darup
# Error bounds for maxout neural network

###### Abstract

Neural network (NN) approximations of model predictive control (MPC) are a versatile approach if the online solution of the underlying optimal control problem (OCP) is too demanding and if an exact computation of the explicit MPC law is intractable. The drawback of such approximations is that they typically do not preserve stability and performance guarantees of the original MPC. However, such guarantees can be recovered if the maximum error with respect to the optimal control law and the Lipschitz constant of that error are known. We show in this work how to compute both values exactly when the control law is approximated by a maxout NN. We build upon related results for ReLU NN approximations and derive mixed-integer (MI) linear constraints that allow a computation of the output and the local gain of a maxout NN by solving an MI feasibility problem. Furthermore, we show theoretically and experimentally that maxout NN exist for which the maximum error is zero.

## 1 Introduction

Model predictive control (MPC) (see, e.g, (Rawlings et al., 2017)) has become a standard tool for the control of dynamical systems with state and input constraints and has been successfully applied in different industrial fields (see, e.g., (Qin and Badgwell, 2003) for an overview). In the classical setup, MPC requires to solve an optimization problem (OP) in every time step. For systems with a short sampling period as, e.g., in power electronics (Karamanakos et al., 2020), this can be challenging because the OP may be too complex to be solved within the sampling period. If we consider a linear discrete-time prediction model in combination with a quadratic cost function, the resulting OP is a quadratic program (QP). In principle, we can compute the solution of the parametric QP offline for all feasible states. This results in an explicit control law with a piecewise affine (PWA) input-output relation defined on a polyhedral partition of the state space (Bemporad et al., 2002). Given the explicit control law, the online computational effort reduces to the evaluation of the PWA function. However, the number of polyhedral regions may grow exponentially with the state dimension and the number of constraints in the OP. Hence, exactly computing the explicit MPC law becomes untractable for complex systems. As a consequence, various techniques have been developed to approximate the control law in MPC (Jones and Morari, 2009; Bemporad and Filippi, 2003). In this context, neural networks (NN) are very popular (Chen et al., 2018; Karg and Lucia, 2020; Chen et al., 2022) since they can approximate a large class of functions, including PWA functions, with arbitrary accuracy (Hornik et al., 1989). In addition, besides their computational demanding offline training, NN are typically fast to evaluate online, which is essential if they are used as controllers. Moreover, some types of NN share the PWA structure of the control law (Schulze Darup, 2020; Hanin, 2017; Arora et al., 2016), making them the perfect choice for approximating MPC. Unfortunately, in general, stability cannot be guaranteed for the approximated controllers. One possibility to recover stability and recursive feasibility is to project the output of the NN onto a suitable set as, e.g., in (Paulson and Mesbah, 2020; Chen et al., 2018). Alternatively, the output of the NN can be used as an initial guess for a solver and not directly for control (Chen et al., 2022). The drawback of both approaches is that they require an additional optimization-based computation step online. Ideally, the NN can be used as a controller without additional online computation, while still providing stability guarantees. In (Fabiani and Goulart, 2022), it is proven that this is possible if the maximum error with respect to the optimal control law and the Lipschitz constant of the corresponding error function are known. Both values can indeed be computed exactly presupposed the output and the local gain of the used NN can be computed by solving a mixed-integer (MI) feasibility problem. This is known to be possible for NN using rectified linear units (ReLU) as activation functions (Fabiani and Goulart, 2022, Thm. 6.1).

In the work at hand, we will extend the result of (Fabiani and Goulart, 2022) by showing that the maximum error and the Lipschitz constant of the error can also be computed exactly for a controller approximation based on a maxout NN. Since maxout NN include other NN with PWA input-output relation such as, e.g., ReLU and leaky ReLU, as a special case, the results provide a generalization to a broader class of PWA NN. Furthermore, we use the PWA structure of maxout NN to compute NN that exactly describe MPC control laws and validate experimentally that these exact maxout NN indeed lead to a maximum error and Lipschitz constant of zero.

The paper is organized as follows. In the remainder of this section, we introduce relevant notation. In Section 2, we summarize some basics on MPC as well as PWA NN and describe concepts
for approximating MPC in more depth. Section 3 is devoted to our main result, i.e., the computation of the maximum error and the Lipschitz constant of the error function related to maxout NN approximating MPC. The obtained method is applied to various maxout NN approximations of MPC laws in Section 4. Finally, conclusions and an outlook are given in Section 5.

### Notation

We will denote the index set containing \(p_{i}\in\mathbb{N}\) integers starting at \(p_{i}(l-1)+1,\;l\in\mathbb{N}\) by

\[\mathcal{A}_{l}^{(i)}:=\{p_{i}(l-1)+1,\ldots,p_{i}l\}.\]

For vectors \(\mathbf{x}\in\mathbb{R}^{n}\) we denote the \(i\)-th element by \(\mathbf{x}_{i}\) and the elements between the indices \(a_{1}\) and \(a_{2}>a_{1}\) by \(\mathbf{x}_{a_{1}:a_{2}}\). For matrices \(\mathbf{K}\in\mathbb{R}^{m\times n}\) we denote the element in the \(i\)-th row and \(j\)-th column by \(\mathbf{K}_{i,j}\), the \(i\)-row and \(j\)-th column by \(\mathbf{K}_{i,\cdot}\) and \(\mathbf{K}_{\cdot,j}\), respectively. If we only write \(\mathbf{K}_{i}\) then we refer to the \(i\)-th row. A block diagonal matrix is defined as

\[\mathrm{diag}(\mathbf{\alpha}_{1},\ldots,\mathbf{\alpha}_{N}):=\begin{pmatrix}\mathbf{ \alpha}_{1}&\mathbf{0}&\ldots&\mathbf{0}\\ \mathbf{0}&\ddots&&\vdots\\ \vdots&&&\mathbf{0}\\ \mathbf{0}&\ldots&\mathbf{0}&\mathbf{\alpha}_{w_{i}}\end{pmatrix},\]

with \(\mathbf{\alpha}\in\mathbb{R}^{w_{i}\times p_{i}}\). A continuous function \(\mathbf{F}(\mathbf{x}):\mathcal{P}\subset\mathbb{R}^{n}\to\mathbb{R}^{m}\) of the form

\[\mathbf{F}(\mathbf{x})=\left\{\begin{array}{ccc}\mathbf{G}^{(1)}\mathbf{x}+\mathbf{g}^{(1)}& \text{if}&\mathbf{x}\in\mathcal{P}^{(1)},\\ \vdots&&\vdots\\ \mathbf{G}^{(s)}\mathbf{x}+\mathbf{g}^{(s)}&\text{if}&\mathbf{x}\in\mathcal{P}^{(r)},\end{array}\right. \tag{1}\]

with a polyhedral partition \(\mathcal{P}=\cup_{i=1}^{r}\mathcal{P}^{(i)}\) and \(\mathrm{int}(\mathcal{P}^{(i)})\cap\mathrm{int}(\mathcal{P}^{(j)})=\emptyset \forall i\neq j\) is denoted as piecewise affine (PWA) function. We further define the local gain \(\mathbf{K}(\mathbf{x}):\cup_{i=1}^{r}\mathrm{int}(\mathcal{P}^{(i)})\to\mathbb{R}^{m}\) of a PWA function as

\[\mathbf{K}(\mathbf{x}):=\left\{\begin{array}{ccc}\mathbf{G}^{(1)}&\text{if}&\mathbf{x}\in \mathrm{int}(\mathcal{P}^{(1)}),\\ \vdots&&\vdots\\ \mathbf{G}^{(s)}&\text{if}&\mathbf{x}\in\mathrm{int}(\mathcal{P}^{(r)}).\end{array}\right. \tag{2}\]

## 2 Fundamentals of MPC and NN

### Model predicitive control

Model predictive control (MPC) for linear discrete-time systems builds on solving an optimal control problem (OCP) of the form

\[V_{N}(\mathbf{x}):= \min_{\hat{\mathbf{u}}(0),\ldots,\hat{\mathbf{u}}(N)\atop\hat{\mathbf{u}}(0 ),\ldots,\hat{\mathbf{u}}(N-1)}\varphi(\hat{\mathbf{x}}(N))+\sum_{\kappa=0}^{N-1} \ell(\hat{\mathbf{x}}(\kappa),\hat{\mathbf{u}}(\kappa))\] (3) s.t. \[\hat{\mathbf{x}}(0) =\mathbf{x},\] \[\hat{\mathbf{x}}(\kappa+1) =\mathbf{A}\,\hat{\mathbf{x}}(\kappa)+\mathbf{B}\hat{\mathbf{u}}(\kappa), \forall\kappa\in\{0,...,N-1\},\] \[(\hat{\mathbf{x}}(\kappa),\hat{\mathbf{u}}(\kappa)) \in\mathcal{X}\times\mathcal{U}, \forall\kappa\in\{0,...,N-1\},\] \[\hat{\mathbf{x}}(N) \in\mathcal{T}\]

in every time step \(k\in\mathbb{N}\) for the current state \(\mathbf{x}=\mathbf{x}(k)\). Here, \(N\in\mathbb{N}\) refers to the prediction horizon and

\[\varphi(\mathbf{x}):=\mathbf{x}^{\top}\mathbf{P}\mathbf{x}\quad\text{and}\quad\ell(\mathbf{x},\bm {u}):=\mathbf{x}^{\top}\mathbf{Q}\mathbf{x}+\mathbf{u}^{\top}\mathbf{Ru} \tag{4}\]

denote the terminal and stage cost, respectively, where the weighting matrices \(\mathbf{P}\), \(\mathbf{Q}\), and \(\mathbf{R}\) are positive (semi-) definite. The dynamics of the linear prediction model are described by \(\mathbf{A}\in\mathbb{R}^{n\times n}\) and \(\mathbf{B}\in\mathbb{R}^{n\times m}\). State and input constraints can be incorporated via the polyhedral sets \(\mathcal{X}\) and \(\mathcal{U}\). Finally, the terminal set \(\mathcal{T}\) allows to enforce closed-loop stability (see (Mayne et al., 2000) for details). The resulting control law \(\mathbf{\pi}:\mathcal{F}_{N}\to\mathcal{U}\) is defined as

\[\mathbf{\pi}(\mathbf{x}):=\hat{\mathbf{u}}^{*}(0), \tag{5}\]

where \(\mathcal{F}_{N}\) denotes the feasible set of (3) and where \(\hat{\mathbf{u}}^{*}(0)\) refers to the first element of the optimal input sequence. For the considered setup it is well known that \(\mathbf{\pi}(\mathbf{x})\) is a PWA function (Bemporad et al., 2002, Thm. 4) of the form (1) with \(\mathbf{G}^{(i)}=\mathbf{K}^{(i)},\mathbf{g}^{(i)}=\mathbf{b}^{(i)}\), \(r=r_{\text{MPC}}\), polyhedral sets \(\mathcal{P}^{(i)}=\mathcal{R}^{(i)}\;\forall i\in\{1,\ldots,r_{\text{MPC}}\}\) and local gain \(\mathbf{K}_{\text{MPC}}(\mathbf{x})\).

### Neural networks with piecewise affine activations

In general, a feed-forward-NN with \(\ell\in\mathbb{N}\) hidden layers and \(w_{i}\) neurons in layer \(i\) can be written as a composition of the form

\[\mathbf{\Phi}(\mathbf{x})=\mathbf{f}^{(\ell+1)}\circ\mathbf{g}^{(\ell)}\circ\mathbf{f}^{(\ell)} \circ\cdots\circ\mathbf{g}^{(1)}\circ\mathbf{f}^{(1)}(\mathbf{x}). \tag{6}\]

Here, the functions \(\mathbf{f}^{(i)}:\mathbb{R}^{w_{i-1}}\to\mathbb{R}^{p_{i}w_{i}}\) for \(i\in\{1,\ldots,\ell\}\) refer to preactivations, where the parameter \(p_{i}\in\mathbb{N}\) allows to consider "multi-channel" preactivations as required for maxout (see (Goodfellow et al., 2013)). Moreover, \(\mathbf{g}^{(i)}:\mathbb{R}^{p_{i}w_{i}}\to\mathbb{R}^{w_{i}}\) stand for activation functions and \(\mathbf{f}^{(\ell+1)}:\mathbb{R}^{w_{\ell}}\to\mathbb{R}^{w_{\ell+1}}\) reflects postactivation. The functions \(\mathbf{f}^{(i)}\) are typically affine, i.e.,

\[\mathbf{f}^{(i)}(\mathbf{y}^{(i-1)})=\mathbf{W}^{(i)}\mathbf{y}^{(i-1)}+\mathbf{b}^{(i)}, \tag{7}\]

where \(\mathbf{W}^{(i)}\in\mathbb{R}^{p_{i}w_{i}\times w_{i-1}}\) is a weighting matrix, \(\mathbf{b}^{(i)}\in\mathbb{R}^{p_{i}w_{i}}\) is a bias vector, and \(\mathbf{y}^{(i-1)}\) denotes the output of the previous layer with \(\mathbf{y}^{(0)}:=\mathbf{x}\in\mathbb{R}^{n}\).

Now, various activation functions have been proposed. As already stated in the introduction, we here focus on PWA activation functions, i.e., we consider the ReLU activation function

\[\mathbf{g}^{(i)}_{\text{ReLU}}(\mathbf{z}^{(i)})=\max\left\{\mathbf{0},\mathbf{z}^{(i)}\right\}: =\begin{pmatrix}\max\left\{0,\mathbf{z}^{(i)}_{1}\right\}\\ \vdots\\ \max\left\{0,\mathbf{z}^{(i)}_{w_{i}}\right\}\end{pmatrix} \tag{8}\]

and the maxout activation function

\[\mathbf{g}^{(i)}_{\text{max}}(\mathbf{z}^{(i)})=\begin{pmatrix}\max\limits_{1\leq j\leq p _{i}}\left\{\mathbf{z}^{(i)}_{j}\right\}\\ \vdots\\ \max\limits_{p_{i}:(w_{i}-1)+1\leq j\leq p_{i}w_{i}}\left\{\mathbf{z}^{(i)}_{j} \right\}\end{pmatrix}, \tag{9}\]

where we use the shorthand notation

\[\max\limits_{1\leq j\leq p_{i}}\left\{\mathbf{z}^{(i)}_{j}\right\}:=\max\left\{\mathbf{z}^ {(i)}_{1},\ldots,\mathbf{z}^{(i)}_{p_{i}}\right\}.\]
We will refer to the resulting NN as ReLU NN and maxout NN, respectively. The proof of (Goodfellow et al., 2013, Thm. 4.3) shows that maxout NN are PWA functions of the form (1) with \(\mathbf{G}^{(i)}=\mathbf{K}^{(i)}_{\text{NN}}\), \(\mathbf{g}^{(i)}=\mathbf{b}^{(i)}_{\text{NN}}\), \(r=r_{\text{NN}}\), polyhedral sets \(\mathcal{P}^{(i)}=\mathcal{R}^{(i)}_{\text{NN}}\;\forall i\in\{1,\ldots,r_{ \text{NN}}\}\) and local gain \(\mathbf{K}_{\text{NN}}(\mathbf{x})\). The number of parameters needed to describe a NN is

\[\#_{p}:=\sum_{i=1}^{\ell}(w_{i-1}+1)p_{i}w_{i}+(w_{\ell}+1)w_{\ell+1},\]

where \(\ell\), \(p_{i}\) with \(i\in\{1,\ldots,\ell\}\) and \(w_{i}\) with \(i\in\{1,\ldots,\ell+1\}\) describe the topology of the NN.

### Approximate MPC

The use of approximate MPC is particularly useful, if the OCP (3) is too complex to be solved online and the explicit solution has too many regions and thus a large memory footprint ((Kvasnica and Fikar, 2012)). In this case, the exact control law may be approximated by a function that is fast to evaluate and has a small memory footprint. A promising candidate for such a function is an NN, as it combines both required properties and is a universal function approximator (Hornik et al., 1989, Thm. 2.4). In addition, due to the common PWA structure of some NN (Hanin, 2017, Thm. 2), (Arora et al., 2016, Thm. 2.1) and the control law (Bemporad et al., 2002, Thm. 4), they seem to be a natural choice. In fact, for a suitable choice of the weighting matrices and bias vectors, ReLU NN can represent the control law exactly (Schulze Darup, 2020, Thm. 1), (Karg and Lucia, 2020, Thm. 1). Unfortunately, despite the ability of NN to exactly represent the control law, the approximated version of the control law typically does not preserve desirable properties of the MPC, such as stability and performance. Crucial for preserving these properties is the error function

\[\mathbf{e}(\mathbf{x}):=\mathbf{\pi}(\mathbf{x})-\mathbf{\Phi}(\mathbf{x}). \tag{10}\]

More precisely we have to compute the maximum error

\[\overline{e}_{\alpha}:=\max_{x\in\mathcal{X}}\left\|\mathbf{e}(\mathbf{x})\right\|_{\alpha} \tag{11}\]

and the \(\alpha\)-Lipschitz constant of the error

\[\mathcal{L}_{\alpha}(\mathbf{e},\mathcal{X}):=\sup_{\mathbf{x}\neq\mathbf{y}\in\mathcal{X }}\frac{\left\|\mathbf{e}(\mathbf{x})-\mathbf{e}(\mathbf{y})\right\|_{\alpha}}{\left\|\mathbf{x}- \mathbf{y}\right\|_{\alpha}}.\]

According to (Gorokhovic et al., 1994, Prop. 3.4), the \(\alpha\)-Lipschitz constant of a PWA function is equal to the maximum \(\alpha\)-norm of the local gain. If we consider a PWA NN, the error is as difference of two PWA functions also PWA (Gorokhovic et al., 1994, Prop. 1.1) and the \(\alpha\)-Lipschitz constant is thus

\[\mathcal{L}_{\alpha}(\mathbf{e},\mathcal{X})=\max_{\mathbf{x}\in\mathcal{X}}\left\|\bm {K}_{\text{MPC}}(\mathbf{x})-\mathbf{K}_{\text{NN}}(\mathbf{x})\right\|_{\alpha}. \tag{12}\]

Now, if \(e_{\alpha}\) and \(\mathcal{L}_{\alpha}(\mathbf{e},\mathcal{T})\) are below certain values, specified in (Fabiani and Goulart, 2022, Eq. (24)-(25)), then the closed-loop system with the approximated NN controller converges exponentially to the origin according to (Fabiani and Goulart, 2022, Thm. 3.4). The problem at this point is that during the training of a NN, the error (10) is only evaluated at discrete samples \((\mathbf{x}_{i}^{\top}\;\mathbf{\pi}(\mathbf{x}_{i})^{\top})\) of the control law, and the parameters of the NN are chosen such that the mean squared error (MSE)

\[\hat{e}^{2}:=\frac{1}{D}\sum_{i=1}^{D}\left\|\mathbf{\pi}(\mathbf{x}_{i})-\mathbf{\Phi}( \mathbf{x}_{i})\right\|_{2}^{2} \tag{13}\]

over the \(D\in\mathbb{N}\) training samples is minimized. We thus do not have any guarantees for the error at points \(\mathbf{x}\in\mathcal{F}_{N}\) not included in the training samples. Moreover, a low MSE does not necessarily mean that the values of (11) and (12) are low. Therefore, to certify stability of the closed-loop system with a pre-trained NN we need a way to compute these values exactly. For \(\alpha=\{1,\infty\}\) this is possible by solving a mixed-integer linear program (MILP) (Fabiani and Goulart, 2022, Thm. 6.1) if both the output \(\mathbf{\Phi}(\mathbf{x})\) and the local gain \(\mathbf{K}_{\text{NN}}(\mathbf{x})\) of the NN can be computed by solving an MI feasibility problem. Which is proven for ReLU NN in (Fabiani and Goulart, 2022, Thm. 6.1). In the remainder of this paper, we will extend the results to maxout NN.

## 3 Maxout neural networks for approximate MPC

In most of the recent work where the control law (5) is approximated by an NN, ReLU NN are used as function approximators (see, e.g., (Chen et al., 2022; Drummond et al., 2022; Karg and Lucia, 2020)). Maxout NN are rarely considered in this context, although they offer a number of advantages. First, the fact that ReLU NN can represent every PWA function exactly is often used as justification for their use to approximate the PWA control law. However, most ReLU NN that allow an exact description are based on the representation of PWA functions as a sum of the type

\[\hat{F}(\mathbf{x})=\sum_{i=1}^{M}\sigma_{i}\max_{1\leq j\leq J}\{\mathbf{\beta}^{(i)}_ {j}\mathbf{x}+\gamma^{(i)}_{j}\}. \tag{14}\]

From (Wang and Sun, 2005) and (Kripfganz and Schulze, 1987) it is known for which \(M\) and \(J\) we can find parameters \(\sigma_{i}\), \(\mathbf{\beta}_{j}\) and \(\gamma_{j}\) such that \(\hat{F}(\mathbf{x})=F(\mathbf{x})\) holds for \(m=1\). Since (14) is a maxout NN with \(\ell=1\), \(w_{1}=M\), \(p_{1}=J\), \(w_{2}=1\), these results can be used directly to find suitable topologies for maxout NN that allow an exact description of the control law. For ReLU NN these results are not directly applicable. Therefore, (14) is decomposed in, e.g., (Hanin, 2017; Arora et al., 2016) to find a ReLU topology that allows an exact description. Such a topology is used in, e.g., (Karg and Lucia, 2020, Thm. 1) to represent the control law. The decomposition step typically leads to a more conservative ReLU topology in terms of number of layers \(\ell\) and neurons per layer \(w_{i}\) compared to a maxout NN that directly represents (14). Another advantage of the maxout activation is that it trivially include the ReLU activation as a special case. In fact, a ReLU activation is a maxout activation with \(p_{i}=2\) where every second affine segment is set to zero (cf. (8) and
(9)). Thus an approach that allow the computation of (11) and (12) for maxout NN is also applicable to ReLU NN and hence extends the known results.

### Exact maxout neural networks

**Corollary 1**.: _Let \(F(\mathbf{x})\) be an arbitrary PWA function of the form (1) with one dimensional output, i.e, \(m=1\). Then for a maxout NN \(\Phi(\mathbf{x})\) with \(\ell=1\), \(w_{2}=1\) and \(\mathbf{b}^{(2)}=0\) there exist parameters_

1. \(p_{1}\!\in\!\mathbb{N}\)_,_ \(\mathbf{W}^{(1)}\!\in\!\mathbb{R}^{2p_{1}\times n}\)_,_ \(\mathbf{b}^{(1)}\!\in\!\mathbb{R}^{2p_{1}\times 1}\) _and_ \(\mathbf{W}^{(2)}\!\in\!\mathbb{R}^{1\times 2}\)__
2. \(w_{1}\!\in\!\mathbb{N}\)_,_ \(\mathbf{W}^{(1)}\!\in\!\mathbb{R}^{w_{1}(n+1)\times n}\)_,_ \(\mathbf{b}^{(1)}\!\in\!\mathbb{R}^{w_{1}(n+1)\times 1}\) _and_ \(\mathbf{W}^{(2)}\!\in\!\mathbb{R}^{1\times w_{1}}\)__

_with \(\Phi(\mathbf{x})=F(\mathbf{x})\)._

Proof.: Since the resulting maxout NN are of the form (14) with \(J=p_{1}\), \(M=2\) for (i) and \(J=n+1\), \(M=w_{1}\) for (ii), the proof follows from (Kripfganz and Schulze, 1987, Lem. 1) and (Wang and Sun, 2005, Thm. 1), respectively. 

Corollary 1 provides two different topologies (i) and (ii) for maxout NN with one hidden layer that can represent every PWA function with one dimensional output exactly, if \(p_{1}\) and \(w_{1}\), respectively are chosen large enough. Although the results are formulated for \(m=1\) they can also be applied to PWA functions with \(m>1\) by applying Corollary 1 to every dimension of the output individually. Thus maxout NN can represent the PWA control law (5) for arbitrary state and input dimension.

### Maxout neural networks as MILP

In this section, we will derive our main result, which is the computation of the output and the local gain of NN with maxout activation by solving an MI feasibility problem. This result then allows to compute (11) and (12) for the case where the control law (5) is approximated by a maxout NN. Both values are according to the descriptions in Section 2.3 sufficient to prove stability of the closed-loop system. Moreover, these values provide a more profound way to evaluate the success of the training than by just considering the error at the training samples (13) as in, e.g., (Karg and Lucia, 2020; Teichrib and Schulze Darup, 2021). We will derive our results based on the observation that the output of a maxout NN with \(\ell\) hidden layers can be modeled by the recursion

\[\mathbf{y}^{(0)} =\mathbf{x},\] \[\mathbf{y}^{(i)} =\Delta^{(i)}(\mathbf{W}^{(i)}\mathbf{y}^{(i-1)}+\mathbf{b}^{(i)}),\;1\leq i \leq\ell,\] \[\mathbf{\Phi}(\mathbf{x}) =\mathbf{W}^{(\ell+1)}\mathbf{y}^{(\ell)}+\mathbf{b}^{(\ell+1)} \tag{15}\]

where \(\Delta^{(i)}\) is a block diagonal matrix of the form

\[\Delta^{(i)}:=\operatorname{diag}(\mathbf{\delta}^{(i)}_{1:p_{i}},\ldots,\mathbf{ \delta}^{(i)}_{p_{i}(w_{i}-1)+1:p_{i}:w_{i}}) \tag{16}\]

with binary variables \(\mathbf{\delta}^{(i)}\in\mathbb{R}^{1\times p_{i}w_{i}}\). The matrix \(\Delta^{(i)}\) is such that for all elements \(\mathbf{\delta}^{(i)}_{j}\) the logical implication

\[[\mathbf{\delta}^{(i)}_{k_{s}}=1]\Longleftrightarrow[\mathbf{W}^{(i)}_{j} \mathbf{y}^{(i-1)}+\mathbf{b}^{(i)}_{j}\leq\mathbf{W}^{(i)}_{k_{s}}\mathbf{y}^{(i-1)}+\mathbf{b}^{ (i)}_{k_{s}},\] \[\qquad\qquad\qquad\forall j\in\mathcal{A}^{(i)}_{s}\setminus k_ {s},\;k_{s}\in\mathcal{A}^{(i)}_{s}],\] \[\forall s\in\{1,\ldots,w_{i}\},\;\forall i\in\{1,\ldots,\ell\}. \tag{17}\]

holds. Thus, for every neuron in layer \(i\) the matrix \(\Delta^{(i)}\) selects the largest affine segment among all \(p_{i}\) segments. With the results from (Fischetti and Jo, 2018, Sec. 2) we can model the logical implication (17) by the following MI linear constraints

\[\mathbf{q}^{(i)}_{s} \leq\mathbf{W}^{(i)}_{j}\mathbf{q}^{(i-1)}+\mathbf{b}^{(i)}_{j}+\overline{ \overline{b}}^{(i)}(1-\mathbf{\delta}^{(i)}_{j}),\] \[-\mathbf{q}^{(i)}_{s} \leq-\mathbf{W}^{(i)}_{j}\mathbf{q}^{(i-1)}-\mathbf{b}^{(i)}_{j}-\varepsilon( 1-\mathbf{\delta}^{(i)}_{j}),\] \[\mathbf{q}^{(0)} =\mathbf{x},\] \[\sum_{j\in\mathcal{A}^{(i)}_{s}}\mathbf{\delta}^{(i)}_{j}=1,\] \[\qquad\forall j\in\mathcal{A}^{(i)}_{s},\;\forall s\in\{1, \ldots,w_{i}\},\;\forall i\in\{1,\ldots,\ell\}, \tag{18}\]

with a constant upper bound \(\overline{\overline{b}}^{(i)}\in\mathbb{R}\) and a small \(\varepsilon\geq 0\). The variables \(\mathbf{q}^{(i)}\in\mathbb{R}^{w_{i}}\) and \(\mathbf{\delta}^{(i)}\) are real and binary optimization variables, respectively.

**Lemma 2**.: _Let \(\mathbf{q}^{(i)}\) and \(\mathbf{\delta}^{(i)}\) be such that the constraints (18) with \(\varepsilon=0\) hold. Then, the output of the maxout NN (6) is given by_

\[\mathbf{\Phi}(\mathbf{x})=\mathbf{W}^{(\ell+1)}\mathbf{q}^{(\ell)}+\mathbf{b}^{(\ell+1)}. \tag{19}\]

Proof.: If we can show that

\[\mathbf{q}^{(i)}=\mathbf{y}^{(i)} \tag{20}\]

holds for all \(i\in\{0,\ldots,\ell\}\) then (19) holds according to (6) and (7). We prove this by induction. The base case \(i=0\) is true by assumption since we have \(\mathbf{q}^{(0)}=\mathbf{x}=\mathbf{y}^{(0)}\). Moreover, the constraints (18) are such that for every index set \(\mathcal{A}^{(i)}_{s}\) there exists exactly one \(k_{s}\in\mathcal{A}^{(i)}_{s}\) with \(\mathbf{\delta}^{(i)}_{k_{s}}=1\) and \(\mathbf{\delta}^{(i)}_{j}=0,\forall j\in\mathcal{A}^{(i)}_{s}\setminus k_{s}\). Thus if we assume that the induction hypothesis (20) is true for one \(i=t\) we obtain for \(i=t+1\) the constraints

\[\mathbf{q}^{(t+1)}_{s} =\mathbf{W}^{(t+1)}_{k_{s}}\mathbf{y}^{(t)}+\mathbf{b}^{(t+1)}_{k_{s}},\;k_{s} \in\mathcal{A}^{(t+1)}_{s},\] \[\mathbf{q}^{(t+1)}_{s} \leq\mathbf{W}^{(t+1)}_{j}\mathbf{y}^{(t)}+\mathbf{b}^{(t+1)}_{j}+\overline{ \overline{b}}^{(t+1)},\;\forall j\in\mathcal{A}^{(t+1)}_{s}\setminus k_{s},\] \[\mathbf{q}^{(t+1)}_{s} \geq\mathbf{W}^{(t+1)}_{j}\mathbf{y}^{(t)}+\mathbf{b}^{(t+1)}_{j},\;\forall j \in\mathcal{A}^{(t+1)}_{s}\setminus k_{s},\] \[\forall s \in\{1,\ldots,w_{t+1}\}.\]

The first two constraints always hold for a large \(\overline{b}^{(i)}\). Hence the constraints (18) imply

\[\mathbf{q}^{(t+1)}_{s}=\mathbf{W}^{(t+1)}_{k_{s}}\mathbf{y}^{(t)}+\mathbf{b}^{(t+ 1)}_{k_{s}}\geq\mathbf{W}^{(t+1)}_{j}\mathbf{y}^{(t)}+\mathbf{b}^{(t+1)}_{j},\] \[\forall j\in\mathcal{A}^{(t+1)}_{s}\setminus k_{s},\;\forall s\in\{1, \ldots,w_{t+1}\}.\]

Which is exactly the relation (17) for \(i=t+1\), i.e., the largest affine segment of the \(s\)-th neuron is the \(k_{s}\)-th segment. Thus the
relation

\[\left(\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
Thus the solution exists and is unique for \(\mathbf{x}\in\mathbb{R}^{n}\setminus\overline{\mathbf{B}}\). Next, we prove that (24) holds for feasible \(\mathbf{x}\), starting with the case \(\mathbf{\delta}_{h}^{(i)}=0\), where we have

\[\tilde{\mathbf{\xi}}_{h,r}^{(i)}=0\text{ and }\underline{w}_{h,r}^{(i)}\leq\tilde{\mathbf{W} }_{h,r}^{(i)}\leq\overline{w}_{h,r}^{(i)}\]

and for \(\mathbf{\delta}_{h}^{(i)}=1\)

\[\underline{w}_{h,r}^{(i)}\leq\tilde{\mathbf{\xi}}_{h,r}^{(i)}\leq \overline{w}_{h,r}^{(i)}\text{ and }\tilde{\mathbf{\xi}}_{h,r}^{(i)}=\tilde{\mathbf{W}}_{h,r}^{(i)},\] \[\forall h\in\{1,\ldots,p_{i}w_{i}\},\;\forall r\in\{1,\ldots,n\}, \;\forall i\in\{1,\ldots,\ell\}.\]

We can now rewrite the equality constraint in the third line of (23) as follows

\[\mathbf{\xi}_{1,r}^{(i)}=\sum_{h=1}^{p_{i}}\mathbf{\delta}_{h}^{(i)} \tilde{\mathbf{W}}_{h,r}^{(i)}=\Delta_{1,:}^{(i)}\tilde{\mathbf{W}}_{:,r}^{(i)},\] \[\vdots\] \[\mathbf{\xi}_{w_{i},r}^{(i)}=\sum_{h=p_{i}(w_{i}-1)+1}^{p_{i}w_{i}} \mathbf{\delta}_{h}^{(i)}\tilde{\mathbf{W}}_{h,r}^{(i)}=\Delta_{w_{i},:}^{(i)}\tilde{ \mathbf{W}}_{:,r}^{(i)},\] \[\forall r\in\{1,\ldots,n\},\;\forall i\in\{1,\ldots,\ell\}.\]

This can be written in a more compact form as matrix multiplication

\[\mathbf{\xi}^{(i)}=\Delta^{(i)}\tilde{\mathbf{W}}^{(i)}=\Delta^{(i)}\mathbf{W}^{(i)}\mathbf{ \xi}^{(i-1)},\;\forall i\in\{1,\ldots,\ell\}.\]

Starting with \(\mathbf{\xi}^{(0)}=\mathbf{I}\), this recursion leads to

\[\mathbf{\xi}^{(\ell)}=\prod_{i=1}^{\ell}\Delta^{(i)}\mathbf{W}^{(i)}.\]

By substituting the former relation in (21) we can show that the local gain is indeed given by (24), which completes the proof. 

The constraints (18) and (23) can now be used to model a maxout NN by MI linear constraints and to compute the output (19) and the local gain (24) by solving an MI feasibility problem. If we further combine the results of Lemma 2 and 3, we can compute the maximum error and the \(\alpha\)-Lipschitz constant of the error.

**Theorem 4**.: _Let \(\mathbf{\Phi}(\mathbf{x})\) be a maxout NN and \(\alpha=\{1,\infty\}\). Then the maximum error (11) and the \(\alpha\)-Lipschitz constant of the error (12) can be computed by solving an MILP._

Proof.: The proof follows from (Fabiani and Goulart, 2022, Thm. 6.1), where it is shown that both values can be computed by solving an MILP, if the output \(\mathbf{\Phi}(\mathbf{x})\) and the local gain \(\mathbf{K}_{\text{NN}}(\mathbf{x})\) of the NN can be computed by solving an MI feasibility problem. This is always possible for maxout NN according to Lemma 2 and 3. Thus, by replacing the MI linear constraints for the ReLU NN in (Fabiani and Goulart, 2022, Thm. 6.1) with the constraints (18) and (23), we can compute (11) and (12) for a maxout NN by solving an MILP. 

## 4 Numerical examples

We consider two simple examples to highlight our theoretical findings. In both examples, the maximum error (11) and the \(\alpha\)-Lipschitz constant of the error (12) are computed by solving the MILP according to Theorem 4 with the MOSEK optimization toolbox for MATLAB (see (MOSEK ApS, 2022)). The NN are implemented in Python with Keras (Chollet et al., 2015) and Tensorflow (Abadi et al., 2015). During the training, the MSE (13) is minimized with respect to the weighting matrices and bias vectors of the NN using a stochastic gradient descent.

### Example system with \(n=1\)

We consider the system from (Schulze Darup and Cannon, 2016, Ex. 2) with the dynamics

\[x(k+1)=\tfrac{6}{5}x(k)+u(k)\]

and the constraints \(\mathcal{X}=[-10,10]\) and \(\mathcal{U}=[-1,1]\). As in (Schulze Darup and Cannon, 2016), we choose \(Q=19/5\), \(R=1\), \(P=5\), and \(\mathcal{T}=[-1,1]\). Finally, we select \(N=2\) for illustration purposes here. Explicitly solving the OCP (3) then leads to the control law

\[\pi(x)=\left\{\begin{array}{ccc}-1&\text{if}&x\in\left[-\tfrac{20}{9},-1 \right],\\ -x&\text{if}&x\in\left[-1,1\right],\\ 1&\text{if}&x\in\left[1,\tfrac{20}{9}\right].\end{array}\right. \tag{27}\]

For a maxout NN with the topology from Corollary 1 (i) with \(p_{1}=2\),

\[\mathbf{W}^{(1)} =\begin{pmatrix}-1&0&-1&0\end{pmatrix}^{\top},\,\mathbf{b}^{(1)}= \begin{pmatrix}0&-1&-1&0\end{pmatrix}^{\top}\text{ and }\] \[\mathbf{W}^{(2)} =\begin{pmatrix}1&-1\end{pmatrix} \tag{28}\]

we have

\[\Phi(x)=\max\{-x,-1\}-\max\{-x-1,0\}=\pi(x). \tag{29}\]

We implemented the constraints (18) and (23) with \(\overline{b}^{(i)}=\overline{w}_{h,r}^{(i)}=-\underline{w}_{h,r}^{(i)}=10^{4}\), \(\varepsilon=10^{-5}\) and used them to compute (11) and (12) with \(\alpha=\infty\) according to Theorem 4. For a maxout NN with the parameters (28), which exactly represents the control law, we computed \(\overline{e}_{\infty}=4.5\times 10^{-16}\) and \(\mathcal{L}_{\infty}(e,\mathcal{T})=0\). This can be seen as an experimental evidence that the implementation computes correctly the maximum error and the \(\infty\)-Lipschitz constant of the error.

To obtain the data in Table 1, we sampled randomly \(1000\) points from the control law (27) and trained different maxout NN for \(1000\) epochs.

The results show that a maxout NN with \(w_{1}=2\) and \(p_{1}=2\), which can represent the control law exactly (cf. (29)), is sufficient to get a good approximation, in terms of maximum error and \(\infty\)-Lipschitz constant of the error. We only get an approximation with a relatively high error for the first and last topology of Table 1. This is not surprising since for \(w_{1}=1\), \(p_{1}=4\) the maxout NN is a convex PWA function and for \(w_{1}=4\), \(p_{1}=1\) it is an affine function. In both cases, it is not possible to find an approximation of (27) with a maximum error close to zero.
### Example system with \(n=2\)

We consider the OCP (3) with

\[\mathbf{A} =\begin{pmatrix}1&1\\ 0&1\end{pmatrix},\quad\mathbf{B}=\begin{pmatrix}0.5\\ 1\end{pmatrix},\] \[\mathcal{X} =\{\mathbf{x}\in\mathbb{R}^{2}\mid|\mathbf{x}_{1}|\leq 25,\,|\mathbf{x}_{2}| \leq 5\},\] \[\mathcal{U} =\{u\in\mathbb{R}\mid|u|\leq 1\},\]

and choose \(\mathbf{Q}=\mathbf{I}\), \(R=1\), \(N=3\), \(\mathbf{P}\) as the solution of the discrete-time algebraic Riccati equation and \(\mathcal{T}\) as the maximal output admissible set (see (Gilbert and Tan, 1991) for details). Explicitly solving the OCP leads to a control law with \(29\) regions. Using the procedure described in (Kripfganz and Schulze, 1987, Sec. 1), we find a maxout NN of the type (i) from Corollary 1 with \(p_{1}=38\) and \(231\) parameters, that exactly describes the control law. For this NN we computed \(\overline{e}_{\infty}=2.68\times 10^{-12}\) and \(\mathcal{L}_{\infty}(e,\mathcal{T})=3.11\times 10^{-5}\). A ReLU NN that can exactly represent the control law has according to (Karg and Lucia, 2020, Thm. 1) \(457\) parameters. This shows that, as already stated in Section 3, ReLU NN for the exact representation of PWA functions are more conservative compared to maxout NN.

Table 2 summarizes the training results for different maxout NN trained for \(1000\) epochs on \(10^{4}\) samples of the control law. The values of (11) and (12) are computed as in the first example, except that in this example \(\varepsilon=10^{-3}\) is chosen.

The experimental data indicate that the first maxout NN is the best choice. Note that this is a maxout NN with the topology (i), which theoretically allows an exact representation of the control law, trained on samples of the control law. The seventh maxout NN is of the type (ii), where \(w_{1}\) is chosen such that the number of parameters is equal to that of the first maxout NN. It has a comparable maximum error but a \(\infty\)-Lipschitz constant that is about \(183\)-times higher. This observation may be explained by the fact that maxout NN with more neurons represent PWA functions with more regions compared to maxout NN with less neurons (see (Montufar et al., 2021, Thm. 3.6)). This makes it more likely that there exits one region for which the deviation between the local gains of the NN and the control law is larger, resulting in a larger \(\infty\)-Lipschitz constant (cf. (12)). The increasing \(\infty\)-Lipschitz constant with the number of neurons \(w_{1}\) in Table 2 supports this assumption.

## 5 Conclusions and outlook

We presented a method to compute two values, i.e., the maximum error (11) with respect to the optimal control law and the Lipschitz constant (12) of the error function (10), which are crucial to certify stability of the closed-loop system with a maxout NN controller that approximates the optimal control law of MPC. Our results are derived by showing how the output and the local gain of maxout NN can be computed by solving an MI feasibility problem (cf. Lems. 2 and 3). The combination of both lemmas leads to Theorem 4, which states that the maximum error and the \(\alpha\)-Lipschitz constant of the error function can be computed by solving an MILP. The computation of both values has been successfully applied to a number of realizations, including maxout NN which exactly describe the control law.

An interesting direction for future research is the combination of the proposed method with the results from (Teichrib and Schulze Darup, 2022), where maxout NN are designed that allow an exact description of the piecewise quadratic optimal value function in MPC. Such a combination may provide a new method to analyze maxout NN approximations of the optimal value function.
