# Hamilton-Jacobi Reachability Analysis for Hybrid Systems with Controlled and
Forced Transitions

Javier Borquez1, Shuang Peng1, Yiyu Chen1, Quan Nguyen1, Somil Bansal1
1Authors are with the University of Southern California, LA, USA.

###### Abstract

Hybrid dynamical systems with non-linear dynamics are one of the most general
modeling tools for representing robotic systems, especially contact-rich
systems. However, providing guarantees regarding the safety or performance of
such hybrid systems can still prove to be a challenging problem because it
requires simultaneous reasoning about continuous state evolution and discrete
mode switching. In this work, we address this problem by extending classical
Hamilton-Jacobi (HJ) reachability analysis, a formal verification method for
continuous non-linear dynamics in the presence of bounded inputs and
disturbances, to hybrid dynamical systems. Our framework can compute reachable
sets for hybrid systems consisting of multiple discrete modes, each with its
own set of non-linear continuous dynamics, discrete transitions that can be
directly commanded or forced by a discrete control input, while still
accounting for control bounds and adversarial disturbances in the state
evolution. Along with the reachable set, the proposed framework also provides
an optimal continuous and discrete controller to ensure system safety. We
demonstrate our framework in simulation on an aircraft collision avoidance
problem, as well as on a real-world testbed to solve the optimal mode planning
problem for a quadruped with multiple gaits.

##  I Introduction

Hybrid dynamical system is a popular and versatile tool to model robotic
systems that exhibit both continuous and discrete dynamics [1, 2]. For
example, for a legged robot, the swinging of a leg has continuous dynamics,
whereas contacts with ground are well-modeled as discrete events. However, as
with many advanced modeling tools, their very complexity and richness presents
challenges, particularly when it comes to ensuring safety and performance for
such systems. Designing safe controllers for hybrid dynamical systems demands
simultaneous reasoning about the continuous evolution of states and the
discrete mode transitions, a task that can quickly become computationally
intensive and conceptually challenging.

Reachability analysis is a powerful and effective approach for analyzing and
controlling hybrid dynamical systems. It provides a comprehensive
understanding of the system’s behavior by characterizing the Backward
Reachable Tube (BRT) of the system – the set of all possible initial states
that can eventually reach some set of target states under optimal control
actions. If the target set represents the undesirable states for the system,
the BRT represents the set of states that are potentially unsafe for the
system and should be avoided. The converse of the BRT thus provides a safe
operation region for the system. Reachability analysis for hybrid dynamical
systems has been studied using a variety of different frameworks, such as
zonotopes [3, 4, 5, 6], hybrid denotable sets [7, 8], Taylor models[9, 10,
11], satisfiability modulo theory [12, 13], support functions [14, 15, 16],
timed automata [17] and linear hybrid automata [17, 18, 19, 20], among others.
However, these methods typically impose restrictive assumptions on the system,
such as limiting the analysis to linear continuous dynamics, not accounting
for discrete control switches, bounding the number of mode changes, not
allowing continuous controls, or not considering disturbances and dynamics
uncertainty. Moreover, most of this methods only provide an over-approximation
of the BRT.

Another approach for computing BRT for hybrid dynamical systems is via
Hamilton-Jacobi (HJ) Reachability analysis. Its advantages include
compatibility with general non-linear system dynamics, formal treatment of
bounded disturbances, and the ability to deal with state and input constraints
[21]. Several works have addressed the control and safety analysis of hybrid
dynamical systems through HJ Reachability [22, 23, 24, 25, 26]. However, these
methods often rely on hand and hard-coded heuristics to account for discrete
mode switching. Moreover, the BRT computation relies on an iterative
refinement that may not necessarily converge. One particular work for the
computation of region-of-attraction and stabilizing controllers for walking
robots has at its core a generalization of the HJ reachability framework to
account for discontinuous state changes originating from state resets [27].
However, generalizing the HJ reachability framework for multiple discrete
modes and accounting for discrete controls still remains challenging.

In this work, we propose a general HJ reachability framework for hybrid
dynamical systems. Our framework formulates the BRT computation for hybrid
systems as a robust optimal control problem. Using the principle of dynamic
programming, we simultaneously reason about optimal continuous and discrete
control inputs for this optimal control problem, resulting in a generalized
version of Hamilton-Jacobi-Bellman equation. We demonstrate that the existing
numerical algorithms to compute BRTs can easily be extended to compute this
generalized value function, thereby providing both the BRT as well as the safe
discrete and continuous control inputs for the hybrid system. Our framework
can handle hybrid systems consisting of multiple discrete modes with nonlinear
dynamics, discrete control switches that allow or force transitions between
different discrete modes, (discontinuous) state resets upon discrete switches,
control bounds, as well as, potentially adversarial, disturbances. We
demonstrate the efficacy of our framework in simulation and on a real-world
quadruped testbed.

##  II Problem Formulation

Consider a hybrid dynamical system with a finite set of discrete modes
$Q=\\{q_{1},q_{2},\ldots,q_{N}\\}$, connected through controlled or forced
switches. We also refer to $q_{i}$ as the discrete state of the system.

Figure 1: Hybrid dynamical system with controlled and forced transitions.

Let $x\in X\subset\mathbb{R}^{n_{x}}$ be the continuous state of the system,
$u\in U\subset\mathbb{R}^{n_{u}}$ be the continuous control input, and $d\in
D\subset\mathbb{R}^{n_{d}}$ be the continuous disturbance in the system. In
each discrete mode $q_{i}\in Q$, the continuous state evolves according to the
dynamics: $\dot{x}=f_{i}(x,u,d)$, as long as $x\in S_{i}\subset X$. Here,
$S_{i}$ is the valid operation domain for mode $q_{i}$. $q_{i}$ also has
discrete control switches $\sigma_{ij}$ that allow controlled transition into
another discrete mode $q_{j}$, where
$j\in\\{j_{1},j_{2},\ldots,j_{n}\\}\subset\\{1,2,\ldots,N\\}$. Note that the
number of discrete modes the system can transition into (that is, $n$) may
vary across different $q_{i}$. As $x$ approaches ${S_{i}}^{C}$ (denoted
$x\to{S_{i}}^{C}$), the system must take one of the forced control switches,
$\varsigma_{ik}$. This leads to a forced transition into another discrete
modes $q_{k}$ where
$k\in\\{k_{1},k_{2},...,k_{m}\\}\subset\\{1,2,\ldots,N\\}$. Each discrete
transition from $q_{i}$ to $q_{j}$, whether controlled or forced, might lead
to a state reset $x^{+}=R_{ij}(x)$, where $x$ is the state before the
transition, $x^{+}$ is the state after the transition.

Our key objective in this work is to compute the Backward Reachable Tube (BRT)
of this hybrid dynamical system, defined as the set of initial discrete modes
and continuous states $(x,q)$ of the system for which the agent acting
optimally and under worst-case disturbances will eventually reach a target set
$\mathcal{L}$ within the time horizon $[t,T]$. Mathematically, the BRT is
defined as:

$\mathcal{V}(t)=\\{(x_{0},q_{0}):\forall d(\cdot),\exists
u(\cdot),\sigma(\cdot),\varsigma(\cdot),\exists
s\in[t,T],x(s)\in\mathcal{L}\\},$

where $x(s)$ denotes the (continuous) state of the system at time $s$,
starting from state $x_{0}$ and discrete mode $q_{0}$ under continuous control
profile $u(\cdot)$, discrete control switch profile $\sigma(\cdot)$111If there
is no discrete control being leveraged at time $s$, $\sigma(s)$ or
$\varsigma(s)$ can be thought of a “dummy” discrete control that keeps the
system in the current discrete mode., forced control switch profile
$\varsigma(\cdot)$, and disturbance $d(\cdot)$. Along with the BRT, we are
interested in finding the optimal discrete and continuous control inputs that
steer the system to the target set.

Conversely, if $\mathcal{L}$ represents the set of undesirable states for the
system (such as obstacles for a navigation robot), we are interested in
computing the BRT with the role of control and disturbance switched, i.e.,
finding the set of initial states from which getting into $\mathcal{L}$ is
unavoidable, despite best control effort.

Running example (Dog1D): To illustrate our approach, we will use the following
running example throughout this paper. We consider very simplified, high-level
dynamics for a robot quadruped moving in one dimension in a room with an
obstacle (table). Figure 2: A simple 1D longitudinal navigation scenario. The
robot goal is to reach the marked area on the other side of the table as shown
in Fig. 2. The hybrid system representation of this system is shown in Fig. 3:
Figure 3: Dog1D hybrid system formulation. Here, the continuous state
$x\in\mathbb{R}$ is the position of the robot and $q_{i}$ with
$i\in\\{1,2,3\\}$ are the discrete modes of the system. $u\in[0,1]$ is the
continuous control of the system. The target set is given as the set of
positions $9\leq x\leq 10$. In this example, $q_{1}$ represents a walking
gait, $q_{3}$ represents frozen dynamics when the walking robot comes in
contact with the table obstacle, and $q_{2}$ is a crawling gait that allows
the robot to move under the table but at a slower rate compared to walking.
The BRT for this example corresponds to all the combinations of continuous
states and discrete modes $(x,q)$ from which the quadruped can reach the
target set within a $5$ second time horizon.

##  III  Background: Hamilton-Jacobi Reachability

One way to compute the BRT for continuous dynamical systems is through
Hamilton-Jacobi (HJ) reachability analysis. To illustrate HJ reachability, we
assume that the system has only one discrete mode, i.e., there are no discrete
transitions.

In HJ reachability, the BRT computation is formulated as a zero-sum game
between control and disturbance. This results in a robust optimal control
problem that can be solved using the dynamic programming principle. First, a
target function $l(x)$ is defined whose sub-zero level set is the target set
$\mathcal{L}$, i.e. $\mathcal{L}=\\{x:l(x)\leq 0\\}$. The BRT seeks to find
all states that could enter $\mathcal{L}$ at any point within the time
horizon. This is computed by finding the minimum distance to $\mathcal{L}$
over time:

$J(x,t,u(\cdot),d(\cdot))=\min_{\tau\in[t,T]}l\left(x(\tau)\right)$ (1)

Our goal is to capture this minimum distance for optimal system trajectories.
Thus, we compute the optimal control that minimizes this distance (drives the
system towards the target) and the worst-case disturbance signal that
maximizes the distance. The value function corresponding to this robust
optimal control problem is:

$V(x,t)=\adjustlimits{\sup}_{d(\cdot)}{\inf}_{u(\cdot)}\\{J(x,t,u(\cdot),d(\cdot))\\}.$
(2)

The value function in (2) can be computed using dynamic programming, which
results in the following final value Hamilton-Jacobi-Isaacs Variational
Inequality (HJI-VI):

$\begin{array}[]{c}\min\left\\{D_{t}V(x,t)+\mathcal{H}(x,t),l(x)-V(x,t)\right\\}=0\\\
V(x,T)=l(x)\end{array}$ (3)

$D_{t}$ and $\nabla$ represent the time and spatial gradients of the value
function. $\mathcal{H}$ is the Hamiltonian, which optimizes over the inner
product between the spatial gradients of the value function and the dynamics
to compute the optimal control and disturbance:

$\mathcal{H}(x,t)=\max_{d\in D}\min_{u\in U}\nabla V(x,t)\cdot f(x,u,d).$ (4)

The term $l(x)-V(x,t)$ in (3) restricts system trajectories that enter and
leave the target set, enforcing that any trajectory with a negative distance
at any time will continue to have a negative distance for the rest of the time
horizon. Once the value function is obtained, the BRT is given as the sub-zero
level set of the value function:

$\mathcal{V}(t)=\\{x:V(x,t)\leq 0\\}$ (5)

The corresponding optimal control can be derived as:

$u^{*}(x,t)=\operatorname*{arg\,min}_{u\in U}\max_{d\in D}\nabla V(x,t)\cdot
f(x,u,d)$ (6)

The system can guarantee reaching the target set as long as it starts inside
the BRT and applies the optimal control in (6) at the BRT boundary. The
optimal adversarial disturbance can be similarly obtained as (6).

This formulation can be used to provide safety guarantees by switching the
roles of the control and disturbance in (4). In that case, the BRT represents
the states that will be driven into the target by optimal disturbance even if
we apply optimal control, so safety can be maintained as long as the system
stays outside the BRT and applies optimal control at the boundary.

However, one of the key limitations of HJ reachability analysis is that it
assumes continuous dynamics and does not account for discrete mode switching
or state resets. Overcoming this limitation is the core focus of this work.

Running example (Dog1D): Computing the BRT for the running example using the
classical Hamilton-Jacobi formulation with no optimal discrete mode switches
gives the results shown in Fig. 4. For the crawling gait it can be observed
that the BRT (the blue shaded area) propagates through the obstacle and the
BRT includes $x=4m$ because the system can move at $1m/s$ reaching the
boundary of the target in exactly $5s$. For the walking wait, even though the
system is allowed to move faster, the BRT stops at the boundary of the
obstacle as the walking gait gets stuck, so only states that start from the
right side of the obstacle can reach the target set. If we imagine the system
being able to transition optimally between the discrete modes, it would only
crawl for states in the obstacle and walk elsewhere, allowing to reach states
farther to the left. As we demonstrate later, the proposed framework can
reason about such optimal transitions. Figure 4:  (Top) Value function for
crawling gait mode. (Bottom) Value function for walking gait mode. Obstacle is
shown in shaded red. Blue shade shows BRT as the subzero level of the value
function. Green shade shows target set as the subzero level of $l(x)$.

##  IV HJ Reachability for General Hybrid Systems

The core contribution of this work is Theorem 1, an extension of the classical
HJ reachability framework to hybrid dynamical systems with controlled and
forced transitions, as well as state resets.

###### Theorem 1

For a hybrid dynamical system, the value function $V(x,q_{i},t)$ that captures
the BRT for an implicit target function $l(x)$ is given by solving the
following constrained VI. If $x\in S_{i}$:

$\min\\{l(x)-V(x,q_{i},t),\min_{\sigma_{ij}}V(R_{ij}(x),q_{j},t)-V(x,q_{i},t),\\\
D_{t}V(x,q_{i},t)+\max_{d\in D}\min_{u\in U}\nabla V(x,q_{i},t)\cdot
f_{i}(x,u,d)\\}=0,$

and if $x\in S_{i}^{C}$,
$V(x,q_{i},t)=\min_{\varsigma_{ik}}V(R_{ik}(x),q_{k},t)$, with terminal time
condition:

$V(x,q_{i},T)=l(x)$ (7)

Intuitively, Theorem 1 updates the value function for each discrete mode
$q_{i}$ with the optimal value across the discrete modes it can transition to.
For states inside the operation domain $S_{i}$, this corresponds to taking one
of the control switches $\sigma_{ij}$ with their corresponding reset maps
$R_{ij}(x)$ or staying in the operation mode $q_{i}$ itself, which will lead
to a continuous flow of the value function, similar to (3). In contrast, for
states outside the operation mode domain $\delta S_{i}^{C}$, the system has to
take a forced switch $\varsigma_{ik}$ with their corresponding reset maps
$R_{ik}(x)$. Thus, the optimal value function is given by the one that takes
the system to the minimum value upon the discrete transition. The proof of
Theorem 1 formalizes this intuition using the Bellman principle of optimality,
the detail of the proof can be found on the Appendix. Once we compute the
value function, the BRT for the hybrid system is given as:

$\mathcal{V}(t)=\\{(x,q):V(x,q,t)\leq 0\\}$ (8)

Finally, the optimal discrete and continuous controls at state $(x,q_{i})$ at
time $t$ are given as:

$\displaystyle u^{*}(x,q_{i},t)$ $\displaystyle=\operatorname*{arg\,min}_{u\in
U}\max_{d\in D}\nabla V(x,q_{i},t)\cdot f_{i}(x,u,d)$ (9)
$\displaystyle\sigma_{ij^{*}}(x,q_{i},t)$
$\displaystyle=\operatorname*{arg\,min}_{\sigma_{ij}}V(R_{ij}(x),q_{j},t)\text{~{}if~{}}x\in
S_{i}$ $\displaystyle\varsigma_{ik^{*}}(x,q_{i},t)$
$\displaystyle=\operatorname*{arg\,min}_{\varsigma_{ik}}V(R_{ik}(x),q_{k},t)\text{~{}if~{}}x\in
S_{i}^{C},$

where, for completeness, we add a “dummy” discrete control $\sigma_{ii}$ that
keeps the system in the same discrete mode.

###### Remark 1

Note that the proposed framework can easily be extended to scenarios where the
forced transitions cannot be controlled and represent adversarial or uncertain
transitions instead. In this case, we can use $\max$ instead of $\min$ over
$\varsigma_{ik}$ in Theorem 1 to account for the worst-case behavior.

Numerical implementation: We now present a numerical algorithm that can be
used to calculate the value function in Theorem 1. It builds upon the value
function calculation for the classical HJI-VI in (2), which is solved using
currently available level set methods [28]. Specifically, the value function
is computed over a discretized state-space grid and is propagated in time
using a small timestep $\delta$. After each propagation step, the value
function is updated for all $(x,q)$ using Theorem 1, until the time horizon
$T$ is reached.

Input: $l(x),T,\delta$

Output: $V(x,q_{i},t)\forall i$

Initialization: $V(x,q_{i},T)=l(x)\forall i$; $t=T$

while  _$(t >0)$_ do

foreach  _$(q_{i}\in Q)$_ do

$V(x,q_{i},t-\delta)=V(x,q_{i},t)+\newline \hskip
50.00008pt\displaystyle\max_{d}\min_{u}\nabla V(x,q_{i},t)\cdot
f_{i}(x,u,d)\delta$

foreach  _$(q_{i}\in Q)$_ do

if  _$(x\in S_{i}$)_ then

$\displaystyle V(x,q_{i},t-\delta)=\min\\{l(x),V(x,q_{i},t-\delta),\newline
\hskip 68.00012pt\min_{\sigma_{ij}}V(R_{ij}(x),q_{j},t-\delta)\\}$

else if  _$(x\in S_{i}^{C}$))_ then

$\displaystyle
V(x,q_{i},t-\delta)=\min_{\varsigma_{ik}}V(R_{ik}(x),q_{k},t-\delta)$

$t=t-\delta$

Algorithm 1 Value function computation for hybrid dynamical system

It is important to stress the remarkable similarity of this new hybrid
reachability algorithm to its continuous counterpart in (2). Indeed, there are
only two key differences: (1) we now propagate $N$ value functions
simultaneously (corresponding to $N$ discrete modes), as opposed to just one
value function. (2) We “adjust” the value function for each discrete mode to
account for discrete switches (the second for-loop in the algorithm). As a
result, the proposed hybrid reachability algorithm leads to a quadratic (in
$N$) growth in the worst-case computational complexity of its continuous
system counterpart.

##  V Case Studies

###  V-A Running example (Dog1D):

We now apply the proposed method to compute the BRT for the running example.
To implement Algorithm 1, we use a modified version of helperOC library and
the level set toolbox, ToolboxLS[28], both of which are used to solve
classical HJI-VI. We use a grid of 301 points over the $x$ dimension for each
of the 3 discrete operation modes. The overall calculation takes 12.125s
running on an Intel Core i5-6200U CPU @ 2.30GHz. Slices of the BRT starting in
mode $q_{1}$ (the walking gait) for different time horizons are shown in Fig.
5:

Figure 5: Value functions for states starting in a walking gait for a $3$ and
$5$ seconds time horizon. Area within the obstacle is shown in shaded red.
Blue shade shows BRT as the subzero level of the value function. Green shade
shows the target set as the subzero level of the implicit target function
$l(x)$.

The intuitive solution for the optimal decision-making in this scenario is to
use the walking gait everywhere except in the obstacle area. The fast walking
gait allows the system to reach the target quicker, while crawling allows to
expand the reachable states beyond the obstacle. The computed solution using
our algorithm indeed aligns with this intuition. For example, if we consider
the top value function for the 3-second time horizon in Fig. 5, we can observe
that the limit of the BRT is $7m$ away from the target set, which coincides
with the intuitive solution of walking for $1s$ covering $3m$, then crawling
under the table for $1s$ covering its $1m$ width, to finally go back to
walking for $1s$ covering other $3m$.

The bottom value function in Fig. 5 shows the BRT corresponding to a $5s$ time
horizon. Compared to Fig. 4, where the BRT stopped growing beyond the obstacle
boundary, we can see how the proposed algorithm can optimally leverage
discrete transitions along with continuous control to reach the target set
from a wider set of initial states.

###  V-B Two-Aircraft Conflict Resolution:

We next consider the two-aircraft conflict resolution example presented in
[26]. Here, two aircraft flying at a fixed altitude are considered. Each
aircraft controls its speed; we assume that the first aircraft is
controllable, while the other aircraft’s speed is uncertain (within a known
interval) and hence modeled as an adversarial disturbance to ensure safety.
Using relative coordinates, the continuous dynamics of the system are
described by:

$\dot{x}_{r}=-u+d\cos\psi_{r}+\omega_{u}y_{r},~{}\dot{y}_{r}=d\sin\psi_{r}-\omega_{u}x_{r},~{}\dot{\psi}_{r}=\omega_{d}-\omega_{u},$

where state $[x_{r},y_{r},\varphi_{r}]^{T}$ is the relative $xy$ position and
heading between the aircraft. $u$ and $d$ are the velocities of the two
aircraft, and $\omega_{u}$ and $\omega_{d}$ are their turning rates.

The conflict resolution protocol consists of three different operation modes,
shown in Fig. 6. The aircraft begin in mode $q_{1}$ with a straight flight
($\omega_{u}=\omega_{d}=0$), keeping a costant relative heading. In this mode,
the aircraft are on a collision course. At some time, the aircraft can begin
the collision avoidance maneuver by taking the switch $\sigma_{12}$,
transitioning the system into mode $q_{2}$ and an instantaneous heading change
of $\pi/2$. In mode $q_{2}$, the aircraft undergo a circular flight path,
where both aircraft use the same constant turning rate
($\omega_{u}=\omega_{d}=1$). Once the aircraft undergo a heading change of
$\pi$ radians, the aircraft are forced to switch to mode $q_{3}$, making
another instantaneous $\pi/2$ turn and resuming their straight flight
($\omega_{u}=\omega_{d}=0$). The two aircraft are now on a collision-free
path. To keep track of the transition time on configuration $q_{2}$ we add an
additional timer state $z$.

Figure 6: Hybrid formulation for the two-aircraft conflict resolution
protocol.

For this system, we are interested in finding the set of all initial states in
mode $q_{1}$ from which a collision cannot be averted between the two aircraft
under the above protocol, despite the best control effort. Thus, the
complement of the BRT will represent the safe states for the system. For BRT
calculations, we consider a circular collision target set of 5 units around
the controlled aircraft. This corresponds to the two aircraft being in close
proximity of each other, also referred to as “loss of separation”. We use an
initial relative heading of $\varphi_{r}=2\pi/3$ (shown in Fig. 6) and a time
horizon of $t=5s$. For aircraft velocities we consider $u\in[1.5,3]$ and
$d\in[2,4]$. The calculations are carried on a $[x_{r},y_{r},z]$ grid of
$[201,201,101]$ points, $\varphi_{r}$ has null dynamics in all operation modes
and is not considered in the grid.

Results for BRT starting from operation mode $q_{1}$ are shown in Fig.7. They
can be interpreted as follows: if the relative coordinates of the aircraft
belong to a state outside the BRT, the collision avoidance maneuver can be
initiated while ensuring a collision can always be avoided as long as the
optimal discrete and continuous control is applied.

Figure 7: Backward reachable tube for two-aircraft conflict resolution
starting from mode $q_{1}$. (Left) BRT for fixed speed of aircraft. (Right)
BRT for aircraft under optimal control and disturbance.

The left BRT exactly replicates the results of the case presented in [26]
where both aircraft keep their speed constant at their maximum value $u=3$ and
$d=4$, where the inside of the BRT correspond to states where the adversarial
airplane can place itself behind the controlled aircraft and use its higher
velocity to force a collision. However, unlike [26], we do not use any hard-
coded heuristics to account for the discrete transitions during the BRT
computation. Furthermore, the right BRT considers the scenario where, in
addition to the control on the transition, both aircraft control their
velocities optimally to avoid/force a collision with the other aircraft. The
changes in the BRT are reflected in the shaded areas, where the growth in the
BRT close to the target corresponds to cases where the adversarial aircraft
slows down to align itself with the controlled aircraft and then use its
higher velocity to force a collision. The decrease in size near the tail
corresponds to states where the blue airplane slows down to avoid a collision.

##  VI Hardware Experiments

Figure 8: (a) Trajectories of the quadruped safely navigating through the
slope or narrow ground obstacles to reach the goal area (pink cylinder). The
goal area is independent of the $z$ position, so the robot can reach it via
slope or ground route. (b) Overlaid trajectories in the real-world obstacle
setup. (c, d) First-person view along the slope route (c) and the ground route
(d).

We next apply our method for task-based, high-level mode planning on a real-
world quadruped robot, to reach a goal position in a terrain consisting of
various obstacles. The quadruped has different walking modes, such as normal
walking, walking on a slope, etc., each of which is modeled as a discrete mode
in the hybrid system (see Fig. 9).

Figure 9: Hybrid control modes for the quadruped system.

Within each discrete mode, we use simplified continuous dynamics to describe
center-of-mass evolution:

$\dot{x}=k_{x}V_{N}\sin\varphi+d_{x};\quad\dot{y}=k_{y}V_{N}\cos\varphi+d_{y};\quad\dot{\varphi}=k_{\omega}\omega_{N},$

where the state $[x,y,\varphi]^{T}$ represent the $xy$ position and heading
angle of the robot. $V_{N}=0.25m/s$ is the fixed (normalized) linear velocity
and $\omega_{N}\in[-0.5,0.5]rad/s$ is the angular velocity control. $k_{x}$,
$k_{y}$, and $k_{\omega}$ are velocity gains for different operation modes.
$d_{x}$ and $d_{y}$ are bounded disturbances on $xy$ velocity. Also, we use
discrete state $[z,\theta,\xi]^{T}$ to set the body’s height, pitch, and roll
angle in different operation modes for overcoming various obstacles.

Each operation mode matches a specific obstacle or terrain shown in Fig. 8(b).
The quadruped begins in mode $q_{1}$ for fast walking. If the robot reaches
the boundary of the slope area, the system will be forced to switch to mode
$q_{2}$ with a $20^{\circ}$ body pitch angle for slope climbing. At some time,
the control may switch to mode $q_{3}$ to walk slowly with a lower body
height, this allows the robot to make a narrow turn or crawl under obstacles.
If the robot encounters a tilted obstacle, the system will go through a forced
switch into mode $q_{4}$ to walk with a tilted body. Whenever the robot
touches a ground obstacle, the system will stay in mode $q_{5}$ permanently
with frozen dynamics. The quadruped needs to reason about optimally switching
between these different walking modes in order to reach its goal area (the
$xy$ area marked as the pink cylinder in Fig. 8(a, b)) as quickly as possible
without colliding with any obstacles.

BRT for this experiment is calculated on a $[x,y,\varphi]$ grid of
$[40,40,72]$ points, over a time horizon of $t=45s$. Having the BRT, we can
acquire optimal velocity control, operation mode, and the desired discrete
body states at each time. The quadruped is equipped with an Intel T265
tracking camera, allowing it to estimate its global state via visual-inertial
odometry. The high-level optimal velocity control and discrete body states
provided by our framework are then tracked by a low-level MPC controller that
uses the centroidal dynamics of the quadruped [29].

