Since $\hat{R}_n$ in setting (GD) and $\check{R}_n$ in setting (GC) are consistent estimators of $R$, it suffices to check that $M_{R,r}$ is a continuous function of $R$.
To do so, we need to inspect the formulas for $M_1$ and $M_2$ in Theorems <ref> and <ref>.
The crucial point is that the eigenvalues and eigenvectors of the upper and lower diagonal blocks $R_1$ (dimension $p \times p$) and $R_2$ (dimension $q \times q$) depend continuously on $R$, since by assumption these two blocks have $p$ and $q$ distinct eigenvalues, respectively.

§ ADDITIONAL NUMERICAL EXPERIMENTS

§.§ Goodness-of-fit for finite samples for $\dep_2$

The results of Section <ref> were only presented for $\dep_1$ in view of space considerations. Here we conduct the exact same simulations using $\dep_2$. The results are presented in Figure <ref>.

[] P-P plots for 3000 repetitions of the centred and (empirically) rescaled estimator of $\dep_2$ in the three settings from Section <ref> for increasing  sample sizes (from left to right). The results are presented for an empirical correlation matrix in the case of Gaussian data.

§.§ Goodness-of-fit for rank-based estimation of the correlation matrix

We now repeat the simulations in the same settings as those of Section <ref> for the Gaussian copula case, that is when the estimated correlation matrix is $\check{R}_{n}$. The results for $\dep_1$ and $\dep_2$ are shown in Figures <ref> and <ref>, respectively.

[] P-P plots for 3000 repetitions of the centred and (empirically) rescaled estimator of $\dep_2$ in the three settings from Section <ref> for increasing  sample sizes (from left to right). The results are presented for a Gaussian copula relying on $\check{R}_{n}$.

[] P-P plots for 3000 repetitions of the centred and (empirically) rescaled estimator of $\dep_2$ for the three settings from Section <ref> for increasing  sample sizes (from left to right). The results are presented for a Gaussian copula relying on $\check{R}_{n}$.

§.§ Shrinkage evaluation for EEG data

For the EEG case study in Section <ref>, we conducted a preliminary assessment to evaluate whether the shrinkage methods in Section <ref> produce confidence intervals performing as they should.
The sample size and parameter values were taken to match those of the data.
The empirical coverage of the confidence intervals was estimated based on 2000 replications.
The results are presented in Figure <ref>. In plots (a) and (c), the advantage of shrinking the eigenvalues is clearly visible for coefficient $\dep_1$.

Empirical coverage of 95% confidence intervals estimated from 2000 replications in the Gaussian copula setting with sample size and parameters derived from the case study in Section <ref>. MLE refers to no shrinkage while DS1 and DS2 refer to the two shrinkage methods in Section <ref>. (a) Alcoholic group with $\dep_1$. (b) Alcoholic group with $\dep_2$. (c) Control group with $\dep_1$. (d) Control group with $\dep_2$.


# A Compositional Resilience Index for Computationally Efficient Safety
Analysis of Interconnected Systems

Luyao Niu1, Abdullah Al Maruf1, Andrew Clark2, J. Sukarno Mertoguno3, and
Radha Poovendran1 1Luyao Niu, Abdullah Al Maruf, and Radha Poovendran are with
the Network Security Lab, Department of Electrical and Computer Engineering,
University of Washington, Seattle, WA 98195-2500
<EMAIL_ADDRESS>Clark is with the Electrical and Systems
Engineering Department, McKelvey School of Engineering, Washington University
in St. Louis, St. Louis, MO 63130<EMAIL_ADDRESS>Sukarno Mertoguno
is with the School of Cybersecurity and Privacy, Georgia Institute of
Technology, Atlanta, GA 30332<EMAIL_ADDRESS>
###### Abstract

Interconnected systems such as power systems and chemical processes are often
required to satisfy safety properties in the presence of faults and attacks.
Verifying safety of these systems, however, is computationally challenging due
to nonlinear dynamics, high dimensionality, and combinatorial number of
possible faults and attacks that can be incurred by the subsystems
interconnected within the network. In this paper, we develop a compositional
resilience index to verify safety properties of interconnected systems under
faults and attacks. The resilience index is a tuple serving the following two
purposes. First, it quantifies how a safety property is impacted when a
subsystem is compromised by faults and attacks. Second, the resilience index
characterizes the needed behavior of a subsystem during normal operations to
ensure safety violations will not occur when future adverse events occur. We
develop a set of sufficient conditions on the dynamics of each subsystem to
satisfy its safety constraint, and leverage these conditions to formulate an
optimization program to compute the resilience index. When multiple subsystems
are interconnected and their resilience indices are given, we show that the
safety constraints of the interconnected system can be efficiently verified by
solving a system of linear inequalities. We demonstrate our developed
resilience index using a numerical case study on chemical reactors connected
in series.

##  I Introduction

Safety-critical interconnected systems are widely seen in real-world
applications such as power systems [1] and chemical processes [2]. Safety
violations can lead to significant economic losses and severe damage to the
system and/or human operators engaged with the system [1, 3, 4, 5, 6].
Therefore, it is of critical importance to verify safety properties for such
large-scale or even societal-scale systems.

One approach to verify safety is to use reachability analysis. Computing
reachable sets for nonlinear systems is known to be undecidable [7].
Alternatively, solutions to safety verification by ensuring forward invariance
of safety sets [8, 9, 10, 11] or approximating reachable sets [12, 13, 14]
have been developed. However, these approaches do not scale to interconnected
systems of high dimensions. Large-scale systems such as power systems
generally consist of multiple interconnected subsystems, motivating the
development of compositional approaches [15, 16, 17, 18]. These approaches
decompose the safety verification problem into a set of problems of smaller
scales formulated on the subsystems, and thus are more tractable.

The approaches in [15, 16, 17, 18] assume that the systems are operated under
benign environments, making the verified safety properties invalid for systems
under faults and attacks. For interconnected systems, an error from one faulty
or compromised subsystem could propagate and accumulate through
interconnections and impact the safety of other subsystems. A naïve approach
to safety verification for interconnected systems operated under adversarial
environments is to enumerate all possible faulty or compromised subsystems,
and perform safety analysis. However, the number of possible faults or attacks
that can be incurred by the interconnected system is combinatorial. At
present, scalable safety verification of large-scale interconnected systems
under faults and attacks has been less studied.

In this paper, we develop a compositional safety verification approach for
large-scale interconnected systems whose subsystems can be faulty or
compromised by attacks. Each subsystem is subject to a safety constraint. We
derive a set of conditions on the dynamics of a subsystem to guarantee its
safety. We parameterize these conditions using a tuple of real numbers, termed
resilience index. Our resilience index defines the amount of time that the
system can safely remain in a faulty state, the amount of time required to
recover from faults, and constraints on the system dynamics that must be
satisfied during faulty as well as normal operation. The resilience index
allows us to convert the problem of safety verification of large-scale
interconnected systems to a set of algebraic computations, and thus makes
safety verification feasible for large-scale systems. To summarize, this paper
makes the following contributions.

  * •

We formulate a resilient index for a subsystem that experiences faults or
attacks. We prove safety guarantees for a subsystem based on the resilience
index. We develop a sum-of-squares optimization to compute the resilience
index for a subsystem.

  * •

We derive a system of linear inequalities to quantify how the resilience index
of a subsystem changes due to interconnections. Using the derived linear
inequalities, we develop the conditions on the interconnections so that all
subsystems are safe under faults and attacks.

  * •

We demonstrate the proposed resilience indices and their usage for safety
analysis by using a numerical case study on chemical process.

The rest of this paper is organized as follows. Section II presents related
work. Section III describes the problem formulation. Section IV develops the
compositional resilience index for each subsystem. In Section V, we derive the
set of inequalities to compute resilience indices after interconnection.
Section VI demonstrates the proposed approach using a numerical case study. We
conclude the paper in Section VII.

##  II Related Work

Safety verification [19, 11, 8] and safety-critical control [20, 21, 4, 22]
haven been investigated for systems operated in benign environments. To
mitigate faults and attacks against safety-critical systems, various
techniques have been developed. Attack detection and secure state estimation
under attacks have been studied in [23, 24]. Fault-tolerant and resilient
control schemes [25, 26, 27, 9] have been proposed to withstand the attacks
and guarantee system safety. For interconnected systems consisting of multiple
subsystems, the systems are of high dimensionality and the attack surface
grows as interconnected systems involving more subsystems, making safety
verification computationally expensive.

Compositional approaches have been adopted for safety verification of
interconnected systems deployed in the absence of faults or attacks [15, 16,
17, 18]. These approaches have utilized techniques including barrier
certificates [15, 18], small-gain theorem [16], and dissipativity property
[17]. When the system is operated under faulty or adversarial environments,
these approaches become less effective.

The authors of [28, 29] re-configured the control laws of each subsystem and
interconnection topology to guarantee safety. Such approach is computationally
expensive when re-configuring the network topology and control laws.
Furthermore, in applications such as power systems, re-designing
interconnection topology is less desired or even impractical. In this paper,
we develop a compositional resilience index and prove that the safety of each
subsystem interconnected within a network can be analyzed by solving a system
of linear inequalities derived using resilience index. Our developed approach
does not require re-configuring the network topology or control laws, and
hence is more computationally efficient. In [30], the authors decomposed the
dynamics of each subsystem into intrinsic and coupled terms, where the former
term is independent of the other subsystems and the latter one depends on
interconnections. A resilience index was defined for each term of a
compromised subsystem by bounding how fast a subsystem approaches the boundary
of safety set. Such resilience indices were computed by sum-of-squares
optimization, and allowed the synthesis of safe control law under fixed
interconnections. When the interconnections change, the resilience indices and
safe control law in [30] could not always guarantee safety property. In this
paper, we propose a resilience index and derive a system of linear
inequalities that can applied to verify safety when interconnections change.
When the inequalities are feasible, the interconnected system can satisfy the
safety constraints.

##  III Problem Formulation

We first define the notations that will be used throughout this paper. Let
$x\in\mathbb{R}^{n}$. We denote the $k$-th entry of $x$ as $[x]_{k}$. A
continuous function $\alpha:[-b,a)\rightarrow(-\infty,\infty)$ belongs to
extended class $\mathcal{K}$ if it is strictly increasing and $\alpha(0)=0$
for some $a,b>0$. Linear functions $\alpha(x)=zx$ defined over $[-b,a)$ are
extended class $\mathcal{K}$ functions when $z>0$.

