# Reachability of Nonlinear Systems with Unknown Dynamics

Taha Shafa and Melkior Ornik Taha Shafa and Melkior Ornik are with the
Department of Aerospace Engineering and the Coordinated Science Laboratory,
University of Illinois Urbana-Champaign, Urbana, USA (e-mail:
<EMAIL_ADDRESS>mornik@illinois.edu).

###### Abstract

Determining the reachable set for a given nonlinear control system is crucial
for system control and planning. However, computing such a set is impossible
if the system’s dynamics are not fully known. This paper is motivated by a
scenario where a system suffers an adverse event mid-operation, resulting in a
substantial change to the system’s dynamics, rendering them largely unknown.
Our objective is to conservatively approximate the system’s reachable set
solely from its local dynamics at a single point and the bounds on the rate of
change of its dynamics. We translate this knowledge about the system dynamics
into an ordinary differential inclusion. We then derive a conservative
approximation of the velocities available to the system at every system state.
An inclusion using this approximation can be interpreted as a control system;
the trajectories of the derived control system are guaranteed to be
trajectories of the unknown system. To illustrate the practical implementation
and consequences of our work, we apply our algorithm to a simplified model of
an unmanned aerial vehicle.

Notice of Previous Publication. This manuscript substantially improves the
work of [1]. Theory has been generalized to include a class of non-invertible
matrices and improved to provide a larger set of reachable states. All lemmas,
corollaries, and Theorems 1, 2, and 4 are entirely novel. Theorem 3 has been
slightly modified from existing theorems in [1] given our new results.

Index terms— Reachable Set Computation, Nonlinear Control Systems, Uncertain
Systems, Aerospace Systems, Autonomous Systems

##  I Introduction

Damage to a control system can cause significant change to its dynamics. In
order to avoid endangering people located in the vicinity of the system, it is
crucial to understand the system’s remaining capabilities. Motivated by
specific examples like an aircraft losing a wing [2] or a UAV becoming damaged
in an urban environment [3], [4] our goal is to conservatively approximate the
unknown system’s set of reachable states [5], [6] while assuming minimal
knowledge about the system dynamics. We call such a set the guaranteed
reachable set (GRS).

The primary contribution of this paper is to provide a meaningful
underapproximation of the GRS of a control-affine system. We assume the only
available information at the time of computation consists of (i) local
dynamics at a single point, which can be obtained with an arbitrarily small
error from applying test control inputs over a short period of time [7], and
(ii) Lipschitz bounds on the rate of change of the system’s dynamics provided
by prior knowledge of the system design and physical laws. The reachable set
of an unknown system is impossible to compute. Prior work determining the
reachable set often focused on overapproximations [7],[8]. Without discussing
reachable sets, the work in [7], operating under similar assumptions as our
paper, focused on optimistic reachability, i.e., attempting to reach a
particular objective while there exists any chance of reaching it. Conversely,
this paper computes states that are guaranteed to be reachable using
admissible control signals.

Apart from [7], work on reachability under uncertainty considered computation
of reachable sets with dynamics generated by a finite number of uncertain
parameters [9], [10] or having bounded disturbances [11], [12]. Work in
adaptive and robust control [8], [13], [14] assumes more knowledge on the
magnitude or structure of the system dynamics, and classical data-driven
learning methods [2], [15], [16] collect data through repeated system runs. In
contrast, this paper contains substantially fewer assumptions and focuses on
deriving as much information as possible for the GRS based on one system run.

Our approach relies on the interpretation of a control system as a
differential inclusion [17], [18], [19], [20] whose right hand side equals the
set of velocities that the system can achieve at every state in the state
space. While, for an unknown system, exact velocity sets may not be available
anywhere, we can determine the family of all velocity sets that are consistent
with prior knowledge of the local dynamics at a single point and Lipschitz
bounds on the rate of change of the system’s dynamics. The intersection of all
elements of such a family is defined as the guaranteed velocity set (GVS).
Such a set is difficult to express in closed form, however we can analytically
derive its underapproximation to compute an underapproximation of the GRS.

The outline of the paper is as follows: we discuss the problem statement in
greater detail in Section II. We then proceed to derive two simply expressible
sets in Section III, with one set being a ball, and the other being a more
complex convex set, later utilized to derive a polygon such that both the ball
and polygon are contained within the GVS. Next, in Section IV, we use these
sets to derive two classes of simple control-affine systems whose reachable
sets are contained in the GRS. Lastly, we present numerical examples in
Section V with a brief discussion on the implementation of our method.

###  I-A Notation

The set of all matrices with $n$ rows and $m$ columns is denoted by
$\mathbb{R}^{n\times m}$. For any vector $v$, $\|v\|$ denotes its Euclidean
norm and $\|v\|_{1}$ denotes its 1-norm. For any matrix $M$, $M^{T}$ denotes
its transpose and $\|M\|$ denotes its Euclidean norm: $\|M\|=$ max${\|v\|=1}$
$\|Mv\|$. Equivalently, $\|M\|=\sigma_{1}(M)$, where $\sigma_{i}(M)$
represents the $i$-th singular value of $M$. We also let $M^{\dagger}$ denote
the Moore-Penrose pseudoinverse, $\text{Im}(M)$ denote the image (range space)
of $M$, and Ker$(M)$ denote the kernel (null space) of $M$. For matrices $M$
and $N$, we will say $M\in\text{Imm}(N)$ if $M=NP$ for some matrix $P$.
Notation $\mathbb{B}^{n}(a;b)$ denotes a closed ball in $\mathbb{R}^{n}$
centered at $a\in\mathbb{R}^{n}$ with radius $b\geq 0$ under the Euclidean
norm. Set $C_{L}(\mathbb{R}^{n};\mathbb{R}^{m})$ denotes the set of all
functions $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$ with a Lipschitz constant $L$,
i.e., the set of all functions $f$ that satisfy $\|f(x)-f(y)\|\leq
L\|x-y\|\hskip 2.84526pt\text{for all}\hskip 2.84526ptx$,
$y\in\mathbb{R}^{n}$. Notation $a+B\mathcal{X}$ where $a\in\mathbb{R}^{n}$,
$B\in\mathbb{R}^{n\times m}$, and $\mathcal{X}\subseteq\mathbb{R}^{m}$ denotes
the set $a+B\mathcal{X}=\\{a+Bx~{}|~{}x\in\mathcal{X}\\}$.

##  II Problem Statement

Throughout the paper, we attempt to meaningfully underapproximate the
reachable set of a nonlinear control-affine system $\mathcal{M}(f,G)$ defined
by

$\dot{x}(t)=f(x(t))+G(x(t))u(t),\indent x(0)=x_{0},$ (1)

where all $t\geq 0$, $x(t)\in\mathbb{R}^{n}$, admissible inputs
$u(t)\in\mathcal{U}=\mathbb{B}^{m}(0;1)$, which is a common setting in
reachability analysis [21], [22], and functions
$f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, and
$G:\mathbb{R}^{n}\to\mathbb{R}^{n\times m}$ are globally Lipschitz continuous
with Lipschitz constants $L_{f}\geq 0$ and $L_{G}\geq 0$, i.e., $f\in
C_{L_{f}}(\mathbb{R}^{n};\mathbb{R}^{n})$ and $G\in
C_{L_{G}}(\mathbb{R}^{n};\mathbb{R}^{n\times m})$. Cases where $L_{f}=0$ or
$L_{G}=0$ are simple, thus for the remainder of the paper we assume $L_{f}>0$
and $L_{G}>0$. Noting that any sufficiently smooth function is globally
Lipschitz continuous on any compact set, the theory developed in this paper
can also be applied to a system whose states are guaranteed to be bounded,
which is naturally true for a large class of systems [3], [4]. Without loss of
generality, we assume $x_{0}=0$.

###  II-A Assumptions and Technical Requirements

In order to approximate the reachable set of (1), the technical work of [1]
requires full actuation at $x_{0}=0$, i.e., $m=n$ with $G(0)$ being full rank.
We relax these requirements, so the system $\mathcal{M}(f,G)$ need not be
fully actuated at $x_{0}=0$, that is, $m$ does not necessarily equal $n$.
Instead of assuming full actuation, we use the following assumption:

Assumption 1: Functions $f$ and $G$ are of the form $f(x)=Rr(x)$ and
$G(x)=RH(x)$ where $R\in\mathbb{R}^{n\times m}$, $r(x)\in\mathbb{R}^{m}$, and
$H(x)\in\mathbb{R}^{m\times m}$ such that $H(0)$ is invertible.

The case of full actuation in [1] corresponds the case where $R$ = $I_{n}$ in
Assumption 1. Motivated by the online learning technique introduced in [7], we
make the following assumption about the knowledge regarding the system
dynamics.

Assumption 2: Bounds $L_{f}$ and $L_{G}$ are known, as well as values $f(0)$
and $G(0)$ such that $G(0)\neq 0$.

Note we are not assuming any knowledge about matrix $R$; we only assume such
an $R$ exists. We only consider the case where $G(0)\neq 0$ because it is
otherwise impossible to determine anything about guaranteed velocities at
states $x\neq 0$.

It is easily shown that $\text{Im}(G(0))=\text{Im}(R)$ under the conditions of
Assumption 1. For our future results, it is important to determine the set of
$x$ such that $\text{Im}(G(0))=\text{Im}(G(x))$. Crucially, we show that the
images of $G(x)$ and $G(0)$ are equal in some neighborhood of $0$.

Lemma 1: Under Assumptions 1 and 2, if
$\|x\|<\frac{\|G(0)^{\dagger}\|^{-1}}{L_{G}}$, then
$\text{Im}(G(x))=\text{Im}(G(0))$.

###### Proof.

Weyl’s inequality for singular values dictates that singular values as
functions on matrices are uniformly Lipschitz with respect to the operator
norm [23]: $\|\sigma_{s}(G(x))-\sigma_{s}(G(0))\|\leq\|G(x)-G(0)\|\leq
L_{G}\|x\|$, such that $1\leq s\leq r$ with $r=\text{rank}(G(0))$. The
Eckhart-Young-Mirsky theorem [24] along with the singular value decomposition
of $G(0)$ show that $\sigma_{r}(G(0))=\|G(0)^{\dagger}\|^{-1}$ is the smallest
non-zero singular value of $G(0)$. For $x$ that satisfies
$\|x\|<\frac{\|G(0)^{\dagger}\|^{-1}}{L_{G}}$, we thus have
$\|\sigma_{s}(G(x))-\sigma_{s}(G(0))\|<\sigma_{r}(G(0))$, i.e.,
$\sigma_{s}(G(x))>0$. Therefore rank$(G(x))\geq\text{rank}(G(0))$.

We defined $G(x)=RH(x)$, so $\text{Im}(G(x))\subset\text{Im}(R)$. Since
$\text{Im}(G(0))=\text{Im}(R)$, then $\text{Im}(G(x))\subset\text{Im}(G(0))$.
Knowing rank$(G(x))\geq\text{rank}(G(0))$, we conclude
$\text{Im}(G(x))=\text{Im}(G(0))$.

∎

###  II-B Guaranteed Reachable Set

Let us denote a set $\mathcal{D}_{con}\subseteq
C_{L_{f}}(\mathbb{R}^{n};\mathbb{R}^{n})\times
C_{L_{G}}(\mathbb{R}^{n};\mathbb{R}^{n\times m})$ as the set of all pairs $f$
and $G$ consistent with Assumptions 1 and 2. We want to underapproximate the
set of reachable states given solely the knowledge of $\mathcal{D}_{con}$. We
first define the (forward) reachable set
$\mathcal{R}^{\hat{f},\hat{G}}(T,x_{0})=\\{\phi^{\hat{f},\hat{G}}_{u}(t;x_{0})\hskip
2.84526pt|\hskip 2.84526ptu:[0,T]\to\mathcal{U},t\in[0,T]\\}$, where
$\phi^{\hat{f},\hat{G}}_{u}(\cdot;x_{0})$ denotes the controlled trajectory of
the system $\mathcal{M}(\hat{f},\hat{G})$ with control signal $u$ and
$\phi^{\hat{f},\hat{G}}_{u}(0;x_{0})=0$.

Let $T\geq 0$. We describe the guaranteed reachable set (GRS) as:

$\mathcal{R}^{\mathcal{G}}(T,x_{0})=\bigcap_{(\hat{f},\hat{G})\hskip
2.84526pt\in\hskip
2.84526pt\mathcal{D}_{con}}\mathcal{R}^{\hat{f},\hat{G}}(T,x_{0}).$ (2)

The GRS describes the set of all states that are reachable by any system
consistent with our knowledge of the system dynamics.

Problem 1: Determine or underapproximate the GRS.

To solve Problem 1, we first represent ordinary differential equations with
control inputs as an ordinary differential inclusion (ODI). We discuss doing
so in Section III. Given the assumed knowledge of the system dynamics, we
develop underapproximations to the right hand side of this inclusion. In
Section IV, we use these underapproximations to derive two control-affine
systems whose reachable sets are subsets of $\mathcal{R}^{\mathcal{G}}(T,0)$.

##  III Guaranteed Velocities

We follow the classical approach of interpreting ordinary differential
equations with control inputs as inclusions [17], [18], [19], [20]. In this
section, we formally define the guaranteed velocity set of an unknown control
system and determine analytically computable underapproximations of such a
set.

###  III-A Guaranteed Velocity Set

We define the available velocity set of the system $\mathcal{M}(f,G)$ at state
$x$ by $\mathcal{V}_{x}=f(x)+G(x)\mathcal{U}$, and introduce the following
ODI:

$\dot{x}\in\mathcal{V}_{x}=f(x)+G(x)\mathcal{U},\indent x(0)=x_{0}.$ (3)

