\section{Introduction}\label{Sec:intro}

Recently, interest in the study of multirotor air vehicles has been
growing quickly due to their high maneuverability and many
applications, such as inspection and surveillance. In particular, the
quadrotor has been a popular platform with which to conduct UAV
(unmanned aerial vehicle) research. While the quadrotor can
tolerate partial rotor faults (see, e.g., \cite{SMGZ,ZZRT,RK,ZC}), it
becomes uncontrollable if one of its rotors fails completely (though
\cite{FLL} shows that only yaw control need be lost). Thus the
octorotor has been studied as an alternative that is more robust to
rotor failures \cite{MWY,AE,AE2,SLFFSS,SLFSF}.

An approach to addressing flight safety for UAVs in a number of recent
works has been through the use of control barrier functions (CBFs). A
CBF is a scalar-valued function on the state space of a control system
whose support (here understood as the region where the function
attains positive values) is forward invariant under appropriate controls. Thus,
one can (mathematically) guarantee the safety of the system by finding
a CBF with support contained in a safe region of the state space. In
particular, in \cite{WAE}, CBFs are used to ensure that teams of
quadrotors are able to fly in a collision free manner; in \cite{NA},
CBFs are used for obstacle avoidance in quadrotor path planning in a
surveillance scenario; and in \cite{XS}, CBFs are used for flight
safety for quadrotors under some degree of human control. Other works
using barrier functions for UAV safety include
\cite{WS,WTE,KZC}. Other applications of barrier functions include
adaptive cruise control and lane keeping \cite{AXGT,XWPGETGA,XGTA},
bipedal walking robots \cite{NS2,NHGAS}, and collision avoidance for
multirobot systems on land \cite{WAE2}. See \cite{ACENST} for a survey
of control barrier functions.

Instead of utilizing barrier functions within control algorithms to ensure
safety, as in the aforementioned references, the focus of this paper
is to use barrier functions to analyze the safety properties of
conventionally defined controllers. Specifically, our aim is to
provide a method of verifying the safety of controllers using
satisfiability modulo theories (SMT) solving. Our work is in the
spirit of \cite{JV2018,JV2019}, in which controllers for flight
vehicles are designed, and the flight envelope and asymptotic
stability of the controller are verified by the automated theorem
prover MetiTarski. However, in these works, the design of the
controllers is fairly specialized, and is closely tied together with
how the controller's safety is verified.  As above, a principal goal
of our work is to decouple safety analysis from controller design so
that our techniques can be employed even in cases where the underlying
controller can no longer be modified or where it is for other reasons
not feasible to incorporate barrier functions at the design stage.

The goal of this paper is therefore to introduce another approach to the formal
verification of the flight envelope of a UAV, one particularly that
has the potential to be used in a way that is independent of the
design of the controller. Our core method involves the use of barrier
functions, more specifically exponential barrier functions \cite{NS}
of a certain form, to produce an invariant subset of a prescribed safe
region near a desired commanded point in the control system's state
space (e.g., as one might see in connection with gain scheduling).
Our method aims to analyze systems with a controller designed to track
commands and whose dynamics can be approximated reasonably well by the
linearization around a desired operating state. In particular, our
approach lends itself well to analyzing UAV safety in states where the
vehicle is nearly upright or mildly tilting. It also has potential for
use in the verification of gain scheduled controllers that are
currently in common use for many kinds of aircraft.

The specific type of UAV we analyze is an octorotor, which, as
mentioned above, is robust to rotor failures. This allows us to
address the issue of fault tolerance in our approach. We account for
complete rotor failures, in which one or more rotors
stop working entirely and exert zero thrust, and more generally
failures in which rotors become stuck at specific speeds. To mitigate
these failures, we utilize a simple control allocator to attempt to
maintain the ideal octorotor behavior. The main property of the
control allocator that needs to be verified is that the rotors are
always capable of exerting the rotor thrusts commanded by the
allocator. We assume that failures are perfectly known, and do not
address the issue of fault detection, as in \cite{SLFFSS,SLFSF}.

In Section \ref{Sec:model}, we present the equations of motion of the
octorotor model we study. In Section \ref{Sec:ctrl}, we describe the
controller. This includes an explanation of the kinds of rotor
failures considered and how they are accounted for by the control
allocator. Section \ref{Sec:barrier} introduces the general barrier
function framework employed, and Section \ref{Sec:octo_barrier}
describes precisely the form of barrier functions used for the
octorotor. In Section \ref{Sec:ver}, we formulate in detail the
properties that are analyzed in dReal and that confirm the desired
safety properties of the octorotor controller. In Section
\ref{Sec:results}, we describe the results obtained for a specific
controller: we formally verify that the octorotor will follow commands
closely in spite of dynamic disturbances when functioning normally and
under various rotor failures, and we identify some combinations of
rotor failures in which the desired conditions are violated. Finally,
concluding remarks are found in Section \ref{Sec:con}.

\section{Octorotor dynamics model}\label{Sec:model}

In this section, we describe the octorotor dynamics model that we
analyze. We follow \cite{GLL,LLM} in modeling the rigid body dynamics
and some characteristics of the rotors, and our octorotor model is
based on that in \cite{MWY}. See also \cite{MKC} for an introduction
to the modeling of multirotor aerial vehicles. The octorotor consists
of eight identical rotors arranged in the shape of a regular octagon as
illustrated in Figure \ref{Fig:octo_diag}. Each rotor exerts an upward
thrust normal to the plane of the octorotor, and a torque about the
octorotor's center of mass as explained below.

To express the octorotor's dynamics, we first describe the two
reference frames required. First, the world reference
frame, which we assume is inertial, is given by the standard unit
basis vectors $e_1=[1,0,0]^T$, $e_2=[0,1,0]^T$, and
$e_3=[0,0,1]^T$. We assume that $e_1$ points north, $e_2$ points east,
and $e_3$ points down. Next, we assume that the body of the octorotor
has reference frame given by the orthonormal (vertical) vectors
$b_1,b_2,b_3$, where $b_1$ is considered the forward direction with
respect to the octorotor, $b_2$ is the right direction, and $b_3$ is
the down direction. (In Figure \ref{Fig:octo_diag}, $b_3$ points into
the page.) The $3\times3$ matrix $[b_1,b_2,b_3]$ is denoted by $R$,
which is the rotation matrix from the body frame to the world frame.

Now we describe the state space of the system. We denote the
inertial position of the octorotor's center of mass by $r=[x,y,z]^T$;
and the inertial velocity by $v=[v_x,v_y,v_z]^T$.  Another component
of the state of the octorotor is its orientation, which is given by
$R$. Alternatively, we express the orientation in terms of Euler
angles. Letting $\phi,\theta,\psi$ be the roll, pitch, and yaw angles,
respectively, corresponding to the $Z$-$Y$-$X$ sequence for Euler
angles, we have the following relationship between $R$ and
$[\phi,\theta,\psi]$:
\begin{equation*}
\begin{aligned}
R&=
\begin{bmatrix}
c_\psi&-s_\psi&0\\
s_\psi&c_\psi&0\\
0&0&1
\end{bmatrix}
\begin{bmatrix}
c_\theta&0&s_\theta\\
0&1&0\\
-s_\theta&0&c_\theta
\end{bmatrix}
\begin{bmatrix}
1&0&0\\
0&c_\phi&-s_\phi\\
0&s_\phi&c_\phi
\end{bmatrix}\\
&=\begin{bmatrix}
c_\theta c_\psi&s_\phi s_\theta c_\psi-c_\phi s_\psi&c_\phi s_\theta c_\psi+s_\phi s_\psi\\
c_\theta s_\psi&s_\phi s_\theta s_\psi+c_\phi s_\psi&c_\phi s_\theta s_\psi-s_\phi c_\psi\\
-s_\theta&s_\phi c_\theta&c_\phi c_\theta
\end{bmatrix}.
\end{aligned}
\end{equation*}
Here and below, we let $s_*=\sin(*)$, $c_*=\cos(*)$, and
$t_*=\tan(*)$. The final component of the state space is the angular
velocity vector $\Omega=[\Omega_1,\Omega_2,\Omega_3]^T$ of the
octorotor with respect to the body reference frame.

Next, we describe the control space of the octorotor. The control
space consists of the eight rotors, each of which exerts a force and
torque on the octorotor. For simplicity, we assume that the thrust of
each rotor can be directly set, and we do not model the motor
dynamics. Also, we assume that each rotor can rotate only in one
direction as shown in Figure \ref{Fig:octo_diag}, and that the thrust
exerted is in the upward direction $-b_3$ with respect to the
octorotor. Let $f_j$, $j=1,\ldots,8$, denote the forces exerted by the
rotors. We assume that each $f_j$ is restricted to an interval
$[f_{min},f_{max}]$, where $f_{min}$ and $f_{max}$ are the least and
greatest possible forces that each rotor can exert on the octorotor,
respectively. Since we are assuming that the rotors can only rotate in
one direction, we have $f_{min}\geq0$. Let $d$ denote the distance
from the center of mass of the octorotor to the center of each rotor, and
note that the thrust of the $j$th rotor induces a torque of magnitude
$d\cdot f_j$ about an axis in the octorotor plane perpendicular to the
arm containing the rotor. The $j$th rotor also induces a torque about
the $b_3$ axis of magnitude $c\cdot f_j$, where $c$ is the ratio of
induced torque to thrust. In particular, the rotors that rotate
counterclockwise induce a torque in the direction of $b_3$, and those
that rotate clockwise induce a torque in the direction of $-b_3$.

\begin{figure}
\centering
\includegraphics[width=0.7\columnwidth]{octorotor_diagram.png}
\caption{Octorotor diagram}
\label{Fig:octo_diag}
\end{figure}

Let $F$ denote the total thrust exerted by the rotors, and let
$\tau_1,\tau_2,\tau_3$ denote the net torque exerted by the rotors
about the body axes $b_1,b_2,b_3$, respectively. Then letting
$\gamma=\frac{\pi}{8}$, the relationship between
$u=[F,\tau_1,\tau_2,\tau_3]^T$ and $f_{all}=[f_1,f_2,\ldots,f_8]^T$ is
given by $u=\Lambda f_{all}$, where $\Lambda$ is the matrix
\[{\footnotesize
\begin{bmatrix}
1&1&1&1&1&1&1&1\\
ds_\gamma&-ds_\gamma&-dc_\gamma&-dc_\gamma&-ds_\gamma&ds_\gamma&dc_\gamma&dc_\gamma\\
dc_\gamma&dc_\gamma&ds_\gamma&-ds_\gamma&-dc_\gamma&-dc_\gamma&-ds_\gamma&ds_\gamma\\
c&c&-c&-c&c&c&-c&-c
\end{bmatrix}}.\]

\iffalse
\begin{equation}\label{rotortopseudo}
\begin{bmatrix}F\\\tau_1\\\tau_2\\\tau_3\end{bmatrix}=
\begin{bmatrix}
1&1&1&1&1&1&1&1\\
ds_\gamma&-ds_\gamma&-dc_\gamma&-dc_\gamma&-ds_\gamma&ds_\gamma&dc_\gamma&dc_\gamma\\
dc_\gamma&dc_\gamma&ds_\gamma&-ds_\gamma&-dc_\gamma&-dc_\gamma&-ds_\gamma&ds_\gamma\\
c&c&-c&-c&c&c&-c&-c
\end{bmatrix}
\begin{bmatrix}f_1\\f_2\\\vdots\\f_8\end{bmatrix}.
\end{equation}
Let $\Lambda$ denote the $4\times8$ matrix above so that $u=\Lambda f_{all}$.
\fi

Next, let $m\in\mathbb R$ denote the mass of the octorotor,
$g=9.81\,\frac{\mathrm{m}}{\mathrm{s}^2}$ be the gravitational
  acceleration, and
