# Error-free approximation of explicit linear MPC through lattice piecewise
affine expression

Jun Xu<EMAIL_ADDRESS>Yunjiang Lou<EMAIL_ADDRESS>Bart De Schutter
<EMAIL_ADDRESS>Zhenhua Xiong<EMAIL_ADDRESS>School of
Mechanical Engineering and Automation, Harbin Institute of Technology,
Shenzhen  Key Laboratory of System Control and Information Processing,
Ministry of Education, Shanghai  Delft Center for Systems and Control, Delft
University of Technology, The Netherlands  State Key Laboratory of Mechanical
Systems and Vibration, Shanghai Jiao Tong University, Shanghai, China

###### Abstract

In this paper, the disjunctive and conjunctive lattice piecewise affine (PWA)
approximations of explicit linear model predictive control (MPC) are proposed.
The training data are generated uniformly in the domain of interest,
consisting of the state samples and corresponding affine control laws, based
on which the lattice PWA approximations are constructed. Re-sampling of data
is also proposed to guarantee that the lattice PWA approximations are
identical to the explicit MPC control law in the unique order (UO) regions
containing the sample points as interior points. Additionally, if all the
distinct affine functions have been sampled, the disjunctive lattice PWA
approximation constitutes a lower bound while the conjunctive lattice PWA
approximation formulates an upper bound of the original optimal control law.
The equivalence of the two lattice PWA approximations then guarantees that the
approximations are error-free in the domain of interest, which is tested
through a statistical guarantee. The complexity of the entire procedure is
analyzed, which is polynomial with respect to the number of samples. The
performance of the proposed approximation strategy is tested through two
simulation examples, and the results show that with a moderate number of
sample points, we can construct lattice PWA approximations that are equivalent
to the optimal control law of the explicit linear MPC.

###### keywords:

linear MPC; lattice piecewise affine; error-free approximation.

AND red

††thanks: This work was supported in part by the National Natural Science
Foundation of China under Grant U1813224, 62173113, and in part by the Science
and Technology Innovation Committee of Shenzhen Municipality under Grant
GXWD20201230155427003-20200821173613001, JCYJ20200109113412326. The material
in this paper was partially presented at the 60th IEEE Conference on Decision
and Control, December 13-17, 2021, Austin, Texas, USA.

, , ,

##  1 Introduction

Model predictive control (MPC) is currently the most popular control
methodology employed in process control and its impact on industry has been
recognized widely [1]. In MPC, the control action is calculated through
solving a finite-horizon open-loop optimal control problem at each sampling
instant, which is computationally expensive and only suitable for systems with
slow dynamics. Therefore, for fast dynamical systems, in order to use MPC, the
complexity of the online optimization should be reduced. A natural thought is
to move the online optimization offline, which is the idea of explicit MPC.
Explicit MPC was first introduced in [2], in which the linear MPC problem is
formulated as a multi-parametric quadratic programming (mpQP) problem and
solved offline. The optimal control law is proved to be continuous piecewise
affine (PWA) with respect to the state. The subregions as well as the
corresponding affine functions defined on them are recorded. For online
implementation, given the current state, one must only find the subregion in
which the state lies, and the function evaluation of the corresponding affine
function gives rise to the optimal control law.

However, the offline construction of such subregions, the memory required to
store them, and the online search for the right subregion are the main
limitations of explicit MPC [3]. Much work has been done to solve these three
problems. To overcome the complex offline geometric computations,
combinatorial approaches are proposed that are based on implicitly enumerating
all possible combinations of active constraints [4, 5]. To reduce the memory
required to store the subregions as well as the control laws, region-free
explicit MPC is proposed [3, 6]. Moreover, the online search complexity can be
reduced by storing additional information [7, 8], introducing an improved
binary search tree (orthogonal truncated binary search tree) [9], or resorting
to the method of convex lifting [10, 11]. The lattice PWA representation has
also been used to exactly express the explicit MPC control law, resulting in a
much lower storage requirement [12, 13]. As the complexity of solving the
explicit MPC problem increases exponentially with the size of the optimization
problem, all these methods can only alleviate the computational burden to some
extent.

Another idea is to formulate the approximate MPC controller [14, 15] or semi-
explicit MPC controller [16]. In these methods, training data containing the
values of states and corresponding optimal control laws of the MPC problem are
generated, and the approximated controller is constructed using these data. In
general, the samples are required to be distributed sufficiently evenly over
the domain [17]. Different approaches have been used to generate the
approximation, such as the canonical piecewise linear function [18], radial
basis functions [19], wavelets [20], and so on. In addition, reinforcement
learning has also been used to derive a data-driven MPC control law in [21].
In the work of [15], [20], and [22], the approximations are based on
particular partitions of the domain of interest, and the interpolation based
algorithm can be developed [23]. In fact, the partitions of the domain of
interest employed in these works are different from the domain partitions in
the explicit linear MPC law, e.g., the continuous PWA solution of the linear
MPC problem has not been fully explored in the approximation.

To resemble explicit MPC control law to a larger extent, the information of
the local affine functions of the linear explicit MPC is utilized to derive
the lattice PWA approximation. In our previous work, the lattice PWA
_representation_ of the explicit MPC control law was derived, which however
also scales poorly with the dimension of the parameters. In this paper, in
order to handle more complex problems, we present an _error-free
approximation_ that consists of disjunctive and conjunctive lattice PWA
approximations, which coincide with the explicit MPC control in a statistical
sense. Moreover, the offline calculation complexity depends mainly on the
number of sample points, which scales well with the state dimension. A
preliminary thought of the disjunctive lattice PWA approximation of explicit
linear MPC was presented in [24], in which the approximated control law is not
guaranteed to be error-free. However, in this work, under mild assumptions,
the equivalence of the disjunctive and conjunctive approximations guarantees
that the two approximations are identical to the optimal control law in the
domain of interest. The approximation can also be simplified to further lower
the storage requirements and online computational complexity.

The rest of this paper is organized as follows. Section 2 gives the
preliminaries about the explicit linear MPC problem and the lattice PWA
expression. The offline approximations of the explicit linear MPC control law
through the lattice PWA expression are given in detail in Section 3, in which
the sampling and re-sampling procedures, as well as the simplification of the
approximation are provided. In Section 4, the approximation error and the
complexity of the proposed procedure are analyzed. Section 5 provides the
simulation results and the paper ends with conclusions and plans for future
work in Section 6.

##  2 Preliminaries

###  2.1 Explicit linear MPC problem

In particular, MPC for a discrete-time linear time-invariant system can be
cast as the following optimization problem at time step $t$:

$\displaystyle\min\limits_{U}$
$\displaystyle\Bigg{\\{}J(U,\bm{x}_{0})=v_{N_{p}}(\bm{x}_{N_{p}})+\sum\limits_{k=0}^{N_{p}-1}v(\bm{x}_{k},\bm{u}_{k})\Bigg{\\}}$
(1a) s.t.
$\displaystyle\bm{x}_{k+1}=A\bm{x}_{k}+B\bm{u}_{k},k=0,\ldots,N_{p}-1$ (1b)
$\displaystyle(\bm{x}_{k},\bm{u}_{k})\in\mathcal{G},k=0,\ldots,N_{p}-1$ (1c)
$\displaystyle\bm{x}_{N_{p}}\in\mathcal{F}$ (1d)

in which the optimized variable is
$U=[\bm{u}_{0}^{T},\ldots,\bm{u}_{N_{p}-1}^{T}]^{T}$, $N_{p}$ is the
prediction horizon, the variables $\bm{x}_{k}\in\mathbb{R}^{n_{x}}$ and
$\bm{u}_{k}\in\mathbb{R}^{n_{u}}$ denote the predicted state and input at time
step $k$, respectively, using (1b). The terminal penalty is denoted as $v_{N}$
and $v(\cdot,\cdot)$ is the stage cost; $\mathcal{G}$ and $\mathcal{F}$ are
full-dimensional polyhedral sets of appropriate dimensions. In this paper, we
assume a strictly convex cost, i.e.,
$v_{N}=\bm{x}_{N}^{T}Q_{N}\bm{x}_{N},v(\bm{x}_{k},\bm{u}_{k})=\bm{x}_{k}^{T}Q_{k}\bm{x}_{k}+\bm{u}_{k}^{T}Q_{u}\bm{u}_{k}$,
in which $Q_{u}\succ 0,Q_{k},Q_{N}\succeq 0$. After solving the optimization
problem (1), the optimal
$U^{*}=[({\bm{u}_{0}^{*}})^{T},\ldots,({\bm{u}_{N_{p}-1}^{*})}^{T}]^{T}$ is
obtained, and only $\bm{u}_{0}^{*}$ is applied to the system. The optimization
problem is subsequently reformulated and solved at the next time steps
$t=1,2,\ldots$ by updating the given state vector $\bm{x}_{0}$.

It has been proved in [2] that the solution $U^{*}$ is a _continuous PWA
function_ of the state $\bm{x}_{0}$, and we use $\bm{x}$ instead hereafter in
this paper. In fact, this conclusion is obtained through solving an mpQP
problem of the form

$\begin{array}[]{rl}\min\limits_{U}&\frac{1}{2}U^{T}HU+\bm{x}^{T}FU\\\
s.t.&GU\leq\bm{w}+E\bm{x}\end{array}$ (2)

where $U\in\mathbb{R}^{N_{p}\cdot n_{u}}$ is the vector of optimization
variables, the parameter vector is $\bm{x}\in\mathbb{R}^{n_{x}}$, and the
matrices $H,F,G$, and $E$ are calculated through the optimization problem (1)
[2]. Under the assumption that $Q_{k},Q_{N}\succeq 0,Q_{u}\succ 0$, we have
$H\succ 0$.

The definition of a continuous PWA function as well as the lemma concerning
the continuous PWA property of the solution to the mpQP problem is presented
as follows.

######  Definition 1.

[25] A function $f:\Omega\rightarrow\mathbb{R}^{m}$, where
$\Omega\subseteq\mathbb{R}^{n_{x}}$ is convex, is said to be continuous PWA if
it is continuous on the domain $\Omega$ and the following conditions are
satisfied:

  1. 1.

The domain space $\Omega$ is divided into a finite number of nonempty convex
polyhedra, i.e.,
$\Omega=\cup_{i=1}^{\hat{N}}\Omega_{i},~{}\Omega_{i}\neq\emptyset$, the
polyhedra are closed and have non-overlapping interiors,
$\mathrm{int}(\Omega_{i})\cap\mathrm{int}(\Omega_{j})=\emptyset,~{}\forall
i,j\in\\{1,\ldots,\hat{N}\\},i\neq j$. These polyhedra are also called local
regions. The boundaries of the polyhedra are nonempty sets in
($n-1$)-dimensional space.

  2. 2.

In each local region $\Omega_{i}$, $f$ equals a local affine function
$u_{\mathrm{loc}(i)}$:

$f(\bm{x})=u_{\mathrm{loc}(i)}(\bm{x}),~{}\forall x\in\Omega_{i}.$

######  Lemma 2.

[2] Considering the mpQP problem (2) and assuming that $H\succ 0$, then the
set of feasible parameters $\Omega\subset\mathbb{R}^{n_{x}}$ is convex, the
optimizer $U^{*}:\Omega\rightarrow\mathbb{R}^{p\cdot n_{u}}$ is continuous
PWA, and the value function $J^{*}:\Omega\rightarrow\mathbb{R}$ is continuous
convex and piecewise quadratic.

