# Resilient State Estimation for Nonlinear Discrete-Time Systems via Input and
State Interval Observer Synthesis

Mohammad Khajenejad, Zeyuan Jin, Thach Ngoc Dinh and Sze Zheng Yong
M. Khajenejad is with the University of California, San Diego, CA, USA. Z. Jin
is with Arizona State University, Tempe, AZ, USA. T.N. Dinh is with
Conservatoire National des Arts et Métiers (CNAM), CEDRIC-Laetitia, Paris,
France. S.Z. Yong is with Northeastern University, Boston, MA, USA. (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>s.yong@northeastern.edu).This work is partially supported by NSF grant
CNS-2313814.

###### Abstract

This paper addresses the problem of  resilient state estimation and attack
reconstruction for bounded-error nonlinear discrete-time systems with
nonlinear observations/constraints, where both sensors and actuators can be
compromised by false data injection attack signals/unknown inputs.  By
leveraging mixed-monotone decomposition of nonlinear  functions, as well as
affine parallel outer-approximation of the observation functions,  along with
introducing auxiliary states  to cancel out the effect of the attacks/unknown
inputs, our proposed observer recursively computes interval estimates that by
construction, contain the true states and unknown inputs of the system.
Moreover, we provide several semi-definite programs to synthesize observer
gains  to ensure input-to-state  stability of the proposed observer and
optimality of the design in the sense of  minimum $\mathcal{H}_{\infty}$ gain.

##  I Introduction

State estimation and unknown input reconstruction are indispensable in various
engineering applications such as aircraft tracking, fault detection, attack
detection and mitigation in cyber-physical systems (CPS) and urban
transportation [1, 2, 3].  Particularly, set-membership approaches have been
proposed for bounded-error systems to provide hard accuracy bounds, which is
especially useful for obtaining robustness guarantees for safety-critical
systems. Moreover, since attackers may be strategic in adversarial settings,
the ability to simultaneously estimate states and inputs without imposing any
assumptions on the unknown inputs /attack signals is desirable and often
crucial.

_Literature review._ Numerous studies in the literature have investigated
_secure estimation_ , i.e., how to accurately estimate the states of a system
when it is under attack or subject to adversarial signals. For instance,
secure state estimation and control problem was addressed in the presence of
false data injection attacks on both the actuators and sensors in [4], in
which a $\chi^{2}$ detector was proposed to detect malicious attacks. The
research in [5] proposed a sliding-mode observer to simultaneously estimate
system states and attacks, while the work in [6] provided a projected sliding-
mode observer-based estimation approach to reconstruct system states. Further,
the work in [7] reconstructed attack signals from the equivalent output
injection signal using a sliding-mode observer, while in [8], an attack was
considered as an auxiliary state and estimated by employing a robust switching
Luenberger observer assuming sparsity. However, all the aforementioned works
considered stochastic/Gaussian noise and hence do not apply to the bounded-
error setting we consider in this paper, where noise/disturbance signals are
assumed to be distribution-free and bounded.

A related body of literature that could be applied to resilient state
estimation in the bounded-error setting is that of unknown input interval
observers. Particularly, the works in [9, 10, 11] considered the problem of
designing unknown input interval observer s for continuous-time linear
parameter varying (LPV), uncertain linear time-invariant (LTI) and discrete-
time switched linear systems, respectively, where the authors  in [9]
formulated the necessary Metzler property as part of a semi-definite program.
A similar problem was considered for nonlinear continuous-time systems with
linear observations in [12]. However,  these approaches are not suitable for
general discrete-time nonlinear systems  and the unknown input signal  does
not affect the output/measurement equation (needed for representing false data
injection attacks on the sensors) in either of  the works in [9, 10, 11, 12].

On the other hand, while our previous works [13, 14] do consider the design of
state and unknown input interval observers for nonlinear discrete-time systems
with nonlinear observations, no stabilizing gains were synthesized in [13,
14].  We aim to address this shortcoming in this paper.

_Contributions._ By leveraging a combination of mixed-monotone decomposition
of nonlinear functions [15, 16] and parallel affine outer-approximation of
observation functions [17], we synthesize a resilient interval observer, i.e.,
a discrete-time dynamical system that by construction, _simultaneously_
returns interval-valued estimates of states and unknown inputs  (representing
false data injection signals on both the actuators and sensors) for  a broad
range of nonlinear discrete-time systems with nonlinear observations. Our
proposed design is a significant improvement to our previous input and state
interval observer designs [14, 13], in which no stabilizing gains were
considered and so the stability of the  previous observer designs only hinged
upon some dynamical systems properties. Moreover, in contrast to many unknown
input (interval) observer designs in the literature, our design considers
arbitrary unknown input signals with no assumptions of  _a priori_ known
intervals, being stochastic with zero mean (as is often assumed for noise) or
bounded. Further, we provide sufficient conditions for the input-to-state-
stability of the proposed observer, which at the same time ensures the
optimality of the design in the sense of  minimum $\mathcal{H}_{\infty}$ gain
by solving semi-definite programs.

##  II Preliminaries

_Notation._ $\lor$ denotes the logical disjunction (the OR truth-functional
operator). $\mathbb{R}^{n},\mathbb{R}^{n\times
p},\mathbb{D}_{n},\mathbb{N},\mathbb{N}_{n},{\mathbb{R}_{\geq 0}}$ and
$\mathbb{R}_{>0}$ denote the $n$-dimensional Euclidean space and the sets of
$n$ by $p$ matrices, $n$ by $n$ diagonal matrices, natural numbers  (including
0), natural numbers  from 1 to $n$, non-negative and positive real numbers,
respectively, while $\mathbb{M}_{n}$ denotes the set of all $n$ by $n$ Metzler
matrices, i.e., square matrices whose off-diagonal elements are non-negative.
Euclidean norm of a vector $x\in\mathbb{R}^{n}$ is denoted by
$\|x\|_{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\triangleq\sqrt{x^{\top}x}}$.
For $M\in\mathbb{R}^{n\times p}$, $M_{ij}$ denotes $M$’s entry in the $i$’th
row and the $j$’th column, $M^{\oplus}\triangleq\max(M,\mathbf{0}_{n,p})$,
$M^{\ominus}=M^{\oplus}-M$ and $|M|\triangleq M^{\oplus}+M^{\ominus}$, where
$\mathbf{0}_{n,p}$ is the zero matrix in $\mathbb{R}^{n\times p}$,  while
$\textstyle{\mathrm{sgn}}(M)\in\mathbb{R}^{n\times p}$ is the element-wise
sign of $M$ with $\textstyle{\mathrm{sgn}}(M_{ij})=1$ if $M_{ij}\geq 0$ and
$\textstyle{\mathrm{sgn}}(M_{ij})=-1$, otherwise. $M\succ 0$ and $M\prec 0$
(or $M\succeq 0$ and $M\preceq 0$) denote that $M$ is positive and negative
(semi-)definite, respectively. Further, a function
$f:S\subseteq\mathbb{R}^{n}\to\mathbb{R}$, where $0\in S$, is positive
definite if $f(x)>0$ for all $x\in S{\,\setminus\\{0\\}}$, and $f(0)=0$.
Finally, an interval
$\mathcal{I}\triangleq[\underline{{z}},\overline{{z}}]\subset\mathbb{R}^{n}$
is the set of all real vectors ${z\in\mathbb{R}^{n_{z}}}$ that satisfies
$\underline{{z}}\leq{z}\leq\overline{{z}}$ (component-wise), where
$\|\overline{{z}}-\underline{{z}}\|{{}_{\infty}\triangleq\max_{i\in\\{1,\cdots,{n_{z}}\\}}{|z_{i}|}}$
is the interval width of $\mathcal{I}$.

Next, we  review some related results and definitions.

######  Proposition 1 (Jacobian Sign-Stable Decomposition [15, Proposition
2]).

If a mapping $f:{\mathcal{Z}}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{p}$ has
Jacobian matrices satisfying
$J^{f}({x})\in[\underline{J}^{f},\overline{J}^{f}]$,
$\forall{z\in\mathcal{Z}}$, where
$\underline{J}^{f},\overline{J}^{f}\in\mathbb{R}^{p\times n_{z}}$ are known
matrices, then the mapping $f$ can be decomposed into an additive remainder-
form:

$\displaystyle\forall{z\in\mathcal{Z}},f({{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}z}})=H{z}+\mu({z}),$
(1)

where the matrix
$H\in\mathbb{R}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}p\times
n_{z}}}$ satisfies

$\displaystyle\forall(i,j)\in\mathbb{N}_{p}\times\mathbb{N}_{n_{z}},H_{ij}=\underline{J}^{f}_{ij}\
\lor H_{ij}=\overline{J}^{f}_{i,j},$ (2)

and $\mu(\cdot)$ and $H{z}$ are nonlinear and linear Jacobian sign-stable
(JSS) mappings, respectively, i.e., the signs of each element of their
Jacobian matrices do not change within their domains ($J^{\nu}_{ij}(\cdot)\geq
0\ \text{ or }J^{\nu}_{ij}(\cdot)\leq 0$, $\nu(z)\in\\{\mu(z),Hz\\}$).

######  Definition 1 (Mixed-Monotonicity and Decomposition Functions).

[18, Definition 1],[19, Definition 4] Consider the discrete-time dynamical
system $x_{k+1}=g(x_{k})$, with initial state
$x_{0}\in\mathcal{X}_{0}\triangleq[\underline{x}_{0},\overline{x}_{0}]{\subset\mathbb{R}^{n}}$.
Furthermore, ${g}:{\mathcal{X}}\subset\mathbb{R}^{n}\to\mathbb{R}^{{n}}$ is
the vector field, and $\mathcal{X}$ is the entire state space. A function
$g_{d}:{\mathcal{X}\times\mathcal{X}}\to\mathbb{R}^{n}$ is a discrete-time
mixed-monotone decomposition mapping for the vector field $g$ if it satisfies
the following conditions: i) $g_{d}({x,x})=g({x})$, ii)
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g_{d}}$
is monotone increasing in its first argument, i.e.,
$\hat{{x}}\geq{x}\Rightarrow g_{d}(\hat{{x}},{x}^{\prime})\geq
g_{d}(x,x^{\prime})$, and iii) $g_{d}$ is monotone decreasing in its second
argument, i.e., $\hat{x}\geq x\Rightarrow g_{d}({x}^{\prime},\hat{{x}})\leq
g_{d}({x}^{\prime},{x})$.

######  Proposition 2 (Tight and Tractable Decomposition Functions for JSS
Mappings).

[15, Proposition 4 & Lemma 3] Suppose
$\mu:{\mathcal{Z}}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{p}$ is a JSS mapping
on its domain. Then, for each $\mu_{i}$, $i\in\mathbb{N}_{p}$, its tight
decomposition function is:

$\displaystyle\mu_{d,i}({z}_{1},{z}_{2})=\mu_{i}(D^{i}{z}_{1}+(I_{n}-D^{i}){z}_{2}),$
(3)

for any ordered ${z_{1},z_{2}\in\mathcal{Z}}$, with a binary diagonal matrix
$D^{i}$ that is determined by the vertex of the interval $[{z}_{1},{z}_{2}]$
that minimizes the function $\mu_{i}$ (if $z_{1}<z_{2}$) or the vertex of the
interval $[z_{2},z_{1}]$ that maximizes $\mu_{i}$ (if $z_{2}\leq z_{1}$),
i.e.,

$D^{i}=\textstyle{\mathrm{diag}}(\max(\textstyle{\mathrm{sgn}}({\overline{J}^{\mu}_{i}}),\mathbf{0}_{1,{n_{z}}})).$

Moreover, if the JSS mapping $\mu$ is a remainder term of a JSS decomposition
of a function $f$ as discussed in Proposition 1, then for any interval domain
$\underline{z}\leq z\leq\overline{z}$, with
$z,\underline{z},\overline{z}\in\mathcal{Z}$ and
$\varepsilon\triangleq\overline{z}-\underline{z}$, the following inequality
holds:
$\delta^{\mu}_{d}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\triangleq\mu_{d}(\overline{z},\underline{z})-\mu_{d}(\underline{z},\overline{z})}\leq\overline{F}_{\mu}\varepsilon$,
with $\overline{F}_{\mu}\triangleq
2\max(\overline{J}_{f}-H,\mathbf{0}_{p,n_{z}})-\underline{J}_{f}+H$ and
$H\in\mathbb{R}^{p\times n_{z}}$ given in Proposition 1.

Consequently, by applying Proposition 2 to the Jacobian sign-stable
decomposition obtained using Proposition 1, a tight and tractable
decomposition function can be obtained (cf. details in [15]). Furthermore, in
the case that the mapping is not JSS, a tractable algorithm has been
introduced in [20, Algorithm 1] to compute _tight remainder-form decomposition
functions_ for a very broad class of nonlinear functions.

######  Definition 2 (Embedding System).

[16, Definition 6] For a discrete-time dynamical system $x_{k+1}=g(x_{k})$
defined over mapping
${g}:{\mathcal{X}}\subset\mathbb{R}^{n}\to\mathbb{R}^{{n}}$ with a
corresponding decomposition function ${g}_{d}(\cdot)$, its embedding system is
a $2n$-dimensional system with initial condition
$\begin{bmatrix}\overline{x}_{0}^{\top}&\underline{x}_{0}^{\top}\end{bmatrix}^{\top}$
defined as
$\begin{bmatrix}{\underline{x}}^{\top}_{k+1}&{\overline{x}}^{\top}_{k+1}\end{bmatrix}^{\top}=\begin{bmatrix}\underline{g}_{d}^{\top}(\underline{x}_{k},\overline{x}_{k})&\overline{g}_{d}^{\top}(\overline{x}_{k},\underline{x}_{k})\end{bmatrix}^{\top}$.