If a trajectory $\phi(\cdot;x_{0})$ satisfies (3), then it obviously serves as
a solution to the control system (1) for an admissible control input, and vice
versa. Given Assumption 2, set $\mathcal{V}_{x_{0}}=\mathcal{V}_{0}$ is known.
The goal of this section is to provide an underapproximation for set
$\mathcal{V}_{x}$ using sets $\mathcal{D}_{con}$ and $\mathcal{V}_{0}$. We
first define the guaranteed velocity set (GVS) below:

$\mathcal{V}^{\mathcal{G}}_{x}=\bigcap_{(\hat{f},\hat{G})\hskip
2.84526pt\in\hskip
2.84526pt\mathcal{D}_{con}}\hat{f}(x)+\hat{G}(x)\mathcal{U}\subseteq\mathcal{V}_{x}.$
(4)

The GVS $\mathcal{V}^{\mathcal{G}}_{x}$ is the set of all velocities that can
be taken by all systems consistent with the assumed knowledge of the dynamics.
Let us consider the following ODI:

$\dot{x}\in\mathcal{V}^{\mathcal{G}}_{x},\indent x_{0}=0.$ (5)

If $\mathcal{V}^{\mathcal{G}}_{\phi(T;x_{0})}=\emptyset$, we will consider by
convention that the trajectory of (5) ceases to exist at time $T$. The
following proposition then holds directly from (2) and (4).

Proposition 1: Let $T\geq 0$. If a trajectory $\phi$ :
$[0,+\infty)\to\mathbb{R}^{n}$ satisfies (5) at all times $t\leq T$, then
$\phi(T)\in\mathcal{R}^{\mathcal{G}}(T,x_{0})$.

Proposition 1 implies that the reachable set of (5) is a subset of
$\mathcal{R}^{\mathcal{G}}(T,x_{0})$. As briefly described in [1], these sets
are not necessarily equal; establishing conditions for the equality of the
reachable set of (5) and $\mathcal{R}^{\mathcal{G}}(T,x_{0})$ is an open
problem for future work.

###  III-B Ball Underapproximation

Proposition 1 motivates us to underapproximate
$\mathcal{R}^{\mathcal{G}}(T,x_{0})$ by determining the reachable set of (5).
We start by examining the geometry of the ODI. Given $x\in\mathbb{R}^{n}$, our
previous assumptions show that

$\begin{gathered}\\{(\hat{f}(x),\hat{G}(x))\hskip 2.84526pt|\hskip
2.84526pt(\hat{f},\hat{G})\in\mathcal{D}_{con}\\}=\\\
(\mathbb{B}^{n}(f(0);L_{f}\|x\|)\cap\text{Im}(R))\times(\mathbb{B}^{n\times
m}(G(0);L_{G}\|x\|)\cap\text{Imm}(R)).\end{gathered}$ (6)

Given that $\mathcal{U}=\mathbb{B}^{m}(0;1)$, $\mathcal{V}^{\mathcal{G}}_{x}$
is an intersection of infinitely many ellipsoids $a+B\mathcal{U}$, where
$a\in\mathbb{B}^{n}(f(0);L_{f}\|x\|)\cap\text{Im}(R)$ and
$B\in\mathbb{B}^{n\times m}(G(0);L_{G}\|x\|)\cap\text{Imm}(R)$. An
intersection of infinitely many ellipsoids is generally not a geometrically
simple object [25]. Thus, we will determine an underapproximation of
$\mathcal{V}^{\mathcal{G}}_{x}$. Our approach will be to implicitly exploit
convexity of $\mathcal{V}^{\mathcal{G}}_{x}$ and underapproximate the distance
of the boundary of set $\mathcal{V}^{\mathcal{G}}_{x}$ from $0$ in every
direction. In Theorem 1, we calculate one such underapproximation.

Theorem 1: Let $\mathcal{U}$, $L_{f}$, and $L_{G}$ be defined as above. Let
$x\in\mathbb{R}^{n}$ satisfy $(L_{f}+L_{G})\|x\|<\|G(0)^{\dagger}\|^{-1}$.
Define

$\bar{\mathcal{V}}^{\mathcal{G}}_{x}=\mathbb{B}^{n}(f(0);\|G(0)^{\dagger}\|^{-1}-L_{f}\|x\|-L_{G}\|x\|)\hskip
2.84526pt\cap\hskip 2.84526pt\text{Im}(G(0)).$ (7)

Then,
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}\subseteq\mathcal{V}^{\mathcal{G}}_{x}$.

###### Proof.

Set
$\cap_{(\hat{f},\hat{G})\in\mathcal{D}_{con}}\hat{f}(x)-f(0)+\hat{G}(x)\mathcal{U}\text{~{}equals~{}}\mathcal{V}^{\mathcal{G}}_{x}-f(0)$.
On the other hand,
$\cap_{(\hat{f},\hat{G})\in\mathcal{D}_{con}}\hat{f}(x)-f(0)+\hat{G}(x)\mathcal{U}=\cap_{(\tilde{f},\hat{G})\in\tilde{\mathcal{D}}_{con}}\tilde{f}(x)+\hat{G}(x)\mathcal{U}$,
where $\tilde{\mathcal{D}}_{con}$ is defined the same as before, just with the
assumption that $f(0)=0$. Thus, we can assume without loss of generality that
$f(0)=0$.

Let $d\in\text{Im}(G(0))=\text{Im}(R)$ such that $\|d\|=1$. We will prove that
if $|k|\leq\|G(0)^{\dagger}\|^{-1}-L_{f}\|x\|-L_{G}\|x\|$, then equation

$k\cdot d=\hat{f}(x)+\hat{G}(x)u,$ (8)

where $(\hat{f},\hat{G})\in\mathcal{D}_{con}$, admits a solution
$u\in\mathcal{U}=\mathbb{B}^{m}(0;1)$.

We subtract $\hat{f}(x)$ from both sides of (8). Since
$\hat{f}(x)\in\text{Im}(R)$ by Assumption 1 and $kd\in\text{Im}(R)$ by
definition, then $kd-\hat{f}(x)\in\text{Im}(R)$. Also,
$\text{Im}(\hat{G}(x))=\text{Im}(R)$ by Lemma 1. Hence, there exists a vector
$\bar{u}\in\mathbb{R}^{m}$ such that $kd-\hat{f}(x)=\hat{G}(x)\bar{u}$. Now,
through the rank-nullity theorem [24], we can write $\bar{u}=u+u_{2}$ where
$u\in\text{Im}(\hat{G}(x)^{T})$ and $u_{2}\in\text{Ker}(\hat{G}(x))$. Thus,
$\hat{G}(x)\bar{u}=\hat{G}(x)(u+u_{2})=\hat{G}(x)u$; hence,
$kd-\hat{f}(x)=\hat{G}(x)u$.

We multiply both sides of $kd-\hat{f}(x)=\hat{G}(x)u$ on the left by
$\hat{G}(x)^{\dagger}$, resulting in
$\hat{G}(x)^{\dagger}(kd-\hat{f}(x))=\hat{G}(x)^{\dagger}\hat{G}(x)u$. The
term $\hat{G}(x)^{\dagger}\hat{G}(x)u$ results in the projection of $u$ onto
the $\text{Im}(\hat{G}(x)^{T})$ [24]. Given that
$u\in\text{Im}(\hat{G}(x)^{T})$, by definition of a projection,
$\hat{G}(x)^{\dagger}(kd-\hat{f}(x))=\hat{G}(x)^{\dagger}\hat{G}(x)u=u$. Thus,
if we prove that:

$\|\hat{G}(x)^{\dagger}(k\cdot d-\hat{f}(x))\|\leq 1,$ (9)

we will have $\|u\|\leq 1$, i.e., $u\in\mathcal{U}$. Utilizing $\|d\|=1$ along
with the product and triangle inequalities for matrices, we arrive at (10) and
(11):

$\displaystyle\|\hat{G}(x)^{\dagger}(k\cdot
d-\hat{f}(x))\|\leq|k|\|\hat{G}(x)^{\dagger}~{}d\|+\|\hat{G}(x)^{\dagger}\hat{f}(x)\|,$
(10)
$\displaystyle\leq|k|\|\hat{G}(x)^{\dagger}\|+\|\hat{G}(x)^{\dagger}\|\|\hat{f}(x)\|.$
(11)

From (11) it follows that the set of all $k$ that satisfy
$|k|\|\hat{G}(x)^{\dagger}\|+\|\hat{G}(x)^{\dagger}\|\|\hat{f}(x)\|\leq 1$ is
a subset of all $k$ that satisfy (9). In other words, if:

$|k|\leq\|\hat{G}(x)^{\dagger}\|^{-1}-\|\hat{f}(x)\|,$ (12)

then $k$ satisfies (9). We note $\|\hat{G}(x)^{\dagger}\|\neq 0$ from the
definition of the Moore-Penrose pseudoinverse and because $G(x)\neq 0$ from
Lemma 1.

By Weyl’s inequality for singular values [23] and Assumption 2, we obtain the
following inequalities:
$\|\hat{G}(x)^{\dagger}\|^{-1}\geq\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|$ and
$\|\hat{f}(x)\|\leq L_{f}\|x\|$. Thus, since we assumed that $k$ satisfies:

$|k|\leq\|G(0)^{\dagger}\|^{-1}-L_{f}\|x\|-L_{G}\|x\|,$ (13)

it satisfies (8).

∎

We slightly generalize Theorem 1 by also considering all $x$ which satisfy
$(L_{f}+L_{G})\|x\|=\|G(0)^{\dagger}\|^{-1}$.

Corollary 1: If $\|x\|=\frac{\|G(0)^{\dagger}\|^{-1}}{L_{f}+L_{G}}$, then
$f(0)\in\mathcal{V}^{\mathcal{G}}_{x}$.

###### Proof.

Let us take $x$ such that $\|x\|=\frac{\|G(0)^{\dagger}\|^{-1}}{L_{f}+L_{G}}$
and a sequence $x_{1},x_{2},...$ such that
$\|x_{i}\|<\frac{\|G(0)^{\dagger}\|^{-1}}{L_{f}+L_{G}}$ and $\|x_{i}\|\to x$
as $i\to\infty$. Theorem 1 shows that for all $i$, there exists a
$u_{i}\in\mathcal{U}$ such that $f(0)=\hat{f}(x_{i})+\hat{G}(x_{i})u_{i}$.
Thus, because $\mathcal{U}$ is a compact set, there exists a subsequence
$u_{p_{1}},u_{p_{2}},...$ which converges to some $u^{*}\in\mathcal{U}$ [26].
Since $x_{p_{i}}\to x$ and $u_{p_{i}}\to u^{*}$, we have
$f(0)=\hat{f}(x_{p_{i}})+\hat{G}(x_{p_{i}})u_{p_{i}}\to\hat{f}(x)+\hat{G}(x)u^{*}$.
Therefore, $f(0)=\hat{f}(x)+\hat{G}(x)u^{*}$.

∎

###  III-C Advanced Convex Underapproximation

For all $x$ that satisfy
$\|x\|\leq\frac{\|G(0)^{\dagger}\|^{-1}}{L_{f}+L_{G}}$, we now have a set
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ guaranteed to be a subset of
$\mathcal{V}^{\mathcal{G}}_{x}$. Such a set is a projection of a ball onto
$\text{Im}(G(0))$. However, there are instances, particularly when singular
values of $G(0)$ are far apart, where a ball may be a poor underapproximation
of $\mathcal{V}^{\mathcal{G}}_{x}$. Consequently, we derive a new
underapproximated set $\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ in Theorem
2.

Theorem 2: Let $\mathcal{U}$, $L_{f}$, and $L_{G}$ be defined as above. Let
$\mu=1$ if rank($G(0))=m=n$, $\mu=\sqrt{2}$ if rank($G(0))=\text{min}(m,n)$
and $m\neq n$, $\mu=\frac{1+\sqrt{5}}{2}$ if rank($G(0))<\text{min}(m,n)$, and
let $x$ satisfy $(L_{f}+L_{G})\|x\|\leq\|G(0)^{\dagger}\|^{-1}$. If

$\begin{gathered}\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}=\\{f(0)+kd\hskip
2.84526pt|\hskip 2.84526pt\|d\|=1,d\in\text{Im}(R),\hskip 2.84526pt0\leq k\leq
K(d)\\}\\\ \text{s.t.}\hskip
2.84526ptK(d)=\frac{\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|-L_{f}\|x\|}{\|G(0)^{\dagger}d\|(\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|)+\mu\|G(0)^{\dagger}\|L_{G}\|x\|},\end{gathered}$
(14)

then
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}\subseteq\mathcal{V}^{\mathcal{G}}_{x}$.

###### Proof.

As in Theorem 1, we will show that for $k$ and $d$ given in (14), equation (8)
admits a solution $u\in\mathcal{U}$. Like in the proof of Theorem 1, without
loss of generality, we set $f(0)=0$. From inequality (10), it follows that the
set of all $k$ that satisfy

$|k|\|\hat{G}(x)^{\dagger}d\|+\|\hat{G}(x)^{\dagger}\|\|\hat{f}(x)\|\leq 1$
(15)

is a subset of all $k$ that satisfy (9). The term $\|G(x)^{\dagger}\|$ is the
inverse of the smallest nonzero singular value of $G(x)$. Therefore, using
Weyl’s inequality for singular values [23], we have
$\|\hat{G}(x)^{\dagger}\|\leq(\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|)^{-1}$. By
Assumption 2, we also have $\|\hat{f}(x)\|\leq L_{f}\|x\|$. We conclude that
the set of all $k$ that satisfy

$|k|\|\hat{G}(x)^{\dagger}d\|+\frac{L_{f}\|x\|}{\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|}\leq
1$ (16)

is a subset of all $k$ that satisfy (9).