Our experiment results are shown in Fig. 8 and can also be seen in the
accompanying video. In our environment setup, the optimal (fastest) route to
reach the target is via leveraging the slope on the right. The corresponding
trajectory is shown in orange in Fig. 8(a) and (b). Specifically, the robot
remains in the normal walking mode (blue boxes in Fig. 8(a)) and switches to
the slope walking mode once near the slope (magenta boxes). When the robot
approaches the end of the slope, it needs to make a tight left turn to avoid a
collision with the wall ahead. Our framework is able to recognize that and
makes a transition to the slow walking mode to allow for a tighter turn (green
boxes). Once the turn is complete, the robot goes back in the faster walking
mode. Some first-person RGB images along the robot’s path are in Fig. 8(c).

In our second experiment, we put some papers on the slope, making it more
slippery. This is encoded by adding a higher disturbance in the slope mode
$q_{2}$. The new BRT makes the slope path infeasible and hence the robot needs
to reach the target via the (slower) ground route. The new robot trajectory is
shown in light green in Fig. 8(a) and (b). Once again, apart from selecting a
new safe route, the proposed framework is able to switch between different
walking modes along the route to handle tight turns and slanted obstacles, as
shown by the activation of the tilt walking mode (cyan boxes in Fig. 8(a)).

Finally, to test the closed-loop robustness of our method, we add human
disturbance in the ground route experiment. The robot is kicked and dragged by
a human during the locomotion process, as shown by red arrows in Fig. 10.
Specifically, the robot was dragged to face backward, forcing it near the
tilted obstacle. However, the reachability controller is able to ensure safety
by reactivating the tilt walking mode (Fig. 10(b)), demonstrating that the
proposed framework is able to reason about both closed-loop continuous and
discrete control inputs. Nevertheless, it should be emphasized that the
presented algorithm does not solve the entirety of legged locomotion planning,
but rather only provides the top level of a hierarchical planning architecture
that typically include a robust footstep planner, a wholebody controller, and
reliable state estimation [30]. We will further explore integrating our
framework within such architectures in future.

Figure 10: (a) Robot trajectory with human disturbances. (b) Partial
trajectory shows the robot autonomously switches from $q_{3}$ to $q_{4}$ after
being pushed into the tilt obstacle area to ensure safety.

##  VII Discussion And Future Work

We present an extension of the classical HJ reachability framework to hybrid
dynamical systems with controlled and forced transitions, and state resets.
Along with the BRT, the proposed framework provides optimal continuous and
discrete control inputs for the hybrid system. Simulation studies and hardware
experiments demonstrate the proposed method, both to reach goal areas and in
maintaining safety. Our work opens up several exciting future research
directions. First, we rely on grid-based numerical methods to compute the BRT,
whose computational complexity scales exponentially with the number of
continuous states, limiting a direct use of our framework to relatively low-
dimensional systems. We will explore recent advances in learning-based methods
to solve HJI-VI [31, 32] to overcome this challenge. Another interesting
direction would be to extend our framework to the cases involving uncertainty
in the transition surface $S_{i}$ or where controlled switches have different
transition surfaces. Finally, we will apply our framework on a broader range
of robotics applications and systems.

## References

  * [1] R. Goebel, R. G. Sanfelice, and A. R. Teel, “Hybrid dynamical systems,” _IEEE control systems magazine_ , vol. 29, no. 2, pp. 28–93, 2009.
  * [2] A. M. Johnson, S. A. Burden, and D. E. Koditschek, “A hybrid systems model for simple manipulation and self-manipulation systems,” _The International Journal of Robotics Research_ , vol. 35, no. 11, pp. 1354–1392, 2016.
  * [3] M. Althoff, O. Stursberg, and M. Buss, “Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes,” _Nonlinear Analysis: Hybrid Systems_ , vol. 4, no. 2, pp. 233–249, 2010, iFAC World Congress 2008. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1751570X09000442
  * [4] M. Althoff and B. H. Krogh, “Avoiding geometric intersection operations in reachability analysis of hybrid systems,” in _Proceedings of the 15th ACM international conference on Hybrid Systems: Computation and Control_ , 2012, pp. 45–54.
  * [5] M. Maiga, N. Ramdani, L. Travé-Massuyès, and C. Combastel, “A comprehensive method for reachability analysis of uncertain nonlinear hybrid systems,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 9, pp. 2341–2356, 2015.
  * [6] C. Tang and M. Althoff, “Formal verification of robotic contact tasks via reachability analysis,” _arXiv preprint arXiv:2307.13977_ , 2023.
  * [7] L. Benvenuti, D. Bresolin, A. Casagrande, P. Collins, A. Ferrari, E. Mazzi, A. Sangiovanni-Vincentelli, and T. Villa, “Reachability computation for hybrid systems with ariadne,” _IFAC Proceedings Volumes_ , vol. 41, no. 2, pp. 8960–8965, 2008.
  * [8] P. Collins and J. Lygeros, “Computability of finite-time reachable sets for hybrid systems,” in _Proceedings of the 44th IEEE Conference on Decision and Control_ , 2005, pp. 4688–4693.
  * [9] X. Chen, E. Ábrahám, and S. Sankaranarayanan, “Flow*: An analyzer for non-linear hybrid systems,” in _Computer Aided Verification: 25th International Conference, CAV 2013, Saint Petersburg, Russia, July 13-19, 2013\. Proceedings 25_. Springer, 2013, pp. 258–263.
  * [10] X. Chen, E. Abraham, and S. Sankaranarayanan, “Taylor model flowpipe construction for non-linear hybrid systems,” in _2012 IEEE 33rd Real-Time Systems Symposium_. IEEE, 2012, pp. 183–192.
  * [11] N. Kochdumper and M. Althoff, “Reachability analysis for hybrid systems with nonlinear guard sets,” ser. HSCC ’20. Association for Computing Machinery, 2020.
  * [12] S. Kong, S. Gao, W. Chen, and E. Clarke, “dreach: $\delta$-reachability analysis for hybrid systems,” in _Tools and Algorithms for the Construction and Analysis of Systems: 21st International Conference, TACAS 2015, Held as Part of the European Joint Conferences on Theory and Practice of Software, ETAPS 2015, London, UK, April 11-18, 2015, Proceedings 21_. Springer, 2015, pp. 200–205.
  * [13] W. Damm, S. Disch, H. Hungar, S. Jacobs, J. Pang, F. Pigorsch, C. Scholl, U. Waldmann, and B. Wirtz, “Exact state set representations in the verification of linear hybrid systems with large discrete state space,” in _International Symposium on Automated Technology for Verification and Analysis_. Springer, 2007, pp. 425–440.
  * [14] G. Frehse, C. Le Guernic, A. Donzé, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang, and O. Maler, “Spaceex: Scalable verification of hybrid systems,” in _Proc. 23rd International Conference on Computer Aided Verification (CAV)_ , ser. LNCS, S. Q. Ganesh Gopalakrishnan, Ed. Springer, 2011.
  * [15] C. Le Guernic and A. Girard, “Reachability analysis of linear systems using support functions,” _Nonlinear Analysis: Hybrid Systems_ , vol. 4, no. 2, pp. 250–262, 2010.
  * [16] G. Frehse and R. Ray, “Flowpipe-guard intersection for reachability computations with support functions,” _IFAC Proceedings Volumes_ , vol. 45, no. 9, pp. 94–101, 2012.
  * [17] R. Alur, C. Courcoubetis, T. A. Henzinger, and P.-H. Ho, “Hybrid automata: An algorithmic approach to the specification and verification of hybrid systems,” in _International Hybrid Systems Workshop_. Springer, 1991, pp. 209–229.
  * [18] T. A. Henzinger, P.-H. Ho, and H. Wong-Toi, “A user guide to hytech,” in _International Workshop on Tools and Algorithms for the Construction and Analysis of Systems_. Springer, 1995, pp. 41–71.
  * [19] S. Yovine, “Kronos: A verification tool for real-time systems,” _Int. J. Softw. Tools Technol. Transf._ , vol. 1, no. 1-2, pp. 123–133, 1997.
  * [20] O. Maler, A. Pnueli, and J. Sifakis, “On the synthesis of discrete controllers for timed systems,” in _STACS 95: 12th Annual Symposium on Theoretical Aspects of Computer Science Munich, Germany, March 2–4, 1995 Proceedings 12_. Springer, 1995, pp. 229–242.
  * [21] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-Jacobi Reachability: A brief overview and recent advances,” in _IEEE Conference on Decision and Control (CDC)_ , 2017.
  * [22] A. Dhinakaran, M. Chen, G. Chou, J. C. Shih, and C. J. Tomlin, “A hybrid framework for multi-vehicle collision avoidance,” in _Conference on Decision and Control_ , 2017, pp. 2979–2984.
  * [23] J. H. Gillula, G. M. Hoffmann, H. Huang, M. P. Vitus, and C. J. Tomlin, “Applications of hybrid reachability analysis to robotic aerial vehicles,” _The International Journal of Robotics Research_ , vol. 30, no. 3, pp. 335–354, 2011.
  * [24] C. Tomlin, J. Lygeros, and S. Sastry, “Computing controllers for nonlinear hybrid systems,” in _Hybrid Systems: Computation and Control: Second International Workshop, HSCC’99 Berg en Dal, The Netherlands, March 29–31, 1999 Proceedings 2_. Springer, 1999, pp. 238–255.
  * [25] I. Mitchell and C. J. Tomlin, “Level set methods for computation in hybrid systems,” in _International workshop on hybrid systems: Computation and control_. Springer, 2000, pp. 310–323.
  * [26] ——, “Level set methods for computation in hybrid systems,” in _International workshop on hybrid systems: Computation and control_. Springer, 2000, pp. 310–323.
  * [27] J. J. Choi, A. Agrawal, K. Sreenath, C. J. Tomlin, and S. Bansal, “Computation of regions of attraction for hybrid limit cycles using reachability: An application to walking robots,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4504–4511, 2022.
  * [28] I. Mitchell, “A toolbox of level set methods,” _http://www. cs. ubc. ca/mitchell/ToolboxLS/toolboxLS. pdf, Tech. Rep. TR-2004-09_ , 2004.
  * [29] D. Kim, J. D. Carlo, B. Katz, G. Bledt, and S. Kim, “Highly dynamic quadruped locomotion via whole-body impulse control and model predictive control,” 2019\.
  * [30] J. Norby and A. M. Johnson, “Fast global motion planning for dynamic legged robots,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 3829–3836.
  * [31] S. Bansal and C. J. Tomlin, “DeepReach: A deep learning approach to high-dimensional reachability,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2021.
  * [32] A. Lin and S. Bansal, “Generating formal safety assurances for high-dimensional reachability,” 2022. [Online]. Available: https://arxiv.org/abs/2209.12336
  * [33] J. Lygeros, “On reachability and minimum cost optimal control,” _Automatica_ , vol. 40, no. 6, pp. 917–927, 2004.