\[J=\begin{bmatrix}
J_1&0&0\\
0&J_2&0\\
0&0&J_3
\end{bmatrix}\]
be the inertia matrix of the octorotor with respect to its body frame. Then let $\Delta_r=[\Delta_x,\Delta_y,\Delta_z]^T$ and $\Delta_R=[\Delta_{R,1},\Delta_{R,2},\Delta_{R,3}]^T$ denote unstructured force and torque disturbances due to dynamics that are unaccounted for. In particular, we view these disturbances as encompassing aerodynamic effects such as air drag and blade flapping, which, for simplicity, we do not model in detail. (See, e.g., \cite{HHWT, MKC} for more information.) We assume the disturbances are unknown, untracked, and measurable functions of time, and satisfy the inequalities $|\Delta_x|,|\Delta_y|,|\Delta_z|\leq\Delta_{r,max}$, and $|\Delta_{R,1}|,|\Delta_{R,2}|\leq\Delta_{R,12,max}$, $|\Delta_{R,3}|\leq\Delta_{R,3,max}$, where $\Delta_{r,max}$, $\Delta_{R,12,max}$, and $\Delta_{R,3,max}$ are fixed bounds. Then the differential equations describing the motion of the octorotor are
\begin{equation*}
\begin{aligned}
\dot{r}&=v\\
m\dot{v}&=\begin{bmatrix}
0\\ 0\\ mg
\end{bmatrix}
-F
\begin{bmatrix}
c_\phi s_\theta c_\psi+s_\phi s_\psi\\
c_\phi s_\theta s_\psi-s_\phi c_\psi\\
c_\phi c_\theta
\end{bmatrix}
+\Delta_r\\
\begin{bmatrix}
\dot{\phi}\\ \dot{\theta}\\ \dot{\psi}
\end{bmatrix}
&=\begin{bmatrix}
1& s_\phi t_\theta & c_\phi t_\theta \\
0 & c_\phi & -s_\phi \\
0 & \frac{s_\phi}{c_\theta} & \frac{c_\phi}{c_\theta}
\end{bmatrix}
\begin{bmatrix}
\Omega_1\\ \Omega_2\\ \Omega_3
\end{bmatrix}\\
J\dot{\Omega}&=\tau-\Omega\times J\Omega+\Delta_R.
\end{aligned}
\end{equation*}

In this paper, we focus on the dynamics of the components of the state vector
\[s=[v_z,\phi,\theta,\psi,\Omega_1,\Omega_2,\Omega_3]^T,\]
which we consider to be the ``inner loop'' components of the octorotor system. Note that the derivative of these components do not depend on $[x,y,z,v_x,v_y]$, and so we may safely view the components of $s$ as forming the state space of a well-defined control system.

\section{Controller}\label{Sec:ctrl}

As mentioned above, we focus on controlling the dynamics of the octorotor's vertical velocity and orientation. In particular, the controller we use is designed to have the components $(v_z,\phi,\theta,\psi)$ of the octorotor state track the command $(v_{z,d},\phi_d,\theta_d,\psi_d)$. In other words, we wish for $s$ to track $s_d=[v_{z,d},\phi_d,\theta_d,\psi_d,0,0,0]^T$. The block diagram of the full control system is shown in Figure \ref{Fig:block_diagram}.

\begin{figure}
\centering
\includegraphics[width=0.98\columnwidth]{control_system_diagram.png}
\caption{Octorotor system block diagram}
\label{Fig:block_diagram}
\end{figure}

The controller consists of two components. The first, block $K$ in the
diagram, takes as input the difference $s_d-s$ between the current
state and the commanded state, and returns the commanded vector
$u_d=[F_d,\tau_{1,d},\tau_{2,d},\tau_{3,d}]^T$ of net force and
torques. This component disregards whether the octorotor is actually
capable of executing $u_d$, that is, whether there exists a valid
$f_{all}$ which yields the desired command. So $u_d$ is not
necessarily the actual control input
$u=[F,\tau_1,\tau_2,\tau_3]^T$. To find $u$, we include a second
component, the control allocator, which takes $u_d$ as input and
returns a value for $f_{all}$ such that $f_{min}\leq f_j\leq f_{max}$
for all $j$ and the resulting control input $u=\Lambda f_{all}$ is
intended to equal, or be close to, $u_d$. Note additionally that we
use the control allocator to implement rotor failures. Specifically,
we assume that the rotor failures are always known and the control
allocator is defined so that the thrusts $f_j$ corresponding to the
failed rotors match the values at which they are stuck.

While under some circumstances the commanded value $u_d$ is not equal
to $u$, we want to ensure this never happens for the purposes of
verifying safety. In particular, we wish to find an invariant set for
which all the commanded values $u_d$ can be executed by the
octorotor. This allows us to split up our verification procedure into
$2$ separate steps as follows. The first step is to verify that a
chosen candidate invariant set is rendered invariant by the controller
under the assumption that $u$ is always equal to $u_d$. The second
step is to verify that the equality $u=u_d$ is in fact satisfied for
the control allocator for all states in the candidate invariant
set. This step is done for each rotor failure combination analyzed.

\subsection{Controller block $K$}\label{Ssec:K}

The block $K$ controller component has a simple PD structure very
similar to \cite{MWY}, and is defined by the following:
\begin{equation}\label{ctrl}
\begin{aligned}
F_d&=\frac{mg}{c_\phi c_\theta}+K_{dz}(v_z-v_{z,d})\\
\tau_{1,d}&=-K_{p\phi}(\phi-\phi_d)-K_{d\phi}\Omega_1\\
\tau_{2,d}&=-K_{p\theta}(\theta-\theta_d)-K_{d\theta}\Omega_2\\
\tau_{3,d}&=-K_{p\psi}(\psi-\psi_d)-K_{d\psi}\Omega_3.
\end{aligned}
\end{equation}
Note that $mg/c_\phi c_\theta$ is the total rotor thrust needed to
keep the vertical acceleration of the octorotor at $0$. We take the
various coefficients $K_\cdot$ from a linear quadratic regulator that
we compute for the linearized octorotor dynamics at $s=0$, and under
the assumption that $v_{z,d},\phi_d,\theta_d,\psi_d,\Delta_{r,max}$,
$\Delta_{R,12,max},\Delta_{R,3,max}$ are all equal to $0$.

\subsection{Control allocator}\label{Ssec:allocation}

To find the individual rotor thrusts exerted by the octorotor in
response to the commanded control input $u_d$, we use a simple
pseudo-inverse control allocation method (see \cite{ODB}). At this
stage, we account for the possibility that some of the rotors have
failed. So let $W\subseteq\{1,\ldots,8\}$ be the (possibly empty) set
of indices corresponding to rotors that have failed, and for $j\in W$,
let $\bar{f}_j\in[f_{min},f_{max}]$ be thrust at which rotor $j$ is
stuck. Also let $\bar{f}_j=0$ for $j\notin W$ and
$\bar{f}_{all}=[\bar{f}_1,\ldots,\bar{f}_8]^T$. Next, define $\Lambda_W$
to be the matrix with the same entries as $\Lambda$, except that
column $j$ of $\Lambda_W$ is zeroed out for all $j\in W$. Then let
$\Lambda_W^\dagger$ be the pseudo-inverse of $\Lambda_W$. Now we
define $\tilde{f}_{all}=[\tilde{f}_1,\ldots,\tilde{f}_8]^T$ by
\begin{equation}\label{tildefall}
\tilde{f}_{all}=\bar{f}_{all}+\Lambda_W^\dagger(u_d-\Lambda\bar{f}_{all}).
\end{equation}
The vector $\tilde{f}_{all}$ consists of potential rotor values, with
$\tilde{f}_j=\bar{f}_j$ for all $j\in W$ as required. However, it is
possible that $\tilde{f}_j$ falls outside $[f_{min},f_{max}]$ for some
$j\notin W$. We thus define $f_j$ for $j=1,\ldots,8$ by
\[f_j=\begin{cases}
\bar{f}_j & \text{if $j\in W$}\\
\midd\{f_{min},\tilde{f}_j,f_{max}\} & \text{if $j\notin W$,}
\end{cases}\]
where $\midd\{f_{min},\tilde{f}_j,f_{max}\}$ denotes the middle value
of the set $\{f_{min},\tilde{f}_j,f_{max}\}$. Thus
$f_j\in[f_{min},f_{max}]$ for all $j\notin W$, and so we let
$f_{all}=[f_1,\ldots,f_8]^T$. This completes the definition of the
control allocator.

To conclude this section, we wish to describe conditions for which
$u=u_d$. We first note that the matrix $\Lambda_W$ is of rank $4$ if
$|W|<4$, or $|W|=4$ and $W$ is not equal to $\{1,2,3,8\}$,
$\{1,2,4,7\}$, $\{1,2,5,6\}$, or any index set that can be obtained
from one of these sets by increasing each element by a common even
number and taking the remainder upon division by $8$ when
necessary. (This is equivalent to rotating the rotor failure patterns
corresponding to the $3$ sets about the octorotor center by multiples
of $90$ degrees.) When $\Lambda_W$ is of rank $4$, we have
$\Lambda_W\Lambda_W^\dagger=I_{4\times4}$ and
$\Lambda\tilde{f}_{all}=u_d$. So to determine whether $u=u_d$, it
suffices to check whether $\tilde{f}_j\in[f_{min},f_{max}]$ for all
$j\notin W$, and this is the condition we use, since we only consider
scenarios with up to $2$ rotor failures.

\section{Barrier Functions}\label{Sec:barrier}

The safety conditions that we wish to prove rely on the notion of
invariant sets. An invariant set for a dynamical system is a region
$I$ of the state space of the system such that if the state begins in
$I$, then the state remains in $I$ at all future times. Our goal is to
ensure that for all commands $s_d$ one wishes to give the octorotor,
there exists an invariant set of octorotor states $s$ which are close
to $s_d$. Obtaining and verifying such invariant sets allows one to
gauge with certainty how close the controller is able to keep the
octorotor state to the given commands in spite of the disturbances
given in the dynamics model.

To find these invariant sets, we utilize barrier functions. The type
of barrier function we use is a form of exponential barrier function
as introduced by \cite{NS}. (See, in particular, Remark 5, Proposition
1, and the related discussion.) With this concept, we find a set of
linear inequalities on the state space which forms a candidate
invariant region which we can verify. In this section, we provide a
brief overview of the reasoning behind exponential barrier functions,
and then explain the general form of exponential barrier functions we
use.

Let $\dot{x}=f(x)$ define a dynamical system on some set
$X\subseteq\mathbb R^n$, and let $h=h_0:X\rightarrow\mathbb R$ be a differentiable
function. The basic logic underlying exponential barrier functions is
that, for a trajectory $x(t)$ of the dynamical system, $h(x(t))\geq0$
for all $t\geq0$ whenever $h(x(0))\geq0$ and
\begin{equation*}
  \begin{aligned}
    &\frac{d}{dt}(h(x(t)))+p_1h(x(t))\\
    &=(\nabla h)(x(t))\cdot f(x(t))+p_1h(x(t))\geq0
  \end{aligned}
\end{equation*}
for all $t\geq0$, where $p_1>0$ is a constant. Thus, if one seeks an
invariant subset of $\{x\in X:h(x)\geq0\}$, one can, if necessary,
enforce the extra inequality $h_1(x):=(\nabla h)(x)\cdot
f(x)+p_1h(x)\geq0$, and form the set $\{x\in
X:h_0(x),h_1(x)\geq0\}$. If one can then prove that whenever the
initial condition is chosen from this set, the inequality
$h_1(x)\geq0$ holds indefinitely, then $h_0(x)\geq0$ holds
automatically. One can furthermore inductively form a sequence of
functions $h_j$ defined by $h_j(x)=(\nabla h_{j-1})(x)\cdot
f(x)+p_jh_{j-1}(x)$ and examine the sets $\{x\in
X:h_j(x)\geq0,j=0,\ldots,n\}$ to see if any of them are invariant, and
in particular, checking if the inequality $h_n(x)\geq0$ holds
indefinitely when $h_0(x),\ldots,h_n(x)\geq0$ holds at the initial
condition. See \cite{NS} for more details regarding the reasoning
above in the context of controlled dynamical systems.