The details of constructing such continuous PWA function are given as follows.

First, the mpQP problem can be rewritten in the form

$\begin{array}[]{rl}\min\limits_{\bm{z}}&\frac{1}{2}\bm{z}^{T}H\bm{z}\\\
s.t.&G\bm{z}\leq\bm{w}+S\bm{x}\end{array}$ (3)

by letting $\bm{z}=U+H^{-1}F^{T}\bm{x}$ and $S=E+GH^{-1}F^{T}$. Once the
optimal solution $\bm{z}^{*}$ of the optimization problem (3) is available, we
can easily obtain the optimal $U^{*}$ as

$U^{*}=\bm{z}^{*}-H^{-1}F^{T}\bm{x}.$

The optimal solution $\bm{z}^{*}$ for a fixed $\bm{x}$ is fully characterized
by the Karush-Kuhn-Tucker (KKT) conditions:

$\displaystyle
H\bm{z}^{*}+G_{\mathcal{A}^{*}}^{T}\bm{\lambda}^{*}+G_{\mathcal{N}^{*}}^{T}\bm{\mu}^{*}=0$
(4a) $\displaystyle
G_{\mathcal{A}^{*}}\bm{z}^{*}=\bm{w}_{\mathcal{A}^{*}}+S_{\mathcal{A}^{*}}\bm{x}$
(4b) $\displaystyle
G_{\mathcal{N}^{*}}\bm{z}^{*}<\bm{w}_{\mathcal{N}^{*}}+S_{\mathcal{N}^{*}}\bm{x}$
(4c) $\displaystyle\bm{\lambda}^{*}\geq 0$ (4d) $\displaystyle\bm{\mu}^{*}\geq
0$ (4e)
$\displaystyle{\bm{\lambda}^{*}}^{T}(G_{\mathcal{A}^{*}}\bm{z}^{*}-\bm{w}_{\mathcal{A}^{*}}-S_{\mathcal{A}^{*}}\bm{x})=0$
(4f)
$\displaystyle{\bm{\mu}^{*}}^{T}(G_{\mathcal{N}^{*}}\bm{z}^{*}-\bm{w}_{\mathcal{N}^{*}}-S_{\mathcal{N}^{*}}\bm{x})=0$
(4g)

in which (4b) and (4c) are the active and inactive constraints at
$\bm{z}^{*}$, respectively. Assuming that $G\in\mathbb{R}^{p\times N_{p}\cdot
n_{u}},\bm{w}\in\mathbb{R}^{p},S\in\mathbb{R}^{p\times n_{x}}$, and
$G_{i},w_{i}$, and $S_{i}$ denote the $i$-th row of $G,w$, and $S$,
respectively, the active as well as inactive index sets can be written as

$\mathcal{A}^{*}=\\{j\in\\{1,\ldots,p\\}|G_{j}\bm{z}^{*}=\bm{w}_{j}+S_{j}\bm{x}\\}$

and

$\mathcal{N}^{*}=\\{j\in\\{1,\ldots,p\\}|G_{j}\bm{z}^{*}<\bm{w}_{j}+S_{j}\bm{x}\\},$

respectively. It is apparent that
$\mathcal{A}^{*}=\\{1,\ldots,p\\}\setminus\mathcal{N}^{*}$. For an inactive
constraint $j$, we have $\bm{\mu}_{j}^{*}=0$. For a particular
$\mathcal{A}^{*}$, and assume $G_{\mathcal{A}^{*}}$ is full row rank, we have

$\bm{\lambda}^{*}=-(G_{\mathcal{A}^{*}}H^{-1}G_{\mathcal{A}^{*}}^{T})^{-1}(\bm{w}_{\mathcal{A}^{*}}+S_{\mathcal{A}^{*}}\bm{x}),$
(5)

as well as

$\bm{z}^{*}=H^{-1}G_{\mathcal{A}^{*}}^{T}(G_{\mathcal{A}^{*}}H^{-1}G_{\mathcal{A}^{*}}^{T})^{-1}(\bm{w}_{\mathcal{A}^{*}}+S_{\mathcal{A}^{*}}\bm{x}).$
(6)

The local region for which the local affine function (6) is defined is called
_critical region_ , and it can be constructed by the constraints of primal
feasibility (4c) and dual feasibility (4d).

######  Remark 3.

For the case in which $G_{\mathcal{A}^{*}}$ is not full row rank, i.e., the
rows of $G_{\mathcal{A}^{*}}$ are linearly dependent, the linear independence
constraints qualification (LICQ) is violated according to [26], and this is
referred to as primary degeneracy [2] (dual degeneracy cannot occur as $H\succ
0$). Assuming that the rank of $G_{\mathcal{A}^{*}}$ is $r$, we can then
arbitrarily select $r$ independent constraints, and proceed with the new
reduced active index set [27].

To search for all the local affine functions and critical regions, one must
enumerate all possible active index sets $\mathcal{A}^{*}$, apply the KKT
conditions accordingly, and then the continuous PWA control law can be
obtained. In the next subsection, the lattice PWA representation is presented,
which is used to express the resulting continuous PWA control law in our
previous work [13].

###  2.2 Lattice PWA representation

It is stated in [13] that any continuous PWA function can be represented by
the lattice PWA representation.

######  Lemma 4.

[13] Letting $f$ be a continuous PWA function defined in Definition 1, then
$f$ can be represented as

$f(\bm{x})=f_{\mathrm{L,d}}(\bm{x})=\max\limits_{i=1,\ldots,N}\min_{j\in
I_{\geq,i}}u_{j}(\bm{x}),~{}\forall x\in\Gamma,$ (7)

or

$f(\bm{x})=f_{\mathrm{L,c}}(\bm{x})=\min\limits_{i=1,\ldots,N}\max_{j\in
I_{\leq,i}}u_{j}(\bm{x}),~{}\forall x\in\Gamma,$ (8)

in which $I_{\geq,i}=\\{j|u_{j}(\bm{x})\geq
u_{i}(\bm{x}),\forall\bm{x}\in\Gamma_{i}\\},I_{\leq,i}=\\{j|u_{j}(\bm{x})\leq
u_{i}(\bm{x}),\forall\bm{x}\in\Gamma_{i}\\}$, and the expressions $\min_{j\in
I_{\geq,i}}u_{j}(\bm{x})$ and $\max_{j\in I_{\leq,i}}u_{j}(\bm{x})$ are called
terms of $f_{\mathrm{L,d}}$ and $f_{\rm L,c}$, respectively. The affine
function $u_{j}(\bm{x})$ is called a literal. The region $\Gamma_{i}$ is a
unique order (UO) region that is subset of the local region and the order of
the affine functions

$u_{1}(\bm{x}),\ldots,u_{N}(\bm{x}),$ (9)

remains unchanged in the interior of $\Gamma_{i}$. The expressions (7) and (8)
are called full disjunctive and conjunctive lattice PWA representations,
respectively, in which the names “disjunctive” and “conjunctive” come from the
terminology in Boolean algebra.

Considering a two-dimensional continuous PWA function (10) with four affine
pieces, Fig. 1 illustrates the UO region and corresponding PWA function.

######  Example 5.

