# Computing control invariant sets of nonlinear systems: decomposition and
distributed computing

Benjamin Decardi-Nelsona, Jinfeng Liua,

a Department of Chemical & Materials Engineering, University of Alberta,

Edmonton AB T6G 1H9, Canada
Corresponding author: J. Liu. Tel: +1-780-492-1317. Fax: +1-780-492-2881.
Email<EMAIL_ADDRESS>
###### Abstract

In this work, we present a distributed framework based on the graph algorithm
for computing control invariant set for nonlinear cascade systems. The
proposed algorithm exploits the structure of the interconnections within a
process network. First, the overall system is decomposed into several
subsystems with overlapping states. Second, the control invariant set for the
subsystems are computed in a distributed manner. Finally, an approximation of
the control invariant set for the overall system is reconstructed from the
subsystem solutions and validated. We demonstrate the efficacy and convergence
of the proposed method to the centralized graph-based algorithm using several
numerical examples including a six dimensional continuous stirred tank reactor
system.

Keywords: Nonlinear systems; Outer approximation; Invariant sets; Graph
theory; Distributed computation.

##  1 Introduction

Many engineered systems such as chemical processes are usually designed to
operate within a prescribed safe set – usually in the form of state
constraints. However, the ability of the state of these engineered systems to
stay within the prescribed safe set is hindered by limits on the available
control energy used to control these systems [1]. Thus, to avoid unsafe and
catastrophic events during operation, it is very critical that the set of
states for which a control input exist to keep the system’s states within the
state constraint at all times is known prior to operation or for controller
design purposes. This set of states is known as a control invariant set. The
concept of control invariant set has close connections with reachable sets
[2], null controllable sets [3] and viability kernels [4].

Control invariant sets have been found to be important in controller design
and assessment [5]. In particular, it has been found to be useful as a
terminal region constraint for model predictive control strategies [6, 7].
Despite the well known importance of control invariant sets, accurate and
efficient determination of control invariant sets is still an open issue, even
for linear systems. This is because determining a control invariant set for a
given system involves considering all possible trajectories of the system.
Unfortunately, computing all possible system trajectories become challenging
when the system dimension is high, a phenomenon known as curse of
dimensionality. For this reason, many current state-of-the-art algorithms have
inadequate scalability when it comes to the dimension of the system states.

To this end, several algorithms and methods for computing control invariant
sets have been proposed [8, 9, 10, 11, 12]. Linear systems which are usually
based on the well-known dynamic programming method [13] and polytopic
computations, have in particular received considerable attention in literature
[9, 10]. However, these methods for linear time invariant (LTI) systems are
based on vertex enumeration of the set’s full-dimensional polytopic
representation [10]. The number of vertices in the polytope grows
exponentially as the system’s dimension increase, resulting in intractable
computation. For nonlinear systems, only a few results exist with majority of
the algorithms either being grid-based [1, 2, 14] or set-based [8]. Similar
issues which are considerably more challenging exist for nonlinear systems.

A trade-off between set complexity and computational complexity is frequently
made to overcome scalability concerns in determining control invariant sets.
The methods presented in [15, 16] have been successful in computing control
invariant set of high dimensional systems but assume a control invariant set
of prescribed complexity or shape which results in faster computation but
overly conservative sets. A decomposition and reconstruction based method for
exact computation of backward reachable sets of self-contained nonlinear
subsystems using the Halmilton-Jacobi (HJ) formulations was developed in [17].
In [18, 9], decentralized and distributed algorithms were developed to obtain
a family of practical positively invariant sets of linear systems. In [19], a
decomposition based method was proposed to compute the backward reachable sets
of nonlinear systems. In these cases, however, the missing interconnection
information were treated as disturbances which results in conservative
results.

In this paper, we focus on cascade systems and present a system decomposition
method and a distributed approach for computing control invariant sets. The
proposed algorithm exploits the structure of the interconnections within a
process network and decomposes the entire process network into smaller
subsystems. Following the decomposition, a distributed approach is developed
to compute the control invariant set of the entire system. The proposed
approach adopts graph-based algorithms in computing the control invariant sets
[8]. In contrast to other works on distributed computation of control
invariant sets, our proposed approach produces sets that approximates the
largest control invariant set since the missing interconnection information is
not treated as disturbances. The proposed distributed computing algorithm may
be extended to more general nonlinear systems. Some preliminary results of
this work were reported in [20]. Compared with [20], this paper provides
significantly extended explanations and discussion, and significantly more
results exploring various features of the proposed distributed computing
approach.

The remainder of the paper is organized as follows: Section 2 is concerned
with preliminaries and problem formulation as well as a brief introduction to
directed graph representation of dynamical system with controls and control
invariant set investigation of the graph. In Section 3, we present the system
decomposition method with overlapping states and discuss the implications of
computing control invariant set via system decomposition. In Section 4, we
present an algorithm to compute the outer approximation of the largest control
invariant set in a distributed manner. We also briefly discuss some
computational issues in this section. We present a several examples in Section
5 to demonstrate the efficacy of the proposed method and then conclude the
paper in Section 6.

##  2 Problem formulation and background

###  2.1 Notation

Throughout this article, the operator proj${}_{i}(x)$ denotes the projection
of the set or point $x$ onto the subspace of subsystem $i$. The operator $/$
denotes set subtraction such that
$\mathbb{A}/\mathbb{B}=\\{x:x\in\mathbb{A},x\notin\mathbb{B}\\}$. $G=(V,E)$
represents a directed graph with $V$ denoting the set of vertices of the graph
and $E$ denoting the set of ordered pairs of vertices known as edges. The
operator $|\cdot|$ denotes the Euclidean norm of a vector.

###  2.2 Problem formulation

We are concerned with a class of nonlinear systems composed of $N$ subsystems
coupled together in a cascade manner.

Figure 1:  A schematic diagram of a cascade process system

A schematic diagram of the cascade process system is shown in Figure 1. The
dynamics of the overall system is described by

$x^{+}=f(x,u)$ (1)

where $x\in\mathbb{R}^{n}$ is the current state of the system,
$u\in\mathbb{R}^{m}$ is the current control input, and
$x^{+}\in\mathbb{R}^{n}$ denotes the state of the system at the next sampling
time. We assume that the state and the control input of the system are
restricted to be in the compact constraint sets $X\subset\mathbb{R}^{n}$ and
$U\subset\mathbb{R}^{m}$ respectively. Without loss of generality, we also
assume that the vector field $f:X\times U\rightarrow X$ is a sufficiently
smooth vector field in $X$. The coupling between the subsystems is such that
the first subsystem is independent of the other subsystems. Also, any
subsystem $i$, other than the first subsystem, is directly affected by the
upstream subsystem $i-1$ but not the downstream subsystem $i+1$. The structure
of the cascade system considered in this article is represented by

$\displaystyle x_{1}^{+}$ $\displaystyle=f_{1}(x_{1},u_{1})$ (2a)
$\displaystyle x_{2}^{+}$ $\displaystyle=f_{2}(x_{2},u_{2})+g_{2}(x_{1})$ (2b)
$\displaystyle\vdots$ $\displaystyle x_{N}^{+}$
$\displaystyle=f_{N}(x_{N},u_{N})+g_{N}(x_{N-1})$ (2c)

where $f_{i}$ denotes the local dynamics of state $x_{i}$, and $g_{i}$ denotes
the dynamics of the coupling between the states $x_{i}$ and $x_{i-1}$. The
subscript $i=1,2,\cdots,N$ represents the $i$th subsystem. The state and input
constraints for subsystem $i$ is given by $X_{i}=\text{proj}_{i}(X)$ and
$U_{i}=\text{proj}_{i}(U)$ respectively.

Before we begin our discussion, let us introduce the following definitions
which are central and referred to throughout this article.

######  Definition 1 (Forward invariant set [5]).

A set $R\subseteq X$ is said to be a forward or positively invariant set of
the system $x^{+}=f(x)$ if for every $x\in R$, $f(x)\in R$.

######  Definition 2 (Control invariant set [5]).

