# Approximate Dynamic Programming for Constrained Piecewise Affine Systems
with Stability and Safety Guarantees

Kanghui He, Shengling Shi, Ton van den Boom, and Bart De Schutter This paper
is part of a project that has received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (Grant agreement No. 101018826 - CLariNet).All authors are
affiliated with Delft Center for Systems and Control, Delft University of
Technology, Delft, The Netherlands (e-mail: {k.he, s.shi-3, a.j.j.vandenboom,
b.deschutter}@tudelft.nl).

###### Abstract

Infinite-horizon optimal control of constrained piecewise affine (PWA) systems
has been approximately addressed by hybrid model predictive control (MPC),
which, however, has computational limitations, both in offline design and
online implementation. In this paper, we consider an alternative approach
based on approximate dynamic programming (ADP), an important class of methods
in reinforcement learning. We accommodate non-convex union-of-polyhedra state
constraints and linear input constraints into ADP by designing PWA penalty
functions. PWA function approximation is used, which allows for a mixed-
integer encoding to implement ADP. The main advantage of the proposed ADP
method is its online computational efficiency. Particularly, we propose two
control policies, which lead to solving a smaller-scale mixed-integer linear
program than conventional hybrid MPC, or a single convex quadratic program,
depending on whether the policy is implicitly determined online or explicitly
computed offline. We characterize the stability and safety properties of the
closed-loop systems, as well as the sub-optimality of the proposed policies,
by quantifying the approximation errors of value functions and policies. We
also develop an offline mixed-integer linear programming-based method to
certify the reliability of the proposed method. Simulation results on an
inverted pendulum with elastic walls and on an adaptive cruise control problem
validate the control performance in terms of constraint satisfaction and CPU
time.

##  I INTRODUCTION

###  I-A Backgrounds

There has been an increasing interest in control of piecewise affine (PWA)
systems due to their capability of representing hybrid models and
approximating nonlinear dynamics [1]. In PWA systems, the state and input
spaces are jointly partitioned into polyhedral regions where an affine
subsystem is defined in each region. Control of PWA systems has a broad range
of applications, including robots with contact phenomena [2] and
transportation networks [3]. For PWA systems with constraints, model
predictive control (MPC) is widely applied. However, MPC for PWA systems still
faces challenges in computational complexity [4], because it involves solving
a mixed-integer linear programming (MILP) problem, which is NP-hard in general
and the complexity of which grows rapidly with the number of subsystems.
Explicit MPC [4], an offline version of MPC, requires solving a parametric
MILP problem, which is also suffering from computational complexity issues.
These issues make MPC only suitable for slow PWA processes or for low-
dimensional problems with short prediction horizons [2].

In contrast to MPC, reinforcement learning (RL) can learn a policy that
minimizes a finite-/infinite-horizon cost and has a much lower online
computational burden [5] than MPC. In RL, two different methodologies can be
distinguished: policy search and dynamic programming [6]. Policy search
directly searches for an optimal policy, while dynamic programming aims to
solve the Bellman equation first and then determines the policy. Dynamic
programming has the advantage over policy search is that it reduces the policy
optimization problem to a one-step look-ahead problem. When applied to systems
with continuous state and input spaces, approximate dynamic programming (ADP),
which uses sample-based methods accompanied by sophisticated function
approximators, has been developed [7, 8]. In this paper, we consider
approximate value iteration (VI), the most basic and direct way to solve the
Bellman equation. Moreover, we provide comprehensive performance guarantees
for stability, safety and sub-optimality of the developed ADP approach.

###  I-B Related work

To reduce the computational cost of MPC, two different types of approaches
have been extensively studied: approximate MPC and reinforcement learning
under constraints.

_Approximate MPC:_ Approximate MPC parameterizes a policy and then uses
supervised learning or gradient-based methods to mimic a predictive control
policy. The online computational cost is thus significantly reduced because
approximate MPC directly computes control actions based on the learned
parameters, rather than online solving an optimization problem. Some
approximate MPC work considers linear systems, with different focuses on,
e.g., stability [9] and constraint satisfaction [10]. Some approximate MPC
approaches can handle nonlinear control problems with constraints, e.g., by
using constraint tightening [11]. However, as these approaches are developed
for general nonlinear systems, they lack deterministic performance
verification tools.

_RL for constrained control:_ RL for constrained control can be roughly
categorized into two groups: policy-projection RL methods, which impose
constraints _after_ learning, and policy-optimization RL methods, which impose
constraints _during_ learning. Using RL for controller synthesis allows
designers to stop the learning process and get a satisfactory policy whenever
the expected performance is met.

For the first group, a predictive safety filter [12, 13] can be adopted to
modify the learned policy, which can be derived from any RL algorithm. This
group of approaches can in general be applied to constrained nonlinear control
problems. For PWA systems with linear constraints, this group needs to solve
mixed-integer convex problems online [13], and as a result it is not suitable
for large-scale systems or for systems requiring fast computation.
Additionally, this group of methods in general jeopardizes optimality and
stability [14].

For the second group, constrained policy optimization [15, 16, 17] is often
used. Most methods [18] consider constrained Markov Decision Process (MDP)
problems, in which constraints are on expected cumulative costs. Recent
developments have been made to transform instantaneous constraints into
constraints on expected cumulative costs. Among these developments, [19]
considers instantaneous expectation constraints, while [20] contemplates
chance constraints. In [21], hard instantaneous constraints are considered,
but the convergence of the RL algorithm relies on the assumptions of the
finiteness of MDP and the continuous differentiability of the value function,
which is not the case for optimal control of PWA systems. Besides, [21, 19] do
not provide stability or safety analysis of the RL policy.

_Lyapunov stable RL:_ In addition to the online computational advantage, the
stability property of RL controllers has been investigated recently. In [22],
stability analysis of approximate VI for unconstrained nonlinear systems is
conducted, but it results in very restrictive conditions, which require the
approximation errors per iteration to be upper bounded. [23, 24] further relax
these conditions by only focusing on the stopping criterion and the
approximation error in the last iteration. Nevertheless, these references
consider unconstrained problems and do not address the sub-optimality of the
RL policy. Besides, the combination of RL and MPC is proposed in [25, 24, 26,
27] as a theoretically justified methodology to produce safe and stabilizing
policies. However, in [25, 24, 26, 27], the online computational issue of MPC
is not solved by introducing RL.

Based on these observations, for the policy optimization methods, one could
conclude that (i) no work has been done yet for PWA systems, and (ii) good and
comprehensive performance of RL-based controllers regarding online computing
convenience, stability, and constraint fulfillment cannot be achieved
simultaneously. In this paper, we will consider policy optimization methods
for PWA systems subject to linear or union-of-polyhedra (UoP) constraints, and
solve the above two issues.

Meanwhile, most RL/ADP methods ignore the properties of the value function and
the policy, e.g., boundedness, continuity, nonlinearity, so that the design of
function approximators is usually done in an empirical way [28]. For PWA
systems with a linear performance index, fortunately, extant research on
explicit hybrid MPC has shown that the optimal value function and the optimal
policy are both PWA [29, 30]. This provides an opportunity to use PWA function
approximation, as will be investigated in this paper.

_Performance verification of learning-based controllers:_ In addition to
controller design, there is some related work on performance analysis and
verification of RL or any learning-based controllers. In [31, 32], a
learner/verifier framework is proposed to verify offline the asymptotic
stability for an input-constrained PWA system controlled by a PWA neural
network (NN). Using MILP to explicitly compute the set of trajectories of
linear systems controlled by an NN can also be applied to evaluate the safety
and stability [33]. However, for nonlinear systems with NN controllers,
achieving practical stability is much easier than achieving asymptotic
stability [23]. For PWA systems with NN controllers, there is currently no
systematic way to verify different properties including practical and
asymptotic stability as well as state constraint satisfaction.

###  I-C Methods and contributions of this paper

In conclusion, using RL to produce a reliable learning-based controller for
constrained PWA systems with performance guarantees and low online
computational requirements is still an open problem. In this paper, we develop
ADP algorithms under linear and UoP constraints. We propose two formulations
for the inclusion of PWA penalties in dynamic programming, i.e., adding
penalties to the stage cost and integrating penalties into the cost-to-go. For
exact dynamic programming algorithms, we prove the continuity and the PWA
property of the value functions when PWA penalty functions and UoP constraints
are considered. These properties allow us to design PWA function
approximators. We then present two different controllers: an implicit
controller that is obtained online by solving an MILP problem that is much
more simple than the one of implicit hybrid MPC, and an explicit controller
that is learned offline by policy gradient. The explicit policy is then
projected online onto the input constraint set, which corresponds to solving a
single convex quadratic program online.

We provide rigorous analysis on the closed-loop stability, safety, as well as
sub-optimality of the controllers. We establish a systematic, MILP-based
procedure that allows us to certify the reliability, in terms of stability and
safety of the closed-loop system.

The paper contributes the state of the art as follows

  1. 1.

This work is the first research on designing policy optimization RL methods
for constrained PWA systems. Systematic performance analysis on the
feasibility, stability, and sub-optimality of the RL-based controllers is
provided. The analysis suggests several ways to employ the proposed algorithms
in practice.

  2. 2.

Compared to MPC, our method exhibits a notable superiority in terms of online
computational simplicity. The proposed explicit controller is also more
computationally tractable than the policy projection methods of [12, 13].

  3. 3.

We develop a mixed-integer optimization-based framework to exactly verify the
stability and safety of the closed-loop system. The framework extends the
verification techniques of [31, 32, 33, 9]. Specifically, we establish a
comprehensive scheme that addresses both practical and asymptotic stability
properties and the enlargement of stable and safe regions while [31, 32, 33,
9] only investigate asymptotic stability.

  4. 4.

For constraints handling, we explore how the proposed two manners of
incorporating PWA penalties can accelerate the learning process. We discuss
their adaptability in different situations.

##  II Preliminaries

_Notations:_ Let $\mathbb{R}=(-\infty,+\infty)$, $\mathbb{R}_{\geq
0}=[0,+\infty)$, and $\mathbb{R}_{>0}=(0,+\infty)$. The boundary of the set
$\mathcal{S}$ is $\partial\mathcal{S}$, and $\mathrm{int}(\mathcal{S})$ stands
for the interior of $\mathcal{S}$. We utilize $A_{i\cdot}$ to represent the
$i$th row of the matrix $A$. We define the sub-level set $\mathcal{B}(J,S)$
for a continuous function $J:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ and a
compact set $S\subseteq\mathbb{R}^{n}$ as
$\mathcal{B}(J,S)\triangleq\\{x\in\mathbb{R}^{n}|J(x)\leq\rho\\}$, where
$\rho=\min_{x\in\partial S}J(x)$. For $a\in\mathbb{R}$, denote by $\lceil
a\rceil$ the smallest integer larger than or equal to $a$.

###  II-A Optimal control of PWA systems

We consider discrete-time PWA systems of the form

$x_{t+1}=f_{\mathrm{PWA}}(x_{t},u_{t})=A_{i}x_{t}+B_{i}u_{t}+f_{i}\quad\text{
if }\left[\begin{array}[]{l}x_{t}\\\
u_{t}\end{array}\right]\in\mathcal{C}_{i},$ (1)

with input and state constraints $x\in X,\;u\in U$111All the results of this
paper also apply to the case when there is a coupled
constraint:$[x^{T}\;\;u^{T}]^{T}\in D$ with $D$ a polyhedron in
$\mathcal{X}\times\mathcal{U}$, by letting $X$ be the projection of $D$ in
$\mathcal{X}$ and by letting $U$ be a time-varying set depending on $x$.. In
(1), $\\{\mathcal{C}_{i}\\}^{s}_{i=1}$ is a polyhedral partition of the state-
input space $\mathcal{X}\times\mathcal{U}$. We want to design a control policy
$u_{t}=\pi(x_{t}),\;t=0,1,\dots$ with $\pi:\mathcal{X}\rightarrow\mathcal{U}$
to simultaneously satisfy the constraints and minimize the infinite-horizon
cost

$J\left({{x_{0}},{\mathbf{u}}}\right)=J_{\pi}(x_{0})=\sum\limits_{t=0}^{\infty}{l\left({{x_{t}},\pi\left(x_{t}\right)}\right)},$
(2)

where
${\mathbf{u}}=\\{u_{t}\\}_{t=0}^{\infty}=\\{\pi(x_{t})\\}_{t=0}^{\infty}$
stands for the infinite-horizon input sequence and the function
$l:\mathbb{R}^{n_{x}}\times\mathbb{R}^{n_{u}}\rightarrow\mathbb{R}_{\geq 0}$
is the stage cost. Throughout the paper, we assume that the dynamics,
constraints, and stage costs satisfy the following assumption.

###### Assumption 1

  * •

_Dynamics:_ The function
$f_{\mathrm{PWA}}\left(\cdot,\cdot\right):\mathcal{X}\times\mathcal{U}\rightarrow\mathcal{X}$
is a continuous PWA function, with $\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$
and $\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$. Besides,
$f_{\mathrm{PWA}}\left(0,0\right)=0$.

  * •

_Constraints:_ The state constraint set $X=\bigcup_{i=1}^{r_{0}}X^{(i)}$ is a
UoP, where $X^{(i)}$ is a polyhedron for each $i=1,\dots,r_{0}$ and $r_{0}$ is
the number of polyhedra. The input constraint set $U$ is a polyhedron.
Moreover, $X\times
U\subset\bigcup_{i=1}^{s}{\mathcal{C}}_{i}=\mathcal{X}\times\mathcal{U}$.

  * •

_Stage cost:_ The stage cost is based on the 1-/$\infty$-norm:
$l\left({{x_{t}},{u_{t}}}\right)=\|Qx_{t}\|_{q_{1}}+\|Ru_{t}\|_{q_{2}}$, where
$q_{1},q_{2}\in\\{1,\infty\\}$ and $Q,R$ have full column rank.

Similar assumptions can be found in other papers on control of PWA systems,
e.g., [29, 34]. In most literature [29, 34, 35], it is usually assumed that
$X$ is a polyhedron, while we generalize it to a UoP. It should be mentioned
that many practical nonlinear and non-convex state constraints, such as
collision avoidance encountered in robotics, can be modeled or outer-
approximated by constraint sets that are UoPs [36].

By combining (1) and (2), our control objective is to solve the constrained
infinite-horizon optimal control problem

$\displaystyle{J^{*}}\left({{x_{0}}}\right)=\mathop{\min}\limits_{\pi,\mathbf{u},\mathbf{x}}$
$\displaystyle J_{\pi}(x_{0})$ s.t.
$\displaystyle{x_{t+1}}={f_{{\text{PWA}}}}({x_{t}},{u_{t}}),\;u_{t}=\pi\left(x_{t}\right),$
$\displaystyle x_{t}\in X,\;u_{t}\in U,\;t=0,1,...,$ (3)

where ${\mathbf{x}}=\\{x_{t}\\}_{t=0}^{\infty}$, and
$J^{*}(\cdot):\mathcal{X}\to[0,\infty]$ is the optimal value function. An
optimal policy, denoted by $\pi^{*}(\cdot)$, minimizes
$J_{\pi}\left(x_{0}\right)$ subject to the constraints in (II-A) for any
initial state $x_{0}$ that makes (II-A) feasible. In this paper, we denote by
$\bar{X}$ the set of feasible initial states $x_{0}$ that make $J^{*}(x_{0})$
finite. In [30], $\bar{X}$ is called the maximal stabilizable set for (1).

###### Assumption 2

The set $\bar{X}$ is non-empty. Furthermore, for any $x_{0}\in\bar{X}$, there
exists a policy $\pi(\cdot)$ such that the system (1) with $u_{t}=\pi(x_{t})$,
starting from $x_{0}$, will reach the origin in a finite number of time steps.

Assumption 2 is a standard stabilizability assumption for discrete-time
systems. Similar assumptions for PWA systems can be found, e.g., in [35, 29].