##  VIII Appendix

###  VIII-A Proof of Theorem 1

We start with an analogous formulation to HJ reachability with a target
function $l(x)$ such that the target set is defined as its sub-zero level set,
$\mathcal{L}=\\{(x):l(x)\leq 0\\}$. The BRT seeks to find all states that
could enter $\mathcal{L}$ at any point within the time horizon. This is
computed by finding the minimum distance to $\mathcal{L}$ over time.

The value function is defined as the cost incurred under optimal discrete and
continuous controls that minimizes distance to the target and optimal
disturbances that maximizes it. We first consider the case when the state
belongs to the interior of the domain $S_{i}$. The value function is given as:

$V(x,q_{i},t)=\max_{d\in D}\min_{u\in
U}\min_{\sigma_{ij}}\min_{s\in[t,T]}l(x(s))$ (10)

With terminal time condition given by:

$V(x,q_{i},T)=l(x)$ (11)

The optimal decision making for discrete transitions is pictorially shown in
the figure below.

The discrete transitions are shown in green arrows and continuous control flow
are indicated by purple trajectories. We also add a “dummy” discrete
transition, $\sigma_{ii}$, that keeps the system in the current discrete mode.
As per (10), the overall value function is given by the minimum cost across
all these possible trajectories of the system. From here, we consider two
different cases:

Case 1: If the system switches to a new discrete mode $j^{*}$, the continuous
state immediately transitions to $R_{ij^{*}}$. In this case, the optimal cost
incurred by the system from the new state $(R_{ij^{*}}(x),q_{j^{*}})$ is, by
definition, given as $V(R_{ij^{*}}(x),q_{j^{*}},t)$. Thus, the optimal cost
across all possible discrete transitions is given as:

$V(x,q_{i},t)=\min_{\sigma_{ij}}V(R_{ij}(x),q_{j},t)$ (12)

Case 2: The system evolves in the same mode $q_{i}$, i.e., it takes the dummy
switch $\sigma_{ii}$. In this case, for a small time step $\delta$ and using
using the dynamic programming principle for the cost function in (10), we
have:

$\displaystyle V(x,q_{i},t)$ $\displaystyle=\max_{d\in D}\min_{u\in
U}\min\\{\min_{s\in[t,t+\delta]}l(x(s)),$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad
V\left(x(t+\delta),q_{i},t+\delta\right)\\}$ $\displaystyle\approx\max_{d\in
D}\min_{u\in U}\min\\{l(x(t)),V\left(x(t+\delta),q_{i},t+\delta\right)\\}$
(13)

Approximating the value function at the next time step with a first order
Taylor expansion:

$\displaystyle V\left(x(t+\delta),q_{i},t+\delta\right)\approx$ $\displaystyle
V(x,q_{i},t)+$ $\displaystyle D_{t}V(x,q_{i},t)\delta+\nabla
V(x,q_{i},t)\cdot\delta x,$ (14)

where the change in the state $\delta x$ can be approximated as
$f_{i}(x,u,d)\delta$. Ignoring the higher order terms and plugging the Taylor
expansion in (VIII-A):

$V(x,q_{i},t)=\min\\{l(x(t)),\\\
V(x,q_{i},t)+D_{t}V(x,q_{i},t)\delta+\max_{d\in D}\min_{u\in U}\nabla
V(x,q_{i},t)\cdot f_{i}(x,u,d)\delta\\}$ (15)

The optimal value function is given by the minimum across the two cases:

$V(x,q_{i},t)=\min\\{l(x),\\\ V(x,q_{i},t)+D_{t}V(x,q_{i},t)\delta+\max_{d\in
D}\min_{u\in U}\nabla V(x,q_{i},t)\cdot f_{i}(x,u,d)\delta,\\\
\min_{\sigma_{ij}}V(R_{ij}(x),q_{ij},t)\\}$ (16)

Subtracting the value function $V(x,q_{i},t)$ from both sides:

$0=\min\\{l(x(t))-V(x,q_{i},t),\\\ \delta(D_{t}V(x,q_{i},t)+\max_{d\in
D}\min_{u\in U}\nabla V(x,q_{i},t)\cdot f_{i}(x,u,d)),\\\
\min_{\sigma_{ij}}V(R_{ij}(x),q_{ij},t)-V(x,q_{i},t)\\}$ (17)

Since the above hold for all $\delta>0$ we must have:

$0=\min\\{l(x(t))-V(x,q_{i},t),\\\ D_{t}V(x,q_{i},t)+\max_{d\in D}\min_{u\in
U}\nabla V(x,q_{i},t)\cdot f_{i}(x,u,d),\\\
\min_{\sigma_{ij}}V(R_{ij}(x),q_{j},t)-V(x,q_{i},t)\\}$ (18)

Note that here we have presented an informal proof based on the Taylor series
expansion of the value function; a more rigorous proof of HJI-VI for
continuous evolution can be found in [33].

For states outside the domain of the discrete mode $q_{i}$ (i.e., $x\in
S_{i}^{C}$), our discrete decision is over the set of available forced
switches $\varsigma_{ik}$. For this case the proof is analogous to the Case 1
previously considered, but without the option to stay in $q_{i}$. The value
function is then given by:

$V(x,q_{i},t)=\min_{\varsigma_{ik}}V(R_{ik}(x),q_{k},t)$ (19)


redrgb1,0,0 MIDL 2020Medical Imaging with Deep Learning 2020 MIDL 2020 – Short
Paper Maxime Chamberland1, Sila Genc1, Erika P. Raven1, Greg D. Parker1, Adam
Cunningham2, Joanne Doherty1,2, Marianne van den Bree2, Chantal M. W. Tax1,
Derek K. Jones1
1 Cardiff University Brain Research Imaging Centre (CUBRIC), Cardiff, United
Kingdom
2 MRC Centre for Neuropsychiatric Genetics and Genomics, Cardiff, United
Kingdom

# Tractometry-based Anomaly Detection for Single-subject White Matter Analysis

###### Abstract

There is an urgent need for a paradigm shift from group-wise comparisons to
individual diagnosis in diffusion MRI (dMRI) to enable the analysis of rare
cases and clinically-heterogeneous groups. Deep autoencoders have shown great
potential to detect anomalies in neuroimaging data. We present a framework
that operates on the manifold of white matter (WM) pathways to learn normative
microstructural features, and discriminate those at genetic risk from controls
in a paediatric population.

###### keywords:

Diffusion MRI, Tractography, Tractometry, Anomaly Detection, Autoencoder

##  1 Introduction

Considerable effort has gone into designing methods for group comparisons in
dMRI (Jones and Cercignani, 2010) (i.e., N patients vs M controls) and as
such, single-subject analysis frameworks (i.e., 1 patient vs M controls) are
currently lacking. For clinically-heterogeneous groups, normative modeling
(Marquand et al., 2016) has been proposed, but often relies on voxel-based
methods, and hence is suboptimal for dMRI since WM tracts offer a more
suitable manifold. In this work, we investigate individual differences in WM
microstructure in children with copy number variants (CNVs) at high genetic
risk of neurodevelopmental and psychiatric disorders (Chawner et al., 2019),
which are relatively rare and challenging to recruit for research imaging
studies (Villalón-Reina et al., 2019). We propose the following unsupervised
framework for anomaly detection: First, we learn a normative set of features
derived from microstructural tract profiles obtained from typically developing
(TD) children. Second, we apply the framework to unseen subjects, to determine
whether these deviate from controls (based on the hypothesis that deviations
will stand out from the normative distribution).

##  2 Methods

###  2.1 Data acquisition & preprocessing

Diffusion data from 90 TD (age 8-18 years) and 8 children with a CNV and no
apparent WM lesions (age 8-15 years) were acquired on a Siemens 3T Connectom
MRI scanner with 14 b0 images, 30 directions at b = 500, 1200 s/mm2, 60
directions at b = 2400, 4000, 6000 s/mm2 and 2$\times$2$\times$2mm3 voxels.
Data were preprocessed as in Chamberland et al. (2019) and rotationally-
invariant spherical harmonic (RISH0, Mirzaalian et al. (2015)) features were
derived for each subject using the b = 6000 s/mm2 shell. Automated WM tract
segmentation was performed using TractSeg (Wasserthal et al., 2018) to obtain
20 bundles of interest (Fig. LABEL:fig:pipeline, left). For each bundle,
Tractometry (Bells et al., 2011; Cousineau et al., 2017) was performed
(sampling at 20 locations, Chamberland et al. (2019)) and the resulting 20
tract profiles were concatenated to form a feature vector (n = 20 tracts
$\times$ 20 locations = 400 features) for each subject.

fig:pipeline

Figure 1: Left: RISH0 feature mapped over 20 WM bundles with Tractometry.
Right: The proposed unsupervised anomaly detection framework trained on
healthy data only.

###  2.2 Anomaly detection

Our autoencoder architecture consists of five fully connected symmetric
layers. The input and output layers have the same number of nodes as the
length of the tract features vector. The inner layers consecutively apply a
compression ratio of 2 by reducing the number of nodes by half, up to the
bottleneck hidden layer (activation: ReLU, epochs: 25, batch size: 24,
learning rate: 1.0e-3, optimiser: Adam, loss: mean squared error).

A validation set (n = 16) was generated and held-out by combining the
individuals with a CNV (n = 8) with a random subset of TD (n = 8). The rest of
the TD (n = 82) data was used to establish a normative distribution (Fig.
LABEL:fig:pipeline, right). 10% of the TD data was held out for testing during
the training phase where the goal was to generate an output ($\hat{x}$)
similar to the input ($x$) by minimising the reconstruction error and
computing the mean absolute error over all features as anomaly score. Age
regression and feature normalization were performed on the training set and
subsequently applied to the validation set. We repeated this entire process 50
times to assess variations within the TD population and derived a mean anomaly
score for each subject. Using the subject labels, we report the mean ROC area
under the curve (AUC) across the iterations. In addition, we compared our
results with two previously-reported approaches: 1) univariate z-score
distribution (Yeatman et al., 2012); and 2) multivariate PCA combined with the
Mahalanobis distance (Yeatman et al., 2018; Sarica et al., 2017; Taylor et
al., 2020) using the aforementioned bootstrapped approach.

##  3 Results

The autoencoder approach was better at identifying CNV subjects as outliers
(AUC: 0.86 $\pm$ 0.06) compared with z-score (AUC: 0.53 $\pm$ 0.06) and PCA
(AUC: 0.61 $\pm$ 0.09). Fig. LABEL:fig:error shows the reconstructed features
of a random CNV subject which highlights significant discrepancies along
various association pathways, suggesting subject-specific differences in
microstructural attributes.