A set $R\subseteq X$ is said to be a control invariant set (CIS) of system (1)
if for every $x\in R$, there exist a feedback control law $u=\mu(x)\in U$ such
that $R$ is forward invariant for the closed-loop system $f(x,u)$.

######  Definition 3 (Largest control invariant set [21]).

A set $R_{X}\subseteq X$ is said to be the largest (with respect to inclusion)
control invariant set of system (1) if $R_{X}$ is control invariant and
contains all other control invariant sets contained in $X$.

In general, the state constraint $X$ for a given system is not control
invariant. However, one may wish to find the largest control invariant set
$R_{X}$ contained in $X$ for controller design and assessment purposes. The
goal of this article is to present a framework for computing an approximation
of the largest control invariant set $R_{X}$ when the state dimension $n$ is
too large to make the standard GIS algorithm tractable [8]. At present, this
happens when $n>4$. This is because the number of cells generated in the
standard GIS algorithm grows exponentially in the state dimension.

Our solution to computing an approximation of $R_{X}$ is to decompose system
(1) into $M$ lower dimensional subsystems with overlapping states.

###  2.3 Graph-based invariant set algorithm

In this section, we briefly introduce the concept of graph-based invariant set
(GIS) algorithm. The reader may refer to the works of Decardi-Nelson and Liu
[8] and Osipenko [22] for a more detailed discussion.

Given a dynamical system in the form of Equation (1), the idea is to cast the
goal of finding $R_{X}$ as a graph theoretical problem. To achieve this, we
first rewrite the function $f$ as set a set-valued map parameterized by the
input

$F(x):=f(x,U)=\\{f(x,u)\\}_{\cup_{u\in U}}$ (3)

This makes system (1) a difference inclusion of the form

$x^{+}\in F(x)$ (4)

The function $F$ associates the state $x$ with the subset $F(x)$ of all
feasible next states. To construct a digraph representation (also known as
symbolic image) of system (1), the state constraint $X$ first need to be
quantized into a collection of cells
$\mathcal{C}=\\{B_{1},B_{2},\cdots,B_{l}\\}$ such that

$R=\cup_{B_{i}\in\mathcal{C}}$ (5)

is an outer covering of $R_{X}$. The diameter $d$ of the collection is given
as

$d=\text{diam}(\mathcal{C}):=\max_{B_{i}\in\mathcal{C}}\text{diam}(B_{i})$

where diam($B_{i}$) = $\text{sup}\\{|x-y|:x,y\in B_{i}\\}$.

Following the quantization of the state constraint, the dynamics of system (1)
is approximated using a directed graph $G=(V,E)$ such that

$\displaystyle V$ $\displaystyle=\mathcal{C}~{}and$ (6) $\displaystyle E$
$\displaystyle=\\{(B_{i},B_{j})\in\mathcal{C}\times\mathcal{C}~{}|~{}F(B_{i})\cap
B_{j}\neq\emptyset\\}$ (7)

where $V$ is the vertices of the directed graph and $E$ denotes the set of
ordered pair of vertices known as edges. From Equations (6) and (7), it is
easy to see that the finer the quantization (the smaller the diameter of the
collection), the more accurate the approximation of the dynamics of system (1)
using a directed graph $G$.

The cells that outer approximate $R_{X}$ may be obtained by analyzing the
resulting directed graph $G$. Specifically, the goal is to split the vertices
of $G$ into two groups namely, leaving and non-leaving cells. The leaving
cells are cells that have only finite paths passing through them whereas the
nonleaving cells are the cells that have infinite admissible paths passing
through them. Let us recall the following definitions in graph theory.

######  Definition 4 (Strongly connected graph).

A directed graph $G=(V,E)$ is said to be strongly connected if there is an
admissible path in both directions between each pair of vertices of the graph.

######  Definition 5 (Strongly connected component (SCC)).

A strongly connected component of a directed graph $G=(V,E)$ is a maximal
strongly connected subgraph of $G$.

Algorithmically, finding the non-leaving cells involves finding the SCC of $G$
and any vertex that is not in the SCC but has a path to any vertex in the SCC.
A graphical illustration of the shown in Figure 2.

Figure 2:  A graphical depiction of the GIS algorithm. (a) Quantization of the
state constraint. (b) Graph construction and analysis. The green cells
constitute the SCC; the yellow cells have an admissible path to any vertex of
the SCC; the red cells do not belong to the SCC nor have a path to the SCC (c)
The final selected cells with infinite admissible paths passing through them.
The remaining or leaving cells are discarded.

The following proposition summarizes how an outer approximation of $R_{X}$ may
be obtained using graph-based tools.

######  Proposition 1 (Non-leaving cells [8]).

Let G=(V,E) be a directed graph approximation of the dynamics of system (1)
with respect to the collection $\mathcal{C}$. Then

  1. i.

The vertices $V_{s}$ of the SCC of G has infinite admissible paths passing
through them.

  2. ii.

any element of $V$ but not in $V_{s}$, that is $V/V_{s}$, with a path to at
least one path to any element of the SCC also has an infinite admissible path
passing through it.

  3. iii.

the union of the elements in (i) and (ii), denoted as $I^{+}(G)$, is a closed
neighborhood of the largest control invariant set $R_{X}$ contained in $X$,
that is

$R_{X}\subseteq I^{+}(G)$

Furthermore, it can be shown that $R$ converges to the largest control
invariant set $R_{X}$ as $d\rightarrow 0$ [8].

##  3 System decomposition and set invariance

In trying to alleviate the exponential cell growth in the GIS algorithm, we
propose to decompose the overall system into smaller subsystems. This makes it
computationally tractable for control invariant set approximation using the
GIS algorithm. System decomposition is a typical way of addressing the
computational challenges associated with control and state estimation of large
scale dynamical systems. This is also the case for control invariant set
calculation. The goal is to divide the original system into many small
subsystems thus making the computations tractable. Usually, the decomposition
is achieved by dividing the system in such a way that the subsystems have weak
to no coupling or interconnection.

In this section, we investigate the system structures suitable for
decomposition and control invariant set computation. In particular, we
consider simple parallel and series/cascade system structures. Throughout this
section, we restrict the discussion to two and three dimensional system
structures namely, series and parallel system structures to make the ideas
presented here easier to follow. The ideas presented here can easily be
generalized to much higher dimensional systems.

We begin this section by analyzing two types of system structures, their
decomposition and the ability to reconstruct the solution from the solution of
the subsystems. Thereafter we described the method of overlapping
decomposition, which is a precursor for the proposed distributed algorithm.

###  3.1 System structures and invariance

Consider the disjoint system structure

$x^{+}=\hat{f}(x)=\begin{bmatrix}\hat{f}_{1}(x_{1})\\\
\hat{f}_{2}(x_{2})\end{bmatrix}=\begin{bmatrix}A_{11}&0\\\
0&A_{22}\end{bmatrix}x$ (8)

where $x=[x_{1}~{}x_{2}]^{T}$ subject to the state constraint $X$. A graphical
depiction of the system is shown in Figure 3.

Figure 3:  Parallel or independent system structure

Since the evolution of the system states are independent of each other, the
system equations can be intuitively decomposed into the following subsystems

$S_{1}:x_{1}^{+}=\hat{f}_{1}(x_{1})=A_{11}x_{1}\text{ and
}S_{2}:x_{2}^{+}=\hat{f}_{2}(x_{2})=A_{22}x_{2}$ (9)

with subspaces $X_{1}$ = proj1($X$) and $X_{2}$ = proj2($X$) respectively. To
construct a graph representation of the dynamics of each subsystem, the
subspaces $X_{1}$ and $X_{2}$ need to be quantized into $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$. In what follows, we show how the control invariant set can
be computed for the full system from the subsystem information.

######  Definition 6 (Graph Cartesian product [23]).

The Cartesian product of the directed graphs $G_{1}=(V_{1},E_{1})$ and
$G_{2}=(V_{2},E_{2})$, denoted $G_{1}\times G_{2}$, is a graph $G=(V,E)$ such
that $V=V_{1}\times V_{2}$ and for any two points $u=(u_{1},u_{2})$ and
$v=(v_{1},v_{2})$ in $V$, the directed edge $(u,v)\in E$ whenever
$u_{1}=v_{1}$ and $(u_{2},v_{2})\in E_{2}$, or $u_{2}=v_{2}$ and
$(u_{1},v_{1})\in E_{1}$.