Note that according to [20, Proposition 3], the embedding system in Definition
2 with decomposition function $g_{d}$ corresponding to the dynamics
$x_{k+1}=g(x_{k})$ has a _state framer property_ , i.e., its solution is
guaranteed to frame the unknown state trajectory $x_{k}$, i.e.,
$\underline{x}_{k}\leq x_{k}\leq\overline{x}_{k}$ for all $k\in\mathbb{N}$.

Next, we will briefly restate our previous result in [17], tailoring it
specifically for intervals to help with computing affine bounding functions
for our functions.

######  Proposition 3.

[17, Affine Outer-Approximation] Consider the  function
$g(.):\mathcal{B}\subset\mathbb{R}^{n}\to\mathbb{R}^{m}$, where $\mathcal{B}$
is an interval with $\overline{x},\underline{x},\mathcal{V}_{\mathcal{B}}$
being its maximal, minimal and set of vertices, respectively. Suppose
$\overline{A}_{\mathcal{B}},\underline{A}_{\mathcal{B}},\overline{e}_{\mathcal{B}},\underline{e}_{\mathcal{B}},\theta_{\mathcal{B}}$
is a solution of the following  linear program (LP):

$\displaystyle\min_{\theta,\overline{A},\underline{A},\overline{e},\underline{e}}{\theta}$
(4) $\displaystyle\quad\quad s.t\
\underline{A}{x}_{s}+\underline{e}+\sigma\leq
g({x}_{s})\leq\overline{A}{x}_{s}+\overline{e}-\sigma,$
$\displaystyle\quad\quad\quad\
(\overline{A}-\underline{A}){x}_{s}+\overline{e}-\underline{e}-2\sigma\leq\theta\mathbf{1}_{m},\
\forall x_{s}\in\mathcal{V}_{\mathcal{B}},$

where $\mathbf{1}_{m}\in\mathbb{R}^{m}$ is a vector of ones and $\sigma$ can
be computed via  [17, Proposition 1] for different function classes. Then,
$\underline{A}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{B}}}{x}+\underline{e}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{B}}}\leq
g(x)\leq\overline{A}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{B}}}{x}+\overline{e}_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{B}}},\forall
x\in\mathcal{B}$.

######  Corollary 1.

By taking the average of upper and lower affine abstractions and
adding/subtracting half of the maximum distance, it is straightforward to
“parallelize” the above upper and lower abstractions as
${A}_{g}{x}+\underline{\epsilon}\leq g(x)\leq A_{g}{x}+\overline{\epsilon}$,
or equivalently
$g(x)=A_{g}x+\epsilon,\epsilon\in[\underline{\epsilon},\overline{\epsilon}]$,
where $A_{g}\triangleq(1/2)(\overline{A}+\underline{A})$,
$\underline{\epsilon}\triangleq(1/2)(\overline{e}+\underline{e}-\theta\mathbf{1}_{m})$
and
$\overline{\epsilon}\triangleq(1/2)(\overline{e}+\underline{e}+\theta\mathbf{1}_{m})$.
We call $A_{g}$ and $\epsilon$ the parallel affine outer-approximation slope
and outer-approximation error of function $g$ on $\mathcal{B}$, respectively.

##  III Problem Formulation

_System Assumptions._ Consider the nonlinear discrete-time system with unknown
inputs and bounded noise

$\displaystyle\begin{array}[]{rl}x_{k+1}&=f(x_{k})+Ww_{k}+Gd_{k},\\\
y_{k}&=h(x_{k})+Vv_{k}+Hd_{k},\end{array}$ (7)

where at time  $k\in\mathbb{N}$, $x_{k}\in\mathcal{X}\subset\mathbb{R}^{n}$,
$d_{k}\in\mathbb{R}^{p}$ and $y_{k}\in\mathbb{R}^{l}$ are the state vector,
unknown input vector, and measurement vector, respectively. The process and
measurement noise signals $w_{k}\in\mathbb{R}^{n}$ and
$v_{k}\in\mathbb{R}^{l}$ are assumed to be bounded,  i.e.,
$w_{k}\in\mathcal{W}\triangleq[\underline{w},\overline{w}\\}$,
$v_{k}\in\mathcal{V}\triangleq[\underline{v},\overline{v}]$ with known lower
and upper bounds, $\underline{w}$, $\overline{w}$ and $\underline{v}$,
$\overline{v}$, respectively. We also assume  that lower and upper bounds for
the initial state, $\underline{x}_{0}$ and $\overline{x}_{0}$, are available,
i.e., $\underline{x}_{0}\leq x_{0}\leq\overline{x}_{0}$. The functions
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$,
$h:\mathbb{R}^{n}\rightarrow\mathbb{R}^{l}$ and matrices $W$, $V$, $G$ and $H$
are known and of appropriate dimensions, where $G$ and $H$ encode the
_locations_ at which the unknown input  (or attack) signal can affect the
system dynamics and measurements. Note that no assumption is made on $H$ to be
either the zero matrix (no direct feedthrough), or to have full column rank
when there is direct feedthrough (in contrast  to [13]).

_Unknown Input (or Attack) Signal Assumptions._ The unknown inputs $d_{k}$
(representing false data injection attack signals) are not constrained to
follow any model nor to be a signal of any type (random or strategic), hence
no prior ‘useful’ knowledge of the dynamics of $d_{k}$ is available
(independent of $\\{d_{\ell}\\}$ $\forall k\neq\ell$, $\\{w_{\ell}\\}$ and
$\\{v_{\ell}\\}$ $\forall\ell$). We also do not assume that $d_{k}$ is bounded
or has known bounds and thus, $d_{k}$ is suitable for representing adversarial
attack signals.

Next, we briefly introduce a similar system transformation as in [3], which
will be used later in our observer structure.

_System Transformation._ Let $p_{H}\triangleq{\rm rk}(H)$. Similar to [3], by
applying singular value decomposition, we have
$H=\begin{bmatrix}U_{1}&U_{2}\end{bmatrix}\begin{bmatrix}\Xi&0\\\
0&0\end{bmatrix}\begin{bmatrix}E_{1}^{\,\top}\\\ E_{2}^{\,\top}\end{bmatrix}$
with $E_{1}\in\mathbb{R}^{p\times p_{H}}$,
$E_{2}\in\mathbb{R}^{p\times(p-p_{H})}$, $\Xi\in\mathbb{R}^{p_{H}\times
p_{H}}$ (a diagonal matrix of full rank; so we can define
$S\triangleq\Xi^{-1}$), $U_{1}\in\mathbb{R}^{l\times p_{H}}$ and
$U_{2}\in\mathbb{R}^{l\times(l-p_{H})}$. Then, since
$D\triangleq\begin{bmatrix}E_{1}&E_{2}\end{bmatrix}$ is unitary:

$\displaystyle d_{k}=E_{1}d_{1,k}+E_{2}d_{2,k},\ d_{1,k}=E_{1}^{\top}d_{k},\
d_{2,k}=E_{2}^{\top}d_{k}.$ (8)

Finally, by defining $T_{1}\triangleq U^{\top}_{1},T_{2}\triangleq
U^{\top}_{2}$, the output equation can be decoupled, by which system (7) can
be rewritten as:

$\displaystyle\begin{array}[]{rl}x_{k+1}&=f(x_{k})+Ww_{k}+G_{1}d_{1,k}+G_{2}d_{2,k},\\\
z_{1,k}&=h_{1}(x_{k})+V_{1}v_{k}+\Xi d_{1,k},\\\
z_{2,k}&={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h_{2}(x_{k})+V_{2}v_{k}},\end{array}$
(12)

where $h_{i}(x_{k})=T_{i}h(x_{k})$, $\forall i\in\\{1,2\\}$ and
$K_{i}\triangleq T_{i}K_{i},\forall K\in\\{V,G\\},\forall i\in\\{1,2\\}$.

Moreover,  we assume the following, which is satisfied for a broad range of
nonlinear functions [21]:

######  Assumption 1.

Functions $f,h$ have bounded Jacobians over the state space $\mathcal{X}$ with
known/computable Jacobian bounds.

######  Assumption 2.

The JSS decomposition of $h_{2}(x_{k})$ via Proposition 1 given by
$h_{2}(x_{k})=C_{2}x_{k}+\psi_{2}(x_{k})$ is such that $\psi_{2}$ is JSS and
further, $C_{2}G_{2}$ has full column rankaaa In the special case that $G=0$,
we would require $G_{2}$ to be empty (and this does happen when $H$ has full
rank), in which case $C_{2}G_{2}$ being full rank is satisfied by assumption..
Consequently, there exists $M_{2}\triangleq(C_{2}G_{2})^{\dagger}$ such that
$M_{2}C_{2}G_{2}=I$.

######  Assumption 3.

(Only needed when the observations are nonlinear, i.e., if
$\psi_{2}(x_{k})\neq 0$)  The entire state space
$\mathcal{X}\subset\mathbb{R}^{n}$ is bounded. Moreover, $A_{g}$ is
invertible, where $A_{g}\in\mathbb{R}^{n\times n}$ is the parallel affine
outer-approximation slope (cf. Proposition 3 and Corollary 1) of the function
$g(x)\triangleq
x+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{2}}M_{2}\psi_{2}(x)$
over the entire state space.

Further, we formally define the notions of _framers_ , _correctness_ and
_stability_ that are used throughout the paper.

######  Definition 3 (Interval Framers).

Given the nonlinear plant (7) (equivalently (12)), the sequences
$\\{\overline{x}_{k},\underline{x}_{k}\\}_{k=0}^{\infty}\subset\mathbb{R}^{n}$
and
$\\{\overline{d}_{k},\underline{d}_{k}\\}_{k=0}^{\infty}\subset\mathbb{R}^{p}$
are called upper and lower framers for the states and inputs of the system in
(7), respectively, if

$\displaystyle\forall
k\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{N}},{\forall
w_{k}\in\mathcal{W},\forall v_{k}\in\mathcal{V}},\
\underline{\nu}_{k}\leq\nu_{k}\leq\overline{\nu}_{k},\forall\nu\in\\{x,d\\}.$

In other words, starting from the initial interval $\underline{x}_{0}\leq
x_{0}\leq\overline{x}_{0}$, the true state of the system in (7), $x_{k}$, and
the unknown input $d_{k}$ are guaranteed to evolve within the interval flow-
pipe $[\underline{x}_{k},\overline{x}_{k}]$ and bounded within the interval
$[\underline{d}_{k},\overline{d}_{k}]$, for all
$(k,{w_{k},v_{k})}\in{\mathbb{N}\times\mathcal{W}\times\mathcal{V}}$,
respectively. Finally, any dynamical system (i.e., tractable algorithm) that
returns upper and lower framers for the states and unknown inputs of system 7
is called a resilient interval framer for (7).

######  Definition 4 (Framer Error).

Given state and input framers
$\\{\underline{x}_{k}\leq\overline{x}_{k}\\}_{k=0}^{\infty}$ and
$\\{\underline{d}_{k}\leq\overline{d}_{k}\\}_{k=1}^{\infty}$, the sequences
$\\{e^{x}_{k}\triangleq\overline{x}_{k}-\underline{x}_{k}\\}_{k=0}^{\infty}$
and
$\\{e^{d}_{k}\triangleq\overline{d}_{k}-\underline{d}_{k}\\}_{k=1}^{\infty}$
are called the state and input framer errors, respectively. It easily  follows
from Definition 3 that $e^{\nu}_{k}\geq 0,\forall
k\in{\mathbb{N}},\forall\nu\in\\{x,d\\}.$

######  Definition 5 (Input-to-State Stability and Interval Observer).

An interval framer is input-to-state stable (ISS), if the framer state error
(cf. Definition 4) is bounded as follows:

$\displaystyle\forall k\in\mathbb{N},\
\|e^{x}_{k}\|_{2}\leq\beta(\|e^{x}_{0}\|_{2},k)+\alpha(\|\delta\|_{\ell_{\infty}}),$
(13)

where $\delta\triangleq[(\delta^{w})^{\top}\
(\delta^{v})^{\top}]^{\top}\triangleq[(\overline{w}-\underline{w})^{\top}\
(\overline{v}-\underline{v})^{\top}]^{\top}$, $\beta$ and $\alpha$ are
functions of classesbbbA function $\alpha:\mathbb{R}_{\geq
0}\to\mathbb{R}_{\geq 0}$ is of class $\mathcal{K}$ if it is continuous,
positive definite, and strictly increasing and is of class
$\mathcal{K}_{\infty}$ if it is also unbounded. Moreover,
$\lambda:\mathbb{R}_{\geq 0}\to\mathbb{R}_{\geq 0}$ is of class $\mathcal{KL}$
if for each fixed $t\geq 0$, $\lambda(\cdot,t)$ is of class $\mathcal{K}$ and
for each fixed $s\geq 0$, $\lambda(s,t)$ decreases to zero as $t\to\infty.$
$\mathcal{KL}$ and $\mathcal{K}_{\infty}$, respectively, and
$\|\delta\|_{\ell_{\infty}}\triangleq\sup_{k\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathbb{N}}}\|\delta_{k}\|_{2}=\|\delta\|_{2}$
is the $\ell_{\infty}$ signal norm. An ISS  resilient interval framer is
called  a resilient interval observer.

######  Definition 6 ($\mathcal{H}_{\infty}$-Optimal Resilient Interval
Observer).