For any $x\in\mathcal{X}$, according to Bellman’s Principle of Optimality
[37], the value function $J^{*}(\cdot)$ and the optimal policy
$\pi^{*}(\cdot)$ satisfy the following equations:

$\displaystyle J^{*}(x)=\Gamma J^{*}(x)\triangleq\min_{u\in
U}l(x,u)+J^{*}(f_{\mathrm{PWA}}(x,u)),$
$\displaystyle\pi^{*}(x)\in\underset{u\in
U}{\arg\min}\;l(x,u)+J^{*}(f_{\mathrm{PWA}}(x,u)),$ (4)

where $\Gamma$ is called the Bellman operator [37]. In (II-A), the domain of
$J^{*}(\cdot)$ is the whole state space $\mathcal{X}$, which means that the
value of $J^{*}(\cdot)$ goes to infinity outside $\bar{X}$. In general, the
equation for $J^{*}(\cdot)$ in (II-A) may have multiple solutions.
Nevertheless, it follows from [38, Proposition 1] that $J^{*}(\cdot)$ can be
the unique solution that satisfies $J^{*}(0)=0$ under Assumption 2.

###  II-B Exact value iteration

Solving the Bellman equations (II-A) is in general computationally prohibitive
for nonlinear systems in the setting of continuous input and state spaces. For
the PWA system, the equations in (II-A), on the other hand, can be solved by
using an exact value iteration (VI) method [29], which solves multiple multi-
parametric linear programs (mp-LPs). To motivate our ADP methods, we summarize
the exact VI method and discuss its limitations in this section. The exact VI
algorithm starts from an initial value function $J_{0}(\cdot)$ that is either
zero in $\mathcal{X}$ (case 1) or a control Lyapunov function defined on a
subset of $\mathcal{X}$ (case 2). In case 2, the following assumption should
be satisfied.

###### Assumption 3

A continuous and PWA control Lyapunov function
$J_{\mathrm{CL}}(\cdot):X_{\mathrm{CI}}\to\mathbb{R}_{\geq 0}$ in a polyhedral
control-invariant set $X_{\mathrm{CI}}$ is available. In other words,
$\min_{u\in U,f_{\mathrm{PWA}}(x,u)\in
X_{\mathrm{CI}}}\;l(x,u)+J_{\mathrm{CL}}(f_{\mathrm{PWA}}(x,u))-J_{\mathrm{CL}}(x)\leq
0,\;\forall x\in X_{\mathrm{CI}}$.

Assumption 3 frequently appears in the stability analysis of MPC [30], where
$J_{\mathrm{CL}}(\cdot)$ is chosen as the terminal cost and $X_{\mathrm{CI}}$
is specified as the terminal constraint. To satisfy Assumption 3, it is
sufficient to compute a stabilizing piecewise linear feedback law on
$X_{\mathrm{CI}}$, and then $J_{\mathrm{CL}}(\cdot)$ can be computed by
solving some nonlinear inequalities that contain 1-/$\infty$-norm of some
linear functions [35].

With the initialization $X_{0}=X$ and $J_{0}(x)=0,\;\forall x\in X_{0}$ (case
1), or $X_{0}=X_{\mathrm{CI}}$ and $J_{0}(x)=J_{\mathrm{CL}}(x),\;\forall x\in
X_{0}$ (case 2), the exact VI method iterates as follows:

$\displaystyle J_{k}(x)=\min_{u\in U,f_{\mathrm{PWA}}(x,u)\in
X_{k-1}}l(x,u)+J_{k-1}(f_{\mathrm{PWA}}(x,u))$ (5)

for $k=1,2,..$. Here, $X_{k}=\mathrm{Pre}(X_{k-1})\cap X_{0}$, where
$\operatorname{Pre}(\mathcal{S})=\left\\{x\in\mathcal{X}|\exists u\in U\text{
s.t. }f_{\mathrm{PWA}}(x,u)\in\mathcal{S}\right\\}$ is the backward-reachable
set to a set $S$. Even for PWA systems with polytopic state constraints, the
backward-reachable set to a polyhedral set can be a non-convex UoP because of
the nonlinear dynamics (1), which means that $X_{k},k=1,2,\dots$ can be non-
convex UoPs [4].

The resulting $J^{*}(\cdot)$ and $\pi^{*}(\cdot)$ are both PWA functions
sharing the same polyhedral partition of the feasible region $\bar{X}$ [29].
However, the complexity (i.e., the number of polyhedral regions or affine
functions) of $J^{*}(\cdot)$ and $\pi^{*}(\cdot)$ is exponential in both the
dimension of the system and the number of constraints in (II-A) [39], so that
storing the affine functions and regions of $J^{*}(\cdot)$ and
$\pi^{*}(\cdot)$ needs a huge amount of memory. Secondly, the number of mp-LPs
that need to be solved per iteration also grows exponentially with the problem
dimension. Moreover, the online implementation of [29] needs to search which
polyhedron the measured state belongs to. For high-dimensional systems, these
regions may have complex representations. The growing complexities of the
explicit controller’s structure, and the offline and online computations limit
the applicability of exact VI to small-scale systems [2].

Based on these observations, it is thus necessary to simplify both the
procedure of solving the Bellman equation and the control policy, by using
some approximation methods.

However, the VI formulation in (5) is not suitable for approximation because
the probably non-convex constraint $f_{\mathrm{PWA}}(x,u)\in X_{k-1}$ leads to
too complex optimization problems when using sample-based approaches.

##  III Value Iteration with Penalty Functions

To deal with this issue, we consider soft state constraints by defining a
penalty function $P(\cdot,\cdot)$. Suppose that each $X^{(i)}$ of the UoP
constraint set $X$ is characterized by the linear inequality $E^{(i)}_{X}x\leq
g^{(i)}_{X}$, where $E^{(i)}_{X}\in\mathbb{R}^{m^{(i)}_{x}\times n_{x}}$ and
$g^{(i)}_{X}\in\mathbb{R}^{m^{(i)}_{x}}$, we design the penalty function
$P(\cdot,\cdot)$ as the following min-max forms:

$\displaystyle P(x,X)$
$\displaystyle=p\min_{i}\max\left\\{0,(E^{(i)}_{X})_{1\cdot}x-(g^{(i)}_{X})_{1},\dots,\right.$
$\displaystyle\quad\left.(E^{(i)}_{X})_{m^{(i)}_{x}\cdot}x-(g^{(i)}_{X})_{m^{(i)}_{x}}\right\\}\;\text{or}$
$\displaystyle P(x,X)$
$\displaystyle=p\min_{i}\sum_{j=1}^{m^{(i)}_{x}}\max\left\\{0,(E^{(i)}_{X})_{j\cdot}x-(g^{(i)}_{X})_{j}\right\\},$
(6)

where $(g^{(i)}_{X})_{j}$ is the $j$th element of the vector $g^{(i)}_{X}$,
$m^{(i)}_{x}$ is the number of rows of $E^{(i)}_{X}$, and $p>0$ is the
constraint violation penalty weight. When $X$ reduces to a polyhedron, i.e.,
$r_{0}=1$, the minimum operator in (III) will be removed.

An important property of $P(\cdot,\cdot)$ is that $P(\cdot,\cdot)$ is a PWA
function w.r.t. its first argument. This means that adding such a penalty
function into the cost function in (II-A) will not change the PWA properties
of the optimal value function and the optimal policy. Note that we avoid
barrier functions such as the logarithmic barrier function because they can go
to infinity in any compact constraint sets and will deprive the value function
of the PWA property. Besides, we do not penalize the input constraint
violation, because the input constraints are single polyhedral constraints
that can be readily handled.

In case 2 of (5), in addition to enforcing a soft penalty for $X$, we need to
reconstruct the initial value function since $J_{\mathrm{CL}}(\cdot)$ is
undefined outside $X_{\mathrm{CI}}$. To achieve this, we need to penalize
$J_{{\text{CL}}}\left(x\right)$ for $x$ outside $X_{\mathrm{CI}}$ by finite
values:

$J_{0}^{{\text{soft}}}(x)=\left\\{\begin{gathered}{J_{{\text{CL}}}}\left(x\right),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;x\in
X_{\mathrm{CI}},\hfill\\\
{J_{{\text{CL}}}}\left({\bar{z}}\right)+P(x,X_{\mathrm{CI}}),\;x\notin
X_{\mathrm{CI}},\hfill\\\ \end{gathered}\right.$ (7)

where $\bar{z}(\cdot)$ is one of the optimizers of the following mp-LP:

$\bar{z}(x)\in\arg\min_{z\in X_{\mathrm{CI}}}||z-x||_{\infty}.$ (8)

In (7), $J_{0}^{{\text{soft}}}(\cdot)$ is continuous on $\mathcal{X}$, which
will be proven in Theorem 1. Based on the defined penalty function, a VI
algorithm with penalty is developed as follows

Algorithm 1 Value iteration with penalty

1:Output: A value function
$J^{\mathrm{soft}}_{k-1}(\cdot):\mathcal{X}\to\mathbb{R}_{\geq 0}$.

2:Initialize the value function (option (a))
$J^{\mathrm{soft}}_{0}(x)\leftarrow 0,\;\forall x\in\mathcal{X}$, or (option
(b)) by (7).

3:for $k=1,2,\dots$ do

4: the value iteration
$J_{k}^{\mathrm{soft}}(x)\leftarrow\Gamma_{\mathrm{p,\alpha}}J_{k-1}^{\mathrm{soft}}(x),\;\forall
x\in\mathcal{X}$, where $\alpha=1$ if option 1 is chosen, or $\alpha=2$ if
option 2 is chosen, and $\Gamma_{\mathrm{p,\alpha}}$ is defined in (9) and
(III).

5: If $J_{k}^{\mathrm{soft}}(x)=J_{k-1}^{\mathrm{soft}}(x),\;\forall
x\in\mathcal{X}$, break.

6:end for

In Algorithm 1, we consider two options for the VI, in which we define two
Bellman operators for $J:\mathcal{X}\to\mathbb{R}$. The first one used in
option 1 is

$\Gamma_{\mathrm{p,1}}J(x)\triangleq\min_{u\in
U}\;l_{\mathrm{p}}(x,u)+J\left(f_{\mathrm{PWA}}(x,u)\right),x\in\mathcal{X},$
(9)

where $l_{\mathrm{p}}(x,u)=l(x,u)+P(x,X)$. The second one used in option 2 is

$\displaystyle\Gamma_{\mathrm{p,2}}J(x)\\!\triangleq\\!\min_{u\in
U}\;l(x,u)\\!+\\!J\left(f_{\mathrm{PWA}}(x,u)\right)$
$\displaystyle\\!+\\!P_{k-1}\left(f_{\mathrm{PWA}}(x,u)\right),$
$\displaystyle x\in\mathcal{X},$ (10)

where $P_{0}(x)=0,\forall x\in\mathcal{X}$ and $P_{k}(x)=P(x,X),\forall
x\in\mathcal{X}$ and $\forall k>0$. Since the state constraints are removed,
the working region of VI is the whole state space $\mathcal{X}$. Besides, we
also consider two options (options (a) and (b)) for the initialization of the
value function. The combinations of above options result in four different
options: options 1(a), 1(b), 2(a), and 2(b). Note that adding a penalty into
the stage cost, which is applied in option 1 of Algorithm 1, is common in
existing constrained ADP methods [40]. In comparison, in option 2 of the
algorithm, we propose a novel scheme in which we add a penalty into the cost-
to-go $J_{k-1}^{\mathrm{soft}}\left(f_{\mathrm{PWA}}(x,u)\right)$. In the
following theorem, we will analyze the PWA property and continuity of each
$J_{k}^{\mathrm{soft}}(\cdot)$ as well as the convergence of the sequence
$\\{J_{k}^{\mathrm{soft}}(\cdot)\\}^{\infty}_{k=0}$ to a fixed optimal value
function in all options.

###### Theorem 1

Considering Algorithm 1, if Assumptions 1-2 hold in option (a) and Assumptions
1-3 hold in option (b), each $J_{k}^{\mathrm{soft}}(\cdot),\;k<\infty$ is a
continuous PWA function on $\mathcal{X}$ and the value function sequence
$\\{J_{k}^{\mathrm{soft}}(x)\\}^{\infty}_{k=0}$ converges point-wise to

$\displaystyle{J^{\mathrm{soft}*}}\left({{x}}\right)=\mathop{\min}\limits_{\\{u_{i},x_{i}\\}^{\infty}_{i=0}}$
$\displaystyle\sum_{i=0}^{\infty}l_{p}\left(x_{i},u_{i}\right)+P(x_{i},X)$
$\displaystyle\mathrm{s.t.}\;\;\;\;$
$\displaystyle{x_{i+1}}={f_{{\text{PWA}}}}({x_{i}},{u_{i}}),\;u_{i}\in
U,\;i=0,1,...$ $\displaystyle x_{0}=x$

for all $x\in X$.

In the proof of Theorem 1 (see Appendix A-A), it is shown that options 1 and 2
lead to the same $J_{k}^{\mathrm{soft}}(x)$ for $x\in X$, even if their ways
of iteration are different. Detailed comparison between options 1 and 2 is
given in Section IV. C.

If the PWA system reduces to a linear time-invariant system and $X$ is a
polyhedron, it is possible to make $J_{k}^{\mathrm{soft}}(x)=J_{k}(x),\forall
x\in X_{k}$, by choosing a large enough $p$, according to the exact penalty
theorem [41, Proposition 5.4.5].

##  IV Constrained ADP algorithm

###  IV-A Algorithm design

Continuity of the value functions, established in Theorem 1, is desired since
it enables a universal approximation capability [42]. With Algorithm 1, a
tractable ADP approach can be developed to approximate each
$J_{k}^{\mathrm{soft}}(\cdot)$. In particular, a function approximator
(critic) $\hat{J}_{k}\left(\cdot,\theta_{k}\right)$, which is parameterized by
$\theta_{k}$, is constructed to replace $J_{k}^{\mathrm{soft}}(\cdot)$. In
each iteration $k$, a set $X_{\mathrm{s}}=\\{x^{(i)}\\}^{N_{x}}_{i=1}$ of
state samples is collected from a compact region of interest
$\Omega_{k}\subseteq\mathcal{X}$, according to some strategies such as
sampling from a uniform grid and random sampling [6]. Here, $N_{x}$ is the
number of samples. The update of the parameter $\theta_{k}$ minimizes
$\sum_{i=1}^{N_{x}}[\Gamma_{\mathrm{p},\alpha}\hat{J}_{k-1}(x,{\theta_{k-1}})|_{x=x^{(i)}}-{\hat{J}_{k}}({x^{(i)},{\theta_{k}}})]^{2}$.
The iterative procedure stops when the difference between $\theta_{k}$ and
$\theta_{k-1}$ is small enough. The detailed procedure is given in Algorithm
2.

Algorithm 2 Constrained approximate value iteration

1:Output: A value function $\hat{J}_{k-1}\left(\cdot,\theta_{k}\right)$

2:Option (a): Initialize the value function
$\hat{J}_{0}(\cdot,\theta_{0})\leftarrow 0,\forall x\in\Omega_{0}$, where
$\Omega_{0}=\mathcal{X}$ in option 1(a) or $\Omega_{0}=X$ in option 2(a).

3:Option (b): Initialize the value function $\hat{J}_{0}(\cdot,\theta_{0})$ by
${\theta_{0}}\leftarrow\mathop{\arg\min}\limits_{\theta}\sum\limits_{x^{(i)}\in
X_{\mathrm{s}}\cap\Omega_{0}}{{{\rho_{v}(x^{(i)})({J_{0}^{\mathrm{soft}}(x^{(i)})-{{\hat{J}}_{0}}({x^{(i)}},{\theta})})}^{2}}},$
(11) where $J_{0}^{\mathrm{soft}}(\cdot)$ is from (7), and
$\Omega_{0}=\mathcal{X}$ in option 1(b) or $\Omega_{0}=X_{\mathrm{CI}}$ in
option 2(b).

4:for $k=1,2,\dots$ do

5: if option 1 is chosen, then let $\Omega_{k}\leftarrow\mathcal{X}$ and
$\alpha\leftarrow 1$.

6: end if

7: if option 2 is chosen, then let
$X_{k}\leftarrow\operatorname{Pre}\left(X_{k-1}\right)\cap X$,
$\Omega_{k}\leftarrow X_{k}$, and $\alpha\leftarrow 2$.

8: end if

9: Obtain the target value $v^{(i)}_{k}$ by
$v^{(i)}_{k}\leftarrow\Gamma_{\mathrm{p,\alpha}}\hat{J}_{k-1}(x,\theta_{k-1})|_{x=x^{(i)}},\forall
x^{(i)}\in X_{\mathrm{s}}\cap\Omega_{k}$ (12)

10: Find $\theta_{k}$ such that
${\theta_{k}}\\!\leftarrow\\!\mathop{\arg\min}\limits_{\theta}\\!\sum\limits_{x^{(i)}\\!\in\\!X_{\mathrm{s}}\\!\cap\Omega_{k}}\\!{{{\rho_{v}(x^{(i)})({v_{k}^{(i)}\\!-\\!{{\hat{J}}_{k}}({x^{(i)}},{\theta})})}\\!^{2}}}.$
(13)

11: if
$|\hat{J}_{k}(x,\theta_{k})-\hat{J}_{k-1}(x,\theta_{k-1})|\leq\epsilon(x),\;\forall
x\in\Omega_{k}\cap\Omega_{k-1}$, then break.

12: end if

13:end for

In (13) of step 8 of Algorithm 2,
$\rho_{v}(\cdot):\mathcal{X}\to\mathbb{R}_{>0}$ is the state relevance
weighting function. In step 9, $\epsilon(\cdot):\mathcal{X}\to\mathbb{R}_{\geq
0}$ is a tolerance function, determining whether
$\hat{J}_{k-1}\left(\cdot,\theta_{k-1}\right)$ is satisfactory. Both of
$\rho_{v}(\cdot)$ and $\epsilon(\cdot)$ will be designed later in Section V.A.
Besides, in practice one would also need a limit on the maximum number of
iterations.

With $\hat{J}_{k-1}\left(\cdot,\theta_{k-1}\right)$ available, a sub-optimal
control policy $\hat{\pi}^{\mathrm{im}}(x)$ can be implicitly determined by

$\hat{\pi}^{\mathrm{im}}(x)\in\arg\min_{u\in
U}\;l(x,u)+\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,u),\theta_{k-1}),\;\forall
x\in\Omega_{k}$ (14)

