# Computation of Regions of Attraction for Hybrid Limit Cycles Using
Reachability: An Application to Walking Robots

Jason J. Choi, Ayush Agrawal, Koushil Sreenath, Claire J. Tomlin, and Somil
Bansal This research is supported in part by NSF Grants CMMI-1931853,
CMMI-1944722, an ONR Basic Research Challenge Program on Multibody Systems,
NASA ULI program (through Stanford) on Safe Aviation Autonomy, and the DARPA
Assured Autonomy Program. The work of Jason Choi received the support of a
fellowship from Kwanjeong Educational Foundation. Jason Choi, Ayush Agrawal,
Koushil Sreenath, and Claire Tomlin are with the University of California,
Berkeley, and Somil Bansal is with the University of Southern California, Los
Angeles. Contact info: {jason.choi, ayush.agrawal, koushils,
<EMAIL_ADDRESS><EMAIL_ADDRESS>Project website:
https://hybridreachability.github.io/website-hybrid-roa/

###### Abstract

Contact-rich robotic systems, such as legged robots and manipulators, are
often represented as hybrid systems. However, the stability analysis and
region-of-attraction computation for these systems are often challenging
because of the discontinuous state changes upon contact (also referred to as
state resets). In this work, we cast the computation of region-of-attraction
as a Hamilton-Jacobi (HJ) reachability problem. This enables us to leverage HJ
reachability tools that are compatible with general nonlinear system dynamics,
and can formally deal with state and input constraints as well as bounded
disturbances. Our main contribution is the generalization of HJ reachability
framework to account for the discontinuous state changes originating from
state resets, which has remained a challenge until now. We apply our approach
for computing region-of-attractions for several underactuated walking robots
and demonstrate that the proposed approach can (a) recover a bigger region-of-
attraction than state-of-the-art approaches, (b) handle state resets,
nonlinear dynamics, external disturbances, and input constraints, and (c) also
provides a stabilizing controller for the system that can leverage the state
resets for enhancing system stability.

##  I Introduction

Practical legged locomotion is one of the fundamental problems in robotics.
Legged robots can exhibit efficient, versatile, and highly dynamic gaits that
allow them to operate in highly unstructured environments. Mathematically,
legged robots are often represented as hybrid dynamical systems. Hybrid
dynamical systems may exhibit both continuous and discrete dynamics, wherein
swinging of a leg has continuous dynamics, whereas contacts with ground are
well-modeled as discrete events [1]. Under hybrid models, walking gaits are
represented as periodic hybrid limit cycles, i.e., a limit cycle involving
discontinuous jumps in the system state when the leg impacts the ground.

To obtain a stable walking behavior, it is often desirable to compute the
region-of-attraction (RoA) for a walking gait – the robot configurations from
which it can eventually converge to the corresponding hybrid limit cycle and
follow the gait. Other than analyzing the stability of a gait, RoA estimation
can also be useful for designing different candidate control laws for walking
and switching between different gaits [2]. Despite significant progress over
the last few years, stability analysis and RoA computation for legged robots
is still quite challenging because the dynamics are highly nonlinear and
intrinsically hybrid, requiring controlling through contact.

In this work, we cast the RoA computation as a Hamilton-Jacobi (HJ)
reachability problem. In reachability analysis, one computes the Backward
Reachable Tube (BRT) of a dynamical system – the set of states such that the
trajectories that start from this set will eventually reach some given target
set despite the worst case disturbance [3]. If we use the limit cycle as the
target set, the corresponding BRT is the RoA of the limit cycle. The advantage
of casting RoA computation as a reachability problem is the availability of
numerical tools that can compute BRT for a variety of nonlinear systems in the
presence of disturbances and control bounds. Even though computationally
intensive compared to other existing methods that can solve for the
underapproximations of RoAs, the reachability-based methods recover the
maximal RoA without imposing any restriction on its shape.

However, one problem that still remains is that, so far, we lack a general
framework for the reachability analysis of dynamical systems that involve
discontinuous state changes (also called state resets). To overcome this
challenge, we extend the reachability framework to dynamics with discontinuous
state resets. Our key insight is that since HJ reachability analysis
ultimately corresponds to solving a dynamic programming problem, we can reason
about the effect of state resets on the BRT by applying the Bellman principle
of optimality at the switching surface. We demonstrate that this corresponds
to a value remapping between post and pre reset states during the BRT
computation. Moreover, the existing numerical algorithms to compute BRTs can
easily be extended to account for this remapping without incurring any
additional computation cost. Thus, the obtained value function, the
corresponding BRT, as well as the stabilizing controller implicitly reason
about the effect of state resets on the stability of the system.

We apply the proposed approach on various “dynamic walking” systems and
demonstrate that it can obtain a bigger RoA than most state-of-the-art
approaches, while accounting for input constraints, state resets, and external
disturbances at the same time. Finally, the proposed approach also provides a
stabilizing controller for the system that can leverage state resets for a
faster stabilization.

##  II Related Work

_Poincaré Map-based methods._ Stability analysis of limit cycles can be
reduced to the problem of analyzing stability of the fixed point of the
corresponding Poincaré map. In particular, if the fixed point of the Poincaré
map is stable, then there exists a stable limit cycle for the original system
that passes through this fixed point [4, 5]. However, finding an analytical
expression for the Poincaré map is challenging; hence, a common approach is to
compute a numerical approximation of its linearization [1, 6]. The local
stability and RoA of the limit cycle are then computed using this discrete-
time linear system [2]. This approach has also been used to analyze stability
in the presence of external disturbances [7].

_Lyapunov-based methods._ One of the limitations of Poincaré analysis is that
the verifiable RoA is restricted to the set of states on the chosen Poincaré
section. The works in [8, 9, 10] use linearized transverse dynamics to
estimate the RoA by constructing a (time-varying) Lyapunov function via Sum-
of-Squares (SOS) programs. As noted in [8], a limitation of this approach is
that the size of the computed RoA depends on the choice of the transverse
coordinates. The work in [11] also proposes a Lyapunov based approach to
compute RoA for nonlinear hybrid systems with state resets using the
_complementarity_ formulation of contact.

In several of the aforementioned methods, Lyapunov functions are computed
using SOS programming. The key advantage of SOS-based methods is that they can
leverage convex optimization to compute fast underapproximations of the RoA,
leading to significant computation gains over reachability-based methods.
However, as a tradeoff, the verified RoA can be a significant
underapproximation of the true RoA. First, the SOS methods restrict the
Lyapunov functions to the class of (piecewise) SOS polynomials. Second,
although one can typically design a convex SOS program when the feedback
controller is fixed, solving for the controller and the Lyapunov function
together to maximize the verifiable RoA results in a nonconvex program which
induces another source of underapproximation when we obtain a local optimum.
These limitations are exacerbated further in the presence of state resets as
the Lyapunov condition at the switching surface leads to a nonconvex program
in general [8]. Ultimately, this results in recovering a smaller RoA. In
contrast, the proposed approach does not impose any specific structure on the
value function and hence can recover the maximal RoA (subjected to numerical
errors), even in the presence of state resets.

_Stabilizing controllers for hybrid systems with state resets._ The _Hybrid
Zero Dynamics_ (HZD) framework [12, 13, 1] was one of the earliest works on
creating stabilizable hybrid limit cycles and stabilizing controllers for
walking robots with discrete impacts. The method of Poincaré can also be
incorporated into the HZD to develop _event-based_ controllers to stabilize or
switch between hybrid limit cycles [14, 15, 16]. Building on the HZD
framework, [17] proposes a _rapidly exponentially stabilizing_ control
Lyapunov function based quadratic program (CLF-QP) to exponentially stabilize
hybrid limit cycles. There is also a rich literature on contact implicit
trajectory optimization [18, 19] and model predictive control [20, 21] that
can simultaneously handle multiple contact phases. Some of these methods have
also leveraged value remapping to account for the state resets. In this work,
we extend this principle to reachability analysis and focus on the RoA
computation, as opposed to trajectory optimization.

##  III Problem Setup

We consider the following class of hybrid systems:

$\displaystyle\dot{x}$ $\displaystyle=f(x,u,d),$ $\displaystyle x\notin S$ (1)
$\displaystyle x^{+}$ $\displaystyle=\Delta\left(x^{-}\right),$ $\displaystyle
x^{-}\in S,$

where $x\in\mathbb{R}^{n}$ is the state, $u\in\mathcal{U}$ is the control
input, and $d\in\mathcal{D}$ is the disturbance. $S$ indicates the switching
surface where the reset map $\Delta$ is applied; under the reset event at time
$t$, the state $x^{-}:=\lim_{\tau\nearrow t}x(\tau)$ instantaneously shifts to
$x^{+}$. Therefore, the resulting trajectory of this hybrid system is composed
of a set of continuous trajectories along the vector flow $f$, and
discontinuous jumps between them whenever the state hits the switching surface
$S$.

Hybrid limit cycle: Let $x^{*}(\cdot)$ represents a non-trivial $T$-periodic
solution of the system that undergoes $N$ impacts at times
$\\{t_{1},t_{2},\ldots,t_{N}\\}+kT$ for non-negative integer $k$. We call
$x^{*}(\cdot)$ a hybrid limit cycle of the system.

We assume that this limit cycle is given along with a baseline control law
$\pi_{0}(\cdot):\mathbb{R}^{n}\rightarrow\mathcal{U}$ that makes the limit
cycle forward invariant and stable in the small neighborhood around it. There
are many existing works that study approaches for designing such stable hybrid
limit cycles [13, 12] as well as the locally-stabilizing control laws [17, 1].

The overall objective is twofold. First, we seek to compute a region of state
space $\Omega(x^{*})\subseteq\mathbb{R}^{n}$ around $x^{*}(\cdot)$, from which
there exists an admissible control that can asymptotically stabilize the
system to $x^{*}(\cdot)$. We call $\Omega(x^{*})$ the (asymptotically)
stabilizable region or RoA for $x^{*}(\cdot)$. Second, we also want to design
a corresponding control law that stabilizes the system trajectory to
$x^{*}(\cdot)$.

###### Remark 1

Note that depending on the baseline controller, the RoA under $\pi_{0}(\cdot)$
can be very small. Our goal in this work is to find a stabilizing control law
that is beyond the capability of the baseline controller and provides us with
the maximal RoA $\Omega(x^{*})$.

Running example (Teleporting Dubins Car).

Figure 1: Teleporting Dubins car system. The car “teleports” to the other side
of the x-axis whenever it hits the x-axis from a positive y-direction.

We now introduce a simple hybrid system example that we will use throughout
the paper to illustrate our method. We introduce state resets to a well-known
Dubins car system which makes the car “teleport” (Fig. 1). It is chosen
because the key principles of our method are well exposed in its design and
results; more practical examples will be introduced in Sec VI. The dynamics
are given by

$\dot{p}_{x}=v\cos\theta,\quad\dot{p}_{y}=v\sin\theta,\quad\dot{\theta}=u,$
(2)

where $x=[p_{x},p_{y},\theta]$ is the state of the system consisting of car’s
position $(p_{x},p_{y})$ and its heading $\theta$. $v$ is the constant
velocity of the car, and $u\in[-\bar{\omega},\bar{\omega}]$ is the control
input, which is the angular speed. Next, we introduce a switching surface and
the reset map as the following:

$\displaystyle S$ $\displaystyle=\\{x\;|\;p_{y}=0,\sin\theta<0\\}$ (3)
$\displaystyle\quad x^{+}$
$\displaystyle=\Delta(x^{-})=\left[-R\cdot(|p_{x}|/R)^{\alpha}\cdot\text{sign}(p_{x}),~{}~{}p_{y},~{}~{}\pi+\theta\right]$

Intuitively, the car undergoes a discrete jump whenever it hits the $x$-axis
from a positive $y$-direction. Upon the reset, the car’s $x$-position and
heading change discontinuously. The parameter $\alpha>0$ in the reset map
dictates the nonlinearity of the reset map. For instance, $\alpha=1$ makes
$p_{x}^{+}=-p_{x}$, i.e., the car position is mirrored around the y-axis. For
$\alpha<1$, the reset map has a contraction property–everytime the state
resets, the value of $|p_{x}|$ gets closer to $R$. On the other hand, value of
$\alpha>1$ induces an expanding reset map.

The desired hybrid limit cycle is a semi-circular trajectory of radius $R$
centered at the origin in the counterclockwise direction. Therefore, the orbit
can be expressed as

$x^{*}(\cdot)=\left\\{x~{}\middle|~{}\|(p_{x},p_{y})\|\\!=\\!R,\;\theta\\!=\\!\frac{\pi}{2}+\tan^{-1}(\frac{p_{y}}{p_{x}}),\frac{\pi}{2}\leq\theta<\frac{3\pi}{2}\right\\}.$

We are interested in computing the RoA for $x^{*}(\cdot)$, i.e., the states
from which the car can be stabilized to follow the semi-circular trajectory
eventually, as well as the corresponding stabilizing controller. A baseline
stabilizing controller $\pi_{0}(\cdot)$ is designed by feedback linearization,
with the output defined as $y=\|(p_{x},p_{y})\|-R$.

##  IV Background

###  IV-A Hamilton-Jacobi Reachability

We present an overview of Hamilton-Jacobi (HJ) reachability analysis, first
for continuous systems without the disturbance: the dynamics of the system is
given by $\dot{x}=f(x,u)$. Let $\xi_{x,t}^{u}(\tau)$ denote the state achieved
at time $\tau$ by starting at initial state $x$ and initial time $t$, and
applying input functions $u(\cdot)$ over $[t,\tau]$. Given a target set
$\mathcal{L}\subset\mathbb{R}^{n}$, we are interested in computing the
Backward Reachable Tube (BRT) of $\mathcal{L}$:

Backward Reachable Tube (BRT) is the set of initial states from which the
agent acting optimally will eventually reach the target set $\mathcal{L}$
within the time horizon $[t,T]$:

$\small{\mathcal{V}(t)=\\{x\;|\;\exists
u(\cdot):[t,T]\rightarrow\mathcal{U},\exists\tau\in[t,T],\xi_{x,t}^{u}(\tau)\in\mathcal{L}\\}}.\vspace{-0.25em}$
(4)

In HJ reachability, the computation of BRT is formulated as an optimal control
problem, which can be solved by using the principle of dynamic programming.

First, a Lipschitz continuous target function $l(x)$ is defined whose zero
sub-level set is the target set $\mathcal{L}$, i.e. $\mathcal{L}=\\{x:l(x)\leq
0\\}$. Typically, $l(x)$ is defined as a signed distance function to
$\mathcal{L}$. The BRT seeks to find all states that could enter $\mathcal{L}$
at any point in the time horizon. This is computed by finding the minimum
distance to $\mathcal{L}$ over time:

$J(x,t,u(\cdot))=\min_{\tau\in[t,T]}l(\xi_{x,t}^{u}(\tau)).\vspace{-0.25em}$
(5)

Our goal is to capture this minimum distance for optimal trajectories of the
system. Thus, we compute the optimal control that minimizes this distance (and
eventually drives the system into the target set when $l(\cdot)\leq 0$). The
value function corresponding to this optimal control problem is:

$\small{V(x,t)=\inf_{u(\cdot)}\Big{\\{}J\Big{(}x,t,u(\cdot)\Big{)}\Big{\\}}.}\vspace{-0.25em}$
(6)

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

$\displaystyle\min\Big{\\{}D_{t}V(x,t)+H(x,t),l(x)-V(x,t)\Big{\\}}=0,$ (7)

with the terminal value function $V(x,T)=l(x)$. $D_{t}$ and $\nabla$ represent
the time and spatial gradients of the value function. $H$ is the Hamiltonian
that encodes the role of dynamics and the optimal control inputs:

$\displaystyle H(x,t)=\min_{u}$ $\displaystyle\langle\nabla
V(x,t),f(x,u)\rangle.$ (8)

Intuitively, the term $l(x)-V(x,t)$ in (7) enforces the value function to
“memorize” the best record (closest instance to the target set) of the optimal
trajectories. Once $V(x,t)$ is obtained, the BRT is given as the zero sub-
level set of the value function $\mathcal{V}(t)=\\{x\;|\;V(x,t)\leq 0\\}$. The
corresponding optimal control for reaching the target set $\mathcal{L}$ is
derived as

$\displaystyle u^{*}(x,t)=\arg\min_{u}$ $\displaystyle\langle\nabla
V(x,t),f(x,u)\rangle.$ (9)

Finally, note that every state that is contained in a finite-time BRT of
$\mathcal{L}$, represents an initial condition from which the system can be
controlled to $\mathcal{L}$ within finite time.

###  IV-B Casting RoA computation as a HJ reachability problem

We now explain how HJ reachability can be used to estimate the RoA for a
hybrid limit cycle, $\Omega(x^{*})$. Let $h(x)$ denote distance to the limit
cycle $x^{*}(\cdot)$, i.e., $h(x):=\inf_{y\in x^{*}(\cdot)}\|x-y\|$. We define
the target function as $l(x):=h(x)-\epsilon$ for some $\epsilon>0$. Thus, the
target set $\mathcal{L}$ is an $\epsilon$-tube around the orbit
$x^{*}(\cdot)$. Note that for a stable hybrid limit cycle, we can pick
$\epsilon$ such the baseline controller $\pi_{0}$ can stabilize to the limit
cycle from anywhere inside $\mathcal{L}$.

Given the value function, $V(x,t)$ and the BRT, $\mathcal{V}(t)$,
corresponding to this target set, we can stabilize to the limit cycle from any
state within the BRT using the following control law:

$u(x,t)=\begin{cases}\arg\min_{u}\langle\nabla V(x,t),f(x,u)\rangle,&\text{if
}x\notin\mathcal{L}\\\ \pi_{0}(x),&\text{otherwise}\end{cases}$ (10)

Intuitively, we can use the controller provided by HJ reachability to steer
the system to $\mathcal{L}$ from any state within the BRT; upon reaching
$\mathcal{L}$, we can switch to the baseline controller. Therefore, for every
state that can be verified as an element of a finite-time BRT of
$\mathcal{L}$, we can conclude that it is an element of $\Omega(x^{*})$. In
other words, $\mathcal{V}(t)$ can be used as an estimation of $\Omega(x^{*})$.
The corresponding stabilizing controller is given by (10).

###### Remark 2

Note that if we pick $\epsilon=0$, only the states that can achieve finite-
time convergence to $x^{*}(\cdot)$ can be verified from HJ reachability
analysis. However, this would exclude many states that are not finite-time
stabilizable but still asymptotically stabilizable to $x^{*}(\cdot)$. In fact,
even for a state where stronger notion like exponential stabilizability [17]
is achievable, its time-to-reach to the orbit is not necessarily a finite
value. Therefore, by using $\epsilon>0$ together with the baseline policy
$\pi_{0}$ which makes sure that every state in $\mathcal{L}$ is stabilizable
to the actual orbit, we can capture the states that are not necessarily
finite-time but are still asymptotically stabilizable to $x^{*}(\cdot)$ and
are inside the finite-time BRT of $\mathcal{L}$. For more detailed analysis
and a constructive way to choose $\epsilon$, please refer to [22]. Finally,
defining $\mathcal{L}$ as the $\epsilon$-neighborhood of $x^{*}(\cdot)$ is
also necessary for numerical methods to compute $V(x,t)$ and to tackle the
presence of disturbance, which will be explained in the next subsection.

Running example (Teleporting Dubins Car). The target function $l(x)$ for the
running example is $l(x):=\|(p_{x},p_{y})\|-R-\epsilon$. The Hamiltonian can
be computed analytically in this case (see [23]):

$H(x,t)=\frac{\partial V}{\partial p_{x}}v\cos\theta+\frac{\partial
V}{\partial p_{y}}v\sin\theta-\bar{\omega}|\frac{\partial
V}{\partial\theta}|.$ (11)

###  IV-C Robustness to bounded disturbance

For systems with bounded disturbance $d\in\mathcal{D}$, HJ reachability can
again be used to synthesize optimal controllers that are robust to the
disturbance. The definition of BRT is now the set of initial states for which,
under worst-case disturbances, the agent acting optimally will eventually
reach the target set $\mathcal{L}$ within time $[t,T]$:

$\small{\mathcal{V}(t)=\\{x:\exists u(\cdot),\forall
d(\cdot):[t,T]\rightarrow\mathcal{D},\exists\tau\in[t,T],\xi_{x,t}^{u,d}(\tau)\in\mathcal{L}\\}}.\vspace{-0.25em}$
(12)

Computation of the BRT can be formulated as a zero-sum game between the
control and the disturbance, and can be solved in a similar way by applying
the dynamic programming [23]. $V(x,t)$ is governed by the same variational
inequality as (7), where the Hamiltonian is as follows:

$\displaystyle H(x,t)=\max_{u}\min_{d}$ $\displaystyle\langle\nabla
V(x,t),f(x,u,d)\rangle.$ (13)

Finally, the optimal control signal that is robust to the disturbance can be
synthesized by the following rule:

$\displaystyle u^{*}(x,t)=\arg\max_{u}\min_{d}$ $\displaystyle\langle\nabla
V(x,t),f(x,u,d)\rangle.$ (14)

For the specific problem of stabilization to the hybrid limit cycle, it must
be noted that any non-zero disturbance signal can break the forward-invariance
of the baseline controller $\pi_{0}(\cdot)$ on the limit cycle $x^{*}(\cdot)$.
Therefore, the notion of $x^{*}(\cdot)$ being robustly asymptotically stable
is defined as an existence of an $\epsilon$-neighborhood of the limit cycle to
which the system can be asymptotically stabilized under the worst-case
disturbance [24]. Under a mild assumption that the baseline controller
$\pi_{0}(\cdot)$ itself renders the set $\mathcal{L}$ locally robustly stable,
any elements of the finite-time BRT for $\mathcal{L}$ can be verified as
robustly asymptotically stabilizable to $x^{*}(\cdot)$.

##  V Hamilton-Jacobi Reachability for Dynamics With State Resets

We now present the main contribution of this paper, which is an extension of
the reachability framework above (for continuous dynamics) to account for
discontinuous state resets. First, the following theorem proves that the value
function in the presence of state resets can be obtained by solving a
constrained version of HJI-VI in (7).

###### Theorem 1

$V(x,t)$ in the presence of state resets is given by the solution to the
following constrained VI:

$\displaystyle\min\Big{\\{}D_{t}V(x,t)+H(x,t),l(x)-V(x,t)\Big{\\}}=0,$
$\displaystyle x\notin S,$ (15) $\displaystyle
V(x,t)=V(\Delta\left(x\right),t),$ $\displaystyle x\in S,$

with the terminal value function

$\displaystyle V(x,T)$ $\displaystyle=l(x),$ $\displaystyle x\notin S,$ (16)
$\displaystyle V(x,T)$ $\displaystyle=l(\Delta\left(x\right)),$ $\displaystyle
x\in S.$

Intuitively, Theorem 1 states that the value function can be obtained by
solving the usual HJI-VI for the states that are not on the switching surface.
For the states that are on the switching surface, the value is given by that
of the corresponding post reset state. This is not surprising because if the
state is on the switching surface, it will instantaneously change to the post
reset state. Consequently, how “close” the state is to the target set is
completely characterized by how close the post reset state is to the target
set. The proof of Theorem 1 formalizes this intuition using the Bellman
principle of optimality and is presented in Appendix A in the supplementary
material.

The interesting aspect about the HJI-VI in (15) is that since it reasons about
the state resets, the obtained value function and the stabilizing controller
implicitly exploit the reset map to reach the target set as quickly as
possible. We will demonstrate this aspect further in Sec. VI-A.

Numerical Implementation. We now present a numerical method to compute the
value function based on the result in Theorem 1. The baseline of the numerical
algorithm is a standard level set-based method that computes the value
function over a discretized computational domain of a compact subset
$\mathcal{X}\subset\mathbb{R}^{n}$ [25]. Let $\\{p(x)\\}_{G}$ denote the set
of values of a function $p(x)$ evaluated over a state space grid $G$. Since
our computation will proceed backward in time, we will numerically solve the
variational inequality in (15) for each discrete time step in a finite
interval $[0,T]$ using the three-step update rule described in Algorithm 1.

Input : Target function $\\{l(x)\\}_{G}$, Time horizon $T$

Output : Value function $\\{V(x,t)\\}_{G}$ for $t=0,\delta,2\delta,\cdots,T$

Initialization: $\\{V(x,T)\\}_{G\setminus S}=\\{l(x)\\}_{G\setminus S}$,

$\\{V(x,T)\\}_{S}=l(x^{+})$,

$t=T$;

while  _$t >0$_ do

$\text{For all }x\in G\setminus S,$ $\displaystyle V(x,t-\delta)=$
$\displaystyle V(x,t)$ (17)
$\displaystyle+\min_{u\in\mathcal{U}}\max_{d\in\mathcal{D}}D_{x}V(x,t)\cdot
f(x,u)\;\delta;$ $V(x,t-\delta)=\min\\{V(x,t-\delta),l(x)\\}$;

$\text{For all }x\in S,$ $\displaystyle V(x,t-\delta)=V(\Delta(x),t-\delta);$
(18)

$t=t-\delta$;

end while

Algorithm 1 Value Function Computation for Dynamics with State Resets.

The numerical algorithm is implemented by means of the computational methods
provided in [25]. It is important to stress the remarkable computational
similarity of this new algorithm to its non-reset counterpart in (7). Indeed,
the only computational overhead is introduced by the step in (18). As a
result, as will be demonstrated in the following section, our method can
compute the BRT at essentially no additional cost compared to the case where
there are no state resets.

###### Remark 3

Note that Theorem 1 and the numerical method in Algorithm 1 are not specific
to the RoA computation problem, which is the main focus of this paper. In
general, these can be applied to any kind of target set $\mathcal{L}$ that can
be expressed as a zero-sublevel set of a Lipschitz continuous function $l(x)$
and has a non-empty volume.

##  VI Case Studies

In this section, we will demonstrate how our extended HJ reachability
framework can be used to compute RoA for various hybrid systems, including a
simple biped walker.

###  VI-A Running example: Teleporting Dubins car

We first apply our method to the running example. We use the following
parameters for our simulations:
$v\\!=\\!$1\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$,
$R\\!=\\!$2\text{\,}\mathrm{m}$$,
$\bar{\omega}\\!=\\!$1\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$$,
$\alpha\\!=\\!0.5$, which induces a contracting reset map, and $\epsilon=0.2$.

To confirm that our new method is able to capture states that exploit discrete
jumps to reach the limit cycle, we compare our method to a simple baseline –
we compute the set of states that are stabilizable to the limit cycle only by
following the continuous part of the dynamics. Thus, for this baseline, we
“freeze” the dynamics on the switching surface, i.e. set $f(x,u)\equiv 0$
whenever $x\in S$.

Fig. 2 visualizes the BRTs computed for $T\\!=\\!$6.3\text{\,}\mathrm{s}$$
from our method (blue), and under the frozen dynamics (brown). It shows the 2D
slices of the BRTs at a fixed heading angle
$\theta=-$3\pi\text{/}4\text{\,}\mathrm{rad}$$. Note that the BRT computed by
freezing the dynamics is strictly contained in the BRT computed by our method.
The gap between the two sets uncovers the states that can reach the target set
within $6.3\text{\,}\mathrm{s}$ only via leveraging discrete jumps. Fig. 3
displays the resulting optimal trajectory (blue) from a particular initial
state that is contained in this gap (marked with the red cross). It exploits
two state jumps to quickly reach the target set after
$4.3\text{\,}\mathrm{s}$; at every discrete jump, the distance to the orbit
reduces significantly. In contrast, the feedback linearization-based baseline-
controller (grey), although doing its best to get closer to the orbit by
heading towards the radial direction, is able to reach the target set only
after $11.9\text{\,}\mathrm{s}$. Furthermore, the discrete jump near the limit
cycle is not “intended” by the baseline controller, since the control law is
not aware of the reset map. In contrast, the proposed method is explicitly
able to reason that since the reset map is a contraction in this case, it can
be exploited to reduce the distance to the orbit and stabilize to the limit
cycle faster.

Figure 2: $p_{x}$-$p_{y}$ slices of the Backward Reachable Tubes for the
teleporting Dubins car at $\theta=-3\pi/4$, computed by value remapping
(blue), and by freezing the dynamics (brown) on the switching surface
($T\\!=\\!$6.3\text{\,}\mathrm{s}$$, $\alpha=0.5$). The gap between the two
sets captures the states that can reach the target set within
$6.3\text{\,}\mathrm{s}$ only by using the “teleports” (reset map). We
evaluate an initial state included in the BRT (cross) in Fig.3.  Figure 3: The
optimal trajectory simulated for $6.3\text{\,}\mathrm{s}$ which reaches the
target set after $4.3\text{\,}\mathrm{s}$ by exploiting two discrete
transitions. For comparison, the feedback-linearization-based baseline
controller is simulated for $14\text{\,}\mathrm{s}$, which reaches the target
set after $11.9\text{\,}\mathrm{s}$.  Figure 4: $p_{x}$-$p_{y}$ slices of the
Backward Reachable Tubes at $\theta=-3\pi/4$ (top) and $\theta=\pi/4$
(bottom), with various values of the reset map nonlinearity parameter $\alpha$
under the value remapping, and under dynamics freezing, applied to the
switching surface ($T\\!=\\!$6.3\text{\,}\mathrm{s}$$). Note that the reset
map is contracting when $\alpha<1$, and expanding when $\alpha>1$.

In Fig. 2, we can observe that many states adjacent to the switching surface
are captured as the stabilizable region to the limit cycle. This is because
$\alpha=0.5$ induces the contracting property of the reset map, which is
beneficial for the stabilization. Fig. 4 visualizes BRTs computed for the same
$T$ under three different values of $\alpha$. Note that $\alpha=3.0$ is
unfavorable for the stabilization as it is a divergent reset map; as a result,
most of the region adjacent to the switching surface is not captured in the
BRT. In fact, with even larger $\alpha$, the BRT eventually shrinks to the BRT
computed by freezing the dynamics on the switching surface, i.e., the reset
map is so divergent that it is no longer possible to stabilize the system
after hitting the switching surface.

Computation details: The grid used for the computation is constructed in the
polar coordinate system of the dynamics ($r=\|(p_{x},p_{y})\|$,
$\alpha=\tan^{-1}(p_{y}/p_{x})$, $\theta$); it is a
$81\\!\times\\!41\\!\times\\!81$ grid over the state space of
$[0,16]\\!\times\\![0,\pi]\\!\times\\![-\pi,\pi]$. It takes under an hour to
compute the BRT of $T\\!=\\!$6.3\text{\,}\mathrm{s}$$ on a computer with a
2.6GHz 6‑core 9th‑generation Intel Core i7 processor CPU.

###  VI-B Rimless wheel system

Figure 5: (Left) Rimless wheel system. (Right) Compass-gait walker.

We next apply the proposed approach to the rimless wheel system. The rimless
wheel is a popular model of a passive-dynamic walker. One of the reasons for
its popularity is that RoA can be computed analytically for the rimless wheel
[26], and hence it can be used to analyze the efficacy of a RoA computation
algorithm.

The system consists of a central mass with several spokes extending radially
outward (see Fig. 5); it has no control input. At any given time, one of the
spokes is pinned at the ground, and the system follows the dynamics of a
simple pendulum, $f(\theta,\dot{\theta})=[\dot{\theta},\sin(\theta)]^{T},$
where $x:=(\theta,\dot{\theta})$ is the state of the system. When an unpinned
spoke contacts the ground, an inelastic collision takes place and there is an
instantaneous change in $\dot{\theta}$. Upon collision, the pinned spoke
leaves the ground and the new spoke becomes the pinned one. The switching
surface and the reset map are given as:

$S=\\{x\;|\theta=\alpha+\gamma\\},\quad
x^{+}=\Delta(x^{-})=\left[\begin{array}[]{c}2\gamma-\theta\\\
\cos(2\alpha)\dot{\theta}\end{array}\right].$

Figure 6: Regions of Attraction (RoA) for the rimless wheel limit cycle
(black). (Left) The true RoA is shown in light pink. RoA obtained using SOS
programming is shown in purple. The figure is reproduced from [8]. (Right) RoA
computed using the proposed approach. The proposes approach is able to recover
the entire RoA.

For a sufficiently inclined slope, the rimless wheel exhibits several stable
hybrid limit cycles that are given as

$x^{*}(\cdot)=\\{x\;|\cos(\theta)+\frac{1}{2}\dot{\theta}^{2}=E\\},$

where $E$ is the total initial energy in the system. One such limit cycle is
shown in black in Fig. 6 (left). The figure is reproduced from [8]. The right
boundary of the plot is the switching surface–whenever the system state
reaches the right boundary, it resets back to the left boundary as per the
reset map. The light pink region is the exact RoA of the system [26]. The
purple region is the RoA obtained using a time-varying Lyapunov function,
which is synthesized using SOS programming [8]. As evident from the figure,
using the SOS programming, we are able to recover only a subset of the RoA. In
Fig. 6 (right), we demonstrate the RoA obtained using our method (dark pink).
The proposed approach is able to recover the entire RoA of the system.

Computation details: We use a grid of size $201\times 201$ over the state
space $[-0.2,0.6]\times[-0.6,1.3]$. It takes within a few minutes to compute
the BRT on a standard laptop.

###  VI-C Compass-gait walker

We next consider a compass-gait walker, which consists of two links with an
actuated joint between them. Similar to the rimless wheel example, we consider
a _pinned_ model of the robot, with the configuration variable
$q\coloneqq[q_{1},q_{2}]^{T}$ as illustrated in Fig. 5, We consider a hybrid
model of walking, with alternating phases of a continuous-time single-support
phase followed by an instantaneous inelastic impact of the swing leg with the
ground. The dynamics in single-support can be represented by the standard
manipulator equation, and the reset map $\Delta$ is derived from an inelastic
impact model [1]. By defining the state as $x\coloneqq[q,\dot{q}]^{T}$, the
continuous dynamics can be expressed as control-affine: $\dot{x}=f(x)+g(x)u$.
The switching surface is defined as the set of states where the swing foot
strikes the ground with a negative velocity and the stance leg angle crosses a
predefined threshold $\bar{q}_{1}$,

$\displaystyle
S:=\\{x~{}|~{}q_{1}\leq\bar{q}_{1},2q_{1}+q_{2}=0,2\dot{q}_{1}+\dot{q}_{2}<0\\}.$
(19)

The target hybrid limit cycle we wish to stabilize to is designed as a
reference walking gait [12, 13, 27], which is tracked by a baseline controller
designed based on feedback linearization and CLF-QP [17]. The configuration of
the swing leg on the gait, $q_{2,r}$, can be parametrized as a polynomial
function of the stance leg angle $q_{1}$, and we define the output as
$y:=q_{2}-q_{2,r}$ for feedback linearization. Each walking step takes
$0.73\text{\,}\mathrm{s}$, where $q_{1}$ starts with $q_{1,d}=0.13$ and ends
at $-q_{1,d}$. The phase plot of the reference gait in $q_{1}$-$q_{2}$ space
is visualized in Fig. 7(a) black dashed line.

Figure 7: Backward Reachable Tubes (BRTs) for the compass-gait walker
visualized by their projection on the $q_{1}$–$q_{2}$ space.
($T\\!=\\!$1.5\text{\,}\mathrm{s}$$ for all three cases) (a) BRT for the
closed-loop dynamics under the CLF-QP. (b) BRT with the control bound
$u\in[-4,4]$, without the disturbance. (c) BRT with the control bound
$u\in[-4,4]$, and the disturbance bound $d\in[-0.3,0.75]$. The color indicates
the thickness of the BRT in the 4D space; it gets thicker as
yellow$\rightarrow$blue.

In Fig. 7, the BRTs computed from our method for
$T\\!=\\!$1.5\text{\,}\mathrm{s}$$ is visualized as the projection of the 4D
set to the $q_{1}$-$q_{2}$ space (colored region). Note that the chosen time
horizon encompasses more than two nominal walking steps. First, we apply our
method to the closed loop dynamics of the baseline controller (Fig. 7(a)).
Although it is designed to be locally stabilizing, the verified RoA comprises
only a small neighborhood of the reference gait. By contrast, Fig. 7(b) shows
the verified stabilizable region for the gait by allowing full control
authority defined by the control bound $u\in[-4,4]$\text{\,}\mathrm{N}$$. The
bound is selected such that on the limit cycle, the CLF-QP barely saturates
this bound, so that we do not give additional advantage to the optimal
controller. However, by exploiting full control capacity and the reset map, we
are able to verify that the actual stabilizable region is much bigger than
what can be achieved by the CLF-QP controller. Fig. 8 shows the
$q_{2}$–$\dot{q_{2}}$ slices of the BRT along the reference gait. Note that
the gait itself is actually contained in the BRT, and that the BRT gets
smaller as it gets closer to the end of the step (Fig.8(c)), because of the
imminent reset.

Figure 8: $q_{2}$–$\dot{q_{2}}$ slices of the BRT corresponding to Figure
7-(a), along the hybrid periodic orbit (reference gait). (a) Initial phase:
$[q_{1},\dot{q_{1}}]=[0.13,-0.67]$. (b) Intermediate phase:
$[q_{1},\dot{q_{1}}]=[0,-0.36]$. (c) Terminal phase (imminent of reset):
$[q_{1},\dot{q_{1}}]=[-0.13,-0.91]$. The black dots indicate the corresponding
$[q_{2},\dot{q_{2}}]$ values on the orbit, which are always contained in the
BRT. The color indicates the value of $V(x,T)$; it is yellow when $V(x,T)=0$
and gets blue as the value get smaller. The diagrams on the top right corners
visualize the swing leg (blue) configuration range that is contained in the
BRT, together with the stance leg (magenta) for each phase.

In Fig. 9, we show some trajectories that start at initial states which are
contained in the verified BRT. These initial states have a large deviation
from the reference gait, corresponding to a large perturbation. We start by
applying the reachability-based optimal control (blue), and then switch to the
CLF-QP controller (orange) when we reach the target set. In Fig. 9(b,c), the
optimal controller is exploiting the reset map to converge to the limit cycle
faster. Physically, this means that the robot steps on the ground first, and
use this ground force to converge to the nominal gait at the next walking
step. For comparison, the trajectories from the same initial states where the
CLF-QP is used from the very beginning are visualized together (grey), which
all fail to recover from the ground impact and cannot stabilize to the limit
cycle. In particular, in Fig. 9(c), the CLF-QP fails to overcome a large
momentum, resulting in the robot to fall backwards.

Figure 9: Phase plots of the reachability-based optimal trajectories
(blue$\rightarrow$orange) and the CLF-QP trajectories (grey) in
$q_{1}$–$q_{2}$ space for various initial states from which the CLF-QP fails
to stabilize to the gait. The trajectories are evaluated for 4 walking steps.
Each start of the step is indicated by the index of the step in italic type.
The evolution of time is indicated by the color fading out. The optimal
trajectory switches from HJ optimal control to CLF-QP once it reaches the
target set (indicated by the transition from blue to orange curve).

Finally, we introduce disturbance to the robot’s dynamics: an unmodeled
repulsive or stiction torque applied to the joint between the leg. This can be
expressed as $\dot{x}=f(x)+g(x)u+\dot{q}_{2}g(x)d$. There could be various
source of uncertainty for this disturbance, for instance, unmodeled motor
dynamics or joint friction. The bound of the disturbance is $d\in[-0.75,0.3]$;
under this bound, on the limit cycle, the CLF-QP is still able to robustly
stabilize the gait. The BRT computed under the robust optimal control setting
(as explained in Sec. IV-C) is visualized in Fig. 7.(c). Since now we only
verify states that are robustly stabilizable to the gait, the resulting BRT is
smaller than the non-robust BRT in Fig. 7.(b). However, it still shows that by
applying the robust optimal controller, it is able to robustly stabilize a
larger set of states than the CLF-QP. For instance, the CLF-QP fails to
stabilize a state slightly perturbed from the gait under the presence of
disturbance $d\\!=\\!-0.75$ (Fig.10.(b)), because the CLF-QP cannot reason
about the additional torque required to compensate the effect of the
disturbance. In contrast, the robust stabilizing controller obtained using the
proposed approach is simultaneously able to counter the disturbances and
stabilize to the reference gait (Fig.10.(c)).

Figure 10: Phase plots of the trajectories in $q_{1}$–$q_{2}$ space evaluated
for 4 walking steps for a same initial state: (a) CLF-QP under no model-plant
mismatch. (b) CLF-QP under model-plant mismatch, which fails to stabilize to
the orbit. (c) Optimal control, which is able to reach the orbit under the
same model-plant mismatch. Each start of the step is indicated by the index of
the step in italic type. The evolution of time is indicated by the color
fading out.

Computation details: We use a grid of size $41\times 81\times 81\times 81$
over the state space of
$[-0.52,0.52]\times[-1.04,1.04]\times[-4,4]\times[-8,8]$, taking roughly 36
hours to compute the BRT of $1.5\text{\,}\mathrm{s}$ for this demonstration
using a 3.3GHz 12-core Intel Core i7-5820K CPU.

##  VII Conclusion and Future Work

We present a reachability-based approach to compute regions-of-attraction for
dynamical systems with state resets. We apply our approach to several
underactuated walking robots and demonstrate that the proposed approach can
synthesize a bigger RoA compared to the state-of-the-art approaches and also
provides a stabilizing controller. However, the proposed approach has several
limitations that we will like to address in future. First, the reachability
computation scales exponentially with the number of states, limiting its
direct use to relatively low-dimensional systems. Recently, there have been
some promising results that use Neural PDE solvers for computing reachable set
for high-dimensional systems [28]. We will explore using these methods to
extend the proposed approach to more complex walking robots. For providing
stability guarantees, it would be interesting to leverage recent sampling-
based methods [29] that provide probabilistic guarantees on the learned value
functions. Second, it will be interesting to extend the proposed approach to
other types of reachability problems, such as reaching a target set while
avoiding some undesirable states, e.g., obstacles. Applying the proposed
approach to real-world robotic systems will be another exciting future
direction.

## References

  * [1] Eric R Westervelt et al.  “Feedback control of dynamic bipedal robot locomotion”  CRC press, 2018
  * [2] Mohamad Shafiee Motahar, Sushant Veer and Ioannis Poulakakis  “Composing limit cycles for motion planning of 3D bipedal walkers”  In _CDC_ , 2016
  * [3] S. Bansal, M. Chen, S. Herbert and C.. Tomlin  “Hamilton-Jacobi Reachability: A Brief Overview and Recent Advances”  In _CDC_ , 2017
  * [4] John Guckenheimer and Philip Holmes  “Structurally stable heteroclinic cycles”  In _Mathematical Proceedings of the Cambridge Philosophical Society_ , 1988  Cambridge University Press
  * [5] JW Grizzle  “Remarks on event-based stabilization of periodic orbits in systems with impulse effects”  In _Int. Symposium on Communications, Control and Signal Processing_ , 2006
  * [6] Benjamin Morris and Jessy W Grizzle  “A restricted Poincaré map for determining exponentially stable periodic orbits in systems with impulse effects: Application to bipedal robots”  In _CDC_ , 2005
  * [7] Sushant Veer and Ioannis Poulakakis  “Input-to-state stability of periodic orbits of systems with impulse effects via Poincaré analysis”  In _TAC_ , 2019
  * [8] Ian R Manchester, Mark M Tobenkin, Michael Levashov and Russ Tedrake  “Regions of attraction for hybrid limit cycles of walking robots”  In _IFAC Proceedings Volumes_ 44.1 Elsevier, 2011, pp. 5801–5806
  * [9] Ian R Manchester  “Transverse dynamics and regions of stability for nonlinear hybrid limit cycles”  In _IFAC Proceedings Volumes_ 44.1 Elsevier, 2011, pp. 6285–6290
  * [10] Leonid B Freidovich, Anton S Shiriaev and Ian R Manchester  “Stability analysis and control design for an underactuated walking robot via computation of a transverse linearization”  In _IFAC Proceedings Volumes_ 41.2 Elsevier, 2008, pp. 10166–10171
  * [11] Michael Posa, Mark Tobenkin and Russ Tedrake  “Lyapunov analysis of rigid body systems with impacts and friction via sums-of-squares”  In _HSCC_ , 2013
  * [12] Jessy W Grizzle, Gabriel Abba and Franck Plestan  “Asymptotically stable walking for biped robots: Analysis via systems with impulse effects”  In _TAC_ , 2001
  * [13] Eric R Westervelt, Jessy W Grizzle and Daniel E Koditschek  “Hybrid zero dynamics of planar biped walkers”  In _TAC_ , 2003
  * [14] Kaveh Akbari Hamed, Brian G Buss and Jessy W Grizzle  “Exponentially stabilizing continuous-time controllers for periodic orbits of hybrid systems: Application to bipedal locomotion with ground height variations”  In _IJRR_ , 2016
  * [15] Ching-Long Shih, JW Grizzle and Christine Chevallereau  “From stable walking to steering of a 3D bipedal robot with passive point feet”  In _Robotica_ 30.7 Cambridge University Press, 2012, pp. 1119–1130
  * [16] A. Agrawal and K. Sreenath  “Discrete Control Barrier Functions for Safety-Critical Control of Discrete Systems with Application to Bipedal Robot Navigation.”  In _RSS_ , 2017
  * [17] Aaron D Ames, Kevin Galloway, Koushil Sreenath and Jessy W Grizzle  “Rapidly exponentially stabilizing control lyapunov functions and hybrid zero dynamics”  In _TAC_ , 2014
  * [18] Michael Posa, Cecilia Cantu and Russ Tedrake  “A direct method for trajectory optimization of rigid bodies through contact”  In _IJRR_ , 2014
  * [19] He Li and Patrick M Wensing  “Hybrid systems differential dynamic programming for whole-body motion planning of legged robots”  In _RA-L_ , 2020
  * [20] S. Cleac’h, T. Howell, M. Schwager and Z. Manchester  “Linear Contact-Implicit Model-Predictive Control”  In _arXiv_ , 2021
  * [21] Tobia Marcucci et al.  “Approximate hybrid model predictive control for multi-contact push recovery in complex environments”  In _Humanoids_ , 2017
  * [22] Fabio Camilli, Lars Grüne and Fabian Wirth  “Control Lyapunov Functions and Zubov’s Method”  In _SIAM Journal on Control and Optimization_ , 2008
  * [23] I.M. Mitchell, A.M. Bayen and C.J. Tomlin  “A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games”  In _TAC_ , 2005
  * [24] Randy A Freeman and Petar V Kokotovic  “Inverse optimality in robust stabilization”  In _SIAM journal on control and optimization_ , 1996
  * [25] I. Mitchell  “A toolbox of level set methods”  In _http://www. cs. ubc. ca/mitchell/ToolboxLS/toolboxLS.pdf_ , 2004
  * [26] Michael Jon Coleman  “A stability study of a three-dimensional passive-dynamic model of human gait”  Cornell University, 1998
  * [27] Ayonga Hereid and Aaron D Ames  “Frost: Fast robot optimization and simulation toolkit”  In _IROS_ , 2017
  * [28] Somil Bansal and Claire Tomlin  “DeepReach: A deep learning approach to high-dimensional reachability”  In _ICRA_ , 2021
  * [29] N. Boffi et al.  “Learning stability certificates from data”  In _CoRL_ , 2020

## Appendix A

###  VII-A Proof of Theorem 1

We will first prove the result for a state that is not on the switching
surface. Let $x(t)\notin S$. The result for this case follows from the proof
of standard HJI-VI. Here we present an informal proof based on the Taylor
series expansion of the value function; a more rigorous proof of HJI-VI can be
found in [23].

Consider a small $\delta>0$ such that there is no state reset in the time
interval $[t,t+\delta]$. The dynamic programming principle for the cost
function in (5) implies that

$\displaystyle V(x(t),t)$
$\displaystyle=\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}\min\\{l(x(t)),V(x(t+\delta),t+\delta)\\},$
$\displaystyle=\min\\{l(x(t)),\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}V(x(t+\delta),t+\delta))\\},$
(20)