A resilient interval framer design is $\mathcal{H}_{\infty}$-optimal if the
$\mathcal{H}_{\infty}$ gain of the framer error system $\tilde{\mathcal{G}}$,
i.e.,
$\|\tilde{\mathcal{G}}\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{H}_{\infty}}}$
is minimized, where
$\|\tilde{\mathcal{G}}\|_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{H}_{\infty}}}\triangleq\sup{\\{}\frac{\|e^{x}\|_{\ell_{2}}}{\|\delta\|_{\ell_{2}}},\delta\neq
0{\\}}$, and
$\|s\|_{\ell_{2}}\triangleq\sqrt{\sum_{0}^{\infty}\|s_{k}\|_{2}^{2}}$ is the
$\ell_{2}$ signal norm for $s\in\\{e^{x},\delta\\}$.

Using the above, we aim to address the following problem.

######  Problem 1.

Given the nonlinear system in (7), as well as Assumptions 1–3, synthesize an
ISS and $\mathcal{H}_{\infty}$-optimal resilient interval observer (cf.
Definitions 3–6).

##  IV Resilient Interval Observer Design

In this section, we describe the proposed resilient interval observer as well
as analyze its correctness and ISS properties.

###  IV-A Interval Framer Design

Our strategy for designing resilient interval observers in the presence of
unknown inputs has three steps. First, we obtain an equivalent representation
of the system in (7) by introducing some auxiliary state variables, such that
the equivalent system is not affected by the attack signal. Then, inspired by
our previous work on synthesizing interval observers for nonlinear systems
[15, 16] we will design embedding systems (cf. Definition 2) for the
equivalent system representation, which returns state framers. Finally, we
obtain input framers  (with a one-step delay since $d_{2,k}$ does not appear
in the measurements $z_{1,k}$ and $z_{2,k}$ in (12)) as functions of the
computed state framers.

First, note that from (12) and with $S\triangleq\Xi^{-1}$, $d_{1,k}$ can be
computed as a function of the state at current time as follows:

$\displaystyle d_{1,k}=S(z_{1,k}-h_{1}(x_{k})-V_{1}v_{k}).$ (14)

Next, we introduce an auxiliary state variable as:

$\displaystyle\xi_{k}\triangleq
x_{k}-N(z_{2,k}-V_{2}v_{k}-\psi_{2}(x_{k}))=(I-NC_{2})x_{k},$ (15)

where  the equality follows from (12) and Assumption 2. Moreover,
$N\in\mathbb{R}^{n\times(l-\tilde{p})}$ is a to-be-designed gain to cancel out
the effect of the unknown input in the state equation. This is done through
the following lemma.

######  Lemma 1.

Suppose Assumption 2 holds and let  $N=G_{2}M_{2}=G_{2}(C_{2}G_{2})^{\dagger}$
and $S\triangleq\Xi^{-1}$. Then, the value of the auxiliary state $\xi_{k}$ at
time step $k+1$ can be computed as:

$\displaystyle\xi_{k+1}$
$\displaystyle=(I-NC_{2})(f(x_{k})+G_{1}S(z_{1,k}-h_{1}(x_{k})-V_{1}v_{k})+Ww_{k}).$
(16)

###### Proof.

By plugging $d_{1,k}$ from (14) into (12), we obtain

$\displaystyle
x_{k+1}=f(x_{k})+G_{1}S(z_{1,k}-h_{1}(x_{k})-V_{1}v_{k})+Ww_{k}+G_{2}d_{2,k}.$
(17)

This, together with the second equality in (15) and  the above choice of $N$
such that $(I-NC_{2})G_{2}=0$, returns (16). ∎

The evolution of the auxiliary state $\xi_{k}$ in (16) is independent of the
unknown input and hence, we can compute propagated framers for $\xi_{k}$
leveraging embedding systems  (cf. Proposition 2). However, we do not have a
way of  directly retrieving the propagated framers for the original states,
i.e., $\\{\underline{x}_{k},\overline{x}_{k}\\}$ in terms of
$\\{\underline{\xi}_{k},\overline{\xi}_{k}\\}$ from the second equality of
(15), since $I-NC_{2}=I-G_{2}(C_{2}G_{2})^{\dagger}C_{2}$ can be shown to be
not invertible. To overcome this difficulty, given Assumption 3, we introduce
a new auxiliary state:

$\displaystyle\gamma_{k}\triangleq
x_{k}-\Lambda(N(z_{2,k}-V_{2}v_{k})-\epsilon_{k}),$ (18)

with $\Lambda\triangleq A_{g}^{-1}$, where $A_{g}$ and
$\epsilon_{k}\in[{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underline{\epsilon}},\overline{\epsilon}]$
are parallel affine outer-approximation slope and approximation error of the
mapping $g(x)\triangleq
x+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{2}}M_{2}\psi_{2}(x)$
on the entire space $\mathcal{X}$ (cf. Proposition 3, Corollary 1 and
Assumption 3).

######  Proposition 4.

Given Assumption 3, the two auxiliary states $\gamma_{k}$ and $\xi_{k}$ are
linearly related as:

$\displaystyle\gamma_{k}=\Lambda\xi_{k}.$ (19)

###### Proof.

Computing parallel affine outer-approximation of the mapping
$g(x_{k})=A_{g}x_{k}+\epsilon_{k}$ and applying (15), we obtain

$\displaystyle g(x_{k})$
$\displaystyle\triangleq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}x_{k}}+N\psi_{2}(x)=\xi_{k}+N(z_{2,k}-V_{2}v_{k})$
$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Rightarrow}A_{g}x_{k}$
$\displaystyle=\xi_{k}+N(z_{2,k}-V_{2}v_{k})-\epsilon_{k},\quad\epsilon_{k}\in[\underline{\epsilon},\overline{\epsilon}],$

from which  and given Assumption 3 ( that $A_{g}$ is invertible, with
$\Lambda=A_{g}^{-1}$), we have

$\displaystyle
x_{k}=\Lambda(\xi_{k}+N(z_{2,k}-V_{2}v_{k})-\epsilon_{k}),\quad\epsilon_{k}\in[\underline{\epsilon},\overline{\epsilon}].$
(20)

Plugging $x_{k}$ from (20) into (18) returns the results. ∎

We are now ready to propose an input and state resilient interval framer,
i.e., the following discrete-time dynamical system (27)–(35), which by
construction, outputs /returns framers for the original states
$\\{x_{k}\\}_{k=0}^{\infty}$ and the unknown input signal
$\\{d_{k}\\}_{k=1}^{\infty}$ of system (7). The details of the framer
construction/design will be provided in the proof of Theorem 1. The proposed
resilient interval framer is as follows:

$\displaystyle\begin{array}[]{rlllll}\underline{\gamma}_{k+1}&=(A-LC_{2})^{\oplus}\underline{\gamma}_{k}-(A-LC_{2})^{\ominus}\overline{\gamma}_{k}+\rho_{d}(\underline{x}_{k},\overline{x}_{k})\\\
&+D^{\ominus}\underline{\epsilon}-D^{\oplus}\overline{\epsilon}+L^{\ominus}\psi_{2,d}(\underline{x}_{k},\overline{x}_{k})-L^{\oplus}\psi_{2,d}(\overline{x}_{k},\underline{x}_{k})\\\
&+\hat{V}^{\ominus}\underline{v}-\hat{V}^{\oplus}\overline{v}+\hat{W}^{\ominus}\underline{w}-\hat{W}^{\ominus}\overline{w}+\hat{z}_{k},\\\
\overline{\gamma}_{k+1}&=(A-LC_{2})^{\oplus}\overline{\gamma}_{k}-(A-LC_{2})^{\ominus}\underline{\gamma}_{k}+\rho_{d}(\overline{x}_{k},\underline{x}_{k})\\\
&+D^{\ominus}\overline{\epsilon}-D^{\oplus}\underline{\epsilon}+L^{\ominus}\psi_{2,d}(\overline{x}_{k},\underline{x}_{k})-L^{\oplus}\psi_{2,d}(\underline{x}_{k},\overline{x}_{k})\\\
&+\hat{V}^{\ominus}\overline{v}-\hat{V}^{\oplus}\underline{v}+\hat{W}^{\oplus}\overline{w}-\hat{W}^{\ominus}\underline{w}+\hat{z}_{k},\end{array}$
(27)

$\displaystyle\begin{array}[]{rl}\underline{x}_{k}&=\underline{\gamma}_{k}+\Lambda
Nz_{2,k}+\Lambda^{\ominus}\underline{\epsilon}-\Lambda^{\oplus}\overline{\epsilon}+(\Lambda
NV_{2})^{\ominus}\underline{v}-(\Lambda NV_{2})^{\oplus}\overline{v},\\\
\overline{x}_{k}&=\overline{\gamma}_{k}+\Lambda
Nz_{2,k}+\Lambda^{\ominus}\overline{\epsilon}-\Lambda^{\oplus}\underline{\epsilon}+(\Lambda
NV_{2})^{\ominus}\overline{v}-(\Lambda
NV_{2})^{\oplus}\underline{v},\end{array}$ (30)
$\displaystyle\begin{array}[]{rlll}\underline{d}_{k-1}&=\Phi^{\oplus}\underline{x}_{k}-\Phi^{\ominus}\overline{x}_{k}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa_{d}}(\underline{x}_{k-1},\overline{x}_{k-1})+A_{z}z_{1,k-1}\\\
&+A_{v}^{\oplus}\underline{v}-A_{v}^{\ominus}\overline{v}+\Phi^{\ominus}\underline{w}-\Phi^{\oplus}\overline{w},\\\
\overline{d}_{k-1}&=\Phi^{\oplus}\overline{x}_{k}-\Phi^{\ominus}\underline{x}_{k}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa_{d}}(\overline{x}_{k-1},\underline{x}_{k-1})+A_{z}z_{1,k-1}\\\
&+A_{v}^{\oplus}\overline{v}-A_{v}^{\ominus}\underline{v}+\Phi^{\ominus}\overline{w}-\Phi^{\oplus}\underline{w},\end{array}$
(35)

where $S\triangleq\Xi^{-1}$, $N=G_{2}M_{2}$ and $\Lambda\triangleq
A^{-1}_{g}$. Furthermore, $L\in\mathbb{R}^{n\times(l-p_{H})}$ is an arbitrary
matrix (observer gain) which will be designed later in Theorem 2 to  yield
stability and optimality of the proposed framers. Moreover,
$A\in\mathbb{R}^{n\times n}$ and
$\rho:\mathcal{X}\subset\mathbb{R}^{n}\to\mathbb{R}^{n}$ are obtained by
applying JSS decompositions (cf. Proposition 1) on the mapping
$\tilde{f}(x)\triangleq\Lambda(I-NC_{2})(f(x)-G_{1}Sh_{1}(x))$, while
$\psi_{2,d}$ and $\rho_{d}$ are tight decomposition functions for the JSS
mappings $\rho$ and $\psi_{2}$, respectively, computed through Proposition 1.
Further,

$\displaystyle\begin{array}[]{rlll}\hat{V}&\triangleq(A-LC_{2})\Lambda
NV_{2}+LV_{2}+\Lambda(I-NC_{2})G_{1}SV_{1},\\\ \Phi&\triangleq
E_{2}M_{2}C_{2},A_{v}\triangleq(\Phi
G_{1}-E_{1})SV_{1},A_{z}\triangleq(E_{1}-\Phi G_{1})S,\\\
D&\triangleq(A-LC_{2})\Lambda,\hat{W}\triangleq\Lambda(I-NC_{2})W,\\\
\hat{z}_{k}&\triangleq\Lambda(I-NC_{2})G_{1}Sz_{1,k}+(L+(A-LC_{2})\Lambda
N)z_{2,k},\end{array}$ (40)

and
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa_{d}}$
is the decomposition function of the mapping
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa}(x)\triangleq(\Phi
G_{1}-E_{1})Sh_{1}(x)-\Phi f(x)$, computed via [19, Theorem 1]. Finally,
$A_{g}$ and $\epsilon_{k}\in[\underline{\epsilon},\overline{\epsilon}]$ are
computed  via Corollary 1.

The following theorem formalizes the state and input framer/correctness
property of the proposed resilient interval observer (27)–(35) with respect to
the original system (7).

######  Theorem 1.

Suppose Assumptions 1–3 hold. Then, the sequences
$\\{\underline{x}_{k},\overline{x}_{k}\\}_{k=0}^{\infty}$ and
$\\{\underline{d}_{k},\overline{x}_{k}\\}_{k=1}^{\infty}$ obtained from the
system (27)–(35), construct framers for the states and unknown input signal of
(7), respectively, i.e.,
$\underline{\nu}_{k}\leq\nu_{k}\leq\overline{\nu}_{k},\forall\nu\in\\{x,d\\},\forall
k\in\mathbb{N},\forall w_{k}\in\mathcal{W},\forall v_{k}\in\mathcal{V}$.

###### Proof.

From (16) and (19), we obtain

$\displaystyle\gamma^{+}_{k}=\Lambda(I-NC_{2})(f(x_{x})+G_{1}S(z_{1,k}-h_{1}(x_{k})-V_{1}v_{k})+Ww_{k}).$

Adding the _zero term_ $L(z_{2,k}-C_{2}x_{k}-\psi_{2}(x_{k})-V_{2}v_{k})$ to
the right hand side of  the above equation and applying mixed-monotone
decompositions  on the mapping
$\tilde{f}(x)\triangleq\Lambda(I-NC_{2})(f(x)-G_{1}Sh_{1}(x))$ to decompose it
as $\tilde{f}(x)=Ax+\rho(x)$ (cf. Proposition 1 for more details), yields:

$\displaystyle\gamma_{k+1}=(A-LC_{2})x_{k}+\rho(x_{k})-L\psi_{2}(x_{k})-\tilde{V}v_{k}+\tilde{W}w_{k}+\tilde{z}_{k},$
(41)