In step 1 of Algorithm 2, the function approximator is initialized by
regressing $J_{0}^{\mathrm{soft}}(\cdot)$, provided that the explicit form of
$J_{\mathrm{CL}}(\cdot)$ is known. If the explicit form of
$J_{\mathrm{CL}}(\cdot)$ is not available but a stabilizing and safe piecewise
linear feedback law $\pi_{\mathrm{PWL}}(\cdot):\mathcal{X}\to\mathcal{U}$ on
$X_{\mathrm{CI}}$ is known, we can initialize
$\hat{J}_{0}\left(\cdot,\theta_{0}\right)$ as an approximation of
$J_{\pi_{\mathrm{PWL}}}(\cdot)$, which is also a control Lyapunov function, by
doing a policy evaluation [6].

To carry out the iterative procedure in Algorithm 2 more efficiently, we need
to use a proper function approximator at each iteration. Since each
$J^{\mathrm{soft}}_{k}(\cdot)$ is a PWA function, it is preferable that the
candidate approximator can also output a PWA function. Suitable choices are
thereby NNs with (leaky) rectifier linear units (ReLUs) as activation
functions, difference of two max-affine functions, Min-Max NNs [43], and so
on. Detailed descriptions of these function approximators are given in the
Appendix A-C.

_Computation of $X_{k}$ under the UoP state constraint:_ In option 2, one is
required to compute the reachable sets $X_{k}$. The following lemma states
that we are able to obtain $X_{k}$ by performing polyhedral operations even if
the state constraint is UoP. The proof is given in Appendix A-B.

###### Lemma 1

For the PWA system (1) with a convex polyhedral input constraint $u\in U$ and
a UoP state constraint $x\in X=\bigcup_{i=1}^{r_{0}}X^{(i)}$, the set iterates
$X_{k}=\mathrm{Pre}(X_{k-1})\cap X,\;k=1,2,\dots$ make each $X_{k}$ be a UoP.
Then, if $X_{k-1}$ is in the form of
$X_{k-1}=\bigcup_{i=1}^{r_{k-1}}X^{(i)}_{k-1}$, where each $X^{(i)}_{k-1}$ is
a polyhedron, then $X_{k}$ can be computed by

$X_{k}=\bigcup_{j=1}^{r_{0}}\bigcup_{i=1}^{r_{k-1}}\mathrm{Pre}(X^{(i)}_{k-1})\cap
X^{(j)}$ (15)

_Mixed-integer formulations of problems ( 12) and (14):_ Problems (12) and
(14) have similar forms. They can be transformed into MILP problems since both
the PWA system and PWA function approximators are MILP representable, which is
shown in Appendix A-C. Here, we say a function $J$ is MILP representable if
$J$ can be represented by a set of mixed-integer linear equations and
inequalities containing additional variables. We note this set by
$\mathrm{gr}_{\mathrm{MILP}}(J)$. As a result, we obtain mixed-integer
formulations of problems (12) and (14). Precisely, if, e.g., the first type of
the penalty function in (III) and the infinity norm of the stage cost are
chosen, (12) in option 1 can be equivalently written as the following MILP:

$\displaystyle\min_{\begin{subarray}{c}u,\varepsilon^{(x)},\varepsilon^{(u)},\\\
\varepsilon^{(J)},\varepsilon^{(P)},x^{+}\end{subarray}}\;\varepsilon^{(x)}+\varepsilon^{(u)}+\varepsilon^{(J)}+\varepsilon^{(P)}$
$\displaystyle\mathrm{s.t.}-1_{n_{x}}\varepsilon^{(x)}\leq\pm
Qx,\;-1_{n_{u}}\varepsilon^{(u)}\leq\pm Ru,\;u\in U,$
$\displaystyle\quad(x^{+},\varepsilon^{(J)})\in\mathrm{gr}_{\mathrm{MILP}}(\hat{J}_{k-1}),\;(x,\varepsilon^{(P)})\in\mathrm{gr}_{\mathrm{MILP}}(P(\cdot,X)),$
$\displaystyle\quad([x^{T}\;u^{T}]^{T},x^{+})\in\mathrm{gr}_{\mathrm{MILP}}(f_{\mathrm{PWA}}).$

Similar transformations can be achieved for problem (14), for option 2 of
Algorithm 2, and for the other choice of $P(\cdot,\cdot)$ and
$l(\cdot,\cdot)$.

MILP problems can be effectively solved by using the branch-and-bound approach
[44], which is a global optimization algorithm.

Different from (12), Problem (14) is solved online. Even through it still
belongs to an MILP problem, (14) can be solved more rapidly than a general
hybrid MPC problem with a long horizon, which results in plenty of auxiliary
and binary variables with mixed-integer formulation. That is one of main
benefits of using ADP or RL.

###  IV-B Approximating explicit policies

Since two PWA functions $f_{\mathrm{PWA}}(\cdot,\cdot)$ and
$\hat{J}_{k-1}(\cdot,\theta_{k-1})$ are coupled in (14), (14) may still have
many auxiliary and binary variables, if, e.g., a multiple-layer (deep) NN is
used. As (14) needs to be solved online, the advantage of low computational
complexity brought by ADP is not obvious. To avoid solving complex MILP
problems online, the policy $\hat{\pi}^{\mathrm{im}}(\cdot)$ can also be
represented explicitly, in which case it usually needs to be approximated by a
second function approximator (actor). The actor is also recommended to having
a PWA form since the optimal control policy $\pi^{*}(\cdot)$ is PWA.

As the optimizer $\hat{\pi}^{\mathrm{im}}(\cdot)$ can be discontinuous and not
unique, instead of using supervised learning methods to train the actor, we
can directly construct a parameterized policy
$\hat{\pi}^{\mathrm{ex}}(\cdot,\omega)$ with parameter $\omega$ and update
$\omega$ to minimize the expectation of the objective function in (14) w.r.t.
the sample distribution $d_{s}$ used in Algorithm 2. This results in the
following policy optimization problem

$\displaystyle\begin{aligned}
\omega^{*}\in\arg\min_{\omega}\;\mathbb{E}_{x\sim
d_{s}}&[\rho_{\pi}(x)(l(x,\hat{\pi}^{\mathrm{ex}}(x,\omega))\\\
&+\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{ex}}(x,\omega)),\theta_{k-1}))]\end{aligned}$
$\displaystyle\quad\quad\quad\quad\mathrm{s.t.}\;\mathbb{E}_{x\sim
d_{s}}[\hat{\pi}^{\mathrm{ex}}(x,\omega)]\in U,$ (16)

where $\rho_{\pi}(\cdot):\mathcal{X}\to\mathbb{R}_{>0}$ is another state
relevance weighting function to be specified later in Section V.A. Similar to
the critic, we specify $\hat{\pi}^{\mathrm{ex}}(\cdot,\omega)$ as a PWA
approximator.

To solve (IV-B), the policy gradient method, combined with the Lagrangian
relaxation approach [18] or the augmented Lagrangian approach [19] for
constraints handling, can be employed.

The above procedures are conducted offline. Ideally, if there are no
approximation errors on both the critics and the actor, and the penalty weight
$p$ and the number of iterations are infinite, we have $X_{\infty}=\bar{X}$
and
$\hat{\pi}^{\mathrm{ex}}(\cdot,\omega^{*})=\hat{\pi}^{\mathrm{im}}(\cdot)$.
Consequently, $\hat{\pi}^{\mathrm{ex}}(\cdot,\omega^{*})$ will always make the
system satisfy all the constraints for the initial condition
$x_{0}\in\bar{X}$. However, due to approximation errors and the finite penalty
weight, the policy cannot always satisfy the state and input constraints. As
the input constraints are usually hard constraints, in the online setting we
project $\hat{\pi}^{\mathrm{ex}}(\cdot,\omega^{*})$ onto $U$ when the current
state $x_{t}$ is received. This results in a convex quadratic program:

$\phi(u^{\mathrm{ex}}_{t})=\arg\min_{u\in U}||u-u^{\mathrm{ex}}_{t}||_{2},$
(17)

where
$u^{\mathrm{ex}}_{t}=\hat{\pi}^{\mathrm{ex}}\left(x_{t},\omega^{*}\right)$,
and the function $\phi(\cdot):\mathbb{R}^{n_{u}}\to\mathbb{R}^{n_{u}}$ maps
the output of the actor to its projected value. Problem (17) can be treated as
a parametric quadratic program with the parameter $u_{t}^{\mathrm{ex}}$.
Therefore, the optimizer $\phi(\cdot)$ of (17) is unique, PWA [30, Theorem
6.7], and also MILP representable [9, Lemma 4]. Meanwhile, (17) defines a
projected policy
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)=\phi(\hat{\pi}^{\mathrm{ex}}\left(\cdot,\omega^{*}\right))$
of $\hat{\pi}^{\mathrm{ex}}\left(\cdot,\omega^{*}\right)$.

###  IV-C Discussions on the ADP method

_Comparison between options 1 and 2 of Algorithm 2:_ The main differences
between options 1 and 2 are the training region $\Omega_{k}$ and the way each
$\hat{J}_{k-1}(\cdot,\theta_{k-1})$ iterates. These differences lead to
differences in the adaptability and efficiency of options 1 and 2:

  * •

Compared to option 1, option 2 is more efficient in sampling and can result in
a better approximation accuracy, because the working region $\Omega_{k}$ in
option 2 is in general much smaller than $\mathcal{X}$. This is also the main
advantage of adding penalties into the cost-to-go over adding penalties into
the stage cost. We note that to implement option 1, one should choose a region
of interest for sampling, and the region must be larger than $X$, so that the
constraint violation can be penalized in the critic. However, the states that
can be steered to the origin and have zero constraint violation are all
contained in $\bar{X}$, which is much smaller than $\mathcal{X}$. Accordingly,
in option 2, we concentrate on $X_{k}$, which converges to $\bar{X}$ as $k$
goes to infinity.

  * •

On the other hand, option 2 needs to compute the $k$-step controllable set
$\Omega_{k}$ ($X_{k}$) while option 1 does not. Therefore, for large-scale PWA
systems, for which the exact computation of each $X_{k}$ is computationally
very demanding, option 1 is preferable.

_Comparison to RL with a safety filter:_ In the schemes of [13, 12], RL
policies are projected onto a safe set where both state and input constraints
are considered. That design is motivated from the fact that the policies in
[13, 12] are derived from standard RL algorithms that do not account for
constraints. In comparison, our proposed ADP algorithms already consider the
state constraints during the offline learning process, so we only need to
modify the policy to satisfy the input constraints. Nevertheless, as the
penalty weight is finite, the state constraint satisfaction cannot be strictly
ensured. Therefore, in Section V we provide a tool to offline verify the state
constraint satisfaction.

_Remark:_ We can strictly guarantee the state constraint satisfaction a priori
by projecting the successor state onto a controlled-invariant set [13]. This
approach, however, will make the online computational complexity very large.

##  V Performance analysis and verification

In this section, we will characterize the stability and safety of the closed-
loop system with the policies $\hat{\pi}^{\mathrm{im}}(\cdot)$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$, and also the sub-optimality
properties of these policies. Firstly, we provide general conditions under
which stability and safety hold. These conditions can guide the parameter
tuning of Algorithm 2. Then, we give sub-optimality guarantees, i.e., a bound
on the mismatch between the infinite cost of the policies and real value
functions. Finally, we develop verifiable stability and safety conditions. We
say that a closed-loop system is safe in a set if its states and inputs
satisfy the constraints for all trajectories starting from the set.

###  V-A Stability and safety analysis

First, we state a useful lemma that gives some properties of the value
function $J_{k}^{\mathrm{soft}}(\cdot)$.

###### Lemma 2

Consider Algorithm 1. Suppose that Assumptions 1-2 hold in option (a) and
Assumptions 1-3 hold in option (b).

(i) Then, there exists a positive constant $\gamma<\infty$ such that for all
$k\geq 0$, $J_{k}^{\mathrm{soft}}(x)\leq\gamma l(x,0)$, $\forall x\in\bar{X}$;

(ii) there exists a finite $\bar{k}>0$ such that for all $k\geq\bar{k}$, we
have

$J_{k}^{\mathrm{soft}}(x)-J_{k-1}^{\mathrm{soft}}(x)\leq\beta l(x,0),\;\forall
x\in\bar{X},\;\mathrm{with}\;\beta\in(0,1).$ (18)

The proof of Lemma 2 is in Appendix A-D. Then, we are ready to state the main
result in this subsection.

###### Theorem 2

Consider Algorithm 2 and the proposed policies
$\hat{\pi}^{\mathrm{im}}(\cdot)$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$. Let $\Omega$ be a compact
subset of $X$. Suppose that Assumptions 1-2 hold in option (a) and Assumptions
1-3 hold in option (b). Consider the follow conditions:

(C1): There exist a constant $\zeta\in(0,1)$ and a positive integer $k$ such
that
$\left|{{{{{\hat{J}}_{k-1}}(x)-J_{k-1}^{\mathrm{soft}}(x)}}}\right|\leq\zeta{{J_{k-1}^{\mathrm{soft}}(x)}},\;\forall
x\in\Omega$.

(C2): There exist a constant $e_{\mathrm{p}}>0$ and a positive integer $k$
such that
$\hat{J}_{\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}}(x)-\hat{J}_{\hat{\pi}^{\mathrm{im}}}(x)\leq
e_{\mathrm{p}}l(x,0),\;\forall x\in\Omega$. Here, $\hat{J}_{\pi}(\cdot)$ is
defined by $\hat{J}_{\pi}(x)\triangleq
l(x,\pi(x))+\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,\pi(x)))$.