We consider a collection of subsystems
$\\{\mathcal{S}_{i}\\}_{i\in\mathcal{N}}$, where
$\mathcal{N}=\\{1,\ldots,N\\}$. Each subsystem $\mathcal{S}_{i}$ individually
follows dynamics given as

$\mathcal{S}_{i}:\quad\dot{x}_{i}=f_{i}(x_{i})+g_{i}(x_{i})u_{i},$ (1)

where $x_{i}\in\mathcal{X}_{i}\subseteq\mathbb{R}^{n_{i}}$ is the state of
subsystem $i$, and $u_{i}\in\mathcal{U}_{i}\subseteq\mathbb{R}^{p_{i}}$ is the
external control input $u_{i}$ applied to subsystem $i$. We consider that the
external control input $u_{i}$ will be chosen following a feedback control
policy $\mu_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathcal{U}_{i}$. Given the
control policy $\mu_{i}$ and an initial state $x_{i}(0)$ at time $t=0$, we
denote the trajectory of subsystem $\mathcal{S}_{i}$ as
$x_{i}(t;x_{i}(0),\mu_{i})$. Functions
$f_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}^{n_{i}}$ and
$g_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}^{n_{i}\times p_{i}}$ are
locally Lipschitz continuous.

We consider that each subsystem $\mathcal{S}_{i}$ is required to satisfy a
safety constraint defined over a set
$\mathcal{C}_{i}=\\{x_{i}:h_{i}(x_{i})\geq 0\\}$, i.e.,
$x_{i}(t;x_{i}(0),\mu_{i})\in\mathcal{C}_{i}$ is required to hold for all
$t\geq 0$. Function $h_{i}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}$ is
continuously differentiable. We assume that the subsystem is initially safe,
i.e., $x(0)\in\mathcal{C}_{i}$. For each subsystem $\mathcal{S}_{i}$, we
assume that we are given a control law $\mu_{i}$ such that
$x_{i}(t;x_{i}(0),\mu_{i})\in\mathcal{C}_{i}$ holds for all $t\geq 0$. Such a
safe control law $\mu_{i}$ can be synthesized by using approaches such as
control barrier functions [20].

We assume that each subsystem can be faulty or compromised by an adversary.
When the subsystem is faulty or under attack, the safe control law $\mu_{i}$
becomes offline, and the control input received by subsystem $\mathcal{S}_{i}$
can be arbitrarily altered to some $\tilde{u}_{i}\in\mathcal{U}_{i}$ that
deviates from $\mu_{i}$. To mitigate the persistence of faults and attacks, we
consider that the subsystem recovers control law $\mu_{i}$ after the
occurrence of faults and attacks leveraging fault/ attack detection and
isolation techniques [31].

To capture the fact that faults and attacks cause the control input to deviate
from $\mu_{i}$ to arbitrary $\tilde{u}_{i}$, we represent each subsystem
$\mathcal{S}_{i}$ under faults and attacks as a hybrid system
$\mathcal{H}_{i}=(\mathcal{X}_{i},\mathcal{U}_{i},\mathcal{L},\mathcal{Y}_{i},\mathcal{Y}_{i}^{init},Inv_{i},\mathcal{F}_{i},\allowbreak\Sigma_{i},\mathcal{G}_{i})$,
where

  * •

$\mathcal{X}_{i}\subseteq\mathbb{R}^{n_{i}}$ is the continuous state space of
subsystem $\mathcal{S}_{i}$, and $\mathcal{U}_{i}\subseteq\mathbb{R}^{p_{i}}$
is the set of admissible control inputs.

  * •

$\mathcal{L}=\\{\texttt{offline},\texttt{online}\\}$ is a set of discrete
locations capturing whether control law $\mu_{i}$ is available (online) or not
(offline).

  * •

$\mathcal{Y}_{i}=\mathcal{X}_{i}\times\mathcal{L}$ is the state space of
hybrid system $H$, and $\mathcal{Y}_{i}^{init}\subseteq\mathcal{Y}_{i}$ is the
set of initial states.

  * •

$Inv_{i}:\mathcal{L}\rightarrow 2^{\mathcal{X}_{i}}$ is the invariant that
maps from the set of locations to the power set of $\mathcal{X}_{i}$. That is,
$Inv_{i}(l)\subseteq\mathcal{X}_{i}$ specifies the set of possible continuous
states when the system is at location $l\in\mathcal{L}$.

  * •

$\mathcal{F}_{i}$ is the set of vector fields. For each
$F_{i}\in\mathcal{F}_{i}$ in the form of Eqn. (1), the continuous system state
evolves as $\dot{x}_{i}=F_{i}(x_{i},u_{i},l)$, where $F_{i}$ is jointly
determined by the system dynamics and the availability of control law
$\mu_{i}$, and $\dot{x}_{i}$ represents the time derivative of $x_{i}$.

  * •

$\Sigma_{i}\subseteq\mathcal{Y}_{i}\times\mathcal{Y}_{i}$ is the set of
transitions between the states of the hybrid system. A transition
$\sigma_{i}=((x_{i},l),(x_{i}^{\prime},l^{\prime}))$ models the state
transition from $(x_{i},l)$ to $(x_{i}^{\prime},l^{\prime})$.

In application such as power systems and vehicle platoons, multiple subsystems
are interconnected as a network. The interconnections introduce couplings
among subsystems, leading to the following dynamics for each subsystem

$\dot{x}_{i}=f_{i}(x_{i})+g_{i}(x_{i})\mu_{i}(x_{i})\\\ +\sum_{j\neq
i}W_{ji}(x_{j},x_{i})-\sum_{j\neq i}W_{ij}(x_{i},x_{j}),$ (2)

where $W_{ij}(x_{i},x_{j})$ captures the interconnection between subsystems
$\mathcal{S}_{i}$ and $\mathcal{S}_{j}$. Note that interconnections are not
necessarily symmetric, i.e., $W_{ij}(x_{i},x_{j})$ and $W_{ji}(x_{j},x_{i})$
may not be identical.

We denote the state and joint control input of the interconnected system as
$x=[x_{1}^{\top},\ldots,x_{N}^{\top}]^{\top}$ and
$u=[u_{1}^{\top},\ldots,u_{N}^{\top}]^{\top}$, respectively. The
interconnected system is therefore of high dimension and nonlinear.
Furthermore, when each subsystem can possibly be compromised or faulty, the
number of faults or attacks incurred by the interconnected system is
combinatorial, making safety verification computationally intractable. In this
paper, we investigate the following problem.

######  Problem 1.

Suppose that we are given a collection of subsystems
$\\{\mathcal{S}_{i}\\}_{i=1}^{N}$ and their safe control laws $\mu_{i}$ with
respect to their individual safety set $\mathcal{C}_{i}$, where
$i\in\mathcal{N}$. The subsystems, which are potentially subject to faults and
attacks, are interconnected within a network and each of them follows dynamics
given in Eqn. (2). The goal is to verify whether the interconnected system
satisfies the set of safety constraints defined over $\mathcal{C}_{i}$ for all
$i=1,\ldots,N$.

##  IV Proposed Resilience Index

In this section, we propose a compositional resilience index for each
subsystem to verify safety.

###  IV-A Definition of Resilience Index

We note that the discrete location set $\mathcal{L}$ is uniform to all
subsystems, allowing us to develop a unified index to measure the resilience
of any subsystem under faults and attacks. Our insight is as follows. At
location offline, the safe control law $\mu_{i}$ is unavailable. To avoid
violating the safety constraint, we require the subsystem to stay within a set
$\mathcal{D}_{i}\subseteq\mathcal{C}_{i}$ so that $x_{i}\in\mathcal{C}_{i}$
for all time when the hybrid system is at location offline. When the hybrid
system transitions from offline to online, the control law $\mu_{i}$ becomes
available. Thereafter, the subsystem starts to _recover_ from faults and
attacks. To recover the control law and mitigate potential faults and attacks
in the future, we let the control law $\mu_{i}$ to remain available for at
least $\phi_{i}\geq 0$ amount of time. Furthermore, the system is required to
reach set $\mathcal{D}_{i}$ within $\phi_{i}$ so that safety constraint will
not be violated when attacks or faults occur in the future. Such insight
allows us to define the following resilience index to capture each subsystem’s
resilience under faults and attacks.

######  Definition 1 (Resilience Index of a Subsystem).

Consider a subsystem $\mathcal{S}_{i}$ that uses a feedback control law
$\mu_{i}$ and is under a safety constraint defined on set
$\mathcal{C}_{i}=\\{x_{i}:h_{i}(x_{i})\geq 0\\}$. Let subsystem
$\mathcal{S}_{i}$ be formulated as hybrid system $\mathcal{H}_{i}$ and
$d_{i},\eta_{i}\geq 0$, $\tau_{i},\phi_{i}>0$. We say subsystem
$\mathcal{S}_{i}$ is $(d_{i},\tau_{i},\phi_{i},\eta_{i})$-resilient if the
following conditions hold

  * •

Set $\mathcal{D}_{i}$ defined as
$\mathcal{D}_{i}=\\{x_{i}:h_{i}(x_{i})-d_{i}\geq 0\\}$ is forward invariant if
control law $\mu_{i}$ is used and the hybrid system $\mathcal{H}_{i}$ in at
location _online_.

  * •

After reaching location _offline_ , hybrid system $\mathcal{H}_{i}$ remains at
location _offline_ for at most $\tau_{i}$ amount of time before transition
from _offline_ to _online_ occurs for any $x_{i}$.

  * •

Following the transition from _offline_ to _online_ , hybrid system
$\mathcal{H}_{i}$ remains at location _online_ for at least $\phi_{i}$ amount
of time.

  * •

When $x_{i}$ is at the boundary of $\mathcal{D}_{i}$ and the hybrid system
$\mathcal{H}_{i}$ in at location _online_ , the time derivative of $x_{i}$ is
lower bounded by $\eta_{i}$. Furthermore, given any state
$x_{i}^{\prime}\in\mathcal{C}_{i}\setminus\mathcal{D}_{i}$, the continuous
state $x_{i}$ reaches $\mathcal{D}_{i}$ within $\phi_{i}$ amount of time when
$\mathcal{H}_{i}$ is at location _online_ and control law $\mu_{i}$ is used.