An illustration of the Cartesian product of two fictitious digraphs $G_{1}$
and $G_{2}$ is presented in Figure 4.

Figure 4:  The Cartesian product of two graphs

######  Assumption 1 (Subsystem quantization).

The set $\mathcal{C}$, which is a quantization of $X$, is the cross product of
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$, that is
$\mathcal{C}=\mathcal{C}_{1}\times\mathcal{C}_{2}$.

######  Proposition 2.

Consider the system described by Equation (8) and decomposed in the form of
Equation (9). Let $G_{1}=(V_{1},E_{1})$ and $G_{2}=(V_{2},E_{2})$ be a
directed graph representation of the dynamics of each subsystem $S_{1}$ and
$S_{2}$ respectively based on the sets $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$. Also, let $G=(V,E)$ be the directed graph representation of
system (8) based on $\mathcal{C}$. If Assumption 1 holds, then the full system
solution can be exactly obtained from the subsystem solutions, that is,

$I^{+}(G)=I^{+}(G_{1}\times G_{2})$

###### Proof.

To prove the above assertion, we first need to show that the vertices of the
two graphs $G$ and $G_{1}\times G_{2}$ are equal. Thereafter, we need to show
that if there is an admissible path between any two cells in $G$, then there
is an equivalent admissible path between those same two cells in $G_{1}\times
G_{2}$.

First, we show that the vertices of the two graphs $G$ and $G_{1}\times G_{2}$
are equal. It follows directly from Assumption 1 that
$\mathcal{C}=\mathcal{C}_{1}\times\mathcal{C}_{2}$. Therefore from Definition
6, we have that $V=V_{1}\times V_{2}$.

Now, we show by that if there is a path between any two cells in $G$, then
there exist a path between those same two cells in $G_{1}\times G_{2}$. Let
$B_{i}=(u_{1},u_{2})$ and $B_{j}=(v_{1},v_{2})$ be any two cells in $V$ such
that $u_{1}$ and $v_{1}$ are in $V_{1}$, and $u_{2}$ and $v_{2}$ are cells in
$V_{2}$. From the graph construction procedure in Equation (7), we have that

$\displaystyle E$
$\displaystyle=\\{(B_{i},B_{j})\in\mathcal{C}\times\mathcal{C}:\hat{f}(B_{i})\cap
B_{j}\neq\emptyset\\}$
$\displaystyle=\\{(B_{i},B_{j})\in\mathcal{C}\times\mathcal{C}:\hat{f}_{1}(u_{1})\cap
v_{1}\neq\emptyset\land\hat{f}_{2}(u_{2})\cap v_{2}\neq\emptyset\\}$
$\displaystyle=\\{(B_{i},B_{j})\in\mathcal{C}\times\mathcal{C}:(u_{1},v_{1})\in
E_{1}\land(u_{2},v_{2})\in E_{2}\\}$

From the definition of $G_{1}\times G_{2}$ and $E$, we have that the path
$(u_{1},u_{2})\rightarrow(u_{1},v_{2})\rightarrow(v_{1},v_{2})$ and
$(u_{1},u_{2})\rightarrow(v_{1},u_{2})\rightarrow(v_{1},v_{2})$ exist in
$G_{1}\times G_{2}$. This implies that there is an admissible path from
$(u_{1},u_{2})$ to $(v_{1},v_{2})$ in $G_{1}\times G_{2}$. It then follows
that if an infinite admissible path passes through the cell $B_{i}$ in $G$,
then an infinite admissible path also passes through the cell
$B_{i}=(u_{1},u_{2})$ in $G_{1}\times G_{2}$. Hence we have that

$I^{+}(G)=I^{+}(G_{1}\times G_{2}),$

which completes the proof. ∎

While we have shown that an approximation of $R_{X}$ can be obtained for the
disjoint system by constructing the digraph $G_{1}\times G_{2}$, in most cases
this will require a considerable amount of memory to store the graph for
higher dimensional systems. The following theorem shows how an approximation
of the largest control invariant set can be exactly obtained from the
solutions of the subsystems.

######  Theorem 1.

Consider the system described by Equation (8) and decomposed in the form of
Equation (9). Let $G_{1}=(V_{1},E_{1})$ and $G_{2}=(V_{2},E_{2})$ be a
directed graph representation of the dynamics of each subsystem $S_{1}$ and
$S_{2}$ respectively based on the sets $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$. Also, let $G=(V,E)$ be the directed graph representation of
system (8) based on $\mathcal{C}$. If Assumption 1 holds, then the full system
solution can be exactly obtained from the subsystem solutions, that is,

$I^{+}(G)=I^{+}(G_{1})\times I^{+}(G_{2})$

###### Proof.

We know from Proposition 2 that $I^{+}(G)=I^{+}(G_{1}\times G_{2})$. Let
$B_{i}=(u_{1},u_{2})\in V$ where $u_{1}\in V_{1}$ and $u_{2}\in V_{2}$. It can
be inferred from the definition of $G_{1}\times G_{2}$ that wherever there is
an infinite admissible path passing through $B_{i}$, there is also infinite
admissible paths passing through $u_{1}$ in $G_{1}$ and $u_{2}$ in $G_{2}$.
This implies that

$I^{+}(G)=I^{+}(G_{1})\times I^{+}(G_{2})\qed$

Theorem 1 shows that the graph representation of the subsystems can be
analyzed independently. Then the cells that approximate the largest control
invariant set for full system can be computed by finding the Cartesian product
of the individual subsystem solutions. This is computationally more efficient
since a large graph is not constructed for analysis. To illustrate the ideas
presented so far, let us consider Example 1.

######  Example 1.

Consider the system

$x^{+}=\begin{bmatrix}2&0\\\ 0&2\end{bmatrix}x+\begin{bmatrix}1&0\\\
0&1\end{bmatrix}u$ (10)

where the constraints on the states and input are
$X=\\{x\in\mathbb{R}^{2}:\|x\|_{\infty}\leq 5\\}$ and
$U=\\{u\in\mathbb{R}^{2}:\|u\|_{\infty}\leq 1\\}$ respectively.

Figure 5:  Decomposition and full solution reconstruction for system (10).
Left: Solution of Subsystem 1. Middle: Solution of Subsystem 2. Right: The
reconstructed full system solution from the subsystem solutions (black line)
and the actual solution of the full system (light blue region).

The system in Example 1 can be decomposed into two subsystems in a form
similar to Equation (9). Figure 5 shows the largest control invariant sets
$R_{X_{1}}$ and $R_{X_{2}}$ for each of the subsystems, the largest control
invariant set $R_{X}$ for the full system and the reconstructed solution from
the control invariant sets of the subsystems. It can be seen that the full
system solution can be exactly obtained from the subsystem solution without
incurring any error due to the decomposition. This agrees well with our
earlier arguments in Theorem 1, that $R_{X}=R_{X_{1}}\times R_{X_{2}}$.
Ultimately, this shows that we can freely move between the control invariant
sets of the lower dimension subsystems and the control invariant set for the
full system without incurring an losses. Geometrically, the structure of the
control invariant set is always going to be a box.

While an approximation of the largest control invariant set can be obtained
from the solutions of the subsystems without incurring any errors due to
decomposition for disjoint subsystems, such systems rarely occur in the real
world. In most cases, the subsystems may be interconnected. Let us consider
one of such cases where is coupling between the system states. Consider the
system

$x^{+}=\hat{f}(x)=\begin{bmatrix}A_{11}&0\\\ A_{21}&A_{22}\end{bmatrix}x$ (11)

where $x=[x_{1}~{}x_{2}]^{T}$ subject to the state constraint $X$. A graphical
depiction of the system is shown in Figure 6. Similar to the disjoint
structure in Equation (8), system (11) can be decomposed into two subsystems
such that

