# Computationally Efficient Safe Reinforcement Learning for Power Systems

Daniel Tabas and Baosen Zhang This work is partially supported by the National
Science Foundation Graduate Research Fellowship Program under Grant No.
DGE-1762114 and NSF grant ECCS-1930605. Any opinions, findings, conclusions,
or recommendations expressed in this material are those of the authors and do
not necessarily reflect the views of the National Science Foundation.Authors
are with the department of Electrical and Computer Engineering, University of
Washington, Seattle, WA, United States. {dtabas<EMAIL_ADDRESS>
###### Abstract

We propose a computationally efficient approach to safe reinforcement learning
(RL) for frequency regulation in power systems with high levels of variable
renewable energy resources. The approach draws on set-theoretic control
techniques to craft a neural network-based control policy that is guaranteed
to satisfy safety-critical state constraints, without needing to solve a model
predictive control or projection problem in real time. By exploiting the
properties of robust controlled-invariant polytopes, we construct a novel,
closed-form “safety-filter” that enables end-to-end safe learning using any
policy gradient-based RL algorithm. We then apply the safety filter in
conjunction with the deep deterministic policy gradient (DDPG) algorithm to
regulate frequency in a modified 9-bus power system, and show that the learned
policy is more cost-effective than robust linear feedback control techniques
while maintaining the same safety guarantee. We also show that the proposed
paradigm outperforms DDPG augmented with constraint violation penalties.

##  I INTRODUCTION

Power systems are a quintessential example of safety-critical infrastructure,
in which the violation of operational constraints can lead to large blackouts
with high economic and human cost. As variable renewable energy resources are
integrated into the grid, it becomes increasingly important to ensure that the
system states, such as generator frequencies and bus voltages, remain within a
“safe” region defined by the operators [1].

The design of safe controllers concerns the ability to ensure that an
uncertain dynamical system will satisfy hard state and action constraints
during execution of a control policy [2, 3]. Recently, set-theoretic control
[4] has been applied to a wide range of safety-critical problems in power
system operation [5, 6]. This approach involves computing a _robust
controlled-invariant set_ (RCI) along with an associated control policy which
is guaranteed to keep the system state inside the RCI [4, 7, 8]. If the RCI is
contained in the feasible region of the (safety-critical) state constraints,
then the associated control policy is considered to be safe.

However, the set-theoretic approach requires several simplifying assumptions
for tractability, leading to controllers with suboptimal performance. First,
the disturbances to the system are assumed to be bounded in magnitude but
otherwise arbitrary [4, 6]. Second, the RCIs must be restricted to simple
geometric objects such as polytopes or ellipsoids [9]. Third, many approaches
select an RCI and control policy in tandem, which usually requires the control
policy to be linear and forces a tradeoff between performance and robustness
[10, 11, 12]. Fourth, nonlinear systems must be treated as linear systems plus
an unknown-but-bounded linearization error [5].

Once an RCI is generated using the conservative assumptions listed above,
data-driven approaches can use learning to deal with the true behavior of the
disturbances and nonlinearities without risk of taking unsafe actions [13, 14,
15, 16]. However, these techniques require solving a model predictive control
(MPC) or projection problem each time an action is executed, which may be too
computationally expensive. Several approaches that avoid repeatedly solving an
optimization problem have also been proposed. One such approach involves
tracking the vertices of the set of safe actions, and using a neural network
to specify an action by choosing convex weights on these vertices. However,
this is only possible when the RCI has exceedingly simple geometry [17]. Other
strategies only guarantee safety in expectation, and do not rule out
constraint violations in every situation [18, 19]. Controllers with Lyapunov
stability or robust control guarantees have also been proposed [20, 21, 22],
but stability does not always translate to constraint satisfaction.

In this paper, we present a method to design safe, data-driven, and closed-
form control policies for frequency regulation in power systems. Our approach
combines the advantages of set-theoretic control and learning. In particular,
we use simple linear controllers to find a maximal RCI, and then use
reinforcement learning (RL) to train a neural network-based controller that
improves performance while maintaining safety. The safety of this control
policy is accomplished by constraining the output of the neural network to the
present set of safe actions. By leveraging the structure of polytopic RCIs, we
construct a closed-form safety filter to map the neural network’s output into
the safe action set without solving an MPC or projection problem. The safety
filter is differentiable, allowing end-to-end training of the neural network
using any policy gradient-based RL algorithm. We demonstrate our proposed
control design on a frequency regulation problem in a 9-bus power system model
consisting of several generators, loads, and inverter-based resources (IBRs).
The simulation results demonstrate that our proposed policy maintains safety
and outperforms safe linear controllers without repeatedly solving an
optimization problem in real time.

We focus on applying our algorithm to the problem of primary frequency control
in power systems. Frequency is a signal in the grid that indicates the balance
of supply and demand. Generators typically respond to the change in frequency
by adjusting their power output to bring the frequency back to nominal (e.g.,
60 Hz in the North American system) [23, 24]. For conventional generators,
these responses are limited to be linear (possibly with a dead-band). In
contrast, IBRs such as solar, wind and battery storage can provide almost any
desired response to frequency changes, subject to some actuation constraints
[25]. Currently, however, these resources still use linear responses, largely
because of the difficulty in designing nonlinear control laws. Recently, RL
based methods have been introduced in the literature (see, e.g. [26] and the
references within). However, most approaches treat safety and constraint
satisfaction as soft penalties, and cannot provide any guarantees [26, 27,
28].

The rest of the paper is organized as follows. Section II introduces the power
system model and formulates the problem of safety-critical control from a set-
theoretic perspective. Section III describes the proposed controller design.
Section IV presents simulation results for the modified 9-bus power system.

##  II MODEL AND PROBLEM FORMULATION

###  II-A Model assumptions

In this paper we are interested in a linear system with control inputs and
disturbances. We write the system evolution as

$x_{t+1}=Ax_{t}+Bu_{t}+Ed_{t},$ (1)

where $x_{t}\in\mathbb{R}^{n}$, $u_{t}\in\mathbb{R}^{m}$ and
$d_{t}\in\mathbb{R}^{p}$ are vectors of the state variables, control inputs,
and disturbances at time $t$. We assume the disturbance $d_{t}$ is bounded but
otherwise can take arbitrary values. More precisely, we assume that $d_{t}$
lies in a compact set. This boundedness assumption on $d_{t}$ is fairly
general, since it allows the disturbances to capture uncontrolled input into
the system, model uncertainties in $A$, $B$, and $E$, and linearization error.
For more compact notation, we will sometimes summarize (1) as
$x^{+}=f(x_{t},u_{t},d_{t})$.

The constraints on inputs are $u_{t}\in\textit{{U}}\subset\mathbb{R}^{m}$ and
$d_{t}\in\textit{{D}}\subset\mathbb{R}^{p}$ for all $t$. The sets U and D are
assumed to be polytopes, defined as the bounded intersection of a finite
number of halfspaces or linear inequalities [29]. Specifically, U and D are
defined as

U $\displaystyle=\\{u\in\mathbb{R}^{m}\mid-\bar{u}\leq
V_{u}u\leq\bar{u}\\}\text{ and}$ (2) D
$\displaystyle=\\{d\in\mathbb{R}^{p}\mid-\bar{d}\leq V_{d}d\leq\bar{d}\\}.$
(3)

In safety-critical control problems such as frequency regulation, operators
want to keep the system states within hard constraints. For example,
frequencies are generally kept within a tenth of a hertz of the nominal
frequency and rotor angle deviations are limited for stability considerations
[23]. We use the set

X $\displaystyle=\\{x\in\mathbb{R}^{n}\mid-\bar{x}\leq V_{x}x\leq\bar{x}\\}$
(4)

to denote the constraints that the state $x$ must satisfy in real-time.

###  II-B Safety-critical control

Because of the presence of disturbances, it may not be possible for the system
state to always remain in X. Some states close to the boundary of X could be
pushed out by a disturbance no matter the control action, while for other
states in X, there may exist a control action such that no disturbance would
push the state outside of the prescribed region. This motivates the definition
of a _robust controlled-invariant set_.

######  Definition 1 (Robust controlled-invariant set (RCI) [4]).

An RCI is a set S for which there exists a feedback control policy
$u_{t}=\pi_{0}(x_{t})\in\textit{{U}}$ ensuring that all system trajectories
originating in S will remain in S for all time, under any disturbance sequence
$d_{t}\in\textit{{D}}$.

If S is contained in X, then $\pi_{0}$ is a safe policy. Often, the goal is to
find the policy that maximizes the size of S while being contained in X, since
it corresponds to making most of the acceptable states safe [9]. In general,
this is a difficult problem. Fortunately, if we restrict the policy to be
linear, there are many well-studied techniques that have been shown to be
successful at producing large safety sets [10, 30, 31, 32].

In this paper, we assume that S is a polytope described by $2r$ linear
inequalities, and that $\pi_{0}$ is a linear feedback control policy.
Specifically, we assume

$\displaystyle\textit{{S}}=\\{x\in\mathbb{R}^{n}\mid-\bar{s}\leq
V_{s}x\leq\bar{s}\\}\subseteq\textit{{X}}\text{ and}$ (5)
$\displaystyle\pi_{0}(x)=Kx$ (6)

where $V_{s}\in\mathbb{R}^{r\times n}$, $\bar{s}\in\mathbb{R}^{r}$, and
$Kx\in\textit{{U}}$ for all $x\in\textit{{S}}.$ For robustness, we choose the
largest set satisfying (5). The algorithm used for choosing
$(\textit{{S}},\pi_{0})$ is described in [30]. The algorithm uses a convex
relaxation to find an _approximately_ maximal RCI S and an associated $K$ as
the solution to an SDP. The objective of the SDP is to maximize the volume of
the largest inscribed ellipsoid inside S.

Of course, a linear policy that maximizes the size of S may not lead to
satisfactory control performance. The set S is chosen jointly with the policy
$\pi_{0}$, but there could be many policies (not necessarily linear) that keep
S robustly invariant. We want to optimize over this class of nonlinear
policies to improve the performance of the system. To explore the full range
of safe policies, we define the _safe action set at time $t$_ as

$\displaystyle\Omega(x_{t}):=\\{u_{t}\in\textit{{U}}\mid
x_{t+1}\in\textit{{S}},\ \forall\ d_{t}\in\textit{{D}}\\}$ (7)

where it is assumed that $x_{t}\in\textit{{S}}$. By induction, any policy that
chooses actions from $\Omega(x_{t})$ is a safe policy [4].

We define the set of safe policies with respect to the RCI S as

$\displaystyle\Pi:=\\{\pi:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\mid\pi(x_{t})\in\Omega(x_{t}),\
\forall\ x_{t}\in\textit{{S}}\\}.$ (8)

Given S, we search for a policy by optimizing over $\Pi$:

$\displaystyle\min_{\pi\in\Pi}\mathbb{E}\bigg{[}\frac{1}{T}\sum_{t=1}^{T}J(x_{t},u_{t})\bigg{]}$
(9a) $\displaystyle\text{subject to: }x_{0}\in\textit{{S}}$ (9b)
$\displaystyle x_{t+1}=Ax_{t}+Bu_{t}+Ed_{t}$ (9c) $\displaystyle
u_{t}=\pi(x_{t})$ (9d)

where $\mathbb{E}$ is the expectation with respect to randomness in initial
conditions and in the sequence of disturbances, and $J(x_{t},u_{t})$ is the
cost associated with occupying state $x_{t}$ and taking action $u_{t}$. To
estimate (9a), $d_{t}$ is sampled from D but treated as stochastic, so that
standard RL algorithms can be used to solve (9) [13]. This relaxation does not
require thorough sampling of D to preserve safety, since the constraint
$\pi\in\Pi$ imposes safety for _all_ possible disturbances
$d_{t}\in\textit{{D}}$.

One example of a cost function that can be used in (9) is the classical LQR
cost on state and control [33], but other non-quadratic cost functions can
also be used. For example, for sparsity-promoting controllers, we may set
$J(x_{t},u_{t})=x_{t}^{T}Qx_{t}+c\|u_{t}\|_{1},$ where $Q\succeq 0$ and
$c>~{}0$ [34].

##  III CONTROLLER DESIGN

Figure 1: Policy network architecture.

In this section, we describe how set-theoretic control techniques can be used
to create a safety guarantee for data-driven controllers without solving an
MPC or projection problem in real time. Since $d_{t}$ is unknown, data-driven
approaches for choosing $\pi$ are appropriate if safety guarantees can be
maintained. For control problems with continuous state and action spaces, one
class of RL algorithms involves parameterizing $\pi$ as a neural network or
other function approximator and using stochastic optimization to search over
the parameters of that function class for a (locally) optimal policy.

A common approach to safety-critical control with RL is to combine a model-
predictive controller with a neural network providing an action recommendation
or warm start [13, 14]. However, this makes it difficult to search over $\Pi$
efficiently and leads to control policies with higher computational overhead.
One optimization-free approach involves tracking the vertices of
$\Omega(x_{t})$ and using a neural network to choose convex weights on the
vertices of $\Omega(x_{t})$. However, this is only possible when S has
exceedingly simple geometry [17]. While it is difficult to constrain the
output of a neural network to arbitrary polytopes such as $\Omega(x_{t})$, it
is easy to constrain the output to $\mathbb{B}_{\infty}$, the $\infty$-norm
unit ball in $\mathbb{R}^{m}$, using activation functions like sigmoid or
hyperbolic tangent in the output layer. By establishing a correspondence
between points in $\mathbb{B}_{\infty}$ and points in $\Omega(x_{t}),$ we will
use neural network-based controllers to parameterize $\Pi$.

In particular, we construct a class of safe, differentiable, and closed-form
policies $\pi_{\theta}$, parameterized by $\theta$, that can approximate any
policy in $\Pi$. The policy first chooses a “virtual” action in
$\mathbb{B}_{\infty}$ using a neural network $\psi_{\theta}$. The policy then
uses a novel, closed-form, differentiable “safety filter” to equate
$\psi_{\theta}(x_{t})$ with an action in $\Omega(x_{t})$. Figure 1 illustrates
the way $\psi_{\theta}$, $\Omega$, and $\pi_{0}$ are interconnected using a
novel _gauge map_ in order to form the policy $\pi_{\theta}$. In order to
efficiently map between $\mathbb{B}_{\infty}$ and $\Omega(x_{t})$, we now
introduce the concepts of _C-sets_ and _gauge functions_.

######  Definition 2 (C-set [4]).

A _C-set_ is a set that is convex and compact and that contains the origin as
an interior point.

Any C-set can be used as a “measuring stick” in a way that generalizes the
notion of a vector norm [4]. In particular, the gauge function (or Minkowski
function) of a vector $v\in\mathbb{R}^{m}$ with respect to a C-set
$\textit{{Q}}\subset\mathbb{R}^{m}$ is given by

$\displaystyle\gamma_{\textit{{Q}}}(v)=\inf\\{\lambda\geq 0\mid
v\in\lambda\textit{{Q}}\\}.$ (10)

If Q is a polytopic C-set defined by $\\{w\in\mathbb{R}^{m}\mid F_{i}^{T}w\leq
g_{i},i=1,\ldots,r\\},$ then the gauge function is given by

$\displaystyle\gamma_{\textit{{Q}}}(v)=\max_{i}\Big{\\{}\frac{F_{i}^{T}v}{g_{i}}\Big{\\}},$
(11)

which is easy to compute since it is simply the maximum over $r$ elements.
Equation (11) is derived in Appendix -A. We will use (11) to construct a
closed-form, differentiable bijection between $\mathbb{B}_{\infty}$ and
$\Omega(x_{t})$.

###  III-A Gauge map

We will first show how to use the gauge function to construct a bijection from
$\mathbb{B}_{\infty}$ to any C-set Q, and will then generalize to the case
when Q does not contain the origin as an interior point. For any
$v\in\mathbb{B}_{\infty}$, we define the _gauge map from $\mathbb{B}_{\infty}$
to Q_ as

$\displaystyle
G(v|\textit{{Q}})=\frac{\|v\|_{\infty}}{\gamma_{\textit{{Q}}}(v)}\cdot v.$
(12)

######  Lemma 1.

For any C-set Q, the gauge map $G:\mathbb{B}_{\infty}\rightarrow\textit{{Q}}$
is a bijection. Specifically, $w=G(v|\textit{{Q}})$ if and only if $w$ and $v$
have the same direction and $\gamma_{\textit{{Q}}}(w)=\|v\|_{\infty}.$

The proof of Lemma 1 is provided in Appendix -C. By Lemma 1, choosing a point
in $\mathbb{B}_{\infty}$ is equivalent to choosing a point in Q. The action of
the gauge map is illustrated in Figure 2.

We cannot directly use the gauge map to convert between points in
$\mathbb{B}_{\infty}$ and points in $\Omega(x_{t})$, since $\Omega(x_{t})$ may
not contain the origin as an interior point. Instead, we must temporarily
“shift” $\Omega(x_{t})$ by one of its interior points, making it a C-set.
Lemma 2 provides an efficient way to achieve this.

Figure 2: Action of the gauge map from $\mathbb{B}_{\infty}$ to randomly
generated Q, with the $\frac{1}{4},\frac{1}{2},$ and $\frac{3}{4}$ level sets
of the respective gauge functions shown in white. For each point
$v\in\mathbb{B}_{\infty}$ and its image $w\in\textit{{Q}}$, $v$ and $w$ have
the same direction and $\gamma_{\textit{{Q}}}(w)=\|v\|_{\infty}$.

######  Lemma 2.

If $\pi_{0}(x)=Kx$ is a policy in $\Pi$, then for any $x_{t}$ in the interior
of S, the set $\hat{\Omega}_{t}:=[\Omega(x_{t})-Kx_{t}]$ is a C-set.

The proof of Lemma 2 is provided in Appendix -D. Figure 3 illustrates the way
the gauge map and Lemma 2 are used in the policy network as a safety filter,
by transforming the output of the policy network from $\mathbb{B}_{\infty}$ to
$\Omega(x_{t})$.

###  III-B Policy architecture

######  Theorem 1.

Assume the system dynamics and constraints are given by (1), (2) and (3), and
there exists a choice of $(\textit{{S}},\pi_{0})$ conforming to (5) and (6).
Let $\psi_{\theta}:\textit{{S}}\rightarrow\mathbb{B}_{\infty}$ be a neural
network parameterized by $\theta$. Then for any $x_{t}$ in the interior of S,
the policy

$\displaystyle\pi_{\theta}(x_{t}):=G(\psi_{\theta}(x_{t})|\hat{\Omega}_{t})+Kx_{t}$
(13)

has the following properties.

  1. 1.

$\pi_{\theta}$ is a safe policy.

  2. 2.

$\pi_{\theta}$ can be computed in closed form.

  3. 3.

$\pi_{\theta}$ is differentiable at $x_{t}$.

  4. 4.

$\pi_{\theta}$ can approximate any policy in $\Pi$.

The proof of Theorem 1 is provided in Appendix -E.

Figure 3: In the policy network, the gauge map is used to map virtual actions
to safe actions.

###  III-C Policy optimization through reinforcement learning

We parametrize the search over $\Pi$ using (13) with parameter $\theta$, and
we optimize $\theta$ using policy gradient RL algorithms. The policy gradient
theorem from reinforcement learning allows one to use past experience to
estimate the gradient of the cost function (9a) with respect to $\theta$ [35].
This is a standard approach for RL in continuous control tasks [36]. Policy
gradient methods require that it be possible to compute the gradient of
$\pi_{\theta}$ with respect to $\theta$. More specifically, $G$ must be
differentiable (Thm. 1, part 3) or else the safety filter would have to be
treated as an uncertain influence whose behavior must be estimated from data.

In addition to being differentiable, $\pi_{\theta}$ has two other noteworthy
attributes. First, under the optimal choice of $\theta$, the controller
$\pi_{\theta}$ performs no worse than $\pi_{0}$. This is because $\pi_{0}$ is
a feasible solution to (9), so the optimal solution to (9) can do no worse.
Second, unlike projection-based methods [15], the structure of $\pi_{\theta}$
facilitates exploration of the interior of the safe action set. This is
because smooth functions such as the sigmoid or hyperbolic tangent can be used
as activation functions in the output layer of $\psi_{\theta}$ to constrain
its output to $\mathbb{B}_{\infty}$. By tuning the steepness of the activation
function, it is possible to bias the output of $\psi_{\theta}$ towards or away
from the boundary of $\mathbb{B}_{\infty}$.

##  IV SIMULATIONS

###  IV-A Power system model

The main application considered in this paper is frequency control in power
systems. We consider a system with $N$ synchronous electric generators. The
standard linearized swing equation at generator $i$ is:

$\displaystyle\dot{\delta}_{i}$ $\displaystyle=\omega_{i}$ (14a)
$\displaystyle M_{i}\dot{\omega}_{i}$
$\displaystyle=-D_{i}\omega_{i}-\sum_{j=1}^{N}K_{ij}(\delta_{i}-\delta_{j})+\sum_{k=1}^{m}b_{ik}u_{k}-\sum_{l=1}^{p}e_{il}d_{l},$
(14b)

where $\delta_{i}$ is the rotor angle, $\omega_{i}$ is the frequency
deviation, and $M_{i}$ and $D_{i}$ are the inertia and damping coefficients of
generator $i$. The coefficients $K_{ij}$, $b_{ik}$, and $e_{il}$ are based on
generator and transmission line parameters taken from a modified IEEE 9-bus
test case, and are computed by solving the DC power flow equations. Thus, the
size of the coefficient measures the influence of each element on the dynamics
of generator $i$. The quantity $u_{k}$ represents controller $k$, an IBR such
as a battery energy storage system or wind turbine [37, 38], where the active
power injections can be controlled in response to a change in system
frequency. The feasible control set $\textit{{U}}\subset\mathbb{R}^{m}$
represents limits on power output for each of the $m$ IBRs.

The disturbance $d_{l}$ captures the uncertainties both in load and in
uncontrolled renewable resources. It is also possible to use $d$ to account
for parameteric uncertainties, linearization error associated with the
linearized swing equation dynamics, or error associated with the DC power flow
approximation, by adding virtual disturbances at every bus in the system [7,
39]. The disturbance set $\textit{{D}}\subset\mathbb{R}^{p}$ is conservatively
estimated from the capacity of the $p$ uncontrolled elements [40].

Discretizing the continuous-time system in (14) and assembling block
components gives a system in the form of (1). Let $\delta$ and
$\omega\in\mathbb{R}^{N}$ be vectors representing the rotor angles and
frequency deviations of all generators in the system, and let the system state
be represented by
$x=\begin{bmatrix}\delta&\omega\end{bmatrix}^{T}\in\mathbb{R}^{n}$ where
$n=2N.$ Using time step $\tau$, the discrete-time system matrices are given by

$\displaystyle A=\begin{bmatrix}I&\tau I\\\ -\tau M^{-1}K&I-\tau
M^{-1}D\end{bmatrix},$ $\displaystyle B=\begin{bmatrix}0\\\
M^{-1}\hat{B}\end{bmatrix},E=\begin{bmatrix}0\\\ M^{-1}\hat{E}\end{bmatrix}$

where $[M]_{ii}=M_{i}$, $[D]_{ii}=D_{i}$, $[K]_{ij}=K_{ij}$,
$[\hat{B}]_{ik}=b_{ik}$, and $[\hat{E}]_{il}=e_{il}$.

We simulate the proposed policy network architecture on a 9-bus power system
consisting of three synchronous electric generators, three controllable IBRs,
and three uncontrolled loads. The time step for discretization is 0.05
seconds. The load is modeled as autoregressive noise defined by

$\displaystyle d_{t+1}=\alpha d_{t}+(1-\alpha)\hat{d}$ (15)

where $\hat{d}\in\mathbb{R}^{p}$ is randomly generated from a uniform
distribution over $\textit{{D}},$ and $\alpha\in(0,1).$ The system is
illustrated in Figure 4.

Figure 4: Illustration of 9-bus power system model.

###  IV-B RL algorithm

To train the policy network, we used the Deep Deterministic Policy Gradient
(DDPG) algorithm [36], an algorithm well-suited for RL in continuous control
tasks. DDPG is an actor-critic algorithm, in which the “actor” or policy
chooses actions based on the state of the system, and the “critic” predicts
the value of state-action pairs in order to estimate the gradient of the cost
function (9a) with respect to $\theta$ (the “policy gradient”). In our
simulations, the cost was given by

$\displaystyle J(x_{t},u_{t})=x^{T}Qx+u^{T}Ru$ (16)

where $Q=\textbf{block diag}\\{1000I_{N},10I_{N}\\}$, $R=5I_{m}$, and $I_{N}$
is the identity matrix in $\mathbb{R}^{N\times N}$. The actor was given by
(13). The function $\psi_{\theta}$ was parameterized by a neural network with
two hidden layers of 256 nodes each, with ReLU activation functions in the
hidden layers. We use sigmoid functions in the last layer to limit the the
outputs to be within $[-1,1]$. The critic, or value network, had the same
hidden layers as $\psi_{\theta}$ but a linear output layer. We trained the
system for 200 episodes of 100 time steps each.

###  IV-C Benchmark comparisons

To demonstrate the advantages of the proposed policy architecture, we compare
against two benchmarks. The first is the linear controller $Kx$, chosen to
maximize the size of the associated RCI. Using the algorithm in [30], we
choose $(\textit{{S}},K)$ by solving the optimization problem

$\displaystyle\max_{\textit{{S}}\in\mathcal{S},K\in\mathbb{R}^{n\times m}}$
$\displaystyle\text{vol}(\textit{{S}})$ (17a) s.t. Invariance:
$\displaystyle(A+BK)\textit{{S}}\oplus E\textit{{D}}\subseteq\textit{{S}}$
(17b) Safety: $\displaystyle\textit{{S}}\subseteq\textit{{X}}$ (17c) Control
bounds: $\displaystyle K\textit{{S}}\subseteq\textit{{U}}$ (17d)

where $\oplus$ denotes Minkowski set addition and $\mathcal{S}$ is a class of
polytopes described by (5). Figure 5 displays the accumulated cost during a
number of test trajectories, showing that $\pi_{\theta}$ is a more cost-
effective controller than $Kx$ when the same S is used for each policy. This
makes sense, since the nonlinear policy is afforded additional flexibility in
balancing performance and robustness. Since $\pi_{\theta}$ and $Kx$ are both
policies in $\Pi$, the learned policy has the same safety guarantees as the
linear policy.

Figure 5: Accumulated cost over several test trajectories with a fixed initial
condition and randomly generated disturbance sequences, showing that the RCI
policy network achieves better performance than the safe linear feedback.

The second benchmark is a policy network that is trained using DDPG augmented
with a soft penalty on constraint violations, in order to incentivize
remaining in X. The policy network for this benchmark consists of two 256-node
hidden layers with ReLU activation, and hyperbolic tangent activation
functions in the output layer that clamp the output to the box-shaped set U.
The soft penalty is the total constraint violation, calculated as

$\displaystyle\lambda\|\max\\{V_{x}x_{t}-\bar{x},0\\}-\min\\{V_{x}x_{t}+\bar{x},0\\}\|_{1}$

where the $\max$ and $\min$ are taken elementwise, and $\lambda>0$.

In Fig. 6, we plot an example of the maximum angle deviation during training.
We place a hard constraint of 0.1 radians on this angle deviation. For the
policy given by (13), the trajectory always stays within this bound, by the
design of the the controller. For a policy trained with a soft penalty,
trajectories initially exit the safe set. With enough training, the
trajectories eventually satisfy the state constraints.

Figure 7 shows that safety in training does not imply safety in testing. The
policy network trained using a soft penalty can still result in constraint
violations, whereas the safe policy network does not. In some sense, this is
not unexpected. Only a limited number of disturbances can be seen during
training, and because of the nonlinearity of the neural network-based policy,
it is difficult to provide guarantees from the cost alone. In addition,
picking the right soft penalty parameter is nontrivial. If the penalty
$\lambda$ is too low, constraint satisfaction will not be incentivized, and if
$\lambda$ is too high, convergence issues may arise [41]. In our experiments,
we tuned $\lambda$ by hand to strike the middle ground, but even automatic,
dynamic tuning of $\lambda$ during training is not guaranteed to prevent
constraint violations in all cases [41].

Figure 6: Maximum angle deviation per training episode for the safe policy
network (blue) and the baseline policy network with soft penalty (orange). The
safe policy network guarantees safety during training, while soft penalties
eventually drive the baseline policy towards constraint satisfaction. Figure
7: Even though soft penalties succeed in driving policies to be safe during
training, they do not necessarily provide safety during testing. In this
example, a policy network that was safe during training (Fig. 6) still
exhibits constraint violations during testing. In contrast, starting from the
same initial conditions and subject to the same disturbance sequence, the
proposed safe policy network guarantees constraint satisfaction.

##  V CONCLUSIONS

In this paper, we propose an efficient approach to safety-critical, data-
driven control. The strategy relies on results from set-theoretic control and
convex analysis to provide provable guarantees of constraint satisfaction.
Importantly, the proposed policy chooses actions without solving an
optimization problem, opening the door to safety-critical control in
applications in which computational power may be a bottleneck. We apply the
proposed controller to a frequency regulation problem in power systems, but
the applications are much more wide-ranging. Future work includes
investigating robustness to changes in power system topology and extending the
proposed technique to decentralized control.

## ACKNOWLEDGMENTS

D. Tabas would like to thank Liyuan Zheng for guidance and Sarah H.Q. Li for
helpful discussions.

## References

  * [1] Y. Chen, J. Anderson, K. Kalsi, A. D. Ames, and S. H. Low, “Safety-critical control synthesis for network systems with control barrier functions and assume-guarantee contracts,” IEEE Transactions on Control of Network Systems, vol. 8, no. 1, pp. 487–499, 2021.
  * [2] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games,” IEEE Transactions on automatic control, vol. 50, no. 7, pp. 947–957, 2005\.
  * [3] J. Lygeros, “On reachability and minimum cost optimal control,” Automatica, vol. 40, no. 6, pp. 917–927, 2004.
  * [4] F. Blanchini and S. Miani, Set-theoretic methods in control.  Birkhauser, 2015.
  * [5] A. El-Guindy, Y. C. Chen, and M. Althoff, “Compositional transient stability analysis of power systems via the computation of reachable sets,” Proceedings of the American Control Conference, pp. 2536–2543, 2017.
  * [6] Y. Zhang, Y. Li, K. Tomsovic, S. Djouadi, and M. Yue, “Review on Set-Theoretic Methods for Safety Verification and Control of Power System,” IET Energy Systems Integration, pp. 2–12, 2020.
  * [7] A. El-Guindy, Y. C. Chen, and M. Althoff, “Compositional transient stability analysis of power systems via the computation of reachable sets,” Proceedings of the American Control Conference, pp. 2536–2543, 2017.
  * [8] A. El-Guindy, K. Schaab, B. Schurmann, O. Stursberg, and M. Althoff, “Formal LPV control for transient stability of power systems,” IEEE Power and Energy Society General Meeting, vol. 2018-Janua, pp. 1–5, 2018.
  * [9] J. N. Maidens, S. Kaynama, I. M. Mitchell, M. M. Oishi, and G. A. Dumont, “Lagrangian methods for approximating the viability kernel in high-dimensional systems,” Automatica, vol. 49, no. 7, pp. 2017–2029, 2013.
  * [10] C. Liu, F. Tahir, and I. M. Jaimoukha, “Full-complexity polytopic robust control invariant sets for uncertain linear discrete-time systems,” International Journal of Robust and Nonlinear Control, vol. 29, no. 11, pp. 3587–3605, 2019.
  * [11] F. Blanchini and A. Megretski, “Robust state feedback control of LTV systems: Nonlinear is better than linear,” IEEE Tranactions on Automatic Control, vol. 44, no. 4, pp. 2347–2352, 1999.
  * [12] T. Nguyen and F. Jabbari, “Disturbance Attenuation for Systems with Input Saturation: An LMI Approach,” Tech. Rep. 4, Department of Mechanical and Aerospace Engineering, University of California, Irvine, 1999.
  * [13] K. P. Wabersich and M. N. Zeilinger, “A predictive safety filter for learning-based control of constrained nonlinear dynamical systems,” Automatica, vol. 129, p. 109597, 2021.
  * [14] A. Aswani, H. Gonzalez, S. S. Sastry, and C. Tomlin, “Provably safe and robust learning-based model predictive control,” Automatica, vol. 49, no. 5, pp. 1216–1226, 2013.
  * [15] S. Gros, M. Zanon, and A. Bemporad, “Safe reinforcement learning via projection on a safe set: How to achieve optimality?,” IFAC-PapersOnLine, vol. 53, no. 2, pp. 8076–8081, 2020.
  * [16] R. Cheng, G. Orosz, R. M. Murray, and J. W. Burdick, “End-to-end safe reinforcement learning through barrier functions for safety-critical continuous control tasks,” in 33rd AAAI Conference on Artificial Intelligence, pp. 3387–3395, 2019.
  * [17] L. Zheng, Y. Shi, L. J. Ratliff, and B. Zhang, “Safe reinforcement learning of control-affine systems with vertex networks,” arXiv preprint: arXiv 2003.09488v1, 2020.
  * [18] J. Achiam, D. Held, A. Tamar, and P. Abbeel, “Constrained policy optimization,” in Proceedings of the 34th International Conference on Machine Learning (D. Precup and Y. W. Teh, eds.), vol. 70 of Proceedings of Machine Learning Research, pp. 22–31, PMLR, 06–11 Aug 2017.
  * [19] M. Yu, Z. Yang, M. Kolar, and Z. Wang, “Convergent policy optimization for safe reinforcement learning,” Advances in Neural Information Processing Systems, vol. 32, no. NeurIPS, 2019.
  * [20] W. Cui and B. Zhang, “Reinforcement Learning for Optimal Frequency Control: A Lyapunov Approach,” arXiv preprint: 2009.05654v3, 2021.
  * [21] W. Cui and B. Zhang, “Lyapunov-regularized reinforcement learning for power system transient stability,” IEEE Control Systems Letters, 2021.
  * [22] P. L. Donti, M. Roderick, M. Fazlyab, and J. Z. Kolter, “Enforcing robust control guarantees within neural network policies,” in International Conference on Learning Representations, pp. 1–26, 2021.
  * [23] P. Kundur, N. Balu, and M. Lauby, Power system stability and control.  New York: McGraw-Hill, 7 ed., 1994.
  * [24] C. Zhao, U. Topcu, N. Li, and S. Low, “Design and stability of load-side primary frequency control in power systems,” IEEE Transactions on Automatic Control, vol. 59, no. 5, pp. 1177–1189, 2014.
  * [25] A. Ademola-Idowu and B. Zhang, “Frequency stability using mpc-based inverter power control in low-inertia power systems,” IEEE Transactions on Power Systems, vol. 36, no. 2, pp. 1628–1637, 2020.
  * [26] X. Chen, G. Qu, Y. Tang, S. Low, and N. Li, “Reinforcement Learning for Decision-Making and Control in Power Systems: Tutorial, Review, and Vision,” arXiv preprint: arXiv 2102.01168, 2021.
  * [27] Z. Yan and Y. Xu, “Data-driven load frequency control for stochastic power systems: A deep reinforcement learning method with continuous action search,” IEEE Transactions on Power Systems, vol. 34, no. 2, pp. 1653–1656, 2018.
  * [28] A. Latif, S. S. Hussain, D. C. Das, and T. S. Ustun, “State-of-the-art of controllers and soft computing techniques for regulated load frequency management of single/multi-area traditional and renewable energy based power systems,” Applied Energy, vol. 266, p. 114858, 2020.
  * [29] S. Boyd and L. Vandenberghe, Convex Optimization.  Cambridge University Press, 2009.
  * [30] C. Liu and I. M. Jaimoukha, “The computation of full-complexity polytopic robust control invariant sets,” Proceedings of the IEEE Conference on Decision and Control, vol. 54rd IEEE, no. Cdc, pp. 6233–6238, 2015.
  * [31] F. Tahir, “Efficient computation of Robust Positively Invariant sets with linear state-feedback gain as a variable of optimization,” 7th International Conference on Electrical Engineering, Computing Science and Automatic Control, 2010, pp. 199–204, 2010.
  * [32] T. B. Blanco, M. Cannon, and B. De Moor, “On efficient computation of low-complexity controlled invariant sets for uncertain linear systems,” International Journal of Control, vol. 83, no. 7, pp. 1339–1346, 2010.
  * [33] J. P. Hespanha, Linear systems theory.  Princeton university press, 2018.
  * [34] F. Dörfler, M. R. Jovanović, M. Chertkov, and F. Bullo, “Sparsity-promoting optimal wide-area control of power networks,” IEEE Transactions on Power Systems, vol. 29, no. 5, pp. 2281–2291, 2014.
  * [35] D. P. Bertsekas, Reinforcement Learning and Optimal Control.  Athena Scientific, 2019.
  * [36] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” 4th International Conference on Learning Representations, ICLR 2016 - Conference Track Proceedings, 2016.
  * [37] B. Kroposki, B. Johnson, Y. Zhang, V. Gevorgian, P. Denholm, B. M. Hodge, and B. Hannegan, “Achieving a 100% Renewable Grid: Operating Electric Power Systems with Extremely High Levels of Variable Renewable Energy,” IEEE Power and Energy Magazine, vol. 15, no. 2, pp. 61–73, 2017.
  * [38] B. Xu, Y. Shi, D. S. Kirschen, and B. Zhang, “Optimal battery participation in frequency regulation markets,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6715–6725, 2018.
  * [39] J. Machowski, J. W. Bialek, and J. R. Bumby, Power System Dynamics: Stability and Control.  Wiley, 2 ed., 2008.
  * [40] Y. C. Chen, X. Jiang, and A. D. Domínguez-García, “Impact of power generation uncertainty on power system static performance,” in NAPS 2011 - 43rd North American Power Symposium, 2011.
  * [41] H. Yoo, V. M. Zavala, and J. H. Lee, “A Dynamic Penalty Function Approach for Constraint-Handling in Reinforcement Learning,” IFAC-PapersOnLine, vol. 54, no. 3, pp. 487–491, 2021.
  * [42] A. Baydin, A. A. Radul, B. A. Pearlmutter, and J. M. Siskind, “Automatic Differentiation in Machine Learning: a Survey,” Journal of Machine Learning Research, vol. 18, pp. 1–43, 2018.
  * [43] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural Networks, vol. 2, no. 5, pp. 359–366, 1989.

###  -A Derivation of (11)

Let $\textit{{Q}}=\\{x\in\mathbb{R}^{n}\mid Fx\leq g\\}$ be a C-set, where
$F\in\mathbb{R}^{r\times n},$ $g\in\mathbb{R}^{r}$, $F_{i}^{T}$ denotes the
$i$th row of $F$, and $g_{i}$ denotes the $i$th element of $g$. The gauge
function $\gamma_{\textit{{Q}}}(v)$ is computed as follows.

$\displaystyle\gamma_{\textit{{Q}}}(v)$ $\displaystyle=\inf\\{\lambda\geq
0\mid v\in\lambda\textit{{Q}}\\}$ (18) $\displaystyle=\inf\\{\lambda\geq
0\mid\frac{1}{\lambda}F_{i}^{T}v\leq g_{i},i=1,\ldots,r\\}$ (19)
$\displaystyle=\inf\\{\lambda\geq
0\mid\lambda\geq\frac{F_{i}^{T}v}{g_{i}},i=1,\ldots,r\\}$ (20)
$\displaystyle=\max\\{0,\max_{i}\\{\frac{F_{i}^{T}v}{g_{i}}\\}\\}.$ (21)

We now argue that $\max_{i}\\{\frac{F_{i}^{T}v}{g_{i}}\\}\geq 0.$ If
$F_{i}^{T}v<0$ for all $i$, then Q is unbounded in the direction of $v$ and Q
cannot be a C-set, a contradiction. Further, since
$0\in\textbf{int}(\textit{{Q}}),$ it must hold that $g_{i}>0$ for all $i$.
Therefore, there exists $i$ such that $\frac{F_{i}^{T}v}{g_{i}}\geq 0.$ ∎

###  -B Additional lemmas

The following lemma will be used in the proofs of Lemma 2 and Theorem 1.

######  Lemma 3.

Under the assumptions of Theorem 1, the safe action set $\Omega(x_{t})$ is a
polytope for all $x_{t}\in\textit{{S}}$.

###### Proof.

Starting from (1), (5), and (7), the safe action set is

$\displaystyle\Omega(x_{t})$
$\displaystyle=\\{u_{t}\in\textit{{U}}\mid-\bar{s}$ $\displaystyle\leq
V_{s}x_{t+1}\leq\bar{s},\ \forall\ d_{t}\in\textit{{D}}\\}$ (22)
$\displaystyle=\\{u_{t}\in\textit{{U}}\mid-\bar{s}_{i}$
$\displaystyle-\min_{d\in\textit{{D}}}V_{s}^{(i)T}Ed$ $\displaystyle\leq
V_{s}^{(i)T}(Ax_{t}+Bu_{t})$
$\displaystyle\leq\bar{s}_{i}-\max_{d\in\textit{{D}}}V_{s}^{(i)T}Ed,$
$\displaystyle\hskip 10.00002pt\forall\ i=1,\ldots,r\\}$ (23)

where $\bar{s}_{i}$ is the $i$th element of $\bar{s}$ and $V_{s}^{(i)T}$ is
the $i$th row of $V_{s}$. Since the $\min$ and $\max$ terms evaluate to
constant scalars for each $i$, and since $x_{t}$ is fixed, (23) is a set of
linear inequalities in $u_{t}$, making $\Omega(x_{t})$ a polytope [4]. ∎

###  -C Proof of Lemma 1

We will prove the more general case in which $\mathbb{B}_{\infty}$ is replaced
by any polytopic C-set. Let P and Q be two polytopic C-sets, and define the
gauge map from P to Q as
$G(v|\textit{{P}},\textit{{Q}})=\frac{\gamma_{\textit{{P}}}(v)}{\gamma_{\textit{{Q}}}(v)}\cdot
v$. We will prove that $G$ is a bijection from P to Q. The proof is then
completed by noting that $\gamma_{\mathbb{B}_{\infty}}$ is the same as the
$\infty$-norm.

To prove injectivity, we fix $v_{1},v_{2}\in\textit{{P}}$ and show that if
$G(v_{1}|\textit{{P}},\textit{{Q}})=G(v_{2}|\textit{{P}},\textit{{Q}})$ then
$v_{1}=v_{2}$. Assume
$G(v_{1}|\textit{{P}},\textit{{Q}})=G(v_{2}|\textit{{P}},\textit{{Q}})$. Then
$v_{1}$ and $v_{2}$ must be nonnegative scalar multiples of each other, i.e.
$v_{2}=\beta v_{1}$ for some $\beta\geq 0$. Making this substitution and
applying positive homogeneity of the gauge function [4] yields

$\displaystyle
G(v_{2}|\textit{{P}},\textit{{Q}})=\frac{\gamma_{\textit{{P}}}(v_{2})}{\gamma_{\textit{{Q}}}(v_{2})}v_{2}=\frac{\gamma_{\textit{{P}}}(v_{1})}{\gamma_{\textit{{Q}}}(v_{1})}v_{2}.$
(24)

Noting that
$G(v_{1}|\textit{{P}},\textit{{Q}})=\frac{\gamma_{\textit{{P}}}(v_{1})}{\gamma_{\textit{{Q}}}(v_{1})}v_{1}$,
we conclude that $v_{1}=~{}v_{2}$.

To prove surjectivity, fix $w\in\textit{{Q}}$. We must find $v\in\textit{{P}}$
such that $G(v|\textit{{P}},\textit{{Q}})=w$. Since P and Q are C-sets, each
set contains an open ball around the origin, thus P and Q each contain all
directions at sufficiently small magnitude. Choose $v$ in the same direction
as $w$ such that $\gamma_{\textit{{P}}}(v)=\gamma_{\textit{{Q}}}(w).$ Since
$w\in\textit{{Q}}$, $v\in\textit{{P}}.$ Then, we have

$\displaystyle G(v|\textit{{P}},\textit{{Q}})$
$\displaystyle=\frac{\gamma_{\textit{{P}}}(v)}{\gamma_{\textit{{Q}}}(v)}v$
(25)
$\displaystyle=\frac{\gamma_{\textit{{Q}}}(w)}{\gamma_{\textit{{Q}}}(v)}v.$
(26)

Since $v$ and $w$ are in the same direction,
$\frac{v}{\gamma_{\textit{{Q}}}(v)}=\frac{w}{\gamma_{\textit{{Q}}}(w)}$.
Making this substitution completes the proof. ∎

###  -D Proof of Lemma 2

Let int and bd denote the interior and boundary of a set, and rewrite (23) as
$\Omega(x)=\\{u\in\mathbb{R}^{m}\mid Hu\leq h,F(Ax+Bu)\leq g\\}$. Fix
$x\in\textbf{int}(\textit{{S}})$ and let $u^{*}=Kx$. By Lemma 3, $\Omega(x)$
is convex and compact. To fulfill the properties of a C-set, it remains to
show that $u^{*}\in\textbf{int}(\Omega(x)).$ Since $\pi_{0}\in\Pi$,
$u^{*}\in\Omega(x)$. Assume for the sake of contradiction that
$u^{*}\in\textbf{bd}(\Omega(x))$. Then either $F_{i}^{T}(A+BK)x=g_{i}$ or
$H_{j}^{T}Kx=h_{j}$ for some $i$ or $j$, where the subscript denotes a row
index. Suppose without loss of generality that the former holds, i.e.
$F_{i}^{T}(A+BK)x=g_{i}$ for some $i$. Since $x\in\textbf{int}(\textit{{S}}),$
there exists $\varepsilon\in(0,1)$ and
$\alpha=[1+\varepsilon\cdot\textbf{sign}(g_{i})]$ such that $y=\alpha x$ is
also in $\textbf{int}(\textit{{S}})$. The set $\Omega(y)$ is contained in the
halfspace $\\{u\mid F_{i}^{T}(Ay+Bu)\leq g_{i}\\}.$ Evaluating this inequality
with $u=Ky$, we have $F_{i}^{T}(A+BK)y=\alpha F_{i}^{T}(A+BK)x=\alpha
g_{i}>g_{i},$ thus $Ky\not\in\Omega(y)$ even though $y\in\textit{{S}}$,
contradicting the assumption that $\pi_{0}\in\Pi$. We conclude that
$u^{*}\not\in\textbf{bd}(\Omega(x)).$ Since $u^{*}\in\Omega(x),$ $u^{*}$ must
be an element of $\textbf{int}(\Omega(x)).$ ∎

###  -E Proof of Theorem 1

  1. 1.

It suffices to show that the gauge map from $\mathbb{B}_{\infty}$ to
$\hat{\Omega}_{t}$ is well-defined on $\textbf{int}(\textit{{S}}).$ This is a
direct result of Lemma 2.

  2. 2.

By Lemmas 2 and 3, $\hat{\Omega}_{t}$ is a polytopic C-set. By (11),
$\gamma_{\hat{\Omega}_{t}}$ (and $\pi_{\theta}$) can be computed in closed
form.

  3. 3.

A standard result from convex analysis shows that the subdifferential of (11)
is defined for all $v\in\mathbb{R}^{m}$, for any polytopic C-set Q [4]. If
$\psi_{\theta}$ is a neural network, automatic differentiation techniques can
be used to compute a subgradient of $\pi_{\theta}$ with respect to $\theta$
[42].

  4. 4.

This is due to the fact that $\psi_{\theta}$ is a universal function
approximator for functions from S to $\mathbb{B}_{\infty}$ [43]. By (13) and
Lemma 1, $\pi_{\theta}$ approximates any function in $\Pi$.


# Automated Model Design and Benchmarking of 3D Deep Learning Models for
COVID-19 Detection with Chest CT Scans

Xin He,1 Shihao Wang,1 Xiaowen Chu ,1 Shaohuai Shi,2 Jiangping Tang,3 Xin
Liu,3 Chenggang Yan,3
Jiyong Zhang,3∗ Guiguang Ding4
1Department of Computer Science, Hong Kong Baptist University, Hong Kong,
China
2Department of Computer Science and Engineering, Hong Kong University of
Science and Technology, Hong Kong, China
3School of Automation, Hangzhou Dianzi University, Hang Zhou, China
4School of Software, Tsinghua University, Beijing, China
Corresponding author

###### Abstract

The COVID-19 pandemic has spread globally for several months. Because its
transmissibility and high pathogenicity seriously threaten people’s lives, it
is crucial to accurately and quickly detect COVID-19 infection. Many recent
studies have shown that deep learning (DL) based solutions can help detect
COVID-19 based on chest CT scans. However, most existing work focuses on 2D
datasets, which may result in low quality models as the real CT scans are 3D
images. Besides, the reported results span a broad spectrum on different
datasets with a relatively unfair comparison. In this paper, we first use
three state-of-the-art 3D models (ResNet3D101, DenseNet3D121, and MC3_18) to
establish the baseline performance on the three publicly available chest CT
scan datasets. Then we propose a differentiable neural architecture search
(DNAS) framework to automatically search for the 3D DL models for 3D chest CT
scans classification with the Gumbel Softmax technique to improve the
searching efficiency. We further exploit the Class Activation Mapping (CAM)
technique on our models to provide the interpretability of the results. The
experimental results show that our automatically searched models (CovidNet3D)
outperform the baseline human-designed models on the three datasets with tens
of times smaller model size and higher accuracy. Furthermore, the results also
verify that CAM can be well applied in CovidNet3D for COVID-19 datasets to
provide interpretability for medical diagnosis. Code: https://github.com/HKBU-
HPML/CovidNet3D.

## Introduction

The Corona Virus Disease 2019 (COVID-19), pandemic is an ongoing pandemic
caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). The
SARS-CoV-2 virus can be easily spread among people via small droplets produced
by coughing, sneezing, and talking. COVID-19 is not only easily contagious but
also a severe threat to human lives. The COVID-19 infected patients usually
present pneumonia-like symptoms, such as fever, dry cough and dyspnea, and
gastrointestinal symptoms, followed by a severe acute respiratory infection.
The usual incubation period of COVID-19 ranges from one to 14 days. Many
COVID-19 patients do not even know that they have been infected without any
symptoms, which would easily cause delayed treatments and lead to a sudden
exacerbation of the condition. Therefore, a fast and accurate method of
diagnosing COVID-19 infection is crucial.