Next, we bound $\|\hat{G}(x)^{\dagger}d\|$ using the product and triangle
inequalities for matrices:

$\begin{gathered}\begin{aligned}
\|\hat{G}(x)^{\dagger}d\|\leq\|G(0)^{\dagger}d\|+\|(\hat{G}(x)^{\dagger}-G(0)^{\dagger})d\|\\\
\leq\|G(0)^{\dagger}d\|+\|\hat{G}(x)^{\dagger}-G(0)^{\dagger}\|.\end{aligned}\end{gathered}$
(17)

By Lemma 1, $\text{Im}(G(0))=\text{Im}(\hat{G}(x))$. Therefore, we can apply
the inequality
$\|\hat{G}(x)^{\dagger}-G(0)^{\dagger}\|\leq\mu\|\hat{G}(x)^{\dagger}\|\|G(0)^{\dagger}\|\|\hat{G}(x)-G(0)\|$
(Theorem 3.3 in [27]). We can now rewrite the upper bound on
$\|\hat{G}(x)^{\dagger}d\|$:

$\|\hat{G}(x)^{\dagger}d\|\leq\|G(0)^{\dagger}d\|+\mu\|\hat{G}(x)^{\dagger}\|\|G(0)^{\dagger}\|\|\hat{G}(x)-G(0)\|.$
(18)

We again use
$\|\hat{G}(x)^{\dagger}\|\leq(\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|)^{-1}$.
According to Weyl’s inequality for singular values, $\|\hat{G}(x)-G(0)\|\leq
L_{G}\|x\|$. Therefore, all $\|\hat{G}(x)^{\dagger}d\|$ that satisfy (18) will
also satisfy

$\|\hat{G}(x)^{\dagger}d\|\leq\|G(0)^{\dagger}d\|+\frac{\mu\|G(0)^{\dagger}\|L_{G}\|x\|}{\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|}.$
(19)

By plugging (19) into (16), we obtain all $k$ that satisfy

$\begin{split}|k|\hskip
2.84526pt\left(\|G(0)^{\dagger}d\|+\frac{\mu\|G(0)^{\dagger}\|L_{G}\|x\|}{\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|}\right)\\\
+\frac{L_{f}\|x\|}{\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|}\leq 1,\end{split}$

i.e.,

$\begin{gathered}|k|\leq\frac{\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|-L_{f}\|x\|}{\|G(0)^{\dagger}d\|(\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|)+\mu\|G(0)^{\dagger}\|L_{G}\|x\|},\end{gathered}$
(20)

also satisfy (9). ∎

By rewriting the denominator of $K(d)$ in Theorem 2 as
$\|G(0)^{\dagger}\|^{-1}\|G(0)^{\dagger}d\|+(\mu\|G(0)^{\dagger}\|-\|G(0)^{\dagger}d\|)L_{G}\|x\|$,
we can easily see that it cannot be negative since $\mu\geq 1$ and
$\|G(0)^{\dagger}\|\geq\|G(0)^{\dagger}d\|$. We now compare the derived sets
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ and
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ in Corollary 2.

Corollary 2: For invertible matrices $G(0)$,
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}\subseteq\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$.

###### Proof.

Given the numerator of $K(d)$ in (14) is identical to the radius of
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ derived in Theorem 1,
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ is contained within
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ if and only if

$\|G(0)^{\dagger}d\|(\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|)+\mu\|G(0)^{\dagger}\|L_{G}\|x\|\leq
1.$ (21)

Inequality (21) is obviously equivalent to

$\|G(0)^{\dagger}d\|(\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|)\leq
1-\mu\|G(0)^{\dagger}\|L_{G}\|x\|.$ (22)

Dividing by $\|G(0)^{\dagger}\|^{-1}-L_{G}\|x\|>0$ results in

$\begin{gathered}\|G(0)^{\dagger}d\|\leq\frac{\|G(0)^{\dagger}\|(1-\mu\|G(0)^{\dagger}\|L_{G}\|x\|)}{1-\|G(0)^{\dagger}\|L_{G}\|x\|}.\end{gathered}$
(23)

We see in the case of invertible $G(0)$ where $\mu=1$, (23) reduces to
$\|G(0)^{\dagger}d\|\leq\|G(0)^{\dagger}\|$, which holds true for all d such
that $\|d\|=1$. Therefore, for invertible $G(0)$,
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ is contained within
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$.

∎

Figure 1: Velocity sets plotted for a system with $f(0)=\begin{bmatrix}0\\\
0\end{bmatrix}$, $G(0)=\begin{bmatrix}10&3\\\ 2&7\end{bmatrix}$, $\|x\|=1$,
and $L_{f}=L_{G}=1$. The blue curves represent the boundaries of the available
velocity sets for all systems such that
$(\hat{f},\hat{G})\in\mathcal{D}_{con}$. The intersection of these sets
produces the guaranteed velocity set $\mathcal{V}^{\mathcal{G}}_{x}$ (white).
Underapproximations $\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ (bounded in red) and
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ (bounded in black) are both
contained in $\mathcal{V}^{\mathcal{G}}_{x}$.

Our calculations result in two derived sets. Corollary 2 proves
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}\subseteq\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$
when $G(0)$ is invertible. Hence, Theorem 2 provides a better
underapproximation than the one obtained in [1], under the assumptions present
in that paper. Figure 1 illustrates a simple example of the GVS
$\mathcal{V}^{\mathcal{G}}_{x}$ and its approximations obtained by Theorems 1
and 2. Consistent with the results of Corollary 2, since $G(0)$ is invertible,
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}\subseteq\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$.

For general cases,
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}\not\subseteq\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$.
Thus, we can take the union of both of sets to determine a larger set of
guaranteed velocities. We denote such a set by

$\hat{\mathcal{V}}^{\mathcal{G}}_{x}=\bar{\mathcal{V}}^{\mathcal{G}}_{x}\cup\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}.$

In Section IV, we show how $\hat{\mathcal{V}}^{\mathcal{G}}_{x}$ can be used
to identify a polygon that can generate a control system with solutions that
satisfy (1).

##  IV Reachable Set

In this section, we aim to utilize the derived sets
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ and
$\hat{\mathcal{V}}^{\mathcal{G}}_{x}$ to determine a set of trajectories
guaranteed to satisfy (1). First, we use $\hat{\mathcal{V}}^{\mathcal{G}}_{x}$
to identify a polygon contained within $\mathcal{V}^{\mathcal{G}}_{x}$.

Lemma 2: Let $\mathcal{S}(x)$ be any finite set of points on the boundary of
$\hat{\mathcal{V}}^{\mathcal{G}}_{x}$. Let $P(\mathcal{S}(x))$ be a convex
hull of $\mathcal{S}(x)$. Every solution to $\dot{x}\in P(\mathcal{S}(x))$ is
a solution to (1).

###### Proof.

The guaranteed velocity set is the intersection of an infinitely many
ellipsoids; such a set is convex [25]. Given
$\hat{\mathcal{V}}^{\mathcal{G}}_{x}\subseteq\mathcal{V}^{\mathcal{G}}_{x}$,
we know $P(\mathcal{S}(x))\subseteq\mathcal{V}^{\mathcal{G}}_{x}$. ∎

While Lemma 2 permits us to chose any set of points on the boundary of
$\hat{\mathcal{V}}^{\mathcal{G}}_{x}$ as $\mathcal{S}(x)$, in the remainder of
the paper, we will choose the points along the singular vectors of $G(0)$. We
do so by noting from inequality (20) that the smaller the magnitude of
$\|G(0)^{\dagger}d\|$, the larger the corresponding $K(d)$. The smallest
magnitude of $\|G(0)^{\dagger}d\|$ for $\|d\|=1$ will be obtained when $d$ is
the singular vector corresponding to the largest singular value of $G(0)$
[24]. We choose other points to be along other singular vectors because of the
orthogonality of singular vectors.

Theorem 1 shows that the reachable set of

$\dot{x}\in\bar{\mathcal{V}}^{\mathcal{G}}_{x},\indent x(0)=x_{0},$ (24)

is a subset of $\mathcal{R}^{\mathcal{G}}_{x}(T,x_{0})$, while Lemma 2 shows
the same for the reachable set of

$\dot{x}\in P(\mathcal{S}(x)),\indent x(0)=x_{0}.$ (25)

###  IV-A Underapproximated Control System – Ball

Analogous to the interpretation of dynamics (1) as an ODI (3), inclusion (24)
can be interpreted as a control system

$\dot{x}=a+g(\|x\|)u,\indent x(0)=x_{0},$ (26)

on $\\{x~{}|~{}\|x\|\leq\|G(0)^{\dagger}\|/(L_{f}+L_{G})\\}$, with $a=f(0)$,
$u\in\mathcal{U}=\mathbb{B}^{m}(0;1)\cap\text{Im}(G(0)^{T})$, and where
$g(s)=\|G(0)^{\dagger}\|^{-1}-(L_{f}+L_{G})s$ if
$s\leq\|G(0)^{\dagger}\|^{-1}/(L_{G}+L_{f})$. We thus obtain the following
result.

Theorem 3: Let $\bar{\mathcal{R}}(T,x_{0})$ be defined as the reachable set of
(26) at time $T$. Then,
$\bar{\mathcal{R}}(T,x_{0})\subseteq\mathcal{R}^{\mathcal{G}}(T,x_{0})$.

###### Proof.

Proposition 1 and Theorem 1 show that
$\bar{\mathcal{R}}(T,x_{0})\subseteq\mathcal{R}^{\mathcal{G}}(T,x_{0})$. ∎

We can expand on Theorem 3 using control system (26) to help determine the
geometric structure of $\bar{\mathcal{R}}(T,x_{0})$.

Corollary 3: Consider a control system of the form $\dot{x}=(b-c\|x\|)u$
defined on some ball $\mathbb{B}\subseteq\mathbb{R}^{n}$, such that
$x\in\mathbb{R}^{n}$, $b,c\in\mathbb{R}$, and $u\in\mathcal{U}$. Then, the
reachable set $\bar{\mathcal{R}}(T,x_{0})$ is a ball in $\mathbb{B}$.

###### Proof.

Let $z=Rx$ such that $R\in\mathbb{R}^{n\times n}$ is any orthonormal matrix.
Obviously, $\dot{z}=R\dot{x}=(b-c\|x\|)Ru$. Note that $\|Rx\|=\|x\|=\|z\|$,
thus we have $\dot{z}=(b+c\|z\|)Ru$. Similarly, let $Ru=v$. Then,
$\|v\|=\|Ru\|=\|u\|$. Therefore, $\dot{z}=(b-c\|z\|)v$ such that
$v\in\mathcal{U}$, so the reachable set of $\dot{x}=(b-c\|x\|)u$ is invariant
to all rotations. Any rotation of any point on any trajectory of the original
system is thus itself on some other trajectory of the original system. Hence,
as the trajectories are continuous, the reachable set of the original system
is a ball. ∎

###  IV-B Underapproximated Control System – Polygon

Analogous to Theorem 3, we can define $\hat{\mathcal{R}}(T,x_{0})$ as the
reachable set of (25) at time $T$, and again
$\hat{\mathcal{R}}(T,x_{0})\subseteq\mathcal{R}^{\mathcal{G}}(T,x_{0})$. We
follow the method above and interpret inclusion (25) as a control system
defined in Theorem 4.

Theorem 4: Let $s\in\mathbb{R}$ and
$g(s)=\|G(0)^{\dagger}\|^{-1}-(L_{G}+L_{f})s$,
$\alpha(s)=\|G(0)^{\dagger}\|^{-1}-L_{G}s$,
$\beta(s)=\mu\|G(0)^{\dagger}\|L_{G}s$ with $\mu$ as defined in Theorem 2. Let
$U\Sigma V^{T}$ be the singular value decomposition of $G(0)$ where
$U=[\eta_{1},...,\eta_{n}]$. Let $r=\text{rank}(G(0))$; we define
$\Lambda(s)=\text{diag}(\lambda_{i}(s))$ such that
$\lambda_{i}(s)=\text{max}\\{\frac{g(s)}{\alpha(s)\|G(0)^{\dagger}\eta_{i}\|+\beta(s)},g(s),0\\}$
for $i=1,...,r$ and $\lambda_{i}(s)=0$ elsewhere.

The reachable set of (25) equals the reachable set of the control system

$\dot{x}=a+U\Lambda(\|x\|)u,\indent x(0)=x_{0},$ (27)

on $\\{x~{}|~{}\|x\|\leq\|G(0)^{\dagger}\|/(L_{f}+L_{G})\\}$, with $a=f(0)$
and $u\in\\{u~{}|~{}\|u\|_{1}\leq 1\\}\cap\text{Im}(G(0)^{T})$. If
$\hat{\mathcal{R}}(T,x_{0})$ denotes the reachable set of (27), then
$\hat{\mathcal{R}}(T,x_{0})\subseteq\mathcal{R}^{\mathcal{G}}(T,x_{0})$.

###### Proof.

Let $\|x\|\leq\|G(0)^{\dagger}\|/(L_{f}+L_{G})\\}$. By Theorem 1, the boundary
of $\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ lies at a magnitude of $g(\|x\|)$
along any direction in $\text{Im}(G(0))$ from $f(0)$. By Theorem 2, the
boundary of $\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ lies at a magnitude of
$g(\|x\|)/(\alpha(\|x\|)\|G(0)^{\dagger}\eta_{i}\|+\beta(\|x\|))$ along the
direction $\eta_{i}$ from $f(0)$, for all $i\leq r$. Since
$\hat{\mathcal{V}}^{\mathcal{G}}_{x}=\bar{\mathcal{V}}^{\mathcal{G}}_{x}\cup\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$,
the vertex of $\mathcal{S}(x)$ in the direction $\eta_{i}$ thus has a
magnitude of $\lambda_{i}(\|x\|)$, for $i\leq r$. In the direction of
$\eta_{i}$ for $i>r$, the distance of boundary points for both
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ and
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ from $f(0)$ equals $0$, as both of
these sets are contained in Im(G(0)). By its construction in Lemma 2, polygon
$P(\mathcal{S}(x))$ is thus given by $f(0)+U\Lambda(\|x\|)Q$, where
$Q=\\{u~{}|~{}\|u\|_{1}\leq 1\\}$. ∎