On the other hand, the Taylor expansion of $V(x(t+\delta),t+\delta)$ can be
written as:

$V(x(t+\delta),t+\delta)=V(x(t),t)+D_{t}V(x(t),t)\delta+\nabla V(x(t),t)\delta
x+\text{h.o.t},$

where $\delta x$ is change in the state and can be approximated as
$f(x,u,d)\delta$. Ignoring the higher order terms and plugging the Taylor
approximation in (VII-A), we have

$\displaystyle V(x(t),t)\approx\min\\{l(x(t)),$ $\displaystyle
V(x(t),t)+D_{t}V(x(t),t)\delta$
$\displaystyle+\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}\nabla
V(x(t),t)\cdot f(x,u,d)\delta\\},$

where we have separated terms that do not depend on $u$ and $d$. Subtracting
$V(x(t),t)$ on both sides, we get

$\displaystyle\min\\{$ $\displaystyle l(x(t))-V(x(t),t),$
$\displaystyle\delta\left[D_{t}V(x(t),t)+\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}\nabla
V(x(t),t)\cdot f(x,u,d)\right]\\}\approx 0,$

Since the above equation holds for all $\delta>0$, we must have

$\displaystyle\min\\{$ $\displaystyle l(x(t))-V(x(t),t),$
$\displaystyle\left[D_{t}V(x(t),t)+\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}\nabla
V(x(t),t)\cdot f(x,u,d)\right]\\}\approx 0,$