In light of the above discussion, we now explain the general barrier
function framework we use in constructing invariant sets for the
octorotor. It relies on a few starting assumptions. First, we assume
that a controller is already in place (namely, one of the form
described in Section \ref{Sec:ctrl}), which seeks to track a fixed
commanded state, and the dynamics are governed by a differential
equation of the form $\dot{x}=f(x,x_d,d)$, where $x\in\mathbb R^n$ is the
state, $x_d\in\mathbb R^d$ is the commanded state, $d\in\mathbb R^m$ denotes any
potential disturbances (all three being vertical vectors), and $f$ is
Lipschitz on any candidate invariant set of interest. Next, we assume
that the dynamics can be approximated sufficiently well by a
linearization of the form $\dot{x}=A(x-x_d)$, $A\in\mathbb R^{n\times n}$. In
the case of the octorotor, we assume $A$ results from the
linearization of the dynamics where $s=0$. Next, we assume that the
conditions we desire to be enforced are of the form
$P_i(x-x_d)+D_i\geq0$, where $P_i$ is a horizontal vector of length
$n$, and $D_i$ is a positive scalar. The aim of these inequalities is
to form a starting set of states $x$ nearby $x_d$ in which to find an
invariant set. The goal is then to confirm the ability of the
controller to keep the state $x$ near $x_d$.

Now we form a sequence of barrier functions that aim to enforce
$P_i(x-x_d)+D_i\geq0$. First, we let $p_{i,j},\delta_{i,j}\geq0$ be
constants for $j=1,\ldots,n_i$ such that
$\sum_{j=1}^{n_i}\delta_{i,j}<D_i$. Then we define the sequence of
functions
\begin{equation*}
\begin{aligned}
\tilde{h}_{i,0}(x)=\,&P_i(x-x_d)+D_i\\
\tilde{h}_{i,1}(x)=\,&P_i(I+p_{i,1}A)(x-x_d)+D_i-\delta_{i,1}\\
\tilde{h}_{i,2}(x)=\,&P_i(I+p_{i,1}A)(I+p_{i,2}A)(x-x_d)\\
\,&+D_i-\delta_{i,1}-\delta_{i,2}\\
\vdots\,&\\
\tilde{h}_{i,n_i}(x)=\,&P_i(I+p_{i,1}A)\cdots(I+p_{i,n_i}A)(x-x_d)\\
\,&+D_i-\delta_{i,1}-\cdots-\delta_{i,n_i}.
\end{aligned}
\end{equation*}
where $I$ is the $n\times n$ identity matrix. Notice that under the
linear dynamics,
$\tilde{h}_{i,j}(x)\leq\tilde{h}_{i,j-1}(x)+p_{i,j}\frac{d}{dt}(\tilde{h}_{i,j-1}(x))$
for $j=1,\ldots,n_i$, and so if the inequality $\tilde{h}_{i,j}\geq0$
holds, then $\tilde{h}_{i,j-1}$ holds automatically. Next, we let $\mu>0$
be a parameter and form our barrier function sequence
\begin{equation*}
\begin{aligned}
h_{i,0,\mu}(x)&=\frac{P_i(x-x_d)}{D_i}+\mu\\
h_{i,1,\mu}(x)&=\frac{P_i(I+p_{i,1}A)(x-x_d)}{D_i-\delta_{i,1}}+\mu\\
h_{i,2,\mu}(x)&=\frac{P_i(I+p_{i,1}A)(I+p_{i,2}A)(x-x_d)}{D_i-\delta_{i,1}-\delta_{i,2}}+\mu\\
&\vdots\\
h_{i,n_i,\mu}(x)&=\frac{P_i(I+p_{i,1}A)\cdots(I+p_{i,n_i}A)(x-x_d)}{D_i-\delta_{i,1}-\cdots-\delta_{i,n_i}}+\mu,
\end{aligned}
\end{equation*}
and notice that $h_{i,j,1}$ is a postive scalar multiple of
$\tilde{h}_{i,j}$. We refer to the number $n_i$ as the depth of the
above sequence. Then $I(\mu)=\{x:h_{i,j,\mu}(x)\geq0,\forall i,\forall
j\}$ forms a parameterized family of candidate invariant sets, and
$I(1)$ in particular is the candidate set that enforces the conditions
$P_i(x-x_d)+D_i\geq0$. The purpose of the parameter $\mu$ is to allow
one to potentially establish a region of attraction for $I(1)$ by
showing that when the initial system state is in $I(\mu_{max})$ for
some $\mu_{max}>1$, the state will enter $I(\mu)$ for smaller and
smaller $\mu$ over time until it reaches $I(1)$.

In order the confirm whether $I(\mu)$ is actually invariant for
$\mu\in[1,\mu_{max}]$, it is sufficient to check that for all $x\in
I(\mu_{max})$, we have that, for all $i$ and $j$, $h_{i,j,\mu}(x)=0$
implies that $\frac{d}{dt}(h_{i,j,\mu}(x))>0$. We shall refer to this
implication as the invariance condition for $h_{i,j,\mu}$. Under this
condition, it is clear that the least $\mu$ satisfying
$h_{i,j,\mu}(x)\geq0$ decreases with respect to time while
$\mu\geq1$.

One may notice that under the linearized dynamics, it is only
necessary to check the invariance condition for $h_{i,n_i,\mu}$ for
all $i$, since in this case, the fact that $h_{i,n_i,\mu}(x)\geq0$
holds implies that $h_{i,j,\mu}(x)\geq0$ holds for all $j\leq
n_i$. This continues to be true if we let $\delta_{i,j}=0$ for all $i$
and $j$. However, since we are interested in verifying invariance for
the original nonlinear dynamics, it is necessary to check the
invariance conditions for all functions.
Also, under the nonlinear dynamics, the fact that
$h_{i,j,\mu}(x)\geq0$ holds may not automatically imply that
$h_{i,j-1,\mu}(x)\geq0$ for $\delta_{i,j}=0$, though this issue can be
mitigated by increasing $\delta_{i,j}$.

To conclude this section, we remark that in the case where $I(\mu)$ is
proven to be invariant for $\mu\in[1,\mu_{max}]$, it is straightforward
to check for an expanded region of attraction for $I(1)$ by verifying the
invariance conditions for $\mu\in[\mu_{max},\tilde{\mu}_{max}]$, where
$\tilde{\mu}_{max}>\mu_{max}$; and similarly one can check to see if
there are smaller invariant sets inside $I(1)$ by verifying the invariance
conditions for $\mu\in[\mu_{min},1]$, for $\mu_{min}\in(0,1)$. One can
then in a systematic way, e.g., through a bisection method,
approximate the exact $\mu$-interval over which the invariant conditions
are satisfied as accurately as desired; though SMT solving time may
limit the accuracy one can feasibly obtain.

\section{Octorotor Barrier Functions}\label{Sec:octo_barrier}

In this section, we describe in more detail the form of the barrier
functions we use for the octorotor system. As mentioned above, the
approximating linear dynamics we use is the linearization around
$s=0$. Conveniently, this approximation decomposes into $4$
subsystems. The first subsystem consists of only the vertical velocity
$v_z$. The nonlinear equation governing $v_z$, without disturbances,
is
\[\dot{v}_z=g-\frac{F}{m}c_\phi c_\theta,\]
which, when combined with the controller defined by \eqref{ctrl}, yields
\[\dot{v}_z=-\frac{K_{dz}c_\phi c_\theta}{m}(v_z-v_{z,d}).\]
(Throughout this section, we assume that $u=u_d$.) Hence, the linearized system for $v_z$ is
\[\dot{v}_z=-\frac{K_{dz}}{m}(v_z-v_{z,d}).\]
One of our goals is to ensure that $|v_z-v_{z,d}|\leq D_{v_z}$ for
some constant $D_{v_z}>0$, assuming the bound $|v_{z,d}|\leq
D_{v_{z,d}}$ with $D_{v_{z,d}}>0$. So we define our barrier function
components for $v_z$ by
\[h_{v_z,\mu}^\pm(s,v_{z,d})=\pm\frac{v_z-v_{z,d}}{D_{v_z}}+\mu,\]
one barrier function corresponding to $+$ and one to $-$. This means
that we do not make use of barrier functions of higher depth to
enforce $|v_z-v_{z,d}|\leq D_{v_z}$.

The next subsystem consists of the roll angle and rate
$(\phi,\Omega_1)$. The nonlinear equations governing these variables,
without disturbances, are
\[\dot{\phi}=\Omega_1+t_\theta(s_\phi\Omega_2+c_\phi\Omega_3),\quad
\dot{\Omega}_1=\frac{\tau_1+\Omega_2\Omega_3(J_2-J_3)}{J_1}.\]
Linearizing and combining with the controller \eqref{ctrl} yields
\[\dot{\phi}=\Omega_1,\quad
\dot{\Omega}_1=\frac{-K_{p\phi}(\phi-\phi_d)-K_{d\phi}\Omega_1}{J_1}.\]
Another goal is to ensure that $|\phi-\phi_d|\leq D_\phi$ for some
$D_\phi>0$, assuming that $|\phi_d|\leq D_{\phi_d}$, $D_{\phi_d}>0$,
and we aim to enforce the first inequality with barrier function
sequences of depth $1$. So for constants
$p_{\phi,1},\delta_{\phi,1}\geq0$ with $\delta_{\phi,1}<D_\phi$, we
define our barrier function components for $(\phi,\Omega_1)$ as the
following:
\begin{equation*}
\begin{aligned}
h_{\phi,0,\mu}^\pm(s,\phi_d)&=\pm\frac{\phi-\phi_d}{D_\phi}+\mu,\\
h_{\phi,1,\mu}^\pm(s,\phi_d)&=\pm\frac{\phi-\phi_d+p_{\phi,1}\Omega_1}{D_\phi-\delta_{\phi,1}}+\mu.
\end{aligned}
\end{equation*}

There are two remaining subsystems, one of which consists of the pitch
angle and rate $(\theta,\Omega_2)$, and the other consists of the yaw
angle and rate $(\psi,\Omega_3)$. Both subsystems are very similar to
that of the roll angle and rate. The nonlinear equations governing the
dynamics of these variables are
\begin{equation*}
\begin{aligned}
\dot{\theta}&=c_\phi\Omega_2-s_\phi\Omega_3, & \dot{\psi}&=\frac{s_\phi\Omega_2+c_\phi\Omega_3}{c_\theta}\\
\dot{\Omega}_2&=\frac{\tau_2+\Omega_1\Omega_3(J_3-J_1)}{J_2}, & \dot{\Omega}_3&=\frac{\tau_3+\Omega_1\Omega_2(J_1-J_2)}{J_3}.
\end{aligned}
\end{equation*}
Linearizing and then combining these equations with the controller \eqref{ctrl} yields
\begin{equation*}
\begin{aligned}
\dot{\theta}&=\Omega_2, & \dot{\Omega}_2&=\frac{-K_{p\theta}(\theta-\theta_d)-K_{d\theta}\Omega_2}{J_2}\\
\dot{\psi}&=\Omega_3, & \dot{\Omega}_3&=\frac{-K_{p\psi}(\psi-\psi_d)-K_{d\psi}\Omega_3}{J_3}.
\end{aligned}
\end{equation*}
Our goals for these variables is to ensure that $|\theta-\theta_d|\leq
D_\theta$ and $|\psi-\psi_d|\leq D_\psi$ for some $D_\theta,D_\psi>0$,
assuming that $|\theta_d|\leq D_{\theta_d}$ and $|\psi_d|\leq
D_{\psi_d}$, $D_{\theta_d},D_{\psi_d}>0$. As with the roll angle and
rate, we utilize barrier function sequences of depth 1. So for
constants
$p_{\theta,1},\delta_{\theta,1},p_{\psi,1},\delta_{\psi,1}\geq0$ with
$\delta_{\theta,1}<D_\theta$ and $\delta_{\psi,1}<D_\psi$, we define
the barrier function components for $(\theta,\Omega_2)$ and
$(\psi,\Omega_3)$ as the following:
\begin{equation*}
\begin{aligned}
h_{\theta,0,\mu}^\pm(s,\theta_d)&=\pm\frac{\theta-\theta_d}{D_\theta}+\mu,\\
h_{\theta,1,\mu}^\pm(s,\theta_d)&=\pm\frac{\theta-\theta_d+p_{\theta,1}\Omega_2}{D_\theta-\delta_{\theta,1}}+\mu\\
h_{\psi,0,\mu}^\pm(s,\psi_d)&=\pm\frac{\psi-\psi_d}{D_\psi}+\mu,\\
h_{\psi,1,\mu}^\pm(s,\psi_d)&=\pm\frac{\psi-\psi_d+p_{\psi,1}\Omega_3}{D_\psi-\delta_{\psi,1}}+\mu.
\end{aligned}
\end{equation*}
Lastly, we also include the following functions to ensure bounds of
the form $|\Omega_1|\leq D_{\Omega_1}$, $|\Omega_2|\leq D_{\Omega_2}$,
and $|\Omega_3|\leq D_{\Omega_3}$, where
$D_{\Omega_1},D_{\Omega_2},D_{\Omega_3}>0$:
\begin{equation*}
\begin{aligned}
h_{\Omega_1,\mu}^\pm(s)&=\pm\frac{\Omega_1}{D_{\Omega_1}}+\mu,\\
h_{\Omega_2,\mu}^\pm(s)&=\pm\frac{\Omega_2}{D_{\Omega_2}}+\mu,\\
h_{\Omega_3,\mu}^\pm(s)&=\pm\frac{\Omega_3}{D_{\Omega_3}}+\mu.
\end{aligned}
\end{equation*}
This completes the description of all the barrier functions.