We note that the dynamics of (26) and (27) are entirely known. Finding
$\bar{\mathcal{R}}(T,x_{0})$ and $\hat{\mathcal{R}}(T,x_{0})$ becomes a
standard problem of determining the reachable set of a nonlinear control
system [28]. In order to exploit previous work on computing reachable sets in
$\mathbb{R}^{n}$, we can follow methods outlined in [1], where (24) and (25)
can be continuously extended to all $\mathbb{R}^{n}$ by defining
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}=\\{f(0)\\}$ and
$P(\mathcal{S}(x))=\\{f(0)\\}$ for all $x$ such that
$\|x\|>\frac{\|G(0)^{\dagger}\|^{-1}}{L_{f}+L_{G}}$. Although there is no
method to analytically determine the exact reachable set for all nonlinear
systems, existing level set methods find the reachable set by determining the
viscosity solution to Hamilton-Jacobi equations [11], [29]. Additional methods
create an overapproximation of the true reachable set by utilizing trajectory
piecewise linearized models [30] or set propagation techniques [31], [32],
[33]. For simplicity, in our numerical examples, we approximate the true
reachable set using a Monte Carlo method by solving ODEs (26) and (27) with
random time-varying inputs.

##  V Numerical Examples

We consider two examples. The first shows novel theory applied to an academic
three-dimensional nonlinear system with initial conditions similar to the
system briefly discussed in Section III, shown in Figure 1. This illustrates
an example where Theorem 1 results in a poorer underapproximaiton of the true
reachable set, while utilizing $\hat{\mathcal{R}}(T,x_{0})$ generates a
significantly better underapproximation. The second example is a control
system with decoupled quadrocoptor dynamics motivated by the scenario of
landing a damaged UAV safely [3], [4]. The goal is to determine a reachable
set of pitch and roll velocities in order to help stabilize the UAV for
landing.

As mentioned at the end of Section IV, we calculate the reachable sets by
numerically solving all ordinary differential equations using the standard
ode45 function in MATLAB which implements a Runge-Kutta method [34] with a
variable time step for efficient computation. Doing so avoids limitations of
current numerical solvers such as CORA [35] that often rely on set propagation
methods to calculate the reachable set and face issues as $\|x\|$ is not
differentiable at $0$. Numerical results are supported through analytical
means and theoretical results derived in Sections II, III, and IV.

###  V-A Three-Dimensional Nonlinear System

(a)

(b)

(c)

Figure 2: True reachable set (blue) with the underapproximations
$\bar{\mathcal{R}}(T,0)$ (red) and $\hat{\mathcal{R}}(T,0)$ (green)
numerically calculated for $T\in\\{0.05,0.2,0.5\\}$ seconds.

We consider a system with dynamics

$\dot{x}(t)=f(x(t))+G(x(t))u(t)=\begin{bmatrix}10u_{1}+(3-x_{2})u_{2}\\\
(2-x_{1})u_{1}+7u_{2}\\\ (2.5+x_{3})u_{3}\end{bmatrix},$ (28)

where

$f(x(t))=\begin{bmatrix}0\\\ 0\\\ 0\end{bmatrix},\indent
G(x(t))=\begin{bmatrix}10&3-x_{2}&0\\\ 2-x_{1}&7&0\\\
0&0&2.5+x_{3}\end{bmatrix},$

with the primary interest of finding the reachable set of $x_{1}$ and $x_{2}$.
Taking the Jacobian of $f(x)$ and $G(x)$ yields $L_{f}=L_{G}=1$ to be
acceptable Lipschitz constants. For simplicity, we set $x_{0}=0$. We remind
the reader that $f$ and $G$ are assumed to be unknown, and only $L_{f}$ and
$L_{G}$ are known.

Results from Figure 2, showing the projection of the reachable sets to the
first two coordinates, illustrate that (i) $\bar{\mathcal{R}}(T,0)$ and
$\hat{\mathcal{R}}(T,0)$ are indeed underapproximations of the true reachable
set, (ii) $\hat{\mathcal{R}}(T,0)$ produces a better underapproximation than
$\bar{\mathcal{R}}(T,0)$, and (iii) the accuracy of the underapproximations
increases as $T\to 0$. Phenomenon (i) validates the results from Theorems 3
and 4. The singular value decomposition of $G(0)$ yields $\sigma_{1}\approx
11.43$, $\sigma_{2}\approx 5.6$, and $\sigma_{3}=2.5$. Because of the large
difference between $\sigma_{1}$ and $\sigma_{3}$, the true reachable set is
not accurately represented by the ball $\bar{\mathcal{R}}(T,0)$. On the other
hand, in Theorem 4 we derive control system (27) using
$\hat{\mathcal{V}}^{\mathcal{G}}_{x}$, which more accurately represents the
complex shape of the GVS. Thus, the resulting set $\hat{\mathcal{R}}(T,0)$ is
larger along singular vectors pertaining to larger singular values,
contributing to result (ii).

Lastly, Figure 2 illustrates how the underapproximations
$\bar{\mathcal{R}}(T,0)$ and $\hat{\mathcal{R}}(T,0)$ become asymptotically
perfect as $T\to 0$. As $T\to\infty$, although comparatively worse, these
underapproximations yield progressively larger reachable sets. In the next
example, we apply this novel theory to decoupled quadrocopter dynamics to help
safely land a damaged UAV. In this case, we will show that, in contrast to the
above example, the polygonal approximation obtained by Theorem 4 produces
worse results than the ball approximation derived from Theorem 3.

###  V-B Decoupled Quadrocopter Dymamics

(a)

(b)

(c)

Figure 3: True reachable set (blue) with the underapproximations
$\bar{\mathcal{R}}(T,0)$ (red) and $\hat{\mathcal{R}}(T,0)$ (green)
numerically calculated for $T\in\\{0.05,0.25,1\\}$ seconds.

We show that the novel theory can be applied to a real system by example of
quadrocopter dynamics. We consider the objective of adjusting the
quadrocopter’s pitch and roll velocities. For a safe landing, a UAV ideally
needs these velocities to equal $0$ [36]. Given the physical dimensions of a
standard UAV, the yaw rate is inconsequential due to the symmetrical shape of
the quadrocopter. We consider the scenario where a UAV collides with an
obstacle, which would result in unwanted high velocity rotations. This
scenario translates to a problem of reachability: given initial conditions of
roll and pitch rates, we aim to determine if it is possible to reach $p=q=0$
without knowing the system’s dynamics after the collision.

The dynamics of a standard UAV are modeled in [37]. The model is comprised of
a solid sphere with mass $M=1\text{kg}$ and radius $R=0.1\text{m}$, which
represents the central frame; it is connected to four point masses
$m=0.1\text{kg}$, each representing one of four propellers at an equidistant
length of $l=0.5\text{m}$ away from the central sphere. The dynamics are shown
below:

$\begin{bmatrix}\dot{p}\\\ \dot{q}\\\
\dot{r}\end{bmatrix}=\begin{bmatrix}\frac{J_{y}-J_{z}}{J_{x}}qr\\\
\frac{J_{z}-J_{x}}{J_{y}}pr\\\
\frac{J_{x}-J_{y}}{J_{z}}pq\end{bmatrix}+\begin{bmatrix}\frac{1}{J_{x}}\tau_{\phi}\\\
\frac{1}{J_{y}}\tau_{\theta}\\\ \frac{1}{J_{z}}\tau_{\psi}\end{bmatrix},$ (29)

where states $p$, $q$, and $r$ correspond to velocities pertaining to the roll
rate, pitch rate, and yaw rate, respectively, and

$J_{x}=\frac{2MR^{2}}{5}+2l^{2}m,\indent J_{y}=J_{x},\indent
J_{z}=\frac{2MR^{2}}{5}+4l^{2}m.$

Inputs $\tau_{\phi}$, $\tau_{\theta}$, and $\tau_{\psi}$ pertain to the
applied torque that directly affects the roll, pitch, and yaw velocities. As
$J_{x}=J_{y}$, the yaw rate can directly be changed by increasing the control
action to $\tau_{\psi}$ without affecting the roll or pitch rates. Including
additional states, such as the roll, pitch, and yaw angles and translational
position, would move the model beyond the requirements of Assumption 1.
Applying the theory to a dynamic structure such as this is a subject for
future work.

We note the system dynamics in (29) are not globally Lipschitz continuous.
However, as previously discussed, given the yaw rate can be directly changed
by increasing control action to $\tau_{\psi}$ without affecting other states,
we can trivially reduce $\dot{r}$ to $0$, causing $r$ to be some constant; we
arbitrarily set $r_{0}=\pi/2$. Let the initial conditions after collision be
$p_{0}=15$, $q_{0}=10$ radians per second. Since novel theoretical results are
derived under the assumption $x_{0}=0$, we perform a simple coordinate
transformation; let $\bar{p}=p-15$ and $\bar{q}=q-10$ such that
$\bar{p}_{0}=\bar{q}_{0}=0$. Obviously, $\dot{\bar{p}}=\dot{p}$ and
$\dot{\bar{q}}=\dot{q}$, thus resulting in the new system

$\begin{bmatrix}\dot{\bar{p}}\\\
\dot{\bar{q}}\end{bmatrix}=\begin{bmatrix}\frac{\pi(J_{y}-J_{z})}{2J_{x}}(\bar{q}+10)\\\
\frac{\pi(J_{z}-J_{x})}{2J_{y}}(\bar{p}+15)\end{bmatrix}+\begin{bmatrix}\frac{1}{J_{x}}\tau_{\phi}\\\
\frac{1}{J_{y}}\tau_{\theta}\end{bmatrix},$ (30)

with initial conditions

$f(0)=\begin{bmatrix}\frac{10\pi(J_{y}-J_{z})}{2J_{x}}\\\
\frac{15\pi(J_{z}-J_{x})}{2J_{y}}\end{bmatrix}\approx\begin{bmatrix}-8.73\\\
13.09\end{bmatrix},\indent G(0)=\begin{bmatrix}\frac{1}{J_{x}}&0\\\
0&\frac{1}{J_{y}}\end{bmatrix},$

such that $J_{x}=J_{y}=0.009$, and $\|G(0)^{\dagger}\|^{-1}\approx 111.11$.
The new dynamics (30) are Lipschitz continuous; we overapproximate the
Lipschitz bounds to be $L_{f}=1$ and $L_{G}=1$ to account for overly
conservative knowledge about the rate of change of the dynamics.

As in the previous sections, we bound inputs $u\in\mathbb{B}^{2}(0;1)$ and
denote $x=[\bar{p},\bar{q}]^{T}\in\mathbb{R}^{2}$ as the system states. The
blue shape in Figure 3 illustrates the true reachable set of the states
pertaining to the roll and pitch velocities, in $(\bar{p},\bar{q})$
coordinates, plotted at $T\in\\{0.05,0.25,1\\}$ seconds. Given the structure
of control systems (26) and (27), it is trivial to see that if point
$(-p_{0},-q_{0})$ is contained within either $\bar{\mathcal{R}}(T,0)$ or
$\hat{\mathcal{R}}(T,0)$, then we know there exists a control $u$ which
guarantees velocities $p$ and $q$ can be reduced to $0$, regardless of the
true system dynamics.

Next, we apply novel theory to solve for $\bar{\mathcal{R}}(T,0)$, resulting
in the set shown in red in Figure 3, which is proven in Sections III and IV to
be a guaranteed underapproximation of the GRS. We begin by calculating
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}=\mathbb{B}^{2}(f(0);\|G(0)^{\dagger}\|^{-1}-(L_{f}+L_{G})\|x\|)=\mathbb{B}^{2}(f(0);111.11-2\|x\|)$.
Control system (26) thus equals

$\dot{x}=\begin{bmatrix}-8.73\\\
13.09\end{bmatrix}+\begin{bmatrix}111.11-2\|x\|&0\\\ 0&111.11-2\|x\|\\\
\end{bmatrix}u$ (31)

such that $u\in\mathbb{B}^{2}(0;1)$. Its reachable set is
$\bar{\mathcal{R}}(T,0)$.

Lastly, we consider control system (27) to determine $\hat{\mathcal{R}}(T,0)$.
Since $G(0)$ is invertible, we know from Corollary 2 that
$\bar{\mathcal{V}}^{\mathcal{G}}_{x}\subseteq\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$.
Thus, when determining the diagonal entries of $\Lambda(\|x\|)$, we know from
Theorem 4 that
$\lambda_{i}(\|x\|)=\frac{g(\|x\|)}{\alpha(\|x\|)\|G(0)^{\dagger}\eta_{i}\|+\beta(\|x\|)}$.
Since $G(0)$ is diagonal, we know the matrix of left singular vectors $U$
equals identity. System (27) therefore equals

$\dot{x}=\begin{bmatrix}-8.73\\\
13.09\end{bmatrix}+\begin{bmatrix}\lambda_{1}(\|x\|)&0\\\
0&\lambda_{2}(\|x\|)\end{bmatrix}u$

such that $\|u\|_{1}\leq 1$.