The quadruple $(d_{i},\tau_{i},\phi_{i},\eta_{i})$ is the resilience index of
$\mathcal{S}_{i}$.

In what follows, we derive a set of conditions to compute the resilience index
of a subsystem.

######  Proposition 1.

Consider a subsystem in Eqn. (1) under attack and a safety set
$\mathcal{C}_{i}$. Let $h_{i}^{d}(x_{i})=h_{i}(x_{i})-d_{i}$ and
$\mathcal{D}_{i}=\\{x_{i}:h_{i}^{d}(x_{i})\geq 0\\}\subseteq\mathcal{C}_{i}$.
Suppose $x_{i}(0)\in\mathcal{D}_{i}$. If there exist constants
$d_{i},\eta_{i}\geq 0$, $\tau_{i},\phi_{i}>0$, and an extended class
$\mathcal{K}$ function $\alpha_{i}(\cdot)$ such that

$\displaystyle\frac{\partial h_{i}^{d}}{\partial
x_{i}}(x_{i})[f_{i}(x_{i})+g_{i}(x_{i})u_{i}]\geq-\frac{d_{i}}{\tau_{i}},~{}\forall(x_{i},u_{i})\in\mathcal{C}_{i}\times\mathcal{U}_{i}$
(3a) $\displaystyle\frac{\partial h_{i}^{d}}{\partial
x_{i}}(x_{i})[f_{i}(x_{i})+g_{i}(x_{i})\mu_{i}(x_{i})]\geq\frac{d_{i}}{\phi_{i}},~{}\forall
x_{i}\in\mathcal{C}_{i}\setminus\mathcal{D}_{i}$ (3b)
$\displaystyle\frac{\partial h_{i}^{d}}{\partial
x_{i}}(x_{i})[f_{i}(x_{i})+g_{i}(x_{i})\mu_{i}(x_{i})]\geq-\alpha_{i}(h_{i}^{d}(x_{i}))+\eta_{i},$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad~{}\forall
x\in\mathcal{D}_{i}$ (3c)

then subsystem $\mathcal{S}_{i}$ is safe with respect to $\mathcal{C}_{i}$.

###### Proof.

Without loss of generality, we assume that the subsystem is compromised by
attack at some arbitrary time $t_{0}\geq 0$. For any
$t\in[t_{0},t_{0}+\tau_{i}]$ and control input $u_{i}\in\mathcal{U}_{i}$, we
have
$h_{i}(x_{i}(t))=h_{i}(x_{i}(t_{0}))+\int_{s=t_{0}}^{t}\dot{h}_{i}\,\textrm{d}{s}\geq
d_{i}-\frac{d_{i}}{\tau_{i}}(t-t_{0})\geq 0$, where $\dot{h}_{i}$ represents
the time derivative of function $h_{i}$, the equality holds by the definition
of $h_{i}(x_{i}(t))$, the first inequality holds by Eqn. (3a) and the
observation that $\dot{h}_{i}^{d}=\dot{h}_{i}$, and the last inequality holds
by $t\in[t_{0},t_{0}+\tau]$. Therefore, hybrid system $\mathcal{H}_{i}$
satisfies $x_{i}(t;x_{i}(t_{0}),u_{i})\in\mathcal{C}_{i}$ for all
$u_{i}\in\mathcal{U}_{i}$ when $x_{i}(t_{0})\in\mathcal{D}_{i}$,
$\mathcal{H}_{i}$ is at location $l=\text{{offline}}$ for at most $\tau_{i}$
amount of time, and Eqn. (3a) holds.

