# Data-driven approximation for feasible regions in nonlinear model predictive
control

Yuanqiang Zhou  Dewei Li  Yugeng Xi  Yunwen Xu  Department of Automation,
Shanghai Jiao Tong University, and Key Laboratory of System Control and
Information Processing, Ministry of Education of China, Shanghai, 200240,
China.
e-mail: {zhouyuanqiang, dwli, ygxi<EMAIL_ADDRESS>
###### Abstract

This paper develops a data-driven learning framework for approximating the
feasible region and invariant set of a nonlinear system under the nonlinear
Model Predictive Control (MPC) scheme. The developed approach is based on the
feasibility information of a point-wise data set using low-discrepancy
sequence. Using kernel-based Support Vector Machine (SVM) learning, we
construct outer and inner approximations of the boundary of the feasible
region and then, obtain the feasible region of MPC for the system.
Furthermore, we extend our approach to the perturbed nonlinear systems using
set-theoretic method. Finally, an illustrative numerical example is provided
to show the effectiveness of the proposed approach.

###### keywords:

Nonlinear model predictive control; support vector machine; set-theoretic
method.

scaletikzpicturetowidth[1]

††thanks: This work was supported by the National Science Foundation of China
(Grant No. 61973214, 61433002, 61590924, 61573239).

##  1 Introduction

When it comes to MPC strategy for continuous-time systems, we periodically
compute an input signal from solving an optimization problem and apply this
signal to the system over a bounded horizon, see e.g. Rawlings and Mayne
(2009); Xi and Li (2019). With the improving of control theory and technology,
the new problems which involve large transients often require Nonlinear MPC
(NMPC), see e.g. Allgöwer and Zheng (2012); Grüne and Pannek (2017). NMPC for
continuous-time systems with input constraints is proposed, see e.g. Chen and
Allgöwer (1998), which considers a nonlinear system given by

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x(t),u(t)),\quad x(0)=\bar{x},\quad
t\in\mathbb{R}_{\geq 0},$ (1) $\displaystyle\quad x(t)\in\mathbb{R}^{n},\quad
u(t)\in\mathbb{U}\subset\mathbb{R}^{m},$ (2)

where ${x}$ and $u$ are the state vector and the control input, respectively;
$\mathbb{U}$ is the input constraint set and the mapping $f(\cdot,\cdot)$
satisfy $f(0,0)=0$. Besides,
$f:\mathbb{R}^{n}\times\mathbb{U}\to\mathbb{R}^{n}$ is assumed to be twice
continuously differentiable, which leads to unique solutions under the piece-
wise right-continuous control inputs, and $\mathbb{U}$ is assumed to be
compact, convex, and contains $0$ in its interior.

The quasi-infinite horizon NMPC developed in Chen and Allgöwer (1998) consists
of designing a control invariant region $\Omega$, which is called as _terminal
region_ , and a control Lyapunov function. The set $\Omega$ must be chosen
such that it is invariant for (1) controlled by a local auxiliary controller.
It can be determined off-line such that invariance property of $\Omega$ holds
and all the constraints are satisfied. The actual implementation of the NMPC
developed in Chen and Allgöwer (1998) would not use the local auxiliary
controller, but rather guarantee recursive feasibility of the closed-loop
system. As a contrast, the dual-mode NMPC developed in Michalska and Mayne
(1993) is actually a receding horizon implementation of a switching control
mechanism; that is, it applies the optimal control input from solving a
optimization problem over a finite horizon so as to drive the state to the
terminal region $\Omega$; once the state enters $\Omega$, it switches the
control input to the local feedback law so as to steer the state to the
origin. We highlight that the terminal region $\Omega$ is a positively
invariant set for (1) (see e.g. Bertsekas and Rhodes (1971)) and the analytic
properties of positively invariant sets have been widely studied, see e.g.
Blanchini and Miani (2008). The design method of the set $\Omega$ developed in
Chen and Allgöwer (1998) is based on the results presented in Michalska and
Mayne (1993), which is calculated off-line by linearizing the system (1). This
method is quite conservative, especially for systems with high nonlinearities.
The terminal region $\Omega$ is very important for stability analysis of the
closed-loop NMPC and affects the feasible region of NMPC. Therefore, this
paper considers the data-driven approximations for feasible regions and the
terminal regions of the nonlinear systems.

In this paper, we develop a data-driven learning framework for approximating
feasible region and invariant set of a perturbed nonlinear systems using
feasibility information of data samples. This work can be seen as an extension
of the results developed in Ong et al. (2006) and Chakrabarty et al. (2016,
2020). Here, we consider the continuous-time nonlinear system with bounded
disturbance and polyhedral state and input constraints. The main contributions
of this paper are as follows: (i) we develop a computation method to obtain
the invariant set for the nonlinear system given by (1); (ii) we develop a
data-driven learning method to obtain an inner and outer approximation of the
feasible region of the NMPC using support vector machine; (iii) we extend the
above two results to perturbed nonlinear system, where disturbance is bounded
and the state and input space are polyhedral.