The same methods for determining the previous two reachable sets numerically
can be applied to find $\hat{\mathcal{R}}(T,0)$. Figure 3 displays
$\hat{\mathcal{R}}(T,0)$ for the roll and pitch velocities $\bar{p}$ and
$\bar{q}$. Notice that $\|G(0)^{\dagger}\eta_{i}\|$ is identical for all $i$.
Thus, $\lambda_{i}(\|x\|)=111.11-2\|x\|$ for all $i$, which is identical to
the diagonal terms in control system (31). Since system (27) considers
$\|u\|_{1}\leq 1$ and system (26) considers $u\in\mathbb{B}^{2}(0;1)$, for
this particular control system,
$\hat{\mathcal{R}}(T,0)\subseteq\bar{\mathcal{R}}(T,0)$. The final result is a
reachable set $\hat{\mathcal{R}}(T,0)$ denoted in green in Figure 3. It
naturally resembles a polyhedron as $T\to 0$, with edges that curve as time
increases.

According to the numerical results in Figure 3, roll and pitch velocities with
initial conditions $p_{0}=15$, $q_{0}=10$ radians per second can be reduced to
$0$ in no more than $0.25$ seconds because $(-15,-10)$ is contained in
$\bar{\mathcal{R}}(0.25,0)$. In fact, for any points that lie within
$\bar{\mathcal{R}}(T,0)$, we can guarantee there exists a control signal which
reaches these states within time $T$. For time $T=0.05$, we see that neither
the true reachable set, nor the underapproximations reach $(-15,-10)$, with
$\bar{\mathcal{R}}(T,0)$ becoming asymptotically perfect as $T\to 0$.
Conversely, at time $T=1$, the theory provides larger, generally worse
underapproximations where the true reachable set, along with
$\bar{\mathcal{R}}(T,0)$ and $\hat{\mathcal{R}}(T,0)$, clearly include
$(-15,-10)$.

Unlike in the first numerical example, set $\bar{\mathcal{R}}(T,0)$ is larger
than $\hat{\mathcal{R}}(T,0)$. Naturally, in general, computing
$\bar{\mathcal{R}}(T,0)\cup\hat{\mathcal{R}}(T,0)$ provides the best
approximation of $\mathcal{R}^{\mathcal{G}}(T,0)$ available from our theory.

##  VI Conclusion and Future Work

This paper provides a novel approach to underapproximating the reachable set
of a system with unknown dynamics. By assuming the nonlinear control-affine
system structure and exploiting solely the knowledge of system dynamics at a
single point and — possibly conservative — Lipschitz bounds on the rate of
change, we are able to determine two underapproximations
$\bar{\mathcal{R}}(T,x_{0})$ and $\hat{\mathcal{R}}(T,x_{0})$ that are
guaranteed to be contained within the guaranteed reachable set
$\mathcal{R}^{\mathcal{G}}(T,x_{0})$. Both underapproximations rely on an
intermediate approximation of the GRS by an ODI
$\dot{x}\in\mathcal{V}^{\mathcal{G}}_{x}$, where the right hand side is a set
of guaranteed velocities for the unknown nonlinear control-affine system. The
two underapproximations differ by the shape of the right hand side set, i.e.,
the approximation of $\mathcal{V}^{\mathcal{G}}_{x}$ is either determined by
balls $\bar{\mathcal{V}}^{\mathcal{G}}_{x}$ or a different shape
$\bar{\bar{\mathcal{V}}}^{\mathcal{G}}_{x}$ that more closely resembles the
shape of an intersection of an infinite many ellipsoids.

A natural area of future work is to focus on creating a larger set
underapproximating the GVS. One possibility is to consider approximations in
norms other than the spectral norm considered in this paper. For example,
potentially by bounding the perturbations of the unknown system’s dynamics
with the Frobenius norm instead of the spectral norm, we could utilize the
Mirsky Inequality [23] to produce new underapproximations of the GVS. However,
there is currently no guarantee that such a bound would produce a more
accurate underapproximation of the GVS. Another possibility includes the
utilization of semi-infinite programming [38, 39] to determine the maximal
extent of $\mathcal{V}^{\mathcal{G}}_{x}$ in every direction. Utilizing first-
and second-order optimality conditions, the problem could potentially be
reduced locally to one with finitely many constraints; similarly, given the
convex structure of the GRS, additional optimization techniques such as
duality may help further simplify the problem.

Another approach to obtaining a larger underapproximation of the GVS is to
increase the knowledge of the system dynamics. In other words, making
additional assumptions on the structure of the dynamics could also help
determine a larger underapproximation of the GVS by reducing the size of the
set of systems consistent with prior knowledge about system dynamics. One
option is to utilize knowledge from multiple system runs instead of solely
dynamics at a single point. Expanding on theoretical results derived from this
paper by incorporated additional knowledge consistent with a large class of
systems could potentially result in substantial progress for the development
of sophisticated safety critical systems.

## References

  * [1] M. Ornik, “Guaranteed reachability for systems with unknown dynamics,” in _59th IEEE Conference on Decision and Control_ , 2020, pp. 2756–2761.
  * [2] N. Nguyen, K. Krishnakumar, J. Kaneshige, and P. Nespeca, “Flight dynamics and hybrid adaptive control of damaged aircraft,” _Journal of Guidance, Control, and Dynamics_ , vol. 31, no. 3, pp. 751–764, 2008.
  * [3] G. Chowdhary, E. N. Johnson, R. Chandramohan, M. S. Kimbrell, and A. Calise, “Guidance and control of airplanes under actuator failures and severe structural damage,” _Journal of Guidance, Control, and Dynamics_ , vol. 36, no. 4, pp. 1093–1104, 2013.
  * [4] D. Jourdan, M. Piedmonte, V. Gavrilets, D. Vos, and J. McCormick, “Enhancing uav survivability through damage tolerant control,” in _AIAA Guidance, Navigation, and Control Conference_ , 2010, pp. 7548–7569.
  * [5] R. W. Brockett, “Nonlinear systems and differential geometry,” _Proceedings of the IEEE_ , vol. 64, no. 1, pp. 61–72, 1976.
  * [6] A. Isidori, _Nonlinear Control Systems_. Springer Science & Business Media, 2013.
  * [7] M. Ornik, S. Carr, A. Israel, and U. Topcu, “Control-oriented learning on the fly,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 11, pp. 4800–4807, 2019.
  * [8] I. M. Mitchell and C. J. Tomlin, “Overapproximating reachable sets by hamilton-jacobi projections,” _Journal of Scientific Computing_ , vol. 19, no. 1, pp. 323–346, 2003.
  * [9] T. F. Filippova, “Estimates of reachable sets of a nonlinear dynamical system with impulsive vector control and uncertainty,” in _14th International Conference: Stability and Oscillations of Nonlinear Control Systems_. IEEE, 2018, pp. 1–4.
  * [10] M. Rungger and M. Zamani, “Accurate reachability analysis of uncertain nonlinear systems,” in _21st International Conference on Hybrid Systems: Computation and Control_ , 2018, pp. 61–70.
  * [11] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games,” _IEEE Transactions on Automatic Control_ , vol. 50, no. 7, pp. 947–957, 2005\.
  * [12] B. Zhang, J. Lam, and S. Xu, “Reachable set estimation and controller design for distributed delay systems with bounded disturbances,” _Journal of the Franklin Institute_ , vol. 351, no. 6, pp. 3068–3088, 2014.
  * [13] G. E. Dullerud and F. Paganini, _A Course in Robust Control Theory: A Convex Approach_. Springer Science & Business Media, 2013.
  * [14] P. A. Ioannou and J. Sun, _Robust Adaptive Control_. Courier Corporation, 2012.
  * [15] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016.
  * [16] Y. Chen, H. Peng, J. Grizzle, and N. Ozay, “Data-driven computation of minimal robust control invariant set,” in _IEEE Conference on Decision and Control_. IEEE, 2018, pp. 4052–4058.
  * [17] J.-P. Aubin and A. Cellina, _Differential Inclusions: Set-Valued Maps and Viability Theory_. Springer Science & Business Media, 2012.
  * [18] A. Bressan and B. Piccoli, _Introduction to the Mathematical Theory of Control_ , 2007.
  * [19] G. V. Smirnov, _Introduction to the Theory of Differential Inclusions_. American Mathematical Society, 2002.
  * [20] A. B. Kurzhanski and P. Varaiya, “Ellipsoidal techniques for reachability analysis,” in _International Workshop on Hybrid Systems: Computation and Control_. Springer, 2000, pp. 202–214.
  * [21] M. Margaliot, “On the reachable set of nonlinear control systems with a nilpotent lie algebra,” in _European Control Conference_. IEEE, 2007, pp. 4261–4267.
  * [22] R. Vinter, “A characterization of the reachable set for nonlinear control systems,” _SIAM Journal on Control and Optimization_ , vol. 18, no. 6, pp. 599–610, 1980.
  * [23] G. W. Stewart, “Perturbation theory for the singular value decomposition,” Tech. Rep., 1998.
  * [24] G. Strang, _Introduction to Linear Algebra_. Wellesley-Cambridge Press, 2016.
  * [25] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex Optimization_. Cambridge University Press, 2004.
  * [26] K. A. Ross, _Elementary Analysis_. Springer, 2013.
  * [27] G. W. Stewart, “On the perturbation of pseudo-inverses, projections and linear least squares problems,” _SIAM Review_ , vol. 19, no. 4, pp. 634–662, 1977\.
  * [28] E. D. Sontag, “Controllability is harder to decide than accessibility,” _SIAM Journal on Control and Optimization_ , vol. 26, no. 5, pp. 1106–1118, 1988.
  * [29] M. Chen, _High Dimensional Reachability Analysis: Addressing the Curse of Dimensionality in Formal Verification_. University of California, Berkeley, 2017.
  * [30] Z. Han and B. H. Krogh, “Reachability analysis of nonlinear systems using trajectory piecewise linearized models,” in _American Control Conference_. IEEE, 2006, pp. 1505–1510.
  * [31] M. Althoff, O. Stursberg, and M. Buss, “Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization,” in _47th IEEE Conference on Decision and Control_. IEEE, 2008, pp. 4042–4048.
  * [32] B. Chachuat, B. Houska, R. Paulen, N. Perić, J. Rajyaguru, and M. E. Villanueva, “Set-theoretic approaches in analysis, estimation and control of nonlinear systems,” _9th IFAC Symposium on Advanced Control of Chemical Processes_ , vol. 48, no. 8, pp. 981–995, 2015.
  * [33] N. Ramdani and N. S. Nedialkov, “Computing reachable sets for uncertain nonlinear hybrid systems using interval constraint-propagation techniques,” _Nonlinear Analysis: Hybrid Systems_ , vol. 5, no. 2, pp. 149–162, 2011.
  * [34] J. R. Dormand and P. J. Prince, “A family of embedded runge-kutta formulae,” _Journal of Computational and Applied Mathematics_ , vol. 6, no. 1, pp. 19–26, 1980.
  * [35] M. Althoff, D. Grebenyuk, and N. Kochdumper, “Implementation of taylor models in CORA 2018,” in _5th International Workshop on Applied Verification for Continuous and Hybrid Systems_ , 2018.
  * [36] K. A. Ghamry, Y. Dong, M. A. Kamel, and Y. Zhang, “Real-time autonomous take-off, tracking and landing of uav on a moving ugv platform,” in _24th Mediterranean Conference on Control and Automation_. IEEE, 2016, pp. 1236–1241.
  * [37] R. W. Beard, “Quadrotor dynamics and control,” 2008.
  * [38] R. Hettich and K. O. Kortanek, “Semi-infinite programming: Theory, methods, and applications,” _SIAM Review_ , vol. 35, no. 3, pp. 380–429, 1993.
  * [39] M. Lopez and G. Still, “Semi-infinite programming,” _European Journal of Operational Research_ , vol. 180, no. 2, pp. 491–518, 2007.


R5G/hBN for $V=-56\,\text{meV}$ in the average interaction scheme. Since $V$
is negative, the lowest conduction bands are polarized towards the hBN layer
and hence strongly experience the moiré hBN coupling at the non-interacting
level. As a result, the pseudophonon energies are significantly larger, and
approach $\simeq 10\,\text{meV}$ for $\kappa_{\text{hBN}}=1$.

In Fig. S37, we compare the collective modes of the $|C|=0$ and $|C|=1$
insulators for $\xi=1$ R5G/hBN (our HF calculations can converge to either
state for the same Hamiltonian). Both insulators have pseudophonon modes that
develop comparable gaps. However, the $|C|=0$ state has a significantly
smaller valley gap. The large difference in the valley gap of $|C|=0$ and
$|C|=1$ insulators has previously been pointed out in twisted TMD bilayers
[56], where the difference was attributed to the underlying electronic
topology in the HF state.

In Fig. S38, we show representative results using the single-particle model of
Ref. [26]. Note that following Ref. [26], we use a 2d interaction and the CN
interaction scheme, and project our calculations using the bare band basis.
The pseudophonon gaps are smaller due to the use of the CN scheme, and the
contrast of the valley gap between $|C|=0$ and $|C|=1$ persists here as well.
In Fig. S39, we further compare the results between the CN scheme and the
$(3+3)$-average scheme (see Sec. B.3). As expected, the pseudophonons develop
a larger gap in the latter calculation. This suggests that the enhancement of
the pseudophonon gap in the average interaction scheme relative to the CN
interaction scheme is not unique to the single-particle model constructed in
Ref. [28].