\section{Verification process}\label{Sec:ver}

We will now formulate the statements about the octorotor system that
we formally analyze in dReal. The dReal SMT solver \cite{GKC} is
capable of analyzing statements in the theory of first-order real
arithmetic with non-linear function symbols (including, crucially,
transcendental functions) in order to determine their satisfiability.
In order to overcome well-known decidability issues with this theory,
dReal incorporates a numerical precision constant that mediates the
granularity of the proof search. Given a statement $\varphi$ in this
theory, dReal will reply with one of the following:
\begin{itemize}
\item An assignment of ranges (intervals) to the (free) variables of
  $\varphi$  that make the statement satisfiable modulo the numerical
  precision; or
\item A proof that there is no assignment that would make $\varphi$ satisfiable.
\end{itemize}
The first case is referred to as a \emph{$\delta$-SAT} result.  The second
case is referred to as an \emph{UNSAT} result. In general, it is not
possible to determine, given the assignment corresponding to a
$\delta$-SAT result whether or not there are values in the range of
the assignment that would make the statement satisfiable.  However,
one can sample from these intervals and evaluate the results in order
to further search for a genuinely satisfiable (\emph{SAT}) result.  In
practice, when running dReal we always carry out a na\"ive search
within $\delta$-SAT ranges by polling the midpoint of the returned
box.  Thus, in the sequel \emph{SAT} means that we have found a
genuine counter-example via this further polling of the variable
ranges returned as part of one of dReal's $\delta$-SATs.  We always
formulate safety properties in such a way that the goal is to obtain
an UNSAT result, i.e., a proof of unsatisfiability, from dReal.

First, we introduce a constant
$\epsilon$, a small positive number that we use for various purposes,
one of which is to make our safety conditions slightly more strict to
account for potential numerical errors. Next, we define a search space
over which the statements we formulate are proven. In particular, we
restrict all of the following variables to a symmetric interval about
$0$:\footnote{The regions could be asymmetric and for other kinds of
  tasks that might be more natural.  For the task at hand, symmetry of
  the regions is most reasonable.}
\[v_z,\phi,\theta,\psi,\Omega_1,\Omega_2,\Omega_3,\Delta_z,\Delta_{R,1},\Delta_{R,2},\Delta_{R,3},v_{z,d},\phi_d,\theta_d,\psi_d.\]
That is, for all variables $\xi$ above, we restrict $\xi$ to lie in
the interval $(-\xi_{max},\xi_{max})$, where $\xi_{max}>0$. We also
restrict the variable $\mu$ to the interval
$(1-\epsilon,\mu_{max}+\epsilon)$. This ensures that dReal conducts a
well bounded search, which avoids unnecessary computational
issues. As mentioned above, all conditions are expressed in terms of a
search for a state which \emph{violates} the desired
safety condition. The aim is then to use dReal to verify that such a
violation cannot exist, which proves the condition.

\subsection{Barrier function support is inside search space}
Now the first thing we wish to ensure is that the common support of
the barrier functions we define is a subset of our search space. This
is because it is necessary to search over the full
boundary of the support when checking invariance. For this property,
we formulate the dReal search as
follows. We search for a state $s$ near the boundary of the search
space ---i.e., satisfying $\xi>\xi_{max}-\epsilon$ or $\xi<-\xi_{max}+\epsilon$
for some $\xi\in\{v_z,\phi,\theta,\psi,\Omega_1,\Omega_2,\Omega_3\}$---
such that there exists a barrier function component $h$ and, if
needed, associated command $c$ such that $s$ is in or near the support
of $h(\cdot,c)$, i.e., $h(s,c)>-\epsilon$. If $h$ has no argument for
a command $c$, then the last inequality should be replaced by
$h(s)>-\epsilon$. If dReal verifies that no such state $s$ exists,
then we can conclude that the full support is inside the search
space.

\subsection{Invariance of barrier function support}
Next, we formulate the conditions for checking that the barrier
function support is invariant, and we do this by checking the
invariance property for each barrier function component. Here we are
assuming that the controller defined by \eqref{ctrl} always works,
i.e., $u=u_d$ and we disregard the rotor thrust limits. Also, we
assume that the states $s$ we search over are in or near the support
of every barrier function component, i.e., for all barrier function
components $h$, we have $h(s,c)>-\epsilon$ for some associated command
$c$, if such a command is needed, or otherwise $h(s)>-\epsilon$. With
these assumptions in mind, we formulate one condition for each barrier
function component. For the component $h_*$, we search for a state $s$
and, if needed, a command $c$ associated with $h_*$ such that $s$ is
on or near the boundary of the support of $h_*(\cdot,c)$ (or
$h_*(s)$), i.e., $-\epsilon<h_*(s,c)<\epsilon$ (or
$-\epsilon<h_*(x)<\epsilon$), and the time derivative of $h_*(s,c)$
(or $h_*(s)$) is almost negative, i.e.,
$\frac{d}{dt}(h_*(s,c))<\epsilon$ (or
$\frac{d}{dt}(h_*(s))<\epsilon$). For the sake of clarity, the
condition we impose on the state $s$ to consider it a violation of the
invariance condition for $h_*$ is the following:
\begin{equation*}
\begin{aligned}
&\,\left({\text{for every b.f.\ component $h$, there exists an}\atop
\text{associated command $c$ such that $h(s,c)>-\epsilon$}}\right)\\
\wedge&\,\left({\text{there exists a command $c$ associated to $h_*$}\atop
\text{such that $h_*(s,c)<\epsilon$ and $\frac{d}{dt}(h_*(s,c))<\epsilon$}}\right).
\end{aligned}
\end{equation*}
For $h$ or $h_*$ that has no argument for a command in the above expression, we disregard the corresponding command $c$ and replace $h(s,c)$ with $h(s)$, or $h_*(s,c)$ with $h_*(s)$. If, for every barrier function component $h_*$, there is no state $s$ that satisfies the above condition, then the barrier function support is invariant under the controller defined by \eqref{ctrl}.

\subsection{Rotor thrust bounds}
The last condition we check is the equality $u=u_d$ we assumed in the previous condition. More specifically, we wish to check whether this equality is true over the support of the given barrier function and under the rotor failure combination being considered. As explained in Section \ref{Ssec:allocation}, it suffices to check whether the thrusts $\tilde{f}_j$ of the non-failed rotors given by the allocation method are valid and stay within the interval $[f_{min},f_{max}]$, that is, as long as $\rank(\Lambda_W)=4$. We formulate a condition for each combination of rotor failures, so let $W\subseteq\{1,\ldots,8\}$ such that $\rank(\Lambda_W)=4$. It is clear we can consider $u_d$ a function of $s$ by definition; and as a result, by \eqref{tildefall} we can view $\tilde{f}_{all}$ as a function of $s$. With this in mind, we can formulate our search as follows. We search for a state $s$ such that for every barrier function component $h$, $h(s,c)>-\epsilon$ for some corresponding command $c$ (or $h(s)>-\epsilon)$, and additionally, for some $j\notin W$, $\tilde{f}_j$ is outside, or nearly outside, of $[f_{min},f_{max}]$, i.e., $\tilde{f}_j-f_{max}>-\epsilon$ or $\tilde{f}_j-f_{min}<\epsilon$. If dReal verifies that no such state $s$ exists, then the octorotor controller behaves as intended for the rotor failure combination being analyzed.

\section{Results}\label{Sec:results}

We apply our verification process to an octorotor model with the characteristics of that in \cite{MWY}. Specifically, the octorotor model has mass $m=1.2\,\text{kg}$, inertia matrix entries $J_1=J_2=7.5\times10^{-3}\,\text{kg}\cdot\text{m}^2$ and $J_3=1.3\times10^{-2}\,\text{kg}\cdot\text{m}^2$, arm length $d=0.4\,\text{m}$, and torque-to-thrust ratio $c=(\frac{7.5}{3.13})\times10^{-2}\,\text{m}$. Additionally, we let the minimum and maximum rotor thrusts be $f_{min}=0$ and $f_{max}=\frac{1}{2}mg$, and let the maximum disturbance magnitudes be $\Delta_{r,max}=\frac{1}{10}mg$, $\Delta_{R,12,max}=(0.6\,\frac{\text{rad}}{\text{s}^2})J_1$, and $\Delta_{R,3,max}=(0.6\,\frac{\text{rad}}{\text{s}^2})J_3$.

Next, to find the coefficients of the controller \eqref{ctrl}, we compute the linear quadratic regulator for the linearization of the octorotor system at $s=0$, while assuming that $s_d=0$ and neglecting the disturbances $\Delta_r,\Delta_R$. Specifically, we find the linear controller that minimizes the objective function
\begin{equation*}
\begin{aligned}
&\int_0^\infty(s^TQs+|u|^2)\,dt,\\
&\text{with } Q=\diag\left(40,\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{8},\frac{1}{8},\frac{1}{8}\right),
\end{aligned}
\end{equation*}
under the linearized dynamics. The resulting coefficients in \eqref{ctrl} are then (approximately) as follows:
\begin{equation}\label{K_coeffs}
\begin{aligned}
&K_{dz}=6.32,\quad K_{p\phi}=K_{p\theta}=K_{p\psi}=0.5,\\
&K_{d\phi}=K_{d\theta}=0.364,\quad K_{d\psi}=0.371.
\end{aligned}
\end{equation}

Next, to define the barrier function components, we make the following assignments:
\begin{equation*}
\begin{gathered}
D_{v_z}=0.25,\enspace D_{v_{z,d}}=1,\enspace D_\phi=D_\theta=D_\psi=0.05,\\
D_{\phi_d}=D_{\theta_d}=0.15,\enspace D_{\psi_d}=\epsilon,\enspace D_{\Omega_1}=D_{\Omega_2}=D_{\Omega_3}=0.09,\\
p_{\phi,1}=p_{\theta,1}=p_{\psi,1}=0.7,\enspace \delta_{\phi,1}=\delta_{\theta,1}=\delta_{\psi,1}=0.017,\\
\mu_{max}=2.
\end{gathered}
\end{equation*}
This means that the command $v_{z,d}$ is allowed to span the interval $(-1,1)$, and $\phi_d$ and $\theta_d$ are allowed to span $(-0.15,0.15)$. Note that there is no loss in letting $D_{\psi_d}=\epsilon$, hence $\psi_d\approx0$, since one can always redefine coordinates so that $\psi_d=0$. Also, for any combination of commands, the candidate invariant set resulting from the barrier functions is the set $I(1)$ of states $s$ such that $|v_z-v_{z,d}|\leq0.25$, $(\phi,\Omega_1)\in S_{\phi_d}$, $(\theta,\Omega_2)\in S_{\theta_d}$, and $(\psi,\Omega_3)\in S_{\psi_d}$, where
\[S_{\alpha_d}=\left\{(\alpha,\omega)\in\mathbb R^2:{|\alpha-\alpha_d|\leq0.05,|\omega|\leq0.09\atop |\alpha-\alpha_d+0.7\omega|\leq0.033}\right\}.\]
The set $S_{\alpha_d}$ is the dark blue region illustrated in Figure \ref{Fig:S_alpha_d}. The light blue region is the dilation of $S_{\alpha_d}$ by a factor of $2$, and corresponds to the region of attraction for the invariant set. Lastly, we let $\epsilon=10^{-8}$.