_Notation_ : We denote $\mathbb{R}^{n}$ for the $n$-dimensional Euclidean
space and $\mathcal{C}^{1}(\mathbb{R}^{n})$ for the space of continuously
differentiable functions on $\mathbb{R}^{n}$. Besides, we let $\mathbb{R}$ and
$\mathbb{R}_{\geq 0}$ denote the real and the non-negative real numbers. For a
set $\mathbb{S}$, we let $\mathrm{Card}(\mathbb{S})$ denote the cardinality of
$\mathbb{S}$, $\mathrm{Vol}(\mathbb{S})$ denote the Lebesgue measure of
$\mathbb{S}$, and $\mathrm{Int}(\mathbb{S})$ denote the set of interior points
of $\mathbb{S}$. Furthermore, a continuous function
$\alpha:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ belongs to class
$\mathcal{K}$ if it is continuous, $\alpha(0)=0$, and it is strictly
increasing. It belongs to class $\mathcal{K}_{\infty}$ if it belongs to class
$\mathcal{K}$ and it is unbounded. Finally, given a point $x\in\mathbb{R}^{n}$
and a bounded and non-empty set $\mathcal{A}\subseteq\mathbb{R}^{n}$, the
distance of $x$ to the set $\mathcal{A}$ is defined as
$|x|_{\mathcal{A}}=\mathrm{dist}(x,\mathcal{A})\triangleq\mathrm{inf}\\{\|x-\xi\|_{2}:\xi\in\mathcal{A}\\}$.
A closed ball centralized at a given point $\xi_{0}\in\mathbb{R}^{n}$ with a
radius of $\rho\geq 0$ is denoted as
$\mathbb{B}^{n}(\xi_{0},\rho)\triangleq\\{\xi\in\mathbb{R}^{n}:\|\xi-\xi_{0}\|\leq\rho\\}$.
If $\xi_{0}=0$, we let $\mathbb{B}^{n}(\rho)$ denote $\mathbb{B}^{n}(0,\rho)$
as a shorthand.

##  2 Problem formulation

We extend (1) to a nonlinear perturbed system described as

$\displaystyle\dot{x}{(t)}$
$\displaystyle={F}\left(x(t),u(t),w(t)\right),\quad x(0)=\bar{x},\quad
t\in\mathbb{R}_{\geq 0},$ (3) $\displaystyle\quad
x(t)\in\mathbb{X}\subset\mathbb{R}^{n},\quad
u(t)\in\mathbb{U}\subset\mathbb{R}^{m},$ (4)

where ${x}\in\mathbb{R}^{n}$, $u\in\mathbb{R}^{m}$, and $w\in\mathbb{R}^{n}$
are the state vector, the control input, and the additive disturbance.
Besides, we confine the state space $\mathbb{X}$ and the input space
$\mathbb{U}$ to be polytopes. We give the following assumptions for (3):

###### Assumption 1

A1

The set $\mathbb{U}\in\mathbb{R}^{m}$ is compact, convex and contains
$0\in\mathbb{R}^{m}$ in its interior. The set $\mathbb{X}$ is a
$n$-dimensional hypercube and contains $0\in\mathbb{R}^{n}$ in its interior.

A2

The disturbance $w(\cdot)\in\mathbb{W}$, where
$\mathbb{W}=\left\\{w:|w|\leq\bar{w}\right\\}$.

A3

$F:\mathbb{X}\times\mathbb{U}\times\mathbb{W}\to\mathbb{X}$ is locally
Lipschitz in $x$ such that $\|F(x_{1},u,w)-F(x_{2},u,0)\|\leq
L_{F}\|x_{1}-x_{2}\|+\beta(|w|_{\infty})$ for all
${x_{1},x_{2}}\in\mathbb{X}$, ${u}\in\mathbb{U}$, and $w\in\mathbb{W}$, where
$L_{F}$ is local Lipschitz constant and $\beta(\cdot)$ is a $\mathcal{K}$
function.

The nominal model of the system (3) is given by

$\displaystyle\dot{\hat{x}}(t)$
$\displaystyle={F}\left(\hat{x}(t),u(t),0\right)\triangleq
f(\hat{x}(t),u(t)),$ $\displaystyle\quad\qquad\qquad\hat{x}(0)=\bar{x},\quad
t\in\mathbb{R}_{\geq 0},$ (5)

and the additive disturbance $w(t),t\geq 0$, in (3) satisfies

$\displaystyle w(t)$ $\displaystyle=\dot{x}(t)-f(x(t),u(t)).$ (6)

Assumption A2 indicates that $w(t)$ is bounded on $\mathbb{W}$. This kind of
model uncertainty has been studied in previous papers about robustness in MPC
as in Zou et al. (2019) and Lu et al. (2019).

Given a state measurement $x_{k}$, a prediction horizon $T$, a state cost
$\Psi(\cdot,\cdot)$, a terminal penalty $\Phi(\cdot)$, and a terminal set
$\Omega$, an optimization problem of NMPC for (3) is defined as follows:

$\displaystyle\hat{u}^{\star}=\underset{u}{\arg\min}~{}\mathcal{J}(x_{k},\hat{u}(\cdot)),$
(7a) with
$\displaystyle\mathcal{J}(x_{k},\hat{u}(\cdot))\buildrel\vartriangle\over{=}\int_{0}^{T}\Psi(\hat{x}(\tau),\hat{u}(\tau))\mathrm{d}\tau+\Phi(\hat{x}(T)),$
(7b) subject to
$\displaystyle\dot{\hat{x}}(\tau)=f(\hat{x}(\tau),\hat{u}(\tau)),\quad
0\leq\tau\leq T,$ (7c)
$\displaystyle\hat{x}(\tau)\in\mathbb{X},\quad\hat{u}(\tau)\in\mathbb{U},\quad
0\leq\tau\leq T,$ (7d)
$\displaystyle\hat{x}(0)=x_{k},\quad\hat{x}(T)\in\Omega.$ (7e)

For (7), we give the following two assumptions.

###### Assumption 2

The stage cost $\Psi(x,u)$ in (7b) satisfies $\Psi(0,0)=0$ and
$\Psi(x,u)\geq\alpha_{\Psi}(\|x\|)+\rho_{\Psi}$, where
$\rho_{\Psi}\in\mathbb{R}$ and $\alpha_{\Psi}(\cdot)$ is a
$\mathcal{K}_{\infty}$ function. Moreover, $\Psi(x,u)$ is locally Lipschitz in
$x$ in $\mathbb{X}\times\mathbb{U}$ with constant $L_{\Psi}$; that is,
$\|\Psi(x_{1},u)-\Psi(x_{2},u)\|\leq L_{\Psi}\|x_{1}-x_{2}\|$ for all
$x_{1},x_{2}\in\mathbb{X}$ and all $u\in\mathbb{U}$.

