# CDC 2018

Amir Khazraei  Amir Khazraei, Hamed Kebriaei∗, and Farzad Rajaei Salmasi This
research was in part supported by a grant from the Institute for Research in
Fundamental Sciences (IPM) (No. CS 1397-4-56).Authors are with the School of
Electrical and Computer Engineering, College of Engineering, University of
Tehran, Tehran, Iran. H. Kebriaei is also with School of Computer Science,
Institute for Research in Fundamental Sciences (IPM), P.O. Box 19395-5746,
Tehran, Iran.∗ Corresponding author: H. Kebriaei (e-mail: kebriaei@ut.ac.ir).

(November 2018)

# An Optimal Linear Dynamic Detection Method for Replay Attack in Cyber-
Physical Systems

Amir Khazraei  Amir Khazraei, Hamed Kebriaei∗, and Farzad Rajaei Salmasi This
research was in part supported by a grant from the Institute for Research in
Fundamental Sciences (IPM) (No. CS 1397-4-56).Authors are with the School of
Electrical and Computer Engineering, College of Engineering, University of
Tehran, Tehran, Iran. H. Kebriaei is also with School of Computer Science,
Institute for Research in Fundamental Sciences (IPM), P.O. Box 19395-5746,
Tehran, Iran.∗ Corresponding author: H. Kebriaei (e-mail: kebriaei@ut.ac.ir).

(November 2018)

###### Abstract

The problem of detecting replay attack to the linear stochastic system with
Kalman filer state estimator and LQG controller is addressed. To this end, a
dynamic attack detector method is proposed which is coupled with the dynamics
of the system. While preserving stability of the main system, conditions on
parameters of the attack detector dynamics are obtained such that the attack
can be revealed by destabilization of a residual trajectory which is the
difference between the estimated and measured output of the system. Using this
method, system operator can adjust the detection rate based on the proposed
scheme by changing the design parameters. Nevertheless, since the exogenous
attack detector signal affects the performance of the closed loop control
system, we propose an optimization problem to determine such a detector with
minimum loss effect. In the simulation results, the proposed dynamical attack
detector approach is compared with the well-known additive white noise
watermarking method and the results confirm the superiority of the new scheme.

##  I Introduction

Cyber-physical systems refer to a new generation of systems which integrates
computational, networking and physical processes. Security of Cyber-Physical
systems from control theory point of view has gained lots of attention during
the last few years [1]. Due to widespread use of computer network in control
of physical systems, the malicious entities have seized an opportunity to
inject attacks on cyber-physical systems [2-3].

The analysis over vulnerabilities of cyber-physical systems has attracted a
tremendous amount of attention in past few years. For instance, false data
injection attack against static state estimation in a power grid has been
studied in [4]. The authors showed that if the injected signal to the sensory
measurements of the system lays within the range space of the observation
matrix, then the attack remains stealthy. Conditions on the existence of
undetectable and unidentifiable attacks has been characterized in [5] using
graph theoretic method. It is assumed that the adversary is capable of
injecting the attack signal into the input and output of the system. In [6],
the cyber-physical attack is modeled based on adversary resources including
disruption, disclosure and information of the system. Covert misappropriation
attack against a networked control system has been considered in [7]. In this
type of attack, attacker needs the complete information of the system and more
resources compared to other types of attack. Other standard types of attacks
like Denial of Service (DoS) and deception attacks have been considered in
[8-10].

The problem of resilient state estimation in a standard form has been
addressed in [11-14]. They all assume that the attacker injects attack vector
to less than half of sensors, since it has been shown that it is impossible to
reconstruct the state of the system whenever more than half of the sensors are
under attack. Reconstructing the state of a nonlinear control system subject
to sensor and actuator attacks has also been studied in [15]. However, it is
assumed that some set of sensors are trustworthy. But, in the problem of
replay attack detection addressed in this paper, we assume that all sensors
can be compromised. It should be noted that our goal is only to detect the
presence of attack rather than estimating the true state of the system during
the attack time.

Replay attack against a centralized control system has been introduced in [16]
and [17] where for discrete time linear time invariant (LTI) control systems
an optimal noisy watermarking signal (independent and identically distributed
noise with zero mean and specific constant covariance) has been proposed in
order to detect the attack. In a recent article of the same authors [18], a
stationary process is utilized such that the defender has more degree of
freedom to design the watermarking signal for a discrete time control system,
which improves the detection rate with the same performance loss. The similar
iid approach has been proposed in [19] while the dynamic of the system is
unknown to the defender.

Detection of replay attack in frequency domain has been introduced in [20]. It
is claimed that the attack can be detected without injecting watermarking
signal. However, a restricting assumption is made so that a white Gaussian
noise signal with enough power should exist into the communication network
system to detect the replay attack. Moreover, it has been shown in [21] that
by including a nonlinear term in the control system in non-regular time
intervals, the attack can be detected. But, the performance loss induced by
incorporating such a nonlinear term into the system is not analyzed. The
problem of replay attack against multi agent systems has been studied in [22]
and a method based on sharing white Gaussian watermarking signal between
agents has been proposed. The authors also initiated the idea of dynamic
attack detector signal for the replay attack to the linear continuous time
systems in [23].

In this paper, we introduce a new method based on injection of dynamical
attack detector signal to detect the replay attack against discrete time LTI
stochastic systems. In the replay attack, it is assumed that there exists an
adversary who is capable of recording the sensory measurements for a period of
time, injecting a disruptive signal into the system and altering the real
sensory measurements of the system with the recorded measurements. It is also
assumed that the defender has the full information of the dynamic of the
system. We provide the conditions on the design parameters of the attack
detector signal so that the replay attack is detected while the system remains
stable in the absent of attack. Our approach yields the following advantages
and contributions with respect to the literature

  1. 1.

The proposed method in [16] and [17] are based on injection of the white and
colored noise to the system, respectively that render a constant probability
of attack detection. But in our method, the attack detector signal is
generated based on a dynamic that uses output estimation error. The
probability of attack detection in our method raises as the time increases and
the attack can be detected in a shorter period of time.

  2. 2.

We provide the condition under which the attack can be revealed using the
proposed attack detector signal (Theorem 1). It is also proved that by
applying the proposed attack detector signal, it is always possible for
physical systems to meet such condition (Theorem 2, proposition 3 and Remark
7).

  3. 3.

The detection rate of replay attack can be adjusted by design parameters in
our method and an optimization problem is proposed so that the performance
loss of the system after injection of the attack detector signal is minimized
subject to detection rate constraint.

  4. 4.

The performance of the proposed method is assessed from two aspects

     * •

It has been shown in simulation results that our proposed method is more
useful when defender wants to detect the attack with greater level of
accuracy, since within a finite time, our method can render more number of
alarms in an arbitrary window time.

     * •

Since we do not inject the noise directly to the system, the control signal in
our proposed attack detector has smaller range and lower energy with respect
to the control signal in white Gaussian noise method for the same detection
time and the same desired number of alarms.

##  II PRELIMINARIES

###  II-A System Model

Consider a control system, described by the following time invariant linear
system

$\begin{array}[]{l}x(t+1)=Ax(t)+Bu(t)+w(t)\\\ y(t)=Cx(t)+v(t)\end{array}$ (1)

where $x\in{\mathbb{R}^{n}}$, $u\in{\mathbb{R}^{m}}$, $y\in{\mathbb{R}^{p}}$
are the state, input and output vectors, respectively. Moreover,
$w\in{\mathbb{R}^{n}}$ and $v\in{\mathbb{R}^{p}}$ are the process and
measurement noises which are assumed to be independent white Gaussian
processes with covariance matrices $Q$ and $R$, respectively. Initial
condition $x(0)$ is a Gaussian random variable with mean $\mu$ and covariance
$\Sigma(0)$ (i.e. $x(0)\sim\mathcal{N}\big{(}\mu,\,\Sigma(0)\big{)}$) is
independent from process and measurement noises.

###  II-B Controller Structure

We assume that the pair $(A,B)$ is controllable and $(A,C)$ is observable.
Therefore, one can design a LQG controller that minimizes the following cost

$J=\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}{E\left\\{{{x^{T}}(t)Wx(t)+{u^{T}}(t)Uu(t)}\right\\}}$
(2)

where $W$ and $U$ are semi positive definite and positive definite weight
matrices with suitable dimensions, respectively. For this system, state
estimation unit based on the steady state Kalman filter is in the following
form

$\begin{split}\hat{x}(t)&=\hat{x}(t|t-1)+L\big{(}y(t)-C\hat{x}(t|t-1)\big{)}\\\
\hat{x}(t+1|t)&=A\hat{x}(t)+Bu(t)\end{split}$ (3)

In which $L={\Sigma_{e}}{C^{T}}{(C{\Sigma_{e}}{C^{T}}+R)^{-1}}$ and
${\Sigma_{e}}$ is stabilizing solution of the following Riccati equation

${\Sigma_{e}}=A{\Sigma_{e}}{A^{T}}-A{\Sigma_{e}}{C^{T}}{(C{\Sigma_{e}}{C^{T}}+R)^{-1}}C{\Sigma_{e}}{A^{T}}+Q$
(4)

and the control input to the system is the LQR gain multiplied by the
estimated state of the system as follows

$u(t)=K\hat{x}(t)$ (5)

where $K=-{({B^{T}}PB+U)^{-1}}{B^{T}}PA$ and $P$ is obtained from the
following Riccati equation

$P={A^{T}}PA+W-{A^{T}}PB{({B^{T}}PB+U)^{-1}}{B^{T}}PA$ (6)

##  III PROBLEM FORMULATION

###  III-A Attack Model

We consider there exists an adversary that disrupts the control system (1)
using replay attack. A replay attack model is defined in [16], where the
attacker targets sensory measurements of the system and meanwhile injects a
destructive signal into the system. The attack procedure is performed in two
phases

  1. 1.

During the first phase, from time $-\tau$ to $-1$, attacker records the
sensory measurements.

  2. 2.

In the second phase, from time zero to $\tau-1$, attacker replaces the sensory
measurements of the system with the recorded ones, as follows

$y(t)=y(t-\tau)\,\,\,\,\,\,for\,\,0\leq t\leq\tau-1$

In addition, during the second phase of the attack, the attacker injects a
disruptive signal ${f}(t)$ through the external input channel $B^{a}$ to the
system as follows

${\begin{array}[]{*{20}{c}}{x(t+1)=Ax(t)+Bu(t)+B^{a}f(t)}\\\
{y(t)=y(t-\tau)={y^{\tau}}(t)}\end{array}}$ (7)

The signal $f(t)$ is chosen by the attacker and this information is not
commonly available to the defender. It is worth mentioning that attacker does
not need to know about the model of the system. Having access to disclosure
resources in output channel and disruption resources in input and output
channels is sufficient to perform the replay attack.
We define $r(t)=y(t)-C\hat{x}(t|t-1)$ as the residue signal of the system.
When the system works on normal condition, the residue signal is an innovation
sequence with zero mean and covariance $H=C{\Sigma_{e}}{C^{T}}+R$ [24]. Let
define the following weighted norm of the residual signal

$g(t)={r^{T}}(t){{H}^{-1}}r(t)$ (8)

which has the Chi-square distribution form. It should be considered that the
defender utilizes $g(t)$ to detect the presence of abnormality in the system,
by comparing $g(t)$ with a scalar $\eta$ as a threshold which is determined by
defender. Anomaly detector will warn the system via an alarm if $g(t)>\eta$,
while the system is considered in its normal operation if $g(t)<\eta$. When
the system is in its normal condition, we define $\alpha=P(g(t)>\eta)$ as the
false alarm rate [25]. To reduce the number of false alarm rates in a period
of time, usually $\eta$ is adjusted in a way that $\alpha$ becomes close to
zero.
It is also assumed that there exists a window time $\mathscr{T}$ that defender
counts the number of alarms in the interval $(t-\mathscr{T},t)$ at each time
step $t$. If the number of alarms in a window time exceeds a threshold
(denoted by ${\beta}$) which is determined by the defender, it is considered
that the attack is occurred at time $t$. Since the expected number of false
alarms in a window time is equal to $\alpha\mathscr{T}$, in order to avoid
false attack detection in normal condition of the system, the desired number
of alarms in a window time should be chosen enough larger than
$\alpha\mathscr{T}$ (i.e. $\alpha\mathscr{T}\ll\beta$). Moreover, we define
the detection time as the required time to reach into the desired number of
alarms ${\beta}$.
Notation: In rest of the paper the superscript $\tau$ shows the delayed
signal. For instance
${x^{\tau}}(t)=x(t-\tau)$,
and superscript $a$ shows the under attack system’s signals. Therefore, from
(7) we get $y^{a}(t)=y^{\tau}(t)$.

######  Remark 1.

The amount of $\tau$ is determined by the attacker, and as much as it becomes
larger, the attacker will get enough time to disrupt the system and also the
defender will have enough time to detect the attack. However, in the case that
$\tau$ is too small, we may not be able to detect the attack for a given
$\beta$. In fact, the attack is detected if the detection time be lower than
$\tau$ for a given $\beta$.

The system from $-\tau$ to zero can be described in the following form