$S_{1}:x_{1}^{+}=\hat{f}_{1}(x_{1})=A_{11}x_{1}\text{ and
}S_{2}:x_{2}^{+}=\hat{f}_{2}(x_{2})+\hat{g}_{2}(\tilde{x}_{1})=A_{22}x_{2}+A_{21}\tilde{x}_{1}$
(12)

with subspaces $X_{1}$ = proj1($X$) and $X_{2}$ = proj2($X$) respectively.
Notice that while Subsystem 1 is independent of the dynamics of Subsystem 2,
Subsystem 2 is dependent on the dynamics of Subsystem 1. The state information
about Subsystem 1 required by Subsystem 2 is denoted by $\tilde{x}_{1}$ in
Equation (12). To work with such a decomposition, the fundamental issues that
needs to be addressed in this scenario are how to obtain the values of
$\tilde{x}_{1}$ and how to use this information in Subsystem 2.

Figure 6:  Series or connected system structures

It is easy to see that the range of values for $\tilde{x}_{1}$ can be obtained
from the solution of Subsystem 1. For the choice of usage of the
$\tilde{x}_{1}$ information, one approach is to use the $\tilde{x}_{1}$
information as an input in Subsystem 2. However, this more often than not
result in an over approximation of the solution of the full system.

######  Theorem 2.

Consider the system described by Equation (11) and decomposed in the form of
Equation (12). Let $G_{1}=(V_{1},E_{1})$ and $G_{2}=(V_{2},E_{2})$ be a
directed graph representation of the dynamics of each subsystem $S_{1}$ and
$S_{2}$ respectively based on the sets $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$. Also, let $G=(V,E)$ be the directed graph representation of
system (11) based on $\mathcal{C}$. If Assumption 1 holds, then

$I^{+}(G)\subseteq I^{+}(G_{1})\times I^{+}(G_{2})$

###### Proof.

This proof follows along the same lines as Theorem 1. We know from Proposition
2 that $I^{+}(G)=I^{+}(G_{1}\times G_{2})$. Let $B_{i}=(u_{1},u_{2})\in V$
where $u_{1}\in V_{1}$ and $u_{2}\in V_{2}$. It can be inferred from the
definition of $G_{1}\times G_{2}$ that wherever there is an infinite
admissible path passing through $B_{i}$, there is also infinite admissible
paths passing through $u_{1}$ in $G_{1}$ and $u_{2}$ in $G_{2}$. This implies
that

$I^{+}(G)\subseteq I^{+}(G_{1})\times I^{+}(G_{2})\qed$

The reverse of Theorem 2, that is $I^{+}(G_{1})\times I^{+}(G_{2})\subseteq
I^{+}(G)$ is necessarily not true. To demonstrate this, let us consider the
following example.

######  Example 2.

Consider the system

$x^{+}=\begin{bmatrix}2&0\\\ 1&2\end{bmatrix}x+\begin{bmatrix}1&0\\\
0&1\end{bmatrix}u$ (13)

where the constraints on the states and input are
$X=\\{x\in\mathbb{R}^{2}:\|x\|_{\infty}\leq 5\\}$ and
$U=\\{u\in\mathbb{R}^{2}:\|u\|_{\infty}\leq 1\\}$ respectively.

As an illustration of our arguments so far, consider the system in Example 2.
System (13) can be decomposed into two subsystems following Equation (12).
Figure 7 shows the solution of the subsystems $R_{X_{1}}$ and $R_{X_{2}}$,
when $\tilde{x}_{1}$ is treated as an input in Subsystem 2, the reconstructed
solution of the full system from the solution of the subsystems, and the
solution of the full system $R_{X}$. It can be seen that $R_{X}\subseteq
R_{X_{1}}\times R_{X_{2}}$. Hence, there exists some states in
$R_{X_{1}}\times R_{X_{2}}$ which are not in $R_{X}$.

Figure 7:  Decomposition and full solution reconstruction for system (13).
Left: Solution of Subsystem 1. Middle: Solution of Subsystem 2 when
$\tilde{x}_{1}$ is treated as an input. Right: The reconstructed full system
solution from the subsystem solutions (black line) and the actual solution of
the full system (light blue region).

To address the issue of over approximation, several studies have proposed to
treat $\tilde{x}_{1}$ as an uncertain input in Subsystem 2 [19, 24]. Let us
introduce the following definition

######  Definition 7 (Robust control invariant set).

A set $R$ is said to be a robust control invariant set of the system
$x^{+}=f(x,u,w)$, if for every $x\in R$, there exist a feedback control law
$u=\mu(x)\in U$ such that $R$ is forward invariant for the closed-loop system
$f(x,u,w)$ for all uncertainty $w$ in the uncertainty set $W$.

By treating $\tilde{x}_{1}$ as an uncertain input, a worse case scenario
analysis can be considered during the computation of the solution of Subsystem
2. This results in a robust control invariant set solution for Subsystem 2.
While this approach works, there are several downsides to it. First, it is
neither trivial nor easy to compute robust control invariant sets. In fact,
more resources are needed to compute robust control invariant sets than to
compute control invariant sets. Second, the resulting reconstructed set is
often conservative or small compared to the control invariant set for the full
system. In some instances, the computation may end up with an empty set for
Subsystem 2 which is undesired. Let us consider the system in Example 2 again.
When $\tilde{x}_{1}$ is treated as an uncertainty in its worst case,
$R_{X_{2}}=\emptyset$, implying that $R_{X_{2}}\times R_{X_{2}}=\emptyset$,
which is not desirable. This also shows that $R_{X_{1}}\times
R_{X_{2}}\subseteq R_{X}$ when the $\tilde{x}_{1}$ is treated as an
uncertainty, which is a conservative result.

In a nutshell, it can be seen that when the subsystems are connected, the full
system solution cannot be easily obtained from the solution of the subsystems.
Further scrutinizing the actual control invariant set for system (13) in
Figure 7, it can be seen that $\tilde{x}_{1}$ actually behaves as an input at
some points in $x_{2}$ and as an uncertainty in other $x_{2}$ locations. This
is the central idea we upon which we develop our distributed algorithm.

###  3.2 Overlapping system decomposition

In the preceding section, we have discussed the implication of decomposition
of set invariance using 2 dimensional systems. But what how do we decompose
connected systems when the system dimension is greater than 2? In this
section, we briefly discuss how to decompose three dimensional cascade
systems.

Let us consider the following system

$x^{+}=\hat{f}(x)=\begin{bmatrix}A_{11}&0&0\\\ A_{21}&A_{22}&0\\\
0&A_{32}&A_{33}\\\ \end{bmatrix}x$ (14)

where $x=[x_{1}~{}x_{2}~{}x_{3}]^{T}$ subject to the state constraint $X$.

Figure 8:  Different decomposition strategies for the cascade system

There are different ways to decompose system (14) into several subsystems. A
simple but naive decomposition is to consider each state in Equation (14)
separately as shown in Figure 8. This way, three subsystems will be obtained
with each subsystem focusing on a single state. While this look simple and is
similar to the parallel decomposition, it require information about the state
in the preceding subsystem. Moreover, the control invariant sets obtained this
way based on the subsystems are not very useful in reconstructing the control
invariant set of the original system since they do not contain the
interconnection information about the neighboring states. The only way to
reconstruct the control invariant set from the subsystem solutions is to find
the Cartesian product of the resulting control invariant set for the each
subsystem. This will ultimately result in a large over approximation of
$R_{X}$.

A more useful way is to decompose system (14) into subsystems such that each
subsystem is chained to its neighboring subsystem as shown in Figure 8. This
way, the interactions between the states can be accounted for. As an example,
system (14) can be decomposed into two subsystems that share a common part.
This can be considered as equivalent to expanding the system to a higher
dimension since the overlapping states are repeated in the state equation as
shown in Equation (15). The overlapping states play an important role in
reconstructing the control invariant set from the distributed computing
results. The role of the overlapping states will be made clear in the later
discussion.

