# An Analytical Framework for Control Synthesis of Cyber-Physical Systems with
Safety Guarantee

Luyao Niu1∗, Abdullah Al Maruf2∗, Andrew Clark1, J. Sukarno Mertoguno3, and
Radha Poovendran2 *Authors contributed equally to this work.1Luyao Niu and
Andrew Clark are with the Department of Electrical and Computer Engineering,
Worcester Polytechnic Institute, Worcester, MA 01609
<EMAIL_ADDRESS>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>Sukarno Mertoguno is with Information and Cyber Sciences Directorate, Georgia
Tech Research Institute, Atlanta, GA 30332<EMAIL_ADDRESS>
###### Abstract

Cyber-physical systems (CPS) are required to operate safely under fault and
malicious attacks. The simplex architecture and the recently proposed cyber
resilient architectures, e.g., Byzantine fault tolerant++ (BFT++), provide
safety for CPS under faults and malicious cyber attacks, respectively.
However, these existing architectures make use of different timing parameters
and implementations to provide safety, and are seemingly unrelated. In this
paper, we propose an analytical framework to represent the simplex, BFT++ and
other practical cyber resilient architectures (CRAs). We construct a hybrid
system that models CPS adopting any of these architectures. We derive
sufficient conditions via our proposed framework under which a control policy
is guaranteed to be safe. We present an algorithm to synthesize the control
policy. We validate the proposed framework using a case study on lateral
control of a Boeing 747, and demonstrate that our proposed approach ensures
safety of the system.

##  I Introduction

Cyber-physical systems (CPS) are subject to random failures and malicious
cyber attacks, which have been reported in applications such as transportation
[1] and power system [2]. Failures and attacks can potentially cause safety
violation of the physical components, which leads to severe harm to the plants
and humans.

Fault tolerant control schemes [3, 4, 5] and architectures such as simplex [6,
7, 8] have been proposed to address random failures. These approaches are
effective in CPS when some components are verified to be fault-free. This
requirement, however, may not be viable for all CPS, especially those subject
to malicious attacks.

A malicious adversary can exploit the vulnerabilities in the cyber subsystem
and intrude into CPS. The adversary can then cause common mode failures across
different components, rendering fault-tolerant schemes designed for random
failures inadequate. A seminal work recently proposed a cyber resilient
architecture, named Byzantine fault tolerant++ (BFT++) [9], for CPS under
malicious cyber attacks. BFT++, which is applied to CPS with redundant
controllers, uses one of the controllers as backup. The other controllers are
engineered to crash upon malicious attack, which triggers automatic and fast
controller recovery using the backup. If the controllers are restored in time,
then the system can guarantee safety. This architecture exploits the fact that
the cyber subsystem operates on a shorter timescale than the physical
subsystem, which has an inherent natural resilience from the physical dynamics
against limited cyber disruptions.

Following BFT++, several other unpublished yet effective approaches have
appeared with different implementations of recovery and backup. In parallel,
alternative approaches are proposed in [10, 11] for CPS without redundancy.
The controller in these approaches is programmed to restart proactively or
periodically to recover the system from malicious attack. The physical
subsystem can then maintain safety by utilizing its natural resilience and
tuning the controller availability.

The aforementioned architectures [9, 10, 11], which we collectively refer as
cyber resilient architectures (CRAs) have found successful applications in
different CPS. While the CRAs can independently provide safety guarantees, the
analyses undertaken are distinct and specific to the systems or architectures.
Hence these analyses may not be readily extended from one CPS to another.
Therefore a common framework which allows a general method of analysis for
these seemingly unrelated yet novel architectures is of key interest. Such a
framework will also enable comparison among different architectures under a
common baseline. Currently, such an analytical framework does not exist.

In this paper, we propose a common framework that models the simplex
architecture and the CRAs. We then present a control policy synthesis with
safety guarantee using our proposed framework, which applies to any of these
architectures. We make the following specific contributions:

  * •

We construct a hybrid system to model CPS implementing the simplex
architecture and CRAs. We propose a common framework that captures these
architectures.

  * •

We derive the sufficient conditions for a control policy to satisfy safety
with respect to any specified budget.

  * •

We propose an algorithm to compute a control policy that satisfies our derived
conditions. Our proposed algorithm converges to a feasible solution, given its
existence, within finite number of iterations.

  * •

We validate our proposed approach using a case study on lateral control of a
Boeing 747. We show that our proposed approach guarantees the safety of Boeing
747 with respect to the given budget constraint.

The reminder of the paper is organized as follows. Section II presents the
related work. Section III introduces the CPS model and presents the problem
statement. Section IV gives our proposed framework. Section V presents our
proposed solution approach. Section VI contains a case study on a Boeing 747.
Section VII concludes the paper.

##  II Related Work

Safety verification [12, 13] and safety controller synthesis [14, 15, 16, 17]
for CPS operated in benign environment have been extensively studied.

Fault tolerant controllers [3, 4, 5] and architectures [6, 7, 8] have been
widely adopted for CPS that may incur faults. One of the well-known fault
tolerant designs is the simplex architecture. This architecture consists of a
main controller which is vulnerable to random failures and a safety controller
which is verifiable and fault-free [6]. Under certain conditions, e.g., the
main controller experiences a fault, a decision module instantaneously
switches to the safety controller. The decision module switches back to the
main controller after the main controller recovers from fault. These fault
tolerant approaches assume that there is no common failure for all components
which may not hold for malicious attack.

There exist two main trends of approaches to address CPS under malicious cyber
attacks. The first category aims at protecting the system from malicious
attacks using control- and game-theoretic approaches [18, 19, 20]. These
approaches detect the attack and then filter its impact. The second body of
literature focuses on designing attack tolerant systems. The CRAs [9, 10, 11,
21] belong to this category.

BFT++ and other variants [9] are applied when CPS have redundant controllers.
One of the redundant controllers is used as backup and equipped with a buffer
storing the time-delayed inputs. The non-backup controllers are deliberately
engineered to crash following a malicious exploit, e.g., by implementing
software diversity [22] or memory/instruction randomization [23]. Sensing the
crash, BFT++ recovers the controllers quickly from the backup whose integrity
is ensured by flushing the buffer.

The CRA proposed in [10] and also the restart-based mechanisms [24, 25, 26,
27, 21] are applicable to CPS that do not have redundancy. These approaches
reset the cyber subsystem to a ‘clean’ state via restart to recover from
attack. The authors of [10] tunes controller availability for the safety of
physical subsystem, whereas the restart-based mechanisms use reachability
analysis [25, 26, 27] for safety guarantee.

##  III System Model and Problem Formulation

We first give some notations before presenting the system model. Then we state
the problem investigated in this paper.

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$. Throughout this paper, we use $\mathbb{R}$,
$\mathbb{R}_{\geq 0}$, $\mathbb{R}_{>0}$, and $\mathbb{Z}_{\geq 0}$ to denote
the set of real numbers, non-negative real numbers, positive real numbers, and
non-negative integers, respectively. Given a vector $x\in\mathbb{R}^{n}$, we
denote its $i$-th entry as $[x]_{i}$, where $i=1,\ldots,n$.

Consider a CPS consisting of a cyber subsystem and a physical subsystem. The
physical subsystem is modeled by a plant that evolves following

$\dot{x}_{t}=f(x_{t})+g(x_{t})u_{t},$ (1)

where $x_{t}\in\mathcal{X}\subset\mathbb{R}^{n}$ is the system state and
$u_{t}\in\mathcal{U}\subset\mathbb{R}^{m}$ is the control input. Functions
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ and
$g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times m}$ are assumed to be
Lipschitz continuous. We also assume that
$\mathcal{U}=\prod_{i=1}^{m}[u_{i,min},u_{i,max}]$ with $u_{i,min}<u_{i,max}$.
The physical plant is normally recommended to be operated within a certain
range $\mathcal{C}=\\{x\in\mathcal{X}:h(x)\geq 0\\}$, where
$h:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a continuously differentiable
function. We assume that set $\mathcal{C}$ is compact. Given the system state
$x$, the actuator signal $u$ is determined by a control policy
$\mu:\mathcal{X}\rightarrow\mathcal{U}$.

Although the physical plant evolves in continuous time, the cyber subsystem
interacts with the physical subsystem following functioning cycles. We assume
that the sensors can directly measure the physical state $x$. At each
functioning cycle $k\in\mathbb{Z}_{\geq 0}$, the cyber subsystem measures
$x_{k\delta}$ and updates the actuator signal $u_{k\delta}=\mu(x_{k\delta})$.
The actuator signal remains constant during each functioning cycle $k$. In the
remainder of this paper, we refer to a functioning cycle as an epoch with
length $\delta>0$.

The system is subject to a malicious attack initiated by an intelligent
adversary. The adversary aims at driving the physical plant outside
$\mathcal{C}$ to damage it. The adversary can exploit the vulnerabilities in
the cyber subsystem and intrude into the system. Once the adversary intrudes
successfully, it gains access to the software, actuators, and other
peripherals. As a consequence, the actuator signal is corrupted by the
adversary and deviates from the desired control policy $\mu(\cdot)$. To
recover the system from attack, the CRAs and other mechanisms have been
proposed, as reviewed in Section I and II. Let $t_{1}\geq 0$ be a time instant
when the adversary corrupts the cyber subsystem. The CRAs eliminate the
adversary from the system at some time $\tilde{t}>t_{1}$. We denote the time
instant when the adversary successfully corrupts the cyber subsystem again for
the first time after $\tilde{t}$ as $t_{2}>\tilde{t}$. We define the interval
$[t_{1},t_{2}]$ as an _attack cycle_. Note that the length $A=t_{2}-t_{1}$ of
each attack cycle varies, and is dependent on the adversary. Later in Section
V, we will compute a lower bound for $A$ to guarantee system safety.

Due to the malicious attack, the physical plant may have to be temporarily
operated outside $\mathcal{C}$. To avoid causing irreversible damage to the
plant, we need to minimize the amount of time that the CPS is operated outside
$\mathcal{C}$ or minimize how far the physical state $x$ deviates from
$\mathcal{C}$. We capture the instantaneous damage incurred by the plant when
operated outside $\mathcal{C}$ as a cost
$L:\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$, defined as

$L(h(x))=\begin{cases}L_{1}(-h(x)),&\mbox{ if }h(x)<0\\\ 0,&\mbox{ if
}h(x)\geq 0\end{cases}$ (2)

where $L_{1}:\mathbb{R}_{>0}\rightarrow\mathbb{R}_{\geq 0}$ is a monotone non-
decreasing function. When $L_{1}(-h(x))=1$, then $\int_{t}L(x)\,\textrm{d}{t}$
is equal to the amount of time such that $x_{t}\notin\mathcal{C}$. When
$L_{1}(-h(x))=-h(x)$ for all $x\notin\mathcal{C}$, Eqn. (2) models the
deviation of the physical plant from the boundary of $\mathcal{C}$. We define
the physical safety with respect to budget $B$ as follows.

######  Definition 1 (Physical Safety with Respect to Budget $B$).

The physical plant is safe with respect to a budget $B$ if the following
relation holds for any attack cycle $[t_{1},t_{2}]$:

$J=\int_{t=t_{1}}^{t_{2}}L(h(x_{t}))\,\textrm{d}{t}\leq B.$ (3)

Eqn. (3) enforces an upper bound on the cost incurred by the system during any
attack cycle. When $B=0$, Definition 1 recovers the strict safety constraint
$x_{t}\in\mathcal{C}$ for all $t\geq 0$ as a special case. Given Definition 1,
the problem of synthesizing a control policy with safety guarantee is stated
as follows:

######  Problem 1.

Synthesize a control policy $\mu:\mathcal{X}\rightarrow\mathcal{U}$ for the
CPS such that Definition 1 is satisfied for a given budget $B$.

##  IV Our Proposed Cyber Resilient Framework

In this section, we first detail the timing behaviors of the CRAs. Then we
construct a hybrid system that models CPS adopting any of these architectures.
The simplex architecture reviewed in Section II is also incorporated in our
framework for completeness. We finally reformulate Problem 1 using the
constructed hybrid system.

###  IV-A Timing Behaviors of the CRAs

In this subsection, we present the timing behaviors of the CRAs [9, 10, 11],
which track the status of the cyber subsystem. Note that the status of the
cyber subsystem are discrete. When the adversary intrudes into the system at
epoch $j$, the cyber subsystem changes from the normal to the corrupted
status, indicating the adversary can arbitrarily manipulate the controllers.
If the CPS have redundant controllers as discussed in [9], the non-backup
controllers will crash by epoch $j+N_{1}$ in the worst-case, which triggers
controller restoration, leading the cyber subsystem to transit from the
corrupted status to the restoration status. In practical implementations of
BFT++, we observe that $N_{1}=2$ and the buffer length is chosen to be greater
than $N_{1}$. Denote the worst-case number of epochs needed for controller
restoration as $N_{2}$. Then the cyber subsystem returns to the normal status
using the backup controller by epoch $j+N_{1}+N_{2}$. To ensure safety of the
physical subsystem, crash delay $N_{1}$ needs to be small enough, and the
restoration time $N_{2}$ needs to be tolerated by the physical subsystem’s
resilience $\Delta(x)$ which is determined by the physical state and system
dynamics.

When CPS have no redundant controllers, cyber subsystem recovery can be
triggered by either the attack or the timer [10, 11]. If the controller
crashes due to attack, which takes at most $N_{3}$ epochs, then the system
reboots and re-initializes the controller to recover it. The worst-case time
needed for such controller recovery, denoted as $N_{4}$, is in general larger
than BFT++, i.e., $N_{4}\geq N_{2}$. We remark that the controller restart is
triggered by the timer when the attack does not crash the controller. If the
recovery is triggered by the timer and there is no attack, then the system
restart is executed every $N_{4}+N_{5}$ epochs, where $N_{5}$ is the number of
epochs elapsed in the normal status during one restart period. Safety of the
physical subsystem is then ensured by tuning the time between two consecutive
restarts and the controller availability.

Figure 1: Hybrid system $H=(\mathcal{X},\mathcal{U},\mathcal{L},\mathcal{Y},\mathcal{Y}_{0},Inv,\mathcal{F},\allowbreak\Sigma,\allowbreak\mathcal{E},\allowbreak\Phi)$ captures the simplex architecture (left part) and the CRAs (right part). Each state $y=(x,(l_{1},j))\in\mathcal{Y}$ of the hybrid system captures the continuous physical state $x\in\mathcal{X}$ and the discrete location $l=(l_{1},j)\in\mathcal{L}$ including the system status $l_{1}$ and the epoch index $j$. Each circle in the figure represents a discrete location in $\mathcal{L}$, with the epoch indices being omitted. The arrows in the figure represent the transitions $\Sigma$ of hybrid system $H$. Each transition is labeled using $e\in\mathcal{E}$, and is enabled when the corresponding clock constraint $\phi\in\Phi$ is satisfied. The detailed labels and clock constraints are given in Table I. The transition in red color is triggered by external event, i.e., the adversarial cyber disruption.  Designs | Starting state $(x,(l_{1},j))$ | Target state $(x^{\prime},(l_{1}^{\prime},j^{\prime}))$ | Label $e=(\gamma,e_{2})$ |  Clock constraint $\phi$
---|---|---|---|---
Simplex architecture | $(x,(normal,j))$ | $(x,(SC,j^{\prime}))$ | $(Safety~{}controller~{}activated,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}=0)$
$(x,(SC,j))$ | $(x,(normal,j^{\prime}))$ | $(Safety~{}controller~{}inactivated,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}=0)$
the CRAs | $(x,(normal,j))$ | $(x^{\prime},(corrupted,j^{\prime}))$ | $(Adversarial~{}cyber~{}disruption,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\geq 0)$
$(x,(corrupted,j))$ | $(x^{\prime},(restoration,j^{\prime}))$ | $(Controller~{}crash,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\leq N_{1})$
$(x,(restoration,j))$ | $(x^{\prime},(normal,j^{\prime}))$ | $(Controller~{}restored,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\leq N_{2})$
$(x,(normal,j))$ | $(x^{\prime},(R\&I,j^{\prime}))$ | $(Timer,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}=N_{5})$
$(x,(corrupted,j))$ | $(x^{\prime},(R\&I,j^{\prime}))$ | $(Controller~{}crash,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\leq N_{3})$
$(x,(corrupted,j))$ | $(x^{\prime},(R\&I,j^{\prime}))$ | $(Timer,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\leq N_{5})$
$(x,(R\&I,j))$ | $(x^{\prime},(normal,j^{\prime}))$ | $(Initialization~{}done,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\leq N_{4})$
| $(x,(\cdot,j))$ | $(x^{\prime},(safety~{}violation,j^{\prime}))$ | $(Maximum~{}tolerance,e_{2})$ |  $(j+e_{2}=j^{\prime})\land(e_{2}\geq\Delta(x))$
TABLE I: This table shows the transitions with their corresponding labels and
clock constraints. The second and third columns give the starting and end
states of a transition, respectively. The fourth column presents the label
$e=(\gamma,e_{2})$ associated with the transition
$\sigma=((x,(l_{1},j)),(x^{\prime},(l_{1}^{\prime},j^{\prime})))$. The trigger
event of the transition is denoted as $\gamma$, and the time elapsed in
$l_{1}$ is denoted using $e_{2}\in\mathbb{Z}_{\geq 0}$. The fifth column gives
the clock constraint $\phi$ that needs to be satisfied by the transition
$\sigma$ and parameter $e_{2}$. Parameters $N_{1},N_{2},N_{3},N_{4}$, and
$N_{5}$ are determined by the architecture design. Parameter $\Delta(x)$
captures the maximum tolerance provided by the physical subsystem.

###  IV-B Our Proposed Framework

In this subsection, we first construct a hybrid system to model the CPS
implementing the simplex architecture or CRAs. We then restate Problem 1 in
the context of the hybrid system. We construct a hybrid system
$H=(\mathcal{X},\mathcal{U},\mathcal{L},\mathcal{Y},\mathcal{Y}_{0},Inv,\mathcal{F},\allowbreak\Sigma,\mathcal{E},\Phi)$,
as shown in Fig. 1, where

  * •

$\mathcal{X}\subseteq\mathbb{R}^{n}$ is the continuous state space modeling
the states of the physical subsystem. $\mathcal{U}\subseteq\mathbb{R}^{m}$ is
the set of admissible control inputs of the physical subsystem.

  * •

$\mathcal{L}=\\{normal,R\&I,SC,restoration,\allowbreak corrupted,\allowbreak
safety~{}violation\\}\times\mathbb{Z}_{\geq 0}$ is a set of discrete
locations111Throughout this paper, we denote $l=reboot\&initialization$ as
$l=R\&I$ and denote $l=safety~{}controller~{}driven$ as $l=SC$ for
simplicity., with each location $l\in\mathcal{L}$ modeling the status of the
system at each epoch index.

  * •

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

  * •

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

  * •

$\mathcal{F}$ is the set of vector fields. For each $F\in\mathcal{F}$, the
continuous system state evolves as $\dot{x}=F(x,u,l)$, where $F$ is jointly
determined by the system dynamics and the status of the cyber subsystem, and
$\dot{x}$ is the time derivative of continuous state $x$.

  * •

$\Sigma\subseteq\mathcal{Y}\times\mathcal{Y}$ is the set of transitions
between the states of the hybrid system. A transition
$\sigma=((x,l),(x^{\prime},l^{\prime}))$ models the state transition from
$(x,l)$ to $(x^{\prime},l^{\prime})$.

  * •

$\mathcal{E}=\Gamma\cup\mathbb{Z}_{\geq 0}$ is a set of labels, where $\Gamma$
is the finite alphabet set. Each $\gamma\in\Gamma$ is labeled on some
transition $\sigma\in\Sigma$ indicating the events that triggers the
transition.

  * •

$\Phi$ is a set of clock constraints, with each $\phi\in\Phi$ being defined as
$\phi:\Sigma\times\mathbb{Z}_{\geq 0}\rightarrow\\{0,1\\}$. Function $\phi$
maps the time elapsed in each discrete location labeled on each transition to
the binary set $\\{0,1\\}$, indicating if the transition is enabled or not.

In Fig. 1, we only label the first element of each location $l=(l_{1},j)$,
where $l_{1}\in\\{normal,R\&I,SC,restoration,\allowbreak corrupted,\allowbreak
safety~{}violation\\}$, and the epoch index $j\in\mathbb{Z}_{\geq 0}$ is
omitted. Variable $l_{1}$ represents the status of the system, as explained in
Section IV-A. Particularly, we use $R\&I$ to represent controller reboot and
initialization for CPS without redundancy, and use $SC$ to represent the
status where the system is driven by the safety controller as suggested in the
simplex architecture. In the remainder of this paper, we refer to $l_{1}$ as
the location of $H$ omitting the epoch index when the context is clear. The
set of vector fields $\mathcal{F}$ captures the dynamics at each discrete
location. For instance, when $l=(R\&I,j)$, we have that
$\dot{x}=F(x,u,l)=f(x)$ for all $j\in\mathbb{Z}_{\geq 0}$ since $u_{t}=0$.

We label each transition $\sigma=((x,l),(x^{\prime},l^{\prime}))\in\Sigma$
using $e=(\gamma,e_{2})\in\mathcal{E}$. The detailed label associated with
each transition can be found in Table I. Here we use $\gamma$ to represent the
event that triggers the transition. For instance, the transition from $normal$
to $corrupted$ is triggered by the adversarial cyber disruption, whereas the
transition from $corrupted$ to $restoration$ is triggered by controller crash.
The element $e_{2}\in\mathbb{Z}_{\geq 0}$ denotes the number of epochs elapsed
in status $l_{1}$ before the occurrence of transition
$\sigma=((x,(l_{1},j)),(x^{\prime},(l_{1}^{\prime},j^{\prime})))$. For
example, at most $N_{1}$ epochs elapse at location $corrupted$ before the
transition from $(x,(corrupted,j))$ to $(x^{\prime},(restoration,j^{\prime}))$
occurs. When $e_{2}=0$, it indicates that the transition occurs
instantaneously.

A transition in hybrid system $H$ is enabled if and only if a clock constraint
$\phi$ associated with the transition is satisfied. Consider a transition
$\sigma=((x,(l_{1},j)),(x^{\prime},(l_{1}^{\prime},j^{\prime})))$ labeled with
$e=(\gamma,e_{2})$. A clock constraint $\phi$ verifies if $j$, $j^{\prime}$,
and $e_{2}$ satisfies $j+e_{2}=j^{\prime}$ with $e_{2}$ being determined by
the architecture. The clock constraint enables hybrid system $H$ to always
track the correct epoch index.