$\begin{split}{x^{\tau}}(t+1)&=A{x^{\tau}}(t)+B{u^{\tau}}(t)+{w^{\tau}}(t)\\\
{y^{\tau}}(t)&=C{x^{\tau}}(t)+{\upsilon^{\tau}}(t)\end{split}$ (9)

where $u^{\tau}(t)=K{\hat{x}^{\tau}}(t)$ and the state estimator of the system
in the absence of attack is as follow

${{\hat{x}}^{\tau}}(t+1|t)=A{{\hat{x}}^{\tau}}(t)+B{u^{\tau}}(t)$ (10)
${{\hat{x}}^{\tau}}(t)={{\hat{x}}^{\tau}}(t|t-1)+L[{y^{\tau}}(t)-C{{\hat{x}}^{\tau}}(t|t-1)]$
(11)

When the system is under attack, the state equation and its estimation are as
follow

$\begin{array}[]{l}{x^{a}}(t+1)={\rm A}{x^{a}}(t)+{\rm
B}{u^{a}}(t)+B^{a}f(t)+{w^{a}}(t)\end{array}$ (12)
${{\hat{x}}^{a}}(t+1|t)=A{{\hat{x}}^{a}}(t)+B{u^{a}}(t)$ (13)
${{\hat{x}}^{a}}(t)={{\hat{x}}^{a}}(t|t-1)+L[{y^{a}}(t)-C{{\hat{x}}^{a}}(t|t-1)]$
(14)

By subtracting (10) from (13), we get

$\begin{array}[]{l}{{\hat{x}}^{a}}(t+1|t)-{{\hat{x}}^{\tau}}(t+1|t)={A}({{\hat{x}}^{a}}(t)-{{\hat{x}}^{\tau}}(t))+\\\
+{B}({u^{a}}(t)-{u^{\tau}}(t))=(A+BK)(I-LC)({{\hat{x}}^{a}}(t|t-1)\\\
-{{\hat{x}}^{\tau}}(t|t-1))+(A+BK)L({y^{a}}(t)-{y^{\tau}}(t))\end{array}$ (15)

Since ${y^{a}}(t)={y^{\tau}}(t)$, (15) can be written as

$\begin{array}[]{l}{{\hat{x}}^{a}}(t+1|t)-{{\hat{x}}^{\tau}}(t+1|t)=({A}+BK)(I-LC)\\\
({{\hat{x}}^{a}}(t|t-1)-{{\hat{x}}^{\tau}}(t|t-1))\end{array}$ (16)

Let define $\Gamma\buildrel\Delta\over{=}({A}+BK)(I-LC)$. We know that
${A}+BK$ and $A(I-LC)$ are Schur individually, but the stability of
$({A}+BK)(I-LC)$ cannot be guaranteed. The residue signal of the system under
attack can be obtained as

$\begin{array}[]{l}{r^{a}}(t)={y^{a}}(t)-C{{\hat{x}}^{a}}(t|t-1)\\\
={y^{\tau}}(t)-C({{\hat{x}}^{\tau}}(t|t-1)+\Gamma{}^{t}\left({{{\hat{x}}^{a}}(0|-1)-{{\hat{x}}^{\tau}}(0|-1)}\right)\,)\\\
={r^{\tau}}(t)-C\Gamma{}^{t}\left({{{\hat{x}}^{a}}(0|-1)-{{\hat{x}}^{\tau}}(0|-1)}\right)\end{array}$
(17)

If $\Gamma$ is Schur, then the second term in (17) will vanish for large
enough $t$, leading to converges of ${r^{a}}(t)$ to ${r^{\tau}}(t)$. Thus, we
get

$\begin{array}[]{l}\mathop{\lim}\limits_{t\to\infty}P({g^{a}}(t)>\eta)=\mathop{\lim}\limits_{t\to\infty}P({r^{a}}^{T}(t){{H}^{-1}}{r^{a}}(t)>\eta)\\\
=\mathop{\lim}\limits_{t\to\infty}P({r^{\tau}}^{T}(t){{H}^{-1}}{r^{\tau}}(t)>\eta)=\alpha\end{array}$
(18)

According to (18), in the case the matrix $\Gamma$ is stable, the available
transient time for convergence of $r^{a}(t)$ to $r^{\tau}(t)$ is not reliable
for detection of the attack. Thus, the number of alarms occurred in a window
time during the attack can converge quickly to the number of alarms in the
normal condition of the system. As a result, it is probable that attack
remains stealthy in the transient time and this motivates the need for
designing an attack detector method, when $\Gamma$ is stable. However, if
$\Gamma$ is unstable, then both ${\hat{x}^{a}}$ and ${r^{a}}(t)$ diverge to
infinity and it is straightforward to show that the attack will get detected.

Note that from (12) to (22), both the state estimation and anomaly detector
units use counterfeit output measurement and not the real output of the
system. Consequently, $f(t)$ has no effect on the state estimation and anomaly
detector units.

##  IV ATTACK DETECTOR SIGNAL DYNAMICS

In the previous section, we studied that in the case matrix $\Gamma$ is
stable, the replay attack remains undetected. Here, we propose a method for
detection of the replay attack which causes the residue signal to become
unstable. Using this method in a finite time, the most achievable number of
alarms with arbitrary size can be obtained for a given window time
$\mathscr{T}$. Moreover, the detection rate can be controlled by the defender
and the design parameters can be adjusted to minimize the performance loss
caused by including the attack detector signal.
Let us assume that using the accessible information of model and input-output
data of the system, defender generates attack detector signal $\xi$ according
to the following process

$\displaystyle\zeta(t+1)=$ $\displaystyle\tilde{\rm A}\zeta(t)+\tilde{\rm
M}(y(t)-C\hat{x}(t))$ (19) $\displaystyle\xi(t)=$ $\displaystyle\tilde{\rm
K}\zeta(t)$

where the matrices $\tilde{\rm A}\in{\mathbb{R}^{n\times n}}$, $\tilde{\rm
M}\in{\mathbb{R}^{n\times p}}$ and $\tilde{\rm K}\in{\mathbb{R}^{m\times n}}$
are determined by the defender as design parameters. Fig. 1 shows a schematic
of under attack control system with the proposed attack detector generation
unit. Moreover, we assume that the defender uses the following input signal
that is composed of control signal and generated attack detector signal

$u(t)=K\hat{x}(t)+\xi(t)$ (20)

After injecting the attack detector signal into the system, the first
requirement is that the stability of the free of attack system is guaranteed.

######  Remark 2.

We assume that the system uses standard LQG controller to regulate the origin
in an optimal manner. Therefore, the gain of control signal and observer
parameters cannot be changed for the sake of attack detection.

######  Remark 3.

The attacker can cancel out the effect of attack detector signal if he or she
knows the exact model of the actual system and the structure of attack
detector dynamic. Moreover, it should have access to input-output measurement,
and the estimated state of the system to construct the negative of this
signal. However, clearly having access to such information is not
straightforward for the attacker. Moreover, we assume that the controller unit
is secure and the attacker is not able to make the attack detector signal
zero.

Figure 1: A schematic of under attack control system with a dynamical attack
detector signal generation unit. Attacker replays the previous output
measurements into the system. Output of the system and estimated state are
used for generating attack detector signal.

######  Proposition 1.

The stability of the system (1) after applying the control input (20) can be
guaranteed, if and only if the matrix $\tilde{\rm A}$ is Schur.

###### Proof.

The stability of the system can be verified by augmenting the dynamic of the
system, estimator error and attack detector signal as follows

$\tilde{x}(t+1)=\Theta\,\tilde{x}(t)+\Phi\,(t)$ (21)

where $e(t)=x(t)-\hat{x}(t)$ and

$\begin{array}[]{l}\tilde{x}(t)\buildrel\Delta\over{=}\left[{\begin{array}[]{*{20}{c}}{x(t)}\\\
{\zeta(t)}\\\ {e(t)}\end{array}}\right],\\\
\Theta\buildrel\Delta\over{=}\,\left[{\begin{array}[]{*{20}{c}}{A+BK}&{B\tilde{\rm
K}}&{-BK}\\\ 0&{\tilde{\rm A}}&{\tilde{\rm M}C}\\\
0&0&{(I-LC)A}\end{array}}\right],\\\
\Phi(t)\buildrel\Delta\over{=}\,\left[{\begin{array}[]{*{20}{c}}{w(t)}\\\
{\tilde{\rm M}\upsilon(t)}\\\
{(I-LC)w(t)-L\upsilon(t+1)}\end{array}}\right]\end{array}$

Since the matrix $\Theta$ is upper triangular, $A+BK$ and $(I-LC)A$ are Schur
therefore, the matrix $\Theta$ is Schur if and only if the matrix $\tilde{\rm
A}$ is Schur. ∎

######  Proposition 2.

Injection of the attack detector signal has no effect on the residue signal of
the system.

###### Proof.

We have assumed the control signal as $u(t)=K\hat{x}(t)+\xi(t)$. Now, let’s
define $\bar{e}(t)\buildrel\Delta\over{=}x(t)-\hat{x}(t|t-1)$. Then

$\displaystyle\bar{e}(t+1)$ $\displaystyle=x(t+1)-\hat{x}(t+1|t)$
$\displaystyle=Ax(t)+Bu(t)+w(t)-A\hat{x}(t)-Bu(t)$
$\displaystyle=A(x-\hat{x}(t))+w(t)$
$\displaystyle=A\left({x(t)-\left({I-LC}\right)\hat{x}(t|t-1)-Ly(t)}\right)+w(t)$
$\displaystyle=A\left({I-LC}\right)\bar{e}(t)+ALv(t)+w(t)$ (22)

and also we have

$\displaystyle r(t)$ $\displaystyle=y(t)-C\hat{x}(t|t-1)=$
$\displaystyle=C\left(x(t)-\hat{x}(t|t-1)\right)+v(t)=C\bar{e}(t)+v(t)$ (23)

Since $A(I-LC)$ is Schur, the expectation of $r(t)$ will be equal to zero in
steady state and also its covariance will be bounded. It can be seen that the
attack detector signal in control input has no effect on the residue signal in
the absence of attack because $u(t)$ is a known signal and it disappears in
the subtraction of $\hat{x}(t+1)$ from $x(t+1)$

∎

###  IV-A Detection of replay attack using the proposed attack detector method

In this section, we want to analyze the problem of attack detection using the
proposed dynamic attack detector. State estimation dynamic of the under attack
system after injection of the attack detector signal is obtained as follows

$\displaystyle{{\hat{x}}^{a}}(t+1|t)=$ $\displaystyle{\rm
A}{{\hat{x}}^{a}}(t)+{\rm B}(K{{\hat{x}}^{a}}(t)+{\xi^{a}}(t))$ (24)
$\displaystyle=$
$\displaystyle(A+BK)\left[{(I-LC){{\hat{x}}^{a}}(t|t-1)+L{y^{a}}(t)}\right]$
$\displaystyle+B\tilde{\rm K}{\zeta^{a}}(t)$

Moreover, the state estimation of the attack-free system is given by

$\displaystyle{{\hat{x}}^{\tau}}(t+1|t)=$ $\displaystyle{\rm
A}{{\hat{x}}^{\tau}}(t)+{\rm B}(K{{\hat{x}}^{\tau}}(t)+{\xi^{\tau}}(t))$ (25)
$\displaystyle=$
$\displaystyle(A+BK)\left[{(I-LC){{\hat{x}}^{\tau}}(t|t-1)+L{y^{\tau}}(t)}\right]$
$\displaystyle+B\tilde{\rm K}{\zeta^{\tau}}(t)$

Rewriting (19) for both under attack and attack free conditions, the
difference of the new equations is as follows

$\begin{array}[]{l}{\zeta^{a}}(t+1)-{\zeta^{\tau}}(t+1)=\tilde{\rm
A}({\zeta^{a}}(t)-{\zeta^{\tau}}(t))\\\ -\tilde{\rm
M}C(I-LC)\left[{{{\hat{x}}^{a}}(t|t-1)-{{\hat{x}}^{\tau}}(t|t-1))}\right]\end{array}$
(26)

Also subtracting (25) from (24) we get

$\begin{array}[]{l}{{\hat{x}}^{a}}(t+1|t)-{{\hat{x}}^{\tau}}(t+1|t)=\Gamma({{\hat{x}}^{a}}(t|t-1)-{{\hat{x}}^{\tau}}(t|t-1))\\\
+B\tilde{\rm K}({\zeta^{a}}(t)-{\zeta^{\tau}}(t))\end{array}$ (27)

where we used the fact that ${y^{a}}(t)={y^{\tau}}(t)$. By defining

$\sigma(t)=\left[{\begin{array}[]{*{20}{c}}{{\sigma_{1}}(t)}\\\
{{\sigma_{2}}(t)}\end{array}}\right]\buildrel\Delta\over{=}\left[{\begin{array}[]{*{20}{c}}{{{\hat{x}}^{a}}(t|t-1)-{{\hat{x}}^{\tau}}(t|t-1)}\\\
{{\zeta^{a}}(t)-{\zeta^{\tau}}(t)}\end{array}}\right]$ (28)

(26) and (27) can be represented as the following state space model