###### Assumption 3

The terminal cost $\Phi(x)$ in (7b) satisfies

  1. 1.

$\underline{\alpha}_{\Phi}(\|x\|)\leq\Phi(x)\leq\overline{\alpha}_{\Phi}(\|x\|)$,
where $\underline{\alpha}_{\Phi}(\cdot)$ and $\overline{\alpha}_{\Phi}(\cdot)$
are $\mathcal{K}_{\infty}$ functions;

  2. 2.

$\Phi(f(x,\kappa(x)))-\Phi(x)\leq-\Psi(x,\kappa(x))$ for all $x\in\Omega$,
where $\kappa(\cdot):\mathbb{R}^{n}\to\mathbb{R}^{m}$ is continuous and
$\kappa(0)=0$;

  3. 3.

$\Phi(\cdot)$ is locally Lipschitz in $\Omega$ with a Lipschitz constant
$L_{\Phi}$.

Assumptions 2 and 3 are fairly standard assumptions to guarantee feasibility
and stability for nonlinear systems under MPC framework, see e.g. Allgöwer and
Zheng (2012); Grüne and Pannek (2017). Specifically, several numerical methods
to obtain $\Omega$ and the local auxiliary controller $\kappa(x)$ satisfying
Assumption 3 were proposed in Michalska and Mayne (1993); Chen and Allgöwer
(1998). The feature of $\kappa(x)$ makes it desirable to construct the
feasible control input and further to prove the closed-loop stability.
Furthermore, several ways to design $\Psi(\cdot,\cdot)$ and $\Phi(\cdot)$
satisfying Assumptions 1 and 2 were proposed in Hashimoto et al. (2016). In
practice, $\Psi(x,u)$ and $\Phi(x)$ can be taken as quadratic functions.

Now, we assume that for (7), Assumptions 1 and 2 hold. Then we introduce the
following definitions.

######  Definition 4 (Feasible State)

A sampling state $x_{k}\in\mathbb{X}$ is feasible if there exists a admissible
control function $\kappa(\cdot)\in\mathbb{U}$, $t\geq 0$, satisfying all the
conditions (7c)–(7e) except the minimization of (7b).

Using Definition (4), we define the feasible region $\mathbb{F}$:

$\mathbb{F}\buildrel\vartriangle\over{=}\big{\\{}x_{k}\in\mathbb{X}:x_{k}~{}\mathrm{is}~{}\mathrm{feasible}~{}\mathrm{for}~{}(7)\big{\\}}.$
(8)

Moreover, we give another two definitions for the $\mathbb{F}$ in (8).

######  Definition 5 (Full Region)

The feasible region $\mathbb{F}$ is full if
$\mathbb{F}=\mathbb{X}\subset\mathbb{R}^{n}$.

######  Definition 6 (Globally Full Region)

The feasible region $\mathbb{F}$ is globally full if
$\mathbb{F}=\mathbb{X}=\mathbb{R}^{n}$.

Ideally we would want that the feasible region $\mathbb{F}$ is a globally full
region. But it is impossible for both theoretical and application aspects. The
best case is that the feasible region $\mathbb{F}$ is a full region such that,
we can execute the receding horizon NMPC all over the state space $\mathbb{X}$
with guarantee of the closed-loop stability. Therefore, the main objective of
this paper is to learn the feasible region of a nonlinear system under the MPC
framework, as well as the nonlinear perturbed case. Our method provides
freedom in the design of the feasible region needed for closed-loop stability
and does not need accurate information of the additively perturbed nonlinear
system. Furthermore, we give some considerations for enlarging the feasible
region so as to make it to be a full region.

##  3 Basic Results

In this section, we use the support vector machine learning method to estimate
the feasible region for the nonlinear system given by (2) based on the
feasibility information of a low-discrepancy sequence. Specifically, using
deterministic sampling data, the SVM learning method is employed to learn the
boundary function of the feasible region.

To that end, we let $\partial\mathbb{F}$ denote the boundary of the feasible
region $\mathbb{F}$ given by (8). Then we make the following assumption.

###### Assumption 7

There exists a continuous function $\phi(\cdot)\in\mathcal{C}^{1}(\mathbb{X})$
such that, the feasible region $\mathbb{F}$ can be represented as the zero
superlevel-set of $\phi(\cdot)$; that is

$\displaystyle\mathrm{Int}(\mathbb{F})=\bigl{\\{}x\in\mathbb{X}:\phi(x)>0\bigr{\\}}.$
(9)

Now using Assumption 7, to characterize the feasible region $\mathbb{F}$, we
just need to construct the function $\phi(\cdot)$.

###  3.1 Point-wise sampling data and SVM learning

To construct the function $\phi(\cdot)$, we need to sample the state space
$\mathbb{X}$ and obtain the point-wise sampling data. To that end, we borrow
the standard notion of low-discrepancy sequences from Niederreiter (1988) and
construct a data set on a multilevel sparse-grid.

######  Definition 8 (Chakrabarty et al. (2016))

The discrepancy of a sequence $\\{x_{k}\\}_{k=1}^{N}\subset\mathbb{X}$ is
defined as

$\displaystyle\mathscr{D}_{N}\left(\\{x_{k}\\}_{k=1}^{N}\right)\triangleq\underset{X\in\mathscr{T}}{\mathrm{sup}}\bigg{|}\frac{\mathrm{Card}_{\mathbb{X}}(\mathbb{S})}{N}-\frac{\mathrm{Vol}(\mathbb{S})}{\mathrm{Vol}(\mathbb{X})}\bigg{|},$