More formally, let $z\in\mathbb{R}^{n_{1}+n_{2}}$ denote the states for the
expanded system for system (11) such that
$z=[x_{1}~{}x_{2}~{}x_{2}~{}x_{3}]^{T}$ is the expanded system state vector
with $z_{1}=[x_{1}~{}x_{2}]^{T}\in\mathbb{R}^{n_{1}}$ and
$z_{2}=[x_{2}~{}x_{3}]^{T}\in\mathbb{R}^{n_{2}}$ being the subsystem states.
The two subsystems are shown in Equation (16). The corresponding decomposition
is depicted in Figure 8.

$z^{+}=\left[\begin{array}[]{cc|cc}A_{11}&0&0&0\\\ A_{21}&A_{22}&0&0\\\
\hline\cr A_{21}&0&A_{22}&0\\\ 0&0&A_{32}&A_{33}\\\ \end{array}\right]z$ (15)

$S_{1}:z_{1}^{+}=f_{1}(z_{1})=\begin{bmatrix}A_{11}&0\\\
A_{21}&A_{22}\end{bmatrix}z_{1}\quad
S_{2}:z_{2}^{+}=f_{2}(z_{2},\tilde{x}_{1})=\begin{bmatrix}A_{22}&0\\\
A_{32}&A_{33}\end{bmatrix}z_{2}+\begin{bmatrix}A_{21}\\\
0\end{bmatrix}\tilde{x}_{1}$ (16)

Notice that while Subsystem 1 is independent of any external state, Subsystem
2 requires information about $x_{1}$ similar to our earlier discussions. We
will refer to the information about $x_{1}$ required by Subsystem 2 as missing
state information, denoted by $\tilde{x}_{1}$, in the next section. How the
missing state information is obtained will be discussed in the subsequent
section. As can be seen in the above equations, the decomposed system has a
higher dimension than the original system i.e. $n_{1}+n_{2}>n$. While the
above example used a linear system, this type of decomposition can be readily
applied to nonlinear systems with exogenous inputs and serves as a precursor
to the proposed algorithm. Unfortunately, there is no general way of
reconstructing the control invariant set for the overall system from the
solution of the subsystems. This is because of the coupling between the two
subsystems and that in Subsystem 2, though $x_{1}$ information is used, the
dynamics is not considered.

##  4 Computation of the largest control invariant set via system
decomposition

In this section, we present an algorithm to compute an approximation of
$R_{X}$ via system decomposition. We focus on the cascade structure presented
in Section 3. Given a cascade system with its associated overlapping
decomposition, the algorithm computes the control invariant set for each of
the subsystems, first in a decentralized manner, followed by a distributed.
Thereafter, the control invariant for the overall system is reconstructed and
subset of the cells tested for control invariance. The central idea employed
here is that the missing state information is treated as both a regular input
and an uncertain input. This way, the two aspects of the missing state is
utilized. As we will show in this section, the use of graph analysis makes it
easier for this to be achieved. Throughout this section, we assume there are
only two subsystems $S_{1}$ and $S_{2}$, that is $M=2$, as described in
Equations (16) with subspaces $X_{1}=\text{proj}_{1}(X)$ and
$X_{2}=\text{proj}_{2}(X)$ respectively.

We begin this section by first describing procedure for the decentralized and
distributed computation. Thereafter, we present the procedure for the
reconstruction and validation of the solution. We end this section by
conducting a brief analyzes of the complexity of the proposed algorithm.

###  4.1 Decentralized and distributed computation

The decentralized and the distributed computation makes use of the standard
GIS algorithm to determine the cells which approximate the largest control
invariant sets ($R_{X_{1}}$ and $R_{X_{2}}$) for each of the subsystems. We
denote the solutions of the subsystems by $R_{1}$ and $R_{2}$ respectively.
During the computations, the directed graphs $G_{1}$ and $G_{2}$ for each
subsystem are constructed based on the collections $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$ of each subsystem respectively. Further analysis on these
graphs is vital for the reconstruction of the solution of the overall system.

In the decentralized step, each subsystem is treated as independent of each
other. This allows for the computation of $R_{1}$ and $R_{2}$ in parallel. The
goal is to quickly obtain an approximation of the largest control invariant
sets and reduce the search space for further computations. To make this
possible, the interval for missing state information for $x_{1}$ in Subsystem
2 is held constant during the computation. In this case, an approximation
which is obtained by projecting the state constraint onto the dimension of
$x_{1}$ may be used, that is

$\tilde{x}_{1}=\text{proj}_{u}X$ (17)

where $u$ denotes the dimension of $x_{1}$. Once the computation begin, there
is no communication between the solutions of the two subsystems. The algorithm
for the decentralized step is summarized in Algorithm 1. The algorithm takes
as input the subsystem equations, the quantized subspaces for each subsystem,
the input constraint and the number of subsystems. As mentioned earlier, this
is fixed since only two subsystems are considered here. The algorithm returns
estimates of the control invariant sets for each of the subsystems as well as
their associated digraphs.

Input: Subsystems in Equation (16) and the collections $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$, $U$, $M=2$

Output: $G_{i}$, $R_{i}$, $i=1,\cdots,M$

$G_{1}\leftarrow$ Construct the directed graph representation of $S_{1}$

$R_{1}\leftarrow$ $I^{+}(G_{1})$

for  _$i=2\cdots M$ _ do

// Can be done in parallel

$\tilde{x}_{i-1}\leftarrow\text{proj}_{i-1}X$

$G_{i}\leftarrow$ Construct the directed graph representation of $S_{i}$
utilizing $\tilde{x}_{i-1}$

$R_{i}\leftarrow$ $I^{+}(G_{i})$

return $R_{i},G_{i},(i=1,\cdots,M)$

Algorithm 1 Decentralized computation

In the distributed step, the control invariant set for each subsystem is
computed sequentially. Thus, the solution for the immediately preceding
subsystem is used to estimate the missing state information for the current
subsystem. While the missing state information for each subsystem is obtained
from the state constraint in the decentralized computation, the missing state
information is dynamically obtained from the neighboring subsystem solution in
the distributed computation. This is the key difference between the
decentralized computation and the distributed computation. It is also the
reason for computing the solution for each subsystem in a sequential fashion.
The procedure for dynamically estimating the missing state information in the
computation of the control invariant set of Subsystem 2 is presented in Figure
9.

Figure 9:  Procedure for dynamically estimating the missing $x_{1}$
information for each cell in Subsystem 2. (a) Solution of Subsystem 1 $R_{1}$.
The yellow cells indicate the cells that correspond to $B_{3}$ in the solution
of Subsystem 2. Merging and projecting these cells onto the $x_{1}$ dimension
produces the range of $x_{1}$ for $B_{3}$. (b) Solution of Subsystem 2
$R_{2}$. The green cell indicate the cell whose missing state information is
being estimated. This procedure is repeated for all other cells in $R_{2}$.

By allowing unidirectional communication through dynamic estimation of the
missing state information, the solution of each subsystem is further refined.
Furthermore, a more realistic graph representation of the dynamics of each
subsystem is obtained. This will be useful during the set reconstruction step
of the algorithm. The algorithm for the distributed computation is presented
in Algorithm 2. It is similar to the algorithm used for the decentralized
computation. As mentioned earlier, the key difference is how the missing state
information $\tilde{x}_{i}$ is obtained.

Input: Subsystems in Equation (16) and the collections $\mathcal{C}_{1}$ and
$\mathcal{C}_{2}$, $U$, $M=2$

Output: $G_{i}$, $R_{i}$, $i=1,\cdots,M$

$G_{1}\leftarrow$ Construct the directed graph representation of $S_{1}$

$R_{1}\leftarrow$ $I^{+}(G_{1})$

for  _$i=2\cdots M$_ do

Estimate the range of $\tilde{x}_{i-1}$ for $S_{i}$ from $R_{i-1}$

$G_{i}\leftarrow$ Construct the directed graph representation of $S_{i}$
utilizing $\tilde{x}_{i-1}$

$R_{i}\leftarrow$ $I^{+}(G_{i})$

return $R_{i},G_{i},(i=1,\cdots,M)$

Algorithm 2 Distributed computation