We remark that for each status $l_{1}$, we do not depict the transitions
$(x,(l_{1},j))$ to $(x,(l_{1},j+1))$ for compactness of the figure. These
transitions do not cause any system status jump, and only track the evolution
of epoch indices. We also define a guard set $\mathcal{G}(l,l^{\prime})$ as
$\mathcal{G}(l,l^{\prime})=\\{x\in\mathcal{X}:((x,l),(x,l^{\prime}))\in\Sigma\\}$
which represents the set of physical states starting from which the system can
transit from location $l$ to $l^{\prime}$ on the hybrid system $H$.

The transitions that end at $l_{1}=safety~{}violation$ are triggered by event
$\gamma=Maximum~{}tolerance$, indicating that the physical subsystem has not
received the correct input in time and has utilized all resilience against the
disruption. In this case, safety violation $J>B$ becomes inevitable, and the
system cannot recover after safety violation (captured via the self-loop in
Fig. 1). We capture the maximum tolerance provided by the physical subsystem
using $e_{2}$ labeled on the transitions. Note that the tolerance depends on
the physical system state $x$ and is denoted as $\Delta(x)$.

We are now ready to translate problem 1 using the context of hybrid system $H$
as follows:
Restatement of Problem 1. Given hybrid system $H$, synthesize a control policy
$\mu:\mathcal{X}\rightarrow\mathcal{U}$ such that hybrid system $H$ never
reaches status $l_{1}=safety~{}violation$.

##  V Analysis of the Proposed Framework

This section presents the proposed solution approach to Problem 1. We first
develop sufficient conditions for the control policy that guarantees safety of
the physical subsystem under a cyber attack. Then we formulate the derived
conditions as sum-of-squares (SOS) constraints and propose an algorithm to
compute a control policy and the corresponding parameters. We finally give the
convergence and complexity of our algorithm.

In the following, we derive the sufficient conditions under which a control
policy $\mu:\mathcal{X}\rightarrow\mathcal{U}$ ensures the system to satisfy
Definition 1 with respect to a given budget $B$. The idea is that if control
policy $\mu$ drives the physical subsystem to $\mathcal{C}_{1}=\\{x:h(x)\geq
c\\}\subseteq\mathcal{C}$ when $l_{1}=normal$ and we can constrain the system
trajectory to remain in a set
$\mathcal{D}=\\{x:h(x)\geq-d\\}\supseteq\mathcal{C}$ for any
$l_{1}\in\\{corrupted,R\&I,restoration\\}$, then we can limit the worst-case
cost incurred during one attack cycle to be bounded by $B$ by tuning choices
$c,d\geq 0$. We denote the worst-case number of epochs when the system is at
some status $l_{1}\in\\{corrupted,R\&I,restoration\\}$ as $N$. We then have
the following conditions:

######  Theorem 1.

Consider hybrid system $H$ and let set $\mathcal{C}$ be defined as in Section
III. Let $h_{c}(x)=h(x)-c$ and $h_{d}(x)=h(x)+d$. We define
$\mathcal{C}_{1}=\\{x:h_{c}(x)\geq 0\\}$ and $\mathcal{D}=\\{x:h_{d}(x)\geq
0\\}$. Consider an arbitrary attack cycle denoted as $[t_{1},t_{2}]$ and
suppose $x_{t_{1}}\in\mathcal{C}_{1}$. If there exist constants $c,d\geq 0$,
$\tau>0$, and a control policy $\mu:\mathcal{X}\rightarrow\mathcal{U}$ such
that

$\displaystyle\frac{\partial h_{d}}{\partial
x}(x)(f(x)+g(x)u)\geq-\frac{c+d}{N\delta},~{}\forall(x,u)\in\mathcal{D}\times\mathcal{U}$
(4a) $\displaystyle\frac{\partial h_{c}}{\partial
x}(x)(f(x)+g(x)\mu(x))\geq\frac{c+d}{\tau},~{}\forall
x\in\mathcal{D}\setminus\mathcal{C}_{1}$ (4b) $\displaystyle\frac{\partial
h_{c}}{\partial x}(x)(f(x)+g(x)\mu(x))\geq-\alpha(h_{c}(x)),~{}\forall
x\in\mathcal{C}_{1}$ (4c)
$\displaystyle\frac{N\delta}{c+d}\int_{s=0}^{d}L_{1}(s)\,\textrm{d}{s}+\frac{\tau}{c+d}\int_{s=0}^{d}L_{1}(s)\,\textrm{d}{s}\leq
B$ (4d)

then system (1) is safe with respect to budget $B$ by taking policy $\mu$ at
$l_{1}=normal$, provided that $A=t_{2}-t_{1}\geq\tau+N\delta$. Furthermore,
$x_{t}\in\mathcal{D}$ for $t\in[t_{1},t_{2}]$ and $x_{t}\in\mathcal{C}_{1}$
for $t\in[t_{1}+\tau+N\delta,t_{2}]$.

###### Proof.

The proof consists of two steps. We first find the guard set for each state of
hybrid system $H$ when the conditions in Eqn. (4) hold. We then prove that the
system is safe with respect to budget $B$ according to Definition 1.

In the first step, we show that if $t_{2}\geq t_{1}+\tau+N\delta$,
$x_{t_{1}}\in\mathcal{C}_{1}$ and the conditions in (4) hold, then
$x_{t}\in\mathcal{D}$ for $t\in[t_{1},t_{2}]$. Additionally we show that
$x_{t}\in\mathcal{C}_{1}$ for $t\in[t_{1}+\tau+N\delta,t_{2}]$ using control
policy $\mu$ at $l_{1}=normal$. As hybrid system $H$ takes $N$ epochs to
transition to $l_{1}=normal$ after being corrupted, therefore $t_{1}+N\delta$
is the time instant when transition $\sigma=(\cdot,normal)$ takes place after
$t_{1}$. Then for any $t^{\prime}\in[t_{1},t_{1}+N\delta]$ and for any
$u\in\mathcal{U}$, we have that

$\displaystyle h(x_{t^{\prime}})=h(x_{t_{1}})$
$\displaystyle+\int_{t=t_{1}}^{t^{\prime}}\dot{h}\,\textrm{d}{t}\geq
c-\frac{c+d}{N\delta}(t^{\prime}-t_{1})\geq-d,$ (5)

where the inequality holds by Eqn. (4a) and the assumptions that
$x_{t_{1}}\in\mathcal{C}_{1}$ and $t^{\prime}\in[t_{1},t_{1}+N\delta]$.
Therefore, we have that $x_{t^{\prime}}\in\mathcal{D}$ for all
$t^{\prime}\in[t_{1},t_{1}+N\delta]$, indicating that
$\mathcal{G}(R\&I,normal),\mathcal{G}(restoration,normal)\subseteq\mathcal{D}$.

Now consider that control policy $\mu$ is applied when the system is at
$l_{1}=normal$. Let $\hat{t}$ be any time when the system is at $l_{1}=normal$
for which the trajectory remains in $\mathcal{D}\setminus\mathcal{C}_{1}$
(i.e. $-d\leq h(x_{\hat{t}})<-c)$. Then we can write

$h(x_{\hat{t}})=h(x_{t_{1}+N\delta})+\int_{t=t_{1}+N\delta}^{\hat{t}}\dot{h}\,\textrm{d}{t}\\\
\geq-d+\frac{c+d}{\tau}(\hat{t}-t_{1}-N\delta),$ (6)

where the inequality holds by Eqn. (4b) and the fact that
$\mathcal{G}(R\&I,normal),\mathcal{G}(restoration,normal)\subseteq\mathcal{D}$.
If $\hat{t}\geq t_{1}+N\delta+\tau$, then Eqn. (6) renders $h(x_{\hat{t}})\geq
c$ and thus $x_{\hat{t}}\in\mathcal{C}_{1}$. Further note that by [14, Thm.
2], $\mathcal{C}_{1}$ is forward invariant using control policy $\mu$. Thus,
$x_{t}\in\mathcal{C}_{1},~{}\forall t\in[t_{1}+N\delta+\tau,t_{2}]$. Since
$\mathcal{C}_{1}\subseteq\mathcal{D}$, we have $x_{t}\in\mathcal{D}$ for all
$t\in[t_{1},t_{2}]$.

In the second step, we quantify the worst-case cost incurred by the system. We
first compute the cost incurred during time $[t_{1},t_{1}+N\delta]$. Suppose
the system reaches the boundary of $\mathcal{C}$ at time instants $z_{1}\leq
z_{2}\leq\ldots\leq z_{K}\leq t_{1}+N\delta$, where $z_{1}\geq t_{1}$ and
$z_{K}\leq t_{1}+N\delta$. We have that

$J_{1}\doteq\int_{t=t_{1}}^{t_{1}+N\delta}L(h(x_{t}))\,\textrm{d}{t}=\int_{t=t_{1}}^{z_{1}}L(h(x_{t}))\,\textrm{d}{t}\\\
+\sum_{j=1}^{K-1}\int_{t=z_{j}}^{z_{j+1}}L(h(x_{t}))\,\textrm{d}{t}+\int_{t=z_{K}}^{t_{1}+N\delta}L(h(x_{t}))\,\textrm{d}{t}$

Since $h(x_{t})\in\mathcal{C}$ for $t\in[t_{1},z_{1}]$, therefore
$L(h(x_{t}))=0$. Additionally, since $L_{1}$ is monotone non-decreasing, by
Eqn. (5) we have that $L(h(x_{t}))=L_{1}(-h(x_{t}))\leq
L_{1}(-c+\frac{c+d}{N\delta}(t-t_{1}))$ for any $t\in[z_{j},z_{j+1}]$ if
$h(x_{t})\notin\mathcal{C}$. Using these arguments, we have that

$\displaystyle J_{1}$
$\displaystyle\leq\sum_{j=1}^{K-1}\int_{t=z_{j}}^{z_{j+1}}L_{1}(-c+\frac{c+d}{N\delta}t)\,\textrm{d}{t}$
$\displaystyle\quad\quad\quad\quad\quad\quad+\int_{t=z_{K}}^{t_{1}+N\delta}L_{1}(-c+\frac{c+d}{N\delta}t)\,\textrm{d}{t}$
$\displaystyle=\int_{t=z_{1}}^{t_{1}+N\delta}L_{1}(-c+\frac{c+d}{N\delta}t)\,\textrm{d}{t}.$

Using Eqn. (5), it follows that $z_{1}\geq
t_{1}+c/(\frac{c+d}{N\delta})=t_{1}+\frac{cN\delta}{c+d}$. Therefore, we have
that

$\displaystyle J_{1}$
$\displaystyle\leq\int_{t=t_{1}+\frac{cN\delta}{c+d}}^{t_{1}+N\delta}L_{1}(-c+\frac{c+d}{N\delta}(t-t_{1}))\,\textrm{d}{t}$
$\displaystyle=\int_{t=0}^{\frac{dN\delta}{c+d}}L_{1}(\frac{c+d}{N\delta}t)\,\textrm{d}{t}=\frac{N\delta}{c+d}\int_{s=0}^{d}L_{1}(s)\,\textrm{d}{s},$

where the above holds by variable substitution and the fact that $L_{1}$ is
non-negative.

We now quantify the worst-case cost incurred during time
$[t_{1}+N\delta,t_{2}]$. By Eqn. (6),
$h(x_{t})\geq-d+\frac{c+d}{\tau}(t-t_{1}-N\delta)$ for
$x_{t}\in\mathcal{D}\setminus\mathcal{C}_{1}$. Note that $h(x_{t})\geq 0$ for
all $t\in[t_{1}+N\delta+\frac{d\tau}{c+d},t_{2}]$ using control policy
$\mu(x)$. Therefore we have that