As a result, we have:

(i) If C1 holds with $k\geq\bar{k}$ and

$(1+\zeta)(1-\beta)>\max(2\zeta\gamma,\;1),$ (19)

where $\bar{k}$, $\beta$, and $\gamma$ come from Lemma 2, the closed-loop
system
$x_{t+1}=f_{\mathrm{PWA}}(x_{t},\hat{\pi}^{\mathrm{im}}(x_{t})),\;t=0,1,...$
is asymptotically stable and safe in
$\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$.

(ii) If C1 and C2 hold with $k\geq\bar{k}$, and

$(1+\zeta)(1-\beta)>\max(2\zeta\gamma+e_{\mathrm{p}},\;1),$ (20)

the closed-loop system
$x_{t+1}=f_{\mathrm{PWA}}(x_{t},\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(x_{t})),\;t=0,1,...$
is asymptotically stable and safe in
$\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$.

Although Theorem 2 provides sufficient conditions for stability, some of them
(e.g., C1 and C2) are difficult to verify, unless the explicit form of
$J_{k-1}^{\mathrm{soft}}(\cdot)$ is known. A method that can verify C1 and C2
in a probabilistic way is reported in [11]. On the other hand, Theorem 2
suggests several ways to design the parameters and function approximators in
Algorithm 2.

  * •

All results in Theorem 2 require (18) to hold. The left-hand side of (18) is
about the residual error of VI. It indicates that a suitable tolerance
function $\epsilon(\cdot)$, which determines the stopping condition at step 9
of Algorithm 2, could be $\epsilon(x)=e_{\mathrm{tole}}l(x,0)$, for some
$e_{\mathrm{tole}}\in(0,1)$.

  * •

Condition C1 limits the mismatch between ${\hat{J}}_{k-1}(\cdot)$ and
$J_{k-1}^{\mathrm{soft}}(\cdot)$, which further limits the approximation error
$\upsilon_{i}(x)\triangleq\hat{J}_{i}(x)-\Gamma_{\mathrm{p},\alpha}\hat{J}_{i-1}(x),\;i=1,...,k$
of VI. To make $\zeta$ as small as possible, which helps to fulfill (19) and
(20), $\rho_{v}(\cdot)$ in (13) could be specified by
$\rho_{v}(x)=1/l^{2}(x,0)$. To understand this, we consider the state
trajectory $x_{0},\;x_{1},...,x_{k}$ that is generated from the closed-loop
system $x_{t+1}=f_{\mathrm{PWA}}(x_{t},\pi_{k-t}(x_{t})),\;t=0,...,k-1$, where
$\pi_{i}(\cdot)$ denotes the optimizer of
$\Gamma_{\mathrm{p},\alpha}J^{\mathrm{soft}}_{i-1}$. Then, we have

$\displaystyle\hat{J}_{k}(x_{0})-J^{\mathrm{soft}}_{k}(x_{0})$
$\displaystyle\leq\hat{J}_{k-1}(x_{1})-J^{\mathrm{soft}}_{k-1}(x_{1})+\upsilon_{k}(x_{0})$
$\displaystyle\leq\hat{J}_{0}(x_{k})-J^{\mathrm{soft}}_{0}(x_{k})+\sum_{i=1}^{k}\upsilon_{i}(x_{k-i}).$

Suppose that $|\hat{J}_{0}(x_{k})-J^{\mathrm{soft}}_{0}(x_{k})|\leq
e_{\mathrm{v}}l(x_{k},0)$ and $|\upsilon_{i}(x_{k-i})|\leq
e_{\mathrm{v}}l(x_{k-i},0),\;i=1,...,k$, where $e_{\mathrm{v}}>0$, we obtain

$\hat{J}_{k}(x_{0})\\!-\\!J^{\mathrm{soft}}_{k}(x_{0})\leq
e_{\mathrm{v}}\sum_{i=0}^{k}l(x_{i},0)\leq
e_{\mathrm{v}}J_{k}^{\mathrm{soft}}(x_{0}).$ (21)

Similar procedures can be applied to upper bound the value of
$J^{\mathrm{soft}}_{k}(x_{0})-\hat{J}_{k}(x_{0})$. From (11) and (13), we see
that compared to letting $\rho_{v}(x)=1$, our choice of $\rho_{v}(\cdot)$ is
more likely to lead to a smaller $e_{\mathrm{v}}$, which can contribute to the
reduction of $\zeta$, according to (21). Moreover, to circumvent the
singularity of $\rho_{v}(\cdot)$ at the origin, we let
$\rho_{v}(x)=1/(l^{2}(x,0)+\rho)$, with a small positive constant $\rho$.

  * •

Condition C2 is related to the approximation performance of the policy
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$. Similar to the design of
$\rho_{v}(\cdot)$, we would select $\rho_{\pi}(x)=1/(l(x,0)+\rho)$ in (IV-B),
to make $e_{\mathrm{p}}$ small.

  * •

C1 indicates that $\hat{J}_{k-1}(x)\geq 0,\;\forall x\in\Omega$ and
$\hat{J}_{k-1}(0)=0$. To fulfill these, a modification of the PWA function
approximator in [32] can be adopted. In particular, for any PWA function
approximator $\bar{J}(\cdot):\mathbb{R}^{n_{x}}\to\mathbb{R}$, we can
construct $\hat{J}_{i}(\cdot)$ via
$\hat{J}_{i}(x)=\bar{J}(x)-\bar{J}(0)+|R_{v}x|_{q}$. Here, $R_{v}$ is a full
column rank matrix that is included in the parameter $\theta$, and
$q\in\\{1,\infty\\}$. Similarly, C2 implies that
$\hat{\pi}^{\mathrm{ex}}(0)=0$. To realize this, for any PWA function
approximator $\bar{\pi}(\cdot):\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{u}}$, we
can let $\hat{\pi}^{\mathrm{ex}}(x)=\bar{\pi}(x)-\bar{\pi}(0)$.

Different from the existing stability results on ADP [23, 24, 22, 45], Theorem
2 considers the effects of state constraints. It implies that constraint
satisfaction can only be guaranteed in the sub-level set
$\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$,
which can be small, even if we have extremely small approximation errors.
Besides, condition C1 allows us to analyze the sub-optimality of
$\hat{\pi}^{\mathrm{im}}(\cdot)$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$, which has not been addressed
in existing research [23, 24, 22, 45].

###  V-B Sub-optimality analysis

Based on Theorem 2, we can compute upper bounds on
$J_{\hat{\pi}^{\mathrm{im}}}(\cdot)$ and
$J_{\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}}(\cdot)$, which are the infinite-
horizon costs of $\hat{\pi}^{\mathrm{im}}(\cdot)$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$, defined in (2).

###### Corollary 1

Consider Algorithm 2 and the proposed policies
$\hat{\pi}^{\mathrm{im}}(\cdot)$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$. Let $\Omega$ be a compact
subset of $X$.

(i) Let the assumptions in (i) of Theorem 2 hold and $2\zeta\gamma<1$. Then,
for any
$x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$,
we have the inequality

$J^{\mathrm{soft*}}(x)\leq
J_{\hat{\pi}^{\mathrm{im}}}(x)\leq\frac{1-\zeta}{1-2\zeta\gamma}J^{\mathrm{soft}}_{k-1}(x).$
(22)

(ii) Let the assumptions in (ii) of Theorem 2 hold and
$2\zeta\gamma+e_{\mathrm{p}}<1$. Then, for any
$x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$,
we have the inequality

$J^{\mathrm{soft*}}(x)\leq
J_{\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}}(x)\leq\frac{1-\zeta}{1-2\zeta\gamma-
e_{\mathrm{p}}}J^{\mathrm{soft}}_{k-1}(x).$ (23)

(22) and (23) confirm the intuition that a smaller approximation error of the
critic (and the actor) leads to tighter sub-optimality guarantees. Namely, as
$\zeta\to 0,\;e_{\mathrm{p}}\to 0$, and $k\to\infty$, we have
$J_{\hat{\pi}^{\mathrm{im}}}(x)\to J^{\mathrm{soft*}}(x)$ and
$J_{\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}}(x)\to J^{\mathrm{soft*}}(x)$ for
any
$x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$.
Moreover, if option (a) of Algorithm 2 is used, $J^{\mathrm{soft}}_{k-1}(x)$
in (22) and (23) can be replaced by $J^{\mathrm{soft*}}(x)$ because
$J^{\mathrm{soft}}_{k-1}(x)\leq J^{\mathrm{soft*}}(x),\;\forall
x\in\mathcal{X}$.

###  V-C Stability and safety verification

As mentioned in the previous subsection, conditions C1 and C2 in Theorem 2 can
only be evaluated statistically. Moreover, the conditions in Theorem 2 are
sufficient conditions for $\hat{J}_{k-1}(\cdot)$ and
$J^{\mathrm{soft}}_{k-1}(\cdot)$ to be Lyapunov functions, and thus these
conditions can be conservative. Besides, sometimes only practical stability
can be ensured for nonlinear systems with neural controllers [23]. In this
section, we propose an offline verification framework to simultaneously verify
the practical stability and safety of the system controlled by the projected
policy $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ in a deterministic
manner, based on MILP. A small adaption that can be used to verify the
asymptotic stability for $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ and
$\hat{\pi}^{\mathrm{im}}(\cdot)$ is provided at the end of this section.

The proposed verification procedure contains 3 steps. Different from [31, 32,
33] that directly verify asymptotic stability, for practical stability we need
to first verify the convergence of the closed-loop system to a neighborhood
containing the origin, and then verify the invariance of the neighborhood.
These will be formulated as two MILP problems. Finally, to enlarge the inner-
estimated region of attraction, which is a sub-level set of
$\hat{J}_{k-1}(\cdot)$, the third MILP problem will be formulated.

To verify the stability in any sub-level set
$\mathcal{B}_{r_{1}}=\\{x\in\mathcal{X}|\hat{J}_{k-1}(x)\leq
r_{1},\;r_{1}>0\\}$ that is contained in $X$, we can directly extend the
verifier in [31], in which a mixed-integer quadratic program is solved to test
a quadratic candidate Lyapunov function for PWA systems with input
constraints. In particular, after implementing Algorithm 2 and (IV-B), we have
the value function $\hat{J}_{k-1}(\cdot)$ and the explicit policy
$\hat{\pi}^{\mathrm{ex}}\left(\cdot,\omega^{*}\right)$ at our disposal. Then,
we formulate the following optimization problem:

$\displaystyle a^{*}_{1}=\max_{x,u}$
$\displaystyle\;\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,u))-\hat{J}_{k-1}(x)+c_{1}l(x,0)$
$\displaystyle\mathrm{s.t.}\;$ $\displaystyle
u=\hat{\pi}_{\mathrm{proj}}^{\mathrm{ex}}(x),\;r_{2}\leq\hat{J}_{k-1}(x)\leq
r_{1},$ (24)

where $c_{1}$ is a small positive parameter, and $r_{2}\in(0,r_{1})$. If
$a^{*}_{1}\leq 0$, we can conclude that the closed-loop system with
$\hat{\pi}_{\mathrm{proj}}^{\mathrm{ex}}(\cdot)$ is safe in
$\mathcal{B}_{r_{1}}$, and that any trajectories starting in
$\mathcal{B}_{r_{1}}$ will enter
$\mathcal{B}_{r_{2}}=\\{x\in\mathcal{X}|\hat{J}_{k-1}(x)\leq r_{2},\\}$ in
finite time. The formal result is included in Theorem 3. If $a^{*}_{1}>0$, we
can reduce the values of $r_{1}$ and $c_{1}$. In this way, the objective
function of problem (V-C) and the feasible region of $x$ become smaller, which
will make $a^{*}_{1}$ smaller.

As all functions in (V-C) are PWA and thus MILP representable, problem (V-C)
can be formulated as an MILP problem.

After the trajectories reach $\mathcal{B}_{r_{2}}$, we need to verify that
they will always stay in $\mathcal{B}_{r_{2}}$, i.e., we need to prove the
positive invariance of $\mathcal{B}_{r_{2}}$. This leads to the second MILP
problem:

$\displaystyle a^{*}_{2}=\max_{x,u}$
$\displaystyle\;\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,u))-c_{2}\hat{J}_{k-1}(x)-r_{2}+r_{2}c_{2}$
$\displaystyle\mathrm{s.t.}\;$ $\displaystyle
u=\hat{\pi}_{\mathrm{proj}}^{\mathrm{ex}}(x),\;0\leq\hat{J}_{k-1}(x)\leq
r_{2},$ (25)

where $c_{2}\in[0,1]$. One can directly take $c_{2}=0$ to minimize
$a^{*}_{2}$. Clearly, $a^{*}_{2}\leq 0$ implies the positive invariance of
$\mathcal{B}_{r_{2}}$, which will be proven in Theorem 3. If $a^{*}_{2}>0$,
similarly we can make $r_{2}$ smaller.

However, the stable and safe region $\mathcal{B}_{r_{1}}$ derived from (V-C)
and (V-C) is usually small, as $\mathcal{B}_{r_{1}}$ is a sub-level set of
$\hat{J}_{k-1}(\cdot)$. In particular, if the weights in $Q$ on different
states vary greatly, the resulting $\mathcal{B}_{r_{1}}$ will be rather narrow
and much smaller than the real region of attraction. It could also happen that
the evolution of the closed-loop system does not make $\hat{J}_{k-1}(\cdot)$
decrease at the beginning, but will drive the state into $\mathcal{B}_{r_{1}}$
in a finite number of time steps. Besides, in most cases we are interested in
the performance of a policy in a polyhedron (or a UoP), rather than a sub-
level set.

Therefore, we further develop the third optimization problem that evaluates
the range of trajectories of the closed-loop system in a finite number of time
steps for all initial states in a polyhedron $X_{\mathrm{in}}$ or a UoP of
interest. The problem is

$\displaystyle\mathrm{Check}\;\mathrm{if}\;x_{t}\in
X,\;t\\!=\\!1,...,N\\!-\\!1,\forall x_{0}\in X_{\mathrm{in}}$
$\displaystyle\mathrm{and}\;\mathrm{if}\;\hat{J}_{k-1}(x_{N})\\!\leq\\!r_{1},\forall
x_{0}\in X_{\mathrm{in}}$
$\displaystyle\mathrm{s.t.}\;u_{t}=\hat{\pi}_{\mathrm{proj}}^{\mathrm{ex}}(x_{t}),\;x_{t+1}=f_{\mathrm{PWA}}(x_{t},u_{t}),\;t=0,...,N-1,$
(26)

where $N$ is the number of time steps and $r_{1}$ is such that it makes
$a^{*}_{1}\leq 0$ in (V-C). If (V-C) returns “Yes”, we can conclude that for
any initial state in $X_{\mathrm{in}}$, the states of the closed-loop system
will satisfy the constraints from $t=0$ to $t=N-1$, and the final state
$x_{N}$ will reach the stable and safe region $\mathcal{B}_{r_{1}}$ computed
in (V-C). Similar to (V-C) and (V-C), (V-C) can be exactly expressed in an
MILP form. If $X_{\mathrm{in}}$ is a UoP, we also need some additional binary
variables to formulate the initial condition $x_{0}\in X_{\mathrm{in}}$ (see
Appendix A-C).

The integration of (V-C), (V-C), and (V-C) constitutes the proposed
verification framework, which computes the exact evolution of the closed-loop
system. It does not need any sampling or statistical testing procedure. The
effectiveness of the proposed verification framework is stated in the
following theorem.

###### Theorem 3

Consider the policy $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ and the
proposed verification framework consisting of (V-C), (V-C), and (V-C). If
$\hat{J}_{k-1}(0)=0$, $a^{*}_{1}\leq 0$, $a^{*}_{2}\leq 0$, and (V-C) returns
“Yes”, then the closed-loop system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ is safe in $X_{\mathrm{in}}$,
and any trajectory starting from $X_{\mathrm{in}}$ will approach
$\mathcal{B}_{r_{1}}$ in at most
$N+\lceil(r_{1}-r_{2})\hat{\gamma}/(c_{1}r_{2})\rceil$ time steps and stay in
$\mathcal{B}_{r_{1}}$ thereafter. Here $\hat{\gamma}$ is a positive constant
independent of the initial condition.