where $\tilde{V}\triangleq\Lambda(I-NC_{2})G_{1}SV_{1}+LV_{2}$,
$\tilde{W}\triangleq\Lambda(I-NC_{2})W$ and
$\tilde{z}_{k}\triangleq\Lambda(I-NC_{2})G_{1}Sz_{1,k}+Lz_{2,k}$. Then, by
computing $x_{k}$ in terms of $\gamma_{k}$ from (18) and plugging it back into
the linear terms in the right-hand side of (41), we obtain

$\displaystyle\small\begin{array}[]{rl}\gamma_{k+1}&=(A-LC_{2})\gamma_{k}+\rho(x_{k})-L\psi_{2}(x_{k})\\\
&-\hat{V}v_{k}+\hat{W}w_{t}-D\epsilon_{k}+\hat{z}_{k},\end{array}\normalsize$
(44)

with $\hat{V},D,\hat{W}$ and $\hat{z}_{k}$ given in (40). Next, by applying
Proposition 2 and [22, Lemma 1], we construct the embedding system (27) for
(44), which implies
$\underline{\gamma}_{k}\leq\gamma_{k}\leq\overline{\gamma}_{k},\forall
k\in\mathbb{N}$, by construction. Further, the results in (30) follow from
applying [22, Lemma 1] on (18) to compute framers of $x_{k}$ in terms of the
framers of $\gamma_{k}$.

To obtain input framers, note that multiplying both sides of (17) by
$M_{2}C_{2}$ together with Assumption 2 yields
$d_{2,k-1}=M_{2}C_{2}(x_{k}-f(x_{k-1})+G_{1}Sh_{1}(x_{k-1})+G_{1}S(V_{1}v_{k-1}-z_{1,k-1})-Ww_{k-1})$.
This, along with (8) and (14),  leads to

$\displaystyle d_{k-1}=\Phi
x_{k}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa(}x_{k-1})+A_{z}z_{1,k-1}+A_{v}v_{k-1}-\Phi
Ww_{k-1}.$ (45)

The input framers in (35) are obtained by leveraging [19, Theorem 1] to
compute a decomposition function for the nonlinear function  $\kappa$, as well
as  applying [22, Lemma 1] to bound the linear terms in the right-hand side of
(45). ∎

###  IV-B ISS and $\mathcal{H}_{\infty}$-Optimal Interval Observer Synthesis

Next, we provide sufficient conditions to guarantee the stability of the
proposed framers, i.e., we  seek to synthesize the observer gain $L$ to ensure
input-to-state stability (ISS) of the observer state error,
$e^{x}_{k}\triangleq\overline{x}_{k}-\underline{x}_{k}$ in the sense of
Definition 4, while ensuring that the design is optimal in the sense of
minimizing the $\mathcal{H}_{\infty}$ gain (cf. Definition 5).

First, we derive the observer error dynamics as follows.

######  Lemma 2.

Consider the nonlinear system (7) and suppose all assumptions in Theorem 2
hold. Then, the state framer error dynamics of the resilient interval observer
(27)–(35) and its nonlinear comparison system are as follows:

$\displaystyle\begin{array}[]{rll}e^{x}_{k+1}&=|A-LC_{2}|e^{x}_{k}+\delta^{\rho}_{k}+|L|\delta^{\psi_{2}}_{k}+|\hat{W}|\delta^{w}\\\
&+(|V_{a}-LV_{b}|-|A-LC_{2}||\Lambda NV_{2}|+|\Lambda NV_{2}|)\delta^{v}\\\
&+(|\Lambda|+|D_{a}-LD_{b}|-|A-LC_{2}||\Lambda|)\delta^{\epsilon}\\\
&\leq(|A-LC_{2}|+\overline{F}_{\rho}+|L|\overline{F}_{\psi_{2}})e^{x}_{k}+|\hat{W}|\delta^{w}\\\
&+(|V_{a}-LV_{b}|-|A-LC_{2}||\Lambda NV_{2}|+|\Lambda NV_{2}|)\delta^{v}\\\
&+(|\Lambda|+|D_{a}-LD_{b}|-|A-LC_{2}||\Lambda|)\delta^{\epsilon},\end{array}$
(52)

where
$\delta^{\zeta}_{k}\triangleq\zeta_{d}(\overline{x}_{k},\underline{x}_{k})-\zeta_{d}(\underline{x}_{k},\overline{x}_{k}),\forall\zeta\in\\{\psi_{2},\rho\\}$,
$\delta^{s}\triangleq\overline{s}-\underline{s},\forall
s\in\\{w,v,\epsilon\\}$, and
$\overline{F}_{\zeta},\forall\zeta\in\\{\psi_{2},\rho\\}$ are computed through
Proposition 2. Moreover,

$\displaystyle\begin{array}[]{rl}V_{a}&\triangleq A\Lambda
NV_{2}+\Lambda(I-NC_{2})G_{1}SV_{1},\\\ V_{b}&\triangleq(C_{2}\Lambda
N-I)V_{2},D_{a}\triangleq A\Lambda,D_{b}\triangleq C_{2}\Lambda.\end{array}$
(55)

###### Proof.

It follows from (27) that  the dynamics of
$e^{\gamma}_{k}\triangleq\overline{\gamma}_{k}-\underline{\gamma}_{k}$ is
given by
$e^{\gamma}_{k+1}=|A-LC_{2}|e^{\gamma}_{k}+\delta^{\rho}_{k}+|L|\delta^{\psi_{2}}_{k}+|\hat{V}|\delta^{v}+|\hat{W}|\delta^{w}+|D|\delta^{\epsilon}$.
This, combined with
$e^{x}_{k}=e^{\gamma}_{k}+|\Lambda|\delta^{\epsilon}+|\Lambda
NV_{2}|\delta^{v}$ (followed from (30)) results in the equality in (52),
which together with the facts that
$\delta^{\zeta}_{k}\leq\overline{F}_{\zeta}e^{x}_{k},\forall\zeta\in\\{\rho,\psi_{2}\\}$
(cf. Proposition 2), yields the inequality in (52). ∎

Further, by leveraging slightly different approaches to derive an upper
_linear_ comparison system for the _nonlinear_ error comparison system (52),
we derive different sets of sufficient conditions to guarantee the ISS
property of the proposed observer, as well as to ensure the optimality of the
design in the sense of  minimum $\mathcal{H}_{\infty}$ gain,  as follows.

######  Theorem 2 (ISS & $\mathcal{H}_{\infty}$-Optimal Resilient Interval
Observer Synthesis).

Consider system (7) (equivalently the transformed system (12)) and suppose
Assumptions 1–3 hold. Moreover, suppose there exist matrices
$\mathbb{R}^{n\times n}\ni
P^{*}\succ\mathbf{0}_{n,n},{\Gamma^{*}\in\mathbb{R}^{n\times(l-p_{H})}_{\geq
0}}$ and $\eta^{*}\in\mathbb{R}_{>0}$ such that $-P^{*}\in\mathbb{M}_{n}$ and
the tuple $(P^{*},\Gamma^{*},\eta^{*})$ solves the following problem:

$\displaystyle\begin{array}[]{rl}&\min\limits_{\\{\eta,P,\Gamma\\}}\eta\\\
&s.t.\begin{bmatrix}P&P\tilde{A}-\Gamma\tilde{C}&P\tilde{B}-\Gamma\tilde{D}&0\\\
*&P&0&I\\\ *&*&\eta I&0\\\ *&*&*&\eta I\end{bmatrix}\succ
0,(P,\Gamma)\in\mathbf{C},\\\ \end{array}$ (58)

where the matrices $\tilde{A},\tilde{B},\tilde{C},\tilde{D}$, as well as the
corresponding additional set of constraints  $\mathbf{C}$ can be either of the
following:

  1. (i)

$\mathbf{C}=\\{(P,\Gamma)\mid
P\begin{bmatrix}{A}&V_{a}&D_{a}\end{bmatrix}-\Gamma\begin{bmatrix}{C}_{2}&V_{b}&D_{b}\end{bmatrix}\geq
0\\}$, if:

$\displaystyle\tilde{A}$ $\displaystyle=A+\overline{F}_{\rho},\
\tilde{C}=C_{2}-\overline{F}_{\psi_{2}},$ $\displaystyle\tilde{B}$
$\displaystyle=\begin{bmatrix}V_{a}+(I-A)|\Lambda
NV_{2}|&|\hat{W}|&D_{a}+(I-A)|\Lambda|\end{bmatrix},$ $\displaystyle\tilde{D}$
$\displaystyle=\begin{bmatrix}V_{b}-C_{2}|\Lambda
NV_{2}|&0&D_{b}-C_{2}|\Lambda|\end{bmatrix}.$

  2. (ii)

$\mathbf{C}=\\{(P,\Gamma)\mid\Gamma\begin{bmatrix}{C}_{2}&V_{b}&D_{b}\end{bmatrix}\geq
0\\}$, if

$\displaystyle\tilde{A}$ $\displaystyle=|A|+\overline{F}_{\rho},\
\tilde{C}=-C_{2}-\overline{F}_{\psi_{2}},$ $\displaystyle\tilde{B}$
$\displaystyle=\begin{bmatrix}|V_{a}|+(I-|A|)|\Lambda
NV_{2}|&|\hat{W}|&(I-|A|)|\Lambda|+|D_{a}|\end{bmatrix},$
$\displaystyle\tilde{D}$ $\displaystyle=\begin{bmatrix}C_{2}|\Lambda
NV_{2}|-V_{b}&0&C_{2}|\Lambda|-D_{b}\end{bmatrix}.$

  3. (iii)

$\mathbf{C}=\\{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(P,\Gamma)\mid\,}PA-\Gamma
C_{2}\geq 0$}, if:

$\displaystyle\tilde{A}$ $\displaystyle=A+\overline{F}_{\rho},\
\tilde{C}=C_{2}-\overline{F}_{\psi_{2}},\
\tilde{D}=\begin{bmatrix}-V_{2}&0&0\end{bmatrix},$ $\displaystyle\tilde{B}$
$\displaystyle=\begin{bmatrix}|\Lambda(I-NC_{2})G_{1}SV_{1}|+|\Lambda
NV_{2}|&|\hat{W}|&|\Lambda|\end{bmatrix}.$

Then, the proposed resilient interval framer (27)–(35) with the corresponding
gain $L=({P}^{*})^{-1}\Gamma^{*}$, is a resilient ISS input and state interval
observer in the sense of Definition 5 and also is
$\mathcal{H}_{\infty}$-optimal (cf. Definition 6).  Finally, in any of the
above cases, the LMI in (58) is feasible only if the linear comparison system
$(\tilde{A},\tilde{B},\tilde{C},\tilde{D})$ is detectable.

###### Proof.

We will show that in each of the cases (i)–(iii), given the corresponding
constraint set  $\mathbf{C}$, a linear comparison system for the observer
state error dynamics (52) can be computed in the following form:

$\displaystyle
e^{x}_{k+1}\leq(\tilde{A}-L\tilde{C})e^{x}_{k}+(\tilde{B}-L\tilde{D})\tilde{w},$
(59)

with
$\tilde{w}\triangleq\begin{bmatrix}\delta^{v\top}&\delta^{w\top}&\delta^{\epsilon\top}\end{bmatrix}^{\top}$,
where the detectability of the pair ($\tilde{A},\tilde{C}$) is a necessary
condition for stabilizing the comparison system. If this can be shown, then
using the results in [23, Section 9.2.3], the solution $(P^{*},\Gamma^{*})$ to
the program in (58) returns the optimal observer gain
$L^{*}=(P^{*})^{-1}\Gamma^{*}$ for the linear comparison system (59), and
hence, for the original error dynamics (52) in the  minimum
$\mathcal{H}_{\infty}$ gain sense with an $\mathcal{H}_{\infty}$ gain of
$\eta^{*}$ (cf. Definition 6). This implies that the above linear comparison
system (59) satisfies the following asymptotic gain (AG) property [24]:

$\displaystyle\small\limsup_{k\to\infty}\|e^{x}_{k}\|_{2}\leq\alpha(\|{\tilde{\delta}}\|_{\ell_{\infty}}),\,\forall
e^{x}_{0},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
0\leq{\tilde{\delta}\leq[(\delta^{w})^{\top}\,(\delta^{v})^{\top}\,(\delta^{\epsilon})^{\top}]^{\top}},}\normalsize$

where $\tilde{\delta}$ is any realization of the augmented noise and outer-
approximation error interval width and $\alpha$ is any class
$\mathcal{K}_{\infty}$ function that is lower bounded by
$\eta^{*}\tilde{\delta}$. On the other hand, by setting $\delta=0$, the LMIs
in (58) reduce to their noiseless counterparts in [15, Eq. (19)]. Hence, by
[15, Theorem 2], the comparison system (59) is 0-stable (0-GAS), which in
addition to the AG property  above is equivalent to the ISS property for (59)
by [24, Theorem 1-e]. Thus, the designed observer is also ISS. So, what
remains to complete the proof is to show that the comparison system (59) can
indeed be computed in each of the cases as follows.

Case (i). Consider the nonlinear comparison system in (52). By satisfying the
constraint set $\mathbf{C}$, we enforce $-P$ to be Metzler, as well as
$P\tilde{A}-\Gamma\tilde{C},PV_{a}-\Gamma V_{b}$ and $PV_{a}-\Gamma V_{b}$ to
be non-negative. Also, $\Gamma$ is non-negative by assumption. Consequently,
since $P$ is positive definite, it becomes a non-singular M-matrix, i.e., a
square matrix whose negation is Metzler and whose eigenvalues have non-
negative real parts, and hence is inverse-positive [25, Theorem 1], i.e.,
$P^{-1}\geq 0$. Therefore, $L=P^{-1}\Gamma\geq 0$, $A-LC_{2}=P^{-1}(PA-\Gamma
C_{2})\geq 0$, $V_{a}-LV_{b}=P^{-1}(PV_{a}-\Gamma V_{b})\geq 0$ and
$D_{a}-LD_{b}=P^{-1}(PD_{a}-\Gamma D_{b})\geq 0$, because they are matrix
products of non-negative matrices.  So,
$|L|=L,|A-LC_{2}|=A-LC_{2},|V_{a}-LV_{b}|=V_{a}-LV_{b}$ and
$|D_{a}-LD_{b}|=D_{a}-LD_{b}$,  which turns (52) into the form of (59).