$\sigma(t+1)=\underbrace{\left[{\begin{array}[]{*{20}{c}}\Gamma&{B\tilde{\rm
K}}\\\ {-\tilde{\rm M}C(I-LC)}&{\tilde{\rm
A}}\end{array}}\right]}_{\Delta}\sigma(t)$ (29)

in which $\Delta$ stands for the attack detection matrix. In what follows, we
study the condition that replay attack using the proposed attack detector
signal can be detected.

######  Theorem 1.

If the matrix $\Delta$ has an eigenvalue with modulus greater than one, then
the attack will be detected.

###### Proof.

the residue signal of the under attack system is

$\begin{split}{r^{a}}(t)=&{y^{a}}(t)-C{{\hat{x}}^{a}}(t\left|{t-1}\right.)={y^{\tau}}(t)\\\
&-C({\sigma_{1}}(t)+{{\hat{x}}^{\tau}}(t\left|{t-1}\right))={r^{\tau}}(t)-C{\sigma_{1}}(t)\end{split}$
(30)

Rewriting (8) for under attack condition we get

$\begin{split}{g^{a}}(t)=&{r^{a}}^{T}(t){H^{-1}}{r^{a}}(t)={\left({{r^{\tau}}(t)-C{\sigma_{1}}(t)}\right)^{T}}{H^{-1}}\\\
&\times\left({{r^{\tau}}(t)-C{\sigma_{1}}(t)}\right)={r^{\tau}}^{T}(t){H^{-1}}{r^{\tau}}(t)\\\
&-\sigma_{1}^{T}(t){C^{T}}{H^{-1}}{r^{\tau}}(t)-{r^{\tau}}^{T}(t){H^{-1}}C{\sigma_{1}}(t)\\\
&+\sigma_{1}^{T}(t){C^{T}}{H^{-1}}C{\sigma_{1}}(t)\end{split}$ (31)

It is straightforward to show that for unstable $\sigma_{1}(t)$ we have
$\mathop{\lim}\limits_{t\to\infty}P({g^{a}}(t)>\eta)=1$ and therefore, the
attack will be detected. Now we show that if the matrix $\Delta$ has at least
one eigenvalue outside the unit circle, then ${\sigma_{1}}(t)$ will be
unstable. Assuming the matrix $\Delta$ has $n$ distinct eigenvalues, the
solution of the system (29) can be obtained as

$\sigma(t)={c_{1}}\lambda_{1}^{t}{p_{1}}+...+{c_{n}}\lambda_{n}^{t}{p_{n}}$
(32)

where ${p_{i}}$ is the $i$-th right eigenvector of $\Delta$ and ${c_{i}}$ is a
constant equal to the inner product of the $i$-th left eigenvector and initial
condition $\sigma(0)$. Assume that ${\lambda_{i}}$ is an eigenvalue outside
the unit circle. It is clear that if at least one of the upper $n$ rows of the
vector ${p_{i}}$ is nonzero, then ${\sigma_{1}}(t)$ will be unstable (note
that since $\sigma(0)$ is taken from a stochastic distribution, ${c_{i}}$ is
almost surely a nonzero constant). Now, assume for the sake of contradiction
the case that the $n$ upper rows of the vector ${p_{i}}$ are zero. According
to eigenvalue definition for ${\lambda_{i}}$ we have

$\left[{\begin{array}[]{*{20}{c}}\Gamma&{B\tilde{\rm K}}\\\ {-\tilde{\rm
M}C(I-LC)}&{\tilde{\rm
A}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}0\\\
{{p_{i,2}}}\end{array}}\right]={\lambda_{i}}\left[{\begin{array}[]{*{20}{c}}0\\\
{{p_{i,2}}}\end{array}}\right]$ (33)

where ${p_{i,2}}$ represents the lower $n$ rows of the vector ${p_{i}}$. From
the above equation we have

$\tilde{\rm A}{p_{i,2}}={\lambda_{i}}{p_{i,2}}$ (34)

From (34), it can be observed that ${\lambda_{i}}$ and ${p_{i,2}}$ are the
eigenvalue and eigenvector associated with the matrix $\tilde{\rm A}$,
respectively. But it contradicts to the fact that all the eigenvalues of
$\tilde{\rm A}$ are inside the unit circle. Therefore, we can conclude that
the upper $n$ rows of eigenvector ${p_{i}}$ associated with the unstable
eigenvalue cannot be all equal to zero. ∎

In [17], the detectability of attack has been verified when the matrix
$\Gamma=(A+BK)(I-LC)$ is unstable. In our method, by adding the attack
detector dynamic as expanding the dynamics of the underlying system, the
matrix $\Gamma$ is extended to matrix $\Delta$. In addition to increasing
probability of attack detection in the proposed dynamic attack detection
method, the extended matrix $\Delta$ contains some designed parameters which
can be used to optimize the performance of detection method.

######  Remark 4.

Even if we have an eigenvalue with some multiplicity, there exists at least
one independent eigenvector associated with that eigenvalue. Therefore, the
discussion on the independent eigenvectors which is presented in proof of
Theorem 1 can be applied since, we do not need to use the other generalized
eigenvectors for this purpose and the proof remains unchanged.

######  Remark 5.

In theorem 1, from $\mathop{\lim}\limits_{t\to\infty}P({g}(t)>\eta)=1$, it can
be deduced that in infinite time the anomaly detector gives the most
achievable number of alarms $\beta_{max}$ (continues alarm) in a given window
time, nevertheless, it is sufficient that the number of alarms exceeds a
threshold to detect the attack. Thus, if this threshold is less than
$\beta_{max}$, the attack can be revealed in a finite time.

In the following theorem we show that according to the degrees of freedom on
choosing the matrices $\tilde{\rm M}$, $\tilde{\rm K}$ and $\tilde{\rm A}$,
the defender can make the matrix $\Delta$ unstable.

######  Theorem 2.

If the matrix $\Gamma$ is full rank, by adjusting the matrices $\tilde{\rm
M}$, $\tilde{\rm K}$ and stable $\tilde{\rm A}$, the matrix $\Delta$ will have
at least an eigenvalue with modulus greater than one.

###### Proof.

We know that the product of all eigenvalues of a matrix is equal to the
determinant of that matrix. If $\det(\Delta)>1$, then there exists at least
one eigenvalue of $\Delta$ that is outside of the unit circle. In what
follows, we prove that if the matrix $\Gamma$ is full rank, by adjusting the
design matrices $\tilde{\rm M}$, $\tilde{\rm K}$ and $\tilde{\rm A}$ one can
make the determinant of the matrix $\Delta$ greater than one. We know that if
$\det(\Gamma)\neq 0$, we have

$\displaystyle\begin{array}[]{l}\det\left({\left[{\begin{array}[]{*{20}{c}}\Gamma&{B\tilde{\rm
K}}\\\ {-\tilde{\rm M}C(I-LC)}&{\tilde{\rm
A}}\end{array}}\right]}\right)=\det\left(\Gamma\right)\\\
\times\det\left({\tilde{\rm A}+\tilde{\rm M}C(I-LC){\Gamma^{-1}}B\tilde{\rm
K}}\right)\end{array}$ (39)

Since $\Gamma$ is Schur, the absolute of its determinant is lower than one.
Assume that the element $(i,j)$ of $\tilde{\rm M}$ and the element $(r,s)$ of
$\tilde{\rm K}$ get two arbitrary values $m$ and $k$, respectively, and the
other elements of $\tilde{\rm M}$ and $\tilde{\rm K}$ are set to zero. It is
obvious that only $i$-th row of $\tilde{\rm M}C(I-LC)$ and $s$-th column of
$B\tilde{\rm K}$ remain nonzero. Therefore, the matrix $\tilde{\rm
M}C(I-LC){\Gamma^{-1}}B\tilde{\rm K}$ will have only one nonzero element in
its $(i,s)$ entry, that can be chosen in a desired way. The determinant of
$\tilde{\rm A}$ along the $i$-th row appears as a linear combination of the
determinants of the minors $M_{ij}$, which cannot be all equal to zero if
$\tilde{\rm A}$ is chosen an invertible matrix. Let ${a_{is}}$ be the nonzero
element associated with such a nonzero determinant, therefore we have

$\displaystyle\det(\tilde{\rm A}+\tilde{\rm
M}C(I-LC){\Gamma^{-1}}B\tilde{K})={\mu_{1}}$ (40)
$\displaystyle+{\mu_{2}}mk\det({M_{is}})$

where $\mu_{1}$ and $\mu_{2}$ are constants. Then, one can choose $m$ and $k$
in a way that $det(\Delta)$ is greater than one. ∎

######  Proposition 3.

The matrix $\Gamma$ is full rank if and only if the matrix $A$ is full rank.

###### Proof.

We know that the multiplication of two full rank matrices is a full rank
matrix and vice versa. Therefore, the matrix $\Gamma$ is full rank if and only
if $({A}+BK)$ and $(I-LC)$ are full rank. First we show the matrix $(I-LC)$ is
always full rank. Also we know that the following equality holds for any
$M\in{\mathbb{R}^{p\times k}}$ and $N\in{\mathbb{R}^{k\times p}}$

$\det({I_{p\times p}}-MN)=\det({I_{k\times k}}-NM)$ (41)

Therefore, using the fact that
$L={\Sigma_{e}}{C^{T}}{(C{\Sigma_{e}}{C^{T}}+R)^{-1}}$, $\det({I_{n\times
n}}-LC)$ can be written as

$\begin{array}[]{l}\det\left({{I_{n\times
n}}-LC}\right)=\det\left({{I_{p\times p}}-CL}\right)=\\\
\det\left({I-C{\Sigma_{e}}{C^{T}}{{\left({C{\Sigma_{e}}{C^{T}}+R}\right)}^{-1}}}\right)=\\\
\det\left({\left({C{\Sigma_{e}}{C^{T}}+R-C{\Sigma_{e}}{C^{T}}}\right){{\left({C{\Sigma_{e}}{C^{T}}+R}\right)}^{-1}}}\right)=\\\
\det\left({R{{\left({C{\Sigma_{e}}{C^{T}}+R}\right)}^{-1}}}\right)=\det\left(R\right)\\\
\times\det\left({{{\left({C{\Sigma_{e}}{C^{T}}+R}\right)}^{-1}}}\right)\end{array}$
(42)

Since the matrix $R$ is full rank and positive definite, and also the matrix
$C{\Sigma_{e}}{C^{T}}$ is a positive semi-definite matrix, therefore,
$\det({I_{n\times n}}-LC)>0$.
For the matrix $A+BK$, we can follow the same procedure. $\det(A+BK)$ can be
written in the following form

$\begin{split}&\det({A+BK})=\det\left({A-B{{\left({{B^{T}}PB+U}\right)}^{-1}}{B^{T}}PA}\right)\\\
&=\det(A)\det\left({{I_{n\times
n}}-B{{\left({{B^{T}}PB+U}\right)}^{-1}}{B^{T}}P}\right)\\\
&=\det(A)\det\left({{I_{m\times
m}}-{B^{T}}PB{{\left({{B^{T}}PB+U}\right)}^{-1}}}\right)=\\\
&\det(A)\det\bigg{(}{\Big{(}{{B^{T}}PB+U-{B^{T}}PB}\Big{)}{{\Big{(}{{B^{T}}PB+U}\Big{)}}^{-1}}}\bigg{)}\\\
&=\det(A)\det\left({U{{\left({{B^{T}}PB+U}\right)}^{-1}}}\right)=\det(A)\det(U)\times\\\
&\,\,\,\,\,\,\,\det\left({{{\left({{B^{T}}PB+U}\right)}^{-1}}}\right)\end{split}$
(43)

Since the matrix $U$ is positive definite and also the ${B^{T}}PB$ is positive
semi definite, therefore,
$\det(U)\det\left({{{\left({{B^{T}}PB+U}\right)}^{-1}}}\right)>0$. As a
result, $A+BK$ is full rank if and only if the matrix $A$ is full rank and the
proof is completed. ∎

######  Remark 6.

The cyber-physical systems which are under study of this paper are continuous
time in reality. It is straightforward to show that the matrix $A$ which can
be obtained after discretization of any continuous time dynamic is full rank
and therefore, the condition in Proposition 3 is always satisfied.

Now, among all the feasible values for design matrices $\tilde{\rm M}$,
$\tilde{\rm K}$ and stable $\tilde{\rm A}$, we seek for those with minimum
performance loss effect in the absence of attack. In what follows we obtain
the performance loss of the system after injection of the proposed attack
detector signal.

###  IV-B Performance Loss

Here, we follow the procedure for obtaining the performance loss (new cost
function) after injection of the attack detector signal.
The input of the system after inserting attack detector signal is as follows

$u(t)=K\hat{x}(t)+\xi(t)=K\hat{x}(t)+\tilde{\rm K}\zeta(t)$ (44)

The performance loss is defined as

$\tilde{J}=\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}{E\left\\{{{x^{T}}(t)Wx(t)+{u^{T}}(t)Uu(t)}\right\\}}$
(45)

Therefore, after injection of the attack detector signal to the system, the
performance loss can be written as