which results into the desired HJI VI.

We now turn our attention to a state $x(t)\in S$. Once again, the dynamic
programming principle for the cost function in (5) implies that

$\displaystyle V(x(t),t)$
$\displaystyle=\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}\min\\{l(x(t)),V(x(t+\delta),t+\delta)\\},$

Since $x(t)\in S$, the system state will instantaneously change to $x^{+}(t)$
under any control and disturbance action. Thus, we have

$\displaystyle V(x(t),t)$
$\displaystyle=\sup_{u\in\mathcal{U}}\inf_{d\in\mathcal{D}}\min\\{l(x^{+}(t)),V(x^{+}(t+\delta),t+\delta)\\},$
$\displaystyle=V(x^{+}(t),t),$

where the second equality follows from the definition of value function for
$x^{+}(t)$.


# Graph-based mutually exciting point processes for modelling event times in
docked bike-sharing systems

Francesco Sanna Passino  Department of Mathematics, Imperial College London
180 Queen’s Gate, SW7 2AZ, London  Yining Che  Department of Mathematics,
Imperial College London
180 Queen’s Gate, SW7 2AZ, London  Carlos Cardoso Correia Perello  Department
of Mathematics, Imperial College London
180 Queen’s Gate, SW7 2AZ, London

###### Abstract