Figure S35: $\mathbf{q}=0$ collective modes of the R5G/hBN $|C|=1$ Chern
insulator. Top (bottom) row shows results for the average (CN) scheme. Left
shows the $\mathbf{q}=0$ collective modes, filtered by whether they carry
valley or spin charge, as a function of the strength of hBN coupling, where
$\kappa_{\text{hBN}}=1$ corresponds to the physical model and the
$\kappa_{\text{hBN}}=0$ limit has continuous translation symmetry. Right shows
the HF band structure for $\kappa_{\text{hBN}}=0,1$ in the spin sector with
finite filling (the other spin sector is fully unoccupied). Parameters have to
chosen so that the interacting gap is similar in the calculations. The HF and
TDHF calculations are performed with $(0+4)$ screened basis projection,
$\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times 12$. Figure S36:
$\mathbf{q}=0$ collective modes of the R5G/hBN $|C|=1$ Chern insulator. Left
shows the $\mathbf{q}=0$ collective modes, filtered by whether they carry
valley or spin charge, as a function of the strength of hBN coupling, where
$\kappa_{\text{hBN}}=1$ corresponds to the physical model and the
$\kappa_{\text{hBN}}=0$ limit has continuous translation symmetry. Right shows
the HF band structure for $\kappa_{\text{hBN}}=0,1$ in the spin sector with
finite filling (the other spin sector is fully unoccupied). The HF and TDHF
calculations are performed with $(0+4)$ screened basis projection, average
interaction scheme, $\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times
12$. Figure S37: $\mathbf{q}=0$ collective modes of the $|C|=0$ and $|C|=1$
states in R5G/hBN. Both HF states are obtained as converged solutions in our
HF calculations. Top (bottom) row shows results for $|C|=0$ ($|C|=1$). Left
shows the $\mathbf{q}=0$ collective modes, filtered by whether they carry
valley or spin charge, as a function of the strength of hBN coupling, where
$\kappa_{\text{hBN}}=1$ corresponds to the physical model and the
$\kappa_{\text{hBN}}=0$ limit has continuous translation symmetry. Right shows
the HF band structure for $\kappa_{\text{hBN}}=0,1$ in the spin sector with
finite filling (the other spin sector is fully unoccupied). The HF and TDHF
calculations are performed with $(0+4)$ screened basis projection, average
interaction scheme, $\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times
12$. Figure S38: $\mathbf{q}=0$ collective modes of the $|C|=0$ and $|C|=1$
states in R5G/hBN using the single-particle model of Ref. [26] and the CN
interaction scheme with bare basis projection. Both HF states are obtained as
converged solutions in our HF calculations. Top (bottom) row shows results for
$|C|=0$ ($|C|=1$). Left shows the $\mathbf{q}=0$ collective modes, filtered by
whether they carry valley or spin charge, as a function of the strength of hBN
coupling, where $\kappa_{\text{hBN}}=1$ corresponds to the physical model and
the $\kappa_{\text{hBN}}=0$ limit has continuous translation symmetry. Right
shows the HF band structure for $\kappa_{\text{hBN}}=0,1$ in the spin sector
with finite filling (the other spin sector is fully unoccupied). The HF and
TDHF calculations are performed with $(0+4)$ bands, $\theta=0.60^{\circ}$, a
2d interaction $V_{\text{2d}}(\mathbf{q})$ with $d_{\text{sc}}=25\,\text{nm}$,
and system size $12\times 12$. Figure S39: $\mathbf{q}=0$ collective modes of
the $|C|=1$ state in R5G/hBN using the single-particle model of Ref. [26] with
bare basis projection. Top row shows results for the $(3+3)$-average, while
bottom row shows results for the CN scheme. Left shows the $\mathbf{q}=0$
collective modes, filtered by whether they carry valley or spin charge, as a
function of the strength of hBN coupling, where $\kappa_{\text{hBN}}=1$
corresponds to the physical model and the $\kappa_{\text{hBN}}=0$ limit has
continuous translation symmetry. Right shows the HF band structure for
$\kappa_{\text{hBN}}=0,1$ in the spin sector with finite filling (the other
spin sector is fully unoccupied). The HF and TDHF calculations are performed
with $(0+4)$ bands, $\theta=0.77^{\circ}$, a 2d interaction
$V_{\text{2d}}(\mathbf{q})$ with $d_{\text{sc}}=25\,\text{nm}$, and system
size $12\times 12$.

###  E.2 Collective mode phase diagrams

In this subsection, we repeat the TDHF calculation for a range of $V$ and
interaction strengths $\epsilon_{r}$. For each parameter, we consider the
energies $\Omega_{\text{phonon,1}}$ and $\Omega_{\text{phonon,2}}$ of the
lowest two intraflavor modes (pseudophonons), and the energy
$\Omega_{\text{valley}}$ of the lowest intervalley mode. We also show ratios
of these energy scales. $\Omega_{\text{phonon},1}/\Omega_{\text{valley},1}>1$
suggests that the stability of the state is more tied to stability against
valley fluctuations rather than translational fluctuations.

In Figs S40 and S41, we show results for $\xi=1$ R5G/hBN for the $|C|=0$ and
$|C|=1$ insulator respectively. The valley gap of the $|C|=0$ state is
significantly lower than that of the $|C|=1$ state throughout the entire phase
diagram (a similar phenomenon was pointed out for twisted TMD homobilayers in
Ref. [56], which attributed the difference to the electronic topology of the
HF state). The pseudophonon gap is smaller for larger $V$, which conforms with
the expectation that a higher displacement field drives the conduction
electrons further away from the hBN, and reduces the moiré charge density
variation in the filled valence bands.

In Figs S42 and S43, we show results for the $|C|=1$ state in R4G/hBN and
R6G/hBN respectively (both for $\xi=1$). The reduced pseudophonon gaps for
R6G/hBN is consistent with the intuition that the conduction bands feel the
moiré potential of the valence bands less strongly due to the increased
vertical spacing between them.

Figure S40: $\mathbf{q}=0$ collective modes of the $|C|=0$ insulator in
R5G/hBN as a function of $V$ and $\epsilon_{r}$.
$\Omega_{\text{phonon,1}},\Omega_{\text{phonon,2}}$ are the energies of the
lowest excitations in the intraflavor channel, and are adiabatically connected
to the two gapless Goldstone modes in the limit of zero moiré potential.
$\Omega_{\text{valley}}$ is the energy of the lowest intervalley excitation.
Grey regions denote parameters where there are no gapped HF states. The HF and
TDHF calculations are performed with $(0+4)$ screened basis projection,
$\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times 12$.  Figure S41:
$\mathbf{q}=0$ collective modes of the $|C|=1$ insulator in R5G/hBN as a
function of $V$ and $\epsilon_{r}$.
$\Omega_{\text{phonon,1}},\Omega_{\text{phonon,2}}$ are the energies of the
lowest excitations in the intraflavor channel, and are adiabatically connected
to the two gapless Goldstone modes in the limit of zero moiré potential.
$\Omega_{\text{valley}}$ is the energy of the lowest intervalley excitation.
Grey regions denote parameters where there are no gapped HF states. The HF and
TDHF calculations are performed with $(0+4)$ screened basis projection,
$\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times 12$.  Figure S42:
$\mathbf{q}=0$ collective modes of the $|C|=1$ insulator in R4G/hBN as a
function of $V$ and $\epsilon_{r}$.
$\Omega_{\text{phonon,1}},\Omega_{\text{phonon,2}}$ are the energies of the
lowest excitations in the intraflavor channel, and are adiabatically connected
to the two gapless Goldstone modes in the limit of zero moiré potential.
$\Omega_{\text{valley}}$ is the energy of the lowest intervalley excitation.
Grey regions denote parameters where there are no gapped HF states. The HF and
TDHF calculations are performed with $(0+4)$ screened basis projection,
$\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times 12$.  Figure S43:
$\mathbf{q}=0$ collective modes of the $|C|=1$ insulator in R6G/hBN as a
function of $V$ and $\epsilon_{r}$.
$\Omega_{\text{phonon,1}},\Omega_{\text{phonon,2}}$ are the energies of the
lowest excitations in the intraflavor channel, and are adiabatically connected
to the two gapless Goldstone modes in the limit of zero moiré potential.
$\Omega_{\text{valley}}$ is the energy of the lowest intervalley excitation.
Grey regions denote parameters where there are no gapped HF states. The HF and
TDHF calculations are performed with $(0+4)$ screened basis projection,
$\xi=1$, $\theta=0.77^{\circ}$, and system size $12\times 12$.

###  E.3 Collective mode wavefunctions and dispersion

In Fig. S44, we show $\mathbf{q}=0$ collective mode wavefunctions for the
$|C|=1$ insulator in $\xi=1$ R5G/hBN for hBN coupling factor
$\kappa_{\text{hBN}}=0,1$. The parameters are identical to the top row of Fig.
S35. In particular, we focus on the lowest four intraflavor modes [two
pseudophonons and two excitons] and the lowest intervalley mode. We plot the
probability density of $X$ (ph) and $Y$ (hp) components of the mode operator
(Eq. 149) across the mBZ. The (pseudo)phonons are predominantly composed of ph
and hp operators at the mBZ boundary. This is consistent with the non-
interacting band structure (Fig. 2) where the lowest conduction bands have
near-degeneracies near the
$\tilde{K}_{M},\tilde{K}^{\prime}_{M},\tilde{M}_{M},\tilde{M}^{\prime}_{M}$,
while the lowest conduction band has a large gap to higher bands at
$\tilde{\Gamma}_{M}$. (Note that the pseudophonon wavefunctions for
$\kappa_{\text{hBN}}=0$ appear to break $C_{3}$-symmetry. This is an artifact
of the fact that the pseudophonons are degenerate at $\kappa_{\text{hBN}}=0$,
so the diagonalization routine yields arbitrary superpositions of the
degenerate eigenfunctions.) The two lowest intraflavor excitons are either
localized at $\tilde{K}_{M}$, or have a node there. Finally the lowest
intervalley mode wavefunction, which has only ph components for our
calculation (see App. D), is nearly uniformly spread out over the mBZ.

In Fig. S45, we plot the dispersion of the lowest intraflavor, intervalley,
and interspin collective modes for the $|C|=1$ insulator in $\xi=1$ R5G/hBN
for hBN coupling factor $\kappa_{\text{hBN}}=0,1$. The parameters are again
identical to the top row of Fig. S35. The lowest interspin mode is a magnon
branch with quadratic dispersion at $\mathbf{q}=0$. The intervalley mode is
gapped, consistent with lack of $SU(2)$ symmetry in valley space, and has a
relativelly narrow bandwidth $\simeq 5\,\text{meV}$. Interestingly, there are
regions in momentum space where the lowest intraflavor mode has negative or
complex eigenvalues. This signals a local instability within the variational
Slater determinant manifold. From the $\kappa_{\text{hBN}}=1$ result, it
appears that the instability corresponds to a wavevector around the
$\tilde{M}_{M}$ and $\tilde{M}^{\prime}_{M}$. We have performed HF
calculations allowing for doubling of the moiré unit cell along one axis, and
indeed find spin-valley polarized HF states with lower energy than the best
translation-invariant solution. This occurs for both the average and CN
interaction schemes. We leave a detailed analysis of such moiré translation-
symmetry broken states to future work.

Figure S44: Momentum-resolved particle-hole (ph) and hole-particle (hp)
content of the $\mathbf{q}=0$ collective modes of the $|C|=1$ insulator in
R5G/hBN for hBN-coupling strength $\kappa_{\text{hBN}}=0$ (top two rows) and
$\kappa_{\text{hBN}}=1$ (bottom two rows). The data here use the same system
parameters as the top row of Fig. S35, i.e. the HF and TDHF calculations are
performed with $(0+4)$ screened basis projection, $\xi=1$,
$\theta=0.77^{\circ}$, $V=64\,\text{meV}$, $\epsilon_{r}=5$, and system size
$12\times 12$. We show results for the lowest four intraflavor modes and the
lowest interavalley mode. $N|X(\mathbf{k})|^{2}$ indicates the total
probability amplitude of $\mathbf{q}=0$ ph operators at momentum $\mathbf{k}$
in the collective mode operator Eq. 149, and has been normalized so that
uniform probability across the mBZ corresponds to a value of 1. Similarly,
$N|Y(\mathbf{k})|^{2}$ indicates the total probability of $\mathbf{q}=0$ hp
operators at momentum $\mathbf{k}$. Note that the intervalley modes do not
have hp content. The mBZ is indicated by the black hexagon, whose center
corresponds to $\tilde{\Gamma}_{M}$. Figure S45: Collective mode dispersions
$\Omega(\mathbf{q})$ of the $|C|=1$ insulator in R5G/hBN for hBN-coupling
strength $\kappa_{\text{hBN}}=0$ (top row) and $\kappa_{\text{hBN}}=1$ (bottom
row). The data here use the same system parameters as the top row of Fig. S35,
i.e. the HF and TDHF calculations are performed with $(0+4)$ screened basis
projection, $\xi=1$, $\theta=0.77^{\circ}$, $V=64\,\text{meV}$,
$\epsilon_{r}=5$, and system size $12\times 12$. We show results for the
lowest energy intraflavor, intervalley and interspin mode. White regions
indicate where the collective mode energies are negative or complex,
signalling an instability. The mBZ is indicated by the black hexagon, whose
center corresponds to $\tilde{\Gamma}_{M}$.


# Turán problems for linear forests and cliques 111This work was supported by
the National Nature Science Foundation of China (Nos.11871040, 12271337)

Tao Fang<EMAIL_ADDRESS>Xiying Yuan  Department of Mathematics, Shanghai
University, Shanghai 200444, P.R. China

###### Abstract