$\displaystyle J_{2}$
$\displaystyle\doteq\int_{t=t_{1}+N\delta}^{t_{2}}L(h(x_{t}))\,\textrm{d}{t}$
$\displaystyle=\int_{t=t_{1}+N\delta}^{t_{1}+N\delta+\frac{d\tau}{c+d}}L(h(x_{t}))\,\textrm{d}{t}+\int_{t=t_{1}+N\delta+\frac{d\tau}{c+d}}^{t_{2}}L(h(x_{t}))\,\textrm{d}{t}$
$\displaystyle\leq\int_{t=N\delta}^{N\delta+\frac{d\tau}{c+d}}L_{1}(d-\frac{c+d}{\tau}(t-N\delta))\,\textrm{d}{t}$
$\displaystyle=\int_{t=0}^{\frac{d\tau}{c+d}}L_{1}(d-\frac{c+d}{\tau}t)\,\textrm{d}{t}=\frac{\tau}{c+d}\int_{s=0}^{d}L_{1}(s)\,\textrm{d}{s}$

where the above holds by (6), $L_{1}$ is monotone non-decreasing and non-
negative, and $L(h(x_{t}))=0$ for $x_{t}\in\mathcal{C}$. Therefore we have
that the total cost is upper bounded by $J_{1}+J_{2}$, which yields condition
(4d). ∎

The above result also encompasses the case when the physical subsystem is
subject to a strict safety constraint, i.e., $x_{t}\in\mathcal{C}$ for all
$t\geq 0$. This case can be captured by letting $B=0$ in Definition 1. The
sufficient conditions for a control policy with strict safety guarantee are
given as follows.

######  Corollary 1.

Consider hybrid system $H$ and a safety set $\mathcal{C}$. Let
$h_{c}(x)=h(x)-c$ and $\mathcal{C}_{1}=\\{x:h_{c}(x)\geq 0\\}$. Consider an
arbitrary attack cycle denoted as $[t_{1},t_{2}]$ and suppose
$x_{0}\in\mathcal{C}_{1}$. If there exist constants $c\geq 0$, $\tau>0$, and a
control policy $\mu:\mathcal{X}\rightarrow\mathcal{U}$ such that

$\displaystyle\frac{\partial h}{\partial
x}(x)(f(x)+g(x)u)\geq-\frac{c}{N\delta},~{}\forall(x,u)\in\mathcal{C}\times\mathcal{U}$
(7a) $\displaystyle\frac{\partial h_{c}}{\partial
x}(x)(f(x)+g(x)\mu(x))\geq\frac{c}{\tau},~{}\forall
x\in\mathcal{C}\setminus\mathcal{C}_{1}$ (7b) $\displaystyle\frac{\partial
h_{c}}{\partial x}(x)(f(x)+g(x)\mu(x))\geq-\alpha(h_{c}(x)),~{}\forall
x\in\mathcal{C}_{1}$ (7c)

$x_{t}\in\mathcal{C}$ for all $t\geq 0$ provided that
$A=t_{2}-t_{1}\geq\tau+N\delta$.

###### Proof.

The corollary can be proved as a special case of Theorem 1 with $d=B=0$, which
yields that $\mathcal{D}=\mathcal{C}$ and thereby
$x_{t}\in\mathcal{C},~{}\forall t\in[t_{1},t_{2}]$. Note that $[t_{1},t_{2}]$
is an arbitrary attack cycle and $x_{0}\in\mathcal{C}_{1}$, rendering
$x_{t}\in\mathcal{C},~{}\forall t\geq 0$. ∎

The above analysis can also be used for the safety controller design of the
simplex architecture. The safety controller, which is invoked when the system
approaches the boundary of $\mathcal{C}$, can be obtained using Theorem 1 by
letting $c=d=B=0$ to guarantee strict safety with respect to $\mathcal{C}$.
Control policy $\mu$ only needs to satisfy Eqn. (4c) in this case.

Now we focus on the computation of control policy $\mu$ as well as parameters
$c,d\geq 0$ and $\tau>0$ so that safety is satisfied according to Definition
1. Our idea is to translate the conditions in Theorem 1 to a set of sum-of-
squares (SOS) constraints. We first make the following assumption.

######  Assumption 1.

We assume that functions $f(x)$, $g(x)$, and $h(x)$ are polynomial in $x$.
Additionally, we assume that function $L_{1}$ is polynomial in $-h(x)$.

When Assumption 1 holds, $L_{1}(-h(x))$ can be written as
$L_{1}(-h(x))=\sum_{i=0}^{k}(-h(x))^{i}a_{i}$, where $a_{i}$ is the
coefficient of $(-h(x))^{i}$ for each $i=0,\ldots,k$. Next we formulate Eqn.
(4) as a set of SOS constraints.

######  Proposition 1.

Suppose there exist parameters $c,d\geq 0$ and $\theta>0$ such that the
following expressions are SOS:

$\displaystyle\frac{\partial h_{d}}{\partial
x}(x)[f(x)+g(x)u]+\frac{c+d}{N\delta}-q(x,u)h_{d}(x)$ (8a)
$\displaystyle-\sum_{i=1}^{m}(w_{i}(x,u)([u]_{i}-[u]_{i,min})+v_{i}(x,u)([u]_{i,max}-[u]_{i})),$
$\displaystyle\frac{\partial h_{c}}{\partial
x}(x)[f(x)+g(x)\lambda(x)]-(c+d)\theta$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad-l(x)h_{d}(x)+p(x)h_{c}(x),$
(8b) $\displaystyle\frac{\partial h_{c}}{\partial
x}(x)[f(x)+g(x)\lambda(x)]+\alpha(h_{c}(x))-r(x)h_{c}(x),$ (8c)
$\displaystyle\lambda_{i}(x)-[u]_{i,min},\quad[u]_{i,max}-\lambda_{i}(x),~{}\forall
i=1,\ldots,m,$ (8d)

and the following inequality holds:

$B(c+d)-(N\delta+\frac{1}{\theta})\sum_{i=1}^{k}\frac{a_{i}d^{i+1}}{i+1}\geq
0$ (9)

where $l(x),p(x),q(x,u),r(x)$ are SOS, $\lambda_{i}(x)$ is a polynomial in $x$
for each $i=1,\ldots,m$, and $w_{i}(x,u)$ and $v_{i}(x,u)$ are SOS for each
$i=1,\ldots,m$. Then
$\mu(x)=\lambda(x)=[\lambda_{1}(x),\ldots,\lambda_{m}(x)]^{\top},c,d$, and
$\tau=\frac{1}{\theta}$ satisfy the conditions in Eqn. (4).

###### Proof.

Consider $x\in\mathcal{D}$ and $[u]_{i,min}\leq[u]_{i}\leq[u]_{i,max}$ for all
$i=1,\ldots,m$. Then we have that $h_{d}(x)\geq 0$, $[u]_{i}-[u]_{i,min}\geq
0$, and $[u]_{i,max}-[u]_{i}\geq 0$. Since expression (8a), $q(x,u)$,
$w_{i}(x,u)$, and $v_{i}(x,u)$ are SOS for all $i=1,\ldots,m$, therefore for
all $(x,u)\in\mathcal{D}\times\mathcal{U}$ we can write

$\frac{\partial h_{d}}{\partial x}(x)[f(x)+g(x)u]+\frac{c+d}{N\delta}\geq
q(x,u)h_{d}(x)\\\
+\sum_{i=1}^{m}(w_{i}(x,u)([u]_{i}-[u]_{i,min})+v_{i}(x,u)([u]_{i,max}-[u]_{i}))\geq
0.$

Thus condition (4a) holds.

Expressions (8b) to (8d) can be proved similarly. Eqn. (9) follows by
computing the integrals in Eqn. (4d). Details are omitted due to space
constraint. ∎

Algorithm 1 Heuristic algorithm for computing $c$, $d$, $\tau$ and control
policy $\mu(x)$

1:Input: $f(x)$, $g(x)$, $B$, $\tau_{max}$, $c_{max}$, $\epsilon_{1}>0$,
$\epsilon_{2}>0$

2:Output: $c$, $d$, $\tau$, $\lambda(x)$

3:Initialization: $c=0$.

4:while $c\leq c_{max}$ do

5: $d=0$

6: while $d\leq d_{max}$ do

7: Maximize ${\theta}$ subject to (8) with $c$ and $d$ fixed.

8: if Eqn. (8) is feasible, (9) is satisfied and $\frac{1}{\theta}<\tau_{max}$
then

9: return $d$, $c$, $\tau=\frac{1}{\theta}$, and $\lambda(x)$

10: else

11: $d=d+\epsilon_{1}$

12: end if

13: end while

14: $c=c+\epsilon_{2}$

15:end while

Simultaneously searching for $\lambda(x),c,d$ and $\theta$ that satisfy
Proposition 1 leads to bilinearity in (8). To this end, we propose an
algorithm to compute $\lambda(x),c,d$ and $\theta$ that satisfy Proposition 1,
as shown in Algorithm 1. Algorithm 1 first initializes parameters $c=d=0$. At
each iteration, the algorithm maximizes $\theta$ using the given $c$ and $d$.
If some $\theta^{*}$ can be found in line 7 which satisfies the conditions in
line 8 ($\tau_{max}=\infty$ if not specified), then the algorithm returns $c$,
$d$, and set $\tau=\frac{1}{\theta^{*}}$ and $\mu(x)=\lambda(x)$. Otherwise,
the algorithm increases the values of parameters $c$ and/or $d$ and repeat the
search process for parameter $\theta$. Algorithm 1 terminates at $c=c_{max}$
and $d=d_{max}$ if no feasible solution to (8) and (9) is found, where
$c_{max}=\sup_{x\in\mathcal{C}}h(x)$ and $d_{max}$ is the maximum value of $d$
that satisfies $\sum_{i=0}^{k}N\delta\frac{a_{i}d^{i+1}}{i+1}\leq(c+d)B$.

Now we briefly characterize the convergence and complexity of Algorithm 1. Our
intuition is that if there exists a feasible solution satisfying Eqn. (8) and
(9) strictly, then this solution must lie within the interior of the feasible
solution set. Thus by choosing $\epsilon_{1}$ and $\epsilon_{2}$ appropriately
small, the convergence of Algorithm 1 can be guaranteed. We formalize this
convergence result in the following proposition.

######  Proposition 2.

Suppose the set $\mathcal{C}$ is compact and $L_{1}$ is polynomial in $-h(x)$
with non-zero degree. Further assume that the feasible solution $(c,d)$ for
(8) and (9) satisfy inequality constraints in (8) and (9) strictly with $0\leq
c\leq c_{max}$, $0\leq d\leq d_{max}$ and $0<\tau\leq\tau_{max}$. Then
Algorithm 1 finds a feasible solution with $0<\tau\leq\tau_{max}$ in finite
number of iterations if $\epsilon_{1}$ and $\epsilon_{2}$ are chosen
appropriately small.

###### Proof.

Since $\mathcal{C}$ is compact, we have
$c_{max}=\sup_{x\in\mathcal{C}}h(x)<\infty$. Also, since the order of $d$ in
$\sum_{i=0}^{k}N\delta\frac{a_{i}d^{i+1}}{i+1}$ is greater than that of
$(c+d)B$ and $L_{1}$ is non-negative, therefore $d_{max}<\infty$. Using
$c_{max},d_{max}<\infty$, we have that Algorithm 1 terminates in finite number
of iterations.