Currently, there are two commonly used methods for COVID-19 diagnosis. One is
viral testing, which uses real-time reverse transcription-prognosis chain
reaction (rRT-PCR) to detect viral RNA fragments. The other is making
diagnoses based on characteristic imaging features on chest X-rays or computed
tomography (CT) scan images. (?) conducted the effectiveness comparison
between the two diagnosis methods and concluded that chest CT has a faster
detection from the initial negative to positive than rRT-PCR. However, the
manual process of analyzing and diagnosing based on CT images highly relies on
professional knowledge and is time-consuming to analyze the features of the CT
images. Therefore, many recent studies have tried to use deep learning (DL)
methods to assist COVID-19 diagnosis with chest X-rays or CT scan images.

However, the reported accuracy of the existing DL-based COVID-19 detection
solutions spans a broad spectrum because they were evaluated on different
datasets, making it difficult to achieve a fair comparison. Besides, most
studies focus on 2D CT datasets (?; ?; ?). However, the real CT scan is
usually the 3D data. Thus it is necessary to use 3D models to classify 3D CT
scan data. To this end, we use three state-of-the-art (SOTA) 3D DL models to
establish the baseline performance on the three open-source 3D chest CT scan
datasets: CC-CCII111We find there are some errors and noises in the original
dataset (Version 1.0). Therefore we built our version based on it. (?),
MosMedData (?) and COVID-CTset (?). The details are shown in Table 2.