With the results in Theorems 2 and 3, we now make practical suggestions on how
to implement Algorithm 2 and the proposed verification framework.

After the explicit policy
$\hat{\pi}^{\mathrm{ex}}\left(\cdot,\omega^{*}\right)$ is obtained, we first
find $r_{2}$ that makes $a^{*}_{2}\leq 0$. Then, one can use (V-C) to compute
the safe and stable region and then use (V-C) to enlarge it. If $a^{*}_{2}\leq
0$ but $a^{*}_{1}>0$ whatever $r_{1}$ is chosen, one needs to apply (V-C) with
$r_{1}$ replaced by $r_{2}$. The cost is that one may have to choose a large
horizon $N$ since $r_{2}$ is small. The complexity of (V-C) grows
exponentially w.r.t. $N$. If $a^{*}_{1}>0$ and $a^{*}_{2}>0$ no matter what
$r_{1}$ and $r_{2}$ are chosen, one have to refine the learning process. To
this end, Theorem 2 implies that one can either (i) increase the number of
iterations by tightening the stopping condition, i.e., making $e_{\text{tole
}}$ smaller, or (ii) improve the approximation quality of function
approximators and restart Algorithm 2.

The proposed verification framework generalizes the verification methods in
[9, 31, 32, 33], and is an extension of [9, 32, 33] from linear systems to PWA
systems and from asymptotic stability to practical stability. Specifically, if
we remove (V-C) and (V-C) and fix $r_{2}=0$, the proposed method verifies
asymptotic stability, which is stronger than the properties in Theorem 3, and
it is then similar to the methods in [31, 32]. For the comparison with [33],
we note that to guarantee state constraint satisfaction, [33] needs to search
for a positively invariant set. In comparison, $X_{\mathrm{in}}$ in (V-C) is
not necessarily positively invariant. Besides, (V-C) includes the case when
the state constraint set $X$ and the searching region $X_{\mathrm{in}}$ are
UoPs. Nevertheless, we should mention that similarly to [9, 33], one should
refine the learning process as mentioned before and re-learn a policy if the
verification scheme fails, while the approaches in [31, 32] are more adaptive
since they actively learn a Lyapunov function that can be different from
$\hat{J}_{k-1}(\cdot)$.

In addition, (V-C) can be adjusted to verify the asymptotic stability and
safety of $\hat{\pi}^{\mathrm{im}}(\cdot)$. With $r_{2}=0$ and $c_{1}l(x,0)$
replaced by $c_{1}l(x,u)$, $a^{*}_{1}\leq 0$ tells that $\hat{J}_{k-1}(\cdot)$
is a control Lyapunov function, which further validates the asymptotic
stability and safety of $\hat{\pi}^{\mathrm{im}}(\cdot)$ in
$\mathcal{B}_{r_{2}}$.

##  VI Case study

We validate the proposed methods in two examples. The simulations are
conducted in MATLAB R2021a on an AMD Core R7-5800H CPU @3.20GHz machine. To
guarantee a fair comparison of CPU times, all algorithms are implemented in
MATLAB code. In particular, all MILP problems are solved by the MATLAB
function “intlinprog” while all nonlinear problems are solved by the MATLAB
function “fmincon” with the active-set algorithm.

###  VI-A Example 1: inverted pendulum with elastic walls

We consider a physical system, in which an active inverted pendulum is placed
between elastic walls to maintain the vertical position (see Fig. 1(a)). The
system parameters are chosen as
$m=1\;\mathrm{kg},\;l=1\;\mathrm{m},\;d_{1}=0.12\;\mathrm{m},\;,d_{2}=d_{3}=0.1\;\mathrm{m},\;k_{1}=300\;\mathrm{N/m},\;k_{2}=300\;\mathrm{N/m}$,
and $k_{3}=500\;\mathrm{N/m}$. The state of the system is
$x=[q\;\;\dot{q}]^{T}$. By linearizing the dynamics around the vertical
configuration $q=\dot{q}=0$ and discretizing the system with a sampling time
$0.05\;\mathrm{s}$, we obtain a discrete-time PWA model with 4 modes.

$\displaystyle A_{1}$ $\displaystyle=\left[\begin{array}[]{cc}1&0.05\\\
-29.5&1\end{array}\right],\quad B_{1}=\left[\begin{array}[]{c}0\\\
0.05\end{array}\right],\quad f_{1}=\left[\begin{array}[]{l}0\\\
-3.3\end{array}\right]$ $\displaystyle\mathcal{C}_{1}$
$\displaystyle=\\{(x,u)\mid[1\;\;0]x\leq-0.12\\}$ $\displaystyle A_{2}$
$\displaystyle=\left[\begin{array}[]{cc}1&0.05\\\
-14.5&1\end{array}\right],\quad B_{2}=\left[\begin{array}[]{c}0\\\
0.05\end{array}\right],\quad f_{2}=\left[\begin{array}[]{l}0\\\
-1.5\end{array}\right]$ $\displaystyle\mathcal{C}_{2}$
$\displaystyle=\\{(x,u)\mid-0.12\leq[1\;\;0]x\leq-0.1\\}$ $\displaystyle
A_{3}$ $\displaystyle=\left[\begin{array}[]{cc}1&0.05\\\
0.5&1\end{array}\right],\quad B_{3}=\left[\begin{array}[]{c}0\\\
0.05\end{array}\right],\quad f_{3}=\left[\begin{array}[]{l}0\\\
0\end{array}\right]$ $\displaystyle\mathcal{C}_{3}$
$\displaystyle=\\{(x,u)\mid-0.1\leq[1\;\;0]x\leq 0.1\\}$ $\displaystyle A_{4}$
$\displaystyle=\left[\begin{array}[]{cc}1&0.05\\\
-24.5&1\end{array}\right],\quad B_{4}=\left[\begin{array}[]{c}0\\\
0.05\end{array}\right],\quad f_{4}=\left[\begin{array}[]{l}0\\\
2.5\end{array}\right]$ $\displaystyle\mathcal{C}_{4}$
$\displaystyle=\\{(x,u)\mid[1\;\;0]x\geq 0.1\\}$

(a) An inverted pendulum with elastic walls.

(b) Closed-loop trajectories of the inverted pendulum. The shaded region is
$\bar{X}$. The curves with different colors represent the state trajectories
starting from different vertices. (c) The learning process of the ADP
algorithm with different options.

(d) Safe and stable regions verified by the verification framework. (e)
Statistical analysis of different methods regarding CPU times. D-max means the
difference of two max-affine functions.

(f) Closed-loop trajectories of the inverted pendulum with UoP state
constraints. The shaded region is $\bar{X}$. The curves with different colors
represent the state trajectories starting from different vertices.

Figure 1: Simulation results.

Due to safety considerations and the limited capability of the torque motor,
the inverted pendulum system is supposed to satisfy the constraints
$[-0.15\;\;-1]^{T}\leq x\leq[0.15\;\;1]^{T}$ and $-4\leq u\leq 4$. The
constrained infinite-horizon control problem is constructed with
$l(x,u)=||\mathrm{diag}([20\;\;1])x||_{\infty}+||u||_{\infty}$. A 61 $\times$
61 state data grid is constructed to cover the region
$\\{x|[-0.17\;\;-1.2]^{T}\leq x\leq[0.17\;\;1.2]^{T}\\}$. These training
points are exploited to train the actor and the critic. Both of them are ReLU
NNs that have two hidden layers of widths 8. Starting from a zero value
function and after 10 iterations, the closed-loop behavior of the system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ is illustrated in Fig. 1(b).
Fig. 1(b) plots the maximal stabilizable set $\bar{X}$ and the trajectories of
the closed-loop system (controlled by
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$) starting from some vertices
of $\bar{X}$. Although states starting from these vertices are in general the
most difficult to regulate to the origin, they converge rapidly under
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$. Meanwhile, we should mention
that the state constraints could be slightly violated: the trajectory depicted
by the purple curves violates the constraints by about 5 percent. To avoid
this, one can tighten the state constraints and restart the ADP algorithm. The
dashed purple curve describes the trajectory starting from the purple vertex
with 10 percent constraint tightening. We can observe that constraint
violation is avoided.

To carry out option (b) of Algorithm 2, we consider the linear subsystem
$(A_{3},\;B_{3})$ and readily get a linear stabilizing feedback law
$u_{0}=[-40\;\;-10]x$, which is used to initialize the value function of
option (b). The optimal value function $J^{*}(\cdot)$ is computed by the MPT3
toolbox [46] in 9.6 hours. The learning processes of different versions of
Algorithm 2 are compared in Fig. 1(c), in which the mean square errors between
$\hat{J}_{k-1}\left(x^{(i)},\theta_{k}\right)$ and $J^{*}(x^{(i)})$ are
depicted. One can see that the value function approximation in option (b) has
a much faster convergence rate than that in option (a), and option 2 results
in lower approximation errors than option 1.

Fig. 1(d) depicts the safe and stable regions that are analyzed by the
proposed verification framework. The blue region represents
$\mathcal{B}_{r_{1}}$, which is computed by (V-C) with
$r_{1}=18,\;c_{1}=0.1,\;r_{2}=3$. We note that for some states in
$\mathcal{B}_{r_{2}}$, which is colored red, the objective function in (V-C)
becomes positive, so the verification method in [32] fails. However, (V-C)
outputs a negative $a^{*}_{2}$ with $c_{2}=0.1$, which means that any
trajectory of the closed-loop system starting from $\mathcal{B}_{r_{1}}$ will
reach in the neighborhood $\mathcal{B}_{r_{2}}$ containing the origin in
finite number of time steps. Furthermore, the safe and stable region
$\mathcal{B}_{r_{1}}$ is enlarged to $X_{\mathrm{in}}$ (the yellow region) by
(V-C) with $N=3$.

We compare the CPU time of running the ADP-based controllers, hybrid MPC,
explicit MPC computed by the MPT3 toolbox, nonlinear MPC that uses the
nonlinear model, and action-governor RL in [13]. With an MPC horizon 8, the
implementation of hybrid MPC requires to solve an MILP problem that contains
94 continuous variables, 28 binary variables, 328 inequality constraints, and
21 equality constraints. In comparison, for $\hat{\pi}^{\mathrm{im}}(\cdot)$
with the value function approximator chosen as the ReLU NN (or the difference
of two max-affine functions that has 10 and 3 terms in the first and second
max blocks), the resulting MILP problem has 29(15) continuous variables, 20(7)
binary variables, 94(73) inequality constraints, and 3(4) equality
constraints. We only compare the time of computation that needs the knowledge
of the current state, including constructing and solving optimization
problems. For statistical analysis, we randomly select 100 initial states and
run the system in 50 time steps222The number of time steps is set based on the
fact that in our numerical experiments most scenarios converge to the origin
after this number of time steps.. The average computation time per time step
of the different methods is shown in Fig. 1(e). The computation of
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ requires the least amount of
CPU time, which is almost less than $10^{-5}$ per time step. Besides,
$\hat{\pi}^{\mathrm{im}}(\cdot)$ performs better than hybrid MPC, explicit
MPC, and action-governor RL, regarding online computation time, both when
using the ReLU NN or the difference of two max-affine functions for value
function approximation.

We further consider the case when the state constraints are a UoP, namely,
$x\in X^{(1)}\cup X^{(2)}$ with $X^{(1)}=\\{x|[-0.15\;\;-1]^{T}\leq
x\leq[0.15\;\;1]^{T}\\}$ and $X^{(2)}=\\{x|[-0.08\;\;-1.5]^{T}\leq
x\leq[0.08\;\;1.5]^{T}\\}$. This means a higher angular velocity is allowed
when the pendulum is close to its vertical configuration. After implementing
option 2(a) of Algorithm 2 in 10 iterations, $\bar{X}$ is depicted in Fig.
1(f). The trajectories of the closed-loop system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ starting from some vertices
of $\bar{X}$ are plotted in Fig. 1(f), from which one can find that the
proposed scheme is still valid even when the state constraints are a UoP.

###  VI-B Example 2: centralized adaptive cruise control

We consider an adaptive cruise control problem in which we drive three
vehicles (followers) to follow a leading vehicle (leader) in a highway
environment. The control objective is to keep a desired distance $d=20$ m
between each two adjacent vehicles.

Figure 2: Adaptive cruise control set up.

_Model:_ The differential equation for the velocity $v_{i}(t),\;i=1,2,3$ of
each following vehicle is $m\dot{v_{i}}(t)+cv^{2}_{i}(t)+\mu mg=bf_{i}(t)$,
where $bf_{i}(t)$ is the input traction/brake force, and where
$c=0.5\;\mathrm{kg/m},\;m=800\;\mathrm{kg},\;\mu=0.01,\;b=3700\;\mathrm{N},\;g=9.8\;\mathrm{m/s^{2}}$.
The maximum allowable velocity for followers is
$v_{\mathrm{max}}=35\;\mathrm{m/s}$.

Following [47], a least-squares PWA approximation of the nonlinear damping
force $V(v)=cv^{2}$ results in the continuous-time PWA dynamics with 2 modes
for each vehicle

$m\dot{v_{i}}(t)+c_{j}v_{i}(t)+a_{j}=bf_{i}(t),\;\mathrm{if}\;v_{i}(t)\in[\alpha_{j-1},\alpha_{j}),j=1,2$
(27)

where
$\alpha_{0}=0,\;\alpha_{1}=v_{\mathrm{max}}/2,\;\alpha_{2}=v_{\mathrm{max}}$,
$c_{1}=3cv_{\mathrm{max}}/8,\;c_{2}=13cv_{\mathrm{max}}/8$, $a_{1}=\mu mg$,
and $a_{2}=-5cv_{\mathrm{max}}^{2}/8+\mu mg$. The leader is assumed to have a
constant velocity $v_{\mathrm{r}}=20\;\mathrm{m/s}$. By defining a new state
$x=[p_{1}-p_{\mathrm{r}}-d\quad v_{1}-v_{\mathrm{r}}\quad p_{2}-p_{1}-d\quad
v_{2}-v_{\mathrm{r}}\quad p_{3}-p_{2}-d\quad v_{3}-v_{\mathrm{r}}]^{T}$, where
$p_{1},\;p_{2},\;p_{3}$ and $p_{\mathrm{r}}$ are the positions of the three
followers and the leader, respectively, the control goal is then making $x$
zero. The dynamics of $x$ thereby have $2^{3}$ different modes. Noticing that
when $x$ is zero, $v_{i}\in[\alpha_{1},\alpha_{2}),i=1,2,3$, so we introduce
an input transformation $u=[f_{1}-\frac{v_{\mathrm{r}}c_{2}+a_{2}}{b}\quad
f_{2}-\frac{v_{\mathrm{r}}c_{2}+a_{2}}{b}\quad
f_{3}-\frac{v_{\mathrm{r}}c_{2}+a_{2}}{b}]^{T}$ to make the origin of $x$-$u$
space the equilibrium of the state equations for $x$. By discretizing the
differential equations for $x$ with the sampling time 1s, we get a discrete-
time PWA model with 6 states, 3 inputs, and 8 modes.

As for constraints, we consider limitations on the distance between each two
adjacent vehicles, the velocity, and the traction/brake force input. In
particular, we require that $10\;\mathrm{m}\leq|p_{i}-p_{i-1}|\leq
30\;\mathrm{m}$, $5\;\mathrm{m/s}\leq v_{i}\leq 35\;\mathrm{m/s}$, and
$|f_{i}|\leq 1\;\mathrm{N}$, where $i=1,2,3$ and $p_{0}=p_{\mathrm{r}}$. All
of these constraints can be readily converted into polyhedral constraints on
$x$ and $u$.