###  4.2 Set reconstruction and validation

Following the distributed computation step, the solutions from the two
subsystems, that is $R_{1}$ and $R_{2}$, together with the directed graphs
$G_{1}$ and $G_{2}$ are used to reconstruct the solution for the overall
system $R$. Notice that the Cartesian product of the two sets will result in a
4 dimensional object which is not accurate. This is because of the overlapping
decomposition. The procedure for reconstructing $R$ therefore involves
projecting the cells in $R_{1}$ onto the $x_{1}$ dimension in a dynamic
fashion similar to the procedure for estimating the missing state information
in the distributed computation. However, the range of $x_{1}$ is not merged in
this case since the focus is reconstruction of the set for the overall system.
The procedure for reconstructing the set is described in Figure 10.

Figure 10:  Procedure for reconstructing the control invariant set for the
overall system from the subsystem solutions $R_{1}$ and $R_{2}$. The procedure
is indicated for Cell $B_{3}$ in $R_{2}$. Because of the overlap, $B_{6}$ in
$R_{2}$ has a connection with with $B_{1}$, $B_{2}$ and $B_{3}$ in $R_{1}$.
These cells are then projected onto the $x_{1}$ dimension to obtain the range
of $x_{1}$ for each cell. The Cartesian product of the $B_{3}$ with the range
of $x_{1}$ produces 3 dimensional cells (in this case 3 cells). This is
repeated for all other cells in $R_{2}$ to obtain the solution $R$ for the
overall system.

As described earlier, reconstructing the solution of the overall system from
the subsystem solutions may often than not result an approximation error due
to the decomposition and coupling between the subsystems. However, compared to
the naive decomposition, the approximation error for the overlapping
decomposition is expected to be smaller. Nonetheless, the approximation error
need to be addressed to ensure that the reconstructed solution from the
solution of the subsystems converges to solution when the full system model is
used.

To address the approximation error due to the decomposition, we propose that a
subset of the cells forming the reconstructed set needs to be validated for
inclusion in the final reconstructed solution. The validation test involves
computing the next feasible states of each cell to be tested using the overall
system model. Thereafter, the next feasible states of the cell is checked for
intersection with the reconstructed set. A cell $B_{i}$ fails the validation
test if its successive states does not intersect the reconstructed set $R$,
that is

$F(B_{i})\cap R=\emptyset$ (18)

This is a direct intuition from the definition of control invariance. Any cell
that fails the validation test is removed from the reconstructed set.

What now remains is how to select the cells from the reconstructed set $R$ for
the validation test. In principle, each cell in the reconstructed set need to
be tested. However, this may result in having to check for a large number of
cells especially when the system dimension is high. An alternative approach is
to further investigate the graph $G_{2}$ from the distributed computation step
to determine the cells that need to be tested. This way, the cells to test
will be significantly reduced. The general idea we use to further investigate
the graph $G_{2}$ is to consider the adversarial aspects of the missing state
information used to construct $G_{2}$. We know by definition that the set
$R_{2}$ is made of cells with infinite admissible paths passing through it,
that is non-leaving cells. Thus, if there is a direct path from a cell in
$R_{2}$ to a cell not in $R_{2}$, then that path could have been caused by the
missing state $\tilde{x}_{1}$. By doing this, we consider the contribution of
the missing state $\tilde{x}_{1}$ as a regular input and as an uncertain
input.

Given the collection $\mathcal{C}_{2}$ for Subsystem 2, let
$\mathcal{C}^{lc}_{2}$ denote the leaving cells. It is easy to see that
$\mathcal{C}_{2}=R_{2}\cup\mathcal{C}^{lc}_{2}$.

######  Definition 8 (Incoming neighbors).

Given the directed graph $G=(V,E)$, the incoming or in neighbors of the vertex
$u\in V$ are the nodes $v\in V$ such that the edges $(v,u)\in E$ exist.

The central idea is to find the incoming neighbors of the leaving cells in the
graph of Subsystem 2. Then we find those incoming neighbors that have an
intersection with the non-leaving cells. Let $\mathcal{C}^{t}_{2}$ denote the
cells that need to be tested for Subsystem 2 and
in$(G_{2},\mathcal{C}^{lc}_{2})$ be the incoming neighbors of the leaving
cells on $G_{2}$. Then the following relationship describes the cells that
need to be flagged for further testing after the reconstruction

$\mathcal{C}^{t}_{2}=\text{in}(G_{2},\mathcal{C}^{lc}_{2})\cap R_{2}$ (19)

By obtaining the cells to be tested according to Equation (19), we ensure that
any cell that has a direct path to any of the leaving cells are verified. The
cells to be tested, that is $\mathcal{C}^{t}$ for the full dimensional problem
can be obtained from $\mathcal{C}^{t}_{2}$ by following the set reconstruction
procedure described earlier for $R_{2}$.

Figure 11:  Procedure for the analysis of the distributed graphs to find the
cells to be tested.

To illustrate how a cell is selected for testing, Figure 11 is presented. From
the figure, it can be seen that the non-leaving cells of Subsystem 2 $R_{2}$
are $\\{B_{2},\cdots,B_{6}\\}$ while the leaving cells $\mathcal{C}^{lc}_{2}$
are $\\{B_{8},B_{9}\\}$. Also, the incoming neighbors of the leaving cells are
$\\{B_{5},B_{6}\\}$. Thus, $\mathcal{C}^{t}_{d,2}=\\{B_{5},B_{6}\\}$.

The algorithm for the set validation is presented in Algorithm 3. The
algorithm takes as input the full system model, the constraint set, the cells
to be tested $\mathcal{C}^{t}$ and the reconstructed solution $R$. It returns
the final validated solution $R^{*}$. A while loop is used in the validation
algorithm since the order to conduct the tests for the cells in the testing
set in not known in advanced. Hence, we check several times until no cell is
removed.

######  Remark 1.

It is worth mentioning that the presence of both finite and infinite path
passing through a particular cell may be caused by the effects of the actual
system inputs and not the missing state. The effects of this two can be
separated. However, in this work we take a conservative approach and treat the
presence of both finite and infinite admissible path on a cell as an effect of
the missing state.

Input: System (1), constraint sets $X$ and $U$, $\mathcal{C}^{t}$, $R$

Output: $R^{*}$

$cont=true$

while  _cont_ do

$tmp\leftarrow\emptyset$

for  _$B_{i}\in\mathcal{C}^{t}$_ do

if  _$F(B_{i})\cap R=\emptyset$_ then

Add $B_{i}$ to $tmp$

if  _$tmp\neq\emptyset$_ then

$R\leftarrow R/tmp$ // Set difference

$\mathcal{C}^{t}\leftarrow\mathcal{C}^{t}/tmp$

$R^{*}\leftarrow R$

return $R^{*}$

Algorithm 3 Set validation

###  4.3 Computational complexity

In this section, we briefly analyze the computational complexity of the
proposed algorithm. In particular, we focus on how the complexity of the
algorithm increases in the system state dimension $n$. Consider that only two
subsystems $S_{1},S_{2}$ with state dimensions $n_{1}<n,n_{2}<n$ have been
created and let $n_{c}$ denote the number of interval divisions per state
dimension. Based on the results of [8], we know that the overall complexity of
the graph-based control invariant set computation algorithm is determined by
the graph construction step which is $\mathcal{O}(n_{c}^{n})$. This implies
that the algorithm increases exponentially in the state dimension $n$.