We next consider some arbitrary time
$t\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ when the hybrid system is at
location online. We show that
$x_{i}(t;x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{C}_{i}$ for all
$t\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ when the hybrid system remains
at location online and control law $\mu_{i}$ is used. We further show that
there exists $t^{\prime}\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ such that
$x_{i}(t^{\prime\prime};x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{D}_{i}$ for
all $t^{\prime\prime}\in[t^{\prime},t_{0}+\tau_{i}+\phi_{i}]$ when the hybrid
system remains at location online and control law $\mu_{i}$ is used.

We prove $x_{i}(t;x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{C}_{i}$ for all
$t\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ when the hybrid system remains
at location online and $\mu_{i}$ is used by contradiction. Suppose that the
subsystem leaves the safety set. Since the trajectory of $\mathcal{S}_{i}$ is
continuous, then there exists time
$t\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ such that $h(x(t))=0$ and
$\dot{h}(x(t))<0$. If such $x(t)\in\mathcal{C}_{i}\setminus\mathcal{D}_{i}$,
we then have contradiction to Eqn. (3b) since $\frac{d_{i}}{\phi_{i}}>0$ and
thus $\dot{h}(x(t))>0$. If such $x(t)\in\mathcal{D}_{i}$, we also have
contradiction since Eqn. (3c) implies that $\dot{h}(x(t))>\eta_{i}\geq 0$.
Therefore, we can claim that
$x_{i}(t;x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{C}_{i}$ for all
$t\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ when the hybrid system remains
at location online and $\mu_{i}$ is used.

We now prove that there exists
$t^{\prime}\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ such that
$x_{i}(t^{\prime\prime};x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{D}_{i}$ for
all $t^{\prime\prime}\in[t^{\prime},t_{0}+\tau_{i}+\phi_{i}]$ when the hybrid
system remains at location online and control law $\mu_{i}$ is used. Suppose
no such $t^{\prime}\in[t_{0}+\tau_{i},t_{0}+\tau_{i}+\phi_{i}]$ exists. We
then have that
$x_{i}(t_{0}+\tau_{i}+\phi_{i})\in\mathcal{C}_{i}\setminus\mathcal{D}_{i}$ and
hence $h_{i}(x_{i}(t_{0}+\tau_{i}+\phi_{i}))<d_{i}$. By Eqn. (3b), we have
$h_{i}(x_{i}(t_{0}+\tau_{i}+\phi_{i}))=h_{i}(x_{i}(t_{0}+\tau_{i}))+\int_{s=t_{0}+\tau_{i}}^{t}\dot{h}_{i}\,\textrm{d}{s}\geq\frac{d_{i}}{\phi_{i}}(t-t_{0}-\tau_{i})\geq
d_{i}$, leading to contradiction. Therefore, such $t^{\prime}$ must exist.
Finally, by Eqn. (3c) and [20], we have that if
$x_{i}(t^{\prime};x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{D}_{i}$, then set
$\mathcal{D}_{i}$ if forward invariant when the hybrid system remains at
location online and control law $\mu_{i}$ is used, indicating that
$x_{i}(t^{\prime\prime};x_{i}(t_{0}+\tau_{i}),\mu_{i})\in\mathcal{D}_{i}$ for
all $t^{\prime\prime}\in[t^{\prime},t_{0}+\tau_{i}+\phi_{i}]$. ∎

###  IV-B Computation of Resilience Index

In the following, we formulate a sum-of-squares (SOS) optimization program to
compute the resilience index for any subsystem $\mathcal{S}_{i}$. Under
certain assumptions on the dynamics (1), function $h_{i}(x)$, and control
input set $\mathcal{U}_{i}$, we formulate the SOS program by converting the
conditions in Eqn. (3) into SOS constraints. We make the following assumption.

######  Assumption 1.

For any subsystem $\mathcal{S}_{i}$, we assume that functions $f_{i}(x_{i})$,
$g_{i}(x_{i})$, and $h_{i}(x_{i})$ are polynomial in $x_{i}$. In addition, we
assume that $\mathcal{U}_{i}=\prod_{k=1}^{p_{i}}[u_{k,min},u_{k,max}]$ with
$u_{k,min}<u_{k,max}$.

In the following, we present the set of SOS constraints. We show that any
$d_{i},\eta_{i}\geq 0$ and $\tau_{i},\phi_{i}>0$ satisfying the SOS
constraints constitute the resilience index of $\mathcal{S}_{i}$.

######  Proposition 2.

Assume that control law $\mu_{i}$ is polynomial in $x_{i}$. Suppose there
exist $d_{i},\eta_{i}\geq 0$, and $\tau_{i},\phi_{i}>0$ such that the
following expressions are SOS:

$\displaystyle\frac{\partial h_{i}^{d}}{\partial
x_{i}}(x_{i})[f_{i}(x_{i})+g_{i}(x_{i})u_{i}]+d_{i}\theta_{i}-q(x_{i},u_{i})h_{i}^{d}(x_{i})$
$\displaystyle\quad\quad\quad-\sum_{k=1}^{p_{i}}(w_{k}(x_{i},u_{i})([u]_{k}-[u]_{k,min})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+v_{k}(x_{i},u_{i})([u]_{k,max}-[u]_{k})),$
(4a) $\displaystyle\frac{\partial h_{i}^{d}}{\partial
x_{i}}(x_{i})[f_{i}(x_{i})+g_{i}(x_{i})\mu_{i}(x_{i})]-d_{i}\beta_{i}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad-l(x_{i})h_{i}^{d}(x)+m(x_{i})h_{i}(x_{i}),$
(4b) $\displaystyle\frac{\partial h_{i}^{d}}{\partial
x_{i}}(x_{i})[f_{i}(x_{i})+g_{i}(x_{i})\mu_{i}(x_{i})]$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad+\alpha(h_{i}^{d}(x))-r(x_{i})h_{i}(x_{i})-\eta_{i},$
(4c)

where $l(x_{i}),m(x_{i}),q(x_{i},u_{i}),r(x_{i})$ are SOS, and
$w_{k}(x_{i},u_{i})$ as well as $v_{k}(x_{i},u_{i})$ are SOS for each
$k=1,\ldots,p_{i}$. Then $d_{i}$, $\tau_{i}=\frac{1}{\theta_{i}}$,
$\phi_{i}=\frac{1}{\beta_{i}}$, and $\eta_{i}$ satisfy Eqn. (3).

###### Proof.

We prove that Eqn. (4a) implies Eqn. (3a). The other SOS constraints can be
proved in a similar manner. Consider $x_{i}\in\mathcal{C}_{i}$ and
$[u]_{k,min}\leq[u]_{k}\leq[u]_{k,max}$ for all $k=1,\ldots,p_{i}$. We thus
have that $h_{i}^{d}(x_{i})\geq 0$ since
$\mathcal{D}_{i}\subseteq\mathcal{C}_{i}$. In addition, we have that
$[u]_{k}-[u]_{k,min}\geq 0$, and $[u]_{k,max}-[u]_{k}\geq 0$. When expression
(4a) is SOS, $q(x_{i},u_{i})$, $w_{k}(x_{i},u_{i})$, and $v_{k}(x_{i},u_{i})$
are SOS for all $k=1,\ldots,p_{i}$, we have that the expression in Eqn. (4a)
is non-negative. Therefore, if expression (4a) is SOS and $\tau_{i}$ is chosen
as $\tau_{i}=1/\theta_{i}$, then Eqn. (3a) holds. ∎

We observe that the SOS constraints derived in Proposition 2 are bilinear (see
the terms $d_{i}\theta_{i}$ and $d_{i}\beta_{i}$). Hence the resilience index
cannot be readily computed by implementing these SOS constraints. We overcome
this challenge by developing an alternating optimization procedure, as shown
in Algorithm 1. In Algorithm 1, parameters $\tau_{max}$ and $\phi_{min}$ are
the upper and lower bounds for $\tau_{i}$ and $\phi_{i}$, respectively. If
there exist no bound for $\tau_{i}$ and $\phi_{i}$, then parameters
$\tau_{max}$ and $\phi_{min}$ can be set as infinity and zero, respectively.

Algorithm 1 Algorithm for computing resilience index

1:Input: $\mathcal{S}_{i}$, $\tau_{max}$, $\phi_{min}$, and $\epsilon>0$

2:Output: $d_{i}$, $\tau_{i}$, $\phi_{i}$, $\eta_{i}$

3:Initialization: $d_{i}=0$

4:while $d_{i}\leq\sup_{x_{i}}\\{h_{i}(x_{i})\\}$ do

5: Solve for $\theta_{i}$, $\beta_{i}$, $\eta_{i}$ such that (4) is feasible
with $d_{i}$ fixed.

6: if Eqn. (4) is feasible, $\frac{1}{\theta_{i}}\leq\tau_{max}$, and
$\frac{1}{\beta_{i}}\geq\tau_{min}$ then

7: return $d_{i}$, $\tau_{i}=\frac{1}{\theta_{i}}$,
$\phi_{i}=\frac{1}{\beta_{i}}$, and $\eta_{i}$

8: else

9: $d_{i}=d_{i}+\epsilon$

10: end if

11:end while

##  V Resilience Index After Interconnection

In this section, we consider a setting where multiple subsystems, with each
being formulated by a hybrid system $\mathcal{H}_{i}$, are interconnected
within a network. Suppose that a collection of subsystems
$\\{\mathcal{S}_{i}\\}_{i=1}^{N}$ are interconnected within a network. In the
network, each subsystem $\mathcal{S}_{i}$ follows the dynamics as given by
Eqn. (2).

We note that the interconnected system can be formulated as a hybrid system
$\mathcal{H}$ as defined in Section IV-A. In this case, the continuous state
space is $\mathcal{X}=\prod_{i=1}^{N}\mathcal{X}_{i}$, and the discrete
location is $\mathcal{L}=\\{\text{{offline}},\text{{online}}\\}^{N}$. We
observe that the continuous state is of dimension $n=\sum_{i=1}^{N}n_{i}$.
Furthermore, the transitions among the discrete locations are the Cartesian
product of discrete transitions of all subsystems, which is combinatorial in
nature to capture all possible faults and attacks that can be incurred by the
subsystems. Therefore, safety verification over the hybrid system
$\mathcal{H}$ is computationally intractable for large-scale interconnected
systems. In what follows, we derive how the resilience index of each subsystem
changes due to interconnections. We further show how our proposed resilience
index can be applied to efficiently verify safety constraints of the
interconnected system.

###  V-A Computation of Resilience Index After Interconnection

In the following, we first characterize the behaviors of any subsystem
$\mathcal{S}_{j}$ when being interconnected. We define

$\delta_{j}=\inf_{x}\Big{\\{}\frac{\partial h_{j}^{d}}{\partial
x_{j}}\big{[}\sum_{i\neq j}W_{ij}(x_{i},x_{j})-\sum_{i\neq
j}W_{ji}(x_{j},x_{i})\big{]}\Big{\\}}.$ (5)

We will show that the resilience index of a subsystem $\mathcal{S}_{j}$ after
being interconnected can be bounded using one of the following two sets of
inequalities

$\displaystyle R_{1}:\begin{cases}0\leq d_{j}^{\prime}\leq
d_{j}\leq\sup_{x_{j}}\\{h_{j}(x_{j})\\}\\\
-\frac{d_{j}^{\prime}}{\tau_{j}^{\prime}}\leq-\frac{d_{j}}{\tau_{j}}+\delta_{j}\\\
\frac{\phi_{j}}{d_{j}+\phi_{j}\delta_{j}}d_{j}^{\prime}-\phi_{j}^{\prime}\leq
0\\\
\eta_{j}^{\prime}\leq\delta_{j}+\min\\{\frac{d_{j}}{\phi_{j}},\eta_{j}+\inf_{x_{j}\in\mathcal{D}_{j}}\\{\alpha_{j}(h_{j}(x_{j})-d_{j}^{\prime})\\\
\quad\quad\quad\quad\quad\quad\quad\quad-\alpha_{j}(h_{j}(x_{j})-d_{j})\\}\\}\end{cases}$
(6) $\displaystyle R_{2}:\begin{cases}0\leq d_{j}\leq
d_{j}^{\prime}\leq\sup_{x_{j}}\\{h_{j}(x_{j})\\}\\\
-\frac{d_{j}^{\prime}}{\tau_{j}^{\prime}}\leq-\frac{d_{j}}{\tau_{j}}+\delta_{j}\\\
\frac{d_{j}^{\prime}}{\phi_{j}^{\prime}}\leq\delta_{j}+\min\\{\frac{d_{j}}{\phi_{j}},\inf_{x_{j}\in\mathcal{D}_{j}}\\{-\alpha_{j}(h_{j}^{d}(x_{j}))\\}+\eta_{j}\\}\\\
\eta_{j}^{\prime}\leq\delta_{j}+\eta_{j}-\sup_{x_{j}\in\mathcal{D}_{j}}\\{\alpha_{j}(h_{j}(x_{j})-d_{j}^{\prime})\\\
\quad\quad\quad\quad\quad\quad\quad\quad-\alpha_{j}(h_{j}(x_{j})-d_{j})\\}\\}\end{cases}$
(7)

We define $\mathcal{D}_{j}^{\prime}=\\{x_{j}:h_{j}(x_{j})-d_{j}^{\prime}\geq
0\\}$. The inequalities in $R_{1}$ and $R_{2}$ specify sets
$\mathcal{D}_{j}^{\prime}$ differently. The inequalities in $R_{1}$ specify
that $\mathcal{D}_{j}^{\prime}\supseteq\mathcal{D}_{j}$, whereas $R_{2}$
defines $\mathcal{D}_{j}^{\prime\prime}\subseteq\mathcal{D}_{j}$, which
further leads to distinct behaviors when
$x_{j}\in\mathcal{C}_{j}\setminus\mathcal{D}_{j}^{\prime}$. In what follows,
we show how the behavior of each subsystem following dynamics in Eqn. (2) can
be characterized by the solutions to $R_{1}$ or $R_{2}$. This allows us to
further verify the safety constraints for the interconnected system.

######  Theorem 1.

Consider that a collection of subsystems $\\{\mathcal{S}_{j}\\}_{j=1}^{N}$ are
interconnected, and each $\mathcal{S}_{j}$ follows dynamics as given in Eqn.
(2) for all $j=1,\ldots,N$. We denote their resilience indices before being
interconnected as $(d_{k},\tau_{k},\phi_{k},\eta_{k})$, where $k=1,\ldots,N$.
Define $\mathcal{D}_{j}^{\prime}=\\{x_{j}:h_{j}(x_{j})-d_{j}^{\prime}\geq
0\\}$. If parameters $d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime}$, and
$\eta_{j}^{\prime}$ render either $R_{1}$ or $R_{2}$ to be feasible, then the
following conditions hold for $\mathcal{S}_{j}$ after being interconnected:

$\displaystyle\frac{\partial h_{j}^{d}}{\partial
x_{j}}(x_{j})[f_{j}(x_{j})+g_{j}(x_{j})u_{j}+\sum_{i\neq
j}W_{ij}(x_{i},x_{j})$ (8a) $\displaystyle\quad-\sum_{i\neq
j}W_{ji}(x_{j},x_{i})]\geq-\frac{d_{j}^{\prime}}{\tau_{j}^{\prime}},~{}\forall(x_{j},u_{j})\in\mathcal{C}_{j}\times\mathcal{U}_{j}$
$\displaystyle\frac{\partial h_{j}^{d}}{\partial
x_{j}}(x_{j})[f_{j}(x_{j})+g_{j}(x_{j})\mu_{j}(x_{j})+\sum_{i\neq
j}W_{ij}(x_{i},x_{j})$ (8b) $\displaystyle\quad-\sum_{i\neq
j}W_{ji}(x_{j},x_{i})]\geq\frac{d_{j}^{\prime}}{\phi_{j}^{\prime}},~{}\forall
x_{j}\in\mathcal{C}_{j}\setminus\mathcal{D}_{j}^{\prime}$
$\displaystyle\frac{\partial h_{j}^{d}}{\partial
x_{j}}(x_{j})[f_{j}(x_{j})+g_{j}(x_{j})\mu_{j}(x_{j})+\sum_{i\neq
j}W_{ij}(x_{i},x_{j})$ (8c) $\displaystyle\quad-\sum_{i\neq
j}W_{ji}(x_{j},x_{i})]\geq-\alpha_{j}(h_{j}(x_{j})-d_{j}^{\prime})+\eta_{j}^{\prime},~{}\forall
x\in\mathcal{D}_{j}^{\prime}$

###### Proof.

We first verify that if $d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime}$,
and $\eta_{j}^{\prime}$ satisfy $R_{1}$, then Eqn. (8) holds. We denote
$\dot{h}_{j}^{d}$ as

$\dot{h}_{j}^{d}=\frac{\partial h_{j}^{d}}{\partial
x_{j}}(x_{j})[f_{j}(x_{j})+g_{j}(x_{j})\mu_{j}(x_{j})\\\ +\sum_{i\neq
j}W_{ij}(x_{i},x_{j})-\sum_{i\neq j}W_{ji}(x_{j},x_{i})].$

Suppose that $d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime}$, and
$\eta_{j}^{\prime}$ yield $R_{1}$ to be feasible. In this case, we have that
$\mathcal{D}_{j}=\\{x_{j}:h_{j}(x_{j})-d_{j}\geq
0\\}\subseteq\mathcal{D}_{j}^{\prime}=\\{x_{j}:h_{j}(x_{j})-d_{j}^{\prime}\geq
0\\}$ due to $d_{j}^{\prime}\leq d_{j}$. When $\mathcal{H}_{j}$ is at location
offline, we have

$\displaystyle\dot{h}_{j}^{d}\geq$ $\displaystyle\frac{\partial
h_{j}^{d}}{\partial x_{j}}(x_{j})[f_{j}(x_{j})+g_{j}(x_{j})u_{j}]+\delta_{j}$
(9) $\displaystyle\geq$
$\displaystyle-\frac{d_{j}}{\tau_{j}}+\delta_{j}\geq-\frac{d_{j}^{\prime}}{\tau_{j}^{\prime}},~{}\forall(x_{j},u_{j})\in\mathcal{C}_{j}\times\mathcal{U}_{j}$
(10)

where inequality (9) holds Eqn. (5), the second inequality holds by Eqn. (3a),
and the last inequality holds by the assumption that $R_{1}$ is feasible under
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$.

Consider the case where hybrid system $\mathcal{H}_{j}$ is at location online
and control law $\mu_{j}(x_{j})$ is available. We have

$\dot{h}_{j}^{d}\geq\frac{d_{j}}{\phi_{j}}+\delta_{j}\geq\frac{d_{j}^{\prime}}{\phi_{j}^{\prime}}$
(11)

holds for all $x_{j}\in\mathcal{C}_{j}\setminus\mathcal{D}_{j}^{\prime}$,
where the first inequality holds by Eqn. (3b) and (5), and the second
inequality holds by
$(\mathcal{C}_{j}\setminus\mathcal{D}_{j}^{\prime})\subseteq(\mathcal{C}_{j}\setminus\mathcal{D}_{j})$
given the feasibility of $R_{1}$.

We finally consider the case where $x_{j}\in\mathcal{D}_{j}^{\prime}$ by
dividing our discussion into two scenarios. When
$x_{j}\in\mathcal{D}_{j}^{\prime}\setminus\mathcal{D}_{j}$, Eqn. (11) yields
$\dot{h}_{j}^{d}\geq\frac{d_{j}}{\phi_{j}}+\delta_{j}\geq\eta_{j}^{\prime}$,
where the last inequality holds by the feasibility of $R_{1}$. When
$x_{j}\in\mathcal{D}_{j}\subseteq\mathcal{D}_{j}^{\prime}$, we have that

$\displaystyle\dot{h}_{j}^{d}\geq-\alpha(h_{j}^{d}(x_{j}))+\eta_{j}+\delta_{j}\geq-\alpha(h_{j}(x_{j})-d_{j}^{\prime})+\eta_{j}^{\prime}$
(12)

holds for all $x_{j}\in\mathcal{D}_{j}$, where the last inequality holds by
the feasibility of $R_{1}$, i.e.,
$\eta_{j}^{\prime}\leq\delta_{j}+\eta_{j}+\inf_{x_{j}\in\mathcal{D}_{j}}\\{\alpha(h_{j}(x_{j})-d_{j}^{\prime})-\alpha(h_{j}(x_{j})-d_{j})\\}$.

We next verify that if $d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime}$,
and $\eta_{j}^{\prime}$ satisfy $R_{2}$, then Eqn. (8) holds. Note that in
this case, $\mathcal{D}_{j}^{\prime}=\\{x_{j}:h_{j}(x_{j})-d_{j}^{\prime}\geq
0\\}\subseteq\mathcal{D}_{j}=\\{x_{j}:h_{j}(x_{j})-d_{j}\geq 0\\}$. When
hybrid system $\mathcal{H}_{j}$ is at location offline, Eqn. (8a) can be
derived using Eqn. (9) and (10) given that
$d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime}$, and $\eta_{j}^{\prime}$
satisfy $R_{2}$. We next consider that hybrid system $\mathcal{H}_{j}$ is at
location online. We discuss two possible scenarios that can occur when
$x_{j}\in\mathcal{C}_{j}\setminus\mathcal{D}_{j}^{\prime}$. If
$x_{j}\in\mathcal{C}_{j}\setminus\mathcal{D}_{j}$, we have

