# Non-Markovian Quantum Control via Model Maximum Likelihood Estimation and
Reinforcement learning

Tanmay Neema
Computer Science Laboratory
SRI International
Menlo Park, CA 94025
<EMAIL_ADDRESS>
&Susmit Jha
Computer Science Laboratory
SRI International
Menlo Park, CA 94025
<EMAIL_ADDRESS>
&Tuhin Sahai
Applied Sciences
SRI International
Menlo Park, CA 94025
<EMAIL_ADDRESS>

###### Abstract

Reinforcement Learning (RL) techniques have been increasingly applied in
optimizing control systems. However, their application in quantum systems is
hampered by the challenge of performing closed-loop control due to the
difficulty in measuring these systems. This often leads to reliance on assumed
models, introducing model bias, a problem that is exacerbated in open quantum
dynamics where Markovian approximations are not valid. To address these
challenges, we propose a novel approach that incorporates the non-Markovian
nature of the environment into a low-dimensional effective reservoir. By
initially employing a series of measurements as a ’dataset’, we utilize
machine learning techniques to learn the effective quantum dynamics more
efficiently than traditional tomographic methods. Our methodology aims to
demonstrates that by integrating reinforcement learning with model learning,
it is possible to devise control policies and models that can counteract
decoherence in a spin-boson system. This approach may not only mitigates the
issues of model bias but also provides a more accurate representation of
quantum dynamics, paving the way for more effective quantum control
strategies.

##  1 Introduction

Quantum computation is an exciting area with the potential to revolutionize
material and drug discovery [1], cryptography [2], optimization and search
[3], and efficient simulation of high-dimensional dynamical systems [4, 5]. By
exploiting quantum mechanical phenomena such as superposition and
entanglement, researchers have constructed novel algorithms that provide
exponential to quadratic speed-up over state-of-the-art approaches. For
example, quantum algorithms provide exponential gains for the task of
computing prime factors for any given integer [6]. However, all these gains
require the development of fault tolerant quantum platforms that are robust to
external noise.

Quantum control and state initialization remain central problems towards the
development of fault tolerant and robust quantum computing as according to
DiVincenzo’s criteria [7]. Generally, current methods of quantum control rely
on either assumed models and their accessible gradients and/or rely on
feedback through measurements [8, 9, 10, 11, 12, 13] . For example, Krotov’s
algorithm minimizes a functional $J$ that represents the update cost of the
control parameters $\\{\epsilon^{i}_{l}(t)\\}$ at the i-th iteration and the
infidelity $J_{T}$ of the final propagated state and target state [14, 10]. In
another example of the Gradient Pulse Engineering algorithm, the fidelity is
maximized across the time discretized control parameters through gradient
ascent [15]. Notably, these methods assume model knowledge in finding
analytical gradient. As a result, optimal control theory may not be robust to
open quantum dynamics where additional dissipative and decohering terms remain
unmodeled. To illustrate this point, many of these models work utilize Master
equations at the limit of the Markovian approximation which assume that the
coupling to the environment is weak or the inital state is uncorrelated with
the environment [16, 17], which for the purpose of state initialization is an
unreasonable assumption. Consequently, these control policies fail to produce
desirable final state fidelities when transfering from simulation to
experiment.

In regards to feedback-based control schemes, the process of measurement of a
quantum system extends duration of control process dependent on interaction
times with the measurement device, often acting as computational bottlenecks.
Strasberg et al. demonstrated that time of measurements scale with the
dimensionality of the system [18]. In addition, the process of measurement,
both strong and weak measurements measurements, necessarily collapse the
system into an eigenstate for a projection-valued measures (PVM) or an unknown
state for positive operator valued measures (POVM), where the probability of a
measurement depends on the trace of the state and the measurement operator
[19]). As a result, feedback dependent methods have undefined protocol
lengths, dependent on the inherent stochasticity of a quantum measurement such
as in a sliding mode controller described in the work byDong et al. [13]. Our
proposed method attempts to mitigate both of these issues, through a model
learning phase and policy learning phase.

Reinforcement learning has been demonstrated to be successful in quantum
control applications [20, 21, 22, 23, 24, 25]. Zhang et al. compared the
performance between Krotov’s, stochastic gradient descent (SGD), tabular Q
learning (TQL), deep Q learning (DQL), and policy gradient(PG) and found that
DQL and PG performed best [21]. These RL methods led to more simplified
control sequences, leading to shorter operation times and higher fidelities.
Niu et al. demonstrated that trusted-region policy optimization via
reinforcement learning was robust to leakage and stochastic gate errors [20].
Ma et al. employed curriculum-based reinforcement learning and demonstrated
that it was superior to gradient and genetic algorithms, even in the increase
in the dimensionality of the system and was capable of perform the protocol in
fewer steps [22]. The robustness of the reinforcement learning in the presence
of noise, its effectiveness in larger systems, and the protocal speed ups
motivate the usage of RL in our proposed method.

Still these methods suffer from model-bias, which can lend to poor application
to experimental setups. As a result, as a part of the model learning phase,
our proposed method aims to learn the open quantum dynamics of the system from
measurements.

For a general dynamical map, the quantum process tomography grows
exponentially with the system size and relies on enough measurements to find
an accurate reconstruction, a problem further exacerbated by the fact that the
evolution superoperators do not form a semigroup. In general, the dynamics are
described by a memory kernel, $\mathcal{K}$, in the Nakajima Zwanzig equation
[16].
$\frac{d}{dt}\rho_{S}(t)=-iL_{s}\rho_{S}(t)+\int_{0}^{t}d\tau\,\mathcal{K}(t-\tau)\rho_{S}(\tau)$

Finding the exact form of the kernel can be accomplished through determining a
series of transfer tensors [26, 27, 28, 29]. The evolution is described by,

$\rho(t_{m})=\begin{pmatrix}T_{1}&T_{2}&\cdots&T_{K}\end{pmatrix}\begin{pmatrix}\rho(t_{m-1})\\\
\rho(t_{m-2})\\\ \vdots\\\ \rho(t_{m-K})\end{pmatrix},$ (1)

where ${T_{i}}$ can be found as a function of the dynamical maps between any
time $t_{0}$ to $t_{i}$ and is truncated at $T_{K}$ as the long-term
correlations approach 0. This too, however in practice, encounters the same
difficulty of quantum tomography methods, i.e., having to learn the series of
dynamical maps. In an attempt to reduce sampling size Torlai et al. recently
utilized unsupervised learning [30], indicating that machine learning has
utility in generalizing quantum dynamics.

Recent work by Luchnikov et al. [31] demonstrated that with a sufficiently
long sequence of measurements of a quantum system, non-Markovian effects, even
in strongly coupled regimes, may be learned, demonstrating a sharp distinction
between traditional tomographic approaches where a series of identical states
are prepared and evolved. Here, by evolving the system with intermediate and
equally spaced projective measurements, the model attempts to maximize the
probability a sequence of measurements $p(\\{E_{i}\\}|\mathcal{L}_{S+ER})$
over an iteratively learned generator $\mathcal{L}_{S+ER}$ that describes the
evolution of the system and environment embedding. Luchnikov et al. found that
the method was able to generate maps with a given error $\epsilon$ in
$\epsilon^{-2}$ samples independent of the specific map. In traditional
tomographic methods, one also needs $\epsilon^{-2}$ samples for given error
$\epsilon$, they require $K$ times more measurements to produce $K$ maps with
the same error (e.g. $\epsilon\propto\sqrt{\frac{K}{n}}$).

Treating the model learning as overhead, the proposed method simulates the
quantum system approximately and utilizes reinforcement learning agent as the
controller. Notably, we demonstrate that that the suggested Markovian
embedding can still be useful in highly decohered system such as in a relaxed
spin boson, i.e., in highly correlated systems and environments. In addition,
we demonstrate that a reinforcement learning agent can actively bring the
system away from thermal equilibrium through controls that apply directly to
the system.

The following paper is structured as follows, section 2 describes the
simulated model of interest and the dynamics of the two state system. Section
3 details the learning of the open dynamics and details the derivation of the
iterative learning of the quantum map that maximizes the probability of a
sequence of the simulated data. Section 4 introduces the actor-critic policy
gradient reinforcement learning algorithm and the choice of reward function.
Finally, Section 5 concludes by describing the results and discusses future
implications.

##  2 Model of Open Quantum System

For two state systems, such as in qubits, spin-boson model has been able to
effectively capture the environment and system interactions. The model
consists of a two state system and its linear interaction with a spectrum of
harmonic oscillators. Inspired by this model, we simulate a system with the
following Hamiltonian,

$H_{\text{total}}=H_{0}+H_{c}(t)+H_{\text{bath}}+H_{\text{int}},$ (2)

where each individual Hamiltonian term are given by:

$\displaystyle H_{0}$ $\displaystyle=-\frac{g}{2}B_{0}\sigma_{z},$ (3)
$\displaystyle H_{c}(t)$ $\displaystyle=-\frac{g}{2}B_{x}(t)\sigma_{x},$ (4)
$\displaystyle H_{\text{bath}}$
$\displaystyle=\sum_{k}\hbar\omega_{k}\hat{a}_{k}^{\dagger}\hat{a}_{k},$ (5)
$\displaystyle H_{\text{int}}$
$\displaystyle=\sigma_{z}\sum_{k}\lambda_{k}(\hat{a}_{k}^{\dagger}+\hat{a}_{k}).$
(6)

Here, $H_{0}$ represents the static part of the Hamiltonian, $H_{c}(t)$ the
control Hamiltonian which may be time-dependent, $H_{\text{bath}}$ the
Hamiltonian for the bath of oscillators (which models the environment) with a
discrete range of frequencies that is characterized by the spectral density
$J(\omega)=\sum_{k}|\lambda_{k}|^{2}\delta(\omega-\omega_{k})$, and
$H_{\text{int}}$ the interaction Hamiltonian between the qubit and the bath.

For simplicity, we assume that the environmental noise contains a single mode,
i.e., the qubit is coupled to a single quantum harmonic oscillator (QHO), with
$J(\omega)=\lambda_{0}^{2}\delta(\omega-\omega_{0})$ and
$\omega_{0}\approx\frac{gB_{0}}{2\hbar}$ so that the QHO frequency is near
resonant with the qubit system, ensuring stronger non-Markovian interactions.

During model learning, we initialize the simulated bath state from thermal
equilibrium at a temperature $T$ such that the energy of the oscillator can be
described by the equipartion theorem i.e. $k_{B}T\gg\hbar\omega_{0}$ (we set
$k_{B}=\hbar=1$). So for the simulated bipartite qubit environment system,
this would be equivalent to
$\rho(0)=\rho_{S}(0)\otimes\rho_{bath}(0)=\ket{e}\bra{e}\otimes\frac{1}{Z}e^{\frac{-H_{bath}}{T}}$.

The control policy being able to navigate the control landscape of the open
quantum system utilizing interactions exclusively with the system (which would
be unable to otherwise due to unitary nature of these transformations) would
indicate a deep understanding of the dynamics gleaned from the model and would
indicate an ability to reverse relaxation and dissipative activity.

##  3 Learning Open Quantum Dynamics

In this section, we describe the model learning algorithm which follows the
description provided by Luchnikov et. al. [31]. We describe an equivalent but
different derivation for greater clarity.