This paper introduces graph-based mutually exciting processes (GB-MEP) to
model event times in network point processes, focusing on an application to
docked bike-sharing systems. GB-MEP incorporates known relationships between
nodes in a graph within the intensity function of a node-based multivariate
Hawkes process. This approach reduces the number of parameters to a quantity
proportional to the number of nodes in the network, resulting in significant
advantages for computational scalability when compared to traditional methods.
The model is applied on event data observed on the Santander Cycles network in
central London, demonstrating that exploiting network-wide information related
to geographical location of the stations is beneficial to improve the
performance of node-based models for applications in bike-sharing systems. The
proposed GB-MEP framework is more generally applicable to any network point
process where a distance function between nodes is available, demonstrating
wider applicability.

††footnotetext:  For the purpose of open access, the authors have applied a
Creative Commons Attribution (CC-BY) licence to any Author Accepted Manuscript
version arising.
CCCP completed this work at Imperial College London, and he is currently in
the Department of Mathematics at the University of Oxford.
Corresponding author: Francesco Sanna Passino –<EMAIL_ADDRESS>
Keywords — bike-sharing systems, dynamic networks, Hawkes processes, mutually
exciting processes.

##  1 Introduction

In recent years, a number of convenient and sustainable alternatives for
short-distance travel have emerged in world cities, with the objective of
reducing traffic congestion, and promoting healthier and more sustainable
urban environments. Bike-sharing systems (BSSs) represent one such solution,
which have become increasingly popular especially after the end of
restrictions to movement due to the COVID-19 pandemic (Heydari et al., 2021).
Two types of BSSs are usually offered in cities: traditional docked systems,
which require users to pick up and return bikes at designated stations, and
dockless systems, which utilise GPS technology, allowing users to locate and
unlock bikes using smartphone applications. From a statistical perspective,
the literature about BSSs is mostly focused on two problems: rebalancing,
consisting in the redistribution of bikes across stations and locations
according to the predicted demand (see, for example, the survey of Vallez et
al., 2021), and demand forecasting (for example, Rennie et al., 2023). In
general, modelling docked bike-sharing systems is still an active research
problem, extensively in the literature (see, for example, Choi et al., 2022;
Deng et al., 2023; Busatto and Stingo, 2023; Paul et al., 2023).