$\dot{h}_{j}^{d}\geq\frac{d_{j}}{\phi_{j}}+\delta_{j}\geq\frac{d_{j}^{\prime}}{\phi_{j}^{\prime}},~{}\forall
x_{j}\in\mathcal{C}_{j}\setminus\mathcal{D}_{j}$ (13)

where the first inequality holds by Eqn. (3b) and the definition of
$\delta_{j}$, and the second inequality holds by the feasibility of $R_{2}$.
If $x_{j}\in\mathcal{D}_{j}\setminus\mathcal{D}_{j}^{\prime}$, we have that

$\dot{h}_{j}^{d}\geq-\alpha(h_{j}^{d}(x_{j}))+\delta_{j}\geq\frac{d_{j}^{\prime}}{\phi_{j}^{\prime}},~{}\forall
x_{j}\in\mathcal{D}_{j}\setminus\mathcal{D}_{j}^{\prime},$ (14)

where the first inequality holds by Eqn. (3c) and the definition of
$\delta_{j}$, and the second inequality holds by the feasibility of $R_{2}$.
Combining Eqn. (13) and (14) yields Eqn. (8b).

We finally consider that $x_{j}\in\mathcal{D}_{j}^{\prime}$. We have that

$\dot{h}_{j}^{d}\geq-\alpha(h_{j}^{d}(x_{j}))+\eta_{j}+\delta_{j}\geq-\alpha(h_{j}(x_{j})-d_{j}^{\prime})+\eta_{j}^{\prime}$

holds for all $x_{j}\in\mathcal{D}_{j}^{\prime}$, where the first inequality
holds by Eqn. (3c) and the definition of $\delta_{j}$, and the second
inequality holds by $\mathcal{D}_{j}^{\prime}\subseteq\mathcal{D}_{j}$ along
with the feasibility of $R_{2}$.

Combining the discussion above completes the proof. ∎

We observe that when function $\alpha_{j}$ is linear, computing the resilience
indices after interconnection reduces to solving a linear system. In the
following, we show that given the resilience indices of $\mathcal{S}_{j}$
before it is interconnected along with its control law $\mu_{i}$, we can
efficiently quantify how its resilience index changes due to interconnections
by solving a set of inequalities given in Eqn. (6) and (7).

######  Theorem 2.

Consider that a collection of subsystems $\\{\mathcal{S}_{j}\\}_{j=1}^{N}$ are
interconnected, and each $\mathcal{S}_{j}$ follows dynamics as given in Eqn.
(2) for all $j=1,\ldots,N$. We denote their resilience indices before being
interconnected as $(d_{k},\tau_{k},\phi_{k},\eta_{k})$, where $k=1,\ldots,N$.
If parameters $d_{j}^{\prime},\eta_{j}^{\prime}\geq 0$ and
$\tau_{j}^{\prime},\phi_{j}^{\prime}>0$ satisfy either $R_{1}$ in Eqn. (6) or
$R_{2}$ in Eqn. (7), then $\mathcal{S}_{j}$ is
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$-resilient
under dynamics (2). Furthermore, $\mathcal{S}_{j}$ is safe with respect to
$\mathcal{C}_{j}$ after being interconnected within the network.

###### Proof.

When parameters $d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime}$, and
$\eta_{j}^{\prime}$ render either $R_{1}$ or $R_{2}$ to be feasible, we have
that Eqn. (8) holds by using Theorem 1. By using Proposition 1 and Definition
1, we have that when $d_{j}^{\prime},\eta_{j}^{\prime}\geq 0$ and
$\tau_{j}^{\prime},\phi_{j}^{\prime}>0$, subsystem $\mathcal{S}_{j}$ is
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$-resilient
and safe with respect to $\mathcal{C}_{j}$. ∎

Computing $\delta_{j}$ for each subsystem $\mathcal{S}_{j}$ requires to solve
an optimization problem over joint system state $x\in\mathbb{R}^{n}$, where
$n=\sum_{i=1}^{N}n_{i}$. To alleviate the computations in high dimensional
state space, we approximate parameter $\delta_{j}$ as

$\delta_{j}\geq\sum_{i\neq
j}\inf_{(x_{i},x_{j})\in\mathcal{C}_{i}\times\mathcal{C}_{j}}\\{W_{ij}(x_{i},x_{j})-W_{ji}(x_{j},x_{i})\\}:=\tilde{\delta}_{j}.$

By replacing $\delta_{j}$ with $\tilde{\delta}_{j}$ in Eqn. (8), we can apply
similar approach to show that Theorem 1 and 2 still hold.

###  V-B Feasibility of Resilience Index for Interconnected System

Consider an interconnected system consisting of $N$ subsystems. We need to
determine whether the inequalities in $R_{1}$ or $R_{2}$ need to be solved to
apply Theorem 1 and 2 for safety verification of the interconnected system.
One approach is to combine the inequalities in $R_{1}$ and $R_{2}$ by using a
set of mixed integer constraints and big M-method [32], where the integer
variable $y\in\\{0,1\\}$ models whether $R_{1}$ or $R_{2}$ is solved. In this
subsection, we show that we can determine whether $R_{1}$ or $R_{2}$ is
feasible given the value of $\delta_{j}$, and hence avoid solving the mixed
integer program.

######  Theorem 3.

Consider a subsystem $\mathcal{S}_{j}$ whose resilience index is given as
$(d_{j},\tau_{j},\phi_{j},\eta_{j})$ before being interconnected. Assume that
$\alpha_{j}(h_{j}(x_{j}))=zh_{j}(x_{j})$ for some coefficient $z>0$. If
$\delta_{j}$ satisfies

$\delta_{j}\geq\max\\{-\frac{d_{j}}{\phi_{j}},-\eta_{j}-zd_{j}\\},$ (15)