The stage cost is chosen as
$l(x,u)=\left\|\operatorname{diag}\left(\left[\begin{array}[]{llllll}1&0.5&1&0.5&1&0.5\end{array}\right]\right)x\right\|_{\infty}+\left\|\operatorname{diag}\left(\left[\begin{array}[]{lll}0.1&0.1&0.1\end{array}\right]\right)u\right\|_{\infty}$.

_Setups:_ We uniformly randomly sample 25000 state points from the state space
$\\{x\in\mathbb{R}^{6}\;|\;[-12\;\;-18\;\;-12\;\;-18\;\;-12\;\;-18]^{T}\leq
x\leq[12\;\;18\;\;12\;\;18\;\;12\;\;18]^{T}\\}$. We compare the performance of
different controllers, including the proposed ADP controllers with different
parameter settings, hybrid MPC with different prediction horizons, and the
safety-filter-based RL in [12]. The performance metrics include the averaged
CPU time per time step, the cumulative stage cost over 15 time steps, and the
percentage that the system trajectory satisfies all constraints. To test the
performance of controllers, we generate 100 initial state scenarios that make
the hybrid MPC with horizon 4 recursively feasible and stabilizing. Since
controllable sets are difficult to compute in this example, we start from a
zero value function and adopt option 2 of Algorithm 2 to get
$\hat{J}_{k-1}\left(\cdot,\theta_{k-1}\right)$. All PWA function approximators
used in the proposed method and the safety-filter-based RL are ReLU NNs with 2
hidden layers.

_Results:_ Table 1 shows the results. In Table 1, $N$ represents the
prediction horizon chosen in hybrid MPC, $M$ refers to the number of units in
each hidden layer, and the last column lists the number of the initial state
scenarios that make the system trajectory satisfy all constraints. As
expected, MPC provides the best performance regarding the total cost and
constraint satisfaction, but it needs too much CPU time for computation.
Compared with MPC, $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ and the method in
[12] can reduce the online computation time, especially for
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$. With the increase of the number of
units, the performance of $\hat{\pi}^{\mathrm{im}}$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ regarding the total cost becomes
better and close to that of MPC. As both $\hat{\pi}^{\mathrm{im}}$ and
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ have acceptable performance
regarding the total cost and constraint satisfaction, using
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ is preferable due to its significant
reduction of CPU time. Besides, Table 1 also indicates that with the same
number of units, the RL method in [12] can be computed faster than
$\hat{\pi}^{\mathrm{im}}$, but the performance of $\hat{\pi}^{\mathrm{im}}$
and $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ is superior to that of the RL
method in [12] in terms of both cumulative cost and constraint violations.

TABLE I: Performance of different controllers Methods |  CPU time (mean/max) |  Total cost |  Constraint satisfaction
---|---|---|---
MPC with $N=4$ |  0.1592 / 0.4593 |  21.7664 |  100%
MPC with $N=5$ |  0.2850 / 0.7703 |  21.6653 |  100%
MPC with $N=6$ |  0.4630 / 1.4619 |  21.5997 |  100%
$\hat{\pi}^{\mathrm{im}}$ with $M=15$ |  0.0605 / 0.1036 |  23.0097 |  98%
$\hat{\pi}^{\mathrm{im}}$ with $M=20$ |  0.0916 / 0.1621 |  22.6613 |  100%
$\hat{\pi}^{\mathrm{im}}$ with $M=25$ |  0.1372 / 0.2991 |  22.3042 |  100%
$\hat{\pi}^{\mathrm{im}}$ with $M=30$ |  0.5487 / 1.3904 |  22.0716 |  100%
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ with $M=15$ |  0.0008 / 0.0028 |  23.4043 |  99%
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ with $M=20$ |  0.0007 / 0.0027 |  22.9568 |  99%
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ with $M=25$ |  0.0007 / 0.0027 |  22.7022 |  99%
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}$ with $M=30$ |  0.0008 / 0.0030 |  22.7562 |  99%
RL in [12] with $M=30$,$N=5$ |  0.0809 / 0.3570 |  29.7561 |  85%

##  VII CONCLUSIONS AND FUTURE WORK

We have proposed an ADP control scheme to deal with infinite-horizon optimal
control of PWA systems subject to linear and union-of-polyhedra (UoP)
constraints, based on MILP. With carefully designed PWA penalty functions, the
probably non-convex UoP constraints during the learning process are removed
while the PWA properties of the value functions are maintained. We have
formally analyzed the PWA properties and continuity of the value function, as
well as the closed-loop stability and safety under the approximation errors.
We have also designed an offline verification tool to make the proposed method
reliable. Simulation results show that the ADP-based policies are near-
optimal, and require much less online computational effort than conventional
hybrid MPC. Topics for future work include eliminating the reliance on value
function approximation, exploring more efficient sampling strategies, and
considering multi-time-step constraints.

##  Appendix A Technical proofs

###  A-A Proof of Theorem 1

The proof contains three parts.

####  A-A1 Convergence of the VI sequence in option 1

_Option 1(a):_ For the nonlinear parametric problems in (9), $U$ is closed,
and the set $\left\\{{u\in
U|\;J_{k}^{\mathrm{soft}}(x)=l(x,u)+P(x,X)+J_{k-1}^{\mathrm{soft}}\left({{f_{{\text{PWA}}}}(x,u)}\right)}\right\\}$
is nonempty for any $x$ in $\mathcal{X}$. Then, according to [48, Theorem
4.2.1], it can be recursively proven that
$J_{k}^{\mathrm{soft}}(\cdot),\;k=0,1,\dots$ is lower-semicontinuous on
$\mathcal{X}$. Hence, all sub-level sets
$\\{x\in\mathcal{X}|J_{k}^{\mathrm{soft}}(x)\leq\lambda\\}$ with
$\lambda\in\mathbb{R}$ are closed. As $l_{p}(\cdot,\cdot)$ and $f_{\text{PWA
}}(\cdot,\cdot)$ are continuous on $\mathcal{X}\times\mathcal{U}$, the set
${U_{k}}\left({x,\lambda}\right)=\\{u\in
U|\;{l}\left({x,u}\right)+P(x,X)+J_{k-1}^{\mathrm{soft}}\left({{f_{{\text{PWA
}}}}(x,u)}\right)\leq\lambda\\}$ is closed and thus compact for all
$x\in\mathcal{X}$, $\lambda\in\mathbb{R}$, and $k\geq 1$. As a result, the
compactness assumption in [38] is satisfied and
$\\{J_{k}^{\mathrm{soft}}(\cdot)\\}^{\infty}_{k=0}$ converges point-wise to
$J^{\mathrm{soft*}}(\cdot)$, according to [38, Proposition 2].

_Option 1(b):_ In this case, we define
$J^{\mathrm{\infty}}_{0}(x)=\begin{cases}J_{\mathrm{CL}}(x),&x\in
X_{\mathrm{CI}}\\\ \infty&x\notin X_{\mathrm{CI}}\end{cases}$. The point-wise
convergence of the VI
$J^{\mathrm{\infty}}_{k}(x)=\Gamma_{\mathrm{p,1}}J^{\mathrm{\infty}}_{k-1}(x)$
to $J^{\mathrm{soft}*}(x)$ is always guaranteed by [38, Proposition 2] since
$J^{\mathrm{\infty}}_{0}(x)\geq J^{\mathrm{soft}*}(x),\;\forall
x\in\mathcal{X}$. By applying the monotonicity of the Bellman operator
$\Gamma_{\mathrm{p,1}}$ [37], we get the convergence of the VI sequence
$\\{J_{k}^{\mathrm{soft}}(\cdot)\\}^{\infty}_{k=0}$ to
$J^{\mathrm{soft*}}(\cdot)$ in option 1(b).

####  A-A2 Continuity and PWA property of the value function in option 1

By recursively iterating (9), it is observed that $J_{k}^{\mathrm{soft}}(x)$
can be derived via the following batch approach:

$\displaystyle J_{k}^{\mathrm{soft}}(x)=$
$\displaystyle\min_{u_{0},\ldots,u_{k-1},\;{x_{0}},\ldots,{x_{k}}}\sum_{i=0}^{k-1}l_{\mathrm{p}}\left(x_{i},u_{i}\right)+J_{0}^{\mathrm{soft}}(x_{k})$
s.t. $\displaystyle x_{i+1}=f_{\text{PWA }}\left(x_{i},u_{i}\right),x_{0}=x,$
$\displaystyle u_{i}\in U,i=0,\ldots,k-1.$ (28)

_Option 1(a):_ The proof of continuity follows from the proof of [30,
Corollary 17.2], because the objective function in (A-A2) is continuous and
there is no state constraint.

Now, we prove the PWA property of $J_{k}^{\mathrm{soft}}(\cdot)$. We denote by
$\left\\{\Xi_{j}\right\\}_{j=1}^{r_{0}^{k}}$ the set of all possible switching
sequences for $P(\cdot,\cdot)$, and by $\Xi_{j}^{(i)}$ the $i$-th element of
the sequence $\Xi_{j}$. Namely, $\Xi_{j}^{(i)}=\vartheta$ if
$P(x_{i},X)=p\max(0,(E^{(\vartheta)}_{X})_{1\cdot}x-(g^{(\vartheta)}_{X})_{1},\dots,(E^{(\vartheta)}_{X})_{m^{(\vartheta)}_{x}\cdot}x-(g^{(\vartheta)}_{X})_{m^{(\vartheta)}_{x}})$.
Here, without loss of generality we assume the first type of $P(\cdot,\cdot)$
in (III) is chosen. We notice that for a fixed switching sequence of
$P(\cdot,\cdot)$, $P(\cdot,\cdot)$ can be converted into the linear inequality
constraints $s^{(\vartheta)}\geq 0$ and $s^{(\vartheta)}\geq
p(E^{(\vartheta)}_{X})_{\ell\cdot}x-p(g^{(\vartheta)}_{X})_{\ell}$,
$\ell=1,...,m^{(\vartheta)}_{x}$ by introducing the slack variable
$s^{(\vartheta)}$. Then, for a fixed switching sequence of the PWA dynamics
and a fixed switching sequence of $P(\cdot,\cdot)$, problem (A-A2) reduces to
a finite-horizon optimal control problem for a linear time-varying system with
linear time-varying constraints, and thereby provides a PWA value function. As
a result, $J_{k}^{\mathrm{soft}}(\cdot)$ is thus determined by comparing the
value functions for all switching sequences and selecting the smallest one, so
it is a PWA function in $\mathcal{X}$.

_Option 1(b):_ We firstly need to prove that the initial value function
$J_{0}^{\mathrm{soft}}(\cdot)$ is continuous and PWA on $\mathcal{X}$. The
continuity and PWA property of $J_{0}^{\mathrm{soft}}(\cdot)$ on
$X_{\mathrm{CI}}$ is clear according to its definition in (7). For problem
(8), there always exists a continuous optimizer $\bar{z}(\cdot)$ on
$\mathcal{X}$, according to [30, Theorem 6.5]. Together with the continuity of
$J_{\mathrm{CL}}(\cdot)$ and $P(\cdot,X_{\mathrm{CI}})$ in their domains, this
ensures that $J_{0}^{\mathrm{soft}}(\cdot)$ is continuous on
$\mathcal{X}\backslash X_{\mathrm{CI}}$. The continuity of
$J_{0}^{\mathrm{soft}}(\cdot)$ on the whole state space $\mathcal{X}$ thus
follows from the fact that $\mathop{\lim}\limits_{x\to
y,x\in\mathcal{X}/X_{\mathrm{CI}}}{J^{\mathrm{soft}}_{0}}\left(x\right)={J_{\mathrm{CL}}}\left(y\right)$
holds for any $y\in\partial X_{\mathrm{CI}}$.

Besides, it is straightforward to prove that $J_{0}^{\mathrm{soft}}(\cdot)$ is
also PWA outside $X_{\mathrm{CI}}$, since the optimizer $\bar{z}(\cdot)$ of
the mp-LP, $J_{\mathrm{CL}}(\cdot)$, and $P(\cdot,X_{\mathrm{CI}})$ are all
PWA [30, Theorem 6.5]. Due to the continuity and PWA property of
$J_{0}^{\mathrm{soft}}(\cdot)$, the remainder of the proof of the continuity
and PWA property of $J_{k}^{\mathrm{soft}}(\cdot)$ is similar to that in
option 1(a). Thus, we have completed the proof of the statements of Theorem 1
in option 1.

####  A-A3 Equivalence of the VI sequences for options 1 and 2

In option 2, by iterating
$J_{k}^{\mathrm{soft}}(x)=\Gamma_{\mathrm{p,2}}J_{k-1}^{\mathrm{soft}}(x),\;\forall
x\in\mathcal{X}$ from $k$ to 0, we can get the expression of
$J_{k}^{\mathrm{soft}}(\cdot)$ via the batch approach:

$\displaystyle{J_{k}^{\mathrm{soft}}}\left({{x}}\right)=\mathrm{the\;optimal\;value\;of}\;\eqref{finite}-P(x,X).$
(29)

From (29), we notice that the difference between the value functions in
options 1 and 2 at the same iteration is $P(x,X)$, which is always continuous
in $x$ and equals zero if $x\in X$. Combining (29) with the first and second
parts of the proof proves the statements of Theorem 1.

###  A-B Proof of Lemma 1

The proof can be constructed by induction. Noticing that $X_{0}=X$ is a UoP,
we assume that $X_{k-1}$ is a UoP in the form of
$X_{k-1}=\bigcup_{i=1}^{r_{k-1}}X^{(i)}_{k-1}$. Then, the backward-reachable
set $\mathrm{Pre}(X_{k-1})$ can be computed as
$\mathrm{Pre}(X_{k-1})=\bigcup_{i=1}^{r_{k-1}}\mathrm{Pre}(X^{(i)}_{k-1})$.
This results in

$\displaystyle X_{k}$ $\displaystyle=\mathrm{Pre}(X_{k-1})\cap
X=\bigcup_{j=1}^{r_{0}}\mathrm{Pre}(X_{k-1})\cap X^{(j)}$
$\displaystyle=\bigcup_{j=1}^{r_{0}}\left(\bigcup_{i=1}^{r_{k-1}}\mathrm{Pre}(X^{(i)}_{k-1})\right)\cap
X^{(j)}$
$\displaystyle=\bigcup_{j=1}^{r_{0}}\bigcup_{i=1}^{r_{k-1}}\mathrm{Pre}(X^{(i)}_{k-1})\cap
X^{(j)}$ (30)

In the last line of (A-B), $\mathrm{Pre}(X^{(i)}_{k-1})\cap X^{(j)}$ is a UoP
for each $i$ and $j$ because $X^{(i)}_{k-1}$ and $X^{(j)}$ are polyhedra. As a
result, $X_{k}$ is a UoP. This proves the statement for $k$. Finally, by
induction, we complete the proof of Lemma 1.

###  A-C PWA functions and their mixed-integer formulations

It has been shown that a ReLU NN can represent a PWA function [49]. A ReLU NN
with $L$ hidden layers can be written by

$f_{\mathrm{NN}}(x,\theta)=[f_{L+1}\circ f_{\mathrm{ReLU}}\circ
f_{L}\circ\cdots\circ f_{\mathrm{ReLU}}\circ f_{1}](x)$ (31)

where each hidden layer contains an affine map
$f_{\ell}\left(\kappa_{\ell-1}\right)=W_{\ell}\kappa_{\ell-1}+b_{\ell}$,
followed by a nonlinear map $f_{\mathrm{ReLU}}(f_{\ell})=\max(0,f_{\ell})$,
$\ell=1,2,\dots,L$. Here, $M_{\ell}$ is the width of each layer, referring to
the number of units in each layer, $\kappa_{\ell-1}\in\mathbb{R}^{M_{\ell-1}}$
is the output of the previous layer with $\kappa_{0}$ the input of the
network, $W_{\ell}\in\mathbb{R}^{M_{\ell}\times M_{\ell-1}}$ and
$b_{\ell}\in\mathbb{R}^{M_{\ell}}$ are the weights and biases, condensed in
the parameter $\theta$.