$f=\left\\{\begin{array}[]{ll}\ell_{1}(\bm{x})=-x_{2}+1&\mathrm{if}~{}\bm{x}\in\Gamma_{1},\\\
\ell_{2}(\bm{x})=-x_{1}+1&\mathrm{if}~{}\bm{x}\in\Gamma_{2},\\\
\ell_{3}(\bm{x})=x_{2}+1&\mathrm{if}~{}\bm{x}\in\Gamma_{3},\\\
\ell_{4}(\bm{x})=x_{1}+1&\mathrm{if}~{}\bm{x}\in\Gamma_{4}.\end{array}\right.$
(10)

The polyhedral regions $\Gamma_{1},\Gamma_{2},\Gamma_{3},\Gamma_{4}$, and the
two-dimensional function $f$ are shown in Fig. 1.

(a) Function.

(b) Regions.

Figure 1: Continuous PWA function in Example 5.

The regions $\Gamma_{1},\Gamma_{2},\Gamma_{3}$, and $\Gamma_{4}$ are local
affine regions, and can be divided into UO regions
$\Gamma_{11},\Gamma_{12},\ldots,\Gamma_{41}$, and $\Gamma_{42}$. Taking the UO
region $\Gamma_{31}$ as an example, the order of affine functions is

$\ell_{3}<\ell_{2}<\ell_{1}<\ell_{4}.$

For this continuous PWA function, as it is concave, both the disjunctive and
conjunctive lattice PWA representations are

$f=\min\\{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\\}.$

According to Lemma 4, we can represent a continuous PWA control law using a
lattice PWA function (either disjunctive or conjunctive). The disjunctive
lattice PWA representation of explicit linear MPC was investigated in [12] and
[13], in which the continuous PWA control law was obtained through the MPT3
toolbox [28] in advance. However, as explained in Section 2.1, for problems
with a large number of constraints and a high-dimensional state, the number of
possible combinations of active constraints increases exponentially and the
derivation of the explicit MPC solution is extremely computationally
expensive. Hence, in this paper, we propose an approximated continuous PWA
control law by sampling only a set of states in the domain of interest. We
show that this approximation utilizes the local affine property of the
original explicit MPC control law and is identical to the original control law
at the sample points and in the UO regions the sample points lie. In addition,
under mild assumptions, the lattice PWA approximations are identical to the
explicit MPC control law in the domain of interest.

##  3 Lattice PWA approximation of explicit linear MPC control law

###  3.1 Generation of sample points in the interior of UO regions

As indicated in Lemma 2, the explicit linear MPC control law $U^{*}$ is a
continuous PWA function with respect to the state $\bm{x}$. Then the first
element of $U^{*}$, which is $\bm{u}_{0}^{*}$, is also a continuous PWA
function of $\bm{x}$, i.e., $\bm{u}_{0}^{*}$ is affine in the local regions
that $\bm{x}$ lies in. For simplicity, we omit the subscript in
$\bm{u}_{0}^{*}$ hereafter in the paper, and use $\bm{u}^{*}$ instead. The
sample points
$(\bm{x}_{i},\bm{u}_{i}(\bm{x}_{i}))\in\mathcal{X}_{1}\times\mathcal{U}_{1}$
are generated in the domain of feasible parameters, in which
$\bm{u}_{i}(\bm{x})$ is the affine function at $\bm{x}_{i}$ such that

$\bm{u}_{i}(\bm{x}_{i})=\bm{u}^{*}(\bm{x}_{i}).$

For simplicity, we consider the case $n_{u}=1$; note however that the proposed
methodology can be easily extended to the case when $n_{u}>1$. Moreover, the
domain of interest is assumed to be a hyperbox.

In this subsection, the training points $\bm{x}_{i}\in\mathcal{X}_{1}$ are
required to be in the interior of UO regions, i.e.,
$\bm{x}_{i}\in\mathrm{int}(\Gamma(\bm{x}_{i}))$, in which $\Gamma(\bm{x}_{i})$
is the corresponding UO region. This means that if
$u_{j_{1}}(\bm{x}_{i})>u_{j_{2}}(\bm{x}_{i})$ holds, then we have

$u_{j_{1}}(\bm{x})>u_{j_{2}}(\bm{x}),\forall\bm{x}\in\Gamma(\bm{x}_{i}).$

As a matter of fact, if there are no affine functions $u_{j_{1}}(\bm{x})$ and
$u_{j_{2}}(\bm{x})$ such that $u_{j_{1}}(\bm{x}_{i})=u_{j_{2}}(\bm{x}_{i})$,
then $\bm{x}_{i}\in\mathrm{int}(\Gamma(\bm{x}_{i}))$. If there is a point
$\bm{x}_{i}$ that is not in the interior of the corresponding UO region, i.e.,
there exist $j_{1},j_{2}$, such that
$u_{j_{1}}(\bm{x}_{i})=u_{j_{2}}(\bm{x}_{i})$, then a perturbation should be
applied, i.e., $\bm{x}_{i}=\bm{x}_{i}+\delta$, such that $\bm{x}_{i}$ is in
the interior of the UO region. This will be shown in Example 12 in Section
3.3. It will be shown in Lemma 7 that the condition that $\bm{x}_{i}$ is in
the interior of UO regions is required for the equivalence of lattice PWA
approximations and the original optimal control law in UO regions.

It should be noted that for the checking of UO regions, only the information
of sampled affine functions is needed, and the global information of the PWA
function is unknown, e.g., whether there is some $u_{j}$ that has not been
sampled or how these sampled affine functions are connected.

Algorithm 1 describes the sampling of training points.

######  Algorithm 1.

Sampling of training points in explicit linear MPC control law.

Input:Linear MPC problem, the number of sample points $N_{1}$, sample domain
$\Omega$. Output:Sample data set $\mathcal{X}_{1}\times\mathcal{U}_{1}$.

1: $\mathcal{X}_{1}=\emptyset$, $\mathcal{U}_{1}=\emptyset$.

2: for $i=1$ to $N_{1}$ do

3: Generate a grid point $\bm{x}_{i}$ in $\Omega$.

4: if $\bm{x}_{i}$ is not an interior point of some UO region then

5: Apply a perturbation until $\bm{x}_{i}$ is an interior point of
$\Gamma(\bm{x}_{i})$.

6: end if

7: Solve corresponding QP problem (3) by letting $\bm{x}=\bm{x}_{i}$.

8: Determine active and inactive index sets $\mathcal{A}^{*}$ and
$\mathcal{N}^{*}$, respectively.

9: Solve KKT conditions (4) to obtain the affine function $\bm{z}_{i}^{*}$
through (6).

10: Obtain optimal input, i.e., $U_{i}(\bm{x}_{i})$ and $u_{i}(\bm{x}_{i})$.

11: end for

For a feasible state $\bm{x}_{i}$, line 7 in Algorithm 1 states that the
optimal $\bm{z}_{i}^{*}=U_{i}^{*}+H^{-1}F^{T}\bm{x}_{i}$ can be obtained
through the QP problem (3), which, together with the information of
$\bm{x}_{i}$, determines the active and inactive constraints (4b) and (4c),
respectively. Therefore, the index set $\mathcal{A}_{i}^{*}$ as well as
$\mathcal{N}_{i}^{*}$ is fixed, and if the matrix $G_{\mathcal{A}_{i}^{*}}$ is
full row rank, the affine function $\bm{z}_{i}(\bm{x}_{i})$ can be calculated
through (6) (the rank-deficient case can be handled as indicated in Remark 3).
Then we have

$U_{i}(\bm{x}_{i})=\bm{z}_{i}^{*}(\bm{x}_{i})-H^{-1}F^{T}\bm{x}_{i},$ (11)

and

$u_{i}(\bm{x}_{i})=\left[\begin{array}[]{cccc}\mathbf{I}_{n_{u}\times
n_{u}}&\bm{0}&\cdots&\bm{0}\end{array}\right]U_{i}(\bm{x}_{i})$ (12)

in which $\mathbf{I}_{n_{u}\times n_{u}}$ is the identity matrix with size
$n_{u}\times n_{u}$.

After evaluating Algorithm 1, we can obtain the sample dataset
$\mathcal{X}_{1}\times\mathcal{U}_{1}$, in which $\mathcal{U}_{1}$ is a set of
affine functions $u_{i}(\bm{x}_{i})$. It is noted that compared with ordinary
sampling, in which only the evaluation of $u_{i}^{*}(\bm{x}_{i})$ is
available, here the corresponding affine function is also recorded, which can
be used for the lattice PWA approximation in Section 3.2.

###  3.2 Lattice PWA approximation based on sample points

We now derive both the disjunctive and conjunctive lattice PWA approximations
based on the sample dateset $\mathcal{X}_{1}\times\mathcal{U}_{1}$.

The disjunctive lattice PWA approximation is constructed via the sample points
and can be expressed as follows:

$\hat{f}_{\mathrm{L,d}}(\bm{x})=\max\limits_{i=1,\ldots,N_{1}}\min_{j\in
J_{\geq,i}}u_{j}(\bm{x}),$ (13)

in which the index set $J_{\geq,i}$ is described as

$J_{\geq,i}=\\{j|u_{j}(\bm{x}_{i})\geq u_{i}(\bm{x}_{i})\\}.$ (14)

Similarly, the conjunctive lattice PWA approximation can be described as
follows:

$\hat{f}_{\mathrm{L,c}}(\bm{x})=\min\limits_{i=1,\ldots,N_{1}}\max_{j\in
J_{\leq,i}}u_{j}(\bm{x}),$ (15)

in which the index set $J_{\leq,i}$ is described as

$J_{\leq,i}=\\{j|u_{j}(\bm{x}_{i})\leq u_{i}(\bm{x}_{i})\\}.$ (16)

Compared with the full disjunctive and conjunctive lattice PWA representations
(7) and (8), respectively, we can see that the lattice PWA approximations (13)
and (15) only consider the order of local affine control laws at each sample
point. Under certain conditions as shown in Assumption 6, the lattice PWA
approximations (13) and (15) coincide with the explicit linear MPC control law
at the sample points.

######  Assumption 6.

The terms in both disjunctive and conjunctive lattice PWA approximations
satisfy,

$\min\limits_{j\in J_{\geq,i}}u_{j}(\bm{x}_{k})\leq u_{k}(\bm{x}_{k}),\forall
i,k\in\\{1,\ldots,N_{1}\\}$ (17)

and

$\max\limits_{j\in J_{\leq,i}}u_{j}(\bm{x}_{k})\geq u_{k}(\bm{x}_{k}),\forall
i,k\in\\{1,\ldots,N_{1}\\}$ (18)

######  Lemma 7.

Assume that the disjunctive and conjunctive lattice PWA approximations are
constructed through (13) and (15), respectively. Supposing that Assumption 6
holds, then we have

$\hat{f}_{\mathrm{L,d}}(\bm{x})=u^{*}(\bm{x}),\forall\bm{x}\in\Gamma(\bm{x}_{i}),\forall\bm{x}_{i}\in\mathcal{X}_{1}$
(19)

and

$\hat{f}_{\mathrm{L,c}}(\bm{x})=u^{*}(\bm{x}),\forall\bm{x}\in\Gamma(\bm{x}_{i}),\forall\bm{x}_{i}\in\mathcal{X}_{1}.$
(20)

We first prove (19) for the disjunctive case.

As all the sample points $\bm{x}_{i}$ are _interior points for some UO
regions_ , according to (17), we have

$\min\limits_{j\in J_{\geq,k}}u_{j}(\bm{x})\leq u_{i}(\bm{x}),\forall
i,k\in\\{1,\ldots,N_{1}\\},\forall\bm{x}\in\Gamma(\bm{x}_{i})$

Besides, the equality holds for $i=k$, i.e.,

$\min\limits_{j\in J_{\geq,i}}u_{j}(\bm{x})=u_{i}(\bm{x}),\forall
i\in\\{1,\ldots,N_{1}\\},\forall\bm{x}\in\Gamma(\bm{x}_{i}).$

Combining with the disjunctive lattice PWA approximation (13), we have

$\hat{f}_{\rm L,d}(\bm{x})=u_{i}(\bm{x}),\forall
k\in\\{1,\ldots,N_{1}\\},\forall\bm{x}\in\Gamma(\bm{x}_{i}).$

As $u^{*}(\bm{x})=u_{i}(\bm{x}),\forall
i\in\\{1,\ldots,N_{1}\\},\forall\bm{x}\in\Gamma(\bm{x}_{i})$, we then have
(19).

The conjunctive case can be proved similarly.

######  Remark 8.

It is noted that the lattice PWA approximation differs from the other
approximations in that the lattice PWA approximation equals the original
control law not only at the sample points, but also in the UO regions
containing the sample points as interior points, as (19) and (20) show. To
achieve this, for each sample point, not only is the value of the
corresponding control law recorded, but the specific affine expression is as
well, as Algorithm 1 shows.

###  3.3 Re-sampling

In general, when we generate a moderate number of sample points, (17) and (18)
hold. However, there are situations when the two inequalities are not valid,
the following gives the resampling method such that both (17) and (18) hold.

####  3.3.1 Guaranteeing the validity of Assumption 6

Taking the disjunctive lattice PWA approximation as an example, if (17) is
violated for some $\bm{x}_{\alpha},\bm{x}_{\beta}\in\mathcal{X}_{1}$, i.e.,

$\min\limits_{j\in
J_{\geq,\alpha}}u_{j}(\bm{x}_{\beta})>u^{*}(\bm{x}_{\beta})=u_{\beta}(\bm{x}_{\beta}),$
(21)

we can add sample points in the line segment

$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})=\lambda\bm{x}_{\alpha}+(1-\lambda)\bm{x}_{\beta},\lambda\in(0,1)$
(22)

such that (17) is satisfied for $\alpha$ and $\bm{x}_{\beta}$, which is shown
in Lemma 9.

######  Lemma 9.

Assuming that there are two points $\bm{x}_{\alpha}$ and $\bm{x}_{\beta}$ such
that (21) holds, then there must be some point
$\bm{x}_{\gamma}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$, which is
defined in (22), and the corresponding control solution $u_{\gamma}$, such
that the inequality

$u_{\gamma}(\bm{x}_{\alpha})\geq
u_{\alpha}(\bm{x}_{\alpha}),u_{\gamma}(\bm{x}_{\beta})\leq
u_{\beta}(\bm{x}_{\beta}),$ (23)

holds.

Furthermore, by adding $\bm{x}_{\gamma}$ to the sample dataset, we have

$\min\limits_{j\in J_{\geq,\alpha}}u_{j}(\bm{x}_{\beta})\leq
u_{\beta}(\bm{x}_{\beta}).$ (24)

As the optimal control solution $u^{*}$ is continuous PWA, it is still
continuous PWA when restricted to the line segment
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$.

Defining an index set $\mathrm{aff}(\bm{x}_{\alpha},\bm{x}_{\beta})$ as

$\mathrm{aff}{(\bm{x}_{\alpha},\bm{x}_{\beta})}=\\{j|\exists\bm{x}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})~{}\mbox{such
that}~{}u^{*}(\bm{x})=u_{j}(\bm{x})\\},$

i.e., the index set $\mathrm{aff}{(\bm{x}_{\alpha},\bm{x}_{\beta})}$ includes
all the indices of affine functions in $u^{*}(x)$ when restricted to
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$.

According to [13], we have

$\min\limits_{j\in S_{\geq,\alpha}}u_{j}(\bm{x})\leq
u_{\beta}(\bm{x}_{\beta}),\forall\bm{x}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta}),$

in which $S_{\geq,\alpha}$ is the index set such that

$S_{\geq,\alpha}=\\{j\in\mathrm{aff}{(\bm{x}_{\alpha},\bm{x}_{\beta})}|u_{j}(\bm{x}_{\alpha})\geq
u_{\alpha}(\bm{x}_{\alpha})\\}.$