then there exist $d_{j}^{\prime},\eta_{j}^{\prime}\geq 0$ and
$\tau_{j}^{\prime},\phi_{j}^{\prime}>0$ such that the inequalities in $R_{1}$
are satisfied. If $\delta_{j}$ satisfies

$\delta_{j}\geq\max\\{-\frac{d_{j}}{\phi_{j}},-\eta_{j}+z(\sup_{x_{j}}\\{h_{j}(x_{j})\\}-d_{j})\\},$
(16)

then there exist $d_{j}^{\prime},\eta_{j}^{\prime}\geq 0$ and
$\tau_{j}^{\prime},\phi_{j}^{\prime}>0$ such that the inequalities in $R_{2}$
are satisfied.

###### Proof.

Suppose that Eqn. (15) holds. We rewrite $R_{1}$ in the matrix form
$A_{j}v_{j}\leq b_{j}$, where
$v_{j}=[d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime}]^{\top}$,

$A=\begin{bmatrix}1&0&0&0\\\ 1&-\frac{d_{j}}{\phi_{j}}+\delta_{j}&0&0\\\
z&0&0&1\\\ -1&0&\frac{d_{j}}{\tau_{j}}-\delta_{j}&0\\\
0&1&1&1\end{bmatrix},b=\begin{bmatrix}d_{j}\\\ 0\\\
\delta_{j}+\eta_{j}+zd_{j}\\\ 0\\\
\delta_{j}+\frac{d_{j}}{\phi_{j}}\end{bmatrix}$

When Eqn. (15) holds, we have that there exists no $r_{j}\geq 0$ such that
$b_{j}^{\top}r_{j}<0$. Using Farkas’ Lemma [33], there must exist some
$v_{j}\geq 0$ such that $A_{j}v_{j}\leq b_{j}$, and thus satisfies $R_{1}$.
Similar proof technique can be applied to show that when Eqn. (16) holds,
there must exist non-negative
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$ such
that $R_{2}$ is satisfied. Noticing that
$\phi_{j}^{\prime},\eta_{j}^{\prime}=0$ will make $A_{j}$ and $b_{j}$ ill-
defined completes our proof. ∎

Using Theorem 3, we can decide whether we need to solve the inequalities given
by $R_{1}$ or $R_{2}$ according to the value of $\delta_{j}$. Therefore, we
mitigate the computational complexity by solving $N$ sets of inequalities. By
observing that
$-\eta_{j}-zd_{j}\leq-\eta_{j}+z(\sup_{x_{j}}\\{h_{j}(x_{j})-d_{j})\\})$, we
further have that if there exist non-negative
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$ such
that $R_{2}$ is satisfied, then there must also exist some non-negative
solution to the inequalities in $R_{1}$. Finally, the sign of $\delta_{j}$ can
be used to reason whether interconnections improve the resilience.

######  Proposition 3.

Consider a subsystem $\mathcal{S}_{j}$ whose resilience index is given as
$(d_{j},\tau_{j},\phi_{j},\eta_{j})$ before being interconnected. If
$\delta_{j}\geq 0$, then there exists a resilience index
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$ such
that the set of inequalities given by $R_{1}$ is feasible. Furthermore, the
interconnections improve the resilience of $\mathcal{S}_{j}$ in the sense that

$d_{j}^{\prime}\leq
d_{j},\quad\tau_{j}^{\prime}=\tau_{j},\quad\phi_{j}^{\prime}\leq\phi_{j},\quad\eta_{j}^{\prime}\geq\eta_{j}.$

###### Proof.

We prove the proposition by giving a choice of non-negative
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$ that
satisfies $R_{1}$. We first note that $\tau_{j}^{\prime}=\tau_{j}\geq 0$ is
valid choice for parameter $\tau_{j}^{\prime}$. In this case, if
$\delta_{j}\geq 0$, we have that $d_{j}^{\prime}$ can be chosen as
$d_{j}^{\prime}=d_{j}+\delta_{j}\geq d_{j}\geq 0$. Given this choice of
$d_{j}^{\prime}$, we have that
$\phi_{j}^{\prime}=\frac{d_{j}^{\prime}}{d_{j}+\phi_{j}\delta_{j}}\phi_{j}\geq
0$. Since $\alpha$ is an extended class $\mathcal{K}$ function, we have that
$\alpha(h_{j}(x_{j})-d_{j}^{\prime})\geq\alpha(h_{j}(x_{j})-d_{j})$ for all
$x_{j}$. Therefore,
$\inf_{x_{j}\in\mathcal{D}_{j}}\\{\alpha(h_{j}(x_{j})-d_{j}^{\prime})-\alpha(h_{j}(x_{j})-d_{j})\\}\geq
0$. We can thus choose $\eta_{j}^{\prime}$ as
$\eta_{j}^{\prime}=\min\\{\frac{d_{j}}{\phi_{j}},\inf_{x_{j}\in\mathcal{D}_{j}}\\{\alpha(h_{j}(x_{j})-d_{j}^{\prime})-\alpha(h_{j}(x_{j})-d_{j})\\}\\}+\delta_{j}\geq
0$. Hence,
$(d_{j}^{\prime},\tau_{j}^{\prime},\phi_{j}^{\prime},\eta_{j}^{\prime})$ is a
valid resilience index for $\mathcal{S}_{j}$ after interconnection. ∎

##  VI Case Study

In this section, we demonstrate how the proposed resilience index can be used
to analyze safety constraints of interconnected systems.

We consider two well-mixed, nonisothermal continuous stirred-tank reactors
(CSTRs), denoted as $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$. We assume that
three parallel elementary irreversible exothermic reactions of the form
$A\xrightarrow{r_{1}}B$, $B\xrightarrow{r_{2}}E$, and $A\xrightarrow{r_{3}}Q$
occurs, where $A$ is the reactant species, $B$ is the desired product, and $E$
as well as $Q$ are undesired byproducts. The states of $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ are denoted as $x_{1}=[T_{1},c_{1}]^{\top}$ and
$x_{2}=[T_{2},c_{2}]^{\top}$, where $T_{i}$ and $c_{i}$ respectively represent
temperature of the reactor and concentration of $\mathcal{S}_{i}$ with
$i\in\\{1,2\\}$. Each CSTR utilizes a jacket to remove or provide heat to the
reactor to control the chemical reaction.

When the CSTRs are interconnected in series, their dynamics are given as

$\displaystyle\dot{x}_{1}$
$\displaystyle=\begin{bmatrix}\frac{F_{e,1}}{V_{1}}(T_{0,1}-T_{1})-\sum_{r=1}^{3}\frac{H_{r}}{\rho
p}R_{r}(c_{1},T_{1})\\\
\frac{F_{e,1}}{V_{i}}(c_{0,1}-c_{1})-\sum_{r=1}^{3}R_{r}(c_{1},T_{1})\end{bmatrix}+\begin{bmatrix}\frac{u_{1}}{\rho
pV_{1}}\\\ 0\end{bmatrix}$ (17) $\displaystyle\dot{x}_{2}$
$\displaystyle=\begin{bmatrix}\frac{F_{e,2}}{V_{2}}(T_{0,2}-T_{2})-\sum_{r=1}^{3}\frac{H_{r}}{\rho
p}R_{r}(c_{2},T_{2})\\\
\frac{F_{e,2}}{V_{2}}(c_{0,2}-c_{2})-\sum_{r=1}^{3}R_{r}(c_{2},T_{2})\end{bmatrix}+\begin{bmatrix}\frac{u_{1}}{\rho
pV_{1}}\\\ 0\end{bmatrix}$
$\displaystyle\quad\quad\quad+\frac{F_{1}}{V_{2}}\begin{bmatrix}T_{1}-T_{2}\\\
c_{1}-c_{2}\end{bmatrix}$ (18)

where $F_{e,i}$ is the flow rate,
$R_{r}(c_{i},T_{i})=k_{i}\exp(-E_{i}/lT_{i})c_{i}$, $V_{i}$ is the volume,
$u_{i}$ is the rate of heat input/ removal. We follow the choices of process
parameters given in [2] and summarize them in Table I. The set of admissible
inputs is set as $|u_{1}|\leq 2.7\times 10^{6}$ KJ/hr and $|u_{2}|\leq
2.8\times 10^{6}$ KJ/hr. In Eqn. (18), we have that

$W_{12}(x_{1},x_{2})=\frac{F_{1}}{V_{2}}\begin{bmatrix}T_{1}-T_{2}\\\
c_{1}-c_{2}\end{bmatrix}.$

Figure 1: This figure presents the temperature in CSTRs $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ after they are interconnected in series. The temperature of
$\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are represented in solid and dashed
lines. The portion plotted in red represents that the CSTR is faulty or
compromised, whereas the portion plotted in green denotes that the desired
rate of heat input is available.

We assume that the safety constraints defined for CSTRs $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ are $\mathcal{C}_{1}=\\{T_{1}:h_{1}(T_{1})\geq 0\\}$ and
$\mathcal{C}_{2}=\\{T_{2}:h_{2}(T_{2})\geq 0\\}$, where
$h_{i}(T_{i})=(T_{i}-300)(400-T_{i})$ for all $i\in\\{1,2\\}$. That is, the
temperature in both CSTR needs to be within range $[300,400]$K. Both CSTRs can
be faulty or compromised, leading to manipulated rate of heat input
$\tilde{u}_{i}$, where $i\in\\{1,2\\}$.

TABLE I: This table presents the values of process parameters used in Eqn. (17) and (18). Parameter Value | Unit
---|---
$F_{e,1}=4.998$, $F_{e,2}=30.0$, $F_{1}=4.998$ | $m^{3}/hr$
$V_{1}=1.0$, $V_{2}=3.0$ | $m^{3}$
$T_{0,1}=300.0$, $T_{0,2}=300.0$ | K
$E_{1}=5.0$, $E_{2}=7.53$, $E_{3}=7.53$ |  $\times 10^{4}$ KJ/kmol
$H_{1}=-5.0$, $H_{2}=-5.2$, $H_{3}=-5.4$ |  $\times 10^{4}$ KJ/kmol
$\rho=1000.0$ | kg/$m^{3}$
$k_{1}=3.0\times 10^{6}$, $k_{2}=3.0\times 10^{5}$, $k_{3}=3.0\times 10^{5}$ | hr-1
$l=8.314$ | KJ/kmol$\cdot$K
$p=0.231$ | KJ/kg$\cdot$ K
$c_{0,1}=4.0$, $c_{0,2}=2.0$ | kmol/$m^{3}$