\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth]{S_alpha_d_new_2.png}
\caption{The region $S_{\alpha_d}$}
\label{Fig:S_alpha_d}
\end{figure}

With all of the above definitions, we used dReal to obtain a proof of
safety in the case of no rotor failures, as well as many combinations
of rotor failures. That is, we used dReal to successfully verify that
the set $I(\mu)$ is invariant for $1\leq\mu\leq2$ under the controller
\eqref{ctrl} with coefficients \eqref{K_coeffs}, and that the
octorotor is capable of executing the controller's actions under
various rotor failure combinations. Table \ref{Tab:1} shows the times
it took to prove the various conditions that ensure safety in the case
where there are no rotor failures. Recall that UNSAT cases correspond
to safety proofs found and SAT cases correspond to concrete
counter-examples found.  (Note that there is a small
discrepancy between the sum of the times of the individual steps and
the total time since the latter is the run time for a single script
that performed all the steps in the table, and a small portion of the
script is not taken into account in the recorded times for the steps.)
Additionally, Table \ref{Tab:2} shows the various times it takes to
check whether the individual rotor thrusts $\tilde{f}_j$ stay in
$[f_{min},f_{max}]$ for various rotor failure combinations, in which
failed rotors exert zero thrust, i.e., they experience complete
failures. In particular, note that the octorotor remains safe for up
to two rotor failures, except when the failures occur on adjacent
rotors that rotate the same direction. Lastly, Table \ref{Tab:3} shows
the times it takes to check the individual rotor thrust bounds for
rotor failures where the thrusts of failed rotors can get stuck at
zero or nonzero values. (Note that $\mu_{max}=2$ for all entries
except where $\mu_{max}$ is specified, in which case $\mu_{max}$
was decreased in order to obtain an UNSAT.) All times were obtained
using $16$ cores.

\section{Conclusion}\label{Sec:con}

This paper introduces a framework for the formal verification of the
safety of control systems using exponential barrier functions, with
the target application of ensuring the faithful command tracking of an
octorotor controller. Our method uses barrier functions to construct
candidate invariant sets for the octorotor dynamics, and invariance of
these sets are then formally checked by the SMT solver dReal. We
account for potential rotor failures through a pseudo-inverse control
allocator, which we also verify to produce valid rotor thrusts via
dReal. Using our approach, we verify that a particular controller
causes the octorotor to follow commands within a certain margin of
error under dynamic disturbances and several types of rotor
failures. Our approach is fairly general, and can potentially be used
to construct safe invariant sets in many different systems and using
different kinds of controllers. One particularly promising application
is in the verification of gain scheduled controllers for aircraft.

\begin{table*}
\caption{Run times for checking safety conditions for no rotor failures}
\label{Tab:1}
\begin{center}
{\footnotesize
\begin{tabular}{ | l | l | l | l | l | l | }
\hline
\multirow{2}{*}{Step} & \multirow{2}{*}{Precision} & \multirow{2}{*}{SAT/UNSAT} & \multicolumn{3}{|c|}{Proof times (with 16 cores)} \\ \cline{4-6}
& & & real & user & sys \\ \hline
Search space contains b.f. support & $10^{-2}$ & UNSAT & 0m0.015s & \multicolumn{1}{|c|}{-} & \multicolumn{1}{|c|}{-} \\ \hline
Invariance, $v_z$ b.f. components & $10^{-2}$ & UNSAT & 0m0.165s & \multicolumn{1}{|c|}{-} & \multicolumn{1}{|c|}{-} \\ \hline
Invariance, $\phi,\Omega_1$ b.f. components & $10^{-2}$ & UNSAT & 253m20.942s & \multicolumn{1}{|c|}{-} & \multicolumn{1}{|c|}{-} \\ \hline
Invariance, $\theta,\Omega_2$ b.f. components & $10^{-2}$ & UNSAT & 253m41.323s & \multicolumn{1}{|c|}{-} & \multicolumn{1}{|c|}{-} \\ \hline
Invariance, $\psi,\Omega_3$ b.f. components & $10^{-2}$ & UNSAT & 0m0.055s & \multicolumn{1}{|c|}{-} & \multicolumn{1}{|c|}{-} \\ \hline
Rotor bounds (no rotor failures) & $10^{-2}$ & UNSAT & 0m0.109s & \multicolumn{1}{|c|}{-} & \multicolumn{1}{|c|}{-} \\ \hline
Total time & & & 507m3.938s & 7576m43.915s & 443m29.505s \\ \hline
\end{tabular}}
\end{center}
\end{table*}

\begin{table*}
\caption{Run times for checking rotor bounds with complete rotor failures}
\label{Tab:2}
\begin{center}
{\footnotesize
\begin{tabular}{ | l | l | l | l | l | l | }
\hline
\multirow{2}{*}{Step: Rotor bounds under failures} & \multirow{2}{*}{Precision} & \multirow{2}{*}{SAT/UNSAT} & \multicolumn{3}{|c|}{Proof times (with 16 cores)} \\ \cline{4-6}
& & & real & user & sys \\ \hline
Rotor $1$ complete failure & $10^{-2}$ & UNSAT & 0m0.925s & 0m5.847s & 0m4.571s \\ \hline
Rotor $1$, $2$ complete failures & $10^{-5}$  & SAT & 0m0.818s & 0m4.810s & 0m4.687s \\ \hline
Rotor $1$, $3$ complete failures  & $10^{-2}$ & UNSAT & 0m0.923s & 0m5.823s & 0m4.594s \\ \hline
Rotor $1$, $4$ complete failures  & $10^{-2}$ & UNSAT & 0m0.910s & 0m5.870s & 0m4.416s \\ \hline
Rotor $1$, $5$ complete failures  & $10^{-2}$ & UNSAT & 0m0.803s & 0m4.179s & 0m3.971s \\ \hline
Rotor $1$, $6$ complete failures  & $10^{-2}$ & UNSAT & 0m0.898s & 0m5.660s & 0m4.343s \\ \hline
Rotor $1$, $7$ complete failures  & $10^{-2}$ & UNSAT & 0m0.893s & 0m5.723s & 0m4.313s \\ \hline
Rotor $1$, $8$ complete failures & $10^{-2}$ & UNSAT & 0m1.035s & 0m6.969s & 0m4.322s \\ \hline
\end{tabular}}
\end{center}
\end{table*}

\begin{table*}
\caption{Run times for checking rotor bounds with stuck rotor failures}
\label{Tab:3}
\begin{center}
{\footnotesize
\begin{tabular}{ | l | l | l | l | l | l | }
\hline
\multirow{2}{*}{Step: Rotor bounds under failures} & \multirow{2}{*}{Precision} & \multirow{2}{*}{SAT/UNSAT} & \multicolumn{3}{|c|}{Proof times (with 16 cores)} \\ \cline{4-6}
& & & real & user & sys \\ \hline
\centering{Rotor $1$ thrust stuck at $\frac{mg}{8}$} & $10^{-2}$ & UNSAT & 0m0.946s & 0m5.245s & 0m4.499s \\ \hline
Rotor $1$ thrust stuck at $\frac{mg}{6}$, $\mu_{max}=1.6$ & $10^{-2}$  & UNSAT & 0m2.267s & 0m25.905s & 0m4.608s \\ \hline
Rotor $1$, $2$ thrusts stuck at $0$, $\frac{mg}{8}$ & $10^{-2}$ & SAT & 0m0.932s & 0m3.660s & 0m3.644s \\ \hline
Rotor $1$, $2$ thrusts stuck at $\frac{mg}{8}$, $\frac{mg}{8}$, $\mu_{max}=1.5$ & $10^{-2}$ & UNSAT & 0m7.261s & 1m45.881s & 0m5.297s \\ \hline
Rotor $1$, $2$ thrusts stuck at $0$, $\frac{mg}{6}$  & $10^{-5}$ & SAT & 0m0.887s & 0m5.530s & 0m4.155s \\ \hline
Rotor $1$, $2$ thrusts stuck at $\frac{mg}{6}$, $\frac{mg}{6}$  & $10^{-4}$ & SAT & 0m0.883s & 0m3.833s & 0m3.820s \\ \hline
Rotor $1$, $3$ thrusts stuck at $0$, $\frac{mg}{8}$, $\mu_{max}=1.6$ & $10^{-2}$ & UNSAT & 0m10.498s & 2m37.854s & 0m4.288s \\ \hline
Rotor $1$, $3$ thrusts stuck at $\frac{mg}{8}$, $\frac{mg}{8}$  & $10^{-2}$ & UNSAT & 0m1.270s & 0m9.844s & 0m3.996s \\ \hline
Rotor $1$, $3$ thrusts stuck at $0$, $\frac{mg}{6}$  & $10^{-6}$ & SAT & 0m0.863s & 0m3.947s & 0m4.064s \\ \hline
Rotor $1$, $3$ thrusts stuck at $\frac{mg}{6}$, $\frac{mg}{6}$, $\mu_{max}=1.3$ & $10^{-2}$ & UNSAT & 0m22.850s & 5m51.932s & 0m6.821s \\ \hline
Rotor $1$, $8$ thrusts stuck at $0$, $\frac{mg}{8}$  & $10^{-2}$ & UNSAT & 0m1.037s & 0m6.810s & 0m4.310s \\ \hline
Rotor $1$, $8$ thrusts stuck at $\frac{mg}{8}$, $\frac{mg}{8}$  & $10^{-2}$ & UNSAT & 0m3.899s & 0m51.320s & 0m4.599s \\ \hline
Rotor $1$, $8$ thrusts stuck at $0$, $\frac{mg}{6}$  & $10^{-2}$ & UNSAT & 0m1.502s & 0m13.834s & 0m4.068s \\ \hline
Rotor $1$, $8$ thrusts stuck at $\frac{mg}{6}$, $\frac{mg}{6}$, $\mu_{max}=1.1$  & $10^{-2}$ & UNSAT & 0m8.274s & 2m0.787s & 0m5.040s \\ \hline
\end{tabular}}
\end{center}
\end{table*}

\bibliographystyle{plain}


\subsubsection*{References}}
\usepackage{algorithm2e}
\usepackage{algorithmic}

\usepackage{todonotes}
\usepackage{comment}

\usepackage{tikz}

\theoremstyle{definition}
\newtheorem{theorem}{Theorem}
\newtheorem{claim}{Claim}

\usepackage{hyperref}

\newcommand{\arabic{algorithm}}{\arabic{algorithm}}

\begin{document}

\title{Learning Large DAGs by Combining Continuous Optimization and Feedback Arc Set Heuristics}

\author{Pierre Gillot and Pekka Parviainen\\University of Bergen, Norway}
\date{}
\maketitle

\begin{abstract}
Bayesian networks represent relations between variables using a directed acyclic graph (DAG). Learning the DAG is an NP-hard problem and exact learning algorithms are feasible only for small sets of variables. We propose two scalable heuristics for learning DAGs in the linear structural equation case. Our methods learn the DAG by alternating between unconstrained gradient descent-based step to optimize an objective function and solving a maximum acyclic subgraph problem to enforce acyclicity. Thanks to this decoupling, our methods scale up beyond thousands of variables.
\end{abstract}

\section{Introduction}

Bayesian networks are probabilistic graphical models that represent joint distributions among random variables. They consist of a structure which is a directed acyclic graph (DAG) representing conditional independencies and parameters that specify local conditional distributions.

Bayesian networks can handle both discrete and continuous variables. In this work, we concentrate on continuous variables. Specifically, we study linear structural equation models (SEMs) where the local conditional distribution in a node is a Gaussian whose mean is linear combination of the values of its parents.