Apparently, there should be some $\gamma\in S_{\geq,{\alpha}}$, such that (23)
is valid.

Therefore, if we add one of these $\bm{x}_{\gamma}$ to the sample point set,
we have

$\gamma\in S_{\geq,\alpha}\subset J_{\geq,\alpha}.$

Hence (24) is valid.

It is noted that for the case when (18) is violated, we have similar results.

######  Corollary 10.

Assuming that there are two points $\bm{x}_{\alpha}$ and $\bm{x}_{\beta}$ such
that the following inequality holds:

$\max\limits_{j\in
J_{\leq,\alpha}}u_{j}(\bm{x}_{\beta})<u^{*}(\bm{x}_{\beta})=u_{\beta}(\bm{x}_{\beta}),$
(25)

then there must be some point
$\bm{x}_{\gamma}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$, and the
corresponding control solution $u_{\gamma}$ satisfies the following
inequality:

$u_{\gamma}(\bm{x}_{\alpha})\leq
u_{\alpha}(\bm{x}_{\alpha}),u_{\gamma}(\bm{x}_{\beta})\geq
u_{\beta}(\bm{x}_{\beta}).$ (26)

Additionally, by adding $\bm{x}_{\gamma}$ to the sample dataset, we have

$\max\limits_{j\in J_{\leq,\alpha}}u_{j}(\bm{x}_{\beta})\geq
u_{\beta}(\bm{x}_{\beta}).$ (27)

When (21) or (25) holds, in order to construct a continuous PWA function and
to ensure the validity of (17) and (18) for every _sample point_ in the line
segment $\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$, the line segment is
recursively partitioned to generate new sample points, as Algorithm 2 shows.

######  Algorithm 2.

Recursive partitioning of line segment
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$ in case that (21) or (25) holds.

Input:Linear MPC problem, initial sample dataset
$\mathcal{X}_{1}\times\mathcal{U}_{1}$, the line segment
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$. Output:Additional sample
dataset $\mathcal{X}_{2}\times\mathcal{U}_{2}$.

1: Initialize $flag=1$, $\mathcal{X}_{2}=\emptyset$,
$\mathcal{U}_{2}=\emptyset$.

2: while flag do

3: $N_{a}=0$;

4: for $\bm{x}_{i}\in\mathcal{X}_{1}\cup\mathcal{X}_{2}$ do

5: Select corresponding
$u_{i}(\bm{x}_{i})\in\mathcal{U}_{1}\cup\mathcal{U}_{2}$ .

6: if
${\mathrm{s}ign}(u_{i}(\bm{x}_{i})-u_{i+1}(\bm{x}_{i}))=\mathrm{sign}(u_{i}(\bm{x}_{i+1})-u_{i+1}(\bm{x}_{i+1}))$
then

7: $N_{a}=N_{a}+1$.

8: Add a point $\bm{x}_{\rm new}=0.5(\bm{x}_{i}+\bm{x}_{i+1})$ to
$\mathcal{X}_{2}$.

9: Calculate corresponding affine function $u^{*}(\bm{x}_{\rm new})$ through
lines 7-9 of Algorithm 1, add to $\mathcal{U}_{2}$.

10: end if

11: end for

12: if Na=0 then

13: flag=0.

14: end if

15: end while

Lemma 11 shows that, if we add points according to Algorithm 2, (17) and (18)
are satisfied for all the sample points in
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$.

######  Lemma 11.

Given the line segment $\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$ such that
(21) or (25) holds, if we add points according to Algorithm 2, then we have

$\begin{array}[]{r}\min\limits_{j\in J_{\geq,i}}u_{j}(\bm{x}_{k})\leq
u_{k}(\bm{x}_{k}),\forall\bm{x}_{i},\bm{x}_{k}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})\cap(\mathcal{X}_{1}\cup\mathcal{X}_{2})\end{array}$
(28)

and

$\max\limits_{j\in J_{\leq,i}}u_{j}(\bm{x}_{k})\geq
u_{k}(\bm{x}_{k}),\forall\bm{x}_{i},\bm{x}_{k}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})\cap(\mathcal{X}_{1}\cup\mathcal{X}_{2}).$
(29)

After evaluating Algorithm 2, the condition

$(u_{i}(\bm{x}_{i})-u_{i+1}(\bm{x}_{i}))\cdot(u_{i}(\bm{x}_{i+1})-u_{i+1}(\bm{x}_{i+1}))\leq
0$ (30)

is satisfied for all points
$\bm{x}_{i}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$, which corresponds
to two cases as shown in Fig. 2. To generalize, here we only consider the case
when $u_{i}\neq u_{i+1}$; for the case $u_{i}=u_{i+1}$, the affine function
$u_{i}$ can connect the two points. Then, we can construct a continuous PWA
function connecting the two points $\bm{x}_{i}$ and $\bm{x}_{i+1}$, i.e.,
$\min\\{u_{i},u_{i+1}\\}$ for case 1 and $\max\\{u_{i},u_{i+1}\\}$ for case 2.

(a) Case 1. (b) Case 2.

Figure 2: Two cases satisfying (30).

Supposing that the constructed disjunctive and conjunctive continuous PWA
functions connecting $\bm{x}_{\alpha}$ and $\bm{x}_{\beta}$ are $\hat{f}_{1}$
and $\hat{f}_{2}$, respectively, then we have

$\hat{f}_{1}(\bm{x}_{k})=\hat{f}_{2}(\bm{x}_{k})=u_{k}(\bm{x}_{k}),\forall\bm{x}_{k}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})\cap(\mathcal{X}_{1}\cup\mathcal{X}_{2}).$

Defining the index sets $\mathrm{aff}_{1}(\bm{x}_{\alpha},\bm{x}_{\beta})$ and
$\mathrm{aff}_{1}(\bm{x}_{\alpha},\bm{x}_{\beta})$ as

$\mathrm{aff}_{1}{(\bm{x}_{\alpha},\bm{x}_{\beta})}=\\{j|\exists\bm{x}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})~{}\mbox{such
that}~{}\hat{f}_{1}(\bm{x})=u_{j}(\bm{x})\\},$

and

$\mathrm{aff}_{2}{(\bm{x}_{\alpha},\bm{x}_{\beta})}=\\{j|\exists~{}\bm{x}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})~{}\mbox{such
that}~{}\hat{f}_{2}(\bm{x})=u_{j}(\bm{x})\\},$

we then have

$\begin{array}[]{l}\min\limits_{j\in
J_{\geq,i}\cap\mathrm{aff}_{1}(\bm{x}_{\alpha},\bm{x}_{\beta})}u_{j}(\bm{x}_{k})\leq
u_{k}(\bm{x}_{k})\\\
\quad\quad\quad\quad\quad\quad\quad\quad\forall\bm{x}_{i},\bm{x}_{k}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})\cap(\mathcal{X}_{1}\cup\mathcal{X}_{2})\end{array}$

and

$\begin{array}[]{l}\max\limits_{j\in
J_{\leq,i}\cap\mathrm{aff}_{2}(\bm{x}_{\alpha},\bm{x}_{\beta})}u_{j}(\bm{x}_{k})\geq
u_{k}(\bm{x}_{k})\\\
\quad\quad\quad\quad\quad\quad\quad\quad\forall\bm{x}_{i},\bm{x}_{k}\in\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})\cap(\mathcal{X}_{1}\cup\mathcal{X}_{2}).\end{array}$

According to the above inequalities, we have (28) and (29).

Algorithm 2 can be run repeatedly until, for all the sample points
$(\bm{x}_{i},u_{i})\in\left(\mathcal{X}_{1}\cup\mathcal{X}_{2}\right)\times\left(\mathcal{U}_{1}\cup\mathcal{U}_{2}\right)$,
we have (17) and (18); thus Assumption 6 is satisfied and the resulting
lattice PWA approximation equals the original control solution at all sample
points and in the UO containing the sample points as interior points, as Lemma
7 shows.

A simple one-dimensional example is used to illustrate the process of
constructing the disjunctive and conjunctive lattice PWA approximations and
the re-sampling procedure.

######  Example 12.

Considering a one-dimensional continuous PWA function as in [13],