Given a graph $T$ and a family of graphs $\mathcal{H}$. The generalized Turán
number of $\mathcal{H}$ is the maximum number of copies of $T$ in an
$\mathcal{H}$-free graph on $n$ vertices, denoted by $ex(n,T,\mathcal{H})$.
Let $ex(n,T,\mathcal{H})$ denote the maximum number of copies of $T$ in an
$n$-vertex $\mathcal{H}$-free graph. Recently, Alon and Frankl
(arXiv2210.15076) determined the exact values of
$\rm{ex}(n,\\{K_{r+1},M_{s+1}\\})$, where $K_{r+1}$ and $M_{s+1}$ are complete
graph on $r+1$ vertices and matching of size $s+1$, respectively. Ma and Hou
(arXiv2301.05625) gave the generalized version of Alon and Frankl’s Theorem,
which determine the exact values of $ex(n,K_{r},\\{K_{k+1},M_{s+1}\\})$. Zhang
determined the exact values of $ex(n,K_{r},\mathcal{L}_{n,s})$, where
$\mathcal{L}_{n,s}$ be the family of all linear forests of order $n$ with $s$
edges. Inspired by the work of Zhang and Ma, in this paper, we determined the
exact number of $ex(n,\\{K_{r+1},\mathcal{L}_{n,s}\\})$.

###### keywords:

Turán number, Linear forest, Clique

###### MSC:

[2010] 05C05 05C35

##  1 Introduction

Given a fixed graph $T$ and a family of graphs $\mathcal{H}$. A graph $G$ is
called $\mathcal{H}$-free if for any $H\in\mathcal{H}$, $G$ contains no copy
of $H$ as subgraph. Write $N(G,T)$ for the number of copies of $T$ in $G$. The
generalized Turán number of $\mathcal{H}$ is the maximum number of copies of
$T$ in an $\mathcal{H}$-free graph on $n$ vertices and it can also be
expressed by the following formula.

$\rm{ex}(n,T,\mathcal{H})=\mbox{max $\\{N(G,T):G$ is an $n$-vertex
$\mathcal{H}$-free graph\\}}.$

For $T=K_{2}$, it reduces to the classical Turán number, that is
$ex(n,\mathcal{H})=ex(n,K_{2},\mathcal{H})$.

A matching in a graph $G$ is a set of non-loop edges with no shared endpoints.
Denote by $M_{k}$ a matching containing $k$ edges. We denote by $\nu(G)$ the
number of edges in a maximum matching of $G$. A linear forest is graph
consisting of vertex disjoint paths or isolated vertices. Denote by
$\mathcal{L}_{n,s}$ the set of all linear forests of order $n$ with exactly
$s$ edges. Let $K_{r}$ denote a complete graph on $r$ vertices. Denote by
$E_{k}$ a set of $k$ vertices in the graph that are not adjacent to each
other. The join of two disjoint graphs $H_{1}$ and $H_{2}$, denoted by
$H_{1}\vee H_{2}$, is the graph whose vertex set is $V(H_{1}\vee
H_{2})=V(H_{1})\cup V(H_{2})$ and edge set is $E(H_{1}\vee H_{2})=E(H_{1})\cup
E(H_{2})\cup\\{xy:x\in V(H_{1}),y\in V(H_{2})\\}$.

There are many results about Turán number and generalized Turán number. For
Turán number, Erdös and Gallai [3] determined the Turán number of $M_{s+1}$.
In [6] , Ning and Wang used the closure technique and the counting technique
to determine the Turán number of $\mathcal{L}_{n,s}$. Recently, Alon and
Frankl [2] use Tutte-Berge Theorem to determine $ex(n,\\{K_{r+1},M_{s+1}\\})$.

######  Theorem 1.1 ([2]).

For $n\geqslant 2s+1$ and $r\geqslant 2$,

$ex(n,\\{K_{r+1},M_{s+1}\\})=\mbox{\rm{max}}\left\\{t(2s+1,r),t(s,r-1)+(n-s)s\right\\}.$

For generalized Turán number, Wang [8] use the shifted method to determine the
generalized Turán number of $M_{s+1}$.

######  Theorem 1.2 ([8]).

For any $r\geqslant 2$ and $n\geqslant 2s+1$,

$ex(n,K_{r},M_{s+1})=\mbox{\rm{max}}\left\\{\binom{2s+1}{r},\binom{s}{r}+(n-s)\binom{s}{r-1}\right\\}.$

In [9], Zhang et al. determined the exact values of
$ex(n,K_{r},\mathcal{L}_{n,s})$.

######  Theorem 1.3 ([9]).

For any $r\geqslant 2$ and $n\geqslant s+1$,

$ex(n,K_{r},\mathcal{L}_{n,s})=\mbox{\rm{max}}\left\\{\binom{s}{r},\binom{\left\lceil\frac{s+1}{2}\right\rceil}{r}+\left(n-\left\lceil\frac{s+1}{2}\right\rceil\right)\binom{\left\lfloor\frac{s-1}{2}\right\rfloor}{r-1}\right\\}.$

For some integers $t\geqslant k\geqslant r$, let
$\Delta_{t,k}^{r}=N(T_{k}(t),K_{r})$. Recently, Ma and Hou [5] determined the
exact values of $ex(n,K_{r},\\{K_{k+1},M_{s+1}\\})$.

######  Theorem 1.4 ([5]).

For $n\geqslant 2s+1$ and $k\geqslant r\geqslant 3$,

$ex(n,K_{r},\\{K_{k+1},M_{s+1}\\})=\mbox{\rm{max}}\\{\Delta_{2s+1,k}^{r},\Delta_{s,k-1}^{r}+(n-s)\Delta_{s,k-1}^{r-1}\\}.$

We were inspired by Zhang and Ma’s article. In this paper, we determine the
Turán number of $\\{K_{r+1},\mathcal{L}_{n,s}\\}$, which generalized the
result of Theorem 1.1.

######  Theorem 1.5.

For $n\geqslant 2s+1$ and $r\geqslant 2$,

$ex(n,\left\\{K_{r+1},\mathcal{L}_{n,s}\right\\})=\mbox{\rm{max}}\left\\{t(s,r),t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,r-1\right)+\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right)\left\lfloor\frac{s-1}{2}\right\rfloor\right\\}.$

##  2 Preliminaries

Our proof is mainly based on the strong-shifting operation and the strong-
closure technique, and we will cover their initial versions first.

Let $V(G),E(G)$ be the vertex set and edge set of the graph $G$ respectively.
For $1\leqslant i<j\leqslant n$ and $e\in E(G)$, we define a shifting
operation $S_{ij}$ on $e$ as follow:

$S_{ij}(e)=\begin{cases}(e-\\{j\\})\cup\\{i\\},&\mbox{if $j\in e$, $i\notin e$
and $(e-\\{j\\})\cup\\{i\\}\notin E(G)$},\\\ e,&\mbox{otherwise}.\end{cases}$

We define $S_{ij}(G)$ to be a shifted graph on vertex set $V(G)$ with edge set
$\\{S_{ij}(e):e\in E(G)\\}$.

In 1981, Kelmans [4] first introduced the shifting operation. In 1985, Akiyama
and Frankl [1] proved the property $\nu(S_{ij}(G))\leqslant\nu(G)$ and used it
to give a short proof of Erdös-Gallai theorem. In 2019, Wang [8] proved the
property $N(S_{ij}(G),K_{r})\geqslant N(G,K_{r})$ and obtained the exact value
of $ex(n,K_{r},M_{s+1})$. Obviously, doing shifting operation on graph $G$
will not change the number of edges of the original graph, i.e.,
$e(S_{ij}(G))=e(G)$. 1n 2022, Zhang [9] proved the following property and
obtained the exact value of $ex(n,K_{r},\mathcal{L}_{n,s})$.

######  Lemma 2.1 ([9]).

Suppose $G$ is a graph on vertex set $[n]$. For any $1\leqslant i<j\leqslant
n$, we have $S_{ij}(G)$ is also $\mathcal{L}_{n,s}$-free if $G$ is
$\mathcal{L}_{n,s}$-free.

We then introduce a variant of the shifting operation, the strong-shifting
operation. Let $V(G)$, $E(G)$ be the vertex set and edge set of the graph $G$
respectively. For any $k\in[n]$, $1\leqslant i<j\leqslant n$ and
$e=\\{j,k\\}\in E(G)$, we define a strong-shifting operation $S^{\prime}_{ij}$
on $e$ as follow:

$S^{\prime}_{ij}(e)=\begin{cases}(e-\\{j\\})\cup\\{i\\},&\mbox{if
$(e-\\{j\\})\cup\\{i\\}\notin E(G)$ and $k,i$ different color},\\\
e,&\mbox{otherwise}.\end{cases}$

Let $S^{\prime}_{ij}(G)$ be a strong-shifted graph on vertex set $V(G)$ with
edge set $\\{S^{\prime}_{ij}(e):e\in E(G)\\}$. Then
$e(S^{\prime}_{ij}(G))=e(G)$.

The following lemma is similar to Lemma 2.1, which introduces that for two
endpoints at the edges of different partitions, the graph will keep
$\\{K_{k+1},\mathcal{L}_{n,s}\\}$-free after doing the strong-shifting
operation.

######  Lemma 2.2.

Suppose $G$ is a graph on vertex set $[n]$. For any $1\leqslant i<j\leqslant
n$, we have $S^{\prime}_{ij}(G)$ is also
$\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free if $G$ is
$\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free.

###### Proof.

Let $G$ be the graph without both $K_{r+1}$ and $\mathcal{L}_{n,s}$. We apply
the strong-shifting operation $S_{ij}^{\prime}$ to $G$ for all $i,j$ for
$1\leqslant i<j\leqslant n$. Finally, we obtain a graph $S_{ij}^{\prime}(G)$.
By the definition of strong-shifting operation, the two endpoints of each new
edge belong to different partitions. Hence, $S_{ij}^{\prime}(G)$ is
$K_{r+1}$-free.

By the definition of shifting operation and strong-shifting operation, the
edges that need to be changed for strong-shifting operation are a subset of
shifting operation. It can be seen from Lemma 2.1 that $S_{ij}^{\prime}(G)$ is
$\mathcal{L}_{n,s}$-free. ∎

A complete multipartite graph is a simple graph whose vertices can be
partitioned into sets so that $\\{u,v\\}$ exists if and only if $u$ and $v$
belong to different sets of the partition. The Turán graph $T(n,r)$ is the
complete $r$-partite graph with $n$ vertices whose partite sets differ in size
by at most one. Let $t(n,r)=|T(n,r)|$ be the Turán number. The famous Turán
Theorem [7] states that $ex(n,K_{r+1})=t(n,r)$.

######  Lemma 2.3 ([7]).

Among the $n$-vertex simple graphs with no $r+1$-clique, $T(n,r)$ has the
maximum number of edges.

Next, we will introduce the closure technique and its variant, which we call
the strong-closure technique.

Suppose $G$ is a graph with order $n$, $P$ is a property that applies to $G$,
and $k$ is a positive integer. A property $P$ is considered to be $k$-stable,
if whenever $G+uv$ has property $P$ and $d_{G}(u)+d_{G}(v)\geqslant k$, then
$G$ itself also has property $P$. We can define the $k$-closure of $G$,
represented by ${cl}_{k}(G)$, as the graph $H$ formed by repeatedly connecting
non-adjacent vertices whose degree sum is at least $k$, until
$d_{H}(u)+d_{H}(v)<k$ for all $uv\notin E(H)$. It is then straightforward to
see that if $P$ is $k$-stable and ${cl}_{k}(G)$ possesses property $P$, then
$G$ also has property $P$.

In [6], Ning and Wang proof that the property $\mathcal{L}_{n,s}$-free is
$s$-stable.

######  Lemma 2.4 ([6]).

Let $G$ be a graph on $n$ vertices. Suppose that $u,v\in V(G)$ with
$d(u)+d(v)\geqslant s$. Then $G$ is $\mathcal{L}_{n,s}$-free if and only if
$G+\\{u,v\\}$ is $\mathcal{L}_{n,s}$-free.

For our problem, we have the following further lemma.

######  Lemma 2.5.

Let $G$ be a graph on $n$ vertices. Suppose that $u,v\in V(G)$ with $u,v$
belong to different partitions and $d(u)+d(v)\geqslant s$. Then $G$ is
$\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free if and only if $G+\\{u,v\\}$ is
$\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free.

###### Proof.

If $G+\\{u,v\\}$ is $\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free, then clearly $G$
is $\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free. Therefore we only need to verify
the other direction.

If $G$ is $\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free. Since the added new edge
$\\{u,v\\}$ is between different partitions, it does not increase the number
of partitions of the graph $G$, $G+\\{u,v\\}$ is $K_{r+1}$-free. And then by
Lemma 2.4, $G+\\{u,v\\}$ is $\mathcal{L}_{n,s}$-free, completing the proof. ∎

##  3 Proof of the main Theorem

###### Proof of Theorem 1.5.

For $n\geqslant 2s$, when $s$ is odd, $T(s,r)$ and $T(\frac{s-1}{2},r-1)\vee
E_{n-\frac{s-1}{2}}$ are $\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free graphs. Then
the number of edges is

$\mbox{\rm{max}}\left\\{t(s,r),t\left(\frac{s-1}{2},r-1\right)+\left(n-\frac{s-1}{2}\right)\frac{s-1}{2}\right\\}.$
(3.1)

When $s$ is even, $T(s,r)$ and $T(\frac{s}{2}-1,r-1)\vee E_{n-\frac{s}{2}+1}$
are $\left\\{K_{r+1},\mathcal{L}_{n,s}\right\\}$-free. Then the number of
edges is

$\mbox{\rm{max}}\left\\{t(s,r),t\left(\frac{s}{2}-1,r-1\right)+\left(n-\frac{s}{2}+1\right)\left(\frac{s}{2}-1\right)\right\\}.$
(3.2)

From (3.1) and (3.2), we can see

$ex(n,\left\\{K_{r+1},\mathcal{L}_{n,s}\right\\})\geqslant\mbox{\rm{max}}\left\\{t(s,r),t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,r-1\right)+\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right)\left\lfloor\frac{s-1}{2}\right\rfloor\right\\}.$