Traditionally, there are two main approaches for learning DAGs. In {\em constraint-based} approach (see, e.g., \cite{pearl00, spirtes00}), one performs conditional independence tests and tries to construct a DAG that expresses the same conditional independencies as the test results. We take the {\em score-based} approach (see, e.g., \cite{Cooper1992, Heckerman1995}) where one tries to find a DAG that maximizes a score. Typically, one uses decomposable scores, that is, the score of a DAG is a sum of local scores for each node-parent set pair. This leads to a combinatorial optimization problem where one picks a parent set for each node while satisfying the constraint that the resulting graph is acyclic.

The combinatorial learning problem is NP-hard \cite{chickering96} and developing scalable methods is challenging. Indeed, state-of-the-art exact learning methods scale only up to few hundred nodes \cite{cussens11} and scalable algorithms for SEMs rely on approaches such as local modifications \cite{Aragam_2019}. A recent breakthrough, NOTEARS \cite{zheng2018dags} circumvents the combinatorial problem by formulating a continuous acyclicity constraint. This enables usage of gradient-based optimization methods. The bottleneck with respect to scalability lies in the cubic complexity for the calculation of the acyclicity function which involves the computation of a matrix exponential. GOLEM \cite{NEURIPS2020_d04d42cd} is similar to NOTEARS but replaces the generalized LASSO objective found in NOTEARS by a log-likelihood-based fitness function. It shares however the same computational bottleneck as NOTEARS due to the acyclicity constraint. Some methods circumvent this bottleneck by finding a sparse graph without the acyclicity constraint and impose acyclicity afterwards \cite{DBLP:journals/corr/abs-2006-03005, yu20}.

Our goal is to develop a fast heuristic for learning linear DAGs. In other words, we trade accuracy for speed. We speed-up learning by decoupling the optimization of the objective function from the acyclicity constraint in a similar fashion as \cite{DBLP:journals/jmlr/ParkK17}\footnote{Differences are discussed in Section~\ref{sec:method}.}. At a general level, we learn by iteratively repeating the following steps:
\begin{enumerate}
    \item Given an acyclic graph, find a graph (possibly cyclic) which is better in terms of the objective function value.
    \item Given a cyclic graph, find an acyclic graph.
\end{enumerate}

The first step can be solved efficiently using state-of-the-art gradient-based solvers. We present two variants for this step. {\em ProxiMAS} uses proximal gradient descent whereas {\em OptiMAS} uses standard automatic differentiation and gradient-based updates.

In the second step, the cyclic solution from the first step is converted into an acyclic one. The quality
of the final solution depends crucially on the quality of this conversion. We solve an instance of maximum acyclic subgraph (MAS) problem which has been previously used to learn DAGs \cite{pmlr-v138-gillot20a}. Intuitively, we prefer keeping arcs whose weights are far from zero. Solving the MAS problem exactly is NP-hard but there exists efficient heuristics for solving its complement, the feedback arc set (FAS) problem \cite{10.14778/3021924.3021930}.

Our experiments (Section~\ref{sec:exp}) show that our methods can quickly find reasonable solutions on datasets with thousands variables, even on modest running time. OptiMAS and ProxiMAS perform well compared to GD, NOTEARS and GOLEM in large-scale learning when resources are limited in terms of processors and memory.

\section{Background}

\subsection{Linear structural equation models and Bayesian network structure learning} \label{sec:sem}

A Bayesian network is a representation of a joint probability distribution. It consists of two parts: a structure and parameters. The structure of a Bayesian network is a DAG $G=(N, A)$ where $N$ is the node set and $A$ is the directed adjacency matrix; we denote the parent set of a node $v$ by ${Pa}_v$. Parameters specify local conditional distributions $ P(v\mid {Pa}_v)$ and the joint probability distribution is represented by a factorization
\[P(N) = \prod_{v\in N} P(v\mid {Pa}_v).\]

We consider linear structural equation models (SEMs) where local conditional distributions are Gaussian distributions whose mean is a linear combination of the values of the parents of the variable.
The structure of a linear SEM is determined by a weight matrix $W\in \mathbb{R}^{d\!\times\!d}$; $W(i, j)$ is non-zero if and only if $A(i,j)=1$, that is, there is an arc from $i$ to $j$. For a $d$-dimensional data vector $x$, we have
\[
x = xW + e,
\]
where $e$ is a $d$-dimensional noise vector. The elements of $e$ are independent of each other.

The goal in Bayesian network structure learning is to find a DAG $G$ that fits the data. We are given a data matrix $X\in \mathbb{R}^{n\!\times\!d}$ with $n$ samples of $d$-dimensional vectors. Our goal is to find a weight matrix $W\in \mathbb{R}^{d\!\times\!d}$ that represents an acyclic graph. To quantify how well the DAG and the weights fit the data, we can use the least-squares loss $\frac{1}{2n}\|XW-X\|_2^2$. Furthermore, we want to induce sparsity and therefore we add a regularization term $g(W)$. This leads to the following optimization problem.
\begin{equation} \label{eq:sem}
\begin{aligned}
\underset{W}{\text{argmin}} \text{ } & \frac{1}{2n}\|XW-X\|_2^2 + \lambda_1 g (W)\\
\text{s.t.} \quad & W \text{ is acyclic.}
\end{aligned}
\end{equation}
In the above formulation, $\lambda_1$ is a user-defined constant that determines the strength of regularization. To induce sparsity, we regularize with $L1$-norm, that is,  $g(W)= \|W\|_1 = \sum_{i,j} |W(i, j)|$.

\subsection{Maximum acyclic subgraph and feedback arc set} \label{sec:mas}

Formally, the maximum acyclic subgraph (MAS) problem is defined as follows. We are given a directed graph $G = (V, E)$ and a weight function $w(e)$ that assigns a weight for each arc $e\in E$. The goal is to find an acyclic graph $G' = (V, E')$ such that $E'\subseteq E$ and $\sum_{e\in E'} w(e)$ is maximized.

The maximum acyclic subgraph problem has a dual (or complementary) problem:  the feedback arc set (FAS) problem. In FAS, we are given a directed graph $G = (V, E)$ and a weight function $w(e)$ just like in MAS. The goal is to find an arc set $E''$ such that $G'' = (V, E\setminus E'')$ is acyclic and $\sum_{e\in E''} w(e)$ is minimized. It is well-known that $E' = E\setminus E''$. Thus, MAS can be solved by first solving FAS and then performing a simple subtraction of sets to obtain the solution to MAS.

Both MAS and FAS are NP-hard \cite{Kar72}. Therefore, exact algorithms are intractable on large graphs. Fortunately, there exists fast heuristics for FAS \cite{berger90, EADES1993319, 10.14778/3021924.3021930}.

\section{Proposed method} \label{sec:method}

A critical difficulty in solving Equation \ref{eq:sem} stems from a combination of two problems:
\begin{itemize}
    \item The quadratic objective function for the linear SEM problem has at most $nd^3$ quadratic terms. Indeed, the quadratic expression $\|XW-X\|_2^2$ is a sum of $n\!\times\!d$ squared expressions $\big((XW)_{i,j}-X_{i,j}\big)^2$, where each $(XW)_{i,j}$ is a linear expression consisting of $d$ terms. As $X$ is a continuous data matrix, one can rarely simplify the quadratic objective function significantly.
    \item Enforcing acyclicity. Standard constraints lead to NP-hard combinatorial problems. In the continuous setting, a smooth function exists that encodes acyclicity but with a prohibitive cubic complexity.
\end{itemize}
The main contribution of this work therefore is to address these two concerns. First, we decompose the quadratic optimization problem into a sequence of easier subproblems using iterative optimization techniques. Second, we separate entirely the quadratic optimization from the acyclicity constraints. Acyclicity is enforced by solving a MAS task as a proxy. The outline of the proposed method is shown in Algorithm \ref{alg:proximas} which iteratively does the following steps at each iteration $k$:
\begin{enumerate}
    \item A new objective function is created based on the acyclic solution $W_{k-1}$ obtained at the end of the previous iteration, which penalizes the original linear SEM objective by the least-square term $\frac{\lambda_2}{2} \|W - W_{k-1} \|_2^2$.
    \item An optimization step is performed on the MAS-penalized problem, leading to a new cyclic solution $\widetilde{W}_k$.
    \item An acyclic projection $W_k$ of the previously obtained cyclic solution $\widetilde{W}_k$ is extracted, based on the squared values of $\widetilde{W}_k$. Formally, we attempt to compute $W_k = \widetilde{W}_k \odot A_k$, where $A_k$ is the solution of the following MAS problem:
\begin{equation} \label{eq:MAS}
\begin{aligned}
A_k &= \underset{A}{\text{ argmax }} {\displaystyle\sum_{i,j}}|\widetilde{W}_k(i,j)|^2 A(i,j)\\
&\quad\quad\quad\text{s.t.} \quad A \in \{0,\!1\}^{d\!\times\!d} \text{ is acyclic.}
\end{aligned}
\end{equation}
    As mentioned in Section~\ref{sec:mas}, finding optimal solutions for MAS is usually too time consuming and one has to resort to heuristics. We use a vectorized version of the approximation algorithm by Eades \cite{EADES1993319} to find the acyclic weighted adjacency matrix $W_k$ (Algorithm \ref{alg:greedy}).
\end{enumerate}
Intuitively, steps 1-2 are designed such that the updated weights matrix $\widetilde{W}_k$ will be constrained to remain in the vicinity of the previously found acyclic solution $W_{k-1}$ returned by the MAS heuristic.  In that sense, we approximate the acyclicity function used in NOTEARS and GOLEM by a projection term toward acyclic solutions which is much easier to compute and differentiate. Step 3 aims to preserve edges that represent the most important dependencies. In other words, we want to keep the weights that are far from zero. Solving MAS using weights that are squares of the original weights is equivalent to minimizing $\|W_k - \widetilde{W}_k\|_2$, which corresponds to finding the acyclic solutions that are closest to the cyclic solutions returned by the iterative optimization process. By alternating between optimization steps and MAS extractions via the repetition of steps 1-3, we aim to navigate through the search space of the original linear SEM problem by ``following the trail'' of a sequence of dynamically generated acyclic solutions.

The GD algorithm introduced in \cite{DBLP:journals/jmlr/ParkK17} follows a similar strategy. It proceeds by repeating the three following steps: 1) make a gradient step for the linear SEM objective, 2) project the current cyclic solution to its MAS solution and 3) fit the linear SEM problem constrained by the newly found acyclic structure; as an optional fourth step, when the progress is too small they resort to a swapping order heuristic. The main difference compared to our method is that we do not perform steps 3 and 4. The GD algorithm is greedier than our method, since we do not attempt to optimize the parameters of every discovered acyclic structure. From a practical standpoint, at each step of the GD algorithm, $d$ LASSO instances have to be solved which becomes intractable for large-scale structure learning. In comparison, by directly plugging in step 1 the MAS projections to the linear SEM objective as dynamically evolving penalization terms, our approach circumvents entirely the need to solve these LASSO instances.

\subsection{Connection with online convex optimization}

Perhaps surprisingly, the dynamic nature of the proposed optimization procedure is not particularly challenging to work with in practice. Algorithm \ref{alg:proximas} can, indeed, be seen as a special case of an online convex optimization (OCO) problem. In his seminal paper \cite{10.5555/3041838.3041955}, Zinkevich introduces this framework which he defines as such:
\begin{itemize}
    \item $F \subset \mathbb{R}^n$ is a feasible set (assumed bounded, closed and non-empty).
    \item $(c_k)_k$ is an infinite sequence of smooth convex functions from $F$ to $\mathbb{R}$, with bounded gradients.
    \item At each step $k$, an element $x_k \in F$ is selected, then assigned the cost $c_k(x_k)$.
\end{itemize}
In OCO, the standard optimization error becomes ill-defined and one seeks to optimize instead the so-called regret defined as
\begin{equation*}
    \text{regret } = \text{ } \sum_{k\leq K} c_k(x_k) - \min_{x \in F} \sum_{k \leq K} c_k(x).