$u(x)=\left\\{\begin{array}[]{lr}\ell_{1}(x)=0.5x+0.5&x\in[0,1],\\\
\ell_{2}(x)=2x-1&x\in[1,1.5],\\\ \ell_{3}(x)=2&x\in[1.5,3.5],\\\
\ell_{4}(x)=-2x+9&x\in[3.5,4],\\\
\ell_{5}(x)=-0.5x+3&x\in[4,5],\end{array}\right.$

the plot of which is shown in Fig. 3.

Figure 3: One-dimensional continuous PWA function.

Supposing that we choose sample points as $(x_{1},u(x_{1}))=(0.5,0.75)$,
$(x_{2},u(x_{2}))=(2.5,2)$, $(x_{3},u(x_{3}))=(3.75,1.5)$, and
$(x_{4},u(x_{4}))=(4.5,0.75)$. Apparently, $(x_{2},u(x_{2}))=(2.5,2)$ is on
the boundary of some UO region, as $\ell_{1}(2.5)=\ell_{5}(2.5)$, hence we
apply a perturbation and move $x_{2}$ to 2.4, i.e., we have
$(x_{2},u(x_{2}))=(2.4,2)$.

then the disjunctive lattice PWA approximation is

$\begin{array}[]{r}\hat{f}_{\rm
L,d}=\max\\{\min\\{\ell_{1},\ell_{3},\ell_{4},\ell_{5}\\},\min\\{\ell_{3},\ell_{4}\\},\\\
\min\\{\ell_{1},\ell_{3},\ell_{4}\\},\min\\{\ell_{1},\ell_{3},\ell_{5}\\}\\}.\end{array}$

The conjunctive lattice PWA approximation is constructed as

$\begin{array}[]{r}\hat{f}_{\rm
L,c}=\min\\{\ell_{1},\max\\{\ell_{3},\ell_{5}\\},\max\\{\ell_{4},\ell_{5}\\},\max\\{\ell_{4},\ell_{5}\\}\\}.\end{array}$

It is apparent that the affine function $\ell_{2}$ has not been sampled, and
we check whether (17), (18) are satisfied and find that

$\min\limits_{j\in
J_{\geq,2}}u_{j}(x_{1})=\min\\{\ell_{3}(x_{1}),\ell_{4}(x_{1})\\}>\ell_{1}(x_{1}),$

and

$\max\limits_{j\in J_{\leq,1}}u_{j}(x_{3})=\ell_{1}(x_{3})<\ell_{3}(x_{3}).$

Hence more sample points should added between $x_{1}$ and $x_{3}$.

Applying Algorithm 2, we obtain the newly sampled points $x_{5}=1.2$ and
$u(x_{5})=2\times 1.2-1=1.4$ and now Assumption 6 is satisfied.

Until now, all of the distinct affine functions $\ell_{1},\ldots,\ell_{5}$
have been sampled, and we obtain the following lattice PWA approximations:

$\begin{array}[]{r}\hat{f}_{\rm
L,d}=\max\\{\min\\{\ell_{1},\ell_{3},\ell_{4},\ell_{5}\\},\min\\{\ell_{2},\ell_{3},\ell_{4},\ell_{5}\\},\\\
\min\\{\ell_{2},\ell_{3},\ell_{4}\\},\min\\{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\\},\min\\{\ell_{1},\ell_{2},\ell_{3},\ell_{5}\\}\\}\end{array}$

and

$\begin{array}[]{r}\hat{f}_{\rm
L,c}=\min\\{\max\\{\ell_{1},\ell_{2}\\},\max\\{\ell_{1},\ell_{2}\\},\\\
\max\\{\ell_{1},\ell_{3},\ell_{5}\\},\max\\{\ell_{4},\ell_{5}\\},\max\\{\ell_{4},\ell_{5}\\}\\}.\end{array}$

For all $x$, as
$\min\\{\ell_{2},\ell_{3},\ell_{4}\\}\geq\min\\{\ell_{2},\ell_{3},\ell_{4},\ell_{5}\\}$
and
$\min\\{\ell_{2},\ell_{3},\ell_{4}\\}\geq\min\\{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\\}$,
the disjunctive approximation $\hat{f}_{\rm L,d}$ can be further expressed as

$\begin{array}[]{r}\hat{f}_{\rm
L,d}=\max\\{\min\\{\ell_{1},\ell_{3},\ell_{4},\ell_{5}\\},\min\\{\ell_{2},\ell_{3},\ell_{4}\\},\\\
\min\\{\ell_{1},\ell_{2},\ell_{3},\ell_{5}\\}\\}.\end{array}$

Similarly, the conjunctive approximation can be rewritten as

$\hat{f}_{\rm
L,c}=\min\\{\max\\{\ell_{1},\ell_{2}\\},\max\\{\ell_{1},\ell_{3},\ell_{5}\\},\max\\{\ell_{4},\ell_{5}\\}\\}.$

Fig. 4 and 5 give the plots of $\hat{f}_{\rm L,d}$ and $\hat{f}_{\rm L,c}$,
respectively. It is apparent that for this simple example, after re-sampling,
the disjunctive approximation $\hat{f}_{\rm L,d}$ equals $u(x)$, and there is
deviation between the conjunctive approximation $\hat{f}_{\rm L,c}$ and
$u(x)$. However, both $\hat{f}_{\rm L,d}$ and $\hat{f}_{\rm L,c}$ are
identical to $u(x)$ at the sample points $x_{1},\ldots,x_{6}$ and in
$\Gamma(x_{1}),\ldots,\Gamma(x_{6})$, confirming Lemma 7.

Figure 4: Disjunctive lattice PWA approximation. Figure 5: Conjunctive lattice
PWA approximation.

###  3.4 Simplification of lattice PWA approximation

In Example 12, duplicated or redundant terms have been removed from the
lattice PWA approximation, which simplifies the approximation. When the number
of sample points $N_{s}$ is large, the evaluation of (13) and (15) are not
easy, and hence the simplification is considered in this section.

The simplification of a disjunctive lattice PWA function was addressed in
[13], for which the detailed subregions of the PWA function are known. In this
paper, the information of the subregions of the PWA function is unknown. It is
also difficult to derive the expression of the subregion polyhedra through
lattice PWA approximation. Hence, in this section, the disjunctive and
conjunctive lattice PWA approximations are simplified according to the
following rule. Assuming that the set $C$ denotes the codomain of affine
functions $u_{1},u_{2},\ldots$, the operations $\bigvee$ and $\bigwedge$ are
defined as follows,

$u_{i}\bigvee u_{j}=\max\\{u_{i},u_{j}\\},u_{i}\bigwedge
u_{j}=\min\\{u_{i},u_{j}\\}.$

It has been shown in [29] that the set $C$, together with the operations
$\bigvee$ and $\bigwedge$, constitutes a distributive lattice, and the
following property holds for all $u_{i},u_{j}\in C$:

$R1:\begin{array}[]{l}u_{i}\bigvee(u_{i}\bigwedge u_{j})=u_{i}\\\
u_{i}\bigwedge(u_{i}\bigvee u_{j})=u_{i}.\end{array}$ (31)

Actually, the above rule has been used in the simplification of lattice PWA
approximations $\hat{f}_{\rm L,d}$ and $\hat{f}_{\rm L,c}$ in Example 12.
After simplification, we get lattice PWA approximations with far less terms,
which are equivalent to the lattice PWA approximations constructed through
(13) and (15).

The process of obtaining disjunctive and conjunctive lattice PWA
approximations and simplifying them can be summarized in Algorithm 3.

######  Algorithm 3.

Construction and simplification of disjunctive and conjunctive lattice PWA
approximations.

Input:Linear MPC problem. Output:Simplified disjunctive and conjunctive
lattice PWA approximations.

1: Generate sample dataset $\mathcal{X}_{1}\times\mathcal{U}_{1}$ according to
Algorithm 1.

2: Generate additional sample dataset $\mathcal{X}_{2}\times\mathcal{U}_{2}$
according to Algorithm 2.

3: for $\bm{x}_{i}\in\mathcal{X}_{1}\cup\mathcal{X}_{2}$ do

4: Calculate the index sets $J_{\geq,i}$ and $J_{\leq,i}$ at the point
$\bm{x}_{i}$ according to Equation (14) and (16), respectively.

5: end for

6: Construct the disjunctive and conjunctive lattice PWA approximations
according to Equation (13) and (15).

7: Simplify the lattice PWA approximations according to the rule R1 (31).

The disjunctive and conjunctive lattice PWA approximations of explicit MPC
control law are illustrated using a small example of a linear discrete-time
system.

######  Example 13.

Considering the discrete-time double integrator example introduced in [30],
the system dynamics can be written as

$\begin{array}[]{rcl}x_{k+1}&=&\left[\begin{array}[]{cc}1&T_{s}\\\
0&1\end{array}\right]+\left[\begin{array}[]{c}T_{s}^{2}\\\
T_{s}\end{array}\right]u_{k}\end{array}$

where the sampling interval $T_{s}$ is 0.3. Considering the MPC problem with
$Q=\mathrm{diag}(1,0)$, $R=1$, and $P$ is the solution of the discrete-time
algebraic Riccati equation. The system constraints are $-1\leq u_{k}\leq 1$
and $-0.5\leq x_{k,2}\leq 0.5$. The region is set to be
$[-2.8,2.8]\times[-0.8,0.8]$.

To derive the lattice PWA approximation, 441 points ($21\times 21$) are
uniformly generated in the region $[-1,1]^{2}$. For all 441 points, 68 points
reach the boundary of some subregions, and we apply perturbations to move the
points to the interior of some UO regions. There are five distinct affine
functions, which are

$\begin{array}[]{l}u_{1}=-0.8082x_{1}-1.1559x_{2};\\\
u_{2}=-3.3333x_{2}-2.6667;\\\ u_{3}=-3.3333x_{2}+2.6667;\\\
u_{4}=-1;u_{5}=1.\end{array}$

We generate 441 terms, and after the simplification of terms, the following
lattice PWA approximations are obtained:

$\hat{f}_{\mathrm{L,d}}(\bm{x})=\max\\{u_{2}(\bm{x}),u_{4}(\bm{x}),\min(u_{1}(\bm{x}),u_{3}(\bm{x}),u_{5}(\bm{x}))\\}$

and

$\hat{f}_{\mathrm{L,c}}(\bm{x})=\min\\{\max(u_{1}(\bm{x}),u_{2}(\bm{x}),u_{4}(\bm{x})),u_{3}(\bm{x}),u_{5}(\bm{x})\\}.$

Readers can verify that the disjunctive and conjunctive approximations are
equivalent in $[-2.8,2.8]\times[-0.8,0.8]$. Fig. 6(a) gives the optimal MPC
controller generated by the MPT3 Toolbox [28], and the linear subregions are
also shown, see Fig. 6(b). For this example, the lattice PWA approximations
are identical to explicit MPC. In the next section, we demonstrate that if all
the affine functions have been sampled and the two lattice PWA approximations
are identical, then both of them equal the optimal MPC control law.

(a) Controller. (b) Region.

Figure 6: Explicit MPC controller in Example 13.

##  4 Approximation error and computational complexity

###  4.1 Deviations between the disjunctive and conjunctive approximations

After evaluating Algorithm 1-2, we obtain the sample point set
$(\mathcal{X}_{1}\cup\mathcal{X}_{2})\times(\mathcal{U}_{1}\cup\mathcal{U}_{2})$,
and use $\mathcal{X}=\mathcal{X}_{1}\cup\mathcal{X}_{2}$ and
$\mathcal{U}=\mathcal{U}_{1}\cup\mathcal{U}_{2}$ to denote the sample point
set. Assuming that $\mathcal{X}=\\{\bm{x}_{1},\ldots,\bm{x}_{N_{s}}\\}$, given
both the disjunctive and conjunctive lattice PWA approximations, under the
following assumption, the deviation between the approximations and the optimal
control law can be derived.

######  Assumption 14.

We assume that all the distinct affine functions have been sampled.

All the distinct affine functions for a PWA MPC optimal controller can be
obtained by collecting all critical regions and corresponding local affine
functions $u_{\rm loc}(x)$ as defined in Definition 1 and selecting the
distinct ones.

Assumption 14 can be explained as follows. Suppose we have 3 regions
$\Omega_{1},\Omega_{2},\Omega_{3}$ and the local affine functions $u_{\rm
loc,1}(\bm{x})=\ell_{1}(\bm{x}),u_{\rm loc,2}(\bm{x})=\ell_{1}(\bm{x}),u_{\rm
loc,3}(\bm{x})=\ell_{2}(\bm{x})$, then the distinct affine functions are
$\ell_{1}(\bm{x})$ and $\ell_{2}(\bm{x})$, and by sampling points in
$\Omega_{1}$ and $\Omega_{3}$, we have sampled all the distinct affine
functions.

Theorem 15 bounds the error between the lattice PWA approximations and the
original optimal control law.

######  Theorem 15.

Supposing that

$\hat{f}_{\rm
L,d}(\bm{x})=\max\limits_{i\in\\{1,\ldots,N_{s}\\}}\min\limits_{j\in
J_{\geq,i}}u_{j}(\bm{x})$

and

$\hat{f}_{\rm
L,c}(\bm{x})=\min\limits_{i\in\\{1,\ldots,N_{s}\\}}\max\limits_{i\in
J_{\leq,i}}u_{i}(\bm{x})$

are the disjunctive and conjunctive approximations of the optimal control law
$u^{*}(\bm{x})$ over the domain $\Omega$, assuming that Assumption 14 holds,
and defining

$\varepsilon=\max\limits_{\bm{x}}\left(\hat{f}_{\rm L,c}(\bm{x})-\hat{f}_{\rm
L,d}(\bm{x})\right),$ (32)

we then have

$-\varepsilon\leq\hat{f}_{\rm L,d}(\bm{x})-u^{*}(\bm{x})\leq 0$ (33)

and

$0\leq\hat{f}_{\rm L,c}(\bm{x})-u^{*}(\bm{x})\leq\varepsilon.$ (34)

Furthermore, if $\varepsilon=0$, we have

$\hat{f}_{\rm L,d}(\bm{x})=\hat{f}_{\rm
L,c}(\bm{x})=u^{*}(x),\forall\bm{x}\in\Omega.$ (35)

Assume $\mathcal{N}$ is the index of all UO regions, according to the
conclusion in [13], we have

$u^{*}(\bm{x})=\max\limits_{i\in\mathcal{N}}\min\limits_{j\in
I_{\geq,i}}u_{j}(\bm{x}),\forall\bm{x}\in\Omega,$ (36)
$u^{*}(\bm{x})=\min\limits_{i\in\mathcal{N}}\max\limits_{j\in
I_{\leq,i}}u_{j}(\bm{x}),\forall\bm{x}\in\Omega.$ (37)

in which the index sets $I_{\geq,i}$ and $I_{\leq,i}$ are defined as,

$I_{\geq,i}=\\{j|u_{j}(\bm{x})\geq
u_{i}(\bm{x}),\forall\bm{x}\in\Gamma(\bm{x}_{i})\\},$
$I_{\leq,i}=\\{j|u_{j}(\bm{x})\leq
u_{i}(\bm{x}),\forall\bm{x}\in\Gamma(\bm{x}_{i})\\}.$

If Assumption 14 holds, i.e., all the distinct affine functions have been
sampled, then for a UO region $\Gamma(\bm{x}_{i})$, as the order of affine
functions remains unchanged in the UO region, the set $I_{\geq,i}$ is
identical to $J_{\geq,i}$ defined in (14). Similarly, the sets $I_{\leq,i}$
and $J_{\leq,i}$ defined in (16) are equivalent. Therefore, for all
$\bm{x}\in\Omega$ and all $i\in\mathcal{N}$, the following inequalities hold:

$\min\limits_{j\in J_{\geq,i}}u_{j}(\bm{x})\leq u^{*}(\bm{x}),$ (38)

and

$\max\limits_{j\in J_{\leq,i}}u_{j}(\bm{x})\geq u^{*}(\bm{x}).$ (39)

As the sampled UO regions are only a subset of all UO regions, i.e.,

$\\{1,\ldots,N_{s}\\}\subset\mathcal{N},$

we have inequalities (38) and (39) for all $\bm{x}\in\Omega$ and all
$i\in\\{1,\ldots,N_{s}\\}$.

Then the following is valid:

$\hat{f}_{\rm L,d}(\bm{x})\leq u^{*}(\bm{x})\leq\hat{f}_{\rm
L,c}(\bm{x}),\forall\bm{x}\in\Omega.$

According to (32), we have (33) and (34).

Furthermore, if $\varepsilon=0$, then both approximations are identical to the
optimal control law in the region $\Omega$, i.e., (35) holds.

As Theorem 15 indicates, in order to guarantee that the lattice PWA
approximations are error-free, two conditions should be satisfied, the first
is the validity of Assumption 14, and the second is that $\varepsilon$ as
defined in (32) is zero.

######  Remark 16.

As a necessary condition for the validity of (35), i.e., the lattice PWA
approximations are error-free, Assumption 14 requires that all the distinct
affine functions have been sampled. It is noted that this is not the same as
the condition that all the critical regions should be identified. As [13]
shows, the number of distinct affine functions in explicit MPC is generally
far more less than the number of critical regions, as some critical regions
can share the same affine function.

####  4.1.1 A necessary condition for the validity of Assumption 14

In general, it is not easy to check whether Assumption 14 is satisfied, which
is vital for the validity of (38) and (39) for all $\bm{x}\in\Omega$ and
$i\in\\{1,\ldots,N_{s}\\}$. Hence we check whether the relaxed condition

$\min_{j\in J_{\geq,i}}u_{j}(\bm{x})\leq\max_{j\in
J_{\leq,i}}u_{j}(\bm{x}),\forall
i\in\\{1,\ldots,N_{s}\\},\forall\bm{x}\in\Omega$ (40)

holds, which is a direct result of (38) and (39).

To check whether (40) is satisfied, for any $i,k\in\\{1,\ldots,N_{s}\\}$, the
following optimization problem is solved:

$\begin{array}[]{rl}\min\limits_{\bm{x}}&\max\limits_{j\in
J_{\leq,i}}u_{j}(\bm{x})-\min\limits_{j\in J_{\geq,k}}u_{j}(\bm{x})\\\
s.t.&\bm{x}\in\Omega.\end{array}$ (41)

If the optimal value for all $i,k\in\\{1,\ldots,N_{s}\\}$ is nonnegative, then
(40) holds.

The optimization problem (41) can be transformed into an equivalent linear
programming (LP) problem:

$\begin{array}[]{rl}\min\limits_{\bm{x},y_{1},y_{2}}&y_{1}+y_{2}\\\
s.t.&\bm{x}\in\Omega,\\\ {}\hfil&u_{j}(\bm{x})\leq y_{1},\forall j\in
J_{\leq,i},\forall i\in\\{1,\ldots,N_{s}\\}\\\ {}\hfil&u_{j}(\bm{x})\geq-
y_{2},\forall j\in J_{\geq,k},\forall k\in\\{1,\ldots,N_{s}\\}\end{array}$
(42)

which is easy to solve as $\Omega$ is a polyhedron. If we find a point
$\bm{x}$ and an index pair $i,k$, such that the cost in (42) is negative,
which means that (40) is violated, we can generate more sample points and get
more distinct affine functions, as Lemma 17 shows.

######  Lemma 17.

Suppose there are some $\bm{x}_{\gamma}$ and indices $\alpha,\beta$ such that
$y_{1}(\bm{x}_{\gamma})+y_{2}(\bm{x}_{\gamma})<0$, then more sample points can
be generated in the line segment $\mathcal{L}(\bm{x}_{i},\bm{x}_{\gamma})$ or
$\mathcal{L}(\bm{x}_{\gamma},\bm{x}_{k})$ according to Algorithm 2, and more
distinct affine functions can be sampled.

As $y_{1}(\bm{x}_{\gamma})+y_{2}(\bm{x}_{\gamma})<0$ for some $\alpha,\beta$,
we have

$\max\limits_{j\in J_{\leq,\alpha}}u_{j}(\bm{x}_{\gamma})<\min\limits_{j\in
J_{\geq,\beta}}u_{k}(\bm{x}_{\gamma}).$ (43)

Then at least one of the inequalities

$\min\limits_{j\in
J_{\geq,\alpha}}u_{j}(\bm{x}_{\gamma})>u^{*}(\bm{x}_{\gamma})$ (44)

and

$\max\limits_{j\in
J_{\leq,\beta}}u_{j}(\bm{x}_{\gamma})<u^{*}(\bm{x}_{\gamma})$ (45)

is valid. This is apparent, since, if both (44) and (45) do not hold, we have

$\min\limits_{j\in J_{\geq,\alpha}}u_{j}(\bm{x}_{\gamma})\leq\max\limits_{j\in
J_{\leq,\beta}}u_{j}(\bm{x}_{\gamma}),$

which contradicts (43).

The checking of (44) and and (45) is very simple, as $\bm{x}_{\gamma}$ and
$u^{*}(\bm{x}_{\gamma})$ are both known, and computing $\min\limits_{j\in
J_{\geq,\alpha}}u_{j}(\bm{x}_{\gamma})$ only needs arithmetic operations,
i.e., comparing the functions values of several affine functions at the point
$\bm{x}_{\gamma}$. If (44) holds, then according to Lemma 11, sample points
can be added to the line segment
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\gamma})$ as in Section 3.3.1 to ensure
that