Table 1: Summary of the existing studies of DL-based methods for COVID-19
detection. NCP indicates the novel coronavirus pneumonia, Non-NCP includes CP
(common pneumonia) and Normal. $\ddagger$: the number of scans. *: the number
of patients.

Paper | Type |  Open- source? | Dataset Statistics | Class Statistics (#slices) |  Size of Test Set |  2D/3D model? |  Accuracy (%)
---|---|---|---|---|---|---|---
#patients | #scans | #slices | NCP | Non-NCP
CP | Normal
(?) | X-ray(2D) | Yes | - | - | 5,941 | 68 | 4,290 | 1,583 | 1,188 | 2D | 88.39
(?) | X-ray(2D) | Yes | - | - | 1,531 | 100 | 1,431 | 764 | 2D | -
(?) | X-ray(2D) | Yes | - | - | 100 | 50 | 50 | 20 | 2D | 98.00
(?) | CT(2D) | No | - | - | 133 | 68 | 65 | 26 | 2D | 93.20
(?) | CT(2D) | No | 194 | - | 1,020 | 510 | 510 | 102 | 2D | 99.63
(?) | CT(2D) | Yes | - | - | 425 | 178 | 247 | 45 | 2D | 98.78
(?) | CT(2D) | Yes | 143 | - | 746 | 349 | 397 | 186 | 2D | 86.00
(?) | CT(2D) | Yes | - | - | 746 | 349 | 397 | 105 | 2D | 87.60
(?) | CT(3D) | Yes | 377 | 526 | 12,058 | 244$\ddagger$ | 282$\ddagger$ | 124$\ddagger$ | 2D | -
(?) | CT(3D) | No | 542 | 630 | - | 313* | 229* | 131* | 3D | 90.10
(?) | CT(3D) | No | 3,322 | 4,356 | - | 1,296$\ddagger$ | 1,735$\ddagger$ | 1,325$\ddagger$ | 427$\ddagger$ | 3D | -
(?) | CT(3D) | Yes | 1,110 | 1,110 | 46,411 | 856$\ddagger$ | 254$\ddagger$ | 331$\ddagger$ | 3D | -
(?) | CT(3D) | Yes | 2,778 | 4,356 | 444,034 | 1,578$\ddagger$ | 1,614$\ddagger$ | 1,164$\ddagger$ | 389$\ddagger$ | 3D | 92.49

In addition, designing a high-quality model for the specific medical image
dataset is a time-consuming task and requires much expertise, which hinders
the development of DL technology in the medical field. Recently, neural
architecture search (NAS) has become a prevalent topic, as it can efficiently
discover high-quality DL models automatically. Many studies have used the NAS
technique to image classification and object detection tasks (?; ?; ?; ?). In
this paper, we present a differentiable neural architecture search (DNAS)
method combined with the Gumbel Softmax (?) technique to search neural
architectures on three chest CT datasets: Clean-CC-CCII (?), MosMedData (?),
and COVID-CTset (?). We represent the search space by a supernet. Using the
Gumbel Softmax technique, we can optimize only one subnetwork of the supernet
at a time; therefore, the searching efficiency can be significantly improved,
and the search stage can be finished in about 2 hours using 4 Nvidia Tesla
V100 GPUs. We name the model searched by DNAS as CovidNet3D. The experimental
results show that CovidNet3D can achieve comparable results to human-designed
SOTA models with a smaller size. Furthermore, medical diagnoses generally
require interpretability of the decision, so we apply Class Activation Mapping
(CAM) (?) techniques to provide interpretability for our CovidNet3D models. In
summary, our contributions are summarized as follows:

  * •

We use three manually designed 3D models to establish the baseline performance
on the three open-source COVID-19 chest CT scan datasets.

  * •

To the best of our knowledge, we are the first to apply the NAS technique to
search for 3D DL models for COVID-19 chest CT scan datasets. Our DNAS
framework can efficiently discover competitive neural architectures that
outperform the baseline models on the three CT datasets.

  * •

We use the Class Activation Mapping (CAM) (?) algorithm to add the
interpretability of our DNAS-designed models, which can help doctors quickly
locate the discriminative lesion areas on the CT scan images.

## Related Work

In recent years, DL techniques have been proven to be effective in diagnosing
diseases with X-ray and CT images (?). To enable DL techniques to be applied
in helping the detection of COVID-19, an increasing number of publicly
available COVID-19 datasets have been proposed, as shown in Table 1.

### Publicly-available Datasets of COVID-19

We separate the publicly available datasets into two different categories: the
pre-pandemic datasets and the post-pandemic datasets which mainly differ in
quality and quantity.

#### Pre-pandemic Datasets

In the pre-pandemic period, gathering datasets for COVID-19 is a tough job as
there is no enough data for collection. Most datasets in this period were
gathered from medical papers or uploaded by the public. IEEE8023 Coivd-
chestxray-dataset (?) is a dataset of COVID-19 cases with chest X-ray and CT
images collected from public sources. But its quality has no guarantee since
the images are not verified by medical experts. Covid-ct-dataset (?) is
another CT dataset of COVID-19, mainly composed of CT images extracted from
COVID-19 research papers, and its quality is low. The dataset only contains 2D
information because each patient has only one to several CT images instead of
a complete 3D scan volume.

#### Post-pandemic Datasets

During the pandemic, the number of confirmed cases of COVID-19 has been rising
rapidly, which brings many high-quality COVID-19 chest CT scan datasets, such
as CC-CCII (?) and COVID-CTset (?). Some of them have annotations by doctors,
e.g., COVID-19-CT-Seg-Dataset (?) and MosMedData (?). The three datasets we
use in this work are all from this category.

### DL-based methods for COVID-19 detection

Much research is conducted on CT images, but the 3D information of CT images
is under-explored, such as the work by (?; ?; ?). These work mainly propose 2D
DL models for COVID-19 detection. (?) benchmarks ten 2D CNNs and compares
their performance in classifying 2D CT images on their private dataset with
102 testing images. On the other hand, the studies in utilizing 3D CT images
are relatively rare, mainly due to the lack of 3D COVID-19 CT scan datasets.
(?; ?) propose 3D CNNs with their private 3D CT datasets. There are also some
other studies conducted on X-ray images. For example, (?) proposes three 2D DL
models for COVID-19 detection. (?) introduces a deep anomaly detection model
for fast and reliable screening. (?) investigates the estimation of
uncertainty and interpretability by Bayesian CNN on the X-ray images. (?) uses
both X-ray images and CT images to do segmentation and detection.

### Neural Architecture Search

In recent years, NAS has created many SOTA results by automatically searching
for neural architectures for many tasks (?; ?). (?; ?) first propose to use
reinforcement learning (RL) to search for neural architectures and achieves
comparable results to SOTA human-designed models. Since then, several types of
NAS methods have been proposed, such as evolutionary algorithm (EA) (?),
surrogate model-based optimization (SMBO) (?), and gradient descent (GD) based
methods (?; ?). (?; ?) combine the GD-based method and the Gumbel Softmax (?)
technique to further improve the searching efficiency.

Due to the success of NAS in natural image recognition (such as ImageNet (?)),
researchers also try to extend it to the medical datasets, such as Magnetic
resonance imaging (MRI) segmentation (?). (?) uses five public datasets,
MESSIDOR, OCT images, HAM 10000, Paediatric images, and CXR images, to search
for and train models by Google Cloud AutoML platform. Their experimental
results demonstrate that AutoML can generate competitive classifiers compared
to manually designed DL models. But to the best of our knowledge, there is no
study applying the NAS technique to search for 3D DL models for COVID-19 chest
CT scan datasets. To this end, we exploit the NAS technique to the three open-
source COVID-19 chest CT scan datasets and successfully discover high-quality
3D models that achieve comparable performance with the human-designed SOTA 3D
models.

## Method

In this section, we first describe our search space for 3D CT scans
classification models. Then, we introduce the differentiable neural
architecture search (DNAS) method combined with the Gumbel Softmax technique
(?; ?).

Figure 1: The overview of our search space. The model is generated by stacking
a predefined number of cells. Each cell contains different number of blocks,
and the block of different cells is different and needs to be searched. Conv3d
1$\times$1$\times$1 denotes 3D convolution with 1$\times$1$\times$1 kernel
size, Dwise3d denotes 3D depthwise convolution, BN3d denotes 3D batch norm,
$D\times H\times W\times F$ denotes tensor shape (depth, height, width,
channel), and MBConv denotes mobile inverted bottleneck convolution.

### Search Space

There are two critical points to be considered before designing the search
space. One is that all datasets we use are composed of 3D CT scans; therefore,
the searched model should be good at extracting the information from three-
dimensional spatial data. The other is that the model should be lightweight,
as the time required to process 3D data is much longer than 2D image data.

Although the cell-based search space (?; ?) is one of the most commonly used
search space, it has several problems: 1) the final model is built by stacking
the same cells, which precludes the layer diversity; 2) many searched cells
are very complicated and fragmented and are therefore inefficient for
inference. MobileNetV2 (?) is a lightweight model manually designed for mobile
and embedded devices for efficient inference. Several NAS studies (?; ?) have
successfully used the layer modules (?) including inverted residuals and
linear bottlenecks to search for neural architectures and achieved SOTA
results on the 2D image datasets. Therefore, we use MobileNetV2 as a reference
to design our 3D search space.