where $\mathscr{T}$ is a set of state dimensional intervals and defined as
$\mathscr{T}\buildrel\vartriangle\over{=}\prod_{k=1}^{M}\bigl{[}a_{k},b_{k}\bigr{)}=\bigl{\\{}x_{k}\in\mathbb{X}:a_{k}\leq
x_{k}<b_{k}\bigr{\\}}\subset\mathbb{X}$, and
$M=\mathrm{Card}_{\mathbb{X}}(\mathbb{S})\triangleq\mathrm{Card}\bigl{\\{}k\in\left\\{1,\ldots,N\right\\}:x_{k}\in\mathbb{S}\bigr{\\}}$.
For a low-discrepancy sequence $\\{x_{k}\\}_{k=1}^{N}$,
$\underset{N\to\infty}{\mathrm{lim}}\mathscr{D}_{N}\bigl{(}\left\\{x_{k}\right\\}_{k=1}^{N}\bigr{)}=0$.

Over the state space $\mathbb{X}$, we generate a point-wise data set

$\displaystyle\varPi\buildrel\vartriangle\over{=}\left\\{x_{k}\right\\}_{k=1}^{N}=\bigl{\\{}x_{1},x_{2},\ldots,x_{N}\bigr{\\}},$
(10)

and assume that it satisfies the following assumption.

###### Assumption 9

The sequence $\left\\{x_{k}\right\\}_{k=1}^{N}$ is a low-discrepancy sequence
on $\mathbb{X}$.

Now, using Assumption 9, we describe the training data in $\varPi$ following
the results presented in Chakrabarty et al. (2016). For each point $x_{k}$,
$k=1,\ldots,N$ in $\varPi$, we solve the optimization problem (7) over a
finite-time horizon $T$. If there is a feasible solution for $x_{k}$, then we
label it as $\mathcal{O}(x_{k})=+1$; otherwise, we label it as
$\mathcal{O}(x_{k})=-1$. Then we have the following two cases:

$\displaystyle\mathcal{O}(x_{k})=\begin{cases}+1,\quad x_{k}\in\mathbb{F};\\\
-1,\quad x_{k}\notin\mathbb{F},\\\ \end{cases}$ (11)

for each point $x_{k}$ in $\varPi$. Note that this is a two-class pattern
recognition problem, which divides the data set $\varPi$ into two classes,
that is,

$\displaystyle\varPi^{+}$
$\displaystyle=\big{\\{}x_{k}\in\varPi:~{}\mathcal{O}(x_{k})=+1\big{\\}};$
$\displaystyle\varPi^{-}$
$\displaystyle=\big{\\{}x_{k}\in\varPi:~{}\mathcal{O}(x_{k})=-1\big{\\}}.$

This problem can be solved using the SVM bi-classifier. Here, we use the
kernel-based SVM classification method, see e.g. Vapnik (2013); Burges (1998).
Now, we can construct a decision function
$\phi_{N}(\cdot)\in\mathcal{C}^{1}(\mathbb{X})$,

$\displaystyle\phi_{N}(x):\varPi\to y_{k}\in\\{+1,-1\\},$ (12)

that accurately classifies an arbitrary state in $\varPi$ as feasible or
infeasible.

Next, to facilitate the classification, we introduce a map:
$\Gamma:\mathbb{X}\to\mathbb{H}$, which maps the data in $\varPi$ to a higher-
dimensional Hilbert space $\mathbb{H}$ where the data is linearly separable,
see e.g. Burges (1998). In this case, we can let
$\phi_{N}(x)=\rho\Gamma(x)+v$, where $\rho$ and $v$ are vector parameters that
determine the orientation of the separating hyperplane. Referring to Ong et
al. (2006), our problem can be formulated as the constrained optimization
problem

$\displaystyle\rho^{\star}=\underset{\rho}{\arg\min}~{}\frac{1}{2}\rho^{\mathrm{T}}\rho+L\sum_{k=1}^{N}s_{k},$
(13a) subject to $\displaystyle\rho\Gamma(x_{k})+v\geq 1-s_{k},\quad s_{k}\geq
0,\quad\forall k\in\varPi^{+},$ (13b)
$\displaystyle\rho\Gamma(x_{k})+v\leq-1+s_{k},\quad\forall k\in\varPi^{-},$
(13c)

where $L$ is a regularization parameter and $s_{k}$ is a slack variable to
relax separability constraints.

The optimization problem (13) is convex and can be solved via its dual
formulation. Letting $\alpha_{i}$ and $\alpha_{j}$ denote the Lagrange
multiplier vector for (13b) and (13c), (13) can reformulated as the dual
problem

$\displaystyle\alpha^{\star}=\underset{\alpha}{\mathrm{arg~{}min}}$
$\displaystyle~{}\frac{1}{2}\sum_{i=1}^{N}\sum_{k=1}^{N}\alpha_{i}\alpha_{k}y_{i}y_{k}\Gamma(x_{i})\Gamma(x_{k})-\sum_{k=1}^{N}\alpha_{k},$
(14a) $\displaystyle\mathrm{subject\;to}~{}~{}$
$\displaystyle\sum_{k=1}^{N}\alpha_{k}y_{k}=0,$ (14b) $\displaystyle
0\leq\alpha_{k}\leq L,\quad k=1,\ldots,N.$ (14c)

After solving (14), we obtain the SVM decision function as

$\phi_{N}(x)=\sum_{k=1}^{N}\alpha_{k}^{\star}y_{k}\Gamma(x_{k})\Gamma(x)+v,$
(15)

and the estimated feasibility region boundary is given by

$\mathcal{S}_{N}\triangleq\left\\{x\in\mathbb{X}:\phi_{N}(x)=0\right\\}.$ (16)