$\min\limits_{j\in J_{\geq,i}}u_{j}(\bm{x}_{k})\leq u_{k}(\bm{x}_{k})$

for all the sample points in $\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\gamma})$.
In this process, more distinct affine functions have been identified.

For the conjunctive case, if (45) is valid, we can also add sample points in
the line segment $\mathcal{L}(\bm{x}_{\beta},\bm{x}_{\gamma})$ according to
Lemma 11 such that

$\max\limits_{j\in J_{\leq,i}}u_{j}(\bm{x}_{k})\geq u_{k}(\bm{x}_{k})$

for all the sample points in $\mathcal{L}(\bm{x}_{\beta},\bm{x}_{\gamma})$.

By doing so, more sample points as well as more distinct affine functions have
been sampled, which is necessary for the validity of Assumption 14.

######  Remark 18.

It is noted that although the relaxed condition (40) is only a necessary
condition for the validity of Assumption 14, the checking of (40) is with
respect to the entire domain $\Omega$, not only the sample points. And we
discovered in our numerical experiments that in most cases, if (40) holds in
the domain of interest, all the distinct affine functions active in that
region have been identified.

####  4.1.2 Checking whether the 2 approximations are identical

Both the disjunctive approximation $\hat{f}_{\rm L,d}$ and conjunctive
approximation $\hat{f}_{\rm L,c}$ are continuous PWA functions, so whether
$\hat{f}_{\rm L,d}=\hat{f}_{\rm L,c}$ can be checked in each linear subregion
of both $\hat{f}_{\rm L,d}$ and $\hat{f}_{\rm L,C}$. When the number of
literals and terms in the lattice PWA approximations are large, it is not easy
to identify all the linear subregions. So here we resort to a statistical
method, i.e., generating a huge number of i.i.d. sample points, which
constitute a validation dataset $\mathcal{X}_{\rm
validate}=\\{\bm{x}_{i},{i=1},\ldots,{N_{v}}\\}$. For each sample point, we
define an indicator function as