$\begin{split}\tilde{J}=&\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}E\bigg{\\{}{x^{T}}(t)Wx(t)+\Big{(}K\big{(}x(t)-e(t)\big{)}\\\
&+\tilde{\rm K}\zeta(t)\Big{)}^{T}U\Big{(}K\big{(}x(t)-e(t)\big{)}+\tilde{\rm
K}\zeta(t)\Big{)}\bigg{\\}}=\\\
=&\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}{E\bigg{\\{}{{\left[{\begin{array}[]{*{20}{c}}{{x^{T}}(t)}&{{\zeta^{T}}(t)}&{{e^{T}}(t)}\end{array}}\right]}}}\\\
=&\underbrace{\begin{bmatrix}W+{K^{T}}UK&{{K^{T}}U\tilde{\rm
K}}&{-{K^{T}}UK}\\\ {\tilde{\rm K}^{T}UK}&{\tilde{\rm K}^{T}U\tilde{\rm
K}}&{-\tilde{\rm K}^{T}UK}\\\ {-{K^{T}}UK}&{-{K^{T}}U\tilde{\rm
K}}&{{K^{T}}UK}\end{bmatrix}}_{G}\begin{bmatrix}{x(t)}\\\ {\zeta(t)}\\\
{e(t)}\end{bmatrix}\bigg{\\}}\end{split}$ (46)
$\begin{split}\tilde{J}=&\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}{E\left\\{{tr\left({{{\tilde{x}}^{T}}(t)G\tilde{x}(t)}\right)}\right\\}}\\\
=&\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}{E\left\\{{tr\left({\tilde{x}(t){{\tilde{x}}^{T}}(t)G}\right)}\right\\}}\\\
=&\mathop{\lim}\limits_{N\to\infty}\frac{1}{N}\sum\limits_{t=1}^{N}{tr\left({G{\mathop{\rm
cov}}(\tilde{x}(t))}\right)}=tr\left({G{\mathop{\rm
cov}}(\tilde{x})}\right)\,\end{split}$ (47)

where according to (21), $cov(\tilde{x})$ can be obtained by solving the
following Lyapunov equation

${\mathop{\rm cov}}(\tilde{x})=\Theta{\mathop{\rm
cov}}(\tilde{x}){\Theta^{T}}+E\\{\Phi\Phi{\,{}^{T}}\\}$ (48)

By taking the expectation values of the elements of $\Phi{\,{}^{T}}\Phi$ we
get

$\begin{array}[]{l}{\mathop{\rm cov}}(\tilde{x})=\Theta{\mathop{\rm
cov}}(\tilde{x}){\Theta^{T}}+\\\
\underbrace{\left[{\begin{array}[]{*{20}{c}}Q&0&{Q{{(I-LC)}^{T}}}\\\
0&{\tilde{\rm M}R{{\tilde{\rm M}}^{T}}}&0\\\
{(I-LC)Q}&0&{(I-LC)Q{{(I-LC)}^{T}}+LR{L^{T}}}\end{array}}\right]}_{\Psi}\end{array}$

(49)

###  IV-C Optimization Problem

The optimal values for $\tilde{\rm M}$, $\tilde{\rm K}$ and $\tilde{\rm A}$
can be found from the following optimization problem

$\begin{array}[]{l}\mathop{\min}\limits_{\tilde{\rm A},\,\,\tilde{\rm
M},\,\tilde{\rm K}}tr\left({G{\mathop{\rm cov}}(\tilde{x})}\right)\\\
st.\,\,\,\,{\mathop{\rm cov}}(\tilde{x})=\Theta{\mathop{\rm
cov}}(\tilde{x}){\Theta^{T}}+\Psi\\\ 1<\delta\leq\max\\{abs(eig(\Delta))\\}\\\
\tilde{\rm A}:Schur\,\,and\,\,full\,rank\end{array}$ (50)

######  Proposition 4.

The optimization problem (50) is feasible.

###### Proof.

Since the matrix $\Theta$ is Schur, according to Proposition 1 and also due to
the fact that the matrix $\Psi$ in the right hand side of (49) is a non-
negative definite matrix, there exists a unique non-negative definite solution
to the Lyapunov equation (49). On the other hand, in Theorem 2, we showed that
one can adjust the absolute value of the maximum eigenvalue of the matrix
$\Delta$ to become greater than any scalar $\delta>1$. Therefore, the feasible
domain of the optimization problem (50) is non-empty. ∎

It should be mentioned that the optimization problem (50) is a non-convex
nonlinear programming. Nevertheless, there are many nonlinear programming
methods for solving non-convex problems such as interior point, trust-region-
reflective, etc.

##  V SIMULATION RESULTS

In this section, DC motor systems that are commonly used in cyber-physical
systems are considered to verify the results of the paper. By assuming that
the magnetic field is constant, the state space model of the system is LTI and
can be represented as follows

$\frac{d}{{dt}}\left[{\begin{array}[]{*{20}{c}}\theta\\\ {\dot{\theta}}\\\
i\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}0&1&0\\\
0&{-\frac{b}{J}}&{\frac{K_{t}}{J}}\\\
0&{-\frac{K_{b}}{L}}&{-\frac{R}{L}}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}\theta\\\
{\dot{\theta}}\\\ i\end{array}}\right]+\left[{\begin{array}[]{*{20}{c}}0\\\
0\\\ {\frac{1}{L}}\end{array}}\right]V(t)$ (51)
$y(t)=\left[{\begin{array}[]{*{20}{c}}1&0&0\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}\theta\\\
{\dot{\theta}}\\\ i\end{array}}\right]$ (52)

The physical parameters for this example are given in Table 1. It is assumed
that the model (51) is discretized with sampling time $T_{s}=.01s$. The cost
matrices for this system, W and U, and covariance matrices, Q and R, are equal
to the identity matrix with appropriate dimensions. The cost of LQG controller
is $j^{*}=1.04$ and we set the scalar $\gamma$ in a way that $\alpha=.01$. We
fix the window time $\mathscr{T}=1s$ (equal to 100 sample time). If the number
of alarms occurred in this window time exceed the desired number of alarms
$\beta$, then the attack will be revealed.

It is assumed that the attacker records output measurements from time interval
$[-15\,\ 0]s$ and replays them exactly in time interval $[0\,\ 15]s$. For this
system, the matrix $\Gamma$ has three eigenvalues inside the unit circle.
Consequently, the attack will remain stealthy and the defender has to use
another countermeasure against the attack.

We have designed our proposed attack detector signal by solving the
optimization problem (50) for $\delta=1.03$. The design matrices are as
follows

$\begin{split}\tilde{\rm
A}&=\left[{\begin{array}[]{*{20}{c}}{.48}&{-.81}&{.02}\\\
{.01}&{.61}&{-.92}\\\ {.89}&{.73}&{-.9}\end{array}}\right],\quad\tilde{\rm
M}=\left[{\begin{array}[]{*{20}{c}}{-.84}\\\ {-.49}\\\
{-.81}\end{array}}\right]\\\ \tilde{\rm
K}&=\left[{\begin{array}[]{*{20}{c}}{.9}&{-.1}&{.35}\end{array}}\right]\end{split}$
(53)

Since the system behaves stochastically, we use expected detection time (the
average time that is needed to reach a desired number of alarms $\beta$)
instead of detection time. In Fig. 2, we have compared the expected detection
time in our proposed method with white Gaussian watermarking signal method in
the case that performance loss is equal for both methods. Each individual
curve is obtained by averaging on 1000 experiments.
In Fig. 2, it can be seen that our proposed method gives lower expected
detection time for all desired number of alarms. However, for $\beta\geq 5$
the White Gaussian is not capable of detecting the presence of attack but the
proposed attack detector method can even reveal the attack for the most
achievable number of alarms $\beta_{max}=100$.

Figure 2: Expected detection time versus number of alarms for our proposed method with $\delta=1.03$ and white Gaussian watermarking signal for the same performance loss. Table I: Physical parameters of DC motor Parameter | Value
---|---
$b$ | $3.5\times{10^{-6}}\,\,N.m.s$
$J$ | $3.23\times{10^{-6}}\,\,kg.m^{2}$
$L$ | 2.75 $\times{10^{-6}}\,\,H$
$K_{b}$ | $0.0274\,\,V/rad/sec$
$K_{t}$ | $0.0274\,\,N.m/Amp$
$R$ | $4\,\,\,Ohm$

Figure 3: Control signal after injection of our proposed (blue line) and
watermarking signal (red line).

Figure 4: The amount of performance loss versus $\delta$. By increasing
$\delta$ the amount of performance loss will also increase.

In Fig. 3 we compare the control signal of our proposed method with control
signal of white Gaussian watermarking method in the case that they give the
same detection time ($.5s$) for the same desired number of alarms $\beta=30$.
As it can be seen, the control signal in our proposed attack detector is
smoother, has smaller range and lower energy with respect to the control
signal in white Gaussian noise method.

Moreover, Fig. 4 examines the relation between the amount of performance loss
and maximum eigenvalue of detection matrix $\Delta$. It is clear that by
increasing $\delta$ the amount of performance loss increases as well.

## References

  * [1] Neuman, Clifford. ”Challenges in security for cyber-physical systems.” DHS Workshop on Future Directions in Cyber-Physical Systems Security. 2009.
  * [2] Langner, Ralph. ”To kill a centrifuge: A technical analysis of what stuxnet’s creators tried to achieve.” The Langner Group (2013).
  * [3] S. Karnouskos, ”Stuxnet worm impact on industrial cyber-physical system security,” in Proc. 37th Annu. Conf. IEEE Industrial Electronics Society, Melbourne, Australia, 2011, pp. 4490-4494.
  * [4] Liu, Yao, Peng Ning, and Michael K. Reiter. ”False data injection attacks against state estimation in electric power grids.” ACM Transactions on Information and System Security (TISSEC) 14.1 (2011): 13.
  * [5] Pasqualetti, Fabio, Florian Dorfler, and Francesco Bullo. ”Attack detection and identification in cyber-physical systems.” IEEE Transactions on Automatic Control 58.11 (2013): 2715-2729.
  * [6] Teixeira, Andre, et al. ”A secure control framework for resource-limited adversaries.” Automatica 51 (2015): 135-148.
  * [7] Smith, Roy S. ”Covert misappropriation of networked control systems: Presenting a feedback structure.” IEEE Control Systems 35.1 (2015): 82-92.
  * [8] Yuan, Yuan, et al. ”Resilient control of cyber-physical systems against denial-of-service attacks.” 2013 6th International Symposium on Resilient Control Systems (ISRCS). IEEE, 2013.
  * [9] Pang, Zhong-Hua, and Guo-Ping Liu. ”Design and implementation of secure networked predictive control systems under deception attacks.” IEEE Transactions on Control Systems Technology 20.5 (2011): 1334-1342.
  * [10] Ding, Derui, et al. ”Event-based security control for discrete-time stochastic systems.” IET Control Theory & Applications 10.15 (2016): 1808-1815.
  * [11] Fawzi, Hamza, Paulo Tabuada, and Suhas Diggavi. ”Secure estimation and control for cyber-physical systems under adversarial attacks.” IEEE Transactions on Automatic Control 59.6 (2014): 1454-1467.
  * [12] Shoukry, Yasser, and Paulo Tabuada. ”Event-triggered state observers for sparse sensor noise/attacks.” IEEE Transactions on Automatic Control 61.8 (2016): 2079-2091.
  * [13] Shoukry, Yasser, et al. ”SMT-based observer design for cyber-physical systems under sensor attacks.” ACM Transactions on Cyber-Physical Systems 2.1 (2018): 5.
  * [14] Pajic, Miroslav, Insup Lee, and George J. Pappas. ”Attack-resilient state estimation for noisy dynamical systems.” IEEE Transactions on Control of Network Systems 4.1 (2017): 82-92.
  * [15] Nateghi, Shamila, et al. ”Cyber Attack Reconstruction of Nonlinear Systems via Higher-Order Sliding-Mode Observer and Sparse Recovery Algorithm.” 2018 IEEE Conference on Decision and Control (CDC). IEEE, 2018.
  * [16] Mo, Yilin, and Bruno Sinopoli. ”Secure control against replay attacks.”Communication, Control, and Computing, 2009. Allerton 2009. 47th Annual Allerton Conference on. IEEE, 2009.
  * [17] Mo, Yilin, Rohan Chabukswar, and Bruno Sinopoli. ”Detecting integrity attacks on SCADA systems.” IEEE Transactions on Control Systems Technology 22.4 (2014): 1396-1407.
  * [18] Mo, Yilin, Sean Weerakkody, and Bruno Sinopoli. ”Physical authentication of control systems: designing watermarked control inputs to detect counterfeit sensor outputs.” IEEE Control Systems 35.1 (2015): 93-109.
  * [19] Liu, Hanxiao, et al. ”An On-line Design of Physical Watermarks.” 2018 IEEE Conference on Decision and Control (CDC). IEEE, 2018.
  * [20] Tang, Bixiang, Luis D. Alvergue, and Guoxiang Gu. ”Secure networked control systems against replay attacks without injecting authentication noise.” American Control Conference (ACC), 2015. IEEE, 2015.
  * [21] Hoehn, Andreas, and Ping Zhang. ”Detection of replay attacks in cyber-physical systems.” American Control Conference (ACC), 2016. IEEE, 2016.
  * [22] Khazraei, Amir, Hamed Kebriaei, and Farzad Rajaei Salmasi. ”Replay attack detection in a multi agent system using stability analysis and loss effective watermarking.” American Control Conference (ACC), 2017. IEEE, 2017.
  * [23] Khazraei, Amir, Hamed Kebriaei, and Farzad Rajaei Salmasi. ”A new watermarking approach for replay attack detection in LQG systems.” Decision and Control (CDC), 2017 IEEE 56th Annual Conference on. IEEE, 2017.
  * [24] Simon, Dan. Optimal state estimation: Kalman, H infinity, and nonlinear approaches. John Wiley & Sons, 2006.
  * [25] Willsky, Alan S. ”A survey of design methods for failure detection in dynamic systems.” Automatica 12.6 (1976): 601-611.