\end{equation*}
Zinkevich  was the first to extend the gradient descent algorithm to its online form. It is well known that assuming convexity of the $c_k$ and boundedness of the gradients, online gradient descent achieves $\mathcal{O}(\sqrt{K})$ regret bound and this bound is improved to $\mathcal{O}\big(\log(K)\big)$ assuming strong convexity of the $c_k$ \cite{DBLP:journals/corr/abs-1909-05207}. More general classes of OCO algorithms have been studied \cite{NIPS2009_ec5aa0b7, DBLP:journals/corr/abs-1810-03594}, notably (accelerated) proximal gradient descent algorithms concerned about composite convex functions of the form $\phi_k = f_k + g$ where only the $f_k$ are smooth. Improved regret bounds again hold assuming strong convexity of the $\phi_k$.

The proposed method is therefore theoretically well-behaved: by considering the functions $f_k: W \mapsto \frac{1}{2n}\|XW-X\|_2^2 + \frac{\lambda_2}{2} \|W - W_{k-1} \|_2^2$ and $\phi_k: W \mapsto f_k(W) + \lambda_1 g(W)$ (Algorithm \ref{alg:proximas} line 2), notice that every $\phi_k$ is $\lambda_2$-strongly convex since for each $k$, the function $W \mapsto \phi_k(W) - \tfrac{\lambda_2}{2}\|W - W_{k-1} \|_2^2=\frac{1}{2n}\|XW-X\|_2^2 + \lambda_1 g(W)$ is convex; Algorithm \ref{alg:proximas} therefore inherits aforementioned regret bounds from the OCO setting.

\subsection{Implementation details}

We implemented two variants of the proposed method:
\begin{itemize}
    \item The first implementation, {\bf ProxiMAS}, is designed to take full advantage of the properties of the objective functions $\phi_k$, owing to the decoupling with acyclicity. Recall that we have $f_k: W \mapsto \frac{1}{2n}\|XW-X\|_2^2 + \frac{\lambda_2}{2} \|W - W_{k-1} \|_2^2$ and $\phi_k: W \mapsto f_k(W) + \lambda_1 g(W)$, where $g$ is convex and the $f_k$ are smooth and convex. By smooth, we mean that every $f_k$ is differentiable with its gradient defined as $\nabla f_k(W) = \tfrac{1}{n}X^t X (W - I) + \lambda_2 (W-W_{k-1})$ and using the Cauchy-Schwarz inequality, one can easily show that every $f_k$ has a Lipschitz-continuous gradient with optimal Lipschitz constant $L_k$ upper bounded by $L = \tfrac{1}{n}\|X^tX + n\lambda_2I\|_2$, a value that does not depend on $k$. We can therefore use a proximal gradient descent optimization scheme as a backbone for our implementation, hence the name ProxiMAS. In practice, we use the FISTA algorithm \cite{FISTA}, an accelerated proximal algorithm with $\mathcal{O}(\tfrac{1}{k^2})$ convergence rate (in an offline optimization setting). One should notice that the running time of ProxiMAS does not depend on the number of samples $n$, since the proximal updates depend only on the covariance matrix $X^tX \in \mathbb{R}^{d\!\times\!d}$ which can be pre-computed.
    \item The second implementation, {\bf OptiMAS}, replaces the proximal gradient descent by gradient descent-like steps. The main interest in doing so is that automatic differentiation will handle the optimization using a generic gradient descent-based solver. Despite the linear SEM objective being non-differentiable when the regularization term is the $L1$ norm, automatic differentiation frameworks can in practice optimize such non-smooth objective. Thus, OptiMAS is agnostic to the choice of the optimizer. In principle, one can use any variant of gradient-based optimizers. In our implementation, we have used Adam \cite{DBLP:journals/corr/KingmaB14}.
\end{itemize}
We stress that both variants are taking full advantage of vectorization and are thus GPU accelerated, first because we lifted the need to solve a sequence of LASSO instances at each step, second because the MAS heuristic (Algorithm \ref{alg:greedy}) is efficiently vectorized and runs quasi-linearly with respect to the number of nodes $d$ when a GPU is available.