Case (ii). By applying the triangle inequality, the comparison system in (52)
can get upper bounded again as

$\displaystyle\begin{array}[]{rll}e^{x}_{k+1}&\leq(|A|+|LC_{2}|+\overline{F}_{\rho}+|L|\overline{F}_{\psi_{2}})e^{x}_{k}+|\hat{W}|\delta^{w}\\\
&+(|V_{a}|+|LV_{b}|-|LC_{2}||\Lambda NV_{2}|+(I-|A|)|\Lambda
NV_{2}|)\delta^{v}\\\
&+((I-|A|)|\Lambda|+|D_{a}|+|LD_{b}|-|LC_{2}||\Lambda|)\delta^{\epsilon}.\end{array}$
(63)

By a similar argument as in Case (i), enforcing $-P$ to be  Metzler along with
the constraints set $\mathbf{C}$ results in $|LC_{2}|=LC_{2},|LV_{b}|=LV_{b}$
and $|LD_{b}|=LD_{b}$,  and hence turns (63) into the form of (59).

Case (iii). Note that by the triangle inequality,
$|V_{a}-LV_{b}|=|(A-LC_{2})\Lambda
NV_{2}+LV_{2}+\Lambda(I-NC_{2})G_{1}SV_{1}|\leq|(A-LC_{2})||\Lambda
NV_{2}|+|L||V_{2}|+|\Lambda(I-NC_{2})G_{1}SV_{1}|$, and
$|D_{a}-LD_{b}|=|(A-LC_{2})\Lambda|\leq|(A-LC_{2})||\Lambda|$. These two
combined with (52) yield

$\displaystyle\small\begin{array}[]{rl}e^{x}_{k+1}&\leq(|A-LC_{2}|+\overline{F}_{\rho}+|L|\overline{F}_{\psi_{2}})e^{x}_{k}+|\hat{W}|\delta^{w}+|\Lambda|\delta^{\epsilon}\\\
&+(|L||V_{2}|+|\Lambda
NV_{2}|+|\Lambda(I-NC_{2})G_{1}SV_{1}|)\delta^{v}.\end{array}\normalsize$ (66)

The rest of the proof is to enforce that $A-LC_{2}$ and $L$ are non-negative
to turn (66) into the form of (59), which  is similar to the the proofs of the
previous two cases. ∎

##  V Illustrative Example

We now illustrate the effectiveness of our proposed  resilient observer using
a three-area  power system  [2, Figure 1], where each control area consists of
a generator and load buses with transmission lines between areas. The
nonlinear continuous-time model of the buses is slightly modified based on
[26], with  the subscript $i$ being the bus number:

$\displaystyle\begin{array}[]{rl}\dot{f}_{1}(t)&=-\frac{1}{m_{1}}(\phi_{i}(t)-(P_{M_{1}}(t)+d_{1}(t)))+w_{2,1}(t),\\\
\dot{f}_{i}(t)&=-\frac{1}{m_{i}}(\phi_{i}(t)-P_{M_{i}}(t))+w_{2,i}(t),\
i\in\\{2,3\\},\\\
\dot{\theta}_{i}(t)&=f_{i}(t)+w_{1,i}(t),\quad\quad\quad\quad\quad\quad\
i\in\\{1,2,3\\},\end{array}\normalsize$

with $\phi_{i}(t)\triangleq D_{i}f_{i}(t)+\sum_{l\in
S_{i}}P_{il}(t)+P_{L_{i}}(t))$, where $\theta_{i}$ is the phase angle, $f_{i}$
is the angular frequency, $m_{i}=0.01$, $D_{i}=0.11$, $P_{M_{i}}(t)$ is the
mechanical power (the control input), $P_{L_{i}}(t)$ is a known power demand,
$S_{i}$ is the set of neighboring buses of $i$, and  the nonlinear tie line
power flow equation is as follows:
$P_{il}(t)=-P_{li}(t)=t_{il}\sin(\theta_{i}(t)-\theta_{l}(t))$, with
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}t_{il}=1}$.
Only the actuator of Control Area 1 is attacked and the false data injection
signal is $d_{1}(t)$.

On the other hand, the output equation is given as follows:

$\displaystyle\small y_{i}(t)$ $\displaystyle=[\theta_{i}(t)\
f_{i}(t)]^{\top}+v_{i}(t),\ i\in\\{1,3\\},$ $\displaystyle y_{2}(t)$
$\displaystyle=[\theta_{2}(t)\ f_{2}(t)]^{\top}+d_{2}(t)+v_{2}(t),\normalsize$

where only the sensor $y_{2}(t)$ is injected with a false data signal
$d_{2}(t)$. Thus, the concatenated attack/unknown input signal is
$d(t)=[d_{1}(t)\ d_{2}(t)]^{\top}$ and the $G$ and $H$ matrices in (7)
corresponding to the attack locations are given by
$G=\begin{bmatrix}0&0&0&0&0&0\\\ 0&1&0&0&0&0\end{bmatrix}^{\top}$ and
$H=\begin{bmatrix}0&0&0&0&0&0\\\ 0&0&0&1&0&0\end{bmatrix}^{\top}$.

In our simulations,  the _forward Euler_ method is used to discretize the
system dynamics with a sampling time $dt=0.01s$ and both $P_{M_{i}}(t)$ and
$P_{L_{i}}(t)$ were set to be identically zero. Moreover, for $i=1,\dots,3$,
the process noise $w_{i}(t)$ and the measurement noise $v_{i}(t)$ were assumed
to be bounded within the bounds
$\left[\begin{bmatrix}-50&-50\end{bmatrix}^{\top},\begin{bmatrix}50&50\end{bmatrix}^{\top}\right]$
and
$\left[\begin{bmatrix}-0.5&-0.5\end{bmatrix}^{\top},\begin{bmatrix}0.5&0.5\end{bmatrix}^{\top}\right]$,
respectively.

For the sake of comparison, we first applied our previous input and state
observer [14] that does not have stabilizing gains to the above example, which
we found to not be able to yield stable interval estimates (i.e., the framer
interval width diverges). On the other hand,  when implementing the proposed
observer in (27)–(35), the optimization problem in (58) was solved with the
additional linear constraints in Case (iii), and we obtained the following
observer gain:

$\displaystyle L=\begin{bmatrix}0.70&0&0.27&0&0\\\ 0&0&0.38&0&0\\\
0&0.83&73.19&0&0\\\ -0.0022&0.0084&174.55&0.0056&-0.0001\\\
0&0&0.14&0.70&0.005\\\ 0.0050&0.0098&0.11&0.01&0.62\end{bmatrix}.$

As shown in Figures 1 and 2, all the states and attack signals are bounded by
the framers computed by the proposed observer, demonstrating its correctness
and ability to obtain resilient state estimates and to reconstruct attack
signals. Finally,  as shown in Figure 3, the actual state and input estimation
error sequences  (i.e., the framer interval widths) converge to steady-state
values,  demonstrating the input-to-state stability of the  proposed interval
observer.

Figure 1: States: $\theta_{i},f_{i}$, and their upper and lower framers
$\overline{\theta}_{i},\underline{\theta}_{1},\overline{f}_{i},\underline{f}_{i}$,
returned by the proposed approach.

Figure 2: Attack signals: $d_{1},d_{2}$, and their upper and lower framers
$\overline{d}_{1},\underline{d}_{1},\overline{d}_{2},\underline{d}_{2}$,
returned by the proposed approach.

Figure 3: State and input estimation error sequences.

##  VI Conclusion

In this paper, the problem of  resilient state estimation and attack
reconstruction for nonlinear discrete-time systems with nonlinear
observations/constraints, that are subject to bounded noise signals, was
addressed. In the considered setting, both sensors and actuators could be
affected by attack signals/unknown inputs. By introducing auxiliary states, as
well as taking advantage of mixed-monotone decomposition of nonlinear
functions and affine  parallel outer-approximation of the observation
functions, the proposed observer was shown to be correct, i.e., it recursively
computes interval estimates that by construction, contain the true states and
unknown inputs of the system. Further, several semi-definite programs were
provided to synthesize the proposed observer gains  that guarantee input-to-
state stability  of the observer and  optimality of the  proposed interval
observer design.  Future work will include alternative designs for minimizing
$L_{1}$ gain, similar to [27], as well as an extension to continuous-time
nonlinear systems and hybrid systems.

## References

  * [1] W. Liu and I. Hwang, “Robust estimation and fault detection and isolation algorithms for stochastic linear hybrid systems with unknown fault input,” _IET Control Theory & Applications_, vol. 5, no. 12, pp. 1353–1368, 2011\.
  * [2] S. Z. Yong, M. Zhu, and E. Frazzoli, “Switching and data injection attacks on stochastic cyber-physical systems: Modeling, resilient estimation and attack mitigation,” _ACM Transactions on Cyber-Physical Systems_ , vol. 2, no. 2, p. 9, 2018.
  * [3] S. Z. Yong, “Simultaneous input and state set-valued observers with applications to attack-resilient estimation,” in _American Control Conference (ACC)_. IEEE, 2018, pp. 5167–5174.
  * [4] G. Chen, Y. Zhang, S. Gu, and W. Hu, “Resilient state estimation and control of cyber-physical systems against false data injection attacks on both actuator and sensors,” _IEEE Transactions on Control of Network Systems_ , vol. 9, no. 1, pp. 500–510, 2021.
  * [5] C. Wu, Z. Hu, J. Liu, and L. Wu, “Secure estimation for cyber-physical systems via sliding mode,” _IEEE Transactions on Cybernetics_ , vol. 48, no. 12, pp. 3420–3431, 2018.
  * [6] E. Mousavinejad, F. Yang, Q. Han, and L. Vlacic, “A novel cyber attack detection method in networked control systems,” _IEEE Transactions on Cybernetics_ , vol. 48, no. 11, pp. 3254–3264, 2018.
  * [7] M. Corradini and A. Cristofaro, “Robust detection and reconstruction of state and sensor attacks for cyber-physical systems using sliding modes,” _IET Control Theory & Applications_, vol. 11, no. 11, pp. 1756–1766, 2017\.
  * [8] A. Lu and G. Yang, “Secure state estimation for cyber-physical systems under sparse sensor attacks via a switched Luenberger observer,” _Information sciences_ , vol. 417, pp. 454–464, 2017.
  * [9] N. Ellero, D. Gucik-Derigny, and D. Henry, “An unknown input interval observer for LPV systems under ${L}_{2}$-gain and ${L}_{\infty}$-gain criteria,” _Automatica_ , vol. 103, pp. 294–301, 2019.
  * [10] X. Wang, C. Tan, L. Liu, and Q. Qi, “A novel unknown input interval observer for systems not satisfying relative degree condition,” _International Journal of Robust and Nonlinear Control_ , vol. 31, no. 7, pp. 2762–2782, 2021\.
  * [11] G. Marouani, T. Dinh, T. Raissi, X. Wang, and H. Messaoud, “Unknown input interval observers for discrete-time linear switched systems,” _European Journal of Control_ , vol. 59, pp. 165–174, 2021.
  * [12] L. Wei and W. Yang, “Hybrid observer design for nonlinear systems with unknown inputs,” in _2021 China Automation Congress (CAC)_. IEEE, 2021, pp. 4094–4099.
  * [13] M. Khajenejad and S. Z. Yong, “Simultaneous input and state interval observers for nonlinear systems with full-rank direct feedthrough,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 5443–5448.
  * [14] ——, “Simultaneous input and state interval observers for nonlinear systems with rank-deficient direct feedthrough,” in _2021 European Control Conference (ECC)_. IEEE, 2021, pp. 2311–2316.
  * [15] M. Khajenejad, F. Shoaib, and S. Z. Yong, “Interval observer synthesis for locally Lipschitz nonlinear dynamical systems via mixed-monotone decompositions,” in _2022 American Control Conference (ACC)_. IEEE, 2022, pp. 2970–2975.
  * [16] M. Khajenejad and S. Z. Yong, “$\mathcal{H}_{\infty}$-optimal interval observer synthesis for uncertain nonlinear dynamical systems via mixed-monotone decompositions,” _IEEE Control Systems Letters_ , vol. 6, pp. 3008–3013, 2022.
  * [17] K. Singh, Q. Shen, and S. Z. Yong, “Mesh-based affine abstraction of nonlinear systems with tighter bounds,” in _Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 3056–3061.
  * [18] M. Abate, M. Dutreix, and S. Coogan, “Tight decomposition functions for continuous-time mixed-monotone systems with disturbances,” _IEEE Control Systems Letters_ , vol. 5, no. 1, pp. 139–144, 2020.
  * [19] L. Yang, O. Mickelin, and N. Ozay, “On sufficient conditions for mixed monotonicity,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 12, pp. 5080–5085, 2019.
  * [20] M. Khajenejad and S. Z. Yong, “Tight remainder-form decomposition functions with applications to constrained reachability and guaranteed state estimation,” _IEEE Transactions of Automatic Control_ , pp. 1–16, 2023, early access.
  * [21] L. Yang and N. Ozay, “Tight decomposition functions for mixed monotonicity,” in _Conference on Decision and Control (CDC)_ , 2019, pp. 5318–5322.
  * [22] D. Efimov, T. Raïssi, S. Chebotarev, and A. Zolghadri, “Interval state observer for nonlinear time varying systems,” _Automatica_ , vol. 49, no. 1, pp. 200–205, 2013.
  * [23] G. Duan and H. Yu, _LMIs in control systems: analysis, design and applications_. CRC press, 2013.
  * [24] E. Sontag and Y. Wang, “New characterizations of input-to-state stability,” _IEEE Trans. on Automatic Control_ , vol. 41, no. 9, pp. 1283–1294, 1996\.
  * [25] R. Plemmons, “M-matrix characterizations. I: nonsingular M-matrices,” _Linear Algebra and its Applications_ , vol. 18, no. 2, pp. 175–188, 1977\.
  * [26] H. Kim, P. Guo, M. Zhu, and P. Liu, “Attack-resilient estimation of switched nonlinear cyber-physical systems,” in _2017 American Control Conference (ACC)_. IEEE, 2017, pp. 4328–4333.
  * [27] T. Pati, M. Khajenejad, S. P. Daddala, and S. Z. Yong, “$L_{1}$-robust interval observer design for uncertain nonlinear dynamical systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 3475–3480, 2022.