$I(\bm{x}_{i}):=\left\\{\begin{array}[]{cc}1&\mbox{if}~{}\hat{f}_{\rm
L,d}(\bm{x}_{i})=\hat{f}_{\rm L,c}(\bm{x}_{i})\\\ 0&\mbox{if}~{}\hat{f}_{\rm
L,d}(\bm{x}_{i})\neq\hat{f}_{\rm L,c}(\bm{x}_{i}),\end{array}\right.$

and then the random variables $I(\bm{x}_{i}),i=1,\ldots,N_{v}$ are also i.i.d.
Denoting the probability for $I(\bm{x}_{i}=1)$ as $\mu$, i.e.,
$\mathbb{P}[I(\bm{x}_{i})=1]=\mu$, then according to Hoeffding’s inequality
[31], we have

$\mathbb{P}[|\mu-\bar{I}|\geq\epsilon]\leq 2\exp(-2N_{v}\epsilon^{2}),$

in which $\bar{I}=\frac{1}{N_{v}}\sum\limits_{k=1}^{N_{v}}I(\bm{x}_{k})$.
Therefore, we have

$\mathbb{P}[\mu\geq\bar{I}-\epsilon]>1-2\exp(-2N_{v}\epsilon^{2}),$

meaning that with confidence $1-2\exp(-2N_{v}\epsilon^{2})$, the probability
that the lattice PWA approximations $\hat{f}_{\rm L,d}$ and $\hat{f}_{\rm
L,c}$ are identical is larger than $\bar{I}-\epsilon$. If $\bar{I}=1$, then by
setting a small enough threshold $\epsilon$, we can say that with confidence
$1-2\exp(-2N_{v}\epsilon^{2})$, the lattice PWA approximations $\hat{f}_{\rm
L,d}$ and $\hat{f}_{\rm L,c}$ are almost identical, and thus both equal the
optimal control law. For example, if $\epsilon=10^{-3}$, then $N_{v}\geq
5\times 10^{6}$ can ensure that the confidence is almost 1.

###  4.2 Complexity analysis

####  4.2.1 Online evaluation

Assuming that there are $\tilde{N}$ terms in the final approximation,
according to [12], the worst-case online evaluation complexity is
$O(\tilde{N}^{2})$. In general, we have $\tilde{N}\ll N_{s}$, and hence the
online evaluation is very fast.

####  4.2.2 Storage requirements

Assuming that the disjunctive lattice PWA approximation has $\tilde{N}$ terms
and $M$ literals (i.e. the number of distinct affine functions is $M$), we
must store $(n_{x}+1)\cdot M$ real numbers and
$\sum\limits_{i=1}^{\tilde{N}}|J_{\geq,i}|$ integer numbers, in which
$|J_{\geq,i}|$ is the number of elements in the set $J_{\geq,i}$. As
$|J_{\geq,i}|\leq M$, in total $(n_{x}+1)\cdot M$ real numbers and
$M\cdot\tilde{N}$ integer numbers must be stored.

In many cases, we have $\tilde{N}\ll N_{s}$, and hence the storage requirement
for the disjunctive lattice PWA approximation is very small.

For the conjunctive lattice PWA approximation, we achieve the same result.

####  4.2.3 Offline complexity

The offline time complexity for deriving equivalent disjunctive and
conjunctive lattice PWA approximations can be summarized as follows.

The time complexity consists of two parts. One concerns the training points
sampling and re-sampling in order to obtain lattice PWA approximations and the
other is the complexity of checking whether the lattice PWA approximations are
error-free. Lemma 19 describes the worst-case offline time complexity.

######  Lemma 19.

Assuming that the sample domain is a hyperbox, then the worst-case complexity
of deriving the disjunctive and conjunctive lattice PWA approximations is
$O(N_{v}\cdot N_{s}^{2})$, in which $N_{s}$ and $N_{v}$ are the numbers of
sample points in $\mathcal{X}$ and $\mathcal{X}_{\rm validate}$, respectively.
As indicated previously, the offline complexity comes from evaluating
Algorithm 3 and checking whether the two lattice PWA approximations are
equivalent.

For Algorithm 3, the evaluation of Algorithm 1 and 2, and the simplification
of lattice PWA approximations account for the major part of the computational
complexity.

The number of sample points in Algorithm 1 can be calculated as
$N_{1}=\prod\limits_{i=1}^{n}\frac{b_{i}-a_{i}}{\delta_{i}}$, in which
$[a_{i},b_{i}]$ is the sample range for the $i$-th component, and
correspondingly, $\delta_{i}$ is the length of the sample interval. The
complexity of evaluating Algorithm 1 includes solving $N_{1}$ convex quadratic
programming problems and corresponding KKT conditions, which are basically
solving linear equations. The time complexity of solving $N_{1}$ convex
quadratic programming problems with $N_{p}\cdot n_{u}$ decision variables is
approximately $O(N_{1}\cdot L^{2}(N_{p}\cdot n_{u})^{4})$ by using an
interior-point algorithm [32], in which $L$ is the bit length of the quadratic
programming problem. The dominant algorithmic operation in solving the KKT
conditions is solving $N_{1}$ matrix inversion problems, the worst-case
complexity of which is $O(N_{1}|\mathcal{A}^{*}|^{3})$ using the Gauss-Jordan
elimination algorithm, in which $|\mathcal{A}^{*}|$ is the number of active
constraints. As $|\mathcal{A}^{*}|\leq p$, where $p$ is the number of
constraints in QP (3), the worst-case complexity for solving the KKT
conditions is $O(N_{1}p^{3})$.

We now discuss the worst-case complexity of evaluating Algorithm 2. For two
points $\bm{x}_{\alpha}$ and $\bm{x}_{\beta}$, if (38) is violated, the
evaluation of Algorithm 2 is basically a binary search method for identifying
omitted subregions in the line segment
$\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$ as well as calculating affine
functions for newly generated points. Assume the maximum number of subregions
appearing in $\mathcal{L}(\bm{x}_{\alpha},\bm{x}_{\beta})$ is
$M_{\alpha,\beta}$, and generally $M_{\alpha,\beta}<N_{1}$. The binary
searching of the subregions then yields a worst-case complexity of
$O(\log_{2}N_{1})$ and the affine function calculation yields a worst-case
complexity of $O(p^{3}\log_{2}N_{1})$, in which $p$ is the number of
constraints in QP (3). Supposing that there are $N_{t}$ point pairs such that
(38) is violated, then the worst-case complexity is $O(N_{t}\cdot
p^{3}\log_{2}N_{1})$. In general, the number $N_{t}$ is closely related to the
sample grid size $\delta_{i}$, which then depends on the number of sample
points $N_{1}$ in Algorithm 1. A larger $N_{1}$ will result in a smaller
number of $N_{t}$; hence, the complexity of Algorithm 2 can be decreased by
increasing the complexity of Algorithm 1.

After evaluating Algorithm 1-2, there are $N_{s}$ sample points. The
simplification procedure requires the comparison of the sets $J_{\geq,i}$
($J_{\leq,i}$) for $i=1,\ldots,N_{s}$, which at most yields
$\left({}^{~{}2}_{N_{s}}\right)=\frac{N_{s}(N_{s}-1)}{2}$ times comparisons.
For each comparison, at most $M^{2}$ literals need to be compared. Hence the
worst-case complexity for the simplification is $O(M^{2}N_{s}^{2})$.

After simplification, we have only $\tilde{N}$ terms, which is generally much
less than $N_{s}$, i.e., $\tilde{N}\ll N_{s}$, then for the process of
checking whether the lattice PWA approximations are error-free, one has to
solve at most $\tilde{N}$ LP problems and generate $N_{v}$ validation points.
Assuming that $L$ is the bit length of the LP problem (42), then the worst-
case complexity for solving $\tilde{N}^{2}$ LP problems is
$O(\tilde{N}^{2}\cdot n_{x}^{3.5}\cdot L)$. For the sampling of $N_{v}$ points
in the lattice PWA approximations, as indicated in Section 4.2.1, the function
evaluation process of lattice PWA approximations has a worst-case complexity
of $O(\tilde{N}^{2})$, in which $\tilde{N}<N_{s}$, hence the worst-case
complexity of generating validation points is $O(N_{v}\cdot\tilde{N}^{2})$. As
$N_{v}$ is large, we have $O(\tilde{N}^{2}\cdot n_{x}^{3.5}\cdot L)\ll
O(N_{v}\cdot\tilde{N}^{2})$; hence the worst-case complexity of validation is
$O(N_{v}\cdot\tilde{N}^{2})$.

In general, $N_{1}<N_{s}$, $p\ll N_{1}$, $N_{t}\ll N_{1}$, $M\ll N_{v}$, and
$\tilde{N}\ll N_{s}$, and then the total worst-case complexity is
$O(N_{v}\cdot N_{s}^{2})$.

It is noted that $O(N_{v}\cdot N_{s}^{2})$ is the worst-case complexity. In
the simulation results, we can see that the offline calculation is actually
not time consuming.

##  5 Simulation results

######  Example 20.

Consider an example taken from [18], which is the inverted pendulum on a cart.
The state consists of the angle and the angle speed of the pole, i.e., $\Phi$,
$\dot{\Phi}$, respectively, and the position and speed of the cart, i.e., $s$
and $\dot{s}$, respectively. The constraints for the state are
$\|\bm{x}\|_{\infty}^{T}\leq[1,1.5,0.35,1.0]^{T}$. The input is the force, the
constraint of which is $|u|\leq 1$. The discrete-time dynamic is given by

$A=\left[\begin{array}[]{cccc}1&0.1&0&0\\\ 0&0.9818&0.2673&0\\\ 0&0&1&0.1\\\
0&-0.0455&3.1182&1\end{array}\right],B=\left[\begin{array}[]{c}0\\\ 0.1818\\\
0\\\ 0.4546\end{array}\right].$

The prediction horizon is taken to be $N=10$. The value of matrices in the
cost function is $Q=\rm diag\\{2,2,2,2\\}$, $R=1$, and $P=0$. According to the
MPT3 toolbox, the optimal control solution is a PWA function of the state
$\bm{x}$, with 2,271 polyhedral regions.

In [18], 88,341 samples were generated to train an approximated controller,
which is basically a deep PWL neural network. Here, in order to construct the
disjunctive and conjunctive lattice PWA approximations, only $8^{4}=4096$
samples are generated uniformly in the region
$\Omega=[-0.6~{}0.6]\times[-0.9~{}0.9]\times[-0.21~{}0.21]\times[-0.6~{}0.6]$.
In this case, there are only 13 distinct affine functions, which is due to the
choice of the region. For a larger region, say
$\Omega_{2}=[-0.9~{}0.9]\times[-1.2~{}1.2]\times[-0.25~{}0.25]\times[-0.9~{}0.9]$,
there are 47 distinct affine functions. However, there are a significant
number of infeasible state points in $\Omega_{2}$, and according to the state
trajectory in [18], the region $\Omega$ is enough.

The evaluation of Algorithm 1-2 results in the disjunctive and conjunctive
lattice PWA approximations, both with six terms. All the computations in this
paper are implemented through MatLab 2016b (MathWorks, USA) on an Apple M1 Max
computer.

The entire offline calculation time is 41.25s, i.e., Algorithm 1, 2, 3 takes
40.25s, 0.66s and 0.33s, respectively. The number of parameters for both
approximations is 143, while the number of parameters stored in the MPT3
solution is 104,535. The average online evaluation time for the MPT solution,
the disjunctive as well as conjunctive lattice PWA approximations are 0.0052s,
$8.96\times 10^{-5}$s, and $8.33\times 10^{-5}$s, respectively. It is apparent
that the approximation results in a much lower online computational burden.

It has been tested through $5\times 10^{6}$ test data points that the two
lattice PWA approximations are identical in the region $\Omega$. By setting
$\epsilon=10^{-3}$, it can be concluded that with confidence $\delta=0.9999$
the probability that the approximated lattice PWA control laws equal the
optimal control is larger than 0.999. For the $5\times 10^{6}$ points, the
optimal explicit linear MPC control law is also calculated and it is found
that the lattice PWA approximations are error-free in $\Omega$.

Fig. 7 shows one exemplary closed-loop simulation of the example, and we can
see from the figure that the optimal state trajectory and the trajectory with
the lattice PWA approximations as inputs are identical.

Figure 7: One exemplary closed-loop simulation of Example 20.

######  Example 21.

Consider an example taken from [5], which is an mpQP problem constructed from
the typical MPC setup of the form

$\displaystyle\min\limits_{u_{0},\ldots,u_{N-1}}\bm{x}_{N}^{T}P\bm{x}_{N}+\sum\limits_{k=0}^{N-1}\bm{x}_{k}^{T}Q\bm{x}_{k}+u_{k}^{T}Ru_{k}$
(46a) $\displaystyle~{}~{}~{}~{}s.t.\bm{x}_{k+1}=A\bm{x}_{k}+Bu_{k}$ (46b)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}\bm{x}\in\mathcal{X},u\in\mathcal{U}$
(46c)

with $\bm{x}\in\mathbb{R}^{10}$, $u\in\mathbb{R}$, $P=Q=I_{10},R=1$,
$\mathcal{X}=\\{x|-10\leq x_{i}\leq 10,i=1,\ldots,n_{x}\\}$,
$\mathcal{U}=\\{-1\leq u\leq 1\\}$. The prediction model is obtained by
discretizing the model $1/(s+1)^{10}$ with sampling time of 1 $\mathrm{s}$ and
then converting the discretized model into state space form. The prediction
horizon is taken to be $N=10$.

Apparently this problem is more complex than Example 20, as the dimension is
much higher. As [5] shows, the $mpt_{solve}$ function in MPT3 Toolbox failed
to solve this problem offline, i.e., MATLAB ran out of memory.

Here, in order to construct the disjunctive and conjunctive lattice PWA
approximations, we have to generate sample points in the 10-dimensional
domain, which may suffer from the curse-of-dimensionality if the sample points
are generated uniformly in the domain. Instead, we generate sample points on
several trajectories, which reflect the most possible states of the dynamic
system under optimal MPC control law. Specifically, 300 initial points are
generated in the domain $[-2,2]^{10}$, then for the 300 trajectories we have
7802 sample points.