Let $(c,d)$ satisfy (8) and (9) strictly with $0\leq c\leq c_{max}$, $0\leq
d\leq d_{max}$, and $0<\tau\leq\tau_{max}$. Therefore there exists an interval
$\mathcal{I}\subseteq\mathbb{R}^{2}$ for $(c,d)$ with non-zero measure for
which (8) and (9) are feasible and $c\in[0,c_{max}]$, $d\in[0,d_{max}]$,
$\tau\in(0,\tau_{max}]$. Denote the length of $\mathcal{I}$ in $c$ and $d$ as
$\tilde{c}>0$ and $\tilde{d}>0$, respectively. Let
$\epsilon_{1}\in(0,\tilde{c})$ and $\epsilon_{2}\in(0,\tilde{d})$. Then
Algorithm 1 will always terminate with a feasible solution. Otherwise interval
$\mathcal{I}$ contains some infeasible solutions to Eqn. (8) and (9), which
contradicts its definition. ∎

By the proof of Proposition 2, the computational complexity of Algorithm 1 is
$\lfloor\frac{c_{max}}{\epsilon_{1}}\rfloor\lfloor\frac{d_{max}}{\epsilon_{2}}\rfloor
M$, where $M$ is the computational complexity of line 7 in Algorithm 1.

##  VI Case Studies

This section presents a case study on lateral control of a Boeing 747. The
lateral dynamics for a Boeing 747 (at Mach 0.8 and 40000ft) [28] are given as
$\dot{x}=f(x)+g(x)u$, where

$f(x)=\begin{bmatrix}-0.0558&-0.9968&0.0802&0.0415\\\
0.598&-0.115&-0.0318&0\\\ -3.05&0.388&-0.465&0\\\ 0&0.0805&1&0\end{bmatrix},$

$g(x)=[0.00729,-0.475,0.153,0]^{\top}$, $x\in\mathbb{R}^{4}$ with
$x_{1},x_{2},x_{3}$, and $x_{4}$ representing the side-slip angle, yaw rate,
roll rate, and roll angle, respectively. The cyber subsystem updates the
control signal with frequency $20Hz$. The aircraft aims at maintaining the yaw
rate $x_{2}$ within $\mathcal{C}=\\{x:h(x)\geq 0\\}$ for passengers’ comfort
and minimizing potential damage to baggage, where $h(x)=0.025^{2}-x_{2}^{2}$.
Set $\mathcal{C}$ is represented as the white region in Fig. 2. We set
$x_{0}=[0.01,0.025,0,0]^{\top}$, $B=0.02$, and study the following three
scenarios.

_Scenario I: The aircraft has redundant controllers._ We choose parameter
$N=N_{1}+N_{2}$ epochs during which the aircraft is either corrupted or under
restoration, where $N_{1}=N_{2}=2$ [9]. We consider the system is equipped
with a buffer of length $3$. Using Algorithm 1, we obtain that $c=0$, $d=0.4$,
and $\tau=0.16s$, indicating that we need $4$ epochs for the system to be at
location $l_{1}=normal$ after restoration. We let the length of each attack
cycle be $8$ epochs. The control policy given by Algorithm 1 is $u=Kx$ with
$K=[-9.231\times 10^{-3},0.503,-1.805\times 10^{-3},2.373\times 10^{-5}]$. We
depict the trajectory of the yaw rate over $150$ epochs using the black solid
line in Fig. 2. The non-smoothness in the trajectory is due to switching
between location $normal$ (the controller is available) and locations
$corrupted$ and $restoration$ (the controller is unavailable). The adversary
enforces the yaw rate to exceed $0.025$ from the second to fifth epoch with
cost $0.0038<B$.

Figure 2: The yaw rate of a Boeing 747 adopting different architectures over
$150$ epochs. Set $\mathcal{C}=\\{x_{2}:0.025^{2}-x^{2}\geq 0\\}$ represented
by the white region. The black solid line depicts the yaw rate of an aircraft
with redundancy. The blue dash-dotted line is the yaw rate of an aircraft
without redundancy. The green solid line describes the yaw rate of an aircraft
using the simplex architecture without attack. The safety controller is
invoked when $0.025^{2}-x_{2}^{2}<0$.

_Scenario II: The aircraft has no redundant controller._ In this case, the
aircraft restarts the controller to recover the system. We let $N=N_{3}+N_{4}$
with $N_{3}=2$ and $N_{4}=4$. In this case, we let the adversary attack every
$10$ epochs. Algorithm 1 gives that $c=0$, $d=0.4$, $\tau=0.18s$ (i.e., 4
epochs), and $u=Kx$ with $K=[0.03017,0.05395,-4.753\times 10^{-3},7.513\times
10^{-5}]^{\top}$. The evolution of yaw rate over $150$ epochs is plotted using
blue dash-dotted line in Fig. 2. The adversary enforces the yaw rate to exceed
$0.025$ from the second to eleventh epoch with cost $0.0104<B$.

_Scenario III: The aircraft adopts simplex architecture._ We assume that the
main controller is in a faulty condition and produces random control input
$u_{t}\in\mathcal{U}$ for each epoch. The safety controller is invoked once
$h(x)<0$. Once the yaw rate exceeds $0.025$ (from the first to seventh epoch
in Fig. 2), the safety controller drives the yaw rate to $\mathcal{C}$. Note
that the simplex architecture assumes that there exists no adversary, which is
different with Scenario I and II.

Therefore, the control policy computed using our proposed algorithm ensures
safety of the system with respect to specified budget for any of the CRAs or
the simplex architecture.

##  VII conclusion

In this paper, we studied the problem of developing a common framework that
allows safety analysis and control synthesis of CPS adopting the simplex
architecture or the set of cyber resilient architectures including BFT++. We
presented the models for cyber and physical subsystems, and formulated the
safety property using a budget constraint. Our formulation captures strict
safety constraint as a special case. We constructed a hybrid system that
models CPS implementing any of these architectures. We derived a set of
sufficient conditions for the control policy to satisfy the budget constraint.
We translated the conditions into a set of sum-of-squares constraints, and
proposed an algorithm to compute the control policy. We analyzed the
convergence and complexity of the algorithm. A case study on the lateral
control of a Boeing 747 was presented to demonstrate viability of our proposed
framework.

## References

  * [1] 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/
  * [2] M. R. Lee, J. M. Assante, and T. Conway, “Analysis of the cyber attack on the Ukrainian power grid,” 2016. [Online]. Available: https://media.kasperskycontenthub.com/wp-content/uploads/sites/43/2016/05/20081514/E-ISAC_SANS_Ukraine_DUC_5.pdf
  * [3] Y. Zhang and J. Jiang, “Bibliographical review on reconfigurable fault-tolerant control systems,” _Annual Reviews in Control_ , vol. 32, no. 2, pp. 229–252, 2008.
  * [4] F. Sharifi, M. Mirzaei, B. W. Gordon, and Y. Zhang, “Fault tolerant control of a quadrotor UAV using sliding mode control,” in _Conference on Control and Fault-Tolerant Systems (SysTol)_. IEEE, 2010, pp. 239–244.
  * [5] D. Xu, F. Zhu, Z. Zhou, and X. Yan, “Distributed fault detection and estimation in cyber–physical systems subject to actuator faults,” _ISA Transactions_ , vol. 104, pp. 162–174, 2020.
  * [6] L. Sha, “Using simplicity to control complexity,” _IEEE Software_ , vol. 18, no. 4, pp. 20–28, 2001.
  * [7] S. Bak, D. K. Chivukula, O. Adekunle, M. Sun, M. Caccamo, and L. Sha, “The system-level simplex architecture for improved real-time embedded system safety,” in _15th IEEE Real-Time and Embedded Technology and Applications Symposium_. IEEE, 2009, pp. 99–107.
  * [8] S. Mohan, S. Bak, E. Betti, H. Yun, L. Sha, and M. Caccamo, “S3A: Secure system simplex architecture for enhanced security and robustness of cyber-physical systems,” in _the 2nd ACM International Conference on High Confidence Networked Systems_ , 2013, pp. 65–74.
  * [9] 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.
  * [10] M. A. Arroyo, M. T. I. Ziad, H. Kobayashi, J. Yang, and S. Sethumadhavan, “YOLO: frequently resetting cyber-physical systems for security,” in _Autonomous Systems: Sensors, Processing, and Security for Vehicles and Infrastructure 2019_ , vol. 11009. International Society for Optics and Photonics, 2019, p. 110090P.
  * [11] M. Arroyo, H. Kobayashi, S. Sethumadhavan, and J. Yang, “Fired: frequent inertial resets with diversification for emerging commodity cyber-physical systems,” _arXiv preprint arXiv:1702.06595_ , 2017.
  * [12] 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.
  * [13] M. Pajic, Z. Jiang, I. Lee, O. Sokolsky, and R. Mangharam, “Safety-critical medical device development using the UPP2SF model translation tool,” _ACM Transactions on Embedded Computing Systems (TECS)_ , vol. 13, no. 4s, pp. 1–26, 2014.
  * [14] A. D. Ames, S. Coogan, M. Egerstedt, G. Notomista, K. Sreenath, and P. Tabuada, “Control barrier functions: Theory and applications,” in _18th European Control Conference (ECC)_. IEEE, 2019, pp. 3420–3431.
  * [15] M. H. Cohen and C. Belta, “Approximate optimal control for safety-critical systems with control barrier functions,” in _59th Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 2062–2067.
  * [16] Z. Qin, K. Zhang, Y. Chen, J. Chen, and C. Fan, “Learning safe multi-agent control with decentralized neural barrier certificates,” _arXiv preprint arXiv:2101.05436_ , 2021.
  * [17] S. L. Herbert, M. Chen, S. Han, S. Bansal, J. F. Fisac, and C. J. Tomlin, “Fastrack: A modular framework for fast and guaranteed safe motion planning,” in _56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 1517–1522.
  * [18] 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)_. IEEE, 2014, pp. 163–174.
  * [19] 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.
  * [20] A. A. Cárdenas, S. Amin, and S. Sastry, “Research challenges for the security of control systems.” in _Proceedings of the 3rd Conference on Hot Topics in Security_ , vol. 5. USENIX Association, 2008, p. 15.
  * [21] L. Niu, D. Sahabandu, A. Clark, and P. Radha, “Verifying safety for resilient cyber-physical systems via reactive software restart,” in _To Appear in ACM/IEEE 13th International Conference on Cyber-Physical Systems (ICCPS)_. ACM/IEEE, 2022.
  * [22] P. Larsen, A. Homescu, S. Brunthaler, and M. Franz, “SoK: Automated software diversity,” in _2014 IEEE Symposium on Security and Privacy_. IEEE, 2014, pp. 276–291.
  * [23] G. S. Kc, A. D. Keromytis, and V. Prevelakis, “Countering code-injection attacks with instruction-set randomization,” in _Proceedings of the 10th ACM conference on Computer and communications security_ , 2003, pp. 272–280.
  * [24] F. Abdi, R. Tabish, M. Rungger, M. Zamani, and M. Caccamo, “Application and system-level software fault tolerance through full system restarts,” in _ACM/IEEE 8th International Conference on Cyber-Physical Systems (ICCPS)_. IEEE, 2017, pp. 197–206.
  * [25] 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 _ACM/IEEE 9th International Conference on Cyber-Physical Systems (ICCPS)_. IEEE, 2018, pp. 10–21.
  * [26] R. Romagnoli, P. Griffioen, B. H. Krogh, and B. Sinopoli, “Software rejuvenation under persistent attacks in constrained environments,” _IFAC-PapersOnLine_ , vol. 53, no. 2, pp. 4088–4094, 2020.
  * [27] T. Arauz, J. Maestre, R. Romagnoli, B. Sinopoli, and E. Camacho, “A linear programming approach to computing safe sets for software rejuvenation,” _IEEE Control Systems Letters_ , vol. 6, pp. 1214–1219, 2021.
  * [28] G. F. Franklin, J. D. Powell, A. Emami-Naeini, and J. D. Powell, _Feedback Control of Dynamic Systems_. Prentice hall Upper Saddle River, 2002, vol. 4.