As shown in Fig. 1, we represent the search space by a supernet, which
consists of the stem layer, a fixed number of cells, and a linear layer. The
stem layer performs convolutional operations, and the last linear layer
follows behind a 3D global average pooling operation (?). Each cell is
composed of several blocks. The structures of all blocks need to be searched.
In different cells, the number of channels and the number of blocks are
different and hand-picked empirically. By default, all blocks have a stride of
1. However, if a cell’s input/output resolutions are different, then its first
block has a stride of 2. The blocks within the same cell have the same number
of input/output channels. Inspired by MobileNetV2 (?), each block is a MBConv-
similar module (see Fig. 1). It consists of three sub-modules: 1) a point-wise
($1\times 1\times 1$) convolution; 2) a 3D depthwise convolution with $K\times
K\times K$ kernel size, where $K$ is a searchable parameter; 3) another point-
wise ($1\times 1\times 1$) convolution. All convolutional operations are
followed by a 3D batch normalization and a ReLU6 activation function (?),
which is denoted by Conv3D-BN3D-ReLU6, and the last convolution has no ReLU6
activation. Another searchable parameter is the expansion ratio $e$, which
controls the ratio between the size of the input bottleneck and the inner
size. For example, $5\times 5\times 5$ MBConv6 denotes that the kernel size of
MBConv is $5\times 5\times 5$, and the expansion ratio is 6.