Motivated by the application to docked bike-sharing systems, this work
proposes a statistical network-informed model for event times, generically
called graph-based mutually exciting process (GB-MEP), which could be used to
predict demand across stations. The proposed methodology is tested on publicly
available data from the Santanter Cycles bike-sharing network in central
London. In this work, journeys observed on bike-sharing systems are
interpreted as a realisation of a network point process, with events having a
corresponding duration. In general, demand for bicycles in docked BSSs has
often been modelled as a multivariate point process in the literature (see,
for example, Gervini and Khanal, 2018; Okawa et al., 2019; Foschi, 2023), but
network information, expressing the interaction between different stations, is
often not explicitly incorporated within the existing models. This work
proposes a parsimonious representation for the temporal point processes
observed at each docking station, directly utilising the information provided
by the distance between stations in the network for imposing restrictions on
the model parameters.

In particular, this work considers a particular class of point processes,
called mutually exciting processes, with a focus on Hawkes processes (see, for
example, Hawkes, 1971; Laub et al., 2022). Hawkes processes have been
extensively used in the literature for modelling a number of different
phenomena, such as earthquakes (Ogata, 1988), crime (Mohler et al., 2011),
social media (Zhao et al., 2015; Rizoiu et al., 2017), emails (Fox et al.,
2016; Miscouridou et al., 2018; Sanna Passino and Heard, 2023), and computer
networks (Price-Williams and Heard, 2020; Shlomovich et al., 2022; Sanna
Passino and Heard, 2023). Hawkes processes appear to be a particularly
suitable model for a bike-sharing system, since they capture periods of
increased demand in a data-driven way, avoiding specific parametric
assumptions on the process, such as periodic components or weather conditions.
This property has been proved particularly powerful in cyber-security, where
self-exciting models largely outperform alternatives that model periodic
components parametrically (Price-Williams and Heard, 2020).

The remainder of this article is structured as follows: Section 2 gives
background on point processes on graphs and Hawkes processes, describing
related literature and challenges of existing approaches when applied to
network point processes. Next, Section 3 introduces graph-based mutually
exciting processes, the main contribution of this work, discussing the method,
parameter estimation, and model evaluation procedures. The proposed
methodology is then applied on publicly available data from the Santander
Cycles bike-sharing system in central London (Section 4), followed by a
conclusion and discussion in Section 5.

##  2 Background

Consider a network $\mathbb{G}=(\mathcal{V},\mathcal{E})$ where
$\mathcal{V}=\\{1,\dots,M\\},\ M\in\mathbb{N}$, is the node set, and
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the edge set, such that
$(x,y)\in\mathcal{E}$ if nodes $x$ and $y$ are connected. Within the context
of bike-sharing systems, nodes represent docking stations, and an edge between
two stations can be drawn if at least one trip between the stations is
observed within an observation period. It is further assumed that event data
with a duration are observed on the network, corresponding to a sequence of
quadruples $(x_{k},y_{k},t_{k},t_{k}^{\prime}),\ k=1,2,\dots,N,\
N\in\mathbb{N}$, where $(x_{k},y_{k})\in\mathcal{E}$ is an edge,
$t_{k}\in\mathbb{R}_{+}$ is the start time of the event, and
$t_{k}^{\prime}>t_{k}$ is the end time of the event. This corresponds to
observing events with duration $\delta_{k}=t_{k}^{\prime}-t_{k}$. In practice,
events are ordered by end time, such that
$t_{1}^{\prime}<\cdots<t_{N}^{\prime}$. In bike-sharing systems, the quadruple
$(x_{k},y_{k},t_{k},t_{k}^{\prime})$ denotes a bike journey between stations
$x_{k}$ and $y_{k}$, started at time $t_{k}$ and ended at time
$t_{k}^{\prime}$. From the observed sequence of quadruples, two left-
continuous counting processes $N_{i}:\mathbb{R}_{+}\to\mathbb{N}_{0}$ and
$N_{i}^{\prime}:\mathbb{R}_{+}\to\mathbb{N}_{0}$ could be associated with each
node:

$\displaystyle N_{i}(t)=\sum_{k\geq
1}\mathds{1}_{[0,t)}{(t_{k})}\mathds{1}_{\\{i\\}}{(x_{k})},$ $\displaystyle
N^{\prime}_{i}(t)=\sum_{k\geq
1}\mathds{1}_{[0,t)}{(t_{k}^{\prime})}\mathds{1}_{\\{i\\}}{(y_{k})},$ (2.1)

where $\mathds{1}_{\cdot}(\cdot)$ denotes the indicator function:
$\mathds{1}_{S}(x)=1$ if $x\in S$, and $0$ otherwise. In the BSS application,
the process $N_{i}(t)$ counts the number of journeys started from station $i$
before time $t$, whereas the process $N^{\prime}_{i}(t)$ counts the number of
journeys ended at station $i$ before time $t$. For simplicity, event times for
each node are often denoted with the notation $t_{i,k},\ i\in\mathcal{V},\
k\in\mathbb{N}$, representing the start time of the $k$-th event started from
node $i$, and $t^{\prime}_{i,k}$, corresponding to the ending time of the
$k$-th event ended on node $i$. This notation can be used to obtain a simpler
representation for $N_{i}(t)$ and $N_{i}^{\prime}(t)$:

$\displaystyle N_{i}(t)=\sum_{k\geq 1}\mathds{1}_{[0,t)}{(t_{i,k})},$
$\displaystyle N^{\prime}_{i}(t)=\sum_{k\geq
1}\mathds{1}_{[0,t)}{(t_{i,k}^{\prime})}.$ (2.2)

In general, a counting process could be characterised by its conditional
intensity function (CIF). The CIF $\lambda_{i}(t)$ for $N_{i}(t)$ is defined
as follows:

$\displaystyle\lambda_{i}(t)=\lim_{\mathrm{d}t\to
0}\frac{\mathbb{E}[N_{i}(t+\mathrm{d}t)-N_{i}(t)|\mathcal{H}(t)]}{\mathrm{d}t},$
(2.3)

where $\mathcal{H}(t)$ is the history of the process on the entire network up
to time $t$.

The counting processes $\\{N_{i}(t)\\}_{i\in\mathcal{V}}$ and
$\\{N_{i}^{\prime}(t)\\}_{i\in\mathcal{V}}$ could be considered as
realisations of multivariate temporal point processes. For this class of
processes, a popular model in the literature is the multivariate Hawkes
process (Hawkes, 1971) with scaled exponential kernel (see, for example, Laub
et al., 2022), which models the CIF $\lambda_{i}(t)$ on each node as a
function of all start times of events observed in the entire network before
time $t$:

$\lambda_{i}(t)=\lambda_{i}+\sum_{j=1}^{M}\sum_{t_{j,k}<t}\alpha_{ij}\exp\\{-\beta_{ij}(t-t_{j,k})\\},$
(2.4)

where $\lambda_{i}\in\mathbb{R}+$ is a baseline intensity,
$\alpha_{ij}\in\mathbb{R}_{+}$ is the jump in the CIF of node $i$ caused by an
event on node $j$, and $\beta_{ij}\in\mathbb{R}_{+}$, with
$\alpha_{ij}<\beta_{ij}$, is the corresponding decay rate. The end times of
events could also be considered in the intensity (2.4), obtaining the
following self-and-mutually exciting Hawkes process CIF (see, for example,
Laub et al., 2022):

$\lambda_{i}(t)=\lambda_{i}+\sum_{j=1}^{M}\left[\sum_{t_{j,k}<t}\alpha_{ij}\exp\\{-\beta_{ij}(t-t_{j,k})\\}+\sum_{t_{j,h^{\prime}}<t}\alpha_{ij}^{\prime}\exp\\{-\beta_{ij}^{\prime}(t-t_{j,h}^{\prime})\\}\right],$
(2.5)

where $\alpha_{ij}^{\prime},\beta_{ij}^{\prime}\in\mathbb{R}_{+}$,
$\alpha_{ij}^{\prime}<\beta_{ij}^{\prime}$. The models in Equations (2.4) and
(2.5) have a number of issues when applied to temporal point processes on
networks: first, the number of parameters is $\mathcal{O}(M^{2})$ across the
entire network, implying that it grows quadratically with the number of nodes,
which is unfeasible for estimation in practice, especially when the network
becomes large; second, the network structure is not explicitly incorporated
within the intensity function in Equation (2.4); third, information from
additional processes, such as end times of events, is not included within the
CIF.