To improve accuracy of the SVM learning when effective data in $\varPi^{+}$
are not enough, we iteratively select new points satisfying
$\mathcal{O}(x_{k})=+1$ and expand into the data set $\varPi^{+}$ to ensure
richness of data. This ensures that a non-trivial decision function (15)
exists. Furthermore, for the inner product $\Gamma(x_{k})\Gamma(x)$ in (15),
it can be replaced by using the kernel function. A common choice of kernel
function is the Gaussian kernel:

$\Gamma(x_{k})\Gamma(x)=\frac{{\|x_{k}-x\|}^{2}}{2\sigma^{2}},$

where $\sigma^{2}$ is the kernel variance, see e.g. Steinwart (2001).

###  3.2 Inner and outer approximations of $\partial\mathbb{F}$

The outer and inner approximations of a bounded region were introduced in
Darup and Mönnigmann (2012); Deffuant et al. (2007). Using (9), we define the
strict outer ($\mathbb{F}^{-}$) and inner ($\mathbb{F}^{+}$) region of the
feasible region $\mathbb{F}$ as

$\displaystyle\mathbb{F}^{+}$
$\displaystyle\triangleq\left\\{x\in\mathbb{X}:\phi(x)>\varepsilon\right\\};$
(17) $\displaystyle\mathbb{F}^{-}$
$\displaystyle\triangleq\left\\{x\in\mathbb{X}:\phi(x)<-\varepsilon\right\\}.$
(18)

for some $\varepsilon>0$. In this subsection, we construct strict outer and
inner approximations of the feasible region $\mathbb{F}$, which ensures that
they contain no feasible (infeasible) samples.

Now, based on the super-level sets of $\phi_{N}(x)$, we can find a positive
constant $\varepsilon^{+}>0$ which defines the inner approximation
$\mathbb{F}_{N}^{+}$ by

$\displaystyle\mathbb{F}_{N}^{+}$
$\displaystyle\triangleq\bigl{\\{}x\in\mathbb{X}:\phi_{N}(x)>\varepsilon^{+}\bigr{\\}}.$
(19)

Using the data in $\varPi^{+}$, we determine the parameter $\varepsilon^{+}$
by solving the optimization problem

$\displaystyle\varepsilon^{+}$
$\displaystyle=\underset{\varepsilon>0}{\mathrm{arg~{}min}}~{}\varepsilon,$
(20a) $\displaystyle\mathrm{subject\;to}~{}$
$\displaystyle\phi_{N}(x_{k})>\varepsilon,\quad\forall x_{k}\in\varPi^{+}.$
(20b)

Analogously, there exists a negative constant $\varepsilon^{-}<0$ which
defines the outer approximation $\mathbb{F}_{N}^{-}$ by

$\displaystyle\mathbb{F}_{N}^{-}$
$\displaystyle\triangleq\bigl{\\{}x\in\mathbb{X}:\phi_{N}(x)<\varepsilon^{-}\bigr{\\}}.$
(21)

This constant $\varepsilon^{-}$ can be determined by using the data in
$\varPi^{-}$ and solving the optimization problem

$\displaystyle\varepsilon^{-}$
$\displaystyle=\underset{\varepsilon<0}{\mathrm{arg~{}max}}~{}\varepsilon,$
(22a) $\displaystyle\mathrm{subject\;to}~{}$
$\displaystyle\phi_{N}(x_{k})<\varepsilon,\quad\forall x_{k}\in\varPi^{-}.$
(22b)

We highlight that novel algorithms for constructing strict inner and outer
approximations of the feasible region $\mathbb{F}$ are developed by solving
(20) and (22). Besides, the inner approximation $\mathbb{F}_{N}^{+}$ is based
only on the feasible samples in $\varPi^{+}$ and the outer approximation
$\mathbb{F}_{N}^{-}$ is based only on the infeasible samples in $\varPi^{-}$.

###  3.3 Convergence analysis

In this subsection, we provide convergence analysis for the approximation
error. We show that as the number of samples $N$ from the low-discrepancy
sequence increases, the inner and outer approximations with universal kernels
will converge to a strict approximation of the actual feasible region
$\mathbb{F}$.

First, we introduce definition of universal kernels adopted from Steinwart
(2001) and give a proposition.

###### Definition 10

A continuous kernel $\mathcal{K}(\cdot,\cdot)$ is universal if the space of
all functions induced by $\mathcal{K}(\cdot,\cdot)$ is dense in
$\mathcal{C}(\mathbb{X})$.

###### Proposition 11

Every universal kernel separates all compact subsets of $\mathbb{X}$.

Definition 10 indicates that for each continuous function $\phi$ and
$\varepsilon>0$, there exists a function $\phi_{\mathcal{K}}$ induced by a
universal kernel $\mathcal{K}$ such that
$|\phi-\phi_{\mathcal{K}}|<\varepsilon$. Proposition 11 indicates that the
kernel-based SVM classification method can be adopted to solve the feasibility
boundary learning problem. In this case, if we design

$\displaystyle\mathcal{K}(x_{k},x)=\Gamma(x_{k})\Gamma(x),$ (23)

then we can accurately separate the pairwise disjoint compact sets
$\varPi^{+}$ and $\varPi^{-}$ in $\mathbb{X}$.

We highlight that there is not computational method for determining such a
separating function in Steinwart (2001). By modifying the arguments presented
in Steinwart (2001), we collect the training data from a low-discrepancy
sequence. It guarantees that increasing the number of samples $N$ for the SVM
with a universal kernel results in an increasing accuracy of the feasible
region learning, and reduces the conservativeness of the inner and outer
approximations of $\mathbb{F}$. We summarize our discussions as the following
theorem, which is slightly similar to the results given in Chakrabarty et al.
(2016).

###### Theorem 12