In our experiments, the search space is a fixed macro-architecture supernet
consisting of 6 cells, where each has 4 blocks, but the last cell only has 1
block. We empirically collect the following set of candidate operations:

  * •

$3\times 3\times 3$ MBConv3

  * •

$3\times 3\times 3$ MBConv4

  * •

$3\times 3\times 3$ MBConv6

  * •

$5\times 5\times 5$ MBConv3

  * •

$5\times 5\times 5$ MBConv4

  * •

$7\times 7\times 7$ MBConv3

  * •

$7\times 7\times 7$ MBConv4

  * •

Skip connection

Therefore, it contains $8^{21}\approx 9.2\times 10^{18}$ possible
architectures. Finding an optimal architecture from such a huge search space
is a stupendous task. We will introduce our search strategy in the following.

### Differentiable NAS with Gumbel Softmax

According to (?), gradient descent (GD) based NAS is an efficient method, and
many studies use it to find competitive models with much shorter time and less
computational resources (?; ?) than other NAS methods. Hence, in this paper,
we use the GD-based method and combine it with the Gumbel Softmax (?)
technique to discover models for COVID-19 detection.

#### Preliminary: DARTS

DARTS (?) was one of the first studies to use GD-based method to search for
neural architectures. Each cell is defined as a directed acyclic graph (DAG)
of $N$ nodes, where each node is a network layer, and each edge between node
$i$ and node $j$ indicates a candidate operation (i.e., block structure) that
is selected from the predefined operation space $\mathcal{O}$. To make the
search space continuous, DARTS (?) uses Softmax over all possible operations
to relax the categorical choice of a particular operation, i.e.,