In the application to bike-sharing systems, these issue appear to be
particularly relevant: the number of docking station could be very large
(approximately $800$ in the Santander Cycles BSS in central London), making
estimation of the parameters $\alpha_{ij}$ and $\beta_{ij}$ a high-dimensional
problem. Using the example of $M=800$, the total number of parameters required
for model in Equation (2.4) for the entire network would be
$M+2M^{2}=\numprint{1280800}$. This number is more than the average number of
bike journeys per month in the busiest year in the Santander Cycles BSS (over
$\numprint{950000}$ hires per month in 2022, see Transport for London, 2022).
Additionally, another drawback of the multivariate Hawkes process models in
Equation (2.4) and (2.5) is that information about known, pre-existing network
structure is not incorporated within the CIF: for modelling network point
processes in real-world applications, information about the nodes in the BSS
could be very important. For example, in a bike-sharing system, stations which
are geographically close are expected to have similar characteristics compared
to stations which are further apart, suggesting that geographical information
about nodes should be taken into account within the model.

It must be remarked that the model in Equation (2.4) has been often used
successfully in the literature for modelling multivariate point processes,
including spatio-temporal processes (see, for example, the survey of Reinhart,
2018). In particular, Hawkes processes are in general suitable to discover
existing network structure and causal relationships between nodes (see, for
example, Linderman and Adams, 2014; Etesami et al., 2016; Xu et al., 2016;
Eichler et al., 2017; Yuan et al., 2019), but known characteristics about the
network and nodes are not explicitly incorporated within the model. Using the
approach of Linderman and Adams (2014), the network adjacency matrix, if
known, could be used to set the coefficients $\alpha_{ij}$ to zero where
$(i,j)\not\in\mathcal{E}$, but in general the use of node covariates and known
network information has been overlooked in the literature. Also, it must be
remarked that Hawkes process models on networks are usually applied at two
different levels of resolution: node-based (see, for example, Fox et al.,
2016; Price-Williams and Heard, 2020), or edge-based (see, for example
Blundell et al., 2012; Perry and Wolfe, 2013; Miscouridou et al., 2018; Huang
et al., 2022; Sanna Passino and Heard, 2023).

This work mostly concerns with node-based models, proposing an extension of
the standard multivariate Hawkes process is proposed, exploiting existing
network structure between the nodes to reduce the number of parameters to
estimate to $\mathcal{O}(M)$, providing advantages for scalability to large
networks. Additionally, information from end times is added to the CIF,
providing a more complete modelling framework for event data on networks with
an associated event duration.

##  3 Graph-based mutually exciting processes

The main objective of this work is to propose a graph-based mutually exciting
model (GB-MEP) for the node-specific conditional intensity functions
$\lambda_{i}(t)$ for the temporal point processes of start times of connection
events on a graph, applied to journeys observed on bike-sharing systems. The
GB-MEP intensity function is assumed to take the following generic form:

$\lambda_{i}(t)=\lambda_{i}+\sum_{j=1}^{M}\left\\{\sum_{k=1}^{N_{j}(t)}\omega_{i}\\{\gamma(i,j),t-t_{j,k}\\}+\sum_{h=1}^{N_{j}^{\prime}(t)}\omega_{i}^{\prime}\\{\gamma(i,j),t-t^{\prime}_{j,h}\\}\right\\},$
(3.1)

where $\lambda_{i}\in\mathbb{R}_{+}$ is a baseline rate,
$\omega_{i},\omega_{i}^{\prime}:\mathbb{R}_{+}\times\mathbb{R}_{+}\to\mathbb{R}_{+}$
are spatio-temporal excitation kernel functions, and
$\gamma:\mathcal{V}\times\mathcal{V}\to\mathbb{R}_{+}$ is a distance function
between nodes. The distance function $\gamma$ quantifies how close in space
two nodes are. For example, in bike-sharing systems, the distance function
could be taken to be the haversine metric, which provides an “as-the-crow-
flies” measure of distance between the geographic location of two nodes. Given
the latitude and longitude $(\varphi_{x},\ell_{x})$ and
$(\varphi_{y},\ell_{y})$ of the bike stations $x$ and $y$ (in radians), the
haversine distance between the stations is:

$\gamma(x,y)=2\rho\arcsin\left\\{\sqrt{\sin^{2}\left(\frac{\varphi_{y}-\varphi_{x}}{2}\right)+\cos(\varphi_{x})\cos(\varphi_{y})\sin^{2}\left(\frac{\ell_{y}-\ell_{x}}{2}\right)}\right\\},$
(3.2)

where $\rho$ is the Earth radius. Note that, for more general point processes
on networks where information about the node locations is not present, the
distance $\gamma(x,y)$ between two nodes $x,y\in\mathcal{V}$ could depend on
graph-related statistics, such as the length of the shortest path between the
nodes, or the number of common neighbours. In this case $\gamma$ takes the
interpretation of a dissimilarity metric.

In this work, it is further assumed that the spatial and temporal components
of the excitation kernels $\omega_{i}$ are separable, with the temporal
element of the kernel taking a scaled exponential form. This results in the
following excitation kernel:

$\omega_{i}(\gamma^{\ast},t)=\kappa_{i}(\gamma^{\ast})\cdot\alpha_{i}\exp(-\beta_{i}t),$
(3.3)

where $\alpha_{i},\beta_{i}\in\mathbb{R}_{+}$, $\alpha_{i}<\beta_{i}$, and
$\gamma^{\ast}$ corresponds to a realised value of the distance function
$\gamma$. The function $\kappa_{i}:\mathbb{R}_{+}\to(0,1]$ rescales the node-
specific temporal mutual excitation given by $\alpha_{i}\exp(-\beta_{i}t)$,
using the distance $\gamma$ between the nodes. In the application to bike-
sharing systems, the excitation effect of an event occurring on a node with
small distance from station $i$ is likely to have a larger effect on the
intensity $\lambda_{i}(t)$ than an event occurring on a node located further
away. Therefore, $\kappa_{i}$ could be set to a monotonically decreasing
function such as $\kappa(0)=1$ and
$\lim_{\gamma^{\ast}\to\infty}\kappa_{i}(\gamma^{\ast})=0$. A possible choice,
used in this work, is the truncated exponential function
$\kappa_{i}(\gamma^{\ast})=\exp(-\theta_{i}\gamma^{\ast})\mathds{1}_{[0,\varepsilon)}(\gamma^{\ast})$,
where $\theta_{i}\in\mathbb{R}_{+}$ is a decay parameter and
$\varepsilon\in\mathbb{R}_{+}$ is a fixed constant defining the radius of
neighbouring nodes considered to have non-zero effect on the CIF
$\lambda_{i}(t)$. Note that setting $\varepsilon\to 0$ would result in
$\kappa_{i}(\gamma^{\ast})=\mathds{1}_{\\{0\\}}(\gamma^{\ast})$, which would
only give weight to events involving station $i$, removing the spatial mutual
excitation component, and making the process exclusively self-exciting. A
similar spatio-temporal excitation kernel could also be postulated for
$\omega^{\prime}_{i}(\gamma^{\ast},t)$:

$\omega_{i}^{\prime}(\gamma^{\ast},t)=\exp(-\theta_{i}^{\prime}\gamma^{\ast})\mathds{1}_{[0,\varepsilon)}\\{\gamma^{\ast}\\}\cdot\alpha_{i}^{\prime}\exp(-\beta_{i}^{\prime}t),$
(3.4)

with
$\alpha_{i}^{\prime},\beta_{i}^{\prime},\theta_{i}^{\prime}\in\mathbb{R}_{+}$
and $\alpha_{i}^{\prime}<\beta_{i}^{\prime}$. Overall, this gives the
following extended form for the proposed GB-MEP CIF in Equation (3.1) for
applications to BSSs:

$\lambda_{i}(t)=\lambda_{i}+\sum_{j=1}^{M}\Bigg{\\{}\sum_{k=1}^{N_{j}(t)}\exp\\{-\theta_{i}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}\exp\\{-\beta_{i}(t-t_{j,k})\\}\\\
+\sum_{h=1}^{N_{j}^{\prime}(t)}\exp\\{-\theta_{i}^{\prime}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha^{\prime}_{i}\exp\\{-\beta_{i}^{\prime}(t-t_{j,h}^{\prime})\\}\Bigg{\\}}.$
(3.5)

An example of the intensity functions arising under a GB-MEP process is given
in Figure 1.

Figure 1: Example of CIFs for a GB-MEP with $\mathcal{V}=\\{1,2,3\\}$ and
three realised events expressed as quadruples $(x,y,t,t^{\prime}),\
x,y\in\mathcal{V},\ t,t^{\prime}\in\mathbb{R}_{+},\ t<t^{\prime}$. Events and
their symbols: $(1,2,1.25,2.75)$ corresponds to $\bullet$; $(1,3,4,4.5)$ to
$\blacktriangle$, and $(2,2,2.35,8)$ to $\blacksquare$. Parameters:
$\lambda_{1}=0.2$, $\lambda_{2}=0.3$, $\lambda_{3}=0.15$, $\alpha_{1}=0.8$,
$\alpha_{2}=\alpha_{3}=0.6$, $\beta_{1}=\beta_{2}=\beta_{3}=1$,
$\theta_{1}=\theta_{2}=\theta_{3}=1$,
$\alpha_{1}^{\prime}=\alpha_{2}^{\prime}=0.5$, $\alpha_{3}^{\prime}=0.3$,
$\beta_{1}^{\prime}=\beta_{2}^{\prime}=\beta_{3}^{\prime}=1$,
$\theta_{1}^{\prime}=\theta_{2}^{\prime}=\theta_{3}^{\prime}=1.5$. Node
distances: $\gamma(1,2)=\gamma(2,1)=0.5$, $\gamma(1,3)=\gamma(3,1)=0.75$,
$\gamma(2,3)=\gamma(3,2)=0.25$, $\gamma(1,1)=\gamma(2,2)=\gamma(3,3)=0$.

###  3.1 Parameter estimation

The model parameters
$\bm{\phi}_{i}=(\lambda_{i},\alpha_{i},\beta_{i},\theta_{i},\alpha_{i}^{\prime},\beta_{i}^{\prime},\theta_{i}^{\prime})$
for each node-based process could be learned via maximum likelihood
estimation, by maximising the node-specific log-likelihood function separately
for each station. The log-likelihood function for a temporal point process
takes the following form (see, for example, Daley and Vere-Jones, 2002):

$\displaystyle\ell_{i}(\bm{\phi}_{i})=\sum_{k=1}^{N_{i}(T)}\log[\lambda_{i}(t_{i,k})]-\Lambda_{i}(T),$
(3.6)

where $\Lambda_{i}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ is the compensator,
defined as:

$\displaystyle\Lambda_{i}(t)=\int_{0}^{t}\lambda_{i}(s)\mathrm{d}s.$ (3.7)

For the CIF in Equation (3.5), the compensator takes the following form:

$\Lambda_{i}(t)=\lambda_{i}t\
-\sum_{j=1}^{M}\Bigg{\\{}\frac{\exp\\{-\theta_{i}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}}{\beta_{i}}\sum_{k=1}^{N_{j}(t)}[\exp\\{-\beta_{i}(t-t_{j,k})\\}-1]\\\
-\frac{\exp\\{-\theta_{i}^{\prime}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}^{\prime}}{\beta_{i}^{\prime}}\sum_{h=1}^{N_{i}^{\prime}(t)}[\exp\\{-\beta_{i}^{\prime}(t-t_{i,h}^{\prime})\\}-1]\Bigg{\\}}.$
(3.8)

Note that, because of the double summations in the CIF in Equation (3.5), the
computational complexity for calculating the term
$\sum_{k=1}^{N_{i}(T)}\log[\lambda_{i}(t_{i,k})]$ in Equation (3.6) is
$\mathcal{O}\\{N_{i}(T)\sum_{j=1}^{M}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}[N_{j}(T)+N_{j}^{\prime}(T)]\\}$,
proportional to the sum of the number of events occurring on each node,
multiplied by the total number of events started from the station under
consideration. In practice, this makes model fitting particularly challenging,
and almost unfeasible for large datasets. This issue is largely alleviated by
the scaled exponential excitation kernel, which admits a computationally
cheaper form of the log-likelihood, written in recursive form using the
following expression for the CIF in Equation (3.5):