#  From creep to flow: Granular materials under cyclic shear

Ye Yuan1 Zhikun Zeng1 Yi Xing1 Houfei Yuan1 Shuyang Zhang1 Walter Kob2
<EMAIL_ADDRESS>Yujie Wang1,3,4<EMAIL_ADDRESS>1 School of
Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China
2 Laboratoire Charles Coulomb, University of Montpellier and CNRS, 34095
Montpellier, France
3 State Key Laboratory of Geohazard Prevention and Geoenvironment Protection,
Chengdu University of Technology, Chengdu 610059, China
4 Department of Physics, College of Mathematics and Physics, Chengdu
University of Technology, Chengdu 610059, China

###### Abstract

Granular materials such as sand, powders, and grains are omnipresent in daily
life, industrial applications, and earth-science [1]. When unperturbed, they
form stable structures that resemble the ones of other amorphous solids like
metallic and colloidal glasses [2]. It is commonly conjectured that all these
amorphous materials show a universal mechanical response when sheared slowly,
i.e., to have an elastic regime, followed by yielding [3]. Here we use X-ray
tomography to determine the microscopic dynamics of a cyclically sheared
granular system in three dimensions. Independent of the shear amplitude
$\Gamma$, the sample shows a cross-over from creep to diffusive dynamics,
indicating that granular materials have no elastic response and always yield,
in stark contrast to other glasses. The overlap function [4] reveals that at
large $\Gamma$ yielding is a simple cross-over phenomenon, while for small
$\Gamma$ it shows features of a first order transition with a critical point
at $\Gamma\approx 0.1$ at which one finds a pronounced slowing down and
dynamical heterogeneity. Our findings are directly related to the surface
roughness of granular particles which induces a micro-corrugation to the
potential energy landscape, thus creating relaxation channels that are absent
in simple glasses. These processes must be understood for reaching an
understanding of the complex relaxation dynamics of granular systems.

Introduction

Yielding of amorphous materials is ubiquitous, governing a wide range of
phenomena like mechanical failure of metallic glasses [5], complex rheologies
of soft glasses [6], or geophysical catastrophes [7]. A typical amorphous
solid under quasistatic simple shear will deform elastically at small strain
and flow plastically beyond the yielding strain. Depending on whether the
material is brittle or ductile, the stress-strain curve shows a drop or a
smooth cross-over to a plateau, respectively. On the particle level, yielding
is a cooperative phenomenon of local plastic events [8, 9], and computer
simulations hint that this phenomenon shares many aspects with a first order
phase transition [4, 10, 11]. However, for real systems the precise nature of
this transition has not yet been clarified since such experiments are very
challenging.

In contrast to standard glasses, for granular systems there is at present no
satisfactory understanding of yielding on the level of the particles. It has
been suggested that sheared granular solids are marginally stable, i.e.,
contacts between particles will change irreversibly even under a tiny applied
strain, which implies that such systems have no elastic regime [12]. This view
seems to be at odds with experimental results which found that the stress-
strain curve of granular solids under simple shear does in fact resemble the
one for a glass [13]. One possibility to resolve this discrepancy is to
consider a different type of driving, i.e., cyclic shear, since it allows to
probe directly the reversibility of the particle trajectories and hence the
presence of elasticity in the system. A number of computational studies have
in fact used this setup to investigate the particle dynamics in soft-sphere
jammed states and model glasses [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24].
These works, as well as related experiments on soft glasses [25, 26, 27], have
indeed revealed a reversible-irreversible transition as the cyclic shear
amplitude $\Gamma$ increases, thus supporting the existence of an elastic
behavior if $\Gamma\lesssim 0.1$. However, these findings are in contrast with
experimental results for cyclically sheared granular particles in three
dimensions (3D) which indicate the absence of an elastic regime, since the
particle trajectories are found to be irreversible, at least for the values of
$\Gamma$ considered, i.e., $0.07\leq\Gamma\leq 0.26$ [28, 29].

Although a few experimental studies have probed the microscopic dynamics of
driven 3D granular materials [30, 31, 32], none of them has made an explicit
connection between this dynamics and the phenomenon of yielding. Establishing
this link is important not only for understanding the plasticity of amorphous
materials but also for developing reliable granular constitution laws [33].
The goal of the present work is therefore to identify this connection.

Cyclic shear experiment

In this work, we experimentally investigate a cyclically sheared granular
system using a X-ray tomography technique [28, 29, 34, 35]. The system
contains $\approx$14000 $50:50$ bidisperse plastic spherical beads of
diameters 5 mm and 6 mm and no crystallization was detected. Particles are
placed in a shear cell and cyclically sheared with an amplitude $\Gamma$ [28,
36]. The size of the shear cell at shear $\gamma=0$ is $24d\times 24d\times
24d$, where $d$ is the diameter of the small beads, and which in the following
will be taken as the unit length. The shear rate is small so that we are in
quasistatic conditions with an inertia number less than $10^{-3}$ [37].

Beads are initially placed in the cell, forming a reproducible loose packings,
and then compacted by cyclic shear until the steady state is reached (see
Extended Data Fig. 1). X-ray tomography scans are taken at $\gamma=0$, with a
periodicity between 1 and 10 cycles, from which we extract the microscopic
structure of the system and the dynamics of the particles. To improve the
statistics of the results, we use $3\sim 5$ independent realizations for each
$\Gamma$, and to mitigate finite size effects we exclude the particles located
closer than $6d$ from the boundaries.

Steady-state dynamics

For small $\Gamma$ the dynamics of the particles is basically isotropic
although weak convection caused by gravity emerges for large $\Gamma$ (see
Extended Data Fig. 2). Hence we focus here on the dynamics in the horizontal
$x$-and $y$-directions. Figure 1(a) presents the mean squared displacement
(MSD), $\langle\delta h^{2}(\Delta n)\rangle$, as a function of the number of
cycles $\Delta n$, where $\delta h$ is the horizontal displacement and
$\langle.\rangle$ stands for the average over different particles, starting
configurations, and realizations. At large $\Delta n$ we find the expected
diffusive growth while at small $\Delta n$ the MSD shows a power-law
dependence, with a $\Gamma$-independent exponent $\approx$ 0.65. Such a
universal creep dynamics at small $\Delta n$ demonstrates that our system has
no caging or elastic regime, in agreement with earlier findings [38, 28, 29].
This is also confirmed by the absence of a two-step relaxation in the self-
intermediate scattering function (see Extended Data Fig. 3).

Figure 1:  Bulk dynamics of cyclic shear for different shear amplitudes
$\Gamma$. (a) Mean squared displacement in horizontal directions vs. shear
cycle number $\Delta n$. A universal crossover from sub-diffusion MSD
$\propto\Delta n^{0.65}$ to normal diffusion (dashed lines) is observed,
corresponding to the $\Gamma$-dependent yielding point. (b) Schematics of how
a granular system explores its potential energy landscape. Due to the particle
surface roughness, the PEL has not only a structure on the particle scale $d$,
but also a micro-corrugation. Starting in the left metabasin of the PEL, the
double arrows indicate the back and forth motions of the system during a shear
cycle. Depending on $\Gamma$, the system does/(does not) leave the MB during a
cycle (blue and orange arrows, respectively). Yielding means that the system
has overcome permanently the PEL barrier (on the scale $d$). (c)
$\Gamma$-dependence of the diffusion coefficient $D$ (circles), the yielding
strain $\Delta\gamma_{c}$ (diamonds), and the strain $\Delta\gamma_{M}$ at
which memory is lost (triangles). Note that $\Delta\gamma_{c}\approx
0.4\Delta\gamma_{M}$ for all $\Gamma$. Error bars represent the standard
deviations from $3-5$ independent realizations. (d) Memory $M$, defined in Eq.
(1), as a function of $\Delta\gamma$, from which $\Delta\gamma_{M}$ is
obtained by $M(\Delta\gamma_{M})=0$. For $\Gamma\leq 0.175$,
$\Delta\gamma_{M}$ is estimated from a linear extrapolation of
$M(\Delta\gamma)$. Color codes are the same as in Fig. 1(a).

The absence of caging is related to the fact that granular particles have a
rough surface which engenders to the potential energy landscape (PEL) of the
system a micro-corrugation on a length scale that is much smaller than $d$,
absent in atomic systems or pure hard-spheres, and sketched in Fig. 1(b). This
corrugation permits the particles to accumulate an irreversible displacement
even at the smallest $\Gamma$ hence preventing caging. Starting from a local
minimum in the PEL, the creep dynamics corresponds thus to the exploration of
the local metabasin (MB), i.e., the set of configurations which are on the
particle level very similar to the initial one. This motion, affected strongly
by memory effects (see below), allows the system to slowly reach the boundary
of the MB, i.e., the yielding point, at which many particles will have changed
their neighbors. After yielding the memory is lost and the MSD becomes
diffusive.

In the following we report the dynamics as a function of the accumulated
strain, i.e., $\Delta\gamma=4\Delta n\Gamma$, which takes into account one
part of the expected $\Gamma$-dependence of the dynamics. Figure 1(c) shows
that the diffusion coefficient $D$, obtained from the Einstein relation
$\langle\delta h^{2}\rangle=2D\Delta\gamma$, is basically constant for
$\Gamma\lesssim 0.1$ and starts to grow sharply beyond this threshold. Also
the yielding strain $\Delta\gamma_{c}$, locating the crossover from sub-
diffusive to diffusive regime (see Extended Data Fig. 4), displays a maximum
at $\Gamma=0.1$. These observations indicate a change in the underlying
dynamics as a function of $\Gamma$, directly linked to the way how the system
(i) explores its MB and (ii) yields at large $\Delta n$.

Figure 1(b) illustrates that for $\Gamma\lesssim 0.1$ the system will stay
inside the MB during many cycles before it yields. Since the approach to the
MB boundary is slow if $\Gamma\lesssim 0.1$, crossing the boundary will
involve many particles at a time, i.e., this yielding is a highly collective
process. In contrast to this, if $\Gamma\gtrsim 0.1$ the particles can leave
the MB during each cycle, notably at maximum strain, but due to friction and
particle roughness their trajectories are largely reversible. Since this
reversibility is not perfect, some particles will have left the original MB
and thus the MB is slowly drained during the cycling. In this case the
yielding, occurring once many particles have switched their neighbors, is no
longer a cooperative process since it occurs very gradually. In contrast to
this, simulations for model glasses [20] or soft-sphere jammed states [15,
18], in which the PEL has no micro-corrugation, found that $D$ is only non-
zero above a certain threshold in $\Gamma$, thus a $\Gamma$-dependence that is
very different from the one found here.

Figure 2: Yielding as a phase transition. (a-c) Overlap function distribution
$P(Q)$ as a function of $\Delta\gamma/\Delta\gamma_{c}$ for $\Gamma=0.05$,
$0.1$, and $0.2$. $P(Q)$ shifts to smaller values as
$\Delta\gamma/\Delta\gamma_{c}$ grows (see legends). (d) $P(Q)$ at the
yielding point $\Delta\gamma=\Delta\gamma_{c}$ shows the presence of two
master curves. Inset: The standard deviation of $P(Q)$ versus $\Gamma$ has a
sharp transition at $\Gamma\approx 0.1$. (e) $\langle Q\rangle$ as a function
of $\Delta\gamma-\Delta\gamma_{c}$. The decay is slowest for $\Gamma\approx
0.1$ indicating the presence of a critical slowing down close to a critical
point. Inset: Also $\langle Q\rangle$ versus $\Delta\gamma$ shows a slowing
down at $\Gamma\approx 0.1$. (f) Dynamic phase diagram. For any $\Gamma$, the
system evolves from a creep (sub-diffusion) to flow (diffusion) regime as
$\Delta\gamma$ grows and yields at $\Delta\gamma=\Delta\gamma_{c}$. This
yielding is first-order-like for $\Gamma<0.1$ (purple solid curve) but shows
only a continuous crossover for $\Gamma>0.1$ (red dashed curve). $\Gamma=0.1$
corresponds to the critical point (star). In (d) and (e) the color codes are
the same as in Fig. 1(a).

The mentioned reversibility of the particle motion can be inferred directly
from the smallness of the MSD at $\Delta n=1$ in Fig. 1(a), indicating that
the dynamics is not Markovian but has instead a significant memory [39]. To
quantify this memory we consider the correlation function between the
displacements of particle $i$ at two consecutive intervals of length
$\Delta\gamma$:

$M(\Delta\gamma)=-\langle\delta x_{i}(2\Delta\gamma)\delta
x_{i}(\Delta\gamma)\rangle/\langle\delta x_{i}(\Delta\gamma)^{2}\rangle,$ (1)