$\begin{array}[]{cl}\overline{o}_{i,j}(x)&=\sum_{k=0}^{K}P_{k}o^{k}(x)\\\
s.t.\,\,\,P_{k}&=\frac{\exp\left(\alpha^{k}_{i,j}\right)}{\sum_{l=0}^{K}\exp\left(\alpha^{l}_{i,j}\right)}\end{array},$
(1)

where $o^{k}$ indicates the $k$-th candidate operation performed on input $x$,
$\alpha^{k}_{i,j}$ indicates the weight for the operation $o^{k}$ between a
pair of nodes ($i,j$), and $K$ is the number of predefined candidate
operations. The training and the validation loss are denoted by
$\mathcal{L}_{train}$ and $\mathcal{L}_{val}$, respectively. Therefore, the
task of searching for architectures is transformed into a bilevel optimization
problem of neural architecture $\alpha$ and the weights $\omega_{\alpha}$ of
the architecture:

$\begin{array}[]{cl}{\min_{\alpha}}&{\mathcal{L}_{val}\left(\omega_{\alpha}^{*},\alpha\right)}\\\
{\text{ s.t.
}}&{\omega_{\alpha}^{*}=\operatorname{argmin}_{\omega_{\alpha}}\,\mathcal{L}_{train}(\omega_{\alpha},\alpha)}\end{array}$
(2)

#### Differentiable Model Sampling by Gumbel Softmax

In DARTS, as Fig. 2 (left) shows, the output of each node is the weighted
average of the mixed operation during the whole search stage. It causes a
linear increase in the requirements of computational resources with the number
of candidate operations. To alleviate this problem, we follow the same idea as
(?). Specifically, for each layer, only one operation is sampled and executed
with the sampling probability distribution $P_{\alpha}$ defined in Equation 1.
For example, the probability of being sampled for the three operations in Fig.
2 (left) is 0.1, 0.2, and 0.7, respectively, but only one operation will be
sampled at a time. Therefore, the sampling distribution $P_{\alpha}$ of all
layers is encoded into a one-hot random distribution $Z$, e.g.,
$P_{\alpha}=[0.1,0.2,0.7]\rightarrow Z=[0,0,1]$.

Figure 2: The comparison between two GD-based methods. (Left) Applying a
mixture of all candidate operations, each with different weight. (Right) Only
one operation is sampled at a time. (best viewed in color) Figure 3: The
pipeline of training 3D deep learning models. All CT scans need to be pre-
processed by the slice sampling strategy to make sure that each scan contains
the same number of slices. The input size of network is $bs\times 1\times
d\times h\times w$, where $bs$ is batch size, $d$ is the number of slices, $h$
and $w$ indicate the height and width, respectively.

However, each operation is sampled from a discrete probability distribution
$Z$, so we cannot back-propagate gradients through $Z$ to $\alpha$. To enable
back-propagation, we use a reparameterization trick named Gumbel Softmax (?),
which can be formulated by

$\begin{array}[]{cl}Z_{k}=\frac{\exp\left(\left(\log{\alpha}_{i,j}^{k}+{G}_{i,j}^{k}\right)/\tau\right)}{\sum_{l=0}^{K}\exp\left(\left(\log{\alpha}_{i,j}^{l}+{G}_{i,j}^{l}\right)/\tau\right)}\end{array},$
(3)

where ${G}_{i,j}^{k}=-log(-log(u_{i,j}^{k}))$ is the $k$-th Gumbel sample,
$u_{i,j}^{k}$ is a uniform random variable, and $\tau$ is the softmax
temperature. When $\tau\rightarrow\infty$, the possibility distribution of all
operations between each pair of nodes approximates to the one-hot
distribution. To be noticed, we perform $argmax$ function on Equation 3 during
the forward process but return the gradients according to the Equation 3
during the backward process.

### Class Activation Mapping Algorithm

As mentioned above, the last linear layer follows behind a 3D global average
pooling layer, which enables us to utilize class activation mapping (CAM)
algorithm to generate 3D activation maps for our model. CAM exploits the
global average pooling layer to calculate get the activation map $M_{c}$ for
class $c$, where each spacial element is given by

$M_{c}(x,y,z)=\sum_{k}w_{k}^{c}f_{k}(x,y,z)$ (4)

where in a given image, $f_{k}(x,y,z)$ is the activation of unit $k$ at the
last convolutional layer before global average pooling layer at spatial
location $(x,y,z)$, $w_{k}^{c}$ is the corresponding linear layer weight of
class $c$ for unit $k$. After getting the class activation map, we can simply
upsample it to the size of the input scan images to visualize and identify the
regions most relevant to the specific class.

## Experiments

### Datasets and Pre-processing

In this paper, we use three publicly available datasets: CC-CCII (?),
MosMedData (?) and COVID-CTset (?). The three datasets are all chest CT
volumes. However, since the data format varies from the three datasets, it is
necessary to pre-process each dataset to make them follow a unified way of
reading data.

The original CC-CCII dataset contains a total number 617,775 slices of 6,752
CT scans from 4,154 patients, but it has five main problems (i.e., damaged
data, non-unified data type, repeated and noisy slices, disordered slices, and
non-segmented slices) that would have high negative impacts on the model
performance. To solve these problems, we manually remove the damaged, repeated
and noisy data. Then we segment the lung part for the unsegmented slice image
and convert the whole dataset to PNG format. After addressing the above
problems, we build a clean CC-CCII dataset named Clean-CC-CCII, which consists
of 340,190 slices of 3,993 scans from 2,698 patients.