$\lambda_{i}(t_{i,k})=\lambda_{i}+\sum_{j=1}^{M}\left[\exp\\{-\theta_{i}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}A_{ij}(k)+\exp\\{-\theta_{i}^{\prime}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha^{\prime}_{i}A_{ij}^{\prime}(k)\right].$
(3.9)

In the above expression, $A_{ij}(k)$ and $A_{ij}^{\prime}(k)$ are defined as
follows:

$\displaystyle
A_{ij}(k)=\sum_{l=1}^{N_{j}(t_{i,k})}\exp\\{-\beta_{i}(t_{i,k}-t_{j,l})\\},$
$\displaystyle
A_{ij}^{\prime}(k)=\sum_{l=1}^{N_{j}^{\prime}(t_{i,k})}\exp\\{-\beta_{i}^{\prime}(t_{i,k}-t_{j,l}^{\prime})\\}.$
(3.10)

The initial values $A_{ij}(1)$ and $A^{\prime}_{ij}(1)$ are simply:

$\displaystyle
A_{ij}(1)=\sum_{l=1}^{N_{j}(t_{i,1})}\exp\\{-\beta_{i}(t_{i,1}-t_{j,l})\\},$
$\displaystyle
A^{\prime}_{ij}(1)=\sum_{l=1}^{N_{j}^{\prime}(t_{i,1})}\exp\\{-\beta_{i}^{\prime}(t_{i,1}-t_{j,l}^{\prime})\\}.$
(3.11)

The recursive form of the log-likelihood stems from the following proposition.

######  Proposition 1.

The expressions for $A_{ij}(k)$ and $A_{ij}^{\prime}(k)$ in Equation (3.10)
can be rewritten recursively as follows:

$\displaystyle A_{ij}(k)=\exp\\{-\beta_{i}(t_{i,k}-t_{i,k-1})\\}\cdot
A_{ij}(k-1)+\sum_{l=N_{j}(t_{i,k-1})+1}^{N_{j}(t_{i,k})}\exp\\{-\beta_{i}(t_{i,k}-t_{j,l})\\},$
(3.12) $\displaystyle
A_{ij}^{\prime}(k)=\exp\\{-\beta_{i}^{\prime}(t_{i,k}-t_{i,k-1})\\}\cdot
A_{ij}^{\prime}(k-1)+\sum_{l=N_{j}^{\prime}(t_{i,k-1})+1}^{N_{j}^{\prime}(t_{i,k})}\exp\\{-\beta_{i}^{\prime}(t_{i,k}-t_{j,l}^{\prime})\\}.$
(3.13)

###### Proof.

The result is proved for $A_{ij}(k)$ here; the proof for $A_{ij}^{\prime}(k)$
is analogous. The summation for $A_{ij}(k)$ in Equation (3.10) could be broken
up at $N_{j}(t_{i,k-1})$ as follows:

$\displaystyle
A_{ij}(k)=\sum_{l=1}^{N_{j}(t_{i,k-1})}\exp\\{-\beta_{i}(t_{i,k}-t_{j,l})\\}+\sum_{l=N_{j}(t_{i,k-1})+1}^{N_{j}(t_{i,k})}\exp\\{-\beta_{i}(t_{i,k}-t_{j,l})\\}.$
(3.14)

Pre-multiplying the first element of Equation (3.14) by
$\exp\\{-\beta_{i}(t_{i,k-1})-\beta_{i}(t_{i,k-1})\\}$ gives:

$\displaystyle\sum_{l=1}^{N_{j}(t_{i,k-1})}\exp\\{-\beta_{i}(t_{i,k}-t_{j,l})\\}$
$\displaystyle=\exp\\{-\beta_{i}(t_{i,k-1})-\beta_{i}(t_{i,k-1})\\}\sum_{l=1}^{N_{j}(t_{i,k-1})}\exp\\{-\beta_{i}(t_{i,k}-t_{j,l})\\}$
(3.15)
$\displaystyle=\exp\\{-\beta_{i}(t_{i,k})-\beta_{i}(t_{i,k-1})\\}\sum_{l=1}^{N_{j}(t_{i,k-1})}\exp\\{-\beta_{i}(t_{i,k-1}-t_{j,l})\\}.$
(3.16)

By the definition in Equation (3.10),
$\sum_{l=1}^{N_{j}(t_{i,k-1})}\exp\\{-\beta_{i}(t_{i,k-1}-t_{j,l})\\}=A_{ij}(k-1)$,
which gives the result in Equation (3.12). ∎

It follows that, writing
$B_{ij}(k)=\exp\\{-\theta_{i}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}A_{ij}(k)+\exp\\{-\theta_{i}^{\prime}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}^{\prime}A_{ij}^{\prime}(k)$,
Proposition 1 gives a recursive form for the log-likelihood in Equation (3.6):

$\displaystyle\ell_{i}(\bm{\phi}_{i})=\sum_{k=1}^{N_{i}(T)}\log\left\\{\lambda_{i}+\sum_{j=1}^{M}B_{ij}(k)\right\\}-\Lambda_{i}(T).$
(3.17)

Under this recursive representation, $\ell_{i}(\bm{\phi}_{i})$ can be
evaluated in $\mathcal{O}\\{N_{i}(T)M_{i}^{\varepsilon}\\}$, where
$M_{i}^{\varepsilon}=\sum_{j=1}^{M}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}$
is the number of nodes within distance $\varepsilon$ of node $i$. Note that
this assumes that the update of $B_{ij}(k)$ given $B_{ij}(k-1)$ is
$\mathcal{O}(M_{i}^{\varepsilon})$.

###  3.2 Model evaluation

The quality of the model fit could be assessed by evaluating the distribution
of $p$-values calculated on a training and a test set observation periods,
using the same methodology discussed in Sanna Passino and Heard (2023). By the
time rescaling theorem (see, for example, Brown et al., 2002), the transformed
event times $\Lambda_{i}(t_{i,k}),\ k=1,2,\dots$ are event times of a
homogeneous Poisson process with unit rate, under the null hypothesis of
correct specification of the conditional intensity and its parameters.
Therefore, under the same null hypothesis, the upper tail $p$-values
calculated from the observed events should follow a uniform distribution in
$(0,1)$. Since the inter-arrival times of a unit rate Poisson process are
exponentially distributed with unit rate, the corresponding upper tail
$p$-values for the observed events are:

$p_{i,k}=\exp\left\\{-\int_{t_{i,k-1}}^{t_{i,k}}\lambda_{i}(s)\mathrm{d}s\right\\}=\exp\\{-\Lambda_{i}(t_{i,k})+\Lambda_{i}(t_{i,k-1})\\}.$
(3.18)

Using, the recursive expressions in Proposition 1 and the counting processes
in Equation (2.2), the expression in Equation (3.18) reduces to:

$p_{i,k}=\exp\Bigg{\\{}\sum_{j=1}^{M}\frac{\exp\\{-\theta_{i}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}}{\beta_{i}}[A_{ij}(k)-A_{ij}(k-1)-N_{j}(t_{i,k})+N_{j}(t_{i,k-1})]\\\
+\sum_{j=1}^{M}\frac{\exp\\{-\theta_{i}^{\prime}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}^{\prime}}{\beta_{i}^{\prime}}[A_{ij}^{\prime}(k)-A_{ij}^{\prime}(k-1)-N_{j}^{\prime}(t_{i,k})+N_{j}^{\prime}(t_{i,k-1})]-\lambda_{i}(t_{i,k}-t_{i,k-1})\Bigg{\\}}.$
(3.19)

Considering an observation period $[0,T]$, the model parameters could first be
estimated from a testing period $[0,T^{\ast}]$, $T^{\ast}<T$, using the
procedure described in Section 3.1. The $p$-values can then be calculated
directly by plugging in the estimates into Equation (3.18). The estimates
could also be used to evaluate the performance of the model in a testing
period $(T^{\ast},T]$, by calculating the $p$-values using the same
methodology. In both cases, under a correct model specification and parameter
estimation, the distributions of the $p$-value on the training and test set
events should be approximately uniform in $(0,1)$. The quality of the fit to
the uniform distribution could be evaluated using the Kolmorogov-Smirnov
statistic, denoted here $\mathrm{KS}_{i}$. For the sequence of node-specific
$p$-values $p_{i,k},\ k=1,\dots,N_{i}(T)$ on the training set,
$\mathrm{KS}_{i}$ is defined as follows:

$\mathrm{KS}_{i}=\sup_{p\in[0,1]}\left|\
p-\frac{1}{N_{i}(T)}\sum_{k=1}^{N_{i}(T)}\mathds{1}_{[0,p)}(p_{i,k})\
\right|.$ (3.20)

##  4 Application to the Santander Cycles bike-sharing system

The proposed models were tested on data from the Santander Cycles bike-sharing
network in central London. Transport for London publicly released data about
journeys on the Santander Cycles network, publicly available at
https://cycling.data.tfl.gov.uk/ (powered by TfL). Each record includes
information about the start and end times, IDs of source and destination
stations, journey duration, a bike ID number, and an ID number for the rental.
For the analysis in this section, 32 weeks of data were considered, starting
from 2nd March 2022 until 9th October 2022. A total of $\numprint{4098756}$
journeys recorded between 2nd March 2022 and 21st June 2022 (approximately 16
weeks) were used to construct the training set, whereas the remaining
$\numprint{3995591}$ journeys occurred between 22nd June 2022 and 9th October
2022 (approximately 16 weeks) formed the test set. A total of $M=787$ unique
active docking stations were observed in the data under consideration.
Distances between docking stations were calculated using the haversine
distance formula in Equation (3.2), with the value $\rho=6365.079\mathrm{km}$
(corresponding to the Earth radius at the latitude of London, with no
altitude). A more precise estimate of the cycling distance between two nodes
could be obtained via calls to the TfL Journey Planner API, or any other tool
providing calculations of cycling distances between locations. The main
practical limitation for this approach is that the calculation of the full
distance matrix would require ${M\choose{2}}$ API requests, which normally
largely exceeds the number of permitted requests under a basic plan.

The proposed GB-MEP model in Equation 2.3 has been compared to a number of
alternatives:

  * •

A standard Poisson process model for each station (Poisson):

$\lambda_{i}(t)=\lambda_{i}.$

  * •

A mutually exciting model for each station (MEP):

$\lambda_{i}(t)=\lambda_{i}+\sum_{h=1}^{N_{i}^{\prime}(t)}\alpha_{i}^{\prime}\exp\\{-\beta_{i}^{\prime}(t-t_{i,h}^{\prime})\\}.$

  * •

A self-exciting Hawkes process model for each station (SEP):

$\lambda_{i}(t)=\lambda_{i}+\sum_{k=1}^{N_{i}(t)}\alpha_{i}\exp\\{-\beta_{i}^{\prime}(t-t_{i,k})\\}.$

  * •

A self-and-mutually exciting process for each station (SMEP):

$\lambda_{i}(t)=\lambda_{i}+\sum_{k=1}^{N_{i}(t)}\alpha_{i}\exp\\{-\beta_{i}^{\prime}(t-t_{i,k})\\}+\sum_{h=1}^{N_{i}^{\prime}(t)}\alpha_{i}^{\prime}\exp\\{-\beta_{i}^{\prime}(t-t_{i,h}^{\prime})\\}.$

  * •

A spatial mutually exciting process for each station (SpMEP):

$\lambda_{i}(t)=\lambda_{i}+\sum_{j=1}^{M}\sum_{k=1}^{N_{j}(t)}\exp\\{-\theta_{i}\gamma(i,j)\\}\mathds{1}_{[0,\varepsilon)}\\{\gamma(i,j)\\}\cdot\alpha_{i}\exp\\{-\beta_{i}(t-t_{j,k})\\}.$

Note that all models above can be written as special cases of GB-MEP for
specific choices of the functions $\kappa$ and $\kappa^{\prime}$, cf. Equation
(3.3). The value of $\varepsilon$ in SpMEP and GB-MEP has been set to
$\varepsilon=0.5$. For some stations, $M_{i}^{\varepsilon}<3$ for
$\varepsilon=0.5$; in these cases, $\varepsilon$ was set to the minimum value
such that $M_{i}^{\varepsilon}\geq 3$. All model parameters for GB-MEP and its
competitors were estimated via maximum likelihood, using numerical
optimisation via L-BFGS-B (see, for example, Zhu et al., 1997). Note that the
Poisson model has a closed form solution:
$\hat{\lambda}_{i}=N_{i}(T^{\ast})/T^{\ast}$ for a training period
$[0,T^{\ast})$. The parameters of SEP and MEP were initialised to
$(\lambda_{i},\alpha_{i},\beta_{i})=(e^{-4},e^{-4},2e^{-4})$, and the
corresponding results were used to initialise SMEP. In turn, the results
obtained from SMEP were used to initialise SpMEP and GB-MEP.