# Two-dimensional Penta-NiPS Sheets: Two Stable Polymorphs

Klichchupong Dabsamut  Department of Physics, Faculty of Science, Kasetsart
University, Chatuchak, Bangkok 10900, Thailand  Thanasee Thanasarnsurapong
Department of Physics, Faculty of Science, Kasetsart University, Chatuchak,
Bangkok 10900, Thailand  Intuon Chatratin  Department of Materials Science and
Engineering, University of Delaware, Newark, Delaware 19716, USA  Tosapol
Maluangnont  Electroceramics Research Laboratory, College of Materials
Innovation and Technology, King Mongkut’s Institute of Technology Ladkrabang,
Ladkrabang, Bangkok 10520, Thailand  Sirichok Jungthawan  School of Physics,
Institute of Science, and Center of Excellence in Advanced Functional
Materials, Suranaree University of Technology, Muang, Nakhon Ratchasima,
30000, Thailand  Adisak Boonchun<EMAIL_ADDRESS>Department of Physics,
Faculty of Science, Kasetsart University, Chatuchak, Bangkok 10900, Thailand

###### Abstract

The discovery of new and stable two-dimensional (2D) materials with exotic
properties is essential for technological advancement. Inspired by the
recently reported penta-PdPSe, we proposed penta-NiPS as a new member of the
penta-2D materials based on first-principles calculations. The penta-NiPS
monolayer is stable in two polymorphs including the $\alpha$ phase with
identical structure as penta-PdPSe, and the newly proposed $\beta$ phase with
rotated sublayers. Comprehensive analyses indicated that both phases are
thermodynamically, dynamically, mechanically, and thermally stable. The penta-
NiPS is a soft material with 2D Young’s modulus of $E_{a}$ = 208 Nm-1 and
$E_{b}$ = 178 Nm-1 for the $\alpha$ phase and $E_{a}$ = 184 Nm-1 and $E_{b}$ =
140 Nm-1 for the $\beta$ phase. Interestingly, the $\alpha$-penta-NiPS showed
nearly zero Poisson’s ratios along the in-plane direction, where its
dimensions would be maintained when being extended. For electronic
application, we demonstrated that penta-NiPS is the wide band gap
semiconductors with an indirect band gap of 2.35 eV for $\alpha$ phase, and
2.20 eV for $\beta$ phase.

###### keywords:

American Chemical Society, LaTeX

IR,NMR,UV

##  1 Introduction

Since the discovery of hexacyclic graphene, two-dimensional (2D) materials
have been rapidly advanced in terms of production and materials design. An
increasing demand on their discovery has been driving the theoretical
explorations of stable 2D materials beyond the conventional honeycomb-like
structure, including the understanding of composition-structure-properties
relationships. Recently, Zhang et al.1 theoretically proposed a new pentagonal
structure penta-graphene with carbon atoms uniquely resembling Cairo
pentagonal tiling. Although not yet experimentally synthesized, the properties
of the penta-graphene have been widely explored through first-principles
calculations. The wide band gap of 3.25 eV and high strength could lead to
potential applications in nanoelectronics and nanomechanical devices. Other
prominent examples are SiC2, CN2, BN2, B2C, PdS2, AlN2, PtN2, PdN2, and NiS2
2, 3, 4, 5, 6, 7, 8, 9. In all cases, it is essential that one thoroughly
explores various aspects of materials’ stability including dynamic,
mechanical, and thermal ones.

The investigation of ternary compounds provides a higher degree of freedom in
tuning the properties of penta-2D materials 10. The prominent example is
penta-BCN 10 which is thermodynamically, mechanically, and thermally stable
according to density functional theory (DFT), in addition to possessing a high
piezoelectric response. Later on, the penta-PdPSe with intrinsic in-plane
anisotropy has been experimentally synthesized.11 This material is made up
from two PdPSe sublayers which are uncommon in nature, leading to unusual
electronic and optical properties.12, 13 Nevertheless, penta-PdPSe is the
closely related member of the binary penta-PdSe2 which can be readily
synthesized via mechanical exfoliation from bulk crystals.14 It has been
experimentally demonstrated that the penta-PdSe2 based field-effect
transistors (FETs) exhibited band gap tunability and high electron mobility.14
Accordingly, chemical intuition with appropriate atom substitution is a simple
yet potential method for theoretically discovering novel materials.

Herein, we used first-principles calculations to search for a novel member in
the ternary pentagonal sheet family. Based on the reported atomic
configuration of penta-PdPSe, Pd was replaced with Ni, and Se was replaced
with S in this work. We found that the proposed composition penta-NiPS exists
in two polymorphs which are thermodynamically, dynamically, mechanically, and
thermally stable. The $\alpha$ phase has the structure similar to penta-PdPSe.
Meanwhile, another stable structure, the $\beta$ phase, can be constructed by
rotating the sublayers from the $\alpha$ phase. Notably, the $\alpha$ phase
exhibited zero Poisson’s ratio (ZPR) in the principle $x$ and $y$ directions
(i.e., along the sheet), a unique ability to maintain its shape when being
compressed or extended. This is rather unusual because most solids exhibit a
positive Poisson’s ratio (PPR) in the range of 0.2-0.3,15 which describes the
level of transversally deformation under an axial strain. The electronic and
mechanical properties of two phase have also been compared.

##  2 Computational Method

Our calculations are based on density functional theory (DFT) and the plane-
wave basis projector augmented wave (PAW)16 method, as defined in the Vienna
Ab-initio Simulation Package (VASP)17, 18. The kinetic energy cutoff of 500 eV
is used for the plane-wave expansion. Perdew-Burke-Ernzerhof (PBE)19
functional is imposed for geometry optimization. In order to accurately
describe layered material, the long-range van der Waals (vdW) interactions
were taken into account using the DFT-D2 functional proposed by Grimme.20 The
structures are optimized until total Hellmann-Feynman forces are decreased to
0.02 eV/Åand the energy convergence is set to $10^{-6}$ eV. The centered
k-point grid within $7\times 7\times 1$ is employed for first Brillouin zone
integration within the Monkhorst-Pack21 scheme. A vacuum layer of 16 Å is
included along the direction perpendicular to a crystal layer to avoid the
interaction between the layer and its neighboring image. A PHONOPY package22
over $3\times 3\times 1$ supercell with 500 eV of energy cutoff and a single
k-point at $\Gamma$ is used to calculate phonon spectra to verify the dynamic
stability. Ab initio molecular dynamics (AIMD) simulations is employed to
confirm thermodynamic stability. To perform electronic band structure, since
DFT within the standard approximation severely underestimates band gaps, we
use the Heyd-Scuseria-Ernzerhof (HSE06)23, 24 hybrid functional to correct not
only the band gap but also the band width.

##  3 Results and discussion

Figure 1:  3x3 supercell of pentagonal (a) $\alpha$-NiPS and (b) $\beta$-NiPS
geometrical structures.

We first constructed the penta-NiPS based on the reported atomic configuration
of PdPSe, by replacing the Pd atoms with Ni atoms and Se atom with S atoms.
This is called the $\alpha$-phase which composes of two sublayers. Within one
sublayer, every Ni atom is linked by two P atoms and two S atoms, showing the
pentagon structure, as shown in Fig. 1 (a). The pentagonal primitive unit
cell, shown as the square in Fig. 1 (a), consists of 12 atoms including 4 Ni
atoms, 4 P atoms, and 4 S atoms. The lattice constants of the $\alpha$ phase
penta-NiPS are a = 5.43 Å and b = 5.38 Å. The whole thickness of the layer,
T1, is 3.98 Å, while the distance between two sublayers, T2, is 2.20 Å .

In addition, we discovered another stable phase, namely the $\beta$ phase,
where the two sublayers directly stack on top of each other as shown in Fig.
1(b). The coordinates of the optimized geometry are provided in the Supporting
Information (SI). The 12-atom primitive unit cell is shown in the square in
Fig. 1(b). The lattice constants of the $\beta$ phase are a = 5.45 Å and b =
5.41 Å, which is slightly larger than $\alpha$ phase. Meanwhile, T1 is 3.86 Å,
and T2 is 2.29 Å . The larger T2 of the $\beta$ phase is understandable
because the direct stacking generates a higher coulomb repulsion of the same
atoms between two sublayers, as is commonly observed in other 2D materials.
The fundamental properties of the two phases including bond lengths, bond
angles and band gap are listed in Table. 1.

Table 1: The calculated lattice constants a and b in Å, bond lengths $d_{1}$-$d_{6}$ in Å, bond angles $\theta_{1}$-$\theta_{11}$ in degree, thicknesses T1 and T2 in Å, the cohesive energy per atom in eV/atom, and band gap $E_{g}^{PBE}$ and $E_{g}^{HSE}$ in eV | a | b | d1 | d2 | d3 | d4 | d5 | d6 |
---|---|---|---|---|---|---|---|---|---
$\alpha$-NiPS | 5.43 | 5.38 | 2.18 | 2.18 | 2.18 | 2.17 | 2.14 | 2.23 |
$\beta$-NiPS | 5.45 | 5.41 | 2.17 | 2.17 | 2.17 | 2.18 | 2.12 | 2.29 |
| $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | $\theta_{4}$ | $\theta_{5}$ | $\theta_{6}$ | $\theta_{7}$ | $\theta_{8}$ | $\theta_{9}$
$\alpha$-NiPS | 88.12 | 105.45 | 109.61 | 92.22 | 120.42 | 91.10 | 125.10 | 88.15 | 108.21
$\beta$-NiPS | 88.71 | 107.72 | 110.31 | 92.18 | 122.06 | 90.63 | 126.21 | 88.21 | 109.66
| $\theta_{10}$ | $\theta_{11}$ | T1 | T2 | $E_{coh}$ | $E_{g}^{PBE}$ | $E_{g}^{HSE}$ |  |
$\alpha$-NiPS | 108.45 | 83.48 | 3.98 | 2.20 | -5.44 | 0.93 | 2.35 |  |
$\beta$-NiPS | 110.04 | 90.00 | 3.86 | 2.29 | -5.42 | 0.71 | 2.20 |  |

The electronic properties of $\alpha$\- and $\beta$-NiPS phases are first
elucidated by the calculated band structure as shown in Fig.2. The dash and
solid lines are energy bands calculated from PBE and HSE06 functional,
respectively. By using the PBE functional,19 the electronic band structures
show indirect band gaps of 0.88 and 0.68 eV for $\alpha$ and $\beta$ phase,
respectively. As it is well-known that the standard DFT underestimates band
gap, the HSE06 functional23, 24 was also employed. It is found that a larger
band gap of 2.35 eV was obtained for the $\alpha$ phase, vs 2.20 eV for the
$\beta$ phase, with the indirect band gap character being consistent. As seen
in Fig.2 within HSE06, the valence band maximum (VBM) of $\alpha$ phase is
located along the $M-Y$ path and the conduction band minimum (CBM) at the $M$
point. In the $\beta$ phase, the VBM is located along the $\Gamma-Y$ path and
the CBM is at the $\Gamma$ point.