Consider the system (2), suppose Assumptions 1, 2, 3, 7, and 9 hold, and let
$\mathcal{K}(\cdot,\cdot)$ be a universal kernel on $\mathbb{X}$ in (23) and
$\phi_{N,\mathcal{K}}(x)$ be the kernel-based SVM decision function from (15).
Then for any $L>0$ in the optimization problem (14) and for two compact
subsets $\mathbb{F}^{+}\subset\mathbb{X}$ and
$\mathbb{F}^{-}\subset\mathbb{X}$ in (17) and (18),

$\displaystyle\underset{N\to\infty}{\mathrm{lim}}\phi_{N}(x)=\phi(x),$ (24)
$\displaystyle\underset{N\to\infty}{\mathrm{lim}}\mathbb{F}_{N}^{+}=\mathbb{F}^{+},\quad\underset{N\to\infty}{\mathrm{lim}}\mathbb{F}_{N}^{-}=\mathbb{F}^{-},$
(25)
$\displaystyle\underset{N\to\infty}{\mathrm{lim}}\mathcal{S}_{N}=\partial\mathbb{F}.$
(26)

Theorem 12 indicates that there exists a sufficiently large number of samples
$N_{0}$ such that if $N\geq N_{0}$,

$\begin{cases}\phi(x)>0,&\quad\mathrm{for}\;x\in\mathbb{F}^{+};\\\
\phi(x)<0,&\quad\mathrm{for}\;x\in\mathbb{F}^{-},\end{cases}$

and the strict inner and outer approximations of the feasible region
$\mathbb{F}$ can be obtained by sampling $\mathbb{X}$ using low-discrepancy
sequences. Therefore, the SVM classifier using universal kernels is able to
approximate our feasible region boundaries with sufficiently high accuracy.

##  4 Main Results

In Section 3, we present a method to estimate the feasible region where the
states of the nominal system (1) can be driven to the terminal set $\Omega$
while satisfying the constraints (4). In this section, we extend this result
to the perturbed nonlinear system given by (3). However, the region
$\mathbb{F}$ obtaining from Section 3 may not be feasible for the perturbed
system when there exists unknown disturbance.

Using Assumption A2, (3), and (2), it follows that

$\displaystyle|F(x(t),u(t),w(t))-f(x(t),u(t))|\leq\bar{w}.$ (27)

So for the perturbed system (3), the set $\Omega$ may not be positively
invariant with the local feedback law $u(t)$. In this case, we introduce the
robust control invariant (RCI) set for (3) following the results presented in
Sun et al. (2018).

######  Definition 13 (RCI set)

Consider the perturbed system (3) with constraints $x\in\mathbb{X}$,
$u\in\mathbb{U}$ and $w\in\mathbb{W}$. The set
$\bar{\Omega}\subset\mathbb{X}\subseteq\mathbb{R}^{n}$ is an RCI set, if there
exists an ancillary controller $\kappa(x)$ with
$\bar{u}=u+\kappa(x)\in\mathbb{U}$ such that $\forall\bar{x}\in\bar{\Omega}$
and $\forall w\in\mathbb{W}$, $x(t)\in\bar{\Omega}$.

For the design of the RCI set $\bar{\Omega}$ (terminal region) for (3), the
set-theoretic methods and tube-based MPC scheme are very useful, see e.g.
Blanchini and Miani (2008) and Rakovic et al. (2012). However, the method
developed in Rakovic et al. (2012) can not be applied to our problem, since
the perturbed system given by (3) is continuous-time and nonlinear case; the
authors in Ong et al. (2006) and Chakrabarty et al. (2016) only consider the
nonlinear systems, not for the perturbed system.

###  4.1 RCI set design based on set-theoretic method

Recalling that the disturbance in (6) is actually an additive disturbance and
satisfies

$\displaystyle w(t)=F(x(t),u(t),w(t))-f(x(t),u(t)),~{}t\in\mathbb{R}_{\geq
0}.$ (28)

Besides, we note that the state and control do not depend on the disturbance
$w$.

Now, using Assumption A2, we introduce the following set

$\displaystyle\Omega_{0}\triangleq\bigl{\\{}x\in\Omega:|x|_{\partial\Omega}\geq\bar{w}\bigr{\\}},$
(29)

where $\Omega$ is a terminal region for (2) and $\partial\Omega$ is the
boundary of the set $\Omega$. Then for (29), we have the following lemma.

###### Lemma 14

Consider the set (29) and assume that $\Omega\supset\mathbb{R}^{n}(\bar{w})$
holds. Then, $\Omega_{0}$ is a RCI set for the perturbed system (3).

Note that $\Omega$ is a control invariant set for the nonlinear system (2) and
$\Omega_{0}\subset\Omega$. Since $\Omega\supset\mathbb{R}^{n}(\bar{w})$,
$\Omega\setminus\Omega_{0}\neq\varnothing$. Then using (29), we have that for
any $x\in\Omega_{0}$,
$|x+w|_{\partial\Omega_{0}}\geq|x|_{\partial\Omega}-\bar{w}\geq 0$ for all
$w\in\mathbb{W}$. That is, for any $x\in\Omega_{0}$, $x+w\in\Omega$ for all
$w\in\mathbb{W}$. According to Definition 13, $\Omega_{0}$ is a RCI set for
(3).

Lemma 14 indicates that the RCI set $\bar{\Omega}$ for (3) can be obtained
using the terminal set $\Omega$ for (2), while the set $\Omega$ can be
computed using the results developed in Chen and Allgöwer (1998) for perturbed
nonlinear systems.

###  4.2 Feasible region estimation for robust NMPC

In this subsection, we develop a feasible region estimation method for (3)
under the NMPC scheme using the feasibility information of the low-discrepancy
data samples.

First, using the RCI set $\Omega_{0}$ in (29) for the perturbed system given
by (3), we modify the problem (7) to a new optimization problem of robust NMPC
which is given by