\begin{figure

\begin{minipage}[t]{0.47\textwidth}
\begin{algorithm}[H]
	\caption{Proposed method}
	\label{alg:proximas}
	\begin{center}
	\begin{algorithmic}[1]
	    \REQUIRE Data $X \!\in\! \mathbb{R}^{n\!\times\!d}$, initialization $W_0 \!\in\! \mathbb{R}^{d\!\times\!d}$, number of iterations $K$, $\lambda_1>0$, $\lambda_2 > 0$, optimizer
	    \ENSURE Approximate solution to Equation~\ref{eq:sem}
		\FOR {$1 \leq k \leq K$}
		    \STATE Define $f_k: W \mapsto \frac{1}{2n}\|XW-X\|_2^2 + \frac{\lambda_2}{2} \|W - W_{k-1} \|_2^2$ $\phi_k: W \mapsto f_k(W) + \lambda_1 g(W)$
		    \STATE Make an optimization step on $\phi_k$: $\widetilde{W}_k = \text{step}(\phi_k, \text{ optimizer})$
		    \STATE Project updated weights to their MAS approximation: $W_k = \text{greedy\_MAS}\big(\widetilde{W}_k\big)$
	    \ENDFOR
	    \STATE {\bfseries return} $W_{\text{best}} = \underset{k}{\text{ argmin }} \tfrac{1}{2n} \|XW_k-X\|_2^2 + \lambda_1 g(W_k)$
	    \texttt{\vspace{0.3cm}}
	\end{algorithmic}
	\end{center}
\end{algorithm}
\end{minipage}
\hfill
\begin{minipage}[t]{0.47\textwidth}
\begin{algorithm}[H]
	\caption{Vectorized greedy MAS}
	\label{alg:greedy}
	\begin{center}
	\begin{algorithmic}[1]
	    \REQUIRE $\widetilde{W} \!\in\! \mathbb{R}^{d\!\times\!d}$
	    \ENSURE Approximate solution to Equation \ref{eq:MAS}
        \STATE $\hat{W} = \widetilde{W} \odot \widetilde{W}$
        \STATE $\text{scores} = \hat{W}\!.\text{sum}(\text{dim}\!=\!0)$
        \STATE $\text{order} = \text{zeros}(\text{size}\!=\!d)$
        \STATE $\text{ub} = (d+1)\times \max(\text{scores})$
		\FOR {$0 \leq i < d$}
		    \STATE $\text{node} = \text{ argmin }\text{scores}$
		    \STATE $\text{order}[-(i+1)] = \text{node}$
		    \STATE $\text{scores}[\text{node}] = \text{ub}$
		    \STATE $\text{scores} = \text{scores} - \hat{W}[\text{node},:]$
	    \ENDFOR
	    \STATE $\text{order}^{-1} = \text{ argsort }\text{order}$
	    \STATE $W = \text{lower\_triangular}\big(\widetilde{W}[\text{order},\text{order}]\big)$
	    \STATE $W = W[\text{order}^{-1},\text{order}^{-1}]$
	    \STATE {\bfseries return} $W$
	\end{algorithmic}
	\end{center}
\end{algorithm}
\end{minipage}

\end{figure}

\section{Experiments} \label{sec:exp}

We now present our experimental pipeline. We choose to compare the proposed algorithms (ProxiMAS and OptiMAS) against an iterative method (GD \cite{DBLP:journals/jmlr/ParkK17}) and the current state-of-the-art methods for sparse linear SEM structure recovery (NOTEARS \cite{zheng2018dags} and GOLEM
\cite{NEURIPS2020_d04d42cd}).

\subsection{Data generation}

We adopt a similar setup as in \cite{zheng2018dags, NEURIPS2020_d04d42cd}: we first generate random DAGs based on Erd{\H o}s-R\'enyi ("ER") and scale-free ("SF") models. We consider three sparsity regimes: sampled DAGs have $k\times d$ edges, where $d$ is the number of nodes and $k \in \{1, 2, 4\}$. Graphs are denoted by "ERk" or "SFk" depending on their graph model and sparsity. Then, we generate the weighted adjacency matrices $W$ by assigning random weights uniformly sampled in the range $[-2,-0.5] \cup [0.5, 2]$. Finally, we generate samples $X$ following the linear SEM model: $X=E (I-W)^{-1}$, where $E \in \mathbb{R}^{n\!\times\!d}$ represents $n$ i.i.d. samples from either a Gaussian, exponential or Gumbel distribution in $\mathbb{R}^d$. For all distributions, we investigate both the equal variance ("EV") setting with scale $1.0$ for all variables and the non-equal variance ("NV") setting where every variable has its scale sampled uniformly in the range $[0.5, 1.5]$. Unless stated otherwise, $n$ samples are generated both for the training data and for the validation data, with $n \in \{1000, 10000\}$.

\subsection{Metrics}

In order to evaluate the performance of the different methods, we compute the false negative rate (FNR), false positive rate (FPR) and the normalized structural Hamming distance (SHD) between predicted and groundtruth adjacency matrices. We proceed similarly with the undirected adjacency matrices. The Gaussian negative log-likelihood is also computed on the validation data (unseen during training). Aforementioned metrics are extracted at different thresholding values of the predicted weights matrices. Different methods behave differently at a fixed thresholding value. For example, we observed in our large-scale tests with limited running time that, for any fixed threshold, OptiMAS tend to produce significantly sparser graphs than NOTEARS and GOLEM. Thus, OptiMAS has lower FPR and higher FNR. In order to get a general performance score independent from the choice of a thresholding value, we additionally consider the average precision score as implemented in the scikit-learn package. This metric is robust against strong class imbalance as it occurs in large-scale sparse structure recovery. For brevity, only a fraction of the figures are shown in this paper.

\subsection{Implementation} \label{sec:implementation}

The two proposed methods (ProxiMAS and OptiMAS) are implemented in pytorch 1.8. The GOLEM method comes in two variations GOLEM-EV and GOLEM-NV originally implemented in tensorflow. In order to streamline benchmarking these variations were re-implemented in pytorch. The tensorflow and pytorch implementations were compared at fixed seed and produce nearly identical results given the same data; speedwise, we found the difference between the two implementations to be insignificant for large scale graphs with thousands of nodes. The original implementation of NOTEARS relies on an L-BFGS-B solver implemented in scipy and as mentioned in \cite{NEURIPS2020_d04d42cd} it does not scale to large instances with thousands of variables, thus for fairness we re-implemented it in pytorch as well. The existing implementation of the GD algorithm is written in R and uses the highly optimized package glmnet, thus we did not alter the implementation. All methods have full GPU support, with the exception of GD which relies on the LASSO implementation from the glmnet package and is restricted to CPU.  For equal comparison, all methods are tested in a multi-threaded setting, without GPU.

\subsection{Hyperparameters}

All the tested methods have a hyperparameter $\lambda_1$  that regulates sparsity. An additional hyperparameter $\lambda_2$ exists for the GOLEM, NOTEARS ProxiMAS and OptiMAS methods which enforces "dagness". The values of these two hyperparameters yield different behavior depending on the method. The chosen value of $\lambda_1$ for NOTEARS and the chosen values of $\lambda_1$ and $\lambda_2$ for GOLEM are those recommended by their authors. The original NOTEARS implementation is based on an augmented Lagrangian method and does not use the $\lambda_2$ hyperparameter. We added this hyperparameter to our pytorch implementation of NOTEARS the same way as in the GOLEM implementation. We do not claim to have performed any model selection, but chose values that worked well in our tests. The table below gives the values of the hyperparameters depending on the method:
\begin{center}
\begin{table}[!htbp]
\caption{\label{tab:hyperparam}Sparsity and dagness hyperparameters for each method}
\begin{tabular}{ |c|c|c|c|c|c|c| }
\hline
 & OptiMAS & ProxiMAS & NOTEARS \cite{zheng2018dags} & GOLEM-EV \cite{NEURIPS2020_d04d42cd} & GOLEM-NV \cite{NEURIPS2020_d04d42cd} & GD \cite{DBLP:journals/jmlr/ParkK17} \\
\hline
$\lambda_1$ & $0.1$ & $0.1$ & $0.1$ & $0.02$ & $0.002$ & $0.1$ \\
$\lambda_2$ & $20.0$ & $20.0$ & $5.0$ & $5.0$ & $5.0$ & -\\
\hline
\end{tabular}
\end{table}
\end{center}
Additionally, in all experiments the ProxiMAS and OptiMAS methods are configured to start enforcing acyclicity after $50$ minutes of solving time, whereas NOTEARS and GOLEM enforce acyclicity the entire time as in their original papers. Finally, the methods that rely on automatic differentiation (NOTEARS, GOLEM and OptiMAS) all use the Adam optimizer \cite{DBLP:journals/corr/KingmaB14} as a backbone for gradient descent, with default learning rate $0.001$ as in \cite{NEURIPS2020_d04d42cd}.

\subsection{Benchmarking pipeline}\label{sec:pipeline}

The experiments were run on a cluster with Intel Xeon-Gold 6138 2.0 GHz / 6230 2.1 GHz CPUs. The number of cores and amount of memory used in each experiments are shown in Table~\ref{tab:param}.

We present three different experiments to emphasize the advantageous scaling of the proposed methods comparatively to the state of the art. Table~\ref{tab:param} below shows the different instances we ran for every experiment.

\begin{center}
\begin{table}[!tbph]
\caption{\label{tab:param}Parameters for each experiment}
\begin{tabular}{ |l|l|l|l| }
\hline
 & Experiment 1 & Experiment 2 & Experiment 3\\
\hline
 $d$ & $1000, 5000$ & $5000$ & $5000, 10000, 15000, 20000$\\
 $k$ & $1, 2, 4$ & $1$ & $1$ \\
 $n$ & $1000, 10000$ & $1000, 10000$ & $1000, 10000$\\
 Graph type & ER, SF & ER & ER\\
 Noise type & Gaussian, exponential, Gumbel & Gaussian & Gaussian\\
 Scale type & EV, NV & EV, NV & EV\\
 Repetitions & $10$ & $1$ & $1$ \\
 Total instances & $1440$ & $4$ & $8$ \\
 CPU cores & $4$ & $4$ & $32$\\
 Memory (GB) & $16$ & $16$ & $128$ \\
 Runtime (h) & $1$ & $24$ & $1$\\
 \hline
\end{tabular}
\end{table}
\end{center}

\subsection{Results}

In the first experiment (Experiment~1 in Table~\ref{tab:param}), we generated data with different noise models. We show selected results in Figure~\ref{fig:exp1}.

\begin{center}
\begin{figure}[!htpb]
\begin{tikzpicture}[scale=1]
\node[inner sep=0pt] at (0,0)
    {\includegraphics[scale=0.45]{figs_exp1/exp1_gaussian_EV_average_precision.png}};
\node[inner sep=0pt] at (7.25,0)
    {\includegraphics[scale=0.45]{figs_exp1/exp1_gaussian_NV_average_precision.png}};
\node at (-0.1,3) {Gaussian-EV};
\node at (7.15,3) {Gaussian-NV};
\node[inner sep=0pt] at (0,-6.25)
    {\includegraphics[scale=0.45]{figs_exp1/exp1_exponential_EV_average_precision.png}};
 \node[inner sep=0pt] at (7.25,-6.25)
    {\includegraphics[scale=0.45]{figs_exp1/exp1_gumbel_EV_average_precision.png}};
\node[inner sep=0pt] at (3.625,-9.25)
    {\includegraphics[scale=0.6]{fig/legend_bars.png}};
 \node at (-0.1,-3.25) {Exponential-EV};
 \node at (7.15,-3.25) {Gumbel-EV};
\end{tikzpicture}
\caption{\label{fig:exp1}Average precisions for different noise distributions and data set sizes. $d$= number of nodes, $n$= number of samples. Confidence intervals show the standard deviation over $10$ data sets.}
\end{figure}
\end{center}

Overall, OptiMAS and ProxiMAS overperform the benchmark methods. Generally, GD performs equally good as OptiMAS and ProxiMAS when $d=1000$ or $n=1000$. However, it becomes slow when data sets grow. Especially, GD usually fails to even find a solution when there are lots of samples ($n=10000$). Most of the time, NOTEARS and GOLEM are in par with GD or slightly better. Comparing OptiMAS and ProxiMAS, we notice that their performance is usually similar to each other. The main difference is that OptiMAS performs significantly better on more complex graphs (ER4) with $d=1000$; We suspect that this difference is due to numerical instabilities.

We also wanted to analyze how the available running time affected each method. Therefore, we generated two data sets, one with $1000$ samples and the other with $10000$ samples, from an ER1 model with Gaussian-EV noise (Experiment~2 in Table~\ref{tab:param}) and let all methods run for 24 hours. We recorded a snapshot of the weight matrix $W$ every hour. Average precisions are shown in Figure~\ref{fig:exp2}.

\begin{figure}[!htpb]
\begin{center}
\begin{tikzpicture}[scale=1]
\node[inner sep=0pt] at (0,-4.5)
    {\includegraphics[scale=0.6]{fig/exp2_ER_gaussian_EV_average_precision.png}};
\node[inner sep=0pt] at (0,-7)
    {\includegraphics[scale=0.6]{fig/legend_lines.png}};
\end{tikzpicture}
\caption{\label{fig:exp2}Average precision measured at different time points. Data generated from ER1 with $5000$ nodes and Gaussian-EV noise. Note that on the plot on the left hand side, the curves for OptiMAS and ProxiMAS are overlapping.}
\end{center}
\end{figure}

We observe that both OptiMAS and ProxiMAS find good solutions quickly. However, improvement after the first hour is negligible. GD starts slowly but eventually catches up with OptiMAS and ProxiMAS and often ends up with a slightly better solution. As in Experiment~1, we notice that the scalability of GD suffers from having lots of observations. Initially, NOTEARS is far behind but keeps improving significantly afterwards and after 24 hours it has found a solution that is almost as good as the ones found by OptiMAS and ProxiMAS. GOLEM performs similarly with NOTEARS when there are $10000$ samples but struggles with $1000$ samples.

Next, we study the scalability of the different methods. To this end, we varied the number of nodes between $5000$ and $20000$ and generated either $1000$ or $10000$ data samples from an ER1 model with Gaussian-EV noise (Experiment~3 in Table~\ref{tab:param}). All methods were given 1 hour running time. Average precisions are shown in Figure~\ref{fig:exp3}.
\begin{center}
\begin{figure}[!htpb]
\begin{tikzpicture}[scale=1]
\node[inner sep=0pt] at (0,0)
    {\includegraphics[scale=0.6]{fig/exp3_ER_gaussian_EV_average_precision.png}};
\node[inner sep=0pt] at (0,-2.5)
    {\includegraphics[scale=0.6]{fig/legend_lines.png}};
\end{tikzpicture}
\caption{\label{fig:exp3}Scalability of different methods. Average precision is measured with different number of nodes $d$. Data was generated from ER1 with Gaussian-EV noise.}
\end{figure}
\end{center}
We notice that with $1000$ samples the average precision for OptiMAS is high for all data set sizes and decreases only little when the number of nodes grows. However, with $10000$ samples average precision drops faster when the number of nodes grows. This may seem counter-intuitive as one would expect that more observations would lead to better performance. The likely explanation for this behavior is that, due to the fixed running time, OptiMAS performed fewer iterations and this countered the effect of increase number of observations at this scale. We can contrast this behavior to ProxiMAS whose running time does not depend on the number of observations. With $n=1000$, ProxiMAS starts with nearly as high average precision as OptiMAS but its performance deteriorates quickly after $10000$ nodes. However, with $n=10000$, the drop is less significant and ProxiMAS clearly outperforms OptiMAS when there are $10000$ or more nodes. We also notice that GOLEM and NOTEARS struggle to learn anything within an hour when there are more than $10000$ nodes. GD performs better than GOLEM and NOTEARS when $n=1000$ but when $n=10000$ it only finds a solution for $d=5000$.

Table~\ref{tab:mem_time} shows memory usage of the different methods and their time per iteration (with acyclicity enforced for NOTEARS, GOLEM, OptiMAS and ProxiMAS). We see that OptiMAS is clearly most memory-efficient. ProxiMAS uses more memory than OptiMAS but much less than NOTEARS and GOLEM. The memory consumption of GD is larger than ProxiMAS but smaller than NOTEARS and GOLEM. Time per iteration for GD is very inconsistent due to the order swapping heuristic it uses at certain iterations, thus we omitted it; as a rule, we observed that GD scales unfavorably with respect to both number of nodes and number of samples.

\begin{center}
\begin{table}[!htpb]
\caption{\label{tab:mem_time}Estimation of the memory usage and time per iteration (32 cores, ER1, Gaussian-EV, n=1000)}
\begin{tabular}{ |l|c|c|c|c||c|c|c|c| }
\hline
 & \multicolumn{4}{|c||}{Memory (GB)} & \multicolumn{4}{|c|}{Time/iteration (s)} \\
 \cline{1-9}
 $d$ & $5000$ & $10000$ & $15000$ & $20000$ & $5000$ & $10000$ & $15000$ & $20000$ \\
\hline
 OptiMAS & $1$ & $6$ & $13$ & $24$ & $1$ & $2$ & $4$ & $6$ \\
 ProxiMAS & $1$ & $7$ & $14$ & $25$ & $1$ & $3$ & $7$ & $13$ \\
 NOTEARS \cite{zheng2018dags} & $6$ & $23$ & $52$ & $92$ & $6$ & $40$ & $100$ & $250$ \\
 GOLEM \cite{NEURIPS2020_d04d42cd} & $6$ & $23$ & $53$ & $94$ & $6$ & $45$ & $120$ & $280$ \\
 GD \cite{DBLP:journals/jmlr/ParkK17} & $4$ & $12$ & $27$ & $47$ & $-$ & $-$ & $-$ & $-$ \\
 \hline
 \end{tabular}
 \end{table}
 \end{center}

\section{Discussion} \label{sec:discussion}

We presented two different heuristics (ProxiMAS and OptiMAS) for the structure recovery problem in the linear SEM case, revolving around a decoupling of the acyclicity constraints from the continuous optimization itself. We observed that both methods have excellent scaling (both space and time). OptiMAS scales particularly well when the number of samples $n$ is small. On the contrary, ProxiMAS has invariant scaling with respect to $n$ and scales in practice better than OptiMAS when the number of samples is large.

In our observations, both ProxiMAS and OptiMAS tend to get stuck on local extremum: the sequence of acyclic DAGs returned by the two methods is conditioned by the initial cyclic solution provided to them. This drawback can be alleviated by ``warm-starting'': run the algorithm initially without the MAS penalization and extraction steps (Algorithm \ref{alg:proximas} lines 2 and 4), then add these steps at some point during the execution. This strategy is made viable since a single MAS step is enough to guarantee acyclicity. Our experiments show that in practice, very good DAGs can be found even when most of the running time is dedicated to fitting the model without enforcing acyclicity.

Based on our experiments, OptiMAS and ProxiMAS are most competitive in situations where there is a large number of nodes and limited amount of computational resources. If there are a couple of thousands of nodes or less, the current state of the art is preferred. Similarly, if one can afford to run GD, NOTEARS or GOLEM for a long enough period of time, these algorithms will eventually outperform ProxiMAS/OptiMAS. However, in such a situation one could use OptiMAS or ProxiMAS to find an initial solution and use it to ``warm-start'' GD, NOTEARS or GOLEM.

Another limitation of our methods is that it is unclear at the moment how the theoretical results from online convex optimization translate with respect to the original problem. Currently, we are not aware of any necessary condition for the convergence of the proposed methods. This opens up an avenue of future research: Can we prove anything about the quality of the solutions? Can we say something for a specific type of data? Does the fact that we use a heuristic to find a maximum acyclic subgraph has an effect and would improving MAS also translate in better structure learning?

\subsection*{Acknowledgements}

The computations were performed on resources provided by
UNINETT Sigma2 - the National Infrastructure for High Performance Computing and
Data Storage in Norway.

The authors thank Young Woong Park for providing the R code for GD.

\bibliographystyle{plain}