Figure 2:  Calculated band structure and partial density of state (PDOS) for
(a) $\alpha$-NiPS and (b) $\beta$-NiPS, respectively.

In order to analyze the electronic properties of penta-NiPS, the partial
density of states (PDOS) obtained from the HSE0623, 24 functional were plotted
along with the band structure in Fig.2. For $\alpha$-NiPS, the VBM is
dominated by the P-$p$ and the S-$p$ state, while the CBM has a main
contribution from the Ni-$d$ states. For $\beta$-NiPS, both of VBM and CBM are
mainly contributed by the Ni-$d$ states. The average $d$ band center of Ni in
the $\alpha$ phase is $\sim 1.85\,{\rm eV}$ with respect to VBM, while the
average $d$ band center in the $\beta$ phase is $\sim 1.65\,{\rm eV}$ with
respect to VBM. It is known that the deeper energy of the $d$ states weakens
the $p-d$ coupling, leading to a further increase in band gap 25 as observed
in the present work.

To evaluate the crystal binding stability, the cohesive energy per atom was
calculated, given by the following equation:

$E_{coh}=\frac{E_{tot}(\text{NiPS})-\sum n_{i}E_{i}}{\sum n_{i}}$ (1)

where $E_{tot}(\text{NiPS})$ is the total energy of the penta-NiPS, $n_{i}$ is
the number of Ni, P and S atoms in the unit cell and $E_{i}$ represents the
energies of an isolated single Ni, P and S atom. The calculated cohesive
energy of the penta-NiPS are -5.44 and -5.42 eV/atom for the $\alpha$ and
$\beta$ phase, respectively. The negative values indicate that the penta-NiPS
are more energetically favorable than the isolated atoms. In addition, the
$\alpha$ phase is just slighlty more energetically favorable than the $\beta$
phase. Generally, different polymorphs (e.g., anatase vs rutile for $\rm
TiO_{2}$) could be synthesized by a fine tuning of synthetic conditions such
as temperature, pressure, surrounding gas, among others.26 The selective
synthesis of either $\alpha$ and $\beta$ penta-NiPS remains to be
accomplished.

Figure 3:  Calculated phonon band structure of (a) $\alpha$-NiPS and (b)
$\beta$-NiPS.

We verify the dynamic stability of $\alpha$ and $\beta$ penta-NiPS by
calculating their phonon dispersion. The dynamic stability manifests itself as
the positive frequencies in the phonon spectra throughout the entire Brillouin
zone. On the other hand, if the spectra display imaginary frequency modes
(negative values), the structures reveal non-restorative force against a
displacement of atoms, indicating instability. Fig 3(a) and Fig 3(b) show the
calculated phonon dispersion of penta-NiPS for $\alpha$ and $\beta$ phase,
respectively. There are no imaginary modes throughout the Brillouin zone,
confirming that both are dynamically stable. We note that the imperceptible
negative frequencies at $\Gamma$, (visible only when zoomed in) is a well-
known computational error which can be effectively eliminated using the LDA
functional, for example. These negative frequencies of phonon dispersion are
also observed in penta-PdSe2.14

Mechanical stability is an important property for strain engineering. Such
stability can be verified from the linear elastic constants, C11, C22, C12 and
C66, which are elements in the stiffness tensor. Here, the elastic constants
are obtained by finite difference method. Our calculated elastic constants of
both polymorphs are listed in table 2. Typically, the structures are
considered mechanically stable when the linear elastic constants satisfy the
conditions C11C22-C${}_{12}^{2}>0$ and C${}_{66}>0$. Notably, these conditions
are obeyed for both phases, confirming the mechanical stability of penta-NiPS.

Table 2: The elastic constants $C_{jk}$ (Nm-1), Poisson’s ratio on [100] ($\nu_{a}$) and [010] ($\nu_{b}$) directions and Young’s modulus (Nm-1) on [100] ($E_{a}$) and [010] ($E_{b}$) directions. |  $\alpha$-NiPS |  $\beta$-NiPS
---|---|---
$C_{11}$ | 184.27 | 211.55
$C_{12}$ | 0.52 | 28.52
$C_{22}$ | 140.50 | 182.60
$C_{66}$ | 65.03 | 66.48
$\nu_{a}$ | 0.003 | 0.135
$\nu_{b}$ | 0.004 | 0.156
$E_{a}$ | 184.27 | 207.71
$E_{b}$ | 140.50 | 178.14

Next, we evaluated the 2D Young’s modulus on the [100] and [010] directions
(i.e., along the plane), which is obtained by $E_{a}=$(C${}^{2}_{11}$ \-
C${}^{2}_{12}$)/C11 and $E_{b}=$ (C${}^{2}_{22}$ \- C${}^{2}_{12}$)/C22. In
the $\alpha$ phase, $E_{a}$ = 207.71 Nm-1 and $E_{b}$ = 178.14 Nm-1. These
values are close to those of penta-PdPSe with identical structure, $E_{a}$ =
218 Nm-1 and $E_{b}$ = 152 Nm-1, calculated from DFT calculation.12 Therefore,
the $\alpha$-NiPS are as strong as penta-PdPSe. In $\beta$ phase, the
calculated values are $E_{a}$ = 184.27 Nm-1 and $E_{b}$ = 140.50 Nm-1,
indicating that it is softer than the $\alpha$ phase. Table 2 lists the
calculated Poisson’s ratio. For $\beta$-NiPS, the calculated Poisson’s ratio
on the [100] direction, $\nu_{a}$ = C12/C11 is 0.135 which is slightly lower
than that on the the [010] direction, $\nu_{b}$ = C12/C22 of 0.156.
Interestingly, for $\alpha$-NiPS, the calculated Poisson’s ratios are nearly
zero for both [100] and [010] directions ($0^{\circ}$ and $90^{\circ}$),
$\nu_{a}$ = 0.003 and $\nu_{b}$ = 0.004. The nearly zero Poisson’s ratios was
previously observed in ternary penta-BCN bus it is located near [110]
direction ($\sim 45^{\circ}$) with the value of 0.004.10 Moreover, Hou et al.
27 reported that the ternary penta-BCP exhibited the largest positive
Poisson’s ratio of 1.30 along the [010] direction, which is now the highest
among 2D materials. The near zero Poisson’s ratios of the present $\beta$-NiPS
indicate that $\beta$-NiPS has an ability to maintain its dimensions under
compression/extension along the in-plane directions.

Figure 4:  Variation of (a) in-plane Young’s modulus, and (b) Poisson’s ratio
for $\alpha$-NiPS and $\beta$-NiPS. Note that solid and dash line represent
$\alpha$-NiPS and $\beta$-NiPS, respectively.

To further explore the dependence of mechanical properties on the crystal
orientation of penta-NiPS, we computed the Young’s modulus ($E$) and Poisson’s
ratio ($\nu$) along an arbitrary direction $\theta$ ($\theta$ is the angle
relative to the x direction) using the following formulas:

$\displaystyle
E(\theta)=\frac{C_{11}C_{12}-C_{12}^{2}}{C_{11}s^{4}+C_{22}c^{4}+\left(\frac{C_{11}C_{12}-C_{12}^{2}}{C_{66}}-2C_{12}\right)c^{2}s^{2}}$
(2)
$\displaystyle\nu(\theta)=\frac{\left(C_{11}+C_{12}-\frac{C_{11}C_{12}-C_{12}^{2}}{C_{66}}\right)c^{2}s^{2}-C_{12}(c^{4}+s^{4})}{C_{11}s^{4}\theta+C_{22}c^{4}+\left(\frac{C_{11}C_{12}-C_{12}^{2}}{C_{66}}-2C_{12}\right)c^{2}s^{2}}$
(3)

where $c=\cos\theta$ and $s=\sin\theta$. The calculated results are plotted in
Fig 4. One can see that the maxima of the in-plane Young’s modulus in both
phases are located at the [100] direction with the value of 184.27
($\alpha$-NiPS) and 207.71 ($\beta$-NiPS) ${\rm Nm^{-1}}$, respectively.
Meanwhile, the minima are along the [010] direction (90∘) with the value of
140.50 ($\alpha$-NiPS) and 178.14 ($\beta$-NiPS) ${\rm Nm^{-1}}$. In
$\alpha$-NiPS, the minimum Poisson’s ratios (which is very close to zero) are
located along the [100] ($0^{\circ}$), while the maximum is along [110]
direction (43∘) with values of 0.10. In $\beta$-NiPS, the minimum and maximum
Poisson’s Ratios are located with the same angle as $\alpha$ phase but at
higher positive values. Because the in-plane Young’s modulus and Poisson’s
ratios are dependent on crystal orientation, penta-NiPS exhibit anisotropic
mechanical properties which can be explained by its lattice constant where
$a\neq b$.

Figure 5:  (a) the variation in biaxial strain energy (Es) and (b) its
derivative for $\alpha$-NiPS and $\beta$-NiPS. Blue and red color represent
$\alpha$-NiPS and $\beta$-NiPS, respectively.

For practical application of semiconductors, it is beneficial that one can
modify the electronic properties by a simple method such as strain
engineering. This can be effectively applied to 2D materials via lattice
mismatch on the substrate, thermal expansion, or mechanical loading. Because
2D materials tend to fold up under lateral compression, we focused on the
tensile strain only, specifically the biaxial strain effect. The degree of
strain is defined as $\varepsilon=\Delta l/l_{0}\times 100$, where $\Delta l$
is the difference in a strained length, $l$, and an original, unstrained
length, $l_{0}$. Fig.5 is the plot showing the strain energies, $E_{S}$, and
their derivative ($dE_{S}(\varepsilon)/d\varepsilon$) with respect to the
applied strains from 0 - 20%. For simplicity, the energies of the strain-free
structures are set to zero so that the strain energies, $E_{S}$ is the energy
differences between the total energies of the structure with and without
strain. When the biaxial strain is continuously applied, a crystal expands
until reaching the mechanical failure at $\varepsilon=\varepsilon_{C}$ where
the maximum derivative of strain is observed. The mechanical failure of the
$\alpha$-NiPS and $\beta$-NiPS occurs at $\varepsilon_{C}^{\alpha}=$ 16.90%
and 12.66%, respectively. These values can be compared with 24.3% for penta-
BCN, and 17.5% for hydrogenated penta-BCN (H-BCN).28, 29 Normally, 2D
materials can sustain much larger strains than their bulk structure. Also, 2D
materials with buckling structure can withstand a higher strain than the non-
buckling structure. Here, the resilience of the $\alpha$ phase could be linked
to the bucking height $\rm T_{1}$ compared to that of the $\beta$ phase.

It is known that the band gaps of 2D materials can be potentially modulated by
strain. Thus, we also investigated how the band gaps of $\alpha$-NiPS and
$\beta$-NiPS change with the additional strain, starting from their
equilibrium (i.e., zero strain) up to just prior to the critical point of
strain. In Fig. S1(a) and S1(b) in the supplementary information (SI), the
band gaps of both polymorphs under the strain seem to linearly decrease from
2.35 to 0.94 eV ($\alpha$) and 2.20 to 1.00 eV ($\beta$), respectively. In
$\alpha$-NiPS, the valence band maximum (VBM) shifted from $Y-\Gamma$ path to
$Y-M$ path while that of the conduction band minimum (CBM) changed from $M$
point to $\Gamma$ point. On the other hand, the characteristics of VBM and CBM
of $\beta$-NiPS have not been altered with strain as depicted in Fig. S2(a)
and S2(b) in SI.