fig:error

Figure 2: RISH0 features (20 sections $\times$ 20 tracts) of a CNV subject
before ($x$) and after ($\hat{x}$) reconstruction. Discrepancies were
statistically identified along major association pathways using bootstrapped
permutations (shaded area, p < 0.01). In comparison, features from a
representative TD show no significant anomalies (bottom right).

##  4 Discussion & future work

The framework enabled subject- and tract-specific characterisation of WM
microstructure. By training on only healthy data, our findings revealed that
clinical cases (CNVs) can be classified as outliers. Furthermore, inspection
of discrepancies in tract profiles allow the identification of clinically
relevant differences in microstructural attributes. The framework also
outperformed traditional multivariate outlier detection mostly due to its
ability to handle high-dimensional data non-linearly. This extends the
possibility of using anomaly detection in clinically-heterogeneous groups (Kia
et al., 2018; Wolfers et al., 2020) and extremely rare cases (as little as
N=1), where group comparisons are otherwise impossible. Future work will
include additional cohorts to better assess the generalizability of the
framework and its application to other pathology. Importantly, anomalies
should be linked to clinical findings (Taylor et al., 2020). Our Tractometry-
based anomaly detection framework paves the way to progress from the
traditional paradigm of group-based comparison of patients against controls,
to a personalised medicine approach, and takes us a step closer in
transitioning microstructural MRI from the bench to the beside. The browser-
based tool will be made freely-available to the community via Github.

The authors thank Prof. Maxime Descoteaux and Jean-Christophe Houde
(Sherbrooke Connectivity Imaging Lab) for their helpful discussions and code
sharing. This work was supported by a Wellcome Trust Investigator Award
(096646/Z/11/Z), a Wellcome Trust Strategic Award (104943/Z/14/Z), an EPSRC
equipment grant (EP/M029778/1) and a Sir Henry Wellcome Fellowship
(215944/Z/19/Z).

## References

  * Bells et al. (2011) Sonya Bells, Mara Cercignani, Sean Deoni, Yaniv Assaf, Ofer Pasternak, C John Evans, Alexander Leemans, and Derek K Jones.  Tractometry–comprehensive multi-modal quantitative assessment of white matter along specific tracts.  In _Proc. ISMRM_ , volume 678, 2011.
  * Chamberland et al. (2019) Maxime Chamberland, Erika P. Raven, Sila Genc, Kate Duffy, Maxime Descoteaux, Greg D. Parker, Chantal M.W. Tax, and Derek K. Jones.  Dimensionality reduction of diffusion MRI measures for improved tractometry of the human brain.  _NeuroImage_ , 200:89–100, oct 2019.  10.1016/j.neuroimage.2019.06.020.
  * Chawner et al. (2019) Samuel JRA Chawner, Michael J Owen, Peter Holmans, F Lucy Raymond, David Skuse, Jeremy Hall, and Marianne BM van den Bree.  Genotype–phenotype associations in children with copy number variants associated with high neuropsychiatric risk in the uk (imagine-id): a case-control cohort study.  _The Lancet Psychiatry_ , 6(6):493–505, 2019\.
  * Cousineau et al. (2017) Martin Cousineau, Pierre-Marc Jodoin, Eleftherios Garyfallidis, Marc-Alexandre Cote, Felix C Morency, Verena Rozanski, Marilyn GrandMaison, Barry J Bedell, and Maxime Descoteaux.  A test-retest study on Parkinson’s PPMI dataset yields statistically significant white matter fascicles.  _NeuroImage: Clinical_ , 16:222–233, 2017.
  * Jones and Cercignani (2010) Derek K. Jones and Mara Cercignani.  Twenty-five pitfalls in the analysis of diffusion MRI data.  _NMR in Biomedicine_ , 23(7):803–820, sep 2010\.
  * Kia et al. (2018) Seyed Mostafa Kia, Christian F Beckmann, and Andre F Marquand.  Scalable multi-task gaussian process tensor regression for normative modeling of structured variation in neuroimaging data.  _arXiv preprint arXiv:1808.00036_ , 2018.
  * Marquand et al. (2016) Andre F. Marquand, Iead Rezek, Jan Buitelaar, and Christian F. Beckmann.  Understanding Heterogeneity in Clinical Cohorts Using Normative Models: Beyond Case-Control Studies.  _Biological Psychiatry_ , 80(7):552–561, oct 2016\.
  * Mirzaalian et al. (2015) Hengameh Mirzaalian, Amicie de Pierrefeu, Peter Savadjiev, Ofer Pasternak, Sylvain Bouix, Marek Kubicki, Carl-Fredrik Westin, Martha E Shenton, and Yogesh Rathi.  Harmonizing diffusion MRI data across multiple sites and scanners.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 12–19. Springer, 2015.
  * Sarica et al. (2017) Alessia Sarica, Antonio Cerasa, Paola Valentino, Jason Yeatman, Maria Trotta, Stefania Barone, Alfredo Granata, Rita Nisticò, Paolo Perrotta, Franco Pucci, and Aldo Quattrone.  The corticospinal tract profile in amyotrophic lateral sclerosis.  _Human Brain Mapping_ , 38(2):727–739, feb 2017\.
  * Taylor et al. (2020) Peter Neal Taylor, Nádia Moreira da Silva, Andrew Blamire, Yujiang Wang, and Rob Forsyth.  Early deviation from normal structural connectivity: A novel intrinsic severity score for mild tbi.  _Neurology_ , 2020.
  * Villalón-Reina et al. (2019) Julio E Villalón-Reina, Kenia Martínez, Xiaoping Qu, Christopher RK Ching, Talia M Nir, Deydeep Kothapalli, Conor Corbin, Daqiang Sun, Amy Lin, Jennifer K Forsyth, et al.  Altered white matter microstructure in 22q11. 2 deletion syndrome: a multisite diffusion tensor imaging study.  _Molecular psychiatry_ , pages 1–14, 2019.
  * Wasserthal et al. (2018) Jakob Wasserthal, Peter Neher, and Klaus H. Maier-Hein.  TractSeg - Fast and accurate white matter tract segmentation.  _NeuroImage_ , 183:239–253, dec 2018.
  * Wolfers et al. (2020) Thomas Wolfers, Christian F Beckmann, Martine Hoogman, Jan K Buitelaar, Barbara Franke, and Andre F Marquand.  Individual differences v. the average patient: mapping the heterogeneity in adhd using normative models.  _Psychological Medicine_ , 50(2):314–323, 2020\.
  * Yeatman et al. (2012) Jason D Yeatman, Robert F Dougherty, Nathaniel J Myall, Brian A Wandell, and Heidi M Feldman.  Tract profiles of white matter properties: automating fiber-tract quantification.  _PloS one_ , 7(11):e49790, 2012.
  * Yeatman et al. (2018) Jason D. Yeatman, Adam Richie-Halford, Josh K. Smith, Anisha Keshavan, and Ariel Rokem.  A browser-based tool for visualization and analysis of diffusion MRI data.  _Nature Communications_ , 9(1), dec 2018.