The dynamics of the open quantum system is broken into two parts, an effective
reservoir (ER) and the system (S) and the method attempts to learn the
dynamics of the reservoir interactions through a series of measurements
performed on the system of interest. A dataset of projective value
measurements $\\{E_{i}\\}$ at each time steps $t_{i}$, i.e., the system of
interest (in our case simulation) is evolved under open quantum dynamics for
time $\Delta t$ and then is followed by a measurement within an orthonormal
basis $\\{\ket{\psi_{k}}\bra{\psi_{k}}\\}$. $E_{i}$ is the measurement from
the set $\\{\ket{\psi_{k}}\bra{\psi_{k}}\\}$ associated with an outcome $k$ at
time step $t_{i}$ with a probability
$p_{i}=\text{tr}(\ket{\psi_{k}}\bra{\psi_{k}}\rho_{S}(t_{i})\ket{\psi_{k}}\bra{\psi_{k}}$.
The total probability of the sequence of $\\{E_{i}\\}$,
$p=\Pi_{i=1}^{n}p_{i}$. As $n\rightarrow\infty$, an expectation would be that
probability of said series of measurements would be maximized under the real
system dynamics. Motivated by this intuition, we aim to learn a generator
$\mathcal{L}_{S+ER}$ that maximizes the probability
$p(\\{E_{i}\\}|\mathcal{L})$,

$\mathcal{L}^{*}=\underset{\mathcal{L}}{\arg\max}\,p.$

Evidently, we have no means to directly measure the environment, hence, as a
hyperparameter, we vary the estimated dimensionality of the system $d_{ER}$.
We consider the evolution of a $\rho_{S+ER}$ which once traced out yields the
approximate system dynamics of $\rho_{S}$ given by,

$\dot{\rho}_{S+ER}=\mathcal{L}\rho_{S+ER}.$ (7)

The channel representing the evolution between each time step is given as
$\Phi$,

$\Phi=e^{\mathcal{L}_{S+ER}\Delta t}.$ (8)

Given the channel $\Phi$, we now express the probability of
${E_{1},E_{2},E_{3}...}$ spaced at times $\Delta t$ as:

$p=\text{tr}\left(E_{n}\Phi\left[E_{n-1}\Phi\left[\ldots\left(E_{1}\Phi\left[\rho_{0}\right]E_{1}\right)\ldots\right]E_{n-1}\right]E_{n}\right)$
(9)

An equivalent expression of $p$ can be given by the forward and
backpropagating operators $\tilde{\rho_{i}}$ and $\mathcal{E}_{i}$
respectively,

$p=\text{tr}\left(\tilde{\rho_{i}}\mathcal{E}_{i}\right)\quad\forall
i\in\\{0,1,2,\ldots,n\\},$ (10)

where:

$\displaystyle\tilde{\rho_{i}}$
$\displaystyle=E_{i}\Phi[\tilde{\rho}_{i-1}]E_{i},\quad\text{with}\quad\tilde{\rho_{0}}=\rho_{0},$
(11) $\displaystyle\mathcal{E}_{i}$
$\displaystyle=\Phi^{\dagger}\left[E_{i}\mathcal{E}_{i+1}E_{i}\right],\quad\text{with}\quad
E_{n}=I.$ (12)

Naïvely, one could initialize $\Phi$ with its ${(d_{S}d_{ER})}^{4}$ parameters
and apply an optimization algorithm to learn them, however, since
$\Phi:B(\mathcal{H}_{S+ER})\rightarrow B(\mathcal{H}_{S+ER})$ ($\mathcal{B}$
being the set of bounded linear operators) satisfies the properties of a valid
quantum channel, the set of parameters is constrained.

A valid quantum channel is completely positive and trace-preserving (CPTP).
According to the Stinespring Dilation Theorem, for a CPTP map $\Phi$, there
exists a representation:

$\Phi[\rho_{S+ER}]=\text{tr}_{A}(U(\rho_{S+ER}\otimes\rho_{A})U^{\dagger})$
(13)

with $U=e^{-iH\Delta t}$, $H\in B(\mathcal{H}_{S+ER}\otimes\mathcal{H}_{A})$,
and $\rho_{a}$ is a fixed ancilla system.

The statement of the theorem implies a quantum channel can be represented by a
unitary evolution of a lifted quantum system. Ideally, the dimension of
$d_{A}=d_{S}d_{ER}$ so that the set of learnable parameters still matches the
dimensionality of the channel $(d_{S}d_{ER})^{4}$.

Finding $\mathcal{L}_{S+ER}$ now comes down to learning the joint Hamiltonian
of the system, effective reservoir, and ancilla, given by:

$H^{*}=\underset{H}{\arg\max}\,p$ (14)

To determine this $H^{*}$, we need to find an analytical gradient of the
probability with respect to the Hamiltonian matrix entries. The derivation is
as follows,

$\displaystyle\frac{d\log p}{dH_{\mu\nu}}$
$\displaystyle=\frac{1}{p}\text{tr}\left(\frac{d\tilde{\rho}_{n}}{dH_{\mu\nu}}\mathcal{E}_{n}\right),$
$\displaystyle=\frac{1}{p}\text{tr}\left(\frac{d\Phi}{dH_{\mu\nu}}[\tilde{\rho}_{n-1}]E_{n}\mathcal{E}_{n}E_{n}+\frac{d\tilde{\rho}_{n-1}}{dH_{\mu\nu}}\mathcal{E}_{n-1}\right),$
$\displaystyle=\frac{1}{p}\sum_{i=0}^{n-1}\text{tr}\left(\frac{d\Phi}{dH_{\mu\nu}}[\tilde{\rho}_{i}]E_{i}\mathcal{E}_{i}E_{i}\right),$
$\displaystyle=\frac{1}{p}\sum_{i=0}^{n-1}\text{tr}\left(\left(\frac{dU}{dH_{\mu\nu}}\left[\tilde{\rho}_{i}\otimes\rho_{A}\right]U^{\dagger}+U\left[\tilde{\rho}_{i}\otimes\rho_{A}\right]\frac{dU^{\dagger}}{dH_{\mu\nu}}\right)\left(E_{i}\mathcal{E}_{i}E_{i}\right)\otimes
I_{A}\right).$

Since $p$ is function of a series of linear operators, taking the derivative
can be expanded easily by the product rule, yielding the form shown. We
unravel the expression of the $p$ backwards from $n$ to 0 to achieve the
result. The rest follows trivially.

Next, we find the derivative of the unitary with respect to approximated
Hamiltonian. We perform a first order expansion in $H_{\mu\nu}$ with a
paramter $\epsilon$ at $\epsilon=0$ which is equivalent to taking the
derivative with the matrix element.

$\displaystyle\frac{dU}{dH_{\mu\nu}}$
$\displaystyle=\frac{d}{d\epsilon}U(H+\epsilon\ket{\mu}\bra{\nu})\big{|}_{\epsilon=0}$
$\displaystyle=\frac{d}{d\epsilon}e^{-i(H+\epsilon\ket{\mu}\bra{\nu})\Delta
t}$ $\displaystyle=\frac{d}{d\epsilon}\sum^{\infty}_{n=0}\frac{(-i\Delta
t)^{n}}{n!}(H+\epsilon\ket{\mu}\bra{\nu})^{n}\big{|}_{\epsilon=0}$
$\displaystyle=\frac{d}{d\epsilon}\left(I+-i\Delta
tH+\sum^{\infty}_{n=1}\frac{(-i\Delta
t)^{n}\epsilon}{n!}\sum^{n-1}_{j=0}H^{n-1-j}\ket{\mu}\bra{\nu}H^{j}+\mathcal{O}(\epsilon^{2})\right)\big{|}_{\epsilon=0}$
$\displaystyle=\sum^{\infty}_{n=1}\frac{(-i\Delta
t)^{n}}{n!}\sum^{n-1}_{j=0}H^{n-1-j}\ket{\mu}\bra{\nu}H^{j}$

Since $H$ is hermitian, it can be spectrally decomposed into:

$H=\sum\lambda_{k}P_{k},$

where $P_{k}P_{l}=\delta_{kl}P_{l}$ and are projectors for the eigenstates,
i.e.,

$H^{n}=\sum\lambda_{k}^{n}P_{k}.$

Hence,

$\displaystyle\frac{dU}{dH_{\mu\nu}}$
$\displaystyle=\sum_{k,l}\sum^{\infty}_{n=1}\frac{(-i\Delta
t)^{n}}{n!}\sum^{n-1}_{j=0}\lambda_{k}^{n-1-j}\lambda_{l}^{j}P_{k}\ket{\mu}\bra{\nu}P_{l},$
$\displaystyle=\sum_{k,l}\sum^{\infty}_{n=1}\frac{(-i\Delta
t)^{n}}{n!}\frac{\lambda_{k}^{n-1}\left(1-\left(\frac{\lambda_{l}}{\lambda_{k}}\right)^{n}\right)}{1-\left(\frac{\lambda_{l}}{\lambda_{k}}\right)}P_{k}\ket{\mu}\bra{\nu}P_{l},$
$\displaystyle=\sum_{k,l}\sum^{\infty}_{n=1}\frac{(-i\Delta
t)^{n}}{n!}\frac{\lambda_{k}^{n}-\lambda_{l}^{n}}{\lambda_{k}-\lambda_{l}}P_{k}\ket{\mu}\bra{\nu}P_{l},$
$\displaystyle=\sum_{k,l}\frac{\left(e^{-i\lambda_{k}\Delta
t}-1\right)-\left(e^{-i\lambda_{l}\Delta
t}-1\right)}{\lambda_{k}-\lambda_{l}}P_{k}\ket{\mu}\bra{\nu}P_{l},$
$\displaystyle=\sum_{k,l}\frac{e^{-i\lambda_{k}\Delta
t}-e^{-i\lambda_{l}\Delta
t}}{\lambda_{k}-\lambda_{l}}\braket{\psi_{k}}{\mu}\braket{\nu}{\psi_{l}}\ket{\psi_{k}}\bra{\psi_{l}}.$

Note that for $k=l$, the coefficient reduces to $-i\Delta te^{-i\Delta
t\lambda_{l}}$.

Now that we have the form of the gradient of the total logarithmic
probability, we use optimization methods to compute the maximum value of
$\text{log}(p)$. Specifically, we use the ADAM optimizer [32] over the matrix
entries. The ADAM optimizer adaptively estimates first and second order
moments and has low memory requirements and is appropriate since $p$ is highly
non-convex.

We finish the training when the geometric mean of $\\{p_{i}\\}$,
$p^{\frac{1}{n}}$ stabilizes or when the algorithm has exceeded a prescribed
number of training steps.

Once we have computed $H^{*}$, we can easily determine the channel $\Phi$ and
generator $\mathcal{L^{*}}$. For computational ease, we treat the channel and
generators as matrices over an operator basis in a Liouville space. These
matrices admit the following property,

$vec(\Phi\rho)=M(\Phi)vec(\rho).$ (15)

Consequently, this leads to the representation of $\mathcal{L}_{S+ER}$, which
can be computed by diagonalization,

$M(\mathcal{L}_{S+ER})=\frac{1}{\Delta t}\ln(M(\Phi)).$ (16)

##  4 Actor Critic RL Controller

Reinforcement learning is challenging in when the state space is continuous
and the actor-critic method has emerged as a promising approach in these
settings. Combining the REINFORCE algorithm [33] of policy gradient and
temporal difference learning, the method iteratively updates the policy and
the estimated state value function. The approach tends to be sample efficient,
while mitigating the high bias of value-based methods and high variance of
exclusively policy-based approaches [34, 33].

We initialize an artificial neural network as our policy $\pi$ and state value
approximation $V$, parameterized by $\theta$ and $w$ respectively. The steps
in our algorithm are outlined below,

At every time step, we measure the temporal difference (TD) error (without a
discount) as,

$\delta_{t}=r_{t}+V(s_{t+1})-V(s_{t}),$

where a positive $\delta_{t}$ implies that the state $s_{t}$ is in fact more
valuable than the estimate. In other words, the transition reward being
greater than approximated on that instant. Consequently, we update the $w$
parameters in an attempt to minimize the sum of $\delta_{t}^{2}$ through
bootstrapping,

$w:=w+\alpha_{w}\delta_{t}\nabla_{w}V(S).$

Similarly, we would want to minimize action probabilities that would produce
these TD errors with a good measure of loss given by
$-log(\pi(a_{t}|s_{t}))\delta_{t}$, thus,

$\theta:=\theta-\alpha_{\theta}\delta_{t}\nabla_{\theta}log(\pi(a_{t}|s_{t})).$

Concretely speaking, the actor update is inspired by the Policy Gradient
Theorem with a baseline, but in contrast due to having access to an estimate
of the value function, the algorithm can estimate
$R_{t}-b(s_{t})\approx\delta_{t}$, with $R_{t}$ being the cumulative reward of
the subsequent states. This means that the actor-critic method is online
(independent of a future state), has less variance, and converges faster, as
updates are not performed with respect to that expectation estimate at the
episode termination as in the REINFORCE algorithm.

In our control problem, we define the reward at time $t$, $r_{t}$, as follows,

$r_{t}=F(\text{tr}_{ER}(\rho_{S+ER}(t)),\rho^{*})\mathcal{D}(\rho_{S+ER}(t))-1,$

where $F$ is the fidelity of the state $\text{tr}_{ER}(\rho_{S+ER},\rho^{*})$
given by,

$F(\rho_{S},\rho_{*})=tr\sqrt{\sqrt{\rho^{*}}\rho_{S}\sqrt{\rho^{*}}}.$

We define a metric $\mathcal{D}$ as,

$\mathcal{D}=F(\rho_{S+ER},\text{tr}_{ER}(\rho_{S+ER})\otimes\text{tr}_{S}(\rho_{S+ER})).$
(17)

Here the metric $\mathcal{D}$ captures the measure of separability of the
environment and system. Consequently, its value is equal to 1 in the case of
factorable state. Multiplying the fidelity and ‘separability’, the system is
guided to the desired final state while also being uncorrelated with the
environment.

This reward structure provides non-sparse learning signals by ensuring that a
reward (necessarily negative) is given at each time step based on the product
of $F(\rho_{S+ER},\rho^{*})$ and $\mathcal{D}(\rho_{S+ER})$ minus 1, which
reaches a maximum of 0, when it satisfies the fidelity of 1 and ’separability’
of 1. The episode is terminated at time $T$ or when $r_{t}>\alpha$. The
termination criteria and necessarily negative rewards encourage the policy to
achieve an adequate level of fidelity as quickly as possible, thus, promoting
efficiency in learning. The agent learning can be thought as being divided as
an exploration phase where the agent learns the system dynamics and the
control effects, with the episode terminating after time $T$. This is followed
by the exploitation phase where the agent acts greedily (achieving higher
reward under shorter periods of time) with the episode terminating after
reaching the latter condition.

##  5 Results and Discussion

Figure 1: Short term comparison of the real dynamics and embedded dynamics.
The blue and red graphs correspond to the model and real dynamics
respectively.

Figure 2: Log probabilities of the model (blue) is graphed against the epoch
number. The system converges early, but later gradually approaches the real
system probabilities, reflecting the nonconvex nature of the optimization
problem.

Figure 3: The control policy of the qubit from excited state to ground state
under weak coupling with QHO environment demonstrates ability to navigate
control landscape in 50 control steps.

The model learning phase demonstrated rapid convergence during the initial
epochs, with significant improvements up to epoch 50. Despite a decrease in
the rate of convergence, the model continued to enhance its accuracy
gradually. This trend underscores the model’s capacity to adapt and refine its
predictions over time, even in the face of the complex, nonconvex optimization
landscape inherent in quantum dynamics modeling.

The decision to truncate the Fock state number of the Quantum Harmonic
Oscillator (QHO) to $N=5$ was driven by computational efficiency
considerations. This truncation proved to be most effective when the embedded
dimension $d_{ER}$ was initialized at 5, balancing the model’s fidelity
against computational constraints. This balance was crucial for managing
fluctuations in the Bloch vector, especially for smaller values of $d_{ER}$,
which tended to average out these fluctuations.

For short durations, the model exhibited high efficacy, closely mirroring the
real dynamics. However, for time scales significantly exceeding the duration
of the dataset ($\Delta t\times 10^{5}$ measurements), discrepancies between
the model and actual dynamics became more pronounced. These deviations suggest
that for long-term predictions, a more sophisticated approach to modeling
relaxation dynamics might be warranted, one as the embedded dimension $d_{ER}$
increases. While further parameterized model could potentially offer greater
accuracy, the trade-off in terms of training speed and computational demand
necessitates careful consideration.

In the context of policy learning, the agent demonstrated proficiency in
navigating states under conditions of low system-environment coupling.
However, with an increase in the coupling constant, the purity of the state
was compromised, indicating challenges in mitigating decoherence and reversing
relaxation through non-unital open quantum system dynamics. This suggests that
while the current approach is effective to a degree, exploring strategies such
as model predictive control might offer a more robust solution for actively
managing decoherence, albeit at the risk of overfitting to the learned model.

The initial success of integrating these methodologies lays a promising
foundation for future research. It paves the way for enhancing the accuracy
and scalability of models dealing with larger systems and more intricate bath
dynamics. The insights gained from this study not only contribute to the
theoretical understanding of quantum systems but also hold potential for
practical applications in quantum computing, where managing decoherence is a
critical challenge.

In conclusion, this study has demonstrated the feasibility and potential of
using truncated models and policy gradient methods for simulating and
controlling quantum systems. The ability to accurately predict and influence
quantum system dynamics, even within computational constraints, opens new
avenues for advancing quantum technology applications.

## References

  * Cao et al. [2018] Yudong Cao, Jhonathan Romero, and Alán Aspuru-Guzik.  Potential of quantum computing for drug discovery.  _IBM Journal of Research and Development_ , 62(6):6–1, 2018.
  * Rieffel and Polak [2000] Eleanor Rieffel and Wolfgang Polak.  An introduction to quantum computing for non-physicists.  _ACM Computing Surveys (CSUR)_ , 32(3):300–335, 2000.
  * Marsh and Wang [2020] Samuel Marsh and Jingbo B Wang.  Combinatorial optimization via highly efficient quantum walks.  _Physical Review Research_ , 2(2):023302, 2020.
  * Gnanasekaran et al. [2023] Abeynaya Gnanasekaran, Amit Surana, and Tuhin Sahai.  Efficient quantum algorithms for nonlinear stochastic dynamical systems.  In _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , volume 02, pages 66–75, 2023.  doi: 10.1109/QCE57702.2023.10186.
  * Surana et al. [2022] Amit Surana, Abeynaya Gnanasekaran, and Tuhin Sahai.  Carleman linearization based efficient quantum algorithm for higher order polynomial differential equations.  _arXiv preprint arXiv:2212.10775_ , 2022.
  * Shor [1999] Peter W Shor.  Polynomial-time algorithms for prime factorization and discrete logarithms on a quantum computer.  _SIAM review_ , 41(2):303–332, 1999.
  * Georgescu [2020] Iulia Georgescu.  The divincenzo criteria 20 years on.  _Nature Reviews Physics_ , 2(12):666–666, 2020.
  * Wu et al. [2018] Re-Bing Wu, Bing Chu, David H Owens, and Herschel Rabitz.  Data-driven gradient algorithm for high-precision quantum control.  _Physical Review A_ , 97(4):042122, 2018.
  * Sivak et al. [2022] VV Sivak, A Eickbusch, H Liu, B Royer, I Tsioutsios, and MH Devoret.  Model-free quantum control with reinforcement learning.  _Physical Review X_ , 12(1):011059, 2022.
  * Morzhin and Pechen [2018] Oleg Morzhin and Alexander Pechen.  Krotov method for optimal control in closed quantum systems.  9 2018.  doi: 10.1070/RM9835.
  * Dong and Petersen [2009a] Daoyi Dong and Ian R Petersen.  Quantum control theory and applications: A survey.  10 2009a.  doi: 10.1049/iet-cta.2009.0508.
  * Chen et al. [2021] Binbin Chen, Jing Wang, and Yunsen Zhou.  Quantum control and its application: A brief introduction.  volume 1802. IOP Publishing Ltd, 3 2021.  doi: 10.1088/1742-6596/1802/2/022068.
  * Dong and Petersen [2009b] Daoyi Dong and Ian R Petersen.  Sliding mode control of quantum systems.  _New Journal of Physics_ , 11:105033, 10 2009b.  ISSN 1367-2630.  doi: 10.1088/1367-2630/11/10/105033.  URL https://iopscience.iop.org/article/10.1088/1367-2630/11/10/105033.
  * Yu and Luo [2021] Zara Yu and Da-Wei Luo.  Many-body quantum state control in the presence of environmental noise.  _arXiv preprint arXiv:2112.06330_ , 2021.
  * de Fouquieres et al. [2011] P. de Fouquieres, S. G. Schirmer, S. J. Glaser, and Ilya Kuprov.  Second order gradient ascent pulse engineering.  2 2011.  doi: 10.1016/j.jmr.2011.07.023.
  * Lidar [2019] Daniel A. Lidar.  Lecture notes on the theory of open quantum systems.  2 2019.
  * Breuer and Petruccione [2007] Heinz-Peter Breuer and Francesco Petruccione.  _The Theory of Open Quantum Systems_.  Oxford University PressOxford, 1 2007.  ISBN 0199213909.  doi: 10.1093/acprof:oso/9780199213900.001.0001.
  * Strasberg et al. [2020] Philipp Strasberg, Kavan Modi, and Michalis Skotiniotis.  How long does it take to implement a projective measurement?  11 2020.  doi: 10.1088/1361-6404/ac5a7a.
  * Beneduci [2020] Roberto Beneduci.  Notes on naimark’s dilation theorem.  _Journal of Physics: Conference Series_ , 1638:012006, 10 2020.  ISSN 1742-6588.  doi: 10.1088/1742-6596/1638/1/012006.
  * Niu et al. [2019] Murphy Yuezhen Niu, Sergio Boixo, Vadim N. Smelyanskiy, and Hartmut Neven.  Universal quantum control through deep reinforcement learning.  _npj Quantum Information_ , 5, 12 2019.  ISSN 20566387.  doi: 10.1038/s41534-019-0141-3.
  * Zhang et al. [2019] Xiao Ming Zhang, Zezhu Wei, Raza Asad, Xu Chen Yang, and Xin Wang.  When does reinforcement learning stand out in quantum control? a comparative study on state preparation.  _npj Quantum Information_ , 5, 12 2019.  ISSN 20566387.  doi: 10.1038/s41534-019-0201-8.
  * Ma et al. [2020] Hailan Ma, Daoyi Dong, Steven X. Ding, and Chunlin Chen.  Curriculum-based deep reinforcement learning for quantum control.  12 2020.  URL http://arxiv.org/abs/2012.15427.
  * He et al. [2021] Run-Hong He, Rui Wang, Shen-Shuang Nie, Jing Wu, Jia-Hui Zhang, and Zhao-Ming Wang.  Deep reinforcement learning for universal quantum state preparation via dynamic pulse control.  _EPJ Quantum Technology_ , 8:29, 12 2021.  ISSN 2662-4400.  doi: 10.1140/epjqt/s40507-021-00119-6.
  * Wei et al. [2019] Peng Wei, Na Li, and Zairong Xi.  Open quantum system control based on reinforcement learning.  pages 6911–6916. IEEE, 7 2019.  ISBN 978-9-8815-6397-2.  doi: 10.23919/ChiCC.2019.8865335.
  * Bukov et al. [2017] Marin Bukov, Alexandre G. R. Day, Dries Sels, Phillip Weinberg, Anatoli Polkovnikov, and Pankaj Mehta.  Reinforcement learning in different phases of quantum control.  5 2017.  doi: 10.1103/PhysRevX.8.031086.  URL http://arxiv.org/abs/1705.00565http://dx.doi.org/10.1103/PhysRevX.8.031086.
  * Cerrillo and Cao [2013] Javier Cerrillo and Jianshu Cao.  Non-markovian dynamical maps: Numerical processing of open quantum trajectories.  7 2013.  doi: 10.1103/PhysRevLett.112.110401.
  * Gherardini et al. [2021] Stefano Gherardini, Andrea Smirne, Susana F. Huelga, and Filippo Caruso.  Transfer-tensor description of memory effects in open-system dynamics and multi-time statistics.  1 2021.  doi: 10.1088/2058-9565/ac4422.
  * Pollock and Modi [2017] Felix A. Pollock and Kavan Modi.  Tomographically reconstructed master equations for any open quantum dynamics.  4 2017.  doi: 10.22331/q-2018-07-11-76.
  * Chen et al. [2019] Yu-Qin Chen, Kai-Li Ma, Yi-Cong Zheng, Jonathan Allcock, Shengyu Zhang, and Chang-Yu Hsieh.  Non-markovian noise characterization with the transfer tensor method.  5 2019.
  * Torlai et al. [2023] Giacomo Torlai, Christopher J. Wood, Atithi Acharya, Giuseppe Carleo, Juan Carrasquilla, and Leandro Aolita.  Quantum process tomography with unsupervised learning and tensor networks.  _Nature Communications_ , 14:2858, 5 2023.  ISSN 2041-1723.  doi: 10.1038/s41467-023-38332-9.
  * Luchnikov et al. [2019] I. A. Luchnikov, S. V. Vintskevich, D. A. Grigoriev, and S. N. Filippov.  Machine learning non-markovian quantum dynamics.  2 2019.  doi: 10.1103/PhysRevLett.124.140502.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Sutton and Barto [2018] Richard S. Sutton and Andrew G. Barto.  _Reinforcement Learning: An Introduction_.  The MIT Press, second edition, 2018.  URL http://incompleteideas.net/book/the-book-2nd.html.
  * Wang and Zou [2022] Yue Wang and Shaofeng Zou.  Policy gradient method for robust reinforcement learning.  5 2022.


capbtabboxtable[][] 5pt 0.84(0.08,0.93) ©This work has been submitted to the
IEEE for possible publication. Copyright may be transferred without notice,
after which this version may no longer be accessible.

# Can a wearable device detect airway obstruction?

Kofi Odame, Graham Atkins, Maria Nyamukuru, Katherine Fearon

###### Abstract

Objective: Lung health monitoring may enable early detection and treatment of
exacerbations in respiratory diseases such as chronic obstructive pulmonary
disease (COPD). Our objective is to explore the feasibility of using a
wearable device to continuously monitor lung health. Towards that goal, this
work studies the relationship between wearable device measurements of tidal
breathing and spirometer lung function parameters.

Methods: Data was collected from 25 single visit adult volunteers with a
confirmed or suspected diagnosis of COPD. A respiratory chest belt was used to
measure the fractional inspiratory time, respiratory rate and tidal volume of
subjects during quiet breathing. The subjects also performed standard
spirometry in a pulmonary function testing laboratory under the supervision of
trained clinical staff to produce forced expiratory volume in one second
(fev1), forced vital capacity (fvc) and the ratio fev1/fvc. Two classification
models were built and trained: one to detect the presence of lung airway
obstruction, and another to stratify the severity of airway obstruction.

Results: The classifier detected airway obstruction with a sensitivity of
$95\%$ and a specificity of $80\%$. Severity of airway obstruction was
classified as either mild/moderate or severe/very severe with a sensitivity of
$91\%$ and a specificity of $78\%$.

Conclusion: Tidal breathing parameters that are measured with a wearable
device can be used to detect airway obstruction with a level of accuracy that
is comparable to that of conventional screening tests used in primary care
settings. In addition, these parameters can reliably distinguish between
severe and mild airway obstruction.

###### Index Terms:

Wearable sensors, biomedical telemetry, chronic obstructive pulmonary disease,
COPD

##  I Introduction

The 251 million people worldwide who suffer from chronic obstructive pulmonary
disease (COPD) are at risk of experiencing an acute exacerbation of symptoms,
which reduces quality of life, increases mortality, and accounts for the
majority of COPD-related healthcare costs [20].

Frequent assessment of lung function via spirometry could help detect COPD
exacerbations early and mitigate their impact. Patients will experience a drop
in forced expiratory volume in 1 second (fev1) up to 2 weeks before an
exacerbation [22], but unfortunately, there is no way to measure lung function
reliably in the home environment, as home spirometry is hindered by low
patient adherence and poor patient technique in executing the spirometric
maneuvers [4].

Figure 1: A respiratory chest belt is worn just below the sternum and measures
the rising and falling of the chest as the subject breathes in and out. The
respiratory signal (top right inset) contains information about lung health
that may be clinically relevant for COPD management.

COPD management is further complicated by improper diagnosis of the disease.
The strict diagnosis of COPD requires compatible symptoms plus spirometry that
demonstrates a ratio of forced expiratory volume in 1 second to forced vital
capacity (fev1/fvc) that is less than $70\%$. But COPD is typically first
diagnosed in the primary care office, where spirometry can be a challenge
because of a lack of equipment and appropriately-skilled personnel. A lack of
spirometry leads to widespread underdiagnosis and overdiagnosis of the
disease, which result in increased use of healthcare services, without
properly treating the patients’ true condition [5, 8].

A solution to these problems of patient non-adherence, improper spirometer
technique and access to skilled personnel is automatic lung function
monitoring. This would require neither active engagement from the patient nor
supervision from skilled clinical staff. To that end, this paper studies the
relationship between wearable device measurements of tidal breathing and
clinically-measured spirometer parameters.

##  II Wearable Lung Monitoring

A respiratory chest belt is an instrumented strap of fabric that is worn just
below the breast bone (Fig. 1). It has an embedded sensor (e.g. an
accelerometer, or a resistive or inductive stretch sensor) that measures the
rising and falling motion of the chest as the wearer breathes in and out.

The chest belt measures the basic tidal breathing pattern at rest, from which
we can extract parameters like expiratory phase duration, respiratory rate,
and tidal volume. These measurements require no active user engagement, are
decoupled from the user’s skill or motivation level and therefore promote user
adherence. Despite the passive nature of these measurements, they are still
useful for detecting changes in lung health status, as we describe below.

Figure 2: A typical respiratory signal. The inspiratory phase lasts $T_{\rm
i}$ seconds, and the full respiratory cycle (inspiratory and expiratory
phases) lasts $T_{\rm tot}$ seconds. The amplitude of the signal is measured
from trough to peak as $RA$ Newtons.

###  II-A Expiratory phase duration

COPD has two contributory pathologies, chronic bronchitis and emphysema, and
most patients with COPD have both pathologies. Chronic bronchitis causes
inflammation and excess mucus production in the lungs’ airways, resulting in
impeded airflow. When breathing in (that is, during the “inspiratory phase” of
respiration), COPD patients can overcome this airway obstruction by increasing
the work of the diaphragm and intercostal muscles and breathing in more
forcefully. Airway diameter also increases slightly in inspiration as pleural
pressure becomes more negative than the atmosphere. Breathing out (during the
“expiratory phase”) is a more passive process than breathing in, relying on
elastic recoil rather than muscular action during normal respiration. So, with
no active compensatory mechanism to counter it, airway obstruction will slow
down the rate at which patients can breathe out and thus prolong the
expiratory phase [14].

Small airways in the emphysematous lung can also collapse in exhalation, a
major driver of airflow obstruction in patients with COPD, leading to slow
expiratory phase. Some patients with emphysema will utilize pursed lip
breathing to facilitate maximal exhalation, which also contributes to a long
expiratory phase, but ultimately reduces small airway collapse in expiration
[2].

During spirometry, the ratio of forced expiratory volume in 1 second to forced
vital capacity (fev1/fvc) measures how fast one can forcefully empty out the
lungs after taking a deep breath. A properly-executed forced exhalation
maneuver will make maximal use of the intercostal and abdominal muscles, thus
actively working against any airflow resistance due to airway obstruction.
However, there is a finite amount of force that a person can exert during
exhalation. So, the more obstructed the airways, the slower the lungs can be
emptied and the smaller the fev1/fvc and fev1; this echoes the prolonged
expiratory phase observed in tidal breathing.

###  II-B Tidal volume

At the end of the expiratory phase of tidal breathing, there is an equilibrium
between the elastic recoil pressure of the chest wall and that of the lungs.
This pressure equilibrium will always leave the lungs partially inflated with
some residual amount of air, or functional residual capacity (FRC). But COPD
patients experience a loss of lung elastic recoil, and this causes an increase
in FRC [6, 7]. As such, the amount of additional air that COPD patients can
inhale during tidal breathing (known as “tidal volume”) is decreased [11].

The spirometric analog to tidal volume is forced vital capacity (fvc), which
is the volume of air that can be forcefully exhaled from the lungs after
taking a deep breath. Just as with tidal breathing, the expiratory phase of
forced breathing ends with a pressure equilibrium that leaves some amount of
air, the residual volume (RV), remaining in the lungs. The loss of elastic
recoil means that RV is typically larger than normal in COPD. So, like tidal
volume, fvc is decreased in COPD.

###  II-C Respiratory rate

In addition to diminished elasticity, the air sacs of COPD patients suffer
wall damage, which reduces their surface area and decreases diffusion of
oxygen into, and carbon dioxide out of, the bloodstream. COPD patients
typically overcome this gas exchange inefficiency by breathing faster to
increase their lung ventilation [1, 13]. So, worsening COPD is often marked by
an increased resting respiratory rate.

Gas exchange efficiency is measured during pulmonary function testing as
diffusing capacity of the lung for carbon monoxide (DLCO). But DLCO is not a
routine component of COPD management, and spirometric volume parameters like
fev1 and fvc provide no information about the relative concentrations of
oxygen or carbon dioxide. Despite the lack of direct correspondence, we can
expect some correlation between respiratory rate and the volume parameters, as
they all change with increasing COPD severity.

##  III Methods

To answer the question of whether a wearable device can detect airway
obstruction, we collected and analyzed data from a number of patients at the
pulmonary clinic of the Dartmouth Hitchcock Medical Center (Lebanon, NH). We
obtained informed consent from each subject before their inclusion in the
study. All aspects of the study protocol were reviewed and approved by the
Dartmouth College Institutional Review Board (Committee for the Protection of
Human Subjects-Dartmouth; Protocol Number: 00028641; Date of Continuing Review
Approval: 1/29/2019).

###  III-A Subjects

Twenty-five adults, 10 men and 15 women, with either suspected or previously-
diagnosed COPD, and scheduled to receive a physician-ordered spirometry test,
were recruited for this study. A summary of the subjects’ spirometry test
results and their anthropomorphic data is listed in Table I.

| mean $\pm$ std
---|---
Age (y) | $67.6\pm 11.6$
Height (cm) | $165.3\pm 10.7$
Weight (kg) | $73.6\pm 17.1$
BMI (kg/m2) | $27.0\pm 6.08$
fev1/fvc | $0.52\pm 0.16$
fev1 (L) | $1.50\pm 0.81$
fvc (L) | $2.68\pm 1.09$
TABLE I: Summary of anthropomorphic data and spirometry test results of the
subjects recruited for this study (n=$25$, 10 men).

###  III-B Procedure

Prior to standard-of-care pulmonary function testing, each subject donned a Go
Direct Respiration Belt (Vernier Software & Technology, Beaverton, OR) and
breathed quietly while seated. The Respiration Belt (see Fig. 1) measured the
force produced by the chest as the subject breathed in and out. Figure 2 is a
typical example of the resulting tidal breathing signal waveform that we
collected from each subject. After producing 15-20 breaths (approximately 1
minute of data collection), the subject removed the Respiration Belt before
performing standard spirometry testing.

Each subject performed standard spirometry with a PC-based spirometer under
the supervision of trained clinical staff. The subject performed three
spirometery maneuvers that were acceptable and reproducible according to the
American Thoracic Society criteria. We collected the subject’s best values of
fev1/fvc, fev1 and fvc for subsequent data analysis.

###  III-C Parameter extraction

For each subject’s tidal breathing signal waveform, we identified a region of
at least 6 consecutive respiratory cycles with minimal motion artifact and
baseline wander. From this region of clean data, we calculated an average
value each for fractional inspiratory time, respiratory rate and amplitude.

Fractional inspiratory time (fit) is the ratio of the inspiratory time to the
total respiratory cycle time. Referring to Fig. 2, we calculated this as fit
$=T_{\rm i}/T_{\rm tot}$. Note that the quantity $1-$fit is a measure of the
expiratory phase duration, normalized to the respiratory period.

The respiratory rate is the average number of breaths taken per minute,
calculated as $\textsc{rr}=60/({\rm BMI}\cdot T_{\rm tot})$. We normalize the
respiratory rate with body-mass index (BMI), because they are significantly
correlated [3, 12].

The respiratory amplitude, $RA$, is the difference in force between an
adjacent trough-peak pair (see Fig. 2). This quantity provides an
approximation of the tidal volume as follows. The respiratory amplitude is
proportional to the change in radius of the chest as it expands during the
inspiratory phase. Since BMI is proportional to chest area [9], we can
estimate the tidal volume (modulo constants of proportionality) as
tv$=RA\cdot{\rm BMI}$.

###  III-D Data analysis

A correlation matrix was generated for the relation between the tidal
breathing parameters (fractional inspiratory time, respiratory rate, and
estimated tidal volume) and the spirometric variables (fev1/fvc, fev1, fvc).

We used the results of the correlation analysis to identify which of the tidal
breathing measures would be useful as a predictor variable in estimating each
spirometric parameter. The tidal breathing measures with the highest
correlation coefficients were used to predict fev1/fvc, fev1 and fvc via
multiple regression.

We built a k-nearest neighbour (KNN) classifier to detect the absence or
presence of lung airway obstruction, using fit, rr and tv as the predictor
variables, and applying the GOLD criteria on fev1/fvc to obtain the true class
labels. We chose a KNN classifier because the ‘normal’ class data points are
clustered together in fit-rr-tv space (see Fig. 6). The classifier was
evaluated using leave-one-out cross validation as follows. We held out data
from each of the 25 subjects in turn, built a KNN model on the remaining data,
and then evaluated the model on the held-out data.

The tidal breathing parameters were also used to classify the severity of
airway obstruction. Since this is an ordinal classification problem, we solved
it with a classifier that is based on a regression model (independent
variables: fit, rr, tv; dependent variable: percentage predicted fev1). The
regression model’s estimated fev1 score was used to classify airway
obstruction severity, according to the GOLD criteria [19]. We evaluated the
classifier using leave-one-out cross validation on two tasks: (1) coarse
classification, to distinguish between mild/moderate and severe/very severe
obstruction; (2) fine classification, to distinguish between all four severity
levels of mild, moderate, severe and very severe obstruction.

##  IV Results

Table II is a matrix of correlations between tidal breathing parameters and
spirometric variables. fev1 and fev1/fvc are most strongly correlated with fit
($R^{2}=0.295$, $p=0.005$ and $R^{2}=0.274$, $p=0.007$ respectively), while
fvc is most strongly correlated with tidal volume ($R^{2}=0.329$, $p=0.003$).

| fev1/fvc | fev1 | fvc
---|---|---|---
fit | 0.274 | 0.295 | 0.135
| $p=0.007$ | $p=0.005$ | $p=0.070$
rr | 0.129 | 0.08 | 0.0125
| $p=0.078$ | $p=0.171$ | $p=0.595$
tv | 0.060 | 0.317 | 0.329
| $p=0.238$ | $p=0.003$ | $p=0.003$
TABLE II: Matrix of correlations between tidal breathing parameters and
spirometric variables. TABLE III: fev1/fvc regression model.
fev1/fvc$=0.094+1.57\cdot$fit$-0.227\cdot$rr
---
R2 | p-value | RMSE
0.435 | 0.002 | 0.124
Figure 3: fev1/fvc regression model (Table III) predictions versus true
values. The dashed gray curve is the line of equality. TABLE IV: fev1
regression model. fev1 $=-1.16+5.35\cdot$fit$+0.005\cdot$tv
---
R2 | p-value | RMSE
0.427 | 0.002 | 0.640
Figure 4: fev1 regression model (Table IV) predictions versus true values. The
dashed gray curve is the line of equality. TABLE V: fvc linear regression
model. fvc$=1.55+0.0096\cdot$tv
---
R2 | p-value | RMSE
0.329 | 0.003 | 0.909
Figure 5: fvc regression model (Table V) predictions versus true values. The
dashed gray curve is the line of equality.

Tables III, IV and V show summary information for the fev1/fvc, fev1 and fvc
regression models. Figures 3, 4 and 5 are scatter plots of true spirometer
parameter values versus their corresponding regression model estimates. The
coefficients of determination are fev1$:R^{2}=0.427$; fvc$:R^{2}=0.329$;
fev1/fvc$:R^{2}=0.435$. All three models have a p-value of $p<0.005$.

Figure 6: Detection of airway obstruction: scatter plot of class labels versus
tidal breathing parameter features. Black markers indicate normal airways
(fev1/fvc $>70\%$), while red markers indicate airway obstruction. The
asterisks are the true class labels, and the circles are the classifier’s
inferred output Figure 7: Detection of airway obstruction: classifier
confusion matrix.

The performance of the airway obstruction detection model is summarized in the
confusion matrix of Fig. 7. As Table VI shows, the model performs with a
sensitivity of $95\%$ and a specificity of $80\%$.

Sensitivity | Specificity | Accuracy | F1-Score
---|---|---|---
95% | 80% | 87.5% | 95%
TABLE VI: Performance of classifier for detecting the presence/absence of
airway obstruction.

The performance of the severity classification model is summarized in Table
VII (fine classification task) and Table VIII (coarse classification task).
Note that only data from the subjects with airway obstruction (indicated by
red asterisk markers in Fig. 6) were included in this experiment.

| Mild | Mod. | Severe | V. severe
---|---|---|---|---
Precision | - | 75.0% | 60.0% | 0%
Recall | 0% | 75.0% | 85.7% | 0%
TABLE VII: Precision and recall performance of classifier for 4-class
stratification of airway obstruction severity. The classifier attempts to
distinguish between mild, moderate, severe and very severe airway obstruction.
Cohen’s kappa coefficient is $\kappa=0.38$.

For the fine classification task, a Cohen’s kappa coefficient of $\kappa=0.38$
was achieved, indicating fair agreement between the estimated severity and
ground truth. The model performed considerably better on the coarse
classification task (Table VIII), yielding a sensitivity of $91\%$ and a
specificity of $78\%$. The Cohen’s kappa coefficient for this task was
$\kappa=0.7$, indicating substantial agreement between the estimated severity
and ground truth.

Figure 8: Severity of airway obstruction: scatter plot of class labels versus tidal breathing parameter features. Black markers indicate mild or moderate airway obstruction, while red markers indicate severe or very severe obstruction. The asterisks are the true class labels, and the circles are the classifier’s inferred output. Note: only data from the 20 subjects with airway obstruction (red asterisk markers in Fig. 6) were included in this experiment. Figure 9: Severity of airway obstruction: classifier confusion matrix. Sensitivity | Specificity | Accuracy | F1-Score
---|---|---|---
90.9% | 77.8% | 84.3% | 87.0%
TABLE VIII: Performance of classifier for 2-class stratification of airway
obstruction severity. The classifier attempts to distinguish between
mild/moderate and severe/very severe airway obstruction. The severe/very
severe class is considered the positive case. Cohen’s kappa coefficient,
$\kappa=0.7$.

##  V Discussion

The relatively strong correlation between fev1/fvc (and fev1) and fit in Table
II is to be expected, given that airway obstruction is a common underlying
factor for both parameters. Similarly, fvc showed the strongest correlation
with tidal volume, due to the effect of lung elastic recoil on both these
measures.

Our fev1/fvc model’s accuracy in detecting airway obstruction is comparable to
that of conventional screening tests that are used in primary care settings
[17, 10]. This suggests that our proposed approach deserves further
exploration and refinement as a primary care screening test tool. The benefit
of our approach over conventional tools is that it does not require exertion
from the patient, and does not depend on the training and technical skill of
clinical staff.

Our fev1 model can reliably distinguish between severe and mild airway
obstruction. However, it falls short in making finer distinctions between,
say, severe and very severe airway obstruction. To perform better on the fine
classification task, we could use a more sophisticated approach than our
current linear regression model. A more powerful model would have more
parameters and ultimately require that we collect much more
breathing/spirometry data for training. We could also consider a wider range
of tidal breathing predictor values, like variability measures and airflow
measures, which have been shown to correlate with airflow obstruction [15, 21,
16, 18].

##  VI Conclusion

In this paper, we studied the relationship between wearable device
measurements of tidal breathing and clinically-measured spirometer parameters.
Regression analysis revealed a moderate amount of correlation between the two
sets of parameters. Further, we built and evaluated simple models that used
wearable device measurements to assess lung health. The classifiers accurately
detected airway obstruction, and also provided a coarse stratification of the
severity of airway obstruction. The results of this study provide a strong
premise for further data collection and exploring the concept of wearable
devices for lung health monitoring.

##  VII Acknowledgements

This work was supported in part by the Munck-Pfefferkorn Education and
Research Fund. We would like to thank Jonathan Toledo and Orkan Sezer of
Dartmouth College and Henry Berube, Jeffrey Vonada and Jennifer Hilton-Hancock
of Dartmouth Hitchcock Medical Center for their invaluable support with data
collection. The views and conclusions contained in this document are those of
the authors and should not be interpreted as necessarily representing the
official policies, either expressed or implied, of the sponsors.

## References

  * [1] Andrea Aliverti and Antonio Pedotti  “Mechanics of breathing: new insights from new technologies”  Springer, 2014
  * [2] Cintia Laura Pereira Araujo et al.  “Pursed-lips breathing reduces dynamic hyperinflation induced by activities of daily living test in patients with chronic obstructive pulmonary disease: a randomized cross-over study”  In _Journal of rehabilitation medicine_ 47.10 Medical Journals Limited, 2015, pp. 957–962
  * [3] Mehdi Chlif et al.  “Effects of obesity on breathing pattern, ventilatory neural drive and mechanics”  In _Respiratory physiology & neurobiology_ 168.3 Elsevier, 2009, pp. 198–202
  * [4] Joana Cruz, Dina Brooks and Alda Marques  “Home telemonitoring in COPD: a systematic review of methodologies and patients’ adherence”  In _International journal of medical informatics_ 83.4 Elsevier, 2014, pp. 249–263
  * [5] Nermin Diab et al.  “Underdiagnosis and overdiagnosis of chronic obstructive pulmonary disease”  In _American Journal of Respiratory and Critical Care Medicine_ 198.9 American Thoracic Society, 2018, pp. 1130–1139
  * [6] Bruno-Pierre Dubé, Antoine Guerder, Capucine Morelot-Panzini and Pierantonio Laveneziana  “The clinical relevance of the emphysema-hyperinflated phenotype in COPD”  In _COPD Research and Practice_ 2.1 Springer, 2015, pp. 1
  * [7] Gary T Ferguson  “Why does the lung hyperinflate?”  In _Proceedings of the American Thoracic Society_ 3.2 American Thoracic Society, 2006, pp. 176–179
  * [8] Andrea S Gershon et al.  “Health services burden of undiagnosed and overdiagnosed COPD”  In _Chest_ 153.6 Elsevier, 2018, pp. 1336–1346
  * [9] Brian B Ghoshhajra et al.  “Direct chest area measurement: a potential anthropometric replacement for BMI to inform cardiac CT dose parameters?”  In _Journal of cardiovascular computed tomography_ 5.4 Elsevier, 2011, pp. 240–246
  * [10] Shamil Haroon, Rachel Jordan, Yemisi Takwoingi and Peymane Adab  “Diagnostic accuracy of screening tests for COPD: a systematic review and meta-analysis”  In _BMJ open_ 5.10 British Medical Journal Publishing Group, 2015
  * [11] Kjell Larsson  “Aspects on pathophysiological mechanisms in COPD”  In _Journal of internal medicine_ 262.3 Wiley Online Library, 2007, pp. 311–340
  * [12] Stephen W Littleton  “Impact of obesity on respiratory function”  In _Respirology_ 17.1 Wiley Online Library, 2012, pp. 43–49
  * [13] Stephen H Loring, Mauricio Garcia-Jacques and Atul Malhotra  “Pulmonary characteristics in COPD and mechanisms of increased work of breathing”  In _Journal of applied physiology_ 107.1 American Physiological Society, 2009, pp. 309–314
  * [14] Atul Malhotra et al.  “Research priorities in pathophysiology for sleep-disordered breathing in patients with chronic obstructive pulmonary disease. An official American Thoracic Society research statement”  In _American journal of respiratory and critical care medicine_ 197.3 American Thoracic Society, 2018, pp. 289–299
  * [15] MJ Morris and DJ Lane  “Tidal expiratory flow patterns in airflow obstruction.”  In _Thorax_ 36.2 BMJ Publishing Group Ltd, 1981, pp. 135–142
  * [16] Shayan Motamedi-Fakhr, Rachel C Wilson and Richard Iles  “Tidal breathing patterns derived from structured light plethysmography in COPD patients compared with healthy subjects”  In _Medical Devices (Auckland, NZ)_ 10 Dove Press, 2017, pp. 1
  * [17] Antonius Schneider et al.  “Diagnostic accuracy of spirometry in primary care”  In _BMC pulmonary medicine_ 9.1 Springer, 2009, pp. 31
  * [18] Ville-Pekka Seppä et al.  “Expiratory variability index (EVI) is associated with the severity of acute bronchial obstruction in small children: A proof-of-concept study”  In _Pediatric Allergy and Immunology_ Wiley Online Library, 2020
  * [19] Dave Singh et al.  “Global strategy for the diagnosis, management, and prevention of chronic obstructive lung disease: the GOLD science committee report 2019”  In _European Respiratory Journal_ 53.5 Eur Respiratory Soc, 2019
  * [20] Sean D Sullivan, Scott D Ramsey and Todd A Lee  “The economic burden of COPD”  In _Chest_ 117.2 Elsevier, 2000, pp. 5S–9S
  * [21] M Vitacca et al.  “Differences in spontaneous breathing pattern and mechanics in patients with severe COPD recovering from acute exacerbation”  In _European Respiratory Journal_ 13.2 Eur Respiratory Soc, 1999, pp. 365–370
  * [22] Henrik Watz et al.  “Spirometric changes during exacerbations of COPD: a post hoc analysis of the WISDOM trial”  In _Respiratory research_ 19.1 Springer, 2018, pp. 251


# Optimal Stochastic Coded Computation Offloading
in Unmanned Aerial Vehicles Network

Wei Chong Ng1,2 Wei Yang Bryan Lim1,2 Jer Shyuan Ng1,2 Suttinee Sawadsitang3
Zehui Xiong4 and Dusit Niyato5 _IEEE Fellow_
1Alibaba Group 2Alibaba-NTU JRI 3College of Arts  Media and Technology  Chiang
Mai University
4Singapore University of Technology and Design  Singapore
5School of Computer Science and Engineering  Nanyang Technological University
Singapore

###### Abstract

Today, modern unmanned aerial vehicles (UAVs) are equipped with increasingly
advanced capabilities that can run applications enabled by machine learning
techniques, which require computationally intensive operations such as matrix
multiplications. Due to computation constraints, the UAVs can offload their
computation tasks to edge servers. To mitigate stragglers, coded distributed
computing (CDC) based offloading can be adopted. In this paper, we propose an
Optimal Task Allocation Scheme (OTAS) based on Stochastic Integer Programming
with the objective to minimize energy consumption during computation
offloading. The simulation results show that amid uncertainty of task
completion, the energy consumption in the UAV network is minimized.

###### Index Terms:

Unmanned Aerial Vehicles, Coded Distributed Computing, Stochastic Integer
Programming, Task Allocation

##  I Introduction

In recent years, UAVs are becoming smaller in size and lower in production
cost, thus increasing the popularity of UAVs in the civil and commercial
fields. UAVs’ applications are also expanding in various areas, including
package delivery [1] and traffic monitoring [2]. In sensor networks, the UAVs
can act as an airborne base station to collect and process data from
terrestrial nodes [3]. UAV-based networked airborne computing allows
computation-intensive tasks to be executed on-board of the UAV directly
through networking [4].

Nevertheless, UAVs are faced with challenges in performing computationally
intensive tasks, due to limited battery lifetime and processing power [5]. As
such, the computation tasks can be offloaded to edge servers for processing.
In this paper, we focus mainly on the matrix multiplication task as it is
central to many modern computing applications, including machine learning and
scientific computing [6]. The matrix multiplication in the UAVs can be sped up
by scaling them out across many distributed computing nodes in base stations
(BSs) or edge servers [7] known as the workers. However, there is a major
performance bottleneck which is the latency in waiting for the slowest
workers, or “stragglers” to finish their tasks [7]. Coded distributed
computing is introduced to deal with stragglers in distributed high-
dimensional matrix multiplication. In CDC, the computation strategy for each
worker is carefully designed so that the UAV only needs to wait for the
fastest subset of workers before recovering the output [7]. The minimum number
of workers that the UAV has to wait to recover their results is known as the
recovery threshold. Nevertheless, there is still a probability that the number
of returned successful workers is less than the recovery threshold.

To address this uncertainty, we introduce the Optimal Task Allocation Scheme
(OTAS) in this paper. In OTAS, with the CDC technique such as PolyDot code
[7], the UAVs can choose to compute several copies locally or offload to the
BSs, where each copy is a sub-portion of matrices involved in the matrix
multiplication operation. The UAVs can retrieve the matrix multiplication
output by decoding all the returned copies. However, if the UAV decides to
offload to the BSs, there is an uncertainty that the copies cannot be returned
on time to the UAV due to delays and link failure [8], and the total copies
that the UAV has is less than the recovery threshold. The UAV has to pay a
correction cost to re-compute the number of shortfalls locally to match the
recovery threshold. Since the UAVs have to hover in the sky throughout the re-
computation, the correction cost is more expensive than the local computation
and offloading cost. OTAS can derive the UAV optimal decision through
Stochastic Integer Programming (SIP) formulation with two-stage recourse [9].
Extensive simulations are performed to evaluate the effectiveness of OTAS. The
results show that OTAS can minimize the total cost and the UAVs’ energy
consumption. The contributions of this paper are summarized as follows.

  * •

The proposed OTAS can minimize the overall costs incurred by the UAVs. At the
same time, it can also minimize all the UAVs’ energy consumption.

  * •

The formulated SIP model can find and achieve an optimal solution by taking
into account the workers’ uncertainty. Evidently, this is the first work that
utilizes the capability of SIP to address unpredictable UAV-enabled
distributed computing environments.

  * •

From the performance evaluation, we obtain a few insights such as the optimal
UAV task allocation. The performance comparison among the OTAS and the other
approaches is also presented.

The remainder of the paper is organized as follows: In Section II, we present
the network model. In Section III we formulate the problem. We discuss and
analyze the simulation result in Section IV. Section V concludes the paper.

##  II Network Model

We consider a network (Fig. 1) that consists of a set
$\mathcal{Y}=\\{1,\ldots,y,\ldots,Y\\}$ of UAVs to be deployed within the
coverage of a set $\mathcal{F}=\\{1,\ldots,f,\ldots,F\\}$ of macro BSs with
the height of $H_{f}$. All the UAVs will take off from their respective mobile
charging station which is located at $(\mathbf{x}_{y},\mathbf{y}_{y})$.
$\mathbf{x}_{y}$ and $\mathbf{y}_{y}$ are the x-y coordinates of UAV $y$
mobile charging station. UAV $y$ will take off vertically to the height of
$H_{y}$ and hover in the sky for purposes such as traffic monitoring [2].
$(\mathbf{x}_{y},\mathbf{y}_{y},H_{y})$ and
$(\mathbf{x}_{f},\mathbf{y}_{f},H_{f})$ are the three-dimensional coordinates
of the UAV $y$ and BS $f$ respectively, where $H_{y}>H_{f}$ to maintain a
line-of-sight (LoS) communication link between UAV $y$ and BS $f$. Due to the
hovering capability, we consider only the rotary-wing UAVs [10].

Figure 1: An illustrative example of UAV network with $Y=6$, $F=2$.

After the UAVs reach their respective heights, they can receive and process
tasks that use the data, e.g., captured and collected by their sensors. In
this paper, we consider the task that UAV $y$ compute is the matrix-matrix
product AB involving the two matrices A and B. Matrix multiplication plays an
essential role in various scientific disciplines as it is regarded as the
primary tool for many other computations in areas such as aerodynamic
computations and machine learning [6]. However, the UAV has limited computing
and storage capability [11]. Therefore, the UAV can choose to offload a
portion or the whole matrix multiplication to the independently owned BS [11].

In order to speed up matrix multiplication, massive parallelization has
emerged as a solution. However, it has a computational bottleneck due to
straggling or faults. Coded computation is introduced to make matrix
multiplications resilient to faults and delays, i.e., Polynomial codes [7]. In
PolyDot codes, the system model typically consists of a master node, multiple
worker nodes, and a fusion node. In our proposed OTAS, we consider that the
UAV is the master and the fusion node. BS $f$ can partition its CPU computing
capability into $n_{f}$ portions. Each portion is treated as the worker and
has the computation capability of $\frac{\tau_{f}}{n_{f}}$, where $\tau_{f}$
denotes the CPU computation capability of the BS $f$ (in CPU cycles per
second). Each worker has a storage constraint that limits the worker to store
only $m$ fraction of matrices A and B [6].

In general, three scenarios may occur.

  * •

The UAV can compute all copies locally.

  * •

The UAV can offload a number of copies to the BS.

  * •

The UAV can compute some copies locally and offload some copies to the BS.

$M_{y}^{(l)}$ indicates the number of copies that UAV $y$ computes locally and
$M_{y,f}^{(o)}$ denotes the number of copies that are offloaded to BS $f$. The
final output can be decoded from all the return copies
($M_{y}^{(l)}+M_{y,f}^{(o)}$).

The definition of copy, successful workers, recovery threshold, and shortfall
are given as follows.

Definition 1. [Copy] $m$-th fraction of matrices A and B [6].

Definition 2. [Successful workers] Workers that finish their computation task
successfully and the task is received by the UAV.

Definition 3. [Recovery threshold] The recovery threshold is the worst-case
minimum number of successful workers required by the UAV to complete the
computation successfully [6].

Definition 4. [Shortfall] There exists a shortfall when the total returned
copies from the local computation and from the workers are less than the
recovery threshold.

Following [6], two $N\times N$ square matrices A and B are considered. Note
that our model can be applied to other matrices, e.g., Non square matrices.
Each of the matrices A and B is sliced both horizontally and vertically. For
example, A is sliced into $\frac{N}{t}\times\frac{N}{s}$ matrices and B is
sliced into $\frac{N}{s}\times\frac{N}{t}$. We choose $s$ and $t$ such that it
satisfies $st~{}=m$ [6]. The recovery threshold $k$ is defined [6] as:

$k=t^{2}(2s-1).$ (1)

The processing by the workers may take a longer time when it is currently
occupied with some other tasks. Therefore, the processing in the offloaded
tasks are perceived to have failed if the duration exceeds the threshold time
limit [12]. To recover the computed task, the sum of returned offloaded copies
and locally computed copies must be greater than or equal to the recovery
threshold $k$.

Similar to [6], the UAV uses $d_{enc}=~{}N^{2}(M_{y}^{(l)}+M_{y,f}^{(o)})$
symbols for encoding of matrices and
$d_{dec}=~{}N^{2}t^{2}(2s-1)\log^{2}t^{2}(2s-1)$ to decode the returned
matrices. Each copy contains $m$-th fraction of matrices A and B. UAV will
transmit $d_{comm}^{to}=~{}\frac{N^{2}}{m}$ symbols to each of the worker.
Each copy requires $d_{cmp}=~{}\frac{N^{3}}{mt}$ symbols for computation.
After computation is completed, the worker will send
$d_{comm}^{fr}=~{}\frac{N^{2}}{t^{2}}$ symbols back to the UAV.

As described earlier, the UAV can choose to offload a portion of the task or
the whole task to the workers. For example, when the recovery threshold $k=4$
and the UAV decide to offload two copies of the task for the workers to
compute $M_{y,f}^{(o)}=~{}2$. Therefore, the UAV has to compute at least two
more copies locally to match the recovery threshold $M_{y}^{(l)}=~{}2$. The
UAV will hover in the sky for a threshold time limit $t^{thresh}_{y}$ to wait
for the offloaded copies to return. In this paper, $t^{thresh}_{y}$ is set as
the worst-case scenario, i.e., time required to compute all copies locally by
the UAV. However, there is a probability that the workers may fail, i.e., the
computed task is not returned to the UAV after $t^{thresh}_{y}$. As a result,
the UAV cannot successfully complete the full task if the total returned
copies are less than 4. When the UAV fails to receive a copy, it means that
there exists a shortfall, and hence, the UAV has to recompute the number of
shortfalls locally. In the meantime, the UAV has to continue hover in the sky
when performing the re-computation. In this paper, we use SIP to model the
uncertainty problem to optimize the number of copies to compute locally
$M_{y}^{(l)}$ and offload $M_{y,f}^{(o)}$. At the same time, it will also
minimize the overall cost incurred by the UAVs. In the following, we discuss
the propulsion, computation, and communication related energy costs for the
UAVs.

###  II-A UAV Hovering Energy

The propulsion energy consumption is needed to provide the UAV with sufficient
thrust to support its movement. The propulsion power of a rotary-wing UAV with
speed $V$ can be modeled as follows [10]:

$P(V)=P_{0}(1+\frac{3V^{2}}{U^{2}_{tip}})+P_{1}\biggl{(}\sqrt{1+\frac{V^{4}}{4v_{0}^{4}}}-\frac{V^{2}}{2v_{0}^{2}}\biggr{)}^{\frac{1}{2}}+\frac{1}{2}d_{0}\rho\mathbb{A}V^{3},$
(2)

where

$\displaystyle
P_{0}=\frac{\delta}{8}\rho\mathbf{s}\mathbb{A}\triangle^{3}R^{3},$ (3)
$\displaystyle P_{1}=(1+\mathbf{k})\frac{W^{3/2}}{\sqrt{2\mathbb{A}\rho}}.$
(4)

$P_{0}$ and $P_{1}$ are two constants related to UAV’s weight, rotor radius,
air density, etc. $U_{tip}$ denotes the tip speed of the rotor blade, $v_{0}$
is known as the mean rotor induced velocity in hover, $d_{0}$ and $s$ are the
fuselage drag ratio and rotor solidity, respectively. $\rho$ and $\mathbb{A}$
are the air density and rotor disc area, respectively. $W$ is the UAV weight,
$\delta$ is the profile drag coefficient, and $\triangle$ denotes blade
angular velocity. By substituting $V=0$ into (2) [10], we obtain the power
consumption for hovering status as follows:

$P_{h}=P_{0}+P_{1}.$ (5)

###  II-B Local Computing Model

When one copy is processed locally, the local computation execution time of
UAV $y$ is expressed as [13]:

$t_{y}^{local}=\frac{\mathcal{C}_{y}d(\frac{N^{3}}{mt})}{\tau_{y}},$ (6)

where $\mathcal{C}_{y}$ is the number CPU cycles needed to process a bit,
$\tau_{y}$ denotes the total CPU computing capability of UAV $y$, and
$d(\cdot)$ is a function to translate the number of symbols to the number of
bits for computation. The UAV $y$ takes $t_{y}^{enc}$ seconds to encode one
copy of the matrices, and it is expressed as follows:

$t_{y}^{enc}=\frac{\mathcal{C}_{y}d(N^{2})}{\tau_{y}}.$ (7)

After the UAV $y$ obtains at least $k$ copies, it will take $t_{y}^{dec}$ secs
to decode. $t_{y}^{dec}$ is defined as follows:

$t_{y}^{dec}=\frac{\mathcal{C}_{y}d(N^{2}t^{2}(2s-1)\log^{2}t^{2}(2s-1))}{\tau_{y}}.$
(8)

###  II-C UAV Communication Model

We assume that each UAV is allocated with an orthogonal spectrum resource
block to avoid the co-interference among the UAVs [14]. The transmission rate
from UAV $y$ to BS $f$ can be represented as [5]:

$r_{y,f}=B_{y}\log_{2}(1+P_{y}^{C}h_{y,f}/N_{o}),$ (9)

where the wireless transmission power of the UAV $y$ is expressed as
$P_{y}^{C}$. $h_{y,f}$ is the channel gains, and $N_{0}$ is the variance of
complex white Gaussian noise. UAV to BS communication is most likely to be
dominated by LoS links. Therefore, the air-to-ground channel power gain from
UAV $y$ to BS $f$ can be modeled as follows [15]:

$h_{y,f}=\frac{\beta_{0}}{\mathbb{D}_{y,f}^{2}},$ (10)

where

$\mathbb{D}_{y,f}^{2}=(\mathbf{x}_{y}-\mathbf{x}_{f})^{2}+(\mathbf{y}_{y}-\mathbf{y}_{f})^{2}+(H_{y}-H_{f})^{2},$
(11)

$\mathbb{D}_{y,f}$ denotes the distance between UAV $y$ and BS $f$ and
$\beta_{0}$ represents the reference channel gain at distance $d_{0}=1$m in an
urban area [15]. The transmission time to offload one copy of matrix from UAV
$y$ to BS $f$ can be given as follows:

$t_{y,f}^{to}=\frac{d(\frac{N^{2}}{m})}{r_{y,f}}.$ (12)

The energy $e_{y}$ required by UAV $y$ to receive data from the BS $f$ is
defined as follows [16]:

$\displaystyle e_{y}=P_{y}^{r}\frac{d(\frac{N^{2}}{t^{2}})}{r_{f,y}},$ (13)

where $P_{y}^{r}$ is the receiving power of UAV $y$ and $r_{f,y}$ is similar
to (9).

##  III Problem formulation

This section introduces the Deterministic Integer Programming (DIP) and SIP to
optimize the number of copies to compute locally and offload by minimizing the
UAVs’ total cost.

###  III-A Deterministic Integer Programming System Model

In an ideal case when the number of shortfalls is precisely known, the UAVs
can choose the exact number of copies to compute locally or offload.
Therefore, the correction for the shortfall is not needed, and the correction
cost is zero. Similar to [17], the cost function is proportional to their
offloaded data and to their demand for consuming computation resources. In
total, four types of payments are considered in DIP. It will cost
$\bar{C}_{y}$ for UAV $y$ to compute one copy locally. $\underline{C}_{y}$ is
the cost for UAV $y$ to hover in the sky for $t^{thresh}_{y}$ seconds. Since
$\underline{C}_{y}$ exists in all three scenarios as described in Section II,
we can completely remove this cost from the formulation. The BS $f$ charges
the UAV with $C_{y,f}$ for each copy that is offloaded. $\hat{C_{y}}$ is the
decoding cost to retrieve the final output.

  * •

$\bar{C}_{y}$ denotes the UAV $y$ local computational and encoding cost for
computing of one copy, i.e.,

$\bar{C}_{y}=\alpha_{1}(t_{y}^{local}+t_{y}^{enc}),$ (14)

where $\alpha_{1}$ is the cost coefficient associated to the energy
consumption.

  * •

$C_{y,f}$ denotes the offloading cost and it consists of three parts. The
first part is related to the transmission and encoding delay. The second part
is the UAV $y$ receiving energy cost and the last part $C_{f}$ is the service
cost for BS $f$. It is modeled as follows:

$C_{y,f}=\alpha_{1}(t_{y,f}^{to}+t_{y}^{enc})+\alpha_{2}e_{y}+C_{f},$ (15)

where $\alpha_{2}$ is the cost coefficient with a similar role to
$\alpha_{1}$.

  * •

$\hat{C_{y}}$ denotes the UAV $y$ decoding cost for the returned matrices as
follows:

$\hat{C_{y}}=\alpha_{1}t_{y}^{dec}.$ (16)

A DIP can be formulated and minimize the total cost of the UAVs as follows:

Minimise:

$\displaystyle\sum_{y\in\mathcal{Y}}(M_{y}^{(l)}\bar{C}_{y}+\hat{C}_{y})+\sum_{y\in\mathcal{Y}}\sum_{f\in\mathcal{F}}M_{y,f}^{(o)}C_{y,f},$
(17)

subject to:

$\displaystyle M_{y}^{(l)}+\sum_{f\in\mathcal{F}}M_{y,f}^{(o)}-S_{y}\geq k,$
$\displaystyle\forall y\in\mathcal{Y},$ (18)
$\displaystyle\sum_{y\in\mathcal{Y}}M_{y,f}^{(o)}\leq n_{f},$
$\displaystyle\forall f\in\mathcal{F}.$ (19)

The objective function in (17) is to minimize UAVs’ total cost. The first part
is to minimize the UAVs’ local computation cost and the second part is to
minimize the UAVs’ offloading cost. The constraint in (18) ensures that the
number of copies computed locally and offloaded have to be at least equal to
or larger than the recovery threshold $k$, and $S_{y}$ is the number of
shortfalls. The constraint in (19) ensures that the total number of copies
offloaded to BSs must not exceed the total number of workers.

###  III-B Stochastic Integer Programming System Model

If the number of shortfalls cannot be known precisely, the deterministic
optimization formulation defined in (17) - (19) is no longer applicable.
Therefore, SIP with two-stage recourse is developed. The first stage is the
assignment stage, whereby the UAV decides the number of copies to be computed
locally and offloaded. The decision will be made based on the available cost
information and the probability distribution of the shortfall, which refers to
at least one copy that fails to return to the UAV after the copy has been
offloaded. It can lead to the situation that the total number of copies that
the UAV receives is less than $k$. The second stage is the correction stage.
After the actual shortfall is observed, the UAV will perform a correction
action to compute the number of shortfalls locally so that the total number of
copies can match with $k$. It costs the UAV $\tilde{C}_{y}$ to compute each
copy locally for the shortfall. Similar to DIP, in order to retrieve the final
matrix multiplication output, it will cost UAV $y$ $\hat{C_{y}}$ for decoding.

Let $\omega=(\mathbb{F}_{1},\ldots,\mathbb{F}_{y})$ denotes the scenario of
all UAVs and a set of scenarios is denoted by $\Omega$, i.e.,
$\omega\in\Omega$ [1]. $\mathbb{F}_{y}$ represents a binary parameter of the
shortfall from UAV $y$. When $\mathbb{F}_{y}=1$ means that, from the copies
that UAV $y$ had offloaded, at least $A_{y}$ copy did not return. As a result,
the total number of copies that it currently has is less than $k$, and
$\mathbb{F}_{y}=0$ otherwise. For example, suppose the service provider owns
three UAVs. In this case, the shortfall scenario is denoted by
$\omega=(\mathbb{F}_{1},\mathbb{F}_{2},\mathbb{F}_{3})$ in which (1,0,0)
indicates that UAV 1 have shortfall due to the copy failure to return while
UAV 2 and 3 do not. Let $A=(A_{1},\ldots,A_{Y})$ be a list representing the
number of shortfalls whereby $A_{y}\leq k$.

The objective function of SIP is the same as DIP with an additional correction
cost.

  * •

$\tilde{C}_{y}$ is the correction cost, and it consists of two parts. The
first part denotes the correction cost for one copy if the UAV $y$ has a
shortfall. The second part is the penalty cost for the UAV to stay hover in
the sky when there is a shortfall. The penalty cost can be from delay of
completing the task, additional computation energy and some other opportunity
costs, e.g., UAV cannot process the next task immediately. The correction cost
is given as follows:

$\displaystyle\tilde{C}_{y}=\alpha_{1}(t_{y}^{local}+t_{y}$
${}^{enc})+\alpha_{2}P_{h}t_{y}^{local},$ (20)

We formulate the task allocation as a two-stage SIP model. There are three
decision variables in this model.

  * •

$M_{y}^{(l)}$ indicates the number of copies to compute locally by UAV $y$.

  * •

$M_{y,f}^{(o)}$ indicates the number of copies to be offloaded to the BS $f$
by UAV $y$.

  * •

$M_{y}^{(L)}$ indicates the number of copies to compute locally by UAV $y$
when there is a shortfall.

The objective function given in (III-B) and (22) is to minimize the total cost
of UAVs. The expressions in (III-B) and (22) represent the first stage and
second stage objectives, respectively. In (22), $\mathcal{P}(\omega)$ denotes
the probability if scenario $\omega\in\Omega$ is realized and in reality,
$\mathcal{P}(\omega)$ can be obtained from the historical records. The SIP
formulation can be expressed as follows:

Minimise:

$\displaystyle\sum_{y\in\mathcal{Y}}(M_{y}^{(l)}\bar{C}_{y}+\hat{C}_{y})+\sum_{y\in\mathcal{Y}}\sum_{f\in\mathcal{F}}M_{y,f}^{(o)}$
$\displaystyle C_{y,f}\noindent$ $\displaystyle+\mathbb{E}$
$\displaystyle\biggl{[}\mathcal{Q}(M^{(o)}_{y,f}(\omega))\biggr{]},$ (21)

where

$\displaystyle\mathcal{Q}(M_{y,f}^{(o)}(\omega))=\sum_{y\in\mathcal{Y}}\sum_{\omega\in\Omega}(\mathbb{F}_{y}(\omega)\mathcal{P}(\omega)M_{y}^{(L)}(\omega)\tilde{C}_{y}),$
(22)

subject to:

$\displaystyle M_{y}^{(l)}+\sum_{f\in\mathcal{F}}M_{y,f}^{(o)}\geq k,$
$\displaystyle\forall y\in\mathcal{Y},$ (23) $\displaystyle
M_{y}^{(l)}+\sum_{f\in\mathcal{F}}M_{y,f}^{(o)}+M_{y}^{(L)}(\omega)-A_{y}(\omega)\geq
k,$ $\displaystyle\forall y\in\mathcal{Y},\forall\omega\in\Omega,$ (24)
$\displaystyle\sum_{y\in\mathcal{Y}}M_{y,f}^{(o)}\leq n_{f},$
$\displaystyle\forall f\in\mathcal{F}.$ (25)

The constraint in (23) ensures that the number of copies computed locally and
offloaded should be at least $k$. The constraint in (24) ensures that if the
UAV has a shortfall, the UAV has to compute the shortfall locally to match
$k$. The constraint (25) is the same as constraint (19).

##  IV Simulation result and analysis

In this simulation, we consider the system model with ten UAVs and two BSs. A
conceptual illustration of the system model is shown in Fig. 2. Figure 2 is an
x-y plane that shows the locations of the UAVs and the BSs. They are randomly
allocated in a $1000\times 1000$ $m^{2}$. Each small grid is $25\times 25$
$m^{2}$. We consider the case with $m=2$ [6]. Therefore, we can substitute
$s=\frac{2}{t}$ into equation (1), and we obtain the following:

$k=4t-t^{2},$ (26)

by differentiating equation (26) with respect to the variable $t$ and equating
the result to zero. Then, we are able to obtain $t=2$, $s=1$, and $k=4$. To
solve SIP, We assume that the probability distribution of all scenarios in set
$\Omega$ are known [18], then, the complexity of the problem depends on the
total number of scenarios in stage two [18]. For example in Section III, the
complexity for the formulated two-stage SIP is $|\Omega|$. The simulation
parameters are summarized in Table I and their values are from [5, 10]. For
the presented experiments, we implement the SIP model using GAMS script [19].

Figure 2: x-y coordinates of all the UAVs and BSs.

Figure 3: The optimal solution in a simple UAV network.

Figure 4: Total cost of the network as the number of UAVs is varied.

Figure 5: EVF comparing with SIP.

First, the cost structure of the network is studied. As an illustration, a
primitive UAV network is considered with zero local computation, and
$n_{1}=n_{2}=100$. Furthermore, we first consider 10 UAVs and 1 BS. Since
$k=4$, the number of copies that are offloaded begins with 4. In Fig. 5, the
costs in the first and second stages and the total cost under the different
number of the offloaded copies are presented. As expected, the first stage’s
cost increases as the number of copies that are offloaded increases. However,
the cost in the second stage after knowing the actual shortfall decreases, as
the number of copies that are offloaded increases since the UAV needs to
perform fewer re-computation to match with the shortfall. Figure 5 shows the
optimal solution in this simple network. It can be identified that even in
this simple network, the optimal solution is not trivial to obtain due to the
uncertainty of shortfall. Therefore, SIP formulation is required to guarantee
the minimum cost to the UAVs.

TABLE I: Experiment parameters Parameter | Values
---|---
Length of the matrix, $N$ | 1000
UAV weight in $kg$, $W$ | 10.2
Air density in $kg/m^{3}$, $\rho$ | 1.225
Rotor radius in meter, $R$ | 0.5
Rotor disc area in $m^{2}$, $\mathbb{A}$ | 0.79
Blade angular velocity in radians/second, $\triangle$ | 400
Tip speed of the rotor blade, $U_{tip}$ | 200
Number of blades, $b$ | 4
Blade or aerofoil chord length, $c$ | 0.0196
Rotor solidity, defined as the ratio of the total blade area $bcR$ to the disc area $\mathbb{A}$, $\mathbf{s}$ | 0.05
Fuselage drag ratio, $d_{0}$ | 0.3
Mean rotor induced velocity in hover, $v_{0}$ | 7.2
Profile drag coefficient, $\delta$ | 0.012
Incremental correction factor to induced power, $\mathbf{k}$ | 0.1
UAV hover height in meter, $H_{y}$ | 100
BS height in meter, $H_{f}$ | 20
UAV bandwidth in $MHz$, $B_{y}$ | 2
UAV transmit power in $mW$, $P_{y}^{C}$ | 32
UAV receiving power in $mW$, $P_{y}^{r}$ | 32
White Gaussian noise in $dBm$, $N_{0}$ | -100
UAV computation power in $GHz$, $\tau_{y}$ | 1
BS computation power in $GHz$, $\tau_{f}$ | 20
Number of CPU cycles needed to process a bit, $\mathcal{C}_{y}$ | 20
Channel gains, $\beta_{0}$ | -60dB
Number of workers in $BS_{f}$, $n_{f}$ | 15
Service cost for BS $f$, $C_{f}$ in $ | 0.2
$\alpha_{1}$ | 0.6
$\alpha_{2}$ | 0.0004

Next, we evaluate the effect of number of UAVs in the network. The network is
initialized with only $UAV_{1}$. The total cost of the network is monitored as
the number of UAVs joining the network increases. Figure 5 displays the result
from this set-up. When the number of UAVs in the network is more than 6, we
can observe a sharp increase in this network’s cost far beyond the preceding
trend, as indicated by the dotted line. The total number of workers in each BS
that can work on the matrix multiplication is $n_{f}$. In this case, all the
workers from the BSs are occupied when the number of UAVs exceeded 6. Instead
of offloading to the BSs, the rest of the UAVs have to compute the copy
locally, leading to a sharp increase in cost.

Then, we compare the SIP with Expected-Value Formulation (EVF) [20]. Expected-
value formulation uses the average values of shortfall and solves a DIP. We
vary the price of the correction action to compare the difference between EVF
and SIP. We set $n_{1},n_{2}=~{}100$. Figure 5 depicts the comparison result.
For EVF, the number of copies in the first stage ($M_{y}^{(l)}+M_{y,f}^{(o)}$)
is fixed using the average value of shortfall, an approximation scheme. As
shown in the result, EVF cannot adapt to the change in cost. On the other
hand, SIP can always achieve the optimal solution to reduce the shortfall
cost.

##  V Conclusion

In this paper, we have proposed an OTAS with the objective to minimize UAV
energy consumption in the PolyDot code based computation offloading. To
account for the uncertainty of task completion, we have formulated the OTAS as
the two-stage stochastic programming. The OTAS minimizes the total cost and
the UAVs’ energy consumption. As compared to DIP, OTAS based on SIP can
achieve the optimal solution as it is able to adapt to changes in probability
of task failure.

## Acknowledgement

This research is supported, in part, by Alibaba Group through Alibaba
Innovative Research (AIR) Program and Alibaba-NTU Singapore Joint Research
Institute (JRI), the National Research Foundation, Singapore under the AI
Singapore Programme (AISG) (AISG2-RP-2020-019), WASP/NTU grant M4082187
(4080), Singapore Ministry of Education (MOE) Tier 1 (RG16/20) and SUTD SRG-
ISTD-2021-165.

## References

  * [1] S. Sawadsitang, D. Niyato, P. Tan, and P. Wang.  Joint ground and aerial package delivery services: A stochastic optimization approach.  IEEE Transactions on Intelligent Transportation Systems, 20(6):2241–2254, 2019.
  * [2] Kapseong Ro, Jun-Seok Oh, and Liang Dong.  Lessons learned: Application of small uav for urban highway traffic monitoring.  In 45th AIAA aerospace sciences meeting and exhibit, page 596, 2007\.
  * [3] Haichao Wang, Jinlong Wang, Jin Chen, Yuping Gong, and Guoru Ding.  Network-connected uav communications: Potentials and challenges.  China Communications, 15(12):111–121, 2018.
  * [4] B. Wang, J. Xie, K. Lu, Y. Wan, and S. Fu.  Coding for heterogeneous uav-based networked airborne computing.  In 2019 IEEE Globecom Workshops (GC Wkshps), pages 1–6, 2019.
  * [5] Jienan Chen, Siyu Chen, Siyu Luo, Qi Wang, Bin Cao, and Xiaoqian Li.  An intelligent task offloading algorithm (itoa) for uav edge computing network.  Digital Communications and Networks, 6(4):433–443, 2020.
  * [6] Sanghamitra Dutta, Mohammad Fahim, Farzin Haddadpour, Haewon Jeong, Viveck Cadambe, and Pulkit Grover.  On the optimal recovery threshold of coded matrix multiplication.  IEEE Transactions on Information Theory, 66(1):278–301, 2019.
  * [7] Qian Yu, Mohammad Maddah-Ali, and Salman Avestimehr.  Polynomial codes: an optimal design for high-dimensional coded matrix multiplication.  In Advances in Neural Information Processing Systems, pages 4403–4413, 2017.
  * [8] Baoqian Wang, Junfei Xie, Kejie Lu, Yan Wan, and Shengli Fu.  On batch-processing based coded computing for heterogeneous distributed computing systems.  arXiv preprint arXiv:1912.12559, 2019.
  * [9] John R Birge and Francois Louveaux.  Introduction to stochastic programming.  Springer Science & Business Media, 2011.
  * [10] Yong Zeng, Jie Xu, and Rui Zhang.  Energy minimization for wireless communication with rotary-wing uav.  IEEE Transactions on Wireless Communications, 18(4):2329–2345, 2019\.
  * [11] Qiyu Hu, Yunlong Cai, Guanding Yu, Zhijin Qin, Minjian Zhao, and Geoffrey Ye Li.  Joint offloading and trajectory design for uav-enabled mobile edge computing systems.  IEEE Internet of Things Journal, 6(2):1879–1892, 2018.
  * [12] Sanghamitra Dutta, Viveck Cadambe, and Pulkit Grover.  Coded convolution for parallel and distributed computing within a deadline.  In 2017 IEEE International Symposium on Information Theory (ISIT), pages 2403–2407. IEEE, 2017.
  * [13] Quoc-Viet Pham, Ming Zeng, Rukhsana Ruby, Thien Huynh-The, and Won-Joo Hwang.  Uav communications for sustainable federated learning.  arXiv preprint arXiv:2103.11073, 2021.
  * [14] Zhenyu Zhou, Junhao Feng, Zheng Chang, and Xuemin Shen.  Energy-efficient edge computing service provisioning for vehicular networks: A consensus admm approach.  IEEE Transactions on Vehicular Technology, 68(5):5087–5099, 2019\.
  * [15] Meng Hua, Yi Wang, Qingqing Wu, Haibo Dai, Yongming Huang, and Luxi Yang.  Energy-efficient cooperative secure transmission in multi-uav-enabled wireless networks.  IEEE Transactions on Vehicular Technology, 68(8):7761–7775, 2019\.
  * [16] Jer Shyuan Ng, Wei Yang Bryan Lim, Hong-Ning Dai, Zehui Xiong, Jianqiang Huang, Dusit Niyato, Xian-Sheng Hua, Cyril Leung, and Chunyan Miao.  Joint auction-coalition formation framework for communication-efficient federated learning in uav-enabled internet of vehicles.  IEEE Transactions on Intelligent Transportation Systems, 2020.
  * [17] Giorgos Mitsis, Eirini Eleni Tsiropoulou, and Symeon Papavassiliou.  Data offloading in uav-assisted multi-access edge computing systems: A resource-based pricing and user risk-awareness approach.  Sensors, 20(8):2434, 2020.
  * [18] Martin Dyer and Leen Stougie.  Computational complexity of stochastic programming problems.  Mathematical Programming, 106(3):423–432, December 2006.
  * [19] Debabrata Chattopadhyay.  Application of general algebraic modeling system to power system optimization.  IEEE Transactions on Power Systems, 14(1):15–22, 1999.
  * [20] S. Chaisiri, Bu-Sung Lee, and D. Niyato.  Optimal virtual machine placement across multiple cloud providers.  In 2009 IEEE Asia-Pacific Services Computing Conference (APSCC), pages 103–110, 2009.