# Inhalation and deposition of spherical and pollen particles after middle
turbinate resection in a human nasal
cavity111http://doi.org/10.1016/j.resp.2021.103769

Kiao Inthavong  Yidan Shang  John M. DelGaudio  Sarah K Wise  Thomas S Edwards
Kimberly Bradshaw  Eugene Wong  Murray Smith  Narinder Singh  Mechanical &
Automotive Engineering, School of Engineering, RMIT University, Bundoora,
Victoria 3083, Australia  Department of Otolaryngology-Head and Neck Surgery
Emory University, Atlanta, Georgia USA  Faculty of Medicine & Health, The
University of Sydney, NSW 2006, Australia  Department of Otolaryngology, Head
and Neck Surgery, Westmead Hospital, Westmead, NSW 2145, Australia

###### Abstract

Middle turbinate resection significantly alters the anatomy and redistributes
the inhaled air. The superior half of the main nasal cavity is opened up,
increasing accessibility to the region. This is expected to increase
inhalation dosimetry to the region during exposure to airborne particles. This
study investigated the influence of middle turbinate resection on the
deposition of inhaled pollutants that cover spherical and non-spherical
particles (e.g. pollen). A computational model of the nasal cavity from CT
scans, and its corresponding post-operative model with virtual surgery
performed was created. Two constant flow rates of 5L/min, and 15L/min were
simulated under a laminar flow field. Inhaled particles including pollen (non-
spherical), and a spherical particle with reference density of 1000kg/m3 were
introduced in the surrounding atmosphere. The effect of surgery was most
prominent in the less patent cavity side, since the change in anatomy was
proportionally greater relative to the original airway space. The left cavity
produced an increase in particle deposition at a flow rate of 15L/min. The
main particle deposition mechanisms were inertial impaction, and to a lesser
degree gravitational sedimentation. The results are expected to provide
insight into inhalation efficiency of different aerosol types, and the
likelihood of deposition in different nasal cavity surfaces.

###### keywords:

nasal cavity , deposition , inhalability , pollen , drag coefficient ,
computational fluid dynamics , CFD

††journal: Accepted in Respiratory Physiology & Neurobiology 294:103769

##  1 Introduction

Nasal airway obstruction can lead to significant impairment in quality-of-
life, secondary to a myriad of potential sequelae, including chronic nasal
obstruction, mouth breathing, chronic rhinosinusitis, sleep-disordered
breathing and obstructive sleep apnoea. As the prevalence of nasal airway
obstruction is high, surgery to relieve nasal obstruction is a commonly
performed elective procedure. Post-surgery patients have significantly altered
anatomy and simulating airflow will improve our understanding of how surgical
strategies affect post-surgical airway ventilation, as well as the inhalation
dosimetry caused by exposure to airborne pollutants. The removal of anatomical
structures in the nose opens up the airway, redistributing the airflow and
consequently the particle distribution.

There have been extensive computational studies investigating different nasal
surgeries and their effect on nasal airflow [Xiong et al., 2008, Rhee et al.,
2011, De Wang et al., 2012, Frank et al., 2013a, Hariri et al., 2015, Dayal et
al., 2016, Maza et al., 2019, Inthavong et al., 2019, Siu et al., 2020a],
which showed the ability of virtual surgery to guide surgical interventions to
alleviate nasal airway obstruction. Previous CFD studies of middle turbinate
resection [Zhao et al., 2014, Alam et al., 2019, Lee et al., 2016] have
demonstrated shifts in airflow toward the area of middle turbinate resection
and decreased airflow velocity, decreased wall shear stress, and increased
local air pressure.

Chen et al. [2012] found that following functional endoscopic sinus surgery, a
moderate inspiratory airflow rate and a particle diameter of approximately
10$\upmu$m improved intranasal deposition into the maxillary sinuses. Frank et
al. [2013b] investigated surgical correction of nasal anatomic deformities on
ten patients, and their simulation results showed posterior particle
deposition after surgery increased by 118%. Bahmanzadeh et al. [2015] showed
that sphenoidotomy increased the airflow and particle deposition into the
sphenoid sinus, with the highest deposition occurring for resting breathing
rate with $10\upmu$m. Siu et al. [2020b] studied drug delivery in post-
operative sinonasal geometries of patients who underwent a comprehensive
functional endoscopic sinus surgery and a modified endoscopic Lothrop
procedure and found sinus deposition was more effective with inhalational
transport of low-inertia particles outside of the range produced by many
standard nasal sprays or nebulizer.

A recent outcome from DelGaudio et al. [2019]’s analysis of patients with
middle turbinate resection suggested that the surgery removed the middle
turbinate protective function in preventing inhaled allergen particles from
depositing on the septum. Anatomical changes in the airway aimed at restoring
airway ventilation or increased accessibility to the paranasal sinuses as the
primary aims, may in fact be ignoring the secondary effects of persistent
inhalation exposure to airborne particulates, such as dust [Tian et al.,
2007], dust mites [Zhang et al., 2019], and fibres. Additionally, pollen has
been implicated for initiating the onset of allergic rhinitis (such as hay
fever) and sometimes has exacerbating asthma [Bousquet et al., 2001].

Modelling of non-spherical particles needs to account for the changes in
particle motion due to the different drag coefficients leading to altered
trajectories. While the discrete element method and immersed boundary methods
can resolve more physics such as rotational degrees-of-freedom and inter-
particle interactions, the Lagrangian approach remains popular for its
computational efficiency and ability to easily integrate into finite-volume
CFD codes. Drag correlations for non-spherical particles have been developed
Hölzer and Sommerfeld [2008], Wang et al. [2018], Bagheri and Bonadonna
[2016]. This is in addition to inclusions of lift and torque components
applied to oblate ellipsoidals or fibres [Sanjeevi et al., 2018, Zastawny et
al., 2012, Tian et al., 2012, Tian and Ahmadi, 2013]. Th authors used a
simplified approach in earlier work Inthavong et al. [2013] by using shape
factors and an overall drag coefficient to account for tumbling and rotational
effects, but this does not allow for the interception deposition mechanism.

This study investigated the changes in deposition of inhaled pollutants in the
nasal airways, and the proportion escaping through the nasopharynx towards the
lungs. Spheres with aerodynamic equivalent diameters, $d_{ae}$, were used as a
reference while pollen particles were investigated. A nasal airway model was
used in pre-operative state, followed by a virtual middle turbinectomy surgery
(post-operative state) with two constant inhalation flow rates of 5 and 15
L/min. The airway was separated into different regions allowing a detailed
comparison of regional deposition changes.

##  2 Method

###  2.1 Nasal cavity model

The nasal cavity geometry was reconstructed from high-quality CT scans of a
healthy 25-year-old, Asian female (161cm height, 53kg mass) without any nasal
obstructions. The CT data was scanned with a Siemens Dual Source CT Machine
with parameters of $0.39\times 0.39$mm in pixel size, $512\times 512$ in image
dimension, and 0.5mm in slice thickness. The images were imported into 3D
Slicer, where the entire paranasal sinus model was determined and exported as
a stereolithography (STL) file format for computational modelling. The regions
of the sinonasal cavity included the external face, nostrils, nasal cavity,
each of the paranasal sinuses and the nasopharynx. Fig 1 demonstrates the
surface partitioning of the pre-operative model highlighting the local regions
used to record deposition values. Each image shows outer walls removed
sequentially with the final image depicting the middle turbinate as red:
anterior head of the middle turbinate; green: lateral side; and yellow:medial.
These regions are removed during the virtual middle turbinectomy resection
which was performed in Slicer to create a ‘post-operative’ model.

The segmented models were imported into the commercial CFD software Ansys-
Fluent 2020R1 and meshed with polyhedral cells and six prism layers applied in
near-wall regions. The advantage of poly-hexcore meshing is that it uses fewer
elements, approximate $3.5\times$ fewer than tetrahedral meshing (with the
same size functions) while achieving negligibly different results. The mesh
strategy was based on a size functions evaluated in an earlier study in
[Inthavong et al., 2018] which found that a minimum scoped mesh cell size of
0.45 mm in the nasal cavity region but a much larger, coarse mesh in the outer
ambient air region was optimal. During our mesh independence analysis
(Supplementary Material) we evaluated the change in velocities in coronal
planes resulting in an optimised mesh of 1.58 million cells for the pre-op
model and 1.75 million cells for the post-op model.

Fig 2 shows the internal volume mesh taken in the sagittal planes of the left
and right cavities. The zoomed view of a posterior near wall region depicts
the six prism layers, with polyhedrals connecting to the hex-core cells to
make up the polyhex-core mesh design. The sagittal planes are taken at the
same locations which highlights the opened airway caused by surgery in the
Post-Op models. The difference in the opened airway is most significant in the
right cavity.

###  2.2 Fluid particle simulation

Two constant flow rates were used (5 and 15 L/min), corresponding to shallow
and normal resting breathing rates based on a laminar, incompressible flow
model. The laminar model provides better overall performance for the flowrates
used where the literature [Van Strien et al., 2021, Li et al., 2017] have
demonstrated laminar dominant flow behaviour at flow rates of 15 L/min or less
[Hahn et al., 1993, Doorly et al., 2008].

The flow equations describing the conservation for mass and momentum are
expressed as:

$\displaystyle\frac{\partial u_{i}}{\partial x_{i}}$ $\displaystyle=0$ (1)
$\displaystyle\frac{\partial u_{j}u_{i}}{\partial x_{j}}$
$\displaystyle=-\frac{1}{\rho}\frac{\partial p}{\partial
x_{i}}+\frac{\partial}{\partial x_{j}}\left(\mu\frac{\partial u_{i}}{\partial
x_{j}}\right)$ (2)

where $u_{i}$ is the flow velocity vector and $p$ is the pressure, and $\mu$
is the molecular viscosity. The equations were discretized using the finite
volume approach, with upwind second order schemes for the pressure, and
momentum equations. The pressure-velocity coupling used the ’Coupled’ scheme
and the solution was iterated through the ’pseudo transient’ option and a time
scale factor of 0.01 s. The pseudo-transient method improves stability and
speeds up convergence through a form of implicit relaxation (as opposed to
explicit relaxation). The convergence criteria was set to 1e-5 where the
residual values were calculated based on local scaling.

In this study, only the gravity and the drag force were considered for
tracking the aerodynamic sphere and pollen particles. The particle sizes of
interest were in the micron size range, and therefore forces that affect
submicron particles such as Brownian, and lift forces were excluded as they
are negligible for micron particles. The particle equation of motion is given
by

$\frac{\mathrm{d}u_{p}}{\mathrm{d}t}=F_{d}+\frac{g\left(\rho_{p}-\rho_{g}\right)}{\rho_{p}}$
(3)

where the subscripts $p$ and $g$ represent particle and gas (air) phases,
$\rho$ is the density $g$ is the gravity term and $F_{d}$ is the drag force
per unit mass.

$F_{d}=\frac{18\mu}{\rho_{p}d_{p}^{2}}\frac{C_{d}\mathrm{Re}}{24}(u_{g}-u_{p})=\frac{(u_{g}-u_{p})}{\tau}$
(4)

where $\tau$ is the particle relaxation time, and $d_{p}$ is taken as the
volume equivalent diameter $d_{v}$. For normalising particles with different
densities the aerodynamic diameter is related by

$d_{ae}=d_{g}\sqrt{\rho_{p}/\rho_{0}}$ (5)

where $d_{g}$ is the geometric particle diameter, $\rho_{p}$ is the particle
density, and $\rho_{0}$ is the particle unit density equal to 1 g/cm3. For the
spherical particle, the drag coefficient $C_{d}$ is taken from Morsi and
Alexander [1972] defined by

$C_{d}=a_{1}+\frac{a_{2}}{\mathrm{Re}_{p}}+\frac{a_{3}}{\mathrm{Re}_{p}}$ (6)

where the $a_{1}$, $a_{2}$, and $a_{3}$ are empirical constants for smooth
spherical particles over different ranges of particle Reynolds number.

For pollen, an assumption of a non-spherical particle shape based on the
experimental measurements of Tran-Cong et al. [2004] was used to define the
non-spherical shapes through an aggregation of particles clustered together.
Fig 3 shows scanned electron microscope images of pollen reported in Suárez-
Cervera et al. [2008] and from Dartmouth College. The rough surface non-
spherical particle assumption was used following Tran-Cong et al. [2004]’s
technique (Fig 3c). Two equivalent diameters, volume equivalent diameter
$d_{n}$, and surface equivalent diameter $d_{a}$ and a shape factor called the
‘degree of circularity’ Wadell [1933] were used for the drag correlation given
as,