$\displaystyle\hat{u}^{\star}=\underset{u}{\arg\min}~{}\mathcal{J}(x_{k},\hat{u}(\cdot)),$
(30a) subject to
$\displaystyle\dot{\hat{x}}(\tau)=f(\hat{x}(\tau),\hat{u}(\tau)),\quad
0\leq\tau\leq T,$ (30b)
$\displaystyle\hat{x}(\tau)\in\mathbb{X},\quad\hat{u}(\tau)\in\mathbb{U},\quad
0\leq\tau\leq T,$ (30c)
$\displaystyle\hat{x}(0)=x_{k},\quad\hat{x}(T)\in\Omega_{0},$ (30d)

where $\mathcal{J}(x_{k},\hat{u}(\cdot))$ is given by (7b) and
$\Psi(\cdot,\cdot)$ and $\Phi(\cdot)$ satisfy Assumptions 1 and 3.

Next, using the same procedure developed in Section 3, we generate the data
set and construct an approximate function of the boundary for the feasible
region. Then similar to Theorem 12, we have the following theorem which shows
the approximation of the feasible region $\bar{\mathbb{F}}$ for the perturbed
systems in (3).

###### Theorem 15

Consider the system (3), suppose Assumptions 1, 3, 7, and 9 hold, and let
$\bar{\Omega}$ be the maximal RCI set in $\mathbb{X}$ and $\bar{\mathbb{F}}$
be the maximal feasible region for (3). Then, $\Omega_{0}\subset\bar{\Omega}$
and

$\displaystyle\bar{\mathbb{F}}_{0}\triangleq\bigl{\\{}x\in{\mathbb{F}}:|x|_{\partial{\mathbb{F}}}\geq\bar{w}\bigr{\\}}\subset\bar{\mathbb{F}}.$
(31)

Theorem 15 indicates that based on the feasible region $\mathbb{F}$ for (2),
we can obtain the conservative terminal region $\Omega_{0}$ and the
conservative feasible region $\bar{\mathbb{F}}_{0}$ for the perturbed systems.
Since the disturbance $w$ in (28) is unknown, it is very difficult to compute
the strict feasible region for the perturbed system (2). Therefore, we have
adopted the low-discrepancy sequences, optimal kernel-based SVM
classification, and set-theoretic methods to increase the accuracy of the
data-driven feasible region estimation. Although the feasible region
$\bar{\mathbb{F}}$ can be enlarged with a longer prediction horizon, it needs
more computation resources for solving the optimization problems (30).

##  5 illustrative numerical example

To illustrate the key ideas presented in this paper, we consider a perturbed
nonlinear system representing the cart-and-spring system adopted from Magni et
al. (2003); Li and Shi (2014) given by

$\begin{cases}\dot{x}_{1}(t)=x_{2}(t),\\\
\dot{x}_{2}(t)=-\frac{k_{0}}{M}e^{-x_{1}(t)}x_{1}(t)-\frac{h_{d}}{M}x_{2}(t)+\frac{u(t)}{M}+w(t),\end{cases}$

where $x=\left[x_{1};x_{2}\right]\in\mathbb{R}^{2}$, $u\in\mathbb{R}$, and
$w\in\mathbb{R}$. The parameters are given as: $M=1.8~{}\mathrm{kg}$,
$k_{0}=1.2~{}\mathrm{N/m}$, $h_{d}=0.25~{}\mathrm{N\cdot s/m}$, The
constrained state and control space associated with the system is given by
$\mathbb{X}=\\{x\in\mathbb{R}^{2}:|x|\leq 2\\}$ and
$\mathbb{U}=\\{u\in\mathbb{R}:|u|\leq 3\\}$ and let
$w(t)\in\mathbb{B}^{1}(0.01)$. To satisfy Assumptions 1 and 3, we set the
stage cost as $\Psi(x,u)=\|x\|_{Q}^{2}+\|u\|_{R}^{2}$ with $Q=I_{2}$ and
$R=0.1$ and the terminal cost as $\Phi(x)=\|x\|_{P}^{2}$ with
$P=\left[5.0511,-2.2731;-2.2731,2.4586\right]$. The local controller is
designed as $\kappa(x)=-Kx$ with $K=\left[4.2291~{}4.8551\right]$ using the
procedure presented in Chen and Allgöwer (1998).

Next, the terminal region $\Omega$ is parametrized by
$\Omega(\mu)=\bigl{\\{}x\in\mathbb{R}^{2}:x^{\mathrm{T}}Px\leq\mu\bigr{\\}}$,
where $\mu\in(0,1)$. Then, we use a kernel-based SVM with $\sigma=0.8$ for the
feasible region approximation. Using the feasibility information of data
samples, approximations of the feasible region boundaries with different
prediction horizons $T$ and different $\mu$ are illustrated in Figures. 1–3.
Furthermore, it can be seen from Figures 1 and 2 that larger terminal region
results in larger feasible region. From Figures 1 and 3, we can see that
longer prediction horizon results in larger feasible region.

Figure 1: The boundary of the feasible region $\bar{\mathbb{F}}_{0}$ resulting
from the optimization (30) with the terminal region $\Omega_{0}(0.5)$ and the
prediction horizon $T=1.0~{}s$ in red. The infeasible samples in $\varPi^{-}$
are shown in black; the feasible samples in $\varPi^{+}$ are shown in blue;
the feasible samples selected using (16) are shown in green; the boundary of
the terminal region $\Omega_{0}(0.5)$ is shown in magenta.  Figure 2: The
boundary of the feasible region $\bar{\mathbb{F}}_{0}$ resulting from the
optimization (30) with the terminal region $\Omega_{0}(0.9)$ and the
prediction horizon $T=1.0~{}s$ in red. The other notations are the same as
Fig. 1.  Figure 3: The boundary of the feasible region $\bar{\mathbb{F}}_{0}$
resulting from the optimization (30) with the terminal region
$\Omega_{0}(0.5)$ and the prediction horizon $T=2.0~{}s$ in red. The other
notations are the same as Fig. 1.