Alternatively, [50] has revealed that every continuous PWA scalar function can
be written as the difference of two max-affine functions:

${f_{{\text{D}}-\max}}\left({x,\theta}\right)=\max({W^{(1)}}x+{b^{(1)}})-\max({W^{(2)}}x+{b^{(2)}})$
(32)

where $\max(a)$ selects the maximum element of the vector $a$,
$W^{(1)}\in\mathbb{R}^{m_{1}\times n_{x}}$, and where
$W^{(2)}\in\mathbb{R}^{m_{2}\times n_{x}}$, $b^{(1)}\in\mathbb{R}^{m_{1}}$ and
$b^{(2)}\in\mathbb{R}^{m_{2}}$ constitute the parameter $\theta$. The
regression method that uses the difference of two max-affine functions can
thus be implemented by fixing the numbers $m_{1},\;m_{2}$ of affine functions
in the two max blocks and then updating $\theta$.

_Mixed-integer formulation of PWA function approximators:_ Consider the ReLU
NN in (31), the outputs $\kappa_{\ell+1},\;\kappa_{\ell},\;\ell=0,1,\dots,L$
of any two adjacent layers satisfy
$\kappa_{\ell+1}=f_{\mathrm{ReLU}}(W_{\ell}\kappa_{\ell}+b_{\ell})$. It has
been shown in [51] that the ReLU function has the mixed-integer linear
equivalent representation:

$\displaystyle\mathrm{gr}_{\mathrm{MILP}}(f_{\mathrm{ReLU}})$ $\displaystyle=$
$\displaystyle\left\\{(\kappa_{\ell},\kappa_{\ell+1})\\!\left|\\!\begin{array}[]{l}\kappa_{\ell+1}\geq
W_{\ell}\kappa_{\ell}+b_{\ell}\\\ \kappa_{\ell+1}\leq
W_{\ell}\kappa_{\ell}+b_{\ell}\\!-\\!\operatorname{diag}\left(\underline{M}_{\ell}\right)\left({{\mathbf{1}}_{M\times
1}}\\!-\\!\delta_{\ell}\right)\\\ \kappa_{\ell+1}\geq 0\\\
\kappa_{\ell+1}\leq\operatorname{diag}\left(\bar{M}_{\ell}\right)\delta_{\ell}\end{array}\right.\right\\}$

where $\underline{M}_{\ell}$ and $\bar{M}_{\ell}$ are the lower and upper
bounds of $W_{\ell}\kappa_{\ell}+b_{\ell}$, and
$\delta_{\ell}\in\\{0,1\\}^{M_{\ell}}$ is a binary variable. These bounds
$\underline{M}_{\ell}$ and $\bar{M}_{\ell}$ can be computed by interval bound
propagation or linear programming if the input of the network (the state of
system (1)) is in a compact set [52]. Finally, by applying
$\mathrm{gr}_{\mathrm{MILP}}(f_{\mathrm{ReLU}})$ to all the layers, one can
derive the mixed-integer linear equivalent representation
$\mathrm{gr}_{\mathrm{MILP}}(f_{\mathrm{NN}})$ of the NN (31).

For the approximator (32), since (32) will be contained in the minimization
problems (12) and (14), the first max function can be simply converted into
$m_{1}$ different linear inequalities. The second max function, according to
[51], leads to a mixed-integer formulation. In summary, we have

$\displaystyle\mathrm{gr}_{\mathrm{MILP}}(\max({W^{(2)}}x+{b^{(2)}}))$
$\displaystyle=$
$\displaystyle\left\\{(x,\varpi_{2})\left|\begin{array}[]{l}\varpi_{2}\geq
W^{(2)}_{i,}x+b^{(2)}_{i},\;\mathrm{for}\;i=1,\dots,m_{2}\\\ \varpi_{2}\leq
W^{(2)}_{i,}x+b^{(2)}_{i}+(1-\chi_{i})(\bar{m}_{-i}-\underline{m}_{i}),\\\
\quad\quad\quad\mathrm{for}\;i=1,\dots,m_{2}\\\
\sum\limits_{i=1}^{m_{2}}{{\chi_{i}}}=1\end{array}\right.\right\\}$ (37)

where $\varpi_{2}$ and $\chi=[\chi_{1}\;\cdots\;\chi_{m_{2}}]^{T}$ are
additionally introduced continuous and binary variables, respectively,
$\underline{m}_{i}$ and $\bar{m}_{i}$ are the lower and upper bounds of
$W^{(2)}_{i,}x+b^{(2)}_{i}$ in $X$, and $\bar{m}_{-i}=\max_{j\neq
i}\bar{m}_{j}$.

In addition, some other candidates for PWA approximation are, e.g., max-out
NNs [49], which inherently have equivalent mixed-integer formulations
according to (A-C), and the PWA regression algorithm in [53], in which a
mixed-integer encoding of the algorithm is also provided.

_Mixed-logical dynamical formulation of PWA systems:_ In the traditional MPC
framework for PWA systems [35], the PWA model is usually converted into an
equivalent mixed-logical dynamical form, so that the MPC problem is formulated
as a solvable MILP problem.

By specifying the partition $\left\\{\mathcal{C}_{i}\right\\}_{i=1}^{s}$ as
$\\{{\mathcal{C}}_{i}\\}_{i=1}^{s}\triangleq\\{\\{\left[\begin{array}[]{ll}x^{T}&u^{T}\end{array}\right]^{T}\in\mathbb{R}^{n_{x}+n_{u}}:{\Psi^{(x)}_{i}}x+{\Psi^{(u)}_{i}}u\leq{\psi_{i}}\\}\\}^{s}_{i=1}$,
according to [30], (1) can always be described as a mixed-logical dynamical
(MLD) system. An equivalent representation [54] is given by

$\displaystyle\mathrm{gr}_{\mathrm{MILP}}({f_{{\text{PWA}}}})$
$\displaystyle=$ $\displaystyle\left\\{\left(\left[\begin{array}[]{l}x\\\
u\end{array}\right],y\right)\left|\begin{array}[]{l}y=\sum\limits_{i=1}^{s}{{z_{i}}},\;\sum\limits_{i=1}^{s}{{\Delta_{i}}}=1\hfill\\\
{\text{for}}\;i=1,2,...,s:\hfill\\\ \left\\{\begin{gathered}{z_{i}}\leq
M{\Delta_{i}}\hfill\\\ {z_{i}}\geq m{\Delta_{i}}\hfill\\\
{z_{i}}\leq{A_{i}}x+{B_{i}}u+{f_{i}}-m\left({1-{\Delta_{i}}}\right)\hfill\\\
{z_{i}}\geq{A_{i}}x+{B_{i}}u+{f_{i}}-M\left({1-{\Delta_{i}}}\right)\hfill\\\
{\Psi^{(x)}_{i}}x+{\Psi^{(u)}_{i}}u\leq{\psi_{i}}+M_{i}^{*}[1-{\Delta_{i}}]\hfill\\\
\end{gathered}\right.\hfill\\\ \end{array}\right.\right\\}$ (49)

where $M$, $m$, and $M_{i}^{*}$ are some bounds of affine functions of $x$ and
$u$ over $X\times U$, and $\Delta_{i}\in\\{0,1\\},i=1,2,\dots,s$ are binary
variables. The detailed definitions of the bounds $M$, $m$, and $M_{i}^{*}$
can be found in [54].

The last inequality in (A-C) is used to formulate the logic
$[x^{T}\;\;u^{T}]^{T}\in\cup^{s}_{i=1}\mathcal{C}_{i}$. This can be
generalized to convert any UoP constraint into a mixed-integer constraint,
such as $x_{0}\in X_{\mathrm{in}}$ in (V-C).

_Mixed-integer formulation of penalty functions in ( III):_ As (12) and (14)
are all minimization problems, the maximum operation in $P(\cdot,\cdot)$ can
be eliminated by introducing some linear inequalities and continuous
variables. On the other hand, the minumum operation in $P(\cdot,\cdot)$
results in a mixed-integer formulation similar to (A-C).

###  A-D Proof of Lemma 2

Firstly, the result in (i) is straightforward, based on the fact that both
$J_{k}^{\mathrm{soft}}(\cdot)$ and $l(\cdot,0)$ are continuous PWA functions
on $\bar{X}$ and equal to zero iff $x=0$.

For (ii), we will prove by contradiction the existence of a finite $\bar{k}$
such that
$J_{\bar{k}}^{\mathrm{soft}}(x)-J_{\bar{k}-1}^{\mathrm{soft}}(x)\leq\beta
l(x,0),\;\forall x\in\bar{X}$ holds in option 1. Suppose that it does not
hold, then for any $k>0$, there exists a $\bar{x}\in\bar{X}$ such that
$J_{k}^{\mathrm{soft}}(\bar{x})-J_{k-1}^{\mathrm{soft}}(\bar{x})\geq
l(\bar{x},0)$. According to the contraction property of the Bellman operator
[37], we have
$J_{i}^{\mathrm{soft}}(\bar{x})-J_{i-1}^{\mathrm{soft}}(\bar{x})\geq
l(\bar{x},0),\;\forall i\in\\{1,...,k\\}$. Summing over $i$ yields
$J_{k}^{\mathrm{soft}}(\bar{x})\geq kl(\bar{x},0)$. Letting $k\to\infty$ makes
$J_{\infty}^{\mathrm{soft}}(\bar{x})$ contradict (i) of Lemma 2. Moreover,
again using the contraction property of $\Gamma_{\mathrm{p,1}}$ leads to the
conclusion in (ii). Finally, (29) tells that the difference of
$J_{\bar{k}}^{\mathrm{soft}}(x)$ in options 1 and 2 is $P(x,X)$, which is
independent of $k$. Therefore, (18) also holds for option 2 of Algorithm 1.

###  A-E Proofs of Theorem 2 and Corollary 1

For each $k\geq\bar{k}$, we define the policy $\pi_{k}(\cdot)$ as (one of) the
optimizer(s) of
$\Gamma_{\mathrm{p},\alpha}J_{k-1}^{\mathrm{soft}}(\cdot),\;\alpha=1\;\mathrm{or}\;2$.
For every $x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\subseteq X$,
according to (2),

$J_{k}^{\mathrm{soft}}(x)=\Gamma_{\mathrm{p},\alpha}J_{k-1}^{\mathrm{soft}}(x)\geq
l(x,\pi_{k}(x))+J_{k-1}^{\mathrm{soft}}(f_{\mathrm{PWA}}(x,\pi_{k}(x)))$ (50)

holds in both options. Together with (18), (50) yields

$\displaystyle
J_{k-1}^{\mathrm{soft}}(f_{\mathrm{PWA}}(x,\pi_{k}(x)))-J_{k-1}^{\mathrm{soft}}(x)$
$\displaystyle\leq-(1-\beta)l(x,\pi_{k}(x))$
$\displaystyle\leq-(1-\beta)l(x,0),$ (51)

which means that for every $x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)$,
the policy $\pi_{k}(\cdot)$ will make
$f_{\mathrm{PWA}}(x,\pi_{k}(x))\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)$.

Now, we show that $J^{\mathrm{soft}}_{k-1}(\cdot)$ and $\hat{J}_{k-1}(\cdot)$
are Lyapunov functions for the system with $\hat{\pi}^{\mathrm{im}}(\cdot)$.
In particular, for any
$x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$,
if C1 and (19) hold, we have

$\displaystyle
l(x,\hat{\pi}^{\mathrm{im}}(x))+\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{im}}(x)))$
$\displaystyle\leq$ $\displaystyle
l(x,\pi_{k}(x))+\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,\pi_{k}(x)))$
$\displaystyle\leq$ $\displaystyle
l(x,\pi_{k}(x))+(1+\zeta)J^{\mathrm{soft}}_{k-1}(f_{\mathrm{PWA}}(x,\pi_{k}(x)))$
$\displaystyle\leq$
$\displaystyle(1+\zeta)J^{\mathrm{soft}}_{k-1}(x)+\left(1-(1+\zeta)(1-\beta)\right)l(x,\pi_{k}(x))$
$\displaystyle\leq$ $\displaystyle\hat{J}_{k-1}(x)+2\zeta
J^{\mathrm{soft}}_{k-1}(x)+\left(1-(1+\zeta)(1-\beta)\right)l(x,\pi_{k}(x)).$
(52)

In (A-E), the first inequality is true since $\hat{\pi}^{\mathrm{im}}(\cdot)$
is an optimizer of problem (14); the second and the last inequalities hold
because $x$ and $f_{\mathrm{PWA}}(x,\pi_{k}(x))$ all in $\Omega$, in which C1
holds; and the third inequality is correct owning to (A-E). Since
$1-(1+\zeta)(1-\beta)<0$, considering Lemma 2, (A-E) results in

$\displaystyle\hat{J}_{k-1}(x^{+})-\hat{J}_{k-1}(x)\leq\left(2\zeta\gamma-(1+\zeta)(1-\beta)\right)l(x,0),$
$\displaystyle\mathrm{with}\;x^{+}=f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{im}}(x)).$
(53)

The right-hand side of (A-E) is strictly negative except for $x=0$ according
to (19). Therefore, we can conclude that $\forall
x\in\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$,
we have
$f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{im}}(x))\in\mathcal{B}(\hat{J}_{k-1},\Omega)$.
Meanwhile, condition C1 also implies that $\hat{J}_{k-1}(x)>0,\;\forall
x\in(\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega))\backslash\\{0\\}$
and $\hat{J}_{k-1}(0)=0$. This shows that $J^{\mathrm{soft}}_{k-1}(\cdot)$ is
a Lyapunov function for the closed-loop system
$x_{t+1}=f_{\mathrm{PWA}}(x_{t},\hat{\pi}^{\mathrm{im}}(x_{t}))$.

Similarly to (A-E) and (A-E), we can get the following inequality for
$J^{\mathrm{soft}}_{k-1}(\cdot)$:

$\displaystyle
l(x,\hat{\pi}^{\mathrm{im}}(x))+J^{\mathrm{soft}}_{k-1}(f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{im}}(x)))$
$\displaystyle\\!\leq\\!$
$\displaystyle\frac{(1\\!+\\!\zeta)J^{\mathrm{soft}}_{k-1}(x)\\!+\\!(1\\!-\\!(1\\!+\\!\zeta)(1\\!-\\!\beta))l(x,\\!\pi_{k}(x))\\!-\\!\zeta
l(x,\\!\hat{\pi}^{\mathrm{im}}(x))}{1-\zeta}.$ (54)

With (19) and Lemma 2, (A-E) readily leads to

$\displaystyle
J^{\mathrm{soft}}_{k-1}(x^{+})-J^{\mathrm{soft}}_{k-1}(x)\leq\frac{2\zeta\gamma-(1+\zeta)(1-\beta)}{1-\zeta}l(x,0),$
(55)

which means that $J^{\mathrm{soft}}_{k-1}(\cdot)$ is strictly decreasing from
any
$x\in(\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega))\backslash\\{0\\}$
to the next state. By combining (A-E) and (55), we note that
$\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega)$
is a positively invariant set for the closed-loop system. This together with
the Lyapunov functions $\hat{J}_{k-1}(\cdot)$ and
$J^{\mathrm{soft}}_{k-1}(\cdot)$ leads to the asymptotic stability and safety
of the closed-loop system.

Next, we analyze the behavior of the closed-loop system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$. With C1 and C2,
$\hat{J}_{k-1}(x^{+})-\hat{J}_{k-1}(x)$ with
$x^{+}=f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(x))$ is then
upper bounded by the right-hand side of (A-E) plus $e_{\mathrm{p}}l(x,0)$, and
$J^{\mathrm{soft}}_{k-1}(x^{+})-J^{\mathrm{soft}}_{k-1}(x)$ is also upper
bounded by the right-hand side of (55) plus $e_{\mathrm{p}}l(x,0)$. Together
with (20) this results in the asymptotic stability and safety of the closed-
loop system with policy $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$.