Owing to the decomposition, the computation of the control invariant sets for
the each of the subsystems is $\mathcal{O}(n_{c}^{n_{1}})$ and
$\mathcal{O}(n_{c}^{n_{2}})$ respectively which is less than
$\mathcal{O}(n_{c}^{n})$. Thus, the complexity of the decentralized and the
distributed steps is determined by the subsystem with the largest dimension
i.e. max($\mathcal{O}(n_{c}^{n_{1}}),\mathcal{O}(n_{c}^{n_{2}})$). In both of
the steps, the missing states are determined and stored in a dictionary prior
to the computation and therefore retrieving it has a complexity of
$\mathcal{O}(1)$. The centralized step is probably the most time consuming
part of the algorithm since the control invariant set is reconstructed from
the subsystem solution and validated. However, since only the cells flagged
for testing $\mathcal{C}^{t}$ are tested, three possibilities can occur
depending on the number of cells flagged. Let $n_{t}$ be the number of cells
flagged for testing after the distributed computation step. At worst, all the
cells are flagged and therefore all the cells in the reconstruction step need
to be validated. At the best case, no cell is flagged and therefore the
overall control invariant set can be reconstructed without any validation. On
the average case, not all the cells are flagged for validation. The complexity
in this stage is therefore $\mathcal{O}(n_{t})$ which is on average less than
or at worst equal to $\mathcal{O}(n_{c}^{n})$.

##  5 Examples

In this section we present three numerical examples to demonstrate the
efficacy of the proposed algorithm. The first two are 3 dimensional linear and
nonlinear examples. This is used to demonstrate that the solution converges to
the solution of the full system model for both the linear and nonlinear case.
This is because, it is computationally difficult to obtain an approximation of
$R_{X}$ for higher dimensional systems using the standard GIS algorithm.
Finally, we compute an outer approximation of $R_{X}$ for a 6 dimensional
nonlinear system.

Unless otherwise stated, the number of division in each dimension was fixed at
128.The computation was performed on a laptop computer with Intel i7 CPU at
2.60 GHZ and 16 GB RAM. We refer to the solution from the decomposition as the
distributed solution and that from the full system model as the centralized
solution.

###  5.1 Linear system example

Consider the linear time-invariant cascade system

$x^{+}=Ax+Bu$ (20)

where $A$ and $B$ are given by

$A=\begin{bmatrix}2&0&0\\\ 1&2&0\\\
0&1&2\end{bmatrix}~{}\text{and}~{}B=\begin{bmatrix}1\\\ 0\\\ 0\end{bmatrix}.$

The constraints on the states and input are
$X=\\{x\in\mathbb{R}^{3}:\|x\|_{\infty}\leq 5\\}$ and
$U=\\{u\in\mathbb{R}:\|u\|_{\infty}\leq 1\\}$ respectively.

Figure 12:  Comparison of the convex hull of the cells in the centralized and
distributed solutions. Blue shaded area: centralized solution. Dashed red:
solution of the sets obtained from the distributed computation. Black: Final
solution after reconstructing the solution and validating the cells.

Our goal is to compute the control invariant set for the system. To use the
proposed algorithm, the system is first expanded to a 4-dimensional system by
repeating the state equation for the second state as shown in Equation (21).

$z^{+}=f(z,u)=\left[\begin{array}[]{cc|cc}2&0&0&0\\\ 1&2&0&0\\\ \hline\cr
1&0&2&0\\\ 0&0&1&2\\\ \end{array}\right]z+\left[\begin{array}[]{c}1\\\ 0\\\
\hline\cr 0\\\ 0\\\ \end{array}\right]u$ (21)

where $z=[x_{1}~{}x_{2}~{}x_{2}~{}x_{3}]^{T}\in\mathbb{R}^{4}$. It can be seen
that $n_{z}=4>n_{x}=3$. Equation (21) is then decomposed to obtain two
subsystems namely: subsystem one: $z_{1}=(x_{1},x_{2})$ and subsystem two:
$z_{2}=(x_{2},x_{3})$ as shown in Equations (22) and (23) below.

$z_{1}^{+}=f_{1}(z_{1},u)=\begin{bmatrix}2&0\\\
1&2\end{bmatrix}z_{1}+\begin{bmatrix}1\\\
0\end{bmatrix}u=f_{1}(x_{1},x_{2},u)$ (22)

$z_{2}^{+}=f_{2}(z_{2},x_{1})=\begin{bmatrix}2&0\\\
1&2\end{bmatrix}z_{2}+\begin{bmatrix}1\\\
0\end{bmatrix}x_{1}+\begin{bmatrix}0\\\
0\end{bmatrix}u=f_{2}(x_{2},x_{3},x_{1})$ (23)

The control invariant sets for the two subsystems and ultimately the full
system are computed using the distributed algorithm. The solution from the
distributed algorithm is compared to the projections of the control invariant
set computed using the centralized model. This is presented in Figure 12. It
can be seen that the distributed computation converges to the centralized
computation.

###  5.2 Nonlinear system example

Consider the nonlinear system

$\displaystyle x_{1}^{+}=x_{1}^{2}+u$ (24a) $\displaystyle
x_{2}^{+}=x_{2}^{2}+x_{1}$ (24b) $\displaystyle x_{3}^{+}=x_{3}^{2}+x_{2}$
(24c)

The constraints on the states and input are
$X=\\{x\in\mathbb{R}^{3}:\|x\|_{\infty}\leq 5\\}$ and
$U=\\{u\in\mathbb{R}:\|u\|_{\infty}\leq 1\\}$ respectively.

Our goal is to compute an approximation of the largest control invariant set
for the system. To use the proposed algorithm, the system is first expanded to
a 4-dimensional system by repeating the state equation for the second state as
shown in Equation (25).

$z^{+}=f(z,u)=\left[\begin{array}[]{cc|cc}x_{1}^{2}&0&0&0\\\
x_{1}&x_{2}^{2}&0&0\\\ \hline\cr x_{1}&0&x_{2}^{2}&0\\\ 0&0&x_{2}&x_{3}^{2}\\\
\end{array}\right]+\left[\begin{array}[]{c}1.0\\\ 0\\\ \hline\cr 0\\\ 0\\\
\end{array}\right]u$ (25)

where $z=[x_{1}~{}x_{2}~{}x_{2}~{}x_{3}]^{T}\in\mathbb{R}^{n_{z}}$ with
$n_{z}=4>n_{x}=3$. Thereafter, the expanded system is decomposed into two
subsystems with overlapping states as described in the earlier sections.

The results from the computations using the full system model and that of the
decomposed system model are presented in Figure 13 and 14. Figure 13 shows the
convex hull of the sets in the distributed solution and the centralized
solution. Again, it can be seen that the solution using the decomposed models
converges to that of the full system model after the reconstruction.
Furthermore, Figure 14 shows the scalability of the proposed distributed
algorithm compared to the centralized algorithm. It can be seen that the
algorithm using the decomposed model scales better than when the centralized
model is used.

Figure 13:  Comparison of the convex hull of the cells in the solutions
utilizing the decomposed model and the full system model. Blue shaded area:
Solution from the computation utilizing the full system model. Dashed red:
solution of the sets obtained from the distributed computation. Black: Final
solution after reconstructing the solution and validating the cells.  Figure
14:  The computation vs the number of interval divisions in each dimension. At
256 interval divisions, it took more than an hour for the computation using
the centralized model. An hour has been used for better visualization.

###  5.3 Three Continuously stirred tank reactors in series example

Let us consider an isothermal continuous-stirred tank reactor (CSTR) in which
the following irreversible chained reactions occurs

$A\xrightarrow{k_{1}}B\xrightarrow{k_{2}}C$

The reactor can be described by the following dimensionless modeling equations

$\displaystyle\frac{dx_{1}}{dt^{\prime}}={}$ $\displaystyle-
x_{1}+\text{Da}_{1}x_{1}+u_{1}$ (26a)
$\displaystyle\frac{dx_{2}}{dt^{\prime}}={}$ $\displaystyle-
x_{2}+\text{Da}_{2}x_{2}^{2}-\text{Da}_{1}x_{1}$ (26b)

where Da${}_{1}=1$ and Da${}_{2}=2$ represent the dimensionless Damkholer
number for the reactions 1 and 2, $x$ and $u$ are the dimensionless state and
input vectors and $t^{\prime}$ is dimensionless time. A similar version of the
model was presented in [25]. By connecting three of the CSTR model in series
(shown in Figure 15), the following six dimensional model is obtained.

Figure 15:  A schematic diagram of the three CSTR in series.