In this case, there are 16 distinct affine functions. The evaluation of
Algorithm 1-2 results in the disjunctive and conjunctive lattice PWA
approximations, the number of terms is 9 and 8, respectively. The offline as
well as online complexity of lattice PWA approximations when the prediction
horizon $N$ is 10 are shown in Table 1.

Table 1: Performance of lattice PWA approximations on Example 21, in which $N_{s}$ and $\\#\rm{para}$ are the number of sample points and parameters, and $t_{\rm off}$ and $t_{\rm on}$ denote offline and online computation time. Method | $N$ | $N_{s}$ | #para |  $t_{\rm off}$[s] |  $t_{\rm on}$[s]
---|---|---|---|---|---
Disjunctive | 10 | 7808 | 320 | 84.65 | $1.69\times 10^{-4}$
Conjunctive |  | 7808 | 304 | 84.65 | $1.40\times 10^{-4}$
Online MPC |  |  |  |  | 0.0077
Method in [5] |  |  |  | 72.85 |
Disjunctive | 20 | 10188 | 418 | 229.24 | $2.71\times 10^{-4}$
Conjunctive |  | 10188 | 418 | 229.24 | $2.53\times 10^{-4}$
Online MPC |  |  |  |  | 0.015
Method in [5] |  |  |  |  |

The entire offline calculation time for generating the disjunctive and
conjunctive lattice PWA approximations is 84.65s, i.e., Algorithm 1,2,3
accounts for 80.3016s, 2.8734s, 1.4776s, respectively. The number of
parameters for the disjunctive and conjunctive approximations are 320 and 304,
respectively. The average online evaluation times (over $10^{4}$ trials) for
the online QP, disjunctive as well as conjunctive lattice PWA approximations
are $0.0077$, $1.69\times 10^{-4}$s, and $1.40\times 10^{-4}$s, respectively.
It is apparent that the approximations result in a very low online
computational burden.

It has been tested through $5\times 10^{6}$ test data points that the two
lattice PWA approximations are identical in the region $\Omega$. By setting
$\epsilon=10^{-3}$, it can be concluded that with confidence $\delta=0.9999$
the probability that the approximated lattice PWA control laws equal the
optimal control law is larger than 0.999. For the $5\times 10^{6}$ points, the
optimal linear MPC control law is also calculated and it is found that the
lattice PWA approximations are error-free in $\Omega$.

Fig. 8 shows one exemplary closed-loop simulation of the example, and we can
see from the figure that the optimal state trajectory and the trajectory with
the lattice PWA approximations as inputs are identical.

Figure 8: One exemplary closed-loop simulation of Example 21.

As indicated in [5], the MPT3 Toolbox failed to solve the problem. Moreover
the combinatorial approach proposed in [5] successfully enumerating all
optimal active sets and creating corresponding critical regions, takes 72.85s.
It is noticed that although the calculation time for our procedure is slightly
longer, we obtain the final continuous PWA expression of the optimal control
law in a simplified form, which is easier to be implement online.

To demonstrate more clearly the efficacy of the lattice PWA approximations,
the prediction horizon $N$ is extended to 20, and the complexity of lattice
PWA approximations is also listed in Table 1. Besides, we also generate
$5\times 10^{6}$ validation points to show that with confidence
$\delta=0.9999$ the probability that the approximated lattice PWA control laws
equal the optimal control law is larger than 0.999. When $N=20$, [5] did not
provide a result either. The lattice PWA approximations in this case takes
229.2417s to get error-free approximations in $\Omega=[-2,2]^{10}$, showing
that the lattice PWA approximations scale well with the problem size for this
example.

##  6 Conclusions and Future work

In this paper, we have presented disjunctive and conjunctive lattice PWA
approximations of the explicit linear MPC control law by sampling and
resampling in the state domain. The lattice PWA approximated and exact control
laws are identical for sample points and in UO regions that contain the sample
points as interior points. Furthermore, under the assumption that all the
affine functions have been identified in the domain of interest, the
disjunctive lattice PWA approximation is always smaller than the original
optimal control law, while the conjunctive lattice PWA approximation is always
larger. Then if the disjunctive and conjunctive lattice PWA approximations are
identical, both are equivalent to the optimal control law. The two kinds of
lattice PWA approximations have been simplified to further reduce the storage
and online evaluation complexity. The complexity of the online and offline
approximation as well as the storage requirements, have been analyzed.
Simulation results show that with a moderate number of sample points we can
obtain statistically error-free lattice PWA approximations that are calculated
with relatively small computational cost.

In the future, apart form the information of local affine functions, the
information of critical regions will also be used to construct the lattice PWA
approximation. Moreover, for more general domain of interest, the
approximation error will be considered, which can be used to derive
corresponding feasibility and stability analysis. Besides, practical
applications will also be considered.

## References

  * [1] T. Samad.  A survey on industry impact and challenges thereof.  IEEE Control Systems Magazine, 37(1):17–18, 2017.
  * [2] A. Bemporad, M. Morari, V. Dua, and E. N. Pistikopoulos.  The explicit linear quadratic regulator for constrained systems.  Automatica, 38(1):3–20, 2002.
  * [3] M. Kvasnica, B. Takács, J. Holaza, and S. Di Cairano.  On region-free explicit model predictive control.  In Proceedings of the 2015 IEEE 54th Conference on Decision and Control (CDC), pages 3669–3674, Osaka, Japan, 2015.
  * [4] A. Gupta, S. Bhartiya, and P. Nataraj.  A novel approach to multiparametric quadratic programming.  Automatica, 47(9):2112–2117, 2011.
  * [5] P. Ahmadi-Moshkenani, T. A. Johansen, and S. Olaru.  Combinatorial approach toward multiparametric quadratic programming based on characterizing adjacent critical regions.  IEEE Transactions on Automatic Control, 63(10):3221–3231, 2018\.
  * [6] F. Borrelli, M. Baotić, J. Pekar, and G. Stewart.  On the computation of linear model predictive control laws.  Automatica, 46(6):1035–1041, 2010.
  * [7] M. Herceg, S. Mariéthoz, and M. Morari.  Evaluation of piecewise affine control law via graph traversal.  In Proceedings of the European Control Conference, pages 3083–3088, 2013.
  * [8] F. J. Christophersen, M. Kvasnica, C. N. Jones, and M. Morari.  Efficient evaluation of piecewise control laws defined over a large number of polyhedra.  In Proceedings of the European Control Conference, pages 2360–2367, 2007.
  * [9] F. Bayat, T. A. Johansen, and A. Jalali.  Flexible piecewise function evaluation methods based on truncated binary search trees and lattice representation in explicit MPC.  IEEE Transactions on Control Systems Technology, 20(3):632–640, 2012.
  * [10] N. A. Nguyen, M. Gulan, S. Olaru, and P. Rodriguez-Ayerbe.  Convex lifting: Theory and control applications.  IEEE Transactions on Automatic Control, 63(5):1243–1258, 2017.
  * [11] M. Gulan, G. Takács, N. A. Nguyen, S. Olaru, P. Rodríguez-Ayerbe, and Rohal’-Ilkiv B.  Efficient embedded model predictive vibration control via convex lifting.  IEEE Transactions on Control Systems Technology, 27(1):48–62, 2019\.
  * [12] C. Wen, X. Ma, and B. E. Ydstie.  Analytical expression of explicit MPC solution via lattice piecewise-affine function.  Automatica, 45(4):910–917, 2009.
  * [13] J. Xu, T. J. J. van den Boom, B. De Schutter, and S. Wang.  Irredundant lattice representations of continuous piecewise affine functions.  Automatica, 70:109–120, 2016.
  * [14] A. Bemporad and C. Filippi.  Suboptimal explicit MPC via approximate multiparametric quadratic programming.  In Proceedings of the 40th IEEE Conference on Decision and Control, 2001, volume 5, pages 4851–4856. IEEE, 2001.
  * [15] A. Bemporad, A. Oliveri, T. Poggi, and M. Storace.  Ultra-fast stabilizing model predictive control via canonical piecewise affine approximations.  IEEE Transactions on Automatic Control, 56(12):2883–2897, 2011\.
  * [16] G. Goebel and F. Allgöwer.  Semi-explicit MPC based on subspace clustering.  Automatica, 83:309–316, 2017.
  * [17] A. Chakrabarty, V. Dinh, M. J. Corless, A. E. Rundell, S. H. Żak, and G. T. Buzzard.  Support vector machine informed explicit nonlinear model predictive control using low-discrepancy sequences.  IEEE Transactions on Automatic Control, 62(1):135–148, 2016.
  * [18] B. Karg and S. Lucia.  Efficient representation and approximation of model predictive control laws via deep learning.  IEEE Transactions on Cybernetics, 50(9):3866–3878, 2020.
  * [19] L. Csekő, M. Kvasnica, and B. Lantos.  Explicit MPC-based RBF neural network controller design with discrete-time actual Kalman filter for semiactive suspension.  IEEE Transactions on Control Systems Technology, 23(5):1736–1753, 2015.
  * [20] S. Summers, C. N. Jones, J. Lygeros, and M. Morari.  A multiresolution approximation method for fast explicit model predictive control.  IEEE Transactions on Automatic Control, 56(11):2530–2541, 2011\.
  * [21] S. Gros and M. Zanon.  Data-driven economic NMPC using reinforcement learning.  IEEE Transactions on Automatic Control, 65(2):636–648, 2020.
  * [22] F. Scibilia, S. Olaru, and M. Hovd.  Approximate explicit linear MPC via Delaunay tessellation.  In 2009 European Control Conference (ECC), pages 2833–2838, 2009\.
  * [23] A. Pavlov, I. Shames, and C. Manzie.  Minimax strategy in approximate model predictive control.  Automatica, 111:108649, 2020.
  * [24] J. Xu.  Lattice piecewise affine approximation of explicit linear model predictive control.  In Proceedings of the 60th IEEE Conference on Decision and Control, 2021.
  * [25] L. O. Chua and A. C. Deng.  Canonical piecewise-linear representation.  IEEE Transactions on Circuits and Systems, 35(1):101–111, 1988\.
  * [26] J. Nocedal and S. Wright.  Numerical Optimization.  Springer-Verlag, New York, NY, USA, 2006.
  * [27] F. Borrelli.  Constrained Optimal Control of Linear and Hybrid Systems, volume 290.  New York: Springer, 2003.
  * [28] M. Herceg, M. Kvasnica, C. N. Jones, and M. Morari.  Multi-parametric toolbox 3.0.  In Proceedings of the European Control Conference, pages 502–510, 2013.  http://control.ee.ethz.ch/~mpt.
  * [29] J. M. Tarela, J. M. Perez, and V. Aleixandre.  Minimization of lattice polynomials on piecewise linear functions (Part I).  Mathematics and Computers in Simulation, 17(2):79–85, 1975.
  * [30] M. Johansson.  Piecewise linear control systems: A computational approach.  Springer-Verlag New York, Inc. Secaucus, NJ, USA, 2003.
  * [31] W. Hoeffding.  Probability inequalities for sums of bounded random variables.  In 409-426, editor, The Collected Works of Wassily Hoeffding. Springer, 1994.
  * [32] Y. Ye and E. TSE.  An extension of Karmarkar’s projective algorithm for convex quadratic programming.  Mathematical Programming, 44(1):157–159, 1989.