#### Scan images construction

Each CT scan contains a different number of slices, but DL models require the
same dimensional inputs. To this end, we propose two slice sampling
algorithms: random sampling and symmetrical sampling. Specifically, the random
sampling strategy is applied to the training set, which can be regarded as the
data augmentation, while the symmetrical sampling strategy is performed on the
test set to avoid introducing randomness into the testing results. The
symmetrical sampling strategy refers to sampling from the middle to both sides
at equal intervals. The relative order between slices remains the same before
and after sampling.

Table 2: The statistics of the three CT scan datasets. Dataset [Format] | Classes | #Patients | #Scans
---|---|---|---
|  | Train | Test | Train | Test
Clean- CC-CCII [PNG] | NCP | 726 | 190 | 1213 | 302
CP | 778 | 186 | 1210 | 303
| Normal | 660 | 158 | 772 | 193
| Total | 2164 | 534 | 3195 | 798
MosMedData [PNG] | NCP | 601 | 255 | 601 | 255
Normal | 178 | 76 | 178 | 76
| Total | 779 | 331 | 779 | 331
COVID-CTset [16bit TIFF] | NCP | 80 | 15 | 202 | 42
Normal | 200 | 82 | 200 | 82
| Total | 280 | 97 | 402 | 124

### Benchmarking

We use three manually-designed 3D neural architectures as the baseline
methods: DenseNet3D121 (?), ResNet3D101 (?), and MC3_18 (?). As shown in Fig.
3, after building the scan images by the sampling algorithm, we further apply
transformations to scans, including resize, center-crop, and normalization.
Besides, for the training set, we also perform a random flip operation in the
horizontal or vertical direction. The other implementation details are as
follows: we use the Adam (?) optimizer and the weight decay of 5e-4. We start
the learning rate of 0.001 and anneal it down to 1e-5. All baseline models are
trained for 200 epochs.

### DNAS for CT Scan Classification

To verify the efficiency of the method, we apply the DNAS method combined with
the Gumbel Softmax technique to search for neural architectures on the three
datasets. The pipeline of our method is shown in Fig. 4, which contains two
sequential stages: search stage and evaluation stage.

#### Search stage

In our experiments, the supernet consists of 6 cells with the number of blocks
of $[4,4,4,4,4,1]$. Besides, the blocks within the same cell have the same
number of channels. Here, we test two settings: small-scale and large-scale,
where the number of channels of blocks in the 6 cells is
$[24,40,80,96,192,320]$ and $[32,64,128,256,512,1024]$, respectively. We name
the models searched under the two settings as CovidNet3D-S and CovidNet3D-L,
respectively. The stem block is a Conv3D-BN3D-ReLU6 sequential module with the
number of output channels fixed to 32.

To improve searching efficiency, we set the input resolution to 64$\times$64,
and the number of slices in a scan to 16. We implement three independent
search experiments for the three datasets. During the search stage, we split
the training set into the training set $\mathcal{D_{T}}$ and the validation
set $\mathcal{D_{V}}$. In each step, we first use $\mathcal{D_{V}}$ to update
the architecture parameters $\alpha$, and then use the training set to update
the sampled architecture weights $\omega_{\alpha}$. Besides, the architecture
parameter $\alpha$ is optimized by the Adam (?) optimizer, and the
architecture weights are optimized with the SGD optimizer with a momentum of
3e-4. The initial learning rate for both optimizers is 0.001. Each experiment
is conducted on four Nvidia Tesla V100 GPUs (the 32GB PCIe version) and it can
be finished in about 2 hours. After each epoch, we save the sampled
architecture and its performance (e.g., accuracy). Therefore, we generate 100
neural architectures for each experiment after the search stage.

Figure 4: The pipeline of DNAS consists of two stages: the search and the
evaluation stage.

#### Evaluation stage

As Fig. 4 shows, the search stages records the performance of the sampled
architectures. In the evaluation stage, we select top-10 architectures and
training these architectures with the training set for several batches, then
the best-performing architecture will be retrained for 200 epochs with the
full training set, and then evaluated on the test set. We set different input
resolutions for three datasets to evaluate the generalization of searched
architectures. Besides, since the number of slices contained in CT scans of
different datasets is different, we set the intermediate value for each
dataset, shown in Table 3. Each evaluation experiment uses the same settings
as follows: we use the Adam (?) optimizer with an initial learning rate of
0.001. The cosine annealing scheduler (?) is applied to adjust the learning
rate. We use Cross-entropy as the loss function.

Table 3: The experimental results of manually designed models and DNAS-
designed models.

Dataset | Model |  Model size (MB) | Input size | #Slices |  Accuracy (%) |  Precision (%) |  Sensitivity (%) | F1-score
---|---|---|---|---|---|---|---|---
Clean-CC-CCII | ResNet3D101 | 325.21 | 128$\times$128 | 32 | 85.54 | 89.62 | 77.15 | 0.8292
DenseNet3D121 | 43.06 | 128$\times$128 | 32 | 87.02 | 88.97 | 82.78 | 0.8576
MC3_18 | 43.84 | 128$\times$128 | 32 | 86.16 | 87.11 | 82.78 | 0.8489
CovidNet3D-S | 11.48 | 128$\times$128 | 32 | 88.55 | 88.78 | 91.72 | 0.9023
CovidNet3D-L | 53.26 | 128$\times$128 | 32 | 88.69 | 90.48 | 88.08 | 0.8926
MosMedData | ResNet3D101 | 325.21 | 256$\times$256 | 40 | 81.82 | 81.31 | 97.25 | 0.8857
DenseNet3D121 | 43.06 | 256$\times$256 | 40 | 79.55 | 84.23 | 92.16 | 0.8801
MC3_18 | 43.84 | 256$\times$256 | 40 | 80.4 | 79.43 | 98.43 | 0.8792
CovidNet3D-S | 12.48 | 256$\times$256 | 40 | 81.17 | 78.82 | 99.22 | 0.8785
CovidNet3D-L | 60.39 | 256$\times$256 | 40 | 82.29 | 79.50 | 98.82 | 0.8811
Covid-CTset | ResNet3D101 | 325.21 | 512$\times$512 | 32 | 93.87 | 92.34 | 95.54 | 0.9392
DenseNet3D121 | 43.06 | 512$\times$512 | 32 | 91.91 | 92.57 | 92.57 | 0.9257
MC3_18 | 43.84 | 512$\times$512 | 32 | 92.57 | 90.95 | 94.55 | 0.9272
CovidNet3D-S | 8.36 | 512$\times$512 | 32 | 94.27 | 92.68 | 90.48 | 0.9157
CovidNet3D-L | 62.82 | 512$\times$512 | 32 | 96.88 | 97.50 | 92.86 | 0.9512

## Results and Analysis

### Evaluation Metrics

Our experiment results are summarized in Table 3. We compare our searched
models with SOTA efficient models. We use several commonly used evaluation
metrics to compare the model performance, as follows:

$Precision=\frac{TP}{TP+FP}$ (5)

$Sensitivity(Recall)=\frac{TP}{TP+FN}$ (6)

$F1-score=\frac{2\times(precision\times recall)}{precision+recall}$ (7)

$Accuracy=\frac{TN+TP}{TN+TP+FN+FP}$ (8)

To be noticed, the positive and negative cases are assigned to the COVID-19
class and the non-COVID-19 class, respectively. Specifically, $TP$ and $TN$
indicate the number of correctly classified COVID-19 and non-COVID-19 scans,
respectively. $FP$ and $FN$ indicate the number of wrongly classified COVID-19
and non-COVID-19 scans, respectively. For the Clean-CC-CCII dataset, the non-
COVID-19 class includes both normal and common pneumonia. The accuracy is the
micro-averaging value for all test data to evaluate the overall performance of
the model. Besides, we also take the model size as an evaluation metric to
compare the model efficiency.

### Results on the Three CT Datasets

Table 3 divides the results according to the datasets. We can see that our
searched CovidNet3D models outperform all baseline models on the three
datasets in terms of accuracy. Specifically, CovidNet3D-L models achieve the
highest accuracy of the three datasets. Besides, all CovidNet3D-S models are
with much smaller sizes than the baseline models, but they can also achieve
similar or even better results. For example, CovidNet3D-S (8.36 MB) achieves
94.27% accuracy on Covid-CTset, which is 41$\times$ smaller than ResNet3D101
(325.21 MB) with 0.4% higher accuracy. In summary, the results demonstrate
that our DNAS method can discover well-performing models without inconsistency
on network size, input size or scan depth (the number of slices).

We can also see that the performance of both baseline models and our
CovidNet3D on the MosMedData dataset is not as good as that on the other two
datasets. There are two possible reasons. One is that the MosMedData
datasets’s original data format is NIfTI, but all our models do not converge
when trained with NIfTI files; therefore we convert NIfTI to Portable Network
Graphics (PNG) format, and this process would loss information of the input
files. The other possible reason is that the MosMedData dataset is imbalanced
(shown in Table 2), which increases the difficulty of model training.

We also find that the random seed greatly influences on the training of the
searched CovidNet3D model through experiments. In other words, the results
obtained by using different seeds for the same model would differ
significantly. Hence, how to improve the robustness of NAS-based models is
worthy for further exploring.

Figure 5: The class activation mappings generated on CovidNet3D on a chest CT
scan of the Clean-CC-CCII dataset. Regions colored in red and brighter has
more impact on model’s decision to the class of COVID-19 while blue and darker
region has less.

### Interpretability

Although our model achieves promising result in detecting COVID-19 in CT
images, classification result itself does not help clinical diagnosis without
proving the inner mechanism which leads to the final decision makes sense. To
inspect our CovidNet3D model’s inner mechanism, we apply Class Activation
Mapping (CAM) (?) on it.

CAM is an algorithm that can visualize the discriminative lesion regions that
the model focuses on. In Fig. 5, we apply CAM on each slice of a whole 3D CT
scan volume from Clean-CC-CCII dataset. Regions appear red and brighter have a
larger impact on the model’s decision to classify it to COVID-19. From the
perspective of the scan volume, we can see that some slices have more impacts
on the model’s decision than the others. In terms of a single slice, the areas
that CovidNet3D focuses on has ground-glass opacity, which is proved a
distinctive feature of CT images of COVID-19 Chest CT images (?). CAM enables
the interpretability of our searched models (CovidNet3D), helping doctors
quickly locate the discriminative lesion areas.

## Conclusion

In this work, we introduce the differentiable neural architecture (DNAS)
framework combined with the Gumbel Softmax technique to search for 3D models
on three open-source COVID-19 CT scan datasets. The results show that
CovidNet3D, a family of models discovered by DNAS can achieve comparable
results to the baseline 3D models with smaller size, which demonstrates that
NAS is a powerful tool for assisting in COVID-19 detection. In the future, we
will apply NAS to the task of 3D medical image segmentation to locate the
lesion areas in a more fine-grained manner.

## Acknowledgments

The research was supported by the grant RMGS2019_1_23 from Hong Kong Research
Matching Grant Scheme . We would like to thank the anonymous reviewers for
their valuable comments.