The results are reported in Table 1, Figure 2(a) and Figure 3(a). Table 1
shows the Kolmogorov-Smirnov score for the distribution of $p$-values
calculated across all $M$ stations, calculated against the uniform
distribution. Figure 2(a) reports uniform quantile-quantile (Q-Q) plots for
all stations individually, and for all stations combined. Figure 3(a) contains
the boxplots of KS scores for all six models considered in this section. The
results show that GB-MEP has the best performance across all models, both on
the training and test set, demonstrating that taking into account the
geographical location of nodes within the network is important to improve the
goodness-of-fit of statistical models to real-world network point process
data. GB-MEP is the only model which achieves a median KS score below 0.05
across stations, both in the training and test sets cf. Figure 3(a)).
Additionally, Table 1 shows remarkable agreement between the results on the
training and test sets, which suggests that the models are not overfitting the
data.

Table 1: Kolmogorov-Smirnov scores for the distribution of $p$-values for all stations on training and test sets for six models. | Poisson | MEP | SEP | SMEP | SpMEP | GB-MEP
---|---|---|---|---|---|---
Training set | $0.2186$ | $0.1610$ | $0.0751$ | $0.0326$ | $0.0734$ | $0.0246$
Test set | $0.2217$ | $0.1584$ | $0.0752$ | $0.0350$ | $0.0738$ | $0.0264$
(a) Training set (b) Test set Figure 2: Uniform Q-Q plots of $p$-values,
calculated separately for each station on the training and test sets for six
models.

(a) Training set (b) Test set Figure 3: Boxplots of KS scores across all
station, calculated on the training and test sets for six models.

For most stations, the two best performing models are SMEP and GB-MEP, as
observed in Table 1, Figure 2(a) and Figure 3(a). GB-MEP has a better
performance than SMEP on approximately 73.88% of stations in the test set.
Figure 4(a) analyses more closely the differences between the results of GB-
MEP and SMEP, showing all docking stations in London, coloured by the
difference in KS scores on the test set, and a histogram of differences
between KS scores for SMEP and GB-MEP on the test set.

(a) Plot of station locations, coloured according to the difference in KS
scores between SMEP and GB-MEP on the test set (b) Histogram of differences
between KS scores for SMEP and GB-MEP on the test set Figure 4: Station
locations and histogram of differences in KS scores between SMEP and GB-MEP on
the test set.

The effect of the end times on the intensity of start times of journeys may
seem counterintuitive, but it is consistent with some characteristics of the
Santander Cycles BSS: for example, if a station is empty, a user returning a
bike to that station may imply that there is a higher chance of a journey
starting immediately from that station, with a different user picking up the
newly available bike. Also, until September 2022, users were charged £2 for
unlimited journeys with Santander Cycles in the subsequent 24 hours. For
journeys longer than 30 minutes, users are additionally charged. Therefore,
many users typically returned bikes before the 30 minutes mark, and they would
pick up another bike immediately for a new journey, meaning that returning a
bike to a station could increase the likelihood of a new journey starting from
the same station.

##  5 Conclusion

This article discussed a general strategy to incorporate measures of distance
between nodes on a graph within the context of multivariate Hawkes processes
on graphs, proposing GB-MEP, a graph-based mutually exciting point process.
Scalable inference procedures were discussed, using a recursive form of the
likelihood function under a scaled exponential form of the excitation kernel.
The model was fitted on a large dataset of journeys within the Santander
Cycles bike-sharing system, with over 4 million events in the training set,
and close to $M=800$ nodes. The model performance was evaluated on close to 4
million events comprising the test set, demonstrating superior performance of
GB-MEP over alternative modelling approaches.

A possible limitation of the proposed framework for the application to bike-
sharing systems is that a bound on the capacity of each station is not
included within the model: each station has a limited number of docks, which
could be used to correct the CIF when a station is at capacity. Also,
additional distance functions could be considered for the excitation
component. This work is framed within the context of docked bike-sharing
systems, but it could be more generally applied to mutually exciting processes
on any type of network where a distance metric between nodes is available. If
an adjacency matrix is observed, and not a distance matrix, measures of
similarity between nodes could be used for $\gamma$ within the GB-MEP
framework (for example, the length of the shortest path between nodes).

A possible extension of the version of GB-MEP presented in this work would be
to consider different excitation kernels, such as power laws (see, for
example, Laub et al., 2022). Furthermore, removing the separability assumption
on the excitation kernel could also prove interesting, as it would allow to
capture more complex patterns arising from the interaction between the spatial
and temporal components, which can lead to a better model fit. Unfortunately,
both alternatives would cause an increased computation cost. Finally, adding a
probabilistic mixture model component in the bike-sharing application could
further improve the model fit, reflecting the fact that bikes could be used
for different purposes, such as leisure or commuting.

## Code

A Python library with code and examples to reproduce all results in this
article can be found in the GitHub repository fraspass/gb-mep. The README file
in the repository provides detailed instructions on how to download the data,
preprocess them, and run the experiments.

## References

  * Blundell et al. (2012) Blundell, C., Beck, J., and Heller, K. (2012) Modelling reciprocating relationships with Hawkes processes.  In Advances in Neural Information Processing Systems 25 (eds. F. Pereira, C. Burges, L. Bottou and K. Weinberger), 2600–2608. Curran Associates.
  * Brown et al. (2002) Brown, E., Barbieri, R., Ventura, V., Kass, R., and Frank, L. (2002) The time-rescaling theorem and its application to neural spike train data analysis.  Neural computation, 14, 325–346.
  * Busatto and Stingo (2023) Busatto, C. and Stingo, F. C. (2023) Inference of multiple high-dimensional networks with the Graphical Horseshoe prior.  arXiv e-prints.
  * Choi et al. (2022) Choi, Y., Cho, H., and Son, H. (2022) Capturing usage patterns in bike sharing system via multilayer network fused Lasso.  arXiv e-prints.
  * Daley and Vere-Jones (2002) Daley, D. and Vere-Jones, D. (2002) An Introduction to the Theory of Point Processes – Volume I: Elementary Theory and Methods.  Probability and Its Applications. Springer.
  * Deng et al. (2023) Deng, P., Zhao, Y., Liu, J., Jia, X., and Wang, M. (2023) Spatio-temporal neural structural causal models for bike flow prediction.  Proceedings of the AAAI Conference on Artificial Intelligence, 37, 4242–4249.
  * Eichler et al. (2017) Eichler, M., Dahlhaus, R., and Dueck, J. (2017) Graphical modeling for multivariate Hawkes processes with nonparametric link functions.  Journal of Time Series Analysis, 38, 225–242.
  * Etesami et al. (2016) Etesami, J., Kiyavash, N., Zhang, K., and Singhal, K. (2016) Learning network of multivariate Hawkes processes: A time series approach.  In Proceedings of the Thirty-Second Conference on Uncertainty in Artificial Intelligence, UAI’16, 162–171. AUAI Press.
  * Foschi (2023) Foschi, R. (2023) A point processes approach to bicycle sharing systems’ design and management.  Socio-Economic Planning Sciences, 87, 101608.
  * Fox et al. (2016) Fox, E., Short, M., Schoenberg, F., Coronges, K., and Bertozzi, A. (2016) Modeling e-mail networks and inferring leadership using self-exciting point processes.  Journal of the American Statistical Association, 111, 564–584.
  * Gervini and Khanal (2018) Gervini, D. and Khanal, M. (2018) Exploring patterns of demand in bike sharing systems via replicated point process models.  Journal of the Royal Statistical Society Series C: Applied Statistics, 68, 585–602.
  * Hawkes (1971) Hawkes, A. (1971) Spectra of some self-exciting and mutually exciting point processes.  Biometrika, 58, 83–90.
  * Heydari et al. (2021) Heydari, S., Konstantinoudis, G., and Behsoodi, A. W. (2021) Effect of the COVID-19 pandemic on bike-sharing demand and hire time: Evidence from Santander Cycles in London.  PLOS ONE, 16, 1–16.
  * Huang et al. (2022) Huang, Z., Soliman, H., Paul, S., and Xu, K. S. (2022) A mutually exciting latent space Hawkes process model for continuous-time networks.  In Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence (eds. J. Cussens and K. Zhang), vol. 180 of Proceedings of Machine Learning Research, 863–873. PMLR.
  * Laub et al. (2022) Laub, P., Lee, Y., and Taimre, T. (2022) The Elements of Hawkes Processes.  Springer International Publishing.
  * Linderman and Adams (2014) Linderman, S. W. and Adams, R. P. (2014) Discovering latent network structure in point process data.  In Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32, ICML’14, II–1413–II–1421.
  * Miscouridou et al. (2018) Miscouridou, X., Caron, F., and Teh, Y. W. (2018) Modelling sparsity, heterogeneity, reciprocity and community structure in temporal interaction data.  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, NIPS’18, 2349–2358.
  * Mohler et al. (2011) Mohler, G. O., Short, M. B., Brantingham, P. J., Schoenberg, F. P., and Tita, G. E. (2011) Self-exciting point process modeling of crime.  Journal of the American Statistical Association, 106, 100–108.
  * Ogata (1988) Ogata, Y. (1988) Statistical models for earthquake occurrences and residual analysis for point processes.  Journal of the American Statistical Association, 83, 9–27.
  * Okawa et al. (2019) Okawa, M., Tanaka, Y., Kurashima, T., Toda, H., and Yamada, T. (2019) Marked temporal point processes for trip demand prediction in bike sharing systems.  IEICE Transactions on Information and Systems, E102.D, 1635–1643.
  * Paul et al. (2023) Paul, A., Flynn, K., and Overney, C. (2023) Estimating censored spatial-temporal demand with applications to shared micromobility.  arXiv e-prints.
  * Perry and Wolfe (2013) Perry, P. and Wolfe, P. (2013) Point process modelling for directed interaction networks.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75, 821–849.
  * Price-Williams and Heard (2020) Price-Williams, M. and Heard, N. A. (2020) Nonparametric self-exciting models for computer network traffic.  Statistics and Computing, 30, 209–220.
  * Reinhart (2018) Reinhart, A. (2018) A Review of Self-Exciting Spatio-Temporal Point Processes and Their Applications.  Statistical Science, 33, 299–318.
  * Rennie et al. (2023) Rennie, N., Cleophas, C., Sykulski, A. M., and Dost, F. (2023) Analysing and visualising bike-sharing demand with outliers.  Discover Data, 1, 1.
  * Rizoiu et al. (2017) Rizoiu, M.-A., Xie, L., Sanner, S., et al. (2017) Expecting to be HIP: Hawkes intensity processes for social media popularity.  In Proceedings of the 26th International Conference on World Wide Web, WWW ’17, 735–744.
  * Sanna Passino and Heard (2023) Sanna Passino, F. and Heard, N. A. (2023) Mutually exciting point process graphs for modelling dynamic networks.  Journal of Computational and Graphical Statistics, 32, 116–130.
  * Shlomovich et al. (2022) Shlomovich, L., Cohen, E. A. K., Adams, N., and Patel, L. (2022) Parameter Estimation of Binned Hawkes Processes.  Journal of Computational and Graphical Statistics, 31, 990–1000.
  * Transport for London (2022) Transport for London (2022) London’s iconic Santander Cycles hire scheme sees yet another record-breaking year in 2022.  URL: https://bit.ly/3PPUwbg.
  * Vallez et al. (2021) Vallez, C. M., Castro, M., and Contreras, D. (2021) Challenges and opportunities in dock-based bike-sharing rebalancing: A systematic review.  Sustainability, 13.
  * Xu et al. (2016) Xu, H., Farajtabar, M., and Zha, H. (2016) Learning Granger causality for Hawkes processes.  In Proceedings of The 33rd International Conference on Machine Learning (eds. M. F. Balcan and K. Q. Weinberger), vol. 48 of Proceedings of Machine Learning Research, 1717–1726. New York, New York, USA: PMLR.
  * Yuan et al. (2019) Yuan, B., Li, H., Bertozzi, A. L., Brantingham, P. J., and Porter, M. A. (2019) Multivariate spatiotemporal Hawkes processes and network reconstruction.  SIAM Journal on Mathematics of Data Science, 1, 356–382.
  * Zhao et al. (2015) Zhao, Q., Erdogdu, M. A., He, H. Y., Rajaraman, A., and Leskovec, J. (2015) SEISMIC: A self-exciting point process model for predicting tweet popularity.  In Proceedings of the 21th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, 1513–1522.
  * Zhu et al. (1997) Zhu, C., Byrd, R. H., Lu, P., and Nocedal, J. (1997) Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-Scale Bound-Constrained Optimization.  ACM Transactions on Mathematical Software, 23, 550–560.