Finally, to prove Corollary 1, we note that the first inequality in (22) and
the first inequality in (23) directly follow from the optimality of
$J^{\mathrm{soft*}}(\cdot)$. We derive from (A-E) that

$\frac{1-2\zeta\gamma}{1-\zeta}l(x,\hat{\pi}^{\mathrm{im}}(x))\leq
J^{\mathrm{soft}}_{k-1}(x)-J^{\mathrm{soft}}_{k-1}(f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{im}}(x)))$
(56)

holds for any
$x\in(\mathcal{B}(J^{\mathrm{soft}}_{k-1},\Omega)\cap\mathcal{B}(\hat{J}_{k-1},\Omega))\backslash\\{0\\}$.
Consider the trajectory $x_{0},x_{1},...$ that is generated by applying
$\hat{\pi}^{\mathrm{im}}(x_{t})$ to system (1) at each time step
$t,\;t=0,1,...$ Letting $x=x_{t}$ and summing up both sides of (56) from $t=0$
to $t=\infty$, we get
$J_{\hat{\pi}^{\mathrm{im}}}(x_{0})\leq\frac{1-\zeta}{1-2\zeta\gamma}(J^{\mathrm{soft}}_{k-1}(x_{0})-J^{\mathrm{soft}}_{k-1}(x_{\infty}))$.
The asymptotic stability in (i) of Theorem 2 indicates that
$J^{\mathrm{soft}}_{k-1}(x_{\infty})=0$, so (i) of Corollary 1 is proved.
Similarly, we can upper bound the stage cost when applying
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ by

$\displaystyle\frac{1-2\zeta\gamma-
e_{\mathrm{p}}}{1-\zeta}l(x,\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(x))$
$\displaystyle\leq$ $\displaystyle
J^{\mathrm{soft}}_{k-1}(x)-J^{\mathrm{soft}}_{k-1}(f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(x))).$

(ii) of Corollary 1 will be obtained by summing up the above inequality along
the trajectory controlled by $\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$.

###  A-F Proof of Theorem 3

If (V-C) returns “Yes”, the trajectories of the closed-loop system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ with initial condition
$x_{0}\in X_{\mathrm{in}}$ will be contained in $\mathcal{B}_{r_{1}}$ after
$N$ time steps.

Then, for any initial state $x_{0}\in\mathcal{B}_{r_{1}}$, suppose that the
$t_{\mathrm{f}}-1$-step trajectory $x_{0},\;x_{1},...,\;x_{t_{\mathrm{f}}-1}$
of the closed-loop system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ is not contained in
$\mathcal{B}_{r_{2}}$. Since $a^{*}_{1}\leq 0$, we have

$\hat{J}_{k-1}(x_{t+1})-\hat{J}_{k-1}(x_{t})\leq-
c_{1}l(x_{t},0),\;t=0,...,t_{\mathrm{f}}-1$ (57)

Summing (57) over time yields

$\hat{J}_{k-1}(x_{t_{\mathrm{f}}})\leq\hat{J}_{k-1}(x_{0})-c_{1}\sum_{t=0}^{t_{\mathrm{f}}-1}l(x_{t},0)$
(58)

Meanwhile, since both $\hat{J}_{k-1}(\cdot)$ and $l(x,0)$ are continuous PWA
functions on their domains, similarly to (i) of Lemma 2, there exists a
positive and finite constant $\hat{\gamma}$ such that
$\hat{J}_{k-1}(x)\leq\hat{\gamma}l(x,0)$ for all $x\in X$. As a result, (58)
implies that $\hat{J}_{k-1}(x_{t_{\mathrm{f}}})\leq
r_{1}-\frac{c_{1}r_{2}t_{\mathrm{f}}}{\hat{\gamma}}$. Specifying
$t_{\mathrm{f}}=\lceil(r_{1}-r_{2})\hat{\gamma}/(c_{1}r_{r})\rceil$, which is
finite and does not depend on $x_{0}$, we have
$\hat{J}_{k-1}(x_{t_{\mathrm{f}}})\leq r_{2}$. Combining the above statements,
we can conclude that any trajectory of the closed-loop system with
$\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(\cdot)$ starting from
$X_{\mathrm{in}}$ will reach $\mathcal{B}_{r_{2}}$ in less than
$N+t_{\mathrm{f}}$ time steps. Finally, the positive invariance of
$\mathcal{B}_{r_{2}}$ is straightforward if $a^{*}_{2}\leq 0$, since we have
$\hat{J}_{k-1}(x)\leq
r_{2}\Rightarrow\hat{J}_{k-1}(f_{\mathrm{PWA}}(x,\hat{\pi}^{\mathrm{ex}}_{\mathrm{proj}}(x)))\leq
r_{2}$ from (V-C). This completes the proof of Theorem 3.

## References

  * [1] T. Liu, Y. Gao, and M. Buss, “Adaptive output tracking control of piecewise affine systems with prescribed performance,” _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ , vol. 52, no. 9, pp. 5398–5410, 2021\.
  * [2] S. Sadraddini and R. Tedrake, “Sampling-based polytopic trees for approximate optimal control of piecewise affine systems,” in _2019 International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 7690–7696.
  * [3] N. Mehr, D. Sadigh, R. Horowitz, S. S. Sastry, and S. A. Seshia, “Stochastic predictive freeway ramp metering from signal temporal logic specifications,” in _2017 American Control Conference (ACC)_ , 2017, pp. 4884–4889.
  * [4] F. Borrelli, M. Baotić, A. Bemporad, and M. Morari, “Dynamic programming for constrained optimal control of discrete-time linear hybrid systems,” _Automatica_ , vol. 41, no. 10, pp. 1709–1721, 2005.
  * [5] D. Görges, “Relations between model predictive control and reinforcement learning,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 4920–4928, 2017.
  * [6] L. Busoniu, R. Babuska, B. De Schutter, and D. Ernst, _Reinforcement Learning and Dynamic Programming Using Function Approximators_. CRC Press, 2017.
  * [7] A. Heydari, “Revisiting approximate dynamic programming and its convergence,” _IEEE Transactions on Cybernetics_ , vol. 44, no. 12, pp. 2733–2743, 2014\.
  * [8] Q. Wei, D. Liu, and H. Lin, “Value iteration adaptive dynamic programming for optimal control of discrete-time nonlinear systems,” _IEEE Transactions on Cybernetics_ , vol. 46, no. 3, pp. 840–853, 2015.
  * [9] R. Schwan, C. N. Jones, and D. Kuhn, “Stability verification of neural network controllers using mixed-integer programming,” _arXiv preprint arXiv:2206.13374_ , 2022.
  * [10] S. Chen, K. Saulnier, N. Atanasov, D. D. Lee, V. Kumar, G. J. Pappas, and M. Morari, “Approximating explicit model predictive control using constrained neural networks,” in _2018 Annual American Control Conference (ACC)_ , 2018, pp. 1520–1527.
  * [11] M. Hertneck, J. Köhler, S. Trimpe, and F. Allgöwer, “Learning an approximate model predictive controller with guarantees,” _IEEE Control Systems Letters_ , vol. 2, no. 3, pp. 543–548, 2018.
  * [12] K. P. Wabersich and M. N. Zeilinger, “A predictive safety filter for learning-based control of constrained nonlinear dynamical systems,” _Automatica_ , vol. 129, p. 109597, 2021.
  * [13] Y. Li, N. Li, H. E. Tseng, A. Girard, D. Filev, and I. Kolmanovsky, “Robust action governor for uncertain piecewise affine systems with non-convex constraints and safe reinforcement learning,” _arXiv preprint arXiv:2207.08240_ , 2022.
  * [14] S. Gros, M. Zanon, and A. Bemporad, “Safe reinforcement learning via projection on a safe set: How to achieve optimality?” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 8076–8081, 2020.
  * [15] A. Chakrabarty, R. Quirynen, C. Danielson, and W. Gao, “Approximate dynamic programming for linear systems with state and input constraints,” in _2019 European Control Conference (ECC)_ , 2019, pp. 524–529.
  * [16] L. Beckenbach, P. Osinenko, T. Göhrt, and S. Streif, “Constrained and stabilizing stacked adaptive dynamic programming and a comparison with model predictive control,” in _2018 European Control Conference (ECC)_ , 2018, pp. 1349–1354.
  * [17] J. Duan, Z. Liu, S. E. Li, Q. Sun, Z. Jia, and B. Cheng, “Adaptive dynamic programming for nonaffine nonlinear optimal control problem with state constraints,” _Neurocomputing_ , vol. 484, pp. 128–141, 2022.
  * [18] C. Tessler, D. J. Mankowitz, and S. Mannor, “Reward constrained policy optimization,” _arXiv preprint arXiv:1805.11074_ , 2018.
  * [19] J. Li, D. Fridovich-Keil, S. Sojoudi, and C. J. Tomlin, “Augmented Lagrangian method for instantaneously constrained reinforcement learning problems,” in _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 2982–2989.
  * [20] S. Paternain, M. Calvo-Fullana, L. F. Chamon, and A. Ribeiro, “Learning safe policies via primal-dual methods,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 2019, pp. 6491–6497.
  * [21] D. Yu, H. Ma, S. Li, and J. Chen, “Reachability constrained reinforcement learning,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 25 636–25 655.
  * [22] A. Heydari, “Stability analysis of optimal adaptive control using value iteration with approximation errors,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 9, pp. 3119–3126, 2018.
  * [23] R. Postoyan, M. Granzotto, L. Buşoniu, B. Scherrer, D. Nešsić, and J. Daafouz, “Stability guarantees for nonlinear discrete-time systems controlled by approximate value iteration,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 2019, pp. 487–492.
  * [24] F. Moreno-Mora, L. Beckenbach, and S. Streif, “Predictive control with learning-based terminal costs using approximate value iteration,” _arXiv preprint arXiv:2212.00361_ , 2022.
  * [25] S. Gros and M. Zanon, “Learning for MPC with stability & safety guarantees,” _Automatica_ , vol. 146, p. 110598, 2022.
  * [26] X. Yang, H. Zhang, Z. Wang, H. Yan, and C. Zhang, “Data-based predictive control via multistep policy gradient reinforcement learning,” _IEEE Transactions on Cybernetics_ , vol. 53, no. 5, pp. 2818–2828, 2021.
  * [27] K. He, T. v. d. Boom, and B. De Schutter, “Approximate dynamic programming for constrained linear systems: A piecewise quadratic approximation approach,” _arXiv preprint arXiv:2205.10065_ , 2022.
  * [28] L. Buşoniu, T. de Bruin, D. Tolić, J. Kober, and I. Palunko, “Reinforcement learning for control: Performance, stability, and deep approximators,” _Annual Reviews in Control_ , vol. 46, pp. 8–28, 2018.
  * [29] M. Baoti, F. J. Christophersen, and M. Morari, “Constrained optimal control of hybrid systems with a linear performance index,” _IEEE Transactions on Automatic Control_ , vol. 51, no. 12, pp. 1903–1919, 2006.
  * [30] F. Borrelli, A. Bemporad, and M. Morari, _Predictive Control for Linear and Hybrid Systems_. Cambridge University Press, 2017.
  * [31] S. Chen, M. Fazlyab, M. Morari, G. J. Pappas, and V. M. Preciado, “Learning Lyapunov functions for hybrid systems,” in _Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control_ , 2021, pp. 1–11.
  * [32] H. Dai, B. Landry, L. Yang, M. Pavone, and R. Tedrake, “Lyapunov-stable neural-network control,” _arXiv preprint arXiv:2109.14152_ , 2021.
  * [33] B. Karg and S. Lucia, “Stability and feasibility of neural network-based controllers via output range analysis,” in _2020 59th IEEE Conference on Decision and Control (CDC)_ , 2020, pp. 4947–4954.
  * [34] J. Xu, T. van den Boom, L. Buşoniu, and B. De Schutter, “Model predictive control for continuous piecewise affine systems using optimistic optimization,” in _2016 American Control Conference (ACC)_ , 2016, pp. 4482–4487.
  * [35] M. Lazar, W. Heemels, S. Weiland, and A. Bemporad, “Stabilizing model predictive control of hybrid systems,” _IEEE Transactions on Automatic Control_ , vol. 51, no. 11, pp. 1813–1818, 2006.
  * [36] A. Thirugnanam, J. Zeng, and K. Sreenath, “Safety-critical control and planning for obstacle avoidance between polytopes with control barrier functions,” in _2022 International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 286–292.
  * [37] D. P. Bertsekas, _Reinforcement Learning and Optimal Control_. Athena Scientific, 2019.
  * [38] ——, “Value and policy iterations in optimal control and adaptive dynamic programming,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 28, no. 3, pp. 500–509, 2015.
  * [39] F. Borrelli, M. Baotic, A. Bemporad, and M. Morari, “An efficient algorithm for computing the state feedback optimal control law for discrete time hybrid systems,” in _2003 American Control Conference_ , vol. 6, 2003, pp. 4717–4722.
  * [40] J. Xu, J. Wang, J. Rao, Y. Zhong, and H. Wang, “Adaptive dynamic programming for optimal control of discrete-time nonlinear system with state constraints based on control barrier function,” _International Journal of Robust and Nonlinear Control_ , vol. 32, no. 6, pp. 3408–3424, 2022.
  * [41] D. P. Bertsekas, “Nonlinear programming,” _Journal of the Operational Research Society_ , vol. 48, no. 3, pp. 334–334, 1997.
  * [42] B. Hanin, “Universal function approximation by deep neural nets with bounded width and ReLU activations,” _Mathematics_ , vol. 7, no. 10, p. 992, 2019\.
  * [43] W. Lohmiller, P. Gassert, and J.-J. Slotine, “Deep minmax networks,” in _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 2929–2934.
  * [44] L. A. Wolsey and G. L. Nemhauser, _Integer and Combinatorial Optimization_. John Wiley & Sons, 1999\.
  * [45] D. Liu, Q. Wei, D. Wang, X. Yang, and H. Li, “Finite approximation error-based value iteration ADP,” _Adaptive Dynamic Programming with Applications in Optimal Control_ , pp. 91–149, 2017.
  * [46] M. Herceg, M. Kvasnica, C. Jones, and M. Morari, “Multi-Parametric Toolbox 3.0,” in _2013 European Control Conference (ECC)_ , Zürich, Switzerland, July 17–19 2013, pp. 502–510, https://www.mpt3.org/.
  * [47] D. Corona and B. De Schutter, “Adaptive cruise control for a SMART car: A comparison benchmark for MPC-PWA control methods,” _IEEE Transactions on Control Systems Technology_ , vol. 16, no. 2, pp. 365–372, 2008.
  * [48] B. Bank, J. Guddat, D. Klatte, B. Kummer, and K. Tammer, _Non-linear Parametric Optimization_. Springer, 1983\.
  * [49] G. F. Montufar, R. Pascanu, K. Cho, and Y. Bengio, “On the number of linear regions of deep neural networks,” _Advances in Neural Information Processing Systems_ , vol. 27, 2014.
  * [50] A. Kripfganz and R. Schulze, “Piecewise affine functions as a difference of two convex functions,” _Optimization_ , vol. 18, no. 1, pp. 23–29, 1987\.
  * [51] V. Tjeng, K. Xiao, and R. Tedrake, “Evaluating robustness of neural networks with mixed integer programming,” _arXiv preprint arXiv:1711.07356_ , 2017\.
  * [52] L. Weng, H. Zhang, H. Chen, Z. Song, C.-J. Hsieh, L. Daniel, D. Boning, and I. Dhillon, “Towards fast computation of certified robustness for ReLU networks,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 5276–5285.
  * [53] A. Bemporad, “A piecewise linear regression and classification algorithm with application to learning and model predictive control of hybrid systems,” _IEEE Transactions on Automatic Control_ , 2022.
  * [54] A. Bemporad and M. Morari, “Control of systems integrating logic, dynamics, and constraints,” _Automatica_ , vol. 35, no. 3, pp. 407–427, 1999.