where $\delta x_{i}(t)=x_{i}(t)-x_{i}(t-\Delta\gamma)$. Figure 1(d)
demonstrates that $M$ is significantly positive at small $\Delta\gamma$, i.e.,
the displacements are anti-correlated, and $M$ decays with a
$\Gamma$-dependent rate. For large $\Delta\gamma$ the memory vanishes, i.e.,
after yielding the dynamics becomes Markovian. We determine the strain
$\Delta\gamma_{M}$ at which the memory is lost, $M(\Delta\gamma_{M})=0$, and
find that $\Delta\gamma_{M}(\Gamma)$ tracks $\Delta\gamma_{c}(\Gamma)$ very
well, i.e., $\Delta\gamma_{c}\approx 0.4\Delta\gamma_{M}$, as shown in Fig.
1(c). Hence the creep dynamics is accompanied by strong memory and once the
initial memory has been halved the system yields.

Yielding as a phase transition

The two distinct regimes in the yielding dynamics indicate that the nature of
the corresponding dynamic phase transition, studied in simulations with simple
shear [4, 11], changes with $\Gamma$. One standard approach to probe the
properties of a phase transition is to use as an order parameter the overlap
$Q(\Delta\gamma)\in[0,1]$, which measures the similarity of two configurations
separated by $\Delta\gamma$ (see Methods). One expects that $\langle Q\rangle$
decreases with increasing $\Delta\gamma$ and the distribution $P(Q)$ allows to
identify the nature of the phase transition. Figures 2(a-c) show that $P(Q)$
are close to $1$ for small $\Delta\gamma/\Delta\gamma_{c}$, i.e., most
particles have not yet moved significantly. With increasing
$\Delta\gamma/\Delta\gamma_{c}$, $P(Q)$ shifts to the left before converging
to the random distribution. For all $\Delta\gamma$ considered, the width of
$P(Q)$ is significantly larger for $\Gamma\leq 0.1$ than for $\Gamma=0.2$.
This is quantified in Fig. 2(d) by plotting $P(Q)$ at
$\Delta\gamma=\Delta\gamma_{c}$ for different $\Gamma$, resulting in two
master curves: A wide one for $\Gamma\leq 0.1$ and a narrow one for
$\Gamma>0.1$. Such a $\Gamma$-dependence is supported by the inset of Fig.
2(d), which demonstrates that the standard deviation of $P(Q)$ strongly drops
at $\Gamma\approx 0.1$.

Figure 3:  Self-part of the Van Hove function $G_{s}$ for different $\Gamma$.
The particle displacement is expressed in terms of the normalized distance
$d_{x}=|\delta x|/\sqrt{\langle\delta x^{2}\rangle}$. The color codes are the
same as in Fig. 1(a). Panels (a)-(c) correspond to
$\Delta\gamma=0.25\Delta\gamma_{c}$, $\Delta\gamma_{c}(\Gamma)$, and
$3\Delta\gamma_{c}(\Gamma)$, respectively. Solid and dashed curves are,
respectively, the Gumbel [Eq. (2)] and Gaussian laws, showing that $G_{s}$ is
non-Gaussian even at small $d_{x}$. Insets show the ratio between $G_{s}$ and
the Gumbel law $G^{g}$. The vertical dashed lines at $d_{x}=2$ mark the
distance beyond which $G_{s}$ is no longer described well by the Gumbel law.
For $d_{x}\leq 2$ the Gumbel law is universal for all $\Gamma$. (d) Intensity
of the excess exponential tail obtained by averaging $G_{s}/G^{g}$ in the
interval $d_{x}\in[3.5,4.5]$. These curves show a maximum at
$\Delta\gamma\approx\Delta\gamma_{c}$, i.e., at the yielding point $G_{s}$ has
the most pronounced tail.

While the $P(Q)$ documented in Ref. [4] showed the double peak structure
signalling a first order transition in which the two phases of the system can
coexist, here we do not see clear evidence for this. However, it is well known
that finite size effects do smear out such a structure, leaving only a single
broad peak [4], which might explain the absence of a double peak due to the
moderate particle number that we consider. Alternatively one can argue that
the free energy barrier between the two phases is small since the micro-
corrugation of the PEL allows the existence of many locally stable particle
configurations that cannot be clearly assigned to one of the two phases. In
other words, the surface tension between the two phases is small, and hence
the first order transition is weak. A weak transition hints the vicinity of a
critical point at which the transition becomes second order and the dynamics
shows a critical slowing down. Upon a further increase of the external
parameter (here $\Gamma$), the transition ceases to exist. This scenario is
indeed compatible with our data in that the total strain $\Delta\gamma_{c}$ at
which the yielding occurs shows a maximum, Fig. 1(c), and also the
$\Delta\gamma$-dependence of the order parameter, $\langle
Q\rangle(\Delta\gamma)$, shows a non-monotonic dependence on $\Gamma$, Fig.
2(e). (Also here the divergence expected at the critical point is rounded off
because the dynamics blurs the two phases.) For $\Gamma\geq 0.1$ the
$\Delta\gamma$-dependence of $P(Q)$ indicates that there is no longer a phase
transition but just a smooth crossover between the two phases, Fig. 2(c).
Figure 2(f) summarizes the different behaviors as a function of $\Delta\gamma$
and $\Gamma$.

Dynamical heterogeneity

To obtain a microscopic understanding of the yielding, we investigate the
dynamics on the particle level. Figures 3(a-c) show the distributions of the
particle displacements in the $x$-or $y$-direction, i.e., the self-part of the
Van Hove function $G_{s}(\delta x,\Delta\gamma$) [2], for
$\Delta\gamma=0.25\Delta\gamma_{c}$, $\Delta\gamma_{c}$, and
$3\Delta\gamma_{c}$. Plotting $G_{s}$ for different $\Gamma$ as a function of
the rescaled distance $d_{x}=|\delta x|/\sqrt{\langle\delta x^{2}\rangle}$
roughly results in a master curve. In contrast to thermal systems the shape of
this curve is clearly non-Gaussian even for small $\Delta\gamma$ (dotted
lines), demonstrating the presence of particle motion that is faster than
expected for a Gaussian process.

For $d_{x}<2$ the distribution can be fitted well by a Gumbel law [28]

$G^{g}(d_{x})=B(\lambda)\exp[-d_{x}/\lambda-\exp(-d_{x}/\lambda)]$ (2)

(solid lines), where $\lambda=0.57$ characterizes the shape of the Gumbel
distribution and $B(\lambda)=2.72$ is a normalization constant. If $d_{x}>2$,
$G_{s}$ exceeds the Gumbel law and shows an exponential decay (see insets).
The deviation of this excess tail from the Gumbel law is quantified in Fig.
3(d), which presents the ratio $G_{s}(d_{x})/G^{g}(d_{x})$ averaged over
$d_{x}\in[3.5,4.5]$. This deviation attains a maximum at
$\Delta\gamma\approx\Delta\gamma_{c}$ for all $\Gamma$, signaling that at
yielding the distribution is widest and hence one has maximal dynamical
heterogeneity. The largest excess is found for $\Gamma\approx 0.1$, which is
in line with theoretical arguments that non-Markovian processes with a
significant memory give rise to a pronounced tail in the Van Hove function
[40].

Figure 4: Dynamical heterogeneity as a function of $\Delta\gamma$ and
$\Gamma$. (a) Non-Gaussian parameter $\alpha_{2}$ as a function of
$\Delta\gamma$ for different $\Gamma$. Solid curves are exponential fits
$\alpha_{2}=A\cdot\exp(-\Delta\gamma/\Delta\gamma_{g})$ at large
$\Delta\gamma$. (b) The $\Gamma$-dependence of $\Delta\gamma_{g}$ tracks the
one of $\Delta\gamma_{c}(\Gamma)$, i.e., $\Delta\gamma_{c}\approx
0.15\Delta\gamma_{g}$. Inset shows that $A(\Gamma)$ mildly peaks at
$\Gamma\approx 0.1$. (c) Fraction of particles involved in the largest
connected cluster (two particles are defined as connected if the center
distance is smaller than $1.2$ times their average diameter), among the top
$10\%$ mobile particles as a function of $\Delta\gamma/\Delta\gamma_{c}$. In
(a) and (c) the color codes are the same as in Fig. 1(a).

We quantify the evolution of the shape of $G_{s}$ via the non-Gaussian
parameter $\alpha_{2}=\langle\delta x^{4}\rangle/(3\langle\delta
x^{2}\rangle^{2})-1$. Figure 4(a) shows that at small $\Delta\gamma$,
$\alpha_{2}$ varies little and is significantly larger than the value expected
from the Gumbel law, $\approx 0.5$, in agreement with Fig. 3(d), i.e., the
exponential tail contributes considerably to $\alpha_{2}$. $\alpha_{2}$ starts
to drop quickly for $\Delta\gamma\gtrsim\Delta\gamma_{c}$, and can be
described well by an exponential,
$\alpha_{2}(\Delta\gamma)=A(\Gamma)\exp(-\Delta\gamma/\Delta\gamma_{g}(\Gamma))$,
where $\Delta\gamma_{g}$ is the strain scale for the recovery of Gaussian
dynamics and $A$ is the amplitude of $\alpha_{2}$. Figure 4(b) reveals that
$\Delta\gamma_{g}(\Gamma)$ perfectly tracks $\Delta\gamma_{c}(\Gamma)$,
similar to $\Delta\gamma_{M}$ in Fig. 1(c), and also $A$ peaks mildly at
$\Gamma\approx 0.1$ (Inset). Hence we conclude that the dynamical
heterogeneities are indeed largest at around $\Gamma=0.1$.

As a direct probe of dynamical heterogeneity, we determine the spatial
arrangement of the fastest particles (top $10\%$) by calculating the number of
particles belonging to the largest connected cluster (defined via a nearest
neighbor criterion). Figure 4(c) presents this number, normalized by the total
number of fast particles, $f$, and it peaks at
$\Delta\gamma/\Delta\gamma_{c}\approx 1$, i.e., yielding is accompanied by a
maximal cooperativity, in agreement with Fig. 3(d). A random choice of 10% of
the particles in the sample gives a $f\approx 0.04$, well below the values we
find here, demonstrating that the observed clustering is indeed significant.
(See Extended Data Fig. 5 for the cluster size distribution.) Consistent with
the phase diagram in Fig. 2(f), dynamical heterogeneity peaks at yielding and
$\Gamma\approx 0.1$ as displayed in Figs. 3 and 4.

Conclusions

Under simple shear conditions, yielding is associated with a drop in the
stress-strain curve, signalling the transition from elastic to plastic
behavior. For this type of driving, the mechanical response of granular
materials will be very similar [3] since the PEL’s micro-corrugation is
irrelevant. In contrast to this, the cyclic shear considered here permits to
investigate the highly non-trivial effect of this micro-corrugation and hence
to detect mechanical features of the system inaccessible in a simple shear
setup. In view of the minimal ingredients needed to generate the micro-
corrugation, we expect the surprising creep dynamics and the associated
yielding behavior reported here to be generic features of granular materials
and hence to be important in a multitude of situations, such as small tremors
in geo-sciences or aging of civil engineering structures.

The details of the yielding dynamics, i.e., the nature of the phase
transition, will depend on the micro-roughness, shape, as well as the friction
coefficient of the particles, since all these parameters influence the micro-
corrugation of the PEL and hence the system dynamics. It will thus be
important to study how these quantities affect the fractional diffusion
behavior of granular materials, or the strength of the memory. How this
dependence can be included in the present theoretical approaches is not clear
and thus remains a challenge for the future. We also mention that for simple
shear it is custom to classify the yielding as either ductile or brittle. Our
results show that for granular materials the nature of yielding depends
strongly on the driving protocol, i.e., simple shear vs. cyclic shear, and
also on the shear amplitude $\Gamma$. This dependence, which is absent in more
standard disordered materials, suggests that for granular materials such a
classification might not be possible. Advancing on these points will lead to a
fundamental understanding of granular rheology and thus permit to develop a
new holistic view of the failure of complex materials.