$d_{n}=\sqrt[3]{6V/\pi};\hskip 56.9055ptd_{a}=\sqrt{4A_{proj}/\pi};\hskip
56.9055ptc=\pi\frac{d_{A}}{P_{proj}}$ (7)

where $V$ is the particle volume, and $A_{proj}$ is the projected area of the
sphere. The circularity $c$ is defined as where $P_{proj}$ is the projected
perimeter of the particle in its direction of motion. The empirically defined
correlation for the drag coefficient for the pollen is given as

$C_{D}=\frac{24}{\textrm{Re}_{p}}\frac{d_{a}}{d_{n}}\left[1+\frac{0.15}{\sqrt{c}}\left(\frac{d_{a}}{d_{n}}\textrm{Re}_{p}\right)^{0.687}\right]+\frac{0.42\left(\frac{d_{a}}{d_{n}}\right)^{2}}{\sqrt{c}\left[1+42500\left(\frac{d_{a}}{d_{n}}\textrm{Re}_{p}\right)^{-1.16}\right]}$
(8)

This correlation was implemented into Ansys-Fluent via a user-defined function
allowing a customised drag correlation and has been used in a previous study
Inthavong et al. [2008]. The proposed correlation is valid in the ranges of
variables $0.15<Re<1500$; $0.80<d_{a}/d_{n}<1.50$; and $0.4<c<1.0$. Particles
were individually tracked through the nasal cavity using a one-way coupled
Lagrangian Discrete Phase Model, where the airflow field was first simulated
and subsequently trajectories of each particle determined by integrating the
particle force balance equation while accounting for gravity and drag. A
particle number independence was evaluated which showed results began
converging with less than 5% change in deposition results from 15,000 to
20,000 particles. To ensure particle number indpendence, approximatly 30,000
particles were introduced for each particle diameter tested (e.g. $n=16$
diameters $\times 30,000$ particles). The initial particle locations to
represent outside air exposure under the influence of inhalation flow rates
was set by distributing the particles uniformly on a spherical dome 3 cm from
the face based on the breathing region distance [Shang et al., 2015, Zhu et
al., 2005]. The particles were set with a zero velocity, and the spherical
dome distance from the nostrils allowed the particles to be naturally inhaled
under the influence of the inhalation flow rates.

The morphology of pollen comes in various shapes, diameters and densities,
which is subjected to different aerodynamic properties. For example, Crowder
et al. [2002] reported geometric diameters ranging from $16-27\upmu$m, Kelly
et al. [2002] found a range between $10-60\upmu$m. A density of 550kg/m3 for
pollen was adopted based on Crowder et al. [2002]. The pollen shape exhibited
spiny protrusions that arise from a spherical surface base (see Fig 3). The
geometric diameters for the spherical and pollen particles used in this study
ranged from 4 to 70 microns (4, 6, 8, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40,
50, 60, 70).

##  3 Results

###  3.1 Flow Field

Velocity contours on ten coronal planes from the anterior, main nasal passage
and the posterior regions are shown in Fig 4. The highest velocities were
found in the anterior planes, which exhibited the smallest cross-sectional
areas. The white to red colours in the contours depict the main bulk flow
paths which were located in the medial regions of the planes as the airflow is
directed superiorly from the nostrils. The effect of surgery is evident in
planes $y4$ to $y8$ where the removal of the turbinate bone creates a large
open space in the superior half of the main nasal passage. The flow becomes
more disturbed and has greater penetration into the most superior locations of
the cavity.

Three sagittal planes were taken in both the left and right cavities (for a
total of six planes), and their velocity contours are shown in Fig 5. The
middle turbinate is located by a yellow surface which disrupts the main flow
paths in the Pre-Op model. In contrast, in the Post-Op model the airflow moves
freely towards the ceiling of the nasal cavity (shown in lateral, and middle
planes). In the septal sides of both left and right cavities, the planes show
the septal wall protruding into the anterior airway, and this is confirmed in
plane $y3$ of the anterior coronal planes from Fig 4. It is expected that
inhaled particles will accelerate with the flow field directed at an angle
shown in the velocity contours in the sagittal planes. Furthermore, the
particles have enhanced opportunity to deposit on anteriorly protruding
anterior septal wall, and on the middle turbinate in the Pre-Op model, and
nasal cavity ceiling on the Post-Op model.

###  3.2 Inhalation Efficiency

The inhaled fraction of particles released from an upstream radial position of
3-cm from the nostrils was determined, shown in Fig 6. The radial distance
selected followed the work from the face was based on the breathing influence
and region [Shang et al., 2015, Zhu et al., 2005]. The inhalation fraction
differs from the inhaled efficiency measure (also known as aspiration
efficiency), where the latter is determined from the ratio of aspirated
particle concentration to the ambient freestream particle concentration
[Kennedy and Hinds, 2002, Dai et al., 2006, Li et al., 2016, Tao et al., 2020,
Dai et al., 2006, King Se et al., 2010].

The inhaled fraction of pollen demonstrated greater opportunities of
inhalation compared to aerodynamic spheres. The spherical particles were
treated as aerodynamic diameters with $\rho_{p}=1000$kg/m3 which is
$1.82\times$ greater than pollen which contributes to the spherical particle
inertia, mass, and its sedimentation away from the inhaled air. For the same
diameters, the inhaled fraction of the lighter pollen particles is
significantly greater than the spheres. When the pollen particle diameters are
normalised by it’s aerodynamic diameter (Fig 6b) it’s inhalation fraction
curve moves closer to the spherical particle, but remains greater for the same
aerodynamic diameter.

A secondary influence on the aerodynamic flight of the pollen and spherical
particles is the difference in drag coefficient shown in Fig 6c for Re number
range of $10^{-1}$ to $10^{3}$. The pollen experiences greater drag force
suggesting it adapts to the inhaled air easier than a spherical particle. The
combined effects of lower density and a higher drag coefficient in the pollen
led to greater chances of inhalation than in an aerodynamic spherical
particle. The multi-parameters that affect the inhalation fraction can be
collapsed into a single term which we define as the ‘inhalation parameter’ as:

$\mathrm{IP}=d_{ae}^{2}Q^{-b}$ (9)

where $b$ was a variable tuned to curve fit the deposition data, which was
$b=0.825$. This combines the inhalation rate $Q$ in units of [cm3/s],
aerodynamic equivalent diameter $d_{ae}$ in units of $\upmu$m to conveniently
capture the influence of inhalation rate and particle density. The data
collapsed onto a single trendline for either particle type allowing curve
fitted correlations to predict the inhalation fraction from exposure to
particles released 3 cm from the nostrils.

###  3.3 Local Surface Deposition

The local deposition fraction for spherical and pollen particles on individual
surfaces are shown on Fig 7. In general, the particle diameters achieving peak
deposition shifted to the right (increased diameter sizes) for the right nasal
cavity compared with the left, suggesting the anatomical differences in the
cavities led to different deposition behaviour. In the left cavity, deposition
was found predominantly on the nasal septum, and the effects of surgery were
generally slightly higher. In the right nasal cavity, deposition was mainly
found on the anterior middle turbinate in Pre-Op.

After removal of the middle turbinates, the particles travelled further into
the nasal airway and deposited on the upper lateral walls. Overall deposition
decreased in the Post-Op model at 5L/min, but was relatively constant at
15L/min. Deposition on the nasal septum was relatively similar between Pre-Op
and Post-Op. There was increased deposition for the higher flow rate of
15L/min with more particles depositing in the anterior regions. The right
cavity of the Post-Op model also exhibited increased deposition in the
paranasal sinuses, which was negligible in the Pre-Op models.

The deposition sites on the nasal cavity are shown in Fig 8 where the
particles are coloured by particle diameter). There is a consistent
concentration of particles depositing on the anterior nasal septum that
protrudes into the airway (circled once in the Fig 8c). At 15L/min a second
local concentration site is found at the vestibule. These two regions account
for the high inertial particles (black to red coloured diameters) impacting on
surfaces directly aligned to the inhalation path. At 5L/min more inhaled
particles penetrate the anterior nasal valve, and when the airflow slows in
the main nasal passage, the particles succumb to gravitational sedimentation
and deposition is found along the nasal floor.

Deposition in the main nasal passage appears to scatter across all surfaces
randomly. Fig 9 extracted deposition on the nasal septum before and after
surgery with the middle turbinates were included in Pre-Op model. The head of
the anterior middle turbinate is the frontline region for contact with the
inhaled air and particles, and it is expected to receive high deposition. The
lateral, and medial middle turbinate regions squeeze into the airspace and
deposition in these regions arise from changes in the flow streamlines (i.e.
changes in flow direction). The Post-Op images show the effect of middle
turbinectomy resection, but there are no consistent changes among all cases.
In general, however, there is an increase in deposition in the superior half
of the nasal septum, attributed to the more open passageway allowing more flow
movement the nasal ceiling.

###  3.4 Deposition Trends

Fig 10 collects all the deposition data into a single ‘swarm’ plot (similar to
a box plot). The values spread laterally for repeating deposition fraction
values; thus, a wide section of the plot suggests a common occurrence of the
data. The height represents the occurrence of high deposition fractions in the
specific surface. The results for the Pre-Op model shows that the highest
deposition occurs in the nasal septum ($>0.1$), and this is only achieved for
inhalation rates of 15L/min (red markers only). The next highest deposition
occurs in the anterior middle turbinate (up to a fraction of 0.3). In the
Anterior, and to a lesser extent the Lateral-Upper regions, the higher flow
rate enhanced their deposition. In the Post-Op model (Fig 10b)the middle
turbinate regions (MT-anterior, MT-lateral, MT-medial) are removed, and there
is an increase in the Lateral-Upper, and Sinuses. At the same time, there is
very low deposition in the Lateral-Lower, and Turbinate (superior, and
inferior) regions.

The change in deposition fraction due to surgery is obtained by subtracting
Pre-Op values from Post-Op values shown in Fig 10c. The results indicate that
middle turbinate resection in this patient generally increased deposition in
the Lateral-Upper, Sinuses, and the Septum. The highest increases in
deposition occurred for either particle inhaled at 15L/min. Increased
deposition between 0.0 to 0.03 in the Septum and Sinuses were produced from a
combination of particle type and inhalation rate. A reduction in the
deposition was found in all surfaces except the sinuses, and this was
predominantly caused by inhalation rates of 5L/min. Fig 10d shows the same
change in deposition fraction due to surgery, but this time the markers are
coloured by nasal cavity side. This highlighted the left nasal cavity as the
main contributor for deposition increases in the Sinuses, and Lateral-Upper.
In contrast, the right nasal cavity demonstrated lower values of deposition
decreases in the Lateral-Upper, Lateral-Lower, and Turbinate surfaces.

##  4 Discussion

Pollen particles exhibited higher drag coefficients, and a lower particle
density compared to aerodynamic equivalent spheres. This suggests that pollen
has greater mobility in its aerodynamic flight and greater potential to
penetrate the nasal cavity. Middle turbinate resection in this study created a
larger open space in the superior half of the main nasal passage where the
inhaled air was directed into, thus transporting inhaled particles to the
space. At 5L/min, smaller particles penetrated deeper into the nasal cavity
and leaving through the nasopharynx exit, leading a decrease in nasal
deposition. At 15L/min, the inertial properties of the particles are enhanced
threefold, and this led to increased anterior deposition, particularly the
protruding nasal septum near the nasal valve.

The middle turbinate directly faces the oncoming flow and bifurcates the air
into lateral and septal streams. This suggests that the propensity for
deposition onto either surface is sensitive to the orientation of the anterior
middle turbinate relative to the oncoming air. The left cavity exhibited very
high deposition on the nasal septum before surgery, whereas the right cavity
exhibited more evenly distributed deposition across different surfaces but
predominantly across the anterior, middle turbinate, and the septum. The
larger open space produced from the surgery increased accessibility to the
superior half of the main nasal passage, and this was reflected by the
increase in deposition in the ‘Lateral-Upper’ region, and this was found only
in the left nasal cavity at 15L/min inhalation flow rate. There was an
increase in the nasal septum, but to a lesser extent, however, the changes may
be masked by the initially high deposition in the anterior nasal septum in the
Pre-Op model.

CFD analysis demonstrates that, in the unoperated normal model, a large number
of particles deposit on the middle turbinate, particularly on the anterior
edge and lateral surface of the middle turbinate, along with the septum.
Following resection of the middle turbinatess, particle deposition increases
at the septum and at the upper lateral nasal walls. These deposition trends
demonstrates that middle turbinate excision leads to increased allergen
deposition on the septum, confirming the clinical suspicion that the middle
turbinate plays a protective role in preventing allergen impaction on the
septum. In patients with aspirin exacerbated respiratory disease (AERD),
surgical strategies aimed at preserving the middle turbinate should be
considered.

##  5 Conclusion