$\displaystyle\frac{dx_{1}}{dt^{\prime}}={}$ $\displaystyle-
x_{1}+\text{Da}_{1}x_{1}+u_{1}$ (27a)
$\displaystyle\frac{dx_{2}}{dt^{\prime}}={}$ $\displaystyle-
x_{2}+\text{Da}_{2}x_{2}^{2}-\text{Da}_{1}x_{1}$ (27b)
$\displaystyle\frac{dx_{3}}{dt^{\prime}}={}$ $\displaystyle-
x_{3}+\text{Da}_{1}x_{3}+x_{1}$ (27c)
$\displaystyle\frac{dx_{4}}{dt^{\prime}}={}$ $\displaystyle-
x_{4}+\text{Da}_{2}x_{4}^{2}-\text{Da}_{1}x_{3}+x_{2}$ (27d)
$\displaystyle\frac{dx_{5}}{dt^{\prime}}={}$ $\displaystyle-
x_{5}+\text{Da}_{1}x_{5}+x_{3}$ (27e)
$\displaystyle\frac{dx_{6}}{dt^{\prime}}={}$ $\displaystyle-
x_{6}+\text{Da}_{2}x_{6}^{2}-\text{Da}_{1}x_{5}+x_{4}$ (27f)

The continuous-time model is discretized using a step size of 1 before usage
in GIS algorithm. The constraints on the states and input are
$X=\\{x\in\mathbb{R}^{6}:0\leq x\leq 1\\}$ and $U=\\{u\in\mathbb{R}:0\leq
u\leq 1\\}$ respectively. The results of each reactor is presented in Figures
16–18. The computations were performed on a computing server with single core
Intel Xeon E5 processor with 2.1 GHz frequency and 100 GB of RAM. It took
roughly 4 hours to compute the solution using the distributed approach with
the reconstruction and validation step taking about 70 % of the total
computation time. It is worth mentioning that it is intractable to compute an
approximation of the largest control invariant set for the six dimensional
system using the standard GIS algorithm without system decomposition.

Figure 16:  The solution for Reactor 1.  Figure 17:  The solution for Reactor
2. Red: results from decentralized computation. Blue: results from the
distributed computation. Black: Final solution after validation Figure 18:
The solution for Reactor 3. Red: results from decentralized computation. Blue:
results from the distributed computation. Black: Final solution after
validation

##  6 Concluding remarks

Given a constrained discrete-time time-invariant control system, this paper
obtains outer approximations of the largest control invariant set contained in
the state constraint. For this purpose, we presented a graph-based distributed
algorithm which approximates the dynamics of the control system as a directed
graph allowing for analysis of the system using graph theory. This algorithm
paves a promising way to overcome the “curse of dimensionality” encountered in
control invariant set computation. Simulations using several numerical
examples demonstrate the ability and effectiveness of the proposed method to
approximate the largest control invariant set, just like the centralized
algorithm.

Future work will focus on further improving the computational efficiency of
the proposed algorithm by utilizing recent advances in interdependent graph
networks to analyze the distributed graphs for control invariance. This way
the validation step can be conducted without having to reconstruct the set
from the subsystem solutions. Another research direction is appropriate
representation of the set obtained from the algorithm for MPC applications.
The use of machine learning techniques will be explored.

## References

  * [1] T. Homer and P. Mhaskar, “Constrained control lyapunov function-based control of nonlinear systems,” Systems & Control Letters, vol. 110, pp. 55–61, 2017.
  * [2] I. Mitchell, A. Bayen, and C. Tomlin, “A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games,” IEEE Transactions on Automatic Control, vol. 50, pp. 947–957, jul 2005.
  * [3] T. Homer, M. Mahmood, and P. Mhaskar, “A trajectory-based method for constructing null controllable regions,” International Journal of Robust and Nonlinear Control, vol. 30, no. 2, pp. 776–786, 2020.
  * [4] J. P. Aubin, Viability theory.  Modern Birkhäuser classics, Boston: Birkhäuser, 2009.
  * [5] F. Blanchini, “Set invariance in control,” Automatica, vol. 35, no. 11, pp. 1747–1767, 1999.
  * [6] D. Q. Mayne, “Control of constrained dynamic systems,” European Journal of Control, vol. 7, no. 2-3, pp. 87–99, 2001.
  * [7] M. Cannon, V. Deshmukh, and B. Kouvaritakis, “Nonlinear model predictive control with polytopic invariant sets,” Automatica, vol. 39, no. 8, pp. 1487–1494, 2003.
  * [8] B. Decardi-Nelson and J. Liu, “Computing robust control invariant sets of constrained nonlinear systems: A graph algorithm approach,” Computers & Chemical Engineering, vol. 145, p. 107177, 2021.
  * [9] S. V. Raković, B. Kern, and R. Findeisen, “Practical set invariance for decentralized discrete time systems,” in 49th IEEE Conference on Decision and Control (CDC), pp. 3283–3288, IEEE, 2010.
  * [10] M. Rungger and P. Tabuada, “Computing robust controlled invariant sets of linear systems,” IEEE Transactions on Automatic Control, vol. 62, pp. 3665–3670, July 2017.
  * [11] M. Fiacchini, T. Alamo, and E. Camacho, “On the computation of convex robust control invariant sets for nonlinear systems,” Automatica, vol. 46, pp. 1334–1338, Aug. 2010.
  * [12] J. M. Bravo, D. Limón, T. Alamo, and E. F. Camacho, “On the computation of invariant sets for constrained nonlinear systems: An interval arithmetic approach,” Automatica, vol. 41, no. 9, pp. 1583–1589, 2005.
  * [13] D. Bertsekas, “Infinite time reachability of state-space regions by using feedback control,” IEEE Transactions on Automatic Control, vol. 17, pp. 604–613, October 1972.
  * [14] T. Homer and P. Mhaskar, “Utilizing null controllable regions to stabilize input-constrained nonlinear systems,” Computers & Chemical Engineering, vol. 108, pp. 24–30, 2018.
  * [15] A. B. Kurzhanski and P. Varaiya, “Ellipsoidal techniques for reachability analysis: internal approximation,” Systems & control letters, vol. 41, no. 3, pp. 201–211, 2000.
  * [16] G. Frehse, C. Le Guernic, A. Donzé, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang, and O. Maler, “Spaceex: Scalable verification of hybrid systems,” in International Conference on Computer Aided Verification, pp. 379–395, Springer, 2011.
  * [17] M. Chen, S. L. Herbert, M. S. Vashishtha, S. Bansal, and C. J. Tomlin, “Decomposition of reachable sets and tubes for a class of nonlinear systems,” IEEE Transactions on Automatic Control, vol. 63, no. 11, pp. 3675–3688, 2018.
  * [18] S. Riverso, K. Kouramas, and G. Ferrari-Trecate, “Decentralized and distributed robust control invariance for constrained linear systems,” in 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pp. 5978–5984, IEEE, 2017.
  * [19] A. Li and M. Chen, “Guaranteed-safe approximate reachability via state dependency-based decomposition,” in 2020 American Control Conference (ACC), pp. 974–980, IEEE, 2020.
  * [20] B. Decardi-Nelson and J. Liu, “A distributed control invariant set computing algorithm for nonlinear cascade systems,” in Proceedings of the American Control Conference, accepted, 2022.
  * [21] E. C. Kerrigan, Robust constraint satisfaction: Invariant sets and predictive control.  PhD thesis, University of Cambridge, 2001.
  * [22] G. Osipenko, Dynamical Systems, Graphs, and Algorithms.  No. 1889 in Lecture Notes in Mathematics, Berlin ; New York: Springer, 2007.  OCLC: ocm75927357.
  * [23] W. T. Trotter and P. Erdös, “When the cartesian product of directed cycles is hamiltonian,” Journal of Graph Theory, vol. 2, no. 2, p. 137–142, 1978\.
  * [24] S. Kaynama and M. Oishi, “A modified riccati transformation for decentralized computation of the viability kernel under LTI dynamics,” IEEE Transactions on Automatic Control, vol. 58, p. 2878–2892, Nov 2013.
  * [25] W. Harmon Ray, Advanced Process Control.  Chemical Engineering S., New York, NY: McGraw-Hill, Sept. 1980.