##  6 Conclusions and future work

This paper develops a data-driven learning framework for designing the
terminal region and approximating the feasible region offline using the
feasibility information of low-discrepancy data samples and support vector
machine learning. Our approach provides the freedom in the design of the
feasible region needed for closed-loop stability, when disturbance is subject
to deterministic bounds. Future work will combine the offline feasible region
learning with online NMPC design.

## References

  * Allgöwer and Zheng (2012) Allgöwer, F. and Zheng, A. (2012).  _Nonlinear Model Predictive Control_ , volume 26.  Birkhäuser.
  * Bertsekas and Rhodes (1971) Bertsekas, D. and Rhodes, I. (1971).  On the minimax reachability of target sets and target tubes.  _Automatica_ , 7, 233–247.
  * Blanchini and Miani (2008) Blanchini, F. and Miani, S. (2008).  _Set-Theoretic Methods in Control_.  Springer, London.
  * Burges (1998) Burges, C.J. (1998).  A tutorial on support vector machines for pattern recognition.  _Data mining and knowledge discovery_ , 2(2), 121–167.
  * Chakrabarty et al. (2020) Chakrabarty, A., Danielson, C., Di Cairano, S., and Raghunathan, A. (2020).  Active learning for estimating reachable sets for systems with unknown dynamics.  _IEEE Transactions on Cybernetics_.
  * Chakrabarty et al. (2016) Chakrabarty, A., Dinh, V., Corless, M.J., Rundell, A.E., Żak, S.H., and Buzzard, G.T. (2016).  Support vector machine informed explicit nonlinear model predictive control using low-discrepancy sequences.  _IEEE Transactions on Automatic Control_ , 62(1), 135–148.
  * Chen and Allgöwer (1998) Chen, H. and Allgöwer, F. (1998).  A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability.  _Automatica_ , 34(10), 1205–1217.
  * Darup and Mönnigmann (2012) Darup, M.S. and Mönnigmann, M. (2012).  Low complexity suboptimal explicit nmpc.  _IFAC Proceedings Volumes_ , 45(17), 406–411.
  * Deffuant et al. (2007) Deffuant, G., Chapel, L., and Martin, S. (2007).  Approximating viability kernels with support vector machines.  _IEEE Transactions on Automatic Control_ , 52(5), 933–937.
  * Grüne and Pannek (2017) Grüne, L. and Pannek, J. (2017).  _Nonlinear Model Predictive Control_.  Springer.
  * Hashimoto et al. (2016) Hashimoto, K., Adachi, S., and Dimarogonas, D.V. (2016).  Self-triggered model predictive control for nonlinear input-affine dynamical systems via adaptive control samples selection.  _IEEE Transactions on Automatic Control_ , 62(1), 177–189.
  * Li and Shi (2014) Li, H. and Shi, Y. (2014).  Distributed receding horizon control of large-scale nonlinear systems: Handling communication delays and disturbances.  _Automatica_ , 50(4), 1264–1271.
  * Lu et al. (2019) Lu, J., Xi, Y., and Li, D. (2019).  Stochastic model predictive control for probabilistically constrained markovian jump linear systems with additive disturbance.  _International Journal of Robust and Nonlinear Control_ , 29(15), 5002–5016.
  * Magni et al. (2003) Magni, L., De Nicolao, G., Scattolini, R., and Allgöwer, F. (2003).  Robust model predictive control for nonlinear discrete-time systems.  _International Journal of Robust and Nonlinear Control_ , 13(3-4), 229–246.
  * Michalska and Mayne (1993) Michalska, H. and Mayne, D. (1993).  Robust receding horizon control of constrained nonlinear systems.  _IEEE Transactions on Automatic Control_ , 38(11), 1623–1633.
  * Niederreiter (1988) Niederreiter, H. (1988).  Low-discrepancy and low-dispersion sequences.  _Journal of Number Theory_ , 30(1), 51–70.
  * Ong et al. (2006) Ong, C.J., Sui, D., and Gilbert, E.G. (2006).  Enlarging the terminal region of nonlinear model predictive control using the support vector machine method.  _Automatica_ , 42(6), 1011–1016.
  * Rakovic et al. (2012) Rakovic, S., Kouvaritakis, B., Cannon, M., Panos, C., and Findeisen, R. (2012).  Parameterized tube model predictive control.  _IEEE Transactions on Automatic Control_ , 57(11), 2746–2761.
  * Rawlings and Mayne (2009) Rawlings, J.B. and Mayne, D.Q. (2009).  _Model Predictive Control: Theory and Design_.  Nob Hill Pub.
  * Steinwart (2001) Steinwart, I. (2001).  On the influence of the kernel on the consistency of support vector machines.  _Journal of machine learning research_ , 2(Nov), 67–93.
  * Sun et al. (2018) Sun, T., Pan, Y., Zhang, J., and Yu, H. (2018).  Robust model predictive control for constrained continuous-time nonlinear systems.  _International Journal of Control_ , 91(2), 359–368.
  * Vapnik (2013) Vapnik, V. (2013).  _The nature of statistical learning theory_.  Springer science & business media.
  * Xi and Li (2019) Xi, Y. and Li, D. (2019).  _Predictive Control: Fundamentals and Developments_.  John Wiley & Sons.
  * Zou et al. (2019) Zou, Y., Su, X., Li, S., Niu, Y., and Li, D. (2019).  Event-triggered distributed predictive control for asynchronous coordination of multi-agent systems.  _Automatica_ , 99, 92–98.