An inhalation parameter was introduced to account for the multi-parameter
scenario of inhaled particles from the surrounding air, given as
$d^{2}_{ae}Q^{-0.825}$. The drag coefficient for pollen and its lighter
density led to greater inhalation efficiency compared to spherical particles
at 1000kg/m3. Correlations to curve fit the pollen, and spherical particle was
produced to predict inhalation fractions. Deposition on local surfaces was
obtained for Pre-Op and Post-Op models. The effect of surgery showed an
increase in deposition in the Post-Op at 15L/min but a decrease at 5L/min.
Middle turbinate resection created a more open space in the superior half of
the main nasal passage, increasing accessibility to the region. The increase
in deposition in the Post-Op model occurred in the left cavity in the Lateral-
Upper area. There was also an increase on the nasal septum, but this value may
be masked by the initially high deposition already found in the Pre-Op model.
It is expected that surgical planning of middle turbinate resection should
consider the implications of increased allergen deposition on the septum which
can lead to secondary sequela. This is particularly critical for patients with
aspirin exacerbated respiratory disease (AERD). Future work will investigate
more population samples/models, and the influence of transient inhalation
breathing flow.

##  6 Acknowledgements

The authors gratefully acknowledge the financial support provided by the
Garnett Passe and Rodney Williams Foundation Conjoint Grant 2019-22.

## References

  * Alam et al. [2019] Alam, S., Li, C., Bradburn, K.H., Zhao, K., Lee, T.S., 2019.  Impact of Middle Turbinectomy on Airflow to the Olfactory Cleft: A Computational Fluid Dynamics Study.  American Journal of Rhinology and Allergy doi:10.1177/1945892418816841.
  * Bagheri and Bonadonna [2016] Bagheri, G., Bonadonna, C., 2016\.  On the drag of freely falling non-spherical particles.  Powder Technology doi:10.1016/j.powtec.2016.06.015.
  * Bahmanzadeh et al. [2015] Bahmanzadeh, H., Abouali, O., Faramarzi, M., Ahmadi, G., 2015\.  Numerical simulation of airflow and micro-particle deposition in human nasal airway pre- and post-virtual sphenoidotomy surgery.  Computers in Biology and Medicine doi:10.1016/j.compbiomed.2015.03.015.
  * Bousquet et al. [2001] Bousquet, J., Van Cauwenberge, P., Khaltaev, N., 2001.  Allergic rhinitis and its impact on asthma.  doi:10.1067/mai.2001.118891.
  * Chen et al. [2012] Chen, X.B., Lee, H.P., Chong, V.F., Wang, D.Y., 2012\.  Drug delivery in the nasal cavity after functional endoscopic sinus surgery: A computational fluid dynamics study.  Journal of Laryngology and Otology doi:10.1017/S0022215112000205.
  * Crowder et al. [2002] Crowder, T.M., Rosati, J.A., Schroeter, J.D., Hickey, A.J., Martonen, T.B., 2002.  Fundamental effects of particle morphology on lung delivery: Predictions of Stokes’ law and the particular relevance to dry powder inhaler formulation and development.  doi:10.1023/A:1014426530935.
  * Dai et al. [2006] Dai, Y.T., Juang, Y.J., ying Wu, Y., Breysse, P.N., Hsu, D.J., 2006.  In vivo measurements of inhalability of ultralarge aerosol particles in calm air by humans.  Journal of Aerosol Science doi:10.1016/j.jaerosci.2005.10.002.
  * Dayal et al. [2016] Dayal, A., Rhee, J.S., Garcia, G.J., 2016.  Impact of middle versus inferior total turbinectomy on nasal aerodynamics, in: Otolaryngology - Head and Neck Surgery (United States).  doi:10.1177/0194599816644915.
  * De Wang et al. [2012] De Wang, Y., Lee, H.P., Gordon, B.R., 2012.  Impacts of fluid dynamics simulation in study of nasal airflow physiology and pathophysiology in realistic human three-dimensional nose models.  doi:10.3342/ceo.2012.5.4.181.
  * DelGaudio et al. [2019] DelGaudio, J.M., Levy, J.M., Wise, S.K., 2019.  Central compartment involvement in aspirin-exacerbated respiratory disease: the role of allergy and previous sinus surgery.  International Forum of Allergy and Rhinology doi:10.1002/alr.22367.
  * Doorly et al. [2008] Doorly, D., Taylor, D., Schroter, R., 2008.  Mechanics of airflow in the human nasal airways.  Respiratory physiology & neurobiology 163, 100–110.
  * Frank et al. [2013a] Frank, D.O., A.M., Z., V.H., D., K.A., M., G.M., F., C.S., E., B.A., S., J.S., K., 2013a.  Quantification of airflow into the maxillary sinuses before and after functional endoscopic sinus surgery.  International Forum of Allergy and Rhinology .
  * Frank et al. [2013b] Frank, D.O., Kimbell, J.S., Cannon, D., Rhee, J.S., 2013b.  Computed intranasal spray penetration: Comparisons before and after nasal surgery.  International Forum of Allergy and Rhinology doi:10.1002/alr.21070.
  * Hahn et al. [1993] Hahn, I., Scherer, P.W., Mozell, M.M., 1993.  Velocity profiles measured for airflow through a large-scale model of the human nasal cavity.  Journal of applied physiology 75, 2273–2287.
  * Hariri et al. [2015] Hariri, B.M., Rhee, J.S., Garcia, G.J., 2015.  Identifying patients who may benefit from inferior turbinate reduction using computer simulations.  Laryngoscope doi:10.1002/lary.25367.
  * Hölzer and Sommerfeld [2008] Hölzer, A., Sommerfeld, M., 2008\.  New simple correlation formula for the drag coefficient of non-spherical particles.  Powder Technology doi:10.1016/j.powtec.2007.08.021.
  * Inthavong et al. [2018] Inthavong, K., Chetty, A., Shang, Y., Tu, J., 2018.  Examining mesh independence for flow dynamics in the human nasal cavity.  Computers in Biology and Medicine doi:10.1016/j.compbiomed.2018.09.010.
  * Inthavong et al. [2019] Inthavong, K., Das, P., Singh, N., Sznitman, J., 2019\.  In silico approaches to respiratory nasal flows: A review.  doi:10.1016/j.jbiomech.2019.109434.
  * Inthavong et al. [2013] Inthavong, K., Mouritz, A.P., Dong, J., Tu, J.Y., 2013.  Inhalation and deposition of carbon and glass composite fibre in the respiratory airway.  Journal of Aerosol Science doi:10.1016/j.jaerosci.2013.07.003.
  * Inthavong et al. [2008] Inthavong, K., Wen, J., Tian, Z., Tu, J., 2008.  Numerical study of fibre deposition in a human nasal cavity.  Journal of Aerosol Science doi:10.1016/j.jaerosci.2007.11.007.
  * Kelly et al. [2002] Kelly, J.K., Rasch, A., Kalisz, S., 2002.  A method to estimate pollen viability from pollen size variation.  American Journal of Botany doi:10.3732/ajb.89.6.1021.
  * Kennedy and Hinds [2002] Kennedy, N.J., Hinds, W.C., 2002\.  Inhalability of large solid particles.  Journal of Aerosol Science doi:10.1016/S0021-8502(01)00168-9.
  * King Se et al. [2010] King Se, C.M., Inthavong, K., Tu, J., 2010.  Inhalability of micron particles through the nose and mouth.  Inhalation Toxicology doi:10.3109/08958370903295204.
  * Lee et al. [2016] Lee, K.B., Jeon, Y.S., Chung, S.K., Kim, S.K., 2016\.  Effects of partial middle turbinectomy with varying resection volume and location on nasal functions and airflow characteristics by CFD.  Computers in Biology and Medicine doi:10.1016/j.compbiomed.2016.08.014.
  * Li et al. [2017] Li, C., Jiang, J., Dong, H., Zhao, K., 2017.  Computational modeling and validation of human nasal airflow under various breathing conditions.  Journal of Biomechanics 64, 59–68.
  * Li et al. [2016] Li, X., Inthavong, K., Tu, J., 2016.  Numerical investigation of micron particle inhalation by standing thermal manikins in horizontal airflows.  Indoor and Built Environment doi:10.1177/1420326X14551069.
  * Maza et al. [2019] Maza, G., Li, C., Krebs, J.P., Otto, B.A., Farag, A.A., Carrau, R.L., Zhao, K., 2019\.  Computational fluid dynamics after endoscopic endonasal skull base surgery—possible empty nose syndrome in the context of middle turbinate resection.  International Forum of Allergy and Rhinology doi:10.1002/alr.22236.
  * Morsi and Alexander [1972] Morsi, S.A., Alexander, A.J., 1972\.  An investigation of particle trajectories in two-phase flow systems.  Journal of Fluid Mechanics doi:10.1017/S0022112072001806.
  * Rhee et al. [2011] Rhee, J.S., Pawar, S.S., Garcia, G.J., Kimbell, J.S., 2011\.  Toward personalized nasal surgery using computational fluid dynamics.  Archives of Facial Plastic Surgery doi:10.1001/archfacial.2011.18.
  * Sanjeevi et al. [2018] Sanjeevi, S.K., Kuipers, J.A., Padding, J.T., 2018.  Drag, lift and torque correlations for non-spherical particles from Stokes limit to high Reynolds numbers.  International Journal of Multiphase Flow doi:10.1016/j.ijmultiphaseflow.2018.05.011.
  * Shang et al. [2015] Shang, Y.D., Inthavong, K., Tu, J.Y., 2015.  Detailed micro-particle deposition patterns in the human nasal cavity influenced by the breathing zone.  Computers and Fluids doi:10.1016/j.compfluid.2015.02.020.
  * Siu et al. [2020a] Siu, J., Dong, J., Inthavong, K., Shang, Y., Douglas, R.G., 2020a.  Quantification of airflow in the sinuses following functional endoscopic sinus surgery.  Rhinology doi:10.4193/Rhin19.387.
  * Siu et al. [2020b] Siu, J., Shrestha, K., Inthavong, K., Shang, Y., Douglas, R., 2020b.  Particle deposition in the paranasal sinuses following endoscopic sinus surgery.  Computers in Biology and Medicine doi:10.1016/j.compbiomed.2019.103573.
  * Suárez-Cervera et al. [2008] Suárez-Cervera, M., Vega-Maray, A., Castells, T., Rodríguez-Rajo, F.J., Asturias, J.A., Le Thomas, A., Seoane-Camba, J.A., 2008.  An approach to the knowledge of pollen and allergen diversity through lipid transfer protein localisation in taxonomically distant pollen grains.  Grana doi:10.1080/00173130802513776.
  * Tao et al. [2020] Tao, Y., Yang, W., Inthavong, K., Tu, J., 2020\.  Indoor particle inhalability of a stationary and moving manikin.  Building and Environment doi:10.1016/j.buildenv.2019.106545.
  * Tian and Ahmadi [2013] Tian, L., Ahmadi, G., 2013.  Fiber transport and deposition in human upper tracheobronchial airways.  Journal of Aerosol Science doi:10.1016/j.jaerosci.2013.02.001.
  * Tian et al. [2012] Tian, L., Ahmadi, G., Wang, Z., Hopke, P.K., 2012\.  Transport and deposition of ellipsoidal fibers in low Reynolds number flows.  Journal of Aerosol Science doi:10.1016/j.jaerosci.2011.09.001.
  * Tian et al. [2007] Tian, Z.F., Inthavong, K., Tu, J.Y., 2007.  Deposition of inhaled wood dust in the nasal cavity.  Inhalation Toxicology doi:10.1080/08958370701665525.
  * Tran-Cong et al. [2004] Tran-Cong, S., Gay, M., Michaelides, E.E., 2004.  Drag coefficients of irregularly shaped particles.  Powder Technology doi:10.1016/j.powtec.2003.10.002.
  * Van Strien et al. [2021] Van Strien, J., Shrestha, K., Gabriel, S., Lappas, P., Fletcher, D.F., Singh, N., Inthavong, K., 2021.  Pressure distribution and flow dynamics in a nasal airway using a scale resolving simulation.  Physics of Fluids 33, 011907\.
  * Wadell [1933] Wadell, H., 1933.  Sphericity and roundness of rock particles 41, 310–331.  doi:10.1086/624040.
  * Wang et al. [2018] Wang, Y., Zhou, L., Wu, Y., Yang, Q., 2018.  New simple correlation formula for the drag coefficient of calcareous sand particles of highly irregular shape.  Powder Technology doi:10.1016/j.powtec.2017.12.004.
  * Xiong et al. [2008] Xiong, G., Zhan, J., Zuo, K., Li, J., Rong, L., Xu, G., 2008.  Numerical flow simulation in the post-endoscopic sinus surgery nasal cavity.  Medical and Biological Engineering and Computing doi:10.1007/s11517-008-0384-1.
  * Zastawny et al. [2012] Zastawny, M., Mallouppas, G., Zhao, F., van Wachem, B., 2012\.  Derivation of drag and lift force and torque coefficients for non-spherical particles in flows.  International Journal of Multiphase Flow doi:10.1016/j.ijmultiphaseflow.2011.09.004.
  * Zhang et al. [2019] Zhang, Y., Shang, Y., Inthavong, K., Tong, Z., Sun, B., Zhu, K., Yu, A., Zheng, G., 2019.  Computational investigation of dust mite allergens in a realistic human nasal cavity.  Inhalation Toxicology doi:10.1080/08958378.2019.1647315.
  * Zhao et al. [2014] Zhao, K., Malhotra, P., Rosen, D., Dalton, P., Pribitkin, E.A., 2014.  Computational Fluid Dynamics as Surgical Planning Tool: A Pilot Study on Middle Turbinate Resection.  Anatomical Record doi:10.1002/ar.23033.
  * Zhu et al. [2005] Zhu, S., Kato, S., Murakami, S., Hayashi, T., 2005\.  Study on inhalation region by means of CFD analysis and experiment.  Building and Environment doi:10.1016/j.buildenv.2004.11.009.