We first compute the resilience indices of CSTRs $\mathcal{S}_{1}$ and
$\mathcal{S}_{2}$ when they are not interconnected. Their resilience indices
are given as $(d_{1},\tau_{1},\phi_{1},\eta_{1})=(2100,0.0146,0.308,0)$ and
$(d_{2},\tau_{2},\phi_{2},\eta_{2})=(500,0.0292,0.0222,0)$. Following dynamics
Eqn. (17) and (18), we have that the $\delta_{12}$ defined in Eqn. (5) is
negative, and hence we aim to solve $R_{2}$ to compute the resilience index of
CSTR $\mathcal{S}_{2}$ after interconnection. We have that
$(d_{2}^{\prime},\tau_{2}^{\prime},\phi_{2}^{\prime},\eta_{2}^{\prime})=(800,0.0237,0.1368,0)$.
We observe that by fixing $\tau_{2}=\tau_{2}^{\prime}$,
$d_{2}^{\prime}<d_{2}$, $\phi_{2}^{\prime}<\phi_{2}$, and
$\eta_{2}^{\prime}>\eta_{2}$, which aligns with our result in Proposition 3.
We simulate the temperature in both CSTRs in Fig. 1. We plot the time period
when the control input is compromised in red color, and the time period when
the desired control law is online in green color. We observe that the fault or
attack could manipulate the temperature in both CSTRs by changing the rate of
heat input $u_{i}$. We further demonstrate that safety constraints defined on
the temperature of $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are met. We plot
the boundaries of the safety set $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ using
dash-dotted blue lines. We observe that $T_{1}$ and $T_{2}$ remain within
$[300,400]$K for all time $t\geq 0$, and hence safety constraint is satisfied
if we can find a feasible resilience index, which demonstrates Theorem 2.

##  VII Conclusion

In this paper, we investigated the problem of efficient safety verification
for large-scale interconnected systems under faults and attacks. We developed
a compositional resilience index for each subsystem to characterize its
capability on tolerating faults and attacks without violating safety
constraints. We showed that if a subsystem possessed a resilience index, then
it satisfies the given safety constraint regardless of the faults and attacks.
We formulated a sum-of-squares optimization program to compute the resilience
index. When the resilience index and a safe control law of a subsystem were
given, we proved that the resilience index of the subsystem after being
interconnected could be computed by solving a system of linear inequalities.
We further developed the sufficient conditions over the interconnections to
guarantee the derived linear inequalities to be feasible. We demonstrated the
proposed approach using a case study on interconnected chemical reactors.

## References

  * [1] J. E. Sullivan and D. Kamensky, “How cyber-attacks in Ukraine show the vulnerability of the US power grid,” _The Electricity Journal_ , vol. 30, no. 3, pp. 30–35, 2017.
  * [2] N. H. El-Farra, A. Gani, and P. D. Christofides, “Fault-tolerant control of process systems using communication networks,” _AIChE Journal_ , vol. 51, no. 6, pp. 1665–1682, 2005.
  * [3] A. Greenberg, “Hackers remotely kill a Jeep on the highway–with me in it,” 2015\. [Online]. Available: https://www.wired.com/2015/07/hackers-remotely-kill-jeep-highway/
  * [4] C. Fan, K. Miller, and S. Mitra, “Fast and guaranteed safe controller synthesis for nonlinear vehicle models,” in _International Conference on Computer Aided Verification_. Springer, 2020, pp. 629–652.
  * [5] K. Koscher, S. Savage, F. Roesner, S. Patel, T. Kohno, A. Czeskis, D. McCoy, B. Kantor, D. Anderson, H. Shacham, and S. Savage, “Experimental security analysis of a modern automobile,” in _IEEE Symposium on Security and Privacy_. IEEE, 2010, pp. 447–462.
  * [6] Y. Chen, A. Singletary, and A. D. Ames, “Guaranteed obstacle avoidance for multi-robot operations with limited actuation: A control barrier function approach,” _IEEE Control Systems Letters_ , vol. 5, no. 1, pp. 127–132, 2020\.
  * [7] N. Fijalkow, J. Ouaknine, A. Pouly, J. Sousa-Pinto, and J. Worrell, “On the decidability of reachability in linear time-invariant systems,” in _22nd ACM International Conference on Hybrid Systems: Computation and Control_ , 2019, pp. 77–86.
  * [8] S. Prajna, A. Jadbabaie, and G. J. Pappas, “A framework for worst-case and stochastic safety verification using barrier certificates,” _IEEE Transactions on Automatic Control_ , vol. 52, no. 8, pp. 1415–1428, 2007.
  * [9] L. Niu, D. Sahabandu, A. Clark, and P. Radha, “Verifying safety for resilient cyber-physical systems via reactive software restart,” in _ACM/IEEE 13th International Conference on Cyber-Physical Systems (ICCPS)_. ACM/IEEE, 2022, pp. 104–115.
  * [10] A. Clark, “Verification and synthesis of control barrier functions,” in _60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 6105–6112.
  * [11] Z. Manna and A. Pnueli, _Temporal Verification of Reactive Systems: Safety_. Springer Science & Business Media, 2012.
  * [12] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games,” _IEEE Transactions on Automatic Control_ , vol. 50, no. 7, pp. 947–957, 2005.
  * [13] D. P. Bertsekas and I. B. Rhodes, “On the minimax reachability of target sets and target tubes,” _Automatica_ , vol. 7, no. 2, pp. 233–247, 1971.
  * [14] J. F. Fisac, A. K. Akametalu, M. N. Zeilinger, S. Kaynama, J. Gillula, and C. J. Tomlin, “A general safety framework for learning-based control in uncertain robotic systems,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 7, pp. 2737–2752, 2018.
  * [15] C. Sloth, G. J. Pappas, and R. Wisniewski, “Compositional safety analysis using barrier certificates,” in _Proceedings of the 15th ACM International Conference on Hybrid Systems: Computation and Control_ , 2012, pp. 15–24.
  * [16] Z. Lyu, X. Xu, and Y. Hong, “Small-gain theorem for safety verification of interconnected systems,” _Automatica_ , vol. 139, p. 110178, 2022.
  * [17] S. Coogan and M. Arcak, “A dissipativity approach to safety verification for interconnected systems,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 6, pp. 1722–1727, 2014.
  * [18] M. Anand, A. Lavaei, and M. Zamani, “From small-gain theory to compositional construction of barrier certificates for large-scale stochastic systems,” _arXiv preprint arXiv:2101.06916_ , 2021.
  * [19] E. M. Clarke, “Model checking,” in _International Conference on Foundations of Software Technology and Theoretical Computer Science_. Springer, 1997, pp. 54–56.
  * [20] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 8, pp. 3861–3876, 2016.
  * [21] X. Xu, “Constrained control of input–output linearizable systems using control sharing barrier functions,” _Automatica_ , vol. 87, pp. 195–201, 2018.
  * [22] J. Breeden and D. Panagou, “Guaranteed safe spacecraft docking with control barrier functions,” _IEEE Control Systems Letters_ , vol. 6, pp. 2000–2005, 2021.
  * [23] H. Fawzi, P. Tabuada, and S. Diggavi, “Secure estimation and control for cyber-physical systems under adversarial attacks,” _IEEE Transactions on Automatic control_ , vol. 59, no. 6, pp. 1454–1467, 2014.
  * [24] M. Pajic, J. Weimer, N. Bezzo, P. Tabuada, O. Sokolsky, I. Lee, and G. J. Pappas, “Robustness of attack-resilient state estimators,” in _ACM/IEEE International Conference on Cyber-Physical Systems (ICCPS)_ , 2014, pp. 163–174.
  * [25] L. Sha, “Using simplicity to control complexity,” _IEEE Software_ , vol. 18, no. 4, pp. 20–28, 2001.
  * [26] F. Abdi, C.-Y. Chen, M. Hasan, S. Liu, S. Mohan, and M. Caccamo, “Guaranteed physical security with restart-based design for cyber-physical systems,” in _2018 ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS)_. ACM/IEEE, 2018, pp. 10–21.
  * [27] J. S. Mertoguno, R. M. Craven, M. S. Mickelson, and D. P. Koller, “A physics-based strategy for cyber resilience of CPS,” in _Autonomous Systems: Sensors, Processing, and Security for Vehicles and Infrastructure 2019_ , vol. 11009. International Society for Optics and Photonics, 2019, p. 110090E.
  * [28] H. Yang, B. Jiang, M. Staroswiecki, and Y. Zhang, “Fault recoverability and fault tolerant control for a class of interconnected nonlinear systems,” _Automatica_ , vol. 54, pp. 49–55, 2015.
  * [29] H. Yang, C. Zhang, Z. An, and B. Jiang, “Exponential small-gain theorem and fault tolerant safe control of interconnected nonlinear systems,” _Automatica_ , vol. 115, p. 108866, 2020.
  * [30] A. Al Maruf, L. Niu, A. Clark, J. S. Mertoguno, and R. Poovendran, “A compositional approach to safety-critical resilient control for systems with coupled dynamics,” in _61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 910–917.
  * [31] F. Pasqualetti, F. Dörfler, and F. Bullo, “Attack detection and identification in cyber-physical systems,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 11, pp. 2715–2729, 2013.
  * [32] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex Optimization_. Cambridge university press, 2004.
  * [33] S.-C. Fang and S. Puthenpura, _Linear Optimization and Extensions: Theory and Algorithms_. Prentice-Hall, Inc., 1993\.


# The Laser Gravitational Compass

Ivan S. Ferreira
University of Brasilia, Institute of Physics, Brasilia, DF Brasilia, DF
70910-900<EMAIL_ADDRESS>
C. Frajuca
National Institute for Space Research, Sao Jose dos Campos, 12227-010,
<EMAIL_ADDRESS>
Nadja S. Magalhaes
Physics Department, Sao Paulo Federal University, SP09913-030, Brazil,
<EMAIL_ADDRESS>
M. D. Maia
University of Brasilia, Institute of Physics, Brasilia, DF70910-900,
<EMAIL_ADDRESS>
Claudio M. G. Sousa
Federal University of Para, Santarem, PA 68040-070<EMAIL_ADDRESS>
###### Abstract