Figure 6:  Total energy fluctuation with time during the AIMD simulation at
800 K of (a) $\alpha$-NiPS and (b) $\beta$-NiPS.

The thermal stability of penta-NiPS was examined via ab initio molecular
dynamic simulations. Supercells with a duplication of a $4\times 4\times 1$ of
the unit cell is used. The molecular dynamics simulations were performed for
15 ps with a time step of 1 fs under the temperature from 300 to 800 K. Fig. 6
shows the fluctuation of the potential energy of penta-NiPS throughout the
simulation time at 800 K. For both phases, the respective structure reaches
the steady state and equilibrates around a constant energy after 1,000 fs.
These results clearly indicate the thermal stability of penta-NiPS at least up
to 800 K.

It is true that both polymorphs of NiPS has not been experimentally realized.
Still, the related composition MPS3 (M = Mn, Fe, Co, Ni and Cd) is a classic
example of layered compounds crystallizing in CdCl2-type structure.30 Although
the two structures ($\alpha$\- and $\beta$-penta-NiPS vs CdCl2) are clearly
different, NiPS might be experimentally synthesized. Consider, for example,
PdSe2 which can be successfully synthesized by a self-flux method.14 Moreover,
very recently, penta-PdPS (the same structure with $\alpha$-NiPS) can be
successfully synthesized by exfoliated from bulk crystals via chemical vapor
transport (CVT) method.31 It is premature to explicitly suggest the synthetic
procedure at present, but we hope it could be achievable.

###  3.1 Conclusions

By using DFT calculation, we theoretically discovered penta-NiPS as the new
member in the penta-2D family. The penta-NiPS is stable in two polymorphs. The
$\alpha$-NiPS has the same structure as the previously reported penta-PdPSe,
while the $\beta$-NiPS was obtained by rotating the two sublayers from the
$\alpha$-NiPS. Both phases are dynamically, mechanically, and thermally
stable, as comprehensively verified from phonon dispersion, elastic constants,
and molecular dynamic simulation, respectively. We found that the penta-NiPS
is a semiconductor with an indirect band gap of 2.35 and 2.20 eV, for $\alpha$
phase and $\beta$ phase. They are soft materials with 2D Young’s modulus of
$E_{a}$ = 207.71 Nm-1 and $E_{b}$ = 178.14 Nm-1 for the $\alpha$ phase and
$E_{a}$ = 184.27 Nm-1 and $E_{b}$ = 140.50 Nm-1 for the $\beta$ phase.
Interestingly, the $\alpha$ penta-NiPS exhibited a close to zero Poisson’s
ratio of 0.003-0.004. These values are significantly smaller than those of the
$\beta$ phase, 0.135-0.156. The nearly zero Poisson’s ratios allow
$\beta$-penta-NiPS to maintain its dimensions under extension with potential
applications in wearable electronics.

This project is funded by National Research Council of Thailand (NRCT) and
Kasetsart University : N42A650278. K.D. was supported by the National Research
Council of Thailand (NRCT), No. NRCT5-RGJ63002-028. T.T. was supported by
Graduate Program Scholarship from the Graduate School, Kasetsart University.
S.J. has received funding support from the NSRF via the Program Management
Unit for Human Resources & Institutional Development, Research and Innovation
[grant number B05F640051]. We wish to thank NSTDA Supercomputer Center
(ThaiSC) for providing computing resources for this work.

The coordination of the optimized geometries of $\alpha$ and $\beta$ penta-
NiPS and additional figures as mentioned in the text.

## References

  * Zhang et al. 2015 Zhang, S.; Zhou, J.; Wang, Q.; Chen, X.; Kawazoe, Y.; Jena, P. Penta-graphene: A new carbon allotrope. _Proc. Natl. Acad. Sci. U.S.A._ 2015, _112_ , 2372–2377
  * Berdiyorov and Madjet 2016 Berdiyorov, G. R.; Madjet, M. E.-A. First-principles study of electronic transport and optical properties of penta-graphene, penta-SiC 2 and penta-CN 2\. _RSC Adv._ 2016, _6_ , 50867–50873
  * Zhang et al. 2016 Zhang, S.; Zhou, J.; Wang, Q.; Jena, P. Beyond graphitic carbon nitride: nitrogen-rich penta-CN2 sheet. _J. Phys. Chem. C_ 2016, _120_ , 3993–3998
  * Li et al. 2016 Li, J.; Fan, X.; Wei, Y.; Chen, G. Penta-B x N y sheet: a density functional theory study of two-dimensional material. _Sci. Rep._ 2016, _6_ , 31840
  * Li et al. 2015 Li, F.; Tu, K.; Zhang, H.; Chen, Z. Flexible structural and electronic properties of a pentagonal B2C monolayer via external strain: a computational investigation. _Phys. Chem. Chem. Phys._ 2015, _17_ , 24151–24156
  * Cheng et al. 2021 Cheng, Z.; Zhang, X.; Zhang, H.; Gao, J.; Liu, H.; Yu, X.; Dai, X.; Liu, G.; Chen, G. Pentagonal B2C monolayer with extremely high theoretical capacity for Li-/Na-ion batteries. _Phys. Chem. Chem. Phys._ 2021, _23_ , 6278–6285
  * Lan et al. 2019 Lan, Y.-S.; Chen, X.-R.; Hu, C.-E.; Cheng, Y.; Chen, Q.-F. Penta-PdX 2 (X= S, Se, Te) monolayers: promising anisotropic thermoelectric materials. _J. Mater. Chem. A_ 2019, _7_ , 11134–11142
  * Li et al. 2016 Li, J.; Fan, X.; Wei, Y.; Liu, H.; Li, S.; Zhao, P.; Chen, G. Half-metallicity and ferromagnetism in penta-AlN2 nanostructure. _Sci. Rep._ 2016, _6_ , 1–10
  * Wang and Du 2020 Wang, C.-T.; Du, S. A unique pentagonal network structure of the NiS 2 monolayer with high stability and a tunable bandgap. _Phys. Chem. Chem. Phys._ 2020, _22_ , 7483–7488
  * Zhao et al. 2020 Zhao, K.; Guo, Y.; Shen, Y.; Wang, Q.; Kawazoe, Y.; Jena, P. Penta-BCN: A New Ternary Pentagonal Monolayer with Intrinsic Piezoelectricity. _J. Phys. Chem. Lett._ 2020, _11_ , 3501–3506
  * Li et al. 2021 Li, P.; Zhang, J.; Zhu, C.; Shen, W.; Hu, C.; Fu, W.; Yan, L.; Zhou, L.; Zheng, L.; Lei, H., et al. Penta-PdPSe: A New 2D Pentagonal Material with Highly In-Plane Optical, Electronic, and Optoelectronic Anisotropy. _Adv. Mater._ 2021, _33_ , 2102541
  * Mortazavi et al. 2022 Mortazavi, B.; Shahrokhi, M.; Zhuang, X.; Rabczuk, T.; Shapeev, A. V. Mechanical, thermal transport, electronic and photocatalytic properties of penta-PdPS,-PdPSe and-PdPTe monolayers explored by first-principles calculations. _J. Mater. Chem. C_ 2022, _10_ , 329–336
  * Bafekry et al. 2022 Bafekry, A.; Fadlallah, M.; Faraji, M.; Shafique, A.; Jappor, H.; Sarsari, I. A.; Ang, Y. S.; Ghergherehchi, M. Two-dimensional penta-like PdPSe with a puckered pentagonal structure: a first-principles study. _Phys. Chem. Chem. Phys._ 2022, _24_ , 9990–9997
  * Oyedele et al. 2017 Oyedele, A. D.; Yang, S.; Liang, L.; Puretzky, A. A.; Wang, K.; Zhang, J.; Yu, P.; Pudasaini, P. R.; Ghosh, A. W.; Liu, Z., et al. PdSe2: Pentagonal Two-Dimensional Layers with High Air Stability for Electronics. _J. Am. Chem. Soc._ 2017, _139_ , 14090–14097, PMID: 28873294
  * Greaves et al. 2011 Greaves, G. N.; Greer, A. L.; Lakes, R. S.; Rouxel, T. Poisson’s ratio and modern materials. _Nat. Mater._ 2011, _10_ , 823–837
  * Blöchl 1994 Blöchl, P. E. Projector augmented-wave method. _Phys. Rev. B_ 1994, _50_ , 17953
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. _Comput. Mater. Sci._ 1996, _6_ , 15–50
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. _Phys. Rev. B_ 1996, _54_ , 11169
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. _Phys. Rev. Lett._ 1996, _77_ , 3865
  * Grimme 2006 Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. _J. Comput. Chem._ 2006, _27_ , 1787–1799
  * Monkhorst and Pack 1976 Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. _Phys. Rev. B_ 1976, _13_ , 5188
  * Togo and Tanaka 2015 Togo, A.; Tanaka, I. First principles phonon calculations in materials science. _Scr. Mater._ 2015, _108_ , 1–5
  * Heyd et al. 2003 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. _J. Chem. Phys._ 2003, _118_ , 8207–8215
  * Heyd et al. 2006 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)]. _J. Chem. Phys._ 2006, _124_ , 219906
  * Janotti et al. 2006 Janotti, A.; Segev, D.; Van de Walle, C. G. Effects of cation $d$ states on the structural and electronic properties of III-nitride and II-oxide wide-band-gap semiconductors. _Phys. Rev. B_ 2006, _74_ , 045202
  * Wells 2012 Wells, A. F. _Structural inorganic chemistry_ ; Clarendon Press: Oxford, U.K., 2012
  * Hou et al. 2022 Hou, C.; Shen, Y.; Sun, W.; Chen, Y.; Ni, D.; Wang, Q. A penta-BCP sheet with strong piezoelectricity and a record high positive Poisson’s ratio. _J. Mater. Chem. C_ 2022, _10_ , 10302–10309
  * Dabsamut et al. 2021 Dabsamut, K.; Thanasarnsurapong, T.; Maluangnont, T.; T-Thienprasert, J.; Jungthawan, S.; Boonchun, A. Strain engineering and thermal conductivity of a penta-BCN monolayer: a computational study. _J. Phys. D: Appl. Phys._ 2021, _54_ , 355301
  * Thanasarnsurapong et al. 2021 Thanasarnsurapong, T.; Dabsamut, K.; Maluangnont, T.; T-Thienprasert, J.; Jungthawan, S.; Boonchun, A. Piezoelectric and electronic properties of hydrogenated penta-BCN: A computational study. _J. Appl. Phys._ 2021, _129_ , 095101
  * Ouvrard et al. 1985 Ouvrard, G.; Brec, R.; Rouxel, J. Structural determination of some MPS3 layered phases (M= Mn, Fe, Co, Ni and Cd). _Mater. Res. Bull._ 1985, _20_ , 1181–1189
  * Duan et al. 2022 Duan, R.; He, Y.; Zhu, C.; Wang, X.; Zhu, C.; Zhao, X.; Zhang, Z.; Zeng, Q.; Deng, Y.; Xu, M., et al. 2D Cairo Pentagonal PdPS: Air-Stable Anisotropic Ternary Semiconductor with High Optoelectronic Performance. _Adv. Funct. Mater._ 2022, _32_ , 2113255