Let $G$ be an $\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free graph on vertex set $[n]$
with maximum number of edges. Then we apply the strong-shifting operation
$S_{ij}^{\prime}$ to $G$ for all $i,j$ for $1\leqslant i<j\leqslant n$.
Finally, we obtain the strong-shifted graph $G^{\prime}$. By Lemma 2.2, we
have $e(G)=e(G^{\prime})$ and $G^{\prime}$ is also
$\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free.

Let $D$ be the set of all vertices in $G^{\prime}$ with degree at least
$\left\lceil\frac{s}{2}\right\rceil$. Since $G^{\prime}$ is $K_{r+1}$-free and
Lemma 2.4, $D$ forms a complete $\ell$-partite graph, where $\ell\leqslant r$.
If all vertices with degree greater than $\left\lceil\frac{s}{2}\right\rceil$
are in the same partial set, we make $D$ the empty set. Let $A$ be the set of
vertices that forms a maximal $\ell$-partite graph that contains $D$ in
$G^{\prime}$. Since $G^{\prime}$ is the graph with the maximum number of
edges, and $A$ forms a maximal $\ell$-partite graph, $A$ forms an
$\ell$-partite Turán graph. Denote $a=|A|$. It is easy to see that $a\leqslant
s$, otherwise $G^{\prime}$ contains a linear forest with $s$ edges, which
contradicts the fact that $G^{\prime}$ is $\mathcal{L}_{n,s}$-free.

Let $B=V(G^{\prime})-A$. For any $v\in B$, since $v$ is not in $D$, we have
$d_{G^{\prime}}(v)\leqslant\left\lceil\frac{s}{2}\right\rceil-1$. Let
$B_{1}(\subseteq B)$ be the set of vertices that are not fully connected to
the vertices in $A$, and $B_{2}(\subseteq B)$ the set of vertices that are
fully connected to the vertices in $A$. Denote $b_{1}=|B_{1}|$ and
$b_{2}=|B_{2}|$. Without losing generality, we can number the vertices in $A$
as $[a]$, the vertices in $B_{1}$ as $[b_{1}+a]\backslash[a]$, and the
vertices in $B_{2}$ as $[n]\backslash[b_{1}+a]$, respectively.

If there is an edge $\\{u_{1},u_{2}\\}$ in $G^{\prime}[B_{1}]$ with
$a+1\leqslant u_{1}<u_{2}\leqslant b_{1}+a$ and there exists a vertex $u(\in
A)$ such that $\\{u,u_{1}\\}\notin G^{\prime}$. Obviously, $u_{1}$ and $u_{2}$
are not the same color as the vertices in $A$, otherwise it contradicts the
fact that $A$ is the set that forms a maximum $\ell$-partite graph. Since
$u\leqslant a<u_{2}$ and $G^{\prime}$ is the strong-shifted graph, then
$\\{u,u_{1}\\}\in G^{\prime}$, a contradiction.

If there is an edge $\\{u_{1},u_{2}\\}$ in $G^{\prime}[B_{1},B_{2}]$ with
$a+1\leqslant u_{1}\leqslant b_{1}+a$, $b_{1}+a+1\leqslant u_{2}\leqslant n$
and there exists a vertex $u(\in A)$ such that $\\{u,u_{1}\\}\notin
G^{\prime}$. We have that $u_{1}$ and $u_{2}$ are not the same color as the
vertices in $A$, otherwise it contradicts the fact that $A$ is the set that
forms a maximum $\ell$-partite graph. Since $u\leqslant a<b_{1}+a+1\leqslant
u_{2}$ and $G^{\prime}$ is the strong-shifted graph, then $\\{u,u_{1}\\}\in
G^{\prime}$, a contradiction.

If there is an edge $\\{u_{1},u_{2}\\}$ in $G^{\prime}[B_{2}]$ with
$b_{1}+a+1\leqslant u_{1}<u_{2}\leqslant n$. We will discuss into following
two cases.

Case 1  If $B_{2}$ is not an empty set. Recall that $B_{2}$ is the set of
vertices that are fully connected to the vertices in $A$ and for any $v\in B$,
we have $d_{G^{\prime}}(v)\leqslant\left\lceil\frac{s}{2}\right\rceil-1$
because of $v\notin D$, $d_{A}(v)\leqslant a$. Since for any $v^{\prime}\in
B_{2}$, $d_{G^{\prime}}(v^{\prime})\geqslant a$, we have
$a\leqslant\left\lceil\frac{s}{2}\right\rceil-1$.

The following equality depends on a trick to estimate the edges outside $A$,
which was presented in [6].

$\displaystyle e(B)+e(A,B)$ $\displaystyle=\frac{1}{2}\sum_{v\in
B}d_{B}(x)+\sum_{v\in B}d_{A}(x)$ $\displaystyle=\frac{1}{2}\sum_{v\in
B}\left(d_{G^{\prime}}(x)+d_{A}(x)\right).$

Thus,

$\displaystyle e(G^{\prime})$ $\displaystyle=e(A)+e(B)+e(A,B)$
$\displaystyle=e(A)+\frac{1}{2}\sum_{v\in
B}\left(d_{G^{\prime}}(x)+d_{A}(x)\right)$ $\displaystyle\leqslant
t(a,\ell)+\frac{1}{2}\left(\left\lceil\frac{s}{2}\right\rceil-1+a\right)(n-a).$

Let
$f(a)=t(a,\ell)+\frac{1}{2}\left(\left\lceil\frac{s}{2}\right\rceil-1+a\right)(n-a)$.
Since $n\geqslant 2s$, we have

$\displaystyle f(a+1)-f(a)$
$\displaystyle=t(a+1,\ell)-t(a,\ell)+\frac{1}{2}\left(n-2a-\left\lceil\frac{s}{2}\right\rceil\right)$
$\displaystyle>\frac{1}{2}\left(n-2a-\left\lceil\frac{s}{2}\right\rceil\right)$
$\displaystyle>0,$

which implies $f(a)$ is an increasing function of $a$. Thus we get

$\displaystyle e(G^{\prime})$ $\displaystyle\leqslant f(a)$
$\displaystyle\leqslant f\left(\left\lceil\frac{s}{2}\right\rceil-1\right)$
$\displaystyle=t\left(\left\lceil\frac{s}{2}\right\rceil-1,\ell\right)+\left(\left\lceil\frac{s}{2}\right\rceil-1\right)\left(n-\left\lceil\frac{s}{2}\right\rceil+1\right)$
$\displaystyle\leqslant
t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,r-1\right)+\left\lfloor\frac{s-1}{2}\right\rfloor\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right).$
(3.3)

Case 2  If $B_{2}$ is an empty set, then $B=B_{1}$ is an independent set in
$G^{\prime}$. If there is no edge in $G^{\prime}[A,B]$, then we have

$e(G^{\prime})=e(A)\leqslant t(a,\ell)\leqslant t(s,r).$ (3.4)

If there are some edges in $G^{\prime}[A,B]$ and $\ell\geqslant 3$, then the
largest linear forest $\mathcal{L}$ will pass all the vertex in $A$. Suppose
not, if we assume that the largest linear forest in $G^{\prime}$ does not pass
thought $v(\in A)$, then we can choose an edge $\\{u,w\\}$ in $\mathcal{L}$,
where $u,w\in A$ and $u,w$ are not in the same partition as $v$, then
$\mathcal{L}-\\{u,w\\}+\\{u,v\\}+\\{w,v\\}$ is the largest linear forest in
$G^{\prime}$, a contradiction to the fact that $\mathcal{L}$ is the largest
linear forest in $G^{\prime}$. Since $A$ forms a maximal Turán graph, we can
further assume that the largest linear forest $\mathcal{L}$ pass the vertex
$x$, where $x\in B$. Let $y$ be the vertex of $G^{\prime}[A]$ such that
$\\{x,y\\}\notin G^{\prime}$. Then $G^{\prime}+\\{x,y\\}$ is
$\\{K_{r+1},\mathcal{L}_{n,s}\\}$-free and
$e(G^{\prime}+\\{x,y\\})>e(G^{\prime})$, a contradiction to the fact that
$G^{\prime}$ has maximum number of edges.

If there are some edges in $G^{\prime}[A,B]$ and $\ell=2$. For any $v\in B$,
we want to give an upper bound on $d_{G^{\prime}}(v)$. On one hand, since $v$
is not in $D$, we have
$d_{G^{\prime}}(v)\leqslant\left\lceil\frac{s}{2}\right\rceil-1$. On the other
hand, since $B_{1}$ is the set that are not fully connected to the vertices in
$A$, we have $d_{G^{\prime}}(v)\leqslant a-1$. Consequently,
$d_{G^{\prime}}(v)\leqslant\mbox{\rm{min}}\left\\{\left\lceil\frac{s}{2}\right\rceil-1,a-1\right\\}$.

Subcase 2.1  $a\leqslant\left\lceil\frac{s}{2}\right\rceil-1$.

For any $v\in B$, it follows that $d_{G^{\prime}}(v)\leqslant a-1$. Since
$B=B_{1}$ is an independent set in $G^{\prime}$, the number of edges in
$G^{\prime}$ can be bounded as follows:

$\displaystyle e(G^{\prime})$ $\displaystyle=e(A)+e(A,B)$
$\displaystyle\leqslant t(a,2)+(n-a)(a-1).$

Let $g_{1}(a)=t(a,2)+(n-a)(a-1)$. We have

$\displaystyle g_{1}(a+1)-g_{1}(a)$ $\displaystyle=t(a+1,2)-t(a,2)+n-2a$
$\displaystyle>n-2a$ $\displaystyle>0,$

which implies that $g_{1}(a)$ is an increasing function of $a$. Then

$\displaystyle e(G^{\prime})$ $\displaystyle\leqslant g_{1}(a)$
$\displaystyle\leqslant
g_{1}\left(\left\lfloor\frac{s-1}{2}\right\rfloor\right)$
$\displaystyle=t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,2\right)+\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right)\left(\left\lfloor\frac{s-1}{2}\right\rfloor-1\right)$
$\displaystyle<t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,2\right)+\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right)\left\lfloor\frac{s-1}{2}\right\rfloor.$
(3.5)

Subcase 2.2  $a\geqslant\left\lceil\frac{s}{2}\right\rceil$.

For any $v\in B$, it follows that
$d_{G^{\prime}}(v)\leqslant\left\lceil\frac{s}{2}\right\rceil-1$. Since
$B=B_{1}$ is an independent set in $G^{\prime}$, we have

$\displaystyle e(G^{\prime})$ $\displaystyle=e(A)+e(A,B)$
$\displaystyle\leqslant
t(a,2)+(n-a)\left(\left\lceil\frac{s}{2}\right\rceil-1\right).$

Let $g_{2}(a)=t(a,2)+(n-a)\left(\left\lceil\frac{s}{2}\right\rceil-1\right)$
and
$g_{2}^{\prime}(a)=g_{2}(a+1)-g_{2}(a)=\left\lceil\frac{a}{2}\right\rceil-\left\lceil\frac{s}{2}\right\rceil+1$.
We have

$g_{2}^{\prime}(a+1)-g_{2}^{\prime}(a)=\left\lceil\frac{a+1}{2}\right\rceil-\left\lceil\frac{a}{2}\right\rceil>0.$

This implies $g_{2}(a)$ is convex on
$\left[\left\lceil\frac{s}{2}\right\rceil,s-1\right]$. Thus,

$\displaystyle e(G^{\prime})$ $\displaystyle\leqslant g_{2}(a)$
$\displaystyle\leqslant\mbox{\rm{max}}\left\\{g_{2}\left(\left\lceil\frac{s}{2}\right\rceil\right),g_{2}(s-1)\right\\}$
$\displaystyle=\mbox{\rm{max}}\left\\{t\left(\left\lceil\frac{s}{2}\right\rceil,2\right)+\left(n-\left\lceil\frac{s}{2}\right\rceil\right)\left(\left\lceil\frac{s}{2}\right\rceil-1\right),t(s-1,2)+(n-s+1)\left(\left\lceil\frac{s}{2}\right\rceil-1\right)\right\\}$
$\displaystyle<t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,2\right)+\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right)\left\lfloor\frac{s-1}{2}\right\rfloor.$
(3.6)

By (3.3), (3.4), (3.5) and (3.6), we have

$ex(n,\left\\{K_{r+1},\mathcal{L}_{n,s}\right\\})\leqslant\mbox{\rm{max}}\left\\{t(s,r),t\left(\left\lfloor\frac{s-1}{2}\right\rfloor,r-1\right)+\left(n-\left\lfloor\frac{s-1}{2}\right\rfloor\right)\left\lfloor\frac{s-1}{2}\right\rfloor\right\\}.$

The proof is complete. ∎

## References

  * [1] Akiyama J., Frankl P., On the size of graphs with complete-factors, J. Graph Theory 9 (1985) 197–201.
  * [2] Alon N., Frankl P., Turán graphs with bounded mathching number, arXiv2210.15076.
  * [3] Erdös P., Gallai T., On maximal paths and circuits of graphs, Acta Math. Acad. Sci. Hung. 10 (1959) 337–356.
  * [4] Kelmans A.K., On graphs with randomly deleted edges. Acta. Math. Acad. Sci. Hung. 37 (1981) 77–88
  * [5] Ma Y., Hou X., Generalized Turán problem with bounded matching number, arXiv2301.05625.
  * [6] Ning B., Wang J., The formula for Turán number of spanning linear forests, Discrete Math. 343(8) (2020) 111924.
  * [7] Turán P., On an extremal problem in graph theory, Matematikaiés Fizikai Lapok Hung., 48 (1941) 436–452.
  * [8] Wang J., The shifting method and generalized Turán number of matchings, Eur. J. Combin, 85 (2020) 103057.
  * [9] Zhang L.P., Wang L., Zhou J., The generalized Turán number of spanning linear forest, Graphs Combin., 38(2) (2022) 40.