Using the observational properties of Einstein’s gravitational field it is
shown that a minimum of four non-coplanar mass probes are necessary for the
Michelson and Morley interferometer to detect gravitational waves within the
context of General Relativity. With fewer probes, some alternative theories of
gravitation can also explain the observations. The conversion of the existing
gravitational wave detectors to four probes is also suggested.

##  1 The Observable Gravitational Wave

Einstein’s prediction of gravitational waves (gw) was originally derived from
arbitrarily small perturbations of the Minkowski metric
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}$, such that Einstein’s equations reduce
to a linear wave equation, written in a special (de Donder) coordinate gauge,
as (Except when explicitly stated, Greek indices run from 0 to 3 and small
case Latin indices run from 1 to 3).

$\Box^{2}\Psi_{\mu\nu}=0,\;\;\;\;\Psi_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}h\eta_{\mu\nu},\;\;\;\>h=\eta^{\mu\nu}h_{\mu\nu}.$
(1)

The currently operating laser gw observatories are inspired by the Michelson &
Morley (M&M) interferometers[1, 3, 2], where the data acquired by the Fabri-
Perót interferometer (an etalon) is used to generate a numerical simulation,
thus producing a template from which the most probable source is estimated[4].
The purpose of this note is to show that in order to detect gw described by
General Relativity (GR) with an M&M interferometer requires a minimum of four
non-coplanar mass probes.

The observables of Einstein’s gravitational field are given by the eigenvalues
of the Riemann curvature[5, 6, 7], defined by
$R(U,V)W=[\nabla_{U},\nabla_{V}]W-\nabla_{[U,V]}W$, whose components in any
given basis $\\{e_{\mu}\\}$ are
$R(e_{\mu},e_{\nu})e_{\rho}=R_{\mu\nu\rho\sigma}e^{\sigma}$. Then we find that
there are at most six independent eigenvectors $X_{\mu\nu}$ and six
eigenvalues $\lambda$, solutions of the eigenvalue equation
$R_{\mu\nu\rho\sigma}X^{\rho\sigma}=\lambda X_{\mu\nu}$, including the zero
eigenvalue, corresponding to the absence of gravitation[8, 9, 10]. Thus, using
the language of field theory, Einstein’s gravitation is said to have five non-
trivial observables or degrees of freedom $(dof)$. The spin-statistics theorem
relates the $dof$ to the helicity or the orbital spin of the field as
$s=(dof-1)/2$, so that Einstein’s gravitation is also said to be a spin-2
field.

Alternative theories of gravitation may have distinct definitions of
observables (not necessarily related to curvature) and their gravitational
waves, if they exist, may require different methods of observations. Well
known examples include: the spin-1 gauge theories of gravitation (there are
several of them), characterized by dof=3; topological gravitation in three
dimensions; projective theories of gravitation; F(R) theories; among many
others. Therefore, in order to understand the observation of a gw it is
essential to specify the observables of the theory on which the experiment is
based.

The most general massless spin-2 field $h_{\mu\nu}$ was defined in the
Minkowski space-time by Fierz and Pauli[11], as a trace-free field
$h=h_{\mu}{}^{\mu}=0$, satisfying the field equations $\Box^{2}h_{\mu\nu}=0$.
This is a linear field not to be confused with Einstein’s gravitation. Since
in the case of the present observations of gw, the supporting theory is
Einstein’s gravitation, then the observational signature to be sought is that
of a $dof=5$ or of a spin-2 field, characterized by the observable curvature.

The use of the Fierz-Pauli field $h_{\mu\nu}$ as the perturbation of the
Minkowski metric makes it possible to free Eq. (1) of coordinate gauges so
that its solutions can be written as as a superposition of plane polarized
gravitational waves, characterized by the Traceless-Transverse-Plane-Polarized
(TTPP) gauge conditions[1]:

$h=0,\;\;h_{i0}=0,\;\;\Box^{2}h_{\mu\nu}=0,\;\;h_{\mu\nu;\rho}=0.$ (2)

Then, these conditions are used to simulate a template, from which the source
source of the gw observations is estimated.

##  2 The Equivalence Principle and the M&M gw Detector

The use of an M&M detectors for gw is based on the principle of equivalence of
GR: Given 2 masses A and B with attached mirrors, under the exclusive action
of a known gravitational field, they propagate (or “free fall”) along time-
like geodesics, with unit tangent vectors $T_{A}$ and $T_{B}$ respectively,
satisfying geodesic equations $\nabla_{T_{A}}T_{A}=0$ and
$\nabla_{T_{B}}T_{B}=0$ . Eventually, probe A sends a light signal with
velocity $P$ to particle B along the light geodesic with equation
$\nabla_{P}P=0$. After a while, probe A receives back the reflected signal, so
that these geodesics describe a closed parallelogram, with the closing
condition $\nabla_{T}P=\nabla_{P}T$.

The curvature tensor calculated in that parallelogram is
$R(T,P)T=[\nabla_{T},\nabla_{P}]T=\nabla_{T}(\nabla_{T}P)=\nabla_{T}a$, where
$a=\nabla_{T}P$ is the acceleration of the signal along the fall. Defining a
basis by $\\{e_{0}=T,e_{i}=P\\}$, where $e_{i}$ denotes any space-like
direction, we obtain the geodesic deviation equation in components
$\frac{da_{i}}{cdt}=R_{0i0i}$, where $t$ denotes the time parameter of the
time-like geodesic. As we can see, the motion of the probes generates a
2-dimensional world-sheet, whose curvature $R_{0i0i}$, is translated directly
into the the variation of the fall acceleration.

In the currently operating and some future planned M&M gw detectors, three
mass-probes are used, defining a space-like plane in space-time, whose motion
under the action of a gw generates a 3-dimensional world-volume, whose
curvature is measured by a $2\times 2$ acceleration array $a_{ij}$, obeying
the geodesic deviation equation

$\frac{da_{ij}}{cdt}=\\!\\!R_{0i0j},\;i,j=1,2.$ (3)

Therefore, such detectors are capable to measure only 3 curvature components
$R_{0101},R_{0102},R_{0202}$, so that at most three degrees of freedom of
Einstein’s gravitational field can be obtained. Since the observed gw are very
weak, it has been assumed that the missing degrees of freedom in one detector
may be complemented by the data collected on another detector located
somewhere else. Such understanding is supported by the numerical simulation,
from which an estimate of the wave source capable of reproducing the same
observation in two separate detectors is obtained. Although it is possible to
parallelly transport the curvature tensor from one detector to another, the
left hand side content of Eq. (3) represents a locally measured quantity which
cannot be transported to another detector, under the penalty of breaching the
principle of equivalence.

The definitive solution for the missing dof problem can be obtained by a
direct measure of the curvature tensor using the geodesic deviation equations
Eq. (3), extended to $i,j=1,2,3$. One such detector, the “gravitational
compass” was conceived by P. Szekeres. It consists of four non-coplanar mass
probes connected by six dynamometers, to measure the six curvature
eigenvectors. By comparing three of these eigenvectors with the three
principal directions of curvature, the compass can points directly to the
source[12].

Szekeres’ gravitational compass was improved by F. Pirani[10], where instead
the dynamometers, the four mass-probes had attached mirrors and light sources.
Then, applying Eq. (3) again for $i,j=1,2,3$, all components of the curvature
can be read locally, defining the 4-dimensional gravitational field at the
center of mass of the four probes. The interesting fact is that four non-
coplanar masses define a virtual spheroidal surface in space-time, then the
natural oscillation modes of the spheroid give a direct measure of the
curvature eigenvalues.

The Szekeres-Pirani gravitational compass can be implemented in the presently
operating gw observatories, with the addition of a fourth mass probe not
belonging to the plane of the existing 2-dimensional M&M detectors, located
either in a tower or well, not necessarily having the same length as the
existing arms, thus defining a spheroidal gw detector.
Summarizing, the presently operating gw detectors are hindered by the fact
that they cannot detect directly the all degrees of freedom of Einstein’s
gravitational field at each point, thus leaving open the possibility of
alternative explanations for the gw. For the specific case of gravitation in
the sense of Einstein’s theory, we suggest the addition of a fourth mass to
the existing detectors, so that all degrees of freedom of gravitation are
directly measured from the five fundamental oscillation modes of a spheroid.

Acknowledgments: NSM and CF acknowledge FAPESP for its support to their
research through the thematic project 2013/26258-4.

## References

  * [1] A. Giazotto, Phys. Reports 182, 365 (1989).
  * [2] D. M. Eardley, D. L. Lee and A. P. Lightman, Phys.Rev. Lett. 30, 884 (1973).
  * [3] M. Pitkin, S. Reid, S. Rowan and J. Hough, Living Reviews Relativity 14, 5 (2011).
  * [4] B. P. Abbot et al. Phys. Rev. Lett. 116, 061102 (2016).
  * [5] V. D. Zakharov, Gravitational Waves in Einstein’s Theory (John Willey and Sons, New York, N.Y. 1972).
  * [6] F. A. E. Pirani, Acta Physica Polonica 15, 389 (1956).
  * [7] F. A. E. Pirani, Physical Review 105, 1089 (1957).
  * [8] F. A. E. Pirani, Gravitational radiation, in  Gravitation, An Introduction to Current Research, ed. L. Witten (John Willey and Sons, New York , 1962).
  * [9] R. K. Sachs, Gravitational Waves in General Relativity. Proc. Roy. Soc. A, 270, 103 (1962).
  * [10] F. A. E. Pirani, Introduction to gravitational radiation theory, in Lectures in General Relativity, Brandeis Summer Institute in Theoretical Physics, eds. S. Deser and K. W. Ford (Prentice-Hall, Englewood Cliffs, 1965).
  * [11] M. Fierz and W. Pauli, On Relativistic Wave Equations for Particles of Arbitrary Spin in an Electromagnetic Field. Proc. Roy. Soc. Lond. A 173, 211 (1939).
  * [12] P. Szekeres, The Gravitational Compass, Jour. Math. Phys. 6, 1387 (1965).
  * [13] N. S. Magalhaes, W. W. Johnson, C. Frajuca and O. D. Aguiar, Determination of astrophysical parameters from the spherical gravitational wave detector data. Mon. Not. R. Astron. Soc.274, 670 (1995).
  * [14] N. S. Magalhaes, W. W. Johnson, C. Frajuca and O. D. Aguiar, A geometric method for location of gravitational wave sources. Astrophys. Jour. 475, 462 (1997).