Figure 1: Computational model depicting (a) the surface mesh, and the initial
particle positions (labelled in red) used for inhalation exposure analysis;
(b) the face, anterior region (orange) which consists of the nasal vesitbule,
and nasal valve, paranasal sinuses (green), nasopharynx (blue); (c) the face
removed; (d) the main nasal passage with individual sub-regions; followed by
(e) a zoomed view of the middle turbinates with the outer lateral walls
removed.

Figure 2: Internal volume mesh shown on sagittal planes in the left and right
cavities for the (A) Pre-Op and (B) Post-Op models. A zoomed vieew of the mesh
near the wall depicts the prism layers which are connected to the hex-core
mesh elements through polyhedral cells. (C) The computational domain is shown
with the imposed boundary conditions. An extension was added at the
nasopharynx to prevent backflow occurring at the outlet.

Figure 3: Scanning electron microscope images of pollen grains reproduced
from: (a) Louisa Howard and Charles Daghlian of Dartmouth College – . Zeiss
DSM 962 SEM, Rippel Electron Microscope Facility; and (b) Suárez-Cervera et
al. [2008] demonstrating pollen morphological diversity .

Figure 4: Ten coronal slices (viewed from the face) taken in the anterior,
main nasal passage, and the posterior regions. The anterior and posterior
slices were identical in the Pre-Op and Post-Op models, but differed in the
main nasal passage, where the middle turbinate resection took place. Velocity
contours for the anterior planes (y1 to y3) were identical and therefore only
one set is shown. The main passage (y3 to y8) and posterior planes (y9, y10)
demonstrated different velocity flow fields. The peak velocity was 4.1m/s in
the anterior coronal planes. The contour legend was capped at 3m/s to provide
better contrast in colours in the main and posterior planes.

(a) Sagittal plane locations

(b) Sagittal planes in the right nasal cavity

(c) Sagittal planes in the left nasal cavity

Figure 5: (a) Top view, looking down at the nasal cavity showing sagittal
planes in the left and right nasal cavities. (b) Velocity contours in the
right nasal cavity, and; (c) velocity contours in the left nasal airway for
Pre-Op and Post-Op models.

(a) Inhalation fraction for sphere and pollen

(b) Inhalation fraction normalised by aerodynamic diameter

(c) Drag coefficient for a sphere and pollen particle

(d) Curve fit of inhalation fraction for sphere, pollen.

Figure 6: (a) Inhalation fraction for sphere and pollen particles at 5L/min
and 15L/min inhalation rates. (b) Inhalation fraction for sphere and pollen
normalised by aerodynamic diameter. (c) Drag coefficient values for pollen
(Eqn 6) and sphere (Eqn 6). (d) Inhalation fraction as a function of an
inhalation parameter for pollen and spherical particles released at 3-cm
radius from nostril inlets. The curve fit correlation provides a method for
normalising the inhaled fraction
$IF=0.85\exp(-0.001053x^{2})+0.15\exp(-0.02x)$

Figure 7: Deposition fraction, defined as a ratio of the number of particles
depositing on a surface to the number introduced from the breathing region
dome. Therefore, the surface deposition excludes particles escaping through
the nasopharynx or not inhaled. (a) Spherical particles inhaled at 5L/min. (b)
Pollen particles inhaled at 5L/min. (c) Spherical particles inhaled at
15L/min. (d) Pollen particles inhaled at 15L/min. Each subfigure depicts the
deposition for PreOp (top panels), and PostOp (lower panels) models, separated
in to left and right cavities.

Figure 8: Deposition sites of inhaled particles. The paranasal sinuses have
been excluded for improved visualisation of the main nasal passage deposition
which was dominant. (a) Spherical particles inhaled at 5L/min. (b) Pollen
particles inhaled at 5L/min. (c) Spherical particles inhaled at 15L/min. (d)
Pollen particles inhaled at 15L/min. Each subfigure depicts the deposition for
Pre-Op (upper locations), and Post-Op (lower locations) models, separated in
to left and right cavities

.

Figure 9: Deposition on the septum and middle turbinate (surface is coloured
in yellow). Local deposition sites are coloured by surface: black - septum;
blue - lateral middle turbinate; red - anterior middle turbinate; green -
medial middle turbinate. (a) Spherical particles inhaled at 5L/min. (b) Pollen
particles inhaled at 5L/min. (c) Spherical particles inhaled at 15L/min. (d)
Pollen particles inhaled at 15L/min. Each subfigure depicts the deposition for
Pre-Op (upper locations), and Post-Op (lower locations) models, separated in
to left and right cavities.

(a) Deposition efficiency in Pre-Op model

(b) Deposition efficiency in Post-Op model

(c) Change in deposition (Post-Op minus Pre-Op) coloured by particle type

(d) Change in deposition (Post-Op minus Pre-Op) coloured by nasal cavity side

Figure 10: Swarm plots for total deposition fractions in one cavity. Points on
the plot having the same values of deposition will spread laterally.The height
correlates with the deposition fraction, thus a higher point represents an
occurrence of greater deposition fraction.

## Supplementary Material 1: Mesh Independence Analysis

Mesh independence analysis demonstrating the variation in velocity magnitude
for a coarse mesh of 0.98 million polyhex-core cells, and a more refined mesh
of 1.58 million polyhex-core cells.


Execution Time Program Verification With Tight Bounds

Ana Carolina Silva1 Manuel Barbosa1,2 Mário Florido1,3

A. C. Silva et al.

FCUP, Universidade do Porto INESC TEC LIACC, Universidade do Porto

This paper presents a proof system for reasoning about execution time bounds for a core imperative programming language. Proof systems are defined for three different scenarios: approximations of the worst-case execution time, exact time reasoning, and less pessimistic execution time estimation using amortized analysis. We define a Hoare logic for the three cases and prove its soundness with respect to an annotated cost-aware operational semantics. Finally, we define a verification conditions (VC) generator that generates the goals needed to prove program correctness, cost, and termination. Those goals are then sent to the Easycrypt toolset for validation. The practicality of the proof system is demonstrated with an implementation in OCaml of the different modules needed to apply it to example programs. Our case studies are motivated by real-time and cryptographic software.

§ INTRODUCTION

§ BACKGROUND

§ RELATED WORK

§ COST-AWARE PROGRAM LOGIC

§ AMORTIZED COSTS

§ EXACT COSTS

§ IMPLEMENTATION

§ CONCLUSIONS

[1]
Albert, E., Arenas, P., Genaim, S., Gómez-Zamalloa, M., Puebla, G., Ramírez,
D., Román, G., Zanardini, D.: Termination and cost analysis with costa and
its user interfaces. Electronic Notes in Theoretical Computer Science
258(1),  109–121 (2009).
proceedings of the Ninth Spanish Conference on Programming and Languages
(PROLE 2009)

[2]
Almeida, J., Frade, M., Pinto, J., Sousa, S.: Rigorous Software Development. An
Introduction to Program Verification. Springer Verlag (01 2011).

[3]
Atkey, R.: Amortised resource analysis with separation logic. Logical Methods
in Computer Science  7(2) (jun 2011).

[4]
Avanzini, M., Dal Lago, U.: Automating sized-type inference for complexity
analysis. Proc. ACM Program. Lang.  1(ICFP) (aug 2017).
10.1145/3110287, <https://doi.org/10.1145/3110287>

[5]
Barbosa, M., Barthe, G., Grégoire, B., Koutsos, A., Strub, P.Y.: Mechanized
proofs of adversarial complexity and application to universal composability.
In: Proceedings of the 2021 ACM SIGSAC Conference on Computer and
Communications Security. p. 2541–2563. CCS '21, Association for Computing
Machinery, New York, NY, USA (2021). 10.1145/3460120.3484548,

[6]
Brockschmidt, M., Emmes, F., Falke, S., Fuhs, C., Giesl, J.: Alternating
runtime and size complexity analysis of integer programs. In:
Ábrahám, E., Havelund, K. (eds.) Tools and Algorithms for the
Construction and Analysis of Systems. pp. 140–155. Springer Berlin
Heidelberg, Berlin, Heidelberg (2014)

[7]
Carbonneaux, Q., Hoffmann, J., Ramananandro, T., Shao, Z.: End-to-end
verification of stack-space bounds for c programs. SIGPLAN Not.
49(6),  270–281 (jun 2014). 10.1145/2666356.2594301,

[8]
Carbonneaux, Q., Hoffmann, J., Shao, Z.: Compositional certified resource
bounds. SIGPLAN Not.  50(6),  467–478 (jun 2015).
10.1145/2813885.2737955, <https://doi.org/10.1145/2813885.2737955>

[9]
Cook, S.A.: Soundness and completeness of an axiom system for program
verification. SIAM Journal on Computing  7(1),  70–90 (1978)

[10]
Dijkstra, E.W.: Usability engineering. Englewood Cliffs, N.J. : Prentice-Hall

[11]
Floyd, R.W.: Assigning meanings to programs. In: Proceedings of Symposia in
Applied Mathematics Vol. 19. pp. 19–32. Amer. Math. Soc. (1967)

[12]
Gulwani, S., Mehra, K.K., Chilimbi, T.: Speed: Precise and efficient static
estimation of program computational complexity. In: Proceedings of the 36th
Annual ACM SIGPLAN-SIGACT Symposium on Principles of Programming Languages.
p. 127–139. POPL '09, Association for Computing Machinery, New York, NY,
USA (2009). 10.1145/1480881.1480898,

[13]
Haslbeck, M.P.L., Nipkow, T.: Hoare logics for time bounds. In: Beyer, D.,
Huisman, M. (eds.) Tools and Algorithms for the Construction and Analysis of
Systems. pp. 155–171. Springer International Publishing, Cham (2018)

[14]
Hoare, C.A.R.: An axiomatic basis for computer programming. Commun. ACM
12(10),  576–580 (oct 1969). 10.1145/363235.363259,

[15]
Hoffmann, J., Das, A., Weng, S.C.: Towards automatic resource bound analysis
for ocaml. SIGPLAN Not.  52(1),  359–373 (jan 2017).
10.1145/3093333.3009842, <https://doi.org/10.1145/3093333.3009842>

[16]
Hoffmann, J., Hofmann, M.: Amortized resource analysis with polynomial
potential. In: Gordon, A.D. (ed.) Programming Languages and Systems. pp.
287–306. Springer Berlin Heidelberg, Berlin, Heidelberg (2010)

[17]
Hoffmann, M., Jost, S.: Type-based amortised heap-space analysis. In: Sestoft,
P. (ed.) Programming Languages and Systems. pp. 22–37. Springer Berlin
Heidelberg, Berlin, Heidelberg (2006)

[18]
Kaminski, B.L., Katoen, J.P., Matheja, C., Olmedo, F.: Weakest precondition
reasoning for expected runtimes of randomized algorithms. J. ACM
65(5) (aug 2018). 10.1145/3208102,

[19]
King, J.C.: A program verifier. Carnegie Mellon University (1970)

[20]
Leroy, X.: Formal verification of a realistic compiler. Commun. ACM
52(7),  107–115 (jul 2009). 10.1145/1538788.1538814,

[21]
Manna, Z.: The correctness of programs. Journal of Computer and System Sciences
3(2),  119–127 (1969)

[22]
Nielson, H.R.: A hoare-like proof system for analysing the computation time of
programs. Science of Computer Programming  9(2),  107–136 (1987).

[23]
Nielson, H., Nielson, F.: Semantics with Applications: An Appetizer.
Undergraduate Topics in Computer Science, Springer London (2007),

[24]
Radiček, I., Barthe, G., Gaboardi, M., Garg, D., Zuleger, F.: Monadic
refinements for relational cost analysis. Proc. ACM Program. Lang.
2(POPL) (dec 2017). 10.1145/3158124,

[25]
Serrano, A., López-García, P., Hermenegildo, M.V.: Resource usage analysis of
logic programs via abstract interpretation using sized types. Theory and
Practice of Logic Programming  14 (05 2014).

[26]
Simões, H., Vasconcelos, P., Florido, M., Jost, S., Hammond, K.: Automatic
amortised analysis of dynamic memory allocation for lazy functional programs.
SIGPLAN Not.  47(9),  165–176 (sep 2012).
10.1145/2398856.2364575, <https://doi.org/10.1145/2398856.2364575>

[27]
Tarjan, R.E.: Amortized computational complexity. SIAM Journal on Algebraic
Discrete Methods  6(2),  306–318 (1985)

[28]
Vasconcelos, P.B., Jost, S., Florido, M., Hammond, K.: Type-based allocation
analysis for co-recursion in lazy functional languages. In: Programming
Languages and Systems - 24th European Symposium on Programming, ESOP 2015,
Proceedings. LNCS, vol. 9032, pp. 787–811. Springer (2015)