## References

  * [Ai et al. 2020] Ai, T.; Yang, Z.; Hou, H.; Zhan, C.; Chen, C.; Lv, W.; Tao, Q.; Sun, Z.; and Xia, L.  2020\.  Correlation of chest ct and rt-pcr testing in coronavirus disease 2019 (covid-19) in china: a report of 1014 cases.  Radiology 200642\.
  * [Alom et al. 2020] Alom, M. Z.; Rahman, M. M. S.; Nasrin, M. S.; Taha, T. M.; and Asari, V. K.  2020\.  Covid_MTNet: Covid-19 detection with multi-task deep learning approaches.  arXiv preprint arXiv:2004.03747.
  * [Ardakani et al. 2020] Ardakani, A. A.; Kanafi, A. R.; Acharya, U. R.; Khadem, N.; and Mohammadi, A.  2020\.  Application of deep learning technique to manage COVID-19 in routine clinical practice using CT images: Results of 10 convolutional neural networks.  Computers in Biology and Medicine 121(March):103795.
  * [Bai et al. 2020] Bai, H. X.; Hsieh, B.; Xiong, Z.; Halsey, K.; Choi, J. W.; Tran, T. M. L.; Pan, I.; Shi, L.-B.; Wang, D.-C.; Mei, J.; et al.  2020\.  Performance of radiologists in differentiating covid-19 from non-covid-19 viral pneumonia at chest ct.  Radiology 296(2):E46–E54.
  * [Cohen, Morrison, and Dao 2020] Cohen, J. P.; Morrison, P.; and Dao, L.  2020\.  Covid-19 image data collection.  arXiv preprint 2003.11597.
  * [Deng et al. 2009] Deng, J.; Dong, W.; Socher, R.; Li, L.; Li, K.; and Li, F.  2009\.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR 2009), 20-25 June 2009, Miami, Florida, USA, 248–255.  IEEE Computer Society.
  * [Diba et al. 2017] Diba, A.; Fayyaz, M.; Sharma, V.; Karami, A. H.; Arzani, M. M.; Yousefzadeh, R.; and Van Gool, L.  2017\.  Temporal 3d convnets: New architecture and transfer learning for video classification.  arXiv preprint arXiv:1711.08200.
  * [Dong and Yang 2019] Dong, X., and Yang, Y.  2019\.  Searching for a robust neural architecture in four GPU hours.  In IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2019, Long Beach, CA, USA, June 16-20, 2019, 1761–1770.  Computer Vision Foundation / IEEE.
  * [Elsken, Metzen, and Hutter 2018] Elsken, T.; Metzen, J. H.; and Hutter, F.  2018\.  Neural architecture search: A survey.  arXiv preprint arXiv:1808.05377.
  * [Faes et al. 2019] Faes, L.; Wagner, S.; Fu, D.; Liu, X.; Korot, E.; Ledsam, J.; Back, T.; Chopra, R.; Pontikos, N.; Kern, C.; Moraes, G.; Schmid, M.; Sim, D.; Balaskas, K.; Bachmann, L.; Denniston, A.; and Keane, P.  2019\.  Automated deep learning design for medical image classification by health-care professionals with no coding experience: a feasibility study.  The Lancet Digital Health 1:e232–e242.
  * [Ghoshal and Tucker 2020] Ghoshal, B., and Tucker, A.  2020\.  Estimating Uncertainty and Interpretability in Deep Learning for Coronavirus (COVID-19) Detection.  arXiv preprint arXiv:2003.10769 1–14.
  * [He et al. 2020] He, X.; Yang, X.; Zhang, S.; Zhao, J.; Zhang, Y.; Xing, E.; and Xie, P.  2020\.  Sample-efficient deep learning for covid-19 diagnosis based on ct scans.  medRxiv.
  * [He, Zhao, and Chu 2021] He, X.; Zhao, K.; and Chu, X.  2021\.  Automl: A survey of the state-of-the-art.  Knowledge-Based Systems 212:106622.
  * [Howard et al. 2017] Howard, A. G.; Zhu, M.; Chen, B.; Kalenichenko, D.; Wang, W.; Weyand, T.; Andreetto, M.; and Mobilenets, H. A.  2017\.  Efficient convolutional neural networks for mobile vision applications.  arXiv preprint arXiv:1704.04861.
  * [Jang, Gu, and Poole 2017] Jang, E.; Gu, S.; and Poole, B.  2017\.  Categorical reparameterization with gumbel-softmax.  In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings.  OpenReview.net.
  * [Jun et al. 2020] Jun, M.; Cheng, G.; Yixin, W.; Xingle, A.; Jiantao, G.; Ziqi, Y.; Minqing, Z.; Xin, L.; Xueyuan, D.; Shucheng, C.; Hao, W.; Sen, M.; Xiaoyu, Y.; Ziwei, N.; Chen, L.; Lu, T.; Yuntao, Z.; Qiongjie, Z.; Guoqiang, D.; and Jian, H.  2020\.  Covid-19 ct lung and infection segmentation dataset.  Zenodo.
  * [Kim et al. 2019] Kim, S.; Kim, I.; Lim, S.; Baek, W.; Kim, C.; Cho, H.; Yoon, B.; and Kim, T.  2019\.  Scalable neural architecture search for 3d medical image segmentation.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, 220–228.  Springer.
  * [Kingma and Ba 2015] Kingma, D. P., and Ba, J.  2015\.  Adam: A method for stochastic optimization.  In Bengio, Y., and LeCun, Y., eds., 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings.
  * [Li et al. 2020] Li, L.; Qin, L.; Xu, Z.; Yin, Y.; Wang, X.; Kong, B.; Bai, J.; Lu, Y.; Fang, Z.; Song, Q.; et al.  2020\.  Artificial intelligence distinguishes covid-19 from community acquired pneumonia on chest ct.  Radiology 200905\.
  * [Litjens et al. 2017] Litjens, G.; Kooi, T.; Bejnordi, B. E.; Setio, A. A. A.; Ciompi, F.; Ghafoorian, M.; van der Laak, J. A.; van Ginneken, B.; and Sánchez, C. I.  2017\.  A survey on deep learning in medical image analysis.  Medical Image Analysis 42(1995):60–88.
  * [Liu et al. 2018] Liu, C.; Zoph, B.; Neumann, M.; Shlens, J.; Hua, W.; Li, L.-J.; Fei-Fei, L.; Yuille, A.; Huang, J.; and Murphy, K.  2018\.  Progressive neural architecture search.  Proceedings of the European Conference on Computer Vision (ECCV) 19–34.
  * [Liu, Simonyan, and Yang 2019] Liu, H.; Simonyan, K.; and Yang, Y.  2019\.  DARTS: differentiable architecture search.  In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019.  OpenReview.net.
  * [Loshchilov and Hutter 2017] Loshchilov, I., and Hutter, F.  2017\.  SGDR: stochastic gradient descent with warm restarts.  In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings.  OpenReview.net.
  * [Mobiny et al. 2020] Mobiny, A.; Cicalese, P. A.; Zare, S.; Yuan, P.; Abavisani, M.; Wu, C. C.; Ahuja, J.; de Groot, P. M.; and Van Nguyen, H.  2020\.  Radiologist-Level COVID-19 Detection Using CT Scans with Detail-Oriented Capsule Networks.  arXiv preprint arXiv:2004.07407.
  * [Morozov et al. 2020] Morozov, S.; Andreychenko, A.; Pavlov, N.; Vladzymyrskyy, A.; Ledikhova, N.; Gombolevskiy, V.; Blokhin, I.; Gelezhe, P.; Gonchar, A.; Chernina, V.; and Babkin, V.  2020\.  Mosmeddata: Chest ct scans with covid-19 related findings.  medRxiv.
  * [Narin, Kaya, and Pamuk 2020] Narin, A.; Kaya, C.; and Pamuk, Z.  2020\.  Automatic Detection of Coronavirus Disease (COVID-19) Using X-ray Images and Deep Convolutional Neural Networks.  arXiv preprint arXiv:2003.10849.
  * [Pham et al. 2018] Pham, H.; Guan, M. Y.; Zoph, B.; Le, Q. V.; and Dean, J.  2018\.  Efficient neural architecture search via parameter sharing.  In Dy, J. G., and Krause, A., eds., Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018, volume 80 of Proceedings of Machine Learning Research, 4092–4101.  PMLR.
  * [Rahimzadeh, Attar, and Sakhaei 2020] Rahimzadeh, M.; Attar, A.; and Sakhaei, S. M.  2020\.  A fully automated deep learning-based network for detecting covid-19 from a new and large lung ct scan dataset.  medRxiv.
  * [Real et al. 2019] Real, E.; Aggarwal, A.; Huang, Y.; and Le, Q. V.  2019\.  Regularized evolution for image classifier architecture search.  In The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019, 4780–4789.  AAAI Press.
  * [Sandler et al. 2018] Sandler, M.; Howard, A. G.; Zhu, M.; Zhmoginov, A.; and Chen, L.  2018\.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In 2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018, 4510–4520.  IEEE Computer Society.
  * [Singh et al. 2020] Singh, D.; Kumar, V.; Vaishali; and Kaur, M.  2020\.  Classification of COVID-19 patients from chest CT images using multi-objective differential evolution-based convolutional neural networks.  European journal of clinical microbiology & infectious diseases : official publication of European Society of Clinical Microbiology.
  * [Tan et al. 2019] Tan, M.; Chen, B.; Pang, R.; Vasudevan, V.; Sandler, M.; Howard, A.; and Le, Q. V.  2019\.  Mnasnet: Platform-aware neural architecture search for mobile.  In IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2019, Long Beach, CA, USA, June 16-20, 2019, 2820–2828.  Computer Vision Foundation / IEEE.
  * [Tran et al. 2018] Tran, D.; Wang, H.; Torresani, L.; Ray, J.; LeCun, Y.; and Paluri, M.  2018\.  A closer look at spatiotemporal convolutions for action recognition.  In 2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018, 6450–6459.  IEEE Computer Society.
  * [Wu et al. 2019] Wu, B.; Dai, X.; Zhang, P.; Wang, Y.; Sun, F.; Wu, Y.; Tian, Y.; Vajda, P.; Jia, Y.; and Keutzer, K.  2019\.  Fbnet: Hardware-aware efficient convnet design via differentiable neural architecture search.  In IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2019, Long Beach, CA, USA, June 16-20, 2019, 10734–10742.  Computer Vision Foundation / IEEE.
  * [Zhang et al. 2020a] Zhang, J.; Xie, Y.; Li, Y.; Shen, C.; and Xia, Y.  2020a.  COVID-19 Screening on Chest X-ray Images Using Deep Learning based Anomaly Detection.  arXiv preprint arXiv:2003.12338.
  * [Zhang et al. 2020b] Zhang, K.; Liu, X.; Shen, J.; Li, Z.; Sang, Y.; Wu, X.; Zha, Y.; Liang, W.; Wang, C.; Wang, K.; et al.  2020b.  Clinically applicable AI system for accurate diagnosis, quantitative measurements, and prognosis of covid-19 pneumonia using computed tomography.  Cell.
  * [Zhao et al. 2020] Zhao, J.; He, X.; Yang, X.; Zhang, Y.; Zhang, S.; and Xie, P.  2020\.  Covid-ct-dataset: A ct scan dataset about covid-19.
  * [Zheng et al. 2020] Zheng, C.; bo Deng, X.; Fu, Q.; Zhou, Q.; Feng, J.; Ma, H.; Liu, W.; and Wang, X.  2020\.  Deep learning-based detection for covid-19 from chest ct using weak label.  medRxiv.
  * [Zhou et al. 2016] Zhou, B.; Khosla, A.; Lapedriza, À.; Oliva, A.; and Torralba, A.  2016\.  Learning deep features for discriminative localization.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016, 2921–2929.  IEEE Computer Society.
  * [Zoph and Le 2017] Zoph, B., and Le, Q. V.  2017\.  Neural architecture search with reinforcement learning.  In 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings.  OpenReview.net.
  * [Zoph et al. 2018] Zoph, B.; Vasudevan, V.; Shlens, J.; and Le, Q. V.  2018\.  Learning transferable architectures for scalable image recognition.  In 2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018, 8697–8710.  IEEE Computer Society.