## References

  * Jaeger  _et al._ [1996] H. M. Jaeger, S. R. Nagel, and R. P. Behringer, Granular solids, liquids, and gases, Rev. Mod. Phys. 68, 1259 (1996).
  * Binder and Kob [2011] K. Binder and W. Kob,  _Glassy materials and disordered solids: An introduction to their statistical mechanics_ (World Scientific, 2011).
  * Nicolas  _et al._ [2018] A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Deformation and flow of amorphous solids: Insights from elastoplastic models, Rev. Mod. Phys. 90, 045006 (2018).
  * Jaiswal  _et al._ [2016] P. K. Jaiswal, I. Procaccia, C. Rainone, and M. Singh, Mechanical yield in amorphous solids: A first-order phase transition, Phys. Rev. Lett. 116, 085501 (2016).
  * Wang  _et al._ [2004] W.-H. Wang, C. Dong, and C. Shek, Bulk metallic glasses, Mater. Sci. Eng., R 44, 45 (2004).
  * Bonn  _et al._ [2017] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield stress materials in soft condensed matter, Rev. Mod. Phys. 89, 035005 (2017).
  * Ancey [2007] C. Ancey, Plasticity and geophysical flows: A review, J. Non-Newton. Fluid Mech. 142, 4 (2007).
  * Schall  _et al._ [2007] P. Schall, D. A. Weitz, and F. Spaepen, Structural rearrangements that govern flow in colloidal glasses, Science 318, 1895 (2007).
  * Ghosh  _et al._ [2017] A. Ghosh _et al._ , Direct observation of percolation in the yielding transition of colloidal glasses, Phys. Rev. Lett. 118, 148001 (2017).
  * Parisi  _et al._ [2017] G. Parisi, I. Procaccia, C. Rainone, and M. Singh, Shear bands as manifestation of a criticality in yielding amorphous solids, Proc. Natl. Acad. Sci. U.S.A. 114, 5577 (2017).
  * Ozawa  _et al._ [2018] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Random critical point separates brittle and ductile yielding transitions in amorphous materials, Proc. Natl. Acad. Sci. U.S.A. 115, 6656 (2018).
  * Cates  _et al._ [1998] M. Cates, J. Wittmer, J.-P. Bouchaud, and P. Claudin, Jamming, force chains, and fragile matter, Phys. Rev. Lett. 81, 1841 (1998).
  * Murphy  _et al._ [2019] K. A. Murphy, K. A. Dahmen, and H. M. Jaeger, Transforming mesoscale granular plasticity through particle shape, Phys. Rev. X 9, 011014 (2019).
  * Royer and Chaikin [2015] J. R. Royer and P. M. Chaikin, Precisely cyclic sand: Self-organization of periodically sheared frictional grains, Proc. Natl. Acad. Sci. U.S.A. 112, 49 (2015).
  * Kawasaki and Berthier [2016] T. Kawasaki and L. Berthier, Macroscopic yielding in jammed solids is accompanied by a nonequilibrium first-order transition in particle trajectories, Phys. Rev. E 94, 022615 (2016).
  * Dagois-Bohy  _et al._ [2017] S. Dagois-Bohy, E. Somfai, B. P. Tighe, and M. van Hecke, Softening and yielding of soft glassy materials, Soft Matter 13, 9036 (2017).
  * Nagasawa  _et al._ [2019] K. Nagasawa, K. Miyazaki, and T. Kawasaki, Classification of the reversible–irreversible transitions in particle trajectories across the jamming transition point, Soft Matter 15, 7557 (2019).
  * Das  _et al._ [2020] P. Das, H. Vinutha, and S. Sastry, Unified phase diagram of reversible–irreversible, jamming, and yielding transitions in cyclically sheared soft-sphere packings, Proc. Natl. Acad. Sci. U.S.A. 117, 10203 (2020).
  * Priezjev [2013] N. V. Priezjev, Heterogeneous relaxation dynamics in amorphous materials under cyclic loading, Phys. Rev. E 87, 052302 (2013).
  * Fiocco  _et al._ [2013] D. Fiocco, G. Foffi, and S. Sastry, Oscillatory athermal quasistatic deformation of a model glass, Phys. Rev. E 88, 020301 (2013).
  * Regev  _et al._ [2015] I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Reversibility and criticality in amorphous solids, Nat. Commun. 6, 1 (2015).
  * Leishangthem  _et al._ [2017] P. Leishangthem, A. D. Parmar, and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nat. Commun. 8, 1 (2017).
  * Jin  _et al._ [2018] Y. Jin, P. Urbani, F. Zamponi, and H. Yoshino, A stability-reversibility map unifies elasticity, plasticity, yielding, and jamming in hard sphere glasses, Sci. Adv. 4, eaat6387 (2018).
  * Yeh  _et al._ [2020] W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, Glass stability changes the nature of yielding under oscillatory shear, Phys. Rev. Lett. 124, 225502 (2020).
  * Keim and Arratia [2014] N. C. Keim and P. E. Arratia, Mechanical and microscopic properties of the reversible plastic regime in a 2d jammed material, Phys. Rev. Lett. 112, 028302 (2014).
  * Knowlton  _et al._ [2014] E. D. Knowlton, D. J. Pine, and L. Cipelletti, A microscopic view of the yielding transition in concentrated emulsions, Soft Matter 10, 6931 (2014).
  * Nagamanasa  _et al._ [2014] K. H. Nagamanasa, S. Gokhale, A. Sood, and R. Ganapathy, Experimental signatures of a nonequilibrium phase transition governing the yielding of a soft glass, Phys. Rev. E 89, 062308 (2014).
  * Kou  _et al._ [2017] B. Kou _et al._ , Granular materials flow like complex fluids, Nature 551, 360 (2017).
  * Kou  _et al._ [2018] B. Kou _et al._ , Translational and rotational dynamical heterogeneities in granular systems, Phys. Rev. Lett. 121, 018002 (2018).
  * Pouliquen  _et al._ [2003] O. Pouliquen, M. Belzons, and M. Nicolas, Fluctuating particle motion during shear induced granular compaction, Phys. Rev. Lett. 91, 014301 (2003).
  * Slotterback  _et al._ [2012] S. Slotterback _et al._ , Onset of irreversibility in cyclic shear of granular packings, Phys. Rev. E 85, 021309 (2012).
  * Denisov  _et al._ [2016] D. Denisov, K. Lörincz, J. Uhl, K. A. Dahmen, and P. Schall, Universality of slip avalanches in flowing granular matter, Nat. Commun. 7, 1 (2016).
  * Forterre and Pouliquen [2008] Y. Forterre and O. Pouliquen, Flows of dense granular media, Annu. Rev. Fluid Mech. 40, 1 (2008).
  * Li  _et al._ [2021] Z. Li _et al._ , Microscopic structure and dynamics study of granular segregation mechanism by cyclic shear, Sci. Adv. 7, eabe8737 (2021).
  * Yuan  _et al._ [2021] Y. Yuan _et al._ , Experimental test of the Edwards volume ensemble for tapped granular packings, Phys. Rev. Lett. 127, 018002 (2021).
  * Xing  _et al._ [2021] Y. Xing _et al._ , X-ray tomography investigation of cyclically sheared granular materials, Phys. Rev. Lett. 126, 048002 (2021).
  * MiDi [2004] G. D. R. MiDi, On dense granular flows, Eur. Phys. J. E 14, 341 (2004).
  * Marty and Dauchot [2005] G. Marty and O. Dauchot, Subdiffusion and cage effect in a sheared granular material, Phys. Rev. Lett. 94, 015701 (2005).
  * Metzler  _et al._ [2014] R. Metzler, J.-H. Jeon, A. G. Cherstvy, and E. Barkai, Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking, Phys. Chem. Chem. Phys. 16, 24128 (2014).
  * Burov [2022] S. Burov, private communications, (2022).

Methods

Quasistatic shear The cyclic shear is driven by a step motor attached to the
bottom plate of the shear box. The shear rate is $\dot{\gamma}\approx 0.13$/s,
giving a dimensionless inertia number $I=\dot{\gamma}d\sqrt{\rho/P}\approx
6\times 10^{-4}$, where we estimate the pressure to be given by $P\approx\rho
gL$ and $L\approx 24d$. This value corresponds to a quasistatic shear
condition [37].

Overlap function: For a system with $N$ particles the overlap function, which
quantifies the similarity of two configurations (here separated by a strain
$\Delta\gamma$), is defined as

$Q(\Delta\gamma)=\frac{1}{N}\sum_{i=1}^{N}\Theta(c-|\delta
x_{i}(\Delta\gamma)|),$ (3)

where $\delta x_{i}$ is the particle displacement, $\Theta(\cdot)$ is the
Heaviside step function, and $c$ is a preset threshold. By definition, $0\leq
Q\leq 1$, and $Q$ decreases as the system moves away from its initial
configuration. In practice, we divide the cubic probe space into $2\times
2\times 2$ non-overlapping subsystems, each having $N\approx 500$, to increase
the number of measurements of $Q$. Then, for given $\Gamma$ and
$\Delta\gamma$, $Q$ is sampled from different subsystems, starting
configurations, $x$-and $y$-directions, as well as $3\sim 5$ independent
realizations. To calculate $Q$ we choose $c=1.15\sqrt{\langle\delta
x^{2}(\Delta\gamma_{c})\rangle}$, which makes that $\langle
Q(\Delta\gamma_{c})\rangle=0.5$. This threshold must be chosen to depend on
$\Gamma$, since the MSD in the subdiffusive regime changes strongly with
$\Gamma$, see Fig. 1(a).

Data availability. The data that support the findings of this study are
available from the corresponding authors on reasonable request.

Acknowledgments We are grateful to S. Burov for discussions. This work is
supported by the National Natural Science Foundation of China (No. 11974240),
and the Science and Technology Commission of Shanghai Municipality (No.
22YF1419900). Y. Y. acknowledges support from the fellowship of China
Postdoctoral Science Foundation (No. 2021M702151). W. K. is a senior member of
the Institut Universitaire de France.

Author contributions Y.Y., Y.W. and W.K. designed the research. Y.Y., Z.Z.,
Y.X., H.Y. and S.Z. performed the experiment. Y.Y., W.K., and Y.W. analysed
the data and wrote the paper.

Competing interests The authors declare no competing interests.

Correspondence and requests for materials should be addressed to Yujie Wang
and Walter Kob.

Supplemental Material

Extended Data Fig. 1: Transient compaction process. (a) Packing fraction
$\phi$ as a function of accumulated strain $\gamma_{\rm acc}$ (defined
similarly to $\Delta\gamma$ but now for a dynamics that is not in steady
state) for different $\Gamma$, evolving from initially deposited packings,
i.e., random loose packings (RLP, marked by black dot) with $\phi_{\rm
RLP}\approx 0.608$, towards the steady state. This process can be
characterized by a stretched exponential law $\phi=\phi_{\rm s}+(\phi_{\rm
RLP}-\phi_{\rm s})\exp[-(\gamma_{\rm acc}/\gamma_{\rm s})^{\beta}]$, where
$\beta=0.35$, $\phi_{\rm s}$ is the steady-state packing fraction, and
$\gamma_{\rm s}$ is the compaction strain scale. (b) and (c) show,
respectively, $\phi_{s}$ and $\gamma_{s}$ versus $\Gamma$ and one observes a
crossover at $\Gamma\approx 0.1$. Error bars represent the standard deviations
from different realizations.
Extended Data Fig. 2:  Convection strength. (a) Up/down asymmetry in the
particle motion, characterized by $C(\Delta\gamma)=\langle
z\rangle/\sqrt{\langle\delta x^{2}(\Delta\gamma)\rangle}$, as a function of
$\Delta\gamma$. $C$ grows linearly with $\Delta\gamma$ with slopes that depend
on $\Gamma$, indicating an upward motion, i.e., convection. (b) The associated
slope $C/\Delta\gamma$ as a function of $\Gamma$ is basically zero for
$\Gamma\lesssim 0.1$ and increases significantly for larger $\Gamma$. This is
further evidence that the dynamics changes at $\Gamma\approx 0.1$ allowing the
convection to set in for large $\Gamma$. Error bars represent the standard
deviations from different realizations. Extended Data Fig. 3: Self-
intermediate scattering function $F_{s}(q,\Delta\gamma)$. (a)
$F_{s}(q,\Delta\gamma)$ for $q=3.5$, i.e., the first peak in static structure
factor, and different $\Gamma$. Solid curves indicate a fit with a stretched
exponential
$F_{s}(q,\Delta\gamma)=\exp[-(\Delta\gamma/\Delta\gamma_{\alpha})^{\beta_{\alpha}}]$,
where $\Delta\gamma_{\alpha}$ is the relaxation strain scale. Inset: Stretched
exponent $\beta_{\alpha}$ as a function of $q$. (b) $\Delta\gamma_{\alpha}$ as
a function of $q$ for different $\Gamma$. The scaling
$\Delta\gamma_{\alpha}\propto q^{-2}$ for small $q$ (dashed line) indicates
the Gaussian dynamics for large length scale, in agreement with the limit
$\beta_{\alpha}\rightarrow 1$ for $q\rightarrow 0$ shown in the inset of panel
(a). The colors of the curves in (a) and (b) are the same as in Fig. 1(a) of
the main text. (c) For $q=3.5$ one finds that $\Delta\gamma_{\alpha}\propto
D^{-1}$ for different $\Gamma$ (dashed line), indicating that
$F_{s}(q,\Delta\gamma)$ conveys the same information as the MSD presented in
the main text. Error bars represent the standard deviations from different
realizations. Extended Data Fig. 4:  Scaling of the MSD to determine the
yielding strain $\Delta\gamma_{c}$. (a) Power-law exponent
$\varepsilon=d\log(\langle\delta h^{2}\rangle)/d\log(\Delta\gamma)$ as a
function of $\Delta\gamma/\Delta\gamma_{c}$. The yielding strain
$\Delta\gamma_{c}$ is located at $\varepsilon=0.825$ (vertical and horizontal
dashed lines). The evolution of $\varepsilon$ can be described by
$\varepsilon=\varepsilon_{1}+(\varepsilon_{0}-\varepsilon_{1})\exp(-\Delta\gamma/\Delta\gamma_{c})$,
where $\varepsilon_{0}$ and $\varepsilon_{1}$ are fit parameters. (b)
$\Gamma$-dependence of the parameters $\varepsilon_{0}$ and $\varepsilon_{1}$.
This graph confirms that the two exponents in the MSD are independent of
$\Gamma$. Error bars represent the standard deviations from different
realizations.  Extended Data Fig. 5:  Characterizing the clusters of top
$10\%$ fastest particles, i.e., the same threshold used in Fig. 4(c). Color of
symbols are the same as in Fig. 1(a) of the main text. (a) Cluster size
distribution $P(s)$ for different $\Gamma$ at the yielding point, i.e.,
$\Delta\gamma=\Delta\gamma_{c}$. An approximately universal scaling
$P(s)\propto s^{-2}$ (dashed line) is found. $P(s)$ from the randomly chosen
$10\%$ (crosses with line) deviates strongly from this master curve. (b) The
average gyration radius $R_{g}(s)$ as a function of $s$ also shows a universal
scaling $R_{g}(s)\propto s^